diff --git a/Project.toml b/Project.toml index 41568022ab..527451bca2 100644 --- a/Project.toml +++ b/Project.toml @@ -51,6 +51,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" SCCNonlinearSolve = "9dfe8606-65a1-4bb3-9748-cb89d1561431" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SciMLPublic = "431bcebd-1456-4ced-9d72-93c2757fff0b" SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" Serialization = "9e88b42a-f829-5b0c-bbe9-9e923198166b" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" @@ -145,6 +146,7 @@ Reexport = "0.2, 1" RuntimeGeneratedFunctions = "0.5.9" SCCNonlinearSolve = "1.0.0" SciMLBase = "2.91.1" +SciMLPublic = "1.0.0" SciMLStructures = "1.7" Serialization = "1" Setfield = "0.7, 0.8, 1" diff --git a/docs/src/systems/PDESystem.md b/docs/src/API/PDESystem.md similarity index 100% rename from docs/src/systems/PDESystem.md rename to docs/src/API/PDESystem.md diff --git a/docs/src/API/System.md b/docs/src/API/System.md new file mode 100644 index 0000000000..33a4d064ab --- /dev/null +++ b/docs/src/API/System.md @@ -0,0 +1,174 @@ +# [The `System` type](@id System_type) + +ModelingToolkit.jl uses `System` to symbolically represent all types of numerical problems. +Users create `System`s representing the problem they want to solve and `mtkcompile` transforms +them into a format ModelingToolkit.jl can generate code for (alongside performing other +optimizations). + +```@docs +System +``` + +## Utility constructors + +Several utility constructors also exist to easily construct alternative system formulations. + +```@docs +NonlinearSystem +SDESystem +JumpSystem +OptimizationSystem +``` + +## Accessor functions + +Several accessor functions exist to query systems for the information they contain. In general, +for every field `x` there exists a `has_x` function which checks if the system contains the +field and a `get_x` function for obtaining the value in the field. Note that fields of a system +cannot be accessed via `getproperty` - that is reserved for accessing variables, subsystems +or analysis points of the hierarchical system. + +```@docs +ModelingToolkit.has_eqs +ModelingToolkit.get_eqs +equations +equations_toplevel +full_equations +ModelingToolkit.has_noise_eqs +ModelingToolkit.get_noise_eqs +ModelingToolkit.has_jumps +ModelingToolkit.get_jumps +jumps +ModelingToolkit.has_constraints +ModelingToolkit.get_constraints +constraints +ModelingToolkit.has_costs +ModelingToolkit.get_costs +costs +ModelingToolkit.has_consolidate +ModelingToolkit.get_consolidate +ModelingToolkit.has_unknowns +ModelingToolkit.get_unknowns +unknowns +unknowns_toplevel +ModelingToolkit.has_ps +ModelingToolkit.get_ps +parameters +parameters_toplevel +tunable_parameters +ModelingToolkit.has_brownians +ModelingToolkit.get_brownians +brownians +ModelingToolkit.has_iv +ModelingToolkit.get_iv +ModelingToolkit.has_observed +ModelingToolkit.get_observed +observed +observables +ModelingToolkit.has_name +ModelingToolkit.get_name +nameof +ModelingToolkit.has_description +ModelingToolkit.get_description +ModelingToolkit.description +ModelingToolkit.has_defaults +ModelingToolkit.get_defaults +defaults +ModelingToolkit.has_guesses +ModelingToolkit.get_guesses +guesses +ModelingToolkit.has_initialization_eqs +ModelingToolkit.get_initialization_eqs +initialization_equations +ModelingToolkit.has_continuous_events +ModelingToolkit.get_continuous_events +continuous_events +continuous_events_toplevel +ModelingToolkit.has_discrete_events +ModelingToolkit.get_discrete_events +discrete_events_toplevel +ModelingToolkit.has_assertions +ModelingToolkit.get_assertions +assertions +ModelingToolkit.has_metadata +ModelingToolkit.get_metadata +SymbolicUtils.getmetadata(::ModelingToolkit.AbstractSystem, ::DataType, ::Any) +SymbolicUtils.setmetadata(::ModelingToolkit.AbstractSystem, ::DataType, ::Any) +ModelingToolkit.has_is_dde +ModelingToolkit.get_is_dde +ModelingToolkit.is_dde +ModelingToolkit.has_tstops +ModelingToolkit.get_tstops +ModelingToolkit.symbolic_tstops +ModelingToolkit.has_tearing_state +ModelingToolkit.get_tearing_state +ModelingToolkit.does_namespacing +toggle_namespacing +ModelingToolkit.iscomplete +ModelingToolkit.has_preface +ModelingToolkit.get_preface +ModelingToolkit.preface +ModelingToolkit.has_parent +ModelingToolkit.get_parent +ModelingToolkit.has_initializesystem +ModelingToolkit.get_initializesystem +ModelingToolkit.is_initializesystem +``` + +## `getproperty` syntax + +ModelingToolkit allows obtaining in a system using `getproperty`. For a system `sys` with a +subcomponent `inner` containing variable `var`, `sys.inner.var` will obtain the appropriately +namespaced version of `var`. Note that this can also be used to access subsystems (`sys.inner`) +or analysis points. + +!!! note + + By default, top-level systems not marked as `complete` will apply their namespace. Systems + marked as `complete` will not do this namespacing. This namespacing behavior can be toggled + independently of whether the system is completed using [`toggle_namespacing`](@ref) and the + current namespacing behavior can be queried via [`ModelingToolkit.does_namespacing`](@ref). + +```@docs +Base.getproperty(::ModelingToolkit.AbstractSystem, ::Symbol) +``` + +## Functions for querying system equations + +```@docs +has_diff_eqs +has_alg_eqs +get_diff_eqs +get_alg_eqs +has_diff_equations +has_alg_equations +diff_equations +alg_equations +is_alg_equation +is_diff_equation +``` + +## String parsing + +ModelingToolkit can parse system variables from strings. + +```@docs +ModelingToolkit.parse_variable +``` + +## Dumping system data + +```@docs +ModelingToolkit.dump_unknowns +ModelingToolkit.dump_parameters +``` + +```@docs; canonical = false +ModelingToolkit.dump_variable_metadata +``` + +## Debugging utilities + +```@docs +debug_system +``` diff --git a/docs/src/API/codegen.md b/docs/src/API/codegen.md new file mode 100644 index 0000000000..9cc7c7af4f --- /dev/null +++ b/docs/src/API/codegen.md @@ -0,0 +1,43 @@ +# Code generation utilities + +These are lower-level functions that ModelingToolkit leverages to generate code for +building numerical problems. + +```@docs +ModelingToolkit.generate_rhs +ModelingToolkit.generate_diffusion_function +ModelingToolkit.generate_jacobian +ModelingToolkit.generate_tgrad +ModelingToolkit.generate_hessian +ModelingToolkit.generate_W +ModelingToolkit.generate_dae_jacobian +ModelingToolkit.generate_history +ModelingToolkit.generate_boundary_conditions +ModelingToolkit.generate_cost +ModelingToolkit.generate_cost_gradient +ModelingToolkit.generate_cost_hessian +ModelingToolkit.generate_cons +ModelingToolkit.generate_constraint_jacobian +ModelingToolkit.generate_constraint_hessian +ModelingToolkit.generate_control_jacobian +ModelingToolkit.build_explicit_observed_function +``` + +For functions such as jacobian calculation which require symbolic computation, there +are `calculate_*` equivalents to obtain the symbolic result without building a function. + +```@docs +ModelingToolkit.calculate_tgrad +ModelingToolkit.calculate_jacobian +ModelingToolkit.jacobian_sparsity +ModelingToolkit.jacobian_dae_sparsity +ModelingToolkit.calculate_hessian +ModelingToolkit.hessian_sparsity +ModelingToolkit.calculate_massmatrix +ModelingToolkit.W_sparsity +ModelingToolkit.calculate_W_prototype +ModelingToolkit.calculate_cost_gradient +ModelingToolkit.calculate_cost_hessian +ModelingToolkit.cost_hessian_sparsity +ModelingToolkit.calculate_control_jacobian +``` diff --git a/docs/src/API/model_building.md b/docs/src/API/model_building.md new file mode 100644 index 0000000000..40072bc2c4 --- /dev/null +++ b/docs/src/API/model_building.md @@ -0,0 +1,253 @@ +# [Model building reference](@id model_building_api) + +This page lists functionality and utilities related to building hierarchical models. It is +recommended to read the page on the [`System`](@ref System_type) before this. + +## Hierarchical model composition + +The `System` data structure can represent a tree-like hierarchy of systems for building models +from composable blocks. The [`ModelingToolkit.get_systems`](@ref) function can be used for +querying the subsystems of a system. The `@component` macro should be used when writing +building blocks for model composition. + +```@docs +@component +``` + +Every constructor function should build either a component or a connector. Components define +the dynamics of the system. Connectors are used to connect components together and propagate +information between them. See also [`@connector`](@ref). + +### Scoping of variables + +When building hierarchical systems, is is often necessary to pass variables from a parent system +to the subsystems. If done naively, this will result in the child system assuming it "owns" the +variables passed to it and any occurrences of those variables in the child system will be +namespaced. To prevent this, ModelingToolkit has the concept of variable scope. The scope allows +specifying which system a variable belongs to relative to the system in which it is used. + +```@docs +LocalScope +ParentScope +GlobalScope +``` + +Note that the scopes must be applied to _individual variables_ and not expressions. For example, +`ParentScope(x + y)` is incorrect. Instead, `ParentScope(x) + ParentScope(y)` is the correct usage. +Applying the same scope (more generally, the same function) to all variables in an expression is a +common task, and ModelingToolkit exposes a utility for the same: + +```@docs +ModelingToolkit.apply_to_variables +``` + +It is still tedious to manually use `apply_to_variables` on any symbolic expression passed to a +subsystem. The `@named` macro automatically wraps all symbolic arguments in `ParentScope` and +uses the identifier being assigned as the name of the system. + +```@docs +@named +``` + +### Exploring the tree structure + +The `System` type implements the `AbstractTrees` interface. This can be used to explore the +hierarchical structure. + +```@docs +hierarchy +``` + +### [Connection semantics](@id connect_semantics) + +ModelingToolkit implements connection semantics similar to those in the [Modelica specification](https://specification.modelica.org/maint/3.6/connectors-and-connections.html). +We do not support the concept of `inner` and `outer` elements or `expandable` connectors. +Connectors in ModelingToolkit are systems with the appropriate metadata added via the `@connector` +macro. + +```@docs +connect +domain_connect +@connector +``` + +Connections can be expanded using `expand_connections`. + +```@docs +expand_connections +``` + +Similar to the `stream` and `flow` keyword arguments in the specification, ModelingToolkit +allows specifying how variables in a connector behave in a connection. + +```@docs +Equality +Flow +Stream +``` + +These are specified using the `connect` metadata. ModelingToolkit also supports `instream`. +Refer to the Modelica specification on [Stream connectors](https://specification.modelica.org/maint/3.6/stream-connectors.html) +for more information. + +```@docs +instream +``` + +### System composition utilities + +```@docs +extend +compose +substitute_component +``` + +### Flattening systems + +The hierarchical structure can be flattened. This operation is performed during simplification. + +```@docs +flatten +``` + +## System simplification + +`System`s can be simplified to reformulate them in a way that enables it to be solved numerically, +and also perform other optimizations. This is done via the `mtkcompile` function. Connection expansion +and flattening are preprocessing steps of simplification. + +```@docs +mtkcompile +@mtkcompile +``` + +It is also possible (though not always advisable) to build numerical problems from systems without +passing them through `mtkcompile`. To do this, the system must first be marked as "complete" via +the `complete` function. This process is used to indicate that a system will not be modified +further and allows ModelingToolkit to perform any necessary preprocessing to it. `mtkcompile` +calls `complete` internally. + +```@docs +complete +``` + +### Exploring the results of simplification + +Similar to how [`full_equations`](@ref) returns the equations of a system with all variables +eliminated during `mtkcompile` substituted, we can perform this substitution on an arbitrary +expression. + +```@docs +ModelingToolkit.substitute_observed +ModelingToolkit.empty_substitutions +ModelingToolkit.get_substitutions +``` + +### Experimental simplification + +ModelingToolkit may have a variety of experimental simplification passes. These are not +enabled by default, but can be used by passing to the `additional_passes` keyword argument +of `mtkcompile`. + +```@docs +ModelingToolkit.IfLifting +``` + +## Event handling + +Time-dependent systems may have several events. These are used to trigger discontinuities +in the model. They compile to standard callbacks from `DiffEqCallbacks.jl`. + +```@docs +ModelingToolkit.SymbolicContinuousCallback +ModelingToolkit.SymbolicDiscreteCallback +``` + +The affect functions for the above callbacks can be symbolic or user-defined functions. +Symbolic affects are handled using equations as described in the [Events](@ref events) +section of the documentation. User-defined functions can be used via `ImperativeAffect`. + +```@docs +ModelingToolkit.AffectSystem +ModelingToolkit.ImperativeAffect +``` + +## Modelingtoolkitize + +ModelingToolkit can take some numerical problems created non-symbolically and build a +symbolic representation from them. + +```@docs +modelingtoolkitize +``` + +## Using FMUs + +ModelingToolkit is capable of importing FMUs as black-box symbolic models. Currently only +a subset of FMU features are supported. This functionality requires importing `FMI.jl`. + +```@docs +ModelingToolkit.FMIComponent +``` + +## Model transformations + +ModelingToolkit exposes a variety of transformations that can be applied to models to aid in +symbolic analysis. + +```@docs +liouville_transform +stochastic_integral_transform +Girsanov_transform +change_independent_variable +add_accumulations +noise_to_brownians +convert_system_indepvar +``` + +## Hybrid systems + +Hybrid systems are dynamical systems involving one or more discrete-time subsystems. These +discrete time systems follow clock semantics - they are synchronous systems and the relevant +variables are only defined at points where the clock ticks. + +While ModelingToolkit is unable to simplify, compile and solve such systems on its own, it +has the ability to represent them. Compilation strategies can be implemented independently +on top of [`mtkcompile`](@ref) using the `additional_passes` functionality. + +!!! warn + + These operators are considered experimental API. + +```@docs; canonical = false +Sample +Hold +SampleTime +``` + +ModelingToolkit uses the clock definition in SciMLBase + +```@docs +SciMLBase.TimeDomain +SciMLBase.Clock +SciMLBase.SolverStepClock +SciMLBase.Continuous +``` + +### State machines + +While ModelingToolkit has the capability to represent state machines, it lacks the ability +to compile and simulate them. + +!!! warn + + This functionality is considered experimental API + +```@docs +initial_state +transition +activeState +entry +ticksInState +timeInState +``` diff --git a/docs/src/API/problems.md b/docs/src/API/problems.md new file mode 100644 index 0000000000..5549224a09 --- /dev/null +++ b/docs/src/API/problems.md @@ -0,0 +1,121 @@ +# Building and solving numerical problems + +Systems are numerically solved by building and solving the appropriate problem type. +Numerical solvers expect to receive functions taking a predefeined set of arguments +and returning specific values. This format of argument and return value depends on +the function and the problem. ModelingToolkit is capable of compiling and generating +code for a variety of such numerical problems. + +## Dynamical systems + +```@docs +ODEFunction(::System, args...) +ODEProblem(::System, args...) +DAEFunction(::System, args...) +DAEProblem(::System, args...) +SDEFunction(::System, args...) +SDEProblem(::System, args...) +DDEFunction(::System, args...) +DDEProblem(::System, args...) +SDDEFunction(::System, args...) +SDDEProblem(::System, args...) +JumpProblem(::System, args...) +BVProblem(::System, args...) +DiscreteProblem(::System, args...) +ImplicitDiscreteProblem(::System, args...) +``` + +## Nonlinear systems + +```@docs +NonlinearFunction(::System, args...) +NonlinearProblem(::System, args...) +SCCNonlinearProblem(::System, args...) +NonlinearLeastSquaresProblem(::System, args...) +SteadyStateProblem(::System, args...) +IntervalNonlinearFunction(::System, args...) +IntervalNonlinearProblem(::System, args...) +ModelingToolkit.HomotopyContinuationProblem +HomotopyNonlinearFunction(::System, args...) +``` + +## Optimization and optimal control + +```@docs +OptimizationFunction(::System, args...) +OptimizationProblem(::System, args...) +ODEInputFunction(::System, args...) +JuMPDynamicOptProblem(::System, args...) +InfiniteOptDynamicOptProblem,(::System, args...) +CasADiDynamicOptProblem(::System, args...) +DynamicOptSolution +``` + +## The state vector and parameter object + +Typically the unknowns of the system are present as a `Vector` of the appropriate length +in the numerical problem. The state vector can also be constructed manually without building +a problem. + +```@docs +ModelingToolkit.get_u0 +ModelingToolkit.varmap_to_vars +``` + +By default, the parameters of the system are stored in a custom data structure called +`MTKParameters`. The internals of this data structure are undocumented, and it should +only be interacted with through defined public API. SymbolicIndexingInterface.jl contains +functionality useful for this purpose. + +```@docs +MTKParameters +ModelingToolkit.get_p +``` + +The following functions are useful when working with `MTKParameters` objects, and especially +the `Tunables` portion. For more information about the "portions" of `MTKParameters`, refer +to the [`SciMLStructures.jl`](https://docs.sciml.ai/SciMLStructures/stable/) documentation. + +```@docs +reorder_dimension_by_tunables! +reorder_dimension_by_tunables +``` + +## Initialization + +```@docs +generate_initializesystem +InitializationProblem +``` + +## Linear analysis + +```@docs +linearization_function +LinearizationProblem +linearize +CommonSolve.solve(::LinearizationProblem) +linearize_symbolic +``` + +There are also utilities for manipulating the results of these analyses in a symbolic context. + +```@docs +ModelingToolkit.similarity_transform +ModelingToolkit.reorder_unknnowns +``` + +### Analysis point transformations + +Linear analysis can also be done using analysis points to perform several common +workflows. + +```@docs +get_sensitivity_function +get_sensitivity +get_comp_sensitivity_function +get_comp_sensitivity +get_looptransfer_function +get_looptransfer +open_loop +``` diff --git a/docs/src/basics/Variable_metadata.md b/docs/src/API/variables.md similarity index 63% rename from docs/src/basics/Variable_metadata.md rename to docs/src/API/variables.md index 43d0a7f2ea..0ff2e7799d 100644 --- a/docs/src/basics/Variable_metadata.md +++ b/docs/src/API/variables.md @@ -1,8 +1,22 @@ -# [Symbolic Metadata](@id symbolic_metadata) +# Symbolic variables and variable metadata -It is possible to add metadata to symbolic variables, the metadata will be displayed when calling help on a variable. +ModelingToolkit uses [Symbolics.jl](https://docs.sciml.ai/Symbolics/stable/) for the symbolic +manipulation infrastructure. In fact, the `@variables` macro is defined in Symbolics.jl. In +addition to `@variables`, ModelingToolkit defines `@parameters`, `@independent_variables`, +`@constants` and `@brownians`. These macros function identically to `@variables` but allow +ModelingToolkit to attach additional metadata. -The following information can be added (note, it's possible to extend this to user-defined metadata as well) +```@docs +Symbolics.@variables +@independent_variables +@parameters +@constants +@brownians +``` + +Symbolic variables can have metadata attached to them. The defaults and guesses assigned +at variable construction time are examples of this metadata. ModelingToolkit also defines +additional types of metadata. ## Variable descriptions @@ -39,6 +53,11 @@ help?> u Symbolics.VariableSource: (:variables, :u) ``` +```@docs +hasdescription +getdescription +``` + ## Connect Variables in connectors can have `connect` metadata which describes the type of connections. @@ -61,6 +80,16 @@ hasconnect(i) getconnect(k) ``` +```@docs +hasconnect +getconnect +``` + +```@docs; canonical = false +Flow +Stream +``` + ## Input or output Designate a variable as either an input or an output using the following @@ -78,16 +107,20 @@ isinput(u) isoutput(y) ``` +```@docs +isinput +isoutput +ModelingToolkit.setinput +ModelingToolkit.setoutput +``` + ## Bounds Bounds are useful when parameters are to be optimized, or to express intervals of uncertainty. -```@example metadata -@variables u [bounds = (-1, 1)] +```@repl metadata +@variables u [bounds = (-1, 1)]; hasbounds(u) -``` - -```@example metadata getbounds(u) ``` @@ -95,51 +128,45 @@ Bounds can also be specified for array variables. A scalar array bound is applie element of the array. A bound may also be specified as an array, in which case the size of the array must match the size of the symbolic variable. -```@example metadata -@variables x[1:2, 1:2] [bounds = (-1, 1)] +```@repl metadata +@variables x[1:2, 1:2] [bounds = (-1, 1)]; hasbounds(x) -``` - -```@example metadata getbounds(x) -``` - -```@example metadata getbounds(x[1, 1]) -``` - -```@example metadata getbounds(x[1:2, 1]) -``` - -```@example metadata -@variables x[1:2] [bounds = (-Inf, [1.0, Inf])] +@variables x[1:2] [bounds = (-Inf, [1.0, Inf])]; hasbounds(x) -``` - -```@example metadata getbounds(x) -``` - -```@example metadata getbounds(x[2]) +hasbounds(x[2]) ``` -```@example metadata -hasbounds(x[2]) +```@docs +hasbounds +getbounds ``` ## Guess -Specify an initial guess for custom initial conditions of an `ODESystem`. +Specify an initial guess for variables of a `System`. This is used when building the +[`InitializationProblem`](@ref). -```@example metadata -@variables u [guess = 1] +```@repl metadata +@variables u [guess = 1]; hasguess(u) +getguess(u) ``` -```@example metadata -getguess(u) +```@docs +hasguess +getguess +``` + +When a system is constructed, the guesses of the involved variables are stored in a `Dict` +in the system. After this point, the guess metadata of the variable is irrelevant. + +```@docs; canonical = false +guesses ``` ## Mark input as a disturbance @@ -151,6 +178,10 @@ Indicate that an input is not available for control, i.e., it's a disturbance in isdisturbance(u) ``` +```@docs +isdisturbance +``` + ## Mark parameter as tunable Indicate that a parameter can be automatically tuned by parameter optimization or automatic control tuning apps. @@ -160,20 +191,32 @@ Indicate that a parameter can be automatically tuned by parameter optimization o istunable(Kp) ``` +```@docs +istunable +ModelingToolkit.isconstant +``` + +!!! note + + [`@constants`](@ref) is a convenient way to create `@parameters` with `tunable = false` + metadata + ## Probability distributions A probability distribution may be associated with a parameter to indicate either uncertainty about its value, or as a prior distribution for Bayesian optimization. -```julia -using Distributions -d = Normal(10, 1) -@parameters m [dist = d] +```@repl metadata +using Distributions; +d = Normal(10, 1); +@parameters m [dist = d]; hasdist(m) +getdist(m) ``` -```julia -getdist(m) +```@docs +hasdist +getdist ``` ## Irreducible @@ -187,6 +230,10 @@ it can be accessed in [callbacks](@ref events) isirreducible(important_value) ``` +```@docs +isirreducible +``` + ## State Priority When a model is structurally simplified, the algorithm will try to ensure that the variables with higher state priority become states of the system. A variable's state priority is a number set using the `state_priority` metadata. @@ -196,18 +243,24 @@ When a model is structurally simplified, the algorithm will try to ensure that t state_priority(important_dof) ``` +```@docs +state_priority +``` + ## Units Units for variables can be designated using symbolic metadata. For more information, please see the [model validation and units](@ref units) section of the docs. Note that `getunit` is not equivalent to `get_unit` - the former is a metadata getter for individual variables (and is provided so the same interface function for `unit` exists like other metadata), while the latter is used to handle more general symbolic expressions. -```@example metadata -using DynamicQuantities -@variables speed [unit = u"m/s"] +```@repl metadata +using DynamicQuantities; +@variables speed [unit = u"m/s"]; hasunit(speed) +getunit(speed) ``` -```@example metadata -getunit(speed) +```@docs +hasunit +getunit ``` ## Miscellaneous metadata @@ -215,13 +268,23 @@ getunit(speed) User-defined metadata can be added using the `misc` metadata. This can be queried using the `hasmisc` and `getmisc` functions. -```@example metadata -@variables u [misc = :conserved_parameter] y [misc = [2, 4, 6]] +```@repl metadata +@variables u [misc = :conserved_parameter] y [misc = [2, 4, 6]]; hasmisc(u) +getmisc(y) ``` -```@example metadata -getmisc(y) +```@docs +hasmisc +getmisc +``` + +## Dumping metadata + +ModelingToolkit allows dumping the metadata of a variable as a `NamedTuple`. + +```@docs +ModelingToolkit.dump_variable_metadata ``` ## Additional functions @@ -250,25 +313,52 @@ lb, ub = getbounds(p) # operating on a vector, we get lower and upper bound vect b = getbounds(sys) # Operating on the system, we get a dict ``` -See also: [`ModelingToolkit.dump_variable_metadata`](@ref), [`ModelingToolkit.dump_parameters`](@ref), -[`ModelingToolkit.dump_unknowns`](@ref). +See also: + +```@docs; canonical = false +tunable_parameters +ModelingToolkit.dump_unknowns +ModelingToolkit.dump_parameters +``` + +## Symbolic operators -## Index +ModelingToolkit makes heavy use of "operators". These are custom functions that are applied +to symbolic variables. The most common operator is the `Differential` operator, defined in +Symbolics.jl. -```@index -Pages = ["Variable_metadata.md"] +```@docs +Symbolics.Differential ``` -## Docstrings +ModelingToolkit also defines a plethora of custom operators. -```@autodocs -Modules = [ModelingToolkit] -Pages = ["variables.jl"] -Private = false +```@docs +Pre +Initial +Shift ``` +While not an operator, `ShiftIndex` is commonly used to use `Shift` operators in a more +convenient way when writing discrete systems. + ```@docs -ModelingToolkit.dump_variable_metadata -ModelingToolkit.dump_parameters -ModelingToolkit.dump_unknowns +ShiftIndex +``` + +### Sampled time operators + +The following operators are used in hybrid ODE systems, where part of the dynamics of the +system happen at discrete intervals on a clock. While ModelingToolkit cannot yet simulate +such systems, it has the capability to represent them. + +!!! warn + + These operators are considered experimental API. + +```@docs +Sample +Hold +SampleTime +sampletime ``` diff --git a/docs/src/basics/AbstractSystem.md b/docs/src/basics/AbstractSystem.md index d1707f822f..62f842b2b2 100644 --- a/docs/src/basics/AbstractSystem.md +++ b/docs/src/basics/AbstractSystem.md @@ -152,7 +152,7 @@ a lower level in the system. ## Namespacing By default, unsimplified systems will namespace variables accessed via `getproperty`. -Systems created via `@mtkbuild`, or ones passed through `structural_simplify` or +Systems created via `@mtkcompile`, or ones passed through `mtkcompile` or `complete` will not perform this namespacing. However, all of these processes modify the system in a variety of ways. To toggle namespacing without transforming any other property of the system, use `toggle_namespacing`. diff --git a/docs/src/basics/Composition.md b/docs/src/basics/Composition.md index d3de71d696..f6abe97d38 100644 --- a/docs/src/basics/Composition.md +++ b/docs/src/basics/Composition.md @@ -42,7 +42,7 @@ equations(connected) # Differential(t)(decay1₊x(t)) ~ decay1₊f(t) - (decay1₊a*(decay1₊x(t))) # Differential(t)(decay2₊x(t)) ~ decay2₊f(t) - (decay2₊a*(decay2₊x(t))) -simplified_sys = structural_simplify(connected) +simplified_sys = mtkcompile(connected) equations(simplified_sys) ``` @@ -84,7 +84,7 @@ example, let's say there is a variable `x` in `unknowns` and a variable `x` in `subsys`. We can declare that these two variables are the same by specifying their equality: `x ~ subsys.x` in the `eqs` for `sys`. This algebraic relationship can then be simplified by transformations -like `structural_simplify` which will be described later. +like `mtkcompile` which will be described later. ### Numerics with Composed Models @@ -169,7 +169,7 @@ parameters(level3) In many cases, the nicest way to build a model may leave a lot of unnecessary variables. Thus one may want to remove these equations -before numerically solving. The `structural_simplify` function removes +before numerically solving. The `mtkcompile` function removes these trivial equality relationships and trivial singularity equations, i.e. equations which result in `0~0` expressions, in over-specified systems. @@ -227,7 +227,7 @@ values. The user of this model can then solve this model simply by specifying the values at the highest level: ```@example compose -sireqn_simple = structural_simplify(sir) +sireqn_simple = mtkcompile(sir) equations(sireqn_simple) ``` @@ -251,7 +251,7 @@ sol[reqn.R] ## Tearing Problem Construction Some system types (specifically `NonlinearSystem`) can be further -reduced if `structural_simplify` has already been applied to them. This is done +reduced if `mtkcompile` has already been applied to them. This is done by using the alternative problem constructors (`BlockNonlinearProblem`). In these cases, the constructor uses the knowledge of the strongly connected components calculated during the process of simplification diff --git a/docs/src/basics/Debugging.md b/docs/src/basics/Debugging.md index 6e2d471461..5bc509acfd 100644 --- a/docs/src/basics/Debugging.md +++ b/docs/src/basics/Debugging.md @@ -13,7 +13,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D eqs = [D(u1) ~ -√(u1), D(u2) ~ -√(u2), D(u3) ~ -√(u3)] defaults = [u1 => 1.0, u2 => 2.0, u3 => 3.0] @named sys = ODESystem(eqs, t; defaults) -sys = structural_simplify(sys) +sys = mtkcompile(sys) ``` This problem causes the ODE solver to crash: @@ -38,7 +38,7 @@ We could have figured that out ourselves, but it is not always so obvious for mo Suppose we also want to validate that `u1 + u2 >= 2.0`. We can do this via the assertions functionality. ```@example debug -@mtkbuild sys = ODESystem(eqs, t; defaults, assertions = [(u1 + u2 >= 2.0) => "Oh no!"]) +@mtkcompile sys = ODESystem(eqs, t; defaults, assertions = [(u1 + u2 >= 2.0) => "Oh no!"]) ``` The assertions must be an iterable of pairs, where the first element is the symbolic condition and diff --git a/docs/src/basics/Events.md b/docs/src/basics/Events.md index 2808204d61..5d5df0a377 100644 --- a/docs/src/basics/Events.md +++ b/docs/src/basics/Events.md @@ -59,7 +59,7 @@ For example, consider the following system. ```julia @variables x(t) y(t) @parameters p(t) -@mtkbuild sys = ODESystem([x * y ~ p, D(x) ~ 0], t) +@mtkcompile sys = ODESystem([x * y ~ p, D(x) ~ 0], t) event = [t == 1] => [x ~ Pre(x) + 1] ``` @@ -134,7 +134,7 @@ function UnitMassWithFriction(k; name) D(v) ~ sin(t) - k * sign(v)] ODESystem(eqs, t; continuous_events = [v ~ 0], name) # when v = 0 there is a discontinuity end -@mtkbuild m = UnitMassWithFriction(0.7) +@mtkcompile m = UnitMassWithFriction(0.7) prob = ODEProblem(m, Pair[], (0, 10pi)) sol = solve(prob, Tsit5()) plot(sol) @@ -154,8 +154,9 @@ like this root_eqs = [x ~ 0] # the event happens at the ground x(t) = 0 affect = [v ~ -Pre(v)] # the effect is that the velocity changes sign -@mtkbuild ball = ODESystem([D(x) ~ v - D(v) ~ -9.8], t; continuous_events = root_eqs => affect) # equation => affect +@mtkcompile ball = ODESystem( + [D(x) ~ v + D(v) ~ -9.8], t; continuous_events = root_eqs => affect) # equation => affect tspan = (0.0, 5.0) prob = ODEProblem(ball, Pair[], tspan) @@ -174,7 +175,7 @@ Multiple events? No problem! This example models a bouncing ball in 2D that is e continuous_events = [[x ~ 0] => [vx ~ -Pre(vx)] [y ~ -1.5, y ~ 1.5] => [vy ~ -Pre(vy)]] -@mtkbuild ball = ODESystem( +@mtkcompile ball = ODESystem( [ D(x) ~ vx, D(y) ~ vy, @@ -254,7 +255,7 @@ end reflect = [x ~ 0] => (bb_affect!, [v], [], [], nothing) -@mtkbuild bb_sys = ODESystem(bb_eqs, t, sts, par, +@mtkcompile bb_sys = ODESystem(bb_eqs, t, sts, par, continuous_events = reflect) u0 = [v => 0.0, x => 1.0] @@ -299,7 +300,7 @@ injection = (t == tinject) => [N ~ Pre(N) + M] u0 = [N => 0.0] tspan = (0.0, 20.0) p = [α => 100.0, tinject => 10.0, M => 50] -@mtkbuild osys = ODESystem(eqs, t, [N], [α, M, tinject]; discrete_events = injection) +@mtkcompile osys = ODESystem(eqs, t, [N], [α, M, tinject]; discrete_events = injection) oprob = ODEProblem(osys, u0, tspan, p) sol = solve(oprob, Tsit5(); tstops = 10.0) plot(sol) @@ -318,7 +319,7 @@ to ```@example events injection = ((t == tinject) & (N < 50)) => [N ~ Pre(N) + M] -@mtkbuild osys = ODESystem(eqs, t, [N], [M, tinject, α]; discrete_events = injection) +@mtkcompile osys = ODESystem(eqs, t, [N], [M, tinject, α]; discrete_events = injection) oprob = ODEProblem(osys, u0, tspan, p) sol = solve(oprob, Tsit5(); tstops = 10.0) plot(sol) @@ -345,7 +346,7 @@ killing = ModelingToolkit.SymbolicDiscreteCallback( tspan = (0.0, 30.0) p = [α => 100.0, tinject => 10.0, M => 50, tkill => 20.0] -@mtkbuild osys = ODESystem(eqs, t, [N], [α, M, tinject, tkill]; +@mtkcompile osys = ODESystem(eqs, t, [N], [α, M, tinject, tkill]; discrete_events = [injection, killing]) oprob = ODEProblem(osys, u0, tspan, p) sol = solve(oprob, Tsit5(); tstops = [10.0, 20.0]) @@ -374,7 +375,7 @@ killing = ModelingToolkit.SymbolicDiscreteCallback( [20.0] => [α ~ 0.0], discrete_parameters = α, iv = t) p = [α => 100.0, M => 50] -@mtkbuild osys = ODESystem(eqs, t, [N], [α, M]; +@mtkcompile osys = ODESystem(eqs, t, [N], [α, M]; discrete_events = [injection, killing]) oprob = ODEProblem(osys, u0, tspan, p) sol = solve(oprob, Tsit5()) @@ -414,7 +415,7 @@ example: ev = ModelingToolkit.SymbolicDiscreteCallback( 1.0 => [c ~ Pre(c) + 1], discrete_parameters = c, iv = t) -@mtkbuild sys = ODESystem( +@mtkcompile sys = ODESystem( D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [ev]) prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0]) @@ -435,7 +436,7 @@ will be saved. If we repeat the above example with `c` not a `discrete_parameter @variables x(t) @parameters c(t) -@mtkbuild sys = ODESystem( +@mtkcompile sys = ODESystem( D(x) ~ c * cos(x), t, [x], [c]; discrete_events = [1.0 => [c ~ Pre(c) + 1]]) prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0]) @@ -538,7 +539,7 @@ to the system. ```@example events @named sys = ODESystem( eqs, t, [temp], params; continuous_events = [furnace_disable, furnace_enable]) -ss = structural_simplify(sys) +ss = mtkcompile(sys) prob = ODEProblem(ss, [temp => 0.0, furnace_on => true], (0.0, 10.0)) sol = solve(prob, Tsit5()) plot(sol) @@ -651,7 +652,7 @@ We can now simulate the encoder. ```@example events @named sys = ODESystem( eqs, t, [theta, omega], params; continuous_events = [qAevt, qBevt]) -ss = structural_simplify(sys) +ss = mtkcompile(sys) prob = ODEProblem(ss, [theta => 0.0], (0.0, pi)) sol = solve(prob, Tsit5(); dtmax = 0.01) sol.ps[cnt] diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index 10671299c6..e83b1f1336 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -28,7 +28,7 @@ are similarly undocumented. Following is the list of behaviors that should be re - `setindex!(::MTKParameters, value, ::ParameterIndex)` can be used to set the value of a parameter with the given index. - `parameter_values(sys, sym)` will return a `ParameterIndex` object if `sys` has been - `complete`d (through `structural_simplify`, `complete` or `@mtkbuild`). + `complete`d (through `mtkcompile`, `complete` or `@mtkcompile`). - `copy(::MTKParameters)` is defined and duplicates the parameter object, including the memory used by the underlying buffers. @@ -194,7 +194,7 @@ p, replace, alias = SciMLStructures.canonicalize(Tunable(), prob.p) # ERROR: ArgumentError: SymbolicUtils.BasicSymbolic{Real}[xˍt(t)] are missing from the variable map. -This error can come up after running `structural_simplify` on a system that generates dummy derivatives (i.e. variables with `ˍt`). For example, here even though all the variables are defined with initial values, the `ODEProblem` generation will throw an error that defaults are missing from the variable map. +This error can come up after running `mtkcompile` on a system that generates dummy derivatives (i.e. variables with `ˍt`). For example, here even though all the variables are defined with initial values, the `ODEProblem` generation will throw an error that defaults are missing from the variable map. ```julia using ModelingToolkit @@ -206,7 +206,7 @@ eqs = [x1 + x2 + 1 ~ 0 x1 + D(x3) + x4 + 3 ~ 0 2 * D(D(x1)) + D(D(x2)) + D(D(x3)) + D(x4) + 4 ~ 0] @named sys = ODESystem(eqs, t) -sys = structural_simplify(sys) +sys = mtkcompile(sys) prob = ODEProblem(sys, [], (0, 1)) ``` @@ -237,7 +237,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D sts = @variables x1(t) = 0.0 eqs = [D(x1) ~ 1.1 * x1] -@mtkbuild sys = ODESystem(eqs, t) +@mtkcompile sys = ODESystem(eqs, t) prob = ODEProblem{false}(sys, [], (0, 1); u0_constructor = x -> SVector(x...)) ``` diff --git a/docs/src/basics/InputOutput.md b/docs/src/basics/InputOutput.md index 4dc5a3d50f..35a3885dbd 100644 --- a/docs/src/basics/InputOutput.md +++ b/docs/src/basics/InputOutput.md @@ -30,7 +30,7 @@ This function takes a vector of variables that are to be considered inputs, i.e. !!! note "Un-simplified system" - This function expects `sys` to be un-simplified, i.e., `structural_simplify` or `@mtkbuild` should not be called on the system before passing it into this function. `generate_control_function` calls a special version of `structural_simplify` internally. + This function expects `sys` to be un-simplified, i.e., `mtkcompile` or `@mtkcompile` should not be called on the system before passing it into this function. `generate_control_function` calls a special version of `mtkcompile` internally. ### Example: diff --git a/docs/src/basics/Linearization.md b/docs/src/basics/Linearization.md index 1c06ce72d4..27b8ec9903 100644 --- a/docs/src/basics/Linearization.md +++ b/docs/src/basics/Linearization.md @@ -29,7 +29,7 @@ eqs = [u ~ kp * (r - y) # P controller D(x) ~ -x + u # First-order plant y ~ x] # Output equation -@named sys = ODESystem(eqs, t) # Do not call @mtkbuild when linearizing +@named sys = ODESystem(eqs, t) # Do not call @mtkcompile when linearizing matrices, simplified_sys = linearize(sys, [r], [y]) # Linearize from r to y matrices ``` @@ -47,7 +47,7 @@ using ModelingToolkit: inputs, outputs !!! note "Un-simplified system" - Linearization expects `sys` to be un-simplified, i.e., `structural_simplify` or `@mtkbuild` should not be called on the system before linearizing. + Linearization expects `sys` to be un-simplified, i.e., `mtkcompile` or `@mtkcompile` should not be called on the system before linearizing. ## Operating point diff --git a/docs/src/basics/MTKLanguage.md b/docs/src/basics/MTKLanguage.md index ba6d2c34b5..05847581f6 100644 --- a/docs/src/basics/MTKLanguage.md +++ b/docs/src/basics/MTKLanguage.md @@ -1,4 +1,4 @@ -# [ModelingToolkit Language: Modeling with `@mtkmodel`, `@connectors` and `@mtkbuild`](@id mtk_language) +# [ModelingToolkit Language: Modeling with `@mtkmodel`, `@connectors` and `@mtkcompile`](@id mtk_language) ## MTK Model @@ -150,7 +150,7 @@ end - Whenever a parameter or variable has initial value, for example `v(t) = 0.0`, a symbolic variable named `v` with initial value 0.0 and a keyword argument `v`, with default value `nothing` are created.
This way, users can optionally pass new value of `v` while creating a component. ```julia -julia> @mtkbuild model_c1 = ModelC(; v = 2.0); +julia> @mtkcompile model_c1 = ModelC(; v = 2.0); julia> ModelingToolkit.getdefault(model_c1.v) 2.0 @@ -182,7 +182,7 @@ One or more partial systems can be extended in a higher system with `@extend` st ```@example mtkmodel-example using ModelingToolkit: getdefault -@mtkbuild model_c3 = ModelC(; model_a.k_array = [1.0, 2.0]) +@mtkcompile model_c3 = ModelC(; model_a.k_array = [1.0, 2.0]) getdefault(model_c3.model_a.k_array[1]) # 1.0 @@ -529,28 +529,28 @@ end ## Build structurally simplified models: -`@mtkbuild` builds an instance of a component and returns a structurally simplied system. +`@mtkcompile` builds an instance of a component and returns a structurally simplied system. ```julia -@mtkbuild sys = CustomModel() +@mtkcompile sys = CustomModel() ``` This is equivalent to: ```julia @named model = CustomModel() -sys = structural_simplify(model) +sys = mtkcompile(model) ``` -Pass keyword arguments to `structural_simplify` using the following syntax: +Pass keyword arguments to `mtkcompile` using the following syntax: ```julia -@mtkbuild sys=CustomModel() fully_determined=false +@mtkcompile sys=CustomModel() fully_determined=false ``` This is equivalent to: ```julia @named model = CustomModel() -sys = structural_simplify(model; fully_determined = false) +sys = mtkcompile(model; fully_determined = false) ``` diff --git a/docs/src/basics/Precompilation.md b/docs/src/basics/Precompilation.md index 97111f0d6b..0bf9a86653 100644 --- a/docs/src/basics/Precompilation.md +++ b/docs/src/basics/Precompilation.md @@ -22,7 +22,7 @@ using ModelingToolkit @variables x(ModelingToolkit.t_nounits) @named sys = ODESystem([ModelingToolkit.D_nounits(x) ~ -x + 1], ModelingToolkit.t_nounits) -prob = ODEProblem(structural_simplify(sys), [x => 30.0], (0, 100), [], +prob = ODEProblem(mtkcompile(sys), [x => 30.0], (0, 100), [], eval_expression = true, eval_module = @__MODULE__) end diff --git a/docs/src/basics/Validation.md b/docs/src/basics/Validation.md index 79c5d0d214..74715d351e 100644 --- a/docs/src/basics/Validation.md +++ b/docs/src/basics/Validation.md @@ -112,7 +112,7 @@ ps = @parameters s=-1 [unit = u"cm"] c=c [unit = u"cm"] eqs = [D(a) ~ dummycomplex(c, s);] sys = ODESystem( eqs, t, [sts...;], [ps...;], name = :sys, checks = ~ModelingToolkit.CheckUnits) -sys_simple = structural_simplify(sys) +sys_simple = mtkcompile(sys) ``` ## `DynamicQuantities` Literals diff --git a/docs/src/comparison.md b/docs/src/comparison.md index 52d5ab2f70..6691b0023a 100644 --- a/docs/src/comparison.md +++ b/docs/src/comparison.md @@ -12,7 +12,7 @@ - All current Modelica compiler implementations are fixed and not extendable by the users from the Modelica language itself. For example, the Dymola compiler [shares its symbolic processing pipeline](https://www.claytex.com/tech-blog/model-translation-and-symbolic-manipulation/), - which is roughly equivalent to the `dae_index_lowering` and `structural_simplify` + which is roughly equivalent to the `dae_index_lowering` and `mtkcompile` of ModelingToolkit.jl. ModelingToolkit.jl is an open and hackable transformation system which allows users to add new non-standard transformations and control the order of application. @@ -90,7 +90,7 @@ [Dymola symbolic processing pipeline](https://www.claytex.com/tech-blog/model-translation-and-symbolic-manipulation/) with some improvements. ModelingToolkit.jl has an open transformation pipeline that allows for users to extend and reorder transformation passes, where - `structural_simplify` is an adaptation of the Modia.jl-improved alias elimination + `mtkcompile` is an adaptation of the Modia.jl-improved alias elimination and tearing algorithms. - Both Modia and ModelingToolkit generate `DAEProblem` and `ODEProblem` forms for solving with [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/). diff --git a/docs/src/examples/higher_order.md b/docs/src/examples/higher_order.md index fac707525f..95480e283b 100644 --- a/docs/src/examples/higher_order.md +++ b/docs/src/examples/higher_order.md @@ -4,7 +4,7 @@ ModelingToolkit has a system for transformations of mathematical systems. These transformations allow for symbolically changing the representation of the model to problems that are easier to numerically solve. One simple to demonstrate transformation, is -`structural_simplify`, which does a lot of tricks, one being the +`mtkcompile`, which does a lot of tricks, one being the transformation that turns an Nth order ODE into N coupled 1st order ODEs. @@ -32,7 +32,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D D(z) ~ x * y - β * z end end -@mtkbuild sys = SECOND_ORDER() +@mtkcompile sys = SECOND_ORDER() ``` The second order ODE has been automatically transformed to two first order ODEs. @@ -43,7 +43,7 @@ and this syntax extends to `N`-th order. Also, we can use `*` or `∘` to compos `Differential`s, like `Differential(t) * Differential(x)`. Now let's transform this into the `ODESystem` of first order components. -We do this by calling `structural_simplify`: +We do this by calling `mtkcompile`: Now we can directly numerically solve the lowered system. Note that, following the original problem, the solution requires knowing the @@ -54,7 +54,7 @@ but we still have to provide a value for the latter. ```@example orderlowering u0 = [D(sys.x) => 2.0] tspan = (0.0, 100.0) -prob = ODEProblem(sys, u0, tspan, [], jac = true) +prob = ODEProblem(sys, u0, tspan, jac = true) sol = solve(prob, Tsit5()) using Plots plot(sol, idxs = (sys.x, sys.y)) diff --git a/docs/src/examples/modelingtoolkitize_index_reduction.md b/docs/src/examples/modelingtoolkitize_index_reduction.md index b19ea46701..9759dc2081 100644 --- a/docs/src/examples/modelingtoolkitize_index_reduction.md +++ b/docs/src/examples/modelingtoolkitize_index_reduction.md @@ -29,7 +29,7 @@ p = [9.8, 1] tspan = (0, 10.0) pendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p) traced_sys = modelingtoolkitize(pendulum_prob) -pendulum_sys = structural_simplify(dae_index_lowering(traced_sys)) +pendulum_sys = mtkcompile(dae_index_lowering(traced_sys)) prob = ODEProblem(pendulum_sys, [], tspan) sol = solve(prob, Rodas5P(), abstol = 1e-8, reltol = 1e-8) plot(sol, idxs = unknowns(traced_sys)) @@ -157,7 +157,7 @@ numerical solver. Let's try that out: ```@example indexred traced_sys = modelingtoolkitize(pendulum_prob) -pendulum_sys = structural_simplify(dae_index_lowering(traced_sys)) +pendulum_sys = mtkcompile(dae_index_lowering(traced_sys)) prob = ODEProblem(pendulum_sys, Pair[], tspan) sol = solve(prob, Rodas5P()) diff --git a/docs/src/examples/perturbation.md b/docs/src/examples/perturbation.md index 0d1a493cb4..20cef4067c 100644 --- a/docs/src/examples/perturbation.md +++ b/docs/src/examples/perturbation.md @@ -43,7 +43,7 @@ eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2) These are the ODEs we want to solve. Now construct an `ODESystem`, which automatically inserts dummy derivatives for the velocities: ```@example perturbation -@mtkbuild sys = ODESystem(eqs_pert, t) +@mtkcompile sys = System(eqs_pert, t) ``` To solve the `ODESystem`, we generate an `ODEProblem` with initial conditions $x(0) = 0$, and $ẋ(0) = 1$, and solve it: @@ -87,7 +87,7 @@ We follow the same steps as in the previous example to construct the `ODESystem` ```@example perturbation eq_pert = substitute(eq, x => x_series) eqs_pert = taylor_coeff(eq_pert, ϵ, 0:2) -@mtkbuild sys = ODESystem(eqs_pert, t) +@mtkcompile sys = System(eqs_pert, t) ``` We solve and plot it as in the previous example, and compare the solution with $ϵ=0.1$ to the exact solution $x(t, ϵ) = e^{-ϵ t} \sin(\sqrt{(1-ϵ^2)}\,t) / \sqrt{1-ϵ^2}$ of the unperturbed equation: diff --git a/docs/src/examples/remake.md b/docs/src/examples/remake.md index 91dba4d7ae..9d40a83dbf 100644 --- a/docs/src/examples/remake.md +++ b/docs/src/examples/remake.md @@ -14,7 +14,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @variables x(t) y(t) eqs = [D(x) ~ (α - β * y) * x D(y) ~ (δ * x - γ) * y] -@mtkbuild odesys = ODESystem(eqs, t) +@mtkcompile odesys = System(eqs, t) ``` To create the "data" for optimization, we will solve the system with a known set of @@ -24,7 +24,7 @@ parameters. using OrdinaryDiffEq odeprob = ODEProblem( - odesys, [x => 1.0, y => 1.0], (0.0, 10.0), [α => 1.5, β => 1.0, γ => 3.0, δ => 1.0]) + odesys, [x => 1.0, y => 1.0, α => 1.5, β => 1.0, γ => 3.0, δ => 1.0], (0.0, 10.0)) timesteps = 0.0:0.1:10.0 sol = solve(odeprob, Tsit5(); saveat = timesteps) data = Array(sol) diff --git a/docs/src/examples/sparse_jacobians.md b/docs/src/examples/sparse_jacobians.md index 03bd80d432..79a93c3471 100644 --- a/docs/src/examples/sparse_jacobians.md +++ b/docs/src/examples/sparse_jacobians.md @@ -54,7 +54,7 @@ prob = ODEProblem(brusselator_2d_loop, u0, (0.0, 11.5), p) Now let's use `modelingtoolkitize` to generate the symbolic version: ```@example sparsejac -@mtkbuild sys = modelingtoolkitize(prob); +@mtkcompile sys = modelingtoolkitize(prob); nothing # hide ``` diff --git a/docs/src/examples/spring_mass.md b/docs/src/examples/spring_mass.md index 355e5c20b2..8d42592796 100644 --- a/docs/src/examples/spring_mass.md +++ b/docs/src/examples/spring_mass.md @@ -13,13 +13,13 @@ function Mass(; name, m = 1.0, xy = [0.0, 0.0], u = [0.0, 0.0]) ps = @parameters m = m sts = @variables pos(t)[1:2]=xy v(t)[1:2]=u eqs = scalarize(D.(pos) .~ v) - ODESystem(eqs, t, [pos..., v...], ps; name) + System(eqs, t, [pos..., v...], ps; name) end function Spring(; name, k = 1e4, l = 1.0) ps = @parameters k=k l=l @variables x(t), dir(t)[1:2] - ODESystem(Equation[], t, [x, dir...], ps; name) + System(Equation[], t, [x, dir...], ps; name) end function connect_spring(spring, a, b) @@ -43,9 +43,9 @@ g = [0.0, -9.81] eqs = [connect_spring(spring, mass.pos, center) scalarize(D.(mass.v) .~ spring_force(spring) / mass.m .+ g)] -@named _model = ODESystem(eqs, t, [spring.x; spring.dir; mass.pos], []) +@named _model = System(eqs, t, [spring.x; spring.dir; mass.pos], []) @named model = compose(_model, mass, spring) -sys = structural_simplify(model) +sys = mtkcompile(model) prob = ODEProblem(sys, [], (0.0, 3.0)) sol = solve(prob, Rosenbrock23()) @@ -56,18 +56,18 @@ plot(sol) ### Building the components -For each component, we use a Julia function that returns an `ODESystem`. At the top, we define the fundamental properties of a `Mass`: it has a mass `m`, a position `pos` and a velocity `v`. We also define that the velocity is the rate of change of position with respect to time. +For each component, we use a Julia function that returns an `System`. At the top, we define the fundamental properties of a `Mass`: it has a mass `m`, a position `pos` and a velocity `v`. We also define that the velocity is the rate of change of position with respect to time. ```@example component function Mass(; name, m = 1.0, xy = [0.0, 0.0], u = [0.0, 0.0]) ps = @parameters m = m sts = @variables pos(t)[1:2]=xy v(t)[1:2]=u eqs = scalarize(D.(pos) .~ v) - ODESystem(eqs, t, [pos..., v...], ps; name) + System(eqs, t, [pos..., v...], ps; name) end ``` -Note that this is an incompletely specified `ODESystem`. It cannot be simulated on its own, since the equations for the velocity `v[1:2](t)` are unknown. Notice the addition of a `name` keyword. This allows us to generate different masses with different names. A `Mass` can now be constructed as: +Note that this is an incompletely specified `System`. It cannot be simulated on its own, since the equations for the velocity `v[1:2](t)` are unknown. Notice the addition of a `name` keyword. This allows us to generate different masses with different names. A `Mass` can now be constructed as: ```@example component Mass(name = :mass1) @@ -85,7 +85,7 @@ Next, we build the spring component. It is characterized by the spring constant function Spring(; name, k = 1e4, l = 1.0) ps = @parameters k=k l=l @variables x(t), dir(t)[1:2] - ODESystem(Equation[], t, [x, dir...], ps; name) + System(Equation[], t, [x, dir...], ps; name) end ``` @@ -129,7 +129,7 @@ eqs = [connect_spring(spring, mass.pos, center) Finally, we can build the model using these equations and components. ```@example component -@named _model = ODESystem(eqs, t, [spring.x; spring.dir; mass.pos], []) +@named _model = System(eqs, t, [spring.x; spring.dir; mass.pos], []) @named model = compose(_model, mass, spring) ``` @@ -153,10 +153,10 @@ parameters(model) ### Simplifying and solving this system -This system can be solved directly as a DAE using [one of the DAE solvers from DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/). However, we can symbolically simplify the system first beforehand. Running `structural_simplify` eliminates unnecessary variables from the model to give the leanest numerical representation of the system. +This system can be solved directly as a DAE using [one of the DAE solvers from DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/). However, we can symbolically simplify the system first beforehand. Running `mtkcompile` eliminates unnecessary variables from the model to give the leanest numerical representation of the system. ```@example component -sys = structural_simplify(model) +sys = mtkcompile(model) equations(sys) ``` @@ -177,7 +177,7 @@ sol = solve(prob, Rosenbrock23()) plot(sol) ``` -What if we want the timeseries of a different variable? That information is not lost! Instead, `structural_simplify` simply changes unknown variables into `observed` variables. +What if we want the timeseries of a different variable? That information is not lost! Instead, `mtkcompile` simply changes unknown variables into `observed` variables. ```@example component observed(sys) diff --git a/docs/src/examples/tearing_parallelism.md b/docs/src/examples/tearing_parallelism.md index 9540e610bd..4ecc8d2a45 100644 --- a/docs/src/examples/tearing_parallelism.md +++ b/docs/src/examples/tearing_parallelism.md @@ -1,7 +1,7 @@ # Exposing More Parallelism By Tearing Algebraic Equations in ODESystems Sometimes it can be very non-trivial to parallelize a system. In this tutorial, -we will demonstrate how to make use of `structural_simplify` to expose more +we will demonstrate how to make use of `mtkcompile` to expose more parallelism in the solution process and parallelize the resulting simulation. ## The Component Library @@ -16,13 +16,13 @@ using ModelingToolkit: t_nounits as t, D_nounits as D # Basic electric components @connector function Pin(; name) @variables v(t)=1.0 i(t)=1.0 [connect = Flow] - ODESystem(Equation[], t, [v, i], [], name = name) + System(Equation[], t, [v, i], [], name = name) end 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 function ConstantVoltage(; name, V = 1.0) @@ -32,12 +32,12 @@ function ConstantVoltage(; name, V = 1.0) @parameters V = V eqs = [V ~ p.v - n.v 0 ~ p.i + n.i] - compose(ODESystem(eqs, t, [], [V], name = name), p, n) + compose(System(eqs, t, [], [V], name = name), p, n) end @connector function HeatPort(; name) @variables T(t)=293.15 Q_flow(t)=0.0 [connect = Flow] - ODESystem(Equation[], t, [T, Q_flow], [], name = name) + System(Equation[], t, [T, Q_flow], [], name = name) end function HeatingResistor(; name, R = 1.0, TAmbient = 293.15, alpha = 1.0) @@ -51,7 +51,7 @@ function HeatingResistor(; name, R = 1.0, TAmbient = 293.15, alpha = 1.0) 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 @@ -62,7 +62,7 @@ function HeatCapacitor(; name, rho = 8050, V = 1, cp = 460, TAmbient = 293.15) 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 @@ -74,7 +74,7 @@ function Capacitor(; name, C = 1.0) eqs = [v ~ p.v - n.v 0 ~ p.i + n.i D(v) ~ p.i / C] - compose(ODESystem(eqs, t, [v], [C], + compose(System(eqs, t, [v], [C], name = name), p, n) end @@ -88,7 +88,7 @@ function parallel_rc_model(i; name, source, ground, R, C) connect(capacitor.n, source.n, ground.g) connect(resistor.h, heat_capacitor.h)] - compose(ODESystem(rc_eqs, t, name = Symbol(name, i)), + compose(System(rc_eqs, t, name = Symbol(name, i)), [resistor, capacitor, source, ground, heat_capacitor]) end ``` @@ -114,7 +114,7 @@ eqs = [ D(E) ~ sum(((i, sys),) -> getproperty(sys, Symbol(:resistor, i)).h.Q_flow, enumerate(rc_systems)) ] -@named _big_rc = ODESystem(eqs, t, [E], []) +@named _big_rc = System(eqs, t, [E], []) @named big_rc = compose(_big_rc, rc_systems) ``` @@ -122,7 +122,7 @@ Now let's say we want to expose a bit more parallelism via running tearing. How do we do that? ```@example tearing -sys = structural_simplify(big_rc) +sys = mtkcompile(big_rc) ``` Done, that's it. There's no more to it. @@ -175,5 +175,5 @@ so this is better than trying to do it by hand. After performing this, you can construct the `ODEProblem` and set `parallel_form` to use the exposed parallelism in multithreaded function -constructions, but this showcases why `structural_simplify` is so important +constructions, but this showcases why `mtkcompile` is so important to that process. diff --git a/docs/src/internals.md b/docs/src/internals.md index 00b29f1a64..ed83192f21 100644 --- a/docs/src/internals.md +++ b/docs/src/internals.md @@ -18,7 +18,7 @@ and are then used to generate the `observed` equation found in the variable when necessary. In this sense, there is an equivalence between observables and the variable elimination system. -The procedure for variable elimination inside [`structural_simplify`](@ref) is +The procedure for variable elimination inside [`mtkcompile`](@ref) is 1. [`ModelingToolkit.initialize_system_structure`](@ref). 2. [`ModelingToolkit.alias_elimination`](@ref). This step moves equations into `observed(sys)`. diff --git a/docs/src/systems/DiscreteSystem.md b/docs/src/systems/DiscreteSystem.md deleted file mode 100644 index 55a02e5714..0000000000 --- a/docs/src/systems/DiscreteSystem.md +++ /dev/null @@ -1,33 +0,0 @@ -# DiscreteSystem - -## System Constructors - -```@docs -DiscreteSystem -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the discrete system. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the discrete system. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the discrete system. - - `get_iv(sys)`: The independent variable of the discrete system - -## Transformations - -```@docs; canonical=false -structural_simplify -``` - -## Problem Constructors - -```@docs; canonical=false -DiscreteProblem(sys::DiscreteSystem, u0map, tspan) -DiscreteFunction(sys::DiscreteSystem, args...) -``` - -## Discrete Domain - -```@docs; canonical=false -Shift -``` diff --git a/docs/src/systems/ImplicitDiscreteSystem.md b/docs/src/systems/ImplicitDiscreteSystem.md deleted file mode 100644 index d687502b49..0000000000 --- a/docs/src/systems/ImplicitDiscreteSystem.md +++ /dev/null @@ -1,33 +0,0 @@ -# ImplicitDiscreteSystem - -## System Constructors - -```@docs -ImplicitDiscreteSystem -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the implicit discrete system. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the implicit discrete system. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the implicit discrete system. - - `get_iv(sys)`: The independent variable of the implicit discrete system - -## Transformations - -```@docs; canonical=false -structural_simplify -``` - -## Problem Constructors - -```@docs; canonical=false -ImplicitDiscreteProblem(sys::ImplicitDiscreteSystem, u0map, tspan) -ImplicitDiscreteFunction(sys::ImplicitDiscreteSystem, args...) -``` - -## Discrete Domain - -```@docs; canonical=false -Shift -``` diff --git a/docs/src/systems/JumpSystem.md b/docs/src/systems/JumpSystem.md deleted file mode 100644 index 5bd0d50602..0000000000 --- a/docs/src/systems/JumpSystem.md +++ /dev/null @@ -1,33 +0,0 @@ -# JumpSystem - -## System Constructors - -```@docs -JumpSystem -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the jump system. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the jump system. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the jump system. - - `get_iv(sys)`: The independent variable of the jump system. - - `discrete_events(sys)`: The set of discrete events in the jump system. - -## Transformations - -```@docs; canonical=false -structural_simplify -``` - -## Analyses - -## Problem Constructors - -```@docs; canonical=false -DiscreteProblem(sys::JumpSystem, u0map, tspan) -``` - -```@docs -JumpProblem(sys::JumpSystem, prob, aggregator) -``` diff --git a/docs/src/systems/NonlinearSystem.md b/docs/src/systems/NonlinearSystem.md deleted file mode 100644 index 06d587b1b9..0000000000 --- a/docs/src/systems/NonlinearSystem.md +++ /dev/null @@ -1,57 +0,0 @@ -# NonlinearSystem - -## System Constructors - -```@docs -NonlinearSystem -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the nonlinear system. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the nonlinear system. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the nonlinear system. - - `get_u0_p(sys, u0map, parammap)` Numeric arrays for the initial condition and parameters given `var => value` maps. - -## Transformations - -```@docs; canonical=false -structural_simplify -alias_elimination -tearing -``` - -## Analyses - -```@docs; canonical=false -ModelingToolkit.isaffine -ModelingToolkit.islinear -``` - -## Applicable Calculation and Generation Functions - -```@docs; canonical=false -calculate_jacobian -generate_jacobian -jacobian_sparsity -``` - -## Problem Constructors - -```@docs -NonlinearFunction(sys::ModelingToolkit.NonlinearSystem, args...) -NonlinearProblem(sys::ModelingToolkit.NonlinearSystem, args...) -``` - -## Torn Problem Constructors - -```@docs -BlockNonlinearProblem -``` - -## Expression Constructors - -```@docs -NonlinearFunctionExpr -NonlinearProblemExpr -``` diff --git a/docs/src/systems/ODESystem.md b/docs/src/systems/ODESystem.md deleted file mode 100644 index 24e2952fc5..0000000000 --- a/docs/src/systems/ODESystem.md +++ /dev/null @@ -1,74 +0,0 @@ -# ODESystem - -## System Constructors - -```@docs -ODESystem -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the ODE. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the ODE. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the ODE. - - `get_iv(sys)`: The independent variable of the ODE. - - `get_u0_p(sys, u0map, parammap)` Numeric arrays for the initial condition and parameters given `var => value` maps. - - `continuous_events(sys)`: The set of continuous events in the ODE. - - `discrete_events(sys)`: The set of discrete events in the ODE. - - `alg_equations(sys)`: The algebraic equations (i.e. that does not contain a differential) that defines the ODE. - - `get_alg_eqs(sys)`: The algebraic equations (i.e. that does not contain a differential) that defines the ODE. Only returns equations of the current-level system. - - `diff_equations(sys)`: The differential equations (i.e. that contain a differential) that defines the ODE. - - `get_diff_eqs(sys)`: The differential equations (i.e. that contain a differential) that defines the ODE. Only returns equations of the current-level system. - - `has_alg_equations(sys)`: Returns `true` if the ODE contains any algebraic equations (i.e. that does not contain a differential). - - `has_alg_eqs(sys)`: Returns `true` if the ODE contains any algebraic equations (i.e. that does not contain a differential). Only considers the current-level system. - - `has_diff_equations(sys)`: Returns `true` if the ODE contains any differential equations (i.e. that does contain a differential). - - `has_diff_eqs(sys)`: Returns `true` if the ODE contains any differential equations (i.e. that does contain a differential). Only considers the current-level system. - -## Transformations - -```@docs -structural_simplify -ode_order_lowering -dae_index_lowering -change_independent_variable -liouville_transform -alias_elimination -tearing -``` - -## Analyses - -```@docs -ModelingToolkit.islinear -ModelingToolkit.isautonomous -ModelingToolkit.isaffine -``` - -## Applicable Calculation and Generation Functions - -```@docs; canonical=false -calculate_jacobian -calculate_tgrad -calculate_factorized_W -generate_jacobian -generate_tgrad -generate_factorized_W -jacobian_sparsity -``` - -## Standard Problem Constructors - -```@docs -ODEFunction(sys::ModelingToolkit.AbstractODESystem, args...) -ODEProblem(sys::ModelingToolkit.AbstractODESystem, args...) -SteadyStateProblem(sys::ModelingToolkit.AbstractODESystem, args...) -DAEProblem(sys::ModelingToolkit.AbstractODESystem, args...) -``` - -## Expression Constructors - -```@docs -ODEFunctionExpr -DAEFunctionExpr -SteadyStateProblemExpr -``` diff --git a/docs/src/systems/OptimizationSystem.md b/docs/src/systems/OptimizationSystem.md deleted file mode 100644 index bcc8b21de7..0000000000 --- a/docs/src/systems/OptimizationSystem.md +++ /dev/null @@ -1,40 +0,0 @@ -# OptimizationSystem - -## System Constructors - -```@docs -OptimizationSystem -``` - -## Composition and Accessor Functions - - - `get_op(sys)`: The objective to be minimized. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns for the optimization. - - `get_ps(sys)` or `parameters(sys)`: The parameters for the optimization. - - `get_constraints(sys)` or `constraints(sys)`: The constraints for the optimization. - -## Transformations - -## Analyses - -## Applicable Calculation and Generation Functions - -```julia -calculate_gradient -calculate_hessian -generate_gradient -generate_hessian -hessian_sparsity -``` - -## Problem Constructors - -```@docs -OptimizationProblem(sys::ModelingToolkit.OptimizationSystem, args...) -``` - -## Expression Constructors - -```@docs -OptimizationProblemExpr -``` diff --git a/docs/src/systems/SDESystem.md b/docs/src/systems/SDESystem.md deleted file mode 100644 index 5789d2d9cb..0000000000 --- a/docs/src/systems/SDESystem.md +++ /dev/null @@ -1,70 +0,0 @@ -# SDESystem - -## System Constructors - -```@docs -SDESystem -``` - -To convert an `ODESystem` to an `SDESystem` directly: - -``` -ode = ODESystem(eqs,t,[x,y,z],[σ,ρ,β]) -sde = SDESystem(ode, noiseeqs) -``` - -## Composition and Accessor Functions - - - `get_eqs(sys)` or `equations(sys)`: The equations that define the SDE. - - `get_unknowns(sys)` or `unknowns(sys)`: The set of unknowns in the SDE. - - `get_ps(sys)` or `parameters(sys)`: The parameters of the SDE. - - `get_iv(sys)`: The independent variable of the SDE. - - `continuous_events(sys)`: The set of continuous events in the SDE. - - `discrete_events(sys)`: The set of discrete events in the SDE. - - `alg_equations(sys)`: The algebraic equations (i.e. that does not contain a differential) that defines the ODE. - - `get_alg_eqs(sys)`: The algebraic equations (i.e. that does not contain a differential) that defines the ODE. Only returns equations of the current-level system. - - `diff_equations(sys)`: The differential equations (i.e. that contain a differential) that defines the ODE. - - `get_diff_eqs(sys)`: The differential equations (i.e. that contain a differential) that defines the ODE. Only returns equations of the current-level system. - - `has_alg_equations(sys)`: Returns `true` if the ODE contains any algebraic equations (i.e. that does not contain a differential). - - `has_alg_eqs(sys)`: Returns `true` if the ODE contains any algebraic equations (i.e. that does not contain a differential). Only considers the current-level system. - - `has_diff_equations(sys)`: Returns `true` if the ODE contains any differential equations (i.e. that does contain a differential). - - `has_diff_eqs(sys)`: Returns `true` if the ODE contains any differential equations (i.e. that does contain a differential). Only considers the current-level system. - -## Transformations - -```@docs; canonical=false -structural_simplify -alias_elimination -``` - -```@docs -ModelingToolkit.Girsanov_transform -``` - -## Analyses - -## Applicable Calculation and Generation Functions - -```@docs; canonical=false -calculate_jacobian -calculate_tgrad -calculate_factorized_W -generate_jacobian -generate_tgrad -generate_factorized_W -jacobian_sparsity -``` - -## Problem Constructors - -```@docs -SDEFunction(sys::ModelingToolkit.SDESystem, args...) -SDEProblem(sys::ModelingToolkit.SDESystem, args...) -``` - -## Expression Constructors - -```@docs -SDEFunctionExpr -SDEProblemExpr -``` diff --git a/docs/src/tutorials/SampledData.md b/docs/src/tutorials/SampledData.md index c700bae5c2..c2d6c9308d 100644 --- a/docs/src/tutorials/SampledData.md +++ b/docs/src/tutorials/SampledData.md @@ -155,10 +155,9 @@ k = ShiftIndex(clock) function plant(; name) @variables x(t)=1 u(t)=0 y(t)=0 - D = Differential(t) eqs = [D(x) ~ -x + u y ~ x] - ODESystem(eqs, t; name = name) + System(eqs, t; name = name) end function filt(; name) # Reference filter @@ -166,7 +165,7 @@ function filt(; name) # Reference filter a = 1 / exp(dt) eqs = [x(k) ~ a * x(k - 1) + (1 - a) * u(k) y ~ x] - ODESystem(eqs, t, name = name) + System(eqs, t, name = name) end function controller(kp; name) @@ -174,7 +173,7 @@ function controller(kp; name) @parameters kp = kp eqs = [yd ~ Sample(y) ud ~ kp * (r - yd)] - ODESystem(eqs, t; name = name) + System(eqs, t; name = name) end @named f = filt() @@ -187,7 +186,7 @@ connections = [r ~ sin(t) # reference signal Hold(c.ud) ~ p.u # controller output to plant input (zero-order-hold) p.y ~ c.y] # plant output to controller feedback -@named cl = ODESystem(connections, t, systems = [f, c, p]) +@named cl = System(connections, t, systems = [f, c, p]) ``` ```@docs diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md index 751b678dae..b364be4012 100644 --- a/docs/src/tutorials/acausal_components.md +++ b/docs/src/tutorials/acausal_components.md @@ -99,7 +99,7 @@ end end end -@mtkbuild rc_model = RCModel(resistor.R = 2.0) +@mtkcompile rc_model = RCModel(resistor.R = 2.0) u0 = [ rc_model.capacitor.v => 0.0 ] @@ -272,7 +272,7 @@ We can create a RCModel component with `@named`. And using `subcomponent_name.pa the parameters or defaults values of variables of subcomponents. ```@example acausal -@mtkbuild rc_model = RCModel(resistor.R = 2.0) +@mtkcompile rc_model = RCModel(resistor.R = 2.0) ``` This model is acausal because we have not specified anything about the causality of the model. We have @@ -323,7 +323,7 @@ plot(sol) By default, this plots only the unknown variables that had to be solved for. However, what if we wanted to plot the timeseries of a different variable? Do not worry, that information was not thrown away! Instead, transformations -like `structural_simplify` simply change unknown variables into observables which are +like `mtkcompile` simply change unknown variables into observables which are defined by `observed` equations. ```@example acausal diff --git a/docs/src/tutorials/attractors.md b/docs/src/tutorials/attractors.md index 317384b01a..426551b017 100644 --- a/docs/src/tutorials/attractors.md +++ b/docs/src/tutorials/attractors.md @@ -41,10 +41,10 @@ eqs = [ Because our dynamical system is super simple, we will directly make an `ODESystem` and cast it in an `ODEProblem` as in the [`ODESystems` tutorial](@ref programmatically). Since all state variables and parameters have a default value we can immediately write ```@example Attractors -@named modlorenz = ODESystem(eqs, t) -ssys = structural_simplify(modlorenz) +@named modlorenz = System(eqs, t) +ssys = mtkcompile(modlorenz) # The timespan given to the problem is irrelevant for DynamicalSystems.jl -prob = ODEProblem(ssys, [], (0.0, 1.0), []) +prob = ODEProblem(ssys, [], (0.0, 1.0)) ``` This `prob` can be turned to a dynamical system as simply as diff --git a/docs/src/tutorials/bifurcation_diagram_computation.md b/docs/src/tutorials/bifurcation_diagram_computation.md index f15d46e1e4..811ea89e83 100644 --- a/docs/src/tutorials/bifurcation_diagram_computation.md +++ b/docs/src/tutorials/bifurcation_diagram_computation.md @@ -14,7 +14,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @parameters μ α eqs = [0 ~ μ * x - x^3 + α * y, 0 ~ -y] -@mtkbuild nsys = NonlinearSystem(eqs, [x, y], [μ, α]) +@mtkcompile nsys = System(eqs, [x, y], [μ, α]) ``` we wish to compute a bifurcation diagram for this system as we vary the parameter `μ`. For this, we need to provide the following information: @@ -95,7 +95,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @parameters μ eqs = [D(x) ~ μ * x - y - x * (x^2 + y^2), D(y) ~ x + μ * y - y * (x^2 + y^2)] -@mtkbuild osys = ODESystem(eqs, t) +@mtkcompile osys = System(eqs, t) bif_par = μ plot_var = x diff --git a/docs/src/tutorials/callable_params.md b/docs/src/tutorials/callable_params.md index 74e7ea87fa..2500c27015 100644 --- a/docs/src/tutorials/callable_params.md +++ b/docs/src/tutorials/callable_params.md @@ -2,7 +2,7 @@ ModelingToolkit.jl allows creating parameters that represent functions to be called. This is especially useful for including interpolants and/or lookup tables inside ODEs. In this -tutorial we will create an `ODESystem` which employs callable parameters to interpolate data +tutorial we will create an `System` which employs callable parameters to interpolate data inside an ODE and go over the various syntax options and their implications. ## Callable parameter syntax @@ -69,7 +69,7 @@ Tspline = typeof(spline) @variables x(t) @parameters (interp::Tspline)(..) -@mtkbuild sys = ODESystem(D(x) ~ interp(t), t) +@mtkcompile sys = System(D(x) ~ interp(t), t) ``` The derivative of `x` is obtained via an interpolation from DataInterpolations.jl. Note @@ -77,7 +77,7 @@ the parameter syntax. The `(..)` marks the parameter as callable. `(interp::Tspl indicates that the parameter is of type `Tspline`. ```@example callable -prob = ODEProblem(sys, [x => 0.0], (0.0, 1.0), [interp => spline]) +prob = ODEProblem(sys, [x => 0.0, interp => spline], (0.0, 1.0)) solve(prob) ``` @@ -85,7 +85,7 @@ Note that the the following will not work: ```julia ODEProblem( - sys; [x => 0.0], (0.0, 1.0), [interp => LinearInterpolation(0.0:0.1:1.0, 0.0:0.1:1.0)]) + sys; [x => 0.0, interp => LinearInterpolation(0.0:0.1:1.0, 0.0:0.1:1.0)], (0.0, 1.0)) ``` Since the type of the spline doesn't match. diff --git a/docs/src/tutorials/change_independent_variable.md b/docs/src/tutorials/change_independent_variable.md index d55639a669..b5b548a73b 100644 --- a/docs/src/tutorials/change_independent_variable.md +++ b/docs/src/tutorials/change_independent_variable.md @@ -24,8 +24,8 @@ using ModelingToolkit: t_nounits as t, D_nounits as D @parameters g=9.81 v # gravitational acceleration and horizontal velocity eqs = [D(D(y)) ~ -g, D(x) ~ v] initialization_eqs = [D(x) ~ D(y)] # 45° initial angle -M1 = ODESystem(eqs, t; initialization_eqs, name = :M) -M1s = structural_simplify(M1) +M1 = System(eqs, t; initialization_eqs, name = :M) +M1s = mtkcompile(M1) @assert length(equations(M1s)) == 3 # hide M1s # hide ``` @@ -44,7 +44,7 @@ This transformation is well-defined for any non-zero horizontal velocity $v$, so ```@example changeivar M2 = change_independent_variable(M1, x) -M2s = structural_simplify(M2; allow_symbolic = true) +M2s = mtkcompile(M2; allow_symbolic = true) # a sanity test on the 10 x/y variables that are accessible to the user # hide @assert allequal([x, M1s.x]) # hide @assert allequal([M2.x, M2s.x]) # hide @@ -67,7 +67,7 @@ It is straightforward to evolve the ODE for 10 meters and plot the resulting tra ```@example changeivar using OrdinaryDiffEq, Plots -prob = ODEProblem(M2s, [M2s.y => 0.0], [0.0, 10.0], [v => 8.0]) # throw 10 meters +prob = ODEProblem(M2s, [M2s.y => 0.0, v => 8.0], [0.0, 10.0]) # throw 10 meters sol = solve(prob, Tsit5()) plot(sol; idxs = M2.y) # must index by M2.y = y(x); not M1.y = y(t)! ``` @@ -88,16 +88,16 @@ In terms of conformal time $t$, they can be written a = GlobalScope(a) # global var needed by all species function species(w; kw...) eqs = [D(Ω) ~ -3(1 + w) * D(a) / a * Ω] - return ODESystem(eqs, t, [Ω], []; kw...) + return System(eqs, t, [Ω], []; kw...) end @named r = species(1 // 3) # radiation @named m = species(0) # matter @named Λ = species(-1) # dark energy / cosmological constant eqs = [Ω ~ r.Ω + m.Ω + Λ.Ω, D(a) ~ √(Ω) * a^2] initialization_eqs = [Λ.Ω + r.Ω + m.Ω ~ 1] -M1 = ODESystem(eqs, t, [Ω, a], []; initialization_eqs, name = :M) +M1 = System(eqs, t, [Ω, a], []; initialization_eqs, name = :M) M1 = compose(M1, r, m, Λ) -M1s = structural_simplify(M1) +M1s = mtkcompile(M1) ``` Of course, we can solve this ODE as it is: @@ -137,7 +137,7 @@ M3 = change_independent_variable(M2, b, [Da(b) ~ exp(-b), a ~ exp(b)]) We can now solve and plot the ODE in terms of $b$: ```@example changeivar -M3s = structural_simplify(M3; allow_symbolic = true) +M3s = mtkcompile(M3; allow_symbolic = true) prob = ODEProblem(M3s, [M3s.r.Ω => 5e-5, M3s.m.Ω => 0.3], (0, -15), []) sol = solve(prob, Tsit5()) @assert Symbol(sol.retcode) == :Success # surrounding text assumes the solution was successful # hide diff --git a/docs/src/tutorials/discrete_system.md b/docs/src/tutorials/discrete_system.md index 8f6828fde7..d070156076 100644 --- a/docs/src/tutorials/discrete_system.md +++ b/docs/src/tutorials/discrete_system.md @@ -1,7 +1,6 @@ # (Experimental) Modeling Discrete Systems -In this example, we will use the new [`DiscreteSystem`](@ref) API -to create an SIR model. +In this example, we will use the [`System`](@ref) API to create an SIR model. ```@example discrete using ModelingToolkit @@ -23,12 +22,12 @@ recovery = rate_to_proportion(γ * h, δt) * I(k - 1) eqs = [S(k) ~ S(k - 1) - infection * h, I(k) ~ I(k - 1) + infection - recovery, R(k) ~ R(k - 1) + recovery] -@mtkbuild sys = DiscreteSystem(eqs, t) +@mtkcompile sys = System(eqs, t) u0 = [S(k - 1) => 990.0, I(k - 1) => 10.0, R(k - 1) => 0.0] p = [β => 0.05, c => 10.0, γ => 0.25, δt => 0.1] tspan = (0.0, 100.0) -prob = DiscreteProblem(sys, u0, tspan, p) +prob = DiscreteProblem(sys, vcat(u0, p), tspan) sol = solve(prob, FunctionMap()) ``` @@ -39,7 +38,7 @@ the Fibonacci series: ```@example discrete @variables x(t) = 1.0 -@mtkbuild sys = DiscreteSystem([x ~ x(k - 1) + x(k - 2)], t) +@mtkcompile sys = System([x ~ x(k - 1) + x(k - 2)], t) ``` The "default value" here should be interpreted as the value of `x` at all past timesteps. diff --git a/docs/src/tutorials/disturbance_modeling.md b/docs/src/tutorials/disturbance_modeling.md index db8d926498..41c7bd86ee 100644 --- a/docs/src/tutorials/disturbance_modeling.md +++ b/docs/src/tutorials/disturbance_modeling.md @@ -67,7 +67,7 @@ This outer model, `ModelWithInputs`, contains two disturbance inputs, both of ty ```@example DISTURBANCE_MODELING @named model = ModelWithInputs() # Model with load disturbance -ssys = structural_simplify(model) +ssys = mtkcompile(model) prob = ODEProblem(ssys, [], (0.0, 6.0)) sol = solve(prob, Tsit5()) using Plots @@ -108,7 +108,7 @@ y &= g(x, p, t) \end{aligned} ``` -where ``x`` is the state, ``y`` are observed variables, ``p`` are parameters, and ``t`` is time. When using MTK, which variables constitute ``x`` and which are considered part of the output, ``y``, is up to the tool rather than the user, this choice is made by MTK during the call to `@mtkbuild` or the lower-level function `structural_simplify`. +where ``x`` is the state, ``y`` are observed variables, ``p`` are parameters, and ``t`` is time. When using MTK, which variables constitute ``x`` and which are considered part of the output, ``y``, is up to the tool rather than the user, this choice is made by MTK during the call to `@mtkcompile` or the lower-level function `mtkcompile`. If we further consider external inputs to the system, such as controlled input signals ``u`` and disturbance inputs ``w``, we can write the system as @@ -169,7 +169,7 @@ end We demonstrate that this model is complete and can be simulated: ```@example DISTURBANCE_MODELING -ssys = structural_simplify(model_with_disturbance) +ssys = mtkcompile(model_with_disturbance) prob = ODEProblem(ssys, [], (0.0, 10.0)) sol = solve(prob, Tsit5()) using Test @@ -191,7 +191,7 @@ g = ModelingToolkit.build_explicit_observed_function( io_sys, outputs; inputs) op = ModelingToolkit.inputs(io_sys) .=> 0 -x0, _ = ModelingToolkit.get_u0_p(io_sys, op, op) +x0 = ModelingToolkit.get_u0(io_sys, op) p = MTKParameters(io_sys, op) u = zeros(1) # Control input w = zeros(length(disturbance_inputs)) # Disturbance input diff --git a/docs/src/tutorials/domain_connections.md b/docs/src/tutorials/domain_connections.md index d6dc2d8781..f9cd3040d2 100644 --- a/docs/src/tutorials/domain_connections.md +++ b/docs/src/tutorials/domain_connections.md @@ -20,7 +20,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D dm(t), [connect = Flow] end - ODESystem(Equation[], t, vars, pars; name, defaults = [dm => 0]) + System(Equation[], t, vars, pars; name, defaults = [dm => 0]) end nothing #hide ``` @@ -47,7 +47,7 @@ The fluid medium setter for `HydralicPort` may be defined as `HydraulicFluid` wi dm ~ 0 ] - ODESystem(eqs, t, vars, pars; name, defaults = [dm => 0]) + System(eqs, t, vars, pars; name, defaults = [dm => 0]) end nothing #hide ``` @@ -63,7 +63,7 @@ Now, we can connect a `HydraulicFluid` component to any `HydraulicPort` connecto eqs = [port.p ~ p] - ODESystem(eqs, t, [], pars; name, systems) + System(eqs, t, [], pars; name, systems) end @component function FixedVolume(; vol, p_int, name) @@ -89,7 +89,7 @@ end rho ~ port.ρ * (1 + p / port.β) dm ~ drho * vol] - ODESystem(eqs, t, vars, pars; name, systems) + System(eqs, t, vars, pars; name, systems) end nothing #hide ``` @@ -97,7 +97,7 @@ nothing #hide When the system is defined we can generate a fluid component and connect it to the system. Here `fluid` is connected to the `src.port`, but it could also be connected to `vol.port`, any connection in the network is fine. ```@example domain -@component function System(; name) +@component function HydraulicSystem(; name) systems = @named begin src = FixedPressure(; p = 200e5) vol = FixedVolume(; vol = 0.1, p_int = 200e5) @@ -108,17 +108,17 @@ When the system is defined we can generate a fluid component and connect it to t eqs = [connect(fluid, src.port) connect(src.port, vol.port)] - ODESystem(eqs, t, [], []; systems, name) + System(eqs, t, [], []; systems, name) end -@named odesys = System() +@named odesys = HydraulicSystem() nothing #hide ``` -To see how the domain works, we can examine the set parameter values for each of the ports `src.port` and `vol.port`. First we assemble the system using `structural_simplify()` and then check the default value of `vol.port.ρ`, whichs points to the setter value `fluid₊ρ`. Likewise, `src.port.ρ`, will also point to the setter value `fluid₊ρ`. Therefore, there is now only 1 defined density value `fluid₊ρ` which sets the density for the connected network. +To see how the domain works, we can examine the set parameter values for each of the ports `src.port` and `vol.port`. First we assemble the system using `mtkcompile()` and then check the default value of `vol.port.ρ`, whichs points to the setter value `fluid₊ρ`. Likewise, `src.port.ρ`, will also point to the setter value `fluid₊ρ`. Therefore, there is now only 1 defined density value `fluid₊ρ` which sets the density for the connected network. ```@repl domain -sys = structural_simplify(odesys) +sys = mtkcompile(odesys) ModelingToolkit.defaults(sys)[odesys.vol.port.ρ] ``` @@ -151,7 +151,7 @@ If we have a more complicated system, for example a hydraulic actuator, with a s port_a.dm ~ +(port_a.ρ) * dx * area port_b.dm ~ -(port_b.ρ) * dx * area] - ODESystem(eqs, t, vars, pars; name, systems) + System(eqs, t, vars, pars; name, systems) end nothing #hide ``` @@ -174,14 +174,14 @@ A system with 2 different fluids is defined and connected to each separate domai connect(src_a.port, act.port_a) connect(src_b.port, act.port_b)] - ODESystem(eqs, t, [], []; systems, name) + System(eqs, t, [], []; systems, name) end @named actsys2 = ActuatorSystem2() nothing #hide ``` -After running `structural_simplify()` on `actsys2`, the defaults will show that `act.port_a.ρ` points to `fluid_a₊ρ` and `act.port_b.ρ` points to `fluid_b₊ρ`. This is a special case, in most cases a hydraulic system will have only 1 fluid, however this simple system has 2 separate domain networks. Therefore, we can connect a single fluid to both networks. This does not interfere with the mathematical equations of the system, since no unknown variables are connected. +After running `mtkcompile()` on `actsys2`, the defaults will show that `act.port_a.ρ` points to `fluid_a₊ρ` and `act.port_b.ρ` points to `fluid_b₊ρ`. This is a special case, in most cases a hydraulic system will have only 1 fluid, however this simple system has 2 separate domain networks. Therefore, we can connect a single fluid to both networks. This does not interfere with the mathematical equations of the system, since no unknown variables are connected. ```@example domain @component function ActuatorSystem1(; name) @@ -198,7 +198,7 @@ After running `structural_simplify()` on `actsys2`, the defaults will show that connect(src_a.port, act.port_a) connect(src_b.port, act.port_b)] - ODESystem(eqs, t, [], []; systems, name) + System(eqs, t, [], []; systems, name) end @named actsys1 = ActuatorSystem1() @@ -224,7 +224,7 @@ In some cases a component will be defined with 2 connectors of the same domain, eqs = [port_a.dm ~ (port_a.p - port_b.p) * K 0 ~ port_a.dm + port_b.dm] - ODESystem(eqs, t, [], pars; systems, name) + System(eqs, t, [], pars; systems, name) end nothing #hide ``` @@ -245,14 +245,14 @@ Adding the `Restrictor` to the original system example will cause a break in the connect(src.port, res.port_a) connect(res.port_b, vol.port)] - ODESystem(eqs, t, [], []; systems, name) + System(eqs, t, [], []; systems, name) end -@mtkbuild ressys = RestrictorSystem() +@mtkcompile ressys = RestrictorSystem() nothing #hide ``` -When `structural_simplify()` is applied to this system it can be seen that the defaults are missing for `res.port_b` and `vol.port`. +When `mtkcompile()` is applied to this system it can be seen that the defaults are missing for `res.port_b` and `vol.port`. ```@repl domain ModelingToolkit.defaults(ressys)[ressys.res.port_a.ρ] @@ -278,10 +278,10 @@ To ensure that the `Restrictor` component does not disrupt the domain network, t port_a.dm ~ (port_a.p - port_b.p) * K 0 ~ port_a.dm + port_b.dm] - ODESystem(eqs, t, [], pars; systems, name) + System(eqs, t, [], pars; systems, name) end -@mtkbuild ressys = RestrictorSystem() +@mtkcompile ressys = RestrictorSystem() nothing #hide ``` diff --git a/docs/src/tutorials/fmi.md b/docs/src/tutorials/fmi.md index 0e01393652..d17452f612 100644 --- a/docs/src/tutorials/fmi.md +++ b/docs/src/tutorials/fmi.md @@ -76,7 +76,7 @@ initialization semantics. We can simulate this model like any other ModelingToolkit system. ```@repl fmi -sys = structural_simplify(model) +sys = mtkcompile(model) prob = ODEProblem(sys, [sys.mass__s => 0.5, sys.mass__v => 0.0], (0.0, 5.0)) sol = solve(prob, Tsit5()) ``` @@ -105,11 +105,11 @@ constant until the next time the callback triggers. The periodic interval must b `communication_step_size` keyword argument. A smaller step size typically leads to less error but is more computationally expensive. -This model alone does not have any differential variables, and calling `structural_simplify` will lead +This model alone does not have any differential variables, and calling `mtkcompile` will lead to an `ODESystem` with no unknowns. ```@example fmi -structural_simplify(inner) +mtkcompile(inner) ``` Simulating this model will cause the OrdinaryDiffEq integrator to immediately finish, and will not @@ -117,7 +117,7 @@ trigger the callback. Thus, we wrap this system in a trivial system with a diffe ```@example fmi @variables x(t) = 1.0 -@mtkbuild sys = ODESystem([D(x) ~ x], t; systems = [inner]) +@mtkcompile sys = System([D(x) ~ x], t; systems = [inner]) ``` We can now simulate `sys`. @@ -184,7 +184,7 @@ metadata. We can now use this component as a subcomponent of a larger system. ```@repl fmi @variables a(t) b(t) c(t) [guess = 1.0]; -@mtkbuild sys = ODESystem( +@mtkcompile sys = System( [adder.a ~ a, adder.b ~ b, D(a) ~ t, D(b) ~ adder.out + adder.c, c^2 ~ adder.out + adder.value], t; @@ -198,8 +198,9 @@ FMUs in initialization to solve for initial states. As mentioned earlier, we can through an FMU. Thus, automatic differentiation has to be disabled for the solver. ```@example fmi -prob = ODEProblem(sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0], - (0.0, 1.0), [sys.adder.value => 2.0]) +prob = ODEProblem( + sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0, sys.adder.value => 2.0], + (0.0, 1.0)) solve(prob, Rodas5P(autodiff = false)) ``` @@ -217,12 +218,13 @@ fmu = loadFMU( @named adder = ModelingToolkit.FMIComponent( Val(2); fmu, type = :CS, communication_step_size = 1e-3, reinitializealg = BrownFullBasicInit()) -@mtkbuild sys = ODESystem( +@mtkcompile sys = System( [adder.a ~ a, adder.b ~ b, D(a) ~ t, D(b) ~ adder.out + adder.c, c^2 ~ adder.out + adder.value], t; systems = [adder]) -prob = ODEProblem(sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0], - (0.0, 1.0), [sys.adder.value => 2.0]) +prob = ODEProblem( + sys, [sys.adder.c => 2.0, sys.a => 1.0, sys.b => 1.0, sys.adder.value => 2.0], + (0.0, 1.0)) solve(prob, Rodas5P()) ``` diff --git a/docs/src/tutorials/initialization.md b/docs/src/tutorials/initialization.md index ba733e0bfb..cdf6e7e813 100644 --- a/docs/src/tutorials/initialization.md +++ b/docs/src/tutorials/initialization.md @@ -44,7 +44,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D eqs = [D(D(x)) ~ λ * x D(D(y)) ~ λ * y - g x^2 + y^2 ~ 1] -@mtkbuild pend = ODESystem(eqs, t) +@mtkcompile pend = System(eqs, t) ``` While we defined the system using second derivatives and a length constraint, @@ -59,7 +59,7 @@ To see how the system works, let's start the pendulum in the far right position, i.e. `x(0) = 1` and `y(0) = 0`. We can do this by: ```@example init -prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], guesses = [λ => 1]) +prob = ODEProblem(pend, [x => 1, y => 0, g => 1], (0.0, 1.5), guesses = [λ => 1]) ``` This solves via: @@ -93,7 +93,7 @@ sol[λ][1] We can similarly choose `λ = 0` and solve for `y` to start the system: ```@example init -prob = ODEProblem(pend, [x => 1, λ => 0], (0.0, 1.5), [g => 1], guesses = [y => 1]) +prob = ODEProblem(pend, [x => 1, λ => 0, g => 1], (0.0, 1.5); guesses = [y => 1]) sol = solve(prob, Rodas5P()) plot(sol, idxs = (x, y)) ``` @@ -102,7 +102,7 @@ or choose to satisfy derivative conditions: ```@example init prob = ODEProblem( - pend, [x => 1, D(y) => 0], (0.0, 1.5), [g => 1], guesses = [λ => 0, y => 1]) + pend, [x => 1, D(y) => 0, g => 1], (0.0, 1.5); guesses = [λ => 0, y => 1]) sol = solve(prob, Rodas5P()) plot(sol, idxs = (x, y)) ``` @@ -114,7 +114,7 @@ We can also directly give equations to be satisfied at the initial point by usin the `initialization_eqs` keyword argument, for example: ```@example init -prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], guesses = [λ => 0, y => 1], +prob = ODEProblem(pend, [x => 1, g => 1], (0.0, 1.5); guesses = [λ => 0, y => 1], initialization_eqs = [y ~ 0]) sol = solve(prob, Rodas5P()) plot(sol, idxs = (x, y)) @@ -125,7 +125,7 @@ variables and parameters: ```@example init prob = ODEProblem( - pend, [x => 1, D(y) => g], (0.0, 3.0), [g => 1], guesses = [λ => 0, y => 1]) + pend, [x => 1, D(y) => g, g => 1], (0.0, 3.0); guesses = [λ => 0, y => 1]) sol = solve(prob, Rodas5P()) plot(sol, idxs = (x, y)) ``` @@ -141,7 +141,7 @@ conditions = getfield.(equations(pend)[3:end], :rhs) when we initialize with ```@example init -prob = ODEProblem(pend, [x => 1, y => 0], (0.0, 1.5), [g => 1], guesses = [y => 0, λ => 1]) +prob = ODEProblem(pend, [x => 1, y => 0, g => 1], (0.0, 1.5); guesses = [y => 0, λ => 1]) ``` we have two extra conditions to satisfy, `x ~ 1` and `y ~ 0` at the initial point. That gives @@ -149,7 +149,7 @@ we have two extra conditions to satisfy, `x ~ 1` and `y ~ 0` at the initial poin case? ```@example init -prob = ODEProblem(pend, [x => 1], (0.0, 1.5), [g => 1], guesses = [y => 0, λ => 1]) +prob = ODEProblem(pend, [x => 1, g => 1], (0.0, 1.5); guesses = [y => 0, λ => 1]) ``` Here we have 4 equations for 5 unknowns (note: the warning is post-simplification of the @@ -167,7 +167,7 @@ direction, we may have an overdetermined system: ```@example init prob = ODEProblem( - pend, [x => 1, y => 0.0, D(y) => 0], (0.0, 1.5), [g => 1], guesses = [λ => 1]) + pend, [x => 1, y => 0.0, D(y) => 0, g => 1], (0.0, 1.5); guesses = [λ => 1]) ``` Can that be solved? @@ -184,7 +184,7 @@ a `SciMLBase.ReturnCode.InitialFailure`: ```@example init prob = ODEProblem( - pend, [x => 1, y => 0.0, D(y) => 2.0, λ => 1], (0.0, 1.5), [g => 1], guesses = [λ => 1]) + pend, [x => 1, y => 0.0, D(y) => 2.0, λ => 1, g => 1], (0.0, 1.5); guesses = [λ => 1]) sol = solve(prob, Rodas5P()) ``` @@ -217,7 +217,7 @@ y-velocities from a given position. ```@example init prob = ODEProblem( - pend, [x => 1, D(y) => 0], (0.0, 1.5), [g => 1], guesses = [λ => 0, y => 1, x => 1]) + pend, [x => 1, D(y) => 0, g => 1], (0.0, 1.5); guesses = [λ => 0, y => 1, x => 1]) sol1 = solve(prob, Rodas5P()) ``` @@ -312,7 +312,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D # hidden @variables x(t) y(t) @parameters total -@mtkbuild sys = ODESystem([D(x) ~ -x, total ~ x + y], t; +@mtkcompile sys = System([D(x) ~ -x, total ~ x + y], t; defaults = [total => missing], guesses = [total => 1.0]) ``` @@ -380,7 +380,7 @@ with observables, those observables are too treated as initial equations. We can resulting simplified system via the command: ```@example init -isys = structural_simplify(isys; fully_determined = false) +isys = mtkcompile(isys; fully_determined = false) ``` Note `fully_determined=false` allows for the simplification to occur when the number of equations @@ -392,7 +392,7 @@ isys = ModelingToolkit.generate_initializesystem( ``` ```@example init -isys = structural_simplify(isys; fully_determined = false) +isys = mtkcompile(isys; fully_determined = false) ``` ```@example init @@ -503,8 +503,8 @@ eqs = [D(x) ~ α * x - β * x * y D(y) ~ -γ * y + δ * x * y z ~ x + y] -@named sys = ODESystem(eqs, t) -simpsys = structural_simplify(sys) +@named sys = System(eqs, t) +simpsys = mtkcompile(sys) tspan = (0.0, 10.0) ``` diff --git a/docs/src/tutorials/linear_analysis.md b/docs/src/tutorials/linear_analysis.md index 7fcb67f9fb..250dbaa6a7 100644 --- a/docs/src/tutorials/linear_analysis.md +++ b/docs/src/tutorials/linear_analysis.md @@ -58,12 +58,14 @@ The following example builds a simple closed-loop system with a plant $P$ and a ```@example LINEAR_ANALYSIS using ModelingToolkitStandardLibrary.Blocks, ModelingToolkit +using ModelingToolkit: t_nounits as t + @named P = FirstOrder(k = 1, T = 1) # A first-order system with pole in -1 @named C = Gain(-1) # A P controller -t = ModelingToolkit.get_iv(P) + eqs = [connect(P.output, :plant_output, C.input) # Connect with an automatically created analysis point called :plant_output connect(C.output, :plant_input, P.input)] -sys = ODESystem(eqs, t, systems = [P, C], name = :feedback_system) +sys = System(eqs, t, systems = [P, C], name = :feedback_system) matrices_S = get_sensitivity(sys, :plant_input)[1] # Compute the matrices of a state-space representation of the (input)sensitivity function. matrices_T = get_comp_sensitivity(sys, :plant_input)[1] diff --git a/docs/src/tutorials/modelingtoolkitize.md b/docs/src/tutorials/modelingtoolkitize.md index 545879f842..fc74578ea4 100644 --- a/docs/src/tutorials/modelingtoolkitize.md +++ b/docs/src/tutorials/modelingtoolkitize.md @@ -52,7 +52,7 @@ If we want to get a symbolic representation, we can simply call `modelingtoolkit on the `prob`, which will return an `ODESystem`: ```@example mtkize -@mtkbuild sys = modelingtoolkitize(prob) +@mtkcompile sys = modelingtoolkitize(prob) ``` Using this, we can symbolically build the Jacobian and then rebuild the ODEProblem: diff --git a/docs/src/tutorials/nonlinear.md b/docs/src/tutorials/nonlinear.md index 057e856229..13caf96231 100644 --- a/docs/src/tutorials/nonlinear.md +++ b/docs/src/tutorials/nonlinear.md @@ -23,12 +23,12 @@ using ModelingToolkit, NonlinearSolve eqs = [0 ~ σ * (y - x) 0 ~ x * (ρ - z) - y 0 ~ x * y - β * z] -@mtkbuild ns = NonlinearSystem(eqs) +@mtkcompile ns = System(eqs) guesses = [x => 1.0, y => 0.0, z => 0.0] ps = [σ => 10.0, ρ => 26.0, β => 8 / 3] -prob = NonlinearProblem(ns, guesses, ps) +prob = NonlinearProblem(ns, vcat(guesses, ps)) sol = solve(prob, NewtonRaphson()) ``` @@ -38,6 +38,6 @@ Just like with `ODEProblem`s we can generate the `NonlinearProblem` with its ana Jacobian function: ```@example nonlinear -prob = NonlinearProblem(ns, guesses, ps, jac = true) +prob = NonlinearProblem(ns, vcat(guesses, ps), jac = true) sol = solve(prob, NewtonRaphson()) ``` diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md index 0a4cd80803..3113ab0fc2 100644 --- a/docs/src/tutorials/ode_modeling.md +++ b/docs/src/tutorials/ode_modeling.md @@ -35,7 +35,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D end using OrdinaryDiffEq -@mtkbuild fol = FOL() +@mtkcompile fol = FOL() prob = ODEProblem(fol, [], (0.0, 10.0), []) sol = solve(prob) @@ -77,12 +77,12 @@ using ModelingToolkit: t_nounits as t, D_nounits as D end end -@mtkbuild fol = FOL() +@mtkcompile fol = FOL() ``` Note that equations in MTK use the tilde character (`~`) as equality sign. -`@mtkbuild` creates an instance of `FOL` named as `fol`. +`@mtkcompile` creates an instance of `FOL` named as `fol`. After construction of the ODE, you can solve it using [OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs/stable/): @@ -90,7 +90,7 @@ After construction of the ODE, you can solve it using [OrdinaryDiffEq.jl](https: using OrdinaryDiffEq using Plots -prob = ODEProblem(fol, [], (0.0, 10.0), []) +prob = ODEProblem(fol, [], (0.0, 10.0)) plot(solve(prob)) ``` @@ -105,15 +105,15 @@ you likely do not want to write an entirely new `@mtkmodel`. ModelingToolkit supports overwriting the default values: ```@example ode2 -@mtkbuild fol_different_values = FOL(; τ = 1 / 3, x = 0.5) -prob = ODEProblem(fol_different_values, [], (0.0, 10.0), []) +@mtkcompile fol_different_values = FOL(; τ = 1 / 3, x = 0.5) +prob = ODEProblem(fol_different_values, [], (0.0, 10.0)) plot(solve(prob)) ``` Alternatively, this overwriting could also have occurred at the `ODEProblem` level. ```@example ode2 -prob = ODEProblem(fol, [fol.τ => 1 / 3], (0.0, 10.0), [fol.x => 0.5]) +prob = ODEProblem(fol, [fol.x => 0.5, fol.τ => 1 / 3], (0.0, 10.0)) plot(solve(prob)) ``` @@ -147,7 +147,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D end end -@mtkbuild fol = FOL() +@mtkcompile fol = FOL() ``` If you copy this block of code to your REPL, you will not see the above LaTeX equations. @@ -175,7 +175,7 @@ in a simulation result. The intermediate variable `RHS` therefore can be plotted along with the unknown variable. Note that this has to be requested explicitly: ```@example ode2 -prob = ODEProblem(fol, [], (0.0, 10.0), []) +prob = ODEProblem(fol, [], (0.0, 10.0)) sol = solve(prob) plot(sol, idxs = [fol.x, fol.RHS]) ``` @@ -221,7 +221,7 @@ Obviously, one could use an explicit, symbolic function of time: end end -@mtkbuild fol_variable_f = FOL() +@mtkcompile fol_variable_f = FOL() ``` However, this function might not be available in an explicit form. @@ -252,11 +252,11 @@ f_fun(t) = t >= 10 ? value_vector[end] : value_vector[Int(floor(t)) + 1] end end -@mtkbuild fol_external_f = FOLExternalFunction() +@mtkcompile fol_external_f = FOLExternalFunction() ``` ```@example ode2 -prob = ODEProblem(fol_external_f, [], (0.0, 10.0), []) +prob = ODEProblem(fol_external_f, [], (0.0, 10.0)) sol = solve(prob) plot(sol, idxs = [fol_external_f.x, fol_external_f.f]) ``` @@ -292,7 +292,7 @@ end fol_2.f ~ fol_1.x end end -@mtkbuild connected = FOLConnected() +@mtkcompile connected = FOLConnected() ``` Here the total model consists of two of the same submodels (components), @@ -316,7 +316,7 @@ initial unknown and the parameter values can be specified accordingly when building the `ODEProblem`: ```@example ode2 -prob = ODEProblem(connected, [], (0.0, 10.0), []) +prob = ODEProblem(connected, [], (0.0, 10.0)) plot(solve(prob)) ``` @@ -344,13 +344,13 @@ Now have MTK provide sparse, analytical derivatives to the solver. This has to be specified during the construction of the `ODEProblem`: ```@example ode2 -prob_an = ODEProblem(connected, [], (0.0, 10.0), []; jac = true) +prob_an = ODEProblem(connected, [], (0.0, 10.0); jac = true) @btime solve(prob_an, Rodas4()); nothing # hide ``` ```@example ode2 -prob_sparse = ODEProblem(connected, [], (0.0, 10.0), []; jac = true, sparse = true) +prob_sparse = ODEProblem(connected, [], (0.0, 10.0); jac = true, sparse = true) @btime solve(prob_sparse, Rodas4()); nothing # hide ``` diff --git a/docs/src/tutorials/optimization.md b/docs/src/tutorials/optimization.md index 20f1079dcd..9eb72b36ea 100644 --- a/docs/src/tutorials/optimization.md +++ b/docs/src/tutorials/optimization.md @@ -24,7 +24,7 @@ using ModelingToolkit, Optimization, OptimizationOptimJL end @parameters a=1.0 b=1.0 rosenbrock = (a - x)^2 + b * (y - x^2)^2 -@mtkbuild sys = OptimizationSystem(rosenbrock, [x, y], [a, b]) +@mtkcompile sys = OptimizationSystem(rosenbrock, [x, y], [a, b]) ``` Every optimization problem consists of a set of optimization variables. @@ -52,7 +52,7 @@ ModelingToolkit is also capable of constructing analytical gradients and Hessian u0 = [y => 2.0] p = [b => 100.0] -prob = OptimizationProblem(sys, u0, p, grad = true, hess = true) +prob = OptimizationProblem(sys, vcat(u0, p), grad = true, hess = true) u_opt = solve(prob, GradientDescent()) ``` @@ -86,7 +86,7 @@ rosenbrock = (a - x)^2 + b * (y - x^2)^2 cons = [ x^2 + y^2 ≲ 1 ] -@mtkbuild sys = OptimizationSystem(rosenbrock, [x, y], [a, b], constraints = cons) +@mtkcompile sys = OptimizationSystem(rosenbrock, [x, y], [a, b], constraints = cons) prob = OptimizationProblem(sys, [], grad = true, hess = true, cons_j = true, cons_h = true) u_opt = solve(prob, IPNewton()) ``` diff --git a/docs/src/tutorials/parameter_identifiability.md b/docs/src/tutorials/parameter_identifiability.md index 0875a7698c..1731d5795f 100644 --- a/docs/src/tutorials/parameter_identifiability.md +++ b/docs/src/tutorials/parameter_identifiability.md @@ -4,7 +4,7 @@ Ordinary differential equations are commonly used for modeling real-world proces We will start by illustrating **local identifiability** in which a parameter is known up to _finitely many values_, and then proceed to determining **global identifiability**, that is, which parameters can be identified _uniquely_. -The package has a standalone data structure for ordinary differential equations, but is also compatible with `ODESystem` type from `ModelingToolkit.jl`. +The package has a standalone data structure for ordinary differential equations, but is also compatible with `System` type from `ModelingToolkit.jl`. ## Local Identifiability @@ -22,7 +22,7 @@ y_2 = x_5\end{cases}$$ This model describes the biohydrogenation[^1] process[^2] with unknown initial conditions. -### Using the `ODESystem` object +### Using the `System` object To define the ode system in Julia, we use `ModelingToolkit.jl`. @@ -61,7 +61,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D end # define the system -@mtkbuild de = Biohydrogenation() +@mtkcompile de = Biohydrogenation() ``` After that, we are ready to check the system for local identifiability: @@ -179,7 +179,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D end end -@mtkbuild ode = GoodwinOscillator() +@mtkcompile ode = GoodwinOscillator() # check only 2 parameters to_check = [ode.b, ode.c] diff --git a/docs/src/tutorials/programmatically_generating.md b/docs/src/tutorials/programmatically_generating.md index 9fc1db1834..406f65d8d9 100644 --- a/docs/src/tutorials/programmatically_generating.md +++ b/docs/src/tutorials/programmatically_generating.md @@ -18,7 +18,8 @@ as demonstrated in the Symbolics.jl documentation. This looks like: ```@example scripting using ModelingToolkit # reexports Symbolics -@variables t x(t) y(t) # Define variables +@independent_variables t +@variables x(t) y(t) # Define variables D = Differential(t) eqs = [D(x) ~ y D(y) ~ x] # Define an array of equations @@ -43,11 +44,11 @@ using ModelingToolkit: t_nounits as t, D_nounits as D eqs = [D(x) ~ (h - x) / τ] # create an array of equations # your first ODE, consisting of a single equation, indicated by ~ -@named model = ODESystem(eqs, t) +@named model = System(eqs, t) # Perform the standard transformations and mark the model complete # Note: Complete models cannot be subsystems of other models! -fol = structural_simplify(model) +fol = mtkcompile(model) prob = ODEProblem(fol, [], (0.0, 10.0), []) using OrdinaryDiffEq sol = solve(prob) @@ -56,22 +57,22 @@ using Plots plot(sol) ``` -As you can see, generating an ODESystem is as simple as creating an array of equations -and passing it to the `ODESystem` constructor. +As you can see, generating an `System` is as simple as creating an array of equations +and passing it to the `System` constructor. -`@named` automatically gives a name to the `ODESystem`, and is shorthand for +`@named` automatically gives a name to the `System`, and is shorthand for ```@example scripting -fol_model = ODESystem(eqs, t; name = :fol_model) # @named fol_model = ODESystem(eqs, t) +fol_model = System(eqs, t; name = :fol_model) # @named fol_model = ODESystem(eqs, t) ``` -Thus, if we had read a name from a file and wish to populate an `ODESystem` with said name, we could do: +Thus, if we had read a name from a file and wish to populate an `System` with said name, we could do: ```@example scripting namesym = :name_from_file -fol_model = ODESystem(eqs, t; name = namesym) +fol_model = System(eqs, t; name = namesym) ``` ## Warning About Mutation -Be advsied that it's never a good idea to mutate an `ODESystem`, or any `AbstractSystem`. +Be advsied that it's never a good idea to mutate an `System`, or any `AbstractSystem`. diff --git a/docs/src/tutorials/stochastic_diffeq.md b/docs/src/tutorials/stochastic_diffeq.md index 79c71a2e8d..9bc0086d35 100644 --- a/docs/src/tutorials/stochastic_diffeq.md +++ b/docs/src/tutorials/stochastic_diffeq.md @@ -53,7 +53,7 @@ eqs = [D(x) ~ σ * (y - x) + 0.3x * B, D(y) ~ x * (ρ - z) - y + 0.3y * B, D(z) ~ x * y - β * z + 0.3z * B] -@mtkbuild de = System(eqs, t) +@mtkcompile de = System(eqs, t) ``` Even though we did not explicitly use `SDESystem`, ModelingToolkit can still infer this from the equations. @@ -65,7 +65,7 @@ typeof(de) We continue by solving and plotting the SDE. ```@example SDE -prob = SDEProblem(de, [], (0.0, 100.0), []) +prob = SDEProblem(de, [], (0.0, 100.0)) sol = solve(prob, SRIW1()) plot(sol, idxs = [(1, 2, 3)]) ``` @@ -87,8 +87,8 @@ multiple `@brownian` variables have to be declared. eqs = [D(x) ~ σ * (y - x) + 0.3x * Bx, D(y) ~ x * (ρ - z) - y + 0.3y * By, D(z) ~ x * y - β * z + 0.3z * Bz] -@mtkbuild de = System(eqs, t) -prob = SDEProblem(de, [], (0.0, 100.0), []) +@mtkcompile de = System(eqs, t) +prob = SDEProblem(de, [], (0.0, 100.0)) sol = solve(prob, SRIW1()) plot(sol) ``` diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 63468917d1..bac555f67b 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -72,10 +72,9 @@ using RuntimeGeneratedFunctions: drop_expr using Symbolics: degree using Symbolics: _parse_vars, value, @derivatives, get_variables, exprs_occur_in, symbolic_linear_solve, build_expr, unwrap, wrap, - VariableSource, getname, variable, Connection, connect, + VariableSource, getname, variable, NAMESPACE_SEPARATOR, set_scalar_metadata, setdefaultval, - initial_state, transition, activeState, entry, hasnode, - ticksInState, timeInState, fixpoint_sub, fast_substitute, + hasnode, fixpoint_sub, fast_substitute, CallWithMetadata, CallWithParent const NAMESPACE_SEPARATOR_SYMBOL = Symbol(NAMESPACE_SEPARATOR) import Symbolics: rename, get_variables!, _solve, hessian_sparsity, @@ -99,6 +98,7 @@ const DQ = DynamicQuantities import DifferentiationInterface as DI using ADTypes: AutoForwardDiff +import SciMLPublic: @public export @derivatives @@ -122,6 +122,16 @@ for fun in [:toexpr] end end +const INTERNAL_FIELD_WARNING = """ +This field is internal API. It may be removed or changed without notice in a non-breaking \ +release. Usage of this field is not advised. +""" + +const INTERNAL_ARGS_WARNING = """ +The following arguments are internal API. They may be removed or changed without notice \ +in a non-breaking release. Usage of these arguments is not advised. +""" + """ $(TYPEDEF) @@ -157,6 +167,7 @@ include("systems/parameter_buffer.jl") include("systems/abstractsystem.jl") include("systems/model_parsing.jl") include("systems/connectors.jl") +include("systems/state_machines.jl") include("systems/analysis_points.jl") include("systems/imperative_affect.jl") include("systems/callbacks.jl") @@ -296,7 +307,7 @@ export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbanc hasunit, getunit, hasconnect, getconnect, hasmisc, getmisc, state_priority export liouville_transform, change_independent_variable, substitute_component, - add_accumulations, noise_to_brownians + add_accumulations, noise_to_brownians, Girsanov_transform export PDESystem export Differential, expand_derivatives, @derivatives export Equation, ConstrainedEquation @@ -359,4 +370,6 @@ export AbstractDynamicOptProblem, JuMPDynamicOptProblem, InfiniteOptDynamicOptPr CasADiDynamicOptProblem export DynamicOptSolution +@public apply_to_variables + end # module diff --git a/src/structural_transformation/symbolics_tearing.jl b/src/structural_transformation/symbolics_tearing.jl index c3499abebc..3b6e3c7ad6 100644 --- a/src/structural_transformation/symbolics_tearing.jl +++ b/src/structural_transformation/symbolics_tearing.jl @@ -99,48 +99,6 @@ function eq_derivative!(ts::TearingState, ieq::Int; kwargs...) return eq_diff end -function tearing_sub(expr, dict, s) - expr = Symbolics.fixpoint_sub(expr, dict; operator = ModelingToolkit.Initial) - s ? simplify(expr) : expr -end - -function tearing_substitute_expr(sys::AbstractSystem, expr; simplify = false) - empty_substitutions(sys) && return expr - substitutions = get_substitutions(sys) - return tearing_sub(expr, substitutions, simplify) -end - -""" -$(TYPEDSIGNATURES) - -Like `equations(sys)`, but includes substitutions done by the tearing process. -These equations matches generated numerical code. - -See also [`equations`](@ref) and [`ModelingToolkit.get_eqs`](@ref). -""" -function full_equations(sys::AbstractSystem; simplify = false) - isempty(observed(sys)) && return equations(sys) - subs = Dict([eq.lhs => eq.rhs for eq in observed(sys)]) - neweqs = map(equations(sys)) do eq - if iscall(eq.lhs) && operation(eq.lhs) isa Union{Shift, Differential} - return tearing_sub(eq.lhs, subs, simplify) ~ tearing_sub(eq.rhs, subs, - simplify) - else - if !(eq.lhs isa Number && eq.lhs == 0) - eq = 0 ~ eq.rhs - eq.lhs - end - rhs = tearing_sub(eq.rhs, subs, simplify) - if rhs isa Symbolic - return 0 ~ rhs - else # a number - error("tearing failed because the system is singular") - end - end - eq - end - return neweqs -end - function tearing_substitution(sys::AbstractSystem; kwargs...) neweqs = full_equations(sys::AbstractSystem; kwargs...) @set! sys.eqs = neweqs diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl index a0ae5f685a..8b41823250 100644 --- a/src/systems/abstractsystem.jl +++ b/src/systems/abstractsystem.jl @@ -204,7 +204,18 @@ mutable struct Substitutions end Substitutions(subs, deps) = Substitutions(subs, deps, nothing) +""" + $(TYPEDSIGNATURES) + +Obtain the name of `sys`. +""" Base.nameof(sys::AbstractSystem) = getfield(sys, :name) +""" + $(TYPEDSIGNATURES) + +Obtain the description associated with `sys` if present, and an empty +string otherwise. +""" description(sys::AbstractSystem) = has_description(sys) ? get_description(sys) : "" #Deprecated @@ -744,6 +755,9 @@ the global structure of the system. One property to note is that if a system is complete, the system will no longer namespace its subsystems or variables, i.e. `isequal(complete(sys).v.i, v.i)`. + +This namespacing functionality can also be toggled independently of `complete` +using [`toggle_namespacing`](@ref). """ function complete( sys::AbstractSystem; split = true, flatten = true, add_initial_parameters = true) @@ -1010,6 +1024,14 @@ function Base.propertynames(sys::AbstractSystem; private = false) end end +""" + Base.getproperty(sys::AbstractSystem, name::Symbol) + +Access the subsystem, variable or analysis point of `sys` named `name`. To check if `sys` +will namespace the returned value, use `ModelingToolkit.does_namespacing(sys)`. + +See also: [`ModelingToolkit.does_namespacing`](@ref). +""" function Base.getproperty( sys::AbstractSystem, name::Symbol; namespace = does_namespacing(sys)) if has_parent(sys) && (parent = get_parent(sys); parent !== nothing) @@ -1084,8 +1106,18 @@ function Base.setproperty!(sys::AbstractSystem, prop::Symbol, val) """) end -apply_to_variables(f::F, ex) where {F} = _apply_to_variables(f, ex) -apply_to_variables(f::F, ex::Num) where {F} = wrap(_apply_to_variables(f, unwrap(ex))) +""" + $(TYPEDSIGNATURES) + +Apply function `f` to each variable in expression `ex`. `f` should be a function that takes +a variable and returns the replacement to use. A "variable" in this context refers to a +symbolic quantity created directly from a variable creation macro such as +[`Symbolics.@variables`](@ref), [`@independent_variables`](@ref), [`@parameters`](@ref), +[`@constants`](@ref) or [`@brownians`](@ref). +""" +apply_to_variables(f, ex) = _apply_to_variables(f, ex) +apply_to_variables(f, ex::Num) = wrap(_apply_to_variables(f, unwrap(ex))) +apply_to_variables(f, ex::Symbolics.Arr) = wrap(_apply_to_variables(f, unwrap(ex))) function _apply_to_variables(f::F, ex) where {F} if isvariable(ex) return f(ex) @@ -1725,6 +1757,63 @@ function equations_toplevel(sys::AbstractSystem) return get_eqs(sys) end +""" + $(TYPEDSIGNATURES) + +Recursively substitute `dict` into `expr`. Use `Symbolics.simplify` on the expression +if `simplify == true`. +""" +function substitute_and_simplify(expr, dict::AbstractDict, simplify::Bool) + expr = Symbolics.fixpoint_sub(expr, dict; operator = ModelingToolkit.Initial) + simplify ? Symbolics.simplify(expr) : expr +end + +""" + $(TYPEDSIGNATURES) + +Recursively substitute the observed equations of `sys` into `expr`. If `simplify`, call +`Symbolics.simplify` on the resultant expression. +""" +function substitute_observed(sys::AbstractSystem, expr; simplify = false) + empty_substitutions(sys) && return expr + substitutions = get_substitutions(sys) + return substitute_and_simplify(expr, substitutions, simplify) +end + +""" +$(TYPEDSIGNATURES) + +Like `equations(sys)`, but also substitutes the observed equations eliminated from the +equations during `mtkcompile`. These equations matches generated numerical code. + +See also [`equations`](@ref) and [`ModelingToolkit.get_eqs`](@ref). +""" +function full_equations(sys::AbstractSystem; simplify = false) + empty_substitutions(sys) && return equations(sys) + subs = get_substitutions(sys) + neweqs = map(equations(sys)) do eq + if iscall(eq.lhs) && operation(eq.lhs) isa Union{Shift, Differential} + return substitute_and_simplify(eq.lhs, subs, simplify) ~ substitute_and_simplify( + eq.rhs, subs, + simplify) + else + if !_iszero(eq.lhs) + eq = 0 ~ eq.rhs - eq.lhs + end + rhs = substitute_and_simplify(eq.rhs, subs, simplify) + return 0 ~ rhs + end + eq + end + return neweqs +end + +""" + $(TYPEDSIGNATURES) + +Get the flattened jumps of the system. In other words, obtain all of the jumps in `sys` and +all the subsystems of `sys` (appropriately namespaced). +""" function jumps(sys::AbstractSystem) js = get_jumps(sys) systems = get_systems(sys) @@ -1734,6 +1823,12 @@ function jumps(sys::AbstractSystem) return [js; reduce(vcat, namespace_jumps.(systems); init = [])] end +""" + $(TYPEDSIGNATURES) + +Get all of the brownian variables involved in the system `sys` and all subsystems, +appropriately namespaced. +""" function brownians(sys::AbstractSystem) bs = get_brownians(sys) systems = get_systems(sys) @@ -1743,6 +1838,12 @@ function brownians(sys::AbstractSystem) return [bs; reduce(vcat, namespace_brownians.(systems); init = [])] end +""" + $(TYPEDSIGNATURES) + +Recursively consolidate the cost vector of `sys` and all subsystems of `sys`, returning the +final scalar cost function. +""" function cost(sys::AbstractSystem) cs = get_costs(sys) consolidate = get_consolidate(sys) @@ -1772,7 +1873,12 @@ function namespace_constraints(sys) map(cstr -> namespace_constraint(cstr, sys), cstrs) end -function constraints(sys) +""" + $(TYPEDSIGNATURES) + +Get all constraints in the system `sys` and all of its subsystems, appropriately namespaced. +""" +function constraints(sys::AbstractSystem) cs = get_constraints(sys) systems = get_systems(sys) isempty(systems) ? cs : [cs; reduce(vcat, namespace_constraints.(systems))] @@ -1799,6 +1905,11 @@ function initialization_equations(sys::AbstractSystem) end end +""" + $(TYPEDSIGNATURES) + +Get the tstops present in `sys` and its subsystems, appropriately namespaced. +""" function symbolic_tstops(sys::AbstractSystem) tstops = get_tstops(sys) systems = get_systems(sys) @@ -1807,6 +1918,12 @@ function symbolic_tstops(sys::AbstractSystem) return tstops end +""" + $(TYPEDSIGNATURES) + +Obtain the preface associated with `sys` and all of its subsystems, appropriately +namespaced. +""" function preface(sys::AbstractSystem) has_preface(sys) || return nothing pre = get_preface(sys) @@ -2456,10 +2573,46 @@ function component_post_processing(expr, isconnector) end end +""" + $(TYPEDSIGNATURES) + +Mark a system constructor function as building a component. For example, + +```julia +@component function AddOne(; name) + @variables in(t) out(t) + eqs = [out ~ in + 1] + return System(eqs, t, [in, out], []; name) +end +``` + +ModelingToolkit systems are either components or connectors. Components define dynamics of +the model. Connectors are used to connect components together. See the +[Model building reference](@ref model_building_api) section of the documentation for more +information. + +See also: [`@connector`](@ref). +""" macro component(expr) esc(component_post_processing(expr, false)) end +""" + $(TYPEDSIGNATURES) + +Macro shorthand for building and compiling a system in one step. + +```julia +@mtkcompile sys = Constructor(args...; kwargs....) +``` + +Is shorthand for + +```julia +@named sys = Constructor(args...; kwargs...) +sys = mtkcompile(sys) +``` +""" macro mtkcompile(exprs...) expr = exprs[1] named_expr = ModelingToolkit.named_expr(expr) @@ -2584,6 +2737,11 @@ end hierarchy(sys::AbstractSystem; describe = false, bold = describe, kwargs...) Print a tree of a system's hierarchy of subsystems. + +# Keyword arguments + +- `describe`: Whether to also print the description of each subsystem, if present. +- `bold`: Whether to print the name of the system in **bold** font. """ function hierarchy(sys::AbstractSystem; describe = false, bold = describe, kwargs...) print_tree(sys; printnode_kw = (describe = describe, bold = bold), kwargs...) @@ -2626,9 +2784,14 @@ end """ $(TYPEDSIGNATURES) -Extend `basesys` with `sys`. +Extend `basesys` with `sys`. This can be thought of as the `merge` operation on systems. +Values in `sys` take priority over duplicates in `basesys` (for example, defaults). + By default, the resulting system inherits `sys`'s name and description. +The `&` operator can also be used for this purpose. `sys & basesys` is equivalent to +`extend(sys, basesys)`. + See also [`compose`](@ref). """ function extend(sys::AbstractSystem, basesys::AbstractSystem; @@ -2675,14 +2838,31 @@ function extend(sys::AbstractSystem, basesys::AbstractSystem; return T(args...; kwargs...) end +""" + $(TYPEDSIGNATURES) + +Extend `sys` with all systems in `basesys` in order. +""" function extend(sys, basesys::Vector{T}) where {T <: AbstractSystem} foldl(extend, basesys, init = sys) end +""" + $(TYPEDSIGNATURES) + +Syntactic sugar for `extend(sys, basesys)`. + +See also: [`extend`](@ref). +""" function Base.:(&)(sys::AbstractSystem, basesys::AbstractSystem; kwargs...) extend(sys, basesys; kwargs...) end +""" + $(TYPEDSIGNATURES) + +Syntactic sugar for `extend(sys, basesys)`. +""" function Base.:(&)( sys::AbstractSystem, basesys::Vector{T}; kwargs...) where {T <: AbstractSystem} extend(sys, basesys; kwargs...) @@ -2691,8 +2871,11 @@ end """ $(SIGNATURES) -Compose multiple systems together. The resulting system would inherit the first -system's name. +Compose multiple systems together. This adds all of `systems` as subsystems of `sys`. +The resulting system inherits the name of `sys` by default. + +The `∘` operator can also be used for this purpose. `sys ∘ basesys` is equivalent to +`compose(sys, basesys)`. See also [`extend`](@ref). """ @@ -2714,9 +2897,22 @@ function compose(sys::AbstractSystem, systems::AbstractArray; name = nameof(sys) @set! sys.ps = unique!(vcat(get_ps(sys), collect(newparams))) return sys end +""" + $(TYPEDSIGNATURES) + +Syntactic sugar for adding all systems in `syss` as the subsystems of `first(syss)`. +""" function compose(syss...; name = nameof(first(syss))) compose(first(syss), collect(syss[2:end]); name = name) end + +""" + $(TYPEDSIGNATURES) + +Syntactic sugar for `compose(sys1, sys2)`. + +See also: [`compose`](@ref). +""" Base.:(∘)(sys1::AbstractSystem, sys2::AbstractSystem) = compose(sys1, sys2) function UnPack.unpack(sys::ModelingToolkit.AbstractSystem, ::Val{p}) where {p} diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl index 2ce6356aa7..12185c943f 100644 --- a/src/systems/analysis_points.jl +++ b/src/systems/analysis_points.jl @@ -176,7 +176,7 @@ function renamespace(sys, ap::AnalysisPoint) end # create analysis points via `connect` -function Symbolics.connect(in, ap::AnalysisPoint, outs...; verbose = true) +function connect(in, ap::AnalysisPoint, outs...; verbose = true) return AnalysisPoint() ~ AnalysisPoint(in, ap.name, collect(outs); verbose) end @@ -217,11 +217,11 @@ typically is not (unless the model is an inverse model). - `verbose`: Warn if an input is connected to an output (reverse causality). Silence this warning if you are analyzing an inverse model. """ -function Symbolics.connect(in::AbstractSystem, name::Symbol, out, outs...; verbose = true) +function connect(in::AbstractSystem, name::Symbol, out, outs...; verbose = true) return AnalysisPoint() ~ AnalysisPoint(in, name, [out; collect(outs)]; verbose) end -function Symbolics.connect( +function connect( in::ConnectableSymbolicT, name::Symbol, out::ConnectableSymbolicT, outs::ConnectableSymbolicT...; verbose = true) allvars = (in, out, outs...) diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl index ba41b2b011..9761d18493 100644 --- a/src/systems/connectors.jl +++ b/src/systems/connectors.jl @@ -1,4 +1,46 @@ -using Symbolics: StateMachineOperator +""" + $(TYPEDEF) + +Struct used to represent a connection equation. A connection equation is an `Equation` +where the LHS is an empty `Connection(nothing)` and the RHS is a `Connection` containing +the connected connectors. + +For special types of connections, the LHS `Connection` can contain relevant metadata. +""" +struct Connection + systems::Any +end + +Base.broadcastable(x::Connection) = Ref(x) +Connection() = Connection(nothing) +Base.hash(c::Connection, seed::UInt) = hash(c.systems, (0xc80093537bdc1311 % UInt) ⊻ seed) +Symbolics.hide_lhs(_::Connection) = true + +""" + $(TYPEDSIGNATURES) + +Connect multiple connectors created via `@connector`. All connected connectors +must be unique. +""" +function connect(sys1::AbstractSystem, sys2::AbstractSystem, syss::AbstractSystem...) + syss = (sys1, sys2, syss...) + length(unique(nameof, syss)) == length(syss) || error("connect takes distinct systems!") + Equation(Connection(), Connection(syss)) # the RHS are connected systems +end + +function Base.show(io::IO, c::Connection) + print(io, "connect(") + if c.systems isa AbstractArray || c.systems isa Tuple + n = length(c.systems) + for (i, s) in enumerate(c.systems) + str = join(split(string(nameof(s)), NAMESPACE_SEPARATOR), '.') + print(io, str) + i != n && print(io, ", ") + end + end + print(io, ")") +end + isconnection(_) = false isconnection(_::Connection) = true """ @@ -20,6 +62,37 @@ function get_connection_type(s) getmetadata(s, VariableConnectType, Equality) end +""" + $(TYPEDSIGNATURES) + +Mark a system constructor function as building a connector. For example, + +```julia +@connector function ElectricalPin(; name, v = nothing, i = nothing) + @variables begin + v(t) = v, [description = "Potential at the pin [V]"] + i(t) = i, [connect = Flow, description = "Current flowing into the pin [A]"] + end + return System(Equation[], t, [v, i], []; name) +end +``` + +Since connectors only declare variables, the equivalent shorthand syntax can also be used: + +```julia +@connector Pin begin + v(t), [description = "Potential at the pin [V]"] + i(t), [connect = Flow, description = "Current flowing into the pin [A]"] +end +``` + +ModelingToolkit systems are either components or connectors. Components define dynamics of +the model. Connectors are used to connect components together. See the +[Model building reference](@ref model_building_api) section of the documentation for more +information. + +See also: [`@component`](@ref). +""" macro connector(expr) esc(component_post_processing(expr, true)) end @@ -63,6 +136,15 @@ is_domain_connector(s) = isconnector(s) && get_connector_type(s) === DomainConne get_systems(c::Connection) = c.systems +""" + $(TYPEDSIGNATURES) + +`instream` is used when modeling stream connections. It is only allowed to be used on +`Stream` variables. + +Refer to the [Connection semantics](@ref connect_semantics) section of the docs for more +information. +""" instream(a) = term(instream, unwrap(a), type = symtype(a)) SymbolicUtils.promote_symtype(::typeof(instream), _) = Real @@ -148,7 +230,7 @@ var1 ~ var3 # ... ``` """ -function Symbolics.connect(var1::ConnectableSymbolicT, var2::ConnectableSymbolicT, +function connect(var1::ConnectableSymbolicT, var2::ConnectableSymbolicT, vars::ConnectableSymbolicT...) allvars = (var1, var2, vars...) validate_causal_variables_connection(allvars) diff --git a/src/systems/if_lifting.jl b/src/systems/if_lifting.jl index 130a4bdab4..c2eba035ff 100644 --- a/src/systems/if_lifting.jl +++ b/src/systems/if_lifting.jl @@ -410,6 +410,11 @@ If lifting converts (nested) if statements into a series of continuous events + Lifting proceeds through the following process: * rewrite comparisons to be of the form eqn [op] 0; subtract the RHS from the LHS * replace comparisons with generated parameters; for each comparison eqn [op] 0, generate an event (dependent on op) that sets the parameter + +!!! warn + + This is an experimental simplification pass. It may have bugs. Please open issues with + MWEs for any bugs encountered while using this. """ function IfLifting(sys::System) if !is_time_dependent(sys) diff --git a/src/systems/state_machines.jl b/src/systems/state_machines.jl new file mode 100644 index 0000000000..347f92e6f8 --- /dev/null +++ b/src/systems/state_machines.jl @@ -0,0 +1,155 @@ +_nameof(s) = nameof(s) +_nameof(s::Union{Int, Symbol}) = s +abstract type StateMachineOperator end +Base.broadcastable(x::StateMachineOperator) = Ref(x) +Symbolics.hide_lhs(_::StateMachineOperator) = true +struct InitialState <: StateMachineOperator + s::Any +end +Base.show(io::IO, s::InitialState) = print(io, "initial_state(", _nameof(s.s), ")") +initial_state(s) = Equation(InitialState(nothing), InitialState(s)) + +Base.@kwdef struct Transition{A, B, C} <: StateMachineOperator + from::A = nothing + to::B = nothing + cond::C = nothing + immediate::Bool = true + reset::Bool = true + synchronize::Bool = false + priority::Int = 1 + function Transition(from, to, cond, immediate, reset, synchronize, priority) + cond = unwrap(cond) + new{typeof(from), typeof(to), typeof(cond)}(from, to, cond, immediate, + reset, synchronize, + priority) + end +end +function Base.:(==)(transition1::Transition, transition2::Transition) + transition1.from == transition2.from && + transition1.to == transition2.to && + isequal(transition1.cond, transition2.cond) && + transition1.immediate == transition2.immediate && + transition1.reset == transition2.reset && + transition1.synchronize == transition2.synchronize && + transition1.priority == transition2.priority +end + +""" + transition(from, to, cond; immediate::Bool = true, reset::Bool = true, synchronize::Bool = false, priority::Int = 1) + +Create a transition from state `from` to state `to` that is enabled when transitioncondition `cond` evaluates to `true`. + +# Arguments: +- `from`: The source state of the transition. +- `to`: The target state of the transition. +- `cond`: A transition condition that evaluates to a Bool, such as `ticksInState() >= 2`. +- `immediate`: If `true`, the transition will fire at the same tick as it becomes true, if `false`, the actions of the state are evaluated first, and the transition fires during the next tick. +- `reset`: If true, the destination state `to` is reset to its initial condition when the transition fires. +- `synchronize`: If true, the transition will only fire if all sub-state machines in the source state are in their final (terminal) state. A final state is one that has no outgoing transitions. +- `priority`: If a state has more than one outgoing transition, all outgoing transitions must have a unique priority. The transitions are evaluated in priority order, i.e., the transition with priority 1 is evaluated first. +""" +function transition(from, to, cond; + immediate::Bool = true, reset::Bool = true, synchronize::Bool = false, + priority::Int = 1) + Equation( + Transition(), Transition(; from, to, cond, immediate, reset, + synchronize, priority)) +end +function Base.show(io::IO, s::Transition) + print(io, _nameof(s.from), " → ", _nameof(s.to), " if (", s.cond, ") [") + print(io, "immediate: ", Int(s.immediate), ", ") + print(io, "reset: ", Int(s.reset), ", ") + print(io, "sync: ", Int(s.synchronize), ", ") + print(io, "prio: ", s.priority, "]") +end + +function activeState end +function entry end +function ticksInState end +function timeInState end + +for (s, T) in [(:timeInState, :Real), + (:ticksInState, :Integer), + (:entry, :Bool), + (:activeState, :Bool)] + seed = hash(s) + @eval begin + $s(x) = wrap(term($s, x)) + SymbolicUtils.promote_symtype(::typeof($s), _...) = $T + function SymbolicUtils.show_call(io, ::typeof($s), args) + if isempty(args) + print(io, $s, "()") + else + arg = only(args) + print(io, $s, "(", arg isa Number ? arg : nameof(arg), ")") + end + end + end + if s != :activeState + @eval $s() = wrap(term($s)) + end +end + +@doc """ + timeInState() + timeInState(state) + +Get the time (in seconds) spent in a state in a finite state machine. + +When used to query the time spent in the enclosing state, the method without arguments is used, i.e., +```julia +@mtkmodel FSM begin + ... + @equations begin + var(k+1) ~ timeInState() >= 2 ? 0.0 : var(k) + end +end +``` + +If used to query the residence time of another state, the state is passed as an argument. + +This operator can be used in both equations and transition conditions. + +See also [`ticksInState`](@ref) and [`entry`](@ref) +""" timeInState + +@doc """ + ticksInState() + ticksInState(state) + +Get the number of ticks spent in a state in a finite state machine. + +When used to query the number of ticks spent in the enclosing state, the method without arguments is used, i.e., +```julia +@mtkmodel FSM begin + ... + @equations begin + var(k+1) ~ ticksInState() >= 2 ? 0.0 : var(k) + end +end +``` + +If used to query the number of ticks in another state, the state is passed as an argument. + +This operator can be used in both equations and transition conditions. + +See also [`timeInState`](@ref) and [`entry`](@ref) +""" ticksInState + +@doc """ + entry() + entry(state) + +When used in a finite-state machine, this operator returns true at the first tick when the state is active, and false otherwise. + +When used to query the entry of the enclosing state, the method without arguments is used, when used to query the entry of another state, the state is passed as an argument. + +This can be used to perform a unique action when entering a state. +""" +entry + +@doc """ + activeState(state) + +When used in a finite state machine, this operator returns `true` if the queried state is active and false otherwise. +""" activeState diff --git a/src/systems/system.jl b/src/systems/system.jl index 5da5d9db79..872cf5fc95 100644 --- a/src/systems/system.jl +++ b/src/systems/system.jl @@ -11,50 +11,237 @@ end const MetadataT = Base.ImmutableDict{DataType, Any} +""" + $(TYPEDEF) + +A symbolic representation of a numerical system to be solved. This is a recursive +tree-like data structure - each system can contain additional subsystems. As such, +it implements the `AbstractTrees.jl` interface to enable exploring the hierarchical +structure. + +# Fields + +$(TYPEDFIELDS) +""" struct System <: AbstractSystem + """ + $INTERNAL_FIELD_WARNING + A unique integer tag for the system. + """ tag::UInt + """ + The equations of the system. + """ eqs::Vector{Equation} # nothing - no noise # vector - diagonal noise # matrix - generic form # column matrix - scalar noise + """ + The noise terms for each equation of the system. This field is only used for flattened + systems. To represent noise in a hierarchical system, use brownians. In a system with + `N` equations and `K` independent brownian variables, this should be an `N x K` + matrix. In the special case where `N == K` and each equation has independent noise, + this noise matrix is diagonal. Diagonal noise can be specified by providing an `N` + length vector. If this field is `nothing`, the system does not have noise. + """ noise_eqs::Union{Nothing, AbstractVector, AbstractMatrix} + """ + Jumps associated with the system. Each jump can be a `VariableRateJump`, + `ConstantRateJump` or `MassActionJump`. See `JumpProcesses.jl` for more information. + """ jumps::Vector{JumpType} + """ + The constraints of the system. This can be used to represent the constraints in an + optimal-control problem or boundary-value differential equation, or the constraints + in a constrained optimization. + """ constraints::Vector{Union{Equation, Inequality}} + """ + The costs of the system. This can be the cost in an optimal-control problem, or the + loss of an optimization problem. Scalar loss values must also be provided as a single- + element vector. + """ costs::Vector{<:Union{BasicSymbolic, Real}} + """ + A function which combines costs into a scalar value. This should take two arguments, + the `costs` of this system and the consolidated costs of all subsystems in the order + they are present in the `systems` field. It should return a scalar cost that combines + all of the individual values. This defaults to a function that simply sums all cost + values. + """ consolidate::Any + """ + The variables being solved for by this system. For example, in a differential equation + system, this contains the dependent variables. + """ unknowns::Vector + """ + The parameters of the system. Parameters can either be variables that parameterize the + problem being solved for (e.g. the spring constant of a mass-spring system) or + additional unknowns not part of the main dynamics of the system (e.g. discrete/clocked + variables in a hybrid ODE). + """ ps::Vector + """ + The brownian variables of the system, created via `@brownians`. Each brownian variable + represents an independent noise. A system with brownians cannot be simulated directly. + It needs to be compiled using `mtkcompile` into `noise_eqs`. + """ brownians::Vector + """ + The independent variable for a time-dependent system, or `nothing` for a time-independent + system. + """ iv::Union{Nothing, BasicSymbolic{Real}} + """ + Equations that compute variables of a system that have been eliminated from the set of + unknowns by `mtkcompile`. More generally, this contains all variables that can be + computed from the unknowns and parameters and do not need to be solved for. Such + variables are termed as "observables". Each equation must be of the form + `observable ~ expression` and observables cannot appear on the LHS of multiple + equations. Equations must be sorted such that every observable appears on + the left hand side of an equation before it appears on the right hand side of any other + equation. + """ observed::Vector{Equation} parameter_dependencies::Vector{Equation} + """ + $INTERNAL_FIELD_WARNING + A mapping from the name of a variable to the actual symbolic variable in the system. + This is used to enable `getproperty` syntax to access variables of a system. + """ var_to_name::Dict{Symbol, Any} + """ + The name of the system. + """ name::Symbol + """ + An optional description for the system. + """ description::String + """ + Default values that variables (unknowns/observables/parameters) should take when + constructing a numerical problem from the system. These values can be overridden + by initial values provided to the problem constructor. Defaults of parent systems + take priority over those in child systems. + """ defaults::Dict + """ + Guess values for variables of a system that are solved for during initialization. + """ guesses::Dict + """ + A list of subsystems of this system. Used for hierarchically building models. + """ systems::Vector{System} + """ + Equations that must be satisfied during initialization of the numerical problem created + from this system. For time-dependent systems, these equations are not valid after the + initial time. + """ initialization_eqs::Vector{Equation} + """ + Symbolic representation of continuous events in a dynamical system. See + [`SymbolicContinuousCallback`](@ref). + """ continuous_events::Vector{SymbolicContinuousCallback} + """ + Symbolic representation of discrete events in a dynamica system. See + [`SymbolicDiscreteCallback`](@ref). + """ discrete_events::Vector{SymbolicDiscreteCallback} + """ + $INTERNAL_FIELD_WARNING + If this system is a connector, the type of connector it is. + """ connector_type::Any + """ + A map from expressions that must be through throughout the solution process to an + associated error message. By default these assertions cause the generated code to + output `NaN`s if violated, but can be made to error using `debug_system`. + """ assertions::Dict{BasicSymbolic, String} + """ + The metadata associated with this system, as a `Base.ImmutableDict`. This follows + the same interface as SymbolicUtils.jl. Metadata can be queried and updated using + `SymbolicUtils.getmetadata` and `SymbolicUtils.setmetadata` respectively. + """ metadata::MetadataT + """ + $INTERNAL_FIELD_WARNING + Metadata added by the `@mtkmodel` macro. + """ gui_metadata::Any # ? + """ + Whether the system contains delay terms. This is inferred from the equations, but + can also be provided explicitly. + """ is_dde::Bool + """ + Extra time points for the integrator to stop at. These can be numeric values, + or expressions of parameters and time. + """ tstops::Vector{Any} + """ + The `TearingState` of the system post-simplification with `mtkcompile`. + """ tearing_state::Any + """ + Whether the system namespaces variables accessed via `getproperty`. `complete`d systems + do not namespace, but this flag can be toggled independently of `complete` using + `toggle_namespacing`. + """ namespacing::Bool + """ + Whether the system is marked as "complete". Completed systems cannot be used as + subsystems. + """ complete::Bool + """ + $INTERNAL_FIELD_WARNING + For systems simplified or completed with `split = true` (the default) this contains an + `IndexCache` which aids in symbolic indexing. If this field is `nothing`, the system is + either not completed, or completed with `split = false`. + """ index_cache::Union{Nothing, IndexCache} + """ + $INTERNAL_FIELD_WARNING + Connections that should be ignored because they were removed by an analysis point + transformation. The first element of the tuple contains all such "standard" connections + (ones between connector systems) and the second contains all such causal variable + connections. + """ ignored_connections::Union{ Nothing, Tuple{Vector{IgnoredAnalysisPoint}, Vector{IgnoredAnalysisPoint}}} + """ + `SymbolicUtils.Code.Assignment`s to prepend to all code generated from this system. + """ preface::Any + """ + After simplification with `mtkcompile`, this field contains the unsimplified system + with the hierarchical structure. There may be multiple levels of `parent`s. The root + parent is used for accessing variables via `getproperty` syntax. + """ parent::Union{Nothing, System} + """ + A custom initialization system to use if no initial conditions are provided for the + unknowns or observables of this system. + """ initializesystem::Union{Nothing, System} + """ + Whether the current system is an initialization system. + """ is_initializesystem::Bool + """ + $INTERNAL_FIELD_WARNING + Whether the system has been simplified by `mtkcompile`. + """ isscheduled::Bool + """ + $INTERNAL_FIELD_WARNING + The `Schedule` containing additional information about the simplified system. + """ schedule::Union{Schedule, Nothing} function System( @@ -114,6 +301,22 @@ function default_consolidate(costs, subcosts) return reduce(+, costs; init = 0.0) + reduce(+, subcosts; init = 0.0) end +""" + $(TYPEDSIGNATURES) + +Construct a system using the given equations `eqs`, independent variable `iv` (`nothing`) +for time-independent systems, unknowns `dvs`, parameters `ps` and brownian variables +`brownians`. + +## Keyword Arguments + +- `discover_from_metadata`: Whether to parse metadata of unknowns and parameters of the + system to obtain defaults and/or guesses. +- `checks`: Whether to perform sanity checks on the passed values. + +All other keyword arguments are named identically to the corresponding fields in +[`System`](@ref). +""" function System(eqs::Vector{Equation}, iv, dvs, ps, brownians = []; constraints = Union{Equation, Inequality}[], noise_eqs = nothing, jumps = [], costs = BasicSymbolic[], consolidate = default_consolidate, @@ -207,10 +410,23 @@ function System(eqs::Vector{Equation}, iv, dvs, ps, brownians = []; initializesystem, is_initializesystem; checks) end +""" + $(TYPEDSIGNATURES) + +Create a time-independent [`System`](@ref) with the given equations `eqs`, unknowns `dvs` +and parameters `ps`. +""" function System(eqs::Vector{Equation}, dvs, ps; kwargs...) System(eqs, nothing, dvs, ps; kwargs...) end +""" + $(TYPEDSIGNATURES) + +Create a time-dependent system with the given equations `eqs` and independent variable `iv`. +Discover variables, parameters and brownians in the system by parsing the equations and +other symbolic expressions passed to the system. +""" function System(eqs::Vector{Equation}, iv; kwargs...) iv === nothing && return System(eqs; kwargs...) @@ -299,6 +515,13 @@ function System(eqs::Vector{Equation}, iv; kwargs...) eqs, iv, collect(allunknowns), collect(new_ps), collect(brownians); kwargs...) end +""" + $(TYPEDSIGNATURES) + +Create a time-independent system with the given equations `eqs`. Discover variables and +parameters in the system by parsing the equations and other symbolic expressions passed to +the system. +""" function System(eqs::Vector{Equation}; kwargs...) eqs = collect(eqs) @@ -329,6 +552,11 @@ function System(eqs::Vector{Equation}; kwargs...) return System(eqs, nothing, collect(allunknowns), collect(new_ps); kwargs...) end +""" + $(TYPEDSIGNATURES) + +Create a `System` with a single equation `eq`. +""" System(eq::Equation, args...; kwargs...) = System([eq], args...; kwargs...) function gather_array_params(ps) @@ -465,6 +693,12 @@ function _check_if_dde(eqs, iv, subsystems) return is_dde end +""" + $(TYPEDSIGNATURES) + +Flatten the hierarchical structure of a system, collecting all equations, unknowns, etc. +into one top-level system after namespacing appropriately. +""" function flatten(sys::System, noeqs = false) systems = get_systems(sys) isempty(systems) && return sys @@ -635,11 +869,23 @@ function Base.hash(sys::System, h::UInt) return h end +""" + $(TYPEDSIGNATURES) + +Get the metadata associated with key `k` in system `sys` or `default` if it does not exist. +""" function SymbolicUtils.getmetadata(sys::AbstractSystem, k::DataType, default) meta = get_metadata(sys) return get(meta, k, default) end +""" + $(TYPEDSIGNATURES) + +Set the metadata associated with key `k` in system `sys` to value `v`. This is an +out-of-place operation, and will return a shallow copy of `sys` with the appropriate +metadata values. +""" function SymbolicUtils.setmetadata(sys::AbstractSystem, k::DataType, v) meta = get_metadata(sys) meta = Base.ImmutableDict(meta, k => v)::MetadataT @@ -653,6 +899,13 @@ function check_complete(sys::System, obj) iscomplete(sys) || throw(SystemNotCompleteError(obj)) end +""" + $(TYPEDSIGNATURES) + +Convert a time-dependent system `sys` to a time-independent system of nonlinear +equations that solve for the steady state of the system where `D(x)` is zero for +each continuous variable `x`. +""" function NonlinearSystem(sys::System) if !is_time_dependent(sys) throw(ArgumentError("`NonlinearSystem` constructor expects a time-dependent `System`")) @@ -678,6 +931,11 @@ end # Utility constructors ######## +""" + $(METHODLIST) + +Construct a [`System`](@ref) to solve an optimization problem with the given scalar cost. +""" function OptimizationSystem(cost; kwargs...) return System(Equation[]; costs = [cost], kwargs...) end @@ -694,6 +952,11 @@ function OptimizationSystem(cost::Array, dvs, ps; kwargs...) return System(Equation[], nothing, dvs, ps; costs = vec(cost), kwargs...) end +""" + $(METHODLIST) + +Construct a [`System`](@ref) to solve a system of jump equations. +""" function JumpSystem(jumps, iv; kwargs...) mask = isa.(jumps, Equation) eqs = Vector{Equation}(jumps[mask]) @@ -708,6 +971,11 @@ function JumpSystem(jumps, iv, dvs, ps; kwargs...) return System(eqs, iv, dvs, ps; jumps, kwargs...) end +""" + $(METHODLIST) + +Construct a system of equations with associated noise terms. +""" function SDESystem(eqs::Vector{Equation}, noise, iv; is_scalar_noise = false, kwargs...) if is_scalar_noise if !(noise isa Vector) diff --git a/src/systems/systems.jl b/src/systems/systems.jl index 13a62f7915..63585541e0 100644 --- a/src/systems/systems.jl +++ b/src/systems/systems.jl @@ -9,10 +9,18 @@ end """ $(SIGNATURES) -Structurally simplify algebraic equations in a system and compute the -topological sort of the observed equations in `sys`. +Compile the given system into a form that ModelingToolkit can generate code for. Also +performs a variety of symbolic-numeric enhancements. For ODEs, this includes processes +such as order reduction, index reduction, alias elimination and tearing. A subset of the +unknowns of the system may be eliminated as observables, eliminating the need for the +numerical solver to solve for these variables. + +Does not rely on metadata to identify variables/parameters/brownians. Instead, queries +the system for which symbolic quantites belong to which category. Any variables not +present in the equations of the system will be removed in this process. + +# Keyword Arguments -### Optional Keyword Arguments: + When `simplify=true`, the `simplify` function will be applied during the tearing process. + `allow_symbolic=false`, `allow_parameter=true`, and `conservative=false` limit the coefficient types during tearing. In particular, `conservative=true` limits tearing to only solve for trivial linear systems where the coefficient has the absolute value of ``1``. + `fully_determined=true` controls whether or not an error will be thrown if the number of equations don't match the number of inputs, outputs, and equations. @@ -147,7 +155,7 @@ function __mtkcompile(sys::AbstractSystem; simplify = false, is_scalar_noise = false end - noise_eqs = StructuralTransformations.tearing_substitute_expr(ode_sys, noise_eqs) + noise_eqs = substitute_observed(ode_sys, noise_eqs) ssys = System(Vector{Equation}(full_equations(ode_sys)), get_iv(ode_sys), unknowns(ode_sys), parameters(ode_sys); noise_eqs, name = nameof(ode_sys), observed = observed(ode_sys), defaults = defaults(sys), diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 16551008e9..4b01917bcf 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -1,11 +1,11 @@ using DataStructures -using Symbolics: linear_expansion, unwrap, Connection +using Symbolics: linear_expansion, unwrap using SymbolicUtils: iscall, operation, arguments, Symbolic using SymbolicUtils: quick_cancel, maketerm using ..ModelingToolkit import ..ModelingToolkit: isdiffeq, var_from_nested_derivative, vars!, flatten, value, InvalidSystemException, isdifferential, _iszero, - isparameter, + isparameter, Connection, independent_variables, SparseMatrixCLIL, AbstractSystem, equations, isirreducible, input_timedomain, TimeDomain, InferredTimeDomain, diff --git a/src/utils.jl b/src/utils.jl index 4a183fd83f..4dea0e6cb6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -750,10 +750,21 @@ end isarray(x) = x isa AbstractArray || x isa Symbolics.Arr +""" + $(TYPEDSIGNATURES) + +Check if any variables were eliminated from the system as part of `mtkcompile`. +""" function empty_substitutions(sys) isempty(observed(sys)) end +""" + $(TYPEDSIGNATURES) + +Get a dictionary mapping variables eliminated from the system during `mtkcompile` to the +expressions used to calculate them. +""" function get_substitutions(sys) Dict([eq.lhs => eq.rhs for eq in observed(sys)]) end diff --git a/src/variables.jl b/src/variables.jl index eb2e54a268..887765a192 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -91,8 +91,50 @@ end ### Connect abstract type AbstractConnectType end +""" + $(TYPEDEF) + +Flag which is meant to be passed to the `connect` metadata of a variable to affect how it +behaves when the connector it is in is part of a `connect` equation. `Equality` is the +default value and such variables when connected are made equal. For example, electric +potential is equated at a junction. + +For more information, refer to the [Connection semantics](@ref connect_semantics) section +of the docs. + +See also: [`connect`](@ref), [`@connector`](@ref), [`Flow`](@ref), +[`Stream`](@ref). +""" struct Equality <: AbstractConnectType end # Equality connection +""" + $(TYPEDEF) + +Flag which is meant to be passed to the `connect` metadata of a variable to affect how it +behaves when the connector it is in is part of a `connect` equation. `Flow` denotes that +the sum of marked variable in all connectors in the connection set must sum to zero. For +example, electric current sums to zero at a junction (assuming appropriate signs are used +for current flowing in and out of the function). + +For more information, refer to the [Connection semantics](@ref connect_semantics) section +of the docs. + +See also: [`connect`](@ref), [`@connector`](@ref), [`Equality`](@ref), +[`Stream`](@ref). +""" struct Flow <: AbstractConnectType end # sum to 0 +""" + $(TYPEDEF) + +Flag which is meant to be passed to the `connect` metadata of a variable to affect how it +behaves when the connector it is in is part of a `connect` equation. `Stream` denotes that +the variable is part of a special stream connector. + +For more information, refer to the [Connection semantics](@ref connect_semantics) section +of the docs. + +See also: [`connect`](@ref), [`@connector`](@ref), [`Equality`](@ref), +[`Flow`](@ref). +""" struct Stream <: AbstractConnectType end # special stream connector """