A continuation-based method for finding laminated composite stacking sequences

A method of recovering laminate ply stacking sequences from a set of up to twelve lamination parameters using polynomial homotopy continuation techniques is presented. The ply angles are treated as continuous variables, and are allowed to take any value between −90◦ and +90◦. The individual plies are assumed to be orthotropic and have constant stiffness. The method is fully deterministic, and does not rely on an optimisation process to establish the stacking sequence. Polyhedral continuation methods are used to limit the solution space in which the stacking sequences are sought. The method can reliably find every stacking sequence solution that exists to achieve a precisely specified set of lamination parameter “targets”, with the number of real solutions to a feasible combination of target properties found to vary from 1 to over 100. The same method is also demonstrated to be able to find stacking sequences to satisfy a set of specified ABD stiffness matrix terms, as might be required following a direct-stiffness modelling design process.


Introduction and background
Every year, thousands of undergraduate students in Engineering learn how to predict the elastic behaviour of thin laminated plates constructed from any number of individually orthotropic laminae using Classical Lamination Theory (CLT). Given the thickness, material properties, location within the lamina stack, and the angle between the principal axes of each lamina (or ply) and a pre-defined axis system in the plane of the plate, it is possible to produce what are known as the extensional, coupling, and bending stiffness matrices [1], usually combined together into the 6 × 6 ABD matrix. The teaching of this process is seldom accompanied by any discussion of the far more difficult "inverse" problem of finding an appropriate set of plies and an accompanying arrangement thereof needed to produce a laminate with a desired set of elastic properties.
The aim of this paper is to capitalise on the advances made possible in the field of kinematics by a polynomial root finding method called polynomial homotopy continuation (PHC) by adapting the method to the (in many ways analogous) task of laminate stacking sequence design. The question the paper seeks to answer is: given a set of either lamination parameters or laminate stiffness values, how many different laminates having exactly these properties can be constructed, and what are their stacking sequences?
Besides improvements in specific strength and stiffness, one of the main advantages of employing non-isotropic composite materials is the possibility of tailoring the structural response, such as bend-twist coupling in composite turbine blades [2]. In principle, any property of the laminate could be altered to achieve the desired result, including the individual ply material properties and thicknesses. However, the majority of the literature regarding laminate design is concerned with the varying of ply angles and the sequence in which they are stacked to produce laminates with particular stiffness, buckling, dynamic, and thermal response behaviour.
Traditionally, only a discrete set of fibre orientation angles have been considered in the laminate design process. The angles 0 • , ±45 • and 90 • are particularly common, and stem from a time when the lay-up process was very manual. This restriction has bred a family of laminate modelling methods in which the ply fibre angles are treated as discrete, rather than continuous variables. Genetic algorithms have become popular in discrete fibre angle stacking sequence optimisation [3]. More recent manufacturing techniques such as automated fibre placement (AFP) allow laminates with any combination of ply fibre angles to be manufactured with a high degree of placement accuracy. These "nonconventional" laminates [4] allow greater flexibility in the design of laminated composites [5,6]. Techniques that treat ply angles as continuous variables, such as the one presented in this paper, do now have a role to play in the search for practical stacking sequences.
The use of lamination parameters in composite plate design is a two-step process. In the first, an optimisation objective function for the plate's stiffness, buckling strength, natural frequency, or other property or properties is constructed using the lamination parameters as optimisation variables. In the second step, the resulting parameters are used to establish a ply stacking sequence whose properties closely approximate those of the optimised plate. Two of the main advantages of using lamination parameters are that the laminate stiffnesses can be expressed as linear combinations of the twelve lamination parameters [7], and the feasible region of lamination parameters is convex [8]. These features make lamination parameters well suited to optimisation problems involving composite plates [9]. However, while the feasible space of lamination parameters is known to be convex, defining the bounds on each individual parameter remains challenging. Much of the time and effort invested in the research of lamination parameters in the last 20 years has been given to determining computationally efficient means of defining feasible domains, both with [10] and without [11] restrictions being placed on the possible ply orientations. A further challenge lies in the fact that once the optimal lamination parameters have been found, a viable stacking sequence must still be obtained. Again, most authors limit themselves to discrete ply angles, with some approaches that have been demonstrated including simple least-squares fitting, various forms of genetic algorithms [12,13,14], ant colony and particle swarm optimisation [15], fractal branch and bound methods [16], and logic-based branch and bound techniques [17]. For smaller sets of parameters it may even be possible to find explicit expressions for stacking sequences. Lipton [18] was able to find such expressions for three-layered laminates derived from the four in-plane lamination parameters. However, for larger sets of lamination parameters, the problem becomes less tractable.
While the majority of stacking sequence recovery algorithms use a fairly restricted set of ply angles, it is not unheard of for nonconventional angle ply laminates to be allowed when extracting stacking sequences from polar parameters [19], although in [19] the ply angles are restricted to integers between −90 • and +90 • .
Besides lamination parameters, the stiffness terms contained within a laminate ABD matrix are sometimes also used in the laminate design process in lieu of using the ply angles themselves. Maintaining feasibility in this "direct stiffness" modelling process is also challenging, and the task of establishing a stacking sequence to obtain the desired ABD values also waits at the end. This paper presents a method of extracting stacking sequences given a set of feasible stiffness requirements or lamination parameters using PHC [20]. The method requires the ply properties to be treated as continuous variables. In principle, any of the ply properties could be treated as unknowns in the process, but attention will be restricted here to the ply angles. It will also be assumed that the stiffness or lamination parameter targets used are already known to be feasible, as the topic of maintaining feasibility is addressed sufficiently elsewhere [18,21,22,10,23,24]. The key advantage of continuation, or homotopy, approaches to solving systems of polynomial equations is that they can be shown to reliably find all the solutions that exist (counting multiplicity) [25], something no optimisation-based approach or root finding technique like Newton's Method can be guaranteed to do. Once the full set of possible laminate stacking sequences has been found, a designer can then exercise their engineering judgement in selecting the most appropriate one. The metrics employed in making such a decision are also not within the scope of this paper.
Attempts to efficiently solve multivariate polynomial systems using computational methods began in the 1970s, and have led to the development of a range of new approaches and accompanying implementations in software. The clear potential for these methods to solve engineering problems has meant they have been brought to bear in a number of applications outside of the field of algebraic geometry [26]. Such applications include chemical equilibrium [27], the synthesis of both planar [28] and spatial mechanisms [29,30,31], and even the study of the buckling of composite shells of certain geometries [32]. In addition to continuation there are, of course, a number of other methods that have been developed for solving multivariate polynomial systems, such as Gröebner bases and resultant matrices [33], and many have been applied to systems in which the variables represent unknown angles [34]. The adaptability of continuation to systems with larger numbers of unknowns, and perhaps more importantly, the ease with which a non-specialist can quickly start to solve polynomial based problems with freely available and mature software packages makes continuation the most appealing tool for this task.
Continuation methods can be applied to solving a variety of families of nonlinear equations, but polynomial systems have received much of the attention. In order to find the roots of a system of multivariate polynomial equations (the target system), PHC methods begin with the construction of a structurally similar system (the start system) for which a full set of solutions is already known. A numerical process called a coefficient homotopy then follows in which the coefficients of the start system are gradually changed into those of the target system. This process occurs over a number of discrete increments, with the solutions of each interim system being computed using the solutions from those at the previous increment in a process called path following. If the start system has been chosen wisely, it is possible to locate every solution of the target system using a coefficient homotopy. The process can present some numerical challenges, with the target system solutions being regular (either real or appearing in complex conjugate pairs), singular, solutions at infinity, or belonging to positive dimensional solution sets [35]. One advantage of PHC methods is that they naturally lend themselves to parallelisation, and the path following of each separate start system solution can be carried out in an independent thread.
The various published procedures for writing PHC algorithms differ mainly in their approaches to constructing start systems. The choice of start system is extremely important as it establishes the number of solutions that must be tracked in the path following process. The total degree of multivariate polynomial solutions grows very quickly with number of equations and the degree of each equation. It also happens that the majority of the solutions for many target systems are not finite, and not of practical interest. An appropriate choice of start system can reduce the number of solutions that need to be tracked by excluding those that would diverge to infinity during the coefficient homotopy step. Some software packages generate efficient multi-homogeneous start systems [36], while others use polyhedral methods [37,38]. A more thorough discussion of the advantages of each approach is beyond the scope of this paper, however, only polyhedral methods will be used here because of the tight upper bound on the number of finite solutions they provide, and because they can be made more efficient in cases in which some of the target system polynomial equations have the same structure [39,40]. The numerical examples presented in this paper were solved using a combination of PHCPack [38], and Matlab [41] code based on the method described in [40]. Matlab files for generating PHCpack input files for sets of lamination parameter or plate stiffness targets, including the scripts used in generating the numerical examples in Section 3 below, accompany this paper. A link is provided in the Data Availability Section at the end.

