[DICE: GAMS -> JuMP] How to debug JuMP model errors (find the culpit)?

Specific Domains Optimization (Mathematical)
, , , ,
1

Hello, this post il dual :slight_smile:
It’s at the same time a request for help for a specific case, but also a much more general question.
I am trying to port the Nordhaus’s DICE climate change model v2023 to Julia, in its simplest possible form.

This is my attempt in Julia/JuMP 
# Implementation of the DICE 2023 model
# This version has all the par values and code elements in a single file
# (just for testing.. know is not good!!)
# Based on DICE2023-b-4-3-10.gms and included files
# (Nonco2-b-4-3-1.gms and FAIR-beta-4-3-1.gms)

using Pkg
Pkg.activate(@__DIR__) # Using the attached Project.toml and Manifest-v1.11.toml guarantees reproducibility
using JuMP, Ipopt # JuMP is the algeabric modeling language, Ipopt is the solver

######################################################################
# Raw (exogenous) parameters
######################################################################

tstep   = 5  # Years per period
ntsteps = 81 # Number of time periods

# --------------------------------------------------------------------
# Population and technology

gama    = 0.300   # Capital elasticity in production function
pop1    = 7752.9  # Initial world population 2020 (millions)
popadj  = 0.145   # Growth rate to calibrate to 2050 population projection
popasym = 10825.0 # Asymptotic population (millions)
dk      = 0.100   # Depreciation rate on capital (per year)
q1      = 135.7   # Initial world output 2020 (trill 2019 USD)
al1     = 5.84    # Initial level of total factor productivity
ga1     = 0.066   # Initial growth rate for TFP per 5 years
dela    = 0.0015  # Decline rate of TFP per 5 years

# --------------------------------------------------------------------
# Emissions parameters and Non-CO2 GHG with sigma = emissions/output

gsigma1   = -0.015 # Initial growth of sigma (per year)
delgsig   = 0.96   # Decline rate of gsigma per period
asymgsig  = -0.005 # Asymptotic gsigma
e1        = 37.56  # Industrial emissions 2020 (GtCO2 per year)
miu1      = 0.05   # Emissions control rate historical 2020
fosslim   = 6000.0 # Maximum cumulative extraction fossil fuels (GtC)
cumemiss0 = 633.5  # Cumulative emissions 2020 (GtC)

# --------------------------------------------------------------------
# Climate damage parameters

a1     = 0.0      # Damage intercept
a2base = 0.003467 # Damage quadratic term
a3     = 2.0      # Damage exponent

# --------------------------------------------------------------------
# Abatement cost

expcost2  = 2.6    # Exponent of control cost function
pback2050 = 515.0  # Cost of backstop in 2019$ per tCO2 (2050)
gback     = -0.012 # Initial cost decline of backstop cost per year
cprice1   = 6.0    # Carbon price in 2020 (2019$ per tCO2)
gcprice   = 0.025  # Growth rate of base carbon price per year

# --------------------------------------------------------------------
# Limits on emissions controls

limmiu2070 = 1.0  # Emission control limit from 2070
limmiu2120 = 1.1  # Emission control limit from 2120
limmiu2200 = 1.05 # Emission control limit from 2220
limmiu2300 = 1.0  # Emission control limit from 2300
delmiumax  = 0.12 # Emission control delta limit per period

# --------------------------------------------------------------------
# Preferences, growth uncertainty, and timing

betaclim = 0.5   # Climate beta
elasmu   = 0.95  # Elasticity of marginal utility of consumption
prstp    = 0.001 # Pure rate of social time preference
pi_val   = 0.05  # Capital risk premium (renamed to avoid conflict with Julia's pi)
k0       = 295.0 # Initial capital stock (10^12 2019 USD)
siggc1   = 0.01  # Annual standard deviation of consumption growth

# --------------------------------------------------------------------
# Scaling so that MU(C(1)) = 1 and objective function = PV consumption

srf    = 1000000.0  # Scaling factor for discounting
scale1 = 0.00891061 # Multiplicative scaling coefficient
scale2 = -6275.91   # Additive scaling coefficient

# --------------------------------------------------------------------
# Parameters for non-industrial emissions

eland0         = 5.9     # Carbon emissions from land 2015 (GtCO2 per year)
deland         = 0.1     # Decline rate of land emissions (per period)
f_misc2020     = -0.054  # Non-abatable forcings 2020
f_misc2100     = 0.265   # Non-abatable forcings 2100
f_ghgabate2020 = 0.518   # Forcings of abatable non-CO2 GHG
f_ghgabate2100 = 0.957   # Forcings of abatable non-CO2 GHG

eco2eghgb2020  = 9.96    # Emissions of abatable non-CO2 GHG (GtCO2e) in 2020
eco2eghgb2100  = 15.5    # Emissions of abatable non-CO2 GHG (GtCO2e) in 2100
emissrat2020   = 1.40    # Ratio of CO2e to industrial CO2 in 2020
emissrat2100   = 1.21    # Ratio of CO2e to industrial CO2 in 2100
fcoef1         = 0.00955 # Coefficient of non-CO2 abateable emissions
fcoef2         = 0.861   # Coefficient of non-CO2 abateable emissions

# --------------------------------------------------------------------
# Parameters for the DFAIR model

yr0      = 2020.0    # Calendar year that corresponds to model year zero

emshare0 = 0.2173    # Carbon emissions share into Reservoir 0
emshare1 = 0.224     # Carbon emissions share into Reservoir 1
emshare2 = 0.2824    # Carbon emissions share into Reservoir 2
emshare3 = 0.2763    # Carbon emissions share into Reservoir 3
tau0     = 1000000.0 # Decay time constant for Reservoir 0
tau1     = 394.4     # Decay time constant for Reservoir 1
tau2     = 36.53     # Decay time constant for Reservoir 2
tau3     = 4.304     # Decay time constant for Reservoir 3

teq1     = 0.324     # Thermal equilibration parameter for box 1
teq2     = 0.44      # Thermal equilibration parameter for box 2
d1       = 236.0     # Thermal response timescale for deep ocean
d2       = 4.07      # Thermal response timescale for upper ocean

irf0     = 32.4      # Pre-industrial IRF100
irc      = 0.019     # Increase in IRF100 with cumulative carbon uptake
irt      = 4.165     # Increase in IRF100 with warming
fco22x   = 3.93      # Forcings of equilibrium CO2 doubling

# --------------------------------------------------------------------
# Initial conditions to be calibrated to history calibration

mat0   = 886.5128014 # Initial concentration in atmosphere in 2020 (GtC)

res00  = 150.093     # Initial concentration in Reservoir 0 in 2020 (GtC)
res10  = 102.698     # Initial concentration in Reservoir 1 in 2020 (GtC)
res20  = 39.534      # Initial concentration in Reservoir 2 in 2020 (GtC)
res30  = 6.1865      # Initial concentration in Reservoir 3 in 2020 (GtC)

mateq  = 588.0       # Equilibrium concentration in atmosphere (GtC)
tbox10 = 0.1477      # Initial temperature box 1 change in 2020 (°C)
tbox20 = 1.099454    # Initial temperature box 2 change in 2020 (°C)
# TODO: check if this should be computed instead (as I think)..
# Changed to prevent numerical instability (GAMS original: 1.24715)...
tatm0  = 1.247154    # Initial atmospheric temperature change in 2020 (°C)

######################################################################
# Computed parameters
######################################################################

# --------------------------------------------------------------------
# Preferences, growth uncertainty, and timing

# Time periods sequence
times = 0:tstep:(ntsteps*tstep)-1  # 0,5,10,...,400
# Time periods index sequences
tidx  = 1:ntsteps                 # 1,2,3,...,81
t0idx = 0:ntsteps-1               # 0,1,2,...,80

# Risk-adjusted rate of time preference
rartp = exp(prstp + betaclim*pi_val)-1
 
# --------------------------------------------------------------------
# Precautionary dynamic parameters

