blog

ナーススケジューリング(2)

Published:

By nob

Category: Posts

Tags: 数理最適化 Python-MIP Python

前回「ナーススケジューリング」 の続き。

看護師の勤務表をプログラムで作成する。

  • 日数: 30日
  • 看護師: 28名
  • シフト: 2交替制（日勤・夜勤・その他セミナー等・休日）

前回うまくスケジュール出来なかったのでモデルを見直す。

制約条件を読み込む

制約条件の詳細は 「ナーススケジューリング」 を参照のこと。 

import pandas as pd
from IPython.display import display

pd.options.display.max_columns = 40
pd.options.display.max_rows = 100

shift_constraints = pd.read_csv("data/pism/53-2-231-table-11.csv")
shift_constraints["day"] = shift_constraints["day"] - 1
nurse_constraints = pd.read_csv("data/pism/53-2-231-table-12.csv")
nurse_constraints["nurse_no"] = nurse_constraints["nurse_no"] - 1
last_shift_table_temp = pd.read_csv("data/pism/53-2-231-table-13-1.csv")
last_shift_table_temp["nurse_no"] = last_shift_table_temp["nurse_no"] - 1
last_shift_table = (
    last_shift_table_temp.set_index("nurse_no").stack().reset_index()
)
last_shift_table.columns = ["staff", "day", "shift"]
# 当月初日をday0とする。前月は31日まである。
last_shift_table["day"] = last_shift_table["day"].astype(int) - 32
last_shift_table = pd.get_dummies(
    last_shift_table,
    columns=["shift"],
    dtype=int,
    prefix="",
    prefix_sep="",
)
last_shift_table = (
    last_shift_table.set_index(["staff", "day"]).stack().reset_index()
)
last_shift_table.columns = ["staff", "day", "shift", "var"]
desired_shifts_temp = pd.read_csv("data/pism/53-2-231-table-13-2.csv")
desired_shifts_temp["nurse_no"] = desired_shifts_temp["nurse_no"] - 1
desired_shifts = (
    pd.DataFrame(desired_shifts_temp.drop("nurse_no", axis=1).unstack(level=1))
    .dropna()
    .reset_index()
)
desired_shifts.columns = ["day", "staff", "shift"]
desired_shifts["day"] = desired_shifts["day"].astype(int) - 1
desired_shifts = desired_shifts.reindex(
    columns=["staff", "day", "shift"]
).sort_values(["staff", "day", "shift"])

モデルを作成する

制約条件からモデルを作成する。目的関数は前回同様シフトの実現度とした。

変数を作成する

モデルの変数を表形式で作成する。

ナースNo. 日付 シフト その日のシフトかどうか(0/1の2値変数)
0 0 - 0
0 0 N 0
0 0 n 0
0 0 / 0
0 0 + 1
0 1 - 0
... ... ... ...

計算を容易にする為、モデルの変数に制約条件の一部を結合しておく。 

import psutil
from mip import Model, OptimizationStatus, maximize, xsum

model = Model(name="nurse_scheduling", solver_name="cbc")
model.threads = psutil.cpu_count(logical=False)

staffs = nurse_constraints["nurse_no"]
staffs.name = "staff"
shift_days = shift_constraints["day"]
shift_days.name = "day"
shifts = pd.Series(["-", "N", "n", "/", "+"], name="shift")

num_staffs = staffs.shape[0]
num_days = shift_days.shape[0]
num_shifts = shifts.shape[0]

# 当月のみの制約条件を作成する為に用いる
shift_table = pd.merge(pd.merge(staffs, shift_days, "cross"), shifts, "cross")
shift_table["var"] = model.add_var_tensor(
    (len(shift_table),), "var", var_type="B"
)

shift_table = shift_table.join(
    nurse_constraints.set_index("nurse_no"),
    on="staff",
    how="left",
)

# 前月・当月を通して制約条件を作成する為に用いる
total_shift_table = pd.concat(
    [
        last_shift_table.join(
            nurse_constraints.set_index("nurse_no"),
            on="staff",
            how="left",
        ),
        shift_table,
    ],
    ignore_index=True,
)

shift_table = shift_table.join(
    shift_constraints.set_index("day"),
    on="day",
    how="left",
)

制約条件を作成する

全員にシフトが割り当てられていること。 

for _, shift in shift_table.groupby(["staff", "day"]):
    model += xsum(shift["var"]) == 1

禁止パターンをアサインしないこと。 

forbidden_patterns = [
    # 夜勤の翌日は夜勤明けでなければならない
    ["N", "-"],
    ["N", "/"],
    ["N", "+"],
    # 夜勤明けの翌日は休みでなければならない
    ["n", "-"],
    ["n", "n"],
    ["n", "N"],
    # 全ての拘束条件を満たしているとされる勤務表の例(表1)を見ると
    # ナース#25の5日からのシフトでは夜勤明け, その他の勤務が許容されている。
    # 勤務の希望(表13)を見るとこれは本人の希望のようである。
    ["n", "+"],
    ["n", "+", "-"],
    ["n", "+", "N"],
    ["n", "+", "n"],
    ["n", "+", "+"],
    ["-", "n"],
    ["/", "n"],
    ["+", "n"],
    # 日勤, 休み, 日勤, 休み, 日勤というパターンは許されていない
    ["-", "/", "-", "/", "-"],
    # 4日連続の日勤は許されていない
    ["-", "-", "-", "-"],
]

total_shift_days = total_shift_table["day"].unique()
total_shift_days.sort()


def sliding_window(alist, window):
    for i in range(0, len(alist) - window + 1, 1):
        yield alist[i : i + window]


# 前月・当月通しての制約
for _, groupby_staff in total_shift_table.groupby("staff"):

    # 1週間に1回は休みが入らなくてはならない
    for days_in_window in sliding_window(total_shift_days, 7):
        model += (
            xsum(
                groupby_staff[
                    (groupby_staff["day"].isin(days_in_window))
                    & (groupby_staff["shift"] == "/")
                ]["var"]
            )
            >= 1
        )

    # 禁止パターン
    for forbidden_pattern in forbidden_patterns:
        window = len(forbidden_pattern)
        for days_in_window in sliding_window(total_shift_days, window):
            model += (
                xsum(
                    [
                        groupby_staff[
                            (groupby_staff["day"] == day)
                            & (groupby_staff["shift"] == shift)
                        ]["var"].values[0]
                        for day, shift in zip(
                            days_in_window, forbidden_pattern
                        )
                    ]
                )
                <= window - 1
            )

    # 夜勤と夜勤の間は少なくとも3日あけなければならない
    # GOOD: N n / - - N n /
    # NG: N n / - N n /
    # 夜勤明けを含めて5日間で計算する
    for days_in_window in sliding_window(total_shift_days, 5):
        model += (
            xsum(
                groupby_staff[
                    (groupby_staff["day"].isin(days_in_window))
                    & (groupby_staff["shift"] == "N")
                ]["var"]
            )
            <= 1
        )

ナース毎の日勤・夜勤・休日の日数 

# 当月のみ、ナース毎の制約
for _, groupby_staff in shift_table.groupby("staff"):
    # 日勤の上限・下限
    for shift_name, shift_code in zip(
        ["day", "night", "off"], ["-", "N", "/"]
    ):
        lower_bound = groupby_staff.at[
            groupby_staff.index[0], "{}_lower_bound".format(shift_name)
        ]
        upper_bound = groupby_staff.at[
            groupby_staff.index[0], "{}_upper_bound".format(shift_name)
        ]
        staff_sum = xsum(
            groupby_staff[groupby_staff["shift"] == shift_code]["var"]
        )
        model += staff_sum >= lower_bound
        model += staff_sum <= upper_bound

日毎に日勤・夜勤のそれぞれについて以下の制約を追加する。

  • 全体の人数
  • ベテランの人数
  • ベテランまたは2年目の人数
  • 新人の人数
  • チーム毎の人数 
team_shift_skill_levels = [
    [
        "all_day",
        ["A", "B", "C"],
        ["skilled", "second_year", "recently_qualified"],
        "-",
    ],
    [
        "all_night",
        ["A", "B", "C"],
        ["skilled", "second_year", "recently_qualified"],
        "N",
    ],
    ["s_day", ["A", "B", "C"], ["skilled"], "-"],
    ["s_night", ["A", "B", "C"], ["skilled"], "N"],
    ["rq_day", ["A", "B", "C"], ["recently_qualified"], "-"],
    ["rq_night", ["A", "B", "C"], ["recently_qualified"], "N"],
    ["a_ss_day", ["A"], ["skilled", "second_year"], "-"],
    ["a_ss_night", ["A"], ["skilled", "second_year"], "N"],
    ["b_ss_day", ["B"], ["skilled", "second_year"], "-"],
    ["b_ss_night", ["B"], ["skilled", "second_year"], "N"],
    ["c_ss_day", ["C"], ["skilled", "second_year"], "-"],
    ["c_ss_night", ["C"], ["skilled", "second_year"], "N"],
    ["a_day", ["A"], ["skilled", "second_year", "recently_qualified"], "-"],
    ["a_night", ["A"], ["skilled", "second_year", "recently_qualified"], "N"],
    ["b_day", ["B"], ["skilled", "second_year", "recently_qualified"], "-"],
    ["b_night", ["B"], ["skilled", "second_year", "recently_qualified"], "N"],
    ["c_day", ["C"], ["skilled", "second_year", "recently_qualified"], "-"],
    ["c_night", ["C"], ["skilled", "second_year", "recently_qualified"], "N"],
]

for day, groupby_day in shift_table.groupby("day"):
    for shift_name, team, skill_level, shift_code in team_shift_skill_levels:
        lower_bound = groupby_day.at[
            groupby_day.index[0], "{}_lower_bound".format(shift_name)
        ]
        upper_bound = groupby_day.at[
            groupby_day.index[0], "{}_upper_bound".format(shift_name)
        ]
        num_assigned_staffs = xsum(
            groupby_day[
                (groupby_day["team"].isin(team))
                & (groupby_day["skill_level"].isin(skill_level))
                & (groupby_day["shift"] == shift_code)
            ]["var"]
        )
        model += num_assigned_staffs >= lower_bound
        model += num_assigned_staffs <= upper_bound

各ナース, 少なくとも1回は週末2連休を必要とする。(2-3日, 3-4日, 9-10日, 16-17日, 23-24日の中から選ぶ.4日は祭日)

週末2連休かどうかを表す変数( holiday_constraints )を追加する。

ナースNo. 2連休の候補 2連休かどうか(0/1の2値変数)
0 0 (2-3) 0
0 1 (3-4) 0
0 2 (9-10) 1
0 3 (16-17) 0
0 4 (23-24) 0
1 0 (2-3) 1
1 1 (3-4) 0
... ... ...

週末2連休かどうかを表す変数( holiday_constraints )をシフトの変数( shift_table )と関連付けることを考える。

少なくとも1回は週末2連休なので、その2日間について休日かどうか(0/1の2値変数)を合計すると以下のパターンとなる。

  • 2連休である: 合計は2
  • 2連休でない: 合計は0以上(どちらが1日が休日である可能性がある) 
holiday_in_rows = [
    [day - 1 for day in holiday_in_row]
    for holiday_in_row in [
        [2, 3],
        [3, 4],
        [9, 10],
        [16, 17],
        [23, 24],
    ]
]

holiday_constraints = model.add_var_tensor(
    (
        num_staffs,
        len(holiday_in_rows),
    ),
    "holiday_constraints",
    var_type="B",
)

# 少なくとも1回は週末2連休である
for staff, groupby_staff in shift_table.groupby("staff"):
    model += xsum(holiday_constraints[staff]) >= 1

# 週末2連休となった場合、その2日はどちらも休日である
for staff, groupby_staff in shift_table.groupby("staff"):
    for i, holiday_in_row in enumerate(holiday_in_rows):
        model += (
            xsum(
                groupby_staff[
                    (groupby_staff["day"].isin(holiday_in_row))
                    & (groupby_staff["shift"] == "/")
                ]["var"]
            )
            # どちらか1日が休日である可能性がある
            >= 2 * holiday_constraints[staff, i]
        )

目的関数を作成する 

# シフトの希望が最大限叶うようにする
shift_satisfaction = desired_shifts.join(
    shift_table.set_index(["day", "staff", "shift"]),
    on=["day", "staff", "shift"],
    how="left",
)

# 夜勤以外を希望 => 日勤・夜勤明け・休日
for index, row in shift_satisfaction[
    shift_satisfaction["shift"].isin(["*"])
].iterrows():
    shift_satisfaction = pd.concat(
        [
            shift_satisfaction,
            shift_table[
                (shift_table["day"] == row["day"])
                & (shift_table["staff"] == row["staff"])
                & (shift_table["shift"].isin(["-", "n", "+", "/"]))
            ],
        ],
        ignore_index=True,
    )

# 日勤以外を希望 => 夜勤・夜勤明け・休日
for index, row in shift_satisfaction[
    shift_satisfaction["shift"].isin(["x"])
].iterrows():
    shift_satisfaction = pd.concat(
        [
            shift_satisfaction,
            shift_table[
                (shift_table["day"] == row["day"])
                & (shift_table["staff"] == row["staff"])
                & (shift_table["shift"].isin(["N", "n", "+", "/"]))
            ],
        ],
        ignore_index=True,
    )

shift_satisfaction = shift_satisfaction[
    shift_satisfaction["shift"].isin(["-", "N", "n", "/", "+"])
]

model.objective = maximize(xsum(shift_satisfaction["var"]))

解を求める 

model.optimize()

Welcome to the CBC MILP Solver 
Version: devel 
Build Date: Aug  4 2024
Starting solution of the Linear programming relaxation problem using Primal Simplex

Coin0506I Presolve 15669 (-1707) rows, 4187 (-153) columns and 58302 (-5006) elements
Clp0030I 32 infeas 0.21554593, obj 81.968868 - mu 0.01691816, its 52, 3968 interior
Clp0030I 50 infeas 0.00064545043, obj 81.952026 - mu 2.3193432e-05, its 62, 4055 interior
Clp1000I sum of infeasibilities 5.20571e-05 - average 3.3223e-09, 132 fixed columns
Coin0506I Presolve 14721 (-948) rows, 4033 (-154) columns and 55380 (-2922) elements
Clp0006I 0  Obj 81.951915 Primal inf 1.0300511e-05 (17) Dual inf 7.7e+13 (142)
Clp0029I End of values pass after 4033 iterations
Clp0014I Perturbing problem by 0.001% of 1.5170141 - largest nonzero change 2.9904382e-05 ( 0.0014952191%) - largest zero change 2.9819115e-05
Clp0000I Optimal - objective value 82
Clp0000I Optimal - objective value 82
Coin0511I After Postsolve, objective 82, infeasibilities - dual 0 (0), primal 0 (0)
Clp0006I 0  Obj 82 Dual inf 2.3333303 (3)
Clp0014I Perturbing problem by 0.001% of 1.5095009 - largest nonzero change 2.9958224e-05 ( 0.0014979112%) - largest zero change 2.9693113e-05
Clp0000I Optimal - objective value 82
Clp0000I Optimal - objective value 82
Clp0000I Optimal - objective value 82
Coin0511I After Postsolve, objective 82, infeasibilities - dual 0 (0), primal 0 (0)
Clp0032I Optimal objective 82 - 0 iterations time 2.272, Presolve 0.04, Idiot 2.22

Starting MIP optimization
threads was changed from 0 to 14
maxSavedSolutions was changed from 1 to 10
Continuous objective value is 82 - 0.011495 seconds
Cgl0002I 116 variables fixed
Cgl0003I 59 fixed, 0 tightened bounds, 12098 strengthened rows, 1482 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 11623 strengthened rows, 26 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 7298 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 4285 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 1175 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 477 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 210 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 117 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 74 strengthened rows, 0 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 57 strengthened rows, 0 substitutions
Cgl0004I processed model has 8157 rows, 3387 columns (3387 integer (3387 of which binary)) and 55561 elements
Coin3009W Conflict graph built in 0.019 seconds, density: 0.072%
Cgl0015I Clique Strengthening extended 8 cliques, 25 were dominated
After applying Clique Strengthening continuous objective value is 82
Cbc0045I Nauty did not find any useful orbits in time 0.905898
Cbc0038I Initial state - 530 integers unsatisfied sum - 137.353
Cbc0038I Pass   1: (26.08 seconds) suminf.   37.89683 (126) obj. 80.4722 iterations 4573
Cbc0038I Pass   2: (26.83 seconds) suminf.   23.33333 (74) obj. 79.6667 iterations 981
Cbc0038I Pass   3: (27.64 seconds) suminf.   15.66667 (53) obj. 79.5 iterations 985
Cbc0038I Pass   4: (28.21 seconds) suminf.    9.83333 (35) obj. 79 iterations 691
Cbc0038I Pass   5: (28.29 seconds) suminf.    9.83333 (35) obj. 79 iterations 81
Cbc0038I Pass   6: (28.52 seconds) suminf.    8.66667 (31) obj. 79 iterations 233
Cbc0038I Pass   7: (28.67 seconds) suminf.    8.66667 (31) obj. 79 iterations 163
Cbc0038I Pass   8: (28.87 seconds) suminf.    9.41667 (31) obj. 79 iterations 198
Cbc0038I Pass   9: (29.11 seconds) suminf.    8.66667 (31) obj. 79 iterations 265
Cbc0038I Pass  10: (29.44 seconds) suminf.    5.33333 (14) obj. 79 iterations 346
Cbc0038I Pass  11: (29.71 seconds) suminf.    4.33333 (12) obj. 79 iterations 274
Cbc0038I Pass  12: (29.96 seconds) suminf.    4.33333 (12) obj. 79 iterations 269
Cbc0038I Pass  13: (30.18 seconds) suminf.    4.33333 (12) obj. 79 iterations 224
Cbc0038I Pass  14: (30.39 seconds) suminf.    4.33333 (12) obj. 79 iterations 215
Cbc0038I Pass  15: (31.67 seconds) suminf.    4.33333 (12) obj. 76 iterations 1296
Cbc0038I Pass  16: (33.00 seconds) suminf.    4.33333 (12) obj. 76 iterations 1113
Cbc0038I Pass  17: (33.09 seconds) suminf.    4.33333 (12) obj. 76 iterations 51
Cbc0038I Pass  18: (34.70 seconds) suminf.    4.33333 (12) obj. 74 iterations 1323
Cbc0038I Pass  19: (36.23 seconds) suminf.    3.66667 (10) obj. 74 iterations 1199
Cbc0038I Pass  20: (36.32 seconds) suminf.    3.66667 (10) obj. 74 iterations 54
Cbc0038I Pass  21: (36.45 seconds) suminf.    3.66667 (10) obj. 74 iterations 89
Cbc0038I Pass  22: (36.78 seconds) suminf.    3.66667 (10) obj. 74 iterations 233
Cbc0038I Pass  23: (38.79 seconds) suminf.    5.00000 (14) obj. 73 iterations 1588
Cbc0038I Pass  24: (40.50 seconds) suminf.    4.33333 (12) obj. 73 iterations 1335
Cbc0038I Pass  25: (40.83 seconds) suminf.    3.66667 (10) obj. 73 iterations 219
Cbc0038I Pass  26: (41.07 seconds) suminf.    3.66667 (10) obj. 73 iterations 165
Cbc0038I Pass  27: (41.13 seconds) suminf.    3.66667 (10) obj. 73 iterations 27
Cbc0038I Pass  28: (41.20 seconds) suminf.    3.66667 (10) obj. 73 iterations 37
Cbc0038I Pass  29: (41.26 seconds) suminf.    3.66667 (10) obj. 73 iterations 32
Cbc0038I Pass  30: (43.10 seconds) suminf.    4.33333 (12) obj. 73.6667 iterations 1468
Cbc0038I No solution found this major pass
Cbc0038I Before mini branch and bound, 2706 integers at bound fixed and 21 continuous
Cbc0038I Full problem 8140 rows 3387 columns, reduced to 995 rows 309 columns
Cbc0038I Mini branch and bound did not improve solution (43.31 seconds)
Cbc0038I Full problem 8141 rows 3387 columns, reduced to 8137 rows 3387 columns - too large
Cbc0038I After 47.20 seconds - Feasibility pump exiting - took 24.27 seconds
Cbc0031I 20 added rows had average density of 18.65
Cbc0013I At root node, 20 cuts changed objective from 82 to 82 in 7 passes
Cbc0014I Cut generator 0 (Probing) - 44 row cuts average 2.2 elements, 0 column cuts (0 active)  in 0.417 seconds - new frequency is 1
Cbc0014I Cut generator 1 (Gomory) - 66 row cuts average 323.8 elements, 0 column cuts (0 active)  in 0.944 seconds - new frequency is 1
Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.252 seconds - new frequency is -100
Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.093 seconds - new frequency is -100
Cbc0014I Cut generator 4 (OddWheel) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.901 seconds - new frequency is -100
Cbc0014I Cut generator 5 (MixedIntegerRounding2) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.102 seconds - new frequency is -100
Cbc0014I Cut generator 6 (FlowCover) - 0 row cuts average 0.0 elements, 0 column cuts (0 active)  in 0.007 seconds - new frequency is -100
Cbc0014I Cut generator 7 (TwoMirCuts) - 193 row cuts average 160.4 elements, 0 column cuts (0 active)  in 0.886 seconds - new frequency is -100
Cbc0014I Cut generator 8 (ZeroHalf) - 22 row cuts average 29.1 elements, 0 column cuts (0 active)  in 3.348 seconds - new frequency is -100
Cbc0010I After 0 nodes, 1 on tree, -1e+50 best solution, best possible 82 (63.47 seconds)
Cbc0010I After 1 nodes, 2 on tree, -1e+50 best solution, best possible 82 (67.14 seconds)
Cbc0010I After 2 nodes, 1 on tree, -1e+50 best solution, best possible 82 (70.08 seconds)
Cbc0010I After 3 nodes, 2 on tree, -1e+50 best solution, best possible 82 (71.94 seconds)
Cbc0010I After 4 nodes, 2 on tree, -1e+50 best solution, best possible 82 (75.10 seconds)
Cbc0010I After 5 nodes, 2 on tree, -1e+50 best solution, best possible 82 (75.87 seconds)
Cbc0010I After 6 nodes, 1 on tree, -1e+50 best solution, best possible 82 (78.62 seconds)
Cbc0010I After 8 nodes, 2 on tree, -1e+50 best solution, best possible 82 (80.38 seconds)
Cbc0010I After 9 nodes, 1 on tree, -1e+50 best solution, best possible 82 (81.91 seconds)
Cbc0010I After 11 nodes, 2 on tree, -1e+50 best solution, best possible 82 (84.32 seconds)
Cbc0010I After 14 nodes, 1 on tree, -1e+50 best solution, best possible 82 (85.85 seconds)
Cbc0010I After 17 nodes, 1 on tree, -1e+50 best solution, best possible 82 (87.07 seconds)
Cbc0010I After 18 nodes, 2 on tree, -1e+50 best solution, best possible 82 (89.11 seconds)
Cbc0010I After 19 nodes, 2 on tree, -1e+50 best solution, best possible 82 (90.30 seconds)
Cbc0010I After 21 nodes, 3 on tree, -1e+50 best solution, best possible 82 (91.17 seconds)
Cbc0010I After 25 nodes, 5 on tree, -1e+50 best solution, best possible 82 (92.52 seconds)
Cbc0010I After 30 nodes, 7 on tree, -1e+50 best solution, best possible 82 (94.01 seconds)
Cbc0010I After 31 nodes, 8 on tree, -1e+50 best solution, best possible 82 (95.19 seconds)
Cbc0010I After 32 nodes, 9 on tree, -1e+50 best solution, best possible 82 (95.91 seconds)
Cbc0010I After 34 nodes, 10 on tree, -1e+50 best solution, best possible 82 (96.66 seconds)
Cbc0010I After 35 nodes, 10 on tree, -1e+50 best solution, best possible 82 (97.55 seconds)
Cbc0010I After 38 nodes, 12 on tree, -1e+50 best solution, best possible 82 (98.50 seconds)
Cbc0010I After 44 nodes, 16 on tree, -1e+50 best solution, best possible 82 (99.85 seconds)
Cbc0010I After 47 nodes, 19 on tree, -1e+50 best solution, best possible 82 (101.03 seconds)
Cbc0010I After 48 nodes, 20 on tree, -1e+50 best solution, best possible 82 (102.02 seconds)
Cbc0010I After 50 nodes, 21 on tree, -1e+50 best solution, best possible 82 (102.81 seconds)
Cbc0010I After 51 nodes, 21 on tree, -1e+50 best solution, best possible 82 (103.58 seconds)
Cbc0010I After 55 nodes, 24 on tree, -1e+50 best solution, best possible 82 (105.06 seconds)
Cbc0010I After 61 nodes, 26 on tree, -1e+50 best solution, best possible 82 (107.07 seconds)
Cbc0010I After 64 nodes, 28 on tree, -1e+50 best solution, best possible 82 (108.25 seconds)
Cbc0010I After 65 nodes, 28 on tree, -1e+50 best solution, best possible 82 (109.40 seconds)
Cbc0010I After 67 nodes, 29 on tree, -1e+50 best solution, best possible 82 (110.12 seconds)
Cbc0010I After 71 nodes, 32 on tree, -1e+50 best solution, best possible 82 (110.85 seconds)
Cbc0010I After 76 nodes, 36 on tree, -1e+50 best solution, best possible 82 (112.25 seconds)
Cbc0010I After 79 nodes, 38 on tree, -1e+50 best solution, best possible 82 (114.00 seconds)
Cbc0010I After 81 nodes, 40 on tree, -1e+50 best solution, best possible 82 (115.06 seconds)
Cbc0010I After 84 nodes, 41 on tree, -1e+50 best solution, best possible 82 (115.82 seconds)
Cbc0010I After 86 nodes, 43 on tree, -1e+50 best solution, best possible 82 (116.63 seconds)
Cbc0010I After 88 nodes, 44 on tree, -1e+50 best solution, best possible 82 (117.38 seconds)
Cbc0010I After 91 nodes, 46 on tree, -1e+50 best solution, best possible 82 (118.54 seconds)
Cbc0010I After 95 nodes, 48 on tree, -1e+50 best solution, best possible 82 (119.36 seconds)
Cbc0010I After 96 nodes, 48 on tree, -1e+50 best solution, best possible 82 (120.75 seconds)
Cbc0010I After 99 nodes, 51 on tree, -1e+50 best solution, best possible 82 (121.49 seconds)
Cbc0010I After 103 nodes, 54 on tree, -1e+50 best solution, best possible 82 (122.53 seconds)
Cbc0010I After 105 nodes, 55 on tree, -1e+50 best solution, best possible 82 (123.32 seconds)
Cbc0010I After 106 nodes, 55 on tree, -1e+50 best solution, best possible 82 (124.37 seconds)
Cbc0010I After 109 nodes, 57 on tree, -1e+50 best solution, best possible 82 (125.22 seconds)
Cbc0010I After 111 nodes, 57 on tree, -1e+50 best solution, best possible 82 (125.97 seconds)
Cbc0010I After 117 nodes, 60 on tree, -1e+50 best solution, best possible 82 (128.56 seconds)
Cbc0010I After 121 nodes, 63 on tree, -1e+50 best solution, best possible 82 (129.84 seconds)
Cbc0010I After 124 nodes, 65 on tree, -1e+50 best solution, best possible 82 (130.71 seconds)
Cbc0010I After 127 nodes, 67 on tree, -1e+50 best solution, best possible 82 (131.65 seconds)
Cbc0010I After 131 nodes, 70 on tree, -1e+50 best solution, best possible 82 (132.67 seconds)
Cbc0010I After 133 nodes, 71 on tree, -1e+50 best solution, best possible 82 (133.94 seconds)
Cbc0010I After 137 nodes, 74 on tree, -1e+50 best solution, best possible 82 (135.18 seconds)
Cbc0010I After 141 nodes, 77 on tree, -1e+50 best solution, best possible 82 (136.41 seconds)
Cbc0010I After 144 nodes, 79 on tree, -1e+50 best solution, best possible 82 (137.40 seconds)
Cbc0010I After 147 nodes, 81 on tree, -1e+50 best solution, best possible 82 (138.74 seconds)
Cbc0010I After 150 nodes, 82 on tree, -1e+50 best solution, best possible 82 (139.48 seconds)
Cbc0010I After 154 nodes, 85 on tree, -1e+50 best solution, best possible 82 (140.29 seconds)
Cbc0010I After 155 nodes, 86 on tree, -1e+50 best solution, best possible 82 (141.63 seconds)
Cbc0010I After 158 nodes, 89 on tree, -1e+50 best solution, best possible 82 (142.73 seconds)
Cbc0010I After 162 nodes, 92 on tree, -1e+50 best solution, best possible 82 (143.61 seconds)
Cbc0010I After 165 nodes, 93 on tree, -1e+50 best solution, best possible 82 (144.75 seconds)
Cbc0010I After 166 nodes, 93 on tree, -1e+50 best solution, best possible 82 (145.84 seconds)
Cbc0010I After 171 nodes, 96 on tree, -1e+50 best solution, best possible 82 (147.21 seconds)
Cbc0010I After 176 nodes, 100 on tree, -1e+50 best solution, best possible 82 (147.94 seconds)
Cbc0010I After 178 nodes, 100 on tree, -1e+50 best solution, best possible 82 (148.88 seconds)
Cbc0010I After 181 nodes, 102 on tree, -1e+50 best solution, best possible 82 (149.71 seconds)
Cbc0010I After 184 nodes, 103 on tree, -1e+50 best solution, best possible 82 (150.64 seconds)
Cbc0010I After 186 nodes, 103 on tree, -1e+50 best solution, best possible 82 (151.97 seconds)
Cbc0010I After 189 nodes, 106 on tree, -1e+50 best solution, best possible 82 (153.21 seconds)
Cbc0010I After 193 nodes, 109 on tree, -1e+50 best solution, best possible 82 (153.93 seconds)
Cbc0010I After 197 nodes, 111 on tree, -1e+50 best solution, best possible 82 (155.04 seconds)
Cbc0010I After 201 nodes, 114 on tree, -1e+50 best solution, best possible 82 (155.77 seconds)
Cbc0010I After 203 nodes, 116 on tree, -1e+50 best solution, best possible 82 (156.67 seconds)
Cbc0010I After 206 nodes, 117 on tree, -1e+50 best solution, best possible 82 (157.57 seconds)
Cbc0010I After 210 nodes, 120 on tree, -1e+50 best solution, best possible 82 (158.43 seconds)
Cbc0010I After 212 nodes, 120 on tree, -1e+50 best solution, best possible 82 (159.16 seconds)
Cbc0010I After 213 nodes, 121 on tree, -1e+50 best solution, best possible 82 (160.09 seconds)
Cbc0010I After 219 nodes, 121 on tree, -1e+50 best solution, best possible 82 (160.95 seconds)
Cbc0010I After 222 nodes, 123 on tree, -1e+50 best solution, best possible 82 (161.79 seconds)
Cbc0010I After 226 nodes, 126 on tree, -1e+50 best solution, best possible 82 (162.70 seconds)
Cbc0010I After 228 nodes, 128 on tree, -1e+50 best solution, best possible 82 (163.93 seconds)
Cbc0010I After 230 nodes, 129 on tree, -1e+50 best solution, best possible 82 (164.64 seconds)
Cbc0010I After 235 nodes, 133 on tree, -1e+50 best solution, best possible 82 (165.80 seconds)
Cbc0010I After 240 nodes, 135 on tree, -1e+50 best solution, best possible 82 (167.14 seconds)
Cbc0010I After 243 nodes, 136 on tree, -1e+50 best solution, best possible 82 (167.86 seconds)
Cbc0010I After 249 nodes, 138 on tree, -1e+50 best solution, best possible 82 (169.33 seconds)
Cbc0010I After 252 nodes, 140 on tree, -1e+50 best solution, best possible 82 (170.70 seconds)
Cbc0010I After 256 nodes, 142 on tree, -1e+50 best solution, best possible 82 (171.46 seconds)
Cbc0010I After 260 nodes, 144 on tree, -1e+50 best solution, best possible 82 (172.35 seconds)
Cbc0010I After 262 nodes, 146 on tree, -1e+50 best solution, best possible 82 (173.10 seconds)
Cbc0010I After 266 nodes, 147 on tree, -1e+50 best solution, best possible 82 (174.95 seconds)
Cbc0010I After 267 nodes, 147 on tree, -1e+50 best solution, best possible 82 (175.96 seconds)
Cbc0010I After 268 nodes, 148 on tree, -1e+50 best solution, best possible 82 (178.15 seconds)
Cbc0010I After 271 nodes, 151 on tree, -1e+50 best solution, best possible 82 (182.61 seconds)
Cbc0010I After 272 nodes, 151 on tree, -1e+50 best solution, best possible 82 (185.63 seconds)
Cbc0010I After 274 nodes, 153 on tree, -1e+50 best solution, best possible 82 (186.60 seconds)
Cbc0010I After 276 nodes, 153 on tree, -1e+50 best solution, best possible 82 (188.24 seconds)
Cbc0010I After 277 nodes, 154 on tree, -1e+50 best solution, best possible 82 (189.25 seconds)
Cbc0010I After 279 nodes, 156 on tree, -1e+50 best solution, best possible 82 (189.96 seconds)
Cbc0010I After 280 nodes, 156 on tree, -1e+50 best solution, best possible 82 (192.65 seconds)
Cbc0010I After 282 nodes, 157 on tree, -1e+50 best solution, best possible 82 (194.53 seconds)
Cbc0010I After 284 nodes, 159 on tree, -1e+50 best solution, best possible 82 (195.81 seconds)
Cbc0010I After 286 nodes, 161 on tree, -1e+50 best solution, best possible 82 (199.00 seconds)
Cbc0010I After 287 nodes, 162 on tree, -1e+50 best solution, best possible 82 (200.20 seconds)
Cbc0010I After 289 nodes, 163 on tree, -1e+50 best solution, best possible 82 (203.35 seconds)
Cbc0010I After 290 nodes, 163 on tree, -1e+50 best solution, best possible 82 (206.61 seconds)
Cbc0010I After 292 nodes, 164 on tree, -1e+50 best solution, best possible 82 (207.79 seconds)
Cbc0010I After 293 nodes, 165 on tree, -1e+50 best solution, best possible 82 (208.85 seconds)
Cbc0010I After 296 nodes, 166 on tree, -1e+50 best solution, best possible 82 (209.88 seconds)
Cbc0010I After 297 nodes, 166 on tree, -1e+50 best solution, best possible 82 (211.23 seconds)
Cbc0010I After 298 nodes, 166 on tree, -1e+50 best solution, best possible 82 (212.30 seconds)
Cbc0010I After 300 nodes, 168 on tree, -1e+50 best solution, best possible 82 (213.20 seconds)
Cbc0010I After 301 nodes, 169 on tree, -1e+50 best solution, best possible 82 (215.88 seconds)
Cbc0010I After 303 nodes, 171 on tree, -1e+50 best solution, best possible 82 (217.89 seconds)
Cbc0010I After 305 nodes, 173 on tree, -1e+50 best solution, best possible 82 (219.32 seconds)
Cbc0010I After 307 nodes, 174 on tree, -1e+50 best solution, best possible 82 (222.02 seconds)
Cbc0010I After 308 nodes, 174 on tree, -1e+50 best solution, best possible 82 (223.21 seconds)
Cbc0010I After 311 nodes, 176 on tree, -1e+50 best solution, best possible 82 (224.88 seconds)
Cbc0010I After 312 nodes, 176 on tree, -1e+50 best solution, best possible 82 (225.86 seconds)
Cbc0010I After 314 nodes, 176 on tree, -1e+50 best solution, best possible 82 (226.59 seconds)
Cbc0004I Integer solution of 82 found after 138086 iterations and 316 nodes (227.37 seconds)
Cbc0012I Integer solution of 82 found by heuristic after 138085 iterations and 315 nodes (227.37 seconds)
Cbc0030I Thread 0 used 27 times,  waiting to start 10.592023,  144 locks, 0.056802034 locked, 0.020977974 waiting for locks
Cbc0030I Thread 1 used 29 times,  waiting to start 7.7063391,  153 locks, 0.034678221 locked, 0.00045466423 waiting for locks
Cbc0030I Thread 2 used 25 times,  waiting to start 14.162436,  133 locks, 0.032462835 locked, 0.0035967827 waiting for locks
Cbc0030I Thread 3 used 26 times,  waiting to start 23.025289,  136 locks, 0.045662403 locked, 0.0012564659 waiting for locks
Cbc0030I Thread 4 used 22 times,  waiting to start 26.228491,  114 locks, 0.024322748 locked, 0.0039508343 waiting for locks
Cbc0030I Thread 5 used 22 times,  waiting to start 16.022214,  116 locks, 0.031300783 locked, 8.058548e-05 waiting for locks
Cbc0030I Thread 6 used 24 times,  waiting to start 32.863837,  125 locks, 0.092724323 locked, 0.00086402893 waiting for locks
Cbc0030I Thread 7 used 27 times,  waiting to start 35.985136,  137 locks, 0.041479349 locked, 0.00010442734 waiting for locks
Cbc0030I Thread 8 used 23 times,  waiting to start 34.618367,  119 locks, 0.033367634 locked, 0.0040180683 waiting for locks
Cbc0030I Thread 9 used 22 times,  waiting to start 30.993783,  112 locks, 0.027930021 locked, 0.0015795231 waiting for locks
Cbc0030I Thread 10 used 18 times,  waiting to start 30.511933,  93 locks, 0.020904779 locked, 0.00078582764 waiting for locks
Cbc0030I Thread 11 used 23 times,  waiting to start 34.515154,  119 locks, 0.030738831 locked, 0.0013768673 waiting for locks
Cbc0030I Thread 12 used 23 times,  waiting to start 31.716897,  117 locks, 0.032862902 locked, 8.893013e-05 waiting for locks
Cbc0030I Thread 13 used 18 times,  waiting to start 38.871023,  92 locks, 0.023838043 locked, 6.5803528e-05 waiting for locks
Cbc0030I Main thread 163.56113 waiting for threads,  809 locks, 0.23262024 locked, 0.020417452 waiting for locks
Cbc0001I Search completed - best objective 82, took 144753 iterations and 329 nodes (243.16 seconds)
Cbc0032I Strong branching done 6468 times (635061 iterations), fathomed 0 nodes and fixed 2 variables
Cbc0035I Maximum depth 29, 0 variables fixed on reduced cost
Cuts at root node changed objective from 82 to 82
Probing was tried 141 times and created 752 cuts (1.08733 seconds)
Gomory was tried 141 times and created 1016 cuts (2.80744 seconds)
Knapsack was tried 105 times and created 0 cuts (0.252274 seconds)
Clique was tried 105 times and created 0 cuts (0.0928044 seconds)
OddWheel was tried 105 times and created 0 cuts (0.901251 seconds)
MixedIntegerRounding2 was tried 105 times and created 0 cuts (0.101727 seconds)
FlowCover was tried 105 times and created 0 cuts (0.00676491 seconds)
TwoMirCuts was tried 105 times and created 2895 cuts (0.88564 seconds)
ZeroHalf was tried 105 times and created 330 cuts (3.34814 seconds)

