¿Cómo realizar una optimización no lineal con scipy / numpy o sympy?

Estoy tratando de encontrar la solución óptima para el siguiente sistema de ecuaciones en Python:

(x-x1)^2 + (y-y1)^2 - r1^2 = 0 (x-x2)^2 + (y-y2)^2 - r2^2 = 0 (x-x3)^2 + (y-y3)^2 - r3^2 = 0 

Dados los valores, un punto (x, y) y un radio (r):

 x1, y1, r1 = (0, 0, 0.88) x2, y2, r2 = (2, 0, 1) x3, y3, r3 = (0, 2, 0.75) 

¿Cuál es la mejor manera de encontrar la solución óptima para el punto (x, y)? Utilizando el ejemplo anterior, sería:
~ (1, 1)

Si entiendo su pregunta correctamente, creo que esto es lo que está buscando:

 from scipy.optimize import minimize import numpy as np def f(coord,x,y,r): return np.sum( ((coord[0] - x)**2) + ((coord[1] - y)**2) - (r**2) ) x = np.array([0, 2, 0]) y = np.array([0, 0, 2]) r = np.array([.88, 1, .75]) # initial (bad) guess at (x,y) values initial_guess = np.array([100,100]) res = minimize(f,initial_guess,args = [x,y,r]) 

Cuyos rendimientos:

 >>> print res.x [ 0.66666666 0.66666666] 

También puede probar el método de mínimos cuadrados que espera una función objective que devuelve un vector. Quiere minimizar la sum de los cuadrados de este vector. Usando mínimos cuadrados, su función objective se vería así:

 def f2(coord,args): x,y,r = args # notice that we're returning a vector of dimension 3 return ((coord[0]-x)**2) + ((coord[1] - y)**2) - (r**2) 

Y lo minimizarías así:

 from scipy.optimize import leastsq res = leastsq(f2,initial_guess,args = [x,y,r]) 

Cuyos rendimientos:

 >>> print res[0] >>> [ 0.77961518 0.85811473] 

Esto es básicamente lo mismo que usar minimize y reescribir la función objective original como:

 def f(coord,x,y,r): vec = ((coord[0]-x)**2) + ((coord[1] - y)**2) - (r**2) # return the sum of the squares of the vector return np.sum(vec**2) 

Esto produce:

 >>> print res.x >>> [ 0.77958326 0.8580965 ] 

Tenga en cuenta que los args se manejan de forma un poco diferente con la leastsq , y que las estructuras de datos devueltas por las dos funciones también son diferentes. Consulte la documentación para scipy.optimize.minimize y scipy.optimize.leastsq para obtener más detalles.

Consulte la documentación de scipy.optimize para obtener más opciones de optimización.

Noté que el código en la solución aceptada ya no funciona … Creo que quizás scipy.optimize ha cambiado su interfaz desde que se publicó la respuesta. Podría estar equivocado. En cualquier caso, scipy.optimize la sugerencia de usar los algoritmos en scipy.optimize , y la respuesta aceptada demuestra cómo (o lo hizo al mismo tiempo, si la interfaz ha cambiado).

Estoy agregando una respuesta adicional aquí, simplemente para sugerir un paquete alternativo que use los algoritmos scipy.optimize en el núcleo, pero es mucho más robusto para la optimización restringida. El paquete es mystic . Una de las grandes mejoras es que mystic ofrece una optimización global restringida.

Primero, aquí está su ejemplo, hecho de manera muy similar a la forma scipy.optimize.minimize , pero utilizando un optimizador global.

 from mystic import reduced @reduced(lambda x,y: abs(x)+abs(y)) #choice changes answer def objective(x, a, b, c): x,y = x eqns = (\ (x - a[0])**2 + (y - b[0])**2 - c[0]**2, (x - a[1])**2 + (y - b[1])**2 - c[1]**2, (x - a[2])**2 + (y - b[2])**2 - c[2]**2) return eqns bounds = [(None,None),(None,None)] #unnecessary a = (0,2,0) b = (0,0,2) c = (.88,1,.75) args = a,b,c from mystic.solvers import diffev2 from mystic.monitors import VerboseMonitor mon = VerboseMonitor(10) result = diffev2(objective, args=args, x0=bounds, bounds=bounds, npop=40, \ ftol=1e-8, disp=False, full_output=True, itermon=mon) print result[0] print result[1] 

Con resultados parecidos a esto:

 Generation 0 has Chi-Squared: 38868.949133 Generation 10 has Chi-Squared: 2777.470642 Generation 20 has Chi-Squared: 12.808055 Generation 30 has Chi-Squared: 3.764840 Generation 40 has Chi-Squared: 2.996441 Generation 50 has Chi-Squared: 2.996441 Generation 60 has Chi-Squared: 2.996440 Generation 70 has Chi-Squared: 2.996433 Generation 80 has Chi-Squared: 2.996433 Generation 90 has Chi-Squared: 2.996433 STOP("VTRChangeOverGeneration with {'gtol': 1e-06, 'target': 0.0, 'generations': 30, 'ftol': 1e-08}") [ 0.66667151 0.66666422] 2.99643333334 

