Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
@hdesing
hdesing authored Nov 3, 2021
1 parent 81aeb1b commit a554ae6
Show file tree
Hide file tree
Showing 10 changed files with 750 additions and 0 deletions.
Binary file added 2018 - Global energy demand.xlsx
View file
Binary file not shown.
Binary file added CO2 emissions and remaining budget.xlsx
View file
Binary file not shown.
162 changes: 162 additions & 0 deletions Historic_and_projected_emission.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,162 @@
%% Historical and projected emissions

tic
clear variables

%% loading and preparing data
% time series for energy demand and CO2 emissions

time_frame = xlsread('CO2 emissions and remaining budget.xlsx','CO2 time series','A4:A172');

P_el_historic_by_source = xlsread('CO2 emissions and remaining budget.xlsx','CO2 time series','X4:AA172');
P_el_historic_by_source (isnan(P_el_historic_by_source))=0;

dot_m_CO2_historic_by_source = xlsread('CO2 emissions and remaining budget.xlsx','CO2 time series','B4:F172');
dot_m_CO2_historic_by_source (isnan(dot_m_CO2_historic_by_source))=0;

m_CO2_cumulative_by_source = cumsum(dot_m_CO2_historic_by_source,1);

IEA_energy_scenario_data = xlsread('CO2 emissions and remaining budget.xlsx','IEA','A21:I23');
IEA_emission_scenario_data = xlsread('CO2 emissions and remaining budget.xlsx','IEA','A27:D29');

IEA_range = 2018:1:2040;
for i=1:size(IEA_emission_scenario_data,2)-1
IEA_emission_scenario_fit = fit(IEA_emission_scenario_data(:,1),IEA_emission_scenario_data(:,i+1),'linearinterp');
IEA_emission_scenario(:,i) = IEA_emission_scenario_fit(IEA_range);
end

IEA_emission_scenario_cumulated = m_CO2_cumulative_by_source(end-1,1:3) + cumsum(IEA_emission_scenario,1);

IPCC_emission_scenario_data = xlsread('CO2 emissions and remaining budget.xlsx','IPCC','C7:F19')';
IPCC_emission_scenario_sigma = xlsread('CO2 emissions and remaining budget.xlsx','IPCC','K8:N13')';

IPCC_range = 2018:1:2100;
for i=1:12
IPCC_emission_scenario_fit = fit(IPCC_emission_scenario_data(:,1),IPCC_emission_scenario_data(:,i+1),'linearinterp');
IPCC_emission_scenario(:,i) = IPCC_emission_scenario_fit(IPCC_range);
end
for i=1:6
IPCC_emission_scenario_FF(:,i) = IPCC_emission_scenario(:,(i*2-1));
IPCC_emission_scenario_total(:,i) = IPCC_emission_scenario(:,(i*2));
end

for i=1:6
IPCC_sigma_fit = fit(IPCC_emission_scenario_data(:,1),IPCC_emission_scenario_sigma(:,i),'linearinterp');
IPCC_sigma(:,i) = IPCC_sigma_fit(IPCC_range);
end

IPCC_emission_scenario_cumulated_total = zeros(size(IPCC_range,2),6);
IPCC_emission_scenario_cumulated_FF = zeros(size(IPCC_range,2),6);
for i=1:6
IPCC_emission_scenario_cumulated_FF(:,i) = sum(m_CO2_cumulative_by_source(end-1,1:4),2) + cumsum(IPCC_emission_scenario(:,(i*2-1)),1);
IPCC_emission_scenario_cumulated_total(:,i) = sum(m_CO2_cumulative_by_source(end-1,:),2) + cumsum(IPCC_emission_scenario(:,(i*2)),1);
end


beta_disp = 0.1;
n_alpha_disp = 1;

%% Uncertainty of historic emissions

n_runs = 10^5;

mu_ff = sum(dot_m_CO2_historic_by_source(:,1:4),2); % fossil and industry emissions
sigma_ff = 0.05 .* mu_ff; % sigma is proportional to mu (Friedlingstein et al. 2019)

mu_l = dot_m_CO2_historic_by_source(:,5); % land use emissions
sigma_l = 2.56; % Gt/a (Friedlingstein et al. 2019)

dot_m_rnd = zeros(size(mu_ff,1),2,n_runs);
for i=1:size(mu_ff,1)
dot_m_rnd(i,1,:) = normrnd(mu_ff(i,1),sigma_ff(i,1),n_runs,1);
dot_m_rnd(i,2,:) = normrnd(mu_l(i,1),sigma_l,n_runs,1);
end

m_cum_rnd = cumsum(dot_m_rnd,1);


