Using SDP to solve an optimization problem subject to LMIs

[Allez-y]

Dear all,

I have worked on a specific optimization problem subject to LMIs with strict constraints. I think SDP is the most recommended solution. To test my code, I have used an example whose results are given. Although I have read some YALMIP's tutorials and examples, I have made mistakes. Could anybody help me?

To clarify my problem, I attached an image, where I present (from left to right) my problem formulation, the corresponding optimization problem, and a numerical example with expected results. The optimization problem can be solved using the LMI (A) or (B), since they are equivalent. Anyway, an incovenient appears, to represent strict constraints using YALMIP. For matrices P>0, the SDP tutorial recommends to use P >= I or enforce trace(P)==1. However, these led to infeasible solutions using FMINCON as basis solver. I thought suitable to enforce trace(P)=n, where "n" is the square matrix order, since trace(I)=n. So doing, for each LMI I have found different and wrong results. My code is presented as follows.

%--------------------------- Begin ---------------------------%

clear all

clc

% Discrete-time linear system:

% xk = A*xk_1 + E*wk_1

% yk_1 = C*xk_1 + F*vk_1

% Test example

A = [0.7846 -0.0154 ; 0.6056 0.9983];

C = eye(2);

E = [-85e-4 -6e-4 ; -0.0603 2e-4];

F = 0.01*eye(2);

n = 2; % number of states

m = 2; % number of measurements

% Solving the optimization problem

%%%% Variables %%%%

gamma = sdpvar(1);

T = sdpvar(n,n,'full');

N = sdpvar(n,m,'full');

L = sdpvar(n,m,'full');

P = sdpvar(n,n,'full');

% Discrete-time error system:

% ek = Ae*ek_1 + Be*dk_1

Ae = T*A - L*C;

Be = [T*E -L*F -N*F];

%%% Constraints %%%

% I want to approximate gamma>0, P>0 and LMI<0 using SDP, since strict constraints are not possible

% First option: LMI (A)

