Budget-cut: introduction to a budget based cutting-plane algorithm for capacity expansion models

We present an algorithm to solve capacity extension problems that frequently occur in energy system optimization models. Such models describe a system where certain components can be installed to reduce future costs and achieve carbon reduction goals; however, the choice of these components requires the solution of a computationally expensive combinatorial problem. In our proposed algorithm, we solve a sequence of linear programs that serve to tighten a budget—the maximum amount we are willing to spend towards reducing overall costs. Our proposal finds application in the general setting where optional investment decisions provide an enhanced portfolio over the original setting that maintains feasibility. We present computational results on two model classes, and demonstrate computational savings up to 96% on certain instances.


Background
Governments worldwide are pushing towards an increasing use of renewable energy technologies. In line with emission targets set by the European Commission-as part of the 2050 Low Carbon Economy roadmap to reduce emissions to 80% below the levels in the year 1990 [6]-the German federal government plans to increase the percentage of energy derived from renewable sources to 80% by the year 2050 [3]. This rampant increase in the share of renewables requires important investment decisions that guide future energy policies, e.g., the expansion of existing energy capacity infrastructure to accommodate the needs and demands of future energy production and supply. Expanding existing capacity to include photovoltaics, hybrid generation systems, and storage devices are a few examples of such investment decisions that can potentially lead to lower greenhouse gas emissions [18].
Optimization models have a rich history for both operations of installed energy systems [1,33], as well as extending existing infrastructure; see, e.g., so-called capacity expansion models [2,25], models for integrating renewable sources with existing fossil fuels [29], and planning for expanding transmission networks [24]. Typically, mixed-integer programs (MIPs) are developed and employed to inform these decisions as well as optimize operations. Two examples of such energy system models are MARKAL [22] and TIMES [23]; see, also [20] and references within. The Framework for Integrated Energy System Assessment (FINE) is an open source python package that provides a framework for the modeling, optimization and analysis of energy systems [7,31]. The goal of such optimization models is to minimize the total annual costs for deploying and operating energy supply infrastructure, subject to the technical and operational constraints of a multi-commodity flow energy-system problem. FINE provides the functionalities to set up energy system models with multiple regions, commodities and time steps. Examples of commodities include electricity, hydrogen, and methane gas. Time steps can also be aggregated to reduce the complexity of the model [12]. In addition to existing infrastructural costs, costs are also incurred by building new components and increasing their capacities. The work in this article arises from a collaboration between mathematicians and energy-and-climate researchers at the Friedrich-Alexander-Universität Erlangen-Nürnberg and the Forschungszentrum Jülich, with the goal to improve the performance of the FINE package as well as other energy system models.
Against the backdrop of decreasing prices of renewable energy sources [13], rising CO 2 emission costs [5,6], a transforming energy demand and new options for energy storage and conversion [32], the consideration of novel technologies offers opportunities for further cost reduction of existing energy systems. We consider the problem of determining which components of an energy capacity infrastructure to install such that total annual costs are further reduced. Installing a new component incurs a fixed cost; however, overall costs are potentially reduced by allowing new components access to larger ranges of energy supplies thereby leading to a more efficient utilization of the entire system. Neumann and Brown consider a similar problem to expand transmission while minimizing total annual costs [27]. Such problems also find application in several other contexts within energy systems that seek to minimize total annual costs; see, e.g., [30] for a model that extends natural gas transmission networks. For an overview of transmission expansion planning problems, see, Mahdavi et al. [26].

