NaNs in your constraints!

[Jasper Meyer]

Hi,

I know this problem has been issued multiple times, but I am not able to implement the solution to my model. I get the error code: "Error using compileinterfacedata (line 1014)

You have NaNs in your constraints!. Read more: http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Extra.NANInModel". Below you can find the model I am trying to solve (already shortened).

If I get it right,  s(n).C_OM(t,1) is an sdpvar and can't be implemented in the cost function. So I tried to solve it by building a concatanation (s(n).C_OM = [s(n).C_OM; p_f * s(n).sr_OM_var * s(n).E_ex_out(t,1)]). But now the problem seems to be infeasible (at least I won't get it solved after a couple of minutes), which it is not, if I don't include C_OM in the Cost function. Any suggestions how to proceed or any mistakes spotted?

clear

clc

constraints = [];

T = 50; %Totally simulated timesteps

t = 1;

cih = 1; %hours per timestep

p_el = 10; %price for electricity

p_f = 20; %fuel price for the variable operation costs

%create random demand for energy (storing and discharging)

E_dem = randn(T,1)*200;

for t=1:T 

       if E_dem(t,1) < 55

          if E_dem(t,1) > -55

             E_dem(t,1) = 0;

          end

       end

end

E_dem_in = zeros(T,1);

E_dem_out = zeros(T,1);

E_dem(1,1) = 0;

for t=2:T 

       if E_dem(t,1) == 0

           E_dem_in(t,1) = 0;

           E_dem_out(t,1) = 0;

       elseif E_dem(t,1) > 0

           E_dem_in(t,1) = E_dem(t,1);

           E_dem_out(t,1) = 0;

       else

           E_dem_in(t,1) = 0;

           E_dem_out(t,1) = abs(E_dem(t,1));

       end

end

storage_count = 1; %number of simulated storages

% Speicher 1

s(1).name = 'PHES'; %Pumped hydro energy storage

s(1).eta = 0.8; %roundtrip efficiency

s(1).E_0 = 0; %initially stored energy

s(1).sr_OM = 0.5; %variable operation costs depending on the discharged MWh

% Definition der Werte für Einspeichermodul

s(1).mod(1).name = 'charging power unit'; 

s(1).mod(1).cap = 500; %charging power capacity [MW]

s(1).mod(1).sc = 10; %specific costs for the installation of 1MW of charging power capacity

% Definition der Werte für Einspeichermodul

s(1).mod(2).name = 'discharging power unit'; 

s(1).mod(2).cap = 500; %discharging power capacity [MW]

s(1).mod(2).sc = 10; %specific costs for the installation of 1MW of discharging power capacity

% Definition der Werte für Speichereinheit

s(1).mod(3).name = 'storage unit';

s(1).mod(3).cap = 5000; %storage capacity [MWh]

s(1).mod(3).sc = 2; %specific costs for the installation of 1MWh of storage capacity

%Calculate incoming and outgoing efficiencies and investment costs

for n=1:storage_count

    

    s(n).eta_in = s(n).eta^0.5; 

    s(n).eta_out = s(n).eta^0.5;

    

    s(n).C_inv = 0;

    for m = 1:3

        s(n).C_inv = s(n).C_inv + s(n).mod(m).cap * s(n).mod(m).sc; 

    end

    

end

%storage variables

for n=1:storage_count

   

    %storageprocess variables

    s(n).E = sdpvar(T,1); %storage level of each timestep 

    s(n).E_ex_in = sdpvar(T,1); %incoming energy from the grid

    s(n).E_ex_out = sdpvar(T,1); %outgoing energy to the grid

    s(n).E_in = sdpvar(T,1); %incoming energy of each timestep

    s(n).E_out = sdpvar(T,1); %outgoing energy of each timestep 

    s(n).C_OM = sdpvar(T,1); %operation and maintenance costs

end

%Systemvariables

E_dem_in_rest = sdpvar(T,1);

E_dem_out_rest = sdpvar(T,1);

CE_ex_in = sdpvar(T,1); %Was haben die Speicher in dem Zeitschritt insgesamt aufgenommen

CE_ex_out = sdpvar(T,1); %Was haben die Speicher in dem Zeitschritt insgesamt aufgenommen

%initial values (t=1)

CE_ex_in(1,1) = 0;

CE_ex_out(1,1) = 0;

E_dem_in_rest(1,1) = 0;

E_dem_out_rest(1,1) = 0;

C_OM = 0;

CE_ex_out_ges = 0;

%initial values (t=1) of the storage

for n=1:storage_count

    

s(n).E(1,1) = s(n).E_0;

s(n).E_ex_in(1,1) = 0;

s(n).E_ex_out(1,1) = 0;

s(n).E_in(1,1) = 0;

s(n).E_out(1,1) = 0;

s(n).C_OM = [0];

end

