Disturbance and input modeling modeling

Disturbances are often seen as external factors that influence a system. Modeling and simulation of such external influences is common in order to ensure that the plant and or control system can adequately handle or suppress these disturbances. Disturbance modeling is also integral to the problem of state estimation, indeed, modeling how disturbances affect the evolution of the state of the system is crucial in order to accurately estimate this state.

This tutorial will show how to model disturbances in ModelingToolkit as disturbance inputs. This involves demonstrating best practices that make it easy to use a single model to handle both disturbed and undisturbed systems, and making use of the model for both simulation and state estimation.

A flexible component-based workflow

We will consider a simple system consisting of two inertias connected through a flexible shaft, such as a simple transmission system in a car. We start by modeling the plant without any input signals:

using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, Test
using ModelingToolkitStandardLibrary.Mechanical.Rotational
using ModelingToolkitStandardLibrary.Blocks
t = ModelingToolkit.t_nounits
D = ModelingToolkit.D_nounits

@mtkmodel SystemModel begin
    @parameters begin
        m1 = 1
        m2 = 1
        k = 10 # Spring stiffness
        c = 3  # Damping coefficient
    end
    @components begin
        inertia1 = Inertia(; J = m1, phi = 0, w = 0)
        inertia2 = Inertia(; J = m2, phi = 0, w = 0)
        spring = Spring(; c = k)
        damper = Damper(; d = c)
        torque = Torque(use_support = false)
    end
    @equations begin
        connect(torque.flange, inertia1.flange_a)
        connect(inertia1.flange_b, spring.flange_a, damper.flange_a)
        connect(inertia2.flange_a, spring.flange_b, damper.flange_b)
    end
end
<< @example-block not executed in draft mode >>

Here, we have added a torque component that allows us to add a torque input to drive the system, but we have not connected any signal to it yet. We have not yet made any attempts at modeling disturbances, and this is deliberate, we will handle this later in order to make the plant model as generically useful as possible.

In order to simulate this system in the presence of disturbances, we must 1. Reason about how disturbances may affect the system, and 2. attach disturbance inputs and disturbance signals to the model. We distinguish between an input and a signal here, where we by input mean an attachment point (connector) to which we may connect a signal, i.e., a time-varying function.

We create a new model that includes disturbance inputs and signals, and attach those to the already defined plant model. We assume that each of the two inertias can be affected by a disturbance torque, such as due to friction or an unknown load on the output inertia.

@mtkmodel ModelWithInputs begin
    @components begin
        input_signal = Blocks.Sine(frequency = 1, amplitude = 1)
        disturbance_signal1 = Blocks.Step(height = -1, start_time = 2) # We add an input signal that equals zero by default so that it has no effect during normal simulation
        disturbance_signal2 = Blocks.Step(height = 2, start_time = 4)
        disturbance_torque1 = Torque(use_support = false)
        disturbance_torque2 = Torque(use_support = false)
        system_model = SystemModel()
    end
    @equations begin
        connect(input_signal.output, :u, system_model.torque.tau)
        connect(disturbance_signal1.output, :d1, disturbance_torque1.tau) # When we connect the input _signals_, we do so through an analysis point. This allows us to easily disconnect the input signals in situations when we do not need them. 
        connect(disturbance_signal2.output, :d2, disturbance_torque2.tau)
        connect(disturbance_torque1.flange, system_model.inertia1.flange_b)
        connect(disturbance_torque2.flange, system_model.inertia2.flange_b)
    end
end
<< @example-block not executed in draft mode >>

This outer model, ModelWithInputs, contains two disturbance inputs, both of type Torque. It also contains three signal specifications, one for the control input and two for the corresponding disturbance inputs. Note how we added the disturbance torque inputs at this level of the model, but the control input was added inside the system model. This is a design choice that is up to the modeler, here, we consider the driving torque to be a fundamental part of the model that is always required to make use of it, while the disturbance inputs may be of interest only in certain situations, and we thus add them when needed. Since we have added not only input connectors, but also connected input signals to them, this model is complete and ready for simulation, i.e., there are no unbound inputs.

@named model = ModelWithInputs() # Model with load disturbance
ssys = mtkcompile(model)
prob = ODEProblem(ssys, [], (0.0, 6.0))
sol = solve(prob, Tsit5())
using Plots
plot(sol)
<< @example-block not executed in draft mode >>

A thing to note in the specification of ModelWithInputs is the presence of three analysis points, :u, :d1, and :d2. When signals are connected through an analysis point, we may at any time linearize the model as if the signals were not connected, i.e., as if the corresponding inputs were unbound. We may also use this to generate a julia function for the dynamics on the form $f(x,u,p,t,w)$ where the input $u$ and disturbance $w$ may be provided as separate function arguments, as if the corresponding input signals were not present in the model. More details regarding this will be presented further below, here, we just demonstrate how we could linearize this system model from the inputs to the angular velocity of the inertias

using ControlSystemsBase, ControlSystemsMTK # ControlSystemsMTK provides the high-level function named_ss and ControlSystemsBase provides the bodeplot function
P = named_ss(model, [ssys.u, ssys.d1, ssys.d2],
    [ssys.system_model.inertia1.w, ssys.system_model.inertia2.w])
bodeplot(P, plotphase = false)
<< @example-block not executed in draft mode >>

It's worth noting at this point that the fact that we could connect disturbance outputs from outside of the plant-model definition was enabled by the fact that we used a component-based workflow, where the plant model had the appropriate connectors available. If the plant model had modeled the system using direct equations without connectors, this would not have been possible and the model would thus be significantly less flexible.

We summarize the findings so far as a number of best practices:

Best practices
  • Use a component-based workflow to model the plant
  • If possible, model the plant without explicit disturbance inputs to make it as generic as possible
  • When disturbance inputs are needed, create a new model that includes the plant model and the disturbance inputs
  • Only add input signals at the top level of the model, this applies to both control inputs and disturbance inputs.
  • Use analysis points to connect signals to inputs, this allows for easy disconnection of signals when needed, e.g., for linearization or function generation.

Modeling for state estimation

In the example above, we constructed a model for simulation of a disturbance affecting the system. When simulating, we connect an input signal of specified shape that simulates the disturbance, above, we used Blocks.Step as input signals. On the other hand, when performing state estimation, the exact shape of the disturbance is typically not known, we might only have some diffuse knowledge of the disturbance characteristics such as "varies smoothly", "makes sudden step changes" or "is approximately periodic with 24hr period". The encoding of such knowledge is commonly reasoned about in the frequency domain, where we specify a disturbance model as a dynamical system with a frequency response similar to the approximate spectrum of the disturbance. For more details around this, see the in-depth tutorial notebook "How to tune a Kalman filter". Most algorithms for state estimation, such as a Kalman-filter like estimators, assume that disturbances are independent and identically distributed (i.i.d.). While seemingly restrictive at first glance, when combined with an appropriate disturbance models encoded as dynamical systems, this assumption still allows for a wide range of non i.i.d. disturbances to be modeled.

When modeling a system in MTK, we essentially (without considering algebraic equations for simplicity in exposition) construct a model of a dynamical system

\[\begin{aligned} \dot x &= f(x, p, t) \\ y &= g(x, p, t) \end{aligned}\]

where $x$ is the state, $y$ are observed variables, $p$ are parameters, and $t$ is time. When using MTK, which variables constitute $x$ and which are considered part of the output, $y$, is up to the tool rather than the user, this choice is made by MTK during the call to @mtkcompile or the lower-level function mtkcompile.

If we further consider external inputs to the system, such as controlled input signals $u$ and disturbance inputs $w$, we can write the system as

\[\begin{aligned} \dot x &= f(x, u, p, t, w) \\ y &= g(x, u, p, t) \end{aligned}\]

To make use of the model defined above for state estimation, we may want to generate a Julia function for the dynamics $f$ and the output equations $g$ that we can plug into, e.g., a nonlinear version of a Kalman filter or a particle filter, etc. MTK contains utilities to do this, namely, ModelingToolkit.generate_control_function and ModelingToolkit.build_explicit_observed_function (described in more details in "Input output"). These functions take keyword arguments disturbance_inputs and disturbance_argument, that indicate which variables in the model are considered part of $w$, and whether or not these variables are to be added as function arguments to $f$, i.e., whether we have $f(x, u, p, t)$ or $f(x, u, p, t, w)$. If we do not include the disturbance inputs as function arguments, MTK will assume that the $w$ variables are all zero, but any dynamics associated with these variables, such as disturbance models, will be included in the generated function. This allows a state estimator to estimate the state of the disturbance model, provided that this state is observable from the measured outputs of the system.

