I. INTRODUCTION

Many real systems, such as solar thermal central receivers, economic systems, and manufacturing systems, are subject to random abrupt changes in their inputs, internal variables and other system parameters^[1]. It is difficult to describe the above systems using conventional linear time-invariant systems. However, Markovian jump linear systems (MJLSs)^{[2, 3, 4]} have been put forward to successfully solve this problem. Therefore, in the past decades MJLSs have received considerable interest, and we can find examples in ^{[5, 6, 7]}. Since it is often difficult to identify exact model parameters, the control of MJLSs with uncertainty in model parameters has been taken an important consideration of ^[1]. With the increasing development of computer science and technology, discrete-time MJLSs has been applied in a receding horizon control scheme^[8].

Model predictive control (MPC), which is also called receding horizon control, has drawn extensive attention from both theoretical and industrial communities with its advantages in processing large-scale, constraint, multi-variable processes. Nowadays, the systems to be handled are increasingly complex^[9]. It is pretty difficult to control these systems with a centralized MPC control strategy due to the computational complexity and communication bandwidth limitations^[10]. Recently, the distributed MPC framework where each subsystem is controlled by an independent controller has been increasingly focused on for its flexibility of system structure, error-tolerance, less computational efforts and no requirements for global information^[11]. In distributed MPC algorithm, we divide the global system into several subsystems and each subsystem exchanges information with others via network. In ^[12], Jia et al. provided a communication based on MPC scheme where each subsystem optimized a local cost function and made coordination via communication. Mercangöz and Doyle^[13] put forward an iterative implementation of a resembling distributed MPC scheme and applied it to a four-tank system. These communication based on distributed MPC algorithm only achieved a Nash-optimization^[14]. A robust distributed MPC algorithm was proposed in ^{[15, 16]}, in which the model uncertain problem was converted into solving a linear matrix inequality (LMI) optimization problem. In ^[17], Ding addressed the distributed MPC of a set of polytopic uncertain local systems with decoupled dynamics.

In this paper, we propose a distributed MPC algorithm for a class of uncertain discrete-time MJLSs subject to actuator saturation. We use state feedback to ensure the stability of system. Thus, we assume that system states are measurable. We decompose the global system into several subsystems and a min-max distributed MPC strategy is proposed for the polytopic uncertain distributed systems. In order to better reflect the practical system under the networked environment, we consider the MJLSs subject to actuator saturation. At each time step an MPC algorithm optimizes a cost function, which is associated with the future states and manipulated variables^[9], to obtain a control input. In next time step, a new cost function is formulated and solved based on the new measurements. An iterative algorithm is developed to get sub-optimizations of the neighbor subsystems.

The rest of this paper is organized as follows. Section I describes the problem statement. In Section II, based on the state feedback control law, a distributed MPC algorithm for input-saturated polytopic uncertain systems is proposed. The effectiveness of the proposed approach is illustrated in Section III with the help of numerical examples. Finally, we give a conclusion to our work in Section IV.

${\bf Notations.}$ The following notations are used throughout the paper. ${{\bf R}}^n$ and ${{\bf R}}^{n\times m}$ denote the $n$-dimensional Euclidean space and the set of all $n\times m$ real matrices, respectively. ${x_i}(k + n|k)$ denotes the state of subsystem $i$ at time $k+n$, predicted based on the measurements at time $k. u_i(k+n|k)$ refers to the predicted control move at time step $k+n$ and $u_i(k|k)$ is the control move to be implemented at time step $k$. $\|x\|^2$ denotes the Euclidean norm of the vector $x$. $I$ and 0 denote the identity matrix and zero matrix of compatible dimension, respectively. ${\rm E}\{x\}$ stands for the expectation of stochastic variable $x$. $G>0$ means matrix $G$ is real symmetric and positive definite. The symbol "$*$" denotes the symmetric part in a symmetric matrix.

II. PROBLEM STATEMENT

Consider a certain discrete-time MJLS subject to actuator saturation described by

\begin{align}\label{1} \begin{cases} x(k+1)=A^{(\xi)}(r_k)x(k)+B^{(\xi)}(r_k)\sigma(u(k)), \\[1mm] y(k)=C^{(\xi)}(r_k)x(k), \end{cases} \end{align}

(1)

where $x(k)\in {\bf R}^ {n_x}$ is the system state; $u(k)\in {\bf R}^{n_u}$ is the control input; $y(k)\in {\bf R}^{n_y}$ is the controlled output; $r_k$ is the system mode. The initial state is $x_0$, and the initial mode is $r_0$. According to the definition in ^[6], we can see the true "state" involves in two parts: the continuous parts, i.e., $x(k)$, and the discrete part, i.e., the mode $r_k$. For a fixed mode $r_k$, the MJLS is a linear time-varying system in essence. In another word, mode $r_k$, i.e., the present mode, is known, while the next mode $r_{k+1}$ is stochastic and uncertain. Let $r_k$ $(k\in [0, N])$ be a Markov chain taking values in a finite state space ${M}=\{1$, $2, \ldots, S\}$ with transition probability given by

\begin{align}\label{2} p(g, h)={\rm prob}\left\{r_{k+1}=h|r_k=g\right\}, \ \ \forall g, h\in {M}, \end{align}

(2)

where $p(g, h)\ge0$ $(g, h\in{{M}})$ is the transition probability from $g$ to $h$ and $\sum_{h=1}^S p(g, h)=1$, $\forall g\in{M}$.

It is assumed that for $\forall r_k \in {M}$, $A^{(\xi)}(r_k), B^{(\xi)}(r_k)$ and $C^{(\xi)}(r_k)$ are unknown matrices which are involved in a convex polyhedral set $\Omega(r_k)$ described by $L$ vertices

\begin{align*} \Omega(r_k)\in \text{Co}\Big\{&\left[A^{(1)}(r_k)~~B^{(1)}(r_k)~~C^{(1)}(r_k)\right]\notag\\[1mm] &\left[A^{(2)}(r_k)~~B^{(2)}(r_k)~~C^{(2)}(r_k)\right]\notag\\[1mm] &\qquad\qquad\qquad\vdots\notag\\[1mm] &\left[A^{(L)}(r_k)~~B^{(L)}(r_k)~~C^{(L)}(r_k)\right]\Big\}, \notag \end{align*}

where Co refers to the convex hull, that is,

\begin{align} &\left[A^{(\xi)}(r_k)~~B^{(\xi)}(r_k)~~C^{(\xi)}(r_k)\right]=\notag\\[1mm] &\qquad \sum_{l=1}^L\lambda_l\left[A^{(l)}(r_k)~~ B^{(l)}(r_k)~~ C^{(l)}(r_k)\right], \end{align}

(3)

where $\sum_{l=1}^L\lambda_l=1$, $\lambda_l\geq0$.

${\bf Remark 1.}$ The function $\sigma:{\bf R}^{m}\rightarrow {\bf R}^{m}$ is the standard saturation function of appropriate dimensions defined as $\sigma(u(k))$ $=$ $[\sigma(u_1(k)) ~\sigma(u_2(k))~\ldots~\sigma(u_m(k))]$ where $u_i(k)$ $=$ ${\rm sgn}{(u)}{\min}{\{1, |u|\}}$, where the notation of "sgn" denotes the signum function. To be convenient, we define $\bar{u}_i(k)$ $=$ $\sigma{(u_i(k))}$.

We decompose the global system into $M$ subsystems and design the corresponding distributed MPC controllers for them. In the following, we give some assumptions for our discussion.

${\bf Assumption 1.}$ The vector of states $x^{\rm T}_i=[x^{\rm T}_{11}, \ldots, x^{\rm T}_{ii}, $ $\ldots, $ $x^{\rm T}_{MM}]^{\rm T}$, $x_{ii}\in {\bf R}^{n_{x_{ii}}}$ includes all the states of the system, that is, the local state $x_{ii}$ that can be measured or estimated, while the other subsystems' states could be obtained by communication via network. Thus, the states of subsystems are the same to the global system, that is to say $x_1=x_2=\cdots$ $=$ $x$, we only discuss the inputs decomposition in this paper. For simplicity, we denote the state value of $i$-th subsystem as $x_i$ to stress the vector to be used for computing the local control input $u_i$.

${\bf Assumption 2.}$ In this paper, we neglect the effect of packet loss, packet disorder, time delay, etc.

Thus, the distributed polytopic uncertain system models with actuator saturation by decomposition can be described by

\begin{align}\label{4} \begin{cases} x_i(k+1)=A_i^{(\xi)}(r_k)x_i(k)+B_i^{(\xi)}(r_k)\sigma(u_i(k))\\ ~~~~~~~~~~~~~~ +\sum\limits_{j=1, j\neq i}^M B_j^{(\xi)}(r_k)\sigma(u_j(k)), \\[5mm] y_i(k)=C_i^{(\xi)}(r_k)x_i(k), ~~~~i=1, 2, \ldots, M, \\ \end{cases} \end{align}

(4)

where $x_i\in {\bf R}^{n_x}$ is the system state of subsystem $i$; $u_i\in {\bf R}^{n_{u_i}}$ is the control input of subsystem $i$; and $y_i\in {\bf R}^{n_{y_i}}$ is the controlled output of subsystem $i$.

\begin{align}\label{5} &[A_i^{(\xi)}(r_k)~~B_i^{(\xi)}(r_k)~~C_i^{(\xi)}(r_k)]\notag\\[1mm] &\qquad =\sum_{l=1}^L\lambda_l\left[A^{(l)}_i(r_k)~~ B^{(l)}_i(r_k)~~ C^{(l)}_i(r_k)\right]. \end{align}

(5)

${\bf Remark 2.}$ It is important to reduce the computation complexity to apply the MPC algorithm to the practical system. In ^[18], by employing a scalar $a$, an MPC algorithm subject to certain LMIs was presented. In this paper, by decomposing the control input, that is, $n_u=\sum_{i=1}^M n_{u_i}$, system (4) becomes a distributed MPC problem. Compared with the centralized MPC, the distributed MPC strategies can reduce the computational complexity and relax communication bandwidth restrictions. Besides, with the inevitable parameter uncertainties taken into consideration, the problem is also a robust MPC.

${\bf Remark 3.}$ Owing to the limitation of the communication bandwidth, the large information-carry packets may be saturated at actuator node. The actuator saturation is an essential issue in the control community, which leads to nonlinearity of the system and has being focused on in recent years such as in ^[20]. In this paper, we discuss the distributed MPC problem for MJLs with actuator saturation.

The aim of this paper is to find a set of state feedback laws as follows:

\begin{align}\label{6} u_i(k+n|k)=F_i(k, r_{k+n|k})x_i(k+n|k), \end{align}

(6)

in order to stabilize the system.

For a given $F_i\in {\bf R}^{m_i\times {n_x}}$, let $\phi(F_i)=\{x_i\in {\bf R}^{n_x}:|f^{\rm T}_{i, \alpha}|$ $\le$ $1, $ $\alpha=1, \ldots, m_i\}$, where $f^{\rm T}_{i, \alpha}$ is the $\alpha$-th row of matrix $F_i$. In this way, it can ensure the robust stability of certain discrete-time MJLS with actuator saturation. In the following, we provide three lemmas that are used in this paper.

$\textbf{Lemma 1.}$ ^[19] Let $P_i\in {\bf R}^{n\times n}$, $P_i>0$ and $\rho_i(k)>0$. An ellipsoid $\Omega(P_i, \rho_i(k))=\left\{x_i\in {\bf R}^{n_x}:x_i^{\rm T}P_ix_i\le\rho_i(k)\right\}$ is inside $\phi{(F_i)}$ if and only if

