Data Reconciliation

Overview

Plant measurements are noisy, and taken at face value they contradict the model: nominations do not close, balances do not balance, pressure drops do not match flows. Data reconciliation finds the smallest statistically weighted adjustment to the measurements that makes them satisfy the model equations. Because the model carries information the sensors do not, the reconciled estimates are more precise than the raw measurements — often by a factor of two or three.

difflow.reconciliation provides:

• Constrained weighted least squares against any differentiable residual function \(F(x, \theta) = 0\) — a flowsheet, a gas network, or an equation set you write yourself.

• Joint parameter estimation: mark a variable unmeasured and it is estimated rather than reconciled, with a standard error, in the same solve.

• Covariance of the estimates, from the inverse KKT matrix, and the sensitivity \(\partial \hat x/\partial y\) by automatic differentiation.

• Gross error detection: the global \(\chi^2\) test, the measurement test, and serial elimination.

• Observability and redundancy classification, run before the solve so an ill-posed problem raises a named error instead of returning NaN.

• Sensor placement: what a proposed meter would buy, before anyone buys it.

Everything that makes this work — the constraint Jacobian, the covariance, the sensitivities — is derivative information. A conventional flowsheet package hand-derives it per unit operation; a differentiable one gets it from jax.jacobian.

The module is domain-agnostic. difflow_gas.residuals shows how a plugin supplies its equation set; see examples/28_data_reconciliation.ipynb for the worked gas-network case.

---

Mathematical Formulation

The problem

\[\min_x \; (x - y)^T W (x - y) \quad \text{subject to} \quad F(x, \theta) = 0\]

where \(y\) are the measurements, \(W = \operatorname{diag}(1/\sigma_i^2)\), and \(F\) are the model equations. Entries with \(\sigma_i = \infty\) get \(W_{ii} = 0\): they are estimated, not reconciled. A finite \(\sigma\) on an unmetered variable acts as a Bayesian prior, which is the graceful way to handle a weakly identified parameter.

The KKT system

The first-order conditions are

\[\begin{split}\begin{bmatrix} W & A^T \\ A & 0 \end{bmatrix} \begin{bmatrix} \Delta x \\ \lambda \end{bmatrix} = \begin{bmatrix} -W(x - y) \\ -F(x) \end{bmatrix}, \qquad A = \frac{\partial F}{\partial x},\end{split}\]

iterated to convergence. This is Gauss–Newton: it drops the \(\sum_k \lambda_k \nabla^2 F_k\) term of the exact Newton Jacobian, which is where a pipe law’s non-smooth \(q|q|\) second derivative would enter. Exact gradients are recovered afterwards from one implicit-function-theorem correction that does use the full Jacobian, so jax.grad through reconcile is exact with respect to \(y\), \(\sigma\) and \(\theta\).

Solvability is observability

For \(W \succeq 0\), the KKT matrix is nonsingular iff \(A\) has full row rank and \(Z^T W Z \succ 0\) on a basis \(Z\) of \(\ker A\). With \(W\) diagonal and zero exactly on the unmeasured entries, the second condition collapses to

the unmeasured columns \(A_U\) must have full column rank

which is precisely the classical observability condition. classify runs this test first and raises ReconciliationStructureError naming the culprits.

Two consequences worth knowing:

• Boundary flows must be state variables, not fixed parameters. With them fixed, the node-balance block of \(A\) is the incidence matrix, whose rank is only \(n_{\text{nodes}} - 1\), and the KKT matrix is singular for purely structural reasons.

• Ranks come from svd(A), never from the eigenvalues of \(A^T A\). Squaring the matrix squares its condition number, and a structurally zero singular value reappears near \(\sqrt{\varepsilon}\,\sigma_{\max}\) — an unobservable system reported as full rank.

Covariance

\[\Sigma_{\hat x} = [K^{-1}]_{11}\]

For a fully measured problem this equals the textbook projection \(\Sigma - \Sigma A^T (A\Sigma A^T)^{-1} A \Sigma\), but unlike it, it stays well defined when a variable is unmeasured — so the standard error of an estimated parameter comes from the same expression as the reconciled variance of a metered flow. The adjustment covariance is \(\Sigma_{\text{adj}} = \Sigma - \Sigma_{\hat x}\).

