#### Abstract

The optimal design of damping parameters for passive vibration control remains a challenge for both research and industrial applications. Here, we introduce a design methodology based on limit cycle analysis in concert with design optimization. A state-space representation is used to model the vibrational behavior converged to its limit cycle. The design approach is outlined and applied to mechanical systems undergoing periodic forces. This method is applicable to both vibration mitigation and energy harvesting, and examples of both are shown. We conclude with a summary of the results and an outlook for future developments and applications.

#### 1. Introduction

Proper design for the control of vibrations in mechanical systems is a problem found in both fields of mechanical [1, 2] and civil engineering [3]. Here, we address the passive control of vibrations in which so-called passive elements, i.e., masses, springs, and dampers, are optimally designed. Systems are assessed considering only their converged limit cycle while neglecting the transient part of the response. Specifically, design optimization problems of mechanical systems related to energy harvesting and vibration mitigation are treated. For both types of problems, energy dissipation is the driving physical phenomenon. A general design optimization formulation is provided in which the method for calculation is based on the state-space representation and limit cycles of dynamic oscillators.

In the literature, several design and design optimization formulations are proposed for this topic. The formulations can be generally categorized in the domain in which the problems are solved and further decomposed into subcategories of the specific solving approach. This proposed categorization is as follows:(1)Time domain:(a)Formulation as a first-order differential equation via state-space representation (e.g., [4, 5])(b)Time integration of the second-order differential equations, which retain the structure of Euler–Lagrange equations (e.g., [6])(2)Frequency domain:(a)Modal analysis (e.g., [7–9])(b)Frequency response and transfer functions (e.g., [10])(c)Inverse eigenvalue problems (e.g., [11–13])

Optimal design problems of structures under dynamic aspects are handled in the literature of structural design optimization, and definitions and formulations are handled generally in [14–17]. These problems are typically solved using a frequency-domain formulation of the vibrational problem [18–20]. Natural frequencies are typically treated in these cases as lower-bounded constraints (e.g. [7]) though a maximization of the lowest eigenfrequency may also be carried out (e.g., [9]). Formulations of constrained eigenfrequency problems can be found in [21], which have been further developed in [8, 22].

The design problem can also be defined that a certain eigenfrequency is desired, and this is then formulated as an inverse eigenvalue problem. This method is typically referred to as *dynamic structural modification* to differentiate it from structural optimization, and examples of this can be found in [13, 23–25].

In this paper, the first category of the list above is of interest and specifically within the subcategory of state-space approaches: modeling with the first-order differential equation via state-space representation and solving in time domain. A design optimization approach in concert with state-space representation to model linear mechanical systems is proposed. Using asymptotic analysis, we are able to achieve the benefits of both state-space and frequency-domain methods. Frequency-domain methods address the steady state of the system, neglecting the transient state. Objective and constraint functions can, though, only be formulated in the frequency domain [19, 20]. On the contrary, time-domain methods can address arbitrary cost functions [4], but they do not address naturally the steady state of the system. Here, the optimization formulation is based on the objective function P, which is an infinite-horizon passive control problem, i.e., minimization of the integral of the following form:where *L* is an arbitrary function called the Lagrangian, is the state variable vector, and *t* is the time.

The model is based on asymptotic convergence under certain assumptions of a linear viscoelastic mechanical system under a periodic force having a fundamental period , which leads to a limit cycle with the same fundamental period. This allows us to formulate arbitrary objective and constraint functions in the steady state of the system. The limit cycle is expressed explicitly in terms of the design variables and the time *t*, allowing us to explicitly formulate the objective function (1) in terms of the design variables.

#### 2. Mathematical Model for Vibrating Mechanical Systems

The governing equation for a linear mechanical system is the following second-order differential equation:where is the () mass matrix, is the () damping matrix, is the () stiffness matrix, is external force vector, is the position vector of the system in generalized coordinates, is the velocity vector, is the acceleration vector, and *n* is the number of degrees of freedom. The nomenclature showing dependence of refers here to the dependency of the system matrices on the design variable vector . It should be noted that, in this work, underlined symbols are vectors, double underlined symbols are matrices, and the symbols without an underline are scalars. It should be further noted that the optimization parameters (of general and nonmechanical nature) use Sans Serif symbols.

