Convex optimization techniques in compliant assembly simulation

A special class of quadratic programming (QP) problems is considered in this paper. This class emerges in simulation of assembly of large-scale compliant parts, which involves the formulation and solution of contact problems. The considered QP problems can have up to 20,000 unknowns, the Hessian matrix is fully populated and ill-conditioned, while the matrix of constraints is sparse. Variation analysis and optimization of assembly process usually require massive computations of QP problems with slightly different input data. The following optimization methods are adapted to account for the particular features of the assembly problem: an interior point method, an active-set method, a Newton projection method, and a pivotal algorithm for the linear complementarity problems. Equivalent formulations of the QP problem are proposed with the intent of them being more amenable to the considered methods. The methods are tested and results are compared for a number of aircraft assembly simulation problems.


Introduction
In the last decade a new modeling approach has been developed and applied to variation simulation and assembly optimization in aerospace and automotive industry (see Lupuleac et al. 2010Lupuleac et al. , 2011Lupuleac et al. , 2019bDahlström and Lindkvist 2007;Lindau et al. 2016;Yang et al. 2016). This approach considers the contact interaction between compliant parts. By using the variational formulation (Galin 1961;Tu and Gazis 1964;Lions and Stampacchia 1967;Kinderlehrer and Stampacchia 1980) and the substructuring (Turner et al. 1956;Guyan 1965;Wriggers 2006;Petukhova et al. 2014), contact detection is reduced to a quadratic programming problem that allows for large-scale computations of contact problems during variation simulation and assembly optimization. Since the accuracy of the assembly model influences the simulation results, it is important to use fine FEM meshes that can reflect the most essential features of assemblies. Implementation of highly refined computational meshes and need for massive computations of problems with slightly different input data require choosing appropriately tailored solvers for the QP problem; that is the primary objective of the present paper.
The paper is organized as follows: in Sect. 2 the contact problem arising in assembly simulation is presented and its features are described; task-level parallelization is also discussed. Section 3 is dedicated to the adaptation of numerical methods for solving the described contact problem. In Sect. 4 computational results for the adapted methods are provided and discussed for several problems of aircraft assembly simulation. Conclusions are presented in Sect. 5.

Formulation of contact problem
A finite element model of several parts fixed in an assembly jig is schematically shown in Fig. 1. This kind of assembly is of sandwich type, i.e., at least three parts are lap-joint. The assemblies with no intermediate parts are called simple (for instance, the assemblies consisting of two parts). The area where the parts overlap is called the junction area. The finite element meshes of the parts are assumed to be conformal. The nodes of the finite element model that are located in the junction area are denoted as computational nodes. These nodes are highlighted in Fig. 1. The degrees of freedom in the computational nodes considered in the analysis are arranged into the vector x ∈ ℝ n , where n is the dimension of the problem. Although where f ∈ ℝ n is the vector of loads (for example, loads from fastening elements), K ∈ ℝ n×n is the reduced stiffness matrix computed using finite element analysis, and A ∈ ℝ n×m ( n ≥ m ) is a linear operator that defines the contact pairs.
As a rule, contact problems arising in assembly simulation have a number of unknowns that varies from 1000 to 20,000. The mechanical properties of the structure, the influence of assembly jig type, and the configuration of the joining parts are described by the stiffness matrix K and the constraint matrix A . When performing a series of computations, these matrices do not change. This feature is exploited to define some auxiliary matrices based on matrices K and A once for every assembly model. Matrix K is symmetric, positive-definite and ill-conditioned. It has blockdiagonal structure with fully populated blocks. The number and the size of the blocks are defined by the number of assembled parts and the number of considered degrees of freedom in computational nodes in each part. Matrix A , which defines Fig. 2 Cross section of the assembly the pairs of nodes that can come in contact, is sparse in the considered problems: every column of the matrix contains 1 or 2 non-zero elements.
The analysis of assembly quality involves the simulation of shape variations of the joined parts. The shape variations are caused by tolerances in manufacturing, positioning, fixating, etc. In industrial applications (e.g., see Lupuleac et al. 2019b) the variation simulation is realized by generation of the large number of initial gap vectors g (cloud of gaps) and subsequent massive solving of contact problem (1). Such an approach can be considered as generalization of the Method of Influence Coefficients proposed in Liu and Hu 1997.

