A column generation algorithm for solving energy system planning problems

Energy system optimization models are typically large models which combine sub-models which range from linear to very nonlinear. Column generation (CG) is a classical tool to generate feasible solutions of sub-models, defining columns of global master problems, which are used to steer the search for a global solution. In this paper, we present a new inner approximation method for solving energy system MINLP models. The approach is based on combining CG and the Frank Wolfe algorithm for generating an inner approximation of a convex relaxation and a primal heuristic for computing solution candidates. The features of this approach are: (i) no global branch-and-bound tree is used, (ii) sub-problems can be solved in parallel to generate columns, which do not have to be optimal, nor become available at the same time to synchronize the solution, (iii) an arbitrary solver can be used to solve sub-models, (iv) the approach (and the implementation) is generic and can be used to solve other nonconvex MINLP models. We perform experiments with decentralized energy supply system models with more than 3000 variables. The numerical results show that the new decomposition method is able to compute high-quality solutions and has the potential to outperform state-of-the-art MINLP solvers.


Energy system optimization
Energy systems involve the conversion, distribution, and storage of energy. They represent complex structures of interconnected subsystems, including different technologies, actors, and markets, depending on the problem under consideration and the defined system boundary (Möst et al. 2009). Energy system studies range from technology comparisons for domestic applications (Ashouri et al. 2013) to the design of long-term transnational transformation strategies (Brown et al. 2018). When mathematical optimization methods are applied to the analysis, it is the modeler's responsibility to select an appropriate level of detail for the investigated problem to ensure the computability of the model and the significance of the obtained results (Brooks and Tobias 1996). Hence, the modeler has to address the trade-off between model accuracy and the required solution effort. In the optimization-based process synthesis, Chen and Grossmann distinguish between aggregated models, shortcut models, and rigorous models (Chen and Grossmann 2017).
Aggregated models typically represent the underlying physical relationships in a simplified way. The problems under investigation are often embedded in large energy systems, and usually, the plant operation is considered in more detail. The resulting models are typically linear (LP) or mixed-integer linear (MILP), and commercial and non-commercial software is available to solve them effectively, e.g., (Gurobi Optimization 2020), (CPLEX IBM 2020), (CBC Forrest et al. 2020) or (SCIP Gamrath et al. 2020). Examples include design concepts for central plants to supply district heating networks (Bruche and Tsatsaronis 2019), industrial parks with electricity, heating, and cooling demands (Voll et al. 2013), and decentralized urban areas with consideration of heat network expansion (Rieder et al. 2014). Rigorous models, in contrast, feature a significantly higher level of detail of the modeled technical devices and their fundamental principles. They represent the plant behavior at the component level and take into account, for example, chemical equilibria or heat transfer phenomena. The models are typically developed using special simulation software, and for optimization, derivative-free black-box methods such as evolutionary algorithms are frequently employed. These often time-consuming techniques enable investigation of complex models but have the drawback that they do not provide quality proof of the results and may produce only local optima. Examples are the design optimization for supercritical steam power plants (Wang et al. 2014) and natural gas-fired oxyfuel plants (Teichgraeber et al. 2017).
The shortcut models are a compromise between those mentioned above. They can either be large and thus closer to the aggregated models, but take into account, for instance, more advanced (nonlinear) functions for describing the component performance or costs, e.g., (Elsido et al. 2017;Goderbauer et al. 2016). Alternatively, the models can be more comparable to rigorous models with reduced complexity to enable deterministic solving methods to handle them. Model simplification can be achieved by considering only a few characteristic load points to limit the number of variables, and introduce model equations with lower mathematical complexity, e.g., to describe the thermodynamic properties (Ahadi-Oskui et al. 2010;Jüdes et al. 2009). In both cases, the results are commonly nonlinear mathematical models (NLP). Most of them are also nonconvex and have additional integer variables to state, for example, if a component is present or absent in a design. In this case, the problem class is mixed-integer nonlinear (MINLP). There are several deterministic methods available for solving nonconvex MINLPs. They can be subdivided into (i) branch-and-bound (Burer and Letchford 2012;Sahinidis 2020), (ii) MIP approximation Goderbauer et al. (2016), and (iii) decomposition methods (Nowak et al. 2019;Muts et al. 2020b). Most of the available commercial and open-source solvers, e.g., (BARON Sahinidis 2020) or (SCIP Gamrath et al. 2020) apply the branch-and-bound method. However, finding (high quality) solutions for large-scale MINLP problems with branch-and-bound-based solvers is often difficult with limited time and memory resources, as the size of the branch and bound tree to be managed may grow with the problem dimension.
Hence, a second approach in dealing with the challenging nonlinearities in optimization models is introducing piecewise-linear MIP approximations for the nonlinear functions. The MIP approximations are usually applied manually during the model preprocessing. The variety of applications includes combined cooling, heat, and power plants (Bischi et al. 2014), energy systems with heat demands at different temperature levels (Yokoyama et al. 2002), and electrochemical conversion technologies (Gabrielli et al. 2018). However, a MIP approximation solution is not an exact solution of the nonlinear model and might even be infeasible in the original problem. Goderbauer et al. addressed this challenge with a problem-tailored adaptive algorithm for the synthesis of decentralized energy supply systems (DESS) . Their iterative two-step process combines the identification of approximate solutions with discretized variables and the solving of NLPs with fixed selected variables. Between two iterations, the discretization grid is reduced and shifted in the direction of the previously calculated solution. With this approach, Goderbauer et al. generated high-quality MINLP solutions in short computation time and outperformed the commercial solver BARON. However, we should take into account that unlike BARON, the method is not a general purpose MINLP solver, but it has been tuned to solve specifically DESS problems.
A third deterministic method for solving MINLPs is decomposition. The authors have gained very good experience with problem decomposition via Column Generation (CG) in solving huge network problems (Borndörfer et al. 2013;Nowak 2014). Motivated by CG's excellent performance, we present a new MINLP CG method. It is based on combining CG and the Frank-Wolfe method (FW) for generating an inner approximation (IA) of a convex relaxation, and a primal heuristic for computing solution candidates. The characteristics of this approach are: 1 3 (i) it does not require a global branch-and-bound tree, (ii) sub-problems can be solved in parallel to generate columns, which do not have to be optimal, nor become available at the same time to synchronize the solution, (iii) an arbitrary solver can be used to solve sub-models.
The new method can be applied to general nonconvex or convex sparse real-world MINLPs, e.g. instances of the MINLPLib (Vigerske 2018). These models can be reformulated as block-separable problems, defined by low-dimensional sub-problems coupled by a moderate number of global linear constraints. In this work, we focus on the characteristics of the superstructure-based synthesis of decentralized energy supply system presented in Goderbauer et al. (2016).

