pythonoptimizationor-toolseuclidean-distancecp-sat

Maximizing the individual distance between 2D points with OR tools CP-SAT


I am trying to write write an algorithm which maximizes the sum of Euclidian distance between 2D points. Based on (and thanks to) this discussion, I want to put the whole thing in for loop so that this can be done for multiple points.

Here is what I attempted to do:

from ortools.sat.python import cp_model
model = cp_model.CpModel()

distance_range = 20
number_of_points = 3

a = []    # list to collect all points
# create set of required (x, y) points 
for i in range(number_robots):
    x_i = model.NewIntVar(1, distance_range , 'x_%i' %i)
    y_i = model.NewIntVar(1, distance_range , 'y_%i' %i)
    a.append((x_i, y_i))

# find difference 
x_diff = model.NewIntVar(-distance_range , distance_range , 'x_diff')
y_diff = model.NewIntVar(-distance_range , distance_range , 'y_diff')
difference = []
for i in range(2):
    for j in range(2):
         model.AddAbsEquality(x_diff, a[i][0] - a[j][0])
         model.AddAbsEquality(y_diff, a[i][1] - a[j][1])
         difference.append((x_diff, y_diff))

# square the difference value
x_diff_sq = model.NewIntVar(0, distance_range^2 , '')
y_diff_sq = model.NewIntVar(0, distance_range^2 , '')
difference_sq = []

for i in range(2):
    model.AddMultiplicationEquality(x_diff_sq, a_diff[i][0], a_diff[i][0])
    model.AddMultiplicationEquality(y_diff_sq, a_diff[i][1], a_diff[i][1])
    difference_sq.append((x_diff_sq, y_diff_sq))

# sum of square distances, (x_i^2 + y_i^2)
distance_sq= model.NewIntVar(0, 2*distance_range , '')
for i in range(2):
    model.AddAbsEquality(distance_sq, difference_sq[i][0]+ difference_sq[i][1])

#Objective
model.Maximize(distance_sq)

# Creates a solver and solves the model.
solver = cp_model.CpSolver()
status = solver.Solve(model)


if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
    for i in range(number_of_points):
        print("x_%i" %i, solver.Value(a[i][0]))
        print("y_%i" %i, solver.Value(a[i][1]))
else:
    print('No solution found.')

I tried executing the above given code but I never get the distance maximized. All the points have the same coordinate values. So I tried putting the below given constraint for 3 points

model.AllDifferent([a[0][0], a[1][0], a[2][0])    # all values of x should be different
model.AllDifferent([a[0][1], a[1][1], a[2][1])    # all values of y should be different

Then I am not getting any optimal solution. How can I tackle this?

I was able to use the above code to maximize/minimize the distance between a fixed point and multiple variables (decision variables) but few of the points are overlapping and that is why I am trying to enforce this maximization objective between the points.


Solution

  • Here is a working version with the square distance:

    from ortools.sat.python import cp_model
    
    model = cp_model.CpModel()
    
    distance_range = 20
    number_of_points = 3
    
    # create set of required (x, y) points
    x = [
        model.NewIntVar(1, distance_range, f'x_{i}')
        for i in range(number_of_points)
    ]
    y = [
        model.NewIntVar(1, distance_range, f'y_{i}')
        for i in range(number_of_points)
    ]
    
    # Build intermediate variables and get the list of square differences on each dimension.
    objective_terms = []
    for i in range(number_of_points - 1):
        for j in range(i + 1, number_of_points):
            x_diff = model.NewIntVar(-distance_range, distance_range, f'x_diff{i}')
            y_diff = model.NewIntVar(-distance_range, distance_range, f'y_diff{i}')
            model.Add(x_diff == x[i] - x[j])
            model.Add(y_diff == y[i] - y[j])
            x_diff_sq = model.NewIntVar(0, distance_range ** 2, f'x_diff_sq{i}')
            y_diff_sq = model.NewIntVar(0, distance_range ** 2, f'y_diff_sq{i}')
            model.AddMultiplicationEquality(x_diff_sq, x_diff, x_diff)
            model.AddMultiplicationEquality(y_diff_sq, y_diff, y_diff)
            objective_terms.append(x_diff_sq)
            objective_terms.append(y_diff_sq)
    
    #Objective
    model.Maximize(sum(objective_terms))
    
    # Creates a solver and solves the model.
    solver = cp_model.CpSolver()
    solver.parameters.log_search_progress = True
    solver.parameters.num_workers = 8
    status = solver.Solve(model)
    
    # Prints the solution.
    if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:
        for i in range(number_of_points):
            print(f'x_{i}={solver.Value(x[i])}')
            print(f'y_{i}={solver.Value(y[i])}')
    else:
        print('No solution found.')
    

    So many wrong things in your code:

    So, because we are using square of distances, the solution is degenerate (2 points in a corner, the other in the opposite corner.

    x_0=1
    y_0=1
    x_1=20
    y_1=20
    x_2=1
    y_2=1
    

    To fix the problem, let's introduce an integer variable that will be less that the square root of the distance between 2 robots. To increase precision, we will add a scaling factor.

    here is the intermediate code:

    objective_terms = []
    scaling = 100
    for i in range(number_of_points - 1):
        for j in range(i + 1, number_of_points):
            x_diff = model.NewIntVar(-distance_range, distance_range, f'x_diff{i}')
            y_diff = model.NewIntVar(-distance_range, distance_range, f'y_diff{i}')
            model.Add(x_diff == x[i] - x[j])
            model.Add(y_diff == y[i] - y[j])
            x_diff_sq = model.NewIntVar(0, distance_range ** 2, f'x_diff_sq{i}')
            y_diff_sq = model.NewIntVar(0, distance_range ** 2, f'y_diff_sq{i}')
            model.AddMultiplicationEquality(x_diff_sq, x_diff, x_diff)
            model.AddMultiplicationEquality(y_diff_sq, y_diff, y_diff)
    
            # we create a dist variable such that dist * dist < x_diff_sq + y_diff_sq.
            # To improve the precision, we scale x_diff_sq and y_diff_sq by a constant.
            scaled_distance = model.NewIntVar(0, distance_range * 2 * scaling, '')
            scaled_distance_square = model.NewIntVar(
                0, 2 * distance_range * distance_range * scaling * scaling, '')
            model.AddMultiplicationEquality(scaled_distance_square,
                                            scaled_distance, scaled_distance)
            model.Add(scaled_distance_square <= x_diff_sq * scaling +
                      y_diff_sq * scaling)
            objective_terms.append(scaled_distance)
    

    which gives a reasonable answer:

    x_0=1
    y_0=1
    x_1=20
    y_1=1
    x_2=1
    y_2=20
    
    

    I also rewrote the example to maximize the minimum euclidian distance between two robots. The code is available here.