andreasvarga / matrixequations.jl Goto Github PK

Solution of Lyapunov, Sylvester and Riccati matrix equations using Julia

License: MIT License

Julia 100.00%

riccati-equations lyapunov-equations matrix-equations sylvester-operators control-systems lyapunov-operators condition-estimation sylvester-equations

matrixequations.jl's People

Contributors

Stargazers

Watchers

Forkers

mortenpi mateuszbaran juliatagbot tangwang-ustc putianyi889 baggepinnen blegat erikqqy

matrixequations.jl's Issues

Routine Package Test Failure

Hello,

After installing and executing the tests in package mode, I received the following error message:

I also would like to mention that my particular use case does not depend on this test, but I was not sure if this error was encountered by others. Can someone else reproduce this error? For what it's worth, my current version of Julia is 1.9. Thank you in advance for reviewing my comment.

possible test failure in upcoming Julia version 1.5

A PkgEval run for a PR which changes the generated numbers for randn! indicates that the tests of this package might fail in Julia 1.5 (and on Julia current master). Apologies if this is a false positive.
cf. https://github.com/JuliaCI/NanosoldierReports/blob/7de24e455342298cbef56826b5827f0d7640d2c1/pkgeval/by_hash/b89e35c_vs_098ef24/logs/MatrixEquations/1.5.0-DEV-71a4a114c2.log

Condition Matrix(T)' = Matrix(T') not fulfilled for all defined operators

Let T be, for example, a Lyapunov operator T:X -> Y, which maps the upper triangular parts of symmetric matrices X to the upper triangular parts of the symmetric matrices Y computed as Y := AX+XA'. Similarly, the transposed operator Tt =T' maps X to Y := A'X+XA. Consider the example:

julia> A = [11 12; 21 22]
2×2 Array{Int64,2}:
 11  12
 21  22

julia> T = lyapop(A,her=true)
Linear operator
  nrow: 3
  ncol: 3
  eltype: Int64
  symmetric: false
  hermitian: false

julia> Tt = lyapop(A',her=true)
Linear operator
  nrow: 3
  ncol: 3
  eltype: Int64
  symmetric: false
  hermitian: false

julia> Matrix(T)
3×3 Array{Int64,2}:
 22  24   0
 21  33  12
  0  42  44

julia> Matrix(T')
3×3 Array{Int64,2}:
 22  42   0
 12  33  21
  0  24  44

julia> Matrix(Tt)
3×3 Array{Int64,2}:
 22  42   0
 12  33  21
  0  24  44

julia> Matrix(T)' == Matrix(T')
false

julia> Matrix(Tt)' == Matrix(Tt')
false

However we also have, as expected,

julia> Matrix(T') == Matrix(Tt)
true

julia> Matrix(Tt') == Matrix(T)
true

Fixing this inconsistent behaviour would be very desirable in a future release.

Documentation error of `arec`

In the README.md,

arec Solution of the continuous Riccati equations AX+XA'-XRX+Q = 0 and A'X+XA-(XB+S)R^(-1)(B'X+S')+Q = 0.

should be

arec Solution of the continuous Riccati equations A'X+XA-XRX+Q = 0 and A'X+XA-(XB+S)R^(-1)(B'X+S')+Q = 0.

according to the latest docs.

Performance Benchmarks

Hi, it is common for Julia packages to share performance benchmarks in the Github readme.
See ParallelKmeans.jl or FixedEffectModels.jl.
Would this be possible for this package at some point?

Issue with LapackUtil

Hi Andreas,

first, let me express I am very pleased to see that you cannot abandon your lifetime topic and decided to contribute to Julia ecosystem. This might be a major reinforcement for the part of Julia community interested in control-related algorithms.

Second, let me inform you that I am having troubles to even start using your package. This is what I get after adding MatrixEquations.jl to my Julia 1.1:

julia> using MatrixEquations
[ Info: Precompiling MatrixEquations [99c1a7ee-ab34-5fd5-8076-27c950a045f4]
ERROR: LoadError: ArgumentError: Package MatrixEquations does not have LapackUtil in its dependencies:
- If you have MatrixEquations checked out for development and have
  added LapackUtil as a dependency but haven't updated your primary
  environment's manifest file, try `Pkg.resolve()`.
- Otherwise you may need to report an issue with MatrixEquations
Stacktrace:
 [1] require(::Module, ::Symbol) at ./loading.jl:836
 [2] include at ./boot.jl:326 [inlined]
 [3] include_relative(::Module, ::String) at ./loading.jl:1038
 [4] include(::Module, ::String) at ./sysimg.jl:29
 [5] top-level scope at none:2
 [6] eval at ./boot.jl:328 [inlined]
 [7] eval(::Expr) at ./client.jl:404
 [8] top-level scope at ./none:3
in expression starting at /home/hurak/.julia/packages/MatrixEquations/T51DM/src/MatrixEquations.jl:4
ERROR: Failed to precompile MatrixEquations [99c1a7ee-ab34-5fd5-8076-27c950a045f4] to /home/hurak/.julia/compiled/v1.1/MatrixEquations/1uOBF.ji.
Stacktrace:
 [1] compilecache(::Base.PkgId, ::String) at ./loading.jl:1197
 [2] _require(::Base.PkgId) at ./loading.jl:960
 [3] require(::Base.PkgId) at ./loading.jl:858
 [4] require(::Module, ::Symbol) at ./loading.jl:853

This must be somehow related to the package design and implementation. I am not yet sufficiently fluent with these things to advise.

Economics Example

Hi and thank you for this package.

Solving optimal control problems is at the heart of economics.
(including LQR/LQG etc)
I'm working to see what packages in the Julia ecosystem can be helpful.
See my note on discourse.

I want to replicate that (simple) firm investment problem with your package.

Or consider a slightly more general economics HW problem from Ken Kasa:
Continuous Time

Discrete Time

PS: here are my closed form solutions

Note in the Readme/Docs

It may help to put something like the following in the Readme/Docs

This can help many users bc there are different formulations of LQR with slightly different notation.
For example, I lost hours today because in my formulation the objective had x'Su instead of 2x'Su...

Move Alternating-direction Implicit (ADI) Method solver here?

We have a quick implementation of ADI:

https://github.com/dlfivefifty/AlternatingDirectionImplicit.jl

This could live in this package.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Algebraic equation in matrix form

Hi, I wanted to solve a Algebraic equation in matrix form that looks like:

The equation has the form: AX=0
So I need to solve det(A)=0 with respect to k. Actually it is a Algebraic equation in matrix form, is it possible to solve it with MatrixEquations.jl ?
Note: It is easy to do it by running:
syms k; solution = double(solve(det(A)));
Any suggestion will be helpful, Thank you very much!

usage of MatrixEquations.jl in Manifolds.jl

Thank you for a this useful package! In Manifolds.jl we recently started using MatrixEquations.jl for projection on the Stiefiel manifold: JuliaManifolds/Manifolds.jl#49 . Since we are supporting Julia 1.0 I'm considering whether we should conditionally use your package through Requires.jl or maybe you could support Julia 1.0? All tests in MatrixEquations.jl pass for me on Julia 1.0 if I define isnothing (which is really simple: https://github.com/JuliaLang/julia/blob/bc82c359d4fe08750915f5db84fa83bbef063629/base/some.jl#L64 ).

Derivative of cholesky factorization

Given a Cholesky factorization Q = U'U and a symmetric matrix dQ, I would like to find the matrix dU such that dQ = U' * dU + dU' * U. The equation is similar to the ones of this package but does not seem to match any, or does it ?

Orthogonal transformation between commuting matrices

In SumOfSquares, we have pairs of commuting matrices A, B and need to compute an orthogonal transformation X such that X * A = B * X. This is implemented in orthogonal_transformation_to. It basically calls schur(A) and schur(B) but then it's a bit subtle since the schur decomposition is not unique so we need to canonicalize in a similar way for both. It looks a bit like some of the equations of this package so I was wondering if it could fit into one of them.

Adjoint handling

From what I understand, the adj field is typically computed from type information, like isa(A, Adjoint) and things like that. That means that the value of adj can be computed at compile time. Once the objects are generated, this type information is lost and only a Bool value is kept. This makes it unreachable for the compiler and in particular for multiple dispatch.

I'd propose to consider to move the adj information into a type parameter, for instance like:

struct GeneralizedLyapunovMap{T,TA <: AbstractMatrix{T},TE <: AbstractMatrix{T},Adj} <: LyapunovMatrixEquationsMaps{T}
   A::TA
   E::TE
   disc::Bool
   her::Bool
   adj::Bool
   function GeneralizedLyapunovMap(A::AbstractMatrix{T}, E::AbstractMatrix{T}; disc=false, her=false) where {T}
      n = LinearAlgebra.checksquare(A)
      n == LinearAlgebra.checksquare(E) ||
               throw(DimensionMismatch("E must be a square matrix of dimension $n"))
      return new{T,typeof(A),typeof(E),false}(A, E, disc, her)
   end
   function GeneralizedLyapunovMap(A::Adjoint{T}, E::Adjoint{T}; disc=false, her=false) where {T}
      n = LinearAlgebra.checksquare(A)
      n == LinearAlgebra.checksquare(E) ||
               throw(DimensionMismatch("E must be a square matrix of dimension $n"))
      return new{T,typeof(A),typeof(E),true}(A, E, disc, her)
   end
end

You could then dispatch on types ::GeneralizedLyapunovMap{T,TA,TE,true} and ::GeneralizedLyapunovMap{T,TA,TE,false} if needed. In this specific case, I don't even see why you would need that information, at least not in the meoperators part. Anyway, I don't think it would lead to a major performance increase, but it could help clean up the code a little bit, i.e., have less branches.

`gared` crash on Julia v1.10.0-beta2

Hi,

Thanks for all the great works! It's great to have an OO alternative to MATLAB.

On Julia v.1.10.0-beta2 (you can install it with UpdateJulia.jl or juliaup), ared function crash.

If I do:

using MatrixEquations
A=[0.5;;]
B=[1.0;;]
Q=[1.0;;]
R=[1.0;;]
ared(A, B, R, Q)

I get this error:

ERROR: MethodError: no method matching *(::LinearAlgebra.AdjointQ{Float64, LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}}, ::LinearMaps.BlockMap{Float64, Tuple{LinearMaps.WrappedMap{Float64, Matrix{Float64}}, LinearMaps.UniformScalingMap{Bool}, Vararg{LinearMaps.WrappedMap{Float64, Matrix{Float64}}, 4}}, Tuple{Int64, Int64}})

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:585
  *(::LinearMaps.ScaledMap, ::LinearMaps.LinearMap)
   @ LinearMaps ~/.julia/packages/LinearMaps/gMfcQ/src/scaledmap.jl:49
  *(::LinearMaps.CompositeMap, ::LinearMaps.LinearMap)
   @ LinearMaps ~/.julia/packages/LinearMaps/gMfcQ/src/composition.jl:147
  ...

Stacktrace:
 [1] gared(A::Matrix{…}, E::LinearAlgebra.UniformScaling{…}, B::Matrix{…}, R::Matrix{…}, Q::Matrix{…}, S::Matrix{…}; as::Bool, rtol::Float64)
   @ MatrixEquations ~/.julia/packages/MatrixEquations/vBvfe/src/riccati.jl:833
 [2] gared
   @ MatrixEquations ~/.julia/packages/MatrixEquations/vBvfe/src/riccati.jl:774 [inlined]
 [3] #ared#38
   @ MatrixEquations ~/.julia/packages/MatrixEquations/vBvfe/src/riccati.jl:693 [inlined]
 [4] ared
   @ MatrixEquations ~/.julia/packages/MatrixEquations/vBvfe/src/riccati.jl:690 [inlined]
 [5] ared(A::Matrix{Float64}, B::Matrix{Float64}, R::Matrix{Float64}, Q::Matrix{Float64})
   @ MatrixEquations ~/.julia/packages/MatrixEquations/vBvfe/src/riccati.jl:690
 [6] top-level scope
   @ ~/Documents/Personnel/ModelPredictiveControl.jl/example/Untitled-2.jl:6
Some type information was truncated. Use `show(err)` to see complete types.

and the complete types with show(err) are:

julia> show(err)
1-element ExceptionStack:
LoadError: MethodError: no method matching *(::LinearAlgebra.AdjointQ{Float64, LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}}, ::LinearMaps.BlockMap{Float64, Tuple{LinearMaps.WrappedMap{Float64, Matrix{Float64}}, LinearMaps.UniformScalingMap{Bool}, Vararg{LinearMaps.WrappedMap{Float64, Matrix{Float64}}, 4}}, Tuple{Int64, Int64}})

Closest candidates are:
  *(::Any, ::Any, ::Any, ::Any...)
   @ Base operators.jl:585
  *(::LinearMaps.ScaledMap, ::LinearMaps.LinearMap)
   @ LinearMaps ~/.julia/packages/LinearMaps/gMfcQ/src/scaledmap.jl:49
  *(::LinearMaps.CompositeMap, ::LinearMaps.LinearMap)
   @ LinearMaps ~/.julia/packages/LinearMaps/gMfcQ/src/composition.jl:147
  ...

Stacktrace:
  [1] gared(A::Matrix{Float64}, E::LinearAlgebra.UniformScaling{Bool}, B::Matrix{Float64}, R::Matrix{Float64}, Q::Matrix{Float64}, S::Matrix{Float64}; as::Bool, rtol::Float64)
    @ MatrixEquations ~/.julia/packages/MatrixEquations/vBvfe/src/riccati.jl:833
  [2] gared
    @ MatrixEquations ~/.julia/packages/MatrixEquations/vBvfe/src/riccati.jl:774 [inlined]
  [3] #ared#38
    @ MatrixEquations ~/.julia/packages/MatrixEquations/vBvfe/src/riccati.jl:693 [inlined]
  [4] ared
    @ MatrixEquations ~/.julia/packages/MatrixEquations/vBvfe/src/riccati.jl:690 [inlined]
  [5] ared(A::Matrix{Float64}, B::Matrix{Float64}, R::Matrix{Float64}, Q::Matrix{Float64})
    @ MatrixEquations ~/.julia/packages/MatrixEquations/vBvfe/src/riccati.jl:690
  [6] top-level scope
    @ ~/Documents/Personnel/ModelPredictiveControl.jl/example/Untitled-2.jl:6
  [7] eval
    @ Base ./boot.jl:383 [inlined]
  [8] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:2070
  [9] invokelatest(::Any, ::Any, ::Vararg{Any}; kwargs::@Kwargs{})
    @ Base ./essentials.jl:887
 [10] invokelatest(::Any, ::Any, ::Vararg{Any})
    @ Base ./essentials.jl:884
 [11] inlineeval(m::Module, code::String, code_line::Int64, code_column::Int64, file::String; softscope::Bool)
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.51.1/scripts/packages/VSCodeServer/src/eval.jl:261
 [12] (::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.51.1/scripts/packages/VSCodeServer/src/eval.jl:181
 [13] withpath(f::VSCodeServer.var"#67#72"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams}, path::String)
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.51.1/scripts/packages/VSCodeServer/src/repl.jl:274
 [14] (::VSCodeServer.var"#66#71"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.51.1/scripts/packages/VSCodeServer/src/eval.jl:179
 [15] hideprompt(f::VSCodeServer.var"#66#71"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.51.1/scripts/packages/VSCodeServer/src/repl.jl:38
 [16] (::VSCodeServer.var"#65#70"{Bool, Bool, Bool, Module, String, Int64, Int64, String, VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.51.1/scripts/packages/VSCodeServer/src/eval.jl:150
 [17] with_logstate(f::Function, logstate::Any)
    @ Base.CoreLogging ./logging.jl:515
 [18] with_logger
    @ VSCodeServer ./logging.jl:627 [inlined]
 [19] (::VSCodeServer.var"#64#69"{VSCodeServer.ReplRunCodeRequestParams})()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.51.1/scripts/packages/VSCodeServer/src/eval.jl:253
 [20] #invokelatest#2
    @ Base ./essentials.jl:887 [inlined]
 [21] invokelatest(::Any)
    @ Base ./essentials.jl:884
 [22] (::VSCodeServer.var"#62#63")()
    @ VSCodeServer ~/.vscode/extensions/julialang.language-julia-1.51.1/scripts/packages/VSCodeServer/src/eval.jl:34
in expression starting at /home/francis/Documents/Personnel/ModelPredictiveControl.jl/example/Untitled-2.jl:6

Francis Gagnon

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

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.