Hongguang Pan ,Weimin Zhong *,Zaiying Wang

1 School of Electrical and Control Engineering,Xi'an University of Science and Technology,Xi'an 710071,China

2 Key Laboratory of Advanced Control and Optimization for Chemical Processes,Ministry of Education,East China University of Science and Technology,Shanghai 200237,China

1.Introduction

Model predictive control(MPC)refers to a class of computer control algorithms thatutilize an explicitprocess modelto predictthe future response of a plant.Due to the advantage of controlling the constrained multivariable processes,MPC has been widely and successfully implemented in the process industries over recent years[1-3].

In lots oflarge scale enterprises,a universalhierarchy shown in Fig.1 for optimization and control is generally adopted.With this hierarchical framework,in time scales,the re finement is successive from top to bottom;in the spatialscales,itis also varying from plant-wide optimization at the top to regulatory control via single PID loops at the bottom[4,5].Specially,a real time optimization(RTO)is mainly implemented to improve economic pro fit and give optimal set points(optimal inputs and outputs of the system)to the double-layered M-PC,which has a steady-state target calculation(SSTC)layer and a dynamic control layer;and then the SSTC layer can give steady-state targets which are tracked in the dynamic layer.

The RTO level and the double-layered MPC level are linked via the optimal set points,which are solved in the RTO level and directly translated to the double-layered MPC,i.e.,to the SSTC layer.The SSTC layer is mainly implemented to correct the optimal set points from the RTO level[5,6].Because the steady-state targets may be changed by the disturbance or interference caused by environment and operators atany instant,the SSTC is necessarily inserted on top of dynamic control layer[7-9].Generally,the SSTC layer adopts a linear steady-state model,which is a steady-state version of a dynamic model used in the dynamic control layer.From Fig.1,we can get that the optimization problem of the SSTC layer is calculated at the same frequency as the dynamic control layer.

The steady-state targets can be calculated in a single or in several optimization problems,named the single priority rank method or the multi-priority rank method,respectively.The single priority rank method may be classified into three types,i.e.,the “self-optimizing”,the target tracking and the weighting method according to various objectives of the cost function[9-11].The multi-priority rank method adopted in industrialpractices can distinguish the setpoints in different importance levels[5,11,12].

The stability of the double-layered MPC has always been a hard topic since the double-layered MPC is adopted.However,there have been lots of systematical results on the stability of the MPC when a synthesis approach is adopted[13],and these results can be borrowed to do some research on the stability of the doublelayered MPC in some extent.Generally,a region of attraction(ROA)of the dynamic controller is inconsistent with the feasible domain of the SSTC layer,and the inconsistency between the doublelayered MPC will result in the in feasibility of the dynamic control layer during the control process.Another advantage of adopting the synthesis approach in dynamic control layer is that the ROA can be conveniently calculated.

Fig.1.Control structure with RTO and double layered MPC.

Based on the above consideration,the synthesis approach is adopted in this paper and an on-line constraintsoftening strategy is presented to guarantee the feasibility of the dynamic control layer.Specially,a series ofsoftened constraints are given via relaxing the softconstraints several times,and a series of ROA centering at the steady-state targets is calculated at each instant;then,the ROA is chosen on-line considering the following two conditions:1)the ROA should contain the current state at each instant and 2)the relaxation should be as small as possible.When the ROA is fixed,the correspondingly softened constraint is adopted in dynamic control layer and the feasibility can be guaranteed.

The paper is organized as follows.Section 2 gives preliminary results of the double layered MPC and analyzes the reason of the inconsistency between the SSTC and the dynamic control;in Section 3,the calculation process of the approximate ROA and the presented algorithm are introduced;Section 4 gives two simulation examples and Section 5 concludes the paper.

Notation:Rndenotes the n-dimensional Euclidean space,x*the optimal value of x;y(x,u)denotes the system output y∈Rny(state x∈Rnx,input u∈Rnu),ys(xs,us)the steady-state output(state,input),yt(ut)the set points of the output(input),yˉs(xˉs,uˉs,yˉt,uˉt)the upper bound of ys(xs,us,yt,ut);y(k+i|k)(x(k+i|k),u(k+i|k))denotes the future value of u(k+i)(x(k+i),x(k+i))predicted at instant k;Qs,Rs,Q,R,P denote the weighting matrices,which are positive definite;⊕denotes the Minkowski sum,i.e.,A⊕B={a+b|?a∈A,b∈B};Indenotes the n-dimensional identity matrix,N the control horizon;||x||denotes the 2-norm of x and||x||Q2denotes xTQx.

