from pyomo.environ import *
import networkx as nx

# Initialize a new directed graph for the road network (using networkx for clarity)
G = nx.Graph()

# Add edges and weights to the graph
edges = [
    (1, 2, 720), (1, 6, 347), (1, 37, 875), (2, 3, 81), (2, 9, 652),
    (3, 7, 138), (4, 5, 52), (4, 8, 139), (5, 8, 239), (5, 6, 96),
    (6, 21, 220), (7, 8, 176), (7, 9, 284), (8, 11, 313), (9, 10, 135),
    (10, 11, 84), (10, 14, 409), (11, 12, 195), (12, 13, 201), (12, 19, 353),
    (13, 14, 67), (13, 15, 260), (14, 15, 329), (15, 16, 146), (16, 17, 53),
    (15, 30, 135), (16, 18, 119), (17, 29, 243), (17, 18, 17), (18, 19, 231),
    (18, 26, 48), (19, 20, 79), (20, 21, 37), (20, 23, 138), (21, 22, 210),
    (22, 23, 55), (22, 24, 57), (23, 25, 59), (24, 36, 149), (24, 25, 49),
    (25, 26, 99), (25, 34, 129), (26, 27, 161), (27, 28, 157), (27, 33, 112),
    (28, 29, 265), (28, 31, 187), (29, 30, 46), (31, 32, 88), (32, 33, 211),
    (32, 40, 68), (33, 34, 16), (34, 35, 79), (35, 36, 38), (35, 39, 116),
    (36, 37, 52), (37, 38, 155), (38, 39, 75), (38, 41, 252), (39, 40, 163),
    (40, 41, 263)
]

# Add edges and weights to the graph
for u, v, weight in edges:
    G.add_edge(u, v, weight=weight)

# Initialize Pyomo model
model = ConcreteModel()

# Define sets for edges and nodes
model.V = Set(initialize=list(G.nodes()))
model.E = Set(initialize=[(i, j) for i, j in G.edges()] + [(j, i) for i, j in G.edges()])

# Decision variables x_ij where x_ij = 1 if edge (i,j) is traversed, otherwise 0
model.x = Var(model.E, domain=Binary)

# Connectivity variable z_ij where z_ij = 1 if edge is part of the main path
model.z = Var(model.E, domain=Binary)

# Additional variables u_i for MTZ subtour elimination
model.u = Var(model.V, domain=NonNegativeReals)

# Edge weights (distances)
weights = {(i, j): G[i][j]['weight'] for i, j in G.edges()}
for i, j in G.edges():
    weights[j, i] = weights[i, j]  # Symmetric since it's an undirected graph

# Objective: Maximize the total distance covered (maximize the sum of weights of chosen edges)
def objective_rule(model):
    return sum(weights[i, j] * model.x[i, j] for (i, j) in model.E)
model.objective = Objective(rule=objective_rule, sense=maximize)

# Constraint: No road is traveled in both directions
def direction_constraint(model, i, j):
    return model.x[i, j] + model.x[j, i] <= 1
model.direction_constraint = Constraint(model.E, rule=direction_constraint)

# Flow conservation constraint: The number of edges entering and leaving any node should be even (Eulerian path condition)
def flow_conservation_rule(model, v):
    return sum(model.x[i, v] for i in model.V if (i, v) in model.E) == sum(model.x[v, j] for j in model.V if (v, j) in model.E)
model.flow_conservation = Constraint(model.V, rule=flow_conservation_rule)

# Modified subtour elimination constraint to allow loops connected to the main path
def modified_subtour_elimination_rule(model, i, j):
    if i != 1 and j != 1:  # Skip node 1 (start node) to avoid trivial constraints
        # Allow loops only if connected to the main path
        return model.u[i] - model.u[j] + len(G.nodes()) * model.x[i, j] <= len(G.nodes()) - 1 + len(G.nodes()) * (1 - model.z[i, j])
    else:
        return Constraint.Skip
model.subtour_elimination = Constraint(model.E, rule=modified_subtour_elimination_rule)

# Ensure that all loops are connected to the main path via at least one edge
def loop_connection_rule(model, i, j):
    return model.z[i, j] <= sum(model.z[k, l] for k, l in model.E if (i == k or i == l or j == k or j == l)) 
model.loop_connection = Constraint(model.E, rule=loop_connection_rule)

# Ensure that every cycle is connected to the main path
def main_path_constraint(model, i, j):
    return model.z[i, j] >= model.x[i, j]  # z ensures that selected edge (i,j) is part of the main path if it's traversed
model.main_path_constraint = Constraint(model.E, rule=main_path_constraint)

# Fix u[1] to 0 to define the starting point
model.u[1].fix(0)

# Solve the model using CBC solver
solver = SolverFactory('cbc')  # Change from 'glpk' to 'cbc'
solver.solve(model)

# Print results
print("Longest closed path distance:", value(model.objective))
edges_traveled = [(i, j) for (i, j) in model.E if value(model.x[i, j]) > 0.5]
print("Edges in the path:", edges_traveled)


Modify the subtour elimination to include the cycles connected to my main graph