I would like to add an additional constraint to the traveling salesman problem: that a given city is visited within a given distance (say `100`

) from start. Is there a way to do this? My question is related to this unanswered CS question.

I have a mixed integer program using the R package `CVXR`

that find the shortest route without subroutes (see below). The city order is represent in the vector `node_order`

. The strategy I’ve pursued so far is:

- Re-organize
`node_order`

so that the index is the order and the value is the city id - Look up the associated distances in
`distances`

- Compute a vector with the cumulative sum of these distances.
- Add the constraint that city
`i`

must occur before the first index in (3) exceeding the distance constraint for that city.

The issues I’ve encountered with this approach is that I have not found a way to include finding-index-by-value into the optimization in `CVXR`

. This is needed in both step (1) and (4) above. Maybe this is possible after all, or there is another approach? I am willing to use other packages than `CVXR`

and other software than `R`

.

My current program:

`library(CVXR) # Make distances N = 10 distances = matrix(1:(N*N), ncol = N) # Flag 1 iff we travel that path. 0 otherwise do_transition = Variable(N, N, boolean = TRUE) # Minimize the total duration of the traveled paths. objective = Minimize(sum(do_transition * distances)) # Only go one tour. Order is 1:(N-1) node_order = Variable(N-1, integer = TRUE) ii = t(node_order %*% matrix(1, ncol = N - 1)) # repeat as N rows jj = node_order %*% matrix(1, ncol = N - 1) # repeat as N cols # Constraints constraints = list( do_transition * diag(N) == 0, # Disallow transitions to self (diagonal elements) sum_entries(do_transition, 1) == rep(1, N), # Exactly one entrance to each node sum_entries(do_transition, 2) == rep(1, N), # Exactly one exit from each node (jj - ii) + N * do_transition[2:N, 2:N] <= N - 1, # One tour constraint (no subtours) node_order >= 1, # This interval represents order as ranks (1 to N-1) node_order <= N-1 ) # Find optimum solution = solve(Problem(objective, constraints)) `

A bit of code pertaining to my current (unsuccessful) attempts:

`# Get tour order #tour = order(c(NA, result$ getValue(node_order))) # R solution tour = rep(NA, N-1) tour[result$ getValue(node_order)] = 2:N # Get tour distances distances_optim = diag(distances[tour, tour[2:N]]) # Tour cumulative distances distances_cumul = cumsum_axis(distances_optim) `