基于OR-Tools的设施选址MIP实现问题：转译SCIP方案后出现全零解及索引错误

Let's walk through the problems in your code and fix them—both the modeling logic that's causing the full-zero solution, and the index error when extracting results.

1. Why You're Getting a Full-Zero Solution #

The core issue is incorrect variable definitions and constraints that don't properly model the facility location problem:

• Wrong x[i][j] variable type: You defined x[i][j] as a boolean variable, but in a capacitated facility location problem, x[i][j] should either be:
  • A continuous non-negative variable representing how much of customer i's demand is served by facility j, or
  • A binary variable indicating whether customer i is fully served by facility j (for single-source FLP)
    Your current boolean x paired with Sum(x[i][j]) <= demands[i] makes no sense—since demands[i] is likely larger than 1, this constraint is always satisfied, so the solver picks all zeros to minimize cost (no setup or transport costs).
• Demand constraint is too weak: Instead of allowing demand to be under-served (<=), you need to enforce full demand satisfaction with ==. Depending on your variable type:
  • For continuous x: Sum(x[i][j] for j in facility_count) == demands[i]
  • For binary x: Sum(x[i][j] for j in facility_count) == 1 (each customer has exactly one serving facility)
• Redundant constraint: The x[i][j] <= demands[i] * y[j] line is unnecessary for continuous x (the capacity constraint already links x to y). For binary x, it should simplify to x[i][j] <= y[j] (you can't assign a customer to a closed facility).

2. Fixing the Index Error #

The IndexError happens for two reasons:

• Full-zero solution: Since the solver returns all zeros, numpy.where(x_val[:, j] == 1) returns an empty array—trying to access [0] on that fails. This will go away once the model is fixed to find a valid solution.
• Incorrect loop range: You're iterating over range(len(customers)) but accessing x_val[:, j]—x_val is shaped (number_of_customers, number_of_facilities), so j should iterate over facilities if you're checking columns, but actually, you want to loop over each customer i and find which facility j serves them (checking rows of x_val).

3. Corrected Code #

Here's a revised version using continuous x[i][j] variables (the standard approach for capacitated FLP, allowing split demand if needed):

def or_tools_scip_mine(facilities, customers, time_limit=None):
    import numpy
    from ortools.linear_solver import pywraplp
    if time_limit is None:
        time_limit = 1000 * 60  # 1 minute (in milliseconds)
    solver = pywraplp.Solver.CreateSolver('SCIP')
    customer_count = range(len(customers))
    facility_count = range(len(facilities))
    
    # Define variables
    y = [solver.BoolVar(f"y({j})") for j in facility_count]  # 1 if facility j is open
    x = [[solver.NumVar(0, solver.infinity(), f"x({i},{j})") for j in facility_count] 
         for i in customer_count]  # Demand from i served by j
    
    # Extract problem parameters
    facility_capacities = [facilities[i][2] for i in facility_count]
    facility_setup_costs = [facilities[i][1] for i in facility_count]
    demands = [customers[i][1] for i in customer_count]
    c = dist_matrix(facilities, customers)  # Assume this returns a 2D cost matrix
    
    # Add constraints
    # 1. All customer demand must be satisfied
    for i in customer_count:
        solver.Add(solver.Sum(x[i][j] for j in facility_count) == demands[i], f"Demand({i})")
    
    # 2. Facility capacity can't be exceeded if the facility is open
    for j in facility_count:
        solver.Add(solver.Sum(x[i][j] for i in customer_count) <= facility_capacities[j] * y[j], 
                   f"Capacity({j})")
    
    # Set objective: minimize total setup + transportation cost
    setup_cost = solver.Sum(facility_setup_costs[j] * y[j] for j in facility_count)
    transport_cost = solver.Sum(c[i][j] * x[i][j] for i in customer_count for j in facility_count)
    solver.Minimize(setup_cost + transport_cost)
    
    # Solve the model
    solver.SetTimeLimit(time_limit)
    result_status = solver.Solve()
    
    # Print solver status for debugging
    status_map = {
        solver.OPTIMAL: "Optimal solution found",
        solver.FEASIBLE: "Feasible solution found (time limit reached)",
        solver.INFEASIBLE: "No feasible solution exists",
        solver.UNBOUNDED: "Unbounded problem"
    }
    print(f"Solver Status: {status_map.get(result_status, 'Unknown status')}")
    
    # Extract solution
    solution = []
    x_val = numpy.array([[x[i][j].solution_value() for j in facility_count] for i in customer_count])
    
    for i in customer_count:
        # Find which facility serves this customer (use epsilon to handle floating point precision)
        serving_facilities = numpy.where(x_val[i] > 1e-6)[0]
        if serving_facilities.size == 0:
            solution.append(None)  # Shouldn't happen if model is correct
        else:
            solution.append(serving_facilities[0])  # Pick first serving facility (adjust if split demand is allowed)
    
    return solver.Objective().Value(), solution

4. Additional Tips #

• Single-Source FLP: If you need each customer to be served by exactly one facility, switch x[i][j] to BoolVar, change the demand constraint to Sum(x[i][j]) == 1, and adjust the capacity constraint to Sum(demands[i] * x[i][j] for i in customer_count) <= facility_capacities[j] * y[j].
• Floating Point Precision: Always use a small epsilon (like 1e-6) when checking for non-zero continuous variables—solvers often return tiny values instead of exact zeros.
• Check Solver Status: Always verify the result_status to understand if the solver found an optimal solution, a feasible one, or if the problem is infeasible.

内容的提问来源于stack exchange，提问作者Asael Bar Ilan