Solver not applicable (cplex & gurobi)

[Jasper Meyer]

Hi,

I'm writing my master thesis about the economics of energy storages. Now I have the Problem, that for the two solver (cplex & gurobi) which I was told to use I get the error that they are not applicable for my problem of calculating the levelized cost of energy (objective possibility #3 at the end of the code), costs divided by the energy leaving the storage). But on the other hand both solvers can optimize the two input terms 'investment cost' and 'energy' individually (objective possibilities 1&2). I simplified the code so it is easier to read and hope one of you has an idea what to do as I am pretty inexperienced using matlab. 

Thanks in advance for anybody who can help me with this.

(please note: time steps are 1 hour so the power of a period equals the energy)

--------------------------------------------------

clear

clc

ops=sdpsettings('solver','cplex');

T = 10; % # of periods 

E_ex = randn(T,1)*100; % create random demand function (E_ex > 0: store, E_ex < 0: discharge).

eta = 0.9; % roundtrip efficiency

eta_in = eta^0.5; % efficiency of storing process (square root of roundtrip efficiency)

eta_out = eta^0.5; % efficiency of discharging process (square root of roundtrip efficiency)

sc_Pin = 10; % specific costs of charging-power ($/kW)

sc_Pout = 10; % specific costs of discharging-power ($/kW)

sc_E = 11; % specific costs of storage capacity ($/kWh)

E = sdpvar(T,1); % stored energy of each time step t. E(1,1) = initial stored Energy E_0

E_ex_rest = sdpvar(T,1); % energy that the storage can't store / demand the storage can't meet 

E_out = sdpvar(T,1); % discharged energy in step t

E_in = sdpvar(T,1); % stored energy in step t

E_ex_eff_in = sdpvar(T,1); % the energy that the storage is absorbing from the grid.

E_ex_eff_out = sdpvar(T,1); % the energy that the storage is emitting to the grid.

CE_out_ex = sdpvar(1,1); %Sum of the emitted energy between t and T to calculate LCOE (Levelized Cost of Energy)

P_in = sdpvar(1,1); %charging-power

P_out = sdpvar(1,1); %discharging-power

E_max = sdpvar(1,1); %storage capacity

E_min = 0; %minimal storage level

C_inv = sdpvar(1,1); %investment costs

LCOE = sdpvar(1,1); %levelized cost of energy

t = 1;

constraints = [];

%Definittion of initial values

E(1,1) = 0;

E_in(1,1)=0;

E_out(1,1)=0;

E_ex_eff_out(1,1) %set initial value one to enable the emitted energy of the other time steps to become zero (otherwise maybe error due to LCOE = C_inv / sum(E_ex_eff_out))

for t = 2:T  

 

   % No demand from grid

   if E_ex(t,1) == 0                                                          

       E_ex_eff_in(t,1) = 0;

       E_ex_eff_out(t,1) = 0;

       E(t,1) = E(t-1,1);                                   

       E_out(t,1) = 0;

       E_in(t,1) = 0;

       E_ex_rest(t,1) = 0;

        

   % grid offers energy -> store

   elseif E_ex(t,1) > 0

        E_ex_eff_in(t,1) = E_ex(t,1); %as the storage should be able to store all the energy offered, the absorbed energy is equal to the offered

        E_in(t,1) = E_ex_eff_in(t,1) * eta_in; %Energy "arriving" in the storage

        E(t,1) = E(t-1,1) + E_in(t,1); %New energy in the storage

        E_ex_rest(t,1) = E_ex(t,1) - E_ex_eff_in(t,1); %residual energy that the storage can't store (in this case of full coverage = 0)

        

        

   % demand for energy from grid -> discharge 

   else 

       E_ex_eff_out(t,1) = min([abs(E_ex(t,1)),P_out,(E(t-1,1)-E_min)*eta_out]); %the demand which the storage can meet can be defined by three possible 'bottlenecks': the external demand, the discharging power or the energy in the storage (technical efficiency is considered).

       E_out(t,1) = E_ex_eff_out(t,1) / eta_out; % due to the efficiency the discharged energy is higher than the emitted 

       E(t,1) = E(t,1) - E_out(t,1); %new energy in the storage

       E_ex_rest(t,1) = E_ex(t,1) + E_ex_eff_out(t,1); %energy demand that the storage can't meet

       

   end