CI = 0.99; % confidence intervall (can be chosen freely)
m_cum(:,2) = quantile(sum(m_cum_rnd,2),0.5,3);
m_cum(:,1) = quantile(sum(m_cum_rnd,2),(1-CI)/2,3);
m_cum(:,3) = quantile(sum(m_cum_rnd,2),1-(1-CI)/2,3);


mu_IPCC = IPCC_emission_scenario(:,[2 4 6 8 10 12]); % total emissions
sigma_IPCC = IPCC_sigma;

dot_m_rnd_IPCC = zeros(size(IPCC_range,2),6,n_runs);
for i=1:size(IPCC_range,2)
for j=1:6
dot_m_rnd_IPCC(i,j,:) = normrnd(mu_IPCC(i,j),sigma_IPCC(i,j),n_runs,1);
end
end

m_cum_rnd_IPCC = sum(m_cum_rnd(end-1,:,:),2) + cumsum(dot_m_rnd_IPCC,1);

m_cum_IPCC(:,:,2) = quantile(m_cum_rnd_IPCC,0.5,3);
m_cum_IPCC(:,:,1) = quantile(m_cum_rnd_IPCC,(1-CI)/2,3);
m_cum_IPCC(:,:,3) = quantile(m_cum_rnd_IPCC,1-(1-CI)/2,3);

%% Calculating negative emissions in IPCC pathways

CDR_IPCC = IPCC_emission_scenario_total - IPCC_emission_scenario_FF;
for i=1:size(IPCC_range,2)
for j=1:6
if CDR_IPCC(i,j) > 0
CDR_IPCC(i,j) = 0;
end
end
end
CDR_cum = cumsum(CDR_IPCC,1);

Net_negative = max(IPCC_emission_scenario_cumulated_total,[],1) - IPCC_emission_scenario_cumulated_total(end,:);

%% Calculation of target violation

% uncertainty distribution of 1.5°C and 2°C targets

m_CO2_2017 = 2313;
mu_1_5 = 2.774; %data from regression (see "Transition model")
sigma_1_5 = 0.3967;

mu_2 = 3.196;
sigma_2 = 0.3267;

m_target_1_5 = lognrnd(log(10^(mu_1_5)),log(10^(sigma_1_5)),n_runs,1) + m_CO2_2017;
m_target_2 = lognrnd(log(10^(mu_2)),log(10^(sigma_2)),n_runs,1) + m_CO2_2017;


[IPCC_max, i_t_max] = max(IPCC_emission_scenario_cumulated_total,[],1);
t_max = IPCC_range(1,i_t_max);

%% probability of violation
n_scenarios = 6;

for i=1:n_scenarios
n_2100=0;
m_2100=0;
n_max=0;
m_max=0;
for j=1:n_runs
if m_cum_rnd_IPCC(end,i,j)>= m_target_1_5(j,1);
n_2100=n_2100+1;
end
if m_cum_rnd_IPCC(end,i,j)>= m_target_2(j,1);
m_2100=m_2100+1;
end
if m_cum_rnd_IPCC(i_t_max(1,i),i,j)>= m_target_1_5(j,1);
n_max=n_max+1;
end
if m_cum_rnd_IPCC(i_t_max(1,i),i,j)>= m_target_2(j,1);
m_max=m_max+1;
end
end
P_v_1_5_2100(i) = n_2100/n_runs;
P_v_2_2100(i) = m_2100/n_runs;
P_v_1_5_max(i) = n_max/n_runs;
P_v_2_max(i) = m_max/n_runs;
end


%%
toc
163 changes: 163 additions & 0 deletions PV_growth_fitting_to_transition_model.m
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,163 @@
% Maximum historic growth for PV
% Fit to logistic curve

clear variables
close all

tic

%% loading and preparing data
% time series for energy demand and CO2 emissions

EROI = 20; %- (energy return on investment)
t_L = 30; %a (lifetime)
dt = .01; %a (time step)
EPBT = t_L/EROI; % a (energy payback time)
alpha_step=0.01;
alpha = 0:alpha_step:1; % - (continuous replacement of fossil in the economy)
n_alpha = size(alpha,2);
tau(1,1,:) = t_L./((EROI).*(1-alpha)); % adjusted doubling time constant

