Bilevel optimization for joint scheduling of production and energy systems

Energy-intensive production sites are often supplied with energy by on-site energy systems. Commonly, the scheduling of the systems is performed sequentially, starting with the scheduling of the production system. Often, the on-site energy system is operated by a different company than the production system. In consequence, the production and the energy system schedule their operation towards misaligned objectives leading in general to suboptimal schedules for both systems. To reflect the independent optimization with misaligned objectives, the scheduling problem of the production system can be formulated as a bilevel problem. We formulate the bilevel problem with mixed-integer decision variables in the upper and the lower level, and propose an algorithm to solve this bilevel problem based on the deterministic and global algorithm by Djelassi, Glass and Mitsos (J Glob Optim 75:341–392, 2019. https://doi.org/10.1007/s10898-019-00764-3) for bilevel problems with coupling equality constraints. The algorithm works by discretizing the independent lower-level variables. In the scheduling problem considered herein, the only coupling equality constraints are energy balances in the lower level. Since an intuitive distinction is missing between dependent and independent variables, we specialize the algorithm and add a procedure to identify independent variables to be discretized. Thereby, we preserve convergence guarantees. The performance of the algorithm is demonstrated in two case studies. In the case studies, the production system favors different technologies for the energy supply than the energy system. By solving the bilevel problem, the production system identifies an energy demand, which leads to minimal cost. Additionally, we demonstrate the benefits of solving the bilevel problem instead of solving the common integrated or sequential problem.


Variables
Size of batch started on production unit ( Big-M to fulfill energy balance in discretization point OC stor s Cost of storing product p CHP,prod Subsidies for on-site consumed electricity from combined-heatand-power engines p CHP,sell Subsidies for electricity sold to the grid from combined-heat-andpower engines p el,buy Electricity price for purchase p el ES∕grid Electricity price paid by the production system to the energy system/grid p el,sell Electricity price for selling p heat Price for heat p gas Price for gas Δt Length of time steṗ V min/max u,e Minimal/maximal load of energy conversion unit u u Optimality tolerance Sets and elements u ∈ d free k,f Free energy conversion units in discretization point k e ∈ E Energy forms i ∈ I Tasks j ∈ J Production equipment k ∈ K Discretization points s ∈ S Products t ∈ T Time steps u ∈ U Energy conversion units (u, e 1 , e 2 ) ∈ C Energy conversion units connecting energy forms u ∈ D max/min

Introduction
Scheduling of production systems is an active research field and important for almost all kinds of industries (Harjunkoski et al. 2014;Castro et al. 2018). Recently, research in scheduling increasingly focuses on energy aspects (Gahm et al. 2016), where new challenges arise from more volatile electricity prices and more energy supply from renewable energy sources (Merkert et al. 2015;Mitsos et al. 2018).
Large energy-intensive production systems are often supplied by on-site energy systems. Still, the on-site energy system is commonly not considered during the scheduling of the production system. Scheduling of production and energy systems is commonly performed sequentially: First, the production system is scheduled, and, second, the energy system is scheduled for the fixed energy demand. The fixed energy demand during energy system optimization can result in unnecessarily high energy cost because the production system does not account for the performance of the energy system. Integrated scheduling of production and energy system can reduce total cost. Agha et al. (2010) consider integrated scheduling of production and energy systems for operation. Leenders et al. (2020) extend the integrated scheduling of production and energy systems by the provision of control reserve, whereas Leenders et al. (2019b) also consider the synthesis of both systems. Thereby, the integrated scheduling assumes that the production system has full control over the energy system as well as complete information of the energy systems optimization problem.
However, in practice, systems are often operated by different companies. Thus, distributed optimization is performed. For this purpose, Allman and Zhang (2020) propose a solution algorithm to coordinate the optimization between different companies, i.e., an industrial process and its customers. The proposed algorithm considers that the different systems do not share their internal optimization problems with the other systems. The algorithm uses the method of Alternating Direction Method of Multipliers for distributed optimization (Boyd et al. 2010). Wenzel et al. (2020) propose an algorithm to coordinate between different production plants in a large production site. Shared resources couple the production plants. The algorithm iteratively coordinates between the systems to meet the network balance of shared resources. The algorithm in Wenzel et al. (2020) is based on extensive earlier work of the authors, e.g., in Wenzel et al. (2016) an algorithm is proposed to adjust the prices for shared resources. Maxeiner and Engell (2020) compare different dual-based methods for the distributed optimization and apply the methods to a case study with coupled semi-batch reactors. The reactors are coupled by a maximum for the combined feed flow rate.
Distributed optimization is important to take confidential issues into account while optimizing large industrial sites. In practice, each system often optimizes to its own objective and the decisions are taken sequentially. In this case, the overall problem cannot be considered as a single-level problem but should rather be considered as a hierarchical problem. This class of problem can be modeled as a Stackelberg game where the leader decides first, and subsequently, the follower decides according to the leader's decision. Here, we consider a production site where the production system (leader) decides first by optimizing its objective, and the energy system (follower) satisfies the resulting energy demand while optimizing its own objective.
In a prior contribution (Leenders et al. 2019a(Leenders et al. , 2021, we proposed a method to coordinate the scheduling of production and energy systems with different objectives without requiring the subsystems to share their internal optimization problems. Since only incomplete information is shared, the method does not guarantee to identify the optimal solution in terms of the leader's objective. Also other authors solved Stackelberg games for energy systems with exchange of incomplete information, e.g., between a utility company and customers (Maharjan et al. 2013;Soliman and Leon-Garcia 2014;Yu and Hong 2016). In contrast to our previous work, in this work, we let the production system consider the energy system's response with complete information. Thus, we consider complete information as in the integrated optimization but the production system does not have the full control over the energy system. Hence, we model the optimization problem as a bilevel problem with mixed-integer decision variables in the upper and lower level. The bilevel problem enables the identification of an optimal solution point for the production system in the given Stackelberg game situation. With the bilevel problem, we consider the sequential decision making between energy and production systems which is the common way to schedule their operation (Agha et al. 2010;Zhao et al. 2014;Zulkafli and Kopanos 2016). As an alternative, one can think of setting up a contract for demand-side management between the systems. Depending on the design of the demand-side management contract, an integrated optimization or a bilevel optimization should be applied. If the demand-side management contract is designed so that the energy system cannot decide anything anymore on its own, an integrated optimization will provide the best solution. If the demand-side management contract allows for independent decisions of the energy system, the demand-side management contract will result in additional constraints in the bilevel problem. Thus, a Stackelberg game remains, and the bilevel problem would provide the optimal solution for the production system.
In bilevel problems, misaligned optimization problems are considered by an upper level (leader) and a lower level (follower). This hierarchical structure is realized by embedding the lower-level optimization problem in the constraints of the upper-level optimization problem.
To the best of our knowledge, this work is the first to model the scheduling of production and energy systems as a bilevel problem. In the bilevel problem, the production system is considered as the upper level, and the energy system is the lower level. Here, we can model a bilevel problem since the production system has complete information on the energy system. The assignment of production system to the upper level and energy system to the lower level is chosen, since we assume that the production system announces its schedule first, and subsequently, the energy system fulfills the energy demand resulting from the production schedule. This is the same order of decision making as in the sequential optimization, which is the most used in industry (Pablos et al. 2021). Compared to the sequential optimization, the production system considers the sequential decision making in the bilevel optimization and benefits from this knowledge.
In literature, others looked at scheduling based on a bilevel problem. Avraamidou and Pistikopoulos (2019) proposed an algorithm based on multi-parametric programming to solve a MILBP. In their case study, the authors design and schedule a process with two production stages for three products. Yokoyama et al. (2019) solve a MILBP between a central power energy system and a distributed co-generation system. The problem is solved by an algorithm that uses a Karush-Kuhn-Tucker reformulation. Kostarelou and Kozanidis (2021) propose a heuristic solution algorithm to solve a MILBP for the optimal price-bidding in day-ahead electricity markets. In the MILBP, the upper level is the electricity producer and the lower level is a electricity consumer. The heuristic solution algorithm is compared to an exact solution algorithm and shows less computational requirements. Recently, Wogrin et al. (2020) present a review on applications of bilevel optimization in energy and electricity markets. As a challenge for further research, the authors identify lowerlevel problems including unit-commitment constraints with binary variables.
The production system considered in this paper results in a MILBP with coupling equality constraints, where the lower-level problem has unit commitment constraints with binary variables, i.e., binary variables if a unit is idle or running and binary variables for modeling part-load performance. For its solution, we build on our previous work ) that provides an extension of the methods in Mitsos et al. (2008) and Mitsos (2010). The extension allows the presence of coupling equality constraints. We use our previously proposed algorithm since this algorithm precedes other algorithms, which allowed us to collect lots of experience with this algorithm. From this experience and the performance in applications of industrialscale problems, we expect a good performance of the algorithm. Furthermore, the original algorithm can solve bilevel problems involving non-linearities. In the algorithm proposed in this paper, we could model non-linearities in the upper-level problem.
The algorithms from our previous work (Mitsos et al. 2008;Mitsos 2010) solve bilevel programs by successively approximating the response of the lower-level program by a discretization of the lower-level feasible set. In , we cope with the presence of coupling equality constraints by separating the lower-level variables into independent and dependent variables. Then, the independent variables are discretized as per the original algorithms, while the dependent variables are represented as variables which are fixed by the coupling equality constraints.
On the one hand, the bilevel problem of production system scheduling with an energy system as follower is a special case to the one considered in , which allows us to skip certain subproblems of the algorithm. On the other hand, in  we assume a known set of dependent lower-level variables, which are fixed uniquely by the coupling equality constraints. In our case, this assumption is violated in the sense that there exists no fixed set of dependent variables but rather a known number of dependent variables. For this purpose, we propose a procedure that determines the set of dependent variables. The choice of the dependent variable changes based on the energy system's optimal solution. The identification of the dependent variables is part of an augmented algorithm. This paper is organized as follows. In Sect. 2, we present the bilevel problem to be solved, i.e., scheduling of production and energy system. In Sect. 3, the algorithm Bilevel optimization for joint scheduling of production and… is presented as well as the optimization problems used in the algorithm. In Sect. 4, we apply the algorithm to two case studies. The case studies are based on the production scheduling example from Kondili et al. (1993) and the energy system model from Voll et al. (2013) and are extended to a bilevel problem. With the case studies, we present the performance of the algorithm, and demonstrate the merits to solve bilevel problems instead of integrated and sequential optimization problems. For comparison, we also apply the method from Leenders et al. (2019a) for incomplete information exchange to the case studies. In Sect. 5, conclusions are presented.

Bilevel problem for scheduling production and energy systems
In this section, first, we introduce the bilevel problem to be solved, i.e., scheduling of production and energy systems. Second, we analyze the properties of the lowerlevel problem. With the properties, we devise a method for generating discretization points.

Formulation of the bilevel problem
The bilevel problem considers scheduling the production system in the upper level and scheduling the energy system in the lower level. The decision variables of the upper level, i.e., production systems, are the continuous variables ∈ ℝ n x and the binary variables ∈ {0, 1} n y . includes the energy demand Ė demand t,e in each time step t and for each energy form e. We denote the joint upper-level variables as ∶= ( , ) . The decision variables of the lower level, i.e., energy system, have the same notation as it is common in energy system modeling and are the power supplied by the energy conversion units ̇ � ∈ ℝ n V and the operational state of the energy conversion units � ∈ {0, 1} n o . As it is the common notation in bilevel modeling, the decision variables of the lower level are written with a prime symbol ′ . For the representation of the lower-level variables in the upper-level problem, the prime symbol ′ is omitted. The mixed-integer linear bilevel problem (MILBP) considered in this paper has the following form: The upper-level optimization problem (scheduling of the production system) is presented in Eqs.
(1) and (2) and the lower-level optimization problem (scheduling of energy system) in Eqs.
(3)-(6). The upper-level objective f u ( ,̇ , ) considers production cost PS T and energy cost to be paid by the production system U,ES,V Ṫ + U,ES,o T . The vectors PS T , U,ES,V T and U,ES,o T describe the cost depending on the upper-and lower-level variables. are the upper-level variables, and thus, the variables of the production system. ̇ and represent the lower-level variables in the upper level. Constraints of the upper level consider the production recipe, task duration, etc. Here, we model the production system as a mixed-integer linear problem and summarize the constraints in matrix form with PS and PS as parameters of the production system. In the case study, we model the production system as a batch production system and use the State-Task-Network formulation (Kondili et al. 1993). However, also other models of the production system could be used in the upper level.
The lower-level objective f l (̇ � , � ) considers cost of the energy system L,ES,V Ṫ ′ + L,ES,o T � to supply the production system with energy. L,ES,V T and L,ES,o T are vectors for the cost depending on the lower-level variables. The lower-level constraints consider energy balances, part-load performance, the operational limits of energy conversion units etc. ̇ ′ is a vector of the power V ′ t,u,e of energy form e provided by energy conversion unit u in time step t. ′ is a vector of the operational state o ′ t,u of energy conversion unit u in time step t. o ′ t,u is a binary variable and equals 1 if energy conversion unit u is operated in time step t. ̇ ′ and ′ are the decision variables of the energy system.
The coupling equality constraints between lower and upper level, i.e., energy balances, are stated separately (Eq. (4)). In the coupling equality constraints, variables of the upper level are the energy demands of the production system Ė demand t,e for each energy form e in time step t. To clarify that this is an upper level decision variable, we also write Ė demand In general, the behavior of energy conversion units is nonlinear, e.g., by a variable efficiency (Goderbauer et al. 2016). As it is common in energy system modeling, we linearize this nonlinear unit behavior by piecewise-affine linear functions (Voll et al. 2013). For an easier separation into dependent and independent variables in the proposed algorithm (c.f. Sect. 3.2), we model each segment of the piecewiseaffine linear function as a separate energy conversion unit. Their operational states are connected via constraints. In the case study, we use the energy system model from Voll et al. (2013). However, also other energy system models could be used in the lower-level problem. Still, the part-load performance needs to be modeled by piecewise-affine functions, and all equations need to be linear as we assume this in our proof (Sect. 2.2 and Appendix B).
, ∀t ∈ T, (u, e 1 , e 2 ) ∈ C Here, all inequality constraints of the lower level are reformulated as parametric upper and lower bounds of the entries of ̇ ′ . Thereby, V − t,u,e (o � t,u ) and V + t,u,e (o � t,u ) are affine linear functions (Eq. (5)). In Eq. (6), the connection between two energy forms is modeled. These connections appear only for specific energy conversion units, e.g., combined-heat-and-power engines which supply electricity and heat simultaneously. These energy conversion units u with the corresponding two energy forms e 1 and e 2 are summarized in the set C. a u and b u are parameters to describe the connection between the energy forms.

Properties of the lower-level problem
The purpose of the energy system model (lower-level problem) is to identify optimal schedules with a medium time resolution, e.g., one hour. Thus, we neglect dynamics such as ramping constraints. Furthermore, we model an energy system without storage. Hence, no constraints to connect time steps need to be considered within the energy system model. Still, the algorithm is suitable to consider storage units, ramping constraints, and other constraints that connect time steps. For this purpose, the discretization points introduced later need to be defined for a specific energydemand curve and not for single energy demands. We expect that this extension would result in increased computational times.
In this paper, the lower-level decisions at one time step are independent of the decisions at all other time steps. Therefore and by linearity of the lower-level objective, the lower-level problem can be decomposed in a lower-level problem for each time step. The optimization problem of the lower level for a single time step t is: The decision variables of the lower level are the output power of the energy conversion units ̇ t and the operational state of the energy conversion units t . In Sect. 3, we present the algorithm to solve the bilevel problem from Eqs. (1)-(6). In the algorithm, we use the property that there always exists an optimal solution of the lower-level problem where all except one energy conversion unit per energy form e are operated at one of their operational limits, which we prove in the following two Propositions.
Proposition 1 Let the lower-level problem be given for a single energy form e and one time step t as Therein, ̇ t,e is the vector of power supplied by the energy conversion units and t is the vector of the operational state of each energy conversion unit. Note that the constraints to connect the energy forms (Eq. (10)) are omitted because only one energy form is considered.
For this problem, there exists an optimal solution in which at most one energy conversion unit is not operated at one of its operating limits. More precisely, let Problem (11) be feasible and let ( * t ,̇ * t,e ) be globally optimal in Problem (11). Then, there exists a solution ( * t ,̃̇ t,e ) that is also globally optimal in Problem (11) and for which there exists at most one u ∈ U such that Proof Consider the case that ( * t ,̇ * t,e ) is given such that for zero or one u ∈ U the following equation hold: Then, the result is proven immediately with ̃̇ t,e =̇ * t,e . Now consider instead that u 1 , u 2 ∈ U, u 1 ≠ u 2 be given such that and Then, it follows from construction of Eq. (11) that all elements of the set are feasible in Eq. (11). ̇ * t,e is an optimal solution of an linear program on a facet of the feasible set. Accordingly, all points on that facet (points in M) are optimal in Eq. (11).
The point Bilevel optimization for joint scheduling of production and… lies within M and satisfies either ̃V t,u=u 1 ,e =V + t,u=u 1 ,e ( * t ) or ̃V t,u=u 2 ,e =V − t,u=u 2 ,e ( * t ) , proving the desired property. Note that in an optimal solution where two energy conversion units are not operated at their operational limits the cost factors c L,ES,V t,u,e for these energy conversion units are equal. Otherwise, a change in the supplied power by the energy conversion units could result in a lower objective function value.
Finally, if there are more than two , the above construction can be applied successively until the same result is reached. ◻ We proved that there exists a solution for which only one energy conversion unit is not operated at one of its operational limits if there is only one form of energy e. More generally, in Appendix B, we show that the number of energy conversion units not operated at one of their operational limits can be reduced to the number of energy forms or below.

Algorithm to solve the bilevel problem
In Sects. 3.1-3.3, we first present the algorithm and second, in Sect. 3.4, we present the optimization problems considered in the algorithm, i.e., lower-bounding problem, lower-level problem, and auxiliary problem.

Algorithm
The proposed algorithm is based our algorithm from . In Djelassi et al. (2019) we use discretization points from optimal lower-level solutions for representing the lower-level objective in the upper-level problem. We assume that the lower-level variables are split into sets of dependent and independent variables. These sets are treated then differently in the discretization: independent variables are fixed and dependent variables are determined by the equality constraints.
In the bilevel problem in this paper, the choices for dependent and independent variables change for each discretization point. In the previous section and in the appendix, we proved that there are at most as many energy conversion units not operated at their operational limits as there are energy forms e considered. Thus, we can always set the power of the energy conversion units not operated at their operational limits as the dependent variables. The power supply by the energy conversion units operated at their operational limits are then the independent variables. Hence, there are at most as many dependent variables in an optimal solution point of the lower-level problem as energy forms are supplied by the energy system (c.f. Sect. 2.2).
We use these properties to identify discretization points from lower-level solutions, which means we identify independent and dependent variables. Thereby, the algorithm can solve the mixed-integer linear bilevel problem of production and energy system scheduling stated in Sect. 2.1.
In the following Sect. 3.2, we introduce the discretization points used in the algorithm. Afterward, we present the steps of the algorithm in Sect. 3.3.

Discretization points
A discretization point defines the independent and dependent variables for an optimal lower-level solution. Therein, the independent variables are fixed and the dependent variables are determined by the equality constraints. Discretization points are used to represent optimal lower-level solutions in the upper-level problem.
In the bilevel problem of scheduling production and energy systems, the energy balances are the only coupling constraints between the upper-and the lower-level problem. The energy balance is given in Eq. (4) and repeated here: The energy demands Ė demand t,e are the upper-level variables in the energy balance and the powers V t,u,e provided by each energy conversion unit u are the lower-level variables. Thus, the only variable from the upper-level problem in the lower-level problem are the energy demands Ė demand t,e . Since discretization points are used to define optimal lower-level solutions in the upper level, discretization points are defined for given energy demands Ė demand t,e . Thus, discretization points indicate the power supplied by each energy conversion unit if a particular energy demand occurs.
In Sect. 2.2 and Appendix B, we proved that there exists an optimal solution in which the number of energy conversion units not operated at one of their operational limits is less or equal to the number of energy forms e considered. We identify this optimal solution in the auxiliary problem (Sect. 3.4.3). This optimal solution is only valid for a specific energy demand. However, suppose we vary this energy demand, keep the power of the energy conversion units operated at their operational limits fixed, and compensate for the variation by adjusting the power of the energy conversion units not operated at their operational limits. In that case, the solution is still feasible until an energy conversion unit exceeds its operational limits. In this range of validity, the lower-level objective function value identified by the discretization point provides a lower-level feasible solution and, therefore, an upper bound of the lower-level objective function. Still, within the range of validity of a discretization point, another discretization point can give a tighter upper bound. We call the energy conversion units not operated at their operational limits free energy conversion units. These units define the dependent variables. The energy conversion units operated at their operational limits are named discretized energy conversion units and are the independent variables. For discretized energy conversion units, the operation is fixed in the discretization point k. These energy conversion units are either idle (index-set: D zero k ) or operated at their minimal or maximal load (index-sets: D min k and D max k ). Thus, the energy V D k,t,u=d,e supplied by these discretized energy conversion units is fixed in the discretization point k: Here, V min u,e and V max u,e are the minimal and maximal load of energy conversion unit u, respectively.
For each discretization point k, time point t, and energy form e, the power of one energy conversion unit u is not fixed, and thus, is free. These free energy conversion units are summarized by the set d free k,f . The index f counts the free energy conversion units for each discretization point k.
To fulfill each energy demand, the power of the free energy conversion units outside the operational limits of the free energy conversion unit. Outside the operational limits, the upper bound on the lower-level objective function is not considered anymore (c.f. Appendix A, Eq. (22) and Eq. (23)).

Steps of the algorithm
Algorithm 1 solves the MILBP by a lower-bounding problem, lower-level problem and auxiliary problem. The lower-bounding problem is a finite (discretization-based) relaxation of the MILBP. The lower-level problem is solved for fixed upper-level variables. The lower-level optimal solution is used to evaluate the upper-level objective to generate an upper bound. This upper bound is valid in our case since there are no upper-level constraints depending on the lower-level variables. Bilevel optimality of the upper-bounding solution is assessed based on the lower bound provided by the lower-bounding problem. In an auxiliary problem, discretization points are identified from the solutions of the lower-level problem. The algorithm repeatedly solves these problems while the lower-bounding problem is tightened by adding discretization points until the optimality tolerance u between upper and lower bound is reached.
Our algorithm exploits the structure of the lower-level problem. In the lower-level problem, all constraints are written for each time step and no constraint connect the time steps. The objective function only sums the cost of each time step. Thus, the lower-level problem can be decomposed into individual optimization problems for each time step t. Each of these individual optimization problems is parameterized only by its energy demand Ė demand t,e . The solution obtained at one time step constitutes a valid discretization point for all time steps because discretization points are defined for energy demands in a single time step which can occur in multiple time steps. Thus, probably fewer iterations are needed to solve the bilevel problem. In general, the algorithm might also be suitable to consider constraints that connect time steps in the lower-level constraints. In this case, each discretization point defines the lower-level solution for a specific energy-demand curve of the whole time horizon.
The algorithm terminates in finite time because there is a finite number of possible discretization points. The number is finite because there is a maximum number of combinations of free energy conversion units and discretized energy conversion units. In the worst case, an enumeration of all the discretization points would be necessary. However, in our applications of the algorithm, we observe that the algorithm does much better in practice and only needs a few iterations.

Subproblems of the algorithm
In this section, we present the optimization problems considered in the algorithm: lower-bounding problem, lower-level problem and auxiliary problem. To make the explanation more specific, we present the algorithm for the supply of two energy forms e, i.e., heat and electricity. These energy forms are assumed to be supplied by boilers, combined-heat-and-power engines, and electricity grids.

Lower-bounding problem (LBD)
The lower-bounding problem is the upper-level optimization problem (production system) while considering the constraints of the lower-level problem (energy system) and the lower-level objective by discretization points.
In the 1st iteration and, thus, without discretization points, the lower-bounding problem is composed from Eqs. (1), (2), (4), (5) and (6). Thus, in the 1st iteration, the lower-bounding problem corresponds to the integrated scheduling, where the production system has full control over the energy system.
After the 1st iteration, the optimal solutions of lower-level problems are considered by discretization points K. A discretization point represents the optimal solution of the lower-level problem for a particular energy demand. The lower level's objective resulting from a discretization point f l,Discr. Point  In Appendix A, equations are stated to identify if the power of the free energy conversion units is within the operational limits.

Lower-level problem (LLP)
The lower-level problem is the operation optimization of the energy system for a given energy demand. The lower-level problem can be solved independently for each time step t because there are no constraints that connect time steps in the lowerlevel. The lower-level problem was already stated in Eqs. (7)-(10). Compared to the algorithm in , here we do not have to evaluate in the upperbounding problem if the lower-level solution is feasible for the upper level. The evaluation is not necessary, since there are no lower-level variables in the constraints of the upper level. Hence, we can directly use the lower-level problem solutions to evaluate the upper-level objective.

Auxiliary problem
The auxiliary problem identifies the discretization points from each lower-level solution. In discretization points, as few as possible energy conversion units are not operated at their operational limits. Thus, the objective of the auxiliary problem is to minimize the number of free energy conversion units (units that are not operated at their operational limits) while retaining the objective of the lower-level solution.
The auxiliary problem considers the constraints of the lower-level problem. Because no constraints connecting time steps are considered in the lower-level problem, an auxiliary problem can be solved for each time step. For each energy conversion unit, we define the 4 binary variables min t,u , max t,u , zero t,u and free t,u . min t,u and max t,u identify if the energy conversion unit is operated at its minimal or maximal operational limit, respectively. zero t,u identifies if the energy conversion unit is idle. free t,u identifies if the energy conversion unit is operated between its operational limits.
The objective function of the auxiliary problem is the minimization of the number of free energy conversion units in time step t.
In the constraints, the objective function value of the lower-level problem calculated with the variables of the auxiliary problem f l,aux t needs to be equal or lower than the objective function value from the lower-level problem f l,LLP t : The constraints of the lower-level problem need to hold in the auxiliary problem (Eqs. (8)-(10)). Additional constraints define the 4 binary variables min t,u , max t,u , zero t,u and free t,u . An energy conversion unit u is either operated at its operational limits, operated between its operational limits (free), or idle. Thus, for the sum of the binary variables we can state: As stated at the beginning of Sect. 3.4, we show the equations to identify the energy conversion units at their operational limits for electricity grids, boilers and combined-heat-and-power engines. Electricity grids are an example of energy conversion units without an operational limit. Boilers and combined-heat-and-power engines are examples for energy conversion units with upper and lower operational limits. Combined-heat-and-power engines are also an example for energy conversion units connecting two energy forms.
We model two electricity grids, one for the supply with electricity and one for the feed-in of electricity into the grid. For the electricity grids, the 4 binary variables are defined in the following: The electricity grids have no upper operational 1 3 Bilevel optimization for joint scheduling of production and… limits and zero is the lower operational limit of the electricity grids. Thus, max t,u and zero t,u are defined to be zero: D Grid is a set containing all electricity grids. We identify if the electricity grid is at the lower operational limit ( min t,u = 1 ) by:

M max
e=el is a sufficiently large number such that the electricity demand from the electricity grid or the electricity supply to the electricity grid never exceeds this number. Thus, a electricity grid is identified as the free energy conversion unit for electricity ( free t,u = 1) if the power V aux t,u,e=el is greater than 0. The constraints to define the 4 binary variables for boilers and combined-heatand-power engines are given in the following: heat and electricity outputs are connected for a combined-heat-and-power engine. Thus, if a combined-heat-andpower engine is operated at its operational limit in heat supply, also the operational limit in the electricity supply is reached. Consequently, we identify if a combined-heat-and-power engine is operated at its operational limit only for one energy form, i.e., heat. We identify if the power of a combined-heat-and-power engine or a boiler is at its lower operational limit by: If min t,u equals 0, the supplied power V aux t,u,e=heat needs to be lower or equal than the heat supply in maximal load V max u,e=heat . If min t,u equals 1, the supplied power V aux t,u,e=heat needs to be lower or equal than the heat supply in minimal part-load V min u,e=heat . Since the other constraints require that if a unit is operated the supplied power is greater or equal than the minimal part-load, the unit is operated at its minimal operational limit. D CHP and D Boiler are sets with all combined-heat-and-power engine and boilers, respectively.
If a combined-heat-and-power engine or a boiler is operated at its upper operational limit V max u,e=heat is identified by: A similar explanation as for Eq. (30) applies for the upper operational limit. Whether a combined-heat-and-power engine or a boiler is idle is identified by checking if the power is within the upper and lower operational limit: Thus, if the heat supply V aux t,u,e=heat is between the operational limits, the binary variables min t,u , max t,u and zero t,u are zero and Eq. (26) sets free t,u to 1. In this case, the corresponding unit is identified as a free energy conversion unit for heat.

Integration of auxiliary problem in lower-level problem
It should be mentioned that the auxiliary problem can also be integrated into the lower-level problem if the number of free energy conversion units is known. For this purpose, the following constraints should be added: Again, free t,u identifies if an energy conversion unit u is not operated at its operational limits ( free t,u =1). The sum over free t,u is constrained by the number of allowed free energy conversion units M free . For the considered energy system in this paper, the energy system can supply 2 energy forms, and the maximum number of free energy conversion units M free is 2 (c.f. Sect. 2.2).
In this paper, we use the auxiliary problem and present a more general formulation in which an identification of the number of free energy conversion units is not necessary.

Case studies
The algorithm is applied to solve the bilevel problem for scheduling production and energy system to two case studies. The production system is modeled based on the State-Task-Network formulation for batch production systems (Kondili et al. 1993). The production system in the case studies is based on the case study from Kondili et al. (1993) (Fig. 1). In both case studies, we use the same energy system model from Voll et al. (2013). The models of both systems are described in detail in our previous publication (Leenders et al. 2019a).
In the two case studies, the energy systems differ: in the first case study (CHP subsidies), we consider: 3 boilers ( 4 MW , 1.5 MW , 0.5 MW ) and 1 combined-heatand-power engine ( 3.5 MW ) which are connected to the electricity grid and the gas grid. In the second case study (Grid alternative), we consider: 3 boilers ( 4 MW , 1.5 MW , 0.5 MW ) and 2 combined-heat-and-power engines ( 1.5 MW , 1.5 MW ) which are connected to the electricity grid and the gas grid. The big-M parameters in Eqs. (22) and (29) cannot be calculated exactly. Thus, the values need to be defined big enough to properly model the reformulation as well as small enough to prevent bad numerical behavior. In the case studies, we didn't experience any bad behavior from the choice of our big-M parameters.
Bilevel optimization for joint scheduling of production and… For both case studies, we compare the bilevel optimization with three other optimization approaches: integrated optimization, sequential optimization and the method for incomplete information exchange from Leenders et al. (2019a). In the method from Leenders et al. (2019a), the energy system responds approximated demand-dependent energy cost to the production system. The production system uses this information as basis for optimization. Thereby, the production system takes the misaligned objectives into account.
We recall that the integrated optimization of production and energy system optimizes both systems to the objective of the production system (upper-level objective). For the production system, the integrated optimization is the ideal benchmark leading to the lowest cost, since, the energy system is operated in favor of the production system.
The sequential optimization is the most common optimization approach. In the sequential optimization, the production system is first minimized to the production cost. Therein, the production system has no information on the energy system. Thus, just the upper-level problem is solved. The result is the production schedule and the energy demand. Subsequently, the energy system is optimized to the objective of the energy system. The energy demand is given by the previous production system optimization.
If we consider the situation of the production system as the leader and the energy system as the follower, both the integrated and the sequential optimization yield suboptimal solutions for the production system. The solutions are suboptimal because the response of the energy system is too optimistic or not anticipated, respectively. Although the solution of the integrated optimization promises the lowest costs for the production system, the energy system would follow its objective instead of the objective of the production system. To model this behavior, we fix the solution of Fig. 1 Production system of the case studies. The production system is a batch production system and based on the example of Kondili et al. (1993). The illustration is adapted from Leenders et al. (2019b). The circles represent the states of the products. Thus, S1, S2, and S3 are raw materials, and S7 and S10 are final products. The rectangles represent the tasks. As additional information, the form of the energy demand as well as the equipment that can be used to perform the task is given. The numbers beside the arrows give the production and consumption fraction of the tasks. E: Equipment; S: State; El.: Electricity the integrated optimization for the production system (upper level) variables and optimize subsequently the energy system (lower level). The cost increase for the production system is called regret since it expresses the excess costs that arise if a cooperative energy system is wrongly assumed. In contrast, the bilevel formulation assumes an independent energy system. Thus, the solution of the bilevel formulation provides the minimal realizable cost.
The algorithm and the benchmarks are executed on an Intel(R) Xeon(R) CPU E5-1660 v3 with 3 GHz running on openSUSE Tumbleweed with Kernelversion 4.17.13. The algorithm is implemented and the subproblems are written in C++ based on an early version of libALE . The algorithm works in interaction with GAMS such that the subproblems are automatically formulated with their input data as GAMS code. The subproblems are then solved with GAMS 25.1.3 (GAMS Development 2018) using CPLEX 12.8.0.0 (IBM Corporation 2017) applying 8 threads. The sequential optimization and the integrated optimization problems are also formulated in GAMS 25.1.3 (GAMS Development 2018) and solved with CPLEX 12.8.0.0 (IBM Corporation 2017) applying 8 threads. The time limit for all optimization problems is set to 1000 s , working memory is set to 30 GB, and the allowed absolute gap is 10 −3 . The optimality tolerance of the upper-level objective function u is set to 0.01. We choose such tight tolerances, since the optimization problems are solved fast and, thereby, we can identify the exact benefits of the bilevel solution compared to the common sequential and integrated optimization.

CHP subsidies
In this case study, the energy system gets subsidies for electricity produced by the combined-heat-and-power engines, as it is the case, e.g., in Germany. The subsidies are not forwarded to the production system resulting in misaligned objectives. In Sect. 4.1.1, we describe the detailed setup and objectives. In Sect. 4.1.2, we present and discuss the results of the case study.

Formulation
In the bilevel problem for CHP subsidies, the production system is scheduled in the upper-level problem for minimizing its production and energy cost. In the lowerlevel problem, the operation of the energy system is optimized. In the lower-level problem, the binary variables are: for each energy conversion unit and time step, a binary variable to decide if the energy conversion unit is operated, and a binary variable that allows to only buy electricity from the grid or either sell electricity to the grid at the same time.
The energy system gets subsidies for electricity produced by the combined-heatand-power engines, but the subsidies are not forwarded to the production system. The amount of subsidies are the German subsidies paid for combined-heat-andpower engine with a capacity over 2 MW . The subsidies and energy prices are provided in Table 1. The prices for energy from the grid are not equal in the case studies to distinguish the case studies further.
In the following, we present the detailed formulation of the objective functions: The objective function of the production system f u is the production cost PS T and the energy cost to be paid by the production system ( U,ES,V Ṫ + U,ES,o T ): The production cost PS T consider fix cost OC fix i,j and variable cost OC var i,j for running task i on production unit j. W t,i,j is a binary variable and equals 1 only in time step t when task i started on production unit j. B t,i,j equals the batch size only in time step t when task i started on production unit j. Furthermore, variable cost OC stor s is considered for storing amount V t,s of product s in time step t.
The energy cost U,ES,V Ṫ + U,ES,o T considers cost for purchasing gas, purchasing electricity, and revenues for selling electricity. The cost for electricity results from purchasing electricity V t,u=grid buy ,e=el for a price of p el,buy . The revenues from selling electricity result from selling electricity V t,u=grid sell ,e=el for a price of p el,sell . The cost for gas results from purchasing gas U gas t for the price of p gas . The amount of consumed gas is calculated by an affine function depending on V t,u,e=heat and o t,u . For the cost to be paid to the energy system, a profit margin of the energy system could be added. Here, we assume only a profit margin for the energy system by subsidies for running combined-heat-and-power engines.
The objective function of the energy system is the cost to be paid by the energy system:

3
These cost result from purchasing gas and electricity, revenues for selling electricity, and subsidies for running combined-heat-and-power engines. The cost for gas is again calculated by the amount of gas U gas t and the price of gas p gas . The cost for electricity is calculated by the bought electricity V t,u=grid buy ,e=el and the price for electricity p el,buy . The revenues for selling electricity result from selling electricity V t,u=grid sell ,e=el for a price of p el,sell . The revenues for selling electricity are increased by the subsidies p CHP,sell for electricity produced by the combined-heat-and-power engines. Further subsidies p CHP,prod are gathered for electricity produced by the combined-heat-and-power engines and not sold to the grid.
The case study CHP subsidies considers a time horizon of 12 h. Each time step has the same length Δt of 1 h. In the last time step of the time horizon, the production system has supply 56 t of S7 and 108 t of S10. The parameters for equipment size, storage size, and operational cost of the production system are the same as in Leenders et al. (2019a). The parameters of the energy conversion units are taken from Voll et al. (2013).

Results
The optimal bilevel solution ( f u,bilvl = 4936 €) saves 3.3 % of the cost for the production system compared to the sequential optimization ( f u,seq = 5106 €, c.f. Figure 2). The method from Leenders et al. (2019a) results in f u,inc.inf . = 5077 € cost for the production system and saves 0.6 % compared to the sequential optimization. The integrated optimization results in even lower cost ( f u,int = 4637 €) than the bilevel optimization and saves 9.2 % compared to sequential optimization. However, as mentioned previously, the integrated optimization assumes full cooperation of the energy system, which is not given in this case study. If we subsequently optimize the energy system to calculate the regret, the cost for the energy system decreases, but cost increases for the production system ( f u,corr.int = 5117€). The cost increase for the production system (regret) is 9.4 % . The regret results in even higher cost for the production system than in the sequential optimization ( 0.2 % ). Thus, the bilevel optimization provides the minimal, realizable cost for the production system (Fig. 2).
The substantial cost benefits arise in the integrated optimization because the combined-heat-and-power engines are not operated (c.f. Fig. 3). Thus, no heat and electricity are supplied by the combined-heat-and-power engines and instead by the boilers and the electricity grid. Electricity from the electricity grid is cheaper for the production system, but for the energy system, the subsidies make the combined-heat-and-power engines more beneficial. Thus, the production system favors that the energy system supplies electricity from the grid and the energy system favors to supply the electricity by the combined-heat-and-power engines. In the sequential optimization, 45.9 % of the electricity is covered by the combined-heatand-power engines. In hours where the energy system can produce more electricity than demanded, the energy system sells additional electricity to the grid.
In the bilevel optimization, 20.1 % of the electricity demand is covered by combined-heat-and-power engines because the energy system significantly benefits when using combined-heat-and-power engines. The production system chooses a demand profile which results in low supply of electricity by the combined-heat-andpower engines due to the relatively high cost of electricity covered by combinedheat-and-power engines.
The subproblems of the bilevel algorithm, the method from Leenders et al. (2019a), the integrated optimization, and the sequential optimization are solved within the time limit and, thus, are solved to the predefined gap.
The sequential optimization and integrated optimization are solved each within 2 s . The method from Leenders et al. (2019a) solves the problem in 23 s . The proposed algorithm solves the bilevel optimization problem in 262 s and needs 4 iterations. The lower-level problem and the auxiliary problem are always solved well below 1 s . The solution time of the lower-bounding problem increases with each iteration with < 1 s , 12 s , 72 s and 178 s . The solution time increases since discretization points are added in each iteration. The algorithm identifies 13 discretization points in total. From the first iteration, 7 discretization points are identified, from the second iteration, 5 additional discretization points are identified, and finally in the third iteration, 1 additional discretization point is identified. With these 13 discretization points in the fourth iteration, the lower-bounding problem reaches the same result as the evaluation of the lower-level solution in the upper-level objective  Fig. 2 In the case study CHP subsidies, cost from the bilevel and the alternative optimization approaches differ. In the integrated optimization, the cost of the production system is minimized in a single-level optimization. The cost of regret is added after a subsequent optimization of the energy system. In the sequential optimization, first, the production system is optimized without considering energy cost. Subsequently, the energy system is optimized. The method from Leenders et al. (2019a) uses incomplete information and reduces the cost slightly (Incomplete inf.) and the algorithm terminates with the optimal bilevel solution. As expected, the cost resulting from the lower-bounding problem (LBD) increase in each iteration (Fig. 4), because in each iteration we add tightening constraints by the discretization points. No such trend exists for the optimal solution of the sum of the cost from the lower-bounding problem and the additional cost from evaluation of the lower-level solution in the upper-level objective. Similar to the integrated optimization, these additional cost are named regret because they also arise from a subsequent optimization of the energy system. Since the lowest cost can only be reached by the optimal solution of the bilevel problem, the cost from the lower-bounding problem and the regret are never lower than the optimal solution of the bilevel problem (Fig. 4).  Fig. 4 Total cost for the production system in the case study CHP subsidies from the lower-bounding problem (LBD) and additional cost by evaluation of the lower-level solution in the upper-level objective (regret). The cost for the production system from the lower-bounding problem (LBD) increase with each iteration

Grid alternative
In this case study, the energy system is only connected to the gas grid. Thus, the energy system provides heat and electricity, while it provides electricity only by combined-heat-and-power engines. The production system can buy additional electricity directly from the electricity grid. First, we describe the detailed setup and the objectives (Sect. 4.2.1). Second, we present the results (Sect. 4.2.2).

Formulation
In the bilevel problem of Grid alternative, the production system is scheduled in the upper level to minimize its production and energy cost. The energy system is scheduled in the lower level. In the lower-level problem, we have the binary variables as in the case study CHP subsidies (Sect. 4.1.1).
In this case study, the production system pays a predefined price for heat and electricity to the energy system. Additionally, the production system is directly connected to the grid and can purchase electricity from the grid. Thus, the energy system has to cover the heat demand but not the electricity demand. The energy prices are given in Table 2.
The objective function of the upper level (production system) is f u ( ,̇ , ) and considers production cost PS T and energy cost ( U,ES,V Ṫ + U,ES,o T ) to be paid by the production system.: The production cost PS T is the same as in CHP subsidies (Eq. (35)). The energy cost U,ES,V Ṫ + U,ES,o T consider cost for purchasing heat and electricity from the energy system as well as cost for electricity purchased directly from the grid. The cost for heat is calculated by the heat demand E demand t,e=heat and the price for heat p heat . The cost for electricity from the energy system is calculated by the price for electricity from the energy system p el ES and the difference of electricity demand E demand t,e=el and electricity already bought from the grid V t,u=grid buy ,e=el . The cost for electricity from the grid is calculated by the price for electricity from the grid p el grid and the amount of electricity bought from the grid V t,u=grid buy ,e=el . (37) u=grid buy ,e=el + p el grid ⋅ V t,u=grid buy ,e=el . Table 2 Energy price of case study Grid alternative. p heat and p el ES are the prices for heat and electricity to be paid by the production system to the energy system, respectively. p el grid is the price for the production system for electricity from the grid. p gas is the gas price and p el,sell is the electricity price for the energy system

3
The objective function of the lower level (energy system) is to maximize the profit (equal to minimization of negative profit). The profit of the energy system P L,ES considers revenues from selling heat and electricity to the production system, cost for purchasing gas and revenues for selling electricity to the grid: The revenues from heat are calculated by the heat demand E demand t,e=heat and the price for heat p heat . The revenues from selling electricity to the production system are calculated by the price for electricity p el ES and the difference between electricity demand E demand t,e=el and electricity already bought from the grid V t,u=grid buy ,e=el . The cost of purchasing gas results from purchasing amount of gas U gas t to run the energy conversion units and the price for gas p gas . The amount of consumed gas is calculated by an affine function depending on V t,u,e=heat and o t,u . The revenues from selling electricity result from selling amount of electricity V t,u=grid buy ,e=el for a price of p el,sell .
In the last time step of the time horizon, the production system has supply the same amount as in the case study CHP subsidies, i.e., 56 t of S7 and 108 t of S10. All other parameters are the same as in CHP subsidies.

Results
Again, the bilevel optimization provides the minimal realizable cost for the production system ( f u,bilvl = 6689€). The bilevel optimization saves 3.2 % compared to sequential optimization ( f u,seq = 6910 €, Fig. 5). In this case study, the method from Leenders et al. (2019a) reaches the same cost ( f u,inc.inf . = 6689 €) as the

Fig. 5
In the case study grid alternative, cost from the bilevel and alternative optimization approaches differ. In the integrated optimization, the cost of the production system is minimized in a single-level optimization. The cost of regret is added after a subsequent optimization of the energy system. In the sequential optimization, first, the production system is optimized without considering energy cost. Subsequently, the energy system is optimized. In this case study, the method from Leenders et al. (2019a) using incomplete information reaches the cost from the bilevel optimization (Incomplete inf.) bilevel optimization. Still, as shown in the previous case study CHP subsidies, the identification of the bilevel solution is not guaranteed. The integrated optimization ( f u,int = 5170 €) would result in 25.2 % cost savings for the production system compared to the sequential optimization. As in CHP subsidies, the integrated optimization results in the best objective for the production system, but the solution from the integrated optimization is again not applicable in practice. The cost for the integrated optimization is increased after a subsequent lower-level optimization ( f u,corr.int = 6989€). The regret is 26.3 % and consequently, the cost for the production system is even 1.1 % higher than from the sequential optimization. The cost benefits by the integrated optimization arise from energy supply by the combined-heat-and-power engines (c.f. Fig. 6). The combined-heat-andpower engines are operated such that electricity demand is totally supplied by the combined-heat-and-power engines, because, in this case study, electricity from the energy system is cheaper than electricity from the grid for the production system. In the sequential optimization and in the bilevel optimization, only 4 % and 4.8 % of electricity demand is covered by combined-heat-and-power engines respectively. This low coverage of electricity by the combined-heat-and-power engines is caused because the energy system does not prefer to use combinedheat-and-power engines to cover the electricity demand. This difference in the use of the combined-heat-and-power engines results in the large cost differences between the integrated optimization and the other approaches.
The sequential and integrated optimizations are solved within 2 s . The method from Leenders et al. (2019a) solves the problem within 30 s . The proposed algorithm solves the bilevel optimization problem in 153 s and needs 3 iterations. The lower-level problem and the auxiliary problem are always solved well below 1 s . The solution time of the lower-bounding problem increases with each iteration with < 1 s , 84 s and 68 s . The algorithm identifies 7 discretization points. Thus, fewer discretization points than time steps are needed, meaning that some discretization points are valid and used for multiple time steps. In the first iteration, 4 discretization points are identified, and in the second iteration, 3 discretization points are identified. The algorithm terminates in the third iteration and the optimal solution of the bilevel problem is identified.
As expected, the cost resulting from the lower-bounding problem increase in each iteration (Fig. 7) and the cost resulting from the evaluation of the lower-level solution in the upper-level objective does not show a trend.
The case studies show that the proposed algorithm solves the bilevel problem efficiently with only a few iterations. Furthermore, the identified bilevel solution considers the misaligned objectives of production and energy system and if complete information on the energy system is available, the bilevel solution allows for the best solution in practical application.

Conclusions
The scheduling of production systems with on-site energy systems is commonly performed sequentially. An integrated optimization of both systems to an overall objective seems beneficial, but often, in practice, both systems have conflicting objectives. Bilevel optimization problems consider conflicting objectives between decision-makers, i.e., a leader and a follower, and thus, lead to a more realistic modeling.
In this paper, we formulate a bilevel problem for production and energy system scheduling. To solve this bilevel problem, we select the relevant parts from the algorithm in  and add a procedure to identify dependent and independent variables. Thereby, we can solve the bilevel problem of production and energy system scheduling.  Fig. 7 Total cost for the production system in the case study grid alternative from the lower-bounding problem (LBD) and from the evaluation of the lower-level solution in the upper-level objective (regret). The cost for the production system from the lower-bounding problem (LBD) increase with each iteration The algorithm iteratively solves the bilevel problem, while 3 optimization problems are solved in each iteration: lower-bounding problem, lower-level problem, and auxiliary problem. In the lower-bounding problem, discretization points bound the lower-level objective function by solutions from the lower-level problem in previous iterations.
The algorithm is successfully applied to two case studies based on literature examples for scheduling a production system and an energy system. The solution of the bilevel problem reaches cost savings of 3.3 % and 3.2 % compared to the sequential optimization, and 3.5 % and 4.3 % compared to an integrated optimization including regret. Furthermore, cost savings are realized compared to a previous method from the authors based on incomplete information exchange.
Thus, on the one hand, the proposed algorithm solves the bilevel problem. On the other hand, the solution of the bilevel problem provides the best realizable solution for scheduling a production system which is supplied by an on-site energy systems.

Additional constraints of the discretization points in the lower-bounding problem
In this section, we provide the additional constraints of discretization points in the lower-bounding problem. By these additional constraints, we identify if the discretization point is valid, i.e., if the free energy conversion units are operated within their operational bounds.
For every discretization point k and time step t, the heat demand Ė demand t,e=heat of the production system needs to be covered by the heat supply V D k,t,u,e=heat of the energy conversion units: The electricity demand Ė demand t,e=el needs to be covered by the electrical supply V D k,t,u,e=el of the energy conversion units: The discretization points are defined in Sect. 3.2. In our case study, we considered two energy forms heat and electricity. In the energy system, we modeled boilers, combined-heat-and-power engines, and the electricity grid. Thus, we define operational limits for energy conversion units supplying heat and electricity. The operational limits for boilers and combined-heat-and-power engines are defined for the energy form heat. The heat and electricity supply of a combined-heat-and-power engine is connected. Thus, the operational limits of the electricity supply can be defined by the operational limit of the heat supply. To identify if the operational limits of the free energy conversion unit hold, big-M formulations are used in the lower-bounding problem. For heat-supplying energy conversion units, the lower operational limit is identified by: The binary variables lower k,t,f identify if the power of the free energy conversion unit V D k,t,u=d free k,f ,e is lower than the minimal load V min u=d free k,f ,e . The big-M parameter M lower k,t,f needs to be chosen such that the heat balance can be fulfilled for every heat demand that can be supplied by the energy system. M lower k,t,f is defined as the maximum heat M max e=heat that can be supplied by the energy system reduced by the heat already supplied by the other energy conversion units V D k,t,u=d,e=heat and the minimal load of the free energy conversion unit V min u=d free k,f ,e=heat : To identify if the upper operational limit holds for heat-supplying energy conversion units, the following equation is stated: The binary variable needs to be chosen such that the heat balance can be fulfilled for every heat demand that can be supplied by the energy system. Therefore, M upper k,t,f is defined as the minimal heat that can be supplied by the energy system increased by the heat already supplied by the other energy conversion units, the maximal heat load of the largest combined-heat-and-power engine V max e=heat and the maximal load of the free energy conversion unit V max u=d free k,f ,e=heat : The parameter M min e=heat is the minimal heat that can be supplied and is commonly 0. The maximal heat load of the largest combined-heat-and-power engine V max e=heat is required for the case a combined-heat-and-power engine unit is selected as free energy conversion unit for electricity. Thus, with the parameter M upper k,t,f , the heat power of the free energy conversion unit can be reduced such that the minimal heat supply can be reached in the heat balance. In our model, we consider two types of electricity grids: electricity purchase and electricity sale. If an electricity grid is the free energy conversion unit, we also need to check if the operational limits hold. The lower operational limit of the electricity grids is 0. Whether the lower operational limit holds is identified by: The big-M parameter M max e=el is the maximum electricity that the production system could demand from the energy system. For the lower operational limit, Eq. (45) checks if V D k,t,u=d free k,f ,e=el is non-negative. The upper operational limit is never exceeded for electricity grids:

Extended proof for multiple energy forms
In this section, we extend the proof to multiple energy forms.
Proposition 2 Let the lower-level problem be given for multiple energy forms e and one time step t as Therein, ̇ t is the vector of power supplied by the energy conversion units u in all energy forms e ∈ E and t is the vector of the operational state of each energy conversion unit u.
For this problem, there exists an optimal solution in which the number of energy conversion units not operated at one of their operating limits is less or equal to the number of energy forms e considered. More precisely, let Problem (47) be feasible and let ( * t ,̇ * t,e ) be globally optimal in Problem (47). Then there exists a solution ( * t ,̃̇ t,e ) that is also globally optimal in Problem (47) and for which there exists at most |E| elements u ∈ U such that Note that a connecting energy conversion unit possesses two power flows that are subject to operational limits. Nevertheless, each connecting energy conversion unit is only counted once for the purposes of this proposition.
Proof Consider the case that ( * t ,̇ * t,e ) is given such that for at most |E| elements of u ∈ U the following equation holds: Then, the result is proven immediately with ̃̇ t,e =̇ * t,e . In the following, we prove that if we consider two energy forms e 1 and e 2 , the maximum number of energy conversion units at their operational limits can be reduced to at maximum two in an optimal solution. This result can then be generalized to an arbitrary number of energy forms e if we iteratively consider two energy forms and reduce the number of the energy conversion units not at their operational limit to at maximum two.
For the case of the two energy forms e 1 and e 2 , energy conversion units can either output energy form e 1 or e 2 , or if they are a connecting energy conversion unit, they output both energy forms e 1 and e 2 . The energy forms e output by an energy conversion unit u is given by the set E u . Let u 1 , u 2 , u 3 ∈ U, u 1 ≠ u 2 ≠ u 3 be given such that and To prove Proposition 2, we distinguish 4 cases that can occur.
• Case 1: All three energy conversion units output a single energy form • Case 2: One energy conversion unit is a connecting energy conversion unit and the other two energy conversion units output a single energy form • Case 3: Two energy conversion units are a connecting energy conversion unit and the other energy conversion unit outputs only a single energy form • Case 4: All three energy conversion units are connecting energy conversion units In the following, we use the fact that if a connecting energy conversion unit reaches the operational limit in one energy form, also the operational limit in the other energy form is reached. Thus, the operational limits of a connecting energy conversion unit u outputting energy forms e 1 and e 2 can be written in terms of the limits in one energy form using Eq. (10): , proving the desired property. Then, only two energy conversion units are not operated at their operational limits.
Case 3 Energy conversion unit u 1 outputs energy form e 1 , u 2 and u 3 are connecting energy conversion units and output energy forms e 1 and e 2 .
From construction of Eq. (47), it follows that all elements of the set are feasible in Eq. (47). ̇ * t,e is an optimal solution of an linear program on a facet of the feasible set. Accordingly, all points on that facet (points in M) are optimal in Eq. (47). Note that again in the last equation of M the constant parts a u=u 2 and a u=u 3 are eliminated. Following the definition of Eq. (54), we define the change in energy form e 1 until an energy conversion unit reaches its operational limits: With ΔV Case3 e=e 1 , we can manipulate the supplied power by the energy conversion units such that the point ( * t ,̃̇ t,e ) with ̃V t,u,e =V * t,u,e , ∀u ∈ U ⧵ {u 1 , u 2 , u 3 }, e ∈ E u and (55)V ( * t ), ∀u ∈ {u 1 , u 2 , u 3 } ∧ a u ⋅ o * t,u = b u ⋅V t,u,e=e 1 −V t,u,e=e 2 , ∀u ∈ {u 2 , u 3 } ∧V t,u=u 1 ,e=e 1 +V t,u=u 2 ,e=e 1 +V t,u=u 3 ,e=e 1 = V * t,u=u 1 ,e=e 1 +V * t,u=u 2 ,e=e 1 +V * t,u=u 3 ,e=e 1 ∧ b u=u 2 ⋅V t,u=u 2 ,e=e 1 + b u=u 3 ⋅V t,u=u 3 ,e=e 1 = b u=u 2 ⋅V * t,u=u 2 ,e=e 1 + b u=u 3 ⋅V * t,u=u 3 ,e=e 1 } (57) t,u=u 2 ,e=e 2 =V − t,u=u 2 ,e=e 2 ( * t ) or ̃V t,u=u 3 ,e=e 2 =V + t,u=u 3 ,e=e 2 ( * t ) , proving the desired property. Note that if b u=u 2 =b u=u 3 , then u 1 does not change its power supply.
With these changes in the supplied power, only two energy conversion units are not operated at their operational limits.

Case 4
This case only occurs if 3 connecting energy conversion units are not at their operational limits and do not have the same conversion factor b u . If at least 2 of the 3 have the same conversion factor b u , the special case from Case 3 can be applied and only the two equal energy conversion units change their power supply.
Energy conversion unit u 1 , u 2 and u 3 are connecting energy conversion units and output energy forms e 1 and e 2 . From construction of Eq. (47), it follows that all elements of the set are feasible in Eq. (47). ̇ * t,e is an optimal solution of an linear program on a facet of the feasible set. Accordingly, all points on that facet (points in M) are optimal in Eq. (47). Note that again in the last equation of M the constant parts a u=u 1 , a u=u 2 and (58)V t,u=u 1 ,e=e 1 =V * t,u=u 1 ,e=e 1 + ΔV Case3 e=e 1 V t,u=u 2 ,e=e 1 =V * t,u=u 2 ,e=e 1 M = {( * t ,̇ t,e )|V t,u,e =V * t,u,e , ∀u ∈ U ⧵ {u 1 , u 2 , u 3 } ∧V − t,u,e=e 1 ( * t ) ≤V t,u,e=e 1 ≤V + t,u,e=e 1 ( * t ), ∀u ∈ {u 1 , u 2 , u 3 } ∧ a u ⋅ o * t,u = b u ⋅V t,u,e=e 1 −V t,u,e=e 2 , ∀u ∈ {u 1 , u 2 , u 3 } ∧V t,u=u 1 ,e=e 1 +V t,u=u 2 ,e=e 1 +V t,u=u 3 ,e=e 1 = V * t,u=u 1 ,e=e 1 +V * t,u=u 2 ,e=e 1 +V * t,u=u 3 ,e=e 1 ∧ b u=u 1 ⋅V t,u=u 1 ,e=e 1 + b u=u 2 ⋅V t,u=u 2 ,e=e 1 + b u=u 3 ⋅V t,u=u 3 ,e=e 1 = b u=u 1 ⋅V * t,u=u 1 ,e=e 1 + b u=u 2 ⋅V * t,u=u 2 ,e=e 1 + b u=u 3 ⋅V * t,u=u 3 ,e=e 1 } a u=u 3 are shortened. Following the definition of the Eq. (54), we define the change in energy form e 1 until an energy conversion unit reaches its operational limits: With ΔV Case4 e=e 1 we can manipulate the supplied power by the energy conversion units such that the point ( * t ,̃̇ t,e ) with ̃V t,u,e =V * t,u,e , ∀u ∈ U ⧵ {u 1 , u 2 , u 3 }, e ∈ E u and lies within M and satisfies either ̃V t,u=u 1 ,e=e 1 =V + t,u=u 1 ,e=e 1 ( * t ) , ̃V t,u=u 2 ,e=e 2 =V − t,u=u 2 ,e=e 2 ( * t ) or ̃V t,u=u 3 ,e=e 2 =V + t,u=u 3 ,e=e 2 ( * t ) , proving the desired property. With these changes in the supplied power, only two energy conversion units are not operated at their operational limits.
Finally, if there are more than three energy conversion units u ∈ U with V − t,u,e ( * t ) <V * t,u,e <V + t,u,e ( * t ) , the above constructions can be applied successively to three energy conversion units until the same result is reached. ◻ Acknowledgements Hatim Djelassi and Alexander Mitsos acknowledge the funding by Réseau de transport d'électricité (RTE, France) as part of the project Hierarchical Optimization for Worst-case Analysis of Power Grids. Ludger Leenders and André Bardow acknowledge the funding by the Swiss Federal Office of Energy's SWEET programme as part of the project PATHFNDR.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interests
The authors declare that there is no known conflict of interest regarding the publication of this paper (60)

ΔV
Case4 e=e 1 = min{V + t,u=u 1 ,e=e 1 Bilevel optimization for joint scheduling of production and… Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long 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/.