\begin{align}\notag f_{i, \alpha}^{\rm T}P_i^{-1}f_{i, \alpha}\le\rho_i^{-1}(k), ~~~\alpha =1, \ldots, m_i. \end{align}

Let $\Xi_i$ be the set of $m_i\times m_i$ diagonal matrices whose diagonal elements are either $0$ or $1$. Then, there are $2^{m_i}$ elements in $\Xi_i$. Suppose that elements of $\Xi_i$ are labeled as $E_{i, \beta}$ with $\beta$ $\in$ $\{1, \ldots, 2^{m_i}\}$. Denote $E^-_{i, \beta}=I-E_{i, \beta}$, obviously $E^-_{i, \beta}$ is also an element of $\Xi_i$ if $E_{i, \beta}\in \Xi_i$.

$\textbf{Lemma 2.}$ ^[20] Let $F_i, H_i\in {\bf R}^{m_i\times {n_x}}$ be given. Supposing that $|h^{\rm T}_{i, \alpha}x_i|\leq 1$ for all $\alpha=1, \ldots, m_i$, then we have

\begin{align}\notag \sigma(F_ix_i)\in \text{Co}\left\{E_{i, \beta}F_ix_i+E^-_{i, \beta}H_ix_i, ~~~\beta=1, \ldots, 2^{m_i}\right\}, \end{align}

where $h^{\rm T}_{i, \alpha}$ is the $\alpha$th row of the matrix $H_i$. In this way, the saturated feedback $\sigma(F_ix_i)$ can be transformed into the convex hull of a group of linear feedbacks.

$\textbf{Lemma 3.}$ ^[11] Let $A(\lambda), B(\lambda)$ be polytopic uncertain as $A(\lambda)=\sum_{i=1}^L\lambda_iA_i$, $B(\lambda)=\sum_{i=1}^L\lambda_iB_i$, and let $F$ be a constant vector. $A(\lambda), B(\lambda)$ and $F$ are of appropriate dimensions. Then, $S(\lambda)=A(\lambda)+B(\lambda)F$ is polytopic uncertain.

III. MAIN RESULTS

A. Controller Design

In this section, the stability of the distributed polytopic certain discrete-time MJLS subject to actuator saturation is discussed, and the corresponding distributed MPC controllers are designed. To discuss a distributed MPC problem under the state feedback control, an upper bound on a robust performance objective is minimized for every subsystem defined by (4)-(6). The min-max problem on the objective is solved as follows:

\begin{align}\label{7} \min_{\substack{ u_i(k+n|k)\\ n\geq0\\ i=1, \ldots, M\\ }}\quad\max_{\substack{ [A^{(\xi)}_i(r_{k+n|k}), B^{(\xi)}_i(r_{k+n|k}), \\ C^{(\xi)}_i(r_{k+n|k})]\in{\Omega(r_k)}\\ i=1, \ldots, M, ~r_{k+n|k}\in{ M}}}~J_i(k), \end{align}

(7)

with

\begin{align*}J_i(k)={{\rm E}}&\bigg[\sum_{n=0}^\infty\Big(\| x_i(k+n|k)\|^2_{Q_i{(r_{k+n|k})}}\\[1mm] &+\sum_{i=1}^M\|{\bar{u}}_i(k+n|k)\|^2_{R_i(r_{k+n|k)}}\Big)\bigg], \end{align*}

where $Q_i> 0$ and $R_i > 0$ are symmetric weighting matrices. With the overall cost function taken into consideration in (7), it becomes a cooperative distributed MPC algorithm subject to actuator saturation. Applying the state feedback law (6) to (4) results in the following closed-loop system:

\begin{align*} &x_i(k+n+1|k)=A^{(\xi)}_i(r_{k+n|k})x_i(k+n|k)\notag\\[1mm] &\qquad + B^{(\xi)}_i(r_{k+n|k})\sigma(F_i(k, r_{k+n|k})x_i(k+n|k))\notag \end{align*} \begin{align}\label{8} & +\sum_{j=1, j\neq i}^M\sigma(F_j(k, r_{k+n|k})x_j(k+n|k)), \end{align}

(8)

where the term $\sum_{j=1, j\neq i}^M\sigma(F_j(k, r_{k+n|k})x_j(k+n|k))$ represents the effects of neighbor subsystems. The coupled inputs are saturated and cannot be handled directly. In order to derive the closed-loop system, we introduce an unsaturated feedback law $\Gamma_j(k, r_{k+n|k})$, which is restructured as follows:

\begin{align}\label{9} \Gamma_{j, \alpha}(k, r_{k+n|k})= \begin{cases} F_{j, \alpha}(k, r_{k+n|k}), \\[1mm] ~~~~~~~ |F_{j, \alpha}(k, r_{k+n|k})x_j(k+n|k)|\le1, \\[4mm] \dfrac{F_{j, p}(k, r_{k+n|k})}{\left|F_{j, p}(k, r_{k+n|k})x_j(k)\right|}, \\[1mm] ~~~~~~~ |F_{j, \alpha}(k, r_{k+n|k})x_j(k+n|k)|>1, \end{cases} \end{align}

(9)

where $\alpha=1, \ldots, m_i$. $\Gamma_{j, \alpha}(k, r_{k+n|k})$ stands for the $\alpha$-th row of $\Gamma_{j}(k, r_{k+n|k})$, $F_{j, \alpha}(k, r_{k+n|k})$ refers to the $\alpha$-th row of $F_{j}(k, r_{k+n|k})$. Then, we have

