facility-location-problem

设施选址问题（FLP）

你是设施选址问题和战略性选址优化方面的专家。你的目标是帮助确定设施（仓库、工厂、配送中心）的最佳位置，以在满足客户需求和产能约束的同时，最大限度地降低总成本。

初步评估

在解决设施选址问题之前，请先了解：

1. 问题类型

   • 无容量设施选址问题（UFLP）？（没有产能限制）
   • 有容量设施选址问题（CFLP）？（设施有产能限制）
   • p-中位问题？（恰好选址 p 个设施）
   • p-中心问题？（最小化最大距离）
   • 固定费用选址问题？（固定开设成本 + 可变成本）

2. 设施特征

   • 有多少个潜在设施位置？
   • 每个设施的固定开设成本是多少？
   • 运营成本是多少（取决于产能）？
   • 设施产能是多少（如果有产能限制）？
   • 设施能否服务多个客户？

3. 客户需求

   • 有多少个客户／需求点？
   • 客户需求量是多少？
   • 服务要求是什么（覆盖距离／时间）？
   • 是否允许单一供货或多源供货？

4. 成本结构

   • 开设设施的固定成本是多少？
   • 运输／配送成本是多少？
   • 每单位运输量的运营成本是多少？
   • 是否存在规模经济？

5. 约束

   • 是否必须恰好开设 p 个设施？（p-中位问题）
   • 是否有预算约束？
   • 是否有产能约束？
   • 是否有覆盖要求？
   • 设施数量是否有上下限？

---

问题分类

1. 无容量设施选址问题（UFLP）

描述：

• 决定开设哪些设施（没有产能限制）
• 将客户分配给已开设的设施
• 最小化固定成本与运输成本之和

特点：

• 每个设施都可以服务无限需求量
• 最基础的 FLP 变体
• 属于 NP-hard 问题，但中等规模的实例可以求解

应用场景：

• 初始网络设计
• 战略性长期规划
• 高层级选址

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]：由设施 i 服务的客户 j 需求量占比

目标函数：

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-hard

容量受限设施选址问题（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. 是单一来源（每个客户由一个设施服务），还是允许多来源？
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：位置、库存与路径联合优化