Robust Quadratic Optimal Control for Discrete-Time Linear Systems with Non-Stochastic Noises

Abstract

In this paper, the quadratic optimal control problem is investigated for the discrete-time linear systems with process and measurement noises which belong to specified ellipsoidal sets. As the noises are non-stochastic, the traditional Kalman filtering and Dynamic Bellman Equation are not applicable for the proposed control problem. To obtain the optimal control, we firstly converted the multi-step quadratic global optimal control problem to multiple one-step quadratic local approximate optimal control problems. For each one-step quadratic optimal control problem, considering that the system states are not fully available, the set-membership filtering is applied to estimate the true state feasible set. Then based on robust optimization, a robust state feedback control strategy can be obtained by solving a certain semidefinite programming (SDP) problem. The method can not only achieve the optimal control, but also estimate the system states more accurately. Finally, the simulation results verify the effectiveness of the proposed algorithm.

1. Introduction

Linear quadratic (LQ) optimal control problems with indeterminacy factors have been widely studied in theoretical and practical perspectives. In most cases, these indeterminacy factors are modeled as stochastic processes. Many achievements in the study of stochastic control have been made, such as minimum variance control [1], Linear-quadratic-Gaussian (LQG) [2,3,4] and probability density function (PDF) control [5]. Actually, due to lack of sufficient information and data, the statistical properties of the uncertainties in some special systems are very difficult to obtain. Taking the mobile system as an example, because the terrain conditions cannot be reused, the stochastic model can not be modeled with the insufficient data [6]. Moreover, there are several uncertainties which are not essentially stochastic.

Non-stochastic uncertainties as an alternative which only consider the range of uncertain information have attracted the interests of scholars. Terra proposed a recursive robust linear quadratic regulator based on the robust weighted least-squares method for discrete-time systems with norm-bounded parametric uncertainties [7]. Huang investigated the robust quadratic optimal control based on robust optimization for the systems in the presence of ellipsoidal bounded disturbance [8]. Zhang studied the robust optimal control with adjustable uncertainty sets which treated the uncertainty sets as additional decision variables [9]. However, the methods above were with the assumption that the system states were observed completely.

Actually, the states in the real systems are usually not measurable completely and accurately. Set-membership filtering can provide an estimated feasible set of true states for systems with non-stochastic noises, and so far, many approaches have been proposed to obtain the solutions [10,11,12,13]. Based on the set-membership estimation theory, Qian studied the output tracking control problem for systems with ellipsoidal uncertainties in system parameters and noise, and an adaptive robust tracking control algorithm was proposed [14]. However, the performance index was given with $l_1$ norm of the output deviation over the finite time horizon, whose transient performance is poorer than the linear quadratic optimal control. Bertsimas studied the linear quadratic control problems with norm-bounded uncertainties [15]. However, for the imperfect state information case, the state estimation error set was assumed in fixed bounded set which neglected the learning feature of the state estimation. Ping investigated the output feedback robust predictive control with ellipsoidal estimation error sets [16,17]. However, the methods updated the estimation state error set by solving one convex optimal problem each time.

Motivated by the above discussions, the LQ control problem with non-stochastic noises is considered in this paper. Specifically, the non-stochastic process and measurement noises are confined to ellipsoidal bounded sets. Ellipsoidal bounded sets as more realistic and flexible non-stochastic uncertain sets, on the one hand, are of good tractability and well-motivated from a stochastic setting, where they naturally occur when a normal distribution is cut off at a certain level of probability. On the one hand, they can describe kinds of bounded uncertainties as well, such as 2-norm bounded (by selecting the identity matrix as the size matrix for the ellipsoid set) [18]. In addition, compared with the traditional stochastic assumption, the bound of the noises can be obtained more easily. As the noises are bounded, the traditional Kalman filtering and Dynamic Bellman Equation are not applicable the proposed control problem. For the LQ control methods with ellipsoidal bounded noises, most of them considered the states to be completely observed and for the output feedback control methods, most of them need to solve one complex optimization problem at each time to update the estimation state error set. In this paper, we converted the multi-step quadratic global optimal control problem to multiple one-step quadratic local approximate optimal control problems. Then set-membership filtering and robust optimization are adopted to design the state feedback control strategy which ensures the optimality of control performance (the quadratic cost function about the system states and the control law).

The remaining part of this paper is organized as follows: In Section 2, the problem formulation is presented. Then, in Section 3, we present the robust control algorithm based on the set-membership filtering and robust optimization. To show the effectiveness of our proposed method, simulation results are presented in Section 4. Finally, conclusions and future research directions are presented in Section 5.

Notations. Throughout the paper, $S\ge 0(S>0)$ means that the real symmetric matrix S is non-negative definite (positive definite). I is the identity matrix with appropriate dimensions. The ellipsoid $\Omega (s,S)$ is denoted by

$$\Omega (s,S)\triangleq \{l\in R^n:{(l-s)}^T{S}^{-1}(l-s)\le 1\},$$

