|
| 1 | +using ModelingToolkit: get_eqs, vars, @set!, get_iv |
| 2 | + |
| 3 | +Base.@kwdef mutable struct AnalysisPoint |
| 4 | + in = nothing |
| 5 | + out = nothing |
| 6 | + name::Symbol |
| 7 | +end |
| 8 | + |
| 9 | +""" |
| 10 | + AnalysisPoint(in, out, name::Symbol) |
| 11 | + AnalysisPoint(in, out; name::Symbol) |
| 12 | + AnalysisPoint(name::Symbol) |
| 13 | +
|
| 14 | +Create an AnalysisPoint for linear analysis. Analysis points can also e created automatically by calling |
| 15 | +``` |
| 16 | +connect(in, :ap_name, out) |
| 17 | +``` |
| 18 | +
|
| 19 | +# Arguments: |
| 20 | +- `in`: A connector of type [`RealOutput`](@ref). |
| 21 | +- `out`: A connector of type [`RealInput`](@ref). |
| 22 | +- `name`: The name of the analysis point. |
| 23 | +
|
| 24 | +See also [`get_sensitivity`](@ref), [`get_comp_sensitivity`](@ref) |
| 25 | +
|
| 26 | +# Example |
| 27 | +```julia |
| 28 | +using ModelingToolkitStandardLibrary.Blocks |
| 29 | +@named P = FirstOrder(k=1, T=1) |
| 30 | +@named C = Gain(-1) |
| 31 | +t = ModelingToolkit.get_iv(P) |
| 32 | +eqs = [ |
| 33 | + connect(P.output, C.input) |
| 34 | + connect(C.output, :plant_input, P.input) # Connect with an automatically created analysis point |
| 35 | +] |
| 36 | +sys = ODESystem(eqs, t, systems=[P,C], name=:feedback_system) |
| 37 | +
|
| 38 | +matrices_S, _ = get_sensitivity(sys, :plant_input) # Compute the matrices of a state-space representation of the (input) sensitivity funciton. |
| 39 | +matrices_T, _ = get_comp_sensitivity(sys, :plant_input) |
| 40 | +``` |
| 41 | +Continued linear analysis and design can be performed using ControlSystemsBase.jl. |
| 42 | +Create `ControlSystemsBase.StateSpace` objects using |
| 43 | +```julia |
| 44 | +using ControlSystemsBase, Plots |
| 45 | +S = ss(matrices_S...) |
| 46 | +T = ss(matrices_T...) |
| 47 | +bodeplot([S, T], lab=["S" "T"]) |
| 48 | +``` |
| 49 | +The sensitivity functions obtained this way should be equivalent to the ones obtained with the code below |
| 50 | +``` |
| 51 | +using ControlSystemsBase |
| 52 | +P = tf(1.0, [1, 1]) |
| 53 | +C = 1 # Negative feedback assumed in ControlSystems |
| 54 | +S = sensitivity(P, C) # or feedback(1, P*C) |
| 55 | +T = comp_sensitivity(P, C) # or feedback(P*C) |
| 56 | +``` |
| 57 | +""" |
| 58 | +AnalysisPoint(in, out; name) = AnalysisPoint(in, out, name) |
| 59 | +AnalysisPoint(name) = AnalysisPoint(; name) |
| 60 | + |
| 61 | +Base.show(io::IO, ap::AnalysisPoint) = show(io, MIME"text/plain"(), ap) |
| 62 | +function Base.show(io::IO, ::MIME"text/plain", ap::AnalysisPoint) |
| 63 | + if get(io, :compact, false) |
| 64 | + print(io, "AnalysisPoint($(ap.in.u), $(ap.out.u); name=$(ap.name))") |
| 65 | + else |
| 66 | + print(io, "AnalysisPoint(") |
| 67 | + printstyled(io, ap.name, color = :cyan) |
| 68 | + if ap.in !== nothing && ap.out !== nothing |
| 69 | + print(io, " from ") |
| 70 | + printstyled(io, ap.in.u, color = :green) |
| 71 | + print(io, " to ") |
| 72 | + printstyled(io, ap.out.u, color = :blue) |
| 73 | + end |
| 74 | + print(io, ")") |
| 75 | + end |
| 76 | +end |
| 77 | + |
| 78 | +""" |
| 79 | + connect(in, ap::AnalysisPoint, out) |
| 80 | + connect(in, ap_name::Symbol, out) |
| 81 | +
|
| 82 | +Connect `in` and `out` with an [`AnalysisPoint`](@ref) inbetween. |
| 83 | +The incoming connection `in` is expected to be of type [`RealOutput`](@ref), and vice versa. |
| 84 | +
|
| 85 | +# Arguments: |
| 86 | +- `in`: A connector of type [`RealOutput`](@ref) |
| 87 | +- `out`: A connector of type [`RealInput`](@ref) |
| 88 | +- `ap`: An explicitly created [`AnalysisPoint`](@ref) |
| 89 | +- `ap_name`: If a name is given, an [`AnalysisPoint`](@ref) with the given name will be created automatically. |
| 90 | +""" |
| 91 | +function ModelingToolkit.connect(in, ap::AnalysisPoint, out) |
| 92 | + ap.in = in |
| 93 | + ap.out = out |
| 94 | + return 0 ~ ap |
| 95 | +end |
| 96 | + |
| 97 | +function ModelingToolkit.connect(in, ap_name::Symbol, out) |
| 98 | + return 0 ~ AnalysisPoint(in, out, ap_name) |
| 99 | +end |
| 100 | + |
| 101 | +function ModelingToolkit.vars(ap::AnalysisPoint; op = Differential) |
| 102 | + vars(connect(ap.in, ap.out); op) |
| 103 | +end |
| 104 | + |
| 105 | +""" |
| 106 | + find_analysis_point(sys, name::Symbol) |
| 107 | +
|
| 108 | +Find and return the analysis point in `sys` with the specified `name`. If no matching [`AnalysisPoint`](@ref) is found, `nothing` is returned. |
| 109 | +""" |
| 110 | +function find_analysis_point(sys, name) |
| 111 | + for eq in equations(sys) |
| 112 | + eq.rhs isa AnalysisPoint && eq.rhs.name == name && (return eq.rhs) |
| 113 | + end |
| 114 | + nothing |
| 115 | +end |
| 116 | + |
| 117 | +""" |
| 118 | + expand_analysis_points(sys) |
| 119 | +
|
| 120 | +Replace analysis points with the identity connection connect(ap.in, ap.out). This is called before a system containing analysis points is simulated, in which case analysis points have no effect. |
| 121 | +""" |
| 122 | +function expand_analysis_points(sys) |
| 123 | + new_eqs = map(get_eqs(sys)) do eq |
| 124 | + eq.rhs isa AnalysisPoint || (return eq) |
| 125 | + ap = eq.rhs |
| 126 | + connect(ap.in, ap.out) |
| 127 | + end |
| 128 | + @set! sys.eqs = new_eqs |
| 129 | + sys |
| 130 | +end |
| 131 | + |
| 132 | +""" |
| 133 | + get_sensitivity(sys, ap::AnalysisPoint; kwargs) |
| 134 | + get_sensitivity(sys, ap_name::Symbol; kwargs) |
| 135 | +
|
| 136 | +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`. |
| 137 | +
|
| 138 | +# Arguments: |
| 139 | +- `kwargs`: Are sent to `ModelingToolkit.linearize` |
| 140 | +
|
| 141 | +See also [`get_comp_sensitivity`](@ref). |
| 142 | +""" |
| 143 | +function get_sensitivity(sys, ap::AnalysisPoint; kwargs...) |
| 144 | + t = get_iv(sys) |
| 145 | + @variables d(t) = 0 # Perturbantion serving as input to sensivity transfer function |
| 146 | + new_eqs = map(get_eqs(sys)) do eq |
| 147 | + eq.rhs == ap || (return eq) |
| 148 | + ap.out.u ~ ap.in.u + d # This assumes that the connector as an internal vaiable named u |
| 149 | + end |
| 150 | + @set! sys.eqs = new_eqs |
| 151 | + @set! sys.states = [states(sys); d] |
| 152 | + ModelingToolkit.linearize(sys, [d], [ap.out.u]; kwargs...) |
| 153 | +end |
| 154 | + |
| 155 | +""" |
| 156 | + get_comp_sensitivity(sys, ap::AnalysisPoint; kwargs) |
| 157 | + get_comp_sensitivity(sys, ap_name::Symbol; kwargs) |
| 158 | +
|
| 159 | +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`. |
| 160 | +
|
| 161 | +# Arguments: |
| 162 | +- `kwargs`: Are sent to `ModelingToolkit.linearize` |
| 163 | +
|
| 164 | +See also [`get_sensitivity`](@ref). |
| 165 | +""" |
| 166 | +function get_comp_sensitivity(sys, ap::AnalysisPoint; kwargs...) |
| 167 | + t = get_iv(sys) |
| 168 | + @variables d(t) = 0 # Perturbantion serving as input to sensivity transfer function |
| 169 | + new_eqs = map(get_eqs(sys)) do eq |
| 170 | + eq.rhs == ap || (return eq) |
| 171 | + ap.out.u + d ~ ap.in.u # This assumes that the connector as an internal vaiable named u |
| 172 | + end |
| 173 | + @set! sys.eqs = new_eqs |
| 174 | + @set! sys.states = [states(sys); d] |
| 175 | + ModelingToolkit.linearize(sys, [d], [ap.in.u]; kwargs...) |
| 176 | +end |
| 177 | + |
| 178 | +# Add a method to get_sensitivity that accepts the name of an AnalysisPoint |
| 179 | +for f in [:get_sensitivity, :get_comp_sensitivity] |
| 180 | + @eval function $f(sys, ap_name::Symbol, args...; kwargs...) |
| 181 | + $f(sys, find_analysis_point(sys, ap_name), args...; kwargs...) |
| 182 | + end |
| 183 | +end |
0 commit comments