Building and solving numerical problems

ModelingToolkit.SemilinearODEFunction(sys::System; kwargs...)
ModelingToolkit.SemilinearODEFunction{iip}(sys::System; kwargs...)
ModelingToolkit.SemilinearODEFunction{iip, specialize}(sys::System; kwargs...)

Create a ModelingToolkit.SemilinearODEFunction from the given sys. iip is a boolean indicating whether the function should be in-place. specialization is a SciMLBase.AbstractSpecalize subtype indicating the level of specialization of the ModelingToolkit.SemilinearODEFunction.

This is a special form of an ODE which uses a SplitFunction internally. The equations are separated into linear, quadratic and general terms and phrased as matrix operations. See calculate_semiquadratic_form for information on how the equations are split. This formulation allows leveraging split ODE solvers such as KenCarp4 and is useful for systems where the stiff and non-stiff terms can be separated out in such a manner. Typically the linear part of the equations is the stiff part, but the keywords stiff_linear, stiff_quadratic and stiff_nonlinear can be used to control which parts are considered as stiff.

Beyond the arguments listed below, this constructor accepts all keyword arguments supported by the DifferentialEquations.jl solve function. For a complete list and detailed descriptions, see the DifferentialEquations.jl solve documentation.

Keyword arguments

• u0: The u0 vector for the corresponding problem, if available. Can be obtained using ModelingToolkitBase.get_u0.

• p: The parameter object for the corresponding problem, if available. Can be obtained using ModelingToolkitBase.get_p.

• t: The initial time for the corresponding problem, if available.

• eval_expression: Whether to compile any functions via eval or RuntimeGeneratedFunctions.

• eval_module: If eval_expression == true, the module to eval into. Otherwise, the module in which to generate the RuntimeGeneratedFunction.

• checkbounds: Whether to enable bounds checking in the generated code.

• simplify: Whether to simplify any symbolically computed jacobians/hessians/etc.

• cse: Whether to enable Common Subexpression Elimination (CSE) on the generated code. This typically improves performance of the generated code but reduces readability.

• sparse: Whether to generate jacobian/hessian/etc. functions that return/operate on sparse matrices. Also controls whether the mass matrix is sparse, wherever applicable.

• check_compatibility: Whether to check if the given system sys contains all the information necessary to create a ModelingToolkit.SemilinearODEFunction and no more. If disabled, assumes that sys at least contains the necessary information.

• expression: Val{true} to return an Expr that constructs the corresponding problem instead of the problem itself. Val{false} otherwise. Constructing the expression does not support callbacks

• jac: Whether to symbolically compute and generate code for the jacobian function.

• sparsity: Whether to provide symbolically compute and provide sparsity patterns for the jacobian/hessian/etc.

• stiff_linear: Whether the linear part of the equations should be part of the stiff function in the split form. Has no effect if the equations have no linear part.

• stiff_quadratic: Whether the quadratic part of the equations should be part of the stiff function in the split form. Has no effect if the equations have no quadratic part.

• stiff_nonlinear: Whether the non-linear non-quadratic part of the equations should be part of the stiff function in the split form. Has no effect if the equations have no such non-linear non-quadratic part.

• kwargs...: Additional keyword arguments passed to the solver

All other keyword arguments are forwarded to the ModelingToolkit.SemilinearODEFunction struct constructor. Note that all three of stiff_linear, stiff_quadratic, stiff_nonlinear cannot be identical, and at least two of A, B, C returned from calculate_semiquadratic_form must be non-nothing. In other words, both of the functions in the split form must be non-empty.

SciMLBase.ModelingToolkit.SemilinearODEProblem(sys::System, op, tspan::NTuple{2}; kwargs...)
SciMLBase.ModelingToolkit.SemilinearODEProblem{iip}(sys::System, op, tspan::NTuple{2}; kwargs...)
SciMLBase.ModelingToolkit.SemilinearODEProblem{iip, specialize}(sys::System, op, tspan::NTuple{2}; kwargs...)

Build a ModelingToolkit.SemilinearODEProblem given a system sys and operating point op and timespan tspan. iip is a boolean indicating whether the problem should be in-place. specialization is a SciMLBase.AbstractSpecalize subtype indicating the level of specialization of the ModelingToolkit.SemilinearODEFunction. The operating point should be an iterable collection of key-value pairs mapping variables/parameters in the system to the (initial) values they should take in ModelingToolkit.SemilinearODEProblem. Any values not provided will fallback to the corresponding default (if present).

ModelingToolkitBase will build an initialization problem where all initial values for unknowns or observables of sys (either explicitly provided or in defaults) will be constraints. To remove an initial condition in the defaults (without providing a replacement) give the corresponding variable a value of nothing in the operating point. The initialization problem will also run parameter initialization. See the Initialization documentation for more information.

This is a special form of an ODE which uses a SplitFunction internally. The equations are separated into linear, quadratic and general terms and phrased as matrix operations. See calculate_semiquadratic_form for information on how the equations are split. This formulation allows leveraging split ODE solvers such as KenCarp4 and is useful for systems where the stiff and non-stiff terms can be separated out in such a manner. Typically the linear part of the equations is the stiff part, but the keywords stiff_linear, stiff_quadratic and stiff_nonlinear can be used to control which parts are considered as stiff.

Keyword arguments

• eval_expression: Whether to compile any functions via eval or RuntimeGeneratedFunctions.

• eval_module: If eval_expression == true, the module to eval into. Otherwise, the module in which to generate the RuntimeGeneratedFunction.

• guesses: The guesses for variables in the system, used as initial values for the initialization problem.

• warn_initialize_determined: Warn if the initialization system is under/over-determined.

• initialization_eqs: Extra equations to use in the initialization problem.

• fully_determined: Override whether the initialization system is fully determined.

• use_scc: Whether to use SCCNonlinearProblem for initialization if the system is fully determined.

• check_initialization_units: Enable or disable unit checks when constructing the initialization problem.

• tofloat: Passed to varmap_to_vars when building the parameter vector of a non-split system.

• u0_eltype: The eltype of the u0 vector. If nothing, finds the promoted floating point type from op.

• u0_constructor: A function to apply to the u0 value returned from varmap_to_vars. to construct the final u0 value.

• p_constructor: A function to apply to each array buffer created when constructing the parameter object.

• warn_cyclic_dependency: Whether to emit a warning listing out cycles in initial conditions provided for unknowns and parameters.

• circular_dependency_max_cycle_length: Maximum length of cycle to check for. Only applicable if warn_cyclic_dependency == true.

• circular_dependency_max_cycles: Maximum number of cycles to check for. Only applicable if warn_cyclic_dependency == true.

• substitution_limit: The number times to substitute initial conditions into each other to attempt to arrive at a numeric value.

• missing_guess_value: An instance of MissingGuessValue which indicates what happens when the initialization problem is missing guess values for variables.

• initsys_mtkcompile_kwargs: A NamedTuple of keyword arguments to pass to mtkcompile when it is called on the initialization system.

• callback: An extra callback or CallbackSet to add to the problem, in addition to the ones defined symbolically in the system.

• check_compatibility: Whether to check if the given system sys contains all the information necessary to create a ModelingToolkit.SemilinearODEProblem and no more. If disabled, assumes that sys at least contains the necessary information.

• expression: Val{true} to return an Expr that constructs the corresponding problem instead of the problem itself. Val{false} otherwise. Constructing the expression does not support callbacks

• stiff_linear: Whether the linear part of the equations should be part of the stiff function in the split form. Has no effect if the equations have no linear part.

• stiff_quadratic: Whether the quadratic part of the equations should be part of the stiff function in the split form. Has no effect if the equations have no quadratic part.

• stiff_nonlinear: Whether the non-linear non-quadratic part of the equations should be part of the stiff function in the split form. Has no effect if the equations have no such non-linear non-quadratic part.

All other keyword arguments are forwarded to the ModelingToolkit.SemilinearODEFunction constructor. Note that all three of stiff_linear, stiff_quadratic, stiff_nonlinear cannot be identical, and at least two of A, B, C returned from calculate_semiquadratic_form must be non-nothing. In other words, both of the functions in the split form must be non-empty.

Extended docs

The following API is internal and may change or be removed without notice. Its usage is highly discouraged.

• build_initializeprob: If false, avoids building the initialization problem.
• check_length: Whether to check the number of equations along with number of unknowns and length of u0 vector for consistency. If false, do not check with equations. This is forwarded to check_eqs_u0.
• return_operating_point: Whether to also return the updated operating point.
• time_dependent_init: Whether to build a time-dependent initialization for the problem. A time-dependent initialization solves for a consistent u0, whereas a time-independent one only runs parameter initialization.
• algebraic_only: Whether to build the initialization problem using only algebraic equations.
• allow_incomplete: Whether to allow incomplete initialization problems.

SciMLBase.SciMLBase.SCCNonlinearProblem(sys::System, op; kwargs...)
SciMLBase.SciMLBase.SCCNonlinearProblem{iip}(sys::System, op; kwargs...)
SciMLBase.SciMLBase.SCCNonlinearProblem{iip, specialize}(sys::System, op; kwargs...)

Build a SciMLBase.SCCNonlinearProblem given a system sys and operating point op . iip is a boolean indicating whether the problem should be in-place. specialization is a SciMLBase.AbstractSpecalize subtype indicating the level of specialization of the SciMLBase.NonlinearFunction. The operating point should be an iterable collection of key-value pairs mapping variables/parameters in the system to the (initial) values they should take in SciMLBase.SCCNonlinearProblem. Any values not provided will fallback to the corresponding default (if present).

