*3.3. Classical Decentralized MPC (DeMPC)*

In Figure 6, the conceptual representation of the DeMPC is presented. When comparing with the distributed strategy from Figure 5, it can be seen that the main difference is given by the fact that the controllers do not exchange information, although the physical coupling remains.

**Figure 6.** A conceptual representation of the decentralized 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).

Hence, the local cost function to be minimized by each controller is:

$$\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{15}$$

$$\text{with} \begin{cases} \mathbf{H}\_i \mathbf{=} \mathbf{G}\_{ii}^T \mathbf{G}\_{ii} \mathbf{f}\_i = -\mathbf{G}\_{ii}^T (\mathbf{R}\_i - \overline{\mathbf{Y}}\_i) \\\ c\_i = \left(\mathbf{R}\_i - \overline{\mathbf{Y}}\_i\right)^T (\mathbf{R}\_i - \overline{\mathbf{Y}}\_i) \end{cases}$$

which is derived from (13), by removing the coupling influence between sub-loops.

### *3.4. Multiple Objective Distributed Model Predictive Control (MODiMPC)*

Nowadays, setpoint tracking is not the only target for the control system. For some fast dynamic systems, the computing speed of the control strategies has an important influence. In order to improve the computing speed, a small loss in tracking performance is made to realize the fast computing speed. The scheme of the MODiMPC is shown in Figure 7. There are three layers of structure, in which the priority is shown as safety > tracking performance > energy.

**Figure 7.** Flow chart for the MODiMPC.

The algorithm starts from initial zero conditions and computes the optimal control effort as an unconstrained solution of the optimization problem that aims to minimize the tracking error.

If the predicted inputs or outputs are not safe (generally out of the hard constraints), the constraints is included in the optimization to ensure safety. If the variables are in the safety interval, then according to the tracking error condition, two options are available, namely: (i) Error > ε, the focus is on performance, and the control effort is kept δ**U***iter <sup>i</sup>* , (i.e., the control effort is kept to the one which minimizes the cost function with (11)); or (ii) Error < ε, the focus is on energy, and the control effort is kept δ**U***iter <sup>i</sup>* = 0 (i.e., the actuator does not need to do any change), where ε is a tolerance error, and in this paper, it is chosen as 1% of the upper bound of the corresponding output. According to the end conditions (δ**U***iter*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> <sup>δ</sup>**U***iter <sup>i</sup>* ≤ ε*<sup>i</sup>* ∨ *iter* + 1 > *iter*) of the DiMPC, the procedure stops or continues to obtain a new result.

Due to the hydraulic cylinder being linked with the valve in the steam/water loop, the frequent changes in valves mean frequent changes in the hydraulic cylinder, which results in a large energy costs. In this sense, the energy is saved if the valves do not need to do any change under the condition Error < ε.

In the traditional MPC, the influence from constraints has always been considered to obtain the optimal inputs for the system. In our study, the quadratic programming is applied. However, the setpoint does not always change during the operation of the system, and most of the time, the system operates at a stable operating point. During this kind of period, the only thing to be considered is energy, in which the control effort is always kept the same as the last sampling time. Hence, no optimization process exists anymore, and there is a huge reduction in computing time.

#### *3.5. Convergence Issue*

In order to analyze the stability of the optimal solution of the distributed control system, first the convergence issue is discussed. A standard MPC formulation is written in a form of a series of static optimization problems shown as follows:

$$SP\_k: \min\_{\mathbb{S}} J(\mathbb{S}) \tag{16}$$

$$\begin{array}{l} \text{s.t.} \, M(\mathbb{S}) = 0 \\ \mathbb{C}(S) \le 0 \end{array}$$

where *S* is the vector of the decision variables, including the state variables *X* and the control variables *U*. *M(S)* is the prediction model and *C(S)* denotes the constraints. Although our process model has an input-output formulation, it can be easily translated into a state-space definition.

In the DiMPC, Equation (16) is decompositioned into subproblems. For the sub-loop *i*:

$$SP\_{ki}:\min\_{\mathcal{S}\_i} I\_i(\mathcal{S}\_{i\prime}\mathcal{S}\_i^{nri})\tag{17}$$

$$\begin{array}{l}\text{s.t. } M\_{\bar{l}}(\mathcal{S}\_{\bar{l}\prime}\mathcal{S}\_{\bar{l}}^{nri}) = 0\\C\_{\bar{l}}(\mathcal{S}\_{\bar{l}\prime}\mathcal{S}\_{\bar{l}}^{nri}) \le 0\end{array}$$