Below, we demonstrate

  1. How to add an integrating disturbance model
  2. how to generate the functions $f$ and $g$ for a typical nonlinear state estimator with explicit disturbance inputs
@mtkmodel IntegratingDisturbance begin
    @variables begin
        x(t) = 0.0
        w(t) = 0.0, [disturbance = true, input = true]
    end
    @components begin
        input = RealInput()
        output = RealOutput()
    end
    @equations begin
        D(x) ~ w
        w ~ input.u
        output.u ~ x
    end
end

@mtkmodel SystemModelWithDisturbanceModel begin
    @components begin
        input_signal = Blocks.Sine(frequency = 1, amplitude = 1)
        disturbance_signal1 = Blocks.Constant(k = 0)
        disturbance_signal2 = Blocks.Constant(k = 0)
        disturbance_torque1 = Torque(use_support = false)
        disturbance_torque2 = Torque(use_support = false)
        disturbance_model = Blocks.Integrator()
        system_model = SystemModel()
    end
    @equations begin
        connect(input_signal.output, :u, system_model.torque.tau)
        connect(disturbance_signal1.output, :d1, disturbance_model.input)
        connect(disturbance_model.output, disturbance_torque1.tau)
        connect(disturbance_signal2.output, :d2, disturbance_torque2.tau)
        connect(disturbance_torque1.flange, system_model.inertia1.flange_b)
        connect(disturbance_torque2.flange, system_model.inertia2.flange_b)
    end
end

@named model_with_disturbance = SystemModelWithDisturbanceModel()
<< @example-block not executed in draft mode >>

We demonstrate that this model is complete and can be simulated:

ssys = mtkcompile(model_with_disturbance)
prob = ODEProblem(ssys, [], (0.0, 10.0))
sol = solve(prob, Tsit5())
using Test
@test SciMLBase.successful_retcode(sol)
<< @example-block not executed in draft mode >>

but we may also generate the functions $f$ and $g$ for state estimation:

inputs = [ssys.u]
disturbance_inputs = [ssys.d1, ssys.d2]
P = ssys.system_model
outputs = [P.inertia1.phi, P.inertia2.phi, P.inertia1.w, P.inertia2.w]

(f_oop, f_ip), x_sym,
p_sym,
io_sys = ModelingToolkit.generate_control_function(
    model_with_disturbance, inputs, disturbance_inputs; disturbance_argument = true)

g = ModelingToolkit.build_explicit_observed_function(
    io_sys, outputs; inputs)

op = ModelingToolkit.inputs(io_sys) .=> 0
x0 = ModelingToolkit.get_u0(io_sys, op)
p = MTKParameters(io_sys, op)
u = zeros(1) # Control input
w = zeros(length(disturbance_inputs)) # Disturbance input
@test f_oop(x0, u, p, t, w) == zeros(5)
@test g(x0, u, p, 0.0) == [0, 0, 0, 0]

# Non-zero disturbance inputs should result in non-zero state derivatives. We call `sort` since we do not generally know the order of the state variables
w = [1.0, 2.0]
@test sort(f_oop(x0, u, p, t, w)) == [0, 0, 0, 1, 2]
<< @example-block not executed in draft mode >>

Input signal library

The Blocks module in ModelingToolkitStandardLibrary contains several predefined input signals, such as Sine, Step, Ramp, Constant etc., a few of which were used in the examples above. If you have an input signal represented as a sequence of samples, you may use an Interpolation block, e.g., as src = Interpolation(ConstantInterpolation, data, timepoints), see the docstring for a complete example.

Disturbance-model library

There is no library explicitly constructed for disturbance modeling. Standard blocks from the Blocks module in ModelingToolkitStandardLibrary, such as Integrator, TransferFunction, StateSpace, can model any disturbance with rational spectrum. Examples of this includes disturbance models such as constants, piecewise constant, periodic, highpass, lowpass, and bandpass. For help with filter design, see ControlSystems.jl: Filter-design and the interface package ControlSystemsMTK.jl. In the example above, we made use of Blocks.Integrator, which is a disturbance model suitable for disturbances dominated by low-frequency components, such as piecewise constant signals or slowly drifting signals.

Further reading

To see full examples that perform state estimation with ModelingToolkit models, see the following resources:

Index

    ModelingToolkit.connectMethod
    connect(output_connector, ap_name::Symbol, input_connector; verbose = true)
    connect(output_connector, ap::AnalysisPoint, input_connector; verbose = true)

    Connect output_connector and input_connector with an AnalysisPoint inbetween. The incoming connection output_connector is expected to be an output connector (for example, ModelingToolkitStandardLibrary.Blocks.RealOutput), and vice versa.

    PLEASE NOTE: The connection is assumed to be causal, meaning that

    @named P = FirstOrder(k = 1, T = 1)
    @named C = Gain(; k = -1)
    connect(C.output, :plant_input, P.input)

    is correct, whereas

    connect(P.input, :plant_input, C.output)

    typically is not (unless the model is an inverse model).

    Arguments

    • output_connector: An output connector
    • input_connector: An input connector
    • ap: An explicitly created AnalysisPoint
    • ap_name: If a name is given, an AnalysisPoint with the given name will be created automatically.

    Keyword arguments

    • verbose: Warn if an input is connected to an output (reverse causality). Silence this warning if you are analyzing an inverse model.
    source
    ModelingToolkit.get_comp_sensitivityFunction
    get_comp_sensitivity(sys, ap::AnalysisPoint; kwargs)
    get_comp_sensitivity(sys, ap_name::Symbol; kwargs)

    Compute the complementary sensitivity function in analysis point ap. The complementary sensitivity function is obtained by introducing an infinitesimal perturbation d at the output of ap, linearizing the system and computing the transfer function between d and the input of ap.

    Arguments:

    • kwargs: Are sent to ModelingToolkit.linearize

    See also get_sensitivity, get_looptransfer.

    source
    ModelingToolkit.get_comp_sensitivity_functionMethod
    get_comp_sensitivity_function(
        sys::ModelingToolkit.AbstractSystem,
        aps;
        kwargs...
    ) -> Tuple{ModelingToolkit.LinearizationFunction{DI, AI, _A, P, _B, _C, J1, J2, J3, J4, IA, @NamedTuple{abstol::Float64, reltol::Float64, nlsolve_alg::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{Missing, NonlinearSolveFirstOrder.GenericTrustRegionScheme{NonlinearSolveFirstOrder.RadiusUpdateSchemes.__Simple, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, NonlinearSolveBase.Dogleg{NonlinearSolveBase.NewtonDescent{Nothing}, NonlinearSolveBase.SteepestDescent{Nothing}}, Nothing, Nothing, Nothing, Val{false}}}} where {DI<:AbstractVector{Int64}, AI<:AbstractVector{Int64}, _A, P<:ODEProblem, _B, _C, J1<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing}}, ModelingToolkit.var"#uff#953"{var"#1488#fun"}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Any}} where {T<:(ForwardDiff.Tag{F} where F<:ModelingToolkit.var"#uff#953"), _A, _B}), var"#1488#fun", _A}}, J2<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing}}, _A, _B, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:ForwardDiff.JacobianConfig, _A, _B}}, J3<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#pff#954"{var"#1489#fun", var"#1490#setter"}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Any}} where {T<:(ForwardDiff.Tag{F} where F<:ModelingToolkit.var"#pff#954"), _A, _B}), var"#1489#fun", var"#1490#setter", _A}}, J4<:(ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#hpf#952"{fun, setter}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Any}} where {T<:(ForwardDiff.Tag{F} where F<:ModelingToolkit.var"#hpf#952"), _A, _B}), fun, setter, _A}), IA<:Union{SciMLBase.NoInit, SciMLBase.OverrideInit{Nothing, Nothing, Nothing}}}, Any}
    

    Return the complementary sensitivity function for the analysis point(s) aps, and the modified system simplified with the appropriate inputs and outputs.

    Keyword Arguments

    - `loop_openings`: A list of analysis points whose connections should be removed and
      the outputs set to the input as a part of the linear analysis.
    
    - `system_modifier`: A function taking the transformed system and applying any
      additional transformations, returning the modified system. The modified system
      is passed to `linearization_function`.

    All other keyword arguments are forwarded to linearization_function.

    source
    ModelingToolkit.get_looptransferFunction
    get_looptransfer(sys, ap::AnalysisPoint; kwargs)
    get_looptransfer(sys, ap_name::Symbol; kwargs)

    Compute the (linearized) loop-transfer function in analysis point ap, from ap.out to ap.in.

    Negative feedback

    Feedback loops often use negative feedback, and the computed loop-transfer function will in this case have the negative feedback included. Standard analysis tools often assume a loop-transfer function without the negative gain built in, and the result of this function may thus need negation before use.

    Arguments:

    • kwargs: Are sent to ModelingToolkit.linearize

    See also get_sensitivity, get_comp_sensitivity, open_loop.

    source
    ModelingToolkit.get_looptransfer_functionMethod
    get_looptransfer_function(
        sys::ModelingToolkit.AbstractSystem,
        aps;
        kwargs...
    ) -> Tuple{ModelingToolkit.LinearizationFunction{DI, AI, _A, P, _B, _C, J1, J2, J3, J4, IA, @NamedTuple{abstol::Float64, reltol::Float64, nlsolve_alg::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{Missing, NonlinearSolveFirstOrder.GenericTrustRegionScheme{NonlinearSolveFirstOrder.RadiusUpdateSchemes.__Simple, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, NonlinearSolveBase.Dogleg{NonlinearSolveBase.NewtonDescent{Nothing}, NonlinearSolveBase.SteepestDescent{Nothing}}, Nothing, Nothing, Nothing, Val{false}}}} where {DI<:AbstractVector{Int64}, AI<:AbstractVector{Int64}, _A, P<:ODEProblem, _B, _C, J1<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing}}, ModelingToolkit.var"#uff#953"{var"#1488#fun"}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Any}} where {T<:(ForwardDiff.Tag{F} where F<:ModelingToolkit.var"#uff#953"), _A, _B}), var"#1488#fun", _A}}, J2<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing}}, _A, _B, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:ForwardDiff.JacobianConfig, _A, _B}}, J3<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#pff#954"{var"#1489#fun", var"#1490#setter"}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Any}} where {T<:(ForwardDiff.Tag{F} where F<:ModelingToolkit.var"#pff#954"), _A, _B}), var"#1489#fun", var"#1490#setter", _A}}, J4<:(ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#hpf#952"{fun, setter}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Any}} where {T<:(ForwardDiff.Tag{F} where F<:ModelingToolkit.var"#hpf#952"), _A, _B}), fun, setter, _A}), IA<:Union{SciMLBase.NoInit, SciMLBase.OverrideInit{Nothing, Nothing, Nothing}}}, Any}
    

    Return the loop-transfer function for the analysis point(s) aps, and the modified system simplified with the appropriate inputs and outputs.

    Keyword Arguments

    - `loop_openings`: A list of analysis points whose connections should be removed and
      the outputs set to the input as a part of the linear analysis.
    
    - `system_modifier`: A function taking the transformed system and applying any
      additional transformations, returning the modified system. The modified system
      is passed to `linearization_function`.

    All other keyword arguments are forwarded to linearization_function.

    source
    ModelingToolkit.get_sensitivityFunction
        get_sensitivity(sys, ap::AnalysisPoint; kwargs)
        get_sensitivity(sys, ap_name::Symbol; kwargs)

    Compute the sensitivity function in analysis point ap. The sensitivity function is obtained by introducing an infinitesimal perturbation d at the input of ap, linearizing the system and computing the transfer function between d and the output of ap.

    Arguments:

    • kwargs: Are sent to ModelingToolkit.linearize

    See also get_comp_sensitivity, get_looptransfer.

    source
    ModelingToolkit.get_sensitivity_functionMethod
    get_sensitivity_function(
        sys::ModelingToolkit.AbstractSystem,
        aps;
        kwargs...
    ) -> Tuple{ModelingToolkit.LinearizationFunction{DI, AI, _A, P, _B, _C, J1, J2, J3, J4, IA, @NamedTuple{abstol::Float64, reltol::Float64, nlsolve_alg::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithm{Missing, NonlinearSolveFirstOrder.GenericTrustRegionScheme{NonlinearSolveFirstOrder.RadiusUpdateSchemes.__Simple, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, NonlinearSolveBase.Dogleg{NonlinearSolveBase.NewtonDescent{Nothing}, NonlinearSolveBase.SteepestDescent{Nothing}}, Nothing, Nothing, Nothing, Val{false}}}} where {DI<:AbstractVector{Int64}, AI<:AbstractVector{Int64}, _A, P<:ODEProblem, _B, _C, J1<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing}}, ModelingToolkit.var"#uff#953"{var"#1488#fun"}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Any}} where {T<:(ForwardDiff.Tag{F} where F<:ModelingToolkit.var"#uff#953"), _A, _B}), var"#1488#fun", _A}}, J2<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing}}, _A, _B, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:ForwardDiff.JacobianConfig, _A, _B}}, J3<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#pff#954"{var"#1489#fun", var"#1490#setter"}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Any}} where {T<:(ForwardDiff.Tag{F} where F<:ModelingToolkit.var"#pff#954"), _A, _B}), var"#1489#fun", var"#1490#setter", _A}}, J4<:(ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#hpf#952"{fun, setter}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Any}} where {T<:(ForwardDiff.Tag{F} where F<:ModelingToolkit.var"#hpf#952"), _A, _B}), fun, setter, _A}), IA<:Union{SciMLBase.NoInit, SciMLBase.OverrideInit{Nothing, Nothing, Nothing}}}, Any}
    

    Return the sensitivity function for the analysis point(s) aps, and the modified system simplified with the appropriate inputs and outputs.

    Keyword Arguments

    - `loop_openings`: A list of analysis points whose connections should be removed and
      the outputs set to the input as a part of the linear analysis.
    
    - `system_modifier`: A function taking the transformed system and applying any
      additional transformations, returning the modified system. The modified system
      is passed to `linearization_function`.

    All other keyword arguments are forwarded to linearization_function.

    source
    ModelingToolkit.linearization_ap_transformMethod
    sys, input_vars, output_vars =
    linearization_ap_transform(
        sys,
        inputs::Union{Vector{Symbol}, Vector{AnalysisPoint}, Symbol, AnalysisPoint},
        outputs,
        loop_openings
    ) -> Tuple{Any, Vector{Any}, Vector{Any}}
    

    Apply analysis-point transformations to prepare a system for linearization.

    Returns

    • sys: The transformed system.
    • input_vars: A vector of input variables corresponding to the input analysis points.
    • output_vars: A vector of output variables corresponding to the output analysis points.
    source
    ModelingToolkit.open_loopMethod
    open_loop(
        sys,
        ap::Union{Symbol, AnalysisPoint};
        system_modifier
    ) -> Tuple{Any, Tuple{Any, Any}}
    

    Apply LoopTransferTransform to the analysis point ap and return the result of apply_transformation.

    Keyword Arguments

    • system_modifier: a function which takes the modified system and returns a new system with any required further modifications performed.
    source
    ModelingToolkit.AnalysisPointType
    struct AnalysisPoint
    AnalysisPoint(input, name::Symbol, outputs::Vector)

    Create an AnalysisPoint for linear analysis. Analysis points can be created by calling

    connect(out, :ap_name, in...)

    Where out is the output being connected to the inputs in.... All involved connectors (input and outputs) are required to either have an unknown named u or a single unknown, all of which should have the same size.

    See also get_sensitivity, get_comp_sensitivity, get_looptransfer, open_loop

    Fields

    • input::Any: The input to the connection. In the context of ModelingToolkitStandardLibrary.jl, this is a RealOutput connector.
    • name::Symbol: The name of the analysis point.
    • outputs::Union{Nothing, Vector{Any}}: The outputs of the connection. In the context of ModelingToolkitStandardLibrary.jl, these are all RealInput connectors.

    Example

    using ModelingToolkit
    using ModelingToolkitStandardLibrary.Blocks
    using ModelingToolkit: t_nounits as t
    
    @named P = FirstOrder(k = 1, T = 1)
    @named C = Gain(; k = -1)
    t = ModelingToolkit.get_iv(P)
    
    eqs = [connect(P.output, C.input)
           connect(C.output, :plant_input, P.input)]
    sys = System(eqs, t, systems = [P, C], name = :feedback_system)
    
    matrices_S, _ = get_sensitivity(sys, :plant_input) # Compute the matrices of a state-space representation of the (input) sensitivity function.
    matrices_T, _ = get_comp_sensitivity(sys, :plant_input)

    Continued linear analysis and design can be performed using ControlSystemsBase.jl. Create ControlSystemsBase.StateSpace objects using

    using ControlSystemsBase, Plots
    S = ss(matrices_S...)
    T = ss(matrices_T...)
    bodeplot([S, T], lab = ["S" "T"])

    The sensitivity functions obtained this way should be equivalent to the ones obtained with the code below

    using ControlSystemsBase
    P = tf(1.0, [1, 1])
    C = 1                      # Negative feedback assumed in ControlSystems
    S = sensitivity(P, C)      # or feedback(1, P*C)
    T = comp_sensitivity(P, C) # or feedback(P*C)
    source