Why do I get Error "Out of Memory caused by failure in user supplied fitness function" in Problem-Based Optimization

Hello,

I am trying to minimize costs of running a heat pump taking spot market prices as reference, but I receive following error:

"Error using optim.problemdef.OptimizationProblem/solve

Out of memory.

Error in Optimierung_WP_Zeitbegrenzung_v2 (line 102)

[sol,fval,exitflag,output] = solve(problem,initialPoint,'Options',options);

Caused by:

Failure in user-supplied fitness function evaluation. Cannot continue.

Failure in initial user-supplied fitness function evaluation. GA cannot continue."

The household I am considering uses a rooftop PV. Therefore my costs do not only depend on the spot market prices, but also on the amount of electricity I can use from the PV and the electricity injected into the grid. The heat demand is given. The optimization should take place for every minute of the day for a whole year. Simplified assumptions to develop a solid code that can be adjusted later on are: Heat pump has to run 15min when started and needs to rest 15min after switched off, therefore start or stop time is assumed to be at each quarter hour. It can run with max power or 2/3 of max power. To test the code I run it over 3 days (4320 minutes), however, the optimization takes place stepwise in 6 hour setps to reduce calculation amount. This is my code

% Input data
x_Ergebnis = zeros(72,1); % variable to save on/off times of Heat pump after each day (15 minute steps)
y_Ergebnis = zeros(4320,1); % variable to save times of grid injection or grid consumption after each day
P_HP_Ergebnis = zeros(72,1); % variable to save power level decision after each day
Leistung_HP_Ergebnis = zeros(4320,1); % variable to save power level in kW after each day
Zustand_HP_Ergebnis = zeros(4320,1); % variable to save on/off times of heat pump after each day (1 minute steps)
s_Ergebnis = zeros(4320,1); % variable to save storage level after each day
storage_capacity = 13.25; % thermal energy storage capacity 
s_Ergebnis(1)=0.8*storage_capacity; % initial heat SOC
for i = 1:12
    Summenprofil = load('Summenprofil.mat');
    COP = 4; % COP of the heat pump
    P_max = 92/COP; % maximum thermal power of the heat pump in kW
    P_max_house = ones(360,1).*86; % maximum total power consumption in kW
    storage_capacity = 13.25; % thermal energy storage capacity
    c_pv_self = ones(360,1).*10; % self-consumed PV electricity price in ct/kWh
    c_pv_feed = ones(360,1).*7; % fed-in PV electricity price in ct/kWh
    diff = optimexpr(360,1); % variable describing load difference
    Self_Consumption_Netzbezug = optimexpr(360,1); % Self-Consumption variable when Grid consumption
    Self_Consumption_Netzeinspeisung = optimexpr(360,1); % Self-Consumption varible when Grid injection
    Grid_Injection = optimexpr(360,1); % amount of grid injection
    Grid_Consumption = optimexpr(360,1); % amount of grid consumption
    Leistung_gesamt = optimexpr(360,1); % total power measured from the household
    s = optimexpr(360,1); % variable describing SOC of heat storage
    z = optimexpr(360,1); % binary variable that is 1 when there is overproduction of PV
    Zustand_HP = optimexpr(360,1); % binary variable (heat pump switched on/off)
    Leistung_HP = optimexpr(360,1); % determines power of heat pump
    m = i*360;
    l = 1+(i-1)*360;
    q = i*24;
    r = 1+(i-1)*24;
    PV_prod = Summenprofil.Summenprofil11.PV(l:m)./1000; 
    Heat_demand = Summenprofil.Summenprofil11.Wrmebedarf(l:m)./1000;
    House_demand = Summenprofil.Summenprofil11.Haushalt3(l:m)./1000;
    c_grid = Summenprofil.Summenprofil11.Kosten(l:m); % grid electricity price in ct/kWh
    % Create optimization variables
    x = optimvar("x",24,"Type","integer","LowerBound",0,"UpperBound",1); % heat pump switched on/off
    P_HP = optimvar("P_HP",24,"Type","integer", "LowerBound",0,"UpperBound",1); % power level of heat pump
    y = optimvar("y",360,"Type","integer","LowerBound",0,"UpperBound",1); % grid injection (yes or no)
    % Set initial starting point for the solver
    initialPoint.x = ones(size(x));
    initialPoint.P_HP = ones(size(P_HP));
    initialPoint.y = ones(size(y));
    for o = 1:24
        if o == 1
            Zustand_HP(1:15) = x(1);
        else
            Zustand_HP((o-1)*15+1:o*15) = x(o); % Zustand_HP = x but with Zustand_HP describing every minute instead of 15 minute time slots
        end
    end
    for t = 1:24
        if t == 1
            Leistung_HP(1:15) = P_HP(1)*(1/3)*P_max+(2/3)*P_max;
        else
            Leistung_HP((t-1)*15+1:t*15) = P_HP(t)*(1/3)*P_max+(2/3)*P_max;
        end
    end
    % Determine values of the optimization expressions for every minute
    for j = 1:360
        diff(j) = (House_demand(j)/60)+(Zustand_HP(j)*Leistung_HP(j)/60)-PV_prod(j)/60; % determines load difference in every minute
        Grid_Consumption(j) = diff(j); % positive "diff" means Grid Consumption
        Self_Consumption_Netzbezug(j) = PV_prod(j)/60; 
        Self_Consumption_Netzeinspeisung(j) = (House_demand(j)/60)+(Zustand_HP(j)*Leistung_HP(j)/60); 
        Grid_Injection(j) = -diff(j); % negative "diff" means Grid Injection
        z(j) = 1-y(j); % z always opposite of y
        Leistung_gesamt(j) = Zustand_HP(j)*Leistung_HP(j)-PV_prod(j)+House_demand(j); % total power in each minute
        % Calculating the SOC   
        if j == 1
            s(1) = s_Ergebnis(l) + (Zustand_HP(1)*Leistung_HP(1) - Heat_demand(1))/60;
        else
            s(j) = s(j-1) + (Zustand_HP(j)*Leistung_HP(j) - Heat_demand(j))/60; 
        end
    end
    % Create problem
    problem = optimproblem;
    % Define problem objective
    problem.Objective = sum(((c_grid./100).*Grid_Consumption.*y)... % cost of grid consumption
        + ((c_pv_self./100).*(Self_Consumption_Netzbezug.*y+Self_Consumption_Netzeinspeisung.*z))... % cost for Self-Consumption
        - ((c_pv_feed./100).*Grid_Injection.*z));
    problem.Constraints.constraint1 = s >= 0; % Minimum SOC
    problem.Constraints.constraint2 = s <= storage_capacity; % Max SOC
    problem.Constraints.constraint3 = diff.*(1-2.*y) <= 0; % y = 1 when Grid Consumption and no Grid Injection
    problem.Constraints.constraint4 = Leistung_gesamt <= P_max_house; % total household power needs to be below maximum
    % Display problem information
    show(problem);
    options = optimoptions(problem,'Display','iter');
    % Solve problem
    [sol,fval,exitflag,output] = solve(problem,initialPoint,'Options',options);
    % Display results
    sol;
    exitflag;
    fval;
 
 % save the solutions
 x_Ergebnis(r:q) = sol.x(1:24); 
 y_Ergebnis(l:m) = sol.y(1:360);
 P_HP_Ergebnis(r:q) = sol.P_HP(1:24);
 for p = 1:24
   if p == 1
        Zustand_HP_Ergebnis(l:l+14) = sol.x(1);
        Leistung_HP_Ergebnis(l:l+14) = sol.P_HP(1)*(1/3)*P_max+(2/3)*P_max;
   else
        Zustand_HP_Ergebnis((p-1)*15+l:p*15+l-1) = sol.x(p); 
        Leistung_HP_Ergebnis((p-1)*15+l:p*15+l-1) = sol.P_HP(p)*(1/3)*P_max+(2/3)*P_max; 
   end
 end

 for k = 1:360
    if k == 1 && i == 1
     s_Ergebnis(l+k-1) = s_Ergebnis(1) + (Zustand_HP_Ergebnis(k)*Leistung_HP_Ergebnis(k) - Heat_demand(k))/60;    
    else
     s_Ergebnis(l+k-1) = s_Ergebnis(l+k-2) + (Zustand_HP_Ergebnis(l+k-1)*Leistung_HP_Ergebnis(k) - Heat_demand(k))/60; % Speichervariable befüllen
    end
 end
 clearvars -except Summenprofil s x_Ergebnis y_Ergebnis P_HP_Ergebnis s_Ergebnis l m Zustand_HP_Ergebnis Leistung_HP_Ergebnis

 if l == 1
  s_Ergebnis(l) = s_Ergebnis(1);
 else
  s_Ergebnis(m+1) = s_Ergebnis(m); % Set current SOC to starting SOC for next time interval
 end

