message_ix icon indicating copy to clipboard operation
message_ix copied to clipboard
verified publisher icon iiasa

Investigate PRICE_COMMODITY spikes with non-uniform duration_period

Open Wegatriespython opened this issue 4 months ago • 10 comments

What happened?

Image

For the baseline scenario under model SSP_SSP2_v6.1, price commodity for electricity shows a spike at 2060.

This has important downstream implications as any technology portfolio with electricity inputs in some small neighbourhood will see a basis flip to the cheaper option under the price spike barring other constraints. This happens for instance in the water model, where groundwater extraction is massively reduced in favour of a pre-dominantly surface water strategy which is problematic.

Now from a purely economic perspective such a price shock should be arbitraged away by the central planner under perfect foresight mode. Since the model doesn't do this, the only implication is that its cheaper for the model to have the price shock than to arbitrage it away. This however would invite the external validity Lucas critique as real world agents would be able to exploit an arbitrage the model cannot, breaking the policy invariance of the model.

The fact that this happens at the shift in the model timestep from 5-10 year periods is also suspicious.

  • This means that either there is some misspecification in how the model is discounting Present Value at the period length change at 2060
  • Or there are some discrete constraints creating a lumpy choice that forces this behaviour.

Code Sample

What did you expect to happen?

A strict 0 arbitrage price trajectory.

Versions

ixmp:              3.11.2.dev44+g1730afd24
              git: 1730afd (HEAD -> main, origin/main, origin/HEAD) Merge pull request #587 from iiasa/fix/gamsmodel-init
message_ix:        3.11.2.dev67+g00bec749
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
message_ix_models: (not installed)
message_data:      (not installed)

click:             8.2.1
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
dask:              2025.5.1
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
genno:             1.28.2
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
graphviz:          0.21
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
ixmp4:             0.10.0
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
jpype:             1.5.2
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
     Java VM path: /usr/lib/jvm/java-21-openjdk-amd64/lib/server/libjvm.so
openpyxl:          3.1.5
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
pandas:            2.2.3
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
pint:              0.24.4
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
xarray:            2025.4.0
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
yaml:              6.0.2
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock

iam_units:         2023.9.12
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
jupyter:           5.8.1
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
matplotlib:        3.10.3
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
plotnine:          0.15.0
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock
pyam:              3.0.0
              git: 00bec749 (HEAD -> main, origin/main, origin/HEAD) Merge pull request #451 from iiasa/material_stock

GAMS:              46.5.0
       system dir: /home/raghunathan/opt/gams/gams46.5_linux_x64_64_sfx

python:            3.12.3 (main, Jun 18 2025, 17:59:45) [GCC 13.3.0]
python-bits:       64
OS:                Linux
OS-release:        6.6.87.1-microsoft-standard-WSL2
machine:           x86_64
processor:         x86_64
byteorder:         little
LC_ALL:            None
LANG:              C.UTF-8
LOCALE:            ('C', 'UTF-8')

Additional Context

No response

Wegatriespython avatar Aug 21 '25 08:08 Wegatriespython

Thanks for filing this. I've removed the 'timeslice' label, as that pertains to "sub-annual time slices", and this appears to occur in a scenario that does not use them.

In Slack, someone suggested that the culprit may be df_period, which is used in the calculation of PRICE_COMMODITY: https://github.com/iiasa/message_ix/blob/cede75ca1f1143024333428c08ea37d185158ea5/message_ix/model/MESSAGE/model_solve.gms#L53-L55 and computed here: https://github.com/iiasa/message_ix/blob/3742426d378e45eb6aa2a1851184ced407cc61a5/message_ix/model/includes/period_parameter_assignment.gms#L90-L97

@Wegatriespython I believe you shared specific df_period values in Slack. Can you please also include them in the issue description, and ideally as text rather than a screenshot image?

khaeru avatar Aug 25 '25 07:08 khaeru

The value calculation for df_period is correct as per the formula. Initially I had thought the value for 2060 was suspicious as it was lower than the one for 2070 but that is a consequence of 2070 having a longer period. At 5% the length of the period outweighs the later time.

Which means that df_period does not have an implementation issue. The problem is that then the change in duration period affects the actual economic logic of the model.

