寻找最佳 ODE 参数¶

本例演示了为描述化学反应链的常微分方程（ODE）寻找最佳参数的工作流程。我们将找到参数，并评估这些参数能在多大程度上再现系统的真实动态。

该问题使用基于问题的优化工作流程解决。

在本例中，我们将使用 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(["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));

描述一个定义$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 ，用于描述化学反应链的 ODE：

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_true 中的$r_1$ 、$r_2$ 和$r_3$ 值是优化算法应尝试重新创建的基准值。

为$С_1 - С_6$ 设置响应模拟持续时间和初始值：

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

使用函数ODEProblem 创建一个 ODE 问题：

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() ，该求解器推荐适用于大多数非刚性问题。

将获得的真解结果绘制成图：

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 ，因此必须事先在函数register 中注册，以便正确集成到库JuMP 的优化环境中。在括号中指定输入参数的数量 ($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 创建一个 ODE 问题：

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

求解 ODE 问题：

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$ 的动态完全相同。

结论¶

在本例中，我们使用[基于问题的优化]工作流 (https://engee.com/community/catalogs/projects/algoritm-resheniia-zadach-optimizatsii) 解决了寻找描述多阶段化学反应链的 ODE 的最优参数的问题。我们模拟了有限观测条件下的情况，并比较了系统的动态变化。