A hierarchical MPC scheme for interconnected systems

TL;DR

This paper proposes a hierarchical MPC control scheme for interconnected systems, combining a robust, long-term optimized higher layer with fast local regulators for subsystem refinement, demonstrated through building temperature control simulations.

Contribution

It introduces a novel hierarchical control architecture that integrates reduced-order and full-order models for improved interconnected system management.

Findings

Effective temperature regulation in building simulation.

Hierarchical MPC improves control performance.

Combines robustness with fast local adjustments.

Abstract

This paper describes a hierarchical control scheme for interconnected systems. The higher layer of the control structure is designed with robust Model Predictive Control (MPC) based on a reduced order dynamic model of the overall system and is aimed at optimizing long-term performance, while at the lower layer local regulators acting at a higher frequency are designed for the full order models of the subsystems to refine the control action. A simulation experiment concerning the control of the temperature inside a building is reported to witness the potentialities of the proposed approach.

Equations

$${\mathcal{B}}_{{\rho }_{{\epsilon }_i}}(\bar{x}):=\{x\in {\mathbb{R}}^{dim}:||x-\bar{x}||\le {\rho }_{{\epsilon }_i}\}$$

$$\Sigma_{i}:\ \left\{\begin{array}[]{lcl}x_{i}(h+1)&=&A_{\rm\scriptscriptstyle L}^{ii}x_{i}(h)+B_{\rm\scriptscriptstyle L}^{ii}u_{i}(h)+E_{\rm\scriptscriptstyle L}^{i}s_{i}(h)\\ z_{i}(h)&=&C_{\rm\scriptscriptstyle L}^{zi}x_{i}(h),\end{array}\right.$$

$$s_i(h)=\sum _{j=1}^M{L}_{ij}{z}_j(h)$$

$$\Sigma :x(h+1)=A_{\mathrm{L}}x(h)+B_{\mathrm{L}}u(h)$$

$$\bar{\Sigma}:\ \begin{array}[]{c}\bar{x}(h+1)=A_{\rm\scriptscriptstyle H}\bar{x}(h)+B_{\rm\scriptscriptstyle H}\bar{u}(h)\end{array}$$

$$\bar{\Sigma}^{{[N_{\rm\scriptscriptstyle L}]}}:\ \begin{array}[]{l}\bar{x}^{{[N_{\rm\scriptscriptstyle L}]}}(k+1)=A_{{\rm{\rm\scriptscriptstyle H}}}^{N_{\rm\scriptscriptstyle L}}\bar{x}^{{[N_{\rm\scriptscriptstyle L}]}}(k)+B_{{\rm{\rm\scriptscriptstyle H}}}^{{[N_{\rm\scriptscriptstyle L}]}}\bar{u}^{{[N_{\rm\scriptscriptstyle L}]}}(k)\end{array}$$

$${\bar{x}}_i^{[N_{\mathrm{L}}]}(k)={\beta }_i{x}_i(k{N}_{\mathrm{L}})$$

$${\bar{x}}^{[N_{\mathrm{L}}]}(k)=\beta x(k{N}_{\mathrm{L}})$$

$$\bar{\Sigma}_{w}^{{[N_{\rm\scriptscriptstyle L}]}}:\ \begin{array}[]{l}\bar{x}^{{[N_{\rm\scriptscriptstyle L}]}}(k+1)=A_{{\rm{\rm\scriptscriptstyle H}}}^{N_{\rm\scriptscriptstyle L}}\bar{x}^{{[N_{\rm\scriptscriptstyle L}]}}(k)+B_{{\rm{\rm\scriptscriptstyle H}}}^{{[N_{\rm\scriptscriptstyle L}]}}\bar{u}^{{[N_{\rm\scriptscriptstyle L}]}}(k)+\bar{w}(k)\end{array}$$

$$\bar{\Sigma}_{w}^{{[N_{\rm\scriptscriptstyle L}]},o}:\ \begin{array}[]{l}\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(k+1)=A_{{\rm{\rm\scriptscriptstyle H}}}^{N_{\rm\scriptscriptstyle L}}\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)+B_{{\rm{\rm\scriptscriptstyle H}}}^{{[N_{\rm\scriptscriptstyle L}]}}\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)\end{array}$$

$$\bar{\Sigma}_{w}^{{[N_{\rm\scriptscriptstyle L}]},e}:\ \begin{array}[]{l}\bar{e}(k+1)=F_{\rm\scriptscriptstyle H}\bar{e}(k)+\bar{w}(k)\end{array}$$

$$\begin{array}[]{l}\min_{\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(t),\overrightarrow{\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}}{(t:t+N_{\rm\scriptscriptstyle H}-1)}}{J_{{\rm\scriptscriptstyle H}}\big{(}\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(t),\overrightarrow{\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}}{(t:t+N_{\rm\scriptscriptstyle H}-1)})}\\ \mbox{subject to:}\\ \bullet\,\,\mbox{the unperturbed model dynamics}~{}\eqref{Eqn:C HLS model_unp}\\ \bullet\,\,\mbox{the initial constraint}~{}\beta x(tN_{\rm\scriptscriptstyle L})-\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(t)\in\mathcal{Z}\\ \bullet\,\,\mbox{the terminal constraint}~{}\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(t+N_{\rm\scriptscriptstyle H})\in\bar{\mathcal{X}}_{F}\\ \bullet\,\,\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)\in{\bar{\mathcal{U}}\ominus\bar{K}_{\rm\scriptscriptstyle H}\mathcal{Z}},\,k=t,\dots,t+N_{\rm\scriptscriptstyle H}-1,\end{array}$$

$$J_{\mathrm{H}}=\sum _{k=t}^{t+N_{\mathrm{H}}-1}\Vert {\bar{x}}^{[N_{\mathrm{L}}],o}(k){\Vert }_{Q_{\mathrm{H}}}^2+\Vert {\bar{u}}^{[N_{\mathrm{L}}],o}(k){\Vert }_{R_{\mathrm{H}}}^2+\Vert {\bar{x}}^{[N_{\mathrm{L}}],o}(t+N_{\mathrm{H}}){\Vert }_{P_{\mathrm{H}}}^2$$

$$F_{\mathrm{H}}^T{P}_{\mathrm{H}}{F}_{\mathrm{H}}-P_{\mathrm{H}}=-(Q_{\mathrm{H}}+{\bar{K}}_{\mathrm{H}}^T{R}_{\mathrm{H}}{\bar{K}}_{\mathrm{H}})$$

$${\bar{u}}^{[N_{\mathrm{L}}]}(t)={\bar{u}}^{[N_{\mathrm{L}}],o}(t|t)+{\bar{K}}_{\mathrm{H}}(\beta x(t{N}_{\mathrm{L}})-{\bar{x}}^{[N_{\mathrm{L}}],o}(t|t))$$

$$u_i(h)={\bar{u}}_i^{[N_{\mathrm{L}}]}(\lfloor h/N_{\mathrm{L}}\rfloor )+\delta u_i(h)$$

$$\hat{\Sigma}_{i}:\ \left\{\begin{array}[]{lcl}\hat{x}_{i}(h+1)&=&A_{\rm\scriptscriptstyle L}^{ii}\hat{x}_{i}(h)+B_{\rm\scriptscriptstyle L}^{ii}\bar{u}_{i}^{{[N_{\rm\scriptscriptstyle L}]}}(\lfloor h/N_{\rm\scriptscriptstyle L}\rfloor)+E_{\rm\scriptscriptstyle L}^{i}\hat{s}_{i}(h)\\ \hat{s}_{i}(h)&=&\sum_{j=1}^{M}L_{ij}\hat{z}_{j}(h)\\ \hat{z}_{i}(h)&=&C_{\rm\scriptscriptstyle L}^{zi}\hat{x}_{i}(h)\\ \hat{x}_{i}(kN_{\rm\scriptscriptstyle L})&=&x_{i}(kN_{\rm\scriptscriptstyle L})\end{array}\right.$$

$$\varDelta\Sigma_{i}:\ \left\{\begin{array}[]{lcl}\delta x_{i}(h+1)&=&A_{\rm\scriptscriptstyle L}^{ii}\delta x_{i}(h)+B_{\rm\scriptscriptstyle L}^{ii}\delta u_{i}(h)+E_{\rm\scriptscriptstyle L}^{i}\delta s_{i}(h)\\ \delta{s}_{i}(h)&=&\sum_{j=1}^{M}L_{ij}\delta{z}_{j}(h)\\ \delta z_{i}(h)&=&C_{\rm\scriptscriptstyle L}^{zi}\delta x_{i}(h)\\ \delta x_{i}(kN_{\rm\scriptscriptstyle L})&=&0\end{array}\right.$$

$$\varDelta\hat{\Sigma}_{i}:\ \left\{\begin{array}[]{lcl}\delta\hat{x}_{i}(h+1)&=&A_{\rm\scriptscriptstyle L}^{ii}\delta\hat{x}_{i}(h)+B_{\rm\scriptscriptstyle L}^{ii}\delta\hat{u}_{i}(h)\\ \delta\hat{x}_{i}(kN_{\rm\scriptscriptstyle L})&=&0\end{array}\right.$$

$$\delta u_i(h)=\delta {\hat{u}}_i(h)+K_i(\delta x_i(h)-\delta {\hat{x}}_i(h))$$

$${\beta }_i{x}_i(k{N}_{\mathrm{L}}+N_{\mathrm{L}})={\bar{x}}_i^{[N_{\mathrm{L}}]}(k+1|k)$$

$${\beta }_i\delta x_i(k{N}_{\mathrm{L}}+N_{\mathrm{L}})={\bar{x}}_i^{[N_{\mathrm{L}}]}(k+1|k)-{\beta }_i{\hat{x}}_i(k{N}_{\mathrm{L}}+N_{\mathrm{L}})$$

$${\beta }_i\delta {\hat{x}}_i(k{N}_{\mathrm{L}}+N_{\mathrm{L}})={\bar{x}}_i^{[N_{\mathrm{L}}]}(k+1|k)-{\beta }_i{\hat{x}}_i(k{N}_{\mathrm{L}}+N_{\mathrm{L}})$$

$$\begin{array}[]{l}\min_{\overrightarrow{\delta\hat{u}}_{i}{(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1)}}J_{{\rm\scriptscriptstyle L}}\big{(}\delta\hat{x}_{i}(kN_{\rm\scriptscriptstyle L}),\overrightarrow{\delta\hat{u}}_{i}{(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1)})\\ \mbox{subject to:}\\ \bullet\,\,\mbox{the dynamics}~{}\eqref{Eqn:LL dndyn}\\ \bullet\,\,\mbox{the terminal constraint}~{}\eqref{eq:terminaln}\\ \bullet\,\,\delta\hat{u}_{i}(kN_{\rm\scriptscriptstyle L}+j)\in{\varDelta\hat{\mathcal{U}}}_{i},\,j=0,\dots,N_{\rm\scriptscriptstyle L}-1,\end{array}$$

$$J_{\mathrm{L}}=\sum _{j=0}^{N_{\mathrm{L}}-1}\Vert \delta {\hat{x}}_i(k{N}_{\mathrm{L}}+j){\Vert }_{Q_i}^2+\Vert \delta {\hat{u}}_i(k{N}_{\mathrm{L}}+j){\Vert }_{R_i}^2$$

$$\kappa (N_{\mathrm{L}})=\Vert \mathcal{B}(N_{\mathrm{L}})\Vert$$

$$\mathcal{B}(N_{\mathrm{L}})=\sum _{j=1}^{N_{\mathrm{L}}}{A}_{\mathrm{H}}^{N_{\mathrm{L}}-j}{B}_{\mathrm{H}}-\beta \sum _{j=1}^{N_{\mathrm{L}}}{A}_{\mathrm{L}}^{N_{\mathrm{L}}-j}{B}_{\mathrm{L}}$$

$${\mathcal{H}}_i(N_{\mathrm{L}})={\beta }_i{\mathcal{R}}_i(N_{\mathrm{L}})\in {\mathbb{R}}^{{\bar{n}}_i\times N{m}_i}$$

$${\rho }_{\delta {\hat{u}}_i}>\frac{\kappa (N_{\mathrm{L}}){\rho }_{\bar{u}}}{\sqrt{N_{\mathrm{L}}}{\underline{\sigma }}_{{\mathcal{H}}_i(N_{\mathrm{L}})}}$$

$${\chi }_i(k{N}_{\mathrm{L}})=\frac{\sqrt{N_{\mathrm{L}}}{\varrho }_u\Vert \mathcal{R}(N_{\mathrm{L}})\Vert \Vert \mathcal{A}(N_{\mathrm{L}})\Vert }{(1-\Vert A_{\mathrm{L}}^{N_{\mathrm{L}}}\Vert )(\sqrt{N_{\mathrm{L}}}{\underline{\sigma }}_{{\mathcal{H}}_i(N_{\mathrm{L}})}{\rho }_{\delta {\hat{u}}_i}-\kappa (N_{\mathrm{L}}){\rho }_{\bar{u}})}\le 1$$

Full text

1 Introduction and main idea

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 represents a problem of increasing difficulty and requires innovative solutions. A classical approach to deal with this challenge consists of resorting to hierarchical control structures, where at the higher layer of the hierarchy simplified models are used to predict and control the long term behavior of the overall system, while at the lower layer local control actions are designed to compensate for model inaccuracies, disturbances, or parametric variations. Along this line, many hierarchical control methods have been described in the past, see e.g. Adetola and Guay (2010), Kadam and Marquardt (2007) in the context of Real Time Optimization (RTO), or Amrit et al. (2011); Grüne (2013); Diehl et al. (2011) in the emerging area of economic MPC.

In view of the potentialities of multilayer control structures, this paper describes a novel approach to the design of a hierarchical control structure for large scale systems composed by interconnected subsystems. The scheme of the proposed solution is sketched in Figure 1: the system under control $\Sigma$ is composed of $M$ interconnected subsystems $\Sigma_{1},...,\Sigma_{M}$. A reduced order model $\bar{\Sigma}_{i}$, $i=1,...,M$ is computed for each subsystem, and the overall reduced order model $\bar{\Sigma}$ is obtained; typically $\bar{\Sigma}_{i}$ and $\bar{\Sigma}$ represent low-frequency approximations of the corresponding systems. At a slow sampling rate, a centralized MPC regulator $R_{\rm\scriptscriptstyle H}$ is designed for $\bar{\Sigma}$ to consider the long-term behavior of the controlled system and to compute the control variables $\bar{u}_{i}$, $i=1,...,M$. Then, local regulators $R_{Li}$, $i=1,...,M$, working at a faster time scale, are designed for each subsystem $\Sigma_{{}_{i}}$: their scope is to compute the control contributions $\delta u_{i}$ compensating for the inaccuracies in the high layer design due to the mismatch between $\Sigma$ and $\bar{\Sigma}$. This structure has already been studied in Picasso et al. (2016) where, however, only independent systems $\Sigma_{i}$ with joint output constraints were considered. The advantage of the approach here proposed is twofold: first, at the slower time scale the optimization problem underlying the MPC solution is of reduced dimension and can minimize a global cost function over a long horizon with a limited computational cost; second, also the local regulators designed for the local subsystems involve the solution to optimization problems whose complexity only depends on the order of the local submodels.

The paper is organized as follows. In Section 2 the models considered at the two layers of the control structure are introduced. Section 3 describes the MPC algorithms adopted at the two layers, while Section 4 is presents the main feasibility and convergence results as well as a summary of the main steps to be performed in the algorithm implementation. Section 5 describes a 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 a set of variables $z_{i}\in{\mathbb{R}}^{q_{i}}$, $i=1,2,\dots,M$, we define the vector whose vector-components are $z_{i}$ in the following compact form: $(z_{1},z_{2},\cdots,z_{M})=[\,z_{1}^{T}\ z_{2}^{T}\ \cdots\ z_{M}^{T}\,]^{T}\in{\mathbb{R}}^{q}$, where $q=\sum_{i=1}^{M}q_{i}$. The symbols $\oplus$/$\ominus$ denote the Minkowski sum/difference. We denote with $\|\cdot\|$ the Euclidean norm. Finally, a ball with radius $\rho_{\varepsilon_{i}}$ and centered at $\bar{x}$ in the ${\mathbb{R}}^{dim}$ space is defined as follows

[TABLE]

2 Models for the two-layer control scheme

In this section we present the model of the complex system under study and the simplified one used for high-level control.

2.1 Large-scale system model

In line with Lunze (1992), we assume that the overall system $\Sigma$ is composed by $M$ discrete-time, linear, interacting subsystems described by

[TABLE]

$i=1,2,\dots,M$, where $x_{i}\subseteq{\mathbb{R}}^{n_{i}}$, $u_{i}\in{\mathcal{U}}_{i}\subseteq{\mathbb{R}}^{m_{i}}$ are the state and input vectors, and where the interconnections among the $\Sigma_{i}^{{}^{\prime}}s$ are represented by the coupling input and output vectors $s_{i}\in{\mathbb{R}}^{p_{si}}$ and $z_{i}\in{\mathbb{R}}^{p_{zi}}$ linked through the following expression

[TABLE]

with $L_{ii}=0$, $i=1,...,M$. The sets ${\mathcal{U}}_{i}$ are closed and convex sets containing the origin.

Collecting all the subsystems (1), the overall dynamical model of $\Sigma$ is

[TABLE]

where $x=(x_{1},\ldots,x_{M})\in\mathbb{R}^{n}$, $n=\sum_{i=1}^{M}n_{i}$, $u=(u_{1},\ldots,u_{M})\in\mathbb{R}^{m}$, $m=\sum_{i=1}^{M}m_{i}$. The diagonal blocks of $A_{L}$ are state transition matrices $A_{\rm\scriptscriptstyle L}^{ii}$, whereas the coupling terms among the subsystems correspond to the non-diagonal blocks of $A_{\rm\scriptscriptstyle L}$, i.e., $A_{\rm\scriptscriptstyle L}^{ij}=E_{\rm\scriptscriptstyle L}^{i}L_{ij}C_{\rm\scriptscriptstyle L}^{zj}$, with $j\neq i$. The collective input matrix is $B_{\rm\scriptscriptstyle L}=$diag$(B_{\rm\scriptscriptstyle L}^{11},...,\,B_{\rm\scriptscriptstyle L}^{MM})$. We define ${\mathcal{U}}=\prod_{i=1}^{M}{\mathcal{U}}_{i}\subseteq\mathbb{R}^{m}$, which is convex by the convexity of $\mathcal{U}_{i}$.

Concerning systems (1) and (3), the following standing assumption is introduced:

Assumption 1

the state $x_{i}$ is measurable, for each $i=1,\dots,M$; 2. 2.

$A_{\rm\scriptscriptstyle L}$* is Schur stable;* 3. 3.

the pair $(A_{\rm\scriptscriptstyle L}^{ii},\,B_{\rm\scriptscriptstyle L}^{ii})$ is reachable, for each $i=1,\dots,M$. $\square$

Note that Assumption 1.1 can be relaxed, using suitable observers - e.g., distributed ones - and accounting for state estimation error in a rigorous way.

2.2 Reduced order models

Associated with each subsystem $\Sigma_{i},\,i=1,...,\,M$, consider a reduced order model $\bar{\Sigma}_{i},\,i=1,...,\,M$, with state $\bar{x}_{i}\in\mathbb{R}^{\bar{n}_{i}}$, $\bar{n}_{i}\leq n_{i}$, and input $\bar{u}_{i}\in{\mathcal{\bar{U}}}_{i}\subseteq{\mathcal{U}}_{i}$. In a collective form, these systems $\bar{\Sigma}_{i}$ define the overall reduced order model

[TABLE]

where $\bar{x}=(\bar{x}_{1},\ldots,\bar{x}_{M})\in\mathbb{R}^{\bar{n}}$, $\bar{n}=\sum_{i=1}^{M}\bar{n}_{i}$, and $\bar{u}=(\bar{u}_{1},\ldots,\bar{u}_{M})\in{\mathcal{\bar{U}}}=\prod_{i=1}^{M}{\mathcal{\bar{U}}}_{i}\subseteq\mathbb{R}^{m}$.

The reduced order models $\bar{\Sigma}_{i}$ can be defined according to different criteria. In any case, it is required that the stability properties of system $\Sigma$ are inherited by $\bar{\Sigma}$. Moreover, it is assumed that, for each subsystem $i=1,\dots,M$, there exists a state projection $\beta_{i}:\mathbb{R}^{{n}_{i}}\rightarrow\mathbb{R}^{\bar{n}_{i}}$, $i=1,...,M$, that allows to establish a connection between the states $x_{i}(h)$ of the original models and the states of the reduced models $\bar{x}_{i}(h)$. Collectively, we define $\beta=\text{diag}(\beta_{1},...,\beta_{M})$. In principle, the ideal case would be to verify $\bar{x}(h)=\beta x(h)$ for all $h\geq 0$. However, due to model reduction approximations, this ideal assumption must be relaxed; instead, we just ask that $\bar{x}(h)=\beta x(h)$ at least in steady-state conditions. Overall, we require the following standing assumption to be satisfied.

Assumption 2

$A_{\rm\scriptscriptstyle H}$* is Schur stable;* 2. 2.

$\beta_{i}$* is full rank, for each $i=1,\dots,M$;* 3. 3.

letting $\hat{G}_{\rm\scriptscriptstyle L}(z)=\beta(zI-A_{\rm\scriptscriptstyle L})^{-1}B_{\rm\scriptscriptstyle L}$ and $G_{\rm\scriptscriptstyle H}(z)=(zI-A_{\rm\scriptscriptstyle H})^{-1}B_{\rm\scriptscriptstyle H}$, it holds that $\hat{G}_{\rm\scriptscriptstyle L}(1)={G}_{\rm\scriptscriptstyle H}(1)$. $\square$

An algorithm to compute the projections $\beta_{i}$ and the matrices of $\bar{\Sigma}$ can be devised along the lines of Picasso et al. (2016).

3 Design of the hierarchical control structure

In this section the regulators at the two layers of the hierarchical control structure are designed.

3.1 Design of the high level regulator

The high level regulator, designed to work at a lower frequency, is based on the reduced order model (4) sampled with period $N_{\rm\scriptscriptstyle L}$ under the assumption that, $\forall\,k\in{\mathbb{N}}$, the $\bar{u}_{i}^{\prime}s$ are held constant over the interval $h\in[kN_{\rm\scriptscriptstyle L},kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1]$. Denoting by $\bar{u}_{i}^{{[N_{\rm\scriptscriptstyle L}]}}(k)$ these constant values and by $\bar{u}^{{[N_{\rm\scriptscriptstyle L}]}}(k)$ the overall input vector, the reduced order model in the slow timescale is

[TABLE]

where $B_{{\rm{\rm\scriptscriptstyle H}}}^{{[N_{\rm\scriptscriptstyle L}]}}=\sum_{j=0}^{N_{\rm\scriptscriptstyle L}-1}A_{{\rm{\rm\scriptscriptstyle H}}}^{j}B_{{\rm{\rm\scriptscriptstyle H}}}$. In order to feedback a value of $\bar{x}^{[N_{\rm\scriptscriptstyle L}]}$ related to the real state $x$ of the system, the projected value $\beta_{i}x_{i}(kN_{\rm\scriptscriptstyle L})$ must be used, so that the reset

[TABLE]

must be applied. In collective form (6) becomes

[TABLE]

The reset (6) at time $k$ may force $\bar{x}^{[N_{\rm\scriptscriptstyle L}]}(k+1)$ to assume a value different from the one computed based on the dynamics of (5) and the applied input $\bar{u}^{[N_{\rm\scriptscriptstyle L}]}(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 $\bar{w}(k)$, i.e.,

[TABLE]

The size of $\bar{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 $\bar{w}(k)\in\mathcal{W}$, where $\mathcal{W}$ is a compact set containing the origin. The characteristics of $\mathcal{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 Mayne et al. (2005). To this end, we first need to define the “unperturbed” prediction model

[TABLE]

and the control gain matrix $\bar{K}_{\rm\scriptscriptstyle H}$ such that, at the same time

• •

$F_{\rm\scriptscriptstyle H}=A_{{\rm{\rm\scriptscriptstyle H}}}^{N_{\rm\scriptscriptstyle L}}+B_{{\rm{\rm\scriptscriptstyle H}}}^{{[N_{\rm\scriptscriptstyle L}]}}\bar{K}_{\rm\scriptscriptstyle H}$ is Schur stable.

• •

$F_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}=A_{\rm\scriptscriptstyle L}^{N_{\rm\scriptscriptstyle L}}+B_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}\bar{K}_{\rm\scriptscriptstyle H}\beta$ is Schur stable, where $B_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}=\sum_{j=0}^{N_{\rm\scriptscriptstyle L}-1}A_{{\rm{\rm\scriptscriptstyle L}}}^{j}B_{\rm\scriptscriptstyle L}$.

We define $\bar{e}(k)=\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)-\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)$ and we let $\mathcal{Z}$ be a robust positively invariant (RPI) set - minimal, if possible - for the autonomous but perturbed system

[TABLE]

Denoting by $\overrightarrow{\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}}{(t:t+N_{\rm\scriptscriptstyle H}-1)}$ the sequence $\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}(t)$, $\dots$, $\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}(t+N_{\rm\scriptscriptstyle H}-1)$, at each slow time-step $t$ the following optimization problem is solved:

