Merge pull request #21 from JuliaComputing/ADeffects

add Sampler with AD effects

Author: baggepinnen

docs/src/tutorials/noise.md

@@ -117,6 +117,7 @@ z = ShiftIndex(Clock(0.1))
     @equations begin
         connect(input.output, quant.input)
         D(x) ~ 0 # Dummy equation
+        # D(x) ~ Hold(quant.y) # Dummy equation
     end
 end
 @named m = QuantizationModel()
@@ -147,3 +148,38 @@ plot(u, [y_mr y_mt], label=["Midrise" "Midtread"], xlabel="Input", ylabel="Outpu
 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.
+
+
+## Sampling with AD effects
+The block [`SampleWithADEffects`](@ref) combines an ideal [`Sampler`](@ref), a [`NormalNoise](@ref) and a [`Quantization`](@ref) block to simulate the undesirable but practically occurring effects of sampling, noise and quantization in an AD converter. The block has the connectors `input` and `output`, where the input is the continuous-time signal to be sampled, and the output is the quantized, noisy signal. Example:
+
+```@example QUANT
+@mtkmodel PracticalSamplerModel begin
+    @components begin
+        input = Sine(amplitude=1.2, frequency=1, smooth=false)
+        sampling = SampleWithADEffects(; dt=0.03, bits=3, y_min = -1, y_max = 1, sigma = 0.1, noisy = true, quantized=true, midrise=true)
+    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, sampling.input)
+        D(x) ~ Hold(sampling.y) # Dummy equation
+    end
+end
+@named m = PracticalSamplerModel()
+m = complete(m)
+ssys = structural_simplify(IRSystem(m))
+prob = ODEProblem(ssys, [m.sampling.noise.y(z-1) => 0], (0.0, 2.0))
+sol = solve(prob, Tsit5())
+plot(sol, idxs=m.input.output.u)
+plot!(sol, idxs=m.sampling.y, label="AD converted output")
+```
+
+Both quantization and noise addition are optional and turned off by default. In the example above, we turn them on with keywords `noisy = true` and `quantized = true`. The noise is Gaussian white noise with standard deviation `sigma`, and the quantization is a 3-bit midrise quantizer (8 output levels) with limits `y_min` and `y_max`. Limits have to be provided when quantization is used. The `dt` parameter is the sampling time, if left unspecified, it will be inferred from context.
+
+Things to notice in the plot:
+- The sampled signal is saturated at the quantization limits ±1.
+- The noise is added to the signal before quantization, which means that the sampled signal has ``2^\text{bits}`` distinct output levels only.
+- 0 is not a possible output value. In situations where 0 is an important value (such as in the presence of integration of a quantized value that is expected to be close to 0), the mid-tread quantizer should be used instead by passing `midrise = false`.
+

src/ModelingToolkitSampledData.jl

@@ -5,7 +5,7 @@ using StableRNGs
 
 export get_clock
 export DiscreteIntegrator, DiscreteDerivative, Delay, Difference, ZeroOrderHold, Sampler,
-       ClockChanger,
+       ClockChanger, SampleWithADEffects,
        DiscretePIDParallel, DiscretePIDStandard, DiscreteStateSpace,
        DiscreteTransferFunction, NormalNoise, UniformNoise, Quantization,
        DiscreteSlewRateLimiter, ExponentialFilter

src/discrete_blocks.jl

@@ -974,11 +974,9 @@ This block is not differentiable, the derivative is zero everywhere exect for at
     @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"]
+        bits = 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, midrise), u(z))
     end
@@ -1080,4 +1078,58 @@ Discrete-time On-Off controller with hysteresis. The controller switches between
             y(z) ~ k*(2*s(z) - 1)
         end
     end
+end
+
+"""
+    SampleWithADEffects(quantized = false, noisy = false)
+
+A sampler with additional effects that appear in practical systems, such as measurement noise and quantization.
+
+The operations occur in the order
+1. Sampling
+2. Noise addition
+3. Quantization
+
+# Structural parameters:
+- `quantized`: If true, the output is quantized. When this option is used, the output is quantized to the number of bits specified by the `bits` parameter. The quantization is midrise if `midrise = true`, otherwise it is midtread. The output is also limited to the range `[y_min, y_max]`.
+- `noisy`: If true, the output is corrupted by additive white Gaussian noise with standard deviation `sigma` (defaults to 0.1). If `noisy = false`, the noise block is a unit gain.
+- `dt`: Sample interval of the sampler. If not specified, the sample interval is inferred from the clock of the system.
+- `clock`: Clock signal of the system. If not specified, the sample interval is inferred from the clock of the system. If `clock` is specified, the parameter `dt` has no effect.
+
+# Parameters:
+- `y_min`: Lower limit of output, defaults to -1. Only used if `quantized = true`.
+- `y_max`: Upper limit of output, defaults to 1. Only used if `quantized = true`.
+- `bits`: Number of bits of quantization, defaults to 8 (256 output levels between `y_min` and `y_max`). Only used if `quantized = true`.
+- `sigma`: Standard deviation of the additive Gaussian noise, defaults to 0.1. Only used if `noisy = true`.
+"""
+@mtkmodel SampleWithADEffects begin
+    @extend input, output = siso = SISO()
+    @structural_parameters begin
+        dt = nothing
+        clock = (dt === nothing ? InferredDiscrete() : Clock(dt))
+        midrise = true
+        quantized = false
+        noisy = false
+    end
+    @parameters begin
+        y_min = -1, [description = "Lower limit of output"]
+        y_max = 1, [description = "Upper limit of output"]
+        bits = 8, [description = "Number of bits of quantization"]
+        sigma = 0.1, [description = "Standard deviation of the additive noise"]
+    end
+    @components begin
+        sampler = Sampler(; clock)
+        if noisy
+            noise = NormalNoise(; sigma, additive = true)
+        else
+            noise = Gain(; k = 1)
+        end
+        quantization = Quantization(; bits, y_min, y_max, midrise, quantized)
+    end
+    @equations begin
+        connect(input, sampler.input)
+        connect(sampler.output, noise.input)
+        connect(noise.output, quantization.input)
+        connect(quantization.output, output)
+    end
 end

test/test_discrete_blocks.jl

@@ -538,4 +538,31 @@ end
     @test sol(0.999, idxs=m.filter.y) == 0
     @test sol(1.1, idxs=m.filter.y) > 0
 
+end
+
+@testset "sampling with AD effects" begin
+    @info "Testing sampling with AD effects"
+    z = ShiftIndex()
+    @mtkmodel PracticalSamplerModel begin
+        @components begin
+            input = Sine(amplitude=1.2, frequency=1, smooth=false)
+            sampling = SampleWithADEffects(; dt=0.03, bits=3, y_min = -1, y_max = 1, sigma = 0.1, noisy = true, quantized=true, midrise=true)
+        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, sampling.input)
+            D(x) ~ Hold(sampling.y) # Dummy equation
+        end
+    end
+    @named m = PracticalSamplerModel()
+    m = complete(m)
+    ssys = structural_simplify(IRSystem(m))
+    prob = ODEProblem(ssys, [m.sampling.noise.y(z-1) => 0], (0.0, 2.0))
+    sol = solve(prob, Tsit5())
+    
+    @test length(unique(sol[m.sampling.y])) == 8
+    @test maximum(abs, sol(0:0.03:1, idxs=m.sampling.y) - sol(0:0.03:1, idxs=m.sampling.u)) < 0.3
+
 end