Parameter and Topology Optimization of Structures in the Frequency Domain under Nevanlinna–Pick Interpolation Constraints

Abstract

The design optimization of structures can be conducted in either the time domain or the frequency domain. The frequency domain approach is advantageous compared to its time domain counterpart, especially if the degree of freedom is large, the objectives and/or constraints are formulated in the frequency domain, or the structure is subject to random loading. In this paper, an attempt is undertaken to obtain feasible optimal solutions by implementing the Nevanlinna–Pick (NP) interpolation theory across multi-objective structural optimization problems in the frequency domain. The NP equations introduce a trade-off that originates from the interpolation theory for complex variables. According to the NP theory, a complex function cannot have an independent amplitude from its derivative at a certain frequency. Consequently, the frequency response of a physical system cannot be shaped arbitrarily at discrete frequencies. Our objectives within this paper include calculating the weight, natural frequency, fatigue life, frequency domain response, and its derivative. To illustrate our claims, sample parameter and topology optimization problems were formulated and solved, both with and without the NP constraints. It was found that the inclusion of NP constraints induced a considerable improvement in the optimal solutions, while also causing the convergence to the optimal solution to become smoother.

1. Introduction

It is no longer sufficient to design a technical system that fulfills its function; it has become important to achieve the highest performance possible with the smallest quantity of material, i.e., optimization. Many studies into structural optimization in the past few decades have responded to this need. Although the optimization of structures can be realized in either the time or frequency domain, frequency domain analysis techniques are widely employed in engineering structures on account of their advantages over their time domain counterparts. For instance, it is more challenging to obtain the time domain response of a system analytically than that of the frequency domain, especially in the case of random loading. In addition, the difficulty increases as the order of the system increases. Furthermore, the sensitivity of a system to parameter variations can be measured or evaluated more conveniently in the frequency domain [1]. The frequency domain provides great convenience, enabling analyses to be performed over a wide frequency range under random loading. Consequently, the optimum design of structures, especially those exposed to harmonic and random loading, is preferably studied in the frequency domain. The frequency domain approach was first predominantly deployed in control theory and has since been used in the analysis of mechanical systems due to its convenience. In chronological order, the relevant texts include [2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18].

For the purpose of optimizing a mechanical system in the frequency domain, various problem formulations have been proposed in the literature. These can be categorized according to their objective functions as follows: response (vibration) control, eigenfrequency maximization, and weight minimization. Although vibration control can be accomplished using an active control approach, deploying a passive control approach during the design phase of mechanical systems is more economical [19]. Vibration control is mainly achieved by reducing the structural response by either directly defining response amplitudes as objective functions [20] or by using transfer function norms (i.e., ${\mathit{H}}_2$ and ${\mathit{H}}_{\mathit{\infty }}$ norms) as objective functions [6,17,21]. Additionally, the root mean square (RMS) or weighted RMS of structural responses such as acceleration and displacement responses are other objective functions that have been used in several studies [4,15,22]. Another approach to vibration control is to formulate an objective function by utilizing different compliance expressions for the purpose of optimizing the frequency response. In [23], an attempt was made to minimize the absolute value (modulus) of the mean compliance and its integration, while in [9], dynamic compliance was used to optimize the dynamic response of the structure. In another study, the maximization of the imaginary part of the complex dynamic compliance of damping materials was investigated in order to reduce the resonance peak response [13]. Moreover, research has utilized both mean compliance formulations [24] and new approaches to compliance minimizations for topology optimization [25].

Additionally, a considerable amount of study has been conducted on the multi-objective optimization of structures over the past few years. This is of especial importance given that mechanical systems must be designed in a manner that considers multiple phenomena [2,22,26,27]. The objective functions, which constitute the cost function, can conflict with each other, as in the cases of weight and displacement [28,29,30] or weight and fatigue/damage [16]. Cost functions with multiple objectives, such as weight, damage, and response amplitude, are used for the optimization problems within this study. Three main methods are reported for solving these multi-objective optimization problems: a priori, a posteriori, and interactive preference methods [31]. In this study, the weighting method, which is an a posteriori method, was employed to solve multi-objective optimization problems.

Furthermore, some specific objective functions were used to optimize structural components, such as the integral of the weighted relative frequency response [5], the frequency response deviation [32], the weighted sum of the difference between the current and initial thrust coefficients and the torque coefficient [10], and the simultaneous minimization of stress and stiffness for topology optimization [33].

Although there are various objective functions involved in the structural optimization problems mentioned so far, the ultimate goal, or at least a constraint that cannot be ignored, in designing and optimizing a mechanical system is its durability. For this purpose, fatigue life has been considered a constraint [34,35,36,37,38,39,40,41] or an objective function in different forms in the time domain optimization of structures in various studies. These expressions have included damage [42], a specific function expression [43], multiple functions derived from the statistical properties of fatigue life (standard deviation and mean) [44], and impact damage area [27]. Conversely, it is noteworthy that only several studies considering fatigue life in frequency domain optimization [12,16,45] exist in the literature. In this respect, this paper contributes to frequency domain optimization research with numerical simulations defining objective functions, including the fatigue damage to a structure.