Result - Optimal solution found
Objective value:                82
Enumerated nodes:               329
Total iterations:               144753
Time (CPU seconds):             1590.88
Time (Wallclock seconds):       245.763
Total time (CPU seconds):       1590.88




<OptimizationStatus.OPTIMAL: 0>

結果を表示する 

if model.status == OptimizationStatus.OPTIMAL:

    model.write("data/pism/ikegami-2shift-data1.lp")

    print("希望シフト")
    display(
        desired_shifts_temp.fillna(" ").rename(columns={"nurse_no": "staff"})
    )

    shift_satisfaction["val"] = shift_satisfaction["var"].astype(float)

    satisfaction = (
        shift_satisfaction.groupby(["staff", "day"])[["val"]]
        .sum()
        .reset_index()
    )
    not_satisfied = satisfaction[satisfaction["val"] < 1][["staff", "day"]]

    print("希望か叶わなかったシフト")
    display(
        pd.merge(
            not_satisfied,
            desired_shifts,
            left_on=["staff", "day"],
            right_on=["staff", "day"],
        )
    )

    print(
        "シフトの実現度: {}% ({}/{})".format(
            int(100 * model.objective_value / desired_shifts.shape[0]),
            int(model.objective_value),
            desired_shifts.shape[0],
        )
    )

    # 最適化された勤務表を取り出す
    shift_table["val"] = shift_table["var"].astype(float)
    optimized_shift_table = shift_table[shift_table["val"] > 0.5]
    optimized_shift_table = optimized_shift_table.pivot(
        columns="day", index="staff", values="shift"
    ).reset_index()

    # 希望が叶わなかったシフトをハイライトする
    def highlight_not_satisfied(row, target, props):
        staff = row.name
        style = [
            (
                props
                if any(t[0] == staff and t[1] == day for t in target)
                else ""
            )
            for day, shift in row.items()
        ]
        return style

    styled = optimized_shift_table.style.apply(
        highlight_not_satisfied,
        target=not_satisfied.to_numpy(),
        props="background-color: #ffffa0; color: #c00000; font-weight: bold;",
        axis=1,
    )

    # Styler未適用のDataFrameと表示を合わせる
    def add_div(func):
        def wrapper(*args, **kwargs):
            return '<div class="dataframe">{}</div>'.format(
                func(*args, **kwargs)
            )

        return wrapper

    styled.to_html = add_div(styled.to_html)

    display(styled)
