GitHub

System construction

Internal API

The content on this page describes internal implementation details of ModelingToolkit. These APIs may change without notice.

A System is the fundamental data structure in MTK. It is a hierarchical data structure which contains all symbolic information about a particular numerical problem. A System in general does not constrain what sorts of information can and cannot be provided about the system. It merely enforces invariants required for the construction of a valid system. Different problem constructors validate that the system they receive is valid for generating code for that type of numerical problem.

Time-dependent and independent systems

The most broad categorization of a system is time-dependence. A system may represent dynamics that vary over time (or any independent variable in general) or a time-independent problem. Time-dependent systems include those for ODEs, SDEs, Jump processes, etc. Time-independent ones include nonlinear systems, optimization systems, and so on. Programmatically, this distinction is made via the iv::Union{SymbolicT, Nothing} field of the system. Union-splitting or type-asserting this field access is highly recommended. For example, if a function knows it operates on time-dependent systems, it should use get_iv(sys)::SymbolicT when accessing the field to improve type-stability and the quality of the code generated by the Julia compiler. In general, these sorts of type-assertions on unstable fields are very useful for performant symbolic-numeric code.

Contracts of the type

The inner constructor for System takes the fields in order. While the order of fields doesn't really matter, it is useful to group them along with similar fields in the struct. All fields should be given as concrete of a type as possible without making System parametric.

The system should never store wrapped types. Any symbolic expression anywhere in the system should always be of type SymbolicT. Storing Num, Arr{T, N}, etc. is invalid. In general, SymbolicT is a useful type-alias. Due to the Const variant of BasicSymbolic, it is often useful to store values of various different types and sizes as a SymbolicT.

All fields should have a docstring. Fields which are only intended for internal access should be annotated with $INTERNAL_FIELD_WARNING at the top of the docstring. Note that "internal" here doesn't mean only ModelingToolkitBase. Internal fields can also be used by ModelingToolkitTearing and ModelingToolkit.

Immutability and updating fields

Immutability of `System`