P_demand = 6; %TW (in 2018 without renewable fraction, which doesn't need to be replaced)

beta_step = 0.001;
beta = (beta_step:beta_step:0.4)';
n_i = size(beta,1);


time_frame = xlsread('CO2 emissions and remaining budget.xlsx','RE','A6:A60'); %a
t_ext = 1980:1:2100;

P_PV = xlsread('CO2 emissions and remaining budget.xlsx','RE','K6:K60'); %TW
P_PV (isnan(P_PV))=0;

ATP_PV = 70; %TW

%% Define functions and perform non-linear regression
% logistic
logistic_curve = @(b_log,t) ATP_PV./(1+exp(-b_log(1).*(t-b_log(2))));

PV_fit = fitnlm(time_frame,P_PV, logistic_curve,[1 2020]);
b_log_estimates = PV_fit.Coefficients.Estimate;
R2_PV = PV_fit.Rsquared.Adjusted;

% exponential
exponential_curve = @(b,t) b(2) .* exp(b(3).*(t-b(1)));

PV_fit_exp = fitnlm(time_frame,P_PV, exponential_curve,[1980 10^(-7) 0.25]);
b_estimates_exp = PV_fit_exp.Coefficients.Estimate;
R2_PV_exp = PV_fit_exp.Rsquared.Adjusted;



%% Search for alpha and beta produce historically fitted growth rate

P_invest = beta.*P_demand; %W
t_PV = 2021:dt:2100; %a
n_t = size(t_PV,2);

P = zeros(n_i,n_t,n_alpha);

for i=1:n_i
P(i,1,:) = 0;
P(i,2,:) = P_invest(i,1) * (t_PV(1,2)-t_PV(1,1)) / EPBT;
for j=3:n_t
for k=1:n_alpha
P(i,j,k) = P(i,j-1,k) + P(i,2,k) .* 2.^((t_PV(1,j-1)-t_PV(1,1))./tau(1,1,k));

end
end
end

% correction with EoL PV
P(:,(t_L/dt+1):n_t,:) = P(:,(t_L/dt+1):n_t,:) - P(:,1:(n_t - t_L/dt),:);

P_rep = zeros(n_i,n_t,n_alpha);
for i=1:n_i
for j=1:n_t
for k=1:n_alpha
if P(i,j,k)>ATP_PV
P(i,j,k)=ATP_PV;
end
if alpha(1,k) .* P (i,j,k) >= P_demand
P_rep(i,j,k) = P_demand;
else
P_rep(i,j,k) = alpha(1,k) .* P (i,j,k);
end
end
end
end


exp_curve_cap = exponential_curve(b_estimates_exp,t_PV);
for j=1:n_t
if exp_curve_cap(1,j)>ATP_PV
exp_curve_cap(1,j)=ATP_PV;
end
end

%% Finding best fit

t_fit = [2021:dt:2037];
n_t_fit = size(t_fit,2);

err_logistic_P = zeros(n_i, n_alpha);
err_exp_P = zeros(n_i, n_alpha);
for i=1:n_i
for k=1:n_alpha
err_logistic_P(i,k) = immse((logistic_curve(b_log_estimates,t_fit)-logistic_curve(b_log_estimates,t_PV(1,1))),P(i,1:n_t_fit,k));
err_exp_P(i,k) = immse((exp_curve_cap(1,1:n_t_fit)-exp_curve_cap(1,1)),P(i,1:n_t_fit,k));
end
end
[err_log_best_fit_beta,n_beta_fit_log] = min(err_logistic_P,[],1);
[err_log_best_fit,n_alpha_fit_log] = min(err_log_best_fit_beta,[],2);
beta_fit_log = beta(n_beta_fit_log(n_alpha_fit_log));
alpha_fit_log = alpha(n_alpha_fit_log);

[err_exp_best_fit_beta,n_beta_fit_exp] = min(err_exp_P,[],1);
[err_exp_best_fit,n_alpha_fit_exp] = min(err_exp_best_fit_beta,[],2);
beta_fit_exp = beta(n_beta_fit_exp(n_alpha_fit_exp));
alpha_fit_exp = alpha(n_alpha_fit_exp);


%% Calculating CO2 emissions for best fit (log)

dot_m_CO2_current = 35; %Gt/a



for k=1:n_alpha_fit_log
for i=1:n_beta_fit_log(n_alpha_fit_log)
for j=1:n_t
if P(i,j,k) > 2 * P_demand
t_transition_fit = t_PV(1,j);
break
end

end
if P(i,end,k) <= 2 * P_demand
t_transition_fit = t_PV(1,end);
end
end
end


% CO2 emissions
for i=1:n_beta_fit_log(n_alpha_fit_log)
for j=1:n_t
for k=1:n_alpha_fit_log
if t_PV(1,j) > t_transition_fit
dot_m_CO2(1,j) = 0;
else
dot_m_CO2(1,j) = dot_m_CO2_current * (1 + beta(i,1) - P_rep (i,j,k)/ P_demand);
end
end
end
end


m_CO2 = sum(dot_m_CO2 .* dt,2) + dot_m_CO2_current/0.8246 * 4; % also land use and other emissions stay constant during time delay


%%
toc
Loading

0 comments on commit a554ae6

Please sign in to comment.