2
$\begingroup$

I am working on a Pyomo optimization model that aims to select the cheapest fuel for production at different facilities over several time periods. The model should select Mode1 for Tech1 when Fuel1 is cheaper and Mode2 when Fuel2 is cheaper. However, the current implementation consistently selects Mode2 regardless of the fuel prices.

I'm new to using Big M technique to linearise non-linearities so I expect the problem lies in there somewhere.

Model Setup:

Time Periods: 2024, 2025, 2026 
Facilities: Facility1, Facility2 
Technologies: Tech1 
Fuels: Fuel1, Fuel2 
Operating Modes: Mode1 (Fuel1), Mode2 (Fuel2) 
Fuel Prices: 2024: Fuel1 = 50, Fuel2 = 70 2025: Fuel1 = 50, Fuel2 = 70 2026: Fuel1 = 70, Fuel2 = 50 

Desired Behavior:

2024, 2025: Select Mode1 for Tech1 (Fuel1 is cheaper) 
2026: Select Mode2 for Tech1 (Fuel2 is cheaper)

Issue: The model incorrectly selects Mode2 in all years. Here is the output showing the operating modes selected and the auxiliary variable values:

OperatingModeActive:
(2024, 'Facility1', 'Tech1', 'Mode1') 0.00
(2024, 'Facility1', 'Tech1', 'Mode2') 1.00
(2024, 'Facility2', 'Tech1', 'Mode1') 0.00
(2024, 'Facility2', 'Tech1', 'Mode2') 1.00
(2025, 'Facility1', 'Tech1', 'Mode1') 0.00
(2025, 'Facility1', 'Tech1', 'Mode2') 1.00
(2025, 'Facility2', 'Tech1', 'Mode1') 0.00
(2025, 'Facility2', 'Tech1', 'Mode2') 1.00
(2026, 'Facility1', 'Tech1', 'Mode1') 0.00
(2026, 'Facility1', 'Tech1', 'Mode2') 1.00
(2026, 'Facility2', 'Tech1', 'Mode1') 0.00
(2026, 'Facility2', 'Tech1', 'Mode2') 1.00

AuxVar:
(2024, 'Facility1', 'Tech1', 'Mode1') 5.00
(2024, 'Facility1', 'Tech1', 'Mode2') 0.00
(2024, 'Facility2', 'Tech1', 'Mode1') 5.00
(2024, 'Facility2', 'Tech1', 'Mode2') 0.00
(2025, 'Facility1', 'Tech1', 'Mode1') 5.00
(2025, 'Facility1', 'Tech1', 'Mode2') 0.00
(2025, 'Facility2', 'Tech1', 'Mode1') 5.00
(2025, 'Facility2', 'Tech1', 'Mode2') 0.00
(2026, 'Facility1', 'Tech1', 'Mode1') 0.00
(2026, 'Facility1', 'Tech1', 'Mode2') 5.00
(2026, 'Facility2', 'Tech1', 'Mode1') 0.00
(2026, 'Facility2', 'Tech1', 'Mode2') 5.00

Here is the code (I've simplified a lot).

##main.py
from pathlib import Path
from pyomo import environ as pyo
from pyomo.core import value

SOLVER = {
    "name": "glpk",
    "path": Path.home() / "glpk-4.65/w64/glpsol",
}

BIG_M = 1e6  # Big-M constant for linearization
HOURS_PER_YEAR = 8760

def get_fuel_prices_data():
    fuel_prices_data = {
        (2024, 'Country1', 'Fuel1'): 50,
        (2025, 'Country1', 'Fuel1'): 50,
        (2026, 'Country1', 'Fuel1'): 70,
        (2024, 'Country1', 'Fuel2'): 70,
        (2025, 'Country1', 'Fuel2'): 70,
        (2026, 'Country1', 'Fuel2'): 50,
    }
    return fuel_prices_data

def dump_variable(variable, variable_name="Variable"):
    print(f"\n{variable_name}:")
    for i in variable:
        val = value(variable[i])
        val = f"{i} {val:,.2f}"
        print(val)

def get_facility_country(facility):
    return "Country1"

def initialise_constraints(model):
    def facility_production_rule(model, y, f):
        production_per_operational_unit = 1
        production = sum(model.AuxVar[y, f, t, m] * production_per_operational_unit for t in model.Technologies for m in model.OperatingModes)
        return production

    model.FacilityProduction = pyo.Expression(model.TimePeriods, model.Facilities, rule=facility_production_rule)

    def facility_production_target_constraint_rule(model, y, f):
        facility_production_target = model.ProductionTarget[y, f]
        facility_production = model.FacilityProduction[y, f]
        return facility_production >= facility_production_target

    model.ProductionTargetConstraint = pyo.Constraint(
        model.TimePeriods, model.Facilities, rule=facility_production_target_constraint_rule
    )

    def mode_selection_constraint_rule(model, y, f, t):
        """
        Ensure that exactly one operating mode is active for each combination of time period, facility,
        and technology.

        Args:
            model: The Pyomo model instance.
            y: The time period.
            f: The facility.
            t: The technology.

        Returns:
            A constraint ensuring only one operating mode is active.
        """
        return sum(model.OperatingModeActive[y, f, t, m] for m in model.OperatingModes) == 1

    model.OperatingModeSelectionConstraint = pyo.Constraint(
        model.TimePeriods, model.Facilities, model.Technologies, rule=mode_selection_constraint_rule
    )

    def aux_var_upper_bound_rule1(model, y, f, t, m):
        """
        Ensure that the auxiliary variable (AuxVar) is less than or equal to the operational units when the
        operating mode is active.

        Args:
            model: The Pyomo model instance.
            y: The time period.
            f: The facility.
            t: The technology.
            m: The operating mode.

        Returns:
            A constraint enforcing the upper bound of AuxVar.
        """
        return model.AuxVar[y, f, t, m] <= model.OperationalUnits[y, f, t]

    model.AuxVarUpperBound1 = pyo.Constraint(
        model.TimePeriods, model.Facilities, model.Technologies, model.OperatingModes,
        rule=aux_var_upper_bound_rule1
    )

    def aux_var_upper_bound_rule2(model, y, f, t, m):
        """
        Ensure that the auxiliary variable (AuxVar) is less than or equal to a large constant (BIG_M) multiplied by
        the binary variable indicating if the operating mode is active.

        Args:
            model: The Pyomo model instance.
            y: The time period.
            f: The facility.
            t: The technology.
            m: The operating mode.

        Returns:
            A constraint enforcing the second upper bound of AuxVar.
        """
        return model.AuxVar[y, f, t, m] <= BIG_M * model.OperatingModeActive[y, f, t, m]

    model.AuxVarUpperBound2 = pyo.Constraint(
        model.TimePeriods, model.Facilities, model.Technologies, model.OperatingModes,
        rule=aux_var_upper_bound_rule2
    )

    def aux_var_lower_bound_rule(model, y, f, t, m):
        """
        Ensure that the auxiliary variable (AuxVar) is greater than or equal to the operational units minus a large
        constant (BIG_M) multiplied by one minus the binary variable indicating if the operating mode is active.

        Args:
            model: The Pyomo model instance.
            y: The time period.
            f: The facility.
            t: The technology.
            m: The operating mode.

        Returns:
            A constraint enforcing the lower bound of AuxVar.
        """
        return model.AuxVar[y, f, t, m] >= model.OperationalUnits[y, f, t] - BIG_M * (
                1 - model.OperatingModeActive[y, f, t, m])

    model.AuxVarLowerBound = pyo.Constraint(
        model.TimePeriods, model.Facilities, model.Technologies, model.OperatingModes,
        rule=aux_var_lower_bound_rule
    )

    return model

def initialise_objective(model):
    model.DiscountRate = pyo.Param(initialize=0.1, domain=pyo.NonNegativeReals)

    def energy_opex_rule(model):
        return sum(
            model.FuelsPrices[y, get_facility_country(f), model.TechnologyFuel[t, m]] *
            model.AuxVar[y, f, t, m]
            for y in model.TimePeriods
            for f in model.Facilities
            for t in model.Technologies
            for m in model.OperatingModes
        )

    model.EnergyOpex = pyo.Expression(rule=energy_opex_rule)

    def npv_energy_opex_rule(model):
        base_year = min(model.TimePeriods)
        return sum(
            model.EnergyOpex / ((1 + model.DiscountRate) ** (y - base_year))
            for y in model.TimePeriods
        )

    model.NPV_EnergyOpex = pyo.Expression(rule=npv_energy_opex_rule)
    model.Objective = pyo.Objective(expr=model.NPV_EnergyOpex, sense=pyo.minimize)

    return model

def create_model():
    model = pyo.ConcreteModel()

    # Define sets
    model.TimePeriods = pyo.Set(initialize=[2024, 2025, 2026])
    model.Facilities = pyo.Set(initialize=['Facility1', 'Facility2'])
    model.Technologies = pyo.Set(initialize=['Tech1'])
    model.Countries = pyo.Set(initialize=['Country1'])
    model.Fuels = pyo.Set(initialize=['Fuel1', 'Fuel2'])
    model.OperatingModes = pyo.Set(initialize(['Mode1', 'Mode2'])

    # Define parameters
    model.FuelsPrices = pyo.Param(model.TimePeriods, model.Countries, model.Fuels, initialize=get_fuel_prices_data(),
                                  within=pyo.NonNegativeReals, mutable=True)
    model.ProductionTarget = pyo.Param(model.TimePeriods, model.Facilities, initialize=5)

    technology_fuel_data = {
        ('Tech1', 'Mode1'): 'Fuel1',
        ('Tech1', 'Mode2'): 'Fuel2',
    }
    model.TechnologyFuel = pyo.Param(model.Technologies, model.OperatingModes, initialize=technology_fuel_data)

    # Define variables
    model.OperatingModeActive = pyo.Var(model.TimePeriods, model.Facilities, model.Technologies, model.OperatingModes,
                                        domain=pyo.Binary)
    model.OperationalUnits = pyo.Var(model.TimePeriods, model.Facilities, model.Technologies,
                                     domain=pyo.NonNegativeIntegers)

    # Auxiliary variable for linearization
    model.AuxVar = pyo.Var(model.TimePeriods, model.Facilities, model.Technologies, model.OperatingModes,
                           domain=pyo.NonNegativeReals)

    return model

Test code:

##test.py
import pytest
from pyomo import environ as pyo
from main import SOLVER, create_model, initialise_constraints, initialise_objective, dump_variable

def mock_get_fuel_prices_data():
    return {
        (2024, 'Country1', 'Fuel1'): 50,
        (2025, 'Country1', 'Fuel1'): 50,
        (2026, 'Country1', 'Fuel1'): 70,
        (2024, 'Country1', 'Fuel2'): 70,
        (2025, 'Country1', 'Fuel2'): 70,
        (2026, 'Country1', 'Fuel2'): 50,
    }

def test_cheapest_fuel_selected(monkeypatch):
    monkeypatch.setattr('cnz.optimization.scratch.main.get_fuel_prices_data', mock_get_fuel_prices_data)

    # Create the model
    model = create_model()
    model = initialise_constraints(model)
    model = initialise_objective(model)

    # Print fuel prices to confirm they are set correctly
    print("Fuel Prices:")
    dump_variable(model.FuelsPrices, "FuelsPrices")

    # Print technology-fuel mapping
    print("\nTechnology to Fuel Mapping:")
    for t in model.Technologies:
        for m in model.OperatingModes:
            print(f"Technology {t}, Mode {m} uses Fuel: {model.TechnologyFuel[t, m]}")

    # Solve the model
    solver = pyo.SolverFactory(SOLVER["name"], executable=SOLVER["path"])
    results = solver.solve(model, tee=True)

    # Check if the model solved successfully
    assert results.solver.termination_condition == pyo.TerminationCondition.optimal

    # Print the operating mode active variable to see which mode is selected
    dump_variable(model.OperatingModeActive, "OperatingModeActive")
    dump_variable(model.OperationalUnits, "OperationalUnits")
    dump_variable(model.AuxVar, "AuxVar")

    # Check if the cheapest fuel mode is selected
    for y in model.TimePeriods:
        for f in model.Facilities:
            for t in model.Technologies:
                mode1_value = model.OperatingModeActive[y, f, t, 'Mode1'].value
                mode2_value = model.OperatingModeActive[y, f, t, 'Mode2'].value
                print(f"Checking {y}, {f}, {t}: Mode1 - {mode1_value}, Mode2 - {mode2_value}")
                if y in [2024, 2025]:
                    assert mode1_value == 1.0, f"Failed at {y}, {f}, {t}: Expected Mode1 to be 1.0 but was {mode1_value}"
                elif y == 2026:
                    assert mode2_value == 1.0, f"Failed at {y}, {f}, {t}: Expected Mode2 to be 1.0 but was {mode2_value}"

if __name__ == "__main__":
    pytest.main()
$\endgroup$
2
  • 6
    $\begingroup$ Big M of 1e6 is huge. Even if the model is "correct" (I haven't checked), because of integrality tolerance, it's possible there is trickle flow; i.e., the intended logic is not enforced by the constraints. Ideally, Big M should be just large enough, not a huge number; that is to prevent trickle flow, as well as result in tighter relaxations, and therefore faster solve times. See or.stackexchange.com/search?q=trickle+flow $\endgroup$ Commented May 22 at 21:33
  • 1
    $\begingroup$ Wow, thank you! I set it to 1e3 and it works right away $\endgroup$
    – Jamie Bull
    Commented May 22 at 21:34

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.