Equivalent formulations of the QP problem
Considering that the matrices K and A are usually the same for the series of computations during aircraft assembly simulation and optimization, problem (1) can be reformulated in a form that is more attractive from a computational point of view. Specifically, the size of the QP problem can be reduced and the form of the constraint matrix A can be simplified. Consider two different reformulations of problem (1): dual problem and relative problem.
The dual formulation of problem (1) is given by where ∈ ℝ m is the vector of Lagrange multipliers corresponding to contact forces in the junction area, Q = A T K −1 A ∈ ℝ m×m is a symmetric positive-definite fully populated matrix, p = A T K −1 f − g ∈ ℝ m and s = − 1 2 f T K −1 f ∈ ℝ 1 (Lupuleac et al. 2013).
The relative formulation of problem (1) is based on a vector of relative displacements in the junction area defined as u = A T x ∈ ℝ m and may be obtained through the minimax form of the problem. First, an additional relation should be derived from optimality conditions. By multiplying the stationarity condition Kx − f + A = 0 with A T K −1 we get Excluding the relation that binds variable x to u is obtained: Then the minimax form of the problem is Calculating the maximum, we obtain a relative formulation of problem (1): All the listed formulations of the QP problem are equivalent and their solutions are connected by the following formulas: Two more formulations arise when applying the Newton projection method and pivotal algorithm for aircraft assembly problems. The Newton projection method is only applicable if the constraint matrix is the identity (see Sect. 3.3 for details). To solve the primal problem (1) with general constraints A T x − g ≤ 0 a change of variables is proposed. If the total number of constraints is less than the number of variables, extra rows are added to the matrix A T to make it square. The proposed idea is to add unit vector rows in such a way that the new matrix becomes nondegenerate. This technique allows to keep the constraint matrix sparse. Thus, a new variable , … , a n T x is introduced and a change of variables in the func- ) is made. Finally, fictitious constraints a m+1 , … , a n T x ≤ +∞ are introduced, and problem (1) is reformulated as follows: The pivotal algorithm considered in the paper is used to solve linear complementarity problems LCP (q, M) that consists in finding the vectors w, z ∈ ℝ N such that for given a vector q ∈ ℝ N and a matrix M ∈ ℝ N×N or in concluding that no such pair of vectors w, z exists; N denotes the dimension of the LCP (q, M ). How the problems (1)-(3) can be formulated as LCP the reader may find in Sect. 3.4.
The basic information about the formulations (1)-(4) such as the number of variables, the number of constraints, type of constraints and structure of Hessian, is presented in Table 1. The details of the LCP formulations are described in Sect. 3.4. (3)

Specifics of parallelization for considered problems
The most important computational feature of the contact problem arising in variation simulation is the need to carry out many similar computations for the same large-sized assembly model; namely, to calculate the displacements for different sets of forces F = {f j } p and for different variants of the initial gap G = {g i } k . The total number of the serial computations can reach up to 10 6 runs during variation simulation and assembly optimization for each assembly model. Usually, both variation simulation and assembly optimization can be parallelized by tasks (Pogarskaia et al. 2018): each task refers to solving a particular QP problem (see Fig. 3). Thus, parallelization of QP solvers is not needed and is therefore not considered in this paper.
Task parallelization gives a rise to another computational feature. Namely, the effect of the dependence of computation time on input data for different solvers is enhanced. For some numerical methods the number of active constraints at the optimum n act affects computation time significantly, thus becoming an essential characteristic of optimization problem when computation time is analysed.  Table 2 presents the dependence of the number of active constraints n act for problem (1) on input data. An increase in the number of installed fasteners leads to an increase in the number of active constraints n act . Variation of the vector of initial gap generally does not affect n act much. Note that in the remainder of the text n act is the number of active constraints for problem (1), (3) and (4) and for problem (2) is equal to m − n act .

