我正试图使用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])