Re: Keep getting NaN

[Johan Löfberg]

Have you even looked at the diagnostic from optimize (does not look like that in your code)? The problem is not being solved (as told by the diagnostic code) as you don't have any solver as you have a bilinear SDP as F1 is bilinear which comes from B1 being bilinear which comes from F being bilinear which comes from the bilinear term E*S. To me though, it looks like you think you've made a variable change (which is possible here), i.e. easily fixed if you read the https://yalmip.github.io/badsdps with your mistake in "That is not a variable change" 

torsdag 10 juni 2021 kl. 23:06:55 UTC+2 skrev yu.ti...@gmail.com:

Dear Johan,

i'm keep getting NaN and don't really know why. The code is as follow:

clear all;

clc;

option = sdpsettings('solver','mosek','verbose',0);

%some defition about the system

Ap=[-1 0.3; 0 0.5];

Bpu=[0.5;-10];

Bpw=[1;0];

Cpy=[0.4  0.8];

Cpz=[1 0;0 1;0 0];

Dpz=[0;0;1];

Ac=[-4 1;0 -8];

Bc=[-1;-0.5];

Cc=[0.3 -2];

Dc=[1];

H=4.0928;

n=2;

m=1;

p=2;

delta=4;

%end 

W=sdpvar(n+m+p,n+m+p);

S=diag(sdpvar(2,1,'full'));

E=sdpvar(n,2,'full');

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

F=E*S;

mu=sdpvar(1,1);

v0=eye(2,2);

gama=sdpvar(1,1);

A=[Ap Bpu zeros(n);

    H*Dc*Cpy  -H  H*Cc;

    Bc*Cpy  zeros(2,1) Ac];

B11=[Bpu  zeros(2,1);

    zeros(1)    eye(1)];

B1=[B11*S;F];

B2=[Bpw;zeros(3,1)];

C2=[Cpz  Dpz  zeros(3,2)];

D2=[-Dpz  zeros(3,1)];

G_hat=[zeros(1,2)  eye(1)  zeros(1,2);

    H*Dc*Cpy  -H  H*Cc];

F1=[W*A'+A*W    -B1+Y'   B2  W*C2';

    -B1'+Y    -2*S    zeros(2,1) S*D2';

    B2'     zeros(1,2)     -eye(1) zeros(1,3);

    C2*W    D2*S        zeros(3,1) -gama*eye(3,3)]<=0;

h=(G_hat*W-Y);

