some questions about "implies" operator

[刘浩翔]

dear professor,

I wanted to model this "if" struct below using "implies" operator:

finally I got NaN result in my code .However,after removing the "implies" constraints I can get a reasonable result . So I wonder what is wrong about my code below.thank you so much!

 

clear

s=100;

N=100;

ds=s/N;

v=sdpvar(1,N);

Te=sdpvar(1,N-1);

Tb=sdpvar(1,N-1);

soc=sdpvar(1,N);

d=binvar(2,1);

socdot=sdpvar(1,N-1);

 

opts=sdpsettings;

opts.solver = '';

opts.verbose = 0;

 

cons1=[];

cons2=[];

cons3=[];

cons4=[];

cons5=[];

cons6=[];

cons9=[];

obj1=[];

 

for i=1:N-1

vdot=[(Te(i)*19-Tb(i)-0.5*4000*9.8*0.02-0.5*0.5*2.5*0.3*1.29*v(i)^2)/(4000*0.5)];

cons1=[cons1,v(i+1)*v(i)==v(i)^2+ds*vdot];

cons2=[cons2,-500<=Te(i),Te(i)<=500];

cons3=[cons3,0<=Tb(i),Tb(i)<=0.5*4000*9.8*0.6];

cons4=[cons4,0<=v(i),v(i)<=10];

cons9=[cons9,sum(d) == 1,

implies(d(1), [Te(i) <= 0, socdot(i) == -(400-(160000-1444*v(i)*Te(i))^(0.5))/6200*v(i)]);

implies(d(2), [Te(i) >= 0, socdot(i) == -(400-(160000-1600*v(i)*Te(i))^(0.5))/6200*v(i)])];

% socdot= [-(400-(160000-1600*v(i)*Te(i))^(0.5))/6200*v(i)];

cons5=[cons5,soc(i+1)==soc(i)+socdot(i)*ds];

cons6=[cons6,0<=soc(i),soc(i)<=1];

obj1=[obj1,ds*cpower(v(i),-1)];

end

cons7=[v(1)==1,v(N)==1];

cons8=[soc(1)==1];

cons=cons1+cons2+cons3+cons4+cons5+cons6+cons7+cons8+cons9;

 

obj=sum(obj1);

sol = optimize(cons,obj,opts);

 

[Johan Löfberg]

You will get an absolutely awful mixed-integer nonlinear program, and there are no well-working solvers for that

You can reformulate it to a mixed-integer quadratic problem (and thus use gurobi) as you can new variable s and z, with constraints

socdot*V == s

(6200*s + 400)^2 == 160000-z

implies(d1, [Te >= 0, z == 1600*V*Te])

implies(d2, [Te <= 0, z == 1444*V*Te])

BTW, your model does not make sense logically as your switches d1 and d2 are fixed over the horizon

[刘浩翔]

Thank you professor! I have removed the above equations from my problem.
BTW, I wonder whether I can get the dual variables of the "cons1" above? I want the lagrange multiplier of the "v" term.  

[Johan Löfberg]

There is no such thing as duals for a mixed-integer program

[刘浩翔]

I mean the problem without that "if" struct ,which doesn't include the mixed-integer program. 

When I run "dual(cons1(i))" ,it gets NaN. Can I get the dual variables of the "cons1" or how can I correct my code? Code is below.

Thank you so much

clear

s=100;

N=100;

ds=s/N;

v=sdpvar(1,N);

Te=sdpvar(1,N-1);

Tb=sdpvar(1,N-1);

opts=sdpsettings('solver','gurobi','verbose',0);

cons1=[];

cons2=[];

cons3=[];

cons4=[];

obj1=[];

for i=1:N-1

vdot=[(Te(i)*19-Tb(i)-0.5*4000*9.8*0.02-0.5*0.5*2.5*0.3*1.29*v(i)^2)/(4000*0.5)];

cons1=[cons1,v(i+1)*v(i)==v(i)^2+ds*vdot];

cons2=[cons2,-500<=Te(i),Te(i)<=500];

cons3=[cons3,0<=Tb(i),Tb(i)<=0.5*4000*9.8*0.6];

cons4=[cons4,0<=v(i),v(i)<=10];

obj1=[obj1,ds*cpower(v(i),-1)];

end

cons5=[v(1)==1,v(N)==1];

cons=cons1+cons2+cons3+cons4+cons5;

obj=sum(obj1);

sol = optimize(cons,obj,opts);

[Johan Löfberg]

duals are not extracted for models solved by nonlinear solvers

On Tuesday 9 January 2024 at 02:23:19 UTC+1 cokewit...@gmail.com wrote:

[刘浩翔]

I got it, but why can't I extract the duals when solving some linear models?

[Johan Löfberg]

Linear models you can always extract the dual. Hence have something in the model which makes it outside a linear program although you think it is linear

[刘浩翔]

Thank you so much ,professor

So what is wrong about my model below, I think it is linear but dual(cons1(i)) returns NaN

clear

s=100;

N=100;

