Steam Turbine Rotor Stress Control through Nonlinear Model Predictive Control

The current flexibility of the energy market requires operating steam turbines that have challenging operation requirements such as variable steam conditions and higher number of startups. This article proposes an advanced control system based on the Nonlinear Model Predictive Control (NMPC) technique, which allows to speed up the start-up of steam turbines and increase the energy produced while maintaining rotor stress as a constraint variable. A soft sensor for the online calculation of rotor stress is presented together with the steam turbine control logic. Then, we present how the computational cost of the controller was contained by reducing the order of the formulation of the optimization problem, adjusting the scheduling of the optimizer routine, and tuning the parameters of the controller itself. The performance of the control system has been compared with respect to the PI Controller architecture fed by the soft sensor results and with standard pre-calculated curves. The control architecture was evaluated in a simulation exploiting actual data from a Concentrated Solar Power Plant. The NMPC technique shows an increase in performance, with respect to the custom PI control application, and encouraging results.


Introduction
In the current energy market, renewable sources and flexible power plants are increasingly exploited in order to meet the energy demand of increasingly connected cities [1], where smart grids are going to become a consolidated reality [2,3]. In this context, the design of new power plants must aim at increasing the level of flexibility, with a resulting multitude of challenges in terms of degradation of the components, restrictions related to the environmental impact, and required level of workforce specialization. For these reasons, the design of Steam Turbines (STs) nowadays must take into account frequent discontinuous operations. For instance, in the case of Concentrated Solar Power Plants (CSPP), the inlet steam highly depends on the weather conditions [4], while in industrial power generation, the discontinuity is related to the scheduling of the electricity production, optimized according to the prices of energy media (natural gas and electricity) that vary on an hourly or daily basis. One of the most important aspects related to discontinuous operation is the need to enforce protection against stress induced by thermal transients, which, due to the continuous changes in the steam conditions, can lead to faster aging of the turbine and, in extreme cases, to component failure due to nucleation and propagation of cracks. While almost every ST component is affected by thermal stress aging, transient thermal stresses in rotors often represent the factor determining the minimum allowable start-up time and the maximum allowable production during start-up. The start-up optimization is carried out through several different steps. The first one usually takes place during the design phase of the machine itself, e.g., through a proper design of the rotor critical locations, which are mostly affected by peak stress. During the design phase, load curves are also optimized by exploiting accurate process models evaluating the specific plant components to which they are connected. The last phase for the start-up design and optimization is represented by the control logic design, which is typically implemented through Proportional Integral Derivative (PID) controllers. In more advanced control systems, which are described later on, the control approaches are based on predictive methodologies considering the machine thermal dynamics or on particular heuristics.
Thermal stress phenomena have been extensively studied in literature, mostly by focusing on their effects and on measurement procedures for early identification of aging phenomena. These are of crucial importance especially for critical components, such as rotors, where cracks after nucleation can rapidly propagate due to vibration phenomena that in normal conditions do not affect rotor integrity. In addition, low cycle fatigue damage cannot be directly detected with non-destructive methodologies. Indirect measurements of the most severe effects of thermal stress are often associated with vibrational phenomena on the machine axis. Most often, the effects of excessive machine stress can only be verified after dismantling the machine. Non-destructive tests are not always capable to highlight issues on the rotor and often cracks can be found only after the rotor is taken out of service.
The management of the machine thermal stress is a rather complex operation and, from the control point of view, is faced through a sequence of steps. In more recent years, these operations have been automated, thanks to a wider use of sensors, which have led to in-depth studies of their behavior through increasingly reliable models. An example of automated pre-warm is described by Badami et al. [5], in which a fuzzy logic-based supervisor has been designed to speed up the rotor warming phase, while taking into account the rotor stress conditions before the cold start of the machine.
The literature highlights two main complementary methods for the control of the rotor thermal stress, which concern either the limitation of the inlet steam flow-common in every steam plant-or the regulation of the temperature of the inlet steam, available only in specific plants. In Combined Cycle (CC) systems, the outlet steam temperature of the Heat Recovery Steam Generator (HRSG) can be adjusted by the direct control of the gas generator. However, this kind of operation shows the disadvantage of forcing the gas turbine to lower efficiency operating conditions. De-superheating of inlet steam by attemperators is also a possible solution for limiting the thermal stress in STs, at the cost of increasing installation costs. For these reasons, a common control strategy aims at adjusting the Heat Transfer Coefficient (HTC) or the inlet steam mass flow. As a first approximation, we can divide the work related to the start-up control into three categories: The first one refers to the modeling of the thermal behavior of the machines and/or connected systems, the second one is linked to the offline calculation of turbomachinery and related system load curves, while the third one concerns the startup online control with/without thermal stress related logic.
Some examples of works related to modelling topics are presented by Carnero et al. [6], which introduce a soft-sensor application aimed at monitoring the ST stress, by exploiting heat transfer equation to model the thermal dynamics in the rotor. Sun et al. in their paper [7] present a soft sensor that models the rotor temperature through transfer function and the damage due to thermal stress through polynomials. Rusin et al. [8] propose two modelling methodologies based on Green's function, Duhamel's integral and Neural Networks (NN). These approaches can be used to solve the inverse problem aimed at optimizing the operation of steam machines and are also exploited in a similar way by Banaszkiewicz et al. [9,10]. Dominiczak et al. [11] apply NNs for the prediction of temperature and stress in the ST critical points.
Offline optimization approaches have been proposed by Casella  Dymola [12]. Faille et al. [13] followed a similar approach with the objective of minimizing the start-up time.
A brief review of online control applications aimed at speeding up the startup of STs show three main classes of methodologies: (i) PID controllers, (ii) Artificial Intelligence (AI)-based techniques, and (iii) Model Predictive Control (MPC). Sun et al. [14] present a controller based on polynomials and PID aimed at extending the life of the machine, limiting its cumulative damage. Boulos et al. [15] propose three different control methods based on pure PID controller, Fuzzy-PID, and a Fuzzy Logic controller, applied to the startup of a CCPP by acting on the load rate. Matsumoto et al. [16] introduce a complex system based Fuzzy Logic and NNs to control the schedule of the power plant load curve as a function of the current status of the machines and the HRSG. MPC architectures have been widely studied for startup applications. Excellent examples are proposed by Nakai et al. [17], who describe a MPC based on linear programming and models based on heat conduction equations. Gallestey et al. [18] present a linear MPC strategy that integrates Fatigue Crack-Growth Models to the aim of extending the residual lifetime of the machine with consequent relevant economic benefits for a power plant application. Shirakawa et al. [19] describe a multi-objective MPC, which integrates several components of a power plant of 600 MW by including the thermal dynamics of the involved STs and by exploiting a complex system of NN models and thermal equations as well as a Genetic Algorithm (GA) for optimization. Larsson et al. extend the MPC application by using different nonlinear components models of a CCPP [20], with a simulation approach implemented in Modelica. D'Amato et al. [21,22], present several excellent works, in which MPC is applied to the startup of turbomachines, by exploiting heat equations and linearization procedures to implement an MPC strategy in Programmable Logic Controllers (PLCs). Rua et al. compared a linear and a nonlinear formulation of MPC to optimize the startup of a combined cycle while taking into account stress limits on the walls of the high-pressure steam drum and the steam turbine rotor [23].
Once the turbine has been correctly started, the management of the variability of the inlet steam conditions and the tracking of the power setpoint are often performed through PI control techniques, fuzzy PID [24], robust control [25], and in some cases through MPC techniques [26].
An in-depth analysis of the aforementioned works shows two possible methodologies for implementing an effective control system for the ST startup, PID controllers and MPCs. The former ones ensure good performance and usually lower modeling and implementation efforts, with some shortcomings related to the inability of the PID controllers to anticipate any control actions to avoid stress overshoots and problematic operating points. On the other hand, the MPC-based strategies naturally manage all the constraints of the machine and, thanks to the model-based approach, can predict its dynamic behavior. This capability makes the control actions effective by allowing anticipation and avoidance of numerous undesired phenomena. Although MPC appears the best-performing strategy, it shows some heavy limitations related to the online calculation of the control action, which requires calculation times that are often not compatible with PLC platforms.
The present research work is part of this context, as it aims at identifying and overcoming the most relevant bottlenecks in the application of MPC for controlling the ST power startup through a local point of view of the phenomenon of the thermal transient, in which the controller focuses on the steam turbine. In particular, this article presents a Nonlinear MPC (NMPC), which aims at maximizing the energy produced during the power startup of STs, by integrating a series of nonlinear models that describe the temperature dynamics in different sections of the rotor as well as the thermodynamic and mechanical behavior of the machine. The modelling scheme implements hybrid methods based on heat equations and empirical numerical models for predicting specific internal quantities. The NMPC optimizes the power load curve during the startup, through a tracking formulation of the optimization strategy, in which the maximum power reference is received from the Distributed Control System (DCS) of the plant. In particular, this article presents several novelties. In the first place, it compares the standard baseline curves for steam turbine start-up with the application of a PI architecture including the soft sensor calculating the rotor stress. Later this PI architecture is compared with the NLMPC application showing the benefits of this control procedure. It also has to be remarked that while other works aim at the application of MPC strategies to the whole plant, here the focus lies on the control of the single steam turbine. The application of these procedures to a CSP case through real data is of particular interest.
The paper is organized as follows. Section 2 describes the modeling methodologies and validation procedures for the models obtained. Section 3 provides a brief introduction on PID and NMPC control and describes the structure of each implemented controller and the associated technical solutions. Section 4 discusses the results of the simulation studies for the implemented controllers by comparing their performances and providing guidelines for tuning the NMPC control parameters. Section 5 provide final comments to the obtained results and some concluding remarks.

Steam Turbine Modelling Approach and Rotor Stress Monitoring
The first fundamental step to define an effective control strategy for optimizing the ST start-up is to design a tool capable of monitoring the thermal state of the machine. As previously mentioned, the stress can only be measured using particularly invasive hardware measuring instruments, which must be designed to work at rather high temperatures and pressures, with all the consequent technical difficulties (e.g., pyrometers, which also require a temperature/stress conversion system while strain gauges are usually unsuitable for rotating equipment). These issues have been overcome through an approach called "soft-sensing", which allows extrapolating, through suitable algorithms, detailed information from the sensors that are normally installed on the machine. This section describes an instrument called the Rotor Stress Monitor (RSM) [27], depicted in Figure 1, which is used to this purpose and whose calculation steps are described in detail in the next sections. tracking formulation of the optimization strategy, in which the maximum power reference is received from the Distributed Control System (DCS) of the plant. In particular, this article presents several novelties. In the first place, it compares the standard baseline curves for steam turbine start-up with the application of a PI architecture including the soft sensor calculating the rotor stress. Later this PI architecture is compared with the NLMPC application showing the benefits of this control procedure. It also has to be remarked that while other works aim at the application of MPC strategies to the whole plant, here the focus lies on the control of the single steam turbine. The application of these procedures to a CSP case through real data is of particular interest. The paper is organized as follows. Section 2 describes the modeling methodologies and validation procedures for the models obtained. Section 3 provides a brief introduction on PID and NMPC control and describes the structure of each implemented controller and the associated technical solutions. Section 4 discusses the results of the simulation studies for the implemented controllers by comparing their performances and providing guidelines for tuning the NMPC control parameters. Section 5 provide final comments to the obtained results and some concluding remarks.

Steam Turbine Modelling Approach and Rotor Stress Monitoring
The first fundamental step to define an effective control strategy for optimizing the ST start-up is to design a tool capable of monitoring the thermal state of the machine. As previously mentioned, the stress can only be measured using particularly invasive hardware measuring instruments, which must be designed to work at rather high temperatures and pressures, with all the consequent technical difficulties (e.g., pyrometers, which also require a temperature/stress conversion system while strain gauges are usually unsuitable for rotating equipment). These issues have been overcome through an approach called "soft-sensing", which allows extrapolating, through suitable algorithms, detailed information from the sensors that are normally installed on the machine. This section describes an instrument called the Rotor Stress Monitor (RSM) [27], depicted in Figure 1, which is used to this purpose and whose calculation steps are described in detail in the next sections. Figure 1. Rotor Stress Monitor, Matlab/Simulink implementation. Cyan represents the measured variables (pressures, temperatures, speed, generator breaker state, and shutdown command), orange represents the thermodynamic parameters, magenta represents the tuning parameters of calculation core, and red represents the mechanical stress parameters. The RSM algorithm allows calculating the current stress, its current limit, and the temperatures of the rotor based on: • Derivation of boundary conditions for the thermal model.

Thermal Model
In order to correctly evaluate thermal stresses, the model must satisfy the minimum requirements in terms of accuracy. However, it is also extremely important for the model to take into account maximum computational cost requirements set by common control industrial hardware e.g., PLCs. Therefore, standard detailed calculation procedures (such as the one used during rotor design phase) based on the 2D and 3D Finite Element Model (FEM) implementing rotor detailed geometrical description are not suitable for control purposes. Such models, which are used to explore the effect of complex thermal transient scenarios, must be simplified into Reduced Order Models (ROM).
In the present case, the selected approach is based on the assumption that the rotor geometry can be approximated through a cylindrical axisymmetric model so that the thermal model can be described (in cylindrical coordinates) by the following partial differential equation: Subjected to the following boundary conditions: where T is the rotor temperature, function of radius r and time t, HTC(t) is the HTC as a function of time, T bulk is the steam temperature, ρ is the rotor density, C p is the rotor local heat capacity (function of temperature), and K is the rotor local thermal conductivity (function of temperature). It can be observed that: • The main assumption is that the rotor geometry can be described by an infinitely extended cylinder excluding the dependency from the axial coordinate; • Local variations of material properties are taken into account leading to a non-linear heat transfer equation with a non-homogeneous material; • Since HTC and T bulk are time-dependent the model is suitable to complex heat transient scenarios.
For the solution of the problem, an explicit finite difference scheme (backward Euler) has been selected with N radial locations and a spatial discretization of ∆r = r i+1 − r i . As a consequence: Energies 2021, 14, 3998 6 of 30 As well known the stability of this numerical scheme depends on the selected time discretization ∆t. Partial derivatives in space and time are discretized as follows: The choice of an explicit method is driven by the fact that in this way, a non-linear solver is not necessary to evaluate rotor radial temperature distribution from step to step. In any case, this poses a limit on the maximum allowable timestep length that has been explored considering the thermal diffusivity of the material of interest for this kind of analysis. In addition, the limit to the timestep length poses challenges for the development of the control scheme (which are later explained). Noticeably, after numerical discretization, the nonlinear partial differential equation can be reconducted to a general control scheme: A (j) ∈ R n×n is a tri-diagonal matrix wherein the non-zero elements are expressed as follows: where α (j) = k (j) is the thermal diffusivity evaluated at the node j. ϕ is the heat flow transferred by convection to the cylinder: The previous expressions highlight a linear dependence of the cylinder temperature state on the heat flow (such as also highlighted by Nakai et al. [17]). The creation of a purely linear control system based on this model is basically compromised by the fact that heat flow is a non-linear combination of disturbance (bulk temperature), system state (last node temperature), and control (HTC). In addition, as the rotor surface temperature approaches bulk temperature, the effect of a control action on HTC becomes less relevant.