[TABLE]

where

[TABLE]

and $\bar{\mathcal{X}}_{F}$ is a positively invariant terminal set for the unperturbed system (9) controlled with the stabilizing auxiliary control law $\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)=\bar{K}_{\rm\scriptscriptstyle H}\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)$, with $\bar{K}_{\rm\scriptscriptstyle H}\bar{\mathcal{X}}_{F}\subseteq{\bar{\mathcal{U}}\ominus\bar{K}_{\rm\scriptscriptstyle H}\mathcal{Z}}$. Note that it is implicitly assumed that ${\bar{\mathcal{U}}\supset\bar{K}_{\rm\scriptscriptstyle H}\mathcal{Z}}$: this can always made possible by reducing $\bar{K}_{\rm\scriptscriptstyle H}$ and set $\mathcal{W}$ and - in turn - $\mathcal{Z}$; as it will be discuss in the following, the latter can be reduced, for example, by increasing $N_{\rm\scriptscriptstyle L}$. The positive definite weighting matrices ${Q_{\rm\scriptscriptstyle H}}$, ${R_{\rm\scriptscriptstyle H}}$ are free design parameters, while ${P_{\rm\scriptscriptstyle H}}$ is computed as the solution to the Lyapunov equation

[TABLE]

Letting $\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(t|t),\overrightarrow{\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}}{(t:t+N_{\rm\scriptscriptstyle H}-1|t)}$ be the solution to the optimization problem (11), the control variable applied at time $t$ is defined as

