#### Abstract

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

#### 1. 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 efficiency improves greatly. Recently, predictive control has been successfully applied to thousands of process control systems and achieved great economic efficiency.

In most modern industrial process control techniques, steady-state target calculation (SSTC) is incorporated into MPC to find the optimal set-point automatically; thus, a double-layered MPC is established [3–6]. In [7], the double-layered 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 different rates. Moreover, the offset-free control is achieved even in the presence of disturbance [9, 10]. Therefore, the double-layered control has become the most effective 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 difficult for MPC to seek a satisfactory solution within an acceptable computational burden. Usually, the large-scale system is first decomposed into low order interconnected subsystems, and then a distributed control scheme is designed based on the subsystems [11–16]. In the distributed architecture, the interconnection influences among subsystems should be considered [11, 13, 17–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 large-scale 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–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 large-scale 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 real-time 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 real-time environments. Therefore, by adopting this strategy, the steady-state target of each subsystem in the distributed double-layered DMC can be solved effectively.

This 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. The 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. The 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.

#### 2. Distributed Open-Loop Prediction Model

Consider a large-scale system with outputs and inputs, its dynamics is described by the following finite step response (FSR) model:where is the modeling horizon. and are the controlled variable and manipulated variable, respectively. is the disturbance variable. , where denotes the value of at the equilibrium. For simplicity, is omitted. and are the step response coefficient matrices of the manipulated variable and disturbance variable, respectively. For all , , .

##### 2.1. Decomposed Systems

Let system equation (1) be decomposed into interconnected subsystems, where each subsystem has as input vector, i.e., and . Similarly, it has as output vector, i.e., , and . In the decomposition, step response coefficient matrices of subsystems are diagonal blocks of . While the nondiagonal blocks of (i.e., ) represent the interconnection influences among subsystems. means that the input of subsystem affects the dynamics of subsystem . Therefore, the FSR model of each subsystem is expressed as follows:where are the output and input vectors of subsystem , respectively. denotes the modeling horizon for the large-scale system. denotes the control horizon. denotes the neighbors of subsystem , which excludes , namely, denoting the collection of the subsystem ’s dynamic neighbor subsystems. represents the interconnection influence on subsystem from neighbor subsystem . is the disturbance variable.

##### 2.2. Distributed Open-Loop Prediction

A large-scale system is decomposed into dynamically coupled subsystems, and the distributed control is established based on the subsystems. The following symbols are used to distinguish different prediction information: is the real prediction of input is the assumed predictive value of input is the optimal predictive value of input

For interconnected subsystems, the open-loop prediction output of subsystem is affected not only by the input of local subsystem itself but also by some information stemming from its neighbors. For free prediction, the local input of each subsystem does not change from time onwards, that is, for all . However, the inputs to the neighbor subsystems keep changing. Due to communication delay, the inputs from its neighbors cannot be obtained immediately. The assumed value is used. Here, the assumed predictive input of each subsystem at time is chosen by the following formula:

Suppose that the inputs from time onwards do not change, , the corresponding prediction of is called free prediction, denoted as , where the superscript means free. In the distributed structure, the free prediction of each subsystem is established by assuming that the input of local subsystem is unchanged from time , whereas the inputs from neighbor subsystems may vary. Thus, the free prediction of subsystem can be expressed as follows:

The free prediction at time is

Suppose that is measurable and will remain unchanged in the future. Then, one has

Suppose that the control horizon of neighbor subsystem is . coupling inputs from neighbor subsystem will keep changing from time 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 and then the real prediction for subsystem can be calculated by

Seen from equation (8), at time , the real output prediction is determined by only the effect of subsystem itself. Thus, the coupling influences from neighbors can be neglected in the subsequent controller design. Notice that, in equations (6) and (8), the prediction step satisfies . By taking , it is allowed to take . Then, we have .

Feedback correction is added to equation (8) and then the open-loop prediction is obtained. Here, the superscript “” denotes open-loop. Based on equation (6), the open-loop prediction of subsystem can be calculated by

At time , the actual input of each subsystem is measurable. Suppose is measurable. The open-loop prediction error is obtained by

Let , then . Let , and can be considered as a bounded disturbance. The feedback correction is added into the open-loop prediction equation to compensate the prediction errors. Let , the submodel of subsystem can be written in vector form as follows:

Notice that the assumption is utilized in deducing equation (11).

For a stable process, the feedback correction is time-invariant, that is, , for . Then, the open-loop steady-state prediction of each subsystem can be expressed as

When replacing with 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. Thus, the open-loop prediction equation can predict the dynamics of each subsystem effectively.

#### 3. 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, steady-state 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 can be calculated bywhere denotes the increment variable of the steady-state target of for subsystem , . denotes the steady-state target value of controlled variable for subsystem .

##### 3.1. Handling the Constraints in SSTC

In practice, some operation requirements and physical restrictions need to be met. The 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.

The steady-state target of each manipulated variable needs to satisfy the corresponding hard magnitude constraints:where denotes the lower bound of and denotes the upper bound of .

In dynamic control module, the increment of manipulated variable has rate constraint and then should further satisfywhere denotes steady-state target increment of at time . denotes the maximum allowable absolute change rate of .

Combining equations (14) and (15), the hard constraint of iswhere

The 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 denotes the lower engineering limit and denotes the upper engineering limit.

Besides, the controlled variable target of each subsystem has to satisfy the operating limit:where denotes lower operating limit of and denotes upper operating limit of .

The engineering limits are more stringent than the operating limits and .

Substituting the steady-state prediction equation (13) into equation (19), we yield

Then, we havewhere

The ETs of and are denoted as and , respectively. , and represent the corresponding set of all the ETs of and . Denote the expected allowable range of as . The limits of ET values for can be taken as follows:

Then, for , check the following soft constraint:

Denote the expected allowable range of as . The limits of ET values for are given by

If the external target of is considered, the following soft constraints need to be further satisfied, that is,

Unifying all the ETs into the constraint of , the corresponding constraint for is denoted as

For , the constraint is

##### 3.2. Feasibility Analysis

For each subsystem , 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)):

