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

Authors: Li Wang; Yuanli Cai; Xin Zan
Journal: Mathematical Problems in Engineering
 (2022)

---

## 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.


---

## Body

## 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 Section2 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 withny outputs and nu inputs, its dynamics is described by the following finite step response (FSR) model:(1)∇yk=∑t=1N−1StuΔuk−t+SNu∇uk−N+∑t=1N−1StfΔfk−t+SNf∇fk−N,where N is the modeling horizon. y∈ℝny and u∈ℝnu are the controlled variable and manipulated variable, respectively. f∈ℝnf is the disturbance variable. ∇y=y−yeq, where yeq denotes the value of y at the equilibrium. For simplicity, ∇ is omitted. Siu and Stf are the step response coefficient matrices of the manipulated variable and disturbance variable, respectively. For all t≥0, SN+tu=SNu, SN+tf=SNf.


### 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 ua∈ℝna as input vector, i.e., u=u1,…,uM and ∑a=1Mnua=nu. Similarly, it has ya∈ℝny as output vector, i.e., y=y1,…,uM, and ∑a=1Mnya=ny. In the decomposition, step response coefficient matrices Su11∈ℝn1×n1,…,SuMM∈ℝnM×nM of M subsystems are diagonal blocks of Su. While the nondiagonal blocks of Su (i.e., Sual,a≠l) represent the interconnection influences among subsystems. Sual≠0 means that the input of subsystem l affects the dynamics of subsystem a. Therefore, the FSR model of each subsystem a∈M is expressed as follows:(2)∇yak=∑t=1N−1StuaΔuak−t+SNuauak−N+∑l∈Na∑t=1N−1StualΔulk−t+SNualulk−N′+∑t=1N−1StfΔfk−t+SNffk−N,where ya∈ℝnya,ua∈ℝnua are the output and input vectors of subsystem a, respectively. N denotes the modeling horizon for the large-scale system. Nc denotes the control horizon. l∈Na denotes the neighbors of subsystem a, which excludes a, namely, Na=1,2,…,a−1,a+1,M denoting the collection of the subsystem a’s dynamic neighbor subsystems. Sual represents the interconnection influence on subsystem a from neighbor subsystem l. f∈ℝnf is the disturbance variable.


### 2.2. Distributed Open-Loop Prediction

A large-scale system is decomposed intoM dynamically coupled subsystems, and the distributed control is established based on the subsystems. The following symbols are used to distinguish different prediction information:uk+t|k is the real prediction of inputu^k+t|k is the assumed predictive value of inputu∗k+t|k is the optimal predictive value of inputFor interconnected subsystems, the open-loop prediction output of subsystema 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, Δuat|k=0 for all t≥k. However, the inputs ulto 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 Δu^lt|k of each subsystem l at time k is chosen by the following formula:(3)Δu^lt|k=Δul∗t|k−1,t=k,…,N−1,Δult,t=1,2,…,k−1.Suppose that the inputs from timek onwards do not change, Δuk+t−1|k=01≤t≤p, the corresponding prediction of yk+p is called free prediction, denoted as yfrk+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. Thus, the free prediction of subsystem a can be expressed as follows:(4)yafrk+p|k=∑t=1N+Nc−1StuaΔuak−t+p|k+SN+Ncuauak−N−Nc+p|k+∑l∈Na∑t=1N+Nc−1StualΔu^lk−i+p|kSN+Ncualulk−N−Nc+p|k+∑t=1N−1StfΔfk−t+SNffk−N.The free prediction at timek−1 is(5)yafrk+p|k−1=∑t=1N+Nc−1StuaΔuak−t+p|k−1+SN+Ncuauak−N−Nc+p|k−1+∑l∈Na∑t=1N+Nc−1StualΔu^lk−t+p|k−1+SN+Ncualu^lk−N−Nc+p|k−1+∑t=1N−1StfΔfk−t+p+SNffk−N+p.Suppose thatfk is measurable and will remain unchanged in the future. Then, one has(6)yafrk+p|k−yafrk+p|k−1=Sp+1uaΔuak−1+∑l∈Na∑t=1Nc−1StualΔu^lk−t+p|k−Δu^lk−t+p|k−1.Suppose that the control horizon of neighbor subsysteml is Nc. M coupling inputs from neighbor subsystem l,l∈Na 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(7)∑t=1Nc−1StualΔu^lk−t+p|k−Δu^lk−i+p|k−1=∑t=1pStualΔul∗k−t+p|k−1−Δul∗k−t+p||k−2+S1ualΔulk−1−Δul∗k−1|k−2.As shown in equation (4), the influences of interconnection inputs from neighbors are considered in free prediction yafrk+p|k and then the real prediction for subsystem a can be calculated by(8)yak+p|k=yafrk+p|k+∑t=1pStuaΔuak−t+p|k.Seen from equation (8), at time k, the real output prediction is determined by only the effect of subsystem a 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 p satisfies p≤0.001. By taking S0f=0, it is allowed to take p≤0.001. Then, we have yfrk|k=yk|k≠yk.Feedback correction is added to equation (8) and then the open-loop prediction yaol is obtained. Here, the superscript “ol” denotes open-loop. Based on equation (6), the open-loop prediction of subsystem a can be calculated by(9)yaolk|k=yaolk|k−1+S1uaΔuak−1+∑l∈NaS1ualΔulk−1−Δul∗k−1|k−2.At timek>0, the actual input Δuak−1 of each subsystem a is measurable. Suppose yak is measurable. The open-loop prediction error is obtained by(10)ϵak=yak−yaolk|k,=yak−yaolk|k−1−S1uaΔuk−1,−∑l∈NaS1ualΔulk−1,−Δul∗k−1|k−2.Letεik=Δuik−Δu^ik, then ∑l∈NaS1ualΔulk−1−Δul∗k−1|k−2=∑l∈NaS1ualεlk−1. Let ei∈Ei=∑i∈NaSiualεl, and ei can be considered as a bounded disturbance. The feedback correction ϵak is added into the open-loop prediction equation to compensate the prediction errors. Let N′=N+Nc, the submodel of subsystem a can be written in vector form as follows:(11)YN′,ak|k=AYN′,aolk|k−1+S1uaS2ua⋮SN′uaΔuak−1+∑l∈NaStual×ΔU^l+S1fS2f⋮SN′fΔfak+eak+1|keak+2|k⋮eak+N′|k,A=0I0⋯0⋮⋱⋱⋱⋮0⋯0I00⋯00I0⋯00I,Stual=S1ualS2ual⋯SNc+1ualS2ualS3ual⋯SNc+2ual⋮⋮⋱⋮SN′ualSN′+1ual⋯SN′+Ncual,ΔU^l=Δu^lk−1|k−Δu^lk−1|k−1Δu^lk|k−Δu^lk||k−1⋮Δu^lk+Nc−1|k−Δu^lk+Nc−1|k−1.Notice that the assumptionyaolk+N′|k−1=yaolk+N′−1|k−1 is utilized in deducing equation (11).For a stable process, the feedback correction is time-invariant, that is,eik+t|k=ϵik, for t≥0. Then, the open-loop steady-state prediction of each subsystem can be expressed as(12)ya,ssol=yaolk+N′+t|k,t≥0.When replacingk 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. Thus, the open-loop prediction equation can predict the dynamics of each subsystem effectively.


## 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 ua∈ℝna as input vector, i.e., u=u1,…,uM and ∑a=1Mnua=nu. Similarly, it has ya∈ℝny as output vector, i.e., y=y1,…,uM, and ∑a=1Mnya=ny. In the decomposition, step response coefficient matrices Su11∈ℝn1×n1,…,SuMM∈ℝnM×nM of M subsystems are diagonal blocks of Su. While the nondiagonal blocks of Su (i.e., Sual,a≠l) represent the interconnection influences among subsystems. Sual≠0 means that the input of subsystem l affects the dynamics of subsystem a. Therefore, the FSR model of each subsystem a∈M is expressed as follows:(2)∇yak=∑t=1N−1StuaΔuak−t+SNuauak−N+∑l∈Na∑t=1N−1StualΔulk−t+SNualulk−N′+∑t=1N−1StfΔfk−t+SNffk−N,where ya∈ℝnya,ua∈ℝnua are the output and input vectors of subsystem a, respectively. N denotes the modeling horizon for the large-scale system. Nc denotes the control horizon. l∈Na denotes the neighbors of subsystem a, which excludes a, namely, Na=1,2,…,a−1,a+1,M denoting the collection of the subsystem a’s dynamic neighbor subsystems. Sual represents the interconnection influence on subsystem a from neighbor subsystem l. f∈ℝnf is the disturbance variable.


