**1. Introduction**

The steam/water loop is an important part of a steam power plant, which plays a role in feed water supply and recycling processes. It is a highly complex and constrained system with multiple variables and interactions [1]. Meanwhile, due to the harsh and challenging operating environment (sea winds, sea waves and sea currents) and various operating modes (automatic start-up, reverse, stop, setting speed, emergency stop and reduction of revolutions) [2], there are difficulties to design a controller that delivers satisfactory performance for the steam/water loop. In order to design an effective approach overcoming the difficulties mentioned above and improving their liability and flexibility of the system, the feasibility of a distributed model predictive control is studied in this paper.

Nowadays, the major concern of the steam power plant is not only the tracking performance, but also other performances such as the consumed energy or safety in terrible conditions. Apart from realizing the load tracking performance, the controllers should also fulfill the flexible requirements for each sub-loop. A general way to improve the flexibility is to apply distributed controllers in the system [3]. Also, the multiple controllers improve the reliability of the system [4]. Concurrently, in order to deal with the constraints in the steam/water loop, model predictive control is selected to maintain

good performance for each sub-loop. In fact, such methods are readily matured for various applications in the Industry 4.0 paradigm [5]. For the particular class of complex processes with varying time delays, the simple loop control has been proposed and successfully tested in real life processes [6]. However, the system under investigation in this paper is a multivariable process. Hence, the single loop control does not suffice to guarantee good performance of all sub-system loops.

In recent years, many advanced control strategies have been developed to improve the performance of steam power plants, such as: Advanced PID control [7–10]; backstepping control [11–13]; sliding mode control [14–18]; robust control [19,20]; adaptive disturbance rejection control [21–23]. As it has the advantage in dealing with the inputs and outputs constraints, a considerable amount of research has been conducted in the application of model predictive control (MPC) for complex systems in general and for steam power plants, in particular. To overcome the control problems arising from load disturbances and intrinsic nonlinearity, an extended state observer based fuzzy MPC was proposed for the ultra-supercritical boiler-turbine unit [24]. By incorporating both the plant-wide economic process optimization and regulatory process control into a hierarchical control structure, a hierarchical MPC was proposed [25]. A stable fuzzy model predictive controller was designed to solve the superheater steam temperature control problem in a power plant [26], in which the effect of modeling mismatched and unknown plant behavior variations were overcome by the use of a disturbance term and steady-state target calculator. Liu proposed an economic model predictive controller by directly utilizing the economic index of the boiler-turbine system as the cost function. The method realized the economic optimization as well as the dynamic tracking [27].

However, the methods mentioned above are mainly concerned with the loop control of drum steam pressure, electric power and steam–water density in the boiler-turbine system. Most of the concerns concern the tracking performance of the system installed on land. The steam power plant in this paper is installed on ships, and the steam/water loop is depicted in Figure 1 with five sub-loops: Drum water level control loop; deaerator water level control loop; deaerator pressure control loop; exhaust manifold pressure control loop; and condenser water level loop. Also, the steam/water loop in large scale ships has special characteristics compared with those installed on land: (i) Receiving more disturbances from the ocean waves; (ii) smaller capacity; (iii) multiple operation points [28]. Hence, effective control strategies are required to ensure the performance of the steam/water loop.

A linear extended predictive self-adaptive control centralized MPC framework was studied in our previous work. The effects of the predictive horizon and the control horizon on the control system were summarized [28,29]. It was concluded that the predictive horizon be selected according to the system dynamics and the control horizon be selected to ensure a good trade-off between the tracking performance and the computational complexity.

