A hierarchical multirate MPC scheme for interconnected systems - $\textit{extended version}$

TL;DR

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

Effective coordination of hierarchical controllers demonstrated in simulation.

Potential for improved control performance in interconnected systems.

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.

Tables

Table 1. Table 1: On-line computation time comparison

Approach			Optimization activated at	Computation time (s)
Proposed one with $N_L=28$	HL regulator		$h=k{N}_L$	1.1
	LL regulator	$1^{st}$ one	$h=k{N}_L$	0.94
		$2^{nd}$ one	$h=k{N}_L$	0.45
Proposed one with $N_L=84$	HL regulator		$h=k{N}_L$	1.07
	LL regulator	$1^{st}$ one	$h=k{N}_L$	2.16
		$2^{nd}$ one	$h=k{N}_L$	3.64
Centralized stabilizing MPC			each fast time instant $h$	21.9

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)\\ y_{i}(h)&=&C_{\rm\scriptscriptstyle L}^{ii}x_{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:\,\left\{\begin{array}[]{lcl}x(h+1)&=&A_{\rm\scriptscriptstyle L}x(h)+B_{\rm\scriptscriptstyle L}u(h)\\ y(h)&=&C_{\rm\scriptscriptstyle L}x(h)\end{array}\right.$$

$$\bar{\Sigma}:\,\left\{\begin{array}[]{clc}\bar{x}(h+1)&=&A_{\rm\scriptscriptstyle H}\bar{x}(h)+B_{\rm\scriptscriptstyle H}\bar{u}(h)\\ \bar{y}(h)&=&C_{\rm\scriptscriptstyle H}\bar{x}(h)\end{array}\right.$$

$$\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{x}_{\rm\scriptscriptstyle S}\oplus\bar{\mathcal{X}}_{\rm\scriptscriptstyle F}\\ \bullet\,\,\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)\in{\bar{u}_{\rm\scriptscriptstyle S}\oplus\bar{\mathcal{U}_{\rm\scriptscriptstyle S}}\ominus\bar{K}_{\rm\scriptscriptstyle H}\mathcal{Z}},\,k=t,\dots,t+N_{\rm\scriptscriptstyle H}-1,\end{array}$$

$$\begin{array}[]{cll}J_{{\rm\scriptscriptstyle H}}&=\sum_{k=t}^{t+N_{\rm\scriptscriptstyle H}-1}\|\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)-\bar{x}_{\rm\scriptscriptstyle S}\|_{Q_{\rm\scriptscriptstyle H}}^{2}+\|\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)-\bar{u}_{\rm\scriptscriptstyle S}\|_{R_{\rm\scriptscriptstyle H}}^{2}\\ &+\|\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(t+N_{\rm\scriptscriptstyle H})-\bar{x}_{\rm\scriptscriptstyle S}\|_{P_{\rm\scriptscriptstyle H}}^{2}\end{array}$$

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

$$\varDelta\hat{\Sigma}_{i}^{[\zeta_{i}]}:\ \left\{\begin{array}[]{c}\begin{array}[]{l}\delta\hat{x}_{i}^{[\zeta_{i}]}(l_{i}+1)=({A_{\rm\scriptscriptstyle L}^{ii}})^{\zeta_{i}}\delta\hat{x}_{i}^{[\zeta_{i}]}(l_{i})+{B_{\rm\scriptscriptstyle L}^{ii}}^{[\zeta_{i}]}\delta\hat{u}^{[\zeta_{i}]}_{i}(l_{i})\end{array}\\ \delta\hat{x}_{i}^{[\zeta_{i}]}(kN_{i})=\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+1)N_{\mathrm{L}})={\bar{x}}_i^{[N_{\mathrm{L}}]}(k+1|k)$$

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

$$\begin{array}[]{rl}&\beta_{i}\delta\hat{x}_{i}^{[\zeta_{i}]}((k+1)N_{i})=\beta_{i}\delta\hat{x}_{i}((k+1)N_{\rm\scriptscriptstyle L})=\\[2.84544pt] &=\bar{x}_{i}^{[N_{\rm\scriptscriptstyle L}]}(k+1|k)-\beta_{i}\hat{x}_{i}((k+1)N_{\rm\scriptscriptstyle L})\end{array}$$

$$\begin{array}[]{l}\min_{\overrightarrow{\delta\hat{u}^{[\zeta_{i}]}_{i}}{(kN_{i}:(k+1)N_{i}-1)}}J_{{\rm\scriptscriptstyle L}}\left(\overrightarrow{\delta\hat{u}^{[\zeta_{i}]}_{i}}{(kN_{i}:(k+1)N_{i}-1)}\right)\\ \mbox{subject to:}\\ \bullet\,\,\mbox{the dynamics}~{}\eqref{Eqn:LL dndynopt}\\ \bullet\,\,\mbox{the terminal constraint}~{}\eqref{eq:terminaln}\\ \bullet\,\,\delta\hat{u}^{[\zeta_{i}]}_{i}(kN_{i}+l_{i})\in{\Delta\hat{\mathcal{U}}_{i}},\,l_{i}=0,\dots,N_{i}-1,\end{array}$$

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

$$\begin{array}[]{rl}&\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}))\end{array}$$

$$\delta {\hat{u}}_i(k{N}_{\mathrm{L}}+j|k{N}_{\mathrm{L}})=\delta {\hat{u}}_i^{[{\zeta }_i]}(k{N}_i+\lfloor j/{\zeta }_i\rfloor |k{N}_{\mathrm{L}})$$

$$\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}}$$

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 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 $\Sigma$ is assumed to be composed by $M$ interconnected subsystems $\Sigma_{1},...,\Sigma_{\rm\scriptscriptstyle M}$. A reduced order model $\bar{\Sigma}$ is first computed and a centralized MPC regulator $R_{\rm\scriptscriptstyle H}$ working at a slow rate is designed to consider the long-term behavior of the system and to compute the control variables $\bar{u}_{i}$, $i=1,...,M$. Then, faster local regulators $R_{{\rm\scriptscriptstyle L}i}$, $i=1,...,M$, are designed for each subsystem $\Sigma_{i}$ to compute the control contributions $\delta u_{i}$ compensating for the inaccuracies of the high layer design due to the mismatch between $\Sigma$ and $\bar{\Sigma}$. 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 $\Sigma_{i}$ with joint output constraints were considered, no multirate implementations were allowed, and the problem was to coordinate the $\Sigma_{i}^{\prime}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 $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_{\rm\scriptscriptstyle M})=[\,z_{1}^{T}\ z_{2}^{T}\ \cdots\ z_{\rm\scriptscriptstyle M}^{T}\,]^{T}\in{\mathbb{R}}^{q}$, where $q=\sum_{i=1}^{M}q_{i}$. We use $\mathcal{N}_{+}$ to denote the set of positive integer numbers. The symbols $\oplus$/$\ominus$ denote the Minkowski sum/difference. We denote by $\|\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 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 $\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{\mathbb{R}}^{m_{i}}$, and $y_{i}\in{\mathbb{R}}^{p_{i}}$ 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 $\Sigma_{i}$ are represented by the coupling input and output vectors $s_{i}\in{\mathbb{R}}^{p_{si}}$ and $z_{i}\in{\mathbb{R}}^{p_{zi}}$, respectively, where

[TABLE]

with $L_{ii}=0$, $i=1,...,M$.

From (1) and (2), the overall dynamical model $\Sigma$ is

[TABLE]

