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

**URL:** https://discourse.julialang.org/t/dice-gams-jump-how-to-debug-jump-model-errors-find-the-culpit/125578
**Category:** Optimization (Mathematical)
**Tags:** question, jump, debug, climate, dice
**Created:** [February 5, 2025, 12:38pm UTC](https://discourse.julialang.org/t/dice-gams-jump-how-to-debug-jump-model-errors-find-the-culpit/125578 "2025-02-05T12:38:30Z")
**Posts on this page:** 7
**Page:** 1

<div class="post-metadata">

### Author: ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [February 5, 2025, 12:38pm UTC](https://discourse.julialang.org/t/dice-gams-jump-how-to-debug-jump-model-errors-find-the-culpit/125578/1 "2025-02-05T12:38:31Z")

</div>

Hello, this post il dual 🙂  
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**
>
> ```julia
> # 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…

---

<div class="post-metadata">

### Author: ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [February 5, 2025, 12:38pm UTC](https://discourse.julialang.org/t/dice-gams-jump-how-to-debug-jump-model-errors-find-the-culpit/125578/2 "2025-02-05T12:38:54Z")

</div>

> **This is the GAMS code (with merged included files)**
>
> ```gams
> $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 ;
> 
> ```

---

<div class="post-metadata">

### Author: ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [February 5, 2025, 1:58pm UTC](https://discourse.julialang.org/t/dice-gams-jump-how-to-debug-jump-model-errors-find-the-culpit/125578/3 "2025-02-05T13:58:18Z")

</div>

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…

---

<div class="post-metadata">

### Author: ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [February 5, 2025, 8:05pm UTC](https://discourse.julialang.org/t/dice-gams-jump-how-to-debug-jump-model-errors-find-the-culpit/125578/4 "2025-02-05T20:05:37Z")

</div>

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](https://github.com/jump-dev/JuMP.jl/issues/3664)

You have a bunch of variables that are defined like:

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

```

and then you have a constraint:

```julia
@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:

```julia
@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:

```Julia
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
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:

```julia
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:

```julia
emshare = [0.2173, 0.224, 0.2824, 0.2763]

```

and work it into:

```julia
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,
)

```

---

<div class="post-metadata">

### Author: ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [February 5, 2025, 9:40pm UTC](https://discourse.julialang.org/t/dice-gams-jump-how-to-debug-jump-model-errors-find-the-culpit/125578/5 "2025-02-05T21:40:26Z")

</div>

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 ?

---

<div class="post-metadata">

### Author: ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [February 5, 2025, 9:50pm UTC](https://discourse.julialang.org/t/dice-gams-jump-how-to-debug-jump-model-errors-find-the-culpit/125578/6 "2025-02-05T21:50:47Z")

</div>

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](https://github.com/jump-dev/JuMP.jl/issues/3664) 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`.

---

<div class="post-metadata">

### Author: ![sylvaticus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sylvaticus/32/203883_2.png) [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [February 15, 2025, 4:15pm UTC](https://discourse.julialang.org/t/dice-gams-jump-how-to-debug-jump-model-errors-find-the-culpit/125578/7 "2025-02-15T16:15:00Z")

</div>

This page from the GAMS documentation could be interesting too for debugging a problematic model, at least for the _kind_ of issue that could arise…

[https://www.gams.com/latest/docs/UG\_ExecErrPerformance.html](https://www.gams.com/latest/docs/UG_ExecErrPerformance.html)

I am impressed with the sentence: “_solvers report the type of arithmetic problem encountered and the problematic equation, but do not identify the particular offending variable or the labels in the index of an equation that cause the error._”