6.1 Optimization Input File¶

This chapter describes FEBio's parameter optimization module. This module tries to estimate model parameters by solving an inverse finite element problem, where the “solution” of the problem is known, and the problem's model parameters are sought. For instance, the solution can be an experimentally determined reaction force curve, and the unknown parameters are the material parameters that will recreate the reaction force curve by solving a forward FE problem. Another common application is determining the model parameters that achieve a desired load, displacement, or other goal.

In any case, the optimization module tries to minimize an objective function of the form,

\[ \varphi\left(\mathbf{a}\right)=\sum\limits^{n}_{i=1}\left[y_{i}-f\left(x_{i};\mathbf{a}\right)\right]^{2}. \]

Here, the \(\left(x_{i},y_{i}\right)\) are user-defined data pairs and \(f\left(x;\mathbf{a}\right)\)is the function that extracts the corresponding data from the model. The optimization module tries to find the model parameters a that minimize the function \(\varphi\). It does this by repeatedly evaluating the function \(f\), which will usually call FEBio to solve a forward FE problem.

The optimization module requires a separate input file that describes all the information needed for solving an optimization problem. This input file is described next.

Like all FEBio input files, the parameter optimization input file is an XML-formatted text file. The following sections are defined:

Task optional section defining the task that is executed to solve the FE problem.

Options optional section defining the optimization control parameters.

Objective defines the objective function that will be minimized.

Parameters defines the model parameters that are to be determined.

Constraints defines linear constraints for the model parameters (requires a solution method that supports linear constraints, e.g. constrained levmar).

In the following paragraphs each section will be explained in detail.

Task Section¶

The optional Task section can be used to define the task that is executed by FEBio when solving the FE model. By default, FEBio will solve the standard forward problem, but if needed, a custom task can be defined.

Options Section¶

This section defines the control parameters for the optimization. The options that can be defined will depend on the chosen optimization method, which is set with an attribute of the Options section. The following table shows the currently supported optimization algorithms.

type	Description
levmar	Levenberg-Marquardt
constrained levmar	Levenberg-Marquardt with linear constraints (1)
powell	Powell's method
scan	Parameter domain scanning (2)

Comments:

The “constrained levmar” uses the third-party levmar libraryhttp://users.ics.forth.gr/~lourakis/levmar/. For additional information please see the library's home page.
The “scan” method will iterate over the parameter range in increments defined by the fourth parameter value defined in the Parameter section. It returns the parameters that resulted in the lowest objective value of visited points, but it may not be the actual minimum.

Note that the Options section is optional. When omitted, the levmar method is used with default values for the control parameters.

The following parameters can be defined.

Parameter	Description	Default
obj_tol	Convergence tolerance for objective function. (1)	0.001
f_diff_scale	Forward difference scale factor (2)	0.001
log_level	Sets the amount of output that is generated on screen and in the logfile by the FEBio solver (3)	LOG_NEVER
tau	Step size scale factor (4)	1.0e-3
print_level	Sets the amount of output generated by the optimization module (5)	PRINT_ITERATIONS

For example,

<Options type="constrained levmar">
  <obj_tol>0.001</obj_tol>
  <f_diff_scale>0.001</f_diff_scale>
  <print_level>PRINT_ITERATIONS</print_level>
</Options>

Comments:

The objective function that is to be minimized is a function of the form:

\[ \varphi\left(\mathbf{a}\right)=\sum\limits^{n}_{i=1}\left[y_{i}-f\left(x_{i};\mathbf{a}\right)\right]^{2}. \]

Here, \(f\left(x;\mathbf{a}\right)\) is the function that describes the model, \(\mathbf{a}\) is a vector with the (unknown) model parameters and the \(\left(x_{i},y_{i}\right)\) are the user-defined data that approximates the ideal model.

The optimization method currently implemented requires the calculation of the gradient of the model function \(f\) with respect to the model parameters \(\mathbf{a}\). Since this gradient is not known, it will be approximated using forward differences. For example, the \(k\)-th component of the gradient is approximated as follows.

\[ \frac{\partial f}{\partial a_{k}}\approx\frac{1}{\delta a_{k}}\left[f\left(a_{1},\cdots,a_{k}+\delta a_{k},\cdots,a_{m}\right)-f\left(a_{1},\cdots,a_{k},\cdots,a_{m}\right)\right] \]

The value for \(\delta a_{k}\)is determined from the following formula.

\[ \delta a_{k}=\varepsilon\left(1+a_{k}\right) \]

where, \(\varepsilon\)is the forward difference scale factor which can be set by the user with the _fdiff_scale _option.

The _log_level _allows the user to control exactly how much output from the time iterations of each optimization iteration is written to the screen and logfile. The following values are allowed:

Value	Description
LOG_DEFAULT	Use default settings (may be omitted)
LOG_NEVER	Don't generate any output
LOG_FILE_ONLY	Save to the logfile only
LOG_SCREEN_ONLY	Only print to the screen
LOG_FILE_AND_SCREEN	Save to the logfile and print to the screen

The step size scale factor tau can be set as an input only for the Constrained Levenberg-Marquardt method. It scales the initial guess for the damping parameter.
The _print_level _sets how much output is generated during the optimizations. In particular, it controls if the objective function values are printed are not. The following values are allowed:

Value	Description
PRINT_ITERATIONS	Print minimal iterations info
PRINT_VERBOSE	Print a lot of information, incl. objective function values

Parameters Section¶

This section defines the material parameters that are to be determined. Each parameter is defined using the param element:

<Parameters>
  <param name="[name]">[guess], [min], [max], [scale]</param>
