fix liouville stuff

Author: ChrisRackauckas

src/ModelingToolkit.jl

@@ -207,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

@@ -5,12 +5,24 @@ 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 and should be used with a zero initial condition.
+
+
+
+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) ~ -tr(calculate_jacobian(sys))
+      neweq = D(trJ) ~ trJ*-tr(calculate_jacobian(sys))
       neweqs = [equations(sys);neweq]
       vars = [states(sys);trJ]
       ODESystem(neweqs,t,vars,parameters(sys))

test/basic_transformations.jl

@@ -1,38 +1,33 @@
 using ModelingToolkit, OrdinaryDiffEq, Test
 
-@parameters t σ ρ β
-@variables x(t) y(t) z(t)
+@parameters t α β γ δ
+@variables x(t) y(t)
 @derivatives D'~t
 
-eqs = [D(D(x)) ~ σ*(y-x),
-       D(y) ~ x*(ρ-z)-y,
-       D(z) ~ x*y - β*z]
+eqs = [D(x) ~ α*x - β*x*y,
+       D(y) ~ -δ*y + γ*x*y]
 
 sys = ODESystem(eqs)
-sys = ode_order_lowering(sys)
 
-u0 = [D(x) => 2.0,
-      x => 1.0,
-      y => 0.0,
-      z => 0.0]
+u0 = [x => 1.0,
+      y => 1.0]
 
-p  = [σ => 28.0,
-      ρ => 10.0,
-      β => 8/3]
+p  = [α => 1.5,
+      β => 1.0,
+      δ => 3.0,
+      γ => 1.0]
 
-tspan = (0.0,100.0)
-prob = ODEProblem(sys,u0,tspan,p,jac=true)
+tspan = (0.0,10.0)
+prob = ODEProblem(sys,u0,tspan,p)
 sol = solve(prob,Tsit5())
 
-sys2 = louiville_transform(sys)
+sys2 = liouville_transform(sys)
 @variables trJ
 
-u0 = [D(x) => 2.0,
-      x => 1.0,
-      y => 0.0,
-      z => 0.0,
-      trJ => 0.0]
+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] ≈ 366.66666666
+@test sol[end,end] ≈ 1.0742818931017244