Introduction

Suspension bridges are among the most advantageous bridge types for long spans exceeding a kilometer, and their structural performance is highly sensitive to layout and dimension parameters1. The design of most suspension bridges follows a traditional approach, where designers propose a few design alternatives on the basis of their experience for comparison during project planning2. The feasibility of these solutions is often influenced by subjective factors rather than the objective conditions of the project. Therefore, determining the optimal layout and dimension parameters of suspension bridges to enhance structural performance and reduce costs has always been a key focus for designers.

Traditional engineering design optimization is usually based on a series of parameter or sensitivity analyses to evaluate the structural response to different design variables3. However, when applied to suspension bridges, this approach makes it difficult to determine the global optimal solution. The development of modern optimization theory has addressed these challenges and has been applied in various ways to cable-stayed bridges and suspension bridges, including size and geometric optimization4,5,6,7,8, cable tension optimization9,10,11,12 and dynamic performance optimization under seismic conditions13,14,15,16.

Metaheuristic optimization algorithms have also achieved significant research results in the structural optimization of suspension bridges. Nieto et al.17 proposed a gradient-based method to optimize the bridge section of suspension bridges. Cao et al.4 presented an enhanced PSO and coupled modeling method to optimize the layout and size parameters of suspension bridges. Chen et al.18 introduced an improved PSO and applied it to shape-finding analysis during the installation of suspension bridges. Wei et al.19 used PSO to solve for the final cable configuration of the suspension bridge main cables. These studies focus primarily on single-objective optimization. However, as large and complex structures, suspension bridges require consideration of both cost and structural performance in structural optimization, necessitating multi-objective optimization in their design. Furthermore, existing structural optimization research on suspension bridges has predominantly used steel box girder suspension bridges as case studies, with a lack of studies focusing on the structural optimization of steel truss suspension bridges.

The modeling of suspension bridges requires shape-finding analysis, where iterative shape-finding is used to determine key parameters such as the initial strain of the main cables and the positions of nodes, which are essential for establishing a FE model of the suspension bridge20,21. Therefore, for the optimization of suspension bridges, a shape-finding method with high computational efficiency and accuracy is crucial. In recent years, global search multi-objective optimization algorithms, such as NSGA-II22,23,24 and MOPSO25,26, have been widely applied in engineering. However, their low convergence speed and tendency to become trapped in local optima significantly limit their application in large-scale structures27. In order to address the shortcomings of multi-objective optimization algorithms, many scholars have conducted extensive research. Wang et al.28 proposed an improved dual-archive algorithm, and experiments showed that the algorithm has good convergence, diversity, and complexity. Kumar et al.29 proposed a multi-objective multi-node optimization algorithm (MOMVO), which is based on two archive concepts that separately focus on convergence and diversity, called MOMVO2arc, and it demonstrates good diversity and convergence. Panagant et al.30 proposed two success history-based adaptive multi-objective differential evolution algorithms (SHAMODE) and success history-based whale optimization adaptive multi-objective differential evolution algorithms (SHAMODE-WO). Comparative results showed that both algorithms rank among the top in truss structure optimization design problems. Wang et al.31 proposed an improved MOEA/D-OMDEA algorithm, and experimental results demonstrated that this algorithm is an advanced multi-objective optimization algorithm. It is evident that dual-archive is an effective method for improving algorithm performance. Therefore, to better perform global optimization for suspension bridges, this paper combines NSGA-II and MOPSO to propose a new improved algorithm.

To address the limitations of traditional shape-finding methods, which cannot automatically perform shape-finding during the optimization process as the design variables change, this study combines the parabolic method with FE-based methods to propose a BP shape-finding method. Additionally, this BP method is coupled with traditional FE modeling to improve the efficiency of FE model establishment. To increase the effectiveness of the optimization results, a crowding-based dynamic crossover and mutation probability is introduced for NSGA-II, and a linear learning factor is proposed for MOPSO. On the basis, NSGA-II and MOPSO are combined to propose a dual-population multi-objective particle swarm optimization (DPMOPSO) algorithm. Furthermore, to address the tendency of both algorithms to get trapped in local optima, a generation distance (GD)-based adaptive mutation strategy is introduced.

Coupled modeling method for suspension bridges based on the BP method

Main cable shape-finding analysis

