Polyhedral Feasible Set Computation of MPC-Based Optimal Control Problems
 IEEE/CAA Journal of Automatica Sinica  2018, Vol. 5 Issue(4): 765-770 PDF
Polyhedral Feasible Set Computation of MPC-Based Optimal Control Problems
Lantao Xie1, Lei Xie1, Hongye Su1, Jingdai Wang2
1. State Key Laboratory of Industrial Control Technology, Institute of Cyber-Systems and Control, Zhejiang University, Hangzhou 310027, China;
2. College of Chemical and Biological Engineering, Zhejiang University, Hangzhou 310027, China
Abstract: Feasible sets play an important role in model predictive control (MPC) optimal control problems (OCPs). This paper proposes a multi-parametric programming-based algorithm to compute the feasible set for OCP derived from MPC-based algorithms involving both spectrahedron (represented by linear matrix inequalities) and polyhedral (represented by a set of inequalities) constraints. According to the geometrical meaning of the inner product of vectors, the maximum length of the projection vector from the feasible set to a unit spherical coordinates vector is computed and the optimal solution has been proved to be one of the vertices of the feasible set. After computing the vertices, the convex hull of these vertices is determined which equals the feasible set. The simulation results show that the proposed method is especially efficient for low dimensional feasible set computation and avoids the non-unicity problem of optimizers as well as the memory consumption problem that encountered by projection algorithms.
Key words: Feasible set computation     linear matrix inequalities (LMIs)     model predictive control (MPC)     polyhedron
Ⅰ. INTRODUCTION

Model predictive control (MPC) is a widely used method in process industries and fields such as automobile, energy, environment, aerospace and medical treatment, etc [1], [2] because of its ability to effectively handle the complex dynamics of systems with multiple inputs and outputs, system constraints, and conflicting control objectives. A key role played in a MPC framework is the so-called feasible set, i.e. the largest subset of the state space for which there exists a control action satisfying all the constraints [3]. MPC-based controllers will finally result in solving a standard optimal control problems (OCPs) such as linear programming (LP), quadratic programming (QP), semidefinite programming (SDP), dynamic programming (DP), etc. If the corresponding OCP is designed to be stable or recursively feasible, once the initial state is in the feasible set, there always exists an optimal control action at each following sample time.

As stated in [4], a larger feasible set usually means that the corresponding MPC algorithm is less conservative. Thus, when comparing conservativeness of different MPC algorithms, feasible sets are often compared, such as in [5], [6]. Feasible sets can also be used to modify initial work states to ensure feasibility, and its shape description plays an important role in the effectiveness of many of the approaches for approximate explicit MPC [3]. By restricting the state of the second step in the prediction horizon in the feasible set and computing each set iteratively, the algorithm proposed in [7] obtains the recursively feasible set and guarantees recursive feasibility. Even with the consideration of techniques used for stability, feasible sets are closely related with the prediction horizon and system constraints; generally longer horizons and looser system constraints result in larger feasible sets.

Oval-shaped and polyhedral feasible sets are two commonly used convex sets in MPC algorithms. Ellipsoid can be described as a weighted 2-norm and its size is related with some measure of the weighted matrix. By performing a maximum volume optimization in terms of the weighted matrix, one can obtain the largest ellipsoid, i.e., the oval-shaped feasible set, such as that shown in [8]. A polyhedron is the intersection of finitely many halfspaces. Computation of a polyhedron is much more complicated than of a ellipsoid. The standard approach to compute polyhedral feasible set is orthogonal projection [9]. However, the orthogonal projection needs the polyhedron to be described by its vertices (V-representation) or the intersection of halfspaces (H-representation), which cannot be applied when the polyhedron is written in linear matrix inequality (LMI) form. Moreover, orthogonal projection often turns out to be computationally intractable in high spatial dimensions. This often happens when the feasible set is computed for MPC-based OCPs with long prediction horizons or with large number of slack variables such as that in some robust MPC [10] and stochastic MPC algorithms [7]. Alternatives to compute feasible set involve parametric linear programming (PLP) methods [11] and set relation methods [3]. However, PLP methods may lead to numerical difficulties due to the non-unicity of the optimizers [12], while set relation methods are only efficient for OCPs without slack variables.