## 2.2. Distributed Open-Loop Prediction

A large-scale system is decomposed intoM dynamically coupled subsystems, and the distributed control is established based on the subsystems. The following symbols are used to distinguish different prediction information:uk+t|k is the real prediction of inputu^k+t|k is the assumed predictive value of inputu∗k+t|k is the optimal predictive value of inputFor interconnected subsystems, the open-loop prediction output of subsystema 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, Δuat|k=0 for all t≥k. However, the inputs ulto 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 Δu^lt|k of each subsystem l at time k is chosen by the following formula:(3)Δu^lt|k=Δul∗t|k−1,t=k,…,N−1,Δult,t=1,2,…,k−1.Suppose that the inputs from timek onwards do not change, Δuk+t−1|k=01≤t≤p, the corresponding prediction of yk+p is called free prediction, denoted as yfrk+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. Thus, the free prediction of subsystem a can be expressed as follows:(4)yafrk+p|k=∑t=1N+Nc−1StuaΔuak−t+p|k+SN+Ncuauak−N−Nc+p|k+∑l∈Na∑t=1N+Nc−1StualΔu^lk−i+p|kSN+Ncualulk−N−Nc+p|k+∑t=1N−1StfΔfk−t+SNffk−N.The free prediction at timek−1 is(5)yafrk+p|k−1=∑t=1N+Nc−1StuaΔuak−t+p|k−1+SN+Ncuauak−N−Nc+p|k−1+∑l∈Na∑t=1N+Nc−1StualΔu^lk−t+p|k−1+SN+Ncualu^lk−N−Nc+p|k−1+∑t=1N−1StfΔfk−t+p+SNffk−N+p.Suppose thatfk is measurable and will remain unchanged in the future. Then, one has(6)yafrk+p|k−yafrk+p|k−1=Sp+1uaΔuak−1+∑l∈Na∑t=1Nc−1StualΔu^lk−t+p|k−Δu^lk−t+p|k−1.Suppose that the control horizon of neighbor subsysteml is Nc. M coupling inputs from neighbor subsystem l,l∈Na 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(7)∑t=1Nc−1StualΔu^lk−t+p|k−Δu^lk−i+p|k−1=∑t=1pStualΔul∗k−t+p|k−1−Δul∗k−t+p||k−2+S1ualΔulk−1−Δul∗k−1|k−2.As shown in equation (4), the influences of interconnection inputs from neighbors are considered in free prediction yafrk+p|k and then the real prediction for subsystem a can be calculated by(8)yak+p|k=yafrk+p|k+∑t=1pStuaΔuak−t+p|k.Seen from equation (8), at time k, the real output prediction is determined by only the effect of subsystem a 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 p satisfies p≤0.001. By taking S0f=0, it is allowed to take p≤0.001. Then, we have yfrk|k=yk|k≠yk.Feedback correction is added to equation (8) and then the open-loop prediction yaol is obtained. Here, the superscript “ol” denotes open-loop. Based on equation (6), the open-loop prediction of subsystem a can be calculated by(9)yaolk|k=yaolk|k−1+S1uaΔuak−1+∑l∈NaS1ualΔulk−1−Δul∗k−1|k−2.At timek>0, the actual input Δuak−1 of each subsystem a is measurable. Suppose yak is measurable. The open-loop prediction error is obtained by(10)ϵak=yak−yaolk|k,=yak−yaolk|k−1−S1uaΔuk−1,−∑l∈NaS1ualΔulk−1,−Δul∗k−1|k−2.Letεik=Δuik−Δu^ik, then ∑l∈NaS1ualΔulk−1−Δul∗k−1|k−2=∑l∈NaS1ualεlk−1. Let ei∈Ei=∑i∈NaSiualεl, and ei can be considered as a bounded disturbance. The feedback correction ϵak is added into the open-loop prediction equation to compensate the prediction errors. Let N′=N+Nc, the submodel of subsystem a can be written in vector form as follows:(11)YN′,ak|k=AYN′,aolk|k−1+S1uaS2ua⋮SN′uaΔuak−1+∑l∈NaStual×ΔU^l+S1fS2f⋮SN′fΔfak+eak+1|keak+2|k⋮eak+N′|k,A=0I0⋯0⋮⋱⋱⋱⋮0⋯0I00⋯00I0⋯00I,Stual=S1ualS2ual⋯SNc+1ualS2ualS3ual⋯SNc+2ual⋮⋮⋱⋮SN′ualSN′+1ual⋯SN′+Ncual,ΔU^l=Δu^lk−1|k−Δu^lk−1|k−1Δu^lk|k−Δu^lk||k−1⋮Δu^lk+Nc−1|k−Δu^lk+Nc−1|k−1.Notice that the assumptionyaolk+N′|k−1=yaolk+N′−1|k−1 is utilized in deducing equation (11).For a stable process, the feedback correction is time-invariant, that is,eik+t|k=ϵik, for t≥0. Then, the open-loop steady-state prediction of each subsystem can be expressed as(12)ya,ssol=yaolk+N′+t|k,t≥0.When replacingk 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. 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 a can be calculated by(13)ya,ssk=SN′uaδua,ssk+ya,ssolk,where δua,ss denotes the increment variable of the steady-state target of u for subsystem a, δua,ssk=ua,ssk−uak−1. ya,ss denotes the steady-state target value of controlled variable y for subsystem a.


