diff --git a/CHANGELOG.md b/CHANGELOG.md index bdfe1ce8..77cea349 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,7 +1,10 @@ # Change Log -## [4.18.17RC] - 2026-06-03 Unreleased in PyPI +## [4.18.17RC] - 2026-06-04 Unreleased in PyPI +- [CHANGED] update the architecture diagram (`doc/img/openTEPES_architecture.svg` and the rendered `.png`) to show the new `openTEPES_ProblemSolvingStageIter.py` as an implemented (green) box in the solver layer, alongside `ProblemSolving`, `Tuning`, `DualExtraction`, `Persistent` and `Benders`; `resolve.py` stays white as the one solver module still planned. The seven boxes are re-spaced to fit the layer's existing width — no other layer changes. +- [ADDED] a test (`tests/test_direct_run.py`) that runs every `openTEPES_*.py` module as a script (via `runpy`, with `__name__ == "__main__"`) to check that each one still imports cleanly that way — i.e. the `try`/`except ImportError` relative-import guard's fallback works. A normal `import openTEPES` only ever takes the relative-import path, so until now nothing checked the fallback; a broken guard, a missing import, or a circular import in any module would only surface when a user runs the file directly (VS Code "Run Python File"). The test solves no model, so it runs in the fast CI job; it excludes `openTEPES_Main.py`, which has a real `__main__` block that would launch the command-line interface. +- [CHANGED] move the per-stage solve loop out of `openTEPES_run` into a new module `openTEPES_ProblemSolvingStageIter.py`. The `(period, scenario, stage)` loop — which activates one stage's load levels, builds that stage's operation constraints, and calls `ProblemSolving` — now lives in a function `StageIterativeSolving`, together with the post-loop work that rebuilds the full stage and load-level sets, reactivates every constraint, and restores the scenario probabilities. `openTEPES.py` still builds the objective and the investment constraints, sets `First_st` / `Last_st` and the empty `pDuals`, and then calls the new driver. The two solve paths are unchanged: a deterministic solve per scenario when there are no expansion decisions and no system emission or RES-energy limit, and one joint stochastic solve at the last period-scenario otherwise. This is a pure move — the loop body is copied unchanged — so results do not change: the whole solve test suite passes with the same locked costs, and a full-output run of case 9n gives the same total cost (252.20132998 MEUR, matching to 13 significant figures; only the per-unit dispatch tables move, by the amount HiGHS already varies run-to-run when it picks among equally optimal storage schedules). The new module uses the same `try`/`except ImportError` relative-import guard as the other split modules, so "Run Python File" still works on it. This continues the solver-layer split (after the Persistent, Tuning and DualExtraction modules) and gives the later re-solve and decomposition drivers a single place to call the stage loop from. - [CHANGED] split the ~2800-line `openTEPES_OutputResults.py` into per-concern modules at the package root — the output-side counterpart of the earlier `openTEPES_Input*` and `openTEPES_ProblemSolving*` splits. The result writers now live in `openTEPES_OutputResultsInvestment.py`, `...Generation.py`, `...Storage.py`, `...Hydrogen.py` (hydrogen network), `...Heat.py` (heat network), `...Network.py` (electricity network and map), `...Economic.py` (marginal, cost-summary, economic), `...Summary.py` (system summary, flexibility, reliability), and `...RawDump.py` (the raw parameter/variable/constraint DuckDB dump). Shared pieces sit in `openTEPES_OutputResultsCommon.py` (the output-directory helper and the three Altair plot builders) and `openTEPES_OutputResultsMapCommon.py` (the flow-series and snapshot-period helpers the three network maps had each copied). The old `openTEPES_OutputResults.py` is removed: `__init__.py` and `openTEPES.py` import the result functions from the concern modules directly, the same way the Input and ProblemSolving splits are wired (no re-export facade). Functions are grouped by topic; the dispatch order is still controlled by `OUTPUT_REGISTRY` in `openTEPES.py` and is unchanged. Every cross-module import uses the same `try`/`except ImportError` relative-import guard as the other split modules, so "Run Python File" works. Verified bit-identical with the output parity tool on case 9n (`PYTHONHASHSEED=0`, all 87 result tables; total cost 164.382043867 MEUR). Two small cleanups made while moving the code: removed 34 no-op `"".join([f"..."])` wrappers (each is just the f-string) and stopped `ESSOperationResults` from writing its inventory-utilization scaling back into `mTEPES.pMaxStorage` (the denominator is now computed locally, so the scaled value no longer leaks to later reads of the parameter). - [ADDED] a `9n_H2` example case (the hydrogen counterpart of `9n_heat`): the minimal 9-node electricity system plus two electrolyzers, hydrogen demand at three nodes, and a hydrogen pipe network. Added to the single-stage CI solve suite in `tests/test_run.py` (expected cost 259.547 MEUR under the 7-day HiGHS fixture, verified deterministic across reruns), so CI exercises the hydrogen result writer (`NetworkH2OperationResults`) directly on a small fast case rather than only through `sSEP`. Output verified bit-identical before and after the sector-coupling module split, for both CSV and DuckDB inputs. - [FIXED] two operating-reserve guards in the marginal and economic results read a leaked loop variable. The down-reserve-marginal guard and the up-reserve-revenue guard loop over storage units (`for ar,es in mTEPES.ar*mTEPES.es`) but checked `pIndOperReserveGen[nr]` / `pIndOperReserveCon[nr]` — `nr` only existed as a leftover from the earlier `for ar,nr ...` loop that builds the area-to-generator map. They now check the loop's own `es`, matching the two sibling guards (up-reserve-marginal, down-reserve-revenue) that were already correct. The down-reserve-revenue guard also read `pIndOperReserveGen[nr]` twice (`Gen ... or ... Gen`) where its siblings read `Gen ... or ... Con`; the second is now `Con`. On case 9n the result is unchanged; on a case where the leftover unit's reserve flags differ from the storage units' these guards now correctly decide which operating-reserve-revenue tables are written. diff --git a/doc/img/openTEPES_architecture.png b/doc/img/openTEPES_architecture.png index bf78baec..d2b0dc1f 100644 Binary files a/doc/img/openTEPES_architecture.png and b/doc/img/openTEPES_architecture.png differ diff --git a/doc/img/openTEPES_architecture.svg b/doc/img/openTEPES_architecture.svg index a1a81978..6f48c43e 100644 --- a/doc/img/openTEPES_architecture.svg +++ b/doc/img/openTEPES_architecture.svg @@ -250,33 +250,37 @@ openTEPES_Problem- Solving*.py - - ProblemSolving.py - orchestrator entry + + ProblemSolving.py + orchestrator entry - - …Tuning.py - solver options + + …StageIter.py + stage loop - - …DualExtraction.py - duals via re-solve + + …Tuning.py + solver options - - …Persistent.py - appsi persistent + + …DualExtraction.py + duals via re-solve - - …Benders.py - L-shaped decomp + + …Persistent.py + appsi persistent - - resolve.py - ★ hot-swap (Mode C) + + …Benders.py + L-shaped decomp + + + resolve.py + ★ hot-swap (Mode C) - + Mode C: Param.store_values() + re-solve diff --git a/openTEPES/__init__.py b/openTEPES/__init__.py index 92814d50..17fda951 100644 --- a/openTEPES/__init__.py +++ b/openTEPES/__init__.py @@ -24,8 +24,9 @@ from .openTEPES_InputDuckDBSource import * from .openTEPES_InputData import * from .openTEPES_ModelFormulation import * -from .openTEPES_ProblemSolving import * -from .openTEPES_ProblemSolvingBenders import * +from .openTEPES_ProblemSolving import * +from .openTEPES_ProblemSolvingStageIter import * +from .openTEPES_ProblemSolvingBenders import * from .openTEPES_OutputResultsCommon import * from .openTEPES_OutputResultsRawDump import * from .openTEPES_OutputResultsInvestment import * diff --git a/openTEPES/openTEPES.py b/openTEPES/openTEPES.py index 7cc00864..1bd6c81d 100644 --- a/openTEPES/openTEPES.py +++ b/openTEPES/openTEPES.py @@ -1,5 +1,5 @@ """ -Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - June 3, 2026 +Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - June 4, 2026 """ # import dill as pickle @@ -17,8 +17,8 @@ try: from .openTEPES_InputData import InputData, DataConfiguration, SettingUpVariables from .openTEPES_InputSource import open_source - from .openTEPES_ModelFormulation import TotalObjectiveFunction, InvestmentElecModelFormulation, InvestmentHydroModelFormulation, InvestmentH2ModelFormulation, InvestmentHeatModelFormulation, GenerationOperationModelFormulationObjFunct, GenerationOperationElecModelFormulationInvestment, GenerationOperationHeatModelFormulationInvestment, GenerationOperationModelFormulationDemand, GenerationOperationModelFormulationStorage, GenerationOperationModelFormulationReservoir, NetworkH2OperationModelFormulation, NetworkHeatOperationModelFormulation, GenerationOperationModelFormulationCommitment, GenerationOperationModelFormulationRampMinTime, NetworkSwitchingModelFormulation, NetworkOperationModelFormulation, NetworkCycles, CycleConstraints - from .openTEPES_ProblemSolving import ProblemSolving + from .openTEPES_ModelFormulation import TotalObjectiveFunction, InvestmentElecModelFormulation, InvestmentHydroModelFormulation, InvestmentH2ModelFormulation, InvestmentHeatModelFormulation + from .openTEPES_ProblemSolvingStageIter import StageIterativeSolving from .openTEPES_OutputResultsRawDump import OutputResultsParVarCon from .openTEPES_OutputResultsInvestment import InvestmentResults from .openTEPES_OutputResultsGeneration import GenerationOperationResults, GenerationOperationHeatResults @@ -33,8 +33,8 @@ sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) from openTEPES.openTEPES_InputData import InputData, DataConfiguration, SettingUpVariables from openTEPES.openTEPES_InputSource import open_source - from openTEPES.openTEPES_ModelFormulation import TotalObjectiveFunction, InvestmentElecModelFormulation, InvestmentHydroModelFormulation, InvestmentH2ModelFormulation, InvestmentHeatModelFormulation, GenerationOperationModelFormulationObjFunct, GenerationOperationElecModelFormulationInvestment, GenerationOperationHeatModelFormulationInvestment, GenerationOperationModelFormulationDemand, GenerationOperationModelFormulationStorage, GenerationOperationModelFormulationReservoir, NetworkH2OperationModelFormulation, NetworkHeatOperationModelFormulation, GenerationOperationModelFormulationCommitment, GenerationOperationModelFormulationRampMinTime, NetworkSwitchingModelFormulation, NetworkOperationModelFormulation, NetworkCycles, CycleConstraints - from openTEPES.openTEPES_ProblemSolving import ProblemSolving + from openTEPES.openTEPES_ModelFormulation import TotalObjectiveFunction, InvestmentElecModelFormulation, InvestmentHydroModelFormulation, InvestmentH2ModelFormulation, InvestmentHeatModelFormulation + from openTEPES.openTEPES_ProblemSolvingStageIter import StageIterativeSolving from openTEPES.openTEPES_OutputResultsRawDump import OutputResultsParVarCon from openTEPES.openTEPES_OutputResultsInvestment import InvestmentResults from openTEPES.openTEPES_OutputResultsGeneration import GenerationOperationResults, GenerationOperationHeatResults @@ -182,11 +182,11 @@ def openTEPES_run(DirName, CaseName, SolverName, pIndOutputResults, pIndLogConso idxDict['y' ] = 1 #%% model declaration - mTEPES = ConcreteModel('Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - Version 4.18.17RC - June 3, 2026') + mTEPES = ConcreteModel('Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - Version 4.18.17RC - June 4, 2026') # In DuckDB-input mode _path may not exist on disk (the case lives in # the DB, not in a directory). Ensure the version-log target exists. os.makedirs(_path, exist_ok=True) - print( 'Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - Version 4.18.17RC - June 3, 2026', file=open(f'{_path}/openTEPES_version_{CaseName}.log','w')) + print( 'Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - Version 4.18.17RC - June 4, 2026', file=open(f'{_path}/openTEPES_version_{CaseName}.log','w')) if _input_source is not None: mTEPES.pInputSource = _input_source @@ -233,274 +233,10 @@ def openTEPES_run(DirName, CaseName, SolverName, pIndOutputResults, pIndLogConso # initialize parameter for dual variables mTEPES.pDuals = {} - # initialize the set of load levels up to the current stage - mTEPES.na = Set(initialize=[]) - - # iterative model formulation for each stage of a year - for p,sc,st in mTEPES.ps*mTEPES.stt: - - # control for not repeating the stages in case of several ones - mTEPES.NoRepetition = 0 - - # activate only load levels to formulate - mTEPES.del_component(mTEPES.st) - mTEPES.del_component(mTEPES.n ) - mTEPES.del_component(mTEPES.n2) - mTEPES.st = Set(doc='stages', initialize=[stt for stt in mTEPES.stt if st == stt if mTEPES.pStageWeight[stt] and sum(1 for (p,sc,st,nn) in mTEPES.s2n)]) - mTEPES.n = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if (p,sc,st,nn) in mTEPES.s2n ]) - mTEPES.n2 = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if (p,sc,st,nn) in mTEPES.s2n ]) - - # load levels multiple of cycles for each ESS/generator - mTEPES.nesc = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [es] == 0] - mTEPES.necc = [(n,ec) for n,ec in mTEPES.n*mTEPES.ec if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [ec] == 0] - mTEPES.neso = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pOutflowsTimeStep[es] == 0] - mTEPES.ngen = [(n,g ) for n,g in mTEPES.n*mTEPES.g if mTEPES.n.ord(n) % mTEPES.pEnergyTimeStep [g ] == 0] - if mTEPES.pIndHydroTopology: - mTEPES.nhc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2h) == 0] - if sum(1 for h,rs in mTEPES.p2r): - mTEPES.np2c = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) == 0] - else: - mTEPES.np2c = [] - if sum(1 for rs,h in mTEPES.r2p): - mTEPES.npc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) == 0] - else: - mTEPES.npc = [] - mTEPES.nrsc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] - mTEPES.nrcc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rn if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] - mTEPES.nrso = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pWaterOutTimeStep [rs] == 0] - - if mTEPES.st: - - # load levels up to the current stage for emission and RES energy constraints - # if (p != mTEPES.p.first() or sc != mTEPES.sc.first()) and st == mTEPES.First_st: - # mTEPES.del_component(mTEPES.na) - if st == mTEPES.First_st: - mTEPES.del_component(mTEPES.na) - mTEPES.na = Set(initialize=mTEPES.n) - else: - temp_na = mTEPES.na | mTEPES.n # Create a temporary set - mTEPES.del_component(mTEPES.na) # Delete the existing set to avoid overwriting - mTEPES.na = Set(initialize=temp_na) # Assign the new set - - print(f'Period {p}, Scenario {sc}, Stage {st}') - - # operation model objective function and constraints by stage - GenerationOperationModelFormulationObjFunct (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - GenerationOperationElecModelFormulationInvestment (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - if mTEPES.pIndHeat: - GenerationOperationHeatModelFormulationInvestment(mTEPES, mTEPES, pIndLogConsole, p, sc, st) - GenerationOperationModelFormulationDemand (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - GenerationOperationModelFormulationStorage (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - if mTEPES.pIndHydroTopology: - GenerationOperationModelFormulationReservoir (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - if mTEPES.pIndHydrogen: - NetworkH2OperationModelFormulation (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - if mTEPES.pIndHeat: - NetworkHeatOperationModelFormulation (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - GenerationOperationModelFormulationCommitment (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - GenerationOperationModelFormulationRampMinTime (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - NetworkSwitchingModelFormulation (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - NetworkOperationModelFormulation (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - # introduce cycle flow formulations - if pIndCycleFlow: - if st == mTEPES.First_st: - NetworkCycles ( mTEPES, pIndLogConsole ) - CycleConstraints (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - - # No generator investments and no generator retirements, and no line investments - if ( (sum(1 for pp,eb in mTEPES.peb if pp <= p) == 0 or mTEPES.pIndBinGenInvest() == 2) - and (sum(1 for pp,gd in mTEPES.pgd if pp <= p) == 0 or mTEPES.pIndBinGenRetire() == 2) - and (sum(1 for pp,rc in mTEPES.prc if pp <= p) == 0 or mTEPES.pIndBinRsrInvest() == 2) - and (sum(1 for pp,ni,nf,cc in mTEPES.plc if pp <= p) == 0 or mTEPES.pIndBinNetElecInvest() == 2) - and (sum(1 for pp,ni,nf,cc in mTEPES.ppc if pp <= p) == 0 or mTEPES.pIndBinNetH2Invest() == 2) - and (sum(1 for pp,ni,nf,cc in mTEPES.phc if pp <= p) == 0 or mTEPES.pIndBinNetHeatInvest() == 2)): - - # No minimum RES requirements and no emission limit - if (max([mTEPES.pRESEnergy[p,ar] for ar in mTEPES.ar]) == 0 and (min([mTEPES.pEmission [p,ar] for ar in mTEPES.ar]) == math.inf or sum(mTEPES.pEmissionRate[nr] for nr in mTEPES.nr) == 0)): - - if (p,sc) == mTEPES.ps.last() and st == mTEPES.Last_st and mTEPES.NoRepetition == 0: - mTEPES.NoRepetition = 1 - - # Writing LP file - if pIndLogConsole: - StartTime = time.time() - mTEPES.write(f'{_path}/openTEPES_{CaseName}_{p}_{sc}_{st}.lp', io_options={'symbolic_solver_labels': True}) - WritingLPFileTime = time.time() - StartTime - StartTime = time.time() - print('Writing LP file ... ', round(WritingLPFileTime), 's') - - # there are no expansion decisions, or they are ignored (it is an operation model), or there are no system emission nor minimum RES constraints - mTEPES.pScenProb[p,sc] = 1.0 - ProblemSolving(DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, mTEPES.p.ord(p)*mTEPES.sc.ord(sc)*mTEPES.st.ord(st)) - mTEPES.pScenProb[p,sc] = 0.0 - - # deactivate the constraints of the previous period and scenario - for c in mTEPES.component_objects(pyo.Constraint, active=True): - if c.name.find(str(p)) != -1 and c.name.find(str(sc)) != -1: - c.deactivate() - - # Minimum RES requirements or emission limit - elif (max([mTEPES.pRESEnergy[p,ar] for ar in mTEPES.ar]) > 0 or (min([mTEPES.pEmission [p,ar] for ar in mTEPES.ar]) < math.inf and sum(mTEPES.pEmissionRate[nr] for nr in mTEPES.nr) > 0)): - - if st == mTEPES.Last_st and mTEPES.NoRepetition == 0: - - if (p,sc) == mTEPES.ps.last(): - mTEPES.NoRepetition = 1 - - mTEPES.del_component(mTEPES.st) - mTEPES.del_component(mTEPES.n ) - mTEPES.del_component(mTEPES.n2) - mTEPES.st = Set(doc='stages', initialize=[stt for stt in mTEPES.stt if mTEPES.pStageWeight[stt] and sum(1 for p,sc,stt,nn in mTEPES.s2n)]) - mTEPES.n = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) - mTEPES.n2 = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) - - # load levels multiple of cycles for each ESS/generator - mTEPES.nesc = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [es] == 0] - mTEPES.necc = [(n,ec) for n,ec in mTEPES.n*mTEPES.ec if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [ec] == 0] - mTEPES.neso = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pOutflowsTimeStep[es] == 0] - mTEPES.ngen = [(n,g ) for n,g in mTEPES.n*mTEPES.g if mTEPES.n.ord(n) % mTEPES.pEnergyTimeStep [g ] == 0] - if mTEPES.pIndHydroTopology: - mTEPES.nhc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2h) == 0] - if sum(1 for h,rs in mTEPES.p2r): - mTEPES.np2c = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) == 0] - else: - mTEPES.np2c = [] - if sum(1 for rs,h in mTEPES.r2p): - mTEPES.npc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) == 0] - else: - mTEPES.npc = [] - mTEPES.nrsc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] - mTEPES.nrcc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rn if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] - mTEPES.nrso = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pWaterOutTimeStep [rs] == 0] - - # Writing an LP file - if pIndLogConsole: - StartTime = time.time() - mTEPES.write(f'{_path}/openTEPES_{CaseName}_{p}_{sc}_{st}.lp', io_options={'symbolic_solver_labels': True}) - WritingLPFileTime = time.time() - StartTime - StartTime = time.time() - print('Writing LP file ... ', round(WritingLPFileTime), 's') - - mTEPES.pScenProb[p,sc] = 1.0 - # there are system emission or RES requirement constraints - ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, mTEPES.p.ord(p)*mTEPES.sc.ord(sc)*mTEPES.st.ord(st)) - # if pIndSectorDecomposition and len(mTEPES.psnel): - # ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, 1) - # SectorDecomposition(DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st ) - # else: - # ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, mTEPES.p.ord(p)*mTEPES.sc.ord(sc)*mTEPES.st.ord(st)) - mTEPES.pScenProb[p,sc] = 0.0 - - # deactivate the constraints of the previous period and scenario - for c in mTEPES.component_objects(pyo.Constraint, active=True): - if c.name.find(str(p)) != -1 and c.name.find(str(sc)) != -1: - c.deactivate() - - # generator investments or generator retirements, or line investments - elif ( (sum(1 for pp,eb in mTEPES.peb if pp <= p) > 0 and mTEPES.pIndBinGenInvest() < 2) - or (sum(1 for pp,gd in mTEPES.pgd if pp <= p) > 0 and mTEPES.pIndBinGenRetire() < 2) - or (sum(1 for pp,rc in mTEPES.prc if pp <= p) > 0 and mTEPES.pIndBinRsrInvest() < 2) - or (sum(1 for pp,ni,nf,cc in mTEPES.plc if pp <= p) > 0 and mTEPES.pIndBinNetElecInvest() < 2) - or (sum(1 for pp,ni,nf,cc in mTEPES.ppc if pp <= p) > 0 and mTEPES.pIndBinNetH2Invest() < 2) - or (sum(1 for pp,ni,nf,cc in mTEPES.phc if pp <= p) > 0 and mTEPES.pIndBinNetHeatInvest() < 2)): - - if (p,sc) == mTEPES.ps.last() and st == mTEPES.Last_st and mTEPES.NoRepetition == 0: - mTEPES.NoRepetition = 1 - - mTEPES.del_component(mTEPES.st) - mTEPES.del_component(mTEPES.n ) - mTEPES.del_component(mTEPES.n2) - mTEPES.st = Set(doc='stages', initialize=[stt for stt in mTEPES.stt if mTEPES.pStageWeight[stt] and sum(1 for p,sc,stt,nn in mTEPES.s2n)]) - mTEPES.n = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) - mTEPES.n2 = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) - - # load levels multiple of cycles for each ESS/generator - mTEPES.nesc = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [es] == 0] - mTEPES.necc = [(n,ec) for n,ec in mTEPES.n*mTEPES.ec if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [ec] == 0] - mTEPES.neso = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pOutflowsTimeStep[es] == 0] - mTEPES.ngen = [(n,g ) for n,g in mTEPES.n*mTEPES.g if mTEPES.n.ord(n) % mTEPES.pEnergyTimeStep [g ] == 0] - if mTEPES.pIndHydroTopology: - mTEPES.nhc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2h) == 0] - if sum(1 for h,rs in mTEPES.p2r): - mTEPES.np2c = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) == 0] - else: - mTEPES.np2c = [] - if sum(1 for rs,h in mTEPES.r2p): - mTEPES.npc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) == 0] - else: - mTEPES.npc = [] - mTEPES.nrsc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] - mTEPES.nrcc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rn if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] - mTEPES.nrso = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pWaterOutTimeStep [rs] == 0] - - # Writing an LP file - if pIndLogConsole: - StartTime = time.time() - mTEPES.write(f'{_path}/openTEPES_{CaseName}_{p}_{sc}_{st}.lp', io_options={'symbolic_solver_labels': True}) - WritingLPFileTime = time.time() - StartTime - StartTime = time.time() - print('Writing LP file ... ', round(WritingLPFileTime), 's') - - # there are investment decisions (it is an expansion and operation model), or there are system emission constraints - ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, 1) - # if pIndSectorDecomposition and len(mTEPES.psnel): - # ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, 1) - # SectorDecomposition(DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st ) - # else: - # ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, 1) - - mTEPES.del_component(mTEPES.st) - mTEPES.del_component(mTEPES.n ) - mTEPES.del_component(mTEPES.n2) - mTEPES.st = Set(doc='stages', initialize=[stt for stt in mTEPES.stt if mTEPES.pStageWeight[stt] and sum(1 for p,sc,stt,nn in mTEPES.s2n)]) - mTEPES.n = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) - mTEPES.n2 = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) - - # load levels multiple of cycles for each ESS/generator - mTEPES.nesc = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [es] == 0] - mTEPES.necc = [(n,ec) for n,ec in mTEPES.n*mTEPES.ec if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [ec] == 0] - mTEPES.neso = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pOutflowsTimeStep[es] == 0] - mTEPES.ngen = [(n,g ) for n,g in mTEPES.n*mTEPES.g if mTEPES.n.ord(n) % mTEPES.pEnergyTimeStep [g ] == 0] - if mTEPES.pIndHydroTopology: - mTEPES.nhc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2h) == 0] - if sum(1 for h,rs in mTEPES.p2r): - mTEPES.np2c = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) == 0] - else: - mTEPES.np2c = [] - if sum(1 for rs,h in mTEPES.r2p): - mTEPES.npc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) == 0] - else: - mTEPES.npc = [] - mTEPES.nrsc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] - mTEPES.nrcc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rn if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] - mTEPES.nrso = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pWaterOutTimeStep [rs] == 0] - - # activate the constraints of all the periods and scenarios - for c in mTEPES.component_objects(pyo.Constraint): - c.activate() - - # Output convention follows the two intended modes: - # - Case 2 (no expansion decisions): the per-scenario solve loop above - # used pScenProb=1.0 for each scenario individually. Reset all - # probabilities to 1.0 so cost-summary writers aggregate the - # independent deterministic solves as a sum. - # - Case 1 (with expansion decisions): the joint solve at the last - # (p,sc) used the input probabilities to minimise probability- - # weighted expected cost. Preserve those probabilities so the - # output reports that same expected cost — otherwise the writers - # would over-count by ~N_scenarios. - _hasExpansion = ( - (sum(1 for _ in mTEPES.peb) > 0 and mTEPES.pIndBinGenInvest() < 2) - or (sum(1 for _ in mTEPES.pgd) > 0 and mTEPES.pIndBinGenRetire() < 2) - or (sum(1 for _ in mTEPES.prc) > 0 and mTEPES.pIndBinRsrInvest() < 2) - or (sum(1 for _ in mTEPES.plc) > 0 and mTEPES.pIndBinNetElecInvest() < 2) - or (sum(1 for _ in mTEPES.ppc) > 0 and mTEPES.pIndBinNetH2Invest() < 2) - or (sum(1 for _ in mTEPES.phc) > 0 and mTEPES.pIndBinNetHeatInvest() < 2) - ) - if not _hasExpansion: - # assign probability 1 to all the periods and scenarios - for p,sc in mTEPES.ps: - mTEPES.pScenProb[p,sc] = 1.0 + # iterative formulation and solve for every stage of the year. The per-stage operation model and the two + # solve paths (deterministic per scenario, or one joint stochastic solve) live in + # openTEPES_ProblemSolvingStageIter; this is a pure extraction, so results are unchanged. + StageIterativeSolving(mTEPES, DirName, CaseName, SolverName, pIndLogConsole, _path, pIndCycleFlow) # pickle the case study data # with open(dump_folder+f'/oT_Case_{CaseName}.pkl','wb') as f: diff --git a/openTEPES/openTEPES_ProblemSolvingStageIter.py b/openTEPES/openTEPES_ProblemSolvingStageIter.py new file mode 100644 index 00000000..17728d5f --- /dev/null +++ b/openTEPES/openTEPES_ProblemSolvingStageIter.py @@ -0,0 +1,315 @@ +""" +openTEPES.openTEPES_ProblemSolvingStageIter — stage-by-stage formulate-and-solve driver. + +``StageIterativeSolving`` runs the model's per-stage loop: for every ``(period, scenario, stage)`` it activates +only that stage's load levels, builds the operation constraints for the stage, and calls ``ProblemSolving`` once +the stage's part of the model is in place. Two solve paths are preserved exactly as they were inside +``openTEPES_run``: + + * **no expansion / no system limits** — each scenario is solved on its own as a deterministic operation model + (``pScenProb`` is set to 1.0 for that single scenario, the solve runs, then its constraints are deactivated); + * **with expansion decisions or a system emission / RES-energy limit** — one joint stochastic solve is run at + the last ``(period, scenario)`` over all stages, using the input scenario probabilities. + +After the loop it rebuilds the full stage / load-level sets, reactivates every constraint, and restores the +scenario probabilities so the output writers report the intended (per-scenario sum, or expected-value) cost. + +This is a pure move out of ``openTEPES_run`` — the loop body is unchanged, so results are identical. It sets up +the seam that later orchestration drivers (Mode C re-solve, sector / Benders decomposition) build on. +""" +from __future__ import annotations + +import math +import os +import time + +import pyomo.environ as pyo +from pyomo.environ import Set + +# Support running this file directly (e.g. VS Code "Run Python File"), where __package__ is empty and the +# relative imports below have no parent package; fall back to absolute package imports in that case. +try: + from .openTEPES_ModelFormulation import GenerationOperationModelFormulationObjFunct, GenerationOperationElecModelFormulationInvestment, GenerationOperationHeatModelFormulationInvestment, GenerationOperationModelFormulationDemand, GenerationOperationModelFormulationStorage, GenerationOperationModelFormulationReservoir, NetworkH2OperationModelFormulation, NetworkHeatOperationModelFormulation, GenerationOperationModelFormulationCommitment, GenerationOperationModelFormulationRampMinTime, NetworkSwitchingModelFormulation, NetworkOperationModelFormulation, NetworkCycles, CycleConstraints + from .openTEPES_ProblemSolving import ProblemSolving +except ImportError: + import sys + sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__)))) + from openTEPES.openTEPES_ModelFormulation import GenerationOperationModelFormulationObjFunct, GenerationOperationElecModelFormulationInvestment, GenerationOperationHeatModelFormulationInvestment, GenerationOperationModelFormulationDemand, GenerationOperationModelFormulationStorage, GenerationOperationModelFormulationReservoir, NetworkH2OperationModelFormulation, NetworkHeatOperationModelFormulation, GenerationOperationModelFormulationCommitment, GenerationOperationModelFormulationRampMinTime, NetworkSwitchingModelFormulation, NetworkOperationModelFormulation, NetworkCycles, CycleConstraints + from openTEPES.openTEPES_ProblemSolving import ProblemSolving + + +def StageIterativeSolving(mTEPES, DirName, CaseName, SolverName, pIndLogConsole, _path, pIndCycleFlow): + """Run the per-stage formulate-and-solve loop in place on ``mTEPES``. + + ``mTEPES`` must already carry the objective, the investment constraints, ``First_st`` / ``Last_st`` and an + initialised ``pDuals`` dict (all set up by ``openTEPES_run`` before this call). ``_path`` is the case path + used for optional LP-file dumps; ``pIndCycleFlow`` toggles the cycle-flow network formulation. + """ + # initialize the set of load levels up to the current stage + mTEPES.na = Set(initialize=[]) + + # iterative model formulation for each stage of a year + for p,sc,st in mTEPES.ps*mTEPES.stt: + + # control for not repeating the stages in case of several ones + mTEPES.NoRepetition = 0 + + # activate only load levels to formulate + mTEPES.del_component(mTEPES.st) + mTEPES.del_component(mTEPES.n ) + mTEPES.del_component(mTEPES.n2) + mTEPES.st = Set(doc='stages', initialize=[stt for stt in mTEPES.stt if st == stt if mTEPES.pStageWeight[stt] and sum(1 for (p,sc,st,nn) in mTEPES.s2n)]) + mTEPES.n = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if (p,sc,st,nn) in mTEPES.s2n ]) + mTEPES.n2 = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if (p,sc,st,nn) in mTEPES.s2n ]) + + # load levels multiple of cycles for each ESS/generator + mTEPES.nesc = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [es] == 0] + mTEPES.necc = [(n,ec) for n,ec in mTEPES.n*mTEPES.ec if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [ec] == 0] + mTEPES.neso = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pOutflowsTimeStep[es] == 0] + mTEPES.ngen = [(n,g ) for n,g in mTEPES.n*mTEPES.g if mTEPES.n.ord(n) % mTEPES.pEnergyTimeStep [g ] == 0] + if mTEPES.pIndHydroTopology: + mTEPES.nhc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2h) == 0] + if sum(1 for h,rs in mTEPES.p2r): + mTEPES.np2c = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) == 0] + else: + mTEPES.np2c = [] + if sum(1 for rs,h in mTEPES.r2p): + mTEPES.npc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) == 0] + else: + mTEPES.npc = [] + mTEPES.nrsc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] + mTEPES.nrcc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rn if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] + mTEPES.nrso = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pWaterOutTimeStep [rs] == 0] + + if mTEPES.st: + + # load levels up to the current stage for emission and RES energy constraints + # if (p != mTEPES.p.first() or sc != mTEPES.sc.first()) and st == mTEPES.First_st: + # mTEPES.del_component(mTEPES.na) + if st == mTEPES.First_st: + mTEPES.del_component(mTEPES.na) + mTEPES.na = Set(initialize=mTEPES.n) + else: + temp_na = mTEPES.na | mTEPES.n # Create a temporary set + mTEPES.del_component(mTEPES.na) # Delete the existing set to avoid overwriting + mTEPES.na = Set(initialize=temp_na) # Assign the new set + + print(f'Period {p}, Scenario {sc}, Stage {st}') + + # operation model objective function and constraints by stage + GenerationOperationModelFormulationObjFunct (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + GenerationOperationElecModelFormulationInvestment (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + if mTEPES.pIndHeat: + GenerationOperationHeatModelFormulationInvestment(mTEPES, mTEPES, pIndLogConsole, p, sc, st) + GenerationOperationModelFormulationDemand (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + GenerationOperationModelFormulationStorage (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + if mTEPES.pIndHydroTopology: + GenerationOperationModelFormulationReservoir (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + if mTEPES.pIndHydrogen: + NetworkH2OperationModelFormulation (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + if mTEPES.pIndHeat: + NetworkHeatOperationModelFormulation (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + GenerationOperationModelFormulationCommitment (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + GenerationOperationModelFormulationRampMinTime (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + NetworkSwitchingModelFormulation (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + NetworkOperationModelFormulation (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + # introduce cycle flow formulations + if pIndCycleFlow: + if st == mTEPES.First_st: + NetworkCycles ( mTEPES, pIndLogConsole ) + CycleConstraints (mTEPES, mTEPES, pIndLogConsole, p, sc, st) + + # No generator investments and no generator retirements, and no line investments + if ( (sum(1 for pp,eb in mTEPES.peb if pp <= p) == 0 or mTEPES.pIndBinGenInvest() == 2) + and (sum(1 for pp,gd in mTEPES.pgd if pp <= p) == 0 or mTEPES.pIndBinGenRetire() == 2) + and (sum(1 for pp,rc in mTEPES.prc if pp <= p) == 0 or mTEPES.pIndBinRsrInvest() == 2) + and (sum(1 for pp,ni,nf,cc in mTEPES.plc if pp <= p) == 0 or mTEPES.pIndBinNetElecInvest() == 2) + and (sum(1 for pp,ni,nf,cc in mTEPES.ppc if pp <= p) == 0 or mTEPES.pIndBinNetH2Invest() == 2) + and (sum(1 for pp,ni,nf,cc in mTEPES.phc if pp <= p) == 0 or mTEPES.pIndBinNetHeatInvest() == 2)): + + # No minimum RES requirements and no emission limit + if (max([mTEPES.pRESEnergy[p,ar] for ar in mTEPES.ar]) == 0 and (min([mTEPES.pEmission [p,ar] for ar in mTEPES.ar]) == math.inf or sum(mTEPES.pEmissionRate[nr] for nr in mTEPES.nr) == 0)): + + if (p,sc) == mTEPES.ps.last() and st == mTEPES.Last_st and mTEPES.NoRepetition == 0: + mTEPES.NoRepetition = 1 + + # Writing LP file + if pIndLogConsole: + StartTime = time.time() + mTEPES.write(f'{_path}/openTEPES_{CaseName}_{p}_{sc}_{st}.lp', io_options={'symbolic_solver_labels': True}) + WritingLPFileTime = time.time() - StartTime + StartTime = time.time() + print('Writing LP file ... ', round(WritingLPFileTime), 's') + + # there are no expansion decisions, or they are ignored (it is an operation model), or there are no system emission nor minimum RES constraints + mTEPES.pScenProb[p,sc] = 1.0 + ProblemSolving(DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, mTEPES.p.ord(p)*mTEPES.sc.ord(sc)*mTEPES.st.ord(st)) + mTEPES.pScenProb[p,sc] = 0.0 + + # deactivate the constraints of the previous period and scenario + for c in mTEPES.component_objects(pyo.Constraint, active=True): + if c.name.find(str(p)) != -1 and c.name.find(str(sc)) != -1: + c.deactivate() + + # Minimum RES requirements or emission limit + elif (max([mTEPES.pRESEnergy[p,ar] for ar in mTEPES.ar]) > 0 or (min([mTEPES.pEmission [p,ar] for ar in mTEPES.ar]) < math.inf and sum(mTEPES.pEmissionRate[nr] for nr in mTEPES.nr) > 0)): + + if st == mTEPES.Last_st and mTEPES.NoRepetition == 0: + + if (p,sc) == mTEPES.ps.last(): + mTEPES.NoRepetition = 1 + + mTEPES.del_component(mTEPES.st) + mTEPES.del_component(mTEPES.n ) + mTEPES.del_component(mTEPES.n2) + mTEPES.st = Set(doc='stages', initialize=[stt for stt in mTEPES.stt if mTEPES.pStageWeight[stt] and sum(1 for p,sc,stt,nn in mTEPES.s2n)]) + mTEPES.n = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) + mTEPES.n2 = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) + + # load levels multiple of cycles for each ESS/generator + mTEPES.nesc = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [es] == 0] + mTEPES.necc = [(n,ec) for n,ec in mTEPES.n*mTEPES.ec if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [ec] == 0] + mTEPES.neso = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pOutflowsTimeStep[es] == 0] + mTEPES.ngen = [(n,g ) for n,g in mTEPES.n*mTEPES.g if mTEPES.n.ord(n) % mTEPES.pEnergyTimeStep [g ] == 0] + if mTEPES.pIndHydroTopology: + mTEPES.nhc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2h) == 0] + if sum(1 for h,rs in mTEPES.p2r): + mTEPES.np2c = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) == 0] + else: + mTEPES.np2c = [] + if sum(1 for rs,h in mTEPES.r2p): + mTEPES.npc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) == 0] + else: + mTEPES.npc = [] + mTEPES.nrsc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] + mTEPES.nrcc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rn if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] + mTEPES.nrso = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pWaterOutTimeStep [rs] == 0] + + # Writing an LP file + if pIndLogConsole: + StartTime = time.time() + mTEPES.write(f'{_path}/openTEPES_{CaseName}_{p}_{sc}_{st}.lp', io_options={'symbolic_solver_labels': True}) + WritingLPFileTime = time.time() - StartTime + StartTime = time.time() + print('Writing LP file ... ', round(WritingLPFileTime), 's') + + mTEPES.pScenProb[p,sc] = 1.0 + # there are system emission or RES requirement constraints + ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, mTEPES.p.ord(p)*mTEPES.sc.ord(sc)*mTEPES.st.ord(st)) + # if pIndSectorDecomposition and len(mTEPES.psnel): + # ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, 1) + # SectorDecomposition(DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st ) + # else: + # ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, mTEPES.p.ord(p)*mTEPES.sc.ord(sc)*mTEPES.st.ord(st)) + mTEPES.pScenProb[p,sc] = 0.0 + + # deactivate the constraints of the previous period and scenario + for c in mTEPES.component_objects(pyo.Constraint, active=True): + if c.name.find(str(p)) != -1 and c.name.find(str(sc)) != -1: + c.deactivate() + + # generator investments or generator retirements, or line investments + elif ( (sum(1 for pp,eb in mTEPES.peb if pp <= p) > 0 and mTEPES.pIndBinGenInvest() < 2) + or (sum(1 for pp,gd in mTEPES.pgd if pp <= p) > 0 and mTEPES.pIndBinGenRetire() < 2) + or (sum(1 for pp,rc in mTEPES.prc if pp <= p) > 0 and mTEPES.pIndBinRsrInvest() < 2) + or (sum(1 for pp,ni,nf,cc in mTEPES.plc if pp <= p) > 0 and mTEPES.pIndBinNetElecInvest() < 2) + or (sum(1 for pp,ni,nf,cc in mTEPES.ppc if pp <= p) > 0 and mTEPES.pIndBinNetH2Invest() < 2) + or (sum(1 for pp,ni,nf,cc in mTEPES.phc if pp <= p) > 0 and mTEPES.pIndBinNetHeatInvest() < 2)): + + if (p,sc) == mTEPES.ps.last() and st == mTEPES.Last_st and mTEPES.NoRepetition == 0: + mTEPES.NoRepetition = 1 + + mTEPES.del_component(mTEPES.st) + mTEPES.del_component(mTEPES.n ) + mTEPES.del_component(mTEPES.n2) + mTEPES.st = Set(doc='stages', initialize=[stt for stt in mTEPES.stt if mTEPES.pStageWeight[stt] and sum(1 for p,sc,stt,nn in mTEPES.s2n)]) + mTEPES.n = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) + mTEPES.n2 = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) + + # load levels multiple of cycles for each ESS/generator + mTEPES.nesc = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [es] == 0] + mTEPES.necc = [(n,ec) for n,ec in mTEPES.n*mTEPES.ec if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [ec] == 0] + mTEPES.neso = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pOutflowsTimeStep[es] == 0] + mTEPES.ngen = [(n,g ) for n,g in mTEPES.n*mTEPES.g if mTEPES.n.ord(n) % mTEPES.pEnergyTimeStep [g ] == 0] + if mTEPES.pIndHydroTopology: + mTEPES.nhc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2h) == 0] + if sum(1 for h,rs in mTEPES.p2r): + mTEPES.np2c = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) == 0] + else: + mTEPES.np2c = [] + if sum(1 for rs,h in mTEPES.r2p): + mTEPES.npc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) == 0] + else: + mTEPES.npc = [] + mTEPES.nrsc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] + mTEPES.nrcc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rn if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] + mTEPES.nrso = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pWaterOutTimeStep [rs] == 0] + + # Writing an LP file + if pIndLogConsole: + StartTime = time.time() + mTEPES.write(f'{_path}/openTEPES_{CaseName}_{p}_{sc}_{st}.lp', io_options={'symbolic_solver_labels': True}) + WritingLPFileTime = time.time() - StartTime + StartTime = time.time() + print('Writing LP file ... ', round(WritingLPFileTime), 's') + + # there are investment decisions (it is an expansion and operation model), or there are system emission constraints + ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, 1) + # if pIndSectorDecomposition and len(mTEPES.psnel): + # ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, 1) + # SectorDecomposition(DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st ) + # else: + # ProblemSolving (DirName, CaseName, SolverName, mTEPES, mTEPES, pIndLogConsole, p, sc, st, 1) + + mTEPES.del_component(mTEPES.st) + mTEPES.del_component(mTEPES.n ) + mTEPES.del_component(mTEPES.n2) + mTEPES.st = Set(doc='stages', initialize=[stt for stt in mTEPES.stt if mTEPES.pStageWeight[stt] and sum(1 for p,sc,stt,nn in mTEPES.s2n)]) + mTEPES.n = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) + mTEPES.n2 = Set(doc='load levels', initialize=[nn for nn in mTEPES.nn if sum(1 for p,sc,st in mTEPES.ps*mTEPES.st if (p,sc,st, nn) in mTEPES.s2n)]) + + # load levels multiple of cycles for each ESS/generator + mTEPES.nesc = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [es] == 0] + mTEPES.necc = [(n,ec) for n,ec in mTEPES.n*mTEPES.ec if mTEPES.n.ord(n) % mTEPES.pStorageTimeStep [ec] == 0] + mTEPES.neso = [(n,es) for n,es in mTEPES.n*mTEPES.es if mTEPES.n.ord(n) % mTEPES.pOutflowsTimeStep[es] == 0] + mTEPES.ngen = [(n,g ) for n,g in mTEPES.n*mTEPES.g if mTEPES.n.ord(n) % mTEPES.pEnergyTimeStep [g ] == 0] + if mTEPES.pIndHydroTopology: + mTEPES.nhc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2h) == 0] + if sum(1 for h,rs in mTEPES.p2r): + mTEPES.np2c = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (h,rs) in mTEPES.p2r) == 0] + else: + mTEPES.np2c = [] + if sum(1 for rs,h in mTEPES.r2p): + mTEPES.npc = [(n,h ) for n,h in mTEPES.n*mTEPES.h if sum(1 for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) and mTEPES.n.ord(n) % sum(mTEPES.pReservoirTimeStep[rs] for rs in mTEPES.rs if (rs,h) in mTEPES.r2p) == 0] + else: + mTEPES.npc = [] + mTEPES.nrsc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] + mTEPES.nrcc = [(n,rs) for n,rs in mTEPES.n*mTEPES.rn if mTEPES.n.ord(n) % mTEPES.pReservoirTimeStep[rs] == 0] + mTEPES.nrso = [(n,rs) for n,rs in mTEPES.n*mTEPES.rs if mTEPES.n.ord(n) % mTEPES.pWaterOutTimeStep [rs] == 0] + + # activate the constraints of all the periods and scenarios + for c in mTEPES.component_objects(pyo.Constraint): + c.activate() + + # Output convention follows the two intended modes: + # - Case 2 (no expansion decisions): the per-scenario solve loop above + # used pScenProb=1.0 for each scenario individually. Reset all + # probabilities to 1.0 so cost-summary writers aggregate the + # independent deterministic solves as a sum. + # - Case 1 (with expansion decisions): the joint solve at the last + # (p,sc) used the input probabilities to minimise probability- + # weighted expected cost. Preserve those probabilities so the + # output reports that same expected cost — otherwise the writers + # would over-count by ~N_scenarios. + _hasExpansion = ( + (sum(1 for _ in mTEPES.peb) > 0 and mTEPES.pIndBinGenInvest() < 2) + or (sum(1 for _ in mTEPES.pgd) > 0 and mTEPES.pIndBinGenRetire() < 2) + or (sum(1 for _ in mTEPES.prc) > 0 and mTEPES.pIndBinRsrInvest() < 2) + or (sum(1 for _ in mTEPES.plc) > 0 and mTEPES.pIndBinNetElecInvest() < 2) + or (sum(1 for _ in mTEPES.ppc) > 0 and mTEPES.pIndBinNetH2Invest() < 2) + or (sum(1 for _ in mTEPES.phc) > 0 and mTEPES.pIndBinNetHeatInvest() < 2) + ) + if not _hasExpansion: + # assign probability 1 to all the periods and scenarios + for p,sc in mTEPES.ps: + mTEPES.pScenProb[p,sc] = 1.0 diff --git a/tests/test_direct_run.py b/tests/test_direct_run.py new file mode 100644 index 00000000..6247fa70 --- /dev/null +++ b/tests/test_direct_run.py @@ -0,0 +1,38 @@ +"""Each package module must import cleanly when run as a script. + +When a file is run directly (VS Code "Run Python File", ``python path/to/file.py``) it executes as ``__main__`` +with an empty ``__package__``, so a top-level ``from .openTEPES_X import ...`` relative import fails. The modules +that import sibling modules at load time guard against this with a ``try``/``except ImportError`` fallback that +puts the repository root on ``sys.path`` and retries the import absolutely. + +Nothing else in the test suite exercises that fallback: a normal ``import openTEPES`` only ever takes the +relative-import path. This test runs every ``openTEPES_*.py`` module as ``__main__`` (via ``runpy``) so the +fallback path is checked on every CI run, catching a broken guard, a missing import, or a circular import before +a user hits it in an editor. No model is solved, so it runs in the fast (non-solve) CI job. + +``openTEPES_Main.py`` is excluded because it has an active ``if __name__ == "__main__"`` block that would launch +the command-line interface rather than just import. +""" +import runpy +from pathlib import Path + +import pytest + +PKG_DIR = Path(__file__).resolve().parent.parent / "openTEPES" + +# Modules with an active __main__ block that does real work (not just imports) — running them as __main__ would +# execute that work, not test the import guard. +EXCLUDE = {"openTEPES_Main.py"} + +MODULES = sorted(p.name for p in PKG_DIR.glob("openTEPES_*.py") if p.name not in EXCLUDE) + + +def test_modules_discovered(): + """Guard against the glob silently matching nothing (e.g. a path change).""" + assert MODULES, f"no openTEPES_*.py modules found under {PKG_DIR}" + + +@pytest.mark.parametrize("module", MODULES) +def test_module_runs_as_main(module): + """Running the module as a script must not raise (the relative-import guard's fallback works).""" + runpy.run_path(str(PKG_DIR / module), run_name="__main__")