diff --git a/README.md b/README.md index cccf5a8..f476c5a 100644 --- a/README.md +++ b/README.md @@ -34,9 +34,9 @@ Since there are a number of perfectly capable open source non-linear solvers in **In testing phase** - DICE-CJL +- v2016R3 **Planned** -- v2016R3 - van der Ploeg safe carbon budget Suggestions and additions welcomed. diff --git a/src/2016.jl b/src/2016.jl index 086f35f..34a970d 100644 --- a/src/2016.jl +++ b/src/2016.jl @@ -70,6 +70,15 @@ end ψ₁₀::Float64 #Initial damage intercept end +@extend struct OptionsV2016R3 <: Options + e₀::Float64 #Industrial emissions 2015 (GtCO2 per year) + μ₀::Float64 #Initial emissions control rate for base case 2015 + tnopol::Float64 #Period before which no emissions controls base + cprice₀::Float64 #Initial base carbon price (2010$ per tCO2) + gcprice::Float64 #Growth rate of base carbon price per year + ψ₁₀::Float64 #Initial damage intercept +end + @extend struct VariablesV2016 <: Variables Eind::Array{JuMP.VariableRef,1} # Industrial emissions (GtCO2 per year) @@ -89,3 +98,4 @@ include("Results2016.jl") include("2016R.jl") include("2016R2.jl") +include("2016R3.jl") diff --git a/src/2016R3.jl b/src/2016R3.jl new file mode 100644 index 0000000..ddd491a --- /dev/null +++ b/src/2016R3.jl @@ -0,0 +1,332 @@ +struct V2016R3 <: Version end + +Base.show(io::IO, v::V2016R3) = print(io, "v2016R3") + +""" + v2016R3() + +Identifier for the 2016R3 version of the model, which extends the v2016R2 version with small alterations discussed in a footnote of Nordhaus, William. 2019. "Climate Change: The Ultimate Challenge for Economics." American Economic Review, 109 (6): 1991-2014. doi:10.1257/aer.109.6.1991 + +# Examples +```jldoctest +julia> v2016R3() +v2016R3 +``` +""" +function v2016R3() + V2016R3() +end + +export v2016R3 + +function options(version::V2016R3; + N::Int = 100, #Number of years to calculate (from 2015 onwards) + tstep::Int = 5, #Years per Period + α::Float64 = 1.45, #Elasticity of marginal utility of consumption + ρ::Float64 = 0.015, #Initial rate of social time preference per year + γₑ::Float64 = 0.3, #Capital elasticity in production function + pop₀::Int = 7403, #Initial world population 2015 (millions) + popadj::Float64 = 0.134, #Growth rate to calibrate to 2050 pop projection + popasym::Int = 11500, #Asymptotic population (millions) + δk::Float64 = 0.1, #Depreciation rate on capital (per year) + q₀::Float64 = 126.6, #Initial world gross output 2018 (trill 2018 USD) + k₀::Float64 = 267.6, #Initial capital value 2018 (trill 2018 USD) + a₀::Float64 = 5.115, #Initial level of total factor productivity + ga₀::Float64 = 0.076, #Initial growth rate for TFP per 5 years + δₐ::Float64 = 0.005, #Decline rate of TFP per 5 years + gσ₁::Float64 = -0.0152, #Initial growth of sigma (continuous per year) + δσ::Float64 = -0.001, #Decline rate of decarbonization per period + eland₀::Float64 = 2.6, #Carbon emissions from land 2015 (GtCO2 per year) + deland::Float64 = 0.115, #Decline rate of land emissions (per period) + e₀::Float64 = 35.85, #Industrial emissions 2015 (GtCO2 per year) + μ₀::Float64 = 0.03, #Initial emissions control rate for base case 2015 + mat₀::Float64 = 851.0, #Initial Concentration in atmosphere 2015 (GtC) + mu₀::Float64 = 460.0, #Initial Concentration in upper strata 2015 (GtC) + ml₀::Float64 = 1740.0, #Initial Concentration in lower strata 2015 (GtC) + mateq::Float64 = 588.0, #Equilibrium concentration atmosphere (GtC) + mueq::Float64 = 360.0, #Equilibrium concentration in upper strata (GtC) + mleq::Float64 = 1720.0, #Equilibrium concentration in lower strata (GtC) + ϕ₁₂::Float64 = 0.12, #Carbon cycle transition matrix coefficient + ϕ₂₃::Float64 = 0.007, #Carbon cycle transition matrix coefficient + t2xco2::Float64 = 3.1, #Equilibrium temp impact (oC per doubling CO2) + fₑₓ0::Float64 = 0.5, #2015 forcings of non-CO2 GHG (Wm-2) + fₑₓ1::Float64 = 1.0, #2100 forcings of non-CO2 GHG (Wm-2) + tocean₀::Float64 = 0.0068, #Initial lower stratum temp change (C from 1900) + tatm₀::Float64 = 0.85, #Initial atmospheric temp change (C from 1900) + ξ₁::Float64 = 0.1005, #Climate equation coefficient for upper level + ξ₃::Float64 = 0.088, #Transfer coefficient upper to lower stratum + ξ₄::Float64 = 0.025, #Transfer coefficient for lower level + η::Float64 = 3.6813, #Forcings of equilibrium CO2 doubling (Wm-2) + ψ₁₀::Float64 = 0.0, #Initial damage intercept + ψ₁::Float64 = 0.0, #Damage intercept + ψ₂::Float64 = 0.00236, #Damage quadratic term + ψ₃::Float64 = 2.0, #Damage exponent + θ₂::Float64 = 2.6, #Exponent of control cost function + pback::Float64 = 660.0, #Cost of backstop 2018$ per tCO2 2015 + gback::Float64 = 0.025, #Initial cost decline backstop cost per period + limμ::Float64 = 1.0, #Upper limit on control rate after 2150 + tnopol::Float64 = 40.0, #Period before which no emissions controls base + cprice₀::Float64 = 2.4, #Initial base carbon price (2018$ per tCO2) + gcprice::Float64 = 0.02, #Growth rate of base carbon price per year + fosslim::Float64 = 6000.0, #Maximum cumulative extraction fossil fuels (GtC) + scale1::Float64 = 0.0302455265681763, #Multiplicative scaling coefficient + scale2::Float64 = -10993.704) #Additive scaling coefficient + OptionsV2016R3(N,tstep,α,ρ,γₑ,pop₀,popadj,popasym,δk,q₀,k₀,a₀,ga₀,δₐ,gσ₁,δσ,eland₀,deland,mat₀,mu₀,ml₀,mateq,mueq,mleq,ϕ₁₂,ϕ₂₃,t2xco2,fₑₓ0,fₑₓ1,tocean₀,tatm₀,ξ₁,ξ₃,ξ₄,η,ψ₁,ψ₂,ψ₃,θ₂,pback,gback,limμ,fosslim,scale1,scale2,e₀,μ₀,tnopol,cprice₀,gcprice,ψ₁₀) +end + +function Base.show(io::IO, ::MIME"text/plain", opt::OptionsV2016R3) + println(io, "Options for 2016R3 version"); + println(io, "Time step"); + println(io, "N: $(opt.N), tstep: $(opt.tstep)"); + println(io, "Preferences"); + println(io, "α: $(opt.α), ρ: $(opt.ρ)"); + println(io, "Population and Technology"); + println(io, "γₑ: $(opt.γₑ), pop₀: $(opt.pop₀), popadj: $(opt.popadj), popasym: $(opt.popasym), δk: $(opt.δk)"); + println(io, "q₀: $(opt.q₀), k₀: $(opt.k₀), a₀: $(opt.a₀), ga₀: $(opt.ga₀), δₐ: $(opt.δₐ)"); + println(io, "Emissions Parameters"); + println(io, "gσ₁: $(opt.gσ₁), δσ: $(opt.δσ), eland₀: $(opt.eland₀)"); + println(io, "deland: $(opt.deland), e₀: $(opt.e₀), μ₀: $(opt.μ₀)"); + println(io, "Carbon Cycle"); + println(io, "mat₀: $(opt.mat₀), mu₀: $(opt.mu₀), ml₀: $(opt.ml₀)"); + println(io, "mateq: $(opt.mateq), mueq: $(opt.mueq), mleq: $(opt.mleq)"); + println(io, "Flow Parameters"); + println(io, "ϕ₁₂: $(opt.ϕ₁₂), ϕ₂₃: $(opt.ϕ₂₃)"); + println(io, "Climate Model Parameters"); + println(io, "t2xco2: $(opt.t2xco2), fₑₓ0: $(opt.fₑₓ0), fₑₓ1: $(opt.fₑₓ1)"); + println(io, "tocean₀: $(opt.tocean₀), tatm₀: $(opt.tatm₀), ξ₁: $(opt.ξ₁)"); + println(io, "ξ₃: $(opt.ξ₃), ξ₄: $(opt.ξ₄), η: $(opt.η)"); + println(io, "Climate Damage Parameters"); + println(io, "ψ₁₀: $(opt.ψ₁₀), ψ₁: $(opt.ψ₁), ψ₂: $(opt.ψ₂), ψ₃: $(opt.ψ₃)"); + println(io, "Abatement Cost"); + println(io, "θ₂: $(opt.θ₂), pback: $(opt.pback), gback: $(opt.gback), limμ: $(opt.limμ)"); + println(io, "tnopol: $(opt.tnopol), cprice₀: $(opt.cprice₀), gcprice: $(opt.gcprice)"); + println(io, "Fossil Fuel Availability"); + println(io, "fosslim: $(opt.fosslim)"); + println(io, "Scaling Parameters"); + print(io, "scale1: $(opt.scale1), scale2: $(opt.scale2)"); +end + +function generate_parameters(c::OptionsV2016R3, model::JuMP.Model) + optlrsav::Float64 = (c.δk + .004)/(c.δk + .004c.α + c.ρ)*c.γₑ; # Optimal savings rate + ϕ₁₁::Float64 = 1 - c.ϕ₁₂; # Carbon cycle transition matrix coefficient + ϕ₂₁::Float64 = c.ϕ₁₂*c.mateq/c.mueq; # Carbon cycle transition matrix coefficient + ϕ₂₂::Float64 = 1 - ϕ₂₁ - c.ϕ₂₃; # Carbon cycle transition matrix coefficient + ϕ₃₂::Float64 = c.ϕ₂₃*c.mueq/c.mleq; # Carbon cycle transition matrix coefficient + ϕ₃₃::Float64 = 1 - ϕ₃₂; # Carbon cycle transition matrix coefficient + σ₀::Float64 = c.e₀/(c.q₀*(1-c.μ₀)); # Carbon intensity 2015 (kgCO2 per output 2018 USD) + λ::Float64 = c.η/c.t2xco2; # Climate model parameter + + @NLparameter(model, ψ₂ == c.ψ₂); + + # Backstop price + pbacktime = Array{Float64}(undef, c.N); + # Growth rate of productivity from 0 to N + gₐ = Array{Float64}(undef, c.N); + # Emissions from deforestation + Etree = Array{Float64}(undef, c.N); + # Average utility social discount rate + rr = Array{Float64}(undef, c.N); + # Carbon price in base case + cpricebase = Array{Float64}(undef, c.N); + + for i in 1:c.N + pbacktime[i] = c.pback*(1-c.gback)^(i-1); + gₐ[i] = c.ga₀*exp.(-c.δₐ*5*(i-1)); + Etree[i] = c.eland₀*(1-c.deland)^(i-1); + rr[i] = 1/((1+c.ρ)^(c.tstep*(i-1))); + cpricebase[i] = c.cprice₀*(1+c.gcprice)^(5*(i-1)); + end + + # Initial conditions and offset required + # Level of population and labor + L = Array{Float64}(undef, c.N); + L[1] = c.pop₀; + # Level of total factor productivity + A = Array{Float64}(undef, c.N); + A[1] = c.a₀; + # Change in sigma (cumulative improvement of energy efficiency) + gσ = Array{Float64}(undef, c.N); + gσ[1] = c.gσ₁; + # CO2-equivalent-emissions output ratio + σ = Array{Float64}(undef, c.N); + σ[1] = σ₀; + + # Cumulative from land + cumtree = Array{Float64}(undef, c.N); + cumtree[1] = 100.0; + + for i in 1:c.N-1 + L[i+1] = L[i]*(c.popasym/L[i])^c.popadj; + A[i+1] = A[i]/(1-gₐ[i]); + gσ[i+1] = gσ[i]*((1+c.δσ)^c.tstep); + σ[i+1] = σ[i]*exp.(gσ[i]*c.tstep); + cumtree[i+1] = cumtree[i]+Etree[i]*(5/3.666); + end + + # Adjusted cost for backstop + θ₁ = Array{Float64}(undef, c.N); + # Exogenous forcing for other greenhouse gases + fₑₓ = Array{Float64}(undef, c.N); + + for i in 1:c.N + θ₁[i] = pbacktime[i]*σ[i]/c.θ₂/1000.0; + fₑₓ[i] = if i < 18 + c.fₑₓ0+(1/17)*(c.fₑₓ1-c.fₑₓ0)*(i-1) + else + c.fₑₓ1 + end; + end + ParametersV2016(ϕ₁₁,ϕ₂₁,ϕ₂₂,ϕ₃₂,ϕ₃₃,σ₀,λ,gₐ,Etree,L,A,gσ,σ,θ₁,fₑₓ,pbacktime,cpricebase,rr,optlrsav,ψ₂,cumtree) +end + +#NOTE: CCATOT and the Tₐₜ upper bound is the only difference to the 2013R models. +function model_vars(version::V2016R3, model::JuMP.Model, N::Int64, cca_ubound::Float64, μ_ubound::Array{Float64,1}, cprice_ubound::Array{Float64,1}) + # Variables # + @variable(model, 0.0 <= μ[i=1:N] <= μ_ubound[i]); # Emission control rate GHGs #TODO: CHANGE: Allow negative emissions. + @variable(model, FORC[1:N]); # Increase in radiative forcing (watts per m2 from 1900) + @variable(model, 0.0 <= Tₐₜ[1:N] <= 20.0); # Increase temperature of atmosphere (degrees C from 1900) -- Changed from 2016R + @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 2018 US dollars per year) + @variable(model, K[1:N] >= 1.0); # Capital stock (trillions 2018 US dollars) + @variable(model, CPC[1:N] >= 0.01); # Per capita consumption (thousands 2018 USD per year) + @variable(model, I[1:N] >= 0.0); # Investment (trillions 2018 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 2018 USD per year) + @variable(model, YGROSS[1:N] >= 0.0); # Gross world product GROSS of abatement and damages (trillions 2018 USD per year) + @variable(model, YNET[1:N]); # Output net of damages equation (trillions 2018 USD per year) + @variable(model, DAMAGES[1:N]); # Damages (trillions 2018 USD per year) + @variable(model, Ω[1:N]); # Damages as fraction of gross output + @variable(model, Λ[1:N]); # Cost of emissions reductions (trillions 2018 USD per year) + @variable(model, MCABATE[1:N]); # Marginal cost of abatement (2018$ per ton CO2) + @variable(model, CCA[1:N] <= cca_ubound); # 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] <= cprice_ubound[i]); # Carbon price (2018$ per ton of CO2) + @variable(model, CEMUTOTPER[1:N]); # Period utility + @variable(model, UTILITY); # Welfare function + VariablesV2016(μ,FORC,Tₐₜ,Tₗₒ,Mₐₜ,Mᵤₚ,Mₗₒ,E,C,K,CPC,I,S,RI,Y,YGROSS,YNET,DAMAGES,MCABATE,CCA,PERIODU,UTILITY,Eind,Ω,Λ,CPRICE,CEMUTOTPER,CCATOT) +end + +function model_eqs(scenario::Scenario, model::JuMP.Model, config::OptionsV2016R3, params::ParametersV2016, vars::VariablesV2016, isMumps::Bool) + N = config.N; + # Equations # + # Emissions Equation + eeq = @constraint(model, [i=1:N], vars.E[i] == vars.Eind[i] + params.Etree[i]); + # Industrial Emissions + @constraint(model, [i=1:N], vars.Eind[i] == params.σ[i] * vars.YGROSS[i] * (1-vars.μ[i])); + # Cumulative total carbon emissions + @constraint(model, [i=1:N], vars.CCATOT[i] == vars.CCA[i] + params.cumtree[i]); + # Radiative forcing equation + @NLconstraint(model, [i=1:N], vars.FORC[i] == config.η * (log(vars.Mₐₜ[i]/588.0)/log(2)) + params.fₑₓ[i]); + # Equation for damage fraction + @NLconstraint(model, [i=1:N], vars.Ω[i] == config.ψ₁*vars.Tₐₜ[i]+params.ψ₂*vars.Tₐₜ[i]^config.ψ₃); + # Damage equation + @constraint(model, [i=1:N], vars.DAMAGES[i] == vars.YGROSS[i]*vars.Ω[i]); + # Cost of exissions reductions equation + @NLconstraint(model, [i=1:N], vars.Λ[i] == vars.YGROSS[i] * params.θ₁[i] * vars.μ[i]^config.θ₂); + # Equation for MC abatement + @NLconstraint(model, [i=1:N], vars.MCABATE[i] == params.pbacktime[i] * vars.μ[i]^(config.θ₂-1)); + # Carbon price equation from abatement + @NLconstraint(model, [i=1:N], vars.CPRICE[i] == params.pbacktime[i] * vars.μ[i]^(config.θ₂-1)); + # Output gross equation + @NLconstraint(model, [i=1:N], vars.YGROSS[i] == params.A[i]*(params.L[i]/1000.0)^(1-config.γₑ)*vars.K[i]^config.γₑ); + # Output net of damages equation + @constraint(model, [i=1:N], vars.YNET[i] == vars.YGROSS[i]*(1-vars.Ω[i])); + # Output net equation + @constraint(model, [i=1:N], vars.Y[i] == vars.YNET[i] - vars.Λ[i]); + # Consumption equation + cc = @constraint(model, [i=1:N], vars.C[i] == vars.Y[i] - vars.I[i]); + # Per capita consumption definition + @constraint(model, [i=1:N], vars.CPC[i] == 1000.0 * vars.C[i] / params.L[i]); + # Savings rate equation + @constraint(model, [i=1:N], vars.I[i] == vars.S[i] * vars.YGROSS[i]); # Altered from 2016R + # Period utility + @constraint(model, [i=1:N], vars.CEMUTOTPER[i] == vars.PERIODU[i] * params.L[i] * params.rr[i]); + # Instantaneous utility function equation + @NLconstraint(model, [i=1:N], vars.PERIODU[i] == ((vars.C[i]*1000.0/params.L[i])^(1-config.α)-1)/(1-config.α)-1); + + # Equations (offset) # + # Cumulative carbon emissions + @constraint(model, [i=1:N-1], vars.CCA[i+1] == vars.CCA[i] + vars.Eind[i]*(5/3.666)); + # Atmospheric concentration equation + @constraint(model, [i=1:N-1], vars.Mₐₜ[i+1] == vars.Mₐₜ[i]*params.ϕ₁₁ + vars.Mᵤₚ[i]*params.ϕ₂₁ + vars.E[i]*(5/3.666)); + # Lower ocean concentration + @constraint(model, [i=1:N-1], vars.Mₗₒ[i+1] == vars.Mₗₒ[i]*params.ϕ₃₃ + vars.Mᵤₚ[i]*config.ϕ₂₃); + # Shallow ocean concentration + @constraint(model, [i=1:N-1], vars.Mᵤₚ[i+1] == vars.Mₐₜ[i]*config.ϕ₁₂ + vars.Mᵤₚ[i]*params.ϕ₂₂ + vars.Mₗₒ[i]*params.ϕ₃₂); + # Temperature-climate equation for atmosphere + @constraint(model, [i=1:N-1], vars.Tₐₜ[i+1] == vars.Tₐₜ[i] + config.ξ₁*((vars.FORC[i+1]-params.λ*vars.Tₐₜ[i])-(config.ξ₃*(vars.Tₐₜ[i]-vars.Tₗₒ[i])))); + # Temperature-climate equation for lower oceans + @constraint(model, [i=1:N-1], vars.Tₗₒ[i+1] == vars.Tₗₒ[i] + config.ξ₄*(vars.Tₐₜ[i]-vars.Tₗₒ[i])); + # Capital balance equation + @constraint(model, [i=1:N-1], vars.K[i+1] <= (1-config.δk)^config.tstep * vars.K[i] + config.tstep*vars.I[i]); + # Interest rate equation + @NLconstraint(model, [i=1:N-1], vars.RI[i] == (1+config.ρ)*(vars.CPC[i+1]/vars.CPC[i])^(config.α/config.tstep)-1); + + #CHANGE: μ limit is 20% more than last step + @constraint(model, [i=2:N], vars.μ[i] <= vars.μ[i-1] + vars.μ[i-1]*0.2); + + # Savings rate for asympotic equilibrium + for i=N-9:N + JuMP.fix(vars.S[i], params.optlrsav; force=true); + end + # Initial conditions + JuMP.fix(vars.CCA[1], 400.0; force=true); + JuMP.fix(vars.Mₐₜ[1], config.mat₀; force=true); + JuMP.fix(vars.Mᵤₚ[1], config.mu₀; force=true); + JuMP.fix(vars.Mₗₒ[1], config.ml₀; force=true); + JuMP.fix(vars.Tₗₒ[1], config.tocean₀; force=true); + + if isMumps && typeof(scenario) <: OptimalPriceScenario + @constraint(model, vars.Tₐₜ[1] == config.tatm₀); + JuMP.fix(vars.K[1], config.k₀; force=true); + elseif isMumps + JuMP.fix(vars.Tₐₜ[1], config.tatm₀; force=true); + @NLconstraint(model, vars.K[1] == config.k₀); + else + JuMP.fix(vars.K[1], config.k₀; force=true); + JuMP.fix(vars.Tₐₜ[1], config.tatm₀; force=true); + end + + @constraint(model, vars.UTILITY == config.tstep * config.scale1 * sum(vars.CEMUTOTPER[i] for i=1:N) + config.scale2); + + # Objective function + @objective(model, Max, vars.UTILITY); + + assign_scenario(scenario, model, config, params, vars); + + EquationsV2016(eeq,cc) +end + +function solve(scenario::Scenario, version::V2016R3; + config::OptionsV2016R3 = options(version), + linear_solver::ipoptLinearSolver=ma97, + optimizer::JuMP.OptimizerFactory = with_optimizer(Ipopt.Optimizer, print_level=5, max_iter=99900,print_frequency_iter=250,sb="yes",linear_solver=selectLinearSolver(linear_solver))) + + # Generate a solver test to implement DICE.jl#35 hacks. + isMumps = optimizer.kwargs[:linear_solver] == "mumps"; + + model = Model(optimizer); + + params = generate_parameters(config, model); + + # Rate limit + μ_ubound = fill(config.limμ, config.N); + cprice_ubound = fill(Inf, config.N); + + variables = model_vars(version, model, config.N, config.fosslim, μ_ubound, cprice_ubound); + + equations = model_eqs(scenario, model, config, params, variables, isMumps); + + optimize!(model); + + results = model_results(model, config, params, variables, equations); + + DICENarrative(config,params,model,scenario,version,variables,equations,results) +end