[TABLE]

3.2 Design of the low level regulators

Recall that (see again Figure 1) the overall control action has components generated by both the high-level and the low-level controllers, i.e.,

[TABLE]

Indeed, the low level regulators are in charge of computing the local control corrections $\delta u_{i}\in{\mathcal{U}}_{i}\ominus{\bar{\mathcal{U}}}_{i}$ compensating for the effect of the model inaccuracies at the high level expressed by the term $\bar{w}(k)$ in (8). To this end, first define the auxiliary system $\hat{\Sigma}_{i}$ given by

[TABLE]

Note that $\hat{\Sigma}_{i}$ can be simulated in a centralized way in the time interval $[kN_{\rm\scriptscriptstyle L},\,kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L})$ once the high level controller has computed $\bar{u}_{i}^{{[N_{\rm\scriptscriptstyle L}]}}(k)$.

Also denote by $\varDelta\Sigma_{i}$ the model given by the difference of the system (1), with (2), (15), and (16) in the form

[TABLE]

where $\delta x_{i}(h)=x_{i}(h)-\hat{x}_{i}(h)$ , $\delta z_{i}(h)=z_{i}(h)-\hat{z}_{i}(h)$ and $\delta s_{i}(h)=s_{i}(h)-\hat{s}_{i}(h)$.

The difference state $\delta x_{i}$ is available at each time instant $h$ since $x_{i}$ is measurable and $\hat{x}_{i}$ can be computed from the available control $\bar{u}_{i}^{{[N_{\rm\scriptscriptstyle L}]}}(\lfloor h/N_{\rm\scriptscriptstyle L}\rfloor)$. However, the difference dynamical system $\varDelta\Sigma_{i}$ is not yet useful for decentralized prediction since it depends upon the interconnection variables $\delta s_{i}(h)$ that, in turn, depend upon the variables $\delta x_{j}(h)$, $j\neq i$, which are not known in advance in the future prediction horizon. For this reason, we define a decentralized (approximated) dynamical system $\varDelta\hat{\Sigma}_{i}$ by discarding all interconnection inputs and with input $\delta\hat{u}_{i}(h)$ (which will be defined as the result of a suitable optimization problem), i.e.,

[TABLE]

For all $i=1,\dots,M$, the input $\delta u_{i}(h)$ to the real model (17) is computed based on $\delta\hat{u}_{i}(h)$, $\delta x_{i}(h)$, and $\delta\hat{x}_{i}(h)$ using a standard state-feedback policy, i.e.,

[TABLE]

and where $K_{i}$ is designed in such a way that the matrix $F_{\rm\scriptscriptstyle L}=A_{\rm\scriptscriptstyle L}+B_{\rm\scriptscriptstyle L}K$ is Schur stable, being $K=$diag$(K_{1},\dots,K_{M})$.

Assume now to be at time $h=kN_{\rm\scriptscriptstyle L}$ and to have run the high level controller, so that both $\bar{u}_{i}^{{[N_{\rm\scriptscriptstyle L}]}}(k)$ and the predicted value $\bar{x}^{[N_{\rm\scriptscriptstyle L}]}(k+1|k)=A_{\rm\scriptscriptstyle H}^{N_{\rm\scriptscriptstyle L}}\bar{x}^{[N_{\rm\scriptscriptstyle L}]}(k)+B_{\rm\scriptscriptstyle H}^{[N_{\rm\scriptscriptstyle L}]}\bar{u}^{[N_{\rm\scriptscriptstyle L}]}(k)$ are available. Therefore, in order to remove the effect of the mismatch at the high level represented by $\bar{w}(k)$ in (8), the low level controller working in the interval $[kN_{\rm\scriptscriptstyle L},kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-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 (20) can only be formulated in an approximated way with reference to its state $\delta\hat{x}_{i}$ as follows:

[TABLE]

Note however that the fulfillment of (21) does not imply that (20) is satisfied due to the neglected interconnections in (18) which make the term $\bar{w}(k)$ in (5) not identically equal to zero, although it contributes to its reduction.

Letting $\overrightarrow{\delta\hat{u}}_{i}{(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1)}=(\delta\hat{u}_{i}(kN_{\rm\scriptscriptstyle L}),\dots,\delta\hat{u}_{i}(kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1))\in({\mathbb{R}}^{m_{i}})^{N_{\rm\scriptscriptstyle L}-1}$, the low level control action is computed, at time instant $h=kN_{\rm\scriptscriptstyle L}$, based on the solution to the following optimization problem:

[TABLE]

where

[TABLE]

and where a discussion on how to select the set ${\varDelta\hat{\mathcal{U}}}_{i}$ is deferred to Appendix A.1.

Finally, at each (fast) time instant, the control component $\delta u_{i}(kN_{\rm\scriptscriptstyle L}+j)$ is given by $\delta u_{i}(kN_{\rm\scriptscriptstyle L}+j)=\delta\hat{u}_{i}(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})+K_{i}(\delta{x}_{i}(kN_{\rm\scriptscriptstyle L}+j)-\delta\hat{x}_{i}(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L}))$.

4 Properties and algorithm implementation

The recursive feasibility and robust convergence properties of the optimization problems stated at the high and low levels are now established. To this end, define

[TABLE]

where

[TABLE]

Also, let ${\mathcal{A}}(N_{\rm\scriptscriptstyle L})=A_{{\rm{\rm\scriptscriptstyle H}}}^{N_{\rm\scriptscriptstyle L}}\beta-\beta A_{{\rm L}}^{N_{\rm\scriptscriptstyle L}}\in{\mathbb{R}}^{\bar{n}\times n}$, $I_{si}=[0_{\bar{n}_{i},\bar{n}_{1}}\ldots I_{\bar{n}_{i},\bar{n}_{i}}\ldots\\ 0_{\bar{n}_{i},\bar{n}_{M}}]$, {\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{]}, $\rho_{{u}}$ be such that ${\mathcal{{U}}}\subseteq\mathcal{B}_{\rho_{{u}}}(0)$. We now introduce the following technical assumption.

Assumption 3

$\|A_{\rm\scriptscriptstyle L}^{N_{\rm\scriptscriptstyle L}}\|<1$; 2. 2.

for each $i=1,...,M$, letting ${\mathcal{R}_{i}}(N_{\rm\scriptscriptstyle L})=\begin{bmatrix}(A_{{\rm{\rm\scriptscriptstyle L}}}^{ii})^{N_{\rm\scriptscriptstyle L}-1}B_{{\rm{\rm\scriptscriptstyle L}}}^{ii}&\dots&B_{{\rm{\rm\scriptscriptstyle L}}}^{ii}\end{bmatrix}$ be the reachability matrix in $N_{\rm\scriptscriptstyle L}$ steps associated to $(A_{{\rm{\rm\scriptscriptstyle L}}}^{ii},B_{{\rm{\rm\scriptscriptstyle L}}}^{ii})$, matrix

[TABLE]

is full-rank with minimum singular value $\underline{\sigma}_{\mathcal{H}_{i}(N_{\rm\scriptscriptstyle L})}>0$; 3. 3.

letting $\rho_{\bar{u}}$ and $\rho_{\delta\hat{u}_{i}}$ be such that ${\mathcal{\bar{U}}}\subseteq\mathcal{B}_{\rho_{\bar{u}}}(0)$ and ${\Delta\hat{\mathcal{U}}_{i}}\supseteq\mathcal{B}_{\rho_{\delta\hat{u}_{i}}}(0)$, respectively, for any $i=1,...,M$ it holds that

[TABLE] 4. 4.

for each $i=1,\dots,M$

[TABLE] 5. 5.

Defining $\Delta\bar{\mathcal{U}}_{i}=\Delta\mathcal{U}_{i}(N_{\rm\scriptscriptstyle L}-1)$, and $\Delta\mathcal{U}_{i}(j)=\Delta\hat{\mathcal{U}}_{i}\oplus\mathcal{B}_{\rho_{\Delta{u}_{i}}(j)}(0)$ where $\rho_{\Delta{u}_{i}}(j)=\sum_{r=2}^{j}\|K_{i}I_{si}F_{\rm\scriptscriptstyle L}^{j-r}(A_{\rm\scriptscriptstyle L}-A_{\rm\scriptscriptstyle L}^{D})\|\\ \rho_{\delta\hat{x}}(r-1)$ for all $j=2,\dots,N_{\rm\scriptscriptstyle L}$, $\rho_{\Delta{u}_{i}}(j)=0$ for $j=0,1$ we require that

[TABLE]

$\square$**

It is now possible to specify the size of the uncertainty set $\mathcal{W}$ to be considered in the high level design. Specifically, let

[TABLE]

where $\rho_{w}=\sum_{j=2}^{N_{\rm\scriptscriptstyle L}}\|\beta F_{\rm\scriptscriptstyle L}^{N_{\rm\scriptscriptstyle L}-j}(A_{\rm\scriptscriptstyle L}-A_{\rm\scriptscriptstyle L}^{D})\|\rho_{\delta\hat{x}}(j-1)$, $\rho_{\delta\hat{x}}(j)=\sqrt{\sum_{i=1}^{M}\rho^{2}_{\delta\hat{x}_{i}}(j)}$,

[TABLE]

and where $A_{\rm\scriptscriptstyle L}^{D}=$diag$(A_{\rm\scriptscriptstyle L}^{11},\dots,A_{\rm\scriptscriptstyle L}^{MM})$.

The following result can be proved.

Theorem 1

Under Assumption 3, if $x(0)$ is such that the problem (11) is feasible at $k=0$ and, for all $i=1,\dots,M$

[TABLE]

*then

(i) $\bar{w}(k)\in\mathcal{W}$ and problems (11) and (22) are feasible for all $k\geq 0$;

(ii) for all $h\geq 0$*

[TABLE]

(iii) the state of the slow time-scale reduced model $\bar{\Sigma}^{[N_{\rm\scriptscriptstyle L}]}$ enjoys robust convergence properties, i.e.,

[TABLE]

(iv) the state of the large scale model $\Sigma$ enjoys robust convergence properties, i.e., for a computable positive constant $\rho_{x}$

[TABLE]

$\square$**

Theorem 1 establishes two important facts. First, it shows that, if the initial state lies in a suitable set (and Assumption 3 holds), the joint feasibility properties of the two control layers can be guaranteed in a recursive fashion. Secondly, it ensures convergence of the state of the small-scale slow system considered by the higher control layer to a set.

Regarding the main technical Assumption 3, note that it involves quantities, that are all functions of the number of steps $N_{\rm\scriptscriptstyle L}$. It is worth now analysing their dependence upon it. Indeed, the following facts can be proved.

• •

In view of Assumption 2.2, $\beta_{i}$ are full rank and, in view of Assumption 1.3, the pairs $(A_{{\rm{\rm\scriptscriptstyle L}}}^{ii},B_{{\rm{\rm\scriptscriptstyle L}}}^{ii})$ are reachable, so that Assumption 3.2 is fulfilled by taking $N_{\rm\scriptscriptstyle L}$ sufficiently large, i.e., by making the upper layer slower without modifying the rate of the lower layer.

• •

In view of Assumptions 1.2, 2.3, and 2.1, $\kappa(N_{\rm\scriptscriptstyle L})=\|\sum_{j=0}^{N_{\rm\scriptscriptstyle L}-1}A_{\rm\scriptscriptstyle H}^{j}B_{\rm\scriptscriptstyle H}-G_{\rm\scriptscriptstyle H}(1)-(\beta\sum_{j=0}^{N_{\rm\scriptscriptstyle L}-1}A_{\rm\scriptscriptstyle L}^{j}B_{\rm\scriptscriptstyle L}-\hat{G}_{\rm\scriptscriptstyle L}(1))\|$ and $G_{\rm\scriptscriptstyle H}(1)=\sum_{j=0}^{+\infty}A_{\rm\scriptscriptstyle H}^{j}B_{\rm\scriptscriptstyle H}$, $\hat{G}_{\rm\scriptscriptstyle L}(1)=\beta\sum_{j=0}^{+\infty}A_{\rm\scriptscriptstyle L}^{j}B_{\rm\scriptscriptstyle L}$. Therefore $\kappa(N_{\rm\scriptscriptstyle L})\leq\|\sum_{j=N_{\rm\scriptscriptstyle L}}^{+\infty}A_{\rm\scriptscriptstyle H}^{j}B_{\rm\scriptscriptstyle H}\|+\|\beta\sum_{j=N_{\rm\scriptscriptstyle L}}^{+\infty}A_{\rm\scriptscriptstyle L}^{j}B_{\rm\scriptscriptstyle L}\|\leq\|A_{\rm\scriptscriptstyle H}^{N_{\rm\scriptscriptstyle L}}\|\|G_{\rm\scriptscriptstyle H}(1)\|\\ +\|A_{\rm\scriptscriptstyle L}^{N_{\rm\scriptscriptstyle L}}\|\|\beta\|\|G_{\rm\scriptscriptstyle L}(1)\|$, where ${G}_{{\rm{\rm\scriptscriptstyle L}}}(z)=(zI-A_{{\rm{\rm\scriptscriptstyle L}}})^{-1}B_{{\rm{\rm\scriptscriptstyle L}}}$. Therefore $\kappa(N_{\rm\scriptscriptstyle L})\rightarrow 0$ exponentially as $N_{\rm\scriptscriptstyle L}\rightarrow+\infty$. This shows that also Assumption 3.25 can be fulfilled by taking $N_{\rm\scriptscriptstyle L}$ sufficiently large.

• •

Equivalently to Proposition 2.3 in Picasso et al., for any $i=1,...,M$ it can be proved that

[TABLE]

• •

The above considerations also show that, by setting a sufficiently large low-level prediction horizon $N_{\rm\scriptscriptstyle L}$, it is always possible to allow for arbitrarily small input constraint sets $\Delta\hat{\mathcal{U}}_{i}$. This, in turn, allows to obtain an arbitrarily small high-level disturbance set $\mathcal{W}$ and, in turn, a small robust positively invariant set $\mathcal{Z}$ which, eventually, allows to define the input sets in such a way that (27) can be verified.

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 $A_{\rm\scriptscriptstyle H}$, $B_{\rm\scriptscriptstyle H}$, and $\beta_{i}$, $i=1,...,M$, such that Assumption (2) is satisfied;

• •

design of $\bar{K}_{\rm\scriptscriptstyle H}$ such that both $F_{\rm\scriptscriptstyle H}=A_{{\rm{\rm\scriptscriptstyle H}}}^{N_{\rm\scriptscriptstyle L}}+B_{{\rm{\rm\scriptscriptstyle H}}}^{{[N_{\rm\scriptscriptstyle L}]}}\bar{K}_{\rm\scriptscriptstyle H}$ and $F_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}=A_{\rm\scriptscriptstyle L}^{N_{\rm\scriptscriptstyle L}}+B_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}\bar{K}_{\rm\scriptscriptstyle H}\beta$ are Schur stable;

• •

design of $K=diag(K_{1},\dots,K_{M})$ such that $F_{\rm\scriptscriptstyle L}=A_{\rm\scriptscriptstyle L}+B_{\rm\scriptscriptstyle L}K$ is Schur stable;

• •

computation of $\rho_{\delta\hat{u}_{i}}$, $\rho_{\bar{u}_{i}}$ (see the procedure proposed in Appendix A.1) and of the sets ${\bar{\mathcal{U}}}_{i}$, $\Delta\hat{\mathcal{U}}_{i}$;

• •

computation of $\mathcal{W}$ according to (28) and (29);

• •

computation of $\bar{\mathcal{X}}_{\rm\scriptscriptstyle F}$, $\mathcal{Z}$, see Rawlings and Mayne (2009), and $P_{\rm\scriptscriptstyle H}$ with (13).

5 Simulation example

Consider the problem of regulating temperatures of two apartments depicted in Figure 2. The first apartment is constituted by rooms $A_{1}$, $B_{1}$, $C_{1}$, $D_{1}$ and $E_{1}$, while the second one by $A_{2}$, $B_{2}$, $C_{2}$, $D_{2}$ and $E_{2}$. Each apartment is equipped with a radiator supplying heats $q_{i}$, $i=1,\,2$. Heat exchange coefficient between neighbouring rooms of different apartment, i.e., $E_{2}$ and $C_{1}$, is $k_{1}^{t}=1$ W/m2K, the one between adjacent rooms inside each apartment is $k_{2}^{t}=2.5$ W/m2K, and the one between the rooms and the external environment is $k_{e}^{t}=0.5$ W/m2K. The external temperature is $T_{E}=0^{\circ}\mathrm{C}$ and, for simplicity, we neglect solar radiation. Furthermore, the height of the walls is $H=4$ m. Air density and heat capacity are $\rho=1.225$ kg/m3 and $c=1005$ J/kgK, respectively. The overall model is made by dynamic energy balance equations of each room. The variables $q_{1},\,q_{2}$ are expressed in Watts, while all the temperature variables are expressed in ${}^{\circ}\mathrm{C}$. The considered equilibrium point is: $\bar{q}=(\bar{q}_{1},\bar{q}_{2})=(354.2,320.8)$, with $\bar{T}=(\bar{T}_{A_{1}},\bar{T}_{B_{1}},\bar{T}_{C_{1}},\bar{T}_{D_{1}},\bar{T}_{E_{1}},\bar{T}_{A_{2}},\bar{T}_{B_{2}},\bar{T}_{C_{2}},\bar{T}_{D_{2}},\bar{T}_{E_{2}})=(19.6,20.3,20.2,\\ 21.7,18.2,17.2,21.2,21.7,19.6,19.4)$. Let $\delta T_{j_{i}}=T_{j_{i}}-\bar{T}_{j_{i}}$ and $\delta q_{i}=q_{i}-\bar{q}_{i}$, for $j=A,\,B,\,C,\,D$ and $i=1,\,2$. In this way, $x_{i}=(\delta T_{A_{i}},\,\delta T_{B_{i}},\,\delta T_{C_{i}},\,\delta T_{D_{i}},\,\delta T_{E_{i}})$ and $u_{i}=\delta q_{i}$ are the state and input variables of the $i$-th subsystem , i.e., $n_{i}=5$ and $m_{i}=1$, with $i=1,\,2$.

The control variables are limited, i.e., $-50\leqslant u_{1},u_{2}\leqslant 50$.

The two subsystems’ continuous-time models have been sampled using the algorithm described in Farina et al. (2013) with $\Delta t=90\rm s$ to obtain their discrete-time counterpart in the fast time scales. The eigenvalues of the first subsystem are $\{0.73,\,0.97,\,0.9,\,0.85,\,0.88\}$, and the eigenvalues of the second one are $\{0.97,\,0.76,\,0.82,\,0.91,\,0.87\}$. Then the procedure described in the Section (1) has been used to compute the discrete-time reduced order model, with $A_{\rm\scriptscriptstyle H}=\rm diag(0.97,0.97)$, i.e., $\bar{n}=2$, as well as the transformation matrices $\beta_{1}=\begin{bmatrix}0&0&-1&0&0\end{bmatrix}$ and $\beta_{2}=\begin{bmatrix}0&0&0&0&-1\end{bmatrix}$. The plant model in the slow time scale has been constructed with $N_{\rm\scriptscriptstyle L}=20$.

Tube-based Robust MPC has been designed at the high level according to the algorithm described in Mayne et al. (2005) with prediction horizon, $N_{\rm\scriptscriptstyle H}=10$, state and input penalty, $Q_{\rm\scriptscriptstyle H}=I_{\bar{n}}$ and $R_{\rm\scriptscriptstyle H}=0.1I_{\bar{m}}$.

At the low level, the finite-horizon optimization algorithms described in (22) have been implemented with $Q_{1}=I_{n_{1}}$, $Q_{2}=I_{n_{2}}$, $R_{1}=R_{2}=10$.

The hierarchical control scheme has been simulated starting from $x(0)=(x_{1}(0),x_{2}(0))=(-2,\dots,-2)$ and $\bar{x}(0)=\beta x(0)$. The transients of the state and control variables are reported in Figures 3- 5. These results show the effectiveness of the proposed algorithm.

6 Conclusions

A two layer control scheme for systems made by interconnected subsystems has been presented. Its performance has been tested in simulation, and its properties of recursive feasibility and convergence to a set have been established. Current research is focusing on the extension of the analysis to guarantee convergence to the origin and to deal with tracking problems, as well as the application of the proposed approach to other systems with large dimensions.

APPENDIX

A.1 Computation of the input constraint sets

