Issue with Implementing Storage Feature in MESSAGEix Framework

[shreeyashn20]

I am encountering challenges while building a national-level model using the MESSAGEix framework. Specifically, I am attempting to integrate the storage feature into the model. To accomplish this, I have referred to test_feature_storage.py, as well as the Central Asia and Brazil models. Despite these references, I have been unable to implement the storage feature correctly.

The issue is that the model is not showing any activity related to the battery, charging, or discharging processes. I have modified test_feature_storage.py from using 4 time slices to 2 time slices, dividing the year into peak and off-peak periods. However, the model does not seem to store energy during the off-peak slice as expected. If I try to force bound activity on power plants, either the bounds get relaxed or model gives infeasibility.

I have attached the relevant code and output results below for reference. Any insights or suggestions for resolving this issue would be greatly appreciated. Thank you in advance.

Library

import pandas as pd
import numpy as np
import ixmp
import message_ix
import itertools
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot')

from io import StringIO
from message_ix import make_df
#from message_ix_models.util import broadcast
from message_ix.utils import make_df

pd.set_option('display.max_rows', None)
pd.set_option('display.max_columns', None)
pd.set_option('display.width', None)
pd.set_option('display.max_colwidth', None)

## Model Platform, Scenario, Spatial, Horizon, Interest Rate, Mode

mp = ixmp.Platform()

scenario = message_ix.Scenario(mp, 
                               model="Storage Model", 
                               scenario="baseline", 
                               version="new")

scenario.add_spatial_sets({"country":"XYZ"})

model_horizon = [2020]
scenario.add_horizon(
    year= model_horizon,
    firstmodelyear=model_horizon[0]
)


scenario.add_set("year", model_horizon)
scenario.add_set("type_year", model_horizon)

year_df = scenario.vintage_and_active_years()
vintage_years, act_years = year_df["year_vtg"], year_df["year_act"]
year_vtg_df=year_df[["year_vtg"]]
year_act_df=year_df[["year_act"]]

scenario.add_par("interestrate", model_horizon, value=0.05, unit='-')

scenario.add_set("mode", ["standard"])

scenario.add_set("time",["peak","offpeak"])
scenario.add_set("lvl_temporal","division")
scenario.add_set("map_temporal_hierarchy", ["division", "peak", "year"])
scenario.add_set("map_temporal_hierarchy", ["division", "offpeak", "year"])
scenario.add_par("duration_time", "peak", 0.4, "%")
scenario.add_par("duration_time", "offpeak", 0.6, "%")

# All parameters with at least one sub-annual time index
parameters = [p for p in scenario.par_list() if "time" in scenario.idx_sets(p)]

# Those parameters with time index that are not empty in our model
[p for p in parameters if not scenario.par(p).empty]

# A function for adding sub-annual data to a parameter
def yearly_to_season(scenario, parameter, data, filters=None):
    if filters:
        old = scenario.par(parameter, filters)
    else:
        old = scenario.par(parameter)
    scenario.remove_par(parameter, old)

    # Finding time related indexes
    time_idx = [x for x in scenario.idx_names(parameter) if "time" in x]
    for h in data.keys():
        new = old.copy()
        for time in time_idx:
            new[time] = h
        new["value"] = data[h] * old["value"]
        scenario.add_par(parameter, new)

scenario.add_set("level","storage")
scenario.add_set("technology",["battery","charger","discharger","dummer"])
scenario.add_set("commodity", ["elec_stored","dummy1","dummy2"])
# Specifying storage reservoir technology
scenario.add_set("storage_tec", "battery")
# Specifying storage level
scenario.add_set("level_storage", "storage")
# Adding time sequence
scenario.add_par("time_order",["division","peak"],1,"-")
scenario.add_par("time_order",["division","offpeak"],2,"-")

# Adding mapping for storage and charger/discharger technologies
for tec in ["charger", "discharger"]:
        scenario.add_set(
            "map_tec_storage",
            ["XYZ", tec, "standard", "battery", "standard", "storage", "elec_stored", "division"],
        )


data_storage_self_discharge = """
node,technology,mode,level,commodity,time,value,unit
XYZ,battery,standard,storage,elec_stored,peak,0.05,-
XYZ,battery,standard,storage,elec_stored,offpeak,0.05,-
"""

df_storage_self_discharge = pd.read_csv(StringIO(data_storage_self_discharge))