Based on the final OCP derived from MPC-based algorithms, this paper proposes a novel multi-parametric programming problem where the final OCP involves both spectrahedron (represented by LMI) and polyhedral (represented by a set of inequalities) constraints. According to the geometrical meaning of an inner product of vectors, the maximum length of the projection vector from the feasible set to a unit spherical coordinates vector is computed and the optimal solution is proven to be one of the vertices of the feasible set. After computing all the vertices of the feasible set, the convex hull of these vertices is computed which is equivalent to the feasible set. The simulation results show that the proposed method is especially efficient for low dimensional feasible set computation and avoids the non-unicity problem of optimizers, as well as the memory consumption problem encountered by projection algorithms.

$\mathit{Notations:}$ The sets of reals is denoted by $\mathbb{R}$. Given two integers $a$ and $b$ where $a < b$, $N_{[a, b]}=[a, a+1, a+2, \ldots, b]$ and given two sets $\mathbb{X} \subseteq \mathbb{R}^n$ and $\mathbb{Y} \subseteq \mathbb{R}^n$, the Minkowski set addition is defined by $\mathbb{X}\oplus\mathbb{Y}=\{x+y|x\in \mathbb{X}, y \in \mathbb{Y} \}$ and the Pontryagin set difference is then defined by $\mathbb{X}\ominus \mathbb{Y}=\{x|x\oplus \mathbb{Y} \subseteq \mathbb{X} \}$. Given set $\mathbb{Z} \subseteq \mathbb{R}^{n+m}$, its projection onto $\mathbb{R}^n$ is defined by $Proj_{x}(\mathbb{Z})=\{x\in \mathbb{R}^n| \exists y \in \mathbb{R}^m \text{ such that } (x, y) \in \mathbb{Z} \}$. Given matrix $Q_i$, diag$(Q_1, \ldots, Q_n)$ denotes the block diagonal matrix with matrix $Q_i$ on the main diagonal. ${E}(\cdot)$ denotes expectation. $|\vec{a}|$ denotes the length of vector $\vec{a}$. $[x]_i$ is the $i$-th element of $x$. $\| x \|$ denotes the Euclidean norm. The semidefinite matrix A is denoted by $A\succeq0$. Given two vectors $\vec{a}, \vec{b}$, the angle between $\vec{a}$ and $\vec{b}$ is expressed as $\angle(\vec{a}, \vec{b})$. conv$\{v_1, v_2, \ldots, v_n\}$ denotes the convex hull of points $(v_1, v_2, \ldots, v_n)$.

Ⅱ. PROBLEM SETUP

Consider the controlled system to be described as

 $$$\label{sys} x_{k+1}=(f+\triangle f)(x_k, u_k, d_k, w_k)$$$ (1)

where $x_k$ is system state and $u_k$ is system input, $\triangle f$ is model uncertainty, $d_k$ is deterministic noise and $w_k$ represents stochastic noise with proper dimensions. Constraints to the systems include expectation constraints which can be expressed as $H{E}(x)+G{E}(u)\ge 0$ [13], probability constraints (or chance constraints) ($Pro\{Hx+Gu\ge 0\}\ge1-\epsilon$) [14], and hard constraints ($Hx+Gu\ge 0$) [15]. By using linear matrix inequalities (LMI) or constraints tightening methods, the optimal control problems of different algorithms can be finally converted to a uniform uncertainty-free form:

 $\begin{eqnarray} &&J(x_k, \theta_k)=\min\limits_{\theta_k}: V(x_k, \theta_k) \nonumber\\ &&{\rm s.t.}\, \, H(x_k, \theta_k)\succeq0, G(x_k, \theta_k)\ge0 \end{eqnarray}$ (2)

where $\succeq$ denotes partial order [16] (i.e. $A\succeq B$ if and only if $A-B$ is positive semidefinite) and $H(x_k, \theta_k)$ is a square matrix. The symbol "$\ge$" represents the vector inequality or componentwise inequality [17] when $G(x_k, \theta_k)$ is a vector or a matrix, i.e., $A\ge B$ means that $a_{ij}\ge b_{ij}$, for $\forall i, j$. $x_k$ is the initial state and $\theta_k$ represents the inputs as well as the necessary slack variables depending on different systems and different MPC-based algorithms.

