scipy最大限度地减少找不到解决方案的情况



我正试图使用scipy.minimize求解一组方程,但我没有得到令人满意的结果,所以可能我做错了什么。我想解下面的方程组。

12.25 * (x + y * 2.2 + z * 4.84) - 8.17437483750257 = 0
12.25 * (x + y * 3.1 + z * 9.61) - 21.9317236606432 = 0
12.25 * (x + y * 4   + z * 16)   - 107.574834524443 = 0

使用Wolfram Alpha我得到答案

x=22.626570068753, y=-17.950683342597, z=3.6223614029055

它确实解决了方程组,给出了的残差

9.407585821463726e-12

现在使用scipy.minimize我做:

import numpy as np
from scipy.optimize import fsolve
from scipy.optimize import minimize
def my_func(p):
points = [8.17437483750257, 21.9317236606432, 107.574834524443]
h1 = abs(12.25 * (p[0] + p[1] * 2.2 + p[2] * 4.84) - points[0])
h2 = abs(12.25 * (p[0] + p[1] * 3.1 + p[2] * 9.61) - points[1])
h3 = abs(12.25 * (p[0] + p[1] * 4   + p[2] * 16)   - points[2])
return h1 + h2 + h3
ini = np.array([22, -15, 5])   # Initial points close to solution
res = minimize(my_func, ini)
print(res)

fun: 1.4196640741924451
hess_inv: array([[ 20.79329103, -14.63447889,   2.36145776],
[-14.63447889,  10.30037625,  -1.66214485],
[  2.36145776,  -1.66214485,   0.26822135]])
jac: array([ 12.25      ,  60.02499545, 254.43249989])
message: 'Desired error not necessarily achieved due to precision loss.'
nfev: 261
nit: 8
njev: 64
status: 2
success: False
x: array([ 21.39197235, -17.08623345,   3.48344393])

第一,它说成功=错误,第二,它找到了不是最佳的解决方案。

为什么具有接近最优解的初始值却无法找到这些解。

优化器的定义有问题吗?

尝试运行它,给出[0,0,0]的初始值,但它只会给出糟糕的结果

ini = np.array([0, 0, 0])   # Initial points close to solution
res = minimize(my_func, ini)
print(res)
fun: 73.66496363902732
hess_inv: array([[ 0.98461683, -0.04223651, -0.1207056 ],
[-0.04223651,  0.88596592, -0.31885642],
[-0.1207056 , -0.31885642,  0.13448927]])
jac: array([ 12.25      ,  15.92499924, -18.98750019])
message: 'Desired error not necessarily achieved due to precision loss.'
nfev: 164
nit: 1
njev: 40
status: 2
success: False
x: array([0.02901304, 0.08994042, 0.29448233])

注意:我不想使用fsolve来找到解决方案,而是使用minimize。原因是我真正的问题涉及到有更多的方程而不是未知数,所以最后我想要一个将所有方程的误差最小化的解决方案。然而,由于它没有给出好的结果,我想首先测试一个简单的问题,这个问题有一个确切的解决方案。但即使在这种情况下,它也不起作用。一旦我使它适用于这个问题,我会扩展它,添加更多的方程。

。。。我真正的问题是有更多的方程而不是未知数,所以最后我想要一个解决方案,最大限度地减少所有这些方程的误差

这听起来很像广义矩量法(GMM(中解决的问题,在GMM中,方程比未知数多。

这类问题通常使用最小二乘法来解决。假设你的整个系统是这样的:

h1(x, y, z) = 0
h2(x, y, z) = 0
h3(x, y, z) = 0
h4(x, y, z) = 0

它有3个未知数和4个方程。那么你的目标函数将是:

F(x, y, z) = H(x, y, z)' * W * H(x, y, z)
  • H(x, y, z)是以上所有hj(x, y, z)的向量
  • H(x, y, z)'是它的转置
  • W是加权矩阵

如果W是单位矩阵,则得到最小二乘目标函数。那么,F(x, y, z)是一个二次型(基本上是多维的抛物线(,它应该很容易优化,因为它是凸的和光滑的。

您的代码使用像h1 = abs(12.25 * (p[0] + p[1] * 2.2 + p[2] * 4.84) - points[0])这样的绝对值,但abs可能很难在原点附近进行区分,但这正是您的最佳值所在,因为您本质上希望h1等于零。

你可以通过对误差进行平方来近似绝对值函数:

h1 =(12.25 * (p[0] + p[1] * 2.2 + p[2] * 4.84) - points[0])**2

这导致了与GMM(或最小二乘法(基本相同的方法,并为您提供了一个易于优化的函数,因为正方形在原点附近是平滑的。

优化问题(和求解者(通常受益于表现良好(平滑(的"优化表面";。当您使用abs函数时;波涛汹涌的";曲面上,带有点的导数是不连续的。

如果不使用abs,而是使用二次函数(具有相同效果(,则可以得到接近预期的解决方案。只需将my_func更改为:

def my_func(p):
points = [8.17437483750257, 21.9317236606432, 107.574834524443]
h1 = (12.25 * (p[0] + p[1] * 2.2 + p[2] * 4.84) - points[0])**2
h2 = (12.25 * (p[0] + p[1] * 3.1 + p[2] * 9.61) - points[1])**2
h3 = (12.25 * (p[0] + p[1] * 4   + p[2] * 16)   - points[2])**2
return h1 + h2 + h3

我得到的是:

fun: 8.437863292878727e-10
hess_inv: array([[ 0.64753863, -0.43474506,  0.06909179],
[-0.43474506,  0.29487798, -0.04722923],
[ 0.06909179, -0.04722923,  0.00761762]])
jac: array([-6.26698693e-12,  6.22490927e-10, -5.11716516e-10])
message: 'Optimization terminated successfully.'
nfev: 55
nit: 7
njev: 11
status: 0
success: True
x: array([ 22.62653789, -17.95066124,   3.62235782])

相关内容

最新更新