Monitoring Station Design

This tutorial runs a multi-objective trade study over instrument settings for a remote monitoring station. The goal is to recover a known environmental signal as faithfully as possible while keeping station cost low.

The underlying signal is fixed. Designs vary the instrument configuration used to observe it.

Run it yourself

The full runnable script is at examples/monitoring_study.py.

uv run python examples/monitoring_study.py

The problem

A reference environmental signal is measured at a remote station. The instrument settings include sample rate, ADC resolution, sensor quality, and the number of co-located sensors. These affect how closely the reference signal is recovered and how much the station costs.

Reference signal

A synthetic 24-hour signal with three components (numpy only):

• Diurnal cycle: smooth sinusoid with a 24-hour period.
• Short-period oscillations: two modes at periods of 5 minutes and 20 minutes, chosen to illustrate finite sampling resolution.
• Correlated noise: an AR(1) process (\(\varphi = 0.995\), decorrelation time \(\approx 200\) s) representing natural variability.

The signal is generated at 1-second resolution, with \(N = 86\,400\) samples for a full day.

Observation model

The simplified observation model applies four transformations:

1. Anti-alias filter. Brick-wall low-pass at the Nyquist frequency of the configured sample rate.
2. Decimation. Subsample to the configured interval.
3. Quantization. Finite ADC resolution maps continuous values to \(2^b\) discrete levels.
4. Sensor noise. Additive Gaussian noise whose standard deviation depends on sensor grade, reduced by \(\sqrt{n}\) averaging when multiple co-located sensors are used. The reduction assumes independent sensor errors.

Why the objectives conflict

• Faster sampling, more bits, and additional sensors change reconstruction error and the example's assumed station cost.
• Coarse sampling removes variation above the Nyquist frequency, while limited ADC resolution introduces quantization error.
• A budget constrains which combinations of sampling and sensor settings are feasible.

Reference signal generation

The plot below shows a 1-hour window of the latent truth, revealing the diurnal trend, the short-period oscillations, and the correlated noise component:

DURATION = 86400.0  # 24 hours [s]
DT_TRUE = 1.0  # true signal resolution [s]
N_TRUE = int(DURATION / DT_TRUE)
T_TRUE = np.arange(N_TRUE) * DT_TRUE  # time axis [s]

RNG = np.random.default_rng(42)


def generate_reference_signal() -> np.ndarray:
    """Build a synthetic 24-hour environmental signal.

    Components:
      - Diurnal cycle (period = 24 h)
      - Short-period oscillations (periods ≈ 5 min, 20 min)
      - Correlated noise (AR(1), φ = 0.995)

    Returns:
        1-D array of length N_TRUE.
    """
    t_hours = T_TRUE / 3600.0

    # Diurnal
    diurnal = 5.0 * np.sin(2.0 * np.pi * t_hours / 24.0)

    # Short-period oscillations (5 min and 20 min)
    t_min = T_TRUE / 60.0
    short = 1.2 * np.sin(2.0 * np.pi * t_min / 5.0) + 0.8 * np.cos(
        2.0 * np.pi * t_min / 20.0
    )

    # Correlated noise  (AR(1), φ = 0.995 → decorrelation ~ 200 s)
    phi = 0.995
    noise = np.empty(N_TRUE)
    noise[0] = RNG.normal()
    for i in range(1, N_TRUE):
        noise[i] = phi * noise[i - 1] + RNG.normal(scale=np.sqrt(1 - phi**2))

    return diurnal + short + noise


REFERENCE = generate_reference_signal()

Observation model

The observe() function implements the four-stage signal chain described above. Each stage is a few lines of numpy:

SENSOR_NOISE: dict[str, float] = {
    "field": 1.0,
    "lab": 0.1,
    "reference": 0.01,
}


