Optimisation of the Operation of an Industrial Power Plant under Steam Demand Uncertainty

: The operation of on-site power plants in the chemical industry is typically determined by the steam demand of the production plants. This demand is uncertain due to deviations from the production plan and ﬂuctuations in the operation of the plants. The steam demand uncertainty can result in an inefﬁcient operation of the power plant due to a surplus or deﬁciency of steam that is needed to balance the steam network. In this contribution, it is proposed to use two-stage stochastic programming on a moving horizon to cope with the uncertainty. In each iteration of the moving horizon scheme, the model parameters are updated according to the new information acquired from the plants and the optimisation is re-executed. Hedging against steam demand uncertainty results in a reduction of the fuel consumption and a more economic generation of electric power, which can result in signiﬁcant savings in the operating cost of the power plant. Moreover, unplanned load reductions due to lack of steam can be avoided. The application of the new approach is demonstrated for the on-site power plant of INEOS in Köln, and signiﬁcant savings are reported in exemplary simulations.


Introduction
Combined Heat and Power (CHP) plants produce heat and power at the same time. This technology is considered to be sustainable and economic due to its potential in reducing the fuel demand and the green-house gas emissions compared to the conventional generation of power through steam condensation [1,2]. The efficiency of a CHP plant depends on its operation due to the considerable change of its performance in partial load regimes [3]. In most situations, a CHP plant has a significantly higher total efficiency factor compared to the traditional separate production of heat and power [4]. As an example, the efficiency of a gas turbine for power production is usually in the range of 36% to 40%, whereas the total efficiency factor reaches above 90%, when the rest heat is utilised as well [5]. In industrial sites, central utility plants convert steam to mechanical work to manage steam distribution at different pressure levels required for the operation of the plants, while co-generating electricity [6]. The role of CHP plants in the transition of the European energy systems has also been recognised by the European Union through the CHP directive (Directive 2004/8/EC, 2004), which promotes the construction and operation of co-generation plants [7]. In a more recent work, Jimenez-Navarro et al. [8] investigated the possible gains through the conversion of existing thermal power plants into CHP plants in Europe. Their result suggests that this transition improves the efficiency of the European energy system and results in a reduction of its operating costs and CO 2 emissions.
The steam and electricity that are required for the operation of industrial plants are often supplied by CHP plants which are either owned by the producing company or by an energy provider in the same industrial park. For a typical chemical production site, the power plant is of a high importance, as it incinerates off-gases produced by the production plants in addition to providing them with steam and electricity. The planning of the operation of CHP plants of such production complexes depends on the steam and the electricity demand of the production units and is typically heat-driven.
The optimisation of the operation of CHP plants has been an on-going research topic for many years. Tina and Passarello [9] developed a model for the CHP plant of a petrochemical complex with plants of several companies and used it to maximise the profit of the CHP plant, considering time-variant electricity prices. In this work, legislative regulations as defined by the European Directive 2004/8/EC are also taken into consideration. Marshman et al. [10] optimised the operation of a CHP plant that has to cover the steam demand of a pulp and paper mill production facility and can sell electricity to the power grid. Mitra et al. [11] presented a model for an industrial CHP plant which models the plant on a component level and accounts for the restrictions of the transitions between different modes of operation. They formulated the optimisation as a Mixed-Integer Linear Program (MILP) and implemented it to optimise the operation considering time sensitive electricity prices. Pablos et al. [12] presented a non-linear gray-box model of a CHP plant that is coupled to a sugar factory where the production rates of the whole system are optimised from an economic perspective, taking dynamic electricity prices into consideration. Bindlish [13] used a non-linear model for the scheduling of a CHP plant in a chemical site and used linear Model Predictive Control (MPC) to implement the optimised operation. In their work, the model parameters of the scheduling model are continually updated using plant data to increase the accuracy of the results. Furthermore, there are several works on optimising the operation of CHP plants by taking the balancing of the electricity grid and day-ahead spot markets into consideration [14,15].
For an industrial CHP plant, the presence of uncertainties in the heat and power demands makes the planning problem challenging, as an imbalance in the utility networks either results in wasted resources or incurs high losses for the other production plants due to shortcomings of steam. Salgado and Pedrero [16] presented the results of a survey on the short-term scheduling of CHP plants. Their work concludes that stochastic models should be developed and applied to CHP scheduling and that the uncertainty in certain parameters, e.g., the efficiency curves, needs to be addressed as well. The work of Palsson and Ravn [17] is an early study on short-term planning of the operation of CHP plants under uncertainty. Their work discusses the operation of a CHP plant that is coupled to a heat storage tank. Under the assumption of an uncertain power production, they solve a stochastic optimisation problem using the Progressive Hedging Algorithm (PHA) [18]. In a recent contribution, Leo and Engell [19] presented an approach for the day ahead energy procurement under uncertain equipment breakdown and electricity prices for a facility with an on-site power production. Kia et al. [20] applied linearisation techniques to the model equations of the units of a CHP plant and formulated a Stochastic Mixed Integer Linear Program (SMILP) to optimise the operation of CHP systems that have thermal and electrical storage capacities. Bornapour et al. [21] proposed a stochastic framework for CHP scheduling in micro-grids, taking into account uncertainties like electricity price, solar radiation and the speed of the wind.
Steam at various pressure levels is essential for the operation of many chemical plants. Their actual steam consumption is influenced by several factors, e.g., the plant utilisation or the ambient temperature. There are several works on developing models that represent the efficient energy performance of chemical processes for given operating conditions with different level of details. These approaches vary from linear regression models proposed by ISO 50006:2017 [22] to non-linear data-based models based on historical data [23] and detailed models based on physico-chemical relations [24]. The actual performance of the plants however deviates from the predictions of these models due to external influences and sub-optimal operation of the plants. The prediction of the mid-term future performance is even more challenging, considering the additional uncertainties of the factors that influence the performance of a process (e.g., variations of load and ambient temperature). The planning of the operation of a power plant in a chemical park, which provides steam to the process plants, requires the knowledge of the demand profiles of the plants. However, some decisions have to be taken well before the exact steam demands are known and have a long-term effect, e.g., the shutdown of a unit, which limits the flexibility of the power plant to react to demand variations. A deviation of the realised demand profile from the planned profile causes an imbalance of the steam network which can result in an inefficient operation of the power plant and can even impair the production of the other plants. As the operation of industrial CHP plants is predominantly heat-driven, it is of high importance to take the uncertainty of the steam demand into account in the planning. This increases the resource efficiency of the production of steam and helps to keep the steam network stable and meet the demand of the plants at all times. To the best knowledge of the authors there is no previous work done that addresses the challenges of scheduling of the operation of industrial CHP plants under steam demand uncertainties. This work presents a new approach to consider this aspect in the planning and closes the current gap in the literature.
In the next section, an overview of the CHP plant that is considered in this investigation is given and the challenges faced due to the steam demand uncertainties are discussed in detail.