The feasible set of OCP (2) is the largest subset of the state space for which there exists a $\theta_k$ satisfying all the constraints that can be defined as

 $\begin{eqnarray} {F}&\hspace{-0.25cm}=&\hspace{-0.25cm}\{x| \exists \theta_k \text{ such that$J(x, \theta_k)$is feasible}\}\nonumber\\ &\hspace{-0.25cm}=&\hspace{-0.25cm}\{x|H(x, \theta_k)\succeq0, G(x, \theta_k)\ge0\}. \end{eqnarray}$ (3)

Definition 1 (Polytope [18]): A bounded convex polyhedron is a polytope or, equivalently, a set ${P}$ is called a polytope if it can be expressed as the convex hull of finitely many points, i.e., ${P}={\text {conv}}\{x_1, x_2, \ldots, x_p\}$. A $k$-dimensional polytope is called a $k$-polytope. This means that some $(k + 1)$-subfamily of $(x_1, x_2, \ldots, x_p)$ is affinely independent, but no such $(k + 2)$-subfamily is affinely independent. The convex hull of $k +1$ affinely independent points (i.e., $k$-polytope) is a $k$-simplex. A simplex is a $k$-simplex if and only if it has $k + 1$ vertices.

From this definition we can see that a 1-simplex is a closed segment, a 2-simplex is a triangle and a 3-simplex is a tetrahedron. We adopt the convention that the empty set is a polytope of dimension $-1$.

$\mathit{Assumption\ 1:}$ The feasible set of OCP (2) is a $(p-1)$-simplex.

Thus, ${F}$ can be represented as [19]:

 $$$\label{VRepFS} {F}=\{x | x = \lambda_1v_1+\lambda_2v_2+\cdots+\lambda_pv_p\ \text{with}\ \lambda_i\ge0, \sum\limits_{i=1}^{p}\lambda_i=1\}$$$ (4)

where $(v_1, v_2, \ldots, v_p)$ are the $p$ vertices of ${F}$. Note that the description of a feasible set by a set of inequalities (H-representation) in (3) and the convex linear combinations (V-representation) in (4) are equivalent. The goal of this paper is to find all the vertices of ${F}$.

Ⅲ. FEASIBLE SET COMPUTATION

As the value of inner product between two vectors $\vec{a}, \vec{b}$ can be treated as the length of the projection vector from $\vec{b}$ to $\vec{a}$ if $\vec{a}$ is a unit vector, the multi-parametric programming problem (mPPP) for computing the vertices of feasible set then can be expressed as:

 $\begin{eqnarray} &&J(x, c)=\max\limits_{x}: c^{T} x \nonumber\\ &&{\rm s.t.}\, \, H(x, \theta_k)\succeq0, G(x, \theta_k)\ge0 \end{eqnarray}$ (5)

where $\|c\|=1$. In fact, the restriction to the length of vector $\vec{c}$ is not necessary as long as the cost function is to compute the maximum length of the projection from $\vec{x}$ to a vector in some direction. However, if we choose $\vec{c}$ to be unit, the cost function of mPPP (5) would be more meaningful as $c^{T} x =\vec{c}\vec{x}=|\vec{p}|$ where $\vec{p}$ is the projection from $\vec{x}$ to $\vec{c}$.

A. 2-D Case

When $x\in \mathbb{R}^2$, the feasible set can be shown in a plane. For convenience, we assume the $p$ vertices $(v_1, v_2, \ldots, v_p)$ of ${F}$ are ordered. Here, "ordered" means $v_i$ is adjacent to $v_{i+1}$, $v_p$ is adjacent to $v_1$, and the index increases anticlockwise. Let $c=c_\alpha=(\cos\alpha, \sin\alpha)$, where $c$ is the unit vector of the polar coordinator. Let $\bar{{F}}$ denote the boundary of ${F}$ and let $\text{int}({F})$ denote the interior of ${F}$, where ${F}=\bar{{F}}\oplus\text{int}({F})$.

