Transient Adjoint Sensitivity Analysis

Authors: Jakob Scheffels, Viet Ha Dinh

Supervisors: Jonas Hillenbrand, Sebastian Schopper

Introduction

Many gradient-based optimization problems require design sensitivity analyses to quantify the changes in performance of the design with respect to design parameters. In structural dynamics, it is common that optimization problems depend on a large amount of design parameters [3], and so computing the sensitivities – in the form of gradients – with the finite difference approach is extremely inefficient and time-consuming. Another approach for computing the gradients is the adjoint method, which is very efficient for a large amount of parameters. This project aims to implement the adjoint method in sensitivity analysis for transient problems and validate the implementation with practical examples such as the performance of a tuned mass damper with respect to its design parameters. The implementation of the adjoint sensitivity analysis is done using with MATLAB with the package bm_mfem as the basis.

Sensitivity Analysis

Optimization problem

First an optimization problem is set up.

Consider a single degree of freedom system consisting of a mass, a spring and a damper. A harmonic load $\begin{array}{l}F(t) = F_0\cos(\omega t)\end{array}$ is applied to the mass element.

We have:

Design variables : s

e.g.: $\begin{array}{l}[𝑴&𝑪&𝑲]^𝑇\end{array}$

State variables: $\begin{array}{l}𝒖(𝒔, 𝑡)\end{array}$ with t : time

Primal problem

Equation of motion for transient problem :

$\begin{array}{l}𝑹(𝒔, 𝒖, 𝒖 ̇, 𝒖 ̈,𝑡)=𝑴(𝒔) 𝒖 ̈+𝑪(𝒔) 𝒖 ̇+𝑲(𝒔)𝒖−𝑭(𝒕) \equiv 0 \\ 𝒖(𝑡=0)= 𝒖_0 \\ 𝒖 ̇(𝑡=0)= 𝒖 ̇_0\end{array}$

The above equation represents the primal problem in which the displacement $\begin{array}{l}u\end{array}$ , the velocity $\begin{array}{l}\dot{u}\end{array}$ , the acceleration $\begin{array}{l}\ddot{u}\end{array}$ are solved iteratively by the Newmark method.

Response function: $\begin{array}{l}𝑱(𝒔, 𝒖, 𝒖 ̇ )= \int_0^T 𝒋(𝒔, 𝒖, 𝒖 ̇ )dt\end{array}$

In structural dynamics, a common response function is the energy dissipation function : $\begin{array}{l}𝑾_𝒄= −\int_𝟎^𝑻 𝒖 ̇^𝑇 𝑪 𝒖 ̇ 𝑑𝑡\end{array}$

The sensitivity analysis is applied to assess the response function of the system with respect to its design parameters. This can be done by computing the derivatives of response function $\begin{array}{l}J\end{array}$ to design parameter $\begin{array}{l}s_i\end{array}$ for $\begin{array}{l}i\end{array}$ as the index of the respected element in the system.

Local sensitivity:

$\begin{array}{l}\displaystyle \frac{dJ}{ds_i} = \int_o^T \frac{dj}{ds_i}dt = \int_o^T (\frac{dj}{ds_i} + \frac{dj}{du} \cdot \frac{du}{ds_i} + \frac{dj}{d \dot{u}} \cdot \frac{d \dot{u}}{ds_i} + \frac{dj}{d \ddot{u}} \cdot \frac{d \ddot{u}}{ds_i} )dt\end{array}$

Adjoint problem

We introduce the Lagrange equation by adding the product of the lambda multiplier and the equation of motion to the response function. For a system in equilibrium, R = 0, therefore, adding it will not deviate the sensitivity results.

$\begin{array}{l}\displaystyle L(u,\dot{u}, \ddot{u}, s) = J + \int_0^T \lambda^T Rdt = \int_0^T({j + \lambda^TR})dt\end{array}$

Differentiate the Lagrange equation, we have the sensitivity:

$\begin{array}{l}\displaystyle \frac{d{J}}{d{s_i}} = \frac{d{L}}{d{s_i}} = \int_0^T \left ( \frac{\partial{j}}{\partial{s_i}} + {\lambda}^T \cdot \frac{\partial{R}}{\partial{s_i}} \right)dt\end{array}$

Throughout further derivation, specifically inserting the derivative of R to the sensitivity and integrating by parts (see Appendix), we achieve a second-order differential equation where lambda is the unknown as below.