Overview of the Power Plant at INEOS in Köln
In this contribution, the power plant of INEOS in Köln is considered as a typical example. INEOS in Köln operates a large integrated petrochemical production site with more than 20 production plants, including naphtha and gas crackers and a central CHP plant. Some of the by-products of the production plants cannot be processed further or be sold. These streams are sent to the CHP plant to be incinerated. Many of the production plants of the site require steam for their operations. The CHP plant utilises the heat that is produced through waste incineration for the generation of steam. However, the energy content of the waste streams is insufficient for satisfying the steam demand of the site, and therefore Natural Gas (NG) is imported as an additional fuel. There are four steam networks at INEOS in Köln with different pressure headers (60, 30, 15 and 5 bar). In order to match the requirements of the steam networks, the high pressure steam from the boilers is sent through several steam turbines which convert part of the thermal energy into electricity. In fact, the electric power that is generated by the CHP plant of INEOS in Köln is a by-product of producing steam on these different pressure levels. In addition, the power plant has a condensation turbine which produces electric power by consuming 5 bar steam. The cost-effectiveness of producing electricity in this turbine depends mainly on the prices of electricity and natural gas. The electricity production of the CHP plant covers approximately 30% of the electric power demand of the production plants and additional electricity is purchased from the public grid.
The steam network of INEOS in Köln is connected to an energy provider at the chemical park, from which the power plant can purchase steam at specific pressure levels within some predefined limits. Figure 1 depicts the schematic of the CHP plant of INEOS in Köln. It consists of five boilers and several multi-and single-stage turbines. It is possible to put any of these units in and out of operation and to vary their loads depending on the steam and electricity demand profiles and the dynamic electricity prices. There are upper and lower bounds on the load of each unit. A switch in the operation mode of some of the units is restricted by a minimum stay time in the new mode. e.g., if Boiler 3 (B3) is turned on, it will have to remain on for at least 24 h. A list of the mode transition constraints can be found in Table 1. The goal of the work is to calculate the optimal operating plan for the operation of the power plant for a fixed horizon in the future.    Unit  B1  B2  B3  B4  B5  T1  T2  T3  T4  T5  T6   On  24  24  24  24  24  2  2  2  2  2  -Off  4  4  4  4  4  2  2  2 2 --

Challenges
The planning of the operation of the power plant requires the knowledge of the steam demand profiles of each of the four pressure levels. Figure 2 presents the steam demand of one pressure header of one of the production plants of INEOS in Köln against the utilisation of the plant. As it can be seen, the steam demand varies even for the same value of the plant utilisation. A deterministic scheduling of the operation of the CHP plant can be done based on average steam demand profiles, which are called the "nominal demand" profiles hereafter. However, the actual steam demand is uncertain and can deviate from the nominal demand profile due to various reasons, e.g.,: • Influence of the operators: In many plants, the operators have a significant influence on the operation, as the more strategic decisions are not automated. • External influences: The model that is used for the calculation of the steam demand does not consider all the possible influences, e.g., the influence of the ambient temperature is often not described. • Deviations from the production plan: The actual production plan of the plants of the site deviates from the values that are used to calculate the steam demands, which results in the deviation of the actual steam demands from the nominal ones. • Disturbances of the operation: There are frequent disturbances of the operation, due to e.g., failure of pieces of equipment or lack of raw materials or other resources.
The power plant has a limited flexibility to react to these short-term disturbances. A deviation of the realised steam demand profile from the planned profile results in an imbalance of the steam network. If the amount of the available steam is higher than the site demand and depending on the electricity price of the power grid, the steam can be sent to the condensation turbine. If the market price of the electric power is negative due to an imbalance of the power grid, or the turbine has reached 100% utilization, the excessive amount of the steam must be vented, which is a waste of resources.
A shortcoming of steam on the other hand requires purchasing steam from the energy provider of the chemical park. In extreme cases, where the shortcoming is above the capacity of the energy provider, the production rates of the plants may have to be reduced. This results in a loss of production and potentially incurs high losses. It also results in imbalances in other networks at the site.
The long-term effect of some of the decisions about the operation of the CHP plant creates an additional complexity of this task. As mentioned earlier, a change in the operation mode of some of the units is subject to minimum stay times in the new mode. Therefore, a deviation from the nominal steam demand in the future can result in inefficient operation of the power plant. For example, after starting up a boiler it is bound to stay on for 24 h. If it is decided to turn a boiler on based on the nominal demand and then the realised demand is critically lower, the imbalance in the network has to be compensated by venting the steam, which is inefficient.

Problem Statement
The problem that is investigated in this contribution can be formulated as follows: the goal is to plan the operation of a CHP plant in an industrial park for a given time horizon in the future (the next few days). It is assumed that the electricity demands of the production plants are known and fixed for this period. Furthermore, it is assumed that the electricity prices are known for the planning horizon. The steam demand is assumed to be known for the next hours, but its future trajectory is not known exactly for the rest of the planning horizon. The following variables are degrees of freedom for the optimisation: • The operation modes of all of the boilers and the turbines; • The distribution of the fuels to each of the boilers; • The flow rates of all of the water/steam streams.
This contribution presents a stochastic programming approach to the optimal scheduling of the operation of a CHP plant under steam demand uncertainty. The proposed framework combines the advantages of two different scheduling approaches, preventive (robust) and reactive scheduling (see Section 2), in order to overcome the effects of the uncertainties.
The paper is structured as follows: Section 2 presents the theoretical background of the methods that are used in this publication. This is followed by the model of the CHP plant and the formulation of the optimisation problem in Section 3. Section 4 discusses the approach to modelling the uncertainties of the steam demand. In Section 5 the application of the approach to the CHP plant of INEOS in Köln is presented and the results for two characteristic scenarios are discussed.

