INVERSE TRUSS DESIGN AS A CONIC MATHEMATICAL PROGRAM WITH EQUILIBRIUM CONSTRAINTS

. We formulate an inverse optimal design problem as a Mathematical Programming problem with Equilibrium Constraints (MPEC). The equi- librium constraints are in the form of a second-order conic optimization problem. Using the so-called Implicit Programming technique, we reformulate the bilevel optimization problem as a single-level nonsmooth nonconvex problem. The major part of the article is devoted to the computation of a subgradient of the resulting composite objective function. The article is concluded by numerical examples demonstrating, for the ﬁrst time, that the Implicit Programming technique can be eﬃciently used in the numerical solution of MPECs with conic constraints on the lower level.


Introduction.
A truss is an assemblage of m pin-jointed uniform straight bars. The bars are subjected to only axial tension and compression when the truss is loaded at the joints; the load is denoted by f ∈ R n . With a given load and a given set of joints at which the truss is fixed, the goal of the designer is to find a truss that is as light as possible and satisfies the equilibrium conditions. In the simplest, yet meaningful, approach, the number of the joints (nodes) and their position are kept fixed. The design variables are the bar cross-sectional areas assembled in a vector a ∈ R m and the only constraints are the equilibrium equation and an upper bound on the weighted sum of the displacements of loaded nodes, so-called compliance. This model (or its equivalent reformulations) has been extensively analyzed in the mathematical and engineering literature (see, e.g., [1,2,17] and the references therein).
Assume that we solved, in the past, the truss optimization problem for some load vector f * ; the solution (bar cross-sectional areas) is denoted by a * . We are now facing the question of solving the inverse problem: given an optimal design a * ∈ R m , we want to find a load vector f ∈ R n such that a * is a solution of the truss optimization problem for this f ∈ R n . This might be useful, for instance, for existing constructions, to determine the loading (position and magnitude of forces) under which the structure performs best.
The resulting problem is a typical Mathematical Program with Equilibrium Constraints (MPEC) where on the upper level we minimize, with respect to the load vector f , the distance of the solution a from a given a * and on the lower level we solve the truss optimization problem for a given f . One of the established techniques for the solution of a class of MPECs is the so-called implicit programming (ImP) technique; see [22,28]. This technique replaces the initial bi-level problem by a single-level one, by using the solution map implicitly defined via the lowerlevel problem. The resulting single-level problem (in variable f only) is typically nonsmooth and one has to use a nonsmooth optimization algorithm for its solution.
For the implicit programming technique to work, some assumptions have to be satisfied. In particular, the lower-level optimization problem (in our case the truss design problem) must satisfy a certain constraint qualification which enables us to work with respective optimality conditions. Further, its solution map (assigning a to f ) must be at least locally unique and locally Lipschitzian. In case of truss design, the problem is usually formulated as a nonlinear nonconvex optimization problem that, unfortunately, does not satisfy the required constraint qualifications. The problem can be reformulated as a convex quadratically constrained quadratic problem which, however, does not include a as a (primal) variable. We thus use yet another reformulation of the truss design problem as a conic optimization problem with a Cartesian product of Lorentz cones. We will show that this problem (in variable a) satisfies the required constraint qualification and arrive in this way at a mathematical program with constraints in the form of a conic optimization problem, or conic MPEC. Under the assumption that its solution of the lower level conic problem is unique we will apply the ImP technique. The conic MPECs have been rarely studied in the literature, due to their technical complexity. The main difficulty lies in finding a subgradient of the (nonsmooth) composite objective that includes the solution map to the conic optimization problem. For that we will use recent results based on the generalized differentiable calculus of B. Mordukhovich ([23, 25]) and on the coderivatives of the projections onto the Lorentz cone ( [30]).
To the best of our knowledge, we present one of the first approaches to the solution of a conic MPEC that leads to a computationally approachable problem. In standard (not conic) MPECs one typically replaces the lower-level equilibrium problem by the corresponding KKT conditions (if possible), in which the (difficult) complementarity condition is either relaxed or added to the upper-level objective via a suitable (exact) penalty. In this way one eventually arrives at a nonlinear optimization problem in all primal variables and, additionally, the multipliers; see [16,32]. Alternatively, under additional assumptions, the ImP approach mentioned above can be used; see [19] for the comparison of these two techniques. The advantage of the first technique is the fact that one does not need the often complicated formula for the subdifferential of the composite function, as in the ImP technique. On the other hand, the resulting nonlinear optimization problem may be difficult to solve numerically because it is nonconvex, despite the original lower-level problem being convex or even linear. Moreover, to guarantee theoretical properties of this technique, one often needs similar assumptions as in ImP. The disadvantages of the first technique are magnified when considering a linear conic optimization problem on the lower level, as in this paper. While avoiding the complex analysis and technical formulas presented in Section 4, it leads to a nonlinear nonconvex conic optimization problem, for which there is no available software so far. Compared to that, the ImP approach used here requires solutions of linear second-order conic problems and of a nonconvex nonsmooth optimization problem in the control variable only. Numerical results presented in Section 5 show that, indeed, with this approach we can successfully find local (and often even global) optima in the considered MPECs. The reader is kindly asked to look at the practical relevance of our application from this point of view, as a demonstration of solvability of conic MPECs.
The plan of the paper is as follows. After introducing necessary notation, in Section 3 we define the truss design problem, give its various reformulations and show that the final formulation as a convex conic optimization problem could amenable for the solution by the ImP technique. In this section we also define our main problem, the inverse truss design, as a conic MPEC. Section 4 presents details of the implicit programming technique, in particular, the computation of the subgradient information for the composite objective function. In Section 5 we show results of two numerical examples.
2. Notation and preliminaries. For a Lipschitz continuous function f : R n → R, ∂f (x) denotes the Clarke subdifferential at x, defined bȳ where Ω f is the set of points at which f is differentiable.
In Section 4 we use the following notions of variational analysis. Given a closed set A ⊂ R n and a pointx ∈ A, we denote by T A (x) the Bouligand tangent cone to A atx, formed by vectors h ∈ R n such that there exist sequences h k → h and t k ↓ 0 for whichx + t k h k ∈ A for all k.
We further denote by N A (x) the Fréchet (regular) normal cone to A atx, defined as the (negative) polar cone to T A (x), i.e., by The limiting (Mordukhovich) normal cone to A at x, denoted N A (x), is defined by where "Lim sup" is the Kuratowski-Painlevé outer limit of sets (see [31]). If A is convex, then N A (x) = N A (x) amounts to the classic normal cone in the sense of convex analysis. We say thatx is a (normally) regular point of A, provided N A (x) = N A (x).
Let K ∈ R m be a closed convex cone with non-empty interior. For a ∈ R m we will use the notation a K 0 for a ∈ K and a K 0 for a ∈ int (K). Consider the following linear conic optimization problem with equality constraints: subject to Definition 2.1. We say that problem (1) satisfies the Conic Mangasarian-Fromovitz constraint qualification (C-MFCQ) if A has linearly independent rows and there exists a vector ξ ∈ R n such that C-MFCQ is equivalent to the Robinson constraint qualification for a general optimization problem; see, e.g., [5,15]. In the absence of equality constraints, it reduces to the Slater constraint qualification.
Definition 2.2. The cone L n ∈ R n defined by L n = (x 1 , . . . , x n−1 , x n ) x 2 1 + · · · + x 2 n−1 ≤ x n is called the n-dimensional Lorentz cone or the second-order cone.
Linear conic problems (1) with the K being the Lorentz cone are called secondorder conic programming problems (SOCP). The rest of our notation is standard; in particular, gph f denotes the graph of a map f , Bd (Ω) is the boundary of a set Ω, and Proj K (·) denotes the metric projection onto a convex set K ⊂ R n . Finally, for x ∈ R, x + := max{0, x}.
3. Optimum truss design. By truss we understand an assemblage of pin-jointed uniform straight bars in two or three spatial dimensions. The bars can only carry axial tension and compression when the truss is loaded at the joints. A truss is determined by positions of nodes and cross-section areas of bars, i.e., by vectors y ∈ Rñ and a ∈ R m . The bar lengths are denoted by i . The material properties of bars are characterized by their Young's moduli E i . The response of the truss to the vector of nodal forces f ∈ Rñ is measured by the vector of nodal displacements u ∈ Rñ. Some of the nodes are assumed to be fixed, that is, some components of u are forced to be zero; the simplest way how to treat these components is to omit them and work with reduced vectors u ∈ R n and f ∈ R n , n <ñ. We will always assume that the nodal coordinates y are fixed.
We now denote by the axial force in i-th bar and introduce bar stiffness matrices and assemble them in the global stiffness matrix of the truss and introduce the equilibrium equation Here r i denotes the vectors of directional cosines of the i-th bar; see [1]. To simplify We first recall the standard truss topology design problem. The design variables are the bar areas a i and we try to find such a design that the truss of limited volume is as stiff as possible, with respect to the given load f . This is formulated as the following standard nonlinear program (NLP), the so-called minimum compliance problem Proof. This follows directly from the KKT conditions of (3).
The next theorem shows that the seemingly straightforward optimization problem (3) does not satisfy (standard) Mangasarian-Fromovitz constraint qualification (MFCQ).
Theorem 3.2. Let (a, u) be feasible for problem (3), and K(a) be singular and such that some of its rows and columns are zero. Then MFCQ is violated at (a, u).
Proof. MFCQ is equivalent to the fact that the Lagrangian multipliers associated with the constraints of the optimization problem are bounded [5]. We will show that the multipliers in problem (3) may be unbounded.
We claim that MFCQ is violated when some components of a are equal to zero. Hence, without loss of generality, we will assume that a = ∞, so that the first constraint reduces to a i ≥ 0. Let µ ∈ R n , ν ∈ R m , λ ∈ R be the Lagrangian multipliers associated with the constraints in (3). The Lagrangian of (3) is then and the Karush-Kuhn-Tucker (KKT) conditions can be written as The first condition shows that µ solves the equilibrium equation for optimal a. However, we have assumed that K(a) is singular, in particular, that some of the rows and columns of K(a) are zero (see the example below the proof). That is, after a proper re-ordering, the stiffness matrix is of the form K(a) = K 1:n ,1:n 0 n ,n 0 n ,n 0 n ,n with some 0 < n , n < n such that n + n = n. The last n components of vector µ may be then arbitrary and thus unbounded. Even unbounded, vector µ must still satisfy the second KKT condition for all bars i = 1, . . . , m. Recall that K(a) = m i=1 a i K i and assume that the last n components of vector µ are indeed unbounded. Then, for every i = 1, . . . , m, either (K i ) n +1:n,n +1:n is a zero matrix and then µ K i u is finite, or (K i ) n +1:n,n +1:n is not zero. In the latter case, µ K i u is unbounded; however, the corresponding a i must be then equal to zero (as (K(a)) n +1:n,n +1:n = 0 and the diagonal elements of K i are nonnegative), the complementarity condition a i ν i = 0 is trivially satisfied and thus also ν i can be (nonnegative) unbounded for the second KKT condition to hold (see Example 1). Figure 1 with V = 1. The element (2) (1) stiffness matrices of the two bars and the right-hand side vector are

