convex-concave-scatterers.py

import numpy as np
import cvxpy as cp
import matplotlib.pyplot as plt

# The idea is that we have two circles, in opposite corners
# that need to move past each other during optimization (in order to minimize
# some simple objective function). The important thing to guarantee
# is that they do not intersect at any point during any iteration.

# Initial/desired positions of circles, radius = 1/4.
x_1 = np.array([0, .4])
x_2 = np.array([1, .6])

x_1_desired = np.array([1, .5])
x_2_desired = np.array([0, .5])

r = 1/4

all_positions_1 = [x_1]
all_positions_2 = [x_2]


for curr_iter in range(10):
    pos_1 = cp.Variable(2)
    pos_2 = cp.Variable(2)

    # Can, of course, be any objective w.r.t. the scatterer positions ! Would say that
    # this could be projected gradient descent, e.g.
    obj = cp.Minimize(cp.sum_squares(x_1_desired - pos_1) + cp.sum_squares(x_2_desired - pos_2))

    # Uncomment line and change iterations from 10 -> 100 to see more "continuous" but slower version, that makes the avoidance clearer in the plot:
    # obj = cp.Minimize(cp.sum_squares(x_1_desired - pos_1) + cp.sum_squares(x_2_desired - pos_2) + 10*cp.sum_squares(pos_1 - x_1) + 10*cp.sum_squares(pos_2 - x_2))

    # We can write the hyperplane lying in between x_1 and x_2 as
    # the set of all y such that ||y - x_1|| <= ||y - x_2||,
    # which is equivalently ||x_1||^2 - ||x_2||^2 <= 2 * (x_1 - x_2)' * y.
    # We want to add a margin of at least r, so this is just:
    # ||x_1||^2 - ||x_2||^2 + r^2 <= 2 * (x_1 - x_2)' * y.
    cons = [
        cp.sum_squares(x_1) - cp.sum_squares(x_2) + r*r <= 2 * (x_1 - x_2) @ pos_1,
        cp.sum_squares(x_1) - cp.sum_squares(x_2) - r*r >= 2 * (x_1 - x_2) @ pos_2
    ]

    prob = cp.Problem(obj, cons)
    print(f'Current objective : {prob.solve(solver=cp.ECOS):.3f}')

    x_1 = pos_1.value
    x_2 = pos_2.value

    all_positions_1.append(x_1)
    all_positions_2.append(x_2)

all_pos_1 = np.array(all_positions_1)
all_pos_2 = np.array(all_positions_2)

plt.figure()
plt.scatter(all_pos_1[:, 0], all_pos_1[:, 1])
plt.scatter(all_pos_2[:, 0], all_pos_2[:, 1])
plt.show()