Theoretical Background
The scheduling of a task can be divided into two main steps: the generation of an initial schedule and the revision of the generated schedule [25]. The generation of a schedule is a proactive procedure to determine the short-term production plan for the different units of a production facility in order to reach the production goals. This plan is created prior to the start of the production process, based on models that consider the constraints of the system. The revision of a schedule is a reactive element where the execution of the schedule and the unfolding of events is monitored, and corrective actions are taken to deal with the deviations from the planned execution.
Preventive scheduling approaches consider the uncertainties during the generation of the schedule. These approaches take the possible realisations of the uncertainty into consideration, and calculate a good solution for all of the cases. Reactive scheduling methods are applied to modify a nominal schedule, as a reaction to the realisation of the uncertainties, to unexpected events or updated parameters of the system during the production process [26].
Stochastic optimisation is a mathematical optimisation technique that includes uncertain parameters. Li and Ierapetritou [27] claimed in their review paper that stochastic optimisation is the most commonly used approach in the literature for preventive scheduling and divide stochastic programming models into two main categories: two-stage or multi-stage stochastic programming approaches and chance constraint programming based approaches. Two-stage stochastic optimisation was introduced by Dantzig [28] for linear models in the 50s, and later Birge and Louveaux [29] extended it to multi-stage. In multi-stage stochastic optimisation, a stage denotes a period of time, for which some values of the uncertain parameters are unveiled at its beginning, while for the future periods, these parameters are uncertain. At each stage, a set of decisions has to be made based on the available information and these decisions are implemented. The future decisions can be adapted to the observed evolution of the uncertainty.
A moving horizon scheduling strategy is a reactive scheduling approach which solves an optimisation problem iteratively over a limited look-ahead horizon into the future, implements a set of decisions and then moves the optimisation horizon forward, updating the parameters at the beginning of each iteration.
In this work, we consider the combination of two-stage stochastic optimisation and moving horizon scheduling, similar to [30] for production scheduling.

Two-Stage Stochastic Optimisation
In two-stage stochastic optimisation, the decision variables (degrees of freedom) are divided into two sets: the first-stage (here and now) decisions and the second-stage (wait and see, or recourse) decisions. The first-stage decisions are the ones that have to be made based on the current information. The second-stage decisions are taken after the realisation of the uncertainty and as a reaction to their effects. As the result of the stochastic optimisation, the first stage decisions are taken such that after the adaptation of the recourse decisions the solution is optimal with respect to the expected value of the cost function and the constraints are met for all realizations of the uncertainty. For the linear case, the standard formulation of two-stage stochastic problems is given by [29,31]: Here, X and Y are the polyhedral domains of the first and second-stage variables, x and y(ω). Equation (2) represents the constraints that incorporate the variables of the first stage. In the second stage, a number of discrete realisations of the uncertain parameters ω are considered. The vectors h(ω) and q(ω), and the matrices W(ω) and T(ω), which correspond to the second stage, depend on the uncertain parameter ω. Equation (4) relates the variables of both stages and has to hold for all realisations of the random parameters. The objective function, Equation (1), is composed of a deterministic term c T x and the expected value of the objective of the second stage, q(ω) T y(ω), over all the possible realisations of ω.
A two-stage Stochastic Mixed Integer Program is a problem in which some of decision variables have to satisfy integer restrictions. If a linear model is assumed and the uncertainties are represented by a set of discrete scenarios S , the SMIP can be reformulated as a deterministic equivalent in the form of an MILP as follows [32]: where s is the index of the scenarios and p denotes the probabilities of the scenarios. A class of SMIP problems, which are called problems with continuous recourse, consider integer restrictions only for X. If Y includes integer restrictions, the SMIP is called a problem with integer recourse. Such problems are challenging to solve if the problem size becomes large, as the well-known decomposition methods, e.g., Benders decomposition, are not directly applicable [33].

Representation of the Scenario Tree
The scenario tree of the stochastic problem introduced in Section 2.1 is presented in Figure 3a. The first-stage decisions are made prior to the realisation of the uncertainties and are therefore shared among all of the scenarios. A complicating feature of this representation is that the first-stage variables appear in all scenarios. Using this representation, direct decomposition approaches like Langrangean decomposition or progressive hedging are not applicable [34]. Ruszczyński [35] proposed an alternative representation of the scenario tree as shown in Figure 3b. In this approach, a separate set of variables is defined for the first stage of each scenario. However, the variables that correspond to the first stage of different scenarios cannot take different values. This is shown by the vertical lines at the first stage of the scenario tree. These lines represent the so called non-anticipativity constraints which must be incorporated into the optimisation problem as additional constraints to ensure the equality of the first-stage variables for the different scenarios. Figure 3. A scenario tree of a two-stage stochastic problem and its alternative representation according to Ruszczyński [35]. (a) Original representation of scenario tree. (b) Alternative representation of scenario tree.
In this work, the alternative form of the tree as shown in Figure 3b is used and the formulation of the non-anticipativity constraints is discussed in Section 3.8.

Optimisation on a Moving Horizon
In real-world problems, the uncertain events unfold over time and the decisions have to be adapted to react to these realisations. In reality, all real-world planning problems over a finite horizon can be considered as multi-stage problems, where at each stage decisions have to be made for the future based on the information that is available at that stage. Ideally, an infinite horizon should be taken into account, but this is not feasible. Therefore, in the moving horizon approach the optimisation problem is solved over a limited horizon which is iteratively shifted into the future and updated. Solving the optimisation problems consecutively on a moving horizon has the advantage of being able to regularly update the stochastic optimisation with the new information that becomes available, while keeping the problem size manageable [30].
The two-stage stochastic optimisation formulation presented in Section 2.1 is an approximation of multi-stage problems by assuming that all uncertainties for the short planning horizon are realised after the end of the first stage. This is a simplification of the real decision-making problem. In combination with the moving horizon formulation, the gain of information is considered by the update of the problem for the new horizon. In each iteration, as shown in Figure 4, the problem parameters are updated based on the new information from the plants to reduce the model uncertainty. Furthermore, in our approach, the scenario tree of the two-stage stochastic problem is updated at the beginning of each rolling horizon based on the current information, which is explained in more detail in Section 4.1. Afterwards, the optimisation is executed and the optimal values of the first-stage variables are implemented. The optimisation horizon is shifted by the length of the first stage and the procedure is repeated.