Laminate equations in polynomial form
Both lamination parameters and the ABD stiffness terms from CLT can be expressed in terms of trigonometric functions of the ply angles, θ k , for layers 1 to N within the laminate. Clearly the laminate equations are not in pure polynomial form, but this is a familiar situation to anyone accustomed to modelling multi-link linkages, such as might be used in a robotic arm. The traditional approach is to simply replace any trigonometric terms with new variables: and to append a new coupling equation to the system for each unknown angle: Another approach involves using the tan half-angle formulae to replace the trigonometric functions with a single unknown for each angle. This has the advantage of not requiring additional coupling equations, but comes at the expense of raising the degree of each existing equation. The choice of substitution variables can have a significant impact on the efficiency of the continuation process that follows. In this case, it is suggested that the following variables be used in converting the laminate equations to polynomial form: Because of the symmetry that is present in many laminate design problems, it is not uncommon for only the u k terms to be present in the system. In cases in which the t k terms are required, the two new variables may be linked with an additional coupling equation: While PHC methods can be used to identify and explore systems possessing more unknowns than equations, as well as systems containing positive dimensional solutions (lines, surfaces, or hypersurfaces in the space formed by the polynomial unknowns that satisfy the equations) [42], this paper is concerned only with square systems (same number of equations as unknowns), and the regular finite, or zero dimensional, solutions thereof. This has some limitations, as the maximum number of independent laminate stiffness equations that can be assembled is twelve, the same as the maximum number of lamination parameters. This means that, at most, twelve unknowns can be sought in any one calculation. However, it is often surprising how many different solutions to such square systems can be found. The twelve unknowns could consist of any combination of non-discrete characteristic of the plies (such as thickness, volume fraction, or ply angle), but in the interest of clarity, only the use of ply angles as unknowns will be addressed here.
In terms of the variable substitutions in Equation 1, the four trigonometric expressions encountered in both lamination parameters and the ABD stiffness matrix elements are: Using these substitutions, a system of laminate equations in polynomial form can be constructed: where p i are the laminate equations to be solved and n is the number of unknown ply angles. Each of the n equations represents a separate laminate property "target", which could be a lamination parameter value, an ABD stiffness matrix term, or some combination of each. Clearly this system will be highly nonlinear in the general case, but its polynomial form allows the efficient finding of all possible solutions using continuation methods.
The following assumptions will be assumed to hold in the remainder of the paper: • all plies are of equal thickness, are orthotropic, and have the same material properties; • at least one stacking sequence capable of achieving the targets exists (the target combination is known to be feasible); • ply angles are continuous variables, with any angle between −90 • and +90 • being achievable; • for simplicity, only symmetric and antisymmetric laminates with an even number of plies are considered, with unsymmetric laminates allowed to have any number of plies.
Note that laminates with an odd number of plies could be treated by splitting the central ply into two of half the original's thickness.

