Efficient Sensitivity-Based Cooperation Concept for Hierarchical Multilayer Process Automation of Steam-Powered Plants

Hierarchical optimization architectures are typically employed to manage industrial steam processes efficiently. The key challenge today is to find a near global optimum despite these subdivided automation structures. This paper proposes a novel cooperation concept between hierarchical layers. In the upper layer, optimal static setpoints are computed for economic operation schedules of the entire plant. In the lower layer, a model predictive controller realizes these schedules, regulates the dynamic plant parts, and treats occurring disturbances optimally in a stochastic manner. To adequately overcome limitations due to mismatches between the optimization layers (caused by model mismatches or disturbances), it is essential to establish an efficient cooperation concept. Therefore, the presented concept exchanges specific sensitivity information between the optimization layers, where it is skillfully exploited to the benefit of global automation objectives, resulting in optimal expected operating costs. The novel concept is demonstrated via simulation studies calibrated with industrial measurements of a chipboard manufacturer. A performance analysis shows that the proposed cooperation concept outperforms alternative approaches, leading to the best possible trade-off between additional fuel and steam demand violation costs.


I. INTRODUCTION A. MOTIVATION
The International Energy Agency (IEA) annually publishes the World Energy Outlook, in which various scenarios provide an overview of tomorrow's energy sector [21]. Of particular interest are two scenarios that envisage compliance with the Paris climate targets [45], or even more ambitious zero net emissions by 2050. Both scenarios rely on renewable energy sources.
The goal of achieving a climate-neutral energy sector through the integration of renewable energies requires greater flexibility on the supply and demand sides of industrial processes while also increasing their efficiency. In energyintensive industries, as in the chipboard manufacturing plant The associate editor coordinating the review of this manuscript and approving it for publication was Ton Duc Do. under investigation, steam plays a major role in converting heat to power or distributing thermal energy throughout the plant to drive its mechanical and chemical processes. Additionally, industrial plants are facing new challenges in the transition to climate-friendly energy systems, such as the task of feeding power into electricity grids at peak times or supplying district heating to their neighborhoods. These challenges require the use of automation methods (e.g., advanced process control methods and operation optimization) to operate industrial steam processes efficiently and sustainably.

B. INDUSTRIAL PROCESS AUTOMATION
A brief overview of advanced process control and its role in the industry is given in [20], [50]. Basic process controls, such as classical PID-based control methods, are designed and implemented within the process components themselves to provide the foundation for operation as well as automation. Advanced process controls are usually added at a higher level, often later, to exploit optimization potential by combining process knowledge with control techniques.
Numerous control methods, from basic to advanced schemes, have been applied in literature for the control of steam-powered systems. A review using both conventional and advanced methods was conducted by [51] and includes academic studies as well as industrial practice. A well-known advanced method is the use of a model predictive controller (MPC), see, e.g., [38]. The MPC computes optimal control actions over a prediction horizon in order to regulate a system as close as possible to optimal setpoints while taking constraints into account.
In addition to process control, operation optimization determines optimal strategies for operating industrial processes while taking advantage of their economic characteristics. The optimal commitment-which defines if a unit is on or off-and the loading level of relevant units of the process need to be determined. This optimization problem is typically called Unit Commitment (UC) problem and originated in 1949 in power systems research [2]. Since then, the UC problem has been the focus of many studies and is also widely applied to thermal processes [1].
To solve the UC problem, a variety of methods and algorithms have been proposed. For example, they include priority listing, dynamic programming, Lagrangian relaxation, simulated annealing, fuzzy systems, artificial neural networks, genetic algorithms, and integer and linear programming. Out of these methods, mixed-integer linear programming (MILP) is an operation research method in which specific variables are integers [1]. This solution technique yields feasible solutions, however, large optimization problems lead to high computational complexity and long run times [46]. Nevertheless, the MILP formulation has become one of the state-of-the-art solutions for the UC problem and is still widely applied.
A concept that goes one step further than traditional operation optimization approaches is the Energy Hub (EH), introduced in [12] and reviewed by [40]. The EH is a concept for optimal energy management of multi-energy-carrier systems, including electrical and thermal energy, water, and gas. It describes a system where multiple energy carriers can be transformed, conditioned, and stored. Thus, this concept emphasized a global view of energy systems, not taking into account one energy carrier in one process, but multiple energy carriers, possibly in a network of energy systems [32]. Although the EH concept is more far-reaching than traditional optimization approaches, it can basically be seen as an extension of the UC problem. Thus, for the formulation and solution of EHs, also methods and algorithms of the UC problem can be used, MILP amongst others [17].
The question arising is whether and how process control and operation optimization can be efficiently combined.

C. HIERARCHICAL PROCESS AUTOMATION ARCHITECTURES
In order to operate industrial steam processes efficiently and combine process control and operation optimization strategies, hierarchical process automation architectures are typically employed [11], [35], [41], [43]. In these concepts, hierarchical process automation is decomposed into several parts (layers) to abstract the complex optimization problem into manageable subproblems. The hierarchical process automation problem is usually decomposed according to the different dynamics of states and disturbances. The sampling time and optimization horizon decreases in the lower layers due to faster process dynamics, disturbances and real-time computation requirements. Typically, real-time optimization (RTO) in the upper layer exploits the economic characteristics of the process while accounting for significant changes in operating conditions, such as price factors and slowly varying disturbances that strongly affect plant profitability. It provides either stationary setpoints or dynamic trajectories as well as additional information on economically optimal operation to the control layer beneath. The control layer for ongoing operations applies the information provided to match the output variables as closely as possible to the optimal setpoints/trajectories in the presence of disturbances or model mismatches [43], [48]. The following RTO/MPC strategies can be distinguished based on their temporal characteristics.
The static RTO/dynamic MPC strategy uses a steady-state plant-wide model for its economic optimization in the upper layer and dynamic models of particular sub-systems for control in the lower layer. The static RTO in the upper layer can only provide steady-state setpoints and is thus less frequently executed. In general, this approach does not consider the economic performance of the plant during the transition from one steady state to another, leading to suboptimal operation for processes that have slow dynamics, experience frequent transitions, or are continuously affected by fast disturbances [22], [43]. The plant's economic performance can suffer since the RTO and MPC solve two different optimization problems, each with different models, objectives, and on different time scales, while accounting for disturbances differently. The conflict between these optimization layers can lead to both infeasibility and unreachability of the economic setpoints in the control layer, resulting in poor economic performance. Therefore, it is essential and inevitable to check the consistency between the different layers even despite their different levels of abstraction. Additionally, an appropriate strategy to unify the (possibly competing) objectives is highly desirable [19].
To address this shortcoming, a new optimization layer was added between the RTO and MPC referred to as steadystate target optimizer [19], [22]. This layer computes feasible setpoint updates for the MPC to which the system must be stabilized, using a steady-state model that is consistent in its formulation as well as frequency with the dynamic MPC, and taking into account information from the RTO [13], [29], [34]. However, this approach is limited with respect to transient operations or fast disturbances.
In other approaches, instead of RTO based on steadystate models, dynamic real-time optimization (DRTO) is performed, which takes into account disturbances and process dynamics and provides target trajectories (rather than steady-state setpoints) to the lower layer [6], [24]. Tosukhowong et al. [44] propose a DRTO approach for the upper optimization layer, which is designed to capture only the dominant slow modes that describe the effective plantwide dynamics but are not affected by local fast disturbances. Kadam et al. [23] decompose the operation of chemical processes into a two-layer architecture, where in the upper layer a DRTO optimizes the dynamic trajectory, which is tracked by the controller of the lower layer. Instead of performing the trajectory optimization at a fixed interval, the DRTO computation in the upper layer is performed only when external disturbances exceed a certain limit. Würth et al. [52] implement a neighboring-extremal control strategy in the lower layer control problem, designed to update the economically optimal control trajectory based on sensitivity information computed at the DRTO level. Vega et al. [48] present a benchmarking of different hierarchical control structures consisting of RTO, DRTO and MPC for wastewater treatment plants. Jamaludin and Swartz [22] developed a closedloop DRTO-MPC formulation along with approximations thereof to manage the computational load and tested their approaches on a polystyrene reactor. Ramesh et al. [36] proposed a DRTO-MPC strategy to handle systems around unstable operating points and tested it on a continuous stirred tank reactor use case. Other strategies aim to approximate the RTO problem setting in the lower-layer MPC, which is considered as economic MPC (EMPC) [3], [10], [37]. Also combinations of RTO/EMPC strategies exist [53].
Nevertheless, strategies such as DRTO-MPC remain computationally expensive, and the use of plant-wide implementation for many different interacting plant components with many states and a large optimization horizon does not appear reasonable as it is too complex and too demanding for real-time requirements. For a reasonable implementation of hierarchical automation architectures, though, efficient cooperation between optimization levels is essential to adequately account for the effects of modeling mismatches between layers (static and dynamic models), disturbances, or plant constraints. In this work, a sensitivity-based stochastic method is developed to do so in a highly efficient manner.

D. MULTI-LAYER OPERATION WITH SENSITIVITY-BASED COOPERATION CONCEPT
In the investigated use case, the main process is steam generation that supplies heat to a chipboard production. In addition, a steam turbine generates electricity, and heat can be fed into a district heating network that supplies a hospital. An implementation of an DRTO-MPC or an economic MPC approach for optimal process operation and control would not be reasonable due to computational complexity. Therefore, a two-layer operation optimization and control architecture for industrial steam process is developed, including a wellsuited cooperation concept.
The operation optimization is formulated as MILP-UC problem based on the EH concept, forming the upper process automation layer. In this automation layer, the entire chipboard manufacturing plant is optimized, considering heat, electricity, interactions with the district heating network as well as the electrical grid, and economic criteria. In addition to the traditionally used EH concept, the product streams and production units of the chipboard manufacturing plant are included in the optimization, according to [16], [17]. To this end, this upper automation layer is further referred to as energy hub optimization (EHO) and is built on static models.
An MPC is chosen for the lower automation layer to track the setpoints computed by the EHO and control the system. The MPC employs a nonlinear dynamic plant model of the crucial parts which is successively linearized around the current trajectory, see, e.g., [25], [54]. In order to account for disturbances, stochastic MPC methods as in [18] are adopted.
In the proposed hierarchical operation optimization and control architecture, the EHO is introduced as an offline economic optimization scheduling problem, while the controller (MPC) is an online optimization problem. Due to stochastic fluctuations in the combustion process of the chipboard manufacturing plant, the steam demand computed by the EHO cannot always be met exactly. There are several ways to handle the resulting problem and bridge the time domain difference between the two automation layers. One option would be to calculate the action integral over time to ensure that the total amount of energy demanded by the EHO is met by the MPC over a given time horizon. In this work, a different approach is taken, which ensures that the energy demand is met according to its importance at a given point in time. In order to have sufficient steam available at all times, excess steam must be generated and thus more fuel must be consumed. An efficient cooperation concept between the EHO and the MPC is developed to optimally solve the tradeoff between fuel consumption and steam demand satisfaction, in which stochastic information and true costs are optimally exploited. Thus, a complex coupled plant operation can be better managed in terms of overall costs and also emissions.
To do so, the EHO passes the (time-varying) sensitivities of its economic cost function to a cooperation concept. In the cooperation concept, optimal expected operating costs (risk costs) determine permitted violation probabilities of the EHO's setpoints. The MPC uses the desired setpoints of the EHO and their permitted violation probabilities computed by the cooperation concept to efficiently produce and distribute the steam over the entire plant, even without knowledge of all plant-wide components and the global economic optimization goal. The occurring process disturbances are handled by the MPC in a stochastic manner, and their corresponding covariances are estimated online. As a result of the sensitivity-based cooperation concept, the plant can be operated at optimal expected costs. To the best of the authors' knowledge, no such process automation framework for industrial steam processes has been presented so far.
The investigated use case is based on plant models that are calibrated with real industrial measurement data. The simulation studies are performed with these validated models.

E. MAIN CONTRIBUTIONS
The main contributions of this paper are as follows: • A novel multi-layer automation framework for industrial steam processes is proposed, which employs a hierarchical structure and leads to minimal expected operating costs (risk costs) while considering stochastic disturbances.
• A developed cooperation strategy between two optimization layers ensures efficient cooperation between the high-layer EHO and the low-layer MPC by computing a permitted violation probability of the EHO's setpoints based on its sensitivity and stochastic information.
• The underlying static and dynamic models are presented and their parameters are identified from real industrial measurement data of a chipboard manufacturing plant.

