WIP julia1.0 compatability

Author: mfalt

LICENSE.md

@@ -1,6 +1,7 @@
 ControlSystems.jl is licensed under the MIT License:
 
-> Copyright (c) 2014: Jim Crist and other contributors:
+> Copyright (c) 2014-2018: Jim Crist, Mattias Fält, Fredrik Bagge Carlson
+> and other contributors:
 >
 > Permission is hereby granted, free of charge, to any person obtaining
 > a copy of this software and associated documentation files (the

README.md

@@ -46,7 +46,7 @@ Some of the available commands are:
 ##### Constructing systems
 ss, tf, zpk, ss2tf
 ##### Analysis
-pole, tzero, norm, norminf, ctrb, obsv, gangoffour, margin, markovparam, damp, dampreport, zpkdata, dcgain, covar, gram, sigma, sisomargin
+pole, tzero, sysnorm, norminf, ctrb, obsv, gangoffour, margin, markovparam, damp, dampreport, zpkdata, dcgain, covar, gram, sigma, sisomargin
 ##### Synthesis
 care, dare, dlyap, lqr, dlqr, place, pid, leadlink, laglink, leadlinkat, rstd, rstc, dab
 ##### Time and Frequency response

REQUIRE

@@ -4,3 +4,6 @@ Polynomials
 LaTeXStrings
 OrdinaryDiffEq
 IterTools
+LinearAlgebra
+Printf
+Colors

docs/src/examples/example.md

@@ -7,13 +7,14 @@ end
 
 # LQR design
 ```julia
+using LinearAlgebra # For identity matrix I
 h       = 0.1
 A       = [1 h; 0 1]
 B       = [0;1]
 C       = [1 0]
 sys     = ss(A,B,C,0, h)
-Q       = eye(2)
-R       = eye(1)
+Q       = I
+R       = I
 L       = dlqr(A,B,Q,R) # lqr(sys,Q,R) can also be used
 
 u(t,x)  = -L*x + 1.5(t>=2.5)# Form control law (u is a function of t and x), a constant input disturbance is affecting the system from t≧2.5
@@ -28,13 +29,14 @@ plot(t,x, lab=["Position", "Velocity"]', xlabel="Time [s]")
 
 # LQR design
 ```julia
+using LinearAlgebra # For identity matrix I
 h       = 0.1
 A       = [1 h; 0 1]
 B       = [0;1]
 C       = [1 0]
 sys     = ss(A,B,C,0, h)
-Q       = eye(2)
-R       = eye(1)
+Q       = I
+R       = I
 L       = dlqr(A,B,Q,R) # lqr(sys,Q,R) can also be used
 
 u(t,x)  = -L*x + 1.5(t>=2.5)# Form control law (u is a function of t and x), a constant input disturbance is affecting the system from t≧2.5

docs/src/lib/analysis.md

@@ -16,7 +16,7 @@ gangofseven
 gram
 margin
 markovparam
-norm
+sysnorm
 obsv
 pole
 sigma

docs/src/makeplots.jl

@@ -12,8 +12,8 @@ A = [1 h; 0 1]
 B = [0;1]
 C = [1 0]
 sys = ss(A,B,C,0, h)
-Q = eye(2)
-R = eye(1)
+Q = Matrix{Float64}(I,2,2n)
+R = Matrix{Float64}(I,1,1)
 L = dlqr(A,B,Q,R) # lqr(sys,Q,R) can also be used
 
 u(x,t)  = -L*x + 1.5(t>=2.5)# Form control law (u is a function of t and x), a constant input disturbance is affecting the system from t≧2.5

src/ControlSystems.jl

@@ -66,14 +66,19 @@ export  LTISystem,
         nyquist,
         sigma,
         # utilities
+        num,    #Deprecated
+        den,    #Deprecated
+        numvec,
+        denvec,
         numpoly,
         denpoly
 
 
 # QUESTION: are these used? LaTeXStrings, Requires, IterTools
-using Polynomials, OrdinaryDiffEq, Plots, LaTeXStrings
-import Base: +, -, *, /, (./), (==), (.+), (.-), (.*), (!=), isapprox, convert, promote_op
-import Base.LinAlg: BlasFloat
+using Polynomials, OrdinaryDiffEq, Plots, LaTeXStrings, LinearAlgebra
+import Base: +, -, *, /, (==), (!=), isapprox, convert, promote_op
+import LinearAlgebra: BlasFloat
+import Printf, Colors
 
 abstract type AbstractSystem end
 
@@ -120,7 +125,10 @@ include("pid_design.jl")
 
 include("plotting.jl")
 
+@deprecate num numvec
+@deprecate den denvec
 
+# @deprecate norm sysnorm
 
 # The path has to be evaluated upon initial import
 const __CONTROLSYSTEMS_SOURCE_DIR__ = dirname(Base.source_path())

src/analysis.jl

@@ -72,15 +72,15 @@ Display a report of the poles, damping ratio, natural frequency, and time
 constant of the system `sys`"""
 function dampreport(io::IO, sys::LTISystem)
     Wn, zeta, ps = damp(sys)
-    t_const = 1./(Wn.*zeta)
+    t_const = 1 ./ (Wn.*zeta)
     header =
     ("|     Pole      |   Damping     |   Frequency   | Time Constant |\n"*
     "|               |    Ratio      |   (rad/sec)   |     (sec)     |\n"*
     "+---------------+---------------+---------------+---------------+")
     println(io, header)
     for i=eachindex(ps)
         p, z, w, t = ps[i], zeta[i], Wn[i], t_const[i]
-        @printf(io, "|  %-13.3e|  %-13.3e|  %-13.3e|  %-13.3e|\n", p, z, w, t)
+        Printf.@printf(io, "|  %-13.3e|  %-13.3e|  %-13.3e|  %-13.3e|\n", p, z, w, t)
     end
 end
 dampreport(sys::LTISystem) = dampreport(STDOUT, sys)
@@ -131,7 +131,7 @@ function tzero(A::Matrix{T}, B::Matrix{T}, C::Matrix{T}, D::Matrix{T}) where T<:
         nf = size(A, 1)
         m = size(D, 2)
         Af = ([A B] * W)[1:nf, 1:nf]
-        Bf = ([eye(nf) zeros(nf, m)] * W)[1:nf, 1:nf]
+        Bf = ([I zeros(nf, m)] * W)[1:nf, 1:nf]
         zs = eigvals(Af, Bf)
         _fix_conjugate_pairs!(zs) # Generalized eigvals does not return exact conj. pairs
     else
@@ -180,10 +180,10 @@ function reduce_sys(A::Matrix{T}, B::Matrix{T}, C::Matrix{T}, D::Matrix{T}, meps
         end
         # Update System
         n, m = size(B)
-        Vm = [V zeros(T, n, m); zeros(T, m, n) eye(T, m)]
+        Vm = [V zeros(T, n, m); zeros(T, m, n) I]
         if sigma > 0
             M = [A B; Cbar Dbar]
-            Vs = [V' zeros(T, n, sigma) ; zeros(T, sigma, n) eye(T, sigma)]
+            Vs = [V' zeros(T, n, sigma) ; zeros(T, sigma, n) I]
         else
             M = [A B]
             Vs = V'
@@ -223,15 +223,19 @@ If `full` return also `fullPhase`
 """
 function margin(sys::LTISystem, w::AbstractVector{S}; full=false, allMargins=false) where S<:Real
     ny, nu = size(sys)
-    vals = (:wgm, :gm, :wpm, :pm, :fullPhase)
+
     if allMargins
-        for val in vals
-            eval(:($val = Array{Array{eltype(sys),1}}($ny,$nu)))
-        end
+        wgm         = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
+        gm          = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
+        wpm         = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
+        pm          = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
+        fullPhase   = Array{Array{numeric_type(sys),1}}(undef, ny,nu)
     else
-        for val in vals
-            eval(:($val = Array{eltype(sys),2}($ny,$nu)))
-        end
+        wgm         = Array{numeric_type(sys),2}(undef, ny, nu)
+        gm          = Array{numeric_type(sys),2}(undef, ny, nu)
+        wpm         = Array{numeric_type(sys),2}(undef, ny, nu)
+        pm          = Array{numeric_type(sys),2}(undef, ny, nu)
+        fullPhase   = Array{numeric_type(sys),2}(undef, ny, nu)
     end
     for j=1:nu
         for i=1:ny
@@ -259,7 +263,7 @@ function sisomargin(sys::LTISystem, w::AbstractVector{S}; full=false, allMargins
     wgm, = _allPhaseCrossings(w, phase)
     gm = similar(wgm)
     for i = eachindex(wgm)
-        gm[i] = 1./abs(freqresp(sys,[wgm[i]])[1][1])
+        gm[i] = 1 ./ abs(freqresp(sys,[wgm[i]])[1][1])
     end
     wpm, fi = _allGainCrossings(w, mag)
     pm = similar(wpm)

src/connections.jl

@@ -25,10 +25,10 @@ function append(systems::StateSpace...)
     if !all(s.Ts == Ts for s in systems)
         error("Sampling time mismatch")
     end
-    A = blkdiag([s.A for s in systems]...)
-    B = blkdiag([s.B for s in systems]...)
-    C = blkdiag([s.C for s in systems]...)
-    D = blkdiag([s.D for s in systems]...)
+    A = blockdiag([s.A for s in systems]...)
+    B = blockdiag([s.B for s in systems]...)
+    C = blockdiag([s.C for s in systems]...)
+    D = blockdiag([s.D for s in systems]...)
     return StateSpace(A, B, C, D, Ts)
 end
 
@@ -37,7 +37,7 @@ function append(systems::TransferFunction...)
     if !all(s.Ts == Ts for s in systems)
         error("Sampling time mismatch")
     end
-    mat = blkdiag([s.matrix for s in systems]...)
+    mat = blockdiag([s.matrix for s in systems]...)
     return TransferFunction(mat, Ts)
 end
 
@@ -53,9 +53,9 @@ function Base.vcat(systems::StateSpace...)
     if !all(s.Ts == Ts for s in systems)
         error("Sampling time mismatch")
     end
-    A = blkdiag([s.A for s in systems]...)
+    A = blockdiag([s.A for s in systems]...)
     B = vcat([s.B for s in systems]...)
-    C = blkdiag([s.C for s in systems]...)
+    C = blockdiag([s.C for s in systems]...)
     D = vcat([s.D for s in systems]...)
 
     return StateSpace(A, B, C, D, Ts)
@@ -87,8 +87,8 @@ function Base.hcat(systems::StateSpace...)
     if !all(s.Ts == Ts for s in systems)
         error("Sampling time mismatch")
     end
-    A = blkdiag([s.A for s in systems]...)
-    B = blkdiag([s.B for s in systems]...)
+    A = blockdiag([s.A for s in systems]...)
+    B = blockdiag([s.B for s in systems]...)
     C = hcat([s.C for s in systems]...)
     D = hcat([s.D for s in systems]...)
 
@@ -156,16 +156,12 @@ end
 #     end
 # end
 
-# TODO: could use cat([1,2], mats...) instead
-# Empty definition to get rid of warning
-Base.blkdiag() = []
-function Base.blkdiag(mats::Matrix...)
+
+blockdiag(mats::Matrix...) = blockdiag(promote(mats...)...)
+
+function blockdiag(mats::Matrix{T}...) where T
     rows = Int[size(m, 1) for m in mats]
     cols = Int[size(m, 2) for m in mats]
-    T = eltype(mats[1])
-    for ind=1:length(mats)
-        T = promote_type(T, eltype(mats[ind]))
-    end
     res = zeros(T, sum(rows), sum(cols))
     m = 1
     n = 1

src/discrete.jl

@@ -17,24 +17,24 @@ function c2d(sys::StateSpace, Ts::Real, method::Symbol=:zoh)
     ny, nu = size(sys)
     nx = nstates(sys)
     if method == :zoh
-        M = expm([A*Ts  B*Ts;
+        M = exp([A*Ts  B*Ts;
             zeros(nu, nx + nu)])
         Ad = M[1:nx, 1:nx]
         Bd = M[1:nx, nx+1:nx+nu]
         Cd = C
         Dd = D
-        x0map = [eye(nx) zeros(nx, nu)]
+        x0map = [I zeros(nx, nu)]
     elseif method == :foh
-        M = expm([A*Ts B*Ts zeros(nx, nu);
-            zeros(nu, nx + nu) eye(nu);
+        M = exp([A*Ts B*Ts zeros(nx, nu);
+            zeros(nu, nx + nu) I;
             zeros(nu, nx + 2*nu)])
         M1 = M[1:nx, nx+1:nx+nu]
         M2 = M[1:nx, nx+nu+1:nx+2*nu]
         Ad = M[1:nx, 1:nx]
         Bd = Ad*M2 + M1 - M2
         Cd = C
         Dd = D + C*M2
-        x0map = [eye(nx) -M2]
+        x0map = [I -M2]
     elseif method == :tustin || method == :matched
         error("NotImplemented: Only `:zoh` and `:foh` implemented so far")
     else
@@ -230,9 +230,9 @@ function lsima(sys::StateSpace, t::AbstractVector, r::AbstractVector{T}, control
         dsys, x0map = c2d(sys, dt, :foh)
     end
     n = size(t, 1)
-    x = Array{T}(size(sys.A, 1), n)
-    u = Array{T}(n)
-    y = Array{T}(n)
+    x = Array{T}(undef, size(sys.A, 1), n)
+    u = Array{T}(undef, n)
+    y = Array{T}(undef, n)
     for i=1:n
         x[:,i] = x0
         y[i] = (sys.C*x0 + sys.D*u[i])[1]

src/freqresp.jl

@@ -15,7 +15,7 @@ function freqresp(sys::LTISystem, w_vec::AbstractVector{S}) where {S<:Real}
     end
     Tsys = numeric_type(sys)
     T = promote_type(typeof(zero(Tsys)/norm(one(Tsys))), Complex64, S)
-    sys_fr = Array{T}(length(w_vec), noutputs(sys), ninputs(sys))
+    sys_fr = Array{T}(undef, length(w_vec), noutputs(sys), ninputs(sys))
 
     if isa(sys, StateSpace)
         sys = _preprocess_for_freqresp(sys)
@@ -76,7 +76,7 @@ end
 
 function evalfr(G::TransferFunction{<:SisoTf{T0}}, s::Number) where {T0}
     T = promote_type(T0, typeof(one(T0)*one(typeof(s))/(one(T0)*one(typeof(s)))))
-    fr = Array{T}(size(G))
+    fr = Array{T}(undef, size(G))
     for j = 1:ninputs(G)
         for i = 1:noutputs(G)
             fr[i, j] = evalfr(G.matrix[i, j], s)
@@ -110,7 +110,7 @@ function (sys::TransferFunction)(z_or_omegas::AbstractVector, map_to_unit_circle
     @assert !iscontinuous(sys) "It makes no sense to call this function with continuous systems"
     vals = sys.(z_or_omegas, map_to_unit_circle)# evalfr.(sys,exp.(evalpoints))
     # Reshape from vector of evalfr matrizes, to (in,out,freq) Array
-    out = Array{eltype(eltype(vals)),3}(length(z_or_omegas), size(sys)...)
+    out = Array{eltype(eltype(vals)),3}(undef, length(z_or_omegas), size(sys)...)
     for i in 1:length(vals)
         out[i,:,:] .= vals[i]
     end
@@ -150,7 +150,7 @@ frequencies `w`
 function sigma(sys::LTISystem, w::AbstractVector)
     resp = freqresp(sys, w)
     nw, ny, nu = size(resp)
-    sv = Array{Float64}(nw, min(ny, nu))
+    sv = Array{Float64}(undef, nw, min(ny, nu))
     for i=1:nw
         sv[i, :] = svdvals(resp[i, :, :])
     end
@@ -162,7 +162,7 @@ function _default_freq_vector(systems::Vector{T}, plot::Symbol) where T<:LTISyst
     min_pt_per_dec = 60
     min_pt_total = 200
     nsys = length(systems)
-    bounds = Array{Float64}(2, nsys)
+    bounds = Array{Float64}(undef, 2, nsys)
     for i=1:nsys
         # TODO : For now we ignore the feature information. In the future,
         # these can be used to improve the frequency vector near features.