Optimization result changes as soon as parts of the solution are inserted into the optimization

Hi Matlab-Team,

I am currently working on an optimization of a PV-battery system. The result of the optimization is a 24hr dispatch strategy for the battery. Input of the optimization are load and pv profiles as well as the system characteristics.

I want to simulate the operation of the program and, therefore, rerun the optimization every hour to handle deviations between forecast and occured data. To build the code without mistakes, I started with the easiest case, where the forecast data equal the occured data. The operation schedule of the battery will be an input of the repeating optimization so that the optimization will consider the charge/discharge made by the battery during the lapsed time.

Here starts the problem: Because there is no deviation between forecast and actual data (because they are equal in this case), the whole input doesn't change. The only thing which changes is the additional dispatch strategy of the battery, which equals the optimization results of the previous optimization runs.

When I am thinking logically, the result of the optimizatin should be equal through every time step bc the input is not changing. But in my case it is.

I added the code and a picture (y-axis = State of charge battery, x-axis = time steps) of the result to show the problem.

function [solution] = demand_limit_optimization_with_thresholdSOC...
    (load_profile, pv_generation, C_bat, SOC_0, P_batMin, P_batMax, SOC_min, SOC_max, ...
    eta_bat, threshold, hour, SOC_trajectory_day, dischargeLimitation, chargeLimitation)
%%Optimization for demand limit minimization

n = length(pv_generation);

prob = optimproblem('ObjectiveSense','minimize');

%Define variables, length of array (n) and upper/lower bounds
P_system = optimvar('P_system', n); %P_PV + P_bat
P_bat = optimvar('P_bat',n,'LowerBound',P_batMin,'UpperBound', P_batMax); %Max and Min power of the battery
SOC = optimvar('SOC', n,'LowerBound',SOC_min,'UpperBound',SOC_max); %SOC level
P_PV = optimvar('P_PV', n); %Forecasted PV generation
P_dmd = optimvar('P_dmd', n); %Forecasted load
threshold = optimvar('threshold', n, 'UpperBound', threshold); %added load threshold which cant be crossed

%%Define objective function
f = 0;

for t=1:n
    if (pv_generation(t) < load_profile(t)) && (load_profile(t) >= 0)
        f = f + (P_dmd(t) - P_system(t)); %P_dmd - (P_PV - P_PVcharge + discharge)
    else
        f = 0;
    end
end

prob.Objective = f;

%%Define constraints
P_system_calculate = optimconstr(n,1); %Constraint to calculate the Power of the PV + Battery system
SOC_level = optimconstr(n,1); % Constraint to calculate the SOC level
initialCharge = optimconstr(n,1); %Constraint to define the inital SOC level
threshold_limitation = optimconstr(n,1);
PV_feed = optimconstr(n,1); %Constraint to implement input PV
load_profile_feed = optimconstr(n,1); %Constraint to implement input load
pvDischarge_limitation  = optimconstr(n,1); %Constraint to battery discharge into grid
pvCharge_limitation = optimconstr(n,1);%Constraint to battery charge from grid
SOC_trajectory = optimconstr(n,1);%Constraint to assure, the optimization includes already made steps

%Time depending constraints (which need to be calculated in every time
%step)
for t=1:n
    %1: Calculating the energy which flows from/into the PV-battery system
    %P_system(t) = P_PV(t) + P_bat(t)
    P_system_calculate(t) = P_PV(t) + P_bat(t) == P_system(t);

      %2: Calculation of the time depending charging level SOC
      if t>1
          SOC_level(t) = SOC(t) == SOC(t-1) - ((P_bat(t) / 4) * eta_bat) / (C_bat / 4); %Devided by 4 bc 15 minutes
      else
          initialCharge = SOC(1) == SOC_0; %Initial charge level of the batter: SOC(1) = SOC_0;
      end

      %Threshold limitation
      threshold_limitation(t) = threshold(t) >= P_dmd(t) - P_system(t);

      %3&4: Setting the input data of the generation and the load profiles
      PV_feed(t) = P_PV(t) == pv_generation(t);
      load_profile_feed(t) = P_dmd(t) == load_profile(t);

      %Constraint to forbid feed in charge
      pvDischarge_limitation(t) = P_bat(t) <= P_dmd(t) - P_PV(t);
      pvCharge_limitation(t) = - P_bat(t) <= P_PV(t);
  end

for i=1:hour
    SOC_trajectory((i-1)*4 + 1) = SOC((i-1)*4 + 1) == SOC_trajectory_day((i-1)*4 + 1);
    SOC_trajectory((i-1)*4 + 2) = SOC((i-1)*4 + 2) == SOC_trajectory_day((i-1)*4 + 2);
    SOC_trajectory((i-1)*4 + 3) = SOC((i-1)*4 + 3) == SOC_trajectory_day((i-1)*4 + 3);
    SOC_trajectory((i-1)*4 + 4) = SOC((i-1)*4 + 4) == SOC_trajectory_day((i-1)*4 + 4);
end

%Add constraints to problem
prob.Constraints.P_system_calculate = P_system_calculate;
prob.Constraints.SOC_level = SOC_level;
prob.Constraints.PV_feed = PV_feed;
prob.Constraints.load_profile_feed = load_profile_feed;
prob.Constraints.initialCharge = initialCharge;
prob.Constraints.threshold_limitation = threshold_limitation;
prob.Constraints.SOC_trajectory = SOC_trajectory;

if dischargeLimitation
    prob.Constraints.pvDischarge_limitation = pvDischarge_limitation;
end

if chargeLimitation
    prob.Constraints.pvCharge_limitation = pvCharge_limitation;
end

%%solve optimization problem
solution = solve(prob);
%showproblem(prob) %See full optimization problem with all constraints and
%boundaries

%

[Matt J]

rerun the optimization every hour to handle deviations between forecast and occured data.

If you are re-running the optimization, shouldn't you be calling solve() within some kind of loop? That doesn't seem to be the case in your code.

Regardless, it doesn't appear that it is a Matlab issue, but rather an issue with the problem data that you have provided to Matlab. You question why you are getting a particular solution, but it is the problem data that you have provided which determines that.

[Uwe]

This optimization is just a function which will be repeated within a script. After the optimization, the whole code reacts if deviations exist and feeds this information in the optimization an re-runs it again to find a new optimal solution. There is a bunch of code behind all this.

Nevertheless, the data doesn't change in the second/third/any run. The only thing, which changes, is the input data of the already made steps of the battery (here called SOC_trajectory_day). But this new input equals the result of the first optimization and therefore follows the "recommended" path.

I added a new code which can be runned without the other parts of the code. The issue is the same though nothing changes.

load_profile = [1,2,2,3,6,4,5,8,3,2,2];
pv_generation = [0,0,0,1,2,4,6,4,3,1,0];
C_bat = 5;
SOC_0 = 0.5;
SOC_trajectory_day = [SOC_0,0,0,0,0,0,0,0,0,0,0];
P_batMin = -3; P_batMax = 3;
SOC_min = 0.1;
SOC_max = 1;
eta_bat = 0.9;
threshold = 6;
dischargeLimitation = false;
chargeLimitation = false;

for j = 1:10
    hour = j;
    optimizationResult = demand_limit_optimization_with_thresholdSOC...
    (load_profile, pv_generation, C_bat, SOC_0, P_batMin, P_batMax, SOC_min, SOC_max, ...
    eta_bat, threshold, hour, SOC_trajectory_day, dischargeLimitation, chargeLimitation);

      SOC_trajectory_day = optimizationResult.SOC;

      plot(SOC_trajectory_day); hold on;
  end

function [solution] = demand_limit_optimization_with_thresholdSOC...
    (load_profile, pv_generation, C_bat, SOC_0, P_batMin, P_batMax, SOC_min, SOC_max, ...
    eta_bat, threshold, hour, SOC_trajectory_day, dischargeLimitation, chargeLimitation)
%%Optimization for demand limit minimization

n = length(pv_generation);

prob = optimproblem('ObjectiveSense','minimize');

%Define variables, length of array (n) and upper/lower bounds
P_system = optimvar('P_system', n); %P_PV + P_bat
P_bat = optimvar('P_bat',n,'LowerBound',P_batMin,'UpperBound', P_batMax); %Max and Min power of the battery
SOC = optimvar('SOC', n,'LowerBound',SOC_min,'UpperBound',SOC_max); %SOC level
P_PV = optimvar('P_PV', n); %Forecasted PV generation
P_dmd = optimvar('P_dmd', n); %Forecasted load
threshold = optimvar('threshold', n, 'UpperBound', threshold); %added load threshold which cant be crossed

%%Define objective function
f = 0;

for t=1:n
    if (pv_generation(t) < load_profile(t)) && (load_profile(t) >= 0)
        f = f + (P_dmd(t) - P_system(t)); %P_dmd - (P_PV - P_PVcharge + discharge)
    else
        f = 0;
    end
end

prob.Objective = f;

%%Define constraints
P_system_calculate = optimconstr(n,1); %Constraint to calculate the Power of the PV + Battery system
SOC_level = optimconstr(n,1); % Constraint to calculate the SOC level
initialCharge = optimconstr(n,1); %Constraint to define the inital SOC level
threshold_limitation = optimconstr(n,1);
PV_feed = optimconstr(n,1); %Constraint to implement input PV
load_profile_feed = optimconstr(n,1); %Constraint to implement input load
pvDischarge_limitation  = optimconstr(n,1); %Constraint to battery discharge into grid
pvCharge_limitation = optimconstr(n,1);%Constraint to battery charge from grid
SOC_trajectory = optimconstr(n,1);%Constraint to assure, the optimization includes already made steps

