Singularities in Structural Optimization of the Ziegler Pendulum

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 nonconservative optimization problems only numerically optimized designs have been reported. The proof of optimality in non-conservative optimization problems is a mathematical challenge related to multiple eigenvalues, singularities in the stability domain, and non-convexity of the merit functional. We present here a study of optimal mass distribution in a classical Ziegler 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 an 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 the problem of the 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 semi-simple 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, the 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, non-conservative systems 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][5][6][7][8][9][10][11][12][13][14][15][16][17][18].For example, Ringertz [9] reported an 838 % increase of the critical flutter load for the Beck 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 a critical flutter load that exceeds the load for a uniform beam by 823 %.We note that although the very notion of the follower forces was debated in [20][21][22], the Beck column [19] as well as its discrete analogues [23][24][25] including the Ziegler pendulum [26], remain popular models for investigating 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 column or of a free-free rod moving under follower thrust have an almost vanishing cross-section, e.g., at the free end, which means vanishing mass of a finite element in the corresponding discretization [9-12, 17, 18].
Another intriguing feature of optimizing non-conservative systems is the 'wandering' critical frequency at the optimal critical load.During optimization the eigenvalue branches experience numerous mutual overlappings and veerings [4,6,7,[27][28][29] with a 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 awaits an 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 multiple eigenvalues with the Jordan block exist [11,12].In order for the last phenomenon to happen, we need at least three modes [11,12], meaning that two-mode approximations [5,8] are unable to detect such optima.This reflects a general question on model reduction and the validity of lowdimensional approximations in non-conservative problems, already discussed by Bolotin [30] and Gasparini et al. [23] and recently raised again in the context of friction-induced vibrations by Butlin and Woodhouse [31].
The above-mentioned phenomena make rigorous proofs of optimality in non-conservative optimization problems substantially more difficult than in conservative ones.To the best of our knowledge, no rigorously proven optimal solutions in optimization problems for distributed circulatory systems have been found.Although the situation is not much better in the finite-dimensional case, it seems reasonable to try to understand the nature of the observed difficulties of optimization on 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) forces.Equation (1) typically originates after linearization and discretization of stability problems for structures under follower loads, in 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 algorithm [33].In the case when m = 2, it reads where λ is an eigenvalue that determines the stability of the trivial solution.
The coefficients of matrix K usually contain the loading parameter, say p, that we need 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 increase substantially.However, 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 finite-dimensional circulatory system (1) in order to increase the critical flutter load, looks promising for explaining the peculiarities of optimizing distributed non-conservative structures.However, it makes sense to tackle first not the most general system (1) but rather a particular finite-dimensional nonconservative system with m degrees of freedom.
In this paper, we propose to take an m-link Ziegler pendulum [23][24][25][26] as a toy model for an investigation of the optimal mass and stiffness distributions that give an extremum to the critical flutter load.It appears that even the classical two-link Ziegler pendulum has rarely been studied from this point of view, in contrast to its continuous analogue -the Beck column.The only example of such a study known to the author is contained in the book by Gajewski and Zyczkowski [1].
This paper is organized as follows.In the next Section, we first consider optimization of an undamped two-link Ziegler pendulum.We find an explicit expression for the critical flutter load as a function of the mass or stiffness distribution, and demonstrate that in the space of the two mass coefficients and the flutter load as well as in the space of the two stiffness coefficients and the flutter load, the critical flutter load forms a surface with a self-intersection and with the Whitney umbrella singular point.We consider the problem of optimal mass redistribution and find that the only 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 only minimum, the ratio of the masses is equal to the ratio of the stiffness coefficients.The maxima are attained at the singular cuspidal points of the stability domain and are characterized by a dramatic increase in the critical frequency of the vibrations.Then, we write down the equations of motion of an m-link undamped Ziegler pendulum and consider the case m = 3, in which we again find that the optimal mass distributions maximizing the critical flutter load correspond to vanishing of some of the three point masses.Other types of local extrema are also found that correspond to points where the flutter boundary has a vertical tangent, such as cuspidal points with triple eigenvalues, cf.[11,12], or 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 pendulum and discuss some hypotheses on plausible optimal solutions and their expected properties.

Structural optimization of the Ziegler pendulum
Let us consider the classical Ziegler 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 stiffness coefficient c 1 and damping  [26,32].The point masses m 1 and m 2 are located at the second revolute joint and at the free end of the second rod, respectively.The second rod is subjected to a tangential follower load P [26,32].