### 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:(14)u¯a≤ua,ssk≤u¯a,k≥0,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Δuak+j|k≤Δu¯0≤j≤M−1 and then δua,ss should further satisfy(15)δua,ssk≤MΔu¯a,k≥0,where δussk=ussk−uk−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 δua,ss is(16)u¯a′≤δua,ssk≤u¯a′k,k≥0,where(17)u¯a′=maxu¯a−uak−1,−MΔu¯a,u¯a′=minu¯a−uak−1,MΔu¯a.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:(18)y¯0,h≤ya,ssk≤y¯0,h,k≥0,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:(19)y¯0≤ya,ssk≤y¯0,k≥0,where y¯0 denotes lower operating limit of y and y¯0 denotes upper operating limit of y.The engineering limits are more stringent than the operating limitsy¯0,h≤y¯0 and y¯0,h≥y¯0.Substituting the steady-state prediction equation (13) into equation (19), we yield(20)y¯0≤SN′uaδua,ssk+ya,ssolk≤y¯0,k≥0.Then, we have(21)y¯k≤SN′uaδua,ssk≤y¯k,where(22)y¯k=minmaxy¯0−yssolk,y¯hk,y¯hk,y¯k=maxminy¯0−yssolk,y¯hk,y¯hk.The ETs ofui,ss and yj,ss are denoted as ui,t and yj,t, respectively. i∈It,j∈Jt, It and Jt represent the corresponding set of all the ETs of ui,ss and yj,ss. Denote the expected allowable range of ui,tk as ui,ss,range. The limits of ET values for ui,ss can be taken as follows:(23)u¯i,ssk=ui,tk+12ui,ss,range−uik−1,i∈It,u¯i,ssk=ui,tk−12ui,ss,range−uik−1,i∈It.Then, foruj,ssk, check the following soft constraint:(24)u¯i,ssk≤δui,ssk≤u¯i,ssk,i∈It.Denote the expected allowable range ofyj,tk as yj,ss,range. The limits of ET values for yj,ss are given by(25)y¯j,ssk=yj,tk+12yj,ss,range−yj,ssolk,j∈Jt,y¯j,ssk=yj,tk−12yj,ss,range−yj,ssolk,j∈Jt.If the external target ofyj,ssk is considered, the following soft constraints need to be further satisfied, that is,(26)SN′uδussk≤y¯k,j∈Jt,SN′uδussk≥y¯k,j∈Jt.Unifying all the ETs into the constraint ofδuss, the corresponding constraint for ui,tk is denoted as(27)δui,ssk=ui,tk−uik−1,i∈It.Foryj,tk, the constraint is(28)SN′,juδussk=yj,tk−yj,ssolk,j∈Jt.


