Merge pull request #1011 from JuliaControl/realplace

make place return real matrix when expected

Author: baggepinnen

docs/src/examples/example.md

@@ -74,33 +74,33 @@ using LinearAlgebra
 using Plots
 
 # Create system - same as LQR example
-Ts      = 0.1
-A       = [1 Ts; 0 1]
-B       = [0; 1]
-C       = I(2)
-P       = ss(A,B,C,0,Ts)
+Ts = 0.1
+A  = [1 Ts; 0 1]
+B  = [0; 1]
+C  = I(2)
+P  = ss(A,B,C,0,Ts)
 
 # Design controller using pole placement
 # Choose desired closed-loop poles (well-damped, faster than original system)
-desired_poles_cont = [-2+0.5im, -2-0.5im]  # Continuous-time poles
-desired_poles = exp.(Ts .* desired_poles_cont)  # Discrete-time poles inside unit circle
+desired_poles_cont = [-2+0.5im, -2-0.5im]       # Continuous-time poles
+desired_poles = exp.(Ts .* desired_poles_cont)  # Discrete-time poles
 
 # Design state feedback gain using pole placement
-L = place(P, desired_poles) |> real
+L = place(P, desired_poles)
 
-# Design observer with poles 5x faster (closer to origin for discrete time)
+# Design observer with poles 5x faster
 observer_poles = exp.(Ts*5 .* desired_poles_cont)
-K = place(P, observer_poles, :o) |> real # Note the :o for observer design
+K = place(P, observer_poles, :o) # Note the :o for observer design
 
-# Create observer-controller using the observer_controller function
+# Create controller system
 controller = observer_controller(P, L, K)
 
 # Form closed-loop system and analyze
 T_cl = feedback(P * controller)
 
 r(x,t)  = [1.5(t>=2.5); 0] # Form control law (r is a function of t and x), change reference to 1.5 at t≧2.5
-t = 0:Ts:5              # Time vector
-x0 = [1, 0, 0, 0]               # Initial condition
+t = 0:Ts:5          # Time vector
+x0 = [1.0, 0, 0, 0] # Initial condition (plant state followed by controller state)
 res = lsim(T_cl, r, t; x0)
 plot(res, lab=["Position" "Velocity"], layout=1, sp=1)
 ```

lib/ControlSystemsBase/src/synthesis.jl

@@ -150,7 +150,7 @@ Please note that this function can be numerically sensitive, solving the placeme
 function place(A, B, p, opt=:c; direct = false, kwargs...)
     n = length(p)
     n != size(A,1) && error("Must specify as many poles as the state dimension")
-    if opt === :c
+    L = if opt === :c
         direct && error("direct = true only applies to observer design")
         n != size(B,1) && error("A and B must have same number of rows")
         if size(B,2) == 1
@@ -172,6 +172,11 @@ function place(A, B, p, opt=:c; direct = false, kwargs...)
     else
         error("fourth argument must be :c or :o")
     end
+    if isreal(A) && isreal(B) && is_self_conjugate(p)
+        @assert all(abs(imag(l)) .< 1e-6 for l in L) "Expected real coefficient in feedback gain, got complex: $L"
+        return real(L)
+    end
+    return L
 end
 function place(sys::AbstractStateSpace, p, opt=:c; direct = false, kwargs...)
     if opt === :c

lib/ControlSystemsBase/src/types/SisoTfTypes/SisoZpk.jl

@@ -129,6 +129,16 @@ function pairup_conjugates!(x::AbstractVector)
     end
     return true
 end
+function is_self_conjugate(x::AbstractVector)
+    length(x) == 0 && return true # Empty vector is self-conjugate
+    # Check that all complex numbers have a conjugate in the vector
+    for i = eachindex(x)
+        if findfirst(≈(conj(x[i])), x) === nothing
+            return false
+        end
+    end
+    return true
+end
 
 function evalfr(f::SisoZpk{T1,TR}, s::T2) where {T1<:Number, TR<:Number, T2<:Number}
     T0 = promote_type(T2, TR)