From bc2c1287c24fd80ce12aa5068c6924c3395665c7 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 19 Dec 2024 13:37:32 +0530
Subject: [PATCH 01/23] feat: add analysis points

---
 src/ModelingToolkit.jl         |   4 +-
 src/systems/analysis_points.jl | 649 +++++++++++++++++++++++++++++++++
 src/systems/connectors.jl      |   1 +
 3 files changed, 653 insertions(+), 1 deletion(-)
 create mode 100644 src/systems/analysis_points.jl

diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 10aba4d8a9..486ef63460 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -145,6 +145,7 @@ include("systems/index_cache.jl")
 include("systems/parameter_buffer.jl")
 include("systems/abstractsystem.jl")
 include("systems/model_parsing.jl")
+include("systems/analysis_points.jl")
 include("systems/connectors.jl")
 include("systems/imperative_affect.jl")
 include("systems/callbacks.jl")
@@ -238,7 +239,8 @@ export SteadyStateProblem, SteadyStateProblemExpr
 export JumpProblem
 export NonlinearSystem, OptimizationSystem, ConstraintsSystem
 export alias_elimination, flatten
-export connect, domain_connect, @connector, Connection, Flow, Stream, instream
+export connect, domain_connect, @connector, Connection, AnalysisPoint, Flow, Stream,
+       instream
 export initial_state, transition, activeState, entry, ticksInState, timeInState
 export @component, @mtkmodel, @mtkbuild
 export isinput, isoutput, getbounds, hasbounds, getguess, hasguess, isdisturbance,
diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
new file mode 100644
index 0000000000..605eb32749
--- /dev/null
+++ b/src/systems/analysis_points.jl
@@ -0,0 +1,649 @@
+"""
+    $(TYPEDEF)
+    $(TYPEDSIGNATURES)
+
+Create an AnalysisPoint for linear analysis. Analysis points can be created by calling
+
+```
+connect(in, :ap_name, out...)
+```
+
+Where `in` is the input to the connection, and `out...` are the outputs. In the context of
+ModelingToolkitStandardLibrary.jl, `in` is a `RealOutput` connector and `out...` are all
+`RealInput` connectors. All involved connectors are required to either have an unknown named
+`u` or a single unknown, all of which should have the same size.
+"""
+struct AnalysisPoint
+    input::Any
+    name::Symbol
+    outputs::Union{Nothing, Vector{Any}}
+end
+
+AnalysisPoint() = AnalysisPoint(nothing, Symbol(), nothing)
+AnalysisPoint(name::Symbol) = AnalysisPoint(nothing, name, nothing)
+
+Base.nameof(ap::AnalysisPoint) = ap.name
+
+Base.show(io::IO, ap::AnalysisPoint) = show(io, MIME"text/plain"(), ap)
+function Base.show(io::IO, ::MIME"text/plain", ap::AnalysisPoint)
+    if ap.input === nothing
+        print(io, "0")
+        return
+    end
+    if get(io, :compact, false)
+        print(io,
+            "AnalysisPoint($(ap_var(ap.input)), $(ap_var.(ap.outputs)); name=$(ap.name))")
+    else
+        print(io, "AnalysisPoint(")
+        printstyled(io, ap.name, color = :cyan)
+        if ap.input !== nothing && ap.outputs !== nothing
+            print(io, " from ")
+            printstyled(io, ap_var(ap.input), color = :green)
+            print(io, " to ")
+            if length(ap.outputs) == 1
+                printstyled(io, ap_var(ap.outputs[1]), color = :blue)
+            else
+                printstyled(io, "[", join(ap_var.(ap.outputs), ", "), "]", color = :blue)
+            end
+        end
+        print(io, ")")
+    end
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Convert an `AnalysisPoint` to a standard connection.
+"""
+function to_connection(ap::AnalysisPoint)
+    return connect(ap.input, ap.outputs...)
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Namespace an `AnalysisPoint` by namespacing the involved systems and the name of the point.
+"""
+function renamespace(sys, ap::AnalysisPoint)
+    return AnalysisPoint(
+        ap.input === nothing ? nothing : renamespace(sys, ap.input),
+        renamespace(sys, ap.name),
+        ap.outputs === nothing ? nothing : map(Base.Fix1(renamespace, sys), ap.outputs)
+    )
+end
+
+# create analysis points via `connect`
+function Symbolics.connect(in, ap::AnalysisPoint, outs...)
+    return AnalysisPoint() ~ AnalysisPoint(in, ap.name, collect(outs))
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Create an `AnalysisPoint` connection connecting `in` to `outs...`.
+"""
+function Symbolics.connect(in::AbstractSystem, name::Symbol, out, outs...)
+    return AnalysisPoint() ~ AnalysisPoint(in, name, [out; collect(outs)])
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Return all the namespaces in `name`. Namespaces should be separated by `.` or
+`$NAMESPACE_SEPARATOR`.
+"""
+namespace_hierarchy(name::Symbol) = map(
+    Symbol, split(string(name), ('.', NAMESPACE_SEPARATOR)))
+
+"""
+    $(TYPEDSIGNATURES)
+
+Remove all `AnalysisPoint`s in `sys` and any of its subsystems, replacing them by equivalent connections.
+"""
+function remove_analysis_points(sys::AbstractSystem)
+    eqs = map(get_eqs(sys)) do eq
+        eq.lhs isa AnalysisPoint ? to_connection(eq.rhs) : eq
+    end
+    @set! sys.eqs = eqs
+    @set! sys.systems = map(remove_analysis_points, get_systems(sys))
+
+    return sys
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Given a system involved in an `AnalysisPoint`, get the variable to be used in the
+connection. This is the variable named `u` if present, and otherwise the only
+variable in the system. If the system does not have a variable named `u` and
+contains multiple variables, throw an error.
+"""
+function ap_var(sys)
+    if hasproperty(sys, :u)
+        return sys.u
+    end
+    x = unknowns(sys)
+    length(x) == 1 && return renamespace(sys, x[1])
+    error("Could not determine the analysis-point variable in system $(nameof(sys)). To use an analysis point, apply it to a connection between causal blocks which have a variable named `u` or a single unknown of the same size.")
+end
+
+"""
+    $(TYPEDEF)
+
+The supertype of all transformations that can be applied to an `AnalysisPoint`. All
+concrete subtypes must implement `apply_transformation`.
+"""
+abstract type AnalysisPointTransformation end
+
+"""
+    apply_transformation(tf::AnalysisPointTransformation, sys::AbstractSystem)
+
+Apply the given analysis point transformation `tf` to the system `sys`. Throw an error if
+any analysis points referred to in `tf` are not present in `sys`. Return a tuple
+containing the modified system as the first element, and a tuple of the additional
+variables added by the transformation as the second element.
+"""
+function apply_transformation end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Given a namespaced subsystem `target` of root system `root`, return a modified copy of
+`root` with `target` modified according to `fn` alongside any extra variables added
+by `fn`.
+
+`fn` is a function which takes the instance of `target` present in the hierarchy of
+`root`, and returns a 2-tuple consisting of the modified version of `target` and a tuple
+of the extra variables added.
+"""
+modify_nested_subsystem(fn, root::AbstractSystem, target::AbstractSystem) = modify_nested_subsystem(
+    fn, root, nameof(target))
+"""
+    $(TYPEDSIGNATURES)
+
+Apply the modification to the system containing the namespaced analysis point `target`.
+"""
+modify_nested_subsystem(fn, root::AbstractSystem, target::AnalysisPoint) = modify_nested_subsystem(
+    fn, root, @view namespace_hierarchy(nameof(target))[1:(end - 1)])
+"""
+    $(TYPEDSIGNATURES)
+
+Apply the modification to the nested subsystem of `root` whose namespaced name matches
+the provided name `target`. The namespace separator in `target` should be `.` or
+`$NAMESPACE_SEPARATOR`. The `target` may include `nameof(root)` as the first namespace.
+"""
+modify_nested_subsystem(fn, root::AbstractSystem, target::Symbol) = modify_nested_subsystem(
+    fn, root, namespace_hierarchy(target))
+
+"""
+    $(TYPEDSIGNATURES)
+
+Apply the modification to the nested subsystem of `root` where the name of the subsystem at
+each level in the hierarchy is given by elements of `hierarchy`. For example, if
+`hierarchy = [:a, :b, :c]`, the system being searched for will be `root.a.b.c`. Note that
+the hierarchy may include the name of the root system, in which the first element will be
+ignored. For example, `hierarchy = [:root, :a, :b, :c]` also searches for `root.a.b.c`.
+An empty `hierarchy` will apply the modification to `root`.
+"""
+function modify_nested_subsystem(
+        fn, root::AbstractSystem, hierarchy::AbstractVector{Symbol})
+    # no hierarchy, so just apply to the root
+    if isempty(hierarchy)
+        return fn(root)
+    end
+    # ignore the name of the root
+    if nameof(root) == hierarchy[1]
+        hierarchy = @view hierarchy[2:end]
+    end
+
+    # recursive helper function which does the searching and modification
+    function _helper(sys::AbstractSystem, i::Int)
+        # we reached past the end, so everything matched and
+        # `sys` is the system to modify.
+        if i > length(hierarchy)
+            sys, vars = fn(sys)
+        else
+            cur = hierarchy[i]
+            idx = findfirst(subsys -> nameof(subsys) == cur, get_systems(sys))
+            idx === nothing &&
+                error("System $(join([nameof(root); hierarchy[1:i-1]], '.')) does not have a subsystem named $cur.")
+
+            newsys, vars = _helper(get_systems(sys)[idx], i + 1)
+            @set! sys.systems[idx] = newsys
+        end
+        if i != 1
+            vars = ntuple(Val(length(vars))) do i
+                renamespace(sys, vars[i])
+            end
+        end
+        return sys, vars
+    end
+
+    return _helper(root, 1)
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Given a system `sys` and analysis point `ap`, return the index in `get_eqs(sys)`
+containing an equation which has as it's RHS an analysis point with name `nameof(ap)`.
+"""
+analysis_point_index(sys::AbstractSystem, ap::AnalysisPoint) = analysis_point_index(
+    sys, nameof(ap))
+"""
+    $(TYPEDSIGNATURES)
+
+Search for the analysis point with the given `name` in `get_eqs(sys)`.
+"""
+function analysis_point_index(sys::AbstractSystem, name::Symbol)
+    name = namespace_hierarchy(name)[end]
+    findfirst(get_eqs(sys)) do eq
+        eq.lhs isa AnalysisPoint && nameof(eq.rhs) == name
+    end
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Create a new variable of the same `symtype` and size as `var`, using `name` as the base
+name for the new variable. `iv` denotes the independent variable of the system. Prefix
+`d_` to the name of the new variable if `perturb == true`. Return the new symbolic
+variable and the appropriate zero value for it.
+"""
+function get_analysis_variable(var, name, iv; perturb = true)
+    var = unwrap(var)
+    if perturb
+        name = Symbol(:d_, name)
+    end
+    if Symbolics.isarraysymbolic(var)
+        T = Array{eltype(symtype(var)), ndims(var)}
+        pvar = unwrap(only(@variables $name(iv)::T))
+        pvar = setmetadata(pvar, Symbolics.ArrayShapeCtx, Symbolics.shape(var))
+        default = zeros(symtype(var), size(var))
+    else
+        T = symtype(var)
+        pvar = unwrap(only(@variables $name(iv)::T))
+        default = zero(T)
+    end
+    return pvar, default
+end
+
+#### PRIMITIVE TRANSFORMATIONS
+
+"""
+    $(TYPEDEF)
+
+A transformation which breaks the connection referred to by `ap`. If `add_input == true`,
+it will add a new input variable which connects to the outputs of the analysis point.
+`apply_transformation` returns the new input variable (if added) as the auxiliary
+information. The new input variable will have the name `Symbol(:d_, nameof(ap))`.
+
+Note that this transformation will remove `ap`, causing any subsequent transformations
+referring to it to fail.
+
+The added variable, if present, will have a default of zero, of the appropriate type and
+size.
+
+## Fields
+
+$(TYPEDFIELDS)
+"""
+struct Break <: AnalysisPointTransformation
+    """
+    The analysis point to break.
+    """
+    ap::AnalysisPoint
+    """
+    Whether to add a new input variable connected to all the outputs of `ap`.
+    """
+    add_input::Bool
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+`Break` the given analysis point `ap` without adding an input.
+"""
+Break(ap::AnalysisPoint) = Break(ap, false)
+
+function apply_transformation(tf::Break, sys::AbstractSystem)
+    modify_nested_subsystem(sys, tf.ap) do breaksys
+        # get analysis point
+        ap_idx = analysis_point_index(breaksys, tf.ap)
+        ap_idx === nothing &&
+            error("Analysis point $(nameof(tf.ap)) not found in system $(nameof(sys)).")
+        breaksys_eqs = copy(get_eqs(breaksys))
+        @set! breaksys.eqs = breaksys_eqs
+
+        ap = breaksys_eqs[ap_idx].rhs
+        deleteat!(breaksys_eqs, ap_idx)
+
+        tf.add_input || return sys, ()
+
+        ap_ivar = ap_var(ap.input)
+        new_var, new_def = get_analysis_variable(ap_ivar, nameof(ap), get_iv(sys))
+        for outsys in ap.outputs
+            push!(breaksys_eqs, ap_var(outsys) ~ new_var)
+        end
+        defs = copy(get_defaults(breaksys))
+        defs[new_var] = new_def
+        @set! breaksys.defaults = defs
+        unks = copy(get_unknowns(breaksys))
+        push!(unks, new_var)
+        @set! breaksys.unknowns = unks
+
+        return breaksys, (new_var,)
+    end
+end
+
+"""
+    $(TYPEDEF)
+
+A transformation which returns the variable corresponding to the input of the analysis
+point.
+
+`apply_transformation` returns the variable as auxiliary information.
+
+## Fields
+
+$(TYPEDFIELDS)
+"""
+struct GetInput <: AnalysisPointTransformation
+    """
+    The analysis point to add the output to.
+    """
+    ap::AnalysisPoint
+end
+
+function apply_transformation(tf::GetInput, sys::AbstractSystem)
+    modify_nested_subsystem(sys, tf.ap) do ap_sys
+        # get the analysis point
+        ap_idx = analysis_point_index(ap_sys, tf.ap)
+        ap_idx === nothing &&
+            error("Analysis point $(nameof(tf.ap)) not found in system $(nameof(sys)).")
+        # get the anlysis point
+        ap_sys_eqs = copy(get_eqs(ap_sys))
+        ap = ap_sys_eqs[ap_idx].rhs
+
+        # input variable
+        ap_ivar = ap_var(ap.input)
+        return ap_sys, (ap_ivar,)
+    end
+end
+
+"""
+    $(TYPEDEF)
+
+A transformation that creates a new input variable which is added to the input of
+the analysis point before connecting to the outputs. The new variable will have the name
+`Symbol(:d_, nameof(ap))`.
+
+If `with_output == true`, also creates an additional new variable which has the value
+provided to the outputs after the above modification. This new variable has the same name
+as the analysis point.
+
+`apply_transformation` returns a 1-tuple of the perturbation variable if
+`with_output == false` and a 2-tuple of the perturbation variable and output variable if
+`with_output == true`.
+
+Removes the analysis point `ap`, so any subsequent transformations requiring it will fail.
+
+The added variable(s) will have a default of zero, of the appropriate type and size.
+
+## Fields
+
+$(TYPEDFIELDS)
+"""
+struct PerturbOutput <: AnalysisPointTransformation
+    """
+    The analysis point to modify
+    """
+    ap::AnalysisPoint
+    """
+    Whether to add an additional output variable.
+    """
+    with_output::Bool
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Add an input without an additional output variable.
+"""
+PerturbOutput(ap::AnalysisPoint) = PerturbOutput(ap, false)
+
+function apply_transformation(tf::PerturbOutput, sys::AbstractSystem)
+    modify_nested_subsystem(sys, tf.ap) do ap_sys
+        # get analysis point
+        ap_idx = analysis_point_index(ap_sys, tf.ap)
+        ap_idx === nothing &&
+            error("Analysis point $(nameof(tf.ap)) not found in system $(nameof(sys)).")
+        # modified quations
+        ap_sys_eqs = copy(get_eqs(ap_sys))
+        @set! ap_sys.eqs = ap_sys_eqs
+        ap = ap_sys_eqs[ap_idx].rhs
+        # remove analysis point
+        deleteat!(ap_sys_eqs, ap_idx)
+
+        # add equations involving new variable
+        ap_ivar = ap_var(ap.input)
+        new_var, new_def = get_analysis_variable(ap_ivar, nameof(ap), get_iv(sys))
+        for outsys in ap.outputs
+            push!(ap_sys_eqs, ap_var(outsys) ~ ap_ivar + new_var)
+        end
+        # add variable
+        unks = copy(get_unknowns(ap_sys))
+        push!(unks, new_var)
+        @set! ap_sys.unknowns = unks
+        # add default
+        defs = copy(get_defaults(ap_sys))
+        defs[new_var] = new_def
+        @set! ap_sys.defaults = defs
+
+        tf.with_output || return ap_sys, (new_var,)
+
+        # add output variable, equation, default
+        out_var, out_def = get_analysis_variable(
+            ap_ivar, nameof(ap), get_iv(sys); perturb = false)
+        defs[out_var] = out_def
+        push!(ap_sys_eqs, out_var ~ ap_ivar + new_var)
+        push!(unks, out_var)
+
+        return ap_sys, (new_var, out_var)
+    end
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+A transformation which adds a variable named `name` to the system containing the analysis
+point `ap`. The added variable has the same type and size as the input of the analysis
+point.
+"""
+struct AddVariable <: AnalysisPointTransformation
+    ap::AnalysisPoint
+    name::Symbol
+end
+
+AddVariable(ap::AnalysisPoint) = AddVariable(ap, nameof(ap))
+
+function apply_transformation(tf::AddVariable, sys::AbstractSystem)
+    modify_nested_subsystem(sys, tf.ap) do ap_sys
+        # get analysis point
+        ap_idx = analysis_point_index(ap_sys, tf.ap)
+        ap_idx === nothing &&
+            error("Analysis point $(nameof(tf.ap)) not found in system $(nameof(sys)).")
+        ap_sys_eqs = copy(get_eqs(ap_sys))
+        ap = ap_sys_eqs[ap_idx].rhs
+
+        # add equations involving new variable
+        ap_ivar = ap_var(ap.input)
+        new_var, new_def = get_analysis_variable(
+            ap_ivar, tf.name, get_iv(sys); perturb = false)
+        # add variable
+        unks = copy(get_unknowns(ap_sys))
+        push!(unks, new_var)
+        @set! ap_sys.unknowns = unks
+        return ap_sys, (new_var,)
+    end
+end
+
+#### DERIVED TRANSFORMATIONS
+
+"""
+    $(TYPEDSIGNATURES)
+
+A transformation enable calculating the sensitivity function about the analysis point `ap`.
+`apply_transformation` returns a 2-tuple `du, u` as auxiliary information.
+
+Removes the analysis point `ap`, so any subsequent transformations requiring it will fail.
+
+The added variables will have a default of zero, of the appropriate type and size.
+"""
+SensitivityTransform(ap::AnalysisPoint) = PerturbOutput(ap, true)
+
+"""
+    $(TYPEDEF)
+
+A transformation to enable calculating the complementary sensitivity function about the
+analysis point `ap`. `apply_transformation` returns a 2-tuple `du, u` as auxiliary
+information.
+
+Removes the analysis point `ap`, so any subsequent transformations requiring it will fail.
+
+The added variables will have a default of zero, of the appropriate type and size.
+"""
+struct ComplementarySensitivityTransform <: AnalysisPointTransformation
+    ap::AnalysisPoint
+end
+
+function apply_transformation(cst::ComplementarySensitivityTransform, sys::AbstractSystem)
+    sys, (u,) = apply_transformation(GetInput(cst.ap), sys)
+    sys, (du,) = apply_transformation(AddVariable(cst.ap, Symbol(:comp_sens_du)), sys)
+    sys, (_du,) = apply_transformation(PerturbOutput(cst.ap), sys)
+
+    # `PerturbOutput` adds the equation `input + _du ~ output`
+    # but comp sensitivity wants `output + du ~ input`. Thus, `du ~ -_du`.
+    eqs = copy(get_eqs(sys))
+    @set! sys.eqs = eqs
+    push!(eqs, du ~ -_du)
+
+    defs = copy(get_defaults(sys))
+    @set! sys.defaults = defs
+    defs[du] = -_du
+    return sys, (du, u)
+end
+
+struct LoopTransferTransform <: AnalysisPointTransformation
+    ap::AnalysisPoint
+end
+
+function apply_transformation(tf::LoopTransferTransform, sys::AbstractSystem)
+    sys, (u,) = apply_transformation(GetInput(tf.ap), sys)
+    sys, (du,) = apply_transformation(Break(tf.ap, true), sys)
+    return sys, (du, u)
+end
+
+### TODO: Move these
+
+canonicalize_ap(ap::Symbol) = [AnalysisPoint(ap)]
+canonicalize_ap(ap::AnalysisPoint) = [ap]
+canonicalize_ap(ap) = [ap]
+function canonicalize_ap(aps::Vector)
+    mapreduce(canonicalize_ap, vcat, aps; init = [])
+end
+
+function handle_loop_openings(sys::AbstractSystem, aps)
+    for ap in canonicalize_ap(aps)
+        sys, (outvar,) = apply_transformation(Break(ap, true), sys)
+        if Symbolics.isarraysymbolic(outvar)
+            push!(get_eqs(sys), outvar ~ zeros(size(outvar)))
+        else
+            push!(get_eqs(sys), outvar ~ 0)
+        end
+    end
+    return sys
+end
+
+function get_sensitivity_function(
+        sys::AbstractSystem, aps; system_modifier = identity, loop_openings = [], kwargs...)
+    sys = handle_loop_openings(sys, loop_openings)
+    aps = canonicalize_ap(aps)
+    dus = []
+    us = []
+    for ap in aps
+        sys, (du, u) = apply_transformation(SensitivityTransform(ap), sys)
+        push!(dus, du)
+        push!(us, u)
+    end
+    linearization_function(system_modifier(sys), dus, us; kwargs...)
+end
+
+function get_comp_sensitivity_function(
+        sys::AbstractSystem, aps; system_modifier = identity, loop_openings = [], kwargs...)
+    sys = handle_loop_openings(sys, loop_openings)
+    aps = canonicalize_ap(aps)
+    dus = []
+    us = []
+    for ap in aps
+        sys, (du, u) = apply_transformation(ComplementarySensitivityTransform(ap), sys)
+        push!(dus, du)
+        push!(us, u)
+    end
+    linearization_function(system_modifier(sys), dus, us; kwargs...)
+end
+
+function get_looptransfer_function(
+        sys, aps; system_modifier = identity, loop_openings = [], kwargs...)
+    sys = handle_loop_openings(sys, loop_openings)
+    aps = canonicalize_ap(aps)
+    dus = []
+    us = []
+    for ap in aps
+        sys, (du, u) = apply_transformation(LoopTransferTransform(ap), sys)
+        push!(dus, du)
+        push!(us, u)
+    end
+    linearization_function(system_modifier(sys), dus, us; kwargs...)
+end
+
+for f in [:get_sensitivity, :get_comp_sensitivity, :get_looptransfer]
+    @eval function $f(sys, ap, args...; kwargs...)
+        lin_fun, ssys = $(Symbol(f, :_function))(sys, ap, args...; kwargs...)
+        ModelingToolkit.linearize(ssys, lin_fun; kwargs...), ssys
+    end
+end
+
+function open_loop(sys, ap::Union{Symbol, AnalysisPoint}; kwargs...)
+    if ap isa Symbol
+        ap = AnalysisPoint(ap)
+    end
+    tf = LoopTransferTransform(ap)
+    return apply_transformation(tf, sys)
+end
+
+function linearization_function(sys::AbstractSystem,
+        inputs::Union{Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}},
+        outputs; loop_openings = [], system_modifier = identity, kwargs...)
+    sys = handle_loop_openings(sys, loop_openings)
+
+    inputs = canonicalize_ap(inputs)
+    outputs = canonicalize_ap(outputs)
+
+    input_vars = []
+    for input in inputs
+        sys, (input_var,) = apply_transformation(PerturbOutput(input), sys)
+        push!(input_vars, input_var)
+    end
+    output_vars = []
+    for output in outputs
+        if output isa AnalysisPoint
+            sys, (output_var,) = apply_transformation(GetInput(output), sys)
+            push!(output_vars, output_var)
+        else
+            push!(output_vars, output)
+        end
+    end
+
+    return linearization_function(system_modifier(sys), input_vars, output_vars; kwargs...)
+end
diff --git a/src/systems/connectors.jl b/src/systems/connectors.jl
index 227b4624bf..f0f49acb85 100644
--- a/src/systems/connectors.jl
+++ b/src/systems/connectors.jl
@@ -481,6 +481,7 @@ end
 
 function expand_connections(sys::AbstractSystem, find = nothing, replace = nothing;
         debug = false, tol = 1e-10, scalarize = true)
+    sys = remove_analysis_points(sys)
     sys, (csets, domain_csets) = generate_connection_set(sys, find, replace; scalarize)
     ceqs, instream_csets = generate_connection_equations_and_stream_connections(csets)
     _sys = expand_instream(instream_csets, sys; debug = debug, tol = tol)

From c43b3e9cd0f94658bf29f8cd01490a8e40082e6b Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 19 Dec 2024 13:37:40 +0530
Subject: [PATCH 02/23] test: test analysis points

---
 test/analysis_points.jl | 23 +++++++++++++++++++++++
 test/runtests.jl        |  1 +
 2 files changed, 24 insertions(+)
 create mode 100644 test/analysis_points.jl

diff --git a/test/analysis_points.jl b/test/analysis_points.jl
new file mode 100644
index 0000000000..b6a9c8c21f
--- /dev/null
+++ b/test/analysis_points.jl
@@ -0,0 +1,23 @@
+using ModelingToolkit, ModelingToolkitStandardLibrary
+using Test
+using ModelingToolkit: t_nounits as t, D_nounits as D
+
+@testset "AnalysisPoint is lowered to `connect`" begin
+    @named P = FirstOrder(k = 1, T = 1)
+    @named C = Gain(; k = -1)
+    
+    ap = AnalysisPoint(:plant_input)
+    eqs = [connect(P.output, C.input)
+           connect(C.output, ap, P.input)]
+    sys_ap = ODESystem(eqs, t, systems = [P, C], name = :hej)
+    sys_ap2 = @test_nowarn expand_connections(sys)
+
+    @test all(eq -> !(eq.lhs isa AnalysisPoint), equations(sys_ap2))
+
+    eqs = [connect(P.output, C.input)
+           connect(C.output, P.input)]
+    sys_normal = ODESystem(eqs, t, systems = [P, C], name = :hej)
+    sys_normal2 = @test_nowarn expand_connections(sys)
+
+    @test isequal(sys_ap2, sys_normal2)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 4a2ec80e6a..af83ffccd8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -84,6 +84,7 @@ end
             @safetestset "print_tree" include("print_tree.jl")
             @safetestset "Constraints Test" include("constraints.jl")
             @safetestset "IfLifting Test" include("if_lifting.jl")
+            @safetestset "Analysis Points Test" include("analysis_points.jl")
         end
     end
 

From eba5dc77d855bca64e4c94e0280b1535009b1f64 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 24 Dec 2024 11:24:54 +0530
Subject: [PATCH 03/23] refactor: format Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 7b0997823c..e4eba3cfad 100644
--- a/Project.toml
+++ b/Project.toml
@@ -73,8 +73,8 @@ MTKBifurcationKitExt = "BifurcationKit"
 MTKChainRulesCoreExt = "ChainRulesCore"
 MTKDeepDiffsExt = "DeepDiffs"
 MTKHomotopyContinuationExt = "HomotopyContinuation"
-MTKLabelledArraysExt = "LabelledArrays"
 MTKInfiniteOptExt = "InfiniteOpt"
+MTKLabelledArraysExt = "LabelledArrays"
 
 [compat]
 AbstractTrees = "0.3, 0.4"

From f0bab4a58daff1d4b6e984173631a1ba75c72fa3 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 30 Dec 2024 00:41:06 +0530
Subject: [PATCH 04/23] feat: allow accessing `AnalysisPoint` via `getproperty`
 syntax

---
 src/systems/abstractsystem.jl | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 168260ae69..9b5ff7183d 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -1162,6 +1162,14 @@ function getvar(sys::AbstractSystem, name::Symbol; namespace = !iscomplete(sys))
         end
     end
 
+    if has_eqs(sys)
+        for eq in get_eqs(sys)
+            if eq.lhs isa AnalysisPoint && nameof(eq.rhs) == name
+                return namespace ? renamespace(sys, eq.rhs) : eq.rhs
+            end
+        end
+    end
+
     throw(ArgumentError("System $(nameof(sys)): variable $name does not exist"))
 end
 

From 04ca326ca985bbc28a66b02f0f1ed89143919a77 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 30 Dec 2024 00:43:03 +0530
Subject: [PATCH 05/23] test: add tests from MTKStdlib

---
 test/analysis_points.jl | 242 +++++++++++++++++++++++++++++++++++++++-
 1 file changed, 239 insertions(+), 3 deletions(-)

diff --git a/test/analysis_points.jl b/test/analysis_points.jl
index b6a9c8c21f..6eae30314e 100644
--- a/test/analysis_points.jl
+++ b/test/analysis_points.jl
@@ -1,11 +1,14 @@
-using ModelingToolkit, ModelingToolkitStandardLibrary
+using ModelingToolkit, ModelingToolkitStandardLibrary.Blocks
+using OrdinaryDiffEq, LinearAlgebra
 using Test
-using ModelingToolkit: t_nounits as t, D_nounits as D
+using ModelingToolkit: t_nounits as t, D_nounits as D, AnalysisPoint, get_sensitivity,
+                       get_comp_sensitivity, get_looptransfer
+using Symbolics: NAMESPACE_SEPARATOR
 
 @testset "AnalysisPoint is lowered to `connect`" begin
     @named P = FirstOrder(k = 1, T = 1)
     @named C = Gain(; k = -1)
-    
+
     ap = AnalysisPoint(:plant_input)
     eqs = [connect(P.output, C.input)
            connect(C.output, ap, P.input)]
@@ -21,3 +24,236 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
 
     @test isequal(sys_ap2, sys_normal2)
 end
+
+# also tests `connect(input, name::Symbol, outputs...)` syntax
+@testset "AnalysisPoint is accessible via `getproperty`" begin
+    @named P = FirstOrder(k = 1, T = 1)
+    @named C = Gain(; k = -1)
+
+    eqs = [connect(P.output, C.input), connect(C.output, :plant_input, P.input)]
+    sys_ap = ODESystem(eqs, t, systems = [P, C], name = :hej)
+    ap2 = @test_nowarn sys_ap.plant_input
+    @test nameof(ap2) == Symbol(join(["hej", "plant_input"], NAMESPACE_SEPARATOR))
+    @named sys = ODESystem(Equation[], t; systems = [sys_ap])
+    ap3 = @test_nowarn sys.hej.plant_input
+    @test nameof(ap3) == Symbol(join(["sys", "hej", "plant_input"], NAMESPACE_SEPARATOR))
+    sys = complete(sys)
+    ap4 = sys.hej.plant_input
+    @test nameof(ap4) == Symbol(join(["hej", "plant_input"], NAMESPACE_SEPARATOR))
+end
+
+### Ported from MTKStdlib
+
+@named P = FirstOrder(k = 1, T = 1)
+@named C = Gain(; k = -1)
+
+ap = AnalysisPoint(:plant_input)
+eqs = [connect(P.output, C.input), connect(C.output, ap, P.input)]
+sys = ODESystem(eqs, t, systems = [P, C], name = :hej)
+@named nested_sys = ODESystem(Equation[], t; systems = [sys])
+
+@testset "simplifies and solves" begin
+    ssys = structural_simplify(sys)
+    prob = ODEProblem(ssys, [P.x => 1], (0, 10))
+    sol = solve(prob, Rodas5())
+    @test norm(sol.u[1]) >= 1
+    @test norm(sol.u[end]) < 1e-6 # This fails without the feedback through C
+end
+
+@testset "get_sensitivity - $name" for (name, sys, ap) in [
+    ("inner", sys, sys.plant_input),
+    ("nested", nested_sys, nested_sys.hej.plant_input),
+    ("inner - nonamespace", sys, :plant_input),
+    ("inner - Symbol", sys, nameof(sys.plant_input)),
+    ("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
+]
+    matrices, _ = get_sensitivity(sys, ap)
+    @test matrices.A[] == -2
+    @test matrices.B[] * matrices.C[] == -1 # either one negative
+    @test matrices.D[] == 1
+end
+
+@testset "get_comp_sensitivity - $name" for (name, sys, ap) in [
+    ("inner", sys, sys.plant_input),
+    ("nested", nested_sys, nested_sys.hej.plant_input),
+    ("inner - nonamespace", sys, :plant_input),
+    ("inner - Symbol", sys, nameof(sys.plant_input)),
+    ("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
+]
+    matrices, _ = get_comp_sensitivity(sys, ap)
+    @test matrices.A[] == -2
+    @test matrices.B[] * matrices.C[] == 1 # both positive or negative
+    @test matrices.D[] == 0
+end
+
+#=
+# Equivalent code using ControlSystems. This can be used to verify the expected results tested for above.
+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)
+=#
+
+@testset "get_looptransfer - $name" for (name, sys, ap) in [
+    ("inner", sys, sys.plant_input),
+    ("nested", nested_sys, nested_sys.hej.plant_input),
+    ("inner - nonamespace", sys, :plant_input),
+    ("inner - Symbol", sys, nameof(sys.plant_input)),
+    ("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
+]
+    matrices, _ = get_looptransfer(sys, ap)
+    @test matrices.A[] == -1
+    @test matrices.B[] * matrices.C[] == -1 # either one negative
+    @test matrices.D[] == 0
+end
+
+#=
+# Equivalent code using ControlSystems. This can be used to verify the expected results tested for above.
+using ControlSystemsBase
+P = tf(1.0, [1, 1])
+C = -1
+L = P*C
+=#
+
+@testset "open_loop - $name" for (name, sys, ap) in [
+    ("inner", sys, sys.plant_input),
+    ("nested", nested_sys, nested_sys.hej.plant_input),
+    ("inner - nonamespace", sys, :plant_input),
+    ("inner - Symbol", sys, nameof(sys.plant_input)),
+    ("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
+]
+    open_sys, (du, u) = open_loop(sys, ap)
+    matrices, _ = linearize(open_sys, [du], [u])
+    @test matrices.A[] == -1
+    @test matrices.B[] * matrices.C[] == -1 # either one negative
+    @test matrices.D[] == 0
+end
+
+# Multiple analysis points
+
+eqs = [connect(P.output, :plant_output, C.input)
+       connect(C.output, :plant_input, P.input)]
+sys = ODESystem(eqs, t, systems = [P, C], name = :hej)
+@named nested_sys = ODESystem(Equation[], t; systems = [sys])
+
+@testset "get_sensitivity - $name" for (name, sys, ap) in [
+    ("inner", sys, sys.plant_input),
+    ("nested", nested_sys, nested_sys.hej.plant_input),
+    ("inner - nonamespace", sys, :plant_input),
+    ("inner - Symbol", sys, nameof(sys.plant_input)),
+    ("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
+]
+    matrices, _ = get_sensitivity(sys, ap)
+    @test matrices.A[] == -2
+    @test matrices.B[] * matrices.C[] == -1 # either one negative
+    @test matrices.D[] == 1
+end
+
+@testset "linearize - $name" for (name, sys, inputap, outputap) in [
+    ("inner", sys, sys.plant_input, sys.plant_output),
+    ("nested", nested_sys, nested_sys.hej.plant_input, nested_sys.hej.plant_output)
+]
+    @testset "input - $(typeof(input)), output - $(typeof(output))" for (input, output) in [
+        (inputap, outputap),
+        (nameof(inputap), outputap),
+        (inputap, nameof(outputap)),
+        (nameof(inputap), nameof(outputap)),
+        (inputap, [outputap]),
+        (nameof(inputap), [outputap]),
+        (inputap, [nameof(outputap)]),
+        (nameof(inputap), [nameof(outputap)])
+    ]
+        if input isa Symbol
+            # broken because MTKStdlib defines the method for
+            # `input::Union{Symbol, Vector{Symbol}}` which we can't directly call
+            @test_broken linearize(sys, input, output)
+            linfun, ssys = @invoke linearization_function(sys::AbstractSystem,
+                input::Union{Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}},
+                output::Any)
+            matrices = linearize(ssys, linfun)
+        else
+            matrices, _ = linearize(sys, input, output)
+        end
+        # Result should be the same as feedpack(P, 1), i.e., the closed-loop transfer function from plant input to plant output
+        @test matrices.A[] == -2
+        @test matrices.B[] * matrices.C[] == 1 # both positive
+        @test matrices.D[] == 0
+    end
+end
+
+@testset "linearize - variable output - $name" for (name, sys, input, output) in [
+    ("inner", sys, sys.plant_input, P.output.u),
+    ("nested", nested_sys, nested_sys.hej.plant_input, sys.P.output.u)
+]
+    matrices, _ = linearize(sys, input, [output])
+    @test matrices.A[] == -2
+    @test matrices.B[] * matrices.C[] == 1 # both positive
+    @test matrices.D[] == 0
+end
+
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
+using ModelingToolkitStandardLibrary.Mechanical.Rotational
+using ModelingToolkitStandardLibrary.Blocks: Sine, PID, SecondOrder, Step, RealOutput
+using ModelingToolkit: connect
+
+@testset "Complicated model" begin
+    # Parameters
+    m1 = 1
+    m2 = 1
+    k = 1000 # Spring stiffness
+    c = 10   # Damping coefficient
+    @named inertia1 = Inertia(; J = m1)
+    @named inertia2 = Inertia(; J = m2)
+    @named spring = Spring(; c = k)
+    @named damper = Damper(; d = c)
+    @named torque = Torque()
+
+    function SystemModel(u = nothing; name = :model)
+        eqs = [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)]
+        if u !== nothing
+            push!(eqs, connect(torque.tau, u.output))
+            return ODESystem(eqs, t;
+                systems = [
+                    torque,
+                    inertia1,
+                    inertia2,
+                    spring,
+                    damper,
+                    u
+                ],
+                name)
+        end
+        ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
+    end
+
+    @named r = Step(start_time = 0)
+    model = SystemModel()
+    @named pid = PID(k = 100, Ti = 0.5, Td = 1)
+    @named filt = SecondOrder(d = 0.9, w = 10)
+    @named sensor = AngleSensor()
+    @named er = Add(k2 = -1)
+
+    connections = [connect(r.output, :r, filt.input)
+                   connect(filt.output, er.input1)
+                   connect(pid.ctr_output, :u, model.torque.tau)
+                   connect(model.inertia2.flange_b, sensor.flange)
+                   connect(sensor.phi, :y, er.input2)
+                   connect(er.output, :e, pid.err_input)]
+
+    closed_loop = ODESystem(connections, t, systems = [model, pid, filt, sensor, r, er],
+        name = :closed_loop, defaults = [
+            model.inertia1.phi => 0.0,
+            model.inertia2.phi => 0.0,
+            model.inertia1.w => 0.0,
+            model.inertia2.w => 0.0,
+            filt.x => 0.0,
+            filt.xd => 0.0
+        ])
+
+    sys = structural_simplify(closed_loop)
+    prob = ODEProblem(sys, unknowns(sys) .=> 0.0, (0.0, 4.0))
+    sol = solve(prob, Rodas5P(), reltol = 1e-6, abstol = 1e-9)
+end

From 8a4ca21724b6631d271450b458d4a041db2ab5de Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 30 Dec 2024 00:43:23 +0530
Subject: [PATCH 06/23] feat: export `AnalysisPoint`-related tools

---
 src/ModelingToolkit.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 486ef63460..29288a5742 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -293,4 +293,8 @@ export MTKParameters, reorder_dimension_by_tunables!, reorder_dimension_by_tunab
 
 export HomotopyContinuationProblem
 
+export AnalysisPoint, Break, PerturbOutput, SampleInput, SensitivityTransform,
+       ComplementarySensitivityTransform, LoopTransferTransform, apply_transformation,
+       get_sensitivity_function, get_comp_sensitivity_function, get_looptransfer_function,
+       get_sensitivity, get_comp_sensitivity, get_looptransfer, open_loop
 end # module

From 0750e62f415e7e66cc4abc1fd9b85a3a2b4a764f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 30 Dec 2024 00:43:47 +0530
Subject: [PATCH 07/23] test: add downstream analysis points tests

---
 test/downstream/analysis_points.jl | 233 +++++++++++++++++++++++++++++
 test/runtests.jl                   |   1 +
 2 files changed, 234 insertions(+)
 create mode 100644 test/downstream/analysis_points.jl

diff --git a/test/downstream/analysis_points.jl b/test/downstream/analysis_points.jl
new file mode 100644
index 0000000000..5dc2b8957c
--- /dev/null
+++ b/test/downstream/analysis_points.jl
@@ -0,0 +1,233 @@
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, ControlSystemsBase
+using ModelingToolkitStandardLibrary.Mechanical.Rotational
+using ModelingToolkitStandardLibrary.Blocks
+using ModelingToolkit: connect, AnalysisPoint, t_nounits as t, D_nounits as D
+import ControlSystemsBase as CS
+
+@testset "Complicated model" begin
+    # Parameters
+    m1 = 1
+    m2 = 1
+    k = 1000 # Spring stiffness
+    c = 10   # Damping coefficient
+    @named inertia1 = Inertia(; J = m1)
+    @named inertia2 = Inertia(; J = m2)
+    @named spring = Spring(; c = k)
+    @named damper = Damper(; d = c)
+    @named torque = Torque()
+
+    function SystemModel(u = nothing; name = :model)
+        eqs = [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)]
+        if u !== nothing
+            push!(eqs, connect(torque.tau, u.output))
+            return ODESystem(eqs, t;
+                systems = [
+                    torque,
+                    inertia1,
+                    inertia2,
+                    spring,
+                    damper,
+                    u
+                ],
+                name)
+        end
+        ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
+    end
+
+    @named r = Step(start_time = 0)
+    model = SystemModel()
+    @named pid = PID(k = 100, Ti = 0.5, Td = 1)
+    @named filt = SecondOrder(d = 0.9, w = 10)
+    @named sensor = AngleSensor()
+    @named er = Add(k2 = -1)
+
+    connections = [connect(r.output, :r, filt.input)
+                   connect(filt.output, er.input1)
+                   connect(pid.ctr_output, :u, model.torque.tau)
+                   connect(model.inertia2.flange_b, sensor.flange)
+                   connect(sensor.phi, :y, er.input2)
+                   connect(er.output, :e, pid.err_input)]
+
+    closed_loop = ODESystem(connections, t, systems = [model, pid, filt, sensor, r, er],
+        name = :closed_loop, defaults = [
+            model.inertia1.phi => 0.0,
+            model.inertia2.phi => 0.0,
+            model.inertia1.w => 0.0,
+            model.inertia2.w => 0.0,
+            filt.x => 0.0,
+            filt.xd => 0.0
+        ])
+
+    sys = structural_simplify(closed_loop)
+    prob = ODEProblem(sys, unknowns(sys) .=> 0.0, (0.0, 4.0))
+    sol = solve(prob, Rodas5P(), reltol = 1e-6, abstol = 1e-9)
+
+    matrices, ssys = linearize(closed_loop, AnalysisPoint(:r), AnalysisPoint(:y))
+    lsys = ss(matrices...) |> sminreal
+    @test lsys.nx == 8
+
+    stepres = ControlSystemsBase.step(c2d(lsys, 0.001), 4)
+    @test Array(stepres.y[:])≈Array(sol(0:0.001:4, idxs = model.inertia2.phi)) rtol=1e-4
+
+    matrices, ssys = get_sensitivity(closed_loop, :y)
+    So = ss(matrices...)
+
+    matrices, ssys = get_sensitivity(closed_loop, :u)
+    Si = ss(matrices...)
+
+    @test tf(So) ≈ tf(Si)
+end
+
+@testset "Analysis points with subsystems" begin
+    @named P = FirstOrder(k = 1, T = 1)
+    @named C = Gain(; k = 1)
+    @named add = Blocks.Add(k2 = -1)
+
+    eqs = [connect(P.output, :plant_output, add.input2)
+           connect(add.output, C.input)
+           connect(C.output, :plant_input, P.input)]
+
+    # eqs = [connect(P.output, add.input2)
+    #        connect(add.output, C.input)
+    #        connect(C.output, P.input)]
+
+    sys_inner = ODESystem(eqs, t, systems = [P, C, add], name = :inner)
+
+    @named r = Constant(k = 1)
+    @named F = FirstOrder(k = 1, T = 3)
+
+    eqs = [connect(r.output, F.input)
+           connect(F.output, sys_inner.add.input1)]
+    sys_outer = ODESystem(eqs, t, systems = [F, sys_inner, r], name = :outer)
+
+    # test first that the structural_simplify works correctly
+    ssys = structural_simplify(sys_outer)
+    prob = ODEProblem(ssys, Pair[], (0, 10))
+    @test_nowarn solve(prob, Rodas5())
+
+    matrices, _ = get_sensitivity(sys_outer, sys_outer.inner.plant_input)
+    lsys = sminreal(ss(matrices...))
+    @test lsys.A[] == -2
+    @test lsys.B[] * lsys.C[] == -1 # either one negative
+    @test lsys.D[] == 1
+
+    matrices_So, _ = get_sensitivity(sys_outer, sys_outer.inner.plant_output)
+    lsyso = sminreal(ss(matrices_So...))
+    @test lsys == lsyso || lsys == -1 * lsyso * (-1) # Output and input sensitivities are equal for SISO systems
+end
+
+@testset "multilevel system with loop openings" begin
+    @named P_inner = FirstOrder(k = 1, T = 1)
+    @named feedback = Feedback()
+    @named ref = Step()
+    @named sys_inner = ODESystem(
+        [connect(P_inner.output, :y, feedback.input2)
+         connect(feedback.output, :u, P_inner.input)
+         connect(ref.output, :r, feedback.input1)],
+        t,
+        systems = [P_inner, feedback, ref])
+
+    P_not_broken, _ = linearize(sys_inner, AnalysisPoint(:u), AnalysisPoint(:y))
+    @test P_not_broken.A[] == -2
+    P_broken, _ = linearize(sys_inner, AnalysisPoint(:u), AnalysisPoint(:y),
+        loop_openings = [AnalysisPoint(:u)])
+    @test P_broken.A[] == -1
+    P_broken, _ = linearize(sys_inner, AnalysisPoint(:u), AnalysisPoint(:y),
+        loop_openings = [AnalysisPoint(:y)])
+    @test P_broken.A[] == -1
+
+    Sinner = sminreal(ss(get_sensitivity(sys_inner, :u)[1]...))
+
+    @named sys_inner = ODESystem(
+        [connect(P_inner.output, :y, feedback.input2)
+         connect(feedback.output, :u, P_inner.input)],
+        t,
+        systems = [P_inner, feedback])
+
+    @named P_outer = FirstOrder(k = rand(), T = rand())
+
+    @named sys_outer = ODESystem(
+        [connect(sys_inner.P_inner.output, :y2, P_outer.input)
+         connect(P_outer.output, :u2, sys_inner.feedback.input1)],
+        t,
+        systems = [P_outer, sys_inner])
+
+    Souter = sminreal(ss(get_sensitivity(sys_outer, :sys_inner_u)[1]...))
+
+    Sinner2 = sminreal(ss(get_sensitivity(
+        sys_outer, :sys_inner_u, loop_openings = [:y2])[1]...))
+
+    @test Sinner.nx == 1
+    @test Sinner == Sinner2
+    @test Souter.nx == 2
+end
+
+@testset "sensitivities in multivariate signals" begin
+    A = [-0.994 -0.0794; -0.006242 -0.0134]
+    B = [-0.181 -0.389; 1.1 1.12]
+    C = [1.74 0.72; -0.33 0.33]
+    D = [0.0 0.0; 0.0 0.0]
+    @named P = Blocks.StateSpace(A, B, C, D)
+    Pss = CS.ss(A, B, C, D)
+
+    A = [-0.097;;]
+    B = [-0.138 -1.02]
+    C = [-0.076; 0.09;;]
+    D = [0.0 0.0; 0.0 0.0]
+    @named K = Blocks.StateSpace(A, B, C, D)
+    Kss = CS.ss(A, B, C, D)
+
+    eqs = [connect(P.output, :plant_output, K.input)
+           connect(K.output, :plant_input, P.input)]
+    sys = ODESystem(eqs, t, systems = [P, K], name = :hej)
+
+    matrices, _ = Blocks.get_sensitivity(sys, :plant_input)
+    S = CS.feedback(I(2), Kss * Pss, pos_feedback = true)
+
+    @test CS.tf(CS.ss(matrices...)) ≈ CS.tf(S)
+
+    matrices, _ = Blocks.get_comp_sensitivity(sys, :plant_input)
+    T = -CS.feedback(Kss * Pss, I(2), pos_feedback = true)
+
+    # bodeplot([ss(matrices...), T])
+    @test CS.tf(CS.ss(matrices...)) ≈ CS.tf(T)
+
+    matrices, _ = Blocks.get_looptransfer(
+        sys, :plant_input)
+    L = Kss * Pss
+    @test CS.tf(CS.ss(matrices...)) ≈ CS.tf(L)
+
+    matrices, _ = linearize(sys, :plant_input, :plant_output)
+    G = CS.feedback(Pss, Kss, pos_feedback = true)
+    @test CS.tf(CS.ss(matrices...)) ≈ CS.tf(G)
+end
+
+@testset "multiple analysis points" begin
+    @named P = FirstOrder(k = 1, T = 1)
+    @named C = Gain(; k = 1)
+    @named add = Blocks.Add(k2 = -1)
+
+    eqs = [connect(P.output, :plant_output, add.input2)
+           connect(add.output, C.input)
+           connect(C.output, :plant_input, P.input)]
+
+    sys_inner = ODESystem(eqs, t, systems = [P, C, add], name = :inner)
+
+    @named r = Constant(k = 1)
+    @named F = FirstOrder(k = 1, T = 3)
+
+    eqs = [connect(r.output, F.input)
+           connect(F.output, sys_inner.add.input1)]
+    sys_outer = ODESystem(eqs, t, systems = [F, sys_inner, r], name = :outer)
+
+    matrices, _ = get_sensitivity(sys_outer, [:, :inner_plant_output])
+
+    Ps = tf(1, [1, 1]) |> ss
+    Cs = tf(1) |> ss
+
+    G = CS.ss(matrices...) |> sminreal
+    Si = CS.feedback(1, Cs * Ps)
+    @test tf(G[1, 1]) ≈ tf(Si)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index af83ffccd8..d38aa86bb2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -111,6 +111,7 @@ end
         @safetestset "Linearization Tests" include("downstream/linearize.jl")
         @safetestset "Linearization Dummy Derivative Tests" include("downstream/linearization_dd.jl")
         @safetestset "Inverse Models Test" include("downstream/inversemodel.jl")
+        @safetestset "Analysis Points Test" include("downstream/analysis_points.jl")
     end
 
     if GROUP == "All" || GROUP == "Extensions"

From 469eda8837e7f7db533763e108e578fb9938bfe9 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 2 Jan 2025 23:33:41 +0530
Subject: [PATCH 08/23] feat: handle array variables as inputs/outputs of
 `linearization_function`

---
 src/systems/abstractsystem.jl | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 9b5ff7183d..3790cb8eb0 100644
--- a/src/systems/abstractsystem.jl
+++ b/src/systems/abstractsystem.jl
@@ -2374,6 +2374,12 @@ function linearization_function(sys::AbstractSystem, inputs,
     op = Dict(op)
     inputs isa AbstractVector || (inputs = [inputs])
     outputs isa AbstractVector || (outputs = [outputs])
+    inputs = mapreduce(vcat, inputs; init = []) do var
+        symbolic_type(var) == ArraySymbolic() ? collect(var) : [var]
+    end
+    outputs = mapreduce(vcat, outputs; init = []) do var
+        symbolic_type(var) == ArraySymbolic() ? collect(var) : [var]
+    end
     ssys, diff_idxs, alge_idxs, input_idxs = io_preprocessing(sys, inputs, outputs;
         simplify,
         kwargs...)

From ccfed451bfb8ba124fc4564eb27d3fd5111642f7 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 2 Jan 2025 23:34:02 +0530
Subject: [PATCH 09/23] fix: fix array variable handling in
 `get_analysis_variable`

---
 src/systems/analysis_points.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
index 605eb32749..ef20377147 100644
--- a/src/systems/analysis_points.jl
+++ b/src/systems/analysis_points.jl
@@ -255,11 +255,11 @@ function get_analysis_variable(var, name, iv; perturb = true)
     if perturb
         name = Symbol(:d_, name)
     end
-    if Symbolics.isarraysymbolic(var)
+    if symbolic_type(var) == ArraySymbolic()
         T = Array{eltype(symtype(var)), ndims(var)}
         pvar = unwrap(only(@variables $name(iv)::T))
         pvar = setmetadata(pvar, Symbolics.ArrayShapeCtx, Symbolics.shape(var))
-        default = zeros(symtype(var), size(var))
+        default = zeros(eltype(symtype(var)), size(var))
     else
         T = symtype(var)
         pvar = unwrap(only(@variables $name(iv)::T))
@@ -429,7 +429,7 @@ function apply_transformation(tf::PerturbOutput, sys::AbstractSystem)
         ap_ivar = ap_var(ap.input)
         new_var, new_def = get_analysis_variable(ap_ivar, nameof(ap), get_iv(sys))
         for outsys in ap.outputs
-            push!(ap_sys_eqs, ap_var(outsys) ~ ap_ivar + new_var)
+            push!(ap_sys_eqs, ap_var(outsys) ~ ap_ivar + wrap(new_var))
         end
         # add variable
         unks = copy(get_unknowns(ap_sys))
@@ -446,7 +446,7 @@ function apply_transformation(tf::PerturbOutput, sys::AbstractSystem)
         out_var, out_def = get_analysis_variable(
             ap_ivar, nameof(ap), get_iv(sys); perturb = false)
         defs[out_var] = out_def
-        push!(ap_sys_eqs, out_var ~ ap_ivar + new_var)
+        push!(ap_sys_eqs, out_var ~ ap_ivar + wrap(new_var))
         push!(unks, out_var)
 
         return ap_sys, (new_var, out_var)
@@ -526,11 +526,11 @@ function apply_transformation(cst::ComplementarySensitivityTransform, sys::Abstr
     # but comp sensitivity wants `output + du ~ input`. Thus, `du ~ -_du`.
     eqs = copy(get_eqs(sys))
     @set! sys.eqs = eqs
-    push!(eqs, du ~ -_du)
+    push!(eqs, du ~ -wrap(_du))
 
     defs = copy(get_defaults(sys))
     @set! sys.defaults = defs
-    defs[du] = -_du
+    defs[du] = -wrap(_du)
     return sys, (du, u)
 end
 

From f9b62e2b3b16cf019fc355bbeec231e951e54409 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 2 Jan 2025 23:34:59 +0530
Subject: [PATCH 10/23] fix: fix naming of new variable in
 `ComplementarySensitivityTransform`

---
 src/systems/analysis_points.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
index ef20377147..e71e67c691 100644
--- a/src/systems/analysis_points.jl
+++ b/src/systems/analysis_points.jl
@@ -519,7 +519,10 @@ end
 
 function apply_transformation(cst::ComplementarySensitivityTransform, sys::AbstractSystem)
     sys, (u,) = apply_transformation(GetInput(cst.ap), sys)
-    sys, (du,) = apply_transformation(AddVariable(cst.ap, Symbol(:comp_sens_du)), sys)
+    sys, (du,) = apply_transformation(
+        AddVariable(
+            cst.ap, Symbol(namespace_hierarchy(nameof(cst.ap))[end], :_comp_sens_du)),
+        sys)
     sys, (_du,) = apply_transformation(PerturbOutput(cst.ap), sys)
 
     # `PerturbOutput` adds the equation `input + _du ~ output`

From ed5fa70d33f4b3c8e950a85bcddd71b63f3695ed Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 2 Jan 2025 23:35:41 +0530
Subject: [PATCH 11/23] fix: handle `loop_openings`, `system_modifier` kwargs
 in `get_*` linear analysis functions

---
 src/systems/analysis_points.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
index e71e67c691..b8ec48eb68 100644
--- a/src/systems/analysis_points.jl
+++ b/src/systems/analysis_points.jl
@@ -611,8 +611,10 @@ function get_looptransfer_function(
 end
 
 for f in [:get_sensitivity, :get_comp_sensitivity, :get_looptransfer]
-    @eval function $f(sys, ap, args...; kwargs...)
-        lin_fun, ssys = $(Symbol(f, :_function))(sys, ap, args...; kwargs...)
+    @eval function $f(
+            sys, ap, args...; loop_openings = [], system_modifier = identity, kwargs...)
+        lin_fun, ssys = $(Symbol(f, :_function))(
+            sys, ap, args...; loop_openings, system_modifier, kwargs...)
         ModelingToolkit.linearize(ssys, lin_fun; kwargs...), ssys
     end
 end

From c8a4127b6580c089f12e51a435268f5b62fb869f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 2 Jan 2025 23:36:09 +0530
Subject: [PATCH 12/23] fix: fix handling of `loop_openings` in
 `linearization_function`

---
 src/systems/analysis_points.jl | 20 ++++++++++++++------
 1 file changed, 14 insertions(+), 6 deletions(-)

diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
index b8ec48eb68..180b730a98 100644
--- a/src/systems/analysis_points.jl
+++ b/src/systems/analysis_points.jl
@@ -630,25 +630,33 @@ end
 function linearization_function(sys::AbstractSystem,
         inputs::Union{Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}},
         outputs; loop_openings = [], system_modifier = identity, kwargs...)
-    sys = handle_loop_openings(sys, loop_openings)
-
+    loop_openings = Set(map(nameof, canonicalize_ap(loop_openings)))
     inputs = canonicalize_ap(inputs)
     outputs = canonicalize_ap(outputs)
 
     input_vars = []
     for input in inputs
-        sys, (input_var,) = apply_transformation(PerturbOutput(input), sys)
+        if nameof(input) in loop_openings
+            delete!(loop_openings, nameof(input))
+            sys, (input_var,) = apply_transformation(Break(input, true), sys)
+        else
+            sys, (input_var,) = apply_transformation(PerturbOutput(input), sys)
+        end
         push!(input_vars, input_var)
     end
     output_vars = []
     for output in outputs
         if output isa AnalysisPoint
-            sys, (output_var,) = apply_transformation(GetInput(output), sys)
-            push!(output_vars, output_var)
+            sys, (output_var,) = apply_transformation(AddVariable(output), sys)
+            sys, (input_var,) = apply_transformation(GetInput(output), sys)
+            push!(get_eqs(sys), output_var ~ input_var)
         else
-            push!(output_vars, output)
+            output_var = output
         end
+        push!(output_vars, output_var)
     end
 
+    sys = handle_loop_openings(sys, collect(loop_openings))
+
     return linearization_function(system_modifier(sys), input_vars, output_vars; kwargs...)
 end

From 0610d67ce61fc08f9fbb9bdb4f9151826144da82 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 2 Jan 2025 23:36:35 +0530
Subject: [PATCH 13/23] test: fix analysis point tests

---
 test/analysis_points.jl | 72 ++---------------------------------------
 1 file changed, 3 insertions(+), 69 deletions(-)

diff --git a/test/analysis_points.jl b/test/analysis_points.jl
index 6eae30314e..c888261db2 100644
--- a/test/analysis_points.jl
+++ b/test/analysis_points.jl
@@ -2,7 +2,7 @@ using ModelingToolkit, ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEq, LinearAlgebra
 using Test
 using ModelingToolkit: t_nounits as t, D_nounits as D, AnalysisPoint, get_sensitivity,
-                       get_comp_sensitivity, get_looptransfer
+                       get_comp_sensitivity, get_looptransfer, open_loop, AbstractSystem
 using Symbolics: NAMESPACE_SEPARATOR
 
 @testset "AnalysisPoint is lowered to `connect`" begin
@@ -13,14 +13,14 @@ using Symbolics: NAMESPACE_SEPARATOR
     eqs = [connect(P.output, C.input)
            connect(C.output, ap, P.input)]
     sys_ap = ODESystem(eqs, t, systems = [P, C], name = :hej)
-    sys_ap2 = @test_nowarn expand_connections(sys)
+    sys_ap2 = @test_nowarn expand_connections(sys_ap)
 
     @test all(eq -> !(eq.lhs isa AnalysisPoint), equations(sys_ap2))
 
     eqs = [connect(P.output, C.input)
            connect(C.output, P.input)]
     sys_normal = ODESystem(eqs, t, systems = [P, C], name = :hej)
-    sys_normal2 = @test_nowarn expand_connections(sys)
+    sys_normal2 = @test_nowarn expand_connections(sys_normal)
 
     @test isequal(sys_ap2, sys_normal2)
 end
@@ -191,69 +191,3 @@ end
     @test matrices.B[] * matrices.C[] == 1 # both positive
     @test matrices.D[] == 0
 end
-
-using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra
-using ModelingToolkitStandardLibrary.Mechanical.Rotational
-using ModelingToolkitStandardLibrary.Blocks: Sine, PID, SecondOrder, Step, RealOutput
-using ModelingToolkit: connect
-
-@testset "Complicated model" begin
-    # Parameters
-    m1 = 1
-    m2 = 1
-    k = 1000 # Spring stiffness
-    c = 10   # Damping coefficient
-    @named inertia1 = Inertia(; J = m1)
-    @named inertia2 = Inertia(; J = m2)
-    @named spring = Spring(; c = k)
-    @named damper = Damper(; d = c)
-    @named torque = Torque()
-
-    function SystemModel(u = nothing; name = :model)
-        eqs = [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)]
-        if u !== nothing
-            push!(eqs, connect(torque.tau, u.output))
-            return ODESystem(eqs, t;
-                systems = [
-                    torque,
-                    inertia1,
-                    inertia2,
-                    spring,
-                    damper,
-                    u
-                ],
-                name)
-        end
-        ODESystem(eqs, t; systems = [torque, inertia1, inertia2, spring, damper], name)
-    end
-
-    @named r = Step(start_time = 0)
-    model = SystemModel()
-    @named pid = PID(k = 100, Ti = 0.5, Td = 1)
-    @named filt = SecondOrder(d = 0.9, w = 10)
-    @named sensor = AngleSensor()
-    @named er = Add(k2 = -1)
-
-    connections = [connect(r.output, :r, filt.input)
-                   connect(filt.output, er.input1)
-                   connect(pid.ctr_output, :u, model.torque.tau)
-                   connect(model.inertia2.flange_b, sensor.flange)
-                   connect(sensor.phi, :y, er.input2)
-                   connect(er.output, :e, pid.err_input)]
-
-    closed_loop = ODESystem(connections, t, systems = [model, pid, filt, sensor, r, er],
-        name = :closed_loop, defaults = [
-            model.inertia1.phi => 0.0,
-            model.inertia2.phi => 0.0,
-            model.inertia1.w => 0.0,
-            model.inertia2.w => 0.0,
-            filt.x => 0.0,
-            filt.xd => 0.0
-        ])
-
-    sys = structural_simplify(closed_loop)
-    prob = ODEProblem(sys, unknowns(sys) .=> 0.0, (0.0, 4.0))
-    sol = solve(prob, Rodas5P(), reltol = 1e-6, abstol = 1e-9)
-end

From bb23e6216360b625314e5491e28196220ffdc900 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Thu, 2 Jan 2025 23:36:48 +0530
Subject: [PATCH 14/23] test: fix and port over remaining analysis point
 downstream tests

---
 test/downstream/analysis_points.jl | 58 +++++++++++++++++++++++++-----
 1 file changed, 50 insertions(+), 8 deletions(-)

diff --git a/test/downstream/analysis_points.jl b/test/downstream/analysis_points.jl
index 5dc2b8957c..85d271a547 100644
--- a/test/downstream/analysis_points.jl
+++ b/test/downstream/analysis_points.jl
@@ -1,7 +1,8 @@
 using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, ControlSystemsBase
 using ModelingToolkitStandardLibrary.Mechanical.Rotational
 using ModelingToolkitStandardLibrary.Blocks
-using ModelingToolkit: connect, AnalysisPoint, t_nounits as t, D_nounits as D
+using ModelingToolkit: connect, AnalysisPoint, t_nounits as t, D_nounits as D,
+                       get_sensitivity, get_comp_sensitivity, get_looptransfer, open_loop
 import ControlSystemsBase as CS
 
 @testset "Complicated model" begin
@@ -154,10 +155,10 @@ end
         t,
         systems = [P_outer, sys_inner])
 
-    Souter = sminreal(ss(get_sensitivity(sys_outer, :sys_inner_u)[1]...))
+    Souter = sminreal(ss(get_sensitivity(sys_outer, sys_inner.u)[1]...))
 
     Sinner2 = sminreal(ss(get_sensitivity(
-        sys_outer, :sys_inner_u, loop_openings = [:y2])[1]...))
+        sys_outer, sys_inner.u, loop_openings = [:y2])[1]...))
 
     @test Sinner.nx == 1
     @test Sinner == Sinner2
@@ -183,23 +184,23 @@ end
            connect(K.output, :plant_input, P.input)]
     sys = ODESystem(eqs, t, systems = [P, K], name = :hej)
 
-    matrices, _ = Blocks.get_sensitivity(sys, :plant_input)
+    matrices, _ = get_sensitivity(sys, :plant_input)
     S = CS.feedback(I(2), Kss * Pss, pos_feedback = true)
 
     @test CS.tf(CS.ss(matrices...)) ≈ CS.tf(S)
 
-    matrices, _ = Blocks.get_comp_sensitivity(sys, :plant_input)
+    matrices, _ = get_comp_sensitivity(sys, :plant_input)
     T = -CS.feedback(Kss * Pss, I(2), pos_feedback = true)
 
     # bodeplot([ss(matrices...), T])
     @test CS.tf(CS.ss(matrices...)) ≈ CS.tf(T)
 
-    matrices, _ = Blocks.get_looptransfer(
+    matrices, _ = get_looptransfer(
         sys, :plant_input)
     L = Kss * Pss
     @test CS.tf(CS.ss(matrices...)) ≈ CS.tf(L)
 
-    matrices, _ = linearize(sys, :plant_input, :plant_output)
+    matrices, _ = linearize(sys, AnalysisPoint(:plant_input), :plant_output)
     G = CS.feedback(Pss, Kss, pos_feedback = true)
     @test CS.tf(CS.ss(matrices...)) ≈ CS.tf(G)
 end
@@ -222,7 +223,8 @@ end
            connect(F.output, sys_inner.add.input1)]
     sys_outer = ODESystem(eqs, t, systems = [F, sys_inner, r], name = :outer)
 
-    matrices, _ = get_sensitivity(sys_outer, [:, :inner_plant_output])
+    matrices, _ = get_sensitivity(
+        sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output])
 
     Ps = tf(1, [1, 1]) |> ss
     Cs = tf(1) |> ss
@@ -230,4 +232,44 @@ end
     G = CS.ss(matrices...) |> sminreal
     Si = CS.feedback(1, Cs * Ps)
     @test tf(G[1, 1]) ≈ tf(Si)
+
+    So = CS.feedback(1, Ps * Cs)
+    @test tf(G[2, 2]) ≈ tf(So)
+    @test tf(G[1, 2]) ≈ tf(-CS.feedback(Cs, Ps))
+    @test tf(G[2, 1]) ≈ tf(CS.feedback(Ps, Cs))
+
+    matrices, _ = get_comp_sensitivity(
+        sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output])
+
+    G = CS.ss(matrices...) |> sminreal
+    Ti = CS.feedback(Cs * Ps)
+    @test tf(G[1, 1]) ≈ tf(Ti)
+
+    To = CS.feedback(Ps * Cs)
+    @test tf(G[2, 2]) ≈ tf(To)
+    @test tf(G[1, 2]) ≈ tf(CS.feedback(Cs, Ps)) # The negative sign appears in a confusing place due to negative feedback not happening through Ps
+    @test tf(G[2, 1]) ≈ tf(-CS.feedback(Ps, Cs))
+
+    # matrices, _ = get_looptransfer(sys_outer, [:inner_plant_input, :inner_plant_output])
+    matrices, _ = get_looptransfer(sys_outer, sys_outer.inner.plant_input)
+    L = CS.ss(matrices...) |> sminreal
+    @test tf(L) ≈ -tf(Cs * Ps)
+
+    matrices, _ = get_looptransfer(sys_outer, sys_outer.inner.plant_output)
+    L = CS.ss(matrices...) |> sminreal
+    @test tf(L[1, 1]) ≈ -tf(Ps * Cs)
+
+    # Calling looptransfer like below is not the intended way, but we can work out what it should return if we did so it remains a valid test
+    matrices, _ = get_looptransfer(
+        sys_outer, [sys_outer.inner.plant_input, sys_outer.inner.plant_output])
+    L = CS.ss(matrices...) |> sminreal
+    @test tf(L[1, 1]) ≈ tf(0)
+    @test tf(L[2, 2]) ≈ tf(0)
+    @test sminreal(L[1, 2]) ≈ ss(-1)
+    @test tf(L[2, 1]) ≈ tf(Ps)
+
+    matrices, _ = linearize(
+        sys_outer, [sys_outer.inner.plant_input], [sys_inner.plant_output])
+    G = CS.ss(matrices...) |> sminreal
+    @test tf(G) ≈ tf(CS.feedback(Ps, Cs))
 end

From 89ed0cf57e6187c0a68c41a68f899aebf5eee516 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 3 Jan 2025 00:27:01 +0530
Subject: [PATCH 15/23] docs: add and update docstrings for analysis points and
 transformations

---
 src/systems/analysis_points.jl | 378 +++++++++++++++++++++++++++------
 1 file changed, 313 insertions(+), 65 deletions(-)

diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
index 180b730a98..9ff06109e5 100644
--- a/src/systems/analysis_points.jl
+++ b/src/systems/analysis_points.jl
@@ -1,25 +1,85 @@
 """
     $(TYPEDEF)