measurement_sensitivity computes \(S = \partial \hat x/\partial y\) by differentiating the solver. For linear constraints \(S \Sigma S^T = \Sigma_{\hat x}\) exactly. For nonlinear ones the two differ by the curvature term the covariance formula drops: they agree to machine precision when the data are consistent (the multipliers vanish) and diverge in proportion to how inconsistent the data are. The classical formula is itself a linearization; the discrepancy is a useful diagnostic of how nonlinear the model is over the range the adjustments span.

Scaling

The KKT matrix mixes \(W\) (units of 1/variable²) with \(A\) (units of residual/variable), so its conditioning depends on the unit system — the same gas network posed in Pa and Pa² rather than bar and bar² is many orders worse conditioned, past what float64 resolves. Variables are scaled by \(d_i = \sigma_i\), which makes the scaled weight matrix a 0/1 mask (so \(1/\sigma^2\) is never evaluated and an infinite \(\sigma\) cannot produce a NaN), and residual rows are equilibrated to unit 2-norm. This is on by default; pass scaling=False only to observe what it prevents.

Gross error detection

The global test statistic is the optimal objective itself, distributed \(\chi^2\) on the degrees of redundancy \(m - \operatorname{rank}(A_U)\). The often-quoted form \(F(y)^T (A\Sigma A^T)^{-1} F(y)\) is a special case that cannot be evaluated at all when a variable is unmeasured, and carries the wrong degrees of freedom.

The measurement test standardizes each adjustment, \(z_i = (\hat x_i - y_i)/\sqrt{\Sigma_{\text{adj},ii}}\), which is standard normal under the null hypothesis. A measurement nothing checks has \(\Sigma_{\text{adj},ii} = 0\) and is reported as untestable rather than given a spurious \(z\). Because that and the redundancy classification are read from the same \(\Sigma_{\hat x}\), they cannot disagree.

Least squares smears a gross error across neighbouring measurements, so identification is reliable only where redundancy is high; serial_elimination discards the prime suspect and re-tests until the data are clean.

---

API Reference

reconcile

def reconcile(
    residual_fn,                    # F(x, params) -> (m,) Array
    y: Array,
    sigma: Array,                   # inf entry => unmeasured
    *,
    params=None, names=None, x0=None,
    unmeasured_init=None, unmeasured_scale=None,
    scaling: Scaling | bool = True,
    max_steps: int = 20, tol: float = 1e-9,
    method: str = "gauss_newton",
    check_structure: bool = True, rank_tol=None,
) -> ReconcileResult

Parameters:

• residual_fn — the model equations, JAX-traceable.

• y — measurements. Entries with infinite sigma are ignored and may be nan.

• sigma — standard deviations; inf marks a variable to estimate.

• params — extra argument threaded to residual_fn; jax.grad with respect to it answers how would the reconciled state move if this fixed parameter changed — a different question from estimating it.

• check_structure — run the observability test first. Disable only inside a jit/vmap sweep whose structure you have already validated.

Returns a ReconcileResult with x, x_named, adjustment, objective, covariance, std, converged, structure and scaling, plus a summary() table.

Other entry points

classify(residual_fn, x, sigma, *, scaling, names=None) -> StructureReport
reconciled_covariance(residual_fn, x, sigma, *, scaling) -> Array
measurement_sensitivity(residual_fn, y, sigma, *, x0, scaling) -> Array

global_test(result, alpha=0.05) -> GlobalTestResult
measurement_test(result, alpha=0.05, bonferroni=True) -> MeasurementTestResult
serial_elimination(residual_fn, y, sigma, *, alpha=0.05, max_removed=3) -> list

sensor_value(residual_fn, x, sigma, *, target, candidate, candidate_sigma) -> dict
sensor_ranking(residual_fn, x, sigma, *, target, candidates, candidate_sigma) -> list

monitor(residual_fn, measurements, sigma, *, names=None, alpha=0.05) -> MonitorResult
blame_concentration(suspects, window=15) -> (fraction, culprit)
reconcile_multi(residual_fn, measurements, sigma, *, shared, names=None) -> MultiReconcileResult

StructureReport.classes assigns every variable one of measured-redundant, measured-just-determined, unmeasured-observable or unmeasured-unobservable, and summary() prints the table.

---

Worked Example

