Sensitivity, uncertainty and calibration

Which parameters matter, how much a result spreads when they are uncertain, which values make the model match a measurement, and which value hits a design target: the workbench's Analysis mode and four diokle commands answer these on any model, with no scripting.

Only free parameters can be varied: public parameters with a literal value, the ones that are editable in the variable browser. Parameters computed from others, protected ones and those evaluated at translation give an error that lists the free ones.

#In the workbench

Simulate a model, switch to Analysis in the toolbar and tick the parameters to vary. Each tab has its own settings and Run; a progress bar shows while it runs and Cancel stops it.

The Calibrate tab loads a CSV, maps its columns to model variables and sets the bounds of each tuner. It plots the measured data against the model before and after the fit, lists the fitted values with standard errors and correlations, and Copy as modifier gives the result as a Modelica modification.

#The math in brief

Write the model as $\dot{x} = f(t, x, p)$, $y = h(x, p)$, with parameters $p$.

Sensitivities are the partial derivatives $S_{ij} = \partial y_i / \partial p_j$: how much an output moves when a parameter moves. The exact way is the forward sensitivity equation, which Diokle solves together with the model,

When the exact method fails, a central difference is used instead, $S_{ij} \approx \bigl(y_i(p + \delta e_j) - y_i(p - \delta e_j)\bigr) / (2\delta)$.

Adjoint gradients are for one scalar $J = \int_0^T g(x, p)\,dt$ and many parameters. One backward solve of the adjoint equations gives every $\mathrm{d}J/\mathrm{d}p_j$ at once,

The formula above is for parameters that enter the model equations; a parameter in the initial state adds one more term, $\lambda(0)^T\,\partial x_0/\partial p$. Its cost does not grow with the number of parameters, so it suits a model with hundreds of them.

Calibration minimizes the weighted squared error between measurements $d_k$ and the model,

by Levenberg–Marquardt steps $\bigl(J^T J + \mu I\bigr)\,\Delta p = -J^T r$, where $r$ holds the residuals and $J$ the sensitivities of the outputs. The --tune name=lo:hi bounds are the projection that keeps every step inside the box.

Design of experiments runs the model at chosen factor settings, then fits a quadratic response surface. With factors $x$ coded to $[-1, 1]$ (x = (p - center) / half-range),

and the reported $R^2$ tells how much of the variation of each quantity the surface explains. Three factors give 10 terms, as in the example below.

#Sensitivity

diokle sensitivity library/Diokle -m Diokle.Examples.QuarterCar.BumpComparison \
    --params linear.m_body,linear.c_spring,linear.c_tyre,bump.h --outputs linear.a_body

sensitivities of Diokle.Examples.QuarterCar.BumpComparison (central finite differences), 2001 points from t = 0 to 2
  note: finite differences used because the exact method failed: step size too small at t = 0.008507344142682349
  note: central differences, relative step 1e-4, tolerance 1e-8

ranking (time average of |p/y dy/dp|, summed over outputs):
   1. linear.m_body                   2.60659e0   (p = 400)
   2. linear.c_spring                 1.74835e0   (p = 30000)
   3. bump.h                          1.20881e0   (p = 0.05)
   4. linear.c_tyre                  5.00852e-1   (p = 250000)

--sens-method chooses how:

• exact integrates the forward sensitivity equations ds/dt = df/dx s + df/dp next to the state, with exact tangents of the compiled model, algebraic loops included. It needs a model without when clauses and delay(), and a run without events.
• fd uses two runs per parameter with a relative step of 1e-4 at a tightened tolerance.
• auto (the default) uses exact where it applies and otherwise says why it used fd, as above, where the hydraulic damper's fast start stops the exact method.

The output holds the time series of dy/dp and of the normalized p/y dy/dp, and -o writes them as JSON or CSV.

#Adjoint gradients

With many parameters and one scalar to differentiate, --adjoint computes the whole gradient for about the cost of one forward and one backward solve, however many parameters there are:

