Shape and topology optimization of plane frames via quantile-based sequential optimization and reliability assessment with global stability constraint

Abstract

In this paper a framework of quantile-based sequential optimization and reliability assessment (SORA) is extended to consider the global stability constraint in the optimization of plane frames. Uncertainty is considered on the structural side without any prior assumptions on their distribution information, and two novel stopping criteria with a reduction coefficient are employed to smoothly shift the constraint boundary for the next iteration. Force density method is introduced for the shape optimization of plane frames to avoid the existence of the melting nodes, and the geometrical stiffness matrix is also penalized to exclude pseudo local buckling modes. The numerical examples illustrate that with the help of reduction coefficients, the shifts of different constraint boundaries in SORA become smoother and the convergence of sequential optimization is improved, and due to shape optimization, the reliability of structural stability can be satisfied with limited increase of structural volume of the one without considering stability. Moreover, it is also shown in the cantilever beam and bridge examples that global structural stability can be enhanced by applying nodal displacement constraints with higher reliability.

Graphical abstract

Keywords

quantile function / sample L-moment / sequential optimization and reliability assessment / force density method / stability constraint / maximum entropy method

1 Introduction

With the emergence of high performance computing and the development of optimization algorithms, structural optimization has become a useful tool for practical design and engineering. Among all types of structures, the optimization of skeletal structures, including trusses and frames, is a popular topic due to its wide applicability to structures in various fields of engineering, and enormous literatures have contributed to this research field [1–4]. When the nodal positions remain fixed, achieving the optimal structural layout involves eliminating redundant members from the array of potential node connections, which is called the ground structure method (GSM) [5] for topology optimization. Although GSM has been extensively studied and explored due to its simplicity, it is recognized that the optimal shape of the structure can further improve the structural performance by changing the nodal locations and member connectivity [6,7]. However, when the locations of nodes are considered as design variables and are allowed to widely move, one of the main challenges is the existence of melting nodes, which leads to the existence of very short members in the structure, causing the stiffness matrix to become singular [8–10]. To alleviate this difficulty, force density method (FDM) is introduced to determine the nodal locations of pin-jointed trusses for shape and topology optimization [11,12]. Shen and Ohsaki [13] extended this method to rigidly-jointed frames by using an auxiliary truss structure to which the FDM is applied. However, due to the random nature in the material and geometrical properties of real-world structures, uncertainty should be considered in the optimization problem.

As evidenced by a large number of literatures [14–16], structural optimization considering uncertainty has increasingly played an important role in practical engineering. Numerous methods have been developed to address uncertainty in the optimization process, and among them one of the popular approaches is the so-called reliability-based design optimization (RBDO), which aims to find the optimal solution while satisfying the target reliability constraints [17,18]. While valuable efforts have been made to enhance the performances of the methods for solving RBDO problem, most of them calculate the structural reliability by searching the most probable point (MPP) [19,20], and it may face with a convergence issue when multiple MPPs exist [21,22]. To alleviate this difficulty, Moustapha et al. [23] proposed a double-loop quantile-based RBDO problem formulation, where the structural reliability is evaluated by using Kriging surrogate model and is nested in the structural optimization procedure [24]. To improve the computational efficiency, Li et al. [25] proposed a quantile-based sequential optimization and reliability assessment (SORA), in which the original RBDO problem is reformulated as a sequence of deterministic optimization problems [20]. Ghasemi et al. [26,27] proposed a sequential optimization method to solve RBDO, where in each iteration a double-stage optimization procedure is implemented to obtain an optimal solution with satisfying reliability constraints. However, due to the lack of knowledge about the exact uncertainty distribution in real-world scenarios, one often has to rely on random samples to assess structural reliability [28]. To reduce the influence of sample size and variability on probability estimation, Pandey [29] approximated the quantile function of a nonnegative random variable using the maximum entropy method (MEM), subject to constraints on sample probability weighted moments. Based on this idea, Shen et al. [30] recently proposed a quantile-based SORA for shape and topology optimization of frame structures using FDM, where the structural reliability associated with nodal displacement is estimated using the MEM, subject to constraints on sample L-moments, while structural stability is not considered in the formulation of optimization problem.

The optimization process for trusses and frames is commonly observed to converge toward structures that exhibit inadequate stability, and the need to include stability constraint into optimization problem has been discussed in many studies [31,32]. Although formulations with different stability constraints, such as nodal stability [33], global stability (also known as linear buckling) [34] and nonlinear buckling [35], have been proposed to ensure structural stability, the global stability constraint is widely used to obtain a stable structure because of its easy implementation [36]. To capture the local and global buckling during optimization, Torii et al. [37] simulated the truss-like structures using frame elements, and formulated the problem with only one global stability constraint. Ohsaki et al. [38] utilized semidefinite programming to solve the truss topology optimization problem with a frequency constraint, and Kanno and Ohsaki [39] further extended this method to frame optimization involving a prescribed linear buckling load factor, with the constraint formulated using its reciprocal. With this reciprocal expression of linear buckling load factor, Shen et al. [40] proposed a robust optimization of plane frames with global stability constraint, where the uncertainties in nodal locations are introduced such that the frame member is no longer straight, and a penalization scheme is used to exclude the pseudo local buckling mode caused by the geometrical stiffness of thin elements. However, since uncertainty is incorporated during the optimization procedure, it is challenging for large-scale problems due to its intensive computational efforts.

Therefore, based on the optimization framework of Ref. [30], in this study we extended our quantile-based SORA method to incorporate stability constraint. Uncertainty is only considered on the structural side, which is assumed to exist in nodal locations, cross-sectional areas and Young’s modulus with unknown distribution, and to clearly address the effect of displacement constraint on structure stability, the optimization problem is formulated as minimizing the structural volume under displacement and global stability constraints. The reliability constraint is formulated based on quantile of structural performance, which is estimated by MEM subject to sample L-moments. Since both displacement constraints and stability constraints are included, the solution spaces among subsequent iterations of quantile-based SORA could be very different and induce drastic shift of the reliability constraint, which will regrade the convergence property. Therefore, a reduction coefficient is employed to reduce the shift of constraint when it is larger than the value of constraint boundary of previous iteration. In addition, the FDM in Ref. [13]. is introduced for shape optimization of plane frames, and a penalization scheme is applied to element geometrical stiffness matrix to accurately calculate the linear buckling load factor and its sensitivity coefficients.

The structure of the paper is as follows. In Section 2 the FDM is briefly introduced, and Section 3 gives the formulation of quantile-based SORA with global stability constraint, together with the estimation of quantile function and the definition of uncertainty in member stiffness. Section 4 summarizes the penalization method for element geometrical stiffness matrix and the corresponding sensitivity analysis. Section 5 examines three numerical examples to illustrate the efficacy of the proposed method, as well as the new findings. Finally, Section 6 summarizes the conclusions of the study.

2 Shape optimization using force density method

In this section, the FDM is applied to an auxiliary pin-jointed truss, in which only axial force exists, to determine the nodal locations of frame for shape optimization. Therefore, it is important to note that the structure considered in this section is a pin-jointed truss, and the nodal layout of the frame structure to be optimized is defined by applying FDM to the auxiliary pin-jointed structure, where the auxiliary truss structure has independent loading and supporting conditions and is irrelevant to the frame to be optimized. Further information can be found in Ref. [13].

Let m denote the number of members and n denote the number of nodes in the auxiliary pin-jointed truss, which are the same as those of the frame to be optimized. Assume that the ith member connects two nodes i₁ and i₂. Then the entry of m-by-n connectivity matrix C is defined as Refs. [13,41]

