【最適解を求めるのは諦める】整数計画問題を線形緩和してソルバーで解いて、緩和問題の解を確率と解釈して許容解を求められる？

import random
import pulp

def generate_set_cover_instance(n=10, m=20, p=0.3, cost_min=1, cost_max=5, seed=None):
    """
    【集合被覆問題のインスタンス生成】
    - n: ユニバース（要素）の数
    - m: 集合の数
    - p: 各集合が特定の要素を含む確率
    - cost_min, cost_max: 各集合のコストを生成する際の最小値・最大値
    - seed: 乱数シード（再現性が必要な場合に使用）

    戻り値: (U, S, cost)
      U: 要素を格納したリスト (例: [0,1,2,...])
      S: 各集合（部分集合）のリスト（S[j] は j番目の部分集合）
      cost: 各集合のコストを格納したリスト
    """
    if seed is not None:
        random.seed(seed)
    
    # ユニバース（要素）のリスト
    U = list(range(n))
    # 各集合を保持するリスト
    S = []
    # 各集合のコストを保持するリスト
    cost = []

    # 乱数を使って集合とコストを生成
    for j in range(m):
        subset = [i for i in U if random.random() < p]
        S.append(subset)
        cost.append(random.randint(cost_min, cost_max))
    
    return U, S, cost


def solve_lp_relaxation_and_get_xstar(U, S, cost):
    """
    【集合被覆問題の線形緩和をPuLPで解き、連続解 x_star を取得する関数】
    モデル：
      min ∑(cost[j] * x_j)
      s.t. ∑(x_j : i ∈ S_j) >= 1, すべての要素 i
           0 <= x_j <= 1 （連続変数）

    戻り値: 
      x_star: 各変数 x_j の連続解（list of float）
      obj_lp: 線形緩和での目的関数値（float）
      status: ソルバーの終了ステータス名（文字列）
    """
    n = len(U)   # 要素数
    m = len(S)   # 集合の個数
    
    # PuLPで問題を定義（最小化問題）
    prob = pulp.LpProblem("SetCoverLP", pulp.LpMinimize)
    
    # 変数 x_j を連続変数として定義 (0 <= x_j <= 1)
    x = pulp.LpVariable.dicts('x', range(m), lowBound=0, upBound=1, cat=pulp.LpContinuous)
    
    # 目的関数: Σ (コスト_j * x_j)
    prob += pulp.lpSum(cost[j] * x[j] for j in range(m)), "Objective"
    
    # 被覆制約: すべての要素 i が少なくとも1つの集合に含まれるようにする
    for i in range(n):
        prob += pulp.lpSum(x[j] for j in range(m) if i in S[j]) >= 1, f"CoverConstraint_{i}"
    
    # ソルバーで解を求める（CBCがデフォルト）
    result_status = prob.solve(pulp.PULP_CBC_CMD(msg=0))
    
    # 線形緩和の目的関数値
    obj_lp = pulp.value(prob.objective)
    # 連続解 x_star
    x_star = [pulp.value(x[j]) for j in range(m)]
    
    return x_star, obj_lp, pulp.LpStatus[result_status]


def randomized_rounding(x_star, cost, S):
    """
    【乱択丸め法】
    - 線形緩和で得た x_star[j] を採用確率とみなし、
      random.random() < x_star[j] により 0 or 1 を決定。
    - 得られた二値解における目的値と、
      「被覆できなかった要素（制約違反）の数」を返す。

    戻り値:
      obj_rounded: 丸め後の目的値
      violations: カバーできなかった要素の個数
      rounded_sol: 二値解（リスト）
    """
    m = len(S)
    # 要素数（Sが空の場合を回避するための簡易的な取得方法）
    if m == 0:
        n = 0
    else:
        # 最大の要素番号+1を要素数とみなす（厳密には要注意）
        n = max((max(s) for s in S if len(s) > 0), default=-1) + 1
    
    # 乱択丸めを実行
    rounded_sol = []
    for j in range(m):
        if random.random() < x_star[j]:
            rounded_sol.append(1)
        else:
            rounded_sol.append(0)
    
    # 丸め後の目的値を計算
    obj_rounded = sum(cost[j] * rounded_sol[j] for j in range(m))
    
    # 制約違反（カバーされなかった要素の数）をカウント
    violations = 0
    for i in range(n):
        # i を含む集合が一つも選ばれていなければカバー失敗
        covered = sum(rounded_sol[j] for j in range(m) if i in S[j])
        if covered < 1:
            violations += 1
    
    return obj_rounded, violations, rounded_sol

