Solver not applicable (gurobi) - on/off variables

[Tommaso De Marco]

Dear Professor,

I'm trying to model an Energy Hub (photovoltaic + energy storage) and to optimize its scheduling. I'm having problems when modeling the electrical energy storage: in particular I'm modelling the maximum energy that can be charged/discharged as a variable and also the non simultaneity of charge and discharge using two binary variables. So the upper bound of the charging flow, for example, end up being the product of a variable and a binary variable. I get a warning stating "Solver not applicable (gurobi) when running the optimization. I'm pretty inexperienced in the field of optimization and solvers but I read the "Unit commitment" tutorial on github realizing that this product might be a problem (as a matter of fact in the tutorial the max/min values of the variables were constants). Is this the problem in this case? Do you have suggestions on how to fix this problem and on how to model the non simultaneity of charge and discharge keeping the model a MILP so that Gurobi can solve it?

(The complete model is attached).

Thank you very much for your attention!

%EES VARIABLES AND CONSTRAINTS

EES_capacity = 50;              %energy capacity of the battery [kWh]

EES_selfdis = 0.001;            %battery selfdischarge (%/100/h)

EES_loss_ch = 0.1;               %charging loss

EES_loss_dis = 0.1;              %dicharging loss

%%%%%%flows variables

Pphub_ees = sdpvar(1,Horizon);      %variable "from power hub to battery"

Pees_phub = sdpvar(1,Horizon);      %variable "from battery to power hub"

Pees_store = sdpvar(1,Horizon+1);   %variable "energy stored" (the 25th step relates to the 24th step of other variables)

%%%%%%binary variables

delta_ees_ch = diag(binvar(Horizon,1));     %battery input flow bin var

delta_ees_dis = diag(binvar(Horizon,1));    %battery output bin var 

ees_binvar_con = (delta_ees_ch + delta_ees_dis)<=1; %constraint on bin vars (charge/discharge not at same time)

%%%%%%in constraints

Pphub_ees_ch_max_start = EES_capacity - Pees_store(1);  %maximum energy that can be charged at the first timestep

Pphub_ees_ch_max = EES_capacity - Pees_store(2:end-1);  %maximum energy that can be charged for timesteps from 2 to 24 regarding the input flow

Pphub_ees_con = [Pphub_ees>=0,Pphub_ees(1)<=Pphub_ees_ch_max_start*delta_ees_ch(1,1),Pphub_ees(2:end)<=Pphub_ees_ch_max*delta_ees_ch([2:end],[2:end])];  %charge                                           constraint: positivity/=0 + smaller than the energy that can fit in the battery + non simultaneity

%%%%%%out constraints

Pees_phub_con = [Pees_phub>=0,Pees_phub<=Pees_store(1:end-1)*delta_ees_dis];   %discharge constraint: positivity/=0 and smaller than the previously stored energy until that moment

%%%%%%store constraints

Pees_store_con = [Pees_store>=0,Pees_store<=EES_capacity];  %energy stored has to be positive/=0 and smaller than the capacity

Pees_store_cont =[Pees_store(2:end)==(1-EES_selfdis)*Pees_store(1:end-1)+(1-EES_loss_ch)*Pphub_ees - (1/(1-EES_loss_dis))*Pees_phub,Pees_store(1)==0];    %continuity: energy at timestep t+1 is equal to the energy stored at timestep t plus all charges and minus discharges

%%%%%%overall constraint

EES_constraints = [Pphub_ees_con,Pees_phub_con,Pees_store_con,Pees_store_cont,ees_binvar_con];  %EES constraints collection

 

[Johan Löfberg]

You are conducting the cardinal sin of multiplying binaries with continuous variables. Products a <= x*y  leads to nonconvex models and cplex/gurobi etc don't support this This is explained in the unitcommitment example asthe first tip. Never ever do this

If y is binary, you have a <= w, implies(y == 0, w == 0) + implies(y == 1, w == x), it easily linearizable

If you have a <= x*y where both x and y are continuous, you are in deep troubles, as that isn't linearizable, and intrincically nonconvex and hard

[Tommaso De Marco]

Dear Professor,

I've modified my code using your suggestions. So I introduced two new variables to avoid the product of continuos variable and binary variable (these two new variable represent the upper bounds of the charging flow and discharging flow). I then tried to use the "implies" in order  to model the non simultaneity of charge and discharge. I also introduced boundaries for these two new variables since I was getting a warning saying "Warning: You have unbounded variables in IMPLIES leading to a lousy big-M relaxation". 

