feat: add unified System type

src/ModelingToolkit.jl

@@ -127,7 +127,6 @@ TODO
 abstract type AbstractSystem end
 abstract type AbstractTimeDependentSystem <: AbstractSystem end
 abstract type AbstractTimeIndependentSystem <: AbstractSystem end
-abstract type AbstractODESystem <: AbstractTimeDependentSystem end
 abstract type AbstractMultivariateSystem <: AbstractSystem end
 abstract type AbstractOptimizationSystem <: AbstractTimeIndependentSystem end
 abstract type AbstractDiscreteSystem <: AbstractTimeDependentSystem end
@@ -162,26 +161,16 @@ include("systems/codegen_utils.jl")
 include("systems/problem_utils.jl")
 include("linearization.jl")
 
-include("systems/optimization/constraints_system.jl")
-include("systems/optimization/optimizationsystem.jl")
 include("systems/optimization/modelingtoolkitize.jl")
 
-include("systems/nonlinear/nonlinearsystem.jl")
 include("systems/nonlinear/homotopy_continuation.jl")
-include("systems/diffeqs/odesystem.jl")
-include("systems/diffeqs/sdesystem.jl")
 include("systems/diffeqs/abstractodesystem.jl")
 include("systems/nonlinear/modelingtoolkitize.jl")
 include("systems/nonlinear/initializesystem.jl")
 include("systems/diffeqs/first_order_transform.jl")
 include("systems/diffeqs/modelingtoolkitize.jl")
 include("systems/diffeqs/basic_transformations.jl")
 
-include("systems/discrete_system/discrete_system.jl")
-include("systems/discrete_system/implicit_discrete_system.jl")
-
-include("systems/jumps/jumpsystem.jl")
-
 include("systems/pde/pdesystem.jl")
 
 include("systems/sparsematrixclil.jl")
@@ -225,7 +214,7 @@ const D = Differential(t)
 PrecompileTools.@compile_workload begin
     using ModelingToolkit
     @variables x(ModelingToolkit.t_nounits)
-    @named sys = ODESystem([ModelingToolkit.D_nounits(x) ~ -x], ModelingToolkit.t_nounits)
+    @named sys = System([ModelingToolkit.D_nounits(x) ~ -x], ModelingToolkit.t_nounits)
     prob = ODEProblem(structural_simplify(sys), [x => 30.0], (0, 100), [], jac = true)
     @mtkmodel __testmod__ begin
         @constants begin
@@ -261,16 +250,14 @@ export AbstractTimeDependentSystem,
        AbstractTimeIndependentSystem,
        AbstractMultivariateSystem
 
-export ODESystem,
-       ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system,
+export ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system,
        add_accumulations, System
 export DAEFunctionExpr, DAEProblemExpr
-export SDESystem, SDEFunction, SDEFunctionExpr, SDEProblemExpr
+export SDEFunction, SDEFunctionExpr, SDEProblemExpr
 export SystemStructure
-export DiscreteSystem, DiscreteProblem, DiscreteFunction, DiscreteFunctionExpr
-export ImplicitDiscreteSystem, ImplicitDiscreteProblem, ImplicitDiscreteFunction,
+export DiscreteProblem, DiscreteFunction, DiscreteFunctionExpr
+export ImplicitDiscreteProblem, ImplicitDiscreteFunction,
        ImplicitDiscreteFunctionExpr
-export JumpSystem
 export ODEProblem, SDEProblem
 export NonlinearFunction, NonlinearFunctionExpr
 export NonlinearProblem, NonlinearProblemExpr
@@ -279,7 +266,6 @@ export IntervalNonlinearProblem, IntervalNonlinearProblemExpr
 export OptimizationProblem, OptimizationProblemExpr, constraints
 export SteadyStateProblem, SteadyStateProblemExpr
 export JumpProblem
-export NonlinearSystem, OptimizationSystem, ConstraintsSystem
 export alias_elimination, flatten
 export connect, domain_connect, @connector, Connection, AnalysisPoint, Flow, Stream,
        instream

src/inputoutput.jl

@@ -161,7 +161,7 @@ has_var(ex, x) = x ∈ Set(get_variables(ex))
 
 """
     (f_oop, f_ip), x_sym, p_sym, io_sys = generate_control_function(
-            sys::AbstractODESystem,
+            sys::System,
             inputs             = unbound_inputs(sys),
             disturbance_inputs = nothing;
             implicit_dae       = false,
@@ -193,7 +193,7 @@ t = 0
 f[1](x, inputs, p, t)
 ```
 """
