A new family of Constitutive Artificial Neural Networks towards automated model discovery

In plain words

A plain neural network is good at fitting a curve, but it knows nothing about physics. Fitted to how a rubber-like material stretches and pulls back, it can produce a stress curve that breaks basic rules: it may store negative energy, or flatten out where further stretching should make the material stiffer. It also tends to fit a few test points perfectly and then go wrong beyond them.

This paper builds a small network whose structure forces the learned energy to follow those rules. The inputs are two standard measures of how much the material is stretched, and the output is the stored energy. The stress is then computed from that energy by a standard thermodynamic relation, rather than learned directly. Because each part of the network matches a known textbook term, its few weights can be read as material constants.

Trained on rubber test data from three loading modes, the network selected a small set of terms that matched known models, and with data from several modes it gave one unique result. A classical network fits the same data well inside the tested range but does not extend sensibly beyond it.

TL;DR

• Constitutive Artificial Neural Networks (CANNs) model an isotropic, incompressible rubber as a free energy of two invariants; stress is derived from it, so the stated thermodynamic and physical constraints hold by design §5.
• Each of the eight network terms uses a custom activation; after training most of the twelve effective weights are zero §5 §7.
• Trained on uniaxial, equibiaxial and pure-shear rubber data, the CANN recovers one model and parameter set for each temperature §7 eq. 48 eq. 49.
• The paper reports no numeric fit error; fit quality and extrapolation are judged from graphs §7.
• Scope is a single material class and one proof-of-concept network; compressible, anisotropic and history-dependent cases are future work §8.

Problem

Classical neural networks approximate the stress as a function of the deformation gradient. They interpolate the training data well, but they do not predict beyond the training regime and they overfit sparse data §1 §4. They also do not enforce thermodynamic consistency, so the learned stress may violate the second law §1 §3.

Physics-informed networks add physics equations as loss terms §1. This paper instead changes the network input, output, architecture and activation functions so that the constraints hold by construction §1 §5. The scope is hyperelastic materials; the isotropic, perfectly incompressible case is the worked example §5 §8.

Method

Kinematics and invariants

The deformation map gives the deformation gradient F, with Jacobian J = det(F) greater than 0, and the right and left Cauchy-Green tensors C = Fᵗ·F and b = F·Fᵗ §2 eq. 1 eq. 2. For an isotropic material the principal invariants are I1, I2 and I3 = J², defined in terms of F or of C and b §2 eq. 3 eq. 4. For nearly incompressible materials, an isochoric split F = J^(1/3)·F̄ gives Ī1, Ī2 and J §2 eq. 5 eq. 6 eq. 7 eq. 8. A transverse-isotropy invariant I4 = λn² is defined with a preferred direction n §2 eq. 9. The fits in §7 use only the isotropic, incompressible case.

Constraints and the design choice each one drives

• Thermodynamic consistency: the Piola stress is the derivative of a free energy ψ with respect to F, P = ∂ψ/∂F §3 eq. 10 eq. 11. The network therefore outputs ψ, not P §3.
• Material objectivity: ψ depends on F only through C §3 eq. 12. The network input is therefore built from invariants, not from the nine components of F §3.
• Material symmetry (isotropy): ψ is a function of I1, I2, I3 §3 eq. 13. The network inputs are the shifted invariants [I1 − 3], [I2 − 3], [I3 − 1], so the reference state gives zero §3.
• Incompressibility: with I3 = 1, the input reduces to I1 and I2 §3 eq. 15.
• Physical restrictions: ψ ≥ 0 eq. 16; ψ = 0 with zero stress in the reference state eq. 17; ψ grows without bound as J → 0 or J → ∞ eq. 18. Activations must therefore be monotone, continuous, smooth at zero, zero at zero, and unbounded §5 fig. 4.
• Polyconvexity: ψ is written as a sum of separate functions of I1, I2 and I3 eq. 19, so the invariants are not mixed in the network. This implies non-negative weights §3 §5.

Isotropic incompressible CANN (the network used in §7)

