Piecewise-Linear Correction for Thermal Oscillator Drift

A single affine model — one slope, one intercept, fit once against a reference clock — is the default treatment in clock drift correction, and it is correct as long as the edge device’s crystal oscillator holds a stable temperature. It does not on a gateway bolted to a weld cell enclosure, mounted above an induction furnace, or living in an unconditioned cabinet that swings 20–40°C between a cold start and full production. Quartz crystal frequency is a function of temperature, so the clock’s rate error is not constant across a shift — it rises as the cabinet heats through the first two hours, plateaus during steady production, and falls again as the line idles. A model that assumes one slope for the whole window either over-corrects the middle of the shift or under-corrects the edges, and either failure shifts a downtime event across a boundary that shift boundary logic is relying on being exact.

Offset versus time under a thermal cycle: a single linear fit diverges where a piecewise fit tracks the curve Measured clock offset rises, plateaus, then falls as an enclosure heats and cools over an eight-hour shift. A single straight-line fit from the first to the last point misses the mid-shift bulge and the late decline. A three-segment, continuous piecewise-linear fit with two breakpoints tracks the measured curve closely throughout. device time into shift (hours) → offset = reference − device (µs) 0h 8h measured offset · heats through hour 2, holds, cools after hour 6 single linear fit residual: line under-predicts offset by ≈110 µs at hour 4 breakpoint 1 breakpoint 2 piecewise-linear fit (continuous, 3 segments)

Why a constant slope fails under thermal cycling

Quartz oscillator frequency versus temperature follows a roughly cubic curve for the AT-cut crystals common in industrial gateways, with a turnover point near 25–30°C where the rate error is minimum and steeper error away from it. An enclosure that starts at ambient and heats past that turnover as the equipment around it runs introduces a rate error that itself changes shape over the shift, not just magnitude. Refitting the two-point affine model from clock drift correction — slope and intercept from a single start/end heartbeat pair — captures only the average rate error across the whole window. As the diagram shows, the true offset bulges upward through the heating phase and sags back down as the cabinet cools late in the shift; a straight line drawn between the endpoints misses both, leaving residual error of a hundred microseconds or more exactly where recalibration would be needed to keep gap-filling algorithms working against a trustworthy time index.

import numpy as np


def residual_of_single_fit(
    device_hours: np.ndarray, offset_us: np.ndarray
) -> np.ndarray:
    """Quantify how badly one affine fit tracks a thermally-cycled offset curve."""
    slope, intercept = np.polyfit(device_hours, offset_us, deg=1)
    predicted = slope * device_hours + intercept
    return offset_us - predicted  # non-zero mid-window residual signals thermal cycling


hours = np.array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0])
offsets_us = np.array([20.0, 90.0, 150.0, 175.0, 185.0, 180.0, 160.0, 110.0, 70.0])
residual = residual_of_single_fit(hours, offsets_us)
# residual peaks near hour 4 at roughly +100us: the single slope cannot see the bulge

Segmented fitting with breakpoint detection

Rather than one global slope, fit a sequence of short affine segments and re-anchor whenever the residual against the current segment exceeds a tolerance tied to the downstream error budget — typically the same offset budget used when configuring NTP and PTP sync at the edge. This greedy sequential-segmentation approach is deterministic and reproducible under replay, unlike a global change-point search: it only ever looks backward from the current point, so re-running it on the same raw heartbeats always produces the same breakpoints.

from dataclasses import dataclass, field

import numpy as np


@dataclass
class Segment:
    """One continuous affine piece of the piecewise-linear drift model."""
    t_start: float          # device hours at the anchor of this segment
    y_anchor: float          # offset (us) at t_start, fixed for continuity
    slope_us_per_hr: float = 0.0
    points: list = field(default_factory=list)


def fit_piecewise(
    device_hours: np.ndarray,
    offset_us: np.ndarray,
    max_residual_us: float = 40.0,
) -> list[Segment]:
    """Greedy, continuity-constrained piecewise-linear fit over a thermal cycle.

    Extends the current segment while its residual against a slope-only fit
    (anchored at the segment start, so no intercept is free) stays within
    `max_residual_us`. On breach, the segment closes at the last good point
    and a new one opens there — never introducing a jump.
    """
    segments: list[Segment] = []
    seg = Segment(t_start=device_hours[0], y_anchor=offset_us[0])
    seg.points.append((device_hours[0], offset_us[0]))

    for t, y in zip(device_hours[1:], offset_us[1:]):
        tau = np.array([p[0] - seg.t_start for p in seg.points] + [t - seg.t_start])
        y_vals = np.array([p[1] for p in seg.points] + [y])
        # Closed-form slope through the fixed anchor (y0 at tau=0):
        # minimizes sum((y0 + slope*tau - y)^2) over slope alone.
        candidate_slope = float(np.dot(tau, y_vals - seg.y_anchor) / np.dot(tau, tau))
        predicted = seg.y_anchor + candidate_slope * tau
        if np.max(np.abs(predicted - y_vals)) > max_residual_us:
            seg.slope_us_per_hr = _closed_form_slope(seg)
            segments.append(seg)
            seg = Segment(t_start=seg.points[-1][0], y_anchor=seg.points[-1][1])
            seg.points.append((t, y))
        else:
            seg.points.append((t, y))

    seg.slope_us_per_hr = _closed_form_slope(seg)
    segments.append(seg)
    return segments


