Piecewise linear functions with multiple indices variables Pyomo

Question

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!

Answer 1

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.

Code

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)

Partial Output:

(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