diff --git a/openTEPES/openTEPES.py b/openTEPES/openTEPES.py index e9c23895..1ea6de26 100644 --- a/openTEPES/openTEPES.py +++ b/openTEPES/openTEPES.py @@ -9,387 +9,350 @@ import time import pyomo.environ as pyo -from pyomo.environ import ConcreteModel, Set - -from .openTEPES_InputData import InputData, DataConfiguration, SettingUpVariables -from .openTEPES_ModelFormulation import TotalObjectiveFunction, InvestmentModelFormulation, GenerationOperationModelFormulationObjFunct, GenerationOperationModelFormulationInvestment, GenerationOperationModelFormulationDemand, GenerationOperationModelFormulationStorage, GenerationOperationModelFormulationReservoir, NetworkH2OperationModelFormulation, NetworkHeatOperationModelFormulation, GenerationOperationModelFormulationCommitment, GenerationOperationModelFormulationRampMinTime, NetworkSwitchingModelFormulation, NetworkOperationModelFormulation, NetworkCycles, CycleConstraints -from .openTEPES_ProblemSolving import ProblemSolving -from .openTEPES_OutputResults import OutputResultsParVarCon, InvestmentResults, GenerationOperationResults, GenerationOperationHeatResults, ESSOperationResults, ReservoirOperationResults, NetworkH2OperationResults, NetworkHeatOperationResults, FlexibilityResults, NetworkOperationResults, MarginalResults, OperationSummaryResults, ReliabilityResults, CostSummaryResults, EconomicResults, NetworkMapResults +from pyomo.environ import ConcreteModel, Set + + +def _assert_single_objective_active(model): + """Defensive check: ensure exactly one active Objective before solve.""" + objs = [o for o in model.component_objects(pyo.Objective, descend_into=True) if o.active] + if len(objs) != 1: + names = [o.name for o in objs] + raise RuntimeError(f"Expected exactly one active Objective, found {len(objs)}: {names}") + + +from .openTEPES_InputData import InputData, DataConfiguration, SettingUpVariables +from .openTEPES_ModelFormulation import ( + TotalObjectiveFunction, + InvestmentModelFormulation, + GenerationOperationModelFormulationObjFunct, + GenerationOperationModelFormulationInvestment, + GenerationOperationModelFormulationDemand, + GenerationOperationModelFormulationStorage, + GenerationOperationModelFormulationReservoir, + NetworkH2OperationModelFormulation, + NetworkHeatOperationModelFormulation, + GenerationOperationModelFormulationCommitment, + GenerationOperationModelFormulationRampMinTime, + NetworkSwitchingModelFormulation, + NetworkOperationModelFormulation, + NetworkCycles, + CycleConstraints +) +from .openTEPES_ProblemSolving import ProblemSolving +from .openTEPES_OutputResults import ( + OutputResultsParVarCon, + InvestmentResults, + GenerationOperationResults, + GenerationOperationHeatResults, + ESSOperationResults, + ReservoirOperationResults, + NetworkH2OperationResults, + NetworkHeatOperationResults, + FlexibilityResults, + NetworkOperationResults, + MarginalResults, + OperationSummaryResults, + ReliabilityResults, + CostSummaryResults, + EconomicResults, + NetworkMapResults +) + +# Static lists of formulation and output functions +_STAGE_FORMULATIONS = [ + ('pIndGenerationOperationModelFormulationObjFunct', GenerationOperationModelFormulationObjFunct), + ('pIndGenerationOperationModelFormulationInvestment', GenerationOperationModelFormulationInvestment), + ('pIndGenerationOperationModelFormulationDemand', GenerationOperationModelFormulationDemand), + ('pIndGenerationOperationModelFormulationStorage', GenerationOperationModelFormulationStorage), + ('pIndGenerationOperationModelFormulationReservoir', GenerationOperationModelFormulationReservoir), + ('pIndNetworkH2OperationModelFormulation', NetworkH2OperationModelFormulation), + ('pIndNetworkHeatOperationModelFormulation', NetworkHeatOperationModelFormulation), + ('pIndGenerationOperationModelFormulationCommitment', GenerationOperationModelFormulationCommitment), + ('pIndGenerationOperationModelFormulationRampMinTime', GenerationOperationModelFormulationRampMinTime), + ('pIndNetworkSwitchingModelFormulation', NetworkSwitchingModelFormulation), + ('pIndNetworkOperationModelFormulation', NetworkOperationModelFormulation), + ('pIndNetworkCycles', NetworkCycles), + ('pIndCycleConstraints', CycleConstraints) +] + +_OUTPUT_FUNCS = [ + ('pIndDumpRawResults', OutputResultsParVarCon, ('DirName', 'CaseName', 'mTEPES', 'mTEPES')), + ('pIndInvestmentResults', InvestmentResults, + ('DirName', 'CaseName', 'mTEPES', 'mTEPES', 'pIndTechnologyOutput', 'pIndPlotOutput')), + ('pIndGenerationOperationResults', GenerationOperationResults, + ('DirName', 'CaseName', 'mTEPES', 'mTEPES', 'pIndTechnologyOutput', 'pIndAreaOutput', 'pIndPlotOutput')), + ('pIndESSOperationResults', ESSOperationResults, + ('DirName', 'CaseName', 'mTEPES', 'mTEPES', 'pIndTechnologyOutput', 'pIndAreaOutput', 'pIndPlotOutput')), + ('pIndReservoirOperationResults', ReservoirOperationResults, + ('DirName', 'CaseName', 'mTEPES', 'mTEPES', 'pIndTechnologyOutput', 'pIndPlotOutput')), + ('pIndNetworkH2OperationResults', NetworkH2OperationResults, ('DirName', 'CaseName', 'mTEPES', 'mTEPES')), + ('pIndNetworkHeatOperationResults', NetworkHeatOperationResults, ('DirName', 'CaseName', 'mTEPES', 'mTEPES')), + ('pIndFlexibilityResults', FlexibilityResults, ('DirName', 'CaseName', 'mTEPES', 'mTEPES')), + ('pIndReliabilityResults', ReliabilityResults, ('DirName', 'CaseName', 'mTEPES', 'mTEPES')), + ('pIndNetworkOperationResults', NetworkOperationResults, ('DirName', 'CaseName', 'mTEPES', 'mTEPES')), + ('pIndNetworkMapResults', NetworkMapResults, ('DirName', 'CaseName', 'mTEPES', 'mTEPES')), + ('pIndOperationSummaryResults', OperationSummaryResults, ('DirName', 'CaseName', 'mTEPES', 'mTEPES')), + ('pIndCostSummaryResults', CostSummaryResults, ('DirName', 'CaseName', 'mTEPES', 'mTEPES')), + ('pIndMarginalResults', MarginalResults, ('DirName', 'CaseName', 'mTEPES', 'mTEPES', 'pIndPlotOutput')), + ('pIndEconomicResults', EconomicResults, + ('DirName', 'CaseName', 'mTEPES', 'mTEPES', 'pIndAreaOutput', 'pIndPlotOutput')) +] + + +def _initialize_indices(): + zero_map = dict(zip([0, 0.0, 'No', 'NO', 'no', 'N', 'n'], [0] * 7)) + one_map = dict(zip(['Yes', 'YES', 'yes', 'Y', 'y'], [1] * 5)) + return {**zero_map, **one_map} + + +def _build_model(desc: str) -> ConcreteModel: + model = ConcreteModel() + model.name = desc + return model + + +def _reading_data(dir, case, model, log_flag): + InputData(dir, case, model, log_flag) + + +def _configure_basic_components(dir_, case, model, log_flag): + DataConfiguration(model) + SettingUpVariables(model, model) + model.First_st, model.Last_st = next(iter(model.st)), None + for st in model.st: + model.Last_st = st + model.pDuals, model.na = {}, Set(initialize=[]) + model.NoRepetition = 0 + model._has_run_last_stage = False + + +def _configure_formulation_flags(model): + flags = dict( + pIndGenerationOperationModelFormulationObjFunct=1, + pIndGenerationOperationModelFormulationInvestment=1, + pIndGenerationOperationModelFormulationDemand=1, + pIndGenerationOperationModelFormulationStorage=1, + pIndGenerationOperationModelFormulationReservoir=1, + pIndNetworkH2OperationModelFormulation=1, + pIndNetworkHeatOperationModelFormulation=1, + pIndGenerationOperationModelulationCommitment=0, + pIndGenerationOperationModelulationRampMinTime=0, + pIndNetworkSwitchingModelulation=0, + pIndNetworkOperationModelulation=0, + pIndNetworkCycles=0, + pIndCycleConstraints=0 + ) + model._formulation_flags = flags + return flags + + +def _configure_output_flags(output_flag: int, init_flags: int) -> dict: + """Generate the output control flags dict.""" + base = dict( + pIndTechnologyOutput=2, + pIndAreaOutput=1, + pIndPlotOutput=0 + ) + detailed = {key: int(output_flag == init_flags) for key, _func, _ in _OUTPUT_FUNCS} + return {**base, **detailed} + + +def _check_conditions(model, flag_name: str) -> bool: + enabled = model._formulation_flags.get(flag_name, 1) + cond_map = { + 'pIndGenerationOperationResults': True, + 'pIndNetworkH2OperationResults': getattr(model, 'pa', False) and getattr(model, 'pIndHydrogen', 0) == 1, + 'pIndNetworkHeatOperationResults': getattr(model, 'ha', False) and getattr(model, 'pIndHeat', 0) == 1, + 'pIndReservoirOperationResults': getattr(model, 'rs', False) and getattr(model, 'pIndHydroTopology', 0) == 1, + 'pIndESSOperationResults': getattr(model, 'es', False), + 'pIndGenerationOperationModelFormulationReservoir': getattr(model, 'rs', False) and getattr(model, + 'pIndHydroTopology', + 0) == 1, + 'pIndNetworkH2OperationModelFormulation': getattr(model, 'pa', False) and getattr(model, 'pIndHydrogen', + 0) == 1, + 'pIndNetworkHeatOperationModelFormulation': getattr(model, 'ha', False) and getattr(model, 'pIndHeat', 0) == 1, + 'pIndNetworkCycles': getattr(model, 'pIndCycleFlow', 0) == 1, + 'pIndCycleConstraints': getattr(model, 'pIndCycleFlow', 0) == 1, + } + return bool(enabled and cond_map.get(flag_name, True)) + + +def _rebuild_stage_sets(model, p, sc, st): + model.del_component(model.st) + model.del_component(model.n) + model.del_component(model.n2) + model.st = Set(doc='stages', initialize=[ + stt for stt in model.stt + if stt == st + and model.pStageWeight[stt] + and any((p, sc, stt, nn) in model.s2n for nn in model.nn) + ]) + model.n = Set(doc='load levels', initialize=[ + nn for nn in model.nn + if (p, sc, st, nn) in model.s2n + ]) + model.n2 = Set(doc='load levels', initialize=model.n) + + +def _compute_cycle_lists(model): + model.nesc = [(n, es) for n, es in model.n * model.es if model.n.ord(n) % model.pStorageTimeStep[es] == 0] + model.necc = [(n, ec) for n, ec in model.n * model.ec if model.n.ord(n) % model.pStorageTimeStep[ec] == 0] + model.neso = [(n, es) for n, es in model.n * model.es if model.n.ord(n) % model.pOutflowsTimeStep[es] == 0] + model.ngen = [(n, g) for n, g in model.n * model.g if model.n.ord(n) % model.pEnergyTimeStep[g] == 0] + if getattr(model, 'pIndHydroTopology', 0) == 1: + model.nhc = [(n, h) for n, h in model.n * model.h + if + model.n.ord(n) % sum(model.pReservoirTimeStep[rs] for rs in model.rs if (rs, h) in model.r2h) == 0] + model.np2c = [(n, h) for n, h in model.n * model.h + if any((h, rs) in model.p2r for rs in model.rs) + and model.n.ord(n) % sum( + model.pReservoirTimeStep[rs] for rs in model.rs if (h, rs) in model.p2r) == 0] + model.npc = [(n, h) for n, h in model.n * model.h + if any((rs, h) in model.r2p for rs in model.rs) + and model.n.ord(n) % sum( + model.pReservoirTimeStep[rs] for rs in model.rs if (rs, h) in model.r2p) == 0] + model.nrsc = [(n, rs) for n, rs in model.n * model.rs if model.n.ord(n) % model.pReservoirTimeStep[rs] == 0] + model.nrcc = [(n, rs) for n, rs in model.n * model.rn if model.n.ord(n) % model.pReservoirTimeStep[rs] == 0] + model.nrso = [(n, rs) for n, rs in model.n * model.rs if model.n.ord(n) % model.pWaterOutTimeStep[rs] == 0] + + +def _accumulate_na(model, st): + if st == model.First_st: + model.del_component(model.na) + model.na = Set(initialize=model.n) + else: + temp = set(model.na) | set(model.n) + model.del_component(model.na) + model.na = Set(initialize=temp) + + +def _apply_formulations(model, log_flag, p, sc, st): + for flag_name, func in _STAGE_FORMULATIONS: + if _check_conditions(model, flag_name): + func(model, model, log_flag, p, sc, st) + + +def _solve_branch(dir_, case, solver, model, log_flag, p, sc, st): + no_invest = ( + (sum(1 for rec in model.peb if rec[0] <= p) == 0 or model.pIndBinGenInvest() == 2) and + (sum(1 for rec in model.pgd if rec[0] <= p) == 0 or model.pIndBinGenRetire() == 2) and + (sum(1 for rec in model.prc if rec[0] <= p) == 0 or model.pIndBinRsrInvest() == 2) and + (sum(1 for rec in model.plc if rec[0] <= p) == 0 or model.pIndBinNetElecInvest() == 2) and + (sum(1 for rec in model.ppc if rec[0] <= p) == 0 or model.pIndBinNetH2Invest() == 2) and + (sum(1 for rec in model.phc if rec[0] <= p) == 0 or model.pIndBinNetHeatInvest() == 2) + ) + min_res = ( + max(model.pRESEnergy[p, ar] for ar in model.ar) > 0 or + (min(model.pEmission[p, ar] for ar in model.ar) < math.inf and sum( + model.pEmissionRate[nr] for nr in model.nr) > 0) + ) + no_res = ( + max(model.pRESEnergy[p, ar] for ar in model.ar) == 0 and + (min(model.pEmission[p, ar] for ar in model.ar) == math.inf or sum( + model.pEmissionRate[nr] for nr in model.nr) == 0) + ) + + def write_lp(): + if log_flag == 1: + t0 = time.time() + model.write(f"{dir_}/{case}/openTEPES_{case}_{p}_{sc}_{st}.lp", + io_options={'symbolic_solver_labels': True}) + print('Writing LP file ...', round(time.time() - t0), 's') + + if no_invest and no_res: + write_lp() + model.pScenProb[p, sc] = 1.0 + ProblemSolving(dir_, case, solver, model, model, log_flag, p, sc, st) + model.pScenProb[p, sc] = 0.0 + elif no_invest and min_res: + write_lp() + ProblemSolving(dir_, case, solver, model, model, log_flag, p, sc, st) + else: + write_lp() + ProblemSolving(dir_, case, solver, model, model, log_flag, p, sc, st) + + for c in model.component_objects(pyo.Constraint, active=True): + nm = c.name + if all(tag in nm for tag in (str(p), str(sc), str(st))): + c.deactivate() + + +def _process_stage(dir_, case, solver, model, log_flag, p, sc, st): + _rebuild_stage_sets(model, p, sc, st) + _compute_cycle_lists(model) + _accumulate_na(model, st) + _apply_formulations(model, log_flag, p, sc, st) + _solve_branch(dir_, case, solver, model, log_flag, p, sc, st) + + if st == model.Last_st and not model._has_run_last_stage: + model.NoRepetition = 1 + model._has_run_last_stage = True + + +def process_stage_loop(dir_, case, solver, model, log_flag): + for p, sc in model.ps: + model._has_run_last_stage = False + for st in model.stt: + _process_stage(dir_, case, solver, model, log_flag, p, sc, st) + + +def finalize_and_output(dir_, case, model, flags): + """Run output routines based on computed flags.""" + context = { + 'DirName': dir_, + 'CaseName': case, + 'OptModel': model, + 'mTEPES': model, + 'pIndTechnologyOutput': flags['pIndTechnologyOutput'], + 'pIndAreaOutput': flags['pIndAreaOutput'], + 'pIndPlotOutput': flags['pIndPlotOutput'] + } + + for flag_name, func, args in _OUTPUT_FUNCS: + if not (flags.get(flag_name) and _check_conditions(model, flag_name)): + continue + # build a list of actual values, not raw strings + r_args = [context[a] if isinstance(a, str) else a for a in args] + # unpack them into separate parameters + func(*r_args) def openTEPES_run(DirName, CaseName, SolverName, pIndOutputResults, pIndLogConsole): - - InitialTime = time.time() - _path = os.path.join(DirName, CaseName) - - #%% replacing string values by numerical values - idxDict = dict() - idxDict[0 ] = 0 - idxDict[0.0 ] = 0 - idxDict[0.0 ] = 0 - idxDict['No' ] = 0 - idxDict['NO' ] = 0 - idxDict['no' ] = 0 - idxDict['N' ] = 0 - idxDict['n' ] = 0 - idxDict['Yes'] = 1 - idxDict['YES'] = 1 - idxDict['yes'] = 1 - idxDict['Y' ] = 1 - 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.6 - July 22, 2025') - print( 'Open Generation, Storage, and Transmission Operation and Expansion Planning Model with RES and ESS (openTEPES) - Version 4.18.6 - July 22, 2025', file=open(f'{_path}/openTEPES_version_{CaseName}.log','w')) - - pIndOutputResults = [j for i,j in idxDict.items() if i == pIndOutputResults][0] - pIndLogConsole = [j for i,j in idxDict.items() if i == pIndLogConsole ][0] - - # introduce cycle flow formulations in DC and AC load flow - pIndCycleFlow = 0 - - # Reading sets and parameters - InputData(DirName, CaseName, mTEPES, pIndLogConsole) - - # Define sets and parameters - DataConfiguration(mTEPES) - - # Define variables - SettingUpVariables(mTEPES, mTEPES) - - # first/last stage - FirstST = 0 - for st in mTEPES.st: - if FirstST == 0: - FirstST = 1 - mTEPES.First_st = st - mTEPES.Last_st = st - - # objective function and investment constraints - TotalObjectiveFunction (mTEPES, mTEPES, pIndLogConsole) - InvestmentModelFormulation(mTEPES, mTEPES, pIndLogConsole) - - # 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 == 1: - 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) - GenerationOperationModelFormulationInvestment (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - GenerationOperationModelFormulationDemand (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - GenerationOperationModelFormulationStorage (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - if mTEPES.pIndHydroTopology == 1: - GenerationOperationModelFormulationReservoir(mTEPES, mTEPES, pIndLogConsole, p, sc, st) - if mTEPES.pIndHydrogen == 1: - NetworkH2OperationModelFormulation (mTEPES, mTEPES, pIndLogConsole, p, sc, st) - if mTEPES.pIndHeat == 1: - 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 == 1: - 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,lc in mTEPES.plc if pp <= p) == 0 or mTEPES.pIndBinNetElecInvest() == 2 - and sum(1 for pp,pc in mTEPES.ppc if pp <= p) == 0 or mTEPES.pIndBinNetH2Invest() == 2 - and sum(1 for pp,hc 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 == 1: - 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.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 == 1: - 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 == 1: - 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) - - # 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,lc in mTEPES.plc if pp <= p) > 0 and mTEPES.pIndBinNetElecInvest() < 2 - or sum(1 for pp,pc in mTEPES.ppc if pp <= p) > 0 and mTEPES.pIndBinNetH2Invest() < 2 - or sum(1 for pp,hc 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 == 1: - 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 == 1: - 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) - - 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 == 1: - 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() - - # assign probability 1 to all the periods and scenarios - for p,sc in mTEPES.ps: - mTEPES.pScenProb[p,sc] = 1.0 - - # pickle the case study data - # with open(dump_folder+f'/oT_Case_{CaseName}.pkl','wb') as f: - # pickle.dump(mTEPES, f, pickle.HIGHEST_PROTOCOL) - - # output results only for every unit (0), only for every technology (1), or for both (2) - pIndTechnologyOutput = 2 - - # output results just for the system (0) or for every area (1). Areas correspond usually to countries - pIndAreaOutput = 1 - - # output plot results - pIndPlotOutput = 1 - - # indicators to control the number of output results - if pIndOutputResults == 1: - pIndDumpRawResults = 0 - pIndInvestmentResults = 1 - pIndGenerationOperationResults = 1 - pIndESSOperationResults = 1 - pIndReservoirOperationResults = 1 - pIndNetworkH2OperationResults = 1 - pIndNetworkHeatOperationResults = 1 - pIndFlexibilityResults = 1 - pIndReliabilityResults = 1 - pIndNetworkOperationResults = 1 - pIndNetworkMapResults = 1 - pIndOperationSummaryResults = 1 - pIndCostSummaryResults = 1 - pIndMarginalResults = 1 - pIndEconomicResults = 1 - else: - pIndDumpRawResults = 0 - pIndInvestmentResults = 1 - pIndGenerationOperationResults = 0 - pIndESSOperationResults = 0 - pIndReservoirOperationResults = 0 - pIndNetworkH2OperationResults = 0 - pIndNetworkHeatOperationResults = 0 - pIndFlexibilityResults = 0 - pIndReliabilityResults = 0 - pIndNetworkOperationResults = 0 - pIndNetworkMapResults = 0 - pIndOperationSummaryResults = 1 - pIndCostSummaryResults = 1 - pIndMarginalResults = 0 - pIndEconomicResults = 0 - - # output parameters, variables, and duals to CSV files - if pIndDumpRawResults: - OutputResultsParVarCon(DirName, CaseName, mTEPES, mTEPES) - - if pIndInvestmentResults == 1: - InvestmentResults (DirName, CaseName, mTEPES, mTEPES, pIndTechnologyOutput, pIndPlotOutput) - if pIndGenerationOperationResults == 1: - GenerationOperationResults (DirName, CaseName, mTEPES, mTEPES, pIndTechnologyOutput, pIndAreaOutput, pIndPlotOutput) - if mTEPES.ch and mTEPES.pIndHeat == 1: - GenerationOperationHeatResults(DirName, CaseName, mTEPES, mTEPES, pIndTechnologyOutput, pIndAreaOutput, pIndPlotOutput) - if pIndESSOperationResults == 1 and mTEPES.es: - ESSOperationResults (DirName, CaseName, mTEPES, mTEPES, pIndTechnologyOutput, pIndAreaOutput, pIndPlotOutput) - if pIndReservoirOperationResults == 1 and mTEPES.rs and mTEPES.pIndHydroTopology == 1: - ReservoirOperationResults (DirName, CaseName, mTEPES, mTEPES, pIndTechnologyOutput, pIndPlotOutput) - if pIndNetworkH2OperationResults == 1 and mTEPES.pa and mTEPES.pIndHydrogen == 1: - NetworkH2OperationResults (DirName, CaseName, mTEPES, mTEPES) - if pIndNetworkHeatOperationResults == 1 and mTEPES.ha and mTEPES.pIndHeat == 1: - NetworkHeatOperationResults (DirName, CaseName, mTEPES, mTEPES) - if pIndFlexibilityResults == 1: - FlexibilityResults (DirName, CaseName, mTEPES, mTEPES) - if pIndReliabilityResults == 1: - ReliabilityResults (DirName, CaseName, mTEPES, mTEPES) - if pIndNetworkOperationResults == 1: - NetworkOperationResults (DirName, CaseName, mTEPES, mTEPES) - if pIndNetworkMapResults == 1: - NetworkMapResults (DirName, CaseName, mTEPES, mTEPES) - if pIndOperationSummaryResults == 1: - OperationSummaryResults (DirName, CaseName, mTEPES, mTEPES) - if pIndCostSummaryResults == 1: - CostSummaryResults (DirName, CaseName, mTEPES, mTEPES) - if pIndMarginalResults == 1: - MarginalResults (DirName, CaseName, mTEPES, mTEPES, pIndPlotOutput) - if pIndEconomicResults == 1: - EconomicResults (DirName, CaseName, mTEPES, mTEPES, pIndAreaOutput, pIndPlotOutput) + start = time.time() + path = os.path.join(DirName, CaseName) + + idx = _initialize_indices() + pOut, pLog = idx[pIndOutputResults], idx[pIndLogConsole] + + model_name = ( + 'Open Generation, Storage, and Transmission Operation and Expansion Planning Model ' + 'with RES and ESS (openTEPES) - Version 4.18.6 - July 22, 2025' + ) + mTEPES = _build_model(model_name) + with open(os.path.join(path, f'openTEPES_version_{CaseName}.log'), 'w') as logf: + print(model_name, file=logf) + + _reading_data(DirName, CaseName, mTEPES, pLog) + _configure_basic_components(DirName, CaseName, mTEPES, pLog) + TotalObjectiveFunction(mTEPES, mTEPES, pLog) + InvestmentModelFormulation(mTEPES, mTEPES, pLog) + + _assert_single_objective_active(mTEPES) + _configure_formulation_flags(mTEPES) + mTEPES.pWriteLP = False + + try: + process_stage_loop(DirName, CaseName, SolverName, mTEPES, pLog) + # show results or not + init_flags_value = 0 + flags = _configure_output_flags(pOut, init_flags_value) + # ensure marginal always on in else-case + flags['pIndMarginalResults'] = 1 + finalize_and_output(DirName, CaseName, mTEPES, flags) + finally: + for c in mTEPES.component_objects(pyo.Constraint): + c.activate() + for p, sc in mTEPES.ps: + mTEPES.pScenProb[p, sc] = 1.0 return mTEPES