Challenges
Such investment optimization problems are frequently formulated as MIPs that employ binary variables to inform "yes" or "no" decisions for utilizing new technologies; see, e.g., [9,20]. The choice of such discrete decision variables governing whether a new technology "is-built" leads to a significant increase in the computational effort to solve the MIP. Models covering multiple spatial regions impose further computational challenges as additional binary variables are required for each region. For an introduction to the challenges and model simplifications, see, e.g., [8,20]. However, as we mention in Sect. 1.1, integrating novel technologies into existing energy systems leads to potentially reduced overall costs and is also advantageous in keeping abreast with the dynamically changing energy sector. To this end, we provide a method that offers scope for reduced runtimes, both for systems with several existing technologies as well as new technologies to choose from.
This work is related to previous works on a so-called "district model" that includes six multi-family houses and households in each building [15], and also to a model with 16 transmission zones within Great Britain [28]. In previous improvements to the FINE package, Kannengießer et al. use a time-series aggregation method and present a temporally aggregated simplification of the model [15]. The authors use a multi-regional district model and fix the binary design decisions for all components. The corresponding optimization model is thus reduced to a linear program (LP); this model is solved easily, however the solutions are suboptimal. In contrast the approach we present is an iterative heuristic that is guaranteed to converge to the optimal solution.  Used unit capacity of component c on edge (l 1 , l 2 ) [kWh] In the above notation, bin denotes a binary decision variable, cap and op denote continuous decision variables for the capacity and operation, respectively. The parameters CapMin and CapMax denote bounds on the cap variable, while TAC denotes total annual cost. We provide details in Sect. 2.2.

Optimization model
The models build with the FINE package allow the inclusion of two types of components: (i) "optional" components that are modeled with additional investment costs independent of the installed unit size, and (ii) already existing components that are modeled with an installation cost contribution that is linearly dependent on the installed unit size. We have five choices for these optional and existing components that are relevant to the discussion in this work: Source, Sink, Storage, Conversion, and Transmission. The complete FINE package includes other components as well, and the user specifications decide the components that form part of the optimization model. We then have a graph where the nodes include combinations of the first four component types, while each Transmission component represents an edge connecting two nodes. Let the index c ∈ C denote the components, index l ∈ L denote the locations of a node, and tuple (l 1 , l 2 ) with l 1 , l 2 ∈ L, l 1 = l 2 denote the start and end locations for an edge, respectively.
Further, let C N ⊂ C denote the set of four components that form nodes, and C E denote the set of components that form edges; i.e., C N = {C Source , C Sink , C Storage , C Conversion } and C E = {C Transmission }. Here, the source components include nodes that provide commodities to the graph, while the sink components withdraw commodities from the graph. The storage components are nodes that connect the different time steps with each other by incorporating a so-called state of charge (SoC) variable; these components operate as a sink (that increases the SoC) or a source (that decreases the SoC) "state variable" (the SoC) at a given time step. The conversion components are nodes that connect multiple "commodity-subgraphs" to each other by converting one commodity to another; e.g., an electrolyzer consumes electricity and converts it into hydrogen. Finally, let L c denote the location(s) of a component c. Although the decision variables and parameters within the FINE package include other indices (such as time and commodity) as well, we suppress these indices as they are not relevant to the discussion in this article; see, the complete model description in the appendix, and further details in [31].
FINE includes binary variables to inform whether new optional components are built, and calls these binary variables as "IsBuilt-variables" [31]. A value of 1 denotes a component is built at a given location. If a component already exists, it does not have a IsBuilt-variable; we assume that there is at least one such component, else there is no configuration to start with. Building an optional component c at location l incurs a fixed total annual cost (TAC) of TAC bin c,l . For the sake of notation, we assume existing components also have associated binary variables that are always 1 with TAC bin c,l = 0. Further, the decision to build is governed by its corresponding capacity variable, cap c,l . The value of a capacity variable corresponds to the scale of the installed component, e.g., for a photovoltaic Source component it corresponds to the area of solar panels installed [31]. If any capacity is installed, the corresponding IsBuilt-variable takes the value 1. All optional components have minimum capacity thresholds that are informed by data; we take this threshold as 0 if data is unspecified. In other words, cap c,l is a semi-continuous variable that is either 0 or lower bounded by the minimum capacity. Analogously, maximum capacity thresholds are available as well, and we take these as +∞ if unspecified. Equations (1) and (2) summarize this discussion for the nodes and edges, respectively. (2) are analogous to equations (1). Further, we have operational variables, that include time indices, corresponding to the components that we denote by op c,l ; these denote how much of the installed capacity is actually used.
Below is the optimization model we consider in this article. The objective function in Eq. (3a) abbreviates the following quantity: Here, cap, bin, and op denote three vectors corresponding to the capacity, IsBuilt, and operational variables, respectively. The coefficients -TAC cap , TAC bin , and TAC op -denote vectors of appropriate size for the TAC corresponding to installing capacity, building new components, and operating these components, respectively. That is, the objective function includes cost contributions that scale with the size of the installed capacities of the components linearly (TAC cap ), cost contributions that are independent of the size of the installed capacity but occur if a component is built (TAC bin ), as well as cost contributions that are connected to the operation of the built components (TAC op ); see, the appendix for details. Then, the objective function seeks to minimize the sum of the TACs of the entire system, with both optional and existing components. The bounding constraints in Eq. (3c) and (3d) enforce further limits on the capacity and operation variables with and without time indices, respectively. The component-linking-constraints of Eq. (3e) include commodity balances, annual commodity inflow and outflow limits as well as shared potentials; they also serve to define limits via the capMax parameter.
The bidirectional and symmetric nature of the components along edges, as well as constraints that model the optimal power flow using the standard linearized DC formulation, is expressed via the transmission constraints of Eq. (3f). The storage constraints of Eq. (3g) express the charging and discharging status via the state of the charge for the different components. For a detailed description of these constraints, see [31]; we provide a complete model description in the appendix. Importantly, constraints (3c)-(3g) do not contain any binary variables.

