Yalmip not working with Gurobi

[jialin liu]

I was using matlab solver (sdpsettings('solver','linprog')) to solve a problem. Now I added gurobi to my working path. 

>> which gurobi

/Users/AndersonLab/Dropbox/Research/P3/Code/gurobi604/mac64/matlab/gurobi.mexmaci64

I also changed the setting to sdpsettings('solver','gurobi'). I did not do any changes except for sdpsettings in terms of code. It did not get faster to solve the problem. Why? Did I miss anything?

Thanks!

[Johan Löfberg]

Perhaps the problem is so small and trivial that it doesn't make a difference which solver you are usng

[jialin liu]

Hi Professor Löfberg,

Thanks for the reply.  Sorry, I did not explain my problem clearly. I think Yalmip is still using the default solver instead of Gurobi, because the time to solve and results are the same as the default solver's when I use Gurobi. I do not have Cplex but I tried sdpsettings('solver','cplex'). It worked and gave me the same thing. 

Also the problem is not really small, it takes a few minutes to solve. Right now, the quality of the solution is really bad. What could I do?

Thanks,
Jialin

[Johan Löfberg]

Gurobi would always be selected by YALMIP if it was available. Hence you have two errors

1. Gurobi is not on the matlab path

2. You are not using the options structure, you only create it

[Johan Löfberg]

and show me the output from

sdpvar x
optimize(x>=0,x)

and

sdpvar x
optimize(x>=0,x,sdpsettings('solver','gurobi'))

[jialin liu]

The results from running your suggested code are :

sdpvar x
optimize(x>=0,x)

Optimize a model with 1 rows, 1 columns and 1 nonzeros

Coefficient statistics:

  Matrix range    [1e+00, 1e+00]

  Objective range [1e+00, 1e+00]

  Bounds range    [0e+00, 0e+00]

  RHS range       [0e+00, 0e+00]

Presolve removed 1 rows and 1 columns

Presolve time: 0.00s

Presolve: All rows and columns removed

Iteration    Objective       Primal Inf.    Dual Inf.      Time

       0    0.0000000e+00   0.000000e+00   0.000000e+00      0s

Solved in 0 iterations and 0.00 seconds

Optimal objective  0.000000000e+00

ans = 

    yalmiptime: 1.0762

    solvertime: 0.0071

          info: 'Successfully solved (GUROBI-GUROBI)'

       problem: 0

sdpvar x
optimize(x>=0,x,sdpsettings('solver','gurobi'))

Optimize a model with 1 rows, 1 columns and 1 nonzeros

Coefficient statistics:

  Matrix range    [1e+00, 1e+00]

  Objective range [1e+00, 1e+00]

  Bounds range    [0e+00, 0e+00]

  RHS range       [0e+00, 0e+00]

Presolve removed 1 rows and 1 columns

Presolve time: 0.00s

Presolve: All rows and columns removed

Iteration    Objective       Primal Inf.    Dual Inf.      Time

       0    0.0000000e+00   0.000000e+00   0.000000e+00      0s

Solved in 0 iterations and 0.00 seconds

Optimal objective  0.000000000e+00

ans = 

    yalmiptime: 0.1146

    solvertime: 0.0035

          info: 'Successfully solved (GUROBI-GUROBI)'

       problem: 0

The way I solve my problem is :

%% Solve Problem

    sdpsettings('solver','gurobi');

%     opts = sdpsettings('verbose',0,'saveduals',0);

    solopt=optimize(constraints,Objective);

It did not indicate Gurobi in the results. My results are like:

Optimize a model with 10488 rows, 5376 columns and 57192 nonzeros

Model has 1800 quadratic objective terms

Coefficient statistics:

  Matrix range    [3e-01, 8e+01]

  Objective range [1e+00, 7e+00]

  Bounds range    [1e+00, 1e+00]

  RHS range       [2e+00, 2e+02]

Presolve removed 8016 rows and 3384 columns

Presolve time: 0.09s

Presolved: 2472 rows, 1992 columns, 12096 nonzeros

Presolved model has 1800 quadratic objective terms

Variable types: 1848 continuous, 144 integer (144 binary)

