用开源求解器和Pyomo实现灵活的班次安排

前段时间讲课时发现简单的优化问题求解还是有些业务需求的，但是无论使用商用还是开源求解器，手工建模仍然需要一定工作量。所以这几天写了个基于pyomo和开源求解器glpk的小排程程序，可以实现多人、自定义班次和排程规则的最优化班次排程(shift assignment)，优化目标包含最低成本和员工最大偏好，并且支持自定义的目标权重。程序支持excel或txt表格形式的标准化输入和输出。

业务需求示例

业务部门要求将n位员工对应至m种班次，在满足每周需求，并尽量满足员工工作时间偏好的同时使人力成本最小。排程单位为小时，planning horizon为一周（168小时）。

业务部门有预先定义的每日班次记录在表格shift里，每个班次预设时长在8-12小时之间，但理论上可以用小于24小时的任何班次设置，也可以跨不同天。

业务部门有n名员工，每位员工每周的工作时间必须满足最小工时/班次数和最大工时/班次数的要求。每位员工的时薪和输出都不同。另外每位员工可以提出对一周内所有班次的偏好，记录在表格preference里。preference为0的班次原则上不安排给员工。每位员工一周内每个小时的availability另外记录在表格availability里。这张表的值为0的时间段也不会安排员工工作（比如休假）。

业务部门定义了每小时的人力需求记录在表格demand里。但每小时的需求不需要严格满足；最终只要求排程满足每周的总体需求，而每小时的需求只需要满足minDemand和maxDemand区间：这是为了保证每小时capacity的上下限。

表格regulation记录所有其他关于班次的要求，目前包含：

每人每周最大/最小工作小时数/班次数
每人的相邻两班次最少相隔小时数
每周最少相同相邻班次个数（保证工作时段稳定性）
如果有夜班，每周最少和最多相邻夜班个数（保证工作时段稳定性）
属于夜班班次的时间范围
基础时薪，以及夜班、周末、法定节假日加班薪酬比例

此外还有关于多目标权重的设定：

排程考虑员工偏好的权重
排程考虑运营成本的权重

输入格式

Demand:

字段	示例	描述
publicHoliday	0	1为公共假日，0为非公共假日
weeks	1	周index，默认只有一周排程
days	1	1-7,周一为1，周日为7
hours	1	从1开始递增，表示排程中的第n小时
hourPerDay	1	1-24,1代表第0-1小时，24代表第23-24小时
demand	2	绝对需求，以小时为单位
min	1	最小capacity
max	4	最大capacity

排程时考虑每周完成绝对demand的总量，每小时保证在最低和最高capacity之间。

Shift：

字段	示例	描述
weeks	1	周index，默认只有一周排程
days	1	1-7,周一为1，周日为7
hours	1	从1开始递增，表示排程中的第n小时
hourPerDay	1	1-24,1代表第0-1小时，24代表第23-24小时
shiftA	1	ID为shiftA的占用时间，1为占用，0为不占用
shiftB	0	ID为shiftB的占用时间，1为占用，0为不占用

只包含1天内的班次，默认一周内每天的可能班次相同。在排程时会考虑跨半夜和跨周的班次连接。另外根据demand，每个班次可能排多人，也可能不会使用。

24小时shift示例：