The 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. The 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 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 rank of priorities, the constraint handled at the th priority can be denoted as

Constraint equations (30) and (31) are treated as hard constraints in the th priority. The following two cases are considered:(i)If the th rank of priority is the desired value of external target, relax the th equality constraint and then the considered constraints are as follows:(ii)If the th rank of priority is the desired upper bound or lower bound of target value, relax the th rank of inequality constraint and then the considered constraints are as follows: where and 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. Then, an equal concern error is assigned to each slack variable . For the upper bound constraints of controlled variables, , and for the lower bound constraints, , and for the ET constraints 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 and can be obtained by solving the following LP: where indicates the th element of and denotes the dimension of .

##### 3.3. 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 and . Suppose the steady-state change of input variable is bounded by

Then, minimizing the increment variable is equivalent to minimizing its constraint bound . 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 is written aswhere is the weighting coefficient.

Notice that 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 of each subsystem is not necessarily optimal for the whole large-scale system. Define the steady-state target vector as . According to the economic optimization problem of each subsystem, the global optimization problem can be chosen aswhere represents the performance index of each subsystem. Their cost functions are assumed to be Lipschitz continuous. Seen from (39), it is a typical multiobjective optimization problem. The Pareto optimality strategy is utilized to handle the multiobjective problem with conflicting cost functions in [29]. The Pareto optimal solution is obtained by Definition 1.

*Definition 1. *(steady-state Pareto solution [29]). A feasible solution of multiobjective optimization problem (39) is said to be Pareto optimal if and only if there exists no other feasible solution such that , and for at least one index .

The 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. The utopia-tracking MPC [28] is utilized to handle this problem.

*Definition 2. *(steady-state Utopia point [30]). The steady-state utopia point is the ideal optimal solution under the condition of , with . is given by solving the following problem:where the utopia cost of each subsystem is denoted as .

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. The compromise solution is defined as , which can be easily calculated by the following distance problem:Equation (41) is regarded as the steady-state utopia-tracking problem [28]. The 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, .

