我有三个未知数4非线性方程X
,Y
和Z
我想要解决。等式的形式为:
F(m) = X^2 + a(m)Y^2 + b(m)XYcosZ + c(m)XYsinZ
…其中a
,b
和c
是取决于F
四个方程式中每个值的常数。
解决此问题的最佳方法是什么?
有两种方法可以做到这一点。
因此,据我所知,您知道在4个不同点处的F,a,b和c,并且想要对模型参数X,Y和Z求逆。我们有3个未知数和4个观测数据点,因此这个问题太确定了。因此,我们将在最小二乘意义上进行求解。
在这种情况下,使用相反的术语更为常见,因此让我们翻转方程式。代替:
F_i = X^2 + a_i Y^2 + b_i X Y cosZ + c_i X Y sinZ
让我们写:
F_i = a^2 + X_i b^2 + Y_i a b cos(c) + Z_i a b sin(c)
我们知道F
,X
,Y
,并Z
在4个不同点(例如F_0, F_1, ... F_i
)。
我们只是在更改变量的名称,而不是方程式本身。(这比我想的要容易得多。)
实际上有可能使该方程线性化。您可以轻松地解决a^2
,b^2
,a b cos(c)
,和a b sin(c)
。为了使操作更简单,让我们再次重新标记:
d = a^2
e = b^2
f = a b cos(c)
g = a b sin(c)
现在,方程式变得简单得多:F_i = d + e X_i + f Y_i + g Z_i
。这很容易做到最小二乘线性反演d
,e
,f
,和g
。然后a
,我们可以从中获取b
,和c
:
a = sqrt(d)
b = sqrt(e)
c = arctan(g/f)
好吧,让我们用矩阵形式写出来。我们将翻译4个观测值(我们将编写的代码将采用任意数量的观测值,但现在让我们保持具体):
F_i = d + e X_i + f Y_i + g Z_i
进入:
|F_0| |1, X_0, Y_0, Z_0| |d|
|F_1| = |1, X_1, Y_1, Z_1| * |e|
|F_2| |1, X_2, Y_2, Z_2| |f|
|F_3| |1, X_3, Y_3, Z_3| |g|
或:(F = G * m
我是地球物理学家,所以我们将其G
用于“格林函数”和m
“模型参数”。通常,我们也将d
“数据”而不是F
)。
在python中,这将转换为:
def invert(f, x, y, z):
G = np.vstack([np.ones_like(x), x, y, z]).T
m, _, _, _ = np.linalg.lstsq(G, f)
d, e, f, g = m
a = np.sqrt(d)
b = np.sqrt(e)
c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
return a, b, c
您也可以使用scipy.optimize
@Joe建议使用来解决此问题。其中最易访问的函数scipy.optimize
是scipy.optimize.curve_fit
默认情况下使用Levenberg-
Marquardt方法。
Levenberg-
Marquardt是一种“爬坡”算法(在这种情况下,它会走下坡路,但是无论如何都使用该术语)。从某种意义上说,您可以初步估计模型参数(所有参数,默认情况下为scipy.optimize
),然后沿observed - predicted
参数空间的坡度向下倾斜至底部。
警告:
选择正确的非线性反演方法,初步猜测并调整方法的参数在很大程度上是一门“黑手艺”。您只能通过做来学习它,并且在很多情况下,事情将无法正常进行。如果您的参数空间相当平滑(应该是这样),那么Levenberg-
Marquardt是一个很好的通用方法。在其他情况下,还有很多其他方法(包括遗传算法,神经网络等,除了更常见的方法(如模拟退火)之外)。我在这里不打算深入探讨这一部分。
有一个常见的陷阱,一些优化工具包试图纠正这一错误,而scipy.optimize
没有解决。如果模型参数的大小不同(例如a=1, b=1000, c=1e-8
),则需要重新缩放比例,以使它们的大小相似。否则scipy.optimize
,“爬山”算法(如LM)将无法准确计算局部梯度的估算值,并且会给出非常不准确的结果。现在,我假设a
,b
以及c
具有相对类似的量值。另外,请注意,基本上所有非线性方法都需要您进行初始猜测,并且对此猜测很敏感。我将其保留在下面(将其作为p0
kwarg传递给它curve_fit
),因为默认值a, b, c = 1, 1, 1
是对的相当准确的猜测a, b, c = 3, 2, 1
。
避免了警告,curve_fit
希望传递一个函数,进行观测的一组点(作为单个ndim x npoints
数组)和观测值。
因此,如果我们这样编写函数:
def func(x, y, z, a, b, c):
f = (a**2
+ x * b**2
+ y * a * b * np.cos(c)
+ z * a * b * np.sin(c))
return f
在将其传递给之前,我们需要将其包装为接受稍有不同的参数curve_fit
。
简而言之:
def nonlinear_invert(f, x, y, z):
def wrapped_func(observation_points, a, b, c):
x, y, z = observation_points
return func(x, y, z, a, b, c)
xdata = np.vstack([x, y, z])
model, cov = opt.curve_fit(wrapped_func, xdata, f)
return model
为了给您一个完整的实现,下面是一个示例
import numpy as np
import scipy.optimize as opt
def main():
nobservations = 4
a, b, c = 3.0, 2.0, 1.0
f, x, y, z = generate_data(nobservations, a, b, c)
print 'Linear results (should be {}, {}, {}):'.format(a, b, c)
print linear_invert(f, x, y, z)
print 'Non-linear results (should be {}, {}, {}):'.format(a, b, c)
print nonlinear_invert(f, x, y, z)
def generate_data(nobservations, a, b, c, noise_level=0.01):
x, y, z = np.random.random((3, nobservations))
noise = noise_level * np.random.normal(0, noise_level, nobservations)
f = func(x, y, z, a, b, c) + noise
return f, x, y, z
def func(x, y, z, a, b, c):
f = (a**2
+ x * b**2
+ y * a * b * np.cos(c)
+ z * a * b * np.sin(c))
return f
def linear_invert(f, x, y, z):
G = np.vstack([np.ones_like(x), x, y, z]).T
m, _, _, _ = np.linalg.lstsq(G, f)
d, e, f, g = m
a = np.sqrt(d)
b = np.sqrt(e)
c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
return a, b, c
def nonlinear_invert(f, x, y, z):
# "curve_fit" expects the function to take a slightly different form...
def wrapped_func(observation_points, a, b, c):
x, y, z = observation_points
return func(x, y, z, a, b, c)
xdata = np.vstack([x, y, z])
model, cov = opt.curve_fit(wrapped_func, xdata, f)
return model
main()
我正在研究一个非线性微分方程求解器。我能得到一般的解决方案,但不能得到具体的解决方案。当我试图找到集成常量时,我得到了错误:E_x不可调用的,我的解决方案被归类为一个列表,所以我不能在其中替换任何东西。 这是我的代码: 我在上得到了可调用错误,在我在书中的例子(数值Python)解决了一个线性ODE,但没有得到这个错误。有没有一种不同的方法来解决非线性DE,或者将E定义为可调用的,而将E_sol定
我正试图找到一个用Python解决一个非线性超定系统的好方法。我在这里查看了优化工具http://docs.cipy.org/doc/scipy/reference/optimize.nonlin.html,但我不知道如何使用它们。到目前为止我所掌握的是 你能帮我一下吗?
我正在用pyomo编程求解非线性优化问题(使用ipopt求解器)。稍后,我想在模型中添加随机元素。我知道在Pyomo中,可以使用复数形式来处理随机规划,但复数形式只能处理线性规划、混合整数规划和二次规划。 一般非线性随机规划问题有求解器吗?如果没有,我们如何使用现有的求解器来处理它?
问题内容: 我想用三个或更多变量求解线性方程。python中有一个好的库吗? 问题答案: 参见http://sympy.org/和http://numpy.scipy.org/。 具体来说,http://docs.scipy.org/doc/numpy/reference/routines.linalg.html 和http://docs.sympy.org/0.7.0/tutorial.html
问题内容: 假设具有以下功能: 有了这个功能跨越的和。 这是通过以下代码确定的: 使用(或任何其他Python根查找器)一次调用找到两个根的正确方法是什么?有其他功能吗? 参考 问题答案: 定义函数,使其可以使用标量或numpy数组作为参数: 然后将参数向量传递给。
我很难理解非线性优化中的性能如何受到求解器引擎接口的特定方式的影响。 我们有一个优化模型,在它的第一个版本中,是用GAMS编写的。IPOPT(一个常见的FOOS非线性求解器引擎)在IPOPT(无函数评估)中为每个优化返回1.4 CPU秒的执行时间,在函数评估中返回0.2 CPU秒的执行时间。 当我们将模型转换为C(以便更好地考虑模型的非优化组件)并通过其C API接口IPOPT时(使用ADOL-C