2.Preliminary Results of the Double-layered MPC

In this paper,a numerical algorithm for the subspace state space system identification(N4SID)method ofsubspace identification method(SIM)is adopted to get the state space model,which is adopted in the double-layered MPC.The reason is that,compared with the traditional approaches,SIM does not need nonlinear search(thence avoiding local minima and convergence problems)and canonical parameterization,which is intrinsically suitable for the identification of multivariable system[14].Because the N4SIDmethod is introduced in many literatures[14,15],it is unnecessary to give details here.

The dynamic state space model{A;B;C}with nuinputs,nxstates and nyoutputs is obtained via N4SID:

which subjects to the following constraints on the outputs,states and inputs:

where Y,X and U denote the admissible set,and are assumed to be nonempty compact convex polyhedrons containing the origin in their interior,f∈ Rnu×nu.The above model is assumed to satisfy the following assumptions:

Assumption 1.The pair(A,B)is stabilizable.

Correspondingly,the steady-state model is as follows:

and the steady-state states and inputs are in the following sets:

Generally,Ys=Y，Xs=X，Us=Y.

2.1.SSTC

In practical industrial applications,almost all MPC software contains SSTC operation before implementing dynamic control.In this paper,we choose the targets tracking cost function for the SSTC.For the set points ytand ut,any ys(k),xs(k)and us(k)satisfying Eq.(6)are on-line calculated.The SSTC optimization problem is

Generally,the solution of problems(7a)and(7b)is unique by minimizing the deviation between the set points and the steady-state targets in least squares sense[16].If(yt,ut)are not reachable,which means that they cannot be reached by any controller,the solution will converge to a solution(ys(k),us(k))which has the minimal distance from(yt,ut),and Js(k)is not equal to zero.With the changing of the set points,the steady-state targets can locate at any point in the admissible set.

In this paper,the setpoints foroutput(yt)and input(ut)are selected(or calculated)through the RTO layer(see Fig.1 in Section 1).The set points remain unchanged when the RTO layer does not give new set points to the double layered MPC.Considering this papermainly focuses on double layered MPC and the feasibility guaranteeing of dynamic layer,the values of set points are given artificially in this paper.

2.2.Dynamic control

In the dynamic control layer,the model predictive controller is designed by the synthesis approach,and the following assumption should be considered.

Assumption 2.

1)K ∈ Rnu×nxis a local stabilizing control gain such that(A+BK)is Schur;

2)P ∈ Rnx×nx,Q ∈ Rnx×nx,R ∈ Rnu×nuare positive definite matrices,where P ∈ Rnx×nxis such that(A+BK)TP(A+BK)-P=-(Q+KTRK);

3)For tracking for systems(1a)and(1b),Xfis a terminal invariant set subject to Eqs.(3)and(4)and K.

When the steady-state targets(ys(k),xs(k),us(k))(where ys(k)=Cxs(k))are calculated in the SSTC layer,the dynamic control optimization problem is shown as follows:

where i=0,1,…,N-1,u′=[u(k|k),u(k+1|k),…,u(k+N-1|k)];the current state can be chosen as(k),which is estimated though a Kalman filter.The Kalman filer is brie fly introduced in Appendix.

In optimization problems(8a)-(8f),all of the states are used in control.Indeed,although the cost functionin Eqs.(7a)and(7b)does not include the steady-state,the steady-state state xs(k)can be calculated as a decision variable(ys(k),xs(k),us(k)are the three decision variables in Eqs.(7a)and(7b)).

Generally,the synthesis method includes the following“three ingredients”,i.e.,the terminal cost function||x(k+N|k)-xs(k)||P2,the terminal invariant set Xfand the local control gain K.The definition of ROA is given here according to[17,18],i.e.,the ROA is a set of states,which can be steered to the terminal invariant set in N steps or less,while assuring that no input and state constraints are violated.According to this definition,the ROA crucially re flects the control ability.