In the scheme proposed in this paper, the dimensions of the input constraint sets $\bar{\mathcal{U}}_{i}$ and $\Delta\hat{\mathcal{U}}_{i}$ are key tuning knobs. They must be tuned in order to satisfy, at the same time, the inequalities (25), for all $i=1,\dots,M$ and (27). 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 $\Delta\hat{\mathcal{U}}_{i}=\mathcal{B}_{\rho_{\delta\hat{u}_{i}}}(0)$ and $\bar{\mathcal{U}}_{i}=\mathcal{B}_{\rho_{\bar{u}_{i}}}(0)$. Under this assumption, the tuning knobs are the vectors $\overrightarrow{\rho}_{\delta\hat{u}}=(\rho_{\delta\hat{u}_{1}},\dots,\rho_{\delta\hat{u}_{M}})$ and $\overrightarrow{\rho}_{\bar{u}}=(\rho_{\bar{u}_{1}},\dots,\rho_{\bar{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 (25), to be verified for all $i=1,\dots,M$. Here the constant $\rho_{\bar{u}}$ appears, defined in such a way that $\bar{\mathcal{U}}=\prod_{i=1}^{M}\mathcal{B}_{\rho_{\bar{u}_{i}}}(0)\subseteq\mathcal{B}_{\rho_{\bar{u}}}(0)$. We can define, for example, $\rho_{\bar{u}}=\sqrt{\sum_{i=1}^{M}\rho_{\bar{u}_{i}}^{2}}\leq\sum_{i=1}^{M}\rho_{\bar{u}_{i}}$. Therefore, to fulfill (25) it is sufficient to verify the following matrix inequality

[TABLE]

where $\mathds{1}_{M\times M}$ is the $M\times M$ matrix whose entries are all equal to $1$. The second main inclusion to be fulfilled is (27), which is verified if, for all $i=1,\dots,M$,

[TABLE]

By definition, $\Delta\bar{\mathcal{U}}_{i}=\Delta\hat{\mathcal{U}}_{i}\oplus\mathcal{B}_{\rho_{\Delta u_{i}(N_{\rm\scriptscriptstyle L}-1)}}(0)$, where

[TABLE]

where $\lambda_{ij}=\sum_{r=2}^{N_{\rm\scriptscriptstyle L}-1}\|K_{i}I_{si}F_{\rm\scriptscriptstyle L}^{N_{\rm\scriptscriptstyle L}-r-1}(A_{\rm\scriptscriptstyle L}-A_{\rm\scriptscriptstyle L}^{D})\|\sum_{k=1}^{r-1}\|(A_{\rm\scriptscriptstyle L}^{jj})^{r-1-k}\\ B_{\rm\scriptscriptstyle L}^{jj}\|$. This implies that $\Delta\bar{\mathcal{U}}_{i}=\mathcal{B}_{\rho_{\delta\hat{u}_{i}}+\sum_{j=1}^{M}\lambda_{ij}\rho_{\delta\hat{u}_{j}}}(0)$. Therefore, to verify (33) it is sufficient to enforce the constraint

[TABLE]

where ${\Lambda}$ is the $M\times M$ matrix whose entries are $\lambda_{ij}$, $i,j=1,\ldots,M$, while $\overrightarrow{\rho}_{{u}}=(\rho_{u_{1}},\dots,\rho_{u_{M}})$, where $\mathcal{B}_{\rho_{i}}(0)\subseteq\mathcal{U}_{i}$ for all $i=1,\dots,M$. Eventually, a suitable choice of $\overrightarrow{\rho}_{\delta\hat{u}}$ and $\overrightarrow{\rho}_{\bar{u}}$ is obtained as the solution to the following linear programming problem:

[TABLE]

where $J_{\rho}=\gamma_{1}\mathds{1}_{1\times M}\overrightarrow{\rho}_{\delta\hat{u}}+\gamma_{2}\mathds{1}_{1\times M}\overrightarrow{\rho}_{\bar{u}}$, where $\gamma_{1}$, $\gamma_{2}$ are arbitrary positive weighting constants.

A.2 Proof of Theorem 1

The proof of Theorem 1 lies on the intermediary results stated below.

Proposition 1

A) Under Assumption 3 and if $\bar{x}^{[N_{\rm\scriptscriptstyle L}]}(k)=\beta x(kN_{\rm\scriptscriptstyle L})$, then for any initial condition $\hat{x}(kN_{\rm\scriptscriptstyle L})=x(kN_{\rm\scriptscriptstyle L})$ such that, for all $i=1,\dots,M$

[TABLE]

*and for any $\bar{u}^{{[N_{\rm\scriptscriptstyle L}]}}\in{\mathcal{\bar{U}}}$ there exists a feasible sequence

$\overrightarrow{\delta\hat{u}}_{i}{(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1|kN_{\rm\scriptscriptstyle L})}\in{\varDelta\hat{\mathcal{U}}_{i}}^{N_{\rm\scriptscriptstyle L}}$ such that the terminal constraint (21) is satisfied.

B) if $x(kN_{\rm\scriptscriptstyle L})$ satisfies condition (37), $\|A_{\rm\scriptscriptstyle L}^{N_{\rm\scriptscriptstyle L}}\|<1$, and, for all $i=1,\dots,M$, (26) is verified, then recursive feasibility of the terminal constraint (21) is guaranteed.*

Proof of Proposition 1

A) Consider the constraint (21) and first note that, since $\delta\hat{x}_{i}(kN_{\rm\scriptscriptstyle L})=0$,

[TABLE]

Moreover, in view of (5)

[TABLE]

Analogously, from (16) written in collective form

[TABLE]

In view of (38), (39), (40), and the definition ${\mathcal{A}}(N_{\rm\scriptscriptstyle L})$, ${\mathcal{B}}(N_{\rm\scriptscriptstyle L})$, and $I_{si}$, the constraint (21) can be written as

[TABLE]

From this expression, the definitions of $\underline{\sigma}_{\mathcal{H}_{i}(N_{\rm\scriptscriptstyle L})}$, $\rho_{\bar{u}}$, $\rho_{\delta\hat{u}_{i}}$, and in view of (24), it can be concluded that a feasible sequence $\overrightarrow{\delta\hat{u}}_{i}{(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1|kN_{\rm\scriptscriptstyle L})}$ can be computed provided that

[TABLE]

from which the result follows.

B) From (3) it holds that

[TABLE]

Therefore

[TABLE]

and, in view of (37)

[TABLE]

for all $i=1,2,\dots,M$. From this expression and the definition of $\chi_{i}(kN_{\rm\scriptscriptstyle L})$ through (26) it turns out that

[TABLE]

for all $i=1,2,\dots,M$ and the result follows. $\square$

Proposition 2

*If Problem (22) is feasible, then

A) For all $k\geq 0$*

[TABLE]

B) For all $k\geq 0$,

[TABLE]

Also it holds that

[TABLE]

for all $j=0,\dots,N_{\rm\scriptscriptstyle L}-1$. $\square$

Proof of Proposition 2

A) Defining the collective vectors $\hat{x}=(\hat{x}_{1},\dots,\hat{x}_{M})$, $\delta x=(\delta x_{1},\dots,\delta x_{M})$, $\delta\hat{x}=(\delta\hat{x}_{1},\dots,\delta\hat{x}_{M})$, and $\varepsilon(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})=\delta x(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})-\delta\hat{x}(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})$, we have that $\bar{w}(k)=\beta x(kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L})-\bar{x}^{[N_{\rm\scriptscriptstyle L}]}(k+1|k)=\beta\hat{x}(kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L})+\beta\delta x(kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L})-\bar{x}^{[N_{\rm\scriptscriptstyle L}]}(k+1|k)=(\beta\hat{x}(kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L})+\beta\delta\hat{x}(kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L})-\bar{x}^{[N_{\rm\scriptscriptstyle L}]}(k+1|k))+\beta\varepsilon(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})=\beta\varepsilon(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})$. The latter equality holds in view of the fact that Problem (22) is feasible, and therefore equality (21) is verified. From (17), (18), (19), we collectively have that

[TABLE]