def observe(
    signal: np.ndarray,
    sample_interval: int,
    adc_bits: int,
    sensor_grade: str,
    n_sensors: int,
    rng: np.random.Generator,
) -> np.ndarray:
    """Apply the instrument signal chain to a reference signal.

    Steps:
      1. Anti-alias low-pass filter (brick-wall in frequency domain).
      2. Decimation to the configured sample interval.
      3. Quantization to *adc_bits* resolution.
      4. Additive sensor noise, reduced by √n_sensors averaging.

    Args:
        signal: High-resolution reference signal.
        sample_interval: Seconds between samples (decimation factor).
        adc_bits: ADC bit depth.
        sensor_grade: One of "field", "lab", "reference".
        n_sensors: Number of co-located sensors to average.
        rng: Numpy random generator for reproducibility.

    Returns:
        Degraded observation array at the decimated sample rate.
    """
    # 1. Anti-alias filter (brick-wall at Nyquist of decimated rate)
    spectrum = np.fft.rfft(signal)
    freqs = np.fft.rfftfreq(len(signal), d=DT_TRUE)
    nyquist = 0.5 / sample_interval
    spectrum[freqs > nyquist] = 0.0
    filtered = np.fft.irfft(spectrum, n=len(signal))

    # 2. Decimate
    decimated = filtered[::sample_interval]

    # 3. Quantize
    sig_range = float(np.ptp(signal)) * 1.1  # 10 % headroom
    n_levels = 2**adc_bits
    step = sig_range / n_levels
    mid = float(np.mean(signal))
    quantized = np.round((decimated - mid) / step) * step + mid

    # 4. Sensor noise (averaged over n_sensors)
    sigma = SENSOR_NOISE[sensor_grade] / np.sqrt(n_sensors)
    return quantized + rng.normal(scale=sigma, size=len(quantized))

Simulator and scorer

The Simulator applies the observation model with a SHA-256-derived, per-configuration seed that is stable across Python processes. The Scorer computes three objectives:

Observable	Direction	What it measures
RMSE	minimize	RMS error after interpolating back to the reference grid
Spectral fidelity	maximize	Fraction of spectral power recovered in the 2.5–30 min band
Station cost	minimize	Additive cost model over sample rate, ADC, grade, and sensor count

class StationSimulator:
    """Simulator that applies the instrument signal chain.

    Returns the full-resolution reference signal as 'truth' and
    the degraded, decimated observations as 'observations'.
    """

    def generate(self, config: dict[str, Any]) -> tuple[Any, Any]:
        """Generate truth and degraded observations.

        Args:
            config: Must contain sample_interval, adc_bits,
                sensor_grade, and n_sensors.

        Returns:
            Tuple of (reference_signal, degraded_observations).
        """
        seed = stable_config_seed(config, stream="monitoring-observation")
        rng = np.random.default_rng(seed)
        obs = observe(
            REFERENCE,
            sample_interval=config["sample_interval"],
            adc_bits=config["adc_bits"],
            sensor_grade=config["sensor_grade"],
            n_sensors=config["n_sensors"],
            rng=rng,
        )
        return REFERENCE, {"observations": obs, "interval": config["sample_interval"]}


def _spectral_fidelity(truth: np.ndarray, reconstructed: np.ndarray) -> float:
    """Fraction of spectral power recovered in the 2.5-30 min band.

    Args:
        truth: Full-resolution reference signal.
        reconstructed: Signal interpolated back to the reference grid.

    Returns:
        Power ratio clamped to [0, 1].
    """
    ref_psd = np.abs(np.fft.rfft(truth)) ** 2
    rec_psd = np.abs(np.fft.rfft(reconstructed)) ** 2
    freqs = np.fft.rfftfreq(len(truth), d=DT_TRUE)
    band = (freqs >= 1 / 1800) & (freqs <= 1 / 150)  # 2.5 min - 30 min
    ref_power = float(np.sum(ref_psd[band]))
    rec_power = float(np.sum(rec_psd[band]))
    return min(rec_power / ref_power, 1.0) if ref_power > 0 else 0.0


def _station_cost(config: dict[str, Any]) -> float:
    """Additive cost model over instrument settings.

    Args:
        config: Instrument configuration dict.

    Returns:
        Total station cost in arbitrary units.
    """
    rate_cost = {1: 50, 5: 30, 15: 15, 60: 5, 300: 1}
    bits_cost = {8: 1, 12: 5, 16: 20, 24: 80}
    grade_cost = {"field": 10, "lab": 50, "reference": 200}
    return float(
        rate_cost[config["sample_interval"]]
        + bits_cost[config["adc_bits"]]
        + grade_cost[config["sensor_grade"]] * config["n_sensors"]
    )


class StationScorer:
    """Score observation quality and station cost."""

    def score(
        self,
        truth: Any,
        observations: Any,
        config: dict[str, Any],
    ) -> dict[str, float]:
        """Compute RMSE, spectral fidelity, and station cost.

        Args:
            truth: Full-resolution reference signal.
            observations: Dict with degraded signal and sample interval.
            config: Instrument configuration.

        Returns:
            Scores for rmse, spectral_fidelity, and station_cost.
        """
        obs = observations["observations"]
        interval = observations["interval"]

        # Interpolate observations back to the reference grid
        obs_times = np.arange(len(obs)) * interval
        reconstructed = np.interp(T_TRUE, obs_times, obs)

        # RMSE
        rmse = float(np.sqrt(np.mean((truth - reconstructed) ** 2)))

        # Spectral fidelity: power ratio in the short-period band
        fidelity = _spectral_fidelity(truth, reconstructed)

        # Station cost (arbitrary units)
        cost = _station_cost(config)

        return {
            "rmse": rmse,
            "spectral_fidelity": fidelity,
            "station_cost": float(cost),
        }

