Authors: Roman Podieiablonskii , Hyoju Lim , Enes Çınar
Supervisor: Felix Schneider
1. Introduction
In computational engineering, the accuracy of simulation models depends on various assumptions and simplifications which define their underlying mathematical framework. These assumptions are necessary to make complex problems solvable, but they can introduce inconsistencies between simulated results and real-world measurements. As a result, predictions obtained from these models may show different results from experimental or field data, which affects their reliability in practical applications.
One of the main reasons for the error is the inaccurate selection of model parameters. Material properties, boundary conditions, and structural characteristics often need to be fine-tuned to ensure the model represents physical behavior accurately. However, due to uncertainties in material properties, manufacturing tolerances, and external environmental factors, it's challenging to determine the most appropriate parameter values. If the chosen parameters do not reflect the actual system characteristics, the model's predictive capability is compromised, leading to unreliable results.
To solve these problems, Sensitivity-Based Finite Element Model Updating (FEMU) approaches have been developed as a framework to enhance model accuracy. These approaches use measurement data such as displacements, eigenmodes, or strain responses to update the model parameters, reducing inconsistency between numerical simulations and real-world observations. By integrating optimization techniques with experimental validation, FEMU approaches provide a robust methodology for improving computational models, which makes them more reliable for modelling of complex phenmena and their analysis.
2. FEMU Approach
2.1. General Idea of Sensitivity-based FEMU Approach
Model parameters can be updated using measurements, such as displacements or eigenmodes, to obtain a more accurate solution.
The update procedure follows the least squares method, a parameter estimation technique that minimizes the sum of the squares of the residuals. The residual \bm{\varepsilon}_z is the difference between measured data (denoted by the vector \bm{z_m} ) and the values predicted by a model (denoted by the vector \bm{z(\theta)} ), which depends on vector \bm{\theta} , containing model parameters to be updated.
\begin{equation} \bm{\varepsilon_z} = \bm{z_m} - \bm{z}(\bm{\theta}) \approx \bm{r}_i - \bm{G}_i (\bm{\theta} - \bm{\theta}_i) \end{equation} |
\begin{equation} \bm{r_i} = \bm{z_m} - \bm{z_i} \end{equation} |
Sensitivity-based FEMU approach builds upon the linearization of the generally non-linear residual function, so it’s developed from a Taylor series expansion of \bm{\varepsilon}_z truncated after the linear term, where \( \bm{G}_i\) denotes the sensitivity matrix which entries are the derivatives of model functions with respect to the model parameters. For each iteration of the update sensitivity matrix \( \bm{G}_i\) is computed at the current value of the complete vector of parameters \( \bm{\theta}_i \) . Then, using this linearized formulation of the residual, the ordinary least squares problem can be formulated, this is denoted by the function J(\bm{\theta}) , and the objective is to minimize this function.
\[ J(\bm{\theta}) = \bm{\varepsilon}_z^T \bm{W}_{\varepsilon} \bm{\varepsilon}_z \] |
Now having introduced the basic idea of sensitivity-based model updating approaches, it is time to dive into the actual algorithm how model parameters are being updated using ordinary least squares method.
2.2. Derivation of the Updating Step
To improve model accuracy, we need to update the parameters iteratively. This requires minimizing the objective function, which measures the difference between predicted and measured values. The objective function is defined as:
\begin{equation} J(\theta) = \boldsymbol{\varepsilon}^T W_{\varepsilon} \boldsymbol{\varepsilon} \end{equation} |
where \bm{W_{\varepsilon}} denotes the weighting matrix. \bm{W_{\varepsilon}} is difficult to estimate, although at the very least this weighting should include scaling to equalize the effect of amplitude and a reasonable choice is
\begin{equation} W_{\varepsilon} = \left[ \text{diag} (z_m) \right]^{-2}. \end{equation} |
In the following, expression for the updating step \bm{\Delta \theta} is derived. For the simplification a single model parameter and one measured value are considered such that vector quantities become scalars $\bm\theta \to \theta$ , \bm{\varepsilon_z} \to \varepsilon} . Taking the derivative with respect to \theta and setting it to zero.
\begin{equation} \frac{\partial J(\theta)}{\partial \theta} = \frac{\partial}{\partial \theta} (\boldsymbol{\varepsilon}^T W_{\varepsilon} \boldsymbol{\varepsilon}) = 0 \end{equation} |
Inserting the Taylor approximation truncated after linear term:
\begin{equation} \frac{\partial J(\theta)}{\partial \theta} = \frac{\partial}{\partial \theta} (W_{\varepsilon} \boldsymbol{\varepsilon}^2) \approx \frac{\partial}{\partial \theta} \left[ W_{\varepsilon} (r_i - G_i (\theta - \theta_i))^2 \right] \end{equation} |
Taking derivative:
\begin{equation} W_{\varepsilon} (r_i - G_i \Delta \theta_i) \cdot 2 \cdot (-G_i) = 0 \end{equation} |
Rearranging:
\begin{equation} G_i W_{\varepsilon} r_i - G_i W_{\varepsilon} G_i \Delta \theta_i = 0 \end{equation} |
Solving for \Delta \theta :
\begin{equation} \Delta \theta_i = (G_i W_{\varepsilon} G_i)^{-1} G_i W_{\varepsilon} r_i \end{equation} |
Here it is of high importance to mention that a necessary (but not necessarily sufficient) condition for obtaining a unique solution is number of measured values should always be larger than number of model parameters to be updated. Otherwise, problem of ill-conditioned sensitivity matrix arrises.
2.3. Generic Algorithm
For optimization, as a first step, we define the model parameter theta, θ. This might be stiffness for a mass-spring system, or geometry for a beam, such as cross-sectional width or thickness. Then, we calculate the residual, which depends on the problem, it might based on forces, eigenvalues, or eigenvectors.
The optimization algorithm includes some steps and these are explained below.
1) Firstly we compute the sensitivity matrix. This matrix shows how small changes in parameters affect outputs.
G_{i} = \left[ \frac{\partial z_{j}}{\partial \theta_{k}} \right]_{\theta = \theta_{i}} |
2) Next, we compute the residual for the current step.
r_i = z_m - z_i |
3) \bm{\theta}_i , which stands for minimizing the objective function, J is computed. To stabilize the solution, and for ill-conditioned problems, we use a regularization parameter, μ. Then, we check convergence by comparing with a predefined tolerance value.
\Delta \theta_i = [G_i^T W_\varepsilon G_i + \mu^2I]^{-1} G_i^T W_\varepsilon r_i |
4) If convergence is not achieved, parameters are updated for the next iteration step, and this whole process is iterated until we get a convergence.
\theta_{i+1} = \theta_i + \Delta \theta_i |
3. Case Analysis
Sensitivity-based model updating approach is applied to three cases to demonstrate the effectiveness of the parameter updating approach. The analysis begins with an initial guess for the model parameters, followed by computing the residuals based on measured and predicted values. Parameters are updated using the derived update equation to minimize the objective function. This case study provides insight into the convergence behaviour, accuracy improvement, and robustness of the sensitivity-based model updating approach.
3.1. Static Case
In the static case, the parameter update process is driven by force-based residuals. Stiffness values k_1, k_2, k_3 are used as model parameters to be updated.
At the initial step (before starting the generic algorithm), the force is computed using the initially guessed stiffness and measured displacement values, defining \bm{f_1} as the initial force vector. Then, the initial residual \bm{r_1} = \bm{f_m} - \bm{f_1} can be constructed. Finally, using \bm{r_1} for the initial iteration the generic algorithm can be started.
With these force-based residuals, a generic algorithm is employed to iteratively refine the model parameters, ensuring the stiffness values are updated to minimize the discrepancy between predicted and measured forces. This force-driven approach is suited for static problems, where equilibrium conditions govern the system's behaviour. The key aspect of this approach is that the residual is calculated based on force differences, obtained by subtracting the computed force vector \bm{f_i} corresponding to the i-th iteration from the measured force vector \bm{f_m} .
3.1.1. Workflow
1) Define initial guess for stiffness values: \begin{bmatrix} k_1 \\ k_2 \\ k_3 \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}
2) Define measured displacements: \begin{equation} u_m = \begin{bmatrix} u_1 \\ u_2 \end{bmatrix} \end{equation}
3) Compute initial force vector: \bm{f_1} = \begin{pmatrix} f_1 \\ f_2 \end{pmatrix} = \begin{bmatrix} k_1 + k_2 & -k_2 \\ -k_2 & k_2 + k_3 \end{bmatrix} \begin{bmatrix} u_{m1} \\ u_{m2} \end{bmatrix}
5) Compute initial residual based on force vectors: \bm{r_1} = \bm{f_m} - \bm{f_1}
6) Use generic algorithm until converged to update parameters
3.1.2. Conditions
We set both underdetermined case and overdetermined case. Underdetermined case have less number of equations than unknowns whilst overdetermined case have more equations than unknowns.
|
Underdetermined Case |
Overdetermined Case |
Initial Guess |
[1, 1, 1] |
|
True Value |
[k1 = 1.2; k2 = 0.8; k3 = 1.0;] |
|
Number of Equations |
2 |
4 |
Number of Unknowns |
3 |
3 |
Measured static displacements(u_m) |
u_1 = [0.60811; 0.27027] |
u_1 = [0.60811; 0.27027] u_2 = [0.27027; 0.67568] |
Measured forces (f_m) |
f_1 = [1; 0] |
f_1 = [1; 0] f_2 = [0; 1] |
3.1.3. Results
1) Underdetermined Case:
|
|
---|---|
Convergence of Estimated Parameters vs True Values | Convergence of Error |
The convergence of the error does not reach zero since there are fewer equations than unknowns in underdetermined case, which makes it impossible to obtain an exact solution. Due to insufficient constraints, multiple parameter sets can satisfy the equations, and algorithm reaches a state where further iterations do not improve the accuracy, as the missing information prevents the system from fully resolving all unknowns. Therefore, the updating process may reduce the error to a certain extent but cannot eliminate it. To solve this issue, additional constraints or regularization techniques are often required to obtain a more accurate and stable solution.
2) Overdetermined Case:
|
|
---|---|
Convergence of Estimated Parameters vs True Values | Convergence of Error |
The number of equations is higher than the number of unknowns in an overdetermined case, and relation between forces (residual) and stiffness (updated parameters) is linear in the static case. Therefore the convergence of the error is achieved in only one iteration, obtaining the solution identical to the measured data. The system in this case is well-defined, therefore the algorithm finds the best-fit parameters that minimize the objective function, leading to full convergence.
3.2. Dynamic Case 1
For the dynamic scenario, a two-degree-of-freedom system is considered, including two masses(\bm{m_1} , \bm{m_2} = 1) connected by two springs. The main objective is to update the stiffness values, \bm{k_1} and \bm{k_2} to their true values of \bm{k_{1,true} = 1.2 N/m} and \bm{k_{2,true}= 0.8 N/m} , respectively. Two methodologies are employed for this purpose; the first one is eigenvalue residual-based updating and the second one is eigenvector residual-based updating. These methods update the stiffness parameters during iteration to minimize the errors between predicted and true modal properties.
3.2.1. Workflow
1) Define initial guess for stiffness values: \bm{\theta = \begin{bmatrix} k_1 \\ k_2 \end{bmatrix}}
2) Define true/measured stiffness values: \begin{equation} \theta_m = \begin{bmatrix} k_{1,true} \\ k_{2,true} \end{bmatrix} \end{equation}
3) Construct stiffness and mass matrices: \mathbf{K} = \begin{bmatrix} k_1 + k_2 & -k_2 \\ -k_2 & k_2 \end{bmatrix}, \mathbf{M} =\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}
4) Solve the eigenvalue problem \begin{equation} det(\bm K - \lambda_i \bm M) = 0\end{equation} to find eigenvectors and eigenvalues with both guessed and true stiffness values separately.
\mathbf{\phi} = \begin{bmatrix} \phi_{1,1} & \phi_{1,2} \\ \phi_{2,1} & \phi_{2,2} \end{bmatrix}, \mathbf{\lambda} =\begin{bmatrix} \lambda_1 \\ \lambda_2 \end{bmatrix}
5) Compute initial residual \bm{r_1} based on eigenvalues or eigenvectors:
i) For the eigenvalue residual method:
\bm{r_1} = \bm{\lambda_m} - \bm{\lambda_{1}}
ii) For the eigenvector residual method:
\bm{r_1} = \bm{\phi_m} - \bm{\phi_{1}}
where \bm{\lambda_{1}} and \bm{\phi_1} denote eigenvalues and eigenvectors, respectively, computed using initial stiffness values
6) Use generic algorithm until converged to update model parameters. For each method, the sensitivity matrix can be computed using the following analytical expression:
i) For the eigenvalue residual method:
G_{jk} = \frac{\partial \lambda_j} {\partial \theta_k} = \bm \phi_j^T \left[ -\lambda_j \frac{\partial \mathbf{M}}{\partial \theta_k} + \frac{\partial \mathbf{K}}{\partial \theta_k} \right] \bm \phi_j
ii) For the eigenvector residual method:
a_{jkh} = \frac{\varphi_h^T \left( -\lambda_j \frac{\partial \mathbf{M}}{\partial \theta_k} + \frac{\partial \mathbf{K}}{\partial \theta_k} \right) \varphi_j}{\left( \lambda_j - \lambda_h \right)} , h \neq j
\frac{\partial \Phi_j}{\partial \theta_k} = \sum_{h=1}^H a_{jkh} \Phi_h;
3.2.2. Results
As a first step, we start with an initial guess \bm{\theta = \begin{bmatrix} k_1 = 1 \\ k_2 = 1 \end{bmatrix}} , and we choose \bm{\mu = 0.001} , and as we expected both methods successfully update the parameters as seen in the below figures. The eigenvalue residual method achieves convergence within five iterations. On the other hand, the eigenvector residual method shows a faster convergence rate (4 iterations) but results in an oscillatory behaviour.
Convergence Issue
When modifying the initial guesses as \bm{\theta = \begin{bmatrix} k_1 = 0.5 \\ k_2 = 1.3 \end{bmatrix}} , the eigenvector residual method struggles to converge, as evident in the plot below. This issue can be weakened by reducing the regularization parameter to a small value such as \bm{\mu = 1x10^{-9}} through trial-and-error tuning. In contrast, the eigenvalue residual method works perfectly regardless of the starting point or regularization amplitude highlighting its robustness. These results suggest that while the eigenvector residual requires careful parameter tuning, whereas the eigenvalue residual method provides a more robust solution framework.
Model updating based on eigenvector residual method with different regularization parameters.
Model updating with eigenvalue residual method with different regularization parameters.
3.3. Dynamic Case 2
In the second dynamic case, FE-approximation of a free vibrating beam is considered. The problem setup is the following: 2d beam of length L = 0.3 m , cross-sectional thickness \bm{t} , width b , discretized with n = 30 elements, with the fixed-end at x = 0 and free end at x = L . There are no external forces applied since we are interested in the free vibration. The investigated case, its dimensions and corresponding coordinate system are summarized in the figure below.
Every finite element is of length l_e = 0.01m and has two nodes, with 3 DOFs per node: horizontal, vertical and rotational around y-axis. Hence, this results in 6 DOFs per element with 6-by-6 elemental mass and stiffness matrices \bm{M_e} and \bm{K_e} , respectively.
\[ K_e = \begin{bmatrix} \frac{AE}{L} & 0 & 0 & -\frac{AE}{L} & 0 & 0 \\ 0 & \frac{12EI}{L^3} & \frac{6EI}{L^2} & 0 & -\frac{12EI}{L^3} & \frac{6EI}{L^2} \\ 0 & \frac{6EI}{L^2} & \frac{4EI}{L} & 0 & -\frac{6EI}{L^2} & \frac{2EI}{L} \\ -\frac{AE}{L} & 0 & 0 & \frac{AE}{L} & 0 & 0 \\ 0 & -\frac{12EI}{L^3} & -\frac{6EI}{L^2} & 0 & \frac{12EI}{L^3} & -\frac{6EI}{L^2} \\ 0 & \frac{6EI}{L^2} & \frac{2EI}{L} & 0 & -\frac{6EI}{L^2} & \frac{4EI}{L} \end{bmatrix} \] |
Where \[ I = \frac{bt^3}{12} \] and \[ A = bt \]
After assembly process, global \bm{K} and \bm{M} are 93-by-93 matrices, which are reduced to 90-by-90 matrices by implementing the fixed-end boundary condition at x = 0 , so at the first FE-node.
For the second dynamic case the model updating procedure was based on eigenvalue residual.
3.3.1. Workflow
Two experiments were conducted with different number of model paramters to be updated: only beam thickness t and two updated parameters b and t . Now it is of high interest to present the workflow considering the exemplary case of upadating only one model parameter: thickness \bm{t} .
FEMU approaches are used to bring model parameters in accordance with measurements, therefore, in the absence of actual measurement data one needs to obtain "synthetic" measurement data. This was done by defining „true“ beam thickness t_m . Then using t_m „true“ eigenvalues were computed and measurement data contained in vector \bm{z_m} was simulated by adding some noise to „true“ eigenvalues (the noise follows normal distribution with zero mean and 10^6 standard deviation).
1) Define "true" thickness \begin{equation} \theta_m = t_m\end{equation}
2) Solve the eigenvalue problem \begin{equation} det(\bm K - \lambda_i \bm M) = 0\end{equation} where \bm{K} and \bm{M} are global stiffness and mass matrices, respectively, which are computed using \begin{equation} \theta_m = t_m\end{equation} , to find "true" eigenvalues.
3) Construct measurement data vector \begin{equation} \bm{z_m} = \bm \lambda + N(\bm 0, 10^6 \times \bm I)\end{equation} by adding measurement noise (10^6 standard deviation corresponds to the order of magnitude of system's eigenvalues)
Next step is to formulate the initial residual \bm{r_1} by subtracting simulated eigenvalues contained in vector \bm{z_1} from measured eigenvalues contained in vector \bm{z_m} . For that, initial guess of thickness t_1 was taken to compute eigenvalues in vector \bm{z_1} .
4) Define initial thickness t_1
5) Calculate eigenvalues contained in vector \begin{equation} \bm{z_1} = \bm{\lambda}\end{equation} by solving the eigenvalue problem as in step (2) with \bm{K} and \bm{M} computed using initial thickness t_1
6) Build up initial residual \begin{equation} \bm{r_1} = \bm{z_m} - \bm{z_1}\end{equation}
7) Use generic algorithm until converged to update model parameters. Sensitivity matrix entries computed by G_{jk} = \frac{\partial \lambda_j} {\partial \theta_k} = \bm \phi_j^T \left[ -\lambda_j \frac{\partial \mathbf{M}}{\partial \theta_k} + \frac{\partial \mathbf{K}}{\partial \theta_k} \right] \bm \phi_j , where \bm{K} and \bm{M} are global stiffness and mass matrices computed every iteration step i with corresponding updated thichness t_i , and \bm{\phi_j} is the eigenvector corresponding to the j-th eigenvalue \lambda_j
3.3.2. Results
1) One Updated Parameter: Thickness \bm{t}
Here true thickness \begin{equation} t_m\end{equation} was set to 0.02 m and initial guessed thickness t_1 was chosen to be equal to 0.025 m .
Considering figures above it can be concluded that in case of only one updated parameter the algorithm converges relatively fast within only 5 iterations. Convergence criterion was set to \begin{equation}\Delta \theta < 10^{-6} \end{equation} . Here it is of high interest to notice that the eigenvalue residual method works perfectly for this case regardless of the starting point or regularization amplitude highlighting its robustness. This could be explained by the fact that in this case the objective function J(\bm{\theta}) depends only on one parameter t and has one global minimum.
2) Two Updated Parameters: Thickness \bm{t} and Width \bm{b}
Here true thickness \begin{equation} t_m\end{equation} and true width \begin{equation} b_m\end{equation} was set to 0.02 m and 0.03 m , respectively, and initial guessed thickness t_1 and width b_1 both were chosen to be equal to 0.025 m .
Considering figures above it can be concluded that in case of two updated parameters the algorithm based on the eigenvalue residual fails to converge in both parameters to "true" values. Although, thickness t is being updated to true value \begin{equation} t_m\end{equation} , width b is being updated in the "right" direction (towards \begin{equation} b_m\end{equation} ), but stops half-way. Moreover, it is important to notice that in this case algorithm requires significantly more iterations for convergence, making overall sensitivity-based FEMU computationally more expensive. Such behaviour could be explained by consideration of the absolute value of the objective function J(\bm{\theta}) in case of two updated parameters depicted in the figures below:
Considering figure on the right one observes that the objective function does not have one global minimum but instead there are two ridges (or valleys). Therefore, it is almost impossible to obtain optimal minimum in the crossing of two valleys, because once trapped in the valley (meaning that one parameter, for instance t , has already converged to true value \begin{equation} t_m\end{equation} or to the close vicinity of true value), gradient with respect to the second updated parameter becomes almost zero, so it cannot be properly updated. Indeed, updated thicknesses t_i converge faster than updated widths \begin{equation} b_i\end{equation} (can be seen in the plots depicting evolution of beam thickness and beam width over iteration), which results in the fact that algorithm is trapped in the valley corresponding to \begin{equation} t_m\end{equation} far from the optimal minimum, and each subsequent updating step for \begin{equation} b_i\end{equation} becomes smaller and smaller asymptotically converging to zero. Convergence criterion was set to \begin{equation}\Delta \theta < 10^{-6} \end{equation} , hence, it is of high interest to analyze algorithm's behaviour for substantially smaller convergence criterion \begin{equation}\Delta \theta < 10^{-12} \end{equation} , which is depicted in the figures below:
Indeed, considering these plots it can be concluded that algorithm is trapped in the valley of the objective function J(\bm{\theta}) corresponding to \begin{equation} t_m\end{equation} where gradients with respect to b are almost zero resulting in asymtotic convergence of updated beam widths \begin{equation} b_i\end{equation} towards "wrong" value of 0.0276 m .
4. Conclusion
4.1. Project Summary
For the static case, the parameter updating approach successfully refines the model by minimizing residuals based on forces. Since the residuals are computed by the difference between the measured and computed forces, the optimization process adjusts the stiffness parameters during iteration to achieve better agreement between the model and experimental data. This implies that the sensitivity-based approach demonstrates the effectiveness of force-based residual minimization in model updating. It improves simulation accuracy and makes computational predictions align more accurately with real-world measurements. However, in the underdetermined case, where there are more unknowns than equations, the solution may converge to an incorrect value due to insufficient constraints.
For the dynamic case 1, the main goal was to update two stiffness values to their true/measured values. As we see in the results, the eigenvalue residual method proves to be robust and reliable across different initial conditions and regularization parameter, while the eigenvector residual stands more sensitive to the initial point and parameter tuning. These findings highlight the importance of selecting an appropriate updating methodology and in case fine tuning parameters to achieve a good result.
For the dynamic case 2, the main goal was to perform an actual sensitivity-based Finite Element Model Updating based on the FE-approximation of a free vibrating beam with one fixed and one free end. Results for the case of only one updated parameter beam thickness demonstrated that the eigenvalue residual method performed exceptionally well. However, in case of two (independent) model parameters algorithm failed to deliver correct results: beam thickness t did converge to the predefined true value, but beam width b failed to do so, although, overall direction of the update was captured correctly. Such behavior was explained by the nature of the objective function J(\bm{\theta}) , which was non-convex, and did not have one global minimum or even multiple local minima. Therefore, it can be concluded that choice of independent model parameters does not guarantee convergence of the algorithm, and, if possible, topology of the objective function J(\bm{\theta}) should be carefully analyzed beforehand.
4.2. Future Work
For future work, damping effects for dynamics cases can be included and analyzed. Also further investigations might be conducted to improve the stability of eigenvector residual methods, mitigating sensitivity to parameter tuning. Moreover, for the second dynamic case it is of high interest to try out other pairs of independent model parametrs to see, if problem of non-convex objective function remains present.
5. References
- 2018, John E. Mottershead et al. “The sensitivity method in finite element model updating: A tutorial”
- Mottershead, J.E., Link, M., Friswell, M.I., Schedlinski, C. (2020). Model Updating. In: Allemang, R., Avitabile, P. (eds) Handbook of Experimental Structural Dynamics. Springer, New York, NY. https://doi.org/10.1007/978-1-4939-6503-8_18-1