Outline
First, Sect. 2 describes a practical application from the area of energy supply system modeling. Section 3 describes the decomposition view of the problem and sketches the resource-constrained approach. Section 4 presents various algorithms based on IA to solve the problem. The numerical evaluation in Sect. 5 shows the potential of the new decomposition-based approximation approach. We conclude in Sect. 6 outlining steps for future investigation.

Decentralized energy supply systems
A DESS is a complex, integrated system consisting of several energy conversion units and energy supply and demand forms. Usually, the starting point for DESS optimization is a superstructure containing all possible components and their interconnections. The optimization goal involves identifying a plant design and a suitable plant schedule simultaneously to minimize or maximize an objective function value. The objective function is usually of economic nature (e.g., maximizing the net present value of an investment). Binary variables (valid values one or zero) can be used to model the facility selection (existence/non-existence) and the operation states of selected units (on/off). Continuous non-negative variables represent the model component input and output streams and the component costs. Both input-output relationships and component costs can be expressed by nonlinear functions. In this case, the underlying model belongs to the class of mixed-integer-nonlinear-problems (MINLP).

General MINLP model formulation
The DESS model used in this paper was initially presented in Voll et al. (2013) and is described in detail in Goderbauer et al. (2016). For the following brief model description, the notation is adopted from Goderbauer et al. (2016). The superstructure of the model is sketched in Fig. 1.
Every superstructure S is assembled from the components s, namely boilers (B), combined heat and power engines (C), electrically powered compression chillers (referred to as "turbo chillers", T), and heat-utilizing absorption chillers (A). In the base case, each component can only be installed once in a superstructure (see Sect. 2.3 for more details). The DESS is required to meet the heating ( Ė heat ), cooling ( Ė cool ), and electricity demands ( Ė el ) for every load case l. The latter are specified by a set of load cases L. Gas-fired boilers and combined heat and power engines are available to provide thermal energy (gas price: p gas,buy ). The electricity produced in the engines is used to satisfy the electricity demand. Missing or excess electricity can be balanced with the connected electricity grid. Costs of p el,buy are incurred for purchased electricity, and the electricity sale is reimbursed with p el,sell . Electrically powered turbo-chillers and heat-utilizing absorption chillers are available to meet the cooling demand. The maximization of the net present value is selected as the objective function. The annual costs are discounted to the investment time using the interest rate i and the cash flow time CF . The factor m s considers the annual maintenance costs as a percentage of the investment costs I s .
The following variable values are determined during the optimization process: The objective function of the optimization problem is presented in Eq. (1). The net present value to be maximized is calculated from the discounted annual cash flows minus the overall investment costs. The annual cash flow is derived from the revenues for electricity sales minus the expenses for electricity import, purchased fuel, and maintenance costs.
Equations (2)-(4) represent the balance constraints for the three energy demands heating, cooling, and electricity. They need to be met with equality in all load cases l.
All components are available in continuous nominal sizes. Equation (5) bounds their size to the range between a lower and upper limit.
Equations (6)-(8) are responsible to set the output energy flow to zero if the binary operation variable is zero. Otherwise the output flow is limited to values between the minimum part-load ( min s ⋅V N s ) and the nominal size. Eq. (9) assures that every unit of a superstructure can only be operated if it exists. (1) A column generation algorithm for solving energy system planning problems All parameter values, the nonlinear part-load performance curves U s (V s ,V N s ) , and investment cost functions I s (V N s ) are specified in the Appendix. As stated by Goderbauer et al. (2016) the MINLP model is nonconvex due to the investment cost functions and the nonlinear equality constraints (2) and (4).

