diff --git a/.travis.yml b/.travis.yml index 6c886fc50..9b7f679c3 100644 --- a/.travis.yml +++ b/.travis.yml @@ -7,7 +7,7 @@ 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 'Pkg.clone("https://github.com/JuliaControl/ControlExamplePlots.jl.git");' + #- julia -e 'Pkg.clone("https://github.com/JuliaControl/ControlExamplePlots.jl.git");' - julia -e 'Pkg.test("ControlSystems"; coverage=true)' after_success: - julia -e 'cd(Pkg.dir("ControlSystems")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())' diff --git a/README.md b/README.md index b79ebbd08..fba824386 100644 --- a/README.md +++ b/README.md @@ -15,9 +15,27 @@ To install, in the Julia REPL: ```julia Pkg.add("ControlSystems") ``` - Note that the latest version of this package requires Julia 0.5. Users of Julia 0.4 should use [v0.1.4](https://github.com/JuliaControl/ControlSystems.jl/tree/v0.1.4) instead. +## News +### 2018-02-01 +- LTISystem types are now more generic and can hold matrices/vectors of arbitrary type. Examples (partly pseudo-code): +```julia +ss(1) +ss(1.) +ss(1im) +ss(ForwardDiff.Dual(1.)) +ss(GPUArray([1])) +ss(SparseMatrix([1]) +``` +Similar for `tf,zpk` etc. +- Continuous time systems are simulated with continuous time solvers from `OrdinaryDiffEq.jl` +- Breaking: `lsim(sys, u::Function)` syntax has changed from `u(t,x)` to `u(x,t)` to be consistent with `OrdinaryDiffEq` +- Breaking: `feedback(P,C)` no longer returns `feedback(P*C)`. The behavior is changed to `feedback(P1, P2) = P1/(1+P1*P2)`. +- Type `Simulator` provides lower level interface to continuous time simulation. +- Example [autodiff.jl](https://github.com/JuliaControl/ControlSystems.jl/tree/master/example/autodiff.jl) provides an illustration of how the new generic types can be used for automatic differentiation of a cost function through the continuous-time solver, which allows for optimization of the cost function w.r.t. PID parameters. + + ## Documentation All functions have docstrings, which can be viewed from the REPL, using for example `?tf `. diff --git a/REQUIRE b/REQUIRE index 49e794c12..729deaaf6 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,5 +1,6 @@ -julia 0.6.0 +julia 0.6.0 0.7.0 Plots 0.7.4 -Polynomials +Polynomials 0.3.0 LaTeXStrings -Requires +OrdinaryDiffEq +IterTools diff --git a/docs/src/lib/constructors.md b/docs/src/lib/constructors.md index 4f535b905..5991c582d 100644 --- a/docs/src/lib/constructors.md +++ b/docs/src/lib/constructors.md @@ -16,6 +16,5 @@ sminreal ss ss2tf tf -tfg zpk ``` diff --git a/docs/src/man/creatingtfs.md b/docs/src/man/creatingtfs.md index 13da426fa..284a38923 100644 --- a/docs/src/man/creatingtfs.md +++ b/docs/src/man/creatingtfs.md @@ -49,26 +49,6 @@ Continuous-time transfer function model The transfer functions created using this method will be of type `TransferFunction{SisoZpk}`. -## tfg - Generalized Representation -If you want to work with transfer functions that are not rational functions, it is possible to use the `tfg` representation -```julia -tfg(str::String), tfg(str::Expr) -``` -This function will either convert `str` to an expression or directly accept an `Expr` and create a transfer function. -### Example: -```julia -tfg("1/((s+1)*exp(-sqrt(s)))") - -## output - -TransferFunction: -1/((s+1)*exp(-sqrt(s))) - -Continuous-time transfer function model -``` -The transfer functions created using this method will be of type `TransferFunction{SisoGeneralized}`. -This type will work with some functions like `bodeplot, stepplot` but not others ,like `poles`. - ## Converting between types It is sometime useful to convert one representation to another, this is possible using the same functions, for example ```julia diff --git a/example/autodiff.jl b/example/autodiff.jl new file mode 100644 index 000000000..85d6224a9 --- /dev/null +++ b/example/autodiff.jl @@ -0,0 +1,114 @@ +using ControlSystems, OrdinaryDiffEq, NLopt, BlackBoxOptim +p0 = Float64[0.2,0.8,1] # Initial guess +K(kp,ki,kd) = pid(kp=kp, ki=ki, kd=kd) +K(p) = K(p...) + +# Define process model +ζ = 0.1 +ω = 1.; ω² = ω^2 +P = tf(ω²,[1, 2ζ*ω, ω²])*tf(1,[1,1]) + +Ω = logspace(-1,2,150) # Frequency vector to eval constraints +h = 0.1 # Sample time for time-domain evaluation +Tf = 60. # Time horizon +t = 0:h:Tf-h + +Ms = 1.4 # Maximum allowed magnitude of sensitivity function +Mt = 1.4 # Maximum allowed magnitude of complimentary sensitivity function + +p = copy(p0) + +function timedomain(p) + C = K(p[1], p[2], p[3]) + S = 1/(1+P*C) # Sensitivity fun + PS = ss(P*S) # TF from load disturbance to output + s = Simulator(PS, (t,x) -> [1]) # Sim. unit step load disturbance + ty = eltype(p) # So that all inputs to solve have same numerical type (ForwardDiff.Dual) + x0 = zeros(PS.nx) .|> ty + tspan = (ty(0.),ty(Tf)) + sol = solve(s, x0, tspan, Tsit5()) + y = PS.C*sol(t) # y = C*x + C,y +end + +function freqdomain(p) + C = K(p[1], p[2], p[3]) + S = 1/(1+P*C) # Sensitivity fun + T = tf(1.) - S# Comp. Sensitivity fun + Sw = vec(bode(S,Ω)[1]) # Freq. domain constraints + Tw = vec(bode(T,Ω)[1]) # Freq. domain constraints + Sw,Tw +end + +# Evaluates and plots candidate controller +function evalsol(p::Vector) + C,y = timedomain(p) + Sw,Tw = freqdomain(p) + plot(t,y', layout=2, show=false) + plot!(Ω, [Sw Tw] , lab=["Sw" "Tw"], subplot=2, xscale=:log10, yscale=:log10, show=false) + plot!([Ω[1],Ω[end]], [Ms,Ms], c = :black, l=:dash, subplot=2, show=false, lab="Ms") + plot!([Ω[1],Ω[end]], [Mt,Mt], c = :purple, l=:dash, subplot=2, lab="Mt") + gui() + false +end +evalsol(res::BlackBoxOptim.OptimizationResults) = evalsol(best_candidate(res)) + + +function constraintfun(p) + Sw,Tw = freqdomain(p) + [maximum(Sw)-Ms; maximum(Tw)-Mt] +end + + +function costfun(p) + C,y = timedomain(p) + mean(abs,y) # ~ Integrated absolute error IAE +end + + +function runopt(p, costfun, constraintfun; + f_tol = 1e-5, + x_tol = 1e-3, + c_tol = 1e-8, + f_cfg = ForwardDiff.GradientConfig(costfun, p), + g_cfg = ForwardDiff.JacobianConfig(constraintfun, p), + lb = zeros(length(p))) + + c1 = constraintfun(p) + np = length(p) + nc = length(c1) + + function f(p::Vector, grad::Vector) + if length(grad) > 0 + grad .= ForwardDiff.gradient(costfun,p,f_cfg) + end + costfun(p) + end + + function c(result, p::Vector, grad) + if length(grad) > 0 + grad .= ForwardDiff.jacobian(constraintfun,p,g_cfg)' + end + result .= constraintfun(p) + end + + opt = Opt(:LD_SLSQP, np) + lower_bounds!(opt, lb) + xtol_rel!(opt, x_tol) + ftol_rel!(opt, f_tol) + + min_objective!(opt, f) + inequality_constraint!(opt, c, c_tol*ones(nc)) + minf,minx,ret = NLopt.optimize(opt, p) +end + +f_cfg = ForwardDiff.GradientConfig(costfun, p) +g_cfg = ForwardDiff.JacobianConfig(constraintfun, p) +@time minf,minx,ret = runopt(1p0, costfun, constraintfun, x_tol=1e-6, c_tol=1e-12, f_cfg=f_cfg, g_cfg=g_cfg) +evalsol(minx) + +# # Optimize costfun using derivative-free method +# res1 = compare_optimizers(costfun; SearchRange = (0.,2.), NumDimensions = length(p0), MaxTime = 20.0) +# res2 = bboptimize(costfun, NumDimensions = length(p0), MaxTime = 20.0) +# evalsol(res2) +# p = best_candidate(res2) diff --git a/example/gks.svg b/example/gks.svg new file mode 100644 index 000000000..e69de29bb diff --git a/example/symbolic_computations.jl b/example/symbolic_computations.jl new file mode 100644 index 000000000..407a15540 --- /dev/null +++ b/example/symbolic_computations.jl @@ -0,0 +1,86 @@ +using Revise +using ControlSystems +using SymPy + +# Some basic demonstations of working with symbolic LTI systems +# Functionality is rather limited, and for complicated expressions the +# printing is awful. Anyway, let's get started... + + +# Need to modify some functions, to work with the Sym type +# Some of these changes are on their way into the respective packages +Base.complex(::Type{SymPy.Sym}) = Sym +Base.eigvals(a::Matrix{Sym}) = Vector{Sym}(collect(keys(call_matrix_meth(a, :eigenvals)))) + +function Polynomials.roots(p::Polynomials.Poly{SymPy.Sym}) + x = Sym("x") + return SymPy.solve(p(x), x) +end + +function SymPy.simplify(G::TransferFunction{ControlSystems.SisoZpk{Sym,Sym}}) + z, p, k = zpkdata(G) + + return zpk(map(x->simplify.(x), z), map(x->simplify.(x), p), map(x->simplify.(x), k)) +end + +function ControlSystems.impulse(sys::StateSpace{Sym}) + t = symbols("t", real=true) + return simplify.(sys.C * exp(sys.A*t) * sys.B) +end + + +# Define a symbolic parameter +a = symbols("a", real=true) + +# Define a statespace and a trasnfer function +sys = ss([1 a; a 1], [0; 1], [1 0], 0) + +s = tf("s") +G = (s/(5a)) / ((s/a)^2 + s/(5a) + 1) + + +# Simple conversions +@edit zpk(sys) +tf(sys) + +# The coefficient looks a bit complicated, but simplifying gives.. +z, p, k = zpkdata(tf(sys)) +simplify.(k[1]) + +zpkdata(G) +ss(G) + + +# Controllability/observability matrices +Mc = ctrb(sys) +Mo = obsv(sys) + + +## Compute the L∞ norm (or actually L∞ norm squared) for two symbolic systems +w = symbols("w", real=true) + +sys_to_consider = sys # G + +sys_fr = simplify(evalfr(sys_to_consider, im*w)[1]) +sys_fr_mag = simplify(abs(sys_fr)^2) + +n, _ = fraction( diff(sys_fr_mag, w) ) +roots = SymPy.solve(SymPy.Poly(n), w) + +real_roots = roots[SymPy.imag(roots) .== 0] + +maximum([subs(sys_fr_mag, w => r) for r in real_roots]) + + +# Compute the impulse resonse of some systems (on statespace form) +impulse(sys)[1] +simplify(impulse(ss(G))[1]) + +rosenbrock = [1/(s+1) 1/(s+a); 1/(s+1) 1/(s+1)] +ss(rosenbrock) +impulse(ss(rosenbrock)) + + +# Analytic impulse response +sys2 = ss([-2.5 0;1 1.5],[1;3],[1 2],Sym(2.5)) +impulse(sys2)[1] diff --git a/src/ControlSystems.jl b/src/ControlSystems.jl index 926d5f4bb..582c44e99 100644 --- a/src/ControlSystems.jl +++ b/src/ControlSystems.jl @@ -5,10 +5,10 @@ export LTISystem, TransferFunction, ss, tf, - tfg, zpk, ss2tf, LQG, + primitivetype, # Linear Algebra balance, care, @@ -57,6 +57,8 @@ export LTISystem, step, impulse, lsim, + solve, + Simulator, # Frequency Response freqresp, evalfr, @@ -67,27 +69,59 @@ export LTISystem, numpoly, denpoly -using Plots, LaTeXStrings, Requires + +# QUESTION: are these used? LaTeXStrings, Requires, IterTools +using Polynomials, OrdinaryDiffEq, Plots, LaTeXStrings import Base: +, -, *, /, (./), (==), (.+), (.-), (.*), (!=), isapprox, convert, promote_op +import Base.LinAlg: BlasFloat + +abstract type AbstractSystem end + +include("types/Lti.jl") + +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") +include("types/SisoTfTypes/conversion.jl") + +include("types/StateSpace.jl") + +# Convenience constructors +include("types/tf.jl") +include("types/zpk.jl") +include("types/ss.jl") -include("types/lti.jl") -include("types/transferfunction.jl") -include("types/statespace.jl") -include("types/tf2ss.jl") -include("types/lqg.jl") +include("types/lqg.jl") # QUESTION: is it really motivated to have an LQG type? +include("utilities.jl") + +include("types/promotion.jl") +include("types/conversion.jl") include("connections.jl") -include("discrete.jl") + +# Analysis +include("freqresp.jl") +include("timeresp.jl") + include("matrix_comps.jl") include("simplification.jl") -include("synthesis.jl") + +include("discrete.jl") include("analysis.jl") -include("timeresp.jl") -include("freqresp.jl") -include("utilities.jl") -include("plotting.jl") +include("synthesis.jl") + +include("simulators.jl") include("pid_design.jl") +include("plotting.jl") + + + # The path has to be evaluated upon initial import const __CONTROLSYSTEMS_SOURCE_DIR__ = dirname(Base.source_path()) diff --git a/src/analysis.jl b/src/analysis.jl index eb6652a5a..d540ee703 100644 --- a/src/analysis.jl +++ b/src/analysis.jl @@ -2,29 +2,17 @@ Compute the poles of system `sys`.""" -> pole(sys::StateSpace) = eigvals(sys.A) +# TODO wrong for MIMO pole(sys::TransferFunction) = [map(pole, sys.matrix)...;] pole(sys::SisoTf) = error("pole is not implemented for type $(typeof(sys))") @doc """`dcgain(sys)` -Compute the gain of SISO system `sys`.""" -> -function dcgain(sys::StateSpace, zs::Vector=tzero(sys)) - !issiso(sys) && error("Gain only defined for siso systems") - return ( sys.C*(-sys.A\sys.B) + sys.D)[1] -end - -@doc """`dcgain(sys)` - Compute the dcgain of system `sys`. equal to G(0) for continuous-time systems and G(1) for discrete-time systems.""" -> -function dcgain(sys::TransferFunction) - !issiso(sys) && error("Gain only defined for siso systems") - if sys.Ts > 0 - return map(s->sum(s.num.a)/sum(s.den.a), sys.matrix) - end - - return map(s -> s.num[end]/s.den[end], sys.matrix) +function dcgain(sys::LTISystem) + return iscontinuous(sys) ? evalfr(sys, 0) : evalfr(sys, 1) end @doc """`markovparam(sys, n)` @@ -51,39 +39,13 @@ 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) - for j = 1:nu - for i = 1:ny - zs[i, j], ps[i, j], ks[i, j] = _zpk_kern(sys, i, j) - end - end - return zs, ps, ks -end + G = convert(TransferFunction{SisoZpk}, sys) -function _zpk_kern(sys::StateSpace, iy::Int, iu::Int) - A, B, C = struct_ctrb_obsv(sys.A, sys.B[:, iu:iu], sys.C[iy:iy, :]) - D = sys.D[iy:iy, iu:iu] - z = tzero(A, B, C, D) - nx = size(A, 1) - nz = length(z) - k = nz == nx ? D[1] : (C*(A^(nx - nz - 1))*B)[1] - return z, eigvals(A), k -end + zs = getfield.(G.matrix, :z) + ps = getfield.(G.matrix, :p) + ks = getfield.(G.matrix, :k) -function _zpk_kern(sys::TransferFunction, iy::Int, iu::Int) - s = sys.matrix[iy, iu] - return _zpk_kern(s) -end - -function _zpk_kern(s::SisoRational) - return roots(s.num), roots(s.den), s.num[1]/s.den[1] -end - -function _zpk_kern(s::SisoZpk) - return s.z, s.p, s.k + return zs, ps, ks end @doc """`Wn, zeta, ps = damp(sys)` @@ -142,22 +104,27 @@ end # # Note that this returns either Vector{Complex64} or Vector{Float64} tzero(sys::StateSpace) = tzero(sys.A, sys.B, sys.C, sys.D) -function tzero(A::Matrix{Float64}, B::Matrix{Float64}, C::Matrix{Float64}, - D::Matrix{Float64}) +# Make sure everything is BlasFloat +function tzero(A::Matrix{<:Number}, B::Matrix{<:Number}, C::Matrix{<:Number}, D::Matrix{<:Number}) + T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D)) + A2, B2, C2, D2 = promote(A,B,C,D, fill(zero(T)/one(T),0,0)) # If Int, we get Float64 + tzero(A2, B2, C2, D2) +end +function tzero(A::Matrix{T}, B::Matrix{T}, C::Matrix{T}, D::Matrix{T}) where T<:BlasFloat # Balance the system A, B, C = balance_statespace(A, B, C) # Compute a good tolerance - meps = 10*eps()*norm([A B; C D]) + meps = 10*eps(T)*norm([A B; C D]) A, B, C, D = reduce_sys(A, B, C, D, meps) A, B, C, D = reduce_sys(A', C', B', D', meps) - if isempty(A) return Float64[] end + if isempty(A) return complex(T)[] end # Compress cols of [C D] to [0 Df] mat = [C D] # To ensure type-stability, we have to annote the type here, as qrfact # returns many different types. - W = full(qrfact(mat')[:Q], thin=false)::Matrix{Float64} + W = full(qrfact(mat')[:Q], thin=false)::Matrix{T} W = flipdim(W,2) mat = mat*W if fastrank(mat', meps) > 0 @@ -165,25 +132,28 @@ function tzero(A::Matrix{Float64}, B::Matrix{Float64}, C::Matrix{Float64}, m = size(D, 2) Af = ([A B] * W)[1:nf, 1:nf] Bf = ([eye(nf) zeros(nf, m)] * W)[1:nf, 1:nf] - zs = eig(Af, Bf)[1] + zs = eigvals(Af, Bf) + _fix_conjugate_pairs!(zs) # Generalized eigvals does not return exact conj. pairs else - zs = Float64[] + zs = complex(T)[] end return zs end +reduce_sys(A::Matrix{<:BlasFloat}, B::Matrix{<:BlasFloat}, C::Matrix{<:BlasFloat}, D::Matrix{<:BlasFloat}, meps::BlasFloat) = + reduce_sys(promote(A,B,C,D)..., meps) """ Implements REDUCE in the Emami-Naeini & Van Dooren paper. Returns transformed A, B, C, D matrices. These are empty if there are no zeros. """ -function reduce_sys(A::Matrix{Float64}, B::Matrix{Float64}, C::Matrix{Float64}, D::Matrix{Float64}, meps::Float64) +function reduce_sys(A::Matrix{T}, B::Matrix{T}, C::Matrix{T}, D::Matrix{T}, meps::BlasFloat) where T <: BlasFloat Cbar, Dbar = C, D if isempty(A) return A, B, C, D end while true # Compress rows of D - U = full(qrfact(D)[:Q], thin=false)::Matrix{Float64} + U = full(qrfact(D)[:Q], thin=false)::Matrix{T} D = U'*D C = U'*C sigma = fastrank(D, meps) @@ -195,7 +165,7 @@ function reduce_sys(A::Matrix{Float64}, B::Matrix{Float64}, C::Matrix{Float64}, end # Compress columns of Ctilde - V = full(qrfact(Ctilde')[:Q], thin=false)::Matrix{Float64} + V = full(qrfact(Ctilde')[:Q], thin=false)::Matrix{T} V = flipdim(V,2) Sj = Ctilde*V rho = fastrank(Sj', meps) @@ -205,15 +175,15 @@ 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 = Array{Float64,2}(0,0) + A = B = Cbar = Dbar = Matrix{T}(0,0) break end # Update System n, m = size(B) - Vm = [V zeros(n, m); zeros(m, n) eye(m)] + Vm = [V zeros(T, n, m); zeros(T, m, n) eye(T, m)] if sigma > 0 M = [A B; Cbar Dbar] - Vs = [V' zeros(n, sigma) ; zeros(sigma, n) eye(sigma)] + Vs = [V' zeros(T, n, sigma) ; zeros(T, sigma, n) eye(T, sigma)] else M = [A B] Vs = V' @@ -230,10 +200,10 @@ end # Determine the number of non-zero rows, with meps as a tolerance. For an # upper-triangular matrix, this is a good proxy for determining the row-rank. -function fastrank(A::Matrix{Float64}, meps::Float64) +function fastrank(A::AbstractMatrix, meps::Real) n, m = size(A, 1, 2) if n*m == 0 return 0 end - norms = Array{Float64}(n) + norms = Vector{eltype(A)}(n) for i = 1:n norms[i] = norm(A[i, :]) end @@ -256,11 +226,11 @@ function margin{S<:Real}(sys::LTISystem, w::AbstractVector{S}; full=false, allMa vals = (:wgm, :gm, :wpm, :pm, :fullPhase) if allMargins for val in vals - eval(:($val = Array{Array{Float64,1}}($ny,$nu))) + eval(:($val = Array{Array{eltype(sys),1}}($ny,$nu))) end else for val in vals - eval(:($val = Array{Float64,2}($ny,$nu))) + eval(:($val = Array{eltype(sys),2}($ny,$nu))) end end for j=1:nu @@ -326,9 +296,9 @@ margin(system, _default_freq_vector(system, :bode); kwargs...) # Interpolate the values in "list" given the floating point "index" fi function interpolate(fi, list) - fif = floor(Integer, fi) - fic = ceil(Integer, fi) - list[fif]+mod(fi,1).*(list[fic]-list[fif]) + fif = floor.(Integer, fi) + fic = ceil.(Integer, fi) + list[fif]+mod.(fi,1).*(list[fic]-list[fif]) end function _allGainCrossings(w, mag) @@ -337,8 +307,8 @@ end function _allPhaseCrossings(w, phase) #Calculate numer of times real axis is crossed on negative side - n = Array{Float64,1}(length(w)) #Nbr of crossed - ph = Array{Float64,1}(length(w)) #Residual + n = Vector{eltype(w)}(length(w)) #Nbr of crossed + ph = Vector{eltype(w)}(length(w)) #Residual for i = eachindex(w) #Found no easier way to do this n[i], ph[i] = fldmod(phase[i]+180,360)#+180 end @@ -346,8 +316,8 @@ function _allPhaseCrossings(w, phase) end function _findCrossings(w, n, res) - wcross = Array{Float64,1}() - tcross = Array{Float64,1}() + wcross = Vector{eltype(w)}() + tcross = Vector{eltype(w)}() for i in 1:(length(w)-1) if res[i] == 0 wcross = [wcross; w[i]] @@ -375,11 +345,11 @@ function delaymargin(G::LTISystem) if G.nu + G.ny > 2 error("delaymargin only supports SISO systems") end - m = margin(G,allMargins=true) - ϕₘ,i= findmin(m[4]) - ϕₘ *= π/180 - ωϕₘ = m[3][i] - dₘ = ϕₘ/ωϕₘ + m = margin(G,allMargins=true) + ϕₘ, i = findmin(m[4]) + ϕₘ *= π/180 + ωϕₘ = m[3][i] + dₘ = ϕₘ/ωϕₘ if G.Ts > 0 dₘ /= G.Ts # Give delay margin in number of sample times, as matlab does end @@ -410,7 +380,9 @@ function gangoffour(P::TransferFunction,C::TransferFunction) return S, D, N, T end + function gangoffour(P::AbstractVector, C::AbstractVector) + Base.depwarn("Deprecrated use of gangoffour(::Vector, ::Vector), use `broadcast` and `zip` instead", :gangoffour) if P[1].nu + P[1].ny + C[1].nu + C[1].ny > 4 error("gangoffour only supports SISO systems") end @@ -424,10 +396,12 @@ function gangoffour(P::AbstractVector, C::AbstractVector) end function gangoffour(P::TransferFunction, C::AbstractVector) + Base.depwarn("Deprecrated use of gangoffour(::TransferFunction, ::Vector), use `broadcast` and `zip` instead", :gangoffour) gangoffour(fill(P,length(C)), C) end function gangoffour(P::AbstractVector, C::TransferFunction) + Base.depwarn("Deprecrated use of gangoffour(::Vector, ::TransferFunction), use `broadcast` and `zip` instead", :gangoffour) gangoffour(P, fill(C,length(P))) end diff --git a/src/connections.jl b/src/connections.jl index fc5214e5f..7e5156cf3 100644 --- a/src/connections.jl +++ b/src/connections.jl @@ -1,18 +1,18 @@ # Model interconnections @doc """ -`series(s1::LTISystem, s2::LTISystem)` +`series(sys1::LTISystem, sys2::LTISystem)` -Connect systems in series, equivalent to `s2*s1` +Connect systems in series, equivalent to `sys2*sys1` """ -> -series(s1::LTISystem, s2::LTISystem) = s2*s1 +series(sys1::LTISystem, sys2::LTISystem) = sys2*sys1 @doc """ -`series(s1::LTISystem, s2::LTISystem)` +`series(sys1::LTISystem, sys2::LTISystem)` -Connect systems in parallel, equivalent to `s2+s1` +Connect systems in parallel, equivalent to `sys2+sys1` """ -> -parallel(s1::LTISystem, s2::LTISystem) = s1 + s2 +parallel(sys1::LTISystem, sys2::LTISystem) = sys1 + sys2 append() = LTISystem[] @doc """ @@ -22,28 +22,23 @@ Append systems in block diagonal form """ -> function append(systems::StateSpace...) Ts = systems[1].Ts - if !all([s.Ts == Ts for s in systems]) + 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]...) - states = vcat([s.statenames for s in systems]...) - outputs = vcat([s.outputnames for s in systems]...) - inputs = vcat([s.inputnames for s in systems]...) - return StateSpace(A, B, C, D, Ts, states, inputs, outputs) + return StateSpace(A, B, C, D, Ts) end function append(systems::TransferFunction...) Ts = systems[1].Ts - if !all([s.Ts == Ts for s in systems]) + if !all(s.Ts == Ts for s in systems) error("Sampling time mismatch") end mat = blkdiag([s.matrix for s in systems]...) - inputs = vcat([s.inputnames for s in systems]...) - outputs = vcat([s.outputnames for s in systems]...) - return TransferFunction(mat, Ts, inputs, outputs) + return TransferFunction(mat, Ts) end append(systems::LTISystem...) = append(promote(systems...)...) @@ -51,107 +46,117 @@ append(systems::LTISystem...) = append(promote(systems...)...) function Base.vcat(systems::StateSpace...) # Perform checks nu = systems[1].nu - if !all([s.nu == nu for s in systems]) + if !all(s.nu == nu for s in systems) error("All systems must have same input dimension") end Ts = systems[1].Ts - if !all([s.Ts == Ts for s in systems]) + if !all(s.Ts == Ts for s in systems) error("Sampling time mismatch") end A = blkdiag([s.A for s in systems]...) B = vcat([s.B for s in systems]...) C = blkdiag([s.C for s in systems]...) D = vcat([s.D for s in systems]...) - states = vcat([s.statenames for s in systems]...) - outputs = vcat([s.outputnames for s in systems]...) - inputs = systems[1].inputnames - if !all([s.inputnames == inputs for s in systems]) - inputs = fill(String(""),size(inputs, 1)) - end - return StateSpace(A, B, C, D, Ts, states, inputs, outputs) + + return StateSpace(A, B, C, D, Ts) end function Base.vcat(systems::TransferFunction...) # Perform checks nu = systems[1].nu - if !all([s.nu == nu for s in systems]) + if !all(s.nu == nu for s in systems) error("All systems must have same input dimension") end Ts = systems[1].Ts - if !all([s.Ts == Ts for s in systems]) + if !all(s.Ts == Ts for s in systems) error("Sampling time mismatch") end mat = vcat([s.matrix for s in systems]...) - outputs = vcat([s.outputnames for s in systems]...) - inputs = systems[1].inputnames - if !all([s.inputnames == inputs for s in systems]) - inputs = fill(String(""),size(inputs, 1)) - end - return TransferFunction(mat, Ts, inputs, outputs) + return TransferFunction(mat, Ts) end Base.vcat(systems::LTISystem...) = vcat(promote(systems...)...) -function Base.vcat{T<:Real}(systems::Union{VecOrMat{T},T,TransferFunction}...) - if promote_type(map(e->typeof(e),systems)...) <: TransferFunction - vcat(map(e->convert(TransferFunction,e),systems)...) - else - cat(1,systems...) - end -end - function Base.hcat(systems::StateSpace...) # Perform checks ny = systems[1].ny - if !all([s.ny == ny for s in systems]) + if !all(s.ny == ny for s in systems) error("All systems must have same output dimension") end Ts = systems[1].Ts - if !all([s.Ts == Ts for s in systems]) + 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 = hcat([s.C for s in systems]...) D = hcat([s.D for s in systems]...) - states = vcat([s.statenames for s in systems]...) - inputs = vcat([s.inputnames for s in systems]...) - outputs = systems[1].outputnames - if !all([s.outputnames == outputs for s in systems]) - outputs = fill(String(""),size(outputs, 1)) - end - return StateSpace(A, B, C, D, Ts, states, inputs, outputs) + + return StateSpace(A, B, C, D, Ts) end function Base.hcat(systems::TransferFunction...) # Perform checks ny = systems[1].ny - if !all([s.ny == ny for s in systems]) + if !all(s.ny == ny for s in systems) error("All systems must have same output dimension") end Ts = systems[1].Ts - if !all([s.Ts == Ts for s in systems]) + if !all(s.Ts == Ts for s in systems) error("Sampling time mismatch") end mat = hcat([s.matrix for s in systems]...) - inputs = vcat([s.inputnames for s in systems]...) - outputs = systems[1].outputnames - if !all([s.outputnames == outputs for s in systems]) - outputs = fill(String(""),size(outputs, 1)) - end - return TransferFunction(mat, Ts, inputs, outputs) + return TransferFunction(mat, Ts) end Base.hcat(systems::LTISystem...) = hcat(promote(systems...)...) -function Base.hcat{T<:Real}(systems::Union{T,VecOrMat{T},TransferFunction}...) - if promote_type(map(e->typeof(e),systems)...) <: TransferFunction - hcat(map(e->convert(TransferFunction,e),systems)...) - else - cat(2,systems...) - end +function Base.cat_t(::Type{Val{1}}, T::Type{<:LTISystem}, X...) + vcat(convert.(T, X)...) +end + +function Base.cat_t(::Type{Val{2}}, T::Type{<:LTISystem}, X...) + hcat(convert.(T, X)...) end +# function hvcat(rows::Tuple{Vararg{Int}}, systems::Union{Number,AbstractVecOrMat{<:Number},LTISystem}...) +# T = Base.promote_typeof(systems...) +# nbr = length(rows) # number of block rows +# rs = Array{T,1}(nbr) +# a = 1 +# for i = 1:nbr +# rs[i] = hcat(convert.(T,systems[a:a-1+rows[i]])...) +# a += rows[i] +# end +# vcat(rs...) +# end + +# function _get_common_sampling_time(sys_vec::Union{AbstractVector{LTISystem},AbstractVecOrMat{<:Number},Number}) +# Ts = -1.0 # Initalize corresponding to undefined sampling time +# +# for sys in sys_vec +# if !all(s.Ts == Ts for s in systems]) +# error("Sampling time mismatch") +# end +# end +# +# end + + +# function Base.hcat{T<:Number}(systems::Union{T,AbstractVecOrMat{T},TransferFunction}...) +# S = promote_type(map(e->typeof(e),systems)...) # TODO: Should be simplified +# +# idx_first_tf = findfirst(e -> isa(e, TransferFunction), systems) +# Ts = sys_tuple[idx_first_tf].Ts +# +# if S <: TransferFunction +# hcat(map(e->convert(TransferFunction,e),systems)...) +# else +# cat(2,systems...) +# end +# end + +# TODO: could use cat([1,2], mats...) instead # Empty definition to get rid of warning Base.blkdiag() = [] function Base.blkdiag(mats::Matrix...) diff --git a/src/discrete.jl b/src/discrete.jl index 3e18f43d5..2a099bf05 100644 --- a/src/discrete.jl +++ b/src/discrete.jl @@ -13,9 +13,9 @@ function c2d(sys::StateSpace, Ts::Real, method::Symbol=:zoh) if !iscontinuous(sys) error("sys must be a continuous time system") end - A, B, C, D = sys.A, sys.B, sys.C, sys.D + A, B, C, D = ssdata(sys) ny, nu = size(sys) - nx = sys.nx + nx = nstates(sys) if method == :zoh M = expm([A*Ts B*Ts; zeros(nu, nx + nu)]) @@ -40,8 +40,7 @@ function c2d(sys::StateSpace, Ts::Real, method::Symbol=:zoh) else error("Unsupported method: ", method) end - return StateSpace(Ad, Bd, Cd, Dd, Ts, sys.statenames, sys.inputnames, - sys.outputnames), x0map + return StateSpace(Ad, Bd, Cd, Dd, Ts), x0map end @@ -189,13 +188,13 @@ function c2d_poly2poly(p,h) end -function c2d(G::TransferFunction, h;kwargs...) +function c2d(G::TransferFunction{S}, h;kwargs...) where {S} @assert iscontinuous(G) ny, nu = size(G) @assert (ny + nu == 2) "c2d(G::TransferFunction, h) not implemented for MIMO systems" sys = ss(G) - sysd = c2d(sys,h,kwargs...)[1] - return ss2tf(sysd) + sysd = c2d(sys, h, kwargs...)[1] + return convert(TransferFunction, sysd) end diff --git a/src/freqresp.jl b/src/freqresp.jl index db3dc91ba..e646a1701 100644 --- a/src/freqresp.jl +++ b/src/freqresp.jl @@ -5,24 +5,29 @@ Evaluate the frequency response of a linear system `w -> C*((iw*im -A)^-1)*B + D` of system `sys` over the frequency vector `w`.""" -> -function freqresp{S<:Real}(sys::LTISystem, w::AbstractVector{S}) - ny, nu = size(sys) - nw = length(w) +function freqresp(sys::LTISystem, w_vec::AbstractVector{S}) where {S<:Real} # Create imaginary freq vector s if !iscontinuous(sys) Ts = sys.Ts == -1 ? 1.0 : sys.Ts - s = exp.(w.*(im*Ts)) + s_vec = exp.(w_vec*(im*Ts)) else - s = im*w + s_vec = im*w_vec end - #Evil but nessesary type instability here - sys = _preprocess_for_freqresp(sys) - sys_fr = Array{Complex{eltype(w)}}(nw, ny, nu) - for i=1:nw + 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)) + + if isa(sys, StateSpace) + sys = _preprocess_for_freqresp(sys) + end + + + for i=1:length(s_vec) # TODO : This doesn't actually take advantage of Hessenberg structure # for statespace version. - sys_fr[i, :, :] = evalfr(sys, s[i]) + sys_fr[i, :, :] .= evalfr(sys, s_vec[i]) end + return sys_fr end @@ -33,22 +38,25 @@ function _preprocess_for_freqresp(sys::StateSpace) if isempty(sys.A) # hessfact does not work for empty matrices return sys end + Tsys = numeric_type(sys) + TT = promote_type(typeof(zero(Tsys)/norm(one(Tsys))), Float32) A, B, C, D = sys.A, sys.B, sys.C, sys.D F = hessfact(A) - H = F[:H]::Matrix{Float64} + H = F[:H]::Matrix{TT} T = full(F[:Q]) P = C*T Q = T\B - StateSpace(H, Q, P, D, sys.Ts, sys.statenames, sys.inputnames, - sys.outputnames) + StateSpace(H, Q, P, D, sys.Ts) end -function _preprocess_for_freqresp(sys::TransferFunction) - map(sisotf -> _preprocess_for_freqresp(sisotf), sys.matrix) -end -_preprocess_for_freqresp(sys::SisoTf) = sys +#_preprocess_for_freqresp(sys::TransferFunction) = sys.matrix +#function _preprocess_for_freqresp(sys::TransferFunction) +# map(sisotf -> _preprocess_for_freqresp(sisotf), sys.matrix) +#end + +#_preprocess_for_freqresp(sys::SisoTf) = sys @doc """ `evalfr(sys, x)` Evaluate the transfer function of the LTI system sys @@ -56,37 +64,33 @@ at the complex number s=x (continuous-time) or z=x (discrete-time). For many values of `x`, use `freqresp` instead. """ -> -function evalfr(sys::StateSpace, s::Number) - S = promote_type(typeof(s), Float64) +function evalfr(sys::StateSpace{T0}, s::Number) where {T0} + T = promote_type(T0, typeof(one(T0)*one(typeof(s))/(one(T0)*one(typeof(s))))) try R = s*I - sys.A - sys.D + sys.C*((R\sys.B)::Matrix{S}) # Weird type stability issue + sys.D + sys.C*((R\sys.B)::Matrix{T}) # Weird type stability issue catch - fill(convert(S, Inf), size(sys)) + fill(convert(T, Inf), size(sys)) end end -function evalfr(sys::TransferFunction, s::Number) - S = promote_type(typeof(s), Float64) - mat = sys.matrix - ny, nu = size(mat) - res = Array{S}(ny, nu) - for j = 1:nu - for i = 1:ny - res[i, j] = evalfr(mat[i, j], s) +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)) + for j = 1:ninputs(G) + for i = 1:noutputs(G) + fr[i, j] = evalfr(G.matrix[i, j], s) end end - return res + return fr end -evalfr(mat::Matrix, s::Number) = map(sys -> evalfr(sys, s), mat) - @doc """ `F(s)`, `F(omega, true)`, `F(z, false)` Notation for frequency response evaluation. - F(s) evaluates the continuous-time transfer function F at s. -- F(omega,true) evaluates the discrete-time transfer function F at i*Ts*omega +- F(omega,true) evaluates the discrete-time transfer function F at exp(i*Ts*omega) - F(z,false) evaluates the discrete-time transfer function F at z """ -> function (sys::TransferFunction)(s) @@ -102,9 +106,15 @@ function (sys::TransferFunction)(z_or_omega::Number, map_to_unit_circle::Bool) end end -function (sys::TransferFunction)(s::AbstractVector, map_to_unit_circle::Bool) +function (sys::TransferFunction)(z_or_omegas::AbstractVector, map_to_unit_circle::Bool) @assert !iscontinuous(sys) "It makes no sense to call this function with continuous systems" - freqresp(sys,s) + 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)...) + for i in 1:length(vals) + out[i,:,:] .= vals[i] + end + return out end @doc """`mag, phase, w = bode(sys[, w])` @@ -166,10 +176,7 @@ end _default_freq_vector(sys::LTISystem, plot::Symbol) = _default_freq_vector( LTISystem[sys], plot) -_default_freq_vector{T<:TransferFunction{SisoGeneralized}}(sys::Vector{T}, plot::Symbol) = - logspace(-2,2,400) -_default_freq_vector(sys::TransferFunction{SisoGeneralized} , plot::Symbol) = - logspace(-2,2,400) + function _bounds_and_features(sys::LTISystem, plot::Symbol) diff --git a/src/matrix_comps.jl b/src/matrix_comps.jl index cb706f19f..53e85b15d 100644 --- a/src/matrix_comps.jl +++ b/src/matrix_comps.jl @@ -86,9 +86,10 @@ function gram(sys::StateSpace, opt::Symbol) end func = iscontinuous(sys) ? lyap : dlyap if opt == :c - return func(sys.A, sys.B*sys.B') + # TODO probably remove type check in julia 0.7.0 + return func(sys.A, sys.B*sys.B')::Array{numeric_type(sys),2} # lyap is type-unstable elseif opt == :o - return func(sys.A', sys.C'*sys.C) + return func(sys.A', sys.C'*sys.C)::Array{numeric_type(sys),2} # lyap is type-unstable else error("opt must be either :c for controllability grammian, or :o for observability grammian") @@ -102,13 +103,13 @@ Compute the observability matrix for the system described by `(A, C)` or `sys`. Note that checking for observability by computing the rank from `obsv` is not the most numerically accurate way, a better method is checking if `gram(sys, :o)` is positive definite.""" -> -function obsv(A, C) +function obsv(A::AbstractMatrix{T}, C::AbstractMatrix{T}) where {T <: Number} n = size(A, 1) ny = size(C, 1) if n != size(C, 2) error("C must have the same number of columns as A") end - res = zeros(n*ny, n) + res = fill(zero(T), n*ny, n) res[1:ny, :] = C for i=1:n-1 res[(1 + i*ny):(1 + i)*ny, :] = res[((i - 1)*ny + 1):i*ny, :] * A @@ -125,13 +126,13 @@ Compute the controllability matrix for the system described by `(A, B)` or Note that checking for controllability by computing the rank from `ctrb` is not the most numerically accurate way, a better method is checking if `gram(sys, :c)` is positive definite.""" -> -function ctrb(A, B) +function ctrb(A::AbstractMatrix{T}, B::AbstractMatrix{T}) where {T <: Number} n = size(A, 1) nu = size(B, 2) if n != size(B, 1) error("B must have the same number of rows as A") end - res = zeros(n, n*nu) + res = fill(zero(T), n, n*nu) res[:, 1:nu] = B for i=1:n-1 res[:, (1 + i*nu):(1 + i)*nu] = A * res[:, ((i - 1)*nu + 1):i*nu] @@ -260,21 +261,23 @@ function normLinf_twoSteps_ct(sys::StateSpace, tol=1e-6, maxIters=1000, approxim # Because of this tuning for example, the relative precision that we provide on the norm computation # is not a true guarantee, more an order of magnitude # outputs: pair of Float64, namely L∞ norm approximation and frequency fpeak at which it is achieved + T = promote_type(numeric_type(sys), Float32) if sys.nx == 0 # static gain - return (norm(sys.D,2), 0.0) + return (norm(sys.D,2), T(0)) end p = pole(sys) # Check if there is a pole on the imaginary axis pidx = findfirst(map(x->isapprox(x,0.0),real(p))) if pidx > 0 - return (Inf, imag(p[pidx])) + return (T(Inf), imag(p[pidx])) # note: in case of cancellation, for s/s for example, we return Inf, whereas Matlab returns 1 else # Initialization: computation of a lower bound from 3 terms - lb = maximum(svdvals(sys.D)); fpeak = Inf - (lb, idx) = findmax([lb, maximum(svdvals(evalfr(sys,0)))]) + lb = maximum(svdvals(sys.D)) + fpeak = T(Inf) + (lb, idx) = findmax([lb, T(maximum((svdvals(evalfr(sys,0)))))]) #TODO remove T() in julia 0.7.0 if idx == 2 - fpeak = 0 + fpeak = T(0) end if isreal(p) # only real poles omegap = minimum(abs.(p)) @@ -282,7 +285,7 @@ function normLinf_twoSteps_ct(sys::StateSpace, tol=1e-6, maxIters=1000, approxim 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)))]) + (lb, idx) = findmax([lb, T(maximum(svdvals(evalfr(sys, omegap*1im))))]) #TODO remove T() in julia 0.7.0 if idx == 2 fpeak = omegap end @@ -290,21 +293,21 @@ function normLinf_twoSteps_ct(sys::StateSpace, tol=1e-6, maxIters=1000, approxim # Iterations iter = 1; while iter <= maxIters - res = (1+2tol)*lb - R = sys.D'*sys.D - res^2*eye(sys.nu) - S = sys.D*sys.D' - res^2*eye(sys.ny) + res = (1+2*T(tol))*lb + R = sys.D'*sys.D - res^2*eye(T, sys.nu) + S = sys.D*sys.D' - res^2*eye(T, sys.ny) M = sys.A-sys.B*(R\sys.D')*sys.C H = [ M -res*sys.B*(R\sys.B') ; res*sys.C'*(S\sys.C) -M' ] - omegas = eigvals(H) + omegas = eigvals(H) .+ 0im # To make type stable omegaps = imag.(omegas[ (abs.(real.(omegas)).<=approximag) .& (imag.(omegas).>=0) ]) sort!(omegaps) if isempty(omegaps) - return (1+tol)*lb, fpeak + return (1+T(tol))*lb, fpeak else # if not empty, omegaps contains at least two values ms = [(x+y)/2 for x=omegaps[1:end-1], y=omegaps[2:end]] for mval in ms - (lb, idx) = findmax([lb, maximum(svdvals(evalfr(sys,mval*1im)))]) + (lb, idx) = findmax([lb, T(maximum(svdvals(evalfr(sys,mval*1im))))]) #TODO remove T() in julia 0.7.0 if idx == 2 fpeak = mval end @@ -312,27 +315,29 @@ function normLinf_twoSteps_ct(sys::StateSpace, tol=1e-6, maxIters=1000, approxim end iter += 1 end - println("The computation of the H-infinity norm did not converge in $maxIters iterations") + error("In norminf: The computation of the H-infinity norm did not converge in $maxIters iterations") end end # discrete-time version of normHinf_twoSteps_ct above # The value fpeak returned by the function is in the range [0,pi)/sys.Ts (in rad/s) -function normLinf_twoSteps_dt(sys::StateSpace,tol=1e-6,maxIters=1000,approxcirc=1e-8) +function normLinf_twoSteps_dt(sys::StateSpace,tol=1e-6, maxIters=1000, approxcirc=1e-8) + T = promote_type(numeric_type(sys), Float32) if sys.nx == 0 # static gain - return (norm(sys.D,2), 0.0) + return (norm(sys.D,2), T(0)) 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))) if (pidx > 0) - return (Inf, angle(p[pidx])/abs(sys.Ts)) + return (T(Inf), angle(p[pidx])/abs(T(sys.Ts))) else # Initialization: computation of a lower bound from 3 terms - lb = maximum(svdvals(evalfr(sys,1))); fpeak = 0 - (lb, idx) = findmax([lb, maximum(svdvals(evalfr(sys,-1)))]) + lb = T(maximum(svdvals(evalfr(sys,1)))) #TODO remove T() in julia 0.7.0 + fpeak = T(0) + (lb, idx) = findmax([lb, T(maximum(svdvals(evalfr(sys,-1))))]) #TODO remove T() in julia 0.7.0 if idx == 2 - fpeak = pi + fpeak = T(pi) end p = p[imag(p).>0] @@ -340,9 +345,9 @@ function normLinf_twoSteps_dt(sys::StateSpace,tol=1e-6,maxIters=1000,approxcirc= # find frequency of pôle closest to unit circle omegap = angle(p[findmin(abs.(abs.(p)-1))[2]]) else - omegap = pi/2 + omegap = T(pi)/2 end - (lb, idx) = findmax([lb, maximum(svdvals(evalfr(sys, exp(omegap*1im))))]) + (lb, idx) = findmax([lb, T(maximum(svdvals(evalfr(sys, exp(omegap*1im)))))]) #TODO remove T() in julia 0.7.0 if idx == 2 fpeak = omegap end @@ -350,23 +355,24 @@ function normLinf_twoSteps_dt(sys::StateSpace,tol=1e-6,maxIters=1000,approxcirc= # Iterations iter = 1; while iter <= maxIters - res = (1+2tol)*lb - R = res^2*eye(sys.nu) - sys.D'*sys.D + res = (1+2*T(tol))*lb + R = res^2*eye(T, sys.nu) - sys.D'*sys.D 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); - sys.C'*(eye(sys.ny)+sys.D*RinvDt)*sys.C L[1:sys.nx,1:sys.nx]'] - zs = eig(L,M)[1] # generalized eigenvalues + zeros(T, sys.nx,sys.nx) eye(T, sys.nx)] + M = [ eye(T, sys.nx) zeros(T, sys.nx,sys.nx); + sys.C'*(eye(T, sys.ny)+sys.D*RinvDt)*sys.C L[1:sys.nx,1:sys.nx]'] + # +0im to make type stable + zs = eigvals(L,M) .+ 0im # generalized eigenvalues # are there eigenvalues on the unit circle? omegaps = angle.(zs[ (abs.(abs.(zs)-1) .<= approxcirc) .& (imag(zs).>=0)]) sort!(omegaps) if isempty(omegaps) - return (1+tol)*lb, fpeak/sys.Ts + return (1+T(tol))*lb, fpeak/T(sys.Ts) else # if not empty, omegaps contains at least two values ms = [(x+y)/2 for x=omegaps[1:end-1], y=omegaps[2:end]] for mval in ms - (lb, idx) = findmax([lb, maximum(svdvals(evalfr(sys,exp(mval*1im))))]) + (lb, idx) = findmax([lb, T(maximum(svdvals(evalfr(sys,exp(mval*1im)))))]) #TODO remove T() in julia 0.7.0 if idx == 2 fpeak = mval end @@ -374,7 +380,7 @@ function normLinf_twoSteps_dt(sys::StateSpace,tol=1e-6,maxIters=1000,approxcirc= end iter += 1 end - println("The computation of the H-infinity norm did not converge in $maxIters iterations") + error("In norminf: The computation of the H-infinity norm did not converge in $maxIters iterations") end end diff --git a/src/pid_design.jl b/src/pid_design.jl index 768e79159..3bd090b76 100644 --- a/src/pid_design.jl +++ b/src/pid_design.jl @@ -7,12 +7,12 @@ Calculates and returns a PID controller on transfer function form. `C = pid(; kp=0, ki=0; kd=0, time=false, series=false)` """ -function pid(; kp=0, ki=0, kd=0, time=false, series=false) +function pid(; kp=0., ki=0., kd=0., time=false, series=false) s = tf("s") if series - return time ? kp*(1 + 1/(ki*s) + kd*s) : kp*(1 + ki/s + kd*s) + return time ? kp*(one(kp) + one(kp)/(ki*s) + kd*s) : kp*(one(kp) + ki/s + kd*s) else - return time ? kp + 1/(ki*s) + kd*s : kp + ki/s + kd*s + return time ? kp + one(kp)/(ki*s) + kd*s : kp + ki/s + kd*s end end @@ -95,17 +95,17 @@ end @userplot Rlocusplot @deprecate rlocus(args...;kwargs...) rlocusplot(args...;kwargs...) -function getpoles(G,K) - poles = map(k -> pole(feedback(G,tf(k))), K) - cat(2,poles...)' -end +# function getpoles(G,K) +# poles = map(k -> pole(feedback(G,tf(k))), K) +# cat(2,poles...)' +# end + -@require OrdinaryDiffEq begin -function getpoles(G, K) # If OrdinaryDiffEq is installed, we overrides getpoles with an adaptive method +function getpoles(G, K) # If OrdinaryDiffEq is installed, we override getpoles with an adaptive method P = G.matrix[1].num.a |> reverse |> Polynomials.Poly Q = G.matrix[1].den.a |> reverse |> Polynomials.Poly - f = (k,y) -> Complex128.(Polynomials.roots(k[1]*P+Q)) - prob = OrdinaryDiffEq.ODEProblem(f,f(0.,0.),(0.,K[end])) + f = (y,_,k) -> Complex128.(Polynomials.roots(k[1]*P+Q)) + prob = OrdinaryDiffEq.ODEProblem(f,f(0.,0.,0.),(0.,K[end])) integrator = OrdinaryDiffEq.init(prob,OrdinaryDiffEq.Tsit5(),reltol=1e-8,abstol=1e-8) poleout = Vector{Vector{Complex128}}() for i in integrator @@ -114,7 +114,7 @@ function getpoles(G, K) # If OrdinaryDiffEq is installed, we overrides getpoles poleout = hcat(poleout...)' end -end + """ rlocusplot(P::LTISystem, K) diff --git a/src/plotting.jl b/src/plotting.jl index a58643235..fb1bb842e 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -103,7 +103,7 @@ lsimplot error("All systems must have the same input/output dimensions") end ny, nu = size(systems[1]) - layout := (ny,1) + layout --> (ny,1) s2i(i,j) = sub2ind((ny,1),j,i) for (si,s) in enumerate(systems) s = systems[si] @@ -116,7 +116,7 @@ lsimplot xguide --> "Time (s)" yguide --> ytext title --> "System Response" - subplot := s2i(1,i) + subplot --> s2i(1,i) label --> "\$G_\{$(si)\}\$" linestyle --> styledict[:l] linecolor --> styledict[:c] @@ -162,7 +162,7 @@ for (func, title, typ) = ((step, "Step Response", Stepplot), (impulse, "Impulse error("All systems must have the same input/output dimensions") end ny, nu = size(systems[1]) - layout := (ny,nu) + layout --> (ny,nu) titles = fill("", 1, ny*nu) title --> titles s2i(i,j) = sub2ind((ny,nu),i,j) @@ -181,7 +181,7 @@ for (func, title, typ) = ((step, "Step Response", Stepplot), (impulse, "Impulse seriestype := style xlabel --> "Time (s)" ylabel --> ytext - subplot := s2i(i,j) + subplot --> s2i(i,j) label --> "\$G_\{$(si)\}\$" linestyle --> styledict[:l] linecolor --> styledict[:c] @@ -208,7 +208,7 @@ bodeplot @recipe function bodeplot(p::Bodeplot; plotphase=true) systems = p.args[1] if !isa(systems, AbstractArray) - systems = [systems] + systems = typeof(systems)[systems] end if length(p.args) >= 2 w = p.args[2] @@ -220,7 +220,7 @@ bodeplot end ny, nu = size(systems[1]) s2i(i,j) = sub2ind((nu,(plotphase?2:1)*ny),j,i) - layout := ((plotphase?2:1)*ny,nu) + layout --> ((plotphase?2:1)*ny,nu) nw = length(w) xticks --> getLogTicks(w) @@ -250,7 +250,7 @@ bodeplot end xguide --> xlab yguide --> "Magnitude $_PlotScaleStr" - subplot := s2i((plotphase?(2i-1):i),j) + subplot --> s2i((plotphase?(2i-1):i),j) title --> "Bode plot from: u($j)" label --> "\$G_\{$(si)\}\$" linestyle --> styledict[:l] @@ -264,7 +264,7 @@ bodeplot grid --> true xscale --> :log10 yguide --> "Phase (deg)" - subplot := s2i(2i,j) + subplot --> s2i(2i,j) xguide --> "Frequency (rad/s)" label --> "\$G_\{$(si)\}\$" linestyle --> styledict[:l] @@ -339,7 +339,7 @@ nyquistplot ny, nu = size(systems[1]) w = length(p.args) < 2 ? _default_freq_vector(systems, :nyquist) : p.args[2] nw = length(w) - layout := (ny,nu) + layout --> (ny,nu) s2i(i,j) = sub2ind((ny,nu),j,i) # Ensure that `axes` is always a matrix of handles for (si,s) = enumerate(systems) @@ -353,7 +353,7 @@ nyquistplot xlims := (min(max(-20,minimum(redata)),-1), max(min(20,maximum(redata)),1)) title --> "Nyquist plot from: u($j)" yguide --> "To: y($i)" - subplot := s2i(i,j) + subplot --> s2i(i,j) label --> "\$G_\{$(si)\}\$" styledict = getStyleSys(si,length(systems)) linestyle --> styledict[:l] @@ -439,7 +439,7 @@ nicholsplot nw = length(w) # Gain circle functions - angle(x) = unwrap(atan2(imag.(x),real.(x))) + angle(x) = unwrap(atan2.(imag.(x),real.(x))) RadM(m) = @. abs(m/(m^2-1)) CentreM(m) = @. m^2/(1-m^2) Ny(mdb,t) = @. CentreM(10^(mdb/20))+RadM(10^(mdb/20)).*(cosd(t)+im.*sind(t)) @@ -729,12 +729,13 @@ Gang-of-Four plot. `kwargs` is sent as argument to Plots.plot.""" -> function gangoffourplot(P::Union{Vector, LTISystem}, C::Vector, args...; plotphase=false, kwargs...) - S,D,N,T = gangoffour(P,C) - fig = bodeplot(LTISystem[[S[i] D[i]; N[i] T[i]] for i = 1:length(C)], args..., plotphase=plotphase; kwargs...) + # Array of (S,D,N,T) + sys = gangoffour.(P,C) + fig = bodeplot([[sys[i][1] sys[i][2]; sys[i][3] sys[i][4]] for i = 1:length(C)], args..., plotphase=plotphase; kwargs...) titles = fill("", 1, plotphase ? 8 : 4) # Empty titles on phase titleIdx = plotphase ? [1,2,5,6] : [1,2,3,4] - titles[titleIdx] = ["\$S = 1/(1+PC)\$", "\$D = P/(1+PC)\$", "\$N = C/(1+PC)\$", "\$T = PC/(1+PC\$)"] + titles[titleIdx] = ["S = 1/(1+PC)", "D = P/(1+PC)", "N = C/(1+PC)", "T = PC/(1+PC)"] Plots.plot!(fig, title = titles) return fig end diff --git a/src/simplification.jl b/src/simplification.jl index 68e395de5..6e11612be 100644 --- a/src/simplification.jl +++ b/src/simplification.jl @@ -6,8 +6,7 @@ determined to be uncontrollable and unobservable based on the location of 0s in `sys` are removed.""" -> function sminreal(sys::StateSpace) A, B, C, inds = struct_ctrb_obsv(sys) - return StateSpace(A, B, C, sys.D, sys.Ts, sys.statenames[inds], - sys.inputnames, sys.outputnames) + return StateSpace(A, B, C, sys.D, sys.Ts) end # Determine the structurally controllable and observable realization for the system diff --git a/src/simulators.jl b/src/simulators.jl new file mode 100644 index 000000000..a20257623 --- /dev/null +++ b/src/simulators.jl @@ -0,0 +1,56 @@ +abstract type AbstractSimulator end + +# ============================================================================================ + +""" + Simulator +# Fields: + P::StateSpace + f = (x,p,t) -> x + y = (x,t) -> y +""" +struct Simulator <: AbstractSimulator + P::StateSpace + f + y +end + + +""" + Simulator(P::StateSpace, u = (x,t) -> 0) + +Used to simulate continuous-time systems. See function `?solve` for additional info. +# Usage: +``` +using OrdinaryDiffEq +h = 0.1 +Tf = 20 +t = 0:h:Tf +P = ss(tf(1,[2,1])^2) +K = 5 +reference(x,t) = [1.] +s = Simulator(P, reference) +x0 = [0.,0] +tspan  = (0.0,Tf) +sol = solve(s, x0, tspan, Tsit5()) +plot(t, s.y(sol, t)[:], lab="Open loop step response") +``` +""" +function Simulator(P::StateSpace, u = (x,t) -> 0) + @assert iscontinuous(P) "Simulator only supports continuous-time system. See function `lsim` for simulation of discrete-time systems." + @assert all(P.D .== 0) "Can not simulate systems with direct term D != 0" + f = (dx,x,p,t) -> dx .= P.A*x .+ P.B*u(x,t) + y(x,t) = P.C*x #.+ P.D*u(x,t) + y(sol::ODESolution,t) = P.C*sol(t) + Simulator(P, f, y) +end + + +""" + sol = solve(s::AbstractSimulator, x0, tspan, args...; kwargs...) +Simulate the system represented by `s` from initial state `x0` over time span `tspan = (t0,tf)`. +`args` and `kwargs` are sent to the `solve` function from `OrdinaryDiffEq`, e.g., `solve(s, x0, tspan, Tsit5(), reltol=1e-5)` solves the problem with solver [`Tsit5()`](http://docs.juliadiffeq.org/stable/solvers/ode_solve.html) and relative tolerance 1e-5. + +See also `Simulator` `lsim` +""" +DiffEqBase.solve(s::AbstractSimulator, x0, tspan, solver=Tsit5(), args...; kwargs...) = solve(ODEProblem(s.f,x0,tspan), solver, args...; kwargs...) diff --git a/src/synthesis.jl b/src/synthesis.jl index 02c1a23db..1ad1be9ba 100644 --- a/src/synthesis.jl +++ b/src/synthesis.jl @@ -143,10 +143,10 @@ end """ `feedback(L)` Returns L/(1+L) -`feedback(P,C)` Returns PC/(1+PC) +`feedback(P1,P2)` Returns P1/(1+P1*P2) """ feedback(L::TransferFunction) = L/(1+L) -feedback(P::TransferFunction, C::TransferFunction) = feedback(P*C) +feedback(P1::TransferFunction, P2::TransferFunction) = P1/(1+P1*P2) #Efficient implementations function feedback{T<:SisoRational}(L::TransferFunction{T}) @@ -156,19 +156,19 @@ function feedback{T<:SisoRational}(L::TransferFunction{T}) P = numpoly(L) Q = denpoly(L) #Extract polynomials and create P/(P+Q) - tf(P[1][:],(P+Q)[1][:], Ts=L.Ts) + tf(P[1][:],(P+Q)[1][:], L.Ts) end -function ControlSystems.feedback{T<:ControlSystems.SisoZpk}(L::TransferFunction{T}) +function feedback(L::TransferFunction{T}) where {T<:SisoZpk} if size(L) != (1,1) error("MIMO TransferFunction inversion isn't implemented yet") end numer = num(L.matrix[1]) k = L.matrix[1].k - denpol = k*prod(numpoly(L)[1])+prod(denpoly(L)[1]) + denpol = k*prod(numpoly(L)[1])+prod(denpoly(L)[1]) # TODO: Chcek indexing into polynomials kden = denpol[1] #Extract polynomials and create P/(P+Q) - zpk(numer,ControlSystems.roots(denpol), k/kden, Ts=L.Ts) + zpk(numer, roots(denpol), k/kden, L.Ts) end """ diff --git a/src/timeresp.jl b/src/timeresp.jl index 976a140c9..65a04e8bd 100644 --- a/src/timeresp.jl +++ b/src/timeresp.jl @@ -2,80 +2,73 @@ # XXX : `step` is a function in Base, with a different meaning than it has # here. This shouldn't be an issue, but it might be. -@doc """`[y, t, x] = step(sys[, Tf])` or `[y, t, x] = step(sys[, t])` +@doc """`y, t, x = step(sys[, Tf])` or `y, t, x = step(sys[, t])` Calculate the step response of system `sys`. If the final time `Tf` or time vector `t` is not provided, one is calculated based on the system pole locations. `y` has size `(length(t), ny, nu)`, `x` has size `(length(t), nx, nu)`""" -> -function Base.step(sys::StateSpace, t::AbstractVector) +function Base.step(sys::StateSpace, t::AbstractVector; method=:cont) lt = length(t) ny, nu = size(sys) nx = sys.nx - u = ones(size(t)) - x0 = zeros(nx, 1) + u = (x,t)->[one(eltype(t))] + x0 = zeros(nx) if nu == 1 - y, t, x = lsim(sys, u, t, x0=x0, method=:zoh) + y, tout, x, _ = lsim(sys, u, t, x0=x0, method=method) else x = Array{Float64}(lt, nx, nu) y = Array{Float64}(lt, ny, nu) for i=1:nu - y[:,:,i], t, x[:,:,i] = lsim(sys[:,i], u, t, x0=x0, method=:zoh) + y[:,:,i], tout, x[:,:,i],_ = lsim(sys[:,i], u, t, x0=x0, method=method) end end return y, t, x end -function Base.step(sys::TransferFunction{SisoGeneralized}, t::AbstractVector) - lsim(sys, ones(length(t), sys.nu), t) -end -Base.step(sys::LTISystem, Tf::Real) = step(sys, _default_time_vector(sys, Tf)) -Base.step(sys::LTISystem) = step(sys, _default_time_vector(sys)) -Base.step(sys::TransferFunction, t::AbstractVector) = step(ss(sys), t::AbstractVector) +Base.step(sys::LTISystem, Tf::Real; kwargs...) = step(sys, _default_time_vector(sys, Tf); kwargs...) +Base.step(sys::LTISystem; kwargs...) = step(sys, _default_time_vector(sys); kwargs...) +Base.step(sys::TransferFunction, t::AbstractVector; kwargs...) = step(ss(sys), t::AbstractVector; kwargs...) -@doc """`[y, t, x] = impulse(sys[, Tf])` or `[y, t, x] = impulse(sys[, t])` +@doc """`y, t, x = impulse(sys[, Tf])` or `y, t, x = impulse(sys[, t])` Calculate the impulse response of system `sys`. If the final time `Tf` or time vector `t` is not provided, one is calculated based on the system pole locations. `y` has size `(length(t), ny, nu)`, `x` has size `(length(t), nx, nu)`""" -> -function impulse(sys::StateSpace, t::AbstractVector) +function impulse(sys::StateSpace, t::AbstractVector; method=:cont) + T = promote_type(eltype(sys.A), Float64) lt = length(t) ny, nu = size(sys) nx = sys.nx - u = zeros(size(t)) - if iscontinuous(sys) + if iscontinuous(sys) #&& method == :cont + u = (x,i) -> [zero(T)] # impulse response equivalent to unforced response of # ss(A, 0, C, 0) with x0 = B. - imp_sys = ss(sys.A, zeros(nx, 1), sys.C, zeros(ny, 1)) + imp_sys = ss(sys.A, zeros(T, nx, 1), sys.C, zeros(T, ny, 1)) x0s = sys.B else + u = (x,i) -> i == t[1] ? [one(T)]/sys.Ts : [zero(T)] imp_sys = sys - x0s = zeros(nx, nu) - u[1] = 1/sys.Ts + x0s = zeros(T, nx, nu) end - if nu == 1 - y, t, x = lsim(sys, u, t, x0=x0s, method=:zoh) + if nu == 1 # Why two cases # QUESTION: Not type stable? + y, t, x,_ = lsim(sys, u, t, x0=x0s[:], method=method) else - x = Array{Float64}(lt, nx, nu) - y = Array{Float64}(lt, ny, nu) + x = Array{T}(lt, nx, nu) + y = Array{T}(lt, ny, nu) for i=1:nu - y[:,:,i], t, x[:,:,i] = lsim(sys[:,i], u, t, x0=x0s[:,i], method=:zoh) + y[:,:,i], t, x[:,:,i],_ = lsim(sys[:,i], u, t, x0=x0s[:,i], method=method) end end return y, t, x end -function impulse(sys::TransferFunction{SisoGeneralized}, t::AbstractVector) - u = zeros(length(t), sys.nu) - u[1,:] = 1/(t[2]-t[1]) - lsim(sys::TransferFunction{SisoGeneralized}, u, t) -end -impulse(sys::LTISystem, Tf::Real) = impulse(sys, _default_time_vector(sys, Tf)) -impulse(sys::LTISystem) = impulse(sys, _default_time_vector(sys)) -impulse(sys::TransferFunction, t::AbstractVector) = impulse(ss(sys), t) +impulse(sys::LTISystem, Tf::Real; kwags...) = impulse(sys, _default_time_vector(sys, Tf); kwags...) +impulse(sys::LTISystem; kwags...) = impulse(sys, _default_time_vector(sys); kwags...) +impulse(sys::TransferFunction, t::AbstractVector; kwags...) = impulse(ss(sys), t; kwags...) @doc """`y, t, x = lsim(sys, u, t; x0, method])` @@ -86,13 +79,11 @@ a zero vector is used. `y`, `x`, `uout` has time in the first dimension. Initial state `x0` defaults to zero. -Continuous time systems are discretized before simulation. By default, the -method is chosen based on the smoothness of the input signal. Optionally, the -`method` parameter can be specified as either `:zoh` or `:foh`, default depends on system smoothnes. +Continuous time systems are simulated using an ODE solver if `u` is a function. If `u` is an array, the system is discretized before simulation. For a lower level inteface, see `?Simulator` and `?solve` `u` can be a function or a matrix/vector of precalculated control signals. -If `u` is a function, then `u(i,x)` is called to calculate the control signal every iteration. This can be used to provide a control law such as state feedback `u(t,x) = -L*x` calculated by `lqr`. -To simulate a unit step, use `(i,x)-> 1`, for a ramp, use `(i,x)-> i*h`, for a step at `t=5`, use (i,x)-> (i*h >= 5) etc. +If `u` is a function, then `u(x,i)` (`u(x,t)`) is called to calculate the control signal every iteration (time instance used by solver). This can be used to provide a control law such as state feedback `u(x,t) = -L*x` calculated by `lqr`. +To simulate a unit step, use `(x,i)-> 1`, for a ramp, use `(x,i)-> i*h`, for a step at `t=5`, use (x,i)-> (i*h >= 5) etc. Usage example: ```julia @@ -104,27 +95,28 @@ Q = eye(2) R = eye(1) L = lqr(sys,Q,R) -u(t,x) = -L*x # Form control law, +u(x,t) = -L*x # Form control law, t=0:0.1:5 x0 = [1,0] y, t, x, uout = lsim(sys,u,t,x0) plot(t,x, lab=["Position", "Velocity"]', xlabel="Time [s]") ``` -""" -> +""" function lsim(sys::StateSpace, u::AbstractVecOrMat, t::AbstractVector; - x0::VecOrMat=zeros(sys.nx, 1), method::Symbol=_issmooth(u) ? :foh : :zoh) + x0::VecOrMat=zeros(sys.nx), method::Symbol=_issmooth(u) ? :foh : :zoh) ny, nu = size(sys) nx = sys.nx if length(x0) != nx error("size(x0) must match the number of states of sys") - elseif !(size(u) in [(length(t), nu) (length(t),)]) + end + if !(size(u) in [(length(t), nu) (length(t),)]) error("u must be of size (length(t), nu)") end dt = Float64(t[2] - t[1]) if !iscontinuous(sys) || method == :zoh - if iscontinuous(sys) + if iscontinuous(sys) # Looks strange to check iscontinuous again dsys = c2d(sys, dt, :zoh)[1] else if sys.Ts != dt @@ -136,7 +128,7 @@ function lsim(sys::StateSpace, u::AbstractVecOrMat, t::AbstractVector; dsys, x0map = c2d(sys, dt, :foh) x0 = x0map*[x0; u[1:1,:].'] end - x = ltitr(dsys.A, dsys.B, map(Float64,u), map(Float64,x0)) + x = ltitr(dsys.A, dsys.B, Float64.(u), Float64.(x0)) y = (sys.C*(x.') + sys.D*(u.')).' return y, t, x end @@ -145,16 +137,17 @@ end @deprecate lsim(sys, u, t, x0, method) lsim(sys, u, t; x0=x0, method=method) function lsim(sys::StateSpace, u::Function, t::AbstractVector; - x0::VecOrMat=zeros(sys.nx, 1), method::Symbol=:zoh) + x0::VecOrMat=zeros(sys.nx), method::Symbol=:cont) ny, nu = size(sys) nx = sys.nx if length(x0) != nx error("size(x0) must match the number of states of sys") - elseif size(u(1,x0)) != (nu,) && size(u(1,x0)) != (nu,1) + elseif size(u(x0,1)) != (nu,) && size(u(x0,1)) != (nu,1) error("return value of u must be of size nu") end + T = promote_type(Float64, eltype(x0)) - dt = Float64(t[2] - t[1]) + dt = T(t[2] - t[1]) if !iscontinuous(sys) || method == :zoh if iscontinuous(sys) dsys = c2d(sys, dt, :zoh)[1] @@ -164,32 +157,23 @@ function lsim(sys::StateSpace, u::Function, t::AbstractVector; end dsys = sys end + x,uout = ltitr(dsys.A, dsys.B, u, t, T.(x0)) else - dsys, x0map = c2d(sys, dt, :foh) - x0 = x0map*[x0; u(1,x0)] + s = Simulator(sys, u) + sol = solve(s, T.(x0), (t[1],t[end]), Tsit5()) + x = sol(t)' + uout = Array{eltype(x)}(length(t), ninputs(sys)) + for (i,ti) in enumerate(t) + uout[i,:] = u(x[i,:],ti)' + end end - x,uout = ltitr(dsys.A, dsys.B, u, length(t), map(Float64,x0)) y = (sys.C*(x.') + sys.D*(uout.')).' return y, t, x, uout end + lsim(sys::TransferFunction, u, t, args...; kwargs...) = lsim(ss(sys), u, t, args...; kwargs...) -function lsim(sys::TransferFunction{SisoGeneralized}, u, t) - ny, nu = size(sys) - if !any(size(u) .== [(length(t), nu) (length(t),)]) - error("u must be of size (length(t), nu)") - end - ny, nu = size(sys) - y = Array{Float64}(length(t),ny) - for i = 1:nu - for o = 1:ny - dt = Float64(t[2]-t[1]) - y[:,o] += lsimabstract(sys.matrix[o,i], u[:,i], dt, t[end]) - end - end - return y, t -end @doc """`ltitr(A, B, u[,x0])` @@ -198,10 +182,10 @@ end Simulate the discrete time system `x[k + 1] = A x[k] + B u[k]`, returning `x`. If `x0` is not provided, a zero-vector is used. -If `u` is a function, then `u(i,x)` is called to calculate the control signal every iteration. This can be used to provide a control law such as state feedback `u=-Lx` calculated by `lqr`. In this case, an integrer `iters` must be provided that indicates the number of iterations. +If `u` is a function, then `u(x,i)` is called to calculate the control signal every iteration. This can be used to provide a control law such as state feedback `u=-Lx` calculated by `lqr`. In this case, an integrer `iters` must be provided that indicates the number of iterations. """ -> function ltitr{T}(A::Matrix{T}, B::Matrix{T}, u::AbstractVecOrMat{T}, - x0::VecOrMat{T}=zeros(T, size(A, 1), 1)) + x0::VecOrMat{T}=zeros(T, size(A, 1))) n = size(u, 1) x = Array{T}(size(A, 1), n) for i=1:n @@ -212,14 +196,15 @@ function ltitr{T}(A::Matrix{T}, B::Matrix{T}, u::AbstractVecOrMat{T}, end -function ltitr{T}(A::Matrix{T}, B::Matrix{T}, u::Function, iters::Int, - x0::VecOrMat{T}=zeros(T, size(A, 1), 1)) +function ltitr{T}(A::Matrix{T}, B::Matrix{T}, u::Function, t, + x0::VecOrMat{T}=zeros(T, size(A, 1))) + iters = length(t) x = Array{T}(size(A, 1), iters) uout = Array{T}(size(B, 2), iters) for i=1:iters x[:,i] = x0 - uout[:,i] = u(i,x0) + uout[:,i] = u(x0,t[i]) x0 = A * x0 + B * uout[:,i] end return x.', uout.' @@ -253,7 +238,6 @@ function _default_Ts(sys::LTISystem) return Ts end -_default_Ts(sys::TransferFunction{SisoGeneralized}) = 0.07 #TODO a reasonable check _issmooth(u::Function) = false diff --git a/src/types/Lti.jl b/src/types/Lti.jl new file mode 100644 index 000000000..61112ede7 --- /dev/null +++ b/src/types/Lti.jl @@ -0,0 +1,53 @@ +abstract type LTISystem <: AbstractSystem end ++(sys1::LTISystem, sys2::LTISystem) = +(promote(sys1, sys2)...) +-(sys1::LTISystem, sys2::LTISystem) = -(promote(sys1, sys2)...) +*(sys1::LTISystem, sys2::LTISystem) = *(promote(sys1, sys2)...) +/(sys1::LTISystem, sys2::LTISystem) = /(promote(sys1, sys2)...) + + +@doc """`issiso(sys)` + +Returns `true` if `sys` is SISO, else returns `false`.""" -> +function issiso(sys::LTISystem) + return ninputs(sys) == 1 && noutputs(sys) == 1 +end + + +@doc """`iscontinuous(sys)` + +Returns `true` if `sys` is continuous, else returns `false`.""" -> +function iscontinuous(sys::LTISystem) + return sys.Ts == 0 +end + + +@doc """`isstable(sys)` + +Returns `true` if `sys` is stable, else returns `false`.""" -> +function isstable(sys::LTISystem) + if iscontinuous(sys) + if any(real.(pole(sys)).>=0) + return false + end + else + if any(abs.(pole(sys)).>=1) + return false + end + end + return true +end + + + + + +function _check_consistent_sampling_time(systems::AbstractVector{LTISystem}) + if !all(s.Ts == Ts for s in systems) + error("Sampling time mismatch") + end +end +function _check_consistent_sampling_time(sys1::LTISystem, sys2::LTISystem) + if sys1.Ts != sys2.Ts + error("Sampling time mismatch") + end +end diff --git a/src/types/SisoTf.jl b/src/types/SisoTf.jl new file mode 100644 index 000000000..202b79493 --- /dev/null +++ b/src/types/SisoTf.jl @@ -0,0 +1,7 @@ +abstract type SisoTf{T<:Number} end + ++(f1::SisoTf, f2::SisoTf) = +(promote(f1,f2)...) +-(f1::SisoTf, f2::SisoTf) = -(promote(f1,f2)...) + +*(f1::SisoTf, f2::SisoTf) = *(promote(f1,f2)...) +/(f1::SisoTf, f2::SisoTf) = /(promote(f1,f2)...) diff --git a/src/types/SisoTfTypes/SisoRational.jl b/src/types/SisoTfTypes/SisoRational.jl new file mode 100644 index 000000000..7335fd69e --- /dev/null +++ b/src/types/SisoTfTypes/SisoRational.jl @@ -0,0 +1,120 @@ +## User should just use TransferFunction +struct SisoRational{T} <: SisoTf{T} + num::Poly{T} + den::Poly{T} + function SisoRational{T}(num::Poly{T}, den::Poly{T}) where T <: Number + if all(den == zero(den)) + error("Cannot create SisoRational with zero denominator") + elseif all(num == zero(num)) + # The numerator is zero, make the denominator 1 + den = one(den) + end + new{T}(num, den) + end +end +function SisoRational(num::Poly{T1}, den::Poly{T2}) where T1 <: Number where T2 <: Number + T = promote_type(T1,T2) + SisoRational{T}(Poly{T}(num.a), Poly{T}(den.a)) +end +function SisoRational{T}(num::AbstractVector, den::AbstractVector) where T <: Number # NOTE: Typearguemnts on the parameters? + SisoRational{T}(Poly{T}(reverse(num)), Poly{T}(reverse(den))) +end +function SisoRational(num::AbstractVector{T1}, den::AbstractVector{T2}) where T1 <: Number where T2 <: Number + T = promote_type(T1,T2) + SisoRational{T}(num, den) +end +# NOTE: How many of these above are actually needed? +# TODO: Add method for scalar inputs + + +Base.zero(::Type{SisoRational{T}}) where T = SisoRational{T}([zero(T)], [one(T)]) +Base.one(::Type{SisoRational{T}}) where T = SisoRational{T}([one(T)], [one(T)]) + +Base.one(f::SisoRational) = one(typeof(f)) +Base.zero(f::SisoRational) = zero(typeof(f)) + +isproper(f::SisoRational) = (length(f.num) <= length(f.den)) + +function minreal(sys::SisoRational, eps::Real=sqrt(eps())) + return SisoRational(minreal(SisoZpk(sys), eps)) +end + +function print_siso(io::IO, f::SisoRational, var=:s) + # Convert the numerator and denominator to strings + numstr = sprint(printpolyfun(var), f.num) + denstr = sprint(printpolyfun(var), f.den) + + # Figure out the length of the separating line + len_num = length(numstr) + len_den = length(denstr) + dashcount = max(len_num, len_den) + + # Center the numerator or denominator + if len_num < dashcount + numstr = "$(repeat(" ", div(dashcount - len_num, 2)))$numstr" + else + denstr = "$(repeat(" ", div(dashcount - len_den, 2)))$denstr" + end + println(io, numstr) + println(io, repeat("-", dashcount)) + println(io, denstr) +end + +# TODO is this working? +function print_compact(io::Base.IO, f::SisoRational, var) + numstr = sprint(print_poly, f.num) + denstr = sprint(print_poly, f.den) + println(io, "($numstr)/($denstr)") +end + + +Base.num(f::SisoRational) = reverse(coeffs(f.num)) +Base.den(f::SisoRational) = reverse(coeffs(f.den)) + +denpoly(f::SisoRational) = f.den +numpoly(f::SisoRational) = f.num + +tzero(f::SisoRational) = roots(f.num) +pole(f::SisoRational) = roots(f.den) + +function evalfr(f::SisoRational{T}, s::Number) where T + S = promote_op(/, promote_type(T, typeof(s)), promote_type(T, typeof(s))) + den = polyval(f.den, s) + if den == zero(S) + convert(S,Inf) + else + polyval(f.num, s)/den + end +end + +==(f1::SisoRational, f2::SisoRational) = (f1.num * f2.den[1] == f2.num * f1.den[1] && f1.den * f2.den[1] == f2.den * f1.den[1]) # NOTE: Not in analogy with how it's done for SisoZpk + +# We might want to consider alowing scaled num and den as equal +function isapprox(f1::SisoRational, f2::SisoRational; rtol::Real=sqrt(eps()), atol::Real=0) + isapprox(f1.num * f2.den[1], f2.num * f1.den[1], rtol=rtol, atol=atol) && isapprox(f1.den * f2.den[1], f2.den * f1.den[1], rtol=rtol, atol=atol) +end + ++(f1::SisoRational, f2::SisoRational) = SisoRational(f1.num*f2.den + f2.num*f1.den, f1.den*f2.den) ++(f::SisoRational, n::Number) = SisoRational(f.num + n*f.den, f.den) ++(n::Number, f::SisoRational) = f + n +#.+(f::SisoRational, n::Number) = t + n +#.+(n::Number, f::SisoRational) = t + n + +-(f1::SisoRational, f2::SisoRational) = SisoRational(f1.num*f2.den - f2.num*f1.den, f1.den*f2.den) +-(n::Number, f::SisoRational) = SisoRational(n*f.den - f.num, f.den) +-(f::SisoRational, n::Number) = +(f, -n) +#.-(f::SisoRational, n::Number) = -(t, n) +#.-(n::Number, f::SisoRational) = -(n, t) + +-(f::SisoRational) = SisoRational(-f.num, f.den) + +*(f1::SisoRational, f2::SisoRational) = SisoRational(f1.num*f2.num, f1.den*f2.den) +*(f::SisoRational, n::Number) = SisoRational(f.num*n, f.den) +*(n::Number, f::SisoRational) = *(f, n) +#.*(f1::SisoRational, f2::SisoRational) = *(f1, f2) +#.*(f::SisoRational, n::Number) = *(f, n) +#.*(n::Number, f::SisoRational) = *(f, n) + +/(n::Number, f::SisoRational) = SisoRational(n*f.den, f.num) +/(f::SisoRational, n::Number) = f*(1/n) +/(f1::SisoRational, f2::SisoRational) = f1*(1/f2) diff --git a/src/types/SisoTfTypes/SisoZpk.jl b/src/types/SisoTfTypes/SisoZpk.jl new file mode 100644 index 000000000..388597ae5 --- /dev/null +++ b/src/types/SisoTfTypes/SisoZpk.jl @@ -0,0 +1,264 @@ +# FIXME: Add some condition to guarantee that: TR <: Union(T, complex(T)) ! + +# T the numeric type of the transfer function +# TR the type of the roots +struct SisoZpk{T,TR<:Number} <: SisoTf{T} + z::Vector{TR} + p::Vector{TR} + k::T + function SisoZpk{T,TR}(z::Vector{TR}, p::Vector{TR}, k::T) where {T<:Number, TR<:Number} + if k == zero(T) + p = TR[] + z = TR[] + end + if TR <: Complex && T <: Real + @assert check_real(z) "zpk model should be real-valued, but zeros do not come in conjugate pairs." + @assert check_real(p) "zpk model should be real-valued, but poles do not come in conjugate pairs." + end + new{T,TR}(z, p, k) + end +end +function SisoZpk{T,TR}(z::Vector, p::Vector, k::Number) where {T<:Number, TR<:Number} + SisoZpk{T,TR}(Vector{TR}(z), Vector{TR}(p), T(k)) +end +function SisoZpk{T}(z::Vector, p::Vector, k::Number) where T + TR = complex(T) + SisoZpk{T,TR}(Vector{TR}(z), Vector{TR}(p), T(k)) +end +function SisoZpk(z::AbstractVector{TZ}, p::AbstractVector{TP}, k::T) where {T<:Number, TZ<:Number, TP<:Number} # NOTE: is this constructor really needed? + TR = promote_type(TZ,TP) + # Could check if complex roots come with their conjugates, + # i.e., if the SisoZpk corresponds to a real-valued system + + if TR <: Complex && T <: Real + @assert check_real(z) "zpk model should be real-valued, but zeros do not come in conjugate pairs." + @assert check_real(p) "zpk model should be real-valued, but poles do not come in conjugate pairs." + end + SisoZpk{T,TR}(Vector{TR}(z), Vector{TR}(p), k) +end + + + + +Base.zero(::Type{SisoZpk{T}}) where T = SisoZpk{T}(T[], T[], zero(T)) +Base.one(::Type{SisoZpk{T}}) where T = SisoZpk{T}(T[], T[], one(T)) + +Base.zero(::Type{SisoZpk{T,TR}}) where {T, TR} = SisoZpk{T,TR}(TR[], TR[], zero(T)) +Base.one(::Type{SisoZpk{T,TR}}) where {T, TR} = SisoZpk{T,TR}(TR[], TR[], one(T)) + +Base.one(f::SisoZpk) = one(typeof(f)) +Base.zero(f::SisoZpk) = zero(typeof(f)) + + +# tzero is not meaningful for transfer function element? But both zero and zeros are taken... +tzero(f::SisoZpk) = f.z # Do minreal first?, +pole(f::SisoZpk) = f.p # Do minreal first? + +numpoly(f::SisoZpk{<:Real}) = f.k*prod(roots2real_poly_factors(f.z)) +denpoly(f::SisoZpk{<:Real}) = prod(roots2real_poly_factors(f.p)) + +numpoly(f::SisoZpk) = f.k*prod(roots2poly_factors(f.z)) +denpoly(f::SisoZpk) = prod(roots2poly_factors(f.p)) + +num(f::SisoZpk) = reverse(coeffs(numpoly(f))) # FIXME: reverse?! +den(f::SisoZpk) = reverse(coeffs(denpoly(f))) # FIXME: reverse?! + +isproper(f::SisoZpk) = (length(f.z) <= length(f.p)) + +# Will not be able to create complex polynomials without type instability + + +function minreal(sys::SisoZpk{T,TR}, eps::Real) where {T, TR} + if length(sys.p) == 0 + return sys + end + + newZ = copy(sys.z) + newP = Vector{TR}(0) + + pidx = 1 + p = sys.p[pidx] + while true + if isempty(newZ) + push!(newP, p) + else + distance, zidx = findmin(abs.(p-newZ)) + + if distance < eps + if imag(p) == 0 && imag(newZ[zidx]) != 0 + newZ[zidx+1] = real(newZ[zidx+1]) + end + if imag(newZ[zidx]) == 0 && imag(p) != 0 + pidx += 1 + p = real(sys.p[pidx]) + deleteat!(newZ, zidx) + continue + end + deleteat!(newZ, zidx) + else + push!(newP, p) + end + end + + pidx += 1 + if pidx > length(sys.p) + break + end + p = sys.p[pidx] + end + SisoZpk{T, TR}(newZ, newP, sys.k) +end + +# FIXME: Perhaps move to some other file with auxilliary functions, +# the name could also be imporoved. Perhaps this functionality can be found in some other package. +""" If TR is Complex and T is Real, check that every pole is matched to its conjugate + this assumes that the compelx poles are ordered as they are output by the LAPACK + routines that return complex-conjugated values, i.e., (x+iy) is followed by (x-iy)""" +function check_real(r_vec::AbstractVector{<:Complex}) + k = 1 + while k <= length(r_vec) + if isreal(r_vec[k]) + # Do nothing + else + if k == length(r_vec) # Element is complex and there is no more elements to match up with + return false + elseif conj(r_vec[k]) == r_vec[k+1] # Element is complex and succeeding element is its conjugate + k=k+1 + else + return false + end + end + k += 1 + end + return true +end + + + +function evalfr(f::SisoZpk{T1,TR}, s::T2) where {T1<:Number, TR<:Number, T2<:Number} + T0 = promote_type(T2, TR) + T = promote_type(T1, Base.promote_op(/, T0, T0)) + den = reduce(*, one(T0), s-a for a in f.p) # Default to one + if den == zero(T0) + return convert(T, Inf) + else + num = reduce(*, one(T0), s-b for b in f.z) + return f.k*num/den + end +end + + +function poly_factors2string(poly_factors::AbstractArray{<:Poly{T}}, var) where T + if length(poly_factors) == 0 + str = sprint(printpolyfun(var), Poly{T}(one(T))) + elseif length(poly_factors) == 1 + str = sprint(printpolyfun(var), poly_factors[1]) + else + str = reduce(*,"",["("*sprint(printpolyfun(var), z)*")" for z in poly_factors]) + end +end + +""" Heurisitc function that tries to add parentheses when printing the coeffcient + for systems on zpk form. Should at least handle the following types + Measurment, Dual, Sym. """ +function _printcoefficient(nbr::Number) + # Print type information as in 1.0f0 for Float32 + # showcompact might be better, but is not consistent with polynomials + nbr_string = sprint(show,nbr) + if contains(nbr_string, " + ") || contains(nbr_string, " - ") || contains(nbr_string, " ± ") + return "(" * nbr_string * ")" # Add parens + else + return nbr_string + end +end + +function print_siso(io::IO, t::SisoZpk, var=:s) + if typeof(t.k) <: Real + numstr = poly_factors2string(roots2real_poly_factors(t.z), var) + denstr = poly_factors2string(roots2real_poly_factors(t.p), var) + else + numstr = poly_factors2string(roots2poly_factors(t.z), var) + denstr = poly_factors2string(roots2poly_factors(t.p), var) + end + # Figure out the length of the separating line + len_num = length(numstr) + len_den = length(denstr) + dashcount = max(len_num, len_den) + + # Center the numerator or denominator + if len_num < dashcount + numstr = "$(repeat(" ", div(dashcount - len_num, 2)))$numstr" + else + denstr = "$(repeat(" ", div(dashcount - len_den, 2)))$denstr" + end + + gainstr = _printcoefficient(t.k) + #Add spaces to account for gain string + numstr = " "^(length(gainstr))*numstr + denstr = " "^(length(gainstr))*denstr + println(io, numstr) + println(io, gainstr*repeat("-", dashcount)) + println(io, denstr) +end + + + + + +==(f1::SisoZpk, f2::SisoZpk) = (f1-f2).k == 0 +function isapprox(f1::SisoZpk, f2::SisoZpk; rtol::Real=sqrt(eps()), atol::Real=sqrt(eps())) + fdiff = f1 - f2 + isapprox(fdiff.k, 0, atol=atol, rtol=rtol) +end + +function +(f1::SisoZpk{T1,TR1}, f2::SisoZpk{T2,TR2}) where {T1<:Number,T2<:Number,TR1<:Number,TR2<:Number} + numPoly = numpoly(f1)*denpoly(f2) + numpoly(f2)*denpoly(f1) + + TRtmp = promote_type(TR1, TR2) + # Calculating roots can make integers floats + TRnew = Base.promote_op(/, TRtmp, TRtmp) + z = convert(Vector{TRnew}, roots(numPoly)) + #TODO gains could remain integers, but numerical precision inhibits this + Ttmp = promote_type(T1,T2) + Tnew = Base.promote_op(/, Ttmp, Ttmp) + if length(numPoly) > 0 + k = convert(Tnew, numPoly[end]) + p = convert(Vector{TRnew}, [f1.p;f2.p]) + else + k = zero(Tnew) + p = TRnew[] + end + + # FIXME: + # Threshold for pole-zero cancellation should depend on the roots of the system + # Note the difference between continuous and discrete-time systems... + minreal(SisoZpk(z,p,k), sqrt(eps())) +end + + ++(f::SisoZpk, n::T) where {T<:Number} = f + SisoZpk{T,T}(T[],T[],n) ++(n::Number, f::SisoZpk) = f + n +#.+(t::SisoZpk, n::Number) = t + n +#.+(n::Number, t::SisoZpk) = t + n + +-(t1::SisoZpk, t2::SisoZpk) = +(t1,-t2) +-(n::T, f::SisoZpk) where {T<:Number} = SisoZpk{T,T}(T[],T[],n) - f +-(t::SisoZpk, n::Number) = +(t, -n) +#.-(t::SisoZpk, n::Number) = -(t, n) +#.-(n::Number, t::SisoZpk) = -(n, t) + + +-(f::SisoZpk) = SisoZpk(f.z, f.p, -f.k) + + +*(f1::SisoZpk, f2::SisoZpk) = SisoZpk([f1.z;f2.z], [f1.p;f2.p], f1.k*f2.k) + +*(f::SisoZpk, n::Number) = SisoZpk(f.z, f.p, f.k*n) +*(n::Number, f::SisoZpk) = *(f, n) +#.*(f1::SisoZpk, f2::SisoZpk) = *(f1, f2) +#.*(t::SisoZpk, n::Number) = *(t, n) +#.*(n::Number, t::SisoZpk) = *(t, n) + +/(n::Number, f::SisoZpk) = SisoZpk(f.p, f.z, n/f.k) +/(f::SisoZpk, n::Number) = SisoZpk(f.z, f.p, f.k/n) +/(f1::SisoZpk, f2::SisoZpk) = f1*(1/f2) diff --git a/src/types/SisoTfTypes/conversion.jl b/src/types/SisoTfTypes/conversion.jl new file mode 100644 index 000000000..921772719 --- /dev/null +++ b/src/types/SisoTfTypes/conversion.jl @@ -0,0 +1,49 @@ +function Base.convert(::Type{SisoZpk{T,TR}}, f::SisoRational{T2}) where {T<:Number, TR<:Number, T2<:Number} + if length(f.num) == 0 + return SisoZpk(T[],TR[],0) + elseif all(f.den == zero(f.den)) + error("Zero denominator, this should not be possible") # NOTE: Should we keep this? + end + + return SisoZpk{T,TR}(roots(f.num), roots(f.den), f.num[end]/f.den[end]) # NOTE: polynomials are stored with increasing coefficient orders +end + +function Base.convert(::Type{<:SisoZpk}, f::SisoRational{T}) where {T<:Number} + TR = complex(Base.promote_op(/, T, T)) # Type of roots(f.z) + Base.convert(SisoZpk{T,TR}, f) +end + + +function Base.convert(::Type{SisoZpk{T,TR}}, f::SisoZpk) where {T<:Number, TR<:Number} + return SisoZpk{T,TR}(f.z, f.p, f.k) +end + + +function Base.convert(::Type{SisoRational{T}}, f::SisoRational) where {T<:Number} + return SisoRational{T}(convert(Poly{T}, f.num), convert(Poly{T}, f.den)) +end + +function Base.convert(::Type{SisoRational{T1}}, f::SisoZpk{T2,TR}) where {T1<:Number, T2<:Number,TR<:Number} + if T2 <: Real + num = prod(roots2real_poly_factors(f.z))*f.k + den = prod(roots2real_poly_factors(f.p)) + else + num = prod(roots2poly_factors(f.z))*f.k + den = prod(roots2poly_factors(f.p)) + end + return SisoRational{T1}(num, den) +end + +function Base.convert(::Type{SisoRational}, f::SisoZpk{T1,TR}) where {T1<:Number, TR<:Number} + if T1 <: Real + T = promote_type(T1,typeof(real(one(TR)))) + # Alternative version if T1 rich enough + # T = T1 + convert(SisoRational{T}, f) + else + convert(SisoRational{promote_type(T1,TR)}, f) + end +end + + +Base.convert(::Type{S}, b::T2) where {T1, S <: SisoTf{T1}, T2<:Number} = convert(T1, b)*one(S) diff --git a/src/types/SisoTfTypes/polyprint.jl b/src/types/SisoTfTypes/polyprint.jl new file mode 100644 index 000000000..d50b37f6f --- /dev/null +++ b/src/types/SisoTfTypes/polyprint.jl @@ -0,0 +1,39 @@ +# Get show compatible function with var set +printpolyfun(var) = (io, p, mimetype = MIME"text/plain"()) -> printpolydesc(io, p, var, mimetype) + +# Specialized function, copied and adapted from Polynomials.jl +# Ignores var as set in Poly +""" Prints polynomial in descending order, with variable `var` +""" +function printpolydesc(io::IO, p::Poly{T}, var, mimetype = MIME"text/plain"()) where {T} + printpoly2(io, p, var, mimetype) +end + + +# Specialized function, copied and adapted from Polynomials.jl +# Ignores var as set in Poly +function printpoly2(io::IO, p::Poly{T}, var, mimetype) where {T} + first = true + printed_anything = false + for i in reverse(eachindex(p)) + printed = showterm2(io, p, i, first, var, mimetype) + first &= !printed + printed_anything |= printed + end + printed_anything || print(io, zero(T)) + return nothing +end + +#Specialized function, copied and adapted from Polynomials.jl +# Ignores var as set in Poly +function showterm2(io::IO, p::Poly{T}, j, first, var, mimetype) where {T} + pj = p[j] + + pj == zero(T) && return false + + pj = Polynomials.printsign(io, pj, j, first, mimetype) + Polynomials.printcoefficient(io, pj, j, mimetype) + Polynomials.printproductsign(io, pj, j, mimetype) + Polynomials.printexponent(io, var, j, mimetype) + true +end diff --git a/src/types/SisoTfTypes/promotion.jl b/src/types/SisoTfTypes/promotion.jl new file mode 100644 index 000000000..591174080 --- /dev/null +++ b/src/types/SisoTfTypes/promotion.jl @@ -0,0 +1,18 @@ +Base.promote_rule(::Type{SisoRational{T1}}, ::Type{SisoRational{T2}}) where {T1, T2} = SisoRational{promote_type(T1, T2)} +Base.promote_rule(::Type{SisoZpk{T1,TR1}}, ::Type{SisoZpk{T2,TR2}}) where {T1, TR1, T2, TR2} = SisoZpk{promote_type(T1, T2), promote_type(TR1,TR2)} + +function Base.promote_rule(::Type{SisoRational{T1}}, ::Type{SisoZpk{T2,TR2}}) where {T1, T2, TR2<:Number} + Tnew = promote_type(T1, T2) + TRnew = promote_type(TR2, complex(Tnew)) + SisoZpk{Tnew, TRnew} +end + + +Base.promote_rule(::Type{SisoRational{T1}}, ::Type{T2}) where {T1<:Number, T2<:Number} = SisoRational{promote_type(T1, T2)} + +function Base.promote_rule(::Type{SisoZpk{T1,TR1}}, ::Type{T2}) where {T1<:Number, TR1<:Number, T2<:Number} + Tnew = promote_type(T1, T2) + TRnew = promote_type(TR1, complex(Tnew)) + # TODO Not obvious that we want to promote to complex poles? + SisoZpk{Tnew, TRnew} +end diff --git a/src/types/StateSpace.jl b/src/types/StateSpace.jl new file mode 100644 index 000000000..ed9baba2a --- /dev/null +++ b/src/types/StateSpace.jl @@ -0,0 +1,251 @@ +##################################################################### +## Data Type Declarations ## +##################################################################### + +type StateSpace{T, MT<:AbstractMatrix{T}} <: LTISystem + A::MT + B::MT + C::MT + D::MT + Ts::Float64 + nx::Int + nu::Int + ny::Int + function StateSpace{T, MT}(A::MT, B::MT, + C::MT, D::MT, Ts::Float64) where {T, MT <: AbstractMatrix{T}} + nx = size(A, 1) + nu = size(B, 2) + ny = size(C, 1) + + if size(A, 2) != nx && nx != 0 + error("A must be square") + elseif size(B, 1) != nx + error("B must have the same row size as A") + elseif size(C, 2) != nx + error("C must have the same column size as A") + elseif nu != size(D, 2) + error("D must have the same column size as B") + elseif ny != size(D, 1) + error("D must have the same row size as C") + end + + # Validate sampling time + if Ts < 0 && Ts != -1 + error("Ts must be either a positive number, 0 + (continuous system), or -1 (unspecified)") + end + new{T, MT}(A, B, C, D, Ts, nx, nu, ny) + end +end +function StateSpace{T,MT}(A::AbstractArray, B::AbstractArray, C::AbstractArray, D::AbstractArray, Ts::Real) where {T, MT <: AbstractMatrix{T}} + return StateSpace{T,Matrix{T}}(MT(to_matrix(T, A)), MT(to_matrix(T, B)), MT(to_matrix(T, C)), + MT(to_matrix(T, D)), Float64(Ts)) +end +function StateSpace(A::AbstractArray, B::AbstractArray, C::AbstractArray, D::AbstractArray, Ts::Real) + # TODO: change back in 0.7 T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D)) + T = promote_type(promote_type(eltype(A),eltype(B)), promote_type(eltype(C),eltype(D))) + @assert (typeof(to_matrix(T, A)) == typeof(to_matrix(T, B)) == typeof(to_matrix(T, C)) == typeof(to_matrix(T, D))) + return StateSpace{T,Matrix{T}}(to_matrix(T, A), to_matrix(T, B), to_matrix(T, C), + to_matrix(T, D), Float64(Ts)) +end + +# Getter functions +get_A(sys::StateSpace) = sys.A +get_B(sys::StateSpace) = sys.B +get_C(sys::StateSpace) = sys.C +get_D(sys::StateSpace) = sys.D + +get_Ts(sys::StateSpace) = sys.Ts + +ssdata(sys::StateSpace) = get_A(sys), get_B(sys), get_C(sys), get_D(sys) + +# Funtions for number of intputs, outputs and states +ninputs(sys::StateSpace) = size(get_D(sys), 2) +noutputs(sys::StateSpace) = size(get_D(sys), 1) +nstates(sys::StateSpace) = size(get_A(sys), 1) + +##################################################################### +## Math Operators ## +##################################################################### + +## EQUALITY ## +function ==(sys1::StateSpace, sys2::StateSpace) + return all(getfield(sys1, f) == getfield(sys2, f) for f in fieldnames(StateSpace)) +end + +## Approximate ## +function isapprox(sys1::StateSpace, sys2::StateSpace) + fieldsApprox = [:A, :B, :C, :D, :Ts] + return all(sys1.f ≈ sys2.f for f in fieldnames(fieldsApprox)) +end + +## ADDITION ## +function +(s1::StateSpace{T,MT}, s2::StateSpace{T,MT}) where {T, MT} + #Ensure systems have same dimensions + if size(s1) != size(s2) + error("Systems have different shapes.") + elseif s1.Ts != s2.Ts + error("Sampling time mismatch") + end + + A = [s1.A fill(zero(T), nstates(s1), nstates(s2)); + fill(zero(T), nstates(s2), nstates(s1)) s2.A] + B = [s1.B ; s2.B] + C = [s1.C s2.C;] + D = [s1.D + s2.D;] + + return StateSpace(A, B, C, D, s1.Ts) +end + ++(sys::StateSpace, n::Number) = StateSpace(sys.A, sys.B, sys.C, sys.D .+ n, sys.Ts) ++(n::Number, sys::StateSpace) = +(sys, n) + +## SUBTRACTION ## +-(sys1::StateSpace, sys2::StateSpace) = +(sys1, -sys2) +-(sys::StateSpace, n::Number) = +(sys, -n) +-(n::Number, sys::StateSpace) = +(-sys, n) + +## NEGATION ## +-(sys::StateSpace) = StateSpace(sys.A, sys.B, -sys.C, -sys.D, sys.Ts) + +## MULTIPLICATION ## +function *(sys1::StateSpace{T,MT}, sys2::StateSpace{T,MT}) where {T, MT} + #Check dimension alignment + #Note: sys1*sys2 = y <- sys1 <- sys2 <- u + if sys1.nu != sys2.ny + error("sys1*sys2: sys1 must have same number of inputs as sys2 has outputs") + elseif sys1.Ts != sys2.Ts + error("Sampling time mismatch") + end + + A = [sys1.A sys1.B*sys2.C; + fill(zero(T), sys2.nx, sys1.nx) sys2.A] + B = [sys1.B*sys2.D ; sys2.B] + C = [sys1.C sys1.D*sys2.C;] + D = [sys1.D*sys2.D;] + + return StateSpace{T,MT}(A, B, C, D, sys2.Ts) +end + +*(sys::StateSpace, n::Number) = StateSpace(sys.A, sys.B, sys.C*n, sys.D*n, sys.Ts) +*(n::Number, sys::StateSpace) = *(sys, n) + +## DIVISION ## +/(sys1::StateSpace, sys2::StateSpace) = sys1*inv(sys2) + +function /(n::Number, sys::StateSpace) + # Ensure s.D is invertible + A, B, C, D = ssdata(sys) + Dinv = try + inv(D) + catch + error("D isn't invertible") + end + return StateSpace(A - B*Dinv*C, B*Dinv, -n*Dinv*C, n*Dinv, get_Ts(sys)) +end + +Base.inv(sys::StateSpace) = 1/sys +/(sys::StateSpace, n::Number) = StateSpace(sys.A, sys.B, sys.C/n, sys.D/n, sys.Ts) + + + +##################################################################### +## Indexing Functions ## +##################################################################### +Base.ndims(::StateSpace) = 2 # NOTE: Also for SISO systems? +Base.size(sys::StateSpace) = (noutputs(sys), ninputs(sys)) # NOTE: or just size(get_D(sys)) +Base.size(sys::StateSpace, d) = d <= 2 ? size(sys)[d] : 1 +Base.eltype(::Type{S}) where {S<:StateSpace} = S + +function Base.getindex(sys::StateSpace, inds...) + if size(inds, 1) != 2 + error("Must specify 2 indices to index statespace model") + end + rows, cols = index2range(inds...) # FIXME: ControlSystems.index2range(inds...) + return StateSpace(copy(sys.A), sys.B[:, cols], sys.C[rows, :], sys.D[rows, cols], sys.Ts) +end + +##################################################################### +## Display Functions ## +##################################################################### + +# TODO : this is a very hacky way of handling StateSpace printing. +# function _string_mat_with_headers(X::Matrix, cols::Vector{String}, +# rows::Vector{String}) +# mat = [[""] reshape(cols,1,length(cols)); +# rows X] +# p = (io, m) -> Base.showarray(io, m, false, header=false) +# return replace(sprint(p, mat), "\"", " ") +# end + +function _string_mat_with_headers(X::Matrix) + #mat = [[""] reshape(cols,1,length(cols)); + # rows X] + p = (io, m) -> Base.showarray(io, m, false, header=false) + return replace(sprint(p, X), "\"", " ") +end + +Base.print(io::IO, sys::StateSpace) = show(io, sys) + +function Base.show(io::IO, sys::StateSpace) + # Compose the name vectors + #inputs = format_names(s.inputnames, "u", "?") + #outputs = format_names(s.outputnames, "y", "?") + #println(io, "StateSpace:") + println(io, typeof(sys)) + if nstates(sys) > 0 + #states = format_names(s.statenames, "x", "?") + println(io, "A = \n", _string_mat_with_headers(sys.A)) + println(io, "B = \n", _string_mat_with_headers(sys.B)) + println(io, "C = \n", _string_mat_with_headers(sys.C)) + end + println(io, "D = \n", _string_mat_with_headers(sys.D), "\n") + # Print sample time + if sys.Ts > 0 + println(io, "Sample Time: ", sys.Ts, " (seconds)") + elseif sys.Ts == -1 + println(io, "Sample Time: unspecified") + end + # Print model type + if nstates(sys) == 0 + print(io, "Static gain") # NOTE: Not quite...still has a time type + elseif iscontinuous(sys) + print(io, "Continuous-time state-space model") + else + print(io, "Discrete-time state-space model") + end +end + + + + +""" +`minsys = minreal(s::StateSpace, tol=sqrt(eps()))` is implemented via `baltrunc` and returns a system on diagonal form. +""" +function minreal(s::StateSpace, tol=sqrt(eps())) + s = baltrunc(s, atol=tol, rtol = 0)[1] + try + return diagonalize(s) + catch + error("Minreal only implemented for diagonalizable systems.") + end +end + + + +""" +`dsys = diagonalize(s::StateSpace, digits=12)` Diagonalizes the system such that the A-matrix is diagonal. The result is rounded to `digits` decimal points. +""" +function diagonalize(s::StateSpace, digits = 12) + r = x -> round(x,digits) + S,V = eig(s.A) + try + A = V\s.A*V .|> r + B = V\s.B .|> r + C = s.C*V .|> r + D = s.D .|> r + return ss(A,B,C,D) + catch e + error("System not diagonalizable", e) + end +end diff --git a/src/types/TransferFunction.jl b/src/types/TransferFunction.jl new file mode 100644 index 000000000..7796bf9e6 --- /dev/null +++ b/src/types/TransferFunction.jl @@ -0,0 +1,202 @@ +type TransferFunction{S<:SisoTf{T} where T} <: LTISystem + matrix::Matrix{S} + Ts::Float64 + nu::Int + ny::Int + function TransferFunction{S}(matrix::Matrix{S}, Ts::Real=0) where {S} + # Validate size of input and output names + ny, nu = size(matrix) + # Validate sampling time + if Ts < 0 && Ts != -1 + error("Ts must be either a positive number, 0 + (continuous system), or -1 (unspecified)") + end + return new{S}(matrix, Ts, nu, ny) + end +end +function TransferFunction(matrix::Matrix{S}, Ts::Float64=0) where {T<:Number, S<:SisoTf{T}} + return TransferFunction{S}(matrix, Ts) +end + +noutputs(G::TransferFunction) = size(G.matrix, 1) +ninputs(G::TransferFunction) = size(G.matrix, 2) + +# function TransferFunction{SisoRational{T}}(num::Matrix{Vector{T}}, den::Matrix{Vector{T}}, Ts::Real=0) where {T <: Number} +# # Validate input and output dimensions match +# ny, nu = size(num, 1, 2) +# if (ny, nu) != size(den, 1, 2) +# error("num and den dimensions must match") +# end +# +# matrix = Matrix{SisoRational{T}}(ny, nu) # TODO: list comprehension seems suitable here +# for o=1:ny +# for i=1:nu +# matrix[o, i] = SisoRational(num[o, i], den[o, i]) +# end +# end +# return TransferFunction(matrix, Float64(Ts)) +# end + +##################################################################### +## Misc. Functions ## +##################################################################### + +## INDEXING ## +Base.ndims(::TransferFunction) = 2 +Base.size(G::TransferFunction) = size(G.matrix) +Base.size(G::TransferFunction, d) = size(G.matrix, d) +Base.eltype(::Type{S}) where {S<:TransferFunction} = S + +function Base.getindex(G::TransferFunction{S}, inds...) where {S<:SisoTf} + if size(inds, 1) != 2 + error("Must specify 2 indices to index TransferFunction model") + end + rows, cols = index2range(inds...) + mat = Matrix{S}(length(rows), length(cols)) + mat[:, :] = G.matrix[rows, cols] + return TransferFunction(mat, G.Ts) +end + +function Base.copy(G::TransferFunction) + return TransferFunction(copy(G.matrix), G.Ts) +end + + +numpoly(G::TransferFunction) = map(numpoly, G.matrix) +denpoly(G::TransferFunction) = map(denpoly, G.matrix) + +@doc """`tf = minreal(tf::TransferFunction, eps=sqrt(eps()))` + +Create a minimial representation of each transfer function in `tf` by cancelling poles and zeros """ -> +function minreal(G::TransferFunction, eps::Real=sqrt(eps())) + matrix = similar(G.matrix) + for i = eachindex(G.matrix) + matrix[i] = minreal(G.matrix[i], eps) + end + return TransferFunction(matrix, G.Ts) +end + +@doc """`isproper(tf)` + +Returns `true` if the `TransferFunction` is proper. This means that order(den) +\>= order(num))""" -> +function isproper(G::TransferFunction) + return all(isproper(f) for f in G.matrix) +end +##################################################################### +## Math Operators ## +##################################################################### + +## EQUALITY ## +function ==(G1::TransferFunction, G2::TransferFunction) + fields = [:Ts, :ny, :nu, :matrix] + for field in fields + if getfield(G1, field) != getfield(G2, field) + return false + end + end + return true +end + +## Approximate ## +function isapprox(G1::TransferFunction, G2::TransferFunction; kwargs...) + G1, G2 = promote(G1, G2) + fieldsApprox = [:Ts, :matrix] + for field in fieldsApprox + if !(isapprox(getfield(G1, field), getfield(G2, field); kwargs...)) + return false + end + end + return true +end + +function isapprox{T<:SisoTf, S<:SisoTf}(G1::Array{T}, G2::Array{S}; kwargs...) + reduce(&, [isapprox(G1[i], G2[i]; kwargs...) for i in eachindex(G1)]) +end + +## ADDITION ## +function +(G1::TransferFunction, G2::TransferFunction) + if size(G1) != size(G2) + error("Systems have different shapes.") + elseif G1.Ts != G2.Ts + error("Sampling time mismatch") + end + + matrix = G1.matrix + G2.matrix + return TransferFunction(matrix, G1.Ts) +end + ++(G::TransferFunction, n::Number) = TransferFunction(G.matrix + n, G.Ts) ++(n::Number, G::TransferFunction) = +(G, n) + +## SUBTRACTION ## +-(n::Number, G::TransferFunction) = TransferFunction(n - G.matrix, G.Ts) +-(G1::TransferFunction, G2::TransferFunction) = +(G1, -G2) +-(G::TransferFunction, n::Number) = +(G, -n) + +## NEGATION ## +-(G::TransferFunction) = TransferFunction(-G.matrix, G.Ts) + +## MULTIPLICATION ## +function *(G1::TransferFunction, G2::TransferFunction) + # Note: G1*G2 = y <- G1 <- G2 <- u + if G1.nu != G2.ny + error("G1*G2: G1 must have same number of inputs as G2 has outputs") + elseif G1.Ts != G2.Ts + error("Sampling time mismatch") + end + + matrix = G1.matrix * G2.matrix + + return TransferFunction{eltype(matrix)}(matrix, G1.Ts) +end + +*(G::TransferFunction, n::Number) = TransferFunction(n*G.matrix, G.Ts) +*(n::Number, G::TransferFunction) = *(G, n) + +## DIVISION ## +function /(n::Number, G::TransferFunction) + if issiso(G) + matrix = reshape([n/G.matrix[1,1]], 1, 1) + else + error("MIMO TransferFunction inversion isn't implemented yet") + end + return TransferFunction(matrix, G.Ts) +end +/(G::TransferFunction, n::Number) = G*(1/n) +/(G1::TransferFunction, G2::TransferFunction) = G1*(1/G2) + +##################################################################### +## Display Functions ## +##################################################################### + +Base.print(io::IO, G::TransferFunction) = show(io, G) + +function Base.show(io::IO, G::TransferFunction) + # Compose the name vectors + #println(io, "TransferFunction:") + println(io, typeof(G)) + var = iscontinuous(G) ? :s : :z + for i=1:G.nu + for o=1:G.ny + if !issiso(G) + println(io, "Input $i to output $o") + end + print_siso(io, G.matrix[o, i], var) + if !(i == G.nu && o == G.ny) + print(io, "\n") + end + end + end + if iscontinuous(G) + print(io, "\nContinuous-time transfer function model") + else + print(io, "\nSample Time: ") + if G.Ts > 0 + print(io, G.Ts, " (seconds)") + elseif G.Ts == -1 + print(io, "unspecified") + end + print(io, "\nDiscrete-time transfer function model") + end +end diff --git a/src/types/conversion.jl b/src/types/conversion.jl new file mode 100644 index 000000000..02dfd7924 --- /dev/null +++ b/src/types/conversion.jl @@ -0,0 +1,299 @@ +# Base.convert(::Type{<:SisoTf}, b::Real) = Base.convert(SisoRational, b) +# Base.convert{T<:Real}(::Type{<:SisoZpk}, b::T) = SisoZpk(T[], T[], b) +# Base.convert{T<:Real}(::Type{<:SisoRational}, b::T) = SisoRational([b], [one(T)]) +# Base.convert{T1}(::Type{SisoRational{Vector{T1}}}, t::SisoRational) = SisoRational(Poly(T1.(t.num.a)),Poly(T1.(t.den.a))) +# Base.convert(::Type{<:StateSpace}, t::Real) = ss(t) +# + +# +# function Base.convert{T<:AbstractMatrix{<:Number}}(::Type{StateSpace{T}}, s::StateSpace) +# AT = promote_type(T, arraytype(s)) +# StateSpace{AT}(AT(s.A),AT(s.B),AT(s.C),AT(s.D), s.Ts, s.statenames, s.inputnames, s.outputnames) +# end + +# TODO Fix these to use proper constructors +# NOTE: no real need to convert numbers to transfer functions, have addition methods.. +# How to convert a number to either Continuous or Discrete transfer function +# In this case it would be motivated with a "Static" Type +# Base.convert(::Type{<:TransferFunction}, b::Number) = tf([b]) +# Base.convert(::Type{<:TransferFunction{<:SisoRational}}, b::Number) = tf(b) +# Base.convert(::Type{<:TransferFunction{<:SisoZpk}}, b::Number) = zpk(b) + +Base.convert(::Type{<:TransferFunction{<:SisoZpk{T1, TR1}}}, b::AbstractMatrix{T2}) where {T1, TR1, T2<:Number} = zpk(T1.(b)) +Base.convert(::Type{<:TransferFunction{<:SisoRational{T1}}}, b::AbstractMatrix{T2}) where {T1, T2<:Number} = tf(T1.(b)) +function convert(::Type{StateSpace{T,MT}}, D::AbstractMatrix{<:Number}) where {T, MT} + (ny, nu) = size(D) + A = MT(fill(zero(T), (0,0))) + B = MT(fill(zero(T), (0,nu))) + C = MT(fill(zero(T), (ny,0))) + D = convert(MT, D) + Ts = 0.0 + return StateSpace{T,MT}(A,B,C,D,Ts) +end + +Base.convert(::Type{<:TransferFunction{<:SisoRational{T}}}, b::Number) where {T} = tf(T(b)) +Base.convert(::Type{<:TransferFunction{<:SisoZpk{T, TR}}}, b::Number) where {T, TR} = zpk(T(b)) +Base.convert(::Type{StateSpace{T, MT}}, b::Number) where {T, MT} = convert(StateSpace{T, MT}, fill(b, (1,1))) + + +#Base.convert(::Type{<:TransferFunction{<:SisoZpk}}, s::TransferFunction) = zpk(s) +#Base.convert(::Type{<:TransferFunction{<:SisoRational}}, s::TransferFunction) = tf(s) + +# +# function Base.convert{T<:Real,S<:TransferFunction}(::Type{S}, b::VecOrMat{T}) +# r = Matrix{S}(size(b,2),1) +# for j=1:size(b,2) +# r[j] = vcat(map(k->convert(S,k),b[:,j])...) +# end +# hcat(r...) +# end +# + +function convert(::Type{TransferFunction{S}}, G::TransferFunction) where S + Gnew_matrix = Matrix{S}(size(G)) + for i in eachindex(G.matrix) + Gnew_matrix[i] = convert(S, G.matrix[i]) + end + return TransferFunction{S}(Gnew_matrix, G.Ts) +end + +function convert(::Type{StateSpace{T,MT}}, sys::StateSpace) where {T, MT} + return StateSpace{T,MT}(convert(MT, sys.A), convert(MT, sys.B), convert(MT, sys.C), convert(MT, sys.D), sys.Ts) +end + +function Base.convert(::Type{StateSpace}, G::TransferFunction{<:SisoTf{T0}}) where {T0<:Number} + T = Base.promote_op(/,T0,T0) + convert(StateSpace{T,Matrix{T}}, G) +end + + +function Base.convert(::Type{StateSpace{T,MT}}, G::TransferFunction) where {T<:Number, MT<:AbstractArray{T}} + if !isproper(G) + error("System is improper, a state-space representation is impossible") + end + + # TODO : These are added due to scoped for blocks, but is a hack. This + # could be much cleaner. + #T = Base.promote_op(/, T0, T0) + + Ac = Bc = Cc = Dc = A = B = C = D = Array{T}(0, 0) + for i=1:ninputs(G) + for j=1:noutputs(G) + a, b, c, d = siso_tf_to_ss(T, G.matrix[j, i]) + if j > 1 + # vcat + Ac = blkdiag(Ac, a) + Bc = vcat(Bc, b) + Cc = blkdiag(Cc, c) + Dc = vcat(Dc, d) + else + Ac, Bc, Cc, Dc = a, b, c, d + end + end + if i > 1 + # hcat + A = blkdiag(A, Ac) + B = blkdiag(B, Bc) + C = hcat(C, Cc) + D = hcat(D, Dc) + else + A, B, C, D = Ac, Bc, Cc, Dc + end + end + # A, B, C = balance_statespace(A, B, C)[1:3] NOTE: Use balance? + return StateSpace{T,MT}(A, B, C, D, G.Ts) +end + +siso_tf_to_ss(T::Type, f::SisoTf) = siso_tf_to_ss(T, convert(SisoRational, f)) + +# Conversion to statespace on controllable canonical form +function siso_tf_to_ss(T::Type, f::SisoRational) + + num0, den0 = Base.num(f), Base.den(f) + # Normalize the numerator and denominator, + # To allow realization of transfer functions that are proper, but now strictly proper + num = num0 / den0[1] + den = den0 / den0[1] + + N = length(den) - 1 # The order of the rational function f + + # Get numerator coefficient of the same order as the denominator + bN = length(num) == N+1 ? num[1] : 0 + + if N==0 #|| num == zero(Poly{T}) + A = fill(zero(T), 0, 0) + B = fill(zero(T), 0, 1) + C = fill(zero(T), 1, 0) + else + A = diagm(fill(one(T), N-1), 1) + A[end, :] .= -reverse(den)[1:end-1] + + B = fill(zero(T), N, 1) + B[end] = one(T) + + C = fill(zero(T), 1, N) + C[1:min(N, length(num))] = reverse(num)[1:min(N, length(num))] + C[:] -= bN * reverse(den)[1:end-1] # Can index into polynomials at greater inddices than their length + end + D = fill(bN, 1, 1) + + return A, B, C, D +end + + +""" +`A, B, C, T = balance_statespace{S}(A::Matrix{S}, B::Matrix{S}, C::Matrix{S}, perm::Bool=false)` + +`sys, T = balance_statespace(sys::StateSpace, perm::Bool=false)` + +Computes a balancing transformation `T` that attempts to scale the system so +that the row and column norms of [T*A/T T*B; C/T 0] are approximately equal. +If `perm=true`, the states in `A` are allowed to be reordered. + +This is not the same as finding a balanced realization with equal and diagonal observability and reachability gramians, see `balreal` +""" +function balance_statespace(A::AbstractMatrix{P}, B::AbstractMatrix{P}, C::AbstractMatrix{P}, perm::Bool=false) where P <: BlasFloat + nx = size(A, 1) + nu = size(B, 2) + ny = size(C, 1) + + # Compute the transformation matrix + mag_A = abs.(A) + mag_B = max.(abs.(B), zero(P)) + mag_C = max.(abs.(C), zero(P)) + T = balance_transform(mag_A, mag_B, mag_C, perm) + + # Perform the transformation + A = T*A/T + B = T*B + C = C/T + + return A, B, C, T +end + +#TODO Throw a warning here at least? +# Fallback mehod for systems with exotic matrices (i.e. TrackedArrays) +balance_statespace(A::AbstractMatrix, B::AbstractMatrix, C::AbstractMatrix, perm::Bool=false) = A,B,C,I + +# First try to promote and see of we get BlasFloat, otherwise, fall beack on function above +function balance_statespace(A::Matrix{<:Number}, B::Matrix{<:Number}, C::Matrix{<:Number}, D::Matrix{<:Number}, perm::Bool=false) + T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D)) + A2, B2, C2, D2 = promote(A,B,C,D, fill(zero(T)/one(T),0,0)) # If Int, we get Float64 + balance_statespace(A2, B2, C2, perm) +end + +function balance_statespace(sys::StateSpace, perm::Bool=false) + A, B, C, T = balance_statespace(sys.A,sys.B,sys.C, perm) + return ss(A,B,C,sys.D), T +end + +""" +`T = balance_transform{R}(A::Matrix{R}, B::Matrix{R}, C::Matrix{R}, perm::Bool=false)` + +`T = balance_transform(sys::StateSpace, perm::Bool=false) = balance_transform(A,B,C,perm)` + +Computes a balancing transformation `T` that attempts to scale the system so +that the row and column norms of [T*A/T T*B; C/T 0] are approximately equal. +If `perm=true`, the states in `A` are allowed to be reordered. + +This is not the same as finding a balanced realization with equal and diagonal observability and reachability gramians, see `balreal` +See also `balance_statespace`, `balance` +""" +function balance_transform(A::AbstractArray{R}, B::AbstractArray{R}, C::AbstractArray{R}, perm::Bool=false) where {R<:BlasFloat} + nx = size(A, 1) + # Compute a scaling of the system matrix M + T = [A B; C zeros(R, size(C*B))] + size(T,1) < size(T,2) && (T = [T; zeros(R, size(T,2)-size(T,1),size(T,2))]) + size(T,1) > size(T,2) && (T = [T zeros(R, size(T,1),size(T,1)-size(T,2))]) + S = diag(balance(T, false)[1]) + Sx = S[1:nx] + Sio = S[nx+1] + # Compute permutation of x (if requested) + pvec = perm ? balance(A, true)[2] * [1:nx;] : [1:nx;] + # Compute the transformation matrix + T = zeros(R, nx, nx) + T[pvec, :] = Sio * diagm(R(1)./Sx) + return T +end + +balance_transform(sys::StateSpace, perm::Bool=false) = balance_transform(sys.A,sys.B,sys.C,perm) + + +convert(::Type{TransferFunction}, sys::StateSpace) = convert(TransferFunction{SisoRational}, sys) + +function convert(::Type{TransferFunction{SisoRational{T}}}, sys::StateSpace) where {T<:Number} + matrix = Matrix{SisoRational{T}}(size(sys)) + + A, B, C, D = ssdata(sys) + + # The following follows from the matrix inversion lemma: + # det(X + uᵀv) = det(X)(1 + vᵀX⁻¹u), or + # det((sI-A)+BC) = det(sI-A)(1 + C(si-A)⁻¹B) + # C(si-A)⁻¹B) + D = 1/det(sI-A) * (det((sI-A)+BC) - I + D*det(sI-A)) + charpolyA = charpoly(A) + for i=1:ninputs(sys), j=1:noutputs(sys) + num = charpoly(A-B[:,i:i]*C[j:j,:]) - charpolyA + D[j, i]*charpolyA + matrix[j, i] = SisoRational{T}(num, charpolyA) + end + TransferFunction{SisoRational{T}}(matrix, get_Ts(sys)) +end +function convert(::Type{TransferFunction{SisoRational}}, sys::StateSpace{T0}) where {T0<:Number} + T = typeof(one(T0)/one(T0)) + convert(TransferFunction{SisoRational{T}}, sys) +end + + +function convert(::Type{TransferFunction{SisoZpk{T,TR}}}, sys::StateSpace) where {T<:Number, TR <: Number} + matrix = Matrix{SisoZpk{T,TR}}(size(sys)) + + A, B, C, D = ssdata(sys) + + for j=1:noutputs(sys), i=1:ninputs(sys) + z, p, k = siso_ss_to_zpk(sys, i, j) + matrix[i, j] = SisoZpk{T,TR}(z, p, k) + end + TransferFunction{SisoZpk{T,TR}}(matrix, get_Ts(sys)) +end +function convert(::Type{TransferFunction{SisoZpk}}, sys::StateSpace{T0}) where {T0<:Number} + T = typeof(one(T0)/one(T0)) + convert(TransferFunction{SisoZpk{T,complex(T)}}, sys) +end + +function siso_ss_to_zpk(sys, i, j) + A, B, C = struct_ctrb_obsv(sys.A, sys.B[:, j:j], sys.C[i:i, :]) + D = sys.D[i:i, j:j] + z = tzero(A, B, C, D) + nx = size(A, 1) + nz = length(z) + k = nz == nx ? D[1] : (C*(A^(nx - nz - 1))*B)[1] + return z, eigvals(A), k +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) + + return prod(roots2poly_factors(Λ)) # Compute the polynomial factors directly? +end +function charpoly(A::AbstractMatrix{<:Real}) + Λ = eigvals(A) + return prod(roots2real_poly_factors(Λ)) +end + + +# function charpoly(A) +# λ = eigvals(A); +# T = promote_type(primitivetype(A), Float64) +# I = one(T) +# p = reduce(*,Poly([I]), Poly[Poly([I, -λᵢ]) for λᵢ in λ]); +# if maximum(imag.(p[:])./(I+abs.(real.(p[:])))) < sqrt(eps(T)) +# for i = 1:length(p) +# p[i] = real(p[i]) +# end +# else +# error("Characteristic polynomial should be real") +# end +# p +# end diff --git a/src/types/lti.jl b/src/types/lti.jl deleted file mode 100644 index c74a8605c..000000000 --- a/src/types/lti.jl +++ /dev/null @@ -1,6 +0,0 @@ -abstract type LTISystem end - -+(sys1::LTISystem, sys2::LTISystem) = +(promote(sys1, sys2)...) --(sys1::LTISystem, sys2::LTISystem) = -(promote(sys1, sys2)...) -*(sys1::LTISystem, sys2::LTISystem) = *(promote(sys1, sys2)...) -/(sys1::LTISystem, sys2::LTISystem) = /(promote(sys1, sys2)...) diff --git a/src/types/polys.jl b/src/types/polys.jl deleted file mode 100644 index 517e02d66..000000000 --- a/src/types/polys.jl +++ /dev/null @@ -1,231 +0,0 @@ -# Taken and modified from https://github.com/vtjnash/Polynomial.jl -# Polynomial.jl was deprecated in favor of Polynomials.jl, which uses reversed -# indexing. The former was copied to this file, and modified/trimmed to just -# the bare functions required for TransferFunction support. - -Base.eps{T}(::Type{T}) = zero(T) -Base.eps{F<:AbstractFloat}(x::Type{F}) = Base.eps(F) -Base.eps{T}(x::Type{Complex{T}}) = eps(T) - -immutable Poly{T<:Number} - a::Vector{T} - nzfirst::Int #for effiencicy, track the first non-zero index - function Poly{S}(a::Vector{S}) where S - la = length(a) - i = 0 - for i = 1:la - if abs(a[i]) > 2*eps(S) break end - end - new(a, i) - end -end - -Poly{T<:Number}(a::Vector{T}) = Poly{T}(a) - -Base.convert{T}(::Type{Poly{T}}, p::Poly) = Poly(convert(Vector{T}, p.a)) -Base.promote_rule{T, S}(::Type{Poly{T}}, ::Type{Poly{S}}) = Poly{promote_type(T, S)} -Base.eltype{T}(::Poly{T}) = T - -Base.length(p::Poly) = length(p.a) - p.nzfirst + 1 -Base.endof(p::Poly) = length(p) -deg(p::Poly) = length(p) - 1 - -Base.getindex(p::Poly, i) = p.a[i - 1 + p.nzfirst] -Base.getindex(p::Poly, c::Colon) = p.a[p.nzfirst:end] -Base.setindex!(p::Poly, v, i) = (p.a[i - 1 + p.nzfirst] = v) - -Base.copy(p::Poly) = Poly(copy(p.a[p.nzfirst:end])) - -Base.zero{T}(p::Poly{T}) = Poly([zero(T)]) -Base.zero{T}(::Type{Poly{T}}) = Poly([zero(T)]) -Base.one{T}(p::Poly{T}) = Poly([one(T)]) -Base.one{T}(::Type{Poly{T}}) = Poly([one(T)]) - -function Base.show(io::IO, p::Poly) - print(io,"Poly(") - print(io,p) - print(io,")") -end - -Base.print(io::IO, p::Poly) = print_poly(io, p) - -function print_poly{T}(io::IO, p::Poly{T}, var=:x) - n = length(p) - if n == 1 - print(io, p[1]) - else - for j = 1:n - pj = p[j] - magpj = abs(pj) - if magpj > 2*eps(T) - if j == 1 - real(pj) < 0 && print(io, "-") #Prepend - if first and negative - else - real(pj) < 0 ? print(io," - ") : print(io," + ") - end - #Print pj if pj is the last coefficient, or pj is not identically 1 - if j == n || abs(magpj - 1) > 2*eps(T) - print(io, magpj) - end - exp = n-j - if exp > 0 - print(io, var) - if exp > 1 - print(io, '^', exp) - end - end - end - end - end -end - -*(c::Number, p::Poly) = Poly(c * p.a[p.nzfirst:end]) -*(p::Poly, c::Number) = Poly(c * p.a[p.nzfirst:end]) -/(p::Poly, c::Number) = Poly(p.a[p.nzfirst:end] / c) -#./(p::Poly, c::Number) = /(p, c) --(p::Poly) = Poly(-p.a[p.nzfirst:end]) - --(p::Poly, c::Number) = +(p, -c) -+(c::Number, p::Poly) = +(p, c) -function +(p::Poly, c::Number) - if length(p) < 1 - return Poly([c,]) - else - p2 = copy(p) - p2.a[end] += c; - return p2; - end -end -function -(c::Number, p::Poly) - if length(p) < 1 - return Poly([c,]) - else - p2 = -p; - p2.a[end] += c; - return p2; - end -end - -function +{T,S}(p1::Poly{T}, p2::Poly{S}) - R = promote_type(T,S) - n = length(p1) - m = length(p2) - if n > m - a = Array{R}(n) - for i = 1:m - a[n-m+i] = p1[n-m+i] + p2[i] - end - for i = 1:n-m - a[i] = p1[i] - end - else - a = Array{R}(m) - for i = 1:n - a[m-n+i] = p1[i] + p2[m-n+i] - end - for i = 1:m-n - a[i] = p2[i] - end - end - Poly(a) -end - -function -{T,S}(p1::Poly{T}, p2::Poly{S}) - R = promote_type(T,S) - n = length(p1) - m = length(p2) - if n > m - a = Array{R}(n) - for i = 1:m - a[n-m+i] = p1[n-m+i] - p2[i] - end - for i = 1:n-m - a[i] = p1[i] - end - else - a = Array{R}(m) - for i = 1:n - a[m-n+i] = p1[i] - p2[m-n+i] - end - for i = 1:m-n - a[i] = -p2[i] - end - end - Poly(a) -end - -function *{T,S}(p1::Poly{T}, p2::Poly{S}) - R = promote_type(T,S) - n = length(p1) - m = length(p2) - if n == 0 || m == 0 - return Poly(R[]) - end - a = zeros(R, n+m-1) - for i = 1:length(p1) - for j = 1:length(p2) - a[i+j-1] += p1[i] * p2[j] - end - end - Poly(a) -end - -function ==(p1::Poly, p2::Poly) - if length(p1) != length(p2) - return false - else - return p1.a[p1.nzfirst:end] == p2.a[p2.nzfirst:end] - end -end - -function isapprox(p1::Poly, p2::Poly; kwargs...) - if length(p1) != length(p2) - return false - else - return isapprox(p1.a[p1.nzfirst:end], p2.a[p2.nzfirst:end]; kwargs...) - end -end - -function polyval{T}(p::Poly{T}, x::Number) - R = promote_type(T, typeof(x)) - lenp = length(p) - if lenp == 0 - return zero(R) - else - y = convert(R, p[1]) - for i = 2:lenp - y = p[i] + x.*y - end - return y - end -end - -# compute the roots of a polynomial -function roots{T}(p::Poly{T}) - R = promote_type(T, Float64) - num_zeros = 0 - if length(p) == 0 - return zeros(R, 0) - end - while abs(p[end-num_zeros]) <= 2*eps(T) - if num_zeros == length(p)-1 - return zeros(R, 0) - end - num_zeros += 1 - end - n = length(p)-num_zeros-1 - if n < 1 - return zeros(R, length(p)-1) - end - companion = zeros(R, n, n) - a0 = p[end-num_zeros] - for i = 1:n-1 - companion[1,i] = -p[end-num_zeros-i] / a0 - companion[i+1,i] = 1; - end - companion[1,end] = -p[1] / a0 - D,V = eig(companion) - r = zeros(eltype(D),length(p)-1) - r[1:n] = 1./D - return r -end diff --git a/src/types/promotion.jl b/src/types/promotion.jl new file mode 100644 index 000000000..37fec014e --- /dev/null +++ b/src/types/promotion.jl @@ -0,0 +1,123 @@ +# Base.eltype{T}(::Type{Poly{T}}) = T +# Base.eltype{T}(::Type{TransferFunction{T}}) = T +# Base.eltype{T}(::Type{SisoZpk{T}}) = T +# Base.eltype{T}(::Type{SisoRational{T}}) = T +# Base.eltype{T}(::Type{StateSpace{T}}) = T +# + + +# """ +# responsetype(T) +# Return the numerical type of the time-response for composite LTISystem type `T`, e.g., `responsetype(ss(1.)) == Float64` +# If `T` is not a type, return `responsetype(typeof(T))` +# """ +# responsetype{T<:Real}(::Type{T}) = T +# responsetype{T<:Complex}(::Type{T}) = promote_type(T.types...) +# responsetype(x) = typeof(x) # Fallback method +# responsetype{T}(::Type{T}) = T # Catch all + + +#matrix_type{T}(::Type{StateSpace{T, M{T}}}) = T + +# Abstract type pyramid ============================================================= + +# NOTE: mycket som inte verkar användas... +#Base.promote_rule(::Type{StateSpace{T1}}, ::Type{StateSpace{T2}}) where {T1,T2} = StateSpace{promote_type(T1, T2)} + +# NOTE: Is the below thing correct always? +Base.promote_rule(::Type{StateSpace{T1,MT1}}, ::Type{StateSpace{T2,MT2}}) where {T1,T2,MT1,MT2} = StateSpace{promote_type(T1, T2),promote_type(MT1, MT2)} + +function Base.promote_rule(::Type{TransferFunction{S1}}, ::Type{StateSpace{T2,MT}}) where {T1,S1<:SisoTf{T1},T2,MT} + #promote_type(Base.promote_op(/, T1, T1), T2) + StateSpace{promote_type(T1,T2), promote_type(Matrix{T1},MT)} +end +# NOTE: Perhaps should try to keep matrix structure? + +Base.promote_rule(::Type{TransferFunction{S1}}, ::Type{TransferFunction{S2}}) where {S1, S2} = TransferFunction{promote_type(S1, S2)} +#Base.promote_rule(::Type{SisoTf}, ::Type{TransferFunction}) = TransferFunction +#Base.promote_rule(::Type{SisoZpk}, ::Type{TransferFunction}) = TransferFunction +#Base.promote_rule(::Type{SisoRational}, ::Type{TransferFunction}) = TransferFunction + + +#function Base.promote_rule{T<:StateSpace,P<:TransferFunction}(::Type{T}, ::Type{P}) #where T <: StateSpace where P <: TransferFunction +# S = promote_type(primitivereal(P), primitivereal(T)) # TODO: this is not correct for P <: zpk, in which case the StateSpace system gets complex matrices +# StateSpace{Matrix{S}} +#end + + +# Promotion of Number ========================================================== + +#Base.promote_rule{T1<:Number,T2<:Number}(::Type{TransferFunction{S}}, ::T2) = promote_type(T1, T2) + +Base.promote_rule(::Type{TransferFunction{S}}, ::Type{T2}) where {S<:SisoTf, T2<:Number} = + TransferFunction{promote_type(S, T2)} + +# TODO Figure out a better way +function Base.promote_rule(::Type{StateSpace{T1, MT}}, ::Type{T2}) where {T1, MT, T2<:Number} + NewMT = Base.promote_op(*, MT, T2) + StateSpace{eltype(NewMT), NewMT} +end + +Base.promote_rule(::Type{TransferFunction{SisoZpk{T1,TR1}}}, ::Type{M2}) where {T1, TR1, T2, M2<:AbstractMatrix{T2}} = TransferFunction{SisoZpk{T1, promote_type(TR1, T2)}} + +Base.promote_rule(::Type{TransferFunction{SisoRational{T1}}}, ::Type{M2}) where {T1, T2, M2<:AbstractMatrix{T2}} = TransferFunction{SisoRational{promote_type(T1, T2)}} + +#Base.promote_rule{S<:TransferFunction{<:SisoTf}}(::Type{S}, ::Type{<:Real}) = S + +# We want this, but not possible, so hardcode for SisoTypes +#Base.promote_rule(::Type{S}, ::Type{T2}) where {T1, S<:SisoTf{T1}, T2<:Number} = S{promote_type(T1, T2)} + +function Base.promote_rule(::Type{SisoZpk{T1,C2}}, ::Type{T2}) where {T1, C2, T2<:Number} + GainType = promote_type(T1, T2) + return SisoZpk{GainType, complex(GainType)} +end +Base.promote_rule(::Type{SisoRational{T1}}, ::Type{T2}) where {T1, T2<:Number} = SisoRational{promote_type(T1,T2)} +#Base.promote_rule(::Type{StateSpace{T1, MT}}, ::Type{T2}) where {T1, T2<:Number} = SisoRational{promote_type(T1,T2)} + +#Base.promote_rule{T<:SisoZpk}(::Type{T}, ::Type{<:Real}) = T +#Base.promote_rule{T}(::Type{SisoRational{T}}, ::Type{<:Real}) = SisoRational{T} +#Base.promote_rule{T}(::Type{StateSpace{T}}, ::Type{<:Real}) = StateSpace{T} +#Base.promote_rule{S<:TransferFunction{<:SisoTf},T<:Real}(::Type{S}, ::Union{Type{Array{T,2}},Type{Array{T,1}}}) = S + + +# Less abstract promotions +#Base.promote_rule(::Type{TransferFunction{SisoRational}}, ::Type{TransferFunction{SisoZpk}}) = TransferFunction{SisoZpk} # NOTE: Is this what we want? + + + + +# function Base.promote_rule{T1c<:SisoZpk,T2c<:SisoRational}(::Type{T1c}, ::Type{T2c}) +# upper_type = SisoZpk +# inner_type = promote_type(arraytype(T1c), arraytype(T2c)) +# upper_type{inner_type} +# end + + + +# @show siso_type = promote_type(eltype.(systems)...) +# if isleaftype(siso_type) +# siso_type = siso_type.name.wrapper +# end +# # array_type = promote_type(arraytype.(systems)...) +# # mat = hcat([(e->convert(siso_type{array_type}, e)).(s.matrix) for s in systems]...) +# @show promote_type([s.matrix for s in systems]...) + + + + +# + +# +# +# # Promote_op types # NOTE: Needed? +# Base.promote_op{T<:SisoTf}(::Any, ::Type{T}, ::Type{T}) = T +# +# #Just default SisoTf to SisoRational +# SisoTf(args...) = SisoRational(args...) +# +# Base.zero(::Type{<:SisoTf}) = zero(SisoRational) +# Base.zero(::SisoTf) = zero(SisoRational) +# Base.zero(::Type{<:SisoZpk}) = SisoZpk(Float64[],Float64[],0.0) +# Base.zero(::SisoZpk) = Base.zero(SisoZpk) +# Base.zero{T}(::Type{SisoRational{T}}) = SisoRational(zero(Poly{T}), one(Poly{T})) # FIXME: Is the in analogy with how zero for SisoRational is createdzer? +# Base.zero{T}(::SisoRational{T}) = Base.zero(SisoRational{T}) diff --git a/src/types/sisogeneralized.jl b/src/types/sisogeneralized.jl deleted file mode 100644 index 9600a3369..000000000 --- a/src/types/sisogeneralized.jl +++ /dev/null @@ -1,119 +0,0 @@ -ExprLike = Union{Expr,Number,Symbol} - -## User should just use TransferFunction -immutable SisoGeneralized <: SisoTf - expr::ExprLike - function SisoGeneralized(expr::ExprLike) - if isa(expr, Expr) && length(expr.args) == 3 && expr.args[1] == :(+) && expr.args[2] == 0 - #Get rid of the zero - expr = expr.args[3] - end - new(expr) - end -end - -SisoGeneralized(str::AbstractString) = SisoGeneralized(parse(str)) - -function minreal(sys::SisoGeneralized, eps::Real=sqrt(eps())) - warn("minreal is not implemented for generalized transferfunctions, returning same system") - sys -end - -function print_siso(io::IO, t::SisoGeneralized, var=:s) - println(io, t.expr) -end - -Base.promote_rule{T<:Real}(::Type{SisoGeneralized}, ::Type{T}) = SisoGeneralized -Base.convert(::Type{SisoGeneralized}, b::Real) = SisoGeneralized(b) - -Base.zero(::Type{SisoGeneralized}) = SisoGeneralized(0) -Base.zero(::SisoGeneralized) = Base.zero(SisoGeneralized) - -Base.length(t::SisoGeneralized) = error("length is not implemented for generalized transferfunctions") -Base.num(t::SisoGeneralized) = error("num is not implemented for generalized transferfunctions") -Base.den(t::SisoGeneralized) = error("den is not implemented for generalized transferfunctions") -pole(t::SisoGeneralized) = error("pole is not implemented for generalized transferfunctions") -tzero(t::SisoGeneralized) = error("tzero is not implemented for generalized transferfunctions") - -#This makes sure that the function can compile once -#Didn't work in julia 0.6 so we now only have inefficient version -function _preprocess_for_freqresp(sys::SisoGeneralized) - sys -end - -evalfr(f::Function, freq) = f(freq) -#evalfr(sys::SisoGeneralized, freq) = _preprocess_for_freqresp(sys)(freq) - -function evalfr(sys::SisoGeneralized, freq) - _f = eval(:(s -> $(sys.expr))) - eval(Expr(:call, - function() - return _f(freq) - end)) -end - -function lsimabstract(sys::SisoGeneralized, uin, dt, Tend) - #TODO make sure U and P are of equal length, fix input arguments, Tend results in half time, make sure u interp is using Tend - N = round(Int, Tend/dt) + 1 - #T=N*dt - T = Tend - dw = pi/T - omega = linspace(-pi/dt, pi/dt, 2N+1) - u = [uin; zeros(N)] - U = fft(u) - #Pf = _preprocess_for_freqresp(sys) - #P = Complex{Float64}[evalfr(Pf, omega[i]*im) for i in 1:2N] - P = evalfr(sys, omega[1:2N]*im) - y = real(ifft(fftshift(P).*U)) - t = dt*(0:N-1) - y[1:N] -end - -function SisoRational(expr::Expr) - if length(expr.args) == 1 - error("Unexpected operator $(expr.args[1]) in expression") - end - if expr.args[1] == :+ - return reduce(+, map(t -> SisoRational(t), expr.args[2:end])) - elseif expr.args[1] == :* - return reduce(*, map(t -> SisoRational(t), expr.args[2:end])) - elseif length(expr.args) == 2 && expr.args[1] == :- - return - SisoRational(expr.args[2]) - elseif length(expr.args) == 3 - if expr.args[1] == :- - return SisoRational(expr.args[2]) - SisoRational(expr.args[3]) - elseif expr.args[1] == :/ - return SisoRational(expr.args[2]) / SisoRational(expr.args[3]) - elseif expr.args[1] == :^ - return SisoRational(expr.args[2]) ^ expr.args[3] - else - error("Unexpected operator $(expr.args[1]) in expression") - end - else - error("Could not parse \"$(expr.args[1])\" in expression") - end -end -SisoRational(expr::Symbol) = (expr == :s ? SisoRational([1, 0],[1]) : error("invalid symbol \"$exp\" only \"s\" is allowed")) -SisoRational(expr::Number) = isa(expr,Real) ? SisoRational([expr],[1]) : error("Only real numers are allowed in transferfunction") - -==(t1::SisoGeneralized, t2::SisoGeneralized) = (t1.expr == t2.expr) - -isapprox(t1::SisoGeneralized, t2::SisoGeneralized; kwargs...) = error("isapprox not implemented for generalized tf") - -+(t1::SisoGeneralized, t2::SisoGeneralized) = SisoGeneralized(:($(t1.expr) + $(t2.expr))) -+(t::SisoGeneralized, n::Real) = SisoGeneralized(:($(t.expr) + $n)) -+(n::Real, t::SisoGeneralized) = SisoGeneralized(:($n + $(t.expr))) - --(t1::SisoGeneralized, t2::SisoGeneralized) = SisoGeneralized(:($(t1.expr) - $(t2.expr))) --(n::Real, t::SisoGeneralized) = SisoGeneralized(:($n - $(t.expr))) --(t::SisoGeneralized, n::Real) = SisoGeneralized(:($(t.expr) - $n)) - --(t::SisoGeneralized) = SisoGeneralized(:(- $(t.expr))) - -*(t1::SisoGeneralized, t2::SisoGeneralized) = SisoGeneralized(:($(t1.expr) * $(t2.expr))) -*(t::SisoGeneralized, n::Real) = SisoGeneralized(:($(t.expr) * $n)) -*(n::Real, t::SisoGeneralized) = SisoGeneralized(:($n * $(t.expr))) - -/(n::Real, t::SisoGeneralized) = SisoGeneralized(:($n / $(t.expr))) -/(t::SisoGeneralized, n::Real) = SisoGeneralized(:($(t.expr) / $n)) -/(t1::SisoGeneralized, t2::SisoGeneralized) = SisoGeneralized(:($(t1.expr) / $(t2.expr))) diff --git a/src/types/sisotf.jl b/src/types/sisotf.jl deleted file mode 100644 index 802d40091..000000000 --- a/src/types/sisotf.jl +++ /dev/null @@ -1,115 +0,0 @@ -## User should just use TransferFunction -immutable SisoRational <: SisoTf - num::Poly{Float64} - den::Poly{Float64} - function SisoRational(num::Poly{Float64}, den::Poly{Float64}) - if all(den == zero(den)) - error("Zero denominator") - elseif all(num == zero(num)) - # The numerator is zero, make the denominator 1 - den = one(den) - end - new(num, den) - end -end -SisoRational(num::Vector, den::Vector) = SisoRational(Poly(map(Float64,num)), Poly(map(Float64,den))) - -function minreal(sys::SisoRational, eps::Real=sqrt(eps())) - return SisoRational(minreal(SisoZpk(sys), eps)) -end - -function print_siso(io::IO, t::SisoRational, var=:s) - # Convert the numerator and denominator to strings - numstr = sprint(print_poly, t.num, var) - denstr = sprint(print_poly, t.den, var) - - # Figure out the length of the separating line - len_num = length(numstr) - len_den = length(denstr) - dashcount = max(len_num, len_den) - - # Center the numerator or denominator - if len_num < dashcount - numstr = "$(repeat(" ", div(dashcount - len_num, 2)))$numstr" - else - denstr = "$(repeat(" ", div(dashcount - len_den, 2)))$denstr" - end - println(io, numstr) - println(io, repeat("-", dashcount)) - println(io, denstr) -end - -function print_compact(io::Base.IO, t::SisoRational, var=:s) - numstr = sprint(print_poly, t.num, var) - denstr = sprint(print_poly, t.den, var) - println(io, "($numstr)/($denstr)") -end - -Base.promote_rule{T<:Real}(::Type{SisoRational}, ::Type{T}) = SisoRational -Base.convert(::Type{SisoRational}, b::Real) = SisoRational([b], [1]) - -Base.zero(::Type{SisoRational}) = SisoRational(zero(Poly{Float64}), one(Poly{Float64})) -Base.zero(::SisoRational) = Base.zero(SisoRational) - -Base.length(t::SisoRational) = max(length(t.num), length(t.den)) - -function Base.num(t::SisoRational) - lt = length(t) - n = zeros(lt) - n[(lt - length(t.num) + 1):end] = t.num[:] - return n -end - -function Base.den(t::SisoRational) - lt = length(t) - d = zeros(lt) - d[(lt - length(t.den) + 1):end] = t.den[:] - return d -end - -denpoly(G::SisoRational) = Poly(den(G)) -numpoly(G::SisoRational) = Poly(num(G)) - -function evalfr(sys::SisoRational, s::Number) - S = promote_type(typeof(s), Float64) - den = polyval(sys.den, s) - if den == zero(S) - convert(S, Inf) - else - polyval(sys.num, s)/den - end -end - -tzero(t::SisoRational) = roots(t.num) -pole(t::SisoRational) = roots(t.den) - -==(t1::SisoRational, t2::SisoRational) = (t1.num == t2.num && t1.den == t2.den) -# We might want to consider alowing scaled num and den as equal -function isapprox(t1::SisoRational, t2::SisoRational; rtol::Real=sqrt(eps()), atol::Real=0) - isapprox(t1.num,t2.num, rtol=rtol, atol=atol) && isapprox(t1.den, t2.den, rtol=rtol, atol=atol) -end - -+(t1::SisoRational, t2::SisoRational) = SisoRational(t1.num*t2.den + t2.num*t1.den, t1.den*t2.den) -+(t::SisoRational, n::Real) = SisoRational(t.num + n*t.den, t.den) -+(n::Real, t::SisoRational) = t + n -#.+(t::SisoRational, n::Real) = t + n -#.+(n::Real, t::SisoRational) = t + n - --(t1::SisoRational, t2::SisoRational) = SisoRational(t1.num*t2.den - t2.num*t1.den, t1.den*t2.den) --(n::Real, t::SisoRational) = SisoRational(n*t.den - t.num, t.den) --(t::SisoRational, n::Real) = +(t, -n) -#.-(t::SisoRational, n::Real) = -(t, n) -#.-(n::Real, t::SisoRational) = -(n, t) - --(t::SisoRational) = SisoRational(-t.num, t.den) - -*(t1::SisoRational, t2::SisoRational) = SisoRational(t1.num*t2.num, t1.den*t2.den) -*(t::SisoRational, n::Real) = SisoRational(t.num*n, t.den) -*(n::Real, t::SisoRational) = *(t, n) -#.*(t1::SisoRational, t2::SisoRational) = *(t1, t2) -#.*(t::SisoRational, n::Real) = *(t, n) -#.*(n::Real, t::SisoRational) = *(t, n) - -/(n::Real, t::SisoRational) = SisoRational(n*t.den, t.num) -/(t::SisoRational, n::Real) = t*(1/n) -/(t1::SisoRational, t2::SisoRational) = t1*(1/t2) diff --git a/src/types/sisozpk.jl b/src/types/sisozpk.jl deleted file mode 100644 index 957832d29..000000000 --- a/src/types/sisozpk.jl +++ /dev/null @@ -1,193 +0,0 @@ -RealOrComplex = Union{Real,Complex} - -## User should just use TransferFunction -immutable SisoZpk <: SisoTf - z::Vector{Complex{Float64}} - p::Vector{Complex{Float64}} - k::Float64 - function SisoZpk(z::Vector{Complex{Float64}}, p::Vector{Complex{Float64}}, k::Float64) - if k == zero(k) - p = [] - z = [] - end - new(z, p, k) - end -end - -SisoZpk{T<:RealOrComplex,S<:RealOrComplex}(z::AbstractArray{T}, p::AbstractArray{S}, k::Real) = SisoZpk(Complex128[z...], Complex128[p...], Float64(k)) - -# Taking care of empty vectors being of type Any -function SisoZpk(z::AbstractArray, p::AbstractArray, k::Real) - if length(z) > 0 - if !(eltype(z) <: RealOrComplex) - error("Zeros must be real or complex") - end - else - z = Array{Float64}(0) - end - if length(p) > 0 - if !(eltype(p) <: RealOrComplex) - error("poles must be real or complex") - end - else - p = Array{Float64}(0) - end - SisoZpk(z, p, k) -end - -function minreal(sys::SisoZpk, eps::Real) - newZ = copy(sys.z) - newP = Vector{Complex{Float64}}() - doubles = Vector{Int64}() - newZ = copy(sys.z) - for p in sys.p - if !isempty(newZ) - val, zi = findmin(abs.(p-newZ)) - else - val = Inf #Keep looping over p, adding poles - end - if val < eps - deleteat!(newZ, zi) - continue - else - push!(newP, p) - end - end - SisoZpk(newZ, newP, sys.k) -end - -function Base.num(t::SisoZpk) - return copy(t.z) -end - -function Base.den(t::SisoZpk) - return copy(t.p) -end - -tzero(sys::SisoZpk) = num(sys) -pole(sys::SisoZpk) = den(sys) - -function zp2polys(vec) - polys = Array{Poly{Float64},1}(0) - polesiLeft = Set(1:length(vec)) - while length(polesiLeft) > 0 - p = vec[pop!(polesiLeft)] - if abs(imag(p)) < sqrt(eps()) - push!(polys,Poly(float([1, -real(p)]))) - else - polesiLeftVec = [i for i in polesiLeft] - polesTest = Complex128[vec[polesiLeftVec]...] - val, i = findmin(abs.(polesTest-conj(p))) - val > 2*sqrt(eps()) && error("Could not find conjugate to pole") - push!(polys,Poly(float([1, -2*real(p), real(p)^2+imag(p)^2]))) - pop!(polesiLeft,polesiLeftVec[i]) - end - end - polys -end - -function numpoly(G::SisoZpk) - zpolys = zp2polys(G.z) -end - -function denpoly(G::SisoZpk) - ppolys = zp2polys(G.p) -end - -function evalfr(sys::SisoZpk, s::Number) - S = promote_type(typeof(s), Float64) - den = reduce(*, (poly) -> polyval(poly, s), zp2polys(sys.p)) - if den == zero(S) - return convert(S, Inf) - else - num = reduce(*, (poly) -> polyval(poly, s), zp2polys(sys.z)) - return sys.k*num/den - end -end - -function print_siso(io::IO, t::SisoZpk, var=:s) - zpolys = zp2polys(t.z) - ppolys = zp2polys(t.p) - # Convert the numerator and denominator to strings - if length(zpolys) < 2 - numstr = ( length(zpolys) == 0 ) ? "1.0" : sprint(print_poly, zpolys[1], var) - else - numstr = reduce(*,"",["("*sprint(print_poly, z, var)*")" for z in zpolys]) - end - if length(ppolys) < 2 - denstr = ( length(ppolys) == 0 ) ? "1.0" : sprint(print_poly, ppolys[1], var) - else - denstr = reduce(*,"",["("*sprint(print_poly, p, var)*")" for p in ppolys]) - end - # Figure out the length of the separating line - len_num = length(numstr) - len_den = length(denstr) - dashcount = max(len_num, len_den) - - # Center the numerator or denominator - if len_num < dashcount - numstr = "$(repeat(" ", div(dashcount - len_num, 2)))$numstr" - else - denstr = "$(repeat(" ", div(dashcount - len_den, 2)))$denstr" - end - - gainstr = string(t.k) - #Add spaces to account for gain string - numstr = " "^(length(gainstr))*numstr - denstr = " "^(length(gainstr))*denstr - println(io, numstr) - println(io, gainstr*repeat("-", dashcount)) - println(io, denstr) -end - -Base.promote_rule{T<:Real}(::Type{SisoZpk}, ::Type{T}) = SisoZpk -Base.convert(::Type{SisoZpk}, b::Real) = SisoZpk([], [], b) - -Base.zero(::Type{SisoZpk}) = SisoZpk([],[],0.0) -Base.zero(::SisoZpk) = Base.zero(SisoZpk) - -Base.length(t::SisoZpk) = max(length(t.z), length(t.p)) - - -==(t1::SisoZpk, t2::SisoZpk) = (t1-t2).k == 0.0 -function isapprox(t1::SisoZpk, t2::SisoZpk; rtol::Real=sqrt(eps()), atol::Real=sqrt(eps())) - tdiff = t1 - t2 - isapprox(tdiff.k, 0, atol=atol, rtol=rtol) -end - -function +(t1::SisoZpk, t2::SisoZpk) - numPoly = t1.k*prod(zp2polys(t1.z))*prod(zp2polys(t2.p))+t2.k*prod(zp2polys(t2.z))*prod(zp2polys(t1.p)) - z = roots(numPoly) - if length(numPoly) > 0 - k = numPoly[1] - p = [t1.p;t2.p] - else - k = 0 - p = [] - end - SisoZpk(z,p,k) -end - -+(t::SisoZpk, n::Real) = t + SisoZpk([],[],n) -+(n::Real, t::SisoZpk) = SisoZpk([],[],n) + t -#.+(t::SisoZpk, n::Real) = t + n -#.+(n::Real, t::SisoZpk) = t + n - --(t1::SisoZpk, t2::SisoZpk) = +(t1,-t2) --(n::Real, t::SisoZpk) = SisoZpk([],[],n) - t --(t::SisoZpk, n::Real) = +(t, -n) -#.-(t::SisoZpk, n::Real) = -(t, n) -#.-(n::Real, t::SisoZpk) = -(n, t) - --(t::SisoZpk) = SisoZpk(t.z, t.p, -t.k) - -*(t1::SisoZpk, t2::SisoZpk) = SisoZpk([t1.z;t2.z], [t1.p;t2.p], t1.k*t2.k) -*(t::SisoZpk, n::Real) = SisoZpk(t.z, t.p, t.k*n) -*(n::Real, t::SisoZpk) = *(t, n) -#.*(t1::SisoZpk, t2::SisoZpk) = *(t1, t2) -#.*(t::SisoZpk, n::Real) = *(t, n) -#.*(n::Real, t::SisoZpk) = *(t, n) - -/(n::Real, t::SisoZpk) = SisoZpk(t.p, t.z, n/t.k) -/(t::SisoZpk, n::Real) = SisoZpk(t.z, t.p, t.k/n) -/(t1::SisoZpk, t2::SisoZpk) = t1*(1/t2) diff --git a/src/types/ss.jl b/src/types/ss.jl new file mode 100644 index 000000000..736123d28 --- /dev/null +++ b/src/types/ss.jl @@ -0,0 +1,48 @@ +# Convenience constructor for creating StateSpace objects + +@doc """`ss(A,B,C,D[, Ts, statenames=..., inputnames=..., outputnames=...]) -> sys` + +Create a state-space model. +This is a continuous-time model if Ts is omitted or set to 0. +Otherwise, this is a discrete-time model with sampling period Ts. +Set Ts=-1 for a discrete-time model with unspecified sampling period. + +State, input and output names: each can be either a vector of strings (one string per dimension), +or a single string (e.g., "x"). In the latter case, an index is automatically appended to identify +the coordinates for each dimension (e.g. "x1", "x2", ...). + +`sys = ss(D[, Ts, ...])` specifies a static gain matrix D.""" -> +function ss(A::Array, B::Array, C::Array, D::Array, Ts::Real=0) + # Check the kwargs for metadata + nu = size(B, 2) + ny, nx = size(C, 1, 2) + return StateSpace(A, B, C, D, Ts) +end + +# Function for accepting scalars +function ss(A::Union{Number,Array}, B::Union{Number,Array}, C::Union{Number,Array}, D::Union{Number,Array}, Ts::Real=0) + T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D)) + A = to_matrix(T, A) + B = to_matrix(T, B) + C = to_matrix(T, C) + if D == 0 + D = fill(zero(T), size(C,1), size(B,2)) + else + D = to_matrix(T, D) + end + ss(A, B, C, D, Ts) +end + +# Function for creation of static gain +function ss(D::Array{T}, Ts::Real=0) where {T<:Number} + ny, nu = size(D, 1, 2) + A = fill(zero(T), 0, 0) + B = fill(zero(T), 0, nu) + C = fill(zero(T), ny, 0) + + return ss(A, B, C, D, Ts) +end +ss(d::Number, Ts::Real=0; kwargs...) = ss([d], Ts) + +# ss(sys) converts to StateSpace +ss(sys::LTISystem) = convert(StateSpace, sys) diff --git a/src/types/statespace.jl b/src/types/statespace.jl deleted file mode 100644 index 714b4cf99..000000000 --- a/src/types/statespace.jl +++ /dev/null @@ -1,346 +0,0 @@ -##################################################################### -## Data Type Declarations ## -##################################################################### - -type StateSpace <: LTISystem - A::Matrix{Float64} - B::Matrix{Float64} - C::Matrix{Float64} - D::Matrix{Float64} - Ts::Float64 - nx::Int - nu::Int - ny::Int - statenames::Vector{String} - inputnames::Vector{String} - outputnames::Vector{String} - - function StateSpace(A::Matrix{Float64}, B::Matrix{Float64}, - C::Matrix{Float64}, D::Matrix{Float64}, Ts::Float64, - statenames::Vector{String}, inputnames::Vector{String}, - outputnames::Vector{String}) - nx = size(A, 1) - nu = size(B, 2) - ny = size(C, 1) - - if size(A, 2) != nx && nx != 0 - error("A must be square") - elseif size(B, 1) != nx - error("B must have the same row size as A") - elseif size(C, 2) != nx - error("C must have the same column size as A") - elseif nu != size(D, 2) - error("D must have the same column size as B") - elseif ny != size(D, 1) - error("D must have the same row size as C") - end - # Validate names of state, input, and output - if size(statenames, 1) != nx - error("Must have same number of statenames as states") - elseif size(inputnames, 1) != nu - error("Must have same number of inputnames as inputs") - elseif size(outputnames, 1) != ny - error("Must have same number of outputnames as outputs") - end - # Validate sampling time - if Ts < 0 && Ts != -1 - error("Ts must be either a positive number, 0 - (continuous system), or -1 (unspecified)") - end - new(A, B, C, D, Ts, nx, nu, ny, statenames, inputnames, outputnames) - end -end - -function StateSpace(A::Array, B::Array, C::Array, D::Array, Ts::Real, - statenames::Vector{String}, inputnames::Vector{String}, - outputnames::Vector{String}) - return StateSpace(float64mat(A), float64mat(B), float64mat(C), - float64mat(D), map(Float64, Ts), statenames, inputnames, outputnames) -end - -##################################################################### -## Constructor Functions ## -##################################################################### - -@doc """`ss(A,B,C,D[, Ts, statenames=..., inputnames=..., outputnames=...]) -> sys` - -Create a state-space model. -This is a continuous-time model if Ts is omitted or set to 0. -Otherwise, this is a discrete-time model with sampling period Ts. -Set Ts=-1 for a discrete-time model with unspecified sampling period. - -State, input and output names: each can be either a vector of strings (one string per dimension), -or a single string (e.g., "x"). In the latter case, an index is automatically appended to identify -the coordinates for each dimension (e.g. "x1", "x2", ...). - -`sys = ss(D[, Ts, ...])` specifies a static gain matrix D.""" -> -function ss(A::Array, B::Array, C::Array, D::Array, Ts::Real=0; kwargs...) - # Check the kwargs for metadata - nu = size(B, 2) - ny, nx = size(C, 1, 2) - kvs = Dict(kwargs) - statenames = validate_names(kvs, :statenames, nx) - inputnames = validate_names(kvs, :inputnames, nu) - outputnames = validate_names(kvs, :outputnames, ny) - return StateSpace(A, B, C, D, Ts, statenames, inputnames, outputnames) -end - -# Function for accepting scalars -function ss(A::Union{Real,Array}, B::Union{Real,Array}, C::Union{Real,Array}, D::Union{Real,Array}, args...; kwargs...) - A = float64mat(A) - B = float64mat(B) - C = float64mat(C) - if D == 0 - D = zeros(size(C,1),size(B,2)) - else - D = float64mat(D) - end - ss(A, B, C, D, args..., kwargs...) -end - -# Function for creation of static gain -function ss(D::Array, Ts::Real=0; kwargs...) - ny, nu = size(D, 1, 2) - A = zeros(0, 0) - B = zeros(0, nu) - C = zeros(ny, 0) - - return ss(A, B, C, D, Ts, kwargs...) -end -ss(d::Real, Ts::Real=0; kwargs...) = ss([d], Ts, kwargs...) - -# ss(sys) converts to StateSpace -ss(sys::LTISystem) = convert(StateSpace, sys) - -# Create a random statespace system -function rss(nx::Int, nu::Int=1, ny::Int=1, feedthrough::Bool=true) - Q = randn(nx, nx) - A = Q*diagm(-100*abs.(randn(nx)))*Q' - B = randn(nx, nu) - C = randn(ny, nx) - if feedthrough - D = randn(ny, nu) - else - D = zeros(ny, nu) - end - return ss(A, B, C, D) -end - -##################################################################### -## Math Operators ## -##################################################################### - -## EQUALITY ## -function ==(s1::StateSpace, s2::StateSpace) - fields = [:Ts, :nx, :ny, :nu, :A, :B, :C, :D, :inputnames, :outputnames, - :statenames] - for field in fields - if getfield(s1, field) != getfield(s2, field) - return false - end - end - return true -end - -## Approximate ## -function isapprox(s1::StateSpace, s2::StateSpace) - fieldsApprox = [:Ts, :nx, :ny, :nu, :A, :B, :C, :D] - fieldsEqual = [:inputnames, :outputnames, :statenames] - for field in fieldsApprox - if !(getfield(s1, field) ≈ getfield(s2, field)) - return false - end - end - for field in fieldsEqual - if getfield(s1, field) != getfield(s2, field) - return false - end - end - return true -end - -## ADDITION ## -function +(s1::StateSpace, s2::StateSpace) - #Ensure systems have same dimensions - if size(s1) != size(s2) - error("Systems have different shapes.") - elseif s1.Ts != s2.Ts - error("Sampling time mismatch") - end - - A = [s1.A zeros(s1.nx, s2.nx); - zeros(s2.nx, s1.nx) s2.A] - B = [s1.B ; s2.B] - C = [s1.C s2.C;] - D = [s1.D + s2.D;] - - # Naming strategy: If only one sys is named, use that. If the names are the - # same, use them. If the names conflict, then they are ignored, and the - # default "" is used. - statenames = [s1.statenames; s2.statenames] - if all(s1.inputnames .== "") - inputnames = s2.inputnames - elseif all(s2.inputnames .== "") || (s1.inputnames == s2.inputnames) - inputnames = s1.inputnames - else - inputnames = fill(String(""),s1.ny) - end - if all(s1.outputnames .== "") - outputnames = s2.outputnames - elseif all(s2.outputnames .== "") || (s1.outputnames == s2.outputnames) - outputnames = s1.outputnames - else - outputnames = fill(String(""),s1.nu) - end - return StateSpace(A, B, C, D, s1.Ts, statenames, inputnames, outputnames) -end - -+(s::StateSpace, n::Real) = StateSpace(s.A, s.B, s.C, s.D .+ n, s.Ts, - s.statenames, s.inputnames, s.outputnames) -+(n::Real, s::StateSpace) = +(s, n) - -## SUBTRACTION ## --(s1::StateSpace, s2::StateSpace) = +(s1, -s2) --(s::StateSpace, n::Real) = +(s, -n) --(n::Real, s::StateSpace) = +(-s, n) - -## NEGATION ## --(s::StateSpace) = StateSpace(s.A, s.B, -s.C, -s.D, s.Ts, s.statenames, - s.inputnames, s.outputnames) - -## MULTIPLICATION ## -function *(s1::StateSpace, s2::StateSpace) - #Check dimension alignment - #Note: s1*s2 = y <- s1 <- s2 <- u - if s1.nu != s2.ny - error("s1*s2: s1 must have same number of inputs as s2 has outputs") - elseif s1.Ts != s2.Ts - error("Sampling time mismatch") - end - - A = [s1.A s1.B*s2.C; - zeros(s2.nx, s1.nx) s2.A] - B = [s1.B*s2.D ; s2.B] - C = [s1.C s1.D*s2.C;] - D = [s1.D*s2.D;] - - statenames = [s1.statenames; s2.statenames] - return StateSpace(A, B, C, D, s1.Ts, statenames, s2.inputnames, - s1.outputnames) -end - -*(s::StateSpace, n::Real) = StateSpace(s.A, s.B, s.C*n, s.D*n, s.Ts, - s.statenames, s.inputnames, s.outputnames) -*(n::Real, s::StateSpace) = *(s, n) - -## DIVISION ## -/(s1::StateSpace, s2::StateSpace) = s1*inv(s2) - -function /(n::Real, s::StateSpace) - # Ensure s.D is invertible - Dinv = try - inv(s.D) - catch - error("D isn't invertible") - end - return StateSpace(s.A - s.B*Dinv*s.C, s.B*Dinv, -n*Dinv*s.C, n*Dinv, s.Ts, - s.statenames, s.outputnames, s.inputnames) -end - -Base.inv(s::StateSpace) = 1/s -/(s::StateSpace, n::Real) = StateSpace(s.A, s.B, s.C/n, s.D/n, s.Ts, - s.statenames, s.inputnames, s.outputnames) - -##################################################################### -## Indexing Functions ## -##################################################################### -Base.ndims(::StateSpace) = 2 -Base.size(s::StateSpace) = (s.ny, s.nu) -Base.size(s::StateSpace, d) = d <= 2 ? size(s)[d] : 1 - -function Base.getindex(s::StateSpace, inds...) - if size(inds, 1) != 2 - error("Must specify 2 indices to index statespace model") - end - rows, cols = ControlSystems.index2range(inds...) - return StateSpace([s.A;], [s.B[:, cols];], [s.C[rows, :];], [s.D[rows, cols];], - s.Ts, [s.statenames;], [s.inputnames[cols];], [s.outputnames[rows];]) -end - -##################################################################### -## Display Functions ## -##################################################################### - -# TODO : this is a very hacky way of handling StateSpace printing. -function _string_mat_with_headers(X::Matrix, cols::Vector{String}, - rows::Vector{String}) - mat = [[""] reshape(cols,1,length(cols)); - rows X] - p = (io, m) -> Base.showarray(io, m, false, header=false) - return replace(sprint(p, mat), "\"", " ") -end - -Base.print(io::IO, s::StateSpace) = show(io, s) - -function Base.show(io::IO, s::StateSpace) - # Compose the name vectors - inputs = format_names(s.inputnames, "u", "?") - outputs = format_names(s.outputnames, "y", "?") - println(io, "StateSpace:") - if s.nx > 0 - states = format_names(s.statenames, "x", "?") - println(io, "A = \n", _string_mat_with_headers(s.A, states, states)) - println(io, "B = \n", _string_mat_with_headers(s.B, inputs, states)) - println(io, "C = \n", _string_mat_with_headers(s.C, states, outputs)) - end - println(io, "D = \n", _string_mat_with_headers(s.D, inputs, outputs), "\n") - # Print sample time - if s.Ts > 0 - println(io, "Sample Time: ", s.Ts, " (seconds)") - elseif s.Ts == -1 - println(io, "Sample Time: unspecified") - end - # Print model type - if s.nx == 0 - print(io, "Static gain") - elseif iscontinuous(s) - print(io, "Continuous-time state-space model") - else - print(io, "Discrete-time state-space model") - end -end - - - -##################################################################### -## Other Functions ## -##################################################################### - -""" -`minsys = minreal(s::StateSpace, tol=sqrt(eps()))` is implemented via `baltrunc` and returns a system on diagonal form. -""" -function minreal(s::StateSpace, tol=sqrt(eps())) - s = baltrunc(s, atol=tol, rtol = 0)[1] - try - return diagonalize(s) - catch - error("Minreal only implemented for diagonalizable systems.") - end -end - -""" -`dsys = diagonalize(s::StateSpace, digits=12)` Diagonalizes the system such that the A-matrix is diagonal. The result is rounded to `digits` decimal points. -""" -function diagonalize(s::StateSpace, digits = 12) - r = x -> round(x,digits) - S,V = eig(s.A) - try - A = V\s.A*V .|> r - B = V\s.B .|> r - C = s.C*V .|> r - D = s.D .|> r - return ss(A,B,C,D) - catch e - error("System not diagonalizable", e) - end -end diff --git a/src/types/tf.jl b/src/types/tf.jl new file mode 100644 index 000000000..6079e61d4 --- /dev/null +++ b/src/types/tf.jl @@ -0,0 +1,50 @@ +function tf(num::AbstractVecOrMat{<:AbstractVector{T1}}, den::AbstractVecOrMat{<:AbstractVector{T2}}, Ts::Real=0.0) where {T1<:Number, T2<:Number} + # Validate input and output dimensions match + ny, nu = size(num, 1, 2) + if (ny, nu) != size(den, 1, 2) + error("num and den dimensions must match") + end + + T = promote_type(T1,T2) + matrix = Matrix{SisoRational{T}}(ny, nu) + for o=1:ny + for i=1:nu + matrix[o, i] = SisoRational{T}(Vector{T}(num[o, i]), Vector{T}(den[o, i])) + end + end + return TransferFunction{SisoRational{T}}(matrix, Ts) +end +function tf(num::AbstractVector{T1}, den::AbstractVector{T2}, Ts::Real=0.0) where {T1<:Number, T2<:Number} + T = promote_type(T1, T2) + return TransferFunction{SisoRational{T}}(fill(SisoRational{T}(num, den), 1, 1), Ts) +end +tf(num::Number, den::Vector, Ts::Real=0.0) = tf([num], den, Ts) + +# Cases for just static gain +function tf(D::AbstractArray{T}, Ts::Real=0.0) where {T<:Number} + ny, nu = size(D, 1, 2) + + matrix = Matrix{SisoRational{T}}(ny, nu) + for i in eachindex(D) + matrix[i] = SisoRational{T}([D[i]], [one(T)]) + end + return TransferFunction{SisoRational{T}}(matrix, Float64(Ts)) +end +tf(n::Number, Ts::Real=0; kwargs...) = tf([n], Ts; kwargs...) + +tf(sys::StateSpace) = convert(TransferFunction, sys) # NOTE: Would perhaps like to write TransferFunction{SisoRational}, but couldn't get this to work... + +function tf(G::TransferFunction{<:SisoTf{T}}) where {T<:Number} + convert(TransferFunction{SisoRational{T}}, G) +end + +# Function for creation of 's' or 'z' var +function tf(var::AbstractString) + var != "s" && error("var must be 's' for continuous time tf.") + return tf([1, 0], [1]) +end +function tf(var::AbstractString, Ts::Real) + var != "z" && error("var must be 'z' for discrete time tf.") + Ts == 0 && error("Ts must not be 0 for discrete time tf.") + return tf([1, 0], [1], Ts) +end diff --git a/src/types/tf2ss.jl b/src/types/tf2ss.jl deleted file mode 100644 index a513a85d0..000000000 --- a/src/types/tf2ss.jl +++ /dev/null @@ -1,174 +0,0 @@ -function Base.convert(::Type{StateSpace}, t::TransferFunction) - if !isproper(t) - error("System is improper, a state-space representation is impossible") - end - ny, nu = size(t) - mat = t.matrix - # TODO : These are added due to scoped for blocks, but is a hack. This - # could be much cleaner. - Ac = Bc = Cc = Dc = A = B = C = D = Array{Float64}(0, 0) - for i=1:nu - for j=1:ny - a, b, c, d = siso_tf_to_ss(mat[j, i]) - if j > 1 - # vcat - Ac = blkdiag(Ac, a) - Bc = vcat(Bc, b) - Cc = blkdiag(Cc, c) - Dc = vcat(Dc, d) - else - Ac, Bc, Cc, Dc = (a, b, c, d) - end - end - if i > 1 - # hcat - A = blkdiag(A, Ac) - B = blkdiag(B, Bc) - C = hcat(C, Cc) - D = hcat(D, Dc) - else - A, B, C, D = (Ac, Bc, Cc, Dc) - end - end - A, B, C = balance_statespace(A, B, C)[1:3] - return ss(A, B, C, D, t.Ts, inputnames=t.inputnames, outputnames=t.outputnames) -end - -Base.convert(::Type{StateSpace}, t::Real) = ss(t) - -Base.promote_rule{T<:SisoTf}(::Type{StateSpace}, ::Type{TransferFunction{T}}) = StateSpace -Base.promote_rule{T<:Real}(::Type{StateSpace}, ::Type{T}) = StateSpace - -siso_tf_to_ss(t::SisoTf) = siso_tf_to_ss(convert(SisoRational, t)) - -function siso_tf_to_ss(t::SisoRational) - t = normalize_tf(t) - tnum = num(t) - tden = den(t) - len = length(tden) - d = Array{Float64}(1, 1) - d[1] = tnum[1] - - if len==1 || tnum == zero(Poly{Float64}) - a = zeros(0, 0) - b = zeros(0, 1) - c = zeros(1, 0) - else - tden = tden[2:end] - a = [-tden' ; eye(len - 2, len - 1)] - b = eye(len - 1, 1) - c = tnum[2:len]' - d * tden[:]' - end - return float64mat(a), float64mat(b), float64mat(c), d -end - -function normalize_tf(t::SisoRational) - d = t.den[1] - return SisoTf(t.num/d, t.den/d) -end - - -""" -`A, B, C, T = balance_statespace{S}(A::Matrix{S}, B::Matrix{S}, C::Matrix{S}, perm::Bool=false)` - -`sys, T = balance_statespace(sys::StateSpace, perm::Bool=false)` - -Computes a balancing transformation `T` that attempts to scale the system so -that the row and column norms of [T*A/T T*B; C/T 0] are approximately equal. -If `perm=true`, the states in `A` are allowed to be reordered. - -This is not the same as finding a balanced realization with equal and diagonal observability and reachability gramians, see `balreal` -""" -function balance_statespace{S}(A::Matrix{S}, B::Matrix{S}, - C::Matrix{S}, perm::Bool=false) - nx = size(A, 1) - nu = size(B,2) - ny = size(C,1) - - # Compute the transformation matrix - mag_A = abs.(A) - mag_B = maximum([abs.(B) zeros(S, nx, 1)], 2) - mag_C = maximum([abs.(C); zeros(S, 1, nx)], 1) - T = balance_transform(mag_A, mag_B, mag_C, perm) - - # Perform the transformation - A = T*A/T - B = T*B - C = C/T - - return A, B, C, T -end - -function balance_statespace(sys::StateSpace, perm::Bool=false) - A, B, C, T = balance_statespace(sys.A,sys.B,sys.C, perm) - return ss(A,B,C,sys.D), T -end - -""" -`T = balance_transform{R}(A::Matrix{R}, B::Matrix{R}, C::Matrix{R}, perm::Bool=false)` - -`T = balance_transform(sys::StateSpace, perm::Bool=false) = balance_transform(A,B,C,perm)` - -Computes a balancing transformation `T` that attempts to scale the system so -that the row and column norms of [T*A/T T*B; C/T 0] are approximately equal. -If `perm=true`, the states in `A` are allowed to be reordered. - -This is not the same as finding a balanced realization with equal and diagonal observability and reachability gramians, see `balreal` -See also `balance_statespace`, `balance` -""" -function balance_transform{R}(A::Matrix{R}, B::Matrix{R}, C::Matrix{R}, perm::Bool=false) - nx = size(A, 1) - # Compute a scaling of the system matrix M - S = diag(balance([A B; C zeros(size(C*B))], false)[1]) - Sx = S[1:nx] - Sio = S[nx+1] - # Compute permutation of x (if requested) - pvec = perm ? balance(A, true)[2] * [1:nx;] : [1:nx;] - # Compute the transformation matrix - T = zeros(R, nx, nx) - T[pvec, :] = Sio * diagm(1./Sx) - return T -end - -balance_transform(sys::StateSpace, perm::Bool=false) = balance_transform(sys.A,sys.B,sys.C,perm) - - -@doc """`sys = ss2tf(s::StateSpace)`, ` sys = ss2tf(A, B, C, Ts = 0; inputnames = "", outputnames = "")` - -Convert a `StateSpace` realization to a `TransferFunction`""" -> -function ss2tf(s::StateSpace) - return ss2tf(s.A, s.B, s.C, s.D, s.Ts, inputnames=s.inputnames, outputnames=s.outputnames) -end - -function ss2tf(A, B, C, D, Ts = 0; inputnames = "", outputnames = "") - nu,ny = size(B,2),size(C,1) - ubernum = Matrix{Vector}(ny,nu) - uberden = Matrix{Vector}(ny,nu) - for i = 1:nu, j=1:ny - ubernum[j,i],uberden[j,i] = sisoss2tf(A, B[:,i], C[j,:]', D[j,i]) - end - tf(ubernum,uberden, Ts, inputnames=inputnames, outputnames=outputnames) -end - -function sisoss2tf(A, B, C, D) - charpolA = charpoly(A) - numP = charpoly(A-B*C) - charpolA + D*charpolA - denP = charpolA - return numP[1:length(numP)], denP[1:length(denP)] -end - -tf(sys::StateSpace) = ss2tf(sys) -zpk(sys::StateSpace) = zpk(ss2tf(sys)) - -function charpoly(A) - λ = eigvals(A); - p = reduce(*,ControlSystems.Poly([1.]), ControlSystems.Poly[ControlSystems.Poly([1, -λᵢ]) for λᵢ in λ]); - if maximum(imag.(p[:])./(1+abs.(real.(p[:])))) < sqrt(eps()) - for i = 1:length(p) - p[i] = real(p[i]) - end - else - error("Characteristic polynomial should be real") - end - p -end diff --git a/src/types/transferfunction.jl b/src/types/transferfunction.jl deleted file mode 100644 index 6913730d1..000000000 --- a/src/types/transferfunction.jl +++ /dev/null @@ -1,494 +0,0 @@ -abstract type SisoTf end -include("polys.jl") -include("sisotf.jl") -include("sisozpk.jl") -include("sisogeneralized.jl") -##################################################################### -## Data Type Declarations ## -##################################################################### - -type TransferFunction{S<:SisoTf} <: LTISystem - matrix::Matrix{S} - Ts::Float64 - nu::Int - ny::Int - inputnames::Vector{String} - outputnames::Vector{String} - function TransferFunction{T}(matrix::Matrix{T}, Ts::Float64, - inputnames::Vector{String}, outputnames::Vector{String}) where T<:SisoTf - # Validate size of input and output names - ny, nu = size(matrix) - if size(inputnames, 1) != nu - error("Must have same number of inputnames as inputs") - elseif size(outputnames, 1) != ny - error("Must have same number of outputnames as outputs") - end - # Validate sampling time - if Ts < 0 && Ts != -1 - error("Ts must be either a positive number, 0 - (continuous system), or -1 (unspecified)") - end - return new{T}(matrix, Ts, nu, ny, inputnames, outputnames) - end -end -TransferFunction{T<:SisoTf}(matrix::Matrix{T}, args...) = TransferFunction{T}(matrix, args...) - -+{T<:Real}(a::TransferFunction, b::AbstractVecOrMat{T}) = +(promote(a,b)...) - -Base.promote_rule{S<:SisoTf,T<:Real}(::Type{TransferFunction{S}}, ::Type{T}) = TransferFunction{S} -Base.promote_rule{S<:SisoTf,T<:Real}(::Type{TransferFunction{S}}, ::Union{Type{Array{T,2}},Type{Array{T,1}}}) = TransferFunction{S} - -Base.convert{T<:Real}(::Type{TransferFunction}, b::T) = tf([b]) -Base.convert{T<:Real}(::Type{TransferFunction{SisoRational}}, b::T) = tf(b) -Base.convert{T<:Real}(::Type{TransferFunction{SisoZpk}}, b::T) = zpk(b) -Base.convert{T<:Real}(::Type{TransferFunction{SisoGeneralized}}, b::T) = tfg(b) - -Base.convert(::Type{TransferFunction{SisoZpk}}, s::TransferFunction) = zpk(s) -Base.convert(::Type{TransferFunction{SisoRational}}, s::TransferFunction) = tf(s) -Base.convert(::Type{TransferFunction{SisoGeneralized}}, s::TransferFunction) = tfg(s) - -Base.promote_rule(::Type{TransferFunction{SisoRational}}, ::Type{TransferFunction{SisoZpk}}) = TransferFunction{SisoZpk} -Base.promote_rule{T<:SisoTf}(::Type{TransferFunction{T}}, ::Type{TransferFunction{SisoGeneralized}}) = TransferFunction{SisoGeneralized} -Base.promote_rule(::Type{SisoRational}, ::Type{SisoZpk}) = SisoZpk -Base.promote_rule{T<:SisoTf}(::Type{T}, ::Type{SisoGeneralized}) = SisoGeneralized - -function Base.convert{T<:Real}(::Type{TransferFunction}, b::VecOrMat{T}) - r = Array{TransferFunction,2}(size(b,2),1) - for j=1:size(b,2) - r[j] = vcat(map(k->convert(TransferFunction,k),b[:,j])...) - end - hcat(r...) -end - -function Base.convert(::Type{SisoZpk}, sys::SisoRational) - if length(sys.num) == 0 - return SisoZpk([],[],0) - elseif all(sys.den == zero(sys.den)) - error("Zero denominator, this should not be possible") - else - return SisoZpk(roots(sys.num),roots(sys.den),sys.num[1]/sys.den[1]) - end -end - -function Base.convert(::Type{SisoRational}, sys::SisoZpk) - num = prod(zp2polys(sys.z))*sys.k - den = prod(zp2polys(sys.p)) - return SisoRational(num, den) -end - -Base.convert(::Type{SisoGeneralized}, sys::SisoRational) = SisoGeneralized(sprint(print_compact, sys)) -Base.convert(::Type{SisoGeneralized}, sys::SisoZpk) = convert(SisoGeneralized, convert(SisoRational, sys)) - -Base.convert(::Type{SisoRational}, sys::SisoGeneralized) = SisoRational(sys.expr) -Base.convert(::Type{SisoZpk}, sys::SisoGeneralized) = convert(SisoZpk, SisoRational(sys.expr)) - -#Just default SisoTf to SisoRational -SisoTf(args...) = SisoRational(args...) -Base.convert(::Type{ControlSystems.SisoTf}, b::Real) = Base.convert(ControlSystems.SisoRational, b) -Base.zero(::Type{SisoTf}) = zero(SisoRational) -Base.zero(::SisoTf) = zero(SisoRational) - -tzero(sys::SisoTf) = error("tzero is not implemented for type $(typeof(sys))") -##################################################################### -## SisoTf Operations ## -##################################################################### - -#These make sure that the matrix operation below works as expected -#Base.convert(::Type{SisoTf}, b::Real) = Base.convert(SisoRational, b) -*{T<:SisoTf}(a::Array{T}, b::Real) = map(x->x*b,a) -*{T<:SisoTf}(b::Real, a::Array{T}) = map(x->x*b,a) -/{T<:SisoTf}(a::Array{T}, b::Real) = map(x->x/b,a) -+{T<:SisoTf}(a::Array{T}, b::Real) = map(x->x+b,a) -+{T<:SisoTf}(b::Real, a::Array{T}) = map(x->x+b,a) --{T<:SisoTf}(a::Array{T}, b::Real) = map(x->x-b,a) --{T<:SisoTf}(b::Real, a::Array{T}) = map(x->b-x,a) --{T<:SisoTf}(a::Array{T}) = map(x-> -x,a) - -#Operations with different types of Siso functions -*(a::SisoTf, b::SisoTf) = *(promote(a,b)...) -+(a::SisoTf, b::SisoTf) = +(promote(a,b)...) --(a::SisoTf, b::SisoTf) = -(promote(a,b)...) -#.*(a::SisoTf, b::SisoTf) = .*(promote(a,b)...) -#.+(a::SisoTf, b::SisoTf) = .+(promote(a,b)...) -#.-(a::SisoTf, b::SisoTf) = .+(promote(a,b)...) - -==(a::SisoTf, b::SisoTf) = ==(promote(a,b)...) -!=(a::SisoTf, b::SisoTf) = !(a==b) -isapprox(a::SisoTf, b::SisoTf; kwargs...) = isapprox(promote(a,b)...; kwargs...) - -# Promote_op types -Base.promote_op{T<:SisoTf}(::Any, ::Type{T}, ::Type{T}) = T - -##################################################################### -## Constructor Functions ## -##################################################################### - -@doc """ `sys = tf(num, den, Ts=0; kwargs...), sys = tf(gain, Ts=0; kwargs...)` - -Create transfer function as a fraction of polynomials: - -`sys = numerator/denominator` - -`num`: the coefficients of the numerator polynomial. Either scalar or vector to create SISO systems -or an array of vectors to create MIMO system. - -`den`: the coefficients of the denominator polynomial. Either vector to create SISO systems -or an array of vectors to create MIMO system. - -`Ts`: Sample time or `0` for continuous system. - -`kwargs`: `inputnames`, `outputnames`: Arrays of strings representing the inputs and outputs. - -Other uses: - -`tf(sys)`: Convert `sys` to `tf` form. - -`tf("s")`, `tf("z")`: Create the continous transferfunction `s`.""" -> -function tf{T<:Vector, S<:Vector}(num::VecOrMat{T}, den::VecOrMat{S}, Ts::Real=0; kwargs...) - # Validate input and output dimensions match - ny, nu = size(num, 1, 2) - if (ny, nu) != size(den, 1, 2) - error("num and den dimensions must match") - end - matrix = Array{SisoRational}(ny, nu) - for o=1:ny - for i=1:nu - matrix[o, i] = SisoRational(num[o, i], den[o, i]) - end - end - kvs = Dict(kwargs) - inputnames = validate_names(kvs, :inputnames, nu) - outputnames = validate_names(kvs, :outputnames, ny) - return TransferFunction(matrix, Float64(Ts), inputnames, outputnames) -end - -@doc """ `zpk(gain, Ts=0; kwargs...), zpk(num, den, k, Ts=0; kwargs...), zpk(sys)` - -Create transfer function on zero pole gain form. The numerator and denominator are represented by their poles and zeros. - -`sys = k*numerator/denominator` - -`num`: the roots of the numerator polynomial. Either scalar or vector to create SISO systems -or an array of vectors to create MIMO system. - -`den`: the roots of the denominator polynomial. Either vector to create SISO systems -or an array of vectors to create MIMO system. - -`k`: The gain of the system. Obs, this is not the same as `dcgain`. - -`Ts`: Sample time or `0` for continuous system. - -`kwargs`: `inputnames`, `outputnames`: Arrays of strings representing the inputs and outputs. - -Other uses: - -`tf(sys)`: Convert `sys` to `tf` form. - -`tf("s")`: Create the transferfunction `s`.""" -> -function zpk{T<:Vector,S<:Vector}(z::VecOrMat{T}, p::VecOrMat{S}, k::VecOrMat, Ts::Real=0; kwargs...) - # Validate input and output dimensions match - ny, nu = size(z, 1, 2) - if (ny, nu) != size(p, 1, 2) || (ny, nu) != size(k, 1, 2) - error("s, p, and k kdimensions must match") - end - matrix = Array{SisoZpk}(ny, nu) - for o=1:ny - for i=1:nu - matrix[o, i] = SisoZpk(z[o, i], p[o, i], k[o, i]) - end - end - kvs = Dict(kwargs) - inputnames = validate_names(kvs, :inputnames, nu) - outputnames = validate_names(kvs, :outputnames, ny) - return TransferFunction(matrix, Float64(Ts), inputnames, outputnames) -end - -function zpk(tf::TransferFunction) - oldmat = tf.matrix - matrix = Array{SisoZpk}(tf.ny, tf.nu) - for i in eachindex(oldmat) - matrix[i] = convert(SisoZpk, oldmat[i]) - end - return TransferFunction(matrix, tf.Ts, copy(tf.inputnames), copy(tf.outputnames)) -end - -function tf(tf::TransferFunction) - oldmat = tf.matrix - matrix = Array{SisoRational}(tf.ny, tf.nu) - for i in eachindex(oldmat) - matrix[i] = convert(SisoRational, oldmat[i]) - end - return TransferFunction(matrix, tf.Ts, copy(tf.inputnames), copy(tf.outputnames)) -end - -@doc """ `sys = tfg(tf::LTISystem), `tfg(s::AbstractString)`, `tfg(exp::Expr)`, `tfg(::Array)` - -Create generalized transfer function represented by an expression. The variable has to be `s`. - -Example: `tfg("1/exp(-sqrt(s))")`, `tfg(["1/exp(-sqrt(s))"), "1/(s+1)])`, `tfg(:(s+1))` - -Other uses: - -`tfg(sys)`: Convert `sys` to `tfg` form. -""" -> -function tfg(tf::TransferFunction) - oldmat = tf.matrix - matrix = Array{SisoGeneralized}(tf.ny, tf.nu) - for i in eachindex(oldmat) - matrix[i] = convert(SisoGeneralized, oldmat[i]) - end - return TransferFunction(matrix, tf.Ts, copy(tf.inputnames), copy(tf.outputnames)) -end - -zpk(sys::TransferFunction{SisoGeneralized}) = zpk(tf(sys)) - -tf(num::Vector, den::Vector, Ts::Real=0; kwargs...) = - tf(reshape(Vector[num], 1, 1), reshape(Vector[den], 1, 1), Ts; kwargs...) - -tf(num::Real, den::Vector, Ts::Real=0; kwargs...) = tf([num], den, Ts; kwargs...) - -zpk(z::Vector, p::Vector, k::Real, Ts::Real=0; kwargs...) = - zpk(reshape(Vector[z], 1, 1), reshape(Vector[p], 1, 1), reshape([k],1,1), Ts; kwargs...) - -# Function for creation of static gain -function tf(gain::Array, Ts::Real=0; kwargs...) - ny, nu = size(gain, 1, 2) - matrix = Array{SisoRational}(ny, nu) - for i in eachindex(gain) - matrix[i] = SisoRational([gain[i]], [1]) - end - kvs = Dict(kwargs) - inputnames = validate_names(kvs, :inputnames, nu) - outputnames = validate_names(kvs, :outputnames, ny) - return TransferFunction(matrix, Float64(Ts), inputnames, outputnames) -end - -function zpk(gain::Array, Ts::Real=0; kwargs...) - ny, nu = size(gain, 1, 2) - matrix = Array{SisoZpk}(ny, nu) - for o=1:ny - for i=1:nu - matrix[o, i] = SisoZpk([],[], gain[o, i]) - end - end - kvs = Dict(kwargs) - inputnames = validate_names(kvs, :inputnames, nu) - outputnames = validate_names(kvs, :outputnames, ny) - return TransferFunction(matrix, Float64(Ts), inputnames, outputnames) -end - -tf(gain::Real, Ts::Real=0; kwargs...) = tf([gain], Ts; kwargs...) -zpk(k::Real, Ts::Real=0; kwargs...) = zpk([], [], k, Ts; kwargs...) - -# Function for creation of 's' or 'z' var -function tf(var::AbstractString) - var != "s" && error("var must be 's' for continuous time tf.") - return tf([1, 0], [1]) -end -function tf(var::AbstractString, Ts::Real) - var != "z" && error("var must be 'z' for discrete time tf.") - Ts == 0 && error("Ts must not be 0 for discrete time tf.") - return tf([1, 0], [1], Ts) -end - -zpk(var::AbstractString) = zpk(tf(var)) -zpk(var::AbstractString, Ts::Real) = zpk(tf(var, Ts)) - -function tfg(systems::Array, Ts::Real=0; kwargs...) - ny, nu = size(systems, 1, 2) - matrix = Array{SisoGeneralized}(ny, nu) - for o=1:ny - for i=1:nu - matrix[o, i] = SisoGeneralized(systems[o, i]) - end - end - kvs = Dict(kwargs) - inputnames = validate_names(kvs, :inputnames, nu) - outputnames = validate_names(kvs, :outputnames, ny) - return TransferFunction(matrix, Float64(Ts), inputnames, outputnames) -end - -tfg(var::Union{AbstractString,ExprLike}, Ts=0; kwargs...) = tfg([var], Ts; kwargs...) -##################################################################### -## Misc. Functions ## -##################################################################### - -## INDEXING ## -Base.ndims(::TransferFunction) = 2 -Base.size(t::TransferFunction) = (t.ny, t.nu) -Base.size(t::TransferFunction, d) = d <= 2 ? size(t)[d] : 1 - -function Base.getindex(t::TransferFunction, inds...) - if size(inds, 1) != 2 - error("Must specify 2 indices to index TransferFunction model") - end - rows, cols = ControlSystems.index2range(inds...) - T = eltype(t.matrix) - mat = Array{T}(length(rows), length(cols)) - mat[:, :] = t.matrix[rows, cols] - innames = String[t.inputnames[i] for i in cols] - outnames = String[t.outputnames[i] for i in rows] - return TransferFunction(mat, t.Ts, innames, outnames) -end - -function Base.copy(t::TransferFunction) - matrix = copy(t.matrix) - inputnames = copy(t.inputnames) - outputnames = copy(t.outputnames) - return TransferFunction(matrix, t.Ts, inputnames, outputnames) -end - -@doc """`tf = minreal(tf::TransferFunction, eps=sqrt(eps()))` - -Create a minimial representation of each transfer function in `tf` by cancelling poles and zeros """ -> -function minreal(t::TransferFunction, eps::Real=sqrt(eps())) - matrix = similar(t.matrix) - for i = eachindex(t.matrix) - matrix[i] = minreal(t.matrix[i], eps) - end - return TransferFunction(matrix, t.Ts, copy(t.inputnames), copy(t.outputnames)) -end -##################################################################### -## Math Operators ## -##################################################################### - -## EQUALITY ## -function ==(t1::TransferFunction, t2::TransferFunction) - fields = [:Ts, :ny, :nu, :inputnames, :outputnames, :matrix] - for field in fields - if getfield(t1, field) != getfield(t2, field) - return false - end - end - return true -end - -## Approximate ## -function isapprox(t1::TransferFunction, t2::TransferFunction; kwargs...) - t1, t2 = promote(t1, t2) - fieldsApprox = [:Ts, :ny, :nu, :matrix] - fieldsEqual = [:inputnames, :outputnames] - for field in fieldsApprox - if !(isapprox(getfield(t1, field), getfield(t2, field); kwargs...)) - return false - end - end - for field in fieldsEqual - if getfield(t1, field) != getfield(t2, field) - return false - end - end - return true -end - -function isapprox{T<:SisoTf, S<:SisoTf}(t1::Array{T}, t2::Array{S}; kwargs...) - reduce(&, [isapprox(t1[i], t2[i]; kwargs...) for i in eachindex(t1)]) -end - -## ADDITION ## -function +(t1::TransferFunction, t2::TransferFunction) - if size(t1) != size(t2) - error("Systems have different shapes.") - elseif t1.Ts != t2.Ts - error("Sampling time mismatch") - end - - # Naming strategy: If only one sys is named, use that. If the names are the - # same, use them. If the names conflict, then they are ignored, and the - # default "" is used. - if all(t1.inputnames .== "") - inputnames = t2.inputnames - elseif all(t2.inputnames .== "") || (t1.inputnames == t2.inputnames) - inputnames = t1.inputnames - else - inputnames = fill(String(""),t1.ny) - end - if all(t1.outputnames .== "") - outputnames = t2.outputnames - elseif all(t2.outputnames .== "") || (t1.outputnames == t2.outputnames) - outputnames = t1.outputnames - else - outputnames = fill(String(""),t1.nu) - end - t1, t2 = promote(t1, t2) - matrix = t1.matrix + t2.matrix - return TransferFunction(matrix, t1.Ts, inputnames, outputnames) -end - -+(t::TransferFunction, n::Real) = TransferFunction(t.matrix + n, t.Ts, - t.inputnames, t.outputnames) -+(n::Real, t::TransferFunction) = +(t, n) - -## SUBTRACTION ## --(n::Real, t::TransferFunction) = TransferFunction(n - t.matrix, t.Ts, - t.inputnames, t.outputnames) --(t1::TransferFunction, t2::TransferFunction) = +(t1, -t2) --(t::TransferFunction, n::Real) = +(t, -n) - -## NEGATION ## --(t::TransferFunction) = TransferFunction(-t.matrix, t.Ts, t.inputnames, - t.outputnames) - -## MULTIPLICATION ## -function *(t1::TransferFunction, t2::TransferFunction) - # Note: t1*t2 = y <- t1 <- t2 <- u - if t1.nu != t2.ny - error("t1*t2: t1 must have same number of inputs as t2 has outputs") - elseif t1.Ts != t2.Ts - error("Sampling time mismatch") - end - matrix = t1.matrix*t2.matrix - return TransferFunction(matrix, t1.Ts, t2.inputnames, t1.outputnames) -end - -*(t::TransferFunction, n::Real) = TransferFunction(n*t.matrix, t.Ts, - t.inputnames, t.outputnames) -*(n::Real, t::TransferFunction) = *(t, n) - -## DIVISION ## -function /(n::Real, t::TransferFunction) - if issiso(t) - matrix = reshape([n/t.matrix[1,1]], 1, 1) - else - error("MIMO TransferFunction inversion isn't implemented yet") - end - return TransferFunction(matrix, t.Ts, t.outputnames, t.inputnames) -end -/(t::TransferFunction, n::Real) = t*(1/n) -/(t1::TransferFunction, t2::TransferFunction) = t1*(1/t2) - -##################################################################### -## Display Functions ## -##################################################################### - -Base.print(io::IO, t::TransferFunction) = show(io, t) - -function Base.show(io::IO, t::TransferFunction) - # Compose the name vectors - inputs = format_names(t.inputnames, "Input ", "?") - outputs = format_names(t.outputnames, "Output ", "?") - println(io, "TransferFunction:") - tftype = iscontinuous(t) ? :s : :z - for i=1:t.nu - for o=1:t.ny - if !issiso(t) - println(io, inputs[i], " to ", outputs[o]) - end - print_siso(io, t.matrix[o, i], tftype) - if !(i == t.nu && o == t.ny) - print(io, "\n") - end - end - end - if iscontinuous(t) - print(io, "\nContinuous-time transfer function model") - else - print(io, "\nSample Time: ") - if t.Ts > 0 - print(io, t.Ts, " (seconds)") - elseif t.Ts == -1 - print(io, "unspecified") - end - print(io, "\nDiscrete-time transfer function model") - end -end diff --git a/src/types/zpk.jl b/src/types/zpk.jl new file mode 100644 index 000000000..bec9a7227 --- /dev/null +++ b/src/types/zpk.jl @@ -0,0 +1,78 @@ +@doc """ `zpk(gain, Ts=0; kwargs...), zpk(num, den, k, Ts=0; kwargs...), zpk(sys)` + +Create transfer function on zero pole gain form. The numerator and denominator are represented by their poles and zeros. + +`sys = k*numerator/denominator` + +`num`: the roots of the numerator polynomial. Either scalar or vector to create SISO systems +or an array of vectors to create MIMO system. + +`den`: the roots of the denominator polynomial. Either vector to create SISO systems +or an array of vectors to create MIMO system. + +`k`: The gain of the system. Obs, this is not the same as `dcgain`. + +`Ts`: Sample time or `0` for continuous system. + +Other uses: + +`tf(sys)`: Convert `sys` to `tf` form. + +`tf("s")`: Create the transferfunction `s`.""" -> +function zpk(z::VecOrMat{<:AbstractVector{TZ}}, p::VecOrMat{<:AbstractVector{TP}}, k::VecOrMat{T0}, Ts::Real=0.0) where {T0<:Number, TZ<:Number, TP<:Number} + # Validate input and output dimensions match + if !(size(z, 1, 2) == size(p, 1, 2) == size(k, 1, 2)) + error("Dimensions for s, p, and k must match") + end + + TR = promote_type(T0,TZ,TP) + T = promote_type(T0, real(TR)) # Ensuring this avoids many problems, e.g., with SisoZpk{Int64,Complex128} + + matrix = Array{SisoZpk{T, TR}}(size(z,1), size(z,2)) # TODO: make this nicer + for o=1:size(z,1) + for i=1:size(z,2) + matrix[o, i] = SisoZpk{T,TR}(z[o, i], p[o, i], k[o, i]) + end + end + return TransferFunction{SisoZpk{T,TR}}(matrix, Float64(Ts)) +end +function zpk(z::AbstractVector{TZ}, p::AbstractVector{TP}, k::T, Ts::Real=0.0) where {T<:Number, TZ<:Number, TP<:Number} + return zpk(fill(z, 1, 1), fill(p, 1, 1), fill(k, 1, 1), Ts) +end + +function zpk(gain::Matrix{T}, Ts::Real=0; kwargs...) where {T <: Number} + TR = promote_type(Complex128,T) + ny, nu = size(gain) + matrix = [SisoZpk{T, TR}(TR[],TR[], gain[o, i]) for o=1:ny, i=1:nu] + return TransferFunction{SisoZpk{T,TR}}(matrix, Ts) +end + +function zpk(z::AbstractVector, p::AbstractVector, k::T) where {T<:Number} # To be able to send in empty vectors [] of type Any + if eltype(z) == Any && eltype(p) == Any + @assert z == [] + @assert p == [] + return zpk(T[], T[], k) + elseif eltype(z) == Any + @assert z == [] + TR = eltype(p) + return zpk(TR[], p, k) + elseif eltype(p) == Any + @assert p == [] + TR = eltype(z) + return zpk(z, TR[], k) + else + error("Non numeric vectors must be empty.") + end +end + +zpk(k::Real, Ts::Real=0) = zpk(eltype(k)[], eltype(k)[], k, Ts) + +zpk(sys::StateSpace) = zpk(tf(sys)) # FIXME: probably better with direct conversion + +function zpk(G::TransferFunction{S}) where {T0, S<:SisoTf{T0}} + T = typeof(one(T0)/one(T0)) + convert(TransferFunction{SisoZpk{T, complex(T)}}, G) +end + +zpk(var::AbstractString) = zpk(tf(var)) +zpk(var::AbstractString, Ts::Real) = zpk(tf(var, Ts)) diff --git a/src/utilities.jl b/src/utilities.jl index 0266fa0f5..c83a7d630 100644 --- a/src/utilities.jl +++ b/src/utilities.jl @@ -1,87 +1,88 @@ -# Misfit functions +# Are only really needed for cases when we accept general LTISystem +# we should either use them consistently, with a good definition, or remove them +numeric_type(::Type{SisoRational{T}}) where T = T +numeric_type(::Type{<:SisoZpk{T}}) where T = T +numeric_type(sys::SisoTf) = numeric_type(typeof(sys)) -## Predicates ## +numeric_type(::Type{TransferFunction{S}}) where S = numeric_type(S) +numeric_type(::Type{<:StateSpace{T}}) where T = T +numeric_type(sys::LTISystem) = numeric_type(typeof(sys)) -@doc """`isstable(sys)` -Returns `true` if `sys` is stable, else returns `false`.""" -> -function isstable(sys::LTISystem) - if sys.Ts == 0 - if any(real.(pole(sys)).>=0) - return false - end - else - if any(abs.(pole(sys)).>=1) - return false - end - end - return true -end +to_matrix(T, A::Vector) = Matrix{T}(reshape(A, length(A), 1)) +to_matrix(T, A::AbstractMatrix) = Matrix{T}(A) +to_matrix(T, A::Number) = fill(T(A), 1, 1) -@doc """`iscontinuous(sys)` -Returns `true` if `sys` is continuous, else returns `false`.""" -> -function iscontinuous(sys::LTISystem) - return sys.Ts == 0 -end -@doc """`issiso(sys)` -Returns `true` if `sys` is SISO, else returns `false`.""" -> -function issiso(sys::LTISystem) - return sys.ny == 1 && sys.nu == 1 -end +# 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 +# returned by LAPACK routines for eigenvalues. +function roots2real_poly_factors(roots::Vector{cT}) where cT <: Number + T = real(cT) + poly_factors = Vector{Poly{T}}(0) + + for k=1:length(roots) + r = roots[k] + + if isreal(r) + push!(poly_factors,Poly{T}([-real(r),1])) + else + if imag(r) < 0 # This roots was handled in the previous iteration # TODO: Fix better error handling + continue + end -@doc """`isproper(tf)` + if k == length(roots) || r != conj(roots[k+1]) + throw(AssertionError("Found pole without matching conjugate.")) + end -Returns `true` if the `TransferFunction` is proper. This means that order(den) -\>= order(num))""" -> -function isproper(t::TransferFunction) - for s in t.matrix - if length(num(s)) > length(den(s)) - return false + push!(poly_factors,Poly{T}([real(r)^2+imag(r)^2, -2*real(r), 1])) + # k += 1 # Skip one iteration in the loop end end - return true -end - -## Helpers (not public) ## -# Convert the argument to a Matrix{Float64} -function float64mat(A::Vector) - A = reshape(A, size(A, 1), 1) - return float64mat(A) + return poly_factors end -float64mat(A::Matrix) = map(Float64,A) -float64mat(A::Matrix{Float64}) = A -float64mat(A::Real) = float64mat([A]) - -# Ensures the metadata for an LTISystem is valid -function validate_names(kwargs, key, n) - names = get(kwargs, key, "") - if names == "" - return fill(String(""), n) - elseif isa(names, Vector) && eltype(names) <: AbstractString - return String[names[i] for i = 1:n] - elseif isa(names, AbstractString) - return String[names * "$i" for i = 1:n] - else - error("$key must be of type `AbstractString` or Vector{AbstractString}") - end +# This function should hande both Complex as well as symbolic types +function roots2poly_factors(roots::Vector{T}) where T <: Number + return [Poly{T}([-r, 1]) for r in roots] end -# Format the metadata for printing -function format_names(names::Vector{String}, default::AbstractString, unknown::AbstractString) - n = size(names, 1) - if all(names .== "") - return String[default * string(i) for i=1:n] - else - for (i, n) in enumerate(names) - names[i] = String((n == "") ? unknown : n) + +""" Typically LAPACK returns a vector with, e.g., eigenvalues to a real matrix, + they are paired up in exact conjugate pairs. This fact is used in some places + for working with zpk representations of LTI systems. eigvals(A) returns a + on this form, however, for generalized eigenvalues there is a small numerical + discrepancy, which breaks some functions. This function fixes small + discrepancies in the conjugate pairs.""" +function _fix_conjugate_pairs!(v::AbstractVector{<:Complex}) + k = 1 + while k <= length(v) - 1 + if isreal(v[k]) + # Do nothing + else + if isapprox(v[k], conj(v[k+1]), rtol=1e-15) + z = (v[k] + conj(v[k+1]))/2 + v[k] = z + v[k+1] = conj(z) + k += 1 + end end - return names + k += 1 end end +function _fix_conjugate_pairs!(v::AbstractVector{<:Real}) + nothing +end + +# Should probably try to get rif of this function... +poly2vec(p::Poly) = p.a[1:end] + function unwrap!(M::Array, dim=1) alldims(i) = [ n==dim ? i : (1:size(M,n)) for n in 1:ndims(M) ] @@ -100,12 +101,6 @@ end unwrap(m::AbstractArray, args...) = unwrap!(collect(m), args...) unwrap(x::Number) = x - -numpoly(G::TransferFunction) = map(numpoly, G.matrix) -denpoly(G::TransferFunction) = map(denpoly, G.matrix) - -poly2vec(p::Poly) = p.a[1:end] - """ outs = index2range(ind1, ind2) Helper function to convert indexes with scalars to ranges. Used to avoid dropping dimensions diff --git a/test/REQUIRE b/test/REQUIRE index ea5c703f4..e27a56f22 100644 --- a/test/REQUIRE +++ b/test/REQUIRE @@ -1,6 +1,6 @@ julia 0.6 -VisualRegressionTests -GR -Plots +# VisualRegressionTests +# GR +# Plots FactCheck -Gtk +# Gtk diff --git a/test/runtests.jl b/test/runtests.jl index 08d2bce75..96deadd9c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,28 +4,29 @@ import Base.isapprox include("framework.jl") my_tests = ["test_statespace", - "test_lqg", "test_transferfunction", - "test_generalizedtf", "test_zpk", - "test_analysis", + "test_promotion", "test_connections", "test_discrete", + "test_conversion", + "test_complex", "test_linalg", "test_simplification", "test_freqresp", - "test_synthesis", - "test_matrix_comps", "test_timeresp", - "test_conversion"] + "test_analysis", + "test_matrix_comps", + "test_lqg", + "test_synthesis"] -try - Pkg.installed("ControlExamplePlots") - push!(my_tests, "test_plots") -catch - warn("The unregistered package ControlExamplePlots is currently needed to test plots, install using: - Pkg.clone(\"https://github.com/JuliaControl/ControlExamplePlots.jl.git\")") -end +# try +# Pkg.installed("ControlExamplePlots") +# push!(my_tests, "test_plots") +# catch +# warn("The unregistered package ControlExamplePlots is currently needed to test plots, install using: +# Pkg.clone(\"https://github.com/JuliaControl/ControlExamplePlots.jl.git\")") +# end run_tests(my_tests) diff --git a/test/test_analysis.jl b/test/test_analysis.jl index 55a4aa2b5..499232034 100644 --- a/test/test_analysis.jl +++ b/test/test_analysis.jl @@ -95,12 +95,12 @@ ex_11 = ss(A, B, C, D) # Test for multiple zeros, siso tf sys = s*(s + 1)*(s^2 + 1)*(s - 3)/((s + 1)*(s + 4)*(s - 4)) -@test tzero(sys) ≈ [-1.0, -im, im, 3.0, 0.0] +@test tzero(sys) ≈ [3.0, -1.0, im, -im, 0.0] ## POLE ## -@test pole(sys) ≈ [-1.0, 4.0, -4.0] -@test pole([sys sys]) ≈ [-1.0, 4.0, -4.0, -1.0, 4.0, -4.0] -@test pole(ex_11) ≈ eig(ex_11.A)[1] +@test pole(sys) ≈ [4.0, -4.0, -1.0] +@test_broken pole([sys sys]) ≈ [4.0, -4.0, -1.0] # Issue #81 +@test pole(ex_11) ≈ eigvals(ex_11.A) poles = [-3.383889568918823 + 0.000000000000000im -2.199935841931115 + 0.000000000000000im @@ -132,19 +132,17 @@ sol_p = vecarray(Complex128, 2, 2, Complex128[], Complex128[-0.5 - 3.12249899919 Complex128[-5.0 + 0.0im], Complex128[-6.0 + 0.0im]) sol_k = [0.0 3.0; 1.0 2.0] z, p, k = zpkdata(H) -@test_array_vecs_eps z sol_z 2*eps(Complex128) -@test_array_vecs_eps p sol_p 2*eps(Complex128) +@test_array_vecs_eps z sol_z 2*eps(Float64) +@test_array_vecs_eps p sol_p 2*eps(Float64) @test k == sol_k z, p, k = zpkdata(G) -@test_array_vecs_eps z sol_z 10*eps(Complex128) -@test_array_vecs_eps p sol_p 10*eps(Complex128) +@test_array_vecs_eps z sol_z 10*eps(Float64) +@test_array_vecs_eps p sol_p 10*eps(Float64) @test k == sol_k ## GAIN ## #Gain is confusing when referring to zpkdata. Test dcgain instead @test [dcgain(H[1, 1]) dcgain(H[1, 2]); dcgain(H[2, 1]) dcgain(H[2, 2])] ≈ [0 0; 0.2 1/3] @test [dcgain(G[1, 1]) dcgain(G[1, 2]); dcgain(G[2, 1]) dcgain(G[2, 2])] ≈ [0 0; 0.2 1/3] -@test_throws ErrorException dcgain(H) -@test_throws ErrorException dcgain(G) ## MARKOVPARAM ## @test markovparam(G, 0) == [0.0 0.0; 1.0 0.0] @@ -159,19 +157,27 @@ z, p, k = zpkdata(G) ## DAMPREPORT ## @test sprint(dampreport, sys) == ( -"| Pole | Damping | Frequency | Time Constant |\n"* -"| | Ratio | (rad/sec) | (sec) |\n"* -"+---------------+---------------+---------------+---------------+\n"* -"| -1.000e+00 | 1.000e+00 | 1.000e+00 | 1.000e+00 |\n"* -"| 4.000e+00 | -1.000e+00 | 4.000e+00 | -2.500e-01 |\n"* -"| -4.000e+00 | 1.000e+00 | 4.000e+00 | 2.500e-01 |\n") + "| Pole | Damping | Frequency | Time Constant |\n"* + "| | Ratio | (rad/sec) | (sec) |\n"* + "+---------------+---------------+---------------+---------------+\n"* + "| -1.000e+00 | 1.000e+00 | 1.000e+00 | 1.000e+00 |\n"* + "| 4.000e+00 | -1.000e+00 | 4.000e+00 | -2.500e-01 |\n"* + "| -4.000e+00 | 1.000e+00 | 4.000e+00 | 2.500e-01 |\n") @test sprint(dampreport, ex_11) == ( -"| Pole | Damping | Frequency | Time Constant |\n"* -"| | Ratio | (rad/sec) | (sec) |\n"* -"+---------------+---------------+---------------+---------------+\n"* -"| -1.000e+00 | 1.000e+00 | 1.000e+00 | 1.000e+00 |\n"* -"| 1.000e+00 | -1.000e+00 | 1.000e+00 | -1.000e+00 |\n"* -"| 2.000e+00 | -1.000e+00 | 2.000e+00 | -5.000e-01 |\n"* -"| -2.000e+00 | 1.000e+00 | 2.000e+00 | 5.000e-01 |\n"* -"| 3.000e+00 | -1.000e+00 | 3.000e+00 | -3.333e-01 |\n") + "| Pole | Damping | Frequency | Time Constant |\n"* + "| | Ratio | (rad/sec) | (sec) |\n"* + "+---------------+---------------+---------------+---------------+\n"* + "| -1.000e+00 | 1.000e+00 | 1.000e+00 | 1.000e+00 |\n"* + "| 1.000e+00 | -1.000e+00 | 1.000e+00 | -1.000e+00 |\n"* + "| 2.000e+00 | -1.000e+00 | 2.000e+00 | -5.000e-01 |\n"* + "| -2.000e+00 | 1.000e+00 | 2.000e+00 | 5.000e-01 |\n"* + "| 3.000e+00 | -1.000e+00 | 3.000e+00 | -3.333e-01 |\n") + + +# Example 5.5 from http://www.control.lth.se/media/Education/EngineeringProgram/FRTN10/2017/e05_both.pdf +G = [1/(s+2) -1/(s+2); 1/(s+2) (s+1)/(s+2)] +@test_broken length(pole(G)) == 1 +@test length(tzero(G)) == 1 +@test_broken size(minreal(ss(G)).A) == fill(2.0, 1, 1) + end diff --git a/test/test_complex.jl b/test/test_complex.jl new file mode 100644 index 000000000..146bd2bc8 --- /dev/null +++ b/test/test_complex.jl @@ -0,0 +1,35 @@ +@testset "test_complex_systems" begin + +C_1 = zpk([], [-1+im], 1.0+1im) +C_2 = zpk([-1+im], [], 1.0+1im) + +# Basic arittmetic +@test C_1 + C_1 == zpk([], [-1+im], 2+2im) +@test (1.0+im)*C_1 == zpk([], [-1+im], 2.0im) +@test C_2 + C_2 == zpk([-1+im], [], 2+2im) + +@test C_1 * C_1 == zpk([], [-1+im,-1+im], 2.0im) +@test C_2 * C_2 ≈ zpk([-1+im,-1+im], [], 2.0im) + +@test im*ss(1) == ss(im) + + +@test pole(zpk([], [-1+im,-1+im,0], 2.0im)) == [-1+im,-1+im,0] +@test tzero(zpk([-1+im,-1+im,0], [-2], 2.0im)) == [-1+im,-1+im,0] + +@test zpk( tf([1.0, 1+im], [1.0, 2+im]) ) == zpk( [-1-im], [-2-im], 1.0+0im) + +@test minreal(zpk([-1+im], [-1+im,-1+im],1+0im)) == zpk([], [-1+im],1+0im) +@test minreal(zpk([-1+im, -1+im], [-1+im],1+1im)) == zpk([-1+im], [], 1+1im) + + +@test_throws AssertionError zpk([-1+im], [-1+im,-1+im],1) # Given the type of k this should be a real-coefficient system, but poles and zeros don't come in conjugate pairs + +@test zpk([-2+im], [-1+im],1+0im)*zpk([], [-1+im],1+0im) == zpk([-2+im], [-1+im, -1+im], 1+0im) +@test zpk([], [-2], 2) + zpk([], [-1], 1) == zpk([-4/3], [-2,-1], 3) + +@test tf(zpk([-2+im], [-1+im],1+0im)) == tf([1, 2-im], [1, 1-im]) + +@test 1 / ( tf("s") + 1 + im ) == tf([1], [1, 1+im]) + +end diff --git a/test/test_connections.jl b/test/test_connections.jl index 3414169fe..ff6c7d7a2 100644 --- a/test/test_connections.jl +++ b/test/test_connections.jl @@ -52,7 +52,8 @@ Dtf_111 = tf([1, 2], [1, 5], 0.005) Dtf_211 = tf([1, 2, 3], [1, 8, 15], 0.005) Dtf_212 = tf(vecarray(2, 1, [1, 2, 3], [1, 2]), vecarray(2, 1, [1, 8, 15], [1, 8, 15]), 0.005) Dtf_221 = tf(vecarray(1, 2, [1, 2, 3], [1, 2]), vecarray(1, 2, [1, 8, 15], [1, 8, 15]), 0.005) -Dtf_222 = [Dtf_221; Dtf_221]; Dtf_222.Ts = 0.005 +Dtf_222 = [Dtf_221; Dtf_221]; +Dtf_222.Ts = 0.005 Dtf_022 = tf(4*eye(2), 0.005) s = tf("s") @@ -91,6 +92,30 @@ s = tf("s") @test [D_111; Dtf_212] == [D_111; ss(Dtf_212)] @test append(D_111, Dtf_211) == append(D_111, ss(Dtf_211)) + +# hcat and vcat for StateSpace and Matrix +A = [-1.1 -1.2; -1.3 -1.4] +B = [1 2; 3 4] +C = [5 6; 7 8] +D = [1 0; 0 1] +P = ss(A, B, C, D) +@test [P fill(2.5, 2, 1)] == ss(A, [B fill(0, 2, 1)], C, [D fill(2.5, 2, 1)]) +@test [fill(2.5, 2, 1) P] == ss(A, [fill(0, 2, 1) B], C, [fill(2.5, 2, 1) D]) +@test [P; fill(2.5, 1, 2)] == ss(A, B, [C; fill(0, 1, 2)], [D; fill(2.5, 1, 2)]) +@test [fill(2.5, 1, 2); P] == ss(A, B, [fill(0, 1, 2); C], [fill(2.5, 1, 2); D]) + +# hcat and vcat for StateSpace and Number +P = ss(-1.0, 2.0, 3.0, 4.0) +@test [P 2.5] == ss(-1.0, [2.0 0.0], 3.0, [4.0 2.5]) +@test [2.5 P] == ss(-1.0, [0.0 2.0], 3.0, [2.5 4.0]) +@test [P; 2.5] == ss(-1.0, 2.0, [3.0; 0.0], [4.0; 2.5]) +@test [2.5; P] == ss(-1.0, 2.0, [0.0; 3.0], [2.5; 4.0]) + +@test [2.5 P 3.5] == ss(-1.0, [0.0 2.0 0.0], 3.0, [2.5 4.0 3.5]) +@test [2.5; P; 3.5] == ss(-1.0, 2.0, [0.0; 3.0; 0.0], [2.5; 4.0; 3.5]) + + + # Combination tfRational and sisoZpk Czpk_111 = zpk([-2],[-5],1) Czpk_211 = zpk([-1+sqrt(2)im,-1-sqrt(2)im], [-5,-3],1) @@ -100,15 +125,19 @@ Czpk_222 = [Czpk_221; Czpk_221] Czpk_022 = [zpk([],[],4) 0; 0 zpk([],[],4)] #Make sure that we get a vector -arr = Array{TransferFunction{ControlSystems.SisoZpk},1}(2) +arr = Array{typeof(zpk(tf(1))),1}(2) arr[1] = zpk(tf(1)); arr[2] = zpk(2); @test [tf(1), zpk(2)] == arr -arr2 = Array{TransferFunction{ControlSystems.SisoRational},1}(2) +arr2 = Array{typeof(tf(1, 0.1)),1}(2) arr2[1] = tf(1, 0.1); arr2[2] = tf(2, 0.1); @test [tf(1, 0.1), tf(2, 0.1)] == arr2 -arr3 = Array{StateSpace,1}(3) -arr3[1] = ss(0); arr3[2] = ss(1); arr3[3] = ss(2) -@test [0, zpk(1), ss(2)] == arr3 +arr3 = Array{typeof(ss(0.)),1}(3) +arr3[1] = ss(0.); arr3[2] = ss(1.); arr3[3] = ss(2.) +@test [0., zpk(1), ss(2.)] == arr3 + +arr4 = Array{typeof(ss(0)),1}(3) +arr4[1] = ss(0); arr4[2] = ss(1); arr4[3] = ss(2) +@test [0., zpk(1), ss(2)] == arr4 @test Czpk_111 ≈ Ctf_111 @test Czpk_211 ≈ Ctf_211 @@ -121,4 +150,10 @@ arr3[1] = ss(0); arr3[2] = ss(1); arr3[3] = ss(2) #This might fail depending on if minreal is used or not @test (Czpk_211+1) ≈ (Ctf_211+1) + + + +@test_broken [D_111 1.0] # Concatenation of discrete system with constant +@test_broken [D_222 fill(1.5, 2, 2)] # Concatenation of discrete system with matrix + end diff --git a/test/test_conversion.jl b/test/test_conversion.jl index a84714b0d..a788bcce7 100644 --- a/test/test_conversion.jl +++ b/test/test_conversion.jl @@ -1,18 +1,122 @@ @testset "test_conversion" begin -f11 = tf(ss([-1 0;1 1],[1;0],[1 1],0)) -f12 = tf([1,0],[1,0,-1]) +# Easy second order system +sys1 = ss([-1 0;1 1],[1;0],[1 1],0) +G1 = tf([1,0],[1,0,-1]) +H1 = zpk([0.0], [-1, 1.0], 1.0) -#Direct term was previously missing -f21 = tf(ss([-1 0;1 1],[1;0],[1 1],1)) -f22 = tf([1,1,-1],[1,0,-1]) +w = 10.0 .^ range(-2,2,50) + +@test ss(sys1) == sys1 +@test freqresp(ss(G1), w) ≈ freqresp(sys1, w) rtol=1e-15 +@test freqresp(ss(H1), w) ≈ freqresp(sys1, w) rtol=1e-15 + +@test tf(sys1) ≈ G1 rtol=1e-15 +@test tf(G1) == G1 +@test tf(H1) == G1 + +@test zpk(sys1) ≈ H1 rtol=1e-15 +@test zpk(G1) == H1 +@test zpk(H1) == H1 + +# With direct term second order system +sys2 = ss([-2.5 0;1 1.5],[1;3],[1 2],2.5) +G2 = tf([2.5, 9.5, 6.125],[1, 1, -3.75]) +H2 = zpk([-2.97703296142690, -0.82296703857310], [1.5, -2.5], 2.5) + +w = 10.0 .^ range(-2,2,50) + +@test ss(sys2) == sys2 +@test freqresp(ss(G2), w) ≈ freqresp(sys2, w) rtol=1e-15 +@test freqresp(ss(H2), w) ≈ freqresp(sys2, w) rtol=1e-15 + +@test tf(sys2) ≈ G2 rtol=1e-15 +@test tf(G2) == G2 +@test tf(H2) ≈ G2 rtol=1e-15 + +@test zpk(sys2) ≈ H2 rtol=1e-15 +@test zpk(G2) ≈ H2 rtol=1e-15 +@test zpk(H2) == H2 + +# Non-proper system +G3 = tf([4.0, 1, 5],[2, 3]) +H3 = zpk([(-1+im*sqrt(79))/8, (-1-im*sqrt(79))/8], [-3/2], 2) + +@test tf(G3) == G3 +@test tf(H3) ≈ G3 rtol=1e-2 + +@test zpk(G3) ≈ H3 rtol=1e-15 +@test zpk(H3) == H3 + +# Test complex 1 +A = [1.0 + im 1; 0 -2-3im] +B = [0;2] +C = [1 0] +D = zeros(1,1); +sys = ss(A,B,C,D) +sys2 = convert(TransferFunction, sys) +w = 10.0 .^ range(-2,2,50) +@test freqresp(sys, w) ≈ freqresp(sys2, w) +csort = v -> sort(v, lt = (x,y) -> abs(x) < abs(y)) +@test csort(pole(zpk(sys2))) ≈ [1+im, -2.0-3im] + +# Test complex 2 +sys4 = ss([-1.0+im 0;1 1],[1;0],[1 im],im) +tf(sys4) +G4 = tf([1.0im,2,-2],[1,-1im,-1+im]) +H4 = zpk([im*(1 + sqrt(1+2im)), im*(1 - sqrt(1+2im))], [-1+im,1], 1.0im) + +w = 10.0 .^ range(-2,2,50) + +@test ss(sys4) == sys4 +@test freqresp(ss(G4), w) ≈ freqresp(sys4, w) rtol=1e-15 +@test freqresp(ss(H4), w) ≈ freqresp(sys4, w) rtol=1e-15 + +@test tf(sys4) ≈ G4 rtol=1e-15 +@test tf(G4) == G4 +@test tf(H4) ≈ G4 rtol=1e-15 + +@test zpk(sys4) ≈ H4 rtol=1e-15 +@test zpk(G4) ≈ H4 rtol=1e-15 +@test zpk(H4) == H4 + +# Test conversion of Int system with direct term +w = [0.5, 1, 1.5, 2.0] +G = tf([2, 3],[5, 1]) +@test freqresp(ss(G), w) ≈ freqresp(G, w) + + +sys2 = ss([-1 0;1 1],[1;0],[1 1],1) +G2 = tf([1,1,-1],[1,0,-1]) +sys3 = ss([0 1;1 0],[0;1],[1 1],0) +G3 = tf([1,1],[1,0,-1]) #Test ss2tf -@test f11 ≈ f12 atol=1e-15 -@test f21 ≈ f22 atol=1e-15 +@test tf(sys2) ≈ G2 atol=1e-15 +@test tf(sys3) ≈ G3 atol=1e-15 -#Test tf2ss +#Test TransferFunction -> StateSpace #ss is not unique so go back and forth -@test f11 ≈ tf(ss(f12)) atol=1e-15 -@test f21 ≈ tf(ss(f22)) atol=1e-15 +@test G2 ≈ tf(ss(G2)) atol=1e-15 +@test G3 ≈ tf(ss(G3)) atol=1e-15 + + +## Test some rather trivial conversions of numeric types + +b = 1.5 +D22 = [1.0 2.0; 3.0 4.0] + +@test convert(StateSpace{Float64,Matrix{Float64}}, D22) == ss(D22) +@test convert(StateSpace{Float64,Matrix{Float64}}, b) == ss(b) + +@test convert(TransferFunction{ControlSystems.SisoRational{Float64}}, D22) == tf(D22) +@test convert(TransferFunction{ControlSystems.SisoRational{Float64}}, b) == tf(b) + +@test convert(TransferFunction{ControlSystems.SisoZpk{Float64,Complex128}}, D22) == zpk(D22) +@test convert(TransferFunction{ControlSystems.SisoZpk{Float64,Complex128}}, b) == zpk(b) + + +# Error, not proper +@test_throws ErrorException ss(tf([1,0],[1])) + end diff --git a/test/test_discrete.jl b/test/test_discrete.jl index ed220f459..31a6337e2 100644 --- a/test/test_discrete.jl +++ b/test/test_discrete.jl @@ -8,37 +8,48 @@ C_222_d = ss([-5 -3; 2 -9], [1 0; 0 2], eye(2), eye(2)) @test c2d(ss(4*eye(2)), 0.5, :zoh) == (ss(4*eye(2), 0.5), zeros(0, 2)) @test c2d(ss(4*eye(2)), 0.5, :foh) == (ss(4*eye(2), 0.5), zeros(0, 2)) @test_c2d(c2d(C_111, 0.01, :zoh), -ss([0.951229424500714], [0.019508230199714396], [3], [0], 0.01), [1 0]) + ss([0.951229424500714], [0.019508230199714396], [3], [0], 0.01), [1 0]) @test_c2d(c2d(C_111, 0.01, :foh), -ss([0.951229424500714], [0.01902855227625244], [3], [0.029506188017136226], 0.01), -[1 -0.009835396005712075]) + ss([0.951229424500714], [0.01902855227625244], [3], [0.029506188017136226], 0.01), + [1 -0.009835396005712075]) @test_c2d(c2d(C_212, 0.01, :zoh), -ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819], -[0.009466805409666932; 0.019219966830212765], eye(2), [0; 0], 0.01), -[1 0 0; 0 1 0]) + ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819], + [0.009466805409666932; 0.019219966830212765], eye(2), [0; 0], 0.01), + [1 0 0; 0 1 0]) @test_c2d(c2d(C_212, 0.01, :foh), -ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482], -[0.008957940478201584; 0.018468989584974498], eye(2), [0.004820885889482196; -0.009738343195298675], 0.01), [1 0 -0.004820885889482196; 0 1 -0.009738343195298675]) + ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482], + [0.008957940478201584; 0.018468989584974498], eye(2), [0.004820885889482196; + 0.009738343195298675], 0.01), [1 0 -0.004820885889482196; 0 1 -0.009738343195298675]) @test_c2d(c2d(C_221, 0.01, :zoh), -ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819], -[0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465], -[1.0 0.0], [0.0 0.0], 0.01), [1 0 0 0; 0 1 0 0]) + ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819], + [0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465], + [1.0 0.0], [0.0 0.0], 0.01), [1 0 0 0; 0 1 0 0]) @test_c2d(c2d(C_221, 0.01, :foh), -ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482], -[0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387], -[1 0], [0.004917457305816479 -9.657141633428213e-5], 0.01), -[1 0 -0.004917457305816479 9.657141633428213e-5; -0 1 -3.219047211142736e-5 -0.009706152723187249]) + ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482], + [0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387], + [1 0], [0.004917457305816479 -9.657141633428213e-5], 0.01), + [1 0 -0.004917457305816479 9.657141633428213e-5; + 0 1 -3.219047211142736e-5 -0.009706152723187249]) @test_c2d(c2d(C_222_d, 0.01, :zoh), -ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819], -[0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465], -eye(2), eye(2), 0.01), [1 0 0 0; 0 1 0 0]) + ss([0.9509478368863918 -0.027970882212682433; 0.018647254808454958 0.9136533272694819], + [0.009753161420545834 -0.0002863560108789034; 9.54520036263011e-5 0.019124514826586465], + eye(2), eye(2), 0.01), [1 0 0 0; 0 1 0 0]) @test_c2d(c2d(C_222_d, 0.01, :foh), -ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482], -[0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387], -eye(2), [1.0049174573058164 -9.657141633428213e-5; 3.219047211142736e-5 1.0097061527231872], 0.01), -[1 0 -0.004917457305816479 9.657141633428213e-5; 0 1 -3.219047211142736e-5 -0.009706152723187249]) + ss([0.9509478368863921 -0.027970882212682433; 0.018647254808454954 0.913653327269482], + [0.009511049106772921 -0.0005531086285713394; 0.00018436954285711309 0.018284620042117387], + eye(2), [1.0049174573058164 -9.657141633428213e-5; 3.219047211142736e-5 1.0097061527231872], 0.01), + [1 0 -0.004917457305816479 9.657141633428213e-5; 0 1 -3.219047211142736e-5 -0.009706152723187249]) + + +# Test c2d for transfer functions +G = tf([1, 1], [1, 3, 1]) +Gd = c2d(G, 0.2) +@test Gd ≈ tf([0, 0.165883310712090, -0.135903621603238], [1.0, -1.518831946985175, 0.548811636094027], 0.2) rtol=1e-14 + +# c2d on a zpk model should arguably return a zpk model +@test_broken typeof(c2d(zpk(f), 1)) <: TransferFunction{<:ControlSystems.SisoZpk} + + # ERRORS @test_throws ErrorException c2d(ss([1], [2], [3], [4], 0.01), 0.01) # Already discrete diff --git a/test/test_freqresp.jl b/test/test_freqresp.jl index 212b94bbc..c8e854a34 100644 --- a/test/test_freqresp.jl +++ b/test/test_freqresp.jl @@ -4,6 +4,8 @@ H = [tf(0) tf([3, 0],[1, 1, 10]) ; tf([1, 1],[1, 5]) tf([2],[1, 6])] G = ss([-5 0 0 0; 0 -1 -2.5 0; 0 4 0 0; 0 0 0 -6], [2 0; 0 1; 0 0; 0 2], [0 3 0 0; -2 0 0 1], [0 0; 1 0]) + + @test evalfr(H, -6) == [0.0 -0.45; 5.0 Inf] @test evalfr(H, -5) == [0.0 -0.5; Inf 2.0] @test evalfr(H, -1) == [0.0 -0.3; 0.0 0.4] @@ -14,26 +16,71 @@ G = ss([-5 0 0 0; 0 -1 -2.5 0; 0 4 0 0; 0 0 0 -6], [2 0; 0 1; 0 0; 0 2], @test evalfr(G, -1) ≈ [0.0 -0.3; 0.0 0.4] @test evalfr(G, 0) ≈ [0.0 0.0; 0.2 1/3] - -## Freqresp +## Constant system w = logspace(-1,1) -H1 = tf(1) -G1 = ss(1) -@test freqresp(G1, w) == ones(Complex128, length(w), 1, 1) -@test freqresp(H1, w) == ones(Complex128, length(w), 1, 1) +sys1 = ss(1) +G1 = tf(1) +H1 = zpk(1) +resp1 = ones(Complex128, length(w), 1, 1) + +@test evalfr(sys1, im*w[1]) == fill(resp1[1], 1, 1) +@test evalfr(G1, im*w[1]) == fill(resp1[1], 1, 1) +@test evalfr(H1, im*w[1]) == fill(resp1[1], 1, 1) + +@test freqresp(sys1, w) == resp1 +@test freqresp(G1, w) == resp1 +@test freqresp(H1, w) == resp1 + +## First order system + +sys2 = ss(-1, 1, 1, 1) +G2 = tf([1, 2], [1,1]) +H2 = zpk([-2], [-1.0], 1.0) +resp2 = reshape((im*w + 2)./(im*w + 1), length(w), 1, 1) + +@test evalfr(sys2, im*w[1]) ≈ fill(resp2[1], 1, 1) +@test evalfr(G2, im*w[1]) == fill(resp2[1], 1, 1) +@test evalfr(H2, im*w[1]) == fill(resp2[1], 1, 1) + +@test freqresp(sys2, w) ≈ resp2 rtol=1e-15 +@test freqresp(G2, w) == resp2 +@test freqresp(H2, w) == resp2 + + +## Complex-coefficient system +sys3 = ss(-1+im, 1, (1-im), (1-im)) +G3 = (1-im)*tf([1,2-im], [1,1-im]) +H3 = zpk([-2+im], [-1+im], (1-im)) +resp3 = reshape((1-im)*(2 + im*(w-1))./(1 + im*(w-1)), length(w), 1, 1) + +@test evalfr(sys3, im*w[1]) ≈ fill(resp3[1], 1, 1) rtol=1e-16 +@test evalfr(G3, im*w[1]) == fill(resp3[1], 1, 1) +@test evalfr(H3, im*w[1]) == fill(resp3[1], 1, 1) + +@test freqresp(sys3, w) ≈ resp3 rtol=1e-15 +@test freqresp(G3, w) ≈ resp3 rtol=1e-16 +@test freqresp(H3, w) == resp3 + -H2 = tf(1, [1,1]) -G2 = ss(-1, 1, 1, 0) -@test freqresp(G2, w) == reshape(1./(1 + im*w), length(w), 1, 1) -@test freqresp(H2, w) == reshape(1./(1 + im*w), length(w), 1, 1) ## Shortcut notation for evalfr ## F = tf([1],[1,0.5],-1) omega = 2 z = 0.5(1+im) -@test F(omega,true)[1] == 1/(exp(-im*2)+0.5) +# This test is not correct if we dont have a sample time +# @test_throws F(omega,true)[1] == 1/(exp(-im*2)+0.5) +@test F(z,false)[1] == 1/(z+0.5) +@test_throws ErrorException F(z,true) + +F = [tf([1],[1,0.5],2.0) 3*tf([1],[1,0.5],2.0)] +omegas = [1,2] +z = 0.5(1+im) +@test F(omegas[1],true) ≈ [1 3].*1/(exp(im*2)+0.5) atol=1e-14 +@test F(omegas[2],true) == [1 3].*1/(exp(2*im*2)+0.5) +@test F(omegas,true) ≈ [k/(exp(omega*im*2)+0.5) for omega=omegas, o=1:1, k=[1,3]] atol=1e-14 +@test F(omegas,false) ≈ [k/(omega+0.5) for omega=omegas, o=1:1, k=[1,3]] atol=1e-14 @test F(z,false)[1] == 1/(z+0.5) @test_throws ErrorException F(z,true) diff --git a/test/test_generalizedtf.jl b/test/test_generalizedtf.jl deleted file mode 100644 index 191a1faf4..000000000 --- a/test/test_generalizedtf.jl +++ /dev/null @@ -1,16 +0,0 @@ -@testset "test_generalizedtf" begin - -# CONTINUOUS -C_011 = tfg("(s+2)") -C_111 = tfg("(s+1)/(s+2)") -@test C_011*C_111 == tfg("(s+2)*((s+1)/(s+2))") - -#We might want to make evalfr scalar -@test C_011(1im) == reshape([2+1im;],1,1) -@test (C_111*C_011)(im) == reshape([1.0+1im],1,1) - -@test tf(C_111) == tf([1,1],[1,2]) -@test zpk(C_011*C_111) == zpk([-1,-2],[-2],1) - -@test bode(C_111*C_011, logspace(-1,1)) == bode(tfg("(s+2)*((s+1)/(s+2))"), logspace(-1,1)) -end diff --git a/test/test_linalg.jl b/test/test_linalg.jl index 94e98d0bc..7679dde85 100644 --- a/test/test_linalg.jl +++ b/test/test_linalg.jl @@ -35,7 +35,7 @@ s = tf("s") f_C_211 = (s+2)*(s+3)/((s+4)*(s+5)) # biquad passband omega0 = 52.0; Q = 10 -f_C_211_bis = (s/(Q*omega0)) / ((s/omega0)^2 + s/(Q*omega0) + 1 ); +f_C_211_bis = (s/(Q*omega0)) / ((s/omega0)^2 + s/(Q*omega0) + 1) C_22tf = [0 tf([3,0],[1,1,10]);tf([1,1],[1,5]) tf(2,[1,6])] @@ -86,36 +86,38 @@ D_static = ss([0.704473 1.56483; -1.6661 -2.16041], 0.07) tolHinf = 1e-12 @test norm(C_212, Inf, tol=tolHinf) ≈ 0.242535625036333 atol=tolHinf @test norm(C_222, Inf, tol=tolHinf) ≈ 0.242535625036333 atol=tolHinf -ninf, fpeak = norminf(C_732, tol=tolHinf) -@test ninf ≈ 4.899135403568278 atol=(10*tolHinf) -@test fpeak ≈ 6.112977387441163 atol=1e-6 +Ninf, ω_peak = norminf(C_732, tol=tolHinf) +@test Ninf ≈ 4.899135403568278 atol=(10*tolHinf) +@test ω_peak ≈ 6.112977387441163 atol=1e-6 @test norm(f_C_211, Inf, tol=tolHinf) ≈ 1.0 atol=(2*tolHinf) -@test norminf(f_C_211_bis, tol=tolHinf)[2] ≈ 52.0 +Ninf, ω_peak = norminf(f_C_211_bis, tol=tolHinf) +@test Ninf ≈ 1.0 +@test ω_peak ≈ 52.0 @test norm(1/(s-1), Inf, tol=tolHinf) ≈ 1.0 # unstable system -ninf, fpeak = norminf(C_22tf, tol=tolHinf) -@test ninf ≈ 3.014974550173459 atol=(10*tolHinf) -@test fpeak ≈ 3.162123338668049 atol=1e-8 - -ninf, fpeak = norminf(D_221, tol=tolHinf) -@test ninf ≈ 17.794697451669421 atol=(20*tolHinf) -@test fpeak ≈ 0 atol=1e-8 -ninf, fpeak = norminf(D_422, tol=tolHinf) -@test ninf ≈ 3.360351099392252 atol=(10*tolHinf) -@test fpeak ≈ 8.320643111730551 atol=1e-8 -ninf, fpeak = norminf(D_311, tol=tolHinf) -@test ninf ≈ 4.458729529942810 atol=(10*tolHinf) -@test fpeak ≈ 11.878021287349698 atol=1e-6 -ninf, fpeak = norminf(C_static, tol=tolHinf) -@test ninf ≈ 1.760164138376307 atol=1e-12 -@test fpeak ≈ 0 atol=1e-12 -ninf, fpeak = norminf(D_static, tol=tolHinf) -@test ninf ≈ 3.205246234285972 atol=1e-12 -@test fpeak ≈ 0 atol=1e-12 +Ninf, ω_peak = norminf(C_22tf, tol=tolHinf) +@test Ninf ≈ 3.014974550173459 atol=(10*tolHinf) +@test ω_peak ≈ 3.162123338668049 atol=1e-8 + +Ninf, ω_peak = norminf(D_221, tol=tolHinf) +@test Ninf ≈ 17.794697451669421 atol=(20*tolHinf) +@test ω_peak ≈ 0 atol=1e-8 +Ninf, ω_peak = norminf(D_422, tol=tolHinf) +@test Ninf ≈ 3.360351099392252 atol=(10*tolHinf) +@test ω_peak ≈ 8.320643111730551 atol=1e-8 +Ninf, ω_peak = norminf(D_311, tol=tolHinf) +@test Ninf ≈ 4.458729529942810 atol=(10*tolHinf) +@test ω_peak ≈ 11.878021287349698 atol=1e-6 +Ninf, ω_peak = norminf(C_static, tol=tolHinf) +@test Ninf ≈ 1.760164138376307 atol=1e-12 +@test ω_peak ≈ 0 atol=1e-12 +Ninf, ω_peak = norminf(D_static, tol=tolHinf) +@test Ninf ≈ 3.205246234285972 atol=1e-12 +@test ω_peak ≈ 0 atol=1e-12 A = [1 100 10000; .01 1 100; .0001 .01 1] -T, P, B = ControlSystems.balance(A) +T, P, B = balance(A) # The scaling is BLAS dependent. However, the ratio should be the same on all # machines. We just need to check that T == res * constant res_diag = [512, 8, 0.0625] diff --git a/test/test_matrix_comps.jl b/test/test_matrix_comps.jl index aa4d4094a..20c9ecf2c 100644 --- a/test/test_matrix_comps.jl +++ b/test/test_matrix_comps.jl @@ -13,6 +13,10 @@ sysr, G = balreal(sys) sysb,T = ControlSystems.balance_statespace(sys) Ab,Bb,Cb,T = ControlSystems.balance_statespace(A,B,C) +@test Ab*T ≈ T*A +@test Bb ≈ T*B +@test Cb*T ≈ C + @test sysb.A ≈ Ab @test ControlSystems.balance_transform(A,B,C) ≈ ControlSystems.balance_transform(sys) @@ -32,4 +36,26 @@ D2 = eye(2) # Discrete system can have direct term @test covar(ss(A,B,C,D2,0.1),W) ≈ [1.00011010837831 -1.0098377309782909e-5; -1.0098377309782909e-5 1.00011010837831] + +# TODO test in Julia 0.7 to see if supported +# # Test special matrices +# As = sparse(A) +# Bs = sparse(B) +# Cs = sparse(C) +# Asb,Bsb,Csb,Ts = ControlSystems.balance_statespace(As,Bs,Cs) #Error no LAPACK function +# +# @test Abs*Ts ≈ Ts*As +# @test Bbs ≈ Ts*Bs +# @test Cbs*Ts ≈ Cs + +# Test special values +Ar = rationalize.(A) +Br = rationalize.(B) +Cr = rationalize.(C) +Arb,Brb,Crb,Tr = ControlSystems.balance_statespace(Ar,Br,Cr) + +@test Arb*Tr ≈ Tr*Ar +@test Brb ≈ Tr*Br +@test Crb*Tr ≈ Cr + end diff --git a/test/test_promotion.jl b/test/test_promotion.jl new file mode 100644 index 000000000..68b7b637d --- /dev/null +++ b/test/test_promotion.jl @@ -0,0 +1,47 @@ +@testset "test_promotion" begin +G = tf(1) +sys = ss(1) +TG = typeof(G) +Tsys = typeof(sys) +@test promote_type(Tsys,TG) == Tsys + +G = tf(1.) +sys = ss(1.) +TG = typeof(G) +Tsys = typeof(sys) +@test promote_type(Tsys,TG) == Tsys + + +G = tf(1.) +sys = ss(1) +TG = typeof(G) +Tsys = typeof(sys) +@test promote_type(Tsys,TG) == typeof(ss(1.)) + + +G = zpk(1) +sys = ss(1) +TG = typeof(G) +Tsys = typeof(sys) +@test promote_type(Tsys,TG) == typeof(ss(1)) + +G1 = zpk(1) +G2 = zpk(1.) +TG1 = typeof(G1) +TG2 = typeof(G2) +@test promote_type(TG1,TG2) == TG2 + +@test promote_type(typeof(zpk(1)), typeof(zpk(1))) == typeof(zpk(1)) +@test promote_type(typeof(tf(1)), typeof(tf(1))) == typeof(tf(1)) +@test promote_type(typeof(ss(1)), typeof(ss(1))) == typeof(ss(1)) + +@test promote(ss(1), ss(1.)) == (ss(1.), ss(1.)) +@test promote(ss(1), tf(1.)) == (ss(1.), ss(1.)) +@test promote(ss(1), zpk(1.)) == (ss(1.), ss(1.)) +@test promote(zpk(1), tf(1.)) == (zpk(1.), zpk(1.)) +@test promote(tf([1 + im]), tf([1.])) == (tf([1. + im]), tf([1. + 0*im])) +@test promote(zpk(1), tf(1)) == (zpk(1.), zpk(1.)) +@test promote(ss(1), ss([1.*im])) == (ss([1+0*im]),ss([1.*im])) +@test promote(ss(1), tf([1.*im])) == (ss([1+0*im]),ss([1.*im])) + +end diff --git a/test/test_simplification.jl b/test/test_simplification.jl index ceee144bc..1d48073ac 100644 --- a/test/test_simplification.jl +++ b/test/test_simplification.jl @@ -21,6 +21,6 @@ sysmin = minreal(sys) t = 0:0.1:10 y1,x1 = step(sys,t)[[1,3]] y2,x2 = step(sysmin,t)[[1,3]] -@test sum(abs2,y1.-y2) < √(eps()) # Test that the output from the two systems are the same +@test sum(abs2,y1.-y2) < 1e-6 # Test that the output from the two systems are the same end diff --git a/test/test_statespace.jl b/test/test_statespace.jl index 084d802b8..067690cb3 100644 --- a/test/test_statespace.jl +++ b/test/test_statespace.jl @@ -37,8 +37,7 @@ D_222_d = ss(da_2, [1 0; 0 2], eye(2), eye(2), 0.005) D_022 = ss(4*eye(2), 0.005) # Definition of input, output and state names -C_222_d_n = ss(a_2, [1 0; 0 2], eye(2), eye(2), - statenames=["i","u"],inputnames=String("e"),outputnames="theta") +C_222_d_n = ss(a_2, [1 0; 0 2], eye(2), eye(2)) # TESTS # Contstuct with scalars @@ -86,30 +85,27 @@ C_222_d_n = ss(a_2, [1 0; 0 2], eye(2), eye(2), @test C_222[1,1:2] == C_221 @test size(C_222[1,[]]) == (1,0) -# Naming signals -@test C_222_d_n.statenames == String["i","u"] -@test C_222_d_n.inputnames == String["e1","e2"] -@test C_222_d_n.outputnames == String["theta1","theta2"] + +A = [-1.0 -2.0; 0.0 -1.0] +B = [0.0; -2.0] +C = [1.0 1.0] +D = 1.0 +sys = ss(A, B, C, D) + +@test sys + 1.0 == ss(A, B, C, D + 1.0) +@test 2.0 + sys == ss(A, B, C, D + 2.0) + +@test -sys == ss(A, B, -C, -D) + # Printing -res = ("StateSpace:\nA = \n x1 x2 \n x1 -5.0 -3.0 \n x2"* - " 2.0 -9.0 \nB = \n u1 u2 \n x1 1.0 0.0 \n"* - " x2 0.0 2.0 \nC = \n x1 x2 \n y1 1.0 0.0 \n"* - " y2 0.0 1.0 \nD = \n u1 u2 \n y1 0.0 0.0 \n"* - " y2 0.0 0.0 \n\nContinuous-time state-space model") +res = ("ControlSystems.StateSpace{Float64,Array{Float64,2}}\nA = \n -5.0 -3.0\n 2.0 -9.0\nB = \n 1.0 0.0\n 0.0 2.0\nC = \n 1.0 0.0\n 0.0 1.0\nD = \n 0.0 0.0\n 0.0 0.0\n\nContinuous-time state-space model") @test sprint(show, C_222) == res -res = ("StateSpace:\nA = \n x1 x2 \n x1 0.2 -0.8 \n x2"* - " -0.8 0.07 \nB = \n u1 u2 \n x1 1.0 0.0 \n"* - " x2 0.0 2.0 \nC = \n x1 x2 \n y1 1.0 0.0 \n"* - " y2 0.0 1.0 \nD = \n u1 u2 \n y1 0.0 0.0 \n"* - " y2 0.0 0.0 \n\nSample Time: 0.005 (seconds)\n"* - "Discrete-time state-space model") +res = ("ControlSystems.StateSpace{Float64,Array{Float64,2}}\nA = \n 0.2 -0.8 \n -0.8 0.07\nB = \n 1.0 0.0\n 0.0 2.0\nC = \n 1.0 0.0\n 0.0 1.0\nD = \n 0.0 0.0\n 0.0 0.0\n\nSample Time: 0.005 (seconds)\nDiscrete-time state-space model") @test sprint(show, D_222) == res -res = ("StateSpace:\nD = \n u1 u2 \n y1 4.0 0.0 \n y2 "* - " 0.0 4.0 \n\nStatic gain") +res = ("ControlSystems.StateSpace{Float64,Array{Float64,2}}\nD = \n 4.0 0.0\n 0.0 4.0\n\nStatic gain") @test sprint(show, C_022) == res -res = ("StateSpace:\nD = \n u1 u2 \n y1 4.0 0.0 \n y2 "* - " 0.0 4.0 \n\nSample Time: 0.005 (seconds)\nStatic gain") +res = "ControlSystems.StateSpace{Float64,Array{Float64,2}}\nD = \n 4.0 0.0\n 0.0 4.0\n\nSample Time: 0.005 (seconds)\nStatic gain" @test sprint(show, D_022) == res # Errors diff --git a/test/test_synthesis.jl b/test/test_synthesis.jl index 4b85d7e95..52f0091b2 100644 --- a/test/test_synthesis.jl +++ b/test/test_synthesis.jl @@ -1,6 +1,6 @@ @testset "test_synthesis" begin -P = tf(1,[1,1]) -C = tf([1,1],[1,0]) +P = tf(1.,[1.,1]) +C = tf([1.,1],[1.,0]) L = P*C Lsys = ss(L) @@ -11,12 +11,21 @@ R = [1,1] S = [1] T = [1] -@test isapprox(minreal(feedback(P,C),1e-5), minreal(feedback(L),1e-5)) +@test isapprox(minreal(feedback(P,C),1e-5), tf([1,0],[1,2,1]), rtol = 1e-5) @test isapprox(numpoly(minreal(feedback(L),1e-5))[1].a, numpoly(tf(1,[1,1]))[1].a)# This test is ugly, but numerical stability is poor for minreal @test feedback2dof(B,A,R,S,T) == tf(B.*T, conv(A,R) + [0;0;conv(B,S)]) @test feedback2dof(P,R,S,T) == tf(B.*T, conv(A,R) + [0;0;conv(B,S)]) @test isapprox(pole(minreal(tf(feedback(Lsys)),1e-5)) , pole(minreal(feedback(L),1e-5)), atol=1e-5) +Pint = tf(1,[1,1]) +Cint = tf([1,1],[1,0]) +Lint = P*C + +@test_broken isapprox(minreal(feedback(Pint,Cint),1e-5), tf([1,0],[1,2,1]), rtol = 1e-5) +@test isapprox(numpoly(minreal(feedback(Lint),1e-5))[1].a, numpoly(tf(1,[1,1]))[1].a)# This test is ugly, but numerical stability is poor for minreal +@test isapprox(pole(minreal(tf(feedback(Lsys)),1e-5)) , pole(minreal(feedback(L),1e-5)), atol=1e-5) + + @test_throws ErrorException feedback(ss(1),ss(1)) @test_throws ErrorException feedback(ss(eye(2), ones(2,2), ones(1,2),0)) end diff --git a/test/test_timeresp.jl b/test/test_timeresp.jl index bb755666a..f38a3df9a 100644 --- a/test/test_timeresp.jl +++ b/test/test_timeresp.jl @@ -8,11 +8,16 @@ Q = eye(2) R = eye(1) L = lqr(sys,Q,R) -u(i,x) = -L*x # Form control law +u(x,i) = -L*x # Form control law t=0:0.1:50 -x0 = [1,0] -y, t, x, uout = lsim(sys,u,t,x0=x0) -@test sum(abs.(x[end,:])) < eps() +x0 = [1.,0] +y, t, x, uout = lsim(sys,u,t,x0=x0) # Continuous time + +th = 1e-6 +@test sum(abs.(x[end,:])) < th + +y, t, x, uout = lsim(c2d(sys,0.1)[1],u,t,x0=x0) # Discrete time +@test sum(abs.(x[end,:])) < th #Do a manual simulation with uout ym, tm, xm = lsim(sys, uout, t, x0=x0) @@ -31,7 +36,7 @@ yd, td, xd = lsim(sysdfb, zeros(501), t, x0=x0) @test x ≈ xd -#Test step and impulse +#Test step and impulse Continuous t0 = 0:0.05:2 systf = [tf(1,[1,1]) 0; 0 tf([1,-1],[1,2,1])] sysss = ss([-1 0 0; 0 -2 -1; 0 1 0], [1 0; 0 1; 0 0], [1 0 0; 0 1 -1], 0) @@ -42,9 +47,42 @@ xreal = zeros(length(t0), 3, 2) y, t, x = step(systf, t0) yreal[:,1,1] = 1-exp.(-t) yreal[:,2,2] = -1+exp.(-t)+2*exp.(-t).*t -@test y ≈ yreal atol=1e-14 +@test y ≈ yreal atol=1e-4 #Step ss y, t, x = step(sysss, t) +@test y ≈ yreal atol=1e-4 +xreal[:,1,1] = yreal[:,1,1] +xreal[:,2,2] = exp.(-t).*t +xreal[:,3,2] = exp.(-t).*(-t-1) + 1 +@test x ≈ xreal atol=1e-5 + +#Impulse tf +y, t, x = impulse(systf, t) +yreal[:,1,1] = exp.(-t) +yreal[:,2,2] = exp.(-t).*(1 - 2.*t) +@test y ≈ yreal atol=1e-2 +#Impulse ss +y, t, x = impulse(sysss, t) +@test y ≈ yreal atol=1e-2 +xreal[:,1,1] = yreal[:,1,1] +xreal[:,2,2] = -exp.(-t).*t + exp.(-t) +xreal[:,3,2] = exp.(-t).*t +@test x ≈ xreal atol=1e-1 + +#Test step and impulse Discrete +t0 = 0:0.05:2 +systf = [tf(1,[1,1]) 0; 0 tf([1,-1],[1,2,1])] +sysss = ss([-1 0 0; 0 -2 -1; 0 1 0], [1 0; 0 1; 0 0], [1 0 0; 0 1 -1], 0) +yreal = zeros(length(t0), 2, 2) +xreal = zeros(length(t0), 3, 2) + +#Step tf +y, t, x = step(systf, t0, method=:zoh) +yreal[:,1,1] = 1-exp.(-t) +yreal[:,2,2] = -1+exp.(-t)+2*exp.(-t).*t +@test y ≈ yreal atol=1e-14 +#Step ss +y, t, x = step(sysss, t, method=:zoh) @test y ≈ yreal atol=1e-14 xreal[:,1,1] = yreal[:,1,1] xreal[:,2,2] = exp.(-t).*t @@ -52,12 +90,12 @@ xreal[:,3,2] = exp.(-t).*(-t-1) + 1 @test x ≈ xreal atol=1e-14 #Impulse tf -y, t, x = impulse(systf, t) +y, t, x = impulse(systf, t, method=:zoh) yreal[:,1,1] = exp.(-t) yreal[:,2,2] = exp.(-t).*(1 - 2.*t) @test y ≈ yreal atol=1e-14 #Impulse ss -y, t, x = impulse(sysss, t) +y, t, x = impulse(sysss, t, method=:zoh) @test y ≈ yreal atol=1e-14 xreal[:,1,1] = yreal[:,1,1] xreal[:,2,2] = -exp.(-t).*t + exp.(-t) @@ -68,17 +106,32 @@ xreal[:,3,2] = exp.(-t).*t #Step response of discrete system with specified final time G = tf([1], [1; zeros(3)], 1) y, t2, x = step(G, 10) -@test y ≈ [zeros(3); ones(8)] atol=1e-14 -@test t2 ≈ 0:1:10 atol=1e-14 +@test y ≈ [zeros(3); ones(8)] atol=1e-5 +@test t2 ≈ 0:1:10 atol=1e-5 #Impulse response of discrete system to final time that is not mulitple of the sample time G = tf([1], [1; zeros(3)], 0.3) y, t2, x = step(G, 2) -@test y ≈ [zeros(3); ones(4)] atol=1e-14 -@test t2 ≈ 0:0.3:1.8 atol=1e-14 - - +@test y ≈ [zeros(3); ones(4)] atol=1e-5 +@test t2 ≈ 0:0.3:1.8 atol=1e-5 #Make sure t was never changed @test t0 == t + + +@testset "Simulators" begin + + +h = 0.1 +Tf = 20 +t = 0:h:Tf +P = ss(tf(1,[2,1])^2) +reference(x,t) = [1.] +s = Simulator(P, reference) +x0 = [0.,0] +tspan  = (0.0,Tf) +sol = solve(s, x0, tspan, OrdinaryDiffEq.Tsit5()) +@test step(P,t)[3] ≈ reduce(hcat,sol.(t))' +end + end diff --git a/test/test_transferfunction.jl b/test/test_transferfunction.jl index daa78e6de..a3071194b 100644 --- a/test/test_transferfunction.jl +++ b/test/test_transferfunction.jl @@ -32,10 +32,25 @@ z = tf("z", 0.005) @test C_022 == tf([4 0;0 4]) @test D_022 == tf([4 0;0 4], 0.005) +# Test equality +@test tf([1,2], [2,3,4]) == tf(2*[1,2], 2*[2,3,4]) +@test tf([1.0], [2.0,3.0]) == tf(π*[1.0], π*[2.0,3.0]) +@test tf([1.0+2.0im], [2.0+im,3.0]) == tf((π+im)*[1+2.0im], (π+im)*[2.0+im,3.0]) + +# Test approximate equlity +# rtol should just be on the order of ϵ, no particular reason that exactly ϵ +# would work, but apparently it does +ϵ = 1e-14 +@test tf([1,2], [2,3,4]) != tf(2*[1,2], (2+ϵ)*[2,3,4]) +@test tf([1,2], [2,3,4]) ≈ tf(2*[1,2], (2+ϵ)*[2,3,4]) rtol=ϵ +@test tf([1.0], [2.0,3.0]) != tf((π+ϵ)*[1.0], π*[2.0,3.0]) +@test tf([1.0], [2.0,3.0]) ≈ tf((π+ϵ)*[1.0], π*[2.0,3.0]) rtol=ϵ +@test tf([1.0+2.0im], [2.0+im,3.0]) != tf((π+(1+ϵ)im)*[1+2.0im], (π+ϵ+im)*[2.0+(1+ϵ*im),3.0]) +@test tf([1.0+2.0im], [2.0+im,3.0]) ≈ tf((π+(1+ϵ)im)*[1+2.0im], (π+ϵ+im)*[2.0+(1+ϵ)*im,3.0]) rtol=ϵ + # Addition @test C_111 + C_111 == tf([2,14,20], [1,10,25]) -@test C_222 + C_222 == -tf(vecarray(2, 2, [2,20,68,108,90], [0,2,20,62,60], +@test C_222 + C_222 == tf(vecarray(2, 2, [2,20,68,108,90], [0,2,20,62,60], [2,20,68,108,90], [0,2,20,62,60]), vecarray(2, 2, [1,16,94,240,225], [1,16,94,240,225], [1,16,94,240,225], [1,16,94,240,225])) @@ -45,10 +60,9 @@ tf(vecarray(2, 2, [2,20,68,108,90], [0,2,20,62,60], # Subtraction @test C_111 - C_211 == tf([0,3,18,15], [1,13,55,75]) -@test 1 - C_222 == tf(vecarray(2, 2, [0,6,12], [1,7,13], [0,6,12], [1,7,13]), - vecarray(2, 2, [1,8,15], [1,8,15], [1,8,15], [1,8,15])) -@test D_111 - D_211 - tf([0,0.3,-2.55,1.2], [1,-0.7,-0.05,0.075], 0.005) == - tf([0], [1], 0.005) +@test 1 - C_222 == tf(vecarray(2, 2, [0,6,12], [1,7,13], [0,6,12], [1,7,13]), vecarray(2, 2, [1,8,15], [1,8,15], [1,8,15], [1,8,15])) +# We are not doing enough to identify zero numerator here +@test_broken D_111 - D_211 - tf([0,0.3,-2.55,1.2], [1,-0.7,-0.05,0.075], 0.005) == tf([0.0], [1], 0.005) # Multiplication @test C_111 * C_221 == tf(vecarray(1, 2, [1,4,7,6], [0,1,4,4]), @@ -57,7 +71,8 @@ tf(vecarray(2, 2, [2,20,68,108,90], [0,2,20,62,60], vecarray(2, 1, [1,13,55,75], [1,13,55,75])) @test 4*C_222 == tf(vecarray(2, 2, [4,8,12], [0,4,8], [4,8,12], [0,4,8]), vecarray(2, 2, [1,8,15], [1,8,15], [1,8,15], [1,8,15])) -@test D_111 * D_221 - tf(vecarray(1, 2, [1,4,7,6], [0,1,4,4]), + # We are not doing enough to identify zero numerator here +@test_broken D_111 * D_221 - tf(vecarray(1, 2, [1,4,7,6], [0,1,4,4]), vecarray(1, 2, [1,-0.7,-0.05,0.075], [1.0,-0.7,-0.05,0.075]), 0.005) == tf(vecarray(1, 2, [0], [0]), vecarray(1, 2, [1], [1]), 0.005) @@ -75,21 +90,20 @@ tf(vecarray(1, 2, [0], [0]), vecarray(1, 2, [1], [1]), 0.005) @test C_222[1,1:2] == C_221 @test size(C_222[1,[]]) == (1,0) +# Errors +@test_throws ErrorException tf(vecarray(1, 1, [1,7,13,15]), + vecarray(2, 1, [1,10,31,30], [1,10,31,30])) + + # Printing -res = ("TransferFunction:\nInput 1 to Output 1\ns^2 + 2.0s + 3.0\n-----------"* - "------\ns^2 + 8.0s + 15.0\n\nInput 1 to Output 2\ns^2 + 2.0s + 3.0\n-"* - "----------------\ns^2 + 8.0s + 15.0\n\nInput 2 to Output 1\n s + "* - "2.0\n-----------------\ns^2 + 8.0s + 15.0\n\nInput 2 to Output 2\n "* - " s + 2.0\n-----------------\ns^2 + 8.0s + 15.0\n\nContinuous-time"* - " transfer function model") -@test sprint(show, C_222) == res -res = ("TransferFunction:\nInput 1 to Output 1\nz^2 + 2.0z + 3.0\n-----------"* - "------\nz^2 - 0.2z - 0.15\n\nInput 1 to Output 2\nz^2 + 2.0z + 3.0\n-"* - "----------------\nz^2 - 0.2z - 0.15\n\nInput 2 to Output 1\n z + "* - "2.0\n-----------------\nz^2 - 0.2z - 0.15\n\nInput 2 to Output 2\n "* - " z + 2.0\n-----------------\nz^2 - 0.2z - 0.15\n\nSample Time: 0.005 "* - "(seconds)\nDiscrete-time transfer function model") -@test sprint(show, D_222) == res +res = ("TransferFunction:\nInput 1 to Output 1\ns^2 + 2s + 3\n-------------\ns^2 + 8s + 15\n\nInput 1 to Output 2\ns^2 + 2s + 3\n-------------\ns^2 + 8s + 15\n\nInput 2 to Output 1\n s + 2\n-------------\ns^2 + 8s + 15\n\nInput 2 to Output 2\n s + 2\n-------------\ns^2 + 8s + 15\n\nContinuous-time transfer function model") +@test_broken sprint(show, C_222) == res +res = ("TransferFunction:\nInput 1 to Output 1\nz^2 + 2.0z + 3.0\n-----------------\nz^2 - 0.2z - 0.15\n\nInput 1 to Output 2\nz^2 + 2.0z + 3.0\n-----------------\nz^2 - 0.2z - 0.15\n\nInput 2 to Output 1\n z + 2.0\n-----------------\nz^2 - 0.2z - 0.15\n\nInput 2 to Output 2\n z + 2.0\n-----------------\nz^2 - 0.2z - 0.15\n\nSample Time: 0.005 (seconds)\nDiscrete-time transfer function model") +@test_broken sprint(show, D_222) == res + + +@test tf(zpk([1.0 2; 3 4])) == tf([1 2; 3 4]) + # Type stability Continuous-time @test eltype(fill(tf("s"),2)) <: TransferFunction diff --git a/test/test_zpk.jl b/test/test_zpk.jl index 8c7400441..007d5b7da 100644 --- a/test/test_zpk.jl +++ b/test/test_zpk.jl @@ -24,30 +24,42 @@ z = zpk("z", 0.005) # TESTS # Constructors -@test s == zpk([0], [], 1) -@test z == zpk([0], [], 1, 0.005) -@test C_022 == zpk(vecarray(Int64, 2, 2, [], [], [], []), vecarray(Int64, 2, 2, [], [], [], []), [4 0; 0 4]) -@test D_022 == zpk(vecarray(2, 2, [], [], [], []), vecarray(2, 2, [], [], [], []), [4 0; 0 4], 0.005) +@test s == zpk([0], Int64[], 1) +@test z == zpk([0], Int64[], 1, 0.005) +@test C_022 == zpk(vecarray(Int64, 2, 2, Int64[], Int64[], Int64[], Int64[]), vecarray(Int64, 2, 2, Int64[], Int64[], Int64[], Int64[]), [4 0; 0 4]) +@test D_022 == zpk(vecarray(2, 2, Int64[], Int64[], Int64[], Int64[]), vecarray(2, 2, Int64[], Int64[], Int64[], Int64[]), [4 0; 0 4], 0.005) @test C_022 == [zpk(4) 0;0 4] #TODO We might want to fix this #@test D_022 == [zpk(4, 0.005) 0;0 4]) +# Test constructors with empty matrices of type Any +@test zpk([-1.0], [], 1.0) == zpk([-1.0], Float64[], 1.0) +@test zpk([], [-1.0], 1.0) == zpk(Float64[], [-1.0], 1.0) +@test zpk([], [], 1.0) == zpk(Float64[], Float64[], 1.0) +@test zpk([], [1.0+im,1.0-im], 1.0) == zpk(Complex128[], [1.0+im,1.0-im], 1.0) + + + #TODO improve polynomial accuracy se these are equal # Addition @test C_111 + C_111 ≈ zpk([-2,-5], [-5,-5], 2) @test C_222 + C_222 ≈ -zpk(vecarray(2, 2, [-1+sqrt(2)im,-1-sqrt(2)im], [-2], [-1+sqrt(2)im,-1-sqrt(2)im], [-2]), vecarray(2, 2, [-5,-3], [-5,-3], [-5,-3], [-5,-3]), [2 2;2 2]) + zpk(vecarray(2, 2, [-1+sqrt(2)im,-1-sqrt(2)im], [-2], [-1+sqrt(2)im,-1-sqrt(2)im], [-2]), vecarray(2, 2, [-5,-3], [-5,-3], [-5,-3], [-5,-3]), [2 2;2 2]) z1 = [-2.5+sqrt(9-2.5^2)im, -2.5-sqrt(9-2.5^2)im] z2 = [-4.5-sqrt(4.5^2-17), -4.5+sqrt(4.5^2-17)] + + @test C_222 + 1 ≈ zpk(vecarray(2, 2, z1, z2, z1, z2), vecarray(2, 2, [-3, -5], [-3, -5], [-3, -5], [-3, -5]), [2 1; 2 1]) @test D_111 + D_111 ≈ zpk([-2],[0.5],2,0.005) # Subtraction -@test C_111 - C_211 ≈ zpk(tf([0,3,18,15], [1,13,55,75])) +@test C_111 - C_211 ≈ zpk([-1], [-5,-3.0], 3.0) +@test minreal(zpk(tf([0,3.0,18,15], [1,13,55,75])), 1e-6) ≈ zpk([-1.0], [-5,-3.0], 3.0) + @test 1 - C_222 ≈ zpk(tf(vecarray(2, 2, [0,6,12], [1,7,13], [0,6,12], [1,7,13]), vecarray(2, 2, [1,8,15], [1,8,15], [1,8,15], [1,8,15]))) @test D_111 - D_211 - zpk(tf([0,0.3,-2.55,1.2], [1,-0.7,-0.05,0.075], 0.005)) ≈ - zpk(tf([0], [1], 0.005)) + zpk(tf([0], [1], 0.005)) # Multiplication @test C_111 * C_221 ≈ zpk(tf(vecarray(1, 2, [1,4,7,6], [0,1,4,4]), @@ -58,7 +70,7 @@ z2 = [-4.5-sqrt(4.5^2-17), -4.5+sqrt(4.5^2-17)] vecarray(2, 2, [1,8,15], [1,8,15], [1,8,15], [1,8,15]))) @test D_111 * D_221 - zpk(tf(vecarray(1, 2, [1,4,7,6], [0,1,4,4]), vecarray(1, 2, [1,-0.7,-0.05,0.075], [1.0,-0.7,-0.05,0.075]), 0.005)) ≈ -zpk(tf(vecarray(1, 2, [0], [0]), vecarray(1, 2, [1], [1]), 0.005)) + zpk(tf(vecarray(1, 2, [0], [0]), vecarray(1, 2, [1], [1]), 0.005)) # Division @test 1/C_111 ≈ zpk(tf([1,5], [1,2])) @@ -69,19 +81,51 @@ zpk(tf(vecarray(1, 2, [0], [0]), vecarray(1, 2, [1], [1]), 0.005)) # Indexing @test size(C_222) == (2, 2) @test size(C_212) == (2, 1) -@test C_222[1,1] == zpk(tf([1, 2, 3], [1, 8, 15])) -@test C_222[1:1,1] == zpk(tf([1, 2, 3], [1, 8, 15])) +@test C_222[1,1] ≈ zpk(tf([1, 2, 3], [1, 8, 15])) +@test C_222[1:1,1] ≈ zpk(tf([1, 2, 3], [1, 8, 15])) @test C_222[1,1:2] == C_221 -@test size(C_222[1,[]]) == (1,0) +@test size(C_222[1,Int64[]]) == (1,0) # Test that number of poles matter -@test !(zpk([],[1,1],1) == zpk([],[1],1)) +@test !(zpk(Int64[],[1,1],1) == zpk(Int64[],[1],1)) + +# Test SispZpk and SisoTf operations +C_212_tf = tf(C_212) +C_111_tf = tf(C_111) +# Add +@test C_212_tf + C_212 ≈ 2*C_212 +@test C_212 + C_212_tf ≈ 2*C_212 +# Approx +@test C_212_tf + C_212 ≈ 2*C_212_tf +# Minus +@test 2*C_212_tf - C_212 ≈ C_212 +# Multiply +@test C_212_tf*C_111 ≈ C_212*C_111 +@test C_212_tf*C_111 ≈ C_212_tf*C_111_tf + # TODO test printing when it is implemented better +# Tests of minreal +@test minreal(zpk([-5.0], [-5.0, -5.0], 1.0)) == zpk(Float64[], [-5.0], 1.0) +@test minreal(zpk([-1.0, -2.0], [-2.0, -2.0], 1.0)) == zpk([-1.0], [-2.0], 1.0) + +# Tests of minreal (roudning and cancelling complex roots with real roots) +@test minreal(zpk([-2.0+1e-10im,-2.0-1e-10im], [-2.0], 1.0)) == zpk([-2.0], Float64[], 1.0) +@test minreal(zpk([-2.0], [-2.0+1e-10im,-2.0-1e-10im], 1.0)) == zpk(Float64[], [-2.0], 1.0) +@test minreal(zpk([-1.0, -2.0, -2.0+1e-10im,-2.0-1e-10im], [-2.0, -2.0], 1.0)) == zpk([-1.0, -2.0], Float64[], 1.0) +@test minreal(zpk([-1.0, -2.0+1e-10im,-2.0-1e-10im,-2.0], [-2.0, -2.0], 1.0)) == zpk([-1.0, -2.0], Float64[], 1.0) +@test minreal(zpk([-2.0, -2.0], [-1.0, -2.0, -2.0+1e-10im,-2.0-1e-10im], 1.0)) == zpk(Float64[], [-1.0, -2.0], 1.0) +@test minreal(zpk([-2.0, -2.0], [-1.0, -2.0+1e-10im,-2.0-1e-10im,-2.0], 1.0)) == zpk(Float64[], [-1.0, -2.0], 1.0) + +# Test of minreal (systems with no poles / no zeros) +@test minreal(zpk([-1.0, -2.0], Float64[], 2.5)) == zpk([-1.0, -2.0], Float64[], 2.5) +@test minreal(zpk(Float64[], [-1.0, -2.0], 2.5)) == zpk(Float64[], [-1.0, -2.0], 2.5) + # Type stability Continuous and discrete time +# QUESTION: why not use typeof instead of eltype(fill()) @test eltype(fill(zpk("s"),2)) <: TransferFunction @test eltype(fill(zpk([1],[1,1],1),2)) <: TransferFunction -@test eltype(fill(zpk([],[1,1],1.0),2)) <: TransferFunction +@test eltype(fill(zpk(Int64[],[1,1],1.0),2)) <: TransferFunction @test eltype(fill(zpk([1 2; 3 4],0.001),2)) <: TransferFunction @test eltype(fill(zpk(1)+zpk(2),2)) <: TransferFunction @test eltype(fill(zpk(1,0.005)/zpk(2, 0.005),2)) <: TransferFunction @@ -89,6 +133,8 @@ zpk(tf(vecarray(1, 2, [0], [0]), vecarray(1, 2, [1], [1]), 0.005)) @test eltype([tf(1,1), zpk(1,1)]) <: TransferFunction +zpk(tf([1 2; 3 4])) == zpk([1 2; 3 4]) + # Errors @test_throws ErrorException C_111 + C_222 # Dimension mismatch @test_throws ErrorException C_111 - C_222 # Dimension mismatch @@ -107,4 +153,8 @@ D_diffTs = zpk(tf([1], [2], 0.1)) @test_throws ErrorException zpk("z") # z creation can't be continuous # Remove this when inferec is implemented @test_throws ErrorException [z 0] # Sampling time mismatch (inferec could be implemented) + + +@test typeof(zpk(tf([1], [2], 0.1))) == TransferFunction{ControlSystems.SisoZpk{Float64,Complex{Float64}}} +@test typeof(zpk([-0.5], [], 1)) == TransferFunction{ControlSystems.SisoZpk{Float64,Float64}} end