end

E_out(isnan(double(E_out))==1)=0;

E_ex_eff_out(isnan(double(E_ex_eff_out))==1)=0;

E_max = max(E); %installed capacity should be high enough to store all the offered energy 

P_in = max(E_ex_eff_in); %the installed charging power should be high enough to store all the offered energy

CE_out_ex = sum(E_ex_eff_out); %sum of the emitted energy

C_inv = sc_Pin * P_in + sc_Pout * P_out + sc_E * E_max; %investment costs

LCOE = C_inv / CE_out_ex; %levelized cost of energy

constraints = [constraints, E(:)>=0];

constraints = [constraints, P_in>=0];

constraints = [constraints, P_out>=0];

constraint = [constraints, E_ex_eff_out(:)>=0];

constraints = [constraints, E_ex_rest(:)<=0]; % All energy must be stored

%objective = C_inv; %#1

%objective =  -CE_out_ex; %#2

objective =  LCOE; %#3

solve = optimize(constraints, objective,ops)

[Johan Löfberg]

Your objective LCOE is C_inv /CE_out_ex where both depend on decision variables. That will not work with cplex and gurobi ,as they work with linear or convex quadratic representable objectives. Your objective is clearly neither.

You can try a bisection. Introduce a constant t and solve a feasibility problem with the constraint C_inv <= t*CE_out_ex (assuming CE_out_ex>=0). Bisect the value t to optimality (see DecayRate example on the Wiki)

BTW, this

constraint = [constraints, E_ex_eff_out(:)>=0];

should probably be

constraints = [constraints, E_ex_eff_out(:)>=0];

[Jasper Meyer]

Thank you very much for the answer. Is there a different suitable solver for the way how I wrote the code?

[Johan Löfberg]

You can throw it at any nonlinear solver. However, a bisection using a LP solver will work nicely and you have guaranteed global optimality with such an approach (a nonlinear solver would most likely give you that too, but you never know as it probably depends a bit on the form of the model YALMIP will send to the solver)

[Jasper Meyer]

Unfortunately, I am pretty inexperienced with yalmip (& also matlab in general). So when I tried to apply the decay example to my problem, not surprisingly I was unable to implement the bisection. I wasn't even sure where in my code to implement the additional calculation. I just did it at the end. I am pretty sure I made some more mistakes as I wasn't sure what the constraint P>=eye(2) does and if I need some similar constraint in my problem, either. 

