MatrixPencils - a new feature

[andreasvarga]

With the last update of MatrixPencils.jl to v1.9.1, I am able to support in DescriptorSystems.jl transfer functions with uncertainties, via the Measurements.jl package. For example, the following example works (other probably not!):

using DescriptorSystems
using Measurements

 s = rtf("s")

J = 2.0*measurement(1,0.2)
b = 0.04*measurement(1,0.2)
K = 1.0*measurement(1,0.2)
R = 0.08*measurement(1,0.2)
L = 1e-4*measurement(1,0.2)

 P = K/(s*((J*s + b)*(L*s + R) + K^2))
sys = dss(P)

The result is:

DescriptorStateSpace{Measurement{Float64}, UniformScaling{Bool}, Matrix{Measurement{Float64}}}

State matrix A:
3×3 Matrix{Measurement{Float64}}:
 -800.0±230.0  -5000.0±2500.0  -0.0±0.0
    1.0±0.0        0.0±0.0      0.0±0.0
    0.0±0.0        1.0±0.0      0.0±0.0

Input matrix B:
3×1 Matrix{Measurement{Float64}}:
 1.0 ± 0.0
 0.0 ± 0.0
 0.0 ± 0.0

Output matrix C:
1×3 Matrix{Measurement{Float64}}:
 0.0±0.0  0.0±0.0  5000.0±1700.0

Feedthrough matrix D:
1×1 Matrix{Measurement{Float64}}:
 0.0 ± 0.0

Continuous-time state-space model.

I think this should work also in ControlSystems.jl, but something goes wrong:

using ControlSystems
using Measurements

s = tf("s")

J = 2.0*measurement(1,0.2)
b = 0.04*measurement(1,0.2)
K = 1.0*measurement(1,0.2)
R = 0.08*measurement(1,0.2)
L = 1e-4*measurement(1,0.2)

P = K/(s*((J*s + b)*(L*s + R) + K^2))
sys = ss(P)

julia> P = K/(s*((J*s + b)*(L*s + R) + K^2))
TransferFunction{Continuous, ControlSystemsBase.SisoRational{Measurement{Float64}}}
                    1.0 ± 0.2
-------------------------------------------------
0.0002 ± 5.7e-5s^3 + 0.16 ± 0.045s^2 + 1.0 ± 0.4s

Continuous-time transfer function model

julia> sys = ss(P)
┌ Warning: Unable to balance state-space, returning original system
└ @ ControlSystemsBase C:\Users\Andreas\.julia\packages\ControlSystemsBase\NYWcD\src\types\conversion.jl:166
ERROR: MethodError: no method matching Float64(::Measurement{Float64})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:265
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:900
  Float64(::IrrationalConstants.Log2π)
   @ IrrationalConstants C:\Users\Andreas\.julia\packages\IrrationalConstants\vp5v4\src\macro.jl:112
  ...

Stacktrace:
  [1] convert(::Type{Float64}, x::Measurement{Float64})
    @ Base .\number.jl:7
  [2] setindex!(A::Memory{Float64}, x::Measurement{Float64}, i1::Int64)
    @ Base .\genericmemory.jl:243
  [3] unsafe_copyto!(dest::Memory{Float64}, doffs::Int64, src::Memory{Measurement{Float64}}, soffs::Int64, n::Int64)
    @ Base .\genericmemory.jl:153
  [4] unsafe_copyto!
    @ .\genericmemory.jl:133 [inlined]
  [5] _copyto_impl!
    @ .\array.jl:308 [inlined]
  [6] copyto!
    @ .\array.jl:294 [inlined]
  [7] copyto!
    @ .\array.jl:319 [inlined]
  [8] copyto_axcheck!
    @ .\abstractarray.jl:1167 [inlined]
  [9] Matrix{Float64}(x::Matrix{Measurement{Float64}})
    @ Base .\array.jl:626
 [10] convert
    @ .\array.jl:618 [inlined]
 [11] to_matrix
    @ C:\Users\Andreas\.julia\packages\ControlSystemsBase\NYWcD\src\utilities.jl:18 [inlined]
 [12] fix_D_matrix
    @ C:\Users\Andreas\.julia\packages\ControlSystemsBase\NYWcD\src\types\StateSpace.jl:72 [inlined]
 [13] StateSpace{Continuous, Float64}(A::Matrix{Measurement{Float64}}, B::Matrix{Measurement{Float64}}, C::Matrix{Measurement{Float64}}, D::Matrix{Measurement{Float64}}, timeevol::Continuous)
    @ ControlSystemsBase C:\Users\Andreas\.julia\packages\ControlSystemsBase\NYWcD\src\types\StateSpace.jl:90
 [14] convert(::Type{StateSpace{…}}, G::TransferFunction{Continuous, ControlSystemsBase.SisoRational{…}}; balance::Bool, minimal::Bool, contr::Bool, obs::Bool, noseig::Bool, kwargs::@Kwargs{})
    @ ControlSystemsBase C:\Users\Andreas\.julia\packages\ControlSystemsBase\NYWcD\src\types\conversion.jl:111
 [15] convert
    @ C:\Users\Andreas\.julia\packages\ControlSystemsBase\NYWcD\src\types\conversion.jl:70 [inlined]
 [16] #convert#155
    @ C:\Users\Andreas\.julia\packages\ControlSystemsBase\NYWcD\src\types\conversion.jl:62 [inlined]
 [17] convert
    @ C:\Users\Andreas\.julia\packages\ControlSystemsBase\NYWcD\src\types\conversion.jl:59 [inlined]
 [18] StateSpace
    @ C:\Users\Andreas\.julia\packages\ControlSystemsBase\NYWcD\src\types\StateSpace.jl:126 [inlined]
 [19] ss(args::TransferFunction{Continuous, ControlSystemsBase.SisoRational{Measurement{Float64}}})
    @ ControlSystemsBase C:\Users\Andreas\.julia\packages\ControlSystemsBase\NYWcD\src\types\StateSpace.jl:146
 [20] top-level scope
    @ REPL[62]:1
Some type information was truncated. Use `show(err)` to see complete types.

Probably there is a simple fix for this issue.

[baggepinnen]

Cool! I have a documentation page dedicated to representing uncertainty in systems using uncertain numbers
https://juliacontrol.github.io/RobustAndOptimalControl.jl/dev/uncertainty/#Parametric-uncertainty-using-[MonteCarloMeasurements.jl](https://github.com/baggepinnen/MonteCarloMeasurements.jl)

I'll have a look at what's causing the issue you report!

[baggepinnen]

The issue is that the function one on on a measurement return a different type

julia> x
1.0 ± 0.0

julia> one(x)
1.0

I guess it is technically following the interface outlined in the docstring for one, but 1.0 ± 0.0 would also have been a multiplicative identity and the one ControlSystems expect to make promotion work. I seem to have found a fix for it though