Undamped case
Small deviations from the vertical equilibrium for an undamped Ziegler pendulum are described by equation ( 1) with 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 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 the flutter instability and the marginal stability domains Equation ( 5) defines a conical surface in the (c 1 , c 2 , P )-space when m 1,2 and l are fixed: flutter inside the cone, Figure 1(a).However, in the (m 1 , m 2 , P )-space the critical surface (5) looks different and has the 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 (5) we get Defining the new critical flutter load as p := P l/c 2 we come to the more symmetric expression The lower value of the critical load corresponds 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 Figure 1(b).It is seen that the stability boundary has a self-intersection along a ray of the p-axis that starts at the Whitney umbrella singularity with the coordinates (0, 0, 2) in the (m 1 , m 2 , p)-space.Indeed, for small absolute values of m 1 c 2 − m 2 c 1 we can expand the critical load in a series which gives an approximation to the flutter boundary having a canonical form for the Whitney umbrella 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 According to inequality (7), the critical load is always not less than p 0 = 2.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 uniformly distributed mass and stiffness in continuous systems such as the Beck column [19].
In structural optimization problems the uniformly distributed stiffness and mass are usually considered as the initial design that is the starting point in optimization procedures.The critical load of the optimized structure is conventionally compared to that of the same structure with 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 Figure 2, the curves p = p(α) bound the flutter domain shown in 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, mass m 1 is vanishing, then the critical flutter load increases to the value At point B in the stability diagram of Figure 2 the flutter boundary has a vertical tangent, which is a typical phenomenon in non-conservative optimization [27][28][29].
The lower part of the flutter boundary corresponds to designs with a complex conjugate pair of pure imaginary double eigenvalues with the Jordan block; the upper part is the designs with two real double eigenvalues of the same magnitude and different sign, each with the Jordan block.Above the upper flutter boundary lies the domain of statical instability or divergence, with two unstable modes corresponding to two different positive real eigenvalues.At point B, the flutter boundary is tangent to the vertical part of the divergence boundary.
To the right of this vertical line, there is a pair of pure imaginary eigenvalues and a pair of real eigenvalues with the same magnitude and different signs.The transition from the stability boundary to the divergence boundary below point B occurs when a pair of pure imaginary eigenvalues goes out of the origin in the complex plane to infinity, merges there and returns back along the real axis.This happens because at the boundary m 1 = 0, i.e. det M = 0.A similar divergence boundary exists at α = 0, which corresponds to m 2 = 0. Transition through the vertical line above point B is accompanied by another eigenvalue bifurcation at infinity: two real eigenvalues of the same magnitude and different signs go out of the origin in the complex plane in order to merge at infinity and then come back along the imaginary axis.The very point B corresponds to an antagonist of a quadruple zero eigenvalue with the Jordan block, i.e. to a quadruplet of complex eigenvalues that merge at infinity in the complex plane.
To summarize, the popular initial design corresponding to uniformly distributed mass and stiffness turns out to give an absolute minimum to the critical flutter load of the Ziegler pendulum.The critical flutter load attains its local maximum, p B , for m 1 = 0 at the singular cuspidal point B of the stability boundary where the flutter domain has a vertical tangent and touches the boundary of the divergence domain.Note that in [11] a local extremum of the flutter load for the free-free beam carrying a point mass was also found to be at the cuspidal point on the flutter boundary.The global maximum of the critical flutter load for the undamped Ziegler pendulum is at infinity when m 2 = 0.
The global maximum corresponds to a vanishing mass at the free end of the column, which is qualitatively in agreement with the numerically found optimized designs of the Beck column available in the literature [9,[13][14][15][16][17][18].Indeed, all known optimized designs of the Beck column are characterized by the vanishing cross-sections at the free end.Moreover, the gradients of the critical flutter load with respect to the mass or stiffness distribution of the Beck column are large, which is, again, in qualitative agreement with our stability diagram of Figure 2. Most interestingly, with the increase in the critical flutter load the higher and higher modes were reported to be involved into the coupling, which 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.

The m-link Ziegler pendulum
It would be very desirable to extend our study to the case of the multiple-degrees-of-freedom Ziegler pendulum.The corresponding models and recursive schemes for deriving the equations of motion were proposed by Gasparini et al. [23] and by Lobas [25].A three-link Ziegler pendulum was considered by Gallina [24].
We take the linearized equations of Lobas [25], since their model deals with the arbitrary masses and stiffnesses in the joints of an m-link Ziegler pendulum, in contrast to the models of Gasparini and Gallina [23,24].The mass and stiffness matrices of the Ziegler-Lobas model look like For m = 2, matrices ( 12) and ( 13) are reduced to (3).Note that Let us consider the Ziegler-Lobas pendulum with m = 3 links.Now the mass at the free end is m 3 , while the masses m 2 and m 1 are located at the third and second joints, respectively.The first joint connects the first rod with the basement.The length of each of the three rods is equal to l.The stiffness coefficients of the joints are c 3 , c 2 , and c 1 , respectively.
For simplicity we assume that c 1 = c 2 = c 3 = c.Then, the characteristic polynomial has the form with the coefficients 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 for this reason we omit them here.However, numerical experiments evidence that the stability boundary is given by the equation Δ 3 = 0 for the stability condition Δ 3 > 0 implies Δ 2 > 0.
In Figure 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 equal lengths of the links l = 1 and 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 longer 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, which complicates the optimization problem.Nevertheless, such diagrams are convenient for analyzing the geometry of the stability boundary and thus for identifying the 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 Figure 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, Figure 3(a).Both maxima are attained at the cuspidal points of the stability boundary, where it has vertical tangents.However, the stability diagram changes when m 1 = 10, Figure 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.040 347 7, 11.961 144), which corresponds to triple pure imaginary eigenvalues λ ±i1.163 524 3 with the Jordan block of the 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 a 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. Figure 2 and Figure 3(c).This is not surprising because big inertia of the first joint makes the three-link pendulum effectively a two-link pendulum.
The right column in Figure 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 Figure 3(f) with r = 0.1 looks similar to Figure 3(c) and Figure 2.
An increase in 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, Figure 3(d,e).Near these singularities, the stability boundary experiences a sharp turn at points B 1 and B 2 , where the tangent to the boundary is vertical.At such points the eigencurves Imλ(P ) form a crossing that can change either to avoided crossing or to overlapping with the origination of a bubble of complex eigenvalues in dependence on the direction of variation of azimuth angle α.This phenomenon has been observed in many numerical studies of non-conservative optimization problems [4,6,7,9,10,13,14,17,18] and was 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 the finite number of degrees of freedom and the finite number of control parameters the stability boundary of an m-link Ziegler 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, unlike the optimization problems for distributed systems.The possibility to work with singularities related to coalescence of more than two eigenvalues allows us to make a qualitative investigation of the question 'Should low-order models be believed?'[31].Although two-mode approximations work well far from such singularities, in their vicinity we have to take into account higher modes.

