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

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

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

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

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

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

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

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

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

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

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

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

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

Pkg.add(["Ipopt", "DifferentialEquations", "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));

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

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_true:

r_true = [2.5, 1.2, 0.45];

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

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

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

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, содержащую три элемента – , и , и задайте границы от 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. В скобках укажите количество входящих аргументов (, и ) и включите автоматическую дифференциацию.

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

Создайте нелинейную целевую функцию для нахождения минимальных значений , и и примените ранее созданную функцию 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]

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

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

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

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, содержащую три элемента – , и , и задайте границы от 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)

Создайте нелинейную целевую функцию для нахождения минимальных значений , и и примените ранее созданную функцию 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_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())

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

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

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

Заключение

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

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

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

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

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

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

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

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

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

Заключение

Тип

Краткое описание

Связанные материалы

Теги

Категории

Языки

Форматы

Уровень

Статус

Опубликовано

Обновлено

Источник