Skip to content

Commit 34a93aa

Browse files
baggepinnenbaggepinnen
committed
add analysis points
1 parent 56eb837 commit 34a93aa

File tree

3 files changed

+244
-0
lines changed

3 files changed

+244
-0
lines changed

‎src/Blocks/Blocks.jl

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,4 +26,7 @@ export Integrator, Derivative, FirstOrder, SecondOrder, StateSpace
2626
export PI, LimPI, PID, LimPID
2727
include("continuous.jl")
2828

29+
export AnalysisPoint, expand_analysis_points, get_sensitivity, get_comp_sensitivity
30+
include("analysis_points.jl")
31+
2932
end

‎src/Blocks/analysis_points.jl

Lines changed: 179 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,179 @@
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+
ModelingToolkit.vars(ap::AnalysisPoint; op = Differential) = vars(connect(ap.in, ap.out); op)
102+
103+
"""
104+
find_analysis_point(sys, name::Symbol)
105+
106+
Find and return the analysis point in `sys` with the specified `name`. If no matching [`AnalysisPoint`](@ref) is found, `nothing` is returned.
107+
"""
108+
function find_analysis_point(sys, name)
109+
for eq in equations(sys)
110+
eq.rhs isa AnalysisPoint && eq.rhs.name == name && (return eq.rhs)
111+
end
112+
nothing
113+
end
114+
115+
"""
116+
expand_analysis_points(sys)
117+
118+
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.
119+
"""
120+
function expand_analysis_points(sys)
121+
new_eqs = map(get_eqs(sys)) do eq
122+
eq.rhs isa AnalysisPoint || (return eq)
123+
ap = eq.rhs
124+
connect(ap.in, ap.out)
125+
end
126+
@set! sys.eqs = new_eqs
127+
sys
128+
end
129+
130+
"""
131+
get_sensitivity(sys, ap::AnalysisPoint; kwargs)
132+
get_sensitivity(sys, ap_name::Symbol; kwargs)
133+
134+
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`.
135+
136+
# Arguments:
137+
- `kwargs`: Are sent to `ModelingToolkit.linearize`
138+
139+
See also [`get_comp_sensitivity`](@ref).
140+
"""
141+
function get_sensitivity(sys, ap::AnalysisPoint; kwargs...)
142+
t = get_iv(sys)
143+
@variables d(t)=0 # Perturbantion serving as input to sensivity transfer function
144+
new_eqs = map(get_eqs(sys)) do eq
145+
eq.rhs == ap || (return eq)
146+
ap.out.u ~ ap.in.u + d # This assumes that the connector as an internal vaiable named u
147+
end
148+
@set! sys.eqs = new_eqs
149+
@set! sys.states = [states(sys); d]
150+
ModelingToolkit.linearize(sys, [d], [ap.out.u]; kwargs...)
151+
end
152+
153+
"""
154+
get_comp_sensitivity(sys, ap::AnalysisPoint; kwargs)
155+
get_comp_sensitivity(sys, ap_name::Symbol; kwargs)
156+
157+
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`.
158+
159+
# Arguments:
160+
- `kwargs`: Are sent to `ModelingToolkit.linearize`
161+
162+
See also [`get_sensitivity`](@ref).
163+
"""
164+
function get_comp_sensitivity(sys, ap::AnalysisPoint; kwargs...)
165+
t = get_iv(sys)
166+
@variables d(t)=0 # Perturbantion serving as input to sensivity transfer function
167+
new_eqs = map(get_eqs(sys)) do eq
168+
eq.rhs == ap || (return eq)
169+
ap.out.u + d ~ ap.in.u # This assumes that the connector as an internal vaiable named u
170+
end
171+
@set! sys.eqs = new_eqs
172+
@set! sys.states = [states(sys); d]
173+
ModelingToolkit.linearize(sys, [d], [ap.in.u]; kwargs...)
174+
end
175+
176+
# Add a method to get_sensitivity that accepts the name of an AnalysisPoint
177+
for f in [:get_sensitivity, :get_comp_sensitivity]
178+
@eval $f(sys, ap_name::Symbol, args...; kwargs...) = $f(sys, find_analysis_point(sys, ap_name), args...; kwargs...)
179+
end

‎test/Blocks/test_analysis_points.jl

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
using ModelingToolkit
2+
using ModelingToolkitStandardLibrary.Blocks
3+
using OrdinaryDiffEq
4+
using ModelingToolkit: get_eqs, vars, @set!, get_iv
5+
6+
@named P = FirstOrder(k=1, T=1)
7+
@named C = Gain(-1)
8+
t = ModelingToolkit.get_iv(P)
9+
10+
# Test with explicitly created AnalysisPoint
11+
ap = AnalysisPoint(:plant_input)
12+
eqs = [
13+
connect(P.output, C.input)
14+
connect(C.output, ap, P.input)
15+
]
16+
sys = ODESystem(eqs, t, systems=[P,C], name=:hej)
17+
18+
19+
ssys = structural_simplify(expand_analysis_points(sys))
20+
prob = ODEProblem(ssys, [P.x => 1], (0, 10))
21+
sol = solve(prob, Rodas5())
22+
@test norm(sol[1]) >= 1
23+
@test norm(sol[end]) < 1e-6 # This fails without the feedback through C
24+
# plot(sol)
25+
26+
matrices, _ = get_sensitivity(sys, ap)
27+
@test matrices.A[] == -2
28+
@test matrices.B[]*matrices.C[] == -1 # either one negative
29+
@test matrices.D[] == 1
30+
31+
32+
matrices, _ = get_comp_sensitivity(sys, ap)
33+
@test matrices.A[] == -2
34+
@test matrices.B[]*matrices.C[] == 1 # both positive or negative
35+
@test matrices.D[] == 0
36+
37+
#=
38+
# Equivalent code using ControlSystems. This can be used to verify the expected results tested for above.
39+
using ControlSystemsBase
40+
P = tf(1.0, [1, 1])
41+
C = 1 # Negative feedback assumed in ControlSystems
42+
S = sensitivity(P, C) # or feedback(1, P*C)
43+
T = comp_sensitivity(P, C) # or feedback(P*C)
44+
=#
45+
46+
# Test with automatically created analysis point
47+
eqs = [
48+
connect(P.output, C.input)
49+
connect(C.output, :plant_input, P.input)
50+
]
51+
sys = ODESystem(eqs, t, systems=[P,C], name=:hej)
52+
53+
matrices, _ = get_sensitivity(sys, :plant_input)
54+
@test matrices.A[] == -2
55+
@test matrices.B[]*matrices.C[] == -1 # either one negative
56+
@test matrices.D[] == 1
57+
58+
59+
matrices, _ = get_comp_sensitivity(sys, :plant_input)
60+
@test matrices.A[] == -2
61+
@test matrices.B[]*matrices.C[] == 1 # both positive
62+
@test matrices.D[] == 0

0 commit comments

Comments
 (0)