-    $(TYPEDSIGNATURES)
+    AnalysisPoint(input, name::Symbol, outputs::Vector)
 
 Create an AnalysisPoint for linear analysis. Analysis points can be created by calling
 
 ```
-connect(in, :ap_name, out...)
+connect(out, :ap_name, in...)
 ```
 
-Where `in` is the input to the connection, and `out...` are the outputs. In the context of
-ModelingToolkitStandardLibrary.jl, `in` is a `RealOutput` connector and `out...` are all
-`RealInput` connectors. All involved connectors are required to either have an unknown named
+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`](@ref), [`get_comp_sensitivity`](@ref), [`get_looptransfer`](@ref), [`open_loop`](@ref)
+
+# Fields
+
+$(TYPEDFIELDS)
+
+# Example
+
+```julia
+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
+
+```julia
+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
+
+```julia
+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)
+```
 """
 struct AnalysisPoint
+    """
+    The input to the connection. In the context of ModelingToolkitStandardLibrary.jl,
+    this is a `RealOutput` connector.
+    """
     input::Any
+    """
+    The name of the analysis point.
+    """
     name::Symbol
+    """
+    The outputs of the connection. In the context of ModelingToolkitStandardLibrary.jl,
+    these are all `RealInput` connectors.
+    """
     outputs::Union{Nothing, Vector{Any}}
 end
 
 AnalysisPoint() = AnalysisPoint(nothing, Symbol(), nothing)
+"""
+    $(TYPEDSIGNATURES)
+
+Create an `AnalysisPoint` with the given name, with no input or outputs specified.
+"""
 AnalysisPoint(name::Symbol) = AnalysisPoint(nothing, name, nothing)
 
 Base.nameof(ap::AnalysisPoint) = ap.name
@@ -78,9 +138,35 @@ function Symbolics.connect(in, ap::AnalysisPoint, outs...)
 end
 
 """
-    $(TYPEDSIGNATURES)
+    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`](@ref) 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
 
