Calibrate and tune MESSAGE-MACRO

“MESSAGE-MACRO” refers to an iterative algorithm that links MESSAGE and MACRO [9]. This algorithm allows to model demand elasticity: MESSAGE solution data on energy prices and total system costs are provided to MACRO, and MACRO solution data on (endogenous) demand and GDP are provided to MESSAGE. This process continues until a convergence criterion or “equilibrium” is reached—briefly, that the MACRO demand output varies minimally between two iterations ([4], further details can be found here).

The linked models can be activated by calling Scenario.solve() with the argument model=’MESSAGE-MACRO’, or using the GAMS MESSAGE-MACRO_run.gms script directly (see Running a model for details about these two methods).

To solve a MESSAGE scenario using MESSAGE-MACRO, it is first necessary to calibrate MACRO. As described in [4], the calibration process…

is parameterized off of a baseline scenario (which assumes some autonomous rate of energy efficiency improvement, AEEI) and is conducted for all MESSAGE regions simultaneously. Therefore, the demand responses motivated by MACRO are meant to represent the additional (compared to the baseline) energy efficiency improvements and conservation that would occur in each region as a result of higher prices for energy services.

In the calibration process, the user provides exogenous, reference energy prices (price_ref) and reference total energy system cost (cost_ref) that correspond to a reference level of demand (demand_ref) in a particular reference year—generally, the ‘historic’ period that directly precedes the first period in the MESSAGE model horizon for optimization (firstmodelyear). This reference year is a period for which commodity prices and energy system cost are known for a given demand of those commodities.

Using these reference values plus energy prices (PRICE_COMMODITY) and total system cost (COST_NODAL_NET) from the solution of MESSAGE for a given demand time series (demand), the calibration process changes two parameters, namely, the autonomous rate of energy efficiency improvement (aeei) and growth in GDP (grow), such that the output of MACRO (GDP and DEMAND) converges to an initially specified time series trajectory of GDP (gdp_calibrate) and demand (demand), respectively. The scenario used for calibration is usually a baseline scenario, meaning that this scenario does not include any constraints that implement policy targets that affect commodity prices (for instance, long-term climate policy targets). Without the calibration, the output of MACRO (GDP and DEMAND) can be different from the initial exogenous assumptions for GDP and demand (gdp_calibrate and demand) in MESSAGE for a given scenario.

The calibration process is invoked using Scenario.add_macro() on a (baseline) scenario and runs for the entire optimization-time horizon, i.e., for all model periods including and after the firstmodelyear. As mentioned, the required price_ref, cost_ref, and demand_ref inputs refer to a period prior to the model horizon. This is detailed in the next section.

The calibration itself is carried out by the message_ix/model/MACRO/macro_calibration.gms. In this iterative process, max_it is used to specify the number of iterations carried out between MESSAGE and MACRO as part of the calibration process. The default value is set to 100 iterations, which has proven to be sufficient for the calibration of MACRO to MESSAGE reference scenario for various models. Adjustment of GDP growth rates (grow) is carried out during even iterations. Adjustment of AEEI improvement rates (aeei) is carried out during odd iterations.

Note

Note, that no actual check is carried out to see if the calibration process has been successful at the end of iterations.

The information from the calibration process is logged in message_ix/model/MACRO_run.lst. Successful calibration of MESSAGE to MACRO can be identified by looking at the reported values for the “PARAMETER growth_correction” for the last “even” iteration, which should be somewhere around 1e-14 to 1e-16 for positive adjustments or -1e-14 to -1e-16 for negative adjustments. Likewise, the “PARAMETER aeei_correction” can be checked for the last “odd” iteration. Once the calibration process has been completed, the scenario will be populated with additional parameters. As part of the calibration process, a final check will automatically be carried out by solving the freshly calibrated scenario in the MESSAGE-MACRO coupled mode, ensuring that the convergence criteria between solution of MESSAGE and MACRO is met after the first iteration.

Input data for calibration 

For calibration, Scenario.add_macro() requires input data that can be provided as either:

a Python dict that maps item names to pandas.DataFrame objects, or
the path to a file in Microsoft Excel format, in which each item is stored as a separate sheet.

For an example of such input data files, see the files message_ix/tests/data/*_macro_input.xlsx included as part of the message_ix test suite; either in your local installation, or here on GitHub.

This section describes the required contents of each item.

`config`: general configuration 

This table/sheet provides structural information for MACRO calibration and the MESSAGE-MACRO linkage. The table has five columns, each of which is a list of labels/codes for a corresponding ixmp set:

“node”, “year”: these columns can each have any length, depending on the number of nodes and periods to be included in the MACRO calibration process.
“sector”, “commodity”, “level”: these 3 columns must have equal lengths. They describe a one-to-one correspondence between items in the MACRO sector set (entries in the “sector” column) and MESSAGE commodity and level sets (paired entries in the “commodity” and “level” columns).

MACRO parameters 

The remaining tables/sheets each contain data for one MACRO parameter. The required dimensions or symbol of each item are given in the same notation used in the documentation of the MACRO core formulation.

price_ref (\(n, s\)): prices of MACRO sector output in the reference year. These can be constructed from the MESSAGE variable PRICE_COMMODITY, using the config mapping. If not provided, message_ix.macro will identify the reference year and extrapolate reference values using an exponential function fitted to PRICE_COMMODITY values; see macro.extrapolate().
cost_ref (\(n\)): total cost of the energy system in the reference year. These can be constructed from the MESSAGE variable COST_NODAL_NET, including dividing by a factor of 1000. If not provided, message_ix.macro will extrapolate using macro.extrapolate().
demand_ref (\(n, s\)): demand for MACRO sector output in the reference year.
lotol (\(n\)): tolerance factor for lower bounds on MACRO variables.
esub (\(\epsilon_n\)): elasticity of substitution between capital-labor and energy.
drate (\(n\)): social discount rate.
depr (\(\mathrm{depr}_n\)): annual percent depreciation.
kpvs (\(\alpha_n\)): capital value share parameter.
kgdp (\(n\)): initial capital to GDP ratio in base year.
gdp_calibrate (\(n, y\)): trajectory of GDP in optimization years calibrated to energy demand to MESSAGE. In order to compute the growth rates in historical years, values are required for the reference year and one prior period—that is, at least two periods prior to the firstmodelyear.
aeei (\(n, s, y\)): annual potential decrease of energy intensity in sector sector.
MERtoPPP (\(n, y\)): conversion factor of GDP from market exchange rates to purchasing power parity.

Numerical issues 

This section describes how to solve two numerical issues that can occur in large MESSAGEix models.

Oscillation detection in the MESSAGE-MACRO algorithm 

The documentation for the MESSAGE_MACRO class describes the algorithm and its three parameters:

convergence_criterion,
max_adjustment, and
max_iteration.

The algorithm detects ‘oscillation’, which occurs when MESSAGE and MACRO each return slightly different solutions, but these two solutions are each stable.

If the difference between these points is greater than convergence_criterion, the algorithm might jump between these two points infinitely. Instead, the algorithm detects oscillation by comparing model solutions on each iteration to previous values recorded in the iteration log. Specifically, the algorithm checks for three patterns across the iterations.

Does the sign of the max_adjustment parameter change?
Are the maximum-positive and maximum-negative adjustments equal to each other?
Do the solutions jump between two objective functions?

If the algorithm picks up on the oscillation between iterations, then after MACRO has solved and before solving MESSAGE, a log message is printed as follows:

--- Restarting execution
--- MESSAGE-MACRO_run.gms(4986) 625 Mb
--- Reading solution for model MESSAGE_MACRO
--- MESSAGE-MACRO_run.gms(4691) 630 Mb
    +++ Indication of oscillation, increase the scaling parameter (4) +++
--- GDX File c:\repo\message_ix\message_ix\model\output\MsgIterationReport_ENGAGE_SSP2_v4_EN_NPi2020_900.gdx
    Time since GAMS start: 1 hour, 10 minutes
    +++ Starting iteration 14 of MESSAGEix-MACRO... +++
    +++ Solve the perfect-foresight version of MESSAGEix +++
--- Generating LP model MESSAGE_LP

Note

This example is from a particular model run, and the actual message may differ.

Which of the three checks listed above has been invoked is logged in the iteration report in MsgIterationReport_<model_name>_<scenario_name>.gdx under the header “oscillation check”.

The algorithm then gradually reduces max_adjustment from the user-supplied value. This has the effect of reducing the allowable relative change in demands, until the convergence_criterion is met.

If none of the checks have been invoked over the iterations, then MESSAGEix and MACRO converged naturally. A log message as follows is printed:

--- Reading solution for model MESSAGE_MACRO
--- Executing after solve: elapsed 7:42:24.622
--- MESSAGE-MACRO_run.gms(5176) 1116 Mb
    +++ Convergence criteria satisfied after 14 iterations +++
    +++ Natural convergence achieved +++

If in any of the iterations, any of the three oscillation checks were invoked, a log message is printed as follows:

--- Reading solution for model MESSAGE_MACRO
--- Executing after solve: elapsed 7:42:24.622
--- MESSAGE-MACRO_run.gms(5176) 1116 Mb
    +++ Convergence criteria satisfied after 14 iterations +++
    +++ Convergence achieved via oscillation check mechanism; check iteration log for further details +++

Issue 1: Oscillations not detected 

Oscillation detection can fail, especially when the oscillation is very small. When this occurs, MESSAGE-MACRO will iterate until max_iteration (default 50) and then print a message indicating that it has not converged.

For the MESSAGEix-GLOBIOM global model, this issue can be encountered with scenarios which have stringent carbon budgets (e.g. <1000 Gt CO₂ cumulative) and require more aggressive reductions of demands.

Identifying oscillation 

In order to find out whether failure to converge is due to undetected oscillation, check the iteration report. The initial iterations will show the objective function value either decreasing or increasing (depending on the model), but after a number of iterations, the objective function will flip-flop between two very similar values.

Preventing oscillation 

The issue can be resolved by tuning max_adjustment and convergence_criterion from their respective default values of 0.2 (20%) and 0.01 (1%). The general approach is to reduce max_adjustment. Reducing this parameter to half of its default value—i.e. 0.1, or 10%—can help, but it can be reduced further, as low as 0.01 (1%).

This may require further tuning of the other parameters: first, ensure that convergence_criterion is smaller than max_adjustment, e.g. set to 0.009 (0.9%) < 0.01. Second, due to the small change allowed to the model solution each iteration, if the initial MESSAGE solution is not close to the convergence point, numerous iterations could be required. Therefore max_iteration may also need an increase.

These changes can be made in two ways:

Pass the values to MESSAGE_MACRO via keyword arguments to Scenario.solve().
Manually edit the default values in MESSAGE-MACRO_run.gms.

Issue 2: MESSAGE solves optimally with unscaled infeasibilities 

By default, message_ix is configured so that the CPLEX solver runs using the lpmethod option set to 4, selecting the barrier method. Solving models the size of MESSAGEix-GLOBIOM would otherwise take very long with the dual simplex method (lpmethod set to 2); scenarios with stringent constraints can take >10 hours on common hardware. With lpmethod set to 4 the model can solve in under a minute.

The drawback of using the barrier method is that, after CPLEX has solved, it crosses over to a simplex optimizer for verification. As part of this verification step, it may turn out that the CPLEX solution is “optimal with unscaled infeasibilities.”

This issue arises when some parameters in the model are not well-scaled, resulting in numerical issues within the solver. This page (from an earlier, 2002 version of the CPLEX user manual) offers some advice on how to overcome the issues. The most direct solution is to rescale the parameters in the model itself.

When this is not possible, there are some workarounds:

Adjust CPLEX’s scaling parameter; specify scaind = 1. This will result in more “aggressive” scaling.
Adjust CPLEX’s barrier crossover algorithm; specify barcrossalg = 2. By default, CPLEX will choose between either Primal crossover or Dual crossover. Unscaled infeasibilities will result only with Primal crossover, hence forcing CPLEX to use the latter will resolve the issue. This will result in longer solving times, but will guarantee overcoming the issue.

Note

This solution has been implemented as part of the MESSAGE-MACRO iterations process. During the iterations, a check is performed on the solution status of MESSAGE. When solving with unscaled infeasibilities, in GAMS, the modelstat will be 1 (Optimal) and the solvestat will be 4 (Terminated by Solver). In this case, a secondary CPLEX configuration file is used for subsequent solving of the MESSAGE model. The secondary CPLEX configuration file message_ixmodelcplex.op2 is a duplicate of message_ixmodelcplex.opt with the addition of the argument barcrossalg = 2. This secondary CPLEX configuration file is generated together with the primary CPLEX configuration file in message_ixmodels.py. Further information on the status description of GAMS can be found here. These differ from those reported by CPLEX.

Adjust CPLEX’s convergence criterion, epopt (this is distinct from the convergence_criterion of the MESSAGE_MACRO algorithm). In message_ix, DEFAULT_CPLEX_OPTIONS sets this to 1e-6 by default. This approach is delicate, as changing the tolerance may also change the solution by a significant amount. This has not been tested in detail and should be handled with care.
Switch to other methods provided by CPLEX, using e.g. lpmethod = 2. A disadvantage of this approach is the longer runtime, as described above.

The arguments can be passed with the solve command, e.g. scenario.solve(solve_options={“barcrossalg”: “2”}) Alternatively the arguments can be specified either in models.py.

Code documentation 

The functions add_model_data() and calibrate() are used by Scenario.add_macro(). Others are internal; prepare_computer() assembles the following functions into a genno.Computer that then executes the necessary calculations to prepare the model data.

`Structures`(level, node, sector, year)	MACRO structure information.
`aconst`(bconst, demand_ref, gdp0, k0, kpvs, rho)	Calculate production function coefficient of capital and labor.
`add_par`(scenario, data, ym1, *, name)	Add data to the scenario.
`add_structure`(scenario, ...)	Add MACRO structure information to scenario.
`bconst`(demand_ref, gdp0, price_ref, rho)	Calculate production function coefficient.
`demand`(model_demand, demand_ref, ...)	Prepare data for the `demand_MESSAGE` MACRO parameter.
`gdp0`(gdp_calibrate, ym1)	Select GDP reference values from "gdp_calibrate".
`growth`(gdp_calibrate)	Calculate GDP growth rates between model periods (MACRO parameter `grow`).
`macro_periods`(demand, config)	Periods ("years") for the MACRO model.
`mapping_macro_sector`(config)	Data for the MACRO set `mapping_macro_sector`.
`price`(model_price, price_ref, ...)	Prepare data for the `price_MESSAGE` MACRO parameter.
`rho`(esub)	Calculate "rho" based on "esub", elasticity of substitution.
`total_cost`(model_cost, cost_ref, ym1)	Combine model_cost and cost_ref (reference year) data.
`unique_set`(column, df)	A `set` of the unique elements in column of df.
`validate_transform`(name, data, s)	Validate df as input data for name, and transform for further calculation.
`ym1`(df, macro_periods)	Period for MACRO initialization: "year minus-one".

The following diagram visualizes the calculation flow:

Diagram of the calculation flow in the calibration of MACRO.

message_ix.macro.INPUT_DATA = ['aeei', 'cost_ref', 'demand_ref', 'depr', 'drate', 'esub', 'gdp_calibrate', 'kgdp', 'kpvs', 'lotol', 'MERtoPPP', 'price_ref'][source]: MACRO calibration parameters to be verified when reading the input data.

class message_ix.macro.Structures(level: set[str], node: set[str], sector: set[str], year: set[int])[source]

MACRO structure information.

year: set[int][source]: Model years for which MACRO is calibrated.

message_ix.macro.aconst(bconst: pandas.Series, demand_ref: pandas.Series, gdp0: pandas.Series, k0: pandas.Series, kpvs: pandas.Series, rho: pandas.Series) → pandas.Series[source]

Calculate production function coefficient of capital and labor.

This is the MACRO GAMS parameter lakl.

message_ix.macro.add_model_data(base: Scenario, clone: Scenario, data: Mapping | PathLike) → None[source]

Calculate and add MACRO structure and data to clone.

Parameters:

base (Scenario) – Base scenario with a solution.
clone (Scenario) – Clone of base scenario for adding calibration parameters.
data (dict) – Data for calibration.

message_ix.macro.add_par(scenario: Scenario, data: pandas.DataFrame, ym1: int, *, name: str) → None[source]: Add data to the scenario.

message_ix.macro.add_structure(scenario: Scenario, mapping_macro_sector: pandas.DataFrame, s: Structures, ym1: int) → None[source]: Add MACRO structure information to scenario.

message_ix.macro.bconst(demand_ref: pandas.DataFrame, gdp0: pandas.Series, price_ref: pandas.DataFrame, rho: pandas.Series) → pandas.DataFrame[source]

Calculate production function coefficient.

This is the MACRO GAMS parameter prfconst.

message_ix.macro.calibrate(s, check_convergence: bool = True, **kwargs)[source]

Calibrate a MESSAGEix scenario to parameters of MACRO.

Parameters:

s (Scenario) – MESSAGEix scenario with calibration data.
check_convergence (bool, optional) – Test is MACRO-calibrated scenario converges in one iteration.
**kwargs – Keyword arguments passed to meth:message_ix.Scenario.solve.

Raises:

RuntimeError – If calibrated scenario solves in more than one iteration.

Returns:

s – MACRO-calibrated scenario.

Return type:

message_ix.Scenario()

message_ix.macro.clean_model_data(data: genno.Quantity, s: Structures) → pandas.DataFrame[source]

Clean MESSAGE variable data for calibration of MACRO parameters.

Parameters:: data (genno.Quantity) – With short dimension names (c, l, n, y), etc.
Returns:: With full column names and a “value” column. Only the labels in s (levels, nodes, sectors, and years) appear in the respective dimensions.
Return type:: pandas.DataFrame

message_ix.macro.demand(model_demand: pandas.DataFrame, demand_ref: pandas.DataFrame, mapping_macro_sector: pandas.DataFrame, ym1: int) → pandas.DataFrame[source]

Prepare data for the demand_MESSAGE MACRO parameter.

Parameters:

model_demand – Values from the DEMAND MESSAGE variable.
demand_ref – Reference values to use for the period ym1.
mapping_macro_sector – MACRO set of the same name; see mapping_macro_sector().
ym1 – First pre-model period; see ym1().

Returns:

With the dimensions (node, sector, year).

Return type:

pandas.DataFrame

Raises:

RuntimeError – If there zero or missing values in the computed data.

message_ix.macro.extrapolate(model_data: pandas.DataFrame, mapping_macro_sector: pandas.DataFrame, ym1: int) → pandas.Series[source]

Extrapolate model_data to cover period ym1.

The extrapolation is done by fitting \(y = b \times m ^ x\), i.e. with two parameters b and m. This is identical to the GROWTH() function in Microsoft Excel ( https://support.microsoft.com/en-us/office/growth-function-541a91dc-3d5e-437d-b156-21324e68b80d ). Data are grouped on all other dimensions, and fitting/extrapolation is performed for each group.

Returns:: The index does not have a year dimension; the data are implicitly for ym1.
Return type:: pandas.Series

message_ix.macro.gdp0(gdp_calibrate, ym1: int) → pandas.Series[source]: Select GDP reference values from “gdp_calibrate”.

message_ix.macro.growth(gdp_calibrate) → pandas.DataFrame[source]: Calculate GDP growth rates between model periods (MACRO parameter grow).

message_ix.macro.macro_periods(demand: genno.Quantity, config: pandas.DataFrame) → set[int][source]

Periods (“years”) for the MACRO model.

The intersection of those appearing in the config data and in the DEMAND variable of the base scenario.

message_ix.macro.mapping_macro_sector(config: pandas.DataFrame) → pandas.DataFrame[source]: Data for the MACRO set mapping_macro_sector.

message_ix.macro.prepare_computer(base: Scenario, target: Scenario | None = None, data: Mapping | PathLike | None = None) → genno.Computer[source]

Prepare a Reporter to perform MACRO calibration calculations.

Parameters:

base (message_ix.Scenario) – Must have a stored solution.
target (message_ix.Scenario) – Scenario to which to add computed data.
data (dict or os.PathLike) – If PathLike, the path to an Excel file containing parameter data, one per sheet. If dict, a mapping from parameter names to data frames.

Raises:

ValueError – if any of the require parameters for MACRO calibration (INPUT_DATA) is missing.