search

baggepinnen / MonteCarloMeasurements.jl

Propagation of distributions by Monte-Carlo sampling: Real number types with uncertainty represented by samples.
https://baggepinnen.github.io/MonteCarloMeasurements.jl/stable/
MIT License
261 stars 17 forks source link

A few bugs with MonteCarloMeasurements + ControlSystems #95

Closed jonniedie closed 3 years ago

jonniedie commented 3 years ago

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.

  1. Here, trying to call bode on F, or G1 crashes my system sometimes. 
    
using ControlSystems
using GenericLinearAlgebra
using MonteCarloMeasurements
using Plots

unsafe_comparisons(true, verbose=false)

s = tf("s")

C = 100 ss((s + 1) / (0.001s + 1))^3

k = 1.0 ± 0.2 m1 = 1.0 ± 0.2 m2 = 1.0 ± 0.2

G1 = 1 / s^2 / m1 G2 = 1 / s^2 / m2

F = [0; G1] [1 -1] + [1 -1]' [0 G2]

bode(G1) # This crashes bode(F) # This crashes bode(F[2,1]) # This also crashes

This gives me the error

Unreachable reached at 000000003CBF36DF

Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks. Exception: EXCEPTION_ILLEGAL_INSTRUCTION at 0x3cbf36df -- getindex at .\array.jl:809 in expression starting at REPL[2]:1 getindex at .\array.jl:809 iterate at .\array.jl:785 [inlined] iterate at .\array.jl:785 [inlined] iterate at .\generator.jl:44 [inlined] _collect at .\array.jl:699 [inlined] collect_similar at .\array.jl:628 [inlined] map at .\abstractarray.jl:2162 [inlined] zpkdata at C:\Users\jdiegelm.julia\packages\ControlSystems\H3Br2\src\analysis.jl:123 _bounds_and_features at C:\Users\jdiegelm.julia\packages\ControlSystems\H3Br2\src\freqresp.jl:152

146 at C:\Users\jdiegelm.julia\packages\ControlSystems\H3Br2\src\freqresp.jl:138

iterate at .\generator.jl:47 [inlined] _collect at .\array.jl:699 [inlined] collect_similar at .\array.jl:628 unknown function (ip: 000000003CBF3D64) map at .\abstractarray.jl:2162 [inlined] _default_freq_vector at C:\Users\jdiegelm.julia\packages\ControlSystems\H3Br2\src\freqresp.jl:138 _default_freq_vector at C:\Users\jdiegelm.julia\packages\ControlSystems\H3Br2\src\freqresp.jl:145 [inlined] bode at C:\Users\jdiegelm.julia\packages\ControlSystems\H3Br2\src\freqresp.jl:108 unknown function (ip: 000000003CBF3640) jl_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\julia.h:1690 [inlined] do_call at /cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:117 eval_value at /cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:206 eval_stmt_value at /cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:157 [inlined] eval_body at /cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:548 jl_interpret_toplevel_thunk at /cygdrive/d/buildbot/worker/package_win64/build/src\interpreter.c:660 jl_toplevel_eval_flex at /cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:840 jl_toplevel_eval_flex at /cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:790 jl_toplevel_eval at /cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:849 [inlined] jl_toplevel_eval_in at /cygdrive/d/buildbot/worker/package_win64/build/src\toplevel.c:883 eval at .\boot.jl:331 eval_user_input at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:134 repl_backend_loop at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:195 start_repl_backend at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:180

run_repl#37 at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:292

run_repl at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\REPL\src\REPL.jl:288

807 at .\client.jl:399

jfptr_YY.807_38807 at C:\Users\jdiegelm\AppData\Local\Programs\Julia\Julia 1.5.2\lib\julia\sys.dll (unknown line) jl_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\julia.h:1690 [inlined] do_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\builtins.c:655 jl_fapply at /cygdrive/d/buildbot/worker/package_win64/build/src\builtins.c:669 [inlined] jl_fapply_latest at /cygdrive/d/buildbot/worker/package_win64/build/src\builtins.c:705

invokelatest#1 at .\essentials.jl:710 [inlined]

invokelatest at .\essentials.jl:709 [inlined] run_main_repl at .\client.jl:383 exec_options at .\client.jl:313 _start at .\client.jl:506 jfptrstart_36391 at C:\Users\jdiegelm\AppData\Local\Programs\Julia\Julia 1.5.2\lib\julia\sys.dll (unknown line) jl_apply at /cygdrive/d/buildbot/worker/package_win64/build/src\julia.h:1690 [inlined] true_main at /cygdrive/d/buildbot/worker/package_win64/build/ui\repl.c:106 wmain at /cygdrive/d/buildbot/worker/package_win64/build/ui\repl.c:227 tmainCRTStartup at /usr/src/debug/mingw64-x86_64-runtime-7.0.0-1/crt\crtexe.c:334 mainCRTStartup at /usr/src/debug/mingw64-x86_64-runtime-7.0.0-1/crt\crtexe.c:223 BaseThreadInitThunk at C:\Windows\System32\KERNEL32.DLL (unknown line) RtlUserThreadStart at C:\Windows\SYSTEM32\ntdll.dll (unknown line) Allocations: 114178914 (Pool: 114140923; Big: 37991); GC: 110