\begin{align}\label{10} &\bar{u}_j(k+n|k)=\sigma(u_j(k+n|k))\notag\\[1mm] &\qquad=\sigma(F_{j}(k, r_{k+n|k})x_j(k+n|k))\notag\\[1mm] &\qquad=\Gamma_{j}(k, r_{k+n|k})x_j(k+n|k). \end{align}

(10)

In the sequel of the paper, for presentation convenience, we denote each possible $r_{k+n|k}=p~(p\in {M}).$ The closed-loop system (8) by applying unsaturated-input couplings can be rewritten as follows:

\begin{align} x_i(k~+&~n+1|k)\notag\\[1mm] =&\ A^{(\xi)}_i(p)x_i(k+n|k)\notag\\[1mm] & + B^{(\xi)}_i(p)\left(E_{i, \beta}F_i(k, p)+E^-_{i, \beta}H_i(k, p)\right)x_i(k+n|k)\notag\\[1mm] & +\sum_{j=1, j\neq i}^MB^{(\xi)}_j(p)\Gamma_j(k, p)x_i(k+n|k)\notag\\[1mm] =&\ \left[{\hat{A}}^{(\xi)}_i(p)+\hat{B}^{(\xi)}_i(p)\right] \times x_i(k+n|k), \end{align}

(11)

where ${\hat{A}}^{(\xi)}_i(p)=A^{(\xi)}_i(p)+\sum_{j=1, j\neq i}^MB^{(\xi)}_i(p)\Gamma_j(k, p), $ $\hat{B}^{(\xi)}_i(p)$ $=$ $B^{(\xi)}_i(p)(E_{i, \beta}F_i(k, p)+E^-_{i, \beta}H_i(k, p))$. By Lemma 3, we can conclude ${\hat{A}}^{(\xi)}_i(p)$ is a polytopic uncertainty.

Consider the following quadratic function:

\begin{align}\notag V_i(k+n|k)=x^{\rm T}_i(k+n|k)P_i(k, p)x_i(k+n|k), \end{align}

where $i=1, \ldots, M$, $n\ge0$, $p\in{M}$. According to the stochastic switching of the discussed system (4) in certain given modes, the system should be a stochastic one and then the quadratic functions $V_i(k+n|k)$ should be stochastic. Thus, during the stability analysis, we should take the expectation of the quadratic functions $V_i(k+n|k)$. To derive an upper bound on (7), we impose the following robust stability constraints

\begin{align*} &{\rm E}\big[V_i(k+n+1|k)-V_i(k+n|k)\big]\notag\\[1mm] & \qquad\le-{\rm E}\bigg[x^{\rm T}_i(k+n|k)Q_i(p)x_i(k+n|k)\notag \end{align*} \begin{align}\label{12} &+\bar{u}^{\rm T}_i(k+n|k)R_i(p)\bar{u}_i(k+n|k)\notag\\ &+\sum_{j=1, j\neq i}^M\bar{u}^{\rm T}_j(k+n|k)R_j(p)\bar{u}_j(k+n|k)\bigg] \end{align}

(12)

$\textbf{Theorem 1\bf.}$ For each subsystem $i$, robust stability constraints (12) are satisfied if there exist a positive scalar $\rho_i(k)$, positive definite matrices $W_i(k, p)=\rho_i(k)P_i^{-1}(k, p)>0$, and any matrices $Y_i(k, p)=F_i(k, p)W_i(k, p)$, such that the following LMIs

\begin{align}\label{13} &\!\!\!\! {\left[\begin{array}{ccccccc} {-W_i(k, p)} & {* } \\ {\Phi_i(k, p)} &\!\!\!\!\!\! {-p^{-1}(p, 1)W_i(k, 1)} \\ {\vdots} & {\vdots}\\ {\Phi_i(k, p)} & {0} \\ {\Psi_i(p)^{\frac{1}{2}}W_i(k, p)} & {0} \\ {R_i^{\frac{1}{2}}(p)\Big(E_{i, \beta}Y_i(k, p)+E^-_{i, \beta}Z_i(k, p)\Big)} & {0} \end{array}\right. } \notag\\ \notag\\ &\left.\qquad\ \ \ \ \begin{array}{ccccc} \!\!\!\!\!{\cdots} & {*} & {*} & {*} \\ \!\!\!\!\!{\cdots} & {*} & {*} & {*} \\ {\ddots} & {\vdots} & {\vdots} & {\vdots} \\ {\!\!\!\!\!\cdots} & {-p^{-1}(p, S)W_i(k, S)} & {*} & {*} \\ \!\!\!\!\!{\cdots} & {0} & {-\rho_i(k)I} & {*} \\ \!\!\!\!\!{\cdots} & {0} & {0} & {-\rho_i(k)I}\end{array}\right] \notag\\ &<0, \notag\\ \end{align} \begin{align}

(13)

\begin{bmatrix}\label{14} -1&*\\ x_i(k)&-W_i(r_{k})\\ \end{bmatrix}\leq0 \end{align}

(14)

hold. Then the gain matrix

\begin{align*} &Y_i(k, p)=F_i(k, p)W_i(k, p), \\[1mm] &Z_i(k, p)=H_i(k, p)W_i(k, p), \end{align*}

and the upper bound of the optimization problem is $\rho_i(k)$, where

\begin{align*}&\Phi_i(k, p)={\hat{A}}^{(\xi)}_i(p)W_i(k, p)\\[1mm] &\qquad+B^{(l)}_i(p)\left(E_{i, \beta}Y_i(k, p)+ E^-_{i, \beta}Z_i(k, p)\right), \\[2mm] &\Psi_i(p)=Q(p)+\sum_{j=1, j\neq i}^M\Gamma^{\rm T}_j(k, p)R_j(p)\Gamma_j(k, p).\end{align*}

$\textbf{Proof}{\bf.}$ The difference of $V_i(k+n|k)$ along system (11) can be calculated as follows:

\begin{align*} {\rm E}\big[\Delta& V_i(k+n|k)\big]={\rm E}\big\{[V_i(k+n+1|k)]-V_i(k+n|k)\big\}\notag\\[1mm] =&\ {\rm E}\big\{[x_i^{\rm T}(k+n+1|k)P_i(k, r_{k+n+1|k})x_i(k+n+1|k)]\big\}\notag\\[1mm] &-x_i^{\rm T}(k+n|k)P_i(k, p)x_i(k+n|k)\notag\\[1mm] =&\ x_i^{\rm T}(k+n|k)\bigg\{\Big[{\hat{A}}^{(\xi)}_i(p)+B^{(\xi)}_i(p)\big(E_{i, \beta}F_i(k, p)\notag\\[1mm] &+E^-_{i, \beta}H_i(k, p)\big)\Big]^{\rm T}{\rm E}\big[P_i(k, r_{k+n+1|k})\big] \notag\\[1mm] &\times\Big[{\hat{A}}^{(\xi)}_i(p) +B^{(\xi)}_i(p)\big(E_{i, \beta}F_i(k, p)+E^-_{i, \beta}H_i(k, p)\big)\Big]\notag\end{align*}\begin{align} &\hspace{-35mm}-P_i(k, p)\bigg\}x_i(k+n|k). \end{align}

(15)

Furthermore, (12) can be reformulated as follows:

\begin{align}\notag {\rm E}[\Delta V_i&(k+n|k)]\notag\\[1mm] \le& -{\rm E}\bigg[x^{\rm T}_i(k+n|k)Q_i(p)x_i(k+n|k)\notag\\[1mm] & +\bar{u}^{\rm T}_i(k+n|k)R_i(p)\bar{u}_i(k+n|k)\notag\\[1mm] & +\sum_{j=1, j\neq i}^M\bar{u}^{\rm T}_j(k+n|k)R_j(p)\bar{u}_j(k+n|k)\bigg]. \end{align}

(16)

By applying (9) into (16), we can get following inequality:

\begin{align}\label{17} {\rm E}[\Delta& V_i(k+n|k)]\notag\\[1mm] \le&-{\rm E}\bigg\{-x^{\rm T}_i(k+n|k)\Big[Q_i(p)+\big(E_{i, \beta}F_i(k, p)\notag\\[1mm] & +E^-_{i, \beta}H_i(k, p)\big)^{\rm T}R_i(p)\left(E_{i, \beta}F_i(k, p)+E^-_{i, \beta}H_i(k, p)\right)\notag\\[1mm] & +\sum_{j=1, j\neq i}^M\Gamma^{\rm T}_j(k, p)R_j(p)\Gamma_j(k, p)\Big]x_i(k+n|k)\bigg\}. \end{align}

(17)

According to (2), it is clear that

$${\rm E}[P_i(k, r_{k+n+1|k})]=\sum_{h=1}^S p(g, h)P_i(k, h).$$

Then combining (15) and (17), we can obtain the following inequality:

\begin{align}\label{18} &x^{\rm T}_i(k+n|k)\notag\\[1mm] &\quad \times\bigg\{\bigg[{\hat{A}}^{(\xi)}_i(p)+B^{(\xi)}_i(p)\left(E_{i, \beta}F_i(k, p)+ E^-_{i, \beta}H_i(k, p)\right)\bigg]^{\rm T}\notag\\[1mm] &\quad \times\sum_{h=1}^Sp(g, h)P_i(k, h)\notag\\[1mm] &\quad\times\bigg[{\hat{A}}^{(\xi)}_i(p)+B^{(\xi)}_i(p)\left(E_{i, \beta}F_i(k, p)+E^-_{i, \beta}H_i(k, p)\right)\bigg]\notag\\[1mm] &\quad -P_i(k, p)+Q_i(p)+\sum_{j=1, j\neq i}^M\Gamma^{\rm T}_j(k, p)R_j(p)\Gamma_j(k, p)\notag\\[1mm] &\quad+\left(E_{i, \beta}F_i(k, p)+E^-_{i, \beta}H_i(k, p)\right)^{\rm T}R_i(p)\big(E_{i, \beta}F_i(k, p)\notag\\[1mm] &\quad +E^-_{i, \beta}H_i(k, p)\big)\bigg\}x_i(k+n|k)\le 0. \end{align}

(18)

According to (5), the above inequality is satisfied if and only if

\begin{align*} &\bigg[{\hat{A}}^{(l)}_i(p)+B^{(l)}_i(p)\left(E_{i, \beta}F_i(k, p)+E^-_{i, \beta}H_i(k, p)\right)\bigg]^{\rm T}\notag\\[1mm] &\quad\times\sum_{h=1}^Sp(g, h)P_i(k, h)\notag\\[1mm] &\quad\times\bigg[{\hat{A}}^{(l)}_i(p)+B^{(l)}_i(p)\left(E_{i, \beta}F_i(k, p)+E^-_{i, \beta}H_i(k, p)\right)\bigg]\notag\end{align*}\begin{align} &\quad -P_i(k, p)+Q_i(p)+\sum_{j=1, j\neq i}^M\Gamma^T_j(k, p)R_j(p)\Gamma_j(k, p)\notag\\ &\quad +\left(E_{i, \beta}F_i(k, p)+E^-_{i, \beta}H_i(k, p)\right)^{\rm T}R_i(p)\big(E_{i, \beta}F_i(k, p)\notag\\ &\quad +E^-_{i, \beta}H_i(k, p)\big)\le 0. \end{align}

(19)

Pre- and post-multiplying (19) with $P^{-1}_i(k, p)$ result in