Root relaxation: objective 6.269299e+04, 31841 iterations, 2.37 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work

 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 62692.9949    0  101          - 62692.9949      -     -    2s

H    0     0                    62737.175161 62692.9949  0.07%     -    2s

     0     0 62695.5999    0   93 62737.1752 62695.5999  0.07%     -    2s

     0     0 62698.8877    0  105 62737.1752 62698.8877  0.06%     -    3s

     0     0 62699.9401    0  108 62737.1752 62699.9401  0.06%     -    3s

     0     0 62701.8630    0  109 62737.1752 62701.8630  0.06%     -    3s

     0     0 62703.7732    0  109 62737.1752 62703.7732  0.05%     -    3s

     0     0 62705.0183    0  108 62737.1752 62705.0183  0.05%     -    4s

     0     0 62705.2829    0  106 62737.1752 62705.2829  0.05%     -    4s

     0     0 62705.5523    0  103 62737.1752 62705.5523  0.05%     -    4s

     0     0 62705.8342    0  102 62737.1752 62705.8342  0.05%     -    4s

     0     0 62705.8342    0  103 62737.1752 62705.8342  0.05%     -    4s

     0     0 62705.8342    0  102 62737.1752 62705.8342  0.05%     -    5s

     0     0 62705.8342    0   98 62737.1752 62705.8342  0.05%     -    5s

     0     0 62705.8342    0   98 62737.1752 62705.8342  0.05%     -    5s

     0     2 62705.8342    0   98 62737.1752 62705.8342  0.05%     -    6s

H  131   109                    62736.180522 62706.1212  0.05%  86.9    7s

H  140   119                    62733.180630 62706.1212  0.04%  86.6    7s

H  181   141                    62733.179561 62706.1212  0.04%  82.4    8s

H  203   156                    62733.179421 62706.1212  0.04%  80.4    8s

H  214   156                    62732.179455 62706.1212  0.04%  78.7    8s

H  218   155                    62731.179166 62706.1212  0.04%  77.7    8s

H  253   168                    62729.175205 62706.1212  0.04%  73.3    8s

Interrupt request received

Cutting planes:

  Gomory: 18

  Cover: 42

  Implied bound: 47

  MIR: 2

  Flow cover: 85

  Flow path: 3

Explored 447 nodes (73695 simplex iterations) in 9.47 seconds

Thread count was 2 (of 4 available processors)

Solve interrupted

Best objective 6.272917520536e+04, best bound 6.270666260368e+04, gap 0.0359%

[Johan Löfberg]

Gurobi is used on your machine.

Your claims are not true though. You have never used linprog to solve this problem. linprog is for LPs only, and from the display, it is obvious that you are solving a MILP.

[jialin liu]

For your reference, I attached my complete code here. Thanks!

clc;

clear all;

close all;

    

%% System Parameters

    nt=24; % time horizon

    % define cost

    cg1=[2 1.75 1 3.25 3 3]'; % 1st order generator cost

    cg2=[0.02 0.0175 0.0625 0.00834 0.025 0.025]'; % 2nd order generator cost

    crgs=[6 6.75 7 5.25 5 5]';

    Conoff=ones(1,6); % generator on cost

    

    mpc=case30;

    loads=repmat(mpc.bus(:,3),1,nt);

    mg=20*ones(1,nt); 

    mbus=10; 

    wbus=22; 

    wl=0; 

    

%% Create B Matrix

    % make Y/B bus matrix and A matrix in Pline=A*Pinj

    [Ybus, Yf, Yt] = makeYbus(mpc);

    B=full(real(i.*Ybus));

    NB=-B;

    Bred=B(1:29,1:29); % Reduced B Matrix

    

    from = mpc.branch(:,1);

    to = mpc.branch(:,2);

    M = zeros(41,30); % from/to matrix

    lineC = zeros(41,1);

    for i=1:41

        M(i,from(i))=1;

        M(i,to(i))=-1;

        lineC(i)=1*mpc.branch(i,7);

    end

    m = M(:,1:29);

    

    x=mpc.branch(:,4); % resistance

    r=mpc.branch(:,3); % reactance

    b=x./(x.^2+r.^2); % susceptance

    D=diag(b); % diagnol matrix of the susceptance

    