F. PAPER STRUCTURE
The paper is organized as follows: Section II describes the use case and the purpose of the different optimization layers in more detail, followed by the description of the static model and the EHO approach in Section III. The dynamic models of the plant are introduced in Section IV, and their parameters are identified from the chipboard manufacturer's data using parameter sensitivity analysis. Section V describes the control concept for the chipboard production process, including the nonlinear stochastic MPC, the state observer, and the disturbance covariance estimation. The cooperation concept between the EHO and the MPC is stated in Section VI. In Section VII, simulation studies to demonstrate the effectiveness of the proposed framework are performed and a conclusion is drawn in Section VIII.

II. USE CASE A. PLANT DESIGN
The investigated plant of the chipboard manufacturer is operated with steam, see Figure 1. First, the steam is generated and superheated. From the steam generation unit, the steam is distributed throughout the plant via two extraction ports. The steam at the lower enthalpy level is used in a heated press to harden the adhesives in the chipboards. The steam at the higher enthalpy level is expanded in a turbine and generates electricity. The turbine has two extraction ports at different pressure levels. At the intermediate-pressure stage, the steam can be extracted and supplied to a chip dryer or to district heating. At the low-pressure stage, the steam can only be used for district heating (if its temperature level is high enough, depending on the season and turbine settings), or is condensed and returned to the plant's water cycle. The plant is modeled in terms of dynamics at two complexity levels: • The static part consists of models for all components in the entire plant: steam producer, steam turbine, presses, dryers, district heating heat exchangers, air condenser (after the turbine), and two additional smaller steam producers that are only used if required. The models are simplified, do not contain any dynamics, and represent the main behavior as well as couplings of the plant. The static component models are implemented in an EHO to compute cost-optimal setpoints.
• The dynamic part decomposes the steam generation process into multiple detailed (sub)models. These dynamic models are implemented in a controller (MPC) to follow the setpoints computed by the EHO.

B. PLANT OPERATION
The EHO computes optimal static setpoints for the entire plant, exploiting its economic characteristics while accounting for significant changes in operating conditions such as price factors or other slowly varying disturbances. The MPC stabilizes the dynamic plant components and follows the EHO's cost-optimal setpoints by minimizing a quadratic cost function that describes the deviations of the states from their desired values.
In a trivial scenario, where no disturbances occur and the operating conditions do not change, the MPC realizes the setpoints of the EHO and performs optimally. However, during a realistic plant operation situation, disturbances and time-varying operating conditions occur. In such situations, an MPC that only tracks the static setpoints with a quadratic penalty on deviations yields suboptimal economic efficiency.
To achieve a (near) optimal solution, a cooperation concept is developed based on the sensitivities of the economic cost function and stochastic information about the disturbances that occur. The proposed cooperation concept manages the interaction between the EHO and the MPC most efficiently in the presence of disturbances.

III. ENERGY HUB OPTIMIZATION
The upper optimization layer (EHO) is based on the previous work [17], where we developed a MILP optimization framework for industrial (manufacturing) processes and applied it to the chipboard manufacturing plant. Using the MILP optimization framework, relevant units of the use-case and their interactions are formulated by linear constraints and objectives, and the models are aggregated to form the final static plant model. This static plant model includes all relevant units of the chipboard production plant, in contrast to the dynamic control design model that only represents the steam generation process.
After the formulation of the static plant model, the optimization of the plant can be carried out, resulting in optimal trajectories of all units' energy flows. These trajectories are the final results of the EHO layer and specify the setpoints to be followed by the lower control layer. Because the control concept only considers the steam generation process, the relevant setpoints define the steam flows for the press and for the steam turbine, see Figure 1.
As the entire formulation and description of the EHO of the chipboard production plant would exceed the scope of this work, we refer to [17] for a detailed description. However, in this Section, we briefly summarize the used constraints and objectives, the modeled units of the static plant model, and relevant optimization features.

A. OPTIMIZATION PROBLEM
The static plant model of the chipboard manufacturing plant is based on a cost-based UC approach using three binary variables to describe the behavior of one unit, according to [8]. Every unit can have an arbitrary number of flows (energy or product flows) entering or leaving. Using the following constraints, a unit itself and its corresponding flows can be described: • Conversion constraints describe the linear relation between different flows of a device, e.g., between input and output. If required, nonlinear behavior can be approximated by piecewise linear segments, using additional binary variables.
• Start-up constraints define that a unit can either be on or off, and can not be started up and shut down simultaneously.
• Maximum/minimum generation constraints set a maximum/minimum value for a flow of a unit.
• Ramp-up/ramp-down constraints limit the maximal positive/negative gradient of a flow of a unit.
• Minimum up-/ downtime constraints define the minimal up-/downtime of a device after it was started/shut down.
• Storage constraints are implemented to model the use of (limited) storage.
• External requirement constraints define a desired fulfillment of external requirements at every time-step, e.g., district heating demand.
• Production schedule constraints also define a desired fulfillment, however, along a specified time interval. Thus, these constraints can be used to define a production schedule of a manufacturing plant. The EHO's optimization problem is based on an objective function J EHO , which is the sum of all fuel costs C EHO , rewards R EHO , and penalties E EHO over an optimization horizon H EHO (time step index k and sampling time T EHO,s ), The objective function depends on the plant operating schedule p p p = p p p k , . . . , p p p k+H EHO −1 , whose entries comprise the input and output flows of all plant components given by the vector p p p k . Rewards include selling electricity to the grid, while penalties are caused by failure to meet external requirements such as meeting district heating demand or production schedules [17]. By minimizing the objective function, subject to plant constraints as outlined above, the optimal plant operating scheduling is found.

