How to solve a limit boundary value ODE in SciML (possibly using shooting methods)

In many fields it is common to encounter systems of differential equations with boundary conditions that occur at limits (rather than at finite points).

For example, the Neoclassical growth model boils down to a system of two boundary value ODEs:
\dot{C}(t)= \frac{ 1 }{ \gamma } \left( \alpha z (K(t))^{\alpha -1} - \delta - \rho \right) \times C(t)
\dot{K}(t) =z (K(t))^{\alpha} - \delta K(t) - C(t)
K(0)=K_{0} given (Initial Condition)
\lim_{T\to \infty} e^{-\rho T} \lambda(T) K(T) = e^{-\rho T} C(T)^{-\gamma} K(T) \to 0 (Terminal Condition)

This system has two variables, one state variable K(t), one control variable C(t).
The initial condition depends only on the state variable.
The terminal condition depends on both variables.
This is standard in optimal control problems.

AFAIK currently SciML does not allow boundary conditions at limits.
This post will attempt to “hack” a solution.

First, let’s look at the closed form solution (for a subset of parameters where we know C0):

using DifferentialEquations, Plots;

function ngm!(dX, X, θ, t)
    C, K          = X  # variables: Consumption, Capital
    ρ, γ, z, α, δ = θ  # parameters
    KN = abs(K)
    dC = ((α*z*(KN^(α-1.0)) -δ -ρ)/γ)*C # Euler Eqn
    dK = z*KN^(α) - δ*KN - C             # Law of Motion
    dX .= (dC, dK)    # Time derivatives

α=0.20; γ=6.0; δ=0.25; z=1.0; ρ = δ*(α*γ-1.0)
K0 = 1.0; C0 = (1-(1/γ))*z*(K0)^α;
KSS= ((α)/(ρ+δ))^(1/(1-α)); CSS = (KSS)^α - δ*KSS;

X0 = [C0, K0]                         # initial conditions
T = 100.0                              # time horizon
tspan = (0.0, T)                      # time span
θ = [ρ, γ, z, α, δ]                   # parameters
TC(t,c,k;γ=γ,ρ=ρ) = exp(-ρ*t) * (c^(-γ)) * k  # Terminal condition
prob = ODEProblem(ngm!, [C0, K0], tspan, θ); # problem
sol = solve(prob,abstol=1e-10,reltol=1e-10)
sol_C(t) = sol(t, idxs=1)
sol_K(t) = sol(t, idxs=2)
plot!(tt, sol_C, label="C(t)", color=1)
plot!(tt, sol_K, label="K(t)", color=2)
plot!(tt, t -> TC(t, sol_C(t), sol_K(t)), label="TC(t)", color=3)
plot!([CSS KSS 0.0],  seriestype = :hline, lab=["CSS" "KSS" ""], color=[1 2 3], l=:dash)

1: the control variable C(t) goes to its steady state $C_{SS}
2: the state variable K(t) goes to its steady state K_{SS}
3: the terminal condition holds TC(t)\to 0

In most cases we don’t have a closed form solution, so we don’t know C0.
Consider a shooting method, to find the appropriate C0 such that the terminal condition holds TC(t)\to 0.
First consider a grid of C0

grid_c0 = .75*K0 :0.01 : 1.25*K0;
grid_TC =[]
for c0 in grid_c0
    prob = ODEProblem(ngm!, [c0, K0], tspan, θ); # problem
    #sol = solve(prob, Tsit5(), abstol=1e-10,reltol=1e-10, maxiters = 10_000_000)
    sol = solve(prob, Tsit5())
    if sol.retcode == :Success
        _tc = TC(T, sol(T, idxs=1), sol(T, idxs=2))
        _tc = -100    
    push!(grid_TC, _tc)
    println(c0, " ", _tc)
cl(x) = clamp(x,-100,100)
plot(grid_c0, cl.(grid_TC) )

It appears the appropriate C0 (which satisfies the TC) is between 0.8 and 0.9.
This is good bc we know the true C0=0.833333

Another approach is to use a nonlinear solver (or optimized) to find the unique value of C0 at which the terminal condition is satisfied.

using Optim;
function f(x; K0=K0,T=100.0,tspan=(0.0, T),θ=θ)
    prob = ODEProblem(ngm!, [x[1], K0], tspan, θ); # problem
    #sol = solve(prob,Tsit5(),reltol=1e-10,maxiters=100_000,abstol=1e-8)
    sol = solve(prob,Tsit5())
    if sol.retcode == :Success #C, K
        cc(x) = sol(x, idxs=1)
        KK(x) = sol(x, idxs=2)
        _tc1 = TC(T, abs(cc(T)), abs(KK(T)))
        _tc1 = -10_000.0    
    return abs(_tc1) 