def experience(num_trials = 100):
    
    feasible_count = 0     # 制約違反が0（完全被覆）の回数
    total_violations = 0   # 違反要素数の合計
    total_obj = 0          # 丸め後の目的値の合計
    
    for _ in range(num_trials):
        #========================================================================
        # 1. 問題インスタンスを生成（集合被覆問題）
        #========================================================================
        #   n=10個の要素、m=20個の集合、p=0.3で「その要素が集合に含まれる確率」を設定。
        #   seed=0 でランダムを固定化して再現性を確保。
        #========================================================================
        U, S, cost = generate_set_cover_instance(n=10, m=20, p=0.3, seed=None)
        n = len(U)
        m = len(S)
    
        #========================================================================
        # 2. 線形緩和を解き、連続変数 x_star を取得
        #========================================================================
        x_star, obj_lp, status = solve_lp_relaxation_and_get_xstar(U, S, cost)
    
        #========================================================================
        # 3. 乱択丸め法を複数回(例:100回)実行して結果を集計
        #========================================================================
        _, violations, _ = randomized_rounding(x_star, cost, S)
        total_violations += violations
        
        # カバーできなかった要素が0なら「許容解」とみなす
        if violations == 0:
            feasible_count += 1
    
    # 乱択丸め法の結果を統計的にまとめる
    avg_violations = total_violations / (num_trials - feasible_count)
    
    #========================================================================
    # 4. 実験結果を表示
    #========================================================================
    print(f"=== 乱択丸め法を {num_trials} 回試行した結果 ===")
    print(f" - 許容解（制約違反0）の回数 : {feasible_count} / {num_trials}")
    print(f" - 制約違反数（被覆漏れ要素数）の平均値 : {avg_violations:.2f}")

# 実行する
experience(num_trials = 100)

import random
import pulp

def generate_shift_scheduling_instance(numP=10, numD=28, numS=3, unavail_prob=0.2, seed=None):
    """
    【シフトスケジューリング問題のインスタンス生成】
    - numP: 従業員数 (|P|)
    - numD: 日数 (|D|)
    - numS: シフト数 (|S|)
    - unavail_prob: 出勤不可の確率（1なら出勤不可、0なら出勤可能）
    - seed: 乱数シード（再現性のため）
    
    戻り値: (P, D, S, availability)
      P: 従業員のリスト（例: [0,1,2,...]）
      D: 日付のリスト（例: [0,1,...]）
      S: シフトのリスト（例: [0,1,...]）
      availability: 辞書 {(p,d): 1なら従業員pが日dに出勤不可、0なら出勤可能}
    """
    if seed is not None:
        random.seed(seed)
    
    P = list(range(numP))
    D = list(range(numD))
    S = list(range(numS))
    
    availability = {}
    for p in P:
        for d in D:
            # 1: 出勤不可, 0: 出勤可能
            availability[(p, d)] = 1 if random.random() < unavail_prob else 0
    
    return P, D, S, availability

