ON OPTIMUM DESIGN OF FRAME STRUCTURES

. Optimization of frame structures is formulated as a non-convex optimization problem, which is currently solved to local optimality. In this contribution, we investigate four optimization approaches: (i) general non-linear optimization, (ii) optimality criteria method, (iii) non-linear semideﬁ-nite programming, and (iv) polynomial optimization. We show that polynomial optimization solves the frame structure optimization to global optimality by building the (moment-sums-of-squares) hierarchy of convex linear semideﬁnite programming problems, and it also provides guaranteed lower and upper bounds on optimal design. Finally, we solve three sample optimization problems and conclude that the local optimization approaches may indeed converge to local optima, without any solution quality measure, or even to infeasible points. These issues are readily overcome by using polynomial optimization, which exhibits a ﬁnite convergence, at the prize of higher computational demands


polynomial op
imization, which exhibits a finite convergence, at the prize of higher computational demands.

Introduction

Designing economical, efficient, and sustainable structures represents a major challenge of the contemporary society.While structural engineers literally explore designs that must satisfy the requirements of limit states analysis, there is usually an infinite number of such designs.Regardless of the properties of this feasible design space, it is required to select only one, the best design.This design quality is measured by an objective function, which usually approximates the expenses of production.Greatly simplified, one such (most common) criterion considers maximization of structural stiffness (while the amount of available material is limited), o , equivalently, minimization of structural volume (while requiring a certain structural stiffness) [1].

Among the structural optimization problems the greatest progress has been achieved so-far in optimizing (the cross-sectional areas of) truss structures [2].This development can be attributed to convexity of the feasible design space [3,Sections 1.3.5,3.4.3,and 4.8], as the (axial) structural stiffness is an affine function of the cross-sectional areas.On the other hand, when the bending stiffness comes into a consideration, convexity does not hold in general.Quite surprisingly, these problems are also much less studied, and to our knowledge only local optimization algorithms have yet been developed [4,5].

In this contribution, we investigate the problem of optimum design of frame structures.In particular, we develop four different methods in Section 2: (i) general non-linear optimization solved by the interior-point method of fmincon, (ii) the first-order optimality criteria (OC) method [1,5], (iii) reformulation of the problem into a non-linear semidefinite program (NSDP) solved by PenLab [6], and (iv) a suitably modified polynomial optmization (PO) problem (iii) solved globally using polynomial optimization methods [7] and the Mosek [8] optimizer.We show that the latter PO approach generates guaranteed lower and upper bounds on the solution, providing a means of assessing the design quality.Finally, Sectio t of three sample optimization problems to compare the optimization approaches.


Solution techniques to frame optimization


Problem statement

Let us consider the problem of optimum design of frame structures composed a finite number of nodes, n n , and of admissible elements, n e , defining the socalled ground structure [2].The frame elements are attributed with non-negative, and thus possibly zero, given-shaped cross-sectional areas a.For simplicity,
PREPRINT K i ( i ) =            Eiai i 0 0 − Eiai i 0 0 12EiIi 3 i 6EiIi 2 i 0 − 12EiIi 3 i 6EiIi 2 i 4EiIi i i i 0 0 12EiIi 3 i − 6EiIi 2 i 4EiIi i           (3)
we assume, in the following text, that the moments of inertia of individual cross-secti gree two1 , i.e., I i = c I,i a2 i for some c I,i > 0, which includes all cross-sections with given aspect ratios of all their components.Note that in such case, each cross-section is fully determined by its area.Optimizing the cross-sectional areas a we search the maximum stiff structures within the available volume V of a linear-elastic material.

The structural stiffness is measured (inversely) by the compliance c, work done by external u = u T K(a)u,(1)
where u constitutes the generalized displacement vector, and K(a) is the symmetric positive semi-definite stiffness matrix-a polynomial function of a,
K(a) = ne i=1 K i (a i ),(2)
assembled from contributions of individual elements K i (a i ).The lower the compliance, the stiffer is the structure with respect to the external forces.

In this contribution, we use the element stiffness matrix of Euler-Bernoulli frame elements (3), with E i denoting the Young modulus.In (1), f (a) is the external force column vector-a linear function of a-to allow for self-weight, assembled fr )
In the following text, we assume that ∀a > 0 : K(a) 0, i.e., the structure is not a kinematic mechanism and the stiffness matrix is positive definite (which is denoted by " 0") for all positive cross-sectional areas.Since K(a) has therefore the full rank, we also ave f (a) ∈ Im (K(a)), where Im(•) is the image of •.

This optimization problem is formalized as
min a,u f (a) T u (5a) s.t. K(a)u = f (a), (5b) T a ≤ V , (5c)a ≥ 0,(5d)
with being the column vector of the frame elements lengths.Notice that in general, (5) constitutes a nonconvex non-linear optimization problem because of the bilinear objective function (5a) and the polynomial equilibrium equality (5b) with possibly

ingular K(a).

On th
other hand, the volume constraint (5c) and the cross-sectional areas non-negativity constraint (5d) are affine functions of the design variables, u and a, and are in turn convex.The optimization problem (5) can be solved to local optimality using standard numerical optimization techniques.In particular, we will solve this problem using the interior-point method implemented in the fmincon function of Matlab.


Optimality criteria

The inherent difficulty of singularity of K(a) in (5b) can be circumvented by assuming a small positive lower-bound on the cross-sectional areas [1], which is denoted by ε in this study.Consequently, the forme e and the displacement field u can be excluded from the design variables 2 .Notice that in the limit when ε → 0, these problems are equivalent, but the smaller ε the higher the condition number of K(a), and thus it is more difficult to solve (5b).Hence, we assume ε = 10 −6 in this st L(a, u, λ, µ, ν) = f (a) T u + λ T [f (a) − K(a)u] +µ T a − V + ν T and ν; and 1 denoting the vector of all ones.In addition to the primal feasibility (5b)-(5d), the Karush-Kuhn-Tucker conditions require feasibility of the dual and complementary PREPRINT vol.no./ V = 0, (7b) ν ≥ 0, (7c) ν T (ε1 − a) = 0,(7d)
together with the stationarity of the Lagrangian with respect to u, d a.However, because the equation K(a)λ = f (a) in ( 8) possesses a unique solution due to K(a) 0, we have λ = u.Using this observation, the necessary first-order optimality conditions read as
0 = ∂L ∂a i = 2u T ∂f (a) ∂a i − u T ∂K(a) ∂a i u + µ i − ν i . (9)
Conditions ( 9) and (7d) then imply that at the optimum, the frame elements with a i > ε have equal constant energy
µ = 1 i u T ∂K(a) ∂a i u − 2 i u T ∂f (a) ∂a i . (10)
Aiming to satisfy (10), optimality criteria methods [1] build update schemes which increase stif tive process, based on the value of µ.In each i eration, the relative change of the design variables is bounded by the move limit ζ, assumed as ζ = 0.2 in this paper.Consequently, the fix-point update scheme reads as (11) with the tuning parameter η = 0.3 and with b
a (k+1) i = max max (1 − ζ)a (k) i , ε , a (k) i b (k) i η(k) i = u T ∂K(a) ∂a i u − 2u T ∂f (a) ∂a i µ i . (12)
Clearly, if ∀i ∈ {1,

. ., n e } : b
(k) i = 1
, we reac
a local minimum as (10) is satisfied.

Combination of ( 11), ( 12) with (5c) allows us to write the (current) volume V = T a (k+1) as a continuous function of , strict decreasing occurs when εI < a. Consequently, the bisection algorithm is used to find µ such that the volume constraint is satisfied.


Nonlinear semidefinite programming