Rotor Stress Models
The temperature field can be exploited to obtain the cylinder stress approximating the actual rotor stress. For this specific case, a generalized plane strain assumption is suitable, leading to the following expression for the Von Mises stress σ eqv as a function of temperature T(r): where β is the thermal expansion coefficient, E is the Young modulus, and ν is the Poisson ratio. Equation (11) provides an estimate of the equivalent stress on the cylinder external surface and is suitable for the identification of crack nucleation (a similar expression can be obtained also for rotor core stress and can be used for detecting crack nucleation and propagation). The effect of centrifugal stress is not taken into account here, as it does not depend on the temperature field, but can be included in an additional term for control purposes.
The stress integral expression can be numerically approximated according to the following linear relation: where: Moreover, the maximum allowable stress for the rotor can be expressed as a function of the local cylinder temperature in order to take into account the typical degradation of steel strength with temperature. Usually, a polynomial approximation is suitable for this purpose. In addition, the coefficients of the polynomial can be updated from start to start in accordance with the expected start-up severity and frequency (for instance, in CSP applications the cold start-up, being the most severe and also less frequent, is usually set up with a high maximum allowable stress limit).
In any case, an actual rotor geometry holds different critical locations, unlike the considered cylinder model, where stress values are amplified. Amplification factors for stress need to be calculated through a dedicated 2D thermo-mechanical Finite Element Analysis (FEA). These factors can be used to properly tune and align the results of the simplified cylindrical model to those of a real rotor geometry (examples of this kind of procedure are given for example in [9,28]). Typically, critical rotor locations are represented by fillets in the proximity of the first stage and bucket tail grooves, where differential expansion drives complex stress conditions and even local plasticization. An example is shown in Figure 2, where the contour plots of Von Mises equivalent stress are shown for an axis-symmetric simulation, at the load step where the stress peak is reached.
where is the thermal expansion coefficient, is the Young modulus, and is the Poisson ratio. Equation (11) provides an estimate of the equivalent stress on the cylinder external surface and is suitable for the identification of crack nucleation (a similar expression can be obtained also for rotor core stress and can be used for detecting crack nucleation and propagation). The effect of centrifugal stress is not taken into account here, as it does not depend on the temperature field, but can be included in an additional term for control purposes.
The stress integral expression can be numerically approximated according to the following linear relation: where: Moreover, the maximum allowable stress for the rotor can be expressed as a function of the local cylinder temperature in order to take into account the typical degradation of steel strength with temperature. Usually, a polynomial approximation is suitable for this purpose. In addition, the coefficients of the polynomial can be updated from start to start in accordance with the expected start-up severity and frequency (for instance, in CSP applications the cold start-up, being the most severe and also less frequent, is usually set up with a high maximum allowable stress limit).
In any case, an actual rotor geometry holds different critical locations, unlike the considered cylinder model, where stress values are amplified. Amplification factors for stress need to be calculated through a dedicated 2D thermo-mechanical Finite Element Analysis (FEA). These factors can be used to properly tune and align the results of the simplified cylindrical model to those of a real rotor geometry (examples of this kind of procedure are given for example in [9,28]). Typically, critical rotor locations are represented by fillets in the proximity of the first stage and bucket tail grooves, where differential expansion drives complex stress conditions and even local plasticization. An example is shown in Figure 2, where the contour plots of Von Mises equivalent stress are shown for an axissymmetric simulation, at the load step where the stress peak is reached.  Two set of tuning functions are created to enable a good correspondence between Finite Element Models and Rotor Stress Monitor. Many shapes of increasing complexity can be chosen for the tuning functions that adjust the inputs and the outputs of the previously explained Reduced Order Model. The performance of the adjustment procedure on tuning functions can be evaluated by testing campaigns: In other words, collecting the response of the system to a stepwise and/or ramp increase of heat flow is able to give relevant indications on the behavior shown in more complex scenarios. A simple example of the response to a testing mission before and after the tuning procedure is shown in Figure 3, where the stress evaluated by FEA in the most critical location is also plotted. Here it can be seen how the chosen tuning function setup is able to match the stress response to a ramp increase of heat flow within a 5% error. At the same time, we chose to keep a conservative margin (~30%) for the response of the system to a step heat flow increase. relevant indications on the behavior shown in more complex scenarios. A simple example of the response to a testing mission before and after the tuning procedure is shown in Figure 3, where the stress evaluated by FEA in the most critical location is also plotted. Here it can be seen how the chosen tuning function setup is able to match the stress response to a ramp increase of heat flow within a 5% error. At the same time, we chose to keep a conservative margin (~30%) for the response of the system to a step heat flow increase. Figure 3. Testing of the tuning procedure to step and ramp increase of heat flow for the rotor stress monitoring, highlighting the stress trends before and after the tuning with respect to the FEM Simulation © 2020 Baker Hughes Company-All rights reserved.

Power and HTC Calculation
The thermo-mechanical FE transient simulation, used to assess ST components' reliability and strength, involves, during the pre-processing phase, a detailed calculation of HTCs and bulk temperatures. These calculations take into account the thermodynamic behavior and geometry of the specific ST and, in the most complex scenario, may also involve Computational Fluid Dynamics (CFD) evaluation with all the geometrical details of the specific ST. These kinds of calculations need to be included in a simplified form for rotor stress monitoring and control purposes through simplified thermodynamic transfer functions. Such functions basically summarize the thermodynamic behavior of the turbine and enable the HTCs calculation, bulk temperature, or ST output power for the location where the stress has to be evaluated. The need for simplified calculations lies in the computational constraint of a User Control Panel implementation. Typically, polynomial Figure 3. Testing of the tuning procedure to step and ramp increase of heat flow for the rotor stress monitoring, highlighting the stress trends before and after the tuning with respect to the FEM Simulation © 2020 Baker Hughes Company-All rights reserved.

Power and HTC Calculation
The thermo-mechanical FE transient simulation, used to assess ST components' reliability and strength, involves, during the pre-processing phase, a detailed calculation of HTCs and bulk temperatures. These calculations take into account the thermodynamic behavior and geometry of the specific ST and, in the most complex scenario, may also involve Computational Fluid Dynamics (CFD) evaluation with all the geometrical details of the specific ST. These kinds of calculations need to be included in a simplified form for rotor stress monitoring and control purposes through simplified thermodynamic transfer functions. Such functions basically summarize the thermodynamic behavior of the turbine and enable the HTCs calculation, bulk temperature, or ST output power for the location where the stress has to be evaluated. The need for simplified calculations lies in the computational constraint of a User Control Panel implementation. Typically, polynomial transfer functions are capable to model, with sufficient precision, the HTCs variation corresponding to the operation field of the rotor. A typical scenario occurs when the influence of controlled and un-controlled extraction can be considered negligible for the determination of the ST condition during start-up, and HTC, power W, and bulk temperature can be approximated by polynomial functions of ST inlet and outlet pressure, inlet temperature, and wheel chamber pressure (the steam pressure after the group including the control stage and the control valve), as follows:

Rotor Stress Control Strategies
This section focuses on the main aspects of PID and MPC strategies in order to introduce the NMPC approach aimed at optimizing the ST power startup and limiting the related rotor stress.

