False result about linear programming problem from pulp to gurobipy
This is my first time ask question on Stack Overflow. Recently I try to convert gurobipy code to pulp and meet some problem. Data resource:
- A list of employees with their available work hours, and hourly salaries:
- A prescibed minimum number of staff needed to be at work at a given hour (fewer staff are needed at night, more staff are needed during peak hours)
(ref: https://blogs.mathworks.com/loren/2016/01/06/generating-an-optimal-employee-work-schedule-using-integer-linear-programming/)
Goal:
- Minimize the total daily wages the employer must pay out.
Constraints:
- At any given hour, you must meet the minimum staffing requirement.
- An employee can only work one shift a day.
- An employee must work within his available hours.
- If an employee is called for duty, they must work at least a specified minimum number of hours, and no more than a specified maximum number of hours
Original gurobipy Code:
import gurobipy
# Create Data
EMPLOYEE, MIN, MAX, COST, START, END = gurobipy.multidict({
"SMITH" : [6, 8, 30, 6, 20], "JOHNSON": [6, 8, 50, 0, 24], 'WILLIAMS': [6, 8, 30, 0, 24],
'JONES' : [6, 8, 30, 0, 24], 'BROWN': [6, 8, 40, 0, 24], 'DAVIS': [6, 8, 50, 0, 24],
'MILLER' : [6, 8, 45, 6, 18], 'WILSON': [6, 8, 30, 0, 24], 'MOORE': [6, 8, 35, 0, 24],
'TAYLOR' : [6, 8, 40, 0, 24], 'ANDERSON': [2, 3, 60, 0, 6], 'THOMAS': [2, 4, 40, 0, 24],
'JACKSON' : [2, 4, 60, 8, 16], 'WHITE': [2, 6, 55, 0, 24], 'HARRIS': [2, 6, 45, 0, 24],
'MARTIN' : [2, 3, 40, 0, 24], 'THOMPSON': [2, 5, 50, 12, 24], 'GARCIA': [2, 4, 50, 0, 24],
'MARTINEZ': [2, 4, 40, 0, 24], 'ROBINSON': [2, 5, 50, 0, 24]})
REQUIRED = [1, 1, 2, 3, 6, 6, 7, 8, 9, 8, 8, 8, 7, 6, 6, 5, 5, 4, 4, 3, 2, 2, 2, 2]
MODEL = gurobipy.Model("Work Schedule")
# Create variable
x = MODEL.addVars(EMPLOYEE, range(24), range(1, 25), vtype=gurobipy.GRB.BINARY)
MODEL.update()
# Obj
MODEL.setObjective(gurobipy.quicksum((j - i) * x[d, i, j] * COST[d] for d in EMPLOYEE for i in range(24) for j in range(i + 1, 25)), sense=gurobipy.GRB.MINIMIZE)
# Constraints
MODEL.addConstrs(x.sum(d) <= (gurobipy.quicksum(x[d, i, j] for i in range(START[d], END[d] + 1) for j in range(i + 1, END[d] + 1) if MIN[d] <= j - i <= MAX[d])) for d in EMPLOYEE)
MODEL.addConstrs(gurobipy.quicksum(x[d, i, j] for i in range(START[d], END[d] + 1) for j in range(i + 1,END[d] + 1) if MIN[d] <= j - i <= MAX[d]) <= 1 for d in EMPLOYEE)
MODEL.addConstrs(gurobipy.quicksum(x[d, i, j] for d in EMPLOYEE for i in range(24) for j in range(i + 1, 25) if i <= c < j) >= REQUIRED[c] for c in range(24))
# Excute
MODEL.optimize()
My pulp Code:
# Create Data
# EMPLOYEE: 0)MIN, 1)MAX, 2)COST, 3)START, 4)END
md = {'SMITH':[6,8,30,6,20],'JOHNSON':[6,8,50,0,24],'WILLIAMS':[6,8,30,0,24],
'JONES':[6,8,30,0,24],'BROWN':[6,8,40,0,24],'DAVIS':[6,8,50,0,24],
'MILLER':[6,8,45,6,18],'WILSON':[6,8,30,0,24],'MOORE':[6,8,35,0,24],
'TAYLOR':[6,8,40,0,24],'ANDERSON':[2,3,60,0,6],'THOMAS':[2,4,40,0,24],
'JACKSON':[2,4,60,8,16],'WHITE':[2,6,55,0,24],'HARRIS':[2,6,45,0,24],
'MARTIN':[2,3,40,0,24],'THOMPSON':[2,5,50,12,24],'GARCIA':[2,4,50,0,24],
'MARTINEZ':[2,4,40,0,24],'ROBINSON':[2,5,50,0,24]}
REQUIRED = [1, 1, 2, 3, 6, 6, 7, 8, 9, 8, 8, 8, 7, 6, 6, 5, 5, 4, 4, 3, 2, 2, 2, 2]
MIN, MAX, COST, START, END = 0, 1, 2, 3, 4
t = 24
# Create Op
MODEL = pulp.LpProblem("Work Schedule", pulp.LpMinimize)
# Create variable
x = pulp.LpVariable.dicts("x", [(d, i, j) for d in md.keys() for i in range(t) for j in range(1, t+1)], cat=pulp.LpBinary)
# Create Obj
MODEL += pulp.lpSum((j - i) * x[d, i, j] * md[d][COST] for i in range(t) for j in range(i+1, t+1) for d in md.keys())
# Constraints
for c in range(len(REQUIRED)):
MODEL += pulp.lpSum(x[d, i, j] for d in md.keys() for i in range(t) for j in range(i+1, t+1) if i <= c < j) >= REQUIRED[c]
for k, v in md.items():
#print(k, v)
MODEL += pulp.lpSum(x[k, i, j] for i in range(v[START], v[END] + 1) for j in range(i + 1, v[END] + 1) if v[MIN] <= (j - i) <= v[MAX]) <= 1
MODEL += pulp.lpSum(x[k, i, j] for i in range(v[START], v[END] + 1) for j in range(i + 1, t + 1)) <= pulp.lpSum(x[k, i, j] for i in range(v[START], v[END] + 1) for j in range(i + 1, v[END] + 1) if v[MIN] <= (j - i) <= v[MAX])
MODEL.solve()
The gurobipy code result is:
My pulp code result is:
Obviously, my result is incorrect, but when I went through each condition individually, I couldn't find any differences. Does anyone know how to fix this? Thanks a lot, everyone.
Solution
I don't trust your results for either of the methods, because you've omitted ANDERSON's first shift half. (See original code's reference to cyclics and GALLERY).
I still don't think that you should have what amounts to O(n^2) variables. It's easy enough to constrain the employee-hour selection variables to be contiguous without that extra dimension.
import pulp
import pandas as pd
md = pd.DataFrame(
columns=('EmployeeName', 'MinHours', 'MaxHours', 'HourlyWage', 'Start', 'End'),
data=[
( 'SMITH', 6, 8, 30, 6, 20),
( 'JOHNSON', 6, 8, 50, 0, 24),
('WILLIAMS', 6, 8, 30, 0, 24),
( 'JONES', 6, 8, 30, 0, 24),
( 'BROWN', 6, 8, 40, 0, 24),
( 'DAVIS', 6, 8, 50, 0, 24),
( 'MILLER', 6, 8, 45, 6, 18),
( 'WILSON', 6, 8, 30, 0, 24),
( 'MOORE', 6, 8, 35, 0, 24),
( 'TAYLOR', 6, 8, 40, 0, 24),
('ANDERSON', 2, 3, 60, 18, 24 + 6),
( 'THOMAS', 2, 4, 40, 0, 24),
( 'JACKSON', 2, 4, 60, 8, 16),
( 'WHITE', 2, 6, 55, 0, 24),
( 'HARRIS', 2, 6, 45, 0, 24),
( 'MARTIN', 2, 3, 40, 0, 24),
('THOMPSON', 2, 5, 50, 12, 24),
( 'GARCIA', 2, 4, 50, 0, 24),
('MARTINEZ', 2, 4, 40, 0, 24),
('ROBINSON', 2, 5, 50, 0, 24),
],
).set_index(keys='EmployeeName', drop=True).sort_index()
# Hourly staffing requirements, over 24 hours only
hour48 = pd.RangeIndex(name='Hour', start=0, stop=48)
required = pd.Series(
name='Required', index=hour48[:24],
data=(1, 1, 2, 3, 6, 6, 7, 8, 9, 8, 8, 8, 7, 6, 6, 5, 5, 4, 4, 3, 2, 2, 2, 2),
)
# Long-style employee-hour selection variables, multi-indexed by employee and hour
flat = pd.DataFrame(
data=pulp.LpVariable.matrix(name='sel', indices=(md.index, hour48), cat=pulp.LpBinary),
index=md.index, columns=hour48,
).stack()
flat.name = 'Select'
# For each employee and available hour, each selection variable beside the employee data
lp_vars = pd.merge(
left=md, right=flat, on='EmployeeName',
).set_index(flat.index)
hour_idx = lp_vars.index.get_level_values('Hour')
lp_vars = lp_vars[
(hour_idx >= lp_vars['Start']) &
(hour_idx < lp_vars['End'])
]
# The expense for each employee, as affine expressions
md['Duration'] = lp_vars['Select'].groupby('EmployeeName').agg(pulp.lpSum)
md['Expense'] = md['HourlyWage'] * md['Duration']
model = pulp.LpProblem(name='Work_Schedule', sense=pulp.LpMinimize)
model.setObjective(pulp.lpSum(md['Expense']))
for name, row in md.iterrows():
model.addConstraint(
name=f'hours_min_{name}',
constraint=row['Duration'] >= row['MinHours'],
)
model.addConstraint(
name=f'hours_max_{name}',
constraint=row['Duration'] <= row['MaxHours'],
)
hour_idx = lp_vars.index.get_level_values('Hour')
for hour, min_staff in required.items():
hour_vars = lp_vars[hour_idx % 24 == hour]
# For each hour (or the same hour in the next day), enforce minimum staff
model.addConstraint(
name=f'min_staff_{hour}',
constraint=pulp.lpSum(hour_vars['Select']) >= min_staff,
)
for name, selects in lp_vars['Select'].groupby('EmployeeName'):
selects = selects.droplevel('EmployeeName')
# At the beginning and end of the day, within the minimum number of working hours, the
# selector values are simply monotonic
employee = md.loc[name]
nmin = employee['MinHours']
for s0, s1 in zip(selects[:nmin-1], selects[1: nmin]):
model.addConstraint(s0 <= s1)
for s0, s1 in zip(selects[-nmin: -1], selects[1-nmin:]):
model.addConstraint(s0 >= s1)
# In the middle of the day, conditionally enforce monotonicity based on selection
M = 48
for i in range(nmin-1, selects.size-nmin):
model.addConstraint(
name=f'm_{name}_{selects.index[i]}',
constraint=M*selects.iloc[i] <=
M*(selects.iloc[i+1] + 1) -
pulp.lpSum(selects.iloc[i+2: 1-nmin]),
)
print(model)
model.solve()
if model.status != pulp.LpStatusOptimal:
raise ValueError(model.status)
md[['Duration', 'Expense']] = md[['Duration', 'Expense']].map(pulp.LpAffineExpression.value)
lp_vars['Select'] = lp_vars['Select'].map(pulp.LpVariable.value)
pd.options.display.width = 200
pd.options.display.max_columns = 99
print(md)
hours_total = lp_vars['Select'].groupby('Hour').sum()
next_day = hours_total[required.size:]
hours_total[:next_day.size] += next_day.values
print(pd.DataFrame({
'required': required,
'allocated': hours_total[:required.size],
}))
print(
lp_vars['Select']
.astype(int)
.unstack('Hour', fill_value=0)
.sort_index(axis=1)
)
Work_Schedule:
MINIMIZE
60*sel_ANDERSON_18 + 60*sel_ANDERSON_19 + ...
SUBJECT TO
hours_min_ANDERSON: sel_ANDERSON_18 + sel_ANDERSON_19 + sel_ANDERSON_20
+ sel_ANDERSON_21 + sel_ANDERSON_22 + sel_ANDERSON_23 + sel_ANDERSON_24
+ sel_ANDERSON_25 + sel_ANDERSON_26 + sel_ANDERSON_27 + sel_ANDERSON_28
+ sel_ANDERSON_29 >= 2
...
VARIABLES
0 <= sel_ANDERSON_18 <= 1 Integer
0 <= sel_ANDERSON_19 <= 1 Integer
0 <= sel_ANDERSON_20 <= 1 Integer
...
Welcome to the CBC MILP Solver
Version: 2.10.3
Build Date: Dec 15 2019
At line 2 NAME MODEL
At line 3 ROWS
At line 467 COLUMNS
At line 5951 RHS
At line 6414 BOUNDS
At line 6833 ENDATA
Problem MODEL has 462 rows, 418 columns and 4229 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 4700 - 0.01 seconds
...
Result - Optimal solution found
Objective value: 4700.00000000
...
MinHours MaxHours HourlyWage Start End Duration Expense
EmployeeName
ANDERSON 2 3 60 18 30 2.0 120.0
BROWN 6 8 40 0 24 8.0 320.0
DAVIS 6 8 50 0 24 8.0 400.0
GARCIA 2 4 50 0 24 4.0 200.0
HARRIS 2 6 45 0 24 6.0 270.0
JACKSON 2 4 60 8 16 2.0 120.0
JOHNSON 6 8 50 0 24 8.0 400.0
JONES 6 8 30 0 24 8.0 240.0
MARTIN 2 3 40 0 24 3.0 120.0
MARTINEZ 2 4 40 0 24 4.0 160.0
MILLER 6 8 45 6 18 8.0 360.0
MOORE 6 8 35 0 24 8.0 280.0
ROBINSON 2 5 50 0 24 3.0 150.0
SMITH 6 8 30 6 20 8.0 240.0
TAYLOR 6 8 40 0 24 8.0 320.0
THOMAS 2 4 40 0 24 4.0 160.0
THOMPSON 2 5 50 12 24 5.0 250.0
WHITE 2 6 55 0 24 2.0 110.0
WILLIAMS 6 8 30 0 24 8.0 240.0
WILSON 6 8 30 0 24 8.0 240.0
required allocated
Hour
0 1 1.0
1 1 1.0
2 2 2.0
3 3 3.0
4 6 6.0
5 6 6.0
6 7 7.0
7 8 8.0
8 9 9.0
9 8 8.0
10 8 8.0
11 8 8.0
12 7 7.0
13 6 6.0
14 6 6.0
15 5 5.0
16 5 5.0
17 4 4.0
18 4 4.0
19 3 3.0
20 2 2.0
21 2 2.0
22 2 2.0
23 2 2.0
Hour 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
EmployeeName
ANDERSON 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0
BROWN 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0
DAVIS 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
GARCIA 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
HARRIS 0 0 0 0 0 0 0 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
JACKSON 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
JOHNSON 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
JONES 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0
MARTIN 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
MARTINEZ 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
MILLER 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0
MOORE 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0
ROBINSON 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0
SMITH 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
TAYLOR 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
THOMAS 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
THOMPSON 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0
WHITE 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
WILLIAMS 0 0 0 0 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
WILSON 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
- AttributeError: install_layout when attempting to install a package in a virtual environment
- Python list comprehension - want to avoid repeated evaluation
- Hash algorithm for dynamic growing/streaming data?
- matplotlib - making labels for violin plots
- Python How to I check if last element has been reached in iterator tool chain?
- Polars and the Lazy API: How to drop columns that contain only null values?
- Why are my Mean, Var, and Std outputs from NumPy different from what the online grader expects?
- Correlation dataframe convertion from results from pl.corr
- Polars DataFrame transformation
- Discord rate limiting while only sending 1 request per minute
- Check if column contains (/,-,_, *or~) and split in another column - Pandas
- How to draw a rectangle at (x,y) in a PyQt GraphicsView?
- how to calculate correlation between ten columns with polars
- How to set class attribute with await in __init__
- Detect hindi encoding, response received from Facebook API in Python
- Is it possible to write a horizontal if statement with a multi-line body?
- Max length of items in list
- Cannot subclass multiprocessing Queue in Python 3.5
- How can I get notified of updates to Python packages in a unified way?
- Using python AST to traverse code and extract return statements
- merge groups of columns in a polars dataframe to single columns
- Group Pandas DataFrame by Continuous Date Ranges
- Flask login @login_required not working
- Odoo: one2many and many2one? KeyError:'___'
- merge some columns in a Polars dataframe and duplicate the others
- Python: Create table from string mixed with separators using FOR loops
- How do I type hint a method with the type of the enclosing class?
- How can I verify an emails DKIM signature in Python?
- Writing a class that accepts a callback in Python?
- Python Paramiko channel.exec_command not returning output intermittently