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_sensitivityget the sensitivity function (wiki), $S(s)$, as defined in the field of control theory.get_comp_sensitivityget the complementary sensitivity function $T(s) : S(s)+T(s)=1$.get_looptransferget 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)$.linearizecan be called with two analysis points denoting the input and output of the linearized system.open_loopreturn 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)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
using ModelingToolkit: t_nounits as t
@named P = FirstOrder(k = 1, T = 1) # A first-order system with pole in -1
@named C = Gain(-1) # A P controller
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 = System(eqs, t, systems = [P, C], name = :feedback_system)
matrices_S = get_sensitivity(sys, :plant_input)[1] # Compute the matrices of a state-space representation of the (input)sensitivity function.
matrices_T = get_comp_sensitivity(sys, :plant_input)[1](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 modelWe 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 modelwhich is equivalent to the following with ControlSystems
L = P * (-C) # Add the minus sign to build the negative feedback into the controllerControlSystemsBase.StateSpace{ControlSystemsBase.Continuous, Float64}
A =
-1.0
B =
-1.0
C =
1.0
D =
-0.0
Continuous-time state-space modelTo 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 modelObtaining 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 modelGain 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 = [-0.0;;], pm = [180.0;;])(they are infinite for this system). A Nyquist plot can be produced using
nyquistplot(P)Operating Point of Unconnected Inputs After Loop Openings
When a connection is broken via loop_openings, the downstream input variable(s) that were previously driven by the connection become free. The value assigned to these inputs determines the operating point at which the linearization is computed, which is significant for nonlinear systems where the Jacobian depends on the operating point. Semantics: The variable introduced by a loop opening is treated as a parameter of the system, not as an additional input to the linearization. It does not appear as an extra column in the $B$ or $D$ matrices of the linearized state-space model. No default value is automatically propagated from the output side of the broken connection — the user is expected to provide the value explicitly via the operating point (e.g., op = [u => value]). For initialization, variables that become free due to loop openings are treated as solvable parameters: the initialization system always considers them as unknowns. If the user provides a value in the operating point, that value is used. If no value is provided, the initialization system will attempt to determine a consistent value, but the system may be underdetermined and a warning will be issued. Motivation: This design reflects a series of tradeoffs discovered through iteration:
- Defaulting free inputs to zero (the original behavior) is incorrect for nonlinear systems because it changes the linearization point relative to the equilibrium, potentially yielding meaningless results.
- Automatically propagating the output value of the broken connection preserves the correct operating point in some cases, makes it impossible to disconnect input sources (e.g., a Step signal). It also silently determines the operating point in a way that is difficult for the user to inspect or override.
- Making the variable a parameter (rather than a linearization input) prevents it from appearing as an extra input dimension in the linearized system. The user requested a linearization from specific inputs to specific outputs; the broken connection's input is not one of them.
- Requiring the user to explicitly provide the operating point for broken connections makes the linearization well-defined and inspectable. If the user wants the value that would have been present with the loop closed, they can determine this from a simulation or an initialization solution and pass it explicitly.
Index
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 toModelingToolkit.linearize
See also get_sensitivity, get_looptransfer.
ModelingToolkit.get_comp_sensitivity_function — Method
get_comp_sensitivity_function(
sys::ModelingToolkitBase.AbstractSystem,
aps;
kwargs...
) -> Tuple{ModelingToolkit.LinearizationFunction{Vector{Int64}, Vector{Int64}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}, P, H, _A, J1, J2, J3, J4, IA, @NamedTuple{abstol::Float64, reltol::Float64, nlsolve_alg::Nothing}} where {P<:ODEProblem, H<:Union{Expr, ModelingToolkitBase.GeneratedFunctionWrapper}, _A, J1<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing}}, ModelingToolkit.var"#uff#6"{var"#50#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#6"), _A, _B}), var"#50#fun", _A}}, J2<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, P, F, B, ADTypes.AutoForwardDiff{nothing, Nothing}} where {P, F, B}}, J3<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#pff#7"{var"#51#fun", var"#52#setter"}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Vector{T} where T<:(ForwardDiff.Dual{T, _A, _B} where {_B, _A, T})}} where {T<:(ForwardDiff.Tag{F} where F<:(ModelingToolkit.var"#pff#7"{_A, SymbolicIndexingInterface.OOPSetter{_A1, System, var"#s177"}} where {_A, _A1, T, var"#s177"<:AbstractVector{T}})), _A, _B}), var"#51#fun", var"#52#setter", _A}}, J4<:(ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#hpf#5"{fun, setter}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Vector{T} where T<:(ForwardDiff.Dual{T, _A, _B} where {_B, _A, T<:(ForwardDiff.Tag{F} where F<:(ModelingToolkit.var"#hpf#5"{_A, SymbolicIndexingInterface.OOPSetter{_A1, System, var"#s177"}} where {_A, _A1, T, var"#s177"<:AbstractVector{T}}))})}} where {T<:(ForwardDiff.Tag{F} where F<:(ModelingToolkit.var"#hpf#5"{_A, SymbolicIndexingInterface.OOPSetter{_A1, System, var"#s177"}} where {_A, _A1, T, var"#s177"<:AbstractVector{T}})), _A, _B}), fun, setter, _A}), IA<:Union{SciMLBase.NoInit, SciMLBase.OverrideInit{Nothing, Nothing, Nothing}}}, System}
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 tolinearization_function.
All other keyword arguments are forwarded to linearization_function.
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.
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 toModelingToolkit.linearize
See also get_sensitivity, get_comp_sensitivity, open_loop.
ModelingToolkit.get_looptransfer_function — Method
get_looptransfer_function(
sys::ModelingToolkitBase.AbstractSystem,
aps;
kwargs...
) -> Tuple{ModelingToolkit.LinearizationFunction{Vector{Int64}, Vector{Int64}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}, P, H, _A, J1, J2, J3, J4, IA, @NamedTuple{abstol::Float64, reltol::Float64, nlsolve_alg::Nothing}} where {P<:ODEProblem, H<:Union{Expr, ModelingToolkitBase.GeneratedFunctionWrapper}, _A, J1<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing}}, ModelingToolkit.var"#uff#6"{var"#50#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#6"), _A, _B}), var"#50#fun", _A}}, J2<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, P, F, B, ADTypes.AutoForwardDiff{nothing, Nothing}} where {P, F, B}}, J3<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#pff#7"{var"#51#fun", var"#52#setter"}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Vector{T} where T<:(ForwardDiff.Dual{T, _A, _B} where {_B, _A, T})}} where {T<:(ForwardDiff.Tag{F} where F<:(ModelingToolkit.var"#pff#7"{_A, SymbolicIndexingInterface.OOPSetter{_A1, System, var"#s177"}} where {_A, _A1, T, var"#s177"<:AbstractVector{T}})), _A, _B}), var"#51#fun", var"#52#setter", _A}}, J4<:(ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#hpf#5"{fun, setter}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Vector{T} where T<:(ForwardDiff.Dual{T, _A, _B} where {_B, _A, T<:(ForwardDiff.Tag{F} where F<:(ModelingToolkit.var"#hpf#5"{_A, SymbolicIndexingInterface.OOPSetter{_A1, System, var"#s177"}} where {_A, _A1, T, var"#s177"<:AbstractVector{T}}))})}} where {T<:(ForwardDiff.Tag{F} where F<:(ModelingToolkit.var"#hpf#5"{_A, SymbolicIndexingInterface.OOPSetter{_A1, System, var"#s177"}} where {_A, _A1, T, var"#s177"<:AbstractVector{T}})), _A, _B}), fun, setter, _A}), IA<:Union{SciMLBase.NoInit, SciMLBase.OverrideInit{Nothing, Nothing, Nothing}}}, System}
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 tolinearization_function.
All other keyword arguments are forwarded to linearization_function.
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 toModelingToolkit.linearize
See also get_comp_sensitivity, get_looptransfer.
ModelingToolkit.get_sensitivity_function — Method
get_sensitivity_function(
sys::ModelingToolkitBase.AbstractSystem,
aps;
kwargs...
) -> Tuple{ModelingToolkit.LinearizationFunction{Vector{Int64}, Vector{Int64}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}, P, H, _A, J1, J2, J3, J4, IA, @NamedTuple{abstol::Float64, reltol::Float64, nlsolve_alg::Nothing}} where {P<:ODEProblem, H<:Union{Expr, ModelingToolkitBase.GeneratedFunctionWrapper}, _A, J1<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing}}, ModelingToolkit.var"#uff#6"{var"#50#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#6"), _A, _B}), var"#50#fun", _A}}, J2<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, P, F, B, ADTypes.AutoForwardDiff{nothing, Nothing}} where {P, F, B}}, J3<:Union{Nothing, ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#pff#7"{var"#51#fun", var"#52#setter"}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Vector{T} where T<:(ForwardDiff.Dual{T, _A, _B} where {_B, _A, T})}} where {T<:(ForwardDiff.Tag{F} where F<:(ModelingToolkit.var"#pff#7"{_A, SymbolicIndexingInterface.OOPSetter{_A1, System, var"#s177"}} where {_A, _A1, T, var"#s177"<:AbstractVector{T}})), _A, _B}), var"#51#fun", var"#52#setter", _A}}, J4<:(ModelingToolkit.PreparedJacobian{true, DifferentiationInterfaceForwardDiffExt.ForwardDiffTwoArgJacobianPrep{Nothing, C, Tuple{Nothing, Nothing, Nothing}}, ModelingToolkit.var"#hpf#5"{fun, setter}, _A, ADTypes.AutoForwardDiff{nothing, Nothing}} where {C<:(ForwardDiff.JacobianConfig{T, _A, _B, <:Tuple{Any, Vector{T} where T<:(ForwardDiff.Dual{T, _A, _B} where {_B, _A, T<:(ForwardDiff.Tag{F} where F<:(ModelingToolkit.var"#hpf#5"{_A, SymbolicIndexingInterface.OOPSetter{_A1, System, var"#s177"}} where {_A, _A1, T, var"#s177"<:AbstractVector{T}}))})}} where {T<:(ForwardDiff.Tag{F} where F<:(ModelingToolkit.var"#hpf#5"{_A, SymbolicIndexingInterface.OOPSetter{_A1, System, var"#s177"}} where {_A, _A1, T, var"#s177"<:AbstractVector{T}})), _A, _B}), fun, setter, _A}), IA<:Union{SciMLBase.NoInit, SciMLBase.OverrideInit{Nothing, Nothing, Nothing}}}, System}
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 tolinearization_function.
All other keyword arguments are forwarded to linearization_function.
ModelingToolkit.linearization_ap_transform — Method
sys, input_vars, output_vars =linearization_ap_transform(
sys,
inputs::Union{Vector{Symbol}, Vector{AnalysisPoint}, Symbol, AnalysisPoint},
outputs,
loop_openings
) -> Tuple{Any, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}}
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.