Background on PI and MPC Control Strategies
The PID controller is the most widespread control algorithm in industry, exploited in more than 90% industrial controllers [29], due to its simple implementation and its effectiveness. Its main characteristic lies in the use of an error-driven philosophy typical of the Model Free control approach [30]. Its implementations may differ depending on the specific application and its realizations have been widely studied in the literature and applied to industrial cases.
The success of PID control strategy is closely linked to factors related to its simplicity as well as to the effectiveness and large availability of efficient and ready-to-use PID controller products. PID is one of the first techniques developed in the field of control engineering, as it was possible to implement it through mechanisms (such as in the case of the Centrifugal governor, the first ST controller that was based on proportional control) or through simple electromechanical circuits and, afterwards, through analog and digital electronics. Its effectiveness also relies on the immediacy of the control action, a necessary feature in the case of low-level control loops (closer to the machines or basic processes such as pressure, flow, level control, etc.). Last but not least is the effectiveness of the control action on processes that work in almost stable operating points, in which the performances achievable through a PID controller are usually quite satisfactory.
The positive aspects of the PID controller are, however, counterbalanced by a series of disadvantages that make it not always the ideal choice for process control. In the first point, PIDs are a Single Input Single Output (SISO) and can act on a single manipulable variable. Multi Input Multi Output (MIMO) control can be obtained through complex control methodologies such as the decoupling of the system dynamics [31] or other more complex logics [32]. However, the most complex aspect, which is widely analyzed in the literature [33], is the tuning of its parameters. The most common approaches are based on standard Ziegler-Nichols (ZN) methods [34], relays [35], and several analytic rules [36]. Auto-tuning methods and approaches based on offline simulation are particularly studied and deepened in literature [37]. Furthermore, AI through Artificial Neural Networks (ANN), Fuzzy Logic, and Evolutionary Algorithms have been widely explored [38].
The MPC is an advanced control technique developed in the refinery and petrochemical industry context in the 1970s, and subsequently deepened in the academic world starting from an initial diffidence and skepticism that allowed, firstly, the definition and stabilization of the theoretical foundations of this methodology and, subsequently, its consistent evolution and research advancement. The foundations of this control technique lie in three basic ideas: (i) The prediction of future behavior of the controlled process, through adequate mathematical models; (ii) the optimization of the control path aimed at minimizing a specific cost function and satisfying process and actuators operative constraints; and (iii) the action on the controlled process through a receding control horizon that implements a closed loop control strategy by manipulating only the current manipulable variable, such as introduced by Propoi [39].
The ideas behind the first applications, based on Dynamic Matrix Control (DMC) [40] and the Identification-Command (IDCOM) software [41], were extended in two fundamental works of Clarke et al. [42,43] and then unified by Garcia et al. in 1989 in their seminal paper [44]. The numerous works in the literature have shown, over time, from an academic point of view, the effectiveness and limitations of this control approach, which is now considered to be the standard for the control of MIMO systems, not only in the industrial context, in which it is obviously widely used [45][46][47][48]. Some novel interesting examples of MPC methodologies application related to the energy management or power generation have been presented by Camacho et al. in their work [49], in which a hybrid  [50] that presents a nonlinear economic MPC approach to maximize the energy produced by wind turbines over all the operating points in variable speed conditions, and by Darwish et al., who presented an integration of Dual-mode MPC and Generalized Predictive Control (GPC) for optimizing the startup of a steam power plant [51].
In deeper detail, the MPC strategy requires three well-defined elements in order to achieve satisfactory control performances: • A sufficiently accurate process model, i.e., the core of the MPC approach, whose complexity must balance two counteracting issues related to the accuracy of the shortlong term dynamic predictions and the computation burden required to simulate it. In the literature, several realizations of the process models are proposed, starting from linear state space models, linear or nonlinear autoregressive exogenous models (ARX) [52], and Hammerstein Wiener models [53], which are considered the standard for industrial process modelling, up to ANNs [54,55], which are less widespread in control applications but particularly effective in modelling nonlinear behavior. Furthermore, the modelling objective it is not only related to the future behavior prediction, but also to the estimation of the current state of the controlled process, which might not be accessible or measurable. • A cost function describing one or more indexes to be optimized. Such cost function depends on the specific task and in general is related to energy description of the system, economic costs, or particular specific mathematical relationships. From this point of view, the MPC can fundamentally be a single or a multi objective optimization problem. • A set of constraints aimed at describing the operative boundaries of the system, including the process itself as well as actuators and related subsystems.
The combination of these three elements defines the complexity of the MPC and, consequently, the difficulty in solving the optimization strategy in real-time. The term "linear MPC" is used when the operating models and constraints of the process to be controlled are linear, and when the cost function is linear or quadratic. The term "nonlinear MPC" is adopted when one of the first two elements is nonlinear and/or the cost function can be described by a Taylor series of a degree higher than two.
The typical MPC closed loop architecture is depicted in Figure 4, where its different elements are divided in functional blocks: (i) The real controlled process/plant; (ii) the state observer (sometimes included within the MPC); and (iii) the MPC block. More in detail, the MPC block includes a mathematical formulation of the optimization problem, including constraints and process/plant dynamics, an objective function and an optimizer, which includes a set of optimization algorithms and routines tailored for the specific application. which includes a set of optimization algorithms and routines tailored for the specific application. The MPC solves in real-time an optimization problem that in the linear case can be described as: The MPC solves in real-time an optimization problem that in the linear case can be described as: where x ∈ R n , u ∈ R m , and y ∈ R o are, respectively, the state, the inputs, and the outputs of the controlled system; A ∈ R n·n , B ∈ R n·m , C ∈ R o·n , and D ∈ R o·m are the matrices of the system dynamics; x 0 is the initial state at the control step k = 0; and H, V, and W are suitable matrices. N p and N c are, respectively, the prediction horizon and the control horizon. Such concepts are graphically described in Figure 5, where the prediction horizon N p is the length of the time window in which the future process dynamic is simulated, and the control horizon N c is the number of degrees of freedom to solve the optimization problem). The prediction horizon is usually longer than the control horizon N p ≥ N c . As far as the cost function is concerned, in the case of a scalar objective function, the two most widely known formulations are reported in Equations (18) and (19). In particular, Equation (18) describes the cost formulation of a receding horizon for a generic optimization problem: The matrices ∈ • , ∈ • , and ∈ • are assumed to be symmetric positive definite in order to formulate a convex optimization problem.
Equation (19) describes a typical cost formulation of a receding horizon for an output tracking task problem, extensively reviewed by Rawlings [56], in which the output of the controlled system must asymptotically follow an available output reference : As far as the cost function is concerned, in the case of a scalar objective function, the two most widely known formulations are reported in Equations (18) and (19). In particular, Equation (18) describes the cost formulation of a receding horizon for a generic optimization problem: (18) The matrices P ∈ R n·n , Q ∈ R n·n , and S ∈ R m·m are assumed to be symmetric positive definite in order to formulate a convex optimization problem. Equation (19) describes a typical cost formulation of a receding horizon for an output tracking task problem, extensively reviewed by Rawlings [56], in which the output of the controlled system must asymptotically follow an available output reference y r : where e r is the output tracking error, defined as: Furthermore, in this case, the matrices W r ∈ R o·o , W u ∈ R l·l , and W ∆u ∈ R l·l are assumed to be symmetric positive definite in order to formulate a convex optimization problem. Once the cost function is defined, the problem is completed with a set of constraints on the inputs, outputs, the velocity of the output ∆y(k), and other specific constraints The formulated optimization problem is sequentially solved in each control instant for a receding horizon by means of appropriate optimization algorithms. The MPC algorithm steps can be summarized as follows: 1.
At the control timestep t, MPC reads the measurements and the estimate of the current state x 0 = x(t) of the process. The initial state defines the simulation starting point of the controlled process; 2.
The optimization algorithm solves, in real time, the formulated problem for a prediction horizon of N p time steps ahead, generating a control sequence u, which minimize the objective function J(x t , N p , N c ), while respecting the constraints: 3. The MPC acts on the controlled system, by implementing the first calculated control move u(1). The MPC keeps in memory the remaining control sequence (k = 2, · · · N c ) as a warm start of the optimization of the next control step. 4.
In the next control step, the MPC strategy comes back to step 1, by implementing a closed loop control, since the applied action on the controlled process is optimized starting from the measurement or the estimate of its current state.
The MPC formulation can be naturally extended to the more complex nonlinear case, which can be generically summarized as follows: where J(x(k), u(k), k) is a generic nonlinear objective function, x is the process state, u is the vector of the manipulated variables, E x t + N p is the terminal cost, N p expresses the length of the prediction horizon, f and h are nonlinear functions describing the dynamic of the process state and output, g, l, r are nonlinear functions describing the constraints on the system, and r x t + N p ≤ 0 is named terminal constraint.
The nonlinear nature of the optimization problem leads to several issues from the computational burden point of view, which particularly differentiates the NMPC from the linear MPC case. Firstly, nonlinear optimization problems are not always characterized by convexity, and this fact does not ensure either the achievement of a global minimum, or in some cases, the convergence of optimization algorithms to a feasible solution. In the case of "simple" nonlinear problem, the main issue is calculating a solution in real-time, while ensuring a certain grade of optimality. For this reason, NMPC is less widespread and applied with respect to linear MPC. Some applications are found for slow processes, whose dynamics are compatible with calculation times. A good formulation often helps in rapidly solving the optimization problem, otherwise it is still possible to implement sub-optimal strategies [57] where the problem is feasible. A possible trade-off solution is offered by the so-called Real-Time Iteration approach [58], which implements a formulation of the optimal problem and of the optimization routines that are a "middle ground" between the linear and the nonlinear case. In cases where the real-time solution of the optimum problem is in any case not feasible, a further approach is offered by the use of methodologies called Explicit MPC, in which the optimum problem is solved offline on a discretization of the operating space

PI Based Rotor Stress Controller
A preliminary implementation using a rotor stress monitor as a basis for turbine control can be obtained through a PI control strategy. The main idea behind this control approach is using an outer control ring with a first PI to create a stress set-point from turbine power current value and current set point. This set point is limited in accordance with the maximum allowed stress. The stress set-point and the stress current value are used by an inner PI in order to determine the wheel chamber pressure set-point that is compliant with the desired stress limit. This kind of architecture is able to grant good results in terms of both production increment during start-up and stress limit compliance. An example is shown in Figure 6 where a fully cold start-up for a power gen application based on a PI rotor stress controller is compared with a standard start-up curve. In this case, the full steam header conditions allow highlighting the great advantage that can be achieved by exploiting dedicated rotor stress controller with respect to the standard startup procedures based on pre-calculated power references. In the specific case study, it is possible to increase the energy produced during the startup by over 15-20%.
At the same time, stress overshoot might appear as PI controllers are subject to finetuning errors. Since the effect of possible mistuning introduces the need for safety factors on the maximum allowed stress, the need of different control routines is justified (as explained in the next sections).

Nonlinear Model Predictive Rotor Stress Controller
The NMPC approach, shown in Figure 7, has been developed in the Matlab/Simulink environment and consists of three main software components:

•
The rotor stress monitor, i.e., the soft sensor presented in Section 2, which estimates the thermal status, focusing on the temperature of the rotor discretized in N points.
The soft sensor is used as a state observer within the NMPC approach, externally installed in PLC platform as an integral part of the steam turbine governor. The block calculates the state of the machine, starting from some relevant inputs: The internal and external case temperature, the pressure and the temperature in the steam header, and the wheel-chamber and outlet pressures.

•
The NMPC algorithm, which calculates a reference control action for the pressure in the wheel chamber.
• A PI controller, which directly acts on the inlet valve by following the reference from the NMPC algorithm. This particular control architecture aims at overcoming the issues related to the nonlinear behavior of the valve system, without including it in the modelling strategy formulated within the NMPC software.
power current value and current set point. This set point is limited in accordance with the maximum allowed stress. The stress set-point and the stress current value are used by an inner PI in order to determine the wheel chamber pressure set-point that is compliant with the desired stress limit. This kind of architecture is able to grant good results in terms of both production increment during start-up and stress limit compliance. An example is shown in Figure 6 where a fully cold start-up for a power gen application based on a PI rotor stress controller is compared with a standard start-up curve. In this case, the full steam header conditions allow highlighting the great advantage that can be achieved by exploiting dedicated rotor stress controller with respect to the standard startup procedures based on pre-calculated power references. In the specific case study, it is possible to increase the energy produced during the startup by over 15-20%. At the same time, stress overshoot might appear as PI controllers are subject to finetuning errors. Since the effect of possible mistuning introduces the need for safety factors on the maximum allowed stress, the need of different control routines is justified (as explained in the next sections).