%% Wind Forecast and Scenarios

    %wind forecast

    load('ObsDays.mat');

    k=54;

    wf1=wl*ObsDays(k,:,11);

    wf1=wf1(1:nt);

    wf=wf1;

    corr1=corrcoef(ObsDays(:,:,11));

    std_dev = [0.12,0.15,0.18,0.5,0.6,0.67,0.72,0.76,0.79,0.82,0.83,0.8315,0.833,0.835,0.836,0.838,0.839,0.841,0.842,0.844,0.845,0.847,0.848,0.85];

    for j=1:nt

        for k=1:nt

        covarr1(j,k)=corr1(j,k)*std_dev(k)*std_dev(j);

        end

    end

    beta=0.01;

    eulernum=exp(1);

    epsilon=0.05;

    Ndelta=1*nt;%num of uncertainties,one wind each period

    Nneed=ceil((1/epsilon)*(eulernum/(eulernum-1))*(log(1/beta)+4*(Ndelta+1)-1));

    N=Nneed;%Number of scenarios

    %create wind error set

    avg=zeros(nt,1)';

    for i=1:N

        err1(i,:)=mvnrnd(avg,covarr1);

    end

    for j=1:N

        wn1(j,:)=wf1+wf1.*err1(j,:);

        pos=wn1(j,:)>=0;

        winderror1(j,:)=wf1.*err1(j,:).*pos;

        ws1(j,:)=wn1(j,:).*pos;

    end

    %total wind forecast

    wst=ws1;

    %create wind error set

    pr1=90;%upper percentile(throughout)

    pr2=10;%lower percentile(throughout)

    windup1=prctile(winderror1,pr1);

    windup=1*windup1;

    winddown1=prctile(winderror1,pr2);

    for i=1:nt

        if wf1(i)+winddown1(i)<0

            winddown1(i)=-wf1(i);

        end

    end

    winddown=1*winddown1;

    

    phi=wst(1,:)-wf;% forecast error

%% Define Optimization Variables and Their Bounds

    % 3. Generators

    ng=6;

    Gmax=[80*ones(1,nt);80*ones(1,nt);40*ones(1,nt);50*ones(1,nt);30*ones(1,nt);55*ones(1,nt)];

    Gmin=[zeros(6,nt)];

    pg=sdpvar(ng,nt,'full');

    onoff=binvar(ng,nt,'full');

    Rgsmin=-0.5*Gmax;

    Rgsmax=0.5*Gmax;

    rgs=sdpvar(ng,nt,'full');

    % LB & UB 

    constraints = [Gmin<=pg<=Gmax,Rgsmin<=rgs<=Rgsmax];   

