A compliance approximation method applied to variable stiffness composite optimisation

A way to approximate the compliance of composites for optimisation is described. A two-level approximation scheme is proposed inspired by traditional approximation concepts such as force approximations and convex linearisation. In level one, an approximation in terms of the reciprocal in-plane stiffness matrix is made. In level two, either the lamination parameters, or the nodal fibre angle distribution are used as design variables. A quadratic approximation is used to build the approximations in terms of the fibre angles. The method of conservative, convex separable approximations is used for the optimisation. Conservativeness is guaranteed by adding a convex damping function to the approximations. Two numerical examples, one optimisng the compliance of a plate clamped on the left, loaded downwards on the bottom right, another one optimising the compliance of a plate loaded with a shear force and a moment show the computational efficiency of the proposed optimisation algorithm.


Introduction
Composite materials are attractive due to their high stiffness-to-weight and strength-to-weight ratio. Traditionally, fibres within a layer have the same orientation, leading to constant stiffness properties. As manufacturing technology has evolved, for example the advent of automated fibre placement machines, the fibre orientation of a layer can be varied continuously leading to varying stiffness properties that can be best tailored for the applied loads. These composites are called variable stiffness laminates (VSL).
To develop constant thickness, steered laminates, a threestep optimisation approach has been developed (Ijsselmuiden 2011; Ijsselmuiden et al. 2010). In the first step the optimal nodal stiffness distribution in terms of lamination parameters is found, in the second step the optimal fibre angles at the nodes are obtained, and in the third step the optimal fibre paths are retrieved. These steps are done one after the other and are not repeated. The critical numerical part of this optimisation are the approximations used during the optimisation.
One way of approximating structural responses is using response surfaces. This means a set of response surfaces is chosen, which are scaled during optimisation. This method has been shown to be accurate, but choosing the response surfaces is problem-dependant (Venter et al. 1998). To approximate eigenvalues, or buckling factors, the Rayleigh-Ritz method is often used. This method has been used to optimise (variable stiffness) composites for buckling, but the shape functions need to be chosen carefully (Wu et al. 2012).
One structural response which has received a lot of attention over the years is the compliance of a structure. Even though in design, no structure is designed using compliance as objective or constraint, it has received a lot of attention over the years (Rozvany et al. 1994). When optimising composite structures, lamination parameters are often used as design variables due to the convex feasible region and limited number of design variables (Fukunag and Vanderplaats 1991). Directly optimising the compliance, discretised at design points is one example (Setoodeh et al. 2005). Disadvantage of using lamination parameters is that no unique stacking sequence corresponds to a lamination parameter distribution.
An interesting approach is introduced by Hammer et al., who approximate the compliance in terms of the lamination parameters based on the mid-plane strain and curvature of the elements. The distribution of lamination parameters is consequently matched using three plies of which the fibre angle and thickness are determined (Hammer et al. 1997a). This way the lamination parameters can be matched exactly. Whilst this method proved that lamination parameters are a good candidate for optimisation, and fibre angles can be found that match it, the method has disadvantages as well. For example, the obtained result is not physically feasible: a material has a certain ply thickness, one cannot just scale the thickness. Furthermore, the change in fibre angle between different points is not taken into account, which means a fibre placement machine is unable to lay it down. Finally, the convergence plot is not monotonically decreasing, meaning the method is not globally convergent.
Another approach is the discrete material and optimisation method (DMO) (Stegmann and Lund 2005). This method determines the best material choice throughout the structure. The different materials are defined by a 6 × 6 stiffness matrix that can represent a different fibre angle or completely different material. The method optimises different patches, that may consist of multiple elements in the FEA. The choice of the size and location of the patches has a significant influence on the final outcome and computational cost. Each patch can have a different material and different orientation. To make sure the outcome is manufacturable, extra constraints have been posed limiting the rate of change over the structure (Sørensen and Lund 2013).
A method that is related to DMO is proposed by Kennedy and Martins, where ply-selection variables are used (Kennedy and Martins 2013). Furthermore, the design variables are continuous, with additional constraints to obtain discrete values at the optimum. These additional constraints are not adduced as constraints directly, but using an l 1 penalty function, which leads to the same solution provided that the penalty parameter is sufficiently large (Kennedy and Martins 2013). Additionally, manufacturing constraints were added limiting the fibre angle difference between adjacent layers, but not between adjacent design regions, meaning manufacturability within a ply is not guaranteed.
A method that tries to combine analytical and numerical optimisation techniques is the DCOC approach Rozvany 1992, 1993). It is a discrete optimality criterion (DOC), which uses ideas from continuous optimality criteria (COC), hence the abbreviation. The basis of the technique is a dual formulation, but the efficiency is increased by computing the Lagrangian multiplier of the stress constraints explicitly. This reduces the computational cost of optimisation compared to using a dual approach significantly, whilst the cost of an analysis is the same. Whilst more iterations are necessary, the total computational time is still decreasing, and the optimum found is the same (Zhou and Rozvany 1992).
The approximate problem that is solved using the DCOC method in each iteration is formulated such that it decreases the computational effort. One of the features is a constraint selection procedure: only the constraints that are 'close' (the exact value is user-defined) to being active are taken into account during optimisation. Furthermore, the internal forces and displacements are assumed to be invariants in each iteration. Finally, a move limit is defined to reduce the chance that a constraint that was not selected becomes active (Zhou and Rozvany 1993).
Another method that is not formulated for a specific type of response is the mutual energy formulation (Taylor and Bendsøe 2001). In this formulation, a general objective/constraint is approximated using the stiffness and strain distribution of a structure. Both the displacement and adjoint displacement are used in the formulation. It is shown that for compliance optimisation, the general formulation reduces to the general, well known, min-max formulation. However, the paper does not provide any numerical examples, only theoretical formulations (Taylor and Bendsøe 2001).
A method to rewrite the min-max problems to a pure minimisation problem by posing additional constraints had been formulated earlier by the same authors (Taylor and Bendsøe 1984). Furthermore, some constraints are relaxed during the optimisation, leading to the disappearance of some singularities. Hence, by adding and relaxing constraints, the optimisation problem becomes easier to solve, although the number of constraints is increasing. It is shown on relatively easy examples that the methodology works (Taylor and Bendsøe 1984).
A method that is applicable in general is linearisation, using for example a Taylor series. Since stress and displacement are exact linear functions of the reciprocal sizing variables in case of a statically determinate structure, they are often approximated in a reciprocal way. A generalisation of linear and reciprocal approximations is the 'convex linearisation' (ConLin) method introduced by Fleury (1989). Whether a variable is approximated in a linear or reciprocal way is dependant of the sign of the derivative at the approximating point: a positive derivative is linearly approximated, a negative derivative in a reciprocal way.
Vanderplaats recognised that approximating for example frequency or stress constraints using the force in a member captured the non-linearity better, which led to better approximations and a more efficient optimisation. This led to a two-level approximation: first, the stress is approximated in terms of the section properties (e.g., member force), then, the physical properties are used to optimise the section properties. Once the optimal section properties are found, the stress is calculated based on the new properties, and a new approximation in terms of the physical properties is made. This is repeated until the stress converges (Vanderplaats and Thomas 1993;Vanderplaats and Kodiyalam 1990).
In this paper the idea of force approximations is used. In level one approximations the structural responses are approximated in terms of the in-plane stiffness matrix. For level two approximation, a quadratic approximation of approximation is constructed along the lines of the Gauss-Newton method. The two-level approximation should not be confused with the three-step optimisation approach: in this paper, level is used to describe the approximation, step is used in the context of optimisation.
First the derivation is done for a truss made of an elastic material, next, the approximations for a general twodimensional structure are derived. This is done since the physical reasoning behind the approximations is clear when a truss made of elastic material is used. For the general two-dimensional plates the complexity of the equations may hide their physical meaning. The approximations will be developed for compliance only in this paper. Before diving into the approximation strategy, the method of successive approximations, the requirements for the approximations and the definition of the design variables are discussed in Section 2. The derivation for the compliance of a truss is done in Section 3. Next, the compliance approximation for a general composite plate is done in Section 4. Next, the approximations in terms of the lamination parameters, and fibre angles are discussed in Sections 5 and 6 respectively. For all approximations, it will be proven that they have the desired properties. The optimisation strategy is discussed in Section 7, and a numerical example is worked out for the compliance approximation in Section 8, two optimisation problems are solved in Section 9. This paper is concluded in Section 10.

