facility-location-problem

施設配置問題（FLP）

あなたは施設配置問題と戦略的な立地選定の最適化に精通しています。顧客の需要と容量制約を満たしながら総コストを最小化する、施設（倉庫、工場、配送センター）の最適な立地を決定できるよう支援してください。

初期評価

施設配置問題を解く前に、次の点を把握してください。

1. 問題の種類

   • 無容量制約施設配置問題（UFLP）ですか？（容量制限なし）
   • 容量制約付き施設配置問題（CFLP）ですか？（施設に容量制限あり）
   • p-メディアン問題ですか？（施設をちょうど p 個配置）
   • p-センター問題ですか？（最大距離を最小化）
   • 固定費用立地問題ですか？（施設の開設固定費＋変動費）

2. 施設の特性

   • 施設の候補地はいくつありますか？
   • 施設ごとの固定開設費はいくらですか？
   • 運用コスト（容量に応じたコスト）はいくらですか？
   • 施設の容量はいくらですか？（容量制約がある場合）
   • 施設は複数の顧客にサービスを提供できますか？

3. 顧客の要件

   • 顧客／需要地点はいくつありますか？
   • 顧客の需要量はいくらですか？
   • サービス要件（対応距離／時間）はありますか？
   • 単一調達と複数調達のどちらが許可されていますか？

4. コスト構造

   • 施設の開設にかかる固定費はいくらですか？
   • 輸送／配送コストはいくらですか？
   • 出荷単位あたりの運用コストはいくらですか？
   • 規模の経済はありますか？

5. 制約

   • 施設をちょうど p 個開設する必要がありますか？（p-メディアン）
   • 予算制約はありますか？
   • 容量制約はありますか？
   • 対応範囲の要件はありますか？
   • 施設数の最小値／最大値はありますか？

---

問題の分類

1. 無容量制約施設配置問題（UFLP）

説明：

• 開設する施設を決定します（容量制限なし）
• 顧客を開設済みの施設に割り当てます
• 固定費と輸送費の合計を最小化します

特徴：

• 各施設は無制限の需要に対応できます
• FLP の最も基本的な変種です
• NP 困難ですが、中規模の問題は解くことができます

用途：

• 初期ネットワークの設計
• 戦略的な長期計画
• 大局的な立地選定

2. 容量制約付き施設配置問題（CFLP）

説明：

• 施設に容量制約があります
• UFLP より現実に即しています
• 各顧客は複数の施設からサービスを受けられます

特徴：

• 容量制約により問題がより難しくなります
• UFLP より多くの施設が必要になる場合があります
• より複雑な解法が必要です

用途：

• 倉庫ネットワークの設計
• 生産施設の立地選定
• サービスセンターの配置

3. p-メディアン問題

説明：

• 施設をちょうど p 個配置します
• 顧客から施設までの（重み付き）距離の合計を最小化します
• 固定費も容量制約もありません

特徴：

• 施設数は事前に決まっています
• 輸送距離の最小化に重点を置きます
• 立地科学における古典的な問題です

用途：

• 緊急サービスの立地選定
• 小売店舗の配置
• 学校／病院の立地選定

4. p-センター問題

説明：

• 施設を p 個配置します
• すべての顧客から最寄りの施設までの距離の最大値を最小化します
• ミニマックス目的（公平性を重視）

特徴：

• 顧客がサービスから遠く離れすぎないようにします
• 目的関数は p-メディアン問題とは異なります
• 緊急サービスに適しています

用途：

• 救急車の拠点配置
• 消防署の配置
• 緊急対応計画

---

数理モデル

無容量制約施設配置問題（UFLP）

集合：

• I = {1, ..., m}：施設候補地の集合
• J = {1, ..., n}：顧客の集合

パラメーター：

• f_i：立地 i に施設を開設する固定費
• c_{ij}：施設 i から顧客 j にサービスを提供するコスト（輸送費）
• d_j：顧客 j の需要量

決定変数：

• y_i ∈ {0,1}：施設 i を開設する場合は 1、それ以外は 0
• x_{ij} ∈ [0,1]：顧客 j の需要のうち、施設 i が提供する割合

目的関数：

Minimize: Σ_i f_i * y_i + Σ_i Σ_j c_{ij} * d_j * x_{ij}
          \_____________/   \___________________________/
          Fixed costs       Transportation costs

制約条件：

1. Demand satisfaction: Each customer fully served
   Σ_i x_{ij} = 1,  ∀j ∈ J

2. Facility opening: Can only serve from open facilities
   x_{ij} ≤ y_i,  ∀i ∈ I, ∀j ∈ J

3. Binary facility decisions:
   y_i ∈ {0,1},  ∀i ∈ I

4. Assignment variables:
   x_{ij} ≥ 0,  ∀i ∈ I, ∀j ∈ J

計算複雑性： NP 困難

容量制約付き施設配置問題（CFLP）

追加パラメーター：

• Q_i：施設 i の容量

変更後の制約：

1. Demand satisfaction:
   Σ_i x_{ij} = 1,  ∀j ∈ J

2. Capacity constraints:
   Σ_j d_j * x_{ij} ≤ Q_i * y_i,  ∀i ∈ I

3. Opening constraints:
   x_{ij} ≤ y_i,  ∀i ∈ I, ∀j ∈ J

4. Binary variables:
   y_i ∈ {0,1},  ∀i ∈ I
   x_{ij} ≥ 0,  ∀i ∈ I, ∀j ∈ J

p-メディアン問題

目的関数：

Minimize: Σ_i Σ_j c_{ij} * x_{ij}

制約：

1. Each customer assigned to exactly one facility:
   Σ_i x_{ij} = 1,  ∀j ∈ J

2. Exactly p facilities opened:
   Σ_i y_i = p

