如何用 scipy 解决非线性优化问题

Question

我需要解决 Python 中的一个非线性优化问题。我发现 scipy 解决了优化问题，但是我不知道我做错了什么，因为使用一些示例输入它找不到我在 NEOS server solver Knitro AMPL.[=28 中的正确解决方案=]

我的问题是，给定一组点，它必须找到最大的内接椭圆，该椭圆最多接触这些点，并且这些点永远不会包含在其中。

理论

制定优化问题，我有 a 和 b 半轴，phi 旋转，xc 和 yc 中心坐标和 points 每个元素的点列表，形式为 [x, y] -> [0, 1] 索引。

纸面上的问题和约束是这些，a, b, phi, xc, yc是真实的，点是整数：

NEOS

我在 NEOS 中使用的文件是这些：

mod

dat

run

成功的结果（complete）：

xc = 143.012
yc = 262.634
a = 181.489
b = 140.429
phi = 1.43575

Python

所以，我的python代码是这样的，这是我第一次使用scipy进行优化，所以我不排除从文档中理解它是如何工作的错误。

from typing import List
import numpy as np
from scipy.optimize import *


def ellipse_calc(
    points: List[List[int]],
    verbose: bool = False
):
    centre = [0, 0]
    for i in range(len(points)):
        centre[0] += points[i][0]
        centre[1] += points[i][1]
    centre[0] /= len(points)
    centre[1] /= len(points)
    if verbose:
        print(f'centre: {centre[0]:.2f}, {centre[1]:.2f}')
    max_x = max([p[0] for p in points])
    max_y = max([p[1] for p in points])
    min_x = min([p[0] for p in points])
    min_y = min([p[1] for p in points])
    initial_axis = 0.25 * (max_x - min_x + max_y - min_y)
    if verbose:
        print(initial_axis)

    constraints = [
        NonlinearConstraint(lambda x: x[0], 1, np.inf),
        NonlinearConstraint(lambda x: x[1], 1, np.inf),
        NonlinearConstraint(lambda x: x[2], 0, np.inf),
    ]
    for i in range(len(points)):
        constraints += [NonlinearConstraint(
            lambda x:
            (points[i][0] - x[3]) ** 2 * (np.cos(x[2]) ** 2 / x[0]**2 + np.sin(x[2]) ** 2 / x[1]**2) +
            (points[i][1] - x[4]) ** 2 * (np.sin(x[2]) ** 2 / x[0]**2 + np.cos(x[2]) ** 2 / x[1]**2) +
            2 * (points[i][0] - x[3]) * (points[i][1] - x[4]) *
            np.cos(x[2]) * np.sin(x[2]) * (1 / x[1]**2 - 1 / x[0]**2), 1, np.inf)]
    result = minimize(
        lambda x: -np.pi * x[0] * x[1],
        [initial_axis, initial_axis, 0, centre[0], centre[1]],
        constraints=constraints
    )
    print(result)


if __name__ == '__main__':
    points = [[50,44],[91,44],[161,44],[177,44],[44,88],[189,88],[239,88],[259,88],[2,132],[250,132],[2,176],[329,176],[2,220],[289,220],[2,264],[288,264],[2,308],[277,308],[2,352],[285,352],[2,396],[25,396],[35,396],[231,396],[284,396],[298,396],[36,440],[76,440],[106,440],[173,440]]
    ellipse_calc(points, True)

这次尝试，它具有我在 NEOS 上尝试过的相同数据，输出如下：

     fun: -8.992626773255127e+40
     jac: array([-5.68832805e+20, -4.96651566e+20, -0.00000000e+00, -0.00000000e+00,
       -0.00000000e+00])
 message: 'Inequality constraints incompatible' 
    nfev: 54
     nit: 10
    njev: 9
  status: 4
 success: False
       x: array([ 1.58089104e+20,  1.81065104e+20, -1.24564497e+15, -1.55647883e+10,
       -2.76654483e+10])

有谁知道我做错了什么以及如何解决？另外，我真的不知道是否可以用 scipy 解决这个问题，在那种情况下，我正在寻找一个免费的库来解决它，甚至寻找寻找椭圆方程的替代方法

Answer 1

这不是一个完整的答案，但它应该可以帮助您入门。这里有两个提示：

传递变量上的简单框约束作为边界，而不是约束。即使用

bounds = [(1, None), (1, None), (0, None), (None, None), (None, None)]

并通过 bounds 参数传递给 minimize。

在循环内通过 lambda 表达式定义约束时需要非常小心，请参阅 here。您需要通过lambda x, i=i: your_fun捕获循环变量i。否则，您的每个约束都使用 i=29 并因此评估最后一点。通过评估特定值的所有约束可以很容易地观察到这一点。

那么您至少应该得到一个 objective 值为 79384 的可行解决方案。另请注意，您可以通过使用 numpy 函数而不是循环来显着缩短代码。

如何用 scipy 解决非线性优化问题

How to solve non linear optimization problem with scipy

scipy

python-3.x

scipy-optimize