where the vector s is the center of the ellipsoid, the matrix $S(\ge 0)$ characterizes its shape.

2. Problem Formulation

Consider the following uncertain discrete-time linear dynamic system,

$$\begin{array}{ccc}\hfill x(k+1)& =& A\left(k\right)x\left(k\right)+B\left(k\right)u\left(k\right)+w\left(k\right),k=0,1,\dots ,N-1,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill y\left(k\right)& =& C\left(k\right)x\left(k\right)+v\left(k\right),k=1,2,\dots ,N,\hfill \end{array}$$

(2)

where $x\left(k\right)\in R^n$ is the state vector, $u\left(k\right)\in R^m$ is the control vector, $y\left(k\right)\in R^p$ is the measure output vector. $w\left(k\right)\in R^n$ and $v\left(k\right)\in R^p$ are the process noise and observation noise respectively, which are confined to specified ellipsoidal sets,

$$w\left(k\right)\in \Omega (0,W(k\left)\right),$$

$$v\left(k\right)\in \Omega (0,V(k\left)\right),$$

where $W\left(k\right)$ and $V\left(k\right)$ represent the bounds of noises $w\left(k\right)$ and $v\left(k\right)$ respectively. The initial state $x\left(0\right)$ takes value in the following ellipsoidal set,

$$x\left(0\right)\in \Omega (\widehat{x}\left(0\right),P\left(0\right)).$$

Considering the quadratic cost function,

$$\begin{array}{c}\hfill J_o=\sum _{k=0}^N{x}^T\left(k\right)Q\left(k\right)x\left(k\right)+\sum _{k=0}^{N-1}{u}^T\left(k\right)R\left(k\right)u\left(k\right),\end{array}$$

(3)

where matrices $Q\left(k\right)\ge 0$ and $R\left(k\right)>0$ are non-negative definite and positive definite symmetric matrices of appropriate dimensions, respectively.

Our target of this work is to find a feedback control series ${\left\{u\left(k\right)\right\}}_{k=0}^{N-1}$ such that the quadratic cost function $J_o$ is minimized when the noises $w\left(k\right)$ and $v\left(k\right)$ belong to bounded ellipsoidal uncertain sets. Namely,

$$\begin{array}{ccc}\hfill & & \underset{u\left(0\right),\cdots ,u(N-1)}{min}{J}_o\hfill \\ & & s.t.x(k+1)=A\left(k\right)x\left(k\right)+B\left(k\right)u\left(k\right)+w\left(k\right)\hfill \\ & & y\left(k\right)=C\left(k\right)x\left(k\right)+v\left(k\right)\hfill \\ & & w\left(k\right)\in \Omega (0,W(k\left)\right)\hfill \\ & & v\left(k\right)\in \Omega (0,V(k\left)\right).\hfill \end{array}$$

(4)

Due to the noises $w\left(k\right)$ and $v\left(k\right)$ being bounded in this paper, the dynamics Bellman equation cannot be applied to solve the multi-step quadratic optimal control problem. In order to effectively address the proposed optimal control problem, we attempt to find a feasible approximate optimal solution.

Because the value of $x^T\left(0\right)Q\left(0\right)x\left(0\right)$ is not influenced by the control law $u\left(k\right)$, so we simplified the cost function to $J={\sum }_{k=0}^{N-1}[x(k+1)Q(k+1)x(k+1)+u^T\left(k\right)R\left(k\right)u\left(k\right)]$. To ensure the optimality of J for all the possible values in the noises sets $\Omega (0,W(k\left)\right)$ and $\Omega (0,V(k\left)\right)$, the minimization of the worst case of J is considered. To realize the approximate optimality of the multi-step quadratic cost function J, we optimize each one-step quadratic cost function $J\left(k\right)$($k=0,1,\dots ,N-1$) which is defined as follows,

$$\begin{array}{c}\hfill J\left(k\right)=x^T(k+1)Q(k+1)x(k+1)+u^T\left(k\right)R\left(k\right)u\left(k\right).\end{array}$$

(5)

Then, the multi-step quadratic global optimal control problem is converted to multiple one-step quadratic local optimal control problems which can provide a feasible approximate optimal solution. Each one-step quadratic optimal control problem is as follows,

$$\begin{array}{ccc}\hfill & & \underset{u\left(k\right)}{min}\underset{w\left(k\right)\in \Omega (0,W(k\left)\right),v\left(k\right)\in \Omega (0,V(k\left)\right)}{max}J\left(k\right)\hfill \\ & & s.t.x(k+1)=A\left(k\right)x\left(k\right)+B\left(k\right)u\left(k\right)+w\left(k\right)\hfill \\ & & y\left(k\right)=C\left(k\right)x\left(k\right)+v\left(k\right).\hfill \end{array}$$

(6)