Nonlinear Model Predictive Rotor Stress Controller
The NMPC approach, shown in Figure 7, has been developed in the Matlab/Simulink environment and consists of three main software components:

•
The rotor stress monitor, i.e., the soft sensor presented in Section 2, which estimates the thermal status, focusing on the temperature of the rotor discretized in N points. The soft sensor is used as a state observer within the NMPC approach, externally installed in PLC platform as an integral part of the steam turbine governor. The block calculates the state of the machine, starting from some relevant inputs: The internal and external case temperature, the pressure and the temperature in the steam header, and the wheel-chamber and outlet pressures.

•
The NMPC algorithm, which calculates a reference control action for the pressure in the wheel chamber. • A PI controller, which directly acts on the inlet valve by following the reference from the NMPC algorithm. This particular control architecture aims at overcoming the issues related to the nonlinear behavior of the valve system, without including it in the modelling strategy formulated within the NMPC software. The tailored NMPC application consists of 4 main blocks depicted in detail in Figure  8: • A Startup Selector for selecting the parameters related to the maximum allowable stress according to the casing temperature and type of startup (cold, warm, hot).

•
A Disturbance Predictor and a Max Pressure Predictor, which estimate, respectively, the future steam temperatures and the pressure available at the header. The steam temperature acts on the thermal dynamics as a non-controllable variable (disturbance). The predicted future header pressure is necessary to estimate its maximum constraints within the problem formulated in the NMPC. In deeper detail, the future predictions are calculated starting from the related current through a linear extrapolation. The linear extrapolation is saturated between a minimum and a maximum allowable value.

•
The Optimizer, the main block of the NMPC algorithm, which implements the optimization problem, routines, and specific logics for calculating the control action. The The tailored NMPC application consists of 4 main blocks depicted in detail in Figure 8: • A Startup Selector for selecting the parameters related to the maximum allowable stress according to the casing temperature and type of startup (cold, warm, hot).

•
A Disturbance Predictor and a Max Pressure Predictor, which estimate, respectively, the future steam temperatures and the pressure available at the header. The steam temperature acts on the thermal dynamics as a non-controllable variable (disturbance). The predicted future header pressure is necessary to estimate its maximum constraints within the problem formulated in the NMPC. In deeper detail, the future predictions are calculated starting from the related current through a linear extrapolation. The linear extrapolation is saturated between a minimum and a maximum allowable value.

•
The Optimizer, the main block of the NMPC algorithm, which implements the optimization problem, routines, and specific logics for calculating the control action. The optimization problem, formulated for a receding horizon problem, is described by the following equations: where N pred is the prediction horizon in samples, u is the manipulated variable, the pressure in the wheel chamber, P re f is the turbine mechanical power reference received from the DCS of the system that takes into account the current limits of the steam generation system, W P and W ∆u are, respectively, the weight for the reference tracking and the variation of the control action, f 1 and f 2 are nonlinear functions, which describe the relationships between the main thermodynamic variables of the ST and HTC and the generated mechanical power P, respectively, and σ min and σ max are respectively the minimum and maximum limits of the stress, designed in function of the specific mission of the turbine. where is the prediction horizon in samples, is the manipulated variable, the pressure in the wheel chamber, is the turbine mechanical power reference received from the DCS of the system that takes into account the current limits of the steam generation system, and Δ are, respectively, the weight for the reference tracking and the variation of the control action, 1 and 2 are nonlinear functions, which describe the relationships between the main thermodynamic variables of the ST and HTC and the generated mechanical power , respectively, and and are respectively the minimum and maximum limits of the stress, designed in function of the specific mission of the turbine. The control sequence = ⋯ is calculated with a blocking strategy, in which the calculated variables are spatially distributed along the prediction horizon (such as shown in Figure 9) and linearly interpolated in the intermediate points with a sampling time , which also discretizes the dynamics of the thermic rotor model. With respect to standard blocking strategies, which concentrate the control action in the beginning of the control horizon, the presented approach aims at limiting the peak of stress for The control sequence u = u 1 u 2 · · · u N X is calculated with a blocking strategy, in which the N X calculated variables are spatially distributed along the prediction horizon (such as shown in Figure 9) and linearly interpolated in the intermediate points with a sampling time T s , which also discretizes the dynamics of the thermic rotor model. With respect to standard blocking strategies, which concentrate the control action in the beginning of the control horizon, the presented approach aims at limiting the peak of stress for the whole prediction horizon, in order to overcome issues related to the infeasibility of the optimization problem due to a short control horizon. The presented formulation is characterized by a large number of constraints mainly due to the evaluation of machine dynamics every for the aforementioned stability issues. In order to reduce the computational burden, two strategies have been implemented within the NMPC algorithm. The first one concerns the decimation of the number of power and stress constraints evaluation along the prediction horizon. The second one concerns the scheduling of the optimizer with a frequency lower than the frequency relative to the model simulation sampling time.
In detail, the decimation of evaluation of stress constraints is motivated by the fact that the stress dynamics follow the thermal ones and are therefore slow. For this reason, it is possible to evaluate the constraints every seconds, where is the decimation parameter: The power constraints are evaluated in the same way, as their dynamics are bound by the limitation on the wheel chamber pressure, which is in turn limited by stress-derived limits. A well-tuned decimation parameter , allows a considerable reduction of the number of constraints on stress and computational costs.
The second strategy implements the decimation of the frequency of the optimizer respect to the NMPC frequency : The optimizer is scheduled with a frequency multiple of the NMPC scheduling frequency: A well-tuned parameter allows overcoming the limitation concerning the fact that the optimizer computational cost is not compatible with the controller frequency. The presented formulation is characterized by a large number of constraints mainly due to the evaluation of machine dynamics every T s for the aforementioned stability issues. In order to reduce the computational burden, two strategies have been implemented within the NMPC algorithm. The first one concerns the decimation of the number of power and stress constraints evaluation along the prediction horizon. The second one concerns the scheduling of the optimizer with a frequency lower than the frequency relative to the model simulation sampling time.
In detail, the decimation of evaluation of stress constraints is motivated by the fact that the stress dynamics follow the thermal ones and are therefore slow. For this reason, it is possible to evaluate the constraints every k cd T s seconds, where k cd is the decimation parameter: The power constraints are evaluated in the same way, as their dynamics are bound by the limitation on the wheel chamber pressure, which is in turn limited by stress-derived limits. A well-tuned decimation parameter k cd , allows a considerable reduction of the number of constraints on stress and computational costs.
The second strategy implements the decimation of the frequency of the optimizer respect to the NMPC frequency f c : Energies The optimizer is scheduled with a frequency f opt multiple of the NMPC scheduling frequency: A well-tuned k mpc-period parameter allows overcoming the limitation concerning the fact that the optimizer computational cost is not compatible with the controller frequency.
The optimizer routine is based on the SQP algorithm [61]. In particular, in this work, two different implementations have been used: The Matlab fmincon routine used as a benchmark, and a tailored Sequential Quadratic Programming (SQP) approach, which is license-free as it was ad-hoc developed.

Tailored SQP Algorithm
SQP is a state-of-the-art algorithm that allows solving general nonlinear optimization problems in the form: where f (x), b(x), c(x) are generic non-linear functions. Nonlinear optimization is an NP-hard problem, namely it is computationally expensive to solve. The basic idea of SQP is to form quadratic sub-problems that are computationally easier and iteratively solve them in order to converge to a solution of the global problem. The approach followed is similar to the one proposed in [62]. Let: be the Lagrangian of the problem. Let H(x) be the Hessian of L(x, σ, λ) as a function only of x at some fixed values of (σ, λ). Let ∆ f be the gradient of the function f and J b , J c the Jacobians of the functions b, c. The sub-problem to solve at iteration k is formulated as follows: The expression is a quadratic problem in the variable d, positive definite since H(x k ) is kept positive definite. The solution d k to the quadratic sub-problem provides the direction of improvement at the next iteration of the original problem solution: To improve the convergence, the step length α is calculated in order to update the solution: where α is computed by optimizing a merit function. If stopping conditions or constraints violation are satisfied, the current solution is given as output. Otherwise, a new iteration is performed. This algorithm converges to a minimum of the objective functions [62] and performs fast iterations, since the quadratic sub-problem is convex at every iteration: Convex quadratic problems can be solved in polynomial time. For this purpose, an active set method has been exploited; the problem is iteratively solved starting from a feasible initial solution calculated with the simplex method. At every iteration, a sub-problem taken only on some equality (active) constraint is solved through QR matrix decompositions [63]. The solution is then iteratively updated in a way that preserves feasibility. This procedure converges in polynomial time to the optimum of the quadratic sub-problem.
With respect to the fmincon routine, the tailored algorithm manages the tolerance on constraints, solutions, and jacobian in a different way. The algorithm tries to converge starting from a minimum allowed tolerance and increasing it of a factor if it is not respected at the end of the number of allowed iterations.
In order to make this algorithm implementable on a PLC platform in the future, the algorithm has been developed in the Matlab environment and the opportunity of translating it into the PLC native language through the PLC Coder Toolbox has been tested. As many Matlab routines (such as the fmincon function, the linear(quadratic) solvers, etc.) cannot be directly translated by the PLC Coder, in their place, the above tailored algorithms have been coded and used. The SQP implementation has been tested with respect to Matlab fmincon for several problems and showed similar or sometimes even better performances. In particular, the two algorithms behave similarly on instances of the nonlinear problem arising from the NMPC. With respect to the fmincon routine, the tailored SQP algorithm for PLC platforms is currently 8-10 times slower.