Based on the brief introduction on the double layered MPC,we can find that the SSTC layer only requires the steady-state targets xs(k)to satisfy the constraints(i.e.,Eq.(7b)),while the dynamic controllayerrequires the predicted values of state in the predictive horizon to satisfy the constraints(i.e.,Eq.(8e)).According to the definition of ROA,when we adopt the synthesis method,if and only if the current state x(k|k)is contained in the ROA,the optimization problems(8a)-(8f)in the dynamic control layer are feasible;else,Eqs.(8a)-(8f)are infeasible.In the actual process,the state constraints of the controller may be modulated with the relaxed or tightened soft constraints,then,the scale of the ROA will change correspondingly.Based on the above analysis,it is found that,during the actual process,there is no specific measure which can guarantee that,the ROA can contain the current state x(k|k)at any instant.Once the current state is out of the ROA,the optimization problems(8a)-(8f)are infeasible.

3.Presented Algorithm for the Double-layered MPC

In this section,we give a solution,which can guarantee the feasibility of the dynamic control optimization problem.At each k,when the steady-state targets are calculated in the SSTC layer,they are passed to the dynamic control layer,and then,a series of ROA can be calculated centering at xs(k)through softening the soft constraints.These measures can guarantee that,at least one ROA can contain the current state.Next,the details are introduced.

Firstly,the terminal invariant set is calculated,and then the ROA are introduced.Here,we adopt the linear quadratic regulator(LQR)and Riccati equation to calculate the local controller gain K,which can be applied to transform an original system into an autonomous system.The autonomous system is shown as

Then,the method proposed in[19,20]is adopted in this paper to get the maximal terminal invariant set Xf（xs（k））,which centers at xs(k).Moreover,based on the maximal terminal invariant sets shown in Eqs.(10a)and(10b),the ROA O（xs（k））can be calculated[20,21]:

where Ol（xs（k））(l=1,2,…,N)is the ROA for control horizon i.In Eqs.(10a)and(10b),O（xs（k））=ON（xs（k））.

3.1.Choosing a “suitable”ROA containing current state on-line

Based on the analysis in last section,through enlarging the size of ROA,the current state can be contained,and the feasibility of the dynamic controller can be guaranteed.

Here,letus soften the softconstraints L times,and Δxˉis the softened amplitude.After the softening,the enlarged state bounds are xˉ0,xˉ1,…,xˉl,…,xˉL(where xˉ0=xˉis the original bound),and Δxˉ=xˉ1-xˉ0= … =xˉl-xˉ(l-1)=…=xˉL-xˉ(L-1).Correspondingly,the admissible sets associated with the enlarged state bounds are rewritten as X0，X1，…，Xl，…，XL.According to this measure,afterobtaining the new bounds,a series of maximal terminal invariant sets,i.e.,（xs（k））centering at xs(k)can be obtained(note that,Xlf（xs（k））is associated with xˉl,l=0,1,…,L).Further,the ROA Ol（xs（k））can be calculated through Eqs.(10a)and(10b).

Assumption3.By choosing the suitable N,L andΔxˉ,there exists a maximal ROA OL（xs（k））satisfying

where X0=X.

Note that,the “suitable”,in the assumption,means that Ol（xs（k））contains(xs(k))with the smallest softening time l.

The above assumption is rational.The rationality is:in order to satisfy Eq.(11),the ROA can be enlarged enough by relaxing the soft constraints(and increasing the control horizon N can also be used to achieve this intention,but it will increase the computation burden obviously).Hence,the trade-off between the computational burden and the relaxed amplitude(associated with L and Δxˉ)should be carefully considered.

Based on the analysis above,the improved optimization problem in the dynamic layer is rewritten as

where i=0,…,N-1.Through solving this problem,the control inputs are obtained and injected to drive the system to the steady-state targets.

3.2.On-line recursive feasible algorithm

Algorithm 1.off-line stage:

(i)select the predictive horizon N,the weighting matrices Qs,Rs,Q,R,the initial Kalman filter covariance H(0|0)and the initial state x(0|0);

(ii)select L and Δxˉ,and get the admissible sets X0，X1，…，Xl，…，XL;

(iii)obtain the local control gain K and the Lyapunov matrix P;

on-line stage:at each k＞0,implement the following steps,

(i)in SSTC layer,obtain the unique xs(k)and us(k)through solving Eqs.(7a)and(7b),and judge whether xs(k)changes or not;if it changes,then continue,else,go to(iii);

(ii)in SSTC layer,calculate Xf（xs（k））and O0（xs（k）），O1（xs（k）），…，Ol（xs（k）），…，OL（xs（k））;

(iii)in SSTC layer,select the “suitable”O(jiān)l（xs（k））with the smallest l to satisfy x（k）∈ Ol（xs（k））,then transmit xs(k),us(k),Xland Xlf（xs（k））to the dynamic control layer;

(iv)in dynamic control layer,solve Eqs.(12a)-(12g)to obtain the optimal control inputs u*(k|k),…,u*(k+N-1|k),and inject the first control input u*(k|k)into the actual plant.

In this on-line algorithm,although choosing a suitable ROA is added to the SSTC layer,the basic structures of the double-layered MPC are not changed,hence,all advantages of the double-layered MPC are remained.

4.Numerical Example and Simulation

In this paper,we show the effectiveness of the presented algorithm through the following two examples.In the first numerical example,the relationship between the ROA and the current state is demonstrated in two figures.In the second simulation example,we mainly give the control results.

4.1.Numerical example 1

Adopt the following LTI system to show the effectiveness of the proposed algorithm,

whereand the constraint bounds are yˉ=4 and uˉ=0.6,respectively.

According to Algorithm 1,the parameters are given as follows:

(i)choose the predictive horizon N=7,the weighting matrix Qs=I2,Rs=I2,Q=I2,R=I2,and initial state x(0|0)=[-8,8]T.Note that,the controlled system is the same as the state space model in this numerical example,hence,the current state x(k|k)is chosen as x(k)and the Kalman filter is unnecessary;

(ii)choose L=3 and Δyˉ=CΔxˉ=[3,3]T,i.e.,CΔxˉ1=C xˉ0+[3,3]T,CΔxˉ2=CΔxˉ1+[3,3]Tand CΔxˉ3=CΔxˉ2+[3,3]T;

(iii)the localcontroller gain K and the weighting matrix P are calculated by the LQR and the Riccati equation solution as:

We verify the effectiveness through a two-stage simulation.The desired targets are yt=xt=[4,2]T,ut=[0,0]T(k∈[0,30]),yt=xt=[2.5,2]T,ut=[0,0]T(k∈[31,60]),respectively.Note that,only the desired targets are changed in each stage.

In the SSTC layer,the steady-state targets calculated ateach instantfor stages 1 and 2 are ys=xs=[4.00,1.21]T,us=[0.29,-0.33]T(k∈[0,30]),and ys=xs=[-4.00,-1.21]T,us=[-0.29,0.33]T(k∈[31,60]),respectively.We can find that not all of the desired targets can be reached in steady-state.

In the dynamic controllayer,the state evolutions overthe two stages are shown in Figs.2,3 and Table 1.

In Table 1,the ROA containing the current state at each instant k are listed,and we can find that O0（xs（k））at instants k=1-34 do not contain the current state x(k|k).If the presented method is not adopted,the dynamic control optimization problem is infeasible at the above instants.Further,Figs.2 and 3 can also give the same conclusion.

4.2.Simulation example 2

Consider the following transfer matrix of a heavy oil fractionator[22]:

Fig.2.The state evolution of example 1 in stage 1.

Fig.3.The state evolution of example 1 in stage 2.

Table 1 The chosen results of the ROA in example 1

The three inputs u1,u2,and u3are the product draw rate from the top of the column,the product draw rate from the side of the column,and the re flux heat duty for the bottom of the column,respectively.The three outputs y1,y2,and y3are the draw composition from the top of the column,the draw composition from the side of the column and the re flux temperature at the bottom of the column,respectively.The inputs and outputs are constrained between-0.5 and 0.5,i.e.,uˉ=yˉ=[0.5,0.5,0.5]T.