# Variance of per capita consumption 
varpcc           = [min((siggc1^2)*t,(siggc1^2)*tstep*47) for t in times]
# Precautionary rate of return   
rprecaut         = @. -0.5 * varpcc * elasmu^2
# STP factor without precautionary factor
rr1              = @. 1 / ((1+rartp)^times)
# STP with precautionary factor
rr               = @. rr1 * (1+rprecaut) ^ -times
# Level of population and labor
l                = fill(pop1, ntsteps); [l[ti] = l[ti-1] * (popasym / l[ti-1])^popadj for ti in tidx[2:end]]
# Growth rate of Total Factor Productivity
ga               = @. ga1*exp(-dela*times)
# Level of total factor productivity
al               = fill(al1, ntsteps); [al[ti] = al[ti-1] / (1 - ga[ti-1]) for ti in tidx[2:end]]
# Optimal long-run savings rate used for transversality
optlrsav         = (dk + 0.004)/(dk + 0.004*elasmu + rartp)*gama
# Carbon price in base case
cpricebase       = @. cprice1*(1+gcprice)^times
# Backstop price 2019$ per ton CO2
pbacktime        = @. pback2050*exp(-tstep*0.01*(tidx-7))
pbacktime[8:end] = @. pback2050*exp(-tstep*0.001*(tidx[8:end]-7))
# Carbon intensity 2020 kgCO2-output 2020
sig1             = e1/(q1*(1-miu1))

# Change in sigma (rate of decarbonization)
gsig             = @. min(gsigma1*delgsig ^(tidx-1),asymgsig)

# CO2-emissions output ratio
sigma            = fill(sig1,ntsteps)
[sigma[ti] = sigma[ti-1] * exp(tstep*gsig[ti-1]) for ti in tidx[2:end] ]

# --------------------------------------------------------------------
# Emissions limits

# Upper bounds on miu
miuup          = fill(1.0, ntsteps)
miuup[1]       = 0.05; miuup[2] = 0.10
@. miuup[3:8]  =  delmiumax*(tidx[3:8]-1)
@. miuup[9:11] =  0.85+.05*(tidx[9:11]-8)
miuup[12:20]  .= limmiu2070
miuup[21:37]  .= limmiu2120
miuup[38:57]  .= limmiu2200
miuup[58:end] .= limmiu2300

# ------------------------------------------------------------------------------
# Parameters emissions and non-CO2 
 
eland          = @. eland0*(1-deland)^t0idx  
co2e_ghgabateb = eco2eghgb2020 .+ [(eco2eghgb2100-eco2eghgb2020) * min(1,ti/16) for ti in t0idx]
f_misc         = f_misc2020    .+ [(f_misc2100-f_misc2020) * min(1,ti/16) for ti in t0idx]
emissrat       = emissrat2020  .+ [(emissrat2100-emissrat2020) * min(1,ti/16) for ti in t0idx]
sigmatot       = @. sigma * emissrat
cost1tot       = @. pbacktime * sigmatot / expcost2/1000

######################################################################
# Optimization model & computation options
######################################################################

m = Model(Ipopt.Optimizer)

######################################################################
# Variables declaration
######################################################################

@variables m begin

    ECO2[tidx]             # Total CO2 emissions (GtCO2 per year)
    ECO2E[tidx]            # Total CO2e emissions including abateable nonCO2 GHG (GtCO2 per year)
    EIND[tidx]             # Industrial CO2 emissions (GtCO2 per yr)
    # TODO: check this, as in GAMS it is declared without time index, but then it is used with time index!
    F_GHGABATE[tidx]       # Forcings abateable nonCO2 GHG    
    
    MIU[tidx] >= 0.0       # Emission control rate GHGs
    C[tidx] >= 0.0         # Consumption (trillions 2019 US dollars per year)
    K[tidx] >= 0.0         # Capital stock (trillions 2019 US dollars)
    CPC[tidx]              # Per capita consumption (thousands 2019 USD per year)
    I[tidx] >= 0.0         # Investment (trillions 2019 USD per year)
    S[tidx]                # Gross savings rate as fraction of gross world product (control)
    Y[tidx] >= 0.0         # Gross world product net of abatement and damages (trillions 2019 USD per year)
    YGROSS[tidx] >= 0.0    # Gross world product GROSS of abatement and damages (trillions 2019 USD per year)
    YNET[tidx] >= 0.0      # Output net of damages equation (trillions 2019 USD per year)
    DAMAGES[tidx]          # Damages (trillions 2019 USD per year)
    DAMFRAC[tidx]          # Damages as fraction of gross output
    ABATECOST[tidx]        # Cost of emissions reductions  (trillions 2019 USD per year)
    MCABATE[tidx]          # Marginal cost of abatement (2019$ per ton CO2)
    CCATOT[tidx]           # Total carbon emissions (GtC)
    PERIODU[tidx]          # One period utility function
    CPRICE[tidx]           # Carbon price (2019$ per ton of CO2)
    TOTPERIODU[tidx]       # Period utility
    UTILITY                # Welfare function
    RFACTLONG[tidx] >= 0.0 #
    RSHORT[tidx[2:end]]    # Real interest rate with precautionary(per annum year on year)
    RLONG[tidx[2:end]]     # Real interest rate from year 0 to T 

    # --------------------------------------------------------------------
    # Climate variables
    # Note: Stock variables correspond to levels at the END of the period

    FORC[tidx]             # Increase in radiative forcing (watts per m2 from 1765)
    TATM[tidx] >= 0.0      # Increase temperature of atmosphere (degrees C from 1765)     
    TBOX1[tidx]            # Increase temperature of box 1 (degrees C from 1765)
    TBOX2[tidx]            # Increase temperature of box 2 (degrees C from 1765)
    RES0[tidx]             # Carbon concentration in Reservoir 0 (GtC from 1765)
    RES1[tidx]             # Carbon concentration in Reservoir 1 (GtC from 1765)
    RES2[tidx]             # Carbon concentration in Reservoir 2 (GtC from 1765)
    RES3[tidx]             # Carbon concentration in Reservoir 3 (GtC from 1765)
    MAT[tidx] >= 0.0       # Carbon concentration increase in atmosphere (GtC from 1765)
    CACC[tidx]             # Accumulated carbon in ocean and other sinks (GtC)
    IRFT[tidx] >= 0.0      # IRF100 at time t
    ALPHA[tidx] >= 0.0     # Carbon decay time scaling factor (control)
    #SUMALPHA              # Placeholder variable for objective function
end

######################################################################
# Constraints
######################################################################

# --------------------------------------------------------------------
# Climate variables
# Note: initial conditionsand stability bounds are set here as contraints, we can also try as fixed value

# Reservoir 0 law of motion
@constraint(m, res0lom[ti in tidx], RES0[ti] ==  ((ti == 1) ? res00 : (emshare0*tau0*ALPHA[ti]*(ECO2[ti]/3.667))*(1-exp(-tstep/(tau0*ALPHA[ti])))+RES0[ti-1]*exp(-tstep/(tau0*ALPHA[ti]))))

# Reservoir 1 law of motion
@constraint(m, res1lom[ti in tidx], RES1[ti] == ((ti == 1) ? res10 :  (emshare1*tau1*ALPHA[ti]*(ECO2[ti]/3.667))*(1-exp(-tstep/(tau1*ALPHA[ti])))+RES1[ti-1]*exp(-tstep/(tau1*ALPHA[ti]))))

# Reservoir 2 law of motion
@constraint(m, res2lom[ti in tidx], RES2[ti] ==  ((ti == 1) ? res20 : (emshare2*tau2*ALPHA[ti]*(ECO2[ti]/3.667))*(1-exp(-tstep/(tau2*ALPHA[ti])))+RES2[ti-1]*exp(-tstep/(tau2*ALPHA[ti]))))

# Reservoir 3 law of motion 
@constraint(m, res3lom[ti in tidx], RES3[ti] ==  ((ti == 1) ? res30 : (emshare3*tau3*ALPHA[ti]*(ECO2[ti]/3.667))*(1-exp(-tstep/(tau3*ALPHA[ti])))+RES3[ti-1]*exp(-tstep/(tau3*ALPHA[ti]))))

# Atmospheric concentration equation
@constraint(m, matlb[ti in tidx], MAT[ti] >= 10.0)
@constraint(m, mmat[ti in tidx], MAT[ti] == ((ti == 1) ? mat0 : mateq + RES0[ti] + RES1[ti] + RES2[ti] + RES3[ti]))

