Документация Engee
Notebook

Поиск оптимальных параметров ОДУ

В данном примере продемонстрирован рабочий процесс поиска оптимальных параметров для обыкновенного дифференциального уравнения (ОДУ), описывающего цепь химических реакций. Мы найдем параметры и оценим, насколько близко они могут воспроизводить истинную динамику системы.

Данная задача решена с помощью рабочего процесса проблемно-ориентированной оптимизации.

В этом примере мы воспользуемся функциями библиотеки JuMP.jl для формулировки оптимизационной задачи, библиотекой нелинейной оптимизации Ipopt.jl, библиотекой для решения дифференциальных уравнений DifferentialEquations.jl и библиотекой Plots.jl для визуализации результатов.

Описание задачи

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

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

  • одно вещество и одно вещество реагируют, образовывая одно вещество со скоростью ;
  • одно вещество и одно вещество реагируют, образовывая одно вещество со скоростью ;
  • одно вещество и одно вещество реагируют, образовывая одно вещество со скоростью .
Screenshot_20241030_112904.png

Скорость реакции пропорциональна произведению количеств требуемых веществ. Таким образом, если представляет количество вещества , то скорость реакции для получения равна . Аналогично, скорость реакции для получения равна , а скорость реакции для получения равна .

Таким образом, мы можем составить дифференциальное уравнение, описывающее процесс изменения нашей системы:

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

Установка библиотек

Если в вашем окружении не установлена последняя версия пакета JuMP, раскомментируйте и запустите ячейку ниже:

In [ ]: 
Pkg.add(["Ipopt", "DifferentialEquations", "JuMP"])
In [ ]: 
#Pkg.add("JuMP");

Для запуска новой версии библиотеки после завершения установки нажмите на кнопку «Личный кабинет»:

Screenshot_20240710_134852.png Затем нажмите на кнопку «Стоп»: Screenshot_20240710_2.png

Перезапустите сессию нажатием кнопки «Старт Engee»:

Screenshot_20240710_135431.png

Подключение библиотек

Подключите библиотеку JuMP:

In [ ]: 
using JuMP;

Подключите библиотеку нелинейного решателя Ipopt:

In [ ]: 
using Ipopt;

Подключите библиотеку для решения дифференциальных уравнений DifferentialEquations:

In [ ]: 
using DifferentialEquations;

Подключите библиотеку построения графиков Plots:

In [ ]: 
using Plots;

Конфигурация параметров визуализации

Задайте кастомизированный размер графиков и бэкграунд plotly:

In [ ]: 
plotly(size=(1500, 600), xlims=(0, :auto));

Опишите функцию, определяющую границы для и осей для приближения графиков:

In [ ]: 
function check_bounds(min::Number, max::Number)
    if min > max
        error("Значение min должно быть меньше чем значение max")
    elseif max < min
        error("Значение max должно быть больше чем значение min")
    end
end
Out[0]:
check_bounds (generic function with 1 method)

Решение дифференциальных уравнений

Создайте функцию dif_fun, описывающую ОДУ нашей цепи химических реакций:

In [ ]: 
function dif_fun(du, u, p, t)
    r1, r2, r3 = p

    du[1] = -r1 * u[1] * u[2]
    du[2] = -r1 * u[1] * u[2]
    du[3] = -r2 * u[3] * u[4] + r1 * u[1] * u[2]- r3 * u[3] * u[4]
    du[4] = -r2 * u[3] * u[4] - r3 * u[3] * u[4]
    du[5] = r2 * u[3] * u[4]
    du[6] = r3 * u[3] * u[4]
end
Out[0]:
dif_fun (generic function with 1 method)

Задайте истинные значения и сохраните их в переменной r_true:

In [ ]: 
r_true = [2.5, 1.2, 0.45];

Значения , и в переменной r_true являются эталоном, который алгоритм оптимизации должен попытаться воссоздать.

Задайте длительность симуляции реакции и начальные значения для :

In [ ]: 
t_span = (0.0, 5.0);
y_0 = [1.0, 1.0, 0.0, 1.0, 0.0, 0.0];

Создайте задачу ОДУ с помощью функции ODEProblem:

In [ ]: 
prob = ODEProblem(dif_fun, y_0, t_span, r_true);

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

In [ ]: 
sol_true = solve(prob, Tsit5(), saveat=0.1);

Для решения данной задачи используем решатель Рунге-Кутты Ципураса 5 порядка Tsit5(), который является рекомендованным для большинства нежестких задач.

Выведите полученные результаты истинных решений на графики:

In [ ]: 
colors = [:red, :blue, :green, :orange, :purple, :brown]

plot(layout=(3,2))
for i in 1:6
    plot!(subplot=i, sol_true.t, sol_true[i,:], color=colors[i], title="y($i)", label="Истинное y($i)")
end
display(current())

Решение задачи оптимизации

Опишите целевую функцию obj_fun, рассчитывающую разницу между истинными параметрами , и и рассчитанными значениями.

In [ ]: 
function obj_fun(r1, r2, r3)
    r = [r1, r2, r3]
    prob_new = ODEProblem(dif_fun, y_0, t_span, r)
    sol_new = solve(prob_new, Tsit5(), saveat=0.1)
    return sum(sum((sol_new[i,:] .- sol_true[i,:]).^2) for i in 1:6)
end
Out[0]:
obj_fun (generic function with 1 method)

Для нахождения оптимальных решений нам нужно будет минимизировать эту функцию.

Создайте оптимизационную задачу с помощью функции Model и укажите название решателя в скобках:

In [ ]: 
model = Model(Ipopt.Optimizer)
Out[0]:
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Ipopt

