DOI: 10.18419/darus-5370
title: A Modular Framework for Non-Linear Optimization of Linear Viscoelastic Model Parameters using a Second-Order Fractional Approach

description: &lt;p&gt; This repository presents a reproducible workflow for fitting the &lt;strong&gt;Linear Second-order Fractional (LSF)&lt;/strong&gt; model to mastercurves of thermoplastic polymers. 
The dataset is associated with the publication &lt;em&gt;Fauser et al. (2025)&lt;/em&gt; and complements the dataset &lt;a href="https://doi.org/10.18419/DARUS-4318"&gt;Data Sets for a Fractional Moisture-Dependent Viscoelasticity Model for Thermoplastic Polymers&lt;/a&gt;.
The LSF model effectively captures the broad relaxation spectra typical of polymers, but its interdependent parameters require robust non-linear optimization. To ensure clarity and maintainability, the fitting workflow is organized into three modular components: &lt;/p&gt; &lt;ol&gt; &lt;li&gt;&lt;strong&gt;Viscoelastic Models Class&lt;/strong&gt; (&lt;code&gt;models.py&lt;/code&gt;)&lt;/li&gt; &lt;li&gt;&lt;strong&gt;Configuration File&lt;/strong&gt; (&lt;code&gt;configuration_file.py&lt;/code&gt;)&lt;/li&gt; &lt;li&gt;&lt;strong&gt;Model Fitting Notebook&lt;/strong&gt; (&lt;code&gt;main.ipynb&lt;/code&gt;)&lt;/li&gt; &lt;/ol&gt;

&lt;h2&gt;1. Viscoelastic Models Class (&lt;code&gt;models.py&lt;/code&gt;)&lt;/h2&gt; &lt;p&gt; This module serves as the core library of constitutive models. Each model is implemented as a self-contained, object-oriented class. &lt;/p&gt; &lt;h3&gt;Purpose&lt;/h3&gt; &lt;p&gt;Encapsulate mathematical definitions of viscoelastic models, enabling modular extensions and reuse.&lt;/p&gt; &lt;h3&gt;Key Component: The LSF Model&lt;/h3&gt; &lt;p&gt; The &lt;strong&gt;Second-Order Fractional (LSF)&lt;/strong&gt; model is defined by its complex modulus: 

&lt;/p&gt; &lt;pre&gt; E*(ω) = (Ξ₀ + (iω)^α₁ Ξ₁ + (iω)^α₂ Ξ₂) / (Λ₀ + (iω)^α₁ Λ₁ + (iω)^α₂ Λ₂) &lt;/pre&gt; &lt;p&gt; 

Parameters: 
&lt;/p&gt; &lt;ul&gt; &lt;li&gt;Ξ₀, Ξ₁, Ξ₂ → strain-related coefficients&lt;/li&gt; 
&lt;li&gt;Λ₀, Λ₁, Λ₂ → stress-related coefficients&lt;/li&gt; 
&lt;li&gt;α₁, α₂ → fractional exponents&lt;/li&gt; &lt;/ul&gt; 

&lt;h3&gt;Parent Class: &lt;code&gt;ViscoelasticityModel&lt;/code&gt;&lt;/h3&gt; &lt;p&gt; 
Defines methods for computing the real/imaginary parts of the complex modulus, as well as the storage modulus, loss modulus, and loss factor: 

&lt;/p&gt; &lt;pre&gt; class ViscoelasticityModel(torch.nn.Module): 
 """
 Parent class for viscoelasticity models 
 """
 def __init__(self): 
 super().__init__()
 
 def complexModulus(self, w): 
 raise Exception
 
 def forward(self, w):
 comp = self.complexModulus(w)
 return comp.imag / comp.real 

 def storageModulus(self,w):
 comp = self.complexModulus(w) 
 return comp.real

 def lossModulus(self,w): 
 comp = self.complexModulus(w)
 return comp.imag&lt;/pre&gt;

 &lt;h3&gt;Subclass: &lt;code&gt;SecondOrderFractionalDerivativeModel&lt;/code&gt;&lt;/h3&gt; &lt;p&gt; Implements the LSF model, with parameters represented as &lt;code&gt;torch.nn.Parameter&lt;/code&gt; for optimization. A scaling factor normalizes experimental and model data for stability. 