System should only ever be an immutable type with no type-parameters. Introducing type-parameters makes ensuring type-stability of most symbolic passes very difficult. In the future, if the struct grows large enough that passing it to functions is a significant slowdown (evidenced by benchmarks or Julia's lowered IR) it can be made a mutable struct with all const fields. This is imperative for ensuring correct behavior.

Almost all fields of System should be considered immutable. Internal mutation is also disallowed, with the specific exceptions of get_initial_conditions(sys) and get_guesses(sys). Outside of these two fields, no other field may be mutated. To change fields of a system, use Setfield.@set or Setfield.@set!. This may be changed in the future to Accessors.@set and Accessors.@reset respectively, without breaking backward compatibility. For example, to add a parameter p to a system, do

sys # The system
ps = copy(get_ps(sys)) # note the `copy`
push!(ps, p)
@set! sys.ps = ps

In case multiple such fields need to be set in quick succession, using ConstructionBase.setproperties directly generates better code after compilation. This function is implemented for AbstractSystem as an @generated function which performs additional validation/updates and is what @set/@set! eventually call.

Constructors

The inner constructor of System optionally performs some basic validation fundamental to the type, and builds the struct. The checks keyword is used to control what validation is performed. The checks here are not allowed to mutate the fields, and should not be unduly expensive to perform. This constructor is only ever called internally, and is not user-facing.

The main user-facing constructor is System(eqs::Vector{Equation}, iv, dvs, ps, brownians = SymbolicT[]; kws...). Most fields of System are provided as keyword arguments. This constructor performs significant additional validation.

!!! note __legacy_defaults__ This keyword arguments only exists for the now-deprecated @mtkmodel. It should not be used by anything else.

This constructor is responsible for:

  • Unwrapping all wrapped symbolic variables/expressions passed to it.
  • Coercing all containers to the appropriate type.
  • Populating var_to_name.
  • Validating bindings.
  • Parsing bindings/initial conditions/guesses from variable metadata.
  • Converting symbolic event shorthand into SymbolicContinousCallback/SymbolicDiscreteCallback.
Parsing metadata

It is often useful to disable parsing bindings/initial conditions/guesses from variable metadata during system construction. This is usually when a new system is built from an existing one. Passing discover_from_metadata = false disables the parsing.

Other constructors are shorthand that call the main constructor. Notably, System(eqs, t) automatically discovers variables and parameters from the equations (and other information) passed to the constructor. Similar behavior exists for the time-independent constructors. The variable discovery process is very carefully crafted. Delays, BVP constraints, etc. all reference variables in non-standard ways which require precise handling.

System equality

System uses the default fallback of isequal, which is to compare referential equality. This is because checking if two expressions are equivalent is undecidable in general. In addition, it is very expensive to check equality in a way that is agnostic to variable/equation ordering. Base.isapprox is provided to check this sort of equality, in case it is necessary.

Base.isapprox — Method
isapprox(sysa::System, sysb::System) -> Bool

Check if two systems are about equal, to the extent that ModelingToolkitBase.jl supports. Note that if this returns true, the systems are not guaranteed to be exactly equivalent (unless sysa === sysb) but are highly likely to represent a similar mathematical problem. If this returns false, the systems are very likely to be different.

source

Metadata API

Systems allow storing and accessing metadata using the same API as symbolic expressions. SymbolicUtils.setmetadata, SymbolicUtils.getmetadata, SymbolicUtils.hasmetadata are implemented appropriately. Note that setmetadata returns a new system with the added metadata, since Systems are immutable (similar to the behavior for expressions).

In addition, System provides a caching API. This cache is mutable, and storing values here does not return a new system. This cache is currently only for internal use. It is highly useful for caching expensive intermediate results, such as symbolic jacobians.

ModelingToolkitBase.check_mutable_cache — Function
check_mutable_cache(
    sys::System,
    K::DataType,
    _::Type{V},
    default
) -> Any

If the mutable cache of sys contains an entry for key K return it. Otherwise, return default. V is the expected type of the cache entry, if it exists.

source
ModelingToolkitBase.should_invalidate_mutable_cache_entry — Function
should_invalidate_mutable_cache_entry(
    K::DataType,
    patch::NamedTuple
) -> Bool

Return a Bool indicating whether the cached entry associated with key K should be invalidated if the given patch is applied to the system. By default, entries will be invalidated on all patches.

source

These functions are only intended to be used on flattened systems. Behavior is undefined when used on non-flattened systems.

Preparation for codegen

Systems have additional constraints before codegen can be done. For differential-equation-based systems (ODEs, DAEs, SDEs, etc.) the equations must either be differential equations of the form

D(x) ~ rhs

or algebraic equations of the form

0 ~ rhs

Here, rhs is an expression that cannot involve operator applications such as D(x). rhs may involve observed variables, in which case the required observed equations will be present during code generation.

In addition, all systems must be flattened and marked as complete via complete. The function will not validate the above structure of equations. It performs the following operations:

  • Discover all GlobalScope variables in the system.
  • Expand all connect equations.
  • Flatten the hierarchical system.
  • Convert discrete variables (@discretes) to parameters.
  • Construct ParameterBindingsGraph.
  • Remove bound parameters from the list of parameters
  • Validate the absence of any equations involving only parameters.
  • Add Initial parameters to the system.
  • Appropriately prepare all symbolic callbacks.
  • If split = true (recommended) creates the IndexCache which contains indexing and structure information for parameters. Also reorders all parameters as required by the data structure.
  • Disabling namespacing for the top-level system.

Note that this list may change between versions. In general, systems are allowed to perform any transformation and caching during complete which does not modify the intent of the system.

Initial parameters

ModelingToolkit inserts special Initial parameters during complete. These are an important utility for correctly building the initialization system during codegen. They are almost always invisible to the user. They do not show up in parameters(sys), but can be obtained via get_ps(sys). They are stored in the .initials field of MTKParameters.

Initial parameters exist for every unknown and observed. For time dependent systems, Initial parameters are also added for the first derivative of each unknown and observed variable. All floating-point discrete variables (ModelingToolkitBase.is_variable_floatingpoint) also have an Initial parameter, as do floating-point bound parameters.

Observed equations must be topologically sorted. This is also required for full_equations to work correctly.

SDEs must be phrased in terms of noise_eqs and not brownians. No brownian variables should be present in the equations.

The order and number of unknowns should correspond to that of equations. For example, if the ith equation is

D(x) ~ rhs

Then x should be the ith unknown. Algebraic variables often have a weaker correspondance.

AtomicArrayDict

ModelingToolkit stores variable-value mappings using AtomicArrayDict. It is a custom AbstractDict wrapping any other AbstractDict, and disallows scalarized/indexed arrays from being keys. Initial conditions, guesses and bindings are stored as this type (or wrappers thereof).

Several functions exist for manipulating AtomicArrayDict.

ModelingToolkitBase.AtomicArrayDict — Type
struct AtomicArrayDict{V, D<:AbstractDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, V}} <: AbstractDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, V}

Wrapper over an AbstractDict{SymbolicT, V} where {V} which disallows keys that are indexed array symbolics. Specifically, if @variables x[1:4] exists, then x can be a key but x[1] cannot.

source
ModelingToolkitBase.as_atomic_dict_with_defaults — Function
as_atomic_dict_with_defaults(
    dict::AbstractDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}},
    default::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}
) -> ModelingToolkitBase.AtomicArrayDict

Convert the symbolic mapping dict to an AtomicArrayDict. If dict contains keys which are elements of a symbolic array, the returned mappng will have a key for the array, and a value which is a symbolic array where entries specified in dict are present and default otherwise.

source
ModelingToolkitBase.write_possibly_indexed_array! — Function
write_possibly_indexed_array!(
    dd::ModelingToolkitBase.AtomicArrayDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, D} where D<:AbstractDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}},
    k::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal},
    v::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal},
    default::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}
) -> ModelingToolkitBase.AtomicArrayDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, D} where D<:AbstractDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}

