Distributed Double-Layered Dynamic Matrix Control for Large-Scale System

A distributed double-layered dynamic matrix control is proposed for large-scale systems. In the scheme, a large-scale system is rst decomposed into several low-dimensional subsystems, and then a double-layered dynamic matrix control algorithm with distributed structure is developed for these subsystems. e distributed open-loop prediction equation of each subsystem is formed based on the predicted output of each local subsystem and eects of its interconnecting neighbor subsystems. Due to simultaneous optimization, at each prediction, the coupling eects of neighbor subsystems are not available in time. us, the assumed value is utilized instead. In the economic optimization stage, conicts may occur among dierent economic optimization goals. Utopia-tracking strategy is introduced to optimize multiple steady-state targets. en, the obtained steady-state target values are taken as reference values and tracked in subsequent dynamic control. e actual control move for each subsystem is nally calculated. e proposed algorithm is tested on Shell heavy oil fractionator benchmark, and the eectiveness is demonstrated by comparing with the typical double-layered dynamic matrix control algorithm.


Introduction
With the advantage of handling various constraints explicitly [1,2], model predictive control (MPC) has become a standard method for solving the constrained multivariable control problem. Utilizing MPC, the production process can be operated on the boundaries of some important constraints, and the economic e ciency improves greatly. Recently, predictive control has been successfully applied to thousands of process control systems and achieved great economic e ciency.
In most modern industrial process control techniques, steady-state target calculation (SSTC) is incorporated into MPC to nd the optimal set-point automatically; thus, a double-layered MPC is established [3][4][5][6]. In [7], the doublelayered zone predictive control strategy was proposed to overcome the performance degradation caused by the mechanical wear. In [8], a method to guarantee the feasibility of optimization problem was presented for dynamic control tracking the steady-state targets. In the double-layered architecture, the economic optimization and dynamic control are separately working at di erent rates. Moreover, the o set-free control is achieved even in the presence of disturbance [9,10]. erefore, the double-layered control has become the most e ective and promising advanced industrial control technology.
However, it is challenge for the double-layered MPC to be applied to large-scale systems. As the size of the controlled process increases, the inherent computational complexity and communication limitation make it di cult for MPC to seek a satisfactory solution within an acceptable computational burden. Usually, the large-scale system is rst decomposed into low order interconnected subsystems, and then a distributed control scheme is designed based on the subsystems [11][12][13][14][15][16]. In the distributed architecture, the interconnection in uences among subsystems should be considered [11,13,[17][18][19]. Besides, in practical application, many physical limitations such as input saturation, the physical limitation of actuators, and external disturbances also need to be considered. In [15], a distributed optimal control was proposed for the interconnected large-scale system with external disturbances and input constraints, where by utilizing an adaptive dynamic programming technique, the worst-case disturbance policy and constrained-input optimal control policy can be approximate synchronously. An event-triggered feed-forward control policy was proposed in [16] to transform the control of largescale systems into equivalent event-triggered control of isolated subsystems. By designing a novel event-triggering condition, three neural networks are eliminated for each subsystem, so that the computational complexity and resource are reduced. A distributed noncooperative MPC was proposed in [20], where input and state constraints are considered and convergence of the closed-loop control system is proved. According to the optimization problem and the types of information exchange among subsystems, the distributed MPCs are often classified as cooperative [21][22][23] and noncooperative [19,20,24]. In cooperative distributed MPC, all subsystems solve a global optimization problem and exchange information with each other. While in noncooperative MPC, each subsystem solves a local optimization problem and exchanges information iteratively with its neighbor subsystems and finally achieves Nash optimal solution [12,24]. In [24], an efficient distributed MPC based on Nash optimality is presented for the largescale system, which mainly aimed at reducing the computational burden in model predictive control. Moreover, as reviewed in [6], both distributed and hierarchical control architectures are proposed for large-scale systems. In [25], the cooperative distributed MPC is extended to the hierarchical structure so as to reduce the communication frequency among subsystems. A novel distributed MPC algorithm for large-scale systems proposed in [26] is very similar to the MPC algorithms employed in industry, where each subsystem needs the reference trajectories only of the variables of its interconnected neighbors.
Motivated by the above mentioned, in this paper, a distributed double-layered dynamic matrix control (DMC) is proposed for large-scale systems. In the distributed structure, when some economic objectives are in conflict at the economic optimization stage among different subsystems, it is tough to solve the optimal steady-state target values. Multiobjective optimization strategy is introduced to solve this knotty problem. A typical method for solving the multiobjective optimization with conflicting cost functions is to construct its Pareto front and then choose the appropriate point as the optimal value [27]. However, this method is computationally heavy and impractical in realtime application. Besides, an ideal set-point varies with changing operating conditions; thus, the corresponding steady-state targets need to update at each iteration, which makes it more complicate to calculate the steady-state target values. In [28], the steady-state compromise solution was proposed to avoid the computation of Pareto fronts in realtime environments. erefore, by adopting this strategy, the steady-state target of each subsystem in the distributed double-layered DMC can be solved effectively.
is paper is organized as follows Section 2 presents the open-loop prediction model of each subsystem in the distributed scheme, where the estimation of coupling effects of neighbor subsystems is utilized. e error correction between the actual value and the prediction value is included in the open-loop dynamic prediction equation. Section 3 analyzes the feasibility of the local optimization of each subsystem. e optimal steady-state target values are calculated based on the improved Pareto optimality. Section 4 designs the dynamic control for the distributed system and then calculates the control move and sent to the plant and implemented. Section 5 demonstrates the effectiveness of the proposed distributed double-layered DMC algorithm by applying it to Shell heavy oil fractionator benchmark and comparing with the centralized double-layered DMC algorithm.

Distributed Open-Loop Prediction Model
Consider a large-scale system with n y outputs and n u inputs, its dynamics is described by the following finite step response (FSR) model: where N is the modeling horizon. y ∈ R n y and u ∈ R n u are the controlled variable and manipulated variable, respectively. f ∈ R n f is the disturbance variable. ∇y � y − y eq , where y eq denotes the value of y at the equilibrium. For simplicity, ∇ is omitted. S u i and S f t are the step response coefficient matrices of the manipulated variable and disturbance variable, respectively. For all t ≥ 0, S u 2.1. Decomposed Systems. Let system equation (1) be decomposed into M interconnected subsystems, where each subsystem a, a ∈ M, M � 1, 2, . . . , M { } has u a ∈ R n a as input vector, i.e., u � (u 1 , . . . , u M ) and M a�1 n ua � n u . Similarly, it has y a ∈ R n y as output vector, i.e., y � (y 1 , . . . , u M ), and M a�1 n y a � n y . In the decomposition, step response coefficient matrices S u 11 ∈ R n 1 ×n 1 , . . . , S u MM ∈ R n M ×n M of M subsystems are diagonal blocks of S u . While the nondiagonal blocks of S u (i.e., S u al , a ≠ l) represent the interconnection influences among subsystems. S u al ≠ 0 means that the input of subsystem l affects the dynamics of subsystem a.
erefore, the FSR model of each subsystem a ∈ M is expressed as follows: where y a ∈ R n y a , u a ∈ R n u a are the output and input vectors of subsystem a, respectively. N denotes the modeling horizon for the large-scale system. N c denotes the control horizon. l ∈ N a denotes the neighbors of subsystem a, which excludes a, namely, N a � 1, 2, . . . , a − 1, a + 1, M { } denoting the collection of the subsystem a's dynamic neighbor subsystems. S u al represents the interconnection influence on subsystem a from neighbor subsystem l. f ∈ R n f is the disturbance variable.

Distributed Open-Loop Prediction.
A large-scale system is decomposed into M dynamically coupled subsystems, and the distributed control is established based on the subsystems. e following symbols are used to distinguish different prediction information: For interconnected subsystems, the open-loop prediction output of subsystem a is affected not only by the input of local subsystem a itself but also by some information stemming from its neighbors. For free prediction, the local input of each subsystem does not change from time k onwards, that is, Δu a (t|k) � 0 for all t ≥ k. However, the inputs u l to the neighbor subsystems keep changing. Due to communication delay, the inputs from its neighbors cannot be obtained immediately. e assumed value is used. Here, the assumed predictive input Δu l (t|k) of each subsystem l at time k is chosen by the following formula: t � 1, 2, . . . , k − 1. (3) Suppose that the inputs from time k onwards do not change, Δu(k + t − 1|k) � 0 (1 ≤ t ≤ p), the corresponding prediction of y(k + p) is called free prediction, denoted as y fr (k + p|k), where the superscript ″ fr ″ means free. In the distributed structure, the free prediction of each subsystem a is established by assuming that the input of local subsystem is unchanged from time k, whereas the inputs from neighbor subsystems may vary. us, the free prediction of subsystem a can be expressed as follows: e free prediction at time k − 1 is Suppose that f(k) is measurable and will remain unchanged in the future. en, one has Suppose that the control horizon of neighbor subsystem l is N c . M coupling inputs from neighbor subsystem l, l ∈ N a will keep changing from time k onwards. According to equation (3), substituting the assumed values into equation (6), the one-step prediction error is obtained by As shown in equation (4), the influences of interconnection inputs from neighbors are considered in free prediction y fr a (k + p|k) and then the real prediction for subsystem a can be calculated by Seen from equation (8), at time k, the real output prediction is determined by only the effect of subsystem a itself.
us, the coupling influences from neighbors can be neglected in the subsequent controller design. Notice that, in equations (6) and (8), the prediction step p satisfies p ≤ 0.001. By taking S f 0 � 0, it is allowed to take p ≤ 0.001. en, we have y fr (k|k) � y(k|k) ≠ y(k).

Mathematical Problems in Engineering 3
Feedback correction is added to equation (8) and then the open-loop prediction y ol a is obtained. Here, the superscript "ol" denotes open-loop. Based on equation (6), the open-loop prediction of subsystem a can be calculated by e open-loop prediction error is obtained by Let i ε l , and e i can be considered as a bounded disturbance. e feedback correction ϵ a (k) is added into the open-loop prediction equation to compensate the prediction errors. Let N ′ � N + N c , the submodel of subsystem a can be written in vector form as follows: Notice that the assumption y ol a (k + N ′ |k − 1) � y ol a (k + N ′ − 1|k − 1) is utilized in deducing equation (11).
For a stable process, the feedback correction is timeinvariant, that is, en, the open-loop steady-state prediction of each subsystem can be expressed as When replacing k with k − 1 in equation (12), Equation (11) is still highly consistent with equation (12). In this case, the final open-loop dynamic prediction is equivalent to the steady-state prediction.
us, the open-loop prediction equation can predict the dynamics of each subsystem effectively.

Calculating Steady-State Target Values Based on Improved Pareto Optimality
For the control of industrial process, the manipulated variables and controlled variables do not always hit the desired external targets given by the real-time optimization (RTO). In order to ensure that the manipulated variables operate as close to the desired set-point as possible, steadystate optimization is incorporated into the DMC algorithm, thus the double-layered DMC algorithm is formed, in which the optimal set-point is calculated by steady-state target calculation (SSTC). SSTC involves two stages: the first stage aims at to analyze the feasibility of the solution which satisfies various constraints and the second one is economic optimization stage, which is designed to seek the optimal steady-state target values. Based on equation (8), the steady-state prediction equation for each subsystem a can be calculated by where δu a,ss denotes the increment variable of the steadystate target of u for subsystem a, δu a,ss (k) � u a,ss (k) − u a (k − 1). y a,ss denotes the steadystate target value of controlled variable y for subsystem a.

Handling the Constraints in SSTC.
In practice, some operation requirements and physical restrictions need to be met. e constraints imposed on controlled variables and manipulated variables include soft constraints and hard constraints. In general, the hard constraints cannot be violated, while the soft constraint can be relaxed appropriately according to the situations. Besides, some variables have their own external target (ET) values given by the real-time optimization control engineers or operators. e steady-state target of each manipulated variable needs to satisfy the corresponding hard magnitude constraints: where u denotes the lower bound of u and u denotes the upper bound of u.
In dynamic control module, the increment of manipulated variable has rate constraint |Δu a (k + j|k)| ≤ Δu(0 ≤ j ≤ M − 1) and then δu a,ss should further satisfy where δu ss (k) � u ss (k) − u(k − 1) denotes steady-state target increment of u at time k. Δu denotes the maximum allowable absolute change rate of u.
Combining equations (14) and (15), the hard constraint of δu a,ss is where e steady-state target of the controlled variables in each subsystem needs to satisfy the engineering limits, referred to as hard magnitude constraints as follows: where y 0,h denotes the lower engineering limit and y 0,h denotes the upper engineering limit. Besides, the controlled variable target of each subsystem has to satisfy the operating limit: where y 0 denotes lower operating limit of y and y 0 denotes upper operating limit of y. e engineering limits are more stringent than the operating limits y 0,h ≤ y 0 and y 0, h ≥ y 0 .
Substituting the steady-state prediction equation (13) into equation (19), we yield en, we have where e ETs of u i,ss and y j,ss are denoted as u i,t and y j,t , respectively. i ∈ I t , j ∈ J t , I t and J t represent the corresponding set of all the ETs of u i,ss and y j,ss . Denote the expected allowable range of u i,t (k) as u i,ss,range . e limits of ET values for u i,ss can be taken as follows: en, for u j,ss (k), check the following soft constraint: Denote the expected allowable range of y j,t (k) as y j,ss,range . e limits of ET values for y j,ss are given by If the external target of y j,ss (k) is considered, the following soft constraints need to be further satisfied, that is, Unifying all the ETs into the constraint of δu ss , the corresponding constraint for u i,t (k) is denoted as For y j,t (k), the constraint is

Feasibility Analysis.
For each subsystem a, the feasible solution of the economic optimization is determined by its constraint conditions including hard constraints and soft constraints, where the hard constraints are classified into manipulated variable constraints (equation (16)) and hard controlled variable constraints (equation (18)): e soft constraints do not have to be satisfied. In practical application, priority strategy is usually introduced to handle soft constraints. All soft constraints are prioritized according to importance. e rule for handling the soft constraints is that the soft constraints with higher priority are processed first and that with lower priority later. Soft constraints that are previously processed (possibly relaxed) are treated as hard constraints in subsequent priority processing. Some soft constraints are relaxed properly so that all Mathematical Problems in Engineering the constraints at the same priority are compatible. For convenience, at each rank of priority, either equality constraints or inequality constraints are dealt with. For the optimization problem with r rank of priorities, the constraint handled at the (r − 1) th priority can be denoted as Constraint equations (30) and (31) are treated as hard constraints in the r th priority. e following two cases are considered: (i) If the r th rank of priority is the desired value of external target, relax the r th equality constraint and then the considered constraints are as follows: (ii) If the r th rank of priority is the desired upper bound or lower bound of target value, relax the r th rank of inequality constraint and then the considered constraints are as follows: where ε (r) eq (k) and ε (r) (k) are the slack variables. In order to seek a smaller |ε|, either LP or QP is invoked. In each rank of priority, if more than one soft constraint needs to be relaxed, correspondingly, there are multiple ε's need to be optimized. Suppose given them equal importance. en, an equal concern error ε is assigned to each slack variable ε. For the upper bound constraints of controlled variables, ε � y h − y 0 , and for the lower bound constraints, ε � y 0 − y h , and ε � 1/2ET range for the ETconstraints or equality constraints. If some ε � 0, the corresponding soft constraint cannot be relaxed and then the corresponding soft constraint can be removed. Set then the slack variables ε (r) eq+ (k) and ε (r) eq− (k) can be obtained by solving the following LP: where τ indicates the τ th element of ε (r) eq (k) and d (r+1) eq denotes the dimension of ε (r) eq (k).

Economic Optimization Stage of SSTC Based on Pareto
Optimality. In the economic optimization stage, an economic cost function representing the profit/loss is constructed based on the change of u ss (k) and y ss (k). Suppose the steady-state change of input variable is bounded by en, minimizing the increment variable |δu i,ss (k)| is equivalent to minimizing its constraint bound U i (k). Having analyzed the feasibility, all the hard constraints and the relaxed soft constraints are simplified into the following formula: By invoking LP, the economic optimization problem of each subsystem a is written as (19), (24), (26)- (28), where h i is the weighting coefficient.
Notice that δu a,ss calculated by equation (38) is local optimal for each subsystem itself, but not for the whole large-scale system. In this distributed control, if some economic objectives of different subsystems are in conflict, the steady-state target value δu a of each subsystem is not necessarily optimal for the whole large-scale system. Define the steady-state target vector as u a,ss � [u 1,ss , u 2,ss , . . .
where [J 1 (δu 1,ss ), J 2 (δu 2,ss ), . . . , J M (δu a,ss )] represents the performance index of each subsystem. eir cost functions are assumed to be Lipschitz continuous. Seen from (39), it is a typical multiobjective optimization problem. e Pareto optimality strategy is utilized to handle the multiobjective problem with conflicting cost functions in [29]. e Pareto optimal solution is obtained by Definition 1.  e surface of the Pareto solution set is usually called the Pareto frontier. In practice, it is time-consuming to construct the Pareto frontier. Since the optimal point selected along the Pareto frontier will change with the conditions, the ideal Pareto optimum may be unreachable. e utopiatracking MPC [28] is utilized to handle this problem.
Definition 2 (steady-state Utopia point [30]). e steadystate utopia point is the ideal optimal solution (δu L a,ss ) under the condition of J a (δu L a,ss ), with a ∈ M. δu L a,ss is given by solving the following problem: where the utopia cost of each subsystem a is denoted as J L a . Considering the conflict of economic goals among subsystems, the utopia value is not reachable. Here, a compromise solution is used instead, which is the closest point from the Pareto frontier to the utopia value. e compromise solution is defined as δu C a,ss , which can be easily calculated by the following distance problem: s.t. (13), a ∈ M. (41) Equation (41) is regarded as the steady-state utopiatracking problem [28]. e schematic of the utopia-tracking approach is shown in Figure 1.
Notice that for the single-objective case, the compromise solution coincides with the utopia point, that is, J 1 (δu 1,ss ) � J L 1 . e compromise solution δu C a,ss calculated by equation (41) is taken as the steady-state target value of subsystem a. Substituting into equation (13), the set-point of y a,ss is obtained. en, the steady-state value (δu a,ss � δu C a,ss , y a,ss ) is sent to dynamic control.

Distributed Dynamic Control
In the dynamic control, update the actual control move at each time step so that the output predicted value closely tracks the steady-state value Δu C a,ss . For each subsystem a, Y ol a,N (k|k) is given. Combining the real prediction equation (8) and the open-loop dynamic prediction equation (11), the closed-loop prediction equation of each subsystem a is obtained: S a is called the dynamic matrix. N is the prediction horizon, and N c ≤ N ≤ N ′ . Equation (42) is the closed-loop prediction equation.
In order to ensure the output to track the steady-state targets as close as possible, meanwhile suppress the controlled variable increment δu a,ss , the cost function is chosen as where the weighting matrix Λ is given off-line.
Besides, an additional constraint to limit δu C a,ss (k) is added to the dynamic control optimization problem: where L � [I I · · · I]. In addition, the inequality constraints including manipulated variable magnitude, manipulated variable increment, and controlled variable magnitude constraints are considered:

Mathematical Problems in Engineering
where S i denotes the i th row of dynamic matrix S: max y a,0 , y a,ss (k) , , y a,ss (k) .

(47)
At time k, the following optimization problem of each subsystem is firstly solved: If the optimization is infeasible, slack variables ε dc (k) and ε dc (k) are added into the cost function and then the optimization problem is modified to where ε dc (k) and ε dc (k) are slack variables. e related constraints on slack variables are as follows: By adding the constraints equations (50) and (51), the optimization problem with slack variables is expressed as follows:  en, u a (k|k) � u a (k|k − 1) + Δu * a (k|k) is sent to each subsystem and implemented. Remark 1. In this paper, the considered distributed structured matrix dynamic control (DMC) is a heuristic algorithm. e dynamic is described by the finite step response (FSR) model, which is not accurate. e model error is compensated in a feedback correction term. e stability has not been stressed in the study of DMC since it was proposed by Garcia in 1986. e distributed double-layered DMC provides a solution for the predictive control of complex industrial systems with multiobjects and multilayered rates. e authors in this paper focus on the economic efficiency and computation burden, while the ingredients for proving stability such as the terminal constraint set, terminal cost function, and feedback control law are not mentioned. On the other hand, the optimization problem, in this paper, is solved by linear programming (LP) or quadratic programming (QP); therefore, the global optimal can be obtained.

Case Study
e distributed double-layered DMC algorithm is applied to the Shell benchmark problem for simulation investigation. In the operation of Shell heavy oil fractionator, the controlled variables and manipulated variables need to meet various requirements. As shown in Figure 2, there are three product extraction ports and three side cyclic refluxes. e product concentrations coming out of the top and side extraction ports are determined by the economics and operating requirements. For the cyclic refluxes, the product separation is performed by removing heat. e exchange heat flows back to reboil the columns in other parts of the plant. An enthalpy controller is added to regulate heat removal in the bottom while the heat exchanges of the other two refluxes are taken as the disturbances to the column.
A three-input-three-output model was presented in [31] for the heavy oil fractionator: where the manipulated variables are the top extraction u 1 , the side extraction u 2 , and the enthalpy controller for the bottom of the column u 3 . d 1 and d 2 are the middle reflux and the top reflux, respectively. e controlled variables include top product concentration y 1 , side product concentration y 2 , and bottom reflux temperature y 3 . All the variables are normalized, and the unit of time is minute.

Mathematical Problems in Engineering
As seen from Figure 2, by applying the interconnection decomposition, the transfer function G u (s) in equation (53) can be decomposed into the following three subsystems: Subsystem 1: 4.05e − 27s /50s + 1 Subsystem 2: 5.72e − 14s /60s + 1 Subsystem 3: 7.20/19s + 1 In the distributed structure, the coupling effects from neighbor subsystems are considered in the control of each subsystem. e sampling period is 4 minutes. e finite step response model of each subsystem can be obtained by Due to the conflicts of the optimization goals of different subsystems, the steady-state target value of each subsystem δ u a,ss cannot be achieved at the same time.
e utopia- tracking strategy [28] is adopted to solve the multiple optimization problem with conflicts. According to equations (40) and (41)  It can be observed that the proposed distributed dynamic prediction model can predict effectively the output of each subsystem at future time. In the distributed structure, the steady-state values of each subsystem obtained by Paretotracking satisfy all the constraints. When the calculated control move is implemented in each subsystem, even in the presence of disturbance, the corresponding output can be settled to the optimal steady-state value, which indicates that the control is offset-free. Meanwhile, all the outputs and control moves satisfy the physical constraints and operational restrictions.
e results indicate that the proposed distributed double-layered DMC algorithm is effective. Next, we will step into the study of the performance of the proposed algorithm. A comparison with the proposed distributed double-layered DMC algorithm and the centralized algorithm is carried out. In both of these algorithms, constant disturbance is considered, respectively. At time 55 ≤ k ≤ 90, the value of disturbance is [0.15; 0.2]. For the centralized algorithm, by tuning parameters, N � 8, N c � 5, the tracking performance is best. e simulation results are shown in Figure 6, and the performance is quantified in Table 1. As seen from Figure 6, though the output of each subsystem starts with a slight bigger error from the steadystate target, as the iteration goes on, there is an obvious improvement in the tracking performance of the distributed algorithm compared with that of the centralized algorithm.  As indicated by the time-consuming in Table 1, the proposed distributed algorithm requires much less computation time than the centralized one. As seen from Table 1, there is a (7.06 % increase in RMSE) improvement in tracking performance as a whole.

Discussion and Conclusion
A distributed double-layered DMC scheme for the large-scale system is presented, where the coupling effects from interconnected subsystems are considered in the distributed openloop prediction equation of each subsystem. e distributed economic optimization problem is treated as a multiple objective optimization problem. en, the utopia-tracking strategy is adopted to calculate the optimal steady-state target values of each subsystem. e proposed distributed doublelayered DMC algorithm is demonstrated through simulation studies. e results show that the computational burden is significantly reduced, and the tracking performance is improved. is algorithm can also be extended to the control of other industrial processes.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.