dual(constraints(1)) returns NaN array

[Danylo Malyuta]

Here is the complete code for an optimization problem that I'm trying to solve, which finds the optimal force input to move a mass from one point to another:

% Attempt to solve a simple problem of finding the optimal force input to minimize

% transport time of a mass, using the optimizer() function. This is an SCP (Sequential

% Convex Programming) problem.

S = 10; % [m] Path length

delta = 0.1; % [m] Spatial discretization

s = 0:delta:S; % [m] Discretized distances

N = length(s)-1; % Number of discrete points before the final point

Qf = 1000; % [J] Initial energy in battery

% *** YALMIP variables

Ekin = sdpvar(N+1,1);

SOC  = sdpvar(N+1,1);

F    = sdpvar(N,1);

v    = sdpvar(N+1,1);

dtds = sdpvar(N+1,1);

% *** Constraints

constraints = [...

Ekin(2:N+1) == Ekin(1:N)+delta*F

SOC(2:N+1)  == SOC(1:N)-delta/Qf*F

F(:)        <= 10

F(:)        >= -10

Ekin(:)     >= 0.5*v(:).^2

Ekin(1)     == 0.01

Ekin(N+1)   == 0.01

SOC(1)      == 1

cone([v.'+dtds.'

      2*ones(1,N+1)

      v.'-dtds.'])

];

% *** Objective function

cost = sum(dtds);

% *** YALMIP options

yalmip_opts = sdpsettings('verbose',1,'solver','gurobi');%'ecos','ecos.maxit',200);

% *** Solution

sol = optimize(constraints,cost,yalmip_opts);

% *** Extract duals of the first set of equality constraints (kinetic energy dynamics)

Ekin_duals = dual(constraints(1)); % PROBLEM: returns vectors of 100 NaN elements ?????

On the last line, I try to get the duals of the first 100 equality constraints. Using both GUROBI and ECOS as solvers, the command returns a [100x1] array of NaN. What am I doing wrong, how do I get the duals? Thanks a lot!

[Danylo Malyuta]

Correction: change the first three lines to:

% Attempt to solve a simple problem of finding the optimal force input to minimize

% transport time of a mass, using the optimize() function.

My bad for not checking them first.

[Johan Löfberg]

I think you cannot extract duals if YALMIP has been forced to do reformulations. Hence, you have to use the cone operator. Additionally, some solvers don't compute duals for SOCP constraints. With gurobi, I think you have to tell it to compute them using the option qcpdual

[Danylo Malyuta]

But I am already using the cone operator, so where else should I use it? I tried using 

yalmip_opts = sdpsettings('verbose',1,'solver','gurobi','gurobi.QCPDual',1);

but it gives NaN as well for the current problem formulation.

[Johan Löfberg]

You still have 

[Danylo Malyuta]

Got it, it works now. Thanks a lot for the super fast reply!

[Rachel]

Dear Johan 

  I have met the similar question:dual(constraints(1)) returns NaN array.

  It's like a unit commitment without variable onoff. The objective is to minimize the generating costs (as all the generators are on, it's a LP model with only one decision variable,P).When I used the expression as objective=Q*P^2+C*P, I could easily get the dual variables. But when I change the expression of objective into one based on block bidding expression:  (NM means the total numbers of bidding segments;Ci,m is the electricity price during the relative segment;Pi,t,m is the winnnig power in relative bidding segment of unit i at time t. The product of Ci,m and Pi,t,m is the the cost of  power in this bidding segment of unit i at time t  )and modify the code,dual(Constraints(1))returns to NaN array(1*Horizon). The codes are shown below:

Horizon=48;

Nunits=3;

...

costC=[

     0 0 0 600 50 700 100;

     0 0 0 500 50 800 100;

    0 0 0 480 50 850 100;

    ]; %P1-----Price1----P2-----Price2-----P3-----Price3------P4,3 bidding segments

HcostC=size(costC,2);

Objective=0;%Total Cost 

for t=1:Horizon

    for i=1:Nunits

        for m=(HcostC-2):-2:1 %every terminal of bidding range 

            c=max((P(i,t)-costC(i,m)),0)-max((P(i,t)-costC(i,m+2)),0);% extract the winning power during thr relative bidding segment

          Objective=Objective+ (max((P(i,t)-costC(i,m)),0)-max((P(i,t)-costC(i,m+2)),0)) *costC(i,m+1)*24/Horizon*1;    %caculate the cost of  power in this bidding segment of unit i at time t    

        end

    end

end

ops = sdpsettings('verbose',1,'debug',1);

optimize(Constraints,Objective,ops);

stairs(value(P)');

legend('Unit 1','Unit 2','Unit 3');

lameda=dual(Constraints(1))

Result:

Optimal solution found (tolerance 1.00e-04)

Best objective 1.330438192796e+06, best bound 1.330438192796e+06, gap 0.0000%

lameda =

  Column 1 - 16 

   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN

   Column17- 32

   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN

   Column 33- 48 

   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN   NaN

Thanks 

Rachel

[Johan Löfberg]

Your objective appears pwa nonconvex and thus MILP-represented

[Rachel]

Thanks for reply! So how can I change the code to get the Lagrange multiplier of the Constraint(1) (equation)? I really have no idea.

[Johan Löfberg]

you cannot, it is a MILP

for simple milps, you can use this interpretation

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

but it won't work for you as you have no control of the internal binaries used to model the nonconvex max operator use. to do it, you would have to manually model that operator so you have the binary available

[Rachel]

Thank you so much! I use another method to calculate the cost based on block bidding, and only use one max operator. Now it works. Thanks for your advice!

[Md Habib Ullah]

Hi Johan and All,

I have the following optimization code with a quadratic constraint. I am using 'gurobi.QCPDual' to get the dual variable for the quadratic constraint. It gives me results when I solve it without 'gurobi.QCPDual', although it returns NaN for lambda. However, I get NaN for the whole solution when I use 'gurobi.QCPDual'.

Could you please help me to solve this problem. Your response is highly appreciated. Thanks. 

clear all

close all

clc

a = [0.002; 0.0035; 0.003; 0.0045];

b = [3; 6; 5; 9];

pmin = [0; 0; 0; -500];

pmax = [200; 300; 200; 0];

mu = [1; 10; 20; -10];

sdpvar p(4,1,'full');

%%

%% Objective Function

%%

Obj = sum( (a.*(p.^2)) + (b.*p) ) ;

%

%% Constraints

C_eq = [ sum(p) >= 0 ];

C_ineq = [ pmin <= p <= pmax ];

c_tst = [p(1)^2 + p(2)^2 <= 400000];

%

F = [C_eq, C_ineq, c_tst];

% %%

ops = sdpsettings('verbose',1,'solver','gurobi','gurobi.QCPDual',1);

% ops = sdpsettings('verbose',1,'solver','gurobi');

% ops = sdpsettings('solver','gurobi');

optimize(F, Obj, ops);

Object = value(Obj)

p = value(p)

lambda = dual(c_tst)

[Johan Löfberg]

Confirmed. Fixed in the develop branch

(simply change line 548 in solvesdp to if any(strfind(solver.tag,'GUROBI')))

[Md Habib Ullah]

Hi Johan. Thanks for your comment, but I'm not sure what you meant. I don't have line 548 in my code. 

[Johan Löfberg]

solvesdp.m definitely has 548 lines.

If you are unsure about editing YALMIP, simply install the latest develop version which has been fixed

[Md Habib Ullah]

Thanks a lot, Johan. Got it, it works now. I install the latest develop version.