Merge pull request #94 from JuliaControl/BugsAnd0.6

Julia 0.6

Author: mfalt

.travis.yml

@@ -1,12 +1,12 @@
 language: julia
 julia:
-  - 0.5
+  - 0.6
 notifications:
     email: false
 script:
     - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi
     - julia -e 'Pkg.clone(pwd()); Pkg.build("ControlSystems")'
-    - julia -e 'ENV["PYTHON"] = ""; Pkg.clone("PyPlot"); Pkg.build("PyPlot")'
+    #- julia -e 'ENV["PYTHON"] = ""; Pkg.clone("PyPlot"); Pkg.build("PyPlot")'
     - julia -e 'Pkg.clone("https://github.com/JuliaControl/ControlExamplePlots.jl.git");'
     - julia -e 'Pkg.test("ControlSystems"; coverage=true)'
 after_success:

REQUIRE

@@ -1,3 +1,3 @@
-julia 0.5-
-Plots v0.7.4
+julia 0.6-
+Plots 0.7.4
 Polynomials

docs/src/makeplots.jl

@@ -4,7 +4,7 @@ println("Generating plots")
 
 plotsDir = (pwd()[end-3:end] == "docs") ? "build/plots" : "docs/build/plots"
 mkdir(plotsDir)
-Plots.pyplot()
+Plots.gr()
 
 
 # PID design functions
@@ -76,7 +76,7 @@ pidplots(P,:gof,;kps=kp,kis=ki, ω= logspace(-2,2,500))
 Plots.savefig(plotsDir*"/pidplotgof1.svg")
 
 kp = linspace(-1,1,8) # Now try a different strategy, where we have specified a gain crossover frequency of 0.1 rad/s
-ki = sqrt(1-kp.^2)/10
+ki = sqrt.(1-kp.^2)./10
 pidplots(P,:nyquist,;kps=kp,kis=ki)
 Plots.savefig(plotsDir*"/pidplotsnyquist2.svg")
 pidplots(P,:gof,;kps=kp,kis=ki)

src/ControlSystems.jl

@@ -64,7 +64,7 @@ export  LTISystem,
         denpoly
 
 import Plots
-import Base: +, -, *, /, (./), (==), (.+), (.-), (.*), (!=), isapprox, call, convert, promote_op
+import Base: +, -, *, /, (./), (==), (.+), (.-), (.*), (!=), isapprox, convert, promote_op
 
 include("types/lti.jl")
 include("types/transferfunction.jl")

src/analysis.jl