Example 1. Consider the truss shown in
Consider the point The global stiffness matrix is then K(a) = 1 0 0 0 and the equilibrium equation is satisfied by any µ of the form µ = (1, The second KKT condition reads then as Clearly, λ = 1. Now, the latter equation reduces to Without loss of generality, we may assume that u 2 < 0 (this component is arbitrary). The equation then remains valid when ν 2 and µ 2 go both to plus infinity with the same rate.
3.1. Dual problem (QCQP and QP formulations). We now introduce a reformulation of (3) that has a great impact on the numerical solution of truss design problems. It is the following quadratically constrained quadratic program (QCQP): Then next theorem is a compilation of several theorems from [1].
). The NLP and QCQP problems (3) and (4) are equivalent in the following sense: (i) If one problem has a solution, then also the other problem has a solution and min(3) = min(4) .
(ii) Let (u * , α * , ρ * ) be a solution to (4). Let further τ * be the vector of Lagrangian multipliers for the inequality constraints associated with this solution. Then (u * , τ * ) is a solution of (3). (iii) Let (u * , a * ) be a solution of (3). Let further r * be the Lagrangian multipliers associated with the upper bounds on a, respectively, and let α * be the multiplier for the volume constraint. Then (u * , α * , r * ) is a solution of (4).
The QCQP problem (4) can further be rewritten as a standard (linearly constrained) quadratic programming problem. We will use the fact that α is always nonnegative. This follows from the quadratic constraints in (4) and the fact that, by Assumption 2, ρ i = 0 for at least one i.
Proof. We first apply transformation of variables α, ρ in (4) Remark 1. Matrix C in Theorem 3.4 may be indefinite, in general. Indeed, if V is not big enough, relative to a i i , it may have one negative eigenvalue and m non-negative ones. Only when V ≥ m i=1 a i i (and this is not a typical case and contradicts Assumption 2), C will become positive semidefinite, according to the Gershgorin theorem [14]. It was shown in [1] that problem (5) is convex nonetheless, as C is positive semidefinite on the feasible set of (5). However, the consequence is that (5) cannot be readily solved by available software for convex quadratic programming, such as Gurobi [11].
3.2. Primal SOCP formulation. We recall a simple but useful lemma that shows the relation between convex quadratic constraints and second-order conic constraints. The proof follows from direct evaluation of the expression on the right-hand side.
Lemma 3.5. Let x ∈ R n , t ∈ R and s ∈ R, s > 0. Then Using this lemma, QCQP (4) can be immediately re-written as an SOCP problem 3.3. Dual SOCP formulation. The primal SOCP formulation (6) does not include the variable a, the cross-section areas. This will be desirable in the next section, so we have to work with its dual.
Proposition 3.6. The dual problem to (6) can be written as where a are the bar areas, q the bar axial forces and the objective function is equal to the compliance.
An important property of the above reformulation is stated in the next proposition.
Proof. By Assumption 1, the gradients of the equality constraints are linearly independent. We need to find a point (a, τ, q) feasible with respect to the equality constraints and strictly feasible with respect to the conic inequalities. The triple (ã,τ ,q) is the vector satisfying Definition 2.1.