F2=[W  h'; h  mu*(v0)^2]<=0;

F3=[mu-delta]<=0;

F=[];

F=[F,W>=0];

F=[F,S>=0];

F=[F,gama>=0];

F=[F,mu>=0];

F=[F,F1];

F=[F,F2];

F=[F,F3];

sol = optimize(F, [], option);

gamaa=value(gama)

Ww=value(W)

Ss=value(S)

Gg=value(G_hat)

Yy=value(Y)

Ee=value(E)

Will be realy greatfull for any of you help.

Best Regards

Yu

[Johan Löfberg]

As explained on the Wiki article, this is a classical variable-change where you simply define a new variable to represent E*S (since E isn't used anywhere else), solve the problem, and then recover

F = sdpvar(...)

solve

E = value(F)*inv(value(S))

fredag 11 juni 2021 kl. 10:53:02 UTC+2 skrev yu.ti...@gmail.com:

Thanks for your reply! I'm a total beginner.  What do i have to do, if i really have to solve a problem with a term like F=E*S or even F=E*S*S, when E and S are both unknown variables and F is a term of a bigger constraint matrix?

[Johan Löfberg]

note though that you are using the variable name F twice in different context, so adjustments have to be made for the recovery to actually run

Also, the code will still fail as the problem is infeasible. Already this is infeasible

optimize([F1,F2])

  struct with fields:

    yalmipversion: '20210331'

    matlabversion: '9.9.0.1524771 (R2020b) Update 2'

       yalmiptime: 0.1212

       solvertime: 0.4048

             info: 'Infeasible problem (learn to debug) (MOSEK-SDP)'

          problem: 1

fredag 11 juni 2021 kl. 10:53:02 UTC+2 skrev yu.ti...@gmail.com:

Thanks for your reply! I'm a total beginner.  What do i have to do, if i really have to solve a problem with a term like F=E*S or even F=E*S*S, when E and S are both unknown variables and F is a term of a bigger constraint matrix?

在2021年6月11日星期五 UTC+2 上午7:23:35<Johan Löfberg> 写道：

[Johan Löfberg]

and infeasibility is obvious

A is stable, hence A*W+W*A' negative definite which you enforce in F1 implies W is positive definite, which is inconsistent with F2 where you enforce W to be negative definite

[Johan Löfberg]

Well you cannot get (actual) results the problem is infeasible (unless you fixed that bug, I presume F2 is supposed to be >=0)

Your talk about F*S F*S^2 makes no sense. In your model you have E*S, and F is the variable change replacing that expression

fredag 11 juni 2021 kl. 12:28:51 UTC+2 skrev yu.ti...@gmail.com:

That was really big help from you, thanks a lot. I think i'm getting there. The coding are running and i'm getting results, but different than the result from the paper(where the whole LMIs come from.)

In original LMI, there is a term -B1*S+Y' in F1, which include F*S and B11*S*S (B1=[B11*S; F]). I think its might be a typo, because then we have F*S^2 etc. which ist not really LMI. So i ignored this " *s ".

If this is really the reason that im getting different results, how should i solve it?

[Johan Löfberg]

and you are solving a feasibility problem, so no reason to expect that you will get the same solution as reported elsewhere

fredag 11 juni 2021 kl. 12:28:51 UTC+2 skrev yu.ti...@gmail.com:

That was really big help from you, thanks a lot. I think i'm getting there. The coding are running and i'm getting results, but different than the result from the paper(where the whole LMIs come from.)

In original LMI, there is a term -B1*S+Y' in F1, which include F*S and B11*S*S (B1=[B11*S; F]). I think its might be a typo, because then we have F*S^2 etc. which ist not really LMI. So i ignored this " *s ".

If this is really the reason that im getting different results, how should i solve it?

在2021年6月11日星期五 UTC+2 上午11:13:28<Johan Löfberg> 写道：

[Piotr Balik]

I have noticed that running the following reduced code yields no error:

clear,clc
option = sdpsettings('verbose',1,'debug',1);

%some defition about the system
Ap=[-1 0.3; 0 0.5];       Bpu=[0.5;-10];
Bpw=[1;0];                Cpy=[0.4  0.8];
Cpz=[1 0;0 1;0 0];     Dpz=[0;0;1];
Ac=[-4 1;0 -8];         Bc=[-1;-0.5];
Cc=[0.3 -2];             Dc=[1];

H=4.0928;    n=2;
m=1;      p=2;      delta=4;

%switch
% E=sdpvar(n,2,'full');
E=ones(n,2);

Y=sdpvar(2*m,n+m+p,'full');
gama=sdpvar(1,1);
W=sdpvar(n+m+p,n+m+p);
S=diag(sdpvar(2,1,'full'));

A=[Ap Bpu zeros(n);
    H*Dc*Cpy  -H  H*Cc;
    Bc*Cpy  zeros(2,1) Ac];
B11=[Bpu  zeros(2,1);
    zeros(1)    eye(1)];

Fs=E*S;
B1=[B11*S;Fs];

B2=[Bpw;zeros(3,1)];
C2=[Cpz  Dpz  zeros(3,2)];
D2=[-Dpz  zeros(3,1)];

F=[W*A'+A*W    -B1+Y'   B2  W*C2';

    -B1'+Y    -2*S    zeros(2,1) S*D2';
    B2'     zeros(1,2)     -eye(1) zeros(1,3);
    C2*W    D2*S        zeros(3,1) -gama*eye(3,3)]<=0;

sol = optimize(F, [], option)

 yalmipversion: '20210331'
    matlabversion: '9.2.0.538062 (R2017a)'
       yalmiptime: 0.2501
       solvertime: 0.0559
             info: 'Successfully solved (SeDuMi-1.3)'
          problem: 0

However when we "switch" the comment, the variable E somehow degrades the solution to:

Warning: Rank deficient, rank = 29, tol =  3.135768e-11.
Reference to non-existent field 'dinf'.

Where 'dinf' has been mentioned here

But shouldn't the solver come up with feasible tested E=[1 1; 1 1] ?
When debugging, the sedumi actually returns only info.pinf=1, so I assume it jumped out of loop without even creating dinf field
Is that worth investigating?

piątek, 11 czerwca 2021 o 12:59:30 UTC+2 yu.ti...@gmail.com napisał(a):

Alright thanks again for your help!

[Johan Löfberg]

I don't know what you are talking about, because f you activate 

E=sdpvar(n,2,'full');

the code will activate the global solver bmibnb and all bets are off. as you introduce a bilinear expression

Fs=E*S;

despite having been told how to fix this easily

[Johan Löfberg]

however, when running the global solver it appears to work well, testing to minimize gama to see what happens

>> sol = optimize(F,gama)

* Starting YALMIP global branch & bound.

* Upper solver     : fmincon

* Lower solver     : MOSEK

* LP solver        : MOSEK

* -Extracting bounds from model

* -Perfoming root-node bound propagation

* -Calling upper solver (no solution found)

* -Branch-variables : 6

* -More root-node bound-propagation

* -Performing LP-based bound-propagation 

* -And some more root-node bound-propagation

* Warning: 4 branch variables are unbounded from below

* Warning: 6 branch variables are unbounded from above

* Starting the b&b process

 Node       Upper       Gap(%)       Lower     Open   Time

    1 :   9.33099E-01     0.00    9.33099E-01    2     1s  Solution found by heuristics | Added 5 cuts on SDP cone  

* Finished.  Cost: 0.9331 (lower bound: 0.9331, relative gap 5.173e-09%)

* Termination with relative gap satisfied 

* Timing: 75% spent in upper solver (2 problems solved)

*         6% spent in lower solver (13 problems solved)

*         7% spent in LP-based domain reduction (64 problems solved)

*         1% spent in upper heuristics (3 candidates tried)

sol = 

  struct with fields:

    yalmipversion: '20210331'

    matlabversion: '9.9.0.1524771 (R2020b) Update 2'

       yalmiptime: 0.5978

       solvertime: 1.2402

             info: 'Successfully solved (BMIBNB)'

          problem: 0

very bad way to solve a model that trivially can be cast as a linear SDP though

lördag 12 juni 2021 kl. 00:18:40 UTC+2 skrev Piotr Balik:

[Piotr Balik]

Now I see that, what I did is a bit embarrassing. Changing the variables does wonders

Best regards