Lamination parameters
The in-plane, coupling, and out-of-plane lamination parameters are defined as: (cos 2θ(z), sin 2θ(z), cos 4θ(z), sin 4θ(z)) dz (V1D, V2D, V3D, V4D) = 12 1 2 − 1 2z 2 (cos 2θ(z), sin 2θ(z), cos 4θ(z), sin 4θ(z)) dz respectively, where z is the through-thickness direction measured from the laminate mid-surface,z = z/h is the normalised through-thickness coordinate, and h is the laminate thickness. In cases in which the ply angle remains constant within each layer, the integral form above can be replaced with a summation: where the h k are the through-thickness locations of the lamina surfaces. Note that there is no universally accepted notation for lamination parameters. Some definitions also exchange the position of the cos 4θ and sin 2θ terms above.
Using the substitutions in Equation 3, the equations in 5 become: The equations in 6 are in pure polynomial form in terms of the t k and u k . The assumption that all the plies in the laminate have the same properties leads to expressions for the lamination parameters that are independent of any material properties. If, for example, it is required that a laminate have V 4A = 0.3, the equation: must be present in the target system (Equation system 4). The number of plies (N ) may not be the same as the number of independent ply angles (n) as there may be some repetition of the ply angles, such as might occur in a symmetric laminate, meaning n ≤ N . A similar situation can arise if there is a consistent enforced relationship between certain ply angles, as in the sign reversal in mirroring plies in an antisymmetric laminate.

