From d638b61781d46aa54b663a9f32f3013a95428283 Mon Sep 17 00:00:00 2001 From: Joe DeCarolis Date: Fri, 11 Sep 2026 14:11:01 -0400 Subject: [PATCH 01/10] Remove stale mypy type-ignore for pint UnitRegistry pint 0.26.1 changed UnitRegistry's type stub so the type-arg ignore is no longer needed; the unused ignore itself now fails mypy under warn_unused_ignores. Unblocks the dependency canary. Ref #372 --- temoa/model_checking/unit_checking/__init__.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/temoa/model_checking/unit_checking/__init__.py b/temoa/model_checking/unit_checking/__init__.py index 0db3c443f..4cfcd2c1c 100644 --- a/temoa/model_checking/unit_checking/__init__.py +++ b/temoa/model_checking/unit_checking/__init__.py @@ -3,8 +3,7 @@ from pint import UnitRegistry from pint.errors import DefinitionSyntaxError -# UnitRegistry is generic but doesn't require type args at instantiation -ureg: UnitRegistry = UnitRegistry() # type: ignore[type-arg] +ureg: UnitRegistry = UnitRegistry() # Load custom unit definitions from the package resources _resource_path = 'temoa.model_checking.unit_checking/temoa_units.txt' From fc10702d3f6528cb2de731c03d1202f7bb5980d7 Mon Sep 17 00:00:00 2001 From: Joe DeCarolis Date: Fri, 11 Sep 2026 14:20:11 -0400 Subject: [PATCH 02/10] chore: bump pint to 0.26.1 to match unit_checking type fix The prior commit (d638b617) removed the type: ignore on UnitRegistry() assuming pint 0.26.1's updated stub, but the committed lockfile still pinned 0.25.3, breaking CI's mypy check in the other direction. Bump pint so lockfile and source agree. Ref #372 --- requirements-dev.txt | 6 +++--- requirements.txt | 6 +++--- uv.lock | 6 +++--- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/requirements-dev.txt b/requirements-dev.txt index 49a81cce6..33c98884a 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1148,9 +1148,9 @@ pillow==12.3.0 \ --hash=sha256:fe3cca2e4e8a592be0f269a1ca4835c25199d9f3ce815c8491048f785b0a0198 \ --hash=sha256:ffd0c5368496f41b0944be820fcb7a838aa6e623d250b01acf2643939c3f99d7 # via matplotlib -pint==0.25.3 \ - --hash=sha256:27eb25143bd5de9fcc4d5a4b484f16faf6b4615aa93ece6b3373a8c1a3c1b97d \ - --hash=sha256:f8f5df6cf65314d74da1ade1bf96f8e3e4d0c41b51577ac53c49e7d44ca5acee +pint==0.26.1 \ + --hash=sha256:1bbde36eae57a5a289cd05081c6405618a5899814940752064ce351cd0204f71 \ + --hash=sha256:e982b129415c09c63308f314ae44697d83e5c96f253bcc0d5b833fc3329eb6d4 # via temoa platformdirs==4.11.7 \ --hash=sha256:4f41487eeeeeb07f3a6625e61d9bc0ae6809f92d3386dbd74392fbb76108104d \ diff --git a/requirements.txt b/requirements.txt index f4af38140..3da57ca78 100644 --- a/requirements.txt +++ b/requirements.txt @@ -557,9 +557,9 @@ pillow==12.3.0 \ --hash=sha256:fe3cca2e4e8a592be0f269a1ca4835c25199d9f3ce815c8491048f785b0a0198 \ --hash=sha256:ffd0c5368496f41b0944be820fcb7a838aa6e623d250b01acf2643939c3f99d7 # via matplotlib -pint==0.25.3 \ - --hash=sha256:27eb25143bd5de9fcc4d5a4b484f16faf6b4615aa93ece6b3373a8c1a3c1b97d \ - --hash=sha256:f8f5df6cf65314d74da1ade1bf96f8e3e4d0c41b51577ac53c49e7d44ca5acee +pint==0.26.1 \ + --hash=sha256:1bbde36eae57a5a289cd05081c6405618a5899814940752064ce351cd0204f71 \ + --hash=sha256:e982b129415c09c63308f314ae44697d83e5c96f253bcc0d5b833fc3329eb6d4 # via temoa platformdirs==4.11.7 \ --hash=sha256:4f41487eeeeeb07f3a6625e61d9bc0ae6809f92d3386dbd74392fbb76108104d \ diff --git a/uv.lock b/uv.lock index f1913064d..dbf2fca40 100644 --- a/uv.lock +++ b/uv.lock @@ -1474,7 +1474,7 @@ wheels = [ [[package]] name = "pint" -version = "0.25.3" +version = "0.26.1" source = { registry = "https://pypi.org/simple" } dependencies = [ { name = "flexcache" }, @@ -1482,9 +1482,9 @@ dependencies = [ { name = "platformdirs" }, { name = "typing-extensions" }, ] -sdist = { url = "https://files.pythonhosted.org/packages/52/9d/b1379cdbd33a49d17d627bc24e2b63cca06a1c5343b38072d2889499e82e/pint-0.25.3.tar.gz", hash = "sha256:f8f5df6cf65314d74da1ade1bf96f8e3e4d0c41b51577ac53c49e7d44ca5acee", size = 255106, upload-time = "2026-03-19T21:57:08.72Z" } +sdist = { url = "https://files.pythonhosted.org/packages/9f/bc/2c38c32e0fb1f966d3695f4493a3f3ce2cc0cca1bbe7c958b92261b31af9/pint-0.26.1.tar.gz", hash = "sha256:1bbde36eae57a5a289cd05081c6405618a5899814940752064ce351cd0204f71", size = 273631, upload-time = "2026-09-10T21:16:49.712Z" } wheels = [ - { url = "https://files.pythonhosted.org/packages/1b/dd/a9fe6a0a09512da23951c68bf36466aeecd89def3183dc095edbc807ddc5/pint-0.25.3-py3-none-any.whl", hash = "sha256:27eb25143bd5de9fcc4d5a4b484f16faf6b4615aa93ece6b3373a8c1a3c1b97d", size = 307488, upload-time = "2026-03-19T21:57:07.022Z" }, + { url = "https://files.pythonhosted.org/packages/6d/5c/9507ea3c732f8a259f45ddf190a741e9bf7d3fe229b5d473917ee6efbf1e/pint-0.26.1-py3-none-any.whl", hash = "sha256:e982b129415c09c63308f314ae44697d83e5c96f253bcc0d5b833fc3329eb6d4", size = 326478, upload-time = "2026-09-10T21:16:48.316Z" }, ] [[package]] From f4310e6e46d8e1905593c91b006a53774395fb87 Mon Sep 17 00:00:00 2001 From: Yamil Essus Date: Fri, 25 Sep 2026 13:27:46 -0400 Subject: [PATCH 03/10] Added solver configuration interface --- README.md | 4 +- temoa/_internal/run_actions.py | 34 +++---- temoa/_internal/temoa_sequencer.py | 2 + temoa/core/config.py | 82 +++++++++++++---- temoa/core/solver_spec.py | 90 +++++++++++++++++++ temoa/extensions/method_of_morris/morris.py | 7 +- .../method_of_morris/morris_evaluate.py | 7 +- .../method_of_morris/morris_sequencer.py | 13 ++- .../MGA_solver_options.toml | 10 +++ .../mga_sequencer.py | 12 ++- .../monte_carlo/MC_solver_options.toml | 10 +++ temoa/extensions/monte_carlo/mc_sequencer.py | 5 +- temoa/extensions/myopic/myopic_sequencer.py | 6 +- .../single_vector_mga/sv_mga_sequencer.py | 3 + .../stochastics/stochastic_sequencer.py | 10 ++- temoa/tutorial_assets/config_sample.toml | 9 +- 16 files changed, 246 insertions(+), 58 deletions(-) create mode 100644 temoa/core/solver_spec.py diff --git a/README.md b/README.md index 55196dbe4..ea3747c00 100644 --- a/README.md +++ b/README.md @@ -178,7 +178,7 @@ config = TemoaConfig( time_sequencing="seasonal_timeslices", input_database="tutorial_database.sqlite", output_database="tutorial_database.sqlite", - solver_name="appsi_highs", + solver="appsi_highs", output_path=output_path, silent=False, ) @@ -238,7 +238,7 @@ scenario = "tutorial" scenario_mode = "perfect_foresight" input_database = "tutorial_database.sqlite" output_database = "tutorial_database.sqlite" -solver_name = "appsi_highs" +solver = "appsi_highs" ``` ### Configuration Options diff --git a/temoa/_internal/run_actions.py b/temoa/_internal/run_actions.py index dcdc01bd5..d586dfa9a 100644 --- a/temoa/_internal/run_actions.py +++ b/temoa/_internal/run_actions.py @@ -3,12 +3,13 @@ """ import sqlite3 -from collections.abc import Generator, Iterable +from collections.abc import Generator, Iterable, Mapping from contextlib import contextmanager from logging import getLogger from pathlib import Path from sys import version_info from time import perf_counter +from typing import Any from pyomo.environ import ( Constraint, @@ -25,6 +26,7 @@ from temoa._internal.table_writer import TableWriter from temoa.core.config import TemoaConfig from temoa.core.model import TemoaModel +from temoa.core.solver_spec import DEFAULT_SOLVER_OPTIONS from temoa.data_processing.db_to_excel import make_excel logger = getLogger(__name__) @@ -174,6 +176,7 @@ def solve_instance( solver_name: str, silent: bool = False, solver_suffixes: Iterable[str] | None = None, + solver_options: Mapping[str, Any] | None = None, ) -> tuple[TemoaModel, SolverResults]: """ Solve the instance and return a loaded instance @@ -181,6 +184,8 @@ def solve_instance( 'duals' is supported in the Temoa Framework. Some solvers may not support duals. :param silent: Run silently :param solver_name: The name of the solver to request from the SolverFactory + :param solver_options: options to set on the solver (see resolve_solver_options). If None, + the Temoa defaults for the solver are used :param instance: the instance to solve :return: loaded instance """ @@ -202,28 +207,11 @@ def solve_instance( if solver_name == 'neos': raise NotImplementedError('Neos based solve is not currently supported') - # Solver Configuration - if solver_name == 'cbc': - pass - - elif solver_name == 'cplex': - # Note: these parameter values match mip-dev / PyPSA - # (see: https://pypsa-eur.readthedocs.io/en/latest/configuration.html) - optimizer.options['lpmethod'] = 4 # barrier - optimizer.options['solutiontype'] = 2 # non basic solution, ie no crossover - optimizer.options['barrier convergetol'] = 1.0e-3 - optimizer.options['feasopt tolerance'] = 1.0e-4 - - elif solver_name == 'gurobi': - # Note: these parameter values match mip-dev / PyPSA (see: https://pypsa-eur.readthedocs.io/en/latest/configuration.html) - optimizer.options['Method'] = 2 # barrier - optimizer.options['Crossover'] = 0 # non basic solution, ie no crossover - optimizer.options['BarConvTol'] = 1.0e-3 - optimizer.options['FeasibilityTol'] = 1.0e-4 - optimizer.options['BarOrder'] = -1 # auto ordering; 2-4x faster than AMD on large models - - elif solver_name == 'appsi_highs': - pass + # Solver Configuration (defaults live in temoa.core.solver_spec.DEFAULT_SOLVER_OPTIONS) + if solver_options is None: + solver_options = DEFAULT_SOLVER_OPTIONS.get(solver_name, {}) + for option, option_value in solver_options.items(): + optimizer.options[option] = option_value # Suffix Handling solver_suffixes_list: list[str] = [] diff --git a/temoa/_internal/temoa_sequencer.py b/temoa/_internal/temoa_sequencer.py index bd44044e3..d781aad1b 100644 --- a/temoa/_internal/temoa_sequencer.py +++ b/temoa/_internal/temoa_sequencer.py @@ -27,6 +27,7 @@ from temoa.core.config import TemoaConfig from temoa.core.model import TemoaModel from temoa.core.modes import TemoaMode +from temoa.core.solver_spec import resolve_solver_options from temoa.data_io.hybrid_loader import HybridLoader from temoa.extensions.method_of_morris.morris_sequencer import MorrisSequencer from temoa.extensions.modeling_to_generate_alternatives.mga_sequencer import MgaSequencer @@ -254,6 +255,7 @@ def _run_perfect_foresight(self) -> None: self.config.solver_name, silent=self.config.silent, solver_suffixes=suffixes, + solver_options=resolve_solver_options(self.config.solver), ) good_solve, msg = check_solve_status(self.pf_results) if not good_solve: diff --git a/temoa/core/config.py b/temoa/core/config.py index 3e68c9225..26a095618 100644 --- a/temoa/core/config.py +++ b/temoa/core/config.py @@ -1,10 +1,14 @@ import shutil import sys import tomllib +import warnings +from collections.abc import Mapping from logging import getLogger from pathlib import Path +from typing import Any from temoa.core.modes import TemoaMode +from temoa.core.solver_spec import SolverSpec from temoa.extensions.framework import normalize_extension_ids, resolve_extension_specs logger = getLogger(__name__) @@ -42,7 +46,7 @@ def __init__( input_database: Path, output_database: Path, output_path: Path, - solver_name: str, + solver_name: str | None = None, neos: bool = False, save_excel: bool = False, save_duals: bool = False, @@ -73,6 +77,7 @@ def __init__( output_threshold_cost: float | None = None, sqlite: dict[str, object] | None = None, extensions: list[str] | tuple[str, ...] | None = None, + solver: str | Mapping[str, Any] | SolverSpec | None = None, ): if '-' in scenario: raise ValueError( @@ -129,7 +134,20 @@ def __init__( self.neos = neos if self.neos: raise NotImplementedError('Neos is currently not supported.') - self.solver_name = solver_name + + # Validate solver input + if solver_name is not None: + if solver is not None: + raise ValueError("Specify either 'solver' or 'solver_name', not both") + warnings.warn( + "The 'solver_name' argument is deprecated, use 'solver' instead", + DeprecationWarning, + stacklevel=2, + ) + solver = solver_name + if solver is None: + raise SolverNotAvailableError('No solver specified in the configuration.') + self.solver = SolverSpec.parse(solver) self.save_excel = save_excel self.save_duals = save_duals @@ -230,6 +248,17 @@ def __init__( if not self.silent: sys.stderr.write('Warning: ' + msg) + @property + def solver_name(self) -> str: + """The name of the selected solver (shorthand for self.solver.name)""" + return self.solver.name + + @solver_name.setter + def solver_name(self, value: str) -> None: + # retained for backward compatibility. Changing solvers drops any configured options, + # as they are specific to the previous solver + self.solver = SolverSpec.parse(value) + @staticmethod def _check_solver_availability(solver_name: str) -> tuple[bool, str | None]: """ @@ -282,27 +311,41 @@ def build_config(config_file: Path, output_path: Path, silent: bool = False) -> data = tomllib.load(f) if 'solver_name' in data: - is_available, location = TemoaConfig._check_solver_availability(data['solver_name']) - if not is_available: - error_message = ( - f"The specified solver '{data['solver_name']}' was not found.\n" - 'Please ensure the solver is installed and accessible.\n' + if 'solver' in data: + raise ValueError( + "Config specifies both 'solver' and 'solver_name'. Use only 'solver' " + "('solver_name' is deprecated)." ) - if data['solver_name'].lower() in SOLVER_DOC_LINKS: - link = SOLVER_DOC_LINKS[data['solver_name'].lower()] - error_message += f'For installation instructions, refer to: {link}\n' - else: - error_message += ( - "Refer to the solver's official documentation for " - 'installation instructions.' - ) - raise SolverNotAvailableError(error_message) + logger.warning( + "The 'solver_name' config key is deprecated and will be removed in a future " + 'release. Replace it with: solver = "%s"', + data['solver_name'], + ) + data['solver'] = data.pop('solver_name') + + if 'solver' not in data: + raise SolverNotAvailableError('No solver name specified in the configuration.') + solver = SolverSpec.parse(data['solver']) + data['solver'] = solver + + is_available, location = TemoaConfig._check_solver_availability(solver.name) + if not is_available: + error_message = ( + f"The specified solver '{solver.name}' was not found.\n" + 'Please ensure the solver is installed and accessible.\n' + ) + if solver.name.lower() in SOLVER_DOC_LINKS: + link = SOLVER_DOC_LINKS[solver.name.lower()] + error_message += f'For installation instructions, refer to: {link}\n' else: - logger.info('Using solver: %s (%s)', data['solver_name'], location) + error_message += ( + "Refer to the solver's official documentation for installation instructions." + ) + raise SolverNotAvailableError(error_message) else: - raise SolverNotAvailableError('No solver name specified in the configuration.') + logger.info('Using solver: %s (%s)', solver.name, location) - if data.get('solver_name') == 'appsi_highs' and data.get('save_duals', False): + if solver.name == 'appsi_highs' and data.get('save_duals', False): raise ValueError( 'save_duals is not supported with appsi_highs (it does not expose duals via the ' 'APPSI interface). Disable save_duals or choose a different solver.' @@ -354,6 +397,7 @@ def __repr__(self) -> str: msg += spacer msg += '{:>{}s}: {}\n'.format('Selected solver', width, self.solver_name) + msg += '{:>{}s}: {}\n'.format('Solver options', width, dict(self.solver.options)) msg += '{:>{}s}: {}\n'.format('NEOS status', width, self.neos) msg += spacer diff --git a/temoa/core/solver_spec.py b/temoa/core/solver_spec.py new file mode 100644 index 000000000..45c6ca1f3 --- /dev/null +++ b/temoa/core/solver_spec.py @@ -0,0 +1,90 @@ +""" +Solver selection and option resolution. + +The config accepts ``solver`` as either a plain solver name or a table with a ``name`` and a +passthrough ``options`` table. Options that reach the solver are layered as: + + DEFAULT_SOLVER_OPTIONS < [solver.options] < extension-specific options (MGA, MC, ...) + +Any solver name known to pyomo's SolverFactory is accepted. Solvers without an entry in +DEFAULT_SOLVER_OPTIONS simply get no Temoa defaults. +""" + +from __future__ import annotations + +from collections.abc import Mapping +from dataclasses import dataclass, field +from typing import Any + +# Note: these parameter values match mip-dev / PyPSA +# (see: https://pypsa-eur.readthedocs.io/en/latest/configuration.html) +DEFAULT_SOLVER_OPTIONS: dict[str, dict[str, Any]] = { + 'cplex': { + 'lpmethod': 4, # barrier + 'solutiontype': 2, # non basic solution, ie no crossover + 'barrier convergetol': 1.0e-3, + 'feasopt tolerance': 1.0e-4, + }, + 'gurobi': { + 'Method': 2, # barrier + 'Crossover': 0, # non basic solution, ie no crossover + 'BarConvTol': 1.0e-3, + 'FeasibilityTol': 1.0e-4, + 'BarOrder': -1, # auto ordering; 2-4x faster than AMD on large models + }, +} + + +@dataclass(frozen=True, slots=True) +class SolverSpec: + """A solver name plus the user-supplied options from the config (defaults excluded).""" + + name: str + options: Mapping[str, Any] = field(default_factory=dict) + + @classmethod + def parse(cls, raw: str | Mapping[str, Any] | SolverSpec) -> SolverSpec: + """ + Build a SolverSpec from a solver name, a {'name': ..., 'options': {...}} mapping, or an + existing SolverSpec + """ + if isinstance(raw, SolverSpec): + return raw + if isinstance(raw, str): + if not raw: + raise ValueError('Solver name must not be empty') + return cls(raw) + if isinstance(raw, Mapping): + unknown = set(raw) - {'name', 'options'} + if unknown: + raise ValueError( + f'Unrecognized key(s) in solver table: {sorted(unknown)}. Expected "name" ' + 'and (optionally) "options". Solver parameters belong under [solver.options]' + ) + name = raw.get('name') + if not isinstance(name, str) or not name: + raise ValueError('The solver table requires a non-empty "name" entry') + options = raw.get('options', {}) + if not isinstance(options, Mapping): + raise ValueError('The "options" entry of the solver table must be a table/dict') + return cls(name, dict(options)) + raise TypeError(f'solver must be a str or a table/dict, got: {type(raw).__name__}') + + +def resolve_solver_options( + spec: SolverSpec, + extension_options: Mapping[str, Any] | None = None, + *, + include_defaults: bool = True, +) -> dict[str, Any]: + """ + Merge the option layers that are passed to the solver + :param spec: the solver spec from the config + :param extension_options: options from an extension's own source (MGA/MC/Morris toml, + stochastic config), which take precedence over everything else + :param include_defaults: include DEFAULT_SOLVER_OPTIONS as the bottom layer. Solve paths that + historically ran without Temoa defaults (MGA base solve, stochastic) turn this off. + :return: a new dict of solver options + """ + defaults = DEFAULT_SOLVER_OPTIONS.get(spec.name, {}) if include_defaults else {} + return {**defaults, **spec.options, **(extension_options or {})} diff --git a/temoa/extensions/method_of_morris/morris.py b/temoa/extensions/method_of_morris/morris.py index 4ae1f308f..19dc114e9 100644 --- a/temoa/extensions/method_of_morris/morris.py +++ b/temoa/extensions/method_of_morris/morris.py @@ -16,6 +16,7 @@ from temoa._internal import run_actions from temoa._internal.table_writer import TableWriter from temoa.core.config import TemoaConfig +from temoa.core.solver_spec import resolve_solver_options from temoa.data_io.hybrid_loader import HybridLoader seed = 42 @@ -40,7 +41,11 @@ def evaluate( dp = DataPortal(data_dict={None: data}) instance = run_actions.build_instance(loaded_portal=dp, extensions=config.extensions) - mdl, res = run_actions.solve_instance(instance=instance, solver_name=config.solver_name) + mdl, res = run_actions.solve_instance( + instance=instance, + solver_name=config.solver_name, + solver_options=resolve_solver_options(config.solver), + ) status = run_actions.check_solve_status(res) if not status: raise RuntimeError('Bad solve during Method of Morris') diff --git a/temoa/extensions/method_of_morris/morris_evaluate.py b/temoa/extensions/method_of_morris/morris_evaluate.py index 9695b1406..daa3e8064 100644 --- a/temoa/extensions/method_of_morris/morris_evaluate.py +++ b/temoa/extensions/method_of_morris/morris_evaluate.py @@ -43,6 +43,7 @@ def evaluate( config: TemoaConfig, log_queue: Any, log_level: int, + solver_options: dict[str, Any] | None = None, ) -> list[float]: """ Run model for params provided and return objective value and emission value @@ -54,6 +55,7 @@ def evaluate( :param data: Data used to build the Data Portal :param i: indexing number :param config: The config file to pull run data from + :param solver_options: resolved solver options. If None, the Temoa defaults are used :return: list of objective value and CO2 emission value """ # get the logger configured... @@ -80,7 +82,10 @@ def evaluate( extensions=config.extensions, ) mdl, res = run_actions.solve_instance( - instance=instance, solver_name=config.solver_name, silent=True + instance=instance, + solver_name=config.solver_name, + silent=True, + solver_options=solver_options, ) status = run_actions.check_solve_status(res) if not status: diff --git a/temoa/extensions/method_of_morris/morris_sequencer.py b/temoa/extensions/method_of_morris/morris_sequencer.py index 6c3d1bc76..678417439 100644 --- a/temoa/extensions/method_of_morris/morris_sequencer.py +++ b/temoa/extensions/method_of_morris/morris_sequencer.py @@ -22,6 +22,7 @@ from SALib.util import compute_groups_matrix, read_param_file # type: ignore[import-untyped] from temoa._internal.table_writer import TableWriter +from temoa.core.solver_spec import resolve_solver_options from temoa.data_io.hybrid_loader import HybridLoader from temoa.extensions.method_of_morris.morris_evaluate import evaluate @@ -70,11 +71,12 @@ def __init__(self, config: TemoaConfig): with open(path, 'rb') as f: all_options = tomllib.load(f) s_options = all_options.get(self.config.solver_name, {}) - logger.info('Using solver options: %s', s_options) except FileNotFoundError: logger.warning('Unable to find solver options toml file. Using default options.') s_options = {} + self.solver_options = resolve_solver_options(self.config.solver, s_options) + logger.info('Using solver options: %s', self.solver_options) # output handling self.verbose = False # for troubleshooting @@ -189,7 +191,14 @@ def start(self) -> Any: sys.stdout.flush() morris_results = Parallel(n_jobs=self.num_cores)( delayed(evaluate)( - param_names, mm_samples[i, :], data, i, self.config, log_queue, log_level + param_names, + mm_samples[i, :], + data, + i, + self.config, + log_queue, + log_level, + solver_options=self.solver_options, ) for i in range(0, len(mm_samples)) ) diff --git a/temoa/extensions/modeling_to_generate_alternatives/MGA_solver_options.toml b/temoa/extensions/modeling_to_generate_alternatives/MGA_solver_options.toml index 869a53384..2b57d00c6 100644 --- a/temoa/extensions/modeling_to_generate_alternatives/MGA_solver_options.toml +++ b/temoa/extensions/modeling_to_generate_alternatives/MGA_solver_options.toml @@ -1,5 +1,7 @@ # A container for solver options # the top level solver name in brackets should align with the solver name in the config.toml +# These are layered on top of Temoa's solver defaults and the config's [solver.options], and take +# precedence over both. (see temoa/core/solver_spec.py) num_workers = 6 @@ -11,6 +13,7 @@ BarConvTol = 0.01 # Relative Barrier Tolerance primal-dual FeasibilityTol= 1e-2 # pretty loose Crossover= 0 # Disabled TimeLimit= 18000 # 5 hrs +BarOrder = -1 # auto ordering (gurobi default) # regarding BarConvTol: https://www.gurobi.com/documentation/current/refman/barrier_logging.html # note that ref above seems to imply that FeasibilyTol is NOT used when using barrier only...? @@ -19,6 +22,13 @@ TimeLimit= 18000 # 5 hrs # 'LogFile': './my_gurobi_log.log', # 'LPWarmStart': 2, # pass basis +[cplex] +# CPLEX's own defaults, restating them overrides the Temoa cplex defaults for MGA workers +lpmethod = 0 # automatic +solutiontype = 0 # automatic +'barrier convergetol' = 1.0e-8 +'feasopt tolerance' = 1.0e-6 + [cbc] # tbd diff --git a/temoa/extensions/modeling_to_generate_alternatives/mga_sequencer.py b/temoa/extensions/modeling_to_generate_alternatives/mga_sequencer.py index a261ba458..6f3f124d3 100644 --- a/temoa/extensions/modeling_to_generate_alternatives/mga_sequencer.py +++ b/temoa/extensions/modeling_to_generate_alternatives/mga_sequencer.py @@ -33,6 +33,7 @@ from temoa._internal.run_actions import build_instance from temoa._internal.table_writer import TableWriter from temoa.components.costs import total_cost_rule +from temoa.core.solver_spec import resolve_solver_options from temoa.data_io.hybrid_loader import HybridLoader from temoa.extensions.modeling_to_generate_alternatives.manager_factory import get_manager from temoa.extensions.modeling_to_generate_alternatives.mga_constants import MgaAxis, MgaWeighting @@ -93,16 +94,21 @@ def __init__(self, config: TemoaConfig): with open(path, 'rb') as f: all_options = tomllib.load(f) s_options = all_options.get(self.config.solver_name, {}) - logger.info('Using solver options: %s', s_options) except FileNotFoundError: logger.warning('Unable to find solver options toml file. Using default options.') s_options = {} all_options = {} - # get handle on solver instance + # get handle on solver instance. The base solve only receives the user's + # [solver.options] (no Temoa defaults, no worker options) to get a more precise base cost self.opt = pyo.SolverFactory(self.config.solver_name) - self.worker_solver_options = s_options + for option, option_value in resolve_solver_options( + self.config.solver, include_defaults=False + ).items(): + self.opt.options[option] = option_value + self.worker_solver_options = resolve_solver_options(self.config.solver, s_options) + logger.info('Using worker solver options: %s', self.worker_solver_options) # some defaults, etc. self.internal_stop = False diff --git a/temoa/extensions/monte_carlo/MC_solver_options.toml b/temoa/extensions/monte_carlo/MC_solver_options.toml index da91cd56e..7a7471d61 100644 --- a/temoa/extensions/monte_carlo/MC_solver_options.toml +++ b/temoa/extensions/monte_carlo/MC_solver_options.toml @@ -1,5 +1,7 @@ # A container for solver options # the top level solver name in brackets should align with the solver name in the config.toml +# These are layered on top of Temoa's solver defaults and the config's [solver.options], and take +# precedence over both. (see temoa/core/solver_spec.py) num_workers = 11 @@ -11,6 +13,7 @@ BarConvTol = 1.0e-2 # Relative Barrier Tolerance primal-dual FeasibilityTol= 1.0e-2 # pretty loose Crossover= 0 # Disabled TimeLimit= 18000 # 5 hrs +BarOrder = -1 # auto ordering (gurobi default) # regarding BarConvTol: https://www.gurobi.com/documentation/current/refman/barrier_logging.html # note that ref above seems to imply that FeasibilyTol is NOT used when using barrier only...? @@ -19,6 +22,13 @@ TimeLimit= 18000 # 5 hrs # 'LogFile': './my_gurobi_log.log', # 'LPWarmStart': 2, # pass basis +[cplex] +# CPLEX's own defaults, restating them overrides the Temoa cplex defaults for MC workers +lpmethod = 0 # automatic +solutiontype = 0 # automatic +'barrier convergetol' = 1.0e-8 +'feasopt tolerance' = 1.0e-6 + [cbc] primalT = 1e-3 dualT = 1e-3 diff --git a/temoa/extensions/monte_carlo/mc_sequencer.py b/temoa/extensions/monte_carlo/mc_sequencer.py index eb99da89a..840df6030 100644 --- a/temoa/extensions/monte_carlo/mc_sequencer.py +++ b/temoa/extensions/monte_carlo/mc_sequencer.py @@ -18,6 +18,7 @@ from typing import TYPE_CHECKING, Any, cast from temoa._internal.table_writer import TableWriter +from temoa.core.solver_spec import resolve_solver_options from temoa.data_io.hybrid_loader import HybridLoader from temoa.extensions.monte_carlo.mc_run import MCRun, MCRunFactory from temoa.extensions.monte_carlo.mc_worker import MCWorker @@ -83,7 +84,6 @@ def __init__(self, config: TemoaConfig): with open(path, 'rb') as f: all_options = tomllib.load(f) s_options = all_options.get(self.config.solver_name, {}) - logger.info('Using solver options: %s', s_options) except FileNotFoundError: if options_file_path: @@ -95,7 +95,8 @@ def __init__(self, config: TemoaConfig): # worker options pulled from file self.num_workers = all_options.get('num_workers', 1) - self.worker_solver_options = s_options + self.worker_solver_options = resolve_solver_options(self.config.solver, s_options) + logger.info('Using solver options: %s', self.worker_solver_options) # internal records self.solve_count = 0 diff --git a/temoa/extensions/myopic/myopic_sequencer.py b/temoa/extensions/myopic/myopic_sequencer.py index e32ed347f..8d01c0c04 100644 --- a/temoa/extensions/myopic/myopic_sequencer.py +++ b/temoa/extensions/myopic/myopic_sequencer.py @@ -16,6 +16,7 @@ from temoa._internal.table_writer import TableWriter from temoa.core.config import TemoaConfig from temoa.core.model import TemoaModel +from temoa.core.solver_spec import resolve_solver_options from temoa.data_io.hybrid_loader import HybridLoader from temoa.data_processing.db_to_excel import make_excel from temoa.extensions.myopic.myopic_index import MyopicIndex @@ -268,7 +269,10 @@ def start(self) -> None: if not self.config.silent and self.progress_mapper and idx: self.progress_mapper.report(idx, 'solve') model, results = run_actions.solve_instance( - instance=instance, solver_name=self.config.solver_name, silent=True + instance=instance, + solver_name=self.config.solver_name, + silent=True, + solver_options=resolve_solver_options(self.config.solver), ) optimal, status = run_actions.check_solve_status(results) diff --git a/temoa/extensions/single_vector_mga/sv_mga_sequencer.py b/temoa/extensions/single_vector_mga/sv_mga_sequencer.py index 543af6678..d0266c317 100644 --- a/temoa/extensions/single_vector_mga/sv_mga_sequencer.py +++ b/temoa/extensions/single_vector_mga/sv_mga_sequencer.py @@ -19,6 +19,7 @@ from temoa.components.costs import total_cost_rule from temoa.core.config import TemoaConfig from temoa.core.model import TemoaModel +from temoa.core.solver_spec import resolve_solver_options from temoa.data_io.hybrid_loader import HybridLoader from temoa.extensions.single_vector_mga.output_summary import summarize from temoa.model_checking.pricing_check import price_checker @@ -100,6 +101,7 @@ def start(self) -> None: solver_name=self.config.solver_name, silent=self.config.silent, solver_suffixes=suffixes, + solver_options=resolve_solver_options(self.config.solver), ) status = res.solver.termination_condition logger.debug('Termination condition: %s', status.name) @@ -163,6 +165,7 @@ def start(self) -> None: solver_name=self.config.solver_name, silent=self.config.silent, solver_suffixes=suffixes, + solver_options=resolve_solver_options(self.config.solver), ) status = res.solver.termination_condition logger.debug('Termination condition: %s', status.name) diff --git a/temoa/extensions/stochastics/stochastic_sequencer.py b/temoa/extensions/stochastics/stochastic_sequencer.py index 5c904ce30..195cf691a 100644 --- a/temoa/extensions/stochastics/stochastic_sequencer.py +++ b/temoa/extensions/stochastics/stochastic_sequencer.py @@ -5,6 +5,7 @@ import pyomo.environ as pyo +from temoa.core.solver_spec import resolve_solver_options from temoa.extensions.stochastics.stochastic_config import StochasticConfig if TYPE_CHECKING: @@ -50,8 +51,13 @@ def start(self) -> None: from temoa.extensions.stochastics.scenario_creator import scenario_creator - # Merge solver options from stoch_config - solver_options = self.stoch_config.solver_options.get(self.config.solver_name, {}) + # Merge solver options: [solver.options] < stoch_config. Temoa defaults are excluded to + # preserve the behavior of stochastic runs prior to the [solver] table + solver_options = resolve_solver_options( + self.config.solver, + self.stoch_config.solver_options.get(self.config.solver_name, {}), + include_defaults=False, + ) options = { 'solver': self.config.solver_name, diff --git a/temoa/tutorial_assets/config_sample.toml b/temoa/tutorial_assets/config_sample.toml index f903f330d..24c6c60d3 100644 --- a/temoa/tutorial_assets/config_sample.toml +++ b/temoa/tutorial_assets/config_sample.toml @@ -75,8 +75,13 @@ neos = false # solver (Mandatory) # Depending on what client machine has installed. -# [appsi_highs, cbc, gurobi, cplex, ...] -solver_name = "appsi_highs" +# [appsi_highs, cbc, gurobi, cplex, ...] (any solver available through pyomo's SolverFactory) +# Either a solver name, which uses Temoa's default options for that solver: +solver = "appsi_highs" +# or a name plus options passed through to the solver, which are merged over Temoa's defaults +# (use an inline table here, as a [solver] table header would capture the keys that follow it): +# solver = { name = "gurobi", options = { Method = 2, Crossover = 0, BarConvTol = 1.0e-3, FeasibilityTol = 1.0e-4, BarOrder = -1 } } +# Note: 'solver_name' is deprecated but still accepted in place of 'solver' # ------------------------------------ # OUTPUTS From 6869621b1cf72e3af69a9376e0f81716b54225a9 Mon Sep 17 00:00:00 2001 From: Yamil Essus Date: Fri, 25 Sep 2026 14:46:14 -0400 Subject: [PATCH 04/10] Redact sensitive information from the configuration options that get logged --- temoa/core/config.py | 6 ++++-- temoa/core/solver_spec.py | 17 +++++++++++++++++ .../method_of_morris/morris_sequencer.py | 4 ++-- .../mga_sequencer.py | 6 ++++-- temoa/extensions/monte_carlo/mc_sequencer.py | 4 ++-- temoa/tutorial_assets/config_sample.toml | 4 +++- 6 files changed, 32 insertions(+), 9 deletions(-) diff --git a/temoa/core/config.py b/temoa/core/config.py index 26a095618..173b296b7 100644 --- a/temoa/core/config.py +++ b/temoa/core/config.py @@ -8,7 +8,7 @@ from typing import Any from temoa.core.modes import TemoaMode -from temoa.core.solver_spec import SolverSpec +from temoa.core.solver_spec import SolverSpec, redact_solver_options from temoa.extensions.framework import normalize_extension_ids, resolve_extension_specs logger = getLogger(__name__) @@ -397,7 +397,9 @@ def __repr__(self) -> str: msg += spacer msg += '{:>{}s}: {}\n'.format('Selected solver', width, self.solver_name) - msg += '{:>{}s}: {}\n'.format('Solver options', width, dict(self.solver.options)) + msg += '{:>{}s}: {}\n'.format( + 'Solver options', width, redact_solver_options(self.solver.options) + ) msg += '{:>{}s}: {}\n'.format('NEOS status', width, self.neos) msg += spacer diff --git a/temoa/core/solver_spec.py b/temoa/core/solver_spec.py index 45c6ca1f3..7985002e0 100644 --- a/temoa/core/solver_spec.py +++ b/temoa/core/solver_spec.py @@ -71,6 +71,23 @@ def parse(cls, raw: str | Mapping[str, Any] | SolverSpec) -> SolverSpec: raise TypeError(f'solver must be a str or a table/dict, got: {type(raw).__name__}') +# substrings (lowercase) of option names whose values are credentials, e.g. gurobi's WLSSecret, +# CloudSecretKey, CSAPIAccessID, ServerPassword, LicenseID +_SENSITIVE_OPTION_MARKERS = ('secret', 'password', 'accessid', 'licenseid', 'key') + + +def redact_solver_options(options: Mapping[str, Any]) -> dict[str, Any]: + """ + Return a copy of the options that is safe to log or print, with credential values masked + """ + return { + option: '***' + if any(marker in option.lower() for marker in _SENSITIVE_OPTION_MARKERS) + else option_value + for option, option_value in options.items() + } + + def resolve_solver_options( spec: SolverSpec, extension_options: Mapping[str, Any] | None = None, diff --git a/temoa/extensions/method_of_morris/morris_sequencer.py b/temoa/extensions/method_of_morris/morris_sequencer.py index 678417439..3a5baf194 100644 --- a/temoa/extensions/method_of_morris/morris_sequencer.py +++ b/temoa/extensions/method_of_morris/morris_sequencer.py @@ -22,7 +22,7 @@ from SALib.util import compute_groups_matrix, read_param_file # type: ignore[import-untyped] from temoa._internal.table_writer import TableWriter -from temoa.core.solver_spec import resolve_solver_options +from temoa.core.solver_spec import redact_solver_options, resolve_solver_options from temoa.data_io.hybrid_loader import HybridLoader from temoa.extensions.method_of_morris.morris_evaluate import evaluate @@ -76,7 +76,7 @@ def __init__(self, config: TemoaConfig): logger.warning('Unable to find solver options toml file. Using default options.') s_options = {} self.solver_options = resolve_solver_options(self.config.solver, s_options) - logger.info('Using solver options: %s', self.solver_options) + logger.info('Using solver options: %s', redact_solver_options(self.solver_options)) # output handling self.verbose = False # for troubleshooting diff --git a/temoa/extensions/modeling_to_generate_alternatives/mga_sequencer.py b/temoa/extensions/modeling_to_generate_alternatives/mga_sequencer.py index 6f3f124d3..ebdfae220 100644 --- a/temoa/extensions/modeling_to_generate_alternatives/mga_sequencer.py +++ b/temoa/extensions/modeling_to_generate_alternatives/mga_sequencer.py @@ -33,7 +33,7 @@ from temoa._internal.run_actions import build_instance from temoa._internal.table_writer import TableWriter from temoa.components.costs import total_cost_rule -from temoa.core.solver_spec import resolve_solver_options +from temoa.core.solver_spec import redact_solver_options, resolve_solver_options from temoa.data_io.hybrid_loader import HybridLoader from temoa.extensions.modeling_to_generate_alternatives.manager_factory import get_manager from temoa.extensions.modeling_to_generate_alternatives.mga_constants import MgaAxis, MgaWeighting @@ -108,7 +108,9 @@ def __init__(self, config: TemoaConfig): ).items(): self.opt.options[option] = option_value self.worker_solver_options = resolve_solver_options(self.config.solver, s_options) - logger.info('Using worker solver options: %s', self.worker_solver_options) + logger.info( + 'Using worker solver options: %s', redact_solver_options(self.worker_solver_options) + ) # some defaults, etc. self.internal_stop = False diff --git a/temoa/extensions/monte_carlo/mc_sequencer.py b/temoa/extensions/monte_carlo/mc_sequencer.py index 840df6030..c5318e23f 100644 --- a/temoa/extensions/monte_carlo/mc_sequencer.py +++ b/temoa/extensions/monte_carlo/mc_sequencer.py @@ -18,7 +18,7 @@ from typing import TYPE_CHECKING, Any, cast from temoa._internal.table_writer import TableWriter -from temoa.core.solver_spec import resolve_solver_options +from temoa.core.solver_spec import redact_solver_options, resolve_solver_options from temoa.data_io.hybrid_loader import HybridLoader from temoa.extensions.monte_carlo.mc_run import MCRun, MCRunFactory from temoa.extensions.monte_carlo.mc_worker import MCWorker @@ -96,7 +96,7 @@ def __init__(self, config: TemoaConfig): # worker options pulled from file self.num_workers = all_options.get('num_workers', 1) self.worker_solver_options = resolve_solver_options(self.config.solver, s_options) - logger.info('Using solver options: %s', self.worker_solver_options) + logger.info('Using solver options: %s', redact_solver_options(self.worker_solver_options)) # internal records self.solve_count = 0 diff --git a/temoa/tutorial_assets/config_sample.toml b/temoa/tutorial_assets/config_sample.toml index 24c6c60d3..d5f4dc8cd 100644 --- a/temoa/tutorial_assets/config_sample.toml +++ b/temoa/tutorial_assets/config_sample.toml @@ -77,10 +77,12 @@ neos = false # Depending on what client machine has installed. # [appsi_highs, cbc, gurobi, cplex, ...] (any solver available through pyomo's SolverFactory) # Either a solver name, which uses Temoa's default options for that solver: -solver = "appsi_highs" +solver = "gurobi" # or a name plus options passed through to the solver, which are merged over Temoa's defaults # (use an inline table here, as a [solver] table header would capture the keys that follow it): # solver = { name = "gurobi", options = { Method = 2, Crossover = 0, BarConvTol = 1.0e-3, FeasibilityTol = 1.0e-4, BarOrder = -1 } } +# Keep license credentials (e.g. gurobi WLS keys) in the solver's license file (gurobi.lic), not +# in these options. # Note: 'solver_name' is deprecated but still accepted in place of 'solver' # ------------------------------------ From ac5fe15e082317247aa9f4d2165edabc74d0db1d Mon Sep 17 00:00:00 2001 From: Yamil Essus Date: Mon, 28 Sep 2026 10:32:55 -0400 Subject: [PATCH 05/10] code rabbit comments --- temoa/core/solver_spec.py | 2 +- temoa/tutorial_assets/config_sample.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/temoa/core/solver_spec.py b/temoa/core/solver_spec.py index 7985002e0..800810228 100644 --- a/temoa/core/solver_spec.py +++ b/temoa/core/solver_spec.py @@ -73,7 +73,7 @@ def parse(cls, raw: str | Mapping[str, Any] | SolverSpec) -> SolverSpec: # substrings (lowercase) of option names whose values are credentials, e.g. gurobi's WLSSecret, # CloudSecretKey, CSAPIAccessID, ServerPassword, LicenseID -_SENSITIVE_OPTION_MARKERS = ('secret', 'password', 'accessid', 'licenseid', 'key') +_SENSITIVE_OPTION_MARKERS = ('secret', 'password', 'accessid', 'licenseid', 'key', 'token') def redact_solver_options(options: Mapping[str, Any]) -> dict[str, Any]: diff --git a/temoa/tutorial_assets/config_sample.toml b/temoa/tutorial_assets/config_sample.toml index d5f4dc8cd..b262cf598 100644 --- a/temoa/tutorial_assets/config_sample.toml +++ b/temoa/tutorial_assets/config_sample.toml @@ -77,7 +77,7 @@ neos = false # Depending on what client machine has installed. # [appsi_highs, cbc, gurobi, cplex, ...] (any solver available through pyomo's SolverFactory) # Either a solver name, which uses Temoa's default options for that solver: -solver = "gurobi" +solver = "appsi_highs" # or a name plus options passed through to the solver, which are merged over Temoa's defaults # (use an inline table here, as a [solver] table header would capture the keys that follow it): # solver = { name = "gurobi", options = { Method = 2, Crossover = 0, BarConvTol = 1.0e-3, FeasibilityTol = 1.0e-4, BarOrder = -1 } } From 1df3264ddb878eaec18f160597a0525a5711b7c4 Mon Sep 17 00:00:00 2001 From: Davey Elder Date: Mon, 28 Sep 2026 13:50:52 -0400 Subject: [PATCH 06/10] Add test coverage for CLI, myopic, stochastic, and extensions framework Signed-off-by: Davey Elder --- tests/test_cli.py | 192 +++++++++++++ tests/test_evolution_updater.py | 17 ++ tests/test_framework_extension_helpers.py | 326 ++++++++++++++++++++++ tests/test_myopic_progress_mapper.py | 87 ++++++ tests/test_stochastic_sequencer.py | 53 ++++ 5 files changed, 675 insertions(+) create mode 100644 tests/test_evolution_updater.py create mode 100644 tests/test_framework_extension_helpers.py create mode 100644 tests/test_myopic_progress_mapper.py create mode 100644 tests/test_stochastic_sequencer.py diff --git a/tests/test_cli.py b/tests/test_cli.py index 5ba1fa1ee..73d0192dd 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -4,6 +4,7 @@ from pathlib import Path import pytest +import tomlkit from typer.testing import CliRunner from temoa.cli import _is_writable, app @@ -469,3 +470,194 @@ def test_cli_run_fails_if_solver_missing(tmp_path: Path, monkeypatch: pytest.Mon # Use the more robust phrase for checking installation instructions assert 'Please ensure the solver is installed and accessible.' in result.stdout assert (tmp_path / 'temoa-run.log').exists() + + +# ============================================================================= +# Tests for the `tutorial` command +# ============================================================================= + + +def test_cli_tutorial_creates_files(tmp_path: Path, monkeypatch: pytest.MonkeyPatch) -> None: + """Test that `temoa tutorial` creates the config, database, and mc_settings files.""" + monkeypatch.chdir(tmp_path) + + args = ['tutorial', 'my_config', 'my_database'] + result = runner.invoke(app, args, catch_exceptions=False) + + assert result.exit_code == 0, f'CLI crashed with error: {result.exception}\n{result.stdout}' + assert (tmp_path / 'my_config.toml').exists() + assert (tmp_path / 'my_database.sqlite').exists() + assert (tmp_path / 'mc_settings.csv').exists() + assert 'Tutorial Setup Complete!' in result.stdout + + +def test_cli_tutorial_default_names(tmp_path: Path, monkeypatch: pytest.MonkeyPatch) -> None: + """Test that `temoa tutorial` uses its default file names when none are given.""" + monkeypatch.chdir(tmp_path) + + result = runner.invoke(app, ['tutorial'], catch_exceptions=False) + + assert result.exit_code == 0, f'CLI crashed with error: {result.exception}\n{result.stdout}' + assert (tmp_path / 'tutorial_config.toml').exists() + assert (tmp_path / 'tutorial_database.sqlite').exists() + + +def test_cli_tutorial_updates_toml_database_paths( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch +) -> None: + """Test that the generated config file points at the newly created database.""" + monkeypatch.chdir(tmp_path) + + result = runner.invoke(app, ['tutorial', 'cfg', 'db_name'], catch_exceptions=False) + + assert result.exit_code == 0, f'CLI crashed with error: {result.exception}\n{result.stdout}' + doc = tomlkit.parse((tmp_path / 'cfg.toml').read_text()) + assert doc['input_database'] == 'db_name.sqlite' + assert doc['output_database'] == 'db_name.sqlite' + + +def test_cli_tutorial_verbose_output(tmp_path: Path, monkeypatch: pytest.MonkeyPatch) -> None: + """Test that `--verbose` prints the extra progress and guidance messages.""" + monkeypatch.chdir(tmp_path) + + result = runner.invoke(app, ['tutorial', 'cfg', 'db', '--verbose'], catch_exceptions=False) + + assert result.exit_code == 0 + assert 'Copying tutorial resources...' in result.stdout + assert 'Updating database paths in configuration...' in result.stdout + assert 'Tutorial files created successfully' in result.stdout + + +def test_cli_tutorial_existing_files_aborts_without_force( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch +) -> None: + """Test that existing tutorial files trigger a confirmation prompt that can be declined.""" + monkeypatch.chdir(tmp_path) + (tmp_path / 'cfg.toml').write_text('placeholder') + + result = runner.invoke(app, ['tutorial', 'cfg', 'db'], input='n\n') + + # A declined confirmation is a graceful, non-error cancellation. + assert result.exit_code == 0 + assert 'Tutorial files already exist' in result.stdout + assert 'Tutorial setup cancelled' in result.stdout + # The placeholder file should be untouched since the user declined. + assert (tmp_path / 'cfg.toml').read_text() == 'placeholder' + + +def test_cli_tutorial_existing_files_force_overwrite( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch +) -> None: + """Test that `--force` overwrites existing tutorial files without prompting.""" + monkeypatch.chdir(tmp_path) + (tmp_path / 'cfg.toml').write_text('placeholder') + + result = runner.invoke(app, ['tutorial', 'cfg', 'db', '--force'], catch_exceptions=False) + + assert result.exit_code == 0, f'CLI crashed with error: {result.exception}\n{result.stdout}' + assert (tmp_path / 'db.sqlite').exists() + assert (tmp_path / 'cfg.toml').read_text() != 'placeholder' + assert 'Tutorial setup cancelled' not in result.stdout + + +# ============================================================================= +# Tests for the `check-units` command +# ============================================================================= + +VALID_UNITS_DB = Path(__file__).parent / 'testing_outputs' / 'utopia_valid_units.sqlite' +INVALID_CURRENCY_DB = Path(__file__).parent / 'testing_outputs' / 'utopia_invalid_currency.sqlite' + +requires_unit_dbs = pytest.mark.skipif( + not (VALID_UNITS_DB.exists() and INVALID_CURRENCY_DB.exists()), + reason='Test databases not created. Ensure conftest.py setup completed successfully.', +) + + +@requires_unit_dbs +def test_cli_check_units_all_clear(tmp_path: Path) -> None: + """Test `temoa check-units` reports success on a valid database.""" + args = ['check-units', str(VALID_UNITS_DB), '--output', str(tmp_path)] + result = runner.invoke(app, args, catch_exceptions=False) + + assert result.exit_code == 0, f'CLI crashed with error: {result.exception}\n{result.stdout}' + assert 'All unit checks passed' in result.stdout + + +@requires_unit_dbs +def test_cli_check_units_all_clear_silent(tmp_path: Path) -> None: + """Test that `--silent` suppresses the success message.""" + args = ['check-units', str(VALID_UNITS_DB), '--output', str(tmp_path), '--silent'] + result = runner.invoke(app, args, catch_exceptions=False) + + assert result.exit_code == 0 + assert 'All unit checks passed' not in result.stdout + + +@requires_unit_dbs +def test_cli_check_units_detects_issues(tmp_path: Path) -> None: + """Test that `temoa check-units` fails and writes a report for a bad database.""" + args = ['check-units', str(INVALID_CURRENCY_DB), '--output', str(tmp_path)] + result = runner.invoke(app, args, catch_exceptions=False) + + assert result.exit_code != 0 + assert 'Unit check found issues' in result.stdout + assert 'Detailed report saved to' in result.stdout + assert 'Report Summary:' in result.stdout + reports = list(tmp_path.glob('units_check_*.txt')) + assert len(reports) == 1 + + +@requires_unit_dbs +def test_cli_check_units_detects_issues_silent(tmp_path: Path) -> None: + """Test that `--silent` suppresses the issue report summary but still fails and writes it.""" + args = ['check-units', str(INVALID_CURRENCY_DB), '--output', str(tmp_path), '--silent'] + result = runner.invoke(app, args, catch_exceptions=False) + + assert result.exit_code != 0 + assert 'Unit check found issues' not in result.stdout + reports = list(tmp_path.glob('units_check_*.txt')) + assert len(reports) == 1 + + +@requires_unit_dbs +def test_cli_check_units_default_output_dir( + tmp_path: Path, monkeypatch: pytest.MonkeyPatch +) -> None: + """Test that omitting `--output` defaults the report to ./unit_check_reports.""" + monkeypatch.chdir(tmp_path) + + args = ['check-units', str(VALID_UNITS_DB)] + result = runner.invoke(app, args, catch_exceptions=False) + + assert result.exit_code == 0 + assert (tmp_path / 'unit_check_reports').is_dir() + + +def test_cli_check_units_missing_database() -> None: + """Test graceful failure for a missing database file.""" + args = ['check-units', 'non_existent_db.sqlite'] + result = runner.invoke(app, args) + + assert result.exit_code != 0 + assert 'non_existent_db.sqlite' in result.stderr + + +# ============================================================================= +# Tests for `_is_writable` +# ============================================================================= + + +def test_is_writable_true_for_writable_dir(tmp_path: Path) -> None: + """Test that a normal, writable directory is reported as writable.""" + assert _is_writable(tmp_path) is True + + +def test_is_writable_false_on_oserror(monkeypatch: pytest.MonkeyPatch, tmp_path: Path) -> None: + """Test that `_is_writable` returns False when touching the probe file raises OSError.""" + + def _raise_oserror(*_args: object, **_kwargs: object) -> None: + raise OSError('mocked failure') + + monkeypatch.setattr(Path, 'touch', _raise_oserror) + + assert _is_writable(tmp_path) is False diff --git a/tests/test_evolution_updater.py b/tests/test_evolution_updater.py new file mode 100644 index 000000000..6d3efbaa9 --- /dev/null +++ b/tests/test_evolution_updater.py @@ -0,0 +1,17 @@ +"""Tests for the myopic evolution_updater template module.""" + +import logging + +import pytest + +from temoa.extensions.myopic.evolution_updater import iterate +from temoa.extensions.myopic.myopic_index import MyopicIndex + + +def test_iterate_logs_base_year(caplog: pytest.LogCaptureFixture) -> None: + idx = MyopicIndex(base_year=2020, step_year=2025, last_demand_year=2024, last_year=2030) + + with caplog.at_level(logging.INFO): + iterate(idx=idx, prev_base_year=2015, last_instance_status='optimal', db_con=None) + + assert 'base year 2020' in caplog.text diff --git a/tests/test_framework_extension_helpers.py b/tests/test_framework_extension_helpers.py new file mode 100644 index 000000000..4bc965322 --- /dev/null +++ b/tests/test_framework_extension_helpers.py @@ -0,0 +1,326 @@ +"""Tests for the extension-agnostic helper functions in temoa.extensions.framework. + +These cover the plumbing (id normalization, manifest/hook merging, and the +enabled/disabled extension table checks) that isn't exercised by +tests/test_extensions.py, which focuses on each concrete extension's model +components. +""" + +from __future__ import annotations + +import sqlite3 +from typing import TYPE_CHECKING, cast + +import pytest + +from temoa.extensions.framework import ( + ExtensionSpec, + _append_extension_schema, + _table_exists, + _table_has_rows, + append_extension_manifest_items, + apply_model_extension_hooks, + assert_disabled_extension_tables_are_empty, + ensure_enabled_extension_tables_exist, + get_known_extension_specs, + merge_regional_group_tables, + normalize_extension_ids, +) + +if TYPE_CHECKING: + from pathlib import Path + + from temoa.data_io.loader_manifest import LoadItem + + +# ============================================================================= +# normalize_extension_ids +# ============================================================================= + + +def test_normalize_extension_ids_none_returns_empty() -> None: + assert normalize_extension_ids(None) == () + + +def test_normalize_extension_ids_empty_list_returns_empty() -> None: + assert normalize_extension_ids([]) == () + + +def test_normalize_extension_ids_dedupes_and_lowercases_preserving_order() -> None: + result = normalize_extension_ids(['Growth_Rates', ' growth_rates ', 'discrete_capacity']) + assert result == ('growth_rates', 'discrete_capacity') + + +def test_normalize_extension_ids_skips_blank_entries() -> None: + assert normalize_extension_ids([' ', 'growth_rates']) == ('growth_rates',) + + +def test_normalize_extension_ids_rejects_non_string() -> None: + with pytest.raises(TypeError, match='Extension ids must be strings'): + normalize_extension_ids([123]) + + +# ============================================================================= +# merge_regional_group_tables +# ============================================================================= + + +def test_merge_regional_group_tables_merges_specs_into_base() -> None: + spec = ExtensionSpec(extension_id='ext_a', regional_group_tables={'tbl_a': 'field_a'}) + merged = merge_regional_group_tables({'tbl_base': 'field_base'}, [spec]) + assert merged == {'tbl_base': 'field_base', 'tbl_a': 'field_a'} + + +def test_merge_regional_group_tables_allows_identical_duplicate_mapping() -> None: + spec = ExtensionSpec(extension_id='ext_a', regional_group_tables={'tbl_base': 'field_base'}) + merged = merge_regional_group_tables({'tbl_base': 'field_base'}, [spec]) + assert merged == {'tbl_base': 'field_base'} + + +def test_merge_regional_group_tables_conflict_raises() -> None: + spec = ExtensionSpec(extension_id='ext_a', regional_group_tables={'tbl_base': 'field_other'}) + with pytest.raises(ValueError, match='conflicting field mappings'): + merge_regional_group_tables({'tbl_base': 'field_base'}, [spec]) + + +# ============================================================================= +# apply_model_extension_hooks / append_extension_manifest_items +# ============================================================================= + + +def test_apply_model_extension_hooks_calls_each_registered_hook() -> None: + calls: list[object] = [] + spec = ExtensionSpec(extension_id='ext_a', register_model_components=calls.append) + model = object() + + apply_model_extension_hooks(model, [spec]) # type: ignore[arg-type] + + assert calls == [model] + + +def test_apply_model_extension_hooks_skips_specs_without_hook() -> None: + spec = ExtensionSpec(extension_id='ext_a') + # Should not raise even though register_model_components is None. + apply_model_extension_hooks(object(), [spec]) # type: ignore[arg-type] + + +def test_append_extension_manifest_items_merges_in_order() -> None: + # Strings stand in for LoadItem; only list order matters here. + item_a = cast('LoadItem', 'item_a') + item_b = cast('LoadItem', 'item_b') + base_item = cast('LoadItem', 'base_item') + spec_a = ExtensionSpec(extension_id='ext_a', build_manifest_items=lambda _model: [item_a]) + spec_b = ExtensionSpec(extension_id='ext_b', build_manifest_items=lambda _model: [item_b]) + + merged = append_extension_manifest_items( + object(), # type: ignore[arg-type] + [base_item], + [spec_a, spec_b], + ) + + assert merged == [base_item, item_a, item_b] + + +# ============================================================================= +# _table_exists / _table_has_rows +# ============================================================================= + + +def test_table_exists_and_has_rows() -> None: + con = sqlite3.connect(':memory:') + try: + assert _table_exists(con, 'missing_table') is False + assert _table_has_rows(con, 'missing_table') is False + + con.execute('CREATE TABLE populated (id INTEGER)') + con.execute('CREATE TABLE empty_table (id INTEGER)') + con.execute('INSERT INTO populated VALUES (1)') + con.commit() + + assert _table_exists(con, 'populated') is True + assert _table_has_rows(con, 'populated') is True + assert _table_exists(con, 'empty_table') is True + assert _table_has_rows(con, 'empty_table') is False + finally: + con.close() + + +# ============================================================================= +# assert_disabled_extension_tables_are_empty +# ============================================================================= + + +def test_assert_disabled_extension_tables_are_empty_warns_when_populated( + caplog: pytest.LogCaptureFixture, +) -> None: + con = sqlite3.connect(':memory:') + try: + con.execute('CREATE TABLE limit_growth_capacity (region TEXT)') + con.execute("INSERT INTO limit_growth_capacity VALUES ('R1')") + con.commit() + + with caplog.at_level('WARNING'): + assert_disabled_extension_tables_are_empty(con, enabled_specs=()) + + assert any('growth_rates' in record.message for record in caplog.records) + finally: + con.close() + + +def test_assert_disabled_extension_tables_are_empty_silent_when_enabled( + caplog: pytest.LogCaptureFixture, +) -> None: + con = sqlite3.connect(':memory:') + try: + con.execute('CREATE TABLE limit_growth_capacity (region TEXT)') + con.execute("INSERT INTO limit_growth_capacity VALUES ('R1')") + con.commit() + + growth_rates_spec = get_known_extension_specs()['growth_rates'] + with caplog.at_level('WARNING'): + assert_disabled_extension_tables_are_empty(con, enabled_specs=(growth_rates_spec,)) + + assert not caplog.records + finally: + con.close() + + +# ============================================================================= +# ensure_enabled_extension_tables_exist / _append_extension_schema +# ============================================================================= + + +def test_ensure_enabled_extension_tables_exist_noop_when_tables_present() -> None: + con = sqlite3.connect(':memory:') + try: + con.execute('CREATE TABLE owned_table (id INTEGER)') + con.commit() + spec = ExtensionSpec(extension_id='ext_a', owned_tables=('owned_table',)) + + # Should not raise or prompt. + ensure_enabled_extension_tables_exist(con, [spec], input_database='db.sqlite', silent=True) + finally: + con.close() + + +def test_ensure_enabled_extension_tables_exist_no_schema_path_raises() -> None: + con = sqlite3.connect(':memory:') + try: + spec = ExtensionSpec(extension_id='ext_a', owned_tables=('missing_table',)) + with pytest.raises(RuntimeError, match='No schema SQL path is registered'): + ensure_enabled_extension_tables_exist( + con, [spec], input_database='db.sqlite', silent=True + ) + finally: + con.close() + + +def test_ensure_enabled_extension_tables_exist_silent_skips_prompt_and_raises() -> None: + con = sqlite3.connect(':memory:') + try: + spec = ExtensionSpec( + extension_id='ext_a', owned_tables=('missing_table',), schema_sql_path='unused.sql' + ) + # silent=True means the prompt is never asked, so should_apply stays False. + with pytest.raises(RuntimeError, match='Re-run and accept the prompt'): + ensure_enabled_extension_tables_exist( + con, [spec], input_database='db.sqlite', silent=True + ) + finally: + con.close() + + +def test_ensure_enabled_extension_tables_exist_prompt_declined_raises( + monkeypatch: pytest.MonkeyPatch, +) -> None: + con = sqlite3.connect(':memory:') + try: + spec = ExtensionSpec( + extension_id='ext_a', owned_tables=('missing_table',), schema_sql_path='unused.sql' + ) + monkeypatch.setattr('builtins.input', lambda _prompt: 'n') + with pytest.raises(RuntimeError, match='Re-run and accept the prompt'): + ensure_enabled_extension_tables_exist( + con, [spec], input_database='db.sqlite', silent=False + ) + finally: + con.close() + + +def test_ensure_enabled_extension_tables_exist_prompt_accepted_applies_schema( + monkeypatch: pytest.MonkeyPatch, tmp_path: Path +) -> None: + schema_file = tmp_path / 'extra_schema.sql' + schema_file.write_text('CREATE TABLE missing_table (id INTEGER);') + + con = sqlite3.connect(':memory:') + try: + spec = ExtensionSpec( + extension_id='ext_a', + owned_tables=('missing_table',), + schema_sql_path=str(schema_file), + ) + monkeypatch.setattr('builtins.input', lambda _prompt: 'y') + + ensure_enabled_extension_tables_exist(con, [spec], input_database='db.sqlite', silent=False) + + assert _table_exists(con, 'missing_table') is True + finally: + con.close() + + +def test_ensure_enabled_extension_tables_exist_still_missing_after_apply_raises( + monkeypatch: pytest.MonkeyPatch, tmp_path: Path +) -> None: + # Schema file exists but doesn't actually create the owned table. + schema_file = tmp_path / 'noop_schema.sql' + schema_file.write_text('CREATE TABLE unrelated_table (id INTEGER);') + + con = sqlite3.connect(':memory:') + try: + spec = ExtensionSpec( + extension_id='ext_a', + owned_tables=('missing_table',), + schema_sql_path=str(schema_file), + ) + monkeypatch.setattr('builtins.input', lambda _prompt: 'y') + + with pytest.raises(RuntimeError, match='still missing'): + ensure_enabled_extension_tables_exist( + con, [spec], input_database='db.sqlite', silent=False + ) + finally: + con.close() + + +def test_append_extension_schema_no_path_raises() -> None: + con = sqlite3.connect(':memory:') + try: + spec = ExtensionSpec(extension_id='ext_a') + with pytest.raises(RuntimeError, match='no schema SQL path configured'): + _append_extension_schema(con, spec) + finally: + con.close() + + +def test_append_extension_schema_missing_file_raises() -> None: + con = sqlite3.connect(':memory:') + try: + spec = ExtensionSpec(extension_id='ext_a', schema_sql_path='/no/such/file.sql') + with pytest.raises(FileNotFoundError, match='not found'): + _append_extension_schema(con, spec) + finally: + con.close() + + +def test_append_extension_schema_executes_and_commits(tmp_path: Path) -> None: + schema_file = tmp_path / 'schema.sql' + schema_file.write_text('CREATE TABLE new_table (id INTEGER);') + + con = sqlite3.connect(':memory:') + try: + spec = ExtensionSpec(extension_id='ext_a', schema_sql_path=str(schema_file)) + _append_extension_schema(con, spec) + assert _table_exists(con, 'new_table') is True + finally: + con.close() diff --git a/tests/test_myopic_progress_mapper.py b/tests/test_myopic_progress_mapper.py new file mode 100644 index 000000000..db024d520 --- /dev/null +++ b/tests/test_myopic_progress_mapper.py @@ -0,0 +1,87 @@ +"""Tests for MyopicProgressMapper, the console progress visualizer for myopic solves.""" + +import re + +import pytest + +from temoa.extensions.myopic.myopic_index import MyopicIndex +from temoa.extensions.myopic.myopic_progress_mapper import MyopicProgressMapper + +YEARS = [2020, 2025, 2030, 2035] + + +def _index(base_year: int, step_year: int, last_demand_year: int) -> MyopicIndex: + return MyopicIndex( + base_year=base_year, + step_year=step_year, + last_demand_year=last_demand_year, + last_year=YEARS[-1] + 1, + ) + + +def test_init_computes_tag_width_and_positions() -> None: + mapper = MyopicProgressMapper(YEARS) + + assert mapper.years == YEARS + assert mapper.tag_width == max(len(str(y)) for y in YEARS) + 2 * len(mapper.leader) + # Positions are in increasing order, one per year. + assert list(mapper.pos.keys()) == YEARS + assert all(mapper.pos[YEARS[i]] < mapper.pos[YEARS[i + 1]] for i in range(len(YEARS) - 1)) + + +def test_draw_header_prints_years_and_label(capsys: pytest.CaptureFixture[str]) -> None: + mapper = MyopicProgressMapper(YEARS) + mapper.draw_header() + + out = capsys.readouterr().out + assert 'Myopic Progress' in out + assert 'HH:MM:SS' in out + for year in YEARS: + assert str(year) in out + + +def test_timestamp_format() -> None: + mapper = MyopicProgressMapper(YEARS) + assert re.match(r'^Elapsed: \d{2}:\d{2}:\d{2}\s+$', mapper.timestamp()) + + +@pytest.mark.parametrize( + 'status,tag', + [ + ('load', 'LOAD'), + ('solve', 'SOLV'), + ('check', 'CHEK'), + ('evolve', 'EVLV'), + ], +) +def test_report_prints_expected_tag_for_status( + capsys: pytest.CaptureFixture[str], status: str, tag: str +) -> None: + mapper = MyopicProgressMapper(YEARS) + idx = _index(base_year=2020, step_year=2025, last_demand_year=2025) + + mapper.report(idx, status) # type: ignore[arg-type] + + out = capsys.readouterr().out + # One tag per year from base_year through last_demand_year (2020, 2025). + assert out.count(tag) == 2 + assert 'Elapsed:' in out + + +def test_report_status_report_uses_step_year(capsys: pytest.CaptureFixture[str]) -> None: + mapper = MyopicProgressMapper(YEARS) + idx = _index(base_year=2020, step_year=2030, last_demand_year=2025) + + mapper.report(idx, 'report') + + out = capsys.readouterr().out + # One tag per year from base_year up to (not including) step_year: 2020, 2025. + assert out.count('RECD') == 2 + + +def test_report_rejects_invalid_status() -> None: + mapper = MyopicProgressMapper(YEARS) + idx = _index(base_year=2020, step_year=2025, last_demand_year=2025) + + with pytest.raises(ValueError, match='bad status'): + mapper.report(idx, 'bogus') # type: ignore[arg-type] diff --git a/tests/test_stochastic_sequencer.py b/tests/test_stochastic_sequencer.py new file mode 100644 index 000000000..84f506d2a --- /dev/null +++ b/tests/test_stochastic_sequencer.py @@ -0,0 +1,53 @@ +"""Tests for StochasticSequencer's constructor validation of stochastic config files.""" + +from pathlib import Path +from types import SimpleNamespace + +import pytest + +from temoa.extensions.stochastics.stochastic_sequencer import StochasticSequencer + + +def _config(stochastic_config: Path | None) -> SimpleNamespace: + """A minimal stand-in for TemoaConfig with just the attribute the sequencer reads.""" + return SimpleNamespace(stochastic_config=stochastic_config) + + +def test_missing_stochastic_config_raises() -> None: + with pytest.raises(ValueError, match="requires a 'stochastic_config'"): + StochasticSequencer(_config(None)) # type: ignore[arg-type] + + +def test_nonexistent_stochastic_config_path_raises(tmp_path: Path) -> None: + missing = tmp_path / 'does_not_exist.toml' + with pytest.raises(ValueError, match='not found'): + StochasticSequencer(_config(missing)) # type: ignore[arg-type] + + +def test_stochastic_config_path_is_directory_raises(tmp_path: Path) -> None: + with pytest.raises(ValueError, match='is not a file'): + StochasticSequencer(_config(tmp_path)) # type: ignore[arg-type] + + +def test_invalid_toml_content_raises_wrapped_error(tmp_path: Path) -> None: + bad_toml = tmp_path / 'stoch.toml' + bad_toml.write_text('not valid = toml = content [[[') + + with pytest.raises(ValueError, match='Error parsing stochastic config'): + StochasticSequencer(_config(bad_toml)) # type: ignore[arg-type] + + +def test_valid_stochastic_config_loads_successfully(tmp_path: Path) -> None: + good_toml = tmp_path / 'stoch.toml' + good_toml.write_text( + """ + [scenarios] + base = 0.5 + high = 0.5 + """ + ) + + sequencer = StochasticSequencer(_config(good_toml)) # type: ignore[arg-type] + + assert sequencer.stoch_config.scenarios == {'base': 0.5, 'high': 0.5} + assert sequencer.objective_value is None From 2d83cc26b51410ce72a6abfea74cda93fbe765ec Mon Sep 17 00:00:00 2001 From: Davey Elder Date: Mon, 28 Sep 2026 13:51:13 -0400 Subject: [PATCH 07/10] Stop refreshing databases for test collection Signed-off-by: Davey Elder --- tests/conftest.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/tests/conftest.py b/tests/conftest.py index 444aebd50..7cb63eb9a 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -185,8 +185,11 @@ def create_unit_test_dbs() -> None: logger.info('Created unit test DB: %s', db_name) -def pytest_configure(config: Config) -> None: # noqa: ARG001 +def pytest_configure(config: Config) -> None: """Setup test databases before test collection.""" + # Skip for collect-only (e.g. IDE test discovery) so it can't race a real run for DB locks. + if config.getoption('collectonly'): + return refresh_databases() try: create_unit_test_dbs() From 1002c32f6c81e6930748093e40587c88f916a0ab Mon Sep 17 00:00:00 2001 From: Davey Elder Date: Mon, 28 Sep 2026 15:24:31 -0400 Subject: [PATCH 08/10] Specify exact expected exit codes in test_cli Signed-off-by: Davey Elder --- tests/test_cli.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/tests/test_cli.py b/tests/test_cli.py index 73d0192dd..6e1551b9a 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -128,7 +128,7 @@ def test_cli_validate_failure_on_invalid_db(tmp_path: Path) -> None: args = ['validate', str(test_config_path), '--output', str(tmp_path)] result = runner.invoke(app, args) - assert result.exit_code != 0, 'CLI should exit with a non-zero code on failure' + assert result.exit_code == 1, 'validate should exit 1 on a failed validation' assert 'Validation failed' in result.stdout # Check that the log was still created, containing the detailed error assert (tmp_path / 'temoa-run.log').exists() @@ -139,7 +139,7 @@ def test_cli_run_missing_config() -> None: args = ['run', 'non_existent_file.toml'] result = runner.invoke(app, args) - assert result.exit_code != 0 + assert result.exit_code == 2, 'missing config file should be rejected as a bad argument' # Check that the error mentions the missing file (more robust than exact string match) assert 'non_existent_file.toml' in result.stderr @@ -293,7 +293,7 @@ def mock_is_writable_always_false(_path: Path) -> bool: args = ['migrate', str(input_file)] result = runner.invoke(app, args, catch_exceptions=False) - assert result.exit_code != 0, 'Migration should fail with a non-zero exit code' + assert result.exit_code == 1, 'migrate should exit 1 when no writable output location exists' # Normalize whitespace to handle platform-specific line breaks from rich.print() normalized_output = ' '.join(result.stdout.split()) assert 'Error: Neither input directory' in normalized_output @@ -309,7 +309,7 @@ def test_cli_migrate_invalid_file() -> None: args = ['migrate', 'non_existent.sql'] result = runner.invoke(app, args) - assert result.exit_code != 0 + assert result.exit_code == 2, 'missing input file should be rejected as a bad argument' # Typer handles file existence check, so error is in stderr assert 'does not exist' in result.stderr or 'does not exist' in str(result.exception) @@ -321,7 +321,7 @@ def test_cli_migrate_unknown_type(tmp_path: Path) -> None: args = ['migrate', str(unknown_file)] result = runner.invoke(app, args) - assert result.exit_code != 0 + assert result.exit_code == 1, 'migrate should exit 1 for an undeterminable migration type' assert 'Cannot determine migration type' in result.stdout @@ -434,7 +434,7 @@ def test_cli_validate_fails_if_solver_missing( args = ['validate', str(test_config_path), '--output', str(tmp_path)] result = runner.invoke(app, args, catch_exceptions=False) - assert result.exit_code != 0, ( + assert result.exit_code == 1, ( f'Validate should have failed: {result.exception}\n{result.stderr}\n{result.stdout}' ) assert isinstance(result.exception, SystemExit) @@ -460,7 +460,7 @@ def test_cli_run_fails_if_solver_missing(tmp_path: Path, monkeypatch: pytest.Mon args = ['run', str(test_config_path), '--output', str(tmp_path)] result = runner.invoke(app, args, catch_exceptions=False) - assert result.exit_code != 0, ( + assert result.exit_code == 1, ( f'Run should have failed: {result.exception}\n{result.stderr}\n{result.stdout}' ) assert isinstance(result.exception, SystemExit) @@ -599,7 +599,7 @@ def test_cli_check_units_detects_issues(tmp_path: Path) -> None: args = ['check-units', str(INVALID_CURRENCY_DB), '--output', str(tmp_path)] result = runner.invoke(app, args, catch_exceptions=False) - assert result.exit_code != 0 + assert result.exit_code == 1, 'check-units should exit 1 when issues are found' assert 'Unit check found issues' in result.stdout assert 'Detailed report saved to' in result.stdout assert 'Report Summary:' in result.stdout @@ -613,7 +613,7 @@ def test_cli_check_units_detects_issues_silent(tmp_path: Path) -> None: args = ['check-units', str(INVALID_CURRENCY_DB), '--output', str(tmp_path), '--silent'] result = runner.invoke(app, args, catch_exceptions=False) - assert result.exit_code != 0 + assert result.exit_code == 1, 'check-units should exit 1 when issues are found' assert 'Unit check found issues' not in result.stdout reports = list(tmp_path.glob('units_check_*.txt')) assert len(reports) == 1 @@ -638,7 +638,7 @@ def test_cli_check_units_missing_database() -> None: args = ['check-units', 'non_existent_db.sqlite'] result = runner.invoke(app, args) - assert result.exit_code != 0 + assert result.exit_code == 2, 'missing database file should be rejected as a bad argument' assert 'non_existent_db.sqlite' in result.stderr From a676abd4c98917e3f3d019363cd8dd74266f68b0 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 2 Oct 2026 13:15:46 +0000 Subject: [PATCH 09/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- scripts/electricity_prices.py | 166 ++++++++++++++++++++++++---------- scripts/reserve_dual_toy.py | 70 ++++++++++---- 2 files changed, 168 insertions(+), 68 deletions(-) diff --git a/scripts/electricity_prices.py b/scripts/electricity_prices.py index df6262b73..ff3e7e65e 100644 --- a/scripts/electricity_prices.py +++ b/scripts/electricity_prices.py @@ -38,8 +38,23 @@ import pandas as pd # Multipliers to convert a commodity unit to MWh and a cost unit to dollars. -ENERGY_TO_MWH = {'PJ': 1e15 / 3.6e9, 'TJ': 1e12 / 3.6e9, 'GJ': 1 / 3.6, 'TWH': 1e6, 'GWH': 1e3, 'MWH': 1.0} -COST_TO_USD = {'MUSD': 1e6, 'M$': 1e6, 'MDOLLAR': 1e6, 'BUSD': 1e9, 'KUSD': 1e3, 'USD': 1.0, '$': 1.0} +ENERGY_TO_MWH = { + 'PJ': 1e15 / 3.6e9, + 'TJ': 1e12 / 3.6e9, + 'GJ': 1 / 3.6, + 'TWH': 1e6, + 'GWH': 1e3, + 'MWH': 1.0, +} +COST_TO_USD = { + 'MUSD': 1e6, + 'M$': 1e6, + 'MDOLLAR': 1e6, + 'BUSD': 1e9, + 'KUSD': 1e3, + 'USD': 1.0, + '$': 1.0, +} NAME_RE = re.compile(r'^(?P\w+)\[(?P.*)\]$') @@ -52,14 +67,18 @@ def connect(path: Path) -> sqlite3.Connection: return sqlite3.connect(f'file:{path}?mode=ro&immutable=1', uri=True) -def discount_factors(con: sqlite3.Connection, base_year: int | None) -> tuple[pd.DataFrame, float, int]: +def discount_factors( + con: sqlite3.Connection, base_year: int | None +) -> tuple[pd.DataFrame, float, int]: row = con.execute( "SELECT value FROM metadata_real WHERE element = 'global_discount_rate'" ).fetchone() if row is None: sys.exit('global_discount_rate missing from metadata_real (Temoa requires it).') gdr = float(row[0]) - future = [r[0] for r in con.execute("SELECT period FROM time_period WHERE flag = 'f' ORDER BY period")] + future = [ + r[0] for r in con.execute("SELECT period FROM time_period WHERE flag = 'f' ORDER BY period") + ] p0 = base_year if base_year is not None else future[0] rows = [] for p, p_next in zip(future[:-1], future[1:]): @@ -73,13 +92,14 @@ def discount_factors(con: sqlite3.Connection, base_year: int | None) -> tuple[pd return pd.DataFrame(rows), gdr, p0 -def unit_conversion(con: sqlite3.Connection, commodity: str, cost_scale: float | None, - energy_to_mwh: float | None) -> tuple[float, str]: +def unit_conversion( + con: sqlite3.Connection, commodity: str, cost_scale: float | None, energy_to_mwh: float | None +) -> tuple[float, str]: """Return multiplier taking (cost unit / commodity unit) to $/MWh.""" c_units = con.execute('SELECT units FROM commodity WHERE name = ?', (commodity,)).fetchone() c_units = (c_units[0] or '').strip() if c_units else '' cost_units = con.execute( - "SELECT units FROM output_cost WHERE units IS NOT NULL LIMIT 1" + 'SELECT units FROM output_cost WHERE units IS NOT NULL LIMIT 1' ).fetchone() cost_units = (cost_units[0] or '').strip() if cost_units else '' @@ -90,7 +110,9 @@ def unit_conversion(con: sqlite3.Connection, commodity: str, cost_scale: float | if cost_scale is None: cost_scale = COST_TO_USD.get(cost_units.upper().replace(' ', '')) if cost_scale is None: - sys.exit(f'Unknown cost unit {cost_units!r}; pass --cost-scale (dollars per cost unit).') + sys.exit( + f'Unknown cost unit {cost_units!r}; pass --cost-scale (dollars per cost unit).' + ) note = f'{cost_units or "?"}/{c_units or "?"} -> $/MWh (x{cost_scale / energy_to_mwh:g})' return cost_scale / energy_to_mwh, note @@ -98,10 +120,12 @@ def unit_conversion(con: sqlite3.Connection, commodity: str, cost_scale: float | def window_year(dual_scenario: str, scenario: str) -> int: """Myopic runs save each window's duals as '-'; perfect foresight uses '' alone (treated as one window starting before every period).""" - return -1 if dual_scenario == scenario else int(dual_scenario[len(scenario) + 1:]) + return -1 if dual_scenario == scenario else int(dual_scenario[len(scenario) + 1 :]) -def load_duals(con: sqlite3.Connection, scenario: str, commodities: list[str]) -> tuple[pd.DataFrame, pd.DataFrame]: +def load_duals( + con: sqlite3.Connection, scenario: str, commodities: list[str] +) -> tuple[pd.DataFrame, pd.DataFrame]: """Return (per-slice duals for the requested commodities, all commodity-balance duals). For myopic runs, each period's dual is taken from the latest window whose base year is @@ -109,7 +133,7 @@ def load_duals(con: sqlite3.Connection, scenario: str, commodities: list[str]) - clear and rewrite results from their base year on, but output_dual_variable is never cleared, so look-ahead periods of earlier windows are still in it).""" raw = con.execute( - "SELECT scenario, constraint_name, dual FROM output_dual_variable " + 'SELECT scenario, constraint_name, dual FROM output_dual_variable ' "WHERE (scenario = ? OR scenario GLOB ? || '-[0-9][0-9][0-9][0-9]') " "AND (constraint_name LIKE 'commodity_balance_constraint[%' " "OR constraint_name LIKE 'annual_commodity_balance_constraint[%')", @@ -134,8 +158,11 @@ def load_duals(con: sqlite3.Connection, scenario: str, commodities: list[str]) - annual = pd.DataFrame(annual_rows, columns=['window', 'region', 'period', 'commodity', 'dual']) hit = annual[annual.commodity.isin(commodities)] if not hit.empty: - print(f'NOTE: {sorted(hit.commodity.unique())} are annual commodities; ' - 'their duals are annual-balance duals and are not per-slice.', file=sys.stderr) + print( + f'NOTE: {sorted(hit.commodity.unique())} are annual commodities; ' + 'their duals are annual-balance duals and are not per-slice.', + file=sys.stderr, + ) return all_slice[all_slice.commodity.isin(commodities)].copy(), all_slice @@ -143,17 +170,19 @@ def load_weights(con: sqlite3.Connection, scenario: str, commodities: list[str]) """Energy consumed from the commodity node in each region/slice (excludes storage charging and exports, which are recorded under 'A-B' exchange regions).""" q = ( - "SELECT f.region, f.period, f.season, f.tod, f.input_comm AS commodity, SUM(f.flow) AS load " - "FROM output_flow_in f JOIN technology t ON t.tech = f.tech " - f"WHERE f.scenario = ? AND f.input_comm IN ({','.join('?' * len(commodities))}) " + 'SELECT f.region, f.period, f.season, f.tod, f.input_comm AS commodity, SUM(f.flow) AS load ' + 'FROM output_flow_in f JOIN technology t ON t.tech = f.tech ' + f'WHERE f.scenario = ? AND f.input_comm IN ({",".join("?" * len(commodities))}) ' "AND t.flag NOT LIKE '%s%' AND f.input_comm != f.output_comm " - "GROUP BY f.region, f.period, f.season, f.tod, f.input_comm" + 'GROUP BY f.region, f.period, f.season, f.tod, f.input_comm' ) return pd.read_sql_query(q, con, params=[scenario, *commodities]) def segment_fractions(con: sqlite3.Connection) -> pd.DataFrame: - seasons = pd.read_sql_query('SELECT season, segment_fraction AS sf_season FROM time_season', con) + seasons = pd.read_sql_query( + 'SELECT season, segment_fraction AS sf_season FROM time_season', con + ) tods = pd.read_sql_query('SELECT tod, hours FROM time_of_day', con) tods['tod_frac'] = tods.hours / tods.hours.sum() seg = seasons.merge(tods[['tod', 'tod_frac']], how='cross') @@ -167,14 +196,24 @@ def weighted(g: pd.DataFrame, w: str) -> float: def main() -> None: - ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + ap = argparse.ArgumentParser( + description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter + ) ap.add_argument('database', type=Path) ap.add_argument('--scenario', help='scenario name (default: the only one with duals)') - ap.add_argument('--commodity', nargs='+', default=['ELC'], help='commodities to price (default ELC)') + ap.add_argument( + '--commodity', nargs='+', default=['ELC'], help='commodities to price (default ELC)' + ) ap.add_argument('--base-year', type=int, help='override P0 (default: first future period)') - ap.add_argument('--cost-scale', type=float, help='dollars per objective cost unit (default from units)') - ap.add_argument('--energy-to-mwh', type=float, help='MWh per commodity unit (default from units)') - ap.add_argument('--flip-sign', action='store_true', help='negate duals (only if your solver flips them)') + ap.add_argument( + '--cost-scale', type=float, help='dollars per objective cost unit (default from units)' + ) + ap.add_argument( + '--energy-to-mwh', type=float, help='MWh per commodity unit (default from units)' + ) + ap.add_argument( + '--flip-sign', action='store_true', help='negate duals (only if your solver flips them)' + ) ap.add_argument('--spike', type=float, default=500.0, help='flag slice prices above this $/MWh') ap.add_argument('--out-dir', type=Path, default=Path('.'), help='where to write CSVs') args = ap.parse_args() @@ -202,13 +241,18 @@ def main() -> None: dfs, gdr, p0 = discount_factors(con, args.base_year) duals, all_balance = load_duals(con, scenario, args.commodity) if duals.empty: - sys.exit(f'No commodity_balance_constraint duals for {args.commodity} in scenario {scenario}.') + sys.exit( + f'No commodity_balance_constraint duals for {args.commodity} in scenario {scenario}.' + ) # Sign sanity check: across all commodity balances, nonzero duals should be mostly positive. nz = all_balance.dual[all_balance.dual.abs() > 1e-9] if len(nz) and (nz < 0).mean() > 0.5 and not args.flip_sign: - print(f'WARNING: {100 * (nz < 0).mean():.0f}% of nonzero commodity-balance duals are negative; ' - 'your solver may report the opposite sign convention (see --flip-sign).', file=sys.stderr) + print( + f'WARNING: {100 * (nz < 0).mean():.0f}% of nonzero commodity-balance duals are negative; ' + 'your solver may report the opposite sign convention (see --flip-sign).', + file=sys.stderr, + ) duals = duals.merge(dfs, on='period', how='left') if duals.discount_factor.isna().any(): @@ -233,26 +277,44 @@ def main() -> None: # annual summaries rows = [] for (c, r, p), g in prices.groupby(['commodity', 'region', 'period']): - rows.append({ - 'commodity': c, 'region': r, 'period': p, - 'load_weighted': weighted(g, 'load'), - 'time_weighted': weighted(g, 'segment_fraction'), - 'min': g.price.min(), 'median': g.price.median(), 'max': g.price.max(), - 'load_PJ': g['load'].sum(), - 'n_slices': len(g), - 'n_negative': int((g.price < -1e-6).sum()), - 'n_zero': int((g.price.abs() <= 1e-6).sum()), - 'n_zero_with_load': int(((g.price.abs() <= 1e-6) & (g['load'] > 0)).sum()), - 'n_spike': int((g.price > args.spike).sum()), - }) + rows.append( + { + 'commodity': c, + 'region': r, + 'period': p, + 'load_weighted': weighted(g, 'load'), + 'time_weighted': weighted(g, 'segment_fraction'), + 'min': g.price.min(), + 'median': g.price.median(), + 'max': g.price.max(), + 'load_PJ': g['load'].sum(), + 'n_slices': len(g), + 'n_negative': int((g.price < -1e-6).sum()), + 'n_zero': int((g.price.abs() <= 1e-6).sum()), + 'n_zero_with_load': int(((g.price.abs() <= 1e-6) & (g['load'] > 0)).sum()), + 'n_spike': int((g.price > args.spike).sum()), + } + ) annual = pd.DataFrame(rows) args.out_dir.mkdir(parents=True, exist_ok=True) stem = f'{args.database.stem}_{scenario}' - slice_cols = ['commodity', 'region', 'period', 'season', 'tod', 'window', 'dual', - 'discount_factor', 'price', 'load', 'segment_fraction'] - prices.sort_values(keys)[slice_cols].rename(columns={'price': 'price_usd_per_mwh', 'load': 'load_PJ'}) \ - .to_csv(args.out_dir / f'{stem}_prices_by_slice.csv', index=False) + slice_cols = [ + 'commodity', + 'region', + 'period', + 'season', + 'tod', + 'window', + 'dual', + 'discount_factor', + 'price', + 'load', + 'segment_fraction', + ] + prices.sort_values(keys)[slice_cols].rename( + columns={'price': 'price_usd_per_mwh', 'load': 'load_PJ'} + ).to_csv(args.out_dir / f'{stem}_prices_by_slice.csv', index=False) annual.to_csv(args.out_dir / f'{stem}_prices_annual.csv', index=False) pd.set_option('display.width', 200) @@ -260,18 +322,24 @@ def main() -> None: print(dfs.to_string(index=False)) for c, g in annual.groupby('commodity'): print(f'\n{c}: load-weighted average price ($/MWh), region x period') - print(g.pivot(index='region', columns='period', values='load_weighted').round(1).to_string()) + print( + g.pivot(index='region', columns='period', values='load_weighted').round(1).to_string() + ) # diagnostics print('\nDiagnostics') neg = prices[prices.price < -1e-6] - print(f' negative slice prices: {len(neg)} of {len(prices)}' - + (f' (min {neg.price.min():.1f} $/MWh)' if len(neg) else '')) + print( + f' negative slice prices: {len(neg)} of {len(prices)}' + + (f' (min {neg.price.min():.1f} $/MWh)' if len(neg) else '') + ) zero_load = prices[(prices.price.abs() <= 1e-6) & (prices['load'] > 0)] print(f' zero prices in slices with load: {len(zero_load)}') spikes = prices[prices.price > args.spike] - print(f' slice prices > {args.spike:g} $/MWh: {len(spikes)}' - + (f' (max {spikes.price.max():.0f})' if len(spikes) else '')) + print( + f' slice prices > {args.spike:g} $/MWh: {len(spikes)}' + + (f' (max {spikes.price.max():.0f})' if len(spikes) else '') + ) # Degeneracy hint: duals that jump between adjacent slices with nearly identical load, # or many slices sharing one exact price with a few far-off outliers. ratio = annual.load_weighted / annual.time_weighted @@ -279,7 +347,9 @@ def main() -> None: print(f' region-periods where load- and time-weighted averages differ >1.5x: {len(off)}') no_load = annual[annual.load_PJ <= 0] if len(no_load): - print(f' region-periods with no recorded load (load-weighted average undefined): {len(no_load)}') + print( + f' region-periods with no recorded load (load-weighted average undefined): {len(no_load)}' + ) print(f'\nWrote {args.out_dir / (stem + "_prices_by_slice.csv")}') print(f'Wrote {args.out_dir / (stem + "_prices_annual.csv")}') diff --git a/scripts/reserve_dual_toy.py b/scripts/reserve_dual_toy.py index 9aa5900a1..f677b79d9 100644 --- a/scripts/reserve_dual_toy.py +++ b/scripts/reserve_dual_toy.py @@ -1,38 +1,68 @@ """Toy capacity-expansion LP with Temoa's discounting, to show where a 2025 turbine's capital shows up in later reserve-margin duals. One region, one peak slice, gas turbines only.""" -import highspy, numpy as np + +import highspy +import numpy as np g, p0, periods, p_end = 0.02, 2020, [2025, 2030, 2035, 2040, 2045, 2050], 2055 -A, F = 26.77, 20.527 # yearly capital-equivalent and fixed O&M, M$/GW-yr -cc, prm = 0.9, 0.35 # capacity credit, reserve margin +A, F = 26.77, 20.527 # yearly capital-equivalent and fixed O&M, M$/GW-yr +cc, prm = 0.9, 0.35 # capacity credit, reserve margin pa = lambda n: ((1 + g) ** n - 1) / (g * (1 + g) ** n) -DF = {p: pa(5) / (1 + g) ** (p - p0) for p in periods} # same factor as my script +DF = {p: pa(5) / (1 + g) ** (p - p0) for p in periods} # same factor as my script # Temoa: loan payments from vintage v to horizon end, FOM every active period -> both sum DF_p, p>=v cost = {v: (A + F) * sum(DF[p] for p in periods if p >= v) for v in periods} + def solve(load, existing, allow_build): - h = highspy.Highs(); h.setOptionValue('output_flag', False) + h = highspy.Highs() + h.setOptionValue('output_flag', False) n = len(periods) - for v in periods: # N_v >= 0 - h.addVar(0, highspy.kHighsInf if allow_build(v) else 0); h.changeColCost(periods.index(v), cost[v]) - for i, p in enumerate(periods): # cc*(E + sum_{v<=p} N_v) >= (1+prm)*load_p + for v in periods: # N_v >= 0 + h.addVar(0, highspy.kHighsInf if allow_build(v) else 0) + h.changeColCost(periods.index(v), cost[v]) + for i, p in enumerate(periods): # cc*(E + sum_{v<=p} N_v) >= (1+prm)*load_p idx = [j for j, v in enumerate(periods) if v <= p] - h.addRow((1 + prm) * load[p] - cc * existing, highspy.kHighsInf, len(idx), np.array(idx, dtype=np.int32), np.full(len(idx), cc)) + h.addRow( + (1 + prm) * load[p] - cc * existing, + highspy.kHighsInf, + len(idx), + np.array(idx, dtype=np.int32), + np.full(len(idx), cc), + ) h.run() - sol = h.getSolution(); N = sol.col_value; rho = sol.row_dual - print(' period build_GW dual(disc.) dual/DF=M$/GW-yr of requirement $/MWh adder in 57-h slice') + sol = h.getSolution() + N = sol.col_value + rho = sol.row_dual + print( + ' period build_GW dual(disc.) dual/DF=M$/GW-yr of requirement $/MWh adder in 57-h slice' + ) for i, p in enumerate(periods): - d = (1 + prm) * rho[i] # d obj / d (1 GW more peak load) - print(f' {p} {N[i]:7.2f} {d:9.2f} {d / DF[p]:7.2f} {d / DF[p] * 1e6 / 57000:7.0f}') - print(' check: cost of 2025 turbine =', round(cost[2025] * cc / cc, 2), ' = sum of cc*rho_p over periods it is active:', - round(sum(cc * rho[i] for i in range(len(periods))), 2) if N[0] > 1e-9 else '(not built)') + d = (1 + prm) * rho[i] # d obj / d (1 GW more peak load) + print( + f' {p} {N[i]:7.2f} {d:9.2f} {d / DF[p]:7.2f} {d / DF[p] * 1e6 / 57000:7.0f}' + ) + print( + ' check: cost of 2025 turbine =', + round(cost[2025] * cc / cc, 2), + ' = sum of cc*rho_p over periods it is active:', + round(sum(cc * rho[i] for i in range(len(periods))), 2) if N[0] > 1e-9 else '(not built)', + ) + -print(f'DF: ' + ', '.join(f'{p}:{DF[p]:.3f}' for p in periods)) -print(f'(A+F)={A+F:.2f} M$/GW-yr; with reserve factor (1+prm)/cc={(1+prm)/cc:.2f} -> {(A+F)*(1+prm)/cc:.1f} M$/GW-yr of load\n') +print('DF: ' + ', '.join(f'{p}:{DF[p]:.3f}' for p in periods)) +print( + f'(A+F)={A + F:.2f} M$/GW-yr; with reserve factor (1+prm)/cc={(1 + prm) / cc:.2f} -> {(A + F) * (1 + prm) / cc:.1f} M$/GW-yr of load\n' +) base = 10.0 print('CASE 1: peak load grows 1 GW every period -> a new turbine every period') solve({p: base + i + 1 for i, p in enumerate(periods)}, base * (1 + prm) / cc, lambda v: True) print('\nCASE 2: peak steps up 1 GW in 2025 and then stays flat -> only the 2025 turbine') -solve({p: base + 1 for p in periods}, base * (1 + prm) / cc, lambda v: True) -print('\nCASE 3: 2025 has spare capacity, peak steps up 1 GW in 2030, but 2030+ builds are not allowed') -solve({p: base if p == 2025 else base + 1 for p in periods}, base * (1 + prm) / cc, lambda v: v == 2025) +solve(dict.fromkeys(periods, base + 1), base * (1 + prm) / cc, lambda v: True) +print( + '\nCASE 3: 2025 has spare capacity, peak steps up 1 GW in 2030, but 2030+ builds are not allowed' +) +solve( + {p: base if p == 2025 else base + 1 for p in periods}, + base * (1 + prm) / cc, + lambda v: v == 2025, +) From 2c72c7618c5242a95e21de66341f60151ce801c9 Mon Sep 17 00:00:00 2001 From: Joe DeCarolis Date: Fri, 2 Oct 2026 09:28:18 -0400 Subject: [PATCH 10/10] Remove accidentally-included local scratch files from merge commit These were untracked WIP files in the working directory (electricity price estimation from duals) that got swept in by an overly broad git add, not part of the unstable -> main sync. Includes a pre-commit.ci reformat of those same files from while they were briefly in the PR. --- scripts/electricity_prices.py | 358 --------------- scripts/electricity_prices_notes.md | 665 ---------------------------- scripts/reserve_dual_toy.py | 68 --- 3 files changed, 1091 deletions(-) delete mode 100644 scripts/electricity_prices.py delete mode 100644 scripts/electricity_prices_notes.md delete mode 100644 scripts/reserve_dual_toy.py diff --git a/scripts/electricity_prices.py b/scripts/electricity_prices.py deleted file mode 100644 index ff3e7e65e..000000000 --- a/scripts/electricity_prices.py +++ /dev/null @@ -1,358 +0,0 @@ -""" -Convert Temoa commodity-balance duals into electricity prices ($/MWh). - -What the duals are (verified against temoa/components/commodities.py and costs.py): - -* Electricity commodities (flag 'p', not annual, not waste) are balanced per - (region, period, season, tod) by ``commodity_balance_constraint`` as an - equality: produced == consumed. (Annual commodities use - ``annual_commodity_balance_constraint[r,p,c]``; waste commodities use >=.) -* Flows in a slice are energy over that slice for one representative year - (capacity constraint: FO <= CF * C2A * SEG * CAP), so a dual is the change in - the objective from one extra unit of energy consumed in that slice in every - year of the period. No division by segment fraction is needed. -* The objective (minimize) discounts annual costs in period p by - DF_p = annuity_to_pv(GDR, LEN_p) * (1 + GDR) ** -(p - P0) - = sum_{k=1..LEN_p} (1 + GDR) ** -(p - P0 + k) - i.e. end-of-year convention (no mid-year), P0 = first optimized period - (or myopic_discounting_year), LEN_p = next period - p, and the last - period's length is set by the final 'f' row in time_period. -* With Pyomo's convention (dual = d obj / d rhs for body = produced - consumed), - a positive dual is the marginal cost of serving more load, so no sign flip is - needed. The script checks this and warns if the sample looks flipped. -* Duals are stored raw (discounted objective units per commodity unit) in - ``output_dual_variable(scenario, constraint_name, dual)``. - -Usage: - python scripts/electricity_prices.py temoa_power.sqlite [--commodity ELC ELCP] -""" - -from __future__ import annotations - -import argparse -import re -import sqlite3 -import sys -from pathlib import Path - -import pandas as pd - -# Multipliers to convert a commodity unit to MWh and a cost unit to dollars. -ENERGY_TO_MWH = { - 'PJ': 1e15 / 3.6e9, - 'TJ': 1e12 / 3.6e9, - 'GJ': 1 / 3.6, - 'TWH': 1e6, - 'GWH': 1e3, - 'MWH': 1.0, -} -COST_TO_USD = { - 'MUSD': 1e6, - 'M$': 1e6, - 'MDOLLAR': 1e6, - 'BUSD': 1e9, - 'KUSD': 1e3, - 'USD': 1.0, - '$': 1.0, -} - -NAME_RE = re.compile(r'^(?P\w+)\[(?P.*)\]$') - - -def connect(path: Path) -> sqlite3.Connection: - # immutable=1 lets us read WAL-mode databases without write access to the -shm file; - # mode=ro stops sqlite from creating an empty file when the path is wrong - if not path.is_file() or path.stat().st_size == 0: - sys.exit(f'{path} does not exist or is empty.') - return sqlite3.connect(f'file:{path}?mode=ro&immutable=1', uri=True) - - -def discount_factors( - con: sqlite3.Connection, base_year: int | None -) -> tuple[pd.DataFrame, float, int]: - row = con.execute( - "SELECT value FROM metadata_real WHERE element = 'global_discount_rate'" - ).fetchone() - if row is None: - sys.exit('global_discount_rate missing from metadata_real (Temoa requires it).') - gdr = float(row[0]) - future = [ - r[0] for r in con.execute("SELECT period FROM time_period WHERE flag = 'f' ORDER BY period") - ] - p0 = base_year if base_year is not None else future[0] - rows = [] - for p, p_next in zip(future[:-1], future[1:]): - length = int(p_next - p) - if gdr == 0: - df = float(length) - else: - annuity = ((1 + gdr) ** length - 1) / (gdr * (1 + gdr) ** length) - df = annuity / (1 + gdr) ** int(p - p0) - rows.append({'period': p, 'period_length': length, 'discount_factor': df}) - return pd.DataFrame(rows), gdr, p0 - - -def unit_conversion( - con: sqlite3.Connection, commodity: str, cost_scale: float | None, energy_to_mwh: float | None -) -> tuple[float, str]: - """Return multiplier taking (cost unit / commodity unit) to $/MWh.""" - c_units = con.execute('SELECT units FROM commodity WHERE name = ?', (commodity,)).fetchone() - c_units = (c_units[0] or '').strip() if c_units else '' - cost_units = con.execute( - 'SELECT units FROM output_cost WHERE units IS NOT NULL LIMIT 1' - ).fetchone() - cost_units = (cost_units[0] or '').strip() if cost_units else '' - - if energy_to_mwh is None: - energy_to_mwh = ENERGY_TO_MWH.get(c_units.upper()) - if energy_to_mwh is None: - sys.exit(f'Unknown energy unit {c_units!r} for {commodity}; pass --energy-to-mwh.') - if cost_scale is None: - cost_scale = COST_TO_USD.get(cost_units.upper().replace(' ', '')) - if cost_scale is None: - sys.exit( - f'Unknown cost unit {cost_units!r}; pass --cost-scale (dollars per cost unit).' - ) - note = f'{cost_units or "?"}/{c_units or "?"} -> $/MWh (x{cost_scale / energy_to_mwh:g})' - return cost_scale / energy_to_mwh, note - - -def window_year(dual_scenario: str, scenario: str) -> int: - """Myopic runs save each window's duals as '-'; perfect foresight uses - '' alone (treated as one window starting before every period).""" - return -1 if dual_scenario == scenario else int(dual_scenario[len(scenario) + 1 :]) - - -def load_duals( - con: sqlite3.Connection, scenario: str, commodities: list[str] -) -> tuple[pd.DataFrame, pd.DataFrame]: - """Return (per-slice duals for the requested commodities, all commodity-balance duals). - - For myopic runs, each period's dual is taken from the latest window whose base year is - <= that period: the window whose solution was kept in the output tables (later windows - clear and rewrite results from their base year on, but output_dual_variable is never - cleared, so look-ahead periods of earlier windows are still in it).""" - raw = con.execute( - 'SELECT scenario, constraint_name, dual FROM output_dual_variable ' - "WHERE (scenario = ? OR scenario GLOB ? || '-[0-9][0-9][0-9][0-9]') " - "AND (constraint_name LIKE 'commodity_balance_constraint[%' " - "OR constraint_name LIKE 'annual_commodity_balance_constraint[%')", - (scenario, scenario), - ).fetchall() - slice_rows, annual_rows = [], [] - for dual_scenario, name, dual in raw: - w = window_year(dual_scenario, scenario) - m = NAME_RE.match(name) - idx = m['idx'].split(',') - if m['con'] == 'commodity_balance_constraint': - r, p, s, d, c = idx - if w <= int(p): - slice_rows.append((w, r, int(p), s, d, c, dual)) - else: - r, p, c = idx - if w <= int(p): - annual_rows.append((w, r, int(p), c, dual)) - cols = ['window', 'region', 'period', 'season', 'tod', 'commodity', 'dual'] - all_slice = pd.DataFrame(slice_rows, columns=cols) - all_slice = all_slice.sort_values('window').drop_duplicates(cols[1:6], keep='last') - annual = pd.DataFrame(annual_rows, columns=['window', 'region', 'period', 'commodity', 'dual']) - hit = annual[annual.commodity.isin(commodities)] - if not hit.empty: - print( - f'NOTE: {sorted(hit.commodity.unique())} are annual commodities; ' - 'their duals are annual-balance duals and are not per-slice.', - file=sys.stderr, - ) - return all_slice[all_slice.commodity.isin(commodities)].copy(), all_slice - - -def load_weights(con: sqlite3.Connection, scenario: str, commodities: list[str]) -> pd.DataFrame: - """Energy consumed from the commodity node in each region/slice (excludes storage charging - and exports, which are recorded under 'A-B' exchange regions).""" - q = ( - 'SELECT f.region, f.period, f.season, f.tod, f.input_comm AS commodity, SUM(f.flow) AS load ' - 'FROM output_flow_in f JOIN technology t ON t.tech = f.tech ' - f'WHERE f.scenario = ? AND f.input_comm IN ({",".join("?" * len(commodities))}) ' - "AND t.flag NOT LIKE '%s%' AND f.input_comm != f.output_comm " - 'GROUP BY f.region, f.period, f.season, f.tod, f.input_comm' - ) - return pd.read_sql_query(q, con, params=[scenario, *commodities]) - - -def segment_fractions(con: sqlite3.Connection) -> pd.DataFrame: - seasons = pd.read_sql_query( - 'SELECT season, segment_fraction AS sf_season FROM time_season', con - ) - tods = pd.read_sql_query('SELECT tod, hours FROM time_of_day', con) - tods['tod_frac'] = tods.hours / tods.hours.sum() - seg = seasons.merge(tods[['tod', 'tod_frac']], how='cross') - seg['segment_fraction'] = seg.sf_season * seg.tod_frac - return seg[['season', 'tod', 'segment_fraction']] - - -def weighted(g: pd.DataFrame, w: str) -> float: - tot = g[w].sum() - return float((g.price * g[w]).sum() / tot) if tot > 0 else float('nan') - - -def main() -> None: - ap = argparse.ArgumentParser( - description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter - ) - ap.add_argument('database', type=Path) - ap.add_argument('--scenario', help='scenario name (default: the only one with duals)') - ap.add_argument( - '--commodity', nargs='+', default=['ELC'], help='commodities to price (default ELC)' - ) - ap.add_argument('--base-year', type=int, help='override P0 (default: first future period)') - ap.add_argument( - '--cost-scale', type=float, help='dollars per objective cost unit (default from units)' - ) - ap.add_argument( - '--energy-to-mwh', type=float, help='MWh per commodity unit (default from units)' - ) - ap.add_argument( - '--flip-sign', action='store_true', help='negate duals (only if your solver flips them)' - ) - ap.add_argument('--spike', type=float, default=500.0, help='flag slice prices above this $/MWh') - ap.add_argument('--out-dir', type=Path, default=Path('.'), help='where to write CSVs') - args = ap.parse_args() - - con = connect(args.database) - - scenarios = [r[0] for r in con.execute('SELECT DISTINCT scenario FROM output_dual_variable')] - if not scenarios: - flows = [r[0] for r in con.execute('SELECT DISTINCT scenario FROM output_cost')] - found = f'results exist for scenarios {flows}' if flows else 'no model results at all' - sys.exit( - f'output_dual_variable is empty in {args.database.name} ({found}). ' - 'Run Temoa with save_duals = true and a solver that returns duals ' - '(e.g. gurobi; not appsi_highs), then run this script again.' - ) - # strip myopic window suffixes ('-2020', ...) to get the scenario the flows are stored under - bases = sorted({re.sub(r'-\d{4}$', '', s) for s in scenarios}) - scenario = args.scenario or (bases[0] if len(bases) == 1 else None) - if scenario is None: - sys.exit(f'Multiple scenarios with duals: {bases}. Pass --scenario.') - windows = sorted(s for s in scenarios if s != scenario and s.startswith(scenario + '-')) - if windows: - print(f'Myopic run: stitching duals from windows {windows}', file=sys.stderr) - - dfs, gdr, p0 = discount_factors(con, args.base_year) - duals, all_balance = load_duals(con, scenario, args.commodity) - if duals.empty: - sys.exit( - f'No commodity_balance_constraint duals for {args.commodity} in scenario {scenario}.' - ) - - # Sign sanity check: across all commodity balances, nonzero duals should be mostly positive. - nz = all_balance.dual[all_balance.dual.abs() > 1e-9] - if len(nz) and (nz < 0).mean() > 0.5 and not args.flip_sign: - print( - f'WARNING: {100 * (nz < 0).mean():.0f}% of nonzero commodity-balance duals are negative; ' - 'your solver may report the opposite sign convention (see --flip-sign).', - file=sys.stderr, - ) - - duals = duals.merge(dfs, on='period', how='left') - if duals.discount_factor.isna().any(): - sys.exit('Duals found for periods without a discount factor; check time_period.') - sign = -1.0 if args.flip_sign else 1.0 - - pieces = [] - for c, g in duals.groupby('commodity'): - mult, note = unit_conversion(con, c, args.cost_scale, args.energy_to_mwh) - g = g.copy() - g['price'] = sign * g.dual / g.discount_factor * mult - pieces.append(g) - print(f'{c}: units {note}') - prices = pd.concat(pieces) - - weights = load_weights(con, scenario, args.commodity) - seg = segment_fractions(con) - keys = ['region', 'period', 'season', 'tod', 'commodity'] - prices = prices.merge(weights, on=keys, how='left').merge(seg, on=['season', 'tod'], how='left') - prices['load'] = prices['load'].fillna(0.0) - - # annual summaries - rows = [] - for (c, r, p), g in prices.groupby(['commodity', 'region', 'period']): - rows.append( - { - 'commodity': c, - 'region': r, - 'period': p, - 'load_weighted': weighted(g, 'load'), - 'time_weighted': weighted(g, 'segment_fraction'), - 'min': g.price.min(), - 'median': g.price.median(), - 'max': g.price.max(), - 'load_PJ': g['load'].sum(), - 'n_slices': len(g), - 'n_negative': int((g.price < -1e-6).sum()), - 'n_zero': int((g.price.abs() <= 1e-6).sum()), - 'n_zero_with_load': int(((g.price.abs() <= 1e-6) & (g['load'] > 0)).sum()), - 'n_spike': int((g.price > args.spike).sum()), - } - ) - annual = pd.DataFrame(rows) - - args.out_dir.mkdir(parents=True, exist_ok=True) - stem = f'{args.database.stem}_{scenario}' - slice_cols = [ - 'commodity', - 'region', - 'period', - 'season', - 'tod', - 'window', - 'dual', - 'discount_factor', - 'price', - 'load', - 'segment_fraction', - ] - prices.sort_values(keys)[slice_cols].rename( - columns={'price': 'price_usd_per_mwh', 'load': 'load_PJ'} - ).to_csv(args.out_dir / f'{stem}_prices_by_slice.csv', index=False) - annual.to_csv(args.out_dir / f'{stem}_prices_annual.csv', index=False) - - pd.set_option('display.width', 200) - print(f'\nScenario {scenario}: GDR={gdr}, P0={p0}') - print(dfs.to_string(index=False)) - for c, g in annual.groupby('commodity'): - print(f'\n{c}: load-weighted average price ($/MWh), region x period') - print( - g.pivot(index='region', columns='period', values='load_weighted').round(1).to_string() - ) - - # diagnostics - print('\nDiagnostics') - neg = prices[prices.price < -1e-6] - print( - f' negative slice prices: {len(neg)} of {len(prices)}' - + (f' (min {neg.price.min():.1f} $/MWh)' if len(neg) else '') - ) - zero_load = prices[(prices.price.abs() <= 1e-6) & (prices['load'] > 0)] - print(f' zero prices in slices with load: {len(zero_load)}') - spikes = prices[prices.price > args.spike] - print( - f' slice prices > {args.spike:g} $/MWh: {len(spikes)}' - + (f' (max {spikes.price.max():.0f})' if len(spikes) else '') - ) - # Degeneracy hint: duals that jump between adjacent slices with nearly identical load, - # or many slices sharing one exact price with a few far-off outliers. - ratio = annual.load_weighted / annual.time_weighted - off = annual[(ratio > 1.5) | (ratio < 0.67)] - print(f' region-periods where load- and time-weighted averages differ >1.5x: {len(off)}') - no_load = annual[annual.load_PJ <= 0] - if len(no_load): - print( - f' region-periods with no recorded load (load-weighted average undefined): {len(no_load)}' - ) - print(f'\nWrote {args.out_dir / (stem + "_prices_by_slice.csv")}') - print(f'Wrote {args.out_dir / (stem + "_prices_annual.csv")}') - - -if __name__ == '__main__': - main() diff --git a/scripts/electricity_prices_notes.md b/scripts/electricity_prices_notes.md deleted file mode 100644 index f2f99cbc8..000000000 --- a/scripts/electricity_prices_notes.md +++ /dev/null @@ -1,665 +0,0 @@ -# Electricity Shadow Prices from Temoa: What the Duals Mean and How to Use Them - -## Purpose - -Temoa is a linear-programming capacity-expansion model. When it is solved, every constraint -comes with a *dual value*, also called a shadow price. The dual on the electricity balance -constraint is the natural place to look for the marginal cost of electricity supply: how much -total system cost would rise if one more unit of electricity had to be supplied in a given -region, period and time slice. - -The raw duals Temoa writes out are not directly usable as prices. They are expressed in the -units of Temoa's objective function, which discounts every cost back to a base year and adds -up costs over all years of each multi-year period. They also depend on the time-slice -structure, on whether the model was run with perfect foresight or myopically, and on how -capital costs are spread over time in the objective. - -This document explains: - -1. what the raw electricity duals mean in Temoa; -2. what costs they include, in perfect-foresight and myopic runs, using a new gas turbine as a - worked example; -3. how the companion script `scripts/electricity_prices.py` turns them into estimated marginal - prices in $/MWh by region, period and time slice; -4. the caveats to keep in mind when reading those prices; -5. how those marginal prices compare with the average cost of electricity computed from - Temoa's cost outputs. - -Everything here was checked against the Temoa source code (repo commit `fc10702d`) and the -`temoa_power.sqlite` input database. Numerical claims were tested by solving small cases on -copies of Temoa's Utopia tutorial database. A small stand-alone example model is in -`scripts/reserve_dual_toy.py`. - ---- - -## 1. What the raw electricity duals mean - -### 1.1 A quick refresher on duals - -In a cost-minimizing linear program, the dual of a constraint tells you how much the optimal -objective would change if that constraint's right-hand side were loosened or tightened by one -unit, with everything else free to adjust. For a supply–demand balance, "tighten by one unit" -means "serve one more unit of demand." The dual is therefore the system's marginal cost of -supplying that commodity, at that place and time. - -Two features of Temoa determine how to read that number: which constraint is the electricity -balance, and how the objective function measures cost. - -### 1.2 Which constraint balances electricity - -In `temoa_power`, electricity passes through several commodities. Generators produce `ELCP`, -which can be thought of as busbar electricity. A loss technology, `E_ELCTDLOSS` (efficiency -0.953), converts `ELCP` to `ELC`. `ELC` is the regional grid node. Storage (`E_Batt`, -`E_Batt8hr`, `E_HYDPS_R`), interregional transmission (`E_TRANS_N`, `E_TRANS_R`), and the -technologies that deliver electricity to the residential, commercial and transport sectors all -connect at `ELC`. `ELC` is therefore the most natural place to measure the price of grid -electricity, and it is the script's default. - -Both `ELC` and `ELCP` are ordinary physical commodities in Temoa's classification: they are -neither "annual" nor "waste" commodities. Temoa therefore balances them with -**`commodity_balance_constraint[r, p, s, d, c]`** (`temoa/components/commodities.py`). This is -an equality, `produced == consumed`, written separately for every region r, period p, season s -and time of day d. Other balance-type constraints apply to different kinds of commodity: - -- `demand_constraint` applies only to end-use demand commodities. -- `annual_commodity_balance_constraint` applies only to commodities balanced over a whole - year. -- Waste commodities use an inequality instead of an equality. - -None of these applies to grid electricity. - -### 1.3 What a time slice represents - -`temoa_power` has two seasons, covering 62.4% and 37.6% of the year. Each season has 96 -time-of-day entries, and each entry is a single representative hour. That gives 192 slices per -region and period. (The `hours` column in the `time_of_day` table is 0.25 for every entry, but -Temoa uses it only as a relative weight: each hour receives an equal 1/96 share of its season.) -Because 192 representative hours stand in for all 8,760 hours of the year, each slice is -weighted to represent many real hours: about 57 hours per year for each `S1` hour, and about -34 for each `S2` hour. - -Temoa's flow variables measure **energy over the whole slice for one representative year**, -not power. This follows from the capacity constraint, `output ≤ CF × C2A × SEG × capacity`, -where SEG is the slice's fraction of the year. The electricity balance is therefore written -in PJ per slice-year. Its dual is the cost of consuming one extra PJ in that slice, **in every -year of the period**. Because the flows already measure energy over the slice, there is no -need to divide by the slice's length to get a per-energy price. - -### 1.4 Sign - -Temoa's modeling layer (Pyomo) stores the balance as `produced − consumed = 0` and reports each -dual as the change in the objective per unit increase in the right-hand side. Raising the -right-hand side means production must exceed consumption, which is the same as adding load. -So for this cost-minimizing model, **a positive dual means a positive marginal cost**, and no -sign flip is needed. - -This was confirmed on the Utopia test database solved with Gurobi. Fuel balance duals came out -positive, and for a fuel whose marginal source is an import with a fixed cost of 2.0 M$/PJ, the -dual divided by the discount factor (§1.6) came back as exactly 2.0000 in every period. - -Because the balance is an equality, the dual can still be negative in some slices. Section 4 -explains when that happens. - -### 1.5 Where duals are stored, and in what units - -When a run is made with `save_duals = true`, Temoa writes duals to the output table -`output_dual_variable(scenario, constraint_name, dual)`. They are raw solver values, keyed by -strings such as `commodity_balance_constraint[CA,2030,S1,H5,ELC]`. - -The units are those of the objective function divided by those of the commodity. For -`temoa_power` that is **millions of dollars, discounted to 2020, per PJ**. - -Perfect-foresight runs store duals under the scenario name. Myopic runs store them separately -for each window, under `-`; see §3. - -### 1.6 How Temoa's objective discounts costs - -Temoa's objective (`temoa/components/costs.py`) is the total discounted cost of the energy -system. Costs that recur every year, such as fuel, variable and fixed O&M, and emission -charges, are calculated for a single representative year of each period. They are then -multiplied by a **period discount factor**: - -``` -DF_p = (sum over the years of period p of each year's discount factor) - = Σ_{k = 1 .. LEN_p} (1 + r)^−(p − P0 + k) -``` - -where: - -- **r** is the global discount rate, 2% in `temoa_power`. -- **P0** is the base year, 2020: the first optimized period. In myopic runs Temoa uses the - first future period for every window, so P0 is also 2020 there. -- **LEN_p** is the period length: the gap to the next period. Every period from 2020 to 2050 is - five years long. The final 2055 entry only marks the end of the horizon and sets the length - of the last period. - -The convention is **end of year**: the first year of each period is already discounted by one -full year, with no mid-year adjustment. For `temoa_power`: - -| Period | 2020 | 2025 | 2030 | 2035 | 2040 | 2045 | 2050 | -|---|---|---|---|---|---|---|---| -| DF_p | 4.713 | 4.269 | 3.867 | 3.502 | 3.172 | 2.873 | 2.602 | - -Because the dual is the change in this discounted, multi-year objective, it carries the same -factor. **Dividing a dual by DF_p converts it back to an ordinary cost per unit in a single -year of period p, in undiscounted dollars.** This is the key step in turning duals into -prices. - -Capital costs are handled differently. The overnight cost of new capacity is first converted -into loan payments at a technology-specific loan rate. That stream of payments is converted to -a present value and spread evenly over the plant's lifetime at the global discount rate. -**Only the payments that fall before the end of the model horizon are charged.** In perfect -foresight the horizon ends in 2055; in myopic runs it is the end of the current window. This -is the standard way to avoid charging the model for years it cannot see, and it matters for -the discussion of myopic runs below. - ---- - -## 2. What the duals include: a gas turbine example - -### 2.1 The electricity dual already includes every other constraint - -It is natural to ask whether the electricity balance dual captures everything, or whether you -also need the duals on the reserve-margin constraint, emission caps and so on. **The balance -dual already captures everything.** When the model must supply a little more electricity in a -slice, the solver re-optimizes everything to meet it: dispatch, storage, imports, new -investment, and every constraint those decisions touch. The dual is the total cost of that -cheapest response, including all knock-on effects. - -The optimality conditions of the linear program show how the pieces fit together. Suppose a -gas turbine supplies the extra electricity in a peak slice. At the optimum: - -``` -λ_ELC = DF_p × VOM + λ_gas / η + μ_cap + (1 + PRM) × ρ_reserve + EF × π_CO2 + … -``` - -where: - -- **λ_ELC** is the raw dual of the `ELC` balance in that slice. -- **λ_gas** is the dual of the balance for the gas the turbine burns (`E_NGA`, also balanced - per slice): the model's own marginal cost of gas in that slice. It is divided by the - efficiency η (0.351) because one PJ of electricity needs 2.85 PJ of gas. -- **μ_cap** is the dual of the turbine's capacity constraint, which is nonzero only when the - turbine is running flat out. -- **ρ_reserve** is the dual of the planning reserve-margin constraint in that slice. The - reserve requirement is 1 + PRM times the output of reserve-eligible plants, so one more PJ - from the turbine raises the requirement by 1.35 PJ with a 35% margin. -- **π_CO2** is the dual of any binding emissions limit, multiplied by the turbine's emissions - per PJ (EF). - -Every term is in the same discounted units. Only VOM shows an explicit `DF_p`, because the -objective's cost coefficients are discounted while the duals come out discounted -automatically. - -The practical conclusion is that **you should not add the reserve-margin or emission duals to -the electricity dual; that would count them twice.** Those other duals are useful only for -breaking the price into its parts (energy, capacity, carbon, renewable mandates), for example -to explain why a slice is expensive. - -### 2.2 Which of the turbine's costs end up in the price - -Suppose a peak slice needs new capacity, and the cheapest option is a new gas turbine -(`E_NGAACT_N`). The inputs used here are for Texas, 2030 vintage: - -- overnight cost 733 M$/GW; -- loan rate 3.9% over a 15-year loan; plant lifetime 50 years; -- fixed O&M 20.5 M$/GW-yr; variable O&M 1.35 M$/PJ (about 4.85 $/MWh); -- efficiency 0.351; -- capacity credit 0.9; -- a 35% planning reserve margin. - -The loan rate, loan life and lifetime were read from other regions' rows and assumed to apply -to Texas. - -After dividing by the discount factor, the price in that slice breaks down as: - -``` -price = running cost + (yearly capital payment + fixed O&M) × k / (slice-years of energy per GW) -``` - -**Running cost.** This is the variable O&M, plus fuel at the model's own gas price divided by -the efficiency, plus any emission charges. - -**Capital, as a yearly payment.** The overnight cost never enters the price directly. Temoa -converts it into an even yearly payment: here, about **26.8 M$/GW per year**. In a myopic run -with one-period windows, the objective charges exactly five years of these payments, and -dividing by the discount factor recovers one year's payment. Adding fixed O&M gives about -**47.3 M$/GW-yr**. - -**Spreading it over very few hours.** A GW of capacity can deliver only about 57 GWh per year in -an `S1` slice, because the slice covers only about 57 hours. One extra PJ in that slice -therefore needs about 4.9 GW of extra capacity, and the yearly capacity cost is spread over a -small amount of energy: - -- **If the turbine's own capacity limit is what binds (k = 1),** the capacity cost adds about - **830 $/MWh** to the price in that slice. -- **If the reserve margin binds first,** which is likely with a 35% margin, each extra unit of - output needs 1.35 units of reserve, and each GW counts at its 0.9 capacity credit, so - k ≈ 1.35 / 0.9 = 1.5. The capacity cost then adds about **1,245 $/MWh**. In this case the - capacity cost reaches the electricity price through the reserve dual, not the turbine's - capacity dual. (Under Temoa's "dynamic" reserve method, the turbine's availability in that - slice replaces the capacity credit.) - -For comparison, spreading the full overnight cost over one year of that slice would give about -12,900 $/MWh. The dual never does that, because it carries yearly payments, not the lump sum. - -These figures are approximate: they assume a single binding slice and the Texas 2030 inputs. -They do, however, explain why peak-slice prices in the thousands of dollars per MWh are normal -in a model like this. - -### 2.3 What about capacity built in earlier periods? - -A related question is whether the capital cost of plants built in earlier periods should ever -appear in a later period's price. The answer depends **not on when the plant was built, but on -whether its size was a decision in the same model solve that produced the dual.** - -**Capacity that is fixed in the solve never contributes its capital to the price.** This -always includes capacity that existed before the model's first period. In myopic runs it also -includes everything built in earlier windows: Temoa loads those builds as fixed existing -capacity (`_load_existing_capacity` in `temoa/data_io/hybrid_loader.py`). The remaining -payments on fixed capacity don't depend on any decision in the current solve, so including them -would only add a constant to the objective, and a constant cannot change a dual. Such plants -are paid, if at all, through the gap between the price and their running cost. - -**Capacity that is chosen in the same solve can contribute its capital to later periods' -prices.** In a perfect-foresight run, a turbine built in 2025 is a decision variable when the -2030 prices are computed, because the model decided to build it with 2030 in view. The next -section shows how that works. - -Earlier builds also affect prices indirectly, through how much capacity exists: an overbuilt -system has prices near running cost, while a tight one has large scarcity rents. That effect -works through the amount of capacity, not its cost. - -### 2.4 Perfect foresight in detail: how a 2025 turbine's capital reaches the 2030 price - -In the objective, a GW of turbine built in 2025 is charged one yearly payment (capital plus -fixed O&M) for every year from 2025 to the end of the horizon, each discounted to 2020. That -total can be written as a sum of per-period pieces: - -``` -Cost_2025 = (A + F) × (DF_2025 + DF_2030 + … + DF_2050) = 47.30 × 20.29 ≈ 959 M$ per GW -``` - -Let ρ_p be the reserve-margin dual in period p, and cc the capacity credit. A linear program -only builds something when its cost exactly equals the value it provides at the optimal duals. -A GW built in 2025 supplies cc GW of reserve in every period from 2025 onward, so: - -``` -Cost_v = cc·ρ_v + cc·ρ_(v+5) + … + cc·ρ_2050 for every vintage v that is built -Cost_v ≥ cc·ρ_v + cc·ρ_(v+5) + … + cc·ρ_2050 for every vintage v that is not built -``` - -The first line is how the 2025 turbine's capital reaches the 2030 price: ρ_2030 is one of the -terms that must add up to its cost. The equation does not say how big each term is. That is -determined by the alternatives the model had, expressed through the inequalities. - -The example model in `scripts/reserve_dual_toy.py` makes this concrete. It has one region and -one peak slice, uses Temoa's discounting, and has the turbine's costs above. It was solved with -the HiGHS solver under three load patterns. The table shows the resulting capacity contribution -to the price in the peak slice, in $/MWh: - -| Case | What gets built | 2025 | 2030 | 2035–2050 | -|---|---|---|---|---| -| 1. Peak load grows 1 GW every period | a new turbine every period | 1,245 | 1,245 | 1,245 each | -| 2. Peak steps up in 2025, then stays flat | the 2025 turbine only | 5,914 | 0 | 0 | -| 3. Peak steps up in 2030; no building allowed after 2025 | the 2025 turbine only | 0 | 6,530 | 0 | - -In all three cases, the reserve duals over the turbine's years of service add up exactly to -its cost (959.43 M$/GW). - -**Case 1 is the typical situation.** Because the model builds in both 2025 and 2030, both -conditions hold with equality. Subtracting one from the other shows that the 2025 dual carries -exactly one period's yearly payments, and so does the 2030 dual. Those 2030 payments are the -2025 turbine's payments for 2030–2034, but they are also exactly what a turbine built in 2030 -would cost. Since the two numbers are identical, it is not meaningful to ask whose capital is -in the 2030 price. The price is the avoidable cost of one more GW of reserve in 2030: one more -period of yearly payments. - -**Case 3 shows the 2030 price carrying all of the 2025 turbine's capital.** The turbine is -needed only in 2030, but the model is not allowed to build then. It must build in 2025 and -carry an idle plant for five years, and the whole 959 M$/GW, including the 2025–2029 -payments, lands in the 2030 price, roughly five times case 1. In a real model this happens -only when building later is blocked or more expensive: new-capacity or growth limits, a cost -that rises over time, or the technology no longer being available. That other constraint would -also have its own nonzero dual. - -**Case 2 shows that the split can be arbitrary.** The reserve constraint binds in every -period against the same turbine, so many different splits of its cost across periods are -equally valid. Case 1's even split of one yearly payment per period is one of them. The solver -happened to put the whole cost on 2025; another solver, or a slightly different data set, could -put it on 2030 instead. The total payback is fixed, but individual period prices are not. If -period prices jump around with no corresponding change in costs or load, this is the likely -cause. - -### 2.5 Myopic runs - -In a myopic run, Temoa solves the horizon as a series of shorter windows, each blind to the -periods after it. `temoa_power` is typically run with one-period windows (`view_depth = 1`, -`step_size = 1`). Capacity built in earlier windows enters later windows as fixed existing -capacity, so **in myopic mode the capital cost of earlier builds never enters a later period's -price.** With one-period windows, each period's price includes only running costs plus the -yearly capital and fixed O&M of capacity built in that period. - -Applied to the three cases above: - -- **Case 1:** the 2030 price is still about 1,245 $/MWh, because building a new turbine in - 2030 is the marginal option. -- **Case 2:** the 2030 price is zero, with no ambiguity, because no new capacity is needed. -- **Case 3:** the 2025 window has no reason to build, because it cannot see the 2030 peak. With - 2030 building ruled out, the model would have to meet the 2030 peak another way, or could not - meet it at all. - -Two consequences follow. First, **nothing ensures that earlier investments recover their -capital.** If a later window brings a tighter carbon cap or a cheaper technology, earlier -plants can be left stranded, and prices need not cover their remaining payments. This reflects -the myopic planner's limited foresight; it is not an error. Second, **window boundaries can -shift prices.** New investments are still priced consistently, because a plant built in a -window is charged a full yearly payment for each year it operates in that window. But -constraints that span several periods are cut off at the window edge, such as multi-period -emission budgets, growth limits, or storage carried between periods. That can affect duals, -especially in the last period of each window. - -### 2.6 A reporting gap in myopic runs - -Reviewing the myopic code turned up a separate issue: **some capital costs are missing from the -cost outputs of myopic runs.** This affects reported totals, not the model's decisions or its -duals. - -In each window, Temoa charges a new plant only for the loan payments that fall inside that -window. Later windows charge investment only for their own new builds. As a result, the -payments due after the building window ends are never charged anywhere, even when the plant -keeps operating. The table that records costs (`output_cost`, both the discounted and -undiscounted investment columns) uses the same cut-off. So does the myopic total system cost, -which is simply the sum of that table. - -This was measured on the Utopia database by comparing a myopic run (one-period windows) with a -perfect-foresight run. The investment cost recorded per unit of new capacity was: - -| Vintage | Myopic as a share of perfect foresight | Why | -|---|---|---| -| 1990 | 33% | 10 of 30 in-horizon years charged | -| 2000 | 50% (40-yr plants), 67% (15-yr plants) | 10 of 20, or 10 of 15, years charged | -| 2010 (last window) | 100% | both runs end at the same year | - -Fixed O&M for plants carried into later windows is charged in every period, and variable and -emission costs are unaffected. - -For `temoa_power` with five-year windows, a rough estimate: a 30-year plant built in 2020 -would have only about 21% of its discounted capital recorded, a 2045 plant about 52%, and a -2050 plant all of it. Reported total system costs will be understated accordingly. - -The model's decisions are unaffected, because within each window the cost charged matches the -years the plant operates there, and in later windows the earlier capital is correctly treated -as sunk. The duals are unaffected for the reasons in §2.3. The practical consequence is that -**prices should not be compared with `output_cost` capital** to judge whether plants recover -their costs; for that, capital would need to be rebuilt over the full horizon from the built -capacity. - ---- - -## 3. How the script turns duals into prices - -`scripts/electricity_prices.py` reads a solved Temoa database and produces estimated marginal -prices for electricity. For each slice dual on the chosen commodity (`ELC` by default): - -``` -price ($/MWh) = dual / DF_p × (dollars per cost unit) / (MWh per commodity unit) -``` - -**Undoing the discounting.** The script reads the discount rate and periods from the database -and rebuilds DF_p exactly as Temoa's objective does (§1.6). Dividing by it turns each dual into -the undiscounted cost of one more unit in one year of the period. That is the quantity that -behaves like a price. - -**Converting units.** Units are read from the database. For `temoa_power`, millions of dollars -per PJ become dollars per MWh by multiplying by 3.6, since one PJ is 277,778 MWh. If the units -aren't recognized, the script stops and asks for the conversion explicitly. - -**Checking the sign.** No sign flip is applied by default (§1.4). As a safeguard, the script -warns if most nonzero balance duals in the run are negative, which would suggest a solver with -the opposite convention. A `--flip-sign` option is available for that case. - -**Combining myopic windows.** In myopic runs, each window's duals are saved separately, and -the dual table is never cleared between windows. It therefore also contains the look-ahead -results of earlier windows for periods that were later re-solved. For each period, the script -uses the duals from the latest window that starts at or before that period, which is the -solution Temoa kept in its other output tables. The per-slice output records which window each -price came from. - -**Averaging.** A per-slice price alone doesn't summarize a period, so the script reports two -averages for each region and period: - -- a **load-weighted average**, weighting each slice by the electricity drawn from the grid node - (storage charging and exports excluded); -- a **time-weighted average**, weighting each slice by its share of the year. - -**Outputs.** - -- A per-slice CSV: raw dual, discount factor, price, load, slice length, and source window. -- An annual CSV, by region and period: both averages, minimum, median and maximum prices, total - load, and counts of negative, zero and very high (default above 500 $/MWh) slice prices. -- A printed table of load-weighted prices by region and period, plus diagnostics flagging - negative prices, zero prices in slices with load, price spikes, and periods where the two - averages differ sharply. - -**Safety.** The script opens the database read-only and changes nothing in it. - -**Validation.** On copies of the Utopia tutorial database: - -- a fuel with a known 2.0 M$/PJ import cost was priced at exactly 7.20 $/MWh in every period, - in both perfect-foresight and myopic runs (confirming the discount factor, units, sign and - myopic base year); -- in a myopic run with overlapping windows, the script correctly ignored an earlier window's - look-ahead dual for a period that was later re-solved; -- perfect-foresight and myopic electricity prices differed only by amounts explained by their - different investment decisions. - -**Usage.** Solve the model with `save_duals = true` and a solver that returns duals. Gurobi -works; the `appsi_highs` interface does not return duals, and Temoa's Monte Carlo mode disables -them. Then run: - -``` -uv run python scripts/electricity_prices.py temoa_power.sqlite --commodity ELC ELCP --out-dir price_out -``` - ---- - -## 4. Caveats when reading the prices - -**These are marginal values.** A dual is the slope of the cost curve at the optimum, valid for -small changes. One PJ in a 57-hour slice is roughly 4.9 GW of extra load, enough to change -which plant or constraint is marginal. Treat "the cost of one more PJ" as the cost of the first -small amount, extended in a straight line. Where the solution is degenerate (several equally -good solutions), adding load and removing load can have different marginal costs, and the -solver reports one of them. - -**A slice is a representative hour, not a calendar hour.** Each price applies to one -representative hour that stands for about 57 real hours per year in `S1` and 34 in `S2`. Because -capacity costs are concentrated on the few hours of binding slices, peak-slice prices above -1,000 $/MWh are expected when new capacity is needed. - -**Individual peak prices can be arbitrary.** When several slices, or several periods, tie at -the binding peak, how the capacity cost is split among them is not unique (case 2 in §2.4). The -symptoms are one very expensive slice next to cheap neighbours, or prices that jump between -periods without any change in costs or load. In the Utopia test, the winter day and night -slices split into +124 and −49 $/MWh in one period, and +160 and −103 in another, while -all other slices were around 19–26 $/MWh. That pattern reflects a constraint linking the two -slices, such as storage or a fixed day/night ratio, not a real negative price. **Load-weighted -period averages are robust to this; individual peak-slice prices are only indicative.** - -**Negative prices** can be real. Because the balance is an equality, surplus production with -nowhere to go can have a negative value: for example must-run output, renewable mandates, or -fixed output ratios. They can also arise from the degeneracy above. The script counts them; -check whether they line up with a plausible cause. - -**The price depends on where it is measured.** `ELC` is measured after grid losses. The busbar -price at `ELCP` is roughly 0.953 × the `ELC` price. Prices delivered to the residential, -commercial and transport sectors add their delivery costs. The script can price several nodes -at once. - -**Policy costs are included.** If a CO₂ cap or a renewable mandate binds, its cost per MWh is -in the price. That is correct for marginal cost under the policy. Separating out the physical -resource cost requires breaking the price into parts using those constraints' duals (§2.1). - -**This is a planning model's marginal cost, not a market price.** It assumes investment -responds optimally within each solve. Capacity costs appear as very high prices concentrated in -binding slices, rather than as a separate capacity payment. In perfect foresight, prices can -include the capital of plants built in earlier periods (§2.4); in myopic runs they never do -(§2.5). - -**Myopic prices reflect limited foresight.** They are correct for each window as solved, given -the fleet it inherited. Differences from perfect-foresight prices come from what the planner -cannot see and from window-boundary effects, not from the missing costs in the reports (§2.6). -Nothing guarantees that earlier investments recover their capital. - -**Solver settings affect precision.** Temoa sets Gurobi to use the barrier algorithm, stop -without the final "crossover" step to an exact corner solution, and accept a loose convergence -tolerance (`temoa/_internal/run_actions.py`, lines 219–222). The duals are therefore only -approximately optimal: expect small nonzero prices where zero would be exact, and slice prices -that are off by a few percent, especially in a model as large as `temoa_power`. Without the -final step, the solver also tends to spread a tied capacity cost evenly across slices and -periods instead of placing it in one. That makes hourly profiles look smoother, but the split -remains a numerical choice, not something the model determines. **For results that matter, -re-solve once with a tight tolerance and crossover enabled, and compare.** If slice prices move -a lot, report period averages rather than hourly shapes. - -**How the reserve margin works matters.** Temoa's reserve requirement scales with the output of -reserve-eligible plants in each slice, not with load directly. The 1.5× capacity factor in §2.2 -applies when the extra supply comes from a reserve-eligible plant. Temoa's static and dynamic -reserve methods also credit capacity differently, so check which one a run used. - -**Reported costs and prices are different questions.** The understated capital in myopic cost -outputs (§2.6) does not affect the prices, but it does affect any comparison between prices and -costs. - ---- - -## 5. Comparing marginal prices with average cost - -A natural check on the marginal prices is to compare them with the average cost of -electricity: total cost divided by electricity produced. Doing that correctly requires knowing -exactly what Temoa's cost table reports, and the result is informative whichever way the -comparison comes out. - -### 5.1 What the `output_cost` columns measure - -Temoa reports costs in `output_cost`, with a discounted and an undiscounted version of each -cost type. The columns do not all cover the same span of time. - -**Fixed, variable and emission costs cover the whole period.** Temoa calculates each as one -representative year's cost and multiplies it by the period discount factor DF_p (§1.6). So -`d_fixed`, `d_var` and `d_emiss` are the discounted totals for all years of the period. The -undiscounted columns `fixed`, `var` and `emiss` are one year's cost multiplied by the period -length. This was confirmed on the Utopia test database, which has 10-year periods: the -recorded variable cost was exactly 10 times one year's activity times its variable cost, and -the ratio of discounted to undiscounted cost was 0.7722, which is DF_1990 / 10. Fixed costs gave -the same ratio. In other words, `d_x / DF_p` and `x / LEN_p` both give the cost of a single -representative year. - -**Investment cost is a lump sum recorded in the year capacity is built.** `d_invest` appears -under the period equal to the plant's vintage. It is the entire stream of capital payments the -model charges for that capacity, discounted to the base year. That stream covers the plant's -life up to the end of the horizon: 2055 in perfect foresight, or the end of the window in -myopic runs, which is the reporting gap described in §2.6. The undiscounted `invest` column is -the same stream without discounting. A plant built in 2025 therefore shows all of its -2025–2054 payments in its 2025 row, and nothing in later rows. - -### 5.2 Computing an average cost that matches the prices - -To compare like with like, the average cost has to be on the same basis as the marginal prices -produced by the script: undiscounted dollars for one year of the period, per MWh at the `ELC` -node. - -- **Fixed, variable and emission costs:** dividing discounted cost by discounted electricity - for a single period works, because the discount factor cancels. It is equivalent to - `(fixed + var + emiss) / LEN_p`, divided by one year's electricity. -- **Investment: do not use `d_invest` as listed.** Because it places decades of capital in the - build period and none afterwards, it would make average cost spike in build periods and fall - too low in all others. Instead, add up the yearly capital payment (for example about - 26.8 M$/GW-yr for the gas turbine in §2.2) for every vintage still in service in that - period. In myopic runs this also avoids the truncated capital reporting. -- **Use the same scope as the price.** The `ELC` price includes fuel valued at its marginal - cost, grid losses, storage and transmission. The average cost should include upstream fuel - costs (or value fuel at its own dual), and the costs of the grid, storage and transmission - technologies. -- **Use the same energy.** Divide by electricity consumed at `ELC`, the same quantity the - script uses for load-weighting, not by generation at `ELCP`. Generation is larger because of - grid losses, so dividing by it would understate the cost per MWh delivered. - -### 5.3 Should average cost always be lower than the average marginal price? - -No. The comparison can go either way, and which way it goes tells you something about the -system. - -**The benchmark case.** Linear programs satisfy an exact accounting identity: at the optimum, -total cost equals the sum over all constraints of each constraint's dual multiplied by its -right-hand side. Suppose every plant in use is newly built at a constant cost, and nothing else -binds: no capacity or resource limits, no existing capacity, no policy constraints. Then revenue -at the marginal prices exactly covers cost, and **the load-weighted average marginal price -equals the average cost.** Departures from that benchmark identify which of the following -situations applies. - -**When the average marginal price is above average cost.** Marginal prices then pay rents to -something whose cost is missing from, or understated in, `output_cost`: - -- **Existing capacity.** Plants built before 2020 have no investment cost in `output_cost`, and - in myopic runs the capital of earlier windows is only partly recorded. When that capacity is - scarce, prices still reflect the cost of new capacity. -- **Limited low-cost resources**, such as caps on wind, solar or hydro sites, capacity upper - bounds, or resource limits. The cheap units earn the difference between the price and their - cost. -- **Binding CO₂ caps or renewable mandates.** The price includes the constraint's dual times - the emissions per unit of output. Unless the database also sets an explicit emission price - (`cost_emission`), that cost appears nowhere in `output_cost`. -- **Upstream rents.** Fuel enters the price at its marginal value, which can be well above its - average supply cost. - -**When average cost is above the average marginal price.** Costs are then being incurred that -the marginal price does not pay for: - -- **Excess or forced capacity.** Capacity overbuilt relative to later needs, which is common in - myopic runs, as well as minimum-build constraints, minimum-activity constraints, or required - shares. -- **Falling technology costs.** Older plants cost more than today's new build, and today's new - build is what sets the marginal price. -- **Surplus hours.** Zero or negative prices in slices with surplus electricity pull the average - price down. -- **Solver imprecision.** Temoa's loose barrier settings make the duals approximate (§4). - -### 5.4 What to expect for `temoa_power` - -In early periods, the average marginal price will probably be above average cost, because of -existing capacity, resource limits and any binding policy. In later periods, if new builds -dominate and costs are stable, the two should move closer together. In myopic runs, recorded -investment is truncated (§2.6), which biases a measured average cost downward, so the capital -accounting should be rebuilt before drawing conclusions. - -When average cost exceeds the marginal price, look for stranded or forced capacity. When the -marginal price is well above average cost, the rents can be traced to the binding constraints -with large duals. - ---- - -## 6. Quick reference - -| Item | `temoa_power` | -|---|---| -| Electricity balance constraint | `commodity_balance_constraint[r,p,s,d,ELC]`: equality, one per slice | -| Meaning of a positive dual | marginal cost of serving more load (no sign flip needed) | -| Raw dual units | M$ discounted to 2020, per PJ, in `output_dual_variable` | -| Period discount factor | sum of yearly discount factors over the 5-year period, end-of-year convention, 2% rate, base year 2020 | -| Conversion to a price | `price ($/MWh) = dual ÷ DF_p × 3.6` | -| Time slices | 2 seasons × 96 representative hours; each stands for ≈ 57 h/yr (`S1`) or 34 h/yr (`S2`) | -| Myopic duals | stored per window; the script uses the latest window starting at or before each period | -| Earlier-period capital in prices | possible in perfect foresight; never in myopic runs | -| Myopic cost reports | investment costs cut off at each window's end; prices unaffected | -| `output_cost` time span | fixed, variable and emission costs: whole period; investment: lump sum of in-horizon payments in the vintage row | -| Average cost vs. marginal price | equal only in the benchmark case; the direction of the gap points to rents or forced costs (§5) | - -## Appendix: status of the `temoa_power` analysis - -At the time of writing, `temoa_power.sqlite` contains model inputs only; it has not yet been -solved with duals saved. The `temoa_power` figures above (time slices, discount factors, gas -turbine costs) are therefore computed from its inputs, and the tested numerical results come -from the Utopia tutorial database. The next step is to solve `temoa_power` in myopic mode with -`save_duals = true` using Gurobi, then run the script as shown in §3. diff --git a/scripts/reserve_dual_toy.py b/scripts/reserve_dual_toy.py deleted file mode 100644 index f677b79d9..000000000 --- a/scripts/reserve_dual_toy.py +++ /dev/null @@ -1,68 +0,0 @@ -"""Toy capacity-expansion LP with Temoa's discounting, to show where a 2025 turbine's capital -shows up in later reserve-margin duals. One region, one peak slice, gas turbines only.""" - -import highspy -import numpy as np - -g, p0, periods, p_end = 0.02, 2020, [2025, 2030, 2035, 2040, 2045, 2050], 2055 -A, F = 26.77, 20.527 # yearly capital-equivalent and fixed O&M, M$/GW-yr -cc, prm = 0.9, 0.35 # capacity credit, reserve margin -pa = lambda n: ((1 + g) ** n - 1) / (g * (1 + g) ** n) -DF = {p: pa(5) / (1 + g) ** (p - p0) for p in periods} # same factor as my script -# Temoa: loan payments from vintage v to horizon end, FOM every active period -> both sum DF_p, p>=v -cost = {v: (A + F) * sum(DF[p] for p in periods if p >= v) for v in periods} - - -def solve(load, existing, allow_build): - h = highspy.Highs() - h.setOptionValue('output_flag', False) - n = len(periods) - for v in periods: # N_v >= 0 - h.addVar(0, highspy.kHighsInf if allow_build(v) else 0) - h.changeColCost(periods.index(v), cost[v]) - for i, p in enumerate(periods): # cc*(E + sum_{v<=p} N_v) >= (1+prm)*load_p - idx = [j for j, v in enumerate(periods) if v <= p] - h.addRow( - (1 + prm) * load[p] - cc * existing, - highspy.kHighsInf, - len(idx), - np.array(idx, dtype=np.int32), - np.full(len(idx), cc), - ) - h.run() - sol = h.getSolution() - N = sol.col_value - rho = sol.row_dual - print( - ' period build_GW dual(disc.) dual/DF=M$/GW-yr of requirement $/MWh adder in 57-h slice' - ) - for i, p in enumerate(periods): - d = (1 + prm) * rho[i] # d obj / d (1 GW more peak load) - print( - f' {p} {N[i]:7.2f} {d:9.2f} {d / DF[p]:7.2f} {d / DF[p] * 1e6 / 57000:7.0f}' - ) - print( - ' check: cost of 2025 turbine =', - round(cost[2025] * cc / cc, 2), - ' = sum of cc*rho_p over periods it is active:', - round(sum(cc * rho[i] for i in range(len(periods))), 2) if N[0] > 1e-9 else '(not built)', - ) - - -print('DF: ' + ', '.join(f'{p}:{DF[p]:.3f}' for p in periods)) -print( - f'(A+F)={A + F:.2f} M$/GW-yr; with reserve factor (1+prm)/cc={(1 + prm) / cc:.2f} -> {(A + F) * (1 + prm) / cc:.1f} M$/GW-yr of load\n' -) -base = 10.0 -print('CASE 1: peak load grows 1 GW every period -> a new turbine every period') -solve({p: base + i + 1 for i, p in enumerate(periods)}, base * (1 + prm) / cc, lambda v: True) -print('\nCASE 2: peak steps up 1 GW in 2025 and then stays flat -> only the 2025 turbine') -solve(dict.fromkeys(periods, base + 1), base * (1 + prm) / cc, lambda v: True) -print( - '\nCASE 3: 2025 has spare capacity, peak steps up 1 GW in 2030, but 2030+ builds are not allowed' -) -solve( - {p: base if p == 2025 else base + 1 for p in periods}, - base * (1 + prm) / cc, - lambda v: v == 2025, -)