def solve_lp_relaxation_shift_scheduling(P, D, S, availability):
    """
    【シフトスケジューリング問題の線形緩和をPuLPで解いて連続解 x_star を取得する関数】
    モデル：
      変数 x_{d,s}^p ∈ [0,1]（連続変数）
      補助変数 t : 各従業員の総出勤回数の上限（最大出勤回数を t とする）
      
      目的関数：
        min t
      制約条件：
        (1) ∀d∈D, ∀s∈S:  ∑_{p∈P} x_{d,s}^p = 1
        (2) ∀d∈D, ∀p∈P:  ∑_{s∈S} x_{d,s}^p ≤ 1
        (3) ∀p∈P, ∀d∈D, ∀s∈S:  x_{d,s}^p ≤ 1 - a_d^p
        (4) ∀p∈P:  ∑_{d∈D, s∈S} x_{d,s}^p ≤ 10
        (5) ∀p∈P:  ∑_{d∈D, s∈S} x_{d,s}^p ≤ t   （各従業員の出勤回数の上限を t で表す）
    
    戻り値:
      x_star: 辞書 {(p,d,s): 連続解値}
      t_star: LPの目的関数値（tの値）
      status: ソルバー終了ステータス（文字列）
    """
    # 問題の定義（最小化問題）
    model = pulp.LpProblem("ShiftSchedulingLP", pulp.LpMinimize)
    
    # 変数 x_{d,s}^p (連続変数, 0 <= x <= 1)
    x = {}
    for p in P:
        for d in D:
            for s in S:
                x[(p, d, s)] = pulp.LpVariable(f"x_{p}_{d}_{s}", lowBound=0, upBound=1, cat=pulp.LpContinuous)
    
    # 補助変数 t
    t = pulp.LpVariable("t", lowBound=0, cat=pulp.LpContinuous)
    
    # 目的関数： t を最小化（各従業員の出勤回数の最大値を小さくする）
    model += t, "Minimize_maximum_workload"
    
    # 制約1: 各日 d, 各シフト s に対して、必ず1人の従業員が割り当てられる
    for d in D:
        for s in S:
            model += pulp.lpSum(x[(p, d, s)] for p in P) == 1, f"Cover_day{d}_shift{s}"
    
    # 制約2: 各日 d, 各従業員 p に対して、その日に複数のシフトに入らない
    for d in D:
        for p in P:
            model += pulp.lpSum(x[(p, d, s)] for s in S) <= 1, f"OneShiftPerDay_p{p}_d{d}"
    
    # 制約3: 出勤不可日の従業員には勤務を割り当てない
    for p in P:
        for d in D:
            for s in S:
                model += x[(p, d, s)] <= 1 - availability[(p, d)], f"Availability_p{p}_d{d}_s{s}"
    
    # 制約4: 各従業員の総出勤回数が10回以内
    for p in P:
        model += pulp.lpSum(x[(p, d, s)] for d in D for s in S) <= 10, f"MaxShifts_p{p}"
    
    # 制約5: 各従業員の総出勤回数が t 以下（t = 各従業員の出勤回数の上限）
    for p in P:
        model += pulp.lpSum(x[(p, d, s)] for d in D for s in S) <= t, f"Define_t_p{p}"
    
    # LP を解く
    model.solve(pulp.PULP_CBC_CMD(msg=0))
    t_star = pulp.value(model.objective)
    
    # 連続解 x_star を辞書に格納
    x_star = {}
    for p in P:
        for d in D:
            for s in S:
                x_star[(p, d, s)] = pulp.value(x[(p, d, s)])
    
    return x_star, t_star, pulp.LpStatus[model.status]

def randomized_rounding_shift(x_star, P, D, S):
    """
    【乱択丸め法によるシフトスケジューリング問題の解生成】
    - LP の連続解 x_star を採用確率として、各変数 x_{d,s}^p を 0/1 に丸める。
    
    戻り値:
      x_round: 辞書 {(p,d,s): 0 または 1}
    """
    x_round = {}
    for p in P:
        for d in D:
            for s in S:
                x_round[(p, d, s)] = 1 if random.random() < x_star[(p, d, s)] else 0
    return x_round

