polytopic LDI + gevp

[Fred K]

Hi, I am a beginner with yalmip. I just came across a gevp problem but with the addition of elements to the right side of the inequalities due to the s-procedure. Here is the BMI:

I approached it with the decay rate example. I am using sdpt3-4.0 to solve it but I get lambdas NaN when I am searching for the upper bound of gamma before performing the bisection method.

F = [P1>=eye(2), P2>=eye(2),l1>=0,l2>=0,l3>=0,l1>=0,...
A1'*P1+P1*A1<=l1*P2-l1*P1-t_upper*P1,...
A2'*P2+P2*A2<=l2*P1-l2*P2-t_upper*P2,...
A1'*P1+P1*A1<=l3*P1-l3*P2-t_upper*P1,...
A2'*P2+P2*A2<=l4*P2-l4*P1-t_upper*P2];
sol = optimize(F,[],ops);

Instead if I try to assign lambdas to zero and reducing the inequalities to 2, it works. Well but then when I perform the bisection method I get negative values of lambdas even though I forced them to be greater or equal to zero.

Thanks in advance,

Fred

[Johan Löfberg]

If both lambda and P are decision variables, then you cannot use SDPT (as the problem is bilinear in lambda and P). Hence, the sol structure should tell you that the solver isn't applicable

[Johan Löfberg]

To be more precise: You do not have a GEVP, but a general bilinear semidefinite program.

[Fred K]

Thanks for you answer. Sorry for that bad piece of code. Instead of choosing an initial value for lambdas I declared lambdas as variables. Actually what I do is mutually fixing P and lambdas. In particular what I do is:

1- Calculate a suboptimal solution in this way:

A1 = [-1 -1;1 -1];

A2 = [-1 -10; 1/10 -1];

P1=sdpvar(2,2);

P2=sdpvar(2,2);

F = [P1>=eye(2),P2>=eye(2), ...

     A1'*P1+P1*A1<=-eye(2), ...

     A2'*P2+P2*A2<=-eye(2)]; 

sol=optimize(F,trace(P1)+trace(P2);

It founds something but returns Mosek error: MSK_RES_TRM_STALL () but i continued with next steps..

2)  I calculate bounds for gamma with a function in this way:

[t_lower,t_upper]=gamma_bounds(A1,A2,P1,P2) where I try to search for an upper bound:

ops = sdpsettings('solver','mosek');

%P1 and P2 must be LMI variables

% Search for lower bound

t_lower1 = max(eig(double(P1)^-1*(A1'*double(P1)+double(P1)*A1)))/2;

t_lower2 = max(eig(double(P2)^-1*(A2'*double(P2)+double(P2)*A1)))/2;

t_lower = max(abs(t_lower1),abs(t_lower2));

t_upper = t_lower*2;

% fix lambdas

l1=abs(randn(1));

l2=abs(randn(1));

l3=abs(randn(1));

l4=abs(randn(1));

sol.problem=0;

% search for upper bound until problem is unfeasible

while ~(sol.problem==1)

    t_upper = t_upper*2; 

F = [P1>=eye(2), P2>=eye(2),...

     A1'*P1+P1*A1-l1*(P2-P1)<= -t_upper*P1,...

     A2'*P2+P2*A2-l2*(P1-P2)<= -t_upper*P2,...

     A1'*P2+P2*A1-l3*(P1-P2)<= -t_upper*P2,...

     A2'*P1+P1*A2-l4*(P2-P1)<= -t_upper*P1];

    sol = solvesdp(F,trace(P1)+trace(P2),ops);

end

3) Fixing P's and optimizing lambda by performing bisection i get: gamma=3.99  all Ps like 10^-12 order and lambdas all equal to 1... Well looks like all P are almost equal to zero..sounds not so good.

while (t_upper-t_lower)>tol

  t_test = (t_upper+t_lower)/2;

  disp([t_lower t_upper t_test])

  F2=[l1>=0,l2>=0,l3>=0,l4>=0];

  F = F2+ [A1'*P1+P1*A1-l1*(P2-P1)<= -t_upper*P1,...

     A2'*P2+P2*A2-l2*(P1-P2)<= -t_upper*P2,...

     A1'*P2+P2*A1-l3*(P1-P2)<= -t_upper*P2,...

     A2'*P1+P1*A2-l4*(P2-P1)<= -t_upper*P1];

  sol = solvesdp(F,trace(P1)+trace(P2),ops);

  if ~sol.problem==0

    t_upper = t_test;

  else

    t_lower = t_test;

 end

end

gamma=t_test;

4) I tried an extra optimization with gamma fixed and lambdas fixed  but I get an unfeasible problem.

My questions are:

1- is it still applicable the bisection method to find gamma even the existence of those lambdas?

2- Do I have to perform the last optimization step? Until step 3) I get  P's though they are something like 10^-12 or -13

Thanks again.

[Johan Löfberg]

No reason to assume that a random assignment of lambda would lead to a feasible problem, and as expected, when I try to solve feasibility problem with t_lower, it is infeasible

t_lower = max(abs(t_lower1),abs(t_lower2));

t_upper = t_lower;

F = [P1>=eye(2), P2>=eye(2),...

     A1'*P1+P1*A1-l1*(P2-P1)<= -t_upper*P1,...

     A2'*P2+P2*A2-l2*(P1-P2)<= -t_upper*P2,...

     A1'*P2+P2*A1-l3*(P1-P2)<= -t_upper*P2,...

     A2'*P1+P1*A2-l4*(P2-P1)<= -t_upper*P1];

    sol = solvesdp(F,trace(P1)+trace(P2),ops)

    yalmiptime: 0.0596

    solvertime: 0.2284

          info: 'Infeasible problem (MOSEK)'

       problem: 1

Hence, your "bisection" (quotes as you treat the lambdas incorrectly) is doomed from the beginning

Your problem is bilinear in the lambdas, and bisection makes no sense unless you can solve the problem, for fixed gamma, globally over both lambda and P. You cannot do that. Hence, you cannot apply bisection in any kind of safe way. You can do hacks of course, with iterative strategies over lambda and P etc, but you have to be much more careful than what you are now (P being zero means that the problem wasn't solved correctly)

BTW, if you fix P, it makes no sense what so ever to perform bisection, as the problem then is linear in lambda and gamma

[Fred K]

The papers where I'm studying at actually show results based the algorithm on the iterative "path following method" by Hassibi -Boyd paper. Then I have to try it and eventually post here again.

Thanks again.

.

[Fred K]

I am trying the iterative solution. I try to always start with fixing the right side of the optimization as:

A1 = [-1 -1;1 -1];

A2 = [-1 -10; 1/10 -1];

P1=sdpvar(2,2);

P2=sdpvar(2,2);

%%%% Initial sub optimal problem

F = [P1>=eye(2), P2>=eye(2), ...

     A1'*P1+P1*A1<= -eye(2),...

     A2'*P2+P2*A2<= -eye(2)];

sol1=optimize(F,trace(P1)+trace(P2));

P1=double(P1);

P2=double(P2);

Then, whenever I try to fix variable lile gamma I get unfeasible problem since the beginning. The paper doesn't suggest any other constraints to be added as instead suggested in some tutorial online.

attempt#1:  FIX P1,P2, K Unfeasible

sdpvar l1 l2 l3 l4;

k=0;

sol.problem=0;

while sol.problem==0

F = [l1>=0, l2>=0, l3>=0, l4>=0, ...

     A1'*P1+P1*A1-l1*(P2-P1)<= -k*P1,...

     A2'*P2+P2*A2-l2*(P1-P2)<= -k*P2,...

     A1'*P2+P2*A1-l3*(P1-P2)<= -k*P2,...

     A2'*P1+P1*A2-l4*(P2-P1)<= -k*P1]; 

sol=optimize(F,[]);

 if sol.problem==0

     break;

 else

k=k+0.01

 end

 

end

attempt#1:  FIX P1,P2, l1,l2,l3,l4  randomly   Unfeasible

When you say iterative every variable you mean some big innested loops? or single as I am doing now but adding more constraints?

Because so far I am definitely stuck...

Thanks in advance.

[Johan Löfberg]

If you do a simple coordinate fixing strategy, why would you start to fix the right-hand side to -I? That has no connection to the original problem. More reasonable to start by fixing gamma to a small value, and thus start with -gamma*Pj. This gives you a feasible set of P. Now you fix P, and maximize gamma with l free.

[Johan Löfberg]

...or simply download penlab and see if it solves the problem.

[Fred K]

yet unfeasible..hum is it correct leaving free lamdas and gamma free ad i did?

A1 = [-1 -1;1 -1];

a=10;
A2 = [-1 -a; 1/a -1];
%Fix gamma. Optimize P
P1=sdpvar(2,2);
P2=sdpvar(2,2);
gamma=0.01;

F = [P1>=eye(2), P2>=eye(2),...

A1'*P1+P1*A1<= -gamma*P1,...
A2'*P2+P2*A2<= -gamma*P2];
sol=optimize(F,-gamma);
k=1;
while sol.problem==0

F = [P1>=eye(2), P2>=eye(2),...

A1'*P1+P1*A1<= -gamma*P1,...
A2'*P2+P2*A2<= -gamma*P2];
sol=optimize(F,-gamma);
if sol.problem==1
break
else
feasP{k,1}=double(P1);
feasP{k,2}=double(P2);
gamma=gamma+0.1;
k=k+1;
end
end
sol.problem=1;
gamma=2;
%%
% attempt 1 Fixed P's, Optimized gamma, lambdas =UNFEASIBLE.
% For fixed gamma=2.0.
%This gives you a feasible set of P. Now you fix P, and maximize gamma with l free.

sdpvar l1 l2 l3 l4;
k=0;

feas=0;
j=1;
gamma=0.01;
while j<max(size(feasP))
P1=feasP{j,1};
P2=feasP{j,2};

F = [l1>=0, l2>=0, l3>=0, l4>=0,...

A1'*P1+P1*A1-l1*(P2-P1) <= -gamma*P1,...
A2'*P2+P2*A2-l2*(P1-P2) <= -gamma*P2,...
A1'*P2+P2*A1-l3*(P1-P2) <= -gamma*P2,...
A2'*P1+P1*A2-l4*(P2-P1) <= -gamma*P1];
sol=optimize(F,-gamma);
feas2(j)=sol.problem;
j=j+1;
end
% attempt 2 Fixed P's, Fixed lambdas, Optimize gamma =UNFEASIBLE.
%This gives you a feasible set of P. Now you fix P, and maximize gamma with l free.
%%
sdpvar gamma;
l1=abs(randn(1));
l2=abs(randn(1));
l3=abs(randn(1));
l4=abs(randn(1));
k=0;
feas=0;
j=1;
while j<max(size(feasP))-1
P1=feasP{j,1};
P2=feasP{j,2};
F = [gamma>=0,...
A1'*P1+P1*A1-l1*(P2-P1) <= -gamma*P1,...
A2'*P2+P2*A2-l2*(P1-P2) <= -gamma*P2,...
A1'*P2+P2*A1-l3*(P1-P2) <= -gamma*P2,...
A2'*P1+P1*A2-l4*(P2-P1) <= -gamma*P1];
sol=optimize(F,-gamma);
feas2(j)=sol.problem;
j=j+1;
end

[Johan Löfberg]

Lets start with the obvious strangeness

> A1 = [-1 -1;1 -1];
> a=10;
> A2 = [-1 -a; 1/a -1];
> %Fix gamma. Optimize P
> P1=sdpvar(2,2);
> P2=sdpvar(2,2);
> gamma=0.01;
> F = [P1>=eye(2), P2>=eye(2),...
> A1'*P1+P1*A1<= -gamma*P1,...
> A2'*P2+P2*A2<= -gamma*P2];
> sol=optimize(F,-gamma);
> k=1;
> while sol.problem==0
> F = [P1>=eye(2), P2>=eye(2),...
> A1'*P1+P1*A1<= -gamma*P1,...
> A2'*P2+P2*A2<= -gamma*P2];
> sol=optimize(F,-gamma);

gamma is constant, so why send it as objective??

> if sol.problem==1
> break
> else
> feasP{k,1}=double(P1);
> feasP{k,2}=double(P2);
> gamma=gamma+0.1;
> k=k+1;
> end
> end

At this point, P1 and P2 do not contain any feasible solutions. You just exited from the loop since the problem was infeasible, hence no reason to believe that there are any good solutions available in P1 and P2 (if there are any values, they must be infeasible)

[Fred K]

Ok thanks, I adjusted like as follows:

%Fix gamma. Optimize P 

P1=sdpvar(2,2); 

P2=sdpvar(2,2); 

gamma=0; 

F = [P1>=eye(2), P2>=eye(2),... 

A1'*P1+P1*A1<= -gamma*P1,... 

A2'*P2+P2*A2<= -gamma*P2]; 

sol=optimize(F,[]); 

k=1; 

feasP{k,1}=double(P1); 

feasP{k,2}=double(P2);

while sol.problem==0  

        

    gamma=gamma+0.2

    F = [P1>=eye(2), P2>=eye(2),... 

    A1'*P1+P1*A1<= -gamma*P1,...    

    A2'*P2+P2*A2<= -gamma*P2]; 

    sol=optimize(F,[]); 

    if sol.problem==1 

        disp(['Unfeasible with gamma=' num2str(gamma)])

            gamma=gamma-0.2  

            break;

    else 

    feasP{k,1}=double(P1); 

    feasP{k,2}=double(P2); 

    k=k+1; 

      end 

end 

Executing this code it returns about 10x2cells of feasibles P's with maximum value for gamma equal to 2.0.

Well, the rest is unfeasible as soon as I add more constraints with lambdas..

[Johan Löfberg]

Is that a question or are you saying you are done now?

[Fred K]

Err, sorry, actually it is a question. The paper where I took the examples show  feasible solution. What happen to me is that after I execute that part of code I perform this optimization as you kindly suggested: using all the P's found before, i optimize gamma and lamda as follow

%% 

% attempt 1 Fixed P's, Optimized gamma, lambdas =UNFEASIBLE. 

% For fixed gamma=2.0. 

sdpvar gamma l1 l2 l3 l4;  

j=1;  

while j<=max(size(feasP)) 

P1=feasP{j,1}; 

P2=feasP{j,2}; 

F = [l1>=0, l2>=0, l3>=0, l4>=0,... 

A1'*P1+P1*A1-l1*(P2-P1) <= -gamma*P1,... 

A2'*P2+P2*A2-l2*(P1-P2) <= -gamma*P2,... 

A1'*P2+P2*A1-l3*(P1-P2) <= -gamma*P2,... 

A2'*P1+P1*A2-l4*(P2-P1) <= -gamma*P1]; 

sol=optimize(F,-gamma); 

feas2(j)=sol.problem;

feas_gamma(j)=double(gamma);

j=j+1; 

end 

Well, removing constraints over gamma, I get NEGATIVE gamma and FEASIBLE solution.

If I constraint gamma to be positive -> Gamma is negative.

If I remove gamma from the objective function -> Yalmip returns Unbounded objective function (sol.problem == 2)

Should I add extra constraints? The paper doesn't mention it actually..

[Johan Löfberg]

'Does the paper say how they solve this *nonconvex* problem (i.e. exactly, not just the term path-following). If not, there is no way for you to reproduce the results consistently. You are using one way you've invented, and it apparently doesn't work. The Ps you compute in the first limited model are not useable in the second run

Why don't you use penlab? or a more advanced approach where you linearize the model etc, trust regions etc. The two-step method you use now seems very odd, no reason for original P to be feasible when you add two new constraints to the model