Python code that finds the distance between two ellipses?

Question

I'm trying to implement a Python script that finds the distance between two ellipses using CPLEX. An ellipse is defined as a quadric: the set of points (x,y) that satisfy the equation: ax2 +bxy+cy2 +dx+ey+f =0, where b2 − 4ac < 0. The distance between two ellipses is the minimum distance between two points (i.e., (x1, y1) and (x2, y2)) within these two ellipses (i.e., (a1,b1,c1,d1,e1,f1) and (a2,b2,c2,d2,e2,f2)) can be found using the following optimization problem:

minimize z = (x1 − x2)^2 + (y1 − y2)^2

subject to: a1x21 + b1x1y1 + c1y12 + d1x1 + e1y1 + f1 ≤ 0
            a2x2 + b2x2y2 + c2y2 + d2x2 + e2y2 + f2 ≤ 0

I have to represent it by using a .txt file ellipse.txt which contains the coefficients of two ellipses in two rows: ellipse.txt:

9 0 25 -18 100 -116


16 0 9 160 -72 400

I guess I have to define a function first such as but I couldn't figure out where to start. What do you guys think? Thanks.

"where to start" -- Are you asking how to program? How to do numerical analysis to minimize something? Or something else? — Rick James, Jan 13 '21 at 22:41
Check out this SciPy [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) for the `scipy.optimize.minimize()` function — Jacob K, Jan 13 '21 at 22:43
Can you just use [QCQP](https://en.wikipedia.org/wiki/Quadratically_constrained_quadratic_program) solver for Python? Such as https://github.com/cvxgrp/qcqp — fdermishin, Jan 14 '21 at 06:05
why not use parametric equation? then "brute" force the 2 angles with some step (like 10deg) and remember the smallest distance points, then search again around found solution with smaller step (like 1deg) and repeat until wanted accuracy is met ... if [coded right](https://stackoverflow.com/a/36163847/2521214) it could be `O(log^2(n))` where `n` is `360deg/smallest_step` — Spektre, Jan 14 '21 at 08:47

Alex Fleischer · Answer 1 · 2021-01-15T09:00:02.783

you could use docplex or OPL to call CPLEX : most common CPLEX API

docplex python:

from docplex.mp.model import Model

M=10000

a1=9
b1= 0
c1= 25
d1= -18
e1= 100
f1= -116

a2=16
b2= 0
c2= 9
d2= 160
e2= -72
f2= 400

mdl = Model(name='ellipse')
x1 = mdl.continuous_var(name='x1',lb=-M,ub=M)
y1 = mdl.continuous_var(name='y1',lb=-M,ub=M)
x2 = mdl.continuous_var(name='x2',lb=-M,ub=M)
y2 = mdl.continuous_var(name='y2',lb=-M,ub=M)


mdl.add_constraint(a1*x1*x1+ b1*x1*y1 + c1*y1*y1 + d1*x1 + e1*y1 + f1 <= 0)
mdl.add_constraint(a2*x2*x2+ b2*x2*y2 + c2*y2*y2 + d2*x2 + e2*y2 + f2 <= 0)
mdl.minimize((x1-x2)*(x1-x2) +(y1-y2)*(y1-y2))

if (mdl.solve()):
   mdl.print_information()
   for v in mdl.iter_continuous_vars():
       print(v," = ",v.solution_value)
else:
        print("Problem has no solution")

OPL:

int a1=9;
int b1= 0;
int c1= 25;
int d1= -18;
int e1= 100;
int f1= -116;

int a2=16;
int b2= 0;
int c2= 9;
int d2= 160;
int e2= -72;
int f2= 400;

dvar float x1;
dvar float x2;
dvar float y1;
dvar float y2;

minimize (x1-x2)^2 +(y1-y2)^2;

subject to
{
  
a1*x1^2+ b1*x1*y1 + c1*y1^2 + d1*x1 + e1*y1 + f1 <= 0;
a2*x2^2+ b2*x2*y2 + c2*y2^2 + d2*x2 + e2*y2 + f2 <= 0;
}

and we get

x1  =  -2.951172430117367
y1  =  -0.16158545959478943
x2  =  -3.494259698129099
y2  =  0.5403290464128077

and you may check the solution with a display:

from sympy import *

a1=9
b1= 0
c1= 25
d1= -18
e1= 100
f1= -116

a2=16
b2= 0
c2= 9
d2= 160
e2= -72
f2= 400


x1b, y1b = symbols('x y')
x2b, y2b = symbols('x y')
p1=plot_implicit(a1*x1b*x1b+ b1*x1b*y1b + c1*y1b*y1b + d1*x1b + e1*y1b + f1 <= 0,
              (x1b,-10,10),(y1b,-10,10))
p2=plot_implicit(a2*x2b*x2b+ b2*x2b*y2b + c2*y2b*y2b + d2*x2b + e2*y2b + f2 <= 0,
                 (x2b,-10,10),(y2b,-10,10))

p1.extend(p2)


p1.show()

Python code that finds the distance between two ellipses?

1 Answers1