Search code examples
pyomo

How to add piece-wise linear constraint in pyomo


In a task allocation problem, let's assume:

  1. we need to allocate 4 tasks to 4 days.

  2. There is only one worker and has to work at most 8 hours a day.

  3. We want to minimize the total working time.

  4. We can bundle tasks and save time in one day. The number of hours needed is

0 task --> 0 hours

1 task --> 6 hours

2 task --> 8 hours

An illustration for the # of Tasks V.S. # of Hours needed:

enter image description here

Blue dashed line: y = 6*n

Black line: Actual time = min(6*n, 6 + 2 * (n - 1))

Red line: Daily limit

This computation function (for the black line) can be defined in Python as

def compute_hours(n):
    return min(6*n, 6 + 2 * (n - 1))

Using the function defined above in the constraint setting, pyomo give errors when I put this in my constraint list.

> File "<ipython-input-69-57dbda13a71f>", line 23, in compute_hours
>     return min(6*n, 6 + 2 * (n - 1))   File "pyomo\core\expr\logical_expr.pyx", line 183, in
> pyomo.core.expr.logical_expr.InequalityExpression.__bool__
> pyomo.common.errors.PyomoException: Cannot convert non-constant
> expression to bool. This error is usually caused by using an
> expression in a boolean context such as an if statement. For example, 
>     m.x = Var()
>     if m.x <= 0:
>         ... would cause this exception.

The error is also easy to reproduce using the following sample code:

# Allocating tasks to days

from pyomo.environ import *
import pandas as pd

m = ConcreteModel()

# We have 4 tasks in total. Each cost 6 hours
m.TASKS = ('task1', 'task2', 'task3', 'task4')

# There are 4 days to allocate to
m.DAYS = (1, 2, 3, 4)

# Hours per task
HOURS_PER_TASK = 6


def compute_hours(n):
    return min(6*n, 6 + 2 * (n - 1))

# for x in range(0, 4):
#     print(x, compute_hours(x))

# Create a decision array for the allocation
m.ALLOCATION = Var(m.TASKS, m.DAYS, domain=Binary)

# Create an array for total hours computation
m.TOTAL_HOURS = Var(m.DAYS, domain=Integers)

# minimize the number of days that are allocated tasks
m.OBJ = Objective(expr=sum([m.TOTAL_HOURS[day] for day in m.DAYS]), sense=minimize)

m.c = ConstraintList()

# One task is done once only in all days
for task in m.TASKS:
    m.c.add(sum([m.ALLOCATION[task, day] for day in m.DAYS]) == 1)

# Compute Total hours a day
for day in m.DAYS:
    m.c.add(sum([m.ALLOCATION[task, day] for task in m.TASKS])*HOURS_PER_TASK == m.TOTAL_HOURS[day])
    # The following computation does not work
    #m.c.add(compute_hours(sum([m.ALLOCATION[task, day] for task in m.TASKS])*HOURS_PER_TASK) == m.TOTAL_HOURS[day])


# Add max working hours per day
for day in m.DAYS:
    m.c.add(m.TOTAL_HOURS[day] <= 8)

SolverFactory('glpk').solve(m).write()

# Show the results
SCHEDULE_df = pd.DataFrame(index = m.TASKS, columns= m.DAYS)
SCHEDULE_HOURS_df = pd.DataFrame(index = m.TASKS, columns= m.DAYS)

for i in m.TASKS:
    for j in m.DAYS:
        SCHEDULE_df.loc[i,j] = m.ALLOCATION[i,j]()

print('------------------------------------------')
print('ALLOCATION:')
print(SCHEDULE_df)

print('------------------------------------------')
print('Total hours per day:')
print(m.TOTAL_HOURS.get_values())
print('Total hours:')
print(value(m.OBJ))

This example can still give results if we assume each task always takes 6 hours.

------------------------------------------
ALLOCATION:
       1  2  3  4
task1  1  0  0  0
task2  0  1  0  0
task3  0  0  1  0
task4  0  0  0  1
------------------------------------------
Total hours per day:
{1: 6.0, 2: 6.0, 3: 6.0, 4: 6.0}
Total hours:
24.0

If we can realize the bundling, then this worker only needs to work for two days with workload 8 hours/day (instead of working 6 hours/day for 4 days).