$\mathit{Lemma\ 1:}$ Given a point $x\in {F}$ and a vector $\vec{q}$, if the angle between $\vec{x}-\vec{v}_i$ and $\vec{q}$ is equal or less than $90^\circ$ for $\forall i \in N_{[1, p]}$, then $x$ is on the boundary of polytope ${F}$.

$\mathit{Proof:}$ If $p\le2$, it is obvious that $x\in \bar{{F}}$ as $\bar{{F}}={F}$. When $p\ge3$, assume $x$ is a interior point of ${F}$ and $x_q$ is a point on the boundary of ${F}$ with $\vec{x}-\vec{x}_p=\lambda\vec{q}$ and $\lambda > 0$. Assume $x_q$ is on the edge of $\overline{v_jv_{j+1}}$. Then, there exists a $i$ such that $\angle v_{j}xv_{j+1}+\angle v_{j+1}xv_{j+2}+\cdots+\angle v_{j+i-1}xv_{j+i}\le180^\circ$ and $\angle v_{j}xv_{j+1}+\angle v_{j+1}xv_{j+2}+\cdots+\angle v_{j+i}xv_{j+i+1} > 180^\circ$ where the subscript of $v$ represents the remainder when dividing by $p$ if it is larger than $p$. As $\angle(\vec{x}-\vec{v}_{j+i}, \vec{q}), \angle(\vec{x}-\vec{v}_{j+i+1}, \vec{q})\in [0^\circ, 90^\circ]$, which means that $\angle v_{j+i}xq, \angle v_{j+i+1}xq \in[0^\circ, 90^\circ]$, $\angle v_{j+i}xv_{j+i+1}= 360^\circ- \angle v_{j+i}xq -\angle v_{j+i+1}xq \ge180^\circ$. Since $\angle v_{j+i}xv_{j+i+1}$ is a interior angle of $\bigtriangleup v_{j+i}xv_{j+i+1}$ and less than $180^\circ$ and conflicts with $\angle v_{j+i}xv_{j+i+1}\ge180^\circ$, $x$ cannot be a interior point of ${F}$. Thus, $x\in \bar{{F}}$.

$\mathit{Lemma\ 2:}$ The solution $x^*$ to problem (5) is on the boundary of ${F}$, i.e., $x^*\in\bar{{F}}$.

$\mathit{Proof:}$ As $x^*$ is the solution to problem (5), we have $cx^*\ge cv_i, \forall i \in N_{[1, p]}$. That is $\vec{c}(\vec{x}^*-\vec{v}_i)=|\vec{c} ||\vec{x}^*-\vec{v}_i|\cos(\angle(\vec{c}, \vec{x}^*-\vec{v}_i))\ge0, \forall i \in N_{[1, p]}$. Thus, the angle between $\vec{c}$ and $\vec{x}^*-\vec{v}_i$ is equal to or less than $90^\circ$ for $\forall i \in N_{[1, p]}$. According to Lemma 1, $x^*$ would be on the boundary of ${F}$.

$\mathit{Lemma\ 3:}$ There exists a certain $\alpha_i \in [0^\circ, 360^\circ]$ such that the solution $x^*$ to problem (5) would be the vertice $v_i$ of ${F}$ for $\forall i \in N_{[1, p]}$.

$\mathit{Proof:}$ From Lemma 2 we know that $x^*\in\bar{{F}}$. Assume $\exists x_\alpha \in\bar{{F}}$ such that $x_\alpha\neq v_i$ and $c_\alpha x_\alpha\ge c_\alpha v_i$ for $\forall \alpha \in [0^\circ, 360^\circ]$. That is $\exists x_\alpha(\neq v_i) \in\bar{{F}}$ such that $\vec{c}_\alpha(\vec{x}_\alpha-\vec{v}_i) \ge0$ for $\forall \alpha \in [0^\circ, 360^\circ]$. Let $\alpha=\alpha_1$ such that $\angle(\vec{c}_{\alpha_1}, \vec{v}_{i-1}-\vec{v}_i)=90^\circ$ and $\angle(\vec{c}_{\alpha_1}, \vec{v}_{i+1}-\vec{v}_i) < 90^\circ$. Then, let $\alpha_2=180^\circ +\alpha_1+\epsilon$ where $\epsilon$ is a extremely small angle such that $\angle(\vec{c}_{\alpha_2}, \vec{v}_{i-1}-\vec{v}_i) > 90^\circ$ and $\angle(\vec{c}_{\alpha_2}, \vec{v}_{i+1}-\vec{v}_i) > 90^\circ$. Thus, we cannot find a $x_{\alpha_2} \in\bar{{F}}$ such that $\vec{c}_{\alpha_2}(\vec{x}_{\alpha_2}-\vec{v}_i) \ge0$ as Fig. 1 shows a contradiction.