Today's industry processes, for instance, those in heavy industries (e.g., papermaking, steelmaking, petrochemical and power generation), are becoming more and more complex [30], and are generally composed of a couple of interconnected systems, and possess high nonlinear and stochastic dynamics [31]. Due to the interplay between high performance local control and interactions between subsystem to subsystem, the centralized model predictive control (CMPC) is largely viewed as impractical, inflexible and unsuitable for control of such complex processes [32]. For example, as the authors have five sub-loops in the complete steam/water process, it is necessary to ensure the reliability of the system. This may be achieved by localizing control functions near each sub-loop with remote monitoring and supervision (i.e., distributed control) instead of one central controller (i.e., centralized control). Hence, the applicability of the distributed model predictive control (DiMPC) is investigated in this paper. A comparison is conducted between the decentralized model predictive control (DeMPC), the CMPC and the DiMPC, and the results show the effectiveness of the DiMPC.

The convergence issue is discussed to ensure the stability of the DiMPC. It is proved that the DiMPC, in our study, satisfies certain conditions, under which the properties of the CMPC problem are enjoyed by the DiMPC. Then, by the means provided in [33], the stability of the CMPC is proved, meanwhile, the stability of the DiMPC is guaranteed.

The rest of the paper is organized as follows: Section 2 describes the steam/water loop, and the modeling of the system is shown. The details about the DiMPC, CMPC and DeMPC are introduced in Section 3. Section 4 shows the results and analysis. Finally, some conclusions are drawn in Section 5.

**Figure 1.** Scheme of a complete steam/water loop [28] (reproduced with permission from Zhao, S.; Maxim, A.; Liu, S.; De Keyser, R.; and Ionescu, C, Processes; published by MDPI, 2018).

#### **2. Description of the Steam**/**Water Loop**

As shown in Figure 1, the steam/water loop is composed of two main loops, in which one is the water loop indicated by green line and the other is steam loop indicated by red line. The system works as follows. Firstly, the water from the water tank goes to the condenser. Then it is deoxygenated in the deaerator. After being pumped to boiler, the feed water goes into the mud drum due to its high density. The feed water is turned into a mixture of steam and water in the risers. Following this, the steam is separated from the mixture and heated in the superheater. Finally, the steam with a certain pressure and temperature services in the steam turbine. The used steam is sent back to the exhaust manifold and most of the steam is condensed in the condenser, while the remaining part services in the deaerator for deoxygenation [28]. The sub-loops have strong interactions between each other, such as the water level between the deaerator and the condenser, the pressure between the deaerator and the exhaust manifold system. Hence, there are challenges to obtain a desired controller for the steam/water loop.

In order to explore the characteristics of the steam/water loop, staircase experiments are conducted around the operating point on the system. The normalized outputs and corresponding static gains are shown in Figure 2; Figure 3, respectively. In the experiment, every 10% step changes are imposed in one input variable, while keeping the other inputs constant. The results show that the static gains change considerably along with the input changes, which indicates the nonlinearity of the system.

**Figure 2.** The normalized outputs of the steam/water loop with 10% step change in the input variables.

*Processes* **2019**, *7*, 442

**Figure 3.** Static gain under different input changes. (The figures on left hand indicate the static gain for drum water level loop, deaerator water level loop, and condenser water level loop; and on the right hand indicate the static gain for deaerator pressure loop, and exhaust manifold pressure loop).

By linearization around the operating point, the model of the system is obtained as shown in (1), with five inputs and five outputs. The input vector **u** = [*u*1*, u*2*, u*3*, u*4*, u*5] contains the positions of the valves that control the flow rates of feedwater to the drum (*u*1), exhaust steam from the exhaust manifold (*u*2), exhaust steam to the deaerator (*u*3), water from the deaerator (*u*4) and water to the condenser (*u*5), respectively. The output vector **y** = [*y*1*, y*2*, y*3*, y*4*, y*5] contains the values of the water level in drum (*y*1), pressure in exhaust manifold (*y*2), water level (*y*3) and pressure (*y*4) in deaerator, and water level of condenser (*y*5), respectively. The ranges and operating points of the output variables are listed in Table 1, and the operating points are obtained according to a real large-scale ship.

$$
\begin{bmatrix} y\_1 \\ y\_2 \\ \vdots \\ y\_5 \end{bmatrix} = \begin{bmatrix} G\_{11} & G\_{12} & \cdots & G\_{15} \\ G\_{21} & G\_{22} & \cdots & G\_{25} \\ \vdots & \vdots & \ddots & \vdots \\ G\_{51} & G\_{52} & \cdots & G\_{55} \end{bmatrix} \begin{bmatrix} u\_1 \\ u\_2 \\ \vdots \\ u\_5 \end{bmatrix} \tag{1}
$$


**Table 1.** The parameters used in the steam/water loop.


The rates and amplitudes of the five inputs are constrained to:

$$\begin{cases} -0.007 \le \frac{du\_1}{dt} \le 0.007 & 0 \le u\_1 \le 1\\ -0.01 \le \frac{du\_2}{dt} \le 0.01 & 0 \le u\_2 \le 1\\ -0.01 \le \frac{du\_3}{dt} \le 0.01 & 0 \le u\_3 \le 1\\ -0.007 \le \frac{du\_4}{dt} \le 0.007 & 0 \le u\_4 \le 1\\ -0.007 \le \frac{du\_5}{dt} \le 0.007 & 0 \le u\_5 \le 1 \end{cases} \tag{2}$$

The inputs units are normalized percentage values of the valve opening (i.e., 0 represents a fully closed valve, and 1 is completely opened). Additionally, the input rates are measured in percentage per second.

#### **3. Control Strategies**

#### *3.1. Introduction for Centralized MPC (CMPC)*

The following is a short summary of the extended prediction self-adaptive control (EPSAC) and more details are found in [34]. Consider a discrete system described below:

$$y(k) = x(k) + w(k)\tag{3}$$

where *k* is the discrete-time index; *y*(*k*) indicates the measured output sequence of system; *x*(*k*) is the output sequence of model; and *w*(*k*) is the model/process disturbance sequence. The output of the model *x*(*k*) depends on the past outputs and inputs, and can be expressed generically as:

$$\mathbf{x}(k) \;= \; f[\mathbf{x}(k-1), \mathbf{x}(k-2), \dots, \mathbf{u}(k-1), \mathbf{u}(k-2), \dots] \tag{4}$$

In EPSAC, the future input consists of two parts:

$$
\mu(k+l|k) = \mu\_{\text{base}}(k+l|k) + \delta\mu(k+l|k) \tag{5}
$$

where *ubase*(*k* + *l <sup>k</sup>*) indicates *<sup>l</sup>*th predicted basic future control scenario and <sup>δ</sup>*u*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup> <sup>k</sup>*) indicates the optimizing future control actions, and they are all based on the states and inputs of time *k*. Then, following the results of *l*th predicted output is obtained by applying (5) as the control effort.

$$\log(k+l|k) = \lg\_{\text{base}}(k+l|k) + \jmath\_{\text{opt}}(k+l|k) \tag{6}$$

where *ybase*(*k* + *l <sup>k</sup>*) is the effect of base future control and *yopt*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup> <sup>k</sup>*) is the effect of optimizing future control actions δ*u*(*k <sup>k</sup>*), ... , <sup>δ</sup>*u*(*<sup>k</sup>* <sup>+</sup> *Nc* <sup>−</sup> <sup>1</sup> *<sup>k</sup>*). The part of *yopt*(*<sup>k</sup>* <sup>+</sup> *<sup>l</sup> <sup>k</sup>*) can be expressed as a discrete time convolution as follows:

$$y\_{opt}(k+l|k) = h\_l \delta u(k|k) + h\_{l-1} \delta u(k+1|k) + \dots + g\_{l-N\_c+1} \delta u(k+N\_c-1|k) \tag{7}$$

where *h*1, ... *hNp* are impulse response coefficients; *g*1, ... *gNp* are the step response coefficients; *Nc*, *Np* are the control horizon and prediction horizon, respectively. Thus, the following formulation is obtained:

$$
\mathbf{Y} = \overline{\mathbf{Y}} + \mathbf{G}\mathbf{U} \tag{8}
$$