-Create an `AnalysisPoint` connection connecting `in` to `outs...`.
+```julia
+@named P = FirstOrder(k = 1, T = 1)
+@named C = Gain(; k = -1)
+connect(C.output, :plant_input, P.input)
+```
+
+is correct, whereas
+
+```julia
+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`](@ref)
+  - `ap_name`: If a name is given, an [`AnalysisPoint`](@ref) with the given name will be created automatically.
 """
 function Symbolics.connect(in::AbstractSystem, name::Symbol, out, outs...)
     return AnalysisPoint() ~ AnalysisPoint(in, name, [out; collect(outs)])
@@ -198,19 +284,23 @@ function modify_nested_subsystem(
 
     # recursive helper function which does the searching and modification
     function _helper(sys::AbstractSystem, i::Int)
-        # we reached past the end, so everything matched and
-        # `sys` is the system to modify.
         if i > length(hierarchy)
+            # we reached past the end, so everything matched and
+            # `sys` is the system to modify.
             sys, vars = fn(sys)
         else
+            # find the subsystem with the given name and error otherwise
             cur = hierarchy[i]
             idx = findfirst(subsys -> nameof(subsys) == cur, get_systems(sys))
             idx === nothing &&
                 error("System $(join([nameof(root); hierarchy[1:i-1]], '.')) does not have a subsystem named $cur.")
 
+            # recurse into new subsystem
             newsys, vars = _helper(get_systems(sys)[idx], i + 1)
+            # update this system with modified subsystem
             @set! sys.systems[idx] = newsys
         end
+        # only namespace variables from inner systems
         if i != 1
             vars = ntuple(Val(length(vars))) do i
                 renamespace(sys, vars[i])
@@ -270,6 +360,15 @@ end
 
 #### PRIMITIVE TRANSFORMATIONS
 
+const DOC_WILL_REMOVE_AP = """
+    Note that this transformation will remove `ap`, causing any subsequent transformations \
+    referring to it to fail.\
+    """
+
+const DOC_ADDED_VARIABLE = """
+    The added variable(s) will have a default of zero, of the appropriate type and size.\
+    """
+
 """
     $(TYPEDEF)
 