ABD matrix terms
The lamina stiffness matrix terms (Q) for the k th layer of a laminate can be expressed in terms of the invariant stiffnesses U 1 to U 5 as [1]: Again using the substitutions in Equation 3 this becomes: The laminate stiffness matrices can then be assembled in the usual way: If, for example, it is required that a laminate have B 16 = 100 N, the equation: must be present in the target system (Equation system 4). As mentioned above, it is possible to construct the ABD stiffness terms as linear combinations of lamination parameters, but the formulation given in Equation 9 will be used directly in the numerical examples below when such stiffness terms are specified as design targets. Much like differential equation solving algorithms can experience numerical instability if the magnitudes of the coefficients in the differential equations are sufficiently dissimilar, PHC software can also fail in the path tracking or endgame refinement process if the scale of the monomial coefficients are very different from one another. Most PHC software includes some pre-scaling process in which all the coefficients of a particular equation are multiplied by a scalar to bring them into line with the magnitude of coefficients present in the other equations, or sometimes temporary variable substitutions are performed to achieve the same effect. However, when constructing laminate equation target systems for many typical material types it is often better to use units of mm, kN and MPa, rather than m, N and Pa. This minimises the work a pre-scaler must do, especially when mixtures of in-plane and out-of-plane stiffness targets are used.

Non-singular target systems
There are occasionally certain combinations of independent variables for which the Jacobian matrix of the system in Equation 4 becomes rank deficient (known as singular points or singular solutions). This complicates the continuation process as the final stages of the path following algorithm can not necessarily benefit from the quadratic convergence present around regular points. Of greater concern are cases in which the equations in the target system are not independent, effectively leading to an underdetermined system for which the Jacobian matrix is always singular. There are also some targets for which it is certain that no real stacking sequence solutions exist, such as a non-zero target in the coupling stiffness matrix for a symmetric laminate. In the context of selecting a combination of targets that could potentially have at least one real stacking sequence solution, the term "feasible combination" will be used below. Whether the combination is feasible insofar as the numerical values of the targets allow for a real stacking sequence solution to exist is another matter.
The task of identifying feasible combinations of lamination parameters is relatively straightforward in that all twelve are already known to be independent. It simply remains to identify those lamination parameters known to be zero in symmetric and antisymmetric matrices and avoid setting these as nonzero targets. The allowable combinations are shown in Table 1. Any subset of targets from each column in the table can be used to attempt to find a stacking sequence. Symmetric laminates can be characterised by up to eight lamination parameters, antisymmetric laminates by six, while laminates with no symmetry may require all twelve to be fully specified. This means that the largest lamination parameter target system will consist of 24 independent polynomial equations in 24 unknowns (recalling that there are two variables, t k and u k , introduced for each unknown ply angle).
Target systems for ABD stiffness matrix targets may also have up to twelve target stiffness terms. Also like lamination parameters, there are some situations in which the terms are identically zero, such as is the case for the coupling stiffness terms for all symmetric laminates. However, unlike lamination parameters, not all combinations of the stiffness matrix terms are independent. Of the 18 upper triangular elements of the three stiffness matrices, the targets that are not identically zero and the combinations in which they can be used to form independent target systems are given in Table 2. Combinations of target stiffnesses may be assembled from only one set from each column in the table.  26 ] is not as the three diagonal elements of the A matrix are not independent. As a further example, an unsymmetric laminate with a target system based on the stiffnesses [A 26 , B 12 , B 66 , D 11 , D 12 , D 22 ] will not be independent because both the combination of B 12 with B 66 , and the combination of D 11 and D 12 with D 22 will make the Jacobian of the target system rank deficient. Again, symmetric laminates can be characterised by up to eight stiffness values, antisymmetric laminates by six, while laminates with no symmetry may require all twelve to be fully specified.

