Riccati Matrix Equation Solvers
Standard Riccati Matrix Equations
MatrixEquations.arec
— Functionarec(A, G, Q = 0; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, Z, scalinfo)
Compute X
, the hermitian/symmetric stabilizing solution (if as = false
) or anti-stabilizing solution (if as = true
) of the continuous-time algebraic Riccati equation
A'X + XA - XGX + Q = 0,
where G
and Q
are hermitian/symmetric matrices or uniform scaling operators. Scalar-valued G
and Q
are interpreted as appropriately sized uniform scaling operators G*I
and Q*I
. The Schur method of [1] is used.
To enhance the accuracy of computations, a block scaling of matrices G
and Q
is performed, if the default setting scaling = 'B'
is used. This scaling is however performed only if norm(Q) > norm(G)
. A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G'
is selected. An alternative, experimental structure preserving scaling can be performed using the option scaling = 'S'
. A symmetric matrix equilibration based scaling is employed if scaling = 'K'
, for which the underlying vector norm can be specified using the keyword argument nrm = p
, where p = 1
is the default setting. Scaling can be disabled with the choice scaling = 'N'
. If pow2 = true
, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false
).
By default, the lower bound for the 1-norm reciprocal condition number rtol
is n*ϵ
, where n
is the order of A
and ϵ
is the machine epsilon of the element type of A
.
EVALS
is a vector containing the (stable or anti-stable) eigenvalues of A-GX
.
Z = [U; V]
is an orthogonal basis for the stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi
, where Sx
and Sxi
are diagonal scaling matrices contained in the named tuple scalinfo
as scalinfo.Sx
and scalinfo.Sxi
, respectively.
Note: To solve the continuous-time algebraic Riccati equation
A'X + XA - XBR^(-1)B'X + Q = 0,
with R
a hermitian/symmetric matrix and B
a compatible size matrix, G = BR^(-1)B'
must be provided. This approach is not numerically suited when R
is ill-conditioned and/or B
has large norm.
Reference:
[1] Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.
Example
julia> using LinearAlgebra
julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
-6.0 -2.0 1.0
5.0 1.0 -1.0
-4.0 -2.0 -1.0
julia> G = [1. 0. 0.; 0. 5. 0.; 0. 0. 10.]
3×3 Array{Float64,2}:
1.0 0.0 0.0
0.0 5.0 0.0
0.0 0.0 10.0
julia> X, CLSEIG = arec(A,G,2I);
julia> X
3×3 Array{Float64,2}:
0.459589 0.333603 -0.144406
0.333603 0.65916 -0.0999216
-0.144406 -0.0999216 0.340483
julia> A'*X+X*A-X*G*X+2I
3×3 Array{Float64,2}:
2.22045e-16 4.44089e-16 -1.77636e-15
4.44089e-16 6.66134e-16 1.11022e-16
-1.77636e-15 1.11022e-16 -1.33227e-15
julia> CLSEIG
3-element Array{Complex{Float64},1}:
-4.411547592296008 + 2.4222082620381102im
-4.411547592296008 - 2.4222082620381102im
-4.337128244724371 + 0.0im
julia> eigvals(A-G*X)
3-element Array{Complex{Float64},1}:
-4.4115475922960075 - 2.4222082620381076im
-4.4115475922960075 + 2.4222082620381076im
-4.337128244724374 + 0.0im
arec(A, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, orth = false, nrm = 1) -> (X, EVALS, F, Z, scalinfo)
Compute X
, the hermitian/symmetric stabilizing solution (if as = false
) or anti-stabilizing solution (if as = true
) of the continuous-time algebraic Riccati equation
A'X + XA - (XB+S)R^(-1)(B'X+S') + Q = 0,
where R
and Q
are hermitian/symmetric matrices or uniform scaling operators such that R
is nonsingular. Scalar-valued R
and Q
are interpreted as appropriately sized uniform scaling operators R*I
and Q*I
. S
, if not specified, is set to S = zeros(size(B))
. The Schur method of [1] is used.
To enhance the accuracy of computations, a block scaling of matrices R
, Q
and S
is performed, if the default setting scaling = 'B'
is used. This scaling is however performed only if norm(Q) > norm(B)^2/norm(R)
. A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G'
is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'
. A symmetric matrix equilibration based scaling is employed if scaling = 'K'
, for which the underlying vector norm can be specified using the keyword argument nrm = p
, where p = 1
is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D'
or scaling = 'T'
. Scaling can be disabled with the choice scaling = 'N'
. If pow2 = true
, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false
).
By default, the lower bound for the 1-norm reciprocal condition number rtol
is n*ϵ
, where n
is the order of A
and ϵ
is the machine epsilon of the element type of A
.
EVALS
is a vector containing the (stable or anti-stable) eigenvalues of A-BF
.
F
is the stabilizing or anti-stabilizing gain matrix F = R^(-1)(B'X+S')
.
Z = [U; V; W]
is a basis for the relevant stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi
and F = -Sr*(W/U)*Sxi
, where Sx
, Sxi
and Sr
are diagonal scaling matrices contained in the named tuple scalinfo
as scalinfo.Sx
, scalinfo.Sxi
and scalinfo.Sr
, respectively. An orthogonal basis Z
can be determined, with an increased computational cost, by setting orth = true
.
Reference:
[1] Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.
Example
julia> using LinearAlgebra
julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
-6.0 -2.0 1.0
5.0 1.0 -1.0
-4.0 -2.0 -1.0
julia> B = [1. 2.; 2. 0.; 0. 1.]
3×2 Array{Float64,2}:
1.0 2.0
2.0 0.0
0.0 1.0
julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
1.0 0.0
0.0 5.0
julia> X, CLSEIG, F = arec(A,B,R,2I);
julia> X
3×3 Array{Float64,2}:
0.522588 0.303007 -0.327227
0.303007 0.650895 -0.132608
-0.327227 -0.132608 0.629825
julia> A'*X+X*A-X*B*inv(R)*B'*X+2I
3×3 Array{Float64,2}:
-2.66454e-15 -1.55431e-15 8.88178e-16
-1.55431e-15 2.22045e-15 -2.9976e-15
9.99201e-16 -2.9976e-15 4.44089e-16
julia> CLSEIG
3-element Array{Complex{Float64},1}:
-4.37703628399912 + 2.8107164873731247im
-4.37703628399912 - 2.8107164873731247im
-1.8663764577096091 + 0.0im
julia> eigvals(A-B*F)
3-element Array{Complex{Float64},1}:
-4.377036283999118 - 2.8107164873731234im
-4.377036283999118 + 2.8107164873731234im
-1.8663764577096063 + 0.0im
arec(A, B, G, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, orth = false, nrm = 1) -> (X, EVALS, F, Z, scalinfo)
Computes X
, the hermitian/symmetric stabilizing solution (if as = false
) or anti-stabilizing solution (if as = true
) of the continuous-time algebraic Riccati equation
A'X + XA - XGX - (XB+S)R^(-1)(B'X+S') + Q = 0,
where G
, R
and Q
are hermitian/symmetric matrices or uniform scaling operators such that R
is nonsingular. Scalar-valued G
, R
and Q
are interpreted as appropriately sized uniform scaling operators G*I
, R*I
and Q*I
. S
, if not specified, is set to S = zeros(size(B))
. For well conditioned R
, the Schur method of [1] is used. For ill-conditioned R
or if orth = true
, the generalized Schur method of [2] is used.
To enhance the accuracy of computations, a block oriented scaling of matrices G
, R
, Q
and S
is performed using the default setting scaling = 'B'
. This scaling is performed only if norm(Q) > max(norm(G), norm(B)^2/norm(R))
. A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G'
is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'
. A symmetric matrix equilibration based scaling is employed if scaling = 'K'
, for which the underlying vector norm can be specified using the keyword argument nrm = p
, where p = 1
is the default setting. If orth = true
, two experimental scaling procedures can be activated using the options scaling = 'D'
and scaling = 'T'
. Scaling can be disabled with the choice scaling = 'N'
.
By default, the lower bound for the 1-norm reciprocal condition number rtol
is n*ϵ
, where n
is the order of A
and ϵ
is the machine epsilon of the element type of A
.
EVALS
is a vector containing the (stable or anti-stable) eigenvalues of A-BF-GX
.
F
is the stabilizing or anti-stabilizing gain matrix F = R^(-1)(B'X+S')
.
Z = [U; V; W]
is a basis for the relevant stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi
and F = -Sr*(W/U)*Sxi
, where Sx
, Sxi
and Sr
are diagonal scaling matrices contained in the named tuple scalinfo
as scalinfo.Sx
, scalinfo.Sxi
and scalinfo.Sr
, respectively. An orthogonal basis Z
can be determined, with an increased computational cost, by setting orth = true
.
Reference:
[1] Laub, A.J., A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.
[2] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.
MatrixEquations.ared
— Functionared(A, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)
Compute X
, the hermitian/symmetric stabilizing solution (if as = false
) or anti-stabilizing solution (if as = true
) of the discrete-time algebraic Riccati equation
A'XA - X - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,
where R
and Q
are hermitian/symmetric matrices. Scalar-valued R
and Q
are interpreted as appropriately sized uniform scaling operators R*I
and Q*I
. S
, if not specified, is set to S = zeros(size(B))
.
To enhance the accuracy of computations, a block oriented scaling of matrices R,
Q
and S
is performed using the default setting scaling = 'B'
. This scaling is performed only if norm(Q) > norm(B)^2/norm(R)
. A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G'
is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'
. A symmetric matrix equilibration based scaling is employed if scaling = 'K'
, for which the underlying vector norm can be specified using the keyword argument nrm = p
, where p = 1
is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D'
, scaling = 'R'
and scaling = 'T'
. Scaling can be disabled with the choice scaling = 'N'
. If pow2 = true
, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false
).
By default, the lower bound for the 1-norm reciprocal condition number rtol
is n*ϵ
, where n
is the order of A
and ϵ
is the machine epsilon of the element type of A
.
EVALS
is a vector containing the (stable) eigenvalues of A-BF
.
F
is the stabilizing gain matrix F = (R+B'XB)^(-1)(B'XA+S')
.
Z = [U; V; W]
is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi
and F = -Sr*(W/U)*Sxi
, where Sx
, Sxi
and Sr
are diagonal scaling matrices contained in the named tuple scalinfo
as scalinfo.Sx
, scalinfo.Sxi
and scalinfo.Sr
, respectively.
Reference:
[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.
Example
julia> using LinearAlgebra
julia> A = [ 0. 1.; 0. 0. ]
2×2 Array{Float64,2}:
0.0 1.0
0.0 0.0
julia> B = [ 0.; sqrt(2.) ]
2-element Array{Float64,1}:
0.0
1.4142135623730951
julia> R = 1.
1.0
julia> Q = [ 1. -1.; -1. 1. ]
2×2 Array{Float64,2}:
1.0 -1.0
-1.0 1.0
julia> X, CLSEIG, F = ared(A,B,R,Q);
julia> X
2×2 Array{Float64,2}:
1.0 -1.0
-1.0 1.5
julia> A'*X*A-X-A'*X*B*inv(R+B'*X*B)*B'*X*A+Q
2×2 Array{Float64,2}:
0.0 -3.33067e-16
-3.33067e-16 8.88178e-16
julia> CLSEIG
2-element Array{Complex{Float64},1}:
0.4999999999999998 - 0.0im
-0.0 - 0.0im
julia> eigvals(A-B*F)
2-element Array{Float64,1}:
-2.7755575615628914e-16
0.5
Generalized Riccati Matrix Equations
MatrixEquations.garec
— Functiongarec(A, E, G, Q = 0; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, Z, scalinfo)
Compute X
, the hermitian/symmetric stabilizing solution (if as = false
) or anti-stabilizing solution (if as = true
) of the generalized continuous-time algebraic Riccati equation
A'XE + E'XA - E'XGXE + Q = 0,
where G
and Q
are hermitian/symmetric matrices or uniform scaling operators and E
is a nonsingular matrix. Scalar-valued G
and Q
are interpreted as appropriately sized uniform scaling operators G*I
and Q*I
. The generalized Schur method of [1] is used.
To enhance the accuracy of computations, a block scaling of matrices G
and Q
is performed, if the default setting scaling = 'B'
is used. This scaling is however performed only if norm(Q) > norm(G)
. A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G'
is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'
. A symmetric matrix equilibration based scaling is employed if scaling = 'K'
, for which the underlying vector norm can be specified using the keyword argument nrm = p
, where p = 1
is the default setting. Scaling can be disabled with the choice scaling = 'N'
. If pow2 = true
, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false
).
By default, the lower bound for the 1-norm reciprocal condition number rtol
is n*ϵ
, where n
is the order of A
and ϵ
is the machine epsilon of the element type of A
.
EVALS
is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-GXE,E)
.
Z = [U; V]
is an orthogonal basis for the stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E
, where Sx
and Sxi
are diagonal scaling matrices contained in the named tuple scalinfo
as scalinfo.Sx
and scalinfo.Sxi
, respectively.
Note: To solve the continuous-time algebraic Riccati equation
A'XE + E'XA - E'XBR^(-1)B'XE + Q = 0,
with R
a hermitian/symmetric matrix and B
a compatible size matrix, G = BR^(-1)B'
must be provided. This approach is not numerically suited when R
is ill-conditioned and/or B
has large norm.
Reference:
[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.
garec(A, E, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)
Compute X
, the hermitian/symmetric stabilizing solution (if as = false
) or anti-stabilizing solution (if as = true
) of the generalized continuous-time algebraic Riccati equation
A'XE + E'XA - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,
where R
and Q
are hermitian/symmetric matrices such that R
is nonsingular, and E
is a nonsingular matrix. Scalar-valued R
and Q
are interpreted as appropriately sized uniform scaling operators R*I
and Q*I
. S
, if not specified, is set to S = zeros(size(B))
. The generalized Schur method of [1] is used.
To enhance the accuracy of computations, a block oriented scaling of matrices R,
Q
and S
is performed using the default setting scaling = 'B'
. This scaling is performed only if norm(Q) > norm(B)^2/norm(R)
. A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G'
is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'
. A symmetric matrix equilibration based scaling is employed if scaling = 'K'
, for which the underlying vector norm can be specified using the keyword argument nrm = p
, where p = 1
is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D'
or scaling = 'T'
. Scaling can be disabled with the choice scaling = 'N'
. If pow2 = true
, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false
).
By default, the lower bound for the 1-norm reciprocal condition number rtol
is n*ϵ
, where n
is the order of A
and ϵ
is the machine epsilon of the element type of A
.
EVALS
is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF,E)
.
F
is the stabilizing/anti-stabilizing gain matrix F = R^(-1)(B'XE+S')
.
Z = [U; V; W]
is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E
and F = -Sr*(W/U)*Sxi
, where Sx
, Sxi
and Sr
are diagonal scaling matrices contained in the named tuple scalinfo
as scalinfo.Sx
, scalinfo.Sxi
and scalinfo.Sr
, respectively.
Reference:
[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.
Example
julia> using LinearAlgebra
julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
-6.0 -2.0 1.0
5.0 1.0 -1.0
-4.0 -2.0 -1.0
julia> E = [10. 3. 0.; 0. 5. -1.; 0. 0. 10.]
3×3 Array{Float64,2}:
10.0 3.0 0.0
0.0 5.0 -1.0
0.0 0.0 10.0
julia> B = [1. 2.; 2. 0.; 0. 1.]
3×2 Array{Float64,2}:
1.0 2.0
2.0 0.0
0.0 1.0
julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
1.0 0.0
0.0 5.0
julia> X, CLSEIG, F = garec(A,E,B,R,2I);
julia> X
3×3 Array{Float64,2}:
0.0502214 0.0284089 -0.0303703
0.0284089 0.111219 -0.00259162
-0.0303703 -0.00259162 0.0618395
julia> A'*X*E+E'*X*A-E'*X*B*inv(R)*B'*X*E+2I
3×3 Array{Float64,2}:
1.55431e-15 -1.9984e-15 -3.33067e-15
-1.77636e-15 1.33227e-15 -3.33067e-15
-2.88658e-15 -3.21965e-15 1.11022e-15
julia> CLSEIG
3-element Array{Complex{Float64},1}:
-0.6184265391601464 + 0.2913286844595737im
-0.6184265391601464 - 0.2913286844595737im
-0.21613059964451786 + 0.0im
julia> eigvals(A-B*F,E)
3-element Array{Complex{Float64},1}:
-0.6184265391601462 - 0.29132868445957383im
-0.6184265391601462 + 0.2913286844595739im
-0.216130599644518 + 0.0im
garec(A, E, B, G, R, Q, S; scaling = 'B', pw2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)
Compute X
, the hermitian/symmetric stabilizing solution (if as = false
) or anti-stabilizing solution (if as = true
) of the generalized continuous-time algebraic Riccati equation
A'XE + E'XA - E'XGXE - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,
where G
, Q
and R
are hermitian/symmetric matrices such that R
is nonsingular, and E
is a nonsingular matrix. Scalar-valued G
, R
and Q
are interpreted as appropriately sized uniform scaling operators G*I
, R*I
and Q*I
. The generalized Schur method of [1] is used.
To enhance the accuracy of computations, a block oriented scaling of matrices G,
R,
Q
and S
is performed using the default setting scaling = 'B'
. This scaling is performed only if norm(Q) > max(norm(G), norm(B)^2/norm(R))
. A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G'
is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'
. A symmetric matrix equilibration based scaling is employed if scaling = 'K'
, for which the underlying vector norm can be specified using the keyword argument nrm = p
, where p = 1
is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D'
or scaling = 'T'
. Scaling can be disabled with the choice scaling = 'N'
. If pow2 = true
, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false
).
By default, the lower bound for the 1-norm reciprocal condition number rtol
is n*ϵ
, where n
is the order of A
and ϵ
is the machine epsilon of the element type of A
.
EVALS
is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF-GXE,E)
.
F
is the stabilizing/anti-stabilizing gain matrix F = R^(-1)(B'XE+S')
.
Z = [U; V; W]
is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E
and F = -Sr*(W/U)*Sxi
, where Sx
, Sxi
and Sr
are diagonal scaling matrices contained in the named tuple scalinfo
as scalinfo.Sx
, scalinfo.Sxi
and scalinfo.Sr
, respectively.
Reference:
[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.
MatrixEquations.gared
— Functiongared(A, E, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)
Compute X
, the hermitian/symmetric stabilizing solution (if as = false
) or anti-stabilizing solution (if as = true
) of the generalized discrete-time algebraic Riccati equation
A'XA - E'XE - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,
where R
and Q
are hermitian/symmetric matrices, and E
ist non-singular. Scalar-valued R
and Q
are interpreted as appropriately sized uniform scaling operators R*I
and Q*I
. S
, if not specified, is set to S = zeros(size(B))
.
To enhance the accuracy of computations, a block oriented scaling of matrices R,
Q
and S
is performed using the default setting scaling = 'B'
. This scaling is performed only if norm(Q) > norm(B)^2/norm(R)
. A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G'
is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'
. A symmetric matrix equilibration based scaling is employed if scaling = 'K'
, for which the underlying vector norm can be specified using the keyword argument nrm = p
, where p = 1
is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D'
, scaling = 'R'
and scaling = 'T'
. Scaling can be disabled with the choice scaling = 'N'
. If pow2 = true
, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false
).
By default, the lower bound for the 1-norm reciprocal condition number rtol
is n*ϵ
, where n
is the order of A
and ϵ
is the machine epsilon of the element type of A
.
EVALS
is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF,E)
.
F
is the stabilizing/anti-stabilizing gain matrix F = (R+B'XB)^(-1)(B'XA+S')
.
Z = [U; V; W]
is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E
and F = -Sr*(W/U)*Sxi
, where Sx
, Sxi
and Sr
are diagonal scaling matrices contained in the named tuple scalinfo
as scalinfo.Sx
, scalinfo.Sxi
and scalinfo.Sr
, respectively.
Reference:
[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.
Example
julia> using LinearAlgebra
julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
-6.0 -2.0 1.0
5.0 1.0 -1.0
-4.0 -2.0 -1.0
julia> E = [10. 3. 0.; 0. 5. -1.; 0. 0. 10.]
3×3 Array{Float64,2}:
10.0 3.0 0.0
0.0 5.0 -1.0
0.0 0.0 10.0
julia> B = [1. 2.; 2. 0.; 0. 1.]
3×2 Array{Float64,2}:
1.0 2.0
2.0 0.0
0.0 1.0
julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
1.0 0.0
0.0 5.0
julia> X, CLSEIG, F = gared(A,E,B,R,2I);
julia> X
3×3 Array{Float64,2}:
0.065865 -0.0147205 -0.0100407
-0.0147205 0.0885939 0.0101422
-0.0100407 0.0101422 0.0234425
julia> A'*X*A-E'*X*E-A'*X*B*inv(R+B'*X*B)*B'*X*A+2I
3×3 Array{Float64,2}:
-1.33227e-15 -2.48412e-15 1.38778e-16
-2.498e-15 -4.44089e-16 -6.50521e-16
1.80411e-16 -5.91541e-16 -1.33227e-15
julia> CLSEIG
3-element Array{Complex{Float64},1}:
-0.084235615751339 - 0.0im
-0.190533552034239 - 0.0im
-0.5238922629921539 - 0.0im
julia> eigvals(A-B*F,E)
3-element Array{Float64,1}:
-0.5238922629921539
-0.19053355203423886
-0.08423561575133902