@@ -278,11 +377,9 @@ it will add a new input variable which connects to the outputs of the analysis p
 `apply_transformation` returns the new input variable (if added) as the auxiliary
 information. The new input variable will have the name `Symbol(:d_, nameof(ap))`.
 
-Note that this transformation will remove `ap`, causing any subsequent transformations
-referring to it to fail.
+$DOC_WILL_REMOVE_AP
 
-The added variable, if present, will have a default of zero, of the appropriate type and
-size.
+$DOC_ADDED_VARIABLE
 
 ## Fields
 
@@ -340,9 +437,7 @@ end
     $(TYPEDEF)
 
 A transformation which returns the variable corresponding to the input of the analysis
-point.
-
-`apply_transformation` returns the variable as auxiliary information.
+point. Does not modify the system.
 
 ## Fields
 
@@ -350,7 +445,7 @@ $(TYPEDFIELDS)
 """
 struct GetInput <: AnalysisPointTransformation
     """
-    The analysis point to add the output to.
+    The analysis point to get the input of.
     """
     ap::AnalysisPoint
 end
@@ -380,15 +475,12 @@ the analysis point before connecting to the outputs. The new variable will have
 
 If `with_output == true`, also creates an additional new variable which has the value
 provided to the outputs after the above modification. This new variable has the same name
-as the analysis point.
+as the analysis point and will be the second variable in the tuple of new variables returned
+from `apply_transformation`.
 
-`apply_transformation` returns a 1-tuple of the perturbation variable if
-`with_output == false` and a 2-tuple of the perturbation variable and output variable if
-`with_output == true`.
+$DOC_WILL_REMOVE_AP
 
-Removes the analysis point `ap`, so any subsequent transformations requiring it will fail.
-
-The added variable(s) will have a default of zero, of the appropriate type and size.
+$DOC_ADDED_VARIABLE
 
 ## Fields
 
@@ -454,17 +546,33 @@ function apply_transformation(tf::PerturbOutput, sys::AbstractSystem)
 end
 
 """
-    $(TYPEDSIGNATURES)
+    $(TYPEDEF)
 
 A transformation which adds a variable named `name` to the system containing the analysis
-point `ap`. The added variable has the same type and size as the input of the analysis
-point.
+point `ap`. $DOC_ADDED_VARIABLE
+
+# Fields
+
+$(TYPEDFIELDS)
 """
 struct AddVariable <: AnalysisPointTransformation
+    """
+    The analysis point in the system to modify, and whose input should be used as the
+    template for the new variable.
+    """
     ap::AnalysisPoint
+    """
+    The name of the added variable.
+    """
     name::Symbol
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Add a new variable to the system containing analysis point `ap` with the same name as the
+analysis point.
+"""
 AddVariable(ap::AnalysisPoint) = AddVariable(ap, nameof(ap))
 
 function apply_transformation(tf::AddVariable, sys::AbstractSystem)