Mixed volumes of laminate target equations
The mixed volume [43] of a system of polynomials provides a tight upper bound on the number of finite solutions that can be expected to exist. This limit is sometimes called the "Bernstein count" or "BKK bound", and is often substantially smaller than the total degree of the system. The advantage of using polyhedral continuation methods in solving systems of polynomial equations is that a start system with only as many solutions as the integer mixed volume can be constructed (using methods such as those described in [40]), and then only these solutions are tracked across to the target system. Tracking all the solutions from a total-degree start system is likely to be much more computationally expensive.
The set of all the exponents present in an equation is referred to as the equation's "support". The magnitude of the mixed volume of a system is a function of the shapes of the Newton polytopes formed by the supports of the system. This, of course, also means that the mixed volume is a function of the particular substitutions used in converting the laminate equations to polynomial form (Equation 1), and of the parameters chosen as targets. The structure of the laminate equations above are such that the mixed volumes are all powers of two. The number of monomials present in the ABD stiffness target equations tends to result in a larger mixed volume than is present in a lamination parameter target system with the same number of equations.
Since the exact value of the mixed volume does depend on the laminate equation targets chosen, it is difficult to generalise about what to expect. However, it is possible to determine upper and lower bounds. In Table 3 the upper and lower bounds on mixed volume for a lamination parameter target system are given as exponents of two. For example, a full twelve lamination parameter target system will have a mixed volume of 2 18 = 262144, meaning that at most 2 18 stacking sequences may be found for such a system. In practice, the number of real finite solutions will be a small fraction of this number, but the mixed volume still provides a much smaller space in which to search for these real solutions than does the total degree of a twelve lamination parameter target system, which is closer to one billion. See Section 3.3 below for a numerical example. The change in behaviour between six and seven targets occurs because it is possible to construct systems without t k terms if there are six or fewer targets (and hence six or fewer unknown ply angles), meaning the coupling equation in 2 is not required. This same behaviour appears in the bounds for ABD stiffness target-based systems, given in Table 4.

System solution procedure
The solution procedure for recovering stacking sequences using PHC methods is summarised in Figure 1. The details associated with the start system creation and homotopy continuation steps have been largely omitted from this paper, and the reader is referred to a source such as [26] for more information. However, many of the freely available PHC software packages [36,37,38] allow a user to solve large systems of polynomial equations without a detailed knowledge of their workings. Note that a change in the numerical values of the targets changes the coefficients of the target system without altering its structure, allowing the same start system to be reused. Numerical examples illustrating the solution procedure as applied to a number of different types of laminates and targets are given in the following section. A link to access Matlab scripts that can generate PHCpack input files for each of the numerical examples is provided at the end of the paper. As Figure 1 indicates, laminate properties other than ply angles can be used as unknowns, however, as described in Section 2, this paper only explores the use of ply angles as unknowns in detail. thicknesses, material properties, . . .).
Create start system based on above structure using using polyhedral, (or other) methods Assign numerical values to targets Generate target polynomial system Coefficient homotopy from start to target system Assess range of real solutions against secondary considerations Re-run with new target values Check target system will not be singular otherwise. Figure 1: System solution procedure using PHC

