Discrete-time LQR and solutions via LMI

convex optimizationlinear-matrix-inequalitymatricesmatrix equationsschur-complement

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?

Best Answer

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
Related Question