In optimization problem (6), the performance index $J\left(k\right)$ is a function with respect to the system states which are not measurable completely. As a result, the state estimation error directly affects the feedback control performance. If we can yield more accurate true state feasible set, the maximum value of $J\left(k\right)$ can be more accurate. Then by minimizing the maximum value of $J\left(k\right)$, optimal control law $u\left(k\right)$ can be obtained as well.

The following lemmas are used in the designing of the optimal control law.

Lemma 1.

S-procedure [19] Let ${\tilde{F}}_0,\cdots ,{\tilde{F}}_p$ be quadratic functions of the variable $\zeta \in R^m$,

$${\tilde{F}}_i\left(\zeta \right)\triangleq {\zeta }^T{T}_i\zeta +2{\varphi }_i^T\zeta +{\phi }_i,i=0,\cdots ,p,$$

where $T_i=T_i^T$, ${\varphi }_i$ and ${\phi }_i$ are with appropriate dimensions. The following condition on ${\tilde{F}}_0$, ⋯, ${\tilde{F}}_p$,

$${\tilde{F}}_0\left(\zeta \right)\ge 0,$$

for all ζ such that

$${\tilde{F}}_i\left(\zeta \right)\ge 0,i=1,\cdots ,p,$$

holds if there exist

$${\tau }_i\ge 0,$$

such that

$$\left[\begin{array}{cc}{T}_0& {\varphi }_0\\ {\varphi }_0^T& {\phi }_0\end{array}\right]-\sum _{i=1}^p{\tau }_i\left[\begin{array}{cc}{T}_i& {\varphi }_i\\ {\varphi }_i^T& {\phi }_i\end{array}\right]\ge 0.$$

When $p=1$, the converse holds, provided that there is some ${\zeta }_0$ such that ${\tilde{F}}_1\left({\zeta }_0\right)\ge 0$.

Lemma 2.

Schur complement [15] Given the matrices $L=L^T$, $Z=Z^T$ and U with appropriate dimensions. Let

$$H=\left[\begin{array}{cc}L& U^T\\ U& Z\end{array}\right],$$

where $L>0$. Then H is positive (semi) definite if and only if the matrix $Z-U{L}^{-1}{U}^T$ is positive (semi) definite.

3. Robust Optimal Control Strategy

Considering that the states are not measured fully, we firstly adopted set-membership filtering to realize the states estimation. Then, based on the states estimation, taking the negative impact of the estimation error and noises into consideration, we designed the robust optimal control strategy in this paper.

3.1. Estimation of True State Feasible Set

For the unknown-but-bounded noises, set-membership filtering can provide a true state feasible set which contains the true states all the time. In addition, the filtering gain can be obtained off-line with small computation burden.

According to set-membership filtering theory, assuming that at instant $k-1$ the true state feasible set is $\Omega \left(\widehat{x}(k-1),P(k-1)\right)$, where $\widehat{x}(k-1)$ is the center of state feasible set, $P(k-1)$ is used to define its shape.

At instant k, the true states $x\left(k\right)$ belong to the optimal ellipsoidal feasible set $\Omega \left(\widehat{x}\left(k\right),P\left(k\right)\right)$, the center $\widehat{x}\left(k\right)$ and the matrix $P\left(k\right)$ can be updated as follows,

$$\begin{array}{ccc}\hfill \widehat{x}\left(k\right)& =& \widehat{x}\left(k\right|k-1)+K\left(k\right)\delta \left(k\right),\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill \widehat{x}\left(k\right|k-1)& =& A(k-1)\widehat{x}(k-1)+B(k-1)u(k-1),\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill K\left(k\right)& =& P\left(k\right|k-1)C^T\left(k\right){[C\left(k\right)P\left(k\right|k-1)C^T\left(k\right)+q_k^{-1}V\left(k\right)]}^{-1},\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill P\left(k\right|k-1)& =& (p_k^{-1}+1)A(k-1)P(k-1)A^T(k-1)+(p_k+1)W(k-1),\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill P\left(k\right)& =& {\beta }_k\{\left[I-K\left(k\right)C\left(k\right)\right]P\left(k\right|k-1){\left[I-K\left(k\right)C\left(k\right)\right]}^T+q_k^{-1}K\left(k\right)V\left(k\right)K^T\left(k\right)\},\hfill \end{array}$$

(11)

$$\begin{array}{ccc}\hfill {\beta }_k& =& 1+q_k-{\delta }^T\left(k\right){[q_k^{-1}V\left(k\right)+C\left(k\right)P\left(k\right|k-1)C^T\left(k\right)]}^{-1}\delta \left(k\right),\hfill \end{array}$$

(12)

$$\begin{array}{ccc}\hfill \delta \left(k\right)& =& y\left(k\right)-C\left(k\right)\widehat{x}\left(k\right|k-1).\hfill \end{array}$$

(13)