weeks	days	hours	hourPerDay	shiftA	shiftB	shiftC	shiftD	shiftE	shiftF	shiftG	shiftH	shiftI	shiftJ	shiftK
1	1	1	1	1	0	0	1	0	0	1	0	0	0	0
1	1	2	2	1	0	0	1	0	0	1	0	0	0	0
1	1	3	3	1	0	0	1	0	0	1	0	0	0	0
1	1	4	4	1	0	0	1	0	0	1	0	0	1	0
1	1	5	5	1	0	0	1	0	0	1	0	0	1	0
1	1	6	6	1	0	0	1	0	0	1	0	0	1	0
1	1	7	7	0	1	0	1	1	0	1	0	1	1	0
1	1	8	8	0	1	0	1	1	0	1	0	1	1	0
1	1	9	9	0	1	0	0	1	0	1	0	1	1	0
1	1	10	10	0	1	0	0	1	0	1	0	1	1	0
1	1	11	11	0	1	0	0	1	0	1	0	1	1	0
1	1	12	12	0	1	0	0	1	0	1	0	1	0	0
1	1	13	13	0	1	0	0	1	0	0	1	1	0	0
1	1	14	14	0	1	0	0	1	0	0	1	1	0	0
1	1	15	15	0	0	1	0	1	1	0	1	1	0	1
1	1	16	16	0	0	1	0	1	1	0	1	1	0	1
1	1	17	17	0	0	1	0	0	1	0	1	1	0	1
1	1	18	18	0	0	1	0	0	1	0	1	1	0	1
1	1	19	19	0	0	1	0	0	1	0	1	0	0	1
1	1	20	20	0	0	1	0	0	1	0	1	0	0	1
1	1	21	21	0	0	1	0	0	1	0	1	0	0	1
1	1	22	22	0	0	1	0	0	1	0	1	0	0	1
1	1	23	23	1	0	0	1	0	1	0	1	0	0	0
1	1	24	24	1	0	0	1	0	1	0	1	0	0	0

Availability:

字段	示例	描述
publicHoliday	0	1为公共假日，0为非公共假日
weeks	1	周index，默认只有一周排程
days	1	1-7,周一为1，周日为7
hours	1	从1开始递增，表示排程中的第n小时
hourPerDay	1	1-24,1代表第0-1小时，24代表第23-24小时
A	1	ID为A的员工的availability,1代表available,0代表休假
B	1	ID为B的员工的availability,1代表available,0代表休假

Preference 员工期望排程：

员工ID必须对应availability的所有员工。

班次ID必须对应一周内的所有班次，数量为shift表格里班次数量 * 排程天数。

Preference数值越大代表偏好越明显。如果preference设置为0，则排程时视为不可安排。

字段	示例	描述
employees	A	员工ID
employeeHourlyWage	1	员工时薪
employeeHourlyOutput	1	员工效率
employeePriority	7	员工排程优先级，越高越优先满足preference
shiftA1	1	员工对周一班次shiftA的preference，数值越高preference越高。默认取值范围0-3，但没有特殊要求
shiftB1	3	员工对周一班次shiftB的preference
shiftA2	1	员工对周二班次shiftA的preference
shiftB2	3	员工对周二班次shiftB的preference

Regulation 排程规则：

规则名称	示例	单位	描述
minWorkHour	35	hours	每周每人最少工作小时数（如果有假期会自动排除）
maxWorkHour	60	hours	每周每人最大工作小时数
minShiftsPerWeek	3	shifts	每周每人最少工作班次数（如果有假期会自动排除）
maxShiftsPerWeek	5	shifts	每周每人最大工作班次数
minHourBetweenShift	12	hours	每人的相邻两班次最少相隔小时数
minShiftContinuous	2	days	每周最少相同相邻班次个数（保证工作时段稳定性）
minNightShiftContinuous	3	days	如果有夜班，每周最少相邻夜班个数（保证工作时段稳定性）
nightShiftDefinitionStart	24	hour	属于夜班班次的时间范围（起始时间，不包含）。班次中有任何时间段落入范围，整个班次为夜班
nightShiftDefinitionEnd	4	hour	属于夜班班次的时间范围（结束时间，包含）
maxNightShiftContinuous	5	days	如果有夜班，每周最多相邻夜班个数
standardShiftCostPerPersonPerHour	1	unit	基础时薪（以1为单位，1代表100%）
standardShiftPaymentStart	0	hour	基础时薪的时间范围（起始时间，包含）
standardShiftPaymentEnd	24	hour	基础时薪的时间范围（结束时间，包含）
additionalNightShiftCost	1	unit	夜班加薪（基础时薪的百分比）
additionalNightShiftPaymentStart	22	hour	夜班加薪的时间范围（起始时间，不包含），可以不同于夜班班次时间范围，默认大于
additionalNightShiftPaymentEnd	6	hour	夜班加薪的时间范围（结束时间，包含）
additionalWeekendShiftCost	1	unit	周末加薪（基础时薪的百分比）
additionalWeekendDayStart	6	Saturday	周末时间范围（起始日，包含）
additionalWeekendDayEnd	7	Sunday	周末时间范围（结束日，包含）
additionalPublicHolidayCost	1	unit	公共假日加薪（基础时薪的百分比）
weightPreference	0.1	unit	计算排程目标时，考虑员工preference的权重，取值在0-1之间
weightCost	0.9	unit	计算排程目标时，考虑运营成本的权重，取值在0-1之间