DESSLib model instances
The DESS model studied here is published as library DESSLib  including multiple model instances of varying complexity. The instances differ in three properties: • number of similar units of all component types (S4, S8, S12, S16) • number of load cases (L1, L2, L4, L6, L8, L12, L16, L24) • energy demands values (ten sampled variants) For example, instance S8L4-0 consists of maximal two units for each of the four component types (boiler, engine, turbo-chiller, absorption chiller). Both units can have different nominal sizes and can be operated independently. Furthermore, instance S8L4-0 includes four load cases (L4), and uses the data variation with the index zero. The DESSLib is used in Goderbauer et al. (2016) to evaluate their problem-tailored method AdaptDiscAlgo (see explanation in Sect. 1) and to benchmark it with state-of-the-art solvers. The results are also accessible for download at Bahl et al. (2016). Thus, high-quality, near-optimal solutions for the applied MINLP problems are available for all model instances.

MINLP block model formulation
Our model has been developed with the Python-based algebraic modeling language Pyomo. Pyomo provides the component Block, which serves as a container for variables, constraints, and more (Hart et al. 2017). For all energy conversion units of the DESS model, i.e., for each boiler, engine, turbo-chiller, and absorption chiller, individual Pyomo blocks are created. These main blocks contain the component variables for existence, nominal size, and investment costs, and the associated nonlinear investment cost functions. Additionally, one or multiple sub-blocks are created within the main component block to store the variables and constraints with load case dependency. These include operation status variables, input and output flow variables, plus the nonlinear performance constraints, and Eq. (6)-(9). Each sub-block can contain one or more load cases. Thus, the size of the sub-blocks is adjustable and depends on the decision of the model user. A schematic representation of the modeling concept is shown in Fig. 2. Modeling with Pyomo blocks is necessary to support the automatic model decomposition. Further details are provided in Sect. 3.1.

Block-separable formulation
In order to apply the new decomposition method, the DESS model is formulated as a block-separable (or quasi-separable) MINLP problem, where the vector of variables x ∈ ℝ n is partitioned into |K| blocks of dimension n k such that n = ∑ k∈K n k and x k ∈ ℝ n k denotes the variables of block k. The model is described by with global (coupling) linear constraints and local (block) constraints with n k = n k1 + n k2 and The local (block) set X k is defined by linear and nonlinear local constraint functions, g kj ∶ ℝ n k → ℝ , which are assumed to be bounded and continuously differentiable within the set [x k , x k ] . An example of a nonlinear constraint function is depicted in Fig. 3. The linear global constraints in P are defined by . The linear objective function is defined by c T x ∶= ∑ k∈K c T k x k , c k ∈ ℝ n k . Blocks in set K can be detected in an automatic way or can be defined manually. One idea of automatic block detection is based on connected components of the socalled 'sparsity graph', defined by nonzero entries of the Hessian of constraint functions of the original MINLP (Nowak 2005). In this paper, we consider the MINLP problems (DESS models) with manually defined blocks, see Sect. 2.4 for implementation details. Some of these blocks are defined by continuous variables and by linear local constraints. We merge these blocks into one one single linear block which is represented by the following polytope Further details how we deal with linear block (14) are given in Sect. 4.1.1.

Resource-constrained reformulation
In order to define an inner approximation, the MINLP (10) is reformulated as a resource-constrained optimization problem by an affine transformation of the original variables x k , k ∈ K, into resource variables (Muts et al. 2020a) 1 3 The variables w k , k ∈ K, describe how the objective value and the constraint values a T i x are distributed over blocks. Considering the transformed local feasible set the resource-constrained formulation of (10) is defined as with the transformed objective function and the transformed global feasible set Notice that the resource constraint view considers a space with (m + 1)|K| variables instead of the original n variables. For sparse MINLPs with sparse matrices A k the number of resources variables can be greatly reduced, since only a resource variable w ki is needed, if A ki ≠ 0 . For instance, let m = 1499, |K| = 50 and number of nonzero rows of A k , k ∈ K to be 100. Then the overall number of resource variables is 75000. This number is reduced to 5000, if we consider only the nonzero rows of matrix A k .

Convex hull relaxation and inner approximation
A resource-constrained convex hull relaxation of (17) is defined by The quality of relaxation (19) depends strongly on the duality gap, defined by Let S k ∶= {y kj } j∈[S k ] ⊂ X k be a finite set of locally feasible points, called inner approximation of X k , and define the corresponding resources by R k ∶= {A k y ∶ y ∈ S k } , called columns. An inner approximation of (19), called LP-IA, is defined by 1 3 A column generation algorithm for solving energy system planning problems We use the notation [R] ∶= {1, … , |R|} for the index set of the finite set R, i.e.
} . An LP formulation of (21) is given by where We call z kj the weight of column r kj ∈ R, k ∈ K and denote by |R k | ⊂ ℝ |R k | the standard simplex with |R k | vertices in the following way

A MINLP column generation algorithm
In this section we describe a new MINLP CG method. The method, called minlpcg, solves MINLP sub-problems to generate partial feasible solutions x k ∈ X k and w k ∈ W k , which in the end are combined into a global solution which also fulfills the global constraints. The direction for looking for better partial solutions are given by master problems in the resource space. The basic steps of minlpcg are: 1. Initialize column sets R k , k ∈ K. 2. Solve (21) for computing an approximate solution w of (19) . 3. Project w onto the feasible set finding a solution candidate w * . 4. Add new columns to R k , k ∈ K. 5. Repeat from step 2 until stopping criterion.
Since the sub-models of (10) represent difficult MINLP problems, it is not easy to quickly generate many columns. Therefore, we developed several strategies for accelerating the column generation (CG).