For the central planner, then the decision to invest in capacity in 2060 is weighted by the fact that waiting one more period gives 10 years of returns vs 5 years of return for the investment in 2060. Hence the incentive to allow for a price shock in 2060 because the overall ROI is higher by waiting till 2070 for that capital.

This can be tested by setting up a model, one with uniform periods, and the the other with different period lengths. Since the specific price shock comes from the full model a more direct test of change in investment schedules (investment deferral) can proxy the dominant mechanism of issue here.

However since the 5->10year periods are there as a computational trick to improve performance, this behaviour is undesirable(for all the reasons mentioned in the first post on thread). Worth thinking of different approaches to address this.

year_all df_year df_period ratio
1950 1.00000e+00 5.52563e+00 5.52563
1955 7.83526e-01 4.32948e+00 5.52563
1960 6.13913e-01 3.39226e+00 5.52563
1965 4.81017e-01 2.65792e+00 5.52563
1970 3.76889e-01 2.08255e+00 5.52563
1975 2.95303e-01 1.63173e+00 5.52563
1980 2.31377e-01 1.27851e+00 5.52563
1985 1.81292e-01 1.00174e+00 5.52563
1990 1.42046e-01 7.84892e-01 5.52563
1995 1.11297e-01 6.14983e-01 5.52563
2000 8.72037e-02 4.81856e-01 5.52563
2005 6.83264e-02 3.77547e-01 5.52563
2010 5.35355e-02 2.95818e-01 5.52563
2015 4.19465e-02 2.31781e-01 5.52563
2020 3.28662e-02 1.81606e-01 5.52563
2025 2.57515e-02 1.42293e-01 5.52563
2030 2.01772e-02 1.11491e-01 5.52563
2035 1.58092e-02 8.73557e-02 5.52563
2040 1.23869e-02 6.84455e-02 5.52563
2045 9.70547e-03 5.36289e-02 5.52563
2050 7.60449e-03 4.20196e-02 5.52563
2055 5.95832e-03 3.29235e-02 5.52563
2060 4.66850e-03 2.57964e-02 5.52563
2070 2.86605e-03 3.60449e-02 12.5751
2080 1.75951e-03 2.21309e-02 12.5751
2090 1.08018e-03 1.35865e-02 12.5751
2100 6.63144e-04 8.34090e-03 12.5751
2110 4.07114e-04 5.12059e-03 12.5751

Wegatriespython avatar Aug 25 '25 15:08 Wegatriespython

OK, thanks for the update.

Which means that df_period does not have an implementation issue. The problem is that then the change in duration period affects the actual economic logic of the model.

IMO that sensitivity means there is an implementation issue. The model solution should be stable regardless of how time is discretized. In order for this to be true, sometimes it is necessary to use a more 'complicated' expression. Often a 'simple' expression of intended logic will give correct values in some cases, but not all—for instance, at a change in duration_period.

A recent example is #924; see https://gist.github.com/khaeru/59e362941f3eb84afcf08d66d3f36c72. Cell 8 in the notebook shows an expression which is more 'complicated' than what was previously in the GAMS code, but which actually gives correct results for changing period durations.

I think if we went through a similar (careful, first-principles) derivation of df_period and PRICE_COMMODITY, we would be able to reach a symbolic expression and then confirm whether the current code aligns.

khaeru avatar Aug 25 '25 16:08 khaeru

OK, thanks for the update.

Which means that df_period does not have an implementation issue. The problem is that then the change in duration period affects the actual economic logic of the model.

IMO that sensitivity means there is an implementation issue. The model solution should be stable regardless of how time is discretized. In order for this to be true, sometimes it is necessary to use a more 'complicated' expression. Often a 'simple' expression of intended logic will give correct values in some cases, but not all—for instance, at a change in duration_period.

A recent example is #924; see https://gist.github.com/khaeru/59e362941f3eb84afcf08d66d3f36c72. Cell 8 in the notebook shows an expression which is more 'complicated' than what was previously in the GAMS code, but which actually gives correct results for changing period durations.

I think if we went through a similar (careful, first-principles) derivation of df_period and PRICE_COMMODITY, we would be able to reach a symbolic expression and then confirm whether the current code aligns.

Wouldn't that imply we can discretize at no information loss ?

Agreed though with the goal of making the model robust, my intuition is that more that just the duration_period will need to be touched for this.

