Using DifferentialEquations.jl to solve routine optimal control problems

I’m using DifferentialEquations.jl to solve a standard consumption-savings problem.
I wanna see if there is a way to solve these issues before I submit feature-requests.

``````using DifferentialEquations, BoundaryValueDiffEq, Plots
function DE!(du,u,p,t)
c, B   = u
r,ρ,γ,y,B0,BT = p
du[1] = dc = ((r-ρ)/(γ))*(c)
du[2] = dB = r*B + y - c
end
function bc1!(residual, u, p, t)
r,ρ,γ,y,B0,BT = p
residual[1] = u[1][2]   - B0
residual[2] = u[end][2] - BT
end
tspan = (0.0,60.0);
u0=[10, 10];  #initial guess. NOT boundary @ t=0. constant.
function Par(;r=0.03, ρ=0.02, γ=1.0, y=10.0, B0=5.0, BT=0)
(r,ρ,γ,y,B0,BT)
end
p1 =Par(r=0.01);
bvp1 = BVProblem(DE!, bc1!, u0, tspan, p1);
@time sol1 = solve(bvp1, GeneralMIRK4(), dt=0.5);
plot(sol1, lab=["consumption" "wealth"])
``````

The code above works fine.

1 Can I transmit derivatives `du` the same way as variables `u` & parameters `p`?

``````function DE!(du,u,p,t)
dc, dB = du    #does not work
c, B   = u
r,ρ,γ,y,B0,BT = p
dc = ((r-ρ)/(γ))*(c)
dB = r*B + y - c
end
``````

2 In the boundary conditions function `bc1!(residual, u, p, t)` I’m unable to transmit `u`

``````function bc1!(residual, u, p, t)
c, B   = u #does not work!
r,ρ,γ,y,B0,BT = p
residual[1] = B[1]       - B0  # Initial Condition
residual[2] = B[end]     - BT  # Terminal Condition
end
``````

3 How do I include non-constant guesses in `BVProblem(DE!, bc1!, u0, tspan, p1)`?
Currently

``````u0=[10, 10]; #Constant guess for c[t] and B[t]
bvp1 = BVProblem(DE!, bc1!, u0, tspan, p1);
``````

What does this mean? If they are explicit values, then they are explicitly defined so you have to compute them for us anyways. If they are implicitly defined, then you have a DAE and not an ordinary BVP.

That’s not `u`, that’s `sol`. `sol[1]` would be `u0`, and `sol[end]` would be the final `u`.

What do you mean?

Sorry I wasn’t clear.
1 in the original program I store the equations in `DE!(du,u,p,t)`.
Then I “unpack” the variables in `u` and parameters in `p`

``````function DE!(du,u,p,t)
c, B          = u
r,ρ,γ,y,B0,BT = p
du[1] = dc = ((r-ρ)/(γ))*(c)
du[2] = dB = r*B + y - c
end
``````

I would like to similarly unpack the derivatives in `du`

``````function DE!(du,u,p,t)
dc, dB        = du
c, B          = u
r,ρ,γ,y,B0,BT = p
dc = ((r-ρ)/(γ))*(c)
dB = r*B + y - c
end
``````

But this doesn’t work.

You can do something like:

``````function DE!(du,u,p,t)
c, B          = u
r,ρ,γ,y,B0,BT = p
dc = ((r-ρ)/(γ))*(c)
dB = r*B + y - c
du .= (dc, dB)
end
``````

Yep. It works now if I write:

``````function DE!(du,u,p,t)
c, B =u
r,ρ,γ,y,B0,BT = p
dc = ((r-ρ)/(γ))*(c)
dB = r*B + y - c
du .= (dc, dB)
end
``````

But it feels slightly awkward to first “unpack” `u` and `p` before writing the equations, and then to “unpack” `du` after.

I have two unknown variables in this system:
consumption: `c = u[1]`
wealth: `B=u[2]`
I have two boundary conditions:
initial wealth: `B[t=0]=B0`
final wealth: `B[t=T]=BT`
My time-span is (0, 40.0) so I want: `B[1]=B0` and `B[end]=BT` (assuming the first time-index is 0, etc).
Currently, I write:

``````function bc1!(residual, u, p, t)
r,ρ,γ,y,B0,BT = p
residual[1] = u[1][2]   - B0
residual[2] = u[end][2] - BT
end
``````

Here it is `u[time-period][variable index number]` eg `u[end][2]`.
It would be convenient to write the boundary conditions in terms of the mathematical names.
For example `B[end]` instead of `u[end][2]`.

``````function bc1!(residual, u, p, t)
c, B   = u
r,ρ,γ,y,B0,BT = p
residual[1] = B[1]       - B0  # Initial Condition
residual[2] = B[end]     - BT  # Terminal Condition
end
``````

@ChrisRackauckas regarding my 3rd point, the documentation for BVProblem says
“The third argument of `BVProblem` is the initial guess of the solution, which is constant in this example”
I was curious about how to supply non-constant guesses.

Either way, these are mainly cosmetic issues.