B. STATIC PLANT MODEL
In contrast to the dynamic model that only includes the steam generation process with the medium steam, the static plant model contains all relevant parts of the chipboard production plant and the mediums steam, electricity, and product. In addition, two district heating demands need to be met by the plant operators, electricity can be sold to the grid, and a specified production schedule needs to be fulfilled. Thus, in the static plant model, the following units are modeled by the constraints outlined above: steam producer, steam turbine, presses, dryers, district heating heat exchangers, air condenser (after the turbine), and two additional smaller steam producers that are only used if required. The unit models are aggregated by connecting their input and output flows with corresponding other units' flows to form the entire static plant model. The static models were applied to process input data to validate the aggregated plant model, and their results were compared to process output data. The analysis showed that the aggregated static plant model can approximate the behavior of the plant well, showing an average error of only 2.4 % for all energy and product flows, see [17]. Nevertheless, this was tested over a time interval of one month and does not allow accurate conclusions to be drawn about the plant's dynamic short-term behavior.

C. IMPLEMENTATION
The optimization of the static plant model is conducted with the developed EHO framework using the YALMIP [27] toolbox in MATLAB R [30] with the solver GUROBI TM [15]. The horizon of the optimization is set to 12 h using hourly time-steps. The optimization identifies the optimal plant operating schedule to achieve minimal costs-which also includes the best possible fulfillment of the plant's external requirements. As a result, the optimal trajectories of all modeled units' flows are determined, including the info whether these units are on or off. Additionally, sensitivities of the EHO's cost function are exported for the cooperation concept. The trajectories can then be handed to the control concept, which aims to follow the setpoints based on dynamic component models, and the sensitivities are used for the control concept. 4

IV. DYNAMIC MODELING AND PARAMETER IDENTIFICATION
In this section, the dynamic models for the components of the steam generation process are presented and their assembly into a plant model is outlined. The implemented models are mechanistic and their parameters are identified and validated using actual industrial measurement data from a chipboard manufacturer. The models serve as simulated reality of the process and form the basis for the MPC.

A. GENERAL MODEL STRUCTURE
The general model structure of both the dynamic component models and the plant model is given by an implicit nonlinear state space descriptor form and reads wherein x x x ∈ R n x ×1 denotes the state vector, u u u ∈ R n u ×1 the input vector, z z z ∈ R n z ×1 the disturbance vector, Θ Θ Θ ∈ R n Θ ×1 the parameter vector, and M M M the mass matrix. The outputs of the model are given by with y y y ∈ R n y ×1 . In (3) and (4), the function f f f describes the differential equation's right-hand side, and g g g is the function mapping the states and inputs to the outputs. The mass matrix M M M and the function f f f depend on the parameter vector Θ Θ Θ ∈ R n Θ ×1 . This model forms the basis for all model variants described later, and the state meaning/composition remains the same throughout this work.

B. PLANT MODEL
The models are implemented in MATLAB R employing a class structure for each type of component and the assembled plant. The component classes construct the corresponding models and specify their properties, state equations, interface definitions, as well as connection rules to other components. The plant class contains the instances of the component classes as well as information concerning their assembly and is responsible for calling the components' methods. By calling the components' methods, the plant class computes their state equations, connects the components to each other through flow variables (e.g. steam, flue gas), and updates their state vectors. Thereby, the plant class has methods for simulation (via the MATLAB R solver ode15s), linearization, and control. The plant model is composed of the component models described in the following section.

C. COMPONENT MODELS
The main components of the steam generation are illustrated in the lower part of Figure 1. The furnace ''FN'' (symbolized by a flame) generates hot flue gas, which distributes thermal energy throughout the process. Thereby, water evaporates in the drum boiler ''DB'' and is then further heated in the superheaters ''SH''. Finally, the resulting steam is directed to the consumers ''C'' (press and steam turbine). The component models are mechanistic and derived from literature [4], [7]. Basically, there are two flows that drive the steam generation process and connect the components: flue gas and water/steam. The evaluation of these flows is executed through the plant model, which contains the assembled component models.
The flue gas is produced in the furnace and emits heat to the drum boiler ''DB'' and superheaters ''SH'', further also referred to as steam components ''SC''. The temperature of the flue gas at the outlet of the previous steam component corresponds to the inlet temperature of the flue gas of the following steam component.

1) FURNACE
The furnace ''FN'' combusts the fuel, e.g. wood, and releases thermal energy, which is distributed in the process by the flue gas. The initial flue gas temperature (combustion temperature) T FN and the heat capacity c p,FG is affected by the composition of the fuel, e.g. the moisture content of the wood, and thus subject to disturbances. The time response of the furnace is modeled assuming a linear PT1 behavior of the available power P FN , Therein, K FN,PT1 is the time constant of the PT1 behavior, P FN,nom is the nominal power and u FN ∈ [0, 1] the utilization factor of the furnace. The available power P FN determines the mass flowṁ FG of the flue gas ''FG'', As indicated by the chipboard manufacturer's measurement data, the heat flow from the flue gas ''FG'' to a steam component ''SC''Q SC arises primarly due to convection, Therein,α SC is the heat transfer coefficient that determines the heat transfer properties of the corresponding steam component. T SC,log stands for the mean logarithmic temperature difference which depends on the type of heat exchanger, either parallel such as for the drum boiler ''DB'', superheater ''SH2'' and ''SH3'' or countercurrent such as for the superheater ''SH1''. The mean logarithmic temperature difference reads T SC,max = T FGin − T SCin for parallel, T SC,min = T FGout − T SCout for parallel, The properties of the incoming and outgoing flows are denoted by the indices ''in'' and ''out'', respectively. In addition to (8), the energy balance for the flue gas must hold,Q SC =ṁ FG c p,FG (T FGin − T FGout ) .
2) DRUM BOILER In the drum boiler ''DB'', the water circulates in risers connected to the furnace. The water evaporates at saturation temperature and leaves the boiler as steam towards the superheaters. Meanwhile, feed water is constantly supplied to the drum boiler.
To model the dominating drum boiler dynamics, the wellknown and validated second-order model of [4] is adopted. It captures the overall behavior of the steam boiler well, but does not describe the behavior of the drum level. If the latter is of interest, other higher-order models from [4] could be implemented instead. In this work, the second-order model is employed.
The dynamic system equations are derived from the mass and energy balances of the drum boiler, ∂T sat ∂p DB .
The first row of (13) represents the transient mass balance for the water and steam fraction in the drum boiler. The second line results from the instationary energy balance of the drum boiler, for further details see [4]. The symbol V denotes the volume, p the pressure,ṁ the mass flow,Q the heat flow, h the specific enthalpy and ρ the density. The different phases steam and water are described by the indices ''s'' and ''w'', respectively. The influence of the drum boiler's metal parts in (13) depends on its total mass m DB,m (drum, tubes, and risers), the metal's specific heat capacity c DB,m , and the metal temperature, which is strongly correlated with the saturation temperature T sat . The input to the system (13) is the valve position u DB which defines the incoming feed wateṙ wherein K DB,w is the flow constant and p DB,fw the pressure of the feed water supply. The connection variables to other components are the heat flow to the drum boilerQ DB evaluated in (8) and the steam flowṁ DBout,s to the component downstream (5).

