Skip to content

Instantly share code, notes, and snippets.

@Libbum
Last active Feb 28, 2019
Embed
What would you like to do?
Julia 1.1 / JuMP 0.19 Parameter issue
using JuMP
using Ipopt
N = 100 #Number of years to calculate (from 2015 onwards)
tstep = 5 #Years per Period
α = 1.45 #Elasticity of marginal utility of consumption
ρ = 0.015 #Initial rate of social time preference per year ρ
γₑ = 0.3 #Capital elasticity in production function
pop₀ = 7403 #Initial world population 2015 (millions)
popadj = 0.134 #Growth rate to calibrate to 2050 pop projection
popasym = 11500 #Asymptotic population (millions)
δk = 0.1 #Depreciation rate on capital (per year)
q₀ = 105.5 #Initial world gross output 2015 (trill 2010 USD)
k₀ = 223.0 #Initial capital value 2015 (trill 2010 USD)
a₀ = 5.115 #Initial level of total factor productivity
ga₀ = 0.076 #Initial growth rate for TFP per 5 years
δₐ = 0.005 #Decline rate of TFP per 5 years
gσ₁ = -0.0152 #Initial growth of sigma (continuous per year)
δσ = -0.001 #Decline rate of decarbonization per period
eland₀ = 2.6 #Carbon emissions from land 2015 (GtCO2 per year)
deland = 0.115 #Decline rate of land emissions (per period)
e₀ = 35.85 #Industrial emissions 2015 (GtCO2 per year)
μ₀ = 0.03 #Initial emissions control rate for base case 2015
mat₀ = 851.0 #Initial Concentration in atmosphere 2015 (GtC)
mu₀ = 460.0 #Initial Concentration in upper strata 2015 (GtC)
ml₀ = 1740.0 #Initial Concentration in lower strata 2015 (GtC)
mateq = 588.0 #Equilibrium concentration atmosphere (GtC)
mueq = 360.0 #Equilibrium concentration in upper strata (GtC)
mleq = 1720.0 #Equilibrium concentration in lower strata (GtC)
ϕ₁₂ = 0.12 #Carbon cycle transition matrix coefficient
ϕ₂₃ = 0.007 #Carbon cycle transition matrix coefficient
t2xco2 = 3.1 #Equilibrium temp impact (oC per doubling CO2)
fₑₓ0 = 0.5 #2015 forcings of non-CO2 GHG (Wm-2)
fₑₓ1 = 1.0 #2100 forcings of non-CO2 GHG (Wm-2)
tocean₀ = 0.0068 #Initial lower stratum temp change (C from 1900)
tatm₀ = 0.85 #Initial atmospheric temp change (C from 1900)
ξ₁ = 0.1005 #Climate equation coefficient for upper level
ξ₃ = 0.088 #Transfer coefficient upper to lower stratum
ξ₄ = 0.025 #Transfer coefficient for lower level
η = 3.6813 #Forcings of equilibrium CO2 doubling (Wm-2)
ψ₁₀ = 0.0 #Initial damage intercept
ψ₁ = 0.0 #Damage intercept
ψ₂ = 0.00236 #Damage quadratic term
ψ₃ = 2.0 #Damage exponent
θ₂ = 2.6 #Exponent of control cost function
pback = 550.0 #Cost of backstop 2010$ per tCO2 2015
gback = 0.025 #Initial cost decline backstop cost per period
limμ = 1.2 #Upper limit on control rate after 2150
tnopol = 45.0 #Period before which no emissions controls base
cprice₀ = 2.0 #Initial base carbon price (2010$ per tCO2)
gcprice = 0.02 #Growth rate of base carbon price per year
fosslim = 6000.0 #Maximum cumulative extraction fossil fuels (GtC)
scale1 = 0.0302455265681763 #Multiplicative scaling coefficient
scale2 = -10993.704 #Additive scaling coefficient
optlrsav = (δk + .004)/(δk + .004α + ρ)*γₑ; # Optimal savings rate
ϕ₁₁ = 1 - ϕ₁₂; # Carbon cycle transition matrix coefficient
ϕ₂₁ = ϕ₁₂*mateq/mueq; # Carbon cycle transition matrix coefficient
ϕ₂₂ = 1 - ϕ₂₁ - ϕ₂₃; # Carbon cycle transition matrix coefficient
ϕ₃₂ = ϕ₂₃*mueq/mleq; # Carbon cycle transition matrix coefficient
ϕ₃₃ = 1 - ϕ₃₂; # Carbon cycle transition matrix coefficient
σ₀ = e₀/(q₀*(1-μ₀)); # Carbon intensity 2010 (kgCO2 per output 2010 USD 2010)
λ = η/t2xco2; # Climate model parameter
#Temporary. This is usually found by running the optimizer with ψ₂ = 0 first.
photel = [0.0010013987556195469, 0.001274195263657956, 0.0016161080142431494, 0.0020427411864263165, 0.0025727653717606902, 0.0032284505332554046, 0.0040362749478795815, 0.005027619619829979, 0.006239558473165778, 0.007715755523683413, 0.009507481201467405, 0.011674761052498477, 0.014287671194145818, 0.017427796139995046, 0.021189865949953014, 0.02568359110870239, 0.031035715097203386, 0.037392306306965245, 0.04492131272047702, 0.05381540488251007, 0.06429513463072647, 0.07661243943931612, 0.09105452467766423, 0.1079481587541999, 0.1276644189988594, 0.15062392925401671, 0.17730263351273493, 0.20823815358575318, 0.2440367827173673, 0.2853814847112282, 0.3330391360391443, 0.3878715026747601, 0.4508442861794069, 0.5230385027625859, 0.6056626474143237, 0.7000659982244626, 0.8077531583394961, 0.9303999409657192, 1.069870711428674, 1.2282373095937427, 1.4077996859955548, 1.611108395873057, 1.8409891070250737, 2.100569290055184, 2.393307273236377, 2.7230238589720073, 3.0939367147418646, 3.5106977685876255, 3.9784338577028073, 4.502790898645537, 5.089981869189603, 5.746838914975056, 6.480869919020307, 7.300319898918406, 8.214237625258482, 9.232547885572727, 10.366129850962892, 11.626902037508703, 13.027914391509043, 14.583448066294036, 16.309123498209917, 18.222017429308945, 20.340789562126012, 22.685819563509114, 25.279355151422404, 28.14567198351224, 31.31124598070134, 34.80493847353268, 38.658193927780374, 42.90524836913395, 47.58334209164315, 52.73291618924306, 58.397721243919996, 64.62452829833369, 71.4604724255574, 78.91818638992561, 80.2902785176549, 78.29567568071172, 76.33891572404065, 74.43112349373372, 72.57111858664092, 70.7577175077947, 68.98976627793486, 67.26613840730373, 65.58573126077532, 63.947459494543594, 62.35024612429199, 60.793267344947395, 59.2940374105241, 57.76823753067119, 56.35542339380239, 54.94890912810089, 53.576797694253294, 52.23841108060393, 50.932362733226164, 49.65709368708481, 48.41077837915529, 47.185264695299566, 45.76061536153112, 0.6517967340966072];
# Backstop price
pbacktime = Array{Float64}(undef, N);
# Growth rate of productivity from 0 to N
gₐ = Array{Float64}(undef, N);
# Emissions from deforestation
Etree = Array{Float64}(undef, N);
# Average utility social discount rate
rr = Array{Float64}(undef, N);
# Carbon price in base case
cpricebase = Array{Float64}(undef, N);
for i in 1:N
pbacktime[i] = pback*(1-gback)^(i-1);
gₐ[i] = ga₀*exp.(-δₐ*5*(i-1));
Etree[i] = eland₀*(1-deland)^(i-1);
rr[i] = 1/((1+ρ)^(tstep*(i-1)));
cpricebase[i] = cprice₀*(1+gcprice)^(5*(i-1));
end
# Initial conditions and offset required
# Level of population and labor
L = Array{Float64}(undef, N);
L[1] = pop₀;
# Level of total factor productivity
A = Array{Float64}(undef, N);
A[1] = a₀;
# Change in sigma (cumulative improvement of energy efficiency)
= Array{Float64}(undef, N);
gσ[1] = gσ₁;
# CO2-equivalent-emissions output ratio
σ = Array{Float64}(undef, N);
σ[1] = σ₀;
# Cumulative from land
cumtree = Array{Float64}(undef, N);
cumtree[1] = 100.0;
for i in 1:N-1
L[i+1] = L[i]*(popasym/L[i])^popadj;
A[i+1] = A[i]/(1-gₐ[i]);
gσ[i+1] = gσ[i]*((1+δσ)^tstep);
σ[i+1] = σ[i]*exp.(gσ[i]*tstep);
cumtree[i+1] = cumtree[i]+Etree[i]*(5/3.666);
end
# Adjusted cost for backstop
θ₁ = Array{Float64}(undef, N);
# Exogenous forcing for other greenhouse gases
fₑₓ = Array{Float64}(undef, N);
for i in 1:N
θ₁[i] = pbacktime[i]*σ[i]/θ₂/1000.0;
fₑₓ[i] = if i < 18
fₑₓ0+(1/17)*(fₑₓ1-fₑₓ0)*(i-1)
else
fₑₓ1-fₑₓ0
end;
end
optimizer = with_optimizer(Ipopt.Optimizer, print_level=5, max_iter=99900,print_frequency_iter=250,sb="yes");
model = Model(optimizer);
@NLparameter(model, pψ₂ == 0.0);
μ_ubound = [if t < 30 1.0 else limμ end for t in 1:N];
# Variables #
@variable(model, 0.0 <= μ[i=1:N] <= μ_ubound[i]); # Emission control rate GHGs
@variable(model, FORC[1:N]); # Increase in radiative forcing (watts per m2 from 1900)
@variable(model, 0.0 <= Tₐₜ[1:N] <= 12.0); # Increase temperature of atmosphere (degrees C from 1900)
@variable(model, -1.0 <= Tₗₒ[1:N] <= 20.0); # Increase temperatureof lower oceans (degrees C from 1900)
@variable(model, Mₐₜ[1:N] >= 10.0); # Carbon concentration increase in atmosphere (GtC from 1750)
@variable(model, Mᵤₚ[1:N] >= 100.0); # Carbon concentration increase in shallow oceans (GtC from 1750)
@variable(model, Mₗₒ[1:N] >= 1000.0); # Carbon concentration increase in lower oceans (GtC from 1750)
@variable(model, E[1:N]); # Total CO2 emissions (GtCO2 per year)
@variable(model, Eind[1:N]); # Industrial emissions (GtCO2 per year)
@variable(model, C[1:N] >= 2.0); # Consumption (trillions 2010 US dollars per year)
@variable(model, K[1:N] >= 1.0); # Capital stock (trillions 2010 US dollars)
@variable(model, CPC[1:N] >= 0.01); # Per capita consumption (thousands 2010 USD per year)
@variable(model, I[1:N] >= 0.0); # Investment (trillions 2010 USD per year)
@variable(model, S[1:N]); # Gross savings rate as fraction of gross world product
@variable(model, RI[1:N]); # Real interest rate (per annum)
@variable(model, Y[1:N] >= 0.0); # Gross world product net of abatement and damages (trillions 2010 USD per year)
@variable(model, YGROSS[1:N] >= 0.0); # Gross world product GROSS of abatement and damages (trillions 2010 USD per year)
@variable(model, YNET[1:N]); # Output net of damages equation (trillions 2010 USD per year)
@variable(model, DAMAGES[1:N]); # Damages (trillions 2010 USD per year)
@variable(model, Ω[1:N]); # Damages as fraction of gross output
@variable(model, Λ[1:N]); # Cost of emissions reductions (trillions 2010 USD per year)
@variable(model, MCABATE[1:N]); # Marginal cost of abatement (2010$ per ton CO2)
@variable(model, CCA[1:N] <= fosslim); # Cumulative industrial carbon emissions (GTC)
@variable(model, CCATOT[1:N]); # Total carbon emissions (GTC)
@variable(model, PERIODU[1:N]); # One period utility function
@variable(model, CPRICE[i=1:N], start = photel[i]); # Carbon price (2010$ per ton of CO2)
@variable(model, CEMUTOTPER[1:N]); # Period utility
@variable(model, UTILITY); # Welfare function
# Equations #
# Emissions Equation
eeq = @constraint(model, [i=1:N], E[i] == Eind[i] + Etree[i]);
# Industrial Emissions
@constraint(model, [i=1:N], Eind[i] == σ[i] * YGROSS[i] * (1-μ[i]));
# Cumulative total carbon emissions
@constraint(model, [i=1:N], CCATOT[i] == CCA[i] + cumtree[i]);
# Radiative forcing equation
@NLconstraint(model, [i=1:N], FORC[i] == η * (log(Mₐₜ[i]/588.0)/log(2)) + fₑₓ[i]);
# Equation for damage fraction
@NLconstraint(model, [i=1:N], Ω[i] == ψ₁*Tₐₜ[i]+pψ₂*Tₐₜ[i]^ψ₃);
# Damage equation
@constraint(model, [i=1:N], DAMAGES[i] == YGROSS[i]*Ω[i]);
# Cost of exissions reductions equation
@NLconstraint(model, [i=1:N], Λ[i] == YGROSS[i] * θ₁[i] * μ[i]^θ₂);
# Equation for MC abatement
@NLconstraint(model, [i=1:N], MCABATE[i] == pbacktime[i] * μ[i]^(θ₂-1));
# Carbon price equation from abatement
@NLconstraint(model, [i=1:N], CPRICE[i] == pbacktime[i] * μ[i]^(θ₂-1));
# Output gross equation
@NLconstraint(model, [i=1:N], YGROSS[i] == A[i]*(L[i]/1000.0)^(1-γₑ)*K[i]^γₑ);
# Output net of damages equation
@constraint(model, [i=1:N], YNET[i] == YGROSS[i]*(1-Ω[i]));
# Output net equation
@constraint(model, [i=1:N], Y[i] == YNET[i] - Λ[i]);
# Consumption equation
cc = @constraint(model, [i=1:N], C[i] == Y[i] - I[i]);
# Per capita consumption definition
@constraint(model, [i=1:N], CPC[i] == 1000.0 * C[i] / L[i]);
# Savings rate equation
@constraint(model, [i=1:N], I[i] == S[i] * Y[i]);
# Period utility
@constraint(model, [i=1:N], CEMUTOTPER[i] == PERIODU[i] * L[i] * rr[i]);
# Instantaneous utility function equation
@NLconstraint(model, [i=1:N], PERIODU[i] == ((C[i]*1000.0/L[i])^(1-α)-1)/(1-α)-1);
# Equations (offset) #
# Cumulative carbon emissions
@constraint(model, [i=1:N-1], CCA[i+1] == CCA[i] + Eind[i]*(5/3.666));
# Atmospheric concentration equation
@constraint(model, [i=1:N-1], Mₐₜ[i+1] == Mₐₜ[i]*ϕ₁₁ + Mᵤₚ[i]*ϕ₂₁ + E[i]*(5/3.666));
# Lower ocean concentration
@constraint(model, [i=1:N-1], Mₗₒ[i+1] == Mₗₒ[i]*ϕ₃₃ + Mᵤₚ[i]*ϕ₂₃);
# Shallow ocean concentration
@constraint(model, [i=1:N-1], Mᵤₚ[i+1] == Mₐₜ[i]*ϕ₁₂ + Mᵤₚ[i]*ϕ₂₂ + Mₗₒ[i]*ϕ₃₂);
# Temperature-climate equation for atmosphere
@constraint(model, [i=1:N-1], Tₐₜ[i+1] == Tₐₜ[i] + ξ₁*((FORC[i+1]-λ*Tₐₜ[i])-(ξ₃*(Tₐₜ[i]-Tₗₒ[i]))));
# Temperature-climate equation for lower oceans
@constraint(model, [i=1:N-1], Tₗₒ[i+1] == Tₗₒ[i] + ξ₄*(Tₐₜ[i]-Tₗₒ[i]));
# Capital balance equation
@constraint(model, [i=1:N-1], K[i+1] <= (1-δk)^tstep * K[i] + tstep*I[i]);
# Interest rate equation
@NLconstraint(model, [i=1:N-1], RI[i] == (1+ρ)*(CPC[i+1]/CPC[i])^/tstep)-1);
# Savings rate for asympotic equilibrium
@constraint(model, S[N-10:N] .== optlrsav);
# Initial conditions
@constraint(model, CCA[1] == 400.0);
@constraint(model, K[1] == k₀);
@constraint(model, Mₐₜ[1] == mat₀);
@constraint(model, Mᵤₚ[1] == mu₀);
@constraint(model, Mₗₒ[1] == ml₀);
@constraint(model, Tₐₜ[1] == tatm₀);
@constraint(model, Tₗₒ[1] == tocean₀);
@constraint(model, UTILITY == tstep * scale1 * sum(CEMUTOTPER[i] for i=1:N) + scale2);
# Objective function
@objective(model, Max, UTILITY);
# Base
for i in 1:N
if i <= tnopol
JuMP.set_upper_bound(CPRICE[i], max(photel[i],cpricebase[i]));
end
end
optimize!(model);
println(value(UTILITY));
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment