From d7e7d3a0d70307b48dfaa49ca512d76c24a3b3ce Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 22 Jul 2025 07:42:26 +0200
Subject: [PATCH 1/2] handle B matrix with non-full column rank in place

---
 lib/ControlSystemsBase/src/synthesis.jl | 18 ++++++++++++++++--
 1 file changed, 16 insertions(+), 2 deletions(-)

diff --git a/lib/ControlSystemsBase/src/synthesis.jl b/lib/ControlSystemsBase/src/synthesis.jl
index 02d66498f..02beaa8fa 100644
--- a/lib/ControlSystemsBase/src/synthesis.jl
+++ b/lib/ControlSystemsBase/src/synthesis.jl
@@ -242,7 +242,7 @@ 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)
+    QRB = qr(B, ColumnNorm())
     U0, U1 = QRB.Q[:, 1:m], QRB.Q[:, m+1:end] # TODO: check dimension
     Z = QRB.R
     R = svdvals(B)
@@ -259,9 +259,23 @@ function place_knvd(A::AbstractMatrix, B, λ; verbose=false, init=:s, method = 0
         return B\(A - ABF) # Solve for F in (A - BF) = Λ
     end
 
+    mB = size(B, 2)
+    if mB > m
+        # 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, :]'
+        F = place(A, B2, λ; verbose, init, method)
+        return pinv(T)'*F
+    end
+
     S = Matrix{CT}[]
     for j = 1:n
-        qj = qr((U1'*(A- λ[j]*I))')
+        H = (U1'*(A- λ[j]*I))
+        qj = qr(H')
         # Ŝj = qj.Q[:, 1:n-m] # Needed for method 2
         Sj = qj.Q[:, n-m+1:n]
         push!(S, Sj)

From 8a1b3c303412c563c0c35ae3672e39c6e7f6ced8 Mon Sep 17 00:00:00 2001
From: Fredrik Bagge Carlson <baggepinnen@gmail.com>
Date: Tue, 22 Jul 2025 08:08:35 +0200
Subject: [PATCH 2/2] add additional test case

---
 lib/ControlSystemsBase/src/synthesis.jl       | 10 ++++----
 lib/ControlSystemsBase/test/test_synthesis.jl | 24 ++++++++++++++++++-
 2 files changed, 27 insertions(+), 7 deletions(-)

diff --git a/lib/ControlSystemsBase/src/synthesis.jl b/lib/ControlSystemsBase/src/synthesis.jl
index 02beaa8fa..b2451d2ea 100644
--- a/lib/ControlSystemsBase/src/synthesis.jl
+++ b/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
diff --git a/lib/ControlSystemsBase/test/test_synthesis.jl b/lib/ControlSystemsBase/test/test_synthesis.jl
index 239cc6527..f6b22f64f 100644
--- a/lib/ControlSystemsBase/test/test_synthesis.jl
+++ b/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