Solution of a truss topology bilevel programming problem by means of an inexact restoration method

Solution of a truss topology bilevel programming problem by means of an inexact restoration method

Ana Friedlander; Francisco A. M. Gomes

Departamento de Matemática Aplicada, IMECC, Universidade Estadual de Campinas, 13083-859 Campinas, SP, Brazil E-mails: [e-mail] / [e-mail]

ABSTRACT

We formulate a truss topology optimization problem as a bilevel programming problem and solve it by means of a line search type inexact restoration algorithm. We discuss details of the implementation and show results of numerical experiments.

Mathematical subject classification: Primary: 65K05; Secondary: 90C55.

Key words: inexact restoration, bilevel programming, truss topology optimization.

1 Introduction

Bilevel programming problems model two-level hierarchical systems and have been studied since the seventies. In [5, 6, 7], history, applications, algorithms, theoretical questions and almost all relevant references can be found.

In this paper we address the following bilevel programming problem:

where x ∈ Rⁿ, s ∈ R. u, u ∈ R^m , and W, q and g are real functions.

The truss topology optimization problem with frictionless unilateral contact is naturally formulated as the mathematical bilevel problem stated above. In this problem, we wish to maximize the stiffness of a truss, given a bound on the material volume. This volume depends on the cross section area of the bars, that are the design variables of the problem. In the lower level problem, whose variables are the displacements of the nodes in the structure, the potential energy of the truss is minimized, given a fixed area for each bar. The truss structure may come in frictionless unilateral contact with some given supports, which is modeled defining bounds on the displacements.

For details on contact problems see [11]. See [4] for the general truss optimization problem, [12, 13] and references therein for more details on the truss optimization problem with unilateral frictionless contact. The methods used to solve this problem involve essentially non-differentiable techniques for the bilevel problem or SQP algorithms for the reformulation of the problem that results after substituting the lower level problem by its KKT conditions. In the latter case, the fact that the lower level problem is a minimization problem is ignored, increasing the chances to find stationary points that are not minimizers.

In [2] the authors studied the resolution of bilevel programming problems using the inexact restoration (IR) algorithm introduced in [15]. In this paper we follow the approach for the IR algorithm introduced in [8].

The present paper is organized as follows. In Section 2, we describe the IR algorithm and discuss its application for solving bilevel problems. In Section 3, we state and comment the truss topology optimization problem that will be solved. Since the resolution of bilevel problems requires specifications that are problem dependent, we discuss in Section 4 how this IR algorithm can be adapted to our problems. In Section 5, we present numerical experiments and we state some conclusions in Section 6.

2 The line-search IR algorithm

A general nonlinear programming problem can be stated as

where f : Rⁿ -> R and C : Rⁿ -> R^m.

Inexact restoration algorithms for nonlinear programming (see [9, 15]) belong to the class of feasible algorithms because they emphasize the importance of the feasibility of an approximate solution. They were introduced to avoid the excessive computational effort that must be spent to achieve constraint fulfillment at iterations that are still very far from optimality, in problems with strongly nonlinear constraints.

In [8], the authors introduce a line search IR algorithm that captures the essential features of the IR technique. Moreover, its convergence theory is straightforward and allows for considerable freedom in the choice of implementation details.

The basic steps of the algorithm presented in [8] are shown below, assuming that r ∈ [0, 1), β,y, τ > 0 are fixed parameters, h : [l, u] —s- [0, ∞) is a function such that ||C (z)|| < h (z), and Φ(z, p) = f (z) + (1 - p)h(z).

Algorithm 1

• Step 0: Initialization.

Choose l < z ⁰ < u and p₀∈ (0, 1). Set k := 0.

• Step 1: Inexact restoration. Compute y^k ∈ [l, u ] so that

• Step 2: Penalty parameter determination. Determine p_k+₁∈ {2^-lp_k, i = 0, 1, ...} as large as possible so that

• Step 3: Computation of the search direction. Compute d^k so that y^k + d^k ∈ [l, u] and the following inequalities hold for t ∈ [0, τ],