3. Assignment only to open facilities:
   x_{ij} ≤ y_i,  ∀i ∈ I, ∀j ∈ J

4. Binary variables:
   y_i ∈ {0,1},  ∀i ∈ I
   x_{ij} ∈ {0,1},  ∀i ∈ I, ∀j ∈ J

---

厳密解法

1. PuLP による MIP 定式化（UFLP）

from pulp import *
import numpy as np

def solve_uflp(fixed_costs, transport_costs, demands):
    """
    Solve Uncapacitated Facility Location Problem

    Args:
        fixed_costs: list of fixed costs for each facility
        transport_costs: 2D array [facilities x customers] of unit transport costs
        demands: list of customer demands

    Returns:
        dict with optimal solution
    """
    m = len(fixed_costs)  # Number of potential facilities
    n = len(demands)       # Number of customers

    # Create problem
    prob = LpProblem("UFLP", LpMinimize)

    # Decision variables
    # y[i] = 1 if facility i is opened
    y = LpVariable.dicts("facility", range(m), cat='Binary')

    # x[i,j] = fraction of customer j's demand served by facility i
    x = LpVariable.dicts("assign",
                         [(i, j) for i in range(m) for j in range(n)],
                         lowBound=0, upBound=1, cat='Continuous')

    # Objective: Minimize total cost (fixed + transportation)
    prob += (
        lpSum([fixed_costs[i] * y[i] for i in range(m)]) +
        lpSum([transport_costs[i][j] * demands[j] * x[i,j]
               for i in range(m) for j in range(n)]),
        "Total_Cost"
    )

    # Constraints

    # 1. Each customer must be fully served
    for j in range(n):
        prob += (
            lpSum([x[i,j] for i in range(m)]) == 1,
            f"Demand_Customer_{j}"
        )

    # 2. Can only serve from open facilities
    for i in range(m):
        for j in range(n):
            prob += (
                x[i,j] <= y[i],
                f"Open_Facility_{i}_{j}"
            )

    # Solve
    import time
    start_time = time.time()
    prob.solve(PULP_CBC_CMD(msg=1, timeLimit=300))
    solve_time = time.time() - start_time

    # Extract solution
    if LpStatus[prob.status] in ['Optimal', 'Feasible']:
        open_facilities = [i for i in range(m) if y[i].varValue > 0.5]

        assignments = {}
        for j in range(n):
            assignments[j] = []
            for i in range(m):
                if x[i,j].varValue > 0.01:  # Threshold for numerical issues
                    assignments[j].append((i, x[i,j].varValue))

        # Calculate cost breakdown
        fixed_cost_total = sum(fixed_costs[i] for i in open_facilities)
        transport_cost_total = sum(
            transport_costs[i][j] * demands[j] * x[i,j].varValue
            for i in range(m) for j in range(n)
        )

        return {
            'status': LpStatus[prob.status],
            'total_cost': value(prob.objective),
            'fixed_cost': fixed_cost_total,
            'transport_cost': transport_cost_total,
            'open_facilities': open_facilities,
            'num_facilities': len(open_facilities),
            'assignments': assignments,
            'solve_time': solve_time
        }
    else:
        return {
            'status': LpStatus[prob.status],
            'total_cost': None,
            'open_facilities': [],
            'solve_time': solve_time
        }


# Example usage
if __name__ == "__main__":
    # Problem data: 5 potential facilities, 10 customers
    np.random.seed(42)

    # Fixed costs to open each facility
    fixed_costs = [5000, 4500, 6000, 5500, 4800]

    # Transportation costs (facility x customer)
    # Lower cost = closer proximity
    transport_costs = np.array([
        [10, 15, 8, 20, 12, 18, 22, 14, 16, 11],
        [18, 12, 16, 14, 10, 15, 20, 25, 13, 19],
        [14, 20, 18, 10, 16, 12, 15, 18, 22, 14],
        [22, 18, 14, 16, 20, 10, 12, 16, 14, 18],
        [16, 14, 20, 18, 14, 16, 10, 12, 15, 20]
    ])

    # Customer demands
    demands = [100, 150, 80, 120, 90, 110, 130, 95, 105, 125]

    result = solve_uflp(fixed_costs, transport_costs, demands)

    print(f"\n{'='*70}")
    print(f"UNCAPACITATED FACILITY LOCATION PROBLEM - SOLUTION")
    print(f"{'='*70}")
    print(f"Status: {result['status']}")
    print(f"Total Cost: ${result['total_cost']:,.2f}")
    print(f"  Fixed Costs: ${result['fixed_cost']:,.2f}")
    print(f"  Transport Costs: ${result['transport_cost']:,.2f}")
    print(f"\nFacilities Opened: {result['num_facilities']}")
    print(f"Facility IDs: {result['open_facilities']}")
    print(f"\nSolve Time: {result['solve_time']:.2f} seconds")

    print(f"\nCustomer Assignments:")
    for customer_id, assignment in result['assignments'].items():
        print(f"  Customer {customer_id} (demand={demands[customer_id]}):")
        for facility_id, fraction in assignment:
            print(f"    → Facility {facility_id}: {fraction*100:.1f}%")

2. 容量制約付き施設配置問題（CFLP）