Adaptation of QP methods for solving the contact problem
The paper is devoted to adaptation and comparison of various types of quadratic programming methods applied to the contact problem. The following four widely used quadratic programming methods are considered: the interior-point method, the active set method, Newton projection method, and complementary pivot algorithm. The interior-point method is a polynomial-time algorithm that has proved its efficiency in practice (Gondzio and Grothey 2006;Byrd et al. 1999;Mehrotra 1992). It allows to solve a wide range of optimization problems such as linear, convex quadratic-programming, second-order cone programming, and semidefinite programming problems (Wright 2005). This method is known to be one of the main methods used for solving large sparse problems and is implemented in many commercial optimization software packages (MOSEK, CPLEX, IMSL, MATLAB and others).
Another class of optimization methods is active set methods. One of the most commonly used active set methods for solving quadratic programming problems was proposed by Goldfarb and Idnani in (1983). It is a dual active set method which proved to be fast and efficient for small to medium-sized problems (up to ten thousand unknowns), see Burton and Toint (1992); Gaspero et al. (2011);Marron et al. (1997). This paper focuses on the modification made by Powell (1985) that is preferable for ill-conditioned problems (further referred to as ASM). This method is also implemented in some well-known numerical libraries such as IMSL Numerical Library (Rogue Wave Software 2016) and Scilab (2018). Newton projection method is a second order modification of the gradient projection method, which is an extension of the steepest-descent method for constrained minimization problems. The gradient projection method was proposed by Goldstein (1964), Levitin and Polyak (1966) and Newton projection method was thoroughly developed by Bertsekas (1976Bertsekas ( , 1982. It is easily implemented and highly effective if constraints have the form x − g ≤ 0 , as the projection operator is just a minimum operator in this case. Newton projection method works well for large-scale problems, but is not popular with researchers due to the fact that its basic version has a low convergence rate. Convex quadratic problems can be formulated as linear complementarity problems. One of the most commonly used methods for solving problems of this type is complementary pivot algorithm (CPA) or Lemke's method, introduced by Lemke (1968). Lemke's method is widely used because of such features as finiteness, reliability and efficiency (Cottle et al. 1992;Acary and Brogliato 2008). Several numerical packages include implementations of Lemke's method such as Siconos (Acary and Pérignon 2007) or GAMS PATH (Ferris and Munson 2018).
The theoretical estimations of computational complexity do not allow choosing the indisputable leader among these methods for the contact problem solving (see Table 3). Note that Table 3 presents a comparison for the primal formulation of the contact problem (1). For the dual and relative formulations (problem (2) and (3), respectively), the number of blocks p in the stiffness matrix K is 1 and the number of unknowns n is equal to the number of constraints m.
All algorithms have polynomial time on average with a close degree of the polynomial. In particular, computation time depends on the number of active constraints at the optimum n act and the number of conjugate gradient method iterations n cg when iterative method is used for linear system solving. Typically, the number of active constraints n act can run over [0, m] for an assembly model. It can take small values as well as large ones when solving a set of similar contact problems as n act is defined by the vector of applied forces f and the vector of initial gap g ( Table 2). Number of iterations n cg depends on the condition number of the stiffness matrix K and the quality of the applied preconditioner. Typically, the condition number is 10 5 − 10 8 for the considered type of contact problems.
While solving a set of similar contact problems, some of the time-consuming operations can be performed once. These operations include the Cholesky decomposition of the stiffness matrix K , the inversion of the Cholesky matrix, and the reordering of constraints and stiffness matrices. The time required for such data preprocessing is presented in the third column of Table 3. It differs for particular algorithms depending on required operations. The details of implementation and adaptation to the features of contact problem are discussed in the following subsections.

Interior-point method
The primal-dual interior-point method (IPM) reformulates the problem of constrained minimization (1) as a sequence of unconstrained minimization problems using a logarithmic barrier: where is the barrier parameter and a j is the j-th column of A . The solution of problem (6) converges to the solution of problem (1) when goes to zero. At each iteration IPM solves a system of nonlinear equations (7) that represents the first order optimality conditions for problem (6): In the system above, and Y are diagonal matrices with and y = −A T x + g on their diagonals and e = 1 m is a vector of ones. One step of Newton's method is used to find an approximate solution of system of linear equations (7) and then the duality gap is updated. The process is repeated until becomes sufficiently small.
There are two challenging issues related to the implementation of the interior-point method when solving the contact problem: the selection of a starting point and the solution of the linear system for determining the Newton direction. There are two variants of the algorithm, namely feasible and infeasible, which call for the fulfilment of different conditions on the starting point. Feasible IPM requires starting points that satisfy the following conditions:  ( 1 ) whereas for infeasible IPM, a starting point needs only satisfy the inequality conditions (10). So, for the infeasible IPM the starting point search process is simplified considerably. Note, that the theoretical worst case complexity is better for the feasible IPM, although the infeasible IPM is known to have practically more efficient implementations (Gondzio and Terlaky 1995). The choice of starting point for solving contact problem is discussed in Stefanova et al. (2018), where the starting point for feasible IPM is determined based on the physical interpretation of the QP problem (1), and is compared to the infeasible one with the starting point from D'Apuzzo et al. (2010). When applying IPM to assembly problems, a feasible IPM is preferred for assembly models with junction area of simple type, whereas for more complicated sandwich type joints the algorithm described in Stefanova et al. (2018) does not have direct physical justification; an infeasible IPM yields better results in this case (see Table 4). The most time consuming IPM step is the solution of the linear system of equations for determining the Newton direction: where I is the identity matrix, , and is the centering parameter chosen adaptively (Mehrotra 1992). A challenging numerical issue is the increase of condition number as the optimal solution is approached. To solve the system of equations (11), one can use direct or iterative approaches. Direct methods are based on decomposition techniques (Gondzio and Grothey 2006) and can be efficient if the problem is sparse or parallelization is possible. When the problem is large-scale, an iterative approach with the use of an appropriate preconditioner is preferable. Several different approaches for a set of contact problems arising in assembly simulation are discussed further and a preconditioner attractive for the considered type of problems is proposed.
The system of equations (11) is first transformed to an augmented system (12) and then to a normal system (13): Δx and Δy can be computed using the relations Δx = K −1 (s 1 − AΔ ) , Δy = −1 (s 3 − Y ) . The augmented system is symmetric indefinite. For the class of problems considered in the paper, solving the system of equations (12) turned out to be significantly less efficient compared to solving the normal system, therefore it is not considered in the paper. The normal system has smaller dimension; it is symmetric, positive definite and fully populated. The system of equations (13) is solved directly using Cholesky decomposition and iteratively using preconditioned conjugate gradient method (CG).
Three preconditioners are considered and compared in this paper. The first one is a Jacobi preconditioner P 1 = diag(H) that is based on the diagonal elements of matrix H . The second one is an incomplete Cholesky preconditioner P 2 = L H L H T . The incomplete Cholesky factorisation is applied to the matrix H r = H s + I , where is a regularization parameter and H s is a sparse approximation of matrix H . Paper (Lin and Saigal 2000) presents a brief review of several existing approaches to define sparsity pattern for dense matrices based on dropping elements with small magnitude from H . The incomplete Cholesky preconditioner uses more information about matrix structure compared to P 1 ; however its efficiency depends on the parameters such as the number of non-zeros in H s and the value of regularization parameter . While the number of CG iterations is reduced with the increase of the amount of non-zeros in H s , the time required to build the preconditioner increases. By varying the number of non-zeros, the preconditioner with 50% of the largest elements of the original matrix H (and all other elements substituted by zeros) was chosen as the fastest. The regularization parameter is chosen adaptively starting from = 0 ; the regularization parameter is increased if the incomplete Cholesky factorization fails.
The third preconditioner proposed for assembly problems (Stefanova et al.

2018) is a split preconditioner
the Cholesky decomposition of the matrix Q = A T K −1 A . The strong point of the HΔ = p, proposed preconditioner P 3 is the combination of the possibility to be updated quickly at each IPM iteration with a good approximation of H . The matrix L Q does not change from iteration to iteration, so only the diagonal terms √ −1 Y must be modified to update the preconditioner. Moreover, the Cholesky decomposition L Q can be performed once per assembly model.
All of the described approaches are compared in Table 5. The table presents the time necessary to solve contact problem. The computations were done on an Intel Core i5 CPU 3.30GHz computer with 16 GB of RAM. Two sets of problems with different vectors of forces f are considered for wing assembly models with 4386 and 6112 variables. The computational results show that the iterative approach with preconditioner P 3 is considerably faster.
Thus, to solve the system of linear equations the adapted IPM uses CG method with the combined preconditioner P 3 and an appropriate starting point is chosen depending on the assembly model.

Active set method
The main goal of active set methods is to find the set of constraints that are active at the solution point: An initial approximation of the active set I 0 act is chosen and then one constraint is added or removed at each iteration forming the current active set I k act . Since the number of constraints is finite, iteration process finishes after finite number of steps (Boland 1996). The dual Goldfarb-Idnani active set method, considered in this paper, starts at the point of the unconstrained minimum of problem (1), i.e., x 0 = K −1 f , I 0 act = �. If the non-penetration conditions are violated, the iteration procedure starts. At each step, current point x k is forced to satisfy one of the violated constraints in such a way that (A T x k − g) i k = 0 , and the objective function of problem (1) is minimized over the set I k act ∪ i k . The process continues until the point x k is in the admissible set, meaning that the solution is found.
In the case that it is not possible to satisfy any constraint during the iteration step, the rollback is performed: the point is shifted in such a way that one of the current active constraints is violated again.
It is worth mentioning that the iteration step is made in both primary and dual spaces (not only the point is shifted but also the Lagrange multipliers).
In order to speed up the computations several modifications of the described algorithm are proposed that make use of the input data structure given in the paragraph 2 and take the algorithm specifics into account.
1. Storing only the non-zero blocks of the matrix K allows the reduction of the time necessary for the calculation of the objective function and its gradient, because the number of multiplications is decreased from (n ⋅ p) 2 to p ⋅ n 2 , where p is the number of blocks and n is the number of variables in the block. 2. At each iteration step, the next constraint to be satisfied is chosen. In the original version of algorithm, the first most violated constraint is taken into consideration: i next = arg max i=1,m a T i x − g i . A more complicated procedure is introduced for the described algorithm relying on the fact that if a force is applied to a computational node, it means that a fastening element is installed. The gap is then assumed to be closed, and the corresponding constraint must be added to the current active set. 3. As it was mentioned earlier, matrix A T is sparse: each row contains either 1 or 2 elements. Therefore, it is not necessary to process the entire matrix, but only the indices of the non-zero elements, which leads to the decrease of the number of multiplications by a factor of 10: multiplication of matrix A T by a vector takes O(n) operations.
The computation time after the implementation of all the proposed modifications is given in Table 6 for a test problem with 3916 variables. The computations were done on an Intel Core i7 CPU 4.00GHz computer with 32 GB of RAM.
Another important feature that can be implemented with the help of the active set method is a "warm-start" technology (Goswami et al. 2012). When a series of similar problems is solved, it is reasonable to start not from the point of the unconstrained minimum but instead from the solution obtained for the previous problem. Thus the iterations that were already done can be skipped and the speed of the computations is significantly increased. Table 7 illustrates how the computation time changes depending on different starting points. Therefore, the Goldfarb-Idnani active set method can be easily adapted to the peculiarities of the considered problem. Implementation of the modifications listed above allows to reduce computation time drastically. The technology of "warm-start" can be efficiently applied when optimizing the assembly.

Newton projection method
Newton projection method (NPM) is an iterative technique based on gradient descent for solving constrained optimization problems. Each iteration consists of a step in a search direction and a projection of the new point to the feasible region. The process continues until the gradient projection is close enough to zero. An iteration of NPM is defined as where d (k) is the search direction, k is the step size, P(x) is the projection operator.
Constraints with an identity matrix are considered: Relative and dual problems satisfy this condition as well as primal problem (4). The projection operator is simple: The method for choosing the direction d (k) with quadratic rate of convergence was introduced in Bertsekas (1982). It is based on Newton's method. For the algorithm specification the set of indices is introduced where The set I − (x (k) ) includes such indices i that x i would not change if the search was made along the antigradient. These indices are excluded from the search and the Newton's method is used in the subspace defined by the remaining indices to reach the quadratic rate of convergence. Thus, the search direction d (k) can be defined from where Initial estimate x (0) can be set equal to g in order to make the set I − (x (0) ) not empty and reduce the dimension of K (0) P . Solving the Eq. (15) is a computationally expensive part of NPM. Thus, choosing an effective method for solving (15) is a key point in NPM adaptation to assembly problems. Notice that it is possible to avoid storing of zero rows and columns of matrix K (k) P . Thus, the dimension of the system of equations is reduced by n (k) act = | | I − (x (k) ) | | . Three approaches based on Cholesky decomposition, conjugate gradient (CG) method and their combination are compared.
First, Cholesky factorization method is considered. Since every Eq. (15) is a reduction of the full system where some rows and columns are disposed, it is suggested to not recalculate the factorization at each iteration but instead update it from the previous iteration using relations presented in Stewart (1998) and Osborne (2010). This method can significantly reduce the time needed for solving the system of equations (15). The system of equations (15) can also be solved with iterative CG method (van der Vorst 2003). Another proposed idea is to use a combined approach, i.e., to use the CG method with low accuracy at the first iterations of NPM, then increase the accuracy during later iterations, and switch to the Cholesky decomposition when the CG method starts slowing down.
The comparison of the described approaches is presented in Table 8. The computations were done on an Intel Core i5 CPU 3.30GHz computer with 16 GB of RAM. Dimension of the Eq. (15) is n − n (k) act ≈ n − n act . Thus, the total time for obtaining the solution is O (n − n act ) 3 . Cholesky decomposition update reduces NPM solution time, especially for the relative problem. CG method can work either faster or slower than Cholesky update depending on the matrix. Combined method commonly works faster than both of them.
One of the challenges of NPM implementation is the step size selection. There are several common methods for this problem. The first method is a constant step method (Goldstein 1964;Levitin and Polyak 1966), which requires knowledge of the Lipchitz constant and has a poor convergence rate. The second one is Armijo step size rule (Bertsekas 1982). The third method is the one-dimensional functional minimization (McCormick and Tapia 1972): which can be performed by some one-dimensional minimization algorithm. The golden section search method can be used for the minimization. Practically, the minimum is usually either 1 or close to 0, so the golden section method is used for the logarithm of in order to reach values close to zero in less iterations. Another approach is to use the projected line search (PLS) (Šantin et al. 2016), which is actually a piecewise quadratic minimization. The comparison between the Armijo rule the golden section method and PLS is presented in Table 9. PLS works faster than golden section method and much faster than the Armijo step size rule.
The comparison between the primal, relative and dual problems is presented in Table 10. The number of iterations of NPM for primal and relative problems decreases linearly with respect to n act and has a linear dependence on the problem min ∈(0,1] size (Šantin and Havlena 2011). Thus, the number of iterations can be estimated as O(n − n act ) and the time needed to solve the problem as O((n − n act ) 4 ) . Primal and relative problems are solved faster when the number of active constraints n act is large, while dual problems are solved faster when n act is small. Moreover, relative problems are almost always solved much faster than the primal ones as the size of the relative problem is smaller. Unfortunately, the dual problem has a worse convergence rate than primal and relative problems due to the different physical meaning of the variables and the lack of diagonal dominance of the inverse stiffness matrix for the dual problem.  In the implemented method the combination of Cholesky factorization, Cholesky update and CG method is used for solution of (15) and the golden section method is used for the step selection. The dual problem should be solved if the number of active constraints is relatively small (i.e., a small number of fasteners is applied), and the relative problem should be solved otherwise.

Complementary pivot algorithm
Lemke's method or complementary pivot algorithm (CPA) is a finite numerical procedure for solving LCPs. The implementation of the CPA is based on the idea of pivots. First, an artificial variable is introduced to the system to replace one of the basic ones. Then, pivots between the variables are performed until the artificial variable is dropped out of the basis. The resulting vector of basic variables is then considered to be the solution.
Convex quadratic programming problems can be reformulated as LCPs (5). The Karush-Kuhn-Tucker conditions (KKT) for the primal (1) and dual (2) problems can be written with an appropriate transformation of variables as follows: Conditions (16)  Both (17) and (19) LCPs incorporate primal (1) and dual (2) problems. The dimension of the constructed problems N = n + m . The CPA, being a pivoting method, is highly sensitive to the dimension of the problem. Therefore another way of LCP construction is suggested in (Yassine 2008). The quadratic problem without positivity constraints can be transformed into an LCP (q, M) where This LCP can also be obtained if the dual problem (2) is chosen as a base for KKT conditions and hence this LCP is further referred to as dual LCP, while the previous LCPs (LCP (17) and LCP (19)) are referred to as primal. The dimension of the dual LCP (20) N = m . Since the average computation time for the CPA is O(N 3 ) , the dual LCP has an important advantage over the primal ones.
The relative problem (3) can also be reformulated as an LCP. The basic version of the primal LCP (17) and the modified version (19) for the relative problem have a dimension N = 2m , and the dual LCP (20) for the relative problem is the same as for the primal, because the dual quadratic problems for them are identical.
Comparison of computation time of CPA application to primal, dual and relative problems is presented in Table 11. The computations were done on an Intel Core i5 CPU 3.30GHz computer with 16 GB of RAM. Results provided in the table confirm that Lemke's method for the dual LCP (20) works decisively faster than for the relative or primal problems.

Computational results for aircraft assembly
In this section, the computation time of the optimization methods is compared for a set of typical aircraft assembly problems. The wing-to-fuselage assembly is chosen for this purpose. The details of wing-to-fuselage assembly simulation are given in Lupuleac et al. (2018Lupuleac et al. ( , 2019b. The upper and lower outer wing box panels are joined with the corresponding central wing box (CWB) panels (see Fig. 4). The CWB is nested within the central fuselage section. A schematic outline of the wing-to-fuselage assembly is presented in Fig. 5. The CWB is designed to be the most rigid part of the aircraft, and therefore the elements of the assembly joined with CWB can be regarded as fixed. The assembly jig supports the wing at several points. The additional parts (buttstraps) are used to ensure the tightness of the joint. The assemblies for the upper and lower wing box panels are performed independently and interinfluence between these processes is minimal. That is why independent models for the upper and lower wing-to-fuselage assemblies are considered (see Fig. 6).
The upper wing-to-fuselage joint consists of just two parts: the cruciform and the upper wing box panel. The cruciform is already fastened to the CWB, so it can be assumed to be fixed along the edges as it is shown in Fig. 5. The model for the upper wing-to-fuselage joint is of simple type and has one junction area between the wing and the cruciform. The lower wing-to-fuselage joint is more complex compared to the upper one and consists of a lower wing panel, a triform, a part of the CWB, and several buttstraps (Fig. 7). The triform and the buttstraps are simultaneously fastened to the wing and to the CWB panel. Hence, the lower wing-to-fuselage assembly model is of sandwich-type where the lower wing panel and the part of CWB are located above the buttstraps and below the triform (Fig. 5). The triform is joined with the vertical CWB panel, and thus it is fixed along the upper edge. The buttstraps are not fixed on the assembly jig before fastening.
An example of computational results for an upper wing assembly model is shown in Fig. 8. The left figure presents the initial gap field between the cruciform and the upper wing box panel. The right figure shows the residual gap after the installation of 3 fasteners.
Three test problems used for comparison of algorithms are described in Table 12. Two models with 5206 and 11308 unknowns correspond to the upper joint model (UJM). The model with the 6190 unknowns refers to the lower joint model (LJM). The condition number of stiffness matrix K is larger for LJM than for UJM due to the features of the assembly jig fixing. The algorithms are implemented in C++ using MSVS 2013 and Visual C++ 12.0 compiler. All computations presented in this section were done on an Intel Core i5 CPU 3.30GHz computer with 16 GB of RAM. As seen in Sect. 3, each algorithm used its own stopping criteria. However, these criteria are chosen so that the difference between the displacements obtained by the two algorithms does not exceed 10e−7 mm. Figure 9 presents comparison of the optimization methods for UJM and LJM in terms of computation speed. The computation time is shown in logarithmic scale as a function of the number of active constraints. The results are presented only for dual and relative formulations since for the primal formulation the time is much longer due to the size of the problem. Each computation was repeated 10 times and the average value is shown in the Fig. 9. The time deviation varies 4-15% in average.  Let us first consider the results for UJM1 problem. Computation time of NPM and CPA for the dual problem and ASM for the relative problem increases with the increase in the number of active constraints at the optimum, i.e., in the number of installed fastening elements. Time of NPM for the relative problem and ASM for the dual one has the opposite dependency. IPM computation time does not depend strongly on the number of active constraints. Note that this behaviour agrees with the theoretical complexity presented in Table 3. These basic time trends are also present in the results of UJM2 and LJM problems.
The UJM1 and UJM2 differ mainly in the number of variables n. Thus, its influence is clearly observed in the results. With an increase in the number of variables, the time increases significantly for the most of the methods with the exception of IPM. NPM applied to the dual problem is generaly the fastest in solving the contact problem for upper wing models. However, NPM for the relative problem is less time consuming than for the dual one when the number of active constraints is large.
The UJM1 and LJM have similar number of unknowns, but the condition number for the lower wing model is larger. Comparing these two models, ASM is revealed to be the least affected by the condition number of the stiffness matrix. ASM applied to the relative problem is usually the fastest for LJM with the exception of problems with a large number of active constraints, where NPM for the relative problem is the leader.
Finally, let us consider the computation time distribution for variation simulation analysis of an Airbus A350-900 wing-to-fuselage assembly. The aim of the   Fig. 9 Comparison of implemented methods. Solid and dashed lines refer to dual and relative problem formulations respectively analysis is to verify that a given configuration of fastening elements provides a sufficient fastening level to close the gap between the assembled parts (Lupuleac et al. 2019b). Such an analysis requires to solve the contact problems for a set of measured initial gaps. In total about 200 sets of measurements of aircraft produced in 2013-2016 were used for the analysis. The distribution of the number of active constraints for the measured initial gaps is shown in Fig. 10. UJM1 with 50% of installed fasteners was used for variation simulation analysis, i.e the contact problem was solved with the same stiffness matrix, matrix of constraints and vector of forces for all of the measured initial gaps. Figure 11 presents computation time distribution for different solvers. The results presented in Fig. 11 agree with those presented in Fig. 9: the order of the bars in the histograms follows that of the curves for the corresponding number of active constraints. Since the computation time of the CPA for relative problems is considerably greater than for dual ones (see Sect. 3.4), it was only applied to the dual formulation of the test problems.
Note that some problems arising in fastener pattern optimization can require massive computations of up to 10 6 runs and thus task parallelization is needed  (Pogarskaia et al. 2018). Therefore, to reduce the total computation time of assembly simulation, not only the average computation time for one run has to be reduced, but also the maximal computation time. According to the obtained results (see Fig. 11), for both dual and relative problem, NPM is revealed to be the leader both in terms of average and maximal computation time. The CPA works considerably slower than the other implemented methods for the test problems. Moreover, it has the largest deviation between maximal and average time, which creates difficulties for task parallelization. Note that the results in Fig. 11 only represent the variation simulation analysis of one fastener pattern. Different methods can yield the best results for different patterns.

Conclusion
The paper presents and discusses computational aspects related to the solution of a certain class of quadratic programming problems related to variation simulation analysis.
The contact problem arising in assembly simulation is reformulated in dual and relative form. Data preprocessing (e.g., matrix transformations) is only performed once for the assembly model. Therefore the time for each of the multiple problem solving as well as time for the complete assembly analysis is reduced, since the variation simulation analysis requires solving of problems with similar data.
The presented numerical experiments demonstrate the following features of the methods. The computation time of the interior point method is not affected by the number of active constraints at the optimum. Thus, this method is suggested for such problems as optimization of fastening pattern for simple joints where the number of active constraints may vary considerably. The Goldfarb-Idnani active set method is revealed to be the best for middle-scaled ill-conditioned problems such as verification analysis of complex joints. Newton projection method is suitable for problems of various types, to the fullest extent for the large-scale ones and for problems with big number of active constraints. The Complementarity pivot algorithm should be preferred for problems with small number of active constraints. For the example of variation simulation analysis of an assembly simulation problem, Newton projection method was revealed to be the leader in terms of average and maximal computation time.
The adaptation of the discussed numerical methods combined with an appropriate modification of the contact problem leads to a considerable reduction of computation time for the described problems.