3) SUPERHEATER AND SPRAY ATTEMPERATOR
In the superheaters ''SH'', the steam supplied by the drum boiler absorbs the heat from the combustion gases of the furnace. Usually, there are multiple superheater stages in which the steam is brought to the desired temperature and is then delivered for further use, e.g. to a steam turbine. Spray attemperators ''SA'' support the precise control of the steam temperature and are integrated in the corresponding superheater model. The state equations are derived by means of the energy and mass balances [7], (15), as shown at the bottom of the page, with Therein V denotes the total volume, p the pressure, T the temperature,ṁ the mass flow,Q the heat flow, h the specific enthalpy and ρ the density. The different phases steam and water are described by the indices ''s'' and ''w'', respectively. The incoming and outgoing quantities are indicated by the index ''in'' and ''out'', respectively. The input to the system (15) is the valve position u SH,w ∈ [0, 1] of the spray attemperator which determines the incoming spray mass floẇ wherein K SH,w is the flow constant and p SH,w the pressure of the spray attemperator. The connection variables to other components are the heat flow to the superheaterQ SH evaluated in (8), and the incoming and outgoing steam flowṡ m SHin,s andṁ SHout,s (5), respectively.

4) CONSUMER
The consumer model ''C'' represents the press and the steam turbine. It is assumed as static and its properties determine the steam mass flow out of the steam generation process, wherein p SC denotes the pressure in the steam component ''SC'', K C is the flow constant, and p C the pressure of the consumer ''C''. The input to the consumer model is the valve opening u C ∈ [0, 1]. The required steam mass flow to the consumers is determined by the EHO, and the MPC attempts to realize the EHO's setpoints.

D. PARAMETER IDENTIFICATION
The parameters of the models described in Section IV-C are obtained either from data sheets or through parameter estimation of the chipboard manufacturing plant. For the latter, parameter sensitvity analysis is applied utilizing real industrial measurements from the chipboard manufacturer. A viable summary of the steps required for parameter estimation can be found in [39].

1) PARAMETER ESTIMATION
To estimate the unknown parameters not found in data sheets, a suitable quadratic cost function is defined, Here, Tr (· · ·) stands for the trace of a square matrix, representing the sum of its diagonal elements. The vector Θ Θ Θ contains all model parameters. The matrixỸ Y Y ∈ R m×n y contains the n y measurement channels at m measurement samples.  The corresponding simulated model outputs are in the matrix Y Y Y (Θ Θ Θ) ∈ R m×n y . The residuals of the different measurement types are scaled by the weighting matrix Q Q Q y . The optimal parameter estimates are given bŷ The parameter estimations are performed separately for each component, treating them in either a static or a dynamic optimization problem. The static parameter estimation of the flow constants in (5), (14), (16) and (17) yields a standard least-squares regression problem.
In order to solve the dynamic parameter estimation problem, the dynamic component models are simulated forward in time. The optimal parameter estimatesΘ Θ Θ are found using the solver fmincon of MATLAB R and minimizing the cost function J id (Θ Θ Θ) (18). To ensure good quality and significance of the parameter estimates, parameter sensitivity analysis is performed.

2) PARAMETER SENSITIVITY ANALYSIS
Many works have dealt with parameter estimation along with the related problems of optimal experiment design [5], identifiability analysis [14], [31], and measures of estimation quality [47]. In this work, the approach described by [39] for characterizing estimation quality as well as significance is recapitulated and applied.
Given a consistent estimator and the correct model structure, the parameter estimates approach their true values as the number of measurement samples grows,Θ Θ Θ → E (Θ Θ Θ). Consequently, the entries of the covariance matrix approach zero as well, Σ Σ ΣΘ = cov Θ Θ Θ → 0 0 0. Hence, in order to assess the accuracy of the parameter estimatesΘ Θ Θ, the associated parameter covariance matrix, is evaluated. Thereby, E [· · ·] stands for the expected value, and thus E [Θ Θ Θ] represents the actual parameter. In general, however, a numerical approximation directly of equation (20) is not feasible, but a lower bound can be obtained from the Cramér-Rao inequality [9], The Fisher information matrix F F F in (21) is defined as [26], [55]: In (22) Σ Σ Σ e denotes the error covariance matrix and φ φ φ (k, Θ Θ Θ) is seen as the parametric output sensitivity of the n p parameters, The Fisher information matrix is formed by applying finite differences at the optimal parameter estimatesΘ Θ Θ. In order to render the Fisher information matrix comparable for different types of measurements, it is scaled using a diagonal matrixΩ Ω Ω containing the optimal parameter estimates, If F F F s,ii → 0 the m measurement samples of the n y measurement channels do not provide enough information to estimate the parameter Θ i , i ∈ [1, . . . , n p ].

3) PARAMETERS AND MODEL FIT
All parameters and their corresponding scaled Fisher information of the dynamic models are listed in Table 1. The Fisher information indicates an excellent suitability of the available data for estimating the parameters not given by data sheets.
The geometry and material parameters of the drum boiler and the superheaters, namely V DB,t , m DB,t , c DB,t , V SH1 , V SH2 , and V SH3 , are obtained from data sheets. The time constant of the furnace is derived from a step response recorded in the chipboard manufacture.
The flow constants K of the steam flow between two consecutive components (5), the feedwater to the drum boiler (14), the spray attemperator (16) and the outlet (17) are obtained from static parameter estimation.
The other parameters are estimated using dynamic simulations, (18), and MATLAB R 's fmincon solver.
Fits of the model outputs using the optimal parameter estimates are shown in Figure 2. The measurements are historical data from a steam-powered industrial plant for the manufacture of chipboards. The R 2 -fit is 85.4%, 74.0%, 86.2%, and 1% for the models of the drum boiler, the superheater 1, the superheater 2, and the superheater 3, respectively. Even though the R 2 -fit for the superheater 3 model is low, a close examination of the model outputs ( Figure 2) shows that the transient and steady-state behavior of all models is sufficient to be implemented in an MPC.

