Sensitivity Analysis of BVProblems

[seishiroono]

I am trying to implement the "shooting method" to solve a simple differential equation (\ddot{x} = -x*(1-x^2) ) as a boundary value problem. I also need to optimize the initial guess when I solve the differential equation.

When I used DifferentialEquations.jl and Optimization.jl for this problem, I encountered the following error. (Please see optimization_diff in my code attached below.)

ERROR: Incompatible problem+solver pairing.
For example, this can occur if an ODE solver is passed with an SDEProblem.
Solvers are only capable of handling specific problem types. Please double
check that the chosen pairing is capable for handling the given problems.

Problem type: ODEProblem
Solver type: Shooting
Problem types compatible with the chosen solver: nothing

On the other hand, if I solved the differential equation without the optimization (the function diff_eq), I did not enter any problem. Also, when I printed the properties of my problem, I found it is Problem type of differential equation: BVProblem.... So, I wonder if the error comes from the incompatibility between Optimization.jl and DifferentialEquations.jl.

Below is my environment.

MacOS (M2Max) Julia: 1.10.0 Zygote v0.6.68 DifferentialEquations v7.12.0 Optimization v3.20.2 OptimizationOptimJL v0.1.14 SciMLSensitivity v7.51.0

using DifferentialEquations, OptimizationOptimJL, Plots

function diff_eq(u0, initial, final, xspan=(0.0, 1.0))
    function f(du, u, p, x)
        du[1] = u[2]
        du[2] = -u[1] * (1 - u[1]^2)
    end

    function bc!(residual, u, p, x)
        residual[1] = u[1][1] - initial
        residual[2] = u[end][1] - final
    end

    bvp = BVProblem(f, bc!, u0, xspan)
    println("Problem type of differential equation: $bvp")
    sol = solve(bvp, Shooting(Vern7()))
    return sol
end


function optimization_diff(initial_guess, initial, final, xspan=(0.0, 1.0))
    x0 = initial_guess
    function objective(v, p)
        u0 = [initial, v[1]]
        sol = diff_eq(u0, initial, final, xspan)
        maximum(sol)
    end
    optprob = OptimizationFunction(objective, Optimization.AutoZygote())
    prob = OptimizationProblem(optprob, x0, [1.0])
    sol = solve(prob, BFGS())
end

diffeq_test = diff_eq([0.0, 0.0], 0.0, 1.0, (0.0, 1.0))
opt_test = optimization_diff([0.0], 0.0, 1.0, (0.0, 1.0))