Parameter Identification

Simcoon provides a Python-based parameter identification workflow using scipy.optimize and a generic key-based file templating system. This replaces the former built-in C++ genetic algorithm and allows users to choose any optimizer, cost function, or external simulation tool.

Overview

The identification workflow consists of three components:

  1. Key system (Parameter / Constant): generic file templating that replaces placeholders with parameter values in any input file

  2. Forward model: simcoon solver, L_eff homogenization, or any external tool (e.g., fedoo for FEMU)

  3. Optimizer: sim.identification() wraps scipy.optimize.differential_evolution, or call scipy directly for more control

Template files (keys/)          Working files (data/)
┌──────────────────┐    copy    ┌──────────────────┐
│  ... @2p ...     │ ────────>  │  ... @2p ...     │
│  ... @3p ...     │            │  ... @3p ...     │
└──────────────────┘            └──────────────────┘
                                       │ apply
                                       v
                                ┌──────────────────┐
                                │  ... 73000 ...   │  ──> Forward model
                                │  ... 0.22 ...    │      (solver, L_eff, ...)
                                └──────────────────┘

Key System

The key system decouples the optimizer from the simulation tool. Template files contain alphanumeric placeholders (keys) that are replaced at each iteration with the current parameter values.

Parameter class

from simcoon.parameter import Parameter, read_parameters

# Create parameters programmatically
params = [
    Parameter(number=0, bounds=(100, 300), key="@E",
              sim_input_files=["material.json"]),
    Parameter(number=1, bounds=(0.1, 0.4), key="@nu",
              sim_input_files=["material.json"]),
]

# Or read from a file
params = read_parameters("data/parameters.inp")

Parameter attributes:

  • number: parameter index

  • bounds: (min, max) tuple — used as optimizer bounds

  • key: placeholder string in template files (e.g., @E, @0p)

  • sim_input_files: list of files containing this key

  • value: current value (defaults to midpoint of bounds)

Parameters file format (parameters.inp):

#Number  #min     #max     #key  #number_of_files  #files
0        100      300      @E    1                  material.json
1        0.1      0.4      @nu   1                  material.json

A template is any text file, so the JSON inputs of the solver work as they are: keys/material.json with "props": [70000.0, 0.3, 1.0e-5, @sigmaY, @k, @m] is not valid JSON until the substitution and is valid after it. The identification/json_keys_identification.py example runs this chain end to end (KeyTemplate, load_simulation_json, solve, calc_cost, identification).

Constant class

Constants are fixed values (not optimized) that also use the key system:

from simcoon.constant import Constant, read_constants

The Constant class is a NamedTuple with fields: number, key, input_values, value, sim_input_files.

A constant holds one value per test in input_values, e.g. the fibre angle of each specimen of a composite identification. read_constants(n, ...) reads n values per constant from constants.inp:

#Number  #key     #test_0  #test_1  #test_2  #number_of_files  #files
0        @theta   0        45       90       1                 ply.json

value is the first one; KeyTemplate.write(dst, test=i) (below) substitutes the value of test i. A constant with a single value is the same in every test: it is not a constant of the identification but a fixed parameter of the simulation, to write in the template rather than as a key, and write(..., test=i) warns about it.

Writing the working files

A KeyTemplate reads the template files of keys/ once and keeps them in memory. Each write replaces the keys of the parameters, then those of the constants, by their current values and writes every working file once:

from simcoon.parameter import KeyTemplate

template = KeyTemplate("keys", params, consts)   # consts may be omitted

def cost(x):
    for p, v in zip(params, x):   # values set by the optimizer
        p.value = v
    template.write("data")
    ...

With several tests, each one is written in its own directory with the constants of that test:

for i in range(3):                       # one fibre angle per specimen
    template.write(f"data/test_{i}", test=i)

This works with any file format — simcoon input files, FE meshes, JSON configs, etc. The key system is deliberately simple: it performs string replacement, making it compatible with any simulation tool.

identification() — Global Optimization

sim.identification() wraps scipy.optimize.differential_evolution using the bounds from your Parameter objects. After optimization, the identified values are written back to each Parameter.value.

from simcoon.identify import identification
from simcoon.parameter import Parameter

params = [
    Parameter(0, bounds=(10000, 200000), key="@Ef",
              sim_input_files=["job.inp"]),
    Parameter(1, bounds=(0.01, 0.45), key="@nuf",
              sim_input_files=["job.inp"]),
]

result = identification(my_cost_function, params, seed=42, disp=True)
print(f"E_f = {params[0].value:.0f}, nu_f = {params[1].value:.3f}")

job.inp here is a placeholder: a template input deck of an external, file-driven solver, supplied by the user in the keys folder (it is not part of simcoon). sim_input_files is the bridge to such solvers: the keys are substituted into their decks before each evaluation (KeyTemplate.write, called from my_cost_function). Simcoon’s own models no longer need it — pass props, and for mean-field models the phases, straight to solve(), and leave sim_input_files out.

Arguments:

  • cost_fn: callable f(x) -> float where x is a parameter array

  • parameters: list of Parameter (bounds used for search space)

  • **kwargs: forwarded to differential_evolution (maxiter, popsize, tol, seed, polish=True, disp, etc.)

Returns: scipy.optimize.OptimizeResult

Bounded least squares

