Using for loop to define multiple constraints

[benny]

Hi, 

I'm new to YALMIP but I already can see how great it is.

I have a problem that i need to solve which is basically semidefinite problem.

My problem looks like the next form:

for i=1:1:100

       A(i) = ...

       -X(i+1) + X(i) + dt * (X(i) * A(i)' + A(i) * X(i)+L(i)) <0

end

That is true if X,A,L were just variables and I could make vector of number.

However here X,A,L are matrices and their value will be different everytime because of A changing every iteration.

Is there a way to define a structure of LMI's using for loop that eventually i will get the answers for L,X for each iteration?

Thank you very much.

Benny 

[Johan Löfberg]

It is very unclear what you want to do, or what the problem is

[benny]

I had like to make for example 100 constraints from the following form using a loop:

Constraint{1] = [-X{2}+X{1}+dt*(X{1}*A{1}'+A{1}*X{1}+B{1}*L{1})<0, X{1} >= inv(R), X{1} < inv(Gamma)];

Constraint{2} = [-X{3}+X{2}+dt*(X{2}*A{2}'+A{2}*X{2}+B{2}*L{2})<0, X{1} >= inv(R), X{2} < inv(Gamma)];

Constraint{3} = [-X{4}+X{3}+dt*(X{3}*A{3}'+A{3}*X{3}+B{3}*L{3})<0, X{1} >= inv(R), X{3} < inv(Gamma)];

....

and so on.

Basically im going to have many X sdpvar's which everytime depands on the one before.

So im asking if there is a way to do it in for loop because everytime i will have different number of constraints and not constant 100.

Eventually I want to get an LMI of N Constaints and to solve it.

Thank you.

[Johan Löfberg]

constraints = []

for i = 1:100

 constraints = [constraints, morestuff];

end

[benny]

but does it work when defining an sdpvar of the form X{1}, X{2}... as X1,X2 and so on?

[Johan Löfberg]

X = sdpvar(repmat(n,1,M),repmat(n,1,M)) generates a cell of M sdpvars of size nxn

https://yalmip.github.io/tutorial/basics/

[benny]

Thank you Johan for your help. 

I tried what you said and its great.

However i just wonder why do i get different results if i set the following:

option 1:

X = sdpvar(repmat(2,1,T/dt+1),repmat(2,1,T/dt+1));

L = sdpvar(repmat(1,1,T/dt),repmat(2,1,T/dt));

Constraints =[];

for k=1:1:T/dt

Constraints = [Constraints, -X{k+1}+X{k}+dt*(X{k}*A{k}'+A{k}*X{k}+B{k}*L{k})<0, X{k} < inv(Gamma),X{1}>= inv(R)];

end

option 2:

Constraints =[];

for k=1:1:T/dt

X{k}= sdpvar(2,2);

X{k+1}=sdpvar(2,2);

L{k}{sdpvar(1,2);

Constraints = [Constraints, -X{k+1}+X{k}+dt*(X{k}*A{k}'+A{k}*X{k}+B{k}*L{k})<0, X{k} < inv(Gamma),X{1}>= inv(R), X{1}>= inv(R)];

end

[Johan Löfberg]

Because in the second case you create in total 2*(T/dt) X variables (you create two new ones in the loop. Since you redefine them, there is no connection between different k, i.e, X(k+1) will not be equal to X(k) when you increase k. Alternatively explained, the second code is 100% equivalent to

Constraints =[];

for k=1:1:T/dt

C = sdpvar(2,2);

D=sdpvar(2,2);

L{k}{sdpvar(1,2);

Constraints = [Constraints, -D+C+dt*(C*A{k}'+A{k}*C+B{k}*L{k})<0, C < inv(Gamma),C>= inv(R), C>= inv(R)];

end

[benny]

Thank you very much, helped me a lot.

[benny]

I'm sorry but i'm still having problems with defining the constraint structure of this differential LMI.

my code looks as the following:

X = sdpvar(repmat(2,1,T/dt+1),repmat(2,1,T/dt+1));

L = sdpvar(ones(1,T/dt),repmat(2,1,T/dt));

Constraints = [];

for k = 1:T/dt

        

    A{k} = [2/(tf1-t) 0 ; 0 2/(tf2-t)]; 

    B{k} = [-1/(Vc*(tf1-t));-1/(Vc*(tf2-t))];

    

    Constraints = [Constraints, -X{k+1}+X{k}+dt*(X{k}*A{k}'+B{k}*L{k})<=0, X{k} <= inv(Gamma),X{1}>= inv(R)];

    t=t+dt;    

end

options = sdpsettings('Solver','sedumi');

solvesdp(Constraints,[],options);

When inserting the results to check if all the constraint worked, some of them arent and im getting "run into numerical problems" message even that i know there is an answer for the test case im running.

Thank you for your help Johan! 

[Johan Löfberg]

https://yalmip.github.io/faq/solutionviolated/

Note that you don't have to compute eigenvalues manually etc to check if a solution is feasible, simply do use check, check(Constraints)

[benny]

changed my constraints into:

 Constraints = [Constraints, -X{k+1}+X{k}+dt*(X{k}*A{k}'+B{k}*L{k})<=-eye(2)*0.5, X{k} <= inv(Gamma),X{1}>= inv(R)];

I dont get the message of iter problems but still not all of the constraints active in the end when checking them.

What concerns me the most is the in each loop im making 3 constraints while it seems to me like the last one X{1}>= inv(R) basically can be outside of the loop because its const.

Is there any right way to write those constraints to make it work/ change properties of solver?

Thank you.

[Johan Löfberg]

of course it should be outside the loop, you're adding the same constraint multiple times now

why do you want all constraints to be active? if the optimal solution is in the interior on some, so what? redundant constraints simply

[benny]

But if the X solutions im looking for should all be <inv(Gamma) and its not the case because at like half of the iterations X>inv(gamma) [after checking], then isnt it meaning that the X im getting is not the one im looking for?

[Johan Löfberg]

Huh? I'm talking about X{1}>= inv(R) which should be outside the loop (as it  doesn't depend on k)

what does check(Constraints) display

[benny]

once I took this constraint out im having now 2001 constraints and the run into numerical problems message.

Now im taking the X{k},L{k} calculated and inserting them into the constraints to check with logical (0 or 1) to see if it really apply the constraint and draw it.

My graph for the constraints looks like this:

Im thinking all of the figures should be constantly showing 1 (constraint satisfied)

check(constraints) display:

and so on till #2001.

I really appreciate your help.

[Johan Löfberg]

The solver you have used here has failed miserably (and I presume it has told you so). That is a clearly infeasible solution (primal residuals around -1e-8 and such is normal, but you have much larger infeasibilities. A primal residual of -0.5 means the solver has given up completely)

[Johan Löfberg]

Is your data reasonable, i.e., what is inv(Gamma) and inv(R), typical values of A and B etc?

[benny]

The test case im running (I know it has a solution) is using the next data:

and my code looks like this:

[Johan Löfberg]

why don't you post the code instead of a picture of the code...

inv(R) will essentially be 0 for the solver, and it looks like the solver will have to squeeze in X between two very small matrices, leading to a very ill-conditioned problem

[benny]

t=0;

dt = 0.01;

R = [7.8*10^5 0 ; 0 6.3*10^7];

Gamma = [62500 0 ; 0 62500];

X = sdpvar(repmat(2,1,T/dt+1),repmat(2,1,T/dt+1));

L = sdpvar(ones(1,T/dt),repmat(2,1,T/dt));

Constraints = [X{1}>= inv(R)];

for k = 1:T/dt

        

    A{k} = [2/(tf1-t) 0 ; 0 2/(tf2-t)]; 

    B{k} = [-1/(Vc*(tf1-t));-1/(Vc*(tf2-t))];

    

    Constraints = [Constraints, -X{k+1}+X{k}+dt*(X{k}*A{k}'+A{k}*X{k}+B{k}*L{k}+L{k}'*B{k}')<=-eye(2)*0.00005 , X{k} <= inv(Gamma)];

    t=t+dt;    

end

options = sdpsettings('Solver','sedumi');

solvesdp(Constraints,[],options);

Is there any solver that can deal with this small numbers of inv(R) to avoid the ill-conditioned case? all of the numbers will be pretty small in this problem.

[Johan Löfberg]

Undefined function or variable 'T'.

 

Undefined function or variable 'tf1'.

 

[benny]

t=0;

dt = 0.01;

tf1 = 15;

tf2 = 15;

T = 10;

Vc = 5000;

R = [7.8*10^5 0 ; 0 6.3*10^7];

Gamma = [62500 0 ; 0 62500];

X = sdpvar(repmat(2,1,T/dt+1),repmat(2,1,T/dt+1));

L = sdpvar(ones(1,T/dt),repmat(2,1,T/dt));

Constraints = [X{1}>= inv(R)];

for k = 1:T/dt

        

    A{k} = [2/(tf1-t) 0 ; 0 2/(tf2-t)]; 

    B{k} = [-1/(Vc*(tf1-t));-1/(Vc*(tf2-t))];

    

    Constraints = [Constraints, -X{k+1}+X{k}+dt*(X{k}*A{k}'+A{k}*X{k}+B{k}*L{k}+L{k}'*B{k}')<= -eye(2)*0.5 ; X{k} <= inv(Gamma)];

    t=t+dt;    

end

options = sdpsettings('Solver','sedumi','sedumi.maxiter',200);

solvesdp(Constraints,[],options);

[Johan Löfberg]

It's clearly infeasible, already this fails

Constraints = [X{1}>= inv(R)];

for k = 1:2

        

    A{k} = [2/(tf1-t) 0 ; 0 2/(tf2-t)]; 

    B{k} = [-1/(Vc*(tf1-t));-1/(Vc*(tf2-t))];

    

    Constraints = [Constraints, -X{k+1}+X{k}+dt*(X{k}*A{k}'+A{k}*X{k}+B{k}*L{k}+L{k}'*B{k}')<= -eye(2)*0.5 ; X{k} <= inv(Gamma)];

    t=t+dt;    

end

solvesdp(Constraints)

[benny]

from your experience, is it because of defining wrong the constraints of the following problem in my code:

 or because of the numbers used/solver?

[Johan Löfberg]

it's clearly your constraints. it is easily detected as infeasible

back to the drawing board and develop a better theory