For a state-space description of the system equation (2), we introduce the state vector defined as

The system equation (2) is modeled in the state space aswhere the matrix is defined bythe vector isand is the () identity matrix. Using the state-space representation, the system equation (4) has the solution:where are the initial conditions and is the state-transition matrix associated with the homogeneous solution of the linear equation (equation (4)). The latter is defined by

Traditional analysis methods derive the steady-state solution from the linear superposition of the normal modes. Although the solution method introduced here is more complex, it is a general solution of equation (2), i.e., no assumptions must be made. Furthermore, this formulation has the following important advantages over frequency-domain analysis (i.e., modal analysis):

(1)This method allows any arbitrary damping that results in a (symmetric) positive-definite damping matrix; i.e., damping is not restricted to proportional or small damping.(2)It allows the derivation of an explicit solution to system equation (2) for its limit cycle, i.e., the steady state.The former advantage also implies that damping is necessary; i.e., undamped systems cannot be analyzed with this method. Furthermore, this method has the limitation of not being able to properly represent instability problems such as flutter, slip-stick friction, and others related to structure-fluid interactions.

In the following section, the solution to the steady-state behavior of equation (7) is shown in the time domain and without the use of Fourier or Laplace transforms.

#### 3. Limit Cycles for Vibrating Mechanical Systems under Periodic Excitation

The term *limit cycle* is traditionally addressed in the context of autonomous systems as a closed (periodic) trajectory such that its neighboring trajectories are not closed, i.e., isolated [26–28]. In this paper, the authors desire to recognize these kinds of objects in traditional problems of linear vibrations under periodic excitation with the motivation of analyzing the steady state of such systems in state space without making reference to complex-valued transfer functions. To address the problem from a geometric perspective, we will prove that the steady-state solution of equation (4) is closed and isolated. However, we underline that limit cycles do not exist in autonomous linear homogeneous systems. Regardless of the linearity of equation (4) with respect to the state variable, such a system is nonhomogeneous and nonautonomous.

We consider nonautonomous and nonhomogeneous mechanical systems, which are asymptotically stable; i.e., the matrix has eigenvalues with negative real parts. Asymptotic stability is a sufficient condition for the convergence of the system to a limit cycle. A sufficient condition for asymptotic stability is that the system matrices , , and are positive-definite [29]. Furthermore, if the symmetric part of the matrix is negative-definite, the state-transition matrix has the following properties:

The two properties of the state-transition matrix shown in equation (9) and equation (10) allow us to state the following theorem (proof found in [30]):

Theorem 1. *We consider a system equation (4) under a periodic force with a fundamental period . If the state-transition matrix of such a system has the property equations (9) and (10) for every initial condition, the system converges to a closed trajectory of period , described in [30]:where*

In turn, is the () identity matrix and

The solution shown in equation (11) is asymptotically approached from arbitrary initial conditions. Therefore, we confirm that the solution is periodic and isolated, which is the definition of a limit cycle, defined above. This theorem permits the approximation of the objective function equation (1) for cases in which the transitory response is negligible with respect to the steady-state contribution. To formulate the numerical design optimization problem, it is sufficient to introduce a time partition in the limit cycle. The objective function is approximated by a Gauss quadrature:where are the weighting functions.

The constraint function can be mapped to *N* constraints at each point of the time partition:

The derivation of the analytical expression equation (11) is as follows: such equations are defined on the interval , so we can write the periodic force as a Fourier series of an odd function:

We evaluate the integral in equation (11) such that

Substituting equations (12) and (17) in equation (11)) reveals

A generic quadratic objective function is described as

In this specific case, we havewhere

The design optimization problem utilizes both gradient-based and stochastic methods below. The general optimization formulation in this context can be formulated mathematically aswhere is the objective function, is the variable of design variables , is the constraint function (formulated here as an upper-bound constraint based on a response and constraint limit of this response , which will be specifically defined in the constrained design problem), is the lower limit, and is the upper limit of the *i*th design variable.

#### 4. Optimal Design of Damping for an Energy-Harvesting Beam