3.4.
Uniqueness of solution to SOCP (7). It is well known that the optimal solution of problem (7), in particular, the design vector a may be non-unique. This is demonstrated in Figure 2: it shows a 3 × 3 truss with all nodes connected by potential bars. Both trusses presented in this figure and all their convex combinations are optimal solutions of problem (7) (this can easily be confirmed numerically). We have to emphasize, however, that this is a very rare situation requiring a special combination of the nodal positions and the load vector. Indeed, in the above example, a minor change in either the nodal positions or the force direction would result in unique solution a * of (7) . In the MPEC problem (8) below we want to identify some a * , an optimal solution of (7) for a given f * . Let S(·) : R n → R 3m be the solution map assigning a load vector f the (unique) optimal solution (a, τ, q) of (7). Denote by S a a "restriction" of this map satisfying gph S a = {(f, a) | ∃τ, q : (f, a, τ, q) ∈ gph S} . In the following, we will assume that the truss nodal positions and load vector f * are chosen in such a way that S a has a Lipschitzian single-valued localization around (f * , a * ), see [9, p.4]. In our case this means, by virtue of the convex valuedness of S a , that S a (f * ) = a * and there is a neighborhood U of f * such that S a is single-valued and Lipschitz continuous on U. It is readily seen from the proof of Proposition 3.6 and the definition of the axial force q that the single-valuedness and the Lipschitz continuity of S a over U amounts, in fact, the single-valuedness and Lipschitz continuity of S over U. This may be ensured, e.g., via the property of nondegeneracy introduced in the next section and the appropriate strong secondorder sufficient condition, see [6, Thm. 5].
3.5. MPEC for inverse truss design. Assume that we solved, in the past, the dual SOCP formulation (7) of the truss design problem for some load vector f * . We are now facing the question of solving the inverse problem: given an optimal design a * ∈ R m , we want to find a load vector f ∈ R n such that a * is a solution of (7) for this f ∈ R n . As mentioned in the Introduction, this may be useful, for instance, for existing constructions, when we want to determine the loading (positions and magnitudes of forces) under which the structure performs best.
As usual in such problems, we will try to find a vector a ∈ R m that minimizes the 2-norm of the difference between a and a * . The problem can be formulated as follows: subject to a solves (7) with load vector f f ∈ F , where F ⊂ R n is a feasible set for f . Problem (8) is a conic MPEC, as defined in the Introduction. As also mentioned there, one of the established techniques for the solution of MPECs is the so-called implicit programming technique (see [22,28]). The idea of the ImP approach is to replace problem (8) by the following nonsmooth optimization problem, referred to as ImP-MPEC: subject to (9) f ∈ F .
The new problem (9) can now be solved by algorithms and software of nonsmooth optimization. In particular, the ImP technique is known to couple well, e.g., with the bundle algorithm described in [7]. This algorithm requires at each iteration, i.e., for any choice of f ∈ F, the knowledge of the objective function value of (9) and of one arbitrary element of its Clarke subdifferential. While the objective function value is immediately available, the computation of an element of the subdifferential is highly nontrivial and, in case of conic MPEC, has not been conducted in the literature yet. The next section is devoted to this issue.