$\begin{array}{l}M^T\frac{d\dot{\lambda}}{dt} - C^T\dot{\lambda} + K^T\lambda = - \left (\frac{\partial j}{\partial u} \right)^T + \dot{u}^T + \left ( \frac{\partial^2{j}}{\partial u \partial{\dot{u}}} \right) ^T + \ddot{u}^T \left (\frac{\partial^2 j}{\partial{\dot{u^2}}} \right )^T \\ \lambda(T) = 0 \\ M^T\dot{\lambda}(T) = -\left (\frac{\partial j}{\partial \dot{u}} \right)_{t=T}\end{array}$

We can see that the derivatives $\begin{array}{l}\frac{du}{ds_i}, \frac{d \dot{u}}{ds_i}\end{array}$ and $\begin{array}{l}\frac{d \ddot{u}}{ds_i}\end{array}$ are omitted, thus there is no dependency of the design variables in this adjoint problem.

Since the conditions are given at time t = T instead of t = 0, we need to solve the equation backwards in time. Therefore, we introduce the substitution:

$\begin{array}{l}\tau = T - t \\ \lambda(t) = \Lambda(T-t) = \Lambda(\tau), \dot{\lambda}(t) = - \dot{\Lambda}(\tau), \ddot\lambda(t) = \ddot\Lambda(\tau)\end{array}$

Thus, the adjoint problem reads:

$\begin{array}{l}M^T\frac{d\dot{\Lambda}}{dt} + C^T\dot{\Lambda} + K^T\Lambda = - \left (\frac{\partial j}{\partial u} \right)^T + \dot{u}^T + \left ( \frac{\partial^2{j}}{\partial u \partial{\dot{u}}} \right) ^T + \ddot{u}^T \left (\frac{\partial^2 j}{\partial{\dot{u^2}}} \right )^T \\ \Lambda(0) = 0 \\ M^T\dot{\Lambda}(0) = -\left (\frac{\partial j}{\partial \dot{u}} \right)_{t=T}\end{array}$

We solve the adjoint problem for lambda, and then insert lambda back to $\begin{array}{l}\frac{dJ}{ds}\end{array}$ for our sensitivities.

$\begin{array}{l}\displaystyle \frac{dJ}{ds} = \int_0^T {\left ( \frac{\partial{j}}{\partial{s_i}} + \Lambda^T(T-t) \frac{\partial{R}}{\partial{s_i}} \right )dt}\end{array}$

The adjoint problem acquires the use of the transpose matrix of M, C, K and shifting the derivatives to the Lagrange parameters $\begin{array}{l}\lambda\end{array}$ , which is the definition of an 'adjoint operator'. Thus the method is named "adjoint method".

Through this method, we only need to solve for lambda once and then inserting lambda in the integral to compute the sensitivities. Therefore, the adjoint method is very efficient, especially for problems with a large amount of parameters, as we do not need to compute the derivative of u with respect to different parameters.

Examples

Single Degree-of-Freedom System

To show the functionality of the Adjoint Sensitivity Method, we start with a simple mass-spring-damper-system with a single degree of freedom. The system has known design parameters s={m,k,c} and is subjected to a harmonic load with a single excitation frequency $\begin{array}{l}\omega_{ex}\end{array}$ . The system is exposed to an initial displacement $\begin{array}{l}u(0) = u_0\end{array}$ and an initial velocity $\begin{array}{l}v(0) = v_0\end{array}$ . We choose as response function $\begin{array}{l}J = \int_0^T u^2 (t) + v^2 (t) dt\end{array}$ . For this example we chose the following values for the parameters:

$\begin{array}{l}\displaystyle m = 1.0kg, k = 10\frac{N}{m}, c = 2.5 \frac{Ns}{m}, |f| = 2N, \omega_{ex} = 1.0\frac{rad}{s}, T = 10s, dt = 0.01s, u(0) = u_0 = 2m, v(0) = v_0 = -0.4\frac{m}{s}.\end{array}$

Solving our primal problem with a Newmark solver, we get the following displacement of the Dof:

.

Note that our excitation frequency is smaller than the natural eigenfrequency of the system $\begin{array}{l}\omega_n = \sqrt\frac{k}{m} = \sqrt10.\end{array}$ Analyzing the response of the SDOF system in the frequency domain, we get:

From the frequency response, we can deduct an expectation of the sensitivity of the spring-stiffness parameter k. Note that the relation between $\begin{array}{l}\omega_n\end{array}$ and k is proportional to $\begin{array}{l}\sqrt{k}\end{array}$ . Therefore, we can shift the resonance frequency of the system according to k. Since our excitation frequency is smaller than the natural eigenfrequency, we see that a shift of the resonance to the right would lead to a smaller response (equivalent to $\begin{array}{l}\omega_{ex}\end{array}$ becoming smaller), This means that with an increase of k we expect a smaller response. With our choice of J, we expect the sensitivity of J with respect to k to be negative ( $\begin{array}{l}\frac{dJ}{dk} < 0\end{array}$ ).