The network has two inputs (I1 − 3 and I2 − 3) and one output, ψ §5 fig. 5. Each input passes through a first layer that forms its first and second powers, giving four nodes fig. 5 caption. A second layer then applies to each power either the identity or the term exp(α·) − 1, giving eight nodes §5 fig. 5. The output ψ is the sum of the eight node values, so no mixed I1·I2 products arise §5 eq. 24. Each invariant's four terms are a linear term, a linear exponential, a quadratic term and a quadratic exponential §5 fig. 4.

The explicit free energy is written in eq. 24 (order of the weights: the I1 terms use w₁,₁…w₁,₄ and w₂,₁…w₂,₄; the I2 terms use indices 5 to 8). One of each pair of weights is redundant, which leaves twelve effective parameters §5 after eq. 24. The stress is obtained from the derivative of ψ §5 eq. 25 eq. 35.

# Free energy derivatives of the isotropic incompressible CANN  [§5] [eq. 25]
# Input:  dI1 (B,), dI2 (B,) where dI = I - 3 ; w1, w2 (8,) each, >= 0
# Output: dpsi_dI1 (B,), dpsi_dI2 (B,)

def cann_derivatives(dI1, dI2, w1, w2):
    # I1 terms: linear, exponential, quadratic, quadratic exponential  [eq. 24]
    dpsi_dI1 = (
        w2[0] * w1[0]
        + w2[1] * w1[1] * exp(w1[1] * dI1)
        + 2.0 * dI1 * (w2[2] * w1[2] + w2[3] * w1[3] * exp(w1[3] * dI1**2))
    )
    # I2 terms use indices 4 to 7 in the same pattern  [eq. 24] [eq. 25]
    dpsi_dI2 = (
        w2[4] * w1[4]
        + w2[5] * w1[5] * exp(w1[5] * dI2)
        + 2.0 * dI2 * (w2[6] * w1[6] + w2[7] * w1[7] * exp(w1[7] * dI2**2))
    )
    return dpsi_dI1, dpsi_dI2

Uniaxial, equibiaxial and pure-shear stress

For each homogeneous mode the paper gives F, I1, I2, the derivatives with respect to the stretch λ, and the pressure from the zero transverse stress §6 fig. 6. The nominal stress for each mode is written in closed form in terms of ∂ψ/∂I1 and ∂ψ/∂I2:

• Uniaxial tension: F = diag(λ, λ^(−1/2), λ^(−1/2)) eq. 36; I1 = λ² + 2/λ and I2 = 2λ + 1/λ² eq. 37; P1 = 2[∂ψ/∂I1 + (1/λ)∂ψ/∂I2](λ − 1/λ²) eq. 39.
• Equibiaxial tension: F = diag(λ, λ, λ^(−2)) eq. 40; I1 = 2λ² + 1/λ⁴ and I2 = λ⁴ + 2/λ² eq. 41; P1 = P2 = 2[∂ψ/∂I1 + λ²∂ψ/∂I2](λ − 1/λ⁵) eq. 43.
• Pure shear: F = diag(λ, 1, λ^(−1)) eq. 44; I1 = I2 = λ² + 1 + 1/λ² eq. 45; P1 = 2[∂ψ/∂I1 + ∂ψ/∂I2](λ − 1/λ³) and P2 = 2[∂ψ/∂I1 + λ²∂ψ/∂I2](1 − 1/λ²) eq. 47.

# Nominal stress in uniaxial tension, isotropic incompressible CANN  [eq. 39]
# Input:  stretch (B,) ; w1, w2 (8,) each, >= 0
# Output: P1 (B,)

def cann_stress_uniaxial(stretch, w1, w2):
    # Invariants of uniaxial tension  [eq. 37]
    I1 = stretch**2 + 2.0 / stretch
    I2 = 2.0 * stretch + 1.0 / stretch**2
    # Shifted invariants used as network input  [§5]
    dI1 = I1 - 3.0
    dI2 = I2 - 3.0
    # Derivatives of psi with respect to I1 and I2  [eq. 25]
    dpsi_dI1, dpsi_dI2 = cann_derivatives(dI1, dI2, w1, w2)
    # Closed-form nominal stress  [eq. 39]
    return 2.0 * (dpsi_dI1 + dpsi_dI2 / stretch) * (stretch - 1.0 / stretch**2)

Training