Method of successive approximations
The method of successive approximations replaces the optimisation of the problem by a sequence of approximate sub-problems. The first approximate sub-problem is built at a user-defined point. The requirements for the approximation are discussed in Section 2.1. This approximate sub-problem is optimised to find the next iterate. The process continues by building an approximate sub-problem at the new iterate and optimising until convergence is reached (Schmit and Fleury 1980;Bruyneel et al. 2002). A flowchart of the method of successive approximations can be found in Fig. 1.
where f 1 to f n denote structural responses that are optimised and f n+1 up to f m denote structural responses that are constrained. These reponses are all functions of the design variables, denoted by x. The problem is defined as a minimisation problem. Hence if we want, for example, to maximise the buckling load, the inverse buckling load is minimised. Another example is maximising the stiffness, which is formulated as minimum compliance. As final example, the factor of safety is not maximised, instead the failure index is minimised.
Defining the objective as worst case is useful when for example performing a buckling optimisation: by taking multiple modes into account, mode jumping is not a problem (Seyranian et al. 1994). Another example is stress optimisation: the maximum failure index appearing in the structure should be minimised.
The convergence criterion used is a soft convergence criterion. If the improvement of the objective function is less than a certain tolerance and the constraints are satisfied, the optimisation is assumed to have converged. The tolerance is usually a function of the initial value of the objective: an improvement smaller than 10 −3 of the initial value is often used. The exact function is used to determine convergence, not the approximation.