Results
Extensive simulation studies have been performed, through which the performance of the proposed NMPC approach has been compared to the selected benchmark based on the PI Controller as well as to a second benchmark based on the standard procedures for power startup, based on pre-calculated load curves. In order to assess the validity of the NMPC approach, all the control approaches have been simulated by using real data related to a Low-Pressure (LP) ST with a rated power of 45 MW, an inlet steam pressure of 1-22 Bar, an inlet steam temperature of 180-380 • C, and a rated speed 3000 RPM. The external radius is assumed to be equal to 29.5 cm in the equivalent cylindrical model and the rotor is assumed to be made of ferritic-martensitic heat resisting forged steel for power generation applications. Even if material properties cannot be disclosed, similar results can be obtained assuming material properties of A470 class 8 whose data are easily available, for example in the NIMS database. It has to be remarked that these material properties are not fundamental to quantify the goodness of the control strategy (a start-up performed with a worse-performing rotor material determines a longer start-up but does not invalidate the control strategy performance). The machine is discontinuously operated in a CSPP, where the steam characteristics vary as a function of the weather with numerous starts per year. The specific case study is related to a cold startup of a typical winter day, during which the steam does not reach the full header conditions in terms of mass flow and enthalpy with respect to the maximum potential of the plant. The steam characteristics are depicted in Figure 10: The top diagram shows the pressure on the steam header (in blue), the temperature of the inlet steam (red line), and the internal and external temperatures of the LP machine casing (green and orange lines, respectively); the central diagram depicts the rotor speed, while the bottom diagram shows the startup command (blue line) and the state of the generator breaker (orange line). For confidentiality constraints, the time axes are dimensionless.
The control performance has been evaluated through three main indexes, related to the maximum percentage stress overshoot, the generated mechanical energy during the startup, and the computing time. Firstly, the performance of PI controller has been evaluated with respect to the standard startup procedures, starting from a brief discussion related to the controller parameter selection. In a second phase of the work, extensive tests have been performed on the NMPC approach, in order to evaluate the performance sensitivity related to the main controller parameter.

PI Control Results
In order to assess the performances of the PI based controller, it is important to verify the sensitivity of this kind of controller to a possible mistuning of its gain parameters. It has to be remarked that during commissioning the gain parameters of power PI are usually modified to obtain the desired control performances. In this scenario, the control architecture should have pre-set values enabling a faster commissioning. A sensitivity analysis based on a Design of Experiments approach solves this kind of problem with the main objective of identifying the optimal gain values of the internal stress PI. In more detail, the controller has been initially tuned through a proprietary site procedure similar to ZN methodologies, adapted to the constraints of the site itself. In addition, a sensitivity analysis was carried out on the performance, in terms of generated energy and stress overshoot, as the controller parameters varied, in a large neighborhood of the operative point identified by the site procedures, in order to verify the robustness of the methodologies applied to the specific case. The results of the sensitivity analysis are shown in Figure 11 where the energy production and the maximum percentage overshoot are shown. The control performance has been evaluated through three main indexes, related to the maximum percentage stress overshoot, the generated mechanical energy during the startup, and the computing time. Firstly, the performance of PI controller has been evaluated with respect to the standard startup procedures, starting from a brief discussion related to the controller parameter selection. In a second phase of the work, extensive tests have been performed on the NMPC approach, in order to evaluate the performance sensitivity related to the main controller parameter.

PI Control Results
In order to assess the performances of the PI based controller, it is important to verify the sensitivity of this kind of controller to a possible mistuning of its gain parameters. It has to be remarked that during commissioning the gain parameters of power PI are usually modified to obtain the desired control performances. In this scenario, the control architecture should have pre-set values enabling a faster commissioning. A sensitivity analysis based on a Design of Experiments approach solves this kind of problem with the main objective of identifying the optimal gain values of the internal stress PI. In more detail, the controller has been initially tuned through a proprietary site procedure similar to ZN methodologies, adapted to the constraints of the site itself. In addition, a sensitivity analysis was carried out on the performance, in terms of generated energy and stress overshoot, as the controller parameters varied, in a large neighborhood of the operative point identified by the site procedures, in order to verify the robustness of the methodologies applied to the specific case. The results of the sensitivity analysis are shown in Figure 11  The results of the DoE show that the peak production also corresponds to the maximum stress overshoot. The results show that the overshoot peak cannot be fully eliminated unless heavily scarifying the performances through a very large integral coefficient of the PI Controller. The solution requires introducing an additional safety factor on the maximum allowable stress for the turbine rotor, in order to limit the stress reference on The results of the DoE show that the peak production also corresponds to the maximum stress overshoot. The results show that the overshoot peak cannot be fully eliminated unless heavily scarifying the performances through a very large integral coefficient of the PI Controller. The solution requires introducing an additional safety factor on the maximum allowable stress for the turbine rotor, in order to limit the stress reference on the PI controller, with a consequent limited reduction of the energy production during start-up. An example of this behavior is shown in Figures 12 and 13 where the PI intervention can be clearly identified between 48-65% of the simulation time t sim .
In particular, in Figure 12, the top diagram shows the equivalent rotor stress in the core (red), the rotor skin stress (blue), the maximum between core and skin stress (green), and the maximum allowable rotor stress in violet. The bottom figure shows the steam temperature in black, the temperature of the external and internal casing (dashed lines in orange and yellow), and the temperature of the rotor in different point of the radius (from blue in the core to red in the skin). In Figure 13, the top diagram shows the controlled steam pressure in the wheel chamber (blue) and the available pressure in the header (orange). The bottom figure shows the rotor speed in blue and the mechanical power in the turbine shaft normalized with respect to the rated power.
(a) (b) Figure 11. Results of sensitivity analysis for the PI controller parameters: (a) Energy production and (b) maximum percentage overshoot.
The results of the DoE show that the peak production also corresponds to the maximum stress overshoot. The results show that the overshoot peak cannot be fully eliminated unless heavily scarifying the performances through a very large integral coefficient of the PI Controller. The solution requires introducing an additional safety factor on the maximum allowable stress for the turbine rotor, in order to limit the stress reference on the PI controller, with a consequent limited reduction of the energy production during start-up. An example of this behavior is shown in Figures 12 and 13 where the PI intervention can be clearly identified between 48-65% of the simulation time . In particular, in Figure 12, the top diagram shows the equivalent rotor stress in the core (red), the rotor skin stress (blue), the maximum between core and skin stress (green), and the maximum allowable rotor stress in violet. The bottom figure shows the steam temperature in black, the temperature of the external and internal casing (dashed lines in orange and yellow), and the temperature of the rotor in different point of the radius (from blue in the core to red in the skin). In Figure 13, the top diagram shows the controlled steam pressure in the wheel chamber (blue) and the available pressure in the header (orange). The bottom figure shows the rotor speed in blue and the mechanical power in the turbine shaft normalized with respect to the rated power.

DOEs NMPC
This section presents the outcomes of the main tests on the NMPC, which were aimed at evaluating the performances of this approach related to both energy production and startup speed up and viability in terms of computing time. In particular, the extensive test campaign mainly focuses on balancing the control performances and the heavy computation burden of NMPC, with the objective of finding an optimal configuration of its parameters. Among the different tests that have been performed, we present here the main results related to the sensitivity analysis with respect to control, optimizer, and formulation parameters (namely, cost function parameters). In particular, the following parameters are considered: • Time-related parameters , , , in which, when 1 (the optimizer is not scheduled every control step), the control action is kept constant through a sample and hold strategy; • Main optimizer parameters , MaxIter and TolCon; Figure 13. PI controller test, control action, generated power, and turbine speed. Header and wheel chamber pressures (a); rotor speed and generated power (b).