where $x=(x_{1},\dots,x_{\rm\scriptscriptstyle M})\in\mathbb{R}^{n}$, $n=\sum_{i=1}^{\rm\scriptscriptstyle M}n_{i}$, $u=(u_{1},\dots,u_{\rm\scriptscriptstyle M})\in\mathbb{R}^{m}$, $m=\sum_{i=1}^{M}m_{i}$, and $y=(y_{1},\dots,y_{\rm\scriptscriptstyle M})\in\mathbb{R}^{p}$, $p=\sum_{i=1}^{M}p_{i}$. The diagonal blocks of $A_{\rm\scriptscriptstyle L}$ are state transition matrices $A_{\rm\scriptscriptstyle L}^{ii}$, whereas the coupling terms among the $\Sigma_{i}^{\prime}s$ 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 and output matrices are $B_{\rm\scriptscriptstyle L}=$diag$(B_{\rm\scriptscriptstyle L}^{11},...,\,B_{\rm\scriptscriptstyle L}^{\rm\scriptscriptstyle MM})$ and $C_{\rm\scriptscriptstyle L}=$diag$(C_{\rm\scriptscriptstyle L}^{11},...,\,C_{\rm\scriptscriptstyle L}^{\rm\scriptscriptstyle MM})$, respectively.

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.

$m=p$* and the system matrix $S_{\Sigma}=\begin{bmatrix}I-A_{\rm\scriptscriptstyle L}&-B_{\rm\scriptscriptstyle L}\\ C_{\rm\scriptscriptstyle L}&0\end{bmatrix}$ has full rank $n+m$;* 4. 4.

$B_{\rm\scriptscriptstyle L}$* and $C_{\rm\scriptscriptstyle L}$ are full-rank matrices;* 5. 5.

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

2.2 Reduced order models

For each subsystem $\Sigma_{i},\,i=1,...,\,M$, we define 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\mathbb{R}^{m_{i}}$. In a collective form, these systems $\bar{\Sigma}_{i}$ define the overall reduced order model

[TABLE]

where $\bar{x}=(\bar{x}_{1},\dots,\bar{x}_{\rm\scriptscriptstyle M})\in\mathbb{R}^{\bar{n}}$, $\bar{n}=\sum_{i=1}^{M}\bar{n}_{i}$, $\bar{u}=(\bar{u}_{1},\dots,\bar{u}_{\rm\scriptscriptstyle M})\in\mathbb{R}^{m}$, and $\bar{y}\in\mathbb{R}^{p}$. The reduced order models $\bar{\Sigma}_{i}$ can be defined according to different criteria. First, it is necessary that the stability properties of system $\Sigma$ are inherited by $\bar{\Sigma}$. Moreover, we assume 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_{\rm\scriptscriptstyle M})$. In principle, the ideal case would be to verify that, if $\bar{u}(h)={u}(h)$, then $\bar{x}(h)=\beta x(h)$ and $\bar{y}(h)=y(h)$ for all $h\geq 0$ for suitable initial conditions. However, due to model reduction approximations, this ideal assumption must be relaxed; instead, we just ask that $\bar{x}=\beta x$ and that $\bar{y}=y$ in steady-state conditions. We also require that, while matrix $B_{\rm\scriptscriptstyle H}$ can be full, the output matrix $C_{\rm\scriptscriptstyle H}$ preserves the block-diagonal form of $C_{\rm\scriptscriptstyle L}$, i.e., that $C_{\rm\scriptscriptstyle H}=$diag$(C_{\rm\scriptscriptstyle H}^{11},\dots,C_{\rm\scriptscriptstyle H}^{\rm\scriptscriptstyle MM})$, where $C_{\rm\scriptscriptstyle H}^{ii}\in\mathbb{R}^{p_{i}\times\bar{n}_{i}}$ for all $i=1,\dots,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

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

$\beta_{i}$* is full rank and is such that $C_{\rm\scriptscriptstyle H}^{ii}\beta_{i}=C_{\rm\scriptscriptstyle L}^{ii}$, 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 ${G}_{\rm\scriptscriptstyle H}(1)$ is full rank and $\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}$ 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 $\Sigma$ to the reference value ${y}_{\rm\scriptscriptstyle S}$, while respecting suitable input constraints.

Thanks to Assumption 1.(3), it is possible to compute the reference pair $({x}_{\rm\scriptscriptstyle S},{u}_{\rm\scriptscriptstyle S})$ corresponding to ${y}_{\rm\scriptscriptstyle S}$, i.e., such that ${x}_{\rm\scriptscriptstyle S}=A_{\rm\scriptscriptstyle L}{x}_{\rm\scriptscriptstyle S}+B_{\rm\scriptscriptstyle L}{u}_{\rm\scriptscriptstyle S}$ and $C_{\rm\scriptscriptstyle L}{x}_{\rm\scriptscriptstyle S}={y}_{\rm\scriptscriptstyle S}$. Correspondingly, we define $\bar{u}_{\rm\scriptscriptstyle S}={u}_{\rm\scriptscriptstyle S}$ as the steady-state input reference for the reduced-order system $\bar{\Sigma}$, and the corresponding reference steady-state value as $\bar{x}_{\rm\scriptscriptstyle S}=G_{\rm\scriptscriptstyle H}(1)\bar{u}_{\rm\scriptscriptstyle S}=\beta{x}_{\rm\scriptscriptstyle S}$ by Assumption 2.(3).

At the same time, we aim to enforce input constraints of type $u_{i}(h)\in\bar{u}_{{\rm\scriptscriptstyle S},i}\oplus\mathcal{U}_{{\rm\scriptscriptstyle S},i}$ for all $i=1,\dots,M$, where $\bar{u}_{{\rm\scriptscriptstyle S},i}$ is the $i$-th vector component of $\bar{u}_{{\rm\scriptscriptstyle S}}$ and $\mathcal{U}_{{\rm\scriptscriptstyle S},i}$, $i=1,\dots,M$, are closed and convex sets containing the origin in their interiors. Note that, if the reference $\bar{u}_{{\rm\scriptscriptstyle S},i}$ changes, also the set $\mathcal{U}_{{\rm\scriptscriptstyle S},i}$ may vary to enforce absolute input limitations or saturations. At a collective level, the required constraints are $u(h)\in\bar{u}_{\rm\scriptscriptstyle S}\oplus\mathcal{U}_{\rm\scriptscriptstyle S}$, where $\mathcal{U}_{\rm\scriptscriptstyle S}=\prod_{i=1}^{M}\mathcal{U}_{{\rm\scriptscriptstyle S},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 $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},(k+1)N_{\rm\scriptscriptstyle L}-1]$. Therefore, the sampling time of the high-level model is $N_{\rm\scriptscriptstyle L}$ times larger than the basic sampling time, used in the model (1). Denoting by $\bar{u}_{i}^{{[N_{\rm\scriptscriptstyle L}]}}(k)$ the constant values of $\bar{u}_{i}$ in the long sampling period $k$ 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}}}$. To enforce the input constraints specified above, we will require that $\bar{u}^{{[N_{\rm\scriptscriptstyle L}]}}(k)\in\bar{u}_{\rm\scriptscriptstyle S}\oplus\bar{\mathcal{U}}_{\rm\scriptscriptstyle S}$, where $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S}=\prod_{i=1}^{M}\bar{\mathcal{U}}_{{\rm\scriptscriptstyle S},i}$, $\bar{\mathcal{U}}_{{\rm\scriptscriptstyle S},i}\subset{\mathcal{U}}_{{\rm\scriptscriptstyle S},i}$ for each $i=1,\dots,M$, and where the properties of sets $\bar{\mathcal{U}}_{{\rm\scriptscriptstyle S},i}$ are specified later in the paper.

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})$ is used, so that the reset

[TABLE]

must be applied for all $i=1,\dots,M$. 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 (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 [27]. 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}]}}(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]

The prediction horizon for the high-level MPC consists of $N_{\rm\scriptscriptstyle H}$ slow time steps. 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 the input constrained set has been properly tightened in accordance with the used tube-based control approach, and where

[TABLE]

The set $\bar{\mathcal{X}}_{\rm\scriptscriptstyle F}$ is a positively invariant terminal set for the unperturbed system (9) controlled with the stabilizing control law $\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)=\bar{K}_{\rm\scriptscriptstyle H}\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)$, satisfying $\bar{K}_{\rm\scriptscriptstyle H}\bar{\mathcal{X}}_{\rm\scriptscriptstyle F}\subseteq{\bar{\mathcal{U}_{\rm\scriptscriptstyle S}}\ominus\bar{K}_{\rm\scriptscriptstyle H}\mathcal{Z}}$. In view of this, $\bar{x}_{\rm\scriptscriptstyle S}\oplus\bar{\mathcal{X}}_{\rm\scriptscriptstyle F}$ results positively-invariant for (9) controlled with the stabilizing auxiliary control law $\bar{u}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)=\bar{u}_{\rm\scriptscriptstyle S}+\bar{K}_{\rm\scriptscriptstyle H}(\bar{x}^{{[N_{\rm\scriptscriptstyle L}]},o}(k)-\bar{x}_{\rm\scriptscriptstyle S})$.