Requirements of the approximation
An approximation has to have certain properties to be used during the method of successive approximations. Since the method is gradient-based, a first-order approximation is used, meaning (Svanberg 1987): where f denotes the exact function,f the approximation, and x 0 the approximation point. For optimisation purposes, four more properties for approximation are favourable (Boyd and Vandenberghe 2004) -convex: if the approximation is convex, it is guaranteed to have a solution, thus optimising the sub-problem will always give a solution when starting from a feasible point. Mathematically, a function f is convex if for any two points x 1 and x 2 in the feasible domain it holds that where t is any value between 0 and 1. Another sufficient condition for convexity is if the second derivative of a function is positive: -separable: for problems with large number of design variables, like in problems addressed in this work, a separable approximation is desirable. This means that the different design variables do not influence each other. This makes the optimisation computationally efficient. Mathematically, a function is separable if it can be written as a summation of functions of single variables: In this work, separable is interpreted slightly different: x i does not need to be a scalar variable, it can be a (small) vector or a tensor. -conservative: an approximation is conservative if, at each point in the feasible domain, the function that is approximated is lower or equal to the approximation. Mathematically, for a minimisation problem this means where D denotes the feasible domain of f (Boyd and Vandenberghe 2004). f andf denote the exact and approximate function respectively. As we shall see in Section 7, conservativeness plays an important role in guaranteeing global convergence of the total optimisation problem.
-homogeneous: an approximation is homogeneous if the response scales with a certain factor when all design variables are scaled. Mathematically, a function is homogeneous of degree n if This implies a solution can always be found, even if the starting point is infeasible, given that the upper and lower bounds on the design variables allow the required scaling.
These four properties are advantageous for the optimisation, but only convexity is required to use the method of successive approximations. The approximations used in this approach are discussed in the following sections. The approximations themselves are convex, separable, and, if possible, homogeneous and conservative. To render them conservative, an extra part, called damping function in this work, is added to the approximation. This is discussed in Section 7.

Definition of design variables
Before an approximation can be made, the location of the design variables has to be defined. Since the structural responses are calculated using a finite element analysis (FEA), it seems a logical choice to link the design variables to the elements. However, rather than at the elements, the design variables are defined at the nodes of a FE model. This has three advantages. One, the continuity of the design variables is more likely to happen; it cannot be mathematically guaranteed, but in all numerical results presented continuity was preserved. Two, the number of design variables is reduced when triangular elements are used. This in general leads to more elements than nodes. Three, the manufacturing constraint on the minimal turning radius can easily be defined, as is shown in other work by the authors (Peeters et al. 2015). The minimal turning radius is the radius necessary to lay down the fibre when manufacturing using an automated fibre placement. When the design variables are defined in the element, they are next to each other and it is hard to determine the turning radius. When the design variables are at the nodes, the required turning radius can be determined based on the distance between them and the magnitude of the change in angle.
To calculate the in-plane stiffness matrix of an element A e , the stiffness matrices at the nearby nodes A n are interpolated reciprocally. For instance, in triangular element A e is calculated according to For the general case, the stiffness matrices at the different Gauss points can be found using where G denotes the number of nodes that have influence on this Gauss point and N ng denotes the shape function of node n at Gauss point g. The same equations hold for the out-of-plane stiffness matrix D.
Once the A-and D-matrix of the elements are known, the stiffness matrix of the entire model K can be generated. From this point onwards, the standard FEM can be applied. The displacement field u is found using (10) After the displacements and rotations of each node are found, the strain and stress at each Gauss point can be recovered. For a more detailed description of a FEM, the reader is referred to, for example, (Felippa 2001).

Compliance approximation of a truss made of elastic material
First, the compliance of a truss made of an elastic material will be derived to illustrate the procedure. An example of a truss can be seen in Fig. 2. The strain energy of a system is defined as where E denotes the Young's modulus, l the length and the subscript e denotes the element. The cross-sectional area of each element A e are the design variables in this case. The principle of minimum total potential energy is formulated as The Lagrangian for this problem can be written as where σ e are the Lagrangian multipliers of the constraints. The optimum can only be reached if the constraint is satisfied; if the constraint is not satisfied, the Lagrange multiplier will go to either +∞ or −∞, and thus the Rewriting to combine the terms involving the displacement vector u and interchanging the min and max, which is allowed in this case, (13) is rewritten to Observing that u acts as a Lagrangian multiplier of an equality constraint, this can be rewritten to the following equivalent problem, which is the principle of minimum total complementary energy: Using where F e is the internal force of element. Implementing the correct expression for f c , the complementary strain energy is found to be Based on the conservation of energy, the strain energy of the truss U equals the work done by the external force, which equals the compliance of the structure by definition.
For structures of linear elastic continua, strain energy U equals complementary energy U * .
From the principle of total complementary energy, the complementary energy of the structure which satisfies both static and kinematic compatibility, is minimised. Therefore the compliance of the truss is By assuming the internal force F e is constant when optimising the design variables, the compliance can be approximated in terms of the cross-sectional areas as where the superscript (k) denotes the force after the kth iteration, when the current iteration is k + 1. The symbol I f is used since it is a level one approximation.
Therefore, the minimal compliance can be formulated as Observing (20), it is noticed that the compliance is reciprocal in terms of the cross-sectional areas.
The four desirable properties mentioned in Section 2.1 are satisfied: -convex: The second derivative of the compliance with respect to any design variable A e is, Since all terms are strictly positive, the second derivative with respect to any design variable is positive, and thus the approximation is convex. -separable: the approximation is a summation of different functions of A e .
-conservative: as has been shown in (21), the compliance C of a truss with cross-sectional area A is found by minimising the complementary energy: F * , which minimises the complementary energy, is both statically and kinematically admissible. However, the internal force F (k) at the approximation point only satisfies the equilibrium condition in the new loop since this is defined as a constraint. Hence, it is not guaranteed to be kinematically admissible. Therefore, according to the principle of total complementary energy, the compliance obtained from (20) is an upper bound for the exact compliance from (23): Hence, the exact compliance is always lower than or equal to the approximation, meaning the approximation is conservative. -homogeneous: the approximation is homogeneous of order −1.