Most identification costs are sums of squared differences between a model and experiments. method="least_squares" minimizes them with the bounded trust-region reflective algorithm of scipy.optimize.least_squares, which exploits that structure and typically needs 7 to 20 times fewer model evaluations than differential evolution. The cost function then returns the vector of residuals; calc_residuals builds it with the arguments of calc_cost (its sum of squares is the "mse", "nmse" or "nmse_per_response" cost):

from simcoon.identify import identification, calc_residuals

def residuals(x):
    y_num = run_model(x)
    return calc_residuals(y_exp, y_num, metric="nmse_per_response")

result = identification(residuals, params, method="least_squares",
                        n_starts=3, seed=0)
cost = np.sum(result.fun ** 2)   # the calc_cost value at the solution

  • The parameters are mapped onto [0, 1] (scale="auto": logarithmically when the lower bound is positive, linearly otherwise; or "linear"/"log", per parameter if needed). The Jacobian is computed by finite differences in these coordinates (diff_step, relative step, default sqrt(eps)).

  • The first start is Parameter.value where set, the middle of the scaled range otherwise; n_starts - 1 more are drawn by Latin hypercube sampling (seed), and the best result is kept (result.starts lists them all).

  • A model failure is best returned as NaN residuals: the trial step is rejected and the trust region shrunk.

  • Other least_squares options (xtol, ftol, gtol, max_nfev, loss…) are forwarded.

Parallel evaluations

n_threads (default 1, serial) evaluates the cost function in a pool of threads: the finite-difference columns of the Jacobian for "least_squares" (the result does not depend on the number of threads), a whole generation for differential evolution (updating="deferred"). It pays off because solve() and L_eff release the GIL, and requires a thread-safe cost function.

For more control (other optimizers, constraints, custom initialization), call scipy.optimize directly — the Parameter objects provide the bounds and values you need.

calc_cost() — Multi-Level Weighted Cost

sim.calc_cost() computes a weighted cost function with three levels of weights, designed for multi-test identification:

\[C = \text{avg}\left( W^{\text{test}}_i \cdot W^{\text{resp}}_{i,k} \cdot W^{\text{pt}}_{i,k,j} \cdot (y^{\exp}_{i,k,j} - y^{\num}_{i,k,j})^2 \right)\]

Data is organized as a list of 2-D arrays, one per test, each of shape (n_points, n_responses):

from simcoon.identify import calc_cost
import numpy as np

# Two tensile tests, each with force + displacement columns
y_exp = [
    np.column_stack([force_exp_1, disp_exp_1]),  # test 1: (N, 2)
    np.column_stack([force_exp_2, disp_exp_2]),  # test 2: (N, 2)
]
y_num = [
    np.column_stack([force_num_1, disp_num_1]),
    np.column_stack([force_num_2, disp_num_2]),
]

# Simple MSE
cost = calc_cost(y_exp, y_num)

# NMSE per response (balances force in N vs disp in mm)
cost = calc_cost(y_exp, y_num, metric='nmse_per_response')

# Per-test weights (emphasize test 2)
cost = calc_cost(y_exp, y_num, w_test=np.array([1.0, 3.0]))

Weight levels

Three levels of weights are combined multiplicatively:

Level

Argument

Description

Test

w_test

ndarray (n_tests,) — weight per experiment/file. Use to emphasize certain tests over others.

Response

w_response

list of ndarray (n_responses,) — weight per response column. To balance responses of different magnitudes, use metric='nmse_per_response' instead of manual weights.

Point

w_point

list of ndarray (n_points, n_responses) — weight per data point. Use for heterogeneous confidence or to mask outliers.

Metrics

Built-in metrics (numpy only, no extra dependency):

  • "mse" — Mean Squared Error (default)

  • "nmse" — Normalized MSE (divided by variance of all experimental data)

  • "nmse_per_response" — NMSE computed independently per response column, then averaged. Each column is divided by its own sum(y_exp^2), summed over the tests, balancing responses of different magnitudes (e.g., force in N vs displacement in mm). This is the recommended metric for multi-response identification. The tests are pooled within a column, not normalised one by one: a test with larger values weighs more (use w_test to rebalance).

  • "rmse" — Root Mean Squared Error

  • "mae" — Mean Absolute Error

  • "mape" — Mean Absolute Percentage Error (scikit-learn definition)

  • "wmape" — weighted MAPE, sum(w |y_exp - y_num|) / sum(w |y_exp|) (robust to responses crossing zero)

With scikit-learn installed (pip install simcoon[identify]):

  • "r2" — R-squared score

  • "mean_squared_error", "mean_absolute_error" — sklearn wrappers

  • Any sklearn.metrics function that accepts sample_weight

If scikit-learn is not installed and an sklearn metric is requested, a clear error message with install instructions is shown.

Using with External Solvers

The key system works with any simulation tool. For example, with fedoo (Finite Element Model Updating — FEMU):

import subprocess
from simcoon.parameter import KeyTemplate, Parameter
from simcoon.identify import identification, calc_cost

params = [
    Parameter(0, bounds=(100e3, 300e3), key="@E",
              sim_input_files=["material.json"]),
]
template = KeyTemplate("keys", params)

def cost(x):
    params[0].value = x[0]
    template.write("data")

    # Run external solver
    subprocess.run(["python", "run_fedoo_simulation.py"])

    # Load results and compare
    y_num = [np.loadtxt("results/reaction_force.txt")]
    y_exp = [np.loadtxt("exp_data/reaction_force.txt")]
    return calc_cost(y_exp, y_num, metric="nmse_per_response")

result = identification(cost, params, seed=42)