year_act_df['key'] = 0
df_storage_self_discharge['key'] = 0
merged_storage_self_discharge = pd.merge(year_act_df, df_storage_self_discharge, on='key').drop('key', axis=1)
merged_storage_self_discharge = merged_storage_self_discharge.rename(columns={'year_act': 'year'})

scenario.add_par("storage_self_discharge",merged_storage_self_discharge)


data_storage_initial = """
node,technology,mode,level,commodity,time,value,unit
XYZ,battery,standard,storage,elec_stored,peak,0.05,-
XYZ,battery,standard,storage,elec_stored,offpeak,0.05,-
"""

df_storage_initial = pd.read_csv(StringIO(data_storage_initial))

year_act_df['key'] = 0
df_storage_initial['key'] = 0
merged_storage_initial = pd.merge(year_act_df, df_storage_initial, on='key').drop('key', axis=1)
merged_storage_initial = merged_storage_initial.rename(columns={'year_act': 'year'})

scenario.add_par("storage_initial",merged_storage_initial)

scenario.par("storage_initial")

scenario.add_set("commodity", "electricity")
scenario.add_set("level", "final")
scenario.add_set("mode", "standard")
scenario.add_set("technology", ["wind_ppl",
                                "gas_cc_ppl",
                               ])

data_output = """
node_loc,technology,mode,node_dest,commodity,level,time,time_dest,value,unit
XYZ,gas_cc_ppl,standard,XYZ,electricity,final,year,year,1.0,-
XYZ,wind_ppl,standard,XYZ,electricity,final,year,year,1.0,-
XYZ,battery,standard,XYZ,dummy1,final,year,year,1.0,-
XYZ,charger,standard,XYZ,elec_stored,storage,year,year,1.0,-
XYZ,discharger,standard,XYZ,electricity,final,year,year,1.0,-
XYZ,dummer,standard,XYZ,dummy2,final,year,year,1.0,-
"""

df_output = pd.read_csv(StringIO(data_output))

year_df['key'] = 0
df_output['key'] = 0

merged_output = pd.merge(year_df, df_output, on='key').drop('key', axis=1)

scenario.add_par("output",merged_output)

data_input = """
node_loc,technology,mode,node_origin,commodity,level,time,time_origin,value,unit
XYZ,battery,standard,XYZ,dummy2,final,year,year,1.0,-
XYZ,charger,standard,XYZ,electricity,final,year,year,1.0,-
XYZ,discharger,standard,XYZ,elec_stored,storage,year,year,1.0,-
"""

df_input = pd.read_csv(StringIO(data_input))

year_df['key'] = 0
df_input['key'] = 0

merged_input = pd.merge(year_df, df_input, on='key').drop('key', axis=1)

scenario.add_par("input",merged_input)

data_var_cost = """
node_loc,technology,mode,time,value,unit
XYZ,gas_cc_ppl,standard,year,2,USD/GWa
XYZ,wind_ppl,standard,year,0,USD/GWa
XYZ,battery,standard,year,0,USD/GWa
XYZ,charger,standard,year,0.2,USD/GWa
XYZ,discharger,standard,year,0.3,USD/GWa
"""

df_var_cost = pd.read_csv(StringIO(data_var_cost))

year_df['key'] = 0
df_var_cost['key'] = 0

merged_var_cost = pd.merge(year_df, df_var_cost, on='key').drop('key', axis=1)

scenario.add_par("var_cost",merged_var_cost)

data_bound_activity = """
node_loc,technology,year_act,mode,time,value,unit
XYZ,wind_ppl,2020,standard,year,0.05,GWa
XYZ,gas_cc_ppl,2020,standard,year,0.05,GWa
"""
df_bound_activity = pd.read_csv(StringIO(data_bound_activity))

scenario.add_par("bound_activity_up",df_bound_activity)

# Modifying the demand for each season
demand_data = {"peak": 0.70, "offpeak": 0.30}
yearly_to_season(scenario, "demand", demand_data)
bound_data = {"peak": 0.30, "offpeak": 0.70}
# Let's look at demand now
scenario.par("demand")
# Modifying input and output parameters for each season
fixed_data = {"peak": 1, "offpeak": 1}
# Output
yearly_to_season(scenario, "output", fixed_data)
# Input
yearly_to_season(scenario, "input", fixed_data)
# Modifying growth rates for each season
# yearly_to_season(scenario, "growth_activity_up", bound_data)
# Modifying variable cost
yearly_to_season(scenario, "var_cost", fixed_data)
# Modifying variable cost
yearly_to_season(scenario, "bound_activity_up", fixed_data)

