PDAE example

[javier-garcia-tilburg]

Is your feature request related to a problem? Please describe.

It would be nice to have an example of a partial differential algebraic equation (PDAE). There are examples about how to transform a PDE system into an ODE and there are some references to DAEs, but as a newcomer I still find PDAEs elusive, even though some comments imply it is possible.

Describe the solution you’d like

What I have in mind is some elliptic / parabolic problem. In economics it is common to have parabolic equations for the dynamics of the value function J(x, t) and for the first order condition (FOC), which works as an algebraic constraint.

Describe alternatives you’ve considered

I tried embedding the algebraic conditions into a PDESystem but discretize complains that there are no boundary conditions for m(x, t). If I derive boundary conditions from the FOC and include them for m(x, t_min), it raises an exception while checking islinear.

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@parameters x t
@variables J(..) m(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

f(m) = exp(z * m)          # Simplified/dummy example
fm(m) = z * exp(z * m)

eq = [
    # HJB
    0 ~ Dt(J(x, t)) + Dx(J(x, t)) * x * (r + m(x, t) * (μ- r) + f(m(x, t))) + (1 / 2) * Dxx(J(x, t)) * (x^2) * ((m(x, t))^2) * (σ^2),
    # FOC
    0 ~ Dx(J(x, t)) * (μ- r + fm(m(x, t))) + Dxx(J(x, t)) * x * m(x, t) * (σ^2)
]

domains = [x ∈ Interval(x_min, x_max),
    t ∈ Interval(t_min, t_max)]

JT(x, t) = (x^(1 - γ)) / (1 - γ)

# Boundary conditions
bcs = [
    J(x, t_min) ~ JT(x, t_min),
]

@named pdesys = PDESystem(eq, bcs, domains, [x, t], [J(x, t), m(x, t)])

N = 32
order = 2
discretization = MOLFiniteDifference([x => N], t, approx_order=order)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys, discretization)  # Error

Other approaches along the lines PDESystem -> discretize -> ODEProblem -> DAEProblem were also unsuccessful.

[ChrisRackauckas]

Transferring to MethodOfLines as this is more of a MethodOfLines issue. If it generates the right DAE then MTK should be fine.

[sencilla]

I have the same Problem with another system.

Discretize complains that there are no valid boundary conditions

@chris
Even PDAEs can be very tricky, it would be good if we had a simple example that works

[javier-garcia-tilburg]

I came across a simple PDAE problem using Merton's portfolio problem (without consumption and discounting). Since the boundary can be solved explicitly, discretize can easily chop the PDAE into a DAE.

However, all the solvers that I've tried (Rosenbrock23, QNDF, Rodas4P, Rodas5P) fail with MaxIters during the first time step. Perhaps I am doing something wrong, but I just can't spot the mistake.

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Plots