### 3.2. Feasibility Analysis

For each subsystema, 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)):(29)I−ISN′u−SN′uδua,ssk≤u¯′k−u¯′ky¯hk−y¯hk.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 withr rank of priorities, the constraint handled at the r−1 th priority can be denoted as(30)Cr−1δua,ssk≤cr−1k,(31)Ceqr−1δua,ssk≤ceqr−1k.Constraint equations (30) and (31) are treated as hard constraints in the r th priority. The following two cases are considered:(i)
If ther 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:(32a)Cr−1δua,ssk≤cr−1k,(32b)Ceqr−1C˜eqrδua,ssk=ceqr−1kc˜eqrk+εeqrk.(ii)
If ther 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:(33a)Cr−1C˜rδua,ssk≤cr−1kc˜rk+εrk,εrk≥0,(33b)Ceqrδua,ssk=ceqr−1,whereεeqrk and εrk 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, ε¯=y¯h−y¯0, and for the lower bound constraints, ε¯=y¯0−y¯h, and ε=1/2ETrange 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(34)εeqr=εeq+rk−εeq−rk,εeq+rk≥0,εeq−rk≥0,then the slack variablesεeq+rk and εeq−rk can be obtained by solving the following LP:(35)minεeq+r,εeq−rδua,ssk∑τ=1dr+1ε¯eq,τr+1−1×εeq+,τrk+εeq−,τrks.t.32–34,whereτ indicates the τ th element of εeqrk and deqr+1 denotes the dimension of εeqrk.


### 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 ofussk and yssk. Suppose the steady-state change of input variable is bounded by(36)−Uik≤δi,ssk≤Uik.Then, minimizing the increment variableδui,ssk is equivalent to minimizing its constraint bound Uik. Having analyzed the feasibility, all the hard constraints and the relaxed soft constraints are simplified into the following formula:(37)u¯i′k≤δui,ssk≤u′¯ik,i∉Jt,y¯jk≤Sj,N′uδussk≤y¯jk,j∉Jt,Sj,N′uδussk=yj,ssk−yj,ssolk,j∈Jt,δui,ssk=ui,ssk−uik−1,i∈It.By invoking LP, the economic optimization problem of each subsystema is written as(38)minδua,ssk,UikJa=∑i∉Imm∪Ithiδua,ssk+∑i∉ImmhiUik,s.t.19,24,26–28,i∈Imm,where hi is the weighting coefficient.Notice thatδua,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 δua of each subsystem is not necessarily optimal for the whole large-scale system. Define the steady-state target vector as ua,ss=u1,ss,u2,ss,…,uM,ss. According to the economic optimization problem of each subsystem, the global optimization problem can be chosen as(39)minJ∑=min∑a∈MJaδua,ss,where J1δu1,ss,J2δu2,ss,…,JMδua,ss 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 δua,ssp of multiobjective optimization problem (39) is said to be Pareto optimal if and only if there exists no other feasible solution δua,ssf such that Jaδua,ssf≤Jaδua,ssp,∀a∈M, and Jaδua,ssf<Jaδua,ssp for at least one index a∈M.
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 δua,ssL under the condition of Jaδua,ssL, with a∈M. δua,ssL is given by solving the following problem:(40)minδua,ssLJaδua,ss,s.t.13,a∈Mwhere the utopia cost of each subsystem a is denoted as JaL.
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δua,ssC, which can be easily calculated by the following distance problem:(41)minδua,ssJaδua,ss−JaLp,s.t.13,a∈M.
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,J1δu1,ss=J1L.
The compromise solutionδua,ssC calculated by equation (41) is taken as the steady-state target value of subsystem a. Substituting into equation (13), the set-point of ya,ss is obtained. Then, the steady-state value δua,ss=δua,ssC,ya,ss is sent to dynamic control.Figure 1
Schematic of the utopia-tracking approach.