From Lemma 3 we know that, if we can find all the $\alpha_i$, we will obtain all the vertices of ${F}$. Thus, if we traverse $\alpha$ in $[0^\circ, 360^\circ]$ with a proper interval, the feasible set can be computed. This procedure can be described as

1) Let $k=1$.

2) For $\alpha=0^\circ$ to $360^\circ$ with step length $\varepsilon$ (i.e., $\alpha=\alpha+\varepsilon$ for the next step).

3) Compute $x_\alpha$ by solving $J(x, c_\alpha)$ from OCP (5).

4) Let $x_k=x_\alpha$ and $k=k$+1.

5) Endfor.

6) Compute the convex hull of points set $P=(x_1, x_2, \ldots, x_{\text{end}})$ and let ${S}=\text{conv}\{P\}$.

$\mathit{Theorem\ 1:}$ There exists an appropriate $\varepsilon$ such that ${S}={F}$.

$\mathit{Proof:}$ From Lemma 3 we get that there exists $\alpha_i \in [0^\circ, 360^\circ]$ such that the solution $x^*$ to problem (5) would be the vertice $v_i$ of ${F}$. $\varepsilon$ is chosen such that $\alpha_i=N_i \varepsilon$, where $N_i$ is a proper integer for $\forall i \in N_{[1, p]}$. Thus, $v_i, \forall i \in N_{[1, p]}$ would be in the points set $P$ and ${S}=\text{conv}\{P\}={F}$.

From the proof of theorem 1, we know that the optimal $\varepsilon$ would be the maximum to make $\alpha_i/\varepsilon$ an integer so that the total computation time would be reduced to a minimum.

B. 3-D Case

If $x\in \mathbb{R}^3$, let $c=c_{\alpha, \beta}=(\sin \beta, \cos \beta\sin\alpha, \cos \beta\cos\alpha)$ which is a unit vector in spherical coordinates. For vertice $v_i$, we assume its adjacent $ni$ vertices can be expressed as $(v_{i+1}, \ldots, v_{i+ni})$. A similar lemma for the existence of these vertices can be described as:

$\mathit{Lemma\ 4:}$ There exists a certain $\alpha_i \in [0^\circ, 360^\circ]$ and $\beta_j \in [0^\circ, 360^\circ]$ such that the solution $x^*$ to problem (5) would be vertice $v_i$ of ${F}$ for $\forall i \in N_{[1, p]}$.