Dealing with linear block (sub-problem), adapting inner approximation
In the DESS models the dimension n 1 of the linear block (14) is much higher than the dimension n k of other blocks. Preliminary numerical experiments showed that, if we distinguish the linear block instead of running traditional CG over all blocks, we obtain a running time reduction in order of magnitude from 2 days towards 1-2 hours. Mathematically this can be explained by the many vertices that the polytope of a linear block has which should be generated as columns if treated as nonlinear blocks. Therefore, we use the following modified LP-IA master problem (25), for which the linear constraints of X 1 in (14) are directly integrated in the LP-IA: where w k (z k , x 1 ) is defined similarly as in (23) by So we distinguish K ⧵ {1} as the set of nonlinear blocks whereas x 1 ∈ ℝ n 1 are the variables of the linear block, and X 1 is defined by only local linear constraints as in (14).
Lemma 1 Let R 1 = vert(W 1 ) be the set of extreme points of W 1 . Then problems (25) and (21) are equivalent.
Proof By definition (14), X 1 defines a polytope. Hence, W 1 is a polytope defined by a linear transformation of X 1 , and W 1 = conv (R 1 ) . Therefore, A 1 x 1 in (26) can be replaced by w 1 (z 1 ) , which proves the statement. ◻ The procedure solveinnerlp(R) solves (25). It returns a primal solution (x k ) k∈K , where x 1 is the solution of the linear block and x k = x k (z k ) with and S k is the set of generated feasible points of X k related to R k , i.e. r kj = A k y kj . Moreover, the procedure returns the dual values for the global resource constraints.

Generation of columns solving sub-problems
The MINLP CG-algorithm generates columns by solving the following MINLP sub-problems regarding a search direction d ∈ ℝ m+1 , where d is typically defined by a dual solution ∈ ℝ m of the LP-IA (21), i.e. d = (1, T ) . Notice that the result y k corresponds to an extreme point of X k as well as W k and is a so-called supported Pareto point in the resource space, see Muts et al. (2020b).
The procedure solveMinlpSubProbl(d) solves (28) and is used in procedure addCol(d, R k ), described in Algorithm 1, to add a column A k y k to R k . Moreover, the procedure computes the reduced cost which is used later to measure the impact of the procedure. If k < 0 for some k ∈ K , then column A k y k may improve the objective value of (25). In the other case, if = 0 , the objective value of (25) cannot be changed (Nowak 2005) and the column generation algorithm can be stopped.
Algorithm 2 shows a method for generating columns by alternately solving an LP-IA master problem and generating columns using Algorithm 1. The set K ⊆ K is subset of the block set K. Note that columns can be generated in parallel.

Initializing columns
In order to initialize column set R k , we perform a subgradient method using Algorithm 3 for maximizing the dual of (10) regarding the global constraints: We compute the step length p by comparing the values of the function L( ) defined in (30) at different iterations p of Algorithm 3, similarly as in Shor (1985): Note that procedure addcol in Algorithm 3 can be performed in parallel.

Accelerating column generation
In order to accelerate CG without calling an MINLP routine, we developed two methods. The first method generates columns by performing NLP local search from good starting points of an IA-LP. The second method is a Frank-Wolfe algorithm for quickly solving the convex hull relaxation (19).

Fast generation of columns using NLP local search and rounding
The master problem used to provide starting points for the column generation is LP-IA problem (21). Since in the beginning only few columns are available, often LP-IA master problem (21) is infeasible, i.e. H ∩ ∏ k∈K conv (R k ) = � . Therefore, we use an LP-IA master problem (32), that includes slacks and a penalty term ( , s) , similar as in du Merle et al. (1999): where the penalty term for slack variables is min F(w(z, x 1 )) + ( , s) and penalty weights > 0 are sufficiently large. We define w k (z k , x 1 ) similarly to (23) as follows and x 1 ∈ ℝ n 1 are the variables of the linear block, defined in (14). The procedure solveslackinnerlp(R) solves (32)  with e i ∈ ℝ m the coordinate i unit vector.
Since, for the CG algorithm, it is sufficient to compute high-quality local feasible solutions, we present a local search procedure approxsolveminlpsubprobl in Algorithm 4 based on rounding of locally feasible point. The goal of this procedure is to avoid usage of a MINLP solver for solving sub-problems and, therefore, reduce the time for sub-problem solving. The inputs of local search procedure approxsolveminlpsubprobl are the block solution x k as starting point and the direction d or (1, ) as search direction. It starts by running procedure solvenlpsubproblem which computes a local minimizer of the integer relaxed sub-problem starting from the primal solution x k of the LP-IA. Then the procedure round rounds integer variables of block k in ỹ k to obtain x k . Finally, procedure solvefixednlpsubproblem solves again an NLP problem, fixing the rounded integer variables of x k : and using the continuous variable values of x k as starting point. The complete column generation procedure is depicted in Algorithm 5. Note that procedure approxsolveminlpsubprobl in Algorithm 5 can be performed in parallel.

