-
Notifications
You must be signed in to change notification settings - Fork 6
/
pso_model.jl
232 lines (207 loc) · 8.82 KB
/
pso_model.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
using UniqueVectors: UniqueVector
using SparseArrays: spzeros, spdiagm
using LinearAlgebra: I, eigen, svd, diagm
using IterativeSolvers: lobpcg
export PSOModel
export lossless_modes_dense, lossy_modes_dense, lossless_modes_sparse
struct PSOModel{T,
K<:AbstractMatrix{<:Number},
G<:AbstractMatrix{<:Number},
C<:AbstractMatrix{<:Number},
P<:AbstractMatrix{<:Number}} <: AdmittanceModel{T}
K::K
G::G
C::C
P::P
ports::UniqueVector{T}
function PSOModel{T}(k, g, c, p, ports::UniqueVector{T}) where {T}
l = size(k, 1)
@assert all(isreal(m) for m in (k, g, c, p))
@assert size(p) == (l, length(ports))
@assert all(size(y) == (l, l) for y in (k, g, c))
kk, gg, cc, pp = float.((k, g, c, p))
K, G, C, P = typeof.((kk, gg, cc, pp))
return new{T,K,G,C,P}(kk, gg, cc, pp, ports)
end
end
"""
PSOModel(K, G, C, P, ports::AbstractVector{T}) where {T}
Create a positive second order (PSO) model, representing the equations
`(K + G∂ₜ + C∂²ₜ)Φ = Px`, `y = Qᵀ∂ₜΦ`. (If the subscript `t` does not display in
the equations, install a font with greater Unicode coverage.)
"""
PSOModel(K, G, C, P, ports::AbstractVector{T}) where {T} =
PSOModel{T}(K, G, C, P, UniqueVector(ports))
get_Y(pso::PSOModel) = (pso.K, pso.G, pso.C)
get_P(pso::PSOModel) = pso.P
get_ports(pso::PSOModel) = pso.ports
function partial_copy(pso::PSOModel;
Y::Union{AbstractVector, Nothing} = nothing,
P::Union{AbstractMatrix{<:Real}, Nothing} = nothing,
ports::Union{AbstractVector, Nothing} = nothing)
Y = isnothing(Y) ? get_Y(pso) : Y
P = isnothing(P) ? get_P(pso) : P
ports = isnothing(ports) ? get_ports(pso) : ports
return PSOModel(Y..., P, UniqueVector(ports))
end
compatible(psos::AbstractVector{<:PSOModel}) = true
function canonical_gauge(pso::PSOModel)
fullP = collect(pso.P)
F = qr(fullP)
(m, n) = size(fullP)
transform = (F.Q * Matrix(I,m,m)) * [transpose(inv(F.R)) zeros(n, m-n);
zeros(m-n, n) Matrix(I, m-n, m-n)]
return apply_transform(pso, transform)
end
# NOTE we cannot write lossy_modes_sparse because KrylovKit.geneigsolve
# solves Av = λBv but requires that B is positive semidefinite and A is
# symmetric or Hermitian. For lossy modes, we can make B positive semidefinite
# and A nonsymmetric or A symmetric and B not positive semidefinite. If
# KrylovKit gains the ability to solve either type of eigenproblem we can
# write lossy_modes_sparse.
"""
lossy_modes_dense(pso::PSOModel; min_freq::Real=0,
max_freq::Union{Real, Nothing}=nothing, real_tol::Real=Inf, imag_tol::Real=Inf)
Use a dense eigenvalue solver to find the decaying modes of the PSOModel. Each mode consists
of an eigenvalue `λ` and a spatial distribution vector `v`. The frequency of the mode is
`imag(λ)/(2π)` and the decay rate is `-2*real(λ)/(2π)`. `min_freq` and `max_freq` give the
frequency range in which to find eigenvalues. Passivity is enforced using an inverse power
method solver, and `real_tol` and `imag_tol` are the tolerances used in that method to
determine convergence of the real and imaginary parts of the eigenvalue.
"""
function lossy_modes_dense(pso::PSOModel; min_freq::Real=0,
max_freq::Union{Real, Nothing}=nothing, real_tol::Real=Inf, imag_tol::Real=Inf)
K, G, C = get_Y(pso)
Z = zero(K)
M = collect([C Z; Z I]) \ collect([-G -K; I Z])
M_scale = maximum(abs.(M))
# for some reason eigen(Y \ X) seems to perform better than eigen(X, Y)
values, vectors = eigen(M/M_scale)
values *= M_scale
freqs = imag.(values)/(2π)
if max_freq == nothing
max_freq = maximum(freqs)
end
inds = findall(max_freq .>= freqs .>= min_freq)
values, vectors = values[inds], vectors[:, inds]
# enforce passivity using inverse power method
pred(x, y) = ((real(y) <= 0) && (abs(real(x-y)) < real_tol) &&
(abs(imag(x-y)) < imag_tol))
for i in 1:length(values)
if real(values[i]) > 0
λ, v = inv_power_eigen(M, v0=vectors[:, i], λ0=values[i], pred=pred)
values[i] = λ
vectors[:, i] = v
end
end
# the second half of the eigenvector is the spatial part
vectors = vectors[size(K, 1)+1:end, :]
# normalize the columns
return values, vectors./sqrt.(sum(abs2.(vectors), dims=1))
end
"""
lossless_modes_dense(pso::PSOModel; min_freq::Real=0,
max_freq::Union{Real, Nothing}=nothing)
Use a dense svd solver to find the modes of the PSOModel neglecting loss. Each mode consists
of an eigenvalue `λ` and an eigenvector `v`. The frequency of the mode is `imag(λ)/(2π)`
and the decay rate is `-2*real(λ)/(2π)`. `min_freq` and `max_freq` give the
frequency range in which to find eigenvalues.
"""
function lossless_modes_dense(pso::PSOModel; min_freq::Real = 0,
max_freq::Union{Real, Nothing} = nothing)
K, _, C = get_Y(pso)
K = .5 * (transpose(K) + K)
C = .5 * (transpose(C) + C)
sqrt_C, sqrt_K = sqrt(collect(C)), sqrt(collect(K))
M = sqrt_C \ sqrt_K
scale = maximum(abs.(M))
U, S, Vt = svd(M/scale)
freqs = S * scale/(2π)
vectors = real.(sqrt_C \ U)
if max_freq == nothing
max_freq = maximum(freqs)
end
inds = findall(max_freq .>= freqs .>= min_freq)
freqs, vectors = freqs[inds], vectors[:, inds]
inds = sortperm(freqs)
return freqs[inds] * 2π * 1im, vectors[:, inds]
end
"""
lossless_modes_sparse(pso::PSOModel; num_modes = 1, maxiter = 200)
Use a sparse generalized eigenvalue solver to find the `num_modes` lowest frequency modes
of the PSOModel neglecting loss. Each mode consists of an eigenvalue `λ` and an eigenvector
`v`. The frequency of the mode is `imag(λ)/(2π)` and the decay rate is `-2*real(λ)/(2π)`.
"""
function lossless_modes_sparse(pso::PSOModel; num_modes = 1, maxiter = 200)
K, _, C = get_Y(pso)
K = (K + transpose(K))/2
C = (C + transpose(C))/2
K_scale, C_scale = maximum(abs.(K)), maximum(abs.(C))
result = lobpcg(K/K_scale, C/C_scale, false, num_modes, maxiter=maxiter)
values = result.λ * (K_scale/C_scale)
vectors = result.X
return real.(sqrt.(values)) * 1im, vectors
end
"""
port_matrix(c::Circuit, port_edges::Vector{<:Tuple})
Returns a port matrix suitable for constructing a `PSOModel` based on
`c::Circuit`, where the columns of the port matrix correspond to edges between
circuit vertices. The edges are given as tuples of vertex names in `port_edges`.
"""
function port_matrix(c::Circuit, port_edges::Vector{<:Tuple})
coord_matrix = coordinate_matrix(c)
if length(port_edges) > 0
port_indices = [(findfirst(isequal(p[1]), c.vertices),
findfirst(isequal(p[2]), c.vertices)) for p in port_edges]
return hcat([coord_matrix[:, i1] - coord_matrix[:, i2]
for (i1, i2) in port_indices]...)
else
return spzeros(size(coord_matrix, 1), 0)
end
end
"""
circuit_to_pso_matrix(m::AbstractMatrix)
Convert one of the inverse inductance, conductance, or capacitance matrices
characterizing a `Circuit` into the form needed for a `PSOModel`.
"""
circuit_to_pso_matrix(m::AbstractMatrix) =
(-m + spdiagm(0=>sum(m, dims=2)[:,1]))[2:end, 2:end]
"""
PSOModel(circuit::Circuit, port_edges::Vector{<:Tuple},
port_names::AbstractVector)
Create a PSOModel for a circuit with ports on given edges. The spanning tree
where the first vertex is chosen as ground and all other vertices are neighbors
of ground is used.
"""
function PSOModel(circuit::Circuit, port_edges::Vector{<:Tuple},
port_names::AbstractVector)
Y = [circuit_to_pso_matrix(m) for m in matrices(circuit)]
P = port_matrix(circuit, port_edges)
return PSOModel(Y..., P, port_names)
end
"""
coordinate_transform(m::AbstractMatrix)
Construct a coordinate transformation matrix for a `PSOModel` from a matrix
given by e.g. `coordinate_matrix(c::Circuit{T}, tree::SpanningTree{T}) where T`.
"""
coordinate_transform(m::AbstractMatrix) = m[:,2:end] .- m[:,1]
"""
coordinate_transform(coord_from::AbstractMatrix, coord_to::AbstractMatrix)
Construct a coordinate transformation matrix to go from a `PSOModel` constructed
with one spanning tree to a `PSOModel` constructed with another.
"""
function coordinate_transform(coord_from::AbstractMatrix, coord_to::AbstractMatrix)
return coordinate_transform(coord_to) / coordinate_transform(coord_from)
end
"""
PSOModel(circuit::Circuit, port_edges::Vector{<:Tuple},
port_names::AbstractVector, tree::SpanningTree)
Create a PSOModel for a circuit with ports on given edges and using the given
spanning tree.
"""
function PSOModel(circuit::Circuit, port_edges::Vector{<:Tuple},
port_names::AbstractVector, tree::SpanningTree)
pso = PSOModel(circuit, port_edges, port_names)
transform = transpose(coordinate_transform(coordinate_matrix(circuit, tree)))
return apply_transform(pso, transform)
end