\begin{align}\label{20} &\bigg[{\hat{A}}^{(l)}_i(p)P^{-1}_i(k, p)\notag\\[1mm] &~ +B^{(l)}_i(p)\left(E_{i, \beta}F_i(k, p)+E^-_{i, \beta}H_i(k, p)\right)P^{-1}_i(k, p)\bigg]^{\rm T}\notag\\[1mm] &~ \times \sum_{h=1}^Sp(g, h)P_i(k, h)\bigg[{\hat{A}}^{(l)}_i(p)P^{-1}_i(k, p)\notag \\[1mm] &~+ B^{(l)}_i(p)\left(E_{i, \beta}F_i(k, p)+E^-_{i, \beta}H_i(k, p)\right)P^{-1}_i(k, p)\bigg]\notag\\[1mm] &~ -P^{-1}_i(k, p)+P^{-1}_i(k, p)Q_i(p)P^{-1}_i(k, p)\notag\\[1mm] &~+\sum_{j=1, j\neq i}^M\big(\Gamma_j(k, p)P^{-1}_i(k, p)\big)^{\rm T} R_j(p)\big(\Gamma_j(k, p)P^{-1}_i(k, p)\big)\notag\\[1mm] &~+\!\left(E_{i, \beta}F_i(k, p)P^{-1}_i(k, p)+E^-_{i, \beta}H_i(k, p)P^{-1}_i(k, p)\right)^{\rm T}\!R_i(p)\notag\\[1mm] &~\times\big(E_{i, \beta}F_i(k, p)P^{-1}_i(k, p)+E^-_{i, \beta}H_i(k, p)P^{-1}_i(k, p)\big)\le 0. \end{align}

(20)

Then, inequalities (20) can be further simplified by multiplying with $\rho_i(k)$

\begin{align}\label{21} &\bigg[{\hat{A}}^{(l)}_i(p)W_i(k, p)+B^{(l)}_i(p)\left(E_{i, \beta}Y_i(k, p)+E^-_{i, \beta}Z_i(k, p)\right)\bigg]^{\rm T}\notag\\[1mm] &~\times\sum_{h=1}^S p(g, h)W_i^{-1}(k, h)\notag\\[1mm] &~\times\!\bigg[{\hat{A}}^{(l)}_i(p)W_i(k, p)\!+\!B^{(l)}_i(p)\!\left(E_{i, \beta}Y_i(k, p)+E^-_{i, \beta}Z_i(k, p)\right)\!\bigg]\notag\\[1mm] &~-W_i(k, p)+W_i^{\rm T}(k, p)\rho_i^{-1}(k)Q_i(p)W_i(k, p)\notag\\[1mm] &~+W_i^{\rm T}(k, p)\rho_i^{-1}(k)\sum_{j=1, j\neq i}^M\Gamma^{\rm T}_j(k, p)R_j(p)\Gamma_j(k, p)W_i(k, p)\notag\\[1mm] &~+\left(E_{i, \beta}Y_i(k, p)+E^-_{i, \beta}Z_i(k, p)\right)^{\rm T}\rho_i^{-1}(k)R_i(p) \notag\\[1mm] &~\times\left(E_{i, \beta}Y_i(k, p)+E^-_{i, \beta}Z_i(k, p)\right)\bigg\}\le0. \end{align}

(21)

By using the Schur complement, if LMI (13) holds, the above inequality will be satisfied.

In the following, based on the robust stability constraints (12), an upper bound on the worst-case infinite horizon (7) can be obtained. Firstly, taking the expected values on the both sides of (12) conditional on the information available at time instant $k$, the following inequality holds

\begin{align*} &\hspace{-16mm}{\rm E}\Big[V_i(k+n+1|k)-V_i(k+n|k)\Big]\notag \end{align*}\begin{align*} \le& -{\rm E}\Big[x^{\rm T}_i(k+n|k)Q_i(p)x_i(k+n|k)\notag\\[1mm] & +\bar{u}^{\rm T}_i(k+n|k)R_i(p)\bar{u}_i(k+n|k)\notag\\[1mm] & +\sum_{j=1, j\neq i}^M\bar{u}^{\rm T}_j(k+n|k)R_j(p)\bar{u}_j(k+n|k)\Big].\notag \end{align*}

Then, summing up the above inequality from $n=0$ to $\infty$ on both sides, we have

\begin{align*}\notag {\rm E}\Big[V_i&(k)-V_i(\infty|k)\Big]\notag\\[1mm] \ge&\ {\rm E}\bigg\{\sum_{n=0}^\infty\Big[x^{\rm T}_i(k+n|k)Q_i(p)x_i(k+n|k)\notag\\[1mm] & +\bar{u}^{\rm T}_i(k+n|k)R_i(p)\bar{u}_i(k+n|k)\notag\\[1mm] & +\sum_{j=1, j\neq i}^M\bar{u}^{\rm T}_j(k+n|k)R_j(p)\bar{u}_j(k+n|k)\Big]\bigg\}.\notag \end{align*}

Because the robust performance objective function $J_i(k)$ should to be finite, we must have $x_i(\infty|k)=0$, which leads to $V_i(\infty|k)=0$,

\begin{align*}\notag J_i(k)=&\ {{\rm E}}\bigg[\sum_{n=0}^\infty\Big(\| x_i(k+n|k)\|^2_{Q_i{(r_{k+n|k})}}\notag\\[1mm] & +\sum_{i=1}^M\|{\bar{u}}_i(k+n|k)\|^2_{R_i(r_{k+n|k)}}\Big)\bigg]\notag\\[1mm] \le&\ {\rm E}\big[V_i(k)\big], \notag \end{align*}

that is

\begin{align*}\notag J_i(k)\le{\rm E}\big[V_i(k)\big]=x_i^{\rm T}(k)P_i(k, r_k)x_i(k).\notag \end{align*}

By defining an upper bound $x_i^{\rm T}(k)P_i(k, r_k)x_i(k)\le\rho_i(k)$, then by Schur complement it can be written as (14) and we can obtain