CG by solving the convex hull relaxation using a Frank-Wolfe algorithm
In this section, we present a Frank-Wolfe algorithm which is an alternative way to generate columns. It is based on solving convex hull relaxation (19) by a quadratic penalty function approach: where ∈ ℝ m + is a vector of penalty weights. Consider the convex optimization problem Let * be an optimal dual solution of (19) regarding the global constraints w ∈ H and set the penalty  (37).
Algorithm 6 presents a Frank-Wolfe (FW) algorithm for approximately solving the convex penalty problem (37). For acceleration, we use the Nesterov direction update rule (Nesterov 1983), line 17. We set the penalty weight = | | , where is a dual solution of LP-IA (32). One step of the FW-algorithm is performed by approximately solving the problem with a linearized objective which is equivalent to solving the sub-problems

3
A column generation algorithm for solving energy system planning problems The sub-problem (39) is solved with approxsolveminlpsubprobl, depicted in Algorithm 4, in order to compute quickly new columns. The columns can be computed in parallel. Note that the gradient ∇ w k Q(w, ) is defined by Hence, ∇ w k Q(w, ) = (1, (w, ) T ) for all k ∈ K . The quadratic line search problem in step 14 of Algorithm 6 can be easily solved.

A primal heuristic for finding solution candidates
In this section, we present two heuristic procedures for computing solution candidates. The first one computes a feasible solution from the slack LP-IA problem (32) solution. The second one computes high-quality solution candidates.
Algorithm 7 presents the initial primal heuristic, which aims at eliminating slacks in LP-IA master problem (32). It starts with running procedure nlpresourceproject, which performs an NLP local search solution of the following integer relaxed resource-projection NLP master problem where x is the solution of the LP-IA master problem (32).
Using the potentially fractional solution ỹ of (40), the algorithm computes an integer globally feasible solution ŷ by calling the procedure mipproject(ỹ ) which solves MIP-projection master problem The integer globally feasible solution ŷ is then used to perform an NLP local search, where integer variables are fixed starting from ŷ by procedure solvefixednip(ŷ): Algorithm 8 presents a primal heuristic for computing a high-quality solution candidate of MINLP problem (10). The procedure is very similar to Algorithm 7, but it does not use the NLP resource-projection problem (40). Instead, the solution of LP-IA (25) is used directly in MIP-projection problem (41). There is no guarantee that the optimal solution of (41) provides the best primal bound. Moreover, it may be infeasible for the original problem (10). Therefore, we generate a pool Ŷ of feasible solutions of problem (41) provided by the MIP solver. Solution pool Ŷ provides good starting points for an NLP local search over the global space and increases the possibility of improving the quality of solution candidate. Similarly to Algorithm 7, Algorithm 8 starts with computing an inner LP solution x of problem (25) by calling procedure solveinnerlp. Point x is used in a procedure solpoolmipproject(x, N ) to generate solution set (pool) Ŷ of (41) of size N, which also includes the optimal solution. Like in Algorithm 7, those alternative solutions are used to perform an NLP local search over the global space, defined in (42) by fixing the integer valued variables. In order to find better solution candidates, these steps are repeated iteratively with an updated point x . In each iteration, the point x is shifted towards the point x * which corresponds to the best current primal bound of the original problem (10). This is a typical heuristic local search procedure, which aims to generate a different solution pool Ŷ in each iteration of the algorithm. Algorithm 8 terminates when the maximum number of iterations is reached or the best primal bound in the current iteration does not improve the best primal bound in the previous iteration.

Main algorithm
Algorithm 9 describes a MINLP-CG method for computing a solution candidate of (10). The algorithm starts with the initialization of IA with the procedure iainit (Algorithm 3). Since the problem (32) might have nonzero slack values, the algorithm tries to eliminate them by computing a first primal solution. This is done by alternately calling the procedures approxcolgen (Algorithm 5) and findsolutioninit (Algorithm 7). For quickly improving the convex relaxation (19), the algorithm calls FW-based column generation procedure fwcolgen (Algorithm 6).
In the main loop, the algorithm alternately performs colgen (Algorithm 2) and heuristic procedure findsolution (Algorithm 8) for computing solution candidates. The procedure colgen is performed for a subset of blocks K ⊆ K ⧵ {1} , in order to keep the number of solved MINLP sub-problems low. Moreover, focusing on a subset of blocks helps avoiding computing already existing columns. The blocks can be excluded for a while by looking at the value of the reduced cost k , k ∈ K ⧵ {1} , which is computed in line 14, as defined in (29). The reduced block set K contains the blocks where the reduced cost is negative, i.e. k < 0, k ∈ K ⧵ {1} , and is updated at each main iteration by solving the sub-problems for the full set K. Note that if the reduced cost k is nonnegative for all blocks, i.e. k ≥ 0, k ∈ K ⧵ {1} , the Column Generation has converged and the algorithm terminates.

Convergence analysis
Column Generation and Frank-Wolfe algorithms are well-known approaches and their convergence has already been proven. In this section, we discuss the convergence of Column Generation in Algorithm 9 and the Frank-Wolfe method in Algorithm 6.

Convergence of Algorithm 9
The convergence proof of the Column Generation algorithm is due to its equivalence to the dual cutting-plane algorithm, see Lemma 4.10 in Nowak (2005). Note that the proof is not based on the computation of reduced cost k , k ∈ K , defined in (29). However, it can be used for measuring the impact of the new columns and as a criterion for algorithm termination, see Nowak (2005). For the convergence proof, we assume that all LP master problems (21) and MINLP sub-problems (28) are solved to optimality. Since Algorithm 2 is performed regarding the subset of the blocks K , we ensure that main Algorithm 9 converges by performing a standard CG step in line 14 regarding all blocks. Note that the direct integration of the linear block (14) into LP-IA master problem (25) is equivalent to performing the CG algorithm regarding that block until convergence, as shown in Lemma 1.  at the p-th iteration of Algorithm 9 at line 12 and * be the optimal value of convex hull relaxation (19). Then lim p→∞ c T x p = * .

