Exploring the Pyomo Model#

This guide explains the internal structure of the SDOM Pyomo model for advanced users and developers.

Model Architecture#

SDOM uses pyomo to write the optimization model.

When you run

```python
model = initialize_model(
    data,
    n_hours=n_steps,
    with_resilience_constraints=with_resilience_constraints
)
```

the function initialize_model() returns the pyomo instance of the optimization model, which in this case is stores in the variable model.

SDOM leverages on pyomo blocks to separate in different blocks the variables, parameters, expressions, constraints, etc of diferent model components. In this way, in that pyomo instance, SDOM creates the following blocks:

  • Core optimization Blocks (Blocks that include variables, sets and parameters)

    • thermal

    • pv

    • wind

    • hydro

    • storage

  • Blocks containing only parameters (Do not include any decision variables)

    • demand

    • nuclear

    • other_renewables

  • Optional blocks (These are created depending on the configuration provided by the user)

    • imports

    • exports

    • resiliency

Using pprint() function, you’ll be able to explore model objects in a simple way. If you run:

from sdom import load_data, initialize_model

data = load_data('./Data/scenario/')
model = initialize_model(data, n_hours=8760)

# Access model blocks
model.pv.pprint()          # Solar PV block
model.wind.pprint()         # Wind block
model.storage.pprint()      # Storage block
model.thermal.pprint()      # Thermal block
model.hydro..pprint()        # Hydropower block
model.nuclear.pprint()      # Nuclear block (parameters only)
model.demand.pprint()       # Demand block (parameters only)

for example, if you run:

```python
model.thermal.heat_rate.pprint()
```

You can obtain:

```python
heat_rate : Size=13, Index=thermal.plants_set, Domain=Any, Default=None, Mutable=False
Key  : Value
106c :   1.0
147c :   1.0
162c :   1.0
163c :   1.0
167c :   1.0
170c :   1.0
221c :   1.0
223c :   1.0
224c :   1.0
235c :   1.0
241c :   1.0
 83c :   1.0
 98c :   1.0
```

Model Components#

Sets#

Key model sets:

# Hourly time set (1-based indexing)
print(len(model.h))  # e.g., 8760

# VRE plant sets
print(list(model.pv.plants_set))    # ['101', '202', '303', ...]
print(list(model.wind.plants_set))  # ['501', '602', ...]

# Storage technology sets
print(list(model.storage.j))  # All storage techs: ['Li-Ion', 'CAES', 'PHS', 'H2']
print(list(model.storage.b))  # Coupled techs: ['Li-Ion', 'PHS']

# Thermal unit set
print(list(model.thermal.plants_set))  # ['GasCC_1', ...]

# Budget set (if using budget formulation)
if hasattr(model.hydro, 'budget_set'):
    print(list(model.hydro.budget_set))  # [1, 2, ..., 12] for monthly

Parameters#

Access model parameters:

# System-level parameters
print(f"Discount rate: {model.r.value}")
print(f"Carbon-free target: {model.GenMix_Target.value}")

# Time-series parameters
print(f"Hour 1 load: {model.demand.ts_parameter[1]} MW")
print(f"Hour 1 nuclear: {model.nuclear.ts_parameter[1]} MW")

# VRE parameters
solar_plant = list(model.pv.plants_set)[0]
print(f"Solar plant {solar_plant} CAPEX: {model.pv.CAPEX_M[solar_plant]} $/kW")
print(f"Solar plant {solar_plant} capacity: {model.pv.plant_max_capacity[solar_plant]} MW")

# Storage parameters
print(f"Li-Ion efficiency: {model.storage.data['Eff', 'Li-Ion']}")
print(f"Li-Ion E_Capex: {model.storage.data['E_Capex', 'Li-Ion']} $/kWh")

Variables#

Model decision variables:

from sdom.common.utilities import safe_pyomo_value

# VRE installed capacity (continuous variables)
for plant in model.pv.plants_set:
    capacity = safe_pyomo_value(model.pv.plant_installed_capacity[plant])
    print(f"Solar plant {plant}: {capacity:.2f} MW")

# VRE generation (continuous, hourly)
gen_h1 = safe_pyomo_value(model.pv.generation[1])
print(f"Solar generation hour 1: {gen_h1:.2f} MW")

# Storage variables
for tech in model.storage.j:
    pcha = safe_pyomo_value(model.storage.Pcha[tech])
    pdis = safe_pyomo_value(model.storage.Pdis[tech])
    ecap = safe_pyomo_value(model.storage.Ecap[tech])
    print(f"{tech}: Pcha={pcha:.2f} MW, Pdis={pdis:.2f} MW, Ecap={ecap:.2f} MWh")

# Storage state of charge (hourly)
soc_h100 = safe_pyomo_value(model.storage.SOC[100, 'Li-Ion'])
print(f"Li-Ion SOC at hour 100: {soc_h100:.2f} MWh")

# Thermal capacity and generation
thermal_cap = safe_pyomo_value(model.thermal.total_installed_capacity)
thermal_gen_h1 = safe_pyomo_value(model.thermal.generation[1, 'GasCC_1'])
print(f"Thermal capacity: {thermal_cap:.2f} MW")
print(f"Thermal gen hour 1: {thermal_gen_h1:.2f} MW")

Expressions#

Computed expressions (not variables):

# Total installed capacity expressions
total_solar = safe_pyomo_value(model.pv.total_installed_capacity)
total_wind = safe_pyomo_value(model.wind.total_installed_capacity)
print(f"Total solar: {total_solar:.2f} MW")
print(f"Total wind: {total_wind:.2f} MW")