In this section, we describe another approach to eliminate the displacement field variables u from the optimization problem formulation.First, let us rewrite (5) as
min a,c c (13a) s.t. c − f (a) T K(a) † f (a) = 0, (13b) T a ≤ V , (13c) a ≥ 0,(13d)
where K(a) † denotes the Moore-Penrose pseudoinverse of > 0 : K(a) 0, the (possible) singularity of K(a) is caused exclusively by zero rows and columns belonging to the degrees of freedom without any attached finite element (or, equivalently, a i = 0 for all attached eleme allows us to partition the stiffness mat ix into a positive definite principal submatrix K(a) and zero blocks, so that the pseudo-inverse equals
K(a) † = K(a) −1 0 0 T 0 . (14)
Using the same partitioning, the force vector f (a) is split into f (a) T 0 T T , in which the term 0 T appears due to the o (K(a)).Consequently, we see that (13b) can be rewritten to
c − f (a) T K(a) −1 f (a) = 0,(15)
and ( 1 (5).From ( 14) we have K † (a) 0, so that c ≥ 0 based on (13b).Moreover, iff f (a) = 0, we have both f (a) T K(a) † f (a) > 0 and c > 0. Because we minimize c (13a) and f (a) T K(a) † f (a) is bounded from below, (13b) can be simplified ent lemma, e.g., [9,Theorem 16.1], we further need to show that
I − K(a)K(a) † f (a) ) is always satisfied because f (a) ∈ Im(K(a)), so that we can substitute f (a) by K(a)v, where v is a vector of coefficients of the linear combination.Then, (17) is equivalent to
K(a) − K(a)K(a) † K(a) v = 0,(18)
with the term in the square brackets always zero [9, Lemma 14.1], as K(a) is symmetric.Finally, application of the generalized Schur complement lemma to ( 16) and (18) provides us with min
a,c c (19a) s.t. c −f (a)

−f (a) T K(a) 0,(19b)T
≤ V , (19c) a ≥ 0,(19d)s.t. K(a) − sf (a)f (a) T 0, (20b) T a ≤ V , (20c) a ≥ 0,(20d)
which is a useful reformulation when f is constant, i.e., self-weight is not considered.

It shall be noted that due to the polynomial matrix inequalities (19b) and (20b) both the optimization problems are non-convex in general.They can still be solved efficiently (to local optimality) using, e.g., augmented Lagrangian methods [6,10].In this contribution, we solve these problems using the open-source PenLab optimizer [6].


Polynomial o timization

Having introduced three different local approaches to frame structure optimization, it is natural and expected that one asks for a global opt loit the act that although ( 19) is non-convex it is indeed a polynomial optimization (PO) problem at the same time, which allows us to employ modern PO techniques, namely the Lasserre hierarchy [7,11], successively building tighter and tighter convex outer semidefinite programming (SDP) approximations called relaxations [7,Corollary 4.3].

To develop a more efficient formulation suitable for PO, we first recognize that the design vari all bounded both from below and above:
0≤a i ≤ V i , ∀i ∈ {1, . . . , n e },(21a)0≤ c ≤ ĉ. (21b)
Whi sed by the nonnegativity of the cross-sectional areas (19d), the upper bounds arise from the volume constraint (19c): none of the structural elements can occupy larger volume than V .In the compliance case (21b), the lower bound3 is due to K(a) 0, recall the discussion in Section 2.3, and the upper bound is provided by an arbitrary iform cross-sectional ar as
a = V 1 T 1,(22)
which are used for the upper-bound computation i = f (a) T K(a) † f (a),(23)
where K(a) † reduces to K(a) −1 in the c caled compliance c sc instead of the original a and c, where the scaled variables are defined in the [−1, 1] domain.Therefore, we have
a i = V (a sc,i + 1) 2 i (24)
and
c = ĉ(c sc + 1) 2 . (25)
Rewriting ( 21) in terms of a sc,i and c sc as unit ball constraints
a 2 sc,i ≤ 1, ∀i ∈ {1, . . . , n e },(26a)c 2 sc ≤ 1,(26b)
we finally arrive at an equivalent formulation to (19), min asc,csc 0.5ĉ(c sc + 1) (27a)
s.t. 0.5ĉ(c sc + 1) −f (a sc ) T −f (a sc ) T K(a sc ) 0, (27b) 2 − n e − 1 T a sc ≥ 0, (27c) 1 − a 2 sc,i ≥ 0, (27d) 1 − c 2 sc ≥ 0.(27e)
In ( 27), the objective function (27a) an

also the constrai
ts (27b)-(27e) are all polynomial inequalities4 non-negative in the feasible region of (27), which is therefore a semi-algebraic set.

The optimization problem (27) can readily be solved by building the (Lasserre) hie archy of convex linear SDP relaxations, with a monotonously converging objective functions box, and the underlying SDP relaxations are solved by the state-of-the-art Mosek optimizer [8].


Solution process

To simplify the notation, let us now denote the objective function polynomial (27a) by p 0 , an ) by p 1 to p PMI (27b) and let x = c sc a sc,1 . . .a sc,ne T (28) be the vector of the design variables.Moreover, let
b r (x) = 1 x 1 x 2 . . . x ne+1 x 2 1 x 1 x 2 . . . x 1 x ne+1 x 2 2 x 2 x 3 . . . x 2 ne+1 . . . x r 1 . . . x r ne+1 T(29)
PREPRINT vol.no./ On optimum design of frame structures denote the polynomial space basis of the maximum degree r.The , we can express each of the polynomials p j , j ∈ {0, . . ., 3} as a linear combination
p j (x) = |br(x)| β=1 p α,j,β (x α ) β (30)
of monomials
x α = ne+1 m=1 x αm m , ne+1 m=1 α m ≤ ch element in x, and the vector p α,j contains coefficients of the linear combinations of the monomials.Notice tha ne+1 m=1 α m ≤ d j , where d j stands for the degree of the polynomial p j .A similar approach is also applied to the elements of P 4 [12].