%Time depending constraints (which need to be calculated in every time
%step)
for t=1:n
    %1: Calculating the energy which flows from/into the PV-battery system
    %P_system(t) = P_PV(t) + P_bat(t)
    P_system_calculate(t) = P_PV(t) + P_bat(t) == P_system(t);

      %2: Calculation of the time depending charging level SOC
      if t>1
          SOC_level(t) = SOC(t) == SOC(t-1) - ((P_bat(t) / 4) * eta_bat) / (C_bat / 4); %Devided by 4 bc 15 minutes
      else
          initialCharge = SOC(1) == SOC_0; %Initial charge level of the batter: SOC(1) = SOC_0;
      end

      %Threshold limitation
      threshold_limitation(t) = threshold(t) >= P_dmd(t) - P_system(t);

      %3&4: Setting the input data of the generation and the load profiles
      PV_feed(t) = P_PV(t) == pv_generation(t);
      load_profile_feed(t) = P_dmd(t) == load_profile(t);

      %Constraint to forbid feed in charge
      pvDischarge_limitation(t) = P_bat(t) <= P_dmd(t) - P_PV(t);
      pvCharge_limitation(t) = - P_bat(t) <= P_PV(t);
  end

for i=1:hour
    SOC_trajectory((i-1) + 1) = SOC((i-1) + 1) == SOC_trajectory_day((i-1) + 1);
end

%Add constraints to problem
prob.Constraints.P_system_calculate = P_system_calculate;
prob.Constraints.SOC_level = SOC_level;
prob.Constraints.PV_feed = PV_feed;
prob.Constraints.load_profile_feed = load_profile_feed;
prob.Constraints.initialCharge = initialCharge;
prob.Constraints.threshold_limitation = threshold_limitation;
prob.Constraints.SOC_trajectory = SOC_trajectory;

if dischargeLimitation
    prob.Constraints.pvDischarge_limitation = pvDischarge_limitation;
end

if chargeLimitation
    prob.Constraints.pvCharge_limitation = pvCharge_limitation;
end

%%solve optimization problem
solution = solve(prob);
%showproblem(prob) %See full optimization problem with all constraints and
%boundaries

end

[Nergiz Soyak]

Bu kodun tamamını yollayabilirmisiniz çalışır halini

Answer 1

[Uwe]

Hey guys,

Thanks again for the help. I talked with my prof about this topic and his answer was, that it is not of bigger importance if the trajectory stays the same as long as all different solutions have an optimal objective function value. Therefore, no changes in the problem formulation/objective function are needed.

Thanks again!

Best Uwe

Answer 2

[Alan Weiss]

It seems to me that your objective function is a sum of terms, (P_dmd(t) - P_system(t)), and these terms can be negative. In other words, your objective function might well not be minimized when all the terms are equal.

Generally, when you want to minimize deviations, you take the minimum sum of squares, not the minimum sum. Or you could take the minimum sum of absolute values of the terms.

Alan Weiss

MATLAB mathematical toolbox documentation

[Uwe]

Hi Alan,

Thanks for your answer. I thought about the same when I was implementing the algorithm in the first place. Problem is, my work is based on the work of someone else and changing the objective function changes the whole approach. Furthermore, square roots and other changes will change the programming from linear to non linear.

I was just thinking that inserting some steps of the optimization shouldn't change the result, just thinking logically. It's like planning an optimal way to the supermarket, taking some steps on this optimal way and then plan again. If nothing on the way changes, the way should be the same.

Best Uwe

[Matt J]

It should be easy to convert what you've done to a form appropriate for lsqlin(). Just doing

S=prob2struct(prob)

will give you all the constraints in the requisite matrix form.

[Uwe]

I read the information about the recommended approach, thanks for the hint. But as said above, it is not possible for me to just change the kind of programming or the objective function. For this approach it is necessary to implement the objective function as described. The objective is not to just minimize the deviation between demand and system. Negative values are allowed.

The issue, after some testing, lies in the existence of several minima within the problem. Therefore, it is possible for the optimization to achieve the value of the objective function with different paths. Which path will be the "optimal" result of the problem can depend on the matlab version, the hardware architecture etc. etc., as well as small differences of the input data, even though they equal the optimization result (as described by you here: https://www.mathworks.com/matlabcentral/answers/329955-fmincon-produces-different-results-on-different-pcs ).

I think an implementation of a penality value could be the solution if, and only if, the forecasted data equal the occuring data and the optimization still wants to change the optimal solution. This idea, though, jeopardizes the idea of the optimization as the framework is different.

[Nergiz Soyak]

ıt runs the full code.can you send it bank very urgent needed

Answer 3

[Matt J]

I think an implementation of a penality value could be the solution.

Or additional constraints... The constraints must be under-specified in some way, if they don't define a unique global minimum.