|
422 | 422 | psol = solve(pprob, PyomoCollocation("ipopt", LagrangeLegendre(4))) |
423 | 423 | @test psol.sol.u[end] ≈ [π, 0, 0, 0] |
424 | 424 | end |
| 425 | + |
| 426 | +@testset "Parameter estimation - JuMP" begin |
| 427 | + @parameters α = 1.5 β = 1.0 [tunable=false] γ = 3.0 δ = 1.0 |
| 428 | + @variables x(t) y(t) |
| 429 | + |
| 430 | + eqs = [D(x) ~ α * x - β * x * y, |
| 431 | + D(y) ~ -γ * y + δ * x * y] |
| 432 | + |
| 433 | + @mtkcompile sys0 = System(eqs, t) |
| 434 | + tspan = (0.0, 1.0) |
| 435 | + u0map = [x => 4.0, y => 2.0] |
| 436 | + parammap = [α => 1.8, β => 1.0, γ => 6.5, δ => 1.0] |
| 437 | + |
| 438 | + oprob = ODEProblem(sys0, [u0map; parammap], tspan) |
| 439 | + osol = solve(oprob, Tsit5()) |
| 440 | + ts = range(tspan..., length=51) |
| 441 | + data = osol(ts, idxs=x).u |
| 442 | + |
| 443 | + costs = [EvalAt(t)(x)-data[i] for (i, t) in enumerate(ts)] |
| 444 | + consolidate(u, sub) = sum(abs2.(u)) |
| 445 | + |
| 446 | + @mtkcompile sys = System(eqs, t; costs, consolidate) |
| 447 | + |
| 448 | + sys′ = subset_tunables(sys, [γ, α]) |
| 449 | + jprob = JuMPDynamicOptProblem(sys′, u0map, tspan; dt=1/50, tune_parameters=true) |
| 450 | + jsol = solve(jprob, JuMPCollocation(Ipopt.Optimizer, constructTsitouras5())) |
| 451 | + |
| 452 | + @test jsol.sol.ps[γ] ≈ 6.5 rtol=1e-4 |
| 453 | + @test jsol.sol.ps[α] ≈ 1.8 rtol=1e-4 |
| 454 | +end |
0 commit comments