A three-stream splitter whose meters do not close:

import jax.numpy as jnp
from difflow.reconciliation import reconcile, global_test, measurement_test

def balance(x, params=None):
    return jnp.array([x[0] - x[1] - x[2]])       # feed = top + bottom

y     = jnp.array([100.0, 62.0, 40.0])           # out by 2 units
sigma = jnp.array([2.0, 1.0, 1.0])               # the feed meter is worst

res = reconcile(balance, y, sigma, names=["feed", "top", "bottom"])
print(res.summary())
print(global_test(res))

The feed meter, being the least trusted, absorbs most of the adjustment; the reconciled standard deviations are all below the raw sigma. With one equation and three measurements the degree of redundancy is 1.

To estimate a quantity instead of reconciling it, give it an infinite sigma:

sigma = jnp.array([2.0, 1.0, jnp.inf])           # the bottoms meter failed
res = reconcile(balance, y, sigma, names=["feed", "top", "bottom"])
res.x_named["bottom"]                            # inferred from the balance
res.std["bottom"]                                # and its standard error

Ask for one unknown too many and the problem is diagnosed rather than silently failing:

from difflow.reconciliation import ReconciliationStructureError

sigma = jnp.array([2.0, jnp.inf, jnp.inf])
try:
    reconcile(balance, y, sigma, names=["feed", "top", "bottom"])
except ReconciliationStructureError as err:
    print(err)
# ReconciliationStructureError: reconciliation problem is not solvable:
# 1 unmeasured variable(s) cannot be determined from the constraints.
# Unobservable: top, bottom. ...

---

Gas Networks

difflow_gas supplies the equation set for transmission networks. difflow_gas.residuals.network_residuals is the single definition of that set — nodal balances, resistance laws, valve relations, and the compressor relation \(p_{\text{to}} = r\, p_{\text{from}}\). difflow_gas.verify is the reporting layer over it, unflattening the same residual vector into labelled dicts of floats.

verify reports every block except the compressor relation, which a sequential solve satisfies by construction — there is nothing to check. A reconciliation cannot drop it, and on the example network it turns out to be exactly what makes the loop observable.

import jax
import difflow_gas as dg

# A small looped network: src -> a -> (compressor) -> b -> c/d, with a cycle
RATIOS = {"cs1": 1.2}
net = dg.GasNetwork(
    arcs={"p1": ("src", "a", "pipe"), "cs1": ("a", "b", "compressor"),
          "p2": ("b", "c", "pipe"), "p3": ("b", "d", "pipe"),
          "p4": ("c", "d", "pipe")},
    beta={aid: dg.weymouth_beta(length_m=L, diameter_m=0.6, roughness_m=1e-4)
          for aid, L in [("p1", 20e3), ("p2", 40e3), ("p3", 60e3), ("p4", 80e3)]},
    supply_kg_s={"src": 120.0, "c": -50.0, "d": -70.0},
    pressure_bounds_bar={n: (30.0, 80.0) for n in ["src", "a", "b", "c", "d"]},
)
fs, dec = dg.build_network_flowsheet(net, root="src", p_slack_pa=60.0e5, ratios=RATIOS)
streams = fs.solve(tol=1e-12, max_iter=500)
p_true = dg.verify.node_pressures_bar(streams, dec)       # the simulated truth
q_true = dg.verify.arc_flows_kg_s(streams, dec)

layout = dg.gas_state_layout(net, efficiency_arcs=["p3"])   # estimate fouling
x_true = layout.pack(p_true, q_true, net.supply_kg_s, {"eta_p3": 1.0})
sigma  = dg.measurement_sigma(layout)                       # meter accuracies
y      = dg.perturb(x_true, sigma, jax.random.PRNGKey(0))   # simulated data

res = dg.reconcile_network(net, y, sigma, layout, ratios={"cs1": 1.2})
p_bar, q_kg_s, supply = dg.reconciled_values(res, layout)
# the reconciliation was free to move the nominations, so check against those
net_rec = dg.GasNetwork(arcs=net.arcs, beta=net.beta, supply_kg_s=supply)
dg.verify.residuals_from_values(p_bar, q_kg_s, net_rec).ok  # True

