Merge pull request #236 from JuliaControl/release-0.5.4-option2

Release 0.5.4

Author: mfalt

.travis.yml

@@ -2,6 +2,7 @@ language: julia
 julia:
   - 1.0
   - 1.1
+  - 1.2
   - nightly
 matrix:
   allow_failures:

Project.toml

@@ -2,7 +2,7 @@ name = "ControlSystems"
 uuid = "a6e380b2-a6ca-5380-bf3e-84a91bcd477e"
 authors = ["Dept. Automatic Control, Lund University"]
 repo = "https://github.com/JuliaControl/ControlSystems.jl.git"
-version = "0.5.3"
+version = "0.5.4"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
@@ -20,4 +20,4 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-Polynomials = "0.5.3"
+Polynomials = "0.6.0"

README.md

@@ -17,6 +17,9 @@ Pkg.add("ControlSystems")
 ```
 
 ## News
+### 2019-11-03
+- Poles and zeros are "not sorted" as in Julia versions < 1.2, even on newer versions of Julia. This should imply that complex conjugates are kept together.
+
 ### 2019-05-28
 #### Delay systems
 - We now support systems with time delays. Example:

src/ControlSystems.jl

@@ -81,7 +81,9 @@ export  LTISystem,
 
 
 # QUESTION: are these used? LaTeXStrings, Requires, IterTools
-using Polynomials, Plots, LaTeXStrings, LinearAlgebra
+using Plots, LaTeXStrings, LinearAlgebra
+import Polynomials
+import Polynomials: Poly, coeffs, polyval
 using OrdinaryDiffEq, DelayDiffEq
 export Plots
 import Base: +, -, *, /, (==), (!=), isapprox, convert, promote_op
@@ -99,7 +101,6 @@ include("types/SisoTf.jl")
 
 # Transfer functions and tranfer function elemements
 include("types/TransferFunction.jl")
-include("types/SisoTfTypes/polyprint.jl")
 include("types/SisoTfTypes/SisoZpk.jl")
 include("types/SisoTfTypes/SisoRational.jl")
 include("types/SisoTfTypes/promotion.jl")

src/analysis.jl

@@ -1,7 +1,7 @@
 """`pole(sys)`
 
 Compute the poles of system `sys`."""
-pole(sys::AbstractStateSpace) = eigvals(sys.A)
+pole(sys::AbstractStateSpace) = eigvalsnosort(sys.A)
 pole(sys::SisoTf) = error("pole is not implemented for type $(typeof(sys))")
 
 # Seems to have a lot of rounding problems if we run the full thing with sisorational,
@@ -224,7 +224,7 @@ function tzero(A::AbstractMatrix{T}, B::AbstractMatrix{T}, C::AbstractMatrix{T},
         m = size(D_rc, 2)
         Af = ([A_rc B_rc] * W)[1:nf, 1:nf]
         Bf = ([Matrix{T}(I, nf, nf) zeros(nf, m)] * W)[1:nf, 1:nf]
-        zs = eigvals(Af, Bf)
+        zs = eigvalsnosort(Af, Bf)
         _fix_conjugate_pairs!(zs) # Generalized eigvals does not return exact conj. pairs
     else
         zs = complex(T)[]

src/matrix_comps.jl

@@ -47,9 +47,9 @@ function dare(A, B, Q, R)
     if (!isposdef(R))
         error("R must be positive definite.");
     end
-    
+
     n = size(A, 1);
-    
+
     E = [
         Matrix{Float64}(I, n, n) B/R*B';
         zeros(size(A)) A'
@@ -58,10 +58,10 @@ function dare(A, B, Q, R)
         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
 
@@ -292,7 +292,7 @@ function normLinf_twoSteps_ct(sys::AbstractStateSpace, tol=1e-6, maxIters=1000,
         if isreal(p)  # only real poles
             omegap = minimum(abs.(p))
         else  # at least one pair of complex poles
-            tmp = maximum(abs.(imag.(p)./(real.(p).*abs.(p))))
+            tmp = abs.(imag.(p)./(real.(p).*abs.(p)))
             omegap = abs(p[argmax(tmp)])    # TODO This is highly suspicious
         end
         (lb, idx) = findmax([lb, T(maximum(svdvals(evalfr(sys, omegap*1im))))]) #TODO remove T() in julia 0.7.0

src/types/SisoTfTypes/polyprint.jl

@@ -285,12 +285,12 @@ end
 # TODO: Could perhaps be made more accurate. See: An accurate and efficient
 # algorithm for the computation of the # characteristic polynomial of a general square matrix.
 function charpoly(A::AbstractMatrix{<:Number})
-    Λ = eigvals(A)
+    Λ = eigvalsnosort(A)
 
     return prod(roots2poly_factors(Λ)) # Compute the polynomial factors directly?
 end
 function charpoly(A::AbstractMatrix{<:Real})
-    Λ = eigvals(A)
+    Λ = eigvalsnosort(A)
     return prod(roots2real_poly_factors(Λ))
 end
 

src/types/conversion.jl

@@ -16,35 +16,38 @@ to_matrix(T, A::Number) = fill(T(A), 1, 1)
 # Handle Adjoint Matrices
 to_matrix(T, A::Adjoint{R, MT}) where {R<:Number, MT<:AbstractMatrix} = to_matrix(T, MT(A))
 
+# Do no sorting of eigenvalues
+@static if VERSION > v"1.2.0-DEV.0"
+    eigvalsnosort(args...; kwargs...) = eigvals(args...; sortby=nothing, kwargs...)
+    roots(args...; kwargs...) = Polynomials.roots(args...; sortby=nothing, kwargs...)
+else
+    eigvalsnosort(args...; kwargs...) = eigvals(args...; kwargs...)
+    roots(args...; kwargs...) = Polynomials.roots(args...; kwargs...)
+end
 
+""" f = printpolyfun(var)
+`fun` Prints polynomial in descending order, with variable `var`
+"""
+printpolyfun(var) = (io, p, mimetype = MIME"text/plain"()) -> Polynomials.printpoly(io, p, mimetype, descending_powers=true, var=var)
 
 # NOTE: Tolerances for checking real-ness removed, shouldn't happen from LAPACK?
 # TODO: This doesn't play too well with dual numbers..
 # Allocate for maxiumum possible length of polynomial vector?
 #
 # This function rely on that the every complex roots is followed by its exact conjugate,
-# and that the first complex root in each pair has positive real part. This formaat is always
+# and that the first complex root in each pair has positive imaginary part. This format is always
 # returned by LAPACK routines for eigenvalues.
 function roots2real_poly_factors(roots::Vector{cT}) where cT <: Number
     T = real(cT)
     poly_factors = Vector{Poly{T}}()
-    @static if VERSION > v"1.2.0-DEV.0" # Sort one more time to handle GenericLinearAlgebra not being updated
-        sort!(roots, by=LinearAlgebra.eigsortby)
-    end
     for k=1:length(roots)
         r = roots[k]
 
         if isreal(r)
             push!(poly_factors,Poly{T}([-real(r),1]))
         else
-            @static if VERSION > v"1.2.0-DEV.0" # Flipped order in this version
-                if imag(r) > 0 # This roots was handled in the previous iteration # TODO: Fix better error handling
-                    continue
-                end
-            else
-                if imag(r) < 0 # This roots was handled in the previous iteration # TODO: Fix better error handling
-                    continue
-                end
+            if imag(r) < 0 # This roots was handled in the previous iteration # TODO: Fix better error handling
+                continue
             end
 
             if k == length(roots) || r != conj(roots[k+1])

src/utilities.jl

@@ -63,11 +63,11 @@ A = randn(3,3)
 B = randn(3,1)
 p = [3.0,2,1]
 K = ControlSystems.acker(A,B,p)
-@test eigvals(A-B*K) ≈ p
+@test ControlSystems.eigvalsnosort(A-B*K) ≈ p
 
 p = [-1+im, -1-im, -1]
 K = ControlSystems.acker(A,B,p)
-@test eigvals(A-B*K) ≈ p
+@test ControlSystems.eigvalsnosort(A-B*K) ≈ p
 end
 
 end