DOEs NMPC
This section presents the outcomes of the main tests on the NMPC, which were aimed at evaluating the performances of this approach related to both energy production and startup speed up and viability in terms of computing time. In particular, the extensive test campaign mainly focuses on balancing the control performances and the heavy computation burden of NMPC, with the objective of finding an optimal configuration of its parameters. Among the different tests that have been performed, we present here the main results related to the sensitivity analysis with respect to control, optimizer, and formulation parameters (namely, cost function parameters). In particular, the following parameters are considered: • Time-related parameters N x , N pred , k mpc-period , in which, when k mpc-period > 1 (the optimizer is not scheduled every control step), the control action is kept constant through a sample and hold strategy; • Main optimizer parameters N x , MaxIter and TolCon; • Cost function weights, W p and W ∆u , where W and W ∆u have been varied according to the relationships: W ∆u = αW, W p = W; The abovementioned tests are summarized through three different figures. In the figures, the energy is normalized with respect to the maximum energy produced resulted in the specific DoE related to the NMPC application to the case study. Figures 14 and 15 are related to the first test campaign. In particular, Figure 14a shows the maximum calculation time as a function of the prediction horizon and the number of calculated variables N x , which highlights an evident link between the hyperparameters involved and the calculation time. The calculation time appears to be very sensitive to the prediction horizon and slightly less sensitive to the number of variables calculated. Figure 14b shows the produced energy as a function of the varied parameters, in particular highlighting a great sensitiveness to the frequency of the optimizer and an increasing sensitivity with respect to the prediction horizon and the number of calculated freedom degrees when the optimizer is more frequently scheduled (k mpc-period ∈ [1,32]). Figure 15 shows the percentage stress overshoot as a function of the varied parameters. In particular, the graph shows that the overshoot depends only on the optimizer scheduling frequency, due to the fact that a sample and hold strategy for a less frequent optimizer activity is not efficient.
The abovementioned tests are summarized through three different figures. In the figures, the energy is normalized with respect to the maximum energy produced resulted in the specific DoE related to the NMPC application to the case study. Figures 14 and 15 are related to the first test campaign. In particular, Figure 14a shows the maximum calculation time as a function of the prediction horizon and the number of calculated variables , which highlights an evident link between the hyperparameters involved and the calculation time. The calculation time appears to be very sensitive to the prediction horizon and slightly less sensitive to the number of variables calculated. Figure 14b shows the produced energy as a function of the varied parameters, in particular highlighting a great sensitiveness to the frequency of the optimizer and an increasing sensitivity with respect to the prediction horizon and the number of calculated freedom degrees when the optimizer is more frequently scheduled ( ∈ 1, 32 ). Figure 15 shows the percentage stress overshoot as a function of the varied parameters. In particular, the graph shows that the overshoot depends only on the optimizer scheduling frequency, due to the fact that a sample and hold strategy for a less frequent optimizer activity is not efficient.   The abovementioned tests are summarized through three different figures. In the figures, the energy is normalized with respect to the maximum energy produced resulted in the specific DoE related to the NMPC application to the case study. Figures 14 and 15 are related to the first test campaign. In particular, Figure 14a shows the maximum calculation time as a function of the prediction horizon and the number of calculated variables , which highlights an evident link between the hyperparameters involved and the calculation time. The calculation time appears to be very sensitive to the prediction horizon and slightly less sensitive to the number of variables calculated. Figure 14b shows the produced energy as a function of the varied parameters, in particular highlighting a great sensitiveness to the frequency of the optimizer and an increasing sensitivity with respect to the prediction horizon and the number of calculated freedom degrees when the optimizer is more frequently scheduled ( ∈ 1, 32 ). Figure 15 shows the percentage stress overshoot as a function of the varied parameters. In particular, the graph shows that the overshoot depends only on the optimizer scheduling frequency, due to the fact that a sample and hold strategy for a less frequent optimizer activity is not efficient.     Figure 16a,b shows the mean calculation time as a function of the main optimizer parameters, highlighting that it is substantially dependent on N x and less dependent on the number of iterations. As the tolerance on the constraints increases, the calculation time decreases in terms of both average and maximum time, as the optimizer relaxes the constraints more. Figure 17a shows the produced energy in function of the varied parameters, which is dependent on the tolerance on the constraints TolCon. With a number of maximum iterations MaxIter equal to 1, the surface shows a dependence also on N x . Figure 16b shows the percentual stress overshoot in function of the varied parameters, and in particular, that overshoot mainly depends on the tolerance constraints and slightly on N x . constraints increases, the calculation time decreases in terms of both average and maximum time, as the optimizer relaxes the constraints more. Figure 17a shows the produced energy in function of the varied parameters, which is dependent on the tolerance on the constraints TolCon. With a number of maximum iterations MaxIter equal to 1, the surface shows a dependence also on . Figure 16b shows the percentual stress overshoot in function of the varied parameters, and in particular, that overshoot mainly depends on the tolerance constraints and slightly on .
(a) (b)  Figure 18 shows the results related to the third test campaign, which investigates the effect of the cost function weighs on the control performance. Figure 18a depicts the mean computing time in function of the varied parameters, while Figure 18b shows the produced energy. The test shows that both computing time and produced energy highly depend on (and, therefore, on ). constraints increases, the calculation time decreases in terms of both average and maximum time, as the optimizer relaxes the constraints more. Figure 17a shows the produced energy in function of the varied parameters, which is dependent on the tolerance on the constraints TolCon. With a number of maximum iterations MaxIter equal to 1, the surface shows a dependence also on . Figure 16b shows the percentual stress overshoot in function of the varied parameters, and in particular, that overshoot mainly depends on the tolerance constraints and slightly on .
(a) (b)  Figure 18 shows the results related to the third test campaign, which investigates the effect of the cost function weighs on the control performance. Figure 18a depicts the mean computing time in function of the varied parameters, while Figure 18b shows the produced energy. The test shows that both computing time and produced energy highly depend on (and, therefore, on ).  Figure 18 shows the results related to the third test campaign, which investigates the effect of the cost function weighs on the control performance. Figure 18a depicts the mean computing time in function of the varied parameters, while Figure 18b shows the produced energy. The test shows that both computing time and produced energy highly depend on α (and, therefore, on W ∆u ).

Guidelines for Controller Parameters Tuning and Discussion of Results
The analysis clarifies the main dependence relationships between the NMPC parameters and the resulting performance. As mentioned before, three main aspects have to be considered: (i) The performance in terms of produced energy; (ii) the overshoot of stress;

