How to change the graph type

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.

[Steven Lord]

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.

Answer 1

[Star Strider]

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);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Pn1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredn, 3);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

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);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Pn2 = polyfit(FLH_filtereds, LCOH_PEM_filteredn, 3);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

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);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Pn3 = polyfit(FLH_filtereds, LCOH_AEC_filteredn, 3);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

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.

.

[Minhee]

Thank you so much!

But can I add mean value lines on that graph like below graph?

[Star Strider]

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);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Pn1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredn, 3);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Pm1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredm, 3);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

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);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Pn2 = polyfit(FLH_filtereds, LCOH_PEM_filteredn, 3);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Pm2 = polyfit(FLH_filtereds, LCOH_PEM_filteredm, 3);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

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);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Pn3 = polyfit(FLH_filtereds, LCOH_AEC_filteredn, 3);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Pm3 = polyfit(FLH_filtereds, LCOH_AEC_filteredm, 3);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

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]

Hi,

I'm sorry but can I change the below graph like the above one?

I didn't use 3-sigma this time but I don't know how to change the code.

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;

[Star Strider]

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);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

[LCOH_PEMmax,LCOH_PEMmin,LCOH_PEMmean,FLHs] = scatterfilter(FLH,LCOH_PEM);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

[LCOH_AECmax,LCOH_AECmin,LCOH_AECmean,FLHs] = scatterfilter(FLH,LCOH_AEC);

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT.

toc

Elapsed time is 3.266319 seconds.

% % % 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

Elapsed time is 3.853918 seconds.

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]

I apply 3-sigma approach to LCOH, not electricity costs.

How can I plot the same graph with below code?

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;

[Star Strider]

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]

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]);

Answer 2

[Voss]

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;