Solving the adjoint problem in MATLAB for the given parameters, we get the following sensitivities:

We notice that the sensitivity of the stiffness parameter k is positive, which is contradicting our assumption we made based on the frequency analysis. Looking at the displacement plot of the SDOF system on the other hand, we are able to explain this phenomenon: the response of the system is dominated by the initial conditions. Increasing k with fixed u(0) would lead to more strain energy in the spring in the first place (as $\begin{array}{l}E_{init} = 0.5*k*x^2(0)\end{array}$ ). If we run the simulation again with silent initial conditions, we get the following sensitivities:

Here we can see the expected value for the sensitivity $\begin{array}{l}\frac{dJ}{dk} = -0.969 < 0\end{array}$ .

Tower

The implemented adjoint sensitivity method is not limited to simple structures consisting of spring-damper- and concentrated-mass-elements. It can deal with any kind of element and with any amount of degree of freedoms. To show that it's working with more complex models, we build a two-dimensional tower of bar-elements with 22 degrees of freedom. The tower is fixed at the bottom and subjected to a horizontal force acting on the nodes at the top to invoke its first eigenmode.

Not only is the number of Dofs not relevant for the sensitivity analysis, also the number of design parameters $\begin{array}{l}s_i\end{array}$ is irrelevant. In this example, each of the 21 bar elements has three parameters: a Young's modulus (E), its cross-section (A) and the density of the material ( $\begin{array}{l}\rho\end{array}$ ). In this example, we set them all to the same values:

$\begin{array}{l}\displaystyle E = 180 GPa, A = 0.01 m^2, \rho = 8000 \frac{kg}{m^3}\end{array}$

To prevent large shear deformations, we scaled the Young's modulus of the diagonal bar elements by a factor of 100.