# Accumulated carbon in sinks equation
@constraint(m, cacceq[ti in tidx], CACC[ti] == CCATOT[ti]-(MAT[ti]-mateq))

# Radiative forcing equation
@constraint(m, force[ti in tidx], FORC[ti] == fco22x*((log((MAT[ti]/mateq))/log(2))) + f_misc[ti]+F_GHGABATE[ti] )

# Temperature box 1 law of motion
@constraint(m, tbox1eq[ti in tidx], TBOX1[ti] ==  ((ti == 1) ? tbox10 : TBOX1[ti-1]*exp(-tstep/d1)+teq1*FORC[ti]*(1-exp(-tstep/d1))))

# Temperature box 2 law of motion
@constraint(m, tbox2eq[ti in tidx[2:end]], TBOX2[ti] ==  ((ti == 1) ? tbox20 : TBOX2[ti-1]*exp(-tstep/d2)+teq2*FORC[ti]*(1-exp(-tstep/d2))))

# Temperature-climate equation for atmosphere
@constraint(m, tatmlb[ti in tidx], TATM[ti] >= 0.5)
@constraint(m, tatmub[ti in tidx], TATM[ti] <= 20.0)
@constraint(m, tatmeq[ti in tidx], TATM[ti] ==  ((ti == 1) ? tatm0 : TBOX1[ti]+TBOX2[ti]))

# Left-hand side of IRF100 equation
@constraint(m, irfeqlhs[ti in tidx],  IRFT[ti] == (ALPHA[ti]*emshare0*tau0*(1-exp(-100/(ALPHA[ti]*tau0))))+(ALPHA[ti]*emshare1*tau1*(1-exp(-100/(ALPHA[ti]*tau1))))+(ALPHA[ti]*emshare2*tau2*(1-exp(-100/(ALPHA[ti]*tau2))))+(ALPHA[ti]*emshare3*tau3*(1-exp(-100/(ALPHA[ti]*tau3)))))

# Right-hand side of IRF100 equation
@constraint(m, irfeqrhs[ti in tidx],  IRFT[ti] == irf0+irc*CACC[ti]+irt*TATM[ti])

# --------------------------------------------------------------------
# Emissions and Damages

# CO2 Emissions equation
@constraint(m, eco2eq[ti in tidx], ECO2[ti] == (sigma[ti]*YGROSS[ti] + eland[ti])*(1-MIU[ti]))

# Industrial CO2 equation
@constraint(m, eindeq[ti in tidx], EIND[ti] == (sigma[ti]*YGROSS[ti])*(1-MIU[ti]))

# CO2E Emissions equation
@constraint(m, eco2eeq[ti in tidx], ECO2E[ti] == (sigma[ti]*YGROSS[ti] + eland[ti] + co2e_ghgabateb[ti])*(1-MIU[ti]))

# Forcings abateable nonCO2 GHG equation 
@constraint(m, f_ghgabateeq[ti in tidx], F_GHGABATE[ti] == ((ti == 1) ? f_ghgabate2020 : fcoef2*F_GHGABATE[ti-1]+ fcoef1*co2e_ghgabateb[ti-1]*(1-MIU[ti-1])))

# --------------------------------------------------------------------
# Emissions and damage

# Cumulative total carbon emissions
@constraint(m, ccatoteq[ti in tidx], CCATOT[ti] == ((ti==1) ? cumemiss0 : CCATOT[ti-1] +  ECO2[ti-1]*(tstep/3.666)))

# Equation for damage fraction
@constraint(m, damfraceq[ti in tidx], DAMFRAC[ti] == (a1*TATM[ti])+(a2base*TATM[ti]^a3))

# Damage equation
@constraint(m, dameq[ti in tidx], DAMAGES[ti] == YGROSS[ti] * DAMFRAC[ti])

# Cost of emissions reductions equation
@constraint(m, abateeq[ti in tidx], ABATECOST[ti] == YGROSS[ti] * cost1tot[ti] * (MIU[ti]^expcost2))

# Equation for MC abatement
@constraint(m, mcabateeq[ti in tidx], MCABATE[ti] == pbacktime[ti] * MIU[ti]^(expcost2-1))

# Carbon price equation from abatement
@constraint(m, carbpriceeq[ti in tidx], CPRICE[ti] == pbacktime[ti] * MIU[ti]^(expcost2-1))

# --------------------------------------------------------------------
# Economic variables

# Output gross equation
@constraint(m, ygrosseq[ti in tidx], YGROSS[ti] == (al[ti]*(l[ti]/1000)^(1-gama))*(K[ti]^gama))

# Output net of damage equation
@constraint(m, yneteq[ti in tidx], YNET[ti] == YGROSS[ti]*(1-DAMFRAC[ti]))

# Output net equation
@constraint(m, yy[ti in tidx], Y[ti] == YNET[ti] - ABATECOST[ti])

# Consumption equation
@constraint(m, cc[ti in tidx], C[ti] == Y[ti] - I[ti])
@constraint(m, clb[ti in tidx],  C[ti] >= 2)

# Per capita consumption definition
@constraint(m, cpce[ti in tidx], CPC[ti] == 1000 * C[ti] / l[ti])
@constraint(m, cpclb[ti in tidx],  CPC[ti] >= 0.01)

# Savings rate equation
@constraint(m, seq[ti in tidx], I[ti] == S[ti] * Y[ti])

# Capital balance equation
@constraint(m, kk0, K[1] == k0)
@constraint(m, kk[ti in tidx[2:end]],  K[ti] <= (1-dk)^tstep * K[ti-1] + tstep * I[ti-1])
@constraint(m, klb[ti in tidx],  K[ti] >= 1)

# Long interest factor
@constraint(m, rfactlonglb[ti in tidx],  RFACTLONG[ti] >= 0.0001)
@constraint(m, rfactlongeq[ti in tidx],  RFACTLONG[ti] == ((ti ==1) ? 1000000 : srf *(CPC[ti]/CPC[1])^(-elasmu)*rr[ti]))

# Long-run interest rate equation
@constraint(m, rlongeq[ti in tidx[2:end]], RLONG[ti] ==  -log(RFACTLONG[ti]/srf)/(tstep*ti))

# Short-run interest rate equation
@constraint(m, rshorteq[ti in tidx[2:end]], RSHORT[ti] == -log(RFACTLONG[ti]/RFACTLONG[ti-1])/tstep)

# -------------------------------------------------------------------- 
# Utility

# Instantaneous utility function equation
@constraint(m, periodueq[ti in tidx], PERIODU[ti] == ((C[ti]*1000/l[ti])^(1-elasmu)-1)/(1-elasmu)-1)

# Period utility
@constraint(m, totperiodueq[ti in tidx], TOTPERIODU[ti] == PERIODU[ti] * l[ti] * rr[ti])

# Total utility
@constraint(m, utileq, UTILITY == tstep * scale1 * sum(TOTPERIODU[ti] for ti in tidx) + scale2)


# -------------------------------------------------------------------- 
# Other upper and lower bounds for stability

@constraint(m, alphalb[ti in tidx] , ALPHA[ti] >= 0.1)
@constraint(m, alphaub[ti in tidx] , ALPHA[ti] <= 100.0)
@constraint(m, miuub[ti in tidx] , MIU[ti] <= miuup[ti])
@constraint(m, sfix[ti in tidx[38:end]] , S[ti] == 0.28)

#####################################################################
# OBJECTIVE FUNCTION
#####################################################################
@objective(m, Max, UTILITY)

#####################################################################
# OPTIMIZATION
#####################################################################

set_optimizer_attribute(m, "print_level", 5)
optimize!(m)
status= termination_status(m)

The model of course doesn’t work and doesn’t even try to iterate, quickly returning a nice “EXIT: Invalid number in NLP function or derivative detected” (feel like an editor rejection when at least you didn’t waste time…).

In this relatively large model, how can I debug it? How do I find the culprit? Is there an JuMP/Optimizer API that translates the Jacobian entry that erroed back into the original constraint of my JuMP model? I have look at the documentation, but there isn’t even a section for this, but I am pretty confident you can get this information in GAMS…