The compromise solution calculated by equation (41) is taken as the steady-state target value of subsystem . Substituting into equation (13), the set-point of is obtained. Then, the steady-state value is sent to dynamic control.

#### 4. 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 . For each subsystem , is given. Combining the real prediction equation (8) and the open-loop dynamic prediction equation (11), the closed-loop prediction equation of each subsystem is obtained:where is called the dynamic matrix. is the prediction horizon, and . 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 , the cost function is chosen aswhere the weighting matrix is given off-line. , updated online.

Besides, an additional constraint to limit is added to the dynamic control optimization problem:where . In addition, the inequality constraints including manipulated variable magnitude, manipulated variable increment, and controlled variable magnitude constraints are considered:where denotes the th row of dynamic matrix :

At time , the following optimization problem of each subsystem is firstly solved:

If the optimization is infeasible, slack variables and are added into the cost function and then the optimization problem is modified towhere and are slack variables. The 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:

Equations (48) and (52) are standard quadratic problem (QP). By solving equations (48) or (52), is obtained. Then, is sent to each subsystem and implemented.

*Remark 1. *In this paper, the considered distributed structured matrix dynamic control (DMC) is a heuristic algorithm. The dynamic is described by the finite step response (FSR) model, which is not accurate. The model error is compensated in a feedback correction term. The stability has not been stressed in the study of DMC since it was proposed by Garcia in 1986. The distributed double-layered DMC provides a solution for the predictive control of complex industrial systems with multiobjects and multilayered rates. The 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.

#### 5. Case Study

The 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. The 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. The 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:wherewhere the manipulated variables are the top extraction , the side extraction , and the enthalpy controller for the bottom of the column . are the middle reflux and the top reflux, respectively. The controlled variables include top product concentration , side product concentration , and bottom reflux temperature . All the variables are normalized, and the unit of time is minute.

As seen from Figure 2, by applying the interconnection decomposition, the transfer function in equation (53) can be decomposed into the following three subsystems: Subsystem 1: Subsystem 2: Subsystem 3:

In the distributed structure, the coupling effects from neighbor subsystems are considered in the control of each subsystem. The sampling period is 4 minutes. The finite step response model of each subsystem can be obtained bywhere denotes the number of the subsystem. and are the step response coefficient matrices on manipulated variables for local subsystem and neighbor subsystem , respectively. Denote and as the steady-state step response coefficient matrices.

Based on the open-loop prediction of each subsystem , the corresponding steady-state prediction is obtained bywhere denotes the steady-state increment of manipulated variable.

The constraints on manipulated variables and controlled variables are specified by ; ; ; ; ; ; . The external targets of are 0.5. By invoking LP, the steady-state target value of each subsystem is solved by the following optimization problem:

Due to the conflicts of the optimization goals of different subsystems, the steady-state target value of each subsystem cannot be achieved at the same time. The utopia-tracking strategy [28] is adopted to solve the multiple optimization problem with conflicts. According to equations (40) and (41), the compromise solution of steady-state target is calculated. In the dynamic matrix control, by utilizing QP, the actual control move for each subsystem is obtained.

The first investigation is carried out with the performance of the proposed distributed double-layered DMC algorithm with constant disturbance. At time , the value of disturbance is [0.15; 0.2]. In the distributed algorithm, the prediction and the control horizon for distributed double-layered DMC are 8 and 5, respectively. The weighting matrices in dynamic control are selected as , . Take , . The closed-loop system responses of three subsystems are shown in Figures 3–5, respectively.

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 Pareto-tracking 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. The 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 , the value of disturbance is [0.15; 0.2]. For the centralized algorithm, by tuning parameters, , the tracking performance is best. The 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 steady-state 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.

#### 6. 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 open-loop prediction equation of each subsystem. The distributed economic optimization problem is treated as a multiple objective optimization problem. Then, the utopia-tracking strategy is adopted to calculate the optimal steady-state target values of each subsystem. The proposed distributed double-layered DMC algorithm is demonstrated through simulation studies. The results show that the computational burden is significantly reduced, and the tracking performance is improved. This algorithm can also be extended to the control of other industrial processes.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Natural Key Research and Development Program of China under Grant 2018YFB1700100.