After solving the primal problem, we choose the same response function $\begin{array}{l}J = \int_0^T x^2(t)+v^2(t) dt\end{array}$ to solve the adjoint problem and to calculate the sensitivities of our design parameters $\begin{array}{l}s_i\end{array}$ . Even though the corresponding design parameters of the elements have the same value (except for the scaled Young's modulus of the diagonal elements), we expect them to have different sensitivities. After solving the adjoint problem for $\begin{array}{l}\Lambda\end{array}$ once, we are able to get the sensitivity of all parameters by solving the respective integral.

Remember that we have to interpret the sensitivities with respect to our choice of the response function J. Since we chose the same function as in the Single Degree of Freedom example, a negative sensitivity of a design parameter with respect to J means a decrease of the oscillation and velocity value. In this example, we notice that our sensitivities scaled by the parameter value show a trend: while for elements at the bottom of the structure (element 2-6) most sensitivities of the cross section and density are negative, while for the elements at the top (element 17-21) are positive. This can be explained by trying to strengthen the structure at the bottom (increase A and $\begin{array}{l}\rho\end{array}$ ) and reduce the weight at the top (decrease A and $\begin{array}{l}\rho\end{array}$ ) to minimize the response of the structure due to the applied load.

To validate the sensitivities, we adjusted the parameters of the elements according to their sensitivity. Every parameter with a positive sensitivity was scaled by 0.9 and negative sensitivity by 1.1. We ran the primal problem again and compared the value of J for the initial run and the one with the adjusted design parameters and noticed a reduction of J by 14%. We expected a reduction of J as we chose to adjust the parameter in a way to decrease the response function J.

Tuned-Mass-Damper

In the final example we tested the sensitivity analysis on a tuned mass damper attached to a single degree of freedom system with parameters $\begin{array}{l}m = 10 kg, c = 2.0 \frac{Ns}{m}\end{array}$ and $\begin{array}{l}k = 10 \frac{N}{m}\end{array}$ . We chose the ratio of masses $\begin{array}{l}\mu = \frac{m_d}{m} = 0.1\end{array}$ . In literature, there are different criteria to pick 'optimal' design parameters for the tuned mass damper to reduce the movement of the single degree of freedom oscillator [1]. Here we chose two criteria to compare: parameters following Den Hartog criteria and H2 optimization criteria.

$\begin{array}{l}Den Hartog: k_d = k * \frac{\mu}{(1+\mu)^2} = 0.8264 \frac{N}{m}, \\ c_d = \sqrt{k_d * m_d}*\sqrt\frac{3\mu}{(2*(1+\mu)} = 0.3357 \frac{Ns}{m}\\ H2: k_d = \frac{1}{(1+\mu)^2} * \frac{2+\mu}{\mu} * m_d = 17.3554 \frac{N}{m}, \\ c_d = \sqrt{k_d*m_d}*2*\sqrt\frac{\mu*(4+3\mu)}{8*(1+\mu)(2+\mu)} = 1.271 \frac{Ns}{m}\end{array}$

Frequency Analysis

We first perform a frequency analysis of the SDOF system and the two systems with the added tuned mass damper.

As expected, the SDOF system ('Normal') shows large resonance around the natural eigenfrequency of the system $\begin{array}{l}\omega_n = \sqrt\frac{k}{m} = 1\end{array}$ . The two systems with tuned mass damper should reduce the movement of the SDOF oscillator and therefore have a lower response to different excitation frequencies. For the Den Hartog parameters (red curve), we notice a significant drop around the natural eigenfrequency of the SDOF with two peaks shifted to lower and larger frequencies compared to the resonance case of the SDOF system. This behavior can be interpreted as a strong reduction of the movement for excitation frequencies close to the resonance case, but either no different or worse behavior in other frequency ranges. For the H2 optimized parameter values, we see a significant drop at a comparably larger frequency ( $\begin{array}{l}\omega = 4\end{array}$ ) and a slight decrease of the response over a broader frequency range.

With this knowledge, we set up two different cases for the excitation frequency $\begin{array}{l}\omega_{ex}\end{array}$ : first an excitation with a single frequency $\begin{array}{l}\omega_{ex} = 1.0\end{array}$ and therefore the natural eigenfrequency of the SDOF oscillator, and second an excitation with a frequency band containing three different frequencies $\begin{array}{l}\omega_{ex} = [0.75, 1.5, 2.0]\end{array}$ to analyze the performance of the different system and compare the sensitivities of the parameters of the tuned mass damper.

As in the previous examples, we chose the response function $\begin{array}{l}J = \int_0^T x^2(t)+v^2(t) dt\end{array}$ . But compared to the other examples, we are not interested in the movement of the mass (DoF) of the tuned mass damper and therefore only consider the values for the first mass of the single degree of freedom oscillator $\begin{array}{l}J = \int_0^T x_1^2(t) + v_1^2(t)dt\end{array}$ .

Single Frequency Excitation

For the single frequency excitation we get the following displacement for the three different systems:

As we can see, the system with Den Hartog optimized parameters performs better (oscillates less) than the other two. This behavior is expected after looking at the frequency plot of the system and seeing a much lower value at the frequency of interest $\begin{array}{l}\omega = 1.0\end{array}$ . Comparing the value of the response function J, we see that it is significantly lower for the Den Hartog system than the other two ( $\begin{array}{l}J_{DH} = 4.66\end{array}$ , $\begin{array}{l}J_{H2} = 10.30\end{array}$ and $\begin{array}{l}J_{SDOF} = 11.60\end{array}$ ). We solve the adjoint problem of the two systems with tuned mass damper to get the sensitivities of the stiffness and the damping value to interpret the results and validate the sensitivity analysis. Note that we are not interested in the sensitivity of the mass of the attached tuned mass damper $\begin{array}{l}m_d\end{array}$ as a larger mass always leads to smaller movement of the SDOF oscillator and is therefore not optimizable but chosen in the problem settings.

	J	$\begin{array}{l}c_d\end{array}$	$\begin{array}{l}k_d\end{array}$	$\begin{array}{l}\frac{dJ}{dc}\end{array}$	$\begin{array}{l}\frac{dJ}{dk}\end{array}$
Den Hartog	4.66	0.34	0.83	4.96	-5.29
H2	10.30	1.27	17.36	-0.02	0.01

To test the calculated sensitivities, we run the primal problem again with adjusted design parameters for the two different cases to compare a change in the respective response value. Since we are interested in reducing J (reduce oscillation of first mass), we follow the negative sensitivity of the parameters to expect an improvement:

	J	$\begin{array}{l}c_d\end{array}$	$\begin{array}{l}k_d\end{array}$	$\begin{array}{l}\frac{dJ}{dc}\end{array}$	$\begin{array}{l}\frac{dJ}{dk}\end{array}$
Den Hartog	4.66	0.34	0.83	4.96	-5.29
increase $\begin{array}{l}k_d\end{array}$	4.34	0.34	0.9	5.98	-3.4
reduce $\begin{array}{l}c_d\end{array}$	4.12	0.3	0.9	6.26	-3.81

	J	$\begin{array}{l}c_d\end{array}$	$\begin{array}{l}k_d\end{array}$	$\begin{array}{l}\frac{dJ}{dc}\end{array}$	$\begin{array}{l}\frac{dJ}{dk}\end{array}$
H2	10.3016	1.27	17.36	-0.02	0.01
increase $\begin{array}{l}c_d\end{array}$	10.3011	1.3	17.36	-0.02	0.01
reduce $\begin{array}{l}k_d\end{array}$ , increase $\begin{array}{l}c_d\end{array}$	10.2953	1.4	17	-0.02	0.01

As we can deduct from the tables, adjusting the parameters of the tuned mass damper according to their sensitivities, we are able to influence the movement of the first mass. This shows that the sensitivities calculated by the implemented adjoint method work. However, one can notice that they do not point at the same direction/ a global maximum for the two different cases of Den Hartog and H2 optimized parameters. This can be seen by looking at the damping value c for both systems ( $\begin{array}{l}0.34 \frac{Ns}{m}\end{array}$ & $\begin{array}{l}1.27 \frac{Ns}{m}\end{array}$ ) and their respective sensitivities. The sensitivities indicate that in the first case, a decrease of $\begin{array}{l}c_d\end{array}$ leads to a decrease of J ( $\begin{array}{l}\frac{dJ}{dc} > 0\end{array}$ ), whereas in the latter case an increase of $\begin{array}{l}c_d\end{array}$ leads to the decrease. This shows that the sensitivities do not point at a 'best' value for $\begin{array}{l}c_d\end{array}$ in between but they point in different directions/ local minima.

Multiple Excitation Frequencies

For multiple excitation frequencies, namely $\begin{array}{l}\omega_{ex} = [0.75, 1.5, 2.0]\end{array}$ , we get the following displacement plot for the three systems:

We notice that it's hard to determine the least oscillating system by observing the plot, as the response of the three systems is pretty similar. Comparing the value of the response function gives us a value to quantify the performance, but also indicates similar behavior, as the value for the three systems are pretty close together ( $\begin{array}{l}J_{SDOF} = 4.6\end{array}$ , $\begin{array}{l}J_{DH} = 4.67\end{array}$ and $\begin{array}{l}J_{H2} = 5.32\end{array}$ ). According to the value of the response function, the SDOF system actually performs better than the systems with tuned mass damper.

In the same way as for the single frequency excitation, we now focus on the two systems with a tuned mass damper to check their sensitivities and validate them.

	J	$\begin{array}{l}c_d\end{array}$	$\begin{array}{l}k_d\end{array}$	$\begin{array}{l}\frac{dJ}{dc}\end{array}$	$\begin{array}{l}\frac{dJ}{dk}\end{array}$
Den Hartog	4.6656	0.34	0.83	0.0115	3.2197
decrease $\begin{array}{l}c_d\end{array}$	4.6664	0.3	0.83	-0.0593	3.5786
increase $\begin{array}{l}k_d\end{array}$	4.9237	0.3	0.9	-0.7982	3.3769
decrease $\begin{array}{l}k_d\end{array}$	4.228	0.3	0.7	1.3049	3.1722

We notice that by decreasing $\begin{array}{l}c_d\end{array}$ in the first altering simulation, J increases instead of decreasing even though we followed the negative sensitivity (expect a smaller J). This can be explained by looking at the sensitivity $\begin{array}{l}\frac{dJ}{dc}\end{array}$ . In our original setup with parameter values according to Den Hartog, the sensitivity $\begin{array}{l}\frac{dJ}{dc} > 0\end{array}$ , but when adjusted the sensitivity flips its sign and becomes negative. Since the performance of the system is linked to both parameters - $\begin{array}{l}c_d\end{array}$ and $\begin{array}{l}k_d\end{array}$ - we can explain the worse behavior of the adjusted system by the dependence between $\begin{array}{l}c_d\end{array}$ and $\begin{array}{l}k_d\end{array}$ . By decreasing $\begin{array}{l}c_d\end{array}$ and keeping $\begin{array}{l}k_d\end{array}$ constant, we crossed the point, where decreasing the damping constant has a reducing factor on the sensitivity and went too far, which led to an increase of the response. It can also be seen, that for our second trial, we followed the positive sensitivity of the stiffness parameter and increased $\begin{array}{l}k_d\end{array}$ . This led to an increase of J, which is the behavior we expected based on the positive sensitivity $\begin{array}{l}\frac{dJ}{dk} > 0\end{array}$ .

	J	$\begin{array}{l}c_d\end{array}$	$\begin{array}{l}k_d\end{array}$	$\begin{array}{l}\frac{dJ}{dc}\end{array}$	$\begin{array}{l}\frac{dJ}{dk}\end{array}$
H2	5.3211	1.27	17.36	-0.0025	0.0006
increase $\begin{array}{l}c_d\end{array}$	5.3208	1.4	17.36	-0.0024	0.0006
decrease $\begin{array}{l}k_d\end{array}$	5.3206	1.4	17	-0.0025	0.0006
decrease $\begin{array}{l}k_d\end{array}$	5.3202	1.4	16.5	-0.0026	0.0007

For the H2 optimized parameters we show that the system behaves as expected, but changing the parameter values only has a little effect on the response of the system. This can also be seen by the value of the sensitivities, which are very close to 0.

How to Use the Implementation (User's Guide)

Assuming that user knows, how to set up a model and solve the primal problem in bm_mfem - namely creating/importing a model, defining the simulation time T and the time step dt (also in the model by model.setDt(dt)/.setEndTime(T)!), selecting a solver and solving the problem - the first step is to choose a response function J. As part of this project, we implemented two different response functions: the one we used in the examples $\begin{array}{l}J = \int_0^T x^2(t)+v^2(t)dt\end{array}$ as well as the energy dissipation of the system $\begin{array}{l}J = \int_0^T - v^T(t)*C*v(t) dt\end{array}$ . Note that interpreting the sensitivities $\begin{array}{l}\frac{dJ}{ds_i}\end{array}$ depends on the choice of the respected response function J.

In bm_mfem, user can simply create an instance of either response function by calling the respective constructor and giving the model as input variable:

User can get the actual value of J by simply calling the member function 'getValue()' of the instances (i.e. response.getValue()), but that is not necessary for the sensitivity calculation.

After choosing a response function and creating an instance of it, we are ready to solve the adjoint problem. To solve the adjoint problem, we create an instance of a Newmark solver especially implemented to solve for $\begin{array}{l}\Lambda\end{array}$ . This can be done by simply calling the constructor of 'NewmarkSolvingStrategyAdjointProblem' and giving the defined model and the instance of the response as input variables:

Since we set $\begin{array}{l}dt = d\tau\end{array}$ , meaning that we do the same time step to solve the adjoint problem as for solving the primal problem, we have to create a loop running through the steps (starting at 2 until counter = steps + 1, steps = endTime/dt). In the loop, we simply call the member function 'solve(counter)' and increase our counter by 1:

After solving the adjoint problem, we can create an instance of the 'Sensitivity' class with our model and response function as input parameters.

The sensitivity class has two callable functions: first we can either call 'getSensitivity(element, parameter)' with a specific element and a parameter name as parameters (i.e. sensitivity.getSensitivity(springElement,"ELEMENTAL_STIFFNESS")). This function will return the sensitivity of the respective parameter with respect to our response function J.

The second member function 'getAllSensitivities()' calculates the sensitivity of all elements and all respective parameters. It returns an array containing the element id, the parameter name and the value of the sensitivity. With the for loop, user is able to display the sensitivities in the command line.

References

[1] Asami, T., Nishihara, O., and Baz, A. M. (March 26, 2002). "Analytical Solutions to H∞ and H2 Optimization of Dynamic Vibration Absorbers Attached to Damped Linear Systems ." ASME. J. Vib. Acoust. April 2002; 124(2): 284–295. https://doi.org/10.1115/1.1456458

[2] Emma Smith Zbarsky (2023). Mass-Spring-Damper Systems (https://github.com/MathWorks-Teaching-Resources/Mass-Spring-Damper-Systems/releases/tag/v1.0.1), GitHub. Retrieved February 9, 2023.

[3] Fernandez, F., & Tortorelli, D. A. (2018). Semi-analytical sensitivity analysis for nonlinear transient problems. Structural and Multidisciplinary Optimization, 58(6), 2387-2410.

Appendix

Derivation of the adjoint problem

The equation of motion $\begin{array}{l}R\end{array}$ is differentiated with respect to $\begin{array}{l}s_i\end{array}$ .

$\begin{array}{l}\displaystyle \frac{\partial{R}}{\partial{s_i}} = \frac{\partial{R}}{\partial{\ddot{u}}} \frac{\partial{\ddot{u}} }{\partial{s_i}} + \frac{\partial{R}}{\partial{\dot{u}}} \frac{\partial{\dot{u}} }{\partial{s_i}} + \frac{\partial{R}}{\partial{u}} \frac{\partial {u}}{\partial{s_i}} = 0\end{array}$

where

$\begin{array}{l}\displaystyle \frac{\partial{R}}{\partial{\ddot{u}}} = M, \frac{\partial{R}}{\partial{\dot{u}}} = C, \frac{\partial{R}}{\partial{u}} = K\end{array}$

Multiply $\begin{array}{l}\frac{dR}{ds_i}\end{array}$ with Lagrange parameters $\begin{array}{l}\lambda\end{array}$ and add to the sensitivity, results in:

$\begin{array}{l}\displaystyle \frac{d{J}}{d{s_i}} = \frac{d{J}}{d{s_i}} + \int_0^T{\lambda}^T \cdot \frac{dR}{ds_i}dt\end{array}$

Substituting $\begin{array}{l}\frac{dR}{ds_i}\end{array}$ into $\begin{array}{l}\frac{dJ}{ds_i}\end{array}$ and reformulate:

$\begin{array}{l}\displaystyle \frac{dJ}{ds_i} = \int_0^T\left ( (\frac{\partial{j}}{\partial{s_i}} + \lambda^T\frac{\partial{R}}{\partial{s_i}}) + ((\lambda^TK + \frac{\partial{j}}{\partial{u}}) \frac{du}{ds_i} + (\lambda^TC + \frac{\partial{j}}{\partial{\dot{u}}})\frac{d{\dot{u}}}{ds_i} + \lambda^TM\frac{d{\ddot{u}}}{ds_i})\right )dt\end{array}$

where $\begin{array}{l}\lambda\end{array}$ is the arbitrary adjoint vector. Integrating by parts to eliminate the high derivatives :

$\begin{array}{l}\int_0^T(\lambda^TC + \frac{\partial j}{\partial{\dot{u}}})\frac{d{\dot{u}}}{ds_i}}dt = \left [(\lambda^TC + \frac{\partial{j}}{\partial{s_i}})\frac{du}{ds_i} \right ]_0^T - \int_0^T(\dot{\lambda}^TC + \frac{d}{dt}\frac{\partial j}{\partial{\dot{u}}})\frac{du}{ds_i}dt \\ \int_0^T\lambda^TM\frac{d\ddot{u}}{ds_i}dt = \left [ \lambda^TM\frac{d\dot{u}}{ds_i} \right ]_0^T - \left [ \dot{\lambda}^TM\frac{du}{ds_i} \right ]_0^T + \int_0^T \frac{d}{dt}\dot{\lambda}^TM\frac{du}{ds_i}dt\end{array}$

Rearranging the equation, we get:

$\begin{array}{l}\frac{dJ}{ds_i} = \int_0^T\left ( \frac{\partial{j}}{\partial{s_i}} + \lambda^T\frac{\partial{R}}{\partial{s_i}} \right )dt + \left [(\lambda^TC + \frac{\partial{j}}{\partial{s_i}})\frac{du}{ds_i} \right ]_0^T + \left [ \lambda^TM\frac{d\dot{u}}{ds_i} \right ]_0^T - \left [\dot{\lambda}^TM\frac{du}{ds_i} \right ]_0^T \\ { } + \int_0^T \left (\frac{d}{dt}\dot{\lambda}^TM - \dot{\lambda}^TC + \lambda^TK + \frac{\partial{j}}{\partial{u}} - \frac{d}{dt}\frac{\partial j}{\partial{\dot{u}}} \right )\frac{du}{ds_i}dt \\\end{array}$

$\begin{array}{l}\frac{dJ}{ds_i} = \int_0^T\left ( \frac{\partial{j}}{\partial{s_i}} + \lambda^T\frac{\partial{R}}{\partial{s_i}} \right )dt - \left [(\lambda^TC + \frac{\partial{j}}{\partial{s_i}} - \dot{\lambda}^TM)\frac{du}{ds_i}\right]_{t=0} - \left[\lambda^TM\frac{d\dot{u}}{ds_i}\right]_{t=0} \\ { }+ \int_0^T \left (\frac{d}{dt}\dot{\lambda}^TM - \dot{\lambda}^TC + \lambda^TK + \frac{\partial{j}}{\partial{u}} - \frac{d}{dt}\frac{\partial j}{\partial{\dot{u}}} \right )\frac{du}{ds_i}dt + \left [(\lambda^TC + \frac{\partial{j}}{\partial{s_i}} - \dot{\lambda}^TM)\frac{du}{ds_i} \right ]_{t=T} + \left [ \lambda^TM\frac{d\dot{u}}{ds_i} \right ]_{t=T}\end{array}$

Assuming the initial conditions are independent of s_i → $\begin{array}{l}\frac{du}{ds_i}|_{t=0} = 0, \frac{d\dot{u}}{ds_i}|_{t=0} = 0\end{array}$ and setting $\begin{array}{l}\ddot{\lambda}^TM - \dot{\lambda}^TC + \lambda^TK + \frac{\partial{j}}{\partial{u}} - \frac{d}{dt}\frac{\partial j}{\partial{\dot{u}}} = 0\end{array}$ to get rid of $\begin{array}{l}\frac{du}{ds_i}\end{array}$ , we obtain a system of a second-order differential equation with two conditions:

$\begin{array}{l}M^T\ddot\lambda - C^T\dot{\lambda} + K^T\lambda = - \left (\frac{\partial j}{\partial u} \right)^T + \dot{u}^T + \left ( \frac{\partial^2{j}}{\partial u \partial{\dot{u}}} \right) ^T + \ddot{u}^T \left (\frac{\partial^2 j}{\partial{\dot{u^2}}} \right )^T \\ \lambda(T) = 0 \\ M^T\dot{\lambda}(T) = -\left (\frac{\partial j}{\partial \dot{u}} \right)_{t=T}\end{array}$

Since the conditions are given at time t = T instead of t = 0, we need to solve the equation backwards in time. Therefore, we introduce the substitution:

$\begin{array}{l}\tau = T - t \\ \lambda(t) = \Lambda(T-t) = \Lambda(\tau), \dot{\lambda}(t) = - \dot{\Lambda}(\tau), \ddot\lambda(t) = \ddot\Lambda(\tau)\end{array}$

Thus, the adjoint problem reads:

$\begin{array}{l}M^T\ddot\Lambda + C^T\dot{\Lambda} + K^T\Lambda = - \left (\frac{\partial j}{\partial u} \right)^T + \dot{u}^T + \left ( \frac{\partial^2{j}}{\partial u \partial{\dot{u}}} \right) ^T + \ddot{u}^T \left (\frac{\partial^2 j}{\partial{\dot{u^2}}} \right )^T \\ \Lambda(0) = 0 \\ M^T\dot{\Lambda}(0) = -\left (\frac{\partial j}{\partial \dot{u}} \right)_{t=T}\end{array}$

Sensitivity analysis methods

Finite Difference method

The forward finite difference method approximates the derivatives of the response function $\begin{array}{l}J\end{array}$ using Taylor series expansion:

$\begin{array}{l}\displaystyle \frac{dJ}{ds_i} \approx \frac{J(u(s + \in e_i), \dot{u}(s + \in e_i), s + \in e_i, t) - J(s, u, \dot{u})}{\in}\end{array}$

where $\begin{array}{l}e_i = [0,0,...,1,...,0,0]^T\end{array}$ is the unit vector of component $\begin{array}{l}i\end{array}$ , and $\begin{array}{l}\in\end{array}$ is the perturbation.

The finite difference method has the advantages of easy implementation and good order of convergence. However, as seen above, it is necessary to compute $\begin{array}{l}J(u + \in e_i)\end{array}$ for each design variable $\begin{array}{l}s_i\end{array}$ by modifying the finite element model, and so makes the whole process very computationally costly [3].

Direct differentiation method

Recalling the equation of motion $\begin{array}{l}𝑹(𝒔, 𝒖, 𝒖 ̇, 𝒖 ̈,𝑡)=𝑴(𝒔) 𝒖 ̈+𝑪(𝒔) 𝒖 ̇+𝑲(𝒔)𝒖−𝑭(𝒕) \equiv 0\end{array}$

Differentiate both sides of $\begin{array}{l}R\end{array}$ with respect to $\begin{array}{l}s_i\end{array}$ , we have:

$\begin{array}{l}\frac{\partial{R}}{\partial{s_i}} = M \frac{\partial{\ddot{u}} }{\partial{s_i}} + C \frac{\partial{\dot{u}} }{\partial{s_i}} + K \frac{\partial {u}}{\partial{s_i}} \\ \frac{\partial u(0)}{\partial s_{i}} = \frac{\partial {u^0}}{\partial s_{i}} \\ \frac{\partial \dot{u}(0)}{\partial s_{i}} = \frac{\partial \dot{u}^{0}}{\partial s_{i}}\end{array}$

With the direct solver used for the primal problem above, $\begin{array}{l}\frac{\partial{\ddot{u}} }{\partial{s_i}}, \frac{\partial{\dot{u}} }{\partial{s_i}}, \frac{\partial {u}}{\partial{s_i}}\end{array}$ can be obtained. Then, we plug these implicit derivatives into the equation of $\begin{array}{l}\frac{dJ}{ds_i}\end{array}$ to evaluate the sensitivity.

Clearly, this method is computationally inefficient especially for a large amount of design variables, as it solves for the implicit derivatives using $\begin{array}{l}\frac{dR}{ds_i}\end{array}$ for each $\begin{array}{l}s_i\end{array}$ regardless of the number of design parameters. On the other hand, this method provides numerically exact sensitivities [3].

Seitenhierarchie