2
This is the GAMS code (with merged included files) 
$ontext
DICE2023-beta-4-3-10 as of October 16, 2023
$offtext

$title October 12, 2023 (DICE2023-b-4-3-10.gms)

set        t  Time periods (5 years per period)                   /1*81/
PARAMETERS
** If optimal control
        ifopt    Indicator where optimized is 1 and base is 0     /1/
** Population and technology
        gama     Capital elasticity in production function        /.300    /
        pop1     Initial world population 2020 (millions)         /7752.9  /
        popadj   Growth rate to calibrate to 2050 pop projection  /0.145   /
        popasym  Asymptotic population (millions)                 /10825.  /
        dk       Depreciation rate on capital (per year)          /0.100    /
        q1       Initial world output 2020 (trill 2019 USD)       /135.7   /
        AL1       Initial level of total factor productivity       /5.84 /
        gA1      Initial growth rate for TFP per 5 years          /0.066   /
        delA     Decline rate of TFP per 5 years                  /0.0015  /        
** Emissions parameters and Non-CO2 GHG with sigma = emissions/output
        gsigma1   Initial growth of sigma (per year)                   / -0.015 /
        delgsig   Decline rate of gsigma per period                    /.96/
        asymgsig   Asympototic gsigma                                  /-.005/
        e1        Industrial emissions 2020 (GtCO2 per year)           / 37.56  /
        miu1      Emissions control rate historical 2020               / .05    /
        fosslim   Maximum cumulative extraction fossil fuels (GtC)     / 6000   /
        CumEmiss0 Cumulative emissions 2020 (GtC)                      / 633.5/
* Climate damage parameters
        a1        Damage intercept                                    /0      /
        a2base    Damage quadratic term rev 01-13-23                  /0.003467/
        a3        Damage exponent                                     /2.00   /
** Abatement cost
        expcost2  Exponent of control cost function                   / 2.6  /
        pback2050 Cost of backstop 2019$ per tCO2 2050                / 515.  /
        gback     Initial cost decline backstop cost per year         / -.012 /
        cprice1   Carbon price 2020 2019$ per tCO2                    / 6    /
        gcprice   Growth rate of base carbon price per year           /.025   /
** Limits on emissions controls
        limmiu2070    Emission control limit from 2070          /1.0/
        limmiu2120    Emission control limit from 2120          /1.1/
        limmiu2200    Emission control limit from 2220          /1.05/
        limmiu2300    Emission control limit from 2300          /1.0/      
        delmiumax     Emission control delta limit per period   / 0.12/
** Preferences, growth uncertainty, and timing
        betaclim  Climate beta                                      / 0.5  /
        elasmu    Elasticity of marginal utility of consumption     / 0.95 /
        prstp     Pure rate of social time preference               /.001/
        pi        Capital risk premium                              / .05  /
        rartp       Risk-adjusted rate of time preference
        k0        Initial capital stock calibrated (1012 2019 USD)  / 295  /
        siggc1    Annual standard deviation of consumption growth   / .01 /
        sig1      Carbon intensity 2020 kgCO2-output 2020
** Scaling so that MU(C(1)) = 1 and objective function = PV consumption
        tstep       Years per Period                               / 5  /
        SRF         Scaling factor discounting                     /1000000/
        scale1      Multiplicative scaling coefficient             / 0.00891061 /
        scale2      Additive scaling coefficient                   /-6275.91/ ;
** Program control variables
sets     tfirst(t), tlast(t), tearly(t), tlate(t);
PARAMETERS
        L(t)           Level of population and labor
        aL(t)          Level of total factor productivity
        sigma(t)       CO2-emissions output ratio
        sigmatot(t)    GHG-output ratio
        gA(t)          Growth rate of productivity from
        gL(t)          Growth rate of labor and population
        gcost1         Growth of cost factor
        gsig(t)        Change in sigma (rate of decarbonization)
        eland(t)       Emissions from deforestation (GtCO2 per year)
        cost1tot(T)    Abatement cost adjusted for backstop and sigma
        pbacktime(t)   Backstop price 2019$ per ton CO2
        optlrsav       Optimal long-run savings rate used for transversality
        scc(t)         Social cost of carbon
        cpricebase(t)  Carbon price in base case
        ppm(t)         Atmospheric concentrations parts per million
        atfrac2020(t)  Atmospheric share since 2020
        atfrac1765(t)  Atmospheric fraction of emissions since 1765
        abaterat(t)    Abatement cost per net output
        miuup(t)       Upper bound on miu
        gbacktime(t)   Decline rate of backstop price
** Precautionary dynamic parameters
        varpcc(t)         Variance of per capita consumption 
        rprecaut(t)       Precautionary rate of return
        RR(t)             STP with precautionary factor
        RR1(t)            STP factor without precautionary factor;
** Time preference for climate investments and precautionary effect
        rartp           = exp( prstp + betaclim*pi)-1;
        varpcc(t)       =  min(Siggc1**2*5*(t.val-1),Siggc1**2*5*47);       
        rprecaut(t)     = -0.5*varpcc(t)*elasmu**2;
        RR1(t)          = 1/((1+rartp)**(tstep*(t.val-1)));
        RR(t)           = RR1(t)*(1+rprecaut(t))**(-tstep*(t.val-1));
        L("1") = pop1; loop(t, L(t+1)=L(t);); loop(t, L(t+1)  = L(t)*(popasym/L(t))**popadj ;);
        gA(t)           = gA1*exp(-delA*5*((t.val-1))); aL("1") = AL1; loop(t, aL(t+1)=aL(t)/((1-gA(t))););
        optlrsav        =(dk + .004)/(dk + .004*elasmu + rartp)*gama;
        cpricebase(t)   = cprice1*(1+gcprice)**(5*(t.val-1));
        pbacktime(t)    = pback2050*exp(-5*(.01)*(t.val-7));
        pbacktime(t)$(t.val > 7) = pback2050*exp(-5*(.001)*(t.val-7));
        sig1            = e1/(q1*(1-miu1)); sigma("1") = sig1;
        gsig(t)         = min(gsigma1*delgsig **((t.val-1)),asymgsig);
        loop(t, sigma(t+1)  = sigma(t)*exp(5*gsig(t)););
** Emissions limits
        miuup('1')= .05; miuup('2')= .10; miuup(t)$(t.val > 2)  = (delmiumax*(t.val-1));
        miuup(t)$(t.val > 8)  = 0.85+.05*(t.val-8); miuup(t)$(t.val > 11) = limmiu2070;
        miuup(t)$(t.val > 20) = limmiu2120; miuup(t)$(t.val > 37) = limmiu2200; miuup(t)$(t.val > 57) = limmiu2300;       
** Include file for non-CO2 GHGs
* ******************************************************************************
* $include Include\Nonco2-b-4-3-1.gms

* nonco2 Parameters      
Parameters
        CO2E_GHGabateB(t)         Abateable non-CO2 GHG emissions base
        CO2E_GHGabateact(t)       Abateable non-CO2 GHG emissions base (actual)
        F_Misc(t)                 Non-abateable forcings (GHG and other)
        emissrat(t)               Ratio of CO2e to industrial emissions
        sigmatot(t)               Emissions output ratio for CO2e       
        FORC_CO2(t)               CO2 Forcings
;      
** Parameters for non-industrial emission
** Assumes abateable share of non-CO2 GHG is 65%
Parameters
        eland0         Carbon emissions from land 2015 (GtCO2 per year)  / 5.9    /
        deland         Decline rate of land emissions (per period)       / .1     /      
        F_Misc2020     Non-abatable forcings 2020                       /  -0.054    /
        F_Misc2100     Non-abatable forcings 2100                        / .265/       
        F_GHGabate2020 Forcings of abatable nonCO2 GHG                   / 0.518 /
        F_GHGabate2100 Forcings of abatable nonCO2 GHG                   / 0.957 /

        ECO2eGHGB2020  Emis of abatable nonCO2 GHG GtCO2e  2020             /  9.96/     
        ECO2eGHGB2100  Emis of abatable nonCO2 GHG GtCO2e  2100             /  15.5 /     
        emissrat2020   Ratio of CO2e to industrial CO2 2020                 / 1.40 /
        emissrat2100   Ratio of CO2e to industrial CO2 2020                 / 1.21 /       
        Fcoef1         Coefficient of nonco2 abateable emissions            /0.00955/
        Fcoef2         Coefficient of nonco2 abateable emissions            /.861/
        ;