Переведите модель в «молчаливый» режим, чтобы скрыть системные сообщения:

In [ ]: 
set_silent(model)

Создайте переменную r, содержащую три элемента – , и , и задайте границы от 0.1 до 10:

In [ ]: 
@variable(model, 0.1 <= r[1:3] <= 10)
Out[0]:
3-element Vector{VariableRef}:
 r[1]
 r[2]
 r[3]

Так как мы имеем дело с комплексной пользовательской нелинейной функцией obj_fun, для корректной интеграции в оптимизационную среду библиотеки JuMP её предварительно нужно зарегистрировать с помощью функции register. В скобках укажите количество входящих аргументов (, и ) и включите автоматическую дифференциацию.

In [ ]: 
register(model, :obj_fun, 3, obj_fun; autodiff = true)

Создайте нелинейную целевую функцию для нахождения минимальных значений , и и примените ранее созданную функцию obj_fun:

In [ ]: 
@NLobjective(model, Min, obj_fun(r[1], r[2], r[3]))

Решите оптимизационную задачу:

In [ ]: 
optimize!(model);

Сохраните оптимизированные параметры в переменной:

In [ ]: 
optimized_params = value.(r)
Out[0]:
3-element Vector{Float64}:
 2.5000119303476147
 1.200000509788957
 0.4499998880470981

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

In [ ]: 
optimized_params = round.(value.(r), digits=2)
println("Истинные параметры: ", r_true)
println("Оптимизированные параметры: ", optimized_params)

Истинные параметры: [2.5, 1.2, 0.45]
Оптимизированные параметры: [2.5, 1.2, 0.45]

Оптимизация при ограниченных наблюдениях

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

Опишите целевую функцию obj_fun_lim, являющуюся копией функции obj_fun из предыдущего раздела, но возвращающую только результаты для и .

In [ ]: 
function obj_fun_lim(r1, r2, r3)
    r = [r1, r2, r3]
    prob_new = ODEProblem(dif_fun, y_0, t_span, r)
    sol_new = solve(prob_new, Tsit5(), saveat=0.1)
    return sum(sum((sol_new[i,:] .- sol_true[i,:]).^2) for i in 5:6)
end
Out[0]:
obj_fun_lim (generic function with 1 method)

Создайте оптимизационную задачу с помощью функции Model и укажите название решателя в скобках:

In [ ]: 
model_lim = Model(Ipopt.Optimizer)
Out[0]:
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Ipopt

Переведите модель в «молчаливый» режим, чтобы скрыть системные сообщения:

In [ ]: 
set_silent(model_lim)

Создайте переменную r, содержащую три элемента – , и , и задайте границы от 0.1 до 10:

In [ ]: 
@variable(model_lim, 0.1 <= r_lim[1:3] <= 10)
Out[0]:
3-element Vector{VariableRef}:
 r_lim[1]
 r_lim[2]
 r_lim[3]

Зарегистрируйте функцию obj_fun_lim с помощью функции register:

In [ ]: 
register(model_lim, :obj_fun_lim, 3, obj_fun_lim; autodiff = true)

Создайте нелинейную целевую функцию для нахождения минимальных значений , и и примените ранее созданную функцию obj_fun_lim:

In [ ]: 
@NLobjective(model_lim, Min, obj_fun_lim(r[1], r[2], r[3]))

Решите оптимизационную задачу:

In [ ]: 
optimize!(model_lim)

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

In [ ]: 
println("Истинные параметры: ", r_true)
println("Оптимизированные параметры (ограниченные): ", value.(r_lim))

Истинные параметры: [2.5, 1.2, 0.45]
Оптимизированные параметры (ограниченные): [1.8085855441520966, 1.5505236127674895, 0.5814611991591475]

Создайте задачу ОДУ с помощью функции ODEProblem, применив оптимизированные значения r_lim:

In [ ]: 
prob_opt = ODEProblem(dif_fun, y_0, t_span, value.(r_lim));

Решите задачу ОДУ:

In [ ]: 
sol_opt = solve(prob_opt, Tsit5(), saveat=0.1);

Визуализируйте результаты для и :

In [ ]: 
y_min = 0 # @param {type:"slider", min:0, max:1, step:0.01}
y_max = 0.6 # @param {type:"slider", min:0, max:1, step:0.01}
check_bounds(y_min, y_max)
x_min = 0 # @param {type:"slider", min:0, max:5, step:0.1}
x_max = 5 # @param {type:"slider", min:0, max:5, step:0.1}
check_bounds(x_min, x_max)

plot(sol_true.t, [sol_true[5,:] sol_true[6,:]], label=["Истинное y(5)" "Истинное y(6)"], ylimits=(y_min, y_max), xlimits=(x_min, x_max))
plot!(sol_true.t, [sol_opt[5,:] sol_opt[6,:]], label=["Оптимизированное y(5)" "Оптимизированное y(6)"], linestyle=:dash, ylimits=(y_min, y_max), xlimits=(x_min, x_max))
display(current())

Визуализируйте и сравните результаты ранее заданных истинных значений и полученных значений :

In [ ]: 
true_colors = [:darkred, :darkblue, :darkgreen, :darkorange, :darkviolet, :brown]
opt_colors = [:red, :blue, :green, :orange, :violet, :orange]

plot(layout=(3,2))
for i in 1:6
    plot!(subplot=i, sol_true.t, sol_true[i,:], color=true_colors[i], label="Истинное y($i)", title="y($i)")
    plot!(subplot=i, sol_opt.t, sol_opt[i,:], color=opt_colors[i], label="Оптимизированное y($i)", linestyle=:dash)
end
display(current())

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

Заключение

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