dg.monitor_network and dg.reconcile_network_multi are the campaign-scale versions of the same call, and take the same (network, ..., layout, ratios=...) arguments. All three fill in names and unmeasured_scale from the layout, so a gas reconciliation never has to restate them.

Measured nominations are deliberately kept out of GasNetwork, which rejects supplies that do not sum to zero — real nominations do not, and making them close is what the reconciliation is for.

See examples/28_data_reconciliation.ipynb for the full treatment: variance reduction, a biased flow meter found and eliminated, an unmeasured pipe-fouling factor estimated with its standard error, the observability boundary, and sensor placement.

Reconciling data vs. updating the model

Both are the same optimisation — a variable with a finite sigma is a measurement you may move at a cost, one with sigma = inf is a parameter you may move for free — so the interesting question is not how to update a model but when you are entitled to. Letting a parameter float absorbs whatever is wrong, including a broken sensor, and hands back a confident wrong number.

The practical discipline is two clocks: reconcile routinely with parameters fixed, so the global test stays a genuine instrument-health monitor; re-estimate parameters only as a deliberate campaign. The trigger is the pair of tests read together — persistent rejection with the same sensor blamed every day is an instrument fault, while persistent rejection with a wandering suspect is model drift.

The routine clock: monitor

monitor runs the first clock. It reconciles a sequence of data sets against one fixed model and keeps the diagnostics, which is what makes the resulting statistic series readable: the model and the sigmas never move, so every wobble in it comes from the data.

import difflow_gas as dg

layout0 = dg.gas_state_layout(net)               # the plain layout, no eta
sigma0 = dg.measurement_sigma(layout0)
x0_true = layout0.pack(p_true, q_true, net.supply_kg_s)
daily_measurements = [dg.perturb(x0_true, sigma0, jax.random.PRNGKey(k))
                      for k in range(6)]         # six days of simulated data

mon = dg.monitor_network(net, daily_measurements, sigma0, layout0, ratios=RATIOS)
mon.statistic          # the chi-squared series, shape (n_days,)
mon.suspects           # who the measurement test blamed each day
mon.diagnose(window=15)
# model drift: 93% of the last 15 steps reject, blame concentration 40%

monitor(residual_fn, measurements, sigma, names=...) is the domain-agnostic form; dg.monitor_network is the same call with the network’s residual closure and the layout’s names and scales filled in.

diagnose applies the rule above. blame_concentration is the measurement behind it — the fraction of the window blaming the single most-blamed sensor, counted over the whole window including the quiet days, so a campaign that rejects rarely does not read as a concentrated fault just because its few rejections agreed. A verdict of instrument fault names the culprit and means go and calibrate it; model drift is the one verdict that entitles you to the second clock. Both thresholds are arguments, so the rule can be tuned to a plant’s own noise.

A data set whose problem cannot be posed at all records failed on its step rather than aborting the campaign.

The campaign clock: reconcile_multi

One day’s data gives a noisy parameter estimate, so pool a window. reconcile_multi gives each data set its own copy of the plant state while the variables in shared appear once, estimated from all of them in a single solve:

from difflow.reconciliation import global_test

layout = dg.gas_state_layout(net, efficiency_arcs=["p3"])   # carries eta
sigma = dg.measurement_sigma(layout)
plain_layout, week_of_data = layout0, daily_measurements
window = [layout.embed(y, plain_layout) for y in week_of_data]

res = dg.reconcile_network_multi(
    net, window, sigma, layout, shared=["eta_p3"], ratios=RATIOS,
)
res.shared["eta_p3"], res.shared_std["eta_p3"]
res.states[0]                       # day 0's reconciled state
global_test(res)                    # on the pooled problem

Again reconcile_multi(residual_fn, measurements, sigma, shared=...) is the general form. GasStateLayout.embed re-packs measurements taken against a plainer layout into the one carrying the parameter, by name — the added entry is unmeasured, so it is filled with nan.

This is not the same as reconciling each day separately and averaging the estimates, in two ways that matter. The standard error it reports is that of the pooled estimate, roughly \(\sqrt{K}\) tighter, whereas averaging point estimates leaves you holding one day’s error bar for a quantity \(K\) of them informed. And because the parameter is one unknown rather than \(K\) private copies, the degrees of redundancy count it once — pooling recovers the \(K-1\) that separate estimations throw away, and a parameter too weakly identified to be recovered from one data set can still be observable from several. The structure check runs on the stacked problem and reports this.

