如何控制根节点Lazy Constraints数量？PCTSP割生成异常求助

Let's break down why your code isn't adding enough cuts at the root node, and walk through the key fixes needed:

Core Issues in Your Current Code #

1. Wrong Reference Node for Subtour Detection
   You're using node 1 as the depot in your cut-finding functions, but your actual depot is n+1 (defined in your main code). This means you're missing valid subtours that don't include node 1 but do include the real depot, and incorrectly targeting the wrong sets for min-cuts.

2. No Check for Optimal Node Relaxation
   Your MIPNODE callback processes the node even if the relaxation solution isn't optimal. This can lead to using incomplete/invalid solutions to generate cuts, which won't work as intended.

3. Limited Cut Addition per Callback Trigger
   For min-cuts, you only add one cut and stop, but you need to keep adding cuts until the root relaxation has no more violated subtour constraints.

Step-by-Step Fixes #

1. Correct Depot Reference in Cut Functions #

First, update your run_DisconnectedCycle and run_MinCut functions to use the actual depot node (pass it as a parameter instead of hardcoding 1):

def run_DisconnectedCycle(G, nodos, depot, capacity, peso):
    G.remove_edges_from(list(G.edges))
    edges = [(i, j) for (i, j) in capacity if capacity[i, j] > 1e-8]
    G.add_edges_from(edges)
    st_list = list(nx.simple_cycles(G))
    subtour_list = []
    for cycle in st_list:
        # Exclude cycles that include the depot
        if depot not in cycle:
            # Find the node in the cycle with the highest gamma value
            max_gamma_node = max(cycle, key=lambda node: peso[node])
            subtour_list.append((cycle, max_gamma_node))
    return subtour_list

def run_MinCut(G, nodos, depot, capacity, peso):
    G.remove_edges_from(list(G.edges))
    G.add_edges_from(capacity)
    for (i, j) in capacity:
        G[i][j]['capacity'] = capacity[i, j]
    most_violated = None
    max_violation = 0
    for j in nodos:
        if j == depot:
            continue
        # Compute min-cut from depot to node j
        cap, (S, S_c) = nx.minimum_cut(G, depot, j)
        # Check if the subtour constraint is violated: sum(x[delta_out(S_c)]) >= gamma[j]
        violation = peso[j] - cap
        if violation > 1e-8:  # Account for floating point errors
            if violation > max_violation:
                max_violation = violation
                most_violated = (j, S_c)
    return most_violated

2. Update the MIPNODE Callback #

Add a check for optimal relaxation status, and ensure you keep processing cuts until no more violated constraints exist (Gurobi will re-trigger the callback after adding cuts):

elif where == GRB.Callback.MIPNODE:
    # Only process if the node relaxation is optimal
    node_status = m.cbGet(GRB.Callback.MIPNODE_STATUS)
    if node_status != GRB.OPTIMAL:
        return
    
    depth = m.cbGet(GRB.Callback.MIPNODE_NODCNT)
    if depth == 0:  # Only process root node
        x_rec = m.cbGetNodeRel(x)
        gamma_rec = m.cbGetNodeRel(gamma)
        
        # First add disconnected cycle cuts (capacity 0)
        st_DisconnectedCycle = run_DisconnectedCycle(graph.G, nodos, depot, x_rec, gamma_rec)
        for subtour in st_DisconnectedCycle:
            cycle, gamma_node = subtour
            # Add lazy constraint: sum of edges leaving the cycle >= gamma[gamma_node]
            m.cbLazy(quicksum(x[e] for e in delta_out_S(index_2, cycle)) >= gamma[gamma_node])
            cuts.SEC_frac_DC += 1
        
        # If no disconnected cycles, add min-cut cuts until no violations
        if not st_DisconnectedCycle:
            min_cut_result = run_MinCut(graph.G, nodos, depot, x_rec, gamma_rec)
            if min_cut_result:
                j, S_c = min_cut_result
                m.cbLazy(quicksum(x[e] for e in delta_out_S(index_2, S_c)) >= gamma[j])
                cuts.SEC_frac_MC += 1

3. Fix Callback Function Calls #

Make sure you pass the depot parameter when calling the cut functions in your callback.

Additional Notes #

• Floating Point Tolerance: Always use a small epsilon (like 1e-8) when comparing floating point values (e.g., checking if a capacity is non-zero or if a constraint is violated) to avoid issues with numerical precision.
• Lazy Constraints: Gurobi's lazy constraint callback will automatically re-solve the relaxation after you add cuts, so you don't need to manually re-run the solver. The callback will trigger again with the updated solution, allowing you to add more cuts if needed.
• Subtour Constraint Logic: Your current constraint sum(x[delta_out(S)]) >= gamma[k] (where k is the highest gamma node in S) is valid for PCTSP, but you could also use the more standard sum(x[delta_out(S)]) >= sum(gamma[i] for i in S) if you want a tighter relaxation.

Modified Full Code #

Here's the complete code with all fixes applied:

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
from gurobipy import *

def ins_run(n, alpha, t=0):
    np.random.seed(t)
    depot = n + 1
    nodos = [i + 1 for i in range(depot)]
    pos = dict(zip([i for i in nodos],list(map(tuple,np.random.random([n, 2]) * 5))))
    pos[depot] = pos[1]
    c_cam = {(i, j): np.sqrt((pos[i][0] - pos[j][0]) ** 2 + (pos[i][1] - pos[j][1]) ** 2) for i in nodos for j in nodos if i != j}
    c_drone = {e: alpha * c_cam[e] for e in c_cam}
    return pos, c_cam, c_drone

def delta_out(arcos, i):
    return [e for e in arcos if e[0] == i]