** Parameters emissions and non-CO2 
        eland(t) = eland0*(1-deland)**(t.val-1); eland(t) = eland0*(1-deland)**(t.val-1);       
        CO2E_GHGabateB(t)=ECO2eGHGB2020+((ECO2eGHGB2100-ECO2eGHGB2020)/16)*(t.val-1)$(t.val le 16)+((ECO2eGHGB2100-ECO2eGHGB2020))$(t.val ge 17);
        F_Misc(t)=F_Misc2020 +((F_Misc2100-F_Misc2020)/16)*(t.val-1)$(t.val le 16)+((F_Misc2100-F_Misc2020))$(t.val ge 17);
        emissrat(t) = emissrat2020 +((emissrat2100-emissrat2020)/16)*(t.val-1)$(t.val le 16)+((emissrat2100-emissrat2020))$(t.val ge 17); 
        sigmatot(t) = sigma(t)*emissrat(t);
        cost1tot(t) = pbacktime(T)*sigmatot(T)/expcost2/1000;
VARIABLES
        ECO2(t)         Total CO2 emissions (GtCO2 per year)
        ECO2E(t)        Total CO2e emissions including abateable nonCO2 GHG (GtCO2 per year)
        EIND(t)         Industrial CO2 emissions (GtCO2 per yr)
        F_GHGabate      Forcings abateable nonCO2 GHG     
;
Equations
        ECO2eq(t)         CO2 Emissions equation
        ECO2Eeq(t)        CO2E Emissions equation
        EINDeq(t)        Industrial CO2 equation
        F_GHGabateEQ(t)
;  
* ***************************************************************
* Program control definitions
        tfirst(t) = yes$(t.val eq 1);
        tlast(t)  = yes$(t.val eq card(t));
VARIABLES
        MIU(t)          Emission control rate GHGs
        C(t)            Consumption (trillions 2019 US dollars per year)
        K(t)            Capital stock (trillions 2019 US dollars)
        CPC(t)          Per capita consumption (thousands 2019 USD per year)
        I(t)            Investment (trillions 2019 USD per year)
        S(t)            Gross savings rate as fraction of gross world product
        Y(t)            Gross world product net of abatement and damages (trillions 2019 USD per year)
        YGROSS(t)       Gross world product GROSS of abatement and damages (trillions 2019 USD per year)
        YNET(t)         Output net of damages equation (trillions 2019 USD per year)
        DAMAGES(t)      Damages (trillions 2019 USD per year)
        DAMFRAC(t)      Damages as fraction of gross output
        ABATECOST(t)    Cost of emissions reductions  (trillions 2019 USD per year)
        MCABATE(t)      Marginal cost of abatement (2019$ per ton CO2)
        CCATOT(t)       Total carbon emissions (GtC)
        PERIODU(t)      One period utility function
        CPRICE(t)       Carbon price (2019$ per ton of CO2)
        TOTPERIODU(t)   Period utility
        UTILITY         Welfare function
        RFACTLONG(t)
        RSHORT(t)       Real interest rate with precautionary(per annum year on year)
        RLONG(t)        Real interest rate from year 0 to T
;
NONNEGATIVE VARIABLES  MIU, TATM, MAT, MU, ML, Y, YNET, YGROSS, C, K, I, RFACTLONG;
EQUATIONS
*Emissions and Damages
        CCATOTEQ(t)      Cumulative total carbon emissions
        DAMFRACEQ(t)     Equation for damage fraction
        DAMEQ(t)         Damage equation
        ABATEEQ(t)       Cost of emissions reductions equation
        MCABATEEQ(t)     Equation for MC abatement
        CARBPRICEEQ(t)   Carbon price equation from abatement
*Economic variables
        YGROSSEQ(t)      Output gross equation
        YNETEQ(t)        Output net of damages equation
        YY(t)            Output net equation
        CC(t)            Consumption equation
        CPCE(t)          Per capita consumption definition
        SEQ(t)           Savings rate equation
        KK(t)            Capital balance equation
        RSHORTEQ(t)      Short-run interest rate equation
        RLONGeq(t)       Long-run interest rate equation
        RFACTLONGeq(t)   Long interest factor
* Utility
        TOTPERIODUEQ(t)  Period utility
        PERIODUEQ(t)     Instantaneous utility function equation
        UTILEQ           Objective function      ;

** Include file for DFAIR model and climate equations
* *******************************************************************
* $include Include\FAIR-beta-4-3-1.gms

** Equals old FAIR with recalibrated parameters for revised F2xco2 and Millar model.
** Deletes nonnegative reservoirs. See explanation below

sets     tfirst(t), tlast(t);

PARAMETERS
         yr0     Calendar year that corresponds to model year zero         /2020/
        emshare0 Carbon emissions share into Reservoir 0   /0.2173/
        emshare1 Carbon emissions share into Reservoir 1    /0.224/
        emshare2 Carbon emissions share into Reservoir 2    /0.2824/
        emshare3 Carbon emissions share into Reservoir 3    /0.2763/
        tau0    Decay time constant for R0  (year)                            /1000000/
        tau1    Decay time constant for R1  (year)                            /394.4/
        tau2    Decay time constant for R2  (year)       /36.53/
        tau3    Decay time constant for R3  (year) /4.304/

        teq1    Thermal equilibration parameter for box 1 (m^2 per KW)         /0.324/
        teq2    Thermal equilibration parameter for box 2 (m^2 per KW)        /0.44/
        d1      Thermal response timescale for deep ocean (year)               /236/
        d2      Thermal response timescale for upper ocean (year)              /4.07/
 
        irf0    Pre-industrial IRF100 (year)                                        /32.4/
        irC      Increase in IRF100 with cumulative carbon uptake (years per GtC)    /0.019/
        irT      Increase in IRF100 with warming (years per degree K)                /4.165/
        fco22x   Forcings of equilibrium CO2 doubling (Wm-2)                        /3.93/       

** INITIAL CONDITIONS TO BE CALIBRATED TO HISTORY
** CALIBRATION
       mat0   Initial concentration in atmosphere in 2020 (GtC)       /886.5128014/

       res00  Initial concentration in Reservoir 0 in 2020 (GtC)      /150.093 /
       res10  Initial concentration in Reservior 1 in 2020 (GtC)      /102.698 /
       res20  Initial concentration in Reservoir 2 in 2020 (GtC)      /39.534  /
       res30  Initial concentration in Reservoir 3 in 2020 (GtC)      / 6.1865 /

       mateq      Equilibrium concentration atmosphere  (GtC)            /588   /
       tbox10    Initial temperature box 1 change in 2020 (C from 1765)  /0.1477  /
       tbox20    Initial temperature box 2 change in 2020 (C from 1765)  /1.099454/
       tatm0     Initial atmospheric temperature change in 2020          /1.24715 /     
        
;
VARIABLES
*Note: Stock variables correspond to levels at the END of the period
        FORC(t)        Increase in radiative forcing (watts per m2 from 1765)
        TATM(t)        Increase temperature of atmosphere (degrees C from 1765)     
        TBOX1(t)       Increase temperature of box 1 (degrees C from 1765)
        TBOX2(t)       Increase temperature of box 2 (degrees C from 1765)
        RES0(t)        Carbon concentration in Reservoir 0 (GtC from 1765)
        RES1(t)        Carbon concentration in Reservoir 1 (GtC from 1765)
        RES2(t)        Carbon concentration in Reservoir 2 (GtC from 1765)
        RES3(t)        Carbon concentration in Reservoir 3 (GtC from 1765)
        MAT(t)         Carbon concentration increase in atmosphere (GtC from 1765)
        CACC(t)        Accumulated carbon in ocean and other sinks (GtC)
        IRFt(t)        IRF100 at time t
        alpha(t)       Carbon decay time scaling factor
        SumAlpha      Placeholder variable for objective function;

**** IMPORTANT PROGRAMMING NOTE. Earlier implementations has reservoirs as non-negative.
**** However, these are not physical but mathematical solutions.
**** So, they need to be unconstrained so that can have negative emissions.
NONNEGATIVE VARIABLES   TATM, MAT,  IRFt, alpha

EQUATIONS
        FORCE(t)        Radiative forcing equation
        RES0LOM(t)      Reservoir 0 law of motion
        RES1LOM(t)      Reservoir 1 law of motion
        RES2LOM(t)      Reservoir 2 law of motion
        RES3LOM(t)      Reservoir 3 law of motion
        MMAT(t)         Atmospheric concentration equation
        Cacceq(t)       Accumulated carbon in sinks equation
        TATMEQ(t)       Temperature-climate equation for atmosphere
        TBOX1EQ(t)      Temperature box 1 law of motion
        TBOX2EQ(t)      Temperature box 2 law of motion
        IRFeqLHS(t)     Left-hand side of IRF100 equation
        IRFeqRHS(t)     Right-hand side of IRF100 equation
;
** Equations of the model
    res0lom(t+1)..   RES0(t+1) =E=  (emshare0*tau0*alpha(t+1)*(Eco2(t+1)/3.667))*(1-exp(-tstep/(tau0*alpha(t+1))))+Res0(t)*exp(-tstep/(tau0*alpha(t+1)));
    res1lom(t+1)..   RES1(t+1) =E=  (emshare1*tau1*alpha(t+1)*(Eco2(t+1)/3.667))*(1-exp(-tstep/(tau1*alpha(t+1))))+Res1(t)*exp(-tstep/(tau1*alpha(t+1)));
    res2lom(t+1)..   RES2(t+1) =E=  (emshare2*tau2*alpha(t+1)*(Eco2(t+1)/3.667))*(1-exp(-tstep/(tau2*alpha(t+1))))+Res2(t)*exp(-tstep/(tau2*alpha(t+1)));
    res3lom(t+1)..   RES3(t+1) =E=  (emshare3*tau3*alpha(t+1)*(Eco2(t+1)/3.667))*(1-exp(-tstep/(tau3*alpha(t+1))))+Res3(t)*exp(-tstep/(tau3*alpha(t+1)));
    mmat(t+1)..      MAT(t+1)  =E=   mateq+Res0(t+1)+Res1(t+1)+Res2(t+1)+Res3(t+1);
    cacceq(t)..      Cacc(t)   =E=  (CCATOT(t)-(MAT(t)-mateq));
    force(t)..       FORC(t)    =E=  fco22x*((log((MAT(t)/mateq))/log(2))) + F_Misc(t)+F_GHGabate(t);
    
    tbox1eq(t+1)..   Tbox1(t+1) =E=  Tbox1(t)*exp(-tstep/d1)+teq1*Forc(t+1)*(1-exp(-tstep/d1));
    tbox2eq(t+1)..   Tbox2(t+1) =E=  Tbox2(t)*exp(-tstep/d2)+teq2*Forc(t+1)*(1-exp(-tstep/d2));
    tatmeq(t+1)..    TATM(t+1)  =E=   Tbox1(t+1)+Tbox2(t+1);
    irfeqlhs(t)..    IRFt(t)   =E=  ((alpha(t)*emshare0*tau0*(1-exp(-100/(alpha(t)*tau0))))+(alpha(t)*emshare1*tau1*(1-exp(-100/(alpha(t)*tau1))))+(alpha(t)*emshare2*tau2*(1-exp(-100/(alpha(t)*tau2))))+(alpha(t)*emshare3*tau3*(1-exp(-100/(alpha(t)*tau3)))));
    irfeqrhs(t)..    IRFt(t)   =E=  irf0+irC*Cacc(t)+irT*TATM(t);

**  Upper and lower bounds for stability
MAT.LO(t)       = 10;
TATM.UP(t)      = 20;
TATM.lo(t)      = .5;
alpha.up(t) = 100;
alpha.lo(t) = 0.1;

* Initial conditions
MAT.FX(tfirst)    = mat0;
TATM.FX(tfirst)   = tatm0;
Res0.fx(tfirst) = Res00;
Res1.fx(tfirst) = Res10;
Res2.fx(tfirst) = Res20;
Res3.fx(tfirst) = Res30;
Tbox1.fx(tfirst) = Tbox10;
Tbox2.fx(tfirst) = Tbox20;

** Solution options
option iterlim = 99900;
option reslim = 99999;
option solprint = on;
option limrow = 0;
option limcol = 0;
 

* *******************************************************************

**** Equations of the model
**Emissions and Damages
 eco2eq(t)..          ECO2(t)        =E= (sigma(t)*YGROSS(t) + eland(t))*(1-(MIU(t)));
 eindeq(t)..          EIND(t)        =E= (sigma(t)*YGROSS(t))*(1-(MIU(t)));
 eco2Eeq(t)..         ECO2E(t)       =E= (sigma(t)*YGROSS(t) + eland(t) + CO2E_GHGabateB(t))*(1-(MIU(t)));
 F_GHGabateEQ(t+1)..  F_GHGabate(t+1) =E= Fcoef2*F_GHGabate(t)+ Fcoef1*CO2E_GHGabateB(t)*(1-(MIU(t)));
 ccatoteq(t+1)..      CCATOT(t+1)    =E= CCATOT(t) +  ECO2(T)*(5/3.666) ;
 damfraceq(t) ..      DAMFRAC(t)     =E= (a1*(TATM(t)))+(a2base*(TATM(t))**a3) ;
 dameq(t)..           DAMAGES(t)     =E= YGROSS(t) * DAMFRAC(t);
 abateeq(T)..         ABATECOST(T)   =E= YGROSS(T) * COST1TOT(T) * (MIU(T)**EXPCOST2);
 mcabateeq(t)..       MCABATE(t)     =E= pbacktime(t) * MIU(t)**(expcost2-1);
 carbpriceeq(t)..     CPRICE(t)      =E= pbacktime(t) * (MIU(t))**(expcost2-1);
**Economic variables
 ygrosseq(t)..        YGROSS(t)      =E= (AL(t)*(L(t)/1000)**(1-gama))*(K(t)**gama);
 yneteq(t)..          YNET(t)        =E= YGROSS(t)*(1-damfrac(t));
 yy(t)..              Y(t)           =E= YNET(t) - ABATECOST(t);
 cc(t)..              C(t)           =E= Y(t) - I(t);
 cpce(t)..            CPC(t)         =E= 1000 * C(t) / L(t);
 seq(t)..             I(t)           =E= S(t) * Y(t);
 kk(t+1)..            K(t+1)         =L= (1-dk)**tstep * K(t) + tstep * I(t);
 RFACTLONGeq(t+1)..   RFACTLONG(t+1) =E= SRF*(cpc(t+1)/cpc('1'))**(-elasmu)*rr(t+1);
 RLONGeq(t+1)..       RLONG(t+1)     =E= -log(RFACTLONG(t+1)/SRF)/(5*t.val);
 RSHORTeq(t+1)..      RSHORT(t+1)    =E= -log(RFACTLONG(t+1)/Rfactlong(t))/5;
** Welfare functions 
 periodueq(t)..       PERIODU(t)     =E= ((C(T)*1000/L(T))**(1-elasmu)-1)/(1-elasmu)-1;
 totperiodueq(t)..    TOTPERIODU(t)  =E= PERIODU(t) * L(t) * RR(t);
 utileq..             UTILITY        =E= tstep * scale1 * sum(t,  TOTPERIODU(t)) + scale2;

* Various control rate limits, initial and terminal conditions
miu.up(t)       = miuup(t);
K.LO(t)         = 1;
C.LO(t)         = 2;
CPC.LO(t)       = .01;
RFACTLONG.lo(t) =.0001;
*set lag10(t) ;
*lag10(t)                =  yes$(t.val gt card(t)-10);
*S.FX(lag10(t))          = optlrsav;
s.fx(t)$(t.val > 37)    =.28;
ccatot.fx(tfirst)       = CumEmiss0;
K.FX(tfirst)            = k0;
F_GHGabate.fx(tfirst)   = F_GHGabate2020;
RFACTLONG.fx(tfirst)    = 1000000; 