@@ -494,11 +602,12 @@ end
     $(TYPEDSIGNATURES)
 
 A transformation enable calculating the sensitivity function about the analysis point `ap`.
-`apply_transformation` returns a 2-tuple `du, u` as auxiliary information.
+The returned added variables are `(du, u)` where `du` is the perturbation added to the
+input, and `u` is the output after perturbation.
 
-Removes the analysis point `ap`, so any subsequent transformations requiring it will fail.
+$DOC_WILL_REMOVE_AP
 
-The added variables will have a default of zero, of the appropriate type and size.
+$DOC_ADDED_VARIABLE
 """
 SensitivityTransform(ap::AnalysisPoint) = PerturbOutput(ap, true)
 
@@ -506,14 +615,21 @@ SensitivityTransform(ap::AnalysisPoint) = PerturbOutput(ap, true)
     $(TYPEDEF)
 
 A transformation to enable calculating the complementary sensitivity function about the
-analysis point `ap`. `apply_transformation` returns a 2-tuple `du, u` as auxiliary
-information.
+analysis point `ap`. The returned added variables are `(du, u)` where `du` is the
+perturbation added to the outputs and `u` is the input to the analysis point.
+
+$DOC_WILL_REMOVE_AP
+
+$DOC_ADDED_VARIABLE
 
-Removes the analysis point `ap`, so any subsequent transformations requiring it will fail.
+# Fields
 
-The added variables will have a default of zero, of the appropriate type and size.
+$(TYPEDFIELDS)
 """
 struct ComplementarySensitivityTransform <: AnalysisPointTransformation
+    """
+    The analysis point to modify.
+    """
     ap::AnalysisPoint
 end
 
@@ -537,7 +653,25 @@ function apply_transformation(cst::ComplementarySensitivityTransform, sys::Abstr
     return sys, (du, u)
 end
 
+"""
+    $(TYPEDEF)
+
+A transformation to enable calculating the loop transfer function about the analysis point
+`ap`. The returned added variables are `(du, u)` where `du` feeds into the outputs of `ap`
+and `u` is the input of `ap`.
+
+$DOC_WILL_REMOVE_AP
+
+$DOC_ADDED_VARIABLE
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
 struct LoopTransferTransform <: AnalysisPointTransformation
+    """
+    The analysis point to modify.
+    """
     ap::AnalysisPoint
 end
 
@@ -547,8 +681,13 @@ function apply_transformation(tf::LoopTransferTransform, sys::AbstractSystem)
     return sys, (du, u)
 end
 
-### TODO: Move these
+"""
+    $(TYPEDSIGNATURES)
 
+A utility function to get the "canonical" form of a list of analysis points. Always returns
+a list of values. Any value that cannot be turned into an `AnalysisPoint` (i.e. isn't
+already an `AnalysisPoint` or `Symbol`) is simply wrapped in an array.
+"""
 canonicalize_ap(ap::Symbol) = [AnalysisPoint(ap)]
 canonicalize_ap(ap::AnalysisPoint) = [ap]
 canonicalize_ap(ap) = [ap]
@@ -556,6 +695,11 @@ function canonicalize_ap(aps::Vector)
     mapreduce(canonicalize_ap, vcat, aps; init = [])
 end
 
+"""
+    $(TYPEDSIGNATURES)
+
+Given a list of analysis points, break the connection for each and set the output to zero.
+"""
 function handle_loop_openings(sys::AbstractSystem, aps)
     for ap in canonicalize_ap(aps)
         sys, (outvar,) = apply_transformation(Break(ap, true), sys)
@@ -568,63 +712,124 @@ function handle_loop_openings(sys::AbstractSystem, aps)
     return sys
 end
 
-function get_sensitivity_function(
-        sys::AbstractSystem, aps; system_modifier = identity, loop_openings = [], kwargs...)
+const DOC_LOOP_OPENINGS = """
+    - `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.
+"""
+
+const DOC_SYS_MODIFIER = """
+    - `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`. 
+"""
+"""
+    $(TYPEDSIGNATURES)
+
+Utility function for linear analyses that apply a transformation `transform`, which
+returns the added variables `(du, u)`, to each of the analysis points in `aps` and then
+calls `linearization_function` with all the `du`s as inputs and `u`s as outputs. Returns
+the linearization function and modified, simplified system.
+
+# Keyword arguments
+
+$DOC_LOOP_OPENINGS
+$DOC_SYS_MODIFIER
+
+All other keyword arguments are forwarded to `linearization_function`.
+"""
+function get_linear_analysis_function(
+        sys::AbstractSystem, transform, aps; system_modifier = identity, loop_openings = [], kwargs...)
     sys = handle_loop_openings(sys, loop_openings)
     aps = canonicalize_ap(aps)
     dus = []
     us = []
     for ap in aps
-        sys, (du, u) = apply_transformation(SensitivityTransform(ap), sys)
+        sys, (du, u) = apply_transformation(transform(ap), sys)
         push!(dus, du)
         push!(us, u)
     end
     linearization_function(system_modifier(sys), dus, us; kwargs...)
 end
 
-function get_comp_sensitivity_function(
-        sys::AbstractSystem, aps; system_modifier = identity, loop_openings = [], kwargs...)
-    sys = handle_loop_openings(sys, loop_openings)
-    aps = canonicalize_ap(aps)
-    dus = []
-    us = []
-    for ap in aps
-        sys, (du, u) = apply_transformation(ComplementarySensitivityTransform(ap), sys)
-        push!(dus, du)
-        push!(us, u)
-    end
-    linearization_function(system_modifier(sys), dus, us; kwargs...)
+"""
+    $(TYPEDSIGNATURES)
+
+Return the sensitivity function for the analysis point(s) `aps`, and the modified system
+simplified with the appropriate inputs and outputs.
+
+# Keyword Arguments
+
+$DOC_LOOP_OPENINGS
+$DOC_SYS_MODIFIER
+
+All other keyword arguments are forwarded to `linearization_function`.
+"""
+function get_sensitivity_function(sys::AbstractSystem, aps; kwargs...)
+    get_linear_analysis_function(sys, SensitivityTransform, aps; kwargs...)
 end
 
-function get_looptransfer_function(
-        sys, aps; system_modifier = identity, loop_openings = [], kwargs...)
-    sys = handle_loop_openings(sys, loop_openings)
-    aps = canonicalize_ap(aps)
-    dus = []
-    us = []
-    for ap in aps
-        sys, (du, u) = apply_transformation(LoopTransferTransform(ap), sys)
-        push!(dus, du)
-        push!(us, u)
-    end
-    linearization_function(system_modifier(sys), dus, us; kwargs...)
+"""
+    $(TYPEDSIGNATURES)
+
+Return the complementary sensitivity function for the analysis point(s) `aps`, and the
+modified system simplified with the appropriate inputs and outputs.
+
+# Keyword Arguments
+
+$DOC_LOOP_OPENINGS
+$DOC_SYS_MODIFIER
+
+All other keyword arguments are forwarded to `linearization_function`.
+"""
+function get_comp_sensitivity_function(sys::AbstractSystem, aps; kwargs...)
+    get_linear_analysis_function(sys, ComplementarySensitivityTransform, aps; kwargs...)
+end
+
+"""
+    $(TYPEDSIGNATURES)
+
+Return the loop-transfer function for the analysis point(s) `aps`, and the modified
+system simplified with the appropriate inputs and outputs.
+
+# Keyword Arguments
+
+$DOC_LOOP_OPENINGS
+$DOC_SYS_MODIFIER
+
+All other keyword arguments are forwarded to `linearization_function`.
+"""
+function get_looptransfer_function(sys::AbstractSystem, aps; kwargs...)
+    get_linear_analysis_function(sys, LoopTransferTransform, aps; kwargs...)
 end
 
 for f in [:get_sensitivity, :get_comp_sensitivity, :get_looptransfer]
+    utility_fun = Symbol(f, :_function)
     @eval function $f(
             sys, ap, args...; loop_openings = [], system_modifier = identity, kwargs...)
-        lin_fun, ssys = $(Symbol(f, :_function))(
+        lin_fun, ssys = $(utility_fun)(
             sys, ap, args...; loop_openings, system_modifier, kwargs...)
         ModelingToolkit.linearize(ssys, lin_fun; kwargs...), ssys
     end
 end
 
-function open_loop(sys, ap::Union{Symbol, AnalysisPoint}; kwargs...)
+"""
+    $(TYPEDSIGNATURES)
+
+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 peformed.
+"""
+function open_loop(sys, ap::Union{Symbol, AnalysisPoint}; system_modifier = identity)
     if ap isa Symbol
         ap = AnalysisPoint(ap)
     end
     tf = LoopTransferTransform(ap)
-    return apply_transformation(tf, sys)
+    sys, vars = apply_transformation(tf, sys)
+    return system_modifier(sys), vars
 end
 
 function linearization_function(sys::AbstractSystem,
@@ -660,3 +865,46 @@ function linearization_function(sys::AbstractSystem,
 
     return linearization_function(system_modifier(sys), input_vars, output_vars; kwargs...)
 end
+
+@doc """
+        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`](@ref), [`get_looptransfer`](@ref).
+""" get_sensitivity
+
+@doc """
+    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`](@ref), [`get_looptransfer`](@ref).
+""" get_comp_sensitivity
+
+@doc """
+    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`.
+
+!!! info "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`](@ref), [`get_comp_sensitivity`](@ref), [`open_loop`](@ref).
+""" get_looptransfer

From ddceb6b6127c89e78308d04275c0e652b7bc87b3 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Fri, 3 Jan 2025 17:29:22 +0530
Subject: [PATCH 16/23] refactor: use version of Stdlib which disables old
 `AnalysisPoint`

---
 Project.toml                       |  1 +
 test/analysis_points.jl            | 15 ++-------------
 test/downstream/Project.toml       |  3 +++
 test/downstream/analysis_points.jl |  3 +--
 4 files changed, 7 insertions(+), 15 deletions(-)

diff --git a/Project.toml b/Project.toml
index e4eba3cfad..256b1de6fc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -117,6 +117,7 @@ Latexify = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16"
 Libdl = "1"
 LinearAlgebra = "1"
 MLStyle = "0.4.17"
+ModelingToolkitStandardLibrary = "2.19"
 NaNMath = "0.3, 1"
 NonlinearSolve = "4.3"
 OffsetArrays = "1"
diff --git a/test/analysis_points.jl b/test/analysis_points.jl
index c888261db2..9361b5a686 100644
--- a/test/analysis_points.jl
+++ b/test/analysis_points.jl
@@ -1,8 +1,7 @@
 using ModelingToolkit, ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEq, LinearAlgebra
 using Test
-using ModelingToolkit: t_nounits as t, D_nounits as D, AnalysisPoint, get_sensitivity,
-                       get_comp_sensitivity, get_looptransfer, open_loop, AbstractSystem
+using ModelingToolkit: t_nounits as t, D_nounits as D, AnalysisPoint, AbstractSystem
 using Symbolics: NAMESPACE_SEPARATOR
 
 @testset "AnalysisPoint is lowered to `connect`" begin
@@ -164,17 +163,7 @@ end
         (inputap, [nameof(outputap)]),
         (nameof(inputap), [nameof(outputap)])
     ]
-        if input isa Symbol
-            # broken because MTKStdlib defines the method for
-            # `input::Union{Symbol, Vector{Symbol}}` which we can't directly call
-            @test_broken linearize(sys, input, output)
-            linfun, ssys = @invoke linearization_function(sys::AbstractSystem,
-                input::Union{Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}},
-                output::Any)
-            matrices = linearize(ssys, linfun)
-        else
-            matrices, _ = linearize(sys, input, output)
-        end
+        matrices, _ = linearize(sys, input, output)
         # Result should be the same as feedpack(P, 1), i.e., the closed-loop transfer function from plant input to plant output
         @test matrices.A[] == -2
         @test matrices.B[] * matrices.C[] == 1 # both positive
diff --git a/test/downstream/Project.toml b/test/downstream/Project.toml
index b8776e1e4f..ade09e797b 100644
--- a/test/downstream/Project.toml
+++ b/test/downstream/Project.toml
@@ -5,3 +5,6 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
+
+[compat]
+ModelingToolkitStandardLibrary = "2.19"
diff --git a/test/downstream/analysis_points.jl b/test/downstream/analysis_points.jl
index 85d271a547..fabca0f493 100644
--- a/test/downstream/analysis_points.jl
+++ b/test/downstream/analysis_points.jl
@@ -1,8 +1,7 @@
 using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, ControlSystemsBase
 using ModelingToolkitStandardLibrary.Mechanical.Rotational
 using ModelingToolkitStandardLibrary.Blocks
-using ModelingToolkit: connect, AnalysisPoint, t_nounits as t, D_nounits as D,
-                       get_sensitivity, get_comp_sensitivity, get_looptransfer, open_loop
+using ModelingToolkit: connect, t_nounits as t, D_nounits as D
 import ControlSystemsBase as CS
 
 @testset "Complicated model" begin

From de17b878e35ca7fb05c5af2a262f1c9ef799c51e Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 6 Jan 2025 12:31:29 +0530
Subject: [PATCH 17/23] feat: validate causality of analysis points

---
 src/systems/analysis_points.jl | 54 ++++++++++++++++++++++++++++------
 test/analysis_points.jl        | 10 +++++++
 2 files changed, 55 insertions(+), 9 deletions(-)

diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
index 9ff06109e5..a2b7fb7ea8 100644
--- a/src/systems/analysis_points.jl
+++ b/src/systems/analysis_points.jl
@@ -72,6 +72,36 @@ struct AnalysisPoint
     these are all `RealInput` connectors.
     """
     outputs::Union{Nothing, Vector{Any}}
+
+    function AnalysisPoint(input, name::Symbol, outputs; verbose = true)
+        # input to analysis point should be an output variable
+        if verbose && input !== nothing
+            var = ap_var(input)
+            isoutput(var) || ap_warning(1, name, true)
+        end
+        # outputs of analysis points should be input variables
+        if verbose && outputs !== nothing
+            for (i, output) in enumerate(outputs)
+                var = ap_var(output)
+                isinput(var) || ap_warning(2 + i, name, false)
+            end
+        end
+
+        return new(input, name, outputs)
+    end
+end
+
+function ap_warning(arg::Int, name::Symbol, should_be_output)
+    causality = should_be_output ? "output" : "input"
+    @warn """
+    The $(arg)-th argument to analysis point $(name) was not a $causality. This is supported in \
+    order to handle inverse models, but may not be what you intended.
+
+    If you are building a forward mode (causal), you may want to swap this argument with \
+    one on the opposite side of the name of the analysis point provided to `connect`. \
+    Learn more about the causality of analysis points in the docstring for `AnalysisPoint`. \
+    Silence this message using `connect(out, :name, in...; warn = false)`.
+    """
 end
 
 AnalysisPoint() = AnalysisPoint(nothing, Symbol(), nothing)
@@ -133,8 +163,8 @@ function renamespace(sys, ap::AnalysisPoint)
 end
 
 # create analysis points via `connect`