The loss is the mean squared error between model stress and data over the training points §5 eq. 33. The paper uses the ADAM optimiser and constrains all network weights to remain non-negative §5. The paper reports the loss falling by about four orders of magnitude during training (see Results) §7.

# Fit the CANN to uniaxial stress data by gradient descent  [§5] [eq. 33]
# Input:  stretch (N,), stress_data (N,), initial w1, w2 (8,) each
# Output: fitted w1, w2 (8,) each, all >= 0

def fit_cann(stretch, stress_data, w1, w2):
    for step in range(num_steps):
        # Model stress at every training stretch  [eq. 39]
        stress_model = cann_stress_uniaxial(stretch, w1, w2)
        # Mean squared error over the training points  [eq. 33]
        loss = mean((stress_model - stress_data) ** 2)
        # ADAM update of all weights; learning rate not stated  [§5]
        w1, w2 = adam_step(loss, w1, w2)
        # Keep every weight non-negative  [§5]
        w1 = maximum(w1, 0.0)
        w2 = maximum(w2, 0.0)
    return w1, w2

Special cases

The weight combinations that recover known models are given in eq. 26 to eq. 30:

• neo-Hooke: μ = 2·w₁,₁·w₂,₁ eq. 26.
• Blatz–Ko (incompressible form): μ = 2·w₁,₅·w₂,₅ eq. 27.
• Mooney–Rivlin: μ1 = 2·w₁,₁·w₂,₁ and μ2 = 2·w₁,₅·w₂,₅ eq. 28.
• Yeoh: a1 = 2·w₁,₁·w₂,₁, a2 = 2·w₁,₃·w₂,₃, a3 = 0 eq. 29.
• Demiray: a = 2·w₁,₂·w₂,₂ and b = w₁,₂ eq. 30.

The paper also writes a Holzapfel term in I4 eq. 32 and a bulk term for nearly incompressible materials eq. 31, neither of which is fitted in §7.

Hyperparameters

Results

Author-reported. Setup: rubber benchmark data tab. 1 tab. 2; eight-term isotropic incompressible CANN fig. 5; classical tanh networks as the baseline §7. The paper reports no numeric error metric; fit quality is judged from graphs §7.

Learned parameters for the simultaneous multi-mode fit (author-reported):

The printed 50° model writes μ2 with [I1−3]; see Pitfalls.

Paper vs code

Code: https://github.com/LivingMatterLab/CANN at commit a639ea97acfa1064005381645c1a2fda74bcbfcd, linked from the abstract. The repository has no LICENSE file. The CANN code is in RUBBER/; it is written in TensorFlow/Keras and the repo has no classical-network baseline code for §7 . Only the code for the isotropic CANN is checked here.

Checked and matching:

