Add commodity flow of capacity variables

[Jihoon]

This PR cleans up and supersedes #387.

PR checklist

• Continuous integration checks all ✅
• Add or expand tests; coverage checks both ✅
• Add, expand, or update documentation.
• Update release notes.

Things to review

• The formulation and the documentation is explanatory about the new changes that are made.
• Check if there are any issues that can possibly arise from the new formulation. (any cases where it might not work as expected).
• Building and solving a scenario with the new formulation without any erorrs (can be the Westeros tutorial).

Expected review time: One hour or one and a half hour

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 96.1%. Comparing base (9fb0b97) to head (d77ab2b).
⚠️ Report is 47 commits behind head on main.

Additional details and impacted files

@@          Coverage Diff          @@
##            main    #451   +/-   ##
=====================================
  Coverage   96.1%   96.1%           
=====================================
  Files         53      54    +1     
  Lines       5030    5097   +67     
=====================================
+ Hits        4834    4901   +67     
  Misses       196     196           

Files with missing lines	Coverage Δ
message_ix/core.py	97.4% <ø> (ø)
message_ix/models.py	99.2% <100.0%> (+<0.1%)	⬆️
message_ix/report/__init__.py	100.0% <ø> (ø)
message_ix/testing/__init__.py	99.7% <ø> (-0.1%)	⬇️
message_ix/tests/model/message/test_cap_comm.py	100.0% <100.0%> (ø)
...age_ix/tests/test_feature_bound_activity_shares.py	100.0% <100.0%> (ø)
message_ix/tests/test_models.py	100.0% <ø> (ø)
message_ix/tests/test_report.py	100.0% <100.0%> (ø)
message_ix/tests/test_tutorials.py	59.6% <ø> (ø)

... and 1 file with indirect coverage changes

[khaeru]

#452 is modifying the CI setup. Please restart the tests here (manually or via rebase) once that one is merged, hopefully within an hour.

[khaeru]

I edited the PR description to re-add the checklist that's auto inserted by the template. These items, particularly documentation, are essential—even for minor PRs. This one is not minor, because it alters the core model formulation.

At the absolute minimum, when completing the item “Update release notes,” mention the parameters, variables, and equations that are added/affected by the additions here. Then, users might have a thin hope of understanding that changes have been made, and might try to puzzle them out or even use them.

Much better would be to update the inline documentation in the GAMS files, or add a new file in docs/, with plain-language description. For a good example, see #190 by @behnam-zakeri: https://github.com/iiasa/message_ix/pull/190/files, in particular

• additions in model_core.gms and parameter_def.gms.
• how the corresponding item in the release notes links to these descriptions.

As a reminder, if time feels tight, there is always the option to task a specific colleague with some or all of this work.

[khaeru]

Postponed to v3.4 because:

• PR checklist is incomplete.
• 3 requested reviews not done.
• Requested changes from 1 review not completed.

[LauWien]

I've added the WIP tag, since there is still more time/investigation needed to finish up the PR, as discussed in today's MESSAGEix meeting. Due to that and the following

• PR checklist is incomplete,
• 3 requested reviews not done,
• Requested changes from 1 review not completed,

we postpone to v3.5.

[GamzeUnlu95]

map_tec_lifetime issue for the first year scrap release equation:
If we say the first model year is 700 and the vintage is 690:
According to this equation which is valid only for the first model year, based on the conditions specified here
year = 700, year2 = 690. In this case for a technology that has vintage 690 map_tec_lifetime(location,tec,690,690) becomes a condition that needs to be satisfied. And map_tec_lifetime does not have this entry for vintage years. In this case there will be no scrap released for the first model year.

As another example, if we say the vintage year is 680 and the lifetime is 20 years, the only way to make sure this plant will retire in the first model period (700) is to check map_tec_lifetime and see the entry for (680, 690). Otherwise, we don't know when it is supposed to retire. We either need to make this historical vintage year addition to map_tec_lifetime or change the equations to satisfy the necessary condition.

[Jihoon]

Just to add a bit on Gamze's comment above, for the capacity supposed to end its activity in 690 and retire in 700, the historical capacity that existed in 690 cannot be taken into account, because the map_tec_lifetime condition does not pass for the period pair of (year2, year) = (690, 700). map_tec_lifetime(location,tec,vintage,year2) is only defined for year2 >= fmy = 700.

[GamzeUnlu95]

Even if the necessary map_tec_lifetime is added manually in GAMS, the if conditions are not satisfied and no scrap is released in the first year. Not sure why this is happening... The equations printed in the .lst format are attached. The gdx files are in: I:\unlu\power_sector_integration\24.05

MESSAGE_run_test_2405.txt

[GamzeUnlu95]

Reformulation of the first year scrap release equation:
The remaining capacity seems to be a binary value, does it make sense to perform a subtraction such as
( remaining_capacity(node,tec,vintage,year2) - CAP(location,tec,vintage,year) )?

[GamzeUnlu95]

Here we can use something like:

• output by all new capacity of technologies located at 'location' sending to 'node' and 'time' distributed over years of periods

        output_cap_ret(location,tec,vintage,node,commodity,level,time)                                                 
        ( historical_new_capacity(node,tec,vintage) - CAP(location,tec,vintage,year)/ duration_period(year) )  

