search_index.js

var documenterSearchIndex = {"docs":
[{"location":"basics/AbstractSystem/#The-AbstractSystem-Interface-1","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"","category":"section"},{"location":"basics/AbstractSystem/#Overview-1","page":"The AbstractSystem Interface","title":"Overview","text":"","category":"section"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"The AbstractSystem interface is the core of the system level of ModelingToolkit.jl. It establishes a common set of functionality that is used between systems from ODEs and chemical reactions, allowing users to have a common framework for model manipulation and compilation.","category":"page"},{"location":"basics/AbstractSystem/#Constructors-and-Naming-1","page":"The AbstractSystem Interface","title":"Constructors and Naming","text":"","category":"section"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"The AbstractSystem interface has a consistent method for constructing systems. Generally it follows the order of:","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"Equations\nIndependent Variables\nDependent Variables (or States)\nParameters","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"All other pieces are handled via keyword arguments. AbstractSystems share the same keyword arguments, which are:","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"system: This is used for specifying subsystems for hierarchical modeling with reusable components. For more information, see the components page\nDefaults: Keyword arguments like defaults are used for specifying default values which are used. If a value is not given at the SciMLProblem construction time, its numerical value will be the default.","category":"page"},{"location":"basics/AbstractSystem/#Composition-and-Accessor-Functions-1","page":"The AbstractSystem Interface","title":"Composition and Accessor Functions","text":"","category":"section"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"Each AbstractSystem has lists of variables in context, such as distinguishing parameters vs states. In addition, an AbstractSystem also can hold other AbstractSystem types. Direct accessing of the values, such as sys.states, gives the immediate list, while the accessor functions states(sys) gives the total set, which includes that of all systems held inside.","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"The values which are common to all AbstractSystems are:","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"equations(sys): All equations that define the system and its subsystems.\nstates(sys): All the states in the system and its subsystems.\nparameters(sys): All parameters of the system and its subsystems.\nnameof(sys): The name of the current-level system.\nget_eqs(sys): Equations that define the current-level system.\nget_states(sys): States that are in the current-level system.\nget_ps(sys): Parameters that are in the current-level system.\nget_systems(sys): Subsystems of the current-level system.","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"Optionally, a system could have:","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"observed(sys): All observed equations of the system and its subsystems.\nget_observed(sys): Observed equations of the current-level system.\nget_defaults(sys): A Dict that maps variables into their default values.\nindependent_variable(sys): The independent variable of a system.\nget_noiseeqs(sys): Noise equations of the current-level system.","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"Note that there's get_iv(sys), but it is not advised to use, since it errors when the system has no field iv. independent_variable(sys) returns nothing for NonlinearSystems.","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"A system could also have caches:","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"get_jac(sys): The Jacobian of a system.\nget_tgrad(sys): The gradient with respect to time of a system.","category":"page"},{"location":"basics/AbstractSystem/#Transformations-1","page":"The AbstractSystem Interface","title":"Transformations","text":"","category":"section"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"Transformations are functions which send a valid AbstractSystem definition to another AbstractSystem. These are passes, like optimizations (e.g., Block-Lower Triangle transformations), or changes to the representation, which allow for alternative numerical methods to be utilized on the model (e.g., DAE index reduction).","category":"page"},{"location":"basics/AbstractSystem/#Analyses-1","page":"The AbstractSystem Interface","title":"Analyses","text":"","category":"section"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"Analyses are functions on a system which return information about the corresponding properties, like whether its parameters are structurally identifiable, or whether it's linear.","category":"page"},{"location":"basics/AbstractSystem/#Function-Calculation-and-Generation-1","page":"The AbstractSystem Interface","title":"Function Calculation and Generation","text":"","category":"section"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"The calculation and generation functions allow for calculating additional quantities to enhance the numerical methods applied to the resulting system. The calculations, like calculate_jacobian, generate ModelingToolkit IR for the Jacobian of the system, while the generations, like generate_jacobian, generate compiled output for the numerical solvers by applying build_function to the generated code. Additionally, many systems have function-type outputs, which cobble together the generation functionality for a system, for example, ODEFunction can be used to generate a DifferentialEquations-based ODEFunction with compiled version of the ODE itself, the Jacobian, the mass matrix, etc.","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"Below are the possible calculation and generation functions:","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"calculate_tgrad\ncalculate_gradient\ncalculate_jacobian\ncalculate_factorized_W\ncalculate_hessian\ngenerate_tgrad\ngenerate_gradient\ngenerate_jacobian\ngenerate_factorized_W\ngenerate_hessian","category":"page"},{"location":"basics/AbstractSystem/#ModelingToolkit.calculate_tgrad","page":"The AbstractSystem Interface","title":"ModelingToolkit.calculate_tgrad","text":"calculate_tgrad(sys::AbstractSystem)\n\nCalculate the time gradient of a system.\n\nReturns a vector of Num instances. The result from the first call will be cached in the system object.\n\n\n\n\n\n","category":"function"},{"location":"basics/AbstractSystem/#ModelingToolkit.calculate_gradient","page":"The AbstractSystem Interface","title":"ModelingToolkit.calculate_gradient","text":"calculate_gradient(sys::AbstractSystem)\n\nCalculate the gradient of a scalar system.\n\nReturns a vector of Num instances. The result from the first call will be cached in the system object.\n\n\n\n\n\n","category":"function"},{"location":"basics/AbstractSystem/#ModelingToolkit.calculate_jacobian","page":"The AbstractSystem Interface","title":"ModelingToolkit.calculate_jacobian","text":"calculate_jacobian(sys::AbstractSystem)\n\nCalculate the jacobian matrix of a system.\n\nReturns a matrix of Num instances. The result from the first call will be cached in the system object.\n\n\n\n\n\n","category":"function"},{"location":"basics/AbstractSystem/#ModelingToolkit.calculate_factorized_W","page":"The AbstractSystem Interface","title":"ModelingToolkit.calculate_factorized_W","text":"calculate_factorized_W(sys::AbstractSystem)\n\nCalculate the factorized W-matrix of a system.\n\nReturns a matrix of Num instances. The result from the first call will be cached in the system object.\n\n\n\n\n\n","category":"function"},{"location":"basics/AbstractSystem/#ModelingToolkit.calculate_hessian","page":"The AbstractSystem Interface","title":"ModelingToolkit.calculate_hessian","text":"calculate_hessian(sys::AbstractSystem)\n\nCalculate the hessian matrix of a scalar system.\n\nReturns a matrix of Num instances. The result from the first call will be cached in the system object.\n\n\n\n\n\n","category":"function"},{"location":"basics/AbstractSystem/#ModelingToolkit.generate_tgrad","page":"The AbstractSystem Interface","title":"ModelingToolkit.generate_tgrad","text":"generate_tgrad(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)\n\nGenerates a function for the time gradient of a system. Extra arguments control the arguments to the internal build_function call.\n\n\n\n\n\n","category":"function"},{"location":"basics/AbstractSystem/#ModelingToolkit.generate_gradient","page":"The AbstractSystem Interface","title":"ModelingToolkit.generate_gradient","text":"generate_gradient(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; kwargs...)\n\nGenerates a function for the gradient of a system. Extra arguments control the arguments to the internal build_function call.\n\n\n\n\n\n","category":"function"},{"location":"basics/AbstractSystem/#ModelingToolkit.generate_jacobian","page":"The AbstractSystem Interface","title":"ModelingToolkit.generate_jacobian","text":"generate_jacobian(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; sparse = false, kwargs...)\n\nGenerates a function for the jacobian matrix matrix of a system. Extra arguments control the arguments to the internal build_function call.\n\n\n\n\n\n","category":"function"},{"location":"basics/AbstractSystem/#ModelingToolkit.generate_factorized_W","page":"The AbstractSystem Interface","title":"ModelingToolkit.generate_factorized_W","text":"generate_factorized_W(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; sparse = false, kwargs...)\n\nGenerates a function for the factorized W-matrix matrix of a system. Extra arguments control the arguments to the internal build_function call.\n\n\n\n\n\n","category":"function"},{"location":"basics/AbstractSystem/#ModelingToolkit.generate_hessian","page":"The AbstractSystem Interface","title":"ModelingToolkit.generate_hessian","text":"generate_hessian(sys::AbstractSystem, dvs = states(sys), ps = parameters(sys), expression = Val{true}; sparse = false, kwargs...)\n\nGenerates a function for the hessian matrix matrix of a system. Extra arguments control the arguments to the internal build_function call.\n\n\n\n\n\n","category":"function"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"Additionally, jacobian_sparsity(sys) and hessian_sparsity(sys) exist on the appropriate systems for fast generation of the sparsity patterns via an abstract interpretation without requiring differentiation.","category":"page"},{"location":"basics/AbstractSystem/#Problem-Constructors-1","page":"The AbstractSystem Interface","title":"Problem Constructors","text":"","category":"section"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"At the end, the system types have DEProblem constructors, like ODEProblem, which allow for directly generating the problem types required for numerical methods. The first argument is always the AbstractSystem, and the proceeding arguments match the argument order of their original constructors. Whenever an array would normally be provided, such as u0 the initial condition of an ODEProblem, it is instead replaced with a variable map, i.e., an array of pairs var=>value, which allows the user to designate the values without having to know the order that ModelingToolkit is internally using.","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"For the value maps, the parameters are allowed to be functions of each other, and value maps of states can be functions of the parameters, i.e. you can do:","category":"page"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"u0 = [\n  lorenz1.x => 2.0\n  lorenz2.x => lorenz1.x * lorenz1.p\n]","category":"page"},{"location":"basics/AbstractSystem/#Default-Value-Handling-1","page":"The AbstractSystem Interface","title":"Default Value Handling","text":"","category":"section"},{"location":"basics/AbstractSystem/#","page":"The AbstractSystem Interface","title":"The AbstractSystem Interface","text":"The AbstractSystem types allow for specifying default values, for example defaults inside of them. At problem construction time, these values are merged into the value maps, where for any repeats the value maps override the default. In addition, defaults of a higher level in the system override the defaults of a lower level in the system.","category":"page"},{"location":"systems/ReactionSystem/#ReactionSystem-1","page":"ReactionSystem","title":"ReactionSystem","text":"","category":"section"},{"location":"systems/ReactionSystem/#","page":"ReactionSystem","title":"ReactionSystem","text":"A ReactionSystem represents a system of chemical reactions. Conversions are provided to generate corresponding chemical reaction ODE models, chemical Langevin equation SDE models, and stochastic chemical kinetics jump process models. As a simple example, the code below creates a SIR model, and solves the corresponding ODE, SDE, and jump process models.","category":"page"},{"location":"systems/ReactionSystem/#","page":"ReactionSystem","title":"ReactionSystem","text":"using ModelingToolkit, OrdinaryDiffEq, StochasticDiffEq, DiffEqJump\r\n@parameters β γ t\r\n@variables S(t) I(t) R(t)\r\n\r\nrxs = [Reaction(β, [S,I], [I], [1,1], [2])\r\n       Reaction(γ, [I], [R])]\r\nrs  = ReactionSystem(rxs, t, [S,I,R], [β,γ])\r\n\r\nu₀map    = [S => 999.0, I => 1.0, R => 0.0]\r\nparammap = [β => 1/10000, γ => 0.01]\r\ntspan    = (0.0, 250.0)\r\n\r\n# solve as ODEs\r\nodesys = convert(ODESystem, rs)\r\noprob = ODEProblem(odesys, u₀map, tspan, parammap)\r\nsol = solve(oprob, Tsit5())\r\n\r\n# solve as SDEs\r\nsdesys = convert(SDESystem, rs)\r\nsprob = SDEProblem(sdesys, u₀map, tspan, parammap)\r\nsol = solve(sprob, EM(), dt=.01)\r\n\r\n# solve as jump process\r\njumpsys = convert(JumpSystem, rs)\r\nu₀map    = [S => 999, I => 1, R => 0]\r\ndprob = DiscreteProblem(jumpsys, u₀map, tspan, parammap)\r\njprob = JumpProblem(jumpsys, dprob, Direct())\r\nsol = solve(jprob, SSAStepper())","category":"page"},{"location":"systems/ReactionSystem/#System-Constructors-1","page":"ReactionSystem","title":"System Constructors","text":"","category":"section"},{"location":"systems/ReactionSystem/#","page":"ReactionSystem","title":"ReactionSystem","text":"Reaction\r\nReactionSystem","category":"page"},{"location":"systems/ReactionSystem/#ModelingToolkit.Reaction","page":"ReactionSystem","title":"ModelingToolkit.Reaction","text":"struct Reaction{S, T<:Number}\n\nOne chemical reaction.\n\nFields\n\nrate\nThe rate function (excluding mass action terms).\nsubstrates\nReaction substrates.\nproducts\nReaction products.\nsubstoich\nThe stoichiometric coefficients of the reactants.\nprodstoich\nThe stoichiometric coefficients of the products.\nnetstoich\nThe net stoichiometric coefficients of all species changed by the reaction.\nonly_use_rate\nfalse (default) if rate should be multiplied by mass action terms to give the rate law. true if rate represents the full reaction rate law.\n\nconnection_type\ntype: type of the system\n\nExamples\n\nusing ModelingToolkit\n@parameters t k[1:20]\n@variables A(t) B(t) C(t) D(t)\nrxs = [Reaction(k[1], nothing, [A]),            # 0 -> A\n       Reaction(k[2], [B], nothing),            # B -> 0\n       Reaction(k[3],[A],[C]),                  # A -> C\n       Reaction(k[4], [C], [A,B]),              # C -> A + B\n       Reaction(k[5], [C], [A], [1], [2]),      # C -> A + A\n       Reaction(k[6], [A,B], [C]),              # A + B -> C\n       Reaction(k[7], [B], [A], [2], [1]),      # 2B -> A\n       Reaction(k[8], [A,B], [A,C]),            # A + B -> A + C\n       Reaction(k[9], [A,B], [C,D]),            # A + B -> C + D\n       Reaction(k[10], [A], [C,D], [2], [1,1]), # 2A -> C + D\n       Reaction(k[11], [A], [A,B], [2], [1,1]), # 2A -> A + B\n       Reaction(k[12], [A,B,C], [C,D], [1,3,4], [2, 3]),          # A+3B+4C -> 2C + 3D\n       Reaction(k[13], [A,B], nothing, [3,1], nothing),           # 3A+B -> 0\n       Reaction(k[14], nothing, [A], nothing, [2]),               # 0 -> 2A\n       Reaction(k[15]*A/(2+A), [A], nothing; only_use_rate=true), # A -> 0 with custom rate\n       Reaction(k[16], [A], [B]; only_use_rate=true),             # A -> B with custom rate.\n       Reaction(k[17]*A*exp(B), [C], [D], [2], [1]),              # 2C -> D with non constant rate.\n       Reaction(k[18]*B, nothing, [B], nothing, [2]),             # 0 -> 2B with non constant rate.\n       Reaction(k[19]*t, [A], [B]),                                # A -> B with non constant rate.\n       Reaction(k[20]*t*A, [B,C], [D],[2,1],[2])                  # 2A +B -> 2C with non constant rate.\n  ]\n\nNotes:\n\nnothing can be used to indicate a reaction that has no reactants or no products. In this case the corresponding stoichiometry vector should also be set to nothing.\nThe three-argument form assumes all reactant and product stoichiometric coefficients are one.\n\n\n\n\n\n","category":"type"},{"location":"systems/ReactionSystem/#ModelingToolkit.ReactionSystem","page":"ReactionSystem","title":"ModelingToolkit.ReactionSystem","text":"struct ReactionSystem <: ModelingToolkit.AbstractSystem\n\nA system of chemical reactions.\n\nFields\n\neqs\nThe reactions defining the system.\niv\nIndependent variable (usually time).\nstates\nDependent (state) variables representing amount of each species.\nps\nParameter variables.\nobserved\nname\nThe name of the system\nsystems\nsystems: The internal systems\n\nExample\n\nContinuing from the example in the Reaction definition:\n\nrs = ReactionSystem(rxs, t, [A,B,C,D], k)\n\n\n\n\n\n","category":"type"},{"location":"systems/ReactionSystem/#Composition-and-Accessor-Functions-1","page":"ReactionSystem","title":"Composition and Accessor Functions","text":"","category":"section"},{"location":"systems/ReactionSystem/#","page":"ReactionSystem","title":"ReactionSystem","text":"get_eqs(sys) or equations(sys): The reactions that define the system.\nget_states(sys) or states(sys): The set of chemical species in the system.\nget_ps(sys) or parameters(sys): The parameters of the system.\nindependent_variable(sys): The independent variable of the reaction system, usually time.","category":"page"},{"location":"systems/ReactionSystem/#Query-Functions-1","page":"ReactionSystem","title":"Query Functions","text":"","category":"section"},{"location":"systems/ReactionSystem/#","page":"ReactionSystem","title":"ReactionSystem","text":"oderatelaw\r\njumpratelaw\r\nismassaction","category":"page"},{"location":"systems/ReactionSystem/#ModelingToolkit.oderatelaw","page":"ReactionSystem","title":"ModelingToolkit.oderatelaw","text":"oderatelaw(rx; combinatoric_ratelaw=true)\n\nGiven a Reaction, return the symbolic reaction rate law used in generated ODEs for the reaction. Note, for a reaction defined by\n\nk*X*Y, X+Z --> 2X + Y\n\nthe expression that is returned will be k*X(t)^2*Y(t)*Z(t). For a reaction of the form\n\nk, 2X+3Y --> Z\n\nthe expression that is returned will be k * (X(t)^2/2) * (Y(t)^3/6).\n\nNotes:\n\nAllocates\ncombinatoric_ratelaw=true uses factorial scaling factors in calculating the rate   law, i.e. for 2S -> 0 at rate k the ratelaw would be k*S^2/2!. If   combinatoric_ratelaw=false then the ratelaw is k*S^2, i.e. the scaling factor is   ignored.\n\n\n\n\n\n","category":"function"},{"location":"systems/ReactionSystem/#ModelingToolkit.jumpratelaw","page":"ReactionSystem","title":"ModelingToolkit.jumpratelaw","text":"jumpratelaw(rx; rxvars=get_variables(rx.rate), combinatoric_ratelaw=true)\n\nGiven a Reaction, return the symbolic reaction rate law used in generated stochastic chemical kinetics model SSAs for the reaction. Note, for a reaction defined by\n\nk*X*Y, X+Z --> 2X + Y\n\nthe expression that is returned will be k*X^2*Y*Z. For a reaction of the form\n\nk, 2X+3Y --> Z\n\nthe expression that is returned will be k * binomial(X,2) * binomial(Y,3).\n\nNotes:\n\nrxvars should give the Variables, i.e. species and parameters, the rate depends on.\nAllocates\ncombinatoric_ratelaw=true uses binomials in calculating the rate law, i.e. for 2S -> 0 at rate k the ratelaw would be k*S*(S-1)/2. If combinatoric_ratelaw=false then the ratelaw is k*S*(S-1), i.e. the rate law is not normalized by the scaling factor.\n\n\n\n\n\n","category":"function"},{"location":"systems/ReactionSystem/#ModelingToolkit.ismassaction","page":"ReactionSystem","title":"ModelingToolkit.ismassaction","text":"ismassaction(rx, rs; rxvars = get_variables(rx.rate),\n                              haveivdep = any(var -> isequal(get_iv(rs),var), rxvars),\n                              stateset = Set(states(rs)))\n\nTrue if a given reaction is of mass action form, i.e. rx.rate does not depend on any chemical species that correspond to states of the system, and does not depend explicitly on the independent variable (usually time).\n\nArguments\n\nrx, the Reaction.\nrs, a ReactionSystem containing the reaction.\nOptional: rxvars, Variables which are not in rxvars are ignored as possible dependencies.\nOptional: haveivdep, true if the Reaction rate field explicitly depends on the independent variable.\nOptional: stateset, set of states which if the rxvars are within mean rx is non-mass action.\n\n\n\n\n\n","category":"function"},{"location":"systems/ReactionSystem/#Transformations-1","page":"ReactionSystem","title":"Transformations","text":"","category":"section"},{"location":"systems/ReactionSystem/#","page":"ReactionSystem","title":"ReactionSystem","text":"Base.convert\r\nstructural_simplify","category":"page"},{"location":"systems/ReactionSystem/#Base.convert","page":"ReactionSystem","title":"Base.convert","text":"Base.convert(::Type{<:ODESystem},rs::ReactionSystem)\n\nConvert a ReactionSystem to an ODESystem.\n\nNotes:\n\ncombinatoric_ratelaws=true uses factorial scaling factors in calculating the rate\n\nlaw, i.e. for 2S -> 0 at rate k the ratelaw would be k*S^2/2!. If combinatoric_ratelaws=false then the ratelaw is k*S^2, i.e. the scaling factor is ignored.\n\n\n\n\n\nBase.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)\n\nConvert a ReactionSystem to an NonlinearSystem.\n\nNotes:\n\ncombinatoric_ratelaws=true uses factorial scaling factors in calculating the rate\n\nlaw, i.e. for 2S -> 0 at rate k the ratelaw would be k*S^2/2!. If combinatoric_ratelaws=false then the ratelaw is k*S^2, i.e. the scaling factor is ignored.\n\n\n\n\n\nBase.convert(::Type{<:SDESystem},rs::ReactionSystem)\n\nConvert a ReactionSystem to an SDESystem.\n\nNotes:\n\ncombinatoric_ratelaws=true uses factorial scaling factors in calculating the rate\n\nlaw, i.e. for 2S -> 0 at rate k the ratelaw would be k*S^2/2!. If combinatoric_ratelaws=false then the ratelaw is k*S^2, i.e. the scaling factor is ignored.\n\nnoise_scaling=nothing::Union{Vector{Num},Num,Nothing} allows for linear\n\nscaling of the noise in the chemical Langevin equations. If nothing is given, the default value as in Gillespie 2000 is used. Alternatively, a Num can be given, this is added as a parameter to the system (at the end of the parameter array). All noise terms are linearly scaled with this value. The parameter may be one already declared in the ReactionSystem. Finally, a Vector{Num} can be provided (the length must be equal to the number of reactions). Here the noise for each reaction is scaled by the corresponding parameter in the input vector. This input may contain repeat parameters.\n\n\n\n\n\nBase.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true)\n\nConvert a ReactionSystem to an JumpSystem.\n\nNotes:\n\ncombinatoric_ratelaws=true uses binomials in calculating the rate law, i.e. for 2S -> 0 at rate k the ratelaw would be k*S*(S-1)/2. If combinatoric_ratelaws=false then the ratelaw is k*S*(S-1), i.e. the rate law is not normalized by the scaling factor.\n\n\n\n\n\n","category":"function"},{"location":"systems/ReactionSystem/#Analyses-1","page":"ReactionSystem","title":"Analyses","text":"","category":"section"},{"location":"mtkitize_tutorials/sparse_jacobians/#Automated-Sparse-Analytical-Jacobians-1","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"","category":"section"},{"location":"mtkitize_tutorials/sparse_jacobians/#","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"In many cases where you have large stiff differential equations, getting a sparse Jacobian can be essential for performance. In this tutorial we will show how to use modelingtoolkitize to regenerate an ODEProblem code with the analytical solution to the sparse Jacobian, along with the sparsity pattern required by DifferentialEquations.jl's solvers to specialize the solving process.","category":"page"},{"location":"mtkitize_tutorials/sparse_jacobians/#","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"First let's start out with an implementation of the 2-dimensional Brusselator partial differential equation discretized using finite differences:","category":"page"},{"location":"mtkitize_tutorials/sparse_jacobians/#","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"using DifferentialEquations, ModelingToolkit\n\nconst N = 32\nconst xyd_brusselator = range(0,stop=1,length=N)\nbrusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.\nlimit(a, N) = a == N+1 ? 1 : a == 0 ? N : a\nfunction brusselator_2d_loop(du, u, p, t)\n  A, B, alpha, dx = p\n  alpha = alpha/dx^2\n  @inbounds for I in CartesianIndices((N, N))\n    i, j = Tuple(I)\n    x, y = xyd_brusselator[I[1]], xyd_brusselator[I[2]]\n    ip1, im1, jp1, jm1 = limit(i+1, N), limit(i-1, N), limit(j+1, N), limit(j-1, N)\n    du[i,j,1] = alpha*(u[im1,j,1] + u[ip1,j,1] + u[i,jp1,1] + u[i,jm1,1] - 4u[i,j,1]) +\n                B + u[i,j,1]^2*u[i,j,2] - (A + 1)*u[i,j,1] + brusselator_f(x, y, t)\n    du[i,j,2] = alpha*(u[im1,j,2] + u[ip1,j,2] + u[i,jp1,2] + u[i,jm1,2] - 4u[i,j,2]) +\n                A*u[i,j,1] - u[i,j,1]^2*u[i,j,2]\n    end\nend\np = (3.4, 1., 10., step(xyd_brusselator))\n\nfunction init_brusselator_2d(xyd)\n  N = length(xyd)\n  u = zeros(N, N, 2)\n  for I in CartesianIndices((N, N))\n    x = xyd[I[1]]\n    y = xyd[I[2]]\n    u[I,1] = 22*(y*(1-y))^(3/2)\n    u[I,2] = 27*(x*(1-x))^(3/2)\n  end\n  u\nend\nu0 = init_brusselator_2d(xyd_brusselator)\nprob = ODEProblem(brusselator_2d_loop,u0,(0.,11.5),p)","category":"page"},{"location":"mtkitize_tutorials/sparse_jacobians/#","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"Now let's use modelingtoolkitize to generate the symbolic version:","category":"page"},{"location":"mtkitize_tutorials/sparse_jacobians/#","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"sys = modelingtoolkitize(prob_ode_brusselator_2d)","category":"page"},{"location":"mtkitize_tutorials/sparse_jacobians/#","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"Now we regenerate the problem using jac=true for the analytical Jacobian and sparse=true to make it sparse:","category":"page"},{"location":"mtkitize_tutorials/sparse_jacobians/#","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"sparseprob = ODEProblem(sys,Pair[],(0.,11.5),jac=true,sparse=true)","category":"page"},{"location":"mtkitize_tutorials/sparse_jacobians/#","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"Hard? No! How much did that help?","category":"page"},{"location":"mtkitize_tutorials/sparse_jacobians/#","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"using BenchmarkTools\n@btime solve(prob,save_everystep=false) # 51.714 s (7317 allocations: 70.12 MiB)\n@btime solve(sparseprob,save_everystep=false) # 2.880 s (55533 allocations: 885.09 MiB)","category":"page"},{"location":"mtkitize_tutorials/sparse_jacobians/#","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"Notice though that the analytical solution to the Jacobian can be quite expensive. Thus in some cases we may only want to get the sparsity pattern. In this case, we can simply do:","category":"page"},{"location":"mtkitize_tutorials/sparse_jacobians/#","page":"Automated Sparse Analytical Jacobians","title":"Automated Sparse Analytical Jacobians","text":"sparsepatternprob = ODEProblem(sys,Pair[],(0.,11.5),sparse=true)\n@btime solve(sparsepatternprob,save_everystep=false) # 2.880 s (55533 allocations: 885.09 MiB)","category":"page"},{"location":"systems/NonlinearSystem/#NonlinearSystem-1","page":"NonlinearSystem","title":"NonlinearSystem","text":"","category":"section"},{"location":"systems/NonlinearSystem/#System-Constructors-1","page":"NonlinearSystem","title":"System Constructors","text":"","category":"section"},{"location":"systems/NonlinearSystem/#","page":"NonlinearSystem","title":"NonlinearSystem","text":"NonlinearSystem","category":"page"},{"location":"systems/NonlinearSystem/#ModelingToolkit.NonlinearSystem","page":"NonlinearSystem","title":"ModelingToolkit.NonlinearSystem","text":"struct NonlinearSystem <: ModelingToolkit.AbstractSystem\n\nA nonlinear system of equations.\n\nFields\n\neqs\nVector of equations defining the system.\nstates\nUnknown variables.\nps\nParameters.\nobserved\njac\nJacobian matrix. Note: this field will not be defined until calculate_jacobian is called on the system.\n\nname\nName: the name of the system. These are required to have unique names.\n\nsystems\nsystems: The internal systems\n\ndefaults\ndefaults: The default values to use when initial conditions and/or parameters are not supplied in ODEProblem.\n\nstructure\nstructure: structural information of the system\n\nconnection_type\ntype: type of the system\n\nExamples\n\n@variables x y z\n@parameters σ ρ β\n\neqs = [0 ~ σ*(y-x),\n       0 ~ x*(ρ-z)-y,\n       0 ~ x*y - β*z]\nns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β])\n\n\n\n\n\n","category":"type"},{"location":"systems/NonlinearSystem/#Composition-and-Accessor-Functions-1","page":"NonlinearSystem","title":"Composition and Accessor Functions","text":"","category":"section"},{"location":"systems/NonlinearSystem/#","page":"NonlinearSystem","title":"NonlinearSystem","text":"get_eqs(sys) or equations(sys): The equations that define the nonlinear system.\nget_states(sys) or states(sys): The set of states in the nonlinear system.\nget_ps(sys) or parameters(sys): The parameters of the nonlinear system.","category":"page"},{"location":"systems/NonlinearSystem/#Transformations-1","page":"NonlinearSystem","title":"Transformations","text":"","category":"section"},{"location":"systems/NonlinearSystem/#","page":"NonlinearSystem","title":"NonlinearSystem","text":"structural_simplify\r\nalias_elimination\r\ntearing","category":"page"},{"location":"systems/NonlinearSystem/#ModelingToolkit.StructuralTransformations.tearing","page":"NonlinearSystem","title":"ModelingToolkit.StructuralTransformations.tearing","text":"tearing(sys; simplify=false)\n\nTear the nonlinear equations in system. When simplify=true, we simplify the new residual residual equations after tearing.\n\n\n\n\n\n","category":"function"},{"location":"systems/NonlinearSystem/#Analyses-1","page":"NonlinearSystem","title":"Analyses","text":"","category":"section"},{"location":"systems/NonlinearSystem/#","page":"NonlinearSystem","title":"NonlinearSystem","text":"ModelingToolkit.islinear","category":"page"},{"location":"systems/NonlinearSystem/#Applicable-Calculation-and-Generation-Functions-1","page":"NonlinearSystem","title":"Applicable Calculation and Generation Functions","text":"","category":"section"},{"location":"systems/NonlinearSystem/#","page":"NonlinearSystem","title":"NonlinearSystem","text":"calculate_jacobian\r\ngenerate_jacobian\r\njacobian_sparsity","category":"page"},{"location":"systems/NonlinearSystem/#Problem-Constructors-1","page":"NonlinearSystem","title":"Problem Constructors","text":"","category":"section"},{"location":"systems/NonlinearSystem/#","page":"NonlinearSystem","title":"NonlinearSystem","text":"NonlinearProblem","category":"page"},{"location":"systems/NonlinearSystem/#SciMLBase.NonlinearProblem","page":"NonlinearSystem","title":"SciMLBase.NonlinearProblem","text":"function DiffEqBase.NonlinearProblem{iip}(sys::NonlinearSystem,u0map,\n                                          parammap=DiffEqBase.NullParameters();\n                                          jac = false, sparse=false,\n                                          checkbounds = false,\n                                          linenumbers = true, parallel=SerialForm(),\n                                          kwargs...) where iip\n\nGenerates an NonlinearProblem from a NonlinearSystem and allows for automatically symbolically calculating numerical enhancements.\n\n\n\n\n\n","category":"type"},{"location":"systems/NonlinearSystem/#Torn-Problem-Constructors-1","page":"NonlinearSystem","title":"Torn Problem Constructors","text":"","category":"section"},{"location":"systems/NonlinearSystem/#","page":"NonlinearSystem","title":"NonlinearSystem","text":"BlockNonlinearProblem","category":"page"},{"location":"tutorials/ode_modeling/#Composing-Ordinary-Differential-Equations-1","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"","category":"section"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"This is an introductory example for the usage of ModelingToolkit (MTK). It illustrates the basic user-facing functionality by means of some examples of Ordinary Differential Equations (ODE). Some references to more specific documentation are given at appropriate places.","category":"page"},{"location":"tutorials/ode_modeling/#Copy-Pastable-Simplified-Example-1","page":"Composing Ordinary Differential Equations","title":"Copy-Pastable Simplified Example","text":"","category":"section"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"A much deeper tutorial with forcing functions and sparse Jacobians is all below. But if you want to just see some code and run, here's an example:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"using ModelingToolkit\r\n\r\n@variables t x(t) RHS(t)  # independent and dependent variables\r\n@parameters τ       # parameters\r\nD = Differential(t) # define an operator for the differentiation w.r.t. time\r\n\r\n# your first ODE, consisting of a single equation, indicated by ~\r\n@named fol_separate = ODESystem([ RHS  ~ (1 - x)/τ,\r\n                                  D(x) ~ RHS ])\r\n\r\nusing DifferentialEquations: solve\r\nusing Plots: plot\r\n\r\nprob = ODEProblem(structural_simplify(fol_separate), [x => 0.0], (0.0,10.0), [τ => 3.0])\r\nsol = solve(prob)\r\nplot(sol, vars=[x,RHS])","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"(Image: Simulation result of first-order lag element, with right-hand side)","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Now let's start digging into MTK!","category":"page"},{"location":"tutorials/ode_modeling/#Your-very-first-ODE-1","page":"Composing Ordinary Differential Equations","title":"Your very first ODE","text":"","category":"section"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Let us start with a minimal example. The system to be modelled is a","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"first-order lag element:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"dotx = fracf(t) - x(t)tau","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Here, t is the independent variable (time), x(t) is the (scalar) state variable, f(t) is an external forcing function, and tau is a constant parameter. In MTK, this system can be modelled as follows. For simplicity, we first set the forcing function to a constant value.","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"using ModelingToolkit\r\n\r\n@variables t x(t)  # independent and dependent variables\r\n@parameters τ       # parameters\r\nD = Differential(t) # define an operator for the differentiation w.r.t. time\r\n\r\n# your first ODE, consisting of a single equation, indicated by ~\r\n@named fol_model = ODESystem(D(x) ~ (1 - x)/τ)\r\n      # Model fol_model with 1 equations\r\n      # States (1):\r\n      #   x(t)\r\n      # Parameters (1):\r\n      #   τ","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Note that equations in MTK use the tilde character (~) as equality sign. Also note that the @named macro simply ensures that the symbolic name matches the name in the REPL. If omitted, you can directly set the name keyword.","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"After construction of the ODE, you can solve it using DifferentialEquations.jl:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"using DifferentialEquations: solve\r\nusing Plots: plot\r\n\r\nprob = ODEProblem(fol_model, [x => 0.0], (0.0,10.0), [τ => 3.0])\r\nsolve(prob) |> plot","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"(Image: Simulation result of first-order lag element)","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"The initial state and the parameter values are specified using a mapping from the actual symbolic elements to their values, represented as an array of Pairs, which are constructed using the => operator.","category":"page"},{"location":"tutorials/ode_modeling/#Algebraic-relations-and-structural-simplification-1","page":"Composing Ordinary Differential Equations","title":"Algebraic relations and structural simplification","text":"","category":"section"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"You could separate the calculation of the right-hand side, by introducing an intermediate variable RHS:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"@variables RHS(t)\r\n@named fol_separate = ODESystem([ RHS  ~ (1 - x)/τ,\r\n                                  D(x) ~ RHS ])\r\n      # Model fol_separate with 2 equations\r\n      # States (2):\r\n      #   x(t)\r\n      #   RHS(t)\r\n      # Parameters (1):\r\n      #   τ","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"To directly solve this system, you would have to create a Differential-Algebraic Equation (DAE) problem, since besides the differential equation, there is an additional algebraic equation now. However, this DAE system can obviously be transformed into the single ODE we used in the first example above. MTK achieves this by means of structural simplification:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"fol_simplified = structural_simplify(fol_separate)\r\n\r\nequations(fol_simplified)\r\n      # 1-element Array{Equation,1}:\r\n      #  Differential(t)(x(t)) ~ (τ^-1)*(1 - x(t))\r\n\r\nequations(fol_simplified) == equations(fol_model)\r\n      # true","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"You can extract the equations from a system using equations (and, in the same way, states and parameters). The simplified equation is exactly the same as the original one, so the simulation performance will also be the same. However, there is one difference. MTK does keep track of the eliminated algebraic variables as \"observables\" (see Observables and Variable Elimination). That means, MTK still knows how to calculate them out of the information available in a simulation result. The intermediate variable RHS therefore can be plotted along with the state variable. Note that this has to be requested explicitly, though:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"prob = ODEProblem(fol_simplified, [x => 0.0], (0.0,10.0), [τ => 3.0])\r\nsol = solve(prob)\r\nplot(sol, vars=[x, RHS])","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"(Image: Simulation result of first-order lag element, with right-hand side)","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Note that similarly the indexing of the solution works via the names, and so sol[x] gives the timeseries for x, sol[x,2:10] gives the 2nd through 10th values of x matching sol.t, etc. Note that this works even for variables which have been eliminated, and thus sol[RHS] retrieves the values of RHS.","category":"page"},{"location":"tutorials/ode_modeling/#Specifying-a-time-variable-forcing-function-1","page":"Composing Ordinary Differential Equations","title":"Specifying a time-variable forcing function","text":"","category":"section"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"What if the forcing function (the \"external input\") f(t) is not constant? Obviously, one could use an explicit, symbolic function of time:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"@variables f(t)\r\n@named fol_variable_f = ODESystem([f ~ sin(t), D(x) ~ (f - x)/τ])","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"But often there is time-series data, such as measurement data from an experiment, we want to embed as data in the simulation of a PDE, or as a forcing function on the right-hand side of an ODE – is it is the case here. For this, MTK allows to \"register\" arbitrary Julia functions, which are excluded from symbolic transformations but are just used as-is. So, you could, for example, interpolate a given time series using DataInterpolations.jl. Here, we illustrate this option by a simple lookup (\"zero-order hold\") of a vector of random values:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"value_vector = randn(10)\r\nf_fun(t) = t >= 10 ? value_vector[end] : value_vector[Int(floor(t))+1]\r\n@register f_fun(t)\r\n\r\n@named fol_external_f = ODESystem([f ~ f_fun(t), D(x) ~ (f - x)/τ])\r\nprob = ODEProblem(structural_simplify(fol_external_f), [x => 0.0], (0.0,10.0), [τ => 0.75])\r\n\r\nsol = solve(prob)\r\nplot(sol, vars=[x,f])","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"(Image: Simulation result of first-order lag element, step-wise forcing function)","category":"page"},{"location":"tutorials/ode_modeling/#Building-component-based,-hierarchical-models-1","page":"Composing Ordinary Differential Equations","title":"Building component-based, hierarchical models","text":"","category":"section"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Working with simple one-equation systems is already fun, but composing more complex systems from simple ones is even more fun. Best practice for such a \"modeling framework\" could be to use factory functions for model components:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"function fol_factory(separate=false;name)\r\n    @parameters τ\r\n    @variables t x(t) f(t) RHS(t)\r\n\r\n    eqs = separate ? [RHS ~ (f - x)/τ,\r\n                      D(x) ~ RHS] :\r\n                      D(x) ~(f - x)/τ\r\n\r\n    ODESystem(eqs;name)\r\nend","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Such a factory can then used to instantiate the same component multiple times, but allows for customization:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"@named fol_1 = fol_factory()\r\n@named fol_2 = fol_factory(true) # has observable RHS","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Now, these two components can be used as subsystems of a parent system, i.e. one level higher in the model hierarchy. The connections between the components again are just algebraic relations:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"connections = [ fol_1.f ~ 1.5,\r\n                fol_2.f ~ fol_1.x ]\r\n\r\n@named connected = ODESystem(connections; systems=[fol_1,fol_2])\r\n      # Model connected with 5 equations\r\n      # States (5):\r\n      #   fol_1₊f(t)\r\n      #   fol_2₊f(t)\r\n      #   fol_1₊x(t)\r\n      #   fol_2₊x(t)\r\n      # ⋮\r\n      # Parameters (2):\r\n      #   fol_1₊τ\r\n      #   fol_2₊τ","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"All equations, variables and parameters are collected, but the structure of the hierarchical model is still preserved. That is, you can still get information about fol_1 by addressing it by connected.fol_1, or its parameter by connected.fol_1.τ. Before simulation, we again eliminate the algebraic variables and connection equations from the system using structural simplification:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"connected_simp = structural_simplify(connected)\r\n      # Model connected with 2 equations\r\n      # States (2):\r\n      #   fol_1₊x(t)\r\n      #   fol_2₊x(t)\r\n      # Parameters (2):\r\n      #   fol_1₊τ\r\n      #   fol_2₊τ\r\n      # Incidence matrix:\r\n      #   [1, 1]  =  ×\r\n      #   [2, 1]  =  ×\r\n      #   [2, 2]  =  ×\r\n      #   [1, 3]  =  ×\r\n      #   [2, 4]  =  ×\r\n\r\nequations(connected_simp)\r\n      # 2-element Array{Equation,1}:\r\n      #  Differential(t)(fol_1₊x(t)) ~ (fol_1₊τ^-1)*(1.5 - fol_1₊x(t))\r\n      #  Differential(t)(fol_2₊x(t)) ~ (fol_2₊τ^-1)*(fol_1₊x(t) - fol_2₊x(t))","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"As expected, only the two state-derivative equations remain, as if you had manually eliminated as many variables as possible from the equations. As mentioned above, the hierarchical structure is preserved though. So the initial state and the parameter values can be specified accordingly when building the ODEProblem:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"u0 = [ fol_1.x => -0.5,\r\n       fol_2.x => 1.0 ]\r\n\r\np = [ fol_1.τ => 2.0,\r\n      fol_2.τ => 4.0 ]\r\n\r\nprob = ODEProblem(connected_simp, u0, (0.0,10.0), p)\r\nsolve(prob) |> plot","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"(Image: Simulation of connected system (two first-order lag elements in series))","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"More on this topic may be found in Composing Models and Building Reusable Components.","category":"page"},{"location":"tutorials/ode_modeling/#Defaults-1","page":"Composing Ordinary Differential Equations","title":"Defaults","text":"","category":"section"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Often it is a good idea to specify reasonable values for the initial state and the parameters of a model component. Then, these do not have to be explicitly specified when constructing the ODEProblem.","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"function unitstep_fol_factory(;name)\r\n    @parameters τ\r\n    @variables t x(t)\r\n    ODESystem(D(x) ~ (1 - x)/τ; name, defaults=Dict(x=>0.0, τ=>1.0))\r\nend\r\n\r\nODEProblem(unitstep_fol_factory(name=:fol),[],(0.0,5.0),[]) |> solve","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Note that the defaults can be functions of the other variables, which is then resolved at the time of the problem construction. Of course, the factory function could accept additional arguments to optionally specify the initial state or parameter values, etc.","category":"page"},{"location":"tutorials/ode_modeling/#Symbolic-and-sparse-derivatives-1","page":"Composing Ordinary Differential Equations","title":"Symbolic and sparse derivatives","text":"","category":"section"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"One advantage of a symbolic toolkit is that derivatives can be calculated explicitly, and that the incidence matrix of partial derivatives (the \"sparsity pattern\") can also be explicitly derived. These two facts lead to a substantial speedup of all model calculations, e.g. when simulating a model over time using an ODE solver.","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"By default, analytical derivatives and sparse matrices, e.g. for the Jacobian, the matrix of first partial derivatives, are not used. Let's benchmark this (prob still is the problem using the connected_simp system above):","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"using BenchmarkTools\r\nusing DifferentialEquations: Rodas4\r\n\r\n@btime solve(prob, Rodas4());\r\n      # 251.300 μs (873 allocations: 31.18 KiB)","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Now have MTK provide sparse, analytical derivatives to the solver. This has to be specified during the construction of the ODEProblem:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"prob_an = ODEProblem(connected_simp, u0, (0.0,10.0), p; jac=true, sparse=true)\r\n\r\n@btime solve(prob_an, Rodas4());\r\n      # 142.899 μs (1297 allocations: 83.96 KiB)","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"The speedup is significant. For this small dense model (3 of 4 entries are populated), using sparse matrices is counterproductive in terms of required memory allocations. For large, hierarchically built models, which tend to be sparse, speedup and the reduction of memory allocation can be expected to be substantial. In addition, these problem builders allow for automatic parallelism using the structural information. For more information, see the ODESystem page.","category":"page"},{"location":"tutorials/ode_modeling/#Notes-and-pointers-how-to-go-on-1","page":"Composing Ordinary Differential Equations","title":"Notes and pointers how to go on","text":"","category":"section"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Here are some notes that may be helpful during your initial steps with MTK:","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Sometimes, the symbolic engine within MTK is not able to correctly identify the independent variable (e.g. time) out of all variables. In such a case, you usually get an error that some variable(s) is \"missing from variable map\". In most cases, it is then sufficient to specify the independent variable as second argument to ODESystem, e.g. ODESystem(eqs, t).\nA completely macro-free usage of MTK is possible and is discussed in a separate tutorial. This is for package developers, since the macros are only essential for automatic symbolic naming for modelers.\nVector-valued parameters and variables are possible. A cleaner, more consistent treatment of these is work in progress, though. Once finished, this introductory tutorial will also cover this feature.","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Where to go next?","category":"page"},{"location":"tutorials/ode_modeling/#","page":"Composing Ordinary Differential Equations","title":"Composing Ordinary Differential Equations","text":"Not sure how MTK relates to similar tools and packages? Read Comparison of ModelingToolkit vs Equation-Based Modeling Languages.\nDepending on what you want to do with MTK, have a look at some of the other Symbolic Modeling Tutorials.\nIf you want to automatically convert an existing function to a symbolic representation, you might go through the ModelingToolkitize Tutorials.\nTo learn more about the inner workings of MTK, consider the sections under Basics and System Types.","category":"page"},{"location":"comparison/#Comparison-of-ModelingToolkit-vs-Equation-Based-Modeling-Languages-1","page":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","title":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","text":"","category":"section"},{"location":"comparison/#Comparison-Against-Modelica-1","page":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","title":"Comparison Against Modelica","text":"","category":"section"},{"location":"comparison/#","page":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","title":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","text":"Both Modelica and ModelingToolkit.jl are acausal modeling languages.\nModelica is a language with many different implementations, such as Dymola and OpenModelica, which have differing levels of performance and can give different results on the same model. Many of the commonly used Modelica compilers are not open source. ModelingToolkit.jl is a language with a single canonical open source implementation.\nAll 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 which is roughly equivalent to the dae_index_lowering and structural_simplify 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.\nModelica is a declarative programming language. ModelingToolkit.jl is a declarative symbolic modeling language used from within the Julia programming language. Its programming language semantics, such as loop constructs and conditionals, can be used to more easily generate models.\nModelica is an object-oriented single dispatch language. ModelingToolkit.jl, built on Julia, uses multiple dispatch extensively to simplify code.\nMany Modelica compilers supply a GUI. ModelingToolkit.jl does not.\nModelica can be used to simulate ODE and DAE systems. ModelingToolkit.jl has a much more expansive set of system types, including nonlinear systems, SDEs, PDEs, and more.","category":"page"},{"location":"comparison/#Comparison-Against-Simulink-1","page":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","title":"Comparison Against Simulink","text":"","category":"section"},{"location":"comparison/#","page":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","title":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","text":"Simulink is a causal modeling environment, whereas ModelingToolkit.jl is an acausal modeling environment. For an overview of the differences, consult academic reviews such as this one. In this sense, ModelingToolkit.jl is more similar to the Simscape sub-environment.\nSimulink is used from MATLAB while ModelingToolkit.jl is used from Julia. Thus any user defined functions have the performance of their host language. For information on the performance differences between Julia and MATLAB, consult open benchmarks which demonstrate Julia as an order of magnitude or more faster in many cases due to its JIT compilation.\nSimulink uses the MATLAB differential equation solvers while ModelingToolkit.jl uses DifferentialEquations.jl. For a systematic comparison between the solvers, consult open benchmarks which demonstrate two orders of magnitude performance advantage for the native Julia solvers across many benchmark problems.\nSimulink comes with a Graphical User Interface (GUI), ModelingToolkit.jl does not.\nSimulink is a proprietary software, meaning users cannot actively modify or extend the software. ModelingToolkit.jl is built in Julia and used in Julia, where users can actively extend and modify the software interactively in the REPL and contribute to its open source repositories.\nSimulink covers ODE and DAE systems. ModelingToolkit.jl has a much more expansive set of system types, including SDEs, PDEs, optimization problems, and more.","category":"page"},{"location":"comparison/#Comparison-Against-CASADI-1","page":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","title":"Comparison Against CASADI","text":"","category":"section"},{"location":"comparison/#","page":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","title":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","text":"CASADI is written in C++ but used from Python/MATLAB, meaning that it cannot be directly extended by users unless they are using the C++ interface and run a local build of CASADI. ModelingToolkit.jl is both written and used from Julia, meaning that users can easily extend the library on the fly, even interactively in the REPL.\nCASADI includes limited support for Computer Algebra System (CAS) functionality, while ModelingToolkit.jl is built on the full Symbolics.jl CAS.\nCASADI supports DAE and ODE problems via SUNDIALS IDAS and CVODES. ModelingToolkit.jl supports DAE and ODE problems via DifferentialEquations.jl, of which Sundials.jl is <1% of the total available solvers and is outperformed by the native Julia solvers on the vast majority of the benchmark equations. In addition, the DifferentialEquations.jl interface is confederated, meaning that any user can dynamically extend the system to add new solvers to the interface by defining new dispatches of solve.\nCASADI's DAEBuilder does not implement efficiency transformations like tearing which are standard in the ModelingToolkit.jl transformation pipeline.\nCASADI supports special functionality for quadratic programming problems while ModelingToolkit only provides nonlinear programming via OptimizationSystem.\nModelingToolkit.jl integrates with its host language Julia, so Julia code can be automatically converted into ModelingToolkit expressions. Users of CASADI must explicitly create CASADI expressions.","category":"page"},{"location":"comparison/#Comparison-Against-Modia.jl-1","page":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","title":"Comparison Against Modia.jl","text":"","category":"section"},{"location":"comparison/#","page":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","title":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","text":"Modia.jl is a Modelica-like system built in pure Julia. As such, its syntax is a domain-specific language (DSL) specified by macros to mirror the Modelica syntax.\nModia's compilation pipeline is similar to the Dymola symbolic processing pipeline 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 and tearing algorithms.\nModia supports DAE problems via SUNDIALS IDAS. ModelingToolkit.jl supports DAE and ODE problems via DifferentialEquations.jl, of which Sundials.jl is <1% of the total available solvers and is outperformed by the native Julia solvers on the vast majority of the benchmark equations. In addition, the DifferentialEquations.jl interface is confederated, meaning that any user can dynamically extend the system to add new solvers to the interface by defining new dispatches of solve.\nModelingToolkit.jl integrates with its host language Julia, so Julia code can be automatically converted into ModelingToolkit expressions. Users of Modia must explicitly create Modia expressions within its macro.\nModia covers DAE systems. ModelingToolkit.jl has a much more expansive set of system types, including SDEs, PDEs, optimization problems, and more.","category":"page"},{"location":"comparison/#Comparison-Against-Causal.jl-1","page":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","title":"Comparison Against Causal.jl","text":"","category":"section"},{"location":"comparison/#","page":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","title":"Comparison of ModelingToolkit vs Equation-Based Modeling Languages","text":"Causal.jl is a causal modeling environment, whereas ModelingToolkit.jl is an acausal modeling environment. For an overview of the differences, consult academic reviews such as this one.\nBoth ModelingToolkit.jl and Causal.jl use DifferentialEquations.jl as the backend solver library.\nCausal.jl lets one add arbitrary equation systems to a given node, and allow the output to effect the next node. This means an SDE may drive an ODE. These two portions are solved with different solver methods in tandem. In ModelingToolkit.jl, such connections promote the whole system to an SDE. This results in better accuracy and stability, though in some cases it can be less performant.\nCausal.jl, similar to Simulink, breaks algebraic loops via inexact heuristics. ModelingToolkit.jl treats algebraic loops exactly through algebraic equations in the generated model.","category":"page"},{"location":"systems/ControlSystem/#ControlSystem-1","page":"ControlSystem","title":"ControlSystem","text":"","category":"section"},{"location":"systems/ControlSystem/#System-Constructors-1","page":"ControlSystem","title":"System Constructors","text":"","category":"section"},{"location":"systems/ControlSystem/#","page":"ControlSystem","title":"ControlSystem","text":"ControlSystem","category":"page"},{"location":"systems/ControlSystem/#ModelingToolkit.ControlSystem","page":"ControlSystem","title":"ModelingToolkit.ControlSystem","text":"struct ControlSystem <: ModelingToolkit.AbstractControlSystem\n\nA system describing an optimal control problem. This contains a loss function and ordinary differential equations with control variables that describe the dynamics.\n\nFields\n\nloss\nThe Loss function\neqs\nThe ODEs defining the system.\niv\nIndependent variable.\nstates\nDependent (state) variables.\ncontrols\nControl variables.\nps\nParameter variables.\nobserved\nname\nName: the name of the system. These are required to have unique names.\n\nsystems\nsystems: The internal systems\n\ndefaults\ndefaults: The default values to use when initial conditions and/or parameters are not supplied in ODEProblem.\n\nExample\n\nusing ModelingToolkit\n\n@variables t x(t) v(t) u(t)\nD = Differential(t)\n\nloss = (4-x)^2 + 2v^2 + u^2\neqs = [\n    D(x) ~ v\n    D(v) ~ u^3\n]\n\nsys = ControlSystem(loss,eqs,t,[x,v],[u],[])\n\n\n\n\n\n","category":"type"},{"location":"systems/ControlSystem/#Composition-and-Accessor-Functions-1","page":"ControlSystem","title":"Composition and Accessor Functions","text":"","category":"section"},{"location":"systems/ControlSystem/#","page":"ControlSystem","title":"ControlSystem","text":"get_eqs(sys) or equations(sys): The equations that define the system.\nget_states(sys) or states(sys): The set of states in the system.\nget_ps(sys) or parameters(sys): The parameters of the system.\nget_controls(sys) or controls(sys): The control variables of the system","category":"page"},{"location":"systems/ControlSystem/#Transformations-1","page":"ControlSystem","title":"Transformations","text":"","category":"section"},{"location":"systems/ControlSystem/#","page":"ControlSystem","title":"ControlSystem","text":"ModelingToolkit.runge_kutta_discretize\nstructural_simplify","category":"page"},{"location":"systems/ControlSystem/#ModelingToolkit.runge_kutta_discretize","page":"ControlSystem","title":"ModelingToolkit.runge_kutta_discretize","text":"runge_kutta_discretize(sys::ControlSystem,dt,tspan;\n                       tab = ModelingToolkit.constructRadauIIA5())\n\nTransforms a nonlinear optimal control problem into a constrained OptimizationProblem according to a Runge-Kutta tableau that describes a collocation method. Requires a fixed dt over a given timespan. Defaults to using the 5th order RadauIIA tableau, and altnerative tableaus can be specified using the SciML tableau style.\n\n\n\n\n\n","category":"function"},{"location":"systems/ControlSystem/#ModelingToolkit.structural_simplify","page":"ControlSystem","title":"ModelingToolkit.structural_simplify","text":"structural_simplify(sys)\n\n\nStructurally simplify algebraic equations in a system and compute the topological sort of the observed equations.\n\n\n\n\n\n","category":"function"},{"location":"systems/ControlSystem/#Analyses-1","page":"ControlSystem","title":"Analyses","text":"","category":"section"},{"location":"tutorials/optimization/#Modeling-Optimization-Problems-1","page":"Modeling Optimization Problems","title":"Modeling Optimization Problems","text":"","category":"section"},{"location":"tutorials/optimization/#","page":"Modeling Optimization Problems","title":"Modeling Optimization Problems","text":"using ModelingToolkit, GalacticOptim\n\n@variables x y\n@parameters a b\nloss = (a - x)^2 + b * (y - x^2)^2\nsys = OptimizationSystem(loss,[x,y],[a,b])\n\nu0 = [\n    x=>1.0\n    y=>2.0\n]\np = [\n    a => 6.0\n    b => 7.0\n]\n\nprob = OptimizationProblem(sys,u0,p,grad=true,hess=true)\nsolve(prob,Newton())","category":"page"},{"location":"tutorials/optimization/#","page":"Modeling Optimization Problems","title":"Modeling Optimization Problems","text":"Needs more text but it's super cool and auto-parallelizes and sparsifies too. Plus you can hierarchically nest systems to have it generate huge optimization problems.","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#Nonlinear-Optimal-Control-1","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"","category":"section"},{"location":"tutorials/nonlinear_optimal_control/#Note:-this-is-still-a-work-in-progress!-1","page":"Nonlinear Optimal Control","title":"Note: this is still a work in progress!","text":"","category":"section"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"The ControlSystem type is an interesting system because, unlike other system types, it cannot be numerically solved on its own. Instead, it must be transformed into another system before solving. Standard methods such as the \"direct method\", \"multiple shooting\", or \"discretize-then-optimize\" can all be phrased as symbolic transformations to a ControlSystem: this is the strategy of this methodology.","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#Defining-a-Nonlinear-Optimal-Control-Problem-1","page":"Nonlinear Optimal Control","title":"Defining a Nonlinear Optimal Control Problem","text":"","category":"section"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"Here we will start by defining a classic optimal control problem. Let:","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"x^ = u^3(t)","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"where we want to optimize our controller u(t) such that the following is minimized:","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"L(theta) = sum_i Vert 4 - x(t_i) Vert + 2 Vert x^prime(t_i) Vert + Vert u(t_i) Vert","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"where i is measured on (0,8) at 0.01 intervals. To do this, we rewrite the ODE in first order form:","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"beginaligned\nx^prime = v \nv^ = u^3(t) \nendaligned","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"and thus","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"L(theta) = sum_i Vert 4 - x(t_i) Vert + 2 Vert v(t_i) Vert + Vert u(t_i) Vert","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"is our loss function on the first order system.","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"Defining such a control system is similar to an ODESystem, except we must also specify a control variable u(t) and a loss function. Together, this problem looks as follows:","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"using ModelingToolkit\n\n@variables t x(t) v(t) u(t)\n@parameters p[1:2]\nD = Differential(t)\n\nloss = (4-x)^2 + 2v^2 + u^2\neqs = [\n    D(x) ~ v - p[2]*x\n    D(v) ~ p[1]*u^3 + v\n]\n\nsys = ControlSystem(loss,eqs,t,[x,v],[u],p)","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#Solving-a-Control-Problem-via-Discretize-Then-Optimize-1","page":"Nonlinear Optimal Control","title":"Solving a Control Problem via Discretize-Then-Optimize","text":"","category":"section"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"One common way to solve nonlinear optimal control problems is by transforming them into an optimization problem by performing a Runge-Kutta discretization of the differential equation system and imposing equalities between variables in the same steps. This can be done via the runge_kutta_discretize transformation on the ControlSystem. While a tableau tab can be specified, it defaults to a 5th order RadauIIA collocation, which is a common method in the field. To perform this discretization, we simply need to give a dt and a timespan on which to discretize:","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"dt = 0.1\ntspan = (0.0,1.0)\nsys = runge_kutta_discretize(sys,dt,tspan)","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"Now sys is an OptimizationSystem which, when solved, gives the values of x(t), v(t), and u(t). Thus we solve the OptimizationSystem using GalacticOptim.jl:","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"u0 = rand(length(states(sys))) # guess for the state values\nprob = OptimizationProblem(sys,u0,[0.1,0.1],grad=true)\n\nusing GalacticOptim, Optim\nsol = solve(prob,BFGS())","category":"page"},{"location":"tutorials/nonlinear_optimal_control/#","page":"Nonlinear Optimal Control","title":"Nonlinear Optimal Control","text":"And this is missing some nice interfaces and ignores the equality constraints right now so the tutorial is not complete.","category":"page"},{"location":"basics/FAQ/#Frequently-Asked-Questions-1","page":"Frequently Asked Questions","title":"Frequently Asked Questions","text":"","category":"section"},{"location":"basics/FAQ/#Getting-the-index-for-a-symbol-1","page":"Frequently Asked Questions","title":"Getting the index for a symbol","text":"","category":"section"},{"location":"basics/FAQ/#","page":"Frequently Asked Questions","title":"Frequently Asked Questions","text":"Since ordering of symbols is not guaranteed after symbolic transformations, one should normally refer to values by their name. For example, sol[lorenz.x] from the solution. But what if you need to get the index? The following helper function will do the trick:","category":"page"},{"location":"basics/FAQ/#","page":"Frequently Asked Questions","title":"Frequently Asked Questions","text":"indexof(sym,syms) = findfirst(isequal(sym),syms)\nindexof(σ,parameters(sys))","category":"page"},{"location":"basics/FAQ/#Transforming-value-maps-to-arrays-1","page":"Frequently Asked Questions","title":"Transforming value maps to arrays","text":"","category":"section"},{"location":"basics/FAQ/#","page":"Frequently Asked Questions","title":"Frequently Asked Questions","text":"ModelingToolkit.jl allows (and recommends) input maps like [x => 2.0, y => 3.0] because symbol ordering is not guaranteed. However, what if you want to get the lowered array? You can use the internal function varmap_to_vars. For example:","category":"page"},{"location":"basics/FAQ/#","page":"Frequently Asked Questions","title":"Frequently Asked Questions","text":"pnew = varmap_to_vars([β=>3.0, c=>10.0, γ=>2.0],parameters(sys))","category":"page"},{"location":"tutorials/acausal_components/#Acausal-Component-Based-Modeling-the-RC-Circuit-1","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"","category":"section"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"In this tutorial we will build a hierarchical acausal component-based model of the RC circuit. The RC circuit is a simple example where we connect a resistor and a capacitor. Kirchoff's laws are then applied to state equalities between currents and voltages. This specifies a differential-algebraic equation (DAE) system, where the algebraic equations are given by the constraints and equalities between different component variables. We then simplify this to an ODE by eliminating the equalities before solving. Let's see this in action.","category":"page"},{"location":"tutorials/acausal_components/#Copy-Paste-Example-1","page":"Acausal Component-Based Modeling the RC Circuit","title":"Copy-Paste Example","text":"","category":"section"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"using ModelingToolkit, Plots, DifferentialEquations\n\n@parameters t\n\n# Basic electric components\nfunction Pin(;name)\n    @variables v(t) i(t)\n    ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])\nend\n\nfunction Ground(;name)\n    @named g = Pin()\n    eqs = [g.v ~ 0]\n    ODESystem(eqs, t, [], [], systems=[g], name=name)\nend\n\nfunction Resistor(;name, R = 1.0)\n    val = R\n    @named p = Pin()\n    @named n = Pin()\n    @variables v(t)\n    @parameters R\n    eqs = [\n           v ~ p.v - n.v\n           0 ~ p.i + n.i\n           v ~ p.i * R\n          ]\n    ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)\nend\n\nfunction Capacitor(; name, C = 1.0)\n    val = C\n    @named p = Pin()\n    @named n = Pin()\n    @variables v(t)\n    @parameters C\n    D = Differential(t)\n    eqs = [\n           v ~ p.v - n.v\n           0 ~ p.i + n.i\n           D(v) ~ p.i / C\n          ]\n    ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)\nend\n\nfunction ConstantVoltage(;name, V = 1.0)\n    val = V\n    @named p = Pin()\n    @named n = Pin()\n    @parameters V\n    eqs = [\n           V ~ p.v - n.v\n           0 ~ p.i + n.i\n          ]\n    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)\nend\n\nR = 1.0\nC = 1.0\nV = 1.0\n@named resistor = Resistor(R=R)\n@named capacitor = Capacitor(C=C)\n@named source = ConstantVoltage(V=V)\n@named ground = Ground()\n\nfunction connect_pins(ps...)\n    eqs = [\n           0 ~ sum(p->p.i, ps) # KCL\n          ]\n    # KVL\n    for i in 1:length(ps)-1\n        push!(eqs, ps[i].v ~ ps[i+1].v)\n    end\n\n    return eqs\nend\n\nrc_eqs = [\n          connect_pins(source.p, resistor.p)\n          connect_pins(resistor.n, capacitor.p)\n          connect_pins(capacitor.n, source.n, ground.g)\n         ]\n\n@named rc_model = ODESystem(rc_eqs, t,\n                            systems=[resistor, capacitor, source, ground])\nsys = structural_simplify(rc_model)\nu0 = [\n      capacitor.v => 0.0\n      capacitor.p.i => 0.0\n     ]\nprob = ODAEProblem(sys, u0, (0, 10.0))\nsol = solve(prob, Tsit5())\nplot(sol)","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"(Image: )","category":"page"},{"location":"tutorials/acausal_components/#Explanation-1","page":"Acausal Component-Based Modeling the RC Circuit","title":"Explanation","text":"","category":"section"},{"location":"tutorials/acausal_components/#Building-the-Component-Library-1","page":"Acausal Component-Based Modeling the RC Circuit","title":"Building the Component Library","text":"","category":"section"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"For each of our components we use a Julia function which emits an ODESystem. At the top we start with defining the fundamental qualities of an electrical circuit component. At every input and output pin a circuit component has two values: the current at the pin and the voltage. Thus we define the Pin component to simply be the values there:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"function Pin(;name)\n    @variables v(t) i(t)\n    ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])\nend","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Note that this is an incompletely specified ODESystem: it cannot be simulated on its own because the equations for v(t) and i(t) are unknown. Instead this just gives a common syntax for receiving this pair with some default values. Notice that in a component we define the name as a keyword argument: this is because later we will generate different Pin objects with different names to correspond to duplicates of this topology with unique variables. One can then construct a Pin like:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Pin(name=:mypin1)","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"or equivalently using the @named helper macro:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"@named mypin1 = Pin()","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Next we build our ground node. A ground node is just a pin that is connected to a constant voltage reservoir, typically taken to be V=0. Thus to define this component, we generate an ODESystem with a Pin subcomponent and specify that the voltage in such a Pin is equal to zero. This gives:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"function Ground(;name)\n    @named g = Pin()\n    eqs = [g.v ~ 0]\n    ODESystem(eqs, t, [], [], systems=[g], name=name)\nend","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Next we build a resistor. A resistor is an object that has two Pins, the positive and the negative pins, and follows Ohm's law: v = i*r. The voltage of the resistor is given as the voltage difference across the two pins while by conservation of charge we know that the current in must equal the current out, which means (no matter the direction of the current flow) the sum of the currents must be zero. This leads to our resistor equations:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"function Resistor(;name, R = 1.0)\n    val = R\n    @named p = Pin()\n    @named n = Pin()\n    @variables v(t)\n    @parameters R\n    eqs = [\n           v ~ p.v - n.v\n           0 ~ p.i + n.i\n           v ~ p.i * R\n          ]\n    ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)\nend","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Notice that we have created this system with a defaults for the resistor's resistance. By doing so, if the resistance of this resistor is not overridden by a higher level default or overridden at ODEProblem construction time, this will be the value of the resistance.","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Using our knowledge of circuits we similarly construct the Capacitor:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"function Capacitor(; name, C = 1.0)\n    val = C\n    @named p = Pin()\n    @named n = Pin()\n    @variables v(t)\n    @parameters C\n    D = Differential(t)\n    eqs = [\n           v ~ p.v - n.v\n           0 ~ p.i + n.i\n           D(v) ~ p.i / C\n          ]\n    ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)\nend","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Now we want to build a constant voltage electrical source term. We can think of this as similarly being a two pin object, where the object itself is kept at a constant voltage, essentially generating the electrical current. We would then model this as:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"function ConstantVoltage(;name, V = 1.0)\n    val = V\n    @named p = Pin()\n    @named n = Pin()\n    @parameters V\n    eqs = [\n           V ~ p.v - n.v\n           0 ~ p.i + n.i\n          ]\n    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)\nend","category":"page"},{"location":"tutorials/acausal_components/#Connecting-and-Simulating-Our-Electrical-Circuit-1","page":"Acausal Component-Based Modeling the RC Circuit","title":"Connecting and Simulating Our Electrical Circuit","text":"","category":"section"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Now we are ready to simulate our circuit. Let's build our four components: a resistor, capacitor, source, and ground term. For simplicity we will make all of our parameter values 1. This is done by:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"R = 1.0\nC = 1.0\nV = 1.0\n@named resistor = Resistor(R=R)\n@named capacitor = Capacitor(C=C)\n@named source = ConstantVoltage(V=V)\n@named ground = Ground()","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Next we have to define how we connect the circuit. Whenever two Pins in a circuit are connected together, the system satisfies Kirchoff's laws, i.e. that currents sum to zero and voltages across the pins are equal. Thus we will build a helper function connect_pins which implements these rules:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"function connect_pins(ps...)\n    eqs = [\n           0 ~ sum(p->p.i, ps) # KCL\n          ]\n    # KVL\n    for i in 1:length(ps)-1\n        push!(eqs, ps[i].v ~ ps[i+1].v)\n    end\n\n    return eqs\nend","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Finally we will connect the pieces of our circuit together. Let's connect the positive pin of the resistor to the source, the negative pin of the resistor to the capacitor, and the negative pin of the capacitor to a junction between the source and the ground. This would mean our connection equations are:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"rc_eqs = [\n          connect_pins(source.p, resistor.p)\n          connect_pins(resistor.n, capacitor.p)\n          connect_pins(capacitor.n, source.n, ground.g)\n         ]","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Finally we build our four component model with these connection rules:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"@named rc_model = ODESystem(rc_eqs, t,\n                            systems=[resistor, capacitor, source, ground])","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Notice that this model is acasual because we have not specified anything about the causality of the model. We have simply specified what is true about each of the variables. This forms a system of differential-algebraic equations (DAEs) which define the evolution of each state of the system. The equations are:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"equations(rc_model)\n\n16-element Vector{Equation}:\n 0 ~ resistor₊p₊i(t) + source₊p₊i(t)\n source₊p₊v(t) ~ resistor₊p₊v(t)\n 0 ~ capacitor₊p₊i(t) + resistor₊n₊i(t)\n resistor₊n₊v(t) ~ capacitor₊p₊v(t)\n ⋮\n Differential(t)(capacitor₊v(t)) ~ capacitor₊p₊i(t)*(capacitor₊C^-1)\n source₊V ~ source₊p₊v(t) - (source₊n₊v(t))\n 0 ~ source₊n₊i(t) + source₊p₊i(t)\n ground₊g₊v(t) ~ 0","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"the states are:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"states(rc_model)\n\n16-element Vector{Term{Real}}:\n resistor₊p₊i(t)\n source₊p₊i(t)\n source₊p₊v(t)\n resistor₊p₊v(t)\n ⋮\n source₊n₊v(t)\n ground₊g₊v(t)\n resistor₊v(t)\n capacitor₊v(t)","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"and the parameters are:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"parameters(rc_model)\n\n3-element Vector{Any}:\n resistor₊R\n capacitor₊C\n source₊V","category":"page"},{"location":"tutorials/acausal_components/#Simplifying-and-Solving-this-System-1","page":"Acausal Component-Based Modeling the RC Circuit","title":"Simplifying and Solving this System","text":"","category":"section"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"This system could be solved directly as a DAE using one of the DAE solvers from DifferentialEquations.jl. However, let's take a second to symbolically simplify the system before doing the solve. The function structural_simplify looks for all of the equalities and eliminates unnecessary variables to build the leanest numerical representation of the system. Let's see what it does here:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"sys = structural_simplify(rc_model)\nequations(sys)\n\n2-element Vector{Equation}:\n 0 ~ capacitor₊v(t) + resistor₊R*capacitor₊p₊i(t) - source₊V\n Differential(t)(capacitor₊v(t)) ~ capacitor₊p₊i(t)*(capacitor₊C^-1)","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"states(sys)\n\n2-element Vector{Any}:\n capacitor₊v(t)\n capacitor₊p₊i(t)","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"After structural simplification we are left with a system of only two equations with two state variables. One of the equations is a differential equation while the other is an algebraic equation. We can then give the values for the initial conditions of our states and solve the system by converting it to an ODEProblem in mass matrix form and solving it with an ODEProblem mass matrix DAE solver. This is done as follows:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"u0 = [\n      capacitor.v => 0.0\n      capacitor.p.i => 0.0\n     ]\nprob = ODEProblem(sys, u0, (0, 10.0))\nsol = solve(prob, Rodas4())\nplot(sol)","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"(Image: )","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"However, we can also choose to use the \"torn nonlinear system\" to remove all of the algebraic variables from the solution of the system. Note that this requires having done structural_simplify. This is done by using ODAEProblem like:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"u0 = [\n      capacitor.v => 0.0\n      capacitor.p.i => 0.0\n     ]\nprob = ODAEProblem(sys, u0, (0, 10.0))\nsol = solve(prob, Rodas4())\nplot(sol)","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"(Image: )","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"Notice that this solves the whole system by only solving for one variable!","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"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 state variables into observed variables. Let's see what our observed variables are:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"observed(sys)\n\n14-element Vector{Equation}:\n resistor₊p₊i(t) ~ capacitor₊p₊i(t)\n capacitor₊n₊v(t) ~ 0.0\n source₊n₊v(t) ~ 0.0\n ground₊g₊i(t) ~ 0.0\n source₊n₊i(t) ~ capacitor₊p₊i(t)\n source₊p₊i(t) ~ -capacitor₊p₊i(t)\n capacitor₊n₊i(t) ~ -capacitor₊p₊i(t)\n resistor₊n₊i(t) ~ -capacitor₊p₊i(t)\n ground₊g₊v(t) ~ 0.0\n source₊p₊v(t) ~ source₊V\n capacitor₊p₊v(t) ~ capacitor₊v(t)\n resistor₊p₊v(t) ~ source₊p₊v(t)\n resistor₊n₊v(t) ~ capacitor₊p₊v(t)\n resistor₊v(t) ~ -((capacitor₊p₊v(t)) - (source₊p₊v(t)))","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"These are explicit algebraic equations which can then be used to reconstruct the required variables on the fly. This leads to dramatic computational savings because implicitly solving an ODE scales like O(n^3), so making there be as few states as possible is good!","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"The solution object can be accessed via its symbols. For example, let's retrieve the voltage of the resistor over time:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"sol[resistor.v]","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"or we can plot the timeseries of the resistor's voltage:","category":"page"},{"location":"tutorials/acausal_components/#","page":"Acausal Component-Based Modeling the RC Circuit","title":"Acausal Component-Based Modeling the RC Circuit","text":"plot(sol, vars=[resistor.v])","category":"page"},{"location":"basics/Composition/#components-1","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"","category":"section"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"The symbolic models of ModelingToolkit can be composed together to easily build large models. The composition is lazy and only instantiated at the time of conversion to numerical models, allowing a more performant way in terms of computation time and memory.","category":"page"},{"location":"basics/Composition/#Simple-Model-Composition-Example-1","page":"Composing Models and Building Reusable Components","title":"Simple Model Composition Example","text":"","category":"section"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"The following is an example of building a model in a library with an optional forcing function, and allowing the user to specify the forcing later. Here, the library author defines a component named decay. The user then builds two decay components and connects them, saying the forcing term of decay1 is a constant while the forcing term of decay2 is the value of the state variable x.","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"using ModelingToolkit\n\nfunction decay(;name)\n  @parameters t a\n  @variables x(t) f(t)\n  D = Differential(t)\n  ODESystem([\n      D(x) ~ -a*x + f\n    ];\n    name=name)\nend\n\n@named decay1 = decay()\n@named decay2 = decay()\n\n@parameters t\nD = Differential(t)\n@named connected = ODESystem([\n                        decay2.f ~ decay1.x\n                        D(decay1.f) ~ 0\n                      ], t, systems=[decay1, decay2])\n\nequations(connected)\n\n#4-element Vector{Equation}:\n# Differential(t)(decay1₊f(t)) ~ 0\n# decay2₊f(t) ~ decay1₊x(t)\n# Differential(t)(decay1₊x(t)) ~ decay1₊f(t) - (decay1₊a*(decay1₊x(t)))\n# Differential(t)(decay2₊x(t)) ~ decay2₊f(t) - (decay2₊a*(decay2₊x(t)))\n\nsimplified_sys = structural_simplify(connected)\n\nequations(simplified_sys)\n\n#3-element Vector{Equation}:\n# Differential(t)(decay1₊f(t)) ~ 0\n# Differential(t)(decay1₊x(t)) ~ decay1₊f(t) - (decay1₊a*(decay1₊x(t)))\n# Differential(t)(decay2₊x(t)) ~ decay1₊x(t) - (decay2₊a*(decay2₊x(t)))","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"Now we can solve the system:","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"x0 = [\n  decay1.x => 1.0\n  decay1.f => 0.0\n  decay2.x => 1.0\n]\np = [\n  decay1.a => 0.1\n  decay2.a => 0.2\n]\n\nusing DifferentialEquations\nprob = ODEProblem(simplified_sys, x0, (0.0, 100.0), p)\nsol = solve(prob, Tsit5())\nsol[decay2.f]","category":"page"},{"location":"basics/Composition/#Basics-of-Model-Composition-1","page":"Composing Models and Building Reusable Components","title":"Basics of Model Composition","text":"","category":"section"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"Every AbstractSystem has a system keyword argument for specifying subsystems. A model is the composition of itself and its subsystems. For example, if we have:","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"@named sys = ODESystem(eqs,indepvar,states,ps,system=[subsys])","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"the equations of sys is the concatenation of get_eqs(sys) and equations(subsys), the states are the concatenation of their states, etc. When the ODEProblem or ODEFunction is generated from this system, it will build and compile the functions associated with this composition.","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"The new equations within the higher level system can access the variables in the lower level system by namespacing via the nameof(subsys). For example, let's say there is a variable x in states 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.","category":"page"},{"location":"basics/Composition/#Numerics-with-Composed-Models-1","page":"Composing Models and Building Reusable Components","title":"Numerics with Composed Models","text":"","category":"section"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"These composed models can then be directly transformed into their associated SciMLProblem type using the standard constructors. When this is done, the initial conditions and parameters must be specified in their namespaced form. For example:","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"u0 = [\n  x => 2.0\n  subsys.x => 2.0\n]","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"Note that any default values within the given subcomponent will be used if no override is provided at construction time. If any values for initial conditions or parameters are unspecified an error will be thrown.","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"When the model is numerically solved, the solution can be accessed via its symbolic values. For example, if sol is the ODESolution, one can use sol[x] and sol[subsys.x] to access the respective timeseries in the solution. All other indexing rules stay the same, so sol[x,1:5] accesses the first through fifth values of x. Note that this can be done even if the variable x is eliminated from the system from transformations like alias_elimination or tearing: the variable will be lazily reconstructed on demand.","category":"page"},{"location":"basics/Composition/#Variable-scope-and-parameter-expressions-1","page":"Composing Models and Building Reusable Components","title":"Variable scope and parameter expressions","text":"","category":"section"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"In some scenarios, it could be useful for model parameters to be expressed in terms of other parameters, or shared between common subsystems. To fascilitate this, ModelingToolkit supports sybmolic expressions in default values, and scoped variables.","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"With symbolic parameters, it is possible to set the default value of a parameter or initial condition to an expression of other variables.","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"# ...\nsys = ODESystem(\n    # ...\n    # directly in the defauls argument\n    defaults=Pair{Num, Any}[\n    x => u,\n    y => σ,\n    z => u-0.1,\n])\n# by assigning to the parameter\nsys.y = u*1.1","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"In a hierarchical system, variables of the subsystem get namespaced by the name of the system they are in. This prevents naming clashes, but also enforces that every state and parameter is local to the subsystem it is used in. In some cases it might be desirable to have variables and parameters that are shared between subsystems, or even global. This can be accomplished as follows.","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"@variables a b c d\n\n# a is a local variable\nb = ParentScope(b) # b is a variable that belongs to one level up in the hierarchy\nc = ParentScope(ParentScope(c)) # ParentScope can be nested\nd = GlobalScope(d) # global variables will never be namespaced","category":"page"},{"location":"basics/Composition/#Structural-Simplify-1","page":"Composing Models and Building Reusable Components","title":"Structural Simplify","text":"","category":"section"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"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 these trivial equality relationships and trivial singularity equations, i.e. equations which result in 0~0 expressions, in over-specified systems.","category":"page"},{"location":"basics/Composition/#Inheritance-and-Combine-(TODO)-1","page":"Composing Models and Building Reusable Components","title":"Inheritance and Combine (TODO)","text":"","category":"section"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"Model inheritance can be done in two ways: explicitly or implicitly. The explicit way is to shadow variables with equality expressions. For example, let's assume we have three separate systems which we want to compose to a single one. This is how one could explicitly forward all states and parameters to the higher level system:","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"using ModelingToolkit, OrdinaryDiffEq, Plots\n\n## Library code\n\n@parameters t\nD = Differential(t)\n\n@variables S(t), I(t), R(t)\nN = S + I + R\n@parameters β,γ\n\n@named seqn = ODESystem([D(S) ~ -β*S*I/N])\n@named ieqn = ODESystem([D(I) ~ β*S*I/N-γ*I])\n@named reqn = ODESystem([D(R) ~ γ*I])\n\n@named sir = ODESystem([   \n                    S ~ ieqn.S,\n                    I ~ seqn.I,\n                    R ~ ieqn.R,\n                    ieqn.S ~ seqn.S,\n                    seqn.I ~ ieqn.I,\n                    seqn.R ~ reqn.R,\n                    ieqn.R ~ reqn.R,\n                    reqn.I ~ ieqn.I], t, [S,I,R], [β,γ],\n                    systems=[seqn,ieqn,reqn],\n                    default_p = [\n                        seqn.β => β\n                        ieqn.β => β\n                        ieqn.γ => γ\n                        reqn.γ => γ\n                    ])","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"Note that the states are forwarded by an equality relationship, while the parameters are forwarded through a relationship in their default values. The user of this model can then solve this model simply by specifying the values at the highest level:","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"sireqn_simple = structural_simplify(sir)\n\nequations(sireqn_simple)\n\n# 3-element Vector{Equation}:\n#Differential(t)(seqn₊S(t)) ~ -seqn₊β*ieqn₊I(t)*seqn₊S(t)*(((ieqn₊I(t)) + (reqn₊R(t)) + (seqn₊S(t)))^-1)\n#Differential(t)(ieqn₊I(t)) ~ ieqn₊β*ieqn₊I(t)*seqn₊S(t)*(((ieqn₊I(t)) + (reqn₊R(t)) + (seqn₊S(t)))^-1) - (ieqn₊γ*(ieqn₊I(t)))\n#Differential(t)(reqn₊R(t)) ~ reqn₊γ*ieqn₊I(t)\n\n## User Code\n\nu0 = [seqn.S => 990.0,\n      ieqn.I => 10.0,\n      reqn.R => 0.0]\n\np = [\n    β => 0.5\n    γ => 0.25\n]\n\ntspan = (0.0,40.0)\nprob = ODEProblem(sireqn_simple,u0,tspan,p,jac=true)\nsol = solve(prob,Tsit5())\nsol[reqn.R]","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"However, one can similarly simplify this process of inheritance by using combine which concatenates all of the vectors within the systems. For example, we could equivalently have done:","category":"page"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"@named sir = combine([seqn,ieqn,reqn])","category":"page"},{"location":"basics/Composition/#Tearing-Problem-Construction-1","page":"Composing Models and Building Reusable Components","title":"Tearing Problem Construction","text":"","category":"section"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"Some system types, specifically ODESystem and NonlinearSystem, can be further reduced if structural_simplify has already been applied to them. This is done by using the alternative problem constructors, ODAEProblem and BlockNonlinearProblem respectively. In these cases, the constructor uses the knowledge of the strongly connected components calculated during the process of simplification as the basis for building pre-simplified nonlinear systems in the implicit solving. In summary: these problems are structurally modified, but could be more efficient and more stable.","category":"page"},{"location":"basics/Composition/#Automatic-Model-Promotion-(TODO)-1","page":"Composing Models and Building Reusable Components","title":"Automatic Model Promotion (TODO)","text":"","category":"section"},{"location":"basics/Composition/#","page":"Composing Models and Building Reusable Components","title":"Composing Models and Building Reusable Components","text":"In many cases one might want to compose models of different types. For example, one may want to include a NonlinearSystem as a set of algebraic equations within an ODESystem, or one may want to use an ODESystem as a subsystem of an SDESystem. In these cases, the compostion works automatically by promoting the model via promote_system. System promotions exist in the cases where a mathematically-trivial definition of the promotion exists.","category":"page"},{"location":"tutorials/stochastic_diffeq/#Modeling-with-Stochasticity-1","page":"Modeling with Stochasticity","title":"Modeling with Stochasticity","text":"","category":"section"},{"location":"tutorials/stochastic_diffeq/#","page":"Modeling with Stochasticity","title":"Modeling with Stochasticity","text":"All models with ODESystem are deterministic. SDESystem adds another element to the model: randomness. This is a stochastic differential equation which has a deterministic (drift) component and a stochastic (diffusion) component. Let's take the Lorenz equation from the first tutorial and extend it to have multiplicative noise.","category":"page"},{"location":"tutorials/stochastic_diffeq/#","page":"Modeling with Stochasticity","title":"Modeling with Stochasticity","text":"using ModelingToolkit, StochasticDiffEq\n\n# Define some variables\n@parameters t σ ρ β\n@variables x(t) y(t) z(t)\nD = Differential(t)\n\neqs = [D(x) ~ σ*(y-x),\n       D(y) ~ x*(ρ-z)-y,\n       D(z) ~ x*y - β*z]\n\nnoiseeqs = [0.1*x,\n            0.1*y,\n            0.1*z]\n\nde = SDESystem(eqs,noiseeqs,t,[x,y,z],[σ,ρ,β])\n\nu0map = [\n    x => 1.0,\n    y => 0.0,\n    z => 0.0\n]\n\nparammap = [\n    σ => 10.0,\n    β => 26.0,\n    ρ => 2.33\n]\n\nprob = SDEProblem(de,u0map,(0.0,100.0),parammap)\nsol = solve(prob,SOSRI())","category":"page"},{"location":"basics/Validation/#Model-Validation-and-Units-1","page":"Model Validation and Units","title":"Model Validation and Units","text":"","category":"section"},{"location":"basics/Validation/#","page":"Model Validation and Units","title":"Model Validation and Units","text":"ModelingToolkit.jl provides extensive functionality for model validation and unit checking. This is done by providing metadata to the variable types and then running the validation functions which identify malformed systems and non-physical equations.","category":"page"},{"location":"basics/Validation/#Consistency-Checking-1","page":"Model Validation and Units","title":"Consistency Checking","text":"","category":"section"},{"location":"basics/Validation/#","page":"Model Validation and Units","title":"Model Validation and Units","text":"check_consistency","category":"page"},{"location":"basics/Validation/#Unit-and-Type-Validation-1","page":"Model Validation and Units","title":"Unit and Type Validation","text":"","category":"section"},{"location":"basics/Validation/#","page":"Model Validation and Units","title":"Model Validation and Units","text":"ModelingToolkit.validate","category":"page"},{"location":"systems/SDESystem/#SDESystem-1","page":"SDESystem","title":"SDESystem","text":"","category":"section"},{"location":"systems/SDESystem/#System-Constructors-1","page":"SDESystem","title":"System Constructors","text":"","category":"section"},{"location":"systems/SDESystem/#","page":"SDESystem","title":"SDESystem","text":"SDESystem","category":"page"},{"location":"systems/SDESystem/#ModelingToolkit.SDESystem","page":"SDESystem","title":"ModelingToolkit.SDESystem","text":"struct SDESystem <: ModelingToolkit.AbstractODESystem\n\nA system of stochastic differential equations.\n\nFields\n\neqs\nThe expressions defining the drift term.\nnoiseeqs\nThe expressions defining the diffusion term.\niv\nIndependent variable.\nstates\nDependent (state) variables.\nps\nParameter variables.\nobserved\ntgrad\nTime-derivative matrix. Note: this field will not be defined until calculate_tgrad is called on the system.\n\njac\nJacobian matrix. Note: this field will not be defined until calculate_jacobian is called on the system.\n\nWfact\nWfact matrix. Note: this field will not be defined until generate_factorized_W is called on the system.\n\nWfact_t\nWfact_t matrix. Note: this field will not be defined until generate_factorized_W is called on the system.\n\nname\nName: the name of the system\n\nsystems\nSystems: the internal systems. These are required to have unique names.\n\ndefaults\ndefaults: The default values to use when initial conditions and/or parameters are not supplied in ODEProblem.\n\nconnection_type\ntype: type of the system\n\nExample\n\nusing ModelingToolkit\n\n@parameters t σ ρ β\n@variables x(t) y(t) z(t)\nD = Differential(t)\n\neqs = [D(x) ~ σ*(y-x),\n       D(y) ~ x*(ρ-z)-y,\n       D(z) ~ x*y - β*z]\n\nnoiseeqs = [0.1*x,\n            0.1*y,\n            0.1*z]\n\nde = SDESystem(eqs,noiseeqs,t,[x,y,z],[σ,ρ,β])\n\n\n\n\n\n","category":"type"},{"location":"systems/SDESystem/#Composition-and-Accessor-Functions-1","page":"SDESystem","title":"Composition and Accessor Functions","text":"","category":"section"},{"location":"systems/SDESystem/#","page":"SDESystem","title":"SDESystem","text":"get_eqs(sys) or equations(sys): The equations that define the SDE.\nget_states(sys) or states(sys): The set of states in the SDE.\nget_ps(sys)s or parameters(sys): The parameters of the SDE.\nindependent_variable(sys): The independent variable of the SDE.","category":"page"},{"location":"systems/SDESystem/#Transformations-1","page":"SDESystem","title":"Transformations","text":"","category":"section"},{"location":"systems/SDESystem/#","page":"SDESystem","title":"SDESystem","text":"structural_simplify\r\nalias_elimination","category":"page"},{"location":"systems/SDESystem/#Analyses-1","page":"SDESystem","title":"Analyses","text":"","category":"section"},{"location":"systems/SDESystem/#Applicable-Calculation-and-Generation-Functions-1","page":"SDESystem","title":"Applicable Calculation and Generation Functions","text":"","category":"section"},{"location":"systems/SDESystem/#","page":"SDESystem","title":"SDESystem","text":"calculate_jacobian\r\ncalculate_tgrad\r\ncalculate_factorized_W\r\ngenerate_jacobian\r\ngenerate_tgrad\r\ngenerate_factorized_W\r\njacobian_sparsity","category":"page"},{"location":"systems/SDESystem/#Problem-Constructors-1","page":"SDESystem","title":"Problem Constructors","text":"","category":"section"},{"location":"systems/SDESystem/#","page":"SDESystem","title":"SDESystem","text":"SDEFunction\r\nSDEProblem","category":"page"},{"location":"systems/SDESystem/#SciMLBase.SDEFunction","page":"SDESystem","title":"SciMLBase.SDEFunction","text":"function DiffEqBase.SDEFunction{iip}(sys::SDESystem, dvs = sys.states, ps = sys.ps;\n                                     version = nothing, tgrad=false, sparse = false,\n                                     jac = false, Wfact = false, kwargs...) where {iip}\n\nCreate an SDEFunction from the SDESystem. The arguments dvs and ps are used to set the order of the dependent variable and parameter vectors, respectively.\n\n\n\n\n\n","category":"type"},{"location":"systems/SDESystem/#SciMLBase.SDEProblem","page":"SDESystem","title":"SciMLBase.SDEProblem","text":"function DiffEqBase.SDEProblem{iip}(sys::SDESystem,u0map,tspan,p=parammap;\n                                    version = nothing, tgrad=false,\n                                    jac = false, Wfact = false,\n                                    checkbounds = false, sparse = false,\n                                    sparsenoise = sparse,\n                                    skipzeros = true, fillzeros = true,\n                                    linenumbers = true, parallel=SerialForm(),\n                                    kwargs...)\n\nGenerates an SDEProblem from an SDESystem and allows for automatically symbolically calculating numerical enhancements.\n\n\n\n\n\n","category":"type"},{"location":"internals/#Internal-Details-1","page":"Internal Details","title":"Internal Details","text":"","category":"section"},{"location":"internals/#","page":"Internal Details","title":"Internal Details","text":"This is a page for detailing some of the inner workings to help future contributors to the library.","category":"page"},{"location":"internals/#Observables-and-Variable-Elimination-1","page":"Internal Details","title":"Observables and Variable Elimination","text":"","category":"section"},{"location":"internals/#","page":"Internal Details","title":"Internal Details","text":"In the variable \"elimination\" algorithms, what is actually done is that variables are removed from being states and equations are moved into the observed category of the system. The observed equations are explicit algebraic equations which are then substituted out to completely eliminate these variables from the other equations, allowing the system to act as though these variables no longer exist.","category":"page"},{"location":"internals/#","page":"Internal Details","title":"Internal Details","text":"However, as a user may have wanted to interact with such variables, for example, plotting their output, these relationships are stored and are then used to generate the observed equation found in the SciMLFunction interface, so that sol[x] lazily reconstructs the observed variable when necessary. In this sense, there is an equivalence between observables and the variable elimination system.","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#Automated-Index-Reduction-of-DAEs-1","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"","category":"section"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"In many cases one may accidentally write down a DAE that is not easily solvable by numerical methods. In this tutorial we will walk through an example of a pendulum which accidentally generates an index-3 DAE, and show how to use the modelingtoolkitize to correct the model definition before solving.","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#Copy-Pastable-Example-1","page":"Automated Index Reduction of DAEs","title":"Copy-Pastable Example","text":"","category":"section"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"using ModelingToolkit\nusing LinearAlgebra\nusing OrdinaryDiffEq\nusing Plots\n\nfunction pendulum!(du, u, p, t)\n    x, dx, y, dy, T = u\n    g, L = p\n    du[1] = dx\n    du[2] = T*x\n    du[3] = dy\n    du[4] = T*y - g\n    du[5] = x^2 + y^2 - L^2\n    return nothing\nend\npendulum_fun! = ODEFunction(pendulum!, mass_matrix=Diagonal([1,1,1,1,0]))\nu0 = [1.0, 0, 0, 0, 0]\np = [9.8, 1]\ntspan = (0, 10.0)\npendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)\ntraced_sys = modelingtoolkitize(pendulum_prob)\npendulum_sys = structural_simplify(dae_index_lowering(traced_sys))\nprob = ODAEProblem(pendulum_sys, Pair[], tspan)\nsol = solve(prob, Tsit5(),abstol=1e-8,reltol=1e-8)\nplot(sol, vars=states(traced_sys))","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#Explanation-1","page":"Automated Index Reduction of DAEs","title":"Explanation","text":"","category":"section"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#Attempting-to-Solve-the-Equation-1","page":"Automated Index Reduction of DAEs","title":"Attempting to Solve the Equation","text":"","category":"section"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"In this tutorial we will look at the pendulum system:","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"beginaligned\n    x^prime = v_x\n    v_x^prime = Tx\n    y^prime = v_y\n    v_y^prime = Ty - g\n    0 = x^2 + y^2 - L^2\nendaligned","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"As a good DifferentialEquations.jl user, one would follow the mass matrix DAE tutorial to arrive at code for simulating the model:","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"using OrdinaryDiffEq, LinearAlgebra\nfunction pendulum!(du, u, p, t)\n    x, dx, y, dy, T = u\n    g, L = p\n    du[1] = dx; du[2] = T*x\n    du[3] = dy; du[4] = T*y - g\n    du[5] = x^2 + y^2 - L^2\nend\npendulum_fun! = ODEFunction(pendulum!, mass_matrix=Diagonal([1,1,1,1,0]))\nu0 = [1.0, 0, 0, 0, 0]; p = [9.8, 1]; tspan = (0, 10.0)\npendulum_prob = ODEProblem(pendulum_fun!, u0, tspan, p)\nsolve(pendulum_prob,Rodas4())","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"However, one will quickly be greeted with the unfortunate message:","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"┌ Warning: First function call produced NaNs. Exiting.\n└ @ OrdinaryDiffEq C:\\Users\\accou\\.julia\\packages\\OrdinaryDiffEq\\yCczp\\src\\initdt.jl:76\n┌ Warning: Automatic dt set the starting dt as NaN, causing instability.\n└ @ OrdinaryDiffEq C:\\Users\\accou\\.julia\\packages\\OrdinaryDiffEq\\yCczp\\src\\solve.jl:485\n┌ Warning: NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.\n└ @ SciMLBase C:\\Users\\accou\\.julia\\packages\\SciMLBase\\DrPil\\src\\integrator_interface.jl:325","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"Did you implement the DAE incorrectly? No. Is the solver broken? No.","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#Understanding-DAE-Index-1","page":"Automated Index Reduction of DAEs","title":"Understanding DAE Index","text":"","category":"section"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"It turns out that this is a property of the DAE that we are attempting to solve. This kind of DAE is known as an index-3 DAE. For a complete discussion of DAE index, see this article. Essentially the issue here is that we have 4 differential variables (x, v_x, y, v_y) and one algebraic variable T (which we can know because there is no D(T) term in the equations). An index-1 DAE always satisfies that the Jacobian of the algebraic equations is non-singular. Here, the first 4 equations are differential equations, with the last term the algebraic relationship. However, the partial derivative of x^2 + y^2 - L^2 w.r.t. T is zero, and thus the Jacobian of the algebraic equations is the zero matrix and thus it's singular. This is a very quick way to see whether the DAE is index 1!","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"The problem with higher order DAEs is that the matrices used in Newton solves are singular or close to singular when applied to such problems. Because of this fact, the nonlinear solvers (or Rosenbrock methods) break down, making them difficult to solve. The classic paper DAEs are not ODEs goes into detail on this and shows that many methods are no longer convergent when index is higher than one. So it's not necessarily the fault of the solver or the implementation: this is known.","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"But that's not a satisfying answer, so what do you do about it?","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#Transforming-Higher-Order-DAEs-to-Index-1-DAEs-1","page":"Automated Index Reduction of DAEs","title":"Transforming Higher Order DAEs to Index-1 DAEs","text":"","category":"section"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"It turns out that higher order DAEs can be transformed into lower order DAEs. If you differentiate the last equation two times and perform a substitution, you can arrive at the following set of equations:","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"beginaligned\nx^prime = v_x \nv_x^prime = x T \ny^prime = v_y \nv_y^prime = y T - g \n0 = 2 left(v_x^2 + v_y^2 + y ( y T - g ) + T x^2 right)\nendaligned","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"Note that this is mathematically-equivalent to the equation that we had before, but the Jacobian w.r.t. T of the algebraic equation is no longer zero because of the substitution. This means that if you wrote down this version of the model it will be index-1 and solve correctly! In fact, this is how DAE index is commonly defined: the number of differentiations it takes to transform the DAE into an ODE, where an ODE is an index-0 DAE by substituting out all of the algebraic relationships.","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#Automating-the-Index-Reduction-1","page":"Automated Index Reduction of DAEs","title":"Automating the Index Reduction","text":"","category":"section"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"However, requiring the user to sit there and work through this process on potentially millions of equations is an unfathomable mental overhead. But, we can avoid this by using methods like the Pantelides algorithm for automatically performing this reduction to index 1. While this requires the ModelingToolkit symbolic form, we use modelingtoolkitize to transform the numerical code into symbolic code, run dae_index_lowering lowering, then transform back to numerical code with ODEProblem, and solve with a numerical solver. Let's try that out:","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"traced_sys = modelingtoolkitize(pendulum_prob)\npendulum_sys = structural_simplify(dae_index_lowering(traced_sys))\nprob = ODEProblem(pendulum_sys, Pair[], tspan)\nsol = solve(prob, Rodas4())\n\nusing Plots\nplot(sol, vars=states(traced_sys))","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"(Image: )","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"Note that plotting using states(traced_sys) is done so that any variables which are symbolically eliminated, or any variable reorderings done for enhanced parallelism/performance, still show up in the resulting plot and the plot is shown in the same order as the original numerical code.","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"Note that we can even go a little bit further. If we use the ODAEProblem constructor, we can remove the algebraic equations from the states of the system and fully transform the index-3 DAE into an index-0 ODE which can be solved via an explicit Runge-Kutta method:","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"traced_sys = modelingtoolkitize(pendulum_prob)\npendulum_sys = structural_simplify(dae_index_lowering(traced_sys))\nprob = ODAEProblem(pendulum_sys, Pair[], tspan)\nsol = solve(prob, Tsit5(),abstol=1e-8,reltol=1e-8)\nplot(sol, vars=states(traced_sys))","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"(Image: )","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize_index_reduction/#","page":"Automated Index Reduction of DAEs","title":"Automated Index Reduction of DAEs","text":"And there you go: this has transformed the model from being too hard to solve with implicit DAE solvers, to something that is easily solved with explicit Runge-Kutta methods for non-stiff equations.","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize/#Automatically-Accelerating-ODEProblem-Code-1","page":"Automatically Accelerating ODEProblem Code","title":"Automatically Accelerating ODEProblem Code","text":"","category":"section"},{"location":"mtkitize_tutorials/modelingtoolkitize/#","page":"Automatically Accelerating ODEProblem Code","title":"Automatically Accelerating ODEProblem Code","text":"For some DEProblem types, automatic tracing functionality is already included via the modelingtoolkitize function. Take, for example, the Robertson ODE defined as an ODEProblem for DifferentialEquations.jl:","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize/#","page":"Automatically Accelerating ODEProblem Code","title":"Automatically Accelerating ODEProblem Code","text":"using DifferentialEquations\nfunction rober(du,u,p,t)\n  y₁,y₂,y₃ = u\n  k₁,k₂,k₃ = p\n  du[1] = -k₁*y₁+k₃*y₂*y₃\n  du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃\n  du[3] =  k₂*y₂^2\n  nothing\nend\nprob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize/#","page":"Automatically Accelerating ODEProblem Code","title":"Automatically Accelerating ODEProblem Code","text":"If we want to get a symbolic representation, we can simply call modelingtoolkitize on the prob, which will return an ODESystem:","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize/#","page":"Automatically Accelerating ODEProblem Code","title":"Automatically Accelerating ODEProblem Code","text":"sys = modelingtoolkitize(prob)","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize/#","page":"Automatically Accelerating ODEProblem Code","title":"Automatically Accelerating ODEProblem Code","text":"Using this, we can symbolically build the Jacobian and then rebuild the ODEProblem:","category":"page"},{"location":"mtkitize_tutorials/modelingtoolkitize/#","page":"Automatically Accelerating ODEProblem Code","title":"Automatically Accelerating ODEProblem Code","text":"jac = eval(ModelingToolkit.generate_jacobian(sys)[2])\nf = ODEFunction(rober, jac=jac)\nprob_jac = ODEProblem(f,[1.0,0.0,0.0],(0.0,1e5),(0.04,3e7,1e4))","category":"page"},{"location":"systems/PDESystem/#PDESystem-1","page":"PDESystem","title":"PDESystem","text":"","category":"section"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"PDESystem is the common symbolic PDE specification for the SciML ecosystem. It is currently being built as a component of the ModelingToolkit ecosystem,","category":"page"},{"location":"systems/PDESystem/#Vision-1","page":"PDESystem","title":"Vision","text":"","category":"section"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"The vision for the common PDE interface is that a user should only have to specify their PDE once, mathematically, and have instant access to everything as simple as a finite difference method with constant grid spacing, to something as complex as a distributed multi-GPU discrete Galerkin method.","category":"page"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"The key to the common PDE interface is a separation of the symbolic handling from the numerical world. All of the discretizers should not \"solve\" the PDE, but instead be a conversion of the mathematical specification to a numerical problem. Preferably, the transformation should be to another ModelingToolkit.jl AbstractSystem, but in some cases this cannot be done or will not be performant, so a SciMLProblem is the other choice.","category":"page"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"These elementary problems, such as solving linear systems Ax=b, solving nonlinear systems f(x)=0, ODEs, etc. are all defined by SciMLBase.jl, which then numerical solvers can all target these common forms. Thus someone who works on linear solvers doesn't necessarily need to be working on a Discontinuous Galerkin or finite element library, but instead \"linear solvers that are good for matrices A with properties ...\" which are then accessible by every other discretization method in the common PDE interface.","category":"page"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"Similar to the rest of the AbstractSystem types, transformation and analyses functions will allow for simplifying the PDE before solving it, and constructing block symbolic functions like Jacobians.","category":"page"},{"location":"systems/PDESystem/#Constructors-1","page":"PDESystem","title":"Constructors","text":"","category":"section"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"PDESystem","category":"page"},{"location":"systems/PDESystem/#ModelingToolkit.PDESystem","page":"PDESystem","title":"ModelingToolkit.PDESystem","text":"struct PDESystem <: ModelingToolkit.AbstractSystem\n\nA system of partial differential equations.\n\nFields\n\neqs\nThe equations which define the PDE\nbcs\nThe boundary conditions\ndomain\nThe domain for the independent variables.\nindvars\nThe independent variables\ndepvars\nThe dependent variables\nps\nThe parameters\ndefaults\ndefaults: The default values to use when initial conditions and/or parameters are not supplied in ODEProblem.\n\nconnection_type\ntype: type of the system\n\nExample\n\nusing ModelingToolkit\n\n@parameters t x\n@variables u(..)\nDxx = Differential(x)^2\nDtt = Differential(t)^2\nDt = Differential(t)\n\n#2D PDE\nC=1\neq  = Dtt(u(t,x)) ~ C^2*Dxx(u(t,x))\n\n# Initial and boundary conditions\nbcs = [u(t,0) ~ 0.,# for all t > 0\n       u(t,1) ~ 0.,# for all t > 0\n       u(0,x) ~ x*(1. - x), #for all 0 < x < 1\n       Dt(u(0,x)) ~ 0. ] #for all  0 < x < 1]\n\n# Space and time domains\ndomains = [t ∈ (0.0,1.0),\n           x ∈ (0.0,1.0)]\n\npde_system = PDESystem(eq,bcs,domains,[t,x],[u])\n\n\n\n\n\n","category":"type"},{"location":"systems/PDESystem/#Domains-(WIP)-1","page":"PDESystem","title":"Domains (WIP)","text":"","category":"section"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"Domains are specifying by saying indepvar in domain, where indepvar is a single or a collection of independent variables, and domain is the chosen domain type. A 2-tuple can be used to indicate an Interval.  Thus forms for the indepvar can be like:","category":"page"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"t ∈ (0.0,1.0)\r\n(t,x) ∈ UnitDisk()\r\n[v,w,x,y,z] ∈ VectorUnitBall(5)","category":"page"},{"location":"systems/PDESystem/#Domain-Types-(WIP)-1","page":"PDESystem","title":"Domain Types (WIP)","text":"","category":"section"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"Interval(a,b): Defines the domain of an interval from a to b (requires explicit","category":"page"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"import from DomainSets.jl, but a 2-tuple can be used instead)","category":"page"},{"location":"systems/PDESystem/#discretize-and-symbolic_discretize-1","page":"PDESystem","title":"discretize and symbolic_discretize","text":"","category":"section"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"The only functions which act on a PDESystem are the following:","category":"page"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"discretize(sys,discretizer): produces the outputted AbstractSystem or SciMLProblem.\nsymbolic_discretize(sys,discretizer): produces a debugging symbolic description of the discretized problem.","category":"page"},{"location":"systems/PDESystem/#Boundary-Conditions-(WIP)-1","page":"PDESystem","title":"Boundary Conditions (WIP)","text":"","category":"section"},{"location":"systems/PDESystem/#Transformations-1","page":"PDESystem","title":"Transformations","text":"","category":"section"},{"location":"systems/PDESystem/#Analyses-1","page":"PDESystem","title":"Analyses","text":"","category":"section"},{"location":"systems/PDESystem/#Discretizer-Ecosystem-1","page":"PDESystem","title":"Discretizer Ecosystem","text":"","category":"section"},{"location":"systems/PDESystem/#NeuralPDE.jl:-PhysicsInformedNN-1","page":"PDESystem","title":"NeuralPDE.jl: PhysicsInformedNN","text":"","category":"section"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"NeuralPDE.jl defines the PhysicsInformedNN discretizer which uses a DiffEqFlux.jl neural network to solve the differential equation.","category":"page"},{"location":"systems/PDESystem/#DiffEqOperators.jl:-MOLFiniteDifference-(WIP)-1","page":"PDESystem","title":"DiffEqOperators.jl: MOLFiniteDifference (WIP)","text":"","category":"section"},{"location":"systems/PDESystem/#","page":"PDESystem","title":"PDESystem","text":"DiffEqOperators.jl defines the MOLFiniteDifference discretizer which performs a finite difference discretization using the DiffEqOperators.jl stencils. These stencils make use of NNLib.jl for fast operations on semi-linear domains.","category":"page"},{"location":"systems/OptimizationSystem/#OptimizationSystem-1","page":"OptimizationSystem","title":"OptimizationSystem","text":"","category":"section"},{"location":"systems/OptimizationSystem/#System-Constructors-1","page":"OptimizationSystem","title":"System Constructors","text":"","category":"section"},{"location":"systems/OptimizationSystem/#","page":"OptimizationSystem","title":"OptimizationSystem","text":"OptimizationSystem","category":"page"},{"location":"systems/OptimizationSystem/#ModelingToolkit.OptimizationSystem","page":"OptimizationSystem","title":"ModelingToolkit.OptimizationSystem","text":"struct OptimizationSystem <: ModelingToolkit.AbstractSystem\n\nA scalar equation for optimization.\n\nFields\n\nop\nVector of equations defining the system.\nstates\nUnknown variables.\nps\nParameters.\nobserved\nequality_constraints\ninequality_constraints\nname\nName: the name of the system.  These are required to have unique names.\n\nsystems\nsystems: The internal systems\n\ndefaults\ndefaults: The default values to use when initial conditions and/or parameters are not supplied in ODEProblem.\n\nExamples\n\n@variables x y z\n@parameters σ ρ β\n\nop = σ*(y-x) + x*(ρ-z)-y + x*y - β*z\nos = OptimizationSystem(eqs, [x,y,z],[σ,ρ,β])\n\n\n\n\n\n","category":"type"},{"location":"systems/OptimizationSystem/#Composition-and-Accessor-Functions-1","page":"OptimizationSystem","title":"Composition and Accessor Functions","text":"","category":"section"},{"location":"systems/OptimizationSystem/#","page":"OptimizationSystem","title":"OptimizationSystem","text":"get_eqs(sys) or equations(sys): The equation to be minimized.\nget_states(sys) or states(sys): The set of states for the optimization.\nget_ps(sys) or parameters(sys): The parameters for the optimization.","category":"page"},{"location":"systems/OptimizationSystem/#Transformations-1","page":"OptimizationSystem","title":"Transformations","text":"","category":"section"},{"location":"systems/OptimizationSystem/#Analyses-1","page":"OptimizationSystem","title":"Analyses","text":"","category":"section"},{"location":"systems/OptimizationSystem/#Applicable-Calculation-and-Generation-Functions-1","page":"OptimizationSystem","title":"Applicable Calculation and Generation Functions","text":"","category":"section"},{"location":"systems/OptimizationSystem/#","page":"OptimizationSystem","title":"OptimizationSystem","text":"calculate_gradient\r\ncalculate_hessian\r\ngenerate_gradient\r\ngenerate_hessian\r\nhessian_sparsity","category":"page"},{"location":"systems/OptimizationSystem/#Problem-Constructors-1","page":"OptimizationSystem","title":"Problem Constructors","text":"","category":"section"},{"location":"systems/OptimizationSystem/#","page":"OptimizationSystem","title":"OptimizationSystem","text":"OptimizationProblem","category":"page"},{"location":"systems/OptimizationSystem/#SciMLBase.OptimizationProblem","page":"OptimizationSystem","title":"SciMLBase.OptimizationProblem","text":"function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem,\n                                          parammap=DiffEqBase.NullParameters();\n                                          u0=nothing, lb=nothing, ub=nothing,\n                                          grad = false,\n                                          hess = false, sparse = false,\n                                          checkbounds = false,\n                                          linenumbers = true, parallel=SerialForm(),\n                                          kwargs...) where iip\n\nGenerates an OptimizationProblem from an OptimizationSystem and allows for automatically symbolically calculating numerical enhancements.\n\n\n\n\n\n","category":"type"},{"location":"tutorials/higher_order/#Automatic-Transformation-of-Nth-Order-ODEs-to-1st-Order-ODEs-1","page":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","title":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","text":"","category":"section"},{"location":"tutorials/higher_order/#","page":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","title":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","text":"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 the ode_order_lowering transformation that sends an Nth order ODE to a 1st order ODE.","category":"page"},{"location":"tutorials/higher_order/#","page":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","title":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","text":"To see this, let's define a second order riff on the Lorenz equations. We utilize the derivative operator twice here to define the second order:","category":"page"},{"location":"tutorials/higher_order/#","page":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","title":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","text":"using ModelingToolkit, OrdinaryDiffEq\r\n\r\n@parameters t σ ρ β\r\n@variables x(t) y(t) z(t)\r\nD = Differential(t)\r\n\r\neqs = [D(D(x)) ~ σ*(y-x),\r\n       D(y) ~ x*(ρ-z)-y,\r\n       D(z) ~ x*y - β*z]\r\n\r\nsys = ODESystem(eqs)","category":"page"},{"location":"tutorials/higher_order/#","page":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","title":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","text":"Note that we could've used an alternative syntax for 2nd order, i.e. D = Differential(t)^2 and then E(x) would be the second derivative, and this syntax extends to N-th order. Also, we can use * or ∘ to compose Differentials, like Differential(t) * Differential(x).","category":"page"},{"location":"tutorials/higher_order/#","page":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","title":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","text":"Now let's transform this into the ODESystem of first order components. We do this by simply calling ode_order_lowering:","category":"page"},{"location":"tutorials/higher_order/#","page":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","title":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","text":"sys = ode_order_lowering(sys)","category":"page"},{"location":"tutorials/higher_order/#","page":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","title":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","text":"Now we can directly numerically solve the lowered system. Note that, following the original problem, the solution requires knowing the initial condition for x', and thus we include that in our input specification:","category":"page"},{"location":"tutorials/higher_order/#","page":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","title":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","text":"u0 = [D(x) => 2.0,\r\n      x => 1.0,\r\n      y => 0.0,\r\n      z => 0.0]\r\n\r\np  = [σ => 28.0,\r\n      ρ => 10.0,\r\n      β => 8/3]\r\n\r\ntspan = (0.0,100.0)\r\nprob = ODEProblem(sys,u0,tspan,p,jac=true)\r\nsol = solve(prob,Tsit5())\r\nusing Plots; plot(sol,vars=(x,y))","category":"page"},{"location":"tutorials/higher_order/#","page":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","title":"Automatic Transformation of Nth Order ODEs to 1st Order ODEs","text":"(Image: Lorenz2)","category":"page"},{"location":"systems/ODESystem/#ODESystem-1","page":"ODESystem","title":"ODESystem","text":"","category":"section"},{"location":"systems/ODESystem/#System-Constructors-1","page":"ODESystem","title":"System Constructors","text":"","category":"section"},{"location":"systems/ODESystem/#","page":"ODESystem","title":"ODESystem","text":"ODESystem","category":"page"},{"location":"systems/ODESystem/#ModelingToolkit.ODESystem","page":"ODESystem","title":"ModelingToolkit.ODESystem","text":"struct ODESystem <: ModelingToolkit.AbstractODESystem\n\nA system of ordinary differential equations.\n\nFields\n\neqs\nThe ODEs defining the system.\niv\nIndependent variable.\nstates\nDependent (state) variables.\nps\nParameter variables.\nobserved\ntgrad\nTime-derivative matrix. Note: this field will not be defined until calculate_tgrad is called on the system.\n\njac\nJacobian matrix. Note: this field will not be defined until calculate_jacobian is called on the system.\n\nWfact\nWfact matrix. Note: this field will not be defined until generate_factorized_W is called on the system.\n\nWfact_t\nWfact_t matrix. Note: this field will not be defined until generate_factorized_W is called on the system.\n\nname\nName: the name of the system\n\nsystems\nsystems: The internal systems. These are required to have unique names.\n\ndefaults\ndefaults: The default values to use when initial conditions and/or parameters are not supplied in ODEProblem.\n\nstructure\nstructure: structural information of the system\n\nconnection_type\ntype: type of the system\n\nExample\n\nusing ModelingToolkit\n\n@parameters t σ ρ β\n@variables x(t) y(t) z(t)\nD = Differential(t)\n\neqs = [D(x) ~ σ*(y-x),\n       D(y) ~ x*(ρ-z)-y,\n       D(z) ~ x*y - β*z]\n\nde = ODESystem(eqs,t,[x,y,z],[σ,ρ,β])\n\n\n\n\n\n","category":"type"},{"location":"systems/ODESystem/#Composition-and-Accessor-Functions-1","page":"ODESystem","title":"Composition and Accessor Functions","text":"","category":"section"},{"location":"systems/ODESystem/#","page":"ODESystem","title":"ODESystem","text":"get_eqs(sys) or equations(sys): The equations that define the ODE.\nget_states(sys) or states(sys): The set of states in the ODE.\nget_ps(sys) or parameters(sys): The parameters of the ODE.\nindependent_variable(sys): The independent variable of the ODE.","category":"page"},{"location":"systems/ODESystem/#Transformations-1","page":"ODESystem","title":"Transformations","text":"","category":"section"},{"location":"systems/ODESystem/#","page":"ODESystem","title":"ODESystem","text":"structural_simplify\r\node_order_lowering\r\ndae_index_lowering\r\nliouville_transform\r\nalias_elimination\r\ntearing","category":"page"},{"location":"systems/ODESystem/#ModelingToolkit.ode_order_lowering","page":"ODESystem","title":"ModelingToolkit.ode_order_lowering","text":"ode_order_lowering(sys::ODESystem) -> Any\n\n\nTakes a Nth order ODESystem and returns a new ODESystem written in first order form by defining new variables which represent the N-1 derivatives.\n\n\n\n\n\n","category":"function"},{"location":"systems/ODESystem/#ModelingToolkit.StructuralTransformations.dae_index_lowering","page":"ODESystem","title":"ModelingToolkit.StructuralTransformations.dae_index_lowering","text":"dae_index_lowering(sys::ODESystem) -> ODESystem\n\nPerform the Pantelides algorithm to transform a higher index DAE to an index 1 DAE.\n\n\n\n\n\n","category":"function"},{"location":"systems/ODESystem/#ModelingToolkit.liouville_transform","page":"ODESystem","title":"ModelingToolkit.liouville_transform","text":"liouville_transform(sys::Any) -> ODESystem\n\n\nGenerates the Liouville transformed set of ODEs, which is the original ODE system with a new variable trJ appended, corresponding to the -tr(Jacobian). This variable is used for properties like uncertainty propagation from a given initial distribution density.\n\nFor example, if u=p*u and p follows a probability distribution f(p), then the probability density of a future value with a given choice of p is computed by setting the inital trJ = f(p), and the final value of trJ is the probability of u(t).\n\nExample:\n\nusing ModelingToolkit, OrdinaryDiffEq, Test\n\n@parameters t α β γ δ\n@variables x(t) y(t)\nD = Differential(t)\n\neqs = [D(x) ~ α*x - β*x*y,\n       D(y) ~ -δ*y + γ*x*y]\n\nsys = ODESystem(eqs)\nsys2 = liouville_transform(sys)\n@variables trJ\n\nu0 = [x => 1.0,\n      y => 1.0,\n      trJ => 1.0]\n\nprob = ODEProblem(sys2,u0,tspan,p)\nsol = solve(prob,Tsit5())\n\nWhere sol[3,:] is the evolution of trJ over time.\n\nSources:\n\nProbabilistic Robustness Analysis of F-16 Controller Performance: An Optimal Transport Approach\n\nAbhishek Halder, Kooktae Lee, and Raktim Bhattacharya https://abhishekhalder.bitbucket.io/F16ACC2013Final.pdf\n\n\n\n\n\n","category":"function"},{"location":"systems/ODESystem/#Analyses-1","page":"ODESystem","title":"Analyses","text":"","category":"section"},{"location":"systems/ODESystem/#","page":"ODESystem","title":"ODESystem","text":"ModelingToolkit.islinear\r\nModelingToolkit.isautonomous","category":"page"},{"location":"systems/ODESystem/#Applicable-Calculation-and-Generation-Functions-1","page":"ODESystem","title":"Applicable Calculation and Generation Functions","text":"","category":"section"},{"location":"systems/ODESystem/#","page":"ODESystem","title":"ODESystem","text":"calculate_jacobian\r\ncalculate_tgrad\r\ncalculate_factorized_W\r\ngenerate_jacobian\r\ngenerate_tgrad\r\ngenerate_factorized_W\r\njacobian_sparsity","category":"page"},{"location":"systems/ODESystem/#Standard-Problem-Constructors-1","page":"ODESystem","title":"Standard Problem Constructors","text":"","category":"section"},{"location":"systems/ODESystem/#","page":"ODESystem","title":"ODESystem","text":"ODEFunction\r\nODEProblem\r\nSteadyStateFunction\r\nSteadyStateProblem","category":"page"},{"location":"systems/ODESystem/#SciMLBase.ODEFunction","page":"ODESystem","title":"SciMLBase.ODEFunction","text":"function DiffEqBase.ODEFunction{iip}(sys::AbstractODESystem, dvs = states(sys),\n                                     ps = parameters(sys);\n                                     version = nothing, tgrad=false,\n                                     jac = false,\n                                     sparse = false,\n                                     kwargs...) where {iip}\n\nCreate an ODEFunction from the ODESystem. The arguments dvs and ps are used to set the order of the dependent variable and parameter vectors, respectively.\n\n\n\n\n\n","category":"type"},{"location":"systems/ODESystem/#SciMLBase.ODEProblem","page":"ODESystem","title":"SciMLBase.ODEProblem","text":"function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem,u0map,tspan,\n                                    parammap=DiffEqBase.NullParameters();\n                                    version = nothing, tgrad=false,\n                                    jac = false,\n                                    checkbounds = false, sparse = false,\n                                    simplify=false,\n                                    linenumbers = true, parallel=SerialForm(),\n                                    kwargs...) where iip\n\nGenerates an ODEProblem from an ODESystem and allows for automatically symbolically calculating numerical enhancements.\n\n\n\n\n\n","category":"type"},{"location":"systems/ODESystem/#SciMLBase.SteadyStateProblem","page":"ODESystem","title":"SciMLBase.SteadyStateProblem","text":"function DiffEqBase.SteadyStateProblem(sys::AbstractODESystem,u0map,\n                                    parammap=DiffEqBase.NullParameters();\n                                    version = nothing, tgrad=false,\n                                    jac = false,\n                                    checkbounds = false, sparse = false,\n                                    linenumbers = true, parallel=SerialForm(),\n                                    kwargs...) where iip\n\nGenerates an SteadyStateProblem from an ODESystem and allows for automatically symbolically calculating numerical enhancements.\n\n\n\n\n\n","category":"type"},{"location":"systems/ODESystem/#Torn-Problem-Constructors-1","page":"ODESystem","title":"Torn Problem Constructors","text":"","category":"section"},{"location":"systems/ODESystem/#","page":"ODESystem","title":"ODESystem","text":"ODAEProblem","category":"page"},{"location":"tutorials/tearing_parallelism/#Exposing-More-Parallelism-By-Tearing-Algebraic-Equations-in-ODESystems-1","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"","category":"section"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"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 parallelism in the solution process and parallelize the resulting simulation.","category":"page"},{"location":"tutorials/tearing_parallelism/#The-Component-Library-1","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"The Component Library","text":"","category":"section"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"The following tutorial will use the following set of components describing electrical circuits:","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"using ModelingToolkit, OrdinaryDiffEq\n\nfunction connect_pin(ps...)\n    eqs = [\n           0 ~ sum(p->p.i, ps) # KCL\n          ]\n    # KVL\n    for i in 1:length(ps)-1\n        push!(eqs, ps[i].v ~ ps[i+1].v)\n    end\n\n    return eqs\nend\n\nfunction connect_heat(ps...)\n    eqs = [\n           0 ~ sum(p->p.Q_flow, ps)\n          ]\n    for i in 1:length(ps)-1\n        push!(eqs, ps[i].T ~ ps[i+1].T)\n    end\n\n    return eqs\nend\n\n# Basic electric components\n@parameters t\nconst D = Differential(t)\nfunction Pin(;name)\n    @variables v(t) i(t)\n    ODESystem(Equation[], t, [v, i], [], name=name, defaults=Dict([v=>1.0, i=>1.0]))\nend\n\nfunction Ground(;name)\n    @named g = Pin()\n    eqs = [g.v ~ 0]\n    ODESystem(eqs, t, [], [], systems=[g], name=name)\nend\n\nfunction ConstantVoltage(;name, V = 1.0)\n    val = V\n    @named p = Pin()\n    @named n = Pin()\n    @parameters V\n    eqs = [\n           V ~ p.v - n.v\n           0 ~ p.i + n.i\n          ]\n    ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)\nend\n\nfunction HeatPort(;name)\n    @variables T(t) Q_flow(t)\n    return ODESystem(Equation[], t, [T, Q_flow], [], defaults=Dict(T=>293.15, Q_flow=>0.0), name=name)\nend\n\nfunction HeatingResistor(;name, R=1.0, TAmbient=293.15, alpha=1.0)\n    R_val, TAmbient_val, alpha_val = R, TAmbient, alpha\n    @named p = Pin()\n    @named n = Pin()\n    @named h = HeatPort()\n    @variables v(t) RTherm(t)\n    @parameters R TAmbient alpha\n    eqs = [\n           RTherm ~ R*(1 + alpha*(h.T - TAmbient))\n           v ~ p.i * RTherm\n           h.Q_flow ~ -v * p.i # -LossPower\n           v ~ p.v - n.v\n           0 ~ p.i + n.i\n          ]\n    ODESystem(\n        eqs, t, [v, RTherm], [R, TAmbient, alpha], systems=[p, n, h],\n        defaults=Dict(\n            R=>R_val, TAmbient=>TAmbient_val, alpha=>alpha_val,\n            v=>0.0, RTherm=>R_val\n        ),\n        name=name,\n    )\nend\n\nfunction HeatCapacitor(;name, rho=8050, V=1, cp=460, TAmbient=293.15)\n    rho_val, V_val, cp_val = rho, V, cp\n    @parameters rho V cp\n    C = rho*V*cp\n    @named h = HeatPort()\n    eqs = [\n           D(h.T) ~ h.Q_flow / C\n          ]\n    ODESystem(\n        eqs, t, [], [rho, V, cp], systems=[h],\n        defaults=Dict(rho=>rho_val, V=>V_val, cp=>cp_val),\n        name=name,\n    )\nend\n\nfunction Capacitor(;name, C = 1.0)\n    val = C\n    @named p = Pin()\n    @named n = Pin()\n    @variables v(t)\n    @parameters C\n    eqs = [\n           v ~ p.v - n.v\n           0 ~ p.i + n.i\n           D(v) ~ p.i / C\n          ]\n    ODESystem(\n        eqs, t, [v], [C], systems=[p, n],\n        defaults=Dict(v => 0.0, C => val),\n        name=name\n    )\nend\n\nfunction rc_model(i; name, source, ground, R, C)\n    resistor = HeatingResistor(name=Symbol(:resistor, i), R=R)\n    capacitor = Capacitor(name=Symbol(:capacitor, i), C=C)\n    heat_capacitor = HeatCapacitor(name=Symbol(:heat_capacitor, i))\n\n    rc_eqs = [\n              connect_pin(source.p, resistor.p)\n              connect_pin(resistor.n, capacitor.p)\n              connect_pin(capacitor.n, source.n, ground.g)\n              connect_heat(resistor.h, heat_capacitor.h)\n             ]\n\n    rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground, heat_capacitor], name=Symbol(name, i))\nend","category":"page"},{"location":"tutorials/tearing_parallelism/#The-Model-1","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"The Model","text":"","category":"section"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"Assuming that the components are defined, our model is 50 resistors and capacitors connected in parallel. Thus following the acausal components tutorial, we can connect a bunch of RC components as follows:","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"V = 2.0\nsource = ConstantVoltage(name=:source, V=V)\nground = Ground(name=:ground)\nN = 50\nRs = 10 .^range(0, stop=-4, length=N)\nCs = 10 .^range(-3, stop=0, length=N)\nrc_systems = map(1:N) do i\n    rc_model(i; name=:rc, source=source, ground=ground, R=Rs[i], C=Cs[i])\nend\n@variables E(t)\neqs = [\n       D(E) ~ sum(((i, sys),)->getproperty(sys, Symbol(:resistor, i)).h.Q_flow, enumerate(rc_systems))\n      ]\nbig_rc = ODESystem(eqs, t, [E], [], systems=rc_systems, defaults=Dict(E=>0.0))","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"Now let's say we want to expose a bit more parallelism via running tearing. How do we do that?","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"sys = structural_simplify(big_rc)","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"Done, that's it. There's no more to it.","category":"page"},{"location":"tutorials/tearing_parallelism/#What-Happened?-1","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"What Happened?","text":"","category":"section"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"Yes, that's a good question! Let's investigate a little bit more what had happened. If you look at the system we defined:","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"equations(big_rc)\n\n1051-element Vector{Equation}:\n Differential(t)(E(t)) ~ rc10₊resistor10₊h₊Q_flow(t) + rc11₊resistor11₊h₊Q_flow(t) + rc12₊resistor12₊h₊Q_flow(t) + rc13₊resistor13₊h₊Q_flow(t) + rc14₊resistor14₊h₊Q_flow(t) + rc15₊resistor15₊h₊Q_flow(t) + rc16₊resistor16₊h₊Q_flow(t) + rc17₊resistor17₊h₊Q_flow(t) + rc18₊resistor18₊h₊Q_flow(t) + rc19₊resistor19₊h₊Q_flow(t) + rc1₊resistor1₊h₊Q_flow(t) + rc20₊resistor20₊h₊Q_flow(t) + rc21₊resistor21₊h₊Q_flow(t) + rc22₊resistor22₊h₊Q_flow(t) + rc23₊resistor23₊h₊Q_flow(t) + rc24₊resistor24₊h₊Q_flow(t) + rc25₊resistor25₊h₊Q_flow(t) + rc26₊resistor26₊h₊Q_flow(t) + rc27₊resistor27₊h₊Q_flow(t) + rc28₊resistor28₊h₊Q_flow(t) + rc29₊resistor29₊h₊Q_flow(t) + rc2₊resistor2₊h₊Q_flow(t) + rc30₊resistor30₊h₊Q_flow(t) + rc31₊resistor31₊h₊Q_flow(t) + rc32₊resistor32₊h₊Q_flow(t) + rc33₊resistor33₊h₊Q_flow(t) + rc34₊resistor34₊h₊Q_flow(t) + rc35₊resistor35₊h₊Q_flow(t) + rc36₊resistor36₊h₊Q_flow(t) + rc37₊resistor37₊h₊Q_flow(t) + rc38₊resistor38₊h₊Q_flow(t) + rc39₊resistor39₊h₊Q_flow(t) + rc3₊resistor3₊h₊Q_flow(t) + rc40₊resistor40₊h₊Q_flow(t) + rc41₊resistor41₊h₊Q_flow(t) + rc42₊resistor42₊h₊Q_flow(t) + rc43₊resistor43₊h₊Q_flow(t) + rc44₊resistor44₊h₊Q_flow(t) + rc45₊resistor45₊h₊Q_flow(t) + rc46₊resistor46₊h₊Q_flow(t) + rc47₊resistor47₊h₊Q_flow(t) + rc48₊resistor48₊h₊Q_flow(t) + rc49₊resistor49₊h₊Q_flow(t) + rc4₊resistor4₊h₊Q_flow(t) + rc50₊resistor50₊h₊Q_flow(t) + rc5₊resistor5₊h₊Q_flow(t) + rc6₊resistor6₊h₊Q_flow(t) + rc7₊resistor7₊h₊Q_flow(t) + rc8₊resistor8₊h₊Q_flow(t) + rc9₊resistor9₊h₊Q_flow(t)\n 0 ~ rc1₊resistor1₊p₊i(t) + rc1₊source₊p₊i(t)\n rc1₊source₊p₊v(t) ~ rc1₊resistor1₊p₊v(t)\n 0 ~ rc1₊capacitor1₊p₊i(t) + rc1₊resistor1₊n₊i(t)\n rc1₊resistor1₊n₊v(t) ~ rc1₊capacitor1₊p₊v(t)\n ⋮\n rc50₊source₊V ~ rc50₊source₊p₊v(t) - rc50₊source₊n₊v(t)\n 0 ~ rc50₊source₊n₊i(t) + rc50₊source₊p₊i(t)\n rc50₊ground₊g₊v(t) ~ 0\n Differential(t)(rc50₊heat_capacitor50₊h₊T(t)) ~ rc50₊heat_capacitor50₊h₊Q_flow(t)*(rc50₊heat_capacitor50₊V^-1)*(rc50₊heat_capacitor50₊cp^-1)*(rc50₊heat_capacitor50₊rho^-1)","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"You see it started as a massive 1051 set of equations. However, after eliminating redundancies we arrive at 151 equations:","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"equations(sys)\n\n151-element Vector{Equation}:\n Differential(t)(E(t)) ~ -rc10₊capacitor10₊p₊i(t)*(rc10₊source₊V - rc10₊capacitor10₊v(t)) - (rc11₊capacitor11₊p₊i(t)*(rc11₊source₊V - rc11₊capacitor11₊v(t))) - (rc12₊capacitor12₊p₊i(t)*(rc12₊source₊V - rc12₊capacitor12₊v(t))) - (rc13₊capacitor13₊p₊i(t)*(rc13₊source₊V - rc13₊capacitor13₊v(t))) - (rc14₊capacitor14₊p₊i(t)*(rc14₊source₊V - rc14₊capacitor14₊v(t))) - (rc15₊capacitor15₊p₊i(t)*(rc15₊source₊V - rc15₊capacitor15₊v(t))) - (rc16₊capacitor16₊p₊i(t)*(rc16₊source₊V - rc16₊capacitor16₊v(t))) - (rc17₊capacitor17₊p₊i(t)*(rc17₊source₊V - rc17₊capacitor17₊v(t))) - (rc18₊capacitor18₊p₊i(t)*(rc18₊source₊V - rc18₊capacitor18₊v(t))) - (rc19₊capacitor19₊p₊i(t)*(rc19₊source₊V - rc19₊capacitor19₊v(t))) - (rc1₊resistor1₊p₊i(t)*(rc1₊source₊V - rc1₊capacitor1₊v(t))) - (rc20₊capacitor20₊p₊i(t)*(rc20₊source₊V - rc20₊capacitor20₊v(t))) - (rc21₊capacitor21₊p₊i(t)*(rc21₊source₊V - rc21₊capacitor21₊v(t))) - (rc22₊capacitor22₊p₊i(t)*(rc22₊source₊V - rc22₊capacitor22₊v(t))) - (rc23₊capacitor23₊p₊i(t)*(rc23₊source₊V - rc23₊capacitor23₊v(t))) - (rc24₊capacitor24₊p₊i(t)*(rc24₊source₊V - rc24₊capacitor24₊v(t))) - (rc25₊capacitor25₊p₊i(t)*(rc25₊source₊V - rc25₊capacitor25₊v(t))) - (rc26₊capacitor26₊p₊i(t)*(rc26₊source₊V - rc26₊capacitor26₊v(t))) - (rc27₊capacitor27₊p₊i(t)*(rc27₊source₊V - rc27₊capacitor27₊v(t))) - (rc28₊capacitor28₊p₊i(t)*(rc28₊source₊V - rc28₊capacitor28₊v(t))) - (rc29₊capacitor29₊p₊i(t)*(rc29₊source₊V - rc29₊capacitor29₊v(t))) - (rc2₊capacitor2₊p₊i(t)*(rc2₊source₊V - rc2₊capacitor2₊v(t))) - (rc30₊capacitor30₊p₊i(t)*(rc30₊source₊V - rc30₊capacitor30₊v(t))) - (rc31₊capacitor31₊p₊i(t)*(rc31₊source₊V - rc31₊capacitor31₊v(t))) - (rc32₊capacitor32₊p₊i(t)*(rc32₊source₊V - rc32₊capacitor32₊v(t))) - (rc33₊capacitor33₊p₊i(t)*(rc33₊source₊V - rc33₊capacitor33₊v(t))) - (rc34₊capacitor34₊p₊i(t)*(rc34₊source₊V - rc34₊capacitor34₊v(t))) - (rc35₊capacitor35₊p₊i(t)*(rc35₊source₊V - rc35₊capacitor35₊v(t))) - (rc36₊capacitor36₊p₊i(t)*(rc36₊source₊V - rc36₊capacitor36₊v(t))) - (rc37₊capacitor37₊p₊i(t)*(rc37₊source₊V - rc37₊capacitor37₊v(t))) - (rc38₊capacitor38₊p₊i(t)*(rc38₊source₊V - rc38₊capacitor38₊v(t))) - (rc39₊capacitor39₊p₊i(t)*(rc39₊source₊V - rc39₊capacitor39₊v(t))) - (rc3₊capacitor3₊p₊i(t)*(rc3₊source₊V - rc3₊capacitor3₊v(t))) - (rc40₊capacitor40₊p₊i(t)*(rc40₊source₊V - rc40₊capacitor40₊v(t))) - (rc41₊capacitor41₊p₊i(t)*(rc41₊source₊V - rc41₊capacitor41₊v(t))) - (rc42₊capacitor42₊p₊i(t)*(rc42₊source₊V - rc42₊capacitor42₊v(t))) - (rc43₊capacitor43₊p₊i(t)*(rc43₊source₊V - rc43₊capacitor43₊v(t))) - (rc44₊capacitor44₊p₊i(t)*(rc44₊source₊V - rc44₊capacitor44₊v(t))) - (rc45₊capacitor45₊p₊i(t)*(rc45₊source₊V - rc45₊capacitor45₊v(t))) - (rc46₊capacitor46₊p₊i(t)*(rc46₊source₊V - rc46₊capacitor46₊v(t))) - (rc47₊capacitor47₊p₊i(t)*(rc47₊source₊V - rc47₊capacitor47₊v(t))) - (rc48₊capacitor48₊p₊i(t)*(rc48₊source₊V - rc48₊capacitor48₊v(t))) - (rc49₊capacitor49₊p₊i(t)*(rc49₊source₊V - rc49₊capacitor49₊v(t))) - (rc4₊resistor4₊p₊i(t)*(rc4₊source₊V - rc4₊capacitor4₊v(t))) - (rc50₊capacitor50₊p₊i(t)*(rc50₊source₊V - rc50₊capacitor50₊v(t))) - (rc5₊capacitor5₊p₊i(t)*(rc5₊source₊V - rc5₊capacitor5₊v(t))) - (rc6₊capacitor6₊p₊i(t)*(rc6₊source₊V - rc6₊capacitor6₊v(t))) - (rc7₊capacitor7₊p₊i(t)*(rc7₊source₊V - rc7₊capacitor7₊v(t))) - (rc8₊capacitor8₊p₊i(t)*(rc8₊source₊V - rc8₊capacitor8₊v(t))) - (rc9₊capacitor9₊p₊i(t)*(rc9₊source₊V - rc9₊capacitor9₊v(t)))\n 0 ~ rc1₊capacitor1₊v(t) + rc1₊resistor1₊R*rc1₊resistor1₊p₊i(t)*(1 + rc1₊resistor1₊alpha*(rc1₊heat_capacitor1₊h₊T(t) - rc1₊resistor1₊TAmbient)) - rc1₊source₊V\n Differential(t)(rc1₊capacitor1₊v(t)) ~ rc1₊resistor1₊p₊i(t)*(rc1₊capacitor1₊C^-1)\n Differential(t)(rc1₊heat_capacitor1₊h₊T(t)) ~ rc1₊resistor1₊p₊i(t)*(rc1₊heat_capacitor1₊V^-1)*(rc1₊heat_capacitor1₊cp^-1)*(rc1₊heat_capacitor1₊rho^-1)*(rc1₊source₊V - rc1₊capacitor1₊v(t))\n 0 ~ rc2₊resistor2₊R*rc2₊capacitor2₊p₊i(t)*(1 + rc2₊resistor2₊alpha*(rc2₊heat_capacitor2₊h₊T(t) - rc2₊resistor2₊TAmbient)) + rc2₊capacitor2₊v(t) - rc2₊source₊V\n ⋮\n Differential(t)(rc49₊heat_capacitor49₊h₊T(t)) ~ rc49₊capacitor49₊p₊i(t)*(rc49₊heat_capacitor49₊V^-1)*(rc49₊heat_capacitor49₊cp^-1)*(rc49₊heat_capacitor49₊rho^-1)*(rc49₊source₊V - rc49₊capacitor49₊v(t))\n 0 ~ rc50₊capacitor50₊v(t) + rc50₊resistor50₊R*rc50₊capacitor50₊p₊i(t)*(1 + rc50₊resistor50₊alpha*(rc50₊heat_capacitor50₊h₊T(t) - rc50₊resistor50₊TAmbient)) - rc50₊source₊V\n Differential(t)(rc50₊capacitor50₊v(t)) ~ rc50₊capacitor50₊p₊i(t)*(rc50₊capacitor50₊C^-1)\n Differential(t)(rc50₊heat_capacitor50₊h₊T(t)) ~ rc50₊capacitor50₊p₊i(t)*(rc50₊heat_capacitor50₊V^-1)*(rc50₊heat_capacitor50₊cp^-1)*(rc50₊heat_capacitor50₊rho^-1)*(rc50₊source₊V - rc50₊capacitor50₊v(t))","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"That's not all though. In addition, the tearing process has turned the sets of nonlinear equations into separate blocks and constructed a DAG for the dependencies between the blocks. We can use the bipartite graph functionality to dig in and investigate what this means:","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"using ModelingToolkit.BipartiteGraphs\nbig_rc = initialize_system_structure(big_rc)\ninc_org = BipartiteGraphs.incidence_matrix(structure(big_rc).graph)\nblt_org = StructuralTransformations.sorted_incidence_matrix(big_rc, only_algeqs=true, only_algvars=true)\nblt_reduced = StructuralTransformations.sorted_incidence_matrix(sys, only_algeqs=true, only_algvars=true)","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"(Image: )","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"The figure on the left is the original incidence matrix of the algebraic equations. Notice that the original formulation of the model has dependencies between different equations, and so the full set of equations must be solved together. That exposes no parallelism. However, the Block Lower Triangular (BLT) transformation exposes independent blocks. This is then further impoved by the tearing process, which removes 90% of the equations and transforms the nonlinear equations into 50 independent blocks which can now all be solved in parallel. The conclusion is that, your attempts to parallelize are neigh: performing parallelism after structural simplification greatly improves the problem that can be parallelized, so this is better than trying to do it by hand.","category":"page"},{"location":"tutorials/tearing_parallelism/#","page":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","title":"Exposing More Parallelism By Tearing Algebraic Equations in ODESystems","text":"After performing this, you can construct the ODEProblem/ODAEProblem and set parallel_form to use the exposed parallelism in multithreaded function constructions, but this showcases why structural_simplify is so important to that process.","category":"page"},{"location":"basics/DependencyGraphs/#Dependency-Graphs-1","page":"Dependency Graphs","title":"Dependency Graphs","text":"","category":"section"},{"location":"basics/DependencyGraphs/#Types-1","page":"Dependency Graphs","title":"Types","text":"","category":"section"},{"location":"basics/DependencyGraphs/#","page":"Dependency Graphs","title":"Dependency Graphs","text":"BipartiteGraph","category":"page"},{"location":"basics/DependencyGraphs/#ModelingToolkit.BipartiteGraphs.BipartiteGraph","page":"Dependency Graphs","title":"ModelingToolkit.BipartiteGraphs.BipartiteGraph","text":"mutable struct BipartiteGraph{I<:Integer, F<:Array{Array{I<:Integer, 1}, 1}, B<:Union{Array{Array{I<:Integer, 1}, 1}, I<:Integer}, M} <: LightGraphs.AbstractGraph{I<:Integer}\n\nA bipartite graph representation between two, possibly distinct, sets of vertices (source and dependencies). Maps source vertices, labelled 1:N₁, to vertices on which they depend (labelled 1:N₂).\n\nFields\n\nne\nfadjlist\nbadjlist\nmetadata\n\nExample\n\nusing ModelingToolkit\n\nne = 4\nsrcverts = 1:4\ndepverts = 1:2\n\n# six source vertices\nfadjlist = [[1],[1],[2],[2],[1],[1,2]]\n\n# two vertices they depend on\nbadjlist = [[1,2,5,6],[3,4,6]]\n\nbg = BipartiteGraph(7, fadjlist, badjlist)\n\n\n\n\n\n","category":"type"},{"location":"basics/DependencyGraphs/#Utility-functions-for-BiPartiteGraphs-1","page":"Dependency Graphs","title":"Utility functions for BiPartiteGraphs","text":"","category":"section"},{"location":"basics/DependencyGraphs/#","page":"Dependency Graphs","title":"Dependency Graphs","text":"Base.isequal","category":"page"},{"location":"basics/DependencyGraphs/#Base.isequal","page":"Dependency Graphs","title":"Base.isequal","text":"Base.isequal(bg1::BipartiteGraph{T}, bg2::BipartiteGraph{T}) where {T<:Integer}\n\nTest whether two BipartiteGraphs are equal.\n\n\n\n\n\n","category":"function"},{"location":"basics/DependencyGraphs/#Functions-for-calculating-dependency-graphs-1","page":"Dependency Graphs","title":"Functions for calculating dependency graphs","text":"","category":"section"},{"location":"basics/DependencyGraphs/#","page":"Dependency Graphs","title":"Dependency Graphs","text":"equation_dependencies\nasgraph\nvariable_dependencies\nasdigraph\neqeq_dependencies\nvarvar_dependencies","category":"page"},{"location":"basics/DependencyGraphs/#ModelingToolkit.equation_dependencies","page":"Dependency Graphs","title":"ModelingToolkit.equation_dependencies","text":"equation_dependencies(sys::AbstractSystem; variables=states(sys))\n\nGiven an AbstractSystem calculate for each equation the variables it depends on.\n\nNotes:\n\nVariables that are not in variables are filtered out.\nget_variables! is used to determine the variables within a given equation.\nreturns a Vector{Vector{Variable}}() mapping the index of an equation to the variables it depends on.\n\nExample:\n\nusing ModelingToolkit\n@parameters β γ κ η t\n@variables S(t) I(t) R(t)\n\n# use a reaction system to easily generate ODE and jump systems\nrxs = [Reaction(β, [S,I], [I], [1,1], [2]),\n       Reaction(γ, [I], [R]),\n       Reaction(κ+η, [R], [S])]\nrs = ReactionSystem(rxs, t, [S,I,R], [β,γ,κ,η])\n\n# ODEs:\nodesys = convert(ODESystem, rs)\n\n# dependency of each ODE on state variables\nequation_dependencies(odesys)\n\n# dependency of each ODE on parameters\nequation_dependencies(odesys, variables=parameters(odesys))\n\n# Jumps\njumpsys = convert(JumpSystem, rs)\n\n# dependency of each jump rate function on state variables\nequation_dependencies(jumpsys)\n\n# dependency of each jump rate function on parameters\nequation_dependencies(jumpsys, variables=parameters(jumpsys))\n\n\n\n\n\n","category":"function"},{"location":"basics/DependencyGraphs/#ModelingToolkit.asgraph","page":"Dependency Graphs","title":"ModelingToolkit.asgraph","text":"asgraph(eqdeps, vtois)\n\nConvert a collection of equation dependencies, for example as returned by equation_dependencies, to a BipartiteGraph.\n\nNotes:\n\nvtois should provide a Dict like mapping from each Variable dependency in eqdeps to the integer idx of the variable to use in the graph.\n\nExample: Continuing the example started in equation_dependencies\n\ndigr = asgraph(equation_dependencies(odesys), Dict(s => i for (i,s) in enumerate(states(odesys))))\n\n\n\n\n\nasgraph(sys::AbstractSystem; variables=states(sys),\n                                      variablestoids=Dict(convert(Variable, v) => i for (i,v) in enumerate(variables)))\n\nConvert an AbstractSystem to a BipartiteGraph mapping the index of equations to the indices of variables they depend on.\n\nNotes:\n\nDefaults for kwargs creating a mapping from equations(sys) to states(sys) they depend on.\nvariables should provide the list of variables to use for generating the dependency graph.\nvariablestoids should provide Dict like mapping from a Variable to its Int index within variables.\n\nExample: Continuing the example started in equation_dependencies\n\ndigr = asgraph(odesys)\n\n\n\n\n\n","category":"function"},{"location":"basics/DependencyGraphs/#ModelingToolkit.variable_dependencies","page":"Dependency Graphs","title":"ModelingToolkit.variable_dependencies","text":"variable_dependencies(sys::AbstractSystem; variables=states(sys), variablestoids=nothing)\n\nFor each variable determine the equations that modify it and return as a BipartiteGraph.\n\nNotes:\n\nDependencies are returned as a BipartiteGraph mapping variable indices to the indices of equations that modify them.\nvariables denotes the list of variables to determine dependencies for.\nvariablestoids denotes a Dict mapping Variables to their Int index in variables.\n\nExample: Continuing the example of equation_dependencies\n\nvariable_dependencies(odesys)\n\n\n\n\n\n","category":"function"},{"location":"basics/DependencyGraphs/#ModelingToolkit.asdigraph","page":"Dependency Graphs","title":"ModelingToolkit.asdigraph","text":"asdigraph(g::BipartiteGraph, sys::AbstractSystem; variables = states(sys), equationsfirst = true)\n\nConvert a BipartiteGraph to a LightGraph.SimpleDiGraph.\n\nNotes:\n\nThe resulting SimpleDiGraph unifies the two sets of vertices (equations and then states in the case it comes from asgraph), producing one ordered set of integer vertices (SimpleDiGraph does not support two distinct collections of vertices so they must be merged).\nvariables gives the variables that g is associated with (usually the states of a system).\nequationsfirst (default is true) gives whether the BipartiteGraph gives a mapping from equations to variables they depend on (true), as calculated by asgraph, or whether it gives a mapping from variables to the equations that modify them, as calculated by variable_dependencies.\n\nExample: Continuing the example in asgraph\n\ndg = asdigraph(digr)\n\n\n\n\n\n","category":"function"},{"location":"basics/DependencyGraphs/#ModelingToolkit.eqeq_dependencies","page":"Dependency Graphs","title":"ModelingToolkit.eqeq_dependencies","text":"eqeq_dependencies(eqdeps::BipartiteGraph{T}, vardeps::BipartiteGraph{T}) where {T <: Integer}\n\nCalculate a LightGraph.SimpleDiGraph that maps each equation to equations they depend on.\n\nNotes:\n\nThe fadjlist of the SimpleDiGraph maps from an equation to the equations that modify variables it depends on.\nThe badjlist of the SimpleDiGraph maps from an equation to equations that depend on variables it modifies.\n\nExample: Continuing the example of equation_dependencies\n\neqeqdep = eqeq_dependencies(asgraph(odesys), variable_dependencies(odesys))\n\n\n\n\n\n","category":"function"},{"location":"basics/DependencyGraphs/#ModelingToolkit.varvar_dependencies","page":"Dependency Graphs","title":"ModelingToolkit.varvar_dependencies","text":"varvar_dependencies(eqdeps::BipartiteGraph{T}, vardeps::BipartiteGraph{T}) where {T <: Integer} = eqeq_dependencies(vardeps, eqdeps)\n\nCalculate a LightGraph.SimpleDiGraph that maps each variable to variables they depend on.\n\nNotes:\n\nThe fadjlist of the SimpleDiGraph maps from a variable to the variables that depend on it.\nThe badjlist of the SimpleDiGraph maps from a variable to variables on which it depends.\n\nExample: Continuing the example of equation_dependencies\n\nvarvardep = varvar_dependencies(asgraph(odesys), variable_dependencies(odesys))\n\n\n\n\n\n","category":"function"},{"location":"tutorials/nonlinear/#Modeling-Nonlinear-Systems-1","page":"Modeling Nonlinear Systems","title":"Modeling Nonlinear Systems","text":"","category":"section"},{"location":"tutorials/nonlinear/#","page":"Modeling Nonlinear Systems","title":"Modeling Nonlinear Systems","text":"In this example we will go one step deeper and showcase the direct function generation capabilities in ModelingToolkit.jl to build nonlinear systems. Let's say we wanted to solve for the steady state of the previous ODE. This is the nonlinear system defined by where the derivatives are zero. We use (unknown) variables for our nonlinear system.","category":"page"},{"location":"tutorials/nonlinear/#","page":"Modeling Nonlinear Systems","title":"Modeling Nonlinear Systems","text":"using ModelingToolkit, NonlinearSolve\r\n\r\n@variables x y z\r\n@parameters σ ρ β\r\n\r\n# Define a nonlinear system\r\neqs = [0 ~ σ*(y-x),\r\n       0 ~ x*(ρ-z)-y,\r\n       0 ~ x*y - β*z]\r\nns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])\r\n\r\nguess = [x => 1.0,\r\n         y => 0.0,\r\n         z => 0.0]\r\n\r\nps = [\r\n      σ => 10.0\r\n      ρ => 26.0\r\n      β => 8/3\r\n      ]\r\n\r\nprob = NonlinearProblem(ns,guess,ps)\r\nsol = solve(prob,NewtonRaphson())","category":"page"},{"location":"tutorials/nonlinear/#","page":"Modeling Nonlinear Systems","title":"Modeling Nonlinear Systems","text":"We can similarly ask to generate the NonlinearProblem with the analytical Jacobian function:","category":"page"},{"location":"tutorials/nonlinear/#","page":"Modeling Nonlinear Systems","title":"Modeling Nonlinear Systems","text":"prob = NonlinearProblem(ns,guess,ps,jac=true)\r\nsol = solve(prob,NewtonRaphson())","category":"page"},{"location":"systems/JumpSystem/#JumpSystem-1","page":"JumpSystem","title":"JumpSystem","text":"","category":"section"},{"location":"systems/JumpSystem/#System-Constructors-1","page":"JumpSystem","title":"System Constructors","text":"","category":"section"},{"location":"systems/JumpSystem/#","page":"JumpSystem","title":"JumpSystem","text":"JumpSystem","category":"page"},{"location":"systems/JumpSystem/#ModelingToolkit.JumpSystem","page":"JumpSystem","title":"ModelingToolkit.JumpSystem","text":"struct JumpSystem{U<:RecursiveArrayTools.ArrayPartition} <: ModelingToolkit.AbstractSystem\n\nA system of jump processes.\n\nFields\n\neqs\nThe jumps of the system. Allowable types are ConstantRateJump, VariableRateJump, MassActionJump.\n\niv\nThe independent variable, usually time.\nstates\nThe dependent variables, representing the state of the system.\nps\nThe parameters of the system.\nobserved\nname\nThe name of the system. . These are required to have unique names.\nsystems\nThe internal systems.\ndefaults\ndefaults: The default values to use when initial conditions and/or parameters are not supplied in ODEProblem.\n\nconnection_type\ntype: type of the system\n\nExample\n\nusing ModelingToolkit\n\n@parameters β γ t\n@variables S I R\nrate₁   = β*S*I\naffect₁ = [S ~ S - 1, I ~ I + 1]\nrate₂   = γ*I\naffect₂ = [I ~ I - 1, R ~ R + 1]\nj₁      = ConstantRateJump(rate₁,affect₁)\nj₂      = ConstantRateJump(rate₂,affect₂)\nj₃      = MassActionJump(2*β+γ, [R => 1], [S => 1, R => -1])\njs      = JumpSystem([j₁,j₂,j₃], t, [S,I,R], [β,γ])\n\n\n\n\n\n","category":"type"},{"location":"systems/JumpSystem/#Composition-and-Accessor-Functions-1","page":"JumpSystem","title":"Composition and Accessor Functions","text":"","category":"section"},{"location":"systems/JumpSystem/#","page":"JumpSystem","title":"JumpSystem","text":"get_eqs(sys) or equations(sys): The equations that define the jump system.\nget_states(sys) or states(sys): The set of states in the jump system.\nget_ps(sys) or parameters(sys): The parameters of the jump system.\nindependent_variable(sys): The independent variable of the jump system.","category":"page"},{"location":"systems/JumpSystem/#Transformations-1","page":"JumpSystem","title":"Transformations","text":"","category":"section"},{"location":"systems/JumpSystem/#","page":"JumpSystem","title":"JumpSystem","text":"structural_simplify","category":"page"},{"location":"systems/JumpSystem/#Analyses-1","page":"JumpSystem","title":"Analyses","text":"","category":"section"},{"location":"systems/JumpSystem/#Problem-Constructors-1","page":"JumpSystem","title":"Problem Constructors","text":"","category":"section"},{"location":"systems/JumpSystem/#","page":"JumpSystem","title":"JumpSystem","text":"DiscreteProblem\nJumpProblem","category":"page"},{"location":"systems/JumpSystem/#SciMLBase.DiscreteProblem","page":"JumpSystem","title":"SciMLBase.DiscreteProblem","text":"function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan,\n                                    parammap=DiffEqBase.NullParameters; kwargs...)\n\nGenerates a blank DiscreteProblem for a pure jump JumpSystem to utilize as its prob.prob. This is used in the case where there are no ODEs and no SDEs associated with the system.\n\nContinuing the example from the JumpSystem definition:\n\nusing DiffEqBase, DiffEqJump\nu₀map = [S => 999, I => 1, R => 0]\nparammap = [β => .1/1000, γ => .01]\ntspan = (0.0, 250.0)\ndprob = DiscreteProblem(js, u₀map, tspan, parammap)\n\n\n\n\n\nDiscreteProblem(sys::DiscreteSystem, u0map::Any, tspan::Any) -> Any\nDiscreteProblem(sys::DiscreteSystem, u0map::Any, tspan::Any, parammap::Any; eval_module, eval_expression, kwargs...) -> Any\n\n\nGenerates an DiscreteProblem from an DiscreteSystem.\n\n\n\n\n\n","category":"type"},{"location":"systems/JumpSystem/#DiffEqJump.JumpProblem","page":"JumpSystem","title":"DiffEqJump.JumpProblem","text":"function DiffEqBase.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)\n\nGenerates a JumpProblem from a JumpSystem.\n\nContinuing the example from the DiscreteProblem definition:\n\njprob = JumpProblem(js, dprob, Direct())\nsol = solve(jprob, SSAStepper())\n\n\n\n\n\n","category":"type"},{"location":"basics/ContextualVariables/#Contextual-Variable-Types-1","page":"Contextual Variable Types","title":"Contextual Variable Types","text":"","category":"section"},{"location":"basics/ContextualVariables/#","page":"Contextual Variable Types","title":"Contextual Variable Types","text":"ModelingToolkit.jl has a system of contextual variable types which allows for helping the system transformation machinery do complex manipulations and automatic detection. The standard variable definition in ModelingToolkit.jl is the @variable which is defined by Symbolics.jl. For example:","category":"page"},{"location":"basics/ContextualVariables/#","page":"Contextual Variable Types","title":"Contextual Variable Types","text":"@variables x y(x)","category":"page"},{"location":"basics/ContextualVariables/#","page":"Contextual Variable Types","title":"Contextual Variable Types","text":"This is used for the \"normal\" variable of a given system, like the states of a differential equation or objective function. All of the macros below support the same syntax as @variables.","category":"page"},{"location":"basics/ContextualVariables/#Parameters-1","page":"Contextual Variable Types","title":"Parameters","text":"","category":"section"},{"location":"basics/ContextualVariables/#","page":"Contextual Variable Types","title":"Contextual Variable Types","text":"All modeling projects have some form of parameters. @parameters marks a variable as being the parameter of some system, which allows automatic detection algorithms to ignore such variables when attempting to find the states of a system.","category":"page"},{"location":"basics/ContextualVariables/#Flow-Variables-(TODO)-1","page":"Contextual Variable Types","title":"Flow Variables (TODO)","text":"","category":"section"},{"location":"basics/ContextualVariables/#","page":"Contextual Variable Types","title":"Contextual Variable Types","text":"In many engineering systems some variables act like \"flows\" while others do not. For example, in circuit models you have current which flows, and the related voltage which does not. Or in thermal models you have heat flows. In these cases, the connect statement enforces conservation of flow between all of the connected components.","category":"page"},{"location":"basics/ContextualVariables/#","page":"Contextual Variable Types","title":"Contextual Variable Types","text":"For example, the following specifies that x is a 2x2 matrix of flow variables:","category":"page"},{"location":"basics/ContextualVariables/#","page":"Contextual Variable Types","title":"Contextual Variable Types","text":"@variables x[1:2,1:2] type=flow","category":"page"},{"location":"basics/ContextualVariables/#Stream-Variables-1","page":"Contextual Variable Types","title":"Stream Variables","text":"","category":"section"},{"location":"basics/ContextualVariables/#","page":"Contextual Variable Types","title":"Contextual Variable Types","text":"TODO","category":"page"},{"location":"basics/ContextualVariables/#Brownian-Variables-1","page":"Contextual Variable Types","title":"Brownian Variables","text":"","category":"section"},{"location":"basics/ContextualVariables/#","page":"Contextual Variable Types","title":"Contextual Variable Types","text":"TODO","category":"page"},{"location":"#ModelingToolkit.jl:-High-Performance-Symbolic-Numeric-Equation-Based-Modeling-1","page":"Home","title":"ModelingToolkit.jl: High-Performance Symbolic-Numeric Equation-Based Modeling","text":"","category":"section"},{"location":"#","page":"Home","title":"Home","text":"ModelingToolkit.jl is a modeling language for high-performance symbolic-numeric computation in scientific computing and scientific machine learning. It then mixes ideas from symbolic computational algebra systems with causal and acausal equation-based modeling frameworks to give an extendable and parallel modeling system. It allows for users to give a high-level description of a model for symbolic preprocessing to analyze and enhance the model. Automatic transformations, such as index reduction, can be applied to the model before solving in order to make it easily handle equations would could not be solved when modeled without symbolic intervention.","category":"page"},{"location":"#Installation-1","page":"Home","title":"Installation","text":"","category":"section"},{"location":"#","page":"Home","title":"Home","text":"To install ModelingToolkit.jl, use the Julia package manager:","category":"page"},{"location":"#","page":"Home","title":"Home","text":"using Pkg\nPkg.add(\"ModelingToolkit\")","category":"page"},{"location":"#Citation-1","page":"Home","title":"Citation","text":"","category":"section"},{"location":"#","page":"Home","title":"Home","text":"If you use ModelingToolkit in your work, please cite the following:","category":"page"},{"location":"#","page":"Home","title":"Home","text":"@misc{ma2021modelingtoolkit,\n      title={ModelingToolkit: A Composable Graph Transformation System For Equation-Based Modeling},\n      author={Yingbo Ma and Shashi Gowda and Ranjan Anantharaman and Chris Laughman and Viral Shah and Chris Rackauckas},\n      year={2021},\n      eprint={2103.05244},\n      archivePrefix={arXiv},\n      primaryClass={cs.MS}\n}","category":"page"},{"location":"#Feature-Summary-1","page":"Home","title":"Feature Summary","text":"","category":"section"},{"location":"#","page":"Home","title":"Home","text":"ModelingToolkit.jl is a symbolic-numeric modeling package. Thus it combines some of the features from symbolic computing packages like SymPy or Mathematica with the ideas of equation-based modeling systems like the causal Simulink and the acausal Modelica. It bridges the gap between many different kinds of equations, allowing one to quickly and easily transform systems of DAEs into optimization problems, or vice-versa, and then simplify and parallelize the resulting expressions before generating code.","category":"page"},{"location":"#Feature-List-1","page":"Home","title":"Feature List","text":"","category":"section"},{"location":"#","page":"Home","title":"Home","text":"Causal and acausal modeling (Simulink/Modelica)\nAutomated model transformation, simplification, and composition\nAutomatic conversion of numerical models into symbolic models\nComposition of models through the components, a lazy connection system, and tools for expanding/flattening\nPervasive parallelism in symbolic computations and generated functions\nTransformations like alias elimination and tearing of nonlinear systems for efficiently numerically handling large-scale systems of equations\nThe ability to use the entire Symbolics.jl Computer Algebra System (CAS) as part of the modeling process.\nImport models from common formats like SBML, CellML, BioNetGen, and more.\nExtendability: the whole system is written in pure Julia, so adding new functions, simplification rules, and model transformations has no barrier.","category":"page"},{"location":"#","page":"Home","title":"Home","text":"For information on how to use the Symbolics.jl CAS system that ModelingToolkit.jl is built on, consult the Symbolics.jl documentation","category":"page"},{"location":"#Equation-Types-1","page":"Home","title":"Equation Types","text":"","category":"section"},{"location":"#","page":"Home","title":"Home","text":"Ordinary differential equations\nStochastic differential equations\nPartial differential equations\nNonlinear systems\nOptimization problems\nContinuous-Time Markov Chains\nChemical Reactions\nNonlinear Optimal Control","category":"page"},{"location":"#Model-Import-Formats-1","page":"Home","title":"Model Import Formats","text":"","category":"section"},{"location":"#","page":"Home","title":"Home","text":"CellMLToolkit.jl: Import CellML models into ModelingToolkit\nRepository of more than a thousand pre-made models\nFocus on biomedical models in areas such as: Calcium Dynamics, Cardiovascular Circulation, Cell Cycle, Cell Migration, Circadian Rhythms, Electrophysiology, Endocrine, Excitation-Contraction Coupling, Gene Regulation, Hepatology, Immunology, Ion Transport, Mechanical Constitutive Laws, Metabolism, Myofilament Mechanics, Neurobiology, pH Regulation, PKPD, Protein Modules, Signal Transduction, and Synthetic Biology.\nSbmlInterface.jl: Import SBML models into ModelingToolkit\nUses the robust libsbml library for parsing and transforming the SBML\nReactionNetworkImporters.jl: Import various models into ModelingToolkit\nSupports the BioNetGen .net file\nSupports importing networks specified by stoichiometric matrices","category":"page"},{"location":"#Extension-Libraries-1","page":"Home","title":"Extension Libraries","text":"","category":"section"},{"location":"#","page":"Home","title":"Home","text":"Because ModelingToolkit.jl is the core foundation of a equation-based modeling ecosystem, there is a large set of libraries adding features to this system. Below is an incomplete list of extension libraries one may want to be aware of:","category":"page"},{"location":"#","page":"Home","title":"Home","text":"Catalyst.jl: Symbolic representations of chemical reactions\nSymbolically build and represent large systems of chemical reactions\nGenerate code for ODEs, SDEs, continuous-time Markov Chains, and more\nSimulate the models using the SciML ecosystem with O(1) Gillespie methods\nDataDrivenDiffEq.jl: Automatic identification of equations from data\nAutomated construction of ODEs and DAEs from data\nRepresentations of Koopman operators and Dynamic Mode Decomposition (DMD)\nMomentClosure.jl: Automatic transformation of ReactionSystems into deterministic systems\nGenerates ODESystems for the moment closures\nAllows for geometrically-distributed random reaction rates\nReactionMechanismSimulator.jl: simulating and analyzing large chemical reaction mechanisms\nIdeal gas and dilute liquid phases.\nConstant T and P and constant V adiabatic ideal gas reactors.\nConstant T and V dilute liquid reactors.\nDiffusion limited rates. Sensitivity analysis for all reactors.\nFlux diagrams with molecular images (if molecular information is provided).","category":"page"},{"location":"#Compatible-Numerical-Solvers-1","page":"Home","title":"Compatible Numerical Solvers","text":"","category":"section"},{"location":"#","page":"Home","title":"Home","text":"All of the symbolic systems have a direct conversion to a numerical system which can then be handled through the SciML interfaces. For example, after building a model and performing symbolic manipulations, an ODESystem can be converted into an ODEProblem to then be solved by a numerical ODE solver. Below is a list of the solver libraries which are the numerical targets of the ModelingToolkit system:","category":"page"},{"location":"#","page":"Home","title":"Home","text":"DifferentialEquations.jl\nMulti-package interface of high performance numerical solvers for ODESystem, SDESystem, and JumpSystem\nNonlinearSolve.jl\nHigh performance numerical solving of NonlinearSystem\nGalacticOptim.jl\nMulti-package interface for numerical solving OptimizationSystem\nNeuralPDE.jl\nPhysics-Informed Neural Network (PINN) training on PDESystem\nDiffEqOperators.jl\nAutomated finite difference method (FDM) discretization of PDESystem","category":"page"},{"location":"#Contributing-1","page":"Home","title":"Contributing","text":"","category":"section"},{"location":"#","page":"Home","title":"Home","text":"Please refer to the SciML ColPrac: Contributor's Guide on Collaborative Practices for Community Packages for guidance on PRs, issues, and other matters relating to contributing to ModelingToolkit.\nThere are a few community forums:\nThe #diffeq-bridged channel in the Julia Slack\nJuliaDiffEq on Gitter\nOn the Julia Discourse forums (look for the modelingtoolkit tag\nSee also SciML Community page","category":"page"}]
}