def solve_cflp(fixed_costs, transport_costs, demands, capacities):
    """
    Solve Capacitated Facility Location Problem

    Args:
        fixed_costs: list of fixed costs for each facility
        transport_costs: 2D array [facilities x customers]
        demands: list of customer demands
        capacities: list of facility capacities

    Returns:
        dict with optimal solution
    """
    m = len(fixed_costs)
    n = len(demands)

    # Create problem
    prob = LpProblem("CFLP", LpMinimize)

    # Decision variables
    y = LpVariable.dicts("facility", range(m), cat='Binary')
    x = LpVariable.dicts("assign",
                         [(i, j) for i in range(m) for j in range(n)],
                         lowBound=0, cat='Continuous')

    # Objective
    prob += (
        lpSum([fixed_costs[i] * y[i] for i in range(m)]) +
        lpSum([transport_costs[i][j] * demands[j] * x[i,j]
               for i in range(m) for j in range(n)]),
        "Total_Cost"
    )

    # Constraints

    # 1. Each customer fully served
    for j in range(n):
        prob += (
            lpSum([x[i,j] for i in range(m)]) == 1,
            f"Demand_{j}"
        )

    # 2. Capacity constraints at each facility
    for i in range(m):
        prob += (
            lpSum([demands[j] * x[i,j] for j in range(n)]) <= capacities[i] * y[i],
            f"Capacity_{i}"
        )

    # 3. Can only serve from open facilities
    for i in range(m):
        for j in range(n):
            prob += (
                x[i,j] <= y[i],
                f"Open_{i}_{j}"
            )

    # Solve
    start_time = time.time()
    prob.solve(PULP_CBC_CMD(msg=1, timeLimit=300))
    solve_time = time.time() - start_time

    # Extract solution
    if LpStatus[prob.status] in ['Optimal', 'Feasible']:
        open_facilities = [i for i in range(m) if y[i].varValue > 0.5]

        # Calculate utilization for each open facility
        utilization = {}
        for i in open_facilities:
            used_capacity = sum(demands[j] * x[i,j].varValue for j in range(n))
            utilization[i] = (used_capacity / capacities[i]) * 100

        assignments = {}
        for j in range(n):
            assignments[j] = []
            for i in range(m):
                if x[i,j].varValue > 0.01:
                    assignments[j].append((i, x[i,j].varValue))

        return {
            'status': LpStatus[prob.status],
            'total_cost': value(prob.objective),
            'open_facilities': open_facilities,
            'num_facilities': len(open_facilities),
            'utilization': utilization,
            'assignments': assignments,
            'solve_time': solve_time
        }
    else:
        return {
            'status': LpStatus[prob.status],
            'solve_time': solve_time
        }


# Example usage
fixed_costs = [8000, 7500, 9000, 8500, 7800]
transport_costs = np.array([
    [10, 15, 8, 20, 12, 18, 22, 14, 16, 11],
    [18, 12, 16, 14, 10, 15, 20, 25, 13, 19],
    [14, 20, 18, 10, 16, 12, 15, 18, 22, 14],
    [22, 18, 14, 16, 20, 10, 12, 16, 14, 18],
    [16, 14, 20, 18, 14, 16, 10, 12, 15, 20]
])
demands = [100, 150, 80, 120, 90, 110, 130, 95, 105, 125]
capacities = [400, 350, 450, 380, 420]  # Facility capacities

result = solve_cflp(fixed_costs, transport_costs, demands, capacities)

print(f"\n{'='*70}")
print(f"CAPACITATED FACILITY LOCATION PROBLEM - SOLUTION")
print(f"{'='*70}")
print(f"Status: {result['status']}")
print(f"Total Cost: ${result['total_cost']:,.2f}")
print(f"Facilities Opened: {result['num_facilities']}")

print(f"\nFacility Utilization:")
for facility_id in result['open_facilities']:
    print(f"  Facility {facility_id}: {result['utilization'][facility_id]:.1f}% "
          f"(capacity={capacities[facility_id]})")

3. p-メディアン問題

def solve_p_median(distances, demands, p):
    """
    Solve p-Median Problem

    Locate exactly p facilities to minimize total weighted distance

    Args:
        distances: 2D array of distances [facilities x customers]
        demands: customer demands (weights)
        p: number of facilities to open

    Returns:
        optimal solution
    """
    m = len(distances)     # Potential facility locations
    n = len(demands)       # Customers

    prob = LpProblem("p_Median", LpMinimize)

    # Decision variables
    y = LpVariable.dicts("facility", range(m), cat='Binary')
    x = LpVariable.dicts("assign",
                         [(i, j) for i in range(m) for j in range(n)],
                         cat='Binary')

    # Objective: Minimize total weighted distance
    prob += (
        lpSum([distances[i][j] * demands[j] * x[i,j]
               for i in range(m) for j in range(n)]),
        "Total_Weighted_Distance"
    )

    # Constraints

    # 1. Each customer assigned to exactly one facility
    for j in range(n):
        prob += (
            lpSum([x[i,j] for i in range(m)]) == 1,
            f"Assign_{j}"
        )

    # 2. Exactly p facilities opened
    prob += (
        lpSum([y[i] for i in range(m)]) == p,
        "p_Facilities"
    )

    # 3. Assignment only to open facilities
    for i in range(m):
        for j in range(n):
            prob += (
                x[i,j] <= y[i],
                f"Open_{i}_{j}"
            )

    # Solve
    start_time = time.time()
    prob.solve(PULP_CBC_CMD(msg=0))
    solve_time = time.time() - start_time

    if LpStatus[prob.status] in ['Optimal', 'Feasible']:
        open_facilities = [i for i in range(m) if y[i].varValue > 0.5]

        assignments = {}
        for j in range(n):
            for i in range(m):
                if x[i,j].varValue > 0.5:
                    assignments[j] = i
                    break

        return {
            'status': LpStatus[prob.status],
            'total_distance': value(prob.objective),
            'open_facilities': open_facilities,
            'assignments': assignments,
            'solve_time': solve_time
        }
    else:
        return {'status': LpStatus[prob.status]}


# Example: Locate 3 facilities among 8 candidates to serve 12 customers
np.random.seed(42)

# Generate random coordinates
facility_coords = np.random.rand(8, 2) * 100
customer_coords = np.random.rand(12, 2) * 100

# Calculate Euclidean distance matrix
distances = np.zeros((8, 12))
for i in range(8):
    for j in range(12):
        distances[i][j] = np.linalg.norm(facility_coords[i] - customer_coords[j])