数据检查、优化目标和约束

检查【min(员工available时间,规定最长工作时间) * 员工单位产出】是否大于等于【客户需求demand总量】。如果最大产出都不能满足需求总量，在进行优化排程前会调整demand到产出的上线。
排程优化会使用以下硬性约束（i.e.如果不能满足，就没有优化结果）：
- 每周最大/最小工作小时数。如果员工有安排假期，则相应向下调整最小工作小时数
- 每周最大/最小班次数。如果员工有安排假期，则相应向下调整最小班次数
- 相邻班次间最小间隔小时数
- 每周最多连续夜班数
- 员工availability设为0，和员工preference设为0的位置不可以安排班次
- 每周总产出必须大于等于总需求
- 每小时的总产出比如在同一时间的最大和最小需求限制以内（比如如果有产能限制）
另外会考虑逻辑或的约束：
- 每周最少连续相同/相似班次数。或员工完全不安排班次，就视为合理
- 每周最少连续相同/相似夜班次数。或员工完全不安排夜班，就视为合理
注意：非凸优化问题可能产生局部最优解；尤其是逻辑或约束个数很多时，结果可能只部分满足约束条件。
排程优化目标：
- 最小化运营成本（根据标准时薪、夜班、周末、公共节假日的设定，计算不同时间段的成本）
- 最大化员工期望的班次数
两个目标的平衡是通过设定weightPreference 和 weightCost两个参数实现的。在数据处理过程中会对所有目标相关的数据做归一化处理。

如果设置最小班次数和最小工作小时数同时为0，在满足所有硬性约束的条件下，有可能有员工没有被安排任何班次（根据成本优先安排已经有工时的员工）。

如果有要求所有员工都必须满足最小班次数和最小工作小时数，可以通过设置最小班次数和最小工作小时数中的一项。

其他数据处理

shift表格只包含24小时的班次信息。会将24小时的班次数据扩展到7天（168小时）。对于跨0时的班次，程序会自动做班次拼接的处理。
计算互斥班次：程序会计算所有不满足班次最小间隔的互斥班次。在计算过程中会自动考虑跨越周日和周一的情况（i.e. 默认所有周都是一周的重复循环）。
计算最大连续夜班班次个数：会生成所有的连续夜班组合，这些组合都比规定的最大连续夜班个数多1个班次。在排程过程中，员工的排程中间结果会与这些夜班组合逐一比较，并保证不出现与任一夜班组合重合的情况。
计算最小连续相同班次个数：会生成所有的连续相同/相似班次组合，这些组合都恰好包含规定的最小连续班次个数。在排程过程中，员工的排程中间结果会与这些连续相同班次组合逐一比较，并保证至少与其中一个组合重合。

对于相似班次的处理：入参中可以设置considerSimilarShifts=True, 同时设置“相似”班次的标准：similarShiftHours有多少个小时重合的情况下视为相似班次。
计算最小连续相同夜班个数：同4.

建模简述

除了需要逻辑或(Logical-OR)的约束，其他约束和目标都满足线性条件，可以使用标准线性建模方法。

需要逻辑或的约束条件是通过convex hull实现的。建模使用了pyomo中的Generalized Disjunctive Programming (GDP)，具体参考这里。

由于逻辑或是non-convex约束，大量的逻辑或会增加模型的计算量，而且并不保证能够找到可行解。在建模时把逻辑或约束放在其他约束之前，可以在一定程度上提高解的质量。

另外通过createCombinations函数会生成hourlyOutput的不同组合作为模型的超参数输入。不同组合会生成不同的解。程序会根据逻辑或约束的满足程度对解排序，在一定程度上提高解的质量。程序的执行效率通过进程池加速。

hourlyOutput组合是所有满足hourlyOutput小于等于给定output的值的组合。所以给定各个员工的hourlyOutput越高或员工数量越大，组合数量就越大，遍历所有组合需要的时间也越长。之所以选择hourlyOutput作为可变化的超参数，是因为默认所有小于给定hourlyOutput的值也都属于合理范围，等价于牺牲成本换取对约束的满足。

如果仍然出现结果不满足逻辑或约束的情况，可以尝试放宽demand的上下限。

详细建模说明：将业务需求转化为目标和约束

决策变量：

设优化问题的决策变量为 x，这个矩阵的维度是【员工人数n * 班次个数m】，取值为0或1。

优化求解的目的是找到一组x，使目标函数（偏好-成本）达到最大值。

约束：

1… 最大maxshift和最小班次数minshift的约束可以表示为：

∑jmxij>=minshift,i∈1,2,⋯ ,n\sum_j^m x_{ij} >= minshift, i \in {1,2, \cdots, n}j∑mxij>=minshift,i∈1,2,⋯,n

∑jmxij<=maxshift,i∈1,2,⋯ ,n\sum_j^m x_{ij} <= maxshift, i \in {1,2, \cdots, n}j∑mxij<=maxshift,i∈1,2,⋯,n

2… 给定每个班次所在的时间段为向量：

shiftj,j∈1,2,⋯ ,mshift_j, j \in {1,2, \cdots, m}shiftj,j∈1,2,⋯,m

每个向量的长度为 t = 7*24 = 168。

则有每位员工对应的工作时间为：

onDutyik=∑jmxij∗shiftjk,i∈1,2,⋯ ,n,k∈1,2,⋯ ,tonDuty_{ik} = \sum_j^m x_{ij} * shift_{jk}, i \in {1,2, \cdots, n}, k \in {1,2, \cdots, t}onDutyik=j∑mxij∗shiftjk,i∈1,2,⋯,n,k∈1,2,⋯,t

onDuty的取值应为0或1。

3… 最大maxhour工作时长和最小工作时长minhour的约束可以表示为：

∑ktonDutyik>=minhour,i∈1,2,⋯ ,n\sum_k^t onDuty_{ik} >= minhour, i \in {1,2, \cdots, n}k∑tonDutyik>=minhour,i∈1,2,⋯,n

∑ktonDutyik<=maxhour,i∈1,2,⋯ ,n\sum_k^t onDuty_{ik} <= maxhour, i \in {1,2, \cdots, n}k∑tonDutyik<=maxhour,i∈1,2,⋯,n

4… 给定avail为给定的每位员工每周可工作时段，取值为0或1，维度为【员工数n * 每周小时数t】；取值为0的时段为休假时段。

避免在availability为0的时段安排班次：

∑in∑ktonDutyik∗(1−availik)=0\sum_i^n \sum_k^t onDuty_{ik} * (1-avail_{ik}) = 0i∑nk∑tonDutyik∗(1−availik)=0

5… 给定pref为每位员工对每个班次的偏好矩阵，维度为【员工数n * 班次数m】，取值为0的是完全不接受的班次；定义noGoPref为：