Constraints = [P>=eye(2), A'*P+P*A <= -2*t*P];

I also wasn't sure if I should include the other Constraints or only include the C_inv <= b*CE_out_ex constraint (note: I used b instead of t as the variable as t is already used).

So I just added the following code behind the code I posted above. 

b_lower = 0;

b_upper = 1;

sdpvar b

ops = sdpsettings('solver','cplex');

Constraints = [C_inv <= b*CE_out_ex];

Solver = optimizer(Constraints,[],ops,b,LCOE);

[~, errorcode] = Solver{b_upper};

while ~(errorcode==1)

    b_upper = b_upper*2;

    [~, errorcode] = Solver{b_upper};

end

b_works = b_lower;

while (b_upper-b_lower)>tol

    b_test = (b_upper+b_lower)/2;

    disp([b_lower b_upper b_test])

    [LCOEopt, errorcode] = Solver{b_test};

    if errorcode==1

        b_upper = b_test;

    else

        b_lower = b_test;

        b_works = b_test;

        LCOEworks = LCOEopt;

    end

end

F = [C_inv <= b_upper*CE_out_ex];

ops = sdpsettings('verbose',0,'warning',0);

sol = optimize(F,[],ops);

while ~(sol.problem==1)

    b_upper = b_upper*2;

    F = [C_inv <= b_upper*CE_out_ex];

    sol = optimize(F,[],ops);

end

bol = 0.01;

b_works = b_lower;

while (b_upper-b_lower)>bol

  b_test = (b_upper+b_lower)/2;

  disp([b_lower b_upper b_test])

  F = [C_inv <= b_upper*CE_out_ex];

  sol = optimize(F,[],ops);

  if sol.problem==1

    b_upper = b_test;

  else

    b_lower = b_test;

    b_works = b_test;

    C_inv_works = value(C_inv);

    CE_out_ex_works = value(CE_out_ex);

 end

end

When running the script I get the error:

 

Error using yalmip2cplex (line 25)

Nonlinear monomials remain when calling CPLEX. If you are using OPTIMIZER, ensure your model really is solvable by

CPLEX for fixed parameters. If you still think so, please report this and ask for a feature improvement.

Error in export (line 160)

        model = yalmip2cplex(interfacedata);

Error in optimizer (line 200)

    [aux1,aux2,aux3,model] = export((x == repmat(pi,nIn*mIn,1))+Constraints,Objective,options,[],[],0);

Error in calc_1storage_v2_forum_engl_bisec (line 104)

Solver = optimizer(Constraints,[],ops,b,LCOE);

Thanks again for the help so far.

[Johan Löfberg]

Why are you throwing away all the constraints in the model? You seem to be missing the point. To minimize f(x)/g(x), we minimize t subject to f(x)<=t*g(x). Now we fix t and simply solve feasibility problems for various t, but we do this cleverly by bisection to do as few tries as possible

constraints = [constraints, E(:)>=0];
constraints = [constraints, P_in>=0];
constraints = [constraints, P_out>=0];

constraint = [constraints, E_ex_eff_out(:)>=0];

constraints = [constraints, E_ex_rest(:)<=0]; % All energy must be stored

% Find an upper bound on the optimal value. Solve problems until feasible
upper = 1
s = solvesdp([C_inv <= upper*CE_out_ex,constraints]);
while s.problem
    upper = upper*2
    s = solvesdp([C_inv <= upper*CE_out_ex,constraints]);
end
% We can never do better than 0
lower = 0;
while upper-lower > 0.01
    test = (upper + lower)/2;
    s = solvesdp([C_inv <= test*CE_out_ex,constraints]);
    if s.problem == 0
        tworks = test;
        upper = test;
    else
        lower = test;
    end
end
% solve again at the optimal t so you can extract solutions
solvesdp([C_inv <= tworks*CE_out_ex,constraints]);
value(E)
value...

[Johan Löfberg]

The crash in optimizer is weird though. I will have to investigate that further. It appears to be related to some of the mixed-integer parts of the model being hidden when the temporary model with parameterized b is created. Using the optimizer approach is not necessary though, it just speeds up and cleans up the code a bit

[Johan Löfberg]

OK, found out that it is a bug when cplex is used with nonlinearly parameterized object in optimizer. Remove the offending line in yalmip2cplex as a quick fix.

[Jasper Meyer]

Thanks a lot for your support. It is working now. Just to be sure I am understanding it correctly: the bisection leads to the point of where LCOE is optimal when the last 'tworks' is reached. What would you define the 0.01 criteria as? As the 'confidence level'?

And I am still a little bit confused with the optimizer crash you mentioned as it wasn't included in the code you posted. Do I understand this also correct, that the code you posted leads to the optimal solution but using optimizer instead of solvesdp could maybe lead to a speed up?

[Johan Löfberg]

0.01 is the optimality tolerance. It is impossible to generate a solution with the cost lower, but it is possible with upper. Hence, we now the optimal solution has a cost in between. Which tolerance you want is your choice

No, I coded it without to make it very simple and clear what we are doing. Using optimizer removes some YALMIP overhead and is thus faster (and makes the code more compact)

[Jasper Meyer]

Thank you once more very much for your help.