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

Открыть пример в Engee

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

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

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

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

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

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

• одно вещество $C_1$ и одно вещество $C_2$ реагируют, образовывая одно вещество $C_3$ со скоростью $r_1$;
• одно вещество $C_3$ и одно вещество $C_4$ реагируют, образовывая одно вещество $C_5$ со скоростью $r_2$;
• одно вещество $C_3$ и одно вещество $C_4$ реагируют, образовывая одно вещество $C_6$ со скоростью $r_3$.

Скорость реакции пропорциональна произведению количеств требуемых веществ. Таким образом, если $y_i$ представляет количество вещества $C_i$, то скорость реакции для получения $C_3$ равна $r_1y_1y_2$. Аналогично, скорость реакции для получения $C_5$ равна $r_2y_3y_4$, а скорость реакции для получения $C_6$ равна $r_3y_3y_4$.

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

$$ \frac{dy}{dt} = \begin{bmatrix} -r_1y_1y_2 \\ -r_1y_1y_2 \\ -r_2y_3y_4 + r_1y_1y_2 - r_3y_3y_4 \\ -r_2y_3y_4 - r_3y_3y_4 \\ r_2y_3y_4 \\ r_3y_3y_4 \\ \end{bmatrix} $$

Зададим начальные значения $y(0)=[1,1,0,1,0,0]$. Такие начальные значения гарантируют, что все вещества полностью прореагируют, в результате чего $C_1–C_4$ будут приближаться к нулю по мере увеличения времени.

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

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

#Pkg.add("JuMP");

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

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

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

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

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

using JuMP;

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

using Ipopt;

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

using DifferentialEquations;

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

using Plots;

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

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

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

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

function check_bounds(min::Number, max::Number)
    if min > max
        error("Значение min должно быть меньше чем значение max")
    elseif max < min
        error("Значение max должно быть больше чем значение min")
    end
end

check_bounds (generic function with 1 method)

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

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

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

dif_fun (generic function with 1 method)

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

r_true = [2.5, 1.2, 0.45];

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

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

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

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

prob = ODEProblem(dif_fun, y_0, t_span, r_true);

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

sol_true = solve(prob, Tsit5(), saveat=0.1);

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

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

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, рассчитывающую разницу между истинными параметрами $r_1$, $r_2$ и $r_3$ и рассчитанными значениями.

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

obj_fun (generic function with 1 method)

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

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

model = Model(Ipopt.Optimizer)

A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Ipopt

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

set_silent(model)

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

@variable(model, 0.1 <= r[1:3] <= 10)

3-element Vector{VariableRef}:
 r[1]
 r[2]
 r[3]

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

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

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

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

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

optimize!(model);

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

optimized_params = value.(r)

3-element Vector{Float64}:
 2.5000119303476147
 1.200000509788957
 0.4499998880470981

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

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]

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

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

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

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

obj_fun_lim (generic function with 1 method)

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

model_lim = Model(Ipopt.Optimizer)

A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Ipopt

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

set_silent(model_lim)

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

@variable(model_lim, 0.1 <= r_lim[1:3] <= 10)

3-element Vector{VariableRef}:
 r_lim[1]
 r_lim[2]
 r_lim[3]

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

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

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

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

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

optimize!(model_lim)

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

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

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

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

prob_opt = ODEProblem(dif_fun, y_0, t_span, value.(r_lim));

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

sol_opt = solve(prob_opt, Tsit5(), saveat=0.1);

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

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

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

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

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

Заключение

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