iiasa / message_ix

The integrated assessment and energy systems model MESSAGEix
https://docs.messageix.org
Apache License 2.0
117 stars 152 forks source link

parameter `bound_activity_lo` does not bound `ACT` but `CAP` #860

Closed SongminYu closed 2 months ago

SongminYu commented 2 months ago

I am developing a power system model with MESSAGE. There are 31 nodes (provinces in China) and each of them has following technologies: coal_ppl, gas_ppl, nuclear_ppl, hydro_ppl, biomass_ppl, solar_ppl, and wind_ppl. I want to optimize the ACT and CAP pathway of all technologies in all nodes. The firstmodelyear in my model is 2025 and it optimizes until 2060.

The technologies are parametrized with historical_new_capacity from 1990 to 2020. I am following the westeros_historical_new_capacity tutorial to force the use of coal_ppl by using the bound_activity_lo parameter. Because I want to make sure that the coal_ppl capacity are not retired too early despite the relatively higher generation cost, and they should still be used at least at the level of bound_activity_lo. Given the vintage and lifetime, they should last until 2045.

However, as shown in this excel (bound_activity_lo.xlsx), the bound_activity_lo parameter is not limiting the ACT but the CAP. The calculated bound_activity_lo / CAP = capacity_factor area showed that, the CAP is calculated with capacity_factor = 0.95. However, this is the uplimit of the capacity_factor as I defined in the input. The real capacity_factor can be much lower due to the competition with other technologies. As a result, the ACT is significantly lower than the bound_activity_lo.

In summary, I think the bound_activity_lo is not working as the name indicates. The example in the westeros_historical_new_capacity tutorial is not reflecting this, because the gas_ppl in the example is not cheap enough. In my case, the model chooses to investment cheap wind_ppl and hydro_ppl starting from the firstmodelyear and leave coal_ppl working at very low real capacity_factor. The bound_activity_lo is not working as expected.

Versions

Output of message-ix show-versions ``` ```

The command message-ix show-versions shows message_ix: 0.0.0 but this should not be true. I don't know why. It should also be 3.8.0 as shown in the screenshot below.

Screenshot
glatterf42 commented 2 months ago

Hi @SongminYu, thank you for your detailed error description :) Looking at your bound_activity_lo.xlsx, I'm wondering how you got this. Because just looking at the numbers, what you call bound_activity_lo seems to always be greater than or equal to what you call ACT, so it almost seems like these two got mixed up. This would also make sense for the calculation you mention, because indeed, assuming $\mathrm{duration\_time} = 1$, our docs show that $\mathrm{ACT} <= \mathrm{capacity\_factor} \cdot \mathrm{CAP}$, and if there is a minimum required ACT and CAP is not free to build, I guess the equation will in practice often evaluate with $=$.

This is not to say that there's no error on our side, maybe we mix up these parameters when reading them from the model solution file. However, I would expect that one of our internal researchers would have already run into the same issue, too.

SongminYu commented 2 months ago

Hi @glatterf42, thanks for your swift reply!

The bound_activity_lo numbers are my input. The ACT and CAP are results from MESSAGE. I used following code to export them:

reporter = message_ix.Reporter.from_scenario(scenario)
for data_name in ["ACT", "CAP"]:
    reporter.write(reporter.get(data_name), f'{data_name}.csv')

Then, I just selected the numbers of coal_ppl from ACT and CAP, to compare them with my bound_activity_lo input. The duration_time in my model is 5 years.

This would also make sense for the calculation you mention, because indeed, assuming duration_time=1, our docs show that ACT<=capacity_factor⋅CAP, and if there is a minimum required ACT and CAP is not free to build, I guess the equation will in practice often evaluate with =.

Yes, when CAP is not free to build in such case, the = holds, meaning that CAP = bound_activity_lo/capacity_factor. But, ACT is not bound by bound_activity_lo. I assigned coal_ppl the capacity_factor 0.95, but as I understand, this is the uplimit in the optimization. The real operation hours of coal_ppl is usually lower than 8760*0.95=8322 hours.

I haven't digged into the GAMS code as I am less familiar with it. Do you have some suggestions on how to debug this? Thanks!

glatterf42 commented 2 months ago

Looking at your data file again and comparing with our docs for bound_activity_lo, this parameter usually has the dimensions node_loc | tec | year_act | mode | time, but your data file just shows the first three, I presume? Similarly, ACT is defined for node_loc | technology | year_vtg | year_act | mode | time | lvl | mrg, and it's not clear which information you are showing in your data file. Is it summed over all year_vtg? Please provide us with complete data, which you can extract e.g. via scenario.var("ACT") and scenario.par("bound_activity_lo"), or, better yet, with a complete minimal working example of how you add these data to the scenario, how you solve it, and how you extract your results.

SongminYu commented 2 months ago

Hi @glatterf42, thanks! I extracted the data and did a rough pivot analysis: bound_activity_lo_pivot.xlsx

The pivot_analysis sheet summarizes data from the other four sheets. I hope the red notes are clear enough...

One thing I also noticed while preparing the excel file: the bound_activity_lo_param and bound_activity_lo_input are not exactly equal... There is a 2% difference... I don't understand yet. I am looking for the reason. But, it is still clear that the bound_activity_lo is limiting the CAP rather than ACT.

I want to show you the code and my workflow but it is a bit too complicated... I invited you to my private repo. After setting up the env, you can run projects/power_cement_storage_48h/main.py. This should produce the same output as shown above...

glatterf42 commented 2 months ago

Hi @SongminYu, thanks for the additional information and access to your private repo! That looks like it was a lot of work.

Unfortunately, I think the premises of why you put in the work are invalid. In your private message, you said you wanted to wrap message_ix because you're more familiar with PyCharm and because you wanted to use Excel files as input. message_ix does not care about which editor you're using. We have several people in our team using PyCharm with our regular code base and it's working fine. We also already provide the scenario.read_excel() function so that input data can be read from Excel. As far as I know, this is also in frequent use within our team.

By now, your private repository might contain several additions to the code that our regular code does not provide -- and we would be glad to have them! The best way to go about this would be to start working on a fork of message_ix, tweaking all the things that you need to tweak in order to make it work, and then open a pull request to collaborate with us in bringing your changes to the regular code base. What we have now, are several hundred lines of code that are not really documented and do what you want -- almost. To troubleshoot private code like this would require us to spend considerable time understanding all nuances of the things that have changed, before we could being assessing what has gone wrong. Unfortunately, we don't have a magic tool to tell us why things are not working. For our regular code base, we do have a suite of tests that checks that everything is working as expected (which makes me more confident that bound_activity_lo is likely working as intended), but your private repo does not have that.

So I may find some time to study your code in detail, but I don't yet know when. I recommend to still do what would have been best from the beginning: please try to make your changes as part of the regular message code and open a PR; this will get a much quicker review. Plus, it will mean that everyone can benefit from your improvements!

SongminYu commented 2 months ago

Hi @glatterf42, thanks for your reply! Yes, I know it is too much to ask you to study my code. They work for my formatting of the input tables only, and convert them into the message format. If you can run projects/power_cement_storage_48h/main.py, there will be an output folder generated (ignored by git) and you can find all the final input parameters formatted in the message way as shown below.

Screenshot 2024-07-03 at 21 18 10

Then, you can see if some parameters are wrong.

I am not asking you to do it soon... Just letting you know that you can do this before looking at my additional code... I am aware of the scenario.read_excel() function, but I felt it more familiar to prepare the data in my wide format, as I also do for my other models.

Taking the param_bound_activity_lo as an example, the converted parameter dataframe is like below:

Screenshot 2024-07-03 at 21 20 09

Please note that, these "param_xxx" files are exported by the scenario object, after it being initialized. So, the "conversion process" do produce parameter dataframes that can be successfully imported into the message scenario object.

Thanks again for your patient and swift answers!

SongminYu commented 2 months ago

Hi @glatterf42, some good updates: I don't understand why but "refactoring adding timeslices" fixed my problem of bound_activity_lo...

Thanks to @behnam-zakeri for directing me to PR#852. I refactored my previous code for adding timeslices. It looks cleaner and also solved the bug in this thread, interestingly. So, there should be no bug in the parameter bound_activity_lo. Sorry for creating the confusion. I will close this thread now.