How to override defaults with u0 in solve

Specific Domains Modelling & Simulations
1

I have an ODEProblem with a system with default values, and with a starting state prob.u0. I want to be able to solve this problem from an arbitrary starting state u0 to the next timestep using solve(prob, QNDF(); save_start = true, saveat=Ts). The problem is that solve does not use prob.u0, but it rather uses defaults(prob.f.sys) as initial values.

julia> defaults(prob.f.sys)
Dict{Any, Any} with 82 entries:
  (pos(t))[3, 5]     => 24.0023
  (vel(t))[2, 10]    => 0.0
  (pos(t))[2, 3]     => 0.0
  (set_values(t))[1] => 0.0
  (tether_vel(t))[2] => 0.0
  (vel(t))[2, 6]     => 0.0
  (pos(t))[3, 12]    => 53.5791
  ⋮                  => ⋮

julia> println(prob.u0)
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

julia> @time sol = solve(prob, solver; u0=prob.u0, save_start = true, saveat=Ts)
  0.003862 seconds (310 allocations: 127.625 KiB)
retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
 0.0
 0.05
u: 2-element Vector{Vector{Float64}}:
 [12.49482743653575, 0.3989368894537224, 24.002340185174127, 12.49482743653575, -0.3989368894537224, 24.002340185174127, 2.1788935686914535, 0.0, 24.90486745229364, 0.17453292519943295  …  0.0, 0.0, 0.0, 0.0, 57.60663544529114, 57.60663544529114, 50.0, 0.0, 0.0, 0.0]
 [12.08517813866516, 0.423952531012562, 25.969115391426246, 12.085178138665146, -0.4239525310125608, 25.96911539142627, 2.3196742729500457, -1.440317109756495e-17, 25.018474359821806, 0.18357897074375468  …  3.246285727696867, 2.680128140006882, -9.338003043785488e-14, 2.3886546391917896, 59.53250254353384, 59.53250254353384, 51.85721596131688, 67.78422267562938, 67.78422267562944, 66.79836234747967]

sol.u[1] result should be filled with ones, but it is filled with the default values.
I tried creating a dict with u0 values, but the solver still does not start at u0 .= 1.0.

julia> default
Dict{SymbolicUtils.BasicSymbolic{Real}, Float64} with 52 entries:
  (pos(t))[3, 5]        => 1.0
  (vel(t))[2, 10]       => 1.0
  (tether_vel(t))[2]    => 1.0
  (vel(t))[2, 6]        => 1.0
  (pos(t))[3, 12]       => 1.0
  (pos(t))[2, 9]        => 1.0
  (vel(t))[1, 4]        => 1.0
  (pos(t))[3, 10]       => 1.0
  ⋮                     => ⋮

julia> @time sol = solve(prob, solver; u0=default, save_start = true, saveat=Ts)
  0.004908 seconds (6.26 k allocations: 1.112 MiB)
retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
 0.0
 0.05
u: 2-element Vector{Vector{Float64}}:
 [12.49482743653575, 0.3989368894537224, 24.002340185174127, 12.49482743653575, -0.3989368894537224, 24.002340185174127, 2.1788935686914535, 0.0, 24.90486745229364, 0.17453292519943295  …  0.0, 0.0, 0.0, 0.0, 57.60663544529114, 57.60663544529114, 50.0, 0.0, 0.0, 0.0]
 [12.08517813866516, 0.423952531012562, 25.969115391426246, 12.085178138665146, -0.4239525310125608, 25.96911539142627, 2.3196742729500457, -1.440317109756495e-17, 25.018474359821806, 0.18357897074375468  …  3.246285727696867, 2.680128140006882, -9.338003043785488e-14, 2.3886546391917896, 59.53250254353384, 59.53250254353384, 51.85721596131688, 67.78422267562938, 67.78422267562944, 66.79836234747967]

What can I do to override the defaults and use prob.u0?

2

The solution is using the integrator interface: Integrator Interface · DifferentialEquations.jl

But still, I don’t understand why changing prob.u0 with solve does not work.

3

There is an initialization phase in the DAE solve that uses prob.u0 as the guesses. If you don’t change the initialization problem (i.e. you bypass the SII and directly change the pieces of the generated code) then you’re not changing the initial values, just the guesses.

1 Like