In view of the fact that $\varepsilon(kN_{\rm\scriptscriptstyle L}|kN_{\rm\scriptscriptstyle L})=\delta\hat{x}(kN_{\rm\scriptscriptstyle L}|kN_{\rm\scriptscriptstyle L})=0$, then $\bar{w}(k)=\beta\sum_{j=2}^{N_{\rm\scriptscriptstyle L}}F_{\rm\scriptscriptstyle L}^{N_{\rm\scriptscriptstyle L}-j}(A_{\rm\scriptscriptstyle L}-A_{\rm\scriptscriptstyle L}^{D})\delta\hat{x}(kN_{\rm\scriptscriptstyle L}+j-1|kN_{\rm\scriptscriptstyle L})$. From this it follows that

[TABLE]

Since $\delta\hat{u}_{i}$ are bounded for all $i=1,\dots,M$, i.e., scalar $\rho_{\delta\hat{u}_{i}}$ are defined such that $\delta\hat{u}_{i}\in\mathcal{B}_{\rho_{\delta\hat{u}_{i}}}(0)$. In view of this, we compute that $\|\delta\hat{x}(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})\|\leq\rho_{\delta\hat{x}}(j)$, where $\rho_{\delta\hat{x}_{i}}(j)$ is defined in (29). Therefore, $\delta\hat{x}(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})$ are bounded, for all $j=1,\dots,N_{\rm\scriptscriptstyle L}-1$ and more specifically we get that $\|\delta\hat{x}(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})\|\leq\sqrt{\sum_{i=1}^{M}\rho^{2}_{\delta\hat{x}_{i}}(j)}$. Therefore one has (47) for all $k\geq 0$.

B) From (50) we have that $\varepsilon(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})=\sum_{r=2}^{j}F_{\rm\scriptscriptstyle L}^{j-r}(A_{\rm\scriptscriptstyle L}-A_{\rm\scriptscriptstyle L}^{D})\delta\hat{x}(kN_{\rm\scriptscriptstyle L}+r-1|kN_{\rm\scriptscriptstyle L})$ and therefore $\delta u_{i}(kN_{\rm\scriptscriptstyle L}+j)-\delta\hat{u}_{i}(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})=K_{i}\varepsilon_{i}(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})=K_{i}I_{si}\sum_{r=2}^{j}F_{\rm\scriptscriptstyle L}^{j-r}(A_{\rm\scriptscriptstyle L}-A_{\rm\scriptscriptstyle L}^{D})\delta\hat{x}(kN_{\rm\scriptscriptstyle L}+r-1|kN_{\rm\scriptscriptstyle L})$. From this it follows that $\delta u_{i}(kN_{\rm\scriptscriptstyle L}+j)\in\delta\hat{u}_{i}(kN_{\rm\scriptscriptstyle L}+j|kN_{\rm\scriptscriptstyle L})\oplus\mathcal{B}_{\rho_{\Delta{u}_{i}}(j)}(0)$ and therefore $\delta u_{i}(kN_{\rm\scriptscriptstyle L}+j)\in\Delta\mathcal{U}_{i}(j)$. In view of the monotonicity property $\rho_{\Delta{u}_{i}}(j+1)\geq\rho_{\Delta{u}_{i}}(j)$ for all $j$, it holds that $\mathcal{B}_{\rho_{\Delta{u}_{i}}(j+1)}(0)\supseteq\mathcal{B}_{\rho_{\Delta{u}_{i}}(j)}(0)$, which implies (49). $\square$

Proof of Theorem 1

(i) If $\|x(0)\|\leq\lambda_{i}(N_{\rm\scriptscriptstyle L})$ and recalling that Assumption 3 holds, from Proposition 1, recursive feasibility of the optimization problems (22) is guaranteed, i.e., that there exists, for all $k\geq 0$, a feasible sequence $\overrightarrow{\delta\hat{u}}_{i}{(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1|kN_{\rm\scriptscriptstyle L})}\in{\varDelta\hat{\mathcal{U}}_{i}}^{N_{\rm\scriptscriptstyle L}}$ such that the terminal constraint (21) is satisfied.

Also, from Proposition 2.A, it is proved that $\bar{w}(k)\in\mathcal{W}$ for all $k\geq 0$, which allows to apply the recursive feasibility arguments of Mayne et al. (2005), proving that also (11) enjoys recursive feasibility properties.

(ii) It is now possible to conclude that, in view of the feasibility of (11), $\bar{u}^{[N_{\rm\scriptscriptstyle L}]}(k)\in\bar{\mathcal{U}}$; also, from Proposition 2.B it follows that $\delta u_{i}(kN_{\rm\scriptscriptstyle L}+j)\in\Delta\bar{\mathcal{U}}_{i}$ for all $k\geq 0$, $j=0,\dots,N_{\rm\scriptscriptstyle L}-1$, and $i=1,\dots,M$. From this, under (27), the inclusion (30) can also be proved.

(iii) We apply the results in Mayne et al. (2005), which guarantee robust convergence properties. In other words, we have that $\bar{x}^{[N_{\rm\scriptscriptstyle L}],o}(k)\rightarrow 0$ as $k\rightarrow+\infty$, and that $\bar{x}^{[N_{\rm\scriptscriptstyle L}]}(k)$ is asymptotically driven to lie in the robust positively invariant set $\mathcal{Z}$.

(iv) To show robust convergence of the global system state, from (3) we obtain that

[TABLE]

Denoting $\mathcal{B}_{L,N_{\rm\scriptscriptstyle L}}^{C}=\begin{bmatrix}A_{\rm\scriptscriptstyle L}^{N-1}B_{\rm\scriptscriptstyle L}&\dots&B_{\rm\scriptscriptstyle L}\end{bmatrix}$, we obtain that

[TABLE]

Also, recall that $\overrightarrow{\delta{u}}{(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1)}=\overrightarrow{\delta\hat{u}}{(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1)}+{\rm{diag}}(K,\dots,K)\overrightarrow{\varepsilon}{(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1)}$ and that, by defining

[TABLE]

[TABLE]

where $A_{\rm\scriptscriptstyle L}^{C}=A_{\rm\scriptscriptstyle L}-A_{\rm\scriptscriptstyle L}^{D}$ and

[TABLE]

Recalling that $\bar{x}^{[N_{\rm\scriptscriptstyle L}]}(k)=\beta x(kN_{\rm\scriptscriptstyle L})$ and that

[TABLE]

we can rewrite (52) as

[TABLE]

Recall that $\bar{x}^{[N_{\rm\scriptscriptstyle L}],o}(k),\,\bar{u}^{[N_{\rm\scriptscriptstyle L}],o}(k)\rightarrow 0$ as $k\rightarrow+\infty$. Also, we compute that

[TABLE]

Based on this, we define $\kappa_{\delta u}=\|\mathcal{B}_{L,N_{\rm\scriptscriptstyle L}}^{C}(I+{\rm{diag}}(K_{i},\dots,K_{i})\mathcal{F}_{N_{\rm\scriptscriptstyle L}}\\ \cdot{\rm{diag}}(A_{\rm\scriptscriptstyle L}^{C},\dots,A_{\rm\scriptscriptstyle L}^{C})\mathcal{B}_{\rm\scriptscriptstyle L})\|$ and we write (54) as

[TABLE]

where \|w_{\rm\scriptscriptstyle L}(k)\|\leq\kappa_{\delta u}\sqrt{N_{\rm\scriptscriptstyle L}}\max_{h\in\{kN_{\rm\scriptscriptstyle L},\dots,(k+1)N_{\rm\scriptscriptstyle L}-1\}}\|\delta\hat{u}(h)\|$$\leq\kappa_{\delta u}\sqrt{N_{\rm\scriptscriptstyle L}}\sqrt{\sum_{i=1}^{M}\rho^{2}_{\delta\hat{u}_{i}}}. Therefore, since $F_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}=A_{\rm\scriptscriptstyle L}^{N_{\rm\scriptscriptstyle L}}+B_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}\bar{K}_{\rm\scriptscriptstyle H}\beta$ is Schur stable, then the asymptotic result follows, where $\rho_{x}=\kappa_{\delta u}\sqrt{N_{\rm\scriptscriptstyle L}}\sqrt{\sum_{i=1}^{M}\rho^{2}_{\delta\hat{u}_{i}}}$. $\square$