def _closed_form_slope(seg: Segment) -> float:
    tau = np.array([p[0] - seg.t_start for p in seg.points])
    y_vals = np.array([p[1] - seg.y_anchor for p in seg.points])
    denom = float(np.dot(tau, tau))
    return float(np.dot(tau, y_vals) / denom) if denom > 0 else 0.0

max_residual_us should be set from the same error budget the correction feeds — 40 µs is generous for a TSDB write cadence in seconds, tight for a microsecond-class control loop. A shorter tolerance yields more, shorter segments and tracks the thermal curve more closely at the cost of more frequent recalibration heartbeats against the PTP/NTP reference.

Continuity constraints across segment boundaries

The closed-form slope above deliberately fixes each segment’s intercept to the previous segment’s last corrected value (y_anchor) rather than fitting slope and intercept freely. This is what keeps the model C0C^0-continuous: two segments meeting at a breakpoint agree exactly at that point, so applying the correction never produces the kind of instantaneous step that clock drift correction treats as an NTP re-anchor event. An unconstrained per-segment least-squares fit — treating each window as an independent regression — would instead leave small jumps at every breakpoint, and those jumps would silently violate the monotonicity guarantee the correction engine depends on.

def correct_with_segments(
    device_hours: np.ndarray, segments: list[Segment]
) -> np.ndarray:
    """Apply the continuous piecewise model to a raw device-time array."""
    corrected = np.empty_like(device_hours, dtype=float)
    bounds = [s.t_start for s in segments] + [np.inf]

    for seg, hi in zip(segments, bounds[1:]):
        mask = (device_hours >= seg.t_start) & (device_hours < hi)
        tau = device_hours[mask] - seg.t_start
        corrected[mask] = seg.y_anchor + seg.slope_us_per_hr * tau

    return corrected  # continuous by construction: no step at any breakpoint

Because each segment’s slope is a rate error rather than a fixed offset, the correction still grows with distance from the segment’s own anchor — the same principle as the single-window affine model, just re-anchored more often. Validate continuity directly by checking that the corrected value evaluated from the end of one segment equals the corrected value evaluated from the start of the next, to within floating-point tolerance; see handling floating-point drift in sensor readings for why an exact equality check on float64 offsets is the wrong test.

Gotchas & anti-patterns

  • Recalibrating on a fixed wall-clock schedule instead of on residual. A nightly recalibration misses a two-hour thermal transient that happens mid-shift; trigger new segments from the residual test, not a cron job.
  • Free-fitting slope and intercept per segment. This reintroduces the step discontinuity the whole approach exists to avoid — always anchor the new segment’s intercept to the previous segment’s endpoint.
  • Segments too short to estimate a slope reliably. A two-point segment fits the noise, not the trend; require a minimum point count or minimum duration before allowing a breakpoint.
  • Confusing thermal drift with an NTP step. A genuine step-correction from an intermittent NTP sync looks superficially like a sharp breakpoint; distinguish them by magnitude — a single-digit-millisecond discontinuity is a step, a smooth multi-minute bulge is thermal.
  • Ignoring hysteresis. A crystal’s frequency-temperature curve is not perfectly reversible on heating versus cooling; if the enclosure’s thermal cycle is repeatable shift to shift, log segment slopes by time-of-shift so hysteresis becomes a visible, correctable pattern rather than noise.

Quick reference

Regime Model Recalibration trigger
Stable-temperature cabinet (climate-controlled) Single affine fit Fixed interval (e.g. daily), residual rarely exceeds budget
Thermal ramp (cold start through warm-up) Piecewise-linear, short segments Residual against current segment exceeds error budget
Steady-state production (plateaued temperature) Piecewise-linear, long segments Residual test rarely fires; segments naturally lengthen
Cyclical thermal load (batch process, intermittent furnace) Piecewise-linear + logged hysteresis by shift phase Residual test, plus scheduled audit against PTP/NTP reference

Related in this section