BMIBNB, lower solver very quick and accurate. Upper solver very slow and non-converging

[Michiel Jansen]

Hi,

I am making a natural gas model with storage. This is a non-linear non convex problem.

I found out that the BMIBNB gives the most accurate result, however when I increase the size of the system it stops converging.

I can see that the lower solver (CPLEX) really quickly finds an accurate solution, while the upper solver (fmincon) has a lot of trouble converging.

Is there a way to force the upper solver in checking the solution of the lower solver? As this would create an optimization model which would work.

Thanks in advance!

My code is:

% Start model
Pr = sdpvar(nnode,horizon);                                               
Si = sdpvar(nnode,horizon);                                            
shedg=sdpvar(nnode,horizon);                                             
Psto=sdpvar(nnode,horizon);
Pin=sdpvar(nnode,horizon);
Pout=sdpvar(nnode,horizon);

%% Constraints - General
assign(Pr,gas.node(:,9)*ones(1,horizon));                                  % assign initial values
assign(Si,gas.node(:,10)*ones(1,horizon));
assign(shedg,zeros(nnode,horizon));                                        % assign initial values

Constraints = [];                                                          % set constraints to zero

%% Constraints - Branch flow constraint
Constraints = [Constraints, 0.01 <= Gasf*(Pr.^2) <= (Capmax.^2)];   % fix maximum pipe flow

%% Constraints - Equality constraint

Constraints = [Constraints,Si-Pin+Pout+shedg-d-((Cft*(sqrtm(Gasf*(Pr.^2)))))  == 0];% equality constraint 

%Constraints = [Constraints,Si+shedg-d-((Cft*(sign(Gasf*(Pr.^2)).*sqrtm(abs(Gasf*(Pr.^2))))))  == 0];    % equality constraint 

% Constraints - Generation constraint
Constraints = [Constraints,Simin <= Si <= Simax];                          % fix gas generation maximum and minimum

% Constraints - Pressure constraint
Constraints = [Constraints,Prmin <= Pr <= Prmax];                          % fix maximum gas generation
Constraints = [Constraints,Pr(7,:) >= 62];                                 % fix maximum gas generation

% Constraints - Storage
Constraints = [Constraints,Psto(:,2:horizon).*Stovec(:,2:horizon)-Psto(:,1:horizon-1).*Stovec(:,1:horizon-1)==Pin(:,2:horizon)-Pout(:,2:horizon)];
Constraints = [Constraints,Psto(:,1).*Stovec(:,1)==Pinit+Pin(:,1)-Pout(:,1)];
Constraints = [Constraints, 0<=Psto.*Stovec<=Stomax];
Constraints = [Constraints, 0==Psto.*(1-Stovec)];
Constraints = [Constraints,0<=Pin<=10];
Constraints = [Constraints,0<=Pout<=10];

% Constraints - Shedding
Constraints = [Constraints,shedg >= 0];                                    % fix shedding to positive
Constraints = [Constraints,shedg <= d];                                    % fix shedding to maximum demand

% Objective 
Objective = 0;                                                             % set objective to zero

Objective= sum(C1'*Si)+sum(sum(shedg.*shedcost))+sum(sum((Pin.*Stoincost)));                           % objective function

% Optimize
if horizon==0
options= sdpsettings('solver','fmincon', 'verbose',1,'usex0',1,'allownonconvex',1,'savesolveroutput',1);
else
options= sdpsettings('solver','bmibnb','verbose',1,'usex0',1,'savesolveroutput',1,'bmibnb.lowersolver',...
    'CPLEX','bmibnb.uppersolver', 'fmincon','bmibnb.lpsolver', 'CPLEX','bmibnb.lpreduce',1,'bmibnb.roottight',1);
end

options.fmincon.MaxIter=3000;                                            
options.fmincon.MaxFunEvals=1000000;
options.fmincon.TolX=1E-6;
options.fmincon.Algorithm= 'interior-point';
options.bmibnb.vartol=1E-6;

exitflag = optimize(Constraints,Objective,options)   

[Johan Löfberg]

YALMIP always tries to use the relaxation solution to create an upper bound (i.e., a feasible solution). If you claim that the relaxation supplied is feasible, it is very strange (although I don't see how you would be able to see that, as YALMIP doesn't return the relaxation).

The fact that the lower bound is equal to the known globally optimal solution, does not mean that the relaxation is feasible. It could be that it is easy to create a strong relaxation value, but actually finding a solution is hard

Discussed further by email, and two improvements on the model was to work with the squared Pr variable instead, and possibly using a pwa approximation with sos2 models to obtain a MILP

[Michiel Jansen]

Thanks a lot Johan it worked! I will apply adaptive gridding to make it more accurate.