andreasvarga / MatrixEquations.jl

Solution of Lyapunov, Sylvester and Riccati matrix equations using Julia
MIT License
79 stars 8 forks source link

Spread the use of `Discrete` and `Continuous` throughout the package? #17

Open dkarrasch opened 2 years ago

dkarrasch commented 2 years ago

There is an opportunity to make use of the above mentioned types also in other places than the MEoperators. For instance, one could make

lyapd(args...) = lyap(Discrete(), args...)
lyapc(args...) = lyap(Continuous(), args...)

and get away from having a new function for each job, but have one function with two methods. I'm not sure that's helpful/possible, but one potential benefit could be to write ContinuousOrDiscrete-agnostic "metacodes" for both cases that just do the right thing depending on the Discrete or Continuous type. This is issue is meant to discuss if at all and in case yes, then which parts could benefit from this little refactoring. I think this could be made nonbreaking by the above redirection.

andreasvarga commented 2 years ago

This is certainly possible, but it would involve some effort which in this moment I would like to invest to start a new project. This idea would be probably the best way to enhance the definitions of the basic objects in the DescriptorSystems package. I have to mention that embeding the system type into the definition of control related objects is already implemented in the ControlSystems package (in a somehow more involved variant where discrete(Ts) is used instead discrete(), with Ts the sampling time).

Probably you remember that there is already a function lyap in Julia as part of the LinearAlgebra package. The corresponding function in ME is lyapc, which, in contrast to lyap produces hermitian/symmetric solution if the right-hand side is hermitian/symmetric, and, in principle, is the right way to compute the solution. If the right-hand side is non-hermitian/non-symmetric, the two implementations are practically the same (relying on the same LAPACK wrapper). Recently I improved significantly the performance of Lyapunov solvers, which are now (almost) at the level of the performance of Matlab solvers (based on Fortran codes). See also the discussion here. Further enhancements are certainly possible, but as I already mentioned, would involve some additional efforts. It would however be nice to have a lyap function in Julia which covers both continuous and discrete cases. The way you propose to implement it would be very nice.

baggepinnen commented 2 years ago

Expanding upon https://github.com/andreasvarga/DescriptorSystems.jl/issues/8#issuecomment-1046181390

As Andreas says, in ControlSystems.jl, we have gone the route of including the time-domain in the type information, i.e., Continuous/Discrete is part of the type, but the numeric value of the same time is not. It's important to not include the sample time value itself as a type parameter, since then you would incur compilation of all methods for each unique sample time, the sample time is data stored as a field of type T in the type Discrete{T}. The only time our approach is inadequate is in situations when you'd like to promote the I in [I sysd], since you do not know which sample time to give the static system I due to the sample time not being part of the type. One solution to this problem would be to allow the sample time -1 to denote "inherited" sample time, i.e., that the systems that have sample time -1 should adopt the sample times of the surrounding systems. I have yet to implement the "inherited" sample time in ControlSystems.jl, but for all other purposes, our model has worked out well.

We have recently also gone the route of dispatching matrix equations based on these types like @dkarrasch mentions in the OP, see e.g. https://github.com/JuliaControl/ControlSystems.jl/blob/master/src/matrix_comps.jl#L9

We dispatch based on the following union types

const ContinuousType = Union{Continuous, Type{Continuous}}

in order to be able to pass either a type or an instance. This allows us to avoid manual dispatch like

if isdiscrete(sys)
    lyapd(...)
else
    lyapc(...)
end

and instead call

lyap(sys.timeevol, ...)

Edit: We now have a complete implementation also of the machinery for promoting unknown sample times. See the function common_timeevol which we call in the beginning of each function that takes more than one system definition in order to get a consistent TimeEvolution instance and throw an error if there is a mismatch: https://github.com/JuliaControl/ControlSystems.jl/blob/08eee7ecceddcc5b5bd68d309f7f4d8b7df88352/src/types/TimeEvolution.jl#L42