From d6818e4c1fa777740011d0ada3edafbe38c76f27 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Tue, 3 Sep 2024 13:15:29 +0200 Subject: [PATCH 1/5] add stable noise sources and tutorial --- Project.toml | 10 ++ docs/Project.toml | 1 + docs/make.jl | 4 +- docs/src/{ => tutorials}/SampledData.md | 22 +++-- docs/src/tutorials/noise.md | 122 ++++++++++++++++++++++++ src/ModelingToolkitSampledData.jl | 1 + src/discrete_blocks.jl | 100 ++++++++++++++++--- test/test_discrete_blocks.jl | 33 ++++++- 8 files changed, 270 insertions(+), 23 deletions(-) rename docs/src/{ => tutorials}/SampledData.md (93%) create mode 100644 docs/src/tutorials/noise.md diff --git a/Project.toml b/Project.toml index 0f7997b..7c72009 100644 --- a/Project.toml +++ b/Project.toml @@ -4,14 +4,24 @@ authors = ["Fredrik Bagge Carlson"] version = "0.1.0" [deps] +DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" JuliaSimCompiler = "8391cb6b-4921-5777-4e45-fd9aab8cb88d" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" ModelingToolkitStandardLibrary = "16a59e39-deab-5bd0-87e4-056b12336739" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" [compat] +DiffEqCallbacks = "~3.8" julia = "1.10" +StableRNGs = "1" +JuliaSimCompiler = "1.19" +ModelingToolkit = "9" +ModelingToolkitStandardLibrary = "2" +OrdinaryDiffEq = "6.89" +Random = "1" + [extras] ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" diff --git a/docs/Project.toml b/docs/Project.toml index 781c10b..c762060 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +ControlSystemIdentification = "3abffc1c-5106-53b7-b354-a47bfc086282" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" JuliaSimCompiler = "8391cb6b-4921-5777-4e45-fd9aab8cb88d" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" diff --git a/docs/make.jl b/docs/make.jl index 33aa31b..848ac61 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -26,12 +26,14 @@ makedocs(; pages = [ "Home" => "index.md", "Tutorials" => [ - "Getting started with Sampled-Data Systems" => "SampledData.md", + "Getting started with Sampled-Data Systems" => "tutorials/SampledData.md", + "Noise" => "tutorials/noise.md", ], "Examples" => [ "Controlled DC motor" => "examples/dc_motor_pi.md", ], "Block library" => "blocks.md", + "API" => "api.md", ], ) diff --git a/docs/src/SampledData.md b/docs/src/tutorials/SampledData.md similarity index 93% rename from docs/src/SampledData.md rename to docs/src/tutorials/SampledData.md index ec38dd0..12240df 100644 --- a/docs/src/SampledData.md +++ b/docs/src/tutorials/SampledData.md @@ -14,17 +14,22 @@ A sampled-data system contains both continuous-time and discrete-time components ## Clocks, operators and difference equations A clock can be seen as an *event source*, i.e., when the clock ticks, an event is generated. In response to the event the discrete-time logic is executed, for example, a control signal is computed. For basic modeling of sampled-data systems, the user does not have to interact with clocks explicitly, instead, the modeling is performed using the operators - - [`Sample`](@ref) - - [`Hold`](@ref) - - [`ShiftIndex`](@ref) +- [`ModelingToolkit.Sample`](@ref) +- [`ModelingToolkit.Hold`](@ref) +- [`ModelingToolkit.ShiftIndex`](@ref) + + or their corresponding components + +- [`ModelingToolkitSampledData.Sampler`](@ref) +- [`ModelingToolkitSampledData.ZeroOrderHold`](@ref) When a continuous-time variable `x` is sampled using `xd = Sample(x, dt)`, the result is a discrete-time variable `xd` that is defined and updated whenever the clock ticks. `xd` is *only defined when the clock ticks*, which it does with an interval of `dt`. If `dt` is unspecified, the tick rate of the clock associated with `xd` is inferred from the context in which `xd` appears. Any variable taking part in the same equation as `xd` is inferred to belong to the same *discrete partition* as `xd`, i.e., belonging to the same clock. A system may contain multiple different discrete-time partitions, each with a unique clock. This allows for modeling of multi-rate systems and discrete-time processes located on different computers etc. -To make a discrete-time variable available to the continuous partition, the [`Hold`](@ref) operator is used. `xc = Hold(xd)` creates a continuous-time variable `xc` that is updated whenever the clock associated with `xd` ticks, and holds its value constant between ticks. +To make a discrete-time variable available to the continuous partition, the `Hold` operator is used. `xc = Hold(xd)` creates a continuous-time variable `xc` that is updated whenever the clock associated with `xd` ticks, and holds its value constant between ticks. -The operators [`Sample`](@ref) and [`Hold`](@ref) are thus providing the interface between continuous and discrete partitions. +The operators `Sample` and `Hold` are thus providing the interface between continuous and discrete partitions. -The [`ShiftIndex`](@ref) operator is used to refer to past and future values of discrete-time variables. The example below illustrates its use, implementing the discrete-time system +The `ShiftIndex` operator is used to refer to past and future values of discrete-time variables. The example below illustrates its use, implementing the discrete-time system ```math x(k+1) = 0.5x(k) + u(k) @@ -48,7 +53,7 @@ eqs = [ A few things to note in this basic example: - The equation `x(k+1) = 0.5x(k) + u(k)` has been rewritten in terms of negative shifts since positive shifts are not yet supported. - - `x` and `u` are automatically inferred to be discrete-time variables, since they appear in an equation with a discrete-time [`ShiftIndex`](@ref) `k`. + - `x` and `u` are automatically inferred to be discrete-time variables, since they appear in an equation with a discrete-time `ShiftIndex` `k`. - `y` is also automatically inferred to be a discrete-time-time variable, since it appears in an equation with another discrete-time variable `x`. `x,u,y` all belong to the same discrete-time partition, i.e., they are all updated at the same *instantaneous point in time* at which the clock ticks. - The equation `y ~ x` does not use any shift index, this is equivalent to `y(k) ~ x(k)`, i.e., discrete-time variables without shift index are assumed to refer to the variable at the current time step. - The equation `x(k) ~ 0.5x(k-1) + u(k-1)` indicates how `x` is updated, i.e., what the value of `x` will be at the *current* time step in terms of the *past* value. The output `y`, is given by the value of `x` at the *current* time step, i.e., `y(k) ~ x(k)`. If this logic was implemented in an imperative programming style, the logic would thus be @@ -106,11 +111,10 @@ eqs = [ ] ``` -(see also [ModelingToolkitStandardLibrary](https://docs.sciml.ai/ModelingToolkitStandardLibrary/stable/) for a discrete-time transfer-function component.) ## Initial conditions -The initial condition of discrete-time variables is defined using the [`ShiftIndex`](@ref) operator, for example +The initial condition of discrete-time variables is defined using the `ShiftIndex` operator, for example ```julia ODEProblem(model, [x(k-1) => 1.0], (0.0, 10.0)) diff --git a/docs/src/tutorials/noise.md b/docs/src/tutorials/noise.md new file mode 100644 index 0000000..b4d9f37 --- /dev/null +++ b/docs/src/tutorials/noise.md @@ -0,0 +1,122 @@ +# Measurement noise and corruption +Measurement noise is practically always present in signals originating from real-world sensors. In a sampled-data system, analyzing the influence of measurement noise using simulation is relatively straight forward. Below, we add Gaussian white noise to the speed sensor signal in the DC motor example. The noise is added using the [`NormalNoise`](@ref) block. + +This block has two modes of operation +1. If `additive = false` (default), the block has the connector `output` only, and this output is the noise signal. +2. If `additive = true`, the block has the connectors `input` and `output`, and the output is the sum of the input and the noise signal, i.e., the noise is _added_ to the input signal. This mode makes it convenient to add noise to a signal in a sampled-data system. + +## Example: Noise +```@example NOISE +using ModelingToolkit +using ModelingToolkit: t_nounits as t +using ModelingToolkitStandardLibrary.Electrical +using ModelingToolkitStandardLibrary.Mechanical.Rotational +using ModelingToolkitStandardLibrary.Blocks +using ModelingToolkitSampledData +using JuliaSimCompiler +using OrdinaryDiffEq +using Plots + +R = 0.5 # [Ohm] armature resistance +L = 4.5e-3 # [H] armature inductance +k = 0.5 # [N.m/A] motor constant +J = 0.02 # [kg.m²] inertia +f = 0.01 # [N.m.s/rad] friction factor +tau_L_step = -0.3 # [N.m] amplitude of the load torque step +nothing # hide + +z = ShiftIndex() + +@mtkmodel NoisyClosedLoop begin + @components begin + ground = Ground() + source = Voltage() + ref = Blocks.Step(height = 1, start_time = 0, smooth = false) + sampler = Sampler(dt = 0.002) + noise = NormalNoise(sigma = 0.1, additive = true) + pi_controller = DiscretePIDStandard( + K = 1, Ti = 0.035, u_max = 10, with_D = false) + zoh = ZeroOrderHold() + R1 = Resistor(R = R) + L1 = Inductor(L = L) + emf = EMF(k = k) + fixed = Fixed() + load = Torque() + load_step = Blocks.Step(height = tau_L_step, start_time = 1.3) + inertia = Inertia(J = J) + friction = Damper(d = f) + speed_sensor = SpeedSensor() + angle_sensor = AngleSensor() + end + + @equations begin + connect(fixed.flange, emf.support, friction.flange_b) + connect(emf.flange, friction.flange_a, inertia.flange_a) + connect(inertia.flange_b, load.flange) + connect(inertia.flange_b, speed_sensor.flange, angle_sensor.flange) + connect(load_step.output, load.tau) + connect(ref.output, pi_controller.reference) + connect(speed_sensor.w, sampler.input) + connect(sampler.output, noise.input) + connect(noise.output, pi_controller.measurement) + connect(pi_controller.ctr_output, zoh.input) + connect(zoh.output, source.V) + connect(source.p, R1.p) + connect(R1.n, L1.p) + connect(L1.n, emf.p) + connect(emf.n, source.n, ground.g) + end +end + + +@named noisy_model = NoisyClosedLoop() +noisy_model = complete(noisy_model) +ssys = structural_simplify(IRSystem(noisy_model)) # Conversion to an IRSystem from JuliaSimCompiler is required for sampled-data systems + +noise_prob = ODEProblem(ssys, [unknowns(noisy_model) .=> 0.0; noisy_model.pi_controller.I(z-1) => 0; noisy_model.pi_controller.eI(z-1) => 0; noisy_model.noise.y(z-1) => 0], (0, 2.0)) +noise_sol = solve(noise_prob, Tsit5()) + +figy = plot(noise_sol, idxs=noisy_model.noise.y, label = "Measured speed", ) +plot!(noise_sol, idxs=noisy_model.inertia.w, ylabel = "Angular Vel. [rad/s]", + label = "Actual speed", legend=:bottomleft, dpi=600, l=(2, :blue)) +figu = plot(noise_sol, idxs=noisy_model.source.V.u, label = "Control signal [V]", ) +plot(figy, figu, plot_title = "DC Motor with Discrete-time Speed Controller") +``` + +## Linear analysis of noise +Propagation of Gaussian noise through linear time-invariant systems is well understood, the stationary covariance of the output can be computed by solving a Lyapunov equation. Unfortunately, ModelingToolkit models that contain both continuous time and discrete time components cannot yet be linearized and linear analysis is thus made slightly harder. Below, we instead use a data-driven linearization approach where we use recorded signals from the simulation and fit a linear model using subspace-based identification. The function `subspaceid` below is provided by the package [ControlSystemIdentification.jl](https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/). + +We let the angular velocity of the inertia be the output, and the output of the noise block as well as the output of the load disturbance be the inputs. + +```@example NOISE +using ControlSystemIdentification, ControlSystemsBase +Tf = 20 +prob2 = remake(noise_prob, p=Dict(noisy_model.load_step.height=>0.0), tspan=(0.0, Tf)) +noise_sol = solve(prob2, Tsit5()) +tv = 0:0.002:Tf +y = noise_sol(tv, idxs=noisy_model.inertia.w) |> vec +un = noise_sol(tv, idxs=noisy_model.noise.y)-y |> vec +ud = noise_sol(tv, idxs=noisy_model.load_step.output.u) |> vec +d = iddata(y', [un ud]', 0.002) +lsys,_ = newpem(d, 4, focus=:simulation, zeroD=false) +``` +With an LTI model available, we can ask for the theoretical output covariance we should obtain if we feed a white noise signal with covariance matrix ``0.1^2 I`` through the noise input of the system. We compare this to the actual output covariance obtained from the simulation (discarding the initial transient as well as the transient caused by the load disturbance). +```@example NOISE +sqrt(covar(lsys[1,1],0.1^2*I)), std(y[[50:648; 750:end]]) +``` + +## Noise filtering +No discrete-time filter components are available yet. You may, e.g. +- Add exponential filtering using `xf(k) ~ (1-α)xf(k-1) + α*x(k)`, where `α` is the filter coefficient and `x` is the signal to be filtered. +- Add moving average filtering using `xf(k) ~ 1/N sum(i->x(k-i), i=0:N-1)`, where `N` is the number of samples to average over. + +## Colored noise +Colored noise can be achieved by filtering white noise through a filter with the desired spectrum. No components are available for this yet. + +## Internal details +Internally, a random number generator from [StableRNGs.jl](https://github.com/JuliaRandom/StableRNGs.jl) is used to produce reproducible streams of random numbers. Each draw of a random number is seeded by `hash(t, hash(seed))`, where `seed` is a parameter in the noise source component, and `t` is the current simulation time. This ensures that +1. The user can alter the stream of random numbers with `seed`. +2. Multiple calls to the random number generator at the same time step all return the same number. + +## Quantization +Not yet available. \ No newline at end of file diff --git a/src/ModelingToolkitSampledData.jl b/src/ModelingToolkitSampledData.jl index 44cc7c7..58ce736 100644 --- a/src/ModelingToolkitSampledData.jl +++ b/src/ModelingToolkitSampledData.jl @@ -1,6 +1,7 @@ module ModelingToolkitSampledData using ModelingToolkit using JuliaSimCompiler +using StableRNGs export get_clock export DiscreteIntegrator, DiscreteDerivative, Delay, Difference, ZeroOrderHold, Sampler, diff --git a/src/discrete_blocks.jl b/src/discrete_blocks.jl index ff342c6..8d89fc6 100644 --- a/src/discrete_blocks.jl +++ b/src/discrete_blocks.jl @@ -144,35 +144,68 @@ A discrete-time noise source that returns a normally-distributed value at each c # Parameters - `mu = 0`: Mean - `sigma = 1`: Standard deviation +- `seed = 1`: Seed for the random number generator # Structural parameters -- `rng`: A random number generator, defaults to `Random.Xoshiro()`. - `z`: The `ShiftIndex` used to indicate clock partition. # Connectors: - `output` """ @mtkmodel NormalNoise begin + @structural_parameters begin + z = ShiftIndex() + additive = false + end @components begin output = RealOutput() + if additive + input = RealInput() + end end @variables begin y(t), [description = "Output variable"] + if additive + u(t), [description = "Input variable"] + end + n(t), [description = "Internal noise variable"] end @parameters begin mu = 0 sigma = 1 - end - @structural_parameters begin - z = ShiftIndex() - rng = Random.Xoshiro() + seed = 1 end @equations begin - y(z) ~ mu + sigma*Symbolics.term(randn, rng; type=Real) output.u ~ y + n(z) ~ mu + sigma*Symbolics.term(seeded_randn, seed, t; type=Real) + # n(z) ~ mu + sigma*Symbolics.term(randn, rng; type=Real) + if additive + y(z) ~ u(z) + n(z) + 1e-100*y(z-1) # The 0*y(z-1) is a workaround for a bug in the compiler, to force the y variable to be a discrete-time state variable + u ~ input.u + else + y(z) ~ n(z) + 1e-100*y(z-1) + end end end +""" + seeded_randn(seed, t) + +Internal function. This function seeds the seed parameter as well as the current simulation time. +""" +function seeded_randn(seed, t) + rng = StableRNGs.StableRNG(hash(t, hash(seed))) + randn(rng) +end +""" + seeded_rand(seed, t) + +Internal function. This function seeds the seed parameter as well as the current simulation time. +""" +function seeded_rand(seed, t) + rng = StableRNGs.StableRNG(hash(t, hash(seed))) + rand(rng) +end """ UniformNoise() @@ -182,6 +215,7 @@ A discrete-time noise source that returns a uniformly distributed value at each # Parameters - `l = 0`: Lower bound - `u = 1`: Upper bound +- `seed = 1`: Seed for the random number generator # Structural parameters - `rng`: A random number generator, defaults to `Random.Xoshiro()`. @@ -191,26 +225,70 @@ A discrete-time noise source that returns a uniformly distributed value at each - `output` """ @mtkmodel UniformNoise begin + @structural_parameters begin + z = ShiftIndex() + rng = Random.Xoshiro() + additive = false + end @components begin output = RealOutput() + if additive + input = RealInput() + end end @variables begin y(t), [description = "Output variable"] + n(t), [description = "Internal noise variable"] end @parameters begin l = 0 u = 1 - end - @structural_parameters begin - z = ShiftIndex() - rng = Random.Xoshiro() + seed = 1 end @equations begin - y(z) ~ l + (u-l)*Symbolics.term(rand, rng; type=Real) output.u ~ y + n(z) ~ l + (u-l)*Symbolics.term(seeded_rand, seed, t; type=Real) + # y(z) ~ l + (u-l)*Symbolics.term(rand, rng; type=Real) + + if additive + y(z) ~ input.u(z) + n(z) + 1e-100*y(z-1) # The 0*y(z-1) is a workaround for a bug in the compiler, to force the y variable to be a discrete-time state variable + else + y(z) ~ n(z) + 1e-100*y(z-1) + end end end +""" + GenericNoise() + +A discrete-time noise source that at each clock tick returns a random value distributed according to the provided distribution. + +# Structural parameters +- `rng`: A random number generator, defaults to `Random.Xoshiro()`. +- `z`: The `ShiftIndex` used to indicate clock partition. +- `d`: The distribution to sample from.` + +# Connectors: +- `output` +""" +# @mtkmodel GenericNoise begin +# @components begin +# output = RealOutput() +# end +# @variables begin +# y(t), [description = "Output variable"] +# end +# @structural_parameters begin +# z = ShiftIndex() +# rng = Random.Xoshiro() +# d +# end +# @equations begin +# y(z) ~ Symbolics.term(rand, rng, d; type=Real) +# output.u ~ y +# end +# end + """ ZeroOrderHold() diff --git a/test/test_discrete_blocks.jl b/test/test_discrete_blocks.jl index 2871cb8..9707387 100644 --- a/test/test_discrete_blocks.jl +++ b/test/test_discrete_blocks.jl @@ -279,9 +279,10 @@ end ## Noise sources # ============================================================================== using ModelingToolkitStandardLibrary.Blocks +using Statistics @testset "NormalNoise" begin - k = ShiftIndex(Clock(1)) + k = ShiftIndex(Clock(0.01)) @mtkmodel NoiseModel begin @components begin @@ -301,7 +302,35 @@ using ModelingToolkitStandardLibrary.Blocks prob = ODEProblem(ssys, [], (0.0, 10.0)) sol = solve(prob, Tsit5()) @test !all(iszero, sol.u) - # TODO: add tests + tv = 0:k.clock.dt:sol.t[end] + @test std(sol(tv, idxs = m.plant.u)) ≈ 1 rtol=0.1 + @test mean(sol(tv, idxs = m.plant.u)) ≈ 0 atol=0.08 +end + +@testset "UniformNoise" begin + k = ShiftIndex(Clock(0.01)) + + @mtkmodel NoiseModel begin + @components begin + noise = UniformNoise(z = k) + zoh = ZeroOrderHold(z = k) + plant = FirstOrder(T = 1e-4) # Included due to bug with only discrete-time systems + end + @equations begin + connect(noise.output, zoh.input) + connect(zoh.output, plant.input) + end + end + + @named m = NoiseModel() + m = complete(m) + ssys = structural_simplify(IRSystem(m)) + prob = ODEProblem(ssys, [], (0.0, 10.0)) + sol = solve(prob, Tsit5()) + @test !all(iszero, sol.u) + tv = 0:k.clock.dt:sol.t[end] + @test minimum(sol(tv, idxs = m.plant.u)) ≈ 0 atol=0.02 + @test maximum(sol(tv, idxs = m.plant.u)) ≈ 1 atol=0.02 end From 0a98bf476d12a81becebc0a4d1cd7fb88a494427 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Tue, 3 Sep 2024 14:46:06 +0200 Subject: [PATCH 2/5] add quantization --- Project.toml | 2 +- src/ModelingToolkitSampledData.jl | 2 +- src/discrete_blocks.jl | 50 ++++++++++++++++++++++++++++++- test/test_discrete_blocks.jl | 39 ++++++++++++++++++++++++ 4 files changed, 90 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 7c72009..378659a 100644 --- a/Project.toml +++ b/Project.toml @@ -16,7 +16,7 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" DiffEqCallbacks = "~3.8" julia = "1.10" StableRNGs = "1" -JuliaSimCompiler = "1.19" +JuliaSimCompiler = "0.1.19" ModelingToolkit = "9" ModelingToolkitStandardLibrary = "2" OrdinaryDiffEq = "6.89" diff --git a/src/ModelingToolkitSampledData.jl b/src/ModelingToolkitSampledData.jl index 58ce736..c578b54 100644 --- a/src/ModelingToolkitSampledData.jl +++ b/src/ModelingToolkitSampledData.jl @@ -7,7 +7,7 @@ export get_clock export DiscreteIntegrator, DiscreteDerivative, Delay, Difference, ZeroOrderHold, Sampler, ClockChanger, DiscretePIDParallel, DiscretePIDStandard, DiscreteStateSpace, - DiscreteTransferFunction, NormalNoise, UniformNoise + DiscreteTransferFunction, NormalNoise, UniformNoise, Quantization include("discrete_blocks.jl") end diff --git a/src/discrete_blocks.jl b/src/discrete_blocks.jl index 8d89fc6..099df90 100644 --- a/src/discrete_blocks.jl +++ b/src/discrete_blocks.jl @@ -876,4 +876,52 @@ See also [ControlSystemsMTK.jl](https://juliacontrol.github.io/ControlSystemsMTK # push!(eqs, input.u ~ u) # push!(eqs, output.u ~ y) # compose(ODESystem(eqs, t, sts, pars; name = name), input, output) -# end \ No newline at end of file +# end + +""" + Quantization + +A quantization block that quantizes the input signal to a specified number of bits. + +# Parameters: +- `y_max`: Upper limit of output +- `y_min`: Lower limit of output +- `bits`: Number of bits of quantization +- `quantized`: If quantization effects shall be computed. If false, the output is equal to the input, which may be useful for, e.g., linearization. + +# Connectors: +- `input` +- `output` + +# Variables +- `y`: Output signal, equal to `output.u` +- `u`: Input signal, equal to `input.u` +""" +@mtkmodel Quantization begin + @extend u, y = siso = SISO() + @structural_parameters begin + z = ShiftIndex() + end + @parameters begin + y_max = 1, [description = "Upper limit of output"] + y_min = -1, [description = "Lower limit of output"] + bits::Int = 8, [description = "Number of bits of quantization"] + quantized::Bool = true, [description = "If quantization effects shall be computed."] + end + begin + end + @equations begin + y(z) ~ ifelse(quantized == true, quantize(u(z), bits, y_min, y_max), u(z)) + end +end + +function quantize(u, bits, y_min, y_max) + d = y_max - y_min + y1 = clamp(u, y_min, y_max) + y2 = (y1 - y_min) / d # between 0 and 1 + levels = 2^Int(bits)-1 + y3 = round(y2 * levels) / levels # quantized between 0 and 1 + y4 = y3*d + y_min + return y4 +end +@register_symbolic quantize(u, bits, y_min, y_max) diff --git a/test/test_discrete_blocks.jl b/test/test_discrete_blocks.jl index 9707387..333782d 100644 --- a/test/test_discrete_blocks.jl +++ b/test/test_discrete_blocks.jl @@ -360,3 +360,42 @@ end # @test reduce(vcat, sol((0:10) .+ 1e-2))[:]≈[zeros(2); 1; zeros(8)] atol=1e-2 # end* + +@testset "quantization" begin + @info "Testing quantization" + + function test_quant(y_min, y_max, bits) + u = y_min:(1/2^bits):y_max + y = ModelingToolkitSampledData.quantize.(u, bits, y_min, y_max) + uy = unique(y) + @test uy ≈ range(y_min, stop=y_max, length=2^bits) + end + + test_quant(-1, 1, 2) # Symmetric + test_quant(-1, 2, 2) # Not symmetric + test_quant(-1, 1, 3) # Symmetric, uneven number of bits + test_quant(-5, -2, 2) # Only negative + test_quant(5, 12, 2) # Only positive + test_quant(-5, 12, 20) # Large number of bits + + z = ShiftIndex(Clock(0.1)) + @mtkmodel QuantizationModel begin + @components begin + input = Sine(amplitude=1, frequency=1) + quant = Quantization(; z, bits=2) + end + @equations begin + connect(input.output, quant.input) + end + end + @named m = QuantizationModel() + m = complete(m) + ssys = structural_simplify(IRSystem(m)) + prob = ODEProblem(ssys, [], (0.0, 10.0)) + sol = solve(prob, Tsit5(), dtmax=0.01) + y = sol[m.quant.y] + uy = unique(y) + @test length(uy) == 4 +end + + From 1a26c1e01a3bf82c5e5dd6b77974ff996d85bf94 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Tue, 3 Sep 2024 15:00:31 +0200 Subject: [PATCH 3/5] add stats --- Project.toml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 378659a..5324fa1 100644 --- a/Project.toml +++ b/Project.toml @@ -14,18 +14,18 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" [compat] DiffEqCallbacks = "~3.8" -julia = "1.10" -StableRNGs = "1" JuliaSimCompiler = "0.1.19" ModelingToolkit = "9" ModelingToolkitStandardLibrary = "2" OrdinaryDiffEq = "6.89" Random = "1" - +StableRNGs = "1" +julia = "1.10" [extras] ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["ControlSystemsBase", "Test"] +test = ["ControlSystemsBase", "Test", "Statistics"] From d8a093c13e93ffc64d3e5eb73577f3e02e85247d Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Tue, 3 Sep 2024 15:13:22 +0200 Subject: [PATCH 4/5] add quantization tutorial --- docs/src/tutorials/noise.md | 33 ++++++++++++++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/docs/src/tutorials/noise.md b/docs/src/tutorials/noise.md index b4d9f37..35d1980 100644 --- a/docs/src/tutorials/noise.md +++ b/docs/src/tutorials/noise.md @@ -119,4 +119,35 @@ Internally, a random number generator from [StableRNGs.jl](https://github.com/Ju 2. Multiple calls to the random number generator at the same time step all return the same number. ## Quantization -Not yet available. \ No newline at end of file +A signal may be quantized to a fixed number of levels (e.g., 8-bit) using the [`Quantization`](@ref) block. This may be used to simulate, e.g., the quantization that occurs in a AD converter. Below, we have a simple example where a sine wave is quantized to 2 bits (4 levels), limited between -1 and 1: +```@example +using ModelingToolkit, ModelingToolkitSampledData, OrdinaryDiffEq, Plots +using ModelingToolkit: t_nounits as t, D_nounits as D +z = ShiftIndex(Clock(0.1)) +@mtkmodel QuantizationModel begin + @components begin + input = Sine(amplitude=1.5, frequency=1) + quant = Quantization(; z, bits=2, y_min = -1, y_max = 1) + end + @variables begin + x(t) = 0 # Dummy variable to work around a bug for models without continuous-time state + end + @equations begin + connect(input.output, quant.input) + D(x) ~ 0 # Dummy equation + end +end +@named m = QuantizationModel() +m = complete(m) +ssys = structural_simplify(IRSystem(m)) +prob = ODEProblem(ssys, [], (0.0, 2.0)) +sol = solve(prob, Tsit5()) +plot(sol, idxs=m.input.output.u) +plot!(sol, idxs=m.quant.y, label="Quantized output") +``` + +The output of the quantizer is always between `y_min` and `y_max` inclusive, and the number of distinct levels it can take is `2^bits`. The possible values are given by +```@example +bits = 2; y_min = -1; y_max = 1 +collect(range(y_min, stop=y_max, length=2^bits)) +``` From dd692893510732caf14bbe3f4fd22b3a33aba0d4 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Tue, 3 Sep 2024 16:28:42 +0200 Subject: [PATCH 5/5] add quantization mode midtread --- docs/src/tutorials/noise.md | 18 +++++++++++++-- src/discrete_blocks.jl | 25 ++++++++++++++++----- test/test_discrete_blocks.jl | 43 +++++++++++++++++++++++++++++++++--- 3 files changed, 76 insertions(+), 10 deletions(-) diff --git a/docs/src/tutorials/noise.md b/docs/src/tutorials/noise.md index 35d1980..a724d1e 100644 --- a/docs/src/tutorials/noise.md +++ b/docs/src/tutorials/noise.md @@ -120,7 +120,7 @@ Internally, a random number generator from [StableRNGs.jl](https://github.com/Ju ## Quantization A signal may be quantized to a fixed number of levels (e.g., 8-bit) using the [`Quantization`](@ref) block. This may be used to simulate, e.g., the quantization that occurs in a AD converter. Below, we have a simple example where a sine wave is quantized to 2 bits (4 levels), limited between -1 and 1: -```@example +```@example QUANT using ModelingToolkit, ModelingToolkitSampledData, OrdinaryDiffEq, Plots using ModelingToolkit: t_nounits as t, D_nounits as D z = ShiftIndex(Clock(0.1)) @@ -146,8 +146,22 @@ plot(sol, idxs=m.input.output.u) plot!(sol, idxs=m.quant.y, label="Quantized output") ``` -The output of the quantizer is always between `y_min` and `y_max` inclusive, and the number of distinct levels it can take is `2^bits`. The possible values are given by + + +### Different quantization modes +With the default option `midrise = true`, the output of the quantizer is always between `y_min` and `y_max` inclusive, and the number of distinct levels it can take is `2^bits`. The possible values are given by ```@example bits = 2; y_min = -1; y_max = 1 collect(range(y_min, stop=y_max, length=2^bits)) ``` +Notably, these possible levels _do not include 0_. If `midrise = false`, a mid-tread quantizer is used instead. The two options are visualized below: +```@example QUANT +y_min = -1; y_max = 1; bits = 2 +u = y_min:0.01:y_max +y_mr = ModelingToolkitSampledData.quantize_midrise.(u, bits, y_min, y_max) +y_mt = ModelingToolkitSampledData.quantize_midtread.(u, bits, y_min, y_max) +plot(u, [y_mr y_mt], label=["Midrise" "Midtread"], xlabel="Input", ylabel="Output", framestyle=:zerolines, l=2, seriestype=:step) +``` +Note how the default mid-rise quantizer mode has a rise at the middle of the interval, while the mid-tread mode has a flat region (a tread) centered around the middle of the interval. + +The default option `midrise = true` includes both end points as possible output values, while `midrise = false` does not include the upper limit. \ No newline at end of file diff --git a/src/discrete_blocks.jl b/src/discrete_blocks.jl index 099df90..f6683ed 100644 --- a/src/discrete_blocks.jl +++ b/src/discrete_blocks.jl @@ -901,6 +901,7 @@ A quantization block that quantizes the input signal to a specified number of bi @extend u, y = siso = SISO() @structural_parameters begin z = ShiftIndex() + midrise = true end @parameters begin y_max = 1, [description = "Upper limit of output"] @@ -911,17 +912,31 @@ A quantization block that quantizes the input signal to a specified number of bi begin end @equations begin - y(z) ~ ifelse(quantized == true, quantize(u(z), bits, y_min, y_max), u(z)) + y(z) ~ ifelse(quantized == true, quantize(u(z), bits, y_min, y_max, midrise), u(z)) end end -function quantize(u, bits, y_min, y_max) +function quantize_midrise(u, bits, y_min, y_max) d = y_max - y_min y1 = clamp(u, y_min, y_max) y2 = (y1 - y_min) / d # between 0 and 1 - levels = 2^Int(bits)-1 - y3 = round(y2 * levels) / levels # quantized between 0 and 1 + Δ = 2^Int(bits)-1 + y3 = round(y2 * Δ) / Δ # quantized between 0 and 1 y4 = y3*d + y_min return y4 end -@register_symbolic quantize(u, bits, y_min, y_max) + +function quantize_midtread(u, bits, y_min, y_max) + Δ = (y_max - y_min) / (2^Int(bits)-1) + # clamp(Δ * floor(u / Δ + 0.5), y_min, y_max) + k = sign(u) * max(0, floor((abs(u) - Δ/2) / Δ + 1)) + y0 = sign(k) * (Δ/2 + Δ*(abs(k)-1/2)) + y1 = iszero(y0) ? zero(y0) : y0 # remove -0.0 + return clamp(y1, y_min, y_max - Δ/2) +end + +function quantize(u, bits, y_min, y_max, midrise) + midrise ? quantize_midrise(u, bits, y_min, y_max) : quantize_midtread(u, bits, y_min, y_max) +end + +@register_symbolic quantize(u::Real, bits::Real, y_min::Real, y_max::Real, midrise::Bool) diff --git a/test/test_discrete_blocks.jl b/test/test_discrete_blocks.jl index 333782d..18595e4 100644 --- a/test/test_discrete_blocks.jl +++ b/test/test_discrete_blocks.jl @@ -299,7 +299,7 @@ using Statistics @named m = NoiseModel() m = complete(m) ssys = structural_simplify(IRSystem(m)) - prob = ODEProblem(ssys, [], (0.0, 10.0)) + prob = ODEProblem(ssys, [m.noise.y(k-1) => 0], (0.0, 10.0)) sol = solve(prob, Tsit5()) @test !all(iszero, sol.u) tv = 0:k.clock.dt:sol.t[end] @@ -325,7 +325,7 @@ end @named m = NoiseModel() m = complete(m) ssys = structural_simplify(IRSystem(m)) - prob = ODEProblem(ssys, [], (0.0, 10.0)) + prob = ODEProblem(ssys, [m.noise.y(k-1) => 0], (0.0, 10.0)) sol = solve(prob, Tsit5()) @test !all(iszero, sol.u) tv = 0:k.clock.dt:sol.t[end] @@ -366,7 +366,7 @@ end function test_quant(y_min, y_max, bits) u = y_min:(1/2^bits):y_max - y = ModelingToolkitSampledData.quantize.(u, bits, y_min, y_max) + y = ModelingToolkitSampledData.quantize_midrise.(u, bits, y_min, y_max) uy = unique(y) @test uy ≈ range(y_min, stop=y_max, length=2^bits) end @@ -378,6 +378,24 @@ end test_quant(5, 12, 2) # Only positive test_quant(-5, 12, 20) # Large number of bits + + + function test_quant2(y_min, y_max, bits) + u = y_min:(1/2^bits):y_max + y = ModelingToolkitSampledData.quantize_midtread.(u, bits, y_min, y_max) + uy = unique(y) + # @test (2^bits - 2 <= length(uy) <= 2^bits) # This check is not reliable since there might be one ulp difference when the last clamp is applied + @test maximum(y) <= y_max + @test minimum(y) >= y_min + end + + test_quant2(-1, 1, 2) # Symmetric + test_quant2(-1, 2, 2) # Not symmetric + test_quant2(-1, 1, 3) # Symmetric, uneven number of bits + test_quant2(-5, -2, 2) # Only negative + test_quant2(5, 12, 2) # Only positive + test_quant2(-5, 12, 20) # Large number of bits + z = ShiftIndex(Clock(0.1)) @mtkmodel QuantizationModel begin @components begin @@ -396,6 +414,25 @@ end y = sol[m.quant.y] uy = unique(y) @test length(uy) == 4 + + + @mtkmodel QuantizationModel2 begin + @components begin + input = Sine(amplitude=1, frequency=1) + quant = Quantization(; z, bits=2, midrise=false) + end + @equations begin + connect(input.output, quant.input) + end + end + @named m = QuantizationModel2() + m = complete(m) + ssys = structural_simplify(IRSystem(m)) + prob = ODEProblem(ssys, [], (0.0, 10.0)) + sol = solve(prob, Tsit5(), dtmax=0.01) + y = sol[m.quant.y] + uy = unique(y) + @test length(uy) == 4 end