TrackingDifferentiator

Tracking differentiator: position and rate estimates for a scalar measurement, with a bound on the rate of change of the rate estimate.

The estimate is driven toward the measurement by the time-optimal control law for a double integrator whose acceleration is bounded by r, so the rate output changes by at most $r T_s$ between consecutive samples. A discontinuity in the measurement, whether a step or an isolated outlier, therefore produces a bounded transient of finite duration, whereas a linear tracker responds in proportion to the size of the discontinuity: an isolated outlier of magnitude $\Delta$ changes the rate estimate of this block by $r T_s$ and that of an AlphaBetaFilter by $\beta \Delta / T_s$.

Within the bound the update is linear and coincides exactly with the critically damped AlphaBetaFilter: with $p = 1 - T_s / h_0$, the rate output equals, sample for sample, that of AlphaBetaFilter(alpha = 1 - p^2, beta = (1 - p)^2). Both poles are located at $p$, the response has no overshoot, the rate estimate is delayed by approximately $2 h_0$, and a constant acceleration $a$ biases it by approximately $2 h_0 a$. The filter factor h0 therefore determines the trade-off between noise attenuation and lag; r does not. Reducing r below the value the signal requires rate-limits the estimate rather than filtering it. The update remains linear while $|\dot u| \le r h_0$ and $|\ddot u| \le r$. Choosing $r \ge \max(a_{\max}, v_{\max} / h_0)$ with margin ensures that the block departs from the linear tracker only for discontinuities larger than $r h_0^2$.

The defaults h0 = Ts and r = Inf reduce the block to the backward difference $(u_k - u_{k-1}) / T_s$; both the filtering and the bound must be enabled explicitly.

When the measurement is oversampled and the signal is smooth, an AlphaBetaGammaFilter tuned to the same noise level responds faster and has no acceleration bias, at the cost of approximately 40 % overshoot on velocity steps and a longer settling time. This block is preferable when a monotone response or robustness to outliers is required.

The rate state initializes to zero under every initialization variant; the enum selects only the position state. Since the position advances by the previous rate, the first output equals the selected previous position.

Usage

DiscreteComponents.TrackingDifferentiator(r=Inf)

Parameters:

Connectors

• u - This connector represents a real signal as an input to a component (RealInput)
• y - This connector represents a real signal as an output from a component (RealOutput)
• rate - This connector represents a real signal as an output from a component (RealOutput)

Variables

Behavior

using DiscreteComponents #hide
using ModelingToolkit #hide
@variables r #hide
@named sys = DiscreteComponents.TrackingDifferentiator(r=r) #hide
full_equations(sys) #hide

Source