Computation of the Benefits of the Stochastic Optimisation
A simulation framework is used in this contribution to calculate a quantitative indicator of the advantage that the presented stochastic optimisation framework provides. This approach, presented graphically in Figure 5, simulates two operators which optimise the operation of the power plant in parallel on a moving horizon. One operator uses two-stage stochastic optimisation and the other one performs a deterministic optimisation based on the expected value of the different steam demand scenarios.
At the beginning of each optimisation instance, the parameters of the model are updated based on the current information about the state of the power plant, and the scenario tree of the optimisation is adapted for both strategies based on the current realisation of the uncertain parameters. Afterwards, the optimisations are executed, the decisions which correspond to the first stage of the two-stage stochastic optimisation are implemented for both operators and the horizon is shifted to the beginning of the second stage and the procedure is repeated until the final optimisation instance for the planning horizon is reached. A reasonably large number of realisations of the uncertainties are simulated to represent the assumed distribution of the steam demand scenarios.

Model of the Power Plant
The model of the power plant is formulated as a MILP and includes the mass and energy balances of several compartments like pre-heaters, boilers and turbines. In order to avoid non-linearities, the pressure and the temperature of the streams and consequently their enthalpies are not decision variables and are assumed to be constant. These parameters are however updated at the beginning of each optimisation instance based on new measured data from the plant (see Section 2.2). Figure 6 presents a visual illustration of the balancing domains of the compartments of the model. The overall mass balance equates the overall amount of the water/steam that enters and leaves the plant and is discussed in Section 3.1. The water streams that enter the power plant go through different units/section and merge with/split to several water/steam streams at certain sections. These units/sections are defined as internal balance nodes and are discussed in more depth in Section 3.2.

Overall Mass Balance
The overall water inlet to the power plant consists of: •ṁ ts W,in,B,total : The total fresh water sent to the different boilers. In case a boiler has a pre-heating system, this stream enters the first internal balance node of this system. Otherwise, the stream enters the internal balance node of the first pump prior to the boiler. •ṁ ts W,inj,total : The total fresh water sent to the internal balance nodes with injection streams (see Section 3.2). •ṁ ts S,in,P,total : The steam which is sent from the steam networks to some of the preheaters of the boilers. The values of these variables are determined by the energy balance of the respective internal balance node (see Section 3.2). •ṁ ts S,in,K ,total : The steam which is sent from the steam networks to the condensation turbine.
The overall outlet of the power plant consists of steam flows at pressure levels of 60, 30, 15 and 5 bar which are sent to the network of the different steam headers and the condensate from the condensation turbine. The set of steam headers is defined as D = {ST60, ST30, ST15, ST5, STC}, where STC stands for the condensate. The total amount of steam export from the power plant is then calculated as: where T and S are the sets of time steps and scenarios, respectively. C denotes the subset of the internal balance nodes that are connected to the steam network andṁ ts S,c,d,exp denotes the flow rate of the output of node c that is sent to the steam network d. The overall mass balance of the power plant is formulated as: m ts W,in,B,total +ṁ ts W,inj,total +ṁ ts S,in,P,total +ṁ ts S,in,K ,total =ṁ ts S,exp,total ∀t ∈ T , s ∈ S , (10)

Internal Balance Nodes
The water/steam streams pass through or are merged/split in the balance nodes. The set of all nodes (N ) is composed of pumps, pre-heaters, boilers, turbines, valves and stream junctions. For all nodes a mass balance is formulated as: where the sets I n and O n represent the input and output streams of node n. In a subset of the nodes which consists of the pre-heaters and pumps P, boilers B, turbines K and let-down stations of the turbines L , a significant change in the water/steam properties is happening. Consequently, an energy balance is formulated for these nodes in addition to the mass balances as: where h t i denotes the enthalpy of the corresponding streams at time t. The variables P ts in,n and P ts out,n are the additional energy imported to or exported from a balance node n at time t under the realisation of scenario s.

Injections
In a subset of the balance nodes, H ⊂ N , fresh water is injected into an overheated steam stream to adjust its temperature. These streams are denoted asṁ ts W,inj,h . Figure 7 shows the ratio of the injected water to steam for one of the injection points on an hourly basis for two years. As it can be seen, other than at a few points in the range of a low mass flow, which is not a normal operating condition, the ratio of these two streams remains in the same range. The average relative deviation of the measurements from the mean of the data points with a ratio below 0.2 is 6.6%. In the following, this ratio factor is therefore assumed to be constant for each optimisation instance and is updated at the beginning of each iteration based on the plant measurements.

Model of the Boilers
Each boiler b has a set of burners G b which incinerate the fuels and produce the heat for the steam production. The energy released in the burners is calculated based on the Lower Heating Value (LHV) of the fuel which is sent to them. The energy that is transferred from the burners of the boiler to the water is calculated as: P ts whereṁ ts BR, f bg is the mass flow of the fuel f that is sent to the burner g of the boiler b at the time step t for the scenario s, and P ts BR, f bg is the produced heat in the respective burner. LHV t f b is the lower heating value of fuel f for boiler b at time t and η t f b is an overall efficiency factor of fuel f in boiler b. Previous works on modelling of boilers showed, based on empirical data, that for a known design and an operating range of above 50% of the maximum load, the variations of the efficiency factor are very small (e.g., Aguilar et al. [36] and Varbanov et al. [37]). At INEOS in Köln the technical lower limit of the operation of the boilers is above 50% of the maximal load and therefore the efficiency factor is assumed to be constant. The variable P ts BO,b is the total heat transferred from boiler b to the water at time t for scenario s and enters the energy balance of the boiler as the imported energy to the balance node (see Equation (12)).
In each boiler, some burners can incinerate only a specific fuel, whereas others have the possibility to combust different ones. There are limits on the minimum and maximum fuel mass flow that each burner can combust, which depends on the fuel type and the boiler. These restrictions are imposed by Equations (15) and (16), where binary variables y ts BR, f bg are defined for the burners and take the value of 1 if the fuel f is sent to the burner g of the boiler b at the time step t for scenario s. Furthermore, each burner can combust only one fuel at a time which is insured by Equation (17).
where LB BR, f b and UB BR, f b are the upper and lower bounds of the mass flow of fuel f for the burners of boiler b. An additional restriction exists due to the fact that some fuels cannot be fed to specific boilers. Even though this constraint can be imposed on the mass flow of the fuels sent to the burners by setting both LB BR, f b and UB BR, f b to zero, it then does not limit the value that the binary variables of the burners can take. To avoid this, binary variables y ts FL, f b are defined for the fuels, which take the value of 1 if the fuel f is sent to the boiler b and zero otherwise and the following constraints are formulated: Equations (18)- (20) make sure that if a fuel is not allowed to be sent to a boiler, all of the corresponding binaries of the burners are equal to zero. Moreover, if a fuel is assigned to a boiler, at least one of the burners of the boiler must be assigned to it.
Each boiler has a limited range of production of steam when it is in operation. The amount of steam produced in each boiler b ∈ B is equal to the output of the corresponding balancing node, which is calculated as ∑ j∈O bṁ ts W,j (cf. Equation (11)). Considering that each boiler has only one output stream, which is the amount of the steam produced, these streams are represented asṁ ts S,b,out for the sake of better clarity of the representation.
The following constraints are then defined to restrict the operation of each boiler to the allowed range:ṁ where LB b and UB b are the lower and upper bounds of steam production of boiler b. y ts b is a binary variable that is equal to 0 if the boiler is out of operation and 1 otherwise.

