python - 在python中求解非线性方程

Question

我有4个非线性方程，要解决的三个未知数X，Y和Z。等式的形式为:

F(m) = X^2 + a(m)Y^2 + b(m)XYcosZ + c(m)XYsinZ

...其中 a， b和 c是常数，它们取决于四个方程式中 F的每个值。

解决此问题的最佳方法是什么？

Answer 1

有两种方法可以做到这一点。

使用非线性求解器

线性化问题并用最小二乘法解决问题

设定

因此，据我所知，您知道在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)

我们在4个不同的点(例如 F)知道 X， Y， Z和 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()