How to change the graph type
Show older comments
I would like to change a type of graph.
Code is;
num_simulations = 10000;
%Common parameters
Discount_Rate_min = 0.06; % assume 6-8%
Discount_Rate_max = 0.08;
Discount_Rate_values = unifrnd(Discount_Rate_min, Discount_Rate_max, [num_simulations, 1]);
Lifetime = 20; % years
Electricity_Cost_mean = 0.255; %EUR/kWh
Electricity_Cost_std = 0.04;
Electricity_Cost_values = normrnd(Electricity_Cost_mean, Electricity_Cost_std, [num_simulations,1]);
Electricity_Cost_values(Electricity_Cost_values < 0.02) = 0.02;
Electricity_Cost_values(Electricity_Cost_values > 0.9) = 0.9;
FLH_min = 1000;
FLH_max = 8760;
FLH = unifrnd(FLH_min, FLH_max, [num_simulations, 1]);
LHV = 33.33; %kWh/kgH2
%SOEC parameters
CAPEX_System_SOEC_mean = 4200; %$/kW
CAPEX_System_SOEC_std = 500;
CAPEX_System_SOEC_values = normrnd(CAPEX_System_SOEC_mean, CAPEX_System_SOEC_std, [num_simulations,1]);
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values < 2800) = 2800;
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values > 5600) = 5600;
CAPEX_Stack_SOEC_values = 0.5*CAPEX_System_SOEC_values; % 50% of CAPEX system
CAPEX_SOEC_values = (CAPEX_System_SOEC_values + CAPEX_Stack_SOEC_values);
OPEX_SOEC_values = 3; % 3% of CAPEX/a
System_Efficiency_SOEC_mean = 0.775;
System_Efficiency_SOEC_std = 0.05;
System_Efficiency_SOEC_values = normrnd(System_Efficiency_SOEC_mean, System_Efficiency_SOEC_std, [num_simulations,1]);
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values < 0.74) = 0.74;
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values > 0.81) = 0.81;
%PEM parameters
CAPEX_System_PEM_mean = 1450; %$/kW
CAPEX_System_PEM_std = 50;
CAPEX_System_PEM_values = normrnd(CAPEX_System_PEM_mean, CAPEX_System_PEM_std, [num_simulations,1]);
CAPEX_System_PEM_values(CAPEX_System_PEM_values < 1100) = 1100;
CAPEX_System_PEM_values(CAPEX_System_PEM_values > 1800) = 1800;
CAPEX_Stack_PEM_values = 0.35*CAPEX_System_PEM_values; % 35% of CAPEX system
CAPEX_PEM_values = (CAPEX_System_PEM_values + CAPEX_Stack_PEM_values);
OPEX_PEM_values = 3;
System_Efficiency_PEM_mean = 0.58;
System_Efficiency_PEM_std = 0.01;
System_Efficiency_PEM_values = normrnd(System_Efficiency_PEM_mean, System_Efficiency_PEM_std, [num_simulations,1]);
System_Efficiency_PEM_values(System_Efficiency_PEM_values < 0.56) = 0.56;
System_Efficiency_PEM_values(System_Efficiency_PEM_values > 0.6) = 0.6;
%AEC parameters
CAPEX_System_AEC_mean = 950; % $/kW
CAPEX_System_AEC_std = 50;
CAPEX_System_AEC_values = normrnd(CAPEX_System_AEC_mean, CAPEX_System_AEC_std, [num_simulations,1]);
CAPEX_System_AEC_values(CAPEX_System_AEC_values < 500) = 500;
CAPEX_System_AEC_values(CAPEX_System_AEC_values > 1400) = 1400;
CAPEX_Stack_AEC_values = 0.35*CAPEX_System_AEC_values; % 35% of CAPEX system
CAPEX_AEC_values = (CAPEX_System_AEC_values + CAPEX_Stack_AEC_values);
OPEX_AEC_values = 3;
System_Efficiency_AEC_mean = 0.665;
System_Efficiency_AEC_std = 0.05;
System_Efficiency_AEC_values = normrnd(System_Efficiency_AEC_mean, System_Efficiency_AEC_std, [num_simulations,1]);
System_Efficiency_AEC_values(System_Efficiency_AEC_values < 0.63) = 0.63;
System_Efficiency_AEC_values(System_Efficiency_AEC_values > 0.7) = 0.7;
% Calculate SOEC LCOH values
term1_S = LHV ./ (System_Efficiency_SOEC_values);
term2_S = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_S = (OPEX_SOEC_values / 100);
term4_S = CAPEX_SOEC_values ./ FLH;
LCOH_SOEC = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + Electricity_Cost_values);
% Calculate PEM LCOH values
term1_P = LHV ./ (System_Efficiency_PEM_values);
term2_P = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_P = (OPEX_PEM_values / 100);
term4_P = CAPEX_PEM_values ./ FLH;
LCOH_PEM = term1_P .* ((term2_P ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_P) .* term4_P + Electricity_Cost_values);
% Calculate AEC LCOH values
term1_A = LHV ./ (System_Efficiency_AEC_values);
term2_A = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_A = (OPEX_AEC_values / 100);
term4_A = CAPEX_AEC_values ./ FLH;
LCOH_AEC = term1_A .* ((term2_A ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_A) .* term4_A + Electricity_Cost_values);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
% Plot the graph
figure;
scatter(FLH_filtered, LCOH_SOEC_filtered, 'o', 'DisplayName', 'SOEC');
hold on;
scatter(FLH_filtered, LCOH_PEM_filtered, 'x', 'DisplayName', 'PEM');
scatter(FLH_filtered, LCOH_AEC_filtered, '+', 'DisplayName', 'AEC');
xlabel('FLH');
ylabel('LCOH');
title('3-Sigma Approach Graph');
legend('Location', 'best');
grid on;
hold off;
This graph is the result of the code but I need the graph like below one.
Please let me know the code for it.
1 Comment
Steven Lord
on 5 Jan 2024
Just a comment, you have a number of blocks of code that look awfully similar aside from using different variable names for the data. For example:
CAPEX_System_SOEC_mean = 4200; %$/kW
CAPEX_System_SOEC_std = 500;
CAPEX_System_SOEC_values = normrnd(CAPEX_System_SOEC_mean, CAPEX_System_SOEC_std, [num_simulations,1]);
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values < 2800) = 2800;
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values > 5600) = 5600;
Rather than duplicating all that code, I'd consider creating a function that accepts the desired mean and std values, generates the data, and returns it to its caller. [I'm not exactly sure how you generate your lower and upper bounds; I suspect it's some number of standard deviations, but if they vary you could always accept the lower and upper bounds as additional inputs.]
function data = generateData(meanvalue, stdvalue, num_simulations)
LB = meanvalue - 3*stdvalue;
UB = meanvalue + 3*stdvalue;
data = normrnd(meanvalue, stdvalue, [num_simulations,1]);
data = min(max(data, LB), UB);
% or
%{
data(data < LB) = LB;
data(data > UB) = UB;
%}
Now instead of duplicating this code for CAPEX_System_SOEC_values, Electricity_Cost_values, System_Efficiency_SOEC_values, etc.you just call generateData once per variable. It would make your code a bit shorter and IMO easier to understand.
Accepted Answer
I believe that this is reasonably close —
num_simulations = 10000;
%Common parameters
Discount_Rate_min = 0.06; % assume 6-8%
Discount_Rate_max = 0.08;
Discount_Rate_values = unifrnd(Discount_Rate_min, Discount_Rate_max, [num_simulations, 1]);
Lifetime = 20; % years
Electricity_Cost_mean = 0.255; %EUR/kWh
Electricity_Cost_std = 0.04;
Electricity_Cost_values = normrnd(Electricity_Cost_mean, Electricity_Cost_std, [num_simulations,1]);
Electricity_Cost_values(Electricity_Cost_values < 0.02) = 0.02;
Electricity_Cost_values(Electricity_Cost_values > 0.9) = 0.9;
FLH_min = 1000;
FLH_max = 8760;
FLH = unifrnd(FLH_min, FLH_max, [num_simulations, 1]);
LHV = 33.33; %kWh/kgH2
%SOEC parameters
CAPEX_System_SOEC_mean = 4200; %$/kW
CAPEX_System_SOEC_std = 500;
CAPEX_System_SOEC_values = normrnd(CAPEX_System_SOEC_mean, CAPEX_System_SOEC_std, [num_simulations,1]);
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values < 2800) = 2800;
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values > 5600) = 5600;
CAPEX_Stack_SOEC_values = 0.5*CAPEX_System_SOEC_values; % 50% of CAPEX system
CAPEX_SOEC_values = (CAPEX_System_SOEC_values + CAPEX_Stack_SOEC_values);
OPEX_SOEC_values = 3; % 3% of CAPEX/a
System_Efficiency_SOEC_mean = 0.775;
System_Efficiency_SOEC_std = 0.05;
System_Efficiency_SOEC_values = normrnd(System_Efficiency_SOEC_mean, System_Efficiency_SOEC_std, [num_simulations,1]);
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values < 0.74) = 0.74;
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values > 0.81) = 0.81;
%PEM parameters
CAPEX_System_PEM_mean = 1450; %$/kW
CAPEX_System_PEM_std = 50;
CAPEX_System_PEM_values = normrnd(CAPEX_System_PEM_mean, CAPEX_System_PEM_std, [num_simulations,1]);
CAPEX_System_PEM_values(CAPEX_System_PEM_values < 1100) = 1100;
CAPEX_System_PEM_values(CAPEX_System_PEM_values > 1800) = 1800;
CAPEX_Stack_PEM_values = 0.35*CAPEX_System_PEM_values; % 35% of CAPEX system
CAPEX_PEM_values = (CAPEX_System_PEM_values + CAPEX_Stack_PEM_values);
OPEX_PEM_values = 3;
System_Efficiency_PEM_mean = 0.58;
System_Efficiency_PEM_std = 0.01;
System_Efficiency_PEM_values = normrnd(System_Efficiency_PEM_mean, System_Efficiency_PEM_std, [num_simulations,1]);
System_Efficiency_PEM_values(System_Efficiency_PEM_values < 0.56) = 0.56;
System_Efficiency_PEM_values(System_Efficiency_PEM_values > 0.6) = 0.6;
%AEC parameters
CAPEX_System_AEC_mean = 950; % $/kW
CAPEX_System_AEC_std = 50;
CAPEX_System_AEC_values = normrnd(CAPEX_System_AEC_mean, CAPEX_System_AEC_std, [num_simulations,1]);
CAPEX_System_AEC_values(CAPEX_System_AEC_values < 500) = 500;
CAPEX_System_AEC_values(CAPEX_System_AEC_values > 1400) = 1400;
CAPEX_Stack_AEC_values = 0.35*CAPEX_System_AEC_values; % 35% of CAPEX system
CAPEX_AEC_values = (CAPEX_System_AEC_values + CAPEX_Stack_AEC_values);
OPEX_AEC_values = 3;
System_Efficiency_AEC_mean = 0.665;
System_Efficiency_AEC_std = 0.05;
System_Efficiency_AEC_values = normrnd(System_Efficiency_AEC_mean, System_Efficiency_AEC_std, [num_simulations,1]);
System_Efficiency_AEC_values(System_Efficiency_AEC_values < 0.63) = 0.63;
System_Efficiency_AEC_values(System_Efficiency_AEC_values > 0.7) = 0.7;
% Calculate SOEC LCOH values
term1_S = LHV ./ (System_Efficiency_SOEC_values);
term2_S = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_S = (OPEX_SOEC_values / 100);
term4_S = CAPEX_SOEC_values ./ FLH;
LCOH_SOEC = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + Electricity_Cost_values);
% Calculate PEM LCOH values
term1_P = LHV ./ (System_Efficiency_PEM_values);
term2_P = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_P = (OPEX_PEM_values / 100);
term4_P = CAPEX_PEM_values ./ FLH;
LCOH_PEM = term1_P .* ((term2_P ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_P) .* term4_P + Electricity_Cost_values);
% Calculate AEC LCOH values
term1_A = LHV ./ (System_Efficiency_AEC_values);
term2_A = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_A = (OPEX_AEC_values / 100);
term4_A = CAPEX_AEC_values ./ FLH;
LCOH_AEC = term1_A .* ((term2_A ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_A) .* term4_A + Electricity_Cost_values);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
[FLH_filtereds,sidx] = sort(FLH_filtered);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
LCOH_SOEC_filteredx = movmax(LCOH_SOEC_filtered(sidx), 100);
LCOH_SOEC_filteredn = movmin(LCOH_SOEC_filtered(sidx), 100);
Px1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredx, 3);
Pn1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredn, 3);
Fx1 = polyval(Px1, FLH_filtereds);
Fn1 = polyval(Pn1, FLH_filtereds);
LCOH_PEM_filteredx = movmax(LCOH_PEM_filtered(sidx), 100);
LCOH_PEM_filteredn = movmin(LCOH_PEM_filtered(sidx), 100);
Px2 = polyfit(FLH_filtereds, LCOH_PEM_filteredx, 3);
Pn2 = polyfit(FLH_filtereds, LCOH_PEM_filteredn, 3);
Fx2 = polyval(Px2, FLH_filtereds);
Fn2 = polyval(Pn2, FLH_filtereds);
LCOH_AEC_filteredx = movmax(LCOH_AEC_filtered(sidx), 100);
LCOH_AEC_filteredn = movmin(LCOH_AEC_filtered(sidx), 100);
Px3 = polyfit(FLH_filtereds, LCOH_AEC_filteredx, 3);
Pn3 = polyfit(FLH_filtereds, LCOH_AEC_filteredn, 3);
Fx3 = polyval(Px3, FLH_filtereds);
Fn3 = polyval(Pn3, FLH_filtereds);
% Plot the graph
figure;
hs{1} = scatter(FLH_filtered, LCOH_SOEC_filtered, 'o', 'DisplayName', 'SOEC');
hold on;
hs{2} = scatter(FLH_filtered, LCOH_PEM_filtered, 'x', 'DisplayName', 'PEM');
hs{3} = scatter(FLH_filtered, LCOH_AEC_filtered, '+', 'DisplayName', 'AEC');
hold off;
xlabel('FLH');
ylabel('LCOH');
title('3-Sigma Approach Graph');
legend([hs{:}], 'Location', 'best');
grid on;
figure;
hp{1} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx1; flip(Fn1)], 'r', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'SOEC');
hold on;
hp{2} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx2; flip(Fn2)], 'g', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'PEM');
hp{3} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx3; flip(Fn3)], 'b', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'AEC');
hold off;
xlabel('FLH');
ylabel('LCOH');
title('3-Sigma Approach Graph');
legend([hp{:}], 'Location', 'best');
grid on;
I am not certain what result you want. This finds a moving maximum and minimum of every group of data (after sorting the independent variable in order to sort the others), and then does a simple 3-degree polynomial regression to provide a smooth upper and lower edge to the patch plots, and then draws the patch plots.
.
7 Comments
Minhee
on 5 Jan 2024
Thank you so much!
But can I add mean value lines on that graph like below graph?
As always, my pleasure!
Yes.
I added one movmean call for each data set, and then added corresponding polyfit and polyval calls to perform the identical polynomial fits as with the others. I noted the movmean calls as ‘% ADDED’ and the polyfit and polyval results are ‘Pm#’ and ‘Fm#’ respectively, where ‘#’ denotes the numbers of the variables (as I have named them) that correspond with the others. I also added them to the legend call. If you do not want them in the legend, change it to:
legend([hp{1:3}], 'Location', 'best');
and the legend entries for the mean value line plots will not appear. (I tested that to be certain.)
Try this —
num_simulations = 10000;
%Common parameters
Discount_Rate_min = 0.06; % assume 6-8%
Discount_Rate_max = 0.08;
Discount_Rate_values = unifrnd(Discount_Rate_min, Discount_Rate_max, [num_simulations, 1]);
Lifetime = 20; % years
Electricity_Cost_mean = 0.255; %EUR/kWh
Electricity_Cost_std = 0.04;
Electricity_Cost_values = normrnd(Electricity_Cost_mean, Electricity_Cost_std, [num_simulations,1]);
Electricity_Cost_values(Electricity_Cost_values < 0.02) = 0.02;
Electricity_Cost_values(Electricity_Cost_values > 0.9) = 0.9;
FLH_min = 1000;
FLH_max = 8760;
FLH = unifrnd(FLH_min, FLH_max, [num_simulations, 1]);
LHV = 33.33; %kWh/kgH2
%SOEC parameters
CAPEX_System_SOEC_mean = 4200; %$/kW
CAPEX_System_SOEC_std = 500;
CAPEX_System_SOEC_values = normrnd(CAPEX_System_SOEC_mean, CAPEX_System_SOEC_std, [num_simulations,1]);
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values < 2800) = 2800;
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values > 5600) = 5600;
CAPEX_Stack_SOEC_values = 0.5*CAPEX_System_SOEC_values; % 50% of CAPEX system
CAPEX_SOEC_values = (CAPEX_System_SOEC_values + CAPEX_Stack_SOEC_values);
OPEX_SOEC_values = 3; % 3% of CAPEX/a
System_Efficiency_SOEC_mean = 0.775;
System_Efficiency_SOEC_std = 0.05;
System_Efficiency_SOEC_values = normrnd(System_Efficiency_SOEC_mean, System_Efficiency_SOEC_std, [num_simulations,1]);
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values < 0.74) = 0.74;
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values > 0.81) = 0.81;
%PEM parameters
CAPEX_System_PEM_mean = 1450; %$/kW
CAPEX_System_PEM_std = 50;
CAPEX_System_PEM_values = normrnd(CAPEX_System_PEM_mean, CAPEX_System_PEM_std, [num_simulations,1]);
CAPEX_System_PEM_values(CAPEX_System_PEM_values < 1100) = 1100;
CAPEX_System_PEM_values(CAPEX_System_PEM_values > 1800) = 1800;
CAPEX_Stack_PEM_values = 0.35*CAPEX_System_PEM_values; % 35% of CAPEX system
CAPEX_PEM_values = (CAPEX_System_PEM_values + CAPEX_Stack_PEM_values);
OPEX_PEM_values = 3;
System_Efficiency_PEM_mean = 0.58;
System_Efficiency_PEM_std = 0.01;
System_Efficiency_PEM_values = normrnd(System_Efficiency_PEM_mean, System_Efficiency_PEM_std, [num_simulations,1]);
System_Efficiency_PEM_values(System_Efficiency_PEM_values < 0.56) = 0.56;
System_Efficiency_PEM_values(System_Efficiency_PEM_values > 0.6) = 0.6;
%AEC parameters
CAPEX_System_AEC_mean = 950; % $/kW
CAPEX_System_AEC_std = 50;
CAPEX_System_AEC_values = normrnd(CAPEX_System_AEC_mean, CAPEX_System_AEC_std, [num_simulations,1]);
CAPEX_System_AEC_values(CAPEX_System_AEC_values < 500) = 500;
CAPEX_System_AEC_values(CAPEX_System_AEC_values > 1400) = 1400;
CAPEX_Stack_AEC_values = 0.35*CAPEX_System_AEC_values; % 35% of CAPEX system
CAPEX_AEC_values = (CAPEX_System_AEC_values + CAPEX_Stack_AEC_values);
OPEX_AEC_values = 3;
System_Efficiency_AEC_mean = 0.665;
System_Efficiency_AEC_std = 0.05;
System_Efficiency_AEC_values = normrnd(System_Efficiency_AEC_mean, System_Efficiency_AEC_std, [num_simulations,1]);
System_Efficiency_AEC_values(System_Efficiency_AEC_values < 0.63) = 0.63;
System_Efficiency_AEC_values(System_Efficiency_AEC_values > 0.7) = 0.7;
% Calculate SOEC LCOH values
term1_S = LHV ./ (System_Efficiency_SOEC_values);
term2_S = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_S = (OPEX_SOEC_values / 100);
term4_S = CAPEX_SOEC_values ./ FLH;
LCOH_SOEC = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + Electricity_Cost_values);
% Calculate PEM LCOH values
term1_P = LHV ./ (System_Efficiency_PEM_values);
term2_P = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_P = (OPEX_PEM_values / 100);
term4_P = CAPEX_PEM_values ./ FLH;
LCOH_PEM = term1_P .* ((term2_P ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_P) .* term4_P + Electricity_Cost_values);
% Calculate AEC LCOH values
term1_A = LHV ./ (System_Efficiency_AEC_values);
term2_A = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_A = (OPEX_AEC_values / 100);
term4_A = CAPEX_AEC_values ./ FLH;
LCOH_AEC = term1_A .* ((term2_A ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_A) .* term4_A + Electricity_Cost_values);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
[FLH_filtereds,sidx] = sort(FLH_filtered);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
LCOH_SOEC_filteredx = movmax(LCOH_SOEC_filtered(sidx), 100);
LCOH_SOEC_filteredn = movmin(LCOH_SOEC_filtered(sidx), 100);
LCOH_SOEC_filteredm = movmean(LCOH_SOEC_filtered(sidx), 100); % ADDED
Px1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredx, 3);
Pn1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredn, 3);
Pm1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredm, 3);
Fx1 = polyval(Px1, FLH_filtereds);
Fn1 = polyval(Pn1, FLH_filtereds);
Fm1 = polyval(Pm1, FLH_filtereds);
LCOH_PEM_filteredx = movmax(LCOH_PEM_filtered(sidx), 100);
LCOH_PEM_filteredn = movmin(LCOH_PEM_filtered(sidx), 100);
LCOH_PEM_filteredm = movmean(LCOH_PEM_filtered(sidx), 100); % ADDED
Px2 = polyfit(FLH_filtereds, LCOH_PEM_filteredx, 3);
Pn2 = polyfit(FLH_filtereds, LCOH_PEM_filteredn, 3);
Pm2 = polyfit(FLH_filtereds, LCOH_PEM_filteredm, 3);
Fx2 = polyval(Px2, FLH_filtereds);
Fn2 = polyval(Pn2, FLH_filtereds);
Fm2 = polyval(Pm2, FLH_filtereds);
LCOH_AEC_filteredx = movmax(LCOH_AEC_filtered(sidx), 100);
LCOH_AEC_filteredn = movmin(LCOH_AEC_filtered(sidx), 100);
LCOH_AEC_filteredm = movmean(LCOH_AEC_filtered(sidx), 100); % ADDED
Px3 = polyfit(FLH_filtereds, LCOH_AEC_filteredx, 3);
Pn3 = polyfit(FLH_filtereds, LCOH_AEC_filteredn, 3);
Pm3 = polyfit(FLH_filtereds, LCOH_AEC_filteredm, 3);
Fx3 = polyval(Px3, FLH_filtereds);
Fn3 = polyval(Pn3, FLH_filtereds);
Fm3 = polyval(Pm3, FLH_filtereds);
% % % Plot the graph
% % figure;
% % hs{1} = scatter(FLH_filtered, LCOH_SOEC_filtered, 'o', 'DisplayName', 'SOEC');
% % hold on;
% % hs{2} = scatter(FLH_filtered, LCOH_PEM_filtered, 'x', 'DisplayName', 'PEM');
% % hs{3} = scatter(FLH_filtered, LCOH_AEC_filtered, '+', 'DisplayName', 'AEC');
% % hold off;
% % xlabel('FLH');
% % ylabel('LCOH');
% % title('3-Sigma Approach Graph');
% % legend([hs{:}], 'Location', 'best');
% % grid on;
figure;
hp{1} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx1; flip(Fn1)], 'r', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'SOEC');
hold on;
hp{4} = plot(FLH_filtereds, Fm1, '-r', 'LineWidth',2, 'DisplayName','\mu(SOEC)');
hp{2} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx2; flip(Fn2)], 'g', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'PEM');
hp{5} = plot(FLH_filtereds, Fm2, '-g', 'LineWidth',2, 'DisplayName','\mu(PEM)');
hp{3} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx3; flip(Fn3)], 'b', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'AEC');
hp{6} = plot(FLH_filtereds, Fm3, '-b', 'LineWidth',2, 'DisplayName','\mu(AEC)');
hold off;
xlabel('FLH');
ylabel('LCOH');
title('3-Sigma Approach Graph');
legend([hp{:}], 'Location', 'best');
grid on;
.
Minhee
on 9 Jan 2024
Sure! No apologies necessary.
This runs (in less than 1 second on my offline computer), however I cannot get it to run here. I get a red pop-up that says ‘Something went wrong ...’
num_simulations = 10000;
%Common parameters
Discount_Rate_min = 0.06; % assume 6-8%
Discount_Rate_max = 0.08;
Discount_Rate_values = unifrnd(Discount_Rate_min, Discount_Rate_max, [num_simulations, 1]);
Lifetime = 20; % years
Electricity_Cost_values = 0.38;
FLH_min = 1000;
FLH_max = 8760;
FLH = unifrnd(FLH_min, FLH_max, [num_simulations, 1]);
LHV = 33.33; %kWh/kgH2
%SOEC parameters
CAPEX_System_SOEC_mean = 750; %$/kW
CAPEX_System_SOEC_std = 50;
CAPEX_System_SOEC_values = normrnd(CAPEX_System_SOEC_mean, CAPEX_System_SOEC_std, [num_simulations,1]);
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values < 500) = 500;
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values > 1000) = 1000;
CAPEX_Stack_SOEC_values = 0.5*CAPEX_System_SOEC_values; % 50% of CAPEX system
CAPEX_SOEC_values = (CAPEX_System_SOEC_values + CAPEX_Stack_SOEC_values);
OPEX_SOEC_values = 3; % 3% of CAPEX/a
System_Efficiency_SOEC_mean = 0.835;
System_Efficiency_SOEC_std = 0.01;
System_Efficiency_SOEC_values = normrnd(System_Efficiency_SOEC_mean, System_Efficiency_SOEC_std, [num_simulations,1]);
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values < 0.77) = 0.77;
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values > 0.9) = 0.9;
%PEM parameters
CAPEX_System_PEM_mean = 550; %$/kW
CAPEX_System_PEM_std = 50;
CAPEX_System_PEM_values = normrnd(CAPEX_System_PEM_mean, CAPEX_System_PEM_std, [num_simulations,1]);
CAPEX_System_PEM_values(CAPEX_System_PEM_values < 200) = 200;
CAPEX_System_PEM_values(CAPEX_System_PEM_values > 900) = 900;
CAPEX_Stack_PEM_values = 0.35*CAPEX_System_PEM_values; % 35% of CAPEX system
CAPEX_PEM_values = (CAPEX_System_PEM_values + CAPEX_Stack_PEM_values);
OPEX_PEM_values = 3;
System_Efficiency_PEM_mean = 0.705;
System_Efficiency_PEM_std = 0.01;
System_Efficiency_PEM_values = normrnd(System_Efficiency_PEM_mean, System_Efficiency_PEM_std, [num_simulations,1]);
System_Efficiency_PEM_values(System_Efficiency_PEM_values < 0.67) = 0.67;
System_Efficiency_PEM_values(System_Efficiency_PEM_values > 0.74) = 0.74;
%AEC parameters
CAPEX_System_AEC_mean = 450; % $/kW
CAPEX_System_AEC_std = 50;
CAPEX_System_AEC_values = normrnd(CAPEX_System_AEC_mean, CAPEX_System_AEC_std, [num_simulations,1]);
CAPEX_System_AEC_values(CAPEX_System_AEC_values < 200) = 200;
CAPEX_System_AEC_values(CAPEX_System_AEC_values > 700) = 700;
CAPEX_Stack_AEC_values = 0.35*CAPEX_System_AEC_values; % 35% of CAPEX system
CAPEX_AEC_values = (CAPEX_System_AEC_values + CAPEX_Stack_AEC_values);
OPEX_AEC_values = 3;
System_Efficiency_AEC_mean = 0.75;
System_Efficiency_AEC_std = 0.01;
System_Efficiency_AEC_values = normrnd(System_Efficiency_AEC_mean, System_Efficiency_AEC_std, [num_simulations,1]);
System_Efficiency_AEC_values(System_Efficiency_AEC_values < 0.7) = 0.7;
System_Efficiency_AEC_values(System_Efficiency_AEC_values > 0.8) = 0.8;
% Calculate SOEC LCOH values
term1_S = LHV ./ (System_Efficiency_SOEC_values);
term2_S = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_S = (OPEX_SOEC_values / 100);
term4_S = CAPEX_SOEC_values ./ FLH;
LCOH_SOEC = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + Electricity_Cost_values);
% Calculate PEM LCOH values
term1_P = LHV ./ (System_Efficiency_PEM_values);
term2_P = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_P = (OPEX_PEM_values / 100);
term4_P = CAPEX_PEM_values ./ FLH;
LCOH_PEM = term1_P .* ((term2_P ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_P) .* term4_P + Electricity_Cost_values);
% Calculate AEC LCOH values
term1_A = LHV ./ (System_Efficiency_AEC_values);
term2_A = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_A = (OPEX_AEC_values / 100);
term4_A = CAPEX_AEC_values ./ FLH;
LCOH_AEC = term1_A .* ((term2_A ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_A) .* term4_A + Electricity_Cost_values);
% % Plotting LCOH vs FLH for SOEC, PEM, and AEC in one graph
% figure;
% % Plotting for SOEC
% plot(FLH, LCOH_SOEC, 'o', 'DisplayName', 'SOEC');
% hold on;
% % Plotting for PEM
% plot(FLH, LCOH_PEM, 'o', 'DisplayName', 'PEM');
% % Plotting for AEC
% plot(FLH, LCOH_AEC, 'o', 'DisplayName', 'AEC');
% title('LCOH vs FLH for Different Technologies');
% xlabel('FLH (hours/year)');
% ylabel('LCOH ($/kgH2)');
% legend('Location', 'Best');
% grid on;
% hold off;
[LCOH_SOECmax,LCOH_SOECmin,LCOH_SOECmean,FLHs] = scatterfilter(FLH,LCOH_SOEC);
[LCOH_PEMmax,LCOH_PEMmin,LCOH_PEMmean,FLHs] = scatterfilter(FLH,LCOH_PEM);
[LCOH_AECmax,LCOH_AECmin,LCOH_AECmean,FLHs] = scatterfilter(FLH,LCOH_AEC);
toc
% % % Plot the graph
% % figure;
% % hs{1} = scatter(FLH_filtered, LCOH_SOEC_filtered, 'o', 'DisplayName', 'SOEC');
% % hold on;
% % hs{2} = scatter(FLH_filtered, LCOH_PEM_filtered, 'x', 'DisplayName', 'PEM');
% % hs{3} = scatter(FLH_filtered, LCOH_AEC_filtered, '+', 'DisplayName', 'AEC');
% % hold off;
% % xlabel('FLH');
% % ylabel('LCOH');
% % title('3-Sigma Approach Graph');
% % legend([hs{:}], 'Location', 'best');
% % grid on;
figure;
hp{1} = patch([FLHs; flip(FLHs)], [LCOH_SOECmin; flip(LCOH_SOECmax)], 'r', 'EdgeColor','none', 'FaceAlpha',0.5, 'DisplayName', 'SOEC');
hold on;
hp{4} = plot(FLHs, LCOH_SOECmean, '-r', 'LineWidth',2, 'DisplayName','\mu(SOEC)');
hp{2} = patch([FLHs; flip(FLHs)], [LCOH_PEMmax; flip(LCOH_PEMmin)], 'g', 'EdgeColor','none', 'FaceAlpha',0.5, 'DisplayName', 'PEM');
hp{5} = plot(FLHs, LCOH_PEMmean, '-g', 'LineWidth',2, 'DisplayName','\mu(PEM)');
hp{3} = patch([FLHs; flip(FLHs)], [LCOH_AECmax; flip(LCOH_AECmin)], 'b', 'EdgeColor','none', 'FaceAlpha',0.5, 'DisplayName', 'AEC');
hp{6} = plot(FLHs, LCOH_AECmean, '-b', 'LineWidth',2, 'DisplayName','\mu(AEC)');
hold off;
xlabel('FLH');
ylabel('LCOH');
title('3-Sigma Approach Graph');
legend([hp{:}], 'Location', 'best');
grid on;
toc
function [Fmax,Fmin,Fmean,xs] = scatterfilter(x,y)
[xs,sidx] = sort(x,'ascend');
ymax = movmax(y(sidx), 100);
ymin = movmin(y(sidx), 100);
ymean = movmean(y(sidx), 100); % ADDED
Px1 = polyfit(xs, ymax, 3);
Pn1 = polyfit(xs, ymin, 3);
Pm1 = polyfit(xs, ymean, 3);
Fmax = polyval(Px1, xs);
Fmin = polyval(Pn1, xs);
Fmean = polyval(Pm1, xs);
end
It produces this result —
EDIT — (9 Jan 2024 at 16:40)
Returned to check it again and without any changes to the code, it now runs. I have no idea what the problem was previously.
.
Minhee
on 23 Jan 2024
I added the loglog plot just to see what iut would look like. You can delete it if you want to.
Try this —
num_simulations = 10000;
%Common parameters
Discount_Rate_min = 0.06; % assume 6-8%
Discount_Rate_max = 0.08;
Discount_Rate_values = unifrnd(Discount_Rate_min, Discount_Rate_max, [num_simulations, 1]);
Lifetime = 20; % years
electricity_cost_SE = [-0.002,0,0.001,0.001,0.002,0.003,0.003,0.004,0.005,0.006,0.007,0.01,0.013,0.017,0.021,0.026,0.031,0.036,0.04,0.046,0.051,0.057,0.063,0.07,0.079,0.091,0.101];
FLH = [0,100,200,300,400,500,600,700,800,900,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000,7500,8000,8500,8760];
LHV = 33.33; %kWh/kgH2
% Create a table
dataTable_SE = table(FLH', electricity_cost_SE', 'VariableNames', {'FLH', 'electricity_cost_SE'});
%SOEC parameters
CAPEX_System_SOEC_mean = 4200; %$/kW
CAPEX_System_SOEC_std = 1142.7;
CAPEX_System_SOEC_values = normrnd(CAPEX_System_SOEC_mean, CAPEX_System_SOEC_std, [num_simulations,1]);
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values < 2800) = 2800;
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values > 5600) = 5600;
%CAPEX_Stack_SOEC_values = 0.5*CAPEX_System_SOEC_values; % 50% of CAPEX system
CAPEX_SOEC_values = CAPEX_System_SOEC_values;
OPEX_SOEC_values = 3; % 3% of CAPEX/a
System_Efficiency_SOEC_mean = 0.775;
System_Efficiency_SOEC_std = 0.0408;
System_Efficiency_SOEC_values = normrnd(System_Efficiency_SOEC_mean, System_Efficiency_SOEC_std, [num_simulations,1]);
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values < 0.74) = 0.74;
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values > 0.81) = 0.81;
%PEM parameters
CAPEX_System_PEM_mean = 1450; %$/kW
CAPEX_System_PEM_std = 367.42;
CAPEX_System_PEM_values = normrnd(CAPEX_System_PEM_mean, CAPEX_System_PEM_std, [num_simulations,1]);
CAPEX_System_PEM_values(CAPEX_System_PEM_values < 1100) = 1100;
CAPEX_System_PEM_values(CAPEX_System_PEM_values > 1800) = 1800;
%CAPEX_Stack_PEM_values = 0.35*CAPEX_System_PEM_values; % 35% of CAPEX system
CAPEX_PEM_values = CAPEX_System_PEM_values;
OPEX_PEM_values = 3;
System_Efficiency_PEM_mean = 0.58;
System_Efficiency_PEM_std = 0.0163;
System_Efficiency_PEM_values = normrnd(System_Efficiency_PEM_mean, System_Efficiency_PEM_std, [num_simulations,1]);
System_Efficiency_PEM_values(System_Efficiency_PEM_values < 0.56) = 0.56;
System_Efficiency_PEM_values(System_Efficiency_PEM_values > 0.6) = 0.6;
%AEC parameters
CAPEX_System_AEC_mean = 950; % $/kW
CAPEX_System_AEC_std = 387.3;
CAPEX_System_AEC_values = normrnd(CAPEX_System_AEC_mean, CAPEX_System_AEC_std, [num_simulations,1]);
CAPEX_System_AEC_values(CAPEX_System_AEC_values < 500) = 500;
CAPEX_System_AEC_values(CAPEX_System_AEC_values > 1400) = 1400;
%CAPEX_Stack_AEC_values = 0.35*CAPEX_System_AEC_values; % 35% of CAPEX system
CAPEX_AEC_values = CAPEX_System_AEC_values;
OPEX_AEC_values = 3;
System_Efficiency_AEC_mean = 0.665;
System_Efficiency_AEC_std = 0.0271;
System_Efficiency_AEC_values = normrnd(System_Efficiency_AEC_mean, System_Efficiency_AEC_std, [num_simulations,1]);
System_Efficiency_AEC_values(System_Efficiency_AEC_values < 0.63) = 0.63;
System_Efficiency_AEC_values(System_Efficiency_AEC_values > 0.7) = 0.7;
% Calculate SOEC LCOH values
term1_S = LHV ./ System_Efficiency_SOEC_values;
term2_S = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_S = (OPEX_SOEC_values / 100);
term4_S = CAPEX_SOEC_values ./ FLH;
LCOH_SOEC = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + electricity_cost_SE);
% Calculate PEM LCOH values
term1_P = LHV ./ System_Efficiency_PEM_values;
term2_P = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_P = (OPEX_PEM_values / 100);
term4_P = CAPEX_PEM_values ./ FLH;
LCOH_PEM = term1_P .* ((term2_P ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_P) .* term4_P + electricity_cost_SE);
% Calculate AEC LCOH values
term1_A = LHV ./ System_Efficiency_AEC_values;
term2_A = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_A = (OPEX_AEC_values / 100);
term4_A = CAPEX_AEC_values ./ FLH;
LCOH_AEC = term1_A .* ((term2_A ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_A) .* term4_A + electricity_cost_SE);
% Calculate mean and standard deviation of Electricity_Cost_values
LCOH_SOEC_mean = mean(LCOH_SOEC);
LCOH_PEM_mean = mean(LCOH_PEM);
LCOH_AEC_mean = mean(LCOH_AEC);
SOEC_std = 1.9;
PEM_std = 0.77;
AEC_std = 0.72;
% Define the 3-sigma range
lower_limit_S = LCOH_SOEC_mean - 3 * SOEC_std;
upper_limit_S = LCOH_SOEC_mean + 3 * SOEC_std;
lower_limit_P = LCOH_PEM_mean - 3 * PEM_std;
upper_limit_P = LCOH_PEM_mean + 3 * PEM_std;
lower_limit_A = LCOH_AEC_mean - 3 * AEC_std;
upper_limit_A = LCOH_AEC_mean + 3 * AEC_std;
mean_vectors = [LCOH_SOEC_mean; LCOH_PEM_mean; LCOH_AEC_mean];
Lvmv = any(isfinite(mean_vectors),1);
mean_vectors = mean_vectors(:,Lvmv);
limit_vectors = [lower_limit_S; upper_limit_S; lower_limit_P; upper_limit_P; lower_limit_A; upper_limit_A];
Lvlv = any(isfinite(limit_vectors),1);
limit_vectors = limit_vectors(:,Lvlv);
FLHpatch = [FLH(Lvmv) flip(FLH(Lvmv))];
DN = ["SOEC", "PEM", "AEC"];
figure
hold on
cm = eye(3);
for k = 1:size(mean_vectors,1)
patchidx = [1 2] + 2*(k-1);
hp = patch(FLHpatch, [limit_vectors(patchidx(1),:) flip(limit_vectors(patchidx(2),:))], cm(k,:), 'EdgeColor','none', 'FaceAlpha',0.5, 'DisplayName',[DN(k)+" Limits"]);
plot(FLH(Lvmv), mean_vectors(k,:), 'LineWidth',1, 'Color',hp.FaceColor, 'DisplayName',[DN(k)+" Mean"])
end
hold off
xlabel('FLH')
ylabel('LCOH')
legend('Location','best')
figure
hold on
cm = eye(3);
for k = 1:size(mean_vectors,1)
patchidx = [1 2] + 2*(k-1);
hp = patch(FLHpatch, [limit_vectors(patchidx(1),:) flip(limit_vectors(patchidx(2),:))], cm(k,:), 'EdgeColor','none', 'FaceAlpha',0.5, 'DisplayName',[DN(k)+" Limits"]);
plot(FLH(Lvmv), mean_vectors(k,:), 'LineWidth',1, 'Color',hp.FaceColor, 'DisplayName',[DN(k)+" Mean"])
end
hold off
xlabel('FLH')
ylabel('LCOH')
legend('Location','best')
Ax = gca;
Ax.XScale = 'log';
Ax.YScale = 'log';
Make appropriate changes to get the result you want.
.
Minhee
on 6 Feb 2024
Thank you for your answers!
Can I plot this graph like your graph?
The code is;
num_simulations = 1000;
%Common parameters
Discount_Rate_min = 0.06; % assume 6-8%
Discount_Rate_max = 0.08;
Discount_Rate_values = unifrnd(Discount_Rate_min, Discount_Rate_max, [num_simulations, 1]);
Lifetime = 20; % years
% Define electricity costs
FLH_DE = [0,100,200,300,400,500,600,700,800,900,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000];
FLH_FI = [0,100,200,300,400,500,600,700,800,900,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000];
FLH_SE = [0,100,200,300,400,500,600,700,800,900,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000];
FLH_NO = [0,100,200,300,400,500,600,700,800,900,1000,1500,2000,2500,3000,3500,4000,4500,5000,5500,6000,6500,7000,7500,8000,8500,8760];
%electricity_cost_DE_LU = [-0.019,-0.001,0.001,0.006,0.011,0.018,0.024,0.03,0.035,0.04,0.044,0.061,0.075,0.087,0.099,0.11,0.12,0.129,0.138,0.147,0.157,0.167,0.178,0.191,0.206,0.223,0.235];
%electricity_cost_FI = [-0.002,0,0.001,0.002,0.004,0.005,0.006,0.007,0.008,0.008,0.009,0.013,0.018,0.025,0.032,0.04,0.048,0.056,0.064,0.073,0.082,0.091,0.101,0.112,0.125,0.141,0.154];
%electricity_cost_SE = [-0.002,0,0.001,0.001,0.002,0.003,0.003,0.004,0.005,0.006,0.007,0.01,0.013,0.017,0.021,0.026,0.031,0.036,0.04,0.046,0.051,0.057,0.063,0.07,0.079,0.091,0.101];
%electricity_cost_NO = [0,0.005,0.012,0.018,0.023,0.027,0.031,0.035,0.039,0.043,0.046,0.057,0.065,0.071,0.077,0.082,0.086,0.091,0.095,0.098,0.102,0.106,0.111,0.118,0.126,0.138,0.146];
electricity_cost_DE_LU_H = [0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.067,0.075,0.083,0.091,0.1,0.108];
electricity_cost_FI_H = [0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.055,0.054,0.055,0.058,0.062,0.067,0.073,0.079,0.086];
electricity_cost_SE_H = [0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.055,0.053,0.052,0.052,0.053,0.055];
electricity_cost_NO_H = [0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.063,0.067,0.07,0.074,0.077,0.081,0.085,0.089,0.095,0.101,0.11,0.116];
electricity_cost_DE_LU_L = [0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.054,0.056,0.06,0.066,0.072,0.079];
electricity_cost_FI_L = [0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.051,0.045,0.042,0.041,0.042,0.045,0.048,0.053,0.057];
electricity_cost_SE_L = [0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.052,0.047,0.043,0.041,0.04,0.041];
electricity_cost_NO_L = [0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.06,0.055,0.056,0.059,0.062,0.065,0.068,0.072,0.075,0.078,0.082,0.085,0.088];
%electricity_cost_DE_LU_P = [0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,-0.019,-0.001,0.001,0.006,0.011,0.018,0.024,0.03,0.035,0.04,0.044,0.061,0.075,0.087];
%electricity_cost_FI_P = [0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,-0.002,0,0.001,0.002,0.004,0.005,0.006,0.007,0.008,0.008,0.009,0.013,0.018,0.025];
%electricity_cost_SE_P = [0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,-0.002,0,0.001,0.001,0.002,0.003,0.003,0.004,0.005,0.006,0.007,0.01,0.013,0.017];
%electricity_cost_NO_P = [0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0,0.005,0.012,0.018,0.023,0.027,0.031,0.035,0.039,0.043,0.046,0.057,0.065,0.071];
% Create a table
%dataTable_DE_LU = table(FLH', electricity_cost_DE_LU', 'VariableNames', {'FLH', 'electricity_cost_DE_LU'});
%dataTable_FI = table(FLH', electricity_cost_FI', 'VariableNames', {'FLH', 'electricity_cost_FI'});
%dataTable_SE = table(FLH', electricity_cost_SE', 'VariableNames', {'FLH', 'electricity_cost_SE'});
%dataTable_NO = table(FLH', electricity_cost_NO', 'VariableNames', {'FLH', 'electricity_cost_NO'});
dataTable_DE_LU_H = table(FLH_DE', electricity_cost_DE_LU_H', 'VariableNames', {'FLH', 'electricity_cost_DE_LU_H'});
dataTable_FI_H = table(FLH_FI', electricity_cost_FI_H', 'VariableNames', {'FLH', 'electricity_cost_FI_H'});
dataTable_SE_H = table(FLH_SE', electricity_cost_SE_H', 'VariableNames', {'FLH', 'electricity_cost_SE_H'});
dataTable_NO_H = table(FLH_NO', electricity_cost_NO_H', 'VariableNames', {'FLH', 'electricity_cost_NO_H'});
dataTable_DE_LU_L = table(FLH_DE', electricity_cost_DE_LU_L', 'VariableNames', {'FLH', 'electricity_cost_DE_LU_L'});
dataTable_FI_L = table(FLH_FI', electricity_cost_FI_L', 'VariableNames', {'FLH', 'electricity_cost_FI_L'});
dataTable_SE_L = table(FLH_SE', electricity_cost_SE_L', 'VariableNames', {'FLH', 'electricity_cost_SE_L'});
dataTable_NO_L = table(FLH_NO', electricity_cost_NO_L', 'VariableNames', {'FLH', 'electricity_cost_NO_L'});
%dataTable_DE_LU_P = table(FLH', electricity_cost_DE_LU_P', 'VariableNames', {'FLH', 'electricity_cost_DE_LU_P'});
%dataTable_FI_P = table(FLH', electricity_cost_FI_P', 'VariableNames', {'FLH', 'electricity_cost_FI_P'});
%dataTable_SE_P = table(FLH', electricity_cost_SE_P', 'VariableNames', {'FLH', 'electricity_cost_SE_P'});
%dataTable_NO_P = table(FLH', electricity_cost_NO_P', 'VariableNames', {'FLH', 'electricity_cost_NO_P'});
LHV = 33.33; %kWh/kgH2
%AEC parameters
CAPEX_System_AEC_mean = 950; % $/kW
CAPEX_System_AEC_std = 387.3;
CAPEX_System_AEC_values = normrnd(CAPEX_System_AEC_mean, CAPEX_System_AEC_std, [num_simulations,1]);
CAPEX_System_AEC_values(CAPEX_System_AEC_values < 500) = 500;
CAPEX_System_AEC_values(CAPEX_System_AEC_values > 1400) = 1400;
%CAPEX_Stack_AEC_values = 0.35*CAPEX_System_AEC_values; % 35% of CAPEX system
CAPEX_AEC_values = CAPEX_System_AEC_values;
OPEX_AEC_values = 3;
System_Efficiency_AEC_mean = 0.665;
System_Efficiency_AEC_std = 0.0271;
System_Efficiency_AEC_values = normrnd(System_Efficiency_AEC_mean, System_Efficiency_AEC_std, [num_simulations,1]);
System_Efficiency_AEC_values(System_Efficiency_AEC_values < 0.63) = 0.63;
System_Efficiency_AEC_values(System_Efficiency_AEC_values > 0.7) = 0.7;
% Calculate AEC LCOH values
term1_S = LHV ./ (System_Efficiency_AEC_values);
term2_S = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_S = (OPEX_AEC_values / 100);
term4_DE = CAPEX_AEC_values ./ FLH_DE;
term4_FI = CAPEX_AEC_values ./ FLH_FI;
term4_SE = CAPEX_AEC_values ./ FLH_SE;
term4_NO = CAPEX_AEC_values ./ FLH_NO;
%DE_LU = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + electricity_cost_DE_LU);
%FI = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + electricity_cost_FI);
%SE = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + electricity_cost_SE);
%NO = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + electricity_cost_NO);
DE_LU_H = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_DE + electricity_cost_DE_LU_H);
FI_H = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_FI + electricity_cost_FI_H);
SE_H = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_SE + electricity_cost_SE_H);
NO_H = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_NO + electricity_cost_NO_H);
DE_LU_L = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_DE + electricity_cost_DE_LU_L);
FI_L = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_FI + electricity_cost_FI_L);
SE_L = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_SE + electricity_cost_SE_L);
NO_L = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_NO + electricity_cost_NO_L);
%DE_LU_P = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + electricity_cost_DE_LU_P);
%FI_P = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + electricity_cost_FI_P);
%SE_P = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + electricity_cost_SE_P);
%NO_P = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + electricity_cost_NO_P);
%avg_DE_LU = mean(DE_LU);
%avg_FI = mean(FI);
%avg_SE = mean(SE);
%avg_NO = mean(NO);
avg_DE_LU_H = mean(DE_LU_H);
avg_FI_H = mean(FI_H);
avg_SE_H = mean(SE_H);
avg_NO_H = mean(NO_H);
avg_DE_LU_L = mean(DE_LU_L);
avg_FI_L = mean(FI_L);
avg_SE_L = mean(SE_L);
avg_NO_L = mean(NO_L);
%avg_DE_LU_P = mean(DE_LU_P);
%avg_FI_P = mean(FI_P);
%avg_SE_P = mean(SE_P);
%avg_NO_P = mean(NO_P);
% Plotting LCOH for each country based on FLH
figure;
hold on;
%plot(FLH, avg_DE_LU, 'DisplayName', 'DE', 'LineWidth', 2, 'Color', [0.2,0.5,0.2]);
%plot(FLH, avg_FI, 'DisplayName', 'FI', 'LineWidth', 2, 'Color', 'm');
%plot(FLH, avg_SE, 'DisplayName', 'SE', 'LineWidth', 2, 'Color', 'b');
%plot(FLH, avg_NO, 'DisplayName', 'NO', 'LineWidth', 2, 'Color', 'r');
plot(FLH_DE, avg_DE_LU_H, '-', 'DisplayName', 'DE', 'LineWidth', 2, 'Color', [0.2,0.5,0.2]);
plot(FLH_FI, avg_FI_H, '-', 'DisplayName', 'FI', 'LineWidth', 2, 'Color', 'm');
plot(FLH_SE, avg_SE_H, '-', 'DisplayName', 'SE', 'LineWidth', 2, 'Color', 'b');
plot(FLH_NO, avg_NO_H, '-', 'DisplayName', 'NO', 'LineWidth', 2, 'Color', 'r');
plot(FLH_DE, avg_DE_LU_L, '--', 'DisplayName', 'DE', 'LineWidth', 2, 'Color', [0.2,0.5,0.2]);
plot(FLH_FI, avg_FI_L, '--', 'DisplayName', 'FI', 'LineWidth', 2, 'Color', 'm');
plot(FLH_SE, avg_SE_L, '--', 'DisplayName', 'SE', 'LineWidth', 2, 'Color', 'b');
plot(FLH_NO, avg_NO_L, '--', 'DisplayName', 'NO', 'LineWidth', 2, 'Color', 'r');
%plot(FLH, avg_DE_LU_P, ':', 'DisplayName', 'DE_P', 'LineWidth', 2, 'Color', [0.2,0.5,0.2]);
%plot(FLH, avg_FI_P, ':', 'DisplayName', 'FI_P', 'LineWidth', 2, 'Color', 'm');
%plot(FLH, avg_SE_P, ':', 'DisplayName', 'SE_P', 'LineWidth', 2, 'Color', 'b');
%plot(FLH, avg_NO_P, ':', 'DisplayName', 'NO_P', 'LineWidth', 2, 'Color', 'r');
hold off;
% Set axis labels and legend
xlabel('FLH');
ylabel('LCOH');
title('LCOH of Each Country');
legend('Location', 'Best');
grid on;
ylim([0 12]);
More Answers (1)
Maybe fill along with boundary.
num_simulations = 10000;
%Common parameters
Discount_Rate_min = 0.06; % assume 6-8%
Discount_Rate_max = 0.08;
Discount_Rate_values = unifrnd(Discount_Rate_min, Discount_Rate_max, [num_simulations, 1]);
Lifetime = 20; % years
Electricity_Cost_mean = 0.255; %EUR/kWh
Electricity_Cost_std = 0.04;
Electricity_Cost_values = normrnd(Electricity_Cost_mean, Electricity_Cost_std, [num_simulations,1]);
Electricity_Cost_values(Electricity_Cost_values < 0.02) = 0.02;
Electricity_Cost_values(Electricity_Cost_values > 0.9) = 0.9;
FLH_min = 1000;
FLH_max = 8760;
FLH = unifrnd(FLH_min, FLH_max, [num_simulations, 1]);
LHV = 33.33; %kWh/kgH2
%SOEC parameters
CAPEX_System_SOEC_mean = 4200; %$/kW
CAPEX_System_SOEC_std = 500;
CAPEX_System_SOEC_values = normrnd(CAPEX_System_SOEC_mean, CAPEX_System_SOEC_std, [num_simulations,1]);
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values < 2800) = 2800;
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values > 5600) = 5600;
CAPEX_Stack_SOEC_values = 0.5*CAPEX_System_SOEC_values; % 50% of CAPEX system
CAPEX_SOEC_values = (CAPEX_System_SOEC_values + CAPEX_Stack_SOEC_values);
OPEX_SOEC_values = 3; % 3% of CAPEX/a
System_Efficiency_SOEC_mean = 0.775;
System_Efficiency_SOEC_std = 0.05;
System_Efficiency_SOEC_values = normrnd(System_Efficiency_SOEC_mean, System_Efficiency_SOEC_std, [num_simulations,1]);
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values < 0.74) = 0.74;
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values > 0.81) = 0.81;
%PEM parameters
CAPEX_System_PEM_mean = 1450; %$/kW
CAPEX_System_PEM_std = 50;
CAPEX_System_PEM_values = normrnd(CAPEX_System_PEM_mean, CAPEX_System_PEM_std, [num_simulations,1]);
CAPEX_System_PEM_values(CAPEX_System_PEM_values < 1100) = 1100;
CAPEX_System_PEM_values(CAPEX_System_PEM_values > 1800) = 1800;
CAPEX_Stack_PEM_values = 0.35*CAPEX_System_PEM_values; % 35% of CAPEX system
CAPEX_PEM_values = (CAPEX_System_PEM_values + CAPEX_Stack_PEM_values);
OPEX_PEM_values = 3;
System_Efficiency_PEM_mean = 0.58;
System_Efficiency_PEM_std = 0.01;
System_Efficiency_PEM_values = normrnd(System_Efficiency_PEM_mean, System_Efficiency_PEM_std, [num_simulations,1]);
System_Efficiency_PEM_values(System_Efficiency_PEM_values < 0.56) = 0.56;
System_Efficiency_PEM_values(System_Efficiency_PEM_values > 0.6) = 0.6;
%AEC parameters
CAPEX_System_AEC_mean = 950; % $/kW
CAPEX_System_AEC_std = 50;
CAPEX_System_AEC_values = normrnd(CAPEX_System_AEC_mean, CAPEX_System_AEC_std, [num_simulations,1]);
CAPEX_System_AEC_values(CAPEX_System_AEC_values < 500) = 500;
CAPEX_System_AEC_values(CAPEX_System_AEC_values > 1400) = 1400;
CAPEX_Stack_AEC_values = 0.35*CAPEX_System_AEC_values; % 35% of CAPEX system
CAPEX_AEC_values = (CAPEX_System_AEC_values + CAPEX_Stack_AEC_values);
OPEX_AEC_values = 3;
System_Efficiency_AEC_mean = 0.665;
System_Efficiency_AEC_std = 0.05;
System_Efficiency_AEC_values = normrnd(System_Efficiency_AEC_mean, System_Efficiency_AEC_std, [num_simulations,1]);
System_Efficiency_AEC_values(System_Efficiency_AEC_values < 0.63) = 0.63;
System_Efficiency_AEC_values(System_Efficiency_AEC_values > 0.7) = 0.7;
% Calculate SOEC LCOH values
term1_S = LHV ./ (System_Efficiency_SOEC_values);
term2_S = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_S = (OPEX_SOEC_values / 100);
term4_S = CAPEX_SOEC_values ./ FLH;
LCOH_SOEC = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + Electricity_Cost_values);
% Calculate PEM LCOH values
term1_P = LHV ./ (System_Efficiency_PEM_values);
term2_P = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_P = (OPEX_PEM_values / 100);
term4_P = CAPEX_PEM_values ./ FLH;
LCOH_PEM = term1_P .* ((term2_P ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_P) .* term4_P + Electricity_Cost_values);
% Calculate AEC LCOH values
term1_A = LHV ./ (System_Efficiency_AEC_values);
term2_A = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_A = (OPEX_AEC_values / 100);
term4_A = CAPEX_AEC_values ./ FLH;
LCOH_AEC = term1_A .* ((term2_A ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_A) .* term4_A + Electricity_Cost_values);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
% Plot the graph
figure;
hold on
idx = boundary(FLH_filtered, LCOH_SOEC_filtered);
fill(FLH_filtered(idx),LCOH_SOEC_filtered(idx),'g','EdgeColor','none','DisplayName','SOEC','FaceAlpha',0.75)
idx = boundary(FLH_filtered, LCOH_PEM_filtered);
fill(FLH_filtered(idx),LCOH_PEM_filtered(idx),'b','EdgeColor','none','DisplayName','PEM','FaceAlpha',0.75)
idx = boundary(FLH_filtered, LCOH_AEC_filtered);
fill(FLH_filtered(idx),LCOH_AEC_filtered(idx),'y','EdgeColor','none','DisplayName','AEC','FaceAlpha',0.75)
xlabel('FLH');
ylabel('LCOH');
title('3-Sigma Approach Graph');
legend('Location', 'best');
grid on;
hold off;
Categories
Find more on Graphics Object Properties 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!
