Add fundamental HJB example to readme?

[azev77]

In Jan 2018 @MaximilianJHuber posted a question about (one of) the most common PDEs that economists encounter.

Is NeuralNetPDE.jl ready to solve this type of problem?

[ChrisRackauckas]

Is it the example a terminal PDE?

[azev77]

@ChrisRackauckas Here is the simplest example formulated as a boundary value PDE:

Note it's really just one PDE, plugin the second equation c_t into the first. I forgot the parameter T, you can let T=120.

[azev77]

Here is my attempt:

using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux

@parameters t a θ
@variables V(..)
@derivatives Da'~a
@derivatives Dt'~t

# PDE
eq  =
    0.01*V(t,a,θ)  ~
    ( (Da(V(t,a,θ)))^(1.0 - 1.0/2.0) )/(1.0 - 2.0) +
    Da(V(t,a,θ)) * (10.0 + 0.05*a - (Da(V(t,a,θ)))^(-1.0/2.0) ) +
    Dt(V(t,a,θ))
#
# Boundary conditions
bcs = [V(120,a,θ) ~ 0.f0,
       Da(V(120,a,θ)) ~ 0.f0]
# Space and time domains
domains = [t ∈ IntervalDomain(0.0,120.0),
           a ∈ IntervalDomain(0.0,20.0)]
# Discretization
dx = 0.1
# Neural network
dim = 2 # number of dimensions
chain = FastChain(FastDense(dim,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

discretization = PhysicsInformedNN(dx,chain)

pde_system = PDESystem(eq,bcs,domains,[t,a],[V])
prob = discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

res = GalacticOptim.solve(prob, Optim.BFGS(); cb = cb, maxiters=1000) # ERROR

I get

julia> res = GalacticOptim.solve(prob, Optim.BFGS(); cb = cb, maxiters=1000)
ERROR: DomainError with -0.010144570572797796:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

[ChrisRackauckas]

You probably want to do ADAM before BFGS.

[azev77]

Sorry, I don't know how to do that

[azev77]

@ChrisRackauckas if my pde only has first-order derivatives, V_t, V_a, do I only need one boundary condition? (I think so)

[ChrisRackauckas]

Yes. You can overconstrain it if you want though.

[ChrisRackauckas]

Sorry, I don't know how to do that

Things like: https://diffeqflux.sciml.ai/dev/examples/neural_ode_sciml/ . 300 iterations of ADAM and then BFGS. Etc.

[azev77]

Same error w/ ADAM:

julia> res = GalacticOptim.solve(prob, ADAM(); cb = cb, maxiters=300)
ERROR: DomainError with -0.05701087375174788:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

I think the issue is that the pde has terms like

[matthieugomez]

Maybe you can try doing max(eps(), va)^(-1/gamma) to avoid this?

On Fri, Nov 6, 2020 at 2:28 PM azev77 notifications@github.com wrote:

Same error w/ ADAM:

julia> res = GalacticOptim.solve(prob, ADAM(); cb = cb, maxiters=300) ERROR: DomainError with -0.05701087375174788: Exponentiation yielding a complex result requires a complex argument. Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

I think the issue is that the pde has terms like [image: image] https://user-images.githubusercontent.com/7883904/98420404-1ca7cc80-2055-11eb-92df-aa7fb8308dcf.png

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/SciML/NeuralPDE.jl/issues/165#issuecomment-723328027, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABPPPXL5NVSSRQXLYT6DONDSORZ7LANCNFSM4S3WUKKQ .

[azev77]

@matthieugomez thanks! Now it's working but very slow & inaccurate. @ChrisRackauckas I'm able to use ADAM, but not sure how to combine ADAM w/ BFGS yet.

Issue here: https://discourse.julialang.org/t/solve-a-routine-economics-pde-from-an-hjb-w-modelingtoolkit/49718

[azev77]

Btw, if I do: max(0.0, va)^(-1/gamma): it looks good for the first few iter, then I get NaNs (dividing by 0) max(eps(), va)^(-1/gamma): it's slow & doesn't converge

1. I know from theory that V_a >0. Is there a way in modeling toolkit to include this information? (as a constraint maybe?)
2. These types of HJB PDEs are usually tackled w/ viscosity solutions. @ChrisRackauckas are there viscosity solvers in the package?
3. I've encountered these types of issues (dividing by 0) using Mathematica & Matlab. Oren Levintal advises:

Useful tip 1: Since consumption and capital should always be positive, it is better to change these variables into logs. This change of variables avoids division by zero or taking the root of a negative number, which may cause the algorithm to fail. In general, it is recommended to use an appropriate change of variables that rules out undesired arithmetic operations.

[ChrisRackauckas]

I know from theory that V_a >0. Is there a way in modeling toolkit to include this information? (as a constraint maybe?)

It's not MTK, it's the neural network. You want to constrain your neural network to only give positive values.

chain = FastChain(FastDense(dim,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1),(x,p)->x^2)

is one way to do it. Then V(dims...) is a positive function since the outputs are always the square of whatever comes out of a neural network.

[ChrisRackauckas]

We could automate log transforms too. We don't do that right now though.