Modify an atomic array mapping dd to map k to v. If k is an indexed array symbolic, update the array to have value v at the corresponding index. If the array is not a key, create the key and set all other entries to default.

source
ModelingToolkitBase.has_possibly_indexed_key — Function
has_possibly_indexed_key(
    dd::ModelingToolkitBase.AtomicArrayDict,
    k::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}
) -> Any

Check if dd has the key k. If k is indexed, check if dd has the array as a key.

source
ModelingToolkitBase.get_possibly_indexed — Function
get_possibly_indexed(
    dd::ModelingToolkitBase.AtomicArrayDict,
    k::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal},
    default
) -> Any

Equivalent to get(dd, k, default). If k is an indexed array, then return dd[arr][idxs...] for the corresponding array arr and indices, or default if arr does not exist.

source

In addition, AtomicArraySet wraps AtomicArrayDict to expose an AbstractSet which cannot contain scalarized/indexed arrays.

Missing docstring.

Missing docstring for ModelingToolkitBase.AtomicArraySet. Check Documenter's build log for details.

ModelingToolkitBase.push_as_atomic_array! — Function
push_as_atomic_array!(
    x::ModelingToolkitBase.AtomicArraySet,
    item::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}
) -> ModelingToolkitBase.AtomicArraySet

Add item to x. If item is an indexed array, add the array instead.

source
Missing docstring.

Missing docstring for ModelingToolkitBase.contains_possibly_indexed_element. Check Documenter's build log for details.

When using an AtomicArrayDict as a list of substitution rules, it is typically necessary to wrap it in AtomicArrayDictSubstitutionWrapper.

ModelingToolkitBase.AtomicArrayDictSubstitutionWrapper — Type
struct AtomicArrayDictSubstitutionWrapper{D} <: AbstractDict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}

A wrapper around AtomicArrayDict{SymbolicT, D} which makes it more amenable to substitution. This wrapper allows indexing with indexed arrays, invoking get_possibly_indexed, has_possibly_indexed_key and write_possibly_indexed_array! as appropriate. More significantly, when Base.get (and Base.getindex) is called with an indexed array for which the corresponding array is present in the wrapped dict but the value is default, it will return the indexed array. For example, if the wrapped dict has k => [default, other_var] as a key-value pair and this wrapper is accessed at k[1], it will return k[1] instead of default. This is useful since default is used to represent missing values. Substituting something like sin(k[1]) will then not attempt to perform sin(default) but instead return sin(k[1]).

source

Variable utilities

To get all the discrete variables of a system, use ModelingToolkitBase.get_all_discretes. If the system has been marked as complete, use ModelingToolkitBase.get_all_discretes_fast instead.

ModelingToolkitBase.get_all_discretes — Function
get_all_discretes(
    sys::ModelingToolkitBase.AbstractSystem
) -> ModelingToolkitBase.AtomicArraySet{Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, Nothing}}

Find all discretes from symbolic event affects.

source
ModelingToolkitBase.get_all_discretes_fast — Function
get_all_discretes_fast(
    sys::ModelingToolkitBase.AbstractSystem
) -> ModelingToolkitBase.AtomicArraySet{Dict{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, Nothing}}

Identical to get_all_discretes, except it uses the IndexCache as a fast path if possible.

source

For array variables, it is often necessary to obtain the parent array from a scalarized element. This should be done using ModelingToolkitBase.split_indexed_var. To obtain a SymbolicUtils.StableIndex{Int} representing the index the scalarized variable accesses the array at, use ModelingToolkitBase.get_stable_index.

ModelingToolkitBase.split_indexed_var — Function
split_indexed_var(
    x::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}
) -> Tuple{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}, Bool}

Given a symbolic variable x, check whether it is an indexed array symbolic. If it is, return the array and true. Otherwise, return x, false.

source
ModelingToolkitBase.get_stable_index — Function
get_stable_index(
    x::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}
) -> SymbolicUtils.StableIndex{Int64}

Given a symbolic variable x, assume split_indexed_var(x)[2] is true. Return the corresponding SymbolicUtils.StableIndex.

source

Best practices

struct Schedule is defined alongside System for convenience, since it is stored in the system by ModelingToolkitTearing when running mtkcompile and contains useful information about the simplified system.

Always use @unpack to obtain variables/parameters/subsystems from a system you intend to extend.

Always use @named to call component constructors. All component constructors should use the @component or @connector macros as appropriate.

In case of conflicts, initial conditions and guesses from parent systems take priority over those in children.

Any container of symbolic variables/expressions should use SymbolicT as the corresponding type. Functions should be as concretely typed as possible, even to the extent of restricting inputs to exact types such as Vector{SymbolicT}. This helps improve type-stability, reduce invalidations, and improve precompilation. This is all the more true of internal functions. Such practice also helps identify bugs when unsupported types are passed into functions, or the impacts of changing a container type when it is passed downstream.