</Parameters>

The name attribute gives the name of the parameter that is to be determined. The following section describes the format to reference a model parameter. Each parameter takes four values: [guess] is the initial guess for this parameter, [min] and [max] are the minimum and maximum values respectively for this parameter, and [scale] is a representative scale (magnitude) for this parameter (although the precise interpretation of this parameter depends on the particular solver). This value is used to normalize the optimization parameter and improve convergence. If not specified by the user, it defaults to the initial guess.

For example, to optimize the material parameters of a neo-Hookean solid, named mat1, use the following syntax.

<Parameters>
  <param name="fem.material('mat1').E">1.0, 0.0, 5.0, 1.0</param>
  <param name="fem.material('mat1').v">0.1, 0.0, 0.5, 1.0</param>
</Parameters>

This example defines two material parameters. The first component of the name (here mat1) is the name of the material as defined in the model input file. The second component (here E and v) are the names of the material parameters. See Appendix B for more information on how to reference model parameters.

Objective Section¶

This section defines the functional model that is used to evaluate the objective function (i.e. the function \(f\) above). The particular functional model is defined via the type attribute, and the following table lists the available options.

Value	Description
data-fit	Fits data, usually specified as a function of time, to a model output parameter
target	Optimizes until one or more output model parameters achieve a target value.
element-data	The target data is defined as a tabulated list of element values.
node-data	The target data is defined as a tabulated list of nodal values.

These options will be described separately in the following sections.

The data-fit model¶

The data-fit model tries to fit predicted FE data to a user-defined data curve. The data curve usually represents some experimentally obtained result, e.g. reaction force as a function of time. The data-fit model needs a data source, i.e. a function for extracting the model data, and the data curve.

<Objective type="data-fit">
  <fnc type="[enter type]">
    ...
  </fnc>
  <data>
    <pt>0, 0</pt>
    ...
    <pt>1, 1</pt>
  </data>
</Objective>

The data source is defined by the type attribute of the fnc tag. The following table shows the supported values.

type	Description
parameter	The function will track the evolution of a model parameter.
filter_positive_only	The function tracks the evolution of a model parameter, but filters only the positive values.
filter_sum	The function is defined by summing nodal values.

parameter

In the case of parameter, a child element defines the parameter that will act as the data source for the objective function. (See Appendix B on more information on referencing model parameters.)

<fnc type="parameter">
  <param name="fem.rigidbody('rigid').Fz"/>
</fnc>

This defines the data function as a function of time. You can also change the ordinate of the function to create other functional dependencies. For example, the following creates a load-displacement curve for a rigid body.

<fnc type="parameter">
  <ordinate name="fem.rigidbody('rigid').position.x"/>
  <param name="fem.rigidbody('rigid').Fz"/>
</fnc>

As another example, the following creates a stress-strain curve for a single solid element, using the left Hencky (logarithmic) strain component along the z-axis (denoted by hz), and the Cauchy normal stress component along z-axis (denoted by sz) in element number 1,

<fnc type="parameter">
  <ordinate name="fem.element_data('hz',1)"/>
  <param name="fem.element_data('sz',1)"/>
</fnc>

filter_positive_only

This option requires another data source that is needed to extract model data (e.g. parameter), but only positive values are returned. If the model parameter's value is negative, zero will be returned by the function.

<fnc type="filter_positive_only">
  <source type="parameter">
    <param name="fem.rigidbody('rigid').Fz"/>
  </source>
</fnc>

filter_sum

This option will sum up nodal values and use that in the objective function. The summed value is always assumed to be a function of time. The following example illustrates how this option can be used to define an objective function from the sum of nodal reaction forces.

<fnc type="filter_sum">
  <node_data data="Rx" node_set="top surface"/>
</fnc>

The target model¶

The target model tries to optimize model parameters such that one or more output variables achieve a user-defined target value. This model requires a list of variables and target values for each variable.

<Objective type="target">
  <var name="[enter name]">[target value]</var>
  ...
</Objective>

The name refers to a model output parameter. The example below will optimize until the X-euler angle of the rigid body 'Material2' becomes 0.5.

<Objective type="target">
  <var name="fem.rigidbody('Material2').euler.y">0.5</var>
</Objective>

The element-data model¶

The element-data model uses a tabulated list of element values as the target for the optimization. In addition, it requires an output variable that will be matched against the element data.

<Objective type="element-data">
  <var type="[enter variable name]"/>
  <data>
    <elem id="1">0.123</elem>
    ...
  </data>
</Objective>

The variable can be any of the element output data listed in the logfile section (See ElementSet Section).

The node-data model¶

The node-data model uses a tabulated list of nodal values as the target for the optimization. In addition, it requires an output variable that will be matched against the nodal data.

<Objective type="node-data">
  <var type="[enter variable name]"/>
  <data>
    <node id="1">0.123</elem>
    ...
  </data>
</Objective>

The variable can be any of the node output data listed in the logfile section (See Node_Data Class).

Constraints Section¶

The Constrained Levenberg-Marquardt method allows linear constraints on the material parameters. If \(\mathbf{a}\) is the material parameter vector, then the linear constraint is of the form:

\[ c_{1}a_{1}+c_{2}a_{2}+...+c_{n}a_{n}+b=0\quad. \]

The coefficients \(c_{1},c_{2},...,c_{n},b\) are the inputs of the constraint tag. For example, if the linear constraint is \(2a_{1}-a_{2}+3=0\), then the Constraints section would be:

<Constraints>
  <constraint>2, -1, 3</constraint>
</Constraints>