{'fraction_non_alphanumeric': 0.03855185074499424, 'fraction_numerical': 0.02241386671220595, 'mean_word_length': 4.493433395872421, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Structural optimization of the Ziegler's pendulum: singularities and exact optimal solutions

Structural optimization of non-conservative systems with respect to stability criteria is a research area with important applications in fluid-structure interactions, friction-induced instabilities, and civil engineering. In contrast to optimization of conservative systems where rigorously proven optimal solutions in buckling problems have been found, for non-conservative optimization problems only numerically optimized designs were reported. The proof of optimality in the non-conservative optimization problems is a mathematical challenge related to multiple eigenvalues, singularities on the stability domain, and non-convexity of the merit functional. We present a study of the optimal mass distribution in a classical Ziegler's pendulum where local and global extrema can be found explicitly. In particular, for the undamped case, the two maxima of the critical flutter load correspond to a vanishing mass either in a joint or at the free end of the pendulum; in the minimum, the ratio of the masses is equal to the ratio of the stiffness coefficients. The role of the singularities on the stability boundary in the optimization is highlighted and extension to the damped case as well as to the case of higher degrees of freedom is discussed.


Introduction
Structural optimization of conservative and non-conservative systems with respect to stability criteria is a rapidly growing research area with important applications in industry [1,2,3].
Optimization of conservative elastic systems such as a problem of optimal shape of a column against buckling is already non-trivial because some optimal solutions could be multi-modal and thus correspond to a multiple semisimple eigenvalue which creates a conical singularity of the merit functional [3]. Despite these complications, a number of rigorous optimal solutions are known in conservative structural optimization. Nevertheless an increase in the critical divergence load given by the optimal design in such problems is usually not very large in comparison with the initial design [2,3].
In contrast to conservative systems, the non-conservative ones can loose stability both by divergence and by flutter. It is known that mass and stiffness modification can increase the critical flutter load by hundreds percent, which is an order of magnitude higher than typical gains achieved in optimization of conservative systems [4]- [18]. For example, Ringertz [9] reported an 838% increase of the critical flutter load for the Beck's column [19] from 20.05 for a uniform design to 188.1 for an optimized shape. Recently Temis and Fedorov [17] found for a free-free beam moving under the follower thrust an optimal design with the critical flutter load that exceeds that for a uniform beam by 823%. We note that despite the very notion of the follower forces was debated in [20,21,22], the Beck's column [19] as well as its discrete analogues [23,24,25] including the Ziegler's pendulum [26] remain popular models for investigation of mode-coupling instabilities in non-conservative systems and related optimization problems.
In both conservative and non-conservative problems of structural optimization of slender structures, their optimal or optimized shapes often possess places with small or even vanishing cross-sections. The known optimized shapes of the Beck's column or of a free-free rod moving under follower thrust have almost vanishing cross-section, e.g., at the free end, which means vanishing mass of a finite element in the corresponding discretization [9,10,11,12,17,18].
Another intrigue of optimization of non-conservative systems is the 'wandering' critical frequency at the optimal critical load. During the optimization the eigenvalue branches experience numerous mutual overlappings and veerings [4,6,7,27,28,29] with the tendency for the critical frequency to increase and to correspond to higher modes [9,14,15,16,17,30]. This puzzling behavior of the critical frequency still waits for its explanation.
In some problems, such as the optimal placement of the point mass along a uniform free-free rod moving under the follower thrust [5,8], the local maxima were found to correspond to singularities of the flutter boundary such as cuspidal points where the multiple eigenvalues with the Jordan block exist [11,12]. In order for the last phenomenon to happen on needs at least three modes [11,12] meaning that the two-modes approximations [5,8] are unable to detect such optima. This reflects a general question on model reduction and validity of low-dimensional approximations in non-conservative problems discussed already by Bolotin [30] and Gasparini et al. [23] and recently raised up again in the context of friction-induced vibrations by Butlin and Woodhouse [31].
The above mentioned phenomena make rigorous proofs of the optimality in the non-conservative optimization problems substantially more difficult than in conservative ones. To the best of our knowledge, the rigorously proven optimal solutions in optimization problems for distributed circulatory systems were not found. Although in the finite-dimensional case the situation is not much better, it looks reasonable to try to understand the nature of the observed difficulties of optimization on the final dimensional non-conservative systems that depend on a finite number of control parameters.
Let us consider a circulatory system where dot indicates time differentiation, M is a real symmetric m × m mass matrix and K is a real non-symmetric m × m matrix of positional forces that include both potential and non-potential (circulatory) ones. The equation (1) typically originates after linearization and discretization of stability problems for structures under follower loads, in the problems of friction-induced vibrations and even in rotor dynamics when the damping is not taken into account [30,28,32].
The characteristic equation for the circulatory system (1) is given by the modified Leverrier's algorithm [33]. In case when m = 2, it reads det Mλ 4 + (trMtrK − tr(MK))λ 2 + det K = 0, where λ is an eigenvalue that determines stability of the trivial solution.
Usually, the coefficients of the matrix K contain the loading parameter, say p, that one needs to increase by varying the coefficients of, e.g., the mass matrix. During this optimization process some masses can come close to zero so that the mass matrix can degenerate yielding det M = 0. As a consequence, some eigenvalues λ can substantially increase. On the other hand such a singular perturbation of the characteristic polynomial may cause large values of the gradient of the critical load with respect to the mass distribution. We see that a problem of optimal mass distribution in a finitedimensional circulatory system (1) in order to increase the critical flutter load, looks promising for explanation of the peculiarities of optimization of distributed non-conservative structures. However, it makes sense to tackle first not the most general system (1) but rather a particular finite-dimensional non-conservative system with m degrees of freedom.
In this article, we propose to take an m-link Ziegler's pendulum [23,24,25,26] as a toy model for investigation of optimal mass and stiffness distributions that give an extremum to the critical flutter load. It appears that even the classical two-link Ziegler's pendulum had rarely been studied from this point of view in contrast to its continuous analogue -the Beck's column. The only example of such a study known to the author is contained in the book by Gajewski and Zyczkowski [1].
The paper is organized as follows. In the next Section we first consider optimization of an undamped two-link Ziegler's pendulum. We find an explicit that correspond to the points where the flutter boundary has a vertical tangent such as cuspidal points with triple eigenvalues, cf. [11,12], or the points where the flutter boundary experiences a sharp turn, cf. [27,28,29]. Finally we consider the problem of optimal mass distribution for a two-link Ziegler pendulum with dissipation. In conclusion, we formulate an optimization problem for an m-link Ziegler's pendulum and discuss some hypotheses on plausible optimal solutions and their expected properties.

Structural optimization of the Ziegler's pendulum
Let us consider the classical Ziegler's pendulum consisting of two light and rigid rods of equal length l. The pendulum is attached to a firm basement by a viscoelastic revolute joint with the stiffness coefficient c 1 and the damping coefficient d 1 . Another viscoelastic revolute joint with the stiffness coefficient c 2 and the damping coefficient d 2 connects the two rods [26,32].
At the second revolute joint and at the free end of the second rod the point masses m 1 and m 2 are located, respectively. The second rod is subjected to a tangential follower load P [26, 32].

Undamped case
Small deviations from the vertical equilibrium for the undamped Ziegler's pendulum are described by the equation (1) with the mass and stiffness matrices that have the following form [26,32] where x = (θ 1 , θ 2 ) T is the vector consisting of small angle deviations from the vertical equilibrium position.
Calculating the characteristic equation det(Mλ 2 +K) = 0 for the Ziegler's pendulum without dissipation, we find By direct calculation of the roots of the characteristic equation (4) or by using the Gallina criterion [24] we find a critical surface that separates flutter instability and marginal stability domains The equation (5)  However, in the (m 1 , m 2 , P )-space the critical surface (5) looks differently and has a form of a self-intersecting surface with the Whitney umbrella singularity at the (0, 0, 2c 2 /l)-point. Indeed, expressing the critical load P from Defining the new critical flutter load as p := P l/c 2 we come to the more symmetric expression The case when c 1 = c 2 = 1, m 1 = 2 and m 2 = 1 corresponds to the classical result of Ziegler [26] The lower value of the critical load correspond to the boundary between marginal stability and flutter while the higher critical load corresponds to the conventional transition from flutter to divergence through the double zero eigenvalue with the Jordan block.
The critical load (7) as a function of the masses p = p(m 1 , m 2 ) is plotted in Fig. 1(b). It is seen that the stability boundary has a self-intersection along a ray of the p-axis that starts at the a Whitney umbrella singularity with the coordinates (0, 0, 2) in the (m 1 , m 2 , p)-space. Due to the symmetry of the expression (7), the critical load as a function of the stiffness coefficients, p = p(c 1 , c 2 ), forms the identical surface in the (c 1 , c 2 , p)-space. In the following we will study the function p = p(m 1 , m 2 ).

Divergence: both modes
One mode stable, another unstable by divergence One mode stable, another unstable by divergence According to the inequality (7) the critical load is always not less than The minimum is reached when the masses satisfy the constraint Note that the equal stiffness coefficients c 1 = c 2 imply equal masses m 1 = m 2 .
This situation corresponds to a uniformly distributed mass and stiffness in continuous systems such as the Beck's column [19].
Usually, in the structural optimization problems the uniformly distributed stiffness and mass are considered as the initial design that is a starting point in optimization procedures. The critical load of the optimized structure is conventionally compared to that of the same structure with the uniform mass and stiffness distributions [4,5,6,7,8,9,10,13,14,15,16,17,18].
Since p(m 1 , m 2 ) is a ruled surface and thus p effectively depends on the mass ratio only, it is convenient to introduce the azimuth angle α by assuming m 1 = cos α and m 2 = sin α and to plot the critical load as a function of α.
In Fig. 2 the curves p = p(α) bound the flutter domain shown in the light gray. When α tends to zero, which corresponds to the vanishing mass m 2 , the critical load increases to infinity. When α tends to π 2 and, correspondingly, the mass m 1 is vanishing, then the critical flutter load increases to the value At the point B in the stability diagram of Fig. 2 the flutter boundary has a vertical tangent, which is a typical phenomenon in non-conservative optimization [27,28,29]. higher modes were reported to be involved into the coupling that indicates the onset of flutter [9,13,14,15,16,17,18]. Our simple model shows that this phenomenon seems to be natural for the optimal design that causes the degeneracy in the mass matrix that gives rise to the critical frequency that increases without bounds. It would be interesting to derive the optimal solutions for the optimization of the Ziegler's pendulum by means of the Pontryagin's maximum principle, which is used in numerical optimization of the distributed non-conservative systems, see e.g. [10,11].

The m-link Ziegler's pendulum
It would be very much desirable to extend our study to the case of the multiple-degrees-of-freedom Ziegler's pendulum. The corresponding models and recursive schemes for deriving the equations of motions where proposed by Gasparini et al. [23] and by Lobas [25]. A three-link Ziegler's pendulum was considered by Gallina [24].
We take the linearized equations of Lobas [25] for its model deals with the arbitrary masses and stiffnesses in the joints of an m-link Ziegler's pendulum in contrast to the models of Gasparini and Gallina [23,24]. The mass and stiffness matrices of the Ziegler-Lobas model look like −c 2 c 2 + c 3 − P l · · · 0 P l · · · · · · · · · · · · · · · 0 0 · · · c m−1 + c m − P l −c m−1 + P l For m = 2 the matrices (11) and (12) are reduced to (3). Note that  For simplicity we assume that c 1 = c 2 = c 3 = c. Then, the characteristic polynomial has the form a 0 λ 6 + a 1 λ 4 + a 2 λ 2 + a 3 = 0, with the coefficients a 0 = l 6 m 1 m 2 m 3 , Composing the discriminant matrix [24] and calculating the discriminant sequence consisting of the determinants of the three main minors of even order, we find that ∆ 1 = 3l 12 m 2 1 m 2 2 m 2 3 > 0. The expressions for determinants ∆ 2 and ∆ 3 = det S are rather involved and by this reason we omit them here. However, the numerical experiments evidence that the stability boundary is given by the equation ∆ 3 = 0 for the stability condition ∆ 3 > 0 implies ∆ 2 > 0.
In Fig. 3 using the inequality ∆ 3 > 0 we present the stability diagrams in the (α, P )-plane where the azimuth angle α in the (m 2 , m 3 )-plane is introduced by assuming m 2 = r cos α and m 3 = r sin α. We assume the equal lengths of the links l = 1 and the equal stiffness coefficients c 1 = c 2 = c 3 = 1 and vary the radial distance r in the (m 2 , m 3 )-plane and the mass m 1 .
Since for m = 3 the critical surface P (m 2 , m 3 ) is no more a ruled surface as it was in the case m = 2, the pictures in the (α, P )-plane change with the variation of the radial distance r that complicates the optimization problem.
Nevertheless such diagrams are convenient for the analysis of the geometry of the stability boundary and thus for the identification of potential extrema.
Moreover, the critical surface P (m 2 , m 3 ) has a self-intersection along a ray of the P -axis that starts from the singularity Whitney umbrella as in the case of the 2-link pendulum. Therefore, at small values of m 2 and m 3 the critical load can be locally approximated by a ruled surface.
In the left column of Fig. 3 the radial distance r in the (m 2 , m 3 )-plane is fixed to r = 1 while the mass m 1 is increasing. As in the case m = 2, (marginal) stability is possible for α ∈ [0, π/2]. For m 1 = 0, two finite maximal values P A and P B are identified at α = 0 (m 3 = 0) and α = π/2 (m 2 = 0), respectively, Fig. 3(a). Both maxima are attained at the cuspidal points of the stability boundary where it has the vertical tangents. However, the stability diagram changes when m 1 = 10, Fig. 3(b). Again, local extrema exist at the boundary points α = 0 (m 3 = 0) and α = π/2 (m 2 = 0), while the global maximum is at the point of the sharp turn of the flutter boundary with the vertical tangent near the cuspidal point C 1 ≃ (0.0403477, 11.961144) that corresponds to triple pure imaginary eigenvalues λ ≃ ±i1.1635243 with the Jordan block of third order, cf. [11,12] where at such a singular point a maximum of the critical flutter load was found for a free-free beam under the follower thrust.
With the further increase in the first mass up to m 1 = 200, the stability diagram converges to that similar to the diagram of the two-link pendulum, cf. Fig. 2 and Fig. 3(c). This is not surprising because big inertia of the first joint makes the three-link pendulum effectively a two-link one.
The right column in Fig. 3 corresponds to the fixed first mass m 1 = 5 and varied radial distance r in the (m 2 , m 3 )-plane. Small values of r correspond effectively to a two-link pendulum. That is why Fig. 3(f) with r = 0.1 looks similar to Fig. 3(c) and Fig. 2.
Increase of r is accompanied by a complication of the stability diagram.
In particular, two cusp point singularities originate corresponding to triple pure imaginary eigenvalues with the Jordan block, Fig. 3(d,e). Near these singularities, the stability boundary experience a sharp turn at the points B 1 and B 2 where the tangent to the boundary is vertical. At such points the eigencurves Imλ(P ) form crossing that can change either to avoided crossing or to the overlapping with the origination of a bubble of complex eigenvalues in dependence on the direction of variation of the azimuth angle α. This phenomenon was observed in many numerical studies of non-conservative optimization problems [4,6,7,9,10,13,14,17,18] and had been described analytically by Kirillov and Seyranian in [27,28]. Moreover, the critical load at these points can jump to a higher value corresponding to the merging of other (often higher) modes and is thus undefined [4,6,7]. Nevertheless, these points could be local extrema of the merit functional, see [29] where the necessary conditions for that were derived.
We stress that due to finite number of degrees of freedom and finite number of control parameters the stability boundary of an m-link Ziegler's pendulum can be thoroughly analyzed both analytically and numerically.
In particular, the coordinates of the singular points can be calculated with arbitrary precision and thus the issues of high sensitivity of the critical flutter load at the optima could be successfully resolved in this model, contrary to the optimization problems for distributed systems. A possibility to work with the singularities related to coalescence of more than two eigenvalues allows one to investigate qualitatively the question 'Should low-order models be believed?' [31]. Indeed, two-mode approximations work well far from such singularities, in their vicinity one has to take into account higher modes.

Damped case
In the presence of dissipation the equation of the Ziegler's pendulum is with the damping matrix [25,26] For the two-link damped Ziegler's pendulum with m = 2, we find the characteristic equation Applying the Routh-Hurwitz criterion, we find the critical flutter load For c 1 = c 2 = 1, l = 1, m 1 = 2 and m 2 = 1 it was found to be [34] The equation (22) defines a surface with the Whitney umbrella singularity in the (d 1 , d 2 , P )-space which explains Ziegler's destabilization paradox by vanishing dissipation [26] as it was first demonstrated by Bottema already in 1956 [32,35].
In contrast to the previous studies, e.g. [32,34,35,36], we consider the critical flutter load (21) as function of masses P = P (m 1 , m 2 ) for the fixed damping distribution.
In Fig. 4, the stability diagrams for the damped two-link Ziegler's pendulum are shown in the assumption of c 1 = c 2 = 1, l = 1, m 1 = r cos α and m 2 = r sin α. The critical load P (m 1 , m 2 ) does not constitute a ruled surface and thus the function P (α) depends on the radial distance r in the (m 1 , m 2 )-plane. We see that the extrema again correspond to the boundary points m 1 = 0 and m 2 = 0 although the singularity at α = π/2 is an intersection point. Nevertheless, with the increase of the number of degrees of freedom and control parameters new types of singularities will appear. The planar stability diagrams of Fig. 4 depend on the damping distribution but do not tend to that of the undamped pendulum when damping goes to zero (destabilization paradox [32,35,36]).

'Problema novum ad cuius solutionem Mathematici invitantur'
'A new problem that mathematicians are invited to solve' is a translated from Latin title of 1696 work by J. Bernoulli where he proposed the famous Brachistochrone problem [37]. Supporting this good old tradition we would like to formulate the following optimization problem: Given a circulatory system (1) with the matrices M and K defined as in (11) and (12). Find all local extrema, the absolute maximum of the critical flutter load P , and the corresponding extremal mass distributions {m 1 , m 2 , . . . , m m }.
One can consider also the problem of optimal stiffness distribution or even vary both stiffness and mass. The same problems can be formulated for the damped pendulum with the damping matrix (18).
We expect that both in the undamped and in the damped case there exists a class of extrema, corresponding to the distributions {m 1 , m 2 , . . . , m m } with some masses m i = 0. It should be possible to find these optimal mass distributions explicitly perhaps with the use of the Pontryagin's maximum principle. It would be interesting to identify the singularities of the stability boundary that correspond to these extrema; some of them should be related to infinite eigenvalues λ.
On the other hand, some local extrema should exist with the mass distributions that do not contain the vanishing masses m i . It would be interesting to understand at which points -smooth or non-smooth -of the stability boundaries they are attained. Since the system is finite-dimensional and contains the finite number of control parameters with the clear physical meaning, the locations of the singularities corresponding to multiple eigenvalues can easily be found numerically with the high accuracy. In the vicinity of such points where at least three pure imaginary eigenvalues couple, the question 'Should low-order models be believed' [31] makes sense because here a one more degree of freedom is crucial for the correct solution.
Knowledge of the rigorously established optimal solutions at every m should shed light on the behavior of the optimal mass distributions and critical frequencies with the increase of the number of degrees of freedom. As a by-product, such a study will give an insight to the problem of dimension reduction and will serve a nice playground both for application of the existing methods of nonsmooth analysis and eigenvalue optimization [38,39,40] and for their further development in the direction of more tight relation both with the singularity theory and the applications' needs. We expect that the proposed model optimization problem will yield useful recommendations for improvement of the algorithms of optimization of real non-conservative structures in industry.

Conclusions
In the present article we proposed to consider an m-link Ziegler