diokle sensitivity crates/diokle-core/tests/models/adjoint.mo -m Adjoint.Chain --adjoint \
    --objective 'integral(y)' --params k1,k2,k3,k4,k5,k6,g,w,e

gradient of integral(y) for Adjoint.Chain (adjoint: backward solve with reverse-mode tape passes)
  note: adjoint: 1000 backward Radau IIA steps over 1001 mesh nodes; 4000 reverse and 12000 tangent tape passes (forward 0.8 ms, backward 10.1 ms)

  J = 4.0528170837e0

  parameter                             value            dJ/dp    p/J dJ/dp
  k1                               1.000000e0    7.79789943e-2    1.9241e-2
  k5                              6.000000e-1    1.13501154e-1    1.6803e-2
  g                               9.000000e-1    7.33660425e-2    1.6292e-2
  k6                               1.100000e0   -5.09701197e-2   -1.3834e-2
  e                                1.500000e0    3.60191614e-2    1.3331e-2
  k3                              8.000000e-1    6.24062793e-2    1.2319e-2
  k4                               1.700000e0    2.13180931e-2    8.9421e-3
  k2                               1.300000e0    1.97269796e-2    6.3277e-3
  w                                2.000000e0   -5.16754637e-3   -2.5501e-3

Nine parameters, one backward solve of about 10 ms.

• Objectives: final(x), at(x,T), integral(x) or mean(x), with an optional window such as integral(x,T0,T1). max and min are differentiated at the maximum or minimum on the output grid.
• The backward pass solves -dλ/dt = (df/dx)ᵀ λ + (dg/dx)ᵀ with a two-stage Radau IIA step, and each f_pᵀ λ is one reverse pass over the compiled model, algebraic loops included.
• --grad-method auto|adjoint|forward|fd: auto uses the adjoint and falls back to forward sensitivities or differences for models with events, saying why.
• The gradients agree with forward sensitivities and central differences to about 1e-7 relative on the test models (crates/diokle-core/tests/adjoint.rs).

The same gradients drive optimize --algorithm lbfgsb and calibrate --optimizer lbfgsb (bound-constrained L-BFGS-B). On a 6-state chain with 9 parameters, L-BFGS-B converges in about 5 runs plus 5 gradients, where Nelder–Mead needs several hundred runs.

#Uncertainty and global sensitivity

montecarlo samples the parameters from their distributions and reports the mean, standard deviation and the 5, 50 and 95 % bands of each output; gsa adds Sobol indices (first-order and total, in N (k + 2) runs) or Morris screening (--morris). Runs use all cores.

diokle gsa library/Diokle -m Diokle.Examples.QuarterCar.BumpComparison \
    --dist linear.m_body=250:350 --dist linear.c_spring=20e3:40e3 --dist bump.h=normal:0.05:0.01 \
    --qoi 'max(linear.a_body)' -N 128

Sobol indices of Diokle.Examples.QuarterCar.BumpComparison: 640 runs, 0 failed

quantity                             mean          std           5%          50%          95%
max(linear.a_body)              1.64112e1    3.72167e0    1.04179e1    1.62511e1    2.29829e1

max(linear.a_body) (mean 1.64266e1, variance 1.45367e1):
  parameter                              S1               ST
  linear.m_body              0.1981 +-0.088    0.1619 +-0.046
  linear.c_spring            0.0079 +-0.028    0.0111 +-0.003
  bump.h                     0.6979 +-0.191    0.6082 +-0.129

The bump height explains about 70 % of the spread in peak body acceleration and the spring rate almost none. That took 14.5 s for 640 runs on four cores.

• Distributions: p=lo:hi (uniform), p=normal:MEAN:STD, p=lognormal:MEAN:STD; Latin hypercube or random sampling (--sampling), seeded (--seed).
• Quantities: final(x), max(x), min(x), mean(x), integral(x) and at(x,T), each with an optional window such as max(x,T0,T1).

#Calibration

calibrate fits parameters to measured data with a projected Levenberg–Marquardt method. The CSV has a time column and measured columns, matched to variables by name or with --map column=var. Here the measurement is a simulated run with m_body = 360 and c_spring = 36000, with Gaussian noise added, and the fit starts from the model's 400 and 30000. Make the file with a seeded noise step:

diokle simulate library/Diokle -m Diokle.Examples.QuarterCar.BumpComparison \
    -p linear.m_body=360 -p linear.c_spring=36000 --vars linear.a_body,linear.body.s -o meas0.csv
python3 - <<'EOF'
import csv, random
random.seed(1)
rows = list(csv.reader(open('meas0.csv')))
out = [rows[0]] + [[r[0]] + [str(float(x) + random.gauss(0, 0.02 * abs(float(x)) + 1e-3))
                             for x in r[1:]] for r in rows[1:]]
csv.writer(open('meas.csv', 'w', newline='')).writerows(out)
EOF

The fitted values depend on the seed:

diokle calibrate library/Diokle -m Diokle.Examples.QuarterCar.BumpComparison \
    --data meas.csv --tune linear.m_body=250:500 --tune linear.c_spring=15e3:50e3

calibration of Diokle.Examples.QuarterCar.BumpComparison: 4 iterations, 4052 residuals, sensitivities: finite differences
  converged: cost and step no longer decrease
  weighted residual RMS: 6.810141e-1 -> 6.676438e-1

  parameter                       initial         fitted   std. error
  linear.m_body               4.0000000e2    3.6174678e2      3.761e0
  linear.c_spring             3.0000000e4    3.5693377e4      9.340e2

  correlation:
    [  1.000   0.342]
    [  0.342   1.000]

  column                         weight     RMS before      RMS after
  linear.a_body               5.1789e-1    3.021453e-1    4.033338e-2
  linear.body.s                7.2494e1    1.310878e-2    1.302128e-2

modifier: (linear(m_body = 361.74678, c_spring = 35693.377))

Weights default to one over each column's standard deviation (--weight changes them), and --window T0:T1 restricts the fit in time.

#Optimization

optimize runs Nelder–Mead inside the bounds. The objective is the sum of --minimize q, --maximize q and --target q=v terms, with the same quantities as above. Which body mass brings the peak body acceleration over the bump down to 10 m/s²?

diokle optimize library/Diokle -m Diokle.Examples.QuarterCar.BumpComparison \
    --tune linear.m_body=250:600 --target 'max(linear.a_body)=10'

optimization of Diokle.Examples.QuarterCar.BumpComparison (Nelder-Mead): 29 iterations, 58 runs; simplex converged
  objective: 5.751403e0 -> 2.923810e-17
  max(linear.a_body): 1.239821e1 -> 1.000000e1

  parameter                       initial        optimal
  linear.m_body               4.0000000e2    5.0070870e2

modifier: (linear(m_body = 500.7087))

Dynamic optimization, where the inputs vary over time and constraints apply along the trajectory, is the job of the CasADi export.

#Design of experiments

diokle doe runs a designed set of experiments and fits a quadratic response surface to each quantity. Three factors on the quarter car, in a central composite design:

diokle doe library/Diokle -m Diokle.Examples.QuarterCar.BumpComparison \
    --factor linear.m_body=300:500 --factor linear.c_spring=20e3:40e3 --factor bump.h=0.03:0.07 \
    --design ccd --qoi 'max(linear.a_body)' --qoi 'max(linear.body.s)' --pareto min,min

central composite design of Diokle.Examples.QuarterCar.BumpComparison: 15 runs over linear.m_body, linear.c_spring, bump.h, 0 failed
  note: surface coefficients in coded units: x = (p - center) / half-range

max(linear.a_body): R^2 = 0.99959, adjusted R^2 = 0.99884, RMSE = 1.7867e-1 (15 runs, 10 terms)
  term                              coefficient    std error           t
  1                                  1.239930e1    9.6031e-2     1.291e2
  linear.m_body                     -3.124491e0    5.6500e-2    -5.530e1
  linear.c_spring                   4.985366e-1    5.6500e-2     8.824e0
  bump.h                             5.198314e0    5.6500e-2     9.201e1
  linear.m_body*bump.h              -1.249334e0    6.3169e-2    -1.978e1
  linear.m_body^2                   7.431336e-1    1.1142e-1     6.670e0
  ...