&lt;/p&gt; &lt;pre&gt; class SecondOrderFractionalDerivativeModel(ViscoelasticityModel):
 """
 Class for a second order derivative model 
 (G0 + d/dt^alpha G1 + d/dt^beta G2) stress 
 = (E0 + d/dt^alpha E1 + d/dt^beta E2) strain 
 """
 def __init__(self,E0,E1,E2,G0,G1,G2,alpha,beta):
 super().__init__()
 
 self.E0 = torch.nn.Parameter(E0 * torch.ones(1, dtype=torch.float64))
 self.E1 = torch.nn.Parameter(E1 * torch.ones(1, dtype=torch.float64))
 self.E2 = torch.nn.Parameter(E2 * torch.ones(1, dtype=torch.float64))

 # normalize G0 to 1, G0 is not fitted 
 self.G0 = 1 * torch.ones(1, dtype=torch.float64) 
 self.G1 = torch.nn.Parameter(G1 * torch.ones(1, dtype=torch.float64))
 self.G2 = torch.nn.Parameter(G2 * torch.ones(1, dtype=torch.float64))

 self.alpha1 = torch.nn.Parameter(alpha * torch.ones(1, dtype=torch.float64))
 self.alpha2 = torch.nn.Parameter(beta * torch.ones(1, dtype=torch.float64))


 def complexModulus(self, w): 

 w = w[:,None]
 iw = torch.view_as_complex(torch.cat((torch.zeros(w.size()),w),1))

 iw_alpha1 = iw ** self.alpha1
 iw_alpha2 = iw ** self.alpha2

 comp = ((self.E0 + self.E1 * iw_alpha1 + self.E2 * iw_alpha2)
 / (self.G0 + self.G1 * self.E1 * iw_alpha1 + self.G2 * self.E2 * iw_alpha2))
 return comp &lt;/pre&gt;

&lt;h2&gt;2. Configuration File (&lt;code&gt;configuration_file.py&lt;/code&gt;)&lt;/h2&gt; &lt;p&gt; 

This module centralizes all hyperparameters for non-linear optimization. Separating configuration ensures reproducibility, readability, and simplified parameter tuning. 
&lt;/p&gt; &lt;h3&gt;Example&lt;/h3&gt; 
&lt;pre&gt; def get_diaplex_50_optimization_params():
 dict = {}
 dict["lr"] = 1e-3
 dict["n_iter"] = 40000
 dict["data_range"] = [5e-6,1e4]
 dict["w_threshold"] = 0.0032
 dict["w_left"] = 1
 dict["w_right"] = 2
 dict["criterion"] = torch.nn.MSELoss(reduction="sum")
 return dict &lt;/pre&gt; 

&lt;h3&gt;Key Parameters&lt;/h3&gt; &lt;ul&gt; &lt;li&gt;&lt;code&gt;lr&lt;/code&gt; → learning rate for Adam optimizer&lt;/li&gt; &lt;li&gt;&lt;code&gt;n_iter&lt;/code&gt; → number of optimization iterations&lt;/li&gt; &lt;li&gt;&lt;code&gt;data_range&lt;/code&gt; → frequency range for fitting&lt;/li&gt; &lt;li&gt;&lt;code&gt;w_threshold&lt;/code&gt; → frequency of peak loss factor&lt;/li&gt; &lt;li&gt;&lt;code&gt;w_left&lt;/code&gt;, &lt;code&gt;w_right&lt;/code&gt; → weighting around &lt;code&gt;w_threshold&lt;/code&gt;&lt;/li&gt; &lt;li&gt;&lt;code&gt;criterion&lt;/code&gt; → loss function (L2 norm)&lt;/li&gt; &lt;/ul&gt;

&lt;h2&gt;3. Model Fitting (&lt;code&gt;main.ipynb&lt;/code&gt;)&lt;/h2&gt;

&lt;p&gt;
The Jupyter notebook orchestrates the full workflow for fitting &lt;strong&gt;fractional viscoelastic models&lt;/strong&gt; to experimental mechanical data. 
It performs data loading, preprocessing, model definition, and optimization of material parameters using &lt;a href="https://pytorch.org/docs/stable/generated/torch.optim.Adam.html"&gt;Adam optimizer&lt;/a&gt;.
&lt;/p&gt;

&lt;h3&gt;Workflow Overview&lt;/h3&gt;
&lt;ol&gt;
&lt;li&gt;&lt;strong&gt;Load configuration and utilities&lt;/strong&gt;&lt;br&gt;
Imports helper functions and parameters from &lt;code&gt;configuration_file.py&lt;/code&gt;:
&lt;ul&gt;
 &lt;li&gt;&lt;code&gt;get_diaplex_filenames()&lt;/code&gt; — returns paths to experimental CSV data.&lt;/li&gt;
 &lt;li&gt;&lt;code&gt;get_diaplex_50_optimization_params()&lt;/code&gt; — provides hyperparameters (learning rate, iteration count, weighting, etc.).&lt;/li&gt;
&lt;/ul&gt;
Initializes the fractional viscoelastic model &lt;code&gt;SecondOrderFractionalDerivativeModel&lt;/code&gt; for each dataset.
&lt;/li&gt;

