Calibrating Models to Data

objective(prob, alg)

When considering the InverseProblem or JointInverseProblem in the context of calibration, we can define an objective function that associates a number to a vector of possible search space values. This number is a measure of how close is the given vector to a solution of the inverse problem. The objective for an inverse problem is defined as a sum of the individual contributions from each experiment.

\[\sum_{i=1}^{N} loss(i).\]

In the above loss(i) is the individual contribution of the ith experiment and this can be computed via cost_contribution and is determined by the loss passed in the corresponding experiment constructor.

Calling objective(prob, alg) will return a function of the search space values, so in order to use it, it is best to store it in a variable as follows:

cost = objective(prob, alg);

In order to obtain a particular value that corresponds to some search space values, one has to input a NamedTuple or a vector of pairs to the above mentioned function.

As an example, one would call

cost([a => 1.0, b => 3.0])

where the search space in the inverse problem was defined as containing the a and b parameters, i.e. [a => (0.0, 10.0), b => (-5.0, 5.0)].

The above function also accepts calibration results as input, so if we have a calibration result obtained from something like

r = calibrate(prob, alg)

then we can compute the cost value corresponding to the calibration result using cost(r), where cost is the above mentioned variable.

The cost function returned by objective can also be called with 0 arguments, i.e. cost(), in which case the cost will be computed using the initial state of the optimization problem, which uses the default values present in the model for all parameters and initial conditions.

Positional arguments

• prob: the InverseProblem or JointInverseProblem defining the objective
• alg: the calibration algorithm that formulates the objective function

cost_contribution(alg, experiment::AbstractExperiment, invprob, x)

Compute the contribution of the given experiment to the total cost associated with the inverse problem invprob using the given calibration algorithm (alg).

Positional arguments

• alg: the calibration algorithm that formulates the objective function.
• experiment: the experiment for which we compute the contribution.
• invprob: the InverseProblem defining the objective.
• x: the values of the tuned parameters or initial conditions.

calibrate(invprob, alg; adtype = nothing, progress = true)

Find the best parameters to solve the inverse problem invprob using the calibration algorithm given by alg.

Arguments

• invprob: the InverseProblem or JointInverseProblem to solve.
• alg: the calibration algorithm used for building the loss function. This can be JuliaHubJob for launching batch jobs in JuliaHub. In this case, the actual alg is wrapped inside JuliaHubJob object.

Keyword Arguments

• adtype: Automatic differentiation choice, see the Optimization.jl docs for details. nothing (the default) selects AutoForwardDiff() for gradient-based optimizers (wrapped in SecondOrder when the optimizer needs Hessians). AutoEnzyme(; mode = set_runtime_activity(Reverse)) is supported as a first-class reverse-mode alternative for SingleShooting; it requires an adjoint sensealg on the experiment (e.g. GaussAdjoint(autojacvec = EnzymeVJP())) and optimize_tunables = false on the InverseProblem.
• progress: Show the progress of the optimization (current loss value) in a progress bar. Defaults to true.

calibrate(res, alg; adtype = AutoForwardDiff(), progress = true)

Continue calibration from res which is obtained from a previous calibration using the algorithm given by alg.

Arguments

• res: object of type CalibrationResult which is obtained from a previous calibration of an InverseProblem or JointInverseProblem.
• alg: the calibration algorithm used for building the loss function. This can be JuliaHubJob for launching batch jobs in JuliaHub. In this case, the actual alg is wrapped inside JuliaHubJob object.

Keyword Arguments

• adtype: Automatic differentiation choice, see the Optimization.jl docs for details. Defaults to AutoForwardDiff().
• progress: Show the progress of the optimization (current loss value) in a progress bar. Defaults to true.

SingleShooting(;maxiters = nothing, maxtime = nothing, optimizer = MadNLPOptimizer(hessian_approximation = MadNLP.CompactLBFGS), kwargs...)

Single shooting is the simplest transcription method for building the loss function. In this case we solve all the model equation and then we compare the solutions against the data with the user given error functions.

The optimization method will get the cost function as the sum of the individual previously mentioned costs. To change the optimization method, you can use the optimizer keyword argument. The default optimization method is MadNLPOptimizer(hessian_approximation = MadNLP.CompactLBFGS) which can be used from OptimizationMadNLP. The maximum number of iterations for the optimizer is given by maxiters, which is a required keyword argument. Atleast one of maxiters or maxtime is required. Optimization stops when either one of the condition is met. Atleast one of maxiters or maxtime is required. If both of them are given, optimization stops after either of condition is first met.

Keyword Arguments

• maxiters: Maximum numbers of iterations when using the method with a time stepping optimizer.
• maxtime: Maximum time for the optimization process.
• optimizer: An Optimization.jl algorithm for the global optimization. Defaults to MadNLPOptimizer(hessian_approximation = MadNLP.CompactLBFGS).
• kwargs: Keyword arguments passed onto Optimization.solve.

MultipleShooting(;maxiters = nothing, maxtime = nothing, optimizer = MadNLPOptimizer(), trajectories, ensemblealg = EnsembleThreads(), continuitytype = ModelStatesPenalty, continuitylossfun = squaredl2loss, continuitylossweight = 100, initialization = DefaultSimulationInitialization(), kwargs...)

Multiple shooting is a calibration algorithm in which we split the solves of the model equations in several parts, as it can be easier to avoid local minima if we fit the individual segments and them match them up. The number of groups can be controlled via the trajectories keyword argument. The initial conditions between the segments are added to the optimization as extra hidden parameters that have to be fit. The continuity of the resulting solution ensemble is enforced using a continuity penalty between segments.

The continuity loss function can be changed via continuitylossfun, defaulting to squaredl2loss and the penalty used for the continuity is weighted via a the continuitylossweight, which defaults to 100. All segments can be initialized using initialization keyword argument. It defaults to DefaultSimulationInitialization(), i.e. the values are obtained by simulating the model passed in the Experiment with default parameters. For information on other initializations, refer Initialization Methods section in the manual.

The optimization method will get the cost function as the sum of the individual previously mentioned costs. To change the optimization method, you can use the optimizer keyword argument. The default optimization method is MadNLPOptimizer which can be used from OptimizationMadNLP. The maximum number of iterations for the optimizer is given by maxiters, which is a required keyword argument. Atleast one of maxiters or maxtime is required. If both of them are given, optimization stops after either of condition is first met.

Keyword Arguments

• maxiters: Maximum numbers of iterations when using the method with a time stepping optimizer.
• maxtime: Maximum time for the optimization process.
• optimizer: An Optimization.jl algorithm for the global optimization. Defaults to MadNLPOptimizer which can be used from OptimizationMadNLP.
• trajectories: Required. This is number of segments that the whole time span is split into.
• ensemblealg: Parallelization method for the ensemble problem. Defaults to EnsembleThreads().
• continuitytype: Type of continuity enforced. This can be ModelStatesPenalty which adds a term in the loss which corresponds to continuity loss given by continuitylossfun or it can be ConstraintBased where the continuity equations are added as constraints for the optimizer.
• continuitylossfun: Loss function to compute continuity penalty. This is only applicable if continuitytype is ModelStatesPenalty.
• continuitylossweight: Weight multiplied with the continuity loss term in the total loss. This is only applicable if continuitytype is ModelStatesPenalty.
• initialization: Initialization method of the segments. This is only applicable if continuitytype is ModelStatesPenalty.
• kwargs: Keyword arguments passed onto Optimization.solve.

DefaultSimulationInitialization()

This is used for initializing segments for doing calibration using MultipleShooting. The values for all the segments are computed by simulating the model passed in the Experiment with default parameters and using its ODESolution.

DataInitialization(; interpolation = CubicSpline)

This is used for initializing segments for doing calibration using MultipleShooting. It is used for those states where data is present and the values are obtained by fitting using interpolation. For other states, initial condition for that state is used for all the segments.

Keyword Arguments

• interpolation: Interpolation method from DataInterpolations which can be passed to fit data for obtaining initial values for the segments.

RandomInitialization()

This is used for initializing segments for doing calibration using MultipleShooting. The values for all the segments are computed as a random value from 0 to 1 multiplied by the initial condition for that state.

StochasticMultipleShooting(; trajectories, batch_size, block_size = 2,
    sampling = :stratified_pairs,
    inner = Adam(1e-3), inner_kwargs = (; save_best = false, epochs = 5),
    auglag_kwargs = (;),
    maxiters, maxtime = nothing,
    diagnostics = NoDiagnostics(),
    ensemblealg = nothing,
    initialization = DefaultSimulationInitialization(),
    limit_blas_threads = true,
    kwargs...)

Mini-batch variant of MultipleShooting. Splits the full trajectory into trajectories uniform segments, samples batch_size of them per objective evaluation, and enforces continuity at all trajectories - 1 junctions as equality constraints via OptimizationAuglag's method-of-multipliers outer loop. The inner loop is a stochastic optimizer (typically Adam) whose gradient, computed on the sampled batch, is an unbiased Monte-Carlo estimator of the full-batch gradient.

The outer optimizer is fixed to OptimizationAuglag.AugLag since it is currently the only Optimization.jl algorithm that handles equality constraints together with a stochastic, data-iterator-aware inner loop, both of which are required by this method's formulation.

When to use mini-batch

For long trajectories with many states or many parameters (e.g. neural ODEs), a full-batch multiple-shooting gradient can dominate wall-clock time. Sampling a small number of segments per Adam step, while leaving the continuity constraints full-batch, trades a modest Monte-Carlo variance for an order-of-magnitude reduction in per-step cost.

Keyword Arguments

• trajectories: Required. Total number of uniform segments N.
• batch_size: Required. Number of segments sampled per inner objective evaluation. Must be in [1, trajectories] and a multiple of block_size.
• block_size: Size of consecutive-segment blocks used for batching (M in the design notes). M=1 is pure i.i.d. sampling, M=2 (default) pairs each sampled segment with its successor — giving one fully-coordinated continuity term per block — and M=N recovers full-batch multiple shooting. For the stratified sampling modes, trajectories must be a multiple of block_size — a trailing partial block would never be sampled, so its segments' data would never be fitted.
• sampling: One of :stratified_pairs (default), :shuffle, :sequential. Controls how blocks are ordered across epochs.
• inner: Stochastic optimizer used inside each AugLag outer iteration. Defaults to Adam(1e-3). Forwarded as AugLag(; inner = ...).
• inner_kwargs: NamedTuple of keyword arguments forwarded to the inner optimizer's solve (and init) via AugLag's inner_kwargs. The default (; save_best = false, epochs = 5) runs five full passes over the mini-batches per outer AugLag iteration, which gives Adam time to approximately solve the unconstrained subproblem before the multipliers update. User-supplied inner_kwargs is merged over (; save_best = false), so omitting save_best keeps the safe default while user-supplied keys win. save_best = false matters because AugLag's outer ρ/λ update consumes the LAST inner iterate of each sub-problem; reverting θ to best-loss-seen against the (about-to-change) sub-objective sabotages the multiplier schedule. Drop epochs to 1 for cheaper outer iterations (at the cost of looser subproblem convergence) or raise it for tighter convergence guarantees. maxiters is also accepted if you want to set an absolute Adam-step budget instead. Note that one epoch is trajectories ÷ batch_size Adam steps for :shuffle sampling and (trajectories ÷ block_size) ÷ (batch_size ÷ block_size) for the stratified modes.
• auglag_kwargs: NamedTuple splatted into the AugLag constructor to tune the method-of-multipliers outer loop itself (ρ initial value, ρ-growth factor, outer/inner tolerance splits, etc.). See OptimizationAuglag.AugLag for the available knobs. Defaults to (;).
• maxiters: Required. Number of outer AugLag iterations.
• maxtime: Optional outer-loop time budget.
• ensemblealg: Ensemble method for per-segment solves. When left unset (default), picks EnsembleThreads() for batch_size > 1 and falls back to EnsembleSerial() at batch_size = 1 where thread-spawn overhead costs more than it saves. User-supplied ensemblealg is respected regardless of batch_size.
• initialization: How to seed per-segment initial conditions. See DefaultSimulationInitialization, DataInitialization, RandomInitialization.
• limit_blas_threads: If true (default), sets BLAS.set_num_threads(1) around each mini-batch and junction-constraint solve, restoring the prior setting afterwards. Disable if the model has large matmuls or where BLAS threading matters.
• kwargs...: Forwarded to Optimization.solve.

DataShooting(;maxiters = nothing, maxtime = nothing, optimizer = MadNLPOptimizer(), groupsize, ensemblealg = EnsembleThreads(), continuitylossfun = squaredl2loss, continuitylossweight = 100, kwargs...)

DataShooting is a variant of the multiple shooting method (see MultipleShooting for more details), where we use the data for the internal initial conditions between segments instead of having the optimizer find them. If the data is considered noisy, it might be better to let the optimizer find the initial conditions, as we would inevitably introduce errors if we use DataShooting. To specify the group size, which is the number of time points in each trajectory where data points on both ends are taken as initial conditions, we can use groupsize keyword argument. The continuity of the resulting solution ensemble is enforced using a continuity penalty between segments. The continuity losss function can be changed via continuitylossfun, defaulting to squaredl2loss and the penalty used for the continuity is weighted via a the continuitylossweight, which defaults to 100.

The optimization method will get the cost function as the sum of the individual previously mentioned costs. To change the optimization method, you can use the optimizer keyword argument. The default optimization method is MadNLPOptimizer which can be used from OptimizationMadNLP. The maximum number of iterations for the optimizer is given by maxiters, which is a required keyword argument. Atleast one of maxiters or maxtime is required. If both of them are given, optimization stops after either of condition is first met.

Keyword Arguments

• maxiters: Maximum numbers of iterations when using the method with a time stepping optimizer.
• maxtime: Maximum time for the optimization process.
• optimizer: An Optimization.jl algorithm for the global optimization. Defaults to Ipopt.Optimizer.
• groupsize: Required. Number of time points in each trajectory.
• ensemblealg: Parallelization method for the ensemble problem. Defaults to EnsembleThreads().
• continuitytype: Type of continuity enforced. This can be ModelStatesPenalty which adds a term in the loss which corresponds to continuity loss given by continuitylossfun or it can be ConstraintBased where the continuity equations are added as constraints for the optimizer.
• continuitylossfun: Loss function to compute continuity penalty. This is only applicable if continuitytype is ModelStatesPenalty.
• continuitylossweight: Weight multiplied with the continuity loss term in the total loss. This is only applicable if continuitytype is ModelStatesPenalty.
• initialization: Initialization method of the segments. This is only applicable if continuitytype is ModelStatesPenalty.
• kwargs: Keyword arguments passed onto Optimization.solve.

DiscreteFixedGainPEM(alpha)

This model transformation implements the prediction error method as a discrete callback that brings the solution closer to the data at the corresponding timepoints as the integrator steps through with a factor that depends on alpha, which is a scalar value from 0 to 1. If alpha is 0, then this is equivalent with not applying any correction, while alpha=1 means that we start each step from the data.

If there is data for the unknowns, it is directly used. If there is data for observed variables, corresponding equations which are linear with respect to the unknowns are used to estimate other unknowns in those equations for which there is no data. This is only the case if the system of equations are invertible with respect to those other unknowns.

It is also useful for calibrating unstable systems, as pure simulation methods cannot track the true solution over longer time intervals and diverge.