I still get a series of errors that I honestly don't understand. Do you have suggestions on how to fix the problem and on how to deal with the time-varying bounds of charge and discharge together with the non-simultaneity?

(the model is attached)

Thank you very much for your attention.

%EES VARIABLES AND CONSTRAINTS

EES_capacity = 50;              %energy capacity of the battery [kWh]

EES_selfdis = 0.001;            %battery selfdischarge (%/100/h)

EES_loss_ch = 0.1;               %charging efficiency

EES_loss_dis = 0.1;              %dicharging efficiency

SOCmin = 0.15*EES_capacity;     %minimum state of charge [kWh]

%%%%%%flows variables

Pphub_ees = sdpvar(1,Horizon);      %variable "from power hub to battery"

Pees_phub = sdpvar(1,Horizon);      %variable "from battery to power hub"

Pees_store = sdpvar(1,Horizon+1);   %variable "energy stored" (the 25th step relates to the 24th step of other variables)

%%%%%%binary variables

delta_ees_ch = diag(binvar(Horizon,1));     %battery input flow bin var

delta_ees_dis = diag(binvar(Horizon,1));    %battery output bin var 

ees_binvar_con = (delta_ees_ch + delta_ees_dis)<=1; %constraint on bin vars (charge/discharge not at same time)

AvL_capacity = sdpvar(1,Horizon);   %new variable (available capacity/energy that can be stored)

%%%%%%in constraints

AvL_capacity_start = EES_capacity - Pees_store(1);  %maximum energy that can be charged at the first timestep

AvL_capacity_rest = EES_capacity - Pees_store(2:end-1);  %maximum energy that can be charged for timesteps from 2 to 24 regarding the input flow

AvL_capacity_bounds = [AvL_capacity>=0,AvL_capacity<=EES_capacity]; %%%boundaries of the new "available capacity" variable

Pphub_ees_con = [Pphub_ees>=0];     %charge constraint: positivity/=0 