where *Si* is the vector of the decision variables for sub-loop *i*, and *Snei <sup>i</sup>* is the vector of the variables for the sub-loops that have interaction with sub-loop *i*. *Mi* and *Ci* meet the conditions: ∪*iSi* = *S*,∪*iMi* = *M*,∪*iCi* = *C*.

According to [35], if the distributed control methodologies satisfy some conditions, the properties that can be proved for the equivalent CMPC problem are enjoyed by the solution obtained using the DiMPC implementation. Also, the convergence issue of the DiMPC is equal to the CMPC. The optimal results from the DiMPC converge to the global optimal point. The conditions are listed as follows:


However, the conditions 2 and 3 are over strict as many systems are nonconvex and have nonlinearity in reality. Further, [36,37] show that these two conditions can be relaxed to nonconvex optimal problems with nonlinearity.

Moreover, the convergence of the DiMPC is further analyzed using the study given in [33]. Hence, starting from the unconstrained optimal solution of the distributed algorithm, the idea is to rewrite it with a recursive matrix formulation. After some matrix manipulation, a compact description is obtained:

$$\mathbf{U}^\*(k) = \mathbf{F}(k) - \mathbf{\hat{H}U}^\*(k-1) \tag{18}$$

where U<sup>∗</sup> (*k*) consists of the optimal sequences of all the sub-loops, computed at sample time *k*, while U∗ (*k* − 1) are the shifted optimal trajectories computed at the previous sampling time *k* − 1, with the last term doubled, to ensure the dimensions consistencies. The term F(*k*) is variable and computed at each sampling instant using the prediction error, whereas Hˆ is a constant term that is computed off-line in the initialization stage of the algorithm (see [33] for further details).

Note that using Equation (18), the convergence of the local optimal solutions can be checked by verifying that all the eigenvalues from H are inside the unit circle. Additionally, Equation (18) can be ˆ reformulated in the classical system approach as:

$$\mathbf{U}^\*(k) = (\mathbf{I} - q^{-1} \, \hat{\mathbf{H}} \,) \mathbf{F}(k) \tag{19}$$

where *q*−<sup>1</sup> is an operator that shifts the data backward one sampling period, F(*k*) is regarded as the system's input, while U<sup>∗</sup> (*k*) denotes the system's output. It is noteworthy to mention that Equation (19) can be used to analyze the stability of the optimal solution U<sup>∗</sup> (*k*) at sampling time *k*, in the classical linear time-invariant framework by verifying that all the eigenvalues of the system equivalent matrix (I−*q*−<sup>1</sup> Hˆ ) are inside the unit circle. Furthermore, if this condition is satisfied on the equality case, it results in the optimal solution of the distributed algorithm being marginally stable.

Hence, using this simple approach, the evolution of the system is computed using F(*k*) as the system's input, which is calculated using the prediction error from each sampling period. Although this is an analytical approach to recursively place the system's progression in time, it can be straightforwardly computed in an automatic manner, using the simulation tools available for a control engineer. Moreover, all the computations are computed in a distributed manner, since using Equation (18), each sub-loop computes the optimal trajectories of the coupling neighbors, and knowing this information, it computes its own optimal trajectory at each sampling instant.

#### **4. Simulation Results and Analysis**

According to our previous work, the parameter configuration for the EPSAC method is shown in Table 2.


**Table 2.** The parameters applied in various MPC controllers.

Where the *Ts* is the sampling time; *Nc*1, *Nc*2, ... , *Nc*<sup>5</sup> are control horizons; (the control horizons were selected by finding a good trade-off between tracking performance and computation time for each loop), *Np*1, *Np*2, ... , *Np*<sup>5</sup> are prediction horizons of the five loops, respectively (the prediction horizons were selected taking into account the specific transient dynamics for each loop); *Ns* is the number of the samples. The step setpoints are provided in Table 3. In the experiments, the initial condition was set at the operating point of the steam/water loop.

**Table 3.** Step setpoints changes in the experiments.


The simulation results are shown in Figure 8, including the system outputs and the corresponding control efforts. In order to test which case provides the best result, performance indexes in an average value for the five sub-loops were compared including integrated absolute relative error (*IARE*), integral secondary control output (*ISU*), ratio of integrated absolute relative error (*RIARE*), ratio of integral secondary control output (*RISU*) and combined index (*J*). These indexes are calculated with the following expressions:

$$IARE \;= \frac{1}{5} \sum\_{i=1}^{5} \sum\_{k=0}^{N\_s - 1} \left| r\_i(k) - y\_i(k) \right| / r\_i(k) \tag{20}$$

$$ISII = \frac{1}{5} \sum\_{i=1}^{5} \sum\_{k=0}^{N\_s - 1} \left( u\_i(k) - u\_{ssi}(k) \right)^2 \tag{21}$$

$$IRARE(C\_2, C\_1) = \frac{IARE(C\_2)}{IARE(C\_1)}\tag{22}$$

$$IRSI(\mathcal{C}\_2, \mathcal{C}\_1) = \frac{ISL(\mathcal{C}\_2)}{ISL(\mathcal{C}\_1)} \tag{23}$$

$$f(\text{C2}, \text{C1}) = \frac{w\_1 RIARE(\text{C2}, \text{C1}) + w\_2 RISI(\text{C2}, \text{C1})}{w\_1 + w\_2} \tag{24}$$

where *ussi* is the steady state value of the *i*th input; *C*1, *C*<sup>2</sup> are the compared controllers and the weighting factors *w*<sup>1</sup> and *w*<sup>2</sup> in (24) are chosen as *w*<sup>1</sup> = *w*<sup>2</sup> = 0.5.

As depicted in Figure 8, the CMPC had similar performance as the DiMPC, and both outperformed the DeMPC. This conclusion is not only valid for this process, but also for other processes, since the DeMPC strategy does not take the interactions into account, which lead to some severe fluctuations when the setpoint changes in other variables. Although the control efforts in the DiMPC are obtained separately, the performance is still good due to the iteratively communication between the controllers for each sub-loop. In real life operations, where the addition of the effects of noise and stochastic disturbances are needed, perhaps with adding periodic disturbances from sea dynamics, the DeMPC may even lead to instability in the overall system.

**Figure 8.** *Cont*.

**Figure 8.** *Cont*.

**Figure 8.** The responses of the steam/water loop under the DeMPC, CMPC and DiMPC for (**a**) drum water level control loop, (**b**) deaerator water level control loop, (**c**) deaerator pressure control loop, (**d**) condenser water level control loop and (**e**) exhaust manifold pressure control loop (The figures on left hand indicate the outputs, and on the right hand indicate the inputs).

The iterations are shown in Figure 9 during the optimization of DiMPC. In the algorithm, the two conditions to end the iteration are designed as: (i) The difference between consecutive optimal inputs fulfill the condition δ**U***iter*+<sup>1</sup> *<sup>i</sup>* <sup>−</sup> <sup>δ</sup>**U***iter <sup>i</sup>* ≤ 0.002 ; (ii) the maximum iteration time is five, i.e., *iter* > 5.

**Figure 9.** The iteration times when the DiMPC is applied to the steam/ water loop.

The same conclusion is also obtained according to the numerical values shown in Tables 4 and 5. As the index *J* implies, the DiMPC and CMPC have similar results. However, there is only one controller in the CMPC, which means the system may be out of service if there is any problem with the controller. On the contrary, the DiMPC has more ability in fault-tolerance and flexibility without much performance loss compared with the CMPC. As the DiMPC only needs a part of the entire model, it is much easier to find a feasible solution, while CMPC needs the entire model to obtain all the solutions at one time. In this context, the DiMPC has better robust performance than the CMPC. Hence, the system model required for the DiMPC can be less accurate than the CMPC. In an industry context, a staggering 60–70% of the project time is spent on model development, while the rest is claimed for the controller design and validation [38]. Hence, any reduction in identification time requirement greatly diminishes

overall loop maintenance control related costs. The analysis given in this paper provides a trade-off solution, yet with acceptable performance but with great yields for cost reduction in control design and validation time.


**Table 4.** Performance indexes for *IARE* and *ISU.*


The results for the DiMPC and the multiple objective distributed model predictive control (MODiMPC) are shown in Figure 10. The performance indexes are shown in Table 6. The computing time for MODiMPC is 2.81 s and for the DiMPC 29.36 s, respectively.

**Figure 10.** *Cont*.

