Last active
February 28, 2019 10:24
-
-
Save Libbum/194ffdfae7fc169a5e106d26c761c740 to your computer and use it in GitHub Desktop.
Julia 1.1 / JuMP 0.19 Parameter issue
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 | |
# 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) | |
gσ = 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]); # 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 | |
optimize!(model); | |
photel = value.(CPRICE); | |
for i in 1:N | |
if i <= tnopol | |
JuMP.set_upper_bound(CPRICE[i], max(photel[i],cpricebase[i])); | |
end | |
end | |
set_value(pψ₂, ψ₂); | |
optimize!(model); | |
println(value(UTILITY)); |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment