SCIP/PySCIPOpt で日能研の広告問題を解きなおす（非線形制約をそのまま書く）

pip3 install pyscipopt

from pyscipopt import Model, quicksum

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


# モデル初期化（モデル名の指定は任意）
m = Model("Timetable")
if False:
  # ログ非表示（任意）
  m.hideOutput()

# 目的関数はダミー（可行解が欲しいだけ）
# SCIP は目的関数は線形のみ可、定数でもよい
m.setObjective(0)

# 0-1 のバイナリ変数（型指定は省略形 vtype="B" も可）
# 曜日 d の p 時限目に教科 s を割り当てるなら 1 とする
# PuLP のような辞書変数作成の特殊関数は使わない
x = {(d, p, s): m.addVar(vtype="Binary", name=f"x[{d},{p},{s}]")
      for d in days for p in periods for s in subjects}

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


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

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

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

# 4. 水曜日に国語と社会がある
for s in ["国語", "社会"]:
  m.addCons(quicksum(x[2, p, s] for p in periods) >= 1)

# 5. 理科を置くなら 4 限のみ（1〜3 限の理科は禁止）
for d in days:
  for p in [1, 2, 3]:
    m.addCons(x[d, p, "理科"] == 0)

# 6. 国語と社会は、それぞれ 1 日 2 時間以上あったとき、その教科は次の日にはない
# （SCIP の長所による改善）
# 「その日に 2 コマ以上 ⇒ 翌日はなし」の非線形な含意表現を直接表現：
#    a = 日 d の s のコマ数, b = 日 d+1 の s のコマ数
#    (a - 1) * b <= 0 で a >= 2 なら b == 0 を強制 (a, b は非負)
for s in ["国語", "社会"]:
  for d in range(4):
    a = quicksum(x[d, p, s] for p in periods)
    b = quicksum(x[d + 1, p, s] for p in periods)
    m.addCons((a - 1) * b <= 0)  # 非線形制約をそのまま書く

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


# 求解
m.optimize()

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

solutions = []

# 解けなくなるまで繰り返す
while True:
  m.optimize()
  if m.getStatus() != "optimal":
    break

  bs = m.getBestSol()
  sol = {
    (d, p): next(s for s in subjects if bs[x[(d, p, s)]])
      for d in days for p in periods
  }
  solutions.append(sol)

  # 再最適化したいときは、既存の解決をクリアする必要がある
  m.freeTransform()
  # 同じ割当てを禁止する制約を追加
  m.addCons(
    quicksum(x[d, p, s] for (d, p), s in sol.items())
      <= len(days) * len(periods) - 1
  )

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

m.writeParams(
  "params_default.txt",
  onlychanged=False,
)

m.writeParams(
  "params_default.txt",
  onlychanged=True,
)