Further, introducing y, with its components y β corresponding to a monomial in the basis b r (x), we build the Lasserre hierarchy of convex l

ear semidefinite programming r
laxations min For all our test cases, this convergence was always finite.

In (32), the matrix M r (y) is the moment matri of the r-th order, and M r−dj is the (r − d j )-order localization matrix associated with p j s, we refer the reader to [7,12].


Recognizing global optimality

There are (at least) two ways to recognize whether y * r , a r-th order relaxation solution to (32), is globally optimal.

The first (common) way is based a sufficient condition for global optimality of (32) [15,16],
rank M r (y * r ) = rank M r−d (y * r ), (33)
in which y * constitutes the optimal solution in the basis b r (c sc , a sc ), and rank M r (y * ) is less or equal to the number of these solutions [7].The values of the original design variables can be extracted using Cholesky or singular value decompositions [16].

For the optimization pro

e another suffic
ent condition for global optimality.Let ãr be the (optimized) values of the cross-sectional areas and let cr denote the compliance found in the r

h order relaxati
n.Indeed, cr ≤ c * [7].Moreover, because any feasible design a generates an upper bound ĉr , recall Eq. ( 23    Note that (35) is substantially simpler to check than (33).If (35) is not satisfied in the (current) relaxation order r, we have at least the quality measure for the current solution in hand.


Sample problems

In this section, we describe and solve three rather small-scaled structural optimization problems, and compare the four optimization approaches.


Cantilever beam

Consider a (simpl ) cantilever beam design problem, as shown in Fig. 1a.The beam is discretized into n e finite elements of equal lengths, each of them assigned a prismatic square cross-section.The beam is loaded with a skew force at its tip.Further, we assume (dimensionless) E = 1 and V = 0.1.It can be seen from Table 1 that the problem is

latively "easy" to solve
as all of the tested algorithms converge to the unique global optima (shown   (1)  0.3 1429.31443.20 0.103 0.082 0.000 0.029 0.000 0.000 0.121 0.095 0.069 PO (2)  5.2 959.32 959.31 0.070 0.186 0.000 0.043 0.000 0.098 0.000 0.064 0.000 in 1b-1e), proved by the PO approach (Figs.1h-1j).


PREPRINT

In fact, we include this optimization problem mainly to show that PO scales unambiguously worse with the problem dimension than the other methods, and that the optimality criteria (OC) method is the fastest one.

Using OC, we can optimize the discretizations of 150 or 300 elements in 0.4 and 0.8 seconds (Figs.1f and  1h).


10-beam frame structure

Second, we consider the topology optimization problem of designing circular cross-sections of the 1

beam structure shown in Fig. 2a.T
is structure is loaded by two moments of magnitudes 1 and 2, placed at the nodes 2 and 3 , respectively.Note that there is no intersection between the elements 3 and 4 , and between 6 and 7 .We assume the Young modulus E = 1 and the volume bound V = 0.5.

Evaluating the four optimization methods, Table 2, we see that while the local optimization methods converged to nearby local minima of topologies shown in Fig. 2b, PO found the global optimum of a clearly different topology, Fig. 2c.Moreover, the global optimum possesses the compliance lower by 8%, which was reached in the second relaxation, Fig. 2d.


I-shaped girder with self-weight

Last, we consider a simply-supported I-shaped girder beam of the span 20, loaded by a uniform load 1 and self-weight.Due to the problem symmetry, we consider only one half of the problem, Fig. 3a, and discretize it using 5 finite elements of equal length.This (half of the) girder beam is allowed to utilize at most V = 0.2 material of the Young modulus E = 10 4 and of density ρ = 3.The beam cross-sections are parameterized by the parameter t p , which denotes the thicknesses of the flang

and web.The
cross-sectional area equals a(t p ) = 18t 2 p .We assume that the upper surfaces of the top flanges are aligned, so that we have I(t p ) = 696t 4 p = 58/27a(t p ) 2 .Consequently, we optimize a i of individual elements, rather than t p .

Table 3 reveals that for this specific problem the OC converged to the global optimum.On the other hand, the fmincon and non-linear semidefinite programming approaches converged to an infeasible point.Polynomial optimization exhibited convergence to the global optimum in three relaxations, Fig. 3c.Notice, moreover, that the designs found for r = {1, 2} were of very high qualities.


Conclusions

In thi contribution, we investigated four topology optimiz