$$\begin{array}{rclcrcl}\text{with}&\mathbf{Y} &=& \left[y(k+N\_1|k)\dots y(k+N\_p|k)\right]^T & \mathbf{U} &=& \left[\delta u(k|k)\dots \delta u(k+N\_c-1|k)\right]^T\\\mathbf{Y} &=& \left[y\_{\text{base}}(k+N\_1|k)\dots y\_{\text{base}}(k+N\_P|k)\right]^T & \text{and} \end{array}$$

$$\mathbf{G} = \begin{bmatrix} h\_{N\_1} & h\_{N\_1 - 1} & \dots & \mathcal{G}\_{N\_1 - N\_c + 1} \\ h\_{N\_1 + 1} & h\_{N\_1} & \dots & \dots \\ \dots & \dots & \dots & \dots \\ h\_{N\_P} & h\_{N\_P - 1} & \dots & \mathcal{G}\_{N\_P - N\_c + 1} \end{bmatrix} \tag{9}$$

where *N*<sup>1</sup> indicates the time-delay in the system.

The disturbance term *w*(*k*) is defined as a filtered white noise signal [30]. When there is no information concerning the noise, the disturbance model used in (3) is chosen as an integrator, to ensure zero steady-state error in the reference tracking experiment:

$$w(k) = w(k-1) + e(k)\tag{10}$$

where *e(k)* denotes the white noise sequence.

In order to apply EPSAC for a multiple-input and multiple-output (MIMO) system, the individual error of each output is minimized separately. The cost function for the steam/water system with five sub-loops is as follows:

$$J\_i = \sum\_{l=N\_1}^{N\_P} \left[ r\_i(k+l|k) - y\_i(k+l|k) \right]^2 (i = 1, 2, \dots, 5) \tag{11}$$

where *ri* (*i* = 1, 2, ... , 5) are the setpoints for the five loops.

By defining **G***ij* as the influence from *jth* input to *ith* output, (11) is rewritten as:

$$(\mathbf{R}\_i - \mathbf{Y}\_i)^T (\mathbf{R}\_i - \mathbf{Y}\_i) = (\mathbf{R}\_i - \overline{\mathbf{Y}}\_i - \sum\_{j=1}^5 \mathbf{G}\_{ij} \mathbf{U}\_j)^T (\mathbf{R}\_i - \overline{\mathbf{Y}}\_i - \sum\_{j=-1}^5 \mathbf{G}\_{ij} \mathbf{U}\_j) \tag{12}$$

with **R***<sup>i</sup>* denotes the reference for loop *i*, and **Y***<sup>i</sup>* denotes the predicted output for loop *i*.

Taking constraints from inputs and outputs into account, the process to find the minimum cost function becomes an optimization problem which is called quadratic programming.

$$\min\_{\mathbf{U}\_i} \mathbf{J}\_i(\mathbf{U}\_i) = \mathbf{U}\_i^T \mathbf{H}\_i \mathbf{U}\_i + 2\mathbf{f}\_i^T \mathbf{U}\_i + c\_i \text{ subject to } \mathbf{A}\mathbf{U} \le \mathbf{b} \tag{13}$$

$$\text{with}\begin{cases} \mathbf{H}\_{i} \mathbf{=} \mathbf{G}\_{ii}^{T} \mathbf{G}\_{ii} \mathbf{f}\_{i} \mathbf{=} - \mathbf{G}\_{ii}^{T} (\mathbf{R}\_{i} - \overline{\mathbf{Y}}\_{i} - \sum\_{j=1, j \neq i}^{5} \mathbf{G}\_{ij} \mathbf{U}\_{j}) \\\ c\_{i} = (\mathbf{R}\_{i} - \overline{\mathbf{Y}}\_{i} - \sum\_{j=1, j \neq i}^{5} \mathbf{G}\_{ij} \mathbf{U}\_{j})^{T} (\mathbf{R}\_{i} - \overline{\mathbf{Y}}\_{i} - \sum\_{j=1, j \neq i}^{5} \mathbf{G}\_{ij} \mathbf{U}\_{j})) \end{cases}$$