def evaluate_solution(x_round, P, D, S, availability):
    """
    【乱択丸め後の解の評価】
    以下の各制約について、違反している回数を計算する。
    
    制約1: 各日 d, 各シフト s で、割り当てられている従業員数がちょうど1であること
    制約2: 各日 d, 各従業員 p で、その日のシフト割り当てが1回以内であること
    制約3: 出勤不可日の従業員に勤務が割り当てられていないこと
    制約4: 各従業員 p の総出勤回数が10回以内であること
    
    戻り値:
      viol1, viol2, viol3, viol4: 各制約の違反数
      max_workload: 各従業員の出勤回数の最大値（丸め後の目的値）
    """
    # 制約1評価
    viol1 = 0
    for d in D:
        for s in S:
            assigned = sum(x_round[(p, d, s)] for p in P)
            if assigned != 1:
                viol1 += 1

    # 制約2評価
    viol2 = 0
    for d in D:
        for p in P:
            shifts_assigned = sum(x_round[(p, d, s)] for s in S)
            if shifts_assigned > 1:
                viol2 += 1

    # 制約3評価
    viol3 = 0
    for p in P:
        for d in D:
            for s in S:
                if availability[(p, d)] == 1 and x_round[(p, d, s)] == 1:
                    viol3 += 1

    # 制約4評価
    viol4 = 0
    for p in P:
        total_shifts = sum(x_round[(p, d, s)] for d in D for s in S)
        if total_shifts > 10:
            viol4 += 1

    # 各従業員の出勤回数の最大値（目的値）
    max_workload = max(sum(x_round[(p, d, s)] for d in D for s in S) for p in P)
    
    return viol1, viol2, viol3, viol4, max_workload

def experience_shift(num_trials=100):
    """
    【シフトスケジューリング問題の乱択丸め法による試行実験】
    - インスタンスを生成し、LP 線形緩和を解いて連続解 x_star を取得
    - 乱択丸め法により 0-1 解を生成し、各制約違反数と最大出勤回数を評価
    - 複数回試行し、各違反数の平均と許容解（全制約違反0）の回数を出力
    """
    total_viol1 = 0
    total_viol2 = 0
    total_viol3 = 0
    total_viol4 = 0
    total_max_workload = 0
    feasible_count = 0  # 全制約違反がゼロの試行回数

    # 乱択丸め法の試行を複数回実施
    for _ in range(num_trials):
    
        # インスタンス生成（例: 従業員10名、日28日、シフト3種類、出勤不可確率0.2）
        P, D, S, availability = generate_shift_scheduling_instance(numP=10, numD=28, numS=3, unavail_prob=0.2, seed=None)
    
        # LP 線形緩和を解いて連続解を取得
        x_star, t_star, status = solve_lp_relaxation_shift_scheduling(P, D, S, availability)
    
    
        x_round_trial = randomized_rounding_shift(x_star, P, D, S)
        viol1, viol2, viol3, viol4, max_w = evaluate_solution(x_round_trial, P, D, S, availability)
        
        total_viol1 += viol1
        total_viol2 += viol2
        total_viol3 += viol3
        total_viol4 += viol4
        total_max_workload += max_w
        
        if viol1 == 0 and viol2 == 0 and viol3 == 0 and viol4 == 0:
            feasible_count += 1
    
    avg_viol1 = total_viol1 / num_trials
    avg_viol2 = total_viol2 / num_trials
    avg_viol3 = total_viol3 / num_trials
    avg_viol4 = total_viol4 / num_trials
    avg_max_workload = total_max_workload / num_trials
    
    print("=== 乱択丸め法を複数回試行した結果 ===")
    print(f"試行回数: {num_trials}")
    print(f"許容解（全制約違反0）の回数: {feasible_count}")
    print(f"制約1（各日・シフトに1人割り当て）違反の平均: {avg_viol1:.2f}")
    print(f"制約2（各日・従業員は1シフト以内）違反の平均: {avg_viol2:.2f}")
    print(f"制約3（出勤不可日の割り当て禁止）違反の平均: {avg_viol3:.2f}")
    print(f"制約4（各従業員の総出勤回数<=10）違反の平均: {avg_viol4:.2f}")
    print(f"最大出勤回数の平均: {avg_max_workload:.2f}")