(1){\mathit{C}}_{\left(i,j\right)}=\{\begin{array}{ll}1,& j=i_1,\\ -1,& j=i_2,\\ 0,& \mathrm{e}\mathrm{l}\mathrm{s}\mathrm{e},\end{array}(i=1,\mathrm{2...},m;i_1,i_2=1,2,...,n),

where subscript of (i,j) refers to the specific entry of matrix C located at the ith row and jth column. Let t_i be the force density of the ith member which is defined as

(2)t_i=\frac{N_i}{L_i},

where N_i stands for the axial force and L_i refers to the length of the ith member. Let x_free, y_free, x_fix, and y_fix be the x- and y-coordinates of free nodes and fixed nodes, respectively, and \mathit{t}={\left(t_1,t_2,\dots ,t_m\right)}^{\mathrm{T}} is the vector of force density. Note that the nodes with the fixed locations in the optimization process of the frame are included in ‘fixed nodes’ of the auxiliary truss, although their displacements may not be fixed in the frame to be optimized. Rearrange the columns of connectivity matrix as \mathit{C}=\left({\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},{\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}}\right) in which the columns of {\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}} and {\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}} correspond to the free nodes and fixed nodes, respectively, and denote the external loads applied to free nodes and fix nodes as P_x,free, P_y,free, P_x,fix, and P_y,fix, respectively. Then the equilibrium equations corresponding to the free and fixed nodes of the auxiliary pin-jointed truss are then formulated as follows Refs. [13,41].

(3)\begin{array}{c}{\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}{\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}+{\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}}{\mathit{x}}_{\mathrm{f}\mathrm{i}\mathrm{x}}={\mathit{P}}_{x,\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\\ {\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}+{\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}}{\mathit{y}}_{\mathrm{f}\mathrm{i}\mathrm{x}}={\mathit{P}}_{y,\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\\ {\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}{\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}+{\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}}{\mathit{x}}_{\mathrm{f}\mathrm{i}\mathrm{x}}={\mathit{P}}_{x,\mathrm{f}\mathrm{i}\mathrm{x}},\\ {\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}+{\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}}{\mathit{y}}_{\mathrm{f}\mathrm{i}\mathrm{x}}={\mathit{P}}_{y,\mathrm{f}\mathrm{i}\mathrm{x}},\end{array}

where \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right) is the diagonal matrix by placing the components of vector t along its main diagonal. Details of Eqs. (1)–(3) and related illustrative example can be found in Refs. [11,41]. The fixed nodes are designated as the supporting and loading nodes, and therefore the vectors P_x,free and P_y,free are both zero vectors. Based on Eq. (3), the locations of free nodes can be obtained by a set of force densities t as

(4)\begin{array}{c}{\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}=-{\left({\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\right)}^{-1}{\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}}{\mathit{x}}_{\mathrm{f}\mathrm{i}\mathrm{x}},\\ {\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}=-{\left({\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\right)}^{-1}{\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{t}\right){\mathit{C}}_{\mathrm{f}\mathrm{i}\mathrm{x}}{\mathit{y}}_{\mathrm{f}\mathrm{i}\mathrm{x}}.\end{array}

It has been proven in Ref. [42]. that if at least one node is fixed in the structure and t is nonzero, the matrix {\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}^{\mathrm{T}}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathit{q}\right){\mathit{C}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}} in Eq. (4) is invertible and there is a solution for Eq. (4). Moreover, it has been investigated in Ref. [13]. that in Eq. (3) when N_i is nonzero, extremely short members can be avoided by indirectly controlling L_i with bounded t_i.

3 Quantile-based sequential optimization and reliability assessment for frame optimization

Based on the framework in Ref. [30], in this section the RBDO is proposed for plane frames with displacement and stability constraints, and the reduction coefficient is introduced to prevent drastic shift of the constraint boundary. Then calculation of quantile estimation is given, where the MEM is used under constraints on sample L-moments. Uncertainty modeling of geometrical and material parameters is presented in the last part of the section.

3.1 Problem formulation

As discussed in Section 2, because the fixed nodes cannot move throughout the optimization procedure, it is considered that the design variables consist of locations of free nodes x_free and y_free and cross-sectional areas of members A. Let E represents the Young’s modulus of members, and uncertainty is assumed to be the random perturbation on locations of nodes \mathrm{\Delta }\mathit{x} and \mathrm{\Delta }\mathit{y}, cross-sectional areas \mathrm{\Delta }\mathit{A} and Young’s modulus \mathrm{\Delta }\mathit{E}, which are written as \theta =\left(\mathrm{\Delta }\mathit{x},\mathrm{\Delta }\mathit{y},\mathrm{\Delta }\mathit{A},\mathrm{\Delta }\mathit{E}\right), and their details will be given in Subsection 3.3.

Consider an RBDO problem where the objective is to minimize the total volume of the structure subject to displacement and global stability constraints as below:

(5)\begin{array}{rl}& \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}V\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\mathit{A}\right),\\ & \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}Pr\left\{g_j\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\mathit{A};\theta \right)⩽{\overline{g}}_j\right\}⩾R_j,j=1,2,\dots ,n_g,\\ & Pr\left\{\frac{1}{{\lambda }^{\mathrm{c}\mathrm{r}}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\mathit{A};\theta \right)}⩽\frac{1}{{\lambda }_L}\right\}⩾R_{\lambda },\\ & {\underset{\_}{\mathit{x}}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}⩽{\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}⩽{\overline{\mathit{x}}}_{{}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}},{\underset{\_}{\mathit{y}}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}⩽{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}⩽{\overline{\mathit{y}}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\underset{\_}{\mathit{A}}⩽\mathit{A}⩽\overline{\mathit{A}},\end{array}

where V\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\mathit{A}\right) is the total structural volume, and g_j\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\mathit{A};\theta \right) represents the jth displacement component with upper bound {\overline{g}}_j and target probability R_j; {\lambda }^{\mathrm{c}\mathrm{r}}\left(\mathit{x},\mathit{y},\mathit{A};\theta \right) is the linear buckling load factor with lower bound {\lambda }_{\mathrm{L}} and target probability R_{\lambda }; and the lower and upper bars on x_free, y_free and A indicate the corresponding lower and upper bounds, respectively. Note that {\lambda }^{\mathrm{c}\mathrm{r}} is the smallest positive result of the eigenvalue problem below:

(6)\left(\mathit{K}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\mathit{A};\theta \right)-\lambda {\mathit{K}}_{\mathrm{G}}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\mathit{A};\theta \right)\right)\mathit{\Phi }=0,

where K and K_G are global elastic and geometrical stiffness matrices with shape of n_dof -by-n_dof, respectively, and n_dof represents the total number of degrees of freedom (DOFs) of the structure; \mathit{\Phi } is the eigenvector corresponding to the eigenvalue\lambda. In Eq. (5) the global stability constraint is formulated by taking the reciprocals of {\lambda }^{\mathrm{c}\mathrm{r}} and {\lambda }_{\mathrm{L}}, ensuring that {\lambda }^{\mathrm{c}\mathrm{r}} is either greater than {\lambda }_{\mathrm{L}}or negative, as adopted in some other researches [37,39]. Henceforth, for simplicity we denote the reciprocals as {\gamma }^{\mathrm{c}\mathrm{r}}=1/{\lambda }^{\mathrm{c}\mathrm{r}} and \overline{\gamma }=1/{\lambda }_{\mathrm{L}}.

By incorporating Eq. (4) in Section 2, the RBDO problem (5) is written as

(7)\begin{array}{rl}& \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}V\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right),\\ & \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}Pr\left\{g_j\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A};\theta \right)⩽{\overline{g}}_j\right\}⩾R_j,\\ & j=1,2,\dots ,n_g,\\ & Pr\left\{{\gamma }^{\mathrm{c}\mathrm{r}}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A};\theta \right)⩽\overline{\gamma }\right\}⩾R_{\gamma },\\ & \underset{\_}{\mathit{t}}⩽\mathit{t}⩽\overline{\mathit{t}},\underset{\_}{\mathit{A}}⩽\mathit{A}⩽\overline{\mathit{A}}.\end{array}

Furthermore, let Q_R denote the quantile with specified probability R. Then, according to the equivalence between the probability and the quantile structural response, we can define the following n_g + 1 quantiles for the displacement and global stability constraints [30,43]

(8)\begin{array}{rl}& Q_{R_j}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A};\theta \right)\\ =& inf\left\{Q:Pr\left\{g_j\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A};\theta \right)⩽Q\right\}⩾R_j\right\},\\ & j=1,2,\dots ,n_g,\\ & Q_{R_{\gamma }}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A};\theta \right)\\ =& inf\left\{Q:Pr\left\{{\gamma }^{\mathrm{c}\mathrm{r}}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A};\theta \right)⩽Q\right\}⩾R_{\gamma }\right\}\end{array},

as substitutes to the probability constraints in problem (7), which can be restated as

(9)\begin{array}{c}\mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}V\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right),\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{Q}_{R_j}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A};\theta \right)⩽{\overline{g}}_j,j=1,2,\dots ,n_g,\\ Q_{R_{\gamma }}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A};\theta \right)⩽\overline{\gamma },\underset{\_}{\mathit{t}}⩽\mathit{t}⩽\overline{\mathit{t}};\underset{\_}{\mathit{A}}⩽\mathit{A}⩽\overline{\mathit{A}}.\end{array}

The optimal solution to problem (9) can be derived through the resolution of a sequence of deterministic optimization problems [30], which is achieved by shifting the upper bounds of constraint functions to align with the feasible domain. Let the shift of constraint boundaries {\overline{g}}_j and \overline{\gamma } at (k + 1)th iteration be {\overline{c}}_j^{k+1} and {\overline{c}}_{\gamma }^{k+1}, respectively, which are calculated as

(10)\begin{array}{rl}{\overline{c}}_j^{k+1}=& Q_{R_j}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{A}}^k;\theta \right)\\ & -g_j\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{A}}^k\right),j=1,2,\dots ,n_g,\\ {\overline{c}}_{\gamma }^{k+1}=& Q_{R_{\gamma }}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{A}}^k;\theta \right)\\ & -{\gamma }^{\mathrm{c}\mathrm{r}}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{A}}^k\right),\end{array}

where t^k and A^k are the solutions of the deterministic optimization problem at the kth iteration as described below in Eq. (11); g_j({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}({\mathit{t}}^k),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}({\mathit{t}}^k),{\mathit{A}}^k) and {\gamma }^{\mathrm{c}\mathrm{r}}({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}({\mathit{t}}^k),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}({\mathit{t}}^k),A^k) are the values of constraint functions without uncertainty. It is worth noting that in our previous research [30] when uncertainty is taken into account, large shifting value of displacement constraint could happen to satisfy the reliability requirement if the structure is under unstable state, and a narrow feasible domain could be also constructed (see numerical example 3 in Ref. [30]). Since in the present study both displacement and global stability constraints are considered in problem (9), it can imagine that drastic shift in Eq. (10) would be obtained and result in negative bounds on one or some of the constraints for the next iteration, which is an unexpected issue in structural optimization. Therefore, in order to prevent such large shifting value from influencing the convergence of optimization procedure, in this study the reduction coefficients are introduced for calculating the shift of bounds at each iteration, and the quantile-based SORA at (k + 1)th iteration leads to the following deterministic optimization problem

(11)\begin{array}{rl}& \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}V\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right),\\ & \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{g}_j\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A};\theta \right)⩽{\overline{g}}_j-{\epsilon }_j{\overline{c}}_j^{k+1},\\ & j=1,2,\dots ,n_g,\\ & {\gamma }^{\mathrm{c}\mathrm{r}}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A};\theta \right)⩽{\overline{\gamma }}_{\gamma }-{\epsilon }_{\gamma }{\overline{c}}_{{}_{\gamma }}^{k+1},\\ & \underset{\_}{\mathit{t}}⩽\mathit{t}⩽\overline{\mathit{t}};\underset{\_}{\mathit{A}}⩽\mathit{A}⩽\overline{\mathit{A}},\end{array}

where the reduction coefficients {\epsilon }_j and {\epsilon }_{\gamma } are determined as

(12)

Note that for k = 0 (i.e., the initial iteration), there is no solution from the former iteration and thus Eq. (10) is not available. Therefore, the shift of constraint boundaries {\overline{c}}_j^{k+1} and {\overline{c}}_{\gamma }^{k+1}are considered to be 0 in problem (11) at the initial iteration as a start [20,30].

The stopping criteria of the quantile-based SORA are as follows: 1) each reliability constraint is fulfilled; 2) one or more reliability constraint remains equal to prevent a conservative result. These two criteria can be mathematically expressed as follows.

(13)\begin{array}{rl}& Q_{R_j}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{A}}^k;\theta \right)⩽{\overline{\mathrm{g}}}_j,\\ & \mathrm{\forall }j\in \left\{1,2,\dots ,n_g\right\}\mathrm{a}\mathrm{n}\mathrm{d}{Q}_{R_{\lambda }}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{A}}^k;\theta \right)⩽\overline{\gamma },\end{array}

(14)\begin{array}{rl}& Q_{R_{j_e}}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{A}}^k;\theta \right)={\overline{\mathrm{g}}}_{j_e},\\ & \mathrm{\exists }{j}_e\in \left\{1,2,\dots ,n_g\right\}\mathrm{o}\mathrm{r}{Q}_{R_{\lambda }}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{A}}^k;\theta \right)=\overline{\gamma }.\end{array}

3.2 Estimation of quantile

To calculate Eq. (10), one need to calculate the quantiles Q_{R_j}and Q_{R_{\lambda }}. after obtaining the solution t^k and A^k of problem (11) at the kth iteration. Based on Ref. [30], in this study this is achieved by using MEM under the constraint of sample L-moment, and its brief derivation is given below for completeness ofhe paper. See Ref. [30]. for more information.

Take the jth (j = 1,2,…,n_g) displacement constraint as example, and suppose g_j is a continuous random variable under uncertainty, i.e., Z_j^k=g_j\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left({\mathit{t}}^k\right),{\mathit{A}}^k;\theta \right), with unknown cumulative distribution function (CDF) F_j^k\left({\mathrm{z}}_j^k\right) and probability density function (PDF) f_j^k\left({\mathrm{z}}_j^k\right). Let Q_j^k(q) be the quantile function and {Q^{\prime }}_j^k(q) be the corresponding density-quantile function of Z_j^k [44], and H_j^k be the entropy of Z_j^k. By denoting L_{j,r}^k and l_{j,r}^k as the rth order (r ≥ 1) of L-moment and sample L-moment of Z_j^k, respectively, the optimization problem to estimate {Q^{\prime }}_j^k(q) using MEM with n_L sample L-moment constraints is formulated as Ref. [30]:

(15)\begin{array}{rl}& \mathrm{M}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{H}_j^k={\int }_0^1\ln {Q^{\prime }}_j^k\left(q\right)\mathrm{d}q,\\ & \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{L}_{j,r}^k=l_{j,r}^k,r=1,2,\dots ,n_{\mathrm{L}}.\end{array}

Let P_{r-1}^{\ast }\left(q\right) be the (r – 1)th order shifted Legendre polynomial, and define K_r\left(q\right) as

(16)K_r\left(q\right)={\int }_q^1{P}_{r-1}^{\ast }\left(v\right)\mathrm{d}v.

We can calculate the rth order L-moment L_{j,r}^k, as Ref. [30]:

(17)L_{j,r}^k={\int }_0^1{K}_r\left(q\right){Q^{\prime }}_j^k\left(q\right)\mathrm{d}q-{\left[K_r\left(q\right)Q_j^k\left(q\right)\right]}_0^1,

