Two point boundary value problem solution not obeying boundary condition

General Usage
, ,
1

Hello,

I am trying to replicate the solution of an optimal control problem (Redirecting).

The canonical system is given by

\begin{aligned} & \dot{x}(t)=x(t)\left(1-\frac{x(t)}{a}\right)-\frac{1}{(c+\tau)} \frac{x(t)^2}{1+x(t)^2}-u(t) x(t) \\ & \dot{\lambda}(t)=\lambda(t)\left(r-1+\frac{2 x(t)}{a}+\frac{2 x(t)}{(c+\tau)\left(1+x(t)^2\right)^2}\right)-u(t)(p-\lambda(t)) \end{aligned}

where x is the state variable, \lambda the costate, and u the control.

The maximizer is given by u=\frac{x}{2}(p-\lambda).

The boundary conditions are the initial condition x(0) =x_0 and the transversality condition \lim_{t\rightarrow \infty} e^{-rt}\lambda(t)=0.

This system has a hyperbolic saddle point \hat y and I am trying to solve for the stable path from the initial state to the steady state.

Let

\dot y(t) = \begin{pmatrix} \dot x(t)\\ \dot \lambda(t) \end{pmatrix} = g(y(t))

and denote by J(\hat y) the Jacobian matrix at the the steady state.

Let E^s(\hat y) be the eigenspace associated with the negative eigenvalues, which is tangent to the stable manifold at \hat y.

Then, for T large enough, the transversality condition can be approximated by F^\prime(y(T)-\hat y) =0, where F^\prime is a matrix spanning the orthogonal complement to the stable eigenspace.

My try of implementing this so far:

using DifferentialEquations
using Optim
using ForwardDiff
using LinearAlgebra
using Plots

r = 0.1
τ = 0.25
c = 0.375
a = 7.5
pr = 0.1

par = [r, τ, c, a, pr]

T = 500

G(X) = X * (1 - X / a)
P(X) = (1 / (c + τ)) * (X^2 / (1 + X^2))

Getting the Jacobian matrix and F

function f(u)

    fish = u[1] / 2 * (pr - u[2])
    du1 = G(u[1]) - P(u[1]) - u[1] * fish
    du2 = u[2] * (r - 1 + 2 * u[1] / a + 2 * u[1] / ((c + τ) * (1 + u[1]^2)^2)) - fish * (pr - u[2])
    return [du1, du2]
end
Df(u) = ForwardDiff.jacobian(f, u)

function get_steady_state(par)
    r, τ, c, a, pr = par
    fish(x, barλ) = x * (pr - barλ) / 2
    eqx(x, λ) = G(x) - P(x) - x * fish(x, λ)
    eqλ(x, λ) = λ * (r - 1 + 2 * x / a + 2 * x / ((c + τ) * (1 + x^2)^2)) - fish(x, λ) * (pr - λ)
    eqs(x, λ) = [eqx(x, λ), eqλ(x, λ)]
    res = Optim.optimize(x -> sum(abs.(eqs(x[1], x[2]))), [3, 0.5], Optim.Options(g_tol=1e-15))
    x = res.minimizer[1]
    λ = res.minimizer[2]
    return x, λ
end

barx, barλ = get_steady_state(par)

J = Df([barx, barλ])
ev = eigvals(J)
stable_evecs = eigvecs(J)[:, ev.<0]
F = nullspace(stable_evecs')

Boundary conditions

x0 = barx + 0.1

function bc1!(residuala, sola, p)
    par, x0, barx, barλ = p
    residuala[1] = sola[1][1] - x0
end

function bc2!(residualb, solb, p)
    par, x0, barx, barλ = p
    residualb[1] = (F'*(solb[end].-[barx, barλ]))[1]
end

Here I also ran into a problem when trying to call sola(0.0)[1] as in the tutorial (Boundary Value Problems · DifferentialEquations.jl) it returns an error

ERROR: MethodError: objects of type Vector{Float64} are not callable
Use square brackets [] for indexing an Array.
Stacktrace:
  [1] bc1!(residuala::SubArray{…}, sola::Vector{…}, p::Tuple{…})

Solving


strt = [barx, barλ]
prob = TwoPointBVProblem(f!, (bc1!, bc2!), strt, (0.0, 500.0), (par, x0, barx, barλ); bcresid_prototype=(zeros(1), zeros(1)))
@time sol = solve(prob, MIRK4(), dt=0.01)

Which yields
fishery_bvp
where for some reason the solution diverges at the end, violating the boundary condition.

Help would be appreciated!

2

I misunderstood the solution object in the boundary conditions in the two point boundary value problem vs the general boundary value problem. sola and solb only contain the value of the solution at the start and end points respectively, so solb[end] was a scalar.

The correct boundary conditions are

function bc1!(residuala, sola, p)
    par, x0, barx, barλ = p
    residuala[1] = sola[1] - x0
end

function bc2!(residualb, solb, p)
    par, x0, barx, barλ = p
    residualb[1] = (F'*(solb.-[barx, barλ]))[1]
end

which yield the expected solution.