Pareto front of max(linear.a_body) (min) and max(linear.body.s) (min): 1 runs
  run    2:    5.77769e0   6.04319e-1   linear.m_body = 500.000000, linear.c_spring = 20000.000000, bump.h = 0.030000

In coded units, the bump height raises the peak acceleration by 5.2 m/s² over half its range and the body mass lowers it by 3.1 m/s². A heavier body also softens the effect of the bump (the m_body*bump.h interaction).

• Designs: full factorial with p=lo:hi:levels; frac, 2^(k−p) fractions with minimum-aberration generators, your own --generators or a --resolution; ccd (default), with --alpha face|rotatable|orthogonal; bbd, Box–Behnken; and lhs.
• Per quantity: R², adjusted R², RMSE and coefficients with standard errors and t-values. Terms the design cannot estimate are dropped and listed.
• Main effects, two-factor interaction tables and, with --pareto for two quantities, the non-dominated runs.

The workbench's DOE tab picks the design and the quantities and shows main effects, interactions, coefficients and the Pareto front. Use applies a Pareto point's parameters.

#Validated enclosures (prototype)

crates/diokle-core/src/validated.rs encloses a simulation: it returns an interval that must contain the exact end state for every initial state and parameter in a box, rather than one trajectory. Each step checks an a priori box that contains the solution over the step, and then computes the end state by one of two schemes, plain Euler or the mean-value form with rotated coordinates. A third path, affine arithmetic (enclose_affine), keeps linear dependencies between the states and gives the narrowest enclosure on the test oscillator: about 0.43 wide, against 2.10 for Euler and 25.6 for the mean-value form.

The rounding is outward at each operation, so the result is guaranteed under the usual floating point model, but the affine path's error term is a per-operation allowance rather than a full analysis, and the library is a prototype. It is not on the command line.

#Identifiability

crates/diokle-core/src/identifiability.rs takes the sensitivities of a simulated experiment and the measurement noise, and returns the Fisher information of the parameters, their standard errors (from its inverse, the Cramér-Rao bound) and the pairs that cannot be told apart. A pair whose sensitivity vectors are proportional, such as two parameters that only enter as a product, is reported directly as confounded, so the result does not depend on the inverse. The module is a library function for now; it is not on the command line or in the workbench yet.

#Instant surrogates

crates/diokle-core/src/surrogate.rs fits a quadratic response surface to the results of a sweep or a design of experiments, so a new point can be estimated without a simulation. The inputs are scaled to [-1, 1] over their sampled ranges. Every holdout_every-th point is held out, and the root-mean-square error on those points is reported with every prediction.

A prediction outside the sampled box is refused, and a caller can set a limit on the error with predict_checked: above it, the exact model should be run instead. Quadratic surfaces are inaccurate across strongly nonlinear designs, and the held-out error is what shows it. The error is one number for the whole surface, so it does not say where the surface is poor. The module is a library function for now; it is not on the command line yet.

#Fault injection and FMEA

crates/diokle-core/src/faults.rs runs each fault at each operating point and compares the result with the nominal run of that point. A fault is a set of parameter overrides, such as a pump stuck off or a leak that opened. The largest deviation of an output is measured against that output's own nominal range, so outputs in different units can be compared. The effect of a fault is Failed when the deviation reaches the failure threshold, Degraded when it reaches the degraded threshold, and None otherwise. A faulted run that does not finish is also Failed.

The result is a failure mode and effects table (FMEA): one row per fault, one column per operating point, with the deviation in each cell, and to_markdown() writes it as a table. Requirements are not part of this first version; the output deviation is the only criterion. The module is a library function for now; it is not on the command line yet.

#Checked against known answers

crates/diokle-core/tests/analysis.rs compares these analyses with analytic results: exponential decay, a mass-spring-damper, the Ishigami function's Sobol indices, Monte Carlo moments, and parameter recovery from synthetic data.