where {\left[K_r\left(q\right)Q_j^k\left(q\right)\right]}_0^1 represents the difference of K_r\left(q\right)Q_j^k\left(q\right) between q = 0 and 1. Then by incorporating Eq. (17) to the constraints in Eq. (15) and denoting h_{j,r}^k=l_{j,r}^k+{\left[K_r\left(q\right)Q_j^k\left(q\right)\right]}_0^1, problem (15) is restated in terms of {Q^{\prime }}_j^k\left(q\right) as

(18)\begin{array}{rl}& \mathrm{M}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}{\int }_0^1\ln {Q^{\prime }}_j^k\left(q\right)\mathrm{d}q,\\ & \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{h}_{j,r}^k={\int }_0^1{K}_r\left(q\right){Q^{\prime }}_j^k\left(q\right)\mathrm{d}q,r=1,2,\cdots ,n_{\mathrm{L}},\end{array}

where the calculation of sample L-moment l_{j,r}^k in h_{j,r}^k is given in Appendix A in Electronic Supplementary materials, and the details of solving problem (18) are also given in Appendix B in Electronic Supplementary materials. After solving problem (18), the value of Q_j^k\left(q\right) is obtained by integrating {Q^{\prime }}_j^k\left(q\right) from q = 0 as

(19)Q_j^k\left(q\right)=Q_j^k\left(0\right)+{\int }_0^q{Q^{\prime }}_j^k\left(u\right)\mathrm{d}u.

Following Ref. [30], in this study Q_j^k\left(0\right) in Eq. (19) is approximated by the smallest order statistic Z_{j,1:m_s}^k (i.e., the smallest value among the samples) [30,45] as given in Appendix A in Electronic Supplementary materials, and the integral in Eq. (19) is evaluated up to q = 0.9999 in the following examples to avoid numerical difficulty. Accordingly, the quantiles Q_{R_j} and Q_{R_{\lambda }} can be directly obtained by letting q in quantile function Eq. (19) equal to the target probability R_j and R_γ in Eq. (5).

3.3 Uncertainty in member stiffness

As stated in Subsection 3.1, the vector of uncertainty θ. includes the random perturbation in nodal locations, cross-sectional areas and Young’s modulus, which are modeled in the same manner as in Ref. [40]. The main steps of modeling these uncertain parameters are given below for completeness of the paper, and details are referred to Ref. [40].

Consider member i which is connected by two red nodes i1 and i2 in Fig. 1(a) with length L_i and global coordinate system (x, y), and introduce random perturbations (Δx_j, Δy_j) to the x- and y-coordinates nodes i1 and i2 in as

(20)\left({x^{\prime }}_j,{y^{\prime }}_j\right)=\left(x_j,y_j\right)+\left(\mathrm{\Delta }{x}_j,\mathrm{\Delta }{y}_j\right)\mathrm{f}\mathrm{o}\mathrm{r}j=i1,i2,

where \left(x_j,y_j\right) and \left(x_j^{\prime },y_j^{\prime }\right) denote the original coordinate and the coordinate after adding perturbations of node j, and they are mked as red and blue dots in Fig. 1(a), respectively. Then, the member length L_i is adjusted to L_i^{\prime }. A single Euler–Bernoulli beam element will be used to simulate the member if L_i^{\prime } is small enough. Otherwise, the member is discretized into four elements of equal length, and these elements are separated by three intermediate nodes, i3, i4, and i5, as illustrated by the red dots in Fig. 1(b). The eccentricity e of the ith member is specified, and the random perturbations at nodes i3, i4, and i5 are constrained to be perpendicular to the orientation of the ith member. Then the random perturbations \left(\mathrm{\Delta }{x}_j,\mathrm{\Delta }{y}_j\right) are added to the three nodes as

(21)\begin{array}{rl}& \left(x_j^{″},y_j^{″}\right)=\left(x_j^{\prime },y_j^{\prime }\right)+\left(\mathrm{\Delta }{x}_j,\mathrm{\Delta }{y}_j\right)\mathrm{f}\mathrm{o}\mathrm{r}j=i3,i4,i5,\\ & \sqrt{\mathrm{\Delta }{x_j}^2+\mathrm{\Delta }{y_j}^2}⩽e{L}_i^{\prime },\left(\mathrm{\Delta }{x}_j,\mathrm{\Delta }{y}_j\right)\mathrm{\perp }\left(x_{i2}^{\prime }-x_{i1}^{\prime },y_{i2}^{\prime }-y_{i1}^{\prime }\right),\end{array}

where \left(x_j^{″},y_j^{″}\right),\left(j=i3,i4,i5\right) are the uncertain locations of nodes i3, i4, and i5 as shown by the blue dots in Fig. 1(b). Moreover, to account for the dependence of uncertain locations between end nodes and intermediate nodes, the following exponential decay function is used to calculate the correlation coefficient [40].

(22)c_{j1j2}=\exp \left(-\frac{\Vert \left(x_{j1},y_{j1}\right)-\left(x_{j2},y_{j2}\right)\Vert }{L_{\mathrm{e}\mathrm{x}\mathrm{p}}}\right),j1,j2=i1,i2,...,i5,

where c_j1j2 represents the correlation coefficient describing the dependence between nodes j1 and j2, and L_exp defines the range over which these uncertainties are likely to be correlated. In a similar way, the random perturbations of cross-sectional area and Young’s modulus are given as

(23)A_i^{\prime }=A_i+\mathrm{\Delta }{A}_i,E_i^{\prime }=E_i+\mathrm{\Delta }{E}_i,

where A_i, A_i^{\prime }, E_i, and E_i^{\prime } are the original and uncertain cross-sectional areas and Young’s modulus of the ith member, respectively. It is worth noting that in Fig. 1(b) all the beam elements have the same A_i^{\prime } and E_i^{\prime }.

4 Penalization on geometrical stiffness and sensitivity analysis

In this section we use penalization method to alleviate the singularity phenomenon, and the corresponding sensitivity analysis is derived for use in a gradient-based algorithm to solve the deterministic optimization problem (11) at each iteration of the quantile-based SORA method.

4.1 Penalization on element geometrical stiffness

It is well known that the singularity in global stability results in superficial local buckling [36,46]. Therefore, to accurately calculate the global geometrical stiffness, a penalization scheme is applied using an approximated Heaviside function F\left(A_i\right) as

(24){\overline{\mathit{K}}}_{Gi}=F\left(A_i\right){\mathit{K}}_{Gi},i=1,2,\dots ,m_e,

where {\mathit{K}}_{Gi} and {\overline{\mathit{K}}}_{Gi} are the original and penalized ith element geometrical stiffness matrices with a shape of 6-by-6, and m_e is the total number of elements. \left(A_i\right) is defined as

(25)F\left(A_i\right)=\frac{1}{2}\left[1+\mathrm{e}\mathrm{r}\mathrm{f}\left(\frac{A_i-\mu }{\sigma \sqrt{2}}\right)\right],

where the error function erf(·) is used to approximate Heaviside function with parameters \mu and \sigma. Equation (25) is continuous function and differentiable, and μ acts as a threshold for the cross-sectional area. Elements with areas below this threshold will have their geometrical stiffness penalized, and \sigma controls the width of the transition zone around \mu. Figure 2 illustrates the shape of F(A_i) in Eq. (25) for a positive value of \mu and different values of \sigma. Since the geometrical stiffness of thin element is expected to be excluded from the global geometrical stiffness to suppress the superficial local buckling, the values of \mu and \sigma are considered small enough in the numerical examples, and their influence on final optimal solutions are also given in Electronic Supplementary materials Appendix E.

4.2 Sensitivity of linear buckling load factor

In this section only the sensitivity coefficients of the linear buckling load are derived below. Those of the structural volume and displacements are referred to Refs. [4,13]. Let {\mathit{\Phi }}^{\mathrm{c}\mathrm{r}} be the buckling shape of {\lambda }^{\mathrm{c}\mathrm{r}}, which satisfies the following normalization condition [4]:

(26){{\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}}^{\mathrm{T}}{\mathit{K}}_{\mathrm{G}}{\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}=1.

Then based on Eqs. (6) and (26), the sensitivity coefficients of {\lambda }^{\mathrm{c}\mathrm{r}} to the design variables t_i and A_i are calculated as Refs. [4,13]:

(27)\begin{array}{rl}\frac{\mathrm{\partial }{\lambda }^{\mathrm{c}\mathrm{r}}}{\mathrm{\partial }{t}_i}=& \sum _{j=1}^{n_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}}\left(\frac{\mathrm{\partial }{\lambda }^{\mathrm{c}\mathrm{r}}}{\mathrm{\partial }{x}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}\frac{\mathrm{\partial }{x}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}{\mathrm{\partial }{t}_i}+\frac{\mathrm{\partial }{\lambda }^{\mathrm{c}\mathrm{r}}}{\mathrm{\partial }{y}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}\frac{\mathrm{\partial }{y}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}{\mathrm{\partial }{t}_i}\right),\\ & i=1,2,\dots ,m,\\ \frac{\mathrm{\partial }{\lambda }^{\mathrm{c}\mathrm{r}}}{\mathrm{\partial }{A}_i}=& {{\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}}^{\mathrm{T}}\frac{\mathrm{\partial }\mathit{K}}{\mathrm{\partial }{A}_i}{\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}\\ & +{\lambda }^{\mathrm{c}\mathrm{r}}\sum _{k=1}^{m_{\mathrm{e}}}{{\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}}^{\mathrm{T}}\left(\frac{\mathrm{\partial }F\left(A_{k_{\mathrm{e}}}\right)}{\mathrm{\partial }{A}_i}{\mathit{K}}_{\mathrm{G}k}+\frac{\mathrm{\partial }{\overline{\mathit{K}}}_{\mathrm{G}k}}{\mathrm{\partial }{N}_k}\frac{\mathrm{\partial }{N}_k}{\mathrm{\partial }{A}_i}\right){\mathit{\Phi }}^{\mathrm{c}\mathrm{r}},\\ & i=1,2,\dots ,m,\end{array}

where N_k is the axial force of element k, and x_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j} and y_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j} are the jth entries of {\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}} and {\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}, respectively. Note that in Eq. (27) the sizes of {\mathit{K}}_{\mathrm{G}k} and {\overline{\mathit{K}}}_{\mathrm{G}k} are n_dof -by-n_dof, i.e., in a global form, and their nonzero entries corresponds to the DOFs of element k and are obtained from Eq. (24), where the other entries are all zero. By denoting the length of element k as l_k, the derivatives \mathrm{\partial }{\lambda }^{\mathrm{c}\mathrm{r}}/\mathrm{\partial }{x}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j} and \mathrm{\partial }{\lambda }^{\mathrm{c}\mathrm{r}}/\mathrm{\partial }{y}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j} in Eq. (27) are calculated by Refs. [4,13]:

(28)\begin{array}{rl}\frac{\mathrm{\partial }{\lambda }^{\mathrm{c}\mathrm{r}}}{\mathrm{\partial }{x}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}=& {{\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}}^{\mathrm{T}}\frac{\mathrm{\partial }\mathit{K}}{\mathrm{\partial }{x}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}{\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}\\ & +{\lambda }^{\mathrm{c}\mathrm{r}}\sum _{k=1}^{m_e}{{\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}}^{\mathrm{T}}\left(\frac{\mathrm{\partial }{\overline{\mathit{K}}}_{\mathrm{G}k}}{\mathrm{\partial }{l}_k}\frac{\mathrm{\partial }{l}_k}{\mathrm{\partial }{x}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}+\frac{\mathrm{\partial }{\overline{\mathit{K}}}_{\mathrm{G}k}}{\mathrm{\partial }{N}_k}\frac{\mathrm{\partial }{N}_k}{\mathrm{\partial }{x}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}\right){\mathit{\Phi }}^{\mathrm{c}\mathrm{r}},\\ \frac{\mathrm{\partial }{\lambda }^{\mathrm{c}\mathrm{r}}}{\mathrm{\partial }{y}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}=& {{\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}}^{\mathrm{T}}\frac{\mathrm{\partial }\mathit{K}}{\mathrm{\partial }{y}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}{\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}\\ & +{\lambda }^{\mathrm{c}\mathrm{r}}\sum _{k=1}^{m_e}{{\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}}^{\mathrm{T}}\left(\frac{\mathrm{\partial }{\overline{\mathit{K}}}_{\mathrm{G}k}}{\mathrm{\partial }{l}_k}\frac{\mathrm{\partial }{l}_k}{\mathrm{\partial }{y}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}+\frac{\mathrm{\partial }{\overline{\mathit{K}}}_{\mathrm{G}k}}{\mathrm{\partial }{N}_k}\frac{\mathrm{\partial }{N}_k}{\mathrm{\partial }{y}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},j}}\right){\mathit{\Phi }}^{\mathrm{c}\mathrm{r}}\end{array}.

Based on Eqs. (27) and (28), we can obtain the sensitivity coefficients of {\gamma }^{\mathrm{c}\mathrm{r}} to t_i and A_i in problem (11) as {\gamma }^{\mathrm{c}\mathrm{r}}=\left(1/{\lambda }^{\mathrm{c}\mathrm{r}}\right).

5 Numerical examples

In this section we explore three numerical examples to demonstrate the effectiveness of the proposed method. The deterministic optimization problem (11) in each iteration of quantile-based SORA is solved using fmincon in MATLAB [47], which is run on a laptop with 13th Gen Intel Core i9-13950HX, and the sequential quadratic programming (SQP) is selected as the optimization algorithm to solve problem (11), which consists of three main stages of updating the Hessian matrix, solving the quadratic programming problem and initialization line search and objective function evaluation. The parameters for SQP are set to default values in fmincon in MATLAB, and more information about SQP in fmincon can refer to Ref. [47]. For comparison, the results of the quantile-based SORA method at both the initial and final iterations are presented in each example, and they are denoted as solutions I and F, respectively. Note that since SQP can only lead to local optimum for nonconvex optimization problem, and the proposed SORA could fall into an infeasible region where the boundary moves back and forth, in the following examples the solution is the first converged result found by running several trials. The sample size of sample L-moments in Eqs. (15) and (18) is m_s = 50 with the first four orders, and to assess the accuracy of quantile estimation by MEM using sample L-moments, the Monte Carlo simulation (MCS) with 1 × 10⁵ samples is employed in numerical examples. The eccentricity e in Eq. (21) is 0.01, which is the same as the value in Ref. [40], and the vector of random perturbations \theta =\left(\mathrm{\Delta }\mathit{x},\mathrm{\Delta }\mathit{y},\mathrm{\Delta }\mathit{A},\mathrm{\Delta }\mathit{E}\right) is assumed to be uniformly distributed, that is, \mathrm{\Delta }\mathit{x}\in \left[\mathrm{\Delta }{\mathit{x}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}},\mathrm{\Delta }{\mathit{x}}_{\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}\right], \mathrm{\Delta }\mathit{y}\in [\mathrm{\Delta }{\mathit{y}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}},\mathrm{\Delta }{\mathit{y}}_{\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}], \mathrm{\Delta }\mathit{A}\in \left[\mathrm{\Delta }{\mathit{A}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}},\mathrm{\Delta }{\mathit{A}}_{\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}\right] and \mathrm{\Delta }\mathit{E}\in \left[\mathrm{\Delta }{\mathit{E}}_{\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}},\mathrm{\Delta }{\mathit{E}}_{\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}\right], with the subscripts “lower” and “upper” representing the lower and upper limits of the uniform distribution. However, it is worth emphasizing that in Section 3 no prior distribution information of the random perturbation is available to formulate the quantile-based SORA problem, and the uniform distribution is only utilized for generating samples in numerical examples [30,48].

The parameter values used in the following numerical examples are summarized in Table 1, and \overline{\mathit{I}} is the unit vector with same dimension as the corresponding parameter. These parameters are determined according to previous research [40,43]. As illustrated in Fig. 1(b), the uncertainty in the locations of end nodes, along with the eccentricity value e, determine the variation range of the intermediate nodes, and the locations of fixed nodes are assumed to be invariant without uncertainty, only the bounds for the random perturbations of free nodes are given and are denoted as \mathrm{\Delta }{\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}}, \mathrm{\Delta }{\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}, \mathrm{\Delta }{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}}, and \mathrm{\Delta }{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e},\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{e}\mathrm{r}}. Each frame member longer than 0.1 m will be divided into four beam elements, while shorter members will be simulated by a single element. Note that in Table 1 the number of nodes n excludes any intermediate nodes within member. For simplicity and completeness of the paper, some of the results are presented in Appendix C in Electronic Supplementary materials, and a brief comparison of the proposed method with surrogate model-based RBDO method is also given in Appendix D in Electronic Supplementary materials.

5.1 Example 1

The first example is an optimization of a simple plane frame with two square grids in the initial structure as shown in Fig. 3. The structure is supported at nodes 1 and 2 by pinned connections, with a 1000 kN downward load applied at node 5. Following Eq. (11), the corresponding deterministic optimization problem for the (k + 1)th iteration of quantile-based SORA is expressed as

(29)\begin{array}{rl}& \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}V\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\mathit{A}\right),\\ & \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{g}_5\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right)⩽5\times {10}^{-3}-{\epsilon }_5{\overline{c}}_5^{k+1},\\ & {\gamma }^{\mathrm{c}\mathrm{r}}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right)⩽0.05-{\epsilon }_{\gamma }{\overline{c}}_{\gamma }^{k+1},\\ & \underset{\_}{\mathit{t}}⩽\mathit{t}⩽\overline{\mathit{t}},\underset{\_}{\mathit{A}}⩽\mathit{A}⩽\overline{\mathit{A}},\end{array}

where g₅ is defined as the vertical displacement of node 5 in the downward direction, and {\overline{c}}_5^{k+1}, {\overline{c}}_{\gamma }^{k+1} and {\epsilon }_5, {\epsilon }_{\gamma } can be obtained by Eqs. (10) and (12), respectively. The optimization process converged after 7 iterations with CPU time 16.8913 s, and the solutions I and F and their corresponding buckling modes are presented in Figs. 4(a) and 4(b) and Figs. 5(a) and 5(b), respectively. For comparison purpose, the following deterministic shape and topology optimization problem without global stability constraint is also solved:

(30)\begin{array}{rl}& \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}V\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}},\mathit{A}\right),\\ & \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{g}_5\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right)⩽5\times {10}^{-3},\\ & \underset{\_}{\mathit{t}}⩽\mathit{t}⩽\overline{\mathit{t}},\underset{\_}{\mathit{A}}⩽\mathit{A}⩽\overline{\mathit{A}}.\end{array}

The optimization result of problem (30) is presented in Fig. 4(c), and its buckling mode is displayed in Fig. 5(c). Table 2 lists the structural volume, linear buckling load factor and quantiles of solutions I and F with global stability constraints and the result of problem (30), and the numbers in parentheses are the quantiles obtained through MCS. Appendix C in Electronic Supplementary materials provides detailed information of nodal locations, force densities t, cross-sectional areas A and member length L.

It can be seen from Fig. 4 and Table 2 that the result of problem (30) has a different structural shape from solutions I and F with 4.7% and 14% smaller structural volume, respectively, while the linear buckling load factors of solutions I and F are about 5 times as large as that of the result of problem (30). By comparing Fig. 5(c) with Figs. 5(a) and 5(b), it can be found that the buckling modes of solutions I and F are similar, while the global stability of the solution without stability constraint in Fig. 5(c) will be significantly deteriorated by the buckling of long lower chord member. This indicates that when global stability constraint is considered, shape optimization can result in an effective improvement without increasing too much on the structural volume. Moreover, from Table 2 it can be observed that the quantile Q_{R_{\lambda }} of the result of problem (30) decreases about 9% from {\lambda }^{\mathrm{c}\mathrm{r}}, whereas the quantiles Q_{R_{\lambda }} of solution I and F decrease about 7.3% and 7.5% from their corresponding {\lambda }^{\mathrm{c}\mathrm{r}}, respectively, demonstrating that the decrease of linear buckling load factor is also enhanced when global stability constraint is imposed.

On the other hand, although solutions I and F have similar geometries and topologies, to satisfy the reliability requirements, the cross-sectional areas of the members in solution F are larger than those in solution I, resulting in an increased structural volume. It is interesting to note in Table 2 that the difference between {\gamma }^{\mathrm{c}\mathrm{r}} and quantile Q_{R_{\gamma }} for solution I is approximately equal to that of solution F, and both quantiles of constraints are approximately on the boundaries, showing the effectiveness of the proposed method to shift the upper bounds of constraints toward feasible domain of the original RBDO problem.

5.2 Example 2

The second example considers a cantilever beam with pinned supports at nodes 1, 2, and 3, as illustrated in Fig. 6. The beam consists of a 3 × 2 grid, comprising 12 nodes and 27 members, and is subjected to a vertical downward load of 1000 kN at node 11. Accordingly, nodes 1, 2, 3, and 11 are selected as fixed nodes for the FDM. Following Eq. (11), a displacement constraint is enforced at node 11, leading to the formulation of the (k + 1)th optimization problem as:

(31)\begin{array}{rl}& \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}V\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right),\\ & \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{g}_{11}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right)⩽3\times {10}^{-3}-{\epsilon }_{11}{\overline{c}}_{{}_{11}}^{k+1},\\ & {\gamma }^{\mathrm{c}\mathrm{r}}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right)⩽0.05-{\epsilon }_{\gamma }{\overline{c}}_{\gamma }^{k+1},\\ & \underset{\_}{\mathit{t}}⩽\mathit{t}⩽\overline{\mathit{t}},\underset{\_}{\mathit{A}}⩽\mathit{A}⩽\overline{\mathit{A}}.\end{array}

The optimization procedure terminated after 4 iterations with CPU time 794.7956 s, and Fig. 7 illustrates the results from the initial and final iterations. Table 3 summarizes the structural volumes and corresponding quantiles Q_{R_{11}}and Q_{R_{\gamma }}for solutions I and F, and the numbers in parentheses show quantiles from MCS. The corresponding nodal coordinates, force densities, cross-sectional areas and lengths of members are given in Tables C3 and C4 in Electronic Supplementary materials Appendix C. As we can see from Table 3 the quantile Q_{R_{\gamma }} of global stability decreases about 33% from solution I to solution F, however, the displacement quantile Q_{R_{11}} decreases only about 7% from solution I to solution F, indicating that the improvement of nodal displacement constraint will contribute to the improvement of global structural stability. Reference [32] also explores the relation between nodal and global stabilities, and it is found in the next numerical example where only the global stability constraints are active in the solution F.