@parameters x t
@variables J(..) m(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

r::Float64 = 0.03     # bond rate
μ::Float64 = r + 0.07 # expected market return
σ::Float64 = 0.2      # market volatility

γ::Float64 = 2        # risk aversion

x_min::Float64 = 5
x_max::Float64 = 100
t_min::Float64 = 0.0
t_max::Float64 = 40

# Notice that time is reversed
m_analytical(x, t) = (μ - r) / (γ * (σ^2))
J_analytical(x, t) = (exp(t .* (1 .- γ) .* (r + (1 / (2 * γ)) * (((μ - r)^2) / ((σ^2)))))) .* (x .^ (1 .- γ)) ./ (1 .- γ)

JT(x) = (x^(1 - γ)) / (1 - γ)

eq = [
    # HJB
    0 ~ -Dt(J(x, t)) + Dx(J(x, t)) * (x * (r + m(x, t) * (μ - r))) + (1 / 2) * Dxx(J(x, t)) * (x^2) * ((m(x, t))^2) * (σ^2),
    # FOC
    0 ~ Dx(J(x, t)) * x * (μ - r) + Dxx(J(x, t)) * (x^2) * (m(x, t)) * (σ^2),
]

domains = [x ∈ Interval(x_min, x_max),
    t ∈ Interval(t_min, t_max)]

# Boundary conditions
bcs = [
    J(x, t_min) ~ JT(x),
    m(x, t_min) ~ (μ - r) / (γ * (σ^2))
]

@named pdesys = PDESystem(eq, bcs, domains, [x, t], [J(x, t), m(x, t)])

#
# Method of lines discretization
#

Nx::Int64 = 32

order::Int64 = 2 # This may be increased to improve accuracy of some schemes

# Use a Float to directtly specify stepsizes dx.
discretization = MOLFiniteDifference([x => Nx], t, approx_order=order)

# Convert the PDAE problem into an DAE problem
println("Discretization:")
@time prob = discretize(pdesys, discretization)

#
# Solving the problem
#

println("Solve:")
@time sol = solve(prob, Rodas4P(), saveat=0.1)

[ChrisRackauckas]

Is the system structurally simplified? Is there a Jacobian singularity at the start?

[javier-garcia-tilburg]

First of all, I want to thank you @ChrisRackauckas for your time, and warn you that I'm a complete noob when it comes to solving PDEs. So feel free to close the issue if you believe that it becomes sterile.

A singularity seems quite likely indeed, since $m(x,t)$ is constant over the whole domain in that simplified problem. Perhaps a slightly better example would be the PDAE formulation below of Merton's portfolio problem with human capital, in which the optimal $m(x,t)$ varies with $x$ and $t$. However the derivatives $m_x(x,t)$ and $m_t(x,t)$ are still zero at $t=0$, which may still cause trouble. Is it possible to work around it?

PDAE formulation code

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Sundials, ODEInterfaceDiffEq, LSODA, ForwardDiff, Plots

@parameters x t
@variables J(..) m(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

bond_rate::Float64 = 0.03
stock_expected_return::Float64 = bond_rate + 0.07
stock_volatility::Float64 = 0.2

risk_aversion::Float64 = 2
initial_wealth::Float64 = 10
external_portfolio_contribution::Float64 = 8

x_min::Float64 = 5
x_max::Float64 = 100
t_min::Float64 = 0.0
t_max::Float64 = 40

JT(x) = (x^(1 - risk_aversion)) / (1 - risk_aversion)

analytical_A = ((-(1 - risk_aversion) .* (bond_rate + (((stock_expected_return - bond_rate) ./ stock_volatility) .^ 2) ./ (2 .* risk_aversion))) ./ risk_aversion)

analytical_g(t) = exp.(.-analytical_A .* t)

analytical_J(x, t) = ((analytical_g(t)) .^ (risk_aversion)) .* (((x .+ external_portfolio_contribution .* (1 .- exp.(-t .* bond_rate)) ./ (bond_rate)) .^ (1 .- risk_aversion)) ./ (1 .- risk_aversion))

analytical_m(x, t) = ((stock_expected_return .- bond_rate) ./ (stock_volatility .^ 2)) .* ((x .+ external_portfolio_contribution .* (1 .- exp.(.-t .* bond_rate)) / (bond_rate)) ./ (x .* risk_aversion))

eq = [
    # HJB
    0 ~ -Dt(J(x, t)) + Dx(J(x, t)) * (x * (bond_rate + m(x, t) * (stock_expected_return - bond_rate)) + external_portfolio_contribution) + (1 / 2) * Dxx(J(x, t)) * (x^2) * ((m(x, t))^2) * (stock_volatility^2),
    # FOC
    0 ~ Dx(J(x, t)) * x * (stock_expected_return - bond_rate) + Dxx(J(x, t)) * (x^2) * (m(x, t)) * (stock_volatility^2),
]

domains = [x ∈ Interval(x_min, x_max),
    t ∈ Interval(t_min, t_max)]

# Boundary conditions
bcs = [
    J(x, t_min) ~ JT(x),
    m(x, t_min) ~ (stock_expected_return - bond_rate) / (risk_aversion * (stock_volatility^2))
]

@named pdesys = PDESystem(eq, bcs, domains, [x, t], [J(x, t), m(x, t)])

#
# Method of lines discretization
#

Nx::Int64 = 128

order::Int64 = 2 # This may be increased to improve accuracy of some schemes

# Use a Float to directtly specify stepsizes dx.
# discretization = MOLFiniteDifference([x => Nx], t, approx_order=order)
discretization = MOLFiniteDifference([x => Nx], t, approx_order=order, use_ODAE=true)

# Convert the PDE problem into an ODE problem
println("Discretization:")
@time prob = discretize(pdesys, discretization)

#
# Solving the problem
#

println("Solve:")
@time sol = solve(prob, Sundials.IDA(), saveat=0.1)

Another problem with the following PDAE formulation above is that MOLFiniteDifference(..., use_ODAE=true) produces an ODESystem instead of a DAE (I assume it should be a DAE because $m(x,t)$ always appears with order 0). Using Sundials.IDA() throws an error.

The system has an equivalent PDE formulation whenever $m(x,t)$ can be solved explicitly from the FOC. Sundials and ODEInterfaceDiffEq solvers can solve 14 years, while native julia solvers only manage to find the solution for the first 5 years.

The problem with this equivalent PDE formulation is that I'm unable to reliably extract $J_x$ and $J_{xx}$ from the solution, so I can't recover $m(x,t)$.

Equivalent PDE formulation code

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Sundials, ODEInterfaceDiffEq, LSODA, ForwardDiff, FiniteDifferences, Plots

@parameters x t
@variables J(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

bond_rate::Float64 = 0.03
stock_expected_return::Float64 = bond_rate + 0.07
stock_volatility::Float64 = 0.2

risk_aversion::Float64 = 2
initial_wealth::Float64 = 10
external_portfolio_contribution::Float64 = 8

x_min::Float64 = 5
x_max::Float64 = 100
t_min::Float64 = 0.0
t_max::Float64 = 40
# t_max::Float64 = 10

JT(x) = (x^(1 - risk_aversion)) / (1 - risk_aversion)

analytical_A = ((-(1 - risk_aversion) .* (bond_rate + (((stock_expected_return - bond_rate) ./ stock_volatility) .^ 2) ./ (2 .* risk_aversion))) ./ risk_aversion)

analytical_g(t) = exp.(.-analytical_A .* t)

analytical_J(x, t) = ((analytical_g(t)) .^ (risk_aversion)) .* (((x .+ external_portfolio_contribution .* (1 .- exp.(-t .* bond_rate)) ./ (bond_rate)) .^ (1 .- risk_aversion)) ./ (1 .- risk_aversion))

eq = [
    # HJB
    0 ~ -Dt(J(x, t)) + Dx(J(x, t)) * (x * bond_rate + external_portfolio_contribution) - (1 / 2) * ((stock_expected_return - bond_rate)^2) / (stock_volatility^2) * ((Dx(J(x, t))^2) / Dxx(J(x, t))),
]

domains = [x ∈ Interval(x_min, x_max),
    t ∈ Interval(t_min, t_max)]

# Boundary conditions
bcs = [
    J(x, t_min) ~ JT(x)
]

@named pdesys = PDESystem(eq, bcs, domains, [x, t], [J(x, t)])

#
# Method of lines discretization
#

Nx::Int64 = 128

order::Int64 = 2 # This may be increased to improve accuracy of some schemes

# Use a Float to directtly specify stepsizes dx.
# discretization = MOLFiniteDifference([x => Nx], t, approx_order=order)
discretization = MOLFiniteDifference([x => Nx], t, approx_order=order, jac=true) # It seems that jac is not useful at all

# Convert the PDE problem into an ODE problem
println("Discretization:")
# @time prob = discretize(pdesys, discretization)
@time prob = discretize(pdesys, discretization, jac=true) # It seems that jac is not useful at all

#
# Solving the problem
#

println("Solve:")
@time sol = solve(prob, Sundials.CVODE_BDF(), saveat=0.1)
# The solution is actually quite accurate

# It seems that interpolation doesn't include derivatives, since this command fails
sol(20.0, 8.0, deriv=Val{1})
# Trying to obtain the second derivative Jxx using ForwardDiff gives 0, which is wrong
ForwardDiff.derivative(x_val -> ForwardDiff.derivative(x_val -> sol(x_val, 8.0, dv=J(x, t)), x_val), 20.0)
# Using finite differences to approximate the nonconstant component of m(x,t) yields nonsense
forward_fdm(12, 1, max_range=9e-20)(x_val -> sol(x_val, 8.0, dv=J(x, t)), 20.0) / (- 20 * forward_fdm(12, 2, max_range=9e-20)(x_val -> sol(x_val, 8.0, dv=J(x, t)), 20.0))