"""
Tracking differentiator: position and rate estimates for a scalar measurement, with a bound on the rate of change of the rate estimate.

The estimate is driven toward the measurement by the time-optimal control law for a double integrator whose acceleration is bounded by `r`, so the `rate` output changes by at most ``r T_s`` between consecutive samples. A discontinuity in the measurement, whether a step or an isolated outlier, therefore produces a bounded transient of finite duration, whereas a linear tracker responds in proportion to the size of the discontinuity: an isolated outlier of magnitude ``\Delta`` changes the rate estimate of this block by ``r T_s`` and that of an [`AlphaBetaFilter`](@ref) by ``\beta \Delta / T_s``.

Within the bound the update is linear and coincides exactly with the critically damped [`AlphaBetaFilter`](@ref): with ``p = 1 - T_s / h_0``, the `rate` output equals, sample for sample, that of `AlphaBetaFilter(alpha = 1 - p^2, beta = (1 - p)^2)`. Both poles are located at ``p``, the response has no overshoot, the rate estimate is delayed by approximately ``2 h_0``, and a constant acceleration ``a`` biases it by approximately ``2 h_0 a``. The filter factor `h0` therefore determines the trade-off between noise attenuation and lag; `r` does not. Reducing `r` below the value the signal requires rate-limits the estimate rather than filtering it. The update remains linear while ``|\dot u| \le r h_0`` and ``|\ddot u| \le r``. Choosing ``r \ge \max(a_{\max}, v_{\max} / h_0)`` with margin ensures that the block departs from the linear tracker only for discontinuities larger than ``r h_0^2``.

The defaults `h0 = Ts` and `r = Inf` reduce the block to the backward difference ``(u_k - u_{k-1}) / T_s``; both the filtering and the bound must be enabled explicitly.

When the measurement is oversampled and the signal is smooth, an [`AlphaBetaGammaFilter`](@ref) tuned to the same noise level responds faster and has no acceleration bias, at the cost of approximately 40 % overshoot on velocity steps and a longer settling time. This block is preferable when a monotone response or robustness to outliers is required.

The rate state initializes to zero under every `initialization` variant; the enum selects only the position state. Since the position advances by the previous rate, the first output equals the selected previous position.
"""
component TrackingDifferentiator@[input clk extends Discrete]
  "Measured value"
  u = RealInput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": -100, "y1": 450, "x2": 0, "y2": 550, "rot": 0}
      },
      "tags": []
    }
  }
  "Position estimate"
  y = RealOutput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": 1000, "y1": 170, "x2": 1100, "y2": 270, "rot": 0}
      },
      "tags": []
    }
  }
  "Rate estimate [units of u per second]"
  rate = RealOutput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": 1000, "y1": 640, "x2": 1100, "y2": 740, "rot": 0}
      },
      "tags": []
    }
  }
  "Sample interval [s]"
  structural parameter Ts::Real = SampleTime()
  "Filter factor [s], at least Ts. Within the linear region the double pole is located at 1 - Ts/h0 and the rate estimate is delayed by approximately 2*h0; the default gives the unfiltered backward difference"
  structural parameter h0::Real = Ts
  structural parameter initialization::InitialCondition = DiscreteComponents.InitialCondition.InitialOutput(y0 = 0.0)
  "Acceleration bound on the rate estimate [units of u per second squared]; the rate estimate changes by at most r*Ts per sample. Inf disables the bound"
  parameter r::Real = Inf
  "Tracking error predicted h0 ahead, (y - u) + h0*rate, from the previous state and the current sample"
  variable yy::Real
  "Time-optimal rate target; equals rate + yy/h0 within the linear region"
  variable a::Real
  "Acceleration command, bounded by r"
  variable fh::Real
relations
  yy@clk = y@(clk-1) - u@clk + h0 * rate@(clk-1)
  a@clk = ifelse(abs(yy@clk) > r * h0 ^ 2, rate@(clk-1) + (sqrt((r * h0) ^ 2 + 8 * r * abs(yy@clk)) - r * h0) / 2 * sign(yy@clk), rate@(clk-1) + yy@clk / h0)
  fh@clk = ifelse(abs(a@clk) > r * h0, -r * sign(a@clk), -a@clk / h0)
  y@clk = y@(clk-1) + Ts * rate@(clk-1)
  rate@clk = rate@(clk-1) + Ts * fh@clk
  initial rate@(clk-1) = 0
  switch initialization
    case InitialOutput
      initial y@(clk-1) = initialization.y0
    case SteadyState
      initial y@(clk-1) = u@clk
    case InitialState
      initial y@(clk-1) = initialization.x0
  end
metadata {
  "Dyad": {"icons": {"default": "dyad://DiscreteComponents/TrackingDifferentiator.svg"}}
}
end