$\mathit{Proof:}$ Assume $\exists x_{\alpha, \beta} \in {F}$ such that $x_{\alpha, \beta}\neq v_i$ and $c_{\alpha, \beta} x_{\alpha, \beta}\ge c_{\alpha, \beta} v_i$ for $\forall \alpha \in [0^\circ, 360^\circ]$ and $\forall \beta \in [0^\circ, 360^\circ]$. That is $\exists x_{\alpha, \beta}(\neq v_i) \in {F}$ such that $\vec{c}_{\alpha, \beta}(\vec{x}_{\alpha, \beta}-\vec{v}_i) \ge0$ for $\forall \alpha \in [0^\circ, 360^\circ]$ and $\forall \beta \in [0^\circ, 360^\circ]$. Let $\alpha=\alpha_1, \beta=\beta_1$ such that $\angle(\vec{c}_{\alpha_1, \beta_1}, \vec{v}_{i+1}-\vec{v}_{i})=90^\circ$ and $\angle(\vec{c}_{\alpha_1, \beta_1}, \vec{v}_{i+j}-\vec{v}_i) < 90^\circ, j=2, 3, \ldots, ni$. Then, let $\alpha_2=180^\circ +\alpha_1+\epsilon_1, \beta_2=360^\circ -\beta_1+\epsilon_2$ where $\epsilon_1$ and $\epsilon_2$ are extremely small angles such that $\angle(\vec{c}_{\alpha_2, \beta_2}, \vec{v}_{i+j}-\vec{v}_i) > 90^\circ$ for some $j$ (such $j$ always exists; otherwise, ${F}$ will be concave or ${v}_i$ cannot be a vertice). Thus, we cannot find a $x_{\alpha_2, \beta_2} \in{F}$ such that $\vec{c}_{\alpha_2, \beta_2}(\vec{x}_{\alpha_2, \beta_2}-\vec{v}_i) \ge0$ as Fig. 2 shows this contradiction.

Again, if we traverse $\alpha$ and $\beta$ through $0^\circ$ to $360^\circ$, we will find all the needed vertices. This procedure can be described as:

1) Let $k_1=1, k_2=1$.

2) For $\beta=0^\circ$ to $360^\circ$ with step length $\varepsilon_1$

3)   For $\alpha=0^\circ$ to $360^\circ$ with step length $\varepsilon_2$

4)     Compute $x_{\alpha, \beta}$ by solving $J(x, c_{\alpha, \beta})$ from OCP (5).

5)     Let $x_{k_1, k_2}=x_{\alpha, \beta}$ and $k_2=k_2+1$.

6)   Endfor.

7)   $k_1=k_1+1$.

8) Endfor.

9) Compute the convex hull of points set $P=(x_{1, 1}, x_{1, 2},$ $\ldots, x_{\text{end}, \text{end}})$ and let ${S}=\text{conv}\{P\}$.

$\mathit{Theorem\ 2:}$ There exists an appropriate $\varepsilon=(\varepsilon_1, \varepsilon_2)$ such that ${S}={F}$.

$\mathit{Proof:}$ Similar to Theroem 1, omitted here.

Parameters $\varepsilon_1$ and $\varepsilon_2$ can be chosen individually according to prior knowledge of the structure of the feasible set. If the vertices of the feasible set are smooth in some axis, the according $\varepsilon$ can be chosen a little larger. The best $\varepsilon$ would still be the maximum $(\varepsilon_1, \varepsilon_2)$ to make $\alpha_i/\varepsilon_1$ and $\beta_j/\varepsilon_2$ integers so that the number of circulation can be reduced to minimum.

C. Higher Dimensional Case

When $x\in \mathbb{R}^n$ where $n\ge4$ and considering the n-dimensional spherical coordinates, we specify $c=(c_1, c_2, \ldots, c_n)$ where:

 $\begin{eqnarray} c_1&=&\sin \phi_1\nonumber\\ c_2&=&\cos \phi_1 \sin \phi_2\nonumber\\ c_3&=&\cos\phi_1 \cos \phi_2\sin \phi_3\nonumber\\ & \vdots&\nonumber\\ c_{n-1}&=&\cos \phi_1 \cdots\cos\phi_{n-2}\sin \phi_{n-1}\nonumber\\ c_{n}&=&\cos \phi_1 \cdots\cos\phi_{n-2}\cos\phi_{n-1}. \end{eqnarray}$ (6)

Theoretically, if we traverse $(\phi_1, \phi_2, \ldots, \phi_{n-1})$ from $0^\circ$ to $360^\circ$ with appropriate interval $(\varepsilon_1, \varepsilon_2, \ldots, \varepsilon_{n-1})$, we will end up with a convex hull $\text{conv}\{P\}=\text{conv}\{x_{1, \ldots, 1}, \ldots, x_{\text{end}, \ldots, \text{end}}\}$ with all the vertices $v_i \in \text{conv}\{P\}, \forall i\in N_{[1, p]}$. Thus, ${S}=\text{conv}\{P\}$ just like in the 2-D and 3-D cases. However, if the computation time for every loop is $T_0$ with the majority comes from solving the mPPP (5), the total computation time $T_t=(360^\circ/\varepsilon_0)^{n-1}T_0$ where we treat $\varepsilon_0=\varepsilon_1=\cdots=\varepsilon_{n-1}$ will increase so rapidly as the dimension of $x$ increases that the proposed method may not be tractable.