@@ -52,9 +52,9 @@ Compute the zeros, poles, and gains of system `sys`.
 `k` : Matrix{Float64}, (ny x nu)""" ->
 function zpkdata(sys::LTISystem)
     ny, nu = size(sys)
-    zs = Array(Vector{Complex128}, ny, nu)
-    ps = Array(Vector{Complex128}, ny, nu)
-    ks = Array(Float64, ny, nu)
+    zs = Array{Vector{Complex128}}(ny, nu)
+    ps = Array{Vector{Complex128}}(ny, nu)
+    ks = Array{Float64}(ny, nu)
     for j = 1:nu
         for i = 1:ny
             zs[i, j], ps[i, j], ks[i, j] = _zpk_kern(sys, i, j)
@@ -96,11 +96,11 @@ function damp(sys::LTISystem)
         Ts = sys.Ts == -1 ? 1 : sys.Ts
         ps = log(ps)/Ts
     end
-    Wn = abs(ps)
+    Wn = abs.(ps)
     order = sortperm(Wn)
     Wn = Wn[order]
     ps = ps[order]
-    ζ = -cos(angle(ps))
+    ζ = -cos.(angle.(ps))
     return Wn, ζ, ps
 end
 
@@ -205,7 +205,7 @@ function reduce_sys(A::Matrix{Float64}, B::Matrix{Float64}, C::Matrix{Float64},
             break
         elseif nu == 0
             # System has no zeros, return empty matrices
-            A = B = Cbar = Dbar = Float64[]
+            A = B = Cbar = Dbar = Array{Float64,2}(0,0)
             break
         end
         # Update System
@@ -233,7 +233,7 @@ end
 function fastrank(A::Matrix{Float64}, meps::Float64)
     n, m = size(A, 1, 2)
     if n*m == 0     return 0    end
-    norms = Array(Float64, n)
+    norms = Array{Float64}(n)
     for i = 1:n
         norms[i] = norm(A[i, :])
     end
@@ -300,7 +300,7 @@ function sisomargin{S<:Real}(sys::LTISystem, w::AbstractVector{S}; full=false, a
         gm, idx = findmin([gm;Inf])
         wgm = [wgm;NaN][idx]
         fi = [fi;NaN][idx]
-        pm, idx = findmin([abs(pm);Inf])
+        pm, idx = findmin([abs.(pm);Inf])
         wpm = [wpm;NaN][idx]
         if full
             if !isnan(fi) #fi may be NaN, fullPhase is a scalar

src/connections.jl

@@ -51,7 +51,7 @@ append(systems::LTISystem...) = append(promote(systems...)...)
 #This is needed until julia removes deprecated vect()
 function Base.vect(sys::Union{LTISystem, Real}...)
     T = Base.promote_typeof(sys...)
-    copy!(Array(T,length(sys)), sys)
+    copy!(Array{T}(length(sys)), sys)
 end
 
 function Base.vcat(systems::StateSpace...)

src/discrete.jl

@@ -231,9 +231,9 @@ function lsima{T}(sys::StateSpace, t::AbstractVector, r::AbstractVector{T}, cont
         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}(size(sys.A, 1), n)
+    u = Array{T}(n)
+    y = Array{T}(n)
     for i=1:n
         x[:,i] = x0
         y[i] = (sys.C*x0 + sys.D*u[i])[1]

src/freqresp.jl

@@ -17,7 +17,7 @@ function freqresp{S<:Real}(sys::LTISystem, w::AbstractVector{S})
     end
     #Evil but nessesary type instability here
     sys = _preprocess_for_freqresp(sys)
-    resp = Array(Complex128, nw, ny, nu)
+    resp = Array{Complex128}(nw, ny, nu)
     for i=1:nw
         # TODO : This doesn't actually take advantage of Hessenberg structure
         # for statespace version.
@@ -66,7 +66,7 @@ function evalfr(sys::TransferFunction, s::Number)
     S = promote_type(typeof(s), Float64)
     mat = sys.matrix
     ny, nu = size(mat)
-    res = Array(S, ny, nu)
+    res = Array{S}(ny, nu)
     for j = 1:nu
         for i = 1:ny
             res[i, j] = evalfr(mat[i, j], s)
@@ -111,7 +111,7 @@ at frequencies `w`
 `mag` and `phase` has size `(length(w), ny, nu)`""" ->
 function bode(sys::LTISystem, w::AbstractVector)
     resp = freqresp(sys, w)[1]
-    return abs(resp), rad2deg(unwrap!(angle(resp),1)), w
+    return abs.(resp), rad2deg.(unwrap!(angle.(resp),1)), w
 end
 bode(sys::LTISystem) = bode(sys, _default_freq_vector(sys, :bode))
 
@@ -136,7 +136,7 @@ frequencies `w`
 function sigma(sys::LTISystem, w::AbstractVector)
     resp = freqresp(sys, w)[1]
     nw, ny, nu = size(resp)
-    sv = Array(Float64, nw, min(ny, nu))
+    sv = Array{Float64}(nw, min(ny, nu))
     for i=1:nw
         sv[i, :] = svdvals(resp[i, :, :])
     end
@@ -148,7 +148,7 @@ function _default_freq_vector{T<:LTISystem}(systems::Vector{T}, plot::Symbol)
     min_pt_per_dec = 60
     min_pt_total = 200
     nsys = length(systems)
-    bounds = Array(Float64, 2, nsys)
+    bounds = Array{Float64}(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.
@@ -180,7 +180,7 @@ function _bounds_and_features(sys::LTISystem, plot::Symbol)
         zp = [tzero(sys); pole(sys)]
     end
     # Get the frequencies of the features, ignoring low frequency dynamics
-    fzp = log10(abs(zp))
+    fzp = log10.(abs.(zp))
     fzp = fzp[fzp .> -4]
     fzp = sort!(fzp)
     # Determine the bounds on the frequency vector

src/matrix_comps.jl

@@ -54,7 +54,7 @@ function dare(A, B, Q, R)
          -Ait*Q        Ait]
 
     S = schurfact(Z)
