Как выполнить нелинейную оптимизацию с помощью scipy / numpy или sympy?

Я пытаюсь найти оптимальное решение следующей системы уравнений в 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

При заданных значениях точки (x, y) и радиуса (r): 

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

Как лучше всего найти оптимальное решение для точки (x, y) Используя приведенный выше пример, это будет:
~ (1, 1)

4 8
python numpy scipy sympy

4 ответа:

Если я правильно понял ваш вопрос, я думаю, что это то, что вам нужно: 

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

Который дает: 

>>> print res.x
[ 0.66666666  0.66666666]

Вы также можете попробовать метод наименьших квадратов, который ожидает целевую функцию, возвращающую вектор. Он хочет минимизировать сумму квадратов этого вектора. Используя наименьшие квадраты, ваша целевая функция будет выглядеть следующим образом: 

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)

И вы бы минимизировали его следующим образом: 

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

Который дает: 

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

Это в основном то же самое, что использовать minimize и переписать исходную целевую функцию как: 

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)

Это дает: 

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

Обратите внимание, что args обрабатываются немного по-разному с leastsq, и что структуры данных, возвращаемые этими двумя функциями, также различны. Смотрите документацию для scipy.optimize.minimize и еще scipy.optimize.leastsq для более подробной информации.

См.scipy.optimize документация для дополнительных вариантов оптимизации.

14

Я заметил, что код в принятом решении больше не работает... Я думаю, что, возможно, scipy.optimize изменил свой интерфейс с тех пор, как был опубликован ответ. Я могу ошибаться. Несмотря на это, я поддерживаю предложение использовать алгоритмы в scipy.optimize, и принятый ответ действительно демонстрирует, как (или делал когда-то, если интерфейс изменился).

Я добавляю здесь дополнительный ответ, чисто чтобы предложить альтернативный пакет, который использует алгоритмы scipy.optimize в ядре, но гораздо больше робастный для ограниченной оптимизации. Пакет - это mystic. Одним из больших улучшений является то, что mystic дает ограниченную глобальную оптимизацию.

Во-первых, вот ваш пример, выполненный очень похоже на способ scipy.optimize.minimize, но с использованием глобального оптимизатора. 
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]

С результатами, которые выглядят следующим образом: 

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
Как уже отмечалось, выбор lambda в reduced влияет на то, какую точку находит оптимизатор, поскольку фактического решения уравнений не существует.

mystic также обеспечивает возможность преобразования символьных уравнений в функцию, где результирующая функция может использоваться в качестве цели или в качестве штрафной функции. Здесь та же проблема, но с использованием уравнений в качестве штрафа вместо цели.

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]

С результатами: 

[ 0.77958328  0.8580965 ]
3.6473132399e+12
Результаты отличаются от предыдущих, потому что примененное наказание отличается от того, что мы применяли ранее в reduced. В mystic вы можете выбрать, какое наказание вы хотите применить.

Эта точка зрения была высказана что уравнение не имеет решения. Вы можете видеть из приведенного выше результата, что результат сильно наказан, так что это хороший признак того, что нет никакого решения. Однако у mystic есть другой способ, который вы можете увидеть там в отсутствии решения. Вместо применения более традиционного penalty, который наказывает решение, где ограничения нарушаются... mystic обеспечивает constraint, которое по сути является преобразованием ядра, которое удаляет все потенциальные решения, которые не соответствуют константам. 

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]

С результатами: 

[          nan  657.17740835]
0.0

, где nan по существу указывает на отсутствие допустимого решения.

К вашему сведению, я автор, так что у меня есть некоторая предвзятость. Однако mystic существует почти столько же, сколько scipy.optimize, является зрелым и имеет более стабильный интерфейс в течение этого периода времени. Суть в том, что если вам нужен гораздо более гибкий и мощный нелинейный оптимизатор с ограниченными возможностями, я предлагаю mystic.
3

Эти уравнения можно рассматривать как описывающие все точки на окружности трех окружностей в двумерном пространстве. Решение будет заключаться в точках, где круги пересекаются.

Сумма их радиусов окружностей меньше, чем расстояния между их центрами, поэтому окружности не перекрываются. Я нарисовал круги, чтобы масштабировать их ниже: геометрическое решение

Нет точек, удовлетворяющих этой системе уравнений.

1

Я сделал пример сценария следующим образом. Заметим, что в последней строке будет найдено оптимальное решение (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))

Выходы: 

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])
0