Frequency domain optimization approaches use both traditional and modern techniques. For instance, nonlinear [2] and linear programming algorithms [4] and methods of moving asymptotes [9,13,46] exemplify the traditional optimization techniques used for structural optimization. In addition, modern optimization techniques, such as genetic (GA) and evolutionary algorithms (EA) [6,7,16,18], particle swarm optimization (PSO) [8], harmony search [14], and colliding bodies optimization [47], are also used to optimize the frequency domain. In addition to traditional and modern optimization methods, a combination of both procedures, namely a hybrid approach, may also be deployed to achieve more efficient and accurate calculations.

The optimum solutions obtained in the frequency domain exhibit difficulties when applied to physical systems. There remains a need to verify the accuracy of frequency domain results via their comparison with time domain results. Despite the aforementioned advantages of frequency domain analysis, it was demonstrated in [16] that frequency domain analysis may yield infeasible results.

In short, while a number of studies exist on structural optimization in the frequency domain, trade-off equations that limit the achievable best solution have not been investigated in the physical domain. Rectifying this issue constitutes the focus of this paper and is achieved using the Nevanlinna–Pick (NP) interpolation theory [48]. Although NP interpolation problems have been studied in control theory [49,50], they have not been investigated for structural optimization problems. Motivated by this fact, NP constraints are implemented in frequency domain optimization formulations to develop an approach to solving constrained optimization problems by considering reachability constraints. The principal contribution of this study is that the NP equations are imposed as nonlinear constraints, frequency response amplitudes are calculated during optimization iterations, and the numerical results are compared to those obtained without the NP constraint. The design variables that cannot satisfy the NP constraints are eliminated, and the feasible region in the parameter space is used in the search for optimal solutions. In this paper, the development of a formulation that is applicable to multiple load cases and multi-objective optimization problems, while not experiencing convergence difficulties or instability in minimizing stress levels as a result of minimizing fatigue damage indices, is of particular importance. The optimization formulation in this paper does not deploy conventional compliance minimization, as this has been reported to be ineffective in some cases [25]. Since the NP equations include frequency response amplitudes, compliance is indirectly included in the cost functions. In addition, we reduce stress levels to the greatest extent possible through the frequency domain fatigue damage minimization cost function term. It is shown that using the NP equations in frequency domain optimization problems improves the optimality of the solutions by providing smaller cost functions at smaller iteration numbers with smoother convergence. This paper is organized as follows: First, the mathematical formulations of the associated optimization problem are given. Next, numerical results are presented to verify the theoretical derivations. Finally, conclusions are drawn.

2. Implementation of the Nevanlinna–Pick Constraint into Structural Dynamics

According to the boundary NP interpolation theorem [48], the existence of a solution, which is a rational function mapping from the open right half plane ${\pi }^{+}$ into an open unit disk $D$, is related to the positive definiteness of the following Pick matrix ${\mathsf{\Lambda }}_{lm}$

