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
- Activating Anaconda Environment in VsCode
- How to implement the Softmax function in Python?
- Dataframe - Select the minimum set of rows to cover all possible values of each columns
- Non-equi join in polars
- Tkinter and matplotlib - NavigationToolbar2Tk - cursor position precision
- Media Image is not loaded in the output of Django Project
- selenium upload part not working on windows. Path related?
- How to encode UTF8 filename for HTTP headers? (Python, Django)
- Split columns in Csv file
- In the Django admin site, how do I change the display format of time fields?
- Condition in pytest fixture depending on some test-suite results
- How to make a class JSON serializable
- Apache Airflow unit and integration test
- TypeError: 'MultiPolygon' object is not iterable with Plotly Choropleth
- How do I get the current time in milliseconds in Python?
- Python: Decoding base64 encoded strings
- How to fix the error when try to plot data with pandas?
- Python Image - Finding largest branch from image skeleton
- Enumerate is Returning A None Type
- Python finite difference functions?
- How to extract hex color codes from a colormap
- Regex: Get "London, UK" from "Evolution Recruitment (Agency) (London, UK)"
- Visual Studio Code - How to add multiple paths to python path?
- What are the Spark transformations that causes a Shuffle?
- How to find the second lowest lists into a nested list by their second value?
- python-docx: Parse a table to Panda Dataframe
- Deposit and Withdraw on Account
- Why does this code snippet give the error list object is not callable?
- How to find out if argparse argument has been actually specified on command line?
- Print current call stack from a method in code