SciML / ModelingToolkit.jl

An acausal modeling framework for automatically parallelized scientific machine learning (SciML) in Julia. A computer algebra system for integrated symbolics for physics-informed machine learning and automated transformations of differential equations
https://mtk.sciml.ai/dev/
Other
1.42k stars 206 forks source link

Nonlinear problems does not seem to follow the same interface as ODESystems #1074

Closed JKRT closed 3 years ago

JKRT commented 3 years ago

Hi!

A wall of text again from me:

How do I construct a system without differential terms? A regular system can be defined like this:

    using ModelingToolkit
    using DiffEqBase
    using DifferentialEquations
    function Model(tspan = (0.0, 1.0))
        pars = ModelingToolkit.@parameters(begin
            (a, t)
        end)
        vars = ModelingToolkit.@variables(begin 
            (x(t),)
        end)
        der = Differential(t)
        eqs = [begin
            der(x) ~ -a * x
        end]
        nonLinearSystem =
            ModelingToolkit.ODESystem(eqs, t, vars, pars, name = :($(Symbol("Model"))))
        pars = Dict(begin
            a => float(begin
                1.0
            end)
        end, t => tspan[1])
        initialValues = [begin
            x => begin
                1.0
            end
        end]
        firstOrderSystem = ModelingToolkit.ode_order_lowering(nonLinearSystem)
        reducedSystem = ModelingToolkit.dae_index_lowering(firstOrderSystem)
        problem = ModelingToolkit.ODEProblem(reducedSystem, initialValues, tspan, pars)
        return problem
    end
   problem = Model()
    function Simulate(tspan = (0.0, 1.0))
        solve(problem, tspan = tspan)
    end
    function Simulate(tspan = (0.0, 1.0); solver = Tsit5())
        solve(problem, tspan = tspan, solver)
    end
Simulate()

Result:


retcode: Success
Interpolation: specialized 4th order "free" interpolation
t: 5-element Vector{Float64}:
 0.0
 0.10001999200479662
 0.3420842593371956
 0.6553980113778367
 1.0
u: 5-element Vector{Vector{Float64}}:
 [1.0]
 [0.9048193287657775]
 [0.710288370182418]
 [0.5192354411753066]
 [0.3678795934275483]

However, how do I go about defining a nonlinear model? For instance a model on a straight line depending on time. Below I go about the same way as above but using a Nonlinear problem instead, inspo from https://github.com/SciML/ModelingToolkit.jl/issues/969

However, it seems that Nonlinear systems do not quite follow the same interface as the ODE and DAE problems. Nonsense example below:

using ModelingToolkit
using DiffEqBase
using DifferentialEquations
using NonlinearSolve
function Model(tspan = (0.0, 1.0))
    pars = ModelingToolkit.@parameters(begin
        (a, t)
    end)
    vars = ModelingToolkit.@variables(begin 
        (x,)
    end)
    eqs = [begin
        0 ~ -t * x + 100
    end]
    nonLinearSystem =
        ModelingToolkit.NonlinearSystem(eqs, t, vars, pars, name = :($(Symbol("Model"))))
    pars = Dict(begin
        a => float(begin
            1.0
        end)
    end, t => tspan[1])
    initialValues = [begin
        x => begin
            1.0
        end
    end]
    problem = ModelingToolkit.NonlinearProblem(nonLinearSystem, initialValues, tspan, pars)
    return problem
end
problem = Model()
function Simulate(tspan = (0.0, 1.0))
    solve(problem, tspan = tspan)
end
Simulate()

Result:

ERROR: LoadError: MethodError: no method matching NonlinearSystem(::Vector{Equation}, ::Num, ::Vector{Num}, ::Vector{Num}; name=:Model)
Closest candidates are:
  NonlinearSystem(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any) at 
....

I note that I probably have done something wrong, but the documentation about the non-linear systems currently does not involve systems depending on time.

Cheers. John

lamorton commented 3 years ago

the documentation about the non-linear systems currently does not involve systems depending on time.

You're right, nonlinear systems in this context means 'time-independent systems of equations.' ODESystems can have linear or nonlinear differential equations. The 'nonlinear' here is to distinguish from linear systems of time independent equations' (ie, linear algebra). The error is happening b/c you are supplying too many arguments to the constructor.

Have a look at the tutorial for NonlinearSystems for an example.

JKRT commented 3 years ago

the documentation about the non-linear systems currently does not involve systems depending on time.

You're right, nonlinear systems in this context means 'time-independent systems of equations.' ODESystems can have linear or nonlinear differential equations. The 'nonlinear' here is to distinguish from linear systems of time independent equations' (ie, linear algebra). The error is happening b/c you are supplying too many arguments to the constructor.

Have a look at the tutorial for NonlinearSystems for an example.

Right, thanks for clarifying! I went with what I had above after reading Yingbos comment here: https://github.com/SciML/ModelingToolkit.jl/issues/969

Next try, change it to an ODESystem (However, it has only one linear (algebraic) equation, no differential terms. Which makes it a DAE with the set of differential equations = ∅):

using ModelingToolkit
using DiffEqBase
using DifferentialEquations

function Model(tspan = (0.0, 1.0))
    pars = ModelingToolkit.@parameters(begin
        (a, t)
    end)
    vars = ModelingToolkit.@variables(begin 
        (x(t),)
    end)
    eqs = [begin
        0 ~ -x + t + 1
    end]
    nonLinearSystem =
        ModelingToolkit.ODESystem(eqs, t, vars, pars, name = :($(Symbol("Model"))))
    pars = Dict(begin
        a => float(begin
            1.0
        end)
    end, t => tspan[1])
    initialValues = [begin
        x => begin
            1.0
        end
    end]
    problem = ModelingToolkit.ODEProblem(nonLinearSystem, initialValues, tspan, pars)
    return problem
end
problem = Model()
function Simulate(tspan = (0.0, 1.0))
    solve(problem, tspan = tspan)
end
t: 9-element Vector{Float64}:
  0.0
  1.0e-6
  1.1e-5
  0.00011099999999999999
  0.0011109999999999998
  0.011110999999999996
  0.11111099999999996
  1.1111109999999995
 10.0
u: 9-element Vector{Vector{Float64}}:
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]
 [1.0]

This does not behave as I would expect, for instance in Modelica the corresponding model is written as:

model StraightLine
  Real x;
equation 
  0 = -x + time + 1;
end StraightLine;
lamorton commented 3 years ago

The mistaken output here is because you are putting t in the vector of parameters and supplying an initial value for t. When you do that, it freezes time. This is #1063, which is now fixed. Try updating ModelingToolkit to the latest version, you should get an error message at this line:

ModelingToolkit.ODESystem(eqs, t, vars, pars, name = :($(Symbol("Model"))))

Your fix is to do:

    pars = @parameters((a,))
   @parameters t

so that you don't count time as a parameter.

YingboMa commented 3 years ago

Could you file a PR to document this if it's not documented already?

lamorton commented 3 years ago

I think we can close this now that #1091, #1079, and #1064 are merged.