The scalar parameters $p_k>0$ and $q_k>0$ are selected to make the sizes of ellipsoidal sets $\Omega \left(\widehat{x}\left(k\right|k-1),P\left(k\right|k-1)\right)$ and $\Omega \left(\widehat{x}\left(k\right),P\left(k\right)\right)$ to be minimized respectively. There are two common criteria to optimize size of the ellipsoid. One is Minimum-trace criterion, and the other is Minimum-volume criterion [13]. In this paper, we adopt the latter one and $p_k$ satisfies

$$\sum _{{}_{i=1}}^n\frac{1}{{\gamma }_i^{\left(p\right)}+p_k}=\frac{n}{p_k(p_k+1)},$$

where ${\gamma }_i^{\left(p\right)}$ denote the eigenvalues of matrix $A(k-1)P(k-1)A^T(k-1)W^{-1}\left(k\right)$ and n is the dimension of the state. $q_k$ satisfies

$$\sum _{{}_{i=1}}^n\frac{{\gamma }_i^{\left(q\right)}}{1+q_k{\gamma }_i^{\left(q\right)}}=n\frac{{\beta }_k^{{}^{\prime }}\left(q_k\right)}{{\beta }_k\left(q\left(k\right)\right)},$$

where ${\gamma }_i^{\left(q\right)}$ denote the eigenvalues of matrix $P\left(k\right|k-1)C^T\left(k\right)V^{-1}\left(k\right)C\left(k\right)$, ${\beta }_k\left(q_k\right)$ is given by (12) and ${\beta }_k^{{}^{\prime }}\left(q_k\right)$ is the derivative of ${\beta }_k\left(q_k\right)$ with respect to $q_k$.

3.2. Analysis of State Estimation Error

The center $\widehat{x}\left(k\right)$ of true state feasible ellipsoid set is adopted as the estimated state in this paper, then state estimation error can be defined as follows,

$$\eta \left(k\right)=x\left(k\right)-\widehat{x}\left(k\right).$$

Based on the set-membership filtering results, it yields that at instant k the estimation error $\eta \left(k\right)$ belongs to the following ellipsoid set,

$$\eta \left(k\right)\in \Omega (0,P(k\left)\right),$$

namely,

$$\eta \left(k\right)\in \left\{\eta \left(k\right)\right|{\eta }^T\left(k\right)P^{-1}\left(k\right)\eta \left(k\right)\le 1\},$$

the shape matrix $P\left(k\right)$ of error ellipsoid set can be refreshed according to Equations (9)–(13), which avoids solving complex convex optimization problem.

3.3. Robust Control Strategy

According to the definition $\eta \left(k\right)$ and the dynamics system state Equation (1), we have

$$\begin{array}{ccc}\hfill x(k+1)& =& A\left(k\right)x\left(k\right)+B\left(k\right)u\left(k\right)+w\left(k\right)\hfill \\ & =& A\left(k\right)\widehat{x}\left(k\right)+A\left(k\right)\eta \left(k\right)+B\left(k\right)u\left(k\right)+w\left(k\right).\hfill \end{array}$$

(14)

Here, let

$$G\left(k\right)=\left[\begin{array}{cc}I& A\left(k\right)\end{array}\right],$$

$$\rho \left(k\right)=\left[\begin{array}{c}w\left(k\right)\\ \eta \left(k\right)\end{array}\right],$$

then Equation (14) can be rewritten as follows,

$$\begin{array}{c}\hfill x(k+1)=A\left(k\right)\widehat{x}\left(k\right)+B\left(k\right)u\left(k\right)+G\left(k\right)\rho \left(k\right).\end{array}$$

(15)

Substituting Equation (15) into Equation (5), it has

$$\begin{array}{ccc}\hfill J\left(k\right)& =& {\widehat{x}}^T\left(k\right)\tilde{A}\widehat{x}\left(k\right)+2b^{T}u\left(k\right)+u^T\left(k\right)\tilde{B}u\left(k\right)\hfill \\ & & +2c^T\rho \left(k\right)+2u^T\left(k\right)D\rho \left(k\right)+{\rho }^T\left(k\right)\tilde{G}\rho \left(k\right),\hfill \end{array}$$

(16)

where

$$\tilde{A}=A^T\left(k\right)Q(k+1)A\left(k\right),$$

$$b=B^T\left(k\right)Q^T(k+1)A\left(k\right)\widehat{x}\left(k\right),$$

$$\tilde{B}=B^T\left(k\right)Q(k+1)B\left(k\right)+R\left(k\right),$$

$$c=G^T\left(k\right)Q^T(k+1)A\left(k\right)\widehat{x}\left(k\right),$$

$$D=B^T\left(k\right)Q(k+1)G\left(k\right),$$

$$\tilde{G}=G^T\left(k\right)Q(k+1)G\left(k\right).$$

It can be seen from (16) that $J\left(k\right)$ contains two uncertainties, one is the process noise and the other is the estimation state error. To realize the control optimality, how to deal with the uncertainties effectively is important. In this subsection, we introduce two uncertainties into the constrains of optimization problem, based on the lemmas outlined in Section 2, the uncertain problem (6) can be converted into one certain semidefinite programming (SDP) problem by using the robust optimization technique.

The main result is given by the following theorem.

Theorem 1.

Under the bounded ellipsoid uncertain noises, the problem (6) can be solved by the following certain semi-definite programming (SDP) problem,

$$\begin{array}{ccc}\hfill & & min\lambda \hfill \\ & & s.t.\hfill \\ & & \left[\begin{array}{ccc}I& z& F\\ z^T& \lambda -\sum _{i=1}^2{\tau }_i& -h^T\\ F^T& -h& -\tilde{G}+\sum _{i=1}^2{\tau }_i{M}_i\left(k\right)+F^{T}F\end{array}\right]\ge 0\hfill \\ & & {\tau }_i\ge 0,\hfill \end{array}$$

(17)

where λ, ${\tau }_1$, ${\tau }_2$ and z are the decision variables, and

$$F={\tilde{B}}^{-1/2}D,h=c-D^T{\tilde{B}}^{-1}b,$$

$$M_1\left(k\right)=\left[\begin{array}{cc}{W}^{-1}\left(k\right)& 0\\ 0& 0\end{array}\right],$$

$$M_2\left(k\right)=\left[\begin{array}{cc}0& 0\\ 0& P^{-1}\left(k\right)\end{array}\right],$$

the control law is

$$u\left(k\right)=-{\tilde{B}}^{-1}b+{\tilde{B}}^{-1/2}z.$$

Proof of Theorem 1. 

For the system performance index in Equation (16), because $Q(k+1)\ge 0$ and $R\left(k\right)>0$, therefore $B^T\left(k\right)Q(k+1)B\left(k\right)\ge 0$ and $\tilde{B}>0$, respectively. Meanwhile, because $Q(k+1)$ and $R\left(k\right)$ are symmetric matrices, so $\tilde{B}$ is symmetric and positive definite matrix. Therefore, it can yield that ${\tilde{B}}^{-1/2}={\left({\tilde{B}}^{-1/2}\right)}^T$.

Let $u\left(k\right)={\tilde{B}}^{-1/2}z-{\tilde{B}}^{-1}b$, then substituting it into Equation (16), it has

$$\begin{array}{ccc}\hfill J\left(k\right)& =& {\widehat{x}}^T\left(k\right)\tilde{A}\widehat{x}\left(k\right)-b^T{\tilde{B}}^{-1}b+{\rho }^T\left(k\right)\tilde{G}\rho \left(k\right)\hfill \\ & & +z^{T}z+2{(c-D^T{\tilde{B}}^{-1}b)}^T\rho \left(k\right)\hfill \\ & & +2z^T\left({\tilde{B}}^{-1/2}D\right)\rho \left(k\right).\hfill \end{array}$$

At stage k, the first two parts of $J\left(k\right)$ can be considered as a known constant. Let

$$F={\tilde{B}}^{-1/2}D,h=c-D^T{\tilde{B}}^{-1}b,$$

the optimal problem of $J\left(k\right)$ can be converted to be that of $\tilde{J}\left(k\right)$,

$$\tilde{J}\left(k\right)=z^{T}z+2h^T\rho \left(k\right)+2z^{T}F\rho \left(k\right)+{\rho }^T\left(k\right)\tilde{G}\rho \left(k\right).$$

Therefore, the problem (6) is equivalent to the following problem,

$$\begin{array}{c}\hfill \underset{u\left(k\right)}{min}\underset{w\left(k\right)\in \Omega (0,W(k\left)\right),\eta \left(k\right)\in \Omega (0,P(k\left)\right)}{max}\tilde{J}\left(k\right).\end{array}$$

(18)

By introducing an auxiliary variable $\lambda$, problem (18) can be written as follows,

$$\begin{array}{ccc}\hfill min& & \lambda \hfill \\ \hfill s.t.& & \tilde{J}\left(k\right)\le \lambda \hfill \\ & & w^T\left(k\right)W^{-1}\left(k\right)w\left(k\right)\le 1\hfill \\ & & {\eta }^T\left(k\right)P^{-1}\left(k\right)\eta \left(k\right)\le 1,\hfill \end{array}$$

(19)

the constrain in (19)

$$w^T\left(k\right)W^{-1}\left(k\right)w\left(k\right)\le 1,$$

can be rewritten in the following form

$${\left[\begin{array}{c}w\left(k\right)\\ \eta \left(k\right)\end{array}\right]}^T\left[\begin{array}{cc}{W}^{-1}\left(k\right)& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}w\left(k\right)\\ \eta \left(k\right)\end{array}\right]\le 1.$$

Let

$$M_1\left(k\right)=\left[\begin{array}{cc}{W}^{-1}\left(k\right)& 0\\ 0& 0\end{array}\right],$$

it yields

$${\rho }^T\left(k\right)M_1\left(k\right)\rho \left(k\right)\le 1,$$

namely,

$$\begin{array}{c}{\left[\begin{array}{c}1\\ \rho \left(k\right)\end{array}\right]}^T\left[\begin{array}{cc}1& 0\\ 0& -M_1\left(k\right)\end{array}\right]\left[\begin{array}{c}1\\ \rho \left(k\right)\end{array}\right]\ge 0,\end{array}$$

(20)

the constrain in (19)

$${\eta }^T\left(k\right)P^{-1}\left(k\right)\eta \left(k\right)\le 1,$$

can be rewritten in the following form

$${\left[\begin{array}{c}w\left(k\right)\\ \eta \left(k\right)\end{array}\right]}^T\left[\begin{array}{cc}0& 0\\ 0& P^{-1}\left(k\right)\end{array}\right]\left[\begin{array}{c}w\left(k\right)\\ \eta \left(k\right)\end{array}\right]\le 1.$$

Let

$$M_2\left(k\right)=\left[\begin{array}{cc}0& 0\\ 0& P^{-1}\left(k\right)\end{array}\right],$$

it yields

$${\rho }^T\left(k\right)M_2\left(k\right)\rho \left(k\right)\le 1,$$

namely,

$$\begin{array}{c}{\left[\begin{array}{c}1\\ \rho \left(k\right)\end{array}\right]}^T\left[\begin{array}{cc}1& 0\\ 0& -M_2\left(k\right)\end{array}\right]\left[\begin{array}{c}1\\ \rho \left(k\right)\end{array}\right]\ge 0.\end{array}$$

(21)

And the constrain in (19)

$$\tilde{J}\left(k\right)\le \lambda ,$$

namely,

$$\lambda -z^{T}z-2h^T\rho \left(k\right)-2z^{T}F\rho \left(k\right)-{\rho }^T\left(k\right)\tilde{G}\rho \left(k\right)\ge 0,$$

can be rewritten as follows

$$\begin{array}{c}{\left[\begin{array}{c}1\\ \rho \left(k\right)\end{array}\right]}^T\left[\begin{array}{cc}\lambda -z^{T}z& -h^T-z^{T}F\\ -h-F^{T}z& -\tilde{G}\end{array}\right]\left[\begin{array}{c}1\\ \rho \left(k\right)\end{array}\right]\ge 0.\end{array}$$

(22)

According to Lemma 1 (S-procedure), for all $\rho \left(k\right)$ that satisfies (20) and (21), the Equation (22) holds if there exist

$${\tau }_i\ge 0,i=1,2,$$

such that

$$\begin{array}{c}\left[\begin{array}{cc}\lambda -z^{T}z& -h^T-z^{T}F\\ -h-F^{T}z& -\tilde{G}\end{array}\right]-\sum _{i=1}^2{\tau }_i\left[\begin{array}{cc}1& 0\\ 0& -M_i\left(k\right)\end{array}\right]\ge 0,\end{array}$$

namely,

$$\begin{array}{c}\hfill \left[\begin{array}{cc}\lambda -\sum _{i=1}^2{\tau }_i& -h^T\\ -h& -\tilde{G}+\sum _{i=1}^2{\tau }_i{M}_i\left(k\right)+F^{T}F\end{array}\right]-{\left[\begin{array}{cc}z& F\end{array}\right]}^T\left[\begin{array}{cc}z& F\end{array}\right]\ge 0.\end{array}$$

According to Lemma 2 (Schur complement), it yields

$$\begin{array}{c}\hfill \left[\begin{array}{ccc}I& z& F\\ z^T& \lambda -\sum _{i=1}^2{\tau }_i& -h^T\\ F^T& -h& -\tilde{G}+\sum _{i=1}^2{\tau }_i{M}_i\left(k\right)+F^{T}F\end{array}\right]\ge 0.\end{array}$$

This completes the proof. □

The proposed control law $u\left(k\right)=-{\tilde{B}}^{-1}b+{\tilde{B}}^{-1/2}z$, let $\Gamma \left(k\right)={\tilde{B}}^{-1}\left[B^T\left(k\right)Q^T(k+1)A\left(k\right)\right]$, then the control law can be rewiritten as $u\left(k\right)=-\Gamma \left(k\right)\widehat{x}\left(k\right)+{\tilde{B}}^{-1/2}z$, which contains two parts: the estimated states feedback and uncertainties compensation. For the first part, the feedback gain $\Gamma \left(k\right)={[B^T\left(k\right)Q(k+1)B\left(k\right)+R\left(k\right)]}^{-1}\left[B^T\left(k\right)Q^T(k+1)A\left(k\right)\right]$ can be computed off-line. For the second part ${\tilde{B}}^{-1/2}z$, it denotes the compensation for the estimated states error and noises which can be obtained by solving the semi-definite programming problem (17).

4. Simulation Analysis

This section illustrates the features of the non-stochastic robust control law proposed in this paper by the following two examples compared with the output feedback robust model predictive control(OFRMPC) in [16].

The proposed robust optimal control algorithm can be run as follows.