&lt;li&gt;&lt;strong&gt;Load and preprocess experimental data&lt;/strong&gt;&lt;br&gt;
Reads CSV files, extracts and sorts:
&lt;ul&gt;
 &lt;li&gt;Extensional Reduced Frequency&lt;/li&gt;
 &lt;li&gt;Extensional Storage Modulus (&lt;em&gt;E′&lt;/em&gt;)&lt;/li&gt;
 &lt;li&gt;Extensional Loss Modulus (&lt;em&gt;E″&lt;/em&gt;)&lt;/li&gt;
&lt;/ul&gt;
Data are converted to PyTorch tensors for GPU-based optimization.
&lt;/li&gt;

&lt;li&gt;&lt;strong&gt;Define training configuration&lt;/strong&gt;&lt;br&gt;
For each dataset:
&lt;ul&gt;
 &lt;li&gt;Finds the frequency corresponding to the loss peak (&lt;code&gt;tan δ = E″/E′&lt;/code&gt;).&lt;/li&gt;
 &lt;li&gt;Loads optimization settings: &lt;code&gt;criterion&lt;/code&gt;, &lt;code&gt;lr&lt;/code&gt;, &lt;code&gt;n_iter&lt;/code&gt;, &lt;code&gt;w_left&lt;/code&gt;, &lt;code&gt;w_right&lt;/code&gt;, and &lt;code&gt;data_range&lt;/code&gt;.&lt;/li&gt;
 &lt;li&gt;Computes &lt;code&gt;data_scaling&lt;/code&gt; to weight frequency ranges differently around the loss peak.&lt;/li&gt;
&lt;/ul&gt;
&lt;/li&gt;
&lt;/ol&gt;

&lt;h3&gt;Optimization Loop&lt;/h3&gt;
&lt;pre&gt;
for t in tqdm(range(n_iter)):

 pred = models[key](angfreq[key][indices[key]])
 meas = ((lossmod[key])/(stormod[key]))[indices[key]]
 loss = criterion(meas/meas * data_scaling, pred/meas * data_scaling)

 optimizer.zero_grad()
 loss.backward()
 optimizer.step()
&lt;/pre&gt;

&lt;p&gt;
The optimization minimizes the difference between measured and predicted loss factor (&lt;em&gt;tan δ&lt;/em&gt;). 
Weighting (&lt;code&gt;w_left&lt;/code&gt;, &lt;code&gt;w_right&lt;/code&gt;) emphasizes frequencies below or above the loss peak for improved robustness. 
Convergence depends strongly on initial guesses, data range, and learning rate.
&lt;/p&gt;

&lt;h3&gt;Visualization and Output&lt;/h3&gt;
&lt;p&gt;
After convergence, the notebook produces log–log plots comparing model predictions and measurements:
&lt;/p&gt;
&lt;ul&gt;
 &lt;li&gt;&lt;strong&gt;Loss factor&lt;/strong&gt; (&lt;em&gt;tan δ&lt;/em&gt;)&lt;/li&gt;
 &lt;li&gt;&lt;strong&gt;Storage modulus&lt;/strong&gt; (&lt;em&gt;E′&lt;/em&gt;)&lt;/li&gt;
&lt;/ul&gt;
&lt;p&gt;
Figures are saved as &lt;code&gt;*_fig_loss_factor.pdf&lt;/code&gt; and &lt;code&gt;*_fig_storage_modulus.pdf&lt;/code&gt;, ready for publication.
&lt;/p&gt;

&lt;h3&gt;Key Notes&lt;/h3&gt;
&lt;ul&gt;
 &lt;li&gt;&lt;strong&gt;Convergence sensitivity:&lt;/strong&gt; Good initial parameter guesses and suitable step sizes are essential.&lt;/li&gt;
 &lt;li&gt;&lt;strong&gt;Weighting strategy:&lt;/strong&gt; Improves robustness near the loss peak, particularly for noisy data.&lt;/li&gt;
 &lt;li&gt;&lt;strong&gt;Modular design:&lt;/strong&gt; New datasets or model variants can be added easily through configuration files.&lt;/li&gt;
&lt;/ul&gt;


&lt;h2&gt;4. Fitted Parameters&lt;/h2&gt;&lt;p&gt;
In addition to the Python scripts for non-linear optimization, this repository also provides the fitted parameters of each master curve. These parameters, published in the related dataset 
(&lt;a href="https://doi.org/10.18419/DARUS-4318"&gt;Data Sets for a Fractional Moisture-Dependent Viscoelasticity Model for Thermoplastic Polymers&lt;/a&gt;), are available for both the LSF model and the Generalized Maxwell model.
&lt;/p&gt;