Merge pull request #748 from JuliaControl/lsimerror

Add error hint to lsim

Author: baggepinnen

Diff for: .github/workflows/test.yml

@@ -19,7 +19,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        version: ['1.7']
+        version: ['1.8.1']
         os: [ubuntu-latest]
         arch: [x64]
         group:
@@ -47,8 +47,8 @@ jobs:
           GROUP: ${{ matrix.group }}
         # continue-on-error: ${{ matrix.version == 'nightly' }} # Allow nightly to fail and workflow still count as completed
       - uses: julia-actions/julia-processcoverage@v1
-        if: ${{ matrix.version == '1.7' }}
+        if: ${{ matrix.version == '1.8.1' }}
       - uses: codecov/codecov-action@v1
-        if: ${{ matrix.version == '1.7' }}
+        if: ${{ matrix.version == '1.8.1' }}
         with:
           file: lcov.info

Diff for: docs/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 ControlSystems = "a6e380b2-a6ca-5380-bf3e-84a91bcd477e"
+ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
@@ -8,4 +9,4 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-Documenter = "0.27.10"
+Documenter = "0.27.10"

Diff for: docs/src/examples/example.md

@@ -128,7 +128,7 @@ P = tf(B,A)
 
 # output
 
-TransferFunction{Continuous, ControlSystems.SisoRational{Float64}}
+TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
         1.0
 -------------------
 1.0s^2 + 0.4s + 1.0

Diff for: docs/src/lib/nonlinear.md

@@ -42,16 +42,16 @@ plot(step([G; Gu], 5), lab = ["Linear y" "Linear u"])
 plot!(step([Gnl; Gunl], 5), lab = ["Nonlinear y" "Nonlinear u"])
 ```
 
-Since the saturating nonlinearity is common, we provide the constructor [`ControlSystems.saturation`](@ref) that automatically forms the equivalent to `nonlinearity(x->clamp(x, -0.7, 0.7))` while at the same time making sure the function has a recognizable name when the system is printed
+Since the saturating nonlinearity is common, we provide the constructor [`ControlSystemsBase.saturation`](@ref) that automatically forms the equivalent to `nonlinearity(x->clamp(x, -0.7, 0.7))` while at the same time making sure the function has a recognizable name when the system is printed
 ```@example nonlinear
 using ControlSystems: saturation
 saturation(0.7)
 ```
 
-See also [`ControlSystems.ratelimit`](@ref) that saturates the derivative of a signal.
+See also [`ControlSystemsBase.ratelimit`](@ref) that saturates the derivative of a signal.
 
 ### Non-zero operating point
-It's common to linearize nonlinear systems around some operating point. We may make use of the helper constructor [`ControlSystems.offset`](@ref) to create affine functions at the inputs and outputs of the linearized system to, e.g.,
+It's common to linearize nonlinear systems around some operating point. We may make use of the helper constructor [`ControlSystemsBase.offset`](@ref) to create affine functions at the inputs and outputs of the linearized system to, e.g.,
 1. Make sure that simulations result are given in the original coordinates rather than in the coordinates of the linearization point.
 2. Allow nonlinearities that are added back after the linearization (such as saturations) to operate with their original parameters.
 
@@ -172,7 +172,7 @@ plot(step(duffing, 20), title="Duffing oscillator open-loop step response")
 
 We now show how we can make use of the circle criterion to prove stability of the closed loop. The function `circle_criterion` below plots the Nyquist curve of the loop-transfer function and figures out the circle to avoid by finding sector bounds for the static nonlinearity ``f(x) = x^3``. We then choose a controller an check that it stays outside of the circle. To find the sector bounds, we choose a domain to evaluate the nonlinearity over. The function ``f(x) = x^3`` goes to infinity faster than any linear function, and the upper sector bound is thus ∞, but if we restrict the nonlinearity to a smaller domain, we get a finite sector bound:
 ```@example DUFFING
-function circle_criterion(L::ControlSystems.HammersteinWienerSystem, domain::Tuple; N=10000)
+function circle_criterion(L::ControlSystemsBase.HammersteinWienerSystem, domain::Tuple; N=10000)
     fun = x->L.f[](x)/x
     x = range(domain[1], stop=domain[2], length=N)
     0 ∈ x && (x = filter(!=(0), x)) # We cannot divide by zero
@@ -244,9 +244,9 @@ plot(f1,f2, size=(1300,800))
 ## Docstrings
 
 ```@docs
-ControlSystems.nonlinearity
-ControlSystems.offset
-ControlSystems.saturation
-ControlSystems.ratelimit
-ControlSystems.deadzone
+ControlSystemsBase.nonlinearity
+ControlSystemsBase.offset
+ControlSystemsBase.saturation
+ControlSystemsBase.ratelimit
+ControlSystemsBase.deadzone
 ```