Background
MIP solvers typically rely on branch-and-bound strategies to identify the optimal solution. The presence of additional binary variables results in a larger search tree that can create potential computational challenges. In this section, we propose an algorithm that prunes branches that lead to suboptimal solutions. To this end, we derive a sequence of cuts by solving a sequence of LPs. We use these cuts to determine a "budget" that provides a valid inequality for the original problem-that is now reduced in size. We formalize this concept in the Budget-Cut Algorithm, present it in Fig. 1, and explain it in more detail below. We first note that when constraints (1d) and (2d) are reduced to their continuous relaxation, model (3) is an LP. To this end, we solve two models where all the binary IsBuilt-variables are apriori determined. Consider the following optimization models:z and We denote models (5) and (6) as Existing and Extended with optimal objective function values ofz and z, respectively. We assume Existing is feasible; in Sect. 5 we discuss the implications of this assumption. The Existing model determines a baseline where no optional components are built, and the only components used for the solution are those that that do not have costs independent of the size of the installed capacity; i.e., all the IsBuilt-variables are set to 0. Thus, the second term -TAC bin · bin -in the objective function of model (3) is 0 for all the IsBuilt variables. We note that the Existing model still includes the already existing components, and these are the only components we optimize over. The Extended model determines the other extreme where all the optional components are built; i.e., all the IsBuiltvariables are set to 1. The idea of the Extended model is subtle. Although all the IsBuilt-variables are set to 1, we ignore the cost to build them by not including the second term -TAC bin · bin -of the objective function of model (3). By ignoring the cost to build the components, the Extended model focuses on finding the most cost-effective components to use and install capacities on. Further, we do not include the lower thresholds on capacity for optional components in Eqs. (1c) and (2c). Since optional components are built to reduce total costs, the Extended model represents the best we can hope to achieve. Intuitively, the Existing model chooses an optimal solution from the existing capacity infrastructure, while the Extended model chooses the optimal solution from the existing and the optional capacity infrastructure. The following lemma relates the Existing and Extended models with the "true" model (3).
Proof The feasible region represented by constraints (5b)-(5d) is a subset of the feasible region represented by constraints (3b)-(3g). Further, the objective functions of models (3) and (5) are identical, since TAC bin · bin is 0 for model (5). With a minimization objective, z * ≤z follows. Next, we note that constraints (6b) and (6d) restrict the capacity variables between 0 and capMax; while, constraints (3b) enforce the capacity variables are either 0 or bounded between capMin and capMax. Since the other constraints of models (3) and (6) are identical, the feasible region of model (3) is smaller than that of model (6). Further, the objective function in Eq. (3a) is at least that in Eq. (6a). Thus, z ≤ z * follows.  (6), respectively. Here, bin = 0 and bin = 1. Further, we definez − z = b ≥ 0. Next, we show that b determines a "budget"the maximum cost we are willing to spend on building IsBuilt components. In other words, spending more than b exceeds the potential savings offered by the addition of optional components. We use the following corollary to Lemma 1 in the proof.
Proof From the proof of Lemma 1, it follows that the solution (cap * , op * , 1) is feasible for model (6).

Lemma 2 TAC bin · bin ≤ b is a valid inequality for model (3).
Proof From Lemma 1 and the definition of b we have, TAC cap · cap * + TAC op · op * + TAC bin · bin * − (TAC cap · cap + TAC op · op) ≤ b. Then, from Corollary 1, we have TAC bin bin * ≤ b, for any optimal solution for model (5).
For sake of completeness, we rewrite the valid inequality in its full form: as well as model (3) with the valid inequality: (8c)

A budget-cut algorithm
Equation (7) provides a valid inequality for model (3) that takes the form of a typical 0-1 knapsack constraint. For ease of exposition within this section, we represent Eq. (7) as i∈I a i · bin i ≤ b. Although the weights of the knapsack items, a i , are known, the benefit derived by the addition or removal of a single element requires the solution of a new instance of model (3). We could determine this benefit by solving an instance with a subset of IsBuilt-variables fixed to 1, and then compute the possible savings minus the construction costs. However, this requires a solution to 2 |I | problem instances. Instead, we use b to determine elements that are too expensive to construct, remove them apriori via Lemma 2, and re-solve model (3) with the valid inequality (i.e., model (8)). Small values of b provide tighter models. In this section, we seek to further reduce the size of the reduced model (8) by reducing b. We begin by distinguishing three cases.
1. Case 1: b < min i a i It follows from Lemma 2 that bin i = 0, ∀i ∈ I ; then, from Lemma 1 z = z * . 2. Case 2: min i a i ≤ b ≤ max i a i We proceed by first fixing all binaries with a i > b to 0. Then, we recompute the budget; i.e., we solve Extended with these binaries set to 0. This guarantees the updated budget is no more than the previous budget, and we repeat this process. 3. Case 3: max i a i < b In this case, we cannot reduce the budget further. Further, if i∈I a i < b is true, Eq. (7) is redundant. We reflect these three cases in the Budget-Cut Algorithm. The algorithm takes as input an instance of model (3), and a time limit, T I M E. A key prerequisite of the algorithm is that model (5) is feasible, else the algorithm's step 3 fails. In Sect. 5, we provide a discussion on handling infeasible instances. The other assumption of the algorithm is that not all the components are optional. In the absence of this assumption, the initial configuration has no cost and thus there is nothing to do. The Budget-Cut  (3), or the corresponding optimality gap and the best found feasible solution in the time limit. The Budget-Cut Algorithm provides at least three advantages compared to a naive solution method. First, from Lemma 2 we can directly proceed to Step 23; this happens if the algorithm recognizes that the budget does not allow building any optional elements. Second, if the algorithm does proceed to the while loop, we are guaranteed at least one element has a binary variable that is fixed to 0. This ensures a finite termination of the algorithm in at most |I | iterations. Third, the solutions of the LPs serve as warmstarts for model (8). Figure 1 presents a visualization of the budget update. In the next section, we compare the computational performance of the algorithm with a naive solution method.

Setup
To examine the computational performance of the Budget-Cut Algorithm, we conduct a number of computational experiments on different instances of model (3). The FINE package uses a time series aggregation method to reduce the size of the optimization model; i.e., it aggregates the complete considered time horizon of, e.g., 365 days into so-called "typical days". For details on the time series aggregation methods used within FINE, see [12]. We define a model instance with a time horizon of one year using typical days -7, 14, 28, 56, and 112 -and weather years from 1995 to 2000.
The weather year parameter does not affect the size of the optimization model, but determines which input data set of commodity power demands and supplies is used. However, models with a larger number of typical days result in larger optimization models.
All tests in this article are carried out with Pyomo 5.7.1 [11] using Gurobi 9.0.2 [10] on two machines: (i) an Intel Core i7 2.8 GHz processor with 16 GB of memory, Fig. 2 Structure of the Self-sufficient building scenario of Sect. 4.1 from Figure 1 of [19] and (ii) a node on the Jülich Research on Exascale Cluster Architectures (JURECA) supercomputer, with a cluster's batch partition 2x Intel Xeon E5-2680 v3 (Haswell) and a 2.5GHz processor and 128 GB of memory [14]. We refer to these two machines as Machine I and Machine II, respectively. We solve smaller models of up to 56 typical days on Machine I with the Gurobi threads parameter set to 3, and larger models with 112 typical days on Machine II with the threads parameter set to 32. We use T I M E = 900 seconds and T I M E = 15, 000 seconds as time limit for our computational experiments on Machine I and Machine II, respectively.
We use the self-sufficient building scenario (SelfScen) of Kotzur et al. [19] for our experiments; see, also [16]. The SelfScen optimizes the cost of a single household building to construct and operate on its own, thereby being self-sufficient from an energy perspective. The available technologies include rooftop photovoltaic systems and batteries for short-term electricity storage, as well as reversible fuel cells and liquid organic hydrogen storage systems for long-term energy storage, to meet demand for power. There is also a demand for heat, this is fulfilled by a combination of electric boilers, heat pumps and heat storage. SelfScen includes the following commoditiesheat, electricity, hydrogen, liquid organic hydrogen carrier (LOHC), and high temperature heat. Scenario instances include heat and electricity demand for the household, and the maximum power rate that can be generated by the photovoltaic units, for the weather years 1995-2000. To this end, the SelfScen chooses the technologies to install (i.e., bin), the capacities for each component (i.e., cap), and their operation (i.e., op).
Optional modeled components are the heat pump, the reversible fuel cells, and conversion technologies that are required for using the hydrogen storage components. See Fig. 2 for an illustration of the SelfScen, and [17] for the dataset associated with the SelfScen. In Online Appendix A.3 we provide additional computational experiments for another scenario that includes multiple buildings.