V. CONTROL CONCEPT
The lower automation layer controls the steam generation process subject to dynamics and disturbances. The control concept is composed of a stochastic MPC, a nonlinear plant model, a successively linearized plant model, and an observer for state and disturbance covariance estimation, see Figure 3. It computes optimal control actions to fulfill a setpoint sequence r r r with a given permitted violation probability β. The setpoint sequence is a section of the optimal operating schedule calculated by the EHO, r r r ∈ p p p. The permitted violation probability β is obtained by the cooperation concept. In the cooperation concept, stochastic information about the disturbances and the EHO's cost function sensitivities are exploited to compute the permitted violation probability β that leads to optimal expected plant operating costs (risk costs).
In the control concept, the plant model (nonlinear) predicts the future state trajectories online using the current control inputs but assuming no disturbances. Then, the plant model is successively linearized around this trajectory and discretized in time. The stochastic MPC implements the obtained linearized plant model and computes optimal control inputs based on chance constraints to follow desired reference setpoint sequence from the EHO in a probabilistic manner. In order to provide the stochastic MPC with up-todate information, the observer estimates the true states from measurements of a plant (simulated reality) and computes their estimation error covariances. Additionally, the covariance of the disturbance is estimated.
In the following Section, k is the time step index with a chosen time step length t, so t k = k t, and signals with subscript k denote their values at time t = t k = k t.

A. PLANT (SIMULATED REALITY)
The simulated reality is based on the model (3)-(4) from Section IV calibrated against real measurement data. The process equations with disturbances and the measurement output equations with noise have the form with x x x p ∈ R n x ×1 and y y y p ∈ R n y ×1 . The disturbances z z z k and the measurement noise v v v k are both assumed Gaussian and white. The simulated reality describes the steam generation process subject to disturbances.

B. PLANT MODEL (NONLINEAR)
The nonlinear plant model matches the plant (simulated reality) (25)- (26), except no disturbances are assumed, with x x x m ∈ R n x ×1 and y y y m ∈ R n y ×1 . The plant model predicts the future state trajectories based on the current states and inputs (x x x m (t k+1 ) , u u u k ). The predicted states form the basis for linearization and state estimation.

C. LINEARIZED PLANT MODEL
The plant model equations (27)- (28) are discretized and successively linearized (see, e.g., [54]) around the predicted trajectories, where indicates the deviations from the state trajectories.
The system matrices A A A ∈ R n x ×n x , B B B ∈ R n x ×n u , C C C ∈ R n y ×n x , D D D ∈ R n y ×n u , and E E E ∈ R n x ×n z are derived using finite differences and their time index is omitted for reasons of brevity. The system (29)- (30) can be implemented in a standard MPC or an extended Kalman filter (observer) formulation.

D. OBSERVER
In the control concept, observers are used to estimate the true states, their estimation error covariances, and the disturbance covariance from measurements.

1) EXTENDED KALMAN FILTER
A standard extended Kalman filter (EKF), as described in [42], is implemented to estimate the true states and their estimation error covariances. The state update, x x x k+1 = x x x m,k+1 + K K K K y y y p,k+1 − y y y m,k+1 , is determined using measurements y y y p from the simulated reality and the predicted nonlinear plant model outputs y y y m . The Kalman gain K K K K is computed using the Jacobians of (27)- (28), the evolution of the estimation error covariance Σ Σ Σx ,k , and the process as well as measurement noise. The process noise is obtained by estimating the disturbance covari-anceΣ Σ Σ z,k as described in Section V-D2 and the measurement covariance R R R is assumed to be known.

2) DISTURBANCE COVARIANCE ESTIMATION
The disturbance covariance estimation is based on the EKF's assumptions and the definition of the state covariance matrix. First, the differences of the measurements and the model outputs are computed, s s s k+1 = y y y p,k+1 − y y y m,k+1 .
This current ouput residual is collected with n s older residuals in the innovation matrix, with S S S k ∈ R n y ×n s whereby n s is the window size of the estimator. Then, the covariance matrix of the measurements is defined asΣ The disturbance variance's influence on the states is given bŷ and the disturbance error covariance results in where (·)\ and /(·) denote left-and right-pseudo inverses of matrix (·). The quality of estimation significantly depends on the window size n s . A large window size is required for good convergence, whereas a small window size is better suited to account for a (slowly) time-varying disturbance covariance and thus reduces the effective bias in the estimation. A recursive estimation algorithm for disturbance covariance estimation can be found in [28], which, however, cannot consider a time-varying disturbance covariance.

