I am currently working on modeling a facility location-allocation using optimization with Pyomo module. For my model, I'd like to employ Piecewise linear function, with the constrained variables have two indices. I have tried three approaches that can be seen below, nonetheless none of them seems working properly (or having error)
Does anybody have any ideas on how to solve this?
piece_wise = {'ab': {0:1,1:1,2:2}, 'cd': {0:1,1:1,2:3}, 'ef': {0:1,1:2,2:4}}
piece_wise_ab = {0:1,1:1,2:2}
piece_wise_cd = {0:1,1:1,2:3}
piece_wise_ef = {0:1,1:2,2:4}
################# Option 1 - not working
model.y_ab = Var(model.Nodes, model.Vehicles, within = NonNegativeIntegers, bounds = (0, max_chargers))
model.arriv_ab = Var(model.Nodes, model.Vehicles, within = NonNegativeIntegers, bounds = (0, max_chargers))
model.y_cd = Var(model.Nodes, model.Vehicles, within = NonNegativeIntegers, bounds = (0, max_chargers))
model.arriv_cd = Var(model.Nodes, model.Vehicles, within = NonNegativeIntegers, bounds = (0, max_chargers))
model.y_ef = Var(model.Nodes, model.Vehicles, within = NonNegativeIntegers, bounds = (0, max_chargers))
model.arriv_ef = Var(model.Nodes, model.Vehicles, within = NonNegativeIntegers, bounds = (0, max_chargers))
model.piece_wise_ab = Piecewise(model.Nodes, model.Vehicles, model.arriv_ab, model.y_ab,
pw_pts = list(piece_wise_ab.values()), f_rule = list(piece_wise_ab.keys()), pw_repn = 'SOS2', pw_constr_type='EQ',unbounded_domain_var=True)
model.piece_wise_cd = Piecewise(model.Nodes, model.Vehicles, model.arriv_cd, model.y_cd,
pw_pts = list(piece_wise_cd.values()), f_rule = list(piece_wise_cd.keys()), pw_repn = 'SOS2', pw_constr_type='EQ',unbounded_domain_var=True)
model.piece_wise_ef = Piecewise(model.Nodes, model.Vehicles, model.arriv_cd, model.y_cd,
pw_pts = list(piece_wise_ef.values()), f_rule = list(piece_wise_ef.keys()), pw_repn = 'SOS2', pw_constr_type='EQ',unbounded_domain_var=True)
def service_rule (model,i,veh):
try:
if 'ab' in veh:
return sum(model.arrival [o,d,veh] * model.lam[i,o,d,veh] for o,d in model.sp) <= model.y_ab[i]
elif 'cd' in veh:
return sum(model.arrival [o,d,veh] * model.lam[i,o,d,veh] for o,d in model.sp) <= model.y_cd[i]
elif 'ef' in veh:
return sum(model.arrival [o,d,veh] * model.lam[i,o,d,veh] for o,d in model.sp) <= model.y_ef[i]
except:
return Constraint.Skip
model.service = Constraint(model.Nodes, model.Vehicles, rule = service_rule)
#################################### Option 2 - not working
model.piecewise_constraint = ConstraintList()
all_breakpoints = set(breakpoint for breakpoints in piece_wise.values() for breakpoint in breakpoints)
model.y_pc = Var(model.Nodes, model.Vehicles, within = NonNegativeIntegers, bounds = (0, max_chargers) )
model.arriv_pc = Var(model.Nodes, model.Vehicles, all_breakpoints, within = Binary)
for veh in model.Vehicles:
for node in model.Nodes:
breakpoints = list(piece_wise[veh].keys())
values_at_breakpoints = list(piece_wise[veh].values())
print(breakpoints, all_breakpoints)
model.piecewise_constraint.add(model.y_pc[node, veh] == sum(values_at_breakpoints[k] * model.arriv_pc[node, veh, k] for k in breakpoints))
############################## Option 3 - Not working
for node in model.Nodes:
for veh in model.Vehicles:
if "ab" in veh:
model.piecewise_constraint.add (model.y_pc [node, veh] == Piecewise (expr = [(model.arriv_pc [node, veh], model.y_pc [node, veh])],
pw_pts = list(piece_wise_ab.values()), f_rule = list(piece_wise_ab.keys()), pw_repn = 'SOS2',
pw_constr_type='EQ',unbounded_domain_var=True))
elif "cd" in veh:
model.piecewise_constraint.add (model.y_pc [node, veh] == Piecewise (expr = model.arriv_pc [node, veh],
pw_pts = list(piece_wise_cd.values()), f_rule = list(piece_wise_cd.keys()), pw_repn = 'SOS2',
pw_constr_type='EQ',unbounded_domain_var=True))
elif "ef" in veh:
model.piecewise_constraint.add (model.y_pc [node, veh] == Piecewise (expr = model.arriv_pc [node, veh],
pw_pts = list(piece_wise_ef.values()), f_rule = list(piece_wise_ef.keys()), pw_repn = 'SOS2',
pw_constr_type='EQ',unbounded_domain_var=True))
Thank you in Advance!
Here is an example that develops a piecewise linear constraint between variables x and y, both of which are dual-indexed variables.
You aren't following the dox located here for construction of the Piecewise constraint. Specifically, see the guidance on the f_rule.
from pyomo.environ import Objective, NonNegativeIntegers, ConcreteModel, Constraint, Var, Set, SolverFactory, \
sum_product, maximize, Piecewise
m = ConcreteModel('piecewise_example')
piecewise_elements = {0: 1, 5: 2, 10: 20}
m.Nodes = Set(initialize=[1, 2, 3])
m.Vehicles = Set(initialize=['a81', 'b22', 'c43'])
m.x = Var(m.Nodes, m.Vehicles, domain=NonNegativeIntegers, bounds=(0, 10))
m.y = Var(m.Nodes, m.Vehicles, domain=NonNegativeIntegers)
# max the sum of m.y...
m.obj = Objective(expr=sum_product(m.y), sense=maximize)
# subset the "c" cars...
m.c_vehicles = Set(initialize=[v for v in m.Vehicles if 'c' in v])
# limit the x vals based on node
@m.Constraint(m.Nodes, m.Vehicles)
def x_limit(m, node, vehicle):
if node == 1:
return m.x[node, vehicle] <= 4
return m.x[node, vehicle] <= 9
def limit_vals(model, node, vehicle, x_val):
# per the dox, the limit function needs to catch the model and all indices, but does
# not need to use them...
return piecewise_elements[x_val]
m.pw_limit = Piecewise(
m.Nodes,
m.Vehicles,
m.y,
m.x,
pw_pts=list(piecewise_elements.keys()),
f_rule=limit_vals,
pw_constr_type='UB')
m.pprint()
solver = SolverFactory('cbc')
res = solver.solve(m)
print(res)
for idx in m.x.index_set():
print(idx, m.x[idx].value, m.y[idx].value)
(1, 'a81') 4.0 1.0
(1, 'b22') 4.0 1.0
(1, 'c43') 4.0 1.0
(2, 'a81') 9.0 16.0
(2, 'b22') 9.0 16.0
(2, 'c43') 9.0 16.0
(3, 'a81') 9.0 16.0
(3, 'b22') 9.0 16.0
(3, 'c43') 9.0 16.0