else:
    print(model.status)

希望シフト
staff 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 30
0 0 n / + /
1 1 / * / / *
2 2 / / /
3 3 / /
4 4
5 5 /
6 6 + / /
7 7 + + / / x
8 8 + N n / /
9 9 /
10 10 / + / /
11 11 /
12 12 / N n /
13 13 / /
14 14 n / + / N n /
15 15 N n / N n /
16 16 +
17 17 +
18 18 / + +
19 19 /
20 20
21 21 /
22 22 + / x
23 23 +
24 24 n / N n / N n /
25 25 / N n + -
26 26 n / / + / +
27 27 
希望か叶わなかったシフト
staff day shift
0 25 6 + 
シフトの実現度: 98% (82/83)
day staff 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
0 0 n / - - N n / - / / / + N n / - - N n / - - N n / - / - + N
1 1 / N n / - - N n / / - - - N n / - + N n / - - N n / / / - -
2 2 N n / / / N n / - - N n / - - - / - - N n / - - N n / - - /
3 3 / - / / - N n / - - - N n / / - N n / / - N n / - / - - N n
4 4 - / - - / - - / - N n / - - - N n / / - N n / / - - - N n /
5 5 / / / N n / - - N n / - / - N n / - - N n / - - - N n / / -
6 6 - - N n / - + N n / / - N n / / / - / - - / + - - - N n / /
7 7 - N n / - + - + / N n / - - + / / - N n / - / / / N n / - +
8 8 N n / / - - - / - - / + - N n / / / - - - N n / / - N n / -
9 9 / - - - N n / - / / - N n / - - - + - - / - N n / / - N n /
10 10 / - - N n / - - N n / + - - / N n / - - N n / / + / / - N n
11 11 N n / / - + / / - N n / / - N n / - - N n / - - - N n / - -
12 12 / N n / / - N n / / - N n / - - N n / + - - - N n / - - - /
13 13 - - N n / - - N n / - - - / / - - N n / - / / / N n / - + -
14 14 n / / / / - + N n / / - N n / - - - + / - - N n / - - N n /
15 15 - / / - - N n / - - N n / - - N n / - - / N n / - - N n / /
16 16 - / - - N n / - - - / + / N n / / / N n / - - - N n / - / N
17 17 / - - N n / - - N n / - / + N n / / - / N n / / - - - + / N
18 18 - / / - - - + - / - - + - + - / - - / - - / - - / - - / - -
19 19 - - N n / / - - / / - N n / - - N n / - - N n / - - N n / /
20 20 / - - - / - - N n / - - - N n / - / N n / - / / - - - N n /
21 21 - N n / - - N n / - - - N n / / / N n / / - - - / + / - + N
22 22 N n / - - - + / N n / - / - N n / - - N n / / / N n / + - -
23 23 / / / N n / + - - - / - - / - - N n / / - - N n / - - - N n
24 24 n / - / - N n / / / N n / - - N n / - - - / - N n / - + - /
25 25 / / / - N n / - - N n / + - - / - - - / N n / - - N n / - -
26 26 n / - / + - N n / / - + - / / - - N n / + - / N n / / - N n
27 27 - - N n / + - / - - N n / - / / / - - - / / - - + - / - - -

データの出典