如何用scipy/numpy或sympy进行非线性优化？

Question

我试图找到Python中以下方程组的最优解:

(x-x1)^2 + (y-y1)^2 - r1^2 = 0
(x-x2)^2 + (y-y2)^2 - r2^2 = 0
(x-x3)^2 + (y-y3)^2 - r3^2 = 0

给定点(x,y)和半径(r)的值:

x1, y1, r1 = (0, 0, 0.88)
x2, y2, r2 = (2, 0, 1)
x3, y3, r3 = (0, 2, 0.75)

找到点(x,y)的最优解的最佳方法是使用上面的例子:
〜(1,1)

Answer 1

如果我理解你的问题,我认为这就是你所追求的:

from scipy.optimize import minimize
import numpy as np

def f(coord,x,y,r):
    return np.sum( ((coord[0] - x)**2) + ((coord[1] - y)**2) - (r**2) )

x = np.array([0,   2,  0])
y = np.array([0,   0,  2])
r = np.array([.88, 1, .75])

# initial (bad) guess at (x,y) values
initial_guess = np.array([100,100])

res = minimize(f,initial_guess,args = [x,y,r])

产量:

>>> print res.x
[ 0.66666666  0.66666666]

您也可以尝试使用最小二乘法,该方法需要一个返回向量的目标函数.它希望最小化该向量的平方和.使用最小二乘法,您的目标函数将如下所示:

def f2(coord,args):
    x,y,r = args
    # notice that we're returning a vector of dimension 3
    return ((coord[0]-x)**2) + ((coord[1] - y)**2) - (r**2)

而你最小化它是这样的:

from scipy.optimize import leastsq
res = leastsq(f2,initial_guess,args = [x,y,r])

产量:

>>> print res[0]
>>> [ 0.77961518  0.85811473]

这与使用minimize和重写原始目标函数基本相同:

def f(coord,x,y,r):
    vec = ((coord[0]-x)**2) + ((coord[1] - y)**2) - (r**2)
    # return the sum of the squares of the vector
    return np.sum(vec**2)

这会产生:

>>> print res.x
>>> [ 0.77958326  0.8580965 ]

请注意,args处理方式略有不同leastsq,并且两个函数返回的数据结构也不同.查看文档scipy.optimize.minimize和scipy.optimize.leastsq更多的细节.

有关scipy.optimize更多优化选项,请参阅文档.

Answer 2

我注意到接受的解决方案中的代码不再起作用......我认为scipy.optimize自从答案发布以来它的界面可能已经改变了。我可能是错的。无论如何，我支持在 中使用算法的建议scipy.optimize，并且接受的答案确实演示了如何（或曾经做过，如果界面已更改）。

我在这里添加了一个额外的答案，纯粹是为了建议一个使用scipy.optimize核心算法的替代包，但对于约束优化来说更加健壮。包裹是mystic. 其中一项重大改进是mystic提供了受约束的全局优化。

首先，这是您的示例，与scipy.optimize.minimize方式非常相似，但使用全局优化器。

from mystic import reduced

@reduced(lambda x,y: abs(x)+abs(y)) #choice changes answer
def objective(x, a, b, c):
  x,y = x
  eqns = (\
    (x - a[0])**2 + (y - b[0])**2 - c[0]**2,
    (x - a[1])**2 + (y - b[1])**2 - c[1]**2,
    (x - a[2])**2 + (y - b[2])**2 - c[2]**2)
  return eqns

bounds = [(None,None),(None,None)] #unnecessary

a = (0,2,0)
b = (0,0,2)
c = (.88,1,.75)
args = a,b,c

from mystic.solvers import diffev2
from mystic.monitors import VerboseMonitor
mon = VerboseMonitor(10)

result = diffev2(objective, args=args, x0=bounds, bounds=bounds, npop=40, \ 
                 ftol=1e-8, disp=False, full_output=True, itermon=mon)

print result[0]
print result[1]

结果如下所示：

Generation 0 has Chi-Squared: 38868.949133
Generation 10 has Chi-Squared: 2777.470642
Generation 20 has Chi-Squared: 12.808055
Generation 30 has Chi-Squared: 3.764840
Generation 40 has Chi-Squared: 2.996441
Generation 50 has Chi-Squared: 2.996441
Generation 60 has Chi-Squared: 2.996440
Generation 70 has Chi-Squared: 2.996433
Generation 80 has Chi-Squared: 2.996433
Generation 90 has Chi-Squared: 2.996433
STOP("VTRChangeOverGeneration with {'gtol': 1e-06, 'target': 0.0, 'generations': 30, 'ftol': 1e-08}")
[ 0.66667151  0.66666422]
2.99643333334

如前所述，lambdain的选择reduced会影响优化器找到的点，因为方程没有实际解。

mystic还提供了将符号方程转换为函数的能力，其中结果函数可以用作目标或惩罚函数。这是同样的问题，但使用方程作为惩罚而不是目标。

def objective(x):
    return 0.0

equations = """
(x0 - 0)**2 + (x1 - 0)**2 - .88**2 == 0
(x0 - 2)**2 + (x1 - 0)**2 - 1**2 == 0
(x0 - 0)**2 + (x1 - 2)**2 - .75**2 == 0
"""

bounds = [(None,None),(None,None)] #unnecessary

from mystic.symbolic import generate_penalty, generate_conditions
from mystic.solvers import diffev2

pf = generate_penalty(generate_conditions(equations), k=1e12)

result = diffev2(objective, x0=bounds, bounds=bounds, penalty=pf, \
                 npop=40, gtol=50, disp=False, full_output=True)

print result[0]
print result[1]

结果：

[ 0.77958328  0.8580965 ]
3.6473132399e+12

结果与之前不同，因为应用的惩罚与我们之前应用的不同reduced。在 中mystic，您可以选择要应用的惩罚。

指出该方程无解。您可以从上面的结果中看到，结果受到了严重的惩罚，这很好地表明没有解决方案。但是，mystic您可以通过另一种方式在没有解决方案的情况下看到那里。而不是应用更传统的penalty，它会惩罚违反约束的解决方案......mystic提供了一个constraint，它本质上是一个内核转换，它删除了所有不满足常数的潜在解决方案。

def objective(x):
    return 0.0

equations = """
(x0 - 0)**2 + (x1 - 0)**2 - .88**2 == 0
(x0 - 2)**2 + (x1 - 0)**2 - 1**2 == 0
(x0 - 0)**2 + (x1 - 2)**2 - .75**2 == 0
"""

bounds = [(None,None),(None,None)] #unnecessary

from mystic.symbolic import generate_constraint, generate_solvers, simplify
from mystic.symbolic import generate_penalty, generate_conditions    
from mystic.solvers import diffev2

cf = generate_constraint(generate_solvers(simplify(equations)))

result = diffev2(objective, x0=bounds, bounds=bounds, \
                 constraints=cf, \
                 npop=40, gtol=50, disp=False, full_output=True)

print result[0]
print result[1]

结果：

[          nan  657.17740835]
0.0

当nan基本表示没有有效的解决方案。

仅供参考，我是作者，所以我有一些偏见。然而，mystic它已经存在了几乎和 一样长的时间scipy.optimize，并且已经成熟，并且在这段时间内拥有更稳定的界面。关键是，如果您需要一个更加灵活和强大的约束非线性优化器，我建议使用mystic.