set abs and rel tols in simulate...

Project.toml

@@ -1,7 +1,7 @@
 name = "Controlz"
 uuid = "e034abe6-471a-4d54-96dd-ecd1f4022419"
 authors = ["SimonEnsemble <[email protected]>"]
-version = "0.3.5"
+version = "0.3.6"
 
 [deps]
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
@@ -15,13 +15,13 @@ ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 
 [compat]
-DifferentialEquations = "^7.12.0"
-DataFrames = "^1.6.1"
-Interpolations = "^0.14.7"
-Polynomials = "^4.0.6"
+DifferentialEquations = "^7.13.0"
+DataFrames = "^1.7"
+Interpolations = "^0.15.1"
+Polynomials = "^4.0.9"
 julia = "^1.10"
-Roots = "^2.1.2"
-CairoMakie = "^0.11.8"
+Roots = "^2.1.6"
+CairoMakie = "^0.12.2"
 ColorSchemes = "^3.24.0"
 Colors = "^0.12.10"
 

src/closed_loops.jl

@@ -130,7 +130,7 @@ function tf_to_ss(cl::CLTFStandard)
 end
 
 # see sim.jl for doc string
-function simulate(cl::ClosedLoopTransferFunction, final_time::Float64; nb_time_points::Int=100)
+function simulate(cl::ClosedLoopTransferFunction, final_time::Float64; nb_time_points::Int=250)
 	cls = CLTFStandard(cl)
 
     if ! strictly_proper(cls)
@@ -145,7 +145,7 @@ function simulate(cl::ClosedLoopTransferFunction, final_time::Float64; nb_time_p
 	h(p, t) = zeros(order(cls))
     f(x, h, p, t) = A * x - C * h(p, t - cls.ϕ) # RHS of ODE (ignore p for params)
     prob = DDEProblem(f, x0, h, (0.0, final_time))
-    sol = solve(prob, d_discontinuities=[0.0, cls.ϕ, cls.θ, cls.θ + cls.ϕ])
+    sol = solve(prob, d_discontinuities=[0.0, cls.ϕ, cls.θ, cls.θ + cls.ϕ], abstol=1e-8, reltol=1e-8)
 
 	t = vcat([-0.05 * final_time, -1e-5], 
 	range(1e-5, final_time, length=nb_time_points - 2))

src/sim.jl

@@ -50,7 +50,7 @@ end
 
 
 @doc raw"""
-    data = simulate(Y, final_time, nb_time_points=100) # invert Y(s)
+    data = simulate(Y, final_time, nb_time_points=250) # invert Y(s)
 
 simulate the output $y(t)$ of an LTI system, given the Laplace transform of the output, $Y(s)$, `Y`.
 
@@ -89,7 +89,7 @@ first(data, 5)           # show the first 5 rows of the data frame
    5 │  0.247432  0.316671
 ```
 """
-function simulate(Y::TransferFunction, final_time::Union{Float64, Int64}; nb_time_points::Int=100)
+function simulate(Y::TransferFunction, final_time::Union{Float64, Int64}; nb_time_points::Int=250)
     if ! proper(Y)
         error("LTI system is not proper...")
     end
@@ -105,7 +105,7 @@ function simulate(Y::TransferFunction, final_time::Union{Float64, Int64}; nb_tim
 
     f(x, p, t) = A * x # RHS of ODE (ignore p for params)
     prob = ODEProblem(f, x0, (0.0, 1.0 * final_time))
-    sol = solve(prob, d_discontinuities=[0.0])
+    sol = solve(prob, Tsit5(), d_discontinuities=[0.0], saveat=final_time / nb_time_points, reltol=1e-8, abstol=1e-8)
 
     t = vcat([-0.05 * final_time, -1e-5], range(1e-5, final_time, length=nb_time_points - 2))
     y = [NaN for i = 1:nb_time_points]

test/runtests.jl

@@ -279,9 +279,9 @@ end
     # L[cos(at)] = s/(s^2+a^2)
     a = 2.3
     U = s / (s^2+a^2)
-    u_data = simulate(U, 12.0)
+    u_data = simulate(U, 12.0, nb_time_points=500)
     _cosat(t::Float64) = t < 0.0 ? 0.0 : cos(a*t)
-    @test isapprox(u_data[:, :output], _cosat.(u_data[:, :t]), rtol=0.001)
+    @test isapprox(u_data[:, :output], _cosat.(u_data[:, :t]), rtol=0.005)
     # L[shifted step] = e^{-a*s}/s
     U = exp(-a*s) / s
     u_data = simulate(U, 12.0)
@@ -313,7 +313,7 @@ end
     data = simulate(g, 12.0)
     y_truth = K / τ * exp.(-data[:, :t]/τ)
     y_truth[data[:, :t] .< 0.0] .= 0.0
-    @test isapprox(y_truth, data[:, :output], rtol=0.0001)
+    @test isapprox(y_truth, data[:, :output], rtol=0.005)
 
     # first order ramp input
     a = 2.0 # slope of ramp
@@ -332,7 +332,7 @@ end
     y_truth = τ * ω / (1 + (τ*ω)^2) * exp.(-data[:, :t] ./ τ) .+ 1 / sqrt(1+(τ*ω)^2) * sin.(ω*data[:, :t] .+ ϕ)
     y_truth *= K * A
     y_truth[data[:, :t] .< 0.0] .= 0.0
-    @test isapprox(y_truth, data[:, :output], rtol=0.001)
+    @test isapprox(y_truth, data[:, :output], rtol=0.005)
 
     # FOPTD
     M = 3.3
@@ -341,10 +341,10 @@ end
     K = 9.0
     g = K / (τ * s + 1) * exp(-θ * s)
     u = M / s
-    data = simulate(g * u, 10.0)
+    data = simulate(g * u, 10.0, nb_time_points=1000)
     y_truth = K * M * (1.0 .- exp.(-(data[:, :t] .- θ) ./ τ))
     y_truth[data[:, :t] .< θ] .= 0.0
-    @test isapprox(y_truth, data[:, :output], atol=0.001)
+    @test isapprox(y_truth, data[:, :output], atol=0.005)
 
     ##
     # second order, overdamped
@@ -401,8 +401,8 @@ end
     @test isapprox(interpolate(data, 1.2), 1.2*5)
     τ = 3.45
     g = 1 / (τ * s + 1)
-    data = simulate(g / s, 10.0)
-    @test isapprox(interpolate(data, τ), 0.632, atol=0.002)
+    data = simulate(g / s, 10.0, nb_time_points=300)
+    @test isapprox(interpolate(data, τ), 1.0 - exp(-1.0), atol=0.002)
 end
 
 @testset "testing Controls" begin
@@ -467,9 +467,9 @@ end
     # run without a time delay to compare to other simulate function.
     gp = 3 / (s + 2)
     gc = TransferFunction(PIController(1.0, 3.0))
-    data_old = simulate(gp * gc / (1 + gp * gc), 10.0)
-    data_new = simulate(ClosedLoopTransferFunction(gp * gc, gp * gc), 10.0)
-    @test isapprox(data_old[:, :output], data_new[:, :output], atol=0.0001)
+    data_old = simulate(gp * gc / (1 + gp * gc), 10.0, nb_time_points=1000)
+    data_new = simulate(ClosedLoopTransferFunction(gp * gc, gp * gc), 10.0, nb_time_points=1000)
+    @test isapprox(data_old[:, :output], data_new[:, :output], atol=0.005)
 end
 
 @testset "example notebook (also to generate images for docs)" begin