Hi all!
When trying to find interpolated values of a solution to a constrained BVP using SciMLs collocation method, I get an error. Here is a MWE:
using BoundaryValueDiffEq, OptimizationIpopt
function f!(du, u, p, t)
du[1] = u[2]
du[2] = u[1]
end
function bc!(res, u, p, t)
res[1] = u(0.0)[1] - 1
res[2] = u(1.0)[1]
end
tspan = (0.0, 1.0)
u0 = [0.0, 0.0]
bvp1 = BVProblem(BVPFunction(f!, bc!;f_prototype = zeros(2)),u0,tspan; lb = [-10.0,-10.0],ub = [10.0,10.0])
sol1 = solve(bvp1, MIRK4(; optimize = IpoptOptimizer()); dt = 0.01)
bvp2 = BVProblem(f!, bc!, u0, tspan)
sol2 = solve(bvp2, MIRK4(), dt = 0.01)
println(sol2(0.5))
println(sol1(0.5))
This is the output:
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit https://github.com/coin-or/Ipopt
******************************************************************************
This is Ipopt version 3.14.19, running with linear solver MUMPS 5.9.0.
Number of nonzeros in equality constraint Jacobian...: 802
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 0
Total number of variables............................: 202
variables with only lower bounds: 0
variables with lower and upper bounds: 202
variables with only upper bounds: 0
Total number of equality constraints.................: 202
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 0.0000000e+00 1.00e+00 0.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 0.0000000e+00 3.24e-16 0.00e+00 -1.0 1.31e+00 - 8.83e-01 1.00e+00h 1
Number of Iterations....: 1
(scaled) (unscaled)
Objective...............: 0.0000000000000000e+00 0.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 3.2439329000766293e-16 3.2439329000766293e-16
Variable bound violation: 0.0000000000000000e+00 0.0000000000000000e+00
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 3.2439329000766293e-16 3.2439329000766293e-16
Number of objective function evaluations = 2
Number of objective gradient evaluations = 2
Number of equality constraint evaluations = 2
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 2
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 1
Total seconds in IPOPT = 9.040
EXIT: Optimal Solution Found.
[0.4434094419864666, -0.9595173756819017]
ArgumentError: range must be non-empty
Stacktrace:
[1] minimum
@ ./range.jl:863 [inlined]
[2] copyto!(dest::Vector{Float64}, rdest::UnitRange{Int64}, src::Vector{Float64}, rsrc::UnitRange{Int64})
@ LinearAlgebra ~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/blas.jl:2241
[3] sum_stages!(z::Vector{…}, id::BoundaryValueDiffEqMIRK.MIRKInterpolation{…}, cache::BoundaryValueDiffEqMIRK.MIRKCache{…}, w::Vector{…}, i::Int64, τ::Float64, ::Type{…})
@ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/interpolation.jl:116
[4] interpolant!
@ ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/interpolation.jl:79 [inlined]
[5] interpolation
@ ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/interpolation.jl:68 [inlined]
[6] (::BoundaryValueDiffEqMIRK.MIRKInterpolation{…})(tvals::Float64, idxs::Nothing, deriv::Type, p::SciMLBase.NullParameters, continuity::Symbol)
@ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/interpolation.jl:13
[7] AbstractODESolution
@ ~/.julia/packages/SciMLBase/dpRhw/src/solutions/ode_solutions.jl:241 [inlined]
[8] #_#625
@ ~/.julia/packages/SciMLBase/dpRhw/src/solutions/ode_solutions.jl:225 [inlined]
[9] AbstractODESolution
@ ~/.julia/packages/SciMLBase/dpRhw/src/solutions/ode_solutions.jl:218 [inlined]
[10] (::ODESolution{…})(t::Float64)
@ SciMLBase ~/.julia/packages/SciMLBase/dpRhw/src/solutions/ode_solutions.jl:218
[11] top-level scope
@ In[2]:17
[12] eval(m::Module, e::Any)
@ Core ./boot.jl:489
Some type information was truncated. Use `show(err)` to see complete types.
The constrained solution has a success retcode and I can access the saved values which match those of the unconstrained version but I can’t interpolate the solution to other points. Is it a bug on my end? If not, is it possible to construct an interpolation from the saved points?
Thanks!