for t=2:T

       

    % initial values for each timestep

    CE_ex_in(t,1) = 0;

    CE_ex_out(t,1) = 0;

    

    for n=1:storage_count 

        

        %non-negativity constraints

        constraints = [constraints, s(n).E(t,1) >= 0]; %Storage level can't be negativ

        constraints = [constraints, s(n).E_ex_in(t,1) >= 0]; %the absorbed energy from the grid can't be negativ

        constraints = [constraints, s(n).E_ex_out(t,1) >= 0]; %the absorbed energy by the grid can't be negativ

        

        %capacity constraints

        constraints = [constraints, s(n).E(t,1) <= s(n).mod(3).cap]; %storage level can't exceed storage capacity

        constraints = [constraints, s(n).E_ex_in(t,1) <= s(n).mod(1).cap]; %the absorbed energy arriving from the grid can't exceed the power capacity of the storage system

        constraints = [constraints, s(n).E_ex_out(t,1) <= s(n).mod(2).cap]; %the energy arriving at the grid can't exceed the power capacity of the storage system

        constraints = [constraints, s(n).E_out(t,1) <= s(n).E(t-1,1)]; %the outgoing can't exceed the stored energy

        constraints = [constraints, s(n).E_in(t,1) <= s(n).mod(3).cap - s(n).E(t-1,1)]; %the incoming energyof each timestep can't exceed the free capacity left

        

        %Storage level calculation

        constraints = [constraints, s(n).E(t,1) == s(n).E(t-1,1) + s(n).E_in(t,1) - s(n).E_out(t,1)]; %stored energy of each timestep is storage level of previous timestep plus/minus in/-outgoing energy

        constraints = [constraints, s(n).E_in(t,1) == s(n).E_ex_in(t,1) * s(n).eta_in]; %incoming energy is the energy arriving from the grid adjusted by efficiency

        constraints = [constraints, s(n).E_out(t,1) == s(n).E_ex_out(t,1) * (1/s(n).eta_out - s(n).sr_OM)]; %outgoing energy is the energy arriving at the grid adjusted by efficiency

        

        %Operations & Maintenance costs

        s(n).C_OM(t,1) = p_f * s(n).sr_OM * s(n).E_ex_out(t,1);

        

        CE_ex_in(t,1) = CE_ex_in(t,1) + s(n).E_ex_in(t,1);

        CE_ex_out(t,1) = CE_ex_out(t,1) + s(n).E_ex_out(t,1);

        

    end

    

constraints = [constraints, CE_ex_in(t,1) <= E_dem_in(t,1)];

constraints = [constraints, CE_ex_out(t,1) <= E_dem_out(t,1)];

constraints = [constraints, E_dem_in_rest(t,1) >= 0];

constraints = [constraints, E_dem_out_rest(t,1) >= 0];  

constraints = [constraints, E_dem_in_rest(t,1) == E_dem_in(t,1) - CE_ex_in(t,1)];

constraints = [constraints, E_dem_in_rest(t,1) <=0];

CE_ex_out_ges = CE_ex_out_ges + CE_ex_out(t,1);   

end

C_inv = 0;

for n = 1:storage_count

    C_inv = C_inv + s(n).C_inv;

end

        

C_OM = sum(s(n).C_OM);

Energy = CE_ex_out_ges;

Costs = C_inv;% + C_OM;

LCOE = Costs / Energy;

%Find an upper bound on the optimal value. Solve problems until feasible

upper = 1;

b = solvesdp([Costs <= upper*Energy,constraints]);

while b.problem

    upper = upper*2;

    b = solvesdp([Costs <= upper*Energy,constraints]);

end

% We can never do better than 0

lower = 0;

while upper-lower > 0.01

    test = (upper + lower)/2;

    b = solvesdp([Costs <= test*Energy,constraints]);

    if b.problem == 0

        bworks = test;

        upper = test;

    else

        lower = test;

    end

end

% solve again at the optimal b so you can extract solutions

solvesdp([Costs <= bworks*Energy,constraints]);

[Johan Löfberg]

I cannot reproduce your NaNs, but I can see that weird things can happen as your problem isn't guaranteed to be feasible.

You can replace the generation of the initial upper bound by

sol = solvesdp(constraints)
if sol.problem == 1
 error('Dude, it is not even feasible!')
else
 upper = double(Costs)/double(Energy);
end

[Johan Löfberg]

This

"If I get it right,  s(n).C_OM(t,1) is an sdpvar and can't be implemented in the cost function."

I have no idea what you are saying. The objective has to be a (in cplex socp representable or convex quadratic) sdpvar (or empty)

[Jasper Meyer]

Thank you very much for the answer. Which logic makes the problem not guaranteed to be feasible?

[Jasper Meyer]

"This

"If I get it right,  s(n).C_OM(t,1) is an sdpvar and can't be implemented in the cost function."

I have no idea what you are saying. The objective has to be a (in cplex socp representable or convex quadratic) sdpvar (or empty)"

I have read the explanation for NaNs from the Yalmip error code. (http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Extra.NANInModel). It is stated, that you can't assign sdpvar to a double. I thought this might have been the problem (calculating the costs including the OM-part).

[Johan Löfberg]

Obviously it is related to the random demand being too large since if you reduce that it works for any random case

[Jasper Meyer]

I realized, that by mistake, I my code didn't include the OM-Costs in the Cost function and therefor didn't produce the NaN error. Sorry for that.

Costs = C_inv;% + C_OM; needs to be changed to Costs = C_inv + C_OM; . 

[Johan Löfberg]

As you say your self, you cannot sub-assign sdpvars to doubles. Still you do that.

First you declare a double

s(n).C_OM = [0];

and then you try to insert an sdpvar into that double

[Jasper Meyer]

Thanks again. That was the missing correction to make my problem work again.