noGoPrefij={1if prefij=00elsenoGoPref_{ij} = \begin{cases} 1 &\quad \text{if } pref_{ij}=0 \\ 0 &\quad \text{else} \end{cases}noGoPrefij={10if prefij=0else

避免在偏好为0的班次安排工作：

∑in∑jmnoGoPrefij∗xij=0\sum_i^n \sum_j^m noGoPref_{ij} * x_{ij} = 0i∑nj∑mnoGoPrefij∗xij=0

6… 给定每位员工的工作产出为长度为n的向量hourlyOutput，给定需求为长度为168的向量demand，上下限为minDemand/maxDemand：

每个时段的产出在上下限之间：

∑inonDutyik∗hourlyOutputi>=minDemand,k∈1,2,⋯ ,t\sum_i^n onDuty_{ik} * hourlyOutput_{i} >= minDemand, k \in {1,2, \cdots, t}i∑nonDutyik∗hourlyOutputi>=minDemand,k∈1,2,⋯,t

∑inonDutyik∗hourlyOutputi<=maxDemand,k∈1,2,⋯ ,t\sum_i^n onDuty_{ik} * hourlyOutput_{i} <= maxDemand, k \in {1,2, \cdots, t}i∑nonDutyik∗hourlyOutputi<=maxDemand,k∈1,2,⋯,t

每周总产出大于总需求：

∑kt∑inonDutyik∗hourlyOutputi>=∑ktdemandk\sum_k^t \sum_i^n onDuty_{ik} * hourlyOutput_{i} >= \sum_k^t demand_kk∑ti∑nonDutyik∗hourlyOutputi>=k∑tdemandk

7… 相邻班次间隔必须大于12小时：

首先构造不符合要求的所有班次组合。用于构造不符合要求班次组合的函数在DataProblem.createNoGoShifts。输出的结果noGoShiftGroups是维度为【组合个数S * 班次个数m】的给定矩阵，取值为0或1。每一行上同时取值为1的对应班次是不符合规定的班次组合。

避免不符合规定的班次组合：

∑jmnoGoShiftGroupssj∗xij<=1,i∈1,2,⋯ ,n,s∈1,2,⋯ ,S\sum_j^m noGoShiftGroups_{sj} * x_{ij} <= 1, i \in {1,2, \cdots, n}, s \in {1,2, \cdots, S}j∑mnoGoShiftGroupssj∗xij<=1,i∈1,2,⋯,n,s∈1,2,⋯,S

共生成 S * n 个约束，保证每位员工的排程都避免了所有S个非法组合。

8… 同样方法避免不符合最长连续夜班规定的夜班组合：设最大允许连续夜班数为maxNights。构造不符合最大夜班规定的夜班组合的函数在DataProblem.createNoGoNightShifts，生成维度为【组合个数R * 班次个数m】的给定矩阵：

∑jmnoGoShiftGroupsrj∗xij<=maxNights,i∈1,2,⋯ ,n,r∈1,2,⋯ ,R\sum_j^m noGoShiftGroups_{rj} * x_{ij} <= maxNights, i \in {1,2, \cdots, n}, r \in {1,2, \cdots, R}j∑mnoGoShiftGroupsrj∗xij<=maxNights,i∈1,2,⋯,n,r∈1,2,⋯,R

共生成 R * n 个约束，保证每位员工的排程都避免了所有R个非法组合。

logical_OR逻辑或约束：

在pyomo5以前，涉及决策变量的“逻辑或”可以通过bigM或convex hull的方式实现，在求解前需要调用pyomo.gdp.bigm或pyomo.gdp.chull。在pyomo5以后，增加了CMU贡献的gdpopt基于逻辑的求解器。两种方法的建模过程类似。本程序选择使用gdpopt。

9… 最小连续班次个数约束：首先构造符合条件的向量：用于构造最小连续班次的函数在DataProblem.createContinuousShifts.输出为contShifts矩阵，维度为【组合数Q * 班次个数m】，取值为0或1，每行同时取值为1的对应班次为一组可能的组合。

定义最小连续班次个数为minContShifts。

定义一组约束：

continuousiq=∑jmcontShiftsqj∗xij,i∈1,2,⋯ ,n,q∈1,2,⋯ ,Qcontinuous_{iq} = \sum_j^m contShifts_{qj} * x_{ij}, i \in {1,2, \cdots, n}, q \in {1,2, \cdots, Q} continuousiq=j∑mcontShiftsqj∗xij,i∈1,2,⋯,n,q∈1,2,⋯,Q

一组对应的disjuct block的定义为：

continuousiq>=minContShiftscontinuous_{iq} >= minContShiftscontinuousiq>=minContShifts

以及“不安排任何班次”的n个约束：

∑jmxij=0,i∈1,2,⋯ ,n\sum_j^m x_{ij}=0, i \in {1,2, \cdots, n}j∑mxij=0,i∈1,2,⋯,n

这组不等式共有n*Q+n个，每一位员工有一组Q+1个约束，要求每一组Q+1个不等式至少有一个成立（结果为True)，也就是每人每周至少有minContShifts个连续且重复的班次，或者0个班次，保证员工工作时间的稳定性。

10… 最小连续夜班个数约束同理：给定维度为【组合数W * 班次个数m】，取值为0或1的矩阵contNightShifts。设最小连续夜班个数为minContNightShifts。

定义一组约束：

continuousNightiw=∑jmcontNightShiftswj∗xij,i∈1,2,⋯ ,n,w∈1,2,⋯ ,WcontinuousNight_{iw} = \sum_j^m contNightShifts_{wj} * x_{ij}, i \in {1,2, \cdots, n}, w \in {1,2, \cdots, W} continuousNightiw=j∑mcontNightShiftswj∗xij,i∈1,2,⋯,n,w∈1,2,⋯,W

一组对应的disjuct block的定义为：

continuousNightiw>=minContNightShiftscontinuousNight_{iw} >= minContNightShiftscontinuousNightiw>=minContNightShifts

以及“不安排任何夜班”的n个约束：

∑wWcontinuousNightiw=0,i∈1,2,⋯ ,n\sum_w^W continuousNight_{iw}=0, i \in {1,2, \cdots, n}w∑WcontinuousNightiw=0,i∈1,2,⋯,n

这组不等式共有n*W+n个，每一位员工有一组W+1个约束，要求每一组W+1个不等式至少有一个成立（结果为True)，也就是每人每周至少有minContNightShifts个连续且重复的夜班，或者0个夜班，保证员工工作时间的稳定性。