Furthermore, Fig. 8 displays the quantiles of global stability and vertical displacement at node 11, as well as the counterparts obtained by MCS, and the upper limits for both constraints are indicated by the black solid and dashed horizontal lines, respectively. It can be seen from Fig. 8 and Table 3 that for solution I both the probabilities of exceeding the limits {\overline{g}}_{11}=3\times {10}^{-3}\mathrm{m} and \overline{\gamma }=0.5 are over 0.99. Compared to the numerical example 2 in Ref. [30], the probability of displacement constraint at node 11 exceeding the upper limit is higher in this case. This is because in this study the uncertainties are incorporated in the locations of all nodes, not just the two end nodes (as discussed in Subsection 3.2), introducing a greater degree of uncertainty into the frame member to deviate from a perfectly straight shape. However, solution F achieves the target reliability levels for both constraints, showing the effectiveness of the proposed method on handling more uncertainties in structural optimization.

5.3 Example 3

The last example is a shape and topology optimization of bridge frame with 14 nodes and 31 members, as shown in Fig. 9. The bridge frame is supported by a pin at node 1 and a roller at node 13. Five downward vertical loads of 1000 kN each are applied at nodes 3, 5, 7, 9, and 11. The optimization objective is to minimize the structural volume subject to displacement constraints on these 5 loading nodes and global stability, and optimization problem at (k + 1)th iteration is written as:

(32)\begin{array}{rl}& \mathrm{M}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}V\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right),\\ & \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}{g}_j\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right)⩽3\times {10}^{-3}-{\epsilon }_j{\overline{c}}_{{}_j}^{k+1},\\ & j=3,5,7,9,11,\\ & {\gamma }^{\mathrm{c}\mathrm{r}}\left({\mathit{x}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),{\mathit{y}}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\mathit{t}\right),\mathit{A}\right)⩽0.05-{\epsilon }_{\gamma }{\overline{c}}_{{}_{\gamma }}^{k+1},\\ & \underset{\_}{\mathit{t}}⩽\mathit{t}⩽\overline{\mathit{t}},\underset{\_}{\mathit{A}}⩽\mathit{A}⩽\overline{\mathit{A}},\end{array}

The number of optimization iterations is k = 59 with CPU time 7662.6191 s, and solutions I and F are displayed on Fig. 10. The structural volume and quantiles Q_{R_j}\left(j=3,5,7,9,11\right)and Q_{R_{\gamma }}associated with solutions I and F are summarized in Table 4, where the values in parentheses denote the corresponding quantiles obtained from MCS. Figure 11 displays the displacement quantile functions of nodes 3, 5, 7, 9, and 11 and stability quantile function of solutions I and F, and the limits for displacement and stability constraints are represented by black solid and dashed lines, respectively.

Table 4 summarizes that in solution I the displacement quantiles of five loading nodes are greater than the limit 3 × 10⁻³, while Q_{R_{\gamma }} slightly exceeds the limit 0.05. Besides, in Fig. 10(a) due to its connection to three thick members (12, 16, and 17), node 7 experiences less influence from the uncertainties, and consequently, it exhibits the smallest quantile Q_{R_7} among the five loading nodes. Meanwhile, solution F in Fig. 10(b) exhibits a different configuration compared to solution I, and he displacement constraint at node 7 is the only active constraint, while Q_{R_7} is significantly lower than the limit of 0.05. These observations suggest that placing displacement constraints on specific nodes within the frame can indeed improve global stability to a certain degree, resulting a higher reliability for the structure’s global stability compared to the individual displacements of nodes 3, 5, 7, 9, and 11, and all constraints are satisfied except for the displacement constraint at node 7, which holds true as an equality constraint.

Furthermore, Fig. 12 shows the convergence histories of the structural volume, maximum value of five displacement quantiles Q_{R_j}\left(j=3,5,7,9,11\right) and global stability quantile Q_{R_{\gamma }}. The maximum value among displacement quantiles Q_{R_j} is greater than 6 × 10⁻³ m at first iteration. Consequently, Eq. (12) is employed to assign a value of 0.5 to the reduction coefficient ε_j in Eq. (32) to avoid abrupt changes in the constraint boundary. In a similar manner, the global stability quantile rises to approximately 0.17 at 37th iteration, and this triggers a shift in the constraint boundary that exceeds twice the predefined limit of 0.05. To ensure a smooth optimization process and prevent drastic changes, Eq. (12) is then used to set the reduction coefficient ε_γ in Eq. (32) to 0.1. Furthermore, the result at the final iteration achieves a smaller structural volume compared to the solution at k = 1, while it still satisfies all the reliability constraints. This demonstrates the stopping criteria in Eqs. (13) and (14) are effective in preventing excessively conservative solutions.

6 Conclusions

This study introduces a novel reliability-based method for optimizing the shape and topology of plane frames with a global stability constraint. By decoupling the optimization and reliability analysis procedures with robust sample statistics of L-moments, the proposed method can be applied to any types of random parameters, which is different from structural reliability calculation in standards where random parameters are usually assumed to follow normal distribution. By incorporating two stopping criteria an overly conservative solution can be prevented, and the convergence of the SORA optimization process is enhanced by introducing a reduction coefficient. Moreover, to avoid existence of the pseudo buckling mode, a penalization scheme is applied on the element geometrical stiffness matrix to exclude negative contributions of thin elements into global geometrical stiffness matrix.

The effectiveness of the proposed method is demonstrated through three numerical examples. For each example, a sample size of 50 is used to calculate the sample L-moments required in the reliability analysis. Since in practice the distribution of uncertain parameters is usually not known, in all three examples the random parameters are generated without giving prior assumption on their distribution. In the first example it is shown that the upper limits of both displacement and global stability constraints can be effectively shifted toward the feasible domain, and the final solution can be obtained with both constraints approximately on the boundaries. Compared to simply increasing the cross-sectional areas to meet stability requirements, the proposed method can lead to a more effective and stable structure with a smaller overall structural volume. In the second example it is found that the improvement of nodal displacement constraint will indirectly contribute to the improvement of global structural stability, resulting in a greater decrease in global stability quantile than that of displacement quantile in the final solution F. Moreover, in example 3, although all the displacement constraints and global stability constraints are satisfied at the final iteration, the structure is more stable than is required, suggesting that adding displacement constraints at certain points in the structure can improve its overall stability. Since the main computational efforts comes from continually recalculating the linear buckling load factor and displacement for updated shapes and topologies, artificial intelligence, such as operator learning or deep energy method [49], can be utilized as high efficient solvers to enhance the efficiency of the proposed method, which will provide the scope of our future studies.

It is also interesting to note that in example 3 the iteration history of optimization procedure becomes oscillatory at the 37th iteration due to the sudden large jump of the quantile function of global stability. However, it still converges on the final feasible solution with the help of reduction coefficient, which alleviates the drastic shift of constraint boundary in SORA and thus contributes to the convergence of optimization process. Besides, the excessively conservative result is also avoided by using the two stopping criteria, and the reliability of 5 nodal displacements are all close to the requirements, showing that the proposed method can find efficient and reliable optimal structure under different kinds of reliability constraints.

References

[1] Stolpe M . Truss optimization with discrete design variables: A critical review. Structural and Multidisciplinary Optimization, 2016, 53(2): 349–374

[2] Arora J S , Wang Q . Review of formulations for structural and mechanical system optimization. Structural and Multidisciplinary Optimization, 2005, 30(4): 251–272

[3] BendsøeM PSigmundO. Topology Optimization: Theory, Methods, and Applications. Berlin: Springer, 2003

[4] OhsakiM. Optimization of Finite Dimensional Structures. Boca Raton, FL: CRC Press, 2010

[5] Zegard T , Paulino G H . GRAND-Ground structure based topology optimization for arbitrary 2D domains using MATLAB. Structural and Multidisciplinary Optimization, 2014, 50(5): 861–882

[6] Gil L , Andreu A . Shape and cross-section optimization of a truss structure. Computers & Structures, 2001, 79(7): 681–689

[7] Wang D , Zhang W H , Jiang J S . Truss shape optimization with multiple displacement constraints. Computer Methods in Applied Mechanics and Engineering, 2002, 191(33): 3597–3612