Compliance approximation of a general two-dimensional composite plate
For a general two-dimensional plate with area , and boundary , as shown in Fig. 3, the compliance equals the minimum of the complementary energy of the structure, which follows the same logic as the truss.
Since the stiffness A n is estimated at each node from the design variables, using the material law (i.e., N ij = A ij kl · kl , Fig. 3 Example of a two-dimensional plate for i, j, k, l = 1, 2), the density of the complementary energy can be rewritten as where N n is the resultant force at the node. Using the Frobenius product, the approximation of the compliance is rewritten to where N denotes the number of nodes in the finite element model of the plate, and defining φ n as where A n denotes the area of node n. However, this does not work well: the forces at the nodes are a function of multiple elements. Hence, finding an appropriate expression directly in the form of (20) is not straightforward.
Therefore, we need to go back to the continuous model, (26), and follow a similar procedure as described in Section 3 to obtain a discretised version.
The strain energy of the plate is defined as where is the strain vector of the structure. Discretising the plate to be used in the finite element method, the strain energy of the plate can be expressed as a summation at every Gauss point using the Gauss integration scheme. Thus the strain energy can be obtained, where w g is the weight coefficient times the determinant of Jacobian matrix of Gauss point g. The subscript g denotes the variables at the Gauss point. The total potential energy of the plate is The principle of total potential energy leads to min g ,u s.t. g − B g u = 0, where B g is the strain-displacement matrix at the Gauss point. The Lagrangian is found to be where λ is the Lagrangian multiplier. The optimality condition with respect to g gives λ g = −w g A g g . Hence, where N g is the stress resultant vector at Gauss point g. Substituting (32) into the Lagrangian (31), and replacing g with A −1 g N g , (31) is rewritten as Inverting the order of min and max and rearranging terms leads to Thus the equivalent optimisation problem can be written as where the objective of this optimisation problem is the complementary energy of the structure U * : Using (9), the complementary energy of the plate is By changing the order of summation and employing the Frobenius product, this is rewritten to Assuming the resultant force at each Gauss point N g remains constant when the stiffness A n is changing during the optimisation, according to (19), the approximation of compliance of a plate in discrete form is D. Peeters et al.
The expression for φ n is where N (k) g is the internal force that is both statically and kinematically admissible at the kth approximation point, when this is iteration k + 1. Comparing this approximation to the one found in the beginning, in (26), it can be seen that the form is exactly the same.
The minimisation of the compliance is formulated as where A is the constitutive matrix of each node. The four desirable properties mentioned in Section 2.1 are satisfied: -convex: The approximation of the compliance can be rewritten as where λ l (φ n ) is the eigenvalue of the matrix φ n , and v l is the corresponding eigenvector. Since φ n is positive definite by construction, λ l (φ n ) is positive. Geometrically, (42) describes ellipsoids. Since the constitutive matrix, which is the design variable, is guaranteed to be symmetric and positive definite, the formulation is convex according to the theory of convex optimisation (Boyd and Vandenberghe 2004). -separable: observing (39), it is noted that the constitutive matrices A n are the design variables, which do not influence each other, hence the approximation is separable. -conservative: similar to the truss, the internal force N g that minimises the complementary energy is both statically and kinematically admissible. However, the internal force N (k) g in (39) is only statically admissible, and not necessarily kinematically admissible in the (k + 1)th optimisation loop. Therefore, the approximation is equal or larger than the exact compliance in (41). As a result, (39) is conservative.
-homogeneous: the approximation is homogeneous of order −1.

Approximations in terms of the lamination parameters
Whilst only one level one approximation exists, for level two, multiple options exist. In the second level approximation, the stiffness matrices are approximated in terms of physical design variables. During step one of the three-step optimisation approach, the level two approximation is in terms of the lamination parameters. The three-step approach has been briefly discussed during the introduction, it was originally proposed in the work of Ijsselmuiden (2011). Approximating the stiffness matrix in terms of the lamination parameters is not strictly speaking an approximation since the lamination parameters describe the stiffness matrices exact. During step two, physical design variables exist: fibre angles. The approximation in terms of the fibre angles is discussed in the following section. This section only focuses on the approximation in terms of the lamination parameters.
The definition of the lamination parameters is given in Appendix A. Using lamination parameters, the expressions for the A-and D-matrix simplify considerably to where the laminate stiffness matrices are found as functions of the lamination parameters (LPs) and laminate thickness. The stiffness matrices can be rewritten as whereĀ andD are the normalised stiffness matrices. During step one of the three-step optimisation approach, the stiffness is optimised. The terms of the stiffness matrices are linked, hence directly optimising them is not easy. One would have to use a lot of constraints to assure feasibility. By using the lamination parameters, the feasible region can easily be described. Since the lamination parameters describe the feasible region exact, the level two approximation is an explicit function, much like Vanderplaats proposed for the strength approximation.
The form of the approximation of compliance, introduced in (39), hence changes to where I I f is used to denote it is a level two approximation. The definition of φ n does not change, it is still given by (40).
The advantage of using lamination parameters as parametrisation is that, independent of the number of layers, five design variables are used: 4 in-plane LPs, and the laminate thickness. For optimisation of a constant stiffness laminate, one set of in-plane LPs, and a thickness is sufficient. If variable stiffness, or variable thickness laminates are optimised, multiple points across the structure will have a set of LPs. The feasible region only considers the feasibility of a single laminate, not whether the change from one set of LPs at one point to the set at an adjacent point is manufacturable. Disadvantage of using LPs is that the lay-up of the laminate is unknown: a set of LPs does not describe a unique stacking sequence.
Furthermore, the feasible region is convex. A convex feasible region implies the optimum found during step one is the global optimum. The details of the lamination parameter optimisation are not discussed here, the interested reader is referred to the PhD thesis of Ijsselmuiden (2011).
When changing the design variables, the four desirable properties mentioned in Section 2.1 are not necessarily preserved. Hence, they are checked again: -convex: in terms of h andĀ, convexity can be proven relatively straightforward. The second derivative of the compliance with respect to h is which is guaranteed to be positive: the compliance is always positive and the thickness as well. Hence the approximation is convex in terms of the laminate thickness h. In terms ofĀ, the same proof as in the previous section can be used, when definingφ = φ h . Since h is a positive number,φ is positive semi-definite, just as φ, thus the same proof of the previous section can be given usingφ andĀ instead of φ and A. Whether the combined approximation holds as well still needs to be proven. This is done in Appendix B.
-separable: as was explained in Section 2, strictly speaking this function is not separable since the laminate thickness and different lamination parameters influence each other, but since it is still separable per node, meaning per five design variables, the approximation is still regarded as separable. -conservative: since the definition of φ n does not change, the proof of conservativeness is still the same as in the previous section. -homogeneous: the approximation is homogeneous of order −2.

Approximations in terms of the fibre angles
During the second step of the three-step optimisation approach, the fibre angles are optimised. This is done by building a level two approximation, denoted by I I f (θ ). Contrary to the lamination parameters, the fibre angles only represent the stiffness matrices exactly at the approximation point. The approximation is a second-order Taylor series, based on (39): where I f 0 denotes the value, g the gradient and H is an approximation of the Hessian of the first level approximation at the approximation point. The gradient and Hessian approximation can be calculated starting from Deriving this with respect to the fibre angle θ i , the ith term of the gradient is found to be Deriving again with respect to fibre angle θ j , the ij th term of the Hessian is found to be Using the exact Hessian, convexity is not guaranteed. Convexity is ensured by omitting the underlined part of (50), which is not guaranteed to be positive definite, and leaving the positive semi-definite leading term, called the Gauss-Newton part. A first-order approximation only has to have equal function and gradient values at the approximation point as the approximated function. Thus, using only part of the Hessian does not influence the validity of the approximation. The derivative of the A-matrix with respect to θ k , the fibre angle of layer k, is The four desirable properties mentioned in Section 2.1 are checked: -convex: by construction, the Hessian is positive semidefinite, hence, the approximation is convex. -separable: since the approximation in terms of the stiffness was separable at the nodes, this approximation is also separable per node. -conservative: the approximation is not guaranteed to be conservative. How this can be guaranteed is explained in Section 7. -homogeneous: the approximation is not homogeneous.
Hence, only two out of four desirable properties are satisfied in this case. However, the function can be made conservative, as will be explained in Section 7, and homogeneity is not strictly necessary for an approximation to be used in the method of successive approximations.

Optimisation algorithm
The method of conservative convex separable approximations (CCSA) will be used (Svanberg 2002). This means we start from a starting point and use the approximations to find the next point.
As shown in the previous sections, the approximations, either in terms of the stiffness, or the fibre angles, are always convex and separable. The approximation in terms of the stiffness is also conservative, however, the approximation in terms of the fibre angles is not conservative. To render the approximation conservative, a damping function is added.
wheref is the approximation built in Section 6, d is called the damping shape and ρ the damping factor. The damping function chosen for the fibre angle optimisation is where θ is the change in angles from the approximation point of the level two approximation, and H d is a regularisation matrix given by with α given by where is a damping factor, usually chosen to be 1. This is done for the approximation in terms of the fibre angle only since the level one approximation is conservative without the damping function. For a more detailed description of the damping function, and the optimisation procedure, the reader is referred to earlier work by the authors (Peeters et al. 2015).
The solution procedure to determine the fibre angles which starts after step one of the three-step optimisation is shown in Fig. 4, and is explained in Algorithm 1. For manufacturing reasons, a steering constraint is posed to ensure the fibres do not wrinkle whilst being laid down.
The details of the steering constraint are omitted here, the details can be found in earlier work by the authors (Peeters et al. 2015).

Results compliance approximation
In this section numerical cases are solved with the proposed compliance approximation. A square plate with sides of 300 mm is loaded with a shear force of 675 N, divided in a quadratic manner, and an equivalent moment of 303.75 Nm on the left and 506.25 Nm on the right, which is applied as a force in x-direction. To avoid rigid body motion the plate is simply supported in the middle of the left and right edge. A graphical representation is shown in Fig. 5. The material used has the following properties: E 1 = 177 GPa, E 2 = 10.8 GPa, G 12 = 7.6 GPa, and ν 12 = 0.27. To cheque the approximations, the material of the plate is assumed to be quasi-isotropic (QI), which is defined as all lamination parameters equal to zero. The thickness of the plate is 10 mm.
In this section it is shown that the approximation is meshindependent, converges to the exact compliance and that (40) and (27) converge to the same result.

Mesh independence
First, the mesh independency of the sensitivity φ, which is defined in (40) is discussed. φ n is the summation of internal force N g , shape function N ng and w g of the relevant Gauss where A g is the area at the gth Gauss point connected with the nth node. Since the internal force and the shape functions from FEM analysis are mesh independent, the variablesφ n should be mesh independent. Since φ n is a 3×3 matrix, the mesh independency ofφ is illustrated by plotting one term of the matrix at each node, in this work,φ(1, 1) is used. The figures obtained using different mesh sizes are compared visually. Observing Figs. 6, 7 and 8, it is noticed thatφ(1, 1) is mesh independent. Similar phenomenon happens to the other terms ofφ, which is not shown in the paper for brevity.
Observing Figs. 6 to 8, the value ofφ(1, 1) corresponds to the top left corner in Fig. 5 changes merely from 1007 to 1094 as we double the mesh three times, with 7.9% change. Meanwhile, the value related to the top right corner changes only from 383.2 to 401.1, with 4.5% change. The configuration of theφ(1, 1) is the same in these figures. The distribution ofφ(1, 1) is therefore mesh independent.
Similar phenomenon happens to the other terms ofφ. Hence, the compliance approximation in (39) is mesh independent.

Convergence of sensitivity with respect to mesh
Next, the convergence ofφ with respect to the number of elements is checked. Since the distribution is already shown to be mesh independent, checking a single node is  Fig. 6φ(1, 1) distribution using a 20 × 20 mesh sufficient. For convenience, the node in the middle of the plate is chosen. The normalised sensitivity at this node can be decomposed: where the subscript c denotes the centre of the plate, λ i denotes the eigenvalues ofφ c in descending order, and e i is the corresponding eigenvector. From (40), φ c is positive definite and has rank one. Therefore λ 1 should converge to a positive value, whilst λ 2 and λ 3 converge to zero. Hence, λ 1 e 1 e T 1 can be used to cheque the convergence ofφ c . This essentially corresponds to N n · N T n in (27). In order to estimate the 'exact' N , denoted by N * , whilst limiting the computational work, Richardson's extrapolation is employed to evaluate λ * 1 e * 1 , which is the Fig. 7φ(1, 1) distribution using a 40 × 40 mesh Fig. 8φ(1, 1) distribution using a 60 × 60 mesh exact internal force N * . Here the exact internal force is interpreted as: where N 1 (h) denotes the initial value, which can be found using a rough mesh with mesh size h: Every time the mesh size h is refined in Richardson's extrapolation, the order of the error in (58) is increasing. In other words: N * is approximated better with every refinement, whilst taking the information from the rougher meshes into account. From N * , the exact sensitivity can be found usinḡ The relative error betweenφ h c and theφ * c from different meshes size is calculated using where h denotes the mesh size, and e h φ c is the relative error. A plot of this error on a logarithmic scale is shown in Fig. 9. Observing this figure, it is found that the error decreases with a slope of −1, implyingφ c converges linearly as a function of the number of elements. Thus the compliance approximation of a plate in discrete form, (39), has a linear convergence rate.

Convergence of sensitivity
Finally, it is checked whether the definition of φ n in (40) and (27) converge to the same result. The central node is again Fig. 9 Relative error ofφ c with respect to the number elements used to check this. The internal force at the central node is calculated from the summation of the internal force at the surrounding Gauss points weighted by the area correlated to each Gauss point where N cg denotes the internal force at Gauss point g around central node c. This is graphically shown in Fig. 10. As a reference, the internal force N * from Richardson's extrapolation is used. The relative error is calculated as where N h c is the internal force at the central node of the plate with mesh size h. The absolute value is used to   Fig. 11. Observing this plot, it can be seen that N c is linearly converging as a function of the number of elements. Hence, (40) and (27) converge to the same result. Consequently, compliance approximation derived from the continuous form, (39), converges to the one derived from the discrete form, (26).

Plate under point load
The first example is a flat composite plate, clamped on the left, with a vertical force downwards on the bottom right. A graphical representation of this problem is shown in Fig. 12. This problem has been previously solved (Nagy et al. 2010). The material data are as follows: E 11 = 148GP a, E 22 = 9.65GP a, G 12 = 4.55GP a and ν 12 = 0.3. In this examples, 6 design layers are used. Since the final laminate is to be balanced and symmetric, each design layer represents 4 layers in the actual laminate: a negative of the design layer is right next to the layer, and the complete stack is symmetric. The only objective in this case is the compliance, the laminate is assumed to be balanced and symmetric, hence V 2 and V 4 are zero. The dimensions are changed compared to the problem solved by Nagy et al.: a = 0.4m, b = 2m. This is done to make the steering constraint active during the retrieval step. The scaling has no influence on the optimal compliance since the aspect ratio is the same, and the compliance is normalised using the compliance of a quasi-isotropic plate, meaning all lamination parameters equal to zero.
When running the optimisation in terms of the lamination parameters, the obtained result is shown in Fig. 13. Only In-plane lamination parameters after fibre angle retrieval four iterations were necessary to obtain the optimum. The normalised compliance is 0.4300, which is even a bit lower than was found by Nagy et al. However, the difference is not large, so could be due to the finer mesh used in the current work. The distribution is similar to what was found in the paper, proving that the approximation algorithm works well, and that the optimisation is behaving as expected.
Studying Fig. 13 in more detail, it is noticed that in the region of the load introduction V 1 goes towards −1, implying a fibre angle close to 90 • . This is as expected: at the load introduction the fibres should be aligned with the load. In the part above the load introduction, V 1 is staying constant, whilst V 3 is changing, implying the angles change to 45 • , leading the load towards the clamped region. Close to the clamped edge, both V 1 and V 3 tend to 1, implying that the fibre angle goes towards 0 • in this region. This means that the load is led in a straight line from the right to the left, after being changed direction by fibre steering. Hence, this result makes sense, and we have a clear idea of what the fibre angles should look like.
A steering constraint of 400 mm is used during the retrieval. Not all 6 layers are shown, just two design layers are shown in Fig. 14. The normalised compliance after retrieval is 0.443. This already shows that the values found in terms of the fibre angles can match the result found in terms of lamination parameters fairly close. To highlight how close both cases are, the V 1 and V 3 distribution after fibre angle retrieval is shown in Fig. 15. Here it is noticed that the V 3 distribution matches quite good, but the V 1 distribution looks different. This is due to the steering constraint: the angles cannot change very abrupt, so the fast change around the load introduction point cannot be matched in terms of the fibre angles.

Plate under shear and moment
The same example as in Section 8 is used to optimise the compliance of the composite plate. The objective is to minimise the compliance, no constraints are posed. The laminate has 36 layers in total, since the laminate is designed to be balanced and symmetric, 9 design layers are to be optimised.
During step one of the three-step optimisation, the different lamination parameters are optimised. Since the laminate is to be balanced and symmetric, V 2 and V 4 are zero. The optimal compliance is found after only 4 iterations, as can be seen in the convergence plot in Fig. 16. The optimal compliance is normalised with respect to the compliance of the initial quasi-isotropic (QI) plate. The final result has a compliance of 0.5446. The distribution of V 1 and V 3 is shown in Fig. 17.  Studying the lamination parameters distribution more closely, it is seen that at the top and bottom both V 1 and V 3 tend to 1, meaning an angle of 0 deg which makes sense since the top and bottom mainly take the moment, and thus have to be stiff to move as little as possible due to the moment. In the middle, V 1 tends to zero, and V 3 towards −1, which means the angle are ±45 • meaning the shear is taken by these regions. In between the angles are a good compromise between both loads.
During step two of the three-step optimisation, the fibre angle distribution is optimised. A manufacturing constraint is posed: a lower bound on the steering radius of 333 mm Fig. 18 Convergence of the optimisation problem in terms of the fibre angle distribution is posed. This optimisation is done in just four iterations, as can be seen in Fig. 18. The final compliance found is 0.6842, hence relatively close to the theoretical optimum found in terms of the stiffness. The reason for the decrease in performance is the manufacturing constraint, limiting the rate of change in stiffness, which is high in terms of the lamination parameters, as can be seen in Fig. 17. The initial Fig. 19 Optimised in-plane lamination parameters point is defined as [±30 ± 70 ± 40 ± 20 ± 50 ± 65 ± 40 ± 10 ± 20] S , where the balanced angles are next to each other. The starting point has a performance close to one, which is due to the starting point, which is relatively close to quasi-isotropic behaviour.
For the optimised fibre angle distribution step three of the three-step optimisation, fibre path retrieval, is performed. Not all fibre paths are shown, two design layers are selected and shown in Fig. 19. In total, each half consist of 18 layers, leading to a balanced, symmetric laminate of 36 layers. On these plots, the two balanced layers, which are next to each other, are shown on top of each other. The angles are as expected: tending to 0 • on the top and bottom and towards ±45 • in the middle part. Since each layer is slightly different, and the steering constraint is active, the steep change in lamination parameters cannot be matched exactly, but in general the influence on the compliance is small.

Conclusion
An approximation approach for the compliance optimisation of composites is proposed. The three-step optimisation approach for variable stiffness composites is used as a guideline. In step one the optimal stiffness distribution is found, in step two the optimal nodal fibre angle distribution is found, in step three the fibre paths are found. The approximations in this paper focus on step one and two.
Analogous to the force approximation approach, a two-level approximation is proposed. In level one, the compliance is approximated in terms of the inverse in-plane stiffness matrix. In step one, a level two approximation in terms of the lamination parameters is used to optimise the stiffness distribution. In step two, a second-order Taylor approximation in terms of the change in fibre angles is made at the approximation point as a level two approximation.
The conservative, convex separable approximation method is used during the optimisation. Convexity and separability are guaranteed by construction of the approximations. Conservativeness is not guaranteed. A damping function and damping factor are added to the fibre angle approximation to guarantee conservativeness and thus global convergence.
Two example problems are solved: a plate clamped on the left with a load downwards on the bottom right, and the compliance of a flat plate loaded with a shear force and moment were optimised. In both cases the number of iterations necessary in step one and two was low: not more than five iterations were necessary. This means that the whole optimisation only requires about ten finite element evaluations, which is computationally very efficient. This shows that using these approximations good results can be obtained using limited computational resources.
The lamination parameters are defined as cos (4θ(z)) , sin (4θ(z))) dz (W 1 , W 2 , W 3 , W 4 ) = 1 2 − 1 2z 2 (cos (2θ(z)) , sin (2θ(z)) , cos (4θ(z)) , sin (4θ(z))) dz, . (71) The feasible region of the lamination parameters is defined as the region where a stacking sequence can be found that gives the combination of lamination parameters. From their definition in (71), the feasible region of the in-or out-of-plane lamination parameters separately can be found to be (Hammer et al. 1997b): For the combination of in-and out-of-plane lamination parameters, the feasible region does not have an easy definition. It can be found in for example Setoodeh et al. (2006), Bloomfield et al. (2009), or Wu et al. (2013. Observing (71), it can be seen that if the laminate is balanced, meaning for every θ there is a −θ , V 2 and V 4 are equal to zero. The out-of-plane LPs are generally all non-zero, also for balanced laminates.

Appendix B: Proof of convexity of lamination parameter approximation
In Section 5, it was already proven that the approximation in terms of the lamination parameters is convex in terms of the laminate thickness, and in terms of the normalised in-plane stiffness matrix separately. However, both can be changed at the same time, hence also in terms of both convexity has to be proven. To proof this, start from f (x, y) = g(x) · h(y), where g is convex in terms of x, and h is convex in terms of y. It has to be proven that f is convex in terms of the combination of x and y. It is assumed that g(x) and h(y) are both positive. Next, it is assumed g(x) is homogeneous of order m, and h(y) is homogeneous of order n, meaning g(λx) = λ m g (x) h(λy) = λ n h(y), where λ denotes a real number. Since g(x) and h(y) are convex, the second variation is positive: To keep the formulas short, define Finally, the total vector of design variables of the function f is defined as z = [x, y].
The function f (z) is convex if the second variation with respect to z is positive: This needs to be proven. The second variation of f with respect to z is given by This matrix has to be positive definite. Since h is positive and P is positive definite, the upper left part of the matrix clearly is positive definite. This implies that if the Schur complement is positive, the complete matrix is positive definite as is stated in Appendix A of Boyd and Vandenberghe (2004). The Schur complement is given by which is guaranteed to exist since h is positive and P is positive definite and hence invertible. Using Euler's theorem for homogeneous functions, it is known that Deriving both sides with respect to x leads to Substituting the expression of (75) and (76), and rearranging terms: Since P is invertible, we can multiply both sides from the left with a T P −1 : Using (80), this is rewritten as a T · P −1 · a = m m − 1 · g.
In an analogous way, it can be found that Substituting (84) into (79), and dividing by g which is allowed since g is positive, the Schur complement is rewritten as Since Q is positive definite, using the Cholesky decomposition, it is rewritten as Q = L · L T . Using the Cholesky decomposition, the Schur complement is rewritten as This is positive if the part between brackets is positive. This part can be rewritten as where e is normalised, hence λ is thus defined as Using (85), this can be rewritten as Observing (88), it is noticed this is positive if λ < 1. Hence, For this specific case, (45), the approximation is homogeneous of order −1 in terms of both the normalised in-plane stiffness matrix, and the thickness. Since both functions are positive, and conservative in terms of the normalised inplane stiffness matrix and thickness separately, all assumptions used in this derivation hold. Filling in m = n = −1 in (92), it is found that the approximation is indeed convex in terms of the combination of normalised in-plane stiffness matrix and thickness.