Model of the Fuel Inventory
The amount of the overall fuel that is available for the optimisation horizon is defined as the predicted value of the waste products and an unlimited amount of natural gas. The inventory of the fuels is thus modelled as: where Inv ts f is the available amount of fuel f at time t and scenario s. Inv 0 f denotes the amount of the available fuel at the beginning of the optimisation horizon.

Model of the Steam Turbines
The production of the electric power in the turbines is modelled based on the so-called Willan's lines [38]. These are linear equations which relate the electric power output to the mass flow of steam. Analogous to the outlet of the boilers, the outlet of each turbine, ∑ i∈I kṁ ts W,i is denoted byṁ ts S,k,out for the sake of the simplicity of representation and the electric power that is produced in the turbines is modelled as: P ts el,k = a 1,k ·ṁ ts S,k,out + a 2,k ∀k ∈ K , t ∈ T , s ∈ S , where a 1,k and a 2,k are the slope and the intercept of the Willian's line for turbine k.
For some of the turbines, a bypass stream exists through which the steam is sent to a let-down station in case of emergency and when the turbine is out of operation. At such letdown stations, the steam is expanded to the pressure of the turbine outlet. During normal operation, a small amount of steam is also bypassed due to the technical considerations. The flow rate of this stream is controlled to a constant ratio of the outlet stream of the turbine. In order to consider this in the model, a constraint is added which makes sure that the bypass stream ratio is higher than the minimum ratio: where U is the subset of the turbines which have a bypass stream.ṁ ts BP,u is the flow-rate of the bypass stream and the µ u is the minimum technical flow ratio.
Similar to the boilers, each turbine has minimum and maximum steam flow rates. This is incorporated into the model using the following constraints: where LB k and UB k are the lower and upper bounds of the steam inlet flow for turbine k. y ts k is a binary variable that is equal to 0 if turbine k is out of operation and 1 otherwise.

Mode Transitions
The boilers and the turbines, as shown in Sections 3.3 and 3.5, can be either in or out of operation. Therefore, a set of operation modes is defined as M = {on, off }. The binary variable y ts r,m takes the value 1 if unit r is in mode m at time step t for scenario s. At each point of time only one mode must be active for a unit. y ts r,on = y ts r ∀t ∈ T , s ∈ S , r ∈ R, R = K ∪ B ∑ m∈M y ts r,m = 1 ∀t ∈ T , s ∈ S , r ∈ R.
The mode change of units in R is bound to a minimum stay time in the new state (see Table 1). For this purpose the binary variable z ts r,m,m is defined which takes the value 1 if and only if the mode of unit r changes from m to m from time t − 1 to t. The following constraints are formulated which make sure that if a unit has switched mode from m to m , it remains in the new mode for a minimum duration of θ r,m,m after the change [11]:

Demand of the Production Plants
The production plants on the site require steam and electric power for operation. As mentioned in Section 1.1, the power plant must satisfy the steam demand of the site and covers a part of its electricity demand.

Steam Demand
The production plants require steam at various pressure levels for their production. As mentioned in Section 1.2, the steam demands are to some extent uncertain. The parameter D ts S,d denotes the total steam demand at pressure d at time t for scenario s. The power plant can satisfy the steam demand by producing the full amount itself. Additionally, it can purchase steam at 30, 15 and 5 bar from an external energy provider in the chemical park. However, there is a limit on the amount of steam that the energy provider of the chemical park can provide to INEOS in Köln. If the deficiency in the production of steam is higher than this limit, the loads of the production plants have to be reduced, which may result in significant economic losses and imbalances in the other networks of the site which may have a cascading effect.
The balance of the 30 and 15 bar steam headers is formulated as: whereṁ ts S,d,imp denotes the mass flow of the steam at pressure level d that is imported from the energy provider at time t for scenario s.
The steam at the 60 bar network can be sent to the pre-heaters of boilers (B1-B3). Consequently, the balance of this header is formulated as: ∀t ∈ T , s ∈ S .
The condensation turbine is connected to the 5 bar steam header. In case of excess in the header that cannot be sent to any of the plants or to the condensation turbine, steam has to be vented. The balance of this header is formulated as: m ts S,ST5,exp +ṁ ts S,ST5,imp −ṁ ts S,in,K ,total −ṁ ts S,vented = D ts whereṁ ts S,vented is the amount of the vented steam andṁ ts S,ST5,imp is the amount of the imported steam from the energy provider at time t for scenario s.

Electricity Demand
The production plants require electricity for their operation. Even though electricity is produced by the steam turbines of the power plant, the capacity of the electricity production is not high enough to meet the overall electricity demand and reaches only about 30%. The production plants have to purchase the rest of the required electricity from the power grid. The balance of the electricity network of INEOS in Köln is formulated as: where D t el denotes the total electricity demand of the production plants at time t. As the focus here is on how to handle the uncertainty in the steam demand of the plants, it is assumed that the electricity demand of the production plants is certain for the optimisation horizon and is independent from the scenario s. P ts el,imp indicates the amount of the electricity imported from the power grid at time t for scenario s.