The parabolic method assumes that the self-weight of the main cable is small compared with that of other permanent loads, such as the stiffening girder, allowing all permanent loads to be simplified as uniformly distributed along the span. Under this assumption, the geometric shape of the main cable in each span forms a parabola32. The analytical method is based on the elastic catenary equation, where the main cable is divided into segments at concentrated loads, and the static equilibrium of the loading points and cable segments is solved through iterative calculations33,34. The FE method uses an assumed FE model to update the node positions and internal tension of cable elements through nonlinear structural analysis11.

The traditional shape-finding methods mentioned above cannot achieve automatic shape-finding of the main cable during the optimization process. Therefore, this paper proposes the BP method. In this approach, a parabola is used to fit the initial shape of the main cable, followed by establishing a finite element model. The bisection method is then applied to update the node positions and internal tension of the cable elements through nonlinear structural analysis.

BP method

As shown in Fig. 1, xi and xj represent the coordinates of the pylon anchorage points, and xc represents the coordinates of the midpoint of the main cable. Thus, the parabolic equation fitting the main cable is given by:

$$y = \frac{4f}{{l^{2} }}x^{2} + \frac{C - 4f}{l}x + y_{i}$$
(1)

where l denotes the span of the suspension bridge, f represents the sag of the main cable at the midspan, C is the vertical distance between the tops of the two pylons, and yi is the vertical coordinate of the left pylon anchorage point.

Fig. 1
figure 1

The parabola of fitting the main cable.

Full size image

As shown in Fig. 2, the proposed shape-finding method divides the suspension bridge into a pylon-girder system (as shown in Fig. 2b) and a cable system (as shown in Fig. 2c), which are connected through the hanger tension and the forces acting at the top of the pylon. The BP shape-finding method involves three steps: 1) calculating the support reactions in the pylon-girder system; 2) fitting the initial shape of the main cable via a parabola; and 3) analyzing the cable system.

Fig. 2
figure 2

Pylon-girder system and cable system of suspension bridge.

Full size image

As shown in Fig. 2b, Ri is the reaction force of the ith vertical support. Therefore, the force Ni acting on the lower anchorage point of the ith hanger is equal to Ri. The concentrated force acting on the ith node of the main cable is equal to the sum of Ni and the self-weight of the ith hanger.

Structural analysis requires predefined geometric parameters for the bridge layout. As shown in Fig. 2a, Lm represents the length of the main span. \({L}_{a}^{l}\) and \({L}_{a}^{r}\) denote the distances from the bridge ends to the anchorages on the left and right sides, respectively. f represents the sag of the main cable at the mid-span. \({d}_{x}^{i}\) is the horizontal distance between two adjacent hangers. C is the vertical distance between the tops of the two bridge towers. The initial main cable shape is fitted via Eq. (1) for shape-finding analysis, and the shape of the girder can be expressed via the following equation:

$$y_{B} = h_{0} + x*g_{l}$$
(2)

where h0 is the height of the stiffening girder and gl is the longitudinal slope of the girder.

Figure 3 shows a schematic diagram of the main cable and girder of the suspension bridge.

Fig. 3
figure 3

Diagram of main cable and main girder of suspension bridge.

Full size image

As shown in Fig. 3, the length of the ith hanger under load can be expressed by the following equation:

$$h^{i} = y_{i} - B_{i}$$
(3)

The unstressed length of the ith hanger, \({h}_{0}^{i}\), can be derived from the following equation:

$$h_{0}^{i} = h^{i} - N_{i} *h^{i} /2E_{h} A_{2}$$
(4)

The self-weight of the ith hanger, \({G}_{h}^{i}\), is given by:

$$G_{h}^{i} = w_{h} h_{0}^{i}$$
(5)

The concentrated force Fi acting on the ith node of the main cable is given by:

$$F_{i} = N_{i} + G_{h}^{i}$$
(6)

where Eh and A2 represent the elastic modulus and cross-sectional area of the hanger, respectively; wh denotes the weight of the hanger per unit length.

After Fi is calculated, the main cable is fitted in ANSYS35 according to Eq. (1), and Fi is applied to the corresponding nodes of the main cable. The bisection method is then used to perform iterative shape-finding for the main cable. In the suspension bridge case presented in this paper, after 20 iterations, the displacement of the midpoint of the main cable is only 0.009 mm, meeting the accuracy requirements.