The positive definite and symmetric 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 action, applied to system $\bar{\Sigma}_{w}^{[N_{\rm\scriptscriptstyle L}]}$ 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 $\Sigma$ 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 $\delta u_{i}\in{\mathcal{U}}_{{\rm\scriptscriptstyle S},i}\ominus{\bar{\mathcal{U}}}_{{\rm\scriptscriptstyle S},i}$ with the specific goal of 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}$: for $h=kN_{\rm\scriptscriptstyle L},\dots,(k+1)N_{\rm\scriptscriptstyle L}-1$

[TABLE]

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

Also denote by $\varDelta\Sigma_{i}$ the model given by the difference between systems (1) and (16), with (2) and (15).

[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 with (16) from the available value $\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$, not known in advance in the future prediction horizon. For this reason, we define a decentralized (approximated) dynamical system $\varDelta\hat{\Sigma}_{i}$ with input $\delta\hat{u}_{i}(h)$ and discarding all coupling inputs, i.e.,

[TABLE]

The decentralized dynamical system $\varDelta\hat{\Sigma}_{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 $\Sigma_{i}$, define a new sampling period $\zeta_{i}\in\mathcal{N}_{+}$ such that $N_{\rm\scriptscriptstyle L}/\zeta_{i}=N_{i}\in\mathcal{N}_{+}$ and a corresponding time index $l_{i}$. For clarity, the relations among the time scales with indices $h$, $l_{i}$, and $k$ are shown in Figure 2 in a specific case.

In the time scale of $l_{i}$, define the dynamical system $\Delta\hat{\Sigma}_{i}^{[\zeta_{i}]}$ as

[TABLE]

where ${B^{ii}_{\rm\scriptscriptstyle L}}^{[\zeta_{i}]}={\sum}_{j=0}^{\zeta_{i}-1}{(A^{ii}_{{\rm\scriptscriptstyle L}})}^{j}B^{ii}_{\rm\scriptscriptstyle L}$. In the short time scale, with time index $h$, the input $\delta\hat{u}_{i}(h)=\delta\hat{u}_{i}^{[\zeta_{i}]}(\lfloor h/\zeta_{i}\rfloor)$ is piecewise constant for $h\in[l_{i}\zeta_{i},(l_{i}+1)\zeta_{i}-1]$. Its value will be computed as the result of a suitable optimization problem formulated for system (19).

Given $\delta\hat{u}_{i}(h)$, the evolution of $\delta\hat{x}_{i}(h)$ can be computed thanks to the dynamical model (18), however $\delta\hat{x}_{i}(h)$ is in general different from $\delta x_{i}(h)$ due to the neglected interconnections in systems (18) and (19). For this reason, and 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]

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_{\rm\scriptscriptstyle 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 $h\in[kN_{\rm\scriptscriptstyle L},(k+1)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 (21) can only be formulated in an approximated way with reference to the state $\delta\hat{x}_{i}$ of system (18). In turn, if the resampling is used with $\zeta_{i}\neq 1$, the constraint on $\delta\hat{x}_{i}$ must be reformulated in terms of the state $\delta\hat{x}_{i}^{[\zeta_{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 $\bar{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 $\overrightarrow{\delta\hat{u}^{[\zeta_{i}]}_{i}}(kN_{i}:(k+1)N_{i}-1)=(\delta\hat{u}^{[\zeta_{i}]}_{i}(kN_{i}),\dots,\delta\hat{u}^{[\zeta_{i}]}_{i}((k+1)N_{i}-1))\in({\mathbb{R}}^{m_{i}})^{{N_{i}}}$, 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 $J_{{\rm\scriptscriptstyle L}}$ is a positive definite function with arguments $\delta\hat{x}^{[\zeta_{i}]}_{i}(kN_{i}+l_{i})$, $\delta\hat{u}^{[\zeta_{i}]}_{i}(kN_{i}+l_{i})$, $l_{i}=0,\dots,N_{i}-1$, e.g.,

[TABLE]

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

Finally, for $j=0,\dots,N_{\rm\scriptscriptstyle L}-1$, the control component at each (fast) time instant $\delta u_{i}(kN_{\rm\scriptscriptstyle L}+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 $N_{\rm\scriptscriptstyle L}$ fast sampling times. Due to resampling, the optimization problem computed at the beginning of each slow sampling time, i.e., at time $h=kN_{\rm\scriptscriptstyle L}$, has a prediction horizon of $N_{i}$ steps of lenght $\zeta_{i}$, where indeed $N_{i}\zeta_{i}=N_{\rm\scriptscriptstyle L}$. As a result of this optimization problem, the input sequence $\delta\hat{u}^{[\zeta_{i}]}_{i}(kN_{i}|kN_{\rm\scriptscriptstyle L})$, $\dots$, $\delta\hat{u}^{[\zeta_{i}]}_{i}((k+1)N_{i}-1|kN_{\rm\scriptscriptstyle L})$ is obtained, and at each fast sampling time $kN_{\rm\scriptscriptstyle L}+j$, $j=0,\dots,N_{\rm\scriptscriptstyle L}-1$, the real (low-level) input contribution (25) is used in (15). Note that, in this way, $\delta u_{i}(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 $\bar{u}^{[N_{\rm\scriptscriptstyle L}]}(k)$ using (14);

• •

for any system $\Sigma_{i}$, at the beginning of any long sampling time $k$ solve the optimization problem (23) with cost (24) and compute the full sequence $\delta\hat{u}^{[\zeta_{i}]}_{i}(kN_{i}|kN_{\rm\scriptscriptstyle L})$, $\dots$, $\delta\hat{u}^{[\zeta_{i}]}_{i}((k+1)N_{i}-1|kN_{\rm\scriptscriptstyle L})$;

• •

for any system $\Sigma_{i}$, at any fast sampling time $kN_{\rm\scriptscriptstyle L}+j$ $j=0,\dots,N_{\rm\scriptscriptstyle L}-1$, compute the low-level contribution $\delta u_{i}(kN_{\rm\scriptscriptstyle L}+j)$ using (25) and (26);

• •

for any system $\Sigma_{i}$, at any fast sampling time $kN_{\rm\scriptscriptstyle L}+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 $I_{si}=\begin{bmatrix}0_{\bar{n}_{i}\times\bar{n}_{1}}&\cdots&I_{\bar{n}_{i}}&\cdots&0_{\bar{n}_{i}\times\bar{n}_{M}}\end{bmatrix}$, ${\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}$, {\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}_{\rm\scriptscriptstyle S}}\subseteq\mathcal{B}_{\rho_{{u}}}(0)$. We now introduce the following assumption.

Assumption 3

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

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

[TABLE]

is full-rank with minimum singular value $\underline{\sigma}_{\mathcal{H}_{i}(N_{i})}>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.

Define $\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$, and where $A_{\rm\scriptscriptstyle L}^{\rm\scriptscriptstyle D}=$diag$(A_{\rm\scriptscriptstyle L}^{11},\dots,A_{\rm\scriptscriptstyle L}^{\rm\scriptscriptstyle MM})$ and $\rho_{\delta\hat{x}}(j)=\sqrt{\sum_{i=1}^{M}\rho^{2}_{\delta\hat{x}_{i}}(j)}$,

[TABLE]

We require that

[TABLE]

$\square$**

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 $A_{\rm\scriptscriptstyle L}$, and establishes a lower bound for the parameter $N_{\rm\scriptscriptstyle L}$. On the other hand, in view of Assumption 2.(2), $\beta_{i}$ are full rank and the full rank of matrix $\mathcal{R}_{i}(N_{i})$ 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 $N_{i}$ (and so $N_{\rm\scriptscriptstyle L}$) sufficiently large.

Secondly, Assumptions 3.(3)-3.(5) provide a tradeoff for the selection of parameters $\rho_{\delta\hat{u}_{i}}$, $i=1,\dots,M$ and $\rho_{\bar{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 $\delta{u}_{i}$, whose maximal amplitude is related to $\rho_{\delta\hat{u}_{i}}$) and centralized high-level controller (with input $\bar{u}$, whose maximal amplitude is related to $\rho_{\bar{u}}$);

• on the other hand, the amplitude of $\delta{u}_{i}$ 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 $N_{\rm\scriptscriptstyle L}$, as clarified in the following.

• 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.(3) can be fulfilled by taking $N_{\rm\scriptscriptstyle L}$ 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 $\bar{u}$ as much as possible - by taking a sufficiently large value of parameter $N_{\rm\scriptscriptstyle L}$. 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 $\mathcal{W}$ to be considered in the high level design is given by

[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}^{\rm\scriptscriptstyle D})\|\rho_{\delta\hat{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,\dots,M$

[TABLE]

*then

(i) $\bar{w}(k)\in\mathcal{W}$ and problems (11) and (23) 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]

(v) we can define (see (67) in Appendix 7.3) a function $\sigma(N_{\rm\scriptscriptstyle L})$ of $N_{\rm\scriptscriptstyle L}$ such that, if

[TABLE]

then, as $k\rightarrow\infty$, $x(kN_{\rm\scriptscriptstyle L})\rightarrow{x}_{\rm\scriptscriptstyle S}$. $\square$

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 $\sigma(N_{\rm\scriptscriptstyle L})$ is smaller than one, then convergence of the state to the goal is ensured. The definition of $\sigma(N_{\rm\scriptscriptstyle L})$ 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 $\sigma(N_{\rm\scriptscriptstyle L})$: the more the subsystems are interconnected (in a wide sense, regarding both the existence of dependencies between subsystems and their amplitude), the larger $\sigma(N_{\rm\scriptscriptstyle L})$. On the other hand, it is possible to reduce arbitrarily this parameter by increasing the tuning knob $N_{\rm\scriptscriptstyle L}$. This depends of the fact that, as $N_{\rm\scriptscriptstyle L}\rightarrow+\infty$, both $\|{\mathcal{A}}(N_{\rm\scriptscriptstyle L})\|\rightarrow 0$ and $\|{\mathcal{B}}(N_{\rm\scriptscriptstyle L})\|\rightarrow 0$.

A further remark is that, similarly to Proposition 2.3 in [29], for any $i=1,...,M$ it can be proved that $\lim_{N_{\rm\scriptscriptstyle L}\to+\infty}\lambda_{i}(N_{\rm\scriptscriptstyle L})=+\infty$, 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 $N_{\rm\scriptscriptstyle L}$, one can reduce at will the values of $\rho_{\delta\hat{u}_{i}}$, related to the maximum required amplitude of inputs $\delta u_{i}$. This, from (30) and (33), allows to reduce arbitrarily the high-level disturbance set $\mathcal{W}$. This, in turn, allows to reduce at will the corresponding RPI set $\mathcal{Z}$ and to minimize the conservativity of the present control scheme. We remark that, although a fine tuning of the gain matrices $K_{i}$, $i=1,\dots,M$ and $\bar{K}_{\rm\scriptscriptstyle H}$ can also be beneficial for the reduction of the conservativity of the scheme, the most relevant tuning knob indeed results parameter $N_{\rm\scriptscriptstyle L}$, 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 $N_{\rm\scriptscriptstyle L}\rightarrow+\infty$ 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 $N_{\rm\scriptscriptstyle L}$ 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 $N_{\rm\scriptscriptstyle L}\rightarrow+\infty$, 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 $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=\text{diag}(K_{1},\dots,K_{\rm\scriptscriptstyle 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 7.2) and of the sets ${\bar{\mathcal{U}}}_{{\rm\scriptscriptstyle S},i}$, $\Delta\hat{\mathcal{U}}_{i}$;

• •

computation of $\mathcal{W}$ according to (33) and (30);

• •

computation of $\bar{\mathcal{X}}_{\rm\scriptscriptstyle F}$, $\mathcal{Z}$, see [31], and $P_{\rm\scriptscriptstyle H}$ 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 $R_{1}$, $R_{2}$, $R_{3}$, three distillation columns $C_{1}$, $C_{2}$, $C_{3}$, 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 $D_{1}$, $D_{2}$, $D_{3}$ are the top products of the columns, while $B_{1}$, $B_{2}$, $B_{3}$ are their bottom products. The following reactions occur inside the reactors $R_{1}:\,A+B\rightarrow C+D$, $R_{2}:\,D+E\rightarrow F+B$, and $R_{3}:\,D+E\rightarrow F+B$. The system has six control variables, namely the refluxes ($v_{1},v_{3},v_{5}$) and the vapour ($v_{2},v_{4},v_{6}$) flow rates in the columns $C_{1},C_{2},C_{3}$, respectively, and six outputs: the liquid molar fraction ($r_{1}$) of component $A$ at the top product of $C_{1}$, the liquid molar fraction ($r_{2}$) of component $D$ at the bottom product of $C_{1}$, the liquid molar fraction ($r_{3}$) of component $C$ at the top product of $C_{2}$, the liquid molar fraction ($r_{4}$) of component $D$ at the bottom product of $C_{2}$, the liquid molar fraction ($r_{5}$) of component $B$ at the top product of $C_{3}$, and the liquid molar fraction ($r_{6}$) of component $F$ at the bottom product of $C_{3}$. 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 $v_{nominal}=[330\quad 410\quad 283\quad 385\quad 141\quad 282]^{{}^{\prime}}$ lb mol/h, to which it corresponds the vector of outputs $r_{nominal}=[0.942\quad 0.552\quad 0.827\quad 0.941\quad 0.705\quad 0.991]^{{}^{\prime}}$. 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 $R_{1},R_{2}$ and columns $C_{1},C_{2}$, while the second one is made by $R_{3}$ and $C_{3}$.

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 $\delta v_{i}$ and $\delta r_{i}$ the deviations of the inputs $v_{i}$ and the outputs $r_{i}$, respectively, with respect to their nominal values, the first resulting linear reduced order model is of order $n_{1}=25$, with $m_{1}=4$ inputs, i.e., $u_{1}=(\delta v_{1},\cdots,\delta v_{4})$, $p_{1}=4$ outputs, i.e., $y_{1}=(\delta r_{1},\cdots,\delta r_{4})$, and $p_{z1}=3$ coupling outputs, while the second one has $n_{2}=16$ states, $m_{2}=2$ inputs, i.e., $u_{2}=(\delta v_{5},\delta v_{6})$, $p_{2}=2$ outputs, i.e., $y_{2}=(\delta r_{5},\delta r_{6})$, and $p_{z2}=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 $u_{\rm\scriptscriptstyle S,i}$ is the steady state value corresponding to the output set-point values given by $y_{\rm\scriptscriptstyle S}=10^{-3}\cdot(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 $N_{\rm\scriptscriptstyle L}=28$, $N_{\rm\scriptscriptstyle H}=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 $\beta_{1}$ and $\beta_{2}$, and the reduced order model (4) with order $\bar{n}=6$. The dynamic matrix $A_{\rm\scriptscriptstyle H}={\rm 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 $B_{\rm\scriptscriptstyle H}$ has been computed as described in Appendix 7.1. The resulting model has been resampled with $N_{\rm\scriptscriptstyle L}=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 $\zeta_{1}=4$ and $\zeta_{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 $Q_{1}=10^{3}\cdot\beta_{1}^{\rm\scriptscriptstyle T}\beta_{1}$, $Q_{2}=10^{4}\cdot\beta_{2}^{\rm\scriptscriptstyle T}\beta_{2}$, $R_{1}=I_{m_{1}}$ and $R_{2}=I_{m_{2}}$.

• •

The decentralized feedback gains $K_{1}$ and $K_{2}$ guaranteeing that $F_{\rm\scriptscriptstyle L}$ 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 $F_{\rm\scriptscriptstyle L}$ is Schur stable.

• •

The parameters $\rho_{\delta\hat{u}_{1}}$, $\rho_{\delta\hat{u}_{2}}$, $\rho_{\bar{u}_{1}}$ and $\rho_{\bar{u}_{2}}$ have been computed according to (41) in Appendix 7.2. Specifically we used the cost function $J_{\rho}=\gamma_{1}\mathds{1}_{\rm\scriptscriptstyle 1\times M}\overrightarrow{\rho}_{\delta\hat{u}}-(\overrightarrow{\rho}_{\bar{u}}-\gamma_{2}\overrightarrow{\rho}_{u})^{{}^{\prime}}(\overrightarrow{\rho}_{\bar{u}}-\gamma_{2}\overrightarrow{\rho}_{u})$, where $\gamma_{1}$, $\gamma_{2}$ are positive constants selected as $\gamma_{1}=1$ and $\gamma_{2}=0.3$. Note that the choice of $\gamma_{1},\gamma_{2}$ allows one to modify the feasibility region of the high level optimization problem (11) and (12); typically setting $\gamma_{1}=1$, the feasibility properties of (11) grow with $\gamma_{2}$. The computed values are $\rho_{\delta\hat{u}_{1}}=145.2$, $\rho_{\delta\hat{u}_{2}}=115.7$, $\rho_{\bar{u}_{1}}=49.3$, $\rho_{\bar{u}_{2}}=25.2$. The input vectors of each subsystem, namely $\bar{u}_{i}$, $\delta\hat{u}_{i}$, $i=1,2$, have been limited as follows: $-\rho_{\bar{u}_{i}}/\sqrt{m_{i}}\mathds{1}_{m_{i}\times 1}\leq\bar{u}_{i}-\bar{u}_{\rm\scriptscriptstyle S,i}\leq\rho_{\bar{u}_{i}}/\sqrt{m_{i}}\mathds{1}_{m_{i}\times 1}$, $-\rho_{\delta\hat{u}_{i}}/\sqrt{m_{i}}\mathds{1}_{m_{i}\times 1}\leq\delta\hat{u}_{i}\leq\rho_{\delta\hat{u}_{i}}/\sqrt{m_{i}}\mathds{1}_{m_{i}\times 1}$; and the corresponding sets $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,i}$, $\Delta\hat{\mathcal{U}}_{\rm\scriptscriptstyle S,i}$, $i=1,2$ have been obtained.

Note that $\rho_{\bar{u}_{1}}/\sqrt{m_{1}}+\rho_{\delta\hat{u}_{1}}/\sqrt{m_{1}}=97.3$ and $\rho_{\bar{u}_{2}}/\sqrt{m_{2}}+\rho_{\delta\hat{u}_{2}}/\sqrt{m_{2}}=99.6$ and we can write $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,1}\oplus\Delta\hat{\mathcal{U}}_{\rm\scriptscriptstyle S,1}=\{\bar{u}_{1}+\delta\hat{u}_{1}|-97.3\cdot\mathds{1}_{m_{1}\times 1}\leq\bar{u}_{1}+\delta\hat{u}_{1}-\bar{u}_{\rm\scriptscriptstyle S,1}\leq 97.3\cdot\mathds{1}_{m_{1}\times 1}\}$ and $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,2}\oplus\Delta\hat{\mathcal{U}}_{\rm\scriptscriptstyle S,2}=\{\bar{u}_{2}+\delta\hat{u}_{2}|-99.6\cdot\mathds{1}_{m_{2}\times 1}\leq\bar{u}_{1}+\delta\hat{u}_{2}-\bar{u}_{\rm\scriptscriptstyle S,2}\leq 99.6\cdot\mathds{1}_{m_{2}\times 1}\}$. These results show that the size of both the sets $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,1}\oplus\Delta\hat{\mathcal{U}}_{\rm\scriptscriptstyle S,1}$ and $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,2}\oplus\Delta\hat{\mathcal{U}}_{\rm\scriptscriptstyle S,2}$ is close to that of the real input constraint in (36). In addition, the radius of the sets $\mathcal{B}_{\rho_{\Delta u_{i}(N_{\rm\scriptscriptstyle L}-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 $\rho_{\Delta u_{1}(N_{\rm\scriptscriptstyle L}-1)}=4.0$ and $\rho_{\Delta u_{2}(N_{\rm\scriptscriptstyle L}-1)}=0.37$. Therefore, we can also write $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,1}\oplus\Delta\mathcal{U}_{\rm\scriptscriptstyle S,1}=\{\bar{u}_{1}+\delta u_{1}|-99.3\cdot\mathds{1}_{m_{1}\times 1}\leq\bar{u}_{1}+\delta u_{1}-\bar{u}_{\rm\scriptscriptstyle S,1}\leq 99.3\cdot\mathds{1}_{m_{1}\times 1}\}$ and $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,2}\oplus\Delta\mathcal{U}_{\rm\scriptscriptstyle S,2}=\{\bar{u}_{2}+\delta u_{2}|-99.9\cdot\mathds{1}_{m_{2}\times 1}\leq\bar{u}_{1}+\delta u_{2}-\bar{u}_{\rm\scriptscriptstyle S,2}\leq 99.9\cdot\mathds{1}_{m_{2}\times 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 $Q_{\rm\scriptscriptstyle H}=I_{\bar{n}}$ and $R_{\rm\scriptscriptstyle H}=0.1I_{m}$.

• •

The gain matrix $\bar{K}_{\rm\scriptscriptstyle H}$ guaranteeing that both $F_{\rm\scriptscriptstyle H}$ and $F_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}$ 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 $F_{\rm\scriptscriptstyle H}$ and $F_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}$ were Schur stable.

• •

The disturbance set has been obtained according to (33) and (30) with $\rho_{w}=1.25$, and the RPI set has been computed with the algorithms described at pp. 231-233 of [30], i.e., $\mathcal{Z}=\{z|-(2.66,3.61,1.93,2.64,1.30,1.30)\leq z\leq\break(2.66,3.61,1.93,2.64,1.30,1.30)\}$ and $\bar{K}_{\rm\scriptscriptstyle H}\mathcal{Z}=\{\bar{K}_{\rm\scriptscriptstyle H}z|-(1.64,1.63,1.02,1.02,\break 0.10,0.10)\leq\bar{K}_{\rm\scriptscriptstyle H}z\leq(1.64,1.63,1.02,1.02,0.10,0.10)\}$. The terminal penalty $P_{\rm\scriptscriptstyle H}$ is computed with (13), and the terminal set is calculated under nominal input constraints in (11), i.e.,$\bar{\mathcal{X}}_{\rm\scriptscriptstyle F}=\{\bar{x}^{[N_{\rm\scriptscriptstyle L}],o}|(\bar{x}^{[N_{\rm\scriptscriptstyle L}],o}-\bar{x}_{\rm\scriptscriptstyle S})^{T}P_{\rm\scriptscriptstyle H}(\bar{x}^{[N_{\rm\scriptscriptstyle L}],o}-\bar{x}_{\rm\scriptscriptstyle S})\leq 0.997\}$, see [18].

The values taken by $\rho_{w}$, $\mathcal{Z}$ and $\bar{K}_{\rm\scriptscriptstyle H}\mathcal{Z}$ 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., ${\bar{u}_{\rm\scriptscriptstyle S}\oplus\bar{\mathcal{U}_{\rm\scriptscriptstyle S}}\ominus\bar{K}_{\rm\scriptscriptstyle H}\mathcal{Z}}$ 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 $\rho_{\bar{u}_{1}}$ and $\rho_{\bar{u}_{2}}$.

5.2.2 Tuning the hierarchical control scheme with $N_{\rm\scriptscriptstyle L}=84$, $N_{\rm\scriptscriptstyle H}=1$

Along the same line with the previous Section 5.2.1, the computational results corresponding to tuning knob $N_{\rm\scriptscriptstyle L}=84$ and $N_{\rm\scriptscriptstyle H}=1$ are also listed here.

• •

The following parameters $\rho_{\delta\hat{u}_{1}}$, $\rho_{\delta\hat{u}_{2}}$, $\rho_{\bar{u}_{1}}$ and $\rho_{\bar{u}_{2}}$ have been computed: $\rho_{\delta\hat{u}_{1}}=125.6$, $\rho_{\delta\hat{u}_{2}}=98.3$, $\rho_{\bar{u}_{1}}=59.5$, $\rho_{\bar{u}_{2}}=41.9$. The amplitude of both $\rho_{\delta\hat{u}_{1}}$ and $\rho_{\delta\hat{u}_{2}}$ ($\rho_{\bar{u}_{1}}$ and $\rho_{\bar{u}_{2}}$) is smaller (larger) than that with tuning knob $N_{\rm\scriptscriptstyle L}=28$, $N_{\rm\scriptscriptstyle H}=3$.

• •

The sets $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,1}\oplus\Delta\hat{\mathcal{U}}_{\rm\scriptscriptstyle S,1}=\{\bar{u}_{1}+\delta\hat{u}_{1}|-92.6\cdot\mathds{1}_{m_{1}\times 1}\leq\bar{u}_{1}+\delta\hat{u}_{1}-\bar{u}_{\rm\scriptscriptstyle S,1}\leq 92.6\cdot\mathds{1}_{m_{1}\times 1}\}$ and $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,2}\oplus\Delta\hat{\mathcal{U}}_{\rm\scriptscriptstyle S,2}=\{\bar{u}_{2}+\delta\hat{u}_{2}|-99.1\cdot\mathds{1}_{m_{2}\times 1}\leq\bar{u}_{1}+\delta\hat{u}_{2}-\bar{u}_{\rm\scriptscriptstyle S,2}\leq 99.1\cdot\mathds{1}_{m_{2}\times 1}\}$. These results show that the size of both the sets $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,1}\oplus\Delta\hat{\mathcal{U}}_{\rm\scriptscriptstyle S,1}$ and $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,2}\oplus\Delta\hat{\mathcal{U}}_{\rm\scriptscriptstyle S,2}$ has been reduced slightly compared to that with tuning knob $N_{\rm\scriptscriptstyle L}=28$, $N_{\rm\scriptscriptstyle H}=3$, however, is still close to that of the real input constraint in (36). In addition, the radius of the sets $\mathcal{B}_{\rho_{\Delta u_{i}(N_{\rm\scriptscriptstyle L}-1)}}(0)$, $i=1,2$, have been computed with $\rho_{\Delta u_{1}(N_{\rm\scriptscriptstyle L}-1)}=11.3$ and $\rho_{\Delta u_{2}(N_{\rm\scriptscriptstyle L}-1)}=0.91$. Therefore, we can also write $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,1}\oplus\Delta\mathcal{U}_{\rm\scriptscriptstyle S,1}=\{\bar{u}_{1}+\delta u_{1}|-98.2\cdot\mathds{1}_{m_{1}\times 1}\leq\bar{u}_{1}+\delta u_{1}-\bar{u}_{\rm\scriptscriptstyle S,1}\leq 98.2\cdot\mathds{1}_{m_{1}\times 1}\}$ and $\bar{\mathcal{U}}_{\rm\scriptscriptstyle S,2}\oplus\Delta\mathcal{U}_{\rm\scriptscriptstyle S,2}=\{\bar{u}_{2}+\delta u_{2}|-99.8\cdot\mathds{1}_{m_{2}\times 1}\leq\bar{u}_{1}+\delta u_{2}-\bar{u}_{\rm\scriptscriptstyle S,2}\leq 99.8\cdot\mathds{1}_{m_{2}\times 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 $\rho_{w}=1.56$, and the RPI set has been computed i.e., $\mathcal{Z}=\{z|-(1.83,2.11,1.69,1.82,1.58,1.58)\leq z\leq\break(1.83,2.11,1.69,1.82,1.58,1.58)\}$ and $\bar{K}_{\rm\scriptscriptstyle H}\mathcal{Z}=\{\bar{K}_{\rm\scriptscriptstyle H}z|-(0.37,0.29,0.17,0.18,\break 0.01,0.01)\leq\bar{K}_{\rm\scriptscriptstyle H}z\leq(0.37,0.29,0.17,0.18,0.01,0.01)\}$. The terminal penalty $P_{\rm\scriptscriptstyle H}$ is computed with (13), and the terminal set is calculated under nominal input constraints in (11), i.e.,$\bar{\mathcal{X}}_{\rm\scriptscriptstyle F}=\{\bar{x}^{[N_{\rm\scriptscriptstyle L}],o}|(\bar{x}^{[N_{\rm\scriptscriptstyle L}],o}-\bar{x}_{\rm\scriptscriptstyle S})^{T}P_{\rm\scriptscriptstyle H}(\bar{x}^{[N_{\rm\scriptscriptstyle L}],o}-\bar{x}_{\rm\scriptscriptstyle S})\leq 0.98\}$.

Note that, the size of both the sets $\mathcal{Z}$ and $\bar{K}_{\rm\scriptscriptstyle H}\mathcal{Z}$ is smaller than that with the tuning $N_{\rm\scriptscriptstyle L}=28$, $N_{\rm\scriptscriptstyle H}=3$. The amplitude of $\rho_{\bar{u}_{1}}$ and $\rho_{\bar{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 $y_{\rm\scriptscriptstyle S}=10^{-3}\cdot(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}\cdot(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=\text{diag}(Q_{1},\,Q_{2})$, input penalty $R=I_{m}$ and prediction horizon $N=N_{\rm\scriptscriptstyle H}\cdot N_{\rm\scriptscriptstyle L}=84$. The terminal set has been chosen as $\mathcal{X}_{\rm\scriptscriptstyle F}=\{x|(x-x_{\rm\scriptscriptstyle S})^{T}P(x-x_{\rm\scriptscriptstyle S})\leq 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 $[kN_{\rm\scriptscriptstyle L},\,kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L})$ with $N_{\rm\scriptscriptstyle L}=28$ and $N_{\rm\scriptscriptstyle H}=3$ ($N_{\rm\scriptscriptstyle L}=84$ and $N_{\rm\scriptscriptstyle H}=1$), i.e., $2.49$ s ($6.87$ s), is reduced dramatically compared to that of the centralized MPC, i.e., $21.9\cdot N_{\rm\scriptscriptstyle L}$ 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 $N_{\rm\scriptscriptstyle L}=28$, $N_{\rm\scriptscriptstyle H}=3$ and $N_{\rm\scriptscriptstyle L}=84$, $N_{\rm\scriptscriptstyle H}=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 $N_{\rm\scriptscriptstyle L}=84$ and $N_{\rm\scriptscriptstyle H}=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 $\hat{Q}_{1}=0.1I_{n_{1}}$, $\hat{Q}_{2}=0.1I_{n_{2}}$, while the covariances of the output noises have been chosen as $\hat{R}_{1}=0.01I_{p_{1}}$, $\hat{R}_{2}=0.01I_{p_{2}}$. Finally, the covariances of the initial state estimates have been selected as $P_{1}(0)=0.01I_{n_{1}}$, $P_{2}(0)=0.01I_{n_{2}}$.

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 $r_{nominal}$; 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 $\beta_{i}$ and of the reduced order model

A constructive procedure for the computation of the matrices $\beta_{i}$, $i=1,\dots,\,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 $\beta$. 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 $\bar{n}\geq p+$dim$($Im$G_{\rm\scriptscriptstyle L}^{x}\cap$Ker$C_{\rm\scriptscriptstyle L})$) is, in our case, only necessary.

Here we now describe a (possibly conservative) procedure:

a

find a subspace ${\rm Ker}\,\beta_{i}$ of dimension $n_{i}-\bar{n}_{i}$ so that ${\rm Ker}\,\beta_{i}\subseteq{\rm Ker}\,C_{{\rm\scriptscriptstyle L}}^{ii}$ and ${\mathcal{Z}_{ii}}\cap{\rm Ker}\,\beta_{i}=\{0\}$, where $\mathcal{Z}_{ii}={\rm Im}\,\tilde{G}_{\rm\scriptscriptstyle L}^{ii}(1)\cap{\rm Ker}\,C_{\rm\scriptscriptstyle L}^{ii}$, and $\tilde{G}_{\rm\scriptscriptstyle L}^{ii}(z)=(zI-A_{\rm\scriptscriptstyle L}^{ii})^{-1}\begin{bmatrix}B_{\rm\scriptscriptstyle L}^{ii}&E_{i}\end{bmatrix}$;

b

let $\{\kappa_{1},\dots,\kappa_{n_{i}-\bar{n}_{i}}\}$ be a set of independent vectors in $\mathcal{K}_{\beta,i}=\rm Ker\,\beta_{i}$ and complete this set to a basis ${\mathcal{B}}_{i}=\{v_{1},\dots,v_{\bar{n}_{i}},\kappa_{1},\dots,\kappa_{n_{i}-\bar{n}_{i}}\}$ of the whole space ${\mathbb{R}}^{n_{i}}$;

c

let $\{r_{1},\dots,r_{\bar{n}_{i}}\}$ be a basis of ${\mathbb{R}}^{\bar{n}_{i}}$, define

[TABLE]

and $T_{{\rm{\scriptscriptstyle L}}}$ be the matrix whose columns are the vectors in ${\mathcal{B}}_{i}$, then $\beta_{i}=\hat{\beta}_{i}T_{{\rm{L}}}^{-1}$;

d

define collective matrix $\beta=\rm{diag}(\beta_{1},\dots,\,\beta_{\rm\scriptscriptstyle M})$;

e

choose matrix $A_{\rm\scriptscriptstyle H}$ being Schur stable, and let

[TABLE]

where the suitable choices for $A_{\rm\scriptscriptstyle H}$ 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, $\tilde{G}_{\rm\scriptscriptstyle L}^{ii}(z)=(zI-A_{\rm\scriptscriptstyle L}^{ii})^{-1}B_{\rm\scriptscriptstyle L}^{ii}$. However this choice does not guarantee a-priori that the property $\mathcal{Z}\cap\mathcal{K}_{\beta}=\{0\}$, where $\mathcal{Z}=$Im$G_{\rm\scriptscriptstyle L}^{x}\cap$Ker$C_{\rm\scriptscriptstyle L}$, $G_{\rm\scriptscriptstyle L}^{x}=(I-A_{\rm\scriptscriptstyle L})^{-1}B_{\rm\scriptscriptstyle L}$, and $\mathcal{K}_{\beta}=\prod_{i=1}^{M}\mathcal{K}_{\beta}^{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 $\bar{\mathcal{U}}_{{\rm\scriptscriptstyle S},i}$ and $\Delta\hat{\mathcal{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,\dots,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 $\Delta\hat{\mathcal{U}}_{i}=\mathcal{B}_{\rho_{\delta\hat{u}_{i}}}(0)$ and $\bar{\mathcal{U}}_{{\rm\scriptscriptstyle S},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}_{\rm\scriptscriptstyle M}})$ and $\overrightarrow{\rho}_{\bar{u}}=(\rho_{\bar{u}_{1}},\dots,\rho_{\bar{u}_{\rm\scriptscriptstyle 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,\dots,M$. Here the constant $\rho_{\bar{u}}$ appears, defined in such a way that $\bar{\mathcal{U}}_{{\rm\scriptscriptstyle S}}=\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 (28) it is sufficient to verify the following matrix inequality

[TABLE]

where $\mathds{1}_{\rm\scriptscriptstyle M\times M}$ is the $M\times 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,\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}^{\rm\scriptscriptstyle 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}\subseteq\mathcal{B}_{\rho_{\delta\hat{u}_{i}}+\sum_{j=1}^{M}\lambda_{ij}\rho_{\delta\hat{u}_{j}}}(0)$. Therefore, to verify (38) it is sufficient to enforce the constraint

[TABLE]

where ${\Lambda}$ is the $M\times M$ matrix whose entries are $\lambda_{ij}$, $i,j=1,\dots,M$, while $\overrightarrow{\rho}_{{u}}=(\rho_{u_{1}},\dots,\rho_{u_{\rm\scriptscriptstyle M}})$, where $\mathcal{B}_{\rho_{i}}(0)\subseteq\mathcal{U}_{{\rm\scriptscriptstyle S},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 problem:

[TABLE]

where $J_{\rho}$ is a suitable (linear or quadratic, if possible) cost function that allows to maximize the size of the constraint set.

7.3 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\bar{u}_{\rm\scriptscriptstyle S}\oplus{\bar{\mathcal{U}}_{\rm\scriptscriptstyle S}}$ there exists a feasible sequence$\overrightarrow{\delta\hat{u}^{[\zeta_{i}]}_{i}}{(kN_{i}:(k+1)N_{i}-1|kN_{\rm\scriptscriptstyle L})}\in{\Delta\hat{\mathcal{U}}_{i}}^{N_{i}}$ such that the terminal constraint (22) is satisfied.

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

Proof of Proposition 1

A) Note that, since $\delta\hat{x}_{i}^{[\zeta_{i}]}(kN_{i})=\delta\hat{x}_{i}(kN_{\rm\scriptscriptstyle L})=0$,

[TABLE]

Moreover, in view of (5)

[TABLE]

Analogously, from (16) written in collective form

[TABLE]

Therefore, in view of (43), (44a), (45a), and the definition ${\mathcal{A}}(N_{\rm\scriptscriptstyle L})$, ${\mathcal{B}}(N_{\rm\scriptscriptstyle L})$, $I_{si}$, the constraint (22) can be written as

[TABLE]

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

[TABLE]

from which the result follows.

B) From (3) it holds that

[TABLE]

where $\otimes$ is the Kronecker product. Therefore

[TABLE]

and, in view of (42)

[TABLE]

for all $i=1,\dots,M$. From this expression and Assumption 3.(4)

[TABLE]

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

Proposition 2

*If Problem (23) is feasible at time $h=kN_{\rm\scriptscriptstyle L}$, then

[TABLE]

Also it holds that

[TABLE]

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

Proof of Proposition 2

Defining the collective vectors $\hat{x}=(\hat{x}_{1},\dots,\hat{x}_{\rm\scriptscriptstyle M})$, $\delta x=(\delta x_{1},\dots,\delta x_{\rm\scriptscriptstyle M})$,$\delta\hat{x}=(\delta\hat{x}_{1},\dots,\delta\hat{x}_{\rm\scriptscriptstyle 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 (23) is feasible, and therefore equality (22) is verified. From (17), (18), (20), 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)=\break\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}^{\rm\scriptscriptstyle 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 (30). 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\break\sqrt{\sum_{i=1}^{M}\rho^{2}_{\delta\hat{x}_{i}}(j)}$. Therefore one has (52) for all $k\geq 0$.

From (55) 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}^{\rm\scriptscriptstyle 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})=\break K_{i}I_{si}\sum_{r=2}^{j}F_{\rm\scriptscriptstyle L}^{j-r}(A_{\rm\scriptscriptstyle L}-A_{\rm\scriptscriptstyle L}^{\rm\scriptscriptstyle 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 (54). $\square$

Proof of Theorem 1

(i) If $\|x(0)-x_{\rm\scriptscriptstyle S}\|\leq\lambda_{i}(N_{\rm\scriptscriptstyle L})$ 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\geq 0$, a feasible sequence $\overrightarrow{\delta\hat{u}^{[\zeta_{i}]}_{i}}{(kN_{i}:(k+1)N_{i}-1|kN_{\rm\scriptscriptstyle L})}\in{\Delta\hat{\mathcal{U}}_{i}}^{N_{i}}$ such that the terminal constraint (22) is satisfied.

Also, from Proposition 2, it follows that $\bar{w}(k)\in\mathcal{W}$ for all $k\geq 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), $\bar{u}^{[N_{\rm\scriptscriptstyle L}]}(k)\in\bar{u}_{\rm\scriptscriptstyle S}\oplus\bar{\mathcal{U}}_{\rm\scriptscriptstyle S}$; also, from Proposition 2 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 (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 $\bar{x}^{[N_{\rm\scriptscriptstyle L}],o}(k)\rightarrow\bar{x}_{\rm\scriptscriptstyle S}$ as $k\rightarrow+\infty$, and that $\bar{x}^{[N_{\rm\scriptscriptstyle L}]}(k)$ is asymptotically driven to lie in the robust positively invariant set $\bar{x}_{\rm\scriptscriptstyle S}\oplus\mathcal{Z}$.

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

[TABLE]

Denoting $\mathcal{B}_{{\rm\scriptscriptstyle L},N_{\rm\scriptscriptstyle L}}^{\rm\scriptscriptstyle C}=\begin{bmatrix}A_{\rm\scriptscriptstyle L}^{N_{\rm\scriptscriptstyle L}-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)+\break{\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]

then

[TABLE]

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

[TABLE]

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

[TABLE]

we can rewrite (57) as

[TABLE]

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

[TABLE]

Based on this, we define $\kappa_{\delta u}=\|\mathcal{B}_{{\rm\scriptscriptstyle L},N_{\rm\scriptscriptstyle L}}^{\rm\scriptscriptstyle C}(I+{\rm{diag}}(K_{i},\dots,K_{i})\\ \mathcal{F}_{N_{\rm\scriptscriptstyle L}}{\rm{diag}}(A_{\rm\scriptscriptstyle L}^{\rm\scriptscriptstyle C},\dots,A_{\rm\scriptscriptstyle L}^{\rm\scriptscriptstyle C})\mathcal{B}_{\rm\scriptscriptstyle L})\|$ and we write (59) 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}}}$. Since $x_{\rm\scriptscriptstyle S}=A_{\rm\scriptscriptstyle L}^{N_{\rm\scriptscriptstyle L}}x_{\rm\scriptscriptstyle S}+B_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}\bar{u}_{\rm\scriptscriptstyle S}$ and $B_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}\bar{K}_{\rm\scriptscriptstyle H}\bar{x}_{\rm\scriptscriptstyle S}=B_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}\bar{K}_{\rm\scriptscriptstyle H}\beta{x}_{\rm\scriptscriptstyle S}$, we can rewrite (61) as

[TABLE]

Eventually, 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}}}$.

(v) We can reformulate problem (23) as

[TABLE]

where $\mathcal{Q}=$diag$(R,\dots,R)+\mathcal{B}_{\rm\scriptscriptstyle L}^{T}$diag$(Q,\dots,Q)\mathcal{B}_{\rm\scriptscriptstyle L}$, $R=$diag

$(R_{1},\dots,R_{\rm\scriptscriptstyle M})$, $Q=$diag$(Q_{1},\dots,Q_{\rm\scriptscriptstyle M})$, $\delta\bar{x}_{end}=\bar{x}^{[N_{\rm\scriptscriptstyle L}]}(k+1|k)-\beta\hat{x}(kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L})$ and$\mathcal{R}^{\rm\scriptscriptstyle D}(N_{\rm\scriptscriptstyle L})=\begin{bmatrix}(A_{\rm\scriptscriptstyle L}^{\rm\scriptscriptstyle D})^{N_{\rm\scriptscriptstyle L}-1}B_{\rm\scriptscriptstyle L}&\dots&B_{\rm\scriptscriptstyle L}\end{bmatrix}$. Recalling that $b_{in}>0$ elementwise, in view of continuity arguments, there exists a ball $\mathcal{B}_{\rho_{end}}(0)$ for $\delta\bar{x}_{end}$ such that the solution to problem (62) satisfies $A_{in}\overrightarrow{\delta\hat{u}}{(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1)}\leq b_{in}$. If $\delta\bar{x}_{end}\in\mathcal{B}_{\rho_{end}}(0)$, then

[TABLE]

Also, it holds that the optimal constrained solution fulfills

[TABLE]

and therefore $\|\overrightarrow{\delta\hat{u}}{(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1)}\|\leq\sqrt{\frac{J_{\rm\scriptscriptstyle L}^{max}}{\lambda_{min}(\mathcal{Q})}}$; in view of this, for $\|\delta\bar{x}_{end}\|\not\in\mathcal{B}_{\rho_{end}}(0)$, then $\|\overrightarrow{\delta\hat{u}}(kN_{\rm\scriptscriptstyle L}:kN_{\rm\scriptscriptstyle L}+N_{\rm\scriptscriptstyle L}-1)\|\leq\sqrt{\frac{J_{\rm\scriptscriptstyle L}^{max}}{\lambda_{min}(\mathcal{Q})}}\frac{1}{\rho_{end}}\|\delta\bar{x}_{end}\|$. Defining now $\bar{\kappa}=\max\{\left\|\begin{bmatrix}I&0\end{bmatrix}\begin{bmatrix}2\mathcal{Q}&(\beta\mathcal{R}^{\rm\scriptscriptstyle D}(N_{\rm\scriptscriptstyle L}))^{T}\\ \beta\mathcal{R}^{\rm\scriptscriptstyle D}(N_{\rm\scriptscriptstyle L})&0\end{bmatrix}^{-1}\begin{bmatrix}0\\ I\end{bmatrix}\right\|,\sqrt{\frac{J_{\rm\scriptscriptstyle L}^{max}}{\lambda_{min}(\mathcal{Q})}}\frac{1}{\rho_{end}}\}$ then we conclude that

[TABLE]

Therefore, from (60) and (63), we have that

[TABLE]

where

[TABLE]

In view of this, we can rewrite (59) as

[TABLE]

where $\|w_{x}(k)\|\leq\kappa_{x}(N_{\rm\scriptscriptstyle L})\|x(kN_{\rm\scriptscriptstyle L})-x_{\rm\scriptscriptstyle S}\|$, $\kappa_{x}(N_{\rm\scriptscriptstyle L})=\kappa_{u}\bar{\kappa}\|\mathcal{B}_{{\rm\scriptscriptstyle L},N_{\rm\scriptscriptstyle L}}^{\rm\scriptscriptstyle C}\|\\ \|{\mathcal{A}}(N_{\rm\scriptscriptstyle L})+{\mathcal{B}}(N_{\rm\scriptscriptstyle L})\bar{K}_{\rm\scriptscriptstyle H}\beta\|$, $\|w_{o}(k)\|\leq\kappa_{u}^{o}\|\bar{u}^{[N_{\rm\scriptscriptstyle L}],o}(k)-\bar{u}_{\rm\scriptscriptstyle S}\|+\kappa_{x}^{o}\|\bar{x}^{[N_{\rm\scriptscriptstyle L}],o}(k)-\bar{x}_{\rm\scriptscriptstyle S}\|$, $\kappa_{u}^{o}=\|B_{\rm\scriptscriptstyle L}^{[N_{\rm\scriptscriptstyle L}]}\|+\kappa_{u}\bar{\kappa}\|\mathcal{B}_{{\rm\scriptscriptstyle L},N_{\rm\scriptscriptstyle L}}^{\rm\scriptscriptstyle C}\|\|{\mathcal{B}}(N_{\rm\scriptscriptstyle L})\|$, and

[TABLE]

To derive a stability condition, we recast (65) as the following redundant dynamic system

[TABLE]

where the initial conditions for $\delta x_{1}$ and $\delta x_{2}$ coincide and are equal to $x(0)-x_{\rm\scriptscriptstyle S}$, and where $\|w_{\delta x_{i}}\|\leq\kappa_{x}(N_{\rm\scriptscriptstyle L})\|\delta x_{i}\|$, $i=1,2$. Recall also that, as already discussed, $w_{o}$ 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, $\Gamma$ is Schur stable if and only if the inequality (35) is verified. Then, under the latter condition, $x(kN_{\rm\scriptscriptstyle L})\rightarrow x_{\rm\scriptscriptstyle S}$ as $k\rightarrow+\infty$. $\square$