demands = [100, 150, 80, 120, 90, 110, 130, 95, 105, 125, 115, 140]
p = 3  # Open exactly 3 facilities

result = solve_p_median(distances, demands, p)

print(f"\n{'='*70}")
print(f"p-MEDIAN PROBLEM - SOLUTION")
print(f"{'='*70}")
print(f"Status: {result['status']}")
print(f"Total Weighted Distance: {result['total_distance']:,.2f}")
print(f"Facilities Opened (p={p}): {result['open_facilities']}")
print(f"\nCustomer Assignments:")
for customer_id, facility_id in result['assignments'].items():
    dist = distances[facility_id][customer_id]
    print(f"  Customer {customer_id} → Facility {facility_id} "
          f"(distance={dist:.2f}, demand={demands[customer_id]})")

---

貪欲ヒューリスティクス

1. 貪欲追加アルゴリズム（p-メディアン）

def greedy_add_p_median(distances, demands, p):
    """
    Greedy heuristic for p-Median Problem

    Iteratively add facility that gives maximum cost reduction

    Time complexity: O(p * m * n)
    Quality: Typically within 10-20% of optimal

    Args:
        distances: distance matrix [facilities x customers]
        demands: customer demands
        p: number of facilities to open

    Returns:
        heuristic solution
    """
    m, n = distances.shape

    open_facilities = []
    assignments = {}

    # Initialize: assign all customers to closest potential facility
    for j in range(n):
        assignments[j] = np.argmin(distances[:, j])

    # Greedily add p facilities
    for _ in range(p):
        best_facility = None
        best_improvement = 0

        for i in range(m):
            if i in open_facilities:
                continue

            # Calculate improvement if we add facility i
            improvement = 0
            for j in range(n):
                current_dist = distances[assignments[j]][j]
                new_dist = distances[i][j]

                if new_dist < current_dist:
                    improvement += demands[j] * (current_dist - new_dist)

            if improvement > best_improvement:
                best_improvement = improvement
                best_facility = i

        if best_facility is not None:
            open_facilities.append(best_facility)

            # Update assignments
            for j in range(n):
                current_dist = distances[assignments[j]][j]
                new_dist = distances[best_facility][j]

                if new_dist < current_dist:
                    assignments[j] = best_facility

    # Calculate total cost
    total_cost = sum(
        distances[assignments[j]][j] * demands[j]
        for j in range(n)
    )

    return {
        'open_facilities': open_facilities,
        'assignments': assignments,
        'total_distance': total_cost,
        'method': 'Greedy Add'
    }

2. 貪欲削除アルゴリズム

def greedy_drop_p_median(distances, demands, p):
    """
    Greedy drop heuristic for p-Median

    Start with all facilities, iteratively close facility
    with least impact until p remain

    Args:
        distances: distance matrix
        demands: customer demands
        p: number of facilities to keep

    Returns:
        heuristic solution
    """
    m, n = distances.shape

    # Start with all facilities open
    open_facilities = set(range(m))

    # Close facilities until only p remain
    while len(open_facilities) > p:
        best_facility_to_close = None
        min_cost_increase = float('inf')

        for i in open_facilities:
            # Test closing facility i
            test_facilities = open_facilities - {i}

            # Calculate cost with this facility closed
            cost = 0
            for j in range(n):
                min_dist = min(distances[f][j] for f in test_facilities)
                cost += min_dist * demands[j]

            if cost < min_cost_increase:
                min_cost_increase = cost
                best_facility_to_close = i

        open_facilities.remove(best_facility_to_close)

    # Final assignments
    assignments = {}
    for j in range(n):
        best_facility = min(open_facilities,
                          key=lambda f: distances[f][j])
        assignments[j] = best_facility

    total_cost = sum(
        distances[assignments[j]][j] * demands[j]
        for j in range(n)
    )

    return {
        'open_facilities': list(open_facilities),
        'assignments': assignments,
        'total_distance': total_cost,
        'method': 'Greedy Drop'
    }

3. 節約量に基づくヒューリスティクス（UFLP）

def savings_heuristic_uflp(fixed_costs, transport_costs, demands):
    """
    Savings-based heuristic for UFLP

    Compute savings from opening each facility and open
    facilities in order of savings

    Args:
        fixed_costs: fixed facility costs
        transport_costs: transportation costs
        demands: customer demands

    Returns:
        heuristic solution
    """
    m = len(fixed_costs)
    n = len(demands)

    # Calculate savings for each facility
    # Savings = reduction in transport cost - fixed cost
    savings = []

    for i in range(m):
        # Calculate transport cost if facility i serves customers
        transport_saving = 0

        # Assume without facility i, customers served by most expensive option
        for j in range(n):
            # Cost to serve customer j from facility i
            cost_from_i = transport_costs[i][j] * demands[j]

            # Cost to serve from next best facility
            other_costs = [transport_costs[k][j] * demands[j]
                          for k in range(m) if k != i]

            if other_costs:
                min_other_cost = min(other_costs)
                transport_saving += max(0, min_other_cost - cost_from_i)

        net_saving = transport_saving - fixed_costs[i]
        savings.append((net_saving, i))

    # Sort by savings (descending)
    savings.sort(reverse=True)

    # Open facilities with positive savings
    open_facilities = []
    for saving, facility_id in savings:
        if saving > 0:
            open_facilities.append(facility_id)

    # If no positive savings, open at least one facility
    if not open_facilities:
        # Open facility with lowest total cost for all customers
        best_facility = None
        min_cost = float('inf')

        for i in range(m):
            cost = fixed_costs[i] + sum(
                transport_costs[i][j] * demands[j] for j in range(n)
            )
            if cost < min_cost:
                min_cost = cost
                best_facility = i

        open_facilities = [best_facility]

    # Assign customers to nearest open facility
    assignments = {}
    for j in range(n):
        best_facility = min(open_facilities,
                          key=lambda i: transport_costs[i][j])
        assignments[j] = best_facility

    # Calculate total cost
    fixed_cost_total = sum(fixed_costs[i] for i in open_facilities)
    transport_cost_total = sum(
        transport_costs[assignments[j]][j] * demands[j]
        for j in range(n)
    )
    total_cost = fixed_cost_total + transport_cost_total

    return {
        'open_facilities': open_facilities,
        'assignments': assignments,
        'total_cost': total_cost,
        'fixed_cost': fixed_cost_total,
        'transport_cost': transport_cost_total,
        'method': 'Savings Heuristic'
    }