The catch is lag: pooling estimates the average truth over its window, so choose a window short enough that the parameter is genuinely constant across it. A finite sigma on a shared variable is a prior, and it is applied once — \(K\) copies of one prior would count it \(K\) times.

examples/29_model_updating.ipynb works this through on a pipe that fouls over a 45-day campaign, including the case where a free parameter manufactures a fouling estimate out of a biased flow meter.

Tracking a drifting parameter: track_parameters

The two clocks are a discipline a person applies. A digital twin — a model kept current against the plant so that what you optimise, price or plan against is the plant as it is — has nobody to apply it, and the discipline does not survive being automated naively. Two things have to change.

The verdict has to gate the update, not advise it. diagnose() already decides whether a sensor or the model is at fault. Left as a report, nothing stops a scheduled re-estimation from running on a day when the honest answer was “go calibrate the dp meter”, and section 6 of examples/29_model_updating.ipynb shows what that costs: the free parameter absorbs part of the bias, the twin reports a confident fouling estimate for a pipe nobody inspected, and the \(\chi^2\) statistic falls while it happens. The model looks healthier as it gets wronger. No filter gain prevents this — a slow filter reaches the wrong answer gracefully. Only a gate prevents it.

The update has to have a memory. Pooling a window is a rectangular filter: hard edges, a ten-day-old period weighted like this morning’s, recomputed from scratch each time. Run it on a schedule and the estimate lurches as periods fall off the back of the window.

track_parameters is the two clocks written as an update law:

from difflow.reconciliation import (
    TrackerState, drift_std_from_time_constant, track_parameters,
)

def F(x, params):
    """The network residuals, with the tracked parameter passed through `params`."""
    return dg.residuals.network_residuals(
        x, net, layout0, ratios=RATIOS, efficiencies={"p3": params["eta"]},
    )

run = track_parameters(
    F, daily_measurements, sigma0,
    state=TrackerState.initial(["eta"], [1.0], std=[0.02]),
    drift_std=drift_std_from_time_constant(0.05, 30.0),   # 5% over a month
    names=layout0.names, unmeasured_scale=layout0.default_scale,
)

run.final.as_params()      # {'eta': ...} -- feeds a model, a Block, an LP
run.final.std              # {'eta': ...}
print(run.summary())

Each period does four things:

1. reconcile with the parameters frozen at the current estimate and record both gross-error tests — the routine clock, and the only reason the tests mean anything;

2. draw a verdict from the campaign so far and put it to update_gate;

3. time_update the tracker, whatever the verdict said;

4. only if the gate opened, estimate the parameters from this period and fold the result in with measurement_update.

The filter

The parameter is given the random walk

\[\theta_{k+1} = \theta_k + w_k, \qquad \operatorname{cov}(w_k) = Q\,\Delta t,\]

which is the same model difflow.mhe.augment_parameters puts on a drifting parameter, and drift_std = \(\sqrt{\operatorname{diag} Q}\) is the same knob as process_std there, carrying the same warning: too large and the parameter absorbs sensor noise, too small and a genuine drift is rejected. It is a required argument, never a default, because it is the bandwidth of the twin rather than a nuisance.

A rate is hard to have an opinion about; the question an engineer can answer is how far does this move, and over how long? drift_std_from_time_constant(spread, tau) converts one to the other — a random walk accumulates \(\sqrt{Q}\sqrt{t}\), so fouling that costs five percent of duty in a month is drift_std_from_time_constant(0.05, 30.0) on a daily clock.

No new estimator is involved. A single period’s reconciliation already returns both halves of a Kalman measurement update — the estimate, and the block of reconciled_covariance belonging to it — so the update is the combination of two Gaussians and all the physics stays inside reconcile:

\[S = P^- + R, \quad K = P^- S^{-1}, \quad \theta^+ = \theta^- + K(\hat\theta - \theta^-).\]

Four properties of that combination are worth stating, because they are what make it better than a rolling refit:

• The covariance is kept full. Correlated parameters have a difference variance a diagonal covariance gets wrong by a factor of a few, so the prior is a matrix. Pass one to TrackerState.initial(..., covariance=...) — difflow.estimation.predicted_covariance and reconciled_covariance both return the right shape.

