A few bugs with MonteCarloMeasurements + ControlSystems #95
Comments
|
Using MCM for ControlSystems is not trivial due to the many factorizations that often have to be done. I often define function LinearAlgebra.eigvals(p::Matrix{<:AbstractParticles{T,N}}; kwargs...) where {T,N} # Special case to propte types differently
individuals = map(1:length(p[1])) do i
eigvals(getindex.(p,i); kwargs...)
end
PRT = Complex{Particles{T,N}}
out = Vector{PRT}(undef, length(individuals[1]))
for i = eachindex(out)
c = getindex.(individuals,i)
out[i] = complex(Particles{T,N}(real(c)),Particles{T,N}(imag(c)))
end
out
endbut one must to be very careful with interpreting the results of this computation, in particular when the uncertainty is such that the eigenvalues jump around and switch places in the output vector from I don't know how matlab's The "unreachable reached" appears to be some kind of julia bug, unless I have made something super sketchy in MCM. It's probably a good idea to since I could reproduce the error here. |
|
Alright, thanks. Yeah, some people at work are starting to show interest in Julia and one of them asked me if we could replicate MATLAB's Robust Control Toolbox, so I figured I'd go through some of their examples and see how far I could get. Most of what I do with ControlSystems+MCM works though and makes for a pretty cool demo (especially with Plots). |
|
I have looked into this a bit more (it's remarkably close to what I do at work :) |
|
Here are some attempts julia> using ControlSystems
julia> using GenericLinearAlgebra
julia> using MonteCarloMeasurements
julia> using Plots
julia> unsafe_comparisons(true, verbose=false)
false
julia> s = tf("s");
julia> k = 1.0 + (-0.2 .. 0.2)
Particles{Float64, 2000}
1.0 ± 0.116
julia> m1 = 1.0 + (-0.2 .. 0.2)
Particles{Float64, 2000}
1.0 ± 0.116
julia> m2 = 1.0 + (-0.2 .. 0.2)
Particles{Float64, 2000}
1.0 ± 0.116
julia> function LinearAlgebra.eigvals(p::Matrix{<:AbstractParticles{T,N}}; kwargs...) where {T,N} # Special case to propte types differently
individuals = map(1:length(p[1])) do i
eigvals(getindex.(p,i); kwargs...)
end
PRT = Complex{Particles{T,N}}
out = Vector{PRT}(undef, length(individuals[1]))
for i = eachindex(out)
c = getindex.(individuals,i)
out[i] = complex(Particles{T,N}(real(c)),Particles{T,N}(imag(c)))
end
out
end
julia> function ControlSystems.StateSpace(A::Matrix,B,C,D, ::Continuous, ::Int64, ::Int64, ::Int64)
ss(A,B,C,D)
end
julia> function dmm(k,m1,m2)
c0 = c1 = c2 = 0
A = [
0.0 1 0 0
-k/m1 -(c1 + c0)/m1 k/m1 c1/m1
0 0 0 1
k/m2 c1/m2 -k/m2 -(c1 + c2)/m2
]
B = [0, 1, 0, 0]
C = [0 0 1 0]
sys = ss(A, B, C, 0)
end
dmm
julia> P = dmm(k,m1,m2)
StateSpace{Continuous, Particles{Float64, 2000}, Matrix{Particles{Float64, 2000}}}
A =
0.0 1.0 0.0 0.0
-1.01 ± 0.17 0.0 1.01 ± 0.17 0.0
0.0 0.0 0.0 1.0
1.01 ± 0.16 0.0 -1.01 ± 0.16 0.0
B =
0.0
1.0
0.0
0.0
C =
0.0 0.0 1.0 0.0
D =
0.0
Continuous-time state-space model
julia> C = 100 * ss((s + 1) / (0.001*s + 1))^3
StateSpace{Continuous, Float64, Matrix{Float64}}
A =
-1000.0 -999000.0 -9.99e8
0.0 -1000.0 -999000.0
0.0 0.0 -1000.0
B =
1.0e6
1000.0
1.0
C =
-9.99e7 -9.99e10 -9.99e13
D =
1.0e11
Continuous-time state-space model
julia> L = (P * C)
StateSpace{Continuous, Particles{Float64, 2000}, Matrix{Particles{Float64, 2000}}}
A =
0.0 1.0 0.0 0.0 0.0 0.0 0.0
-1.01 ± 0.17 0.0 1.01 ± 0.17 0.0 -9.99e7 -9.99e10 -9.99e13
0.0 0.0 0.0 1.0 0.0 0.0 0.0
1.01 ± 0.16 0.0 -1.01 ± 0.16 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 -1000.0 -999000.0 -9.99e8
0.0 0.0 0.0 0.0 0.0 -1000.0 -999000.0
0.0 0.0 0.0 0.0 0.0 0.0 -1000.0
B =
0.0
1.0e11
0.0
0.0
1.0e6
1000.0
1.0
C =
0.0 0.0 1.0 0.0 0.0 0.0 0.0
D =
0.0
Continuous-time state-space model
julia> T = feedback(L)
StateSpace{Continuous, Particles{Float64, 2000}, Matrix{Particles{Float64, 2000}}}
A =
0.0 1.0 0.0 0.0 0.0 0.0 0.0
-1.01 ± 0.17 0.0 -1.0e11 ± 0.17 0.0 -9.99e7 -9.99e10 -9.99e13
0.0 0.0 0.0 1.0 0.0 0.0 0.0
1.01 ± 0.16 0.0 -1.01 ± 0.16 0.0 0.0 0.0 0.0
0.0 0.0 -1.0e6 0.0 -1000.0 -999000.0 -9.99e8
0.0 0.0 -1000.0 0.0 0.0 -1000.0 -999000.0
0.0 0.0 -1.0 0.0 0.0 0.0 -1000.0
B =
0.0
1.0e11
0.0
0.0
1.0e6
1000.0
1.0
C =
0.0 0.0 1.0 0.0 0.0 0.0 0.0
D =
0.0
Continuous-time state-space model
julia> w = exp10.(LinRange(-1, 1, 600));
julia> function worst_case_gain(P, w)
g,p = bode(P,w)
g = maximum.(g)
wg,wi = findmax(g, dims=1)
@info "The worst-case gain is $wg at freq. $(w[wi]) rad/s"
wg,w[wi]
end
worst_case_gain
julia> worst_case_gain(T, w)
[ Info: The worst-case gain is [1.039958463247935] at freq. [7.352653051278958] rad/s
([1.03996], [7.35265])
julia> Pnom = MonteCarloMeasurements.mean_object(P)
StateSpace{Continuous, Float64, Matrix{Float64}}
A =
0.0 1.0 0.0 0.0
-1.0139854186529447 0.0 1.0139854186529447 0.0
0.0 0.0 0.0 1.0
1.0132002552961126 0.0 -1.0132002552961126 0.0
B =
0.0
1.0
0.0
0.0
C =
0.0 0.0 1.0 0.0
D =
0.0
Continuous-time state-space model
julia> Tnom = MonteCarloMeasurements.mean_object(T)
StateSpace{Continuous, Float64, Matrix{Float64}}
A =
0.0 1.0 0.0 0.0 0.0 0.0 0.0
-1.0139854186529447 0.0 -9.999999999898604e10 0.0 -9.99e7 -9.99e10 -9.99e13
0.0 0.0 0.0 1.0 0.0 0.0 0.0
1.0132002552961126 0.0 -1.0132002552961126 0.0 0.0 0.0 0.0
0.0 0.0 -1.0e6 0.0 -1000.0 -999000.0 -9.99e8
0.0 0.0 -1000.0 0.0 0.0 -1000.0 -999000.0
0.0 0.0 -1.0 0.0 0.0 0.0 -1000.0
B =
0.0
1.0e11
0.0
0.0
1.0e6
1000.0
1.0
C =
0.0 0.0 1.0 0.0 0.0 0.0 0.0
D =
0.0
Continuous-time state-space model
julia> dcgain(Tnom)[]
1.0
julia> dcgain(T)[]
Particles{Float64, 2000}
1.0 ± 1.0e-7
julia> maximum(maximum.(real.(pole(T))))
-0.807577
julia> gm,pm, wg, wp = margin(Tnom)
([NaN], [Inf], [39.044], [157.783])
julia> function robstab(P)
margins = map(1:MonteCarloMeasurements.nparticles(k)) do i
P0 = MonteCarloMeasurements.replace_particles(P, replacer=p->p[i])
margin(tf(P0))
end
worst_gainm = minimum(reduce(vcat, getindex.(margins, 1)))
worst_phasem = minimum(reduce(vcat, getindex.(margins, 2)))
(;worst_gainm, worst_phasem)
end
robstab
julia> robstab(Tnom)
(worst_gainm = 571.722, worst_phasem = 7.62437) |
|
Nice! Thanks for this. |
|
I'll close this issue in favor of JuliaControl/ControlSystems.jl#396 |
|
@jonniedie I'm not sure if you've seen, but RobustAndOptimalControl.jl has seen some developement recently |
|
@baggepinnen Awesome, thanks! I'll check it out. |
Hey, sorry to dump multiple things on here, I can separate these into different issues if you'd like. I'm not even sure if this or ControlSystems.jl is the right place for this kind of thing. I'm trying to use MonteCarloMeasurements.jl + ControlSystems.jl to get something close to MATLAB's Robust Control Toolbox. Trying to replicate this example fails on a few parts.
bodeon F, or G1 crashes my system sometimes.This gives me the error
What's weird is if I put some other stuff before the
bodecalls, it works forG1.T,P, orLtozpkform hits a stack overflow error. This one, I'm sure, is more on the ControlSystems side though. It looks like it's expecting the eltypes of the state space matrices to beAbstractFloats. I tried converting totffirst, but that gives me an error that the iteration limit for a schur decomposition is reached.What's weird here is that I can call
eigvals(T.A)and it works fine.The text was updated successfully, but these errors were encountered: