Discrete-time LQR and solutions via LMI

Question

Having a infinite horizon discrete-time LQR problem

$J^* = \min_u \ \sum_{k=0}^{\infty} x^\top_k Q x_k +u^\top_kRu_k$

subject to

$x_{k+1}= Ax_k+Bu_k, \quad x(0)=x_0$.

With some algebra manipulations, and setting $J^*=x_kPx_k$, with $P=P^\top\succ 0$ the following LMI is obtained:

$J^* = \begin{bmatrix} u_k\\x_k \end{bmatrix}^\top \begin{bmatrix} B^\top PB+R & B^\top PA\\A^\top PB & A^\top PA +Q \end{bmatrix} \begin{bmatrix} u_k\\x_k \end{bmatrix} \prec 0$

Taking the Schur complement, the resulting state feedback controller $u_k=Kx_k$ is

$K = -(R+B^\top PB)^{-1}B^\top PA$

where P is the solution of the Riccati equation

$P = Q + A^\top PA - A^\top PB(R+B^\top PB)^{-1}B^\top PA$

I implemented an example in Matlab and compared the solutions obtained using the command dlqr and the LMI solved with Yalmip, but the values of the obtained (P,K) are not the same.

the code is:

A=[1.1 -0.3; 1 0];
B=[1;0];
Q=eye(2,2);
R=1;

%% inf horizon
[K,S,e] = dlqr(A,B,Q,R)

%% LMI solver:
P = sdpvar(2,2);
F1 = [A'*P*A + Q,A'*P*B; B'*P*A, R+B'*P*B ];
F = [P>=0, F1<=0];
optimize (F,P);
Pfeasible=value(P);
Kfeasible = -inv(R+B'*Pfeasible*B)*B'*Pfeasible*A

Do you have an idea why? is it due to the solver?

Answer 1

reference: Connections between duality in control theory and convex optimization, V Balakrishnan, L Vandenberghe . The LMI should be formulated as follows

$\max \ \text{trace}(P)$

sbj to

$P\succeq 0$

$\begin{bmatrix} A^TPA+Q-P & A^TPB\\ B^TPA& R+B^TPB \end{bmatrix}\succeq 0$

Solved in Yalmip yields the following code

P = sdpvar(2,2);
mat = [A'*P*A+Q-P, A'*P*B; B'*P*A, R+B'*P*B ];
cons = [P>=0, mat>=0];
optimize(cons, -trace(P));
Pf = value(P)
Kf = inv(R+B'*Pf*B)*B'*Pf*A

If you run the code you should find a feasible solutions for $P$ and $K$.

The code in cvx should be

cvx_begin
variable P(2,2) symmetric
maximize( trace(P) )
subject to
    [A'*P*A+Q-P, A'*P*B; B'*P*A, R+B'*P*B ]>=0;
    P>=0;
cvx_end

Answer 2

First, a question like this is much more suitable for the YALMIP google groups support forum.

The LMI you have used is not correct. If you take the known true solution, it is far from feasible in the model you have

>> assign(P,S);
>> check(F)

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   ID|          Constraint|   Primal residual|   Dual residual|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|   #1|   Matrix inequality|            1.0377|             NaN|
|   #2|   Matrix inequality|           -7.5552|             NaN|
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

In fact, the problem you have posed is infeasible (the command optimize is telling you this if you look at the output)

Finally, the objective you use makes no sense. You probably meant $\text{trace}(P)$ or some other scalar measure.

Answer 3

Good day, I had the same issue few days ago (LQR and LMI). Here is how I convinced myself about it. I had to reformulate the problem from the beginning. I hope it helps.

Let $J$ be the cost function, \begin{align} J = \sum_{k=0}^{\infty} x_k^T Qx_k + u_k^T Ru_k = \sum_{k=0}^{\infty}z_k^Tz_k \end{align} with $z_k= C x_k+Du_k$, $C=\begin{bmatrix} Q^{1/2}\\ 0\end{bmatrix}$ and $D=\begin{bmatrix} 0\\R^{1/2}\end{bmatrix}$. Solving the recurrence relation of $x_{k+1}=Ax_{k}+Bu_k$ with $u_k=Kx_k$ and $x\in\mathbb{R}^n$ one obtain, \begin{align} z_k = (C+DK)(A+BK)^kx_0, \end{align} evaluating it in the cost function, we obtain \begin{align} J = \sum_{k=0}^{\infty}\left\lVert z_k\right\lVert_2^2 =\sum_{k=0}^{\infty}\left\lVert (C+DK)(A+BK)^kx_0\right\lVert_2^2. \end{align} As we need to solve for $K$, we have to deal with matrix norms. Therefore, \begin{align} J &= \sum_{n=0}^{\infty} x_0^T ((A+BK)^T)^k(C+DK)^T (C+DK)(A+BK)^kx_0\\&=\sum_{n=0}^{\infty} tr\left( (C+DK)(A+BK)^kx_0x_0^T ((A+BK)^T)^k(C+DK)^T \right) \end{align} Thus, we arrive to the controllability gramian, and the problem becomes, \begin{align} J = tr\left( (C+DK)P(C+DK)^T\right),\\ P-(A+BK)P(A+BK)^T-x_0x_0^T>0. \end{align} Applying Schur complement and other manipulations we get the optimization problem,

\begin{align} &\min \quad tr(W)\\ s.t.& \quad \begin{bmatrix} W & CP+DF\\ *&P\end{bmatrix}>0\\ &\quad \begin{bmatrix} P & AP+BF & x_0\\ *&P & 0\\x_0^T & 0 & \mathbb{I}\end{bmatrix}>0 \end{align} where $F=KP$ and $P=P^T>0$. Testing the system under interest, python cvxpy.

For LQR

P_ric = solve_discrete_are(Ad,Bd,Qd,Rd)
L_ric = -(Rd+Bd.T*P_ric*Bd).I*Bd.T*P_ric*Ad
J_cost = x0.T*P_ric*x0 

it returns $L_{lqr} = \begin{bmatrix}-0.75265207 & 0.22140237\end{bmatrix}$ and $J_{lqr} = 2.81578553$.

For LMI

P = cp.Variable((nx,nx))
W = cp.Variable((nx+1,nx+1))
F = cp.Variable((nu,nx)) 
gamma = cp.Variable((1,1))
Cz = np.append(Qd_sqrt,np.zeros((nu,nx)),axis=0)
Dz = np.append(np.zeros((nx,nu)),Rd_sqrt,axis=0)
zero_aux = np.zeros((nx,1))
gamma = cp.Variable((1,1))
id_aux = np.matrix([[1.]])
LMI0 = cp.trace(W)
LMI1 = cp.bmat([[ W, Cz*P+Dz*F], [(Cz*P+Dz*F).T, P]])
LMI2 = cp.bmat([[P, Ad*P+Bd*F,x0],[ (P*Ad.T+F.T*Bd.T), P,zero_aux],  [x0.T,zero_aux.T ,id_aux]])
cons_a = P == P.T
cons_b = W == W.T
cons_e = gamma >>0
cons0 = LMI0 << gamma
cons1 = LMI1 >> 0
cons2 = LMI2 >> 0
#
list_const= [cons_a,cons_b,cons_c,cons0,cons1,cons2]
optprob = cp.Problem(cp.Minimize(gamma), constraints=list_const)
result = optprob.solve()

it returns $L_{lmi} = \begin{bmatrix}-0.75255735 & 0.22138669\end{bmatrix}$ and $J_{lmi} = 2.81578088$.

And the solution is the same. After all this analysis, the LQR formulation through LMIs seems to be not so obvious. What I believe, but I'm not completely sure, is that as the stationary Riccati equation is solvable for the only stabilizing solution, it may not find a solution through inequalities. But again, I'm not sure about this, and I'm still trying to find other argument.