Optimize FEBio simulations

Purpose

The pyfebiopt.optimize module builds on XPLT reading. Optimization adapters usually open the .xplt generated by a FEBio run, slice the variables needed for one objective, and return arrays to the residual assembler.

Core classes

ExperimentSeries

Stores measured x and y arrays, plus optional weights. These are the reference data used to form residuals.

Parameter and ParameterSpace

Define fitted variables, initial values, bounds, and active/inactive flags. ParameterSpace maps physical parameters \(\theta\) to optimizer-space variables \(\phi\) when reparameterization is enabled.

FebTemplate, ParameterBinding, and EvaluationBinding

Connect Python parameters to XML nodes in a FEBio .feb file. ParameterBinding writes one parameter value directly. EvaluationBinding computes an arbitrary value from the full parameter dictionary, useful for coupled values such as mirrored fiber directions.

SimulationAdapter

Wraps a callable that reads one FEBio result file and returns simulated x and y arrays. Adapters are where .xplt slicing and project-specific post-processing normally live.

SimulationCase

Groups one template, one set of experiments, matching adapters, grid alignment policies, output folder names, and per-case FEBio runtime options.

EngineOptions

Collects configuration for runners, storage, cleanup, Jacobians, logging, monitoring, and optimizer selection.

Engine

Coordinates the full loop: render FEB files, launch FEBio, collect outputs, assemble residuals, log/report progress, and call the selected SciPy optimizer.

Assemble a script

Minimal scripts follow the same shape even when the objective is complex:

import numpy as np

from pyfebiopt.optimize import (
    CleanupOptions,
    Engine,
    EngineOptions,
    ExperimentSeries,
    FebTemplate,
    GridPolicyOptions,
    JacobianOptions,
    MonitorOptions,
    OptimizerOptions,
    ParameterBinding,
    ParameterSpace,
    RunnerOptions,
    SimulationAdapter,
    SimulationCase,
    StorageOptions,
)
from pyfebiopt.xplt import xplt

def read_force_curve(xplt_path):
    xp = xplt(str(xplt_path))
    xp.readAllStates()
    force = xp.results["reaction forces"].nodeset("grip").comp("z")
    time = np.asarray(xp._time, dtype=float)
    return time, np.asarray(force, dtype=float)

exp_data = np.loadtxt("data.txt")
experiment = ExperimentSeries(x=exp_data[:, 0], y=exp_data[:, 1])

parameters = ParameterSpace(xi=2.0)
parameters.add_parameter(name="G", theta0=0.6, bounds=(0.0, 10.0))

template = FebTemplate(
    "model_template.feb",
    bindings=[
        ParameterBinding(theta_name="G", xpath=".//Material/material[@id='1']/G"),
    ],
)

case = SimulationCase(
    template=template,
    subfolder="biaxial",
    experiments={"force": experiment},
    adapters={"force": SimulationAdapter(read_force_curve)},
    grids={"force": GridPolicyOptions(policy="exp_to_sim")},
    omp_threads=4,
)

options = EngineOptions(
    jacobian=JacobianOptions(enabled=True, perturbation=1e-3, parallel=True),
    cleanup=CleanupOptions(remove_previous=True, mode="all"),
    runner=RunnerOptions(jobs=2, command=("febio4", "-i")),
    storage=StorageOptions(mode="disk", root="runs"),
    monitor=MonitorOptions(enabled=True),
    optimizer=OptimizerOptions(
        name="least_squares",
        settings={"ftol": 1e-6, "xtol": 1e-6},
        reparametrize=True,
    ),
)

result = Engine(parameters, [case], options=options).run()
print(result.theta)

Objective

The default least-squares objective is assembled from all cases and experiments:

\[J(\theta) = \sum_{c \in \mathrm{cases}} \sum_{e \in \mathrm{experiments}} w_{c,e} \lVert f_{c,e}(\theta) - y_{c,e} \rVert_2^2\]

where \(f_{c,e}\) is produced by the adapter, \(y_{c,e}\) is the experimental series, and \(w_{c,e}\) is an optional weight.