Merge pull request #717 from SciML/louiville

add the liouville transformation

Author: ChrisRackauckas

docs/src/systems/ODESystem.md

@@ -17,6 +17,7 @@ ODESystem
 
 ```@docs
 ode_order_lowering
+liouville_transform
 ```
 
 ## Applicable Calculation and Generation Functions

src/ModelingToolkit.jl

@@ -176,6 +176,7 @@ include("systems/diffeqs/abstractodesystem.jl")
 include("systems/diffeqs/first_order_transform.jl")
 include("systems/diffeqs/modelingtoolkitize.jl")
 include("systems/diffeqs/validation.jl")
+include("systems/diffeqs/basic_transformations.jl")
 
 include("systems/jumps/jumpsystem.jl")
 
@@ -206,7 +207,7 @@ export SteadyStateProblem, SteadyStateProblemExpr
 export JumpProblem, DiscreteProblem
 export NonlinearSystem, OptimizationSystem
 export ControlSystem
-export ode_order_lowering
+export ode_order_lowering, liouville_transform
 export runge_kutta_discretize
 export PDESystem
 export Reaction, ReactionSystem, ismassaction, oderatelaw, jumpratelaw

src/systems/diffeqs/basic_transformations.jl

@@ -0,0 +1,56 @@
+"""
+$(TYPEDSIGNATURES)
+
+Generates the Liouville transformed set of ODEs, which is the original
+ODE system with a new variable `trJ` appended, corresponding to the
+-tr(Jacobian). This variable is used for properties like uncertainty
+propagation from a given initial distribution density.
+
+For example, if ``u'=p*u`` and `p` follows a probability distribution
+``f(p)``, then the probability density of a future value with a given
+choice of ``p`` is computed by setting the inital `trJ = f(p)`, and
+the final value of `trJ` is the probability of ``u(t)``.
+
+Example:
+
+```julia
+using ModelingToolkit, OrdinaryDiffEq, Test
+
+@parameters t α β γ δ
+@variables x(t) y(t)
+@derivatives D'~t
+
+eqs = [D(x) ~ α*x - β*x*y,
+       D(y) ~ -δ*y + γ*x*y]
+
+sys = ODESystem(eqs)
+sys2 = liouville_transform(sys)
+@variables trJ
+
+u0 = [x => 1.0,
+      y => 1.0,
+      trJ => 1.0]
+
+prob = ODEProblem(sys2,u0,tspan,p)
+sol = solve(prob,Tsit5())
+```
+
+Where `sol[3,:]` is the evolution of `trJ` over time.
+
+Sources:
+
+Probabilistic Robustness Analysis of F-16 Controller Performance: An
+Optimal Transport Approach
+
+Abhishek Halder, Kooktae Lee, and Raktim Bhattacharya
+https://abhishekhalder.bitbucket.io/F16ACC2013Final.pdf
+"""
+function liouville_transform(sys)
+      t = independent_variable(sys)
+      @variables trJ
+      D = ModelingToolkit.Differential(t)
+      neweq = D(trJ) ~ trJ*-tr(calculate_jacobian(sys))
+      neweqs = [equations(sys);neweq]
+      vars = [states(sys);trJ]
+      ODESystem(neweqs,t,vars,parameters(sys))
+end

test/basic_transformations.jl

@@ -0,0 +1,33 @@
+using ModelingToolkit, OrdinaryDiffEq, Test
+
+@parameters t α β γ δ
+@variables x(t) y(t)
+@derivatives D'~t
+
+eqs = [D(x) ~ α*x - β*x*y,
+       D(y) ~ -δ*y + γ*x*y]
+
+sys = ODESystem(eqs)
+
+u0 = [x => 1.0,
+      y => 1.0]
+
+p  = [α => 1.5,
+      β => 1.0,
+      δ => 3.0,
+      γ => 1.0]
+
+tspan = (0.0,10.0)
+prob = ODEProblem(sys,u0,tspan,p)
+sol = solve(prob,Tsit5())
+
+sys2 = liouville_transform(sys)
+@variables trJ
+
+u0 = [x => 1.0,
+      y => 1.0,
+      trJ => 1.0]
+
+prob = ODEProblem(sys2,u0,tspan,p,jac=true)
+sol = solve(prob,Tsit5())
+@test sol[end,end] ≈ 1.0742818931017244