整数計画問題へ入門するために PuLP/CBC で日能研の広告問題を解く

import pulp

days = list(range(5))  # 月〜金曜日
periods = [1, 2, 3, 4]  # 1限〜4限
subjects = ["算数", "国語", "社会", "理科"]


# 目的関数はダミー（可行解が欲しいだけ）
model = pulp.LpProblem("Timetable", pulp.LpStatusOptimal)

# PuLP での変数定義、0-1 のバイナリ変数
# 曜日 d の p 時限目に教科 s を割り当てるなら 1 とする
x = pulp.LpVariable.dicts(
  "x",
  ((d, p, s) for d in days for p in periods for s in subjects),
  cat="Binary",
  lowBound=0,
  upBound=1,
)

# 各コマに必ずちょうど 1 科目を割り当てる（基本ルール）
for d in days:
  for p in periods:
    model += pulp.lpSum(x[d, p, s] for s in subjects) == 1


# 設問の制約条件を追加
# 1. 月曜日の時間割を与えられたもので固定
model += x[0, 1, "算数"] == 1
model += x[0, 2, "国語"] == 1
model += x[0, 3, "算数"] == 1
model += x[0, 4, "国語"] == 1

# 2. 算数は毎日ちょうど 2 コマ
for d in days:
  model += pulp.lpSum(x[d, p, "算数"] for p in periods) == 2

# 3. 国語、理科、社会は、それぞれ 5 日間で 3 時間以上ある
model += pulp.lpSum(x[d, p, "国語"] for p in periods for d in days) >= 3
model += pulp.lpSum(x[d, p, "理科"] for p in periods for d in days) >= 3
model += pulp.lpSum(x[d, p, "社会"] for p in periods for d in days) >= 3

# 4. 水曜日に国語と社会がある
model += pulp.lpSum(x[2, p, "国語"] for p in periods) >= 1
model += pulp.lpSum(x[2, p, "社会"] for p in periods) >= 1

# 5. 理科を置くなら 4 限のみ
model += pulp.lpSum(x[d, p, "理科"] for p in [1, 2, 3] for d in days) == 0

# 6. 国語と社会は、それぞれ 1 日 2 時間以上あったとき、その教科は次の日にはない
y = pulp.LpVariable.dicts(
  # 月〜木でそれら科目が 2 時間以上あるかを示す従属的な補助変数 y を導入
  "y",
  ((d, s) for d in range(4) for s in ["国語", "社会"]),
  cat="Binary",
  lowBound=0,
  upBound=1,
)
z = pulp.LpVariable.dicts(
  # 火〜金でそれら科目が存在するかを示す従属的な補助変数 z を導入
  "z",
  ((d, s) for d in range(4) for s in ["国語", "社会"]),
  cat="Binary",
  lowBound=0,
  upBound=1,
)
for s in ["国語", "社会"]:
  for d in [0, 1, 2, 3]:
    # 補助変数 y を決定する
    model += pulp.lpSum(x[d, p, s] for p in periods) - 1 <= 3 * y[d, s]
    model += pulp.lpSum(x[d, p, s] for p in periods) >= 2 * y[d, s]
    # 補助変数 z を決定する
    model += pulp.lpSum(x[d + 1, p, s] for p in periods) <= 4 * z[d, s]
    model += pulp.lpSum(x[d + 1, p, s] for p in periods) >= z[d, s]
    # 含意 (y ならば z ではない) を表現する
    model += y[d, s] <= 1 - z[d, s]

# 7. 同じ科目は連続禁止
for d in days:
  for p in [1, 2, 3]:
    for s in subjects:
      model += x[d, p, s] + x[d, p + 1, s] <= 1

# 求解
model.solve()

# 解の表示（火曜日の時間割）
schedule = [
  (p, next(s for s in subjects if pulp.value(x[1, p, s]))) for p in periods
]
for p, s in schedule:
  print(f"{p} 限: {s}")

solutions = []

# 解けなくなるまで繰り返す
while model.solve(pulp.PULP_CBC_CMD(msg=False)) == 1:
  sol = {
    (d, p): s for d in days for p in periods for s in subjects
      if pulp.value(x[d, p, s])
  }
  assert len(sol) == len(days) * len(periods)
  solutions.append(sol)

  # 同じ割当てを禁止する制約を追加
  model += pulp.lpSum(
    x[d, p, s] for (d, p), s in sol.items()
  ) <= len(days) * len(periods) - 1

print("可行解数 =", len(solutions))