• Step 4: Line search. Determine t_k ∈ {2^-i, i = 0, 1, ...} as large as possible so that the first inequality in Step 3 holds and

• Step 5: Step update. Set z^k+¹ := y^k + td^k and k := k + 1. Go to Step 1.

The equality constraints that define the feasible region of the bilevel problem (1) are given by

At the inexact restoration step, we need to determine a point that is more feasible than (x, s, u)^T. For q (x, s), this is clearly stated in Step 1 of the algorithm. However, for the low level problem

we have to explain carefully what we mean by a better value of u.

At the beginning of iteration k, the point z^k = (x^k, s^k, u^k, μ^k_Lμ^k_U)^T is available, where μ^k_L and μ^k_U are estimates of the Lagrange multipliers associated with the bounds on u. For a fixed value of the first two variables, we can use u^k as an initial point and apply an adequate minimization algorithm to solve problem (4). The choice of the algorithm is problem dependent. For the truss topology optimization problem presented in Section 3, we discuss this issue in Section 4.

The inexact restoration step of the algorithm does not require full feasibility. Instead, a sufficient decrease of the feasibility error measured by a majorant of the norm of the constraint vector C (z) given in (2) is required. In the definition of C (z), the infeasibility related to (4) is measured by equalities in the KKT optimality conditions of this minimization problem, i.e.

It is very important to realize that the lower lever problem (4) is not replaced by its KKT conditions. For given values of x, we minimize the original problem in u and use the satisfaction of the KKT conditions with an adaptive tolerance as a stopping criterion for this minimization process.

The bound constraints and the non-negativeness of the Lagrange multipliers are forced at each iteration. Other details of the algorithm, such as the choice of h, the search direction in Step 3 and the line search in Step 4 are problem dependent and will be presented in Section 4.

3 The truss topology problem

The simplest truss topology optimization problem consists in finding the stiffest truss for a given volume. Various formulations of this problem can be found in [1] and [3]. More sophisticated problems can be generated if, for example, we include upper limits on the displacements and on the stresses, or multiple loads.

In this work, we consider the optimization of a truss subject to external loads and frictionless contact supports. We restrict our analysis to a rectangular bidimensional domain discretized into a mesh with n_x x n_y nodes. The truss bars are generated connecting nodes from this mesh. The set of potentials bars is called the ground structure. The size of the ground structure depends on the connectivity level adopted, that determines to which other nodes each node can be connected.

If a node can only be connected to its immediate neighbors, the connectivity level, hereafter denoted n_p, is set to one. If a node can also be connected to the nodes in the neighborhood of its immediate neighbors, then n_p = 2. Generalizing this idea, n_p can assume any value between 1 and the maximum of {n_x — 1, n_y — 1}. In the latter case, the ground structure is said to be fully connected. Figure 1 shows all of the bars of a ground structure that connect one node to its neighbors in the NE quadrant, for n p = 4.

Table 1 shows the growth of the number of bars as a function of the number of nodes and the connectivity level, for a square structure. The number of bars of a fully connected base structure with n_x = n_y = 100 is 30398894. Thus, for highly discretized structures, it is important to keep n p small.

3.1 The truss topology problem with frictionless contact as a bilevel programming problem

In the topology optimization of a truss formed by m potential bars that link n potential nodes, we consider two sets of variables: the vector of cross sectional areas of the bars, x ∈ R, and the vector of horizontal and vertical displacements of the nodes, u ∈ R²ⁿ. The objective is to maximize the stiffness ofthe structure, i.e. minimize the compliance, subject to volume and contact constraints.

Following the min-max formulation of the standard truss topology optimization problem presented in [3,4], the problem with frictionless contact conditions can be stated as

where the vector lcontains the lengths of the bars, V is the upper limit for the structure volume, s is the slack variable of the volume constraint, x_U is the vector of upper limits for the cross sectional areas of the bars, u_L and u_U are the lower and upper limits for the nodal displacements and P is the vector of nodal forces. The global stiffness matrix K (x) is given by where E is the Young's modulus of the material and bj is the j-th column of the compatibility matrix B that relates the nodal forces to the bar forces. More information on the formulation of the truss topology problem can be found in [3].