Wegatriespython avatar Aug 25 '25 16:08 Wegatriespython

Wouldn't that imply we can discretize at no information loss ?

To be more precise with the example in the notebook I linked, the code could not (and I guess does not) claim to model any arbitrary "continuous function $n(t)$ which is additions to capacity per unit time". The derived expression there only works in the case where $n(t)$ is piecewise exponential in form (drawn in red on all the plots). We can imagine other functions that are continuous (splines, sinusoids) for which the formulation would not capture the actual ratio of $N_b$ to $N_a$.

In that case (growth_activity_up), our documentation should remind the user that the expression provides a consistent discretization of an upper bound at a certain (piecewise) rate of exponential growth. And the examples in the notebook show that the new formulation gives the same results for that functional form and a few alternate discretizations of time.

So here, by similar reasoning (and to be clear I don't want to do the whole "(careful, first-principles) derivation" in the comments, just sketching), we could imagine a bunch of cases:

  1. 10 periods of uniform duration_period = 1 year. (Without using the "subannual time slice" feature, this is the finest resolution that the MESSAGE formulation can handle.)
  2. 1 period of duration_period = 10 years. (This is the coarsest resolution covering the same model horizon.)
  3. 2 periods of uniform duration_period = 5 years.
  4. 5 periods of uniform duration_period = 2 years.
  5. 2 periods of non-uniform duration_period = (2, 8) years.
  6. Alternate versions of (5).
  7. 3 periods of non-uniform duration_period = (1, 8, 1) years.
  8. etc.

We could then:

  • Express what we expect to be invariant across theses cases—values of CAP, ACT, etc.
  • Observe and understand how the actual marginals of COMMODITY_BALANCE_AUX relate to those.
  • Think about the notion of "commodity price" (in general) and decide how we believe those values should look: in general and in each case.

This should tell us how to derive a consistent expression (including any intermediate quantities) for PRICE_COMMODITY.

khaeru avatar Aug 29 '25 15:08 khaeru

Also, mulling over this I realized why the term "price shock" wasn't sitting right. This is a term used in some kinds of economic modeling, e.g. CGE modeling. Crucially, in such models, prices are endogenous: forcing or constraining a price to be higher or lower than a reference value is a "shock" in the sense that it forces the model to find a different equilibrium.

In MESSAGE, PRICE_COMMODITY is a reporting variable; it does not enter into the objective function and is computed only after the model is solved (per the first code snippet linked in my earlier https://github.com/iiasa/message_ix/issues/973#issuecomment-3219207566 —observe the solve statement several lines earlier). So (apparent) jumps and drops in prices are not the only or chief influence on the model solution. If there is some constraint other than COMMODITY_BALANCE with a larger marginal, the model will still choose to consume more at given indices (n,c,l,y,h) than at other indices (e.g. y) where PRICE_COMMODITY(n,c,l,y,h) is lower.

For the central planner, then the decision to invest in capacity in 2060 is weighted by the fact that waiting one more period gives 10 years of returns vs 5 years of return for the investment in 2060.

As far as I understand, remaining_capacity as derived from technical_lifetime is supposed to prevent this, by ensuring that the number of years of returns a given vintage gets is not affected by the duration_period of the period in which it's constructed. For instance, something constructed with technical_lifetime = 10 years:

  • as of the start of period '2070' = 2061-01-01 to 2070-12-31 —should only be active within that one period.
  • as of the start of period '2060' = 2056-01-01 to 2060-12-31 —should last until 2065-12-31, i.e. halfway through the period '2070', thus should have remaining_capacity(n,t,2065,2070) = 0.5.

Again, this is what I understand to be the intended representation in the GAMS code. If it has some implementation error, the solver could choose to exploit it.

The example cases in my previous comment would be the best way to investigate this.

khaeru avatar Aug 29 '25 15:08 khaeru

https://github.com/Wegatriespython/message_ix/blob/price_spike/%23973/message_ix/tests/test_nonuniform_periods.py

Got a simple reproduction of the issue going.

Structure:

Single region, single commodity (electricity), single technology (gas power plant)

  • Growing demand: 100 * (1.05^t_elapsed) GWa for t_elapsed = year_i - start_year.
  • Gas plant: 800 USD/GWa investment cost, 50 USD/GWa variable cost, 30-year lifetime
  • Capacity factor: 1.0
  • Historical year = 2015

In general price in excess capacity regimes (CAP_NEW = 0): Price = Variable Cost. In capacity investment regimes: Price = Variable cost + Periodized Capital Cost. For the special case TL =1, the periodized capital cost is the full investment cost. Therefore the price spikes are coming from the periodization of capital costs over different period lengths.

Observe the starting price of the all 10 year periods at interest rate 5 &10% (Cases 2 and 3). Now compare them to the 2040 prices in the mixed periods (Case 4 & 5) at the respective rates. They are identical. This makes sense as the first period in all the 10 year cases is exactly after a historical period 5 years in length.

Given the model treats investment cost as a lump sum amount, the periodization of capital costs is grid dependent on the discretization schema used.

To make this time invariant, we could consider using a capital recovery factor based approach. This can give us a grid invariant annualized capital cost.

Image

Difference in Period Structure and Interest Rates

1. All 5-year periods (5%)

  • Years: [2020, 2025, 2030, 2035, 2040, 2045, 2050, 2055, 2060, 2065, 2070]
  • Interest Rate: 0.05
Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 20.000000 100.000000 100.000000 95.142130 -525.720324
2025 127.628156 5.525631 127.628156 127.628156 95.093740 -411.706130
2030 162.889463 7.052261 162.889463 162.889463 95.016513 -322.320552
2035 207.892818 9.000671 207.892818 207.892818 94.893931 -252.220771
2040 265.329771 11.487391 265.329771 265.329771 94.700677 -197.219113
2045 338.635494 14.661145 338.635494 338.635494 95.538169 -155.892900
2050 432.194238 38.711749 432.194238 432.194238 95.407693 -121.979352
2055 551.601537 29.407091 551.601537 551.601537 95.198558 -95.364514
2060 703.998871 37.531728 703.998871 703.998871 94.864787 -74.458618
2065 898.500779 47.901053 898.500779 898.500779 94.334992 -58.014461
2070 1146.739979 61.135230 1146.739979 1146.739979 93.499760 -45.053387

2. All 10-year periods (5%)

  • Years: [2020, 2030, 2040, 2050, 2060, 2070]
  • Interest Rate: 0.05
Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 20.000000 100.000000 100.000000 108.885432 -1369.549264
2030 162.889463 6.288946 162.889463 162.889463 90.182678 -696.366736
2040 265.329771 10.244031 265.329771 265.329771 89.959040 -426.448619
2050 432.194238 21.686447 432.194238 432.194238 90.558142 -263.545991
2060 703.998871 38.469410 703.998871 703.998871 90.175646 -161.110995
2070 1146.739979 54.518142 1146.739979 1146.739979 89.209096 -97.848025

3. All 10-year periods (10%)

  • Years: [2020, 2030, 2040, 2050, 2060, 2070]
  • Interest Rate: 0.1
Year ACT CAP_NEW CAP Demand Price
2020 100.000000 20.000000 100.000000 100.000000 133.403381
2030 162.889463 6.288946 162.889463 162.889463 102.143920
2040 265.329771 10.244031 265.329771 265.329771 102.044335
2050 432.194238 21.686447 432.194238 432.194238 102.412176
2060 703.998871 38.469410 703.998871 703.998871 102.128049
2070 1146.739979 54.518142 1146.739979 1146.739979 100.390349

4. Mixed 5 and 10-year periods (5%)

  • Years: [2020, 2025, 2030, 2035, 2040, 2050, 2060, 2070]
  • Interest Rate: 0.05
Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 20.000000 100.000000 100.000000 95.316399 -526.683273
2025 127.628156 5.525631 127.628156 127.628156 92.925410 -402.318394
2030 162.889463 7.052261 162.889463 162.889463 95.215672 -322.996151
2035 207.892818 9.000671 207.892818 207.892818 92.725978 -246.458519
2040 265.329771 11.487391 265.329771 265.329771 108.530174 -226.019764
2050 432.194238 21.686447 432.194238 432.194238 90.558142 -263.545991
2060 703.998871 36.706344 703.998871 703.998871 90.175646 -161.110995
2070 1146.739979 53.409359 1146.739979 1146.739979 89.209096 -97.848025

5. Mixed 5 and 10-year periods (10%)

  • Years: [2020, 2025, 2030, 2035, 2040, 2050, 2060, 2070]
  • Interest Rate: 0.1
Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 20.000000 100.000000 100.000000 114.573138 -699.480467
2025 127.628156 5.525631 127.628156 127.628156 113.073745 -428.638456
2030 162.889463 7.052261 162.889463 162.889463 114.551885 -269.629974
2035 207.892818 9.000671 207.892818 207.892818 112.993035 -165.140722
2040 265.329771 11.487391 265.329771 265.329771 133.252144 -120.924217
2050 432.194238 21.686447 432.194238 432.194238 102.412176 -93.538238
2060 703.998871 36.706344 703.998871 703.998871 102.128049 -35.962988
2070 1146.739979 53.409359 1146.739979 1146.739979 100.390349 -13.629372

6. Mixed 5 and 10-year periods (0%)

  • Years: [2020, 2025, 2030, 2035, 2040, 2050, 2060, 2070]
  • Interest Rate: 0.0
Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 20.000000 100.000000 100.000000 76.666667 -383.333333
2025 127.628156 5.525631 127.628156 127.628156 76.666667 -383.333333
2030 162.889463 7.052261 162.889463 162.889463 76.666667 -383.333333
2035 207.892818 9.000671 207.892818 207.892818 76.666667 -383.333333
2040 265.329771 11.487391 265.329771 265.329771 76.666667 -383.333333
2050 432.194238 21.686447 432.194238 432.194238 76.666667 -766.666667
2060 703.998871 36.706344 703.998871 703.998871 76.666667 -766.666667
2070 1146.739979 53.409359 1146.739979 1146.739979 76.666667 -766.666667

Large Initial CAP to Neutralize CAP_NEW at Kink Year

This scenario assumes 5000 historical new capacity with a technical lifetime of 30. New capacity (CAP_NEW) is only required post-kink year.

7. Mixed 5 and 10-year periods (5%)

  • Years: [2020, 2025, 2030, 2035, 2040, 2050, 2060, 2070]
  • Interest Rate: 0.05
Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 0.000000 5000.000000 100.000000 50.000000 -276.281563
2025 127.628156 0.000000 5000.000000 127.628156 50.000000 -216.473834
2030 162.889463 0.000000 5000.000000 162.889463 50.000000 -169.612913
2035 207.892818 0.000000 5000.000000 207.892818 50.000000 -132.896155
2040 265.329771 0.000000 5000.000000 265.329771 50.000000 -104.127615
2050 432.194238 43.219424 432.194238 432.194238 90.558142 -263.545991
2060 703.998871 27.180463 703.998871 703.998871 90.175646 -161.110995
2070 1146.739979 44.274111 1146.739979 1146.739979 89.209096 -97.848025

8. All 5-year periods (5%) - Counterfactual

  • Years: [2020, 2025, 2030, 2035, 2040, 2045, 2050, 2055, 2060, 2065, 2070]
  • Interest Rate: 0.05
Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 0.000000 5000.000000 100.000000 50.000000 -276.281563
2025 127.628156 0.000000 5000.000000 127.628156 50.000000 -216.473834
2030 162.889463 0.000000 5000.000000 162.889463 50.000000 -169.612913
2035 207.892818 0.000000 5000.000000 207.892818 50.000000 -132.896155
2040 265.329771 0.000000 5000.000000 265.329771 50.000000 -104.127615
2045 338.635494 67.727099 338.635494 338.635494 95.538169 -155.892900
2050 432.194238 18.711749 432.194238 432.194238 95.407693 -121.979352
2055 551.601537 23.881460 551.601537 551.601537 95.198558 -95.364514
2060 703.998871 30.479467 703.998871 703.998871 94.864787 -74.458618
2065 898.500779 38.900382 898.500779 898.500779 94.334992 -58.014461
2070 1146.739979 49.647840 1146.739979 1146.739979 93.499760 -45.053387

Technical Lifetime = 1

This scenario sets the technical lifetime to 1, forcing per-period CAP_NEW to equal DEMAND.

9. Mixed 5 and 10-year periods (5%)

  • Years: [2020, 2025, 2030, 2035, 2040, 2050, 2060, 2070]
  • Interest Rate: 0.05
Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 100.000000 100.000000 100.000000 850.000000 -4696.786563
2025 127.628156 127.628156 127.628156 127.628156 850.000000 -3680.055170
2030 162.889463 162.889463 162.889463 162.889463 850.000000 -2883.419520
2035 207.892818 207.892818 207.892818 207.892818 850.000000 -2259.234643
2040 265.329771 265.329771 265.329771 265.329771 850.000000 -1770.169459
2050 432.194238 432.194238 432.194238 432.194238 850.000000 -2473.704582
2060 703.998871 703.998871 703.998871 703.998871 850.000000 -1518.640028
2070 1146.739979 1146.739979 1146.739979 1146.739979 850.000000 -932.313241

10. All 5-year periods (5%)

  • Years: [2020, 2025, 2030, 2035, 2040, 2045, 2050, 2055, 2060, 2065, 2070]
  • Interest Rate: 0.05
Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 100.000000 100.000000 100.000000 850.000000 -4696.786563
2025 127.628156 127.628156 127.628156 127.628156 850.000000 -3680.055170
2030 162.889463 162.889463 162.889463 162.889463 850.000000 -2883.419520
2035 207.892818 207.892818 207.892818 207.892818 850.000000 -2259.234643
2040 265.329771 265.329771 265.329771 265.329771 850.000000 -1770.169459
2045 338.635494 338.635494 338.635494 338.635494 850.000000 -1386.974090
2050 432.194238 432.194238 432.194238 432.194238 850.000000 -1086.730492
2055 551.601537 551.601537 551.601537 551.601537 850.000000 -851.481776
2060 703.998871 703.998871 703.998871 703.998871 850.000000 -667.158252
2065 898.500779 898.500779 898.500779 898.500779 850.000000 -522.735948
2070 1146.739979 1146.739979 1146.739979 1146.739979 850.000000 -409.577293

Wegatriespython avatar Aug 30 '25 19:08 Wegatriespython

Using the CRF formula : for r >0 $$\text{CRF}(r, n) = \frac{r (1+r)^n}{(1+r)^n - 1}$$

r = 5%, n = 30,

CRF ~= 0.0654

For Inv cost = 800, annualized cost = 800 * 0.0654 = 52.3

For the same scenario set inv cost = eps (so that membership in inv techs is not lost) and set fixed cost = 52.3

Years: [2020, 2025, 2030, 2035, 2040, 2050, 2060, 2070], Interest Rate: 0.05

Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 20.000000 100.000000 100.000000 102.041148 -563.841757
2025 127.628156 5.525631 127.628156 127.628156 102.041148 -441.784771
2030 162.889463 7.052261 162.889463 162.889463 102.041148 -346.149927
2035 207.892818 9.000671 207.892818 207.892818 102.041148 -271.217526
2040 265.329771 46.013693 265.329771 265.329771 102.041148 -483.723554
2050 432.194238 39.693293 432.194238 432.194238 102.041148 -296.964301
2060 703.998871 70.399887 703.998871 703.998871 102.041148 -182.310320
2070 1146.739979 44.274111 1146.739979 1146.739979 102.041148 -111.922722

Notice that the model is now effectively retiring previous capacity every period, and choosing to start fresh. Exploiting an early retirement arbitrage.

Deactivating early retirement we get :

Years: [2020, 2025, 2030, 2035, 2040, 2050, 2060, 2070], Interest Rate: 0.05

Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 20.000000 100.000000 100.000000 102.041148 -563.841757
2025 127.628156 5.525631 127.628156 127.628156 102.041148 -441.784770
2030 162.889463 7.052261 162.889463 162.889463 102.041148 -346.149927
2035 207.892818 9.000671 207.892818 207.892818 102.041148 -271.217526
2040 265.329771 11.487391 265.329771 265.329771 102.041148 -483.723554
2050 432.194238 21.686447 432.194238 432.194238 102.041148 -296.964301
2060 703.998871 36.706344 703.998871 703.998871 102.041148 -182.310320
2070 1146.739979 53.409359 1146.739979 1146.739979 102.041148 -111.922722

Which is virtually the same new capacity investment schedule as case 4.

So effectively the issue is resolved by removing the lump sum INV cost from model objective, replacing it with the annuity A = CRF(r,n) * INV, and adding an early retirement penalty = PV sum of the remaining annual payments per CRF.

Pre-liminary results -> the price spike at 2060 vanishes.

Image

Wegatriespython avatar Sep 12 '25 11:09 Wegatriespython

This is a solid, systematic investigation of the issue and a potential solution, and I think worth presenting and discussing during the MESSAGE weekly meeting some week.

So effectively the issue is resolved by removing the lump sum INV cost from model objective, replacing it with the annuity A = CRF(r,n) * INV, and adding an early retirement penalty = PV sum of the remaining annual payments per CRF.

Just to confirm I've understood the diagnosis, here's what I hope is the same thing stated in different terms:

  • The objective function value is computed as a sum over the year dimension of COST_NODAL: https://github.com/iiasa/message_ix/blob/33edac8acb80ed22300519b2882a4a6e438fae89/message_ix/model/MESSAGE/model_core.gms#L340-L341
  • COST_NODAL is defined here: https://github.com/iiasa/message_ix/blob/33edac8acb80ed22300519b2882a4a6e438fae89/message_ix/model/MESSAGE/model_core.gms#L389-L470
  • By "the [current] lump sum inv_cost [in the] model objective", you refer to the fact that inv_cost here is multiplied by CAP_NEW in the same/single period in which a particular technology vintage is constructed (the "vintage period", index year).
  • If meeting demand in later periods requires increasing a single value CAP_NEW(node, tec, year), this affects OBJ via an increase in COST_NODAL(node, year) for that vintage period. Conversely, there is no connection at all between CAP_NEW(…, year)/inv_cost(…, year) and future periods' COST_NODAL(…, year+1) etc.
  • The marginal of COMMODITY_BALANCE_AUX.M(node,commodity,level,year,time) (used to compute PRICE_COMMODITY, same dimensionality) reflects the potential reduction in OBJ via that COST_NODAL(…, year), only for that vintage period (year) and not for future periods (year+1 et seq) in which output from that same technology vintage will also be used.
  • Thus the prices in those future periods are artificially low, and in the vintage period artificially high.

khaeru avatar Sep 12 '25 13:09 khaeru

I'm not sure, so the picture as I see it is at the optimal solution the zero profit condition holds and the model balances all its costs across the horizon(allowing our use of crf). This implicitly requires a capital recovery schedule. However since cap new is a lump sum via the objective and the objective is on an uneven grid, the model's capital recovery schedule inherits the approximation error of the grid.

So for instance if we set the resolution at 1 year and same everything else we get an annual price series closer to the calculated CRF value of 102.

Now I can't for sure say why there is a that small difference between the model's endogenous calculation and 102. But the idea is we are giving the model a nicer pre computed method for investment cost annualization than its actual grid. This is what removes the artefact.

Years: [2020, 2021, 2022, 2023, 2024, 2025, 2026, 2027, 2028, 2029, 2030, 2031, 2032, 2033, 2034, 2035, 2036, 2037, 2038, 2039, 2040, 2041, 2042, 2043, 2044, 2045, 2046, 2047, 2048, 2049, 2050, 2051, 2052, 2053, 2054, 2055, 2056, 2057, 2058, 2059, 2060, 2061, 2062, 2063, 2064, 2065, 2066, 2067, 2068, 2069, 2070], Interest Rate: 0.05

Year ACT CAP_NEW CAP Demand Price Dual_COMM_BAL_AUX
2020 100.000000 100.000000 100.000000 100.000000 99.632543 -99.632543
2021 105.000000 5.000000 105.000000 105.000000 99.621985 -94.878081
2022 110.250000 5.250000 110.250000 110.250000 99.610387 -90.349558
2023 115.762500 5.512500 115.762500 115.762500 99.597649 -86.036194
2024 121.550625 5.788125 121.550625 121.550625 99.583662 -81.927725
2025 127.628156 6.077531 127.628156 127.628156 99.568305 -78.014372
..
2067 990.597109 58.085664 990.597109 990.597109 98.061152 -9.899196
2068 1040.126965 60.989947 1040.126965 1040.126965 97.842668 -9.406800
2069 1092.133313 64.039444 1092.133313 1092.133313 97.603731 -8.936980
2070 1146.739979 67.241417 1146.739979 1146.739979 97.342524 -8.488631

Wegatriespython avatar Sep 16 '25 12:09 Wegatriespython