-    S = ordschur(S, abs(S.values).<=1)
+    S = ordschur(S, abs.(S.values).<=1)
     U = S.Z
 
     (m, n) = size(U)
@@ -181,14 +181,14 @@ end
 covar(sys::TransferFunction, W::StridedMatrix) = covar(ss(sys), W)
 
 
-# Note: the H∞ norm computation is probably not as accurate as with SLICOT, 
+# Note: the H∞ norm computation is probably not as accurate as with SLICOT,
 # but this seems to be still reasonably ok as a first step
 @doc """
 `..  norm(sys, p=2; tol=1e-6)`
 
 `norm(sys)` or `norm(sys,2)` computes the H2 norm of the LTI system `sys`.
 
-`norm(sys, Inf)` computes the L∞ norm of the LTI system `sys`. 
+`norm(sys, Inf)` computes the L∞ norm of the LTI system `sys`.
 The H∞ norm is the same as the L∞ for stable systems, and Inf for unstable systems.
 If the peak gain frequency is required as well, use the function `norminf` instead.
 
@@ -202,7 +202,7 @@ The L∞ norm computation implements the 'two-step algorithm' in:
 N.A. Bruinsma and M. Steinbuch, 'A fast algorithm to compute the H∞-norm
 of a transfer function matrix', Systems and Control Letters 14 (1990), pp. 287-293.
 For the discrete-time version, see, e.g.,: P. Bongers, O. Bosgra, M. Steinbuch, 'L∞-norm
-calculation for generalized state space systems in continuous and discrete time', 
+calculation for generalized state space systems in continuous and discrete time',
 American Control Conference, 1991.
 """ ->
 function Base.norm(sys::StateSpace, p::Real=2; tol=1e-6)
@@ -226,11 +226,11 @@ end
 @doc """
 `.. (peakgain, peakgainfrequency) = norminf(sys; tol=1e-6)`
 
-Compute the L∞ norm of the LTI system `sys`, together with the frequency 
+Compute the L∞ norm of the LTI system `sys`, together with the frequency
 `peakgainfrequency` (in rad/TimeUnit) at which the gain achieves its peak value `peakgain`.
 The H∞ norm is the same as the L∞ for stable systems, and Inf for unstable systems.
 
-`tol` is an optional keyword argument representing the desired relative accuracy for 
+`tol` is an optional keyword argument representing the desired relative accuracy for
 the computed L∞ norm (this is not an absolute certificate however).
 
 sys is first converted to a state space model if needed.
@@ -239,7 +239,7 @@ The L∞ norm computation implements the 'two-step algorithm' in:
 N.A. Bruinsma and M. Steinbuch, 'A fast algorithm to compute the H∞-norm
 of a transfer function matrix', Systems and Control Letters 14 (1990), pp. 287-293.
 For the discrete-time version, see, e.g.,: P. Bongers, O. Bosgra, M. Steinbuch, 'L∞-norm
-calculation for generalized state space systems in continuous and discrete time', 
+calculation for generalized state space systems in continuous and discrete time',
 American Control Conference, 1991.
 """ ->
 function norminf(sys::StateSpace; tol=1e-6)