E. CHANCE-CONSTRAINED MODEL PREDICTIVE CONTROL
MPCs are widely used to control multivariable systems under constraints. The key idea is to compute an optimal trajectory of the controlled system model, predicted over a chosen horizon from the current time into the future. The first control signal value is applied and the whole optimization is repeated at the next time step, based on an updated initial state. In order to account for disturbances and ensure sufficient steam quality, stochastic MPC methods as described in [18] are adapted here. The formulations of [18] result in the following deterministic optimal control problem of the linearized plant model (29)- (30), u u u k+i ≤ u u u con , x Therein, the cost function J MPC,N is minimized over the prediction horizon N to compute an optimal input control sequence u u u = {u u u k , . . . , u u u k+N −1 }. The first input u u u k of the control sequence is applied to the system, and then the optimization problem is recomputed at the next time step. In (37)-(44) the subscript k +i+1 denotes the value at the corresponding future time step based on the current knowledge. The expected deterministic value of the predicted state and output vector (no disturbances assumed) is given by x x x k+i+1 and ỹ y y k+i+1 , respectively. Here, indicates deviations from the trajectory computed by the nonlinear plant model, so that the effective predicted state and output vector is the sum of the two.
When the disturbances z z z k in (29) are modeled as random variables, the system dynamics are described in a probalistic sense. In the chance constraint (40), Pr [·] stands for the probability of an event. The vectors r r r k+i+1 and u u u con contain the desired outputs (setpoints) and the input constraints, respectively. In (40), the probability of violating an output constraint j should not exceed a permitted probability β j ∈ (0, 1]. The chance constraints (40) can be expressed in terms of the expected deterministic predicted output variableỹ y y k+i+1,j and the inverse cumulative density function F −1 k+i+1,j for each output constraint j, y y y k+i+1,j ≥ r r r k+i+1, The complexity of computing (45) depends on the underlying probability distribution of the stochastic disturbances z z z k as well as the specific nonlinearity characteristics of the plant behavior. Assuming z z z k as Gaussian and the plant behavior to be approximately linear with respect to disturbance propagation over the prediction horizon N , the state and output covariance can be propagated over the prediction horizon, In (46), Σ Σ Σ z,k is the disturbance covariance matrix, VOLUME 10, 2022 Since the disturbance covariance matrix is unknown and needs to be estimated as described in Section V-D2, the predicted state vector and thus the output vector is modeled via the Student's t-distribution x x x k+i+1|k ∼ T x x x k+i+1|k , Σ Σ Σx ,k+i+1 , ν , whereby ν is the number of degrees of freedom of the Student's t-distribution.
The cumulative density function F and its inverse can be obtained with initial condition Σ Σ Σx ,k = Σ Σ Σx ,k andx x x k =x x x k , both computed by the EKF in Section V-D1.
The cost function J MPC,N in (37) is designed to track the desired setpoint sequence r r r = {r r r k+1 , . . . r r r k+N } with r r r i ∈ R n y ×1 , wherein Q Q Q y and Q Q Q u are weighting matrices. By minimizing equation (49), the MPC is able to optimally track the setpoints given by the EHO with a given permitted violation probability β. The permitted violation probability β of (40) is obtained from the cooperation concept in Section VI using the EHO's cost function sensitivities.

VI. COOPERATION CONCEPT
The basic idea of the following cooperation concept is to determine a permitted violation probability β of the desired setpoints r r r that leads to optimal expected plant operating costs (risk costs). Therefore, the EHO passes its cost function sensitivities to the cooperation concept. In the cooperation concept, the EHO's sensitivities are implemented in a risk cost function that takes into account the stochastic distribution of the process variables. Minimizing the risk cost function leads to an optimal permitted violation probability β, which can be adopted by the stochastic MPC.

A. RISK COST FUNCTION
The risk cost function is based on static considerations of the EHO's cost function and describes the additional costs incurred by non-ideal plant operation due to disturbances. It consists of two terms, one for fuel costs and one for violation costs, Therein, c fuel denotes the fuel costs per additional output unit and c viol the violation costs per missing output, obtained by deriving the cost function of the EHO with respect to the corresponding setpoint These cost factors are multiplied by the difference between the desired setpoint r and the expected deterministic outputỹ as well as by the expected violations ω.
The optimal permitted violation probability β sought for the stochastic MPC is determined by minimizing the risk cost function, The output of the controlled plant y will fulfill the EHO's setpoint r with the probability 1 − β, see Figure 4. The green area shows the overproduced steam and the red area shows the steam shortage caused by disturbances. The tradeoff between steam overproduction and steam shortage is optimally solved by the cooperation concept.

B. EXPECTED OUTPUT AND VIOLATIONS
The expected outputỹ and the expected violations ω depend on the choice of the permitted violation probability β, see (40) and Figure 5. The expected output and violations are determined using the probability density function f X , which is chosen to be Student's t-distributed. Then, the permitted violation probability is given by the cumulative density function and the expected output is the inverse of it. The expected violations, weighted by severity and frequency, are determined by Using these relations, the optimal permitted violation probability can be computed (52).

C. MULTIPLE OUTPUTS
In the case of multiple outputs, the corresponding permitted violation probabilities can be evaluated separately if there are no other constraints. If the maximum available power is limitted, for example due to plant restrictions, the optimization for determining the permitted violation probabilities for multiple outputs must be performed with respect to complementary constraints.

VII. DEMONSTRATION OF THE NOVEL PROCESS AUTOMATION FRAMEWORK
The performance of the developed cooperation concept (sensitivity-based) is demonstrated via simulation studies. First, the settings of the process automation framework are described. Then, the effectiveness of the stochastic control concept is shown. Finally, simulations are performed using the multi-layer control architecture with and without the cooperation concept.

A. SETTINGS 1) SETTINGS OF THE EHO
The economic plant operation is optimized via the EHO from Section III using the YALMIP [27] toolbox in MATLAB R [30] with the solver GUROBI TM [15]. The EHO's optimization horizon is set to 12 h using hourly timesteps, T EHO,s = 1 h. The optimization is performed under real industrial constraints for district heating demand, production schedules, the electricity market and cost factors.

2) SETTINGS OF THE CONTROL CONCEPT
The control concept from Section V follows the EHO's setpoints with the permitted violation probability β computed by the cooperation concept. The sampling time of the control concept is chosen to be T MPC,s = 10 s. The observer for disturbance estimation has a window size of n s = 3000 simulation steps. The MPC controls the water volume of the drum boiler, as well as the enthalpy flow to the steam turbine and the press. The weighting matrices of the MPC are chosen diagonal. The entries of the output weighting matrix Q Q Q y are 10 5 for the water volume and 10 1 for the demanded enthalpy flow. The entries of the input weighting matrix Q Q Q u are 10 4 for all steam/water valves and 10 6 for the utilization factor of the furnace. The degrees of freedom of the Student's t-distribution are ν = n s − 1. Nonlinear simulation results are obtained via the MATLAB R solver ode15s [30].
The optimization problems of the MPC are solved using MOSEK R [33] and YALMIP [27].

3) SETTINGS OF THE SIMULATION
The parameters of the simulated reality are chosen as described in Section IV and V-A. A stochastic disturbance with a standard deviation of σ z = 0.15 is applied to the  input of the furnace u FN ∈ [0, 1], which is caused by a fluctuating fuel quality (e.g., due to the unknown moisture content of wood).