目标：

11… 根据标准班次、夜班、周末班次、节假日班次的设置，给定每小时成本为长度=168的向量cost。给定每个员工的时薪为长度=n的向量hourlyWage。

为了鼓励尽量安排最少人数（优先安排已经有班次的员工），增加惩罚项penalty（为确定值）。

则有成本约束cost per person：

costppi=∑ktonDutyik∗costi∗hourlyWagei+onDutyik∗penalty,i∈1,2,⋯ ,ncostpp_i = \sum_k^t onDuty_{ik} * cost_i * hourlyWage_i + onDuty_{ik} * penalty, i \in {1,2, \cdots, n}costppi=k∑tonDutyik∗costi∗hourlyWagei+onDutyik∗penalty,i∈1,2,⋯,n

12… 偏好约束pref per person:

prefppi=∑jmxij∗prefij,i∈1,2,⋯ ,nprefpp_i = \sum_j^m x_{ij} * pref_{ij}, i \in {1,2, \cdots, n}prefppi=j∑mxij∗prefij,i∈1,2,⋯,n

将多目标用权重转化成单一目标：

obj=∑inprefWeight∗prefppi−costWeight∗costppiobj = \sum_i^n prefWeight * prefpp_i - costWeight * costpp_iobj=i∑nprefWeight∗prefppi−costWeight∗costppi

输出示例：

shift assignment：员工 vs. 班次的透视表，1为安排班次，0为未安排班次。

示例：

employee	shiftB1	shiftC1	shiftG1	shiftH1	shiftI1	shiftK1
A	0	0	0	0	0	1
B	1	0	0	0	0	0
C	1	0	0	0	0	0
D	0	0	0	0	1	0
E	0	1	0	0	0	0
F	0	0	0	1	0	0
G	0	0	1	0	0	0

schedule：员工 vs. 小时的透视表，1为上班，0为休息。

示例：

weeks	days	hours	hourPerDay	B	C	D	G
1	1	1	1	0	0	0	1
1	1	2	2	0	0	0	1
1	1	3	3	0	0	0	1
1	1	4	4	0	0	0	1
1	1	5	5	0	0	0	1
1	1	6	6	0	0	0	1
1	1	7	7	1	1	1	1
1	1	8	8	1	1	1	1
1	1	9	9	1	1	1	1
1	1	10	10	1	1	1	1

代码实现：

数据准备部分过于啰嗦，这里不再赘述。建模部分如下：