Proposition 1 Let x p be the solution of
Proof The proof is equivalent to the proof of Proposition 4.11 in Nowak (2005). ◻

Convergence of Algorithm 6
Algorithm 6 combines the original Frank-Wolfe algorithm (see Algorithm 1 in Jaggi (2013)) with Nesterov update rule (Nesterov 1983). The approach proposed by Nesterov in Nesterov (1983) has a convergence rate O(1∕p 2 ) , whereas the original Frank-Wolfe algorithm has a slower convergence rate of O(1/p), see Theorem 1 in Jaggi (2013). In order to prove the convergence of Algorithm 6, we assume that all subproblems (39) are solved to global optimality. Next, we state that Algorithm 6 has a convergence rate O(1∕p 2 ).
Proposition 2 Let p ∶= Q(w p , ) be the value of the quadratic penalty function (36) at p-th iteration of Algorithm 6 and * be the optimal value of convex hull relaxation (19).
, where * defines an optimal dual solution of (19). Then there exist a constant C such that ∀p ≥ 0 Proof The proof is equivalent to the proof of Theorem 1 in Mouatasim and Farhaoui (2019). ◻

Numerical results
In this section, we evaluate the performance of Algorithm 9 solving several DESS model instances, described in Sect. 2. Algorithm 9 was implemented with Pyomo (Hart et al. 2017), an algebraic modeling language in Python, as part of the parallel MINLP-solver (Decogo Nowak et al. 2018). In the current implementation of Decogo, the sub-problems are not solved in parallel. The solver utilizes SCIP 7.0.0 (Gamrath et al. 2020) for solving MINLP sub-problems, GUROBI 9.0.1 (Gurobi Optimization 2020) for solving MIP/LP master-problems and IPOPT 3.12.13 (Wächter and Lorenz 2006) for performing an NLP local search in master and subproblems. All computational experiments were performed using a computer with AMD Ryzen Threadripper 1950X 16-Core 3.4 GHz CPU and 128 GB RAM. For the experiments with DESS instances, the blocks K were defined manually using Pyomo Block component (Hart et al. 2017), for more details see Sect. 2.4. For the experiments with Algorithm 9, we use the following stopping criteria, unless another is mentioned: • Time and iteration limit in Algorithm 9 is set to 12 hrs and 20 iterations, respectively; • The iteration limit of Algorithm 2 is set to 5; • The iteration limit of the outer and inner loop of Algorithm 6 is set to 5 and 10, respectively; • In Algorithm 8, the iteration limit is set to 5, the pool size to N = 100 and = 0.5 . To generate a solution pool with Gurobi, we used a parameter value of poolsearchmode=1. In this approach, the solver computes a set of N high-quality solutions including an optimal one, see Gurobi Optimization (2020) for more details; • For SCIP we set the maximum number of processed nodes after the last improvement of the primal bound to 1000, since, for CG, it is sufficient to compute good feasible solutions, for more details see Muts et al. (2020a).
Our main research question is how to exploit the decomposable structure of DESS models to generate a new approach which may be an alternative to state-of-the-art software like BARON. During the investigation, we developed several elements to speed up the computation. To evaluate their effect, first in Sect. 5.1 we measure the effect of distinguishing the linear block and adding the fast FW approach for generating promising columns. The idea to use a pool of high-quality MIP solutions is numerically investigated in Sect. 5.2. The pool generates a group of starting points to look for feasible solutions. That does not only affect the quality of the reached best solution, but also provides alternatives for the design problem that may be analyzed by the decision maker. Finally in Sect. 5.3, we compare the outcomes and the solution process of Algorithm 9 to that of the state-of-the-art solver BARON and to the problem-tailored MIP approximation algorithm presented in Goderbauer et al. (2016).

Effect of linear block integration into LP-IA and fast CG with the Frank-Wolfe approach
To illustrate the impact of direct integration of one linear block into the LP IA master problem, as defined in (25), we compare two runs of Algorithm 9 on the same instance: The first run uses the LP-IA formulation given by (25) and the second uses the formulation given by (22). The comparison was done for instance S16L16-1, see Table 6 for more details.
For this particular instance, approximately 37% of all variables belong to the linear block. Whereas the CG algorithm for formulation (25) relies on generating linear block polytope vertices as columns, the second formulation (22) may save this effort. The difference can be observed in Fig. 4; the running time of Algorithm 9 with the direct linear block integration is drastically reduced from two days to approximately two hours. It also shows a significant improvement of the convergence speed of the IA objective value and a significant improvement of its final value.
The next question focuses on the the convergence impact of including procedure fwcolgen (Algorithm 6) into Algorithm 9. To measure the effect, we perform two equivalent runs of Algorithm 9; one of those does not use procedure fwcolgen. We only focus on the convergence of the IA and don't perform procedure findsolution (Algorithm 8). The test instance is S4L4-1, see Table 6 for more details. Figure 5 shows that the IA objective value (25) converges faster with procedure fwcolgen at the initial stage of Algorithm 9. Moreover, it appears that in the beginning of the algorithm, the fwcolgen IA objective value is very close to the final IA objective value, i.e. after 50 seconds, the IA objective value with procedure fwcolgen is approximately 0.5 % worse than the final IA objective value. Fig. 4 Convergence of the IA objective of Algorithm 9 with and without linear block integration for instance S16L16-1 Fig. 5 Convergence of the IA objective value of Algorithm 9 for S4L4-1 with and without the fast FW column generation Figure 5 indicates that Algorithm 9 with procedure fwcolgen converges slower to the final IA objective value than without procedure fwcolgen. The reason for that is that fwcolgen fails to generate high-quality columns in the later stages of the algorithm. This is due to the fact that we have no guarantee that penalty weights , defined in (36), are sufficiently large. Another reason is that the procedure generates columns by solving NLP sub-problems heuristically. However, it has an advantage in speed when comparing column generation by solving MINLP sub-problems.