B. DEMONSTRATION OF THE STOCHASTIC CONTROL CONCEPT
The performance of the stochastic control concept is demonstrated using a stepwise heat demand of the steam turbine and an chosen permitted violation probability of 10%.
In the case of a linear reference simulation, the system was able to accurately meet the steam demand (realized violation probability 9.8%). However, when considering the nonlinear plant, the constraints get violated in 17.3% of the cases. The deviations in the permitted violation probability occur due to nonlinear plant behavior and the successive linearization. Therefore, the chance constraints (45) are modified by adding a nonlinear correction term nonlin , y y y k+i+1,j ≥ r r r k+i+1,j + F −1 k+i+1,j 1 − β j + nonlin . (55) The nonlinear correction term nonlin is obtained from online analyzing the last n nonlin = 300 samples. With the proposed correction term, a satisfying plant performance was achieved and the constraints were violated in 10.2% of the cases, see Figure 6. Tests with differently chosen permitted violation probabilities lead to similar satisfactory results.  Moreover, the observer performs well for the state estimate, and the disturbance covariance is estimated accurately, see Figure 7. The estimations quickly approache the actual values and converge.

C. DEMONSTRATION OF THE COOPERATION CONCEPT
The process automation framework is tested with three different control concepts in the lower layer: (1) a standard MPC corresponding to the stochastic MPC with a permitted violation probability of 50 %, (2) a stochastic MPC with a constant permitted violation probability of 3%, and (3) a stochastic MPC in combination with the cooperation concept resulting in an optimal permitted violation probability. The constant permitted violation probability of 3% is chosen to serve the overall optimum well and to strongly challenge the cooperation concept.
In the upper layer, the EHO optimizes the economic plant operation statically under real industrial constraints for district heating demand, production schedules, the electricity market, and cost factors. The EHO provides setpoints to the MPC and the cost function sensitivities to the cooperation concept. The disturbances act on the input of the furnace u FN caused by fluctuating fuel quality and are only compensated for in the MPC.
The simulation study is performed with 34 runs of randomly generated disturbances to account for the stochastic nature of the process. The simulation study is evaluated with respect to three different cost factors: (1) total additional plant operating costs, which can be broken down into (2) the steam demand violation costs and (3) the additional fuel costs. The additional plant operating costs result from the nonideal operation of the plant in the presence of disturbances. If no disturbances were present, the steam demand could be met exactly and the additional plant operating costs would be zero. Figure 8 shows the total additional plant operating costs caused by the occurrence of disturbances. The MPC with the novel cooperation concept is able to reduce the additional The additional plant operating costs can be divided into those arising from the violation of the required setpoints (steam demand) and those arising from the fuel costs for excess steam production. The steam demand violation costs are slightly higher for the cooperation concept compared to the stochastic MPC with a permitted violation probability of 3%, see Figure 9. The standard MPC case causes the highest steam demand violation costs. The additional fuel costs are minimal if no measures are taken to stochastically compensate for violations, as in the standard MPC case. Fuel costs are significantly higher using the stochastic MPC with the 3 % permitted violation probability than with the optimal permitted violation probability (cooperation concept), see Figure 10. Thus, in total, the additional plant operating costs are the lowest with the developed cooperation concept, optimally solving the trade-off between additional fuel costs and violation costs.

D. RESULTS
As mentioned in the introduction (Section I-C), the implementation of DRTO-MPC methods is challenging in terms of computational complexity [22], [52]. The concept presented in this paper aims to outperform classical RTO-MPC methods while keeping the computational cost lower than with DRTO-MPC methods. The proposed stochastic MPC with the sensitivity-based cooperation concept only slightly increases the computational effort compared to RTO-MPC methods, but is able to realize an increase in system performance of the same order of magnitude as with the DRTO-MPC methods of [48] in the present of disturbances. For future work, it would be interesting to compare different methods with the proposed cooperation concept using the same benchmark process.
In summary, the developed cooperation concept is able to achieve the cost-critical setpoints while efficiently limiting additional fuel consumption and thus emissions to the necessary level.

VIII. CONCLUSION
The developed process automation framework is able to operate an industrial plant at optimal expected costs in the presence of stochastic disturbances. The method is based on a hierarchical structure in which a higher layer (EHO) optimizes the economic costs and provides setpoints to a lower layer (MPC) that compensates for disturbances and considers dynamic performance. A novel cooperation concept between the optimization layers is proposed to overcome limitations due to mismatches caused by occuring disturbances. In the cooperation concept, the EHO's cost function sensitivities and stochastic information about the disturbances are exploited to compute an optimal permitted violation probability of the setpoints. The MPC tracks the desired setpoints with the permitted violation probability, resulting in an optimal trade-off between additional fuel and steam demand violation costs.
The operation optimization framework is tested using models calibrated with industrial measurement data and in simulation scenarios under real industrial constraints for district heating demand, production schedules, the electricity market, and cost factors. It is seen that the stochastic control framework can reliably realize setpoints with a given permitted violation probability. Observers are able to accurately estimate the true system states and the disturbance covariance from measurements. The cooperation concept between the hierarchical optimization layers leads to optimal plant operation in the presence of disturbances reducing the additional plant operating costs by 55.4 % compared to a standard MPC case. A stochastic control concept with a fixed permitted violation probability of 3% is outperformed by 11.1%. Using the cooperation concept, no further expert knowledge is required to control the plant, e.g., to adjust control parameters.
In real industrial processes, random disturbances always occur, caused for example by imperfections of fuel compositions or reactions. The process automation framework developed here can be applied to properly handle these kinds of disturbances and to operate the entire plant in an optimal and efficient way. A next step to improve the developed method lies in the data-driven modeling of the disturbances using generalized Gaussian mixture distributions, see [49]. This will help to better capture the probability distribution of disturbances such as fluctuating fuel quality and will make the framework even more significant for real-world implementation.
DOMINIK PERNSTEINER received the master's degree in mechanical engineering, in 2018, and the doctorate degree, in 2021. In his Diploma thesis, he worked on the numerical modeling of fluidized beds. In his Ph.D. thesis, he developed advanced modeling, control, and observation concepts for thermal energy systems. Since 2018, he has been working as a Research Assistant with the Research Unit for Control Engineering and Process Automation, Institute of Mechanics and Mechatronics, Vienna University of Technology (TU Wien).
VERENA HALMSCHLAGER received the master's degree in process and chemical engineering, in 2017, and the doctorate degree, in 2021. In her Diploma thesis, she worked on the mathematical description and experimental evaluation of carbon dioxide transfer in a minimal ivasive liquid lung. In her Ph.D. thesis, she developed an optimization framework and grey-box modeling concepts for industrial applications.