[8] Achtziger W . On simultaneous optimization of truss geometry and topology. Structural and Multidisciplinary Optimization, 2007, 33(4-5): 285–304

[9] Ohsaki M . Simultaneous optimization of topology and geometry of a regular plane truss. Computers & Structures, 1998, 66(1): 69–77

[10] Guo X , Liu W , Li H Y . Simultaneous shape and topology optimization of truss under local and global stability constraints. Acta Mechanica Solida Sinica, 2003, 16(2): 95–101

[11] Ohsaki M , Hayashi K . Force density method for simultaneous optimization of geometry and topology of trusses. Structural and Multidisciplinary Optimization, 2017, 56(5): 1157–1168

[12] Hayashi K , Ohsaki M . FDMopt: Force density method for optimal geometry and topology of trusses. Advances in Engineering Software, 2019, 133: 12–19

[13] ShenWOhsakiM. Geometry and Topology Optimization of Plane Frames for Compliance Minimization Using Force Density Method for Geometry Model. London: Springer, 2020

[14] Schuëller G I , Jensen H A . Computational methods in optimization considering uncertainties—An overview. Computer Methods in Applied Mechanics and Engineering, 2008, 198(1): 2–13

[15] Ben-TalA. Laurent El Ghaoui, Nemirovski A. Robust Optimization. Princeton University Press, 2009

[16] ElishakoffIOhsakiM. Optimization and Anti-optimization of Structures Under Uncertainty. London: World Scientific, 2010

[17] Frangopol D M . Structural optimization using reliability concepts. Journal of Structural Engineering, 1985, 111(11): 2288–2301

[18] ChoiS KGrandhiR VCanfieldR A. Reliability-based Structural Design. London: Springer, 2007

[19] Tu J , Choi K K , Park Y H . A new study on reliability-based design optimization. Journal of Mechanical Design, 1999, 121(4): 557–564

[20] Du X , Chen W . Sequential optimization and reliability assessment method for efficient probabilistic design. Journal of Mechanical Design, 2004, 126(2): 225–233

[21] HaoPWangY. A new reliability-based design optimization method with multiple-design points using the active learning kriging. In: Proceedings of Asian Congress of Structural and Multidisciplinary Optimization 2020. Seoul: ASSMO, 2020

[22] Torii A J , Lopez R H , Leandro L F . A general RBDO decoupling approach for different reliability analysis methods. Structural and Multidisciplinary Optimization, 2016, 54: 317–332

[23] Moustapha M , Sudret B , Bourinet J M , Guillaume B . Quantile-based optimization under uncertainties using adaptive Kriging surrogate models. Structural and Multidisciplinary Optimization, 2016, 54(6): 1403–1421

[24] Valdebenito M A , Schuëller G I . A survey on approaches for reliability-based optimization. Structural and Multidisciplinary Optimization, 2010, 42(5): 645–663

[25] Li G , Yang H , Zhao G . A new efficient decoupled reliability-based design optimization method with quantiles. Structural and Multidisciplinary Optimization, 2020, 61(2): 635–647

[26] Ghasemi H , Kerfriden P , Bordas S P A , Muthu J , Zi G , Rabczuk T . Probabilistic multiconstraints optimization of cooling channels in ceramic matrix composites. Composites. Part B, Engineering, 2015, 81: 107–119

[27] Ghasemi H , Brighenti R , Zhuang X , Muthu J , Rabczuk T . Optimal fiber content and distribution in fiber-reinforced solids using a reliability and NURBS based sequential optimization approach. Structural and Multidisciplinary Optimization, 2015, 51(1): 99–112

[28] Kanno Y . A data-driven approach to non-parametric reliability-based design optimization of structures with uncertain load. Structural and Multidisciplinary Optimization, 2019, 60(1): 83–97

[29] Pandey M D . A direct approach to the estimation of quantile function using the maximum entropy principle. Structural Safety, 2000, 22(1): 61–79

[30] Shen W , Ohsaki M , Yamakawa M . Quantile-based sequential optimization and reliability assessment for shape and topology optimization of plane frames using L-moments. Structural Safety, 2022, 94: 102153

[31] Zhou M , Rozvany G I N . Difficulties in truss topology optimization with stress, local buckling and system stability constraints. Structural Optimization, 1996, 11(2): 213–217

[32] OhsakiMIkedaK. Stability and Optimization of Structures: Generalized Sensitivity Analysis. Sendai: Springer, 2010

[33] Descamps B , Filomeno Coelho R . The nominal force method for truss geometry and topology optimization incorporating stability considerations. International Journal of Solids and Structures, 2014, 51(13): 2390–2399

[34] Beghini L L , Beghini A , Baker W F , Paulino G H . Integrated discrete/continuum topology optimization framework for stiffness or global stability of high-rise buildings. Journal of Structural Engineering, 2015, 141(8): 04014207

[35] Ohsaki M . Design sensitivity analysis and optimization for nonlinear buckling of finite-dimensional elastic conservative structures. Computer Methods in Applied Mechanics and Engineering, 2005, 194(30–33): 3331–3358

[36] Tugilimana A , Filomeno Coelho R , Thrall A P . Including global stability in truss layout optimization for the conceptual design of large-scale applications. Structural and Multidisciplinary Optimization, 2018, 57(3): 1213–1232

[37] Torii A J , Lopez R H , Miguel L F F . Modeling of global and local stability in optimization of truss-like structures using frame elements. Structural and Multidisciplinary Optimization, 2015, 51(6): 1187–1198

[38] Ohsaki M , Fujisawa K , Katoh N , Kanno Y . Semi-definite programming for topology optimization of trusses under multiple eigenvalue constraints. Computer Methods in Applied Mechanics and Engineering, 1999, 180(1–2): 203–217

[39] Kanno Y , Ohsaki M , Katoh N . Sequential semidefinite programming for optimization of framed structures under multimodal buckling constraints. International Journal of Structural Stability and Dynamics, 2001, 1(4): 585–602

[40] Shen W , Ohsaki M , Yamakawa M . Robust geometry and topology optimization of plane frames using order statistics and force density method with global stability constraint. International Journal for Numerical Methods in Engineering, 2021, 122(14): 3653–3677

[41] Zhang J Y , Ohsaki M . Adaptive force density method for form-finding problem of tensegrity structures. International Journal of Solids and Structures, 2006, 43(18–19): 5658–5673

[42] Kanno Y , Ohsaki M . Minimum principle of complementary energy of cable networks by using second-order cone programming. International Journal of Solids and Structures, 2003, 40(17): 4437–4460

[43] Shen W , Ohsaki M , Yamakawa M . Multiobjective robust shape and topology optimization of plane frames using order statistics. Structural and Multidisciplinary Optimization, 2021, 63(1): 75–94

[44] Hosking J R M . Distributions with maximum entropy subject to constraints on their L-moments or expected order statistics. Journal of Statistical Planning and Inference, 2007, 137(9): 2870–2891

[45] PrescottPArnoldB CBalakrishnanNNagarajaH N. A First Course in Order Statistics. Philadelphia, PA: Society for Industrial and Applied Mathematics, 1993

[46] Evgrafov A . On globally stable singular truss topologies. Structural and Multidisciplinary Optimization, 2005, 29(3): 170–177

[47] Mathworks. Optimization Toolbox User’s Guide R2018a, 2018

[48] Rocchetta R , Crespo L G , Kenny S P . A scenario optimization approach to reliability-based design. Reliability Engineering & System Safety, 2020, 196: 106755

[49] WangYBaiJLinZWangQAnitescuCSunJ. Artificial intelligence for partial differential equations in computational mechanics. 2024, arXiv:2410.19843

RIGHTS & PERMISSIONS

The Author(s). This article is published with open access at link.springer.com and journal.hep.com.cn