## 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:(14)u¯a≤ua,ssk≤u¯a,k≥0,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Δuak+j|k≤Δu¯0≤j≤M−1 and then δua,ss should further satisfy(15)δua,ssk≤MΔu¯a,k≥0,where δussk=ussk−uk−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 δua,ss is(16)u¯a′≤δua,ssk≤u¯a′k,k≥0,where(17)u¯a′=maxu¯a−uak−1,−MΔu¯a,u¯a′=minu¯a−uak−1,MΔu¯a.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:(18)y¯0,h≤ya,ssk≤y¯0,h,k≥0,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:(19)y¯0≤ya,ssk≤y¯0,k≥0,where y¯0 denotes lower operating limit of y and y¯0 denotes upper operating limit of y.The engineering limits are more stringent than the operating limitsy¯0,h≤y¯0 and y¯0,h≥y¯0.Substituting the steady-state prediction equation (13) into equation (19), we yield(20)y¯0≤SN′uaδua,ssk+ya,ssolk≤y¯0,k≥0.Then, we have(21)y¯k≤SN′uaδua,ssk≤y¯k,where(22)y¯k=minmaxy¯0−yssolk,y¯hk,y¯hk,y¯k=maxminy¯0−yssolk,y¯hk,y¯hk.The ETs ofui,ss and yj,ss are denoted as ui,t and yj,t, respectively. i∈It,j∈Jt, It and Jt represent the corresponding set of all the ETs of ui,ss and yj,ss. Denote the expected allowable range of ui,tk as ui,ss,range. The limits of ET values for ui,ss can be taken as follows:(23)u¯i,ssk=ui,tk+12ui,ss,range−uik−1,i∈It,u¯i,ssk=ui,tk−12ui,ss,range−uik−1,i∈It.Then, foruj,ssk, check the following soft constraint:(24)u¯i,ssk≤δui,ssk≤u¯i,ssk,i∈It.Denote the expected allowable range ofyj,tk as yj,ss,range. The limits of ET values for yj,ss are given by(25)y¯j,ssk=yj,tk+12yj,ss,range−yj,ssolk,j∈Jt,y¯j,ssk=yj,tk−12yj,ss,range−yj,ssolk,j∈Jt.If the external target ofyj,ssk is considered, the following soft constraints need to be further satisfied, that is,(26)SN′uδussk≤y¯k,j∈Jt,SN′uδussk≥y¯k,j∈Jt.Unifying all the ETs into the constraint ofδuss, the corresponding constraint for ui,tk is denoted as(27)δui,ssk=ui,tk−uik−1,i∈It.Foryj,tk, the constraint is(28)SN′,juδussk=yj,tk−yj,ssolk,j∈Jt.


## 3.2. Feasibility Analysis

For each subsystema, 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)):(29)I−ISN′u−SN′uδua,ssk≤u¯′k−u¯′ky¯hk−y¯hk.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 withr rank of priorities, the constraint handled at the r−1 th priority can be denoted as(30)Cr−1δua,ssk≤cr−1k,(31)Ceqr−1δua,ssk≤ceqr−1k.Constraint equations (30) and (31) are treated as hard constraints in the r th priority. The following two cases are considered:(i)
If ther 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:(32a)Cr−1δua,ssk≤cr−1k,(32b)Ceqr−1C˜eqrδua,ssk=ceqr−1kc˜eqrk+εeqrk.(ii)
If ther 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:(33a)Cr−1C˜rδua,ssk≤cr−1kc˜rk+εrk,εrk≥0,(33b)Ceqrδua,ssk=ceqr−1,whereεeqrk and εrk 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, ε¯=y¯h−y¯0, and for the lower bound constraints, ε¯=y¯0−y¯h, and ε=1/2ETrange 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(34)εeqr=εeq+rk−εeq−rk,εeq+rk≥0,εeq−rk≥0,then the slack variablesεeq+rk and εeq−rk can be obtained by solving the following LP:(35)minεeq+r,εeq−rδua,ssk∑τ=1dr+1ε¯eq,τr+1−1×εeq+,τrk+εeq−,τrks.t.32–34,whereτ indicates the τ th element of εeqrk and deqr+1 denotes the dimension of εeqrk.