4.
Computation of subgradients of the composite objective. Assume that f belongs to a neighborhood of f * on which S : f → (a, τ, q) generated by the dual SOCP problem (7) is single-valued and Lipschitz continuous, as discussed in the previous section. In the computation of the "subgradient" information needed for the minimization of the respective composite objective function via the used bundle method, we take advantage of the generalized differential calculus of B. Mordukhovich and employ the chain rule [23, Thm. 1.110]. It follows that for an arbitrary continuously differentiable function ϕ of argument a and for a fixed quadruple (f , a, τ , q) ∈ gph S one has the inclusion whenever S is single-valued and Lipschitzian around f . So any element from the right-hand side of (10) provides a workable subgradient information at the considered point. Unfortunately, the currently available tools of variational analysis do not enable us to compute such an element at all points of gph S. In this section we will explain the used technique and discuss the properties of the resulting information which will be supplied to the bundle algorithm as a desired Clarke subgradient of the composite objective function in ImP-MPEC (9). To this purpose, we replace SOCP (7) by a slightly more general convex conic program which will allow us to simplify both the notation and the final formulas.
This conic program, dependent on parameter f ∈ R n , attains the form min y h(y) (11) subject to where y ∈ R o stands for the triple (a, τ, q), g : R n × R o → R r is affine linear and models the map in the equality constraints in (7), Ω ∈ R o is a closed convex set replacing the bound constraints on a and the constraint b(y) ∈ K generalizes the conic constraints in (7). We will assume that h : R o → R is convex and continuously differentiable, K ⊂ R p is a closed cone with vertex at 0 and non-empty interior and b : R o → R p is twice continuously differentiable and K-convex [5, Def. 2.103]. These assumptions are satisfied by the dual SOCP problem (7) and imply, in particular, that (11) is a convex program. Let Γ be the constraint multifunction in (11) defined by where Further, let Σ : R n ⇒ R o be the "solution map" of (11) which assigns each f the corresponding set of minimizers y. Clearly, by virtue of the assumed convexity, for f = f , y is a solution of (11) if and only if it is a solution to the GE 0 ∈ ∇h(y) + N Γ(f ) (y) , and, with Q(f, y) := N Γ(f ) (y), one has Define now the associated perturbation mapping M : According to [12,Thm. 4.1] one has that whenever M has the calmness property at (0, f , y); see [31, p. 399]. This property is fulfilled under the "standard" 2nd-order qualification condition which ensures even the stronger Aubin property of M around (0, f , y) ([31, Def.9.36]). In most applications of the ImP technique to various MPECs, the coderivative D * Q(f , y, −∇h(y)) can be computed and condition (15) can be verified. Even then, however, the subgradient information based on the right-hand side in (14) need not be always correct because it is an upper estimate of N gph Σ (f , y). In our case, the situation is even more complicated because, under the imposed C-MFCQ, we dispose only with an upper approximation of D * Q(f , y, −∇h(y)) derived in [25,Thm. 4.3]. Notice that this upper approximation requires the fulfillment of an additional calmness assumption and knowledge of all Lagrange multipliers associated with the constraint system (12). Fortunately, recent results in [10] indicate that certain conic programs fulfill generically the so-called nondegeneracy condition and under this condition D * Q(f , y, −∇h(y)) can be computed exactly. (More precisely, in [10] this has been proved for nonlinear semidefinite programming (SDP) problems and thus, as every SOCP problem can be formulated as an SDP problem, also for our dual SOCP (7).) To introduce the nondegeneracy in the case of problem (11), notice first that C-MFCQ at (f , y) for this problem amounts to the implication 0 ∈ (∇ y g(f , y)) λ + (∇b(y)) µ + N Ω (y) This follows from [5,Prop. 2.89] because both these conditions are equivalent with the metric regularity of the canonically perturbed conic constraint system around (y, 0). Following [5,Def. 4.70] and [29, p. 804], we may now say that y is a nondegenerate point of (17), provided (16) is strengthened to the form 0 ∈ (∇ y g(f , y)) λ + (∇b(y)) µ + span N Ω (y) It is easy to prove that if y is a nondegenerate point of (17), then there are unique Lagrange multipliers λ ∈ R r , µ ∈ N K (b(y)) and ν ∈ N Ω (y), satisfying the KKT conditions ∇ y L(f , y, λ, µ, ν) = 0 where L(f, u, λ, µ, ν) := ∇h(y) + ∇ y g(f, y) λ + ∇b(y) µ + ν is the Lagrangian associated with problem (11).
We may now invoke the recent results in [24,26] and arrive at the following statement.
Theorem 4.1. Let y be a nondegenerate point of (17) and λ ∈ R r , µ ∈ N K (b(y)) and ν ∈ N Ω (y) be the unique Lagrange multipliers satisfying (19). Then one has with some vectors d ∈ R r , v ∈ ker ∇ y g(f , y) .
Proof. The proof follows essentially from [26,Thm. 3.1], where this type of a secondorder chain rule has been derived under surjectivity of ∇ y Ψ(f , y). The replacement of this requirement by the nondegeneracy has been conducted in [24,Prop. 5.1]. The restriction v ∈ ker ∇ y g(f , y) comes from the fact that so that we have no restrictions on d, but have to require that v ∈ ker ∇ y g(f , y).
Vector a amounts to the single term ∇ f g(f , y) d because ∇ y g(f , y) does not depend on f .
Under the assumptions of Theorem 4.1 we thus immediately obtain that the right-hand side of (14) amounts to with some vectors d ∈ R r , v ∈ ker ∇ y g(f , y) .