• The eight-term free energy and its I1 and I2 blocks match eq. 24: four terms per invariant, with the same exponential and square structure Models.py Models.py.
• The invariants for uniaxial, equibiaxial and pure shear match eq. 37 eq. 41 eq. 45 ContinuumMechHelper.py.
• The closed-form stresses match: uniaxial eq. 39 and [code: RUBBER/ContinuumMechHelper.py#L31-L41]; pure shear eq. 47 and [code: RUBBER/ContinuumMechHelper.py#L44-L54]; equibiaxial eq. 43 and [code: RUBBER/ContinuumMechHelper.py#L67].
• The optimiser is ADAM and the loss is MSE §5 Rubber_CANN.py.

Discrepancies and undocumented choices:

• Learning rate and epochs are not stated in the paper. The code uses Adam with learning rate 0.001 Rubber_CANN.py, batch size 8 Rubber_CANN.py, and 12000 epochs Rubber_CANN.py.
• Early stopping monitors the loss with patience 1000 and restores the best weights Rubber_CANN.py. The paper states no stopping rule.
• The code adds an L2 penalty of 0.001 to every kernel Models.py (applied at Models.py, [L89], [L94], [L97] and [L129]). The paper mentions no regulariser.
• Initialisation: linear-type weights use Glorot normal, and exponential-type weights are drawn uniformly from [0, 1e-5] Models.py Models.py Models.py.
• The non-negativity constraint is a Keras NonNeg projection applied to every layer Models.py Models.py. The paper says only that weights stay non-negative §5.
• Multi-mode data are resampled: the PS and ET curves are interpolated with np.interp onto a linspace of the length of the UT data Rubber_CANN.py. The paper does not describe this.
• Data split by row index: the first 15 rows are UT, the next 15 are ET, and the rest are PS Rubber_CANN.py. The 20° UT column of tab. 2 has 25 rows, so the code's input is not the printed table as-is (inferred from code and table; the data file is missing, see below).
• Stress scaling: the loaded stress is multiplied by 10 with no comment Rubber_CANN.py. The paper says stresses were converted to MPa tab. 1.
• The training data file input/Treloar.pkl is not in the repository, so the training script cannot run as shipped Rubber_CANN.py.
• By default the script does not train. It loads the saved weights from Results/ Rubber_CANN.py Rubber_CANN.py.
• Multi-mode loss (inferred from code, not stated in the paper): one MSE is applied to three outputs, and Keras sums the outputs by default Rubber_CANN.py Rubber_CANN.py.
• Restarts: if the first prediction is NaN the fit is repeated, up to four attempts Rubber_CANN.py.
• Shipped checkpoint for the multi-mode run [code: RUBBER/Results/Rubber_CANN/all/All/Config_file.txt] reports R² of 0.9986 (UT), 0.9960 (PS) and 0.9826 (ET). The paper prints no R². The stored weights give μ1 = 2·w₁,₁·w₂,₁ ≈ 2.99 in the code's units, which does not match μ1 = 0.2370 MPa in eq. 48. The paper does not say which temperature or data set produced the shipped checkpoint, and the repo has no separate 50° run.
• Weights: the code keeps all sixteen raw values, while the paper reports twelve effective parameters after removing redundancy §5.

Pitfalls

• Missing training settings. The paper names no learning rate, batch size, epoch count, initialisation scheme or stopping rule. The code's values are listed above; an implementation that uses different settings may reach a different local minimum (the paper itself reports non-unique fits for single-mode data fig. 12).
• Non-negativity and sparsity. The paper's "most weights train to zero" §7 §8 is reached in the code by clipping at zero and by the L2 penalty, which means exact zeros in the stored weights [code: RUBBER/Results/Rubber_CANN/all/All/Config_file.txt]. Implement the clipping explicitly, or the sparsity will not appear.
• Non-uniqueness needs multi-mode data. With one loading mode the split of stress between the I1 and I2 terms is non-unique fig. 12 §7. The paper attributes this to data richness, not to the network §7.
• Eq. 49 bracket. The 50° model writes μ2 with [I1−3]. The weight mapping in the text makes μ2 = 2·w₁,₅·w₂,₅, which belongs to the I2 linear term of eq. 24 (inferred; the printed bracket is probably a typo for [I2−3]). The code's index for w₁,₅ is the I2 linear term Models.py.
• Eq. 31 bracket. The prose says the bulk term is in [J−1], but eq. 31 shows [J−3]². This term is not used in the §7 fits.
• Extrapolation evidence. All listed data end at λ ≈ 7.7 tab. 1 tab. 2. The claims that the CANN "predicts" beyond the training range §7 rest on smoothness in graphs; no held-out data beyond that range are reported.
• Figure numbers. The printed figure numbers follow the order of appearance in the compiled paper, not the source file names. For example, the source file fig14 is printed as Figure 7 and fig07 as Figure 8 fig. 7 fig. 8.
• Numerics. The code computes the stress in float32 [code: RUBBER/ContinuumMechHelper.py#L35-L36]. Near λ = 1 the shifted invariants are close to zero, so rounding error is relative to small numbers (inferred from the formulas; not discussed in the paper).
• Scope. The worked example is isotropic and perfectly incompressible §5 §8. Other symmetry classes are described but not fitted.

Lineage

• Builds on: the classical rubber models named in §5 (neo-Hooke, Blatz–Ko, Mooney–Rivlin, Yeoh, Demiray, Holzapfel). None are in this wiki as arXiv entries. The earlier two-layer logistic network for rubber §1 is also not in this wiki.
• Compared against: classical tanh networks with one or two layers §4 §7.
• Superseded by: not known from wiki content. The repository README lists later CANN papers (Bayesian, Gaussian, generalised-invariant); none are in this wiki.