f([0.75]) #
result = optimize(f, [0.75])

Optim.jl returns C0=0.8333338128872135 which is close to the true C0=.83333333333.

  1. Does anyone know of a good way to do this in the Julia ecosystem or otherwise
  2. I’ve tried in on larger problems (2 state/2 control) variables w/o success
  3. I arbitrarily choose T=100.0, is there a better way for finding the limit of TC(t) maybe when some tolerance is reached?

I think a nicer way to handle this is to use the SteadyStateProblem, since the value at t->infty is only well defined if the solution is a steady state, it’s equivalent. It adds a callback that automatically checks if certain criteria are obtained to then cause an exit at a determined steady state, which is then the value “at infinity” with some tolerance (if the ODE is autonomous). That also has some nice overloads for making the adjoint differentiation faster (since it doesn’t need to backpropogate through the steps). It does need an overload for a shortcutting forward mode, but the fallback is fine enough and would just differentiate the algorithm.

You should be able to increase the number of states with more success since those are paired with an initial condition rather than a transversality condition.

Unless you go to crazier methods, either shooting or some sort of implicit scheme are the only things I am aware of to solve these. In all cases you have very poorly behaved root finding problems to solve. The choice of T from my experience is also tricky and problem dependent because the rate of convergence towards the steady state varies widely.

At the heart of the difficulty here is that you are trying to find a saddlepath by shooting and divergent trajectories blow up if you aren’t exactly on that saddlepath. It is a very difficult problem for nonlinear solvers to handle effectively and requires a lot of solver babysitting.

The issue is less the calculation of the steady-state (which can be done by hand here) and more the nature of these mixed initial value/boundary value problems. To find the u(0), think of the initial conditions pinning down a subset of the u(0) and then longrun boundary conditions ultimately being required to figure out the others.

The nature of these problems is that they are steady state problems. Though steady state shooting to initial conditions is somewhat ill-defined because the gradient is non-unique.

I think the way a control-theory person would put it is that the boundary conditions are always determined by the steady-state problem which provide sufficient boundary conditions to make transition dynamics well posed. Maybe that is what you mean.

Regardless, I think that Albert’s question is a practical one. How can you setup a problem to calculate these given these “exotic” boundary conditions. I tried a few years ago and didn’t think it was clear - nor do I think it should necessarily be a goal given that control theory tend to have its own algorithms.

It is the opposite direction. The way these work is that you fix some of the initial conditions and shoot from the other initial conditions towards the boundary values at infinity. I believe that the transversality conditions pin down gradients in theory at the saddlepath solution… you just need to find it.

I think you can easily prove that the gradients of the terminal condition with respect to the initial condition you shoot from is either degenerate, infinite, or on the saddle-path as you increase T. This is a big part of why deterministic control-theory has its own methods… With saddlepaths, you are trying to shoot towards the exact path that doesn’t diverge - where every other path blows up in some way. So you try to pick a T which isn’t too large but isn’t too small. @Albert_Zevelev in my mind, this is the core of why these problems are so difficult to solve with standard methods, and why it isn’t obvious how to fit them into the standard SciML solvers.

BTW, this is the reason that LQ works so well for control problems. By analyzing the linearized dynamical system’s spectrum you can find the paths that don’t blow up.

Is there a good paper on that? The saddlepath would more likely be computed by bifurcation analysis tools, but how to do that mapping would be interesting.

Oh yeah, I think it is super-interesting and at the core of why control-theory is so tricky. I am with you on using bifurcation analysis to analyze the steady-state and see if it helps with shooting.

I don’t know the control-theory literature well enough to see if people have looked at how to tackle this problem with bifurcations for general nonnlinear models, but I suspect some OR/EE people may have? One issue is that a lot of them work with finite-time horizon - where this doesn’t occur.

But your bifurcation thing is is basically what LQ control does manually - since the linearized system is especially easy to analyze for saddlepaths since the jacobian is just the linear evoluation matrix. Furthermore, if you figure out the non-explosive solution you no longer need to worry about shooting since there is only one linear evolution equation which fulfills the criteria. After you find it, you can solve transition dynamics as a linear initial value problem everywhere in the state-space. Other perturbation methods do the same thing with higher-order approximations but their saddlepath is solved through the linear term.