Computational experiments
In Table 1, we compare the computational results for the naive solution method and the Budget-Cut Algorithm. Within our time limit, the Budget-Cut Algorithm succeeds in finding the optimal solution in all the instances. The naive solution method, however, fails to find even a feasible solution for all instances with 56 typical days. For the instance with 112 typical days and the year 1995-that we solve on Machine II-the naive solution terminates with a large MIP gap of 45%. The value in the third column of Table 1 is the best known value of the objective function in model (3); the value reported by the algorithm is indeed the optimal in all instances. Next, we note that the naive solution method performs better for smaller instances up to 7 typical days. However, for the larger instances on Machine I the improvements are significant; for instances with 14, 28, and 56 typical days the average improvement is 61.2%, 88.0%, and 82.2%, respectively. For the largest instance with 112 typical days that we solve on Machine II, the runtimes are an order of magnitude lower except for the year 2000; the average improvement here is 81.4%.
All instances in Table 1 are solved without any iterations of the Budget-Cut Algorithm; i.e., the while loop in the Budget-Cut Algorithm is not entered. Then, the valid inequality in Eq. (7) is trivially true with bin i = 1, ∀i ∈ I . In other words, no component is trivially excluded for having too high costs. However, even then the naive solution method is significantly slower as the algorithm benefits from warmstarts derived from the Extended model as well as additional cuts derived by the optimization solver with the addition of the trivially true valid inequality.
To this end, we modify the SelfScen to ensure i∈I a i · bin i > b; thus, the algorithm enters the while loop. For details on how we modify the SelfScen, see Online Appendix A.2. We denote this modified scenario as ModSelfScen. Table 2 presents our computational results for the ModSelfScen, analogous to Table 1. The optimal objective function values for the ModSelfScen are larger than those of SelfScen due to the increased objective function coefficients of the ModSelfScen. All instances in Table 2 are solved to optimality within the time limit, thus we remove the two columns corresponding to "Gap" in Table 1. All instances require two iterations. Trends for the ModSelfScen mirror those of the SelfScen, however the percentage savings are smaller for the former. The average improvements for instances with 14, 28, 56, and 112 typical days are 25.6%, 53.2%, 61.2% and 68.4%, respectively.
In Online Appendix A.3 we provide an additional set of computational results on another class of scenarios. The results again follow the same trends as those we report above. The "Gap" columns provide the MIP gap reported by Gurobi; a blank denotes a gap of 0%, while ∞ denotes no feasible solution is found. Entries marked with T in the "Time" column denote the maximum time limit is reached. The "Improvement" columns provide the percent improvement in run time for the algorithm compared to the naive solution method. All entries are rounded to the first decimal. For details, see Sect. 4

Limitations
At least two limitations of our proposed algorithm offer work for future research. First, as we mention in Sect. 3, the algorithm fails when the Existing model is infeasible. However, checking the usability of Budget-Cut Algorithm for an optimization model requires the solution of a single LP; computationally this is a small effort. Future work could focus on determining feasible start solutions for Budget-Cut Algorithm, as well as solving the Existing and Extended problems in parallel. Next, given a feasible solution to initiate the algorithm, we can determine a valid upper bound for model (3). That being said, determining feasible solutions for a MIP is, in general, hard [4]. However, for certain classes of problems-such as the one we present-feasibility is often maintained when no investment decisions are made [21]. A second limitation of Budget-Cut Algorithm occurs when the budget is too large; i.e., b > a i . In this case, the algorithm skips directly to the final solution of the MIP without a valid inequality. The only benefits of our proposal in this case is the use of a feasible solution from Extended as a warmstart. However, even then, the computational benefits we demonstrate in Sect. 4 with the SelfScen instances are significant.
Finally, we mention that several practical problems include a known budget-the maximum possible expenditure for building new components. Then, the parameter b of Algorithm Budget-Cut Algorithm serves as an input.

Summary
To summarize, we present a simple-to-implement algorithm for reducing runtimes of a capacity extension problem. This problem is computationally demanding when exponentially many choices for installing new components exist. Intuitively, the general class of problems we study addresses the following concern: given a portfolio of potential investments with varying purchase and operation costs, choose the ones that minimize long-term horizon expenditures subject to a given budget. Such situations also find application in the general setting when investment decisions are optional and only provide an enhanced portfolio, thereby the original problem serves as a basecase. We relate this problem to a knapsack problem, and propose an algorithm that determines a valid inequality to cut off suboptimal branches of the branch-and-bound search tree. Our algorithm rests on determining upper and lower bounds by solving two extremes of problems -the first where no optional components are built, and the second where all optional components are built, respectively. By iteratively pruning off suboptimal solutions, we increase the lower bounds obtained by the second of these extremes.
Documentation and data for the FINE code is available under the MIT License at https://github.com/FZJ-IEK3-VSA/FINE. Data for the SelfScen is available under the MIT License at https://data.mendeley.com/datasets/zhwkrc6k93/1 [16].