Skip to content

Commit aa1a83e

Browse files
baggepinnenbaggepinnen
committed
propagate extra
1 parent 690e7f9 commit aa1a83e

File tree

4 files changed

+105
-13
lines changed

4 files changed

+105
-13
lines changed

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
2626
[compat]
2727
ChainRulesCore = "1"
2828
ComponentArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15"
29-
ControlSystemsBase = "1.12"
29+
ControlSystemsBase = "1.17"
3030
DescriptorSystems = "1.2"
3131
Distributions = "0.25"
3232
GenericSchur = "0.5.2"

‎src/named_systems2.jl

Lines changed: 96 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -152,7 +152,42 @@ Return a named tuple `(; x, u)` containing the operating point of the system `P`
152152
"""
153153
operating_point(P::NamedStateSpace) = get_extra(P, :operating_point, (x=zeros(P.nx), u=zeros(P.nu)))
154154
operating_point(P::AbstractStateSpace) = (x=zeros(P.nx), u=zeros(P.nu))
155-
set_operating_point!(P::NamedStateSpace, xu::NamedTuple{(:x,:u)}) = set_extra!(P, :operating_point, xu)
155+
function set_operating_point!(P::NamedStateSpace, xu::NamedTuple{(:x,:u)})
156+
length(xu.x) == P.nx ||
157+
throw(ArgumentError("Operating point x must have length $(P.nx), got length $(length(xu.x))"))
158+
length(xu.u) == P.nu ||
159+
throw(ArgumentError("Operating point u must have length $(P.nu), got length $(length(xu.u))"))
160+
set_extra!(P, :operating_point, xu)
161+
end
162+
has_operating_point(P::NamedStateSpace) = haskey(P.extra, :operating_point)
163+
164+
165+
"""
166+
infer_operating_point(P1, P2, method = :obsv)
167+
168+
Return the operating point of `P1` inferred from the operating point of `P2` and a similarity transform between the two systems. The method for finding the similarity transform can be specified, default is `:obsv`.
169+
170+
```
171+
s1 = named_ss(ssrand(1,1,2), "P", x = :x, u = :u1, y=:y1)
172+
opx = randn(s1.nx)
173+
opu = randn(s1.nu)
174+
RobustAndOptimalControl.set_operating_point!(s1, (x = opx, u = opu))
175+
176+
s2, _ = balance_statespace(s1) # s1 and s2 are similar systems
177+
op2 = RobustAndOptimalControl.infer_operating_point(s2, s1)
178+
@test RobustAndOptimalControl.operating_point(s2) == op2
179+
```
180+
"""
181+
function infer_operating_point(P1, P2, method = :obsv)
182+
T = ControlSystemsBase.find_similarity_transform(P1, P2, method)
183+
op2 = operating_point(P2)
184+
op1 = (x = T*op2.x, u = op2.u)
185+
end
186+
187+
function infer_operating_point!(P1, P2, method = :obsv)
188+
op = infer_operating_point(P1, P2, method)
189+
set_operating_point!(P1, op)
190+
end
156191

157192
const NamedIndex = Union{Symbol, Vector{Symbol}, Colon}
158193

@@ -265,29 +300,43 @@ iterable(v) = v
265300
function Base.getindex(sys::NamedStateSpace{T,S}, i::NamedIndex, j::NamedIndex) where {T,S}
266301
i,j = iterable.((i, j))
267302
ii = i isa Colon ? i : names2indices(i, sys.y)
268-
jj = j isa Colon ? j : names2indices(j, sys.u)
303+
jj = j isa Colon ? j : names2indices(j, sys.u)
269304

270-
return NamedStateSpace{T,S}(
305+
newsys = NamedStateSpace{T,S}(
271306
sys.sys[ii, jj],
272307
sys.x,
273308
sys.u[jj],
274309
sys.y[ii],
275310
sys.name,
311+
sys.extra,
276312
) # |> sminreal
313+
if has_operating_point(sys)
314+
op = operating_point(sys)
315+
op = (; op.x, u=op.u[jj])
316+
set_operating_point!(newsys, op)
317+
end
318+
return newsys
277319
end
278320

279321
function Base.getindex(sys::NamedStateSpace{T,S}, inds...) where {T,S}
280322
if size(inds, 1) != 2
281323
error("Must specify 2 indices to index statespace model")
282324
end
283325
rows, cols = CS.index2range(inds...)
284-
return NamedStateSpace{T,S}(
326+
newsys = NamedStateSpace{T,S}(
285327
sys.sys[rows, cols],
286328
sys.x,
287329
sys.u[cols],
288330
sys.y[rows],
289331
sys.name,
332+
sys.extra,
290333
) # |> sminreal
334+
if has_operating_point(sys)
335+
op = operating_point(sys)
336+
op = (; op.x, u=op.u[cols])
337+
set_operating_point!(newsys, op)
338+
end
339+
return newsys
291340
end
292341

293342
function Base.show(io::IO, G::NamedStateSpace)
@@ -300,6 +349,10 @@ function Base.show(io::IO, G::NamedStateSpace)
300349
length(G.x) < 50 && (print(io, "\nWith state names: "); println(io, join(G.x, ' ')))
301350
length(G.u) < 50 && (print(io, " input names: "); println(io, join(G.u, ' ')))
302351
length(G.y) < 50 && (print(io, " output names: "); println(io, join(G.y, ' ')))
352+
if has_operating_point(G)
353+
op = operating_point(G)
354+
println(io, "Operating point: x = ", op.x, ", u = ", op.u)
355+
end
303356
end
304357

305358
function Base.:-(s1::NamedStateSpace{T,S}) where {T <: CS.TimeEvolution, S}
@@ -309,6 +362,7 @@ function Base.:-(s1::NamedStateSpace{T,S}) where {T <: CS.TimeEvolution, S}
309362
s1.u,
310363
s1.y,
311364
s1.name,
365+
s1.extra,
312366
)
313367
end
314368

@@ -319,6 +373,7 @@ function Base.:+(s1::NamedStateSpace{T,S}, s2::NamedStateSpace{T,S}) where {T <:
319373
s1.u,
320374
s1.y,
321375
"",
376+
s1.extra,
322377
)
323378
end
324379

@@ -336,6 +391,7 @@ function Base.:*(s1::NamedStateSpace{T}, s2::NamedStateSpace{T}) where {T <: CS.
336391
s2.u,
337392
s1.y,
338393
"",
394+
s1.extra,
339395
)
340396
end
341397

@@ -346,6 +402,7 @@ function Base.:*(s1::Number, s2::NamedStateSpace{T, S}) where {T <: CS.TimeEvolu
346402
s2.u,
347403
[Symbol(string(y)*"_scaled") for y in s2.y],
348404
isempty(s2.name) ? "" : s2.name*"_scaled",
405+
s2.extra,
349406
)
350407
end
351408

@@ -356,6 +413,7 @@ function Base.:*(s1::NamedStateSpace{T, S}, s2::Number) where {T <: CS.TimeEvolu
356413
[Symbol(string(u)*"_scaled") for u in s1.u],
357414
s1.y,
358415
isempty(s1.name) ? "" : s1.name*"_scaled",
416+
s1.extra,
359417
)
360418
end
361419

@@ -367,6 +425,7 @@ function Base.:*(s1::AbstractMatrix, s2::NamedStateSpace{T, S}) where {T <: CS.T
367425
s2.u,
368426
[Symbol(string(y)*"_scaled") for y in s2.y],
369427
isempty(s2.name) ? "" : s2.name*"_scaled",
428+
s2.extra,
370429
)
371430
else
372431
return *(promote(s1, s2)...)
@@ -381,6 +440,7 @@ function Base.:*(s1::NamedStateSpace{T, S}, s2::AbstractMatrix) where {T <: CS.T
381440
[Symbol(string(u)*"_scaled") for u in s1.u],
382441
s1.y,
383442
isempty(s1.name) ? "" : s1.name*"_scaled",
443+
s1.extra,
384444
)
385445
else
386446
return *(promote(s1, s2)...)
@@ -405,25 +465,45 @@ function Base.hcat(systems::NamedStateSpace{T,S}...) where {T,S}
405465
x = generate_unique_x_names(systems...)
406466
u = reduce(vcat, getproperty.(systems, :u))
407467
check_unique(u, "u")
468+
if any(has_operating_point, systems)
469+
opx = reduce(vcat, operating_point(sys).x for sys in systems)
470+
opu = reduce(vcat, operating_point(sys).u for sys in systems)
471+
op = (; x = opx, u = opu)
472+
extra = Dict(:operating_point => op)
473+
else
474+
extra = Dict{Symbol, Any}()
475+
end
408476
return NamedStateSpace{T,S}(
409477
hcat(getproperty.(systems, :sys)...),
410478
x,
411479
u,
412480
systems[1].y,
413481
"",
482+
extra
414483
)
415484
end
416485

417486
function Base.vcat(systems::NamedStateSpace{T,S}...) where {T,S}
418487
x = generate_unique_x_names(systems...)
419488
y = reduce(vcat, getproperty.(systems, :y))
420489
check_unique(y, "y")
490+
if any(has_operating_point, systems)
491+
allequal(operating_point(sys).u for sys in systems) ||
492+
throw(ArgumentError("All systems must have the same input operating point u to be concatenated."))
493+
opx = reduce(vcat, operating_point(sys).x for sys in systems)
494+
opu = operating_point(systems[1]).u
495+
op = (; x = opx, u = opu)
496+
extra = Dict(:operating_point => op)
497+
else
498+
extra = Dict{Symbol, Any}()
499+
end
421500
return NamedStateSpace{T,S}(
422501
vcat(getproperty.(systems, :sys)...),
423502
x,
424503
systems[1].u,
425504
y,
426505
"",
506+
extra,
427507
)
428508
end
429509

@@ -443,7 +523,7 @@ function measure(s::NamedStateSpace, names)
443523
C[i, inds[i]] = 1
444524
end
445525
s2 = ss(A,B,C,0, s.timeevol)
446-
sminreal(named_ss(s2; s.x, s.u, y=names, name=s.name))
526+
sminreal(named_ss(s2; s.x, s.u, y=names, name=s.name, s.extra))
447527
end
448528

449529
"""
@@ -541,7 +621,7 @@ function ControlSystemsBase.feedback(s1::NamedStateSpace{T}, s2::NamedStateSpace
541621
@assert sys.nu == length(W1) + length(W2)
542622
@assert sys.ny == length(Z1) + length(Z2)
543623
@assert sys.nx == length(x1)
544-
nsys = NamedStateSpace{T,typeof(sys)}(sys, x1, [s1.u[W1]; s2.u[W2]], [s1.y[Z1]; s2.y[Z2]], "")
624+
nsys = NamedStateSpace{T,typeof(sys)}(sys, x1, [s1.u[W1]; s2.u[W2]], [s1.y[Z1]; s2.y[Z2]], "", Dict{Symbol, Any}())
545625
# sminreal(nsys)
546626
end
547627

@@ -838,7 +918,7 @@ function named_ss(sys::ExtendedStateSpace{T}, name="";
838918
throw(ArgumentError("Length of performance names must match sys.nz ($(sys.nz))"))
839919

840920
sys2 = ss(sys)
841-
NamedStateSpace{T, typeof(sys2)}(sys2, x, [w; u], [z; y], name)
921+
NamedStateSpace{T, typeof(sys2)}(sys2, x, [w; u], [z; y], name, Dict{Symbol, Any}())
842922
end
843923

844924

@@ -893,7 +973,7 @@ end
893973

894974
function CS.c2d(s::NamedStateSpace{Continuous}, Ts::Real, method::Symbol = :zoh, args...;
895975
kwargs...)
896-
named_ss(c2d(s.sys, Ts, method, args...; kwargs...); s.x, s.u, s.y, s.name)
976+
named_ss(c2d(s.sys, Ts, method, args...; kwargs...); s.x, s.u, s.y, s.name, s.extra)
897977
end
898978

899979

@@ -910,8 +990,13 @@ function CS.append(systems::NamedStateSpace...; kwargs...)
910990
x = reduce(vcat, getproperty.(systems, :x))
911991
y = reduce(vcat, getproperty.(systems, :y))
912992
u = reduce(vcat, getproperty.(systems, :u))
913-
914-
return named_ss(systype(A, B, C, D, timeevol); x, y, u, kwargs...)
993+
994+
newsys = named_ss(systype(A, B, C, D, timeevol); x, y, u, kwargs...)
995+
if any(has_operating_point, systems)
996+
op = (; x = reduce(vcat, operating_point(s).x for s in systems), u = reduce(vcat, operating_point(s).u for s in systems))
997+
set_operating_point!(newsys, op)
998+
end
999+
return newsys
9151000
end
9161001

9171002

@@ -941,7 +1026,7 @@ function CS.minreal(sys::NamedStateSpace, args...; kwargs...)
9411026
named_ss(msys; sys.u, sys.y, sys.name)
9421027
end
9431028

944-
for fun in [:baltrunc, :balreal, :balance_statespace]
1029+
for fun in [:baltrunc, :balreal]
9451030
@eval function CS.$(fun)(sys::NamedStateSpace, args...; kwargs...)
9461031
msys, rest... = CS.$(fun)(sys.sys, args...; kwargs...)
9471032
named_ss(msys; sys.u, sys.y, sys.name), rest...

‎src/uncertainty_interface.jl

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -464,7 +464,9 @@ function ControlSystemsBase.poles(G::TransferFunction{<:ControlSystemsBase.TimeE
464464
end
465465

466466
function ControlSystemsBase.tzeros(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T}, D::AbstractMatrix{T}) where T <: AbstractParticles
467-
bymap(tzeros, A, B, C, D)
467+
all = [tzeros(vecindex(A, i), vecindex(B, i), vecindex(C, i), vecindex(D, i)) for i in 1:nparticles(A[1])]
468+
nz = length(all[1])
469+
[Particles(getindex.(all, i)) for i = 1:nz]
468470
end
469471

470472
function ControlSystemsBase.balance(A::AbstractMatrix{<:AbstractParticles}, perm::Bool=true)

‎test/test_named_systems2.jl

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -222,6 +222,11 @@ opu = randn(s1.nu)
222222
RobustAndOptimalControl.set_operating_point!(s1, (x = opx, u = opu))
223223
@test RobustAndOptimalControl.operating_point(s1) == (x = opx, u = opu)
224224

225+
s1b, T = balance_statespace(s1)
226+
opb = RobustAndOptimalControl.infer_operating_point(s1b, s1)
227+
228+
@test RobustAndOptimalControl.operating_point(s1b).x ≈ opb.x
229+
225230
## Promotion and conversion
226231
@testset "Promotion and conversion" begin
227232
@info "Testing Promotion and conversion"

0 commit comments

Comments
 (0)