## 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 ofussk and yssk. Suppose the steady-state change of input variable is bounded by(36)−Uik≤δi,ssk≤Uik.Then, minimizing the increment variableδui,ssk is equivalent to minimizing its constraint bound Uik. Having analyzed the feasibility, all the hard constraints and the relaxed soft constraints are simplified into the following formula:(37)u¯i′k≤δui,ssk≤u′¯ik,i∉Jt,y¯jk≤Sj,N′uδussk≤y¯jk,j∉Jt,Sj,N′uδussk=yj,ssk−yj,ssolk,j∈Jt,δui,ssk=ui,ssk−uik−1,i∈It.By invoking LP, the economic optimization problem of each subsystema is written as(38)minδua,ssk,UikJa=∑i∉Imm∪Ithiδua,ssk+∑i∉ImmhiUik,s.t.19,24,26–28,i∈Imm,where hi is the weighting coefficient.Notice thatδua,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 δua of each subsystem is not necessarily optimal for the whole large-scale system. Define the steady-state target vector as ua,ss=u1,ss,u2,ss,…,uM,ss. According to the economic optimization problem of each subsystem, the global optimization problem can be chosen as(39)minJ∑=min∑a∈MJaδua,ss,where J1δu1,ss,J2δu2,ss,…,JMδua,ss 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 δua,ssp of multiobjective optimization problem (39) is said to be Pareto optimal if and only if there exists no other feasible solution δua,ssf such that Jaδua,ssf≤Jaδua,ssp,∀a∈M, and Jaδua,ssf<Jaδua,ssp for at least one index a∈M.
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 δua,ssL under the condition of Jaδua,ssL, with a∈M. δua,ssL is given by solving the following problem:(40)minδua,ssLJaδua,ss,s.t.13,a∈Mwhere the utopia cost of each subsystem a is denoted as JaL.
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δua,ssC, which can be easily calculated by the following distance problem:(41)minδua,ssJaδua,ss−JaLp,s.t.13,a∈M.
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,J1δu1,ss=J1L.
The compromise solutionδua,ssC calculated by equation (41) is taken as the steady-state target value of subsystem a. Substituting into equation (13), the set-point of ya,ss is obtained. Then, the steady-state value δua,ss=δua,ssC,ya,ss is sent to dynamic control.Figure 1
Schematic of the utopia-tracking approach.


