julia - julia DifferentialEquations 中多个数据集的参数估计

Question

我一直在寻找，但找不到使用 julia 中的 DifferentialEquations 参数估计来拟合多个数据集的直接方法。因此，假设我们有一个带有两个参数的简单微分方程:

f1 = function (du,u,p,t)
  du[1] = - p[1]*p[2] * u[1]
end

我们有 u[1] 与 t 的实验数据集。每个数据集具有不同的 p[2] 值和/或不同的初始条件。 p[1] 是我们要估计的参数。我可以通过在 for 循环中求解微分方程来实现此目的，该循环迭代不同的初始条件和 p[2] 值，将解存储在数组中并根据实验数据创建损失函数。我想知道是否有一种方法可以用更少的代码行来完成此操作，例如使用DiffEqBase.problem_new_parameters来设置每个数据集的条件。在将模型拟合到实验数据时，这是一种非常常见的情况，但我在文档中找不到很好的示例。

提前谢谢您，

致以诚挚的问候

编辑 1

上面表达的情况只是一个简化的示例。为了使其成为一个实际案例，我们可以通过以下代码创建一些虚假的实验数据:

using DifferentialEquations

# ODE function
f1 = function (du,u,p,t)
  du[1] = - p[1]*p[2] * u[1]
end

# Initial conditions and parameter values.
# p1 is the parameter to be estimated.
# p2 and u0 are experimental parameters known for each dataset.
u0 = [1.,2.]
p1 = .3
p2 = [.1,.2]
tspan = (0.,10.)
t=linspace(0,10,5)


# Creating data for 1st experimental dataset.
## Experimental conditions: u0 = 1. ; p[2] = .1
prob1 = ODEProblem(f1,[u0[1]],tspan,[p1,p2[1]])
sol1=solve(prob1,Tsit5(),saveat=t)

# Creating data for 2nd experimental dataset.
## Experimental conditions: u0 = 1. ; p[2] = .2
prob2 = ODEProblem(f1,[u0[1]],tspan,[p1,p2[2]])
sol2=solve(prob2,Tsit5(),saveat=t)

# Creating data for 3rd experimental dataset.
## Experimental conditions: u0 = 2. ; p[2] = .1
prob3 = ODEProblem(f1,[u0[2]],tspan,[p1,p2[1]])
sol3=solve(prob3,Tsit5(),saveat=t)

sol1、sol2 和 sol3 现在是我们的实验数据，每个数据集使用初始条件和 p[2](代表一些实验变量(例如温度、流量...)的不同组合

目标是使用实验数据 sol1、sol2 和 sol3 来估计 p[1] 的值，让 DiffEqBase.problem_new_parameters 或其他替代方案迭代实验条件。

Answer 1

您可以做的是创建一个 MonteCarloProblem 来同时解决所有三个问题:

function prob_func(prob,i,repeat)
  i < 3 ? u0 = [1.0] : u0 = [2.0]
  i == 2 ? p = (prob.p[1],0.2) : p = (prob.p[1],0.1)
  ODEProblem{true}(f1,u0,(0.0,10.0),p)
end
prob = MonteCarloProblem(prob1,prob_func = prob_func)

@time sol = solve(prob,Tsit5(),saveat=t,parallel_type = :none,
                  num_monte = 3)

然后创建一个损失函数，将每个解决方案与 3 个数据集进行比较，并将它们的损失加在一起。

loss1 = L2Loss(t,data1)
loss2 = L2Loss(t,data2)
loss3 = L2Loss(t,data3)
loss(sol) = loss1(sol[1]) + loss2(sol[2]) + loss3(sol[3])

最后，您需要告诉它如何将优化参数与其正在解决的问题联系起来。在这里，我们的 MonteCarloProblem 持有一个 prob，每当它生成问题时，它就会从中提取 p[1]。我们要优化的值是p[1]，所以:

function my_problem_new_parameters(prob,p)
  prob.prob.p[1] = p[1]
  prob
end

现在我们的目标正是将这些部分组合在一起:

obj = build_loss_objective(prob,Tsit5(),loss,
                           prob_generator = my_problem_new_parameters,
                           num_monte = 3,
                           saveat = t)

现在让我们将其交给 Optim.jl 的 Brent 方法:

using Optim
res = optimize(obj,0.0,1.0)

Results of Optimization Algorithm
 * Algorithm: Brent's Method
 * Search Interval: [0.000000, 1.000000]
 * Minimizer: 3.000000e-01
 * Minimum: 2.004680e-20
 * Iterations: 10
 * Convergence: max(|x - x_upper|, |x - x_lower|) <= 2*(1.5e-08*|x|+2.2e-16): true
 * Objective Function Calls: 11

发现总体最佳值为0.3，这是我们用来生成数据的参数。

这是完整的代码:

using DifferentialEquations

# ODE function
f1 = function (du,u,p,t)
  du[1] = - p[1]*p[2] * u[1]
end

# Initial conditions and parameter values.
# p1 is the parameter to be estimated.
# p2 and u0 are experimental parameters known for each dataset.
u0 = [1.,2.]
p1 = .3
p2 = [.1,.2]
tspan = (0.,10.)
t=linspace(0,10,5)

# Creating data for 1st experimental dataset.
## Experimental conditions: u0 = 1. ; p[2] = .1
prob1 = ODEProblem(f1,[u0[1]],tspan,[p1,p2[1]])
sol1=solve(prob1,Tsit5(),saveat=t)
data1 = Array(sol1)

# Creating data for 2nd experimental dataset.
## Experimental conditions: u0 = 1. ; p[2] = .2
prob2 = ODEProblem(f1,[u0[1]],tspan,[p1,p2[2]])
sol2=solve(prob2,Tsit5(),saveat=t)
data2 = Array(sol2)

# Creating data for 3rd experimental dataset.
## Experimental conditions: u0 = 2. ; p[2] = .1
prob3 = ODEProblem(f1,[u0[2]],tspan,[p1,p2[1]])
sol3=solve(prob3,Tsit5(),saveat=t)
data3 = Array(sol3)

function prob_func(prob,i,repeat)
  i < 3 ? u0 = [1.0] : u0 = [2.0]
  i == 2 ? p = (prob.p[1],0.2) : p = (prob.p[1],0.1)
  ODEProblem{true}(f1,u0,(0.0,10.0),p)
end
prob = MonteCarloProblem(prob1,prob_func = prob_func)

# Just to show what this looks like
sol = solve(prob,Tsit5(),saveat=t,parallel_type = :none,
            num_monte = 3)

loss1 = L2Loss(t,data1)
loss2 = L2Loss(t,data2)
loss3 = L2Loss(t,data3)
loss(sol) = loss1(sol[1]) + loss2(sol[2]) + loss3(sol[3])

function my_problem_new_parameters(prob,p)
  prob.prob.p[1] = p[1]
  prob
end
obj = build_loss_objective(prob,Tsit5(),loss,
                           prob_generator = my_problem_new_parameters,
                           num_monte = 3,
                           saveat = t)

using Optim
res = optimize(obj,0.0,1.0)