Compare Julia 0.64 / JuMP 0.17 to Julia 1.1 / JuMP 0.18

dice2016Rbeta.jl

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=3, max_iter=99900,print_frequency_iter=50,sb="yes");
model = Model(optimizer);

@NLparameter(model, pψ₂ == ψ₂);
μ_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[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
# We add the first optimize! since our run is infeasible without it.
optimize!(model);
set_value(pψ₂, 0.0);
optimize!(model);

photel = value.(CPRICE);

set_value(pψ₂, ψ₂);
for i in 1:N
    if i <= tnopol
        JuMP.set_upper_bound(CPRICE[i], max(photel[i],cpricebase[i]));
    end
end

optimize!(model);
optimize!(model);
optimize!(model);