Solution

  • John.

    Pyomo got a Piecewise Linear Expression for these cases. (kernel layer also got the Piecewise Function Library)

    In your case, it is easy, in fact, to use the Piecewise in Pyomo.

    You need to create an helper variable with the same index than m.TOTAL_HOURS (requested by pyomo), let's call it m.TASK_PER_DAY that will account for the number of task assignned for each day. This TASK_PER_DAY var will be then computed (constrained) with m.ALLOCATION using the summation of task each day (Similar as you try to constraint). Then, you can use the pyomo Piecewise to computed the TOTAL_HOURS depending upon TASK_PER_DAY

    Edit: Further explanation on pw_pts and f_rule, that are the source of error in modelung in Piecewise. pw_pts are the break points (independent values, domain) in the Piecewise. f_rule is the known value at given points pw_pts. This may be a function or a couple of points. If you want to linearly interpolate the points (as your image) you just need to define such points pw_pts=[0,1,2,3,...] and f_rule=[0,6,8,10,...], but for constant piecewise functions you need to specify the same point twice in order to ensure that for domain points, the piecewise returns the same value as in a step function. In your modeling, since variable are integers, there is not problem, but is better to clear the issue.

    m.TASK_PER_DAY = Var(m.DAYS, domain=Integers)
    m.PIECEWISE_DAILY_HOURS = Piecewise(
        m.DAYS,
        m.TOTAL_HOURS,
        m.TASK_PER_DAY,
        pw_pts=[0,1,2,3],
        f_rule=[0,6,8,10],
        pw_repn='SOS2',
        pw_constr_type='EQ',
        unbounded_domain_var=True
    )
    #Constraining TASK_PER_DAY
    def Assign_Task(model, day):
        return model.TASK_PER_DAY[day] == sum(model.ALLOCATION[task,day] for task in model.TASKS)
    m.assign_task = Constraint(m.DAYS, rule=Assign_Task)
    

    In this Approach you're not constraining the day to perform the task, therefore, the assignned day may be any pair of days, but each day certainly got 2 task per day. If you want to contraint the model to perform the task as soon (early days) as possible, you need to add that new constraint.

    The following code can reproduce the desired output:

    import pyomo.environ as pyo
    
    m = pyo.ConcreteModel()
    
    #Sets
    # We have 4 tasks in total. Each cost 6 hours
    m.TASKS = ('task1', 'task2', 'task3', 'task4')
    # There are 4 days to allocate to
    m.DAYS = (1, 2, 3, 4)
    # Hours per task
    HOURS_PER_TASK = 6
    
    #Variables
    # Create a decision array for the allocation
    m.ALLOCATION = pyo.Var(m.TASKS, m.DAYS, domain=Binary)
    m.TOTAL_HOURS = pyo.Var(m.DAYS, domain=Integers)
    m.TASK_PER_DAY = pyo.Var(m.DAYS, domain=Integers)
    
    #Piecewise
    m.PIECEWISE_DAILY_HOURS = pyo.Piecewise(
        m.DAYS,
        m.TOTAL_HOURS,
        m.TASK_PER_DAY,
        pw_pts=[0,0,1,1,2,2,3,3],
        f_rule=[0,0,6,6,8,8,10,10],
        pw_repn='SOS2',
        pw_constr_type='EQ',
        unbounded_domain_var=True
    )
    
    #Objective
    # minimize the number of days that are allocated tasks
    m.OBJ = pyo.Objective(expr=sum([m.TOTAL_HOURS[day] for day in m.DAYS]), sense=minimize)
    
    m.c = pyo.ConstraintList()
    
    # One task is done once only in all days
    for task in m.TASKS:
        m.c.add(sum([m.ALLOCATION[task, day] for day in m.DAYS]) == 1)
    
    #Compute the number of tasks per day
    def Assign_Task(model, day):
        return model.TASK_PER_DAY[day] == sum(model.ALLOCATION[task,day] for task in model.TASKS)
    
    m.assign_task = pyo.Constraint(m.DAYS, rule=Assign_Task)
    
    # Add max working hours per day
    for day in m.DAYS:
        m.c.add(m.TOTAL_HOURS[day] <= 8)
    
    pyo.SolverFactory('gurobi').solve(m, tee=True)
    

    Then the solution is reduced to 16 hours.

    >>>m.OBJ()
    16.0
    >>>m.TOTAL_HOURS.display()
    TOTAL_HOURS : Size=4, Index=TOTAL_HOURS_index
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :  None :  -0.0 :  None : False : False : Integers
          2 :  None :   8.0 :  None : False : False : Integers
          3 :  None :  -0.0 :  None : False : False : Integers
          4 :  None :   8.0 :  None : False : False : Integers
    >>>m.TASK_PER_DAY.display()
    TASK_PER_DAY : Size=4, Index=TASK_PER_DAY_index
        Key : Lower : Value : Upper : Fixed : Stale : Domain
          1 :  None :   0.0 :  None : False : False : Integers
          2 :  None :   2.0 :  None : False : False : Integers
          3 :  None :   0.0 :  None : False : False : Integers
          4 :  None :   2.0 :  None : False : False : Integers
    >>>m.ALLOCATION.display()
    ALLOCATION : Size=16, Index=ALLOCATION_index
        Key          : Lower : Value : Upper : Fixed : Stale : Domain
        ('task1', 1) :     0 :  -0.0 :     1 : False : False : Binary
        ('task1', 2) :     0 :   1.0 :     1 : False : False : Binary
        ('task1', 3) :     0 :  -0.0 :     1 : False : False : Binary
        ('task1', 4) :     0 :  -0.0 :     1 : False : False : Binary
        ('task2', 1) :     0 :  -0.0 :     1 : False : False : Binary
        ('task2', 2) :     0 :  -0.0 :     1 : False : False : Binary
        ('task2', 3) :     0 :  -0.0 :     1 : False : False : Binary
        ('task2', 4) :     0 :   1.0 :     1 : False : False : Binary
        ('task3', 1) :     0 :  -0.0 :     1 : False : False : Binary
        ('task3', 2) :     0 :   1.0 :     1 : False : False : Binary
        ('task3', 3) :     0 :  -0.0 :     1 : False : False : Binary
        ('task3', 4) :     0 :  -0.0 :     1 : False : False : Binary
        ('task4', 1) :     0 :  -0.0 :     1 : False : False : Binary
        ('task4', 2) :     0 :  -0.0 :     1 : False : False : Binary
        ('task4', 3) :     0 :  -0.0 :     1 : False : False : Binary
        ('task4', 4) :     0 :   1.0 :     1 : False : False : Binary