diff --git a/Project.toml b/Project.toml index a199654..62b884c 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Controlz" uuid = "e034abe6-471a-4d54-96dd-ecd1f4022419" authors = ["SimonEnsemble "] -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" diff --git a/src/closed_loops.jl b/src/closed_loops.jl index 96e192b..460f3d9 100644 --- a/src/closed_loops.jl +++ b/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)) diff --git a/src/sim.jl b/src/sim.jl index 612acfd..8eb8190 100644 --- a/src/sim.jl +++ b/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] diff --git a/test/runtests.jl b/test/runtests.jl index 86a0391..30d7dba 100644 --- a/test/runtests.jl +++ b/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