Step 1. Set $k=0$, input $\widehat{x}\left(0\right)$ and $P\left(0\right)$.

Step 2. Update true state feasible set $\Omega (\widehat{x}(k+1),P(k+1))$ according to Equations (7)–(13).

Step 3. Calculate control law $u\left(k\right)$ by solving optimal problem (17).

Step 4. Apply $u\left(k\right)$ to the system. If $k=N-1$, stop; Otherwise, set $k=k+1$ and go back to Step 2.

Example 1.

Consider the following two-dimension uncertain linear system,

$$\begin{array}{ccc}\hfill x(k+1)& =& \left[\begin{array}{cc}0.60& 0.70\\ 0.25& 0.59\end{array}\right]x\left(k\right)+\left[\begin{array}{c}3\\ 0.3\end{array}\right]u\left(k\right)+w\left(k\right),k=0,1,\cdots ,29,\hfill \\ \hfill y\left(k\right)& =& \left[\begin{array}{cc}0& 1\end{array}\right]x\left(k\right)+v\left(k\right),k=1,\cdots ,30,\hfill \end{array}$$

where $x\left(k\right)\in R^2$, $u\left(k\right)\in R$, $y\left(k\right)\in R$, $w\left(k\right)\in R^2$, and $v\left(k\right)\in R$ are the state, the control input, the measured output, the process noise and the observation noise, respectively.

The noises are independent mutually and belong to the ellipsoid sets,

$$w\left(k\right)\in \Omega (0,W(k\left)\right),$$

$$v\left(k\right)\in \Omega (0,V(k\left)\right),$$

where

$$W\left(k\right)=\left[\begin{array}{cc}{\left(0.1arctan\left(k\right)\right)}^2& 0\\ 0& {\left(0.1arctan\left(k\right)\right)}^2\end{array}\right],$$

$$V\left(k\right)={\left(0.1arctan\left(k\right)\right)}^2.$$

The original state $x\left(0\right)$ belongs to the following ellipsoid set,

$$x\left(0\right)\in \Omega (\widehat{x}\left(0\right),P\left(0\right)),$$

where

$$\widehat{x}\left(0\right)=\left[\begin{array}{c}1\\ 1\end{array}\right],$$

$$P\left(0\right)=\left[\begin{array}{cc}25& 0\\ 0& 25\end{array}\right].$$

The matrices in performance index are with

$$Q\left(k\right)=\left[\begin{array}{cc}0.8& 0\\ 0& 0.8\end{array}\right],$$

$$R\left(k\right)=1,$$

and they are selected as non-negative definite and positive definite symmetric matrices of appropriate dimensions, respectively. By using Sedumi toolbox, LMI toolbox and Optimization Toolbox in MATLAB, the simulation results are as follows.

Figure 1 and Figure 2 show the comparison of the system state and the estimated states between the proposed method and the OFRMPC method. Figure 3 and Figure 4 display the estimated errors of system states for the proposed method and OFRMPC method. Here, the estimated states are the centers of the estimated state feasible sets. It can be seen that the estimated states with the proposed method are closer to the true state.

Figure 5 displays the evolution of true state feasible sets for the proposed method. The blue line means the bound of estimated state feasible ellipsoidal sets. It shows that the feasible states sets from instant 2 become smaller and smaller. Meanwhile, they contain the true states all the time and this can be seen directly from the partial enlarged detail in the figure.

Figure 6 shows the curve of the control law for the proposed method and the OFRMPC method. It is observed that the control law proposed in this paper tends to zero and the control law for OFRMPC fluctuates around zero. The upper bound $\lambda$ of the cost function at each time is shown in Figure 7 and they become small gradually. And the quadratic cost function of the proposed method is 49.753, while the quadratic cost function for OFRMPC method is 55.324. This implies that the proposed method has better control performance.

Example 2.

For the position control problem of the industrial robot with six degrees of freedom, because most of the motion in the robot processing is at low speed except the singularity, the coupling of the motion of the robot is weak. In addition, a reasonable path can be planned for the singularity to keep the robot away from its own singularity and avoid strong coupling between movements in different directions. In reference [20], the position control model of industrial robot is simplified as an independent but identical second-order system with disturbance in each direction,

$$\left[\begin{array}{c}\dot{X}\\ \dot{Y}\\ \dot{Z}\end{array}\right]=\left[\begin{array}{ccc}{A}_X& 0& 0\\ 0& A_Y& 0\\ 0& 0& A_Z\end{array}\right]\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]+\left[\begin{array}{ccc}{B}_X& 0& 0\\ 0& B_Y& 0\\ 0& 0& B_Z\end{array}\right]\left[\begin{array}{c}{u}_X\\ u_Y\\ u_Z\end{array}\right]+\left[\begin{array}{c}{d}_X\\ d_Y\\ d_Z\end{array}\right],$$

