Mixed-Integer Linear Programming Formulation of Combined Heat and Power Units for the Unit Commitment Problem

The energy system in Europe changes towards a smart energy system with interconnection of heat, power and gas networks, increasing the flexibility and efficiency. Combined heat and power units play an important role in such systems, since they have a high total efficiency and enable the interconnection between the networks. Yet, control concepts are needed to ensure an optimal operation planning of available capacities. The resulting unit commitment problem can be approached with mixed-integer linear programming. This paper presents a generic modelling method of coupled and decoupled combined heat and power units with mixed-integer linear programming. Furthermore, operation data of two steam turbine types are compared with the model performance to demonstrate the feasibility of this modelling approach. Finally, a model of a combined cycle is presented, derived from the presented method.


INTRODUCTION
New challenges arise with the increasing share of renewable energy in the energy mix due to EU's climate target of 80% reduction in annual green house gas emissions in 2050 compared to 1990 [1].This renewable energy consists mostly of wind and solar power and is therefore of volatile nature, while energy balances in the energy grids have to be satisfied.Especially the electricity generation has to meet the demand at all times due to the incapability of storing electric energy in the grid.Nevertheless, the full integration of renewable energy is desired, which leads currently to positive and negative residual energy that has to be provided mostly by flexible fossil-fuelled power plants.
According to Mathiesen et al. [2,3], the integration of these increasing, fluctuating renewable energy sources, with a target of even 100% renewable energy in the system, should be done by a smart energy systems approach, interconnecting electricity, heat and gas networks and increasing economic benefits and flexibility through conversion between different energy forms.The interaction between electricity and heat grids via district heating allows e.g.utilising surplus heat from power plants, industrial processes or waste incineration, while also allowing more renewable energy such as wind, large-scale solar thermal, and geothermal power in the system and thus enabling a fuel and cost efficient integration of fluctuating renewable energy sources [1,4].
Within smart grids there are many actors, such as Combined Heat and Power (CHP) supply systems for energy-intensive industrial applications or District Heating Networks (DHN).In order to optimally use the available capacities, optimal control methods are needed.These can be achieved by solving a unit commitment problem with Mixed-Integer Linear Programming (MILP), which is the main focus of this paper.State-of-the-art control systems or simulation tools can use forecasts and models of the operating system within a certain prediction horizon to react to imminent changes, e.g. the predicted heat demand of an industrial process or a DHN, or electricity prices (in the case of this paper on an hourly basis).In Leobner et al. [5], two different simulation-based concepts to integrate demand-response strategies into energy management systems in the customer domain of the Smart Grid are presented, highlighting demand response as an efficient way to harmonise demand and supply, but also pointing out extensive research demand concerning reduction of modelling effort, development of interfaces between systems supporting demand response and potentials.An overview of various approaches to solve the UC problem can be found in Padhy [6].In this paper, the modelling of a CHP supply system for Mixed-Integer Linear Programming Unit Commitment (MILP-UC)-problems is discussed and a method for classification and generic model generation of different CHP technologies with certain model capabilities and restrictions is presented.The method is validated by comparison of models and sample data of two Steam Turbines (ST) of different types to investigate feasibility and accuracy.

Combined Heat and Power supply systems
CHP generation is achieved by using the waste heat of an electricity generating unit, e.g. a Gas Turbine (GT), to satisfy a heat demand as in industrial processes in the energy-intensive industry like paper and pulp, steel, chemical or cement industry, or a DHN.Modern CHP units have high power-to-heat ratios and total efficiencies of about 90% or higher and contribute to a high primary energy utilisation, e.g. the block-type thermal power plant built in Kiel, Germany [7].Thermal energy storages are commonly used to enable a more efficient operation of the energy supply systems and ensure the satisfaction of the heat balance.DHN often have thermal energy storages available, mostly realised by hot water tanks, the hydraulic networks themselves have as well a certain storage capacity.Some CHP-technologies are e.g.Block-Type Thermal Power (BTTP) station, Combined Cycle (CC) power plant consisting of a GT with Heat Recovery Steam Generation (HRSG) and a ST with Supplementary Firing (SF), GT with HRSG only, biomass-fired boiler with ST, coal-fired boiler with ST or waste incineration with ST.Besides CHP, there are Power-to-Heat (PtH) technologies as heat pumps or electrode boilers transforming electric power to heat.Regarding the ST, there are various types like the Back Pressure Steam Turbine (BPST) or Extraction-Condensing Steam Turbine (ECST).From the MILP-modelling perspective, all these technologies can be separated into two classes: coupled and decoupled CHP generation.