Keyword arguments

• eval_expression: Whether to compile any functions via eval or RuntimeGeneratedFunctions.

• eval_module: If eval_expression == true, the module to eval into. Otherwise, the module in which to generate the RuntimeGeneratedFunction.

• guesses: The guesses for variables in the system, used as initial values for the initialization problem.

• warn_initialize_determined: Warn if the initialization system is under/over-determined.

• initialization_eqs: Extra equations to use in the initialization problem.

• fully_determined: Override whether the initialization system is fully determined.

• use_scc: Whether to use SCCNonlinearProblem for initialization if the system is fully determined.

• check_initialization_units: Enable or disable unit checks when constructing the initialization problem.

• tofloat: Passed to varmap_to_vars when building the parameter vector of a non-split system.

• u0_eltype: The eltype of the u0 vector. If nothing, finds the promoted floating point type from op.

• u0_constructor: A function to apply to the u0 value returned from varmap_to_vars. to construct the final u0 value.

• p_constructor: A function to apply to each array buffer created when constructing the parameter object.

• warn_cyclic_dependency: Whether to emit a warning listing out cycles in initial conditions provided for unknowns and parameters.

• circular_dependency_max_cycle_length: Maximum length of cycle to check for. Only applicable if warn_cyclic_dependency == true.

• circular_dependency_max_cycles: Maximum number of cycles to check for. Only applicable if warn_cyclic_dependency == true.

• substitution_limit: The number times to substitute initial conditions into each other to attempt to arrive at a numeric value.

• missing_guess_value: An instance of MissingGuessValue which indicates what happens when the initialization problem is missing guess values for variables.

• initsys_mtkcompile_kwargs: A NamedTuple of keyword arguments to pass to mtkcompile when it is called on the initialization system.

• check_compatibility: Whether to check if the given system sys contains all the information necessary to create a SciMLBase.SCCNonlinearProblem and no more. If disabled, assumes that sys at least contains the necessary information.
• expression: Val{true} to return an Expr that constructs the corresponding problem instead of the problem itself. Val{false} otherwise.

All other keyword arguments are forwarded to the SciMLBase.NonlinearFunction constructor.

Extended docs

The following API is internal and may change or be removed without notice. Its usage is highly discouraged.

• build_initializeprob: If false, avoids building the initialization problem.
• check_length: Whether to check the number of equations along with number of unknowns and length of u0 vector for consistency. If false, do not check with equations. This is forwarded to check_eqs_u0.
• return_operating_point: Whether to also return the updated operating point.
• time_dependent_init: Whether to build a time-dependent initialization for the problem. A time-dependent initialization solves for a consistent u0, whereas a time-independent one only runs parameter initialization.
• algebraic_only: Whether to build the initialization problem using only algebraic equations.
• allow_incomplete: Whether to allow incomplete initialization problems.

lin_fun, simplified_sys = linearization_function(sys::AbstractSystem, inputs, outputs; simplify = false, initialize = true, initialization_solver_alg = nothing, kwargs...)

Return a function that linearizes the system sys. The function linearize provides a higher-level and easier to use interface.

lin_fun is a function (variables, p, t) -> (; f_x, f_z, g_x, g_z, f_u, g_u, h_x, h_z, h_u), i.e., it returns a NamedTuple with the Jacobians of f,g,h for the nonlinear sys (technically for simplified_sys) on the form

\[\begin{aligned} ẋ &= f(x, z, u) \\ 0 &= g(x, z, u) \\ y &= h(x, z, u) \end{aligned}\]

where x are differential unknown variables, z algebraic variables, u inputs and y outputs. To obtain a linear statespace representation, see linearize. The input argument variables is a vector defining the operating point, corresponding to unknowns(simplified_sys) and p is a vector corresponding to the parameters of simplified_sys. Note: all variables in inputs have been converted to parameters in simplified_sys.

The simplified_sys has undergone mtkcompile and had any occurring input or output variables replaced with the variables provided in arguments inputs and outputs. The unknowns of this system also indicate the order of the unknowns that holds for the linearized matrices.

Arguments:

• sys: A System of ODEs. This function will automatically apply simplification passes on sys and return the resulting simplified_sys.
• inputs: A vector of variables that indicate the inputs of the linearized input-output model.
• outputs: A vector of variables that indicate the outputs of the linearized input-output model.
• simplify: Apply simplification in tearing.
• initialize: If true, a check is performed to ensure that the operating point is consistent (satisfies algebraic equations). If the op is not consistent, initialization is performed.
• initialization_solver_alg: A NonlinearSolve algorithm to use for solving for a feasible set of state and algebraic variables that satisfies the specified operating point.
• autodiff: An ADType supported by DifferentiationInterface.jl to use for calculating the necessary jacobians. Defaults to using AutoForwardDiff()
• kwargs: Are passed on to find_solvables!

