fmincon variable is start value
Show older comments
Dear everyone,
i am currently trying to use fmincon to find a min for the variabel LCOE and the variabels Cutoff, Threshold, Tank and Cap at which the min is.
I use a for loop so i can use three different start values to check my result.
For all three cases the tank size is always the Start value and never changes, so for start values 2500, the tank result will be 2500 and so on the result of LCOE is about the same for all cases though. And is also what i would suspect, i can make an estimate with a surface plot and some set variabels that were part of the result.
Why does the algorithm work that way, what does it mean regarding my results.
Greetings Isabell
clear all, close all
%Startwerte für die Optimierung
Start_CO = [ 5 15 30]; %Euro/MWh
Start_TH = [15 30 50];%Euro/MWh
Start_Tank = [ 2500 10000 20000];%Kg H2
Start_Cap_EL = [5 10 25]; %MW
%Gesammelte Startwerte zum Übergeben
Start = [Start_CO; Start_TH; Start_Tank; Start_Cap_EL];
%Lower and Upper Bound [Cut off Threshold Tanksize Cap_EL]Cur
LB = [5 0 1000 5];
UB = [50 100 40000 40];
for i = 1:3
rng default % For reproducibility
opts = optimoptions(@fmincon,'Algorithm','sqp','FinDiffRelStep',1e-1);
problem = createOptimProblem('fmincon','objective',...
@main_S1,'x0',Start(:,i),'lb',LB,'ub',UB,'options',opts);
gs = GlobalSearch;
[x,f] = run(gs,problem)
% opts = optimoptions(@fmincon,'Algorithm','sqp','FinDiffRelStep',1e-2);
% problem = createOptimProblem('fmincon','objective',...
% @main_S1,'x0',x,'lb',LB,'ub',UB,'options',opts);
% gs = GlobalSearch;
% [x,f] = run(gs,problem)
%Ausgabe der Outputs von verschiedenen Startwerten
LCOE(i) = f;
Cutoff (i) = x(1);
Threshold(i) = x(2);
Tank(i) = x(3);
Cap(i) = x(4);
end
function [LCOE] = main_S1(input)
CUT_OFF = input(1);
THRESHOLD = input(2);
Tank_max = input(3);
Cap_EL = input(4);
%% SP-Steam Plant, EL -ELectrolyser
EC = importdata('Strompreis2020.txt'); %Strompreis in Euro/MWh
%% Parameters Electrolyser
%Cap_EL = 40; %MW Capacity
Capex_EL = 580 * Cap_EL * 1000; % Euro Investment over all
[OPEX] = Opex_Curve(Cap_EL);
%OPEX=1.71
Opex_fix_EL = OPEX/100 * Capex_EL + 84; %Euro yearly
H2_out= 20*Cap_EL; %kgH2/h Hydrogen productivity
LT_EL = 30; %Jahre Lifetime
WACC_EL= 0.04; % % Weight Average Cost / interest rate
%Sp = 50; %Spacing Linspace
%CUT_OFF = 28;
%CUT_OFF = linspace(24,30,Sp);
%CUT_OFF =28*ones(1,Sp); %€/MWh wenn darunter wird produziert
%CUT_OFF = [24 26 28 30];
%% Parameters Power Plant
Cap_SP = 10; %MW Capacity
Capex_SP = 960 * Cap_SP * 1000; % Euro Capital Investment
Opex_fix_SP = 2.5/100 * Capex_SP; %Euro yearly
H2_in_h= 712; %kgH2/h Hydrogen consumption
LT_SP = 30; %Jahre Lifetime
WACC_SP= 0.073; % % Weight Average Cost / interest rate
%LCOH2 = 2.8; %Euro/kg H2
%THRESHOLD = linspace(30,40,Sp); %Startwert Euro/kg
%THRESHOLD = [30 33 36 40];
%THRESHOLD = 30;
%% Cut Off Check
check_EL=zeros(length(EC),length(CUT_OFF));
%EC_h_EL=zeros(length(EC),Sp);
for i= 1:length(EC)% Stunden im Jahr
if EC(i)< CUT_OFF
check_EL(i) = 1;
% EC_h_EL(i,j)=check_EL(i,j)*EC(i)*Cap_EL; %Kosten pro stunde produziert mit der Kapazität
else
check_EL(i) = 0;
% EC_h_EL(i,j)= 0;
end
end
%% Threshold Check
% Checking EC Price against Threshold Price
check = zeros(length(EC),length(THRESHOLD));
for i= 1:length(EC)
if EC(i)> THRESHOLD
check(i) = 1; %immer 1, wenn produziert, sonst 0
%El_h_SP(i,j)=check(i,j)*Cap_SP; %MW produziert
else
check(i) = 0;
%El_h_SP(i,j)= 0;
end
end
%% Energy Storage,
%gefüllt mit kg H2
%Tank_max = 40000;
%Kosten
CAPEX_Tank = 550*Tank_max; %D.D. Papadias umgerechnet in Euro *0.98, Kosten aus Toms PhD 250Euro/kg
%Initialisierung
Tank=zeros(length(EC),1); %
H_SP=zeros(length(EC),1); % H_SP zählt Wie viel wird die Steamplant genutzt
H_EL=zeros(length(EC),1);% H_EL zählt wie viele Stunden wird die Elektrolyse genutzt
Voll=zeros(length(EC),1);% Voll zählt wie viele Stunden war der Tank voll Elektrolyse konnte nicht
Leer=zeros(length(EC),1);% Leer zählt die Stunden in denen nicht produziert werden kann
Kosten_EL=zeros(length(EC),1); %Kosten für Elektrolyse
Kosten_SP=zeros(length(EC),1);%Kosten für SteamPowerPlant
%Nur für 1 Jahr, weil alle Jahre gleich sind.
for A= 1 : LT_EL % Working Years
for i = 1:length(EC) %Stunden im Jahr
if i== 1 % erste Stunde
if A==1 %im ersten Jahr Tank ist Leer
Tank(i,A) = check_EL(i)*H2_out; %Es kann nur gefüllt werden nicht entnommen
%Wie viele Stunden Produziert EL oder SP
if check_EL(i)< check(i) % In der Ersten Stunde ist der Tank noch nicht gefüllt um das Kraftwerk zu nutzen
H_SP(i,A) = 0;
H_EL(i,A) = 0;
elseif check_EL(i)> check(i) %Elektrolyse läuft
H_SP(i,A) = 0;
H_EL(i,A) = 1;
Kosten_EL(i,A) = EC(i)*Cap_EL;
elseif check_EL(i)==1 && check(i)==1 %Beides könnte Laufen, Aber der Tank ist noch leer
H_SP(i,A) = 0;
H_EL(i,A) = 1;
Kosten_EL(i,A) = EC(i)*Cap_EL;
else %check_EL(i,j)== check(i,j)%Beides nicht
H_SP(i,A) = 0;
H_EL(i,A) = 0;
end
else %1 Stunden in all den anderen Jahren
Tank(1,A) = Tank(length(EC),A-1); %Es kann nur gefüllt werden nicht entnommen
% Wie viele Stunden Produziert EL oder SP
if check_EL(i)< check(i) % In der Ersten Stunde ist der Tank noch nicht gefüllt um das Kraftwerk zu nutzen
H_SP(i,A) = 1;
H_EL(i,A) = 0;
Einnahmen_SP(i,A) = EC(i)*Cap_SP;
elseif check_EL(i)> check(i) %Elektrolyse läuft
H_SP(i,A) = 0;
H_EL(i,A) = 1;
Kosten_EL(i,A) = EC(i)*Cap_EL;
elseif check_EL(i)==1 && check(i)==1 %Beides könnte Laufen, Aber der Tank ist noch leer
H_SP(i,A) = 1;
H_EL(i,A) = 1;
Kosten_EL(i,A) = EC(i)*Cap_EL;
Einnahmen_SP(i,A) = EC(i)*Cap_SP;
else %check_EL(i,j)== check(i,j)%Beides nicht
H_SP(i,A) = 0;
H_EL(i,A) = 0;
end
end
%Alle weiteren Stunden
elseif Tank(i-1,A)+ check_EL(i)*H2_out - check(i)*H2_in_h > Tank_max %Tank voll
Leer(i,A)= 0;
Voll(i,A)= 1;
%Wie viele Stunden Produziert EL oder SP
if check_EL(i)< check(i) % Steamplant läuft
H_SP(i,A) = 1;
H_EL(i,A) = 0;
Einnahmen_SP(i,A) = EC(i)*Cap_SP;
Tank (i,A) = Tank(i-1,A) - check(i)*H2_in_h;
elseif check_EL(i)> check(i) %Elektrolyse läuft nicht da Tank voll
H_SP(i,A) = 0;
H_EL(i,A) = 0;
Tank (i,A) = Tank(i-1,A);
elseif check_EL(i)==1 && check(i)==1 %Beides könnte läuft, aber Tank voll also nur Steamplant
H_SP(i,A) = 1;
H_EL(i,A) = 0;
Einnahmen_SP(i,A) = EC(i)*Cap_SP;
%B(i,j,k)=0;
Tank (i,A) = Tank(i-1,A) - check(i)*H2_in_h;
else %Beides nicht
H_SP(i,A) = 0;
H_EL(i,A) = 0;
Tank (i,A) = Tank(i-1,A);
end
elseif Tank(i-1,A)+ check_EL(i)*H2_out - check(i)*H2_in_h <= 0 %Tank leer
Leer(i,A)= 1;
%Wie viele Stunden Produziert EL oder SP
if check_EL(i)< check(i) %SP läuft nicht, da Tank leer
H_SP(i,A) = 0;
H_EL(i,A) = 0;
Tank (i,A) = Tank(i-1,A);
elseif check_EL(i)> check(i) %Elektrolyse läuft
H_SP(i,A) = 0;
H_EL(i,A) = 1;
Kosten_EL(i,A) = EC(i)*Cap_EL;
Tank (i,A) = Tank(i-1,A)+ check_EL(i)*H2_out;
elseif check_EL(i)==1 && check(i)==1 %Beides könnte läufen, Tank aber Leer also nur EL
H_SP(i,A) = 0;
H_EL(i,A) = 1;
Kosten_EL(i,A) = EC(i)*Cap_EL;
Tank (i,A) = Tank(i-1,A)+ check_EL(i)*H2_out;
%B(i,j)=0;
else %Beides nicht
H_SP(i,A) = 0;
H_EL(i,A) = 0;
Tank (i,A) = Tank(i-1,A);
end
else Tank (i,A) = Tank(i-1,A)+ check_EL(i)*H2_out - check(i)*H2_in_h; %Tank wird gefüllt oder geleert
%Wie viele Stunden Produziert EL oder SP
if check_EL(i)< check(i) %Steamplant läuft
H_SP(i,A) = 1;
H_EL(i,A) = 0;
Einnahmen_SP(i,A) = EC(i)*Cap_SP;
elseif check_EL(i)> check(i) %Elektrolyse läuft
H_SP(i,A) = 0;
H_EL(i,A) = 1;
Kosten_EL(i,A) = EC(i)*Cap_EL;
elseif check_EL(i)==1 && check(i)==1 %Beides läuft
H_SP(i,A) = 1;
H_EL(i,A) = 1;
Einnahmen_SP(i,A) = EC(i)*Cap_SP;
Kosten_EL (i,A) = EC(i)*Cap_EL;
% B(i,j)=1;
else %Beides nicht
H_SP(i,A) = 0;
H_EL(i,A) = 0;
end
end
end
end
% end
Einnahmen_a = sum(Einnahmen_SP);
H_EL_a=sum(H_EL);
H_SP_a=sum(H_SP);
%
% Leer_a = sum(Leer); %Wie viele Stunden im Jahr Leer
% Voll_a = sum(Voll); %Wie viele Stunden im Jahr Voll
%% Levelized Cost of Hydrogen for set Parameters
%H2_prod= sum(check_EL)*H2_out; %kgH2 pro Jahr
H2_prod = H_EL_a*H2_out;
H2_prod = squeeze(H2_prod(1,:,:));
EC_a_EL=sum(Kosten_EL); %Electricity Cost per year for Electrolyser
%EC_a_EL= squeeze(EC_a_EL(1,:,:));
% Nenner_EL = zeros(length(CUT_OFF),Sp,LT_EL);
% Zaehler_EL = zeros(length(CUT_OFF),Sp,LT_EL);
% Nenner_l_EL = zeros(length(CUT_OFF),Sp,LT_EL);
% Zaehler_l_EL = zeros(length(CUT_OFF),Sp,LT_EL);
%Laufendekosten über Lifetime
for t=1:LT_EL % Lifetime
%Berechnung:
Nenner_l_EL(t)= (Opex_fix_EL+EC_a_EL(t))/(1+WACC_EL)^t;
Zaehler_l_EL(t)= H2_prod(t)/(1+WACC_EL)^t;
%Eintragen:
if t == 1
Nenner_EL(t) = Nenner_l_EL(t);
Zaehler_EL(t) = Zaehler_l_EL(t);
%Aufaddieren
else
Nenner_EL(t) = Nenner_EL(t-1) + Nenner_l_EL(t);
Zaehler_EL(t) = Zaehler_EL(t-1) + Zaehler_l_EL(t);
end
end
LC_H2= (Capex_EL + Nenner_EL(LT_EL))/Zaehler_EL(LT_EL); %Euro pro Kg H2
LCOH2=LC_H2;
%% LCOE
%h = sum(check); %hours produced anually for different CutOff prices
h = squeeze(H_SP_a(1,:,:)); %hours produced anually for different CutOff prices and Thresholdprices
H2_in_a= h.*H2_in_h; %Hydrogen consumed annual kg for Different CutOff Price
H2C_a = H2_in_a(1,:,:).*LCOH2; %Cost of hydrogen annual Euro for Different CutOff Price
El_a = h.*Cap_SP; %Electricity produced annual MW for Different CutOff Prices
% Nenner = zeros(LT_SP);
% Zaehler = zeros(length(LT_SP);
% Nenner_l = zeros(length(LT_SP);
% Zaehler_l = zeros(length(LT_SP);
for t=1:LT_SP %Working Years
%Berechnung für die Stunde
Nenner_l(t)= (Opex_fix_SP+H2C_a(t))/(1+WACC_SP)^t;
Zaehler_l(t)= El_a(t)/(1+WACC_SP)^t;
%Eintragen
if t == 1
Nenner(t) = Nenner_l(t);
Zaehler(t) = Zaehler_l(t);
%Aufaddieren
else
Nenner(t) = Nenner(t-1) + Nenner_l(t);
Zaehler(t) = Zaehler(t-1) + Zaehler_l(t);
end
end
LCOE = (Capex_SP + CAPEX_Tank + Nenner(LT_SP))/Zaehler(LT_SP); %Euro pro MWh
end
2 Comments
Torsten
on 28 Jul 2022
Why does the algorithm work that way, what does it mean regarding my results.
I don't understand your question.
The results are a consequence of your model, and since we don't know your model, we can't tell why they come out the way they do.
Isabell Hetmann
on 28 Jul 2022
Accepted Answer
You didn't include everything we need to run the code (e.g., Strompreis2020.txt) however here are some guesses:
(1) Your choice of FinDiffRelStep=1e-1 looks very large and may be giving you poor gradient estimates.
(2) Your function may be piecewise constant as a function of the tank size. In particular, I notice that your function computes binary variables H_SP and H_EL to reach its output. If the function contains discretization operations, it is likely to be piecewise flat in places and therefore have no gradient with whcih to evolve the initial point. Below is my favorite example of that,
opts=optimoptions('fmincon',Display='none');
X0=rand(1,5)*9; %initial points
for i=1:numel(X0)
xstart=X0(i)
xfinal=fmincon(@round,xstart,[],[],[],[],0,10,[],opts)
disp ' '
end
5 Comments
Isabell Hetmann
on 28 Jul 2022
Matt J
on 28 Jul 2022
i did include the file now
But not Opex_Curve...
Isabell Hetmann
on 28 Jul 2022
The following code plots main_S1() as a function of the 3rd argument (tank size) in the near neighborhood of the first start point x0=Start(:,1). It does indeed appear to be piecewise constant. Maybe try ga instead.
x0=Start(:,1);
fun= @(z)main_S1(x0+[0;0;z;0]);
dx=linspace(-0.1,0.1,20);
plot(dx, arrayfun(fun,dx));
xlabel dx; ylabel main_S1; axis padded
Isabell Hetmann
on 1 Aug 2022
More Answers (0)
Categories
Find more on Sparse Matrices in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