The optimal solution of (6) contains only a few of the ground structure bars. Therefore many components of vector x are expected to vanish near the optimum. However, to avoid numerical difficulties that arise when the cross sectional areas of the truss bars are allowed to become zero, it is a common practice to define a positive lower level x_min for these areas. Following this practice, we also require that xi ∈ [xmin, xmax], for i = 1, ..., m.

4 The algorithm of inexact restoration for truss optimization with friction-less contact

The crucial steps of Algorithm 1 are the inexact restoration, the definition of the search direction and the line search. In this section, we show how these steps may be efficiently computed for the bilevel problem (6).

4.1 The inexact restoration step

An approximate solution z for (6) is said to be feasible if

where z = (x, s, u, μ_L, μ_U)^T, μ_L and μ_U are the Lagrange multipliers of the lower level problem in (6), s is the slack variable of the volume constraint,

Here, M_L and M_U are the diagonal matrices containing μ_L and μ_U, respectively. One may notice that the first component of C (z) is the volume constraint, while the remaining components are the equalities in the KKT conditions of the lower level problem.

In the inexact restoration phase, we seek a solution that satisfies the conditions stated at Step 1 of the algorithm. In [8] the authors suggest that this can be accomplished defining

and requiring the restored point at iteration k (the point y^k) to satisfy h(y^k) < rh(z^k) and

The bounds on the variables defined by Ω and the linear volume constraint are forced at each iteration. Therefore, the components x^k and s^k ofthe restored step y^k take the same values as in z^k. This implies that the restoration is restricted to find a vector of nodal displacements u^R that solves the problem

After obtaining u^R, the Lagrange multipliers are estimated and forced to verify ||μ_L^R - μ^k_L||∞< βh (z^k) and ||μ_U^R - μ_U^k ||∞< βh (z^k).

If we define a positive lower bound for the cross sectional areas of the bars, (9) turns out to be a box constrained strictly convex quadratic programming problem that can be easily solved.

4.2 The search direction

Algorithm 1 requires the search direction d^k = (d_x, d_s, d_u, d_uL,d_uU ) to satisfy

In [8], the authors show that this can be obtained taking d^k as the projection of -∇ W (x^k, u^R) onto the set of tangent directions

In this case, d^k is the solution ofthe quadratic programming problem

where

and U^R, U_L and U_U are diagonal matrices containing u^R, u_L and u_U, respectively.

Instead of using this sort of Cauchy step, in our algorithm, an improved direction vector is obtained replacing the objective function of (12) by the quadratic approximation of W(x^k + d_x, u^R + d_u), given by

where

Since the solution of this problem may not satisfy (11), we include a trust region constraint, so the step is recomputed whenever it fails to verify this condition. Thus, our search direction is the solution of

With this strategy all the convergence results are essentially preserved as shown in [10].

The solution of (13) may be obtained using any quadratic programming library that works with sparse matrices. However, if the routine requires the decomposition of a matrix formed of some or all of the columns of ∇C (y^k) or ∇C(y^k)^T ∇ C(y^k), it is better to reorder the variables in a way that the factors generated by the decomposition will remain very sparse. A good reordering can be obtained if we replace d^k by d^k = (d_uL, d_uU, d_u, d_s, d_x). Moreover, a good node ordering should also be used, in order to reduce the fill-in produced by the factorization of K(x^k).

Following the usual approach of IR algorithms (see [9] and [15], for example), we also tried to define d^k as the solution of the nonlinear programming problem that consists in minimizing the original objective function subject to a linearization of the constraints, i.e.

However, the direction generated by this problem, besides being more difficult to obtain, did not show a better performance than the one computed solving (13).

4.3 The line search

The backtracking scheme adopted in Algorithm 1 seems to be quite simple to code. However, to compute a step length that satisfies the condition stated at Step 4 of the algorithm, we need to evaluate the objective function at some intermediate points (x+, u+) = (x^k + λd_x, u^R + λd_u). This task may be costly if a naive approach is adopted, since it requires the computation of the stiffness matrix K (x+), that depends on the number of bars. Fortunately, K varies linearly with x, so we may write

