Add open_loop

depends on SciML/ModelingToolkit.jl#1824

Author: baggepinnen

src/Blocks/Blocks.jl

@@ -27,7 +27,7 @@ export PI, LimPI, PID, LimPID
 include("continuous.jl")
 
 export AnalysisPoint, expand_analysis_points, get_sensitivity, get_comp_sensitivity,
-       get_looptransfer
+       get_looptransfer, open_loop
 include("analysis_points.jl")
 
 end

src/Blocks/analysis_points.jl

@@ -188,8 +188,6 @@ See also [`get_sensitivity`](@ref), [`get_comp_sensitivity`](@ref).
 """
 function get_looptransfer(sys, ap::AnalysisPoint; kwargs...)
     t = get_iv(sys)
-    @variables u(t) = 0
-    @variables y(t) = 0
     new_eqs = map(get_eqs(sys)) do eq
         eq.rhs == ap || (return eq)
         0 ~ 0 # we just want to open the connection
@@ -198,8 +196,36 @@ function get_looptransfer(sys, ap::AnalysisPoint; kwargs...)
     ModelingToolkit.linearize(sys, [ap.out.u], [ap.in.u]; kwargs...)
 end
 
+"""
+    open_sys = open_loop(sys, ap::AnalysisPoint; kwargs)
+    open_sys = open_loop(sys, ap_name::Symbol; kwargs)
+
+Open the loop at analysis point `ap` by breaking the connection through `ap`.
+
+`open_sys` will have `ap.out` as input and `ap.in` as output.
+
+# Arguments:
+- `kwargs`: Are sent to `ModelingToolkit.linearize`
+
+See also [`get_sensitivity`](@ref), [`get_comp_sensitivity`](@ref), [`get_looptransfer`](@ref).
+"""
+function open_loop(sys, ap::AnalysisPoint; kwargs...)
+    t = get_iv(sys)
+    @variables u(t)=0 [input = true]
+    @variables y(t)=0 [output = true]
+    new_eqs = map(get_eqs(sys)) do eq
+        eq.rhs == ap || (return [eq])
+        [ap.out.u ~ u
+         ap.in.u ~ y]
+    end
+    new_eqs = reduce(vcat, new_eqs)
+    @set! sys.eqs = new_eqs
+    @set! sys.states = [states(sys); u; y]
+    sys
+end
+
 # Add a method to get_sensitivity that accepts the name of an AnalysisPoint 
-for f in [:get_sensitivity, :get_comp_sensitivity, :get_looptransfer]
+for f in [:get_sensitivity, :get_comp_sensitivity, :get_looptransfer, :open_loop]
     @eval function $f(sys, ap_name::Symbol, args...; kwargs...)
         $f(sys, find_analysis_point(sys, ap_name), args...; kwargs...)
     end

test/Blocks/test_analysis_points.jl

@@ -67,3 +67,13 @@ P = tf(1.0, [1, 1])
 C = -1
 L = P*C
 =#
+
+# Open loop
+open_sys = Blocks.open_loop(sys, :plant_input)
+@unpack u, y = open_sys
+
+# Linearizing the open-loop system should yield the same system as get_looptransfer
+matrices, _ = linearize(open_sys, [u], [y])
+@test matrices.A[] == -1
+@test matrices.B[] * matrices.C[] == -1 # either one negative
+@test matrices.D[] == 0