Non-Anticipativity Constraints
As mentioned in Section 2.1.1, this work uses a representation of the scenario tree that introduces separate first-stage variables for each scenario. Therefore, non-anticipativity constraints are required to guarantee that the same decisions are made for the time periods where the scenarios are identical. These constraints are formulated as: where T 1 ⊂ T presents the time steps that belong to the first stage of the scenario tree. x ts denotes all of the scenario dependent variables that are presented in Section 3.

Objective Function
The objective function of the optimisation problem is defined as the profit of the power plant. The following factors affect the overall profit: • Cost of the fuel that is purchased in addition to the off-gas streams from the production plants; • Cost of the steam that is purchased from the energy provider; • Cost of the water for the vented steam (the cost for the sold steam is implicitly considered in the steam price); • Income from the steam that is sold to the production plants; • Income from the electricity that is produced in the power plant and sold to the production plants.
Furthermore, an additional penalty term with a large value is defined for steam shortcomings above the amount that can be imported from the energy provider, in order to represent the adverse effects of such imbalances in the steam networks which forces the production plants to reduce their loads. As the cost of steam and electricity is only a small part of the production cost and the margins of the plants are high, this penalty is comparatively high. The objective function of the optimisation is formulated as follows: where p s is the probability of the scenario s and: The parameter c S,d is the price of the steam at pressure level d that the power plant sells to the production plants and is assumed to be constant for the planning horizon. The parameter c S,d,imp is the price for importing steam at pressure level d from the energy provider, which is higher than its selling price (c S,d ) to the production plants. c el,t is the hourly electricity price from the power grid, which is assumed to be known for the optimisation horizon. The variableṁ S,d,shr is the amount of shortcomings in the steam network which cannot be covered by the energy provider and c pen,d is the penalty price for the associated losses and operational problems. Finally, c W is the price of fresh water.

Modelling the Uncertainties
The uncertainty of the steam demand is modelled as follows. As mentioned in Section 1.2, when there are no deviations from the production plan, the actions of the operators and the external influences as, e.g., the ambient temperature are two of the main reasons for the deviations of the steam demand trajectory from the nominal demand for a given production plan. In this work, it is assumed that the plants of the site follow the production plan without any deviations. Considering the shift system of the company, in which the operators change every 8 h, and the fact that the external factors like ambient temperature do not deviate from their predicted values for the next 8 h drastically, it is assumed that the steam demand trajectories for the first 8 h are known and certain. Therefore, the first stage is defined as the first 8 h. The goal is to make the operational decisions of the power plant for the next 8 h by taking the partly uncertain demands of the following 2 days into consideration. The demand after the first 8 h is uncertain and deviations from the nominal demand trajectories that are calculated based on the currently available information are possible. Therefore, the second stage is defined as the 48 h after the end of the first stage with two different possible deviations from the nominal demand trajectory for each of the steam networks; a higher demand than the nominal trajectory and a lower demand. These two other potential trajectories are obtained by multiplication of the predicted demand trajectory with a constant factor. Under normal operating conditions of the site, the production plants require 5, 15 and 30 bar steam for their operation. In this work it is assumed that no correlation exists between the realisation of the uncertainties of the different steam networks and therefore all of the possible combinations are considered. Figure 8 provides a qualitative representation of the scenario tree. Each scenario has a specific combination of the realised demand regions of the steam networks, represented by a multiplication factor for the nominal demand α d,i , where i denotes the realised demand region as high (h), nominal or medium (m) and low (l). Taking the uncertainty of all three networks into consideration results in 3 3 = 27 scenarios for the second-stage of the optimisation problem.
, , , , , Figure 8. The scenario tree of the uncertain steam demands. In each scenario the dynamic nominal demand trajectory of steam network d, which is calculated based on the currently available information, is multiplied with the corresponding multiplication factor α d,i .

Adaptation of the Uncertainties
The definition of the scenario tree and the probabilities of different scenario realisations are important elements of calculating representative results of two-stage stochastic optimisation. Therefore, the historical data was investigated to check if any specific dependencies between the current situation of the site and the probability of the scenarios in the future exists.
In this analysis, the plant data presented in Figure 2 is investigated which presents the hourly steam demand of a production plant of INEOS in Köln. As proposed above, two additional demand profiles are defined as deviations around the nominal demand profile, presented in Figure 9. The steam demand is consequently divided into three regions, the boundaries of which are calculated as the average of two consecutive demand profiles and are indicated in Figure 9 as the thin lines between them. Each of the data points is assigned to one of the regions: the points lying below the lower boundary are assigned to the the low demand region, the data points lying between the two region boundaries are assigned to the nominal demand region and the rest of the points to the high demand region. Later, the demand region of each data point is compared to the demand region of the data point that represents the steam demand of the plant after 8 h, the starting point of the second stage of the current and the first stage of the next optimisation interval. Doing so, the frequency of the changes from the current demand region into another region after 8 h is analysed. The results are listed in Table 2, where α d,j,i denote the realised demand region i at interval j for pressure level d. It can be seen that the frequency of remaining in a similar region after 8 h is the highest. In other words, according to demand data, a drastic change of demand in two consecutive time intervals is not very likely. Therefore, at the beginning of each iteration, the probability of the scenario tree is updated based on the information on the current operating region.

Magnitude of the Uncertainties
The overall steam demand of the site from each network is calculated as the aggregated amount of the steam demand of the production plants: where D t S,d,v denotes the uncertain steam demand of plant v from the steam network d at time t and: where σ D S,d,v denotes the uncertainty of the steam demand of plant v from pressure header d. Assuming that the steam demand of the different plants is not correlated, the uncertain demand from each steam network is calculated as: where and In an ideal case the magnitude of the uncertainty of each plant has to be analysed separately. In this work, as a simplifying assumption and based on discussions with the site experts, it is assumed that the steam demand of each plant has a ±15% deviation around its nominal value. The uncertainties of the overall site demand from each steam network are consequently calculated using Equations (43)-(45), by taking the steam demand of the production plants from the steam networks for a representative month as basis. The results are presented in Table 3.

