# Constrain traveling salesman: visit a given city within a given distance from start

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:

1. Re-organize node_order so that the index is the order and the value is the city id
2. Look up the associated distances in distances
3. Compute a vector with the cumulative sum of these distances.
4. 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.

# 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))


# Unsuccessful attempt

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)


node_order[5 - 1] <= 3

node_order[7 - 1] <= 5,