numerica - Resolviendo ecuaciones no lineales en python

sistema de ecuaciones sympy

Question

Tengo 4 ecuaciones no lineales con tres incógnitas X , Y y Z que quiero resolver. Las ecuaciones son de la forma:

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

... donde a , b y c son constantes que dependen de cada valor de F en las cuatro ecuaciones.

¿Cuál es la mejor manera de resolver esto?

Answer 1

Hay dos maneras de hacer esto.

1. Utilice un solver no lineal
2. Linealice el problema y resuélvalo en el sentido de los mínimos cuadrados.

Preparar

Entonces, según entiendo su pregunta, sabe F, a, b y c en 4 puntos diferentes, y desea invertir para los parámetros de modelo X, Y y Z. Tenemos 3 incógnitas y 4 puntos de datos observados, por lo que El problema está sobredeterminado. Por lo tanto, estaremos resolviendo en el sentido de los mínimos cuadrados.

En este caso, es más común usar la terminología opuesta, así que cambiemos la ecuación. En lugar de:

F_i = X^2 + a_i Y^2 + b_i X Y cosZ + c_i X Y sinZ

Vamos a escribir:

F_i = a^2 + X_i b^2 + Y_i a b cos(c) + Z_i a b sin(c)

Donde conocemos F , X , Y y Z en 4 puntos diferentes (por ejemplo F_0, F_1, ... F_i ).

Solo estamos cambiando los nombres de las variables, no la ecuación en sí. (Esto es más para mi facilidad de pensamiento que cualquier otra cosa.)

Solución lineal

En realidad es posible linealizar esta ecuación. Puede resolver fácilmente para a^2 , b^2 , ab cos(c) y ab sin(c) . Para hacer esto un poco más fácil, volvamos a etiquetar las cosas una vez más:

d = a^2
e = b^2
f = a b cos(c)
g = a b sin(c)

Ahora la ecuación es mucho más simple: F_i = d + e X_i + f Y_i + g Z_i . Es fácil hacer una inversión lineal de mínimos cuadrados para d , e , f , y g . Entonces podemos obtener a , b , c de:

a = sqrt(d)
b = sqrt(e)
c = arctan(g/f)

Bien, escribamos esto en forma de matriz. Vamos a traducir 4 observaciones de (el código que escribiremos tendrá cualquier número de observaciones, pero seamos concretos en este momento):

F_i = d + e X_i + f Y_i + g Z_i

Dentro:

|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|

O: F = G * m (soy un geofisista, así que usamos G para "Funciones de Green" y m para "Parámetros del modelo". Normalmente también usaríamos d para "datos" en lugar de F ).

En python, esto se traduciría a:

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

Solución no lineal

También puedes resolver esto usando scipy.optimize , como sugirió @Joe. La función más accesible en scipy.optimize es scipy.optimize.curve_fit que utiliza un método de Levenberg-Marquardt de forma predeterminada.

Levenberg-Marquardt es un algoritmo de "escalada" (bueno, va cuesta abajo, en este caso, pero el término se usa de todos modos). En cierto sentido, hace una estimación inicial de los parámetros del modelo (todos unos, por defecto en scipy.optimize ) y sigue la pendiente de lo observed - predicted en su espacio de parámetros cuesta abajo hasta el final.

Advertencia: elegir el método de inversión no lineal correcto, la estimación inicial y ajustar los parámetros del método es en gran medida un "arte oscuro". Solo lo aprendes haciéndolo, y hay muchas situaciones donde las cosas no funcionan correctamente. Levenberg-Marquardt es un buen método general si su espacio de parámetros es bastante suave (este debería ser). Hay muchos otros (incluidos algoritmos genéticos, redes neuronales, etc., además de métodos más comunes como el recocido simulado) que son mejores en otras situaciones. No voy a ahondar en esa parte aquí.

Hay un problema común que algunos kits de herramientas de optimización intentan corregir para que scipy.optimize no intente manejar. Si los parámetros de su modelo tienen magnitudes diferentes (por ejemplo, a=1, b=1000, c=1e-8 ), deberá volver a escalar las cosas para que sean similares en magnitud. De scipy.optimize contrario, los scipy.optimize de "escalada" de scipy.optimize (como LM) no calcularán con precisión la estimación del gradiente local y darán resultados imprecisos. Por ahora, asumo que a , b y c tienen magnitudes relativamente similares. Además, tenga en cuenta que, en esencia, todos los métodos no lineales requieren que haga una estimación inicial y son sensibles a esa estimación. Lo dejo a continuación (simplemente pásalo como curve_fit kwarg a curve_fit ) porque el valor predeterminado a, b, c = 1, 1, 1 es una estimación bastante precisa para a, b, c = 3, 2, 1 .

Con las advertencias fuera del camino, curve_fit espera que se le pase una función, un conjunto de puntos donde se realizaron las observaciones (como una sola matriz de ndim x npoints ), y los valores observados.

Entonces, si escribimos la función así:

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

Tendremos que envolverlo para aceptar argumentos ligeramente diferentes antes de pasarlo a curve_fit .

En una palabra:

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

Ejemplo independiente de los dos métodos:

Para darle una implementación completa, aquí hay un ejemplo que

1. genera puntos distribuidos al azar para evaluar la función en,
2. evalúa la función en esos puntos (usando los parámetros del modelo establecido),
3. añade ruido a los resultados,
4. y luego invierte los parámetros del modelo utilizando los métodos lineales y no lineales descritos anteriormente.

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()

Answer 2

Probablemente desee utilizar los solucionadores no lineales de scipy, son muy fáciles: http://docs.scipy.org/doc/scipy/reference/optimize.nonlin.html