Here,we adopt the N4SID method to get the state space model,and the identification process is introduced brie fly.The open loop inputs and outputs data are obtained by exciting the system(13)with 3 general binary noise(GBN)signals of magnitude 0.5 as the inputs.2000 sample points are measured as the identification data with the sampling period 4 min.The system order is set as n=6.Through N4SID command in Matlab R2013a,the state space matrices{A,B,C}are identified,and matrix D is set to zero.{A,B,C}are as follows:

Similarly,the parameters for Algorithm 1 are given as follows:

(i)choose the predictive horizon N=5,the weighting matrix Qs=I3,Rs=I3,Q=diag[1,1,2],R=I3,the initial state y(0|0)=Cx(0|0)=[0 0 0]Tand the initial Kalman filter covariance H(0|0)=100,000I6;

(ii)choose L=3 and Δyˉ=CΔxˉ=[0.1,0.1,0.1]T,i.e.,CΔxˉ1=C xˉ0+[0.1,0.1,0.1]T,CΔxˉ2=CΔxˉ1+[0.1,0.1,0.1]Tand CΔxˉ3=CΔxˉ2+[0.1,0.1,0.1]T;

(iii)K and P are

In this example,we also give a two-stage simulation,and the desired targets in the first and second stages are chosen as yt=[0.4,0.3,0.3]T,ut=[0.3,-0.25,-0.25]T(k∈[0,100]),and yt=[0.5,-0.5,0.4]T,ut=[0.2,-0.2,-0.2]T(k∈[101,600]),respectively;then,the steady-state targets are calculated via Eqs.(7a)and(7b)in the SSTC layer,and the results are ys=[0.48,0.36,0.16]T,us=[0.25,-0.01,-0.07]T(k∈[0,100]),and ys=[0.50,-0.36,0.18]T,us=[-0.08,-0.24,0.21]T(k∈[101,600]),respectively.

In the dynamic control layer,the chosen results of the ROA are shown in Table 2,and the control outputs and inputs corresponding to the steady-state targets given through problem are shown in Fig.4.During the control process(k=101),another group of desired targets(yt=[0.5,-0.5,0.4]Tand ut=[0.2,-0.2,0.2]T)are given by the RTO layer,then,the new group of steady-state targets(ys=[0.48,0.36,0.16]Tand us=[0.25,-0.01,-0.07]T)are calculated through Eqs.(7a)and(7b),and tracked by the dynamic layer.In Table 2,wecan find that O0（xs（k））at instants k ∈ [122,149]and[169,195]do not contain the currentstate x(k|k).Ifthe presented method is notadopted,the dynamic control optimization problem is infeasible at the above instants.Further,Fig.4 can also give the same conclusion.

Table 2 The chosen results of the ROA in example 2

Fig.4.The control results of example 2.

In the above two examples,Eqs.(7a)and(7b)are standard quadratic programming problems which are solved adopting the Matlab optimization toolbox;the sets（xs（k））and Ol（xs（k））are computed via Multi-Parametric Toolbox(MPT)in Matlab(see e.g.[23]);and Eqs.(12a)-(12g)are solved through the Linear Matrix Inequality(LMI)toolbox in Matlab.

5.Conclusions

In this paper,we firstly analyze the inconsistency in the doublelayered model predictive control,and obtain the guideline to solve the infeasibility of the dynamic control layer.Further,the on-line constraint softening strategy is given to guarantee the feasibility of the dynamic control layer through on-line choosing a “suitable”region of attraction.This work lays a foundation for solving the stability problem of the double-layered MPC.In the future,the following two works should be continued,1)considering the disturbance and noise,which is more close to the real processes and 2)analyze the stability of the doublelayered MPC based on this work.

Appendix.Brief introduction of Kalman filter

The gain matrix Kkand the current state are updated as follows:

(1)discrete Kalman filter time update equations:

(2)discrete Kalman filter measurement update equations:

where H denotes Kalman filter covariance,W measurement noise covariance,and I unit matrix.When an initial state estimation(0)and an initial covariance H(0|0)are given,the above two steps are repeated to on-line update the filtergain Kkand the current state estimation(k)with each instant.