-function Symbolics.connect(in, ap::AnalysisPoint, outs...)
-    return AnalysisPoint() ~ AnalysisPoint(in, ap.name, collect(outs))
+function Symbolics.connect(in, ap::AnalysisPoint, outs...; verbose = true)
+    return AnalysisPoint() ~ AnalysisPoint(in, ap.name, collect(outs); verbose)
 end
 
 """
@@ -161,15 +191,21 @@ connect(P.input, :plant_input, C.output)
 
 typically is not (unless the model is an inverse model).
 
-# Arguments:
+# Arguments
+
+- `output_connector`: An output connector
+- `input_connector`: An input connector
+- `ap`: An explicitly created [`AnalysisPoint`](@ref)
+- `ap_name`: If a name is given, an [`AnalysisPoint`](@ref) with the given name will be
+  created automatically.
+
+# Keyword arguments
 
-  - `output_connector`: An output connector
-  - `input_connector`: An input connector
-  - `ap`: An explicitly created [`AnalysisPoint`](@ref)
-  - `ap_name`: If a name is given, an [`AnalysisPoint`](@ref) with the given name will be created automatically.
+- `verbose`: Warn if an input is connected to an output (reverse causality). Silence this
+  warning if you are analyzing an inverse model.
 """
-function Symbolics.connect(in::AbstractSystem, name::Symbol, out, outs...)
-    return AnalysisPoint() ~ AnalysisPoint(in, name, [out; collect(outs)])
+function Symbolics.connect(in::AbstractSystem, name::Symbol, out, outs...; verbose = true)
+    return AnalysisPoint() ~ AnalysisPoint(in, name, [out; collect(outs)]; verbose)
 end
 
 """
diff --git a/test/analysis_points.jl b/test/analysis_points.jl
index 9361b5a686..8c93e919b9 100644
--- a/test/analysis_points.jl
+++ b/test/analysis_points.jl
@@ -24,6 +24,16 @@ using Symbolics: NAMESPACE_SEPARATOR
     @test isequal(sys_ap2, sys_normal2)
 end
 
+@testset "Inverse causality throws a warning" begin
+    @named P = FirstOrder(k = 1, T = 1)
+    @named C = Gain(; k = -1)
+
+    ap = AnalysisPoint(:plant_input)
+    @test_warn ["1-th argument", "plant_input", "not a output"] connect(
+        P.input, ap, C.output)
+    @test_nowarn connect(P.input, ap, C.output; verbose = false)
+end
+
 # also tests `connect(input, name::Symbol, outputs...)` syntax
 @testset "AnalysisPoint is accessible via `getproperty`" begin
     @named P = FirstOrder(k = 1, T = 1)

From 50973353bb7da69dc35d838f90177f5ae1b461d1 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 6 Jan 2025 14:06:29 +0530
Subject: [PATCH 18/23] fix: fix `isvarkind` for `Symbolics.Arr`

---
 src/variables.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/variables.jl b/src/variables.jl
index 3b0d64e7ea..c45c4b0f00 100644
--- a/src/variables.jl
+++ b/src/variables.jl
@@ -90,7 +90,7 @@ struct Equality <: AbstractConnectType end # Equality connection
 struct Flow <: AbstractConnectType end     # sum to 0
 struct Stream <: AbstractConnectType end   # special stream connector
 
-isvarkind(m, x::Num) = isvarkind(m, value(x))
+isvarkind(m, x::Union{Num, Symbolics.Arr}) = isvarkind(m, value(x))
 function isvarkind(m, x)
     iskind = getmetadata(x, m, nothing)
     iskind !== nothing && return iskind

From 2b1f525cb036dadec10bb8c34503e480240e4d5c Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 6 Jan 2025 14:55:13 +0530
Subject: [PATCH 19/23] fix: remove ambiguities in namespacing of analysis
 points

---
 src/systems/analysis_points.jl     | 34 ++++++++---------
 test/analysis_points.jl            | 59 ++++++++++++------------------
 test/downstream/analysis_points.jl | 16 ++++----
 3 files changed, 48 insertions(+), 61 deletions(-)

diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
index a2b7fb7ea8..87f278d4dc 100644
--- a/src/systems/analysis_points.jl
+++ b/src/systems/analysis_points.jl
@@ -314,9 +314,10 @@ function modify_nested_subsystem(
         return fn(root)
     end
     # ignore the name of the root
-    if nameof(root) == hierarchy[1]
-        hierarchy = @view hierarchy[2:end]
+    if nameof(root) != hierarchy[1]
+        error("The name of the root system $(nameof(root)) must be included in the name passed to `modify_nested_subsystem`")
     end
+    hierarchy = @view hierarchy[2:end]
 
     # recursive helper function which does the searching and modification
     function _helper(sys::AbstractSystem, i::Int)
@@ -722,13 +723,14 @@ end
 
 A utility function to get the "canonical" form of a list of analysis points. Always returns
 a list of values. Any value that cannot be turned into an `AnalysisPoint` (i.e. isn't
-already an `AnalysisPoint` or `Symbol`) is simply wrapped in an array.
+already an `AnalysisPoint` or `Symbol`) is simply wrapped in an array. `Symbol` names of
+`AnalysisPoint`s are namespaced with `sys`.
 """
-canonicalize_ap(ap::Symbol) = [AnalysisPoint(ap)]
-canonicalize_ap(ap::AnalysisPoint) = [ap]
-canonicalize_ap(ap) = [ap]
-function canonicalize_ap(aps::Vector)
-    mapreduce(canonicalize_ap, vcat, aps; init = [])
+canonicalize_ap(sys::AbstractSystem, ap::Symbol) = [AnalysisPoint(renamespace(sys, ap))]
+canonicalize_ap(sys::AbstractSystem, ap::AnalysisPoint) = [ap]
+canonicalize_ap(sys::AbstractSystem, ap) = [ap]
+function canonicalize_ap(sys::AbstractSystem, aps::Vector)
+    mapreduce(Base.Fix1(canonicalize_ap, sys), vcat, aps; init = [])
 end
 
 """
@@ -737,7 +739,7 @@ end
 Given a list of analysis points, break the connection for each and set the output to zero.
 """
 function handle_loop_openings(sys::AbstractSystem, aps)
-    for ap in canonicalize_ap(aps)
+    for ap in canonicalize_ap(sys, aps)
         sys, (outvar,) = apply_transformation(Break(ap, true), sys)
         if Symbolics.isarraysymbolic(outvar)
             push!(get_eqs(sys), outvar ~ zeros(size(outvar)))
@@ -776,7 +778,7 @@ All other keyword arguments are forwarded to `linearization_function`.
 function get_linear_analysis_function(
         sys::AbstractSystem, transform, aps; system_modifier = identity, loop_openings = [], kwargs...)
     sys = handle_loop_openings(sys, loop_openings)
-    aps = canonicalize_ap(aps)
+    aps = canonicalize_ap(sys, aps)
     dus = []
     us = []
     for ap in aps
