SciML / BoundaryValueDiffEq.jl

Boundary value problem (BVP) solvers for scientific machine learning (SciML)
Other
41 stars 31 forks source link

Sensitivity Analysis of BVProblems #157

Open seishiroono opened 5 months ago

seishiroono commented 5 months ago

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))
ChrisRackauckas commented 5 months ago

The adjoint of a BVP hasn't been added yet. It's actually written down just hasn't been put to the proper repos. We're working this out over the next few months.

seishiroono commented 5 months ago

Thanks for your reply. I understand the situation. I look forward to having such a wonderful update!

avik-pal commented 3 months ago

I also need to optimize the initial guess when I solve the differential equation

Even when you have adjoints, the gradient wrt the initial guess would be zero by IFT (unless I am getting something wrong)