# Using DifferentialEquations.jl to solve routine optimal control problems

**URL:** <https://discourse.julialang.org/t/using-differentialequations-jl-to-solve-routine-optimal-control-problems/34096>\
**Category:** Finance and Economics\
**Tags:** diffeq\
**Created:** [February 2, 2020, 10:43pm UTC](https://discourse.julialang.org/t/using-differentialequations-jl-to-solve-routine-optimal-control-problems/34096 "2020-02-02T22:43:38Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [February 2, 2020, 10:43pm UTC](https://discourse.julialang.org/t/using-differentialequations-jl-to-solve-routine-optimal-control-problems/34096/1 "2020-02-02T22:43:38Z")

</div>

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.

```nohighlight
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`?

```nohighlight
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`

```nohighlight
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

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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 3, 2020, 12:11am UTC](https://discourse.julialang.org/t/using-differentialequations-jl-to-solve-routine-optimal-control-problems/34096/2 "2020-02-03T00:11:42Z")

</div>

> [@Albert\_Zevelev](#):
>
> Can I transmit derivatives `du` the same way as variables `u` & parameters `p` ?

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.

> [@Albert\_Zevelev](#):
>
> In the boundary conditions function `bc1!(residual, u, p, t)` I’m unable to transmit `u`

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

> [@Albert\_Zevelev](#):
>
> How do I include non-constant guesses in `BVProblem(DE!, bc1!, u0, tspan, p1)` ?

What do you mean?

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [February 3, 2020, 12:32am UTC](https://discourse.julialang.org/t/using-differentialequations-jl-to-solve-routine-optimal-control-problems/34096/3 "2020-02-03T00:32:33Z")

</div>

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`

```nohighlight
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`

```nohighlight
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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 3, 2020, 12:34am UTC](https://discourse.julialang.org/t/using-differentialequations-jl-to-solve-routine-optimal-control-problems/34096/4 "2020-02-03T00:34:15Z")

</div>

You can do something like:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [February 3, 2020, 12:49am UTC](https://discourse.julialang.org/t/using-differentialequations-jl-to-solve-routine-optimal-control-problems/34096/5 "2020-02-03T00:49:20Z")

</div>

Yep. It works now if I write:

```nohighlight
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.

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [February 3, 2020, 1:00am UTC](https://discourse.julialang.org/t/using-differentialequations-jl-to-solve-routine-optimal-control-problems/34096/6 "2020-02-03T01:00:33Z")

</div>

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:

```nohighlight
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]`.

```nohighlight
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

```

---

<div class="post-metadata">

**Author:** ![Albert\_Zevelev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albert_zevelev/32/11844_2.png) [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)\
**Post date:** [February 3, 2020, 3:40am UTC](https://discourse.julialang.org/t/using-differentialequations-jl-to-solve-routine-optimal-control-problems/34096/7 "2020-02-03T03:40:12Z")

</div>

@ChrisRackauckas regarding my 3rd point, the documentation for [BVProblem](https://docs.juliadiffeq.org/release-4.6/tutorials/bvp_example.html) 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!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 3, 2020, 8:38pm UTC](https://discourse.julialang.org/t/using-differentialequations-jl-to-solve-routine-optimal-control-problems/34096/8 "2020-02-03T20:38:18Z")

</div>

No worries at all.

> [@Albert\_Zevelev](#):
>
> “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.

Needs work: [Improve initial guess interface · Issue #27 · SciML/BoundaryValueDiffEq.jl · GitHub](https://github.com/JuliaDiffEq/BoundaryValueDiffEq.jl/issues/27)

> [@Albert\_Zevelev](#):
>
> It would be convenient to write the boundary conditions in terms of the mathematical names.  
> For example `B[end]` instead of `u[end][2]` .

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

> [@Albert\_Zevelev](#):
>
> But it feels slightly awkward to first “unpack” `u` and `p` before writing the equations, and then to “unpack” `du` after.

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