(21)
It now follows from the definition of the coderivative that D * Σ(f , y)(y * ) can be approximated by the set of vectors ξ such that (ξ, −y * ) belongs to the set (21). This type of approximation has been used in numerous applications of the ImP technique to MPECs with complicated variational systems on the lower level [3,4,18].
So, given a vector η ∈ R o , let us define the "adjoint system" in variables d ∈ R r , v ∈ R p : η ∈ ∇ y L(y, λ, µ, ν) v + (∇ y g(f , y)) d On the basis of (21) we propose the following three-step procedure for the computation of a point that will be supplied as an element from ∂(h • Σ)(f ).
Indeed, following [31, Thm.6.14] and [23,Thm.4.1], it can be shown that N gph Σ (f , y) contains the set (21), provided the limiting coderivatives D * N K (b(y), µ) and D * N Ω (y, ν) are replaced by their regular counterparts. This fact, together with inclusion (14) produces a chain of inclusions which, under condition (c), become equalities. The nature of conditions (a)-(c) indicates that in the course of the iteration process we meet almost always points where our methods leads to a correct subgradient. This behaviour has been observed also in the case of the MPEC problem (8). Further notice that the algorithm of the used nonsmooth code BTNCLC is sufficiently robust so that it can cope with occasionally incorrect subgradients. Concerning (a),(b), one could attempt to verify implications (15) and (18) at the considered points. This would, however, be extremely time-consuming and definitely not practicable during the iteration process. It is substantially easier to verify the stationarity of the obtained results on the basis of a suitable stationarity concept. Finally notice that in many of our numerical examples presented in the next section, BTNCLC converged to points giving optimal objective values smaller than 10 −7 . Due to the nature of problem, without any further verification we may conclude that these points are (approximations of) global minima.
We conclude this section with a specification of the above procedure to our MPEC problem (8). With respect to the structure of K(= m i=1 L 3 ), it is convenient to put y = (y 1 , y 2 , . . . , y m ) ∈ R 3m with y i = (a i , τ i , q i ) ∈ R 3 , i = 1, . . . , m .
Putting everything together, we arrive at a system of linear equations in 7m + n + 1 variables α = (α 1 together with the relations System (24)- (27) represents a specific form of adjoint system (22). By a suitable choice of selections from the coderivative mappings D * N L 3 (Ay i , µ i ) and D * N [0, ai] (a i , ν i ), we obtain a square linear system in (α, d, γ, v); see [19, pp. 138-139]. The limiting coderivatives of N L 3 and N [0, ai] are well-known (see, e.g., [30,8] and the Appendix) and so the choice of suitable vectors α i , γ i (dependent on the position of considered points in the graph of N L 3 and N [0, ai] , respectively) does not make any difficulties. In step (iii) of Algorithm 1 we then put ξ = −d 2 , where (α, d, γ, v) is a solution of system (24)- (27).
Finally notice that under assumptions (a)-(c) stated above, the system (24)-(27) possesses a solution. This follows from the fact that the set of all vectors −d 2 arising in the solutions of (24)- (27) amounts under (a)-(c) to the limiting (Mordukhovich) subdifferential of ϕ•S at f , where ϕ is the objective in the ImP-MPEC problem (9). Since ϕ • S is Lipschitz around f , the claim is implied by the non-emptiness of this subdifferential, see [23, p. 86].

Numerical experiments.
In this section we present results of numerical examples which, on the one hand, demonstrate the correctness of our formulas and, on the other hand, the difficulties connected with this problem.
Intentionally, we do not specify realistic physical properties, like material elastic modulus E or density ρ. In all examples, they are assumed to be equal to one. This can be done without any loss of generality or applicability due to the linear dependence of the volume, compliance, and stiffness matrices on a, E and ρ. Switching from our default values to realistic physical values (say, E = 2.1 · 10 11 Pa and ρ = 7.8 · 10 3 kg/m 3 ) is a matter of a simple linear scaling of our results. Further, when we speak of an i × j truss, we have in mind a regular equidistant grid of i times j nodes, i in the horizontal direction and j in the vertical direction. Thus the dimensions of the nodal grid are (i − 1)(j − 1).
The truss topology software together with formulas for the subgradients were implemented in MATLAB. The SOCP problems (7) were defined using YALMIP [21] and solved by the conic solver MOSEK [27]. The nonsmooth optimization problem (9) was solved by the Bundle-Trust code BTNCLC [33] for nonsmooth nonconvex optimization with linear constraints. Stopping criterion for BTNCLC was set to ε = 10 −4 .
The examples were solved on a notebook with Intel Core i7 CPU M620 2.67GHz running 64-bit Windows 7.
In all examples, we have first solved the dual SOCP problem (7) with a given load vector f * to get an optimal design a * . Then, by solving MPEC (8), we tried to reconstruct this design, starting with a random load vector. We did not apply any constraints on f , i.e., we set F = R n in the ImP-MPEC (9). Notice that, using Lemma 3.1, the optimal solution a of the NLP problem (3) (and thus SOCP (7)) is invariant to the scaling of the load vector f . This means that problem (8), without any additional constraints on f , has always infinitely many optimal solutions giving the same value of the objective function.
The ImP-MPEC optimization problem (9) is nonconvex, hence BTNCLC may converge to a local minimum or just to a stationary point. The solution obtained (whether local or global) and the iteration process depends significantly on the initial iteration. We have adopted three strategies for choosing the initial point f ini : INI1: random vector of sum one; INI2: random vector of zero sum with norm 0.1 f * added to the original load f * ; INI3: random vector of sum one applied only at nodes with nonzero components of the original load f * ; here we assume that we know the loaded nodes but not the actual load vector.
Comparison with a derivative-free method. Is the complicated computation of the subgradient information (and the theory behind) really necessary? After all, with F = R n the implicit programming problem (9) is a nonsmooth unconstrained optimization problem. We may thus try to solve it by an appropriate derivativefree optimization (DFO) method. In this the examples below, we thus present also results obtained by a DFO algorithm and compare them with the BTNCLC results. Our DFO method of choice was the universal Nelder-Mead "simplex" algorithm and its implementation by Higham [13] Table 1 summarizes the results obtained by BTNCLC and NMSMAX from our three initial points. We present the computed optimal objective value c * and, in the next four columns, number of iterations needed to reach the objective values 10 −3 -10 −6 . The last column gives the number of iterations needed to reach the   Table 2 presents the results in the same format as Table 1. We can see that, in case of INI3, BTNCLC converged to a local optimum with c * = 0.44; the corresponding design is shown in Figure 5-right. The derivative-free code had serious problems in this example and only obtained an acceptable solution with INI3, however, in much bigger number of iterations than BTNCLC. Table 2. Results of Example 3. The last columns show the number of iterations needed to obtain the value of the objective function smaller than 10 −3 -10 −6 and (the last column) to fulfill the stopping criterion. We reiterate that the load vector f leading to the same design vector a is not unique. Not only is the design a invariant with respect to scaling of f , it may differ substantially from the original given load f * . An example is shown in Figure 5-left which presents the optimal solution obtained by BTNCLC from initial point INI1. As seen in Table 2 (and indeed in Figure 5 itself), this load vector leads to the same design as the original load f * depicted in Figure 4. Appendix. To compute vectors α, γ satisfying (26), (27), we make use of the wellknown relationship a ∈ N K (b) ⇔ b = Proj K (a + b) , which is valid for any closed convex set K. This enables us to express the limiting coderivative of N L n (·) in terms of the limiting coderivative of Proj L n (·) and employ the results of [30]. Since the formulas used to compute an element of the limiting coderivative of the metric projection Proj L n (·) onto the second order cone L n are rather complicated, we decided to give them in full detail in this Appendix. For x = (x, x n ) withx ∈ R n−1 , x n ∈ R, we define the spectral decomposition