Guidelines for Controller Parameters Tuning and Discussion of Results
The analysis clarifies the main dependence relationships between the NMPC parameters and the resulting performance. As mentioned before, three main aspects have to be considered: (i) The performance in terms of produced energy; (ii) the overshoot of stress; and (iii) the computation burden.
In first place, the energy produced depends on a decreasing proportion of the frequency of scheduling of the optimizer, the maximum number of iterations MaxIter, the number of control freedom degrees N x , and the length of the prediction horizon N pred . In more detail, the energy is substantially dependent on the frequency of scheduling of the optimizer: The less frequently the optimizer is scheduled, the less the controller is capable of optimizing the starting curve of the machine. The produced energy is also partially dependent on the combination of the maximum number of iterations and the number of Degrees of Freedom (DOFs) of the control. For a low number of maximum iterations (i.e., MaxIter = 1), the produced energy highly depends on the number of control DOFs. The greater N x , the greater the amount of produced energy, since fewer iterations are, however, balanced by a greater number of variables, which allow better optimization of the control action. For a maximum number of iterations greater than 3, it is possible to obtain the maximum possible energy during start-up (if the tuning of the other parameters is optimal). Finally, the prediction horizon N pred slightly affects the produced energy, due to the accuracy of disturbance predictions, which is less accurate for a long prediction horizon. The computational cost mainly depends on the length of the time horizon N pred , on the number of future variables that are controlled N x , and, finally, on the parameter W ∆u , which weights the variation of control action between consecutive control steps. The tolerance on constraints TolCon has a role on the sensitivity analysis but only when this is rather high, resulting in a decreased computational cost. In deeper details, high values of N pred , N x , and W ∆u lead to higher computational costs. In conclusion, these four parameters must be carefully tuned for first. The length of the prediction horizon must cover the rise time of the stress in order to predict and identify the peak of stress, and, at the same time, it is necessary to balance the counteracting calculation time. Ultimately, this important parameter depends on the construction material of the machine and its size. Once the abovementioned parameters are tuned, the required computation burden-in terms of maximum computation time-can be balanced by tuning the frequency of the scheduling of the optimizer. The maximum percentage stress overshoot mainly depends on the scheduling frequency of the optimizer and on the tolerance of the optimization algorithm on constraints. The overshoot is almost linearly dependent on the two parameters.
An example of fine-tuned NMPC controller application is shown in Figures 19-21. More in detail, the simulation is performed with a sampling time T s = 2 s, a prediction horizon of 40 min (N pred = 1200 samples) and a number of control steps equal to N x = 5. The optimizer within the NMPC has been scheduled every k mpc-period = 4 periods, equal to 8 s, to guarantee an adequate margin for the computation time. Figure 19 shows on the top diagram the stress during start-up (blue in the core, red in the skin of the rotor, green the maximum stress) and the stress limit (in purple). The middle diagram depicts the temperature of the steam in the header (black), the temperature of the external and internal casing (dashed lines in orange and yellow), and the rotor temperature in each radial node (from blue in the core to red in the skin). The bottom diagram shows the state of the generator breaker (blue) and the start-up state (black). Figure 20 shows, in the top diagram, the pressure on the header (yellow) and the regulated pressure in the wheel chamber (blue). The middle diagram shows the rotor speed (blue) and the generated mechanical power (orange). Figure 21 depicts the computing time of the scheduled optimizer in the NMPC. perature of the steam in the header (black), the temperature of the external and internal casing (dashed lines in orange and yellow), and the rotor temperature in each radial node (from blue in the core to red in the skin). The bottom diagram shows the state of the generator breaker (blue) and the start-up state (black). Figure 20 shows, in the top diagram, the pressure on the header (yellow) and the regulated pressure in the wheel chamber (blue). The middle diagram shows the rotor speed (blue) and the generated mechanical power (orange). Figure 21 depicts the computing time of the scheduled optimizer in the NMPC. Figure 19. Example of application of NMPC to the startup of a LP turbine, stress conditions, steam and rotor temperatures, and startup state. Stress trends (a); temperature of the steam, wheel chamber and rotor (b); startup and generator breaker commands (c). the maximum stress) and the stress limit (in purple). The middle diagram depicts the temperature of the steam in the header (black), the temperature of the external and internal casing (dashed lines in orange and yellow), and the rotor temperature in each radial node (from blue in the core to red in the skin). The bottom diagram shows the state of the generator breaker (blue) and the start-up state (black). Figure 20 shows, in the top diagram, the pressure on the header (yellow) and the regulated pressure in the wheel chamber (blue). The middle diagram shows the rotor speed (blue) and the generated mechanical power (orange). Figure 21 depicts the computing time of the scheduled optimizer in the NMPC.   The simulation shows several interesting results. The control mission shows two singular behaviors: The first one occurs around 11-17% of simulation time, when a peak in the rotor stress is observed, which is due to the speed startup and is not controllable by the NMPC that acts on the power startup. The second one occurs around 47-53% of simulation time, when the NMPC effectively regulates the wheel chamber pressure in order to limit the thermal stress conditions in the rotor without overshoots on the rotor stress behavior.
The comparison between the NMPC and PI control strategies shows several interesting results. Table 1 summarizes the main numerical results for the three case studies: (i) PI Rotor Stress Controller; (ii) PI Rotor Stress Controller with anti-overshoot; and (iii) welltuned NMPC. All the case studies refer to the control of the steam in CSPP application, in which the steam header condition is dependent on the weather condition, and rarely reach the full header capacity in terms of temperature and pressures. The results are reported in terms of two different indexes: The Normalized Percentual Integral Absolute Error (NPIAE) and the Normalized Production Capability (NPC) that describes the ratio between the produced energy in the case of controlled rotor stress scenario and the maximum theorical energy production. The NPIAE is a normalized version of the IAE index to overcome confidentiality data issues related to the intellectual property.
where ( ) is the power produced at the time for the presented -th control scenario, max is the power produced in a scenario with uncontrolled rotor stress, and , respectively, are the initial power startup time and the instant in which the startup ends, and and max are, respectively, the energy produced in the -th controlled scenario and the uncontrolled scenario. In particular, the NMPC allow one to totally avoid the stress overshoot, which in some control applications, can result in poor or dangerous behaviors, while through the The simulation shows several interesting results. The control mission shows two singular behaviors: The first one occurs around 11-17% of simulation time, when a peak in the rotor stress is observed, which is due to the speed startup and is not controllable by the NMPC that acts on the power startup. The second one occurs around 47-53% of simulation time, when the NMPC effectively regulates the wheel chamber pressure in order to limit the thermal stress conditions in the rotor without overshoots on the rotor stress behavior.
The comparison between the NMPC and PI control strategies shows several interesting results. Table 1 summarizes the main numerical results for the three case studies: (i) PI Rotor Stress Controller; (ii) PI Rotor Stress Controller with anti-overshoot; and (iii) welltuned NMPC. All the case studies refer to the control of the steam in CSPP application, in which the steam header condition is dependent on the weather condition, and rarely reach the full header capacity in terms of temperature and pressures. The results are reported in terms of two different indexes: The Normalized Percentual Integral Absolute Error (NPIAE) and the Normalized Production Capability (NPC) that describes the ratio between the produced energy in the case of controlled rotor stress scenario and the maximum theorical energy production. The NPIAE is a normalized version of the IAE index to overcome confidentiality data issues related to the intellectual property.
where P s i (t) is the power produced at the time t for the presented i-th control scenario, P max is the power produced in a scenario with uncontrolled rotor stress, t o and t end , respectively, are the initial power startup time and the instant in which the startup ends, and E s i and E max are, respectively, the energy produced in the i-th controlled scenario and the uncontrolled scenario. In particular, the NMPC allow one to totally avoid the stress overshoot, which in some control applications, can result in poor or dangerous behaviors, while through the PI controller it is not possible to eliminate the overshoot unless one adequately limits the stress reference with safety factors (as in the second scenario). Safety factors in fact limit the energy production during the startup. Another fundamental aspect concerns the production of energy during the startup. The PI controller application shows encouraging results, considering the simplicity of the control implementation. In conditions of perfectly tuned NMPC and PI controllers, in the presented case study, the NMPC allows one to produce about 2.3% more energy than the PI-based rotor stress controller. In turn, the controller allows one to increase energy production by over 15-20% compared to what can normally be obtained with pre-calculated load curves.
In general, PI tuning is a difficult task in practice, which is performed by field engineers who usually exploit conservative tuning methodologies. NMPC tuning can be performed offline through an automatic simulation campaign, with relevant savings of commissioning time. In this context, NMPC overcomes not only several issues related to overshoot managing, but also the PI controller mistuning issues, resulting in superior performance.

Conclusions
In this paper, an extensive presentation of a Nonlinear Model Predictive Control approach has been presented and compared with respect to the standard control methodologies for speeding up the power startup of steam turbines. The NMPC approach includes a detailed thermodynamical model of the machine that allows to identify with a great accuracy the thermal behavior of its rotor, in particular focusing on the equivalent stress.
A final analysis of the obtained results and the comparison of the proposed advanced control approach to standard ones show that relevant improvements can be achieved by using NMPC architecture. Firstly, the NMPC-when well-tuned with offline proceduresguarantees a total rejection of the stress overshoots and superior performances in terms of produced energy, with respect to both start-up through pre-calculated curves and the PI controller. In order to prove the effectiveness of the proposed approach, the application to a specific case study was presented regarding the startup of steam turbines in concentrated solar power plants, characterized by particularly variable steam conditions. For this specific application, when the PI controller action is tuned to eliminate stress overshoots, the NMPC gains about 2.3% in terms of energy produced during the startup, on top of the 15-20% output power gain that can be obtained through a standard PI control scheme aimed at thermal stress containment with respect to standard start-up curves.
The only bottleneck in applying NMPC methodologies when starting small machines is the calculation time, which currently does not guarantee application on PLC platforms, except through complex distributions of the optimization calculation on consecutive control cycles. Real-Time Iteration architecture may overcome these issues and will be the object of future developments of this topic. Future work will aim at developing an adaptive blocking strategy, based on the current thermal state and the stress of the machine. In particular, the strategy can be effectively adjusted through Fuzzy Logic-based methodologies that can vary both the number of DOFs and their concentration over time in order to favor a better optimization of the starting curve in its initial moments and a further reduction of the calculation times. A further increase in performance can be obtained by improving the disturbance models that forecasts the steam characteristics at the inlet of the turbine, along the prediction and control horizons of the NMPC. An accurate model of the available steam energy at the inlet is necessary to improve the MPC performances, as it allows one to obtain less restrictive and conservative control curves. Furthermore, a broader point of view on the startup issues faced in a global way could allow to further increase the performance. On the one hand, the control system can be extended to the generation of steam itself, with enormous benefit on the synchronization of the machines and the plant, at the cost of a higher computational burden. On the other hand, the modeling work can be extended in the same point of view, with a greater number of details and consequent precision of the predictions of the vapor characteristics. AI techniques, e.g., Recurrent Neural Networks-based architectures, can undoubtedly improve the prediction of future disturbances and make the start-up curve less conservative. Interesting examples are systems based on concentrated solar thermal technology, in which AI methodologies can certainly benefit from the use of AI techniques for prediction of steam characteristics based on weather data and the condition of the solar field plant, in order to improve not only the management of connected STs but also the overall performance. Funding: The work presented in this paper was developed within the project entitled "Smart Turbine Technologies" (STECH-CUP No. 745529), which has been developed in response to the call FAR-FAS for the financing of fundamental research, industrial research and experimental development projects jointly carried out by companies and research bodies in the field of new technologies in energy, photonics, ICT, robotics, and other related enabling technologies.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.