'Solver not applicable (gurobi)'

[mON HC]

Dear Johan,

I am trying to solve an optimisation problem where I want to find the price (f) that will maximise revenues from demmand (Di). I have a piecewise function of the price based on demmand. I did the formulation using implies, the problem I have is that when I want to maximize price*demand this does not seem to be working. It works when I just maximise f but this is not my calculation of my revenues though. I have tried to use FMINCON but it does not work either, unless I remove the implies and just substitute with the nonlinear equivalent equation.

Any help with this error would be highly appreciated.

Here is my code:

%DATA
y1=[28 23.5847];
x1=[0.1721 0.24];
m1 = (y1(1)-y1(2))/(x1(1)-x1(2));

y2=[23.5847 11.1145];
x2=[0.24 0.28];
m2 = (y2(1)-y2(2))/(x2(1)-x2(2));

y3=[11.1145 10];
x3=[0.28 0.3961];
m3 = (y3(1)-y3(2))/(x3(1)-x3(2));


f = sdpvar(1,24);
Di = sdpvar(1,24);

%CONSTRAINTS
Model = [];

%Model = [ Di(1)+Di(2)+Di(3)+Di(4)==0.5, Di(23)+Di(24)==0.5 ]
for t=1:24
d = binvar(3,1);
Model = [Model,
    sum(d) == 1,
implies(d(1), [  0.1721<= Di <= 0.24, f == m1*(Di-x1(1))+y1(1)]);
implies(d(2), [ 0.24 <= Di <= 0.28, f == m2*(Di-x2(1))+y2(1)]);
implies(d(3), [ 0.28 <= Di <= 0.3961,  f == m3*(Di-x3(1))+y3(1)])];
end


Model = [Model, Di(3)+Di(4)+Di(5)==0.9 ]

Objective = sum(f.*Di,2)


options = sdpsettings('solver','gurobi');
   
sol = optimize(Model,-Objective,options);
   
% Analyze error flags
if sol.problem == 0
    % Extract and display value
    solution = value(Di)
    solution2 = value(f)
    else
    display('Hmm, something went wrong!');
    sol.info
    yalmiperror(sol.problem)
    end

NOW CODE WITH FMINCON

%DATA INPUTS FOR MODEL
T=24; %work timeframe for EV charging station(think if it is a parking lot or charging station only, residential charging?)

%IMPORT DATA AND CREATE PARAMETERS
%param=[9.76584915402396,27.9904155442713,0.252459657607282,-39.8453332837601];

   
%DEFINE VARIABLES
M = sdpvar(1,T);
%Di = sdpvar(1,T);
 
Di =  (70368744177664.*log(-(562949953421312.*M - 15757203126890695)./(562949953421312.*M - 5497684326377347)))./(2803866064518675*log(10)) + 4547908879704953/18014398509481984
  
   
   
Constraints = [10<=M(:)<=28, Di(2)+Di(3)+Di(4)==0.5];
   
Objective = sum(Di.*M,2);
   
options = sdpsettings('solver','fmincon');
   
sol = optimize(Constraints);
   
    % Analyze error flags
    if sol.problem == 0
    % Extract and display value
    solution = value(Di)
    solution2 = value(M)
    solution3 = value(Objective)
    else
    display('Hmm, something went wrong!');
    sol.info
    yalmiperror(sol.problem)
    end
  
   

[Johan Löfberg]

You have a bilinear objective fi*Di, hence nonconvex and not solvable using gurobi

fmincon does not support integer variables

YALMIPs built-in global solver bmibnb easily solves the problem though

(you have to add bounds on all variables involved in the logic constraints. At the moment YALMIP warns a lot about that)

[Mark L. Stone]

Also solvable by CPLEX via YALMIP. And solvable by direct call from MATLAB to Gurobi 9.x, even though not via YALMIP.

[Johan Löfberg]

Ah right always forget about that.

Set cplex.optimalitytarget to 3 and it solves easily

[mON HC]

Dear Johan,

I hope you are well. I have updated the number of variables in my model. However I get the an error: simulation exceeds iterations or simulation time. I was wondering if there is a way to define the number of iterations with bmibnb?. I have tried using cplex but it  does not seem to be working. Here is the only modification I did:

T=24

for t=1:T

d = binvar(3,1);
Model = [Model,
    sum(d) == 1,

implies(d(1), [  0.1721<= Di(t) <= 0.24, f(t) == m1*(Di(t)-x1(1))+y1(1)]);
implies(d(2), [ 0.24 <= Di(t) <= 0.28, f(t) == m2*(Di(t)-x2(1))+y2(1)]);
implies(d(3), [ 0.28 <= Di(t) <= 0.3961,  f(t) == m3*(Di(t)-x3(1))+y3(1)])];
end

Thanks a lot!

[Johan Löfberg]

you would have to supply reproducible code

[Johan Löfberg]

ok, saw now that you had that above

But then there is no change to the core situation. You are setting up a nasty nonlinear mixed-integer model

sure, you can change bmibnb.maxiter and bmibnb.maxtime, but it's still effectively intractable and you should spend time on a better model 

Is this really the cost you want (this is what you are maxiizing no)

x = 0.1721:0.001:0.24;plot(x, x.*(m1*(x-x1(1))+y1(1)));hold on

x = 0.24:0.001:0.28;plot(x, x.*(m2*(x-x2(1))+y2(1)));hold on

x = 0.28:0.001:0.3961;plot(x, x.*(m3*(x-x3(1))+y3(1)));hold on

if you really want this, I would go for an approximation using the interp1 operator, using sos2 etc

[mON HC]

Dear Johan,

Thank you so much for the message. Sorry no, these are the lines that I am trying to make as a piecewise function, I actually broke down a sigmoid because I thought this would be easier for the solver, I have also included the original sigmoid I am trying to approximate in the following code:

y1=[28 23.5847];
x1=[0.1721 0.24];
m1 = (y1(1)-y1(2))/(x1(1)-x1(2));
 
y2=[23.5847 11.1145];
x2=[0.24 0.28];
m2 = (y2(1)-y2(2))/(x2(1)-x2(2));
 
y3=[11.1145 10];
x3=[0.28 0.3961];
m3 = (y3(1)-y3(2))/(x3(1)-x3(2));
 
 

x = 0.1721:0.001:0.24;plot(x,(m1*(x-x1(1))+y1(1)));hold on
x = 0.24:0.001:0.28;plot(x,(m2*(x-x2(1))+y2(1)));hold on
x = 0.28:0.001:0.3961;plot(x,(m3*(x-x3(1))+y3(1)));hold on
 
param = [9.7658, 27.9904, 0.2525, -39.8453];
x = linspace(0.17,0.40,200);
fsigm = param(1)+(param(2)-param(1))./(1+10.^((param(3)-x)*param(4)))
y = fsigm+ 0.1*randn(size(x));
 
plot(x,y,'k')
hold on

Do you still suggest using the interp1 operator, using sos2 etc?

Sincerely,

Monica

[Johan Löfberg]

If all you want to do is to approximate that sigmoid (or the product sigmoid(D) * D), you should use the sos2 interp1 approximation, instead of reinventing the approximation manually, or even worse introducing bilinear products and integer variables

[mON HC]

Dear Johan,

Thanks for the message. I really appreciate it!.

Monica