After the shape-finding is completed, the strain \({\varepsilon }_{i}\) of each segment of the main cable is extracted, and then the unstressed length of the ith segment of the main cable, \({l}_{0}^{i}\), is calculated via the following equation:

$$l_{0}^{i} = [(x_{i + 1} - x_{i} )^{2} + (y_{i + 1} - y_{i} )^{2} ]^{1/2} *(1 - \varepsilon_{i} )$$
(7)

Thus, the unstressed length of the main cable in the mid-span, L, is given by:

$$L = \sum\limits_{i = 1}^{2n} {l_{0}^{i} }$$
(8)

After completing the shape-finding for the main cable in the mid-span, the unstressed length of the side-span main cable is calculated on the basis of the equality of horizontal forces in the side and mid-span cables, thus obtaining the unstressed length of the entire main cable.

Coupling of the BP method with traditional FE Modeling

To perform shape-finding via the BP Method, a FE model of the pylon-girder system of the suspension bridge, as shown in Fig. 2b, should first be established in ANSYS. Then, a static analysis is conducted to extract the vertical support reactions. The initial shape of the main cable is fitted via Eq. (1), and the shape-finding of the suspension bridge’s main cable is carried out according to the method described in Section 2.2.2. After the shape-finding is completed, the nodal coordinates of the main cable and the initial strain of each cable segment are extracted, and a FE model of the suspension bridge is established. FE analysis is then conducted in ANSYS to provide the necessary inputs required during the optimization process.

Multi-objective optimization model for suspension bridges

Design variables

Figure 2a shows the layout and geometric parameters of a single-span steel truss suspension bridge with pylons of unequal heights. In addition to the five geometric parameters described in Section 2.2.2 (Lm, \({L}_{a}^{l}\), \({L}_{a}^{r}\), f, C), Hs and Hb represent the length of the shortest hanger at mid-span and the height of the pylon below the bridge deck, respectively, whereas B denotes the width of the bridge. The side view of the pylons is shown in Fig. 4a, and the simplified cross-sections of the pylons and transverse beams are presented in Fig. 4b,c, respectively. Figure 5 shows the truss section that constitutes the steel truss stiffening girder.

Fig. 4
figure 4

Simplified cross sections of the pylon and cross beam.

Full size image
Fig. 5
figure 5

Truss sections of steel truss stiffened girder.

Full size image

Some geometric parameters of the bridge are determined by practical and site conditions, and therefore remain unchanged during the optimization process4. Thus, for a suspension bridge with a fixed total span, the design variables considered in the optimization include the following:

Pylon: the width w1 and height w2 of the pylon, and the thicknesses of the plates T1 and T2;

Transverse beam: the height w3 of the transverse beam, and the thicknesses of plates T3 and T4;

Stiffening girder: as shown in Fig. 5, the stiffening girder consists of nine different ‘H’-shaped sections, with four design variables for each section. Therefore, the design variables are w4–w21 and T5–T22;

Cable system: the sag-to-span ratio (f/Lm), the area of the main cable A1 and hanger A2.

Multi-objective functions

Most designers focus on single-objective optimization in the design of suspension bridges, with material cost typically serving as the objective function. Although such optimization reduces material costs, it does not include a quantitative analysis of the impact on the structural safety of the bridge. Therefore, this paper conducts multi-objective optimization for suspension bridges, with objective functions consisting of two parts: one is the commonly used material cost function, and the other is a function representing the structural safety of the bridge, as shown below:

$$\left\{ {\begin{array}{*{20}l} {\min f_{1} = \sum\limits_{i = 1}^{N} {C_{i} A_{i} L_{i} } } \hfill \\ {\max f_{2} = \frac{{\sigma_{C}^{cr} }}{{\sigma_{C} }} + \frac{{\sigma_{H}^{cr} }}{{\sigma_{H} }} + \frac{{\sigma_{G}^{cr} }}{{\sigma_{G} }} + \frac{{\sigma_{P}^{cr} }}{{\sigma_{P} }}} \hfill \\ \end{array} } \right.$$
(9)

where, f1 represents the material cost. Ci, Ai and Li denote the unit volume material cost, cross-sectional area, and length of component i, respectively. f2 is the safety factor, where \({\sigma }_{C}^{cr}\), \({\sigma }_{H}^{cr}\), \({\sigma }_{G}^{cr}\), and \({\sigma }_{P}^{cr}\) represent the allowable stresses for the main cable, hanger, girder, and pylon, respectively. \({\sigma }_{C}\), \({\sigma }_{H}\), \({\sigma }_{G}\), and \({\sigma }_{P}\) denote the maximum stresses of the main cable, hanger, girder, and pylon calculated in each FE analysis by ANSYS.

Since evolutionary optimization algorithms always seek to minimize the objective function, Eq. (9) is transformed into the following form:

$$\left\{ {\begin{array}{*{20}l} {\min f_{1} = \sum\limits_{i = 1}^{N} {C_{i} A_{i} L_{i} } } \hfill \\ {\min f_{2} = - \left( {\frac{{\sigma_{C}^{cr} }}{{\sigma_{C} }} + \frac{{\sigma_{H}^{cr} }}{{\sigma_{H} }} + \frac{{\sigma_{G}^{cr} }}{{\sigma_{G} }} + \frac{{\sigma_{P}^{cr} }}{{\sigma_{P} }}} \right)} \hfill \\ \end{array} } \right.$$
(10)

Multi-objective particle swarm optimization algorithm based on the dual-population coevolution strategy (DPMOPSO)

Proposal and improvement of the algorithm

NSGA-II has strong search capabilities but is highly dependent on parameters and sensitive to parameter selection. MOPSO converges quickly but tends to become trapped in local optima, failing to find the global optimal solution. Additionally, NSGA-II and MOPSO differ in their information-sharing mechanisms: NSGA-II transmits information through genetic operators, whereas MOPSO uses the global best particle to guide other particles. This leads each algorithm to have its own advantages and disadvantages. Therefore, this paper combines the two algorithms to leverage their respective strengths and proposes DPMOPSO, with improvements to further enhance performance.

Dual-population coevolution strategy

The dual-population coevolution strategy divides the population into two halves on the basis of the non-dominated sorting of NSGA-II. The upper half, with better Pareto rankings, forms the elite population, which leverages NSGA-II’s strong search capabilities to explore other Pareto solution sets within the region and identify non-dominated solutions. The lower half, with poorer Pareto rankings, adopts the MOPSO algorithm for learning within the nonelite population. To prevent premature convergence caused by the MOPSO algorithm, individuals from the elite population are selected through roulette wheel selection as the global best for MOPSO.

Crowding-based dynamic crossover and mutation probability

Crossover and mutation operations in NSGA-II are essential for the exploration of the solution space by population individuals. However, traditional crossover and mutation operations are often blind. To better control the range of crossover and mutation within the population, a dynamic crossover and mutation probability based on population crowding and the number of evolutionary iterations is adopted to enhance population convergence. The specific improvements are as follows:

$$P_{c} = \left\{ {\begin{array}{*{20}l} {P_{cavg} + \frac{{T_{\max } - t}}{{T_{\max } }}*(P_{c\max } - P_{c\min } ), \, d(i) < d_{avg} (i)} \hfill \\ {P_{cavg} , d(i) = d_{avg} (i)} \hfill \\ {P_{cavg} - \frac{t}{{T_{\max } }}*(P_{c\max } - P_{c\min } ), d(i) > d_{avg} (i)} \hfill \\ \end{array} } \right.$$
(11)
$$P_{m} = \left\{ {\begin{array}{*{20}l} {P_{mavg} + \frac{{T_{\max } - t}}{{T_{\max } }}*(P_{m\max } - P_{m\min } ), \, d(i) < d_{avg} (i)} \hfill \\ {P_{mavg} , d(i) = d_{avg} (i)} \hfill \\ {P_{mavg} - \frac{t}{{T_{\max } }}*(P_{m\max } - P_{m\min } ), d(i) > d_{avg} (i)} \hfill \\ \end{array} } \right.$$
(12)

Here, \({P}_{cavg}\), \({P}_{cmax}\), \({P}_{cmin}\), \({P}_{mavg}\), \({P}_{mmax}\), and \({P}_{mmin}\) represent the average, maximum, and minimum values of the crossover probability Pc and mutation probability Pmrespectively. t is the current iteration number, and \({T}_{max}\) is the maximum number of iterations. d(i) denotes the crowding distance of the ith individual, and \({d}_{avg}(i)\) is the average crowding distance of the front to which the ith individual belongs.

Linear learning factor

Particle swarm optimization (PSO) was proposed by Kennedy and Eberhart in 199536. In 2002, Coello et al.37 applied the PSO algorithm to solve multi-objective optimization problems and proposed the multi-objective particle swarm optimization (MOPSO) algorithm. The position update equation and the velocity update equation are shown as follows:

$$v_{id} (t + 1) = wv_{id} (t) + c_{1} r_{1} (pbest_{id} (t) - x_{id} (t)) + c_{2} r_{2} (gbest_{id} (t) - x_{id} (t))$$
(13)
$$x_{id} (t + 1) = x_{id} (t) + v_{id} (t + 1)$$
(14)

where \({v}_{id}\) represents the velocity of the ith particle in the dth dimension; w is the inertia weight; \({pbest}_{id}\) denotes the personal best of the ith particle in the dth dimension, and \({gbest}_{id}\) represents the global best of the ith particle in the dth dimension; and \({x}_{id}\) is the position of the ith particle in the dth dimension. Additionally, \({c}_{1}\) and \({c}_{2}\) are two learning factors, and \({r}_{1}\) and \({r}_{2}\) are random numbers that are uniformly distributed between 0 and 1.

As shown in Eq. (13), the standard MOPSO uses fixed learning factors, which are not conducive to a global search for the optimal solution. Clerc and Kennedy38 pointed out that during the early stages of the algorithm, the individual learning factor \({c}_{1}\) should be set to a larger value, while the social learning factor \({c}_{2}\) should be set to a smaller value. This helps to expand the search range and increase the diversity of the population. In the later stages of the algorithm, when the number of iterations is large, \({c}_{1}\) should be smaller, and \({c}_{2}\) should be larger, which promotes local search and accelerates the convergence speed of the algorithm. On this basis, the learning factors are adjusted, and the improved iteration formula is as follows:

$$\begin{array}{*{20}l} {v_{i} (k + 1) = wv_{i} (k) + c_{1} r_{1} (p_{i} (k) - x_{i} (k)) + c_{2} r_{2} (x_{NSGA - IIr} (k) - x_{i} (k))} \hfill \\ {c_{1} = c_{1\max } - (c_{1\max } - c_{1\min } )*\frac{t}{{T_{\max } }}} \hfill \\ {c_{2} = c_{2\min } + (c_{2\max } - c_{2\min } )*\frac{t}{{T_{\max } }}} \hfill \\ {x_{i} (k + 1) = x_{i} (k) + v_{i} (k + 1)} \hfill \\ \end{array}$$
(15)

In this formula, \({c}_{1max}\), \({c}_{1min}\), \({c}_{2max}\), and \({c}_{2min}\) are the maximum and minimum values of the learning factors \({c}_{1}\) and \({c}_{2}\), respectively. \({x}_{NSGA-IIr}\) is an individual randomly selected from the NSGA-II population.

Generation distance (GD)-based adaptive mutation strategy

Non-dominated sorting is performed on individuals obtained after NSGA-II’s genetic operations and MOPSO-guided optimization, with individuals ranked as non-dominated level 1 containing the current optimal evolutionary information. However, both NSGA-II and MOPSO are prone to getting trapped in local optima, potentially leading to many similar individuals after several generations, which can prevent the algorithm from finding the optimal Pareto front. Therefore, mutation operations are applied to this subset of solutions to increase the ability of the population to escape local optima. Additionally, the dynamic changes in the group of level 1 non-dominated individuals are also important to monitor. GD is used to measure these changes, and the mutation scale is dynamically adjusted accordingly.

Polynomial mutation is used for mutation, where one dimension of the individual is randomly selected for mutation:

$$x_{new}^{d} = x^{d} + (x_{U}^{d} - x_{L}^{d} )\delta *k$$
(16)

where \({x}_{U}^{d}\) and \({x}_{L}^{d}\) are the upper and lower bounds of the dth dimension of the individual, \(\delta\) is the disturbance term, and k is the disturbance ratio coefficient. The disturbance term \(\delta\) is calculated as follows:

$$\delta = \left\{ {\begin{array}{*{20}c} {(2r)^{{\frac{1}{\mu + 1} - 1}} ,\;\quad \;\;} & {r < 0.5} \\ {1 - [2(1 - r)]^{{\frac{1}{\mu + 1}}} ,} & {r \ge 0.5} \\ \end{array} } \right.$$
(17)

where r is a random number between 0 and 1, and \(\mu\) is the mutation distribution index.

GD represents the Euclidean distance between the solutions in the evaluated solution set and the true Pareto front. To measure the dynamic change in the level 1 non-dominated group, the generation distance between the current generation and the previous generation is denoted as \({GD}^{*}\). The difference in \({GD}^{*}\) between two consecutive generations is denoted as \({\Delta GD}^{*}\), and is given by:

$$\Delta GD^{*} (t) = GD^{*} (t) - GD^{*} (t - 1),\;t \ge 4$$
(18)

\({\Delta GD}^{*}\) can measure the convergence rate of the algorithm. When \({\Delta GD}^{*}>0\), the algorithm converges quickly, so the mutation scale should be reduced to improve algorithm efficiency. When \({\Delta GD}^{*}<0\), the mutation scale needs to be increased to enhance the convergence and diversity of the algorithm. The adjustment of the mutation scale is as follows:

$$Q(t + 1) = \left\{ {\begin{array}{*{20}c} {round(0.75Q(t)),} & {\Delta GD^{*} (t) > 0} \\ {Q(t),\quad \quad \quad \quad } & {\Delta GD^{*} (t) = 0} \\ {round(1.25Q(t)),} & {\Delta GD^{*} (t) < 0} \\ \end{array} } \right.$$
(19)

where, round() denotes the rounding function, and Q(t) represents the mutation scale of the tth generation. When t = 1, 2, 3, Q(t) is equal to the current number of individuals with a non-dominated rank of 1. The flowchart of the proposed algorithm is shown in Fig. 6.

Fig. 6
figure 6

Flowchart of the proposed algorithm.

Full size image

Performance evaluation of DPMOPSO

The ZDT1-4 and ZDT6 test functions from the ZDT series39 are selected to evaluate the performance of the proposed optimization algorithm. The ZDT5 function, which describes a deceptive problem, differs from the other test functions and is therefore not used in this study. Additionally, MOPSO, NSGA-II40, SMPSO41, MPSO/D42 and NSGA-III43 are selected as comparison algorithms. The inverted generational distance (IGD)44, hypervolume (HV)45, Spacing46 and Coverage47 are used as evaluation metrics for the algorithms. To ensure fairness in the comparison, the maximum number of iterations for each algorithm is set to 1000, and the population size is set to 200. The source code of the comparison algorithms used in this paper is provided by PlatEMO48.

Tables 1, 2, and 3 show the IGD, HV, and Spacing values of DPMOPSO and the comparison algorithms on the ZDT series test functions, with the optimal values highlighted in bold. From the three tables, it can be observed that DPMOPSO performs the best, achieving the optimal value for each performance evaluation metric, demonstrating the superiority of the proposed algorithm. Table 4 shows the Coverage values between the results obtained by DPMOPSO and the comparison algorithms. From the table, it can be seen that the solutions obtained by DPMOPSO completely cover those obtained by the comparison algorithms, further proving the superiority of the proposed algorithm. To more intuitively demonstrate the superiority of the proposed algorithm, a validation using the Pareto front simulation charts is also provided. Figures 7, 8, 9, 10, and 11 show the Pareto simulation fronts of DPMOPSO and the comparison algorithms on the ZDT series test functions. From these figures, it can be seen that only the results obtained by DPMOPSO coincide with the true Pareto front, while the results of the other comparison algorithms significantly deviate from the true Pareto front. This further confirms the conclusions drawn from Tables 1, 2, 3 and 4: the proposed DPMOPSO has a stronger ability to handle complex problems.

Table 1 IGD values of the DPMOPSO and the comparison algorithms on the test functions.
Full size table
Table 2 HV values of DPMOPSO and the comparison algorithms on the test functions.
Full size table
Table 3 Spacing values of DPMOPSO and the comparison algorithms on the test functions.
Full size table
Table 4 Coverage values of DPMOPSO and the comparison algorithms on the test functions (A \(\succcurlyeq\) B).
Full size table
Fig. 7
figure 7

Simulation results of DPMOPSO and comparison algorithms on ZDT1.

Full size image
Fig. 8
figure 8

Simulation results of DPMOPSO and comparison algorithms on ZDT2.

Full size image
Fig. 9
figure 9

Simulation results of DPMOPSO and comparison algorithms on ZDT3.

Full size image
Fig. 10
figure 10

Simulation results of DPMOPSO and comparison algorithms on ZDT4.

Full size image
Fig. 11
figure 11

Simulation results of DPMOPSO and comparison algorithms on ZDT6.

Full size image

Case study

Description of the suspension bridge and FE model

The feasibility and effectiveness of the proposed optimization method were verified through the optimization design of a single-span steel truss suspension bridge with a total length of 660 m and width of 28 m, accommodating a two-way four-lane roadway. Table 5 lists the materials and properties used for each bridge component. Table 6 lists the fixed geometric and dimensional parameters of the bridge, with a distance of 8 m between adjacent hangers.

Table 5 Materials and properties of bridge components.
Full size table
Table 6 Fixed geometric parameters of the suspension bridge.
Full size table

This example includes the 46 design variables mentioned in section “Design variables”. Table 7 lists the range of each design variable, estimated on the basis of similar bridges. Structural analysis was conducted via the 3D beam-truss model in the FE analysis software ANSYS. The FE model uses beam elements (Beam188) for the main girder, pylon, and transverse beam, and truss elements (Link10) for the cable system. Each FE model consists of 6219 beam elements and 330 truss elements, with a total of 3946 nodes.

Table 7 Range of design variables.
Full size table

Constraints

The constraints include strength and serviceability requirements, as specified below:

Strength requirements: The allowable stresses for the main cable and hanger are 708 MPa and 590 MPa, respectively; the allowable tensile and compressive stresses for the stiffening girder are 115 MPa; and the allowable stresses for the pylon and transverse beam are 32.4 MPa.

Serviceability requirements: The maximum allowable vertical deflection of the girder should be less than \({L}_{m}\)/500, and the maximum displacement of the pylon should not exceed 1/800 of the total pylon height.

Parameter selection of the algorithms

The parameter settings for the proposed optimization algorithm in suspension bridge optimization are as follows: \({N}_{p}={N}_{r}=\) 80, \({c}_{1max}={c}_{2max}=\) 2, \({c}_{1min}={c}_{2min}=\) 0.2, \(w=\) 0.8, \({P}_{cmax}=\) 0.8, \({P}_{cmin}=\) 0.4, \({P}_{mmax}=\) 0.2, \({P}_{mmin}=\) 0.001, \(k=\) 0.1, and \({T}_{max}=\) 20.

In the suspension bridge optimization process, the parameters for MOPSO are set as follows: \({N}_{p}={N}_{r}=\) 80, \({c}_{1}={c}_{2}=\) 1.1, \(w=\) 0.8, and \({T}_{max}=\) 20.

The parameter settings for NSGA-II are set as follows: \({N}_{p}={N}_{r}=\) 80, \({P}_{c}=\) 0.9, \({P}_{m}=\) 0.1, and \({T}_{max}=\) 20.

Figure 12 presents the flowchart of the suspension bridge optimization process using the proposed DPMOPSO and coupled modeling.

Fig. 12
figure 12

Flowchart of suspension bridge optimization program.

Full size image

Results and discussion

All the optimization processes were executed on an Intel Core i5-12490F, 3.00 GHz machine. DPMOPSO was used for multi-objective optimization of the suspension bridge, with MOPSO and NSGA-II as comparison algorithms. Using Spacing, Set Coverage, and Frontier Spread, a comprehensive evaluation was conducted on the solution sets obtained by the three algorithms49,50. Table 8 presents the computational time and the solution set comparisons for each algorithm.

Table 8 Comparison of the solution sets obtained by different algorithms.
Full size table

As shown in Table 8, DPMOPSO requires the longest computation time because it has a greater complexity than do NSGA-II and MOPSO and requires more evaluations of the objective function. Since calculating the objective function requires FE analysis, which is the most time-consuming part of the optimization process, DPMOPSO consequently takes more computation time. DPMOPSO outperforms the other two algorithms in terms of solution quality (distribution uniformity, dominance, and coverage).

Figure 13 shows the solutions obtained by the three algorithms after the optimization of the suspension bridge. This figure shows that although NSGA-II yields more solutions, most of them are dominated by the optimal solutions of the other two algorithms, indicating lower solution quality. The solutions obtained by DPMOPSO can dominate those obtained by the other two algorithms, indicating that its solution quality is the best. This further confirms the robust optimization capability of the proposed DPMOPSO.

Fig. 13
figure 13

The solutions for suspension bridge optimization.

Full size image

Figure 14 shows a trend plot of the two objective functions (f1 and f2) after sorting the solutions obtained by the three algorithms in ascending order of f2. Here, the independent variable is the design scheme index; the dependent variable f1 represents the material cost, measured in yuan; and f2 is the safety factor, which is dimensionless.

Fig. 14
figure 14

The change trend of two objective functions.

Full size image

Figure 14 shows that some design schemes incur higher material costs while achieving lower safety factors. To investigate the reasons behind this phenomenon and facilitate more rational suspension bridge design, the Lasso function51 is used to conduct elastic net regression, analyzing the impact of the selected design variables on the two objective functions. Figure 15 shows the degree of fit between the actual and predicted values of the objective functions (Fig. 15a represents objective function f1, and Fig. 15b represents objective function f2).

Fig. 15
figure 15

The fitting of two objective functions.

Full size image

As shown in Fig. 15, both objective functions are well fitted. The R2 value for f1 is 0.9997, and for f2, it is 0.9544, both greater than 0.9, which confirms the good fit of the objective functions. Thus, the impact of the design variables on the objective functions, as derived from this regression model, is reliable. Fig. 16 shows the influence of partial design variables on the objective functions, and Table 9 presents the 10 design variables with the greatest impact on the objective functions and the magnitude of their influence coefficients.

Fig. 16
figure 16

The degree of influence of partial design variables on the two objective functions.

Full size image
Table 9 Influence coefficients on the objective functions.
Full size table

As shown in Table 9, increasing the sag-to-span ratio (f/Lm) results in the highest increase in material cost. The remaining design variables that have a positive impact on material cost, in descending order of influence coefficients, are the width of the lower chord web of the truss transverse beam (T14), the width of the upper horizonal link web of the truss (T20), and the width of the lower chord flange plate of the truss transverse beam (T13). On the other hand, increasing the width of the upper horizontal flange plate of the truss (T19) actually reduces the material cost of the bridge. For the bridge safety factor (f2), since this objective function is set to its negative value, a positive influence coefficient implies a negative effect on the function, whereas a negative influence coefficient implies a positive effect. Thus, from Table 9, it can be concluded that appropriately increasing the f/Lm is beneficial for enhancing the bridge’s safety factor. Other design variables with a positive impact on the bridge’s safety factor include the width of the lower chord flange plate of the truss transverse beam (T13), the width of the lower horizonal link web of the truss (T22), and the length of the upper chord web of the truss transverse beam (W4). The design variables with a negative impact on f2 include the width of the inner diagonal webs flange plate of the truss transverse beam (T17) and T19.

In summary, for a single-span steel truss suspension bridge with pylons of unequal heights, appropriately increasing the sag-to-span ratio can increase the bridge safety factor. Increasing the width of the lower chord flange plate of the truss transverse beam (T13) is the most economical way to improve the safety factor of the bridge. Although increasing the width of the upper horizontal flange plate of the truss (T19) can reduce the material cost, it also reduces the safety factor of the bridge.

Conclusion

  1. (1)

    The coupled method based on BP shape-finding and the traditional FE modeling strategy meets the requirement of repeatedly generating FE models during the optimization process. Additionally, through shape-finding analysis of the suspension bridge case in this paper, it is demonstrated that the BP shape-finding method is theoretically simple and highly computationally efficient.

  2. (2)

    By combining MOPSO and NSGA-II, a multi-objective particle swarm optimization method based on a dual-population coevolution strategy (DPMOPSO) was proposed. The performance of the proposed algorithm was tested via the ZDT series of test functions. The test results indicate that the proposed DPMOPSO is a more robust multi-objective optimization algorithm.

  3. (3)

    DPMOPSO, MOPSO, and NSGA-II were each applied to the multi-objective optimization design of the suspension bridge case presented in this paper. The solutions obtained from the optimization process reveal that DPMOPSO achieved a solution set with higher quality. However, owing to the extensive improvements made to the proposed algorithm, the computational complexity increased, thereby increasing the computation time.

  4. (4)

    Elastic net regression analysis using the Lasso function on the solutions obtained from the multi-objective optimization of the suspension bridge indicates that appropriately increasing the sag-to-span ratio can increase the bridge’s safety factor. The most economical way to increase the safety factor is to increase the width of the lower chord flange plate of the truss transverse beam. Although increasing the width of the upper horizontal flange plate of the truss can reduce the material cost overall, it simultaneously decreases the safety factor of the bridge.