where $X={(X1,X2)}^T$, $Y={(Y1,Y2)}^T$, $Z={(Z1,Z2)}^T$ are the displacement and velocity for the center point of the end of the manipulator along X, Y, Z direction respectively. The system matrices and input matrices in each direction are as follows,

$$A_X=\left[\begin{array}{cc}0& 1\\ -1172.5·t^{-2}& -46.2·t^{-1}\end{array}\right],B_X=\left[\begin{array}{c}0\\ 1172.5·t^{-2}\end{array}\right],$$

$$A_Y=\left[\begin{array}{cc}0& 1\\ -536.8·t^{-2}& -38.9·t^{-1}\end{array}\right],B_Y=\left[\begin{array}{c}0\\ 536.8·t^{-2}\end{array}\right],$$

$$A_Z=\left[\begin{array}{cc}0& 1\\ -2285.8·t^{-2}& -132.2·t^{-1}\end{array}\right],B_Z=\left[\begin{array}{c}0\\ 2258.8·t^{-2}\end{array}\right],$$

where t is time. The disturbance matrices in each direction are as follows,

$$d_X=\left[\begin{array}{c}{d}_{X1}\\ d_{X2}\end{array}\right],d_Y=\left[\begin{array}{c}{d}_{Y1}\\ d_{Y2}\end{array}\right],d_Z=\left[\begin{array}{c}{d}_{Z1}\\ d_{Z2}\end{array}\right],$$

which contain the system noises, the un-modeled dynamics and the coupling effect of the motion of the robot in X, Y, Z directions and are bounded.

Because the dynamics system in each direction has the same structure, here we take the motion control in X direction as an example to show the effectiveness of the proposed method. The laser tracker is used to measure the position for the center point of the end of the manipulator, the measurement accuracy is $15\mathsf{\mu }\mathrm{m}+5\mathsf{\mu }$m/m. By discretizing the dynamics system in X direction with the discrete sampling time $T=0.01s$, it yields

$$\begin{array}{ccc}\hfill X(k+1)& =& \left[\begin{array}{cc}1& 0.01\\ -11.725·k^{-2}& 1-0.462·k^{-1}\end{array}\right]X\left(k\right)+\left[\begin{array}{c}0\\ 11.725·k^{-2}\end{array}\right]u_X\left(k\right)+d_X\left(k\right),\hfill \\ \hfill y\left(k\right)& =& \left[\begin{array}{cc}1& 0\end{array}\right]X\left(k\right)+v\left(k\right).\hfill \end{array}$$

The noises are independent mutually and belong to the ellipsoid sets,

$$d_X\left(k\right)\in \Omega (0,W\left(k\right)),$$

$$v\left(k\right)\in \Omega (0,V(k\left)\right),$$

where

$$W\left(k\right)=\left[\begin{array}{cc}9& 0\\ 0& 25\end{array}\right],$$

$$V\left(k\right)=4\times {10}^{-4}.$$

The original state $X\left(0\right)$ belongs to the following ellipsoid set,

$$X\left(0\right)\in \Omega (\widehat{X}\left(0\right),P\left(0\right)),$$

where

$$\widehat{X}\left(0\right)=\left[\begin{array}{c}2\\ 10\end{array}\right],$$

$$P\left(0\right)=\left[\begin{array}{cc}100& 0\\ 0& 100\end{array}\right],$$

and the matrices in performance index are with

$$Q\left(k\right)=\left[\begin{array}{cc}0.001& 0\\ 0& 0.001\end{array}\right],$$

$$R\left(k\right)=1.$$

The simulation results are as follows.

Figure 8 and Figure 9 display the estimation errors of displacement and velocity for the center point of the end of the manipulator in X direction for the proposed method and OFRMPC method. It can be seen that the estimation errors of the proposed method are smaller than those of OFRMPC method and tend to zeros after 15 s.

Figure 10 displays the evolution of true displacement and velocity feasible sets for the proposed method. It shows that the feasible sets become smaller and smaller. Meanwhile, they contain the true states all the time.

Figure 11 shows the curve of the control law for the proposed method and the OFRMPC method. The upper bound $\lambda$ of the cost function at each time is shown in Figure 12 and they become small gradually. The quadratic cost function of the proposed method is 224.3114, while the quadratic cost function for OFRMPC method is 281.3348. This implies that the proposed method has better control performance.

5. Conclusions

In this paper, a robust quadratic control algorithm is proposed for the discrete-time linear systems with the ellipsoidal bounded uncertainties based on the set-membership filtering and robust optimization theory. To obtain the optimal control, the multi-step quadratic global optimal control problem is converted to multiple one-step quadratic local approximate optimal control problems. For each one-step optimal control problem, the proposed method takes the influence of the estimated states errors into consideration which ensures the optimality and robustness of the control strategy. Set-membership filtering is used to estimate the system states, which avoids solving one complex convex optimization to obtain estimated states errors sets and reduces the calculation burden. However, it needs to solve a semi-definite programming problem at each time to obtain the control law. Therefore, how to reduce the online computational burden is a key future research direction.