Interpolating solution to constrained BVP

Specific Domains Numerics
, ,
1

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!

2

What’s different between bvp1 and bvp2 is that bvp1 is built with BVPFunction(f!, bc!; f_prototype = zeros(2)), while bvp2 uses f!/bc! directly with no f_prototype. That’s the only structural difference that matters here; the optimize = IpoptOptimizer() solve algorithm itself isn’t the trigger.

Per SciMLBase’s docstring, f_prototype exists for optimal‑control BVPs problems where the unknown vector u is split into (a) true ODE/DAE state variables governed by f, and (b) extra “control” variables that Ipopt is free to choose, not governed by any dynamics. f_prototype’s length tells the solver how many of the M total unknowns fall into bucket (a).

In BoundaryValueDiffEqMIRK/src/interpolation.jl, sum_stages! checks:

has_control = !isnothing(cache.prob.f.f_prototype)
length_z = has_control ? length(cache.prob.f.f_prototype) : length(z)
...
if has_control
    inc = τ / dt .* (id.u[i + 1] .- id.u[i])
    copyto!(z, (length_z + 1):M, inc, (length_z + 1):M)   # interpolation.jl:116
end

You passed f_prototype = zeros(2), whose length (2) equals the full state length M (2) — i.e., zero actual control variables. So (length_z+1):M is 3:2, an empty range. Base.copyto! tolerates that fine, but the specialized copyto!(dest, rdest, src, rsrc) method in LinearAlgebra (the one actually dispatched to here) calls minimum(rdest) unconditionally — and minimum of an empty UnitRange throws ArgumentError: range must be non-empty.

Note this only fires on interpolation after the solve. The solve itself, the retcode, and the saved mesh values are all computed on a different code path (__mirk_loss!/collocation) that never touches this branch, which is why sol1 “looks fine” until you call sol1(0.5).

You don’t really need f_prototype because your problem has no separate control variables, only two coupled ODE states with box constraints. Just drop it:

bvp1 = BVProblem(BVPFunction(f!, bc!), u0, tspan; lb = [-10.0, -10.0], ub = [10.0, 10.0])

lb/ub are independent BVProblem fields (bounds on solution variables) and don’t require f_prototype — that should let interpolation fall through the normal (non-has_control) branch, same as bvp2.

3

Thanks for explaining the issue!

I had tried removing f_prototype = zeros(2) but then the solver does not work. I get this error:

The function `copy` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  copy(::Core.PhiNode)
   @ Base expr.jl:42
  copy(::LinearAlgebra.Transpose{Bool, BitMatrix})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.12.6+0.x64.linux.gnu/share/julia/stdlib/v1.12/LinearAlgebra/src/bitarray.jl:240
  copy(::Base.JuliaSyntax.SyntaxData)
   @ Base /cache/build/builder-amdci5-4/julialang/julia-release-1-dot-12/base/JuliaSyntax/src/syntax_tree.jl:257
  ...


Stacktrace:
  [1] _broadcast_getindex_evalf
    @ ./broadcast.jl:699 [inlined]
  [2] _broadcast_getindex
    @ ./broadcast.jl:672 [inlined]
  [3] _getindex
    @ ./broadcast.jl:696 [inlined]
  [4] _broadcast_getindex
    @ ./broadcast.jl:671 [inlined]
  [5] _getindex
    @ ./broadcast.jl:620 [inlined]
  [6] getindex
    @ ./broadcast.jl:616 [inlined]
  [7] copy
    @ ./broadcast.jl:933 [inlined]
  [8] materialize
    @ ./broadcast.jl:894 [inlined]
  [9] __init(prob::BVProblem{…}, alg::MIRK4{…}; dt::Float64, abstol::Float64, adaptive::Bool, controller::BoundaryValueDiffEqCore.DefectControl{…}, nlsolve_kwargs::@NamedTuple{…}, optimize_kwargs::@NamedTuple{…}, verbose::DiffEqBase.DEVerbosity{…}, kwargs::@Kwargs{})
    @ BoundaryValueDiffEqMIRK ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/mirk.jl:129
 [10] __init
    @ ~/.julia/packages/BoundaryValueDiffEqMIRK/z0Efo/src/mirk.jl:44 [inlined]
 [11] __solve(::BVProblem{…}, ::MIRK4{…}; kwargs::@Kwargs{…})
    @ BoundaryValueDiffEqCore ~/.julia/packages/BoundaryValueDiffEqCore/CwTMw/src/BoundaryValueDiffEqCore.jl:47
 [12] solve_call(_prob::BVProblem{…}, args::MIRK4{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/F2QKS/src/solve.jl:172
 [13] solve_up(prob::BVProblem{…}, sensealg::Nothing, u0::Vector{…}, p::SciMLBase.NullParameters, args::MIRK4{…}; originator::SciMLBase.ChainRulesOriginator, kwargs::@Kwargs{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/F2QKS/src/solve.jl:646
 [14] solve_up
    @ ~/.julia/packages/DiffEqBase/F2QKS/src/solve.jl:619 [inlined]
 [15] #solve#31
    @ ~/.julia/packages/DiffEqBase/F2QKS/src/solve.jl:603 [inlined]
 [16] top-level scope
    @ In[20]:13
 [17] eval(m::Module, e::Any)
    @ Core ./boot.jl:489
Some type information was truncated. Use `show(err)` to see complete types.

I think when the residual is being initialized, f_prototype is expected when there is any constraints including box bounds : BoundaryValueDiffEqMIRK/z0Efo/src/mirk.jl:129 has

vcat(
    [__alloc(bcresid_prototype)], __alloc.(
        copy.(
            [
                f_prototype
                    for _ in 1:Nig
            ]
        )
    )
)
4

Yeah it’s a bug.