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

changing resultoptimization

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

Best Answer

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

Related Solutions

MATLAB: Error using .* : Conversion to optim.problemdef.OptimizationExpression from optim.problemdef.OptimizationInequality is not possible

P_bat.*(P_bat<=0)

A relationship operator between two optimization expressions returns a constraint expression. Constraint expressions do not return binary values and cannot be used arithmetically such as by multiplying them by an optimization expression.

Unfortunately for you, min() and max() are not supported on optimization expressions either.

You are going to need to find a different approach.

MATLAB: Too much time taking for output

I might have messed your code up some, but you will get the idea. Your biggest problem was you were not preallocating. Several for-loops can be eliminated. If you want people to help you, you must format your code properly.

%Modeling of PV system using MATLAB
filename = 'book1.xls';
sheetname = 'data';
I_L=xlsread(filename, sheetname,'D2:D8761'); %Load Current (A)
T= xlsread(filename, sheetname, 'A2:A8761'); %Ambient Temp (°C )
S= xlsread(filename, sheetname, 'B2:B8761'); %Solar Radiation(W/m²)
NOCT=43.6;
Alpha=0.068; %Temperature Coefficient
IPV_M=6.89; %Current at MPP
T_cell=T+S*(NOCT-20)/800; %Cell Temp
I_PV=S/1000+Alpha*(T_cell-25);
I_PV(I_PV<0)=0;
%(Routine 2):Initialization
PV_eff=0.16; %PV module efficiency
Wire_eff=0.98; %Wire efficiency
INV_eff=0.95; %inverter efficiency
V_sys=230; %System voltage (V)
PSH=(mean(S)*12)/1000;
V_B=12; %Battery voltage
DCharge_eff=0.8; %Battery charging efficiency
DOD=0.8; %Allowed Depth of charge (%)
Alpha=0.068;
PV_WP=120;
%(Routine3):Intuitive method to find the initial N of PV modules needed
P_L=I_L*230;
A_PV=1.408*0.56;
E_L=sum(P_L)/length(I_L)*24;
N_PV = ceil(E_L/(PV_eff*INV_eff*Wire_eff*PSH*1000*A_PV));
N_PVmin=ceil(N_PV/3);
N_PVmax=5*N_PV;
nh=2;
C_battery=(E_L*nh)/(DOD*DCharge_eff);
C_batterymin=C_battery/3;
C_batterymax=5*C_battery;
%(Routine 4):Find LLP at each Battery Capacity and PV module trial
SOC=zeros(1,length(I_L));
I_Load=zeros(1,length(I_L));
I_Charge=zeros(1,length(I_L));
I_Discharge=zeros(1,length(I_L));
I_Deficit=zeros(1,length(I_L));
I_Damp=zeros(1,length(I_L));
I_Battery=zeros(1,length(I_L));
xxx=N_PVmin:N_PVmax;
yyy=C_batterymin:500:C_batterymax;
C_batteryf=zeros(length(xxx),length(yyy));
C_PVf=zeros(length(xxx),length(yyy));
LLP_calculated=zeros(length(xxx),length(yyy));
t=1;
Vmp=17.4;
NOCT=43.6;
K=0.8;
D=1e-5;
ns=6; %Number of cells per battery
SOC1=1; %Max. battery SOC=100% of total capacity
SOC2=SOC1;
SOC3=0.2; %Min. battery SOC=20% of total capacity
x=1;
I_net=I_PV-I_L;
for m=xxx %Number of PV modules
    y=1;
    for n=yyy %capacity of battery
        SOCmax=n; %Battery Capacity (Wh)
        SOCmin=SOCmax.*(1-DOD);
        w=0;
        for k=1:length(I_L) %Number of hours
            %(Routine 5):In Case of the battery is not empty
            if SOCmax>0
                %(Routine 5.1):In Case of I_PV=I_Load
                if I_net(k)==0
                    I_Load(k)=I_PV(k);
                    if k==1
                        SOC(k)=SOC1;
                    elseif w==0 %In case of battery is not discharged on the previous stage
                        SOC(k)=SOC1;
                    elseif w==1%In case of battery is discharged on the previous stage
                        SOC(k)=SOC(k-1);
                    end
                    %(Routine 5.2):In Case of I_PV>I_Load
                elseif I_net(k)>0
                    I_Load(k)=I_L(k);
                    if k==1
                        SOC(k)=SOC1;
                        I_Damp(k)=I_net(k);
                    elseif k>1
                        if w==0
                            SOC(k)=SOC1;
                            I_Damp(k)=I_net(k);
                        elseif w==1
                            if SOC(k-1)>=SOC1
                                SOC(k)=SOC(k-1);
                                I_Damp(k)=I_net(k);
                            elseif SOC(k-1)<SOC1
                                I_Charge(k)=I_net(k);
                                %%Charging mode
                                B=SOC2;
                                V1=(2+0.148*B)*ns;
                                R1=((0.758+(0.1309/(1.06-B)))*ns)/SOCmax;
                                %syms v;


                                ee=K*V1*I_net(k)-D*SOC2*SOCmax;%,v,0,t));
                                SOCx=SOC1+ee/SOCmax;
                                SOC2=SOCx;
                                SOC(k)=SOC(k-1)+abs(SOC1-SOC2);
                            end
                        end
                    end
                    %(Routine 5.3):In Case of I_PV<I_Load
                elseif I_net(k)<0
                    if w==0
                        I_Discharge(k)=I_L(k)-I_PV(k);
                        I_Battery(k)=I_Discharge(k);
                        I_Load(k)=I_PV(k)+I_Battery(k);
                        %%Discharging mode
                        B=SOC2;
                        V1=(1.926+0.124*B)*ns;
                        R1=(0.19+(0.1037/(B-0.14)))*(ns/SOCmax);
                        %syms v;
                        ee=K*V1*I_net(k)-D*SOC2*SOCmax;
                        SOCx=SOC1+ee/SOCmax;
                        SOC2=SOCx;
                        SOC(k)=SOC2;
                        w=w+1;
                    elseif w==1
                        if SOC(k-1)>SOC3
                            I_Discharge(k)=I_L(k)-I_PV(k);
                            I_Battery(k)=I_Discharge(k);
                            I_Load(k)=I_L(k);
                            B=SOC2;
                            V1=(1.926+0.124*B)*ns;
                            R1=(0.19+(0.1037/(B-0.14)))*(ns/SOCmax);
                            %syms v;
                            ee=K.*V1.*I_net(k)-D*SOC2*SOCmax;
                            SOCx=SOC1+ee/SOCmax;
                            SOC2=SOCx;
                            SOC(k)=SOC(k-1)-abs(SOC1-SOC2);
                        elseif SOC(k-1)<=SOC3
                            SOC2=SOC3;
                            SOC(k)=SOC3;
                            I_Deficit(k)=I_net(k);
                        end
                    end
                end
                    % In Case of the battery is empty
            else
                % In Case of I_PV=I_Load
                if I_net(n)==0
                    SOC(k)=SOC1;
                    I_Load(n)=I_PV(n);
                    I_Damp(n)=0;
                    I_Charge(n)=0;
                    I_Deficit(n)=0;
                    I_Discharge(n)=0;
                    I_Battery(n)=0;
                    %In Case of I_PV>I_Load
                elseif I_net(n)>0
                    SOC(k)=SOC1;
                    I_Load(n)=I_L(n);
                    I_Damp(n)=I_net(n);
                    I_Charge(n)=0;
                    I_Deficit(n)=0;
                    I_Discharge(n)=0;
                    I_Battery(n)=0;
                    %(Routine 6.3):In Case of I_PV<I_Load
                elseif I_net(n)<0
                    SOC(k)=SOC1;
                    I_Damp(n)=0;
                    I_Charge(n)=0;
                    I_Deficit(n)=I_L(n)-I_PV(n);
                    I_Discharge(n)=0;
                    I_Battery(n)=0;
                    I_Load(n)=I_PV(n);
                end
            end
        end
        E_Excess=sum(I_Damp.*230);
        SOC_per=SOC./SOCmax;
        C_batteryf(x,y)=n; %Battery Capacity for each hour and PV module
        C_PVf(x,y)=m; %Number of PV modules for each hour
        LLP_calculated(x,y)=abs(sum(I_Deficit)/sum(I_L));
        y=y+1;
    end
    x=x+1;