Numerical examples
To illustrate the mechanics of recovering a stacking sequence from a set of lamination parameters or ABD stiffness targets, a number of numerical examples are presented in Sections 3.1-3.5 below. The examples are intended to show different aspects of the method's capabilities, rather than to be representative of any particular engineering problem. In the following numerical examples it is assumed that Ply 1 has the most negative z value, while Ply N has the most positive z value. Stacking sequences will be given in the order [1, 2, . . . , N − 1, N ] T . The stacking sequence angles will be given in degrees to one decimal place. This is not done in any expectation that such fibre placement accuracy is achievable (at least not easily) in practice, but rather to allow numerically similar but geometrically distinct solutions to be distinguished.

Example a -Unsymmetric laminate with four lamination parameter targets
As an initial example, consider a case in which all possible stacking sequences are to be recovered for a laminate in which the four V B terms have been specified, possibly as the result of an earlier optimisation process. The question here is how many laminates with the specific parameters V 1B = −0.15, V 2B = 0.4, V 3B = 0.2 and V 4B = −0.7 exist, neglecting any discrepancy between the laminate mid-plane and its neutral axis? Looking back to Equation 6, it is clear that this system will contain both t k and u k terms. Because of this, it will require additional coupling equations of the form of Equation 2 to link the two sets of variables. With four lamination parameter targets, a laminate consisting of four plies is needed to make the system fully determined. All four of the targets in this example are V B lamination parameters, and the four target system equations they generate can all be taken from the second equation in 6. Each of the four {t k , u k } pairs require their own coupling equation. Because all the plies are taken to be the same thickness, the ply thickness can be taken to be unity without loss of generality, with a laminate thickness of N , and the ply through-thickness locations h k can be calculated as a fraction of the overall thickness using the ply number k. The lamination parameters can then be expressed purely in terms of the unknown ply angles. The resulting system of eight equations in eight unknowns (after a small amount of re-scaling) is given in Equation 10: The first equation is clearly linear (degree 1), the second and third equations are degree 2, with the remainder of the equations being degree 3, bringing the total degree of the system to 2 2 .3 5 = 972. However, this system has a mixed volume of only 64 (= 2 6 ), meaning that there are at most only 64 finite isolated solutions to be found, depending on the particular monomial coefficients. By creating a random complex coefficient start system using the methods described in [40] only 64 solutions must be tracked from the start system to the target system. Following this approach, it is found that the system in Equation 10 has 14 complex solutions (appearing in 7 conjugate pairs), a further 46 solutions that diverge to infinity, and 4 real regular solutions. The problem is summarised in Table 5. The four real-valued solutions are given, along with the corresponding stacking sequences that they represent. In specifying only a subset of the twelve lamination parameters as targets, all control over the others is relinquished. For example, the full set of lamination parameters for the first two stacking sequence solutions in Table 5 are: Clearly all the V B parameters match, but the others vary greatly. When selecting which of the real stacking sequence solutions is the most appropriate for the application at hand, an engineer could apply a set of secondary selection criteria.     Note that it is difficult to predict a priori how many of the roots of a multivariate polynomial system will be real (although some methods of providing a bound do exist [44]). For example, another system of four targets: V 3A = −0.7, V 1B = −0.15, V 3B = 0.2, V 2D = −0.05, has a total degree of 648, a mixed volume of 64, and 12 real solutions.