An alternative is to consider the orthogonal projection if the feasible set is well structured and the projection is computationally tractable. From the definition of orthogonal projection, we know that given ${P}=\{x, \theta|H(x, \theta)\succeq0, G(x, \theta)\ge0\}$, ${F}=Proj_{x}({P})$. However, LMI $H(x, \theta_k)\succeq0$ cannot be handled in projection algorithms in which only V-representation and H-representation are allowed. In fact, the region of LMI ${H}=\{x, \theta_k|H(x, \theta_k)\succeq0\}$ is called a spectrahedron [20]. A lot of work on the connections between polyhedrons and spectrahedrons has been done such as [20], [21]. Recently, [22] proposed a method to convert the spectrahedron to an identical polyhedron. In particular, we have to find a transform matrix $M$ such that $M^{T}H(x, \theta_k)~M={\text{ diag}}~(Q(x, ~\theta_k), ~D(x, ~\theta_k))$ where $D(x, \theta_k))$ is a diagonal matrix map where [22] proved that ${H}=\{x, \theta_k|D(x, \theta_k)\succeq0\}$. Thus, ${H}$ can be expressed in H-representation ${H}=\{x, \theta_k|D_{ii}(x, \theta_k)\ge0\}$. Though [22] proposed an approach to compute $M$ by seeking the joint invariant subspace and its orthonormal basis, this procedure is rather complicated.

$\mathit{Remark\ 1:}$ The proposed method is particularly efficient for low dimensional feasible set computation. When the total computation time $T_t$ becomes unacceptable, the transformed orthogonal projection then should be considered.

Ⅳ. SIMULATION EXAMPLE

$\mathit{Example\ 1}$: Consider the widely used system [10], [23], [24]:

 $$$\label{2D} x_{k+1}= \left( \begin{array}{c c} 1&1\\ 0&1\\ \end{array}\right)x_k+ \left( \begin{array}{c} 0.5\\ 1\\ \end{array}\right)u_k+w_k$$$ (7)

with system constraints $x\in {X}=\{x|-10\le [x]_1\le 2, -10\le [x]_2\le 10\}$, $u\in\mathbb{U}=\{u|-1\le u\le 1\}$ and $w\in{W}=\{w|-0.01\le [w]_1\le 0.01, -0.01\le [w]_2\le 0.01\}$. Assume the applied algorithm is the min-max MPC proposed in [8] with open-loop control law. When the prediction horizon $N=1$, constraints to the mPPP would be

 $\begin{eqnarray} &&Ax+Bu\in {X}\ominus{S}\nonumber\\[2mm] &&x\in {X}\nonumber\\[2mm] &&u\in {U} \end{eqnarray}$ (8)

where $S = \{ x|0.5{[x]_2} \le 0.99, - 0.66{[x]_1} - $$1.33{[x]_2} \le 0.96,0.66{[x]_1} +$$ 1.32{[x]_2} \le 0.96,0.43{[x]_1}$$+ 0.21{[x]_2} \le 0.95, - 0.43{[x]_1} - 0.21{[x]_2} \le 0.95\}$ is the maximum robust control invariant set which can be computed by methods introduced in [4] and $\theta=u$. Thus, the feasible set computed by the proposed method with $\varepsilon=5^\circ$ is shown in Fig. 3 in dark green and set in light green is the set ${P}=\{x, u|Ax+Bu\in {X}\ominus{S}, x\in {X}, u\in {U}\}$. The result is identical to that computed by projection which means $\varepsilon=5^\circ$ works well.

 Download: larger image Fig. 3 Feasible set of Example 1 with $N=1$.

$\mathit{Example\ 2}$: Consider the continuously stirred tank reactor (CSTR) system in [7]:

 $$$\label{LCSTR} x_{k+1}=Ax_k+Bu_k+B_pp_k+B_qq_k$$$ (9)

