Dear All,
I am trying to find a way to solve as fast as possible a two point boundary value problem (TPBVP) related to the optimal control of a storage system (it arises from the associated optimality necessary conditions - Pontryagin’s minimum principle).
The differential equation is (state and costate equations - the state is the level of energy in the storage):
function area_ODE!(du, u, p, t)
a, _ = p
@inbounds @fastmath du[1] = min( a[3](t)::Float64 , max( a[4](t) , -1/a[2]::Float64*u[2] - a[5](t)::Float64 ) )
@inbounds @fastmath du[2] = -a[1]::Float64*u[1]
nothing
end
a contains parameters and functions of time.
The residual function is:
function bc2!(residual, u, p, t)
_, b = p
@inbounds @fastmath residual[1] = u[1][1] - b[1]::Float64
@inbounds @fastmath residual[2] = u[end][2] - b[2]::Float64*u[end][1]
end
The first residual is the initial state condition (initial level of energy equal to zero), while the second one is the transversality condition, defined at at the final time.
I define the TPBVP as:
bvp2 = TwoPointBVProblem(area_ODE!, bc2!, u0, tspan, p)
Now, when I solve it with command:
sol = @time solve(bvp2, GeneralMIRK4(), dt = 0.01)
The solution takes very long, but seems correct (residuals very close to zero and control is as expected):
julia> include("solve_problem.jl")
┌ Warning: Backwards compatability support of the new return codes to Symbols will be deprecated with the Julia v1.9 release. Please see https://docs.sciml.ai/SciMLBase/stable/interfaces/Solutions/#retcodes for more information
└ @ SciMLBase ~/.julia/packages/SciMLBase/szsYq/src/retcodes.jl:354
1023.778347 seconds (24.47 G allocations: 1.718 TiB, 7.64% gc time, 0.42% compilation time: <1% of which was recompilation)
If instead I use Shooting:
sol = @time solve(bvp2, Shooting(AutoTsit5(Rosenbrock23())), dt = 0.000001, abstol = 1e-12, reltol = 1e-12)
a solution is found very fast:
julia> include("solve_problem.jl")
┌ Warning: Backwards compatability support of the new return codes to Symbols will be deprecated with the Julia v1.9 release. Please see https://docs.sciml.ai/SciMLBase/stable/interfaces/Solutions/#retcodes for more information
└ @ SciMLBase ~/.julia/packages/SciMLBase/szsYq/src/retcodes.jl:354
2.176208 seconds (26.81 M allocations: 2.303 GiB, 7.71% gc time)
but the final (second) residual is very far from zero.
My question is, is there a way to get Shooting return fast a correct solution? Or to speed up solve(bvp2, GeneralMIRK4(), dt = 0.01).
I have seen here https://docs.sciml.ai/DiffEqDocs/dev/solvers/bvp_solve/#Recommended-Methods that MultipleShooting(5, FBDF()) could be better than single shooting, but when I try:
sol = @time solve(bvp2, MultipleShooting(5, FBDF()), dt = 0.000001, abstol = 1e-12, reltol = 1e-12)
I get:
julia> include("solve_problem.jl")
ERROR: LoadError: UndefVarError: `MultipleShooting` not defined
Finally, my ODE seems to be a DynamicalODE. I am wandering if it is necessary or not to write it as explained in https://docs.sciml.ai/DiffEqDocs/dev/solvers/dynamical_solve/#Dynamical,-Hamiltonian,-and-2nd-Order-ODE-Solvers.
Any help would be appreciated, sorry for the long post.
Thank you!
Best
Francesco