** Solution options
option iterlim = 99900; option reslim = 99999; option solprint = on; option limrow = 0; option limcol = 0;

** Initial solution
model  CO2 /all/;
ifopt=1;
solve CO2 maximizing UTILITY using nlp ;
 
**** STATEMENTS FOR DEFINITIONS AND PUT STATEMENTS FOR MAJOR SCENARIOS
$include Include\Putlong-4-3-10.gms

*DISPLAY FOR MAJOR VARIABLES
option decimals = 6;
display gpccons, varpcc,RLONG.l, rprecaut,RSHORT.l,cpc.l,rr,miu.l,scc, mat.l,tatm.l,eco2eq.m,miu.l;
display ifopt, elasmu, rartp, pi, k0, betaclim, prstp,k0,siggc1,delA,a2base, utility.l ;
3

Solved the particular issue removing the 3 equations on long and short term interest rate, that can go to post-processing, as they do not influence the optimization.

Still I have the general question… how I could have obtained, instead of manually trying the various equations, an informative message telling me where the problem was…

4

A tool to debug starting points is one of my goals for Tools to test and debug JuMP models · Issue #3664 · jump-dev/JuMP.jl · GitHub

You have a bunch of variables that are defined like:

    MAT[tidx] >= 0.0       # Carbon concentration increase in atmosphere (GtC from 1765)

and then you have a constraint:

@constraint(m, matlb[ti in tidx], MAT[ti] >= 10.0)

This constraint does not set the lower bound of MAT to 10, it adds a new linear constraint 1.0 * MAT[ti] >= 10. The lower bound of MAT is 0.0.

By default, JuMP sets that start value of a variable to 0 projected onto the variable bounds, so we start with MAT[ti] = 0.0.

However, MAT appears in a log of this constraint:

@constraint(m, force[ti in tidx], FORC[ti] == fco22x*((log((MAT[ti]/mateq))/log(2))) + f_misc[ti]+F_GHGABATE[ti] )

which is not defined when MAT = 0.0, hence the error.

Here’s how I would write your model:

model = Model(Ipopt.Optimizer)
@variables(model, begin
    ECO2[tidx]             # Total CO2 emissions (GtCO2 per year)
    ECO2E[tidx]            # Total CO2e emissions including abateable nonCO2 GHG (GtCO2 per year)
    EIND[tidx]             # Industrial CO2 emissions (GtCO2 per yr)
    # TODO: check this, as in GAMS it is declared without time index, but then it is used with time index!
    F_GHGABATE[tidx]       # Forcings abateable nonCO2 GHG
    0 <= MIU[ti in tidx] <= miuup[ti]       # Emission control rate GHGs
    C[tidx] >= 2.0         # Consumption (trillions 2019 US dollars per year)
    K[tidx] >= 0.0         # Capital stock (trillions 2019 US dollars)
    CPC[tidx] >= 0.01      # Per capita consumption (thousands 2019 USD per year)
    I[tidx] >= 0.0         # Investment (trillions 2019 USD per year)
    S[tidx]                # Gross savings rate as fraction of gross world product (control)
    Y[tidx] >= 0.0         # Gross world product net of abatement and damages (trillions 2019 USD per year)
    YGROSS[tidx] >= 0.0    # Gross world product GROSS of abatement and damages (trillions 2019 USD per year)
    YNET[tidx] >= 0.0      # Output net of damages equation (trillions 2019 USD per year)
    DAMAGES[tidx]          # Damages (trillions 2019 USD per year)
    DAMFRAC[tidx]          # Damages as fraction of gross output
    ABATECOST[tidx]        # Cost of emissions reductions  (trillions 2019 USD per year)
    MCABATE[tidx]          # Marginal cost of abatement (2019$ per ton CO2)
    CCATOT[tidx]           # Total carbon emissions (GtC)
    PERIODU[tidx]          # One period utility function
    CPRICE[tidx]           # Carbon price (2019$ per ton of CO2)
    TOTPERIODU[tidx]       # Period utility
    UTILITY                # Welfare function
    RFACTLONG[tidx] >= 0.0 #
    RSHORT[tidx[2:end]]    # Real interest rate with precautionary(per annum year on year)
    RLONG[tidx[2:end]]     # Real interest rate from year 0 to T
    # --------------------------------------------------------------------
    # Climate variables
    # Note: Stock variables correspond to levels at the END of the period
    FORC[tidx]             # Increase in radiative forcing (watts per m2 from 1765)
    0.5 <= TATM[tidx] <= 20.0      # Increase temperature of atmosphere (degrees C from 1765)
    TBOX1[tidx]            # Increase temperature of box 1 (degrees C from 1765)
    TBOX2[tidx]            # Increase temperature of box 2 (degrees C from 1765)
    RES0[tidx]             # Carbon concentration in Reservoir 0 (GtC from 1765)
    RES1[tidx]             # Carbon concentration in Reservoir 1 (GtC from 1765)
    RES2[tidx]             # Carbon concentration in Reservoir 2 (GtC from 1765)
    RES3[tidx]             # Carbon concentration in Reservoir 3 (GtC from 1765)
    MAT[tidx] >= 10.0       # Carbon concentration increase in atmosphere (GtC from 1765)
    CACC[tidx]             # Accumulated carbon in ocean and other sinks (GtC)
    IRFT[tidx] >= 0.0      # IRF100 at time t
    0.1 <= ALPHA[tidx] <= 100.0     # Carbon decay time scaling factor (control)
end)
for ti in tidx[38:end]
    fix(S[ti], 0.28)
