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.
1 件のコメント
Steven Lord
2024 年 1 月 5 日
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.
採用された回答
I believe that this is reasonably close —
num_simulations = 10000;
%Common parameters
Discount_Rate_min = 0.06; % assume 6-8%
Discount_Rate_max = 0.08;
Discount_Rate_values = unifrnd(Discount_Rate_min, Discount_Rate_max, [num_simulations, 1]);
Lifetime = 20; % years
Electricity_Cost_mean = 0.255; %EUR/kWh
Electricity_Cost_std = 0.04;
Electricity_Cost_values = normrnd(Electricity_Cost_mean, Electricity_Cost_std, [num_simulations,1]);
Electricity_Cost_values(Electricity_Cost_values < 0.02) = 0.02;
Electricity_Cost_values(Electricity_Cost_values > 0.9) = 0.9;
FLH_min = 1000;
FLH_max = 8760;
FLH = unifrnd(FLH_min, FLH_max, [num_simulations, 1]);
LHV = 33.33; %kWh/kgH2
%SOEC parameters
CAPEX_System_SOEC_mean = 4200; %$/kW
CAPEX_System_SOEC_std = 500;
CAPEX_System_SOEC_values = normrnd(CAPEX_System_SOEC_mean, CAPEX_System_SOEC_std, [num_simulations,1]);
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values < 2800) = 2800;
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values > 5600) = 5600;
CAPEX_Stack_SOEC_values = 0.5*CAPEX_System_SOEC_values; % 50% of CAPEX system
CAPEX_SOEC_values = (CAPEX_System_SOEC_values + CAPEX_Stack_SOEC_values);
OPEX_SOEC_values = 3; % 3% of CAPEX/a
System_Efficiency_SOEC_mean = 0.775;
System_Efficiency_SOEC_std = 0.05;
System_Efficiency_SOEC_values = normrnd(System_Efficiency_SOEC_mean, System_Efficiency_SOEC_std, [num_simulations,1]);
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values < 0.74) = 0.74;
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values > 0.81) = 0.81;
%PEM parameters
CAPEX_System_PEM_mean = 1450; %$/kW
CAPEX_System_PEM_std = 50;
CAPEX_System_PEM_values = normrnd(CAPEX_System_PEM_mean, CAPEX_System_PEM_std, [num_simulations,1]);
CAPEX_System_PEM_values(CAPEX_System_PEM_values < 1100) = 1100;
CAPEX_System_PEM_values(CAPEX_System_PEM_values > 1800) = 1800;
CAPEX_Stack_PEM_values = 0.35*CAPEX_System_PEM_values; % 35% of CAPEX system
CAPEX_PEM_values = (CAPEX_System_PEM_values + CAPEX_Stack_PEM_values);
OPEX_PEM_values = 3;
System_Efficiency_PEM_mean = 0.58;
System_Efficiency_PEM_std = 0.01;
System_Efficiency_PEM_values = normrnd(System_Efficiency_PEM_mean, System_Efficiency_PEM_std, [num_simulations,1]);
System_Efficiency_PEM_values(System_Efficiency_PEM_values < 0.56) = 0.56;
System_Efficiency_PEM_values(System_Efficiency_PEM_values > 0.6) = 0.6;
%AEC parameters
CAPEX_System_AEC_mean = 950; % $/kW
CAPEX_System_AEC_std = 50;
CAPEX_System_AEC_values = normrnd(CAPEX_System_AEC_mean, CAPEX_System_AEC_std, [num_simulations,1]);
CAPEX_System_AEC_values(CAPEX_System_AEC_values < 500) = 500;
CAPEX_System_AEC_values(CAPEX_System_AEC_values > 1400) = 1400;
CAPEX_Stack_AEC_values = 0.35*CAPEX_System_AEC_values; % 35% of CAPEX system
CAPEX_AEC_values = (CAPEX_System_AEC_values + CAPEX_Stack_AEC_values);
OPEX_AEC_values = 3;
System_Efficiency_AEC_mean = 0.665;
System_Efficiency_AEC_std = 0.05;
System_Efficiency_AEC_values = normrnd(System_Efficiency_AEC_mean, System_Efficiency_AEC_std, [num_simulations,1]);
System_Efficiency_AEC_values(System_Efficiency_AEC_values < 0.63) = 0.63;
System_Efficiency_AEC_values(System_Efficiency_AEC_values > 0.7) = 0.7;
% Calculate SOEC LCOH values
term1_S = LHV ./ (System_Efficiency_SOEC_values);
term2_S = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_S = (OPEX_SOEC_values / 100);
term4_S = CAPEX_SOEC_values ./ FLH;
LCOH_SOEC = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + Electricity_Cost_values);
% Calculate PEM LCOH values
term1_P = LHV ./ (System_Efficiency_PEM_values);
term2_P = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_P = (OPEX_PEM_values / 100);
term4_P = CAPEX_PEM_values ./ FLH;
LCOH_PEM = term1_P .* ((term2_P ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_P) .* term4_P + Electricity_Cost_values);
% Calculate AEC LCOH values
term1_A = LHV ./ (System_Efficiency_AEC_values);
term2_A = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_A = (OPEX_AEC_values / 100);
term4_A = CAPEX_AEC_values ./ FLH;
LCOH_AEC = term1_A .* ((term2_A ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_A) .* term4_A + Electricity_Cost_values);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
[FLH_filtereds,sidx] = sort(FLH_filtered);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
LCOH_SOEC_filteredx = movmax(LCOH_SOEC_filtered(sidx), 100);
LCOH_SOEC_filteredn = movmin(LCOH_SOEC_filtered(sidx), 100);
Px1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredx, 3);
Pn1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredn, 3);
Fx1 = polyval(Px1, FLH_filtereds);
Fn1 = polyval(Pn1, FLH_filtereds);
LCOH_PEM_filteredx = movmax(LCOH_PEM_filtered(sidx), 100);
LCOH_PEM_filteredn = movmin(LCOH_PEM_filtered(sidx), 100);
Px2 = polyfit(FLH_filtereds, LCOH_PEM_filteredx, 3);
Pn2 = polyfit(FLH_filtereds, LCOH_PEM_filteredn, 3);
Fx2 = polyval(Px2, FLH_filtereds);
Fn2 = polyval(Pn2, FLH_filtereds);
LCOH_AEC_filteredx = movmax(LCOH_AEC_filtered(sidx), 100);
LCOH_AEC_filteredn = movmin(LCOH_AEC_filtered(sidx), 100);
Px3 = polyfit(FLH_filtereds, LCOH_AEC_filteredx, 3);
Pn3 = polyfit(FLH_filtereds, LCOH_AEC_filteredn, 3);
Fx3 = polyval(Px3, FLH_filtereds);
Fn3 = polyval(Pn3, FLH_filtereds);
% Plot the graph
figure;
hs{1} = scatter(FLH_filtered, LCOH_SOEC_filtered, 'o', 'DisplayName', 'SOEC');
hold on;
hs{2} = scatter(FLH_filtered, LCOH_PEM_filtered, 'x', 'DisplayName', 'PEM');
hs{3} = scatter(FLH_filtered, LCOH_AEC_filtered, '+', 'DisplayName', 'AEC');
hold off;
xlabel('FLH');
ylabel('LCOH');
title('3-Sigma Approach Graph');
legend([hs{:}], 'Location', 'best');
grid on;
figure;
hp{1} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx1; flip(Fn1)], 'r', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'SOEC');
hold on;
hp{2} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx2; flip(Fn2)], 'g', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'PEM');
hp{3} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx3; flip(Fn3)], 'b', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'AEC');
hold off;
xlabel('FLH');
ylabel('LCOH');
title('3-Sigma Approach Graph');
legend([hp{:}], 'Location', 'best');
grid on;
I am not certain what result you want. This finds a moving maximum and minimum of every group of data (after sorting the independent variable in order to sort the others), and then does a simple 3-degree polynomial regression to provide a smooth upper and lower edge to the patch plots, and then draws the patch plots.
.
7 件のコメント
Minhee
2024 年 1 月 5 日
Thank you so much!
But can I add mean value lines on that graph like below graph?
As always, my pleasure!
Yes.
I added one movmean call for each data set, and then added corresponding polyfit and polyval calls to perform the identical polynomial fits as with the others. I noted the movmean calls as ‘% ADDED’ and the polyfit and polyval results are ‘Pm#’ and ‘Fm#’ respectively, where ‘#’ denotes the numbers of the variables (as I have named them) that correspond with the others. I also added them to the legend call. If you do not want them in the legend, change it to:
legend([hp{1:3}], 'Location', 'best');
and the legend entries for the mean value line plots will not appear. (I tested that to be certain.)
Try this —
num_simulations = 10000;
%Common parameters
Discount_Rate_min = 0.06; % assume 6-8%
Discount_Rate_max = 0.08;
Discount_Rate_values = unifrnd(Discount_Rate_min, Discount_Rate_max, [num_simulations, 1]);
Lifetime = 20; % years
Electricity_Cost_mean = 0.255; %EUR/kWh
Electricity_Cost_std = 0.04;
Electricity_Cost_values = normrnd(Electricity_Cost_mean, Electricity_Cost_std, [num_simulations,1]);
Electricity_Cost_values(Electricity_Cost_values < 0.02) = 0.02;
Electricity_Cost_values(Electricity_Cost_values > 0.9) = 0.9;
FLH_min = 1000;
FLH_max = 8760;
FLH = unifrnd(FLH_min, FLH_max, [num_simulations, 1]);
LHV = 33.33; %kWh/kgH2
%SOEC parameters
CAPEX_System_SOEC_mean = 4200; %$/kW
CAPEX_System_SOEC_std = 500;
CAPEX_System_SOEC_values = normrnd(CAPEX_System_SOEC_mean, CAPEX_System_SOEC_std, [num_simulations,1]);
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values < 2800) = 2800;
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values > 5600) = 5600;
CAPEX_Stack_SOEC_values = 0.5*CAPEX_System_SOEC_values; % 50% of CAPEX system
CAPEX_SOEC_values = (CAPEX_System_SOEC_values + CAPEX_Stack_SOEC_values);
OPEX_SOEC_values = 3; % 3% of CAPEX/a
System_Efficiency_SOEC_mean = 0.775;
System_Efficiency_SOEC_std = 0.05;
System_Efficiency_SOEC_values = normrnd(System_Efficiency_SOEC_mean, System_Efficiency_SOEC_std, [num_simulations,1]);
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values < 0.74) = 0.74;
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values > 0.81) = 0.81;
%PEM parameters
CAPEX_System_PEM_mean = 1450; %$/kW
CAPEX_System_PEM_std = 50;
CAPEX_System_PEM_values = normrnd(CAPEX_System_PEM_mean, CAPEX_System_PEM_std, [num_simulations,1]);
CAPEX_System_PEM_values(CAPEX_System_PEM_values < 1100) = 1100;
CAPEX_System_PEM_values(CAPEX_System_PEM_values > 1800) = 1800;
CAPEX_Stack_PEM_values = 0.35*CAPEX_System_PEM_values; % 35% of CAPEX system
CAPEX_PEM_values = (CAPEX_System_PEM_values + CAPEX_Stack_PEM_values);
OPEX_PEM_values = 3;
System_Efficiency_PEM_mean = 0.58;
System_Efficiency_PEM_std = 0.01;
System_Efficiency_PEM_values = normrnd(System_Efficiency_PEM_mean, System_Efficiency_PEM_std, [num_simulations,1]);
System_Efficiency_PEM_values(System_Efficiency_PEM_values < 0.56) = 0.56;
System_Efficiency_PEM_values(System_Efficiency_PEM_values > 0.6) = 0.6;
%AEC parameters
CAPEX_System_AEC_mean = 950; % $/kW
CAPEX_System_AEC_std = 50;
CAPEX_System_AEC_values = normrnd(CAPEX_System_AEC_mean, CAPEX_System_AEC_std, [num_simulations,1]);
CAPEX_System_AEC_values(CAPEX_System_AEC_values < 500) = 500;
CAPEX_System_AEC_values(CAPEX_System_AEC_values > 1400) = 1400;
CAPEX_Stack_AEC_values = 0.35*CAPEX_System_AEC_values; % 35% of CAPEX system
CAPEX_AEC_values = (CAPEX_System_AEC_values + CAPEX_Stack_AEC_values);
OPEX_AEC_values = 3;
System_Efficiency_AEC_mean = 0.665;
System_Efficiency_AEC_std = 0.05;
System_Efficiency_AEC_values = normrnd(System_Efficiency_AEC_mean, System_Efficiency_AEC_std, [num_simulations,1]);
System_Efficiency_AEC_values(System_Efficiency_AEC_values < 0.63) = 0.63;
System_Efficiency_AEC_values(System_Efficiency_AEC_values > 0.7) = 0.7;
% Calculate SOEC LCOH values
term1_S = LHV ./ (System_Efficiency_SOEC_values);
term2_S = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_S = (OPEX_SOEC_values / 100);
term4_S = CAPEX_SOEC_values ./ FLH;
LCOH_SOEC = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + Electricity_Cost_values);
% Calculate PEM LCOH values
term1_P = LHV ./ (System_Efficiency_PEM_values);
term2_P = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_P = (OPEX_PEM_values / 100);
term4_P = CAPEX_PEM_values ./ FLH;
LCOH_PEM = term1_P .* ((term2_P ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_P) .* term4_P + Electricity_Cost_values);
% Calculate AEC LCOH values
term1_A = LHV ./ (System_Efficiency_AEC_values);
term2_A = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_A = (OPEX_AEC_values / 100);
term4_A = CAPEX_AEC_values ./ FLH;
LCOH_AEC = term1_A .* ((term2_A ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_A) .* term4_A + Electricity_Cost_values);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
[FLH_filtereds,sidx] = sort(FLH_filtered);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
LCOH_SOEC_filteredx = movmax(LCOH_SOEC_filtered(sidx), 100);
LCOH_SOEC_filteredn = movmin(LCOH_SOEC_filtered(sidx), 100);
LCOH_SOEC_filteredm = movmean(LCOH_SOEC_filtered(sidx), 100); % ADDED
Px1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredx, 3);
Pn1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredn, 3);
Pm1 = polyfit(FLH_filtereds, LCOH_SOEC_filteredm, 3);
Fx1 = polyval(Px1, FLH_filtereds);
Fn1 = polyval(Pn1, FLH_filtereds);
Fm1 = polyval(Pm1, FLH_filtereds);
LCOH_PEM_filteredx = movmax(LCOH_PEM_filtered(sidx), 100);
LCOH_PEM_filteredn = movmin(LCOH_PEM_filtered(sidx), 100);
LCOH_PEM_filteredm = movmean(LCOH_PEM_filtered(sidx), 100); % ADDED
Px2 = polyfit(FLH_filtereds, LCOH_PEM_filteredx, 3);
Pn2 = polyfit(FLH_filtereds, LCOH_PEM_filteredn, 3);
Pm2 = polyfit(FLH_filtereds, LCOH_PEM_filteredm, 3);
Fx2 = polyval(Px2, FLH_filtereds);
Fn2 = polyval(Pn2, FLH_filtereds);
Fm2 = polyval(Pm2, FLH_filtereds);
LCOH_AEC_filteredx = movmax(LCOH_AEC_filtered(sidx), 100);
LCOH_AEC_filteredn = movmin(LCOH_AEC_filtered(sidx), 100);
LCOH_AEC_filteredm = movmean(LCOH_AEC_filtered(sidx), 100); % ADDED
Px3 = polyfit(FLH_filtereds, LCOH_AEC_filteredx, 3);
Pn3 = polyfit(FLH_filtereds, LCOH_AEC_filteredn, 3);
Pm3 = polyfit(FLH_filtereds, LCOH_AEC_filteredm, 3);
Fx3 = polyval(Px3, FLH_filtereds);
Fn3 = polyval(Pn3, FLH_filtereds);
Fm3 = polyval(Pm3, FLH_filtereds);
% % % Plot the graph
% % figure;
% % hs{1} = scatter(FLH_filtered, LCOH_SOEC_filtered, 'o', 'DisplayName', 'SOEC');
% % hold on;
% % hs{2} = scatter(FLH_filtered, LCOH_PEM_filtered, 'x', 'DisplayName', 'PEM');
% % hs{3} = scatter(FLH_filtered, LCOH_AEC_filtered, '+', 'DisplayName', 'AEC');
% % hold off;
% % xlabel('FLH');
% % ylabel('LCOH');
% % title('3-Sigma Approach Graph');
% % legend([hs{:}], 'Location', 'best');
% % grid on;
figure;
hp{1} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx1; flip(Fn1)], 'r', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'SOEC');
hold on;
hp{4} = plot(FLH_filtereds, Fm1, '-r', 'LineWidth',2, 'DisplayName','\mu(SOEC)');
hp{2} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx2; flip(Fn2)], 'g', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'PEM');
hp{5} = plot(FLH_filtereds, Fm2, '-g', 'LineWidth',2, 'DisplayName','\mu(PEM)');
hp{3} = patch([FLH_filtereds; flip(FLH_filtereds)], [Fx3; flip(Fn3)], 'b', 'EdgeColor','none', 'FAceAlpha',0.5, 'DisplayName', 'AEC');
hp{6} = plot(FLH_filtereds, Fm3, '-b', 'LineWidth',2, 'DisplayName','\mu(AEC)');
hold off;
xlabel('FLH');
ylabel('LCOH');
title('3-Sigma Approach Graph');
legend([hp{:}], 'Location', 'best');
grid on;
.
Minhee
2024 年 1 月 9 日
Sure! No apologies necessary.
This runs (in less than 1 second on my offline computer), however I cannot get it to run here. I get a red pop-up that says ‘Something went wrong ...’
num_simulations = 10000;
%Common parameters
Discount_Rate_min = 0.06; % assume 6-8%
Discount_Rate_max = 0.08;
Discount_Rate_values = unifrnd(Discount_Rate_min, Discount_Rate_max, [num_simulations, 1]);
Lifetime = 20; % years
Electricity_Cost_values = 0.38;
FLH_min = 1000;
FLH_max = 8760;
FLH = unifrnd(FLH_min, FLH_max, [num_simulations, 1]);
LHV = 33.33; %kWh/kgH2
%SOEC parameters
CAPEX_System_SOEC_mean = 750; %$/kW
CAPEX_System_SOEC_std = 50;
CAPEX_System_SOEC_values = normrnd(CAPEX_System_SOEC_mean, CAPEX_System_SOEC_std, [num_simulations,1]);
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values < 500) = 500;
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values > 1000) = 1000;
CAPEX_Stack_SOEC_values = 0.5*CAPEX_System_SOEC_values; % 50% of CAPEX system
CAPEX_SOEC_values = (CAPEX_System_SOEC_values + CAPEX_Stack_SOEC_values);
OPEX_SOEC_values = 3; % 3% of CAPEX/a
System_Efficiency_SOEC_mean = 0.835;
System_Efficiency_SOEC_std = 0.01;
System_Efficiency_SOEC_values = normrnd(System_Efficiency_SOEC_mean, System_Efficiency_SOEC_std, [num_simulations,1]);
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values < 0.77) = 0.77;
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values > 0.9) = 0.9;
%PEM parameters
CAPEX_System_PEM_mean = 550; %$/kW
CAPEX_System_PEM_std = 50;
CAPEX_System_PEM_values = normrnd(CAPEX_System_PEM_mean, CAPEX_System_PEM_std, [num_simulations,1]);
CAPEX_System_PEM_values(CAPEX_System_PEM_values < 200) = 200;
CAPEX_System_PEM_values(CAPEX_System_PEM_values > 900) = 900;
CAPEX_Stack_PEM_values = 0.35*CAPEX_System_PEM_values; % 35% of CAPEX system
CAPEX_PEM_values = (CAPEX_System_PEM_values + CAPEX_Stack_PEM_values);
OPEX_PEM_values = 3;
System_Efficiency_PEM_mean = 0.705;
System_Efficiency_PEM_std = 0.01;
System_Efficiency_PEM_values = normrnd(System_Efficiency_PEM_mean, System_Efficiency_PEM_std, [num_simulations,1]);
System_Efficiency_PEM_values(System_Efficiency_PEM_values < 0.67) = 0.67;
System_Efficiency_PEM_values(System_Efficiency_PEM_values > 0.74) = 0.74;
%AEC parameters
CAPEX_System_AEC_mean = 450; % $/kW
CAPEX_System_AEC_std = 50;
CAPEX_System_AEC_values = normrnd(CAPEX_System_AEC_mean, CAPEX_System_AEC_std, [num_simulations,1]);
CAPEX_System_AEC_values(CAPEX_System_AEC_values < 200) = 200;
CAPEX_System_AEC_values(CAPEX_System_AEC_values > 700) = 700;
CAPEX_Stack_AEC_values = 0.35*CAPEX_System_AEC_values; % 35% of CAPEX system
CAPEX_AEC_values = (CAPEX_System_AEC_values + CAPEX_Stack_AEC_values);
OPEX_AEC_values = 3;
System_Efficiency_AEC_mean = 0.75;
System_Efficiency_AEC_std = 0.01;
System_Efficiency_AEC_values = normrnd(System_Efficiency_AEC_mean, System_Efficiency_AEC_std, [num_simulations,1]);
System_Efficiency_AEC_values(System_Efficiency_AEC_values < 0.7) = 0.7;
System_Efficiency_AEC_values(System_Efficiency_AEC_values > 0.8) = 0.8;
% Calculate SOEC LCOH values
term1_S = LHV ./ (System_Efficiency_SOEC_values);
term2_S = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_S = (OPEX_SOEC_values / 100);
term4_S = CAPEX_SOEC_values ./ FLH;
LCOH_SOEC = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + Electricity_Cost_values);
% Calculate PEM LCOH values
term1_P = LHV ./ (System_Efficiency_PEM_values);
term2_P = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_P = (OPEX_PEM_values / 100);
term4_P = CAPEX_PEM_values ./ FLH;
LCOH_PEM = term1_P .* ((term2_P ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_P) .* term4_P + Electricity_Cost_values);
% Calculate AEC LCOH values
term1_A = LHV ./ (System_Efficiency_AEC_values);
term2_A = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_A = (OPEX_AEC_values / 100);
term4_A = CAPEX_AEC_values ./ FLH;
LCOH_AEC = term1_A .* ((term2_A ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_A) .* term4_A + Electricity_Cost_values);
% % Plotting LCOH vs FLH for SOEC, PEM, and AEC in one graph
% figure;
% % Plotting for SOEC
% plot(FLH, LCOH_SOEC, 'o', 'DisplayName', 'SOEC');
% hold on;
% % Plotting for PEM
% plot(FLH, LCOH_PEM, 'o', 'DisplayName', 'PEM');
% % Plotting for AEC
% plot(FLH, LCOH_AEC, 'o', 'DisplayName', 'AEC');
% title('LCOH vs FLH for Different Technologies');
% xlabel('FLH (hours/year)');
% ylabel('LCOH ($/kgH2)');
% legend('Location', 'Best');
% grid on;
% hold off;
[LCOH_SOECmax,LCOH_SOECmin,LCOH_SOECmean,FLHs] = scatterfilter(FLH,LCOH_SOEC);
[LCOH_PEMmax,LCOH_PEMmin,LCOH_PEMmean,FLHs] = scatterfilter(FLH,LCOH_PEM);
[LCOH_AECmax,LCOH_AECmin,LCOH_AECmean,FLHs] = scatterfilter(FLH,LCOH_AEC);
toc
% % % Plot the graph
% % figure;
% % hs{1} = scatter(FLH_filtered, LCOH_SOEC_filtered, 'o', 'DisplayName', 'SOEC');
% % hold on;
% % hs{2} = scatter(FLH_filtered, LCOH_PEM_filtered, 'x', 'DisplayName', 'PEM');
% % hs{3} = scatter(FLH_filtered, LCOH_AEC_filtered, '+', 'DisplayName', 'AEC');
% % hold off;
% % xlabel('FLH');
% % ylabel('LCOH');
% % title('3-Sigma Approach Graph');
% % legend([hs{:}], 'Location', 'best');
% % grid on;
figure;
hp{1} = patch([FLHs; flip(FLHs)], [LCOH_SOECmin; flip(LCOH_SOECmax)], 'r', 'EdgeColor','none', 'FaceAlpha',0.5, 'DisplayName', 'SOEC');
hold on;
hp{4} = plot(FLHs, LCOH_SOECmean, '-r', 'LineWidth',2, 'DisplayName','\mu(SOEC)');
hp{2} = patch([FLHs; flip(FLHs)], [LCOH_PEMmax; flip(LCOH_PEMmin)], 'g', 'EdgeColor','none', 'FaceAlpha',0.5, 'DisplayName', 'PEM');
hp{5} = plot(FLHs, LCOH_PEMmean, '-g', 'LineWidth',2, 'DisplayName','\mu(PEM)');
hp{3} = patch([FLHs; flip(FLHs)], [LCOH_AECmax; flip(LCOH_AECmin)], 'b', 'EdgeColor','none', 'FaceAlpha',0.5, 'DisplayName', 'AEC');
hp{6} = plot(FLHs, LCOH_AECmean, '-b', 'LineWidth',2, 'DisplayName','\mu(AEC)');
hold off;
xlabel('FLH');
ylabel('LCOH');
title('3-Sigma Approach Graph');
legend([hp{:}], 'Location', 'best');
grid on;
toc
function [Fmax,Fmin,Fmean,xs] = scatterfilter(x,y)
[xs,sidx] = sort(x,'ascend');
ymax = movmax(y(sidx), 100);
ymin = movmin(y(sidx), 100);
ymean = movmean(y(sidx), 100); % ADDED
Px1 = polyfit(xs, ymax, 3);
Pn1 = polyfit(xs, ymin, 3);
Pm1 = polyfit(xs, ymean, 3);
Fmax = polyval(Px1, xs);
Fmin = polyval(Pn1, xs);
Fmean = polyval(Pm1, xs);
end
It produces this result —
EDIT — (9 Jan 2024 at 16:40)
Returned to check it again and without any changes to the code, it now runs. I have no idea what the problem was previously.
.
Minhee
2024 年 1 月 23 日
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
2024 年 2 月 6 日
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]);
その他の回答 (1 件)
Maybe fill along with boundary.
num_simulations = 10000;
%Common parameters
Discount_Rate_min = 0.06; % assume 6-8%
Discount_Rate_max = 0.08;
Discount_Rate_values = unifrnd(Discount_Rate_min, Discount_Rate_max, [num_simulations, 1]);
Lifetime = 20; % years
Electricity_Cost_mean = 0.255; %EUR/kWh
Electricity_Cost_std = 0.04;
Electricity_Cost_values = normrnd(Electricity_Cost_mean, Electricity_Cost_std, [num_simulations,1]);
Electricity_Cost_values(Electricity_Cost_values < 0.02) = 0.02;
Electricity_Cost_values(Electricity_Cost_values > 0.9) = 0.9;
FLH_min = 1000;
FLH_max = 8760;
FLH = unifrnd(FLH_min, FLH_max, [num_simulations, 1]);
LHV = 33.33; %kWh/kgH2
%SOEC parameters
CAPEX_System_SOEC_mean = 4200; %$/kW
CAPEX_System_SOEC_std = 500;
CAPEX_System_SOEC_values = normrnd(CAPEX_System_SOEC_mean, CAPEX_System_SOEC_std, [num_simulations,1]);
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values < 2800) = 2800;
CAPEX_System_SOEC_values(CAPEX_System_SOEC_values > 5600) = 5600;
CAPEX_Stack_SOEC_values = 0.5*CAPEX_System_SOEC_values; % 50% of CAPEX system
CAPEX_SOEC_values = (CAPEX_System_SOEC_values + CAPEX_Stack_SOEC_values);
OPEX_SOEC_values = 3; % 3% of CAPEX/a
System_Efficiency_SOEC_mean = 0.775;
System_Efficiency_SOEC_std = 0.05;
System_Efficiency_SOEC_values = normrnd(System_Efficiency_SOEC_mean, System_Efficiency_SOEC_std, [num_simulations,1]);
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values < 0.74) = 0.74;
System_Efficiency_SOEC_values(System_Efficiency_SOEC_values > 0.81) = 0.81;
%PEM parameters
CAPEX_System_PEM_mean = 1450; %$/kW
CAPEX_System_PEM_std = 50;
CAPEX_System_PEM_values = normrnd(CAPEX_System_PEM_mean, CAPEX_System_PEM_std, [num_simulations,1]);
CAPEX_System_PEM_values(CAPEX_System_PEM_values < 1100) = 1100;
CAPEX_System_PEM_values(CAPEX_System_PEM_values > 1800) = 1800;
CAPEX_Stack_PEM_values = 0.35*CAPEX_System_PEM_values; % 35% of CAPEX system
CAPEX_PEM_values = (CAPEX_System_PEM_values + CAPEX_Stack_PEM_values);
OPEX_PEM_values = 3;
System_Efficiency_PEM_mean = 0.58;
System_Efficiency_PEM_std = 0.01;
System_Efficiency_PEM_values = normrnd(System_Efficiency_PEM_mean, System_Efficiency_PEM_std, [num_simulations,1]);
System_Efficiency_PEM_values(System_Efficiency_PEM_values < 0.56) = 0.56;
System_Efficiency_PEM_values(System_Efficiency_PEM_values > 0.6) = 0.6;
%AEC parameters
CAPEX_System_AEC_mean = 950; % $/kW
CAPEX_System_AEC_std = 50;
CAPEX_System_AEC_values = normrnd(CAPEX_System_AEC_mean, CAPEX_System_AEC_std, [num_simulations,1]);
CAPEX_System_AEC_values(CAPEX_System_AEC_values < 500) = 500;
CAPEX_System_AEC_values(CAPEX_System_AEC_values > 1400) = 1400;
CAPEX_Stack_AEC_values = 0.35*CAPEX_System_AEC_values; % 35% of CAPEX system
CAPEX_AEC_values = (CAPEX_System_AEC_values + CAPEX_Stack_AEC_values);
OPEX_AEC_values = 3;
System_Efficiency_AEC_mean = 0.665;
System_Efficiency_AEC_std = 0.05;
System_Efficiency_AEC_values = normrnd(System_Efficiency_AEC_mean, System_Efficiency_AEC_std, [num_simulations,1]);
System_Efficiency_AEC_values(System_Efficiency_AEC_values < 0.63) = 0.63;
System_Efficiency_AEC_values(System_Efficiency_AEC_values > 0.7) = 0.7;
% Calculate SOEC LCOH values
term1_S = LHV ./ (System_Efficiency_SOEC_values);
term2_S = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_S = (OPEX_SOEC_values / 100);
term4_S = CAPEX_SOEC_values ./ FLH;
LCOH_SOEC = term1_S .* ((term2_S ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_S) .* term4_S + Electricity_Cost_values);
% Calculate PEM LCOH values
term1_P = LHV ./ (System_Efficiency_PEM_values);
term2_P = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_P = (OPEX_PEM_values / 100);
term4_P = CAPEX_PEM_values ./ FLH;
LCOH_PEM = term1_P .* ((term2_P ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_P) .* term4_P + Electricity_Cost_values);
% Calculate AEC LCOH values
term1_A = LHV ./ (System_Efficiency_AEC_values);
term2_A = Discount_Rate_values .* (1 + Discount_Rate_values).^Lifetime;
term3_A = (OPEX_AEC_values / 100);
term4_A = CAPEX_AEC_values ./ FLH;
LCOH_AEC = term1_A .* ((term2_A ./ ((1 + Discount_Rate_values).^Lifetime - 1) + term3_A) .* term4_A + Electricity_Cost_values);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
% Calculate mean and standard deviation of Electricity_Cost_values
electricity_cost_mean = mean(Electricity_Cost_values);
electricity_cost_std = std(Electricity_Cost_values);
% Define the 3-sigma range
lower_limit = electricity_cost_mean - 3 * electricity_cost_std;
upper_limit = electricity_cost_mean + 3 * electricity_cost_std;
% Filter data based on the 3-sigma range
valid_indices = (Electricity_Cost_values >= lower_limit) & (Electricity_Cost_values <= upper_limit);
FLH_filtered = FLH(valid_indices);
LCOH_SOEC_filtered = LCOH_SOEC(valid_indices);
LCOH_PEM_filtered = LCOH_PEM(valid_indices);
LCOH_AEC_filtered = LCOH_AEC(valid_indices);
% Plot the graph
figure;
hold on
idx = boundary(FLH_filtered, LCOH_SOEC_filtered);
fill(FLH_filtered(idx),LCOH_SOEC_filtered(idx),'g','EdgeColor','none','DisplayName','SOEC','FaceAlpha',0.75)
idx = boundary(FLH_filtered, LCOH_PEM_filtered);
fill(FLH_filtered(idx),LCOH_PEM_filtered(idx),'b','EdgeColor','none','DisplayName','PEM','FaceAlpha',0.75)
idx = boundary(FLH_filtered, LCOH_AEC_filtered);
fill(FLH_filtered(idx),LCOH_AEC_filtered(idx),'y','EdgeColor','none','DisplayName','AEC','FaceAlpha',0.75)
xlabel('FLH');
ylabel('LCOH');
title('3-Sigma Approach Graph');
legend('Location', 'best');
grid on;
hold off;
カテゴリ
ヘルプ センター および File Exchange で Graphics Object Properties についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