\begin{align}\notag J_i(k)\le\rho_i(k). \end{align}

We can see that $\rho_i(k)$ is an upper bound on the expected cost function.

$\square$

Using Lemma 1, condition $|h^{\rm T}_{i, \alpha}(k)x_i(k)|\le 1$ can be written as follows:

\begin{align} \begin{bmatrix}\label{22} 1&z^{\rm T}_{i, \alpha}(k)\\[1mm] z_{i, \alpha}(k)&W_i(k)\\ \end{bmatrix}\ge 0, ~~\alpha=1, \ldots, m_i, \end{align}

(22)

where $z_{i, \alpha}(k)=W_i(k)h_{i, \alpha}(k)$.

Next, we will minimize upper bound to approximately minimize the worst-case infinite horizon expected cost function. Now, we are ready to present the unconstrained distributed MPC of MJLSs with polytopic uncertainties in terms of a minimization problem at each time instant $k$ as follows:

\begin{align}\label{23} \min_{W_i(k, p), Y_i(k, p), p\in M}~\rho_i(k)\quad {\rm s.t.}\quad (13), ~(14), ~(22). \end{align}

(23)

In the following theorem, we will discuss the feasibility and stability of (23).

$\textbf{Theorem 2\bf.}$ Consider the MJLSs with actuator saturation described by (3)-(5). If there is a feasible solution to optimization problem (23) at time instant $k$ for initial state $x(k)$ and initial mode $r_k$, there will also exist a feasible solution at time instant $t \geq k$; and the distributed MPC control law $F_i(k, r_k)$ $=$ $Y_i(k, r_k)W^{-1}_i(k, r_k)$ based on (23) guarantees closed-loop stability of the system in the mean-square sense.

$\textbf {Proof.}$ The proof of this theorem is similar to ^[11] and ^[15]. If the solution to problem (23) at time $k$ is feasible, then it is also feasible at all future sampling steps $k+n$, $n>0$. This is because the only constraint that depends on the states in problem (23) is constraint (14), i.e., $-1+x^{\rm T}_i(k)W_i(r_k)x_i(k)$ $\leq$ $0$ , where the states are given by $x_i(k + n)=$ $[{\hat{A}}^{(\xi)}_i(p)$ $+$ $\hat{B}^{(\xi)}_i(p)] x_i(k+n-1)$. This constraint can be feasible by using the definition of invariant set that is satisfied at time $k$ (see ^[11] and ^[15]).

$\square$

B. Iteration Algorithm

Notice that there are couplings in LMIs (13) when solving problem (23). We propose an iterative distributed MPC algorithm in this section.

Iterative distributed MPC algorithm subject to actuator saturation of $i$-th subsystem is as follows:

${{\bf Step 1.}}$ At time step $k=0$, set an initial feedback law $F_i(0, r_k)$, $i=1, 2, \ldots, M$.

${{\bf Step 2.}}$ At time step $k$, all subsystems exchange their local states measurements and initial feedback law $F_i(k, r_k)$ via communication, set iteration $t=1$ and $F_i(k, r_k)=F_i(0, r_k)$.

${{\bf Step 3.}}$ All subsystem solve LMI problems (23) in parallel to obtain the optimal $Y^{(t)}_i(k, r_k)$ and $W^{(t)}_i(k, r_k)$ to estimate the optimal feedback laws $F^{(t)}_i(k, r_k)=Y^{(t)}_i(k, r_k)$ $\times$ $W^{-1(t)}_i(k, r_k)$. Each subsystem checks the convergence with a specified error tolerance $\eta_i$ for all the feedback laws.

\begin{align*}\notag \|F^{(t)}_i(k, r_k)-F^{(t-1)}_i(k, r_k)\|\le\eta_i, ~~i=1, 2, \ldots, M. \end{align*}

If the convergence condition or $t=t_{\rm max}$ is satisfied, current $F_i^{(t)}(k, r_k)$ is taken as the optimal feedback law; Otherwise, update the initial feedback laws with $F_i(k, r_k)=F_i^{(t)}(k, r_k)$ and set $t=t+1$, exchange the solutions with other subsystems and repeat Step 3; ${\bf Step 4.}$ The optimal scheme $u_i(k)=F_i(k, r_k)x_i(k)$ is applied to the corresponding subsystems. Set the time interval $k=k+1$ and go back to Step 2.

Iteration algorithm is solved in parallel for all subsystems, thus the feedback laws can be obtained at a same time. In this way, the computational time can be saved by distributed MPC algorithm.

$\textbf{Lemma 4.}$^[15] For cooperative distributed MPC algorithm, the performance indices $\rho_1(k), \rho_2(k), \ldots, $ $\rho_M(k), $ of subsystems will converge to the centralized problem with the increase of iteration times, which implies

\begin{align*}\notag \rho_1^{(t)}(k)=\rho_2^{(t)}(k)=\cdots=\rho_M^{(t)}(k)=\rho(k), \notag \end{align*}

where $\rho(k)$ is the upper bound of centralized MPC at time interval $k$.

$\textbf{Proof.}$ The proof is similar to ^[15]. For subsystem $i$,

\begin{align*}\notag \rho^{(t)}_i&= \rho_i(F^{(t-1)}_1, \ldots, F^{(t)}_i, \ldots, F^{(t-1)}_M)\\[1mm] & =\min_{F^{(t)}_i} \rho_i(F^{(t-1)}_1, \ldots, F^{(t-1)}_i, \ldots, F^{(t-1)}_M), \\[1mm] &\qquad\qquad\qquad\qquad\qquad\quad i\in\{1, \ldots, M\}. \end{align*}

Similarly for subsystem $j$,