experience_shift(num_trials=100)

import matplotlib.pyplot as plt
import matplotlib.patches as patches

def plot_gantt(x_round, P, D, S, availability):
    """
    【ガントチャートによるシフトスケジュールの表示】
    各従業員（y軸）に対して、各日（x軸）の勤務をバーで表示する。
    - バーの色はシフトに対応（ここではシフト0,1,2に対して異なる色を設定）
    - 1日に複数シフト割り当てされている場合は、バーを分割し赤枠で表示（制約2違反）
    - 出勤不可日に割り当てがあれば、バー中央に赤い "X" を表示（制約3違反）
    - さらに、1日あたりの割り当てシフト数が |S| と合わない場合（制約1違反）は、
      その日付全体を赤い点線枠で囲む。
    """
    shift_colors = {0: 'skyblue', 1: 'lightgreen', 2: 'lightcoral'}
    bar_height = 0.8  # バーの高さ
    
    fig, ax = plt.subplots(figsize=(15, 6))
    
    # -----------------------------
    # 各従業員ごとのバーを描画
    # -----------------------------
    for p in P:
        for d in D:
            # 従業員 p の日 d に割り当てられたシフトリスト
            assigned_shifts = [s for s in S if x_round.get((p, d, s), 0) == 1]
            if assigned_shifts:
                multiple_assignment = (len(assigned_shifts) > 1)
                num_assigned = len(assigned_shifts)
                for idx, s in enumerate(assigned_shifts):
                    # 1日を単位幅とし、複数割り当ての場合は幅を分割して表示
                    width = 1.0 / num_assigned
                    x_start = d + idx * width
                    rect = patches.Rectangle(
                        (x_start, p - bar_height/2),  # 左下の座標
                        width,                       # 横幅
                        bar_height,                  # 縦幅
                        facecolor=shift_colors.get(s, 'gray'),
                        edgecolor='red' if multiple_assignment else 'black',
                        linewidth=2 if multiple_assignment else 1
                    )
                    ax.add_patch(rect)
                    
                    # 出勤不可日に割り当てがあれば "X" を表示（制約3違反）
                    if availability.get((p, d), 0) == 1:
                        ax.text(x_start + width/2, p, "X", color="red",
                                ha="center", va="center", fontsize=14, fontweight="bold")
    
    # -----------------------------
    # 各日付で割り当てシフト数をチェックし、制約1違反を可視化
    # -----------------------------
    for d in D:
        # 1日あたりの合計シフト割り当て数
        day_assigned = sum(x_round.get((p, d, s), 0) for p in P for s in S)
        # もし所定のシフト数 len(S) と異なれば、日付全体を赤い点線枠で囲む
        if day_assigned != len(S):
            # y方向は 従業員IDの最小 -0.5 から 最大 +0.5 まで
            rect_day = patches.Rectangle(
                (d, min(P) - 0.5),       # 左下の座標 (x, y)
                1,                       # 横幅（1日ぶん）
                len(P),                  # 縦幅（従業員全体）
                facecolor='none',
                edgecolor='red',
                linestyle='--',
                linewidth=2
            )
            ax.add_patch(rect_day)
    
    ax.set_xlim(0, max(D)+1)
    ax.set_ylim(min(P)-0.5, max(P)+0.5)
    ax.set_xlabel("Day")
    ax.set_ylabel("Worker (ID)")
    ax.set_yticks(P)
    ax.set_title("Gantt Chart")
    plt.gca().invert_yaxis()  # 従業員IDが小さい順に上部に表示
    plt.grid(True, axis='x', linestyle='--', alpha=0.5)
    plt.show()