Thus, once we have computed K (d_x), we can avoid building K from scratch each time λ is changed. Moreover, we can also reduce the cost of obtaining W+ = W (x+, u+) if we write

Vector v₁ = K(x^k)u^R is already available, since it was used to calculate W(x^k, u^R). Therefore, after computing vectors u+, v₂ = K(x^k)d_u + K(d_x)u^R and v₃ = K (d_x)d_u, the objective function value may be obtained with two saxpys and two inner products.

After determining the optimal value for λ, the stiffness matrix may also be updated using (15). To keep the roundoff errors under control, K should be recomputed from scratch once in a while. In our algorithm, this recalculation is done at each 50 iterations.

If λ = 1 satisfies the condition defined in Step 4, then the approach presented here is more expensive than computing K (x+) directly. However, usually the step length needs to be reduced, so it pays for computing three extra matrix-vector products and a matrix update.

Finally, it must be noticed that due to the linearity of the volume constraint and the fact that ∇C(y^k)d = 0, vector x+ will always satisfy this constraint whenever it is satisfied by x^k, independently of the value of λ. This property of the problem allows us to perform an inexact restoration based only on the nodal displacements.

5 Numerical results

In this section, we show the results obtained using the algorithm for nonlinear bilevel problems presented above to optimize the topology ofsome plane trusses supported by a frictionless contact foundation.

In all of the examples, the material used in the trusses has elasticity modulus E = 200000N/mm². The cross sectional areas of the bars are required to belong to the interval [0.00001mm², 4mm²]. The parameters used in the inexact restoration step are r = 0.5 and β = 10.0.

5.1 First problem

The ground structure of our first problem is shown in Figure 2. Due to symmetry, only the right half of the frame is considered. The structure has 668 potential bars and 189 potential nodes. The connectivity level is one.

Only one external vertical force of magnitude 10N is applied at node 9. The volume is limited to 2% of the ground structure's volume. The foundation is not parallel to the bottom line of the structure. In fact, the gap between node 91 and the foundation is equal to 10 mm, and this distance is linearly reduced to 0 mm at node 181 (the bottom right node). It must be noticed that no fixed supports are defined for the structure, but only the frictionless foundation that limits downward node displacements.

The solution of the first problem is shown in Figure 3. The truss obtained resembles the MBB beam that is usually obtained solving the topology optimization problem where the lower right corner of the structure is fixed.

Increasing the connectivity level to 2, we obtain the structure shown in Figure 4.

5.2 Second problem

The domain of the second problem is the same of the first problem. However, an external vertical force of magnitude 1 N is applied at node 1, instead of at node 9. The volume is limited to 2% of the ground structure's volume. Also, the frictionless foundation is flat and covers only one fourth of the ground structure basis.

A schematic representation of the problem is presented in Figure 5. The solution obtained after applying Algorithm 1 is shown in Figure 6. Again, the truss obtained is similar to the usual MBB beam with a fixed lower right node. However, the structure is shorter, due to the presence of the flat foundation.

5.3 Third problem

In the third problem, there are two external forces applied at the lower corners of the structure. A small frictionless foundation prevents the downward displacement of one single node of the structure's basis. The volume is limited to 0,05% of the ground structure's volume. The connectivity level is set to 2.

The problem is shown in Figure 7 and the solution obtained after applying Algorithm 1 is given in Figure 8. As we can see, the upper right border of the structure is curved, and is linked to the node that is in contact with the foundation by a group of thin bars.

It ust be noticed that, due to the less stringent boundary conditions adopted for the problems presented in this section, the stiffness matrix is singular. In fact, the global stiffness matrix always has one redundant row, but this row is usually eliminated when the boundary conditions are applied. Since this is not the case here, the problems above are harder to solve than usual truss topology problems. Even though, the structures shown in Figures 3, 4, 6 and 8 are quite reasonable.

6 Conclusions

The IR approach for solving bilevel problems allows for a lot of freedom in the choice of the method used to solve the lower level optimization problem. For particular problems, some special nonlinear programming algorithms may be devised. In this paper, we illustrate this fact dealing with a simple truss topology optimization problem.

