diff --git a/Project.toml b/Project.toml
index 7b0997823c..256b1de6fc 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"
@@ -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/docs/Project.toml b/docs/Project.toml
index a358455503..40fcb16581 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,13 +1,14 @@
 [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"
 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"
@@ -26,11 +27,11 @@ 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"
 ModelingToolkit = "8.33, 9"
+ModelingToolkitStandardLibrary = "2.19"
 NonlinearSolve = "3, 4"
 Optim = "1.7"
 Optimization = "3.9, 4"
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,
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/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/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
+```
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
diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl
index 10aba4d8a9..7462122998 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,
@@ -291,4 +293,7 @@ export MTKParameters, reorder_dimension_by_tunables!, reorder_dimension_by_tunab
 
 export HomotopyContinuationProblem
 
+export AnalysisPoint, get_sensitivity_function, get_comp_sensitivity_function,
+       get_looptransfer_function, get_sensitivity, get_comp_sensitivity, get_looptransfer,
+       open_loop
 end # module
diff --git a/src/systems/abstractsystem.jl b/src/systems/abstractsystem.jl
index 168260ae69..3790cb8eb0 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
 
@@ -2366,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...)
diff --git a/src/systems/analysis_points.jl b/src/systems/analysis_points.jl
new file mode 100644
index 0000000000..87f278d4dc
--- /dev/null
+++ b/src/systems/analysis_points.jl
@@ -0,0 +1,946 @@
+"""
+    $(TYPEDEF)
+    AnalysisPoint(input, name::Symbol, outputs::Vector)
+
+Create an AnalysisPoint for linear analysis. Analysis points can be created by calling
+
+```
+connect(out, :ap_name, in...)
+```
+
+Where `out` is the output being connected to the inputs `in...`. All involved
+connectors (input and outputs) are required to either have an unknown named
+`u` or a single unknown, all of which should have the same size.
+
+See also [`get_sensitivity`](@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}}
+
+    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)
+"""
+    $(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
+
+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...; verbose = true)
+    return AnalysisPoint() ~ AnalysisPoint(in, ap.name, collect(outs); verbose)
+end
+
+"""
+    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
+
+```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.
+
+# Keyword arguments
+
+- `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...; verbose = true)
+    return AnalysisPoint() ~ AnalysisPoint(in, name, [out; collect(outs)]; verbose)
+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]
+        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)
+        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])
+            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 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(eltype(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
+
+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)
+
+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))`.
+
+$DOC_WILL_REMOVE_AP
+
+$DOC_ADDED_VARIABLE
+
+## 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. Does not modify the system.
+
+## Fields
+
+$(TYPEDFIELDS)
+"""
+struct GetInput <: AnalysisPointTransformation
+    """
+    The analysis point to get the input of.
+    """
+    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 and will be the second variable in the tuple of new variables returned
+from `apply_transformation`.
+
+$DOC_WILL_REMOVE_AP
+
+$DOC_ADDED_VARIABLE
+
+## 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 + wrap(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 + wrap(new_var))
+        push!(unks, out_var)
+
+        return ap_sys, (new_var, out_var)
+    end
+end
+
+"""
+    $(TYPEDEF)
+
+A transformation which adds a variable named `name` to the system containing the analysis
+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)
+    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`.
+The returned added variables are `(du, u)` where `du` is the perturbation added to the
+input, and `u` is the output after perturbation.
+
+$DOC_WILL_REMOVE_AP
+
+$DOC_ADDED_VARIABLE
+"""
+SensitivityTransform(ap::AnalysisPoint) = PerturbOutput(ap, true)
+
+"""
+    $(TYPEDEF)
+
+A transformation to enable calculating the complementary sensitivity function about the
+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
+
+# Fields
+
+$(TYPEDFIELDS)
+"""
+struct ComplementarySensitivityTransform <: AnalysisPointTransformation
+    """
+    The analysis point to modify.
+    """
+    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(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`
+    # but comp sensitivity wants `output + du ~ input`. Thus, `du ~ -_du`.
+    eqs = copy(get_eqs(sys))
+    @set! sys.eqs = eqs
+    push!(eqs, du ~ -wrap(_du))
+
+    defs = copy(get_defaults(sys))
+    @set! sys.defaults = defs
+    defs[du] = -wrap(_du)
+    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
+
+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
+
+"""
+    $(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. `Symbol` names of
+`AnalysisPoint`s are namespaced with `sys`.
+"""
+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
+
+"""
+    $(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(sys, 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
+
+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(sys, aps)
+    dus = []
+    us = []
+    for ap in aps
+        sys, (du, u) = apply_transformation(transform(ap), sys)
+        push!(dus, du)
+        push!(us, u)
+    end
+    linearization_function(system_modifier(sys), dus, us; kwargs...)
+end
+
+"""
+    $(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
+
+"""
+    $(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 = $(utility_fun)(
+            sys, ap, args...; loop_openings, system_modifier, kwargs...)
+        ModelingToolkit.linearize(ssys, lin_fun; kwargs...), ssys
+    end
+end
+
+"""
+    $(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)
+    ap = only(canonicalize_ap(sys, ap))
+    tf = LoopTransferTransform(ap)
+    sys, vars = apply_transformation(tf, sys)
+    return system_modifier(sys), vars
+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(sys, loop_openings)))
+    inputs = canonicalize_ap(sys, inputs)
+    outputs = canonicalize_ap(sys, outputs)
+
+    input_vars = []
+    for input in inputs
+        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(AddVariable(output), sys)
+            sys, (input_var,) = apply_transformation(GetInput(output), sys)
+            push!(get_eqs(sys), output_var ~ input_var)
+        else
+            output_var = output
+        end
+        push!(output_vars, output_var)
+    end
+
+    sys = handle_loop_openings(sys, map(AnalysisPoint, collect(loop_openings)))
+
+    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
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)
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
diff --git a/test/analysis_points.jl b/test/analysis_points.jl
new file mode 100644
index 0000000000..9e8e3fb455
--- /dev/null
+++ b/test/analysis_points.jl
@@ -0,0 +1,181 @@
+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
+    @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_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_normal)
+
+    @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)
+    @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
+
+test_cases = [
+    ("inner", sys, sys.plant_input),
+    ("nested", nested_sys, 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 test_cases
+    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 test_cases
+    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 test_cases
+    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])
+
+test_cases = [
+    ("inner", sys, sys.plant_input),
+    ("nested", nested_sys, 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 "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)
+]
+    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),
+        (inputname, outputap),
+        (inputap, outputname),
+        (inputname, outputname),
+        (inputap, [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
+        @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
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
new file mode 100644
index 0000000000..7a433cb504
--- /dev/null
+++ b/test/downstream/analysis_points.jl
@@ -0,0 +1,272 @@
+using ModelingToolkit, OrdinaryDiffEq, LinearAlgebra, ControlSystemsBase
+using ModelingToolkitStandardLibrary.Mechanical.Rotational
+using ModelingToolkitStandardLibrary.Blocks
+using ModelingToolkit: connect, 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, :r, :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, :u, :y)
+    @test P_not_broken.A[] == -2
+    P_broken, _ = linearize(sys_inner, :u, :y, loop_openings = [:u])
+    @test P_broken.A[] == -1
+    P_broken, _ = linearize(sys_inner, :u, :y, loop_openings = [: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_outer.sys_inner.u)[1]...))
+
+    Sinner2 = sminreal(ss(get_sensitivity(
+        sys_outer, 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, _ = 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, _ = 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, _ = get_looptransfer(
+        sys, :plant_input)
+    L = Kss * Pss
+    @test CS.tf(CS.ss(matrices...)) ≈ CS.tf(L)
+
+    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
+
+@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, [sys_outer.inner.plant_input, 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)
+
+    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], [nameof(sys_inner.plant_output)])
+    G = CS.ss(matrices...) |> sminreal
+    @test tf(G) ≈ tf(CS.feedback(Ps, Cs))
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 4a2ec80e6a..d38aa86bb2 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
 
@@ -110,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"