import pyomo.environ as pyo
[docs]class LinearProgrammingModel:
"""Solve unit commitment and economic dispatch problem using linear programming.
:param data: Full path with file name and extention to the input .dat object
:type data: str
"""
def __init__(self, data):
self.data = data
# instantiate model
self.model = pyo.AbstractModel()
# generators by fuel-type
self.model.Coal = pyo.Set()
self.model.Oil = pyo.Set()
self.model.Gas = pyo.Set()
self.model.Hydro = pyo.Set()
self.model.Solar = pyo.Set()
self.model.Wind = pyo.Set()
# all generators
self.model.Thermal = self.model.Coal | self.model.Oil | self.model.Gas
self.model.Dispatchable = self.model.Thermal | self.model.Hydro
self.model.Generators = self.model.Dispatchable | self.model.Solar | self.model.Wind
# transmission sets
self.model.lines = pyo.Set()
self.model.buses = pyo.Set()
# Generator type
self.model.typ = pyo.Param(self.model.Generators, within=pyo.Any)
# Node name
self.model.node = pyo.Param(self.model.Generators, within=pyo.Any)
# Max capacity
self.model.maxcap = pyo.Param(self.model.Generators)
# Min capacity
self.model.mincap = pyo.Param(self.model.Generators)
# Heat rate
self.model.heat_rate = pyo.Param(self.model.Generators)
# Variable O&M
self.model.var_om = pyo.Param(self.model.Generators)
# Fixed O&M cost
self.model.no_load = pyo.Param(self.model.Generators)
# Start cost
self.model.st_cost = pyo.Param(self.model.Generators)
# Ramp rate
self.model.ramp = pyo.Param(self.model.Generators)
# Minimun up time
self.model.minup = pyo.Param(self.model.Generators)
# Minimun down time
self.model.mindn = pyo.Param(self.model.Generators)
self.model.Reactance = pyo.Param(self.model.lines)
self.model.FlowLim = pyo.Param(self.model.lines)
self.model.LinetoBusMap = pyo.Param(self.model.lines, self.model.buses)
self.model.BustoUnitMap = pyo.Param(self.model.Generators, self.model.buses)
######===== Parameters/initial_conditions to run simulation ======#######
# full range of time series information
self.model.SimHours = pyo.Param(within=pyo.PositiveIntegers)
self.model.SH_periods = pyo.RangeSet(1, self.model.SimHours + 1)
self.model.SimDays = pyo.Param(within=pyo.PositiveIntegers)
self.model.SD_periods = pyo.RangeSet(1, self.model.SimDays + 1)
# Operating horizon information
self.model.HorizonHours = pyo.Param(within=pyo.PositiveIntegers)
self.model.HH_periods = pyo.RangeSet(0, self.model.HorizonHours)
self.model.hh_periods = pyo.RangeSet(1, self.model.HorizonHours)
self.model.ramp_periods = pyo.RangeSet(2, 24)
# Demand over simulation period
self.model.SimDemand = pyo.Param(self.model.buses * self.model.SH_periods, within=pyo.NonNegativeReals)
# Horizon demand
self.model.HorizonDemand = pyo.Param(self.model.buses * self.model.hh_periods, within=pyo.NonNegativeReals,
mutable=True)
# variable resources over simulation period
self.model.SimHydro = pyo.Param(self.model.Hydro, self.model.SH_periods, within=pyo.NonNegativeReals)
self.model.SimSolar = pyo.Param(self.model.Solar, self.model.SH_periods, within=pyo.NonNegativeReals)
self.model.SimWind = pyo.Param(self.model.Wind, self.model.SH_periods, within=pyo.NonNegativeReals)
# Variable resources over horizon
self.model.HorizonHydro = pyo.Param(self.model.Hydro, within=pyo.NonNegativeReals, mutable=True)
self.model.HorizonSolar = pyo.Param(self.model.Solar, self.model.hh_periods, within=pyo.NonNegativeReals,
mutable=True)
self.model.HorizonWind = pyo.Param(self.model.Wind, self.model.hh_periods, within=pyo.NonNegativeReals,
mutable=True)
# Must run resources
self.model.Must = pyo.Param(self.model.buses, within=pyo.NonNegativeReals)
######=======================Decision variables======================########
# amount of day-ahead energy generated by each generator at each hour
self.model.mwh = pyo.Var(self.model.Generators, self.model.HH_periods, within=pyo.NonNegativeReals,
initialize=0)
# 1 if unit is on in hour i, otherwise 0
self.model.on = pyo.Var(self.model.Dispatchable, self.model.HH_periods, within=pyo.NonNegativeReals,
initialize=0)
# 1 if unit is switching on in hour i, otherwise 0
self.model.switch = pyo.Var(self.model.Dispatchable, self.model.HH_periods, within=pyo.NonNegativeReals,
initialize=0)
# slack variables
self.model.S = pyo.Var(self.model.buses, self.model.hh_periods, within=pyo.NonNegativeReals, initialize=0)
# transmission line variables
self.model.Flow = pyo.Var(self.model.lines, self.model.hh_periods)
self.model.Theta = pyo.Var(self.model.buses, self.model.hh_periods)
######================Objective function=============########
# minimization of system costs
self.model.SystemCost = pyo.Objective(rule=self.system_cost_obj_function, sense=pyo.minimize)
#####========== Logical Constraint =========#############
#####==========Ramp Rate Constraints =========#############
self.model.RampCon1 = pyo.Constraint(self.model.Thermal, self.model.ramp_periods, rule=self.ramp_constraint_one)
self.model.RampCon2 = pyo.Constraint(self.model.Thermal, self.model.ramp_periods, rule=self.ramp_constraint_two)
# Constraints for Max & Min Capacity of dispatchable resources
self.model.MaxCap = pyo.Constraint(self.model.Dispatchable, self.model.hh_periods,
rule=self.maximum_capacity_constraint)
self.model.MinCap = pyo.Constraint(self.model.Dispatchable, self.model.hh_periods,
rule=self.minimum_capacity_constraint)
# Max capacity constraints on domestic hydropower
self.model.HydroConstraint = pyo.Constraint(self.model.Hydro, self.model.hh_periods,
rule=self.hydropower_max_capacity_constraint)
# Max capacity constraints on solar
self.model.SolarConstraint = pyo.Constraint(self.model.Solar, self.model.hh_periods,
rule=self.solar_max_capacity_constraint)
# Max capacity constraints on wind
self.model.WindConstraint = pyo.Constraint(self.model.Wind, self.model.hh_periods,
rule=self.wind_max_capacity_constraint)
self.model.Node_Constraint = pyo.Constraint(self.model.buses, self.model.hh_periods,
rule=self.nodal_balance_constraint)
self.model.FlowL_Constraint = pyo.Constraint(self.model.lines, self.model.hh_periods, rule=self.energy_flow)
self.model.ThetaB_Constraint = pyo.Constraint(self.model.hh_periods, rule=self.theta_bus_voltage_angle)
self.model.FlowU_Constraint = pyo.Constraint(self.model.lines, self.model.hh_periods,
rule=self.energy_flow_max_line_limit)
self.model.FlowLL_Constraint = pyo.Constraint(self.model.lines, self.model.hh_periods,
rule=self.energy_flow_min_line_limit)
[docs] @staticmethod
def system_cost_obj_function(model):
"""Objective function for system costs.
:param model: Pyomo model instance
"""
fixed = sum(model.no_load[j] * model.on[j, i] for i in model.hh_periods for j in model.Dispatchable)
starts = sum(model.st_cost[j] * model.switch[j, i] for i in model.hh_periods for j in model.Dispatchable)
coal = sum(
model.mwh[j, i] * (model.heat_rate[j] * 2 + model.var_om[j]) for i in model.hh_periods for j in model.Coal)
oil = sum(
model.mwh[j, i] * (model.heat_rate[j] * 10 + model.var_om[j]) for i in model.hh_periods for j in model.Oil)
gas = sum(
model.mwh[j, i] * (model.heat_rate[j] * 4.5 + model.var_om[j]) for i in model.hh_periods for j in model.Gas)
slack = sum(model.S[z, i] * 100000 for i in model.hh_periods for z in model.buses)
return coal + oil + gas + slack + fixed + starts
[docs] @staticmethod
def ramp_constraint_one(model, j, i):
"""Ramp constraint for the thermal generating units.
:param model: Pyomo model instance
:param j: TODO
:parma i: TODO
"""
a = model.mwh[j, i]
b = model.mwh[j, i - 1]
return a - b <= model.ramp[j]
[docs] @staticmethod
def ramp_constraint_two(model, j, i):
"""Ramp constraint for the thermal generating units.
:param model: Pyomo model instance
:param j: TODO
:parma i: TODO
"""
a = model.mwh[j, i]
b = model.mwh[j, i - 1]
return b - a <= model.ramp[j]
[docs] @staticmethod
def maximum_capacity_constraint(model, j, i):
"""Maximum capacity constraint of dispatchable resources"""
return model.mwh[j, i] <= model.on[j, i] * model.maxcap[j]
[docs] @staticmethod
def minimum_capacity_constraint(model, j, i):
return model.mwh[j, i] >= model.on[j, i] * model.mincap[j]
[docs] @staticmethod
def hydropower_max_capacity_constraint(model, j, i):
"""Max capacity constraints on domestic hydropower."""
daily = sum(model.mwh[j, i] for i in model.hh_periods)
return daily <= model.HorizonHydro[j]
[docs] @staticmethod
def solar_max_capacity_constraint(model, j, i):
"""Max capacity constraints on solar"""
return model.mwh[j, i] <= model.HorizonSolar[j, i]
[docs] @staticmethod
def wind_max_capacity_constraint(model, j, i):
"""Max capacity constraints on wind"""
return model.mwh[j, i] <= model.HorizonWind[j, i]
[docs] @staticmethod
def nodal_balance_constraint(model, z, i):
"""Power balance constraint. Creates a balance for each node (supply == demand)"""
power_flow = sum(model.Flow[l, i] * model.LinetoBusMap[l, z] for l in model.lines)
gen = sum(model.mwh[j, i] * model.BustoUnitMap[j, z] for j in model.Generators)
slack = model.S[z, i]
must_run = model.Must[z]
return gen + slack + must_run - power_flow == model.HorizonDemand[z, i]
[docs] @staticmethod
def energy_flow(model, l, i, base_apparent_power=100):
"""Energy flowing through transmission lines in Mw/h.
:param model: TODO
:param l: TODO
:param i: TODO
:param base_apparent_power: TODO
"""
value = sum(model.Theta[z, i] * model.LinetoBusMap[l, z] for z in model.buses)
return base_apparent_power * value == model.Flow[l, i] * model.Reactance[l]
[docs] @staticmethod
def theta_bus_voltage_angle(model, i, reference_bus='MESA_CAL_20', constant_value=0):
"""Volatge angle for the reference bus.
:param model: TODO
:param i: TODO
:param reference_bus: TODO
:type reference_bus: str
:param constant_value: TODO
:type constant_value: int
"""
return model.Theta[reference_bus, i] == constant_value
[docs] @staticmethod
def energy_flow_max_line_limit(model, l, i):
"""Maximum energy flow limit of the transmission line."""
return model.Flow[l, i] <= model.FlowLim[l]
[docs] @staticmethod
def energy_flow_min_line_limit(model, l, i, negative_factor=-1):
"""Minimum energy flow limit of the transmission line."""
return negative_factor * model.Flow[l, i] <= model.FlowLim[l]
[docs] def create_instance(self):
"""Wrapper to create an instance of the model with data."""
return self.model.create_instance(self.data)