搜索最佳ODE参数

此示例演示了为描述化学反应链的常微分方程(ODE)查找最佳参数的工作流程。我们将找到这些参数，并评估它们如何密切地再现系统的真实动态。

使用[面向问题的优化]工作流解决了此任务（https://engee.com/community/catalogs/projects/algoritm-resheniia-zadach-optimizatsii )。

在这个例子中，我们将使用[跳转的功能。jl]图书馆(https://engee.com/helpcenter/stable/ru/julia/JuMP/index.html ）使用非线性优化库[Ipopt.jl]来制定优化问题（https://engee.com/helpcenter/stable/ru/julia/Ipopt/README.html ），求解微分方程的库[DifferentialEquations.jl]（https://engee.com/helpcenter/stable/ru/julia/DifferentialEquations/index.html ）和[情节。jl]图书馆(https://engee.com/helpcenter/stable/ru/julia/juliaplots/index.html ）来可视化结果。

任务说明

让我们考虑一个涉及几种物质的多阶段化学反应链的模型。 . 物质以一定的速率相互反应它们形成新的物质。此任务的真实反应速率未知。但是我们可以在给定的时间内进行一定次数的观察。 . 根据这些观察结果，我们将选择最佳反应速率。 .

模型中有六种物质，从以前其反应如下:

*一种物质和一种物质它们反应形成单一物质速度快 ;
*一种物质和一种物质它们反应形成单一物质速度快 ;
*一种物质和一种物质它们反应形成单一物质速度快 .

反应速率与所需物质的量的乘积成比例。因此，如果表示物质的量，然后反应速率得到等于 . 同样，反应速率得到等于，和反应速率得到等于 .

因此，我们可以创建一个微分方程，描述改变我们系统的过程。:

设置初始值 . 这样的初始值确保所有物质充分反应，导致随着时间的增加，它们将接近零。

安装库

如果您的环境中未安装最新版本的软件包 JuMP，取消注释并运行下面的单元格:

Pkg.add(["Ipopt", "DifferentialEquations", "JuMP"])

#Pkg.add("JuMP");

要在安装完成后启动新版本的库，请单击"个人帐户"按钮:

然后点击"停止"按钮:

通过单击"Start 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];

使用函数创建ODE任务 ODEProblem:

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

使用函数解决问题 solve. 在括号中，指定求解器的名称和参数中观测值的频率。 saveat:

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

为了解决这个问题，我们使用5阶的Runge-Kutta Tsipouras求解器。 Tsit5()，即[рекомендованным]（https://engee.com/helpcenter/stable/ru/julia/DifferentialEquations/solvers/ode_solve.html ）对于大多数非刚性任务。

在图形上显示真实解的结果:

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]

使用函数创建ODE任务 ODEProblem 通过应用优化值 r_lim:

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

解决ODE问题:

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

我们可以发现，随着有限观测得到的新的优化参数，物质和它们被消耗得更慢，而物质较小程度上累积。然而，物质 , 和对于已找到的参数和真实参数，它们具有完全相同的动态。

结论

在这个例子中，我们解决了使用工作流[面向问题的优化]找到描述多阶段化学反应链的ODE的最佳参数的问题（https://engee.com/community/catalogs/projects/algoritm-resheniia-zadach-optimizatsii )。我们用有限的观察对情况进行了建模，并比较了系统的动力学。