Coupled and decoupled Combined Heat and Power generation
Coupled CHP generation means that heat generation Q, electricity generation P and fuel consumption F are linearly dependent, resulting in only one Degree of Freedom (DOF).Only one of these parameters needs to be considered as optimisation variable in the MILP-UC-problem, the other parameters can be derived from the chosen optimisation variable.Decoupled CHP generation has two DOF in either combination of the parameters Q, P and F, requiring two optimisation variables.The third parameter can again be derived from the two chosen optimisation variables.The method for modelling coupled and decoupled CHP units for the MILP-UC-problem will be discussed later.

Back pressure and extraction-condensing steam turbines
With regard to steam extraction from steam turbines for process heat, a general distinction is made between BPST and ECST.In BPST, the steam in the turbine expands to a fixed pressure and satisfies a downstream heat demand, resulting in lower electric efficiencies than pure condensing turbines.The heat and electricity generation of BPST is therefore coupled.In ECST, a part of the steam is extracted at certain turbine stages to satisfy heat demands, while the rest of the steam expands and condenses.Depending on the required steam temperature and mass flow for the downstream process, the extraction ratio can be varied.Remaining steam, which is not used for the heat decoupling, is supplied to the condenser after full expansion in the turbine, therefore the heat and electricity generation is decoupled [8].

Overview of Mixed-Integer Linear Programming
As the name indicates, only linear formulations with integer or continuous variables are allowed in MILP.This restricts the model capabilities, considering that most physical phenomena have non-linear characteristics.However, nonlinearities can be approximated by piece-wise linearisation, if necessary.Furthermore, MILP solvers are highly advanced and ensure global optimality in contrast to nonlinear solvers.In combination with increasing computing power over the last decades, MILP-UC problems can be solved in reasonable computation times.Two of the currently advanced solvers for MILP are Gurobi and CPLEX [9][10][11].
The mathematical formulation for MILP is stated in eq. ( 1): The objective is to minimise the cost function with the weighting vector c.The vector of optimisation variables x consists of the vector of integer variables xint and the vector of continuous variables xcon.Matrix A and vector b express all inequality constraints of the problem.The vectors lb and ub are lower and upper bounds of the variables.
There are several approaches for MILP algorithms such as branch and bound, cutting plane, decomposition or logic-based methods [12].Many solvers use the branch and bound method, which is based on relaxation, branching and bounding.First, a Linear Programming-Relaxation (LP-Relaxation) for the MILP-problem is solved, meaning, that the integrality constraints are ignored.In the branching step, the LP-Relaxation is divided into two branches, choosing a relaxed integer variable and setting the upper integer value of the LP-relaxed solution as lower bound for one branch and the lower integer value as upper bound for the other branch.If a solution is found satisfying all integrality constraints, any other branches with higher objective values can be eliminated (bounding) [13].Size and strength of the LP relaxation as well as the effect of branching have a strong impact on the computation time.The strongest formulation for the LP-relaxation is called sharp or is denoted as convex hull formulation [14].
With regard to the modelling of the unit commitment problem of energy supply systems, mixed-integer variables are well suited to model system characteristics, such as minimum part load or minimum up-or downtime of a unit or a nonlinear characteristic, separated in multiple linear segments.In linear programming, the feasible operating range has to be convex, meaning that for any two points within a convex set every point on the straight line segment between them belongs to the same set [12].With MILP, a non-convex operating range can be split into multiple convex ones.There are Special Ordered Sets (SOS) of type 1 and 2. In type 1 SOS, only one variable of the set is allowed to be nonzero, e.g. for modelling of discrete variables.Type 2 SOS allow two adjacent variables of the set to be nonzero, e.g. to model piecewise-linear functions [15].
MILP problems count among the complexity class Non-deterministic Polynomial (NP) time, which indicates the computational effort to solve the problem.The computational effort of many MILP-problems with n integer variables scales exponentially with approximately 2 n .The amount of integer variables is therefore essential for adequate computation times [15,16].
A general formulation for the MILP-UC is proposed in Carrion and Arroyo [17].There are several further studies, e.g. on tight formulations for minimal up-and downtimes [18], intra-period load gradient changes [19] or different approaches to model start up costs [20].In Morales-España et al. [21] and Fan and Guan [22] MILP-modelling methods for CC power plants are presented.To model a CC power plant there are mainly three methods: aggregate modelling, component modelling and configuration-base modelling.These works only consider the electric generation of CC power plants, not the CHP generation.In this paper, a component-based model of a CHP CC power plant is presented.In Mitra et al. [23], a generalized mode model on component basis for CHP power plants is presented, claiming to improve the profit up to 5-20%, depending on the case.A MILP formulation has been employed in Alipour et al. [24] to model the non-convex feasible operation region of CHP units, satisfying the technical constraints of generation units.