Impact of using the solution pool in Algorithm 8
This section focuses on using the solution pool in the heuristic procedure findsolution (Algorithm 8). Note that using a single integer feasible point, we perform local search with the fixed NLP (42). This may imply that the result is not a feasible solution of the problem. Using a pool of starting points, may alleviate this effect.   Figure 6 shows the objective function values of all pool solutions, generated for instance S16L16-1. The solution with the best primal bound of approximately −5.15 × 10 7 is identified in the seventh iteration of Algorithm 9. The component configuration of this solution is specified in Table 1. In this solution, all four available units of all component types are selected. From Fig. 6, one can notice that other high-quality solutions could be computed in the earlier iterations of the algorithm. Note that the MIP solver was not always able to provide a solution pool of prescribed size of N = 100 . Instead it computed a solution pool of smaller size. Therefore, in Fig. 6 one can observe that the number of generated feasible solutions varies over the iterations. Moreover, Fig. 6 indicates that a solution pool may contain many feasible points with a very similar objective function value. These points are distinct, since GUROBI adds only a feasible point to the solution pool if the value of Fig. 7 DESS configurations for all solutions of the solution pool of instance S16L16-1 within a range of 2% from the solution with the best primal bound its integer variables is different for at least one integer variable Gurobi Optimization (2020).
One advantage of the solution method is that a variety of feasible solutions are generated during the process. These can be stored in the solution pool. Often nearoptimal solutions are valuable for the user to gain a more profound knowledge of the optimization problem Voll et al. (2015). In addition, if they better satisfy requirements not included in the mathematical model (e.g., safety, maintainability, and operability) near-optimal solutions might be finally selected in the design of energy systems instead of an optimal one Bejan et al. (1996). Figure 7 shows the DESS configurations of all solutions of the solution pool of instance S16L16-1 with a primal bound differing at most 2% from the best solution. In the case of instance S16L16-1, this range comprises 95 solutions. In Fig. 7 they are sorted by their primal bound values. Thus the solution from Table 1 is shown on the very left. First, it is apparent that all four available units of each component type are selected in almost all cases. In further investigations, it could be analyzed whether this indeed provides advantages regarding the objective function value or if this behavior results from the the solution procedure. For the boiler and the CHP engine, the total installed capacity is mainly focused on one or two units, and the other units are installed with a smaller or minimal size. For the cooling equipment, the total capacity tends to be more evenly distributed among all units. This may be explained by their distinct part-load behavior. The configuration of most solutions is relatively similar. In the considerably different solution numbers 10 and 12, it is indicated that a shift of the cold supply from heat-driven absorption chillers to a higher proportion of electric turbo chillers is possible with only small reductions in the objective function value. This simultaneously reduces the required boiler capacity and increases the capacity of the CHP engines to cover the increased internal electricity demand.