class ShiftModel():def __init__(self, modelData,solverChoice='gdpopt',mip='glpk',nlp='ipopt',solverpath_exe=None):self.solverChoice = solverChoiceself.mip = mip # mixed integer programming solverself.nlp = nlp # nonlinear programming solverself.solverpath_exe = solverpath_exeself.model = self.buildConcreteModel(modelData)def buildConcreteModel(self, data):model = ConcreteModel()model.Idx_employee = range(len(data.employees))model.Idx_shift = range(data.shift.shape[1])model.Idx_noGoShiftGroups = range(data.noGoShiftGroups.shape[0])model.Idx_noGoNightShifts = range(data.noGoNightShifts.shape[0])model.Idx_hour = range(data.shift.shape[0])model.Idx_demand = range(len(data.demand.demand))model.x = Var(model.Idx_employee, model.Idx_shift, within=Binary, initialize=1)# the minimum continuous shifts per week turns the problem into non-convex # and would require getting many logical-OR's via disjunction or a big-M relaxation,# which are numerically terrible. if np.sum(data.continuousShifts)>0:model.Idx_contShifts = range(data.continuousShifts.shape[0])model.contShifts = Var(model.Idx_employee, model.Idx_contShifts, initialize=0)        def calContShifts_rule(model, Idx_employee, Idx_contShifts):return sum(data.continuousShifts[Idx_contShifts, i] * model.x[Idx_employee, i]for i in model.Idx_shift) == model.contShifts[Idx_employee, Idx_contShifts]model.contShifts_c = Constraint(model.Idx_employee, model.Idx_contShifts, rule=calContShifts_rule)model.disjuncts = list()model.disjuncts_c = list()for Idx_employee in model.Idx_employee:model.disjuncts.append(list())for Idx_contShifts in model.Idx_contShifts:model.disjuncts[Idx_employee].append(Disjunct())model.disjuncts[Idx_employee][Idx_contShifts].c = Constraint(expr=model.contShifts[Idx_employee, Idx_contShifts]>=int(data.regDict['minShiftContinuous']))# the disjunct of not-any-shift-scheduled:model.disjuncts[Idx_employee].append(Disjunct())model.disjuncts[Idx_employee][-1].c = Constraint(expr=sum(model.x[Idx_employee,i] for i in model.Idx_shift)==0)model.disjuncts_c.append(Disjunction(expr=model.disjuncts[Idx_employee]))# minimum continuous night shifts disjuncts are created in a similar way# except that there is the additional disjunct of having no night shifts scheduled at all.if np.sum(data.continuousNightShifts)>0:model.Idx_contNightShifts = range(data.continuousNightShifts.shape[0])model.contNightShifts = Var(model.Idx_employee, model.Idx_contNightShifts, initialize=0)        def calContNightShifts_rule(model, Idx_employee, Idx_contNightShifts):return sum(data.continuousNightShifts[Idx_contNightShifts, i] * model.x[Idx_employee, i]for i in model.Idx_shift) == model.contNightShifts[Idx_employee, Idx_contNightShifts]model.contNightShifts_c = Constraint(model.Idx_employee, model.Idx_contNightShifts, rule=calContNightShifts_rule)        model.nightDisjuncts = list()model.nightDisjuncts_c = list()for Idx_employee in model.Idx_employee:model.nightDisjuncts.append(list())for Idx_contNightShifts in model.Idx_contNightShifts:model.nightDisjuncts[Idx_employee].append(Disjunct())model.nightDisjuncts[Idx_employee][Idx_contNightShifts].c = Constraint(expr=model.contNightShifts[Idx_employee, Idx_contNightShifts]>=int(data.regDict['minNightShiftContinuous']))                model.nightDisjuncts[Idx_employee].append(Disjunct())model.nightDisjuncts[Idx_employee][-1].c = Constraint(expr=sum(model.contNightShifts[Idx_employee,i] for i in model.Idx_contNightShifts)==0)model.nightDisjuncts_c.append(Disjunction(expr=model.nightDisjuncts[Idx_employee]))if np.sum(data.noGoShiftGroups.shape)>0:def noGoShiftGroups_rule(model, Idx_noGoShiftGroups, Idx_employee):return sum(data.noGoShiftGroups[Idx_noGoShiftGroups,i] * model.x[Idx_employee,i] for i in model.Idx_shift) <= 1model.noGoShiftGroups_c = Constraint(model.Idx_noGoShiftGroups, model.Idx_employee, rule=noGoShiftGroups_rule)if np.sum(data.noGoNightShifts)>0:def noGoNightShifts_rule(model, Idx_noGoNightShifts, Idx_employee):return sum(data.noGoNightShifts[Idx_noGoNightShifts,i] * model.x[Idx_employee,i]for i in model.Idx_shift) <= int(data.regDict['maxNightShiftContinuous'])model.noGoNightShifts_c = Constraint(model.Idx_noGoNightShifts, model.Idx_employee, rule=noGoNightShifts_rule)def minMaxShifts_rule(model, Idx_employee):minshift = max([data.regDict['minShiftsPerWeek'] - int(sum(data.nonAvail[Idx_employee]) / 24), 0])return (minshift, sum(model.x[Idx_employee,i] for i in model.Idx_shift), data.regDict['maxShiftsPerWeek']) model.maxShifts_c = Constraint(model.Idx_employee, rule=minMaxShifts_rule)model.onDuty = Var(model.Idx_employee, model.Idx_hour, within=Binary, initialize=1)def calOnDuty_rule(model, Idx_employee, Idx_hour):shift = np.array(data.shift.iloc[Idx_hour,:], dtype=int)return sum(model.x[Idx_employee,i] * shift[i] for i in model.Idx_shift) == model.onDuty[Idx_employee, Idx_hour]model.onDuty_c = Constraint(model.Idx_employee, model.Idx_hour, rule=calOnDuty_rule)def onDuty_rule(model, Idx_employee):minhour = max([data.regDict['minWorkHour'] - sum(data.nonAvail[Idx_employee]), 0])return (minhour,sum(model.onDuty[Idx_employee, i] for i in model.Idx_hour), data.regDict['maxWorkHour'])model.onDuty_bnd = Constraint(model.Idx_employee, rule=onDuty_rule)if np.sum(data.nonAvail)>0:model.nonAvail_c = ConstraintList()for Idx_employee in model.Idx_employee:if sum(data.nonAvail[Idx_employee]) > 0:lhs = sum(model.onDuty[Idx_employee,i] * data.nonAvail[Idx_employee,i] for i in model.Idx_hour)model.nonAvail_c.add(lhs == 0)model.Idx_noGoPref = np.array([idx for idx,i in enumerate(np.sum(data.noGoPreference, axis=1).astype(bool).astype(int).tolist()) if i>0])def noGoPreference_rule(model, Idx_noGoPref):return sum(data.noGoPreference[Idx_noGoPref,i] * model.x[Idx_noGoPref,i] for i in model.Idx_shift)==0if sum(model.Idx_noGoPref)>0:model.noGoPref_c = Constraint(model.Idx_noGoPref, rule=noGoPreference_rule)def demandHour_rule(model, Idx_hour):return (data.demand.minDemand[Idx_hour],sum(model.onDuty[i,Idx_hour] * data.hourlyOutput[i] for i in model.Idx_employee), data.demand.maxDemand[Idx_hour])model.demandHour_c = Constraint(model.Idx_hour, rule=demandHour_rule)model.demandTotal_c = Constraint(expr=sum(model.onDuty[i,j] * data.hourlyOutput[i]for i in model.Idx_employee for j in model.Idx_hour) >= sum(data.demand.demand))model.costpp = Var(model.Idx_employee, bounds=(0,max(data.cost)*sum(data.demand.maxDemand)))def costpp_rule(model, Idx_employee):return sum(model.onDuty[Idx_employee,i] * data.cost[i] * data.hourlyWage[Idx_employee]for i in model.Idx_hour) == model.costpp[Idx_employee]model.costpp_c = Constraint(model.Idx_employee, rule=costpp_rule)model.prefpp = Var(model.Idx_employee, bounds=(0,np.max(data.preference)*sum(data.demand.maxDemand)))def prefpp_rule(model, Idx_employee):return sum(model.x[Idx_employee,i] * data.preference[Idx_employee,i] for i in model.Idx_shift) == model.prefpp[Idx_employee]model.prefpp_c = Constraint(model.Idx_employee, rule=prefpp_rule)def obj_rule(model):return sum(data.regDict['weightPreference'] * model.prefpp[i]- data.regDict['weightCost'] * model.costpp[i] for i in model.Idx_employee)model.obj = Objective(rule=obj_rule, sense=maximize)return modeldef runOptimizer(self,data):self.model.preprocess()if self.solverChoice == 'gdpopt':opt.solve(self.model, mip=self.mip, mip_options={'tmlim':30}, nlp=self.nlp, tol=1E-8, iterlim=30, tee=True)else:opt.solve(self.model, tee=True)x = list(self.model.x.get_values().values())self.result = xreturn x