**Figure 10.** The responses of the steam/water loop under the MODiMPC and DiMPC for (**a**) drum water level control loop, (**b**) deaerator water level control loop, (**c**) deaerator pressure control loop, (**d**) condenser water level control loop and (**e**) exhaust manifold pressure control loop (The figures on left hand indicate the outputs, and on the right hand indicate the inputs).


**Table 6.** Performance indexes for *IARE* and *ISU.*

It can be seen from the results that there is a large improvement in the computing time after the MODiMPC was applied and without too much loss in tracking performance.

As previously mentioned, in order to guarantee the stability of the optimal solution in a DiMPC framework, the convergence of the optimal solution was firstly discussed. Due to the fact that there are only constraints in input variables in our study, the feasibility of the steam/water loop belongs to the trivial case (according to [33], solving the optimization problem the existence of a feasible solution is ensured at each sampling period). The *H<sup>i</sup>* matrices for the five loops are calculated as follows (for more details about convergence issue, please refer to the reference [33]):

*H*<sup>1</sup> = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 0.9613 −0.0352 −0.0308 −0.0639 −0.0352 0.9673 −0.0293 −0.0716 −0.0308 −0.0293 −0.9730 −0.0781 −0.0639 −0.0716 −0.0781 0.2062 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ , *H*<sup>2</sup> = 0.0003, H3 = 0.0025, H4 = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 0.9866 −0.0125 −0.0113 −0.0224 −0.0125 0.9876 0.0114 −0.0298 −0.0113 −0.0114 −0.9888 −0.0369 −0.0224 −0.0298 −0.0369 0.4255 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ , H5 = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 0.9788 −0.0189 −0.0159 −0.0119 −0.0077 0.0145 −0.0189 0.9789 0.0189 −0.0160 −0.0121 0.0106 −0.0159 −0.0189 −0.9788 −0.0190 −0.0160 0.0041 −0.0119 −0.0160 −0.0190 0.9788 −0.0188 −0.0049 −0.0077 −0.0121 −0.0160 −0.0188 0.9792 −0.0158 0.0145 0.0106 0.0041 −0.0049 −0.0158 0.4403 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

By the eig function in MATLAB (R2016b, Mathworks, Natick, MA, USA, 2016), the eigenvalues are calculated in a centralized manner for the steam/water loop, and the maximum value is ρmax < 1, which indicates that the DiMPC is convergent. In the steam/water loop, the five sub-loops cover the full system, and in the DiMPC the information is exchanged iteratively. Hence, the conditions 1 and 3–6 are satisfied. In order to cover the worst circumstances, the sufficient conditions in Section 3.5 tend to be conservative, and in some cases the convexity are not necessary [37]. Hence, it is concluded that the DiMPC has the same convergence as the CMPC, and the convergence of the DiMPC is guaranteed.

#### **5. Conclusions**

Regarding the multiple sub-loops in the steam/water loop, this paper introduced a distributed model predictive control based on the EPSAC framework. Different types of the MPC were applied to the steam/water loop system, including the DeMPC, CMPC and DiMPC. According to the simulation results, the DiMPC had similar performance with the CMPC, and outperformed the DeMPC. Due to the multiple controllers in the DiMPC strategy, the DiMPC had better performance of fault-tolerance and flexibility than the CMPC which improved the reliability of the steam/water loop. By proving equivalence in stability between the DiMPC and the CMPC, and the stability of the CPMC, the stability of the DiMPC is guaranteed. Meanwhile, a multiple objective MPC was proposed, and the computing speed was improved without too much loss in tracking performance.

**Author Contributions:** Methodology, R.D.K., A.M. and S.Z.; software, S.Z., A.M. and S.L.; formal analysis, S.Z.; writing—original draft preparation, S.Z. and A.M.; writing—review and editing, S.Z., A.M., R.D.K., S.L., and C.M.I.; supervision, C.M.I., S.L. and R.D.K.; funding acquisition, S.Z., S.L. and C.M.I.

**Funding:** Shiquan Zhao acknowledges the financial support from Chinese Scholarship Council (CSC) under grant 201706680021 and the Co-funding for Chinese PhD candidates from Ghent University under grant 01SC1918. Part of this project is funded by a special research fund of Ghent University, MIMOPREC STG020-18 (Ionescu). This work is partly supported by the Fundamental Research Funds for the Central Universities under grant 3072019CF0408, 3072019CFT0403.

**Conflicts of Interest:** The authors declare no conflicts of interest.
