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

I came across a simple PDAE problem using <a href="https://en.wikipedia.org/wiki/Merto

First of all, I want to thank you <a class="user-mention notranslate" data-hovercard-t

PDAE example,about sciml/methodoflines.jl

Comments (5)

ChrisRackauckas commented on July 22, 2024

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

from methodoflines.jl.

sencilla commented on July 22, 2024

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

from methodoflines.jl.

javier-garcia-tilburg commented on July 22, 2024

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)

from methodoflines.jl.

ChrisRackauckas commented on July 22, 2024

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

from methodoflines.jl.

javier-garcia-tilburg commented on July 22, 2024

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

from methodoflines.jl.

PDAE example about methodoflines.jl HOT 5 OPEN

Comments (5)

Related Issues (20)

Recommend Projects

React

Vue.js

Typescript

TensorFlow

Django

Laravel

D3

Recommend Topics

javascript

web

server

Machine learning

Visualization

Game

Recommend Org

Facebook

Microsoft

Google

Alibaba

D3

Tencent