Это может быть просто и эффективно, поверьте мне.
Независимо от того, являетесь ли вы ученым, инженером или кем-то еще, поиск корней — это проблема, с которой вы столкнетесь по пути, занимаясь числовыми задачами. Неважно, пытаетесь ли вы решить физическую модель, проводите какие-то операционные исследования для распределения ресурсов, планируете маршрут или выполняете подгонку данных, поиск корней будет в любом случае.
Python сам по себе медленный, но кто вообще решает проблемы на чистом Python? Если вы здесь, то наверняка знаете, что Python — это просто интерфейс для низкоуровневого скомпилированного кода, который выполняется достаточно быстро для обработки алгоритмов, интенсивно использующих ЦП. В противном случае вы забыли сделать какое-то домашнее задание!
Далее мы увидим, как решить простую задачу поиска корня с помощью различных пакетов Python jungle. По ходу обсуждаются преимущества и недостатки каждого из подходов. Идея состоит в том, чтобы сохранить проблему простой, но достаточно репрезентативной для общего нелинейного приложения для поиска корней.
Постановка задачи
В настоящем исследовании мы реализуем решение для стационарного состояния Брюсселятора при известном стабильном наборе параметров. Если после прочтения вы захотите попрактиковаться в других подобных задачах, вы можете проверить уравнения Лотки-Вольтерры или Орегонатор. Я позволяю вам проверить ссылку на Википедию для математической формулировки уравнений.
Примечание: в следующем посте мы будем использовать тот же набор инструментов для динамического моделирования Brusselator, так что подождите, потому что сегодня наша цель — понять только реализации поиска корней.
Поскольку мы хотим совместно использовать некоторые параметры для различных реализаций, которые будут тестироваться, не беспокоясь об изменении значений во всех файлах или чтении какого-либо текстового/JSON-файла, мы выгружаем приведенный ниже класс Params в params .py файл. Таким образом, мы можем извлечь выгоду из системы импорта Python, чтобы легко загружать эти значения.
# -*- coding: utf-8 -*-
class Params:
""" Standard problem parameters. """
# Initial guess for problem.
x0 = [1.0, 1.0]
# Brusselator coefficients
A = 1.0
B = 1.7
Подходы к решению
Пришло время приложить наши руки к работе! Сегодня мы рассмотрим три способа поиска корней Brusselator. Прежде чем мы введем код, несколько слов о пакетах, которые мы будем использовать.
Начнем с SciPy. Если вы уже применяли некоторые численные методы в Python, вы наверняка знакомы с этим пакетом. Это один из основных игроков в экосистеме Python, и он имеет широкий спектр охватываемых алгоритмов. Там вы найдете оптимизацию, интерполяцию, интеграцию, обработку сигналов и многое другое. Хотя он предоставляет возможности некоторой ограниченной оптимизации, это довольно громоздкий интерфейс, и я бы придерживался его только в простых случаях. Здесь мы используем версию 1.10.1.
Второй пакет, который мы рассмотрим, называется CasADi. Если вы не разбираетесь в средствах управления прогнозирующими моделями, вы никогда об этом не слышали. Этот пакет выполняет работу по вычислению якобианской матрицы вашей задачи за вас. Почему это здорово? Нелинейный поиск корня основан на своего рода градиентном спуске, и подпрограммы должны вычислять производные конечных разностей, чтобы найти направление прогресса. Это ОЧЕНЬ подвержено ошибкам, и следование градиенту, вычисляемому на основе аналитического якобиана, обеспечивает большую надежность. Да, вы можете предоставить SciPy проблемный якобиан, но удачи в его получении и реализации для сложных случаев. Здесь использовалась версия 3.5.6.
В завершение мы вводим Pyomo. У этого пакета более широкая аудитория, чем у CasADi, но и более крутая кривая обучения. Pyomo делает вещи формальными с математической точки зрения. Позвольте мне уточнить. Pyomo позволяет вам объявлять переменные и параметры с заданными доменами, границами, инициализацией и т. д., используя правильную терминологию. Кроме того, он поддерживает ограничения неравенства, что сложнее сделать с CasADi. Если вы занимаетесь операционными исследованиями, это то, что вам нужно. Мы запускали скрипт с версией 6.4.1.
В качестве альтернативы вы также можете использовать Gekko, но я воздержусь от его изучения, потому что исполняемый файл APM, распространяемый вместе с ним, не имеет открытого исходного кода.
Известный способ SciPy
Реализация SciPy является самой простой. Вам нужно сформулировать уравнение невязки задачи и вернуть массив, все элементы которого оцениваются как ноль, когда решение найдено. Для решения проблемы используется подпрограмма root из scipy.optimize. Начальное предположение, а также другие параметры функции предоставляются через его интерфейс. Полный набор параметров приведен в его документации. Я уже использовал эту функцию для решения задач с несколькими тысячами уравнений, с лучшими временами выполнения стены, используя крылов и df-sane методы. По умолчанию будет использоваться решатель hybr.
# -*- coding: utf-8 -*-
from scipy.optimize import root
from params import Params
def brusselator(x, a, b):
""" Steady-state Brusselator residual equation. """
eq1 = x[0] * (b - x[0] * x[1])
eq0 = a - x[0] - eq1
return [eq0, eq1]
def solve(brusselator, x0=Params.x0, a=Params.A, b=Params.B):
""" Solve problem with given coefficients and guess. """
return root(brusselator, x0, args=(a, b))
sol = solve(brusselator)
print(sol)
Объект, возвращаемый методом, содержит флаг, указывающий на успех, некоторые важные значения и решение x проблемы.
In [1]: %run gist-scipy.py
message: The solution converged.
success: True
status: 1
fun: [ 4.441e-16 -4.441e-16]
x: [ 1.000e+00 1.700e+00]
nfev: 5
fjac: [[-9.191e-01 -3.939e-01]
[ 3.939e-01 -9.191e-01]]
r: [ 7.616e-01 -5.252e-01 1.313e+00]
qtf: [-8.019e-13 2.004e-12]
Универсальный CasADi
CasADi предоставляет платформу SX для построения графа символьных операций. Это позволяет пакету выполнять свою основную задачу: автоматическое дифференцирование. Вот как ему удается предоставить решателю аналитический якобиан. Интересной особенностью CasADi является то, что вы можете сформулировать задачу и передать ее различным нелинейным решателям, даже коммерческим, используя один и тот же интерфейс. Он поставляется с отличным решателем с открытым исходным кодом Ipopt. Поскольку наша задача состоит только из уравнений ограничения, мы устанавливаем цель f равной единице и предоставляем массив ограничений g . Обратите внимание, что граф операций принимает свободные параметры через массив p в нелинейной постановке задачи, которые могут быть предоставлены только во время решения. Наконец, при вызове решателя мы можем указать допустимые отклонения для ограничений с помощью lbg и ubg, что может быть массивом того же размера, что и g, или просто скаляром. Установка обоих значений на ноль здесь обеспечивает цель поиска корня.
# -*- coding: utf-8 -*-
from pprint import pprint
from casadi import SX
from casadi import nlpsol
from casadi import vertcat
from params import Params
def create_solver():
""" Create steady-state Brusselator solver. """
p = SX.sym("p", 2)
x = SX.sym("x", 2)
a, b = p[0], p[1]
eq1 = x[0] * (b - x[0] * x[1])
eq0 = a - x[0] - eq1
g = vertcat(eq0, eq1)
nlp = {"x": x, "p": p, "f": 1.0, "g": g}
opts = {"ipopt.print_level": 5}
return nlpsol("brusselator", "ipopt", nlp, opts)
def solve(brusselator, x0=Params.x0, a=Params.A, b=Params.B):
""" Solve problem with given coefficients and guess. """
return brusselator(x0=x0, p=[a, b], lbg=0.0, ubg=0.0)
brusselator = create_solver()
sol = solve(brusselator)
pprint(sol)
При решении проблемы с помощью Ipopt решатель отобразит сводку проблемы, проследит процесс решения и подытожит решение. Уровень детализации можно контролировать при создании решателя, см. переменную opts выше в create_solver. Чтобы получить статус решателя, вы можете использовать brusselator.stats()["return_status"], чего нет в примере кода. Лично я использовал этот пакет в ситуациях, когда чисто численный подход SciPy терпел неудачу, и это мой выбор при разработке управления с прогнозированием моделей (MPC) или любой другой оптимизации с ограничениями. При определенных обстоятельствах вы также можете сгенерировать код C для развертывания решателя.
In [1]: %run gist-casadi.py
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.4, running with linear solver MUMPS 5.4.1.
Number of nonzeros in equality constraint Jacobian...: 4
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 2
Total number of variables............................: 2
variables with only lower bounds: 0
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 2
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 1.0000000e+00 7.00e-01 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 1.0000000e+00 0.00e+00 0.00e+00 -1.7 7.00e-01 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 1
(scaled) (unscaled)
Objective...............: 1.0000000000000000e+00 1.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 0.0000000000000000e+00 0.0000000000000000e+00
Number of objective function evaluations = 2
Number of objective gradient evaluations = 2
Number of equality constraint evaluations = 2
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 1
Total seconds in IPOPT = 0.019
EXIT: Optimal Solution Found.
brusselator : t_proc (avg) t_wall (avg) n_eval
nlp_f | 0 ( 0) 3.00us ( 1.50us) 2
nlp_g | 0 ( 0) 8.00us ( 4.00us) 2
nlp_grad | 0 ( 0) 4.00us ( 4.00us) 1
nlp_grad_f | 0 ( 0) 12.00us ( 4.00us) 3
nlp_hess_l | 0 ( 0) 1.00us ( 1.00us) 1
nlp_jac_g | 0 ( 0) 4.00us ( 1.33us) 3
total | 34.00ms ( 34.00ms) 30.02ms ( 30.02ms) 1
{'f': DM(1),
'g': DM([0, 0]),
'lam_g': DM([0, -0]),
'lam_p': DM([-0, 0]),
'lam_x': DM([0, 0]),
'x': DM([1, 1.7])}
Плотный матричный объект DM, представленный в приведенном выше решении, можно преобразовать в массив NumPy для дальнейшей обработки в Python.
Официально с Pyomo
Я говорил вам, что путь изучения Pyomo был долгим, но это окупается. Что касается CasADi, Pyomo также принимает несколько решателей, в том числе проприетарные. Учитывая его более форматный характер, именно во время объявления мы сообщаем модели, является ли что-то параметром или переменной, его математической областью и так далее. При создании Ограничения вы предоставляете правило, которое может быть неравенством, что сразу не поддерживается CasADi. . Объем кода больше, чем в предыдущих пакетах для той же проблемы, но цели яснее. Pyomo — мой предпочтительный пакет при решении задач, не связанных с ступенчатыми функциями (я все еще ищу надежный способ использования функции Хевисайда в Pyomo) или для моделирования смешанных целых чисел, функция, которая не поддерживается другими обсуждаемыми альтернативами. . Он также отлично подходит для решения проблем с DAE, но это тема для другого поста.
# -*- coding: utf-8 -*-
from pathlib import Path
from pprint import pprint
import pyomo.environ as pe
import pyomo.opt as po
from params import Params
apps = Path.home() / "Applications"
ipopt = apps / "Ipopt-3.14.12-win64-msvs2019-md/bin/ipopt.exe"
def create_solver():
""" Create steady-state Brusselator solver. """
model = pe.ConcreteModel(name="brusselator")
model.i = pe.Set(initialize=[0, 1], dimen=1)
model.a = pe.Param(mutable=True)
model.b = pe.Param(mutable=True)
model.x = pe.Var(model.i)
eq1 = model.x[0] * (model.b - model.x[0] * model.x[1])
eq0 = model.a - model.x[0] - eq1
model.eq0 = pe.Constraint(rule=(0.0 == eq0))
model.eq1 = pe.Constraint(rule=(0.0 == eq1))
return model
def solve(brusselator, x0=Params.x0, a=Params.A, b=Params.B):
""" Solve problem with given coefficients and guess. """
brusselator.a.set_value(a)
brusselator.b.set_value(b)
for k, xk in enumerate(x0):
brusselator.x[k] = xk
solver = po.SolverFactory("ipopt", executable=ipopt)
return solver.solve(brusselator)
brusselator = create_solver()
sol = solve(brusselator)
sol.write()
x = [xk.value for xk in brusselator.x.values()]
pprint(x)
После этого распечатывается статистика решателя, что является преимуществом по сравнению с CasADi, поскольку экономит много времени на обмен данными со стандартным выводом.
In [1]: %run gist-pyomo.py # ========================================================== # = Solver Results = # ========================================================== # ---------------------------------------------------------- # Problem Information # ---------------------------------------------------------- Problem: - Lower bound: -inf Upper bound: inf Number of objectives: 1 Number of constraints: 2 Number of variables: 2 Sense: unknown # ---------------------------------------------------------- # Solver Information # ---------------------------------------------------------- Solver: - Status: ok Message: Ipopt 3.14.12\x3a Optimal Solution Found Termination condition: optimal Id: 0 Error rc: 0 Time: 0.08400368690490723 # ---------------------------------------------------------- # Solution Information # ---------------------------------------------------------- Solution: - number of solutions: 0 number of solutions displayed: 0 [1.0, 1.7]
Заключительные мысли
Я надеюсь, что этот пост помог вам сделать осознанный выбор реализации ваших нелинейных задач поиска корней.
Как мы показали, ни один из исследованных пакетов не поддерживает все функции, которые могут потребоваться, и это необходимо учитывать при запуске нового проекта. Еще одной важной особенностью является размер сообщества и команда разработчиков пакетов. И CasADi, и Pyomo являются нишевыми приложениями, и найти примеры на Stack Overflow может быть проблемой.
Понравился мой контент? Пожалуйста, помогите мне производить больше, подписавшись на мой Патреон!