Como se señaló, la elección de la lambda en reduced efectos reduced en qué punto encuentra el optimizador ya que no hay una solución real para las ecuaciones.

mystic también proporciona la capacidad de convertir ecuaciones simbólicas a una función, donde la función resultante se puede usar como un objective, o como una función de penalización. Aquí está el mismo problema, pero usando las ecuaciones como una penalización en lugar del objective.

 def objective(x): return 0.0 equations = """ (x0 - 0)**2 + (x1 - 0)**2 - .88**2 == 0 (x0 - 2)**2 + (x1 - 0)**2 - 1**2 == 0 (x0 - 0)**2 + (x1 - 2)**2 - .75**2 == 0 """ bounds = [(None,None),(None,None)] #unnecessary from mystic.symbolic import generate_penalty, generate_conditions from mystic.solvers import diffev2 pf = generate_penalty(generate_conditions(equations), k=1e12) result = diffev2(objective, x0=bounds, bounds=bounds, penalty=pf, \ npop=40, gtol=50, disp=False, full_output=True) print result[0] print result[1] 

Con resultados:

 [ 0.77958328 0.8580965 ] 3.6473132399e+12 

Los resultados son diferentes a los anteriores porque la penalización aplicada es diferente a la que aplicamos anteriormente en reduced . En mystic , puede seleccionar qué sanción desea aplicar.

Se señaló que la ecuación no tiene solución. Se puede ver en el resultado anterior, que el resultado está fuertemente penalizado, por lo que es una buena indicación de que no hay solución. Sin embargo, el mystic tiene otra forma de ver que no hay solución. En lugar de aplicar una penalty más tradicional, que penaliza la solución donde se violan las restricciones … mystic proporciona una constraint , que es esencialmente una transformación del kernel, que elimina todas las soluciones potenciales que no cumplen con las constantes.

 def objective(x): return 0.0 equations = """ (x0 - 0)**2 + (x1 - 0)**2 - .88**2 == 0 (x0 - 2)**2 + (x1 - 0)**2 - 1**2 == 0 (x0 - 0)**2 + (x1 - 2)**2 - .75**2 == 0 """ bounds = [(None,None),(None,None)] #unnecessary from mystic.symbolic import generate_constraint, generate_solvers, simplify from mystic.symbolic import generate_penalty, generate_conditions from mystic.solvers import diffev2 cf = generate_constraint(generate_solvers(simplify(equations))) result = diffev2(objective, x0=bounds, bounds=bounds, \ constraints=cf, \ npop=40, gtol=50, disp=False, full_output=True) print result[0] print result[1] 

Con resultados:

 [ nan 657.17740835] 0.0 

Donde el nan esencialmente indica que no hay una solución válida.

Para su información, soy el autor, así que tengo un sesgo. Sin embargo, mystic ha existido casi tanto tiempo como scipy.optimize , está maduro y ha tenido una interfaz más estable durante ese tiempo. El punto es que, si necesita un optimizador no lineal restringido mucho más flexible y poderoso, sugiero mystic .

Se puede ver que estas ecuaciones describen todos los puntos en la circunferencia de tres círculos en el espacio 2D. La solución serían los puntos donde se interceptan los círculos.

La sum de sus radios de los círculos es más pequeña que las distancias entre sus centros, por lo que los círculos no se superponen. He trazado los círculos para escalar a continuación: solución geométrica

No hay puntos que satisfagan este sistema de ecuaciones.

Hice un script de ejemplo de la siguiente manera. Tenga en cuenta que la última línea encontrará una solución óptima (a, b):

 import numpy as np import scipy as scp import sympy as smp from scipy.optimize import minimize a,b = smp.symbols('a b') x_ar, y_ar = np.random.random(3), np.random.random(3) x = np.array(smp.symbols('x0:%d'%np.shape(x_ar)[0])) y = np.array(smp.symbols('y0:%d'%np.shape(x_ar)[0])) func = np.sum(a**2+b**2-x*(a+b)+2*y) print func my_func = smp.lambdify((x,y), func) print 1.0/3*my_func(x_ar,y_ar) ab = smp.lambdify((a,b),my_func(x_ar,x_ar)) print ab(1,2) def ab_v(x): return ab(*tuple(x)) print ab_v((1,2)) minimize(ab_v,(0.1,0.1)) 

Las salidas son:

 3*a**2 + 3*b**2 - x0*(a + b) - x1*(a + b) - x2*(a + b) + 2*y0 + 2*y1 + 2*y2 1.0*a**2 - 0.739792011558683*a + 1.0*b**2 - 0.739792011558683*b +0.67394435712335 12.7806239653 12.7806239653 Out[33]: status: 0 success: True njev: 3 nfev: 12 hess_inv: array([[1, 0], [0, 1]]) fun: 3.6178137388030356 x: array([ 0.36989601, 0.36989601]) message: 'Optimization terminated successfully.' jac: array([ 5.96046448e-08, 5.96046448e-08])