I am currently working on a problem in which I want to optimize several objects over a period of time (specifically: the charging profile of several electric vehicles over a certain period of time).
Does it make a difference whether I create a variable (pyo.Var) for each electric vehicle charging profile and index it with the time steps or whether I create a single variable and index it with both the time steps and the electric vehicles? Or to put it another way: Does this difference go beyond personal taste when scaling up the number of electric vehicles or the number of time steps, is there a better/more practical variant in terms of performance?
I'll try to illustrate the difference with the help of an MWE:
The underlying class is the same. This is a very simple model of an electric vehicle that can charge and drive, thereby changing its state of charge:
import pandas as pd
import pyomo.environ as pyo
class EV:
def __init__(self, idx: int, capacity_kwh: float, starting_soc: float,
time_steps: pd.DatetimeIndex, driving_profile_km: list[float]):
self.idx = idx
self.capacity_kwh = capacity_kwh
self.soc = starting_soc
self.driving_profile_km = pd.Series(driving_profile_km, index=time_steps)
def charge(self, power_kw: float):
duration_h = 1.0
self.soc = self.soc + (power_kw * duration_h) / self.capacity_kwh
def drive(self, distance_km: float):
consumption_kwh_per_km = 0.15
self.soc = self.soc - (distance_km * consumption_kwh_per_km) / self.capacity_kwh
What is also the same in both cases is the start of the model creation and the addition of time steps and electric vehicles:
model = pyo.ConcreteModel()
time_steps = pd.date_range("2023-11-15 10:00", "2023-11-15 13:00", freq="h")
model.time_steps = pyo.Set(initialize=time_steps)
ev1 = EV(1, capacity_kwh=40.0, starting_soc=0.5, time_steps=time_steps,
driving_profile_km=[20.0, 0.0, 0.0, 20.0])
ev2 = EV(2, capacity_kwh=30.0, starting_soc=0.5, time_steps=time_steps,
driving_profile_km=[10.0, 5.0, 0.0, 15.0])
ev3 = EV(3, capacity_kwh=20.0, starting_soc=0.5, time_steps=time_steps,
driving_profile_km=[0.0, 30.0, 25.0, 5.0])
model.electric_vehicles = pyo.Set(initialize=[ev1, ev2, ev3])
Now comes the above-mentioned point of creating the decision variable(s). In the first case, I create one variable per electric vehicle and index it with the time steps:
for ev in model.electric_vehicles:
model.add_component(f"charging_profile_ev{ev.idx}",
pyo.Var(model.time_steps, bounds=[0, 11], initialize=0.0))
In the second case, I create only one variable and index it with both the time steps and the electric vehicles:
model.add_component(f"charging_profiles",
pyo.Var(model.time_steps, model.electric_vehicles,
bounds=[0, 11], initialize=0.0))
The creation of the additionally required expressions, constraints and objectives also changes slightly depending on the case. I will briefly show both cases so that the MWE is complete:
Case 1:
def expression_soc(m, ts, ev):
charging_power_kw = m.component(f"charging_profile_ev{ev.idx}")[ts]
driving_distance_km = ev.driving_profile_km.loc[ts]
ev.charge(charging_power_kw)
ev.drive(driving_distance_km)
return ev.soc
model.add_component(f"soc_expression",
pyo.Expression(model.time_steps, model.electric_vehicles,
rule=expression_soc))
def constraint_soc(m, ts, ev):
soc_value = m.component(f"soc_expression")[ts, ev]
return 0.5, soc_value, 1.0
model.add_component(f"soc_constraint",
pyo.Constraint(model.time_steps, model.electric_vehicles,
rule=constraint_soc))
def objective(m):
value = 0
for ev in m.electric_vehicles:
for ts in m.time_steps:
value = value + m.component(f"charging_profile_ev{ev.idx}")[ts]
return value
model.objective = pyo.Objective(rule=objective, sense=pyo.minimize)
results = pyo.SolverFactory("glpk").solve(model)
model.objective.display()
model.charging_profile_ev1.display()
model.charging_profile_ev2.display()
model.charging_profile_ev3.display()
Case 2:
def expression_soc(m, ts, ev):
charging_power_kw = m.component(f"charging_profiles")[ts, ev]
driving_distance_km = ev.driving_profile_km.loc[ts]
ev.charge(charging_power_kw)
ev.drive(driving_distance_km)
return ev.soc
model.add_component(f"soc_expression",
pyo.Expression(model.time_steps, model.electric_vehicles,
rule=expression_soc))
def constraint_soc(m, ts, ev):
soc_value = m.component(f"soc_expression")[ts, ev]
return 0.5, soc_value, 1.0
model.add_component(f"soc_constraint",
pyo.Constraint(model.time_steps, model.electric_vehicles,
rule=constraint_soc))
def objective(m):
value = 0
for ev in m.electric_vehicles:
for ts in m.time_steps:
value = value + m.component(f"charging_profiles")[ts, ev]
return value
model.objective = pyo.Objective(rule=objective, sense=pyo.minimize)
results = pyo.SolverFactory("glpk").solve(model)
model.objective.display()
model.charging_profiles.display()
The corresponding outputs logically look different and you have to access the profiles differently in order to be able to use them in other applications. But are there any other differences or reasons why one or the other approach should be preferred?
100% you want to go with the second route and index everything rather than having a bunch of singleton variables rolling around.
In the case of your model, you have 2 natural sets: time and vehicle, and I would expect most variables to be dual-indexed by those 2 sets.
First off, however (and this may be a long post to show this....sorry!) you are playing with fire in my opinion in how you are treating soc in your model build. You are implicitly relying on the instance variable soc to cover all the time periods, but it is just a single variable. It cannot simultaneously represent the soc in all of your model's periods. You are (kinda) OK in constructing the expression you are doing, but it is fragile. Meaning: It is making correct expressions for the SOC for the model's use, but it is corrupting the instance variable as you are doing it. An example is shown below. Note that the expression for the SOC is what you would expect as it is built sequentially as long as you pass the time steps in the proper order but when I use the ev.soc variable again, it is corrupted in the "bad" constraint because it is still holding the last usage and it doesn't know anything about time. Make sense? You can (if you are careful) just keep using the expression you generated (see the next constraint) because it has captured the correct expressions when it was built. You'd have to use that expression again later to figure out the SOC for similar reasons.
"""
EV charge model
"""
import pyomo.environ as pyo
class SimpleEV:
idx = 1 # rolling counter
def __init__(self):
self.soc = 10 # initial charge
self.idx = SimpleEV.idx
SimpleEV.idx += 1
def add_juice(self, charge):
self.soc += charge
def use_juice(self, charge):
self.soc -= charge
# for cleaner printing in pyomo, provide a __repr__
def __repr__(self):
return f'ev{self.idx}'
m = pyo.ConcreteModel('ev')
m.T = pyo.Set(initialize=range(3)) # time
m.E = pyo.Set(initialize=[SimpleEV()]) # 1 simple EV
m.charge = pyo.Var(m.E, m.T, domain=pyo.NonNegativeReals, doc='charge added to vehicle e at time t')
m.draws = pyo.Param(m.T, initialize={1:5}, default=0)
def expr_generator(m, e: SimpleEV, t):
e.add_juice(m.charge[e, t])
e.use_juice(m.draws[t])
return e.soc
m.soc_expression = pyo.Expression(m.E, m.T, rule=expr_generator)
# this is NOT going to work using the .soc variable:
def bad_min_charge(m, e: SimpleEV, t):
return e.soc >= 2
m.bad_min_charge_constraint = pyo.Constraint(m.E, m.T, rule=bad_min_charge)
def ok_min_charge(m, e: SimpleEV, t):
return m.soc_expression[e, t] >= 2
m.ok_min_charge_constraint = pyo.Constraint(m.E, m.T, rule=ok_min_charge)
m.pprint()
1 Expression Declarations
soc_expression : Size=3, Index=soc_expression_index
Key : Expression
(ev1, 0) : 10 + charge[ev1,0]
(ev1, 1) : 10 + charge[ev1,0] + charge[ev1,1] - 5
(ev1, 2) : 10 + charge[ev1,0] + charge[ev1,1] - 5 + charge[ev1,2]
2 Constraint Declarations
bad_min_charge_constraint : Size=3, Index=bad_min_charge_constraint_index, Active=True
Key : Lower : Body : Upper : Active
(ev1, 0) : 2.0 : 10 + charge[ev1,0] + charge[ev1,1] - 5 + charge[ev1,2] : +Inf : True
(ev1, 1) : 2.0 : 10 + charge[ev1,0] + charge[ev1,1] - 5 + charge[ev1,2] : +Inf : True
(ev1, 2) : 2.0 : 10 + charge[ev1,0] + charge[ev1,1] - 5 + charge[ev1,2] : +Inf : True
ok_min_charge_constraint : Size=3, Index=ok_min_charge_constraint_index, Active=True
Key : Lower : Body : Upper : Active
(ev1, 0) : 2.0 : soc_expression[ev1,0] : +Inf : True
(ev1, 1) : 2.0 : soc_expression[ev1,1] : +Inf : True
(ev1, 2) : 2.0 : soc_expression[ev1,2] : +Inf : True
On to a better idea: I would construct your model similar to below. A few notes:
__hash__ and __eq__ functions.__repr__() in your code so the outputs are more clear. (see the python dox for guidance, I'm cutting a few corners here)"""
improved EV charge model
"""
import pyomo.environ as pyo
class DriveProfile:
def __init__(self, distances):
self.distances = distances # {time:km} dictionary
def get_distance(self, time):
return self.distances.get(time, None) # NONE if no travels at this time block
class EV:
idx = 1 # rolling counter
def __init__(self, max_soc=100, efficiency=2, route: DriveProfile = None):
self.initial_soc = 10 # initial charge
self.max_soc = max_soc # max charge
self.discharge_rate = efficiency # units/km
self.route = route
self.idx = EV.idx
self.soc = None
EV.idx += 1
def load_soc(self, soc_dict: dict):
self.soc = soc_dict
# for cleaner printing in pyomo, provide a __repr__
def __repr__(self):
return f'ev{self.idx}'
# make some fake data...
times = list(range(4))
r1 = DriveProfile({1: 10, 3: 2.4})
r2 = DriveProfile({2: 20})
e1 = EV(route=r1, efficiency=3)
e2 = EV(max_soc=200, route=r2)
# make the model
m = pyo.ConcreteModel('ev')
# SETS
m.T = pyo.Set(initialize=times) # time
m.E = pyo.Set(initialize=[e1, e2]) # ev's
# VARS
m.soc = pyo.Var(m.E, m.T, domain=pyo.NonNegativeReals, doc='state of charge of E at time T')
m.juice = pyo.Var(m.E, m.T, domain=pyo.NonNegativeReals, doc='juice added to E at time T')
# PARAMS
# --- all params are internal to the EV objects. We could capture them here, but not needed. ---
# OBJ: Minimize the juice delivered
m.obj = pyo.Objective(expr=sum(m.juice[e, t] for e in m.E for t in m.T))
# CONSTRAINTS
# max charge
def max_charge(m, e: EV, t):
return m.soc[e, t] <= e.max_soc
m.max_charge = pyo.Constraint(m.E, m.T, rule=max_charge)
def soc_balance(m, e: EV, t):
# find any use (driving)
use = e.route.get_distance(t)
discharge = use * e.discharge_rate if use else 0
if t == m.T.first():
return m.soc[e, t] == e.initial_soc - discharge + m.juice[e, t]
return m.soc[e, t] == m.soc[e, m.T.prev(t)] - discharge + m.juice[e, t]
m.soc_balance = pyo.Constraint(m.E, m.T, rule=soc_balance)
m.pprint()
# SOLVE
solver = pyo.SolverFactory('cbc')
soln = solver.solve(m)
print(soln)
# show the juice profile
m.juice.display()
# pass the soc info back to the ev's
print('The state of charge in the vehicles:')
for ev in m.E:
# we can use items() to get tuples of ( (e, t), soc ) and screen them, then get the .value from the variable...
res = {time: soc.value for (e, time), soc in m.soc.items() if e == ev}
ev.load_soc(res)
print(ev.__repr__(), ev.soc)
soc_balance : Size=8, Index=soc_balance_index, Active=True
Key : Lower : Body : Upper : Active
(ev1, 0) : 0.0 : soc[ev1,0] - (10 + juice[ev1,0]) : 0.0 : True
(ev1, 1) : 0.0 : soc[ev1,1] - (soc[ev1,0] - 30 + juice[ev1,1]) : 0.0 : True
(ev1, 2) : 0.0 : soc[ev1,2] - (soc[ev1,1] + juice[ev1,2]) : 0.0 : True
(ev1, 3) : 0.0 : soc[ev1,3] - (soc[ev1,2] - 7.199999999999999 + juice[ev1,3]) : 0.0 : True
(ev2, 0) : 0.0 : soc[ev2,0] - (10 + juice[ev2,0]) : 0.0 : True
(ev2, 1) : 0.0 : soc[ev2,1] - (soc[ev2,0] + juice[ev2,1]) : 0.0 : True
(ev2, 2) : 0.0 : soc[ev2,2] - (soc[ev2,1] - 40 + juice[ev2,2]) : 0.0 : True
(ev2, 3) : 0.0 : soc[ev2,3] - (soc[ev2,2] + juice[ev2,3]) : 0.0 : True
11 Declarations: T E soc_index soc juice_index juice obj max_charge_index max_charge soc_balance_index soc_balance
Problem:
- Name: unknown
Lower bound: 57.2
Upper bound: 57.2
Number of objectives: 1
Number of constraints: 16
Number of variables: 16
Number of nonzeros: 3
Sense: minimize
Solver:
- Status: ok
User time: -1.0
System time: 0.0
Wallclock time: 0.0
Termination condition: optimal
Termination message: Model was solved to optimality (subject to tolerances), and an optimal solution is available.
Statistics:
Branch and bound:
Number of bounded subproblems: None
Number of created subproblems: None
Black box:
Number of iterations: 2
Error rc: 0
Time: 0.007892131805419922
Solution:
- number of solutions: 0
number of solutions displayed: 0
juice : juice added to E at time T
Size=8, Index=juice_index
Key : Lower : Value : Upper : Fixed : Stale : Domain
(ev1, 0) : 0 : 27.2 : None : False : False : NonNegativeReals
(ev1, 1) : 0 : 0.0 : None : False : False : NonNegativeReals
(ev1, 2) : 0 : 0.0 : None : False : False : NonNegativeReals
(ev1, 3) : 0 : 0.0 : None : False : False : NonNegativeReals
(ev2, 0) : 0 : 30.0 : None : False : False : NonNegativeReals
(ev2, 1) : 0 : 0.0 : None : False : False : NonNegativeReals
(ev2, 2) : 0 : 0.0 : None : False : False : NonNegativeReals
(ev2, 3) : 0 : 0.0 : None : False : False : NonNegativeReals
The state of charge in the vehicles:
ev1 {0: 37.2, 1: 7.2, 2: 7.2, 3: 0.0}
ev2 {0: 40.0, 1: 40.0, 2: 0.0, 3: 0.0}