---

局所探索による改善

1. 施設交換近傍

def facility_swap_local_search(distances, demands, initial_solution, p):
    """
    Local search with facility swap moves for p-Median

    Swap an open facility with a closed one if it improves objective

    Args:
        distances: distance matrix
        demands: customer demands
        initial_solution: starting solution
        p: number of facilities

    Returns:
        improved solution
    """
    m, n = distances.shape

    open_facilities = set(initial_solution['open_facilities'])
    closed_facilities = set(range(m)) - open_facilities

    def calculate_cost(facilities):
        """Calculate total cost for given facility set"""
        cost = 0
        for j in range(n):
            min_dist = min(distances[f][j] for f in facilities)
            cost += min_dist * demands[j]
        return cost

    current_cost = calculate_cost(open_facilities)
    improved = True

    while improved:
        improved = False
        best_swap = None
        best_cost = current_cost

        # Try all swaps
        for facility_in in closed_facilities:
            for facility_out in open_facilities:
                # Test swap
                test_facilities = (open_facilities - {facility_out}) | {facility_in}
                test_cost = calculate_cost(test_facilities)

                if test_cost < best_cost:
                    best_cost = test_cost
                    best_swap = (facility_in, facility_out)
                    improved = True

        if improved:
            facility_in, facility_out = best_swap
            open_facilities = (open_facilities - {facility_out}) | {facility_in}
            closed_facilities = (closed_facilities - {facility_in}) | {facility_out}
            current_cost = best_cost

    # Final assignments
    assignments = {}
    for j in range(n):
        best_facility = min(open_facilities,
                          key=lambda f: distances[f][j])
        assignments[j] = best_facility

    return {
        'open_facilities': list(open_facilities),
        'assignments': assignments,
        'total_distance': current_cost,
        'method': 'Facility Swap Local Search'
    }

2. 顧客の再割り当て

def customer_reassignment(open_facilities, distances, demands):
    """
    Improve solution by optimally reassigning customers

    Given open facilities, assign each customer to nearest facility

    Args:
        open_facilities: list of open facility indices
        distances: distance matrix
        demands: customer demands

    Returns:
        optimal assignments for given facilities
    """
    n = distances.shape[1]

    assignments = {}
    for j in range(n):
        best_facility = min(open_facilities,
                          key=lambda i: distances[i][j])
        assignments[j] = best_facility

    total_cost = sum(
        distances[assignments[j]][j] * demands[j]
        for j in range(n)
    )

    return {
        'assignments': assignments,
        'total_cost': total_cost
    }

---

メタヒューリスティクス

1. p-メディアン問題に対するシミュレーテッドアニーリング

import random
import math

def simulated_annealing_p_median(distances, demands, p,
                                 initial_temp=1000, cooling_rate=0.95,
                                 max_iterations=5000):
    """
    Simulated Annealing for p-Median Problem

    Args:
        distances: distance matrix
        demands: customer demands
        p: number of facilities
        initial_temp: starting temperature
        cooling_rate: temperature reduction factor
        max_iterations: maximum iterations

    Returns:
        best solution found
    """
    m, n = distances.shape

    def calculate_cost(facilities):
        cost = 0
        for j in range(n):
            min_dist = min(distances[f][j] for f in facilities)
            cost += min_dist * demands[j]
        return cost

    # Initial solution: random p facilities
    current_facilities = set(random.sample(range(m), p))
    current_cost = calculate_cost(current_facilities)

    best_facilities = current_facilities.copy()
    best_cost = current_cost

    temperature = initial_temp

    for iteration in range(max_iterations):
        # Generate neighbor: swap one facility
        closed = list(set(range(m)) - current_facilities)

        if not closed:
            break

        facility_out = random.choice(list(current_facilities))
        facility_in = random.choice(closed)

        neighbor_facilities = (current_facilities - {facility_out}) | {facility_in}
        neighbor_cost = calculate_cost(neighbor_facilities)

        delta = neighbor_cost - current_cost

        # Accept or reject
        if delta < 0 or random.random() < math.exp(-delta / temperature):
            current_facilities = neighbor_facilities
            current_cost = neighbor_cost

            if current_cost < best_cost:
                best_facilities = current_facilities.copy()
                best_cost = current_cost

        # Cool down
        temperature *= cooling_rate

        if temperature < 0.1:
            break

    # Final assignments
    assignments = {}
    for j in range(n):
        best_facility = min(best_facilities,
                          key=lambda f: distances[f][j])
        assignments[j] = best_facility

    return {
        'open_facilities': list(best_facilities),
        'assignments': assignments,
        'total_distance': best_cost,
        'method': 'Simulated Annealing'
    }

2. UFLP に対する遺伝的アルゴリズム