## 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Δua,ssC. For each subsystem a, Ya,Nolk|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:(42)Ya,Nk|k=Ya,Nolk|k+SaΔu˜ak|k,a∈Mwhere(43)Ya,Nk|k=yak+1|kyak+2|k⋮yak+N|k,Δu˜ak|k=Δuak|kΔuak+1|k⋮Δuak+Nc−1|k,Sa=S1ua0⋯0S2uaS1ua⋱⋮⋮⋱⋱0SNcua⋯S2uaS1ua⋮⋮⋮SNua⋯SN−Nc+2uaSN−Nc+1ua,Sa is called the dynamic matrix. N is the prediction horizon, and Nc≤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δua,ss, the cost function is chosen as(44)Jak=∑i=1Nyak+i|k−ya,sskQik2+∑j=0Nc−1Δuak+j|kΛ2,where the weighting matrix Λ is given off-line. Qik=diagqi1k2,qi2k2,…,qi,nyk2, updated online.Besides, an additional constraint to limitδua,ssCk is added to the dynamic control optimization problem:(45)ℒΔu˜ak|k=δua,ssCk,where ℒ=II⋯I. In addition, the inequality constraints including manipulated variable magnitude, manipulated variable increment, and controlled variable magnitude constraints are considered:(46)Δuak+j|k≤Δu¯a,0≤j≤M−1,u¯≤uak−1+∑l=0jΔuak+l|k≤u¯,0≤j≤M−1,y¯a,0′≤yaolk+i|k+SiΔu¯ak|k≤y¯0′k,1≤i≤N,where Si denotes the i th row of dynamic matrix S:(47)y¯a,0′k=maxy¯a,0,ya,ssk,y¯a,0′k=miny¯a,0,ya,ssk.At timek, the following optimization problem of each subsystem is firstly solved:(48)minΔu˜ak|kJak,s.t.45,46,a∈M.If the optimization is infeasible, slack variablesε¯dck and ε¯dck are added into the cost function and then the optimization problem is modified to(49)Ja′k=∑i=1Nyak+i|k−ya,sskQik2+∑j=0M−1Δuak+j|kΛ2+ε¯dckω¯2+ε¯dckω¯2,where ε¯dck and ε¯dck are slack variables. The related constraints on slack variables are as follows:(50)ε¯dck≤y¯a,0′k−y¯a,0,h,1≤i≤N,(51)ε¯dck≤y¯a,0,h−y¯a,0′k,1≤i≤N.By adding the constraints equations (50) and (51), the optimization problem with slack variables is expressed as follows:(52)minε¯dck,ε¯dck,Δu˜ak|kJa′k,s.t.45,46,50,51.Equations (48) and (52) are standard quadratic problem (QP). By solving equations (48) or (52), Δu˜a∗k|k is obtained. Then, uak|k=uak|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. 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 Figure2, 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.Figure 2
Shell heavy oil fractionator.A three-input-three-output model was presented in [31] for the heavy oil fractionator:(53)y1y2y3=Gusu1u2u3+Gdsd1d2,where(54)Gus=4.05e−27s50s+11.77e−28s60s+15.88e−27s50s+15.39e−18s50s+15.72e−14s60s+16.90e−15s40s+14.38e−20s33s+14.42e−22s44s+17.2019s+1,Gds=1.20e−27s45s+11.44e−27s40s+11.52e−15s25s+11.83e−15s20s+11.1427s+11.2632s+1,where the manipulated variables are the top extraction u1, the side extraction u2, and the enthalpy controller for the bottom of the column u3. d1andd2 are the middle reflux and the top reflux, respectively. The controlled variables include top product concentration y1, side product concentration y2, and bottom reflux temperature y3. All the variables are normalized, and the unit of time is minute.As seen from Figure2, by applying the interconnection decomposition, the transfer function Gus in equation (53) can be decomposed into the following three subsystems:Subsystem 1:4.05e−27s/50s+1Subsystem 2:5.72e−14s/60s+1Subsystem 3:7.20/19s+1In 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 by(55)yak=∑i=1N−1SiuaΔuak−i+SNuauak−N+∑l≠a∑iN−1SiualΔulk−i+SNualulk−N+∑i=1N−1SadΔdk−i+SNddk−N,where a=1,2,3 denotes the number of the subsystem. Siua and Siual are the step response coefficient matrices on manipulated variables for local subsystem a and neighbor subsystem l, respectively. Denote SNua and SNual as the steady-state step response coefficient matrices.Based on the open-loop prediction of each subsystema, the corresponding steady-state prediction is obtained by(56)ya,ssk=SNuaδua,ssk+ya,ssolk,where δua,ss denotes the steady-state increment of manipulated variable.The constraints on manipulated variables and controlled variables are specified byu¯a=−0.5,−0.2,−0.35; u¯a=0.5,0.2,0.35; Δu¯a=δu¯a,ss=0.1,0.1,0.1; y¯j,h=−0.4,−0.4,−0.4; y¯j,h=0.4,0.4,0.4; y¯j,0=−0.5; y¯j,0=0.5,j=1,2,3. The external targets of y1,ss,y2,ss,andu3,ss are 0.5. By invoking LP, the steady-state target value of each subsystem a is solved by the following optimization problem:(57)minδua,ss,U2kJa=∑i=1,3hiδui,ssk+h2U2k.Due to the conflicts of the optimization goals of different subsystems, the steady-state target value of each subsystemδua,ss 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 time55≤k≤90, 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 Λ=diag1;4;1, Q=diag4.0;0.67;1.0. Take u1−1=u2−1=u3−1=0, y10=y20=y30=0. The closed-loop system responses of three subsystems are shown in Figures 3–5, respectively.Figure 3
Closed-loop responses of subsystem 1.Figure 4
Closed-loop responses of subsystem 2.Figure 5
Closed-loop responses of subsystem 3.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 time55≤k≤90, the value of disturbance is [0.15; 0.2]. For the centralized algorithm, by tuning parameters, N=8,Nc=5, 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.Figure 6
Tracking performance of the two algorithms.Table 1
Performance comparison.
AlgorithmTrack.-dev.Time-consumingCen. DDMC0.381124.0034Dis. DDMC0.35426.572


## 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.



---

*Source: 4650342-2022-05-30.xml*