.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/analysis/hyperelastic_parameter_identification.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_analysis_hyperelastic_parameter_identification.py: Parameter Identification for Mooney-Rivlin Hyperelastic Material ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ This example demonstrates inverse parameter identification for hyperelastic materials using experimental stress-strain data. We identify the Mooney-Rivlin model parameters from Treloar's classical rubber elasticity experiments. **Optimization**: scipy.optimize.differential_evolution (global optimizer) **Cost function**: Mean Squared Error on stress (sklearn.metrics.mean_squared_error) **Forward model**: simcoon hyperelastic stress functions .. GENERATED FROM PYTHON SOURCE LINES 13-22 .. code-block:: Python import numpy as np import pandas as pd import matplotlib.pyplot as plt from scipy.optimize import differential_evolution, fsolve from sklearn.metrics import mean_squared_error from simcoon import simmit as sim import os .. GENERATED FROM PYTHON SOURCE LINES 23-64 Methodology Overview -------------------- **Inverse Problem Formulation** Parameter identification is an inverse problem: given experimental observations (stress-strain data), we seek the material parameters that minimize the discrepancy between model predictions and experiments. The optimization problem is: .. math:: \boldsymbol{\theta}^* = \arg\min_{\boldsymbol{\theta}} \mathcal{L}(\boldsymbol{\theta}) where :math:`\boldsymbol{\theta}` are the material parameters and :math:`\mathcal{L}` is the cost function. **Cost Function: Stress-Based MSE** We use the Mean Squared Error (MSE) computed on the first Piola-Kirchhoff stress :math:`P_{11}`: .. math:: \text{MSE} = \frac{1}{N} \sum_{i=1}^{N} \left( P_{11}^{\text{model}}(\lambda_i; \boldsymbol{\theta}) - P_{11}^{\text{exp}}(\lambda_i) \right)^2 This stress-based metric directly measures the model's ability to predict mechanical response, which is the primary quantity of interest in constitutive modeling. **Why Differential Evolution?** We use ``scipy.optimize.differential_evolution`` because: 1. **Global optimization**: Unlike gradient-based methods, it explores the entire parameter space and avoids local minima 2. **Derivative-free**: No gradients required, robust for non-smooth objectives 3. **Bounded search**: Natural handling of physical parameter constraints 4. **Robustness**: Less sensitive to initial guess compared to local methods For hyperelastic models, the cost function landscape can be non-convex, making global optimization essential for reliable parameter identification. .. GENERATED FROM PYTHON SOURCE LINES 66-95 Mooney-Rivlin Constitutive Model -------------------------------- The Mooney-Rivlin model is a phenomenological hyperelastic model widely used for rubber-like materials. The strain energy function is: .. math:: W = C_{10}(\bar{I}_1 - 3) + C_{01}(\bar{I}_2 - 3) + U(J) where: - :math:`\bar{I}_1, \bar{I}_2` are the first and second isochoric invariants of the left Cauchy-Green tensor :math:`\mathbf{b} = \mathbf{F}\mathbf{F}^T` - :math:`C_{10}, C_{01}` are the material parameters to identify - :math:`U(J) = \kappa(J \ln J - J + 1)` is the volumetric contribution - :math:`\kappa` is the bulk modulus (fixed, assuming near-incompressibility) The Cauchy stress is computed as: .. math:: \boldsymbol{\sigma} = \boldsymbol{\sigma}_{\text{iso}} + \boldsymbol{\sigma}_{\text{vol}} where simcoon provides ``sigma_iso_hyper_invariants`` and ``sigma_vol_hyper`` to compute these contributions from the strain energy derivatives: .. math:: \frac{\partial W}{\partial \bar{I}_1} = C_{10}, \quad \frac{\partial W}{\partial \bar{I}_2} = C_{01}, \quad \frac{\partial U}{\partial J} = \kappa \ln J .. GENERATED FROM PYTHON SOURCE LINES 97-110 Experimental Data: Treloar's Rubber Experiments ----------------------------------------------- L.R.G. Treloar (1944) performed classical experiments on vulcanized rubber under three loading conditions: 1. **Uniaxial Tension (UT)**: :math:`\lambda_1 = \lambda, \lambda_2 = \lambda_3 = \lambda_t` 2. **Pure Shear (PS)**: :math:`\lambda_1 = \lambda, \lambda_2 = \lambda_t, \lambda_3 = 1` 3. **Equibiaxial Tension (ET)**: :math:`\lambda_1 = \lambda_2 = \lambda, \lambda_3 = \lambda_t` where :math:`\lambda_t` is the transverse stretch determined by equilibrium (zero transverse stress). The data provides stretch :math:`\lambda` and the corresponding first Piola-Kirchhoff stress :math:`P_{11}` [MPa]. .. GENERATED FROM PYTHON SOURCE LINES 110-141 .. code-block:: Python def load_treloar_data(filepath: str) -> pd.DataFrame: """ Load Treloar experimental data. The file contains columns for three loading cases: - lambda_1, P1_MPa: Uniaxial tension - lambda_2, P2_MPa: Pure shear - lambda_3, P3_MPa: Equibiaxial tension Parameters ---------- filepath : str Path to Treloar.txt data file Returns ------- pd.DataFrame Experimental data with columns for each loading case """ df = pd.read_csv( filepath, sep=r"\s+", engine="python", names=["lambda_1", "P1_MPa", "lambda_2", "P2_MPa", "lambda_3", "P3_MPa"], header=0, ) return df .. GENERATED FROM PYTHON SOURCE LINES 142-160 Forward Model: Stress Prediction using simcoon ---------------------------------------------- The forward model computes the first Piola-Kirchhoff stress for given material parameters and stretch values. This is the core function that will be called repeatedly during optimization. **Key steps:** 1. For each stretch :math:`\lambda_1`, solve for the transverse stretch :math:`\lambda_t` that satisfies equilibrium (zero transverse stress) 2. Construct the deformation gradient :math:`\mathbf{F}` 3. Compute the left Cauchy-Green tensor :math:`\mathbf{b} = \mathbf{F}\mathbf{F}^T` using ``sim.L_Cauchy_Green(F)`` 4. Compute strain energy derivatives for Mooney-Rivlin 5. Compute Cauchy stress using ``sim.sigma_iso_hyper_invariants`` and ``sim.sigma_vol_hyper`` 6. Convert to first Piola-Kirchhoff stress: :math:`P_{11} = \sigma_{11}/\lambda_1` .. GENERATED FROM PYTHON SOURCE LINES 160-285 .. code-block:: Python def mooney_rivlin_pk1_stress( C10: float, C01: float, kappa: float, lambda_array: np.ndarray, loading_case: str = "UT", ) -> np.ndarray: """ Compute first Piola-Kirchhoff stress for Mooney-Rivlin material. Uses simcoon's hyperelastic stress functions to compute the Cauchy stress from the strain energy derivatives, then converts to PK1 stress. Parameters ---------- C10 : float First Mooney-Rivlin parameter [MPa] C01 : float Second Mooney-Rivlin parameter [MPa] kappa : float Bulk modulus for volumetric response [MPa] lambda_array : np.ndarray Array of principal stretch values in the loading direction loading_case : str Type of loading: "UT" (uniaxial tension), "PS" (pure shear), or "ET" (equibiaxial tension) Returns ------- np.ndarray First Piola-Kirchhoff stress P_11 [MPa] """ PK1_stress = [] lambda_t_guess = 1.0 # Initial guess for transverse stretch for lam in lambda_array: # ----------------------------------------------------------------- # Step 1: Find transverse stretch satisfying equilibrium # ----------------------------------------------------------------- def equilibrium_residual(lambda_t_vec): """Residual function for finding equilibrium transverse stretch.""" lt = lambda_t_vec[0] # Construct deformation gradient based on loading case if loading_case == "UT": # Uniaxial: F = diag(lambda, lambda_t, lambda_t) F = np.diag([lam, lt, lt]) J = lam * lt**2 elif loading_case == "PS": # Pure shear: F = diag(lambda, lambda_t, 1) F = np.diag([lam, lt, 1.0]) J = lam * lt elif loading_case == "ET": # Equibiaxial: F = diag(lambda, lambda, lambda_t) F = np.diag([lam, lam, lt]) J = lam**2 * lt else: raise ValueError(f"Unknown loading case: {loading_case}") # Compute left Cauchy-Green tensor using simcoon b = sim.L_Cauchy_Green(F) # Mooney-Rivlin strain energy derivatives dW_dI1_bar = C10 dW_dI2_bar = C01 dU_dJ = kappa * np.log(J) # Compute Cauchy stress using simcoon functions sigma_iso = sim.sigma_iso_hyper_invariants( float(dW_dI1_bar), float(dW_dI2_bar), b, J, False ) sigma_vol = sim.sigma_vol_hyper(dU_dJ, b, J, False) sigma = sigma_iso + sigma_vol # Return stress component that should be zero at equilibrium if loading_case == "UT": # sigma_22 = sigma_33 = 0 return 0.5 * (sigma[1, 1] + sigma[2, 2]) elif loading_case == "PS": # sigma_22 = 0 return sigma[1, 1] else: # ET # sigma_33 = 0 return sigma[2, 2] # Solve for equilibrium transverse stretch lambda_t = fsolve(equilibrium_residual, lambda_t_guess, full_output=False)[0] lambda_t_guess = lambda_t # Use as starting point for next iteration # ----------------------------------------------------------------- # Step 2: Compute final stress state at equilibrium # ----------------------------------------------------------------- if loading_case == "UT": F = np.diag([lam, lambda_t, lambda_t]) J = lam * lambda_t**2 elif loading_case == "PS": F = np.diag([lam, lambda_t, 1.0]) J = lam * lambda_t else: # ET F = np.diag([lam, lam, lambda_t]) J = lam**2 * lambda_t b = sim.L_Cauchy_Green(F) # Strain energy derivatives dW_dI1_bar = C10 dW_dI2_bar = C01 dU_dJ = kappa * np.log(J) # Compute Cauchy stress sigma_iso = sim.sigma_iso_hyper_invariants( float(dW_dI1_bar), float(dW_dI2_bar), b, J, False ) sigma_vol = sim.sigma_vol_hyper(dU_dJ, b, J, False) sigma = sigma_iso + sigma_vol # Convert Cauchy stress to PK1: P = J * sigma * F^{-T} # For diagonal F: P_11 = sigma_11 / lambda_1 PK1_stress.append(sigma[0, 0] / lam) return np.array(PK1_stress) .. GENERATED FROM PYTHON SOURCE LINES 286-302 Cost Function: Stress-Based Mean Squared Error ---------------------------------------------- The cost function quantifies the discrepancy between model predictions and experimental data. We use MSE computed on the PK1 stress values. **Normalized MSE (NMSE)**: When combining multiple loading cases with different stress magnitudes, we can normalize by the variance to ensure balanced weighting: .. math:: \text{NMSE} = \frac{\text{MSE}}{\text{Var}(P^{\text{exp}})} This ensures that loading cases with smaller stress magnitudes are not dominated by cases with larger stresses during combined fitting. .. GENERATED FROM PYTHON SOURCE LINES 302-383 .. code-block:: Python def cost_function_mse( params: np.ndarray, lambda_exp: np.ndarray, P_exp: np.ndarray, kappa: float, loading_case: str, normalize: bool = False, ) -> float: """ Compute MSE cost between model prediction and experimental stress data. Parameters ---------- params : np.ndarray Material parameters [C10, C01] lambda_exp : np.ndarray Experimental stretch values P_exp : np.ndarray Experimental PK1 stress values [MPa] kappa : float Fixed bulk modulus [MPa] loading_case : str Loading type ("UT", "PS", or "ET") normalize : bool If True, divide MSE by variance of experimental data Returns ------- float MSE or NMSE value """ C10, C01 = params try: # Predict stress using forward model P_model = mooney_rivlin_pk1_stress(C10, C01, kappa, lambda_exp, loading_case) # Compute MSE using sklearn mse = mean_squared_error(P_exp, P_model) if normalize: variance = np.var(P_exp) return mse / variance if variance > 1e-12 else mse return mse except Exception: # Return large cost if computation fails (e.g., numerical issues) return 1e10 def cost_function_combined( params: np.ndarray, data_dict: dict, kappa: float, normalize: bool = True ) -> float: """ Combined cost function for simultaneous fitting to multiple loading cases. Parameters ---------- params : np.ndarray Material parameters [C10, C01] data_dict : dict Dictionary mapping loading case names to (lambda, P) tuples kappa : float Fixed bulk modulus [MPa] normalize : bool If True, use NMSE for balanced weighting across cases Returns ------- float Sum of (N)MSE values over all loading cases """ total_cost = 0.0 for case_name, (lambda_exp, P_exp) in data_dict.items(): cost = cost_function_mse(params, lambda_exp, P_exp, kappa, case_name, normalize) total_cost += cost return total_cost .. GENERATED FROM PYTHON SOURCE LINES 384-396 Parameter Identification using Differential Evolution ----------------------------------------------------- ``scipy.optimize.differential_evolution`` is a stochastic global optimizer based on evolutionary algorithms. Key parameters: - **bounds**: Search space for each parameter - **strategy**: Mutation strategy ('best1bin' works well for most problems) - **popsize**: Population size (larger = more thorough search) - **mutation**: Mutation constant (controls exploration vs exploitation) - **recombination**: Crossover probability - **polish**: If True, refines result with L-BFGS-B (local optimizer) .. GENERATED FROM PYTHON SOURCE LINES 396-557 .. code-block:: Python def identify_parameters( lambda_exp: np.ndarray, P_exp: np.ndarray, kappa: float = 4000.0, loading_case: str = "UT", bounds: list = None, normalize: bool = False, seed: int = 42, verbose: bool = True, ) -> dict: """ Identify Mooney-Rivlin parameters using differential evolution. Parameters ---------- lambda_exp : np.ndarray Experimental stretch values P_exp : np.ndarray Experimental PK1 stress values [MPa] kappa : float Fixed bulk modulus [MPa] (default: 4000 for near-incompressibility) loading_case : str Loading type ("UT", "PS", or "ET") bounds : list Parameter bounds [(C10_min, C10_max), (C01_min, C01_max)] normalize : bool If True, use normalized MSE seed : int Random seed for reproducibility verbose : bool Print optimization progress Returns ------- dict Optimized parameters and optimization results """ if bounds is None: bounds = [(0.01, 2.0), (-1.0, 1.0)] if verbose: print(f"\n{'=' * 60}") print(f"Optimizing for {loading_case} loading case") print(f"{'=' * 60}") print(f"Parameter bounds: C10 in {bounds[0]}, C01 in {bounds[1]}") print(f"Bulk modulus (fixed): kappa = {kappa} MPa") result = differential_evolution( cost_function_mse, bounds=bounds, args=(lambda_exp, P_exp, kappa, loading_case, normalize), strategy="best1bin", maxiter=500, popsize=15, tol=1e-8, mutation=(0.5, 1.0), recombination=0.7, seed=seed, polish=True, # Refine with L-BFGS-B disp=verbose, ) C10_opt, C01_opt = result.x if verbose: print(f"\nOptimization completed!") print(f" C10 = {C10_opt:.6f} MPa") print(f" C01 = {C01_opt:.6f} MPa") print(f" Final MSE = {result.fun:.6e} MPa^2") return { "C10": C10_opt, "C01": C01_opt, "kappa": kappa, "mse": result.fun, "success": result.success, "result": result, } def identify_parameters_combined( data_dict: dict, kappa: float = 4000.0, bounds: list = None, normalize: bool = True, seed: int = 42, verbose: bool = True, ) -> dict: """ Identify parameters by fitting to multiple loading cases simultaneously. Using NMSE (normalized MSE) ensures balanced contribution from each loading case, preventing cases with larger stress magnitudes from dominating the optimization. Parameters ---------- data_dict : dict Dictionary {case_name: (lambda_exp, P_exp)} for each loading case kappa : float Fixed bulk modulus [MPa] bounds : list Parameter bounds [(C10_min, C10_max), (C01_min, C01_max)] normalize : bool If True, use NMSE (recommended for combined fitting) seed : int Random seed verbose : bool Print progress Returns ------- dict Optimized parameters and results """ if bounds is None: bounds = [(0.01, 2.0), (-1.0, 1.0)] if verbose: print(f"\n{'=' * 60}") print(f"Combined optimization for: {list(data_dict.keys())}") print(f"{'=' * 60}") print(f"Parameter bounds: C10 in {bounds[0]}, C01 in {bounds[1]}") print(f"Bulk modulus (fixed): kappa = {kappa} MPa") print(f"Using {'NMSE' if normalize else 'MSE'} for cost function") result = differential_evolution( cost_function_combined, bounds=bounds, args=(data_dict, kappa, normalize), strategy="best1bin", maxiter=500, popsize=15, tol=1e-8, mutation=(0.5, 1.0), recombination=0.7, seed=seed, polish=True, disp=verbose, ) C10_opt, C01_opt = result.x if verbose: print(f"\nOptimization completed!") print(f" C10 = {C10_opt:.6f} MPa") print(f" C01 = {C01_opt:.6f} MPa") print(f" Final combined cost = {result.fun:.6e}") return { "C10": C10_opt, "C01": C01_opt, "kappa": kappa, "cost": result.fun, "success": result.success, "result": result, } .. GENERATED FROM PYTHON SOURCE LINES 558-560 Visualization Functions ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 560-632 .. code-block:: Python def plot_fit_comparison( lambda_exp: np.ndarray, P_exp: np.ndarray, params: dict, loading_case: str, title: str = None, ax=None, ): """ Plot experimental data vs model prediction. Parameters ---------- lambda_exp : np.ndarray Experimental stretch values P_exp : np.ndarray Experimental PK1 stress values [MPa] params : dict Dictionary with 'C10', 'C01', 'kappa' keys loading_case : str Loading type title : str Plot title (optional) ax : matplotlib.axes.Axes Axes to plot on (creates new figure if None) Returns ------- matplotlib.axes.Axes The axes with the plot """ if ax is None: fig, ax = plt.subplots(figsize=(8, 6)) # Model prediction on fine grid lambda_model = np.linspace(1.0, lambda_exp.max() * 1.02, 200) P_model = mooney_rivlin_pk1_stress( params["C10"], params["C01"], params["kappa"], lambda_model, loading_case ) # Experimental data ax.plot( lambda_exp, P_exp, "o", markersize=8, markerfacecolor="red", markeredgecolor="black", label="Treloar experimental", ) # Model prediction ax.plot( lambda_model, P_model, "-", linewidth=2, color="blue", label=f"Mooney-Rivlin (C10={params['C10']:.4f}, C01={params['C01']:.4f})", ) ax.set_xlabel(r"Stretch $\lambda$ [-]", fontsize=12) ax.set_ylabel(r"PK1 Stress $P_{11}$ [MPa]", fontsize=12) ax.set_title(title or f"{loading_case} - Parameter Identification", fontsize=14) ax.legend(loc="upper left", fontsize=10) ax.grid(True, alpha=0.3) return ax .. GENERATED FROM PYTHON SOURCE LINES 633-644 Main Execution -------------- We demonstrate two identification strategies: 1. **Individual fitting**: Optimize parameters separately for each loading case 2. **Combined fitting**: Find a single parameter set that best fits all cases The individual fits provide loading-case-specific parameters (useful for understanding material behavior), while combined fitting gives a universal parameter set for general use. .. GENERATED FROM PYTHON SOURCE LINES 644-829 .. code-block:: Python if __name__ == "__main__": # Locate data file relative to script # Try multiple approaches to find the data file (supports both direct # execution and Sphinx-Gallery which does not define __file__) data_filename = os.path.join("hyperelasticity", "comparison", "Treloar.txt") possible_paths = [] # Try using __file__ if available try: script_dir = os.path.dirname(os.path.abspath(__file__)) possible_paths.append(os.path.join(script_dir, "..", data_filename)) except NameError: pass # Try relative to current working directory (Sphinx-Gallery context) possible_paths.append(os.path.join(os.getcwd(), "..", data_filename)) possible_paths.append(os.path.join(os.getcwd(), "examples", data_filename)) # Find first existing path data_path = None for path in possible_paths: if os.path.exists(path): data_path = path break if data_path is None: raise FileNotFoundError( f"Could not find Treloar.txt data file. Searched: {possible_paths}" ) print("=" * 70) print(" MOONEY-RIVLIN PARAMETER IDENTIFICATION") print(" Using scipy.differential_evolution and sklearn.mean_squared_error") print("=" * 70) # ------------------------------------------------------------------------- # Load and prepare experimental data # ------------------------------------------------------------------------- df = load_treloar_data(data_path) print(f"\nLoaded Treloar data: {len(df)} data points") # Extract data for each loading case (remove NaN values) ut_mask = ~df["lambda_1"].isna() & ~df["P1_MPa"].isna() ps_mask = ~df["lambda_2"].isna() & ~df["P2_MPa"].isna() et_mask = ~df["lambda_3"].isna() & ~df["P3_MPa"].isna() lambda_ut = df.loc[ut_mask, "lambda_1"].values P_ut = df.loc[ut_mask, "P1_MPa"].values lambda_ps = df.loc[ps_mask, "lambda_2"].values P_ps = df.loc[ps_mask, "P2_MPa"].values lambda_et = df.loc[et_mask, "lambda_3"].values P_et = df.loc[et_mask, "P3_MPa"].values print(f"\nData summary:") print( f" Uniaxial Tension (UT): {len(lambda_ut)} points, " f"lambda in [{lambda_ut.min():.2f}, {lambda_ut.max():.2f}]" ) print( f" Pure Shear (PS): {len(lambda_ps)} points, " f"lambda in [{lambda_ps.min():.2f}, {lambda_ps.max():.2f}]" ) print( f" Equibiaxial (ET): {len(lambda_et)} points, " f"lambda in [{lambda_et.min():.2f}, {lambda_et.max():.2f}]" ) # Fixed bulk modulus (near-incompressible rubber) kappa = 4000.0 # ------------------------------------------------------------------------- # Individual fitting for each loading case # ------------------------------------------------------------------------- print("\n" + "=" * 70) print(" INDIVIDUAL FITTING (one parameter set per loading case)") print("=" * 70) params_ut = identify_parameters(lambda_ut, P_ut, kappa, "UT", verbose=True) params_ps = identify_parameters(lambda_ps, P_ps, kappa, "PS", verbose=True) params_et = identify_parameters(lambda_et, P_et, kappa, "ET", verbose=True) # Reference values from literature (Steinmann et al., 2012) print("\n" + "-" * 70) print("Comparison with literature (Steinmann et al., 2012):") print("-" * 70) print( f"{'Case':<8} {'C10 (lit.)':<12} {'C10 (ident.)':<14} {'C01 (lit.)':<12} {'C01 (ident.)':<14}" ) print("-" * 70) print( f"{'UT':<8} {0.2588:<12.4f} {params_ut['C10']:<14.4f} {-0.0449:<12.4f} {params_ut['C01']:<14.4f}" ) print( f"{'PS':<8} {0.2348:<12.4f} {params_ps['C10']:<14.4f} {-0.065:<12.4f} {params_ps['C01']:<14.4f}" ) print( f"{'ET':<8} {0.1713:<12.4f} {params_et['C10']:<14.4f} {0.0047:<12.4f} {params_et['C01']:<14.4f}" ) # ------------------------------------------------------------------------- # Combined fitting # ------------------------------------------------------------------------- print("\n" + "=" * 70) print(" COMBINED FITTING (single parameter set for all loading cases)") print("=" * 70) data_combined = { "UT": (lambda_ut, P_ut), "PS": (lambda_ps, P_ps), "ET": (lambda_et, P_et), } params_combined = identify_parameters_combined( data_combined, kappa=kappa, normalize=True, verbose=True ) # ------------------------------------------------------------------------- # Visualization # ------------------------------------------------------------------------- print("\n" + "=" * 70) print(" GENERATING PLOTS") print("=" * 70) fig, axes = plt.subplots(2, 3, figsize=(15, 10)) # Row 1: Individual fits plot_fit_comparison( lambda_ut, P_ut, params_ut, "UT", title="UT - Individual fit", ax=axes[0, 0] ) plot_fit_comparison( lambda_ps, P_ps, params_ps, "PS", title="PS - Individual fit", ax=axes[0, 1] ) plot_fit_comparison( lambda_et, P_et, params_et, "ET", title="ET - Individual fit", ax=axes[0, 2] ) # Row 2: Combined fit applied to all cases plot_fit_comparison( lambda_ut, P_ut, params_combined, "UT", title="UT - Combined fit", ax=axes[1, 0] ) plot_fit_comparison( lambda_ps, P_ps, params_combined, "PS", title="PS - Combined fit", ax=axes[1, 1] ) plot_fit_comparison( lambda_et, P_et, params_combined, "ET", title="ET - Combined fit", ax=axes[1, 2] ) fig.suptitle( "Mooney-Rivlin Parameter Identification using Treloar Data\n" "(Top: Individual fits, Bottom: Combined fit)", fontsize=14, fontweight="bold", ) plt.tight_layout() # ------------------------------------------------------------------------- # Summary # ------------------------------------------------------------------------- print("\n" + "=" * 70) print(" SUMMARY OF IDENTIFIED PARAMETERS") print("=" * 70) print(f"{'Case':<12} {'C10 [MPa]':>12} {'C01 [MPa]':>12} {'MSE [MPa^2]':>14}") print("-" * 52) print( f"{'UT':<12} {params_ut['C10']:>12.4f} {params_ut['C01']:>12.4f} {params_ut['mse']:>14.2e}" ) print( f"{'PS':<12} {params_ps['C10']:>12.4f} {params_ps['C01']:>12.4f} {params_ps['mse']:>14.2e}" ) print( f"{'ET':<12} {params_et['C10']:>12.4f} {params_et['C01']:>12.4f} {params_et['mse']:>14.2e}" ) print("-" * 52) print( f"{'Combined':<12} {params_combined['C10']:>12.4f} {params_combined['C01']:>12.4f} {params_combined['cost']:>14.2e}" ) print("=" * 70) print("\nNote: The combined fit uses NMSE (normalized MSE) which explains") print("the different cost magnitude compared to individual MSE values.") plt.show() .. image-sg:: /examples/analysis/images/sphx_glr_hyperelastic_parameter_identification_001.png :alt: Mooney-Rivlin Parameter Identification using Treloar Data (Top: Individual fits, Bottom: Combined fit), UT - Individual fit, PS - Individual fit, ET - Individual fit, UT - Combined fit, PS - Combined fit, ET - Combined fit :srcset: /examples/analysis/images/sphx_glr_hyperelastic_parameter_identification_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none ====================================================================== MOONEY-RIVLIN PARAMETER IDENTIFICATION Using scipy.differential_evolution and sklearn.mean_squared_error ====================================================================== Loaded Treloar data: 25 data points Data summary: Uniaxial Tension (UT): 25 points, lambda in [1.00, 7.61] Pure Shear (PS): 14 points, lambda in [1.00, 4.96] Equibiaxial (ET): 17 points, lambda in [1.00, 4.44] ====================================================================== INDIVIDUAL FITTING (one parameter set per loading case) ====================================================================== ============================================================ Optimizing for UT loading case ============================================================ Parameter bounds: C10 in (0.01, 2.0), C01 in (-1.0, 1.0) Bulk modulus (fixed): kappa = 4000.0 MPa differential_evolution step 1: f(x)= 0.4096615268469581 differential_evolution step 2: f(x)= 0.4096615268469581 differential_evolution step 3: f(x)= 0.3978263758640262 differential_evolution step 4: f(x)= 0.3978263758640262 differential_evolution step 5: f(x)= 0.39194853823940684 differential_evolution step 6: f(x)= 0.39194853823940684 differential_evolution step 7: f(x)= 0.390742001678335 differential_evolution step 8: f(x)= 0.390742001678335 differential_evolution step 9: f(x)= 0.38955535356791926 differential_evolution step 10: f(x)= 0.38955535356791926 differential_evolution step 11: f(x)= 0.38955535356791926 differential_evolution step 12: f(x)= 0.3895493503783034 differential_evolution step 13: f(x)= 0.3895134000272561 differential_evolution step 14: f(x)= 0.3895134000272561 differential_evolution step 15: f(x)= 0.3895134000272561 differential_evolution step 16: f(x)= 0.38951161523788896 differential_evolution step 17: f(x)= 0.38951044296197623 differential_evolution step 18: f(x)= 0.3895083336274455 differential_evolution step 19: f(x)= 0.38950807193905634 differential_evolution step 20: f(x)= 0.3895079700689561 differential_evolution step 21: f(x)= 0.3895079700689561 differential_evolution step 22: f(x)= 0.3895078172607434 differential_evolution step 23: f(x)= 0.3895078172607434 differential_evolution step 24: f(x)= 0.3895078172607434 differential_evolution step 25: f(x)= 0.3895078094426372 differential_evolution step 26: f(x)= 0.3895078094426372 differential_evolution step 27: f(x)= 0.3895078012254065 differential_evolution step 28: f(x)= 0.38950780048910943 differential_evolution step 29: f(x)= 0.38950780048910943 differential_evolution step 30: f(x)= 0.38950780044294975 differential_evolution step 31: f(x)= 0.38950780038681637 differential_evolution step 32: f(x)= 0.38950780038681637 Polishing solution with 'L-BFGS-B' Optimization completed! C10 = 0.407989 MPa C01 = -0.751039 MPa Final MSE = 3.895078e-01 MPa^2 ============================================================ Optimizing for PS loading case ============================================================ Parameter bounds: C10 in (0.01, 2.0), C01 in (-1.0, 1.0) Bulk modulus (fixed): kappa = 4000.0 MPa differential_evolution step 1: f(x)= 0.052568489490637596 differential_evolution step 2: f(x)= 0.009493123290881412 differential_evolution step 3: f(x)= 0.009493123290881412 differential_evolution step 4: f(x)= 0.009493123290881412 differential_evolution step 5: f(x)= 0.002238709247196296 differential_evolution step 6: f(x)= 0.002222638903038 differential_evolution step 7: f(x)= 0.002093401662362508 differential_evolution step 8: f(x)= 0.002051757060569482 differential_evolution step 9: f(x)= 0.002051757060569482 differential_evolution step 10: f(x)= 0.0020515205718457806 differential_evolution step 11: f(x)= 0.0020486426587081443 differential_evolution step 12: f(x)= 0.0020486426587081443 differential_evolution step 13: f(x)= 0.002048074478582873 differential_evolution step 14: f(x)= 0.0020480272988842034 differential_evolution step 15: f(x)= 0.0020480272988842034 differential_evolution step 16: f(x)= 0.0020480272988842034 differential_evolution step 17: f(x)= 0.0020480272988842034 differential_evolution step 18: f(x)= 0.0020480272988842034 differential_evolution step 19: f(x)= 0.0020480250885571804 differential_evolution step 20: f(x)= 0.0020480166211757785 differential_evolution step 21: f(x)= 0.0020480166211757785 differential_evolution step 22: f(x)= 0.0020480166211757785 differential_evolution step 23: f(x)= 0.0020480166211757785 differential_evolution step 24: f(x)= 0.0020480166211757785 differential_evolution step 25: f(x)= 0.0020480166211757785 differential_evolution step 26: f(x)= 0.0020480165511174825 differential_evolution step 27: f(x)= 0.0020480162263757373 differential_evolution step 28: f(x)= 0.0020480162263757373 differential_evolution step 29: f(x)= 0.0020480162252322466 differential_evolution step 30: f(x)= 0.0020480162150828224 differential_evolution step 31: f(x)= 0.0020480162150828224 differential_evolution step 32: f(x)= 0.0020480162150828224 differential_evolution step 33: f(x)= 0.0020480162150828224 differential_evolution step 34: f(x)= 0.0020480162144333124 differential_evolution step 35: f(x)= 0.002048016214353775 Polishing solution with 'L-BFGS-B' Optimization completed! C10 = 0.465673 MPa C01 = -0.296036 MPa Final MSE = 2.048016e-03 MPa^2 ============================================================ Optimizing for ET loading case ============================================================ Parameter bounds: C10 in (0.01, 2.0), C01 in (-1.0, 1.0) Bulk modulus (fixed): kappa = 4000.0 MPa differential_evolution step 1: f(x)= 3.7890783827233254 differential_evolution step 2: f(x)= 1.5801552994420787 differential_evolution step 3: f(x)= 1.5801552994420787 differential_evolution step 4: f(x)= 1.5801552994420787 differential_evolution step 5: f(x)= 0.033449646241513344 differential_evolution step 6: f(x)= 0.033449646241513344 differential_evolution step 7: f(x)= 0.0038665433369498774 differential_evolution step 8: f(x)= 0.0038665433369498774 differential_evolution step 9: f(x)= 0.0038665433369498774 differential_evolution step 10: f(x)= 0.0038665433369498774 differential_evolution step 11: f(x)= 0.0038665433369498774 differential_evolution step 12: f(x)= 0.0038665433369498774 differential_evolution step 13: f(x)= 0.0038665433369498774 differential_evolution step 14: f(x)= 0.0033141935848357527 differential_evolution step 15: f(x)= 0.0028921861132308784 differential_evolution step 16: f(x)= 0.0028921861132308784 differential_evolution step 17: f(x)= 0.0028921861132308784 differential_evolution step 18: f(x)= 0.0028236479044486317 differential_evolution step 19: f(x)= 0.0028236479044486317 differential_evolution step 20: f(x)= 0.0028236479044486317 differential_evolution step 21: f(x)= 0.0028236479044486317 differential_evolution step 22: f(x)= 0.0028236479044486317 differential_evolution step 23: f(x)= 0.0028236479044486317 differential_evolution step 24: f(x)= 0.0028236098274586146 differential_evolution step 25: f(x)= 0.002823571015744447 differential_evolution step 26: f(x)= 0.002823571015744447 differential_evolution step 27: f(x)= 0.002823571015744447 differential_evolution step 28: f(x)= 0.002823571015744447 differential_evolution step 29: f(x)= 0.002823571015744447 differential_evolution step 30: f(x)= 0.0028235703402129573 differential_evolution step 31: f(x)= 0.002823569975441282 differential_evolution step 32: f(x)= 0.002823569890218678 differential_evolution step 33: f(x)= 0.0028235698721478884 differential_evolution step 34: f(x)= 0.00282356976826437 differential_evolution step 35: f(x)= 0.00282356976826437 differential_evolution step 36: f(x)= 0.00282356976826437 differential_evolution step 37: f(x)= 0.00282356976826437 differential_evolution step 38: f(x)= 0.002823569764235306 differential_evolution step 39: f(x)= 0.0028235697641588657 differential_evolution step 40: f(x)= 0.0028235697638824744 differential_evolution step 41: f(x)= 0.0028235697638824744 Polishing solution with 'L-BFGS-B' Optimization completed! C10 = 0.171238 MPa C01 = 0.004730 MPa Final MSE = 2.823570e-03 MPa^2 ---------------------------------------------------------------------- Comparison with literature (Steinmann et al., 2012): ---------------------------------------------------------------------- Case C10 (lit.) C10 (ident.) C01 (lit.) C01 (ident.) ---------------------------------------------------------------------- UT 0.2588 0.4080 -0.0449 -0.7510 PS 0.2348 0.4657 -0.0650 -0.2960 ET 0.1713 0.1712 0.0047 0.0047 ====================================================================== COMBINED FITTING (single parameter set for all loading cases) ====================================================================== ============================================================ Combined optimization for: ['UT', 'PS', 'ET'] ============================================================ Parameter bounds: C10 in (0.01, 2.0), C01 in (-1.0, 1.0) Bulk modulus (fixed): kappa = 4000.0 MPa Using NMSE for cost function differential_evolution step 1: f(x)= 35.21222741563679 differential_evolution step 2: f(x)= 9.537265331692677 differential_evolution step 3: f(x)= 0.9974706563134659 differential_evolution step 4: f(x)= 0.8786897719288405 differential_evolution step 5: f(x)= 0.53121067957046 differential_evolution step 6: f(x)= 0.4965791322994794 differential_evolution step 7: f(x)= 0.4965791322994794 differential_evolution step 8: f(x)= 0.46264130322439295 differential_evolution step 9: f(x)= 0.46264130322439295 differential_evolution step 10: f(x)= 0.46120121776724954 differential_evolution step 11: f(x)= 0.46056166365362206 differential_evolution step 12: f(x)= 0.46056166365362206 differential_evolution step 13: f(x)= 0.46056166365362206 differential_evolution step 14: f(x)= 0.46056166365362206 differential_evolution step 15: f(x)= 0.46056166365362206 differential_evolution step 16: f(x)= 0.46055920122048505 differential_evolution step 17: f(x)= 0.460500046967362 differential_evolution step 18: f(x)= 0.4604919190725185 differential_evolution step 19: f(x)= 0.4604333550735898 differential_evolution step 20: f(x)= 0.4604331492392568 differential_evolution step 21: f(x)= 0.4604298987818909 differential_evolution step 22: f(x)= 0.4604298987818909 differential_evolution step 23: f(x)= 0.46042949481992085 differential_evolution step 24: f(x)= 0.46042949481992085 differential_evolution step 25: f(x)= 0.4604292754019745 differential_evolution step 26: f(x)= 0.4604292754019745 differential_evolution step 27: f(x)= 0.4604292663257766 differential_evolution step 28: f(x)= 0.4604292663257766 differential_evolution step 29: f(x)= 0.4604292663257766 differential_evolution step 30: f(x)= 0.46042926572184356 differential_evolution step 31: f(x)= 0.46042926572184356 differential_evolution step 32: f(x)= 0.4604292652772896 Polishing solution with 'L-BFGS-B' Optimization completed! C10 = 0.192048 MPa C01 = 0.003199 MPa Final combined cost = 4.604293e-01 ====================================================================== GENERATING PLOTS ====================================================================== ====================================================================== SUMMARY OF IDENTIFIED PARAMETERS ====================================================================== Case C10 [MPa] C01 [MPa] MSE [MPa^2] ---------------------------------------------------- UT 0.4080 -0.7510 3.90e-01 PS 0.4657 -0.2960 2.05e-03 ET 0.1712 0.0047 2.82e-03 ---------------------------------------------------- Combined 0.1920 0.0032 4.60e-01 ====================================================================== Note: The combined fit uses NMSE (normalized MSE) which explains the different cost magnitude compared to individual MSE values. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 15.567 seconds) .. _sphx_glr_download_examples_analysis_hyperelastic_parameter_identification.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: hyperelastic_parameter_identification.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: hyperelastic_parameter_identification.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: hyperelastic_parameter_identification.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_