Catastrophe in Elastic Tensegrity Frameworks
CCatastrophe in elastic tensegrity frameworks
Alexander Heaton and Sascha TimmeSeptember 29, 2020
Abstract
We discuss elastic tensegrity frameworks made from rigid bars and elastic cables,depending on many parameters. For any fixed parameter values, the stable equilibriumposition of the framework is determined by minimizing an energy function subject toalgebraic constraints. As parameters smoothly change, it can happen that a stableequilibrium disappears. This loss of equilibrium is called catastrophe since the frame-work will experience large-scale shape changes despite small changes of parameters.Using nonlinear algebra we characterize a semialgebraic subset of the parameter space,the catastrophe set , which detects the merging of local extrema from this parametrizedfamily of constrained optimization problems, and hence detects possible catastrophe.Tools from numerical nonlinear algebra allow reliable and efficient computation of allstable equilibrium positions as well as the catastrophe set itself.
Tensegrity structures appear in nature and engineering, scaling in size from nanometers[LHT +
10] to meters [TP03], used on the earth [Mot03, SdO09] and in outer space [Tib02,ZGS + deployable structures [Pel01]. They can significantly changesize and shape, using several different functional configurations during their application.In this article, we discuss elastic tensegrity frameworks (Definition 1) made from rigidbars and elastic cables, similar to those appearing in [LWPQ17, SJM18], but also similarto the tensegrity frameworks defined in [CW96] which are popular in the mathematics andcombinatorics literature. Instead of edge length inequalities as in [CW96], we use Hooke’slaw to introduce an energy function that distinguishes between bars and elastic cables.The configuration is then determined by solving a constrained optimization problem. Thisprovides a large family of simple models which are effectively treated using the theory ofelasticity and energy minimization (see Definition 2). We use numerical nonlinear algebrato calculate all equilibrium positions, in contrast to the more widely-used iterative methods(e.g. Newton-Raphson) which can only find one solution at a time, with no guarantees onfinding them all.Elastic tensegrity frameworks depend on many parameters, e.g., the length of its rigid barsor the fixed position of some nodes. For a given framework we can choose a space of control1 a r X i v : . [ m a t h . M G ] S e p arameters Ω whose values are viewed as the parameters we can manipulate. A path is a mapfrom the unit interval y : [0 , ⊂ R → Ω which describes how the controls y ( t ) vary in time.We use numerical nonlinear algebra to track the changes in stable equilibrium positions of theframework as the control parameters vary. Most importantly, we are interested in a positive-dimensional semialgebraic subset C Ω ⊂ Ω called the catastrophe set (Definition 7). This setrecords those values of control parameters whose crossing could result in a discontinuousjump in the location of the nearest local equilibrium, since the current equilibrium candisappear after crossing C Ω . This loss of equilibrium and the resulting behavior is called a catastrophe . The importance of this set is well-known (see [Arn86] for an overview), but wefind that studying it from the algebraic perspective provides useful benefits.Therefore, the purpose of this article is to show how techniques from numerical nonlinearalgebra can be used to compute the catastrophe set C Ω . For this we introduce an algebraicreformulation in Section 3 that we use to compute a superset D Ω ⊃ C Ω which containsthe relevant information for the original problem (Section 2). This algebraic set D Ω , the catastrophe discriminant , detects the merging of equilibrium solutions from a parametrizedfamily of constrained optimization problems.Hooke’s law provides a simple model which has proven extremely effective in an enormousamount of real-world situations. Also, in the article The Catastrophe Controversy [Guc79]Guckenheimer writes “The application of Catastrophe - Singularity Theory to problems ofelastic stability has been the greatest success of the theory thus far.” Thus catastrophediscriminants are of known importance, but they are very difficult to explicitly compute andthis has limited their usefulness. With the development of efficient techniques in numericalnonlinear algebra, explicit computation of catastrophe discriminants is now within reach.Therefore, another purpose of this article is to explicitly describe these computations for afamily of simple models (elastic tensegrity frameworks) which will be useful in many differentapplications.
Figure 1: Loop crossing the catastrophe set. The black edge is a rigid bar and the green edgesare elastic cables. Square nodes have fixed positions, the cross node is controlled around aloop, and the circular node’s position is determined by minimizing the potential energy inthe green elastic cables.A running example, simple enough to understand yet complicated enough to illustrate theadvantage of knowing C Ω , is Zeeman’s catastrophe machine . Zeeman’s catastrophe machine2onsists of a rigid bar which can rotate freely around one of its endpoints. Attached to thenon-fixed endpoint are two elastic cables. The end of one of the cables is fixed, the other canbe moved freely. The machine and its behavior is depicted in Figure 1 at six discrete-timesnapshots. For more on this example see [PW73], where they give a parametrization of C Ω for a simplified machine, and implicit equations defining C Ω for the actual machine. See also[Arn86, Section 4]. In contrast, we use sample points to encode C Ω not only for Zeeman’smachine but for any elastic tensegrity framework. The basic idea of Zeeman’s machine is tocontrol the free endpoint y ( t ) ∈ Ω (cid:39) R of one cable while the rotating rigid bar settles intoa position of minimum energy. Using numerical nonlinear algebra we can reliably computeall complex solutions to this constrained optimization problem and find among them thereal-valued and stable local minima. In addition, we compute a pseudo-witness set [HS10]for D Ω allowing effective sampling of the catastrophe set C Ω , and therefore easily computableinformation on when catastrophes may occur, and how to avoid them entirely.For those readers new to Zeeman’s machine, consider the behavior depicted in Figure 1.The black bar can rotate around its base, as the green elastic cables pull on its free endpoint.As one of the cable endpoints moves smoothly, the stable equilibrium position of the machinealso moves smoothly... usually. Upon crossing C Ω it can happen that this stable equilibriumdisappears. This forces the machine to rapidly change shape, moving towards some newequilibrium. Without knowledge of C Ω , these behaviors can be very surprising. For example,moving the control point in a small loop does not ensure a return to the original position forthe machine (see Figure 1). Playing with this example, one quickly discovers the advantagesof knowing C Ω . Seemingly random catastrophes become easily predictable.Figure 2: Catastrophe discriminant D Ω (left, degree 72) and catastrophe set C Ω (right) forZeeman’s machine, sampled numerically using homotopy continuation methods.Section 2 gives the basic definitions for elastic tensegrity frameworks. In Section 3 wedescribe an algebraic reformulation of the relevant energy minimization problem. In so doingwe naturally arrive at the equilibrium degree of an elastic tensegrity framework (Definition4), and the catastrophe degree of its catastrophe discriminant (Definition 6). These numbersare intrinsic to the algebraic approach and characterize the algebraic complexity of each3lastic tensegrity framework for a dense set of control parameters. Though the algebraicapproach naturally deals with the algebraic set D Ω , the original problem deals with thesemialgebraic set C Ω (Definition 7). For Zeeman’s machine, both sets are shown in Figure 2.We note that C Ω in Figure 2 is the envelope of a family of curves, each of which is a conchoid ofNicomedes [Kle95, PW73]. Propositions 2 and 3 and Theorem 3 precisely relate the algebraicreformulation with the original setup. In Section 4 we give more details on the requiredcomputations using numerical nonlinear algebra. Finally, in Section 5 we demonstrate ournewly developed tools on a four-bar linkage, which easily becomes an elastic tensegrityframework upon the attachment of two elastic cables (Figure 5). We compute both D Ω and C Ω (Figure 6) and explicitly demonstrate one possible catastrophe (Figure 7). Code thatreproduces all examples in this article can be found at https://doi.org/10.5281/zenodo.4056121 . In this section we formally introduce elastic tensegrity frameworks and the necessary defini-tions and concepts to talk about their equilibrium positions. Let G = ([ n ] , E ) be a graphon [ n ] := { , , . . . , n } nodes and E = B ∪ C edges. Edges are two-element subsets of [ n ].Every ij ∈ B is a rigid bar between nodes i and j and we have (cid:96) ij as its length. Similarly,every ij ∈ C is an elastic cable between nodes i and j that has natural resting length r ij and a constant of elasticity c ij . The graph G is embedded by a map p : [ n ] → R d and wedenote the coordinates of the n nodes of G by p = ( p , . . . , p d ) , . . . , p n ∈ R d and identifythe space of coordinates with R nd . Example 1 (Zeeman’s catastrophe machine) . We illustrate the definitions and concepts ofthis and the next section on Zeeman’s catastrophe machine. Zeeman’s machine is an elastictensegrity framework on n = 4 nodes with edges E = { , , } partitioned as B = { } and C = { , } . See Figure 3 for an illustration. cablebar
14 23
Figure 3: Our setup of Zeeman’s catastrophe machineFor every rigid bar ij ∈ B we define the bar constraint polynomial b ij := (cid:88) k ∈ [ d ] ( p ik − p jk ) − (cid:96) ij (1)and denote by b the polynomial system whose component functions are the b ij for ij ∈ B .4or each elastic cable ij ∈ C we define its potential energy q ij using Hooke’s law q ij := 12 c ij max , (cid:115) (cid:88) k ∈ [ d ] ( p ik − p jk ) − r ij with Q = (cid:88) ij ∈ C q ij . (2)This says that the energy q ij is proportional to the square of the distance the elastic cablehas been stretched past its natural resting length. Though we have only introduced rigidbars and elastic cables, one could easily add compressed elastic edges which want to expandaccording to Hooke’s law. For ease of exposition we proceed with elastic cables and rigidbars, rather than also including compressive struts in our notation.We have introduced several variables. As shorthand we use the symbols p, (cid:96), r, c to referto the variables p ik for i ∈ [ n ] , k ∈ [ d ] (cid:96) ij for ij ∈ Br ij for ij ∈ Cc ij for ij ∈ C. (3)In various examples some of these variables will be viewed as control parameters y ∈ Y whose values we can fix or manipulate at will, while the other variables will be viewed as internal parameters x ∈ X whose values are determined by the controls y and the principleof energy minimization. Often we may fix several control parameters and let others vary insome subset Ω ⊂ Y . Example 2 (Zeeman’s catastrophe machine (continued)) . We continue with Example 1. Wechoose
X, Y as X = { ( p , p ) } (cid:39) R Y = { ( p , p , p , p , p , p , (cid:96) , r , r , c , c ) } (cid:39) R but only consider the subset Ω ⊂ Y as inΩ = (cid:8) (0 , , , − , p , p , , , , . , .
5) : ( p , p ) ∈ R (cid:9) ⊂ Y. In this setup, we have fixed everything except the coordinates of nodes 3 and 4. We control the y = ( p , p ) ∈ Ω and solve for the x = ( p , p ) ∈ X . This means that for a given y = ( p , p ) ∈ Ω we find the coordinates x = ( p , p ) ∈ X which minimize Q ( p , p ) = 14 max (cid:26) , (cid:113) (2 − p ) + ( − − p ) − (cid:27) + 14 max (cid:26) , (cid:113) ( p − p ) + ( p − p ) − (cid:27) restricted to the set { ( x, y ) : b ( x, y ) = 0 }∩ X × Ω. In this case, since B = { } the constraints b ( x, y ) = 0 have only one equation b ( x, y ) = 0 which reads b ( x, y ) = (0 − p ) + (0 − p ) − = 0 . efinition 1. An elastic tensegrity framework is a graph on nodes [ n ] with edges E ⊂ (cid:0) [ n ]2 (cid:1) along with the energy function Q of (2), a partition E = B ∪ C of the edge set into rigidbars and elastic cables, and a partition of variables p, (cid:96), r, c of (3) into internal and controlparameters X and Ω ⊂ Y . A configuration of an elastic tensegrity framework is a tuple( x, y ) ∈ X × Y satisfying the bar constraints (1). Remark 1.
We note that [CW96] used the concept of an energy function as motivationfor their definition of prestress stability. Their definition of a tensegrity framework usesinequalities on edge lengths to distinguish bars from cables and struts. Our definition putsthe energy function at center stage and also allows for a space of control parameters Ω, whichwe need in order to define catastrophe discriminants D Ω ⊂ Ω below.
Definition 2.
We describe the interaction between an elastic tensegrity framework and theenergy function given in (2) with the following definitions.1. Fix a tuple of control parameters y ∈ Y . An elastic tensegrity framework in configu-ration ( x, y ) is stable if the internal parameters x ∈ X are a strict local minimum ofthe energy function Q restricted to the algebraic set { x ∈ X : b ( x, y ) = 0 } of internalparameters satisfying the bar constraints b ( x, y ) = 0 of (1).2. For fixed controls y ∈ Y we collect all strict local minima in the stability set S y ,defined as all internal parameters x ∈ X such that the corresponding elastic tensegrityframework ( x, y ) is stable.3. The stability correspondence SC is the set of pairs ( x, y ) ∈ X × Y such that x ∈ S y .For a given subset Ω ⊂ Y of controls we let SC Ω be all ( x, y ) ∈ X × Ω ⊂ X × Y suchthat x ∈ S y .If we are only interested in a subset of control parameters Ω ⊂ Y these definitions applyverbatim with Ω replacing Y . Example 3 (Zeeman’s catastrophe machine (continued)) . We continue with Example 2.Figure 4 shows Zeeman’s catastrophe machine in a stable configuration. However, for thatspecific value of y the stability set S y contains two points, with the second configurationshown in grey. Since the constraints b ( x, y ) = 0 essentially describe a circle we can also plotthe periodic energy function in Figure 4. For the particular value of the controls y ∈ Ω wechose, there are two local minima, and hence |S y | = 2.In the following, we focus on stable elastic tensegrity frameworks and the behavior whencontrol parameters y ∈ Ω ⊂ Y change. For this, consider a smooth path of control parameters y : [0 , ⊂ R → Ω ⊂ Y (4) t (cid:55)→ y ( t )and an initial condition ( x (0) , y (0)) which is stable according to Definition 2. We are in-terested in the time evolution of the internal parameters x ( t ) determined by minimizing Q constrained by b for the given path y ( t ) of control parameters. In particular, can we iden-tify certain regions where small changes in y ( t ) can cause large changes in the tensegrity6igure 4: Zeeman machine in a stable configuration with |S y | = 2. The second stable positionof node 4 is depicted in gray.framework? In Section 3 we solve an algebraic reformulation of this problem, defining the catastrophe discriminant D Ω ⊂ Ω ⊂ Y and proving that (Theorem 1) as long as y ( t ) / ∈ D Ω we can always lift a smooth path of controls y ( t ) to a smooth path of equilibria, and alsothat (Theorem 2) stable local minima are preserved along this lift. Additionally, we relatethe algebraic reformulation back to the original probem by showing the stronger statement(Theorem 3) that stable local minima are preserved when the smaller, semialgebraic set C Ω ⊂ D Ω is avoided. In this section, we transfer questions about the stability of elastic tensegrity frameworks intoan algebraic problem. The motivation is as follows. The computation of S y is in general avery hard problem since the points in S y are all local minima of a constrained optimizationproblem. Thus, standard optimization methods are not sufficient since they yield in eachrun at most one local minimum and cannot provide guarantees to find all local minima. Incontrast, if we work with systems of polynomial equations we can apply tools from numericalnonlinear algebra to obtain all solutions. This is discussed in more detail in Section 4.In the following, let ([ n ] , E ) be an elastic tensegrity framework with variables p, (cid:96), r, c from (3) partitioned into the internal parameters x ∈ X (cid:39) C m and the control parameters y ∈ Y (cid:39) C m . Compared to the previous section we now work over the complex numbers. LetΩ be an affine subvariety of the control parameters Y we wish to manipulate with controls y ( t ) ∈ Ω. Why an affine subvariety? This allows us, among other things, to consider7ovement of a node constrained to motion in a sphere, perhaps determined by a rigid bar.We denote by Ω R the real part of Ω and assume that the dimension of Ω R and Ω coincide.We introduce variables δ ij for ij ∈ C to eliminate the square roots in the potential energies q ij . For ij ∈ E , let g ij = (cid:40) (cid:96) ij − (cid:80) k ∈ [ d ] ( p ik − p jk ) if ij ∈ Bδ ij − (cid:80) k ∈ [ d ] ( p ik − p jk ) if ij ∈ C and denote by g ( x, δ, y ) : X × C | C | × Y → C | E | the polynomial systems whose component-wiseentries are the g ij . Furthermore, denote by G y the zero set of g for a fixed y ∈ Y G y := { ( x, δ ) ∈ X × C | C | | g ( x, δ, y ) = 0 } . For ij ∈ C let (cid:102) q ij = 12 c ij ( δ ij − r ij ) with (cid:101) Q y = (cid:88) ij ∈ C (cid:102) q ij an algebraic energy function (cid:101) Q y . The subscript emphasizes possible dependency on y ∈ Y .To study the stability set S y we look at the critical points of (cid:101) Q y ( x, δ ) subject to ( x, δ ) ∈ G y .A point ( x, δ ) ∈ G y is a critical point of the energy function (cid:101) Q if the gradient ∇ ˜ Q is orthog-onal to the tangent space of G y at ( x, δ ). If the variety G y is a complete intersection, i.e., thecodimension of G y equals | E | , then we can directly apply the technique of Lagrange multi-pliers to compute the critical points. In the following, we assume for ease of exposition thatthis is the case. However, this is not a critical assumption and the results can be extendedto non-complete intersection by using standard numerical nonlinear algebra techniques forrandomizing overdetermined systems (see [SW05, Chapter 13]). We introduce the variables λ ij for ij ∈ E to act as Lagrange multipliers and let L y ( x, δ, λ ) = (cid:101) Q y + (cid:88) ij ∈ E λ ij g ij ( x, δ, y ) . (5) Definition 3.
Define the polynomial map dL y by letting its component functions be thevarious partial derivatives of L y with respect to x, δ and λ . dL y := ∂L y ∂ ( x, δ, λ ) : X × C | C | × C | E | → X × C | C | × C | E | , ( x, δ, λ, y ) (cid:55)→ dL y (cid:0) x, δ, λ (cid:1) . Let its zero set be the affine algebraic variety denoted L y := dL − y (0) ⊂ X × C | C | × C | E | .Similarly, we define LC := { ( x, δ, λ, y ) | ( x, δ, λ ) ∈ L y } ⊂ X × C | C | × C | E | × Ωand let LC reg denote its open, dense subset of smooth points, LC sing its singular locus, and LC R its real part. Proposition 1.
If the dimension of Ω and LC coincide, then for almost all y ∈ Ω the variety L y is finite and has the same cardinality N . For all y ∈ Ω the variety L y contains at most N isolated points. 8 roof. This is a standard result in algebraic geometry, e.g., [SW05, Theorem 7.1.6].
Definition 4.
Given Ω ⊂ Y we define the equilibrium degree of a framework to be thecardinality of L y for general y ∈ Y . Proposition 1 implies that the equilibrium degree iswell-defined. Example 4 (Zeeman’s catastrophe machine (continued)) . We continue our running ex-ample with edges E = { , , } partitioned as B = { } and C = { , } . Recallfrom Example 2 we had Ω = { (0 , , , − , p , p , , , , . , .
5) : ( p , p ) ∈ R } ⊂ Y ,and X = { ( p , p ) } (cid:39) R . We write x = ( p , p ). The polynomials defining our con-straints are g ( x, δ, y ) = − (0 − p ) − (0 − p ) δ − (2 − p ) − ( − − p ) δ − ( p − p ) − ( p − p ) = . We obtain the Lagrangian of (5) as L ( p ,p ) = 14 ( δ − + 14 ( δ − + (cid:0) − p − p (cid:1) λ + (cid:0) δ − (2 − p ) − ( − − p ) (cid:1) λ + (cid:0) δ − ( p − p ) − ( p − p ) (cid:1) λ . The polynomial system dL y is given by dL ( p ,p ) ( x, δ, λ ) = − λ p − λ (2 − p ) − λ ( p − p ) − λ p − λ ( − − p ) − λ ( p − p ) ( δ −
1) + 2 δ λ ( δ −
1) + 2 δ λ − p − p δ − (2 − p ) − ( − − p ) δ − ( p − p ) − ( p − p ) . The equilibrium degree for this framework is 16. This means that for generic ( p , p ) ∈ Ωthe equations dL ( p ,p ) ( x, δ, λ ) = 0 have 16 isolated solutions.We have particular interest in those parameter values y ∈ Ω where the number of regularisolated solutions |L y | of dL y ( x, δ, λ ) = 0 is less than the equilibrium degree of the framework,since for those parameters local minima can disappear. Definition 5.
Define the catastrophe discriminant D Ω ⊂ Ω as the Zariski closure of the setof critical values of the projection map π : LC → Ω , z = ( x, δ, λ, y ) (cid:55)→ y = π ( z )where the critical values are defined as those π ( z ) ∈ Ω such that either z ∈ LC sing or thereexists a tangent vector v ∈ T z LC in the kernel of dπ z . The catastrophe discriminant is analgebraic subvariety of Ω of codimension 1. Definition 6.
The catastrophe degree of an elastic tensegrity framework is the degree of thealgebraic variety D Ω . 9 xample 5 (Zeeman’s catastrophe machine (continued)) . We continue with Example 4.Refer back to Figure 2 which shows the catastrophe discriminant D Ω ⊂ Ω for Zeeman’smachine with controls Ω defined in Example 2. Note that D Ω does not depend on y ∈ Ωbut just on the choice of X and Ω ⊂ Y itself. Here, D Ω is an algebraic plane curve ofdegree 72. That is, the catastrophe degree is 72. Over the finite field Z the catastrophediscriminant D Ω is the zero set of the 2701-term polynomial p + 13109 p p − p p + 10676 p p + 7407 p p + 4476 p p + 31981 p p +12338 p p − p p + 19319 p p + 4482 p p + . . . − p − p + 540 . Figure 2 also shows the catastrophe set C Ω which we define below. As we move controls y ( t ) ∈ Ω the set C Ω detects changes in the number of local minima, and hence possiblecatastrophe. Definition 7.
We define C Ω := { y ∈ D Ω ∩ Ω R | there exists ( x, δ, λ, y ) ∈ π − ( y ) with δ ≥ } ⊂ D Ω ∩ Ω R to be the catastrophe set . This is the part of the catastrophe discriminant D Ω which relatesto the original problem. We note that the catastrophe set C Ω partitions Ω R into cells withinwhich the number of strict local minima is constant. Figure 2 depicts the number |S y | ofstable local minima for a typical point y in each connected component of the complementΩ R \ C Ω . Look ahead to Figure 6 for another illustration of this phenomenon for the elasticfour-bar linkage discussed in Section 5. Proposition 2.
The catastrophe set C Ω is a semialgebraic set. Proof.
From the definition of C Ω follows that it is the projection of a semialgebraic set whichby the Tarski-Seidenberg principle is again a semialgebraic set.We now begin to prove theorems justifying our interest in D Ω and C Ω . Theorem 1shows that if a smooth path of controls y ( t ) avoids D Ω then there is always a correspondingsmooth path z ( t ) = ( x, δ, λ, y ) ∈ LC . We will combine this with Theorem 2 and Proposition3 to prove Theorem 3, which says that controls y ( t ) avoiding the semialgebraic catastropheset C Ω always correspond to stable local minima, and thus avoid catastrophes where localminima disappear discontinuously. This is called catastrophe since a real-world system wouldbe forced to move rapidly towards the nearest remaining local minima, and since withoutknowledge of C Ω this sudden change in behavior would be very surprising (loss of equilibrium).For the remainder of this section we aim to prove Theorems 1, 2, and 3. Theorem 1.
Let y : [0 , → Ω R with [0 , ⊂ R be a smooth path of control parameters withinitial conditions y (0) ∈ Ω R and z (0) ∈ LC reg such that π ( z (0)) = y (0) and the expecteddimension dim (cid:0) T z (0) LC (cid:1) = dim (cid:0) T y (0) Ω (cid:1) holds. If y ( t ) / ∈ D Ω for all t , then there exists a smooth lifting z : [0 , → LC with π ( z ( t )) = y ( t ) for all t . 10 roof. Since dim (cid:0) T z (0) LC (cid:1) = dim (cid:0) T y (0) Ω (cid:1) and since y (0) / ∈ D Ω we know that the differential dπ z (0) is an isomorphism. By the inverse function theorem, π is a local diffeomorphism at z (0). Hence there is some open neighborhood U of y (0) in Ω mapped diffeomorphically tosome open neighborhood of z (0) in LC . Therefore we can define z ( t ) = π | − U ( y ( t )) for all t such that y ( t ) ∈ U . Since y ( t ) avoids D Ω we know that π | − U ( y ( t )) avoids the singular locusof LC , so that the dimension conditions dim (cid:0) T z ( t ) LC (cid:1) = dim (cid:0) T y ( t ) Ω (cid:1) continue to hold for all t . Since y ( t ) avoids D Ω we also know that dπ z ( t ) will continue to have full rank, since noneof the dim (cid:0) T z ( t ) LC (cid:1) many tangent vectors are in the kernel. This shows we can continuedefining z ( t ) by the local diffeomorphisms that exist by the inverse function theorem. Hencethere exists a lifting z : [0 , → LC with π ( z ( t )) = y ( t ) for all t .We now want to show that avoiding D Ω preserves the stability of the corresponding elastictensegrity framework. For this, we first need a precise definition of what it means to be alocal minimum of (cid:101) Q y ( x, δ ) subject to the constraint G y . Definition 8.
Let d (cid:101) Q and d g ij be the Hessian matrices of (cid:101) Q y and g ij respectively, whenviewed as functions of the variables x and δ . Let dg denote the Jacobian of the constraints g viewed again as functions of the variables x and δ . We say that z = ( x, δ, λ, y ) ∈ LC R is a strict local minimum for the energy (cid:101) Q y subject to constraints g if ( x, δ ) satisfy the sufficientcondition that the projected Hessian V T (cid:34) d (cid:101) Q + (cid:88) ij ∈ E λ ij d g ij (cid:35) V (6)is positive definite. Here V is a real basis of the null space of dg . The term projected refersto the fact that the Hessian is projected onto the tangent space of the constraints g = 0 at( x, δ ). See, e.g., [GMW82, page 81].The next theorem shows that controls y ( t ) avoiding D Ω preserve the stability of thecorresponding elastic tensegrity framework. Theorem 2.
Let y ( t ) be a smooth path of control parameters with initial conditions as inTheorem 1. Furthermore, if the initial condition z (0) is a strict local minimum according toDefinition 8 then all the lifts z ( t ) are strict local minima as well. Proof.
Let V be a matrix whose columns form a basis for Null( dg ), the tangent space of theconstraint variety. Then Null( V T ) = Col( dg T ) is the normal space. Set H := (cid:34) d (cid:101) Q + (cid:88) ij ∈ E λ ij d g ij (cid:35) . As the controls y ( t ) vary smoothly, by Theorem 1 so does the point z ( t ) ∈ LC reg . Thereforethe eigenvalues of the symmetric matrix V T HV also vary smoothly. Since z (0) began as astrict local minimum, V T HV began with all positive eigenvalues. Suppose that at some z ( t )there appears a zero eigenvalue of V T HV . Then there is a null vector V T HV u = 0. Placingparentheses V T ( HV u ) = 0 we see that
HV u ∈ the normal space and by construction V u ∈ w writing HV u in termsof the columns of dg T , and hence ( V u, − w ) ∈ Null( d L ) where d L = (cid:20) H dg T dg (cid:21) is the Hessian of the Lagrangian L y of (5). Note that the null vector ( V u, − w ) of d L extendsto a tangent vector of T z ( t ) LC by appending zeros in the Ω components. This tangent vectorclearly projects to zero by dπ z ( t ) . But this means that y ( t ) ∈ D Ω , completing the proof.We now discuss how our algebraic reformulation relates back to the original problem.In our algebraic reformulation we removed the square roots by introducing the additionalvariables δ ij for ij ∈ C . In the following proposition we assume that all elastic cables are intension since such systems are only structurally stable when self-stress is induced. Proposition 3.
Consider an elastic tensegrity framework in stable configuration ( x, y ) ∈ SC of Definition 2 with y / ∈ C Ω which also satisfies (cid:115) (cid:88) k ∈ [ d ] ( p ik − p jk ) − r ij > ij ∈ C , so that all elastic cables are in tension. Then there exists δ ∈ R | C |≥ and λ ∈ R | E | such that ( x, δ, λ, y ) ∈ LC reg . Proof.
Let V b := { ( x, y ) : b ( x, y ) = 0 } and V g := { ( x, δ, y ) : g ( x, δ, y ) = 0 } . Now considerthe map s : X × Ω → R | C | defined by coordinate functions s ij ( x, y ) := (cid:113)(cid:80) k ∈ [ d ] ( p ik − p jk ) .Restricting this map to V b we have its graph { ( x, s ( x, y ) , y ) : ( x, y ) ∈ V b } ⊂ X × R | C |≥ × Ωwhich provides a local diffeomorphism between V b and V g near any point ( x, y ) ∈ V b satisfying(7). Observe that, by construction, (cid:101) Q y takes values on the image points equal to the valuestaken by Q on the domain V b , provided condition (7) holds. Therefore, if ( x, y ) ∈ V b is astrict local minimum of Q on V b then ( x, s ( x, y ) , y ) ∈ V g is a strict local minimum of (cid:101) Q y on V g . Also, since we assume y / ∈ C Ω we conclude that ( x, s ( x, y ) , y ) ∈ V g is a non-singular pointof V g . Hence, by the Lagrange multipliers condition for local extrema we know that thereexists λ such that ( x, s ( x, y ) , λ, y ) ∈ LC reg , concluding the proof.Finally, we are able to prove that the stability of the corresponding elastic tensegrityframework is preserved by avoiding only C Ω ⊂ Ω. Theorem 3.
Let y : [0 , → Ω with [0 , ⊂ R be a smooth path of control parameters withinitial conditions y (0) ∈ Ω and ( x (0) , y (0)) ∈ SC satisfying the conditions of Proposition 3.If y ( t ) / ∈ C Ω and condition (7) remains satisfied for all t ∈ [0 ,
1] then there exists a smoothlifting z : [0 , → LC with z ( t ) = ( x ( t ) , δ ( t ) , λ ( t ) , y ( t )) for all t such that ( x ( t ) , y ( t )) ∈ SC .12 roof. By Proposition 3, if we have ( x (0) , y (0)) ∈ SC satisfying (7) then there exist δ, λ which we call δ (0) , λ (0) such that ( x (0) , δ (0) , λ (0) , y (0)) ∈ LC reg . But then by Theorems 1and 2 we have lifts z ( t ) = ( x ( t ) , δ ( t ) , λ ( t ) , y ( t )) satisfying Definition 8 as strict local minimafor (cid:101) Q y on V g . Using the graphs of the maps s ( x ( t ) , y ( t )) as in the proof of Proposition 3 wehave local diffeomorphisms for every t such that strict local minima of Q on V b are mapped tostrict local minima of (cid:101) Q y on V g , since (7) is satisfied for all t . But then each ( x ( t ) , y ( t )) ∈ SC as required. In this section we use the algebraic reformulation developed in the previous section to de-scribe numerical nonlinear algebra routines which can be used to answer the following threecomputational problems:1. Given γ ∈ Ω R compute S γ .2. Given an algebraic path y ( t ) : [0 , → Ω R ⊂ Y and an initial configuration x ∈S y (0) compute the path points γ ⊂ y ([0 , C Ω .We start with the first question. Recall that dL y is a polynomial system. We cancompute all isolated solutions of a polynomial system using homotopy continuation methods[SW05]. Homotopy continuation methods work by first solving a compatible but simplerstart system and then keeping track of these solutions as the start system is deformed intothe system we intended to solve originally (the target system). For our computations we usethe software package HomotopyContinuation.jl [BT18]. To compute S γ for a given γ ∈ Ω R we therefore first solve dL γ ( x, δ, λ ) = 0 which results in finitely many complex solutions L γ .Of these complex solutions we then select those solutions whose components are real-valuedand then further select those real-valued solutions where the projected Hessian defined in(6) is positive definite. Note that computing solutions to dL γ ( x, δ, λ ) = 0 usually requiresthat we track many more paths than the equilibrium degree of L γ . If the goal is to compute S y for many different y ∈ Ω R it is more efficient to use a parameter homotopy [MS89, SW05].There, the idea is to first compute L y for a general (complex) y ∈ Ω and then to usethe parameter homotopy H ( z, t ) = dL ty +(1 − t ) y ( z ) to efficiently compute L y . Using thisparameter homotopy approach allows us to track only the minimal numbers of paths thatstill guarantee all solutions mathcalS y are computed correctly.Consider the second question where we are given an algebraic path y ( t ) : [0 , → Ω R ⊂ Y and an initial configuration x ∈ S y (0) . We want to compute the path points γ ⊆ y ([0 , C Ω and y ([0 , D Ω ∩ α where α ⊂ Ω is an algebraic curve containing y ([0 , α (cid:54)⊂ D Ω . The catastrophe discriminant D Ω is given13y π ( H − (0))) with π from Definition 5 and H Ω ( x, δ, λ, y ) = (cid:20) dL y ( x, δ, λ )det d L y ( x, δ, λ ) (cid:21) . Since the evaluation of a determinant is numerically unstable it is better to instead usethe formulation that there exists a v ∈ P n such that d L y ( x, δ, λ ) · v = 0. Consider thecollection { H Ω , π, π − ( α ) , W} where H Ω and π are the polynomial maps defined above and W = π − ( α ) ∩ H − (0) contains finitely many solution points. In the case that α is a line thisis known as a pseudo-witness set [HS10] since it allows us to perform computations on D Ω without knowing its defining polynomials explicitly. Since W is the zero set of a polynomialsystem it can again be computed by using homotopy continuation techniques. If α is a linethen cardinality of W is the catastrophe degree of the tensegrity framework. To compute C Ω ∩ y ([0 , W we have to select from ( x, δ, λ, γ ) ∈ W all those solutions which havereal-valued coordinates, δ >
0, and γ ∈ y ([0 , ⊆ α .We move to the third question and discuss the computation of the catastrophe set C Ω .This is more involved since C Ω is a positive-dimensional set and we have to decide what“compute” means in our context. Since C Ω is a semialgebraic set it can theoretically bedefined by a union of finite lists of polynomial equalities and inequalities. However, comput-ing the describing polynomials is a very challenging computational problem since it requiresGr¨obner bases computations which have exponential complexity. We were able to computethe polynomial defining D Ω in Example 5, but only over a finite field, and larger exampleswill likely fail to terminate. Instead we opt to obtain a sufficiently dense point sample of C Ω . The idea is to apply the previously described technique to compute repeatedly the in-tersection of C Ω and a real line (cid:96) ⊂ Ω R . To proceed we first compute a pseudo-witness set { H Ω , π, π − ( (cid:96) ) , W } for a general (complex) line (cid:96) ⊂ Ω and then we can compute the pseudowitness set { H Ω , π, π − ( (cid:96) ) , W} by utilizing a parameter homotopy. As discussed above, thisis much more efficient for the repeated solution of our equations. Note that even if the reallines (cid:96) ⊂ Ω R are sampled uniformly, this does not guarantee that the obtained sample pointsconverge to a uniform sample of C Ω . If uniform sampling is of interest the procedure can beaugmented with a rejection step as described in [BM20]. The outlined procedure is an effec-tive method to sample the catastrophe discriminant D Ω . Figure 2 depicts the point samplesobtained for Zeeman’s catastrophe machine using this method, while Figure 6 depicts thoseobtained for the elastic four-bar framework of Section 5. We want to demonstrate the developed techniques on another example. For this we considera planar four-bar linkage which is constructed from four bars connected in a loop by fourrotating joints where one link of the chain is fixed. The resulting mechanism has one degreeof freedom. Four-bar linkages are extensively studied in mechanics as well as numerical non-linear algebra [WS11]. Here, we extend a four-bar linkage to an elastic tensegrity frameworkby introducing two nodes which are attached to the two non-fixed joints by elastic cables.Formally, we introduces six nodes with coordinates p , . . . , p ∈ R , bars B = { , , , } and elastic cables C = { , } . See Figure 5 for an illustration of this basic setup.14 cablebar Figure 5: Setup of a four-bar elastic tensegrity framework.The zero set of the bar constraints b ij , ij ∈ B , is a curve of degree 6 which can beparameterized by the plane curve traced out by the motion of the midpoint ( p + p ) /
2. Inkinematics terminology the midpoint is a coupler point and the plane curve is called the coupler curve of the mechanism .The idea is to fix nodes 1, 2 and 5, and control node 6. For our model this means choosing X = { ( p , p , p , p ) } (cid:39) R as internal parameters and Ω = { ( p , p ) } (cid:39) R as controlparameters. Furthermore, we fix nodes p = ( − , p = (1 , p = (4 , l = 3, l = 1, l = 1 .
5, resting lengths r = r = 0 . c = 1, c = 2.In this setup the framework has an equilibrium degree of 64. The resulting catastrophediscriminant D Ω is a curve of degree 288. D Ω and the catastrophe set C Ω are depicted inFigure 6. The typical sizes of the stability set S γ , γ ∈ Ω, are 2, 3 and 4.Figure 6: The catastrophe discriminant (left) and catastrophe set (right) of the elastic fourbar framework. The cardinality of the stability set for points in each chamber of the com-plement of the catastrophe set is shown in the upper right corner.Finally, we also want to give in Figure 7 another concrete example of a catastrophe.There, the control node 5 is depicted by a cross and it is dragged in a straight line between15ts position in the left figure and its position in the right figure. When the control nodecrosses the catastrophe set C Ω , its previously stable position disappears from S y , and theframework “jumps” to a new position. Again, without knowledge of C Ω these catastrophesare extremely surprising. With knowledge of C Ω and Theorem 3 they become avoidable.Figure 7: Left: The elastic four bar framework in a stable configuration. Right: Configuration of theframework after crossing the catastrophe set. The gray dashed line is the coupler curve of the four barlinkage traced out by the coupler point defined as the midpoint of the bar connecting node 2 and 3. Thecoupler curve allows to parameterize all possible four bar positions. The catastrophe set C Ω is depicted inred. At the bottom are the energy landscapes along the coupler curve with the current position depicted ingreen. This article described elastic tensegrity frameworks as a large family of simple models basedon Hooke’s law and energy minimization. For this family we showed how to explicitlycalculate and track all stable equilibrium positions of a given framework. More importantly,we showed how to calculate the catastrophe set C Ω by using pseudo-witness sets to encode asuperset D Ω ⊃ C Ω . To do this we reformulated the problem algebraically to take advantageof tools in numerical nonlinear algebra. Knowing the catastrophe set provides extremelyuseful information, since Theorem 3 shows that paths of control parameters avoiding C Ω willalso avoid discontinuous loss of equilibrium, and hence avoid surprising large-scale shapechanges.In our two illustrative examples, we chose the controls Ω as a two-dimensional space16verlaid with the configuration itself. These choices were made to demonstrate the ideas.However, the calculation and tracking of all stable local minima by parameter homotopyand the encoding of D Ω ⊃ C Ω by pseudo-witness sets apply much more generally. Thecontrol set Ω can be chosen in any way, and all the same methods apply, even if thereare no easy visualizations for the controls desired. Therefore, for more complicated sets ofcontrol parameters Ω it is of interest to develop more efficient local sampling techniquesbased on Monte Carlo methods [ZHCG18], perhaps only sampling C Ω locally near the initialconfiguration ( x (0) , y (0)) or locally near the intended path y ([0 , C Ω nearest to a given initial or current position y ( t ).We would also mention recent work [BHM +
20] which details a sampling scheme whosegoal is to learn the real discriminant of a parametrized polynomial system, as well as thenumber of real solutions on each connected component. They combine homotopy continua-tion methods with k -nearest neighbors and deep learning techniques. For elastic tensegrityframeworks, these techniques might be used to learn D Ω ∩ Ω R .Finally we discuss the potential of our results for use in mechanobiology [IWS14], wherescientists have frequently and successfully used tensegrity to model cell mechanics. Evensmall and simple elastic tensegrity frameworks (e.g. with 6 or 12 rigid bars, plus morecables) have been used to explain and predict experimental results observed in actual cellsand living tissue [DSLB +
11, VVB00, CS03, WTNC + qualitative phenomena observedin actual experiments could be predicted or explained by elastic tensegrity frameworks.Knowing the catastrophe set for a simple tensegrity model with a biology-informed choiceof Ω would give catastrophe predictions that could then be tested experimentally. References [Arn86] Vladimir I. Arnold.
Catastrophe Theory . Springer Verlag, Berlin, 1986.[BHM +
20] Edgar A. Bernal, Jonathan D. Hauenstein, Dhagash Mehta, Margaret H. Re-gan, and Tingting Tang. Machine learning the real discriminant locus, 2020.[BM20] Paul Breiding and Orlando Marigliano. Random points on an algebraic mani-fold.
SIAM Journal on Mathematics of Data Science , 2(3):683–704, 2020.[BT18] Paul Breiding and Sascha Timme. HomotopyContinuation.jl: A Package forHomotopy Continuation in Julia. In
Mathematical Software – ICMS 2018 ,pages 458–465. Springer International Publishing, 2018.[Cal78] Christopher R. Calladine. Buckminster Fuller’s “Tensegrity” structures andMaxwell’s rules for the construction of stiff frames.
International Journal ofSolids and Structures , 14(2):161 – 172, 1978.17CS03] Mark F. Coughlin and Dimitrije Stamenovi´c. A Prestressed Cable NetworkModel of the Adherent Cell Cytoskeleton.
Biophysical Journal , 84(2):1328 –1336, 2003.[CW96] Robert Connelly and Walter Whiteley. Second-order rigidity and prestressstability for tensegrity frameworks.
SIAM J. Discrete Math. , 9(3):453–491,1996.[DSLB +
11] Gianluca De Santis, Alex B. Lennon, Federica Boschetti, Benedict Verhegghe,Pascal Verdonck, and Patrick J. Prendergast. How can cells sense the elasticityof a substrate? An analysis using a cell tensegrity model.
Eur Cell Mater ,22:202–213, Oct 2011.[GMW82] Philip E. Gill, Walter Murray, and Margaret H. Wright.
Practical optimization .Emerald Group Publishing Limited, 1982.[Guc79] John Guckenheimer. The catastrophe controversy.
Math. Intelligencer , 1(1):15–20, 1978/79.[HS10] Jonathan D. Hauenstein and Andrew J. Sommese. Witness sets of projections.
Applied Mathematics and Computation , 217(7):3349 – 3354, 2010.[IWS14] Donald E. Ingber, Ning Wang, and Dimitrije Stamenovi´c. Tensegrity, cellularbiophysics, and the mechanics of living systems.
Rep Prog Phys , 77(4):046603,Apr 2014.[Kle95] Felix Klein.
Vortr¨age ¨uber ausgew¨ahlte Fragen der Elementargeometrie . Teub-ner, Leipzig, 1895.[LHT +
10] Tim Liedl, Bj¨orn H¨ogberg, Jessica Tytell, Donald E. Ingber, and William M.Shih. Self-assembly of three-dimensional prestressed tensegrity structures fromDNA.
Nature Nanotechnology , 5(9):520–524, 2010.[LWPQ17] Ke Liu, Jiangtao Wu, Glaucio H. Paulino, and H. Jerry Qi. Programmabledeployment of tensegrity structures by stimulus-responsive polymers.
ScientificReports , 7(1):3511, 2017.[Mot03] Ren´e Motro.
Tensegrity: Structural systems for the future . Butterworth-Heinemann, Oxford, 2003.[MS89] Alexander P. Morgan and Andrew J. Sommese. Coefficient-parameter poly-nomial continuation.
Applied Mathematics and Computation , 29(2):123–160,1989.[Pel01] Sergio Pellegrino.
Deployable Structures . Springer Verlag, Vienna, 2001.[PW73] Tim Poston and Alexander E. R. Woodcock. Zeeman’s Catastrophe Ma-chine.
Mathematical Proceedings of the Cambridge Philosophical Society ,74(2):217–226, 1973. 18SdO09] R.E. Skelton and M.C. de Oliveira.
Tensegrity Systems . Springer-Verlag US,Boston, 2009.[SJM18] Menachem Stern, Viraaj Jayaram, and Arvind Murugan. Shaping the topologyof folding pathways in mechanical systems.
Nature Communications , 9(1):4303,2018.[SW05] Andrew J. Sommese and Charles W. Wampler.
The Numerical Solution ofSystems of Polynomials Arising in Engineering and Science . World Scientific,2005.[Tib02] Gunnar Tibert.
Deployable Tensegrity Structures for Space Applications . PhDthesis, KTH Royal Institute of Technology, 2002.[TP03] Gunnar Tibert and Sergio Pellegrino. Deployable tensegrity masts. In . American Institute of Aeronautics and Astronautics, 2003.[VVB00] Konstantin Y. Volokh, Oren Vilnay, and M. Belsky. Tensegrity architectureexplains linear stiffening and predicts softening of living cells.
J Biomech ,33(12):1543–1549, Dec 2000.[WS11] Charles W. Wampler and Andrew J. Sommese. Numerical algebraic geometryand algebraic kinematics.
Acta Numerica , 20:469–567, 2011.[WTNC +
02] Ning Wang, Iva Marija Toli´c-Nørrelykke, Jianxin Chen, Srboljub M. Mi-jailovich, James P. Butler, Jeffrey J. Fredberg, and Dimitrije Stamenovi´c. Cellprestress. i. stiffness and prestress are closely associated in adherent contrac-tile cells.
American Journal of Physiology-Cell Physiology , 282(3):C606–C616,2002. PMID: 11832346.[ZGS +
12] V.S. Zolesi, P.L. Ganga, L. Scolamiero, A. Micheletti, P. Podio-Guidugli,G. Tibert, A. Donati, and M. Ghiozzi. On an innovative deployment conceptfor large space structures. In . American Institute of Aeronautics and Astronautics, 2012.[ZHCG18] Emilio Zappa, Miranda Holmes-Cerfon, and Jonathan Goodman. Monte Carloon manifolds: sampling densities and integrating functions.
Communicationson Pure and Applied Mathematics , 71(12):2609–2647, 2018.[ZO15] Jing Yao Zhang and Makoto Ohsaki.
Tensegrity Structures: Form, Stabilityand Symmetry . Springer Japan, Boston, 2015.
Authors’ addresses:
Alex Heaton, Technische Universit¨at Berlin and Max Planck Institute for Mathematics in the Sciences, [email protected], [email protected]
Sascha Timme, Technische Universit¨at Berlin, [email protected]@math.tu-berlin.de