[1]D.Q.Mayne,Model predictive control:recent developments and future promise,Automatica 50(12)(2014)2967-2986.

[2]R.Amrit,J.B.Rawlings,L.T.Biegler,Optimizing process economics online using model predictive control,Comput.Chem.Eng.58(2013)334-343.

[3]M.L.Darby,M.Nikolaou,MPC:current practice and challenges,Control.Eng.Pract.20(4)(2012)328-342.

[4]M.A.Miiller,D.Angeli,F.Allgower,R.Amrit,J.B.Rawlings,Convergence in economic model predictive control with average constraints,Automatica 50(12)(2014)3100-3111.

[5]M.L.Darby,M.Nikolaou,J.Jones,Doug Nicholson,RTO:an overview and assessment of current practice,J.Process Control 21(6)(2011)874-884.

[6]A.G.Marchetti,A.Ferramosca,A.H.Gonzalez,Steady-state target optimization designs for integrating real-time optimization and model predictive control,J.Process Control 24(1)(2014)129-145.

[7]D.E.Kassmann,T.A.Badgwell,Robert B.Hawkins,Robust steady-state target calculation for model predictive control,AIChE J.46(5)(2000)1007-1024.

[8]A.Nikandrov,C.L.E.Swartz,Sensitivity analysis of LP-MPC cascade control systems,J.Process Control 19(1)(2009)16-24.

[9]C.V.Rao,J.B.Rawlings,Steady states and constraints in model predictive control,AIChE J.45(6)(1999)1266-1278.

[10]J.B.Rawlings,R.Amrit,Optimizing process economic performance using model predictive control,Nonlinear Model Predictive Control,Springer 2009,pp.119-138.

[11]H.G.Pan,H.N.Gao,Y.Sun,Y.Zhang,B.C.Ding,The algorithm and software implementation for double layered model predictive control based on multi priority rank steady state optimization,Acta Automat.Sin.40(3)(2014)405-414.

[12]T.Zou,H.Q.Li,X.X.Zhang,Y.Gu,H.Y.Su,Feasibility and soft constraint of steady state target calculation layer in LP-MPC and QP-MPC cascade control systems,Advanced Control of Industrial Processes(ADCONIP),2011 International Symposium on,IEEE 2011,pp.524-529.

[13]D.Q.Mayne,J.B.Rawlings,C.V.Rao,P.O.M.Scokaert,Constrained model predictive control:stability and optimality,Automatica 36(6)(2000)789-814.

[14]B.Huang,R.Kadali,Dynamic Modeling,Predictive Control and Performance Monitoring,Springer,London,2008.

[15]P.V.Overschee,B.L.R.De Moor,Subspace Identification for Linear Systems:Theory,Implementation,Applications,Springer US,New York,1996.

[16]K.R.Muske,Steady-State Target Optimization in Linear Model Predictive Control,American Control Conference.Proceedings of the 1997,vol.6,IEEE 1997,pp.3597-3601.

[17]D.Limon,T.A.,E.F.Camacho,Enlarging the domain of attraction of mpc controllers,Automatica 41(4)(2005)629-635.

[18]F.Blanchini,S.Miani,Any domain of attraction for a linear constrained system is a tracking domain of attraction,SIAM J.Control.Optim.38(3)(2000)971-994.

[19]B.Pluymers,J.A.Rossiter,J.A.K.Suykens,Bart De Moor,The efficient computation of polyhedral invariant sets for linear systems with polytopic uncertainty,American Control Conference.Proceedings of the 2005,IEEE 2005,pp.804-809.

[20]F.Blanchini,Survey paper:set invariance in control,Automatica 35(11)(1999)1747-1767.

[21]E.C.Kerrigan,Robust Constraint Satisfaction:Invariant Sets and Predictive Control(PhD thesis)University of Cambridge,2001.

[22]D.M.Prett,M.Morari,The Shell Process Control Workshop,Elsevier,2013.

[23]M.Kvasnica,P.Grieder,M.Baotic,M.Morari,Multi-Parametric Toolbox(MPT),Hybrid System:Computation and Control,Springer,Berlin 2004,pp.448-462.