end
@constraints(model, begin
    # Reservoir 0 law of motion
    [ti in tidx], RES0[ti] ==  ((ti == 1) ? res00 : (emshare0*tau0*ALPHA[ti]*(ECO2[ti]/3.667))*(1-exp(-tstep/(tau0*ALPHA[ti])))+RES0[ti-1]*exp(-tstep/(tau0*ALPHA[ti])))
    # Reservoir 1 law of motion
    [ti in tidx], RES1[ti] == ((ti == 1) ? res10 :  (emshare1*tau1*ALPHA[ti]*(ECO2[ti]/3.667))*(1-exp(-tstep/(tau1*ALPHA[ti])))+RES1[ti-1]*exp(-tstep/(tau1*ALPHA[ti])))
    # Reservoir 2 law of motion
    [ti in tidx], RES2[ti] ==  ((ti == 1) ? res20 : (emshare2*tau2*ALPHA[ti]*(ECO2[ti]/3.667))*(1-exp(-tstep/(tau2*ALPHA[ti])))+RES2[ti-1]*exp(-tstep/(tau2*ALPHA[ti])))
    # Reservoir 3 law of motion
    [ti in tidx], RES3[ti] ==  ((ti == 1) ? res30 : (emshare3*tau3*ALPHA[ti]*(ECO2[ti]/3.667))*(1-exp(-tstep/(tau3*ALPHA[ti])))+RES3[ti-1]*exp(-tstep/(tau3*ALPHA[ti])))
    # Atmospheric concentration equation
    [ti in tidx], MAT[ti] == ((ti == 1) ? mat0 : mateq + RES0[ti] + RES1[ti] + RES2[ti] + RES3[ti])
    # Accumulated carbon in sinks equation
    [ti in tidx], CACC[ti] == CCATOT[ti]-(MAT[ti]-mateq)
    # Radiative forcing equation
    [ti in tidx], FORC[ti] == fco22x*((log((MAT[ti]/mateq))/log(2))) + f_misc[ti]+F_GHGABATE[ti]
    # Temperature box 1 law of motion
    [ti in tidx], TBOX1[ti] ==  ((ti == 1) ? tbox10 : TBOX1[ti-1]*exp(-tstep/d1)+teq1*FORC[ti]*(1-exp(-tstep/d1)))
    # Temperature box 2 law of motion
    [ti in tidx[2:end]], TBOX2[ti] ==  ((ti == 1) ? tbox20 : TBOX2[ti-1]*exp(-tstep/d2)+teq2*FORC[ti]*(1-exp(-tstep/d2)))
    # Temperature-climate equation for atmosphere
    [ti in tidx], TATM[ti] ==  ((ti == 1) ? tatm0 : TBOX1[ti]+TBOX2[ti])
    # Left-hand side of IRF100 equation
    [ti in tidx], IRFT[ti] == (ALPHA[ti]*emshare0*tau0*(1-exp(-100/(ALPHA[ti]*tau0))))+(ALPHA[ti]*emshare1*tau1*(1-exp(-100/(ALPHA[ti]*tau1))))+(ALPHA[ti]*emshare2*tau2*(1-exp(-100/(ALPHA[ti]*tau2))))+(ALPHA[ti]*emshare3*tau3*(1-exp(-100/(ALPHA[ti]*tau3))))
    # Right-hand side of IRF100 equation
    [ti in tidx], IRFT[ti] == irf0+irc*CACC[ti]+irt*TATM[ti]
    # --------------------------------------------------------------------
    # Emissions and Damages
    # CO2 Emissions equation
    [ti in tidx], ECO2[ti] == (sigma[ti] * YGROSS[ti] + eland[ti]) * (1 - MIU[ti])
    # Industrial CO2 equation
    [ti in tidx], EIND[ti] == (sigma[ti] * YGROSS[ti]) * (1 - MIU[ti])
    # CO2E Emissions equation
    [ti in tidx], ECO2E[ti] == (sigma[ti]*YGROSS[ti] + eland[ti] + co2e_ghgabateb[ti])*(1-MIU[ti])
    # Forcings abateable nonCO2 GHG equation
    [ti in tidx], F_GHGABATE[ti] == ((ti == 1) ? f_ghgabate2020 : fcoef2*F_GHGABATE[ti-1]+ fcoef1*co2e_ghgabateb[ti-1]*(1-MIU[ti-1]))
    # --------------------------------------------------------------------
    # Emissions and damage
    # Cumulative total carbon emissions
    [ti in tidx], CCATOT[ti] == ((ti==1) ? cumemiss0 : CCATOT[ti-1] +  ECO2[ti-1]*(tstep/3.666))
    # Equation for damage fraction
    [ti in tidx], DAMFRAC[ti] == (a1*TATM[ti])+(a2base*TATM[ti]^a3)
    # Damage equation
    [ti in tidx], DAMAGES[ti] == YGROSS[ti] * DAMFRAC[ti]
    # Cost of emissions reductions equation
    [ti in tidx], ABATECOST[ti] == YGROSS[ti] * cost1tot[ti] * (MIU[ti]^expcost2)
    # Equation for MC abatement
    [ti in tidx], MCABATE[ti] == pbacktime[ti] * MIU[ti]^(expcost2-1)
    # Carbon price equation from abatement
    [ti in tidx], CPRICE[ti] == pbacktime[ti] * MIU[ti]^(expcost2-1)
    # --------------------------------------------------------------------
    # Economic variables
    # Output gross equation
    [ti in tidx], YGROSS[ti] == (al[ti]*(l[ti]/1000)^(1-gama))*(K[ti]^gama)
    # Output net of damage equation
    [ti in tidx], YNET[ti] == YGROSS[ti]*(1-DAMFRAC[ti])
    # Output net equation
    [ti in tidx], Y[ti] == YNET[ti] - ABATECOST[ti]
    # Consumption equation
    [ti in tidx], C[ti] == Y[ti] - I[ti]
    # Per capita consumption definition
    [ti in tidx], CPC[ti] == 1000 * C[ti] / l[ti]
    # Savings rate equation
    [ti in tidx], I[ti] == S[ti] * Y[ti]
    # Capital balance equation
    K[1] == k0
    [ti in tidx[2:end]],  K[ti] <= (1-dk)^tstep * K[ti-1] + tstep * I[ti-1]
    [ti in tidx],  K[ti] >= 1
    # Long interest factor
    [ti in tidx],  RFACTLONG[ti] >= 0.0001
    [ti in tidx],  RFACTLONG[ti] == ((ti ==1) ? 1000000 : srf *(CPC[ti]/CPC[1])^(-elasmu)*rr[ti])
    # Long-run interest rate equation
    [ti in tidx[2:end]], RLONG[ti] ==  -log(RFACTLONG[ti]/srf)/(tstep*ti)
    # Short-run interest rate equation
    [ti in tidx[2:end]], RSHORT[ti] == -log(RFACTLONG[ti]/RFACTLONG[ti-1])/tstep
    # --------------------------------------------------------------------
    # Utility
    # Instantaneous utility function equation
    [ti in tidx], PERIODU[ti] == ((C[ti]*1000/l[ti])^(1-elasmu)-1)/(1-elasmu)-1
    # Period utility
    [ti in tidx], TOTPERIODU[ti] == PERIODU[ti] * l[ti] * rr[ti]
    # Total utility
    UTILITY == tstep * scale1 * sum(TOTPERIODU[ti] for ti in tidx) + scale2
end)
@objective(model, Max, UTILITY)
set_optimizer_attribute(model, "print_level", 5)
optimize!(model)
is_solved_and_feasible(model)

It converges to a feasible point:

julia> solution_summary(model)
* Solver : Ipopt

* Status
  Result count       : 1
  Termination status : LOCALLY_SOLVED
  Message from the solver:
  "Solve_Succeeded"

* Candidate solution (result #1)
  Primal status      : FEASIBLE_POINT
  Dual status        : FEASIBLE_POINT
  Objective value    : 6.66490e+03
  Dual objective value : -6.21803e+03

* Work counters
  Solve time (sec)   : 2.55619e+00
  Barrier iterations : 268

I would also invite you to consider how you can use Julia to make this less like GAMS. For example, instead of:

emshare0 = 0.2173    # Carbon emissions share into Reservoir 0
emshare1 = 0.224     # Carbon emissions share into Reservoir 1
emshare2 = 0.2824    # Carbon emissions share into Reservoir 2
emshare3 = 0.2763    # Carbon emissions share into Reservoir 3

you might do:

emshare = [0.2173, 0.224, 0.2824, 0.2763]

and work it into:

res0  = [150.093, 102.698, 39.534, 6.1865]
emshare = [0.2173, 0.224, 0.2824, 0.2763]
tau = [1000000.0, 394.4, 36.53, 4.304]
@variable(model, RES[r in 1:4, ti in tidx])
@constraint(
    model,
    [r in 1:4, ti in tidx],
    RES[r, ti] == if ti == 1
        res0[r]
    else
        emshare[r] * tau[r] * ALPHA[ti] * ECO2[ti] / 3.667 * (1 - exp(-tstep / (tau[r] * ALPHA[ti]))) +
        RES[r, ti-1] * exp(-tstep / (tau[r] * ALPHA[ti]))
    end,
)
1 Like
5

Thank you very much, I will investigate the issue you reported concerning the bonds. The “similarity” to GAMS is intended in this specific case.

Concerning the general question however what you are proposing/implementing seems different.
You are introducing some analysis that find and report common modelling problems. What I am saying instead is something simpler, a way to report/map in the user model space what the solver reported as an error in its space. Is this possible ?

6

Ipopt doesn’t tell us where the issue is. It just says it found an issue.

I assume GAMS does what JuMP would like to do, which is for us to evaluate the various callbacks and look for NaN or Inf. This requires a bit of effort though. It’s non-trivial at the moment.

Three of the points in Tools to test and debug JuMP models · Issue #3664 · jump-dev/JuMP.jl · GitHub are relevant:

Bounds given as constraints

Identify constraints of the form @constraint(model, l <= x <= u). These should instead be specified as @variable(model, l <= x <= u).

Motivation: The former adds a new linear constraint row to the problem. The latter sets a variable’s bounds. Solvers typically have specialised support for variable bounds. If the solver doesn’t have specialised support for variable bounds, JuMP will reformulate into the constraint version.

Prior art: Most solvers have a presolve that does this, but it’s a sign of user-error that could be improved.

Starting point analysis

Starting points which violate domain restrictions like log(x) when x = 0 to start are common. Therefore, we should analyse the feasibility of a given starting point, and potentially that of the first- and second-derivatives as well.

Domain analysis

Nonlinear programs often include terms like log(x) where x >= 0. This is a common cause of domain violations because the solver may try x = 0. It’s usually better to add an explicit lower bound on x of x >= 0.01.

3 Likes