Diff for: docs/src/lib/plotting.md

@@ -6,7 +6,7 @@ Pages = ["plotting.md"]
     All plotting requires the user to manually load the Plots.jl library, e.g., by calling `using Plots`.
 
 !!! note "Time-domain responses"
-    There are no special functions to plot time-domain results, such as step and impulse responses, instead, simply call `plot` on the result structure ([`ControlSystems.SimResult`](@ref)) returned by [`lsim`](@ref), [`step`](@ref), [`impulse`](@ref) etc.
+    There are no special functions to plot time-domain results, such as step and impulse responses, instead, simply call `plot` on the result structure ([`ControlSystemsBase.SimResult`](@ref)) returned by [`lsim`](@ref), [`step`](@ref), [`impulse`](@ref) etc.
 
 # Plotting functions
 

Diff for: docs/src/man/creating_systems.md

@@ -19,7 +19,7 @@ tf([1.0],[1,2,1])
 
 # output
 
-TransferFunction{Continuous, ControlSystems.SisoRational{Float64}}
+TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
         1.0
 -------------------
 1.0s^2 + 2.0s + 1.0
@@ -42,7 +42,7 @@ zpk([-1.0,1], [-5, -10], 2)
 
 # output
 
-TransferFunction{Continuous, ControlSystems.SisoZpk{Float64, Float64}}
+TransferFunction{Continuous, ControlSystemsBase.SisoZpk{Float64, Float64}}
    (1.0s + 1.0)(1.0s - 1.0)
 2.0-------------------------
    (1.0s + 5.0)(1.0s + 10.0)
@@ -71,7 +71,7 @@ function HeteroStateSpace(A,B,C,D,Ts=0,f::F=to_static) where F
     HeteroStateSpace(f(A),f(B),f(C),f(D),Ts)
 end
 HeteroStateSpace(s,f) = HeteroStateSpace(s.A,s.B,s.C,s.D,s.timeevol,f)
-ControlSystemsBase._string_mat_with_headers(a::SizedArray) = ControlSystems._string_mat_with_headers(Matrix(a)); # Overload for printing purposes
+ControlSystemsBase._string_mat_with_headers(a::SizedArray) = ControlSystemsBase._string_mat_with_headers(Matrix(a)); # Overload for printing purposes
 
 nothing # hide
 ```
@@ -90,7 +90,7 @@ tf(zpk([-1], [1], 2, 0.1))
 
 # output
 
-TransferFunction{Discrete{Float64}, ControlSystems.SisoRational{Int64}}
+TransferFunction{Discrete{Float64}, ControlSystemsBase.SisoRational{Int64}}
 2z + 2
 ------
 z - 1
@@ -221,4 +221,4 @@ StateSpace[P_cont, P_disc]
 The type `StateSpace` is abstract, since the type parameters are not specified.
 
 ## Demo systems
-The module `ControlSystems.DemoSystems` contains a number of demo systems demonstrating different kinds of dynamics.
+The module `ControlSystemsBase.DemoSystems` contains a number of demo systems demonstrating different kinds of dynamics.

Diff for: docs/src/man/differences.md

@@ -9,7 +9,7 @@ The rest of this page will list noteworthy differences between ControlSystems.jl
 - Simulation using [`lsim`](@ref), [`step`](@ref), [`impulse`](@ref) etc. return a structure that can be plotted. These functions never plot anything themselves.
 - Functions [`bode`](@ref), [`nyquist`](@ref) etc. never produce a plot. Instead, see [`bodeplot`](@ref), [`nyquistplot`](@ref) etc.
 - In Julia, functionality is often split up into several different packages. You may therefore have to install and use additional packages in order to cover all your needs. See [Ecosystem](@ref) for a collection of control-related packages.
-- In Julia, `1` has a different type than `1.0`, and the types in ControlSystems.jl respect the types chosen by the user. As an example, `tf(1, [1, 1])` is a transfer function with integer coefficients, while `tf(1.0, [1, 1])` will promote all coefficients to `Float64`.
+- In Julia, `1` has a different type than `1.0`, and the types in ControlSystemsBase.jl respect the types chosen by the user. As an example, `tf(1, [1, 1])` is a transfer function with integer coefficients, while `tf(1.0, [1, 1])` will promote all coefficients to `Float64`.
 - In Julia, code can often be differentiated using automatic differentiation. When using ControlSystems.jl, we recommend trying [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl/) for AD. An example making use of this is available [here](https://github.com/JuliaControl/ControlExamples.jl/blob/master/autotuning.ipynb).
 - In Julia, the source code is often very readable. If you want to learn how a function is implemented, you may find the macros `@edit` or `@less` useful to locate the source code.
 - If you run into an issue (bug) with a Julia package, you can share this issue (bug report) on the package's github page and it will often be fixed promptly. To open an issue with ControlSystems.jl, [click here](https://github.com/JuliaControl/ControlSystems.jl/issues/new/choose). Thank you for helping out improving open-source software!

Diff for: docs/src/man/introduction.md

@@ -32,7 +32,7 @@ T = P/(1+P)
 
 # output
 
-TransferFunction{Continuous, ControlSystems.SisoRational{Float64}}
+TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
     1.0s + 1.0
 -------------------
 1.0s^2 + 3.0s + 2.0
@@ -46,7 +46,7 @@ minreal(T)
 
 # output
 
-TransferFunction{Continuous, ControlSystems.SisoRational{Float64}}
+TransferFunction{Continuous, ControlSystemsBase.SisoRational{Float64}}
    1.0
 ----------
 1.0s + 2.0

Diff for: docs/src/man/numerical.md

@@ -47,7 +47,7 @@ w = exp10.(LinRange(-2, 2, 20000))
 # 55.120 ms (517957 allocations: 24.42 MiB)
 @btime bode($G, $w, unwrap=false); # phase unwrapping is slow
 # 3.624 ms (7 allocations: 2.44 MiB)
-ws = ControlSystems.BodemagWorkspace(G, w)
+ws = ControlSystemsBase.BodemagWorkspace(G, w)
 @btime bodemag!($ws, $G, $w);
 # 2.991 ms (1 allocation: 64 bytes)
 ```

Diff for: lib/ControlSystemsBase/src/ControlSystemsBase.jl

@@ -212,4 +212,14 @@ end
 # The path has to be evaluated upon initial import
 const __ControlSystemsBase_SOURCE_DIR__ = dirname(Base.source_path())
 
+function __init__()
+    Base.Experimental.register_error_hint(MethodError) do io, exc, argtypes, kwargs
+        if exc.f ∈ (lsim, step, impulse) && argtypes[1] <: LTISystem{Continuous}
+            print(io, "\n$(exc.f) with continuous-time systems, including delay systems and nonlinear systems, require the user to first ")
+            printstyled(io, "install and load ControlSystems.jl, or pass the keyword method = :zoh", color=:green, bold=true)
+            print(io, " for automatic discretization (applicable to systems without delays or nonlinearities only).")
+        end
+    end
+end
+
 end

Diff for: lib/ControlSystemsBase/src/delay_systems.jl

@@ -44,7 +44,7 @@ i.e. of z^-N where N = τ / Ts.
 τ must be a multiple of Ts.
 """
 function delayd_ss(τ::Number, Ts::Number)
-    n = Int(round(τ / Ts))
+    n = round(Int, τ / Ts)
     if !(τ - n*Ts ≈ 0)
         error("The delay τ must be a multiple of the sample time Ts")
     end

Diff for: lib/ControlSystemsBase/src/timeresp.jl

@@ -41,8 +41,9 @@ Calculate the step response of system `sys`. If the final time `tfinal` or time
 vector `t` is not provided, one is calculated based on the system pole
 locations. The return value is a structure of type `SimResult` that can be plotted or destructured as `y, t, x = result`.
 
-`y` has size `(ny, length(t), nu)`, `x` has size `(nx, length(t), nu)`"""
-function Base.step(sys::AbstractStateSpace, t::AbstractVector; method=:cont, kwargs...)
+`y` has size `(ny, length(t), nu)`, `x` has size `(nx, length(t), nu)`
+"""
+function Base.step(sys::AbstractStateSpace, t::AbstractVector; method=:zoh, kwargs...)
     T = promote_type(eltype(sys.A), Float64)
     ny, nu = size(sys)
     nx = nstates(sys)
@@ -122,7 +123,7 @@ plot(result, plotu=true, plotx=false)
 ```
 `y`, `x`, `u` have time in the second dimension. Initial state `x0` defaults to zero.
 
-Continuous-time systems are simulated using an ODE solver if `u` is a function. If `u` is an array, the system is discretized (with `method=:zoh` by default) before simulation. For a lower-level inteface, see `?Simulator` and `?solve`
+Continuous-time systems are simulated using an ODE solver if `u` is a function (requires using ControlSystems). If `u` is an array, the system is discretized (with `method=:zoh` by default) 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(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`.
@@ -132,6 +133,7 @@ For maximum performance, see function [`lsim!`](@ref), avaialable for discrete-t
 
 Usage example:
 ```julia
+using ControlSystems
 using LinearAlgebra # For identity matrix I
 using Plots
 
@@ -253,7 +255,7 @@ function lsim(sys::AbstractStateSpace, u::Function, t::AbstractVector;
         end
         x,uout = ltitr(simsys.A, simsys.B, u, t, T.(x0))
     else
-        error("Continuous-time simulation requires using ControlSystems")
+        throw(MethodError(lsim, (sys, u, t)))
     end
     y = sys.C*x
     if !iszero(sys.D)

Diff for: lib/ControlSystemsBase/test/test_analysis.jl

@@ -139,7 +139,6 @@ sol_k = [0.0 3.0; 1.0 2.0]
 z, p, k = zpkdata(H)
 
 @test_array_vecs_eps z sol_z 2*eps(Float64)
-@test_broken true == false # order of poles changed below, should probably be consistent
 #@test_array_vecs_eps p sol_p 2*eps(Float64)
 @test k == sol_k
 z, p, k = zpkdata(G)
@@ -183,11 +182,11 @@ sys = s*(s + 1)*(s^2 + 1)*(s - 3)/((s + 1)*(s + 4)*(s - 4))
 
 @test sprint(dampreport, 1/(s+1+2im)/(s+2+3im)) == "|        Pole        |   Damping     |   Frequency   |   Frequency   | Time Constant |\n|                    |    Ratio      |   (rad/sec)   |     (Hz)      |     (sec)     |\n+--------------------+---------------+---------------+---------------+---------------+\n| -1            -2im |  0.447        |  2.24         |  0.356        |  1            |\n| -2            -3im |  0.555        |  3.61         |  0.574        |  0.5          |\n"
 
-# Example 5.5 from http://www.control.lth.se/media/Education/EngineeringProgram/FRTN10/2017/e05_both.pdf
+# Example 5.5 from https://www.control.lth.se/fileadmin/control/Education/EngineeringProgram/FRTN10/2019/e05_both.pdf
 G = [1/(s+2) -1/(s+2); 1/(s+2) (s+1)/(s+2)]
 @test_broken length(poles(G)) == 1
 @test length(tzeros(G)) == 1
-@test_broken size(minreal(ss(G)).A) == (2,2)
+@test minreal(ss(G)).A ≈ [-2]
 
 
 ## MARGIN ##

Diff for: lib/ControlSystemsBase/test/test_timeresp.jl

@@ -13,7 +13,7 @@ L = lqr(sys,Q,R)
 u1(x,i) = -L*x # Form control law
 t=0:0.1:50
 x0 = [1.,0]
-@test_throws ErrorException lsim(sys,u1,t,x0=x0) # Continuous time simulation not loaded
+@test_throws MethodError lsim(sys,u1,t,x0=x0) # Continuous time simulation not loaded
 
 
 th = 1e-6
@@ -194,4 +194,14 @@ sysd = tf(1, [1, 1], 0.01)
 @test ControlSystemsBase._default_dt(sysint) == 0.05
 @test ControlSystemsBase._default_dt(sysd) == 0.01
 
+
+# Test error hints
+if VERSION >= v"1.7"
+    # If ControlSystems is not loaded, the 
+    G = ssrand(1,1,1)
+    @test_throws "install and load ControlSystems.jl" lsim(G, (u,t)->[1], 10)
+    @test_throws "install and load ControlSystems.jl" lsim(G*delay(1), (u,t)->[1], 10)
+    @test_throws "step with continuous-time systems" step(G*delay(1), 10)
+end
+
 end

Diff for: test/runtests.jl

@@ -19,6 +19,8 @@ dev_subpkg("ControlSystemsBase") # Always dev this package to test with the late
 
 const GROUP = get(ENV, "GROUP", "All") # Get the GROUP attribute from the test.yml file, default to "All" for testing locally
 
+@show GROUP
+
 if GROUP ∈ ("ControlSystems", "All")
     include("runtests_controlsystems.jl")
 end
@@ -28,7 +30,7 @@ if GROUP == "All"
         subpkg_path = joinpath(dirname(@__DIR__), "lib", GROUP)
         Pkg.test(PackageSpec(name = GROUP, path = subpkg_path))
     end
-else
+elseif GROUP != "ControlSystems"
     # dev_subpkg(GROUP) # Do this if more sub packages are added, don't forget to avoid doing it if GROUP is CSBase
     subpkg_path = joinpath(dirname(@__DIR__), "lib", GROUP)
     Pkg.test(PackageSpec(name = GROUP, path = subpkg_path))

Diff for: test/test_timeresp.jl

@@ -96,7 +96,7 @@ xreal[3,:,2] = exp.(-t).*t
     x0             = [0.,0]
     tspan          = (0.0,tfinal)
     sol            = solve(s, x0, tspan, OrdinaryDiffEq.Tsit5())
-    @test step(P,t)[3] ≈ reduce(hcat,sol.(t))
+    @test step(P,t)[3] ≈ reduce(hcat,sol.(t)) rtol=1e-4
 end
 
 # test fwd euler