What's weird is if I put some other stuff before the `bode` calls, it works for `G1`.
```julia
using ControlSystems
using GenericLinearAlgebra
using MonteCarloMeasurements
using Plots

unsafe_comparisons(true, verbose=false)

s = tf("s")

k = 1.0 ± 0.2
m1 = 1.0 ± 0.2
m2 = 1.0 ± 0.2

G1 = 1 / (m1 * s^2)
G2 = 1 / (m2 * s^2)

F = [0; G1] * [1 -1] + [1 -1]' * [0 G2]

P = lft(F, ss(k))

C = 100 * ss((s + 1) / (0.001*s + 1))^3

L = P * C

T = feedback(L, 1)

bode(G1)
  1. Trying to convert T, P, or L to zpk form 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 be AbstractFloats. I tried converting to tf first, but that gives me an error that the iteration limit for a schur decomposition is reached. 
    julia> tf(T)
ERROR: ArgumentError: iteration limit 330 reached
Stacktrace:
[1] _schur!(::GenericLinearAlgebra.HessenbergFactorization{Particles{Float64,2000},Array{Particles{Float64,2000},2},GenericLinearAlgebra.Householder{Particles{Float64,2000},S} where S<:(StridedArray{T, 1} where 
T)}; tol::Float64, debug::Bool, shiftmethod::Symbol, maxiter::Int64, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:104
[2] _schur! at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:91 [inlined]
[3] #_schur!#21 at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:169 [inlined]
[4] _schur! at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:169 [inlined]
[5] #_eigvals!#23 at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:246 [inlined]
[6] _eigvals! at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:246 [inlined]
[7] eigvals!(::Array{Particles{Float64,2000},2}; sortby::Nothing, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at C:\Users\jdiegelm\.julia\packages\GenericLinearAlgebra\JoKjw\src\eigenGeneral.jl:260
[8] #eigvals#73 at D:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.5\LinearAlgebra\src\eigen.jl:326 [inlined]
[9] #eigvalsnosort#82 at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\utilities.jl:35 [inlined]
[10] eigvalsnosort at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\utilities.jl:35 [inlined]
[11] charpoly(::Array{Particles{Float64,2000},2}) at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\types\conversion.jl:310
[12] convert(::Type{TransferFunction{Continuous,ControlSystems.SisoRational{Particles{Float64,2000}}}}, ::StateSpace{Continuous,Particles{Float64,2000},Array{Particles{Float64,2000},2}}) at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\types\conversion.jl:261
[13] convert(::Type{TransferFunction{Continuous,ControlSystems.SisoRational}}, ::StateSpace{Continuous,Particles{Float64,2000},Array{Particles{Float64,2000},2}}) at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\types\conversion.jl:270
[14] convert at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\types\conversion.jl:250 [inlined]
[15] tf(::StateSpace{Continuous,Particles{Float64,2000},Array{Particles{Float64,2000},2}}) at C:\Users\jdiegelm\.julia\packages\ControlSystems\H3Br2\src\types\tf.jl:69

    What's weird here is that I can call eigvals(T.A) and it works fine.

baggepinnen commented 3 years ago

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
end

but 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 eigvals.

I don't know how matlab's ureal works, but it could be interesting to look into. I often exactly what you are trying to with success though, so I wonder if there is something on my fork of ControlSystems that is required to make it work. I'll see if I can look into it.

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

Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks.

since I could reproduce the error here.

jonniedie commented 3 years ago

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).

baggepinnen commented 3 years ago

I have looked into this a bit more (it's remarkably close to what I do at work :) We suffer from pole-zero cancellations not being carried out when systems contain uncertain parameters. I'm not sure how to get around this in an effective way. The matlab ureal appears to be some kind of a symbolic variable, so I tried to work for a while with SymPy and SymbolicControlSystems.jl but didn't really get very far. SymPy bogs down very quickly so by the time I had gotten to the expression for T it was enormous. I think we need better tools for simplification in general, our minreal is not really working that well even for Float64 :/

baggepinnen commented 3 years ago

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)
jonniedie commented 3 years ago

Nice! Thanks for this.

baggepinnen commented 3 years ago

I'll close this issue in favor of https://github.com/JuliaControl/ControlSystems.jl/issues/396 I'm not sure how to approach the interface between MCM and CS in general, I have a feeling it would boil down to a new package that explicitly targets the combination of MCM and CS and implements the functions required for robust control. I hope that most tools required to build such a package are present in MCM already, but several functions like robstab above are probably required in order to make robust control synthesis and analysis convenient. Feel free to open an issue to discuss such a package over at ControlSystems (JuliaControl is probably a better place for such a discussion?).

baggepinnen commented 2 years ago

@jonniedie I'm not sure if you've seen, but RobustAndOptimalControl.jl has seen some developement recently https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/ The uncertainty modelling is still somewhat of a weak spot, real parameter uncertainty (ureal) in particular, but we do handle basic μ-analysis (no μ-synthesis) and complex uncertainties as well as disk margins etc.

jonniedie commented 2 years ago

@baggepinnen Awesome, thanks! I'll check it out.