Pphub_ees_con_start_0 = [Pphub_ees(1)<=AvL_capacity(1),implies(delta_ees_ch(1,1)==0,AvL_capacity(1)==0)];   %the upper bound of the charge flow is the available capacity                                                                                                                                                                                               (the AvL_capacity depends on the binvar

Pphub_ees_con_start_1 = [Pphub_ees(1)<=AvL_capacity(1),implies(delta_ees_ch(1,1)==1,AvL_capacity(1)==AvL_capacity_start)];

Pphub_ees_con_rest_0 = [Pphub_ees(2:end)<=AvL_capacity(2:end),implies(delta_ees_ch([2:end],[2:end])==0,AvL_capacity(2:end)==0)];

Pphub_ees_con_rest_1 = [Pphub_ees(2:end)<=AvL_capacity(2:end),implies(delta_ees_ch([2:end],[2:end])==1,AvL_capacity(2:end)==AvL_capacity_rest)];

Pphub_ees_con_overall = [AvL_capacity_bounds,Pphub_ees_con,Pphub_ees_con_start_0,Pphub_ees_con_start_1,Pphub_ees_con_rest_0,Pphub_ees_con_rest_1];  %overall ees constraint

%%%%%%out constraints

AvL_energy = sdpvar(1,Horizon);     %new variable (amount of energy we can discharge)

AvL_energy_bounds = [AvL_energy>=0,AvL_energy<=EES_capacity];   %boundaries of the new dischargeable energy variable (either 0 or equal to the max capacity)

Pees_phub_con = [Pees_phub>=0];   %discharge constraint: positivity/=0 

Pees_phub_con_rest_0 = [Pees_phub<=AvL_energy,implies(delta_ees_dis==0,AvL_energy==0)]; %the upper bound of the discharge flow is the new variable AvL_energy and                                                                                                                                                                        depends on the binvar

Pees_phub_con_rest_1 = [Pees_phub<=AvL_energy,implies(delta_ees_dis==1,AvL_energy==Pees_store(1:end-1))];

Pees_phub_con_overall = [AvL_energy_bounds,Pees_phub_con,Pees_phub_con_rest_0,Pees_phub_con_rest_1];    %overall constraint

%%%%%%store constraints

Pees_store_con = [Pees_store>=0,Pees_store<=EES_capacity];  %energy stored has to be positive/=0 and smaller than the capacity

Pees_store_cont =[Pees_store(2:end)==(1-EES_selfdis)*Pees_store(1:end-1)+(1-EES_loss_ch)*Pphub_ees - (1/(1-EES_loss_dis))*Pees_phub,Pees_store(1)==0];    %continuity: energy at timestep t+1 is equal to the energy stored at timestep t plus all charges and minus discharges

%%%%%%overall constraint

EES_constraints = [Pphub_ees_con_overall,Pees_phub_con_overall,Pees_store_con,Pees_store_cont,ees_binvar_con];  %EES constraints collection

[Johan Löfberg]

First, why are you adding Pphub_ees(1)<=AvL_capacity(1) twice

Second, implies(delta_ees_ch([2:end],[2:end])==0,AvL_capacity(2:end)==0) doesn't make sense as delta_ees_ch([2:end],[2:end]) is matrix (of which almost all elements are zeros) and the other object is a vector . Why are you even defining that as a matrix, when it really only has elements on the diagonal, and then you try to reference precisely that diagonal

When I try to run the code I get dimension errors in an implies operator, but that most likely comes from your weird diagonal matrices

[Johan Löfberg]

similarily you add Pees_phub<=AvL_energy twice. It looks like you think explicit bonds on variables have to be defined in the same list as the logic constraints where the bounds are used. The explicit bounds only have to appear somewhere in the complete model

[Tommaso De Marco]

I added them twice because I thought it was necessary to write the inequality for the bound when using implies.

You're right, using diagonal matrices doesn't make sense.

So, I've modified them so now the binary variables are just vectors (as the other variables). I've also modified the lines containing the implies and the inequalities Pphub_ees(1)<=AvL_capacity(1) and Pphub_ees(2:end)<=AvL_capacity(2:end) so that now they're used just once. I also moved the bounds from the lines containing implies.

I still get the same errors and I don't get why. Do you have any suggestion or ideas on where to look to solve the issue?

(the model is attached)

Thanks for your help!!

%EES VARIABLES AND CONSTRAINTS

EES_capacity = 50;              %energy capacity of the battery [kWh]

EES_selfdis = 0.001;            %battery selfdischarge (%/100/h)

EES_loss_ch = 0.1;               %charging efficiency

EES_loss_dis = 0.1;              %dicharging efficiency

SOCmin = 0.15*EES_capacity;     %minimum state of charge [kWh]

%%%%%%flows variables

Pphub_ees = sdpvar(1,Horizon);      %variable "from power hub to battery"

Pees_phub = sdpvar(1,Horizon);      %variable "from battery to power hub"

Pees_store = sdpvar(1,Horizon+1);   %variable "energy stored" (the 25th step relates to the 24th step of other variables)

%%%%%%binary variables

delta_ees_ch = binvar(1,Horizon);     %battery input flow bin var

delta_ees_dis = binvar(1,Horizon);    %battery output bin var 

ees_binvar_con = (delta_ees_ch + delta_ees_dis)<=1; %constraint on bin vars (charge/dischargege not at same time)

%%%%%%in constraints

AvL_capacity = sdpvar(1,Horizon);   %new variable (available capacity/energy that can be stored)

AvL_capacity_start = EES_capacity - Pees_store(1);  %maximum energy that can be charged at the first timestep

AvL_capacity_rest = EES_capacity - Pees_store(2:end-1);  %maximum energy that can be charged for timesteps from 2 to 24 regarding the input flow

AvL_capacity_bounds = [AvL_capacity>=0,AvL_capacity<=EES_capacity]; %%%boundaries of the new "available capacity" variable

Pphub_ees_con = [Pphub_ees>=0];     %charge constraint: positivity/=0 

Pphub_ees_con_start = [Pphub_ees(1)<=AvL_capacity(1)];   %the upper bound of the charge flow is the available capacity (the AvL_capacity depends on the binvar

Pphub_ees_con_start_bin = [implies(delta_ees_ch(1)==0,AvL_capacity(1)==0),implies(delta_ees_ch(1)==1,AvL_capacity(1)==AvL_capacity_start)];

Pphub_ees_con_rest = [Pphub_ees(2:end)<=AvL_capacity(2:end)];

Pphub_ees_con_rest_bin = [implies(delta_ees_ch(2:end)==0,AvL_capacity(2:end)==0),implies(delta_ees_ch(2:end)==1,AvL_capacity(2:end)==AvL_capacity_rest)];

Pphub_ees_con_overall = [AvL_capacity_bounds,Pphub_ees_con,Pphub_ees_con_start,Pphub_ees_con_rest,Pphub_ees_con_start_bin,Pphub_ees_con_rest_bin];  %overall ees constraint

%%%%%%out constraints

AvL_energy = sdpvar(1,Horizon);     %new variable (amount of energy we can discharge)

AvL_energy_bounds = [AvL_energy>=0,AvL_energy<=EES_capacity];   %boundaries of the new dischargeable energy variable (either 0 or equal to the max capacity)

Pees_phub_con = [Pees_phub>=0];   %discharge constraint: positivity/=0 and smaller than the stored energy until that moment

Pees_phub_con_rest = [Pees_phub<=AvL_energy]; %the upper bound of the discharge flow is the new variable AvL_energy and depends on the binvar

Pees_phub_con_rest_bin = [implies(delta_ees_dis(1:end)==0,AvL_energy(1:end)==0),implies(delta_ees_dis(1:end)==1,AvL_energy==Pees_store(1:end-1))];

Pees_phub_con_overall = [AvL_energy_bounds,Pees_phub_con,Pees_phub_con_rest,Pees_phub_con_rest_bin];    %overall constraint

%%%%%%store constraints

Pees_store_con = [Pees_store>=0,Pees_store<=EES_capacity];  %energy stored has to be positive/=0 and smaller than the capacity

Pees_store_cont =[Pees_store(2:end)==(1-EES_selfdis)*Pees_store(1:end-1)+(1-EES_loss_ch)*Pphub_ees - (1/(1-EES_loss_dis))*Pees_phub,Pees_store(1)==0];    %continuity: energy at timestep t+1 is equal to the energy stored at timestep t plus all charges and minus discharges

%%%%%%overall constraint

EES_constraints = [Pphub_ees_con_overall,Pees_phub_con_overall,Pees_store_con,Pees_store_cont,ees_binvar_con];  %EES constraints collection

[Johan Löfberg]

I don't get any issues about solver not applicable, as it crashes already during compilation

That is due to a bug in implies, where vectorized implications doesn't work for row vectors. You must convert to column vectors to get it to work

fails

x = sdpvar(1,2);y = sdpvar(1,2);optimize(implies(x==0,y==1))

works

x = sdpvar(1,2);y = sdpvar(1,2);optimize(implies(x'==0,y'==1))

[Johan Löfberg]

and the code looks unnecessarily complicated

cap = [AvL_capacity_start AvL_capacity_rest]

Pphub_ees_con_start_bin = [implies(delta_ees_ch==0,AvL_capacity==0),implies(delta_ees_ch==1,AvL_capacity==cap)];

[Johan Löfberg]

Also note that x == 0 implies y == 0 means if all(x==0) then all(y==0), i.e, implies(x==0,y==0) is expanded to implies(x==0,d) + implies(d,y==0) where d is a scalar binary.

You should not use implications from constraints, but do the logical modelling manually using explicit binary variables, to avoid any such ambiguity.

[Tommaso De Marco]

Dear Professor,

I've modified my code considering the first two of your last suggestions so that now I converted to column vectors the vectors present in the implies and I simplified the code as you suggested. I didn't understand your last suggestion though. 

The fact is that I have to model the upper bound of the variable Pphub_ees through the inequality Pphub_ees<=AvL_capacity (I introduced the AvL_capacity to avoid the product of sdpvar and binvar as you suggested previously). The AvL_capacity assumes different values for each timestep so it's a variable (1x24 row vector). Nonetheless there is the binary variable delta_ees_ch (1x24 row vector) that is used to say if the input flow to the battery is present or not (if delta_ees_ch=1 the input flow is present then the output flow is not present so delta_ees_dis=0). If the first element (element associated to the first timestep) of delta_ees_ch is ==0 then the AvL_capacity has to be zero so that no input flow (Pphub_ees) is present for that timestep since in this situation we would have 0<=Pphub_ees<=0; if the first element of delta_ees_ch is equal to 1 then AvL_capacity is equal to the first element of AvL_capacity_cap. So I should have this kind of implication for each timestep (both for input flow: Pphub_ees and for output flow:Pees_phub). So for example, for each timestep (j) I have to have the following condition:

given Pphub(j)<=AvL_capacity(j), if delta_ees_ch(j)==0 then AvL_capacity(j)==0 if delta_ees_ch(j)==1 then AvL_capacity(j)=AvL_capacity_cap(j)

At this point the code doesn't show errors but I don't know if the Pphub_ees_con_avl and Pphub_ees_con_bin (together with Pees_phub_con_rest and Pees_phub_con_rest_bin regarding the output flow) are well written through the implies in my code and actually represent the implications between binary variables, Pphub_ees and Pees_phub for each timestep as I've tried to explain with this message.

Do you think the code correctly represent the variability of AvL_capacity for each timestep and the non simultaneity of Pphub_ees and Pees_phub for each timestep? If not, how do you think the code shoul be modified?

Thank you very much for your attention and help, it's very appreciated!!

(the model is attached)

%EES VARIABLES AND CONSTRAINTS

EES_capacity = 50;              %energy capacity of the battery [kWh]

EES_selfdis = 0.001;            %battery selfdischarge (%/100/h)

EES_loss_ch = 0.1;               %charging efficiency

EES_loss_dis = 0.1;              %dicharging efficiency

SOCmin = 0.15*EES_capacity;     %minimum state of charge [kWh]

%%%%%%flows variables

Pphub_ees = sdpvar(1,Horizon);      %variable "from power hub to battery"

Pees_phub = sdpvar(1,Horizon);      %variable "from battery to power hub"

Pees_store = sdpvar(1,Horizon+1);   %variable "energy stored" (the 25th step relates to the 24th step of other variables)

%%%%%%binary variables

delta_ees_ch = binvar(1,Horizon);     %battery input flow bin var

delta_ees_dis = binvar(1,Horizon);    %battery output bin var 

ees_binvar_con = (delta_ees_ch + delta_ees_dis)<=1; %constraint on bin vars (charge/dischargege not at same time)

%%%%%%in constraints

AvL_capacity = sdpvar(1,Horizon);   %new variable (available capacity/energy that can be stored)

AvL_capacity_start = EES_capacity - Pees_store(1);  %maximum energy that can be charged at the first timestep

AvL_capacity_rest = EES_capacity - Pees_store(2:end-1);  %maximum energy that can be charged for timesteps from 2 to 24 regarding the input flow

AvL_capacity_bounds = [AvL_capacity>=0,AvL_capacity<=EES_capacity]; %%%boundaries of the new "available capacity" variable

AvL_capacity_cap = [AvL_capacity_start,AvL_capacity_rest];

Pphub_ees_con = [Pphub_ees>=0];     %charge constraint: positivity/=0 

Pphub_ees_con_avl = [Pphub_ees<=AvL_capacity];

Pphub_ees_con_bin = [implies(delta_ees_ch(1:end)'==0,AvL_capacity(1:end)'==0),implies(delta_ees_ch(1:end)'==1,AvL_capacity(1:end)'==AvL_capacity_cap(1:end)')];

Pphub_ees_con_overall = [AvL_capacity_bounds,AvL_capacity_cap,Pphub_ees_con,Pphub_ees_con_avl,Pphub_ees_con_bin];  %overall ees constraint

%%%%%%out constraints

AvL_energy = sdpvar(1,Horizon);     %new variable (amount of energy we can discharge)

AvL_energy_bounds = [AvL_energy>=0,AvL_energy<=EES_capacity];   %boundaries of the new dischargeable energy variable (either 0 or equal to the max capacity)

Pees_phub_con = [Pees_phub>=0];   %discharge constraint: positivity/=0 and smaller than the stored energy until that moment

Pees_phub_con_rest = [Pees_phub<=AvL_energy]; %the upper bound of the discharge flow is the new variable AvL_energy and depends on the binvar

Pees_phub_con_rest_bin = [implies(delta_ees_dis(1:end)'==0,AvL_energy(1:end)'==0),implies(delta_ees_dis(1:end)'==1,AvL_energy(1:end)'==Pees_store(1:end-1)')];

[Johan Löfberg]

simply don't use the pattern implies(binary == 1, action), but write it as implies(binary,action) or implies(1-binary == action)

The case with a vector constraint implications is partly ambiguous, and partly buggy