MATERIALS AND METHODS
From the MILP-modelling perspective, the regarded CHP technologies can be classified in coupled and decoupled CHP generation.The parameters for the energy supply system are heat and power generation Q and P as well as fuel consumption F (representing MWh/h, as in Figure 1 to Figure 5 and Figure 7, since these parameters represent transferred energy within a time step, not directly power) and the binary variable B for each unit, being 1 if the plant is running and 0 if shut down.As mentioned above, the MILP formulation is restrictive, as there are only two coefficients available for the conversion from the optimisation variable to another dependent parameter for a coupled CHP unit, representing a straight line in the three dimensional P-Q-F-diagram, or three coefficients for a decoupled CHP unit, representing a plane, as can be seen in the results section.
In eq. ( 2), the formulation of the objective function for a MILP-UC-problem, as addressed in this paper, with Nu units over NP time steps with fuel costs cF and fuel consumption F, fixed operating costs cB and binary commitment variable B, profits from heat and power generation with heat price cQ and heat generation Q, as well as electricity price cP and electricity generation P and start up and shut down costs SU and SD is given:

Coupled Heat and Power generation
For coupled modelling, it is assumed that heat and power are linearly dependent due to the restrictions of MILP.The fuel consumption F is modelled as well, being linearly dependent from both Q and P, resulting in one DOF.Therefore, besides the binary variable B, only one additional parameter is necessary as optimisation variable.In this paper, the parameter Q is chosen as optimisation variable.Parameters P and F can be derived according to eqs. ( 3) and ( 4): The respective terms of P and F in the cost function have to be added proportionally to the objective term of Q and B as shown in eq. ( 5): min ( " ) % * " + % " # % & " '/ min7 ( % 5 ( " * % 6 ( " / " ) % & % 5 & " * % 6 & " " '8 (5) min ̃( " ) % ̃& " '/ This formulation allows a minimum part load and varying efficiency from part load to full load.

Decoupled heat and power generation
In this case, the parameters Q and P are not coupled anymore, therefore both are needed as optimisation variables.The fuel consumption F can be expressed according to eq. ( 6): As for coupled CHP units, the objective function can be simplified according to eq. ( 7): min ( " ) % * " + % " # % & " '/ min 7 ( % 6 ( " / " ) % * % 6 * " " + % & % 6 & " " '8 The coefficients for these equations can be expressed by plant characteristics, such as minimum and maximum heat and power generation as well as fuel consumption.Of course, other parameters such as the thermal and electrical efficiency for full load and part load could be used instead, as well by adaptations of eq. ( 8) and (9).In this paper, the parameters of the minimum part load and the full load operation (Фmin and Фmax, Ф represents an arbitrary parameter) are chosen as defining parameters for the coupled CHP units.For decoupled CHP units, an additional (auxiliary) operating point (Фmin), not on the line segment between the other two points, is needed.Best suited are operating points on the vertices of the feasible operating range.The reason for the additional point for decoupled units is that the operating range is not a straight line anymore, but a plane that requires three points for definition.Table 1 shows potential defining parameters for coupled and decoupled units.
The coefficients for P and F in eq. ( 3) and ( 4) for coupled CHP units are derived from any two operating points (according to Table 1: max…1, min…2) as follows: The coefficients for F in eq. ( 6) for decoupled CHP units are derived from any three operating points (according to Table 1: max …1, min …2, aux …3) as follows: Depending on the necessary specifications of the plant being modelled, the corresponding constraints for the MILP formulation can be comfortably created.In the following, the formulations for ramp up and ramp down constraints as well as minimum uptime and downtime constraints are described.With these equations as well as values for minimum part load and full load for P, Q and F the feasible operating range can be expressed.
Constraints for ramp up with the maximum ramp up speed (RU) and ramp down with the maximum ramp down speed (RD) for any parameter Ф can be implemented with eqs.(10) and (11).The parameters P QR STU and P QI STU are the necessary minimal values of Ф at start up and shut down, the corresponding terms in eqs.(10) and ( 11) are necessary to cross the gap between 0 and the minimal part load Фmin, if RU or RD are smaller than Фmin: The minimum uptime (UT) and downtime (DT) for a unit can be modelled with eq. ( 12) and ( 13), according to Morales-España et al. [18]: For this purpose, as well as for applying start up and shut down costs, two additional binary optimisation variables Y and Z are introduced.For simulations with receding horizon, information about uptimes and downtimes of each unit at the first time step of the current prediction horizon has to be provided from the previously optimised time interval.
The objective of the UC-problem is to optimise the plant operation according to an objective function for a certain prediction horizon, while satisfying all constraints.
The prediction horizon has a length of NP time steps of size Ts, in our case as mentioned one hour, based on the hourly electricity price data from the spot market.Within this horizon the plant behaviour is modelled in discrete time steps, taking into account ramp up and down constraints, minimum up-and downtimes, varying fill level of storages, satisfaction of energy demands, etc.The size of the prediction horizon is limited, since the amount of integer variables increases proportionally with it and the computational effort grows exponentially.To simulate longer periods, e.g. a whole year on hourly basis, the optimisation can be done step-wise with a receding horizon, optimising only a small time interval (prediction horizon of size NP) of the whole simulation period at once with a certain overlap.There is the receding horizon step-size of NS time steps and an overlap of NO time steps with NP = NS + NO.The overlapping is necessary to take into account among others the inertia of the modelled system.The necessary overlapping size depends e.g. on ramp up and down constraints or storage capacities.To optimise the operation of a seasonal storage, the prediction horizon has to consider several months at least, while for a thermal energy storage for load smoothing a horizon of just several hours, up to a day is sufficient.
The model inputs are: • Plant configuration; • Heat demand; • Spot market electricity prices; • Fuel costs and others as Carbon dioxide (CO2) costs; • Subsidies for each unit, e.g. for CHP electricity generation in Germany [25]; • Size of the prediction horizon; • Step sizes; • Overlapping.The MILP method is capable of modelling minimum part load (semi-continuous) and full load constraints for Q, P and F, variable efficiencies that result of the linear dependencies of Q, P and F, minimum up-and downtimes as well as maximum ramp up and down constraints for Q, P and F. Costs can be applied either on the binary variable (fixed operating costs, if unit is running, start up and shut down costs) or Q, P and F (variable operating costs).
The energy balances of the system have to be satisfied, taking into account the heat demand of the system and the thermal energy storage with its maximum charging and discharging power as well as the maximum fill level according to eq. ( 14):

RESULTS AND DISCUSSION
In the following, a model of a BPST with coupled CHP generation and a model of an ECST with decoupled CHP generation are compared with real operation data of a BPST and an ECST in a DHN energy supply plant from four consecutive years (2012-2015).Figure 1 shows the total set of operation data (blue and red) for the BPST.Only operation data representative of the back pressure operation (blue) are taken into account for the modelling of the back pressure behaviour of the BPST.The deviations of the back pressure line can occur due to start up and shut down procedures, sampling error or suboptimal operation.Figure 2 illustrates the data of the ECST.Again, only representative data are used (blue).In case of the ECST the operation is not limited to the back pressure line.Due to the minimum heat extraction, there can be seen a gap between the fully condensing electricity generation and the extraction CHP mode.

Comparison of the coupled Back Pressure Steam Turbine model and data
For the modelling of a BPST the proposed method for coupled systems can be used with the three parameters: steam input S (replacing the fuel consumption F), heat output Q and power output P. Since the parameters are interdependent for coupled systems, any of those may be chosen as an optimisation variable.
The data of the BPST are compared with the proposed modelling method in Figure 3. Using the least squares method for the calculation of the two coefficients, a coefficient of determination (R 2 ) of 0.98537 for P and 0.98861 for S in dependency of Q is achieved.All corresponding coefficients and R 2 for the description of the models of BPST and ECST are listed in Table 2, with the expressed parameter Ф and the corresponding coefficients φP,Q,B.
Though there is an adequate fit of the linear model to the data, a slightly curved trend in the data can be seen, especially in Figure 3a.Such deviations or nonlinearities can be considered by separating the operating range into multiple sections by introducing additional binary variables for more accurate modelling of the real behaviour.For the purpose of the presented model for operation optimisation, this nonlinearity can be neglected, but it shows the limitation of MILP.
As a further illustration of the model accuracy, Figure 3c and 3d shows the varying efficiency η of the BPST kl ?V; Q m according to the model and the sample data.Though there is a spread in the data, the model lies midst within the data, suggesting adequate representation of the real behaviour.The relative deviations of η are higher for low values of heat and power generation, due to the approximately constant absolute deviations, seen in Figure 3a and 3b.

Comparison of the decoupled Extraction-Condensing Steam Turbine model and data
Figure 4 shows the plane of the ECST model operating range, describing the interdependency of the three parameters P, Q and S, as well as the corresponding data.There is an R 2 of 0.98537, therefore the model can be considered sufficiently accurate and the linear formulation holds not only for the BPST, but for the ECST as well.The model parameter S can be described as in eq. ( 6), again substituting F with S with the corresponding coefficients in Table 2. Parameters P and Q are optimisation variables.Figure 5 illustrates the dependency of the total efficiency η of the ECST of the parameters P and Q.The originating surface is not linear, since it is formed by a quotient of model parameters.For better illustration the surface is limited to 1, representing 100% efficiency.It is important to note that the surfaces seen in Figure 4 and Figure 5 are only valid for the model behaviour within the operating range of the modelled system.The mentioned gap, due to the minimum heat extraction, can be modelled by introducing an additional binary variable for setting the operation mode, either full condensing or mixed condensing/extraction.To this end, the additional binary variable would be used to define a minimal part load for the heat output Q. Depending on the desired accuracy of the model this gap can be ignored, as the overall behaviour and results will not change notably in many cases.Additionally, a reduction of binary variables can be crucial for simulations with high computational effort, since MILP problems belong to the complexity class NP, as mentioned above.

Combined cycle with Extraction-Condensing Steam Turbine
The method of the previous section can be used to model a CC consisting of three model components: GT (representing 2x Siemens SGT 800 according to manufacturer data) and SF with integrated HRSG as well as the presented ECST.In Figure 6, a simplified scheme of such a CC is shown, suited for MILP, not taking into account auxiliary components, e.g.feed water pump.The HRSG is not specifically modelled, but is presumed to be part of the GT and SF, since the output parameter of these components is steam.The GT is modelled as a coupled CHP unit with parameters PGT, SGT and FGT, the SF is defined by parameters SSF and FSF.The ECST is modelled as a decoupled CHP unit as above presented with the parameters PST, QST and SST.GT and SF both provide the steam SGT + SSF for the ECST.There is no further heat recovery from the exhaust gases after the HRSG in this example, meaning only the ECST provides heat (QST) for the downstream system, e.g. a DHN.
The behaviour of each component is described by constraints, as well as their relations.Table 3 lists the parameters, chosen optimisation variables and relations between the CC-parameters.As described previously, the coefficients can be calculated by defining operating points, or in the case of the ECST with the least squares method applied to the sample data.The mentioned model capabilities apply for the CC model as well.The CC has a more complex behaviour than the decoupled ST.This is accomplished by combining multiple units, gaining a higher DOF.As long as the described operating range of each unit is convex, the resulting combined operating range is convex as well and the modelling can be done without the use of additional integer variables.
Figure 7 shows the operating area of the modelled CC, having three distinctive operating areas.For better illustration of the three-dimensionality, the faces are projected on the grid.The actually feasible operating range is not only the three depicted faces, but rather a volume above these faces.Yet, the solutions of the MILP solver will mostly lie on these lower faces, since they represent the lowest possible fuel consumption for a certain combination of heat and power generation.This behaviour results from the chosen objective function of maximising profits of which a part is minimising fuel costs.The three faces in Figure 7 represent: • Part load to full load of GT without SF and variable extraction-condensing operation; • Full load of GT with SF and variable extraction-condensing operation; • Part load to full load of GT with SF and back pressure operation.The preferred operation of the presented models of BPST, ECST and CC lie on the vertices and edges of the feasible operating ranges, depending on the required heat generation, which complies with the linear programming theory.Points not on the vertices or edges may occur e.g.due to limited load gradients.With the presented generic MILP modelling method for CHP units, simulations of complex power systems and components like the CC can be conducted for prognoses of future demands.Another application is the investigation of different plant configurations for design optimisation.In succession, it can be used as a basis for optimal control and scheduling of the energy supply system.As was shown for a BPST and an ECST with comparison of model and real operation data, the MILP formulations are well suited for the modelling of these components, providing a compromise of accuracy and computation effort.

CONCLUSIONS
We discussed the need for more flexible and predictive optimal control concepts resulting in a UC problem, which can be solved with MILP.For future energy systems (especially 100% renewable) the interconnection between electricity, heat and gas networks through smart grids is essential.An example of such is an industrial process or a DHN with CHP generation, connecting heating and electricity networks.The modelling approach of coupled and decoupled CHP components for such systems, representing several technologies, e.g.BPST (coupled), ECST (decoupled), was presented.
By comparison of sample data and the MILP-model of a BPST and an ECST, a sufficient accuracy of the presented models was demonstrated.These can be created by use of the least squares method, if operating data is available, or by specifying the

Figure 3 .
Figure 3. Power output in dependency of: heat output (a); steam input in dependency of heat output (b); efficiency in dependency of heat output (c) and efficiency in dependency of power output (d)

Figure 4 .
Figure 4. Operating range for ECST steam consumption in dependency of heat and power generation

Figure 5 .
Figure 5. Operating range for ECST efficiency in dependency of heat and power generation

Figure 6 .
Figure 6.Scheme of a CC with GT, SF and an ECST

Figure 7 .
Figure 7. Operating range of the CC with fuel consumption in dependency of heat and power output

Table 1 .
Potential parameters for calculation of coefficients