Merge pull request #3285 from AayushSabharwal/as/analysis-points

feat: add analysis points

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"

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"

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,

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",

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]

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)

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)

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
 

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

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
+```

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

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), [])

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

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

src/systems/abstractsystem.jl

@@ -1163,6 +1163,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
 
@@ -2372,6 +2380,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...)