Merge 025938e into e821752

src/matrix_comps.jl

@@ -34,37 +34,35 @@ end
 """`dare(A, B, Q, R)`
 
 Compute `X`, the solution to the discrete-time algebraic Riccati equation,
-defined as A'XA - X - (A'XB)(B'XB + R)^-1(B'XA) + Q = 0, where A and R
-are non-singular.
+defined as A'XA - X - (A'XB)(B'XB + R)^-1(B'XA) + Q = 0, where Q>=0 and R>0
 
 Algorithm taken from:
 Laub, "A Schur Method for Solving Algebraic Riccati Equations."
 http://dspace.mit.edu/bitstream/handle/1721.1/1301/R-0859-05666488.pdf
 """
 function dare(A, B, Q, R)
-    G = try
-        B*inv(R)*B'
-    catch
-        error("R must be non-singular.")
+    if (!ishermitian(Q) || minimum(eigvals(real(Q))) < 0)
+        error("Q must be positive-semidefinite.");
     end
-
-    Ait = try
-        inv(A)'
-    catch
-        error("A must be non-singular.")
+    if (!isposdef(R))
+        error("R must be positive definite.");
     end
-
-    Z = [A + G*Ait*Q   -G*Ait;
-         -Ait*Q        Ait]
-
-    S = schur(Z)
-    S = ordschur(S, abs.(S.values).<=1)
-    U = S.Z
-
-    (m, n) = size(U)
-    U11 = U[1:div(m, 2), 1:div(n,2)]
-    U21 = U[div(m,2)+1:m, 1:div(n,2)]
-    return U21/U11
+
+    n = size(A, 1);
+
+    E = [
+        Matrix{Float64}(I, n, n) B/R*B';
+        zeros(size(A)) A'
+    ];
+    F = [
+        A zeros(size(A));
+        -Q Matrix{Float64}(I, n, n)
+    ];
+
+    QZ = schur(F, E);
+    QZ = ordschur(QZ, abs.(QZ.alpha./QZ.beta) .< 1);
+
+    return QZ.Z[(n+1):end, 1:n]/QZ.Z[1:n, 1:n];
 end
 
 """`dlyap(A, Q)`

test/test_linalg.jl

@@ -7,6 +7,33 @@ r = 3
 @test care(a, b, c'*c, r) ≈ [0.5895174372762604 1.8215747248860816; 1.8215747248860823 8.818839806923107]
 @test dare(a, b, c'*c, r) ≈ [240.73344504496302 -131.09928700089387; -131.0992870008943 75.93413176750603]
 
+## Test dare method for non invertible matrices
+A = [0 1; 0 0];
+B0 = [0; 1];
+Q = Matrix{Float64}(I, 2, 2);
+R0 = 1
+X = dare(A, B0, Q, R0);
+# Reshape for matrix expression
+B = reshape(B0, 2, 1)
+R = fill(R0, 1, 1)
+@test norm(A'X*A - X - (A'X*B)*((B'X*B + R)\(B'X*A)) + Q) < 1e-14
+## Test dare for scalars
+A = 1.0;
+B = 1.0;
+Q = 1.0;
+R = 1.0;
+X0 = dare(A, B, Q, R);
+X = X0[1]
+@test norm(A'X*A - X - (A'X*B)*((B'X*B + R)\(B'X*A)) + Q) < 1e-14
+# Test for complex matrices
+I2 = Matrix{Float64}(I, 2, 2)
+@test dare([1.0 im; im 1.0], I2, I2, I2) ≈ [1 + sqrt(2) 0; 0 1 + sqrt(2)]
+# And complex scalars
+@test dare(1.0, 1, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1)
+@test dare(1.0im, 1, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1)
+@test dare(1.0, 1im, 1, 1) ≈ fill((1 + sqrt(5))/2, 1, 1)
+
+## Test gram, ctrb, obsv
 a_2 = [-5 -3; 2 -9]
 C_212 = ss(a_2, [1; 2], [1 0; 0 1], [0; 0])
 C_222 = ss(a_2, [1 0; 0 2], [1 0; 0 1], zeros(2,2))