scenario.commit("Basic Model")
scenario.set_as_default()
scenario.solve(solve_options={'lpmethod': '4'})

scenario.var("OBJ")

scenario.var("ACT")

act_no_stor = scenario.var("ACT", {"technology": "gas_cc_ppl"})["lvl"].sum()

activitydf = scenario.var("ACT")
pivotact = activitydf.pivot_table('lvl',index='technology', columns='year_act',
                    aggfunc= 'sum')

pivotact = pivotact.round(2)

print("ACTIVITY - Unit GWa")
pivotact.round(2)

# Multiply the sum values by 8.76
pivotactTWh = activitydf.pivot_table('lvl',index='technology', columns='year_act',
                    aggfunc= lambda x: sum(x) * 8.76)

pivotactTWh = pivotactTWh.round(2)

# Print the modified pivot table
print("ACTIVITY - Unit TWh")
pivotactTWh.round(2)

mp.close_db()

[glatterf42]

Hi @shreeyashn20, thanks for opening this discussion :)
I don't immediately know what you mean with the "Central Asia" and "Brazil" models. Could you please provide links to both so that we can be sure?

Next, I'm wondering if you have found our tutorial for the storage feature? If not, then working through it may clarify some additional things.

Now, I don't perform a lot of modelling activity myself, so I don't know intuitively which parts of your code might cause the issue (thanks for providing all of that, this is very helpful!). However, I can go through the code, comparing it to other uses of the feature and try to spot differences.

• First of, I'm wondering if you truly intend to build a model just for the year 2020. That does not seem very useful and it doesn't give the model a lot of room to optimize anything. However, since our test_feature_storage model is also just built on one year, this should not be the cause of your issue.
• Next, I notice you're defining duration_time like this: scenario.add_par("duration_time", "peak", 0.4, "%"), but in the tutorial it is explained (and also used similarly in the test) that duration times are given relative to 1 and should add up to 1, so you should likely use scenario.add_par("duration_time", ["peak"], 0.4, "-").
• The way you set up storage_self_discharge seems error-prone to me. Please ensure that your functionality does the same as the corresponding lines in our test, in particular, the unit should be "%" here.
  Essentially, you have already spelled out exactly what you want to add with the scenario.add_par(...) call except for one column: year, which should be 2020 for both rows in your case. So instead of employing lots of pandas to merge with year_act_df, you could simply add this column in your multi-line string that you read with pd.read_csv(...).
• The same is true for storage_initial: please use "%" as the unit and make sure your functionality is the same as in the test. Also, please note that in our test, the storage_initial is only added for the very first time slice a, not for the others. Similarly, if your year starts with peak, you should probably also only add storage_initial for peak.
• For your input and output parameters, please again compare to the test: you probably don't want to add the string "year" for the columns time and time_origin, but the integer 2020 instead. Also, the test is adding two columns equal to the season, (i.e. peak and offpeak in your case) which you don't seem to add. You can use e.g. scenario.idx_names("input") to figure out which index sets your parameters need. If you don't intend to change the values of input and output, please consider using your yearly_to_season() function as described in the tutorial! Finally, please make sure input and output have valid units. You probably don't want to handle dimensionless quantities of electricity.
• Your var_cost also seems to be missing a seasonal index and it's also using "year" for "time", which I still don't understand how it's supposed to work.
• Similarly, your bound_activity_up is also using "year" for "time" but doesn't seem to have a seasonal index as the tests do.
• Lastly, please make sure that you have adapted all parameters in this list to seasonal parameters before solving the model: [p for p in parameters if not scenario.par(p).empty].

Once you have covered all that, the model might already show different results. Please keep in mind that you seem to add a demand equal to 1 (possibly GWa) but limit your power plants to 0.05 GWa of activity, which makes me wonder how the model is able to find a solution at all.
Please further make sure to follow the test setup fully, i.e. to request the variables STORAGE and STORAGE_CHARGE. You will likely not see the storage in the ACT output, but only in those two (or related equations). Please also see the lines below that in the test to see how to inspect those variables.

Hope this helps :)

[shreeyashn20]

Thank you @glatterf42 for your prompt and detailed response. I was referring to the "Central Asia" and "Brazil" models developed by IIASA. Please find the link below. Yes, I have referred to the tutorial you provided. However, the tutorial is limited to seasonality and does not mention how to add storage technologies. Thank you for all the suggestions. I will try to rerun the model incorporating these suggestions. Can you please help me with how to run the test file?
Link 1: https://docs.messageix.org/en/latest/bibliography.html#id13
Link 2: https://github.com/natiweber/BESMM/tree/main

