I have encountered several optimization problems that involve identifying one or more indices in a vector that maximizes or minimizes a cost. Is there a way to identify such indices in linear programming? I'm open to solutions in mathprog, CVXR, CVXPY, or any other API.

For example, identifying an index is needed for change point problems (find the index at which the function changes), putting distance constraints on the traveling salesman problem (visit city X before cumulative distance Y).

As a simple example, suppose we want to identify the location in a vector where the sum on either side is the most equal (their difference is smallest). In this example, the solution is index 5:

x = c(1, 3, 6, 4, 7, 9, 6, 2, 3)

Attempt 1

Using CVXR, I tried declaring split_index and using that as an index (e.g., x[1:split]):

library(CVXR)
split_index = Variable(1, integer = TRUE)
objective = Minimize(abs(sum(x[1:split_index]) - sum(x[(split_index+1):length(x)])))
result = solve(objective)

It errs 1:split_index with NA/NaN argument.

Attempt 2

Declare an explicit index-vector (indices) and do an elementwise logical test whether split_index <= indices. Then element-wise-multiply that binary vector with x to select one or the other side of the split:

indices = seq_along(x)
split_index = Variable(1, integer = TRUE)
is_first = split_index <= indices
objective = Minimize(abs(sum(x * is_first) - sum(x * !is_first)))
result = solve(objective)

It errs in x * is_first with non-numeric argument to binary operator. I suspect that this error arises because is_first is now an IneqConstraint object.

Symbols in red are decision variables and symbols in blue are constants.

R code:

> library(Rglpk)
> library(CVXR)
> 
> x <- c(1, 3, 6, 4, 7, 9, 6, 2, 3)
> n <- length(x)
> delta <- Variable(n, boolean=T)
> y <- Variable(2)
> order <- list()
> for (i in 2:n) {
+     order[[as.character(i)]] <- delta[i-1] <= delta[i]
+ }
> 
> 
> problem <- Problem(Minimize(abs(y[1]-y[2])),
+                    c(order,
+                      y[1] == t(1-delta) %*% x,
+                      y[2] == t(delta) %*%x))
> result <- solve(problem,solver = "GLPK", verbose=T)
GLPK Simplex Optimizer, v4.47
30 rows, 12 columns, 60 non-zeros
      0: obj =  0.000000000e+000  infeas = 4.100e+001 (2)
*     7: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
*     8: obj =  0.000000000e+000  infeas = 0.000e+000 (0)
OPTIMAL SOLUTION FOUND
GLPK Integer Optimizer, v4.47
30 rows, 12 columns, 60 non-zeros
9 integer variables, none of which are binary
Integer optimization begins...
+     8: mip =     not found yet >=              -inf        (1; 0)
+     9: >>>>>  1.000000000e+000 >=  0.000000000e+000 100.0% (2; 0)
+     9: mip =  1.000000000e+000 >=     tree is empty   0.0% (0; 3)
INTEGER OPTIMAL SOLUTION FOUND
> result$getValue(delta)
      [,1]
 [1,]    0
 [2,]    0
 [3,]    0
 [4,]    0
 [5,]    0
 [6,]    1
 [7,]    1
 [8,]    1
 [9,]    1
> result$getValue(y)
     [,1]
[1,]   21
[2,]   20
> 

The absolute value is automatically linearized by CVXR.