Dear Prof. Johan Löfberg,

[Hailong Chen]

I encountered a strange situation.

I'm trying to use the bisection method to solve the SOS program.

I declare three SOS multipliers to construct an SOS program. Once it is solved (the program is feasible), I fix one of the three multipliers and then solve the SOS program again. 

The "output log" of the solver shows that the program is infeasible.

It is strange, right?

I don't know why this happened.

Best regards,

Hailong Chen

[Johan Löfberg]

not necessarily as it can be a tolerance issue, i.e. the solution obtained is slightly infeasible, and then when you fix and try to solve again it is impossible to find a feasible solution

a*x^2+b*y^2 should be sos, and due to numerical tolerances it returns a=1, b = -1e-16, and then you fix b to -1e-16 and try to solve, and solver immediately claims infeasibility

[Hailong Chen]

Hi Prof. Johan Löfberg,

Thank you for your explanation!

I agree with you that it appears to be caused by the tolerance issue. However, when I increase the tolerance, the "infeasible" status doesn't change.

Old setting:

options = sdpsettings('solver','mosek',...

'verbose',1,...%'sos.model',2,...

'mosek.MSK_DPAR_INTPNT_CO_TOL_PFEAS',1e-12,...

'mosek.MSK_DPAR_INTPNT_CO_TOL_DFEAS',1e-12,...

'mosek.MSK_DPAR_INTPNT_CO_TOL_REL_GAP',1e-12);

New setting:

options = sdpsettings('solver','mosek',...

'verbose',1,...%'sos.model',2,...

'mosek.MSK_DPAR_INTPNT_CO_TOL_PFEAS',1e-6,...

'mosek.MSK_DPAR_INTPNT_CO_TOL_DFEAS',1e-6...

'mosek.MSK_DPAR_INTPNT_CO_TOL_REL_GAP',1e-3);

As for the example, I think that the situation doesn't fit into it since I fix the variable as one of the multipliers returned from the feasible SOS program.

Best regards,

Hailong

[Johan Löfberg]

you would have to supply minimal reproducible examples (but for the example I gave, you would have to use tighter precision in first call, and sloppier in second. changing in one does not make any difference)

[Hailong Chen]

Hi Prof. Johan Löfberg,

The following code is an example that I deal with. I'm not sure if I've written the code clearly enough. If not, please let me know.

clc

clear all

clear classes

z1 = sdpvar;

z2 = sdpvar;

z3 = sdpvar;

z4 = sdpvar;

z5 = sdpvar;

z6 = sdpvar;

z = [z1;z2;z3;z4;z5;z6];

dlambda1 =2;

ds1 = 2;

ds2 = 2;

f_z =[z3-z5-z2*z3+z2*z5;

z1*z3-z1*z5;

2.52274979573e-07-63.4491948576*z1+4.85059010021*z2-0.377*z3-60.7556781551*z4-3.03222653963*z6-24.7152*z4^2-2.6522070632*z4*z6-55.4973055463*z1*z4+4.24268080116*z2*z4-39.6636962907*z1*z6+3.03222653963*z2*z6-34.6927691846*z1*z4*z6+2.6522070632*z2*z4*z6;

-0.0291063949887*z1-0.380732506579*z2-0.6634*z4+0.238005518318*z6-0.0181951435876*z1*z6-0.238005518318*z2*z6;

-5.15416831064e-05+77.593075149*z1+45.996277102*z2-40.2317074323*z4-0.5344*z5-78.8848265475*z6-25.1498577473*z4*z6+67.8685775208*z1*z4+40.2317074323*z2*z4+48.5053935496*z1*z6+28.7534360297*z2*z6+42.4263641565*z1*z4*z6+25.1498577473*z2*z4*z6-15.6692*z6^2;

0.206670513732*z1-0.348641275196*z2+0.304947153688*z4-0.7344*z6+0.180769144097*z1*z4-0.304947153688*z2*z4];

Gz = z1^2+z2^2-2*z2;

Num_epsilon = 1e-4;

ell2 = Num_epsilon*(z1^4+z2^4+z3^4+z5^4+z4^4+z6^4);

Vz = -0.0739874424389*z2*z3-0.129043912198*z2*z5+0.0199135250803*z1*z3-0.0431598064656*z1*z5+23.1007622703*z4^2+8.93486277337*z4*z6+3.37904960982*z1*z4+13.1391048185*z2*z4-12.1562240419*z1*z6+9.84207067669*z2*z6+23.0687967514*z6^2+11.4659775908*z1^2+10.2372018995*z2^2+2.72187354354*z1*z2+0.0827356036976*z3^2-0.0244049688384*z3*z4-0.16092294445*z3*z5+0.04842907423*z4*z5+0.0821052573216*z5^2+0.186458174236*z3*z6-0.0823722997125*z5*z6;

c = 0.0062;

dVz = jacobian(Vz,z);

[s1,s1_coef] = polynomial(z,ds1,0);

[s2,s2_coef] = polynomial(z,ds2,0);

[lambda1,lambda1_coef] = polynomial(z,dlambda1,0);

SOS1 = -s1*(c-Vz)-s2*dVz*f_z-lambda1*Gz-ell2;

constr = [sos(s1),sos(s2),sos(SOS1)];

options = sdpsettings('solver','mosek',...

'verbose',1,...%'sos.model',2,...

'mosek.MSK_DPAR_INTPNT_CO_TOL_PFEAS',1e-6,...

'mosek.MSK_DPAR_INTPNT_CO_TOL_DFEAS',1e-6,...

'mosek.MSK_DPAR_INTPNT_CO_TOL_REL_GAP',1e-3);

sol1 = solvesos(constr,0,options,[lambda1_coef;s1_coef;s2_coef]);

s2_j = replace(s2,s2_coef,value(s2_coef));

s2 = s2_j;

dVz = jacobian(Vz,z);

[s1,s1_coef] = polynomial(z,ds1,0);

[lambda1,lambda1_coef] = polynomial(z,dlambda1,0);

SOS1 = -s1*(c-Vz)-s2*dVz*f_z-lambda1*Gz-ell2;

constr = [sos(s1),sos(SOS1)];

sol2 = solvesos(constr,0,options,[lambda1_coef;s1_coef]);

You can see from the code that I try to fix $s_2$ and then solve the program again. However, it is infeasible.

Thank you in advance. If anything needs to be clarified, please let me know.

Best regards,

Hailong

[Johan Löfberg]

You are suffering from these classical numerical issues

s2 is not actually sos

solvesos(sos(replace(s2,s2_coef,value(s2_coef))))

We try to make a small perturbation to make it feasible (according to solver with tolerances still)

t = sdpvar(length(s2_coef),1);

solvesos(sos(replace(s2,s2_coef,t + value(s2_coef))),norm(t,1),[],t)

and the largest perturbation it needed is 10^-9 or so

and that solution happens to be truly sos (according to solver, once again remembering that the solver has tolerances)

solvesos(sos(replace(s2,s2_coef,value(t) + value(s2_coef))))

[Hailong Chen]

Hi Prof. Johan Löfberg,

I understand why this happened this time. Thank you very much! I will pay attention to the tolerance issue.

Best regards,

Hailong