Damped case
In the presence of dissipation, the equation of the Ziegler pendulum is with the damping matrix [25,26] For the two-link damped Ziegler pendulum with m = 2, we find the characteristic equation with the coefficients 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] Equation ( 23) 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 was first demonstrated by Bottema already in 1956 [32,35].
In contrast to earlier studies, e.g.[32,[34][35][36], we consider the critical flutter load (22) as a function of the masses P = P (m 1 , m 2 ) for the fixed damping distribution.In Figure 4, the stability diagrams for the damped two-link Ziegler 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 an increase in the number of degrees of freedom and control parameters new types of singularities will appear.The planar stability diagrams of Figure 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 translation of the Latin title of the work by J. Bernoulli published in 1696 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 matrices M and K defined as in (12) and (13).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 }.
We can also consider 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 (19).
We expect that both in the undamped case 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 could be related to infinite eigenvalues λ.
On the other hand, some local extrema should exist with mass distributions that do not contain 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 a finite number of control parameters with a clear physical meaning, the locations of the singularities corresponding to multiple eigenvalues can easily be found numerically with 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 one more degree of freedom is crucial for the correct solution.
Knowledge of rigorously established optimal solutions at every m should shed light on the behavior of the optimal mass distributions and critical frequencies with an increase in the number of degrees of freedom.As a by-product, such a study will give an insight into the problem of dimension reduction, and will serve as a nice playground both for applying the existing methods of nonsmooth analysis and eigenvalue optimization [38][39][40] and for developing them further in the direction of a tighter relation both with singularity theory and with the needs of applications.We expect that the proposed model optimization problem will yield useful recommendations for improving the algorithms for optimizing real non-conservative structures in industry.