• input by all new capacity of technologies located at 'location' taking from 'node' and 'time' distributed over years of periods

        - input_cap_ret(location,tec,vintage,node,commodity,level,time)                                               
        (historical_new_capacity(node,tec,vintage) - CAP(location,tec,vintage,year)/ duration_period(year))) 

We might need to sum the historical_new_capacity over vintage years (if there are two different vintages with different life times both retiring in the first model year...)

[Jihoon]

Another suspected issue here is:
For the capacity (e.g. vintage=740) ending its activity in e.g. year2=750, because the condition seq_period(750,760) still passes, this part of the equation keeps the term -CAP(location,tec,740,760) (in line 613), which we expect to be 0 by definition. But differently from our expectation, it appears that while solving the LP model, a spurious non-zero value is assigned to this CAP(location,tec,740,760).

[volker-krey]

The issue related to the historical vintage 690 of coal_ppl seems to be a result of the missing technical_lifetime parameter of that vintage in the Westeros tutorial (see #609).

The problem of non-zero CAP values of technologies after the end of their lifetime is due to the fact that there is no equation that would require CAP to be zero after the end of the technology lifetime. The equation CAPACITY_MAINTENANCE just requires that CAP in a period is lower than in a previous period. This usually does not lead to problems, but in this extension there might be a value in delaying commodity release to better balance a "scrap-type" commodity equation that is set to equal. Usually, there is no cost penalty and also no other benefits as other parameters have not been defined for periods beyond the lifetime.

To avoid this behavior, we may need to introduce an additional equation of the following kind:

CAPACITY_MAINTENANCE2(node,inv_tec,vintage,year)$( map_tec(node,inv_tec,vintage) AND NOT map_tec_lifetime(node,inv_tec,vintage,year) 
AND NOT first_period(year) AND year_order(vintage) < year_order(year))..
    CAP(node,inv_tec,vintage,year) =L= 0 ;

[khaeru]

Another test failure:

FAILED message_ix/tests/test_macro.py::test_calc_price_zero[jdbc] - AssertionError: Regex pattern did not match.

This boils down to the following: on main, the PRICE_COMMODITY value in the solved Westeros scenario look like:

    node    commodity      level  year  time            lvl  mrg
Westeros  electricity  secondary   710  year  4.940656e-324  0.0
Westeros  electricity  secondary   720  year  4.940656e-324  0.0
Westeros  electricity      final   700  year  4.940656e-324  0.0
Westeros        light     useful   700  year  4.940656e-324  0.0
Westeros        light     useful   710  year  4.940656e-324  0.0
Westeros        light     useful   720  year   5.000000e+00  0.0

On this branch, we instead see:

    node    commodity      level  year  time            lvl  mrg
Westeros  electricity  secondary   700  year  4.940656e-324  0.0
Westeros  electricity  secondary   710  year  4.940656e-324  0.0
Westeros  electricity      final   700  year  4.940656e-324  0.0
Westeros  electricity      final   710  year  4.940656e-324  0.0
Westeros  electricity      final   720  year  4.940656e-324  0.0
Westeros        light     useful   720  year   5.000000e+00  0.0

The missing values for (light, useful, 700), (light, useful, 710) cause a preliminary calculation to fail, before reaching the main focus of the test.

Some thing I observe:

• In neither case—not on main or this branch—does PRICE_COMMODITY include all of the "active" combinations of (node, commodity, level, year, time)—that is, the indices for which there are non-zero in- or outflows of commodities.

• The descriptive text for PRICE_COMMODITY says

  commodity price (derived from marginals of COMMODITY_BALANCE constraint)

  • This was set in a commit from 2018: https://github.com/iiasa/message_ix/blame/9fb0b97867f9234d98b4e43e158d7af9b30c42f4/message_ix/model/MESSAGE/model_core.gms#L148
  • That actually precedes commit 753e911 that changed COMMODITY_BALANCE from an equation to a macro (the current PR is changing it back).

• It is defined here:

  message_ix/message_ix/model/MESSAGE/model_solve.gms

  Lines 52 to 55 in 9fb0b97

  	* assign auxiliary variables PRICE_COMMODITY and PRICE_EMISSION for reporting
  	PRICE_COMMODITY.l(node,commodity,level,year,time) =
  	( COMMODITY_BALANCE_GT.m(node,commodity,level,year,time) + COMMODITY_BALANCE_LT.m(node,commodity,level,year,time) )
  	/ df_period(year) ;

• I note that the new COMMODITY_BALANCE_AUX equation does have marginal values present for all "active" (n, c, l, y, h) combinations. I see values like:

$ gdxdump X.gdx

Equation COMMODITY_BALANCE_AUX(node,commodity,level,year_all,time) Auxiliary equation for calculation of commodity balance /
'Westeros'.'light'.'useful'.'720'.'year'.M -8.93317663555801,

Equation COMMODITY_BALANCE_GT(*,*,*,*,*) commodity supply greater than or equal demand /
'Westeros'.'light'.'useful'.'720'.'year'.M 8.93317663555801,

Equation COMMODITY_BALANCE_LT(*,*,*,*,*) commodity supply lower than or equal demand / /;

• In other words, the value of COMMODITY_BALANCE_AUX.M appears to be the negative of (COMMODITY_BALANCE_GT.M + COMMODITY_BALANCE_LT.M), the quantity

Based on these observations, I will try to revise PRICE_COMMODITY to refer to COMMODITY_BALANCE_AUX instead of COMMODITY_BALANCE_{GT,LT}. This appears to be what existed before.