Example b -Symmetric laminate with eight lamination parameter targets
In symmetric laminates, the same ply angle variables are used for opposing pairs of plies equidistant from the laminate mid-plane. A maximum of eight lamination parameters may be specified as targets for symmetric laminates, none of which may be coupling (V B ) parameters. As an example, find all possible stacking sequences for a laminate with in-plane lamination parameters V 1A = 0.04, V 2A = 0.14, V 3A = −0.45 and V 4A = 0.17, and out-of-plane parameters V 1D = 0.16, V 2D = 0.53, V 3D = −0.64 and V 4D = 0.01. The details of this eight-target system example are given in Table 6. For brevity the target equation system is not given here, and the list of real solutions has been omitted from the table.   The efficiency of using polyhedral methods in constructing start systems is on display in this example, with the mixed volume being only 0.4% of the total degree of the system. Rather than tracking nearly a million solutions to the target system, only 2 12 (= 4096) need to be checked. From this restricted set of possible solutions, only 10 are found to be real.

Example c -Unsymmetric laminate with twelve lamination parameter targets
In a sense, this is the heaviest lifting the PHC method can ever be expected to do in solving lamination parameter target problems. The full set of twelve lamination parameter targets requires a system of 24 equations (counting the coupling equations) in 24 unknowns to be solved. While the total degree of this system is an impressive 918330048 (= 2 6 .3 15 ), the mixed volume is a large, but more manageable 262144 (or 2 18 ). The full process of solving a twelve-target system using PHCpack, including creating a random coefficient start system and running the coefficient homotopy, took 2 h 15 min running 16 threads on an Intel Xeon W-2145 processor. The system and the results of the continuation process are summarised in Table 7.  The combination of lamination parameter targets given in Table 7 happens to generate 112 real solutions -too many to list in the table -so only ten are listed. Again, polyhedral methods are shown to be very efficient in in constructing start systems. The mixed volume being only 0.03% of the total degree of the system. A total degree homotopy would have taken just under one year to complete on the same processor.

Example d -Antisymmetric laminate with six ABD stiffness targets
Before constructing target systems based on the ABD stiffness terms for a laminate, it is necessary to assign some material properties to the laminae being used. Using the arbitrarily chosen but not atypical (for fibre reinforced polymers) values in Table 8, this example will illustrate the recovery of antisymmetric stacking sequences using six ABD stiffness targets. As was the case in the lamination parameter targets examples above, each target ABD stiffness target requires an extra equation of the form found in 4 to be added to the system. The ABD target equations are constructed using the expressions in 9, using theQ expressions from 8. The invariant stiffness values in 8 for the lamina used in this example are also given in Table 8. Six is the maximum number of targets that can be specified for an antisymmetric laminate, limiting the number of independent ply angles to six in the interest of maintaining a square system. Choosing the in-plane targets A 11 = 140 kN/mm and A 22 = 400 kN/mm, the coupling targets B 16 = −55 kN and B 26 = −110 kN, and the out-of plane targets D 12 = 250 kN.mm and D 22 = 1250 kN.mm leads to a target system in twelve equations and twelve unknowns, with a total degree of 2 4 .3 8 = 104976 and a mixed volume of 2 12 = 4096. The search for real stacking sequences is again quite fruitful, with twelve being found. The results are given in Table 9. Variables: {t1, . . . , t6, u1, . . . , u6} Solutions 3.5. Example e -Symmetric laminate with one buckling and three ABD stiffness targets As a further example, consider the case of a symmetric laminate for which a designer wants to assign the two perpendicular in-plane stiffnesses (A 11 and A 22 ) and the twist term (D 66 ) to specific values, but also wants the onset of firstmode buckling in a simply-supported version of the laminated plate to occur at a specific load. From [45], a plate of length a and width b, simply supported on all sides, has an axial buckling load given by: where N x is an edge loading in N/mm and m is the mode number along the loaded axis. It is assumed that the mode number in the direction perpendicular to the loaded axis is one. This example will use the ply properties in Table 8 again, and the plate geometric properties in Table 10.  where the D ij terms are functions of u k constructed using Equation 9 as usual. The task of ensuring a non-singular target system is slightly more complicated when constructing equations from linear combinations of ABD stiffness targets. At least one of the constituent targets in the combination equation must be independent of all the other equations in the system. The allowable combinations of targets in Table 2 may still be used to guide this checking process.
The details of this problem are shown in Table 11. Note that unlike the previous examples, this system can be constructed entirely in terms of the u k variables, leading to a relatively small total degree and mixed volume. It also means that the solutions do not carry any information about the signs of the ply angles. It can also lead to situations in which real solutions do not have a corresponding physically meaningful layup. In the current example, two of the four finite solutions are real, but the negative sign in one of the real solutions leads to a complex ply angle result. The two solutions themselves are shown in Table 11 to illustrate this, with only the second solution converted to a layup: [66.0/39.8/74.3/18.9] S . The ply angles are given without sign, but [-66.0/39.8/74.3/18.9] S and any other sign permutation would also satisfy the requirements. As always, the number of real feasible solutions varies with the value and combination of the targets and ply properties. The combination of targets A 11 = 200 kN/mm, A 22 = 150 kN/mm, D 66 = 100 kN.mm and N x = −1050 N/mm has four real solutions, of which three produce physically meaningful stacking sequences.

