Add a line on the graph

How can I add the red line following the bar shape on the graph? (like below image)

The 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_values = 0.255;      %EUR/kWh

FLH = 4380;

LHV = 33.33;    %kWh/kgH2

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

% Create histograms for SOEC, AEC, and PEM LCOH values
figure;

% SOEC LCOH
subplot(3, 1, 1);
histogram(LCOH_SOEC, 'Normalization', 'probability', 'EdgeColor', 'w');
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('SOEC LCOH Distribution');
ylim([0, 0.1]);

% AEC LCOH
subplot(3, 1, 2);
histogram(LCOH_AEC, 'Normalization', 'probability', 'EdgeColor', 'w');
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('AEC LCOH Distribution');
ylim([0, 0.1]);

% PEM LCOH
subplot(3, 1, 3);
histogram(LCOH_PEM, 'Normalization', 'probability', 'EdgeColor', 'w');
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('PEM LCOH Distribution');
ylim([0, 0.1]);

% Adjust layout
sgtitle('Monte Carlo Simulation Results for SOEC, AEC, and PEM');

% Show the plot

[Dyuman Joshi]

Should be easy to do with plot.

What seems to be the problem?

What does the red line denote?

Answer 1

[Star Strider]

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_values = 0.255;      %EUR/kWh
FLH = 4380;
LHV = 33.33;    %kWh/kgH2
%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_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);
% Create histograms for SOEC, AEC, and PEM LCOH values
figure;
% SOEC LCOH
subplot(3, 1, 1);
histogram(LCOH_SOEC, 'Normalization', 'probability', 'EdgeColor', 'w');
Ax1 = gca;
drawRedLine(Ax1);
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('SOEC LCOH Distribution');
ylim([0, 0.1]);
% AEC LCOH
subplot(3, 1, 2);
histogram(LCOH_AEC, 'Normalization', 'probability', 'EdgeColor', 'w');
Ax2 = gca;
drawRedLine(Ax2);
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('AEC LCOH Distribution');
ylim([0, 0.1]);
% PEM LCOH
subplot(3, 1, 3);
histogram(LCOH_PEM, 'Normalization', 'probability', 'EdgeColor', 'w');
Ax3 = gca;
drawRedLine(Ax3);
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('PEM LCOH Distribution');
ylim([0, 0.1]);
% Adjust layout
sgtitle('Monte Carlo Simulation Results for SOEC, AEC, and PEM');

% Show the plot

% function drawRedLine(axh)
% %   Draws the red line in each histogram subplot.
%     H = axh.Children;
%     BinCtrs = H.BinEdges(1:end-1)+diff(H.BinEdges)/2;
%     B = polyfit(BinCtrs,H.Values,2) ;
%     RedLine = polyval(B, BinCtrs);
%     FudgeFactor = max(H.Values) - max(RedLine);
%     hold on
%     plot(axh,BinCtrs, RedLine  + FudgeFactor,  '-r', 'LineWidth',1)
%     hold off
% end

function drawRedLine(axh)
%   Draws the red line in each histogram subplot.
    H = axh.Children;
    BinCtrs = H.BinEdges(1:end-1)+diff(H.BinEdges)/2;
    Data = H.Values;
    u = mean(BinCtrs);
    s = std(Data);
    ndpdf = @(b,x) b(1) .* exp(-(x-b(2)).^2 ./ b(3));
    B = fminsearch(@(b)norm(Data-ndpdf(b,BinCtrs)), [max(H.Values); u; s]);
    % text(axh, B(2), max(H.Values), sprintf('$p(x) = %.3fe^{-\\frac{(%.2f-x)^2}{%.3f}}$', B), 'Horiz','center','Vert','bottom','Interpreter','latex')
    RedLine = ndpdf(B, H.BinEdges);
    hold on
    plot(axh,H.BinEdges, RedLine,  '-r', 'LineWidth',1.5)
    hold off
end

EDIT — (17 Jan 2024 at 04:57)

It occurred to me later that the red line could be a normal probability distribution, although exactly what it is has not yet been divulged. If so, my new ‘drawRedLine’ function may be correct. (The previous one was just a guess.) The distribution parameters are the ‘B’ vector, whose contents are not shown here, however would be readily available. I provided a text call that presents them (commented-out) in case that information is desired.

.

[Minhee]

Thank you for answer.

That's exactly what I want!

But how can I plot these three histogram(with red line) on one graph together to compare?

[Star Strider]

As always, my pleasure!