-function generate_control_function(sys::AbstractODESystem, inputs = unbound_inputs(sys),
+function generate_control_function(sys::AbstractSystem, inputs = unbound_inputs(sys),
         disturbance_inputs = disturbances(sys);
         disturbance_argument = false,
         implicit_dae = false,
@@ -344,7 +344,7 @@ The structure represents a model of a disturbance, along with the input variable
 # Fields:
 
   - `input`: The variable affected by the disturbance.
-  - `model::M`: A model of the disturbance. This is typically an `ODESystem`, but type that implements [`ModelingToolkit.get_disturbance_system`](@ref)`(dist::DisturbanceModel) -> ::ODESystem` is supported.
+  - `model::M`: A model of the disturbance. This is typically a `System`, but type that implements [`ModelingToolkit.get_disturbance_system`](@ref)`(dist::DisturbanceModel) -> ::System` is supported.
 """
 struct DisturbanceModel{M}
     input::Any
@@ -354,7 +354,7 @@ end
 DisturbanceModel(input, model; name) = DisturbanceModel(input, model, name)
 
 # Point of overloading for libraries, e.g., to be able to support disturbance models from ControlSystemsBase
-function get_disturbance_system(dist::DisturbanceModel{<:ODESystem})
+function get_disturbance_system(dist::DisturbanceModel{System})
     dist.model
 end
 
@@ -395,7 +395,7 @@ c = 10   # Damping coefficient
 eqs = [connect(torque.flange, inertia1.flange_a)
        connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
        connect(inertia2.flange_a, spring.flange_b, damper.flange_b)]
-model = ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper],
+model = System(eqs, t; systems = [torque, inertia1, inertia2, spring, damper],
                   name = :model)
 model = complete(model)
 model_outputs = [model.inertia1.w, model.inertia2.w, model.inertia1.phi, model.inertia2.phi]
@@ -427,7 +427,7 @@ function add_input_disturbance(sys, dist::DisturbanceModel, inputs = nothing; kw
 
     eqs = [dsys.input.u[1] ~ d
            dist.input ~ u + dsys.output.u[1]]
-    augmented_sys = ODESystem(eqs, t, systems = [dsys], name = gensym(:outer))
+    augmented_sys = System(eqs, t, systems = [dsys], name = gensym(:outer))
     augmented_sys = extend(augmented_sys, sys)
 
     (f_oop, f_ip), dvs, p, io_sys = generate_control_function(augmented_sys, all_inputs,

src/linearization.jl

@@ -19,7 +19,7 @@ The `simplified_sys` has undergone [`structural_simplify`](@ref) and had any occ
 
 # Arguments:
 
-  - `sys`: An [`ODESystem`](@ref). This function will automatically apply simplification passes on `sys` and return the resulting `simplified_sys`.
+  - `sys`: A [`System`](@ref) of ODEs. This function will automatically apply simplification passes on `sys` and return the resulting `simplified_sys`.
   - `inputs`: A vector of variables that indicate the inputs of the linearized input-output model.
   - `outputs`: A vector of variables that indicate the outputs of the linearized input-output model.
   - `simplify`: Apply simplification in tearing.
@@ -640,15 +640,15 @@ function plant(; name)
     @variables u(t)=0 y(t)=0
     eqs = [D(x) ~ -x + u
            y ~ x]
-    ODESystem(eqs, t; name = name)
+    System(eqs, t; name = name)
 end
 
 function ref_filt(; name)
     @variables x(t)=0 y(t)=0
     @variables u(t)=0 [input = true]
     eqs = [D(x) ~ -2 * x + u
            y ~ x]
-    ODESystem(eqs, t, name = name)
+    System(eqs, t, name = name)
 end
 
 function controller(kp; name)
@@ -657,7 +657,7 @@ function controller(kp; name)
     eqs = [
         u ~ kp * (r - y),
     ]
-    ODESystem(eqs, t; name = name)
+    System(eqs, t; name = name)
 end
 
 @named f = ref_filt()
@@ -668,7 +668,7 @@ connections = [f.y ~ c.r # filtered reference to controller reference
                c.u ~ p.u # controller output to plant input
                p.y ~ c.y]
 
-@named cl = ODESystem(connections, t, systems = [f, c, p])
+@named cl = System(connections, t, systems = [f, c, p])
 
 lsys0, ssys = linearize(cl, [f.u], [p.x])
 desired_order = [f.x, p.x]

src/problems/bvproblem.jl

@@ -0,0 +1,90 @@
+"""
+```julia
+SciMLBase.BVProblem{iip}(sys::AbstractSystem, u0map, tspan,
+                         parammap = DiffEqBase.NullParameters();
+                         constraints = nothing, guesses = nothing,
+                         version = nothing, tgrad = false,
+                         jac = true, sparse = true,
+                         simplify = false,
+                         kwargs...) where {iip}
+```
+
+Create a boundary value problem from the [`System`](@ref). 
+
+`u0map` is used to specify fixed initial values for the states. Every variable 
+must have either an initial guess supplied using `guesses` or a fixed initial 
+value specified using `u0map`.
+
+Boundary value conditions are supplied to Systems in the form of a list of constraints.
+These equations  should specify values that state variables should take at specific points,
+as in `x(0.5) ~ 1`). More general constraints that  should hold over the entire solution,
+such as `x(t)^2 + y(t)^2`, should be  specified as one of the equations used to build the
+`System`.
+
+If a `System` without `constraints` is specified, it will be treated as an initial value problem. 
+
+```julia
+    @parameters g t_c = 0.5
+    @variables x(..) y(t) λ(t)
+    eqs = [D(D(x(t))) ~ λ * x(t)
+           D(D(y)) ~ λ * y - g
+           x(t)^2 + y^2 ~ 1]
+    cstr = [x(0.5) ~ 1]
+    @mtkbuild pend = System(eqs, t; constraints = cstrs)
+
+    tspan = (0.0, 1.5)
+    u0map = [x(t) => 0.6, y => 0.8]
+    parammap = [g => 1]
+    guesses = [λ => 1]
+
+    bvp = SciMLBase.BVProblem{true, SciMLBase.AutoSpecialize}(pend, u0map, tspan, parammap; guesses, check_length = false)
+```
+
+If the `System` has algebraic equations, like `x(t)^2 + y(t)^2`, the resulting 
+`BVProblem` must be solved using BVDAE solvers, such as Ascher.
+"""
+@fallback_iip_specialize function SciMLBase.BVProblem{iip, spec}(
+        sys::System, u0map, tspan, parammap = SciMLBase.NullParameters();
+        check_compatibility = true, cse = true, checkbounds = false, eval_expression = false,
+        eval_module = @__MODULE__, guesses = Dict(), callback = nothing, kwargs...) where {
+        iip, spec}
+    check_complete(sys, BVProblem)
+    check_compatibility && check_compatible_system(BVProblem, sys)
+    isnothing(callback) || error("BVP solvers do not support callbacks.")
+
+    # Systems without algebraic equations should use both fixed values + guesses
+    # for initialization.
+    _u0map = has_alg_eqs(sys) ? u0map : merge(Dict(u0map), Dict(guesses))
+    fode, u0, p = process_SciMLProblem(
+        ODEFunction{iip, spec}, sys, _u0map, parammap; guesses,
+        t = tspan !== nothing ? tspan[1] : tspan, check_compatibility = false, cse, checkbounds,
+        time_dependent_init = false, kwargs...)
+
+    dvs = unknowns(sys)
+    stidxmap = Dict([v => i for (i, v) in enumerate(dvs)])
+    u0_idxs = has_alg_eqs(sys) ? collect(1:length(dvs)) : [stidxmap[k] for (k, v) in u0map]
+    fbc = generate_boundary_conditions(
+        sys, u0, u0_idxs, tspan; expression = Val{false}, cse, checkbounds)
+
+    if (length(constraints(sys)) + length(u0map) > length(dvs))
+        @warn "The BVProblem is overdetermined. The total number of conditions (# constraints + # fixed initial values given by u0map) exceeds the total number of states. The BVP solvers will default to doing a nonlinear least-squares optimization."
+    end
+
+    kwargs = process_kwargs(sys; kwargs...)
+    # Call `remake` so it runs initialization if it is trivial
+    return remake(BVProblem{iip}(fode, fbc, u0, tspan[1], p; kwargs...))
+end
+
+function check_compatible_system(T::Type{BVProblem}, sys::System)
+    check_time_dependent(sys, T)
+    check_not_dde(sys)
+    check_no_cost(sys, T)
+    check_has_constraints(sys, T)
+    check_no_jumps(sys, T)
+    check_no_noise(sys, T)
+    check_is_continuous(sys, T)
+
+    if !isempty(discrete_events(sys)) || !isempty(continuous_events(sys))
+        throw(SystemCompatibilityError("BVP solvers do not support events."))
+    end
+end