ds=s/N;

v=sdpvar(1,N);

Te=sdpvar(1,N-1);

Tb=sdpvar(1,N-1);

opts=sdpsettings('solver','gurobi','verbose',0,'savesolveroutput',1);

cons1=[];

cons2=[];

cons3=[];

cons4=[];

obj1=[];

for i=1:N-1

vdot=[(Te(i)*19-Tb(i)-0.5*4000*9.8*0.02-0.5*0.5*2.5*0.3*1.29*v(i))/(4000*0.5)];

cons1=[cons1,v(i+1)==v(i)+ds*vdot];

cons2=[cons2,-500<=Te(i),Te(i)<=500];

cons3=[cons3,0<=Tb(i),Tb(i)<=0.5*4000*9.8*0.6];

cons4=[cons4,0<=v(i),v(i)<=10];

obj1=[obj1,ds*cpower(v(i),-1)];

end

cons5=[v(1)==1,v(N)==1];

cons=cons1+cons2+cons3+cons4+cons5;

obj=sum(obj1);

sol = optimize(cons,obj,opts);

[Johan Löfberg]

it uses cpower, hence it requires second-order cone programming

[刘浩翔]

Thank you professor

Is there any other good way to model "1/v(i)" in the objective? besides the SOCP and direct way

[Johan Löfberg]

No

But if you want the duals you can extract these for SOCPs, if you do the modelling of 1/v using the cone operator manually

[刘浩翔]

Thanks 

I ran the code below ,exactly as the example shown in  https://yalmip.github.io/faq/extractdual/, but still failed to get the duals ,why?

clear

x = sdpvar(2,1);

Model1 = [x'*x <= 1];

options1 = sdpsettings('solver','gurobi','gurobi.qcpdual',1);

sol1=optimize(Model1,sum(x),options1);

dual(Model1(1))

Model2 = [cone([1;x])];

options2 = sdpsettings('solver','gurobi','gurobi.qcpdual',1);

optimize(Model2,sum(x),options2);

dual(Model2(1))

[Johan Löfberg]

Assuming the problem actually was solved (which we don't know, and you don't seem to check that) the second case works as expected. For models where reformulation is required (first) duals are not supported

Iter       Primal          Dual         Primal    Dual     Compl     Time
   0   0.00000000e+00  0.00000000e+00  0.00e+00 3.16e+00  9.26e-01     0s
   1  -5.73848504e-01 -4.40389050e-01  0.00e+00 1.05e+00  2.17e-01     0s
   2  -1.41284816e+00 -1.49452246e+00  0.00e+00 1.16e-06  2.04e-02     0s
   3  -1.41421032e+00 -1.41429588e+00  0.00e+00 3.17e-10  2.14e-05     0s
   4  -1.41421356e+00 -1.41421365e+00  0.00e+00 4.94e-10  2.42e-08     0s

Barrier solved model in 4 iterations and 0.04 seconds (0.00 work units)
Optimal objective -1.41421356e+00

Solving KKT system to obtain QCP duals...


ans =

    1.4142
    1.0000
    1.0000

[刘浩翔]

sorry I don't really get your point

But I found that ,with Gurobi, if I add this condition " 'gurobi.qcpdual,' 1" to a problem which could have been solved normally, it becomes unsolvable instead(it gave me a NaN result). What may be the problem?

[Johan Löfberg]

It might be that you are using gurobi 11 which messes up the logic for options if you explciitly use these in some scenarios vs 10

Hence, has the problem actually been solved or has gurobi crashed. Since you don't seem to look at the diagnostic, you don't know

The line gurobi.TuneResults = -1; must be removed in sdpsettings.m if you use gurobi 11 and alter options

[刘浩翔]

I checked the info and it says "(Unknown problem (<a href="yalmip.github.io/inside/debug">learn to debug</a>) (Error using gurobiGurobi error 10008: Unable to set parameter TuneTimeLimit to value -1 (minimum is 0)))" 

[Johan Löfberg]

As I guessed then, so easy fix

[刘浩翔]

Thank you so much ,it works now.

I have one more question,how should I model the socp constraint below in yalmip?

cone([dtds+v;

    2*ones(1,N);

    dtds-v])

Is this correct?

[Johan Löfberg]

looks right

[Mark L. Stone]

I am using Gurobi 11 with non-default options with YALMIP develop 20240102.

I commented out the line

gurobi.TuneResults = -1;

from sdpsettings.m

per your instructions above.

I still get the error:

Error using gurobi

Gurobi error 10008: Unable to set parameter TuneTimeLimit to value -1 (minimum is 0)

Error in callgurobi (line 20)
result = gurobi(model,model.params);
Error in solvesdp (line 420)
    eval(['output = ' solver.call '(interfacedata);']);
Error in optimize (line 31)
[varargout{1:nargout}] = solvesdp(varargin{:});

[Johan Löfberg]

*TuneTimeLimit* 

[Mark L. Stone]

Thanks. That worked. It was pretty stupid of me not to see and correct your typo.