def delta_in(arcos, i):
    return [e for e in arcos if e[1] == i]

def delta_out_S(arcos, S):
    return [e for i in S for e in delta_out(arcos, i) if e[1] not in S]

def run_DisconnectedCycle(G, nodos, depot, capacity, peso):
    G.remove_edges_from(list(G.edges))
    edges = [(i, j) for (i, j) in capacity if capacity[i, j] > 1e-8]
    G.add_edges_from(edges)
    st_list = list(nx.simple_cycles(G))
    subtour_list = []
    for cycle in st_list:
        if depot not in cycle:
            max_gamma_node = max(cycle, key=lambda node: peso[node])
            subtour_list.append((cycle, max_gamma_node))
    return subtour_list

def run_MinCut(G, nodos, depot, capacity, peso):
    G.remove_edges_from(list(G.edges))
    G.add_edges_from(capacity)
    for (i, j) in capacity:
        G[i][j]['capacity'] = capacity[i, j]
    most_violated = None
    max_violation = 0
    for j in nodos:
        if j == depot:
            continue
        cap, (S, S_c) = nx.minimum_cut(G, depot, j)
        violation = peso[j] - cap
        if violation > 1e-8:
            if violation > max_violation:
                max_violation = violation
                most_violated = (j, S_c)
    return most_violated

if __name__ == '__main__':
    def PCTSP_CB(m, where):
        def get_sol(model):
            x_sol = []
            for e in index_2:
                arco_x = model.cbGetSolution([x[e]])
                if arco_x[0] > 0.5:
                    x_sol.append(e)
            return x_sol

        if where == GRB.Callback.MIPSOL:
            x_sol = get_sol(m)
            G = nx.DiGraph(x_sol)
            subtour_list = list(nx.simple_cycles(G))
            # Remove cycles containing depot
            subtour_list = [st for st in subtour_list if depot not in st]
            for subtour in subtour_list:
                # Find max gamma node in subtour
                max_gamma_node = max(subtour, key=lambda node: gamma[node].X)
                m.cbLazy(quicksum(x[e] for e in delta_out_S(index_2, subtour)) >= gamma[max_gamma_node])
                cuts.SEC_int += 1

        elif where == GRB.Callback.MIPNODE:
            node_status = m.cbGet(GRB.Callback.MIPNODE_STATUS)
            if node_status != GRB.OPTIMAL:
                return
            
            depth = m.cbGet(GRB.Callback.MIPNODE_NODCNT)
            if depth == 0:
                x_rec = m.cbGetNodeRel(x)
                gamma_rec = m.cbGetNodeRel(gamma)
                
                st_DisconnectedCycle = run_DisconnectedCycle(graph.G, nodos, depot, x_rec, gamma_rec)
                for subtour in st_DisconnectedCycle:
                    cycle, gamma_node = subtour
                    m.cbLazy(quicksum(x[e] for e in delta_out_S(index_2, cycle)) >= gamma[gamma_node])
                    cuts.SEC_frac_DC += 1
                
                if not st_DisconnectedCycle:
                    min_cut_result = run_MinCut(graph.G, nodos, depot, x_rec, gamma_rec)
                    if min_cut_result:
                        j, S_c = min_cut_result
                        m.cbLazy(quicksum(x[e] for e in delta_out_S(index_2, S_c)) >= gamma[j])
                        cuts.SEC_frac_MC += 1

    class Graph():
        def __init__(self):
            self.G = nx.DiGraph()
            self.S_track = list()

    class Cuts():
        def __init__(self):
            self.SEC_frac_DC = 0
            self.SEC_frac_MC = 0
            self.SEC_int = 0

    # Data
    n = 30
    seed = 11
    alpha = 0.5
    pos, c, d = ins_run(n, alpha, seed)
    depot = n + 1
    nodos = [i + 1 for i in range(depot)]
    index_2 = [(i, j) for i in nodos for j in nodos if i != j]
    graph = Graph()
    cuts = Cuts()

    model = Model()
    model.Params.LazyConstraints = 1
    x = model.addVars(index_2, vtype=GRB.BINARY, name='x')
    gamma = model.addVars(nodos, vtype=GRB.CONTINUOUS, ub=1, name='z')

    # Flow conservation constraints
    model.addConstrs(quicksum(x[e] for e in delta_in(index_2, i)) == gamma[i] for i in nodos)
    model.addConstrs(quicksum(x[e] for e in delta_out(index_2, i)) == gamma[i] for i in nodos)

    # Depot constraints
    model.addConstr(x[depot, 1] == 1)
    model.addConstr(x[1, depot] == 0)
    model.addConstr(gamma[1] == 1)
    model.addConstr(gamma[depot] == 1)

    # Additional constraints
    model.addConstrs(x[e] <= quicksum(x[r] for r in delta_in(index_2, depot) if r[0] <= e[1]) for e in delta_out(index_2, 1))
    model.addConstr(quicksum(gamma[i] for i in nodos) >= round((len(nodos) - 1) / 2 + 0.1) + 1)

    # Objective
    obj = quicksum(x[e] * c[e] for e in index_2)
    model.setObjective(obj, GRB.MINIMIZE)

    model.optimize(PCTSP_CB)
    print(f"Disconnected cycle cuts added: {cuts.SEC_frac_DC}")
    print(f"Min-cut cuts added: {cuts.SEC_frac_MC}")
    print(f"Integer subtour cuts added: {cuts.SEC_int}")

These changes should ensure that your code correctly identifies and adds all necessary subtour cuts at the root node, leading to a tighter relaxation and better performance in solving your PCTSP instance.

内容的提问来源于stack exchange，提问作者SebaVLL