时间戳

首页 > 解决方案 > 查找两个椭圆之间距离的Python代码？

python - 查找两个椭圆之间距离的Python代码？

问题描述

我正在尝试实现一个 Python 脚本，该脚本使用 CPLEX 查找两个椭圆之间的距离。椭圆被定义为二次曲线：满足方程的点 (x,y) 的集合：ax2 +bxy+cy2 +dx+ey+f =0，其中b2 − 4ac < 0。(x1, y1)两个椭圆之间的距离是这两个椭圆（即和(x2, y2)）内两点（即 和 ）之间的最小距离，(a1,b1,c1,d1,e1,f1)可以(a2,b2,c2,d2,e2,f2)使用以下优化问题求得：

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

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

我必须使用 .txt 文件 ellipse.txt 来表示它，其中包含两行中两个椭圆的系数：ellipse.txt：

9 0 25 -18 100 -116


16 0 9 160 -72 400

我想我必须先定义一个函数，例如，但我不知道从哪里开始。你们有什么感想？谢谢。

标签： pythonoptimizationgeometryellipseoperations-research

解决方案

您可以使用docplex 或 OPL来调用 CPLEX：最常见的 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")

欧莱雅：

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;
}

我们得到

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

您可以通过显示检查解决方案：

2个椭圆

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()

推荐阅读