# Cost expressions
solar_capex = safe_pyomo_value(model.pv.capex_cost_expr)
solar_fom = safe_pyomo_value(model.pv.fixed_om_cost_expr)
print(f"Solar CAPEX: ${solar_capex:,.2f}")
print(f"Solar FOM: ${solar_fom:,.2f}")

# Total generation expressions
total_solar_gen = safe_pyomo_value(model.pv.total_generation)
print(f"Total solar generation: {total_solar_gen:.2f} MWh")

Constraints#

Model constraints are organized by block:

# List all constraints
for constraint in model.component_objects(pyo.Constraint, active=True):
    print(f"Constraint: {constraint.name}")
    if constraint.is_indexed():
        print(f"  Indices: {len(constraint)} constraints")
    else:
        print(f"  Single constraint")

Key constraint types:

  • Energy Balance (system_energy_balance_rule): Supply = Demand every hour

  • Carbon Target (GenMix_constraint): Clean energy ≥ target × total generation

  • VRE Balance: Generation ≤ CF × Capacity for each hour/plant

  • Storage Dynamics: State of charge evolution constraints

  • Storage Bounds: SOC within [0, Energy Capacity]

  • Thermal Bounds: Generation within [Min, Max Capacity]

  • Hydro Budget: Energy generation sums match budget (if applicable)

Inspecting the Model#

Examine Specific Constraints#

import pyomo.environ as pyo

# Energy balance constraint (one per hour)
balance_h1 = model.system_energy_balance[1]
print(f"Balance constraint hour 1: {balance_h1.expr}")

# GenMix constraint
genmix = model.GenMix_constraint
print(f"GenMix constraint: {genmix.expr}")

# Storage SOC dynamics (example for hour 1, Li-Ion)
if (1, 'Li-Ion') in model.storage.storage_charge_discharge_dynamics:
    dynamics = model.storage.storage_charge_discharge_dynamics[1, 'Li-Ion']
    print(f"Storage dynamics constraint: {dynamics.expr}")

Display Variable Bounds#

# Check variable bounds
for tech in model.storage.j:
    var = model.storage.Pcha[tech]
    print(f"{tech} Pcha bounds: [{var.lb}, {var.ub}]")
    
for h in list(model.h)[:5]:  # First 5 hours
    var = model.pv.generation[h]
    print(f"Solar gen hour {h} bounds: [{var.lb}, {var.ub}]")

Model Modifications#

Changing Parameters Between Runs#

# Load data and initialize model once
data = load_data('./Data/scenario/')
model = initialize_model(data, n_hours=8760)

# Change GenMix_Target (mutable parameter)
model.GenMix_Target = 0.90  # 90% clean energy
results_90 = run_solver(model, solver_config)

model.GenMix_Target = 0.95  # 95% clean energy
results_95 = run_solver(model, solver_config)

model.GenMix_Target = 0.99  # 99% clean energy
results_99 = run_solver(model, solver_config)

Adding Custom Constraints#

import pyomo.environ as pyo

# Add minimum wind capacity constraint
def min_wind_capacity_rule(model):
    return model.wind.total_installed_capacity >= 10000  # 10 GW minimum

model.min_wind_constraint = pyo.Constraint(rule=min_wind_capacity_rule)

# Add maximum solar/wind ratio
def max_solar_wind_ratio_rule(model):
    return model.pv.total_installed_capacity <= 2 * model.wind.total_installed_capacity

model.max_ratio_constraint = pyo.Constraint(rule=max_solar_wind_ratio_rule)

Debugging Model Issues#

Check for Infeasibility#

from pyomo.util.infeasible import log_infeasible_constraints

# After failed solve
results_list, best_result, solver_result = run_solver(model, solver_config)

if solver_result.solver.termination_condition != pyo.TerminationCondition.optimal:
    print("Model is infeasible or suboptimal")
    log_infeasible_constraints(model)

Validate Model Construction#

# Check that variables are defined
assert hasattr(model.pv, 'generation'), "Solar generation variable missing"
assert hasattr(model.storage, 'SOC'), "Storage SOC variable missing"

# Check that sets are populated
assert len(model.pv.plants_set) > 0, "No solar plants loaded"
assert len(model.storage.j) > 0, "No storage technologies"

# Check parameter values
assert model.r.value > 0, "Discount rate must be positive"
assert 0 <= model.GenMix_Target.value <= 1, "GenMix_Target out of range"

Advanced: Direct Pyomo Operations#

import pyomo.environ as pyo

# Write model to file
model.write('sdom_model.lp', io_options={'symbolic_solver_labels': True})

# Get solver results details
print(f"Solver status: {solver_result.solver.status}")
print(f"Termination condition: {solver_result.solver.termination_condition}")
print(f"Solve time: {solver_result.solver.time:.2f} seconds")

# Access dual values (shadow prices) if solver supports
if hasattr(solver_result, 'problem'):
    print(f"Objective value: {solver_result.problem.lower_bound}")
    print(f"Best bound: {solver_result.problem.upper_bound}")

Model Formulation Files#

The model formulation is split across multiple files in src/sdom/models/:

  • formulations_system.py: Objective function and system constraints

  • formulations_vre.py: Solar and wind variables/constraints

  • formulations_storage.py: Storage variables/constraints

  • formulations_thermal.py: Thermal unit variables/constraints

  • formulations_hydro.py: Hydropower constraints

  • formulations_imports_exports.py: Cross-border trade

  • models_utils.py: Helper functions for expressions and constraints

Next Steps#