Thiose are difficult to put on one plot and have the plot make sense, however it can be done. Here, I coloured the distribution line plots the same as the bar colours, however I included a commented-out-option that draws all the lines as red. The ‘secret’ to making this straightforward is to change my ‘drawRedLine’ function to return its data as outputs.

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_values = 0.255;      %EUR/kWh
FLH = 4380;
LHV = 33.33;    %kWh/kgH2
%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_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);
% Create histograms for SOEC, AEC, and PEM LCOH values
figure;
% SOEC LCOH
subplot(3, 1, 1);
histogram(LCOH_SOEC, 'Normalization', 'probability', 'EdgeColor', 'w');
Ax1 = gca;
[RL{1},BE{1},BC{1},BD{1}] = drawRedLine(Ax1);
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('SOEC LCOH Distribution');
ylim([0, 0.1]);
% AEC LCOH
subplot(3, 1, 2);
histogram(LCOH_AEC, 'Normalization', 'probability', 'EdgeColor', 'w');
Ax2 = gca;
[RL{2},BE{2},BC{2},BD{2}] = drawRedLine(Ax2);
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('AEC LCOH Distribution');
ylim([0, 0.1]);
% PEM LCOH
subplot(3, 1, 3);
histogram(LCOH_PEM, 'Normalization', 'probability', 'EdgeColor', 'w');
Ax3 = gca;
[RL{3},BE{3},BC{3},BD{3}] = drawRedLine(Ax3);
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('PEM LCOH Distribution');
ylim([0, 0.1]);
% Adjust layout
sgtitle('Monte Carlo Simulation Results for SOEC, AEC, and PEM');

% Show the plot


DN = ["SOEC LCOH" "AEC LCOH" "PEM LCOH"];

figure
hold on
for k = 1:3
    hb = bar(BC{k}, BD{k}, 'DisplayName',[DN(k)+" Histogram"], 'FaceAlpha',0.5);
    % plot(BE{k}, RL{k}, 'LineWidth',3, 'DisplayName',[DN(k)+" Distribution"], 'Color','r')
    plot(BE{k}, RL{k}, 'LineWidth',3, 'DisplayName',[DN(k)+" Distribution"], 'Color',hb.FaceColor)
end
hold off
grid
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('Normal Distribution Plots')
legend('Location','northoutside')

% function drawRedLine(axh)
% %   Draws the red line in each histogram subplot.
%     H = axh.Children;
%     BinCtrs = H.BinEdges(1:end-1)+diff(H.BinEdges)/2;
%     B = polyfit(BinCtrs,H.Values,2) ;
%     RedLine = polyval(B, BinCtrs);
%     FudgeFactor = max(H.Values) - max(RedLine);
%     hold on
%     plot(axh,BinCtrs, RedLine  + FudgeFactor,  '-r', 'LineWidth',1)
%     hold off
% end

function [RL,BE,BC,BD] = drawRedLine(axh)                           % MODIFIED: Add Outputs
%   Draws the red line in each histogram subplot.
    H = axh.Children;
    BinCtrs = H.BinEdges(1:end-1)+diff(H.BinEdges)/2;
    Data = H.Values;
    u = mean(BinCtrs);
    s = std(Data);
    ndpdf = @(b,x) b(1) .* exp(-(x-b(2)).^2 ./ b(3));
    B = fminsearch(@(b)norm(Data-ndpdf(b,BinCtrs)), [max(H.Values); u; s]);
    % text(axh, B(2), max(H.Values), sprintf('$p(x) = %.3fe^{-\\frac{(%.2f-x)^2}{%.3f}}$', B), 'Horiz','center','Vert','bottom','Interpreter','latex')
    RedLine = ndpdf(B, H.BinEdges);
    hold on
    plot(axh,H.BinEdges, RedLine,  '-r', 'LineWidth',1.5)
    hold off
    RL = RedLine;
    BE = H.BinEdges;
    BC = BinCtrs;
    BD = Data;
end

.

[Minhee]

The result of the code is like the image...

[Star Strider]

We must be doing dofferent things.

When I run my previous code, I get this result —

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.255;      %EUR/kWh
FLH = 4380;
LHV = 33.33;    %kWh/kgH2
%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_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);
% Create histograms for SOEC, AEC, and PEM LCOH values
figure;
% SOEC LCOH
subplot(3, 1, 1);
histogram(LCOH_SOEC, 'Normalization', 'probability', 'EdgeColor', 'w');
Ax1 = gca;
[RL{1},BE{1},BC{1},BD{1}] = drawRedLine(Ax1);
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('SOEC LCOH Distribution');
ylim([0, 0.1]);
% AEC LCOH
subplot(3, 1, 2);
histogram(LCOH_AEC, 'Normalization', 'probability', 'EdgeColor', 'w');
Ax2 = gca;
[RL{2},BE{2},BC{2},BD{2}] = drawRedLine(Ax2);
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('AEC LCOH Distribution');
ylim([0, 0.1]);
% PEM LCOH
subplot(3, 1, 3);
histogram(LCOH_PEM, 'Normalization', 'probability', 'EdgeColor', 'w');
Ax3 = gca;
[RL{3},BE{3},BC{3},BD{3}] = drawRedLine(Ax3);
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('PEM LCOH Distribution');
ylim([0, 0.1]);
% Adjust layout
sgtitle('Monte Carlo Simulation Results for SOEC, AEC, and PEM');

% Show the plot


DN = ["SOEC LCOH" "AEC LCOH" "PEM LCOH"];

figure
hold on
for k = 1:3
    hb = bar(BC{k}, BD{k}, 'DisplayName',[DN(k)+" Histogram"], 'FaceAlpha',0.5);
    % plot(BE{k}, RL{k}, 'LineWidth',3, 'DisplayName',[DN(k)+" Distribution"], 'Color','r')
    plot(BE{k}, RL{k}, 'LineWidth',3, 'DisplayName',[DN(k)+" Distribution"], 'Color',hb.FaceColor)