end

Used solver:

Solving problem using ga.

Single objective optimization:

408 Variable(s)

408 Integer variable(s)

1440 Nonlinear inequality constraint(s)

Options:

CreationFcn: @gacreationuniformint

CrossoverFcn: @crossoverlaplace

SelectionFcn: @selectiontournament

MutationFcn: @mutationpower

I would really appreciate if you could help me. Thank you in advance!

[Cris LaPierre]

There is an error with the code you have shared, preventing us from running the code.

[Bastian]

I had shortened the code a bit, to focus on the problem. Now I entered the full code, it should work now.

[Torsten]

I changed

PV_prod = Summenprofil.PV(l:m)./1000; 
Heat_demand = Summenprofil.Wrmebedarf(l:m)./1000;
House_demand = Summenprofil.Haushalt3(l:m)./1000;
c_grid = Summenprofil.Kosten(l:m); % grid electricity price in ct/kWh

to

PV_prod = Summenprofil.Summenprofil11.PV(l:m)./1000; 
Heat_demand = Summenprofil.Summenprofil11.Wrmebedarf(l:m)./1000;
House_demand = Summenprofil.Summenprofil11.Haushalt3(l:m)./1000;
c_grid = Summenprofil.Summenprofil11.Kosten(l:m); % grid electricity price in ct/kWh