Conclusion
It has been shown that, assuming a feasible solution is known to exist, the PHC method can deterministically find all stacking sequences capable of producing a laminate with a specific set of stiffness properties. The particular approach presented here is limited in that it requires the system to be square, and has only been demonstrated with unknowns representing ply angles. While only the ply angles were used as design variables in this paper, the same continuation-based approach could easily be adapted to use other parameters such as ply thickness, volume fraction, or even the constituent component's material properties as unknowns in the stacking sequence recovery process. The target properties could also be expanded from lamination parameters, plate stiffnesses, and basic linear buckling criteria, to include thermal expansion, simple laminate strength requirements, or even bistable behaviour using models based on out-of-plane stiffness properties [46].
Another limitation of the method is that it does not provide any insight into how distant one solution may be from a more suitable stacking sequence obtainable with only a small change in the numerical value of the target stiffness values. Similarly, an engineer may need to assess the sensitivity of the laminate target properties to any small alterations in the ply angles produced by the PHC process necessary to round them to the nearest whole degree, or the nearest achievable discrete value using the manufacturing equipment available.
It is possible that an improvement to the efficiency of the method presented may be achieved through a combination of elimination and continuation techniques. Many of the most successful applications of polynomial-based methods to kinematic systems have involved first reducing the number of variables through a process of elimination, only switching to continuation when further elimination becomes impractical. The V 1x , for example, are linear in the u k variables, making the elimination process straightforward for systems containing these equations, while further elimination could be achieved via successive application of the method of resultants. It may even be possible to reduce some systems to univariate polynomials, for which highly efficient root finding methods are available. However, the degree of such a univariate polynomial may be very large, with the number of finite solutions to the twelve lamination parameter target system in Section 3.3 above suggesting that the degree of a univariate polynomial reduction of this system could be 60000 or more.
It may be that there is some utility to studying non-square systems, in particular those with more unknowns than equations. The traditional approach in such situations is to apply an optimisation technique to obtain the closest approximation to the targets using the available variables, but it is also possible to use continuation methods to capture the behaviour of these systems by using one or more of the extra variables as natural continuation parameters.
The solution process of a twelve-target lamination parameter system with a mixed volume of 2 18 was found to take a matter of hours to solve using a modern multi-core processor (see the example in Section 3.3 above). A similar process for a twelve-target ABD stiffness system with a mixed volume of 2 24 may take substantially longer (on the order of a week with the same processor), but it is certainly not beyond the realms of possibility that a practical attempt could be made at solving it with the aid of more parallel processing.
The number of unique stacking sequences capable of meeting a specific set of laminate stiffness requirements may exceed 100. By using a continuationbased approach to recover these stacking sequences, rather than optimisation or Newton's method which relies on the prudent selection of starting guesses, a designer can be sure of having every possible option to choose from.

Data availability
The Matlab files used to generate the PHCpack input files for the numerical examples above are accessible at: 10.15126/surreydata.9943043.