Comparison to other approaches
In this section, we compare the performance and solution quality of Algorithm 9 to other solvers on 12 selected instances from the DESSLib library . The test instances were selected varying the number of unit components and number of load cases, as described in Sect. 2.3. Instance characteristics are reported in Table 6. The results of Algorithm 9 are compared to results obtained by the state-ofthe-art MINLP solver BARON 20.4.14 (Sahinidis 2020). We also compare results provided by adaptive discretization MINLP algorithm (AdaptDiscAlgo) reported in Goderbauer et al. (2016). We do not compare the performance of Algorithm 9 to AdaptDiscAlgo, since the implementation of AdaptDiscAlgo is not publicly available. Therefore, we analyze the results by looking at the primal bounds of AdaptDis-cAlgo, reported in Bahl et al. (2016). Note that AdaptDiscAlgo does not provide a valid lower bound and cannot be applied to general MINLP problems.
In this analysis, we use a 12 hour time limit for BARON and Algorithm 9. For the quality evaluation of a feasible solution, we define the gap to a reference (base) objective function value b as where a is an objective value of the feasible solution point. Note that a "−" is used before the term in (43), since we are considering a maximization problem. In case of minimization, it is omitted. Table 2 compares the primal solution quality of decomposition Algorithm 9 versus that of BARON and AdaptDiscAlgo using gap value (43). Let CG and B be the primal bound of CG Algorithm 9 and BARON, respectively. For the sake of simplicity, let * denote the best known primal bound among BARON and AdaptDiscAlgo. Table 2 also presents the duality gap of Algorithm 9 and BARON. Let CG denote the dual bound of Algorithm 9 defined by the objective value of the last solution of (25) and B denote the dual bound of BARON given by the lower bound provided by BARON. Table 2 shows that for five instances, Algorithm 9 improves the primal bound of BARON. Even though AdaptDiscAlgo is specially tailored to DESS, the Algorithm 9 result differs at most 3 % from the best known bound and for one instance improves the solution. We also evaluated the solution after performing two main iterations with Algorithm 9. The solution quality is relatively good. This means that with a limited solution time, relatively good solutions can be generated after two iterations. For the instances with more than 500 variables, the duality gap of Algorithm 9 is smaller that of BARON. One can also observe that the duality gap value  of Algorithm 9 depends mostly on the number of unit components (parameter S) and is less sensitive to the number of load cases (parameter L). Table 3 gives the number of iterations N iter performed by Algorithm 9 and the running time of both algorithms, Algorithm 9 to BARON.
Despite the fact that Algorithm 9 does not solve the sub-problems in parallel and it is implemented with Python, Table 3 shows that for several larger instances it requires less time than BARON. Moreover, we can see that for some instances Algorithm 9 did not reach the maximum number of iterations and the time limit. This indicates that apparently the Column Generation has converged.   Table 3 We are interested in the computing time with increasing size of the problem instance without varying the number of iterations. The estimated running time of Algorithm 9 after two main iterations is sketched in Fig. 8. Moreover, it sketches the actual time for initializing the Inner Approximation (all procedures before entering into main iterations of Algorithm 9). Note that the total time for finishing two iterations includes the actual time for the CG initialization. The idea of time estimation for two main iterations is based on the assumption that, for each instance, in each call of findsolution (Algorithm 8), we obtain a solution pool by a MIP solver of size 500 ( 500 = 5 ⋅ 100 , where 5 is the prescribed number of iterations of Algorithm 8 and 100 is a prescribed size of the solution pool). In fact, the solution pool provided by a MIP solver can be smaller than the prescribed number, see Fig. 6. Therefore, the estimated time indicates how long it would take if we would need to call the NLP solver 500 times to solve fixed NLP problem (42) in the first two main iterations of Algorithm 9. This time could be also useful for estimating the time complexity of similar instances with a much larger size. From Fig. 8, one can notice that the time for initialization of Column Generation is relatively small in comparison to the estimated time after performing two main iterations of Algorithm 8. Figure 5 illustrates that the algorithm quickly computes a good CG bound in the initial stage of the algorithm. Figure 8 indicates that findsolution (Algorithm 8) has a big influence on the running time of the entire algorithm.

Conclusions
Energy system optimization models combine potentially difficult sub-models into a global system. Due to their dimension, such systems may be difficult to solve by generic solvers based on a single branch-and-bound tree like BARON. Decomposition appears to be an appropriate concept to apply to such models due to the structure of sub-models and global constraints. Our research question is how to do this in an efficient way.
In our investigation, we looked into the potential of using Column Generation (CG). The idea is to generate feasible solutions of sub-models, defining columns of a global master problem, which is used to steer the search for a global solution. The master problem is based on an easy to solve inner approximation (IA) of a convex hull relaxation. One of our findings is, that it is more efficient to deal with blocks that have a linear structure in a separate way and include them directly in the master problem instead of generating columns for them. For speeding up column generation, we developed a fast Frank-Wolfe (FW) algorithm to generate feasible solutions of sub-problems. Solution candidates are computed by solving NLP problems with fixed integer variables regarding a solution pool of a MIP projection master problem.
Typical features of this approach for solving energy system MINLP models are: (i) no global branch-and-bound tree is used, (ii) sub-problems can be solved in parallel to generate columns, which do not have to be optimal, nor become available at the same time to synchronize the solution, (iii) an arbitrary solver can be used to generate solutions of sub-models, (iv) the approach (and the implementation) is generic and can be used to solve other nonconvex MINLP models, (v) the process generates feasible, not necessarily optimal, solutions during the algorithm run, which may be inspected by the decision maker.
Notice that this way of working provides a generic procedure where the sub-problems may have a black-box structure. In extremis, one may modify the optimization model during the solution process. The generated columns and solutions can also be used for performing a warm-start if the model has been changed slightly, e.g. in a dynamic optimization model with a moving horizon.
Experiments with instances of several hundreds and thousands of variables show that standard software like BARON is faster for smaller problems and reaches good solutions. We should keep in mind that, in the current Python-based implementation, we do not solve the sub-problems in parallel. The interesting notice is that for the largest DESS models we obtained significantly better dual bounds and slightly better primal bounds than BARON. Hence, for larger models the presented decomposition approach is indeed able to reach better solutions for problems with thousands of variables the energy system models usually consist of.

Appendix A Model parameters and equipment constraints
For the sake of completeness, the equations of the employed technical equipment models are given below. These include the non-linear performance maps in Equations (44)-(45) and the investment cost functions in Equations (46)-(52). The applied technical parameters are given in Table 4. Moreover, Table 5 contains the economic parameters of the DESS model. All parameter values and equations are obtained from Goderbauer et al. (2016).
Part-load performance equations: (44)-(48) Boiler (∀s ∈ B) :   as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.