"""
Tracking differentiator: position and rate estimates for a scalar measurement, with a bound on the rate of change of the rate estimate.

The estimate is driven toward the measurement by the time-optimal control law for a double integrator whose acceleration is bounded by `r`, so the `rate` output changes by at most ``r T_s`` between consecutive samples. A discontinuity in the measurement, whether a step or an isolated outlier, therefore produces a bounded transient of finite duration, whereas a linear tracker responds in proportion to the size of the discontinuity: an isolated outlier of magnitude ``\Delta`` changes the rate estimate of this block by ``r T_s`` and that of an [`AlphaBetaFilter`](@ref) by ``\beta \Delta / T_s``.

Within the bound the update is linear and coincides exactly with the critically damped [`AlphaBetaFilter`](@ref): with ``p = 1 - T_s / h_0``, the `rate` output equals, sample for sample, that of `AlphaBetaFilter(alpha = 1 - p^2, beta = (1 - p)^2)`. Both poles are located at ``p``, the response has no overshoot, the rate estimate is delayed by approximately ``2 h_0``, and a constant acceleration ``a`` biases it by approximately ``2 h_0 a``. The filter factor `h0` therefore determines the trade-off between noise attenuation and lag; `r` does not. Reducing `r` below the value the signal requires rate-limits the estimate rather than filtering it. The update remains linear while ``|\dot u| \le r h_0`` and ``|\ddot u| \le r``. Choosing ``r \ge \max(a_{\max}, v_{\max} / h_0)`` with margin ensures that the block departs from the linear tracker only for discontinuities larger than ``r h_0^2``.

The defaults `h0 = Ts` and `r = Inf` reduce the block to the backward difference ``(u_k - u_{k-1}) / T_s``; both the filtering and the bound must be enabled explicitly.

When the measurement is oversampled and the signal is smooth, an [`AlphaBetaGammaFilter`](@ref) tuned to the same noise level responds faster and has no acceleration bias, at the cost of approximately 40 % overshoot on velocity steps and a longer settling time. This block is preferable when a monotone response or robustness to outliers is required.

The rate state initializes to zero under every `initialization` variant; the enum selects only the position state. Since the position advances by the previous rate, the first output equals the selected previous position.
"""
component TrackingDifferentiator
  "Measured value"
  u = RealInput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": -100, "y1": 450, "x2": 0, "y2": 550, "rot": 0}
      },
      "tags": []
    }
  }
  "Position estimate"
  y = RealOutput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": 1000, "y1": 170, "x2": 1100, "y2": 270, "rot": 0}
      },
      "tags": []
    }
  }
  "Rate estimate [units of u per second]"
  rate = RealOutput@[clk]() {
    "Dyad": {
      "placement": {
        "diagram": {"iconName": "default", "x1": 1000, "y1": 640, "x2": 1100, "y2": 740, "rot": 0}
      },
      "tags": []
    }
  }
  "Sample interval [s]"
  structural parameter Ts::Real = SampleTime()
  "Filter factor [s], at least Ts. Within the linear region the double pole is located at 1 - Ts/h0 and the rate estimate is delayed by approximately 2*h0; the default gives the unfiltered backward difference"
  structural parameter h0::Real = Ts
  structural parameter initialization::InitialCondition = DiscreteComponents.InitialCondition.InitialOutput(y0 = 0.0)
  "Acceleration bound on the rate estimate [units of u per second squared]; the rate estimate changes by at most r*Ts per sample. Inf disables the bound"
  parameter r::Real = Inf
  "Tracking error predicted h0 ahead, (y - u) + h0*rate, from the previous state and the current sample"
  variable yy::Real
  "Time-optimal rate target; equals rate + yy/h0 within the linear region"
  variable a::Real
  "Acceleration command, bounded by r"
  variable fh::Real
relations
  yy@clk = y@(clk-1) - u@clk + h0 * rate@(clk-1)
  a@clk = ifelse(abs(yy@clk) > r * h0 ^ 2, rate@(clk-1) + (sqrt((r * h0) ^ 2 + 8 * r * abs(yy@clk)) - r * h0) / 2 * sign(yy@clk), rate@(clk-1) + yy@clk / h0)
  fh@clk = ifelse(abs(a@clk) > r * h0, -r * sign(a@clk), -a@clk / h0)
  y@clk = y@(clk-1) + Ts * rate@(clk-1)
  rate@clk = rate@(clk-1) + Ts * fh@clk
  initial rate@(clk-1) = 0
  switch initialization
    case InitialOutput
      initial y@(clk-1) = initialization.y0
    case SteadyState
      initial y@(clk-1) = u@clk
    case InitialState
      initial y@(clk-1) = initialization.x0
  end
metadata {
  "Dyad": {"icons": {"default": "dyad://DiscreteComponents/TrackingDifferentiator.svg"}}
}
end

Test Cases

No test cases defined.

Related

• Examples
• Experiments
• Analyses