Linear Analysis

Linear analysis refers to the process of linearizing a nonlinear model and analysing the resulting linear dynamical system. To facilitate linear analysis, ModelingToolkit provides the concept of an AnalysisPoint, which can be inserted in-between two causal blocks (such as those from ModelingToolkitStandardLibrary.Blocks sub module). Once a model containing analysis points is built, several operations are available:

get_sensitivity get the sensitivity function (wiki), $S(s)$, as defined in the field of control theory.
get_comp_sensitivity get the complementary sensitivity function $T(s) : S(s)+T(s)=1$.
get_looptransfer get the (open) loop-transfer function where the loop starts and ends in the analysis point. For a typical simple feedback connection with a plant $P(s)$ and a controller $C(s)$, the loop-transfer function at the plant output is $P(s)C(s)$.
linearize can be called with two analysis points denoting the input and output of the linearized system.
open_loop return a new (nonlinear) system where the loop has been broken in the analysis point, i.e., the connection the analysis point usually implies has been removed.

An analysis point can be created explicitly using the constructor AnalysisPoint, or automatically when connecting two causal components using connect:

connect(comp1.output, :analysis_point_name, comp2.input)

A single output can also be connected to multiple inputs:

connect(comp1.output, :analysis_point_name, comp2.input, comp3.input, comp4.input)

Causality

Analysis points are causal, i.e., they imply a directionality for the flow of information. The order of the connections in the connect statement is thus important, i.e., connect(out, :name, in) is different from connect(in, :name, out).

The directionality of an analysis point can be thought of as an arrow in a block diagram, where the name of the analysis point applies to the arrow itself.

┌─────┐         ┌─────┐
│     │  name   │     │
│  out├────────►│in   │
│     │         │     │
└─────┘         └─────┘

This is signified by the name being the middle argument to connect.

Of the above mentioned functions, all except for open_loop return the output of ModelingToolkit.linearize, which is

matrices, simplified_sys = linearize(...)
# matrices = (; A, B, C, D)

i.e., matrices is a named tuple containing the matrices of a linear state-space system on the form

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

Example

The following example builds a simple closed-loop system with a plant $P$ and a controller $C$. Two analysis points are inserted, one before and one after $P$. We then derive a number of sensitivity functions and show the corresponding code using the package ControlSystemBase.jl

using ModelingToolkitStandardLibrary.Blocks, ModelingToolkit
@named P = FirstOrder(k = 1, T = 1) # A first-order system with pole in -1
@named C = Gain(-1)             # A P controller
t = ModelingToolkit.get_iv(P)
eqs = [connect(P.output, :plant_output, C.input)  # Connect with an automatically created analysis point called :plant_output
       connect(C.output, :plant_input, P.input)]
sys = ODESystem(eqs, t, systems = [P, C], name = :feedback_system)

matrices_S = get_sensitivity(sys, :plant_input)[1] # Compute the matrices of a state-space representation of the (input)sensitivity function.
matrices_T = get_comp_sensitivity(sys, :plant_input)[1]

(A = [-2.0;;], B = [-1.0;;], C = [-1.0;;], D = [0.0;;])

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

using ControlSystemsBase, Plots
S = ss(matrices_S...)
T = ss(matrices_T...)
bodeplot([S, T], lab = ["S" "" "T" ""], plot_title = "Bode plot of sensitivity functions",
    margin = 5Plots.mm)

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]) |> ss
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)

ControlSystemsBase.StateSpace{ControlSystemsBase.Continuous, Float64}
A = 
 -2.0
B = 
 1.0
C = 
 1.0
D = 
 0.0

Continuous-time state-space model

We may also derive the loop-transfer function $L(s) = P(s)C(s)$ using

matrices_L = get_looptransfer(sys, :plant_output)[1]
L = ss(matrices_L...)

ControlSystemsBase.StateSpace{ControlSystemsBase.Continuous, Float64}
A = 
 -1.0
B = 
 -1.0
C = 
 1.0
D = 
 0.0

Continuous-time state-space model

which is equivalent to the following with ControlSystems

L = P * (-C) # Add the minus sign to build the negative feedback into the controller

ControlSystemsBase.StateSpace{ControlSystemsBase.Continuous, Float64}
A = 
 -1.0
B = 
 -1.0
C = 
 1.0
D = 
 -0.0

Continuous-time state-space model

To obtain the transfer function between two analysis points, we call linearize

matrices_PS = linearize(sys, :plant_input, :plant_output)[1]

(A = [-2.0;;], B = [1.0;;], C = [1.0;;], D = [0.0;;])

this particular transfer function should be equivalent to the linear system P(s)S(s), i.e., equivalent to

feedback(P, C)

ControlSystemsBase.StateSpace{ControlSystemsBase.Continuous, Float64}
A = 
 -2.0
B = 
 1.0
C = 
 1.0
D = 
 0.0

Continuous-time state-space model

Obtaining transfer functions

A statespace system from ControlSystemsBase can be converted to a transfer function using the function tf:

tf(S)

ControlSystemsBase.TransferFunction{ControlSystemsBase.Continuous, ControlSystemsBase.SisoRational{Float64}}
1.0s + 1.0
----------
1.0s + 2.0

Continuous-time transfer function model

Gain and phase margins

Further linear analysis can be performed using the analysis methods from ControlSystemsBase. For example, calculating the gain and phase margins of a system can be done using

margin(P)

(wgm = [NaN;;], gm = [Inf;;], wpm = [NaN;;], pm = [Inf;;])

(they are infinite for this system). A Nyquist plot can be produced using

nyquistplot(P)

Index

ModelingToolkit.AnalysisPoint
ModelingToolkit.AnalysisPoint
ModelingToolkit.get_comp_sensitivity
ModelingToolkit.get_comp_sensitivity_function
ModelingToolkit.get_looptransfer
ModelingToolkit.get_looptransfer_function
ModelingToolkit.get_sensitivity
ModelingToolkit.get_sensitivity_function
ModelingToolkit.open_loop
Symbolics.connect

ModelingToolkit.get_comp_sensitivity — Function

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

source

ModelingToolkit.get_comp_sensitivity_function — Method

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#816"{var"#1475#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#816"), _A, _B}), var"#1475#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#817"{var"#1476#fun", var"#1477#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#817"), _A, _B}), var"#1476#fun", var"#1477#setter", _A}}, J4<:(ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#hpf#815"{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#815"), _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 zero 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_looptransfer — Function

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

source

ModelingToolkit.get_looptransfer_function — Method

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#816"{var"#1475#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#816"), _A, _B}), var"#1475#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#817"{var"#1476#fun", var"#1477#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#817"), _A, _B}), var"#1476#fun", var"#1477#setter", _A}}, J4<:(ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#hpf#815"{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#815"), _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 zero 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_sensitivity — Function

    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

source

ModelingToolkit.get_sensitivity_function — Method

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#816"{var"#1475#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#816"), _A, _B}), var"#1475#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#817"{var"#1476#fun", var"#1477#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#817"), _A, _B}), var"#1476#fun", var"#1477#setter", _A}}, J4<:(ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#hpf#815"{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#815"), _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 zero 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.open_loop — Method

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

Symbolics.connect — Method

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.AnalysisPoint — Type

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.

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 = ODESystem(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

ModelingToolkit.AnalysisPoint — Method

AnalysisPoint(name::Symbol) -> AnalysisPoint

Create an AnalysisPoint with the given name, with no input or outputs specified.

source