Visit given city before a given cumulative distance in the traveling salesman problem

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.

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)