where

 $\begin{eqnarray*} A&\hspace{-0.25cm}=\hspace{-0.25cm}&\left( \begin{array}{c c c} 9.827\times10^{-1}&-9.803\times10^{-4}& -1.226\times10^{-5} \\ 4.277\times10^{-1} & 1.000\times10^{0} & 2.351\times10^{-2} \\ 4.649\times10^{-2} & 2.043\times10^{-1} & 6.953\times10^{-1} \\ \end{array}\right) \end{eqnarray*}$
 $\begin{eqnarray*} B_p&\hspace{-0.25cm}=\hspace{-0.25cm}&\left( \begin{array}{c c c} 2.828\times10^{0}& 1.713\times10^{-3} &-8.495\times10^{-7}\\ -1.163\times10^{1}& 3.707\times10^{-4}& 1.728\times10^{-3}\\ -1.281\times10^{0}& 2.762\times10^{-5}&1.873\times10^{-4}\\ \end{array}\right) \end{eqnarray*}$
 $\begin{eqnarray*} B&\hspace{-0.25cm}=\hspace{-0.25cm}&\left( \begin{array}{c} 4.961\times10^{-4}\\ -1.468\times10^{0}\\ -9.860\times10^{1} \end{array}\right), B_q=\left( \begin{array}{c} -5.193\times10^{-7}\\ 1.536\times10^{-3}\\ 1.032\times10^{-1}\\ \end{array}\right). \end{eqnarray*}$

Constraints for this CSTR system are $Pr\{-10\leq [x]_1\leq 10\}\ge 90\, \%$, $Pr\{-5\leq[x]_2\leq 5\}\ge 90\, \%$, $Pr\{-2.8\leq [x]_3\leq 2.8\}\ge 90\, \%$, $-2.156\times 10^{-2}\leq u\leq 0.2$. Deterministic disturbances range $-1\leq[p]_2\leq 1$, $-2\leq [p]_3\leq 2$. $[p]_1$ is a persistent deterministic disturbance and $[p]_1=-1\times 10^{-3}$ when $t\geq 30$ mins, otherwise, $[p]_1=0$. Stochastic disturbance $q\sim$ N(0, 10) and $10 \leq q \leq 10$.

Using the proposed algorithms in [7], the final OCP involves both spectrahedral and polyhedral constraints. Choose $\varepsilon_0=\varepsilon_1=\varepsilon_2=3^\circ$ and the feasible set can be computed accurately when $N=12$ as Fig. 4 shows. The dimension of $\theta$ increases when the prediction horizon $N$ increases which in turn increases, the total time to compute the feasible set. We record the total time for the proposed method as well as the orthogonal projection approach based on Fourier method [25] to compute the feasible set. From Fig. 5 we can see that the computation time for the proposed method increases slowly as it depends on the computation time of the basic mPPP once the dimension of the feasible set is decided. However, the computation time as well as memory consumption of orthogonal projection increases rapidly when the dimension of $\theta$ increases, which ultimately leads to the memory exhaustion when $N=10$. Thus, when the dimension of $\theta$ is large and dimension of $x$ is small, the proposed method would be a good choice to compute the feasible set.

 Download: larger image Fig. 4 Feasible set for example 2 with $N$=12.
In this paper, a multi-parametric programming problem has been proposed to compute the feasible set of the final OCP derived from MPC-based algorithms and involving both spectrahedral (represented by LMI) and polyhedral (represented by a set of inequalities) constraints. In fact, the forms or even the convexity and linearity of constraints are not restricted as long as assumption 1 holds. The simulation results showed that the proposed method is especially efficient for low dimensional feasible set computation and avoided the non-unicity problem of optimizers. Because the computation procedure just involves solving the basic mPPP repeatedly, memory consumption is not large. Thus, the proposed algorithm avoids the memory depletion problem encountered by projection algorithms. If we choose a large $\varepsilon$, a glimpse of the feasible set is available as the computed ${S}$ would be a subset of the feasible set. This can be used to find a rough initial state or check the existence of a nonempty feasible set.