See also linearize which provides a higher-level interface.

mutable struct LinearizationProblem{F<:ModelingToolkit.LinearizationFunction, T}

A struct representing a linearization operation to be performed. Can be symbolically indexed to efficiently update the operating point for multiple linearizations in a loop. The value of the independent variable can be set by mutating the .t field of this struct.

(; A, B, C, D), simplified_sys, extras = linearize(sys, inputs, outputs;    t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false, kwargs...)
(; A, B, C, D), extras                 = linearize(simplified_sys, lin_fun; t=0.0, op = Dict(), allow_input_derivatives = false, zero_dummy_der=false)

Linearize sys between inputs and outputs, both vectors of variables. Return a NamedTuple with the matrices of a linear statespace representation on the form

\[\begin{aligned} ẋ &= Ax + Bu\\ y &= Cx + Du \end{aligned}\]

The first signature automatically calls linearization_function internally, while the second signature expects the outputs of linearization_function as input.

op denotes the operating point around which to linearize. If none is provided, the default values of sys are used.

If allow_input_derivatives = false, an error will be thrown if input derivatives ($u̇$) appear as inputs in the linearized equations. If input derivatives are allowed, the returned B matrix will be of double width, corresponding to the input [u; u̇].

zero_dummy_der can be set to automatically set the operating point to zero for all dummy derivatives.

The return value extras is a NamedTuple (; x, p, t) containing the result of the initialization problem that was solved to determine the operating point.

See also linearization_function which provides a lower-level interface, linearize_symbolic and ModelingToolkit.reorder_unknowns.

See extended help for an example.

The implementation and notation follows that of "Linear Analysis Approach for Modelica Models", Allain et al. 2009

Extended help

This example builds the following feedback interconnection and linearizes it from the input of F to the output of P.

  r ┌─────┐       ┌─────┐     ┌─────┐
───►│     ├──────►│     │  u  │     │
    │  F  │       │  C  ├────►│  P  │ y
    └─────┘     ┌►│     │     │     ├─┬─►
                │ └─────┘     └─────┘ │
                │                     │
                └─────────────────────┘

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
function plant(; name)
    @variables x(t) = 1
    @variables u(t)=0 y(t)=0
    eqs = [D(x) ~ -x + u
           y ~ x]
    System(eqs, t; name = name)
end

function ref_filt(; name)
    @variables x(t)=0 y(t)=0
    @variables u(t)=0 [input = true]
    eqs = [D(x) ~ -2 * x + u
           y ~ x]
    System(eqs, t, name = name)
end

function controller(kp; name)
    @variables y(t)=0 r(t)=0 u(t)=0
    @parameters kp = kp
    eqs = [
        u ~ kp * (r - y),
    ]
    System(eqs, t; name = name)
end

@named f = ref_filt()
@named c = controller(1)
@named p = plant()

connections = [f.y ~ c.r # filtered reference to controller reference
               c.u ~ p.u # controller output to plant input
               p.y ~ c.y]

@named cl = System(connections, t, systems = [f, c, p])

lsys0, ssys = linearize(cl, [f.u], [p.x])
desired_order = [f.x, p.x]
lsys = ModelingToolkit.reorder_unknowns(lsys0, unknowns(ssys), desired_order)

@assert lsys.A == [-2 0; 1 -2]
@assert lsys.B == [1; 0;;]
@assert lsys.C == [0 1]
@assert lsys.D[] == 0

## Symbolic linearization
lsys_sym, _ = ModelingToolkit.linearize_symbolic(cl, [f.u], [p.x])

@assert substitute(lsys_sym.A, ModelingToolkit.defaults(cl)) == lsys.A

(; Ã, B̃, C̃, D̃) = similarity_transform(sys, T; unitary=false)

Perform a similarity transform T : Tx̃ = x on linear system represented by matrices in NamedTuple sys such that

à = T⁻¹AT
B̃ = T⁻¹ B
C̃ = CT
D̃ = D

If unitary=true, T is assumed unitary and the matrix adjoint is used instead of the inverse.

reorder_unknowns(sys::NamedTuple, old, new)

Permute the state representation of sys obtained from linearize so that the state unknown is changed from old to new Example:

lsys, ssys = linearize(pid, [reference.u, measurement.u], [ctr_output.u])
desired_order = [int.x, der.x] # Unknowns that are present in unknowns(ssys)
lsys = ModelingToolkit.reorder_unknowns(lsys, unknowns(ssys), desired_order)

See also ModelingToolkit.similarity_transform