@@ -860,9 +862,7 @@ result of `apply_transformation`.
   with any required further modifications peformed.
 """
 function open_loop(sys, ap::Union{Symbol, AnalysisPoint}; system_modifier = identity)
-    if ap isa Symbol
-        ap = AnalysisPoint(ap)
-    end
+    ap = only(canonicalize_ap(sys, ap))
     tf = LoopTransferTransform(ap)
     sys, vars = apply_transformation(tf, sys)
     return system_modifier(sys), vars
@@ -871,9 +871,9 @@ end
 function linearization_function(sys::AbstractSystem,
         inputs::Union{Symbol, Vector{Symbol}, AnalysisPoint, Vector{AnalysisPoint}},
         outputs; loop_openings = [], system_modifier = identity, kwargs...)
-    loop_openings = Set(map(nameof, canonicalize_ap(loop_openings)))
-    inputs = canonicalize_ap(inputs)
-    outputs = canonicalize_ap(outputs)
+    loop_openings = Set(map(nameof, canonicalize_ap(sys, loop_openings)))
+    inputs = canonicalize_ap(sys, inputs)
+    outputs = canonicalize_ap(sys, outputs)
 
     input_vars = []
     for input in inputs
@@ -897,7 +897,7 @@ function linearization_function(sys::AbstractSystem,
         push!(output_vars, output_var)
     end
 
-    sys = handle_loop_openings(sys, collect(loop_openings))
+    sys = handle_loop_openings(sys, map(AnalysisPoint, collect(loop_openings)))
 
     return linearization_function(system_modifier(sys), input_vars, output_vars; kwargs...)
 end
diff --git a/test/analysis_points.jl b/test/analysis_points.jl
index 8c93e919b9..9e8e3fb455 100644
--- a/test/analysis_points.jl
+++ b/test/analysis_points.jl
@@ -2,6 +2,7 @@ using ModelingToolkit, ModelingToolkitStandardLibrary.Blocks
 using OrdinaryDiffEq, LinearAlgebra
 using Test
 using ModelingToolkit: t_nounits as t, D_nounits as D, AnalysisPoint, AbstractSystem
+import ModelingToolkit as MTK
 using Symbolics: NAMESPACE_SEPARATOR
 
 @testset "AnalysisPoint is lowered to `connect`" begin
@@ -69,26 +70,21 @@ sys = ODESystem(eqs, t, systems = [P, C], name = :hej)
     @test norm(sol.u[end]) < 1e-6 # This fails without the feedback through C
 end
 
-@testset "get_sensitivity - $name" for (name, sys, ap) in [
+test_cases = [
     ("inner", sys, sys.plant_input),
     ("nested", nested_sys, nested_sys.hej.plant_input),
-    ("inner - nonamespace", sys, :plant_input),
-    ("inner - Symbol", sys, nameof(sys.plant_input)),
-    ("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
+    ("inner - Symbol", sys, :plant_input),
+    ("nested - Symbol", nested_sys, nameof(sys.plant_input))
 ]
+
+@testset "get_sensitivity - $name" for (name, sys, ap) in test_cases
     matrices, _ = get_sensitivity(sys, ap)
     @test matrices.A[] == -2
     @test matrices.B[] * matrices.C[] == -1 # either one negative
     @test matrices.D[] == 1
 end
 
-@testset "get_comp_sensitivity - $name" for (name, sys, ap) in [
-    ("inner", sys, sys.plant_input),
-    ("nested", nested_sys, nested_sys.hej.plant_input),
-    ("inner - nonamespace", sys, :plant_input),
-    ("inner - Symbol", sys, nameof(sys.plant_input)),
-    ("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
-]
+@testset "get_comp_sensitivity - $name" for (name, sys, ap) in test_cases
     matrices, _ = get_comp_sensitivity(sys, ap)
     @test matrices.A[] == -2
     @test matrices.B[] * matrices.C[] == 1 # both positive or negative
@@ -104,13 +100,7 @@ S = sensitivity(P, C)      # or feedback(1, P*C)
 T = comp_sensitivity(P, C) # or feedback(P*C)
 =#
 
-@testset "get_looptransfer - $name" for (name, sys, ap) in [
-    ("inner", sys, sys.plant_input),
-    ("nested", nested_sys, nested_sys.hej.plant_input),
-    ("inner - nonamespace", sys, :plant_input),
-    ("inner - Symbol", sys, nameof(sys.plant_input)),
-    ("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
-]
+@testset "get_looptransfer - $name" for (name, sys, ap) in test_cases
     matrices, _ = get_looptransfer(sys, ap)
     @test matrices.A[] == -1
     @test matrices.B[] * matrices.C[] == -1 # either one negative
@@ -125,13 +115,7 @@ C = -1
 L = P*C
 =#
 
-@testset "open_loop - $name" for (name, sys, ap) in [
-    ("inner", sys, sys.plant_input),
-    ("nested", nested_sys, nested_sys.hej.plant_input),
-    ("inner - nonamespace", sys, :plant_input),
-    ("inner - Symbol", sys, nameof(sys.plant_input)),
-    ("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
-]
+@testset "open_loop - $name" for (name, sys, ap) in test_cases
     open_sys, (du, u) = open_loop(sys, ap)
     matrices, _ = linearize(open_sys, [du], [u])
     @test matrices.A[] == -1
@@ -146,13 +130,14 @@ eqs = [connect(P.output, :plant_output, C.input)
 sys = ODESystem(eqs, t, systems = [P, C], name = :hej)
 @named nested_sys = ODESystem(Equation[], t; systems = [sys])
 
-@testset "get_sensitivity - $name" for (name, sys, ap) in [
+test_cases = [
     ("inner", sys, sys.plant_input),
     ("nested", nested_sys, nested_sys.hej.plant_input),
-    ("inner - nonamespace", sys, :plant_input),
-    ("inner - Symbol", sys, nameof(sys.plant_input)),
-    ("nested - Symbol", nested_sys, nameof(nested_sys.hej.plant_input))
+    ("inner - Symbol", sys, :plant_input),
+    ("nested - Symbol", nested_sys, nameof(sys.plant_input))
 ]
+
+@testset "get_sensitivity - $name" for (name, sys, ap) in test_cases
     matrices, _ = get_sensitivity(sys, ap)
     @test matrices.A[] == -2
     @test matrices.B[] * matrices.C[] == -1 # either one negative
@@ -163,15 +148,19 @@ end
     ("inner", sys, sys.plant_input, sys.plant_output),
     ("nested", nested_sys, nested_sys.hej.plant_input, nested_sys.hej.plant_output)
 ]
+    inputname = Symbol(join(
+        MTK.namespace_hierarchy(nameof(inputap))[2:end], NAMESPACE_SEPARATOR))
+    outputname = Symbol(join(
+        MTK.namespace_hierarchy(nameof(outputap))[2:end], NAMESPACE_SEPARATOR))
     @testset "input - $(typeof(input)), output - $(typeof(output))" for (input, output) in [
         (inputap, outputap),
-        (nameof(inputap), outputap),
-        (inputap, nameof(outputap)),
-        (nameof(inputap), nameof(outputap)),
+        (inputname, outputap),
+        (inputap, outputname),
+        (inputname, outputname),
         (inputap, [outputap]),
-        (nameof(inputap), [outputap]),
-        (inputap, [nameof(outputap)]),
-        (nameof(inputap), [nameof(outputap)])
+        (inputname, [outputap]),
+        (inputap, [outputname]),
+        (inputname, [outputname])
     ]
         matrices, _ = linearize(sys, input, output)
         # Result should be the same as feedpack(P, 1), i.e., the closed-loop transfer function from plant input to plant output
diff --git a/test/downstream/analysis_points.jl b/test/downstream/analysis_points.jl
index fabca0f493..7a433cb504 100644
--- a/test/downstream/analysis_points.jl
+++ b/test/downstream/analysis_points.jl
@@ -64,7 +64,7 @@ import ControlSystemsBase as CS
     prob = ODEProblem(sys, unknowns(sys) .=> 0.0, (0.0, 4.0))
     sol = solve(prob, Rodas5P(), reltol = 1e-6, abstol = 1e-9)
 
-    matrices, ssys = linearize(closed_loop, AnalysisPoint(:r), AnalysisPoint(:y))
+    matrices, ssys = linearize(closed_loop, :r, :y)
     lsys = ss(matrices...) |> sminreal
     @test lsys.nx == 8
 
@@ -129,13 +129,11 @@ end
         t,
         systems = [P_inner, feedback, ref])
 
-    P_not_broken, _ = linearize(sys_inner, AnalysisPoint(:u), AnalysisPoint(:y))
+    P_not_broken, _ = linearize(sys_inner, :u, :y)
     @test P_not_broken.A[] == -2
-    P_broken, _ = linearize(sys_inner, AnalysisPoint(:u), AnalysisPoint(:y),
-        loop_openings = [AnalysisPoint(:u)])
+    P_broken, _ = linearize(sys_inner, :u, :y, loop_openings = [:u])
     @test P_broken.A[] == -1
-    P_broken, _ = linearize(sys_inner, AnalysisPoint(:u), AnalysisPoint(:y),
-        loop_openings = [AnalysisPoint(:y)])
+    P_broken, _ = linearize(sys_inner, :u, :y, loop_openings = [:y])
     @test P_broken.A[] == -1
 
     Sinner = sminreal(ss(get_sensitivity(sys_inner, :u)[1]...))
@@ -154,10 +152,10 @@ end
         t,
         systems = [P_outer, sys_inner])
 
-    Souter = sminreal(ss(get_sensitivity(sys_outer, sys_inner.u)[1]...))
+    Souter = sminreal(ss(get_sensitivity(sys_outer, sys_outer.sys_inner.u)[1]...))
 
     Sinner2 = sminreal(ss(get_sensitivity(
-        sys_outer, sys_inner.u, loop_openings = [:y2])[1]...))
+        sys_outer, sys_outer.sys_inner.u, loop_openings = [:y2])[1]...))
 
     @test Sinner.nx == 1
     @test Sinner == Sinner2
@@ -268,7 +266,7 @@ end
     @test tf(L[2, 1]) ≈ tf(Ps)
 
     matrices, _ = linearize(
-        sys_outer, [sys_outer.inner.plant_input], [sys_inner.plant_output])
+        sys_outer, [sys_outer.inner.plant_input], [nameof(sys_inner.plant_output)])
     G = CS.ss(matrices...) |> sminreal
     @test tf(G) ≈ tf(CS.feedback(Ps, Cs))
 end

From c93df8fb418af2db96ca215028657ae8f9f73e5f Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 6 Jan 2025 15:13:50 +0530
Subject: [PATCH 20/23] refactor: do not export individual analysis point
 transformations

---
 src/ModelingToolkit.jl | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 29288a5742..7462122998 100644
--- a/src/ModelingToolkit.jl
+++ b/src/ModelingToolkit.jl
@@ -293,8 +293,7 @@ export MTKParameters, reorder_dimension_by_tunables!, reorder_dimension_by_tunab
 
 export HomotopyContinuationProblem
 
-export AnalysisPoint, Break, PerturbOutput, SampleInput, SensitivityTransform,
-       ComplementarySensitivityTransform, LoopTransferTransform, apply_transformation,
-       get_sensitivity_function, get_comp_sensitivity_function, get_looptransfer_function,
-       get_sensitivity, get_comp_sensitivity, get_looptransfer, open_loop
+export AnalysisPoint, get_sensitivity_function, get_comp_sensitivity_function,
+       get_looptransfer_function, get_sensitivity, get_comp_sensitivity, get_looptransfer,
+       open_loop
 end # module

From a7650737cbcfdd99def579c193b59c4dedec04d1 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 6 Jan 2025 15:14:22 +0530
Subject: [PATCH 21/23] docs: port over linear analysis tutorial from MTKStdlib

---
 docs/Project.toml                     |   3 +
 docs/pages.jl                         |   3 +-
 docs/src/tutorials/linear_analysis.md | 152 ++++++++++++++++++++++++++
 3 files changed, 157 insertions(+), 1 deletion(-)
 create mode 100644 docs/src/tutorials/linear_analysis.md

diff --git a/docs/Project.toml b/docs/Project.toml
index a358455503..be80250603 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
+ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -8,6 +9,7 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
@@ -31,6 +33,7 @@ Distributions = "0.25"
 Documenter = "1"
 DynamicQuantities = "^0.11.2, 0.12, 1"
 ModelingToolkit = "8.33, 9"
+ModelingToolkitStandardLibrary = "2.19"
 NonlinearSolve = "3, 4"
 Optim = "1.7"
 Optimization = "3.9, 4"
diff --git a/docs/pages.jl b/docs/pages.jl
index 2af487adf8..8a95de4ac5 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -13,7 +13,8 @@ pages = [
         "tutorials/bifurcation_diagram_computation.md",
         "tutorials/SampledData.md",
         "tutorials/domain_connections.md",
-        "tutorials/callable_params.md"],
+        "tutorials/callable_params.md",
+        "tutorials/linear_analysis.md"],
     "Examples" => Any[
         "Basic Examples" => Any["examples/higher_order.md",
             "examples/spring_mass.md",
diff --git a/docs/src/tutorials/linear_analysis.md b/docs/src/tutorials/linear_analysis.md
new file mode 100644
index 0000000000..7fcb67f9fb
--- /dev/null
+++ b/docs/src/tutorials/linear_analysis.md
@@ -0,0 +1,152 @@
+# 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`](@ref), 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`](@ref) get the [sensitivity function (wiki)](https://en.wikipedia.org/wiki/Sensitivity_(control_systems)), $S(s)$, as defined in the field of control theory.
+  - [`get_comp_sensitivity`](@ref) get the complementary sensitivity function $T(s) : S(s)+T(s)=1$.
+  - [`get_looptransfer`](@ref) 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`](@ref) can be called with two analysis points denoting the input and output of the linearized system.
+  - [`open_loop`](@ref) 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`](@ref), or automatically when connecting two causal components using `connect`:
+
+```julia
+connect(comp1.output, :analysis_point_name, comp2.input)
+```
+
+A single output can also be connected to multiple inputs:
+
+```julia
+connect(comp1.output, :analysis_point_name, comp2.input, comp3.input, comp4.input)
+```
+
+!!! warning "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`](@ref) return the output of [`ModelingToolkit.linearize`](@ref), which is
+
+```julia
+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
+
+```math
+\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
+
+```@example LINEAR_ANALYSIS
+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]
+```
+
+Continued linear analysis and design can be performed using ControlSystemsBase.jl.
+We create `ControlSystemsBase.StateSpace` objects using
+
+```@example LINEAR_ANALYSIS
+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
+
+```@example LINEAR_ANALYSIS_CS
+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)
+```
+
+We may also derive the loop-transfer function $L(s) = P(s)C(s)$ using
+
+```@example LINEAR_ANALYSIS
+matrices_L = get_looptransfer(sys, :plant_output)[1]
+L = ss(matrices_L...)
+```
+
+which is equivalent to the following with ControlSystems
+
+```@example LINEAR_ANALYSIS_CS
+L = P * (-C) # Add the minus sign to build the negative feedback into the controller
+```
+
+To obtain the transfer function between two analysis points, we call `linearize`
+
+```@example LINEAR_ANALYSIS
+using ModelingToolkit # hide
+matrices_PS = linearize(sys, :plant_input, :plant_output)[1]
+```
+
+this particular transfer function should be equivalent to the linear system `P(s)S(s)`, i.e., equivalent to
+
+```@example LINEAR_ANALYSIS_CS
+feedback(P, C)
+```
+
+### Obtaining transfer functions
+
+A statespace system from [ControlSystemsBase](https://juliacontrol.github.io/ControlSystems.jl/stable/man/creating_systems/) can be converted to a transfer function using the function `tf`:
+
+```@example LINEAR_ANALYSIS_CS
+tf(S)
+```
+
+## Gain and phase margins
+
+Further linear analysis can be performed using the [analysis methods from ControlSystemsBase](https://juliacontrol.github.io/ControlSystems.jl/stable/lib/analysis/). For example, calculating the gain and phase margins of a system can be done using
+
+```@example LINEAR_ANALYSIS_CS
+margin(P)
+```
+
+(they are infinite for this system). A Nyquist plot can be produced using
+
+```@example LINEAR_ANALYSIS_CS
+nyquistplot(P)
+```
+
+## Index
+
+```@index
+Pages = ["linear_analysis.md"]
+```
+
+```@autodocs
+Modules = [ModelingToolkit]
+Pages   = ["systems/analysis_points.jl"]
+Order   = [:function, :type]
+Private = false
+```

From 520daf49d589b49c05c172d8a29e353170bf6c05 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Mon, 6 Jan 2025 18:08:13 +0530
Subject: [PATCH 22/23] docs: remove DifferentialEquations.jl from docs
 environment

---
 docs/Project.toml                                 | 2 --
 docs/src/basics/Composition.md                    | 2 +-
 docs/src/examples/perturbation.md                 | 2 +-
 docs/src/examples/sparse_jacobians.md             | 2 +-
 docs/src/examples/spring_mass.md                  | 2 +-
 docs/src/tutorials/acausal_components.md          | 2 +-
 docs/src/tutorials/modelingtoolkitize.md          | 4 ++--
 docs/src/tutorials/ode_modeling.md                | 6 +++---
 docs/src/tutorials/programmatically_generating.md | 2 +-
 9 files changed, 11 insertions(+), 13 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index be80250603..40fcb16581 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -3,7 +3,6 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
@@ -28,7 +27,6 @@ Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 BenchmarkTools = "1.3"
 BifurcationKit = "0.4"
 DataInterpolations = "6.5"
-DifferentialEquations = "7.6"
 Distributions = "0.25"
 Documenter = "1"
 DynamicQuantities = "^0.11.2, 0.12, 1"
diff --git a/docs/src/basics/Composition.md b/docs/src/basics/Composition.md
index 78dec8445b..2e5d4be831 100644
--- a/docs/src/basics/Composition.md
+++ b/docs/src/basics/Composition.md
@@ -56,7 +56,7 @@ x0 = [decay1.x => 1.0
 p = [decay1.a => 0.1
      decay2.a => 0.2]
 
-using DifferentialEquations
+using OrdinaryDiffEq
 prob = ODEProblem(simplified_sys, x0, (0.0, 100.0), p)
 sol = solve(prob, Tsit5())
 sol[decay2.f]
diff --git a/docs/src/examples/perturbation.md b/docs/src/examples/perturbation.md
index 407248709d..0d1a493cb4 100644
--- a/docs/src/examples/perturbation.md
+++ b/docs/src/examples/perturbation.md
@@ -49,7 +49,7 @@ These are the ODEs we want to solve. Now construct an `ODESystem`, which automat
 To solve the `ODESystem`, we generate an `ODEProblem` with initial conditions $x(0) = 0$, and $ẋ(0) = 1$, and solve it:
 
 ```@example perturbation
-using DifferentialEquations
+using OrdinaryDiffEq
 u0 = Dict([unknowns(sys) .=> 0.0; D(y[0]) => 1.0]) # nonzero initial velocity
 prob = ODEProblem(sys, u0, (0.0, 3.0))
 sol = solve(prob)
diff --git a/docs/src/examples/sparse_jacobians.md b/docs/src/examples/sparse_jacobians.md
index 1ed3b1733c..03bd80d432 100644
--- a/docs/src/examples/sparse_jacobians.md
+++ b/docs/src/examples/sparse_jacobians.md
@@ -11,7 +11,7 @@ First, let's start out with an implementation of the 2-dimensional Brusselator
 partial differential equation discretized using finite differences:
 
 ```@example sparsejac
-using DifferentialEquations, ModelingToolkit
+using OrdinaryDiffEq, ModelingToolkit
 
 const N = 32
 const xyd_brusselator = range(0, stop = 1, length = N)
diff --git a/docs/src/examples/spring_mass.md b/docs/src/examples/spring_mass.md
index e733b11724..355e5c20b2 100644
--- a/docs/src/examples/spring_mass.md
+++ b/docs/src/examples/spring_mass.md
@@ -5,7 +5,7 @@ In this tutorial, we will build a simple component-based model of a spring-mass
 ## Copy-Paste Example
 
 ```@example component
-using ModelingToolkit, Plots, DifferentialEquations, LinearAlgebra
+using ModelingToolkit, Plots, OrdinaryDiffEq, LinearAlgebra
 using ModelingToolkit: t_nounits as t, D_nounits as D
 using Symbolics: scalarize
 
diff --git a/docs/src/tutorials/acausal_components.md b/docs/src/tutorials/acausal_components.md
index c91ba29670..096b576d29 100644
--- a/docs/src/tutorials/acausal_components.md
+++ b/docs/src/tutorials/acausal_components.md
@@ -20,7 +20,7 @@ equalities before solving. Let's see this in action.
 ## Copy-Paste Example
 
 ```@example acausal
-using ModelingToolkit, Plots, DifferentialEquations
+using ModelingToolkit, Plots, OrdinaryDiffEq
 using ModelingToolkit: t_nounits as t, D_nounits as D
 
 @connector Pin begin
diff --git a/docs/src/tutorials/modelingtoolkitize.md b/docs/src/tutorials/modelingtoolkitize.md
index e272a2237c..545879f842 100644
--- a/docs/src/tutorials/modelingtoolkitize.md
+++ b/docs/src/tutorials/modelingtoolkitize.md
@@ -33,10 +33,10 @@ to improve a simulation code before it's passed to the solver.
 ## Example Usage: Generating an Analytical Jacobian Expression for an ODE Code
 
 Take, for example, the Robertson ODE
-defined as an `ODEProblem` for DifferentialEquations.jl:
+defined as an `ODEProblem` for OrdinaryDiffEq.jl:
 
 ```@example mtkize
-using DifferentialEquations, ModelingToolkit
+using OrdinaryDiffEq, ModelingToolkit
 function rober(du, u, p, t)
     y₁, y₂, y₃ = u
     k₁, k₂, k₃ = p
diff --git a/docs/src/tutorials/ode_modeling.md b/docs/src/tutorials/ode_modeling.md
index 755e70ef46..0a4cd80803 100644
--- a/docs/src/tutorials/ode_modeling.md
+++ b/docs/src/tutorials/ode_modeling.md
@@ -34,7 +34,7 @@ using ModelingToolkit: t_nounits as t, D_nounits as D
     end
 end
 
-using DifferentialEquations: solve
+using OrdinaryDiffEq
 @mtkbuild fol = FOL()
 prob = ODEProblem(fol, [], (0.0, 10.0), [])
 sol = solve(prob)
@@ -84,10 +84,10 @@ Note that equations in MTK use the tilde character (`~`) as equality sign.
 
 `@mtkbuild` creates an instance of `FOL` named as `fol`.
 
-After construction of the ODE, you can solve it using [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/):
+After construction of the ODE, you can solve it using [OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs/stable/):
 
 ```@example ode2
-using DifferentialEquations
+using OrdinaryDiffEq
 using Plots
 
 prob = ODEProblem(fol, [], (0.0, 10.0), [])
diff --git a/docs/src/tutorials/programmatically_generating.md b/docs/src/tutorials/programmatically_generating.md
index 76d12dba2a..9fc1db1834 100644
--- a/docs/src/tutorials/programmatically_generating.md
+++ b/docs/src/tutorials/programmatically_generating.md
@@ -49,7 +49,7 @@ eqs = [D(x) ~ (h - x) / τ] # create an array of equations
 # Note: Complete models cannot be subsystems of other models!
 fol = structural_simplify(model)
 prob = ODEProblem(fol, [], (0.0, 10.0), [])
-using DifferentialEquations: solve
+using OrdinaryDiffEq
 sol = solve(prob)
 
 using Plots

From 69e05f5b5031acc3902acf9edd5430dad75efa73 Mon Sep 17 00:00:00 2001
From: Aayush Sabharwal <aayush.sabharwal@gmail.com>
Date: Tue, 7 Jan 2025 12:35:57 +0530
Subject: [PATCH 23/23] docs: ignore Scholarpedia link in Documenter linkcheck

---
 docs/make.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/docs/make.jl b/docs/make.jl
index f380867e6c..36a27bf598 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -25,7 +25,12 @@ makedocs(sitename = "ModelingToolkit.jl",
     modules = [ModelingToolkit],
     clean = true, doctest = false, linkcheck = true,
     warnonly = [:docs_block, :missing_docs, :cross_references],
-    linkcheck_ignore = ["https://epubs.siam.org/doi/10.1137/0903023"],
+    linkcheck_ignore = [
+        "https://epubs.siam.org/doi/10.1137/0903023",
+        # this link tends to fail linkcheck stochastically and often takes much longer to succeed
+        # even in the browser it takes ages
+        "http://www.scholarpedia.org/article/Differential-algebraic_equations"
+    ],
     format = Documenter.HTML(;
         assets = ["assets/favicon.ico"],
         mathengine,