end
hold off
grid
xlabel('LCOH (per kWh)');
ylabel('Probability');
title('Normal Distribution Plots')
legend('Location','northoutside')

% function drawRedLine(axh)
% %   Draws the red line in each histogram subplot.
%     H = axh.Children;
%     BinCtrs = H.BinEdges(1:end-1)+diff(H.BinEdges)/2;
%     B = polyfit(BinCtrs,H.Values,2) ;
%     RedLine = polyval(B, BinCtrs);
%     FudgeFactor = max(H.Values) - max(RedLine);
%     hold on
%     plot(axh,BinCtrs, RedLine  + FudgeFactor,  '-r', 'LineWidth',1)
%     hold off
% end

function [RL,BE,BC,BD] = drawRedLine(axh)                           % MODIFIED: Add Outputs
%   Draws the red line in each histogram subplot.
    H = axh.Children;
    BinCtrs = H.BinEdges(1:end-1)+diff(H.BinEdges)/2;
    Data = H.Values;
    u = mean(BinCtrs);
    s = std(Data);
    ndpdf = @(b,x) b(1) .* exp(-(x-b(2)).^2 ./ b(3));
    B = fminsearch(@(b)norm(Data-ndpdf(b,BinCtrs)), [max(H.Values); u; s]);
    % text(axh, B(2), max(H.Values), sprintf('$p(x) = %.3fe^{-\\frac{(%.2f-x)^2}{%.3f}}$', B), 'Horiz','center','Vert','bottom','Interpreter','latex')
    RedLine = ndpdf(B, H.BinEdges);
    hold on
    plot(axh,H.BinEdges, RedLine,  '-r', 'LineWidth',1.5)
    hold off
    RL = RedLine;
    BE = H.BinEdges;
    BC = BinCtrs;
    BD = Data;
end

.

[Minhee]

Even though I copy and paste your code, there is nothing on the graph...

[Star Strider]

I have no idea what the problem could be, since I am not certain how you are running it.

See if creating a new script file, pasting my code into it, and then running that script independently of everything else will work. (That is essentially what I am doing here.)

[Minhee]

Unable to perform assignment because brace indexing is not supported for variables of this type.

Error in Monte_Carlo (line 76)

[RL{1},BE{1},BC{1},BD{1}] = drawRedLine(Ax1);

I have this error. Do you know how to solve it?

[Star Strider]

My ‘drawRedLine’ function creates those. That syntax should work unless you changed the function somehow.

The other possibility is that you already declared those variable names as something other than cell arrays earlier in your code. In that instance, change the variable names in the ‘drawRedLines’ call, and everywhere you use those specific results, for example:

[dRL{1},dBE{1},dBC{1},dBD{1}] = drawRedLine(Ax1);

or something similar, to distinguish them from the previously assigned or declared variables.

.

[Minhee]

Thank you so much.

Finally I got the graph!!

[Star Strider]

As always, my pleasure!

Answer 2

[Hassaan]

% Assuming 'LCOH_SOEC' is the array containing your SOEC LCOH data
% and you've already created a histogram:
histogram(LCOH_SOEC, 'Normalization', 'probability', 'EdgeColor', 'w');
hold on; % Retain the current plot

% Now, let's say you want to add a vertical line at the mean value of the LCOH
mean_value = mean(LCOH_SOEC);
y_limits = ylim; % Gets the current limits of the y-axis

% Plot a vertical line at the mean value
plot([mean_value mean_value], y_limits, 'r-', 'LineWidth', 2);

hold off; % Release the plot

This script will add a vertical red line at the mean value of the SOEC LCOH data on the histogram. If you want to add other types of lines (horizontal, or at a specific x or y value), you would adjust the plot command accordingly.

-----------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

Professional Interests

• Technical Services and Consulting
• Embedded Systems | Firmware Developement | Simulations
• Electrical and Electronics Engineering

It's important to note that the advice and code are based on limited information and meant for educational purposes. Users should verify and adapt the code to their specific needs, ensuring compatibility and adherence to ethical standards.

Feel free to contact me.

[Adam Danz]

> Now, let's say you want to add a vertical line at the mean value of the LCOH

This answer is not helpful. The OP is not asking how to create vertical lines and the red lines in the OP's question are obviously not based on means.

@Muhammad Hassaan Shah please review the Generative AI guidelines for MATLAB Central

https://www.mathworks.com/matlabcentral/content/community-guidelines.html

The use of generative AI tools is allowed but you are responsible for ensuring that your content is relevant and accurate.

[Hassaan]

@Adam Danz Just for my understanding:

• Can you please elaborate how is it 'generative AI' if my understanding is incorrect for such case and the solution is not as per the requirements.
• There are a millions ways to tackle a problem and OP may or may not agree with the solution. How this falls under 'not helpful' or 'generative AI' or flagged as?
• What are the metrics to be classified as a generative content?.

Thank you.