This paper introduces a hierarchical multirate MPC framework for interconnected linear systems, combining a robust centralized high-level controller with local low-level controllers operating at different rates to enhance control performance.
Contribution
It proposes a novel hierarchical multirate MPC scheme that integrates a reduced order centralized controller with local controllers at varying rates for interconnected systems.
Findings
01
Effective coordination of hierarchical controllers demonstrated in simulation.
02
Potential for improved control performance in interconnected systems.
03
Framework adaptable to different system configurations.
Abstract
This paper presents a hierarchical control scheme for interconnected linear systems. At the higher layer of the control structure a robust centralized Model Predictive Control (MPC) algorithm based on a reduced order dynamic model of the overall system optimizes a long-term performance index penalizing the deviation of the state and the control input from their nominal values. At the lower layer local MPC regulators, possibly working at different rates, are designed for the full order models of the subsystems to refine the control action computed at the higher layer. A simulation experiment is presented to describe the implementation aspects and the potentialities of the proposed approach.
Tables1
Table 1. Table 1: On-line computation time comparison
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsAdvanced Control Systems Optimization · Fault Detection and Control Systems · Control Systems and Identification
Full text
A hierarchical multirate MPC scheme for interconnected systems
extended version
M. Farina
X. Zhang
R. Scattolini
Dipartimento di Elettronica, Informazione e Bioingegneria, Politecnico di Milano, Milan 20133, Italy (e-mail: {marcello.farina,xinglong.zhang,riccardo.scattolini}@polimi.it).
Abstract
This paper presents a hierarchical control scheme for interconnected linear systems. At the higher layer of the control structure a robust centralized Model Predictive Control (MPC) algorithm based on a reduced order dynamic model of the overall system optimizes a long-term performance index penalizing the deviation of the state and the control input from their nominal values. At the lower layer local MPC regulators, possibly working at different rates, are designed for the full order models of the subsystems to refine the control action computed at the higher layer. A simulation experiment is presented to describe the implementation aspects and the potentialities of the proposed approach.
Physical and cyber-physical systems are becoming more and more complex, large-scale, and heterogeneous due to the growing opportunities provided by information technology in terms of computing power, transmission of information, and networking capabilities. As a consequence, also the management and control of these systems represent a problem of increasing difficulty and require innovative solutions. A classical approach consists of resorting to decentralized or distributed control structures, see the fundamental book [36] and, in the context of Model Predictive Control (MPC) here considered, the recent review [34] on distributed MPC (DMPC) or the book [28], where the most recent contributions to the design DMPC methods are reported. According to [34], DMPC algorithms can be cooperative, characterized by an intensive transmission load due to the multiple exchange of information among the regulators within one sampling period, or non-cooperative, characterized conservativeness to compensate for the effects of neglected dynamics.
With the aim of reducing the above limitations of DMPC algorithms, in this paper we propose the novel two-layer control scheme shown in Figure 1.
The system under control Σ is assumed to be composed by M interconnected subsystems Σ1,...,ΣM. A reduced order model Σˉ is first computed and a centralized MPC regulator RH working at a slow rate is designed to consider the long-term behavior of the system and to compute the control variables uˉi, i=1,...,M. Then, faster local regulators RLi, i=1,...,M, are designed for each subsystem Σi to compute the control contributions δui compensating for the inaccuracies of the high layer design due to the mismatch between Σ and Σˉ. Notably, the local regulators can be designed and implemented at different rates to cope with subsystems operating at different time scales, as it often happens in many important industrial fields, see the centralized multirate MPC methods reported in [33, 4, 23, 8, 7, 22, 2], or the multirate implementations described in [32, 16, 38].
A similar control structure has been already proposed in [29] where however the framework was different: only independent systems Σi with joint output constraints were considered, no multirate implementations were allowed, and the problem was to coordinate the Σi′s to guarantee an overall output request, so that different technical tools, with respect to the ones here considered, had to be adopted at the two layers of the control structure.
The paper is organized as follows. Section 2 introduces the models considered at the two layers. Section 3 describes the MPC algorithms adopted at the two layers, while Section 4 presents the main feasibility and convergence results as well as a summary of the main implementation aspects. Section 5 describes the simulation example, while in Section 6 some conclusions are drawn. The proofs of the main results are reported in the Appendix.
Notation: for a given set of variables zi∈Rqi,
i=1,2,…,M, we define the vector whose vector-components are zi in the following compact form: (z1,z2,⋯,zM)=[z1Tz2T⋯zMT]T∈Rq, where q=∑i=1Mqi.
We use N+ to denote the set of positive integer numbers.
The symbols ⊕/⊖ denote the Minkowski sum/difference. We denote by ∥⋅∥ the Euclidean norm.
Finally, a ball with radius ρεi and centered at xˉ in the Rdim
space is defined as follows
[TABLE]
2 Models for the two-layer control scheme
In this section we present the model of the overall system and the simplified one used for high-level control.
2.1 Large-scale system model
In line with [25], we assume that the overall system Σ is
composed by M discrete-time, linear, interacting subsystems described by
[TABLE]
i=1,2,…,M, where xi⊆Rni, ui∈Rmi, and yi∈Rpi
are the state, input, and output vectors, while h is the discrete-time index in the basic time scale according to which the models are defined and the low level regulators will be designed.
The interconnections among the subsystems Σi are represented by the coupling input and output vectors si∈Rpsi
and zi∈Rpzi, respectively, where
[TABLE]
with Lii=0, i=1,...,M.
From (1) and (2), the overall dynamical model Σ is
[TABLE]
where x=(x1,…,xM)∈Rn, n=∑i=1Mni,
u=(u1,…,uM)∈Rm, m=∑i=1Mmi, and y=(y1,…,yM)∈Rp, p=∑i=1Mpi. The diagonal blocks of AL
are state transition matrices ALii, whereas
the coupling terms among the Σi′s correspond to the non-diagonal
blocks of AL, i.e., ALij=ELiLijCLzj,
with j=i. The collective input and output matrices are BL=diag(BL11,...,BLMM) and CL=diag(CL11,...,CLMM), respectively.
Concerning systems (1) and (3), the
following standing assumption is introduced:
Assumption 1
The state xi is measurable, for each i=1,…,M;
2. 2.
AL* is Schur stable;*
3. 3.
m=p* and the system matrix
SΣ=[I−ALCL−BL0]
has full rank n+m;*
4. 4.
BL* and CL are full-rank matrices;*
5. 5.
the pair (ALii,BLii) is
reachable, for each i=1,…,M. □
2.2 Reduced order models
For each subsystem Σi,i=1,...,M, we define a reduced order model Σˉi,i=1,...,M, with state
xˉi∈Rnˉi, nˉi≤ni, and input uˉi∈Rmi.
In a collective form, these systems Σˉi
define the overall reduced order model
[TABLE]
where xˉ=(xˉ1,…,xˉM)∈Rnˉ,
nˉ=∑i=1Mnˉi, uˉ=(uˉ1,…,uˉM)∈Rm, and yˉ∈Rp.
The reduced order models Σˉi can be defined according to different criteria. First, it is necessary that the stability properties of system Σ are inherited by Σˉ. Moreover, we assume that for each subsystem i=1,…,M there exists a state projection
βi:Rni→Rnˉi, i=1,...,M, that allows to establish a connection between the states xi(h) of the original models and the states of the reduced models xˉi(h). Collectively, we define β=diag(β1,...,βM). In principle, the ideal case would be to verify that, if uˉ(h)=u(h), then xˉ(h)=βx(h) and yˉ(h)=y(h) for all h≥0 for suitable initial conditions. However, due to model reduction approximations, this ideal assumption must be relaxed; instead, we just ask that xˉ=βx and that yˉ=y in steady-state conditions. We also require that, while matrix BH can be full, the output matrix CH preserves the block-diagonal form of CL, i.e., that CH=diag(CH11,…,CHMM), where CHii∈Rpi×nˉi for all i=1,…,M.
Overall, we require the following standing assumption to be satisfied to guarantee the compatibility of the models used at the two layers.
Assumption 2
AH* is Schur stable;*
2. 2.
βi* is full rank and is such that CHiiβi=CLii, for each i=1,…,M;*
3. 3.
letting G^L(z)=β(zI−AL)−1BL and GH(z)=(zI−AH)−1BH, it holds that GH(1) is full rank and G^L(1)=GH(1). □
An algorithm to compute the projections βi and the matrices of Σˉ is discussed in Appendix 7.1, along the lines of [29].
3 Design of the hierarchical control structure
In this section the regulators at the two layers of the hierarchical control structure are designed for the solution to a tracking control problem, i.e., to drive the output y(h) of the system Σ to the reference value yS, while respecting suitable input constraints.
Thanks to Assumption 1.(3), it is possible to compute the reference pair (xS,uS) corresponding to yS, i.e., such that xS=ALxS+BLuS and CLxS=yS. Correspondingly, we define uˉS=uS as the steady-state input reference for the reduced-order system Σˉ, and the corresponding reference steady-state value as xˉS=GH(1)uˉS=βxS by Assumption 2.(3).
At the same time, we aim to enforce input constraints of type ui(h)∈uˉS,i⊕US,i for all i=1,…,M, where uˉS,i is the i-th vector component of uˉS and US,i, i=1,…,M, are closed and convex sets containing the origin in their interiors. Note that, if the reference uˉS,i changes, also the set US,i may vary to enforce absolute input limitations or saturations. At a collective level, the required constraints are u(h)∈uˉS⊕US, where US=∏i=1MUS,i is a closed and convex set containing the origin in its interior.
3.1 Design of the high level regulator
The high level regulator, designed to work at a low frequency, is based on the reduced order model (4) sampled
with period NL under the assumption that, ∀k∈N,
the uˉi′s are held constant over the interval h∈[kNL,(k+1)NL−1]. Therefore, the sampling time of the high-level model is NL times larger than the basic sampling time, used in the model (1).
Denoting by uˉi[NL](k) the constant
values of uˉi in the long sampling period k and by uˉ[NL](k) the overall
input vector, the reduced order model in the slow timescale is
[TABLE]
where BH[NL]=∑j=0NL−1AHjBH. To enforce the input constraints specified above, we will require that uˉ[NL](k)∈uˉS⊕UˉS, where UˉS=∏i=1MUˉS,i, UˉS,i⊂US,i for each i=1,…,M, and where the properties of sets UˉS,i are specified later in the paper.
In order to feedback a value of xˉ[NL] related to the real
state x of the system, the projected value βixi(kNL) is used, so that the reset
[TABLE]
must be applied for all i=1,…,M. In collective form (6) becomes
[TABLE]
The reset (6) at time k
may force xˉ[NL](k+1) to assume a value different from
the one computed based on the dynamics (5)
and the applied input uˉ[NL](k). This discrepancy, due to the model reduction error and to the actions of the low level controllers,
is accounted for by including in (5) an additive disturbance wˉ(k), i.e.,
[TABLE]
The size of wˉ(k) depends on the action of the low level regulators and its presence requires to resort to a robust MPC method, which is here designed assuming that wˉ(k)∈W, where W is a compact set containing the origin. The characteristics of W will be defined in the following once the low level regulators have been specified (see Section 4).
The robust MPC algorithm is based on the scheme proposed in [27]. To this end, we first need to define the “unperturbed” prediction model
[TABLE]
and the control gain matrix KˉH such that, at the same time
•
FH=AHNL+BH[NL]KˉH is Schur stable.
•
FL[NL]=ALNL+BL[NL]KˉHβ is Schur stable, where BL[NL]=∑j=0NL−1ALjBL.
We define eˉ(k)=xˉ[NL](k)−xˉ[NL],o(k) and we let Z be a robust positively invariant (RPI) set - minimal, if possible - for the autonomous but perturbed system
[TABLE]
The prediction horizon for the high-level MPC consists of NH slow time steps. Denoting by uˉ[NL],o(t:t+NH−1) the sequence uˉ[NL],o(t), …, uˉ[NL],o(t+NH−1), at each slow time-step t the following optimization problem is solved:
[TABLE]
where the input constrained set has been properly tightened in accordance with the used tube-based control approach, and where
[TABLE]
The set XˉF is a positively invariant terminal set for the unperturbed system (9) controlled with the stabilizing control law uˉ[NL],o(k)=KˉHxˉ[NL],o(k), satisfying KˉHXˉF⊆USˉ⊖KˉHZ. In view of this, xˉS⊕XˉF results positively-invariant for (9) controlled with the stabilizing auxiliary control law uˉ[NL],o(k)=uˉS+KˉH(xˉ[NL],o(k)−xˉS).
The positive definite and symmetric weighting matrices QH, RH are free design parameters, while PH is computed as the solution to the Lyapunov equation
[TABLE]
Letting xˉ[NL],o(t∣t),uˉ[NL],o(t:t+NH−1∣t) be the solution to the optimization problem (11), the control action, applied to system Σˉw[NL] at time t, is
[TABLE]
3.2 Design of the low level regulators
Recall that (see again Figure 1) the overall control action to be applied to the real system Σ has components generated by both the high-level and the low-level controllers, i.e.,
[TABLE]
The low level regulators are in charge of computing the local control corrections δui∈US,i⊖UˉS,i with the specific goal of compensating for the effect of the model inaccuracies at the high level expressed by the term wˉ(k) in (8).
To this end, first define the auxiliary system Σ^i: for h=kNL,…,(k+1)NL−1
[TABLE]
and note that Σ^i can be simulated in a centralized way in the time interval [kNL,(k+1)NL)
once the high level controller has computed uˉi[NL](k) at the beginning of the time interval.
Also denote by ΔΣi the model given by the difference between systems (1) and (16), with (2) and (15).
[TABLE]
where δxi(h)=xi(h)−x^i(h) , δzi(h)=zi(h)−z^i(h)
and δsi(h)=si(h)−s^i(h).
The difference state δxi is available at each time instant
h since xi is measurable and x^i can be computed with (16)
from the available value uˉi[NL](⌊h/NL⌋).
However, the difference dynamical system ΔΣi is
not yet useful for decentralized prediction since it depends upon the interconnection variables δsi(h) that, in turn, depend upon the variables δxj(h), j=i, not known in advance in the future prediction horizon. For this reason, we define a decentralized (approximated)
dynamical system ΔΣ^i with input δu^i(h) and discarding all coupling inputs, i.e.,
[TABLE]
The decentralized dynamical system ΔΣ^i is suitable for prediction, since it does not depend on quantities related to other subsystems. However, the dynamics of the subsystems can be very different from each other and resampling can be advisable for the design of the low level regulators. To this end, for any subsystem Σi, define a new sampling period ζi∈N+ such that NL/ζi=Ni∈N+ and a corresponding time index li. For clarity, the relations among the time scales with indices h, li, and k are shown in Figure 2 in a specific case.
In the time scale of li, define the dynamical system ΔΣ^i[ζi] as
[TABLE]
where BLii[ζi]=∑j=0ζi−1(ALii)jBLii. In the short time scale, with time index h, the input δu^i(h)=δu^i[ζi](⌊h/ζi⌋) is piecewise constant for
h∈[liζi,(li+1)ζi−1]. Its value will be computed as the result of a suitable optimization problem formulated for system (19).
Given δu^i(h), the evolution of δx^i(h) can be computed thanks to the dynamical model (18), however δx^i(h) is in general different from δxi(h) due to the neglected interconnections in systems (18) and (19). For this reason, and for all i=1,…,M, the input δui(h) to the real model (17) is computed based on δu^i(h), δxi(h), and δx^i(h) using a standard state-feedback policy, i.e.,
[TABLE]
where Ki is designed in such a way that the matrix FL=AL+BLK is Schur stable, being K=diag(K1,…,KM).
Assume now to be at time h=kNL and to have run the high level controller,
so that both uˉi[NL](k) and the
predicted value xˉ[NL](k+1∣k)=AHNLxˉ[NL](k)+BH[NL]uˉ[NL](k) are available. Therefore,
in order to remove the effect of the mismatch at the high level represented by wˉ(k) in (8),
the low level controller working in the interval h∈[kNL,(k+1)NL−1] should, if possible, aim to fulfill
[TABLE]
or equivalently,
[TABLE]
Since the model used for low-level control design is the decentralized one (i.e., (18)), the constraint (21) can only be formulated in an approximated way with reference to the state δx^i of system (18). In turn, if the resampling is used with ζi=1, the constraint on δx^i must be reformulated in terms of the state δx^i[ζi] of system (19), so that
[TABLE]
Note that the fulfillment of (22) does not
imply that (21) is satisfied due to the neglected
interconnections in (18) and (19) which make the term wˉ(k) in (5) not identically equal to zero, although it contributes to its reduction.
All these considerations lead to the formulation of the following low-level MPC designs. Letting δu^i[ζi](kNi:(k+1)Ni−1)=(δu^i[ζi](kNi),…,δu^i[ζi]((k+1)Ni−1))∈(Rmi)Ni,
the low level control action is computed, at time instant h=kNL, based on the solution to the following optimization problem:
[TABLE]
where JL is a positive definite function with arguments δx^i[ζi](kNi+li), δu^i[ζi](kNi+li), li=0,…,Ni−1, e.g.,
[TABLE]
and where a discussion on how to select the set ΔU^i is deferred to Appendix 7.2.
Finally, for j=0,…,NL−1, the control component at each (fast) time instant δui(kNL+j) is given by
[TABLE]
where
[TABLE]
A further clarification is finally due. At the low level, the controller is designed mostly to compensate for the effects of model inaccuracies during each long sampling time (i.e., that of the high-level controller). Therefore, the prediction horizon at low level coincides with one large sampling time of the slow high-level controller, i.e., corresponding to NL fast sampling times. Due to resampling, the optimization problem computed at the beginning of each slow sampling time, i.e., at time h=kNL, has a prediction horizon of Ni steps of lenght ζi, where indeed Niζi=NL. As a result of this optimization problem, the input sequence δu^i[ζi](kNi∣kNL), …, δu^i[ζi]((k+1)Ni−1∣kNL) is obtained, and at each fast sampling time kNL+j, j=0,…,NL−1, the real (low-level) input contribution (25) is used in (15). Note that, in this way, δui(h) varies at each sampling time.
In summary, the on-line implementation of the hierarchical control scheme here proposed consists of the following steps:
•
at any long sampling time k solve the centralized simplified slow time-scale MPC problem (11) with cost (12) and obtain the value of uˉ[NL](k) using (14);
•
for any system Σi, at the beginning of any long sampling time k solve the optimization problem (23) with cost (24) and compute the full sequence δu^i[ζi](kNi∣kNL), …, δu^i[ζi]((k+1)Ni−1∣kNL);
•
for any system Σi, at any fast sampling time kNL+jj=0,…,NL−1, compute the low-level contribution δui(kNL+j) using (25) and (26);
•
for any system Σi, at any fast sampling time kNL+j the applied control value is obtained as the sum of the low-level and of the high level contributions, as in (15).
4 Properties and algorithm implementation
The main assumptions, the recursive feasibility and convergence properties of the optimization problems stated at the high and low levels are established in this section.
4.1 Main assumptions and remarks
Define
[TABLE]
where
[TABLE]
Also, let Isi=[0nˉi×nˉ1⋯Inˉi⋯0nˉi×nˉM], A(NL)=AHNLβ−βALNL∈Rnˉ×n,
{\mathcal{R}(N_{\rm\scriptscriptstyle L})}=\big{[}\,B_{{\rm{\rm\scriptscriptstyle L}}}\ A_{{\rm{\rm\scriptscriptstyle L}}}B_{{\rm{\rm\scriptscriptstyle L}}}\ \cdots\ (A_{{\rm{\rm\scriptscriptstyle L}}})^{N_{\rm\scriptscriptstyle L}-1}B_{{\rm{\rm\scriptscriptstyle L}}}\,\big{]},
ρu be such that US⊆Bρu(0).
We now introduce the following assumption.
Assumption 3
∥ALNL∥<1;
2. 2.
letting Ri(Ni)=[BLii[ζi](ALii)ζiBLii[ζi]⋯(ALii)ζi(Ni−1)BLii[ζi]]
for each i=1,...,M, be the reachability matrix in Ni steps associated to ((ALii)ζi,BLii[ζi]),
matrix
[TABLE]
is full-rank with minimum singular value σHi(Ni)>0;
3. 3.
letting ρuˉ and ρδu^i be such that
Uˉ⊆Bρuˉ(0)
and ΔU^i⊇Bρδu^i(0),
respectively, for any i=1,...,M it holds that
[TABLE]
4. 4.
for each i=1,…,M
[TABLE]
5. 5.
Define ΔUˉi=ΔUi(NL−1), and ΔUi(j)=ΔU^i⊕BρΔui(j)(0) where ρΔui(j)=∑r=2j∥KiIsiFLj−r(AL−ALD)∥ρδx^(r−1) for all j=2,…,NL, ρΔui(j)=0 for j=0,1, and where ALD=diag(AL11,…,ALMM) and ρδx^(j)=∑i=1Mρδx^i2(j),
[TABLE]
We require that
[TABLE]
□**
Assumption 3 may be viewed, at a first glance, as a list of purely abstract and technical requirements. However, on the one hand, it represents relevant inherent properties of the system under control and, on the other hand, it implicitly includes important design principles, that will be shortly discussed next.
First, note that Assumption 3.(1) can always be verified in the light of the Assumption 1.(2) on stability of the transition matrix AL, and establishes a lower bound for the parameter NL. On the other hand, in view of Assumption 2.(2), βi are full rank and the full rank of matrix Ri(Ni) is guaranteed by the reachability of the local submodels used at low level control guaranteed by Assumption 1.(5): therefore Assumption 3.(2)
is fulfilled by taking Ni (and so NL) sufficiently large.
Secondly, Assumptions 3.(3)-3.(5) provide a tradeoff for the selection of parameters ρδu^i, i=1,…,M and ρuˉ. In fact
on the one hand, as required by Assumption 3.(5), and more specifically by equation (31), the input u(h) must be shared between the low-level controllers (with local inputs δui, whose maximal amplitude is related to ρδu^i) and centralized high-level controller (with input uˉ, whose maximal amplitude is related to ρuˉ);
on the other hand, the amplitude of δui must be sufficiently large to compensate for the model inaccuracies, as expressed by Assumptions 3.(3), 3.(4), and more specifically by equations (28) and (29).
Importantly, note that the constraints (28) and (29) depend upon quantities that are all functions of the number of steps NL, as clarified in the following.
In view of Assumptions 1.(2), 2.(3), and 2.(1), κ(NL)=∥∑j=0NL−1AHjBH−GH(1)−(β∑j=0NL−1ALjBL−G^L(1))∥ and
GH(1)=∑j=0+∞AHjBH, G^L(1)=β∑j=0+∞ALjBL. Therefore
κ(NL)≤∥∑j=NL+∞AHjBH∥+∥β∑j=NL+∞ALjBL∥≤∥AHNL∥∥GH(1)∥+∥ALNL∥∥β∥∥GL(1)∥,
where GL(z)=(zI−AL)−1BL. Therefore κ(NL)→0 exponentially as NL→+∞. This shows that also Assumption 3.(3) can be fulfilled by taking NL sufficiently large.
Similarly to Proposition 2.3 in [29], for any i=1,...,M it can be proved that
[TABLE]
This proves that also Assumption 3.(4) can be fulfilled - i.e., we can increase the maximum amplitude of uˉ as much as possible - by taking a sufficiently large value of parameter NL. This, as a byproduct, allows to minimize the conservativity of the overall control scheme, as discussed in the next section.
4.2 Main results and conservativity of the scheme
The size of the uncertainty set W to be considered in the high level design is given by
[TABLE]
where ρw=∑j=2NL∥βFLNL−j(AL−ALD)∥ρδx^(j−1).
The main result can now be stated.
Theorem 1
Under Assumption 3, if x(0) is such that the problem (11) is feasible at k=0 and, for all i=1,…,M
[TABLE]
*then
(i) wˉ(k)∈W and problems (11) and (23) are feasible for all k≥0;
(ii) for all h≥0*
[TABLE]
(iii) the state of the slow time-scale reduced model Σˉ[NL] enjoys robust convergence properties, i.e.,
[TABLE]
(iv) the state of the large scale model Σ enjoys robust convergence properties, i.e., for a computable positive constant ρx
[TABLE]
(v) we can define (see (67) in Appendix 7.3) a function σ(NL) of NL such that, if
[TABLE]
then, as k→∞, x(kNL)→xS.
□
Theorem 1 establishes three important facts. First, it shows that, if the initial state lies in a suitable set (and if Assumption 3 holds), the joint feasibility properties of the two control layers can be guaranteed in a recursive fashion. Secondly, it ensures robust convergence of the states of both the reduced-scale and the overall systems to a neighborhood of the corresponding steady-state goals.
Finally, if a suitably-defined function σ(NL) is smaller than one, then convergence of the state to the goal is ensured. The definition of σ(NL) is quite involved and for this reason it is given in Appendix 7.3, in particular see equation (67). A general remark, however, is due on parameter σ(NL): the more the subsystems are interconnected (in a wide sense, regarding both the existence of dependencies between subsystems and their amplitude), the larger σ(NL). On the other hand, it is possible to reduce arbitrarily this parameter by increasing the tuning knob NL. This depends of the fact that, as NL→+∞, both ∥A(NL)∥→0 and ∥B(NL)∥→0.
A further remark is that, similarly to Proposition 2.3 in [29], for any i=1,...,M it can be proved that
limNL→+∞λi(NL)=+∞, allowing to increase at will the feasibility region of the low-level problem.
From the discussion in Section 4.1 it has finally become clear that, by tuning the value of the low-level prediction horizon NL, one can reduce at will the values of ρδu^i, related to the maximum required amplitude of inputs δui. This, from (30) and (33), allows to reduce arbitrarily the high-level disturbance set W. This, in turn, allows to reduce at will the corresponding RPI set Z and to minimize the conservativity of the present control scheme. We remark that, although a fine tuning of the gain matrices Ki, i=1,…,M and KˉH can also be beneficial for the reduction of the conservativity of the scheme, the most relevant tuning knob indeed results parameter NL, especially since the dependence upon all other design parameters results rather straightforward and simple.
From the discussion above a further consideration is due. Although the case NL→+∞ allows to verify all the requirements, to obtain the best dynamic performances from the application of our control scheme, it would be beneficial to set NL to an ”average” value, such that the assumptions are verified, but at the same time that guarantees to control the system in a dynamic fashion also at high level. It is nevertheless important to remark that, when NL→+∞, the scheme can be regarded as a two-layer algorithm that at high level consists of a static optimizer based on a simplified system model and, at low level, consists of a dynamic, reactive, decentralized, and multi-rate, optimization-based regulator.
4.3 Design
The implementation of the multilayer algorithm described in the previous section requires a number of off-line computations here listed for the reader’s convenience.
•
design of AH, BH, and βi, i=1,...,M, such that Assumption 2 is satisfied;
•
design of KˉH such that both FH=AHNL+BH[NL]KˉH and
FL[NL]=ALNL+BL[NL]KˉHβ are Schur stable;
•
design of K=diag(K1,…,KM) such that FL=AL+BLK is Schur stable;
•
computation of ρδu^i, ρuˉi (see the procedure proposed in Appendix 7.2) and of the sets UˉS,i, ΔU^i;
computation of XˉF, Z, see [31], and PH with (13).
5 Simulation example and implementation procedures
The hierarchical control algorithm described in the previous sections has been used to control the model of the large-scale chemical plant described in [26], [35].
5.1 Description of the plant and linearized model
The system is composed of three reactors R1, R2, R3, three distillation columns C1, C2, C3, two recycle streams and six chemical components: A,B,C,D,E,F. The flow diagram of the system is reported in Figure 3, where D1, D2, D3 are the top products of the columns, while B1, B2, B3 are their bottom products. The following reactions occur inside the reactors R1:A+B→C+D, R2:D+E→F+B, and R3:D+E→F+B.
The system has six control variables, namely the refluxes (v1,v3,v5) and the vapour (v2,v4,v6) flow rates in the columns C1,C2,C3, respectively, and six outputs: the liquid molar fraction (r1) of component A at the top product of C1, the liquid molar fraction (r2) of component D at the bottom product of C1, the liquid molar fraction (r3) of component C at the top product of C2, the liquid molar fraction (r4) of component D at the bottom product of C2, the liquid molar fraction (r5) of component B at the top product of C3, and the liquid molar fraction (r6) of component F at the bottom product of C3. A detailed description of the model equations and of the model parameters is reported in [35]. The considered nominal operating point is characterized by the vector of inputs vnominal=[330410283385141282]′ lb mol/h, to which it corresponds the vector of outputs rnominal=[0.9420.5520.8270.9410.7050.991]′. The linearized model at this operating condition is of order n=192, and shows strong interactions among the control and controlled variables, see again [35].
In order to apply the hierarchical control structure described in this paper, the system and the corresponding linearized model have been partitioned into two subsystems (i.e. M=2). The first one includes reactors R1,R2 and columns C1,C2, while the second one is made by R3 and C3.
The continuous-time linearized models of the two subsystems have been discretized with
the algorithm described in [15] and basic sampling time T=0.05 h and a standard model order reduction procedure has been used to remove the unreachable states of the two linear subsystems. Denoting with δvi and δri the deviations of the inputs vi and the outputs ri, respectively, with respect to their nominal values, the first resulting linear reduced order model is of order n1=25, with m1=4 inputs, i.e., u1=(δv1,⋯,δv4), p1=4 outputs, i.e., y1=(δr1,⋯,δr4), and pz1=3 coupling outputs, while the second one has n2=16 states, m2=2 inputs, i.e., u2=(δv5,δv6), p2=2 outputs, i.e., y2=(δr5,δr6), and pz2=2 coupling outputs.
These two linear subsystems have been used as the linear models described in (1) for the implementation of the hierarchical control structure. The control variables are limited by
[TABLE]
for i=1,2, where uS,i is the steady state value corresponding to the output set-point values given by yS=10−3⋅(9.4,5.5,8.3,9.4,7.0,9.9).
5.2 Design and implementation procedures
The procedures used to implement the hierarchical control structure previously described and the computational details are now listed.
5.2.1 Tuning the hierarchical control scheme with NL=28, NH=3
Devising the high-level simplified model and the low-level submodels
•
The procedure described in Appendix 7.1 has been used to compute the matrices
β1 and β2, and the reduced order model (4)
with order nˉ=6. The dynamic matrix
AH=diag(0.972,0.984,0.969,0.969,0.874,0.869) has been computed; its parameters have been selected as the maximum singular values of the reachability matrix of each subsystem previously discretized. The matrix BH has been computed as described in Appendix 7.1. The resulting model has been resampled with NL=28 to obtain the model (5) to be used at the high level in the slow time scale.
•
The models (18) have been re-sampled with ζ1=4 and ζ2=1 to obtain the models (19) to be used at the low level in the fast time scale.
Off-line design of the low-level regulators
•
The low level finite-horizon optimization algorithms described in (23) and (24) have been implemented with
state and input penalties Q1=103⋅β1Tβ1, Q2=104⋅β2Tβ2,
R1=Im1 and R2=Im2.
•
The decentralized feedback gains K1 and K2 guaranteeing that FL is Schur stable can in principle be computed according to the algorithm described in [6] and based on the solution of LMI problems. However, in our case, we simply solved two independent LQ problems and we checked that the resulting FL is Schur stable.
•
The parameters ρδu^1, ρδu^2, ρuˉ1 and ρuˉ2 have been computed according to (41) in Appendix 7.2. Specifically we used the cost function Jρ=γ1\mathds11×Mρδu^−(ρuˉ−γ2ρu)′(ρuˉ−γ2ρu), where γ1, γ2 are positive constants selected as γ1=1 and γ2=0.3. Note that the choice of γ1,γ2 allows one to modify the feasibility region of the high level optimization problem (11) and (12); typically setting γ1=1, the feasibility properties of (11) grow with γ2.
The computed values are ρδu^1=145.2, ρδu^2=115.7, ρuˉ1=49.3, ρuˉ2=25.2. The input vectors of each subsystem, namely uˉi, δu^i, i=1,2, have been limited as follows: −ρuˉi/mi\mathds1mi×1≤uˉi−uˉS,i≤ρuˉi/mi\mathds1mi×1, −ρδu^i/mi\mathds1mi×1≤δu^i≤ρδu^i/mi\mathds1mi×1; and the corresponding sets UˉS,i, ΔU^S,i, i=1,2 have been obtained.
Note that ρuˉ1/m1+ρδu^1/m1=97.3 and ρuˉ2/m2+ρδu^2/m2=99.6 and we can write UˉS,1⊕ΔU^S,1={uˉ1+δu^1∣−97.3⋅\mathds1m1×1≤uˉ1+δu^1−uˉS,1≤97.3⋅\mathds1m1×1} and UˉS,2⊕ΔU^S,2={uˉ2+δu^2∣−99.6⋅\mathds1m2×1≤uˉ1+δu^2−uˉS,2≤99.6⋅\mathds1m2×1}. These results show that the size of both the sets UˉS,1⊕ΔU^S,1 and UˉS,2⊕ΔU^S,2 is close to that of the real input constraint in (36).
In addition, the radius of the sets BρΔui(NL−1)(0), i=1,2, resulting from the feedback policy in (20) for compensating the couplings terms between subsystems have been computed according to Assumption 3.(31) with ρΔu1(NL−1)=4.0 and ρΔu2(NL−1)=0.37. Therefore, we can also write
UˉS,1⊕ΔUS,1={uˉ1+δu1∣−99.3⋅\mathds1m1×1≤uˉ1+δu1−uˉS,1≤99.3⋅\mathds1m1×1} and UˉS,2⊕ΔUS,2={uˉ2+δu2∣−99.9⋅\mathds1m2×1≤uˉ1+δu2−uˉS,2≤99.9⋅\mathds1m2×1}.
In view of this, the conservativeness of the algorithm related to the computation of the input constraint sets described in Appendix 7.2 is small, especially for the second subsystem.
Off-line design of the high-level regulator
•
The high-level tube-based robust MPC has been designed according to the algorithm
described in [27] with state and input penalties QH=Inˉ and RH=0.1Im.
•
The gain matrix KˉH guaranteeing that both FH and FL[NL] are Schur stable can be computed according to the algorithm described in [6] and based on the solution of LMI problems. However, in our case, we simply solved a LQ problem and we checked that the resulting FH and FL[NL] were Schur stable.
•
The disturbance set has been obtained according to (33) and (30) with ρw=1.25, and the RPI set has been computed with the algorithms described at pp. 231-233 of [30], i.e., Z={z∣−(2.66,3.61,1.93,2.64,1.30,1.30)≤z≤\break(2.66,3.61,1.93,2.64,1.30,1.30)} and KˉHZ={KˉHz∣−(1.64,1.63,1.02,1.02,\break0.10,0.10)≤KˉHz≤(1.64,1.63,1.02,1.02,0.10,0.10)}. The terminal penalty PH is computed with (13), and the terminal set is calculated under nominal input constraints in (11),
i.e.,XˉF={xˉ[NL],o∣(xˉ[NL],o−xˉS)TPH(xˉ[NL],o−xˉS)≤0.997}, see [18].
The values taken by ρw, Z and KˉHZ reveal that the feasibility region of high-level regulator might be slightly reduced compared to the one of stabilizing MPC due to the use of tightened input constraint set, i.e., uˉS⊕USˉ⊖KˉHZ in optimization problem (11) and the computation formulas of the disturbance set in (33) and (30), but can be enlarged by increasing the tuning knobs ρuˉ1 and ρuˉ2.
5.2.2 Tuning the hierarchical control scheme with NL=84, NH=1
Along the same line with the previous Section 5.2.1, the computational results corresponding to tuning knob NL=84 and NH=1 are also listed here.
•
The following parameters ρδu^1, ρδu^2, ρuˉ1 and ρuˉ2 have been computed: ρδu^1=125.6, ρδu^2=98.3, ρuˉ1=59.5, ρuˉ2=41.9. The amplitude of both ρδu^1 and ρδu^2 (ρuˉ1 and ρuˉ2) is smaller (larger) than that with tuning knob NL=28, NH=3.
•
The sets UˉS,1⊕ΔU^S,1={uˉ1+δu^1∣−92.6⋅\mathds1m1×1≤uˉ1+δu^1−uˉS,1≤92.6⋅\mathds1m1×1} and UˉS,2⊕ΔU^S,2={uˉ2+δu^2∣−99.1⋅\mathds1m2×1≤uˉ1+δu^2−uˉS,2≤99.1⋅\mathds1m2×1}. These results show that the size of both the sets UˉS,1⊕ΔU^S,1 and UˉS,2⊕ΔU^S,2 has been reduced slightly compared to that with tuning knob NL=28, NH=3, however, is still close to that of the real input constraint in (36).
In addition, the radius of the sets BρΔui(NL−1)(0), i=1,2, have been computed with ρΔu1(NL−1)=11.3 and ρΔu2(NL−1)=0.91. Therefore, we can also write
UˉS,1⊕ΔUS,1={uˉ1+δu1∣−98.2⋅\mathds1m1×1≤uˉ1+δu1−uˉS,1≤98.2⋅\mathds1m1×1} and UˉS,2⊕ΔUS,2={uˉ2+δu2∣−99.8⋅\mathds1m2×1≤uˉ1+δu2−uˉS,2≤99.8⋅\mathds1m2×1}.
In view of this, the conservativity of the algorithm related to the computation of the input constraint sets described in Appendix 7.2 is still small, especially for the second subsystem.
•
The disturbance set has been obtained with ρw=1.56, and the RPI set has been computed i.e., Z={z∣−(1.83,2.11,1.69,1.82,1.58,1.58)≤z≤\break(1.83,2.11,1.69,1.82,1.58,1.58)} and KˉHZ={KˉHz∣−(0.37,0.29,0.17,0.18,\break0.01,0.01)≤KˉHz≤(0.37,0.29,0.17,0.18,0.01,0.01)}. The terminal penalty PH is computed with (13), and the terminal set is calculated under nominal input constraints in (11),
i.e.,XˉF={xˉ[NL],o∣(xˉ[NL],o−xˉS)TPH(xˉ[NL],o−xˉS)≤0.98}.
Note that, the size of both the sets Z and KˉHZ is smaller than that with the tuning NL=28, NH=3. The amplitude of ρuˉ1 and ρuˉ2 has been enlarged and the feasibility region of the high-level regulator can be increased significantly.
5.2.3 On-line implementation
The on-line implementation proceduce is depicted in Algorithm 1.
5.3 Simulation results: application to the linearized model
The overall control actions computed by the high and low level controllers
have been applied to the linear system at each fast sampling time. The output reference values for the linear system have been initially maintained at the nominal setpoints yS=10−3⋅(9.4,5.5,8.3,9.4,7.0,9.9); then, at time t=25.2 h, they have been set equal to
10−2⋅(2.43,−1.01,−0.43,0.15,−0.70,0.41). For comparison, a centralized state-feedback stabilizing MPC has been designed with an auxiliary state-feedback control law computed with LQ control, state penalty matrix Q=diag(Q1,Q2), input penalty R=Im and prediction horizon N=NH⋅NL=84. The terminal set has been chosen as XF={x∣(x−xS)TP(x−xS)≤1} where the terminal penalty P has been taken as the steady state solution of the Riccati equation according to the infinite horizon control problem with Q, R. All the simulation tests have been implemented
within MATLAB Yalmip and MPT toolbox, see [24] and [17], in a PC with Intel Core i5-4200U
2.30 GHz and with Windows 10 operating system. The SDPT3 solver has been used for the implementation of the centralized MPC and of the high-level regulator of the proposed approach, while the Matlab QUADPROG solver has been used for the low-level optimization problems. The detailed on-line computational time required for each controller is reported in Table 1. This table shows that the on-line computational time of the proposed hierarchical approach for each interval [kNL,kNL+NL) with NL=28 and NH=3 (NL=84 and NH=1), i.e., 2.49 s (6.87 s), is reduced dramatically compared to that of the centralized MPC, i.e., 21.9⋅NL s.
The evolution of the output and control variables of the controlled linear system is reported
in Figures 4-5 which show that, after an initial transient, inputs and outputs return to their nominal values until the change of the reference occurs, when both the centralized and the two-layer control systems properly react to bring the controlled variables to their reference values, and the performance of the proposed two-layer approach with both tuning cases NL=28, NH=3 and NL=84, NH=1 is close to that of centralized stabilizing MPC.
5.4 Simulation results: application to the nonlinear system
The two-layer control structure has also been applied to the nonlinear chemical plant with NL=84 and NH=1. Since in a realistic scenario the state is unmeasurable, the distributed Kalman filter described in [14] has been used. The covariances of the noises acting on the states have been set equal to Q^1=0.1In1, Q^2=0.1In2, while the covariances of the output noises have been chosen as R^1=0.01Ip1, R^2=0.01Ip2. Finally, the covariances of the initial state estimates have been selected as P1(0)=0.01In1, P2(0)=0.01In2.
Starting from the nominal operating conditions, the overall control actions computed by high and low level controllers
have been applied to the original nonlinear system at each fast time instant. The output reference values for the nonlinear system have been initially maintained at the nominal point rnominal; then, at time t=25.2 h, they have been set equal to
[0.866,0.558,0.791,0.903,0.704,0.948]. The evolution
of the output and control variables of the controlled nonlinear system are reported
in Figure 6-7. These figures show that, after an initial transient due to the state filter, inputs and outputs return to their nominal values until the change of the reference occurs, when the two-layer control system properly reacts to bring the controlled variables to their reference values.
6 Extensions and conclusions
In this paper a two-layer control scheme for systems made by interconnected subsystems has been presented. The algorithm is based on the solution, at the two layers, of MPC problems of reduced size and allows for a multirate implementation, suitable to deal with systems characterized by significantly different dynamics. Its main properties of recursive feasibility and convergence have been established and its performances have been tested in a nontrivial simulation example.
The main rationale of the proposed control architecture is grounded on the use, at the high level, of a reduced and slow-timescale model for centralized control. At the low-level, each subsystem is endowed of a local controller and is in charge of both compensating for the model inaccuracies introduced at the high level and dealing with the distributed nature of the system.
At high level, in this paper we use a tracking control algorithm based on robust tube-based MPC. This choice, however, is somehow arbitrary, since many alternative (robust) control solutions can be used instead, e.g., the offset-free robust MPC scheme proposed in [5].
The reference signals are here assumed to be previously computed, for instance by an additional Real Time Optimization layer (RTO), e.g., based on the current and predicted external conditions, such as prices of energy or costs of row materials, see [9, 19, 1, 20, 21, 13, 10, 37]. It is important to remark that RTO-based structures must be properly designed to guarantee the compatibility of the models used at the two layers, see [10, 37, 13]. Also, dynamic RTO structures may be prone to stability issues, as noted in [12].
To overcome problems related to the use of RTO to generate reference signals, a robust economic MPC approach can be used, see [3], for the design of a high level regulator which, at the same time, computes the optimal reference values for the controlled outputs.
This, in our opinion, would not entail significant differences in the algorithm implementation and in the theoretical results. Future work will be devoted to this extension.
7 Appendix
7.1 Construction of βi and of the reduced order model
A constructive procedure for the computation of the matrices βi, i=1,…,M, and the reduced order model satisfying Assumption 2 is listed here following the same line as in [29]. Note however that in [29] the case of dynamically decoupled subsystems was considered, and full system reduction actually could be carried out subsystem-by-subsystem. On the contrary, in this case, system reduction must be performed at a full system level and structural additional constraints must be satisfied, i.e. the block-diagonality of matrix β. In view of the presence of these constraints, the sufficient and necessary condition given in Proposition 1 in [29] for guaranteeing the fulfillment of Assumption 2 (i.e., that nˉ≥p+dim(ImGLx∩KerCL)) is, in our case, only necessary.
Here we now describe a (possibly conservative) procedure:
a
find a subspace Kerβi of dimension ni−nˉi
so that Kerβi⊆KerCLii
and Zii∩Kerβi={0}, where Zii=ImG~Lii(1)∩KerCLii, and G~Lii(z)=(zI−ALii)−1[BLiiEi];
b
let {κ1,…,κni−nˉi} be a
set of independent vectors in Kβ,i=Kerβi and complete this
set to a basis Bi={v1,…,vnˉi,κ1,…,κni−nˉi}
of the whole space Rni;
c
let {r1,…,rnˉi} be a basis of Rnˉi,
define
[TABLE]
and TL be the matrix whose columns
are the vectors in Bi, then βi=β^iTL−1;
d
define collective matrix β=diag(β1,…,βM);
e
choose matrix AH being Schur stable, and let
[TABLE]
where the suitable choices for AH are those modeling the dominant dynamics of the low-level collective model.
Steps a-c imply that Assumption 2.(2) is fulfilled, while step e guarantees that Assumption 2.(1) and 2.3 are satisfied.
Note that a less conservative choice can be taken, i.e., defining, in step a, G~Lii(z)=(zI−ALii)−1BLii. However this choice does not guarantee a-priori that the property Z∩Kβ={0}, where Z=ImGLx∩KerCL, GLx=(I−AL)−1BL, and Kβ=∏i=1MKβi. This must be verified after the reduction phase has been carried out.
7.2 Computation of the input constraint sets
In the scheme proposed in this paper, the dimensions of the input constraint sets UˉS,i and ΔU^i are key tuning knobs, which must be selected in order to satisfy, at the same time, the inequalities (28), for all i=1,…,M, and (31). To address the design issue, in this appendix we propose a simple and lightweight algorithm based on a linear program. As a simplifying assumption, we set ΔU^i=Bρδu^i(0) and UˉS,i=Bρuˉi(0). Under this assumption, the tuning knobs are the vectors ρδu^=(ρδu^1,…,ρδu^M) and ρuˉ=(ρuˉ1,…,ρuˉM). Note that, in case of need, such assumption can be relaxed, at the price of a slightly different definition of the inequalities below.
First consider inequality (28), to be verified for all i=1,…,M. Here the constant ρuˉ appears, defined in such a way that UˉS=∏i=1MBρuˉi(0)⊆Bρuˉ(0). We can define, for example, ρuˉ=∑i=1Mρuˉi2≤∑i=1Mρuˉi. Therefore, to fulfill (28) it is sufficient to verify the following matrix inequality
[TABLE]
where \mathds1M×M is the M×M matrix whose entries are all equal to 1.
The second main inclusion to be fulfilled is (31), which is verified if, for all i=1,…,M,
[TABLE]
By definition, ΔUˉi=ΔU^i⊕BρΔui(NL−1)(0), where
[TABLE]
where λij=∑r=2NL−1∥KiIsiFLNL−r−1(AL−ALD)∥∑k=1r−1∥(ALjj)r−1−kBLjj∥. This implies that
ΔUˉi⊆Bρδu^i+∑j=1Mλijρδu^j(0). Therefore, to verify (38) it is sufficient to enforce the constraint
[TABLE]
where Λ is the M×M matrix whose entries are λij, i,j=1,…,M, while ρu=(ρu1,…,ρuM), where Bρi(0)⊆US,i for all i=1,…,M.
Eventually, a suitable choice of ρδu^ and ρuˉ is obtained as the solution to the following problem:
[TABLE]
where Jρ is a suitable (linear or quadratic, if possible) cost function that allows to maximize the size of the constraint set.
The proof of Theorem 1 lies on the intermediary results stated below.
Proposition 1
A) Under Assumption 3 and if xˉ[NL](k)=βx(kNL),
then for any initial condition x^(kNL)=x(kNL) such that, for all i=1,…,M
[TABLE]
*and for any uˉ[NL]∈uˉS⊕UˉS
there exists a feasible sequenceδu^i[ζi](kNi:(k+1)Ni−1∣kNL)∈ΔU^iNi
such that the terminal constraint (22) is satisfied.
B) if x(kNL) satisfies condition (42), ∥ALNL∥<1,
and, for all i=1,…,M, (29) is verified, then recursive feasibility of the terminal constraint (22) is guaranteed.*
Therefore, in view of (43), (44a), (45a), and the definition
A(NL), B(NL), Isi, the constraint (22)
can be written as
[TABLE]
From this expression, the definitions of σHi(Ni),
ρuˉ, ρδu^i, and in view of (27),
it can be concluded that a feasible sequence δu^i[ζi](kNi:(k+1)Ni−1∣kNL)
can be computed provided that
Defining the collective vectors x^=(x^1,…,x^M), δx=(δx1,…,δxM),δx^=(δx^1,…,δx^M), and ε(kNL+j∣kNL)=δx(kNL+j∣kNL)−δx^(kNL+j∣kNL), we have that wˉ(k)=βx(kNL+NL)−xˉ[NL](k+1∣k)=βx^(kNL+NL)+βδx(kNL+NL)−xˉ[NL](k+1∣k)=(βx^(kNL+NL)+βδx^(kNL+NL)−xˉ[NL](k+1∣k))+βε(kNL+j∣kNL)=βε(kNL+j∣kNL). The latter equality holds in view of the fact that Problem (23) is feasible, and therefore equality (22) is verified.
From (17), (18), (20), we collectively have that
[TABLE]
In view of the fact that ε(kNL∣kNL)=δx^(kNL∣kNL)=0, then
wˉ(k)=\breakβ∑j=2NLFLNL−j(AL−ALD)δx^(kNL+j−1∣kNL). From this it follows that
[TABLE]
Since δu^i are bounded for all i=1,…,M, i.e., scalar ρδu^i are defined such that δu^i∈Bρδu^i(0). In view of this, we compute that
∥δx^(kNL+j∣kNL)∥≤ρδx^(j), where ρδx^i(j) is defined in (30). Therefore, δx^(kNL+j∣kNL) are bounded, for all j=1,…,NL−1 and more specifically we get that ∥δx^(kNL+j∣kNL)∥≤\break∑i=1Mρδx^i2(j).
Therefore one has (52) for all k≥0.
From (55) we have that ε(kNL+j∣kNL)=∑r=2jFLj−r(AL−ALD)δx^(kNL+r−1∣kNL) and therefore
δui(kNL+j)−δu^i(kNL+j∣kNL)=Kiεi(kNL+j∣kNL)=\breakKiIsi∑r=2jFLj−r(AL−ALD)δx^(kNL+r−1∣kNL). From this it follows that
δui(kNL+j)∈δu^i(kNL+j∣kNL)⊕BρΔui(j)(0) and therefore δui(kNL+j)∈ΔUi(j). In view of the monotonicity property ρΔui(j+1)≥ρΔui(j) for all j, it holds thatBρΔui(j+1)(0)⊇BρΔui(j)(0), which implies (54). □
(i) If ∥x(0)−xS∥≤λi(NL) and recalling that Assumption 3 holds, from Proposition 1, recursive feasibility of the optimization problems (23) is guaranteed, i.e., that there exists, for all k≥0, a feasible sequence δu^i[ζi](kNi:(k+1)Ni−1∣kNL)∈ΔU^iNi
such that the terminal constraint (22) is satisfied.
Also, from Proposition 2, it follows that wˉ(k)∈W for all k≥0, which allows to apply the recursive feasibility arguments of [27], proving that also (11) enjoys recursive feasibility properties.
(ii) It is now possible to conclude that, in view of the feasibility of (11), uˉ[NL](k)∈uˉS⊕UˉS; also, from Proposition 2 it follows that δui(kNL+j)∈ΔUˉi for all k≥0, j=0,…,NL−1, and i=1,…,M. From this, under (31), the inclusion (34) can also be proved.
(iii) We apply the results in [27], which guarantee robust convergence properties. In other words, it holds that xˉ[NL],o(k)→xˉS as k→+∞, and that xˉ[NL](k) is asymptotically driven to lie in the robust positively invariant set xˉS⊕Z.
(iv) To show robust convergence of the global system state, from (3) we obtain that
[TABLE]
Denoting BL,NLC=[ALNL−1BL…BL], we obtain that
[TABLE]
Also, recall that
δu(kNL:kNL+NL−1)=δu^(kNL:kNL+NL−1)+\breakdiag(K,…,K)ε(kNL:kNL+NL−1) and that, by defining
Recall that xˉ[NL],o(k)→xˉS, uˉ[NL],o(k)→uˉS as k→+∞. Also, we compute that
[TABLE]
Based on this, we define κδu=∥BL,NLC(I+diag(Ki,…,Ki)FNLdiag(ALC,…,ALC)BL)∥ and we write (59) as
[TABLE]
where ∥wL(k)∥≤κδuNLmaxh∈{kNL,…,(k+1)NL−1}∥δu^(h)∥≤κδuNL∑i=1Mρδu^i2.
Since xS=ALNLxS+BL[NL]uˉS and BL[NL]KˉHxˉS=BL[NL]KˉHβxS, we can rewrite (61) as
[TABLE]
Eventually, since FL[NL]=ALNL+BL[NL]KˉHβ is Schur stable, then the asymptotic result follows, where
ρx=κδuNL∑i=1Mρδu^i2.
(R1,…,RM), Q=diag(Q1,…,QM), δxˉend=xˉ[NL](k+1∣k)−βx^(kNL+NL) andRD(NL)=[(ALD)NL−1BL…BL]. Recalling that bin>0 elementwise, in view of continuity arguments, there exists a ball Bρend(0) for δxˉend such that the solution to problem (62) satisfies Ainδu^(kNL:kNL+NL−1)≤bin. If
δxˉend∈Bρend(0), then
[TABLE]
Also, it holds that the optimal constrained solution fulfills
[TABLE]
and therefore ∥δu^(kNL:kNL+NL−1)∥≤λmin(Q)JLmax; in view of this, for ∥δxˉend∥∈Bρend(0), then ∥δu^(kNL:kNL+NL−1)∥≤λmin(Q)JLmaxρend1∥δxˉend∥.
Defining now
κˉ=max{[I0][2QβRD(NL)(βRD(NL))T0]−1[0I],λmin(Q)JLmaxρend1} then we conclude that
where ∥wx(k)∥≤κx(NL)∥x(kNL)−xS∥, κx(NL)=κuκˉ∥BL,NLC∥∥A(NL)+B(NL)KˉHβ∥, ∥wo(k)∥≤κuo∥uˉ[NL],o(k)−uˉS∥+κxo∥xˉ[NL],o(k)−xˉS∥, κuo=∥BL[NL]∥+κuκˉ∥BL,NLC∥∥B(NL)∥, and
[TABLE]
To derive a stability condition, we recast (65) as the following redundant dynamic system
[TABLE]
where the initial conditions for δx1 and δx2 coincide and are equal to x(0)−xS, and where ∥wδxi∥≤κx(NL)∥δxi∥, i=1,2. Recall also that, as already discussed, wo is an asymptotically vanishing input.
The stability of the system (66) can be studied using the (ISS) small gain theorem in [11], according to which the interconnected system above enjoys asymptotic stability properties if the matrix
[TABLE]
is Schur stable, where
[TABLE]
Also, Γ is Schur stable if and only if the inequality (35) is verified.
Then, under the latter condition, x(kNL)→xS as k→+∞. □
Bibliography38
The reference list from the paper itself. Each links out to its DOI / PubMed record.
1[1] V. Adetola and M. Guay. Integration of real-time optimization and model predictive control. Journal of Process Control , 20:125 – 133, 2010.
2[2] M. Baldea and P. Daoutidis. Dynamics and Nonlinear Control of Integrated Process Systems . Cambridge University Press, 2012.
3[3] F. A. Bayer, M. A. Müller, and F. Allgöwer. Tube-based robust economic model predictive control. Journal of Process Control , 24(8):1237–1246, 2014.
4[4] B.W. Bequette. Nonlinear predictive control using multi-rate sampling. The Canadian Journal of Chemical Engineering , 69(1):136–143, 1991.
5[5] G. Betti, M. Farina, and R. Scattolini. A robust MPC algorithm for offset-free tracking of constant reference signals. IEEE Transactions on Automatic Control , 58(9):2394–2400, 2013.
6[6] G Betti, M Farina, and R Scattolini. Distributed MPC: A noncooperative approach based on robustness concepts. In Distributed Model Predictive Control Made Easy , pages 421–435. Springer, 2014.
7[7] X. Chen, M. Heidarinejad, J. Liu, and P. D. Christofides. Composite fast slow MPC design for nonlinear singularly perturbed systems. AI Ch E Journal , 58:1082––1811, 2012.
8[8] X. Chen, M. Heidarinejad, J. Liu, D. Muñoz de la Peña, and P. D. Christofides. Model predictive control of nonlinear singularly perturbed systems: application to a large-scale process network. Journal of Process Control , 21:1296––1305, 2011.