$${Ʌ}_{lm}=\left\{\begin{array}{c} \frac{1-{\overline{f}}_l{f}_m}{{\overline{s}}_l+s_m} , if l\ne m\\ \left|f^{\prime }\left(s_l\right)\right| , if l=m\end{array}\right.$$

(1)

where $f_k$ is the given complex function defined by

$$\left|f_k\right|=1$$

(2)

and $s_k=i{\omega }_k$ is the complex Laplace transform variable, $i$ is the imaginary unit, and ${\omega }_k$ is the excitation frequency. In summary, a physical system cannot display an arbitrary frequency response function (FRF) $f\left(i\omega \right)$ and derivative $f^{\prime }\left(i\omega \right)$ at discrete frequencies ${\omega }_k$. The above equations are implemented into finite element equations for sample structural problems, represented by the following equation:

$$M\ddot{x}+C\dot{x}+Kx=F$$

(3)

where M, C, and K are, respectively, mass, damping, and stiffness matrices. F is the applied force vector, x is the displacement vector, and time rates are denoted by superposed dots. If the Laplace transform of Equation (3) is taken, then we obtain the following equation in terms of the Laplace transformed variables:

$$\left(M{s}^2+Cs+K\right)X\left(s\right)=F\left(s\right)$$

(4)

where s is the Laplace transform variable. The corresponding system transfer function is expressed as follows:

$$G\left(s\right)=\frac{F\left(s\right)}{X\left(s\right)}={\left(M{s}^2+Cs+K\right)}^{-1}$$

(5)

and $f_k=G\left({\omega }_k\right)$. For a harmonic force with the frequency $\omega$, i.e., $F\left(t\right)=F_0{e}^{i\omega t}$, the steady-state response of a structure can be expressed as follows [1]:

$$x_p\left(t\right)=x_0\left|G\left(i\omega \right)\right|e^{i(\omega t+\varnothing )}$$

(6)

where $x_p\left(t\right)$ denotes the particular solution, $\left|G\left(i\omega \right)\right|$ is the magnitude of the transfer function of the structure, and $\varnothing$ is the phase angle. The NP equations are thus related to the steady-state response of the structure to a harmonic force.

To derive insight into imposing the NP constraint equations onto structural equations, the physical meaning of the terms in Equation (1) is pursued. To this end, following Equation (5), it can be written that $X\left(i\omega \right)=G\left(i\omega \right)F\left(i\omega \right)$. Suppose that the forcing term in Equation (3) is an impulsive force with the following harmonic expansion: $F\left(t\right)=F_0\delta \left(t\right)=F_0{\sum }_{i=1}^{\infty }{A}_i\sin \left({\omega }_{i}t\right),$ where $F_0$ is a constant vector, $\delta \left(t\right)$ is the impulsive force, and ${\omega }_i$ the are harmonic frequencies. Then, the impulse response $g\left(t\right)$ can be calculated by using the Bode gain and phase equations [1]:

$$g\left(t\right)={\sum }_{i=1}^{\infty }{F}_0\left|G\left(i{\omega }_i\right)\right|A_i\sin \left({\omega }_{i}t+{\phi }_i\right)$$

(7)

where ${\phi }_i=Arg\left(G\left({\omega }_i\right)\right)$ is the phase angle and $G\left(i\omega \right)$ is found using the Laplace transform of $g\left(t\right)$ as follows: $G\left(i\omega \right)=\mathcal{L}\left\{g\left(t\right)\right\}$. On the other hand, since $G^{\prime }\left(i\omega \right)=-\mathcal{L}\left\{tg\left(t\right)\right\}$ and $\left|G^{\prime }\left(i\omega \right)\right|=\left|\mathcal{L}\left\{tg\left(t\right)\right\}\right|$, e.g., see [1], we derive

$$tg\left(t\right)={\sum }_{i=1}^{\infty }{F}_0\left|G^{\prime }\left({\omega }_i\right)\right|A_i\sin \left({\omega }_{i}t+{{\phi }_i}^{*}\right)$$

(8)

where ${{\phi }_i}^{*}=Arg\left(G^{\prime }\left({\omega }_i\right)\right)$ are phase angles. In summary, Equations (7) and (8) indicate that the frequency domain quantities of $\left|G\left({\omega }_i\right)\right|$ and $\left|G^{\prime }\left({\omega }_i\right)\right|$ are related to the impulse response magnitude and time-weighted impulse response magnitude, respectively. As a consequence, the magnitudes of $g\left(t\right)$ and $tg\left(t\right)$ cannot be independent due to Equation (1). Consequently, imposing the NP constraint Equation (1) in optimization costs assumes the dependence between $g\left(t\right)$ and $tg\left(t\right)$. Therefore, the $\left|G\left({\omega }_i\right)\right|$ and $\left|G^{\prime }\left({\omega }_i\right)\right|$ terms were independently imposed in our optimization formulations in the cost functions. Moreover, the NP constraint Equation (1) was also imposed in our numerical examples to establish the dependence between the $\left|G\left({\omega }_i\right)\right|$ and $\left|G^{\prime }\left({\omega }_i\right)\right|$ terms.

This study demonstrates how the NP interpolation theory can be implemented in frequency domain structural optimization formulations. Using this method, an achievable performance in the frequency domain can be determined such that the optimization algorithm will not attempt to search for a physically infeasible objective. Consequently, the NP constraints are implemented in the optimization problem as follows: as the optimization algorithm updates the design parameter vector during iterations, it checks whether the corresponding Pick matrix is positive definite. For example, if the Pick matrix in Equation (1) does not satisfy the NP constraints due to the vibration amplitudes at selected frequencies, the design parameter vector is redefined such that the Pick matrix is positive definite. Therefore, all design parameter vectors should satisfy the positive definite Pick matrix condition (namely, the feasibility condition). In summary, the proposed approach searches for an optimal solution to the optimization formulations by adding the positive definiteness condition of the Pick matrix as a nonlinear constraint. Thus, during the updating step of the design parameter vector, additional sensitivity calculations are required due to the nonlinear NP constraints stemming from the positive definiteness of the Pick matrix.

Practical Considerations

When the NP constraint was implemented into the sample structural systems, we were required to consider certain properties of the Pick matrix. We observed that the positive definiteness of the Pick matrix changed at certain frequencies in a system. As an illustrative example, consider a 1-DOF mass-spring-damper system with a mass of 0.1 kg, a stiffness of 0.1 N/m, and a damping constant of 0.01 N/m/s. Figure 1 shows the corresponding Bode gain and phase plots of the system.

In numerical simulations, the eigenvalue test was used to determine the positive definiteness of the Pick matrix in Equation (1). We found that the Pick matrix was negative definite at frequencies around the natural frequency, while it was positive definite at all other frequencies. In other words, the positive definiteness of the Pick matrix was violated when the frequency range included the natural frequency of the system, which, in this example, was 1 rad/s. It was also observed that both the frequency range in which the Pick matrix was evaluated and the damping constant affected the positive definiteness of the Pick matrix. The positive definiteness of the Pick matrix was violated with damping constants above almost 0.03 N/m/s, in this example.

When the Bode plots corresponding to different damping constants were examined, it was found that the change in the slope of the Bode gain curves resulted in a negative definite Pick matrix. This occurred because of the flattening of the response curve around the peak of the frequency response curve. The Bode gain curve flattened as the damping constant increased, which is associated with a decrease in the derivative of the transfer function. Figure 2 clearly shows that the peak frequency of the response shifts to the left as the damping increases, but that the natural frequency of the system remains constant. As a result, the positive definiteness of the Pick matrix was affected by both the frequency range (related to the natural frequency) and the damping constant.

In short, the positive definiteness of the Pick matrix is influenced by the peak frequency of the response and the damping constant of the system. The peak frequency of the response almost coincides with the natural frequency at lower damping constants, while it shifts significantly away from the natural frequency at higher damping constants. In addition, the Bode gain plot flattens as the damping constant increases. As the damping constant increases, the frequency range in which the Pick matrix is negative definite enlarges around the peak frequency of the response. For this reason, the frequency range under consideration during optimization is divided into two regions: first the left side of the peak frequency of the response, and, second, the right side of the peak frequency of the response. This excludes the non-smooth part of the curve that renders the Pick matrix negative definite.

3. Numerical Examples

In this section, solutions to parameter and topology optimization problems are presented to illustrate the efficiency of the NP constraint in obtaining optimum solutions. All programs were written using Matlab (2021b) software. Figure 3 shows the flowchart of the optimization algorithm utilized to solve the example problems. After calculating the objective functions and constraints, the optimization algorithm checks for convergence and updates the design variables accordingly. Two optimization methods were used in the simulations: while the Matlab function ‘fmincon’ was used to solve the parameter optimization problems by selecting the interior point optimization algorithm for the three examples in this section, the PSO method was also used in the truss system and topology optimization problems.

3.1. Pareto Optimal Solutions for a Mass-Spring-Damper System

In the first example, the aim was to illustrate the utility of the NP constraint in finding the optimum solution. Considering the three objective function terms, the multi-objective optimization problem was formulated as follows:

$$\underset{ m, c, k}{\mathrm{m}\mathrm{i}\mathrm{n}} F={\sum }_{i=1}^n\left\{\alpha \left|G\left({\omega }_i\right)\right|+\beta \left|G^{\prime }\left({\omega }_i\right)\right|+\gamma (1/{\omega }_n)\right\}$$

(9)

$$Subject to \left\{\begin{array}{c}-{Ʌ}_{ij}+\epsilon \le 0\\ {\mathit{h}}_{lb}{\le \mathit{h}}_n\le {\mathit{h}}_{ub}\end{array}\right.$$

(10)

where the weightings were defined as $0\le \alpha \le 1$, $0\le \beta \le 1$, and $\alpha +\beta +\gamma$ = 1. ${Ʌ}_{ij}$ denotes the Pick matrix given in Equation (1), and its dimension is determined by the number of discrete frequency values, ${\omega }_i$. $\mathsf{\epsilon }$ is a very small value ($\mathsf{\epsilon }\ll 1$) (here, $\epsilon ={10}^{-3}$) used to eliminate spurious oscillations in the positive definiteness check, and $\mathit{h}$ is the vector of the structural parameters, i.e., $\mathit{h}=\left\{m c k\right\}$, in this case bounded by the lower and upper limit values ${\mathit{h}}_{lb}$ and ${\mathit{h}}_{ub}$, respectively. The input parameters are listed in

Table 1. Input parameters for optimization problems.

Parameter	Value
${ h}_{lb}$	$\left\{0.01, 0.0001, 0.1\right\}$
${ h}_{ub}$	$\left\{10, 0.1, 50\right\}$
${ h}_{initial}$	$\left\{0.01, 0.0001, 0.1\right\}$
$∆\omega$	$0.1$
${\omega }_{vec1}$	$0.001:∆\omega :({\omega }_d-{\omega }_d/10)$
${\omega }_{vec2}$	$\left({\omega }_d+{\omega }_d/10\right):∆\omega :5{\omega }_n$

, where${ h}_{initial}$ is the initial parameter vector, $∆\omega$ is the frequency increment used in the frequency domain calculations, ${\omega }_{vec1}$ and ${\omega }_{vec2}$ are frequency vectors used for calculation when $\left|G\left({\omega }_i\right)\right|$, and $\left|G^{\prime }\left({\omega }_i\right)\right|$, and ${\omega }_d$ is the damped natural frequency. The lower bound vector ${\mathit{h}}_{lb}$ and the upper bound vector ${\mathit{h}}_{ub}$ were selected such that the mass is in a physically realizable range for experimental setups, the damping coefficient was in the range for metal materials and the stiffness coefficient was in a practical range. The parameter vector range was kept intentionally large to make the differences between the solutions with and without the NP constraint more clear. The interior point method in Matlab was utilized in this example. Note that the terms in Equation (9) are the magnitude of the frequency response $\left|G\left({\omega }_i\right)\right|$, the roll rate of the frequency response $\left|G^{\prime }\left({\omega }_i\right)\right|$, and the reciprocal of the natural frequency $1/{\omega }_n$. Due to Equations (7) and (8), the minimization of $\left|G\left({\omega }_i\right)\right|$ and $\left|G^{\prime }\left({\omega }_i\right)\right|$ provides, respectively, lower displacement and slope $\left|G^{\prime }\left({\omega }_i\right)\right|$ around natural frequencies, while the minimization of $1/{\omega }_n$ provides the maximization of the natural frequency. In short, the minimization of $\left|G\left({\omega }_i\right)\right|$ and $\left|G^{\prime }\left({\omega }_i\right)\right|$ is equivalent to the minimization of compliance and its rate in the frequency domain.

Figure 4 illustrates the Pareto optimal solution sets obtained by increasing the values of $\alpha$, $\beta$, and $\gamma$ by a step size of 0.1, both with and without the NP constraint equations. It can be observed in Figure 4 that the solutions obtained without the NP constraint show more scattering. In the end, using the NP constraints resulted in a total cost function of 108.11, while without using the NP constraints the total cost function was 118.6.

Table 2. Minimum values in Pareto optimal sets for the three cost functions.

	${\mathit{F}}_1={\displaystyle \sum _{\mathit{i}=1}^{\mathit{n}}}\left\{\begin{array}{c}\mathit{\alpha }\left|\mathit{G}\left({\mathit{\omega }}_{\mathit{i}}\right)\right|+\\ \mathit{\beta }\left|{\mathit{G}}^{\mathbf{\prime }}\left({\mathit{\omega }}_{\mathit{i}}\right)\right|+\mathit{\gamma }(1/{\mathit{\omega }}_{\mathit{n}})\end{array}\right\}$		${\mathit{F}}_2={\displaystyle \sum _{\mathit{i}=1}^{\mathit{n}}}\left\{\begin{array}{c}10\mathit{\alpha }\left|\mathit{G}\left({\mathit{\omega }}_{\mathit{i}}\right)\right|+\\ \mathit{\beta }\left|{\mathit{G}}^{\mathbf{\prime }}\left({\mathit{\omega }}_{\mathit{i}}\right)\right|+\mathit{\gamma }(1/{\mathit{\omega }}_{\mathit{n}})\end{array}\right\}$		${\mathit{F}}_3={\displaystyle \sum _{\mathit{i}=1}^{\mathit{n}}}\left\{\begin{array}{c}\mathit{\alpha }\left|\mathit{G}\left({\mathit{\omega }}_{\mathit{i}}\right)\right|+\\ 10\mathit{\beta }\left|{\mathit{G}}^{\mathbf{\prime }}\left({\mathit{\omega }}_{\mathit{i}}\right)\right|+\mathit{\gamma }(1/{\mathit{\omega }}_{\mathit{n}})\end{array}\right\}$
	w/o NP	w/NP	w/o NP	w/NP	w/o NP	w/NP
$min\left(F_i\right)$	118.6	108.11	570.14	224.98	708.21	595.91
${\omega }_n$	1.1024	0.889	1.051	0.333	1.111	1.3932
$\left|G\left(\omega \right)\right|$	52	42	50	15	53	67
$\left|G^{\prime }\left(\omega \right)\right|$	65.696	64.987	69.189	71.975	65.431	52.819
No. of iterations	3693	4539	2508	4942	3396	4700

also shows the optimum solutions calculated by altering the weight of each objective function. We observed that introducing a greater penalty for an objective resulted in lower values for that objective. Table 2 compares the results of three cost functions, which differ from each other due to the scaling factor of each objective: $F_1={\sum }_{i=1}^n\left\{\alpha \left|G\left({\omega }_i\right)\right|+\beta \left|G^{\prime }\left({\omega }_i\right)\right|+\gamma (1/{\omega }_n\right\}$, $F_2={\sum }_{i=1}^n\left\{10\alpha \left|G\left({\omega }_i\right)\right|+\beta \left|G^{\prime }\left({\omega }_i\right)\right|+\gamma (1/{\omega }_n\right\}$, and $F_3={\sum }_{i=1}^n\left\{\alpha \left|G\left({\omega }_i\right)\right|+10\beta \left|G^{\prime }\left({\omega }_i\right)\right|+\gamma (1/{\omega }_n\right\}$. Note that increasing the penalty of the term $\gamma (1/{\omega }_n$) did not cause a significant change in the costs. The costs of ${\omega }_n, \left|G\left({\omega }_i\right)\right|$ and $\left|G^{\prime }\left({\omega }_i\right)\right|$ given in Table 2 were obtained using optimal parameters. Table 2 shows that the NP constraint always gave lower cost functions, $min\left(F_i\right)$, i = 1, 2, 3, than those obtained without the NP constraint.

In brief, using the NP constraint resulted in a lower total cost function compared to the results obtained without its use. Without using the NP constraint, the optimization algorithm established a new parameter update vector in each iteration, resulting in both transfer function gain and its derivative being minimized simultaneously. However, these were not independent developments; therefore, more iterations of the optimization algorithm were performed. Conversely, using the NP constraint ensured that a new parameter vector was defined such that both the transfer function gain and its derivative were minimized while the NP constraint was satisfied. Therefore, the parameter search took place in the physically realizable region. Here, the optimum solutions obtained with and without the use of the NP constraints in the frequency domain differed. The optimization problem using the NP equations found the global minimum in a smaller number of iterations.

3.2. Optimization of a Truss System

As a second example, a randomly loaded truss structure was considered, shown in Figure 5, and a frequency domain fatigue damage index was included in the cost function. The principal objective was to test the effect of the NP constraint on multi-objective optimization problems with stress minimization and multiple load cases.

Using the truss element formulations in [51], the corresponding mass matrix M, stiffness matrix K, and damping matrix C were calculated using finite element method (FEM) equations. Then, the FRFs were used to calculate the Pick matrix in Equation (1). The cost function was formulated as follows:

$$\underset{h_n}{\mathit{Minimize}}f=\alpha f_1+\beta f_2+\gamma f_3+\lambda f_4$$

(11)

$$Subject to\left\{\begin{array}{c}-{Ʌ}_j+\epsilon \le 0\\ h_{lb}{\le h}_n\le h_{ub}\end{array} \right.$$

(12)

where the objective function terms are $f_1={\sum }_{i=1}^n\left|G\left({\omega }_i\right)\right|, f_2={\sum }_{i=1}^n\left|G^{\prime }\left({\omega }_i\right)\right|, f_3={\sum }_{i=1}^m{h}_i$, and $f_4={\overline{D}}^{DK}$. The FRFs $G\left({\omega }_i\right)$ are calculated using Equations (3)–(5). The notation ${\overline{D}}^{DK}$ refers to the Dirlik damage index [52]. Discrete frequencies ${\omega }_i$ are chosen as natural frequencies. Minimizing $f_4$ in Equation (11) corresponds to maximizing the fatigue life in the frequency domain, and $\alpha$, $\beta$, $\gamma$, and $\lambda$ are the weightings, defined as $0\le \alpha \le 1$, $0\le \beta \le 1,0\le \gamma \le 1$, $0\le \lambda \le 1,$ and $\alpha +\beta +\gamma +\lambda =1$. These are changed by an increment of 0.1 in the search for Pareto optimal sets. In addition, $h_n$ denotes the nth design variable, and $h_{lb}$ and $h_{ub}$ are its lower and upper bounds, respectively. The design parameter $h_n$ is the weight of element number n and $f_3$ is the total cross-sectional area of the truss members, which is proportional to the total weight. Note that $\gamma =0$ corresponds to no penalty on the weight; therefore, this was excluded in the Pareto optimal solutions.

The random excitation force was arbitrarily chosen to be $f\left(t\right)={10}^3\left(sin\left(3t\right)+0.2sin\left(7t\right)+0.1sin\left(13t\right)+0.05ran{d}_{uniform}\left(\mathrm{0,1}\right)+0.02ran{d}_{Gaussian}\left(\mathrm{0,1}\right)\right)\left[N\right]$. The system FRFs were calculated in response to this random transient force.

The lower and upper limits for the cross-sectional areas were $h_{lb}=10$ and $h_{ub}=100 \left[m{m}^2\right]$. To calculate the fatigue damage in the frequency domain, the Dirlik damage index ${\overline{D}}^{DK}$ [52] was utilized to maximize the fatigue life of the structure and the weight was minimized under other constraints. In addition, $ran{d}_{uniform}\left(\mathrm{0,1}\right)$ and $ran{d}_{Gaussian}\left(\mathrm{0,1}\right)$ were uniform and Gaussian random numbers were generated in the interval [0, 1], respectively. The material parameters for the S-N curve were given by a fatigue strength constant of A = 1298.5 MPa, inverse slope of b = −0.125, k = 8, and an ultimate tensile strength of UTS = 510 MPa [53].

Using the cost function of Equation (11) and the constraints of Equation (12) with the interior point method, Pareto optimal solution sets were computed with and without the use of the NP constraints, and are shown in Figure 6. It is observed that the Dirlik fatigue damage index ${\overline{D}}^{DK }$ reduces as the optimum weight increases; thus, there is a trade-off between the fatigue damage index ${\overline{D}}^{DK }$ and the weight. In addition, the frequency domain response gain $\left|G\left(\omega \right)\right|$ falls as the fatigue damage index ${\overline{D}}^{DK }$ is reduced, which indicates that heavier structures result in less deformation and a smaller fatigue damage index. The result of the total cost function Equation (11) was 2.6144 e⁻⁴ with the NP constraint, while it was 0.0270 without the NP constraint. Overall, the use of NP constraints provided smaller cost functions for multi-objective optimization runs. If the number of runs were increased to obtain Pareto optimal sets, the minimum cost function values, both with and without NP-constrained solutions, grew closer to each other. Conversely, when one or more of the weightings (e.g., $\alpha$, $\beta$, $\gamma$ or $\lambda$) were chosen to be zero, the problem formulation excluded some cost functions that could lead to misleading optimum solutions. In parallel with the observations made in the previous section, imposing the NP constraint resulted in a smaller cost function.

Subsequently, the NP-constrained frequency domain optimization formulations were solved using stochastic optimization methods. The PSO algorithm was preferred for this purpose over the GA because it was reported in [54] that the PSO possessed computational superiority over the GA. When solving the NP-constrained optimization problems, the constrained PSO algorithm in [55] was adapted to our problem. An unconstrained PSO algorithm in Matlab was used when the optimization problem was solved without the NP constraint. The numerical results of the NP-constrained optimization problem solved by the constrained PSO algorithm are shown in Figure 7. The total cost function Equation (11) result found by the PSO was 5.7149 e⁻⁴, which was higher than that obtained using the interior point method. On the other hand, the CPU time of the constrained PSO algorithm was about four times longer than that of the interior point method when ran for this example on an Intel i7 CPU computer.

3.3. Topology Optimization of a U-Beam

Considering the simplified 3D engine traverse model shown in Figure 8, the NP constraint was imposed on a 3D topology optimization problem using a solid isotropic material with a penalization (SIMP) approach [56]. The following multi-objective optimization problem given in Equations (11) and (12) was solved.

The objective functions were $f_1={\sum }_{i=1}^n\left|G\left({\omega }_i\right)\right|, f_2={\sum }_{i=1}^n\left|G^{\prime }\left({\omega }_i\right)\right|, f_3={\sum }_{e=1}^m{\rho }^e{V}^e$, and $f_4={\overline{D}}^{DK }$. Discrete frequencies ${\omega }_i$ were chosen as the natural frequencies. Note that it is common to use compliance minimization for the cost function in topology optimization problems in the literature. The compliance minimization was achieved equivalently by minimizing the term $f_1={\sum }_{i=1}^n\left|G\left({\omega }_i\right)\right|$ in Equation (11), since our formulations were based on the NP constraint, which is a function of the FRFs $G\left({\omega }_i\right)$. The total volume was calculated using ${\int }_{\Omega }^{}{\rho }^{e}d{\Omega }^e=V$, where the artificial density of the elements is denoted by ${\rho }^e$, bounded by $\underset{\_}{\rho }\le {\rho }^e\le \overline{\rho }, {\rho }^eϵ\Omega$, $\underset{\_}{\rho }=0.01, \mathrm{a}\mathrm{n}\mathrm{d} \overline{\rho }=1$. Note that there were attempts to solve Equation (11) by imposing the volume fraction ratio in addition to the other constraints; however, these efforts caused convergence difficulties. Therefore, the volume fraction ratio was not included in the solutions. Imposing the volume fraction ratio in a stable way will be examined in the future studies.

The dimensions of the block material of the side walls in the x, y, and z directions were 700 mm, 150 mm, and 200 mm, respectively. Loads were applied to the elements on their top faces, as shown in Figure 8, while the elements at the bottom corners were fixed, i.e., there was an unsymmetric loading condition. The material was assumed to be steel with a density of ${\rho }_{steel}=7.7e^{-9} tonnes/{mm}^3$, Young’s modulus of 210,000 MPa, and Poisson’s ratio of 0.29. The S-N curve parameters of the material were the same as the material parameters given in Section 3.2. The applied random force was similar to that in Section 3.2: $f=f\left(t\right)=F_0\left(sin\left(3t\right)+0.2sin\left(7t\right)+0.1sin\left(13t\right)+0.05ran{d}_{uniform}\left(\mathrm{0,1}\right)+0.02ran{d}_{Gaussian}\left(\mathrm{0,1}\right)\right)\left[N\right]$. Its magnitude was $F_0=\frac{5000}{8}N$. The element formulations for 3D elastic solids in [57] were used, and the FRFs $G\left({\omega }_i\right)$ in Equation (11) were calculated according to Equations (3)–(5). The half-symmetry of the model in the y-z plane was used and there were 160 elements in total. The first attempts to solve Equations (11) and (12) were conducted using the conventional gradient-based optimization algorithm of the interior point method; however, CPU times were too high and convergence difficulties were observed. Therefore, the PSO algorithm was used, which constitutes the basis for the numerical solutions presented here. The weightings $\alpha$, $\beta$, $\gamma$, and $\lambda$ in Equation (11) were used to scale the individual costs $f_i$, so that their nominal values were approximately equal to each other.

When the optimization problem was solved with the NP constraint, the constrained PSO algorithm [58] was used. Similar to Section 3.2, when the optimization problem was solved without the NP constraint, the unconstrained PSO algorithm in Matlab was used. Figure 9 shows the optimization solutions obtained both with and without the NP constraint. The minimum total cost function with the NP constraint was 0.0228, obtained at iteration number 1965 with a total artificial density of 84.161, while 0.0229 was obtained at iteration number 3887 with a total artificial density of 96.468. In short, the global minimum was reached in fewer iterations for the NP-constrained case, and its optimum solution indicated thinner elements. Furthermore, it was observed that the solutions obtained without the NP constraint showed greater scatter during the iterations, which is reflected in the cost function plots in Figure 9. The reason for this is as follows: without using the NP constraint, the optimization algorithm attempts to reduce both the transfer function gain and its derivative simultaneously at each iteration; however, they are not independent. Such an attempt results in greater scatter in the cost functions during the iterations. On the other hand, by using the NP constraint we implicitly consider the relationship between the transfer function gain and its derivative; thus, the parameter vector search takes place in the physically realizable region and the global minimum can be achieved in a smaller number of iterations. The distributions of the artificial densities are shown in Figure 10 for both cases, with and without the NP constraint. Note that the solutions obtained with the NP constraint show a smoother density distribution with a lighter design and exhibit lower-density regions with slopes of approximately 45 degrees. Additionally, the optimum topology solutions without the NP constraint did not converge. This phenomenon is illustrated by the sample artificial density distribution in Figure 10b. In fact, Figure 9 shows the scattering behavior in cost terms in Equation (11) without the NP constraint. In Figure 9, all cost terms obtained without the NP constraint show large scattering across iterations and the results did not converge at all.

Similar to the case reported in [59], unsymmetric topological solutions were observed in our optimum topologies in Figure 10a. In [59], the authors claimed that nonlinear and non-convex constraints were the contributors to unsymmetric optimum solutions for the topology optimization problems in their study. We believe that the same reasons also caused our unsymmetric topology solutions when the NP constraint was used. This feature requires further investigation in future studies.

To reveal the reason for the oscillatory behavior of the cost function when solving Equations (11) and (12) without the NP constraint, the following two cases were also solved: first, $\beta =0$ was selected in Equation (11), the results of which are shown in Figure 11a; then, both $\alpha =0$ and $\beta =0$ were selected in Equation (11), the results of which are shown in Figure 11b. In short, only when $\alpha =0$ and $\beta =0$ were selected was convergence achieved without the NP constraint for solving Equation (11); however, the corresponding optimization formulation was the solution to a different problem than the one with non-zero $\alpha$ and $\beta$ values.

4. Conclusions

A new formulation for frequency domain optimization problems was investigated via an analysis of NP interpolation theory. We emphasized the minimization of frequency response magnitudes and their slopes. In addition, additional costs such as the minimization of the fatigue damage index and the maximization of natural frequencies were also considered. Our results showed that the physically realizable solution space in an optimization problem could be determined by considering the NP constraint equations. Solutions to the numerical examples were presented to prove the mathematical derivations. We discovered that resonance frequencies should be handled with special attention on account of the positive definiteness of the Pick matrix changes. Solutions of sample multi-objective optimization problems implied that considering the NP constraint resulted in relatively smaller optimum cost functions and smoother convergence in comparison to conventional optimization formulations. It was found that the stress level minimization achieved in frequency domain fatigue damage minimization caused convergence difficulties, which were overcome by the use of NP constraint equations. For example, significant convergence difficulties were observed in the topology optimization problems, especially when no NP constraint was used.

In the future, the derived formulation will be implemented to resolve transient time response problems. Moreover, the convergence properties of parameters and topology optimization problems using the NP constraint must be further investigated. Convergence difficulties were encountered when the volume fraction ratio constraint was imposed in topology optimization problems. The reason behind this difficulty will be studied in the future. Furthermore, the NP constraint will be utilized in the simultaneous optimization of structures and their controllers to investigate the achievable bounds of controller performance.