Reliability Engineering and System Safety 72 (2001) 165±177
www.elsevier.com/locate/ress
Optimal structural design under stochastic uncertainty by stochastic linear programming methods K. Marti* Aero-Space Engineering and Technology, Federal Armed Forces University Munich, D-85577 Neubiberg/Munich, Germany Received 30 July 2000; accepted 27 December 2000
Abstract In the optimal plastic design of mechanical structures one has to minimize a certain cost function under the equilibrium equation, the yield condition and some additional simple constraints, like box constraints. A basic problem is that the model parameters and the external loads are random variables with a certain probability distribution. In order to get reliable/robust optimal designs with respect to random parameter variations, by using stochastic optimization methods, the original random structural optimization problem must be replaced by an appropriate deterministic substitute problem. Starting from the equilibrium equation and the yield condition, the problem can be described in the framework of stochastic (linear) programming problems with `complete ®xed recourse'. The main properties of this class of substitute problems are discussed, especially the `dual decomposition' data structure which enables the use of very ef®cient special purpose LP-solvers. q 2001 Elsevier Science Ltd. All rights reserved. Keywords: Optimal structural design; Random structural parameters; Stochastic linear programming
1. Optimal plastic design under stochastic uncertainty The limitations of elastic structural theory with respect to the assessment of actual structural resistance, hence of safety, is well known [1]. Indeed, to accept as valid the elastic behaviour up to the critical point of the mechanical structure means to identify the structural failure with ®rst attainment of a critical stress condition in any point of the structure, i.e. to assume a weakest-link type of behaviour. This holds in practice in some cases, such as statically determinant structures or structures made of brittle material, but in most structures a stress redistribution allows the loads to grow further beyond the level corresponding to the ®rst critical condition of elastic theory. Plastic collapse of the structure occurs [31,35] when the structure is converted into a mechanism by providing a suitable number and dispositions of plastic zones. Thus, according to limit analysis, collapse is identi®ed with the development of plastic hinges, i.e. zones of material with stress in the yield limit condition, in such a number and location to allow a movement of the whole or a part of the structure without requiring deformation of the zones with stresses below the yield limit (trans-
* Present address. Fak. LRT, UniBwMunchen, D-85577 Neubiberg/ Munchen, Germany. Tel.: 149-89-6004-2541; fax: 149-89-6004-4092. E-mail address:
[email protected] (K. Marti).
formation of the structure into a plastic collapse mechanism). Limit analysis is concerned [13,14,31] with establishing the strength of a structure, i.e. its capacity for the supporting of loads. Hence, using the plastic ductility of structural materials in improving the design of structures, limit analysis is not concerned with deformation: it cannot therefore provide the load carrying capacity for a structure with elements that have a limited ductility or deformability, nor for a structure which becomes unstable because of the displacements induced by plastic deformation, see Ref. [31]. For elastic-perfectly plastic materials, the ultimate load condition corresponding to complete collapse of the structure can be obtained by means of plastic limit (collapse) analysis through application of either a pair of dual theorems, i.e. the static (safe or lower bound) theorem and the kinematic (unsafe or upper bound) theorem. Considering a structure made of ductile material, collapse does not occur if there exists a statically admissible stress ®eld, i.e. a stress vector in equilibrium with the applied loads and not violating the strength inequality or yield condition at any point of the structure. Conversely, collapse occurs if a collapse mechanism exists such that the work done by the applied loads is larger than the internal plastic work. Since the fundamental static and kinematic relations for limit (collapse) analysis appear [16], naturally as linear equations, and the yield condition is convex, plastic limit
0951-8320/01/$ - see front matter q 2001 Elsevier Science Ltd. All rights reserved. PII: S 0951-832 0(01)00003-5
166
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
analysis and design problems may be formulated by a linear or a convex program. However, using linear approximates of the convex yield condition, plastic limit analysis and optimal design problems may be formulated always (approximatively) in the framework of linear programming (LP), see e.g. [3,15,36]. Having the powerful and most ef®cient LP tools, a remaining major dif®culty in plastic structural analysis and optimal plastic design as well as in other branches of optimal structural design, e.g. optimal material layout, is the fact that the material coef®cients, such as yield stresses in compression and tension, and the applied (external) loads, but also many other model parameters of the analysis or design problem, such as cost, weight coef®cients, etc., are not given, ®xed quantities in practice. Due to random variations of the material and the loadings, due to manufacturing and modelling errors and further random disturbances, these parameters must be modelled through random variables on a certain probability space, where the joint probability distribution of the random quantities can be assumed to be known. Thus, in plastic analysis and in optimal plastic design the given structural optimization, material optimization or topology optimization problem with random parameters must be replaced by an appropriate deterministic substitute problem using stochastic optimization methods, see e.g. [9,10,19,23]. Based on certain decision theoretical principles, substitute problems incorporate the random parameter variations in order to reach robust optimal designs, i.e. optimal designs which are insensitive with respect to random parameter variations. This goal can be achieved by minimizing the expected structural costs, e.g. the structural weight or volume, subject to constraints for the expected costs of structural failure; as a special case, here we have constraints for the probability of failure for the whole structure and/or certain parts of the structure. 1.1. Fundamental theorems from plastic limit analysis From collapse load analysis, the following fundamental theorems for structures having rigid-plastic or elasticperfectly plastic materials are well known [1,3,4,6,12,26]: Static Theorem (ST) (Lower bound or safe theorem). If any stress distribution throughout the structure can be found which is everywhere in equilibrium internally and balances the external loads, and at the same time does not violate the yield conditions, these external loads will be carried safely by the structure. By means of duality theory (Kuhn±Tucker theorem) the equivalence between the above (ST) and the following kinematic theorem can be established [3,6,27]. Kinematic Theorem (KT) (Upper bound or unsafe theorem). Collapse occurs if a collapse mechanism, i.e. a pair ~ exists, ful®lling the compatibility condition, of vectors
u; u
such that the work done by the external loads is larger than the corresponding internal plastic work. In limit analysis and in optimal plastic design of structures [5,7,25±28,32,35±37] we have therefore the following basic conditions (I), (II): (I) The equilibrium equation CF R
l0 R0 :
1
Here, C is the m £ n equilibrium matrix with rank C m , n; R; R0 are external load vectors, and l0 is a load factor. Moreover, F denotes the vector of internal or member loads (axial forces, twisting and bending moments). (II) The yield condition For trusses [11] or if, for plane or spatial frames, no member load interactions are taken into account [11,15], then the yield condition can be represented by the following inequalities FL
s L ; X # F # FU
s U ; X;
2
where the bounds F L, F U for the n-vector F of member loads depend on the r-vector X of design variables Xk ; k 1; ¼; r; and the vectors s L ; s U of yield stresses in compression (,0) and tension (.0). If for plane or spatial frames member load interactions must be taken into account, then for each beam i 1; ¼; B we have to consider the vectors of member loads [33] at the two (`negative', `positive') ends of the beam ! ! ti ti 2 1 Fi ; Fi
for plane frames
3a m2 m1 i i 1 1 0 0 ti ti C C B B B t C B t C C C B B i i C C B B F1 F2 i B 2 C; i B 1 C
for spatial frames
3b B mi C B mi C A A @ @ 2 1 mi mi containing the axial force ti and the bending moments m2 i ; m1 ; respectively, for plane frames, and the axial force t , the i i 2 1 twisting moment ti and two bending moments m2 i ; mi ; mi ; 1 mi ; respectively, for spatial frames. Let ! Npi
s y ; X ;
4a Fpi Fpi
s y ; X U Mpi
s y ; X 1 0 Npi
s y ; X C B B T
s ; X C C B pi y C;
4b Fpi Fpi
s y ; X U B C B y B Mpi
s y ; X C A @ z Mpi
s y ; X respectively, denote the corresponding principal axial, twisting and bending plastic capacities of the ith beam, where X designates again the design r-vector and s y is the vector of yield stresses s yi ; i 1; ¼; B: Then, the
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
yield condition reads 2 F21 pid Fi
[ Ki ;
1 F21 pid Fi
[ Ki ;
i 1; ¼; B;
5a
where Fpid denotes the diagonal matrix having the components of Fpi on its diagonal, F21 pid is the inverse of Fpid, and K1 ; K2 ; ¼; KB
5b
are the closed convex yield domains [8,15,32,35] in R2 or R4 related to the beams i 1; ¼; B; where the origin 0 of R2 ; R4 ; respectively, is an interior point of Ki. Denoting then by z pi
z U inf l . 0 : [ Ki ; z [ R2
R4 ;
6 l the positively homogeneous, subadditive distance or Minkowski functional [17,18] of Ki ; i 1; ¼; B; the yield conditions (5a) and (5b) can be represented also by the inequalities
pi
Fpid
s y ; X21 F2 i # 1; i 1; ¼; B;
7a
pi
Fpid
s y ; X21 F1 i # 1; i 1; ¼; B:
7b
Multiplying Eqs. (7a) and (7b) with Npi . 0; we ®nd the equivalent conditions 2 pi
r21 pid Fi # Npi
s y ; X; i 1; ¼; B;
8a
1 pi
r21 pid Fi # Npi
s y ; X; i 1; ¼; B;
8b
where the vector of ratios rpi is de®ned by ! Fpi Mpi 0 rpi U 1; ; Npi Npi ! y z Fpi Tpi Mpi Mpi 0 rpi U 1; ; ; ; Npi Npi Npi Npi
9a
9b
respectively. In many cases, e.g. for special classes of design vectors X, or in bilevel structural optimization procedures for frames [11], the ratios (8a) and (8b) are ®xed, hence,
10 rpi cpi ; i 1; ¼; B; with given ®xed vectors cpi, 1 # i # B; independent of s y and X. (II.1) Dual decomposition of the yield condition Using the dual representation [18] of pi ;
pi
z max{u 0 z : u [ Ui }; where the convex, compact sets Ui are de®ned by Ui U {u [ R2
R4 : u 0 z # Qi
z; for all z};
11b
then the distance functional pi of Ki can be approximated by
p~ i
z U
max u 0ij z; j[Ji
2
4
z [ R
R ;
position approach, the yield domains Ki ; i 1; ¼; B; are approximated by the convex yield polygons [15,35]. K~ i U {z [ R2
R4 : p~ i
z # 1} {z [ R2
R4 : u 0ij z # 1; j [ Ji }; i 1; ¼; B:
12a
where Ji is a ®nite index set
uJi u , 11; and uij, j [ Ji ; are certain elements of Ui. Thus, by means of this dual decom-
12b
Hence, from Eqs. (12a) and (12b) we obtain the following approximative yield condition: 2 p~ i
r21 pid Fi # Npi
s y ; X; i 1; ¼; B;
13a
1 p~ i
r21 pid Fi # Npi
s y ; X; i 1; ¼; B:
13b
Assuming now that condition (10) holds, then the approximative yield conditions (13a) and (13b) reads h 0pij F2 i # Npi
s y ; X;
j [ Ji ; i 1; ¼; B;
14a
h 0pij F1 i # Npi
s y ; X;
j [ Ji ; i 1; ¼; B;
14b
where the vector hpij ; j [ Ji ; i 1; ¼; B; is de®ned by hpij U c21 pid uij :
14c
Thus, in this case, the yield conditions (7a) and (7b) can be represented approximatively, cf. [15,25,34,35], by the system of linear inequalities HF # F0
s y ; X
PB
15a
forPF, where H is a 2 i1 uJi u £ n matrix and the 2 Bi1 uJi u -vector F0 F0
s y ; X is de®ned by F0
s y ; X U
Np1 ; ¼; Np1 ; Np1 ; ¼; Np1 ; ¼; NpB ; ¼; NpB ; NpB ; ¼; NpB 0
15b Obviously, the yield condition (2) can also be represented by a system of linear inequalities for F of the type Eqs. (15a) and (15b). Besides the equilibrium equation (1) and the yield conditions (2), (5a), (5b), (7a), (7b), (13a), (13b), (15a) and (15b), respectively, in optimal structural design there is still a primary structural (III) Weight, volume or more general cost function G0 G0
g; X
11a
167
16
for the evaluation of the design vector X, where g is a vector of weight, scale or cost coef®cients. A ®rst important example is the volume, weight or material cost function de®ned by G0
g; X U
B X i1
gi Li Ai
X;
17a
where Li is the length and Ai Ai
X denotes the crosssectional area of the ith element. Considering the bending capacities Mpi ; i 1; ¼; B; of a plane frame, we ®nd the
168
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
`cost' criterion [11,37]. G0
g; s y ; X U
B X i1
gi Mpi
s y ; X:
17b
In addition to the relations (I)±(III), for the design r-vector X we have some further simple deterministic constraints, e.g. box constraints, of the type X [ D0 U {X [ Rr : A0 X b0 ; X $ 0};
18
have to formulate now appropriate deterministic substitute problems for programs (20a)±(20e) and (21a)±(21e), taking into account the random variability of the total random parameter vector a(´). Hence, we have to apply stochastic optimization methods [9,10,19]. 2. Stochastic optimization 2.1. Derivation of deterministic substitute problems
where (A0,b0) is a given, ®xed matrix.
Due to random material variability, stochastic uncertainty of the external loadings, manufacturing and modelling errors, etc., the total parameter vector of model parameters and external loads
Violations of the yield conditions (20c) or (21c) cause damages or failure of the structure. Hence, costs of repair or reconstruction and often also further penalty costs are caused by violations of the yield condition [1,30]. These costs can be represented approximatively by the cost model
a U
s y ; s L ; s U ; Ei ; g; R; ¼;
DG DG
h; F; D;
1.2. Decision making under stochastic uncertainty
19a
where Ei is the elasticity modulus of the ith element, must be modelled by a random vector a a
v; v [
V; A; P;
19b
having a certain probability distribution P a
´ : Consequently, from (I)±(III) and Eq. (18) we obtain the following structural optimization problem with random parameters: min G0
g
v; X
20a
s.t. CF R
v
20b
fp
Fp
s y
v; X; F # 1
fp
rp
v; X; F # N
s y
v; X
20c A 0 X b0
20d
X $ 0:
20e
In the random program (20a)±(20e), the vectorial constraint (20c) represents the yield conditions (7a), (7b), (8a) and (8b), respectively. Using Ð under condition (10) Ð the approximative yield conditions (15a) and (15b), we have the random program min G0
g
v; X
21a
s.t. CF R
v
21b
HF # F0
s y
v; X
21c
A 0 X b0
21d
X $ 0:
21e
In order to obtain reliable/robust optimal designs X p, we
22a
where the discrepancy D in Eqs. (20c) and (21c), respectively is de®ned by fp
Fp
s y
v; X; F 1 D 1;
22b
fp
rp
v; X; F 1 D N
s y
v; X;
22c
HF 1 D F0
s y
v; X;
22d
respectively, and F must ful®l the equilibrium Eq. (1); moreover, the vector h contains certain cost, scale or weight factors. In many cases we have that
DG
h; F; D 0 for D $ 0;
CF R:
23
Inserting now the additional costs caused by damage or failure of the structure into the objective function or into the constraints of programs (20a)±(20e) and (21a)±(21e), different substitute problems occur. 2.2. Substitute problems based on failure probabilities Taking 0±1-cost functions DG DG
h; F; D; corresponding to structural survival/failure, in case of the cost models (22a) and (22d), the risk function, i.e. the expected cost function
D G
h
´; F
´; D
´ U EDG
h
v; F
v; D
v;
24
is equal to the failure probability pf pf
X 1 2 ps
X 1 2 P
there is an internal loading F such that CF R
v and HF # F0
s y
v; X:
25 Based on failure probabilities, the following reliabilityoriented substitute problems (RBO) are well known [20,21,30]:
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
Problem type A (expected cost minimization under reliability constraints) min EG0
g
v; X
26a s.t. pf
X # a
maximum failure probability
26b A 0 X b0
26c
X $ 0:
26d
Problem type B (reliability maximization under expected cost constraints) min pf
X
27a
s.t. EG0
g
v; X # Gmax
maximum budget 0
27b
A 0 X b0
27c
X $ 0:
27d
Using the pair of dual linear programs (LP) related to the survival conditions (21b) and (21c) de®ning the probability of survival ps, i.e. the so-called static/kinematic LP [3,6,20], the probability of failure pf 1 2 ps can be represented [20] by pf
X P
Ml
a
v; X , 0 for at least one index l; 1 # l # l0 ;
28a
where Ml Ml
a; X; l 1; ¼; l0 ;
28b
are the limit state functions of the structure. Remark 2.1. Since l0 is equal to the number of extreme points of the feasible domain of the kinematic LP, the integer l0 is very large in general. Solving substitute problems of type A or B, we have to carry out then the following calculations [1,20,21] which are very time-consuming in general: 1. Determine a certain collection of `most relevant' limit state functions Ml Ml
a; X; l [ L; where L , {1; 2; ¼; l0 } is a (small) index set. 2. Approximate pf
X < p~ f
X by using the above collection Ml ; l [ L; and by using certain probability bounds, e.g. Bonferroni inequalities. 3. Compute the approximates p~f
X numerically by means of FORM/SORM, RSM or sampling techniques, etc. 4. Approximate (if necessary) the expected cost function EG0 by means of Taylor expansion methods, RSM, sampling techniques, etc. 5. Finally, solve the resulting approximate nonlinear reliability-based optimization problem of type A or B by applying an appropriate parameter optimization proce-
169
dure, e.g. an SQP method. 2.3. Stochastic linear programming (SLP) techniques An alternative to the numerically very expensive method described brie¯y in Section 2.2 is obtained now by means of stochastic (linear) programming (SLP) techniques. For trusses, or after appropriate assumptions/approximations in more general cases, for the computation of an optimal design point X p we have, as described in the next section, an LP with a special structure according to Fig. 1, cf. [22]. 3. Optimal structural design by SLP In the following we suppose that the yield condition can be represented by Eqs. (15a) and (15b), hence, HF # F0
s y ; X: This holds for trusses, if no member load interactions are taken into account, for special classes D0 of design vectors, or if the ratios rpi ; cf. Eqs. (9a) and (9b), are kept ®xed and varied by steps only in a bilevel optimization procedure for the general case. Representing now the vector of discrepancies D in Eqs. (22a) and (22d) by
D
v D1
v 2 D2
v; D1
v; D2
v $ 0;
29a
and selecting the linear cost function
DG
q
v; F
v; D
v U q2
v 0 D2
v 1 q1
v 0 D1
v;
29b where h
v q
v U
q2
v; q1
v $ 0 is a further vector of cost coef®cients, for the minimization of the total expected costs E
G0 1 DG; we get the following stochastic optimization problem: min E
G0
g
v; X 1 q2
v 0 D2
v 1 q1
v 0 D1
v
30a s.t. CF
v R
v a:s:
with probability1
30b
HF
v 1 D1
v 2 D2
v F0
s y
v; X a:s:
30c
A 0 X b0
30d
D1
v $ 0;
X $ 0;
D2
v $ 0 a:s:
30e
De®ning now T
v; X U
WU
!
0 2F0
s y
v; X
0
0
C
I
2I
H
;
h
v U
R
v 0
!
31a
!
recourse matrix
31b
170
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
Fig. 1. LP with dual decomposition data structure.
0
D1
v
1
B 2 C C y
v U B @ D
v A;
yI
v U
F
v 0
q1
v
D1
v D2
v
area of the ith rod element of a truss, then in case of the 1st stage cost function (17a) we have that
! ;
c0
v U 0; ci
v U gi
vLi ; i 1; ¼; B;
31c
1
B 2 C C q
v U B @ q
v A ; 0 the above problems (30a)±(30e) can be represented also by 0
min E
G0
g
v; X 1 q
v y
v
32a
s.t. T
v; X 1 Wy
v h
v a:s:
32b
A 0 X b0
32c
X $ 0;
yI
v $ 0 a:s:
32d
However, this is a so-called `stochastic program with complete ®xed recourse' [9,10,19,23]. In the case that the function, mapping, respectively, G0
g
v; X U c0
v 1 c
v 0 X
33a
T
v; X U v0
v 1 T
vX
33b
is linear with respect to X, as e.g. for trusses [5,6,11], on special feasible domains D0, or after a linearization step [19], problem (32a)±(32d) is a `stochastic linear program with complete ®xed recourse' [9,10,23]. Remark 3.1.
If Xi Ai
X denotes the cross-sectional
34a
and for v0
v; T
v we ®nd, cf. Eqs. (2), (15a) and (15b), 0 1 0 B C U C v0
v U 0; T
v U B
34b @ 2s d
v A;
s dL
v where s L ; s U are the vectors of yield stresses in compression and tension, and s dL ; s dU denote the related diagonal matrices, cf. Eq. (5a). For stochastic programs with recourse there exists a large amount of theoretical results [9,10,19] and ef®cient numerical solution procedures [2,10,19,23] based on: ² discretization techniques (for Pa
´ and sequential linearization/convex approximation (with respect to X); ² special purpose LP-solvers: special purpose LP-solvers use extreme point methods (special factorization techniques), interior point methods in combination with factorization schemes, column splitting, decomposition approaches (L-shaped methods, decomposition with sampling, stochastic decomposition), Lagrangian-based approaches (Lagrangian dual ascent, augemented Lagrangian method). We still mention here that for the ef®cient computation of the equilibrium matrix C, recently also a special subroutine to the FEM software package `NASTRAN' is available [24].
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
171
In the following, several important properties of the stochastic program (32a)±(32d) are presented.
is the kinematic-type LP [1,20]
3.1. The second-stage linear program
s.t.
For a given design vector X, a so-called 1st stage decision vector, and a given realization
q
v; T
v; X; h
v of the random variables
q
´; T
´; X; h
´; the so-called secondstage linear program is de®ned by
C 0 u H 0 u~
38b
2q1
v # u~ # q2
v:
38c
max R
v 0 u 2 F0
s y
v; X 0 u~
38a
35a
By means of the optimal value functions Q(q,z) of the second-stage LP (35a)±(35c), the stochastic optimization problem (32a)±(32d) can be represented also by
Wy z
z h
v 2 T
v; X
35b
min E
G0
g
v; X 1 Q
q
v; h
v 2 T
v; X
yI $ 0;
35c
min q
v 0 y s.t.
let then
36 Q Q
q
v; z; z [ Rm ; PB m U m 1 2 i1 uJi u; denote the optimal value of the LP (35a)±(35c) for arbitrary right hand sides z [ Rm ; Q
q
v; z is also called `recourse function' [9,10]. Since q
v 0 y $ 0
37a
for feasible points y in Eqs. (35a)±(35c) {Wy : yI $ 0} Rm ;
37b
we have this lemma:
1. LP (35a)±(35c) has an optimal solution yp
D1p ; D2p ; Fp for each z [ Rm : 2. The optimal value function Q Q
q
v; z of LP (35a)± (35c) is a positively homogeneous, subadditive functional with respect to z. Proof.
s.t. A 0 X b0
39b
X $ 0:
39c
This yields the following result, cf. Eqs. (33a) and (33b): Theorem 3.1. If X ! G0
g
v; X is convex and X ! T
v; X is linear, then problem (39a)±(39c) is a convex program. Proof.
See Refs. [9,10].
3.1.1. Selection of the cost factors q 2, q 1 As already indicated by relation (23), for q 1 we may take q1 0:
Lemma 3.1.
See Refs. [10,18].
Obviously, for z U h
v 2 T
v; X; the second-stage LP of static-type yields the most favorable internal loading Fp Fp
q
v; h
v 2 T
v; X of the structure for given design vector X and data
q
v; T
v; ´; h
v: The dual of the static-type LP (35a)±(35c), having also the form
39a
40a
Moreover, if the 1st stage cost function G0
g
v; X is de®ned by Eq. (17a), then an appropriate choice for the components q2 ij ; j [ Ii ; i 1; ¼; B; reads q q2 ij
v U gi
v
Li ; s i
j [ Ii ; i 1; ¼; B;
40b
where s i s i
s yi ; Ei is a certain stress coef®cient, and gqi
v denotes again a scale factor. Since yI {D1 ; D2 } 2 contains force components, cf. Eq. (15b),
1=gqi q2 ij Dij are volume terms evaluating the member displacements of the structure. Working with the 1st stage cost function (17b) for plane frames, q2 ij may be selected by q q2 ij U gi cpi2 ;
j [ Ii ; i 1; ¼; B;
40c
35d
where cpi2 denotes the second component of the vector cpi ; see Eqs. (9a) and (10), and gqi is a scale factor. For more general cases, the cost vectors q 2, q 1 can be selected in similar way.
CF R
v
35e
3.2. Discretely distributed problems
HF 1 D1 2 D2 F0
s y
v; X
35f
Assume here that the total random parameter vector a(´) has a discrete probability distribution
D1 $ 0;
35g
P
a
v a j aj ; j 1; ¼; s;
min q1
v 0 D1 1 q2
v 0 D2 s.t.
D2 $ 0;
41a
172
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
having realizations a j taken with probabilities aj . 0; j 1; ¼; s: This holds if: ² there are only a ®nite number s of material, cost, loading,¼, scenarios a j ; j 1; ¼; s; which may occur with certain probabilities a1 ; ¼; as ; or ² Eq. (41a) results from a discretization of an original (continuous) distribution P
0 a
´ : Assuming that G0
g
v; X; T
v; X can be represented by Eqs. (33a) and (33b), as for trusses, see Remark 3.1, on special feasible domains D0, or after linearization with respect to X, let then denote T
j
; hvj
j
Uh 2
v0j ; q j ; j
1; ¼; s;
41b
the realizations of T
v; hv
v U h
v 2 v0
v; q
v; respectively, related to the realizations (scenarios) a j ; 1 # j # s; of the random parameter vector a
v: A central role in the various numerical solution techniques for the stochastic optimization problem (32a)±(32d) is played by the following result. Theorem 3.2. Assume that the linearity conditions (33a) and (33b) hold and the distribution assumptions (41a) and (41b) are ful®lled. Then the stochastic optimization problem (32a)±(32d) can be represented by a linear program with a dual decomposition data structure according to Fig. 1. Proof.
4. Basic solution techniques for stochastic programs 4.1. Stochastic linear programs having a discrete distribution
42
with a given vector q0. According to program (39a)±(39c) we have then to solve the convex optimization problem
43a
s.t. A 0 X b0 X $ 0;
s X j1
aj Q
q0 ; hvj 2 T j X;
44a
where c U Ec
v:
44b
Introducing auxiliary variables u
u1 ; ¼; uj ; ¼; us 0 ; program (43a)±(43c) is equivalent to min c 0 X 1
s X j1
aj uj
45a
s.t. A0 X b0
45b
uj $ Q
q0 ; hvj 2 T j X; j 1; ¼; s
45c
X $ 0;
uj $ 0;
j 1; ¼; s:
45d
Note that the non-negativity constraint u $ 0 is due to the fact that we may assume that Q $ 0; see LP (35a)±(35c) and (36). The key to the dual decomposition method for solving program (43a)±(43c) lies in the fact that the values Q
q0 ; hvj 2 T j X; 1 # j # s; of the recourse function Q can be represented, using LP-duality theory, by Q
q0 ; hvj 2 T j X max w 0
hvj 2 T j X;
46a
cf. Eqs. (31a)±(31c), (35a)±(35f), (38a)±(38c), (6), (11a) and (11b), where the feasible domain U of the dual LP reads ! ( ) u 0 0 1 2 ~ 2q0 # u~ # q0 : UU w : C u H u;
46b u~ Thus, representing Eq. (45c) by
In this section we assume that the linearity conditions (33a) and (33b) hold and a
v has a discrete distribution as described by Eqs. (41a) and (41b). Furthermore, for simpli®cation suppose ®xed second-stage cost factors
min E
c
v 0 X 1 Q
q0 ; hv
v 2 T
vX
G
X U c 0 X 1
w[U
See Refs. [9,10].
q
v q0 a:s:
given by
43b
43c
where c
v; T
v and hv
v are de®ned according to Eqs. (33a), (33b) and (41b). Based on the above assumptions, the objective function G G
X of problem (43a)±(43c) is
uj $ w 0
hvj 2 T j X;
w [ U; j 1; ¼; s;
45e
Program (45a)±(45d) can be written also as a semi-in®nite LP. Corresponding to the approximation p~ i of the distance or Minkowski functional pi of Ki, see Eqs. (11a), (11b) and (12a), the recourse function Q is approximated iteratively by ~ 0 ; hvj 2 T j X U max w 0
hvj 2 T j X; Q
q0 ; hvj 2 T j X < Q
q w[Ujn
47 where Ujn , U; n 0; 1; ¼; are ®nite sets which are de®ned recursively. In detail, the dual decomposition method proceeds as follows: Step 1: initialization Since Q $ 0; solve the approximative LP, cf. program (45a)±(45d), min c~ 0 X 1
s X j1
aj uj
48a
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
s.t. A 0 X b0
u $ 0:
X $ 0;
48b
the dual decomposition method yields an optimal design point X after ®nitely many steps.
48c
Proof.
^ u^
u^ j denote an optimal solution of this LP. Let
X; Step 2: optimality cuts ^ j 1; ¼; s; solve the second-stage For z j U hvj 2 T j X; program (35a)±(35c) producing a primal, dual optimal solution y^ j ; w^ j ; respectively, j 1; ¼; s: Hence, 0
^ w^ j
hvj 2 T j X; ^ j 1; ¼; s; Q
q0 ; hvj 2 T j X
49
cf. Eqs. (46a) and (46b). Because of Eq. (45e) we consider now the constraints 0
uj $ w^ j
hvj 2 T j X; j 1; ¼; s:
50 ^ ^ Case 2A:
X; u ful®ls Eq. (50) According to Eqs. (49) and (46a), X^ is optimal in program (43a)±(43c), and the algorithms stops. Note that this argument remains true also in later steps, since only constraints of the type Eq. (50) are added to the preceding approximate LP, see Eq. (45e). ^ u^ violates Eq. (50) Case 2B:
X; In this case there is at least one index j0 ; 1 # j0 # s; such that ^ q 00 y^ j0 w^ j0 0
hvj0 2 T j0 X ^ Q
q0 ; hvj0 2 T j0 X u^ R 2 u^~ j0 0
j0
j0 0
^ F0
s yj0 ; X
^ . u^ j $ 0; U 2Ml0
a ; X 0 j0
51 where Ml0 Ml0
a j0 ; X; 1 # l0 # l0 ; denotes the limit state function corresponding to the extreme point w^ j0 of U, see Eq. (28b). This means that X^ lies in a certain failure domain with positive probability, and the related minimum failure costs exceed the actual optimal cost level u^ j0 : Hence, add to the preceding LP the so-called optimality cuts (50), ^ u^ : cutting off the nonoptimal solution
X; Step 3: Update of the approximate LP If D^ denotes the feasible domain of the current approximate LP, rede®ne D~ U {
X; u [ D^ :
X; u fulfils Eq:
50};
52 and solve the updated LP s X aj u j min c 0 X 1 s.t.
53a
j1
~
X; u [ D;
53b ^ u^ U
X; ~ u~ ; ~ u~ : With
X; yielding an optimal solution
X; return to Step 2. Under weak assumptions, this algorithms ®nds an optimal solution X p of program (43a)±(43c) after ®nitely many steps: Theorem 4.1.
173
If {X : A0 X b0 ; X $ 0} is bounded, then
See Refs. [10].
4.2. Sequential discretization For simpli®cation, suppose again that Eq. (42) holds, i.e. the second-stage cost coef®cients are ®xed. For an arbitrary random matrix
T
v; hv
v in program (43a)±(43c), by discretization we can ®nd a sequence
T1
v; h1v
v;
T2
v; h2v
v; ¼;
TN
v; hNv
v; ¼
54a of piecewise constant random matrices
TN
´; hNv
´; N 1; 2; ¼; converging to
T
´; hv
´ in some probabilistic sense. Consequently, we obtain then the approximative objective function GN
X U c 0 X 1 EQ
q0 ; hNv
v 2 TN
vX;
54b
cf Eq. (44a). Under weak assumptions [18] we have that GN
x ! G
X; N ! 1;
54c
uniformly on compact subsets of Rr : Moreover, the approximative program min GN
X
55a
s.t. A 0 X b0
55b
X$0
55c
has the following properties. Theorem 4.2. 1. If
TN
´; hNv
´ is a piecewise constant random matrix having ®nitely many realizations, then program (55a)± (55c) can be represented by a linear program with a dual decomposition data structure as described by Fig. 1. 2. Suppose that D0 {X [ Rr : A0pX b0 ; X $ 0} is a compact set, and let denote XN an optimal solution of the approximate problem (55a)±(55c). Then the optip mal value GN of (55a)±(55c) converges to the optimal value G p of program (43a)±(43c), as N ! 1; and every p accumulation point X p of
XN is an optimal solution of the original problems (43a)±(43c). Proof.
See Refs. [9,18].
Remark 4.1. By special iterative re®ning techniques applied in the sequential discretization
TN
´; hNv
´ !
TN11
´; hN11
´; N 1; 2; ¼; v of
T
´; hv
´; the increase of the size of the resulting (large scale) LPs (problem (55a)±(55c)) can be kept moderately.
174
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
4.3. Convex approximation by partial linearization If X ! G0
g
v; X and/or X ! T
v; X are nonlinear mappings, then by linearization of the function(s) G0
g
v; X < G0
g
v; X0 1 7 x G0
g
v; X0 0
X 2 X0 ;
56a T
v; X < T
v; X0 1
2T
v; X0
X 2 X0 2X
56b
with respect to X at a certain reference point X0, the original problem (43a)±(43c) can be reduced to a convex problem ful®ling Eqs. (33a) and (33b); let X1 denote its optimal solution and put X0 U X1 : Proceeding iteratively, a minimizing sequence X0 ; X1 ; ¼ for the original problem (43a)± (43c) is obtained. For more details, see Ref. [19].
at time t1 of the structure built using the initial design X at time t0, then Eq. (58) may be given by the linear model F 0
s y ; X; j U F 0I
s yI ; X 1 FII0
s yII ; j;
59b
where the vectors F 0I ; FII0 are de®ned again as in Eq. (15b).
Corresponding to program (21a)±(21e) in Section 1, for the optimal selection of the initial design X X
Pa
´ and the second-stage design function j j
a we have then this random program min G 0
g
v; X 1 G1
z
v; j
a
v
60a
5. Adaptive designs
s.t.
In this section we assume that we have the following adaptive decision situation, see e.g. Ref. [29]. Stage I at time t0: An initial or pre-design X must be selected at a time t0 based on limited information on the model parameters and structural loading represented by the probability distribution Pa
´ of the total parameter vector a(´). Stage II at a later time t1 . t0 : The random vector a(´) is realized a a
v at a later time instant t1 . t0 ; which enables then a correction, completion j j
a of the initial design, a strengthening of the present structure by an additional design action j , etc. Besides the costs G 0 G0
g
v; X caused by the inital or pre-design X, we have now additional costs
CF R
v
60b
HF # F 0
s
v; X; j
a
v
60c
A 0 X b0
60d
A1 j
a
v b 1
60e
G1 G1
z
v; j
a; a U
s y ; s L ; s U ; Ei ; g; h; z; R; ¼;
57 caused by the second-stage design action j . Working again with the yield conditions (15a) and (15b), see Section 1, the vector F0 of principal axial plastic capacities in Eq. (15b) depends now F 0 F0
s y ; X; j
58
on the initial design vector X X
P a
´ and on the correction/completion or strengthening design vector j j
a which can be taken at time t 1 . t0 under full information about the realization a a
v of a(´). For the de®nition (58) of F0 in this two-stage decision situation we mention the following models: 1. If j j
a is a vector of design corrections for the initial design X, then F 0
s y ; X; j F 0
s y ; X 1 j;
59a
where the right hand side of Eq. (59a) is de®ned as in Eq. (15b). 2. If j j
a represents certain pre-planned completions
X $ 0;
j L # j # j U;
60f
where j L ; j U ; are given bounds for j which may depend on X, and (A1,b1) is a given matrix. The transformation of program (60a)±(60f) into a twostage stochastic (linear) program can be done now by the stochastic optimization tools described in Sections 3 and 4. To give an example, in case (59a), and assuming then that G 0
g
v; X U c
v 0 X;
61a
G1
z
v; j : d
v 0 j;
61b
F 0
s y
v; X 1 j U S
v
X 1 j S
vX 1 S
vj
61c with a certain matrix S
v S
s y
v; for the minimization of the expected total costs we get a stochastic linear program with recourse of the type (32a)±(32d) having the following main components, cf. Eqs. (31a)±(31c): T
v; X T
vX U
0 2S
vX
! ;
h
v U
R
v 0
! ;
62a
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
Fig. 2. Three-bar truss, loading.
WU
0
0
I
2I
0
D1
v
0
0
2S
v H 1
C B B D2
v C C B C; y
v U B C B B j
v C A @ F
v 1 0 1 q
v C B B q2
v C C B C: q
v U B C B B d
v C A @ 0
!
recourse matrix
0
D1
v
62b
1
B 2 C C y I
v U B @ D
v A;
j
v
62c
6. Numerical examples 6.1. Three-bar truss with random load Consider ®rst the plane three-bar truss according to Fig. 2. In this simple case the equilibrium matrix C is given by ! k 0 2k 1 C with k p :
63 2 2k 21 2k
175
Assume that the allowable stresses are ®xed and given (Fig. 3) by s L 2190 N=mm2 and s U 216 N=mm2 for all members. Moreover, suppose that the external load is a random vector R R
v
a1
v; a2
v 0 such that a1
v has a uniform distribution over a given support 2r; r with radius r $ 0 and a2
v 2105 N for all v . By using discretization techniques [23], the distribution of the load vector R R
v is approximated by a discrete distribution with mean ER
v
0; 2105 N 0 and 100 possible realizations which are taken with probabilities pj < 1=100; j 1; ¼; 100: Let the 1st stage P cost function G0 be de®ned by the volume G0
v; X 3k1 Xk Lk with the design vector A
X U
X1 ; X2 ; X3 0 of cross-sectional areas X1, X2, X3 of the bars. Based on Eqs. (40a) and (40b), the cost factors qL qU q2 are chosen by Eq. (40b) with the ®xed elastic modulus s E 7; 2 £ 10 4 N=mm2 for all members. If D0 is still de®ned by D0 {X [ R3 ; AX # b; X $ 0} with A
1; 1; 1; b 105 mm2 ; then the substitute problem with cost function (30a), see Section 3, has the following LP-representation, cf. Fig. 3: Fig. 4 below shows the dependence of the optimal crosssectional areas on the radius r of the supporting interval of the discrete distribution for load R1
v; where r is proportional to Var
R1
´: In comparison to the solution of the deterministic problem (r 0; where the unknown parameters are replaced by their mean values), the cross-sectional area of bar 1 and 3 are increasing, while the area of bar 2 is decreasing for increasing uncertainty of the load; furthermore, the expected total costs are increasing. Considering discrete distributions with a large number of realizations, the linear programs may become very large. The following Tables 1 and 2 show the dependence of the elapsed time (time) and the number of iterations (it) on the number (nr) of realizations of the random load vector for two classes of LP-solvers, cf. Ref. [23]: 1. General purpose (standard) LP-solver: OSL, CONOPT. 2. Special purpose LP-solver (implemented in SLP-IOR)
Fig. 3. LP-data structure of the example in Section 6.1.
176
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
Fig. 4. Optimal cross-sectional areas.
exploiting the dual decomposition data structure of the problem: SHOR_2, BPMPD, SDECOM. The results below were obtained on a Pentium II workstation with 333 MHz, where the terminal values of the total objective function of each test problem obtained by the various LP-solvers differ at most 0.02%. Comparing Tables 1 and 2 and referring to the above remarks, we ®nd that the series of the test problems may be solved by both classes of LP-solvers with the same accuracy. However, by using special purpose LP-solvers, the computational expenses may be reduced drastically! Furthermore, we observe, that in general the total number of iterations and the total elapsed time is increasing for a growing problem size.
example in Section 6.1. Moreover, the 1st stage cost function is again the volume of the truss, and the design vector X is the vector of cross-sectional areas (Fig. 6). The external load acts on two knots only, one is stochastic with mean ER1
v
105 N; 0 N 0 ; the other, R2
0 N; 2105 N 0 is ®xed. The random load has now a discrete distribution with ®ve possible realizations taken with equal probabilities pj 1=5; j 1; ¼; 5: Fig. 8 below shows the dependence of the optimal cross-sectional areas on the radius r of the supporting interval of the discrete distribution for load R1
v; where r
6.2. Ten-bar truss with random load The second example is a plane ten-bar truss, cf. Fig. 5. The material parameter
s L ; s U ; E are the same as in the Table 1 Standard LP-solvers
Fig. 5. Initial structure ten bar-truss, loading.
OSL
CONOPT
nr
it
Time
it
Time
100 400 900 1600 2500
551 2338 5172 9050 14294
0.69 12.57 79.976 274.77 689.89
146 445 977 1988 2701
0.78 6.37 34.26 128.33 296.74 Fig. 6. Optimal solution of the deterministic problem (- - - deleted bar).
Table 2 Special purpose LP-solvers SHOR_2
BPMPD
SDECOM
nr
it
Time
it
Time
it
Time
100 400 900 1600 2500
98 102 104 108 106
8.84 27.03 55.71 94.56 138.9
12 12 13 13 13
0.54 1.71 4.23 8.01 13.12
2155 2424 2449 2576 2572
42.79 56.96 59.26 63.28 65.74
Fig. 7. Optimal solution of SOP (- - - deleted bar).
K. Marti / Reliability Engineering and System Safety 72 (2001) 165±177
Fig. 8. Optimal cross-sectional areas.
is proportional to Var
R1
´: Starting with the initial ten-bar structure, cf. Fig. 5, the optimal solution of the deterministic problem with R1 U ER1
v yields an optimal design such that the bars 5, 6, 7, 10 may be deleted. In comparison to this solution, in the optimal solution of the substitute problem described above, bar 5 may be deleted only, where the crosssectional areas of bar 1, 2, 3, 6, 7 and 10 are increasing, while the areas of bar 4, 8 and 9 are decreasing for increasing uncertainty 0 r1 , r2 , r3 , r4 , r5 of the load (see Figs. 6 and 7). Furthermore, the expected total costs are increasing. Acknowledgements The author would like to thank Ms Elisabeth LoÈûl for setting the LaTex ®le of the manuscript and Dipl.Math. Gerald StoÈckl, UniBw MuÈnchen, Fak. LRT, for providing the numerical results given in Section 6. References [1] Augusti G, Baratta A, Casciati F. Probabilistic methods in structural engineering. London/New York: Chapman and Hall, 1984. [2] Birge JR. Stochastic programming computation and applications. INFORMS J Computing 1997;9(2):111±33. [3] Charnes A, Lemke CE, Zienkiewicz OC. Virtual work, linear programming and plastic limit analysis. Proc R Soc A 1959;25:110±6. [4] Cocchetti G, Maier G. Static shakedown theorems in piecewise linearized poroplasticity. Arch Appl Mech 1988;68:651±61. [5] Haftka RT, Kamat MP. Elements of structural optimization. The Hague/Boston/Lancaster: Martinus Nijhoff/Kluwer, 1985. [6] Hemp WS. Optimum structures. Oxford: Oxford University Press, 1973. [7] Heyman J. Plastic design of beams and plane frames for minimum material consumption. Q Appl Math 1950/51;8:373±81. [8] Kaliszky S. PlastizitaÈtslehre. DuÈsseldorf: VDI, 1984. [9] Kall P. Stochastic linear programming. Berlin/Heidelberg/New York: Springer, 1976. [10] Kall P, Wallace SW. Stochastic programming. Chichester: Wiley, 1994. [11] Kirsch U. Structural optimization. Berlin/Heidelberg/New York: Springer, 1993.
177
[12] KoÈnig JA. Shakedown of elastic±plastic structures. Amsterdam: Elsevier, 1987. [13] Krabbenhoft K, Damkilde L. Limit analysis based on lower-bound solutions and nonlinear yield criteria. In: Topping BHV, editor. Finite elements: techniques and developments. Edinburgh: Civil-Comp Press, 2000. p. 117±29. [14] Krenk S, Damkilde L, Hoyer O. Limit analysis and optimal design of plates with equilibrium elements. J Engng Mec 1994;120(6):1237± 54. [15] Lawo M. Optimierung im kunstruktiven ingenieurbau. Braunschweig/ Wiesbaden: Vieweg, 1987. [16] Maier G, Munro J. Mathematical programming applications to engineering plastic analysis. Appl Mec Rev 1982;35(12):1631±43. [17] Marti K. Entscheidungsprobleme mit linearem aktionen- und ergebnisraum. Z Wahrscheinlichkeitstheorie verw Geb 1972;23:133±47. [18] Marti K. Approximationen der entscheidungsprobleme mit linearer ergebnisfunktion und positiv homogener, subadditiver verlustfunktion. Z Wahrscheinlichkeitstheorie verw Geb 1975;31:203±33. [19] Marti K. Approximationen stochastischer optimierungsprobleme. KoÈnigstein/Ts: Verlag Anton Hain Meisenheim GmbH, 1979. [20] Marti K. Approximation and derivatives of probabilities of survival in structural analysis and design. Struct Optimization 1997;13:230±43. [21] Marti K. Solving stochastic structural optimization problems by RSM-based stochastic approximation methods-gradient estimation in case of intermediate variables. MMOR 1997;46(3):409±34. [22] Marti K. Optimal design of trusses as a stochastic linear programming problem. In: Nowak AS, editor. Reliability and optimization of structural systems, Ann Arbor, MI: University of Michigan Press, 1998. p. 231±9. [23] Mayer J. Stochastic linear programming algorithms. Amsterdam: Gordon and Breach, 1998. [24] MSC/NASTRAN Module CMG, MacNeal-Schwendler Corp. Munich, November 1998. [25] Nielsen MP. Limit analysis and concrete plasticity. Englewood Cliffs, NJ: Prentice Hall, 1984. [26] Prager W, Shield RT. A general theory of optimal plastic design. J Appl Mec 1967;34(1):184±6. [27] Rozvany GIN. Structural design via optimality criteria. Dordrecht/ Boston/London: Kluwer Academic, 1989. [28] Rozvany GIN. Topology optimization in structural mechanics. Wien/ New York: Springer, 1997. [29] SchuÈler H. Zur Analyse und Bemessung adaptiver Tragwerke aus Stahlbeton unter dynamischen Einwirkungen. Dissertation, Bauhaus-UniversitaÈt Wien, 1997. [30] SchueÈller GI, Gasser M. Some basic principles of reliability-based optimization (RBO) of structures and mechanical components. In: Marti K, Kall P, editors. Stochastic programming Ð numerical techniques and engineering applications, LNEMS Vor. 458. Berlin/ Heidelberg: Springer, 1998. p. 80±103. [31] Smith DL, editor. Mathematical programming methods in structural plasticity CISM course and lectures, vol. 299. Wien/New York: Springer, 1990. [32] Smith DL, Munro J. Plastic analysis and synthesis of frames subject to multiple loadings. Engng Optimization 1976;2:145±57. [33] Spillers WR. Automated structural analysis: an introduction. New York/Toronto/Oxford: Pergamon, 1972. [34] Steffen PN. Elastoplastische dimenionierung von stahlbetonplatten mittels ®niter bemessungselemente und linearer programmierung. Dissertation, ETH ZuÈrich, Nr. 11611, ZuÈrich, 1996. [35] Tin-Loi F. On the optimal plastic synthesis of frames. Engng Optimization 1990;16:91±108. [36] Tin-Loi F. Plastic limit analysis of plane frames and grids using GAMS. Comput & Struct 1995;54:15±25. [37] Zeman P, Irvine HM. Plastic design. An imposed hinge-rotation approach. London/Boston/Sydney: Allen and Unwin, 1986.