|
| 1 | +# [Optimizing a Parameterized ODE](@id optimizing-parameterized-ode) |
| 2 | + |
| 3 | +Let us fit a parameterized ODE to some data. We will use the Lotka-Volterra model as an |
| 4 | +example. We will use Single Shooting to fit the parameters. |
| 5 | + |
| 6 | +```@example parameterized_ode |
| 7 | +using OrdinaryDiffEq, NonlinearSolve, Plots |
| 8 | +``` |
| 9 | + |
| 10 | +Let us simulate some real data from the Lotka-Volterra model. |
| 11 | + |
| 12 | +```@example parameterized_ode |
| 13 | +function lotka_volterra!(du, u, p, t) |
| 14 | + x, y = u |
| 15 | + α, β, δ, γ = p |
| 16 | + du[1] = dx = α * x - β * x * y |
| 17 | + du[2] = dy = -δ * y + γ * x * y |
| 18 | +end |
| 19 | +
|
| 20 | +# Initial condition |
| 21 | +u0 = [1.0, 1.0] |
| 22 | +
|
| 23 | +# Simulation interval and intermediary points |
| 24 | +tspan = (0.0, 2.0) |
| 25 | +tsteps = 0.0:0.1:10.0 |
| 26 | +
|
| 27 | +# LV equation parameter. p = [α, β, δ, γ] |
| 28 | +p = [1.5, 1.0, 3.0, 1.0] |
| 29 | +
|
| 30 | +# Setup the ODE problem, then solve |
| 31 | +prob = ODEProblem(lotka_volterra!, u0, tspan, p) |
| 32 | +sol = solve(prob, Tsit5(); saveat = tsteps) |
| 33 | +
|
| 34 | +# Plot the solution |
| 35 | +using Plots |
| 36 | +plot(sol; linewidth = 3) |
| 37 | +savefig("LV_ode.png") |
| 38 | +``` |
| 39 | + |
| 40 | +Let us now formulate the parameter estimation as a Nonlinear Least Squares Problem. |
| 41 | + |
| 42 | +```@example parameterized_ode |
| 43 | +function loss_function(ode_param, data) |
| 44 | + sol = solve(prob, Tsit5(); p = ode_param, saveat = tsteps) |
| 45 | + return vec(reduce(hcat, sol.u)) .- data |
| 46 | +end |
| 47 | +
|
| 48 | +p_init = zeros(4) |
| 49 | +
|
| 50 | +nlls_prob = NonlinearLeastSquaresProblem(loss_function, p_init, vec(reduce(hcat, sol.u))) |
| 51 | +``` |
| 52 | + |
| 53 | +Now, we can use any NLLS solver to solve this problem. |
| 54 | + |
| 55 | +```@example parameterized_ode |
| 56 | +res = solve(nlls_prob, LevenbergMarquardt(); maxiters = 1000, show_trace = Val(true), |
| 57 | + trace_level = TraceAll()) |
| 58 | +``` |
| 59 | + |
| 60 | +We can also use Trust Region methods. |
| 61 | + |
| 62 | +```@example parameterized_ode |
| 63 | +res = solve(nlls_prob, TrustRegion(); maxiters = 1000, show_trace = Val(true), |
| 64 | + trace_level = TraceAll()) |
| 65 | +``` |
| 66 | + |
| 67 | +Let's plot the solution. |
| 68 | + |
| 69 | +```@example parameterized_ode |
| 70 | +prob2 = remake(prob; tspan = (0.0, 10.0)) |
| 71 | +sol_fit = solve(prob2, Tsit5(); p = res.u) |
| 72 | +sol_true = solve(prob2, Tsit5(); p = p) |
| 73 | +plot(sol_true; linewidth = 3) |
| 74 | +plot!(sol_fit; linewidth = 3, linestyle = :dash) |
| 75 | +savefig("LV_ode_fit.png") |
| 76 | +``` |
0 commit comments