def genetic_algorithm_uflp(fixed_costs, transport_costs, demands,
                          population_size=50, generations=200,
                          mutation_rate=0.15):
    """
    Genetic Algorithm for Uncapacitated Facility Location

    Args:
        fixed_costs: facility fixed costs
        transport_costs: transportation costs
        demands: customer demands
        population_size: GA population size
        generations: number of generations
        mutation_rate: mutation probability

    Returns:
        best solution
    """
    m = len(fixed_costs)
    n = len(demands)

    def decode_chromosome(chromosome):
        """Convert binary chromosome to open facilities"""
        return [i for i in range(m) if chromosome[i] == 1]

    def calculate_fitness(chromosome):
        """Calculate fitness (inverse of cost)"""
        open_facilities = decode_chromosome(chromosome)

        if not open_facilities:
            return 0  # Invalid solution

        # Fixed costs
        fixed_cost = sum(fixed_costs[i] for i in open_facilities)

        # Transportation costs
        transport_cost = 0
        for j in range(n):
            min_cost = min(transport_costs[i][j] for i in open_facilities)
            transport_cost += min_cost * demands[j]

        total_cost = fixed_cost + transport_cost
        return 1.0 / (1.0 + total_cost)

    def create_individual():
        """Create random chromosome"""
        # Ensure at least one facility open
        chromosome = [random.randint(0, 1) for _ in range(m)]
        if sum(chromosome) == 0:
            chromosome[random.randint(0, m-1)] = 1
        return chromosome

    def crossover(parent1, parent2):
        """Uniform crossover"""
        child = []
        for i in range(m):
            child.append(parent1[i] if random.random() < 0.5 else parent2[i])

        # Ensure at least one facility
        if sum(child) == 0:
            child[random.randint(0, m-1)] = 1

        return child

    def mutate(chromosome):
        """Bit flip mutation"""
        for i in range(m):
            if random.random() < mutation_rate:
                chromosome[i] = 1 - chromosome[i]

        # Ensure at least one facility
        if sum(chromosome) == 0:
            chromosome[random.randint(0, m-1)] = 1

        return chromosome

    # Initialize population
    population = [create_individual() for _ in range(population_size)]

    best_chromosome = None
    best_fitness = 0

    for generation in range(generations):
        # Evaluate fitness
        fitnesses = [calculate_fitness(ind) for ind in population]

        # Track best
        max_fitness = max(fitnesses)
        if max_fitness > best_fitness:
            best_fitness = max_fitness
            best_chromosome = population[fitnesses.index(max_fitness)].copy()

        # Selection and reproduction
        new_population = []

        # Elitism
        elite_count = int(0.1 * population_size)
        elite_indices = sorted(range(len(fitnesses)),
                             key=lambda i: fitnesses[i],
                             reverse=True)[:elite_count]
        new_population = [population[i].copy() for i in elite_indices]

        # Create offspring
        while len(new_population) < population_size:
            # Tournament selection
            parent1 = max(random.sample(list(zip(population, fitnesses)), 3),
                         key=lambda x: x[1])[0]
            parent2 = max(random.sample(list(zip(population, fitnesses)), 3),
                         key=lambda x: x[1])[0]

            child = crossover(parent1, parent2)
            child = mutate(child)

            new_population.append(child)

        population = new_population

    # Decode best solution
    open_facilities = decode_chromosome(best_chromosome)

    # Assign customers
    assignments = {}
    for j in range(n):
        best_facility = min(open_facilities,
                          key=lambda i: transport_costs[i][j])
        assignments[j] = best_facility

    # Calculate costs
    fixed_cost = sum(fixed_costs[i] for i in open_facilities)
    transport_cost = sum(
        transport_costs[assignments[j]][j] * demands[j]
        for j in range(n)
    )

    return {
        'open_facilities': open_facilities,
        'assignments': assignments,
        'total_cost': fixed_cost + transport_cost,
        'fixed_cost': fixed_cost,
        'transport_cost': transport_cost,
        'method': 'Genetic Algorithm'
    }

---

完全な施設配置ソルバー

import time

