diff --git a/Project.toml b/Project.toml index a15e2f201..70031fde4 100644 --- a/Project.toml +++ b/Project.toml @@ -27,5 +27,5 @@ IterTools = "1.0" LaTeXStrings = "1.0" OrdinaryDiffEq = "5.2" Plots = "0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 1.0" -Polynomials = "0.6.0, 0.7" +Polynomials = "0.6.0" julia = "1.0" diff --git a/docs/src/man/creating_systems.md b/docs/src/man/creating_systems.md index fc5ab05bc..cfd1a6783 100644 --- a/docs/src/man/creating_systems.md +++ b/docs/src/man/creating_systems.md @@ -10,7 +10,7 @@ The syntax for creating a transfer function is ```julia tf(num, den, Ts=0) ``` -where `num` and `den` are the polinomial coefficients of the numerator and denominator of the polynomial and `Ts` is the sample time. +where `num` and `den` are the polynomial coefficients of the numerator and denominator of the polynomial and `Ts` is the sample time. ### Example: ```julia tf([1.0],[1,2,1]) @@ -29,7 +29,7 @@ Continuous-time transfer function model The transfer functions created using this method will be of type `TransferFunction{SisoRational}`. ## zpk - Pole-Zero-Gain Representation -Sometimes it's better to represent the transferfunction by its poles, zeros and gain, this can be done using +Sometimes it's better to represent the transfer function by its poles, zeros and gain, this can be done using ```julia zpk(zeros, poles, gain, Ts=0) ``` @@ -72,7 +72,7 @@ A state-space system is created using ```julia ss(A,B,C,D,Ts=0) ``` -and they behave similarily to transfer functions. State-space systems with heterogeneous matrix types are also available, which can be used to create systems with static or sized matrices, e.g., +and they behave similarly to transfer functions. State-space systems with heterogeneous matrix types are also available, which can be used to create systems with static or sized matrices, e.g., ```jldoctest HSS using StaticArrays import ControlSystems.HeteroStateSpace diff --git a/src/ControlSystems.jl b/src/ControlSystems.jl index 282ddd919..95c9f2b0f 100644 --- a/src/ControlSystems.jl +++ b/src/ControlSystems.jl @@ -8,7 +8,6 @@ export LTISystem, DelayLtiSystem, Continuous, Discrete, - Static, ss, tf, zpk, @@ -88,10 +87,8 @@ export LTISystem, denvec, numpoly, denpoly, - sampletime, iscontinuous, - isdiscrete, - isstatic + isdiscrete # QUESTION: are these used? LaTeXStrings, Requires, IterTools @@ -112,8 +109,8 @@ abstract type AbstractSystem end include("types/TimeType.jl") ## Added interface: -# sampletime(Lti) -> Number -# sampletype(Lti) -> TimeType +# time(Lti) -> TimeType (not exported) + include("types/Lti.jl") @@ -174,20 +171,6 @@ function covar(D::Union{AbstractMatrix,UniformScaling}, R) D*R*D' end -function Base.getproperty(sys::Union{StateSpace,HeteroStateSpace,TransferFunction}, s::Symbol) - if s === :Ts - # if !isdiscrete(sys) # NOTE this line seems to be breaking inference of isdiscrete - if !isdiscrete(sys.time) # We use this instead - @warn "Getting sampletime 0.0 for non-discrete systems is deprecated. Check `isdiscrete` before trying to access sampletime." - return 0.0 - else - return sampletime(sys) - end - else - return getfield(sys, s) - end -end - # 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 41abc2bd9..299c62385 100644 --- a/src/analysis.jl +++ b/src/analysis.jl @@ -118,7 +118,7 @@ Compute the zeros, poles, and gains of system `sys`. `k` : Matrix{Float64}, (ny x nu)""" function zpkdata(sys::LTISystem) - G = convert(TransferFunction{sampletype(sys),SisoZpk}, sys) + G = convert(TransferFunction{typeof(time(sys)),SisoZpk}, sys) zs = map(x -> x.z, G.matrix) ps = map(x -> x.p, G.matrix) @@ -133,9 +133,8 @@ Compute the natural frequencies, `Wn`, and damping ratios, `zeta`, of the poles, `ps`, of `sys`""" function damp(sys::LTISystem) ps = pole(sys) - if !iscontinuous(sys) - Ts = isstatic(sys) ? 1 : sampletime(sys) - ps = log(ps)/Ts + if isdiscrete(sys) + ps = log(ps)/sys.Ts end Wn = abs.(ps) order = sortperm(Wn; by=z->(abs(z), real(z), imag(z))) @@ -447,7 +446,7 @@ function delaymargin(G::LTISystem) ωϕₘ = m[3][i] dₘ = ϕₘ/ωϕₘ if isdiscrete(G) - dₘ /= sampletime(G) # Give delay margin in number of sample times, as matlab does + dₘ /= G.Ts # Give delay margin in number of sample times, as matlab does end dₘ end diff --git a/src/connections.jl b/src/connections.jl index c8c55b51f..55dd8e701 100644 --- a/src/connections.jl +++ b/src/connections.jl @@ -22,18 +22,18 @@ Append systems in block diagonal form """ function append(systems::(ST where ST<:AbstractStateSpace)...) ST = promote_type(typeof.(systems)...) - Ts = common_sample_time(s.time for s in systems) + time = common_time(systems...) A = blockdiag([s.A for s in systems]...) B = blockdiag([s.B for s in systems]...) C = blockdiag([s.C for s in systems]...) D = blockdiag([s.D for s in systems]...) - return ST(A, B, C, D, Ts) + return ST(A, B, C, D, time) end function append(systems::TransferFunction...) - Ts = common_sample_time(s.time for s in systems) + time = common_time(systems...) mat = blockdiag([s.matrix for s in systems]...) - return TransferFunction(mat, Ts) + return TransferFunction(mat, time) end append(systems::LTISystem...) = append(promote(systems...)...) @@ -62,8 +62,8 @@ function Base.vcat(systems::ST...) where ST <: AbstractStateSpace B = vcat([s.B for s in systems]...) C = blockdiag([s.C for s in systems]...) D = vcat([s.D for s in systems]...) - Ts = common_sample_time(s.time for s in systems) - return ST(A, B, C, D, Ts) + time = common_time(systems...) + return ST(A, B, C, D, time) end function Base.vcat(systems::TransferFunction...) @@ -72,10 +72,10 @@ function Base.vcat(systems::TransferFunction...) if !all(s.nu == nu for s in systems) error("All systems must have same input dimension") end - Ts = common_sample_time(s.time for s in systems) + time = common_time(systems...) mat = vcat([s.matrix for s in systems]...) - return TransferFunction(mat, Ts) + return TransferFunction(mat, time) end Base.vcat(systems::LTISystem...) = vcat(promote(systems...)...) @@ -86,13 +86,13 @@ function Base.hcat(systems::ST...) where ST <: AbstractStateSpace if !all(s.ny == ny for s in systems) error("All systems must have same output dimension") end - Ts = common_sample_time(s.time for s in systems) + time = common_time(systems...) A = blockdiag([s.A for s in systems]...) B = blockdiag([s.B for s in systems]...) C = hcat([s.C for s in systems]...) D = hcat([s.D for s in systems]...) - return ST(A, B, C, D, Ts) + return ST(A, B, C, D, time) end function Base.hcat(systems::TransferFunction...) @@ -101,10 +101,10 @@ function Base.hcat(systems::TransferFunction...) if !all(s.ny == ny for s in systems) error("All systems must have same output dimension") end - Ts = common_sample_time(s.time for s in systems) + time = common_time(systems...) mat = hcat([s.matrix for s in systems]...) - return TransferFunction(mat, Ts) + return TransferFunction(mat, time) end Base.hcat(systems::LTISystem...) = hcat(promote(systems...)...) @@ -201,13 +201,13 @@ function feedback(sys::Union{StateSpace, DelayLtiSystem}) end function feedback(sys1::StateSpace,sys2::StateSpace) - Ts = common_sample_time(sys1.time,sys2.time) + time = common_time(sys1,sys2) !(iszero(sys1.D) || iszero(sys2.D)) && error("There cannot be a direct term (D) in both sys1 and sys2") A = [sys1.A+sys1.B*(-sys2.D)*sys1.C sys1.B*(-sys2.C); sys2.B*sys1.C sys2.A+sys2.B*sys1.D*(-sys2.C)] B = [sys1.B; sys2.B*sys1.D] C = [sys1.C sys1.D*(-sys2.C)] - ss(A, B, C, sys1.D, Ts) + ss(A, B, C, sys1.D, time) end @@ -230,7 +230,7 @@ See Zhou etc. for similar (somewhat less symmetric) formulas. U1=:, Y1=:, U2=:, Y2=:, W1=:, Z1=:, W2=Int[], Z2=Int[], Wperm=:, Zperm=:, pos_feedback::Bool=false) - Ts = common_sample_time(sys1.time,sys2.time) + time = common_time(sys1,sys2) if !(isa(Y1, Colon) || allunique(Y1)); @warn "Connecting single output to multiple inputs Y1=$Y1"; end if !(isa(Y2, Colon) || allunique(Y2)); @warn "Connecting single output to multiple inputs Y2=$Y2"; end @@ -299,7 +299,7 @@ See Zhou etc. for similar (somewhat less symmetric) formulas. s2_D12*R2*s1_D21 s2_D11 + α*s2_D12*R2*s1_D22*s2_D21] end - return StateSpace(A, B[:, Wperm], C[Zperm,:], D[Zperm, Wperm], Ts) + return StateSpace(A, B[:, Wperm], C[Zperm,:], D[Zperm, Wperm], time) end diff --git a/src/delay_systems.jl b/src/delay_systems.jl index 3dc18b2d0..c6cd203b3 100644 --- a/src/delay_systems.jl +++ b/src/delay_systems.jl @@ -161,7 +161,7 @@ function _bounds_and_features(sys::DelayLtiSystem, plot::Val) end # Againm we have to do something for default vectors, more or less a copy from timeresp.jl -function _default_Ts(sys::DelayLtiSystem) +function _default_dt(sys::DelayLtiSystem) if !isstable(sys.P.P) return 0.05 # Something small else diff --git a/src/discrete.jl b/src/discrete.jl index 10cc172e0..617944cf9 100644 --- a/src/discrete.jl +++ b/src/discrete.jl @@ -222,7 +222,7 @@ function lsima(sys::StateSpace, t::AbstractVector, r::AbstractVector, control_si if iscontinuous(sys) dsys = c2d(sys, dt, :zoh)[1] else - if sampletime(sys) != dt + if sys.Ts != dt error("Time vector must match sample time for discrete system") end dsys = sys diff --git a/src/freqresp.jl b/src/freqresp.jl index 19bab97e6..fad4c63d7 100644 --- a/src/freqresp.jl +++ b/src/freqresp.jl @@ -7,11 +7,10 @@ Evaluate the frequency response of a linear system of system `sys` over the frequency vector `w`.""" function freqresp(sys::LTISystem, w_vec::AbstractVector{<:Real}) # Create imaginary freq vector s - if !iscontinuous(sys) - Ts = isstatic(sys) ? 1.0 : sampletime(sys) - s_vec = exp.(w_vec*(im*Ts)) - else + if iscontinuous(sys) s_vec = im*w_vec + else + s_vec = exp.(w_vec*(im*sys.Ts)) end if isa(sys, StateSpace) sys = _preprocess_for_freqresp(sys) @@ -65,7 +64,7 @@ end 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 exp(i*Ts*omega) +- F(omega,true) evaluates the discrete-time transfer function F at exp(im*Ts*omega) - F(z,false) evaluates the discrete-time transfer function F at z """ function (sys::TransferFunction)(s) @@ -75,7 +74,7 @@ end function (sys::TransferFunction)(z_or_omega::Number, map_to_unit_circle::Bool) @assert isdiscrete(sys) "It only makes no sense to call this function with discrete systems" if map_to_unit_circle - isreal(z_or_omega) ? evalfr(sys,exp(im*z_or_omega.*sampletime(sys))) : error("To map to the unit circle, omega should be real") + isreal(z_or_omega) ? evalfr(sys,exp(im*z_or_omega.*sys.Ts)) : error("To map to the unit circle, omega should be real") else evalfr(sys,z_or_omega) end @@ -169,8 +168,8 @@ function _bounds_and_features(sys::LTISystem, plot::Val) w1 = 0.0 w2 = 2.0 end - if !iscontinuous(sys) && !isstatic(sys) # Do not draw above Nyquist freq for disc. systems - w2 = min(w2, log10(π/sampletime(sys))) + if isdiscrete(sys) # Do not draw above Nyquist freq for disc. systems + w2 = min(w2, log10(π/sys.Ts)) end return [w1, w2], zp end diff --git a/src/matrix_comps.jl b/src/matrix_comps.jl index c6b8e7fa5..53abbfd84 100644 --- a/src/matrix_comps.jl +++ b/src/matrix_comps.jl @@ -241,13 +241,8 @@ state space systems in continuous and discrete time', American Control Conferenc See also [`linfnorm`](@ref). """ -function hinfnorm(sys::AbstractStateSpace; tol=1e-6) - if !isdiscrete(sys) # Continuous or Static - return _infnorm_two_steps_ct(sys, :hinf, tol) - else - return _infnorm_two_steps_dt(sys, :hinf, tol) - end -end +hinfnorm(sys::AbstractStateSpace{<:Continuous}; tol=1e-6) = _infnorm_two_steps_ct(sys, :hinf, tol) +hinfnorm(sys::AbstractStateSpace{<:Discrete}; tol=1e-6) = _infnorm_two_steps_dt(sys, :hinf, tol) hinfnorm(sys::TransferFunction; tol=1e-6) = hinfnorm(ss(sys); tol=tol) """ @@ -374,8 +369,8 @@ function _infnorm_two_steps_dt(sys::AbstractStateSpace, normtype::Symbol, tol=1e on_unit_circle = z -> abs(abs(z) - 1) < approxcirc # Helper fcn for readability - T = promote_type(real(numeric_type(sys)), Float64, typeof(true/sampletime(sys))) - Tw = typeof(one(T)/sampletime(sys)) + T = promote_type(real(numeric_type(sys)), Float64, typeof(true/sys.Ts)) + Tw = typeof(one(T)/sys.Ts) if sys.nx == 0 # static gain return (T(opnorm(sys.D)), Tw(0)) @@ -386,7 +381,7 @@ function _infnorm_two_steps_dt(sys::AbstractStateSpace, normtype::Symbol, tol=1e # Check if there is a pole on the unit circle pidx = findfirst(on_unit_circle, pole_vec) if !(pidx isa Nothing) - return T(Inf), Tw(angle(pole_vec[pidx])/sampletime(sys)) + return T(Inf), Tw(angle(pole_vec[pidx])/sys.Ts) end if normtype == :hinf && any(z -> abs(z) > 1, pole_vec) @@ -439,7 +434,7 @@ function _infnorm_two_steps_dt(sys::AbstractStateSpace, normtype::Symbol, tol=1e sort!(θ_vec) if isempty(θ_vec) - return T((1+tol)*lb), Tw(θ_peak/sampletime(sys)) + return T((1+tol)*lb), Tw(θ_peak/sys.Ts) end # Improve the lower bound diff --git a/src/plotting.jl b/src/plotting.jl index dc23a0a86..fc8a84ae1 100644 --- a/src/plotting.jl +++ b/src/plotting.jl @@ -161,15 +161,15 @@ end @userplot Stepplot @userplot Impulseplot """ - stepplot(sys[, Tf[, Ts]]) -Plot step response of `sys` with optional final time `Tf` and discretization time `Ts`. + stepplot(sys[, tfinal[, dt]]) +Plot step response of `sys` with optional final time `tfinal` and discretization time `dt`. If not defined, suitable values are chosen based on `sys`. """ stepplot """ - impulseplot(sys[, Tf[, Ts]]) -Plot step response of `sys` with optional final time `Tf` and discretization time `Ts`. + impulseplot(sys[, tfinal[, dt]]) +Plot step response of `sys` with optional final time `tfinal` and discretization time `dt`. If not defined, suitable values are chosen based on `sys`. """ impulseplot @@ -183,12 +183,12 @@ for (func, title, typ) = ((step, "Step Response", Stepplot), (impulse, "Impulse systems = [systems] end if length(p.args) < 2 - Ts_list, Tf = _default_time_data(systems) + dt_list, tfinal = _default_time_data(systems) elseif length(p.args) == 2 - Ts_list = _default_Ts.(systems) - Tf = p.args[2] + dt_list = _default_dt.(systems) + tfinal = p.args[2] else - Tf, Ts_list = p.args[2:3] + tfinal, dt_list = p.args[2:3] end if !_same_io_dims(systems...) error("All systems must have the same input/output dimensions") @@ -198,8 +198,8 @@ for (func, title, typ) = ((step, "Step Response", Stepplot), (impulse, "Impulse titles = fill("", 1, ny*nu) title --> titles s2i(i,j) = LinearIndices((ny,nu))[i,j] - for (si,(s, Ts)) in enumerate(zip(systems, Ts_list)) - t = 0:Ts:Tf + for (si,(s, dt)) in enumerate(zip(systems, dt_list)) + t = 0:dt:tfinal y = func(s, t)[1] styledict = getStyleSys(si,length(systems)) for i=1:ny @@ -471,7 +471,7 @@ nicholsplot ny, nu = size(systems[1]) if isdiscrete(systems[1]) - w_nyquist = 2π/sampletime(systems[1]) + w_nyquist = 2π/systems[1].Ts w = w[w.<= w_nyquist] end nw = length(w) @@ -706,9 +706,9 @@ function _same_io_dims(systems::LTISystem...) end function _default_time_data(systems::Vector{T}) where T<:LTISystem - sample_times = [_default_Ts(i) for i in systems] - Tf = 100*maximum(sample_times) - return sample_times, Tf + sample_times = [_default_dt(i) for i in systems] + tfinal = 100*maximum(sample_times) + return sample_times, tfinal end _default_time_data(sys::LTISystem) = _default_time_data(LTISystem[sys]) diff --git a/src/timeresp.jl b/src/timeresp.jl index be2bfa48d..d3db3f2c6 100644 --- a/src/timeresp.jl +++ b/src/timeresp.jl @@ -43,14 +43,14 @@ function impulse(sys::StateSpace, t::AbstractVector; method=:cont) lt = length(t) ny, nu = size(sys) nx = sys.nx - if iscontinuous(sys) || isstatic(sys) #&& method == :cont + 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(T, nx, 1), sys.C, zeros(T, ny, 1)) x0s = sys.B else - u = (x,i) -> i == t[1] ? [one(T)]/sampletime(sys) : [zero(T)] + u = (x,i) -> i == t[1] ? [one(T)]/sys.Ts : [zero(T)] imp_sys = sys x0s = zeros(T, nx, nu) end @@ -122,7 +122,7 @@ function lsim(sys::StateSpace, u::AbstractVecOrMat, t::AbstractVector; if !isdiscrete(sys) dsys = c2d(sys, dt, :zoh)[1] else - if sampletime(sys) != dt + if sys.Ts != dt error("Time vector must match sample time for discrete system") end dsys = sys @@ -155,7 +155,7 @@ function lsim(sys::StateSpace, u::Function, t::AbstractVector; if !isdiscrete(sys) dsys = c2d(sys, dt, :zoh)[1] else - if sampletime(sys) != dt + if sys.Ts != dt error("Time vector must match sample time for discrete system") end dsys = sys @@ -217,31 +217,31 @@ end # TODO: This is a poor heuristic to estimate a "good" time vector to use for # simulation, in cases when one isn't provided. -function _default_time_vector(sys::LTISystem, Tf::Real=-1) - Ts = _default_Ts(sys) - if Tf == -1 - Tf = 100*Ts +function _default_time_vector(sys::LTISystem, tfinal::Real=-1) + dt = _default_dt(sys) + if tfinal == -1 + tfinal = 100*dt end - return 0:Ts:Tf + return 0:dt:tfinal end -function _default_Ts(sys::LTISystem) +function _default_dt(sys::LTISystem) if isdiscrete(sys) - Ts = sampletime(sys) + return sys.Ts elseif !isstable(sys) - Ts = 0.05 + return 0.05 else ps = pole(sys) r = minimum([abs.(real.(ps));0]) if r == 0.0 r = 1.0 end - Ts = 0.07/r + return 0.07/r end - return Ts end + #TODO a reasonable check _issmooth(u::Function) = false diff --git a/src/types/DelayLtiSystem.jl b/src/types/DelayLtiSystem.jl index de36e885d..296bce061 100644 --- a/src/types/DelayLtiSystem.jl +++ b/src/types/DelayLtiSystem.jl @@ -16,18 +16,7 @@ struct DelayLtiSystem{T,S<:Real} <: LTISystem # end end -# Fallback since system is always continuous -function getproperty(sys::DelayLtiSystem, s::Symbol) - if s == :Ts - return sys.P.Ts # Will throw deprecation until removed # DEPRECATED - end - return getfield(sys, s) -end - -sampletime(sys::DelayLtiSystem) = sampletime(sys.P) -iscontinuous(sys::DelayLtiSystem) = iscontinuous(sys.P) -isdiscrete(sys::DelayLtiSystem) = isdiscrete(sys.P) -isstatic(sys::DelayLtiSystem) = isstatic(sys.P) +time(sys::DelayLtiSystem) = time(sys.P) # QUESTION: would psys be a good standard variable name for a PartionedStateSpace # and perhaps dsys for a delayed system, (ambigous with discrete system though) @@ -58,8 +47,6 @@ DelayLtiSystem(sys::StateSpace{Continuous,T,MT}) where {T, MT} = DelayLtiSystem{ # From TransferFunction, infer type TODO Use proper constructor instead of convert here when defined DelayLtiSystem(sys::TransferFunction{TimeT,S}) where {TimeT,T,S<:SisoTf{T}} = DelayLtiSystem{T}(convert(StateSpace{Continuous,T, Matrix{T}}, sys)) -# To handle cases where StateSpace is static -DelayLtiSystem(sys::StateSpace{Static,T,MT}, args...) where {T, MT,S} = DelayLtiSystem(StateSpace{Continuous,T,MT}(sys), args...) # TODO: Think through these promotions and conversions Base.promote_rule(::Type{AbstractMatrix{T1}}, ::Type{DelayLtiSystem{T2,S}}) where {T1<:Number,T2<:Number,S} = DelayLtiSystem{promote_type(T1,T2),S} diff --git a/src/types/Lti.jl b/src/types/Lti.jl index 47aea4d99..477cc6cec 100644 --- a/src/types/Lti.jl +++ b/src/types/Lti.jl @@ -41,29 +41,42 @@ function issiso(sys::LTISystem) return ninputs(sys) == 1 && noutputs(sys) == 1 end +"""`time(sys::LTISystem)` +Get the time::TimeType from system `sys`, usually sys.time """ +time(sys::LTISystem) = sys.time """`iscontinuous(sys)` -Returns `true` if `sys` is continuous, else returns `false`.""" -iscontinuous(sys::LTISystem) = iscontinuous(sys.time) +Returns `true` for a continuous-time system `sys`, else returns `false`.""" +iscontinuous(sys::LTISystem) = time(sys) isa Continuous """`isdiscrete(sys)` -Returns `true` if `sys` is discrete, else returns `false`.""" -isdiscrete(sys::LTISystem) = isdiscrete(sys.time) -"""`isstatic(sys)` +Returns `true` for a discrete-time system `sys`, else returns `false`.""" +isdiscrete(sys::LTISystem) = time(sys) isa Discrete + + +function Base.getproperty(sys::LTISystem, s::Symbol) + if s === :Ts + # if !isdiscrete(sys) # NOTE this line seems to be breaking inference of isdiscrete (is there a test for this?) + if isdiscrete(sys) + return time(sys).Ts + else + @warn "Getting time 0.0 for non-discrete systems is deprecated. Check `isdiscrete` before trying to access time." + return 0.0 + end + else + return getfield(sys, s) + end +end + +Base.propertynames(sys::LTISystem, private::Bool=false) = + (fieldnames(typeof(sys))..., (isdiscrete(sys) ? (:Ts,) : ())...) + -Returns `true` if `sys` is static, else returns `false`.""" -isstatic(sys::LTISystem) = isstatic(sys.time) -"""`sampletime(sys)` -Returns the sampletime of a discrete time system, throws error if the system is continuous time or static. -Always ensure `isdiscrete` before using.""" -sampletime(sys::LTISystem) = sampletime(sys.time) +common_time(systems::LTISystem...) = common_time(time(sys) for sys in systems) -"""`sampletype(sys)` -Get the sampletype of system. Usually typeof(sys.time).""" -sampletype(sys) = typeof(sys.time) """`isstable(sys)` diff --git a/src/types/PartionedStateSpace.jl b/src/types/PartionedStateSpace.jl index 8c9fca649..ee267f2f2 100644 --- a/src/types/PartionedStateSpace.jl +++ b/src/types/PartionedStateSpace.jl @@ -11,7 +11,7 @@ u = [u1 u2]^T y = [y1 y2]^T """ -struct PartionedStateSpace{S} +struct PartionedStateSpace{S<:AbstractStateSpace} <: LTISystem P::S nu1::Int ny1::Int @@ -27,8 +27,6 @@ function getproperty(sys::PartionedStateSpace, d::Symbol) if d === :Ts return P.Ts # Will throw deprecation until removed # DEPRECATED - elseif d === :time - return P.time elseif d === :P return P elseif d === :nu1 @@ -58,13 +56,10 @@ function getproperty(sys::PartionedStateSpace, d::Symbol) end end -sampletime(sys::PartionedStateSpace) = sampletime(sys.P) -iscontinuous(sys::PartionedStateSpace) = iscontinuous(sys.P) -isdiscrete(sys::PartionedStateSpace) = isdiscrete(sys.P) -isstatic(sys::PartionedStateSpace) = isstatic(sys.P) +time(sys::PartionedStateSpace) = time(sys.P) function +(s1::PartionedStateSpace, s2::PartionedStateSpace) - Ts = common_sample_time(s1.P.time,s2.P.time) + time = common_time(s1,s2) A = blockdiag(s1.A, s2.A) @@ -76,7 +71,7 @@ function +(s1::PartionedStateSpace, s2::PartionedStateSpace) D = [(s1.D11 + s2.D11) s1.D12 s2.D12; [s1.D21; s2.D21] blockdiag(s1.D22, s2.D22)] - P = StateSpace(A, B, C, D, Ts) # How to handle discrete? + P = StateSpace(A, B, C, D, time) # How to handle discrete? PartionedStateSpace(P, s1.nu1 + s2.nu1, s1.ny1 + s2.ny1) end @@ -88,7 +83,7 @@ end Series connection of partioned StateSpace systems. """ function *(s1::PartionedStateSpace, s2::PartionedStateSpace) - Ts = common_sample_time(s1.P.time,s2.P.time) + time = common_time(s1,s2) A = [s1.A s1.B1*s2.C1; zeros(size(s2.A,1),size(s1.A,2)) s2.A] @@ -104,7 +99,7 @@ function *(s1::PartionedStateSpace, s2::PartionedStateSpace) s1.D21*s2.D11 s1.D22 s1.D21*s2.D12; s2.D21 zeros(size(s2.D22,1),size(s1.D22,2)) s2.D22 ] - P = StateSpace(A, B, C, D, Ts) + P = StateSpace(A, B, C, D, time) PartionedStateSpace(P, s2.nu1, s1.ny1) end @@ -112,7 +107,7 @@ end # QUESTION: What about algebraic loops and well-posedness?! Perhaps issue warning if P1(∞)*P2(∞) > 1 function feedback(s1::PartionedStateSpace, s2::PartionedStateSpace) - Ts = common_sample_time(s1.P.time,s2.P.time) + time = common_time(s1,s2) X_11 = (I + s2.D11*s1.D11)\[-s2.D11*s1.C1 -s2.C1] X_21 = (I + s1.D11*s2.D11)\[s1.C1 -s1.D11*s2.C1] @@ -151,7 +146,7 @@ function feedback(s1::PartionedStateSpace, s2::PartionedStateSpace) #tmp = [blockdiag(s1.D12, s2.D12); blockdiag(s1.D22, s2.D22)] #D[:, end-size(tmp,2)+1:end] .+= tmp - P = StateSpace(A, B, C, D, Ts) + P = StateSpace(A, B, C, D, time) PartionedStateSpace(P, s2.nu1, s1.ny1) end @@ -164,7 +159,7 @@ end """ function vcat_1(systems::PartionedStateSpace...) # Perform checks - Ts = common_sample_time(sys.P.time for sys in systems) + time = common_time(systems...) nu1 = systems[1].nu1 if !all(s.nu1 == nu1 for s in systems) @@ -184,7 +179,7 @@ function vcat_1(systems::PartionedStateSpace...) D21 = reduce(vcat, [s.D21 for s in systems]) D22 = blockdiag([s.D22 for s in systems]...) - sysnew = StateSpace(A, [B1 B2], [C1; C2], [D11 D12; D21 D22], Ts) + sysnew = StateSpace(A, [B1 B2], [C1; C2], [D11 D12; D21 D22], time) return PartionedStateSpace(sysnew, nu1, sum(s -> s.ny1, systems)) end @@ -198,7 +193,7 @@ end """ function hcat_1(systems::PartionedStateSpace...) # Perform checks - Ts = common_sample_time(sys.P.time for sys in systems) + time = common_time(systems...) ny1 = systems[1].ny1 if !all(s.ny1 == ny1 for s in systems) @@ -218,7 +213,7 @@ function hcat_1(systems::PartionedStateSpace...) D21 = blockdiag([s.D21 for s in systems]...) D22 = blockdiag([s.D22 for s in systems]...) - sysnew = StateSpace(A, [B1 B2], [C1; C2], [D11 D12; D21 D22], Ts) + sysnew = StateSpace(A, [B1 B2], [C1; C2], [D11 D12; D21 D22], time) return PartionedStateSpace(sysnew, sum(s -> s.nu1, systems), ny1) end diff --git a/src/types/StateSpace.jl b/src/types/StateSpace.jl index 330afb78a..4bb63449f 100644 --- a/src/types/StateSpace.jl +++ b/src/types/StateSpace.jl @@ -2,7 +2,7 @@ ## Data Type Declarations ## ##################################################################### -function state_space_validation(A,B,C,D,Ts::TimeType) +function state_space_validation(A,B,C,D,time::TimeType) nx = size(A, 1) nu = size(B, 2) ny = size(C, 1) @@ -22,8 +22,9 @@ function state_space_validation(A,B,C,D,Ts::TimeType) nx,nu,ny end -abstract type AbstractStateSpace <: LTISystem end -struct StateSpace{TimeT<:TimeType, T, MT<:AbstractMatrix{T}} <: AbstractStateSpace +abstract type AbstractStateSpace{TimeT<:TimeType} <: LTISystem end + +struct StateSpace{TimeT, T, MT<:AbstractMatrix{T}} <: AbstractStateSpace{TimeT} A::MT B::MT C::MT @@ -32,13 +33,13 @@ struct StateSpace{TimeT<:TimeType, T, MT<:AbstractMatrix{T}} <: AbstractStateSpa nx::Int nu::Int ny::Int - function StateSpace{TimeT, T, MT}(A::MT, B::MT, C::MT, D::MT, Ts::TimeT) where {TimeT<:TimeType, T, MT <: AbstractMatrix{T}} - nx,nu,ny = state_space_validation(A,B,C,D,Ts) - new{TimeT, T, MT}(A, B, C, D, Ts, nx, nu, ny) + function StateSpace{TimeT, T, MT}(A::MT, B::MT, C::MT, D::MT, time::TimeT) where {TimeT<:TimeType, T, MT <: AbstractMatrix{T}} + nx,nu,ny = state_space_validation(A,B,C,D,time) + new{TimeT, T, MT}(A, B, C, D, time, nx, nu, ny) end end -function StateSpace(A::MT, B::MT, C::MT, D::MT, Ts::TimeT) where {TimeT<:TimeType, T, MT <: AbstractMatrix{T}} - StateSpace{TimeT, T, MT}(A, B, C, D, Ts) +function StateSpace(A::MT, B::MT, C::MT, D::MT, time::TimeT) where {TimeT<:TimeType, T, MT <: AbstractMatrix{T}} + StateSpace{TimeT, T, MT}(A, B, C, D, time) end # Constructor for Discrete system function StateSpace(A::MT, B::MT, C::MT, D::MT, Ts::Number) where {T, MT <: AbstractMatrix{T}} @@ -71,9 +72,9 @@ end const AbstractNumOrArray = Union{Number, AbstractVecOrMat} # Explicit construction with different types -function StateSpace{TimeT,T,MT}(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, Ts::TimeType) where {TimeT, T, MT <: AbstractMatrix{T}} +function StateSpace{TimeT,T,MT}(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, time::TimeType) where {TimeT, T, MT <: AbstractMatrix{T}} D = fix_D_matrix(T, B, C, D) - return StateSpace{TimeT, T,Matrix{T}}(MT(to_matrix(T, A)), MT(to_matrix(T, B)), MT(to_matrix(T, C)), MT(D), TimeT(Ts)) + return StateSpace{TimeT, T,Matrix{T}}(MT(to_matrix(T, A)), MT(to_matrix(T, B)), MT(to_matrix(T, C)), MT(D), TimeT(time)) end # Explicit conversion @@ -81,9 +82,9 @@ function StateSpace{TimeT,T,MT}(sys::StateSpace) where {TimeT, T, MT <: Abstract StateSpace{TimeT,T,MT}(sys.A,sys.B,sys.C,sys.D,sys.time) end -function StateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, Ts::TimeType) +function StateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, time::TimeType) A, B, C, D, T = to_similar_matrices(A,B,C,D) - return StateSpace{typeof(Ts),T,Matrix{T}}(A, B, C, D, Ts) + return StateSpace{typeof(time),T,Matrix{T}}(A, B, C, D, time) end # General Discrete constructor StateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, Ts::Number) = @@ -93,19 +94,19 @@ StateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, StateSpace(A, B, C, D, Continuous()) # Function for creation of static gain -function StateSpace(D::AbstractArray{T}, Ts::TimeType) where {T<:Number} +function StateSpace(D::AbstractArray{T}, time::TimeType) where {T<:Number} ny, nu = size(D, 1), size(D, 2) A = fill(zero(T), 0, 0) B = fill(zero(T), 0, nu) C = fill(zero(T), ny, 0) D = reshape(D, (ny,nu)) - return StateSpace(A, B, C, D, Ts) + return StateSpace(A, B, C, D, time) end StateSpace(D::AbstractArray, Ts::Number) = StateSpace(D, Discrete(Ts)) -StateSpace(D::AbstractArray) = StateSpace(D, Static()) +StateSpace(D::AbstractArray) = StateSpace(D, Continuous()) StateSpace(d::Number, Ts::Number; kwargs...) = StateSpace([d], Discrete(Ts)) -StateSpace(d::Number; kwargs...) = StateSpace([d], Static()) +StateSpace(d::Number; kwargs...) = StateSpace([d], Continuous()) # StateSpace(sys) converts to StateSpace @@ -117,7 +118,7 @@ StateSpace(sys::LTISystem) = convert(StateSpace, sys) Create a state-space model `sys::StateSpace{TimeT, T, MT<:AbstractMatrix{T}}` where `MT` is the type of matrixes `A,B,C,D`, `T` the element type and TimeT is Continuous or Discrete. -This is a continuous-time model if Ts is omitted. +This is a continuous-time model if `Ts` is omitted. Otherwise, this is a discrete-time model with sampling period Ts. `sys = ss(D [, Ts])` specifies a static gain matrix D. @@ -125,7 +126,7 @@ Otherwise, this is a discrete-time model with sampling period Ts. ss(args...;kwargs...) = StateSpace(args...;kwargs...) -struct HeteroStateSpace{TimeT<:TimeType, AT<:AbstractMatrix,BT<:AbstractMatrix,CT<:AbstractMatrix,DT<:AbstractMatrix} <: AbstractStateSpace +struct HeteroStateSpace{TimeT, AT<:AbstractMatrix,BT<:AbstractMatrix,CT<:AbstractMatrix,DT<:AbstractMatrix} <: AbstractStateSpace{TimeT} A::AT B::BT C::CT @@ -136,21 +137,21 @@ struct HeteroStateSpace{TimeT<:TimeType, AT<:AbstractMatrix,BT<:AbstractMatrix,C ny::Int end function HeteroStateSpace(A::AT, B::BT, - C::CT, D::DT, Ts::TimeT) where {TimeT<:TimeType,AT<:AbstractMatrix,BT<:AbstractMatrix,CT<:AbstractMatrix,DT<:AbstractMatrix} - nx,nu,ny = state_space_validation(A,B,C,D,Ts) - HeteroStateSpace{TimeT,AT,BT,CT,DT}(A, B, C, D, Ts, nx, nu, ny) + C::CT, D::DT, time::TimeT) where {TimeT<:TimeType,AT<:AbstractMatrix,BT<:AbstractMatrix,CT<:AbstractMatrix,DT<:AbstractMatrix} + nx,nu,ny = state_space_validation(A,B,C,D,time) + HeteroStateSpace{TimeT,AT,BT,CT,DT}(A, B, C, D, time, nx, nu, ny) end # Explicit constructor -function HeteroStateSpace{TimeT,AT,BT,CT,DT}(A, B, C, D, Ts) where {TimeT,AT,BT,CT,DT} - nx,nu,ny = state_space_validation(A,B,C,D,Ts) - HeteroStateSpace{TimeT,AT,BT,CT,DT}(AT(A), BT(B), CT(C), DT(D), TimeT(Ts), nx, nu, ny) +function HeteroStateSpace{TimeT,AT,BT,CT,DT}(A, B, C, D, time) where {TimeT,AT,BT,CT,DT} + nx,nu,ny = state_space_validation(A,B,C,D,time) + HeteroStateSpace{TimeT,AT,BT,CT,DT}(AT(A), BT(B), CT(C), DT(D), TimeT(time), nx, nu, ny) end HeteroStateSpace(s::AbstractStateSpace) = HeteroStateSpace(s.A,s.B,s.C,s.D,s.time) # Base constructor -function HeteroStateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, Ts::TimeType) +function HeteroStateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, time::TimeType) A = to_abstract_matrix(A) B = to_abstract_matrix(B) C = to_abstract_matrix(C) @@ -159,7 +160,7 @@ function HeteroStateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::Abstr else D = to_abstract_matrix(D) end - return HeteroStateSpace{typeof(Ts),typeof(A),typeof(B),typeof(C),typeof(D)}(A, B, C, D, Ts) + return HeteroStateSpace{typeof(time),typeof(A),typeof(B),typeof(C),typeof(D)}(A, B, C, D, time) end HeteroStateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrArray, D::AbstractNumOrArray, Ts::Number) = @@ -169,20 +170,20 @@ HeteroStateSpace(A::AbstractNumOrArray, B::AbstractNumOrArray, C::AbstractNumOrA # Function for creation of static gain -function HeteroStateSpace(D::AbstractArray{T}, Ts::TimeType) where {T<:Number} +function HeteroStateSpace(D::AbstractArray{T}, time::TimeType) where {T<:Number} ny, nu = size(D, 1), size(D, 2) A = fill(zero(T), 0, 0) B = fill(zero(T), 0, nu) C = fill(zero(T), ny, 0) - return HeteroStateSpace(A, B, C, D, Ts) + return HeteroStateSpace(A, B, C, D, time) end HeteroStateSpace(D::AbstractArray{T}, Ts::Number) where {T<:Number} = HeteroStateSpace(D, Discrete(Ts)) -HeteroStateSpace(D::AbstractArray{T}) where {T<:Number} = HeteroStateSpace(D, Static()) +HeteroStateSpace(D::AbstractArray{T}) where {T<:Number} = HeteroStateSpace(D, Continuous()) HeteroStateSpace(d::Number, Ts::Number; kwargs...) = HeteroStateSpace([d], Discrete(Ts)) -HeteroStateSpace(d::Number; kwargs...) = HeteroStateSpace([d], Static()) +HeteroStateSpace(d::Number; kwargs...) = HeteroStateSpace([d], Continuous()) # HeteroStateSpace(sys) converts to HeteroStateSpace HeteroStateSpace(sys::LTISystem) = convert(HeteroStateSpace, sys) @@ -220,7 +221,7 @@ function +(s1::StateSpace{TimeT,T,MT}, s2::StateSpace{TimeT,T,MT}) where {TimeT, if size(s1) != size(s2) error("Systems have different shapes.") end - Ts = common_sample_time(s1.time,s2.time) + time = common_time(s1,s2) A = [s1.A fill(zero(T), nstates(s1), nstates(s2)); fill(zero(T), nstates(s2), nstates(s1)) s2.A] @@ -228,7 +229,7 @@ function +(s1::StateSpace{TimeT,T,MT}, s2::StateSpace{TimeT,T,MT}) where {TimeT, C = [s1.C s2.C;] D = [s1.D + s2.D;] - return StateSpace{TimeT,T,MT}(A, B, C, D, Ts) + return StateSpace{TimeT,T,MT}(A, B, C, D, time) end function +(s1::HeteroStateSpace, s2::HeteroStateSpace) @@ -236,7 +237,7 @@ function +(s1::HeteroStateSpace, s2::HeteroStateSpace) if size(s1) != size(s2) error("Systems have different shapes.") end - Ts = common_sample_time(s1.time,s2.time) + time = common_time(s1,s2) T = promote_type(eltype(s1.A),eltype(s2.A)) A = [s1.A fill(zero(T), nstates(s1), nstates(s2)); fill(zero(T), nstates(s2), nstates(s1)) s2.A] @@ -244,7 +245,7 @@ function +(s1::HeteroStateSpace, s2::HeteroStateSpace) C = [s1.C s2.C;] D = [s1.D + s2.D;] - return HeteroStateSpace(A, B, C, D, Ts) + return HeteroStateSpace(A, B, C, D, time) end +(sys::ST, n::Number) where ST <: AbstractStateSpace = ST(sys.A, sys.B, sys.C, sys.D .+ n, sys.time) @@ -265,14 +266,14 @@ function *(sys1::StateSpace{TimeT,T,MT}, sys2::StateSpace{TimeT,T,MT}) where {Ti if sys1.nu != sys2.ny error("sys1*sys2: sys1 must have same number of inputs as sys2 has outputs") end - Ts = common_sample_time(sys1.time,sys2.time) + time = common_time(sys1,sys2) 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{TimeT,T,MT}(A, B, C, D, Ts) + return StateSpace{TimeT,T,MT}(A, B, C, D, time) end function *(sys1::HeteroStateSpace, sys2::HeteroStateSpace) @@ -281,7 +282,7 @@ function *(sys1::HeteroStateSpace, sys2::HeteroStateSpace) if sys1.nu != sys2.ny error("sys1*sys2: sys1 must have same number of inputs as sys2 has outputs") end - Ts = common_sample_time(sys1.time,sys2.time) + time = common_time(sys1,sys2) T = promote_type(eltype(sys1.A),eltype(sys2.A)) A = [sys1.A sys1.B*sys2.C; fill(zero(T), sys2.nx, sys1.nx) sys2.A] @@ -289,7 +290,7 @@ function *(sys1::HeteroStateSpace, sys2::HeteroStateSpace) C = [sys1.C sys1.D*sys2.C;] D = [sys1.D*sys2.D;] - return HeteroStateSpace(A, B, C, D, Ts) + return HeteroStateSpace(A, B, C, D, time) end *(sys::ST, n::Number) where ST <: AbstractStateSpace = StateSpace(sys.A, sys.B, sys.C*n, sys.D*n, sys.time) @@ -360,12 +361,10 @@ function Base.show(io::IO, sys::AbstractStateSpace) println(io, "D = \n", _string_mat_with_headers(sys.D), "\n") # Print sample time if isdiscrete(sys) - println(io, "Sample Time: ", sampletime(sys), " (seconds)") + println(io, "Sample Time: ", sys.Ts, " (seconds)") end # Print model type - if isstatic(sys) - print(io, "Static gain state-space model") - elseif iscontinuous(sys) + if iscontinuous(sys) print(io, "Continuous-time state-space model") else print(io, "Discrete-time state-space model") diff --git a/src/types/TimeType.jl b/src/types/TimeType.jl index 68385a951..46dfd30c0 100644 --- a/src/types/TimeType.jl +++ b/src/types/TimeType.jl @@ -1,12 +1,12 @@ abstract type TimeType end -const UNDEF_SAMPLETIME = -1 +const UNDEF_SAMPLEPERIOD = -1 # For handling promotion of Matrix to LTISystem struct Discrete{T} <: TimeType Ts::T function Discrete{T}(Ts::T) where T - if Ts <= 0 && Ts != UNDEF_SAMPLETIME + if Ts <= 0 && Ts != UNDEF_SAMPLEPERIOD throw(ErrorException("Creating a continuous time system by setting sample time to 0 is no longer supported.")) end new{T}(Ts) @@ -15,77 +15,49 @@ end Discrete{T}(x) where T = Discrete{T}(T(x)) struct Continuous <: TimeType end -struct Static <: TimeType end # Basic Constructors Discrete(x::T) where T<:Number = Discrete{T}(x) Continuous(x::Continuous) = x -Static(x::Static) = x # Simple converseion Discrete{T}(x::Discrete) where T = Discrete{T}(x.Ts) -Discrete{T}(::Static) where T = Discrete{T}(UNDEF_SAMPLETIME) -Continuous(::Static) = Continuous() -# Default to checkng type -iscontinuous(x::TimeType) = x isa Continuous -isdiscrete(x::TimeType) = x isa Discrete -isstatic(x::TimeType) = x isa Static -# Interface for getting sample time -sampletime(x::Discrete) = x.Ts -sampletime(x::Continuous) = error("Continuous system has no sample time") -sampletime(x::Static) = error("Static system has no sample time") -unknown_time(::Type{Discrete{T}}) where T = Discrete{T}(UNDEF_SAMPLETIME) -unknown_time(::Type{Continuous}) where T = Continuous() -unknown_time(::Type{Static}) where T = Static() +undef_sampleperiod(::Type{Discrete{T}}) where T = Discrete{T}(UNDEF_SAMPLEPERIOD) +undef_sampleperiod(::Type{Continuous}) where T = Continuous() + # Promotion # Fallback to give useful error Base.promote_rule(::Type{<:Discrete}, ::Type{Continuous}) = throw(ErrorException("Sampling time mismatch")) -Base.promote_rule(::Type{T}, ::Type{Static}) where T <: TimeType = T Base.promote_rule(::Type{Discrete{T1}}, ::Type{Discrete{T2}}) where {T1,T2}= Discrete{promote_type(T1,T2)} Base.convert(::Type{Discrete{T1}}, x::Discrete{T2}) where {T1,T2} = Discrete{T1}(x.Ts) # Promoting two or more systems systems should promote sample times -common_sample_time(x::TimeType) = x -common_sample_time(x::TimeType, y::TimeType) = throw(ErrorException("Sampling time mismatch")) -common_sample_time(x::TimeType, y::TimeType, z...) = common_sample_time(common_sample_time(x, y), z...) -common_sample_time(a::Base.Generator) = reduce(common_sample_time, a) - -function common_sample_time(x::Discrete{T1}, y::Discrete{T2}) where {T1,T2} - if x != y && x.Ts != UNDEF_SAMPLETIME && y.Ts != UNDEF_SAMPLETIME - throw(ErrorException("Sampling time mismatch")) +common_time(x::TimeType) = x +common_time(x::TimeType, y::TimeType) = throw(ErrorException("Sampling time mismatch")) +common_time(x::TimeType, y::TimeType, z...) = common_time(common_time(x, y), z...) +common_time(a::Base.Generator) = reduce(common_time, a) + +function common_time(x::Discrete{T1}, y::Discrete{T2}) where {T1,T2} + if x != y && x.Ts != UNDEF_SAMPLEPERIOD && y.Ts != UNDEF_SAMPLEPERIOD + throw(ErrorException("Sampling time mismatch")) end - if x.Ts == UNDEF_SAMPLETIME + if x.Ts == UNDEF_SAMPLEPERIOD return Discrete{promote_type(T1,T2)}(y) else return Discrete{promote_type(T1,T2)}(x) end end -common_sample_time(x::Continuous, ys::Continuous...) = Continuous() -common_sample_time(x::Static, ys::Static...) = Static() - -# Combinations with static -common_sample_time(x::TimeType, ys::Static) = x -common_sample_time(x::Static, y::TimeType) = y +common_time(x::Continuous, ys::Continuous...) = Continuous() -# Check equality and similarity +# Check equality ==(x::TimeType, y::TimeType) = false ==(x::Discrete, y::Discrete) = (x.Ts == y.Ts) ==(::Continuous, ::Continuous) = true -==(::Static, ::Static) = true - -# Equality with static is true -==(::Static, ::TimeType) = true -==(::TimeType, ::Static) = true isapprox(x::TimeType, y::TimeType, args...; kwargs...) = false isapprox(x::Discrete, y::Discrete, args...; kwargs...) = isapprox(x.Ts, y.Ts, args...; kwargs...) isapprox(::Continuous, ::Continuous, args...; kwargs...) = true -isapprox(::Static, ::Static, args...; kwargs...) = true - -# Equality with static is true -isapprox(::Static, ::TimeType, args...; kwargs...) = true -isapprox(::TimeType, ::Static, args...; kwargs...) = true diff --git a/src/types/TransferFunction.jl b/src/types/TransferFunction.jl index a3787e120..5967cdd33 100644 --- a/src/types/TransferFunction.jl +++ b/src/types/TransferFunction.jl @@ -1,16 +1,16 @@ -struct TransferFunction{TimeT<:TimeType, S<:SisoTf{T} where T} <: LTISystem +struct TransferFunction{TimeT, S<:SisoTf{T} where T} <: LTISystem matrix::Matrix{S} time::TimeT nu::Int ny::Int - function TransferFunction{TimeT,S}(matrix::Matrix{S}, Ts::TimeT) where {S,TimeT} + function TransferFunction{TimeT,S}(matrix::Matrix{S}, time::TimeT) where {S,TimeT} # Validate size of input and output names ny, nu = size(matrix) - return new{TimeT,S}(matrix, Ts, nu, ny) + return new{TimeT,S}(matrix, time, nu, ny) end end -function TransferFunction(matrix::Matrix{S}, Ts::TimeT) where {TimeT<:TimeType, T<:Number, S<:SisoTf{T}} - TransferFunction{TimeT, S}(matrix, Ts) +function TransferFunction(matrix::Matrix{S}, time::TimeT) where {TimeT<:TimeType, T<:Number, S<:SisoTf{T}} + TransferFunction{TimeT, S}(matrix, time) end # # Constructor for Discrete time system @@ -22,6 +22,7 @@ end # return TransferFunction(matrix, Continuous()) # end + noutputs(G::TransferFunction) = size(G.matrix, 1) ninputs(G::TransferFunction) = size(G.matrix, 2) @@ -110,9 +111,9 @@ function +(G1::TransferFunction, G2::TransferFunction) if size(G1) != size(G2) error("Systems have different shapes.") end - Ts = common_sample_time(G1.time,G2.time) + time = common_time(G1,G2) matrix = G1.matrix + G2.matrix - return TransferFunction(matrix, Ts) + return TransferFunction(matrix, time) end +(G::TransferFunction, n::Number) = TransferFunction(G.matrix .+ n, G.time) @@ -130,13 +131,12 @@ end function *(G1::TransferFunction, G2::TransferFunction) # Note: G1*G2 = y <- G1 <- G2 <- u - Ts = common_sample_time(G1.time,G2.time) + time = common_time(G1,G2) if G1.nu != G2.ny error("G1*G2: G1 must have same number of inputs as G2 has outputs") end - Ts = common_sample_time(G1.time,G2.time) matrix = G1.matrix * G2.matrix - return TransferFunction(matrix, Ts) + return TransferFunction(matrix, time) end *(G::TransferFunction, n::Number) = TransferFunction(n*G.matrix, G.time) @@ -182,7 +182,7 @@ function Base.show(io::IO, G::TransferFunction) if iscontinuous(G) print(io, "\nContinuous-time transfer function model") elseif isdiscrete(G) - print(io, "\nSample Time: ", sampletime(G), " (seconds)") + print(io, "\nSample Time: ", G.Ts, " (seconds)") print(io, "\nDiscrete-time transfer function model") else print(io, "\nStatic gain transfer function model") diff --git a/src/types/conversion.jl b/src/types/conversion.jl index a144dd9cd..6e4e49665 100644 --- a/src/types/conversion.jl +++ b/src/types/conversion.jl @@ -20,9 +20,9 @@ # Base.convert(::Type{<:TransferFunction{<:SisoZpk}}, b::Number) = zpk(b) # Base.convert(::Type{TransferFunction{TimeT,SisoZpk{T1, TR1}}}, b::AbstractMatrix{T2}) where {TimeT, T1, TR1, T2<:Number} = - zpk(T1.(b), unknown_time(TimeT)) + zpk(T1.(b), undef_sampleperiod(TimeT)) Base.convert(::Type{TransferFunction{TimeT,SisoRational{T1}}}, b::AbstractMatrix{T2}) where {TimeT, T1, T2<:Number} = - tf(T1.(b), unknown_time(TimeT)) + tf(T1.(b), undef_sampleperiod(TimeT)) function convert(::Type{StateSpace{TimeT,T,MT}}, D::AbstractMatrix{<:Number}) where {TimeT,T, MT} (ny, nu) = size(D) @@ -30,14 +30,14 @@ function convert(::Type{StateSpace{TimeT,T,MT}}, D::AbstractMatrix{<:Number}) wh B = MT(fill(zero(T), (0,nu))) C = MT(fill(zero(T), (ny,0))) D = convert(MT, D) - return StateSpace{TimeT,T,MT}(A,B,C,D,unknown_time(TimeT)) + return StateSpace{TimeT,T,MT}(A,B,C,D,undef_sampleperiod(TimeT)) end # TODO We still dont have TransferFunction{..}(number/array) constructors Base.convert(::Type{TransferFunction{TimeT,SisoRational{T}}}, b::Number) where {TimeT, T} = - tf(T(b), unknown_time(TimeT)) + tf(T(b), undef_sampleperiod(TimeT)) Base.convert(::Type{TransferFunction{TimeT,SisoZpk{T,TR}}}, b::Number) where {TimeT, T, TR} = - zpk(T(b), unknown_time(TimeT)) + zpk(T(b), undef_sampleperiod(TimeT)) Base.convert(::Type{StateSpace{TimeT,T,MT}}, b::Number) where {TimeT, T, MT} = convert(StateSpace{TimeT,T,MT}, fill(b, (1,1))) # diff --git a/src/types/tf.jl b/src/types/tf.jl index ef39c8f48..ce1d6256f 100644 --- a/src/types/tf.jl +++ b/src/types/tf.jl @@ -1,4 +1,4 @@ -""" `sys = tf(num, den[, Ts=0]), sys = tf(gain[, Ts])` +""" `sys = tf(num, den[, Ts]), sys = tf(gain[, Ts])` Create as a fraction of polynomials: @@ -40,9 +40,9 @@ tf(num::AbstractVecOrMat{<:AbstractVector{T1}}, den::AbstractVecOrMat{<:Abstract tf(num::AbstractVecOrMat{<:AbstractVector{T1}}, den::AbstractVecOrMat{<:AbstractVector{T2}}) where {T1,T2} = tf(num, den, Continuous()) -function tf(num::AbstractVector{T1}, den::AbstractVector{T2}, Ts::TimeT) where {TimeT<:TimeType,T1<:Number, T2<:Number} +function tf(num::AbstractVector{T1}, den::AbstractVector{T2}, time::TimeT) where {TimeT<:TimeType,T1<:Number, T2<:Number} T = promote_type(T1, T2) - return TransferFunction{TimeT,SisoRational{T}}(fill(SisoRational{T}(num, den), 1, 1), Ts) + return TransferFunction{TimeT,SisoRational{T}}(fill(SisoRational{T}(num, den), 1, 1), time) end tf(num::AbstractVector{T1}, den::AbstractVector{T2}, Ts::Number) where {T1<:Number, T2<:Number} = tf(num, den, Discrete(Ts)) @@ -52,17 +52,17 @@ tf(num::AbstractVector{T1}, den::AbstractVector{T2}) where {T1<:Number, T2<:Numb tf(num::Number, den::Vector, args...) = tf([num], den, args...) # Cases for just static gain -function tf(D::AbstractArray{T}, Ts::TimeT) where {TimeT<:TimeType, T<:Number} +function tf(D::AbstractArray{T}, time::TimeT) where {TimeT<:TimeType, T<:Number} ny, nu = size(D, 1), size(D, 2) matrix = Matrix{SisoRational{T}}(undef, ny, nu) for i in eachindex(D) matrix[i] = SisoRational{T}([D[i]], [one(T)]) end - return TransferFunction{TimeT,SisoRational{T}}(matrix, Ts) + return TransferFunction{TimeT,SisoRational{T}}(matrix, time) end tf(D::AbstractArray{T}, Ts::Number) where T = tf(D, Discrete(Ts)) -tf(D::AbstractArray{T}) where T = tf(D, Static()) +tf(D::AbstractArray{T}) where T = tf(D, Continuous()) tf(n::Number, args...; kwargs...) = tf([n], args...; kwargs...) @@ -84,7 +84,7 @@ function tf(var::AbstractString, Ts::Real) end ## Constructors for polynomial inputs -function tf(num::AbstractArray{PT}, den::AbstractArray{PT}, Ts::TimeT) where {TimeT<:TimeType,T<:Number, PT <: Polynomials.Poly{T}} +function tf(num::AbstractArray{PT}, den::AbstractArray{PT}, time::TimeT) where {TimeT<:TimeType,T<:Number, PT <: Polynomials.Poly{T}} ny, nu = size(num, 1), size(num, 2) if (ny, nu) != (size(den, 1), size(den, 2)) error("num and den dimensions must match") @@ -96,15 +96,15 @@ function tf(num::AbstractArray{PT}, den::AbstractArray{PT}, Ts::TimeT) where {T matrix[o, i] = SisoRational{T}(num[o, i], den[o, i]) end end - return TransferFunction{TimeT,SisoRational{T}}(matrix, Ts) + return TransferFunction{TimeT,SisoRational{T}}(matrix, time) end tf(num::AbstractArray{PT}, den::AbstractArray{PT}, Ts::Number) where {T<:Number, PT <: Polynomials.Poly{T}} = tf(num, den, Discrete(Ts)) tf(num::AbstractArray{PT}, den::AbstractArray{PT}) where {T<:Number, PT <: Polynomials.Poly{T}} = tf(num, den, Continuous()) -function tf(num::PT, den::PT, Ts::TimeT) where {TimeT<:TimeType, T<:Number, PT <: Polynomials.Poly{T}} - tf(fill(num,1,1), fill(den,1,1), Ts) +function tf(num::PT, den::PT, time::TimeT) where {TimeT<:TimeType, T<:Number, PT <: Polynomials.Poly{T}} + tf(fill(num,1,1), fill(den,1,1), time) end tf(num::PT, den::PT, Ts::Number) where {T<:Number, PT <: Polynomials.Poly{T}} = tf(num, den, Discrete(Number)) diff --git a/src/types/zpk.jl b/src/types/zpk.jl index 80d16d764..89ffd974d 100644 --- a/src/types/zpk.jl +++ b/src/types/zpk.jl @@ -86,4 +86,4 @@ zpk(z, p, k, Ts::Number) = zpk(z, p, k, Discrete(Ts)) zpk(z, p, k) = zpk(z, p, k, Continuous()) # Catch all 1(2) argument versions zpk(gain, Ts::Number; kwargs...) where {T <: Number} = zpk(gain, Discrete(Ts)) -zpk(gain; kwargs...) where {T <: Number} = zpk(gain, Static()) +zpk(gain; kwargs...) where {T <: Number} = zpk(gain, Continuous()) diff --git a/test/runtests.jl b/test/runtests.jl index ddfd5f836..83f3be091 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,7 @@ include("framework.jl") eye_(n) = Matrix{Int64}(I, n, n) my_tests = [ + "test_statespace", "test_transferfunction", "test_delayed_systems", diff --git a/test/test_connections.jl b/test/test_connections.jl index de0e46482..6bd7fe646 100644 --- a/test/test_connections.jl +++ b/test/test_connections.jl @@ -156,11 +156,15 @@ arr4[1] = ss(0); arr4[2] = ss(1); arr4[3] = ss(2) @test [1.0 D_111] == ss([1.0], [0.0 2.0], [3.0], [1.0 4.0], 0.005) # Type and sample time @test [D_111 1.0] isa StateSpace{Discrete{Float64},Float64,Array{Float64,2}} -@test sampletime([D_111 1.0]) == 0.005 +@test [D_111 1.0].Ts == 0.005 # Continuous version @test [C_111 1.0] == ss([1.0], [2.0 0.0], [3.0], [4.0 1.0]) @test [1.0 C_111] == ss([1.0], [0.0 2.0], [3.0], [1.0 4.0]) @test [C_111 1.0] isa StateSpace{Continuous,Float64,Array{Float64,2}} +@test [C_111 1.0].Ts == 0.0 +@test_logs (:warn, + "Getting time 0.0 for non-discrete systems is deprecated. Check `isdiscrete` before trying to access time." + ) [C_111 1.0].Ts # Concatenation of discrete system with matrix @test [D_222 fill(1.5, 2, 2)] == [D_222 ss(fill(1.5, 2, 2),0.005)] @test [C_222 fill(1.5, 2, 2)] == [C_222 ss(fill(1.5, 2, 2))] diff --git a/test/test_delayed_systems.jl b/test/test_delayed_systems.jl index 2be92062e..2921ba8c5 100644 --- a/test/test_delayed_systems.jl +++ b/test/test_delayed_systems.jl @@ -132,6 +132,9 @@ w = 10 .^ (-2:0.1:2) @test freqresp(s11, w) ≈ freqresp(f2[1,1], w) rtol=1e-15 +@test propertynames(delay(1.0)) == (:P, :Tau) + + #FIXME: A lot more tests, including MIMO systems in particular # Test step diff --git a/test/test_partitioned_statespace.jl b/test/test_partitioned_statespace.jl index 03dd68a1c..ee771ae7d 100644 --- a/test/test_partitioned_statespace.jl +++ b/test/test_partitioned_statespace.jl @@ -38,3 +38,5 @@ sys2 = ControlSystems.PartionedStateSpace(ss(fill(1.0, 2, 2), fill(2.0, 2, 5), f # TODO: Add some tests for interconnections, implicitly tested through delay system implementations though @test (sys1 + sys1).P[1, 1] == (sys1.P[1,1] + sys1.P[1,1]) + +@test propertynames(sys1) == (:P, :nu1, :ny1) diff --git a/test/test_statespace.jl b/test/test_statespace.jl index 369e8da1e..584882892 100644 --- a/test/test_statespace.jl +++ b/test/test_statespace.jl @@ -98,6 +98,12 @@ @test -sys == SS(A, B, -C, -D) + # Accessing Ts through .Ts + @test D_111.Ts == 0.005 + + # propertynames + propertynames(C_111) = (:A, :B, :C, :D, :time, :nx, :nu, :ny, :Ts) + propertynames(D_111) = (:A, :B, :C, :D, :time, :nx, :nu, :ny) # Printing if SS <: StateSpace @@ -109,7 +115,7 @@ res = ("StateSpace{Discrete{Float64},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") end @test sprint(show, D_222) == res - res = ("StateSpace{Static,Float64,Array{Float64,2}}\nD = \n 4.0 0.0\n 0.0 4.0\n\nStatic gain state-space model") + res = ("StateSpace{Continuous,Float64,Array{Float64,2}}\nD = \n 4.0 0.0\n 0.0 4.0\n\nContinuous-time state-space model") @test sprint(show, C_022) == res res = "StateSpace{Discrete{Float64},Float64,Array{Float64,2}}\nD = \n 4.0 0.0\n 0.0 4.0\n\nSample Time: 0.005 (seconds)\nDiscrete-time state-space model" @test sprint(show, D_022) == res diff --git a/test/test_transferfunction.jl b/test/test_transferfunction.jl index 8741bfb7c..84d103198 100644 --- a/test/test_transferfunction.jl +++ b/test/test_transferfunction.jl @@ -90,6 +90,14 @@ 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) +# Accessing Ts through .Ts +@test D_111.Ts == 0.005 + +# propertynames +@test propertynames(C_111) == (:matrix, :time, :nu, :ny) +@test propertynames(D_111) == (:matrix, :time, :nu, :ny, :Ts) + + # Errors @test_throws ErrorException tf(vecarray(1, 1, [1,7,13,15]), vecarray(2, 1, [1,10,31,30], [1,10,31,30]))