We do not claim that the lower level strategy presented in Section 4 is better than others. Another strategy can be used in this framework. It is worth to mention that the feasibility of the iterates can be explicitly controlled at all iterations. We do not need to solve for feasibility when we are far from optimality. On the other hand, when we judge that we are sufficiently near to the optimal value, we may solve the problems in Step 1 with stronger stopping criteria by taking smaller values of the parameter r. This may be important in practice if these constraints need to be satisfied at an approximate solution.

With our approach, we were able to solve truss topology optimization problems with singular stiffness matrices and frictionless contact conditions. Besides, it would be easy to include other constraints on the bar volumes, on the displacements or even on the stresses. It would also be possible to consider other bilevel formulations of the truss topology problem, such as the one presented in [14]. We plan to extend this approach for truss problems with nonlinear elastic materials, that include nonlinear stiffness matrices. At least theoretically, these problems will pose no extra difficulty. The implementation, of course, has to be done very carefully.

Received: 15/VIII/10. Accepted: 05/I/11

#CAM-318/10.

This work was supported by PRONEX-Optimization (PRONEX-CNPq/FAPERJ E-26/171.164/2003-APQ1), FAPESP (Grant 06/53768-0) and CNPq.

• [1] W. Achtziger, On simultaneous optimization of truss geometry and topology. Structural and Multidisciplinary Optimization, 33 (2007), 285-304.
• [2] R. Andreani, S.L.C. Castro, J. Chela, A. Friedlander and S.A. Santos, An inexact-restoration method for nonlinear bilevel programming problems. Computational Optimization and Applications, 43 (2009), 307-328.
• [3] A. Ben-Tal and M. P. Bends∅e, A new method for optimal truss topology design. SIAM Journal on Optimization, 3(2) (1993), 322-358.
• [4] M. P. Bends∅e, A. Ben-Tal and J. Zowe, Optimization methods for truss geometry and topology design. Structural Optimization, 7 (1994), 141-159.
• [5] S. Dempe, Foundations of Bilevel Programming. Nonconvex Optimization and its Applications Series, 61, Kluwer Academic Publishers, The Netherlands (2002).
• [6] S. Dempe, Annotated bibliography on bilevel programming and mathematical problems with equilibrium constraints. Optimization, 52 (2003), 333-359.
• [7] S. Dempe, Bilevel programming - A survey. Technical Report TU 2003-11, Bergakademie Freiberg (2003).
• [8] A. Fischer and A. Friedlander, A new line search inexact restoration approach for nonlinear programming. Computational Optimization and Applications, 46 (2010), 333-346.
• [9] C.C. Gonzaga, E.W. Karas and M. Vanti, A globally convergent filter method for nonlinear programming. SIAM Journal on Optimization, 14 (2003), 646-669.
• [10] G. Haeser, Condições sequenciais de otimalidade. Doctoral Thesis, Unicamp 2009. Available for download in www.sjc.unifesp.br/docente/haeser/
• [11] N. Kikuchi and J. T. Oden, Contact problems in elasticity- A studyofvariational inequalities and finite element methods. Studies in Applied Mathematics, SIAM, Philadelphia, USA (1988).
• [12] A. Klarbing, Quadratic progras in frictionless contact probles. International Journal of Engineering and Science, 24 (1986), 1207-1217.
• [13] A. Klarbing, J. Petersson andM. Rönnquist, Truss topology including unilateral contact. Journal of Optimization Theory and Applications, 87(1)(1995), 1-31.
• [14] M. Kocvara and J.V. Outrata, Effective reformulations of the truss topology design problem. Optimization and Engineering, 7 (2006), 201-219.
• [15] J.M. Martinez, Inexact-restoration method with Lagrangian tangent decrease and new merit function for nonlinear programming. Journal of Optimization Theory and Applications, 111(1) (2001), 39-58.

Publication Dates

• Publication in this collection
  22 Mar 2011
• Date of issue
  2011

History

• Accepted
  05 Jan 2011
• Received
  15 Aug 2010