Оценка параметров нескольких наборов данных в julia DifferentialEquations

Я искал и не смог найти прямой способ использования оценки параметра DifferentialEquations в julia для соответствия нескольким наборам данных. Итак, допустим, у нас есть это простое дифференциальное уравнение с двумя параметрами:

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

У нас есть экспериментальные наборы данных u [1] vs 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] (который представляет некоторую экспериментальную переменную (например, температуру, расход ...)

Цель состоит в том, чтобы оценить значение p [1], используя экспериментальные данные sol1, sol2 и sol3, позволяя DiffEqBase.problem_new_parameters или другой альтернативе перебирать экспериментальные условия.


julia least-squares mathematical-optimization differentialequations.jl
person JIbab    schedule 13.06.2018    source источник
comment
Поделитесь некоторым примером набора данных   -  person DeshDeep Singh    schedule 13.06.2018

Ответы (1)


arrow_upward
3
arrow_downward

Что вы можете сделать, так это создать задачу Монте-Карло, которая решает сразу все три проблемы:

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

Наконец, вам нужно сказать ему, как связать параметр (ы) оптимизации с проблемой, которую он решает. Здесь наша задача Монте-Карло содержит 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)

Теперь перебросим это в метод Brent Optim.jl:

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)
person Chris Rackauckas    schedule 16.06.2018
comment
Спасибо, это то, что я искал! - person JIbab; 17.06.2018
comment
Возвращаясь к этому вопросу, у меня есть любопытство. Я тестировал, и использование saveat в целевой функции необходимо, когда это проблема Монтекарло, но не для обычной задачи ODE, я предполагаю, что она интерполирует sol (t [i]). При попытке оптимизировать задачу Монте-Карло без использования сохранения, проблема сводится к неправильному решению. Было бы полезно иметь хотя бы один из этих двух вариантов: 1. Укажите saveat для каждой итерации Монте-Карло (поскольку разные наборы данных могут иметь разные временные точки) 2. Использовать интерполяцию в sol [i] (t [j]) - person JIbab; 26.11.2018
comment
Вероятно, мы сможем справиться с этим на стороне библиотеки. Пожалуйста, откройте вопрос. - person Chris Rackauckas; 26.11.2018