end
cc=1;
LLP_ff=zeros(1,size(LLP_calculated,1)*size(LLP_calculated,1));
C_PV_ff=zeros(1,size(LLP_calculated,1)*size(LLP_calculated,1));
C_battery_ff=zeros(1,size(LLP_calculated,1)*size(LLP_calculated,1));
for ii=1:size(LLP_calculated,1)
    for jj=1:size(LLP_calculated,1)
        if LLP_calculated(ii,jj)>=0.0095 && LLP_calculated(ii,jj)<=0.0105
            LLP_ff(cc)=LLP_calculated(ii,jj);
            C_PV_ff(cc)=C_PVf(ii,jj);
            C_battery_ff(cc)=C_batteryf(ii,jj);
            cc=cc+1;
        end
    end
end
%%Initialization of the cost Function
%PV
CC_PV=456; %Capital Cost for one PV
MC_PV=6.5; %Maintenance cost of one PV module per year
Ls=25; %The duration of operation of the system in years
L_PV=25; %The total lifetime period for PV array
%Battery
Ca_battery=1200; %Capacity for one battery (Wh)
CC_batwh=4.8; %Capital Cost for Wh
CC_bat=CC_batwh*Ca_battery; %Capital Cost for one battery
MC_bat=3.4; %Maintenance cost of one battery per year
L_bat=5; %The total lifetime period for battery
Y_bat=(Ls/L_bat)-1; %The expected numbers of the storage battery replacement during the system lifetime
B_rep=50; %Replacement cost for one battery
CC_B_rep=Y_bat*B_rep; %Cost of the storage battery replacement during the system lifetime
%Charge Controller
CC_cc=400; %Capital Cost for one Charge Controller
MC_cc=0; %Maintenance cost of one Charge Controller per year
L_cc=25; %The total lifetime period for Charge Controller
Y_cc=(Ls/L_cc)-1; %The expected numbers of the Charge
%Controller replacement during the system lifetime
N_cc=1; %Number of Charge Controllers during the system lifetime
%Inverter
CC_inv=800; %Capital Cost for one Inverter
MC_inv=0; %Maintenance cost of one Inverter per year
L_inv=25; %The total lifetime period for Inverter
Y_inv=(Ls/L_inv)-1; %The expected numbers of the Inverter replacement during the system lifetime
N_inv=1; %Number of inverters during the system lifetime
%Other Costs
%Circuit Breaker
N_CB=4; %Number of Circuit breaker
C_CB=25; %Cost for one circuit breaker
CC_CB=N_CB*C_CB; %Cost for all circuit breakers
%Support Structure
CC_SS=200;
%Civil work
CC_CW=400;
%Total cost for the Other Costs
CC_OC=CC_CB+CC_SS+CC_CW+CC_B_rep;
ir=0.035; %Real Interest Rate
fr=0.015; %Inflation Rate
ndr=((1+ir)/(1+fr))-1; %Net of discount?inflation rate
% Total life cycle cost calculation
LCCx=zeros(1,length(C_PV_ff));
CC_D=zeros(1,length(C_PV_ff));
LCC=zeros(1,length(C_PV_ff));
for kk=1:length(C_PV_ff)
    LCCx(kk)=(CC_OC/((((1+ndr)^Ls)-1)/((ndr*((1+ndr)^Ls)))))+((C_PV_ff(kk)*(CC_PV+Ls*MC_PV))/L_PV)+(((ceil(C_battery_ff(kk)/Ca_battery)*CC_bat*(1+Y_bat))+(MC_bat*(LsY_bat)))/L_bat)+(((N_cc*CC_cc*(1+Y_cc))+(MC_cc*(Ls-Y_cc)))/L_cc)+(((N_inv*CC_inv*(1+Y_inv))+(MC_inv*(Ls-Y_inv))/L_inv));
    CC_D(kk)=(C_PV_ff(kk)*(CC_PV))+((ceil(C_battery_ff( kk ) / Ca_battery ) * CC_bat * ( 1 + Y_rbat ) ) ) +(N_cc*CC_cc*(1+Y_cc))+(N_inv*CC_inv*(1+Y_inv))+CC_CB+CC_SS;
    LCC(kk) = LCCx(kk)-(((0.13*(CC_D(kk))))/((((1+ndr)^Ls)-1)/((ndr*((1+nd)^Ls)))));
end
%The optimum sizes of PV and battery combination
[MM,II]=min(LCC);
MM;
C_PV_best= C_PV_ff(II);
C_battery_best= C_battery_ff(II);
fprintf('Best PV Moduls Number is: %d\n',C_PV_best);
fprintf('Best Battery Capacity is: %d\n',C_battery_best);

Best Answer

Related Solutions

MATLAB: Error using .* : Conversion to optim.prob​lemdef.Opt​imizationE​xpression from optim.prob​lemdef.Opt​imizationI​nequality is not possible

MATLAB: Too much time taking for output

Related Question

MATLAB: Error using .* : Conversion to optim.problemdef.OptimizationExpression from optim.problemdef.OptimizationInequality is not possible