%% Power Flow Constraints

    % power flow with forecast wind

    Pinj=sdpvar(30,nt,'full'); % bus nodal matrix with forecast wind

    for i=1:30

      if i==1

      constraints=[constraints,Pinj(i,:)==-loads(i,:)+pg(1,:)];

      elseif i==2

      constraints=[constraints,Pinj(i,:)==-loads(i,:)+pg(2,:)];

      elseif i==13

      constraints=[constraints,Pinj(i,:)==-loads(i,:)+pg(3,:)];

      elseif i==22

      constraints=[constraints,Pinj(i,:)==-loads(i,:)+pg(4,:)];

      elseif i==23

      constraints=[constraints,Pinj(i,:)==-loads(i,:)+pg(5,:)];

      elseif i==27

      constraints=[constraints,Pinj(i,:)==-loads(i,:)+pg(6,:)];

      elseif i==wbus

      constraints=[constraints,Pinj(i,:)==-loads(i,:)+wf];

      elseif i==mbus

      constraints=[constraints,Pinj(i,:)==-loads(i,:)+mg];

      else

      constraints=[constraints,Pinj(i,:)==-loads(i,:)];

      end

    end

   

   theta=sdpvar(29,nt,'full'); % bus angle

   constraints = [constraints,theta==inv(Bred)*Pinj(1:29,:)];

   pflow=sdpvar(41,nt,'full'); % line flow

   constraints = [constraints,pflow==D*m*theta];

   constraints = [constraints,-repmat(lineC,1,nt)<=pflow<=repmat(lineC,1,nt)];%  line constraints

   

   % power flow with wind scenario

    Pinjs=sdpvar(30,nt,'full'); % bus nodal matrix with real wind           

    for i=1:30

          if i==1

          constraints=[constraints,Pinjs(i,:)==-loads(i,:)+pg(1,:)+rgs(1,:)];

          elseif i==2

          constraints=[constraints,Pinjs(i,:)==-loads(i,:)+pg(2,:)+rgs(2,:)];

          elseif i==13

          constraints=[constraints,Pinjs(i,:)==-loads(i,:)+pg(3,:)+rgs(3,:)];

          elseif i==22

          constraints=[constraints,Pinjs(i,:)==-loads(i,:)+pg(4,:)+rgs(4,:)];

          elseif i==23

          constraints=[constraints,Pinjs(i,:)==-loads(i,:)+pg(5,:)+rgs(5,:)];

          elseif i==27

          constraints=[constraints,Pinjs(i,:)==-loads(i,:)+pg(6,:)+rgs(6,:)];

          elseif i==wbus

          constraints=[constraints,Pinjs(i,:)==-loads(i,:)+wst(1,:)];

          elseif i==mbus

          constraints=[constraints,Pinjs(i,:)==-loads(i,:)+mg];

          else

          constraints=[constraints,Pinjs(i,:)==-loads(i,:)];

          end

    end

    

   thetas=sdpvar(29,nt,'full'); % bus angle

   constraints = [constraints,thetas==inv(Bred)*Pinjs(1:29,:)];

   pflows=sdpvar(41,nt,'full'); % line flow

   constraints = [constraints,pflows==D*m*thetas];

   constraints = [constraints,-repmat(lineC,1,nt)<=pflows<=repmat(lineC,1,nt)]; %  line constraints

   

%% Constraints

  

  %constraints=[constraints,sum(rlsup)+sum(rgsdn)==pvup,sum(rlsdn)+sum(rgsup)==pvdn]; 

  constraints=[constraints,-sum(rgs)==max(0,phi)-max(0,-phi)]; % secondary reserve constraint

  

% Generator Constraints

  Rdn=0.3*Gmax;

  Rup=0.3*Gmax;

  constraints=[constraints,-Rdn(:,2:nt)<=pg(:,2:nt)-pg(:,1:nt-1)<=Rup(:,2:nt)]; % ramping constraints

  constraints=[constraints,Gmin.*onoff<=pg+rgs<=Gmax.*onoff]; % Generator Bounds with Rgs

  constraints=[constraints,-Rdn(:,2:nt)<=(pg(:,2:nt)+rgs(:,2:nt))-(pg(:,1:nt-1)+rgs(:,1:nt-1))<=Rup(:,2:nt)]; % Ramping Constraints with Rgs

  constraints=[constraints,sum(pg)+wf-sum(loads)+mg==0];  % Power balance equation forecast

  

%% Objective

     Objective=sum(pg')*cg1+sum(pg')*diag(cg2)*sum(pg')'+sum(abs(rgs'))*crgs+sum(onoff')*Conoff';

%% Solve Problem

    sdpsettings('solver','gurobi');

%     opts = sdpsettings('verbose',0,'saveduals',0);

    solopt=optimize(constraints,Objective);

%% Optimization Variables

    pg=value(pg);

    rgs=value(rgs);

    onoff=value(onoff);

    Pinj=value(Pinj);

    Pinjs=value(Pinjs);

    pflow=value(pflow);

    pflows=value(pflows);

[jialin liu]

Thanks. Prof. Löfberg. I will try to install CPLEX to see if I get better quality solution.

[Johan Löfberg]

You have a solution which is within 0.04% of optimal. Your problem is not the solver, but the model, if you get a solution which you don't like.

[jialin liu]

Thanks! Prof. Lofberg. I do find a problem in my model. I really appreciate your help along the way!