@@ -277,9 +277,9 @@ function normLinf_twoSteps_ct(sys::StateSpace, tol=1e-6, maxIters=1000, approxim
             fpeak = 0
         end
         if isreal(p)  # only real poles
-            omegap = minimum(abs(p))
+            omegap = minimum(abs.(p))
         else  # at least one pair of complex poles
-            tmp = maximum(abs(imag(p)./(real(p).*abs(p))))
+            tmp = maximum(abs.(imag.(p)./(real.(p).*abs.(p))))
             omegap = abs(p[indmax(tmp)])
         end
         (lb, idx) = findmax([lb, maximum(svdvals(evalfr(sys, omegap*1im)))])
@@ -297,7 +297,7 @@ function normLinf_twoSteps_ct(sys::StateSpace, tol=1e-6, maxIters=1000, approxim
             H = [         M              -res*sys.B*(R\sys.B') ;
                    res*sys.C'*(S\sys.C)            -M'            ]
             omegas = eigvals(H)
-            omegaps = imag(omegas[ (abs(real(omegas)).<=approximag) & (imag(omegas).>=0) ])
+            omegaps = imag.(omegas[ (abs.(real.(omegas)).<=approximag) .& (imag.(omegas).>=0) ])
             sort!(omegaps)
             if isempty(omegaps)
                 return (1+tol)*lb, fpeak
@@ -324,7 +324,7 @@ function normLinf_twoSteps_dt(sys::StateSpace,tol=1e-6,maxIters=1000,approxcirc=
     end
     p = pole(sys)
     # Check first if there is a pole on the unit circle
-    pidx = findfirst(map(x->isapprox(x,1.0),abs(p)))
+    pidx = findfirst(map(x->isapprox(x,1.0),abs.(p)))
     if (pidx > 0)
         return (Inf, angle(p[pidx])/abs(sys.Ts))
     else
@@ -338,8 +338,8 @@ function normLinf_twoSteps_dt(sys::StateSpace,tol=1e-6,maxIters=1000,approxcirc=
         p = p[imag(p).>0]
         if ~isempty(p)  # not just real poles
             # find frequency of pôle closest to unit circle
-            omegap = angle(p[findmin(abs(abs(p)-1))[2]])
-        else 
+            omegap = angle(p[findmin(abs.(abs.(p)-1))[2]])
+        else
             omegap = pi/2
         end
         (lb, idx) = findmax([lb, maximum(svdvals(evalfr(sys, exp(omegap*1im))))])
@@ -355,11 +355,11 @@ function normLinf_twoSteps_dt(sys::StateSpace,tol=1e-6,maxIters=1000,approxcirc=
             RinvDt = R\sys.D'
             L = [ sys.A+sys.B*RinvDt*sys.C  sys.B*(R\sys.B');
                   zeros(sys.nx,sys.nx)      eye(sys.nx)]
-            M = [ eye(sys.nx)                              zeros(sys.nx,sys.nx); 
+            M = [ eye(sys.nx)                              zeros(sys.nx,sys.nx);
                   sys.C'*(eye(sys.ny)+sys.D*RinvDt)*sys.C  L[1:sys.nx,1:sys.nx]']
             zs = eig(L,M)[1]  # generalized eigenvalues
             # are there eigenvalues on the unit circle?
-            omegaps = angle(zs[ (abs(abs(zs)-1) .<= approxcirc) & (imag(zs).>=0)])
+            omegaps = angle.(zs[ (abs.(abs.(zs)-1) .<= approxcirc) .& (imag(zs).>=0)])
             sort!(omegaps)
             if isempty(omegaps)
                 return (1+tol)*lb, fpeak/sys.Ts
@@ -428,8 +428,8 @@ P = gram(sys, :c)
 Q = gram(sys, :o)
 Q1 = chol(Q)
 U,Σ,V = svd(Q1*P*Q1')
-Σ = sqrt(Σ)
-Σ1 = diagm((sqrt(Σ)))
+Σ .= sqrt.(Σ)
+Σ1 = diagm(sqrt.(Σ))
 T = Σ1\(U'Q1)
 
 Pz = T*P*T'

src/pid_design.jl

@@ -193,7 +193,7 @@ See also `stabregionPID`, `loopshapingPI`, `pidplots`
 """
 function stabregionPID(P, ω = _default_freq_vector(P,:bode); kd=0, doplot = true)
     Pv      = squeeze(freqresp(P,ω)[1],(2,3))
-    r       = abs(Pv)
+    r       = abs.(Pv)
     phi     = angle(Pv)
     kp      = -cos(phi)./r
     ki      = kd*ω.^2 - ω.*sin(phi)./r
@@ -205,7 +205,7 @@ function stabregionPID(P::AbstractString, ω = logspace(-3,1); kd=0, doplot = tr
     Pe      = parse(P)
     Pf(s)   = eval(:(s -> $(Pe)))(s)
     Pv      = Pf(im*ω)
-    r       = abs(Pv)
+    r       = abs.(Pv)
     phi     = angle(Pv)
     kp      = -cos(phi)./r
     ki      = kd*ω.^2 - ω.*sin(phi)./r
@@ -228,7 +228,7 @@ See also `pidplots`, `stabregionPID`
 function loopshapingPI(P,ωp; ϕl=0,rl=0, phasemargin = 0, doplot = false)
 Pw = P(im*ωp)[1]
 ϕp = angle(Pw)
-rp = abs(Pw)
+rp = abs.(Pw)
 
 if phasemargin > 0
     ϕl = deg2rad(-180+phasemargin)