The first example is the optimal design of a bed of dampers, which maximize the energy harvested from a structure: here, a simply supported beam (cf. Figure 1). The beam structure has a constant cross-sectional geometry and material (i.e., Young’s modulus and density). Furthermore, it is pinned at both ends. The loading of the beam is a periodic point force located at (here, at the midpoint of beam). The bed of dampers consists of *m* dampers having damping coefficient located at , .

The structure is modeled as a Euler–Bernoulli beam of 20 elements with the system having a total of 60 degrees of freedom. The equilibrium of this system is stated as the following fourth-order differential equation:where *E* is Young’s modulus, *I* is the second-moment of area about the bending axis, *ρ* is the density, *A* is the cross-sectional area, is the vertical displacement of the neutral axis, is the internal damping of the beam’s material, is the length of the beam, and is the Dirac delta function, which is used to describe the point loading. The notation is simplified by reducing the number of parameters in our system with the following definitions of dimensionless values:

Equation (23) is solved using the Galerkin method, resulting in the approximation of nondimensional displacement as (shown here with the variables defined in equation (24)) follows:where () is a set of orthonormal functions that satisfy the boundary conditions of the problem. It should be noted that the nomenclature differentiates between the normalized general position *ξ*, the specific positions of the dampers , and the position of the force .

The discretized equations of motion in generalized coordinates arewhere ,

In this case, the matrices and are diagonal matrices and have the following components:

For *m* dampers, we define our design variable vector as

The bed of dampers absorbs the power described by the following:

Therefore, the objective function P, including its value in the limit cycle , is a dimensionless quantity proportional to the dissipated energy in a period *T*. This value is to be maximized in the case of energy harvesting.

The design variable vector is bounded in order to preserve the physical meaning of the position and to keep the optimization problem well conditioned. The damping coefficient is to remain positive, and therefore, we set

Furthermore, the positions of the dampers are constrained by

The constraint function of this design problem limits the displacement of the beam with an allowable normalized displacement :

We address this constraint at discrete points *j* where the displacement is assessed, which we denote as . The maximum displacement throughout the limit cycle is constrained:

The inequality constraint is then normalized giving

The value of is numerically evaluated at nine evenly spaced points. This design problem is carried out using an upper bound to the dimensionless displacementor one-tenth of the total length. We load the system with a dimensionless force:

We can now express the complete optimal design problem mathematically as

It should be noted that this is a maximization problem that can also be described as .

For the numerical design optimization, we utilize DesOptPy [31, 32] and specifically the interface to pyOpt [33] and the second-order algorithm NLPQLP [34, 35]. In this work, we do not provide the analytical sensitivities though this can be done in a future work to further reduce the computational effort. This optimization problem converges in 9 iterations and 100 evaluations. The optimization problem—including start design—and its results are shown in Tables 1 and 2.

In Figure 2, we can observe the state-space response of the first mode for both the initial and optimal designs. For this time-domain depiction, we have chosen an initial condition . It is important to note that the optimization above is valid for arbitrary initial conditions. From our chosen starting point, Figure 2(a) shows the convergence to the limit cycle. The limit cycle is shown alone (without the convergence path) in Figure 2(b).

**(a)**

**(b)**

To understand the difference in energy harvested from the initial and optimal designs, we show the time-domain response in Figure 3. It can be seen that the optimal design harvests less energy than the initial, which is due to the infeasible starting design; i.e., displacements 3 through 7 are violated (cf. Table 1). The first mode plot (Figure 2) is a good approximation of these initially violated displacements. The optimal design does not guarantee that the constraints are respected during the transitory portion of the motion of the system, as can be appreciated in Figure 2(a).

To verify that the use of gradient-based optimization is valid, a range of starting designs were investigated. All starting points result in the same optimal design, needing between 12 and 50 iterations (100 to 551 evaluations). This leads the authors to the assumption (though not mathematically verifiable) that the problem is convex or—at the very least—that gradient-based optimization is permissible for the efficient solving of this problem.

#### 5. Optimal Design of Dynamic Impact Absorbers Mitigating Vibration