\begin{align*}\notag \rho^{(t)}_j&=\rho_j(F^{(t-1)}_1, \ldots, F^{(t)}_j, \ldots, F^{(t-1)}_M)\\[1mm] &=\min_{F^{(t)}_j}\rho(F^{(t-1)}_1, \ldots\!, F^{(t-1)}_j, \ldots, F^{(t-1)}_M), \\[1mm] &\qquad\qquad\qquad\quad\ j\in \{1, \ldots M\}, \ j\neq i. \end{align*}

Then from the convexity of problem (13) and for any $i$ and $j$ pair of subsystems, $\rho^{(t)}_i\leq \rho^{(t-1)}_j$ and $\rho^{(t)}_j\leq \rho^{(t-1)}_i$. Thus, the minimization index of worst-case cost function $\rho_i$ continues to decrease with the increase of iteration times until both inequalities become equalities. Since the minimizations given above are convex and each is leading to a global optimum, we can consequently obtain $\rho^{(t)}_i=\rho^{(t)}_j=\rho$, $\forall i, j\in \{1, \ldots, M\}$, $i\neq j$.

$\square$

${\bf Remark 4.}$ Lemma 4 tells us that the performance by distributed MPC algorithm can be as good as the performance by centralized MPC algorithm within several iterations. However, the more computational time will be cost with the increase of iteration times. There is a tradeoff between computational time and the performance. In terms of computational time, the distributed MPC presented in this paper can outperform centralized MPC with little performance loss when terminated after a few iterations.

IV. NUMERICAL EXAMPLES

${\bf Example 1.}$ Consider the system as follows:

\begin{align}\label{24} \begin{cases} x(k+1)=A^{(\xi)}(r_k)x(k)+B^{(\xi)}(r_k)\sigma(u(k)), \\[1mm] y(k)=C^{(\xi)}(r_k)x(k), \end{cases} \end{align}

(24)

where the system matrices are represented by the following two vertices:

Mode 1.

\begin{align*}\notag &A^{(1)}(r1)= \begin{bmatrix} 0.9&0\\[1mm] 0.75&0.4 \end{bmatrix}, \\[2mm] &B^{(1)}(r1)= \begin{bmatrix} 0.04&0.435\\[1mm] 0.0356&0.94 \end{bmatrix}, \quad\notag\\[2mm] &A^{(2)}(r1)= \begin{bmatrix} 0.10&0\\[1mm] 0&0.1 \end{bmatrix}, \\[2mm] &B^{(2)}(r1)= \begin{bmatrix} 0&0.01\\[1mm] 0&0.01 \end{bmatrix}.\notag \end{align*}

Mode 2.

\begin{align*}\notag &A^{(1)}(r2)= \begin{bmatrix} 0.884&0\\[1mm] 0.3&0.1 \end{bmatrix}, \\[2mm] &B^{(1)}(r2)= \begin{bmatrix} 0.092&0\\[1mm] 0&0.09 \end{bmatrix}, \quad\notag\\[2mm] &A^{(2)}(r2)= \begin{bmatrix} 0.501&0\\[1mm] 1.5&0.1 \end{bmatrix}, \end{align*}\begin{align*} &B^{(2)}(r2)= \begin{bmatrix} 0.2&0.01\\[1mm] 0&0.8 \end{bmatrix}.\notag \end{align*}

where $r_1$ and $r_2$ denote Mode 1 and Mode 2, respectively. Assume that the transition probability matrices with polytopic uncertainties can be denoted as

\begin{align}\notag {P}= \begin{bmatrix} 0.95&0.05\\[1mm] 0.6&0.4 \end{bmatrix}. \end{align}

We decompose the global system into two subsystems, each of which has one control input. The actuator saturation is described as $\sigma\left(u_1(k)\right)={\mbox{sgn}}\left(u_1(k)\right)\min \left\{1, |u_1(k)|\right\}$ and $\sigma\left(u_2(k)\right)={\mbox{sgn}}\left(u_2(k)\right)\min \{1, |u_2(k)|\}$. The proposed distributed MPC algorithm is used to control these two subsystems with initial state $x_0=[-1, 1]^{\rm T}$ and weighting matrices $Q_1=Q_2=I_2$, $R_1=R_2=1$. The performance of distributed algorithm is compared with centralized MPC algorithm under the same weighting matrices.

Fig. 1 shows the invariant sets obtained by different algorithms at the first sampling step. It can be seen that the proposed distributed MPC algorithm reaches the same invariant set as centralized MPC algorithm after three iterations, which validates the result of Lemma 4. Fig. 2 (a) presents the upper-bound trends of centralized MPC, subsystems 1 and 2, respectively, from which we find that only after three iterations the upper bound of the subsystems 1 and 2 almost converge to the upper bound of centralized MPC. From Fig. 2 (b) we can see the upper bound either by centralized MPC or by distributed MPC will converge to zero when the time step $k$ tends to infinite. The comparison of upper bound and computational time between centralized MPC (CMPC) and distributed MPC (DMPC) at the first sampling step is shown in Table I. The distributed MPC algorithm is applied with one single iteration, two iterations and three iterations. From Table I, we can see the distributed MPC algorithm only takes 0.92 s after three iterations, which is less than 1.23 s by centralized MPC. Distributed MPC outperform centralized MPC with little performance loss, which validate the results described in Remark 4.

V. CONCLUSION

Thinking of the uncertainty and abrupt changes of system parameters, we present a distributed MPC for MJLSs. By decomposing the system into several subsystems to design the distributed controllers, it can decrease the system complexity and computation cost. In each subsystem, a distributed MPC strategy is proposed. The distributed controller is obtained by solving an LMI optimization problem and an iteration algorithm is developed for making cooperation among subsystems. Due to the consideration of actuator saturation, the main results made in this paper are more practical than those with ideal network connection. In the end, an example illustrates the effectiveness of the proposed distributed MPC algorithm.