Results
In this section, the results of the proposed approach for the industrial real world case study are presented for two instances. As mentioned in Section 4, the first stage of the optimisation is defined as the first 8 h into the future, where the steam demand is certain, and the second stage is defined as the 48 h after the first stage with an uncertain demand. In the presented instances the steam import from the energy provider of the site is limited to the 30 and 15 bar pressure headers and 5 bar steam cannot be imported. In order to prevent a very high number of variables and to have moderate solution times, the discretization is done on a 2 h basis, which results in 4 and 24 time steps for the first and the second stages in each moving horizon iteration.
Moreover, the scenario tree of the two-stage stochastic optimisation problem is defined by considering the demand uncertainty of only two of the three steam networks. This results in nine scenarios for the second stage of the optimisation problems. This simplification is justified due to the fact that at INEOS in Köln the demand of the 30 bar steam network has smaller values compared to the other networks. Moreover, the amount of 30 bar steam that can be imported from the energy provider of the site in case of a shortcoming (see Section 3.7.1) is higher than for the 15 bar steam network, which reduces the chance of an extreme shortcoming in the network. Therefore, the scenario tree of the two-stage stochastic problem is defined by taking only the uncertainties of the 5 and 15 bar steam networks into consideration. The scenario tree and the realisation of the uncertain steam demand are defined according to Section 4 and by defining deviations around the nominal demand of each pressure level of the steam according to Table 3. The analysis of the time-dependent scenario probabilities in Section 4.1 is assumed to be valid for the demand of all steam networks and the probabilities of the scenarios that depend on the current realisation are defined according to Table 2. These probabilities are updated at the beginning of each moving horizon iteration depending on the current realisation together with all of the other model parameters. Two instances are investigated in what follows: • Instance 1: As discussed in Section 1.2, an extreme shortcoming in the steam networks can incur substantial costs due to the resulting forced reduction in the utilisation of the production plants. The first instance intends to show the difference of the performance between the two-stage stochastic and the deterministic approach in preventing such situations from happening. • Instance 2: The second instance is defined to investigate the possible savings that can be realised by using stochastic optimisation in normal daily operation compared to using a deterministic approach.
The optimisation problem is formulated in AIMMS ® 4.60 and solved using the solver CPLEX ® 12.8 to 1% optimality gap. It has around ∼2 × 10 6 constraints and ∼9 × 10 5 variables, from which ∼2 × 10 5 are binaries. The optimisations were run on a 64-Bit-Windows 10 with an INTEL ® CORE ™ i5-6300U CPU clocked at 2.4 GHz and with 8 GB of RAM. Due to the confidentiality of data all of the results in this section are presented in scaled coordinates.

Instance 1
The instance is defined as a period where one of the boilers (B3) is unavailable. The uncertain steam demands of the different pressure headers for two consecutive moving horizon iterations are presented in Figure 10.   The sub-figures include the trajectories of the different scenarios of the uncertain demand, denoted with "Scn i", together with the weighted average demand of all scenarios w.r.t. their probabilities (expected value), denoted by "EV". The demand of the first stage, on the left side of the vertical separating line, is certain and equal for all scenarios. As the probability of the realisation of different scenarios is dependent on the current situation, the calculated expected value at each iteration of the moving horizon is adapted. In the first iteration, shown in Figure 10a, the realised uncertainty for all steam networks is in the low demand region. This fact is reflected in the expected demands of the second stage which is below the nominal demand, because the scenarios with low demand have a higher probability of realisation (see Table 2). The deterministic optimisation solves the problem by considering the expected value as the single possible demand scenario (the black dashed line), while the stochastic optimisation takes all of the scenarios (Scn1-9) into consideration. This difference results in a different decision for the operation mode of boiler B4 at time point 4. Figure 11 shows the different operating strategies that are proposed by the deterministic (EV) and the two-stage stochastic optimisations (2ST). The deterministic optimisation decides to turn this boiler off, which forces the boiler to stay off also in time step 5 due to the minimum stay time constraints (see Table 1). Based on the expected value trajectories, the power plant is able to compensate the shortcoming of the steam networks at this point of time by purchasing steam from the energy provider. In contrast, in the next optimisation iteration, presented in Figure 10b, the uncertainty of all of the steam networks is realised in the high demand region. Since B4 is off, this realisation results in a large shortcoming in the steam network at time step 5 which cannot be compensated due to the limited amount of the steam that can be purchased from the external energy provider. In case of such shortcomings in the real production site, the production plants must reduce their production which incurs a substantial financial loss. The two-stage stochastic optimisation is however aware of the possibility of such extreme changes and therefore decides to keep B4 on. Thus, implementing the solution of stochastic optimisation in the first moving horizon iteration prevents an extreme shortcoming in the future. This also shows the importance of considering a relative long optimisation horizon (of two days): the optimisation decides on the optimal operation taking into account the impact of the different modes of operation for the different scenarios over a significant future time horizon.

