Using DifferentialEquations.jl to solve routine optimal control problems

Specific Domains Finance and Economics
1

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);
2

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?

3

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.

4

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
5

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.

6

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
7

@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.
I love your package!

8

No worries at all.

Needs work: Improve initial guess interface · Issue #27 · SciML/BoundaryValueDiffEq.jl · GitHub

Yeah, I think that it just needs a good ModelingToolkit DSL over it.

Maybe make a macro? Parameters.jl has some nice stuff for pack/unpack.

1 Like