In this example, we have a forced oscillator with one degree of freedom, which we augment with a series of dynamic impact absorbers, i.e., mass-damper-spring systems, to mitigate the vibration (Figure 4). The equation of motion of this system is the same as equation (2) with the following system matrices of size :

The design problem is to find the optimal values of these impact absorbers such that the energy is minimized.

The original system to which the dynamic impact absorbers are added will be referred to as the principal system and denoted with the subscript “0”. The principal system is described by the frequency ratio and the damping coefficient which are defined as follows:where is the natural frequency, is the damping of the principal system, and is the mass of the principal system. The dynamic impact absorbers are also described accordingly:

The definitions above result in normalized parameters as well reducing the number of parameters by one that are necessary to describe each subsystem. The kinetic energy of the vibrational energy is described by the Lagrangian:

The objective function (the value of P in the limit cycle) is again a dimensionless measure, which is proportional to the mean value of the energy in a period T. This value is to be minimized in the case of vibration mitigation. Mathematically, this is described by the following bounded, but unconstrained optimization problem:

This problem is inherently nonconvex and therefore a zeroth-order algorithm is chosen, namely, differential evolution as implemented in PyGMO [36] via DesOptPy. The algorithm options of a population size of 100 for 100 generations, i.e., 10000 system evaluations, were chosen. Due to the efficiency of the solving routine of the mechanical system, even this large number of evaluations takes less than two seconds on a laptop.

The state-space response of the principal system for an arbitrary solution and optimal design is shown in Figure 5. As we are using a stochastic optimization algorithm, we do not have an initial design as a reference and, therefore, use the upper-bound design. For this time-domain depiction, we have again chosen an initial condition as mentioned above. It can be observed from these two figures that the optimal design has decreased the displacement at the cost of a slight increase in velocity. In Figure 6, a drastic reduction in the vibrational energy can be seen with the optimal solution.

**(a)**

**(b)**

As this method is of stochastic nature, 100 complete optimization runs were carried out to test the robustness of the optimization algorithm in concert with the solution approach. In Tables 3 and 4, results show the system without dynamic impact absorbers, at the lower and upper bounds as well as an exemplary optimization. The histogram showing the results can be found in Figure 7. This shows satisfactory results due to its relatively “tight” spread though better results would be found by increasing the number of evaluations with the cost of more computational effort.

#### 6. Conclusion

This paper has presented a method for the optimal design of passive vibration control. This approach is based on the limit cycle of linear mechanical systems under a periodic force. The theory has been devolved and elaborated upon. Sufficient conditions for the existence of such limit cycle and both constrained and unconstrained formulations of the optimization problem have been shown and solved. The application of this design method to two categories of the vibrating system is introduced:

(1)energy harvesting in which high amounts of vibratory energy is desired and therefore maximized,(2)vibration mitigation where energy is to be minimized.The former example has shown to be well conditioned in its formulation of the optimization problem. Therefore, an efficient gradient-based optimization has been applied. The example of vibration mitigation is indeed inherently multimodal, and therefore, less efficient stochastic methods were investigated, resulting in the use of differential evolution. A general optimization framework has been provided; however, it is necessary to properly understand the characteristics of the optimization problem, including convexity.

Although the design optimization of the problems shown was successful, we have not considered the transient response. The presented optimal designs may experience suboptimal behavior and violate constraints while converging to the limit cycle. The authors see this as the main drawback of the introduced method as it cannot be applied to problems, which are strongly characterized by their transient response. This said there are a great number of problems that are indeed characterized solely by their converged limit cycle response, which can be handled with this methodology.

To conclude, we have shown a general procedure that has successfully been applied to typical benchmarks proving its feasibility. This methodology can also be applied to large-scale problems for the proper design of passive control of vibrations. Future work will consider larger industrial application problems; although as long as the system limitations introduced above are met, the authors do not see this as a hurdle to efficiency or effectivity to this methodology for finding optimal dynamic systems.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The research and results reported in this work are part of the project TN2803, “Mech4SME3: Mechatronics for Smart Maintenance and Energy Efficiency Enhancement,” funded by the Free University of Bozen-Bolzano. This work was supported by the Open Access Publishing Fund provided by the Free University of Bozen-Bolzano.