I'm creating a version of the TSP and have followed an article about it to create a version that very nearly matches my desired parameters. I want to add a line like the one currently commented out towards the bottom after the line def capacity_constraint(model,k):, such that it adds to the capacity constraint.
Right now, it checks model.y (2d boolean indicator matrix of node visitation, where i is the node and k is the vehicle) and multiplies it by model.q (1D array of node demand values). I'd like to do the same for my travel arcs, multiplying model.x (3d boolean matrix of i = node departed, j = node arrived at, k = vehicle) by distances (2d matrix of arc costs).
I think the solution is very simple and I just don't understand how the syntax I'm using is wrong, somewhere. I would greatly appreciate any help, and thank you for reading this far!
import time
from itertools import cycle
import numpy as np
from scipy.spatial.distance import pdist, squareform
import matplotlib.pyplot as plt
import matplotlib as mpl
import networkx as nx
import pyomo.environ as pyo
from pyomo.contrib.appsi.solvers.highs import Highs
np.random.seed(88) # sets seed for consistent results; set to 88
N = 13 #default = 20, scenario=6
demands = np.random.randint(8, 9, size=N) #demands is an analog for both distance and speed of travel;
#...default=1,6, scenario=6
#***NOTE: allow user to upload data that indicates the demand load for each ATON
#***NOTE: stratify this between (1) demand associated with travel and (2) demand associated with ATON service
demands[0] = 0 #starting node has zero demand obviously
capacity = 130 #model input; default= 10, scenario=3
#***NOTE: possibly stratify unless distance and ATON service demand are translated into hours
n_vehicles = 8 #model input; default = 6, scenario=3-5
# coordinates = np.random.rand(N, 2) #randomizes 2d coordinates for N nodes
#***NOTE: replace this with the coordinates of vessel Willow servicing the PR area from a specified homeport location
coordinates = [[0.25,0.55],
[0.51,0.19],
[0.58,0.16],
[0.55,0.05],
[0.62,0.19],
[0.64,0.16],
[0.63,0.06],
[0.71,0.19],
[0.73,0.16],
[0.745,0.175],
[0.815,0.155],
[0.838,0.181],
[0.893,0.04]]
coordinates = np.array(coordinates)
#distances = squareform(pdist(coordinates, metric="euclidean")) #generates a list of distances between n des. 2d matrix where i,j is distance from i to j. doesn't have to be symmetic (simulate wind for example)
#distances = np.round(distances, decimals=4) # avoid numerical errors later
distances = [[0,40,40,44,44,44,48,44,44,44,48,48,52],
[40,0,2,4,4,4,8,4,4,4,8,8,12],
[40,2,0,4,4,4,8,4,4,4,8,8,12],
[44,4,4,0,8,8,4,8,8,8,12,12,8],
[44,4,4,8,0,2,12,2,2,2,4,4,8],
[44,4,4,8,2,0,12,2,2,2,4,4,8],
[48,8,8,4,12,12,0,12,12,12,8,8,4],
[44,4,4,8,2,2,12,0,12,12,8,8,4],
[44,4,4,8,2,2,12,12,0,12,8,8,4],
[44,4,4,8,2,2,12,12,12,0,8,8,4],
[48,8,8,12,4,4,8,8,8,8,0,2,4],
[48,8,8,12,4,4,8,8,8,8,2,0,4],
[52,12,12,8,8,8,4,4,4,4,4,4,0]]
distances = np.array(distances)
model = pyo.ConcreteModel()
#variable setup
model.V = pyo.Set(initialize=range(len(demands))) #nodes
model.A = pyo.Set(initialize=[(i, j) for i in model.V for j in model.V if i != j]) #arcs
model.K = pyo.Set(initialize=range(n_vehicles)) #vehicles
model.Q = pyo.Param(initialize=capacity) #capacity ie travel distance constraint
model.q = pyo.Param(model.V, initialize={i: d for (i, d) in enumerate(demands)}) #demand load
model.c = pyo.Param(model.A, initialize={(i, j): distances[i, j] for (i, j) in model.A}) #arc costs
model.x = pyo.Var(model.A, model.K, within=pyo.Binary) #boolean indicator for arc fulfillment
model.y = pyo.Var(model.V, model.K, within=pyo.Binary) #boolean indicator for node fulfillment
#constraints
#rule that each node must be visited
def arcs_in(model, i):
if i == model.V.first():
return sum(model.x[:, i, :]) == len(model.K)
else:
return sum(model.x[:, i, :]) == 1.0
#each node must also be departed
def arcs_out(model, i):
if i == model.V.first():
return sum(model.x[i, :, :]) == len(model.K)
else:
return sum(model.x[i, :, :]) == 1.0
#if a vehicle visits a node it takes its demand onto itself
def vehicle_assignment(model, i, k):
return sum(model.x[:, i, k]) == model.y[i, k]
#every vehicle also has to leave a node
def comp_vehicle_assignment(model, i, k):
return sum(model.x[i, :, k]) == model.y[i, k]
#total demand must not exceed capacity
def capacity_constraint(model,k):
return sum(model.y[i, k] * model.q[i] for i in model.V) <= model.Q
# return sum(model.y[i, k] * model.q[i] for i in model.V) + sum(model.x[i,j,k]*distances[i,j] for (i,j) in model.A)<= model.Q
#subtour constraint. method eliminate_subtours adds to this when it runs
def subtour_elimination(model, S, Sout, h, k):
nodes_out = sum(model.x[i, j, k] for i in S for j in Sout)
return model.y[h, k] <= nodes_out
model.arcs_in = pyo.Constraint(model.V, rule=arcs_in)
model.arcs_out = pyo.Constraint(model.V, rule=arcs_out)
model.vehicle_assignment = pyo.Constraint(model.V, model.K, rule=vehicle_assignment)
model.comp_vehicle_assignment = pyo.Constraint(model.V, model.K, rule=comp_vehicle_assignment)
model.capacity_constraint = pyo.Constraint(model.K, rule=capacity_constraint)
model.subtour_elimination = pyo.ConstraintList()
I've tried a wide variety of alterations to the syntax of both the relevant constraint and the line that instantiates it in the model:
model.capacity_constraint = pyo.Constraint(model.K, rule=capacity_constraint)