Each configuration is evaluated for one synthetic reference day and one seeded observation realization. The example therefore demonstrates API and objective construction; it does not estimate design performance across days, environmental regimes, calibration states, or repeated sensor-noise draws.

Define observables and factors

Observables declare what we measure and which direction is better:

observables = [
    Observable("rmse", Direction.MINIMIZE),
    Observable("spectral_fidelity", Direction.MAXIMIZE),
    Observable("station_cost", Direction.MINIMIZE),
]

Factors define the instrument search space. The full factorial grid is \(5 \times 4 \times 3 \times 4 = 240\) configurations:

factors = [
    Factor("sample_interval", FactorType.DISCRETE, levels=[1, 5, 15, 60, 300]),
    Factor("adc_bits", FactorType.DISCRETE, levels=[8, 12, 16, 24]),
    Factor("sensor_grade", FactorType.DISCRETE, levels=["field", "lab", "reference"]),
    Factor("n_sensors", FactorType.DISCRETE, levels=[1, 2, 4, 8]),
]

Run the sweep

build_grid with method="full" generates every combination. run_grid evaluates them all and returns a ResultsTable:

    grid = build_grid(factors, method="full")
    print(f"Full factorial grid: {len(grid)} configurations")

    results = run_grid(
        world=StationSimulator(),
        scorer=StationScorer(),
        grid=grid,
        observables=observables,
    )

Inspect the Pareto front

The Pareto front retains evaluated configurations for which no other design is at least as good on every objective and strictly better on one. It describes the sampled configurations under this signal and cost model.

    # Pareto front
    directions = [o.direction for o in observables]
    front_idx = extract_front(
        results.scores,
        directions,
    )
    print(f"Pareto front: {len(front_idx)} / {len(grid)} designs\n")

    print(
        f"{'Interval':>8s}  {'Bits':>4s}  {'Grade':>9s}  {'#Sens':>5s}  "
        f"{'RMSE':>6s}  {'Fidelity':>8s}  {'Cost':>6s}"
    )
    print("-" * 60)
    for i in front_idx:
        cfg = results.configs[i]
        rmse_val, fid, cost = results.scores[i]
        print(
            f"{cfg['sample_interval']:8d}  {cfg['adc_bits']:4d}  "
            f"{cfg['sensor_grade']:>9s}  {cfg['n_sensors']:5d}  "
            f"{rmse_val:6.3f}  {fid:8.4f}  {cost:6.0f}"
        )

    # Cheapest design on the front
    cheapest = front_idx[np.argmin(results.scores[front_idx, 2])]
    print(f"\nCheapest Pareto design: {results.configs[cheapest]}")
    print(
        f"  RMSE={results.scores[cheapest, 0]:.4f}  "
        f"fidelity={results.scores[cheapest, 1]:.4f}  "
        f"cost={results.scores[cheapest, 2]:.0f}"
    )

    # Best RMSE on the front
    best = front_idx[np.argmin(results.scores[front_idx, 0])]
    print(f"\nLowest-RMSE Pareto design: {results.configs[best]}")
    print(
        f"  RMSE={results.scores[best, 0]:.4f}  "
        f"fidelity={results.scores[best, 1]:.4f}  "
        f"cost={results.scores[best, 2]:.0f}"
    )

Best vs cheapest Pareto designs

The comparison below overlays the observed signal (coloured) on the ground-truth reference (grey) for the two extreme Pareto members:

Compare the retained short-period variation, quantization, and noise in the two configurations. Their differences depend on the assumed signal and instrument model.

Pareto front scatter

Each dot is one instrument configuration in a 3-D plot of RMSE, spectral fidelity and cost; highlighted points are nondominated within this grid. The three coordinates show the sampled budget trade-off.

Parallel coordinates

Every line is one design, coloured by Pareto membership. The coordinates compare reconstruction error, spectral fidelity, and assumed station cost.

RMSE strip plot

Scores split by Pareto membership. A design with higher RMSE may remain on the front because it has lower cost or higher spectral fidelity.

What to try next

• Use build_grid(factors, method="lhs", n_samples=60) for a Latin hypercube over continuous factor bounds.
• Wrap the sweep in a multi-phase Study. Screen first with a coarse grid, then refine the promising region.
• Add an Annotation for power consumption or maintenance cost from an external costing sheet.
• Swap the brick-wall filter for a Butterworth via scipy.signal to study filter roll-off effects.