where **A** is a matrix; **b** is a vector according to the constraints and **U***<sup>i</sup>* is the input for sub-loop *i*.

Figure 4 shows the conceptual representation of the centralized MPC [35]. To get the optimal solutions for sub-loop *<sup>i</sup>*, the interaction *uj*<sup>∈</sup>*Ni* , *Ni* <sup>=</sup> *j* ∈ *N* : *Gij* - 0 from other sub-loops is taken into account as shown in (13).

Hence, the optimal centralized solution **U** = [**U1 U2 U3 U4 U5**] is obtained by solving the following global cost function:

$$J = \sum\_{i=1}^{5} p\_i I\_i \tag{14}$$

where *Ji* are defined in (11), and *pi* are weighting factors. In our case, the *pi* are chosen as the values which can normalize the cost function *Ji*.

**Figure 4.** A conceptual representation of the centralized model predictive control (MPC) architecture [33] (Reproduced with permission from Maxim, A.; Copot, D.; De Keyser, R.; and Ionescu, C.M., Journal of Process Control; published by Elsevier, 2018).

#### *3.2. Proposed Distributed MPC (DiMPC)*

It is noteworthy to mention that the centralized approach implies that all the information regarding to all the sub-systems (or sub-loops) is gathered in a single controller, as showed in Figure 4. The advantage is straightforward since the cost function (14) has an optimal solution. However, if one sub-loop malfunctions, then the entire steam/water loop collapses, with serious consequences for safety of the large scale ship.

One solution is provided by the distributed MPC (DiMPC) method, that regards all the sub-systems as independent modules which are controlled by an individual controller. Through the communication network, the inherent interactions are considered.

Thus, the same local cost function (13) is locally minimized by each controller, in which the coupling term <sup>5</sup> *j*=1,*ji* **G***ij***U***<sup>j</sup>* is computed with the input trajectory **U***<sup>j</sup>* received from the neighbors, and several iterations are performed until the local optimal solution is reached. For the sake of clarity, conceptual representation of the distributed MPC architecture is shown in Figure 5, and a pseudo-code is provided:

**Algorithm 1** The Iterative DiMPC

Step 1: Sub-loop *I* receives an optimal local control action δ**U***<sup>i</sup>* at the iterative time as *iter* = 0 according to the EPSAC, and the local control action δ**U***<sup>i</sup>* can be rewritten as δ**U***iter <sup>i</sup>* , where δ**U***<sup>i</sup>* indicates the vector of the optimizing future control actions with length of *Nci*;

Step 2: The δ**U***iter <sup>j</sup>* (*<sup>j</sup>* <sup>∈</sup> *Ni*) is communicated to the loop *<sup>i</sup>*, and the <sup>δ</sup>**U***iter*+<sup>1</sup> *<sup>i</sup>* is calculated again with the <sup>δ</sup>**U***iter j* from other loops;

Step 3: If the termination conditions δ**U***iter*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> <sup>δ</sup>**U***iter <sup>i</sup>* ≤ <sup>ε</sup>*<sup>i</sup>* <sup>∨</sup> *iter* <sup>+</sup> <sup>1</sup> <sup>&</sup>gt; *iter* are reached, the **<sup>U</sup>***iter*+<sup>1</sup> *<sup>i</sup>* is adopted, where ε*<sup>i</sup>* is the positive value and *iter* indicates the upper bound of the number of iteration times. Otherwise, the iter is set as *iter* = *iter* + 1, and return to Step 2;

Step 4: Calculate the optimal control effort as **U***<sup>t</sup>* = **U***base* + δ**U***iter*, and the control effort is applied to the system;

Step 5: Set *t* = *t* + 1, return to Step 1.

**Figure 5.** A conceptual representation of the distributed MPC architecture [33] (Reproduced with permission from Maxim, A.; Copot, D.; De Keyser, R.; and Ionescu, C.M., Journal of Process Control; published by Elsevier, 2018).
