add additional test case

Author: baggepinnen

lib/ControlSystemsBase/src/synthesis.jl

@@ -242,11 +242,11 @@ function place_knvd(A::AbstractMatrix, B, λ; verbose=false, init=:s, method = 0
     λ = sort(vec(λ), by=LinearAlgebra.eigsortby)
     length(λ) == size(A, 1) == n || error("Must specify as many poles as the state dimension")
     Λ = diagm(λ)
-    QRB = qr(B, ColumnNorm())
-    U0, U1 = QRB.Q[:, 1:m], QRB.Q[:, m+1:end] # TODO: check dimension
-    Z = QRB.R
     R = svdvals(B)
     m = count(>(100*eps()*R[1]), R) # Rank of B
+    QRB = qr(B, ColumnNorm())
+    U0, U1 = QRB.Q[:, 1:m], QRB.Q[:, m+1:end] # TODO: check dimension
+    Z = (QRB.R*QRB.P')[:, 1:m] 
     if m == n # Easy case, B is full rank
         r = count(e->imag(e) == 0, λ)
         ABF = diagm(real(λ))
@@ -264,10 +264,8 @@ function place_knvd(A::AbstractMatrix, B, λ; verbose=false, init=:s, method = 0
         # several inputs but not full column rank, this case must be handled separately
         # when B does not have full column rank but that rank is not 1. In that case, find B2 and T from rank-revealing QR (qr(B, ColumnNorm())
         verbose && @info "Projecting down to rank of B"
-        # T = ones(mB)
-        # B2 = B*T
         B2 = QRB.Q[:, 1:m]
-        T = Z[1:m, :]'
+        T = QRB.P * QRB.R[1:m, :]'
         F = place(A, B2, λ; verbose, init, method)
         return pinv(T)'*F
     end

lib/ControlSystemsBase/test/test_synthesis.jl

@@ -114,7 +114,7 @@ end
         # @show cond(gram(sys, :c))
         (; A, B) = sys
         p = [-1.0, -2, -3]
-        L = place(A, B, p)
+        L = place(A, B, p, verbose=false)
         @test eltype(L) <: Real
         @test allin(eigvals(A - B*L), p)
 
@@ -177,6 +177,28 @@ end
     @test allin(eigvals(A - B*L), p; tol=0.025) # Tolerance from paper
     # norm(L)
 
+    ## Rank deficient multi-input B
+    P = let
+        tempA = [0.0 1.0; -4.0 -1.2]
+        tempB = [0.0 0.0; 4.0 4.0]
+        tempC = [1.0 0.0]
+        tempD = [0.0 0.0]
+        ss(tempA, tempB, tempC, tempD)
+    end
+    F = place(P, [-2, -2], verbose=true)
+    @test allin(eigvals(P.A - P.B*F), [-2, -2])
+
+
+    P = let
+        tempA = [0.0 1.0; -4.0 -1.2]
+        tempB = randn(2, 1) * randn(1, 10) # Rank deficient
+        tempC = [1.0 0.0]
+        tempD = zeros(1, 10)
+        ss(tempA, tempB, tempC, tempD)
+    end
+    F = place(P, [-2, -2], verbose=true)
+    @test allin(eigvals(P.A - P.B*F), [-2, -2])
+
 end
 
 @testset "LQR" begin