I guess it's what you wanted.

[Bastian]

Thank you for your answer. The mistake you corrected was caused by changing the database necessary to be able to upload the file in here, but unfortunately it does not solve my initial problem. The solver still ends with the error. Any other ideas? Thanks in adance! :)

[Torsten]

I did not intend it to be an answer ; I just wanted to make the code work without syntax errors.

[Cris LaPierre]

How long does it take for you to get the error? I let this code run for over an hour, and it was still going strong.

[Bastian]

I would say 20-30 mins, but I didn't measure the time. My previous code had P_HP as an continuous optimvar with 360 instead of 24 rows, so that the power of the Heat Pump could be changed for every minute. This code worked, although it was going slow. I changed it to the integer because I wanted the power to be constant over 15 minutes when switched on or off. I also thought that it might be easier to run the code because of a reduced number of variables, but this doesn't seem to be the case.

[Torsten]

Did the MATLAB solver used to solve your problem change from linprog to intlinprog after the changes you made ?

[Bastian]

It solves the problem with ga. Before the changes also ga.

I now reduced the timesteps in each optimization iteration from 360 to 96, by switching from optimizing each minute to optimizing in 15 min steps. This also allows me to optimize a whole day (24 hours = 96 * 15 mins). The code runs now.

I guess it was just the amount of variables that lead to this error. Also before some of the constraints couldnt be met, which now can be solved easily.

However, I do not get why this error hasn't occured before, because I thought changing continuous P_HP(360,1) to integer P_HP(24,1) would be easier to solve.

[Torsten]

However, I do not get why this error hasn't occured before, because I thought changing continuous P_HP(360,1) to integer P_HP(24,1) would be easier to solve.

No, integer-constraint problems are far more difficult to solve than problems with continuous variables. That's why I asked about the solver MATLAB used before. I didn't look at your system in detail, but you should try to see whether it's possible to formulate your objective function and constraints linearly. It will save a lot of time because MATLAB will be able to either use linprog or intlinprog to solve.

[Bastian]

I don't see any possibility to do so, especially since I use some constraints to generate some kind of "if-conditions". Eg. I use optimvar y only to determine when variable "diff" is positive or negative.

As I said, my code worked. I tried to extend it with a condition so that the heat pump, when switched on (x=1), has to run for at least 2 time steps (30 mins), but I can't figure out a way to do so, without getting a nonsense solution or error. Do you have any idea to do that, without taking optimvars or optimexpr into "if-conditions"?