feat: add unified System type

examples/electrical_components.jl

@@ -4,13 +4,13 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @connector function Pin(; name)
     sts = @variables v(t) [guess = 1.0] i(t) [guess = 1.0, connect = Flow]
-    ODESystem(Equation[], t, sts, []; name = name)
+    System(Equation[], t, sts, []; name = name)
 end
 
 @component function Ground(; name)
     @named g = Pin()
     eqs = [g.v ~ 0]
-    compose(ODESystem(eqs, t, [], []; name = name), g)
+    compose(System(eqs, t, [], []; name = name), g)
 end
 
 @component function OnePort(; name)
@@ -20,7 +20,7 @@ end
     eqs = [v ~ p.v - n.v
            0 ~ p.i + n.i
            i ~ p.i]
-    compose(ODESystem(eqs, t, sts, []; name = name), p, n)
+    compose(System(eqs, t, sts, []; name = name), p, n)
 end
 
 @component function Resistor(; name, R = 1.0)
@@ -30,7 +30,7 @@ end
     eqs = [
         v ~ i * R
     ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+    extend(System(eqs, t, [], ps; name = name), oneport)
 end
 
 @component function Capacitor(; name, C = 1.0)
@@ -40,7 +40,7 @@ end
     eqs = [
         D(v) ~ i / C
     ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+    extend(System(eqs, t, [], ps; name = name), oneport)
 end
 
 @component function ConstantVoltage(; name, V = 1.0)
@@ -50,7 +50,7 @@ end
     eqs = [
         V ~ v
     ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+    extend(System(eqs, t, [], ps; name = name), oneport)
 end
 
 @component function Inductor(; name, L = 1.0)
@@ -60,12 +60,12 @@ end
     eqs = [
         D(i) ~ v / L
     ]
-    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
+    extend(System(eqs, t, [], ps; name = name), oneport)
 end
 
 @connector function HeatPort(; name)
     @variables T(t) [guess = 293.15] Q_flow(t) [guess = 0.0, connect = Flow]
-    ODESystem(Equation[], t, [T, Q_flow], [], name = name)
+    System(Equation[], t, [T, Q_flow], [], name = name)
 end
 
 @component function HeatingResistor(; name, R = 1.0, TAmbient = 293.15, alpha = 1.0)
@@ -79,7 +79,7 @@ end
            h.Q_flow ~ -v * p.i # -LossPower
            v ~ p.v - n.v
            0 ~ p.i + n.i]
-    compose(ODESystem(eqs, t, [v, RTherm], [R, TAmbient, alpha],
+    compose(System(eqs, t, [v, RTherm], [R, TAmbient, alpha],
             name = name), p, n, h)
 end
 
@@ -90,6 +90,6 @@ end
     eqs = [
         D(h.T) ~ h.Q_flow / C
     ]
-    compose(ODESystem(eqs, t, [], [rho, V, cp],
+    compose(System(eqs, t, [], [rho, V, cp],
             name = name), h)
 end

examples/rc_model.jl

@@ -13,5 +13,5 @@ rc_eqs = [connect(source.p, resistor.p)
           connect(capacitor.n, source.n)
           connect(capacitor.n, ground.g)]
 
-@named rc_model = ODESystem(rc_eqs, t)
+@named rc_model = System(rc_eqs, t)
 rc_model = compose(rc_model, [resistor, capacitor, source, ground])

examples/serial_inductor.jl

@@ -12,7 +12,7 @@ eqs = [connect(source.p, resistor.p)
        connect(source.n, inductor2.n)
        connect(inductor2.n, ground.g)]
 
-@named ll_model = ODESystem(eqs, t)
+@named ll_model = System(eqs, t)
 ll_model = compose(ll_model, [source, resistor, inductor1, inductor2, ground])
 
 @named source = ConstantVoltage(V = 10.0)
@@ -29,5 +29,5 @@ eqs = [connect(source.p, inductor1.p)
        connect(resistor2.n, inductor2.p)
        connect(source.n, inductor2.n)
        connect(inductor2.n, ground.g)]
-@named ll2_model = ODESystem(eqs, t)
+@named ll2_model = System(eqs, t)
 ll2_model = compose(ll2_model, [source, resistor1, resistor2, inductor1, inductor2, ground])

src/ModelingToolkit.jl

@@ -127,10 +127,8 @@ 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
 
 function independent_variable end
 
@@ -158,30 +156,40 @@ include("systems/connectors.jl")
 include("systems/analysis_points.jl")
 include("systems/imperative_affect.jl")
 include("systems/callbacks.jl")
+include("systems/system.jl")
 include("systems/codegen_utils.jl")
+include("systems/codegen.jl")
 include("systems/problem_utils.jl")
 include("linearization.jl")
 
-include("systems/optimization/constraints_system.jl")
-include("systems/optimization/optimizationsystem.jl")
+include("problems/compatibility.jl")
+include("problems/odeproblem.jl")
+include("problems/daeproblem.jl")
+include("problems/ddeproblem.jl")
+include("problems/sdeproblem.jl")
+include("problems/sddeproblem.jl")
+include("problems/bvproblem.jl")
+include("problems/discreteproblem.jl")
+include("problems/implicitdiscreteproblem.jl")
+include("problems/nonlinearproblem.jl")
+include("problems/intervalnonlinearproblem.jl")
+include("problems/sccnonlinearproblem.jl")
+include("problems/initializationproblem.jl")
+include("problems/jumpproblem.jl")
+include("problems/optimizationproblem.jl")
+
+include("modelingtoolkitize/common.jl")
+include("modelingtoolkitize/odeproblem.jl")
+include("modelingtoolkitize/sdeproblem.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 +233,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 +269,14 @@ export AbstractTimeDependentSystem,
        AbstractTimeIndependentSystem,
        AbstractMultivariateSystem
 
-export ODESystem,
-       ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system,
-       add_accumulations, System
+export ODEFunction, ODEFunctionExpr, ODEProblemExpr, convert_system,
+       System, OptimizationSystem, JumpSystem, SDESystem
 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 +285,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
@@ -291,19 +296,19 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc
        hasunit, getunit, hasconnect, getconnect,
        hasmisc, getmisc, state_priority
 export ode_order_lowering, dae_order_lowering, liouville_transform,
-       change_independent_variable, substitute_component
+       change_independent_variable, substitute_component, add_accumulations
 export PDESystem
 export Differential, expand_derivatives, @derivatives
 export Equation, ConstrainedEquation
 export Term, Sym
 export SymScope, LocalScope, ParentScope, GlobalScope
-export independent_variable, equations, controls, observed, full_equations
+export independent_variable, equations, controls, observed, full_equations, jumps
 export initialization_equations, guesses, defaults, parameter_dependencies, hierarchy
 export structural_simplify, expand_connections, linearize, linearization_function,
        LinearizationProblem
 export solve
 
-export calculate_jacobian, generate_jacobian, generate_function, generate_custom_function,
+export calculate_jacobian, generate_jacobian, generate_rhs, generate_custom_function,
        generate_W
 export calculate_control_jacobian, generate_control_jacobian
 export calculate_tgrad, generate_tgrad

src/discretedomain.jl

@@ -7,6 +7,8 @@ SymbolicUtils.promote_symtype(::Type{<:SampleTime}, t...) = Real
 Base.nameof(::SampleTime) = :SampleTime
 SymbolicUtils.isbinop(::SampleTime) = false
 
+function validate_operator(op::SampleTime, args, iv; context = nothing) end
+
 # Shift
 
 """
@@ -68,6 +70,13 @@ Base.hash(D::Shift, u::UInt) = hash(D.steps, hash(D.t, xor(u, 0x055640d6d952f101
 Base.:^(D::Shift, n::Integer) = Shift(D.t, D.steps * n)
 Base.literal_pow(f::typeof(^), D::Shift, ::Val{n}) where {n} = Shift(D.t, D.steps * n)
 
+function validate_operator(op::Shift, args, iv; context = nothing)
+    isequal(op.t, iv) || throw(OperatorIndepvarMismatchError(op, iv, context))
+    op.steps <= 0 || error("""
+    Only non-positive shifts are allowed. Found shift of $(op.steps) in $context.
+    """)
+end
+
 hasshift(eq::Equation) = hasshift(eq.lhs) || hasshift(eq.rhs)
 
 """
@@ -128,6 +137,18 @@ Base.show(io::IO, D::Sample) = print(io, "Sample(", D.clock, ")")
 Base.:(==)(D1::Sample, D2::Sample) = isequal(D1.clock, D2.clock)
 Base.hash(D::Sample, u::UInt) = hash(D.clock, xor(u, 0x055640d6d952f101))
 
+function validate_operator(op::Sample, args, iv; context = nothing)
+    arg = unwrap(only(args))
+    if !is_variable_floatingpoint(arg)
+        throw(ContinuousOperatorDiscreteArgumentError(op, arg, context))
+    end
+    if isparameter(arg)
+        throw(ArgumentError("""
+        Expected argument of $op to be an unknown, found $arg which is a parameter.
+        """))
+    end
+end
+
 """
     hassample(O)
 
@@ -156,6 +177,11 @@ SymbolicUtils.isbinop(::Hold) = false
 
 Hold(x) = Hold()(x)
 
+function validate_operator(op::Hold, args, iv; context = nothing)
+    # TODO: maybe validate `VariableTimeDomain`?
+    return nothing
+end
+
 """
     hashold(O)
 

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]