class FacilityLocationSolver:
    """
    Comprehensive Facility Location Problem Solver

    Supports UFLP, CFLP, p-Median, and multiple solution methods
    """

    def __init__(self):
        self.fixed_costs = None
        self.transport_costs = None
        self.demands = None
        self.capacities = None
        self.problem_type = None

    def load_problem(self, fixed_costs=None, transport_costs=None,
                    demands=None, capacities=None, problem_type='UFLP'):
        """
        Load problem data

        Args:
            fixed_costs: list of facility fixed costs
            transport_costs: 2D array [facilities x customers]
            demands: customer demands
            capacities: facility capacities (for CFLP)
            problem_type: 'UFLP', 'CFLP', or 'p-Median'
        """
        self.fixed_costs = np.array(fixed_costs) if fixed_costs else None
        self.transport_costs = np.array(transport_costs)
        self.demands = np.array(demands)
        self.capacities = np.array(capacities) if capacities else None
        self.problem_type = problem_type

        print(f"Loaded {problem_type} problem:")
        print(f"  Facilities: {len(transport_costs)}")
        print(f"  Customers: {len(demands)}")
        print(f"  Total demand: {sum(demands)}")
        if capacities:
            print(f"  Total capacity: {sum(capacities)}")

    def solve_exact(self, time_limit=300):
        """
        Solve using exact MIP method

        Args:
            time_limit: time limit in seconds

        Returns:
            optimal/best solution
        """
        print(f"\nSolving {self.problem_type} with MIP (exact)...")

        if self.problem_type == 'UFLP':
            return solve_uflp(self.fixed_costs, self.transport_costs, self.demands)
        elif self.problem_type == 'CFLP':
            return solve_cflp(self.fixed_costs, self.transport_costs,
                            self.demands, self.capacities)
        else:
            raise ValueError(f"Exact solve not implemented for {self.problem_type}")

    def solve_heuristic(self, method='greedy', p=None):
        """
        Solve using heuristic method

        Args:
            method: 'greedy_add', 'greedy_drop', 'savings', 'sa', 'ga'
            p: number of facilities (for p-median)

        Returns:
            heuristic solution
        """
        print(f"\nSolving with {method} heuristic...")

        if method == 'greedy_add' and self.problem_type == 'p-Median':
            return greedy_add_p_median(self.transport_costs, self.demands, p)

        elif method == 'greedy_drop' and self.problem_type == 'p-Median':
            return greedy_drop_p_median(self.transport_costs, self.demands, p)

        elif method == 'savings' and self.problem_type == 'UFLP':
            return savings_heuristic_uflp(self.fixed_costs,
                                         self.transport_costs,
                                         self.demands)

        elif method == 'sa':
            if self.problem_type == 'p-Median':
                return simulated_annealing_p_median(
                    self.transport_costs, self.demands, p)
            else:
                raise ValueError("SA only for p-Median currently")

        elif method == 'ga' and self.problem_type == 'UFLP':
            return genetic_algorithm_uflp(self.fixed_costs,
                                         self.transport_costs,
                                         self.demands)

        else:
            raise ValueError(f"Unknown method {method} for {self.problem_type}")

    def compare_methods(self, methods, p=None):
        """
        Compare multiple solution methods

        Args:
            methods: list of method names
            p: number of facilities (for p-median)

        Returns:
            comparison dataframe
        """
        import pandas as pd

        results = []

        for method in methods:
            start_time = time.time()

            try:
                if method == 'exact':
                    solution = self.solve_exact()
                    cost_key = 'total_cost' if 'total_cost' in solution else 'total_distance'
                else:
                    solution = self.solve_heuristic(method, p)
                    cost_key = 'total_cost' if 'total_cost' in solution else 'total_distance'

                solve_time = time.time() - start_time

                results.append({
                    'Method': method,
                    'Cost': solution[cost_key],
                    'Facilities': len(solution['open_facilities']),
                    'Status': solution.get('status', 'Heuristic'),
                    'Time (s)': f"{solve_time:.3f}"
                })

            except Exception as e:
                print(f"  Error with {method}: {e}")

        df = pd.DataFrame(results)

        # Calculate gap from best
        if len(df) > 0:
            best_cost = df['Cost'].min()
            df['Gap %'] = ((df['Cost'] - best_cost) / best_cost * 100).round(2)

        return df

    def visualize_solution(self, solution, facility_coords, customer_coords):
        """
        Visualize facility location solution

        Args:
            solution: solution dictionary
            facility_coords: array of facility coordinates
            customer_coords: array of customer coordinates
        """
        import matplotlib.pyplot as plt

        plt.figure(figsize=(12, 8))

        # Plot all potential facilities (gray)
        for i, coord in enumerate(facility_coords):
            if i not in solution['open_facilities']:
                plt.plot(coord[0], coord[1], 's', color='lightgray',
                        markersize=10, alpha=0.5)

        # Plot open facilities (red)
        for i in solution['open_facilities']:
            coord = facility_coords[i]
            plt.plot(coord[0], coord[1], 's', color='red',
                    markersize=15, label='Open Facility' if i == solution['open_facilities'][0] else '')

        # Plot customers (blue)
        for j, coord in enumerate(customer_coords):
            plt.plot(coord[0], coord[1], 'o', color='blue',
                    markersize=8, alpha=0.7)

        # Draw assignments
        if 'assignments' in solution:
            for customer_id, assignment in solution['assignments'].items():
                if isinstance(assignment, int):
                    facility_id = assignment
                    fraction = 1.0
                else:
                    facility_id, fraction = assignment[0] if isinstance(assignment, list) else (assignment, 1.0)

                customer_coord = customer_coords[customer_id]
                facility_coord = facility_coords[facility_id]

                plt.plot([customer_coord[0], facility_coord[0]],
                        [customer_coord[1], facility_coord[1]],
                        'k-', alpha=0.2, linewidth=0.5)

        plt.xlabel('X Coordinate')
        plt.ylabel('Y Coordinate')
        plt.title(f'Facility Location Solution - {solution.get("method", "Unknown Method")}')
        plt.legend()
        plt.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.show()


# Complete example
if __name__ == "__main__":
    print("="*70)
    print("FACILITY LOCATION PROBLEM - COMPREHENSIVE EXAMPLE")
    print("="*70)

    # Generate problem data
    np.random.seed(42)

    n_facilities = 10
    n_customers = 20

    # Generate coordinates
    facility_coords = np.random.rand(n_facilities, 2) * 100
    customer_coords = np.random.rand(n_customers, 2) * 100

    # Calculate distance matrix
    distances = np.zeros((n_facilities, n_customers))
    for i in range(n_facilities):
        for j in range(n_customers):
            distances[i][j] = np.linalg.norm(
                facility_coords[i] - customer_coords[j]
            )

    # Problem parameters
    fixed_costs = np.random.uniform(5000, 8000, n_facilities)
    transport_costs = distances * 10  # Cost per unit per distance
    demands = np.random.uniform(50, 200, n_customers)
    capacities = np.random.uniform(800, 1200, n_facilities)

    # Solve UFLP
    print("\n" + "="*70)
    print("UNCAPACITATED FACILITY LOCATION PROBLEM")
    print("="*70)

    solver = FacilityLocationSolver()
    solver.load_problem(fixed_costs, transport_costs, demands, problem_type='UFLP')

    # Compare methods
    comparison = solver.compare_methods(['exact', 'savings', 'ga'])
    print("\nMethod Comparison:")
    print(comparison.to_string(index=False))

    # Solve and visualize best solution
    best_solution = solver.solve_exact()
    print(f"\nBest Solution:")
    print(f"  Total Cost: ${best_solution['total_cost']:,.2f}")
    print(f"  Fixed Cost: ${best_solution['fixed_cost']:,.2f}")
    print(f"  Transport Cost: ${best_solution['transport_cost']:,.2f}")
    print(f"  Facilities Opened: {len(best_solution['open_facilities'])}")
    print(f"  Facility IDs: {best_solution['open_facilities']}")

    # Visualize
    solver.visualize_solution(best_solution, facility_coords, customer_coords)