Instance 2
The comparison framework presented in Section 2.3 is used to evaluate the benefit of the stochastic optimisation over its deterministic counterpart. To simulate cases that represent the realisation of the uncertainties in accordance to Table 2, the following steps are taken:
The first realisation is defined for all of the three steam networks as the nominal demand region (α d,1,m ); 3.
The next realisation of each steam network is selected randomly based on the probabilities presented in Table 2; 4.
Step 3 is repeated until the realised demand regions of the steam networks for the last optimisation iteration (40th) have been defined. 5.
The distribution of the created cases is compared to the values in Table 2: (a) If the distribution was reproduced for all of the steam networks within a predefined tolerance ≤ 0.02, the definition of the uncertainty realisation instances is completed; If not, a new planning horizon of 40 optimisation iterations is added and steps 2-5 are repeated.
Following the above-mentioned steps, a total number of 240 (6 × 40) uncertainty realisation instances were defined.
The comparison of the results is done between the realised performances when using different solution approaches. The costs and the incomes which are considered are the values that are computed for the time points for which the decisions are actually implemented, i.e., the decisions that correspond to the first stage of the optimisation problems. Therefore, the actual amount of the steam that is sold to the production plants is identical for the compared solutions. Thus, differences occur only in the cost factors (purchased fuel, purchased steam, vented steam, penalty for extreme shortcomings) and in the amount of the electric power that the power plant has produced and sold to the production plants. Moreover, the operational decisions of course are different.
For the comparison of the results an Overall Difference Indicator (ODI app1,app2 ) for two solutions is defined in Equation (46) which provides a relative comparison of the solutions: ODI app1,app2 = ∑ t∈T (∑ k∈K (P t el,k,app1 − P t el,k,app2 ) · c el,t + (costs t app2 − costs t app1 )) ∑ t∈T (∑ k∈K P t el,k,app2 · c el,t − costs t app2 ) , where the indices app1 and app2 denote the variables related to solution approach −1 and −2 respectively. P t el,k,app i denotes the produced electric power by turbine k at time t for solution approach app i and costs t app i presents the sum of costs (see Equation (40)) at time t for the solution approach app i.
The first comparison is done between the realised performances when using the twostage stochastic optimisation approach (2ST) and the deterministic solution based on the expected value of the steam demands (EV). ODI 2ST,EV is calculated as 1.8%. This represents the overall expected savings that the stochastic solution can realise when compared to its deterministic counterpart. At first sight, this may seem to be a small improvement, but due to the size of the operation this represents a significant financial gain. Moreover, the power plant is a system that has been continuously improved in small steps, so additional gains are difficult to achieve without further investments in new equipment.
The factors that influence ODI 2ST,EV are presented in Figure 12. As none of the solutions caused an extreme shortcoming, this contribution has been left out. The 2ST approach has a higher overall purchasing cost of fuel and has an insignificantly higher amount of vented steam. In contrast, its expenditure for purchasing steam from the energy provider is lower and it realises a higher income from selling electric power to the production plants. As mentioned in the beginning of Section 5, in order to have moderate solution times, the scenario tree of the 2ST problem is defined by taking only the uncertainties of the 5 and 15 bar steam networks into consideration; the demand uncertainty of the 30 bar steam network is neglected and its nominal trajectory is used to define all of the nine scenarios (see Figure 10). However, the uncertainty of all of the three steam networks are included in the presented comparison framework to define the realisations of the uncertainties. As presented in Table 2, the realisation probability of the second-stage trajectories depend on the current situation. Therefore, as a pragmatic and a computationally cheap extension, it is investigated if the incorporation of this dependency into the calculation of a single trajectory for the second-stage of the 30 bar steam network can be advantageous. To do so, instead of using the nominal trajectory of the 30 bar steam network, the weighted average of the three possible demand trajectories w.r.t. their probabilities presented in Table 2 is used, analogous to the expected value which is calculated for the second-stage of the 15 and 5 bar steam networks in the EV problem. An example of the calculated steam demand trajectories for the 30 bar steam network using the two approaches, with (2ST * ) and without (2ST) considering the dependencies, is depicted in Figure 13. In this example, the steam demand trajectories of the second-stage of an optimisation instance are compared under two different assumptions for the current demand region: a high demand region (D h ) and a low demand region (D l ). Using the 2ST approach, the second-stage trajectory remains the same for both of the realisations. However, the 2ST * approach results in a different trajectory for each of the assumptions, 2ST * l or 2ST * h , since the probability of different possible scenarios is adapted based on the current demand region.
The optimisation problems are solved again and with the sole difference that the second-stage trajectory of the 30 bar steam network is calculated with the 2ST * approach as explained above. ODI 2ST * ,EV then results as 2.6% which indicates the overall savings that can be realised by considering the dependency of the demand trajectory of the second-stage of the optimisation problems on the current realisation for the 30 bar steam network.  Figure 14a presents the factors that influence ODI 2ST * ,EV . Similar to the previous example, as none of the solutions caused an extreme shortcoming, this contribution has been left out. The 2ST * approach has a slightly higher purchasing cost of fuels and has a marginally higher amount of vented steam. In contrast, the incurred cost from purchasing steam from the energy provider is lower and it has a higher income from providing the production plants with electric power. Figure 14b depicts the difference between the 2ST * and the 2ST solutions. The 2ST * approach has a lower overall purchasing cost of fuel and realises a higher income from selling electric power to the production plants. The expenditure of the 2ST * approach for importing steam from the energy provider of the site and the cost incurred due to venting steam are slightly higher.
The incorporation of the dependencies of the second-stage trajectory of the 30 bar steam network on the current demand region into the EV problem did not change the results significantly and the improvements were marginal. The average solution time of the optimisations for the moving horizon iterations is approximately 550 s per iteration of the two-stage stochastic approaches and is approximately 14 s per iteration for the deterministic approach.

Conclusions and Outlook
This contribution addresses the optimisation of the operation of an industrial CHP plant in the presence of uncertain steam demands of the production plants as it is the case in real industrial operations. An approach is proposed which solves two-stage stochastic optimisations consecutively on a moving horizon. The optimisation problem is formulated as an MILP which includes the mass and energy balances, the different operating modes and the mode transitions of certain elements of the CHP plant. In order to reduce the inaccuracies of the linear formulation, the model parameters are updated at the beginning of each iteration. Furthermore, it was shown that the probability distribution of the scenario tree of the two-stage stochastic problem depends on the current situation of the site. An approach is proposed for an improved definition of the scenario tree of the two-stage stochastic optimisation based on the historical data. This fact is reflected in the optimisation by updating the probability distribution of the scenarios of the steam demand based on the information about the current situation.
The proposed approach was applied to the power plant of INEOS in Köln. It is able to reduce the likelihood of extreme shortcomings in the steam network which can force the production plants to reduce their utilisation. Furthermore, under normal operation, it can reduce the operational costs of the power plant significantly when compared to a deterministic solution.
A possible improvement for future work is to strengthen the reactive element of the presented optimisation approach by recomputing the optimisation on shorter intervals, e.g., 2 h. This would result in an earlier reaction to the changes in the state of the power plant and in the demands from the different steam networks.
Another aspect that should be considered in future work is the effect of the CO 2 emissions on the planning. As a measure to reduce the emission of the green-house gases to tackle the challenge of climate change, the companies in Europe are obliged to compensate the CO 2 emission caused due to the activity of a certain group of their processes, e.g., power plants, through purchasing the so-called CO 2 -certificates [39]. Meanwhile the price of these certificates has increased significantly, which can affect the optimal planning of the operation of power plants.
Another aspect for further research is an investigation on the possible savings when the uncertainty of other parameters, e.g., electricity price, is considered in addition to the steam demand uncertainty.

Conflicts of Interest:
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.