• The covariance is propagated in Joseph form, \(P^+ = (I-K)P^-(I-K)^T + KRK^T\). The short form \((I-K)P^-\) is algebraically equal and numerically worse: it loses symmetry over a long run and can go indefinite. A twin is a long run.

• A weakly informative period needs no special case. Its \(R\) is large, \(K\) goes to zero, the estimate does not move.

• The time update runs on held periods too. Holding is not knowing: while the twin refuses to move a parameter, its error bar widens at exactly the rate the drift model claims, and the next permitted update takes a correspondingly larger step. max_std bounds that growth if a quiet year would otherwise make the next step a jump.

Innovation.nis is the filter’s own global test — \(v^T S^{-1} v\), which is \(\chi^2\) on \(p\) degrees of freedom when the drift model and the sigmas are right. A series running well above \(p\) says the parameter is moving faster than drift_std admits.

The gate

update_gate reads a MonitorDiagnosis and permits an update only on model drift:

verdict	gate	why
model drift	update	persistent rejection, blame wanders — the signature of a model fault
instrument fault	hold	blame is concentrated; go calibrate the named sensor
consistent	hold	the data gave the model nothing to correct
undiagnosed	hold	persistent rejection with nothing testable to blame

Holding on consistent makes the loop event-triggered: the parameter sits still until the evidence is strong enough to reject, then moves. That deadband is deliberate. A parameter re-estimated every period tracks whatever that period’s noise favoured, and the twin stops being a model.

The policy is an argument (allow=), so it can be widened — deliberately, and at a known cost. tests/test_tracking.py::TestGate::test_an_ungated_loop_would_have_invented_a_fouling_factor is the regression: the same filter on the same biased-meter data, with the gate opened, reports a pipe several percent off clean.

The gate is a statistical rule, not a guarantee, and the rule has two knobs. track_parameters forwards window, rejection_threshold and concentration_threshold to diagnose, because the defaults are not right for every plant: on a small network a few days early in a sensor bias can read as diffuse, slip through as model drift, and move the parameter before the rule settles. When that happens the fix is to tune the rule to the plant’s own noise — lengthen the window, lower the concentration threshold — never to widen allow, which removes the rule instead of sharpening it. Section 7 of examples/29_model_updating.ipynb shows both the leak and the tuning on a real gas network.

Where the parameter lives

The tracked parameters are threaded through params, so the frozen and free problems are the same residual_fn — a twin whose two clocks run different code drifts apart in a second, less interesting way. parameter_measurement appends them to the state vector with sigma = inf and hands the augmented problem to reconcile unchanged.

They are left free there rather than having the prior passed in as a finite sigma, on two counts. sigma is a vector, so a prior smuggled through it would be diagonal and would discard exactly the correlations the filter keeps. And with the prior outside, the reconciliation’s objective stays a test of data against model, uninflated by how confident the twin already was.

The cost is that each period must identify the parameters on its own. If it cannot, the structure check raises ReconciliationStructureError naming them — loudly, rather than returning a NaN. Pool several periods with reconcile_multi and hand its shared estimate and covariance straight to measurement_update instead.

What it does not fix

A filter on parameters assumes the model form is right. Under structural mismatch the parameter converges to something that is not the physical quantity and that shifts with operating point, so the twin’s gradients are wrong even where its values match — which is fatal, since everything downstream of a differentiable flowsheet is a derivative. The symptom is visible in run.monitor.statistic: the gate opens, the parameter moves, and the statistic does not come back down, because no value of the parameter fits. The answer then is difflow.planning.modifiers.update_modifiers, which corrects values and gradients against the plant, not a faster filter.

track_parameters is the offline driver, and it is also a backtest — run it over a recorded campaign to choose drift_std before trusting the loop live. The four steps are public and stateless (update_gate, time_update, parameter_measurement, measurement_update), so an online loop is the same four calls with a TrackerState carried between them.

Section 7 of examples/29_model_updating.ipynb runs the whole loop on the gas network the notebook builds, over both campaigns. The gated loop tracks the fouling pipe to 1.257 against a truth of 1.300 and never moves at all on the biased-meter campaign; wiring the gate open buys 0.03 of that lag back and reports 1.165 — a confident 16% fouling claim — on the pipe that is clean.