LMI = [(Ae'*P*Ae + eye(n) - P) Ae'*P'*Be ; Be'*P*Ae (Be'*P*Be -gamma^2*eye(n + 2*m))];

% Second option: LMI (B)

%LMI = [(eye(n) - P) zeros(n,n+2*m) Ae'*P' ; zeros(n+2*m,n) -gamma^2*eye(n + 2*m) Be'*P' ; P*Ae P*Be -P];

const = [gamma >= 1e-8 , P >= 0 , trace(P) == n , LMI <= 0 , trace(LMI) == -size(LMI,1)]; % constraints

%%%% Objective function %%%%

goal = gamma^2; % to minimize

%%%% Solution %%%%

sol = optimize(const , goal);

%%%% Analyze error flags %%%%

if sol.problem == 0

 % Extract and display value

 value(T)

 value(L)

 value(N)

 value(gamma)

else

 display('Hmm, something went wrong!');

 sol.info

 yalmiperror(sol.problem)

end

%-------------------------------- End ------------------------------------%

Where could I be making mistakes?

[Johan Löfberg]

This is not an LMI as it is nonlinear in gamma. However, you can trivially make it linear by replacing gamma^2 with  new variable t

[Allez-y]

Thanks for expliciting the variable change in need. In this case, could you indicate a way to solve my problem? Should I use another solver, or fix my code?

[Johan Löfberg]

1. As said, don't use gamma^2, simply call that expression t, and minimize t instead. Equivalent problem, but linear.

2. You've failed to create symmetric matrices. If you look at e.g. the variable LMI, it is not symmetric, i.e. trivial bugs in your code. Hence your model is not a semidefinite programming problem at the moment, but a standard nonlinear nonconvex program

3. The model is extremely hard still (after fiing bugs above), as you have polynomial terms such as e.g. Ae'*P*Ae. You've got a lot of work to do to get something you will be able to solve. As you've called the variable LMI, it smells like you don't understand what the L stands for, and how utterly important that is.

[Allez-y]

Dear Johan,

Finally, I got to implement the first example, after applying change of variables, defining t = sdpvar(1) and using it instead of gamma^2 like you said. Thereby, I obtain an LMI, and I can use SDP-SeDuMi to solve my optimization problem. I have adjusted the value gamma to 0.0099, as desired, but other matrices have been returned. This occurs because I may use a different strategy to avoid null matrices, i.e., I handle trace(*) == x, but I can use another strategy to reach the same gamma, and other matrices can be returned. Anyway, the state estimation leads to results very close.

I have two more doubts, related to the new attached image. The black box at left is another way to solve the first example. This representation corresponds to the change of variables in need to avoid nonlinear products in the matrix inequality. I applied this configuration to solve my first example, and I had success. Using the same code structure to solve the second example, I have obtained NaN for the matrix Y. There is a massage that appears in the beginning of the solution:

"The coefficient matrix is not full row rank, numerical problems may occur".

Could you give me a tip to fix my following code?

%%%%%%%%%%% Begin %%%%%%%%%%%%%%%

clear all

close all

clc

%--------- Discrete-time linear system --------%

% xk = A*xk_1 + E*wk_1

% yk_1 = C*xk_1 + F*vk_1

A = [-0.54 0.45 0.36 0 0 ;

    0.63 0.45 0.18 0.36 0 ; 

    0.09 0.45 0.27 0.09 0.18 ;

    0 0 0.25 0.25*sqrt(2) -0.25*sqrt(2);

    0 0 0 0.25*sqrt(2) 0.25*sqrt(2)];

C = [1 0 0 0 0 ; 0 0 0 1 0];

E = [-1 0 0 0 1]';

F = 0.1*eye(2);

n = 5; % number of states

m = 2; % number of measurements

nw = 1; % number of process noise parcel

nv = 2; % number of measurement noise parcel

%------------- Initializing variables --------------%

TT = [eye(n) ; C]; alp1 = [eye(n) ; zeros(m,n)]; alp2 = [zeros(n,m) ; eye(m)];

PTT = (TT'*TT)\TT'; % pseudo inverse

Psi = eye(n+m) - TT*PTT;

t = sdpvar(1);

W = sdpvar(n,m,'full');

Y = sdpvar(n,n+m,'full');

P = sdpvar(n,n,'sym');

O1 = P*PTT*alp1*A + Y*Psi*alp1*A - W*C;

O2 = P*PTT*alp1*E + Y*Psi*alp1*E;

O3 = -W*F;

O4 = -P*PTT*alp2*F - Y*Psi*alp2*F;

%---------------- Constraints -----------------%

LMI = blkdiag(eye(n)-P , -t*eye(nw), -t*eye(nv), -t*eye(nv), -P);

LMI(end+1-n:end,1:n) = O1;

LMI(1:n,end+1-n:end) = O1';

LMI(end+1-n:end,n+1:n+nw) = O2;

LMI(n+1:n+nw,end+1-n:end) = O2';

LMI(end+1-n:end,n+nw+1:n+nw+nv) = O3;

LMI(n+nw+1:n+nw+nv,end+1-n:end) = O3';

LMI(end+1-n:end,n+nw+nv+1:n+nw+2*nv) = O4;

LMI(n+nw+nv+1:n+nw+2*nv,end+1-n:end) = O4';

ops = sdpsettings('solver','sedumi','sedumi.eps',1e-12,'verbose',2);

const = [sqrt(t) >= 1e-5 , P >= 0 , trace(P) == 11.0025 , LMI <= 0 , trace(LMI) == -17.17];

%---------------- Objective function ------------------%

goal = t; % to minimize

%------------------ Solution ----------------%

sol = optimize(const , goal, ops);

%----------------- Analyze error flags -----------------%

if sol.problem == 0

    T = PTT*alp1 + value(P)\value(Y)*Psi*alp1

    N = PTT*alp2 + value(P)\value(Y)*Psi*alp2

    L = value(P)\value(W)

    gamma = sqrt(value(t))

else

    display('Hmm, something went wrong!');

    sol.info

    yalmiperror(sol.problem)

end

%%%%%%%%%% End %%%%%%%%%%%%%%%%%

My last question is related to the change of variables. The matrix inequality at the second attached image is, in fact, the matrix inequality (B), at the first attached image, after doing the change of variables. However, this equivalence is not direct, because the pseudo inverse. I have done a result comparison for the first example solving: (I) the matrix inequality (B) after applying the change of variables directly, i.e., defining Ua = P*L, Ub = P*T, Uc = P*N and t = gamma^2; (II) the inequality described in second image.

The solution provided by (I) is not right, but the solution from (II) is.

Do you know why it occurs? For me, they should be the same solution.

Thanks.

[Allez-y]

Dear all,

I got to solve my all doubts. Really, the true solution is obtained solving the optimization problem with the matrix inequality attached latter (second image), because my state observer design. If I use the other matrix inequalities, the results are close to zero, since I would be considering the matrices T and N as independent variables, but they are not. This is the reason the last matrix inequality to be different, since it includes the dependency between both matrix variables.

About the solver, I have used another toolbox, which treats by itself the LMI by means of SDP. Then, I do not need to handle the symmetric matrices (by means of traces or identity matrices) until to obtain a feasable solution, which is non null.

Thanks for everything Johan, your help was very useful. I am so grateful to you by the discussions.

[Johan Löfberg]

Not sure what you are trying to say or convey.

Writing sqrt(t)>=1e-5 is an unnecessarily complicated way (as it requires an SOCP representation) of saying t >= 1e-10, which effectively in any numerical solver is the same as saying t>=0. However, since you are minimizing t, and the optimal t is far from zero, why even constrain it...

Many of the variables in Y are never used in the optimization problem, and thus they will be reported as having value nan (as no value is computed). You should either assign(Y,0) before solving the problem to set them, or do Y = value(Y);Y(isnan(Y))=0 after solving.