---

ツールとライブラリ

Python ライブラリ

最適化：

• PuLP：MIP モデリング（UFLP／CFLP に推奨）
• Pyomo：高度な最適化モデリング
• OR-Tools (Google)：高速なヒューリスティクスと厳密解法
• Gurobi：商用ソルバー（アカデミックライセンスを利用可能）
• CPLEX：IBM の商用ソルバー

特化型：

• scikit-learn：施設配置のためのクラスタリング
• scipy.optimize：汎用最適化
• geopy：実際の座標を使った距離計算

可視化：

• matplotlib：基本的な可視化
• plotly：インタラクティブな地図
• folium：地理情報の可視化
• seaborn：統計データの可視化

商用ソフトウェア

• Gurobi Optimizer：最先端の MIP ソルバー
• IBM ILOG CPLEX：エンタープライズ向け最適化
• FICO Xpress：最適化スイート
• AMPL：モデリング言語とソルバー

クラウドサービス

• Google Maps Platform：Distance Matrix API
• AWS Location Service：位置情報ベースのサービス
• Azure Maps：Microsoft の地図サービス

---

よくある課題と解決策

課題：大規模問題

問題：

• 数千の施設候補
• 数万の顧客
• 厳密解法では時間がかかりすぎる

解決策：

• ヒューリスティクスを使用する（節約量、貪欲追加／削除）
• ラグランジュ緩和
• 分解法（Benders）
• メタヒューリスティクス（SA、GA、タブーサーチ）
• 前処理による問題規模の縮小
• 求解前に顧客をクラスタリングする

課題：容量制約

問題：

• CFLP は UFLP より難しい
• 多数の施設が必要になる場合がある
• タイトな上下界を得るのが難しい

解決策：

• 定式化を強化する有効不等式
• 分枝価格アルゴリズム
• 列生成
• 特化型ヒューリスティクス
• 許容される場合は複数施設からの調達を検討する

課題：複数の目的

問題：

• コストを最小化し、サービスを最大化する
• 施設数とサービス品質の間にトレードオフがある
• ステークホルダーによって選好が異なる

解決策：

• 多目的最適化
• 加重和法
• ε制約法
• パレートフロンティア分析
• 異なる重みによるシナリオ分析

課題：需要の変動

問題：

• 需要が時間とともに変化する
• 季節変動がある
• 施設に関する意思決定が長期にわたる

解決策：

• ロバスト最適化
• 確率計画法
• シナリオベースの計画
• 柔軟な／モジュール式の施設設計
• ローリングホライズン方式

課題：実世界の距離

問題：

• ユークリッド距離は現実的ではない
• 実際の道路網が必要
• 交通状況を考慮する必要がある

解決策：

• Google Maps Distance Matrix API を使用する
• 実際の距離を事前計算する
• 時間帯による影響を考慮する
• 実際の物流データを使用する

---

出力形式

施設配置問題の解決レポート

問題のインスタンス：

• 問題の種類：容量制約なし施設配置問題
• 施設の候補地：10
• 顧客数：20
• 総需要：2,450 単位
• 最適化：固定費と輸送費の合計を最小化

最適解：

指標	値
総コスト	$87,342
固定費	$23,500
輸送費	$63,842
開設施設数	4
施設あたりの平均需要	612.5 単位
解決時間	3.7 秒
解の状態	最適

開設する施設：

施設 ID	所在地	固定費	対応する顧客数	総需要
2	(34.2, 67.8)	$6,200	6	720 単位
5	(78.1, 23.4)	$5,800	5	590 単位
7	(12.5, 45.9)	$5,900	4	480 単位
9	(56.3, 89.2)	$5,600	5	660 単位

顧客の割り当て：

顧客	需要	割り当て先の施設	距離	輸送費
C01	120	施設 2	12.3 km	$1,476
C02	150	施設 5	8.7 km	$1,305
...	...	...	...	...

コストの内訳：

• 固定費：総コストの 26.9%
• 輸送費：総コストの 73.1%
• 単位あたりの平均コスト：$35.67

---

確認する質問

詳細が必要な場合は、次の点を確認してください。

1. 施設配置問題の種類は何ですか？（UFLP、CFLP、p-Median、p-Center）
2. 施設の候補地はいくつありますか？顧客は何人ですか？
3. 施設の開設にかかる固定費はいくらですか？
4. 施設に容量制約はありますか？
5. 輸送費／サービス費はいくらですか？
6. 座標または距離行列はありますか？
7. 需要は固定ですか、それとも変動／不確実ですか？
8. 施設数に下限または上限はありますか？
9. 単一ソーシング（顧客ごとに施設を 1 か所割り当てる方式）ですか、それとも複数ソーシングも可能ですか？
10. 計画期間はどのくらいですか？（短期の戦術計画か、長期の戦略計画か）
11. 特別な要件はありますか？（カバー距離、サービスレベル）
12. 総資本投資額に予算制約はありますか？

---

関連スキル

• hub-location-problem：ハブアンドスポークネットワーク向け
• warehouse-location-optimization：倉庫に特化した考慮事項
• distribution-center-network：流通ネットワークの設計
• network-design：より広範なサプライチェーンネットワークの設計
• set-covering-problem：カバレッジに基づく施設配置
• network-flow-optimization：フローに基づく施設割り当て
• optimization-modeling：一般的な MIP 定式化手法
• multi-objective-optimization：複数の評価基準による最適化
• stochastic-optimization：需要が不確実なシナリオ
• inventory-routing-problem：配置・在庫・配送ルーティングの統合