[glatterf42]

Thanks for clarifying the models!
To be clear: my suggestion was to use our test suite for inspiration, not to literally run the tests. I don't think this will provide you much benefit. Also, while most points above are still valid, the most important suggestion (I think) is the last one: to formally request the variables you are interested in.

That being said, here's how you could run the test suite:

Note

This all but requires you to have message_ix installed from source. It definitely requires you to use a terminal to execute some commands.

1. Activate the virtual environment where you have installed message_ix.
2. Make sure this environment contains pytest and all other dependencies required to run the tests.
3. Navigate to the directory from which you installed message_ix.
4. Run pytest message_ix -m "not nightly and not tutorial" -rA --verbose --color=yes --durations=20.

For step 2, please note that there are different ways to do that. If you installed message_ix from source using pip, you can simply navigate to your message_ix-directory and run pip install .[tests] (with your virtual environment active). This will automatically resolve to the list above, which contains further lists like message_ix[docs], which resolves to these packages.
If you're using conda, these features are not supported. As far as I know, your only option would be to go through all these nested and linked lists and ensure manually that all packages mentioned are installed in your virtual environment.

Please let us know if the answer above resolved your initial question :)

[shreeyashn20]

Thank you @glatterf42 ! I was able to use storage feature in MESSAGEix Framework. I am sharing modified test code for storage feature which might be helpful for someone.

import logging
from itertools import product

import pandas as pd

from message_ix import Scenario
import ixmp

mp = ixmp.Platform()

import message_ix

scen = message_ix.Scenario(mp, model="Storage Test", scenario="baseline", version="new")

years=[2020]
scen.add_set("node", "node")
scen.add_set("commodity", "electr")
scen.add_set("level", "level")
scen.add_set("year", years)
scen.add_set("type_year", years)
scen.add_set("technology", ["wind_ppl", "gas_ppl"])
scen.add_set("mode", "M1")
output_specs = ["node", "electr", "level", "year", "year"]
# Chronological order of timesteps in the year
time_order = {"a": 1, "b": 2, "c": 3, "d": 4}
time_duration = {"a": 0.25, "b": 0.25, "c": 0.25, "d": 0.25}

var_cost = {"wind_ppl": 0, "gas_ppl": 2}
for year, (tec, cost) in product(years, var_cost.items()):
    scen.add_par("demand", ["node", "electr", "level", year, "year"], 1, "GWa")
    tec_specs = ["node", tec, year, year, "M1"]
    scen.add_par("output", tec_specs + output_specs, 1, "GWa")
    scen.add_par("var_cost", tec_specs + ["year"], cost, "USD/GWa")

scen.add_set("time", sorted(list(set(time_duration.keys()))))
scen.add_set("lvl_temporal", "season")
for h, duration in time_duration.items():
    scen.add_set("map_temporal_hierarchy", ["season", h, "year"])
    scen.add_par("duration_time", [h], duration, "-")

def year_to_time(scen, parname, time_share):
    old = scen.par(parname)
    scen.remove_par(parname, old)
    time_idx = [x for x in scen.idx_names(parname) if "time" in x]
    for h, share in time_share.items():
        new = old.copy()
        for index in time_idx:
            new[index] = h
        new["value"] = share * old["value"]
        scen.add_par(parname, new)


# A function for adding storage technologies and parameterization
# Adding level of storage
scen.add_set("level", "storage")

# Adding storage technologies (reservoir, charger, and discharger)
scen.add_set("technology", ["dam", "pump", "turbine"])

# Adding a storage commodity
scen.add_set("commodity", ["water", "cooling", "dummay"])

# Specifying storage reservoir technology
scen.add_set("storage_tec", "dam")

# Specifying storage level
scen.add_set("level_storage", "storage")

# Adding mapping for storage and charger/discharger technologies
for tec in ["pump", "turbine"]:
     scen.add_set(
            "map_tec_storage",
            ["node", tec, "M1", "dam", "M1", "storage", "water", "season"],
        )

# Adding time sequence
for h in time_order.keys():
    scen.add_par("time_order", ["season", h], time_order[h], "-")

    # Adding input, output, and capacity factor for storage technologies
output_spec = {
        "dam": ["node", "dummay", "level"],
        "pump": ["node", "water", "storage"],
        "turbine": ["node", "electr", "level"],
    }

input_spec = {
        "dam": ["node", "cooling", "level"],
        "pump": ["node", "electr", "level"],
        "turbine": ["node", "water", "storage"],
    }

var_cost = {"dam": 0, "pump": 0.2, "turbine": 0.3}
stor_tecs = ["dam", "pump", "turbine"]
for year, h, tec in product(set(scen.set("year")), time_order.keys(), stor_tecs):
    tec_sp = ["node", tec, year, year, "M1"]
    scen.add_par("output", tec_sp + output_spec[tec] + [h, h], 1, "GWa")
    scen.add_par("input", tec_sp + input_spec[tec] + [h, h], 1, "GWa")
    scen.add_par("var_cost", tec_sp + [h], var_cost[tec], "USD/GWa")

    # Adding storage self-discharge (as %) and initial content
for year, h in product(set(scen.set("year")), time_order.keys()):
    storage_spec = ["node", "dam", "M1", "storage", "water", year, h]
    scen.add_par("storage_self_discharge", storage_spec, 0.05, "%")

    # Adding initial content of storage (optional)
    storage_spec = ["node", "dam", "M1", "storage", "water", year, "a"]
    scen.add_par("storage_initial", storage_spec, 0.5, "%")


# Main function for building a model with storage and testing the functionality
# First, building a simple model and adding seasonality

# Fixed share for parameters that don't change across timesteps
fixed_share = {"a": 1, "b": 1, "c": 1, "d": 1}
year_to_time(scen, "output", fixed_share)
year_to_time(scen, "var_cost", fixed_share)
    # Variable share for parameters that are changing in each timestep
    # share of demand in each season from annual demand
demand_share = {"a": 0.15, "b": 0.2, "c": 0.4, "d": 0.25}
year_to_time(scen, "demand", demand_share)

for h in time_duration.keys():
        scen.add_par(
            "bound_activity_up", ["node", "wind_ppl", 2020, "M1", h], 0.15, "GWa"
        )
for h in time_duration.keys():
    scen.add_par(
            "bound_activity_up", ["node", "gas_ppl", 2020, "M1", h], 0.15, "GWa"
        )

scen.add_set("technology", "cooler")
df = scen.par("output", {"technology": "turbine"})
df["technology"] = "cooler"
df["commodity"] = "cooling"
scen.add_par("output", df)
    # Changing input of dam from 1 to 1.2 to test commodity balance
df = scen.par("input", {"technology": "dam"})
df["value"] = 1.2
scen.add_par("input", df)

scen.commit("storage data added")
scen.solve(case="with_storage", quiet=True)
act = scen.var("ACT")
act

	node_loc	technology	year_vtg	year_act	mode	time	lvl	mrg
0	node	wind_ppl	2020	2020	M1	a	0.150000	0.000000e+00
1	node	wind_ppl	2020	2020	M1	b	0.150000	0.000000e+00
2	node	wind_ppl	2020	2020	M1	c	0.150000	0.000000e+00
3	node	wind_ppl	2020	2020	M1	d	0.150000	0.000000e+00
4	node	gas_ppl	2020	2020	M1	a	0.005540	0.000000e+00
5	node	gas_ppl	2020	2020	M1	b	0.150000	0.000000e+00
6	node	gas_ppl	2020	2020	M1	c	0.150000	0.000000e+00
7	node	gas_ppl	2020	2020	M1	d	0.100000	0.000000e+00
8	node	dam	2020	2020	M1	a	0.005540	0.000000e+00
9	node	dam	2020	2020	M1	b	0.105263	0.000000e+00
10	node	dam	2020	2020	M1	c	0.000000	0.000000e+00
11	node	dam	2020	2020	M1	d	0.000000	0.000000e+00
12	node	pump	2020	2020	M1	a	0.005540	0.000000e+00
13	node	pump	2020	2020	M1	b	0.100000	0.000000e+00
14	node	pump	2020	2020	M1	c	0.000000	5.000000e-01
15	node	pump	2020	2020	M1	d	0.000000	1.100000e-01
16	node	turbine	2020	2020	M1	a	0.000000	5.000000e-01
17	node	turbine	2020	2020	M1	b	0.000000	5.000000e-01
18	node	turbine	2020	2020	M1	c	0.100000	0.000000e+00
19	node	turbine	2020	2020	M1	d	0.000000	3.900000e-01
20	node	cooler	2020	2020	M1	a	0.006648	0.000000e+00
21	node	cooler	2020	2020	M1	b	0.126316	0.000000e+00
22	node	cooler	2020	2020	M1	c	0.000000	4.940656e-324
23	node	cooler	2020	2020	M1	d	0.000000	4.940656e-324

mp.close_db()