# Nonpenalty Machine Learning Constraint Handling Using PSO-SVM for Structural Optimization

Authors: Marco M. Rosso; Raffaele Cucuzza; Fabio Di Trapani; Giuseppe C. Marano
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

Firstly formulated to solve unconstrained optimization problems, the common way to solve constrained ones with the metaheuristic particle swarm optimization algorithm (PSO) is represented by adopting some penalty functions. In this paper, a new nonpenalty-based constraint handling approach for PSO is implemented, adopting a supervised classification machine learning method, the support vector machine (SVM). Because of its generality, constraint handling with SVM appears more adaptive both to nonlinear and discontinuous boundary. To preserve the feasibility of the population, a simple bisection algorithm is also implemented. To improve the search performances of the algorithm, a relaxation function of the constraints is also adopted. In the end part of this paper, two numerical literature benchmark examples and two structural examples are discussed. The first structural example refers to a homogeneous constant cross-section simply supported beam. The second one refers to the optimization of a plane simply supported Warren truss beam. The obtained results in terms of objective function demonstrate that this new approach represents a valid alternative to solve constrained optimization problems even in structural optimization field. Furthermore, as demonstrated by the Warren truss beam problem, this new algorithm provides an optimal structural design which represents also a good solution from the technical point of view with a trivial rounding-off that does not jeopardize significantly the optimization design process.


---

## Body

## 1. Introduction

Many optimum design engineering problems can be modelled as continuous nonlinear optimization problems, in which the search space is reduced due to the existence of various constraints:(1)minx∈Ωfx,s.t.gqx≤0∀q=1,…,nq,hrx=0∀r=1,…,nr,in which x=x1,…,xnT is the design variable vector, fx is the objective function (OF) to be minimized, and Ω is a box-type search space, defined as Cartesian products among the admissible intervals within a lower and an upper bound for each design variable Ω=x1l,x1u×⋯×xnl,xnu. The equality constraints can also be rewritten in order to have inequality constraints only, i.e., gpx≤0, for p=1,…,nq,nq+1,…,np, with np=nq+2nr. The modern approaches to solve them are related to the Soft Computing techniques; in particular, the particle swarm optimization (PSO) algorithm [1] is adopted in this paper. The swarm exhibits an intelligent search behavior because of the social and cognitive interactions among the individuals of the population. This mechanism is inspired by natural aggregations such as bird flocking, fish schooling, or swarming of insects in nature. The PSO was successfully adopted in structural optimization fields and one can refer to, e.g., [2–6] for truss problems, to, e.g., [7, 8] for composite structures, to, e.g., [9], for shape optimization of arch dams under earthquake loading, or to, e.g., [10] for optimization of steel frames. The starting positions and velocities are randomly sampled using the Latin Hypercube Sampling (LHS) technique to generate an initial population with minimum correlation between samples [11]. In the Newtonian dynamics-based algorithm, each ith particle (with i=1,…,N, in which N is the population size) is completely defined by the position xik and the velocity vik at the kth generation. The update rules for the next generation are the following:(2)vik+1=vik+c1k+1r1i∗xiPbk−xik+c2k+1r2i∗xiPbk−xik,(3)xik+1=xik+τvik+1τ=1,in which the symbol ∗ denotes the term-by-term vector multiplication (Hadamard product, [12]). Each particle is an intelligent agent with memory of its best visited position until the current iteration in terms of minimal OF xiPbk (denoted as pbest), whereas the term xGbk (denoted as gbest) is usually the best individual of the entire swarm in terms of minimal OF (fully connected neighborhood topology). It is possible to adopt also a local strategy limiting the sharing of information about the gbest to a certain number of neighbor particles [13–15]. The velocity update rule (2) depends on a inertia term related to the previous velocity vector, on a cognitive term related to the pbest attraction, and on a social term related to the gbest, which are randomly weighted by r1ik+1,k+1r2i=rand0,1. The positive scalar acceleration factors c1 and c2 are, respectively, called cognitive and social parameters. To protect the cohesion of the swarm, each velocity component is limited by a maximum value typically assumed as vmax=γxu−xl/τ, where τ=1 is a time-related parameter, whereas γ∈0.1,1 (generally set to 0.5) defines how far a particle can move starting from its current position [13]. After the update phase, the algorithm checks if the particles are inside the box-search space. A stopping criterion is usually set through a maximum number of iterations L; however, the strictly required ones to reach the optimum are not know a priori; therefore, one can refer to the suggestions of [16]. Later, to improve the exploration capacity of the swarm, Shi and Eberhart [17, 18] introduced an inertia weight term wk multiplied to the previous velocity in (2), which changes during iterations from an initial value w0 to a final one wL with a linearly decreasing law. Some other inertia weight laws were proposed during the years [14]. The performance of the algorithm is strongly affected by the choice of the parameters such as the swarm size N, usually set in a range of 20,100 with n≤30, or the acceleration factors, which are usually assumed statically fixed to c1=c2=2 [13]. In this study, it is assumed that all of them are constant values equal to c1=c2=2, w0=0.90, and wL=0.40 [2].Despite several studies (e.g., [19]) demonstrating that this optimizer has a good convergence rate, the PSO, such as many other Soft Computing techniques, was originally proposed to solve unconstrained optimization problems. Therefore, the implementation of some effective constraint-handling mechanisms is considered a crucial issue for all biological inspired optimizers [20–23]. These methods have been classified by different authors into certain categories (see, for instance, the state-of-the-art review by [21, 24, 25]): penalty functions-based methods; methods based on special operators and representations; methods based on repair algorithms; methods based on the separation between OF and constraints; hybrid methods. The most adopted method due to its simplicity is the exterior penalty approach which allows to convert the problem in an unconstrained one [26, 27]. Many different approaches such as the death, static, dynamic, or adaptive penalty functions have been proposed in time, e.g., see [26]. A proper choice of the constraint-handling mechanism affects the performance of the algorithm, and one of the critical issues to take into account is the preservation of the diversity of the population [28]. The brutal elimination of the unfeasible particles, such as in the death penalty rule, can jeopardize the exploration performances due to a loss of information [21, 29]. Moreover, it has been verified that several traditional penalty-based approaches may not be adequate to deal with highly complex search spaces, especially for problems in which several constraints are actively close to the optimum [30]. Specific constraint handling operators have been implemented in PSO, as presented in [26], and can be classified in four main groups: penalty-based mechanisms; separatist mechanisms; hybrid mechanisms; other constraint handling mechanisms (such as Del Valle’s ranking approach). Many variant of the PSO have been proposed in the years, e.g., the ILS-PSO by [31] which adopts a local search operator to deal with equality constraints or the PSO + by [27] which is based on the preservation of the feasibility. In the last decades, to overcome the drawbacks of a single heuristic approach, many hybridization have been proposed such as with the simplex method [32], with sequential quadratic programming [6], with radial basis function approximations [33], with fuzzy logic [9, 34], with neural networks [35], with parallel computing techniques [36], with genetic algorithm GA-PSO [12], with differential evolution algorithm DEPSO [37], with ant colony optimization [5], and with many other algorithms such as exposed in [14]. In the following, a new nonpenalty constraint handling approach based on the support vector machine (SVM) is reported in detail, and, after that, two numerical problems are discussed to verify the convergence properties of the new PSO-SVM algorithm. Subsequently, the proposed PSO-SVM is adopted into the structural optimization field analyzing two examples, respectively, related to the size optimization of a simply supported beam and to the size and shape optimization of a Warren truss beam.


## 2. Constraint Handling in the Proposed Framework: The PSO-SVM Approach

In this paper, a new nonpenalty-based constraint handling approach is implemented using the SVM. This machine learning approach is based on statistical learning theory [38]. During the learning process, the machine learns from examples contained in a known training dataset which could be composed of both input and output data (supervised) or only input data (unsupervised) [39]. The SVM is able to map the input data to another space (feature space) usually with a higher dimension [40] where the data are linearly separable and it will search the optimal separating hyperplane according to the principle of the maximal margin [41]. The transformation is performed by a Kernel Function which represents the inner product in the feature space (the so-called kernel trick procedure [39]). Then, after the remapping of the optimal separating hyperplane in the starting design space, the result of this process is the definition of a predictive model that could map the output of other new input data. If the output is a variable from a finite set, which represents the class (or category or label) of the input data, the problem is called classification (or pattern recognition) problem [39, 42]. After a training phase on the current positions of the swarm, the new trial positions are labelled as feasible or unfeasible using the trained predictive model. In this case, since the algorithm has to separate feasible positions from unfeasible ones, it is performing a binary classification problem. The SVM is based on the inner product in feature space, so it is possible to avoid explicitly writing the transformation of training data into the feature space, but it is possible to directly operate defining a proper kernel function, as stated in [39]. This approach is called in [40] as implicit mapping or implicit embedding. According to [43], typical widely spread kernel functions are polynomial kernels and Gaussian kernel [39, 40] which is the isotropic form of the radial basis kernels, RBF [42]. In the proposed approach, the Gaussian kernel is adopted. It can occur that training data are nonseparable, i.e., which cannot be linearly divided without committing a certain error [41]. In this case, the SVM operates with a soft margin defining a hyperplane minimizing the error. This is governed by a penalty coefficient C>0 also know as regularization parameter which controls the tolerance on training classification errors and the margin position related to the complexity of the model [39, 42, 43]. For the proposed algorithm, it was found that a good trade-off between accuracy and computational effort is obtained using C between 100 and 1000.The new proposed approach is implemented in Matlab® code, and its aim is to integrate the PSO algorithm with SVM in order to separate feasible positions from the unfeasible ones and reduce the search space.In Figure1(a), the flowchart of the proposed algorithm is shown. The initial random positions, which are obtained by the LHS technique, are labelled with respect to the constraints, and the SVM is able to construct the nonlinear boundary of the feasible region, thanks to the kernel trick. In order to get very accurate results still from the beginning, it would be necessary to sample a very huge swarm size. However, since the particles move in another configuration at each iteration, new available positions can be used to train a new SVM model at each iteration; therefore, this incremental approach is adopted in the PSO-SVM. Since the constraints are defined as gq≤0, if at least one of them is greater than zero the entire n-tuple (which is a single individual candidate solution) is labelled as unfeasible (yi=−1), otherwise it is labelled as feasible (yi=+1). Unfortunately, sometimes, the SVM really struggles to work properly, e.g., with a very narrow search space. In fact, the preservation of the feasible points only, into a wide unfeasible space, leads to instabilities such as overfitting or underfitting problems of the SVM. To improve the performance of the algorithm, a relax constraint function ψik is defined, wherein the subscript i refers to ith constraint and the k refers to the current iteration. This approach leads to enlarge the real feasible space to a fictitious wider one because of a relaxation of the actual constraints. This means that the real constraints are “moved” from their actual position to a fictitious one through a proper choice of ψik which is not trivial to generalize. This procedure acts as substituting the original inequality constraints with the following relaxed constraints:(4)gi,kx≤ψik⟶gi,kx−ψik≤0.Figure 1
Flowchart of the proposed PSO-SVM algorithm. (a) General flowchart of the proposed algorithm; (b) flowchart of “Max Position Correction” block.
(a)(b)Theψik factor is tuned on the standard deviation σu,k of the amount of all unfeasible points detected in all generations until the current one. This approach allows a better exploration during the first generations and, subsequently, it could be forced to zero implementing a decreasing function, which directly depends on the current generation k and the number of maximum generations kmax. Defining a further reduction factor λ of the standard deviation directly chosen by the user, the following relaxation functions are implemented:(i)
Constant relax:ψik=λσu,k(ii)
Piecewise constant relax:ψik=λσu,k,ifk≤kmax/20,otherwise(iii)
Linear decreasing relax:ψik=λσu,k−λσu,k/kmaxk(iv)
Piecewise linear decreasing relax:ψik=λσu,k−2λσu,k/kmaxk,ifk≤kmax/20,otherwise(v)
Parabolic decreasing relax:ψik=λσu,k−λσu,k/kmax2k2(vi)
Piecewise parabolic decreasing relax:ψik=λσu,k−4λσu,k/kmax2k2,ifk≤/20,otherwiseThe proposed algorithm implements the SVM with soft margin with regularization parameterC=100. Indeed, with respect to a hard margin, to speed up the algorithm, it is numerically convenient to adopt a high box-constraint value, e.g., between 100 and 1000, allowing a certain misclassification rate but assuring all the possible solutions lie in the feasible space boundary respecting the constraints with a certain engineering tolerance. After the training phase, the unfeasible points of the initial population have to be randomly resampled in the feasible region adopting the LHS technique, and, therefore, they are discarded and sampled again until their label becomes yi=+1. At this point, the initial velocity of the particles is randomly sampled always using LHS, and then, the evolutionary phase of the algorithm can begin. The velocity for the next generation is calculated as in equation (2), but it is necessary to check if it respects the allowable maximum velocity (“Max Velocities Correction” block in Figure 1(a)). The next block “Max Position Correction” is a fundamental step because it is the check of the feasibility of the position given by (3), and, for greater clarity, it is expanded in Figure 1(b). First of all, a temporary final position is computed with equation (3) for each particle, and the algorithm checks if it lie inside the box-search space or moves it to the nearest boundary of Ω if it lies outside. If the temporary position label is yi=+1, the new temporary position lies in the feasible region and it is stored, otherwise a simply bisection approach is performed. The velocity vector which leaded to the temporary position is firstly divided by p=2 and a new temporary final position is computed. If the label of this new temporary position becomes positive, so this new position is accepted and stored for the next generation. Instead, if it is still negative, the algorithm will redouble p, getting a new velocity vector and a new temporary position to label and so on. When all the elements of temporary velocity vector tend to zero, it is possible to leave the original position for the next generation. In this way, the particle can move only if the landing position is into the feasible region, otherwise it stands still. After that, new positions are available to increase the database set and train a new refined SVM model getting a new boundary of the feasible region. It is possible that some of the points which were inside the previous boundary may lie outside of the new boundary. In this case, it is necessary to resample them as before using LHS until their label becomes positive. Despite this, the latter naive approach may slow the algorithm, and it is beneficial for the diversity of the population, allowing new starting search positions. The update of the gbest is performed considering only the real feasible points respect to the actual constraints and not the relaxed ones. This procedure ensures that the algorithm gets good and feasible results also using a constant relax constraint function during the generations. The stopping criterion can be set to the achievement of a maximum number of iterations kmax. A further strategy to improve the behaviour of the PSO is reducing the maximum possible velocity range of the particles, governed by γ, updating this latter during generations. For example, it is possible to set γ/2 starting from kmax/3 and to the minimum value (γ=0.1) in the last generations to improve the exploitation. In the following, the PSO-SVM is firstly tested on some numerical mathematical literature problems and statements are reported in Appendix. To make some comparison of the results, in these last examples, the OF value is compared with the genetic algorithm of Matlab ® and the PSO with penalty approach provided from the code proposed by [44].


## 3. Numerical Example 1: Sickle Problem [45]

The first numerical example (from literature [45]) is the benchmark test 1 stated in the Appendix also known as sickle function. Since this optimization problem has two design variables, it is possible to follow the optimization process of the PSO-SVM in the design variables’ plane. In Figure 2(a), the objective function and the constraints are graphically represented as a three-dimensional graph. The projection on the design variables’ plane of the feasible region is really narrow because it is formed by the space between the two constraints’ parabolas. For this problem, a population size of 50 individuals is set, and in Figure 2(b), the first random generation is shown and the support vectors are emphasized. After that, the unfeasible points are resampled until all the population falls inside the SVM-based boundary. At each generation, new points are added to the SVM training data and the boundary is improved. Since the search space is really narrow, running the algorithm with no relax constraint function results in poor performance in defining the SVM boundary (black solid line) with respect to the actual one (dashed lines), as shown in Figure 2. It is possible to improve the performance in terms of the OF decreasing history adopting the relax of constraints, e.g., the constant relax function, as shown in Figure 3. In this latter example, the reduction factor of the standard deviation of the unfeasible point is chosen as λ=0.5. The PSO-SVM performances using constant relax function is compared with other algorithms in Table 1. The comparison is carried out in terms of the objective function value to the optimal solution after 50 runs with kmax=100 each run and collecting mean value μ, standard deviations σ, best OF, and worst OF. As one can check, the convergence of the new proposed method is satisfied since the objective function value is getting close to the global optimum with a little standard deviation as the other existing methods.Figure 2
Example 1 (sickle problem [45]): case no relax constraints’ function. (a) 3D graph of sickle problem design space. (b) Generation 1. (c) Generation 2. (d) Generation 50. (e) Generation 100. (f) Objective function history.
(a)(b)(c)(d)(e)(f)Figure 3
Example 1 (sickle problem [45]): case constant relax constraints’ function. (a) Generation 2. (b) Generation 50. (c) Generation 100. (d) Objective function history.
(a)(b)(c)(d)Table 1
Numerical Example 1 (sickle problem [45]): comparison of PSO-SVM, GA, and PSO-penalty.
PSO-SVMPSO-PenaltyGAμx10.84330.87210.8502σx10.00041.121e − 150.0172μx214.095214.109114.0986σx20.00028.972e − 150.0087μOF−6.9614e + 03−6.9291e + 03−6.9537e + 03σOF0.46094.59e − 1219.2968Best OF−6.9595e + 03−6.93e + 03−6.9618e + 03Worst OF−6.96e + 03−6.93e + 03−6.8547e + 03


## 4. Numerical Example 2: Five Design Variables’ Optimization Problem

The second numerical example, in which the statement is reported in Appendix, is a literature benchmark test optimization problem with five design variables and six constraints. In this example, a comparison between different histories of the objective function using different relax constraint functions is performed. As shown in Figure4, all the examples give a good result and tend to converge to the exact solution with different decreasing rates. In general, one can notice that, in piecewise functions, in the first iterations, the algorithm generally boosts the exploration instead of the exploitation, which is usually enhanced in the second half with a zero relax coefficient. This feature is important because it can affect the performance of the proposed algorithm with different kinds of problems and the user needs to try different relax functions in order to find the most suitable for the specific problem to solve. In Tables 2 and 3, a comparison with the GA and the PSO-Penalty is performed in terms of OF value and optimal design points, running the codes 50 times with kmax=100 each run and collecting mean value μ, standard deviation σ, worst, and best. Also, with this more complex optimization problem, the convergence of the new PSO-SVM method is satisfied, getting an OF value close to the global optimum with a little standard deviation as the other existing methods.Figure 4
Numerical example 2: objective value history comparison among different relax constraint functions for a single run.Table 2
Numerical example 2: comparison among PSO-SVM with different relax constraint functions (check Appendix and [46]).
Comparison among constraint relax functions in PSO-SVMDesign Var.No relaxConst. relaxPiecewise Const.Lin. relaxPiecewise lin.Parabolic relaxPiecewise Par.μx178.000478.008678.002678.000978.000678.002778.0077σx10.00130.02430.01200.00380.00410.01000.0526μx233.006233.004033.005633.010433.010733.074033.0052σx20.01580.01830.01430.04380.03220.22740.0158μx330.002729.998530.003730.000930.003830.036630.0027σx30.01010.01050.01250.02350.01780.12170.0160μx444.656944.236044.364344.985844.800244.081444.6662σx41.59862.96302.00270.06171.07343.28621.3233μx536.898137.075537.013636.767136.836037.044636.8931σx50.64441.19560.80360.06520.43101.38430.5324μOF−3.0655e + 04−3.0644e + 04−3.0647e + 04−3.0664e + 04−3.0659e + 04−3.0634e + 04−3.0655e + 04σOF43.018778.920553.77733.991828.993187.996736.2059Best OF−3.0666e + 04−3.0666e + 04−3.0666e + 04−3.0666e + 04−3.0666e + 04−3.0666e + 04−3.0666e + 04Worst OF−3.0375e + 04−3.0185e + 04−3.0449e + 04−3.0643e + 04−3.0468e + 04−3.0186e + 04−3.0452e + 04Table 3
Numerical example 2: results from PSO-SVM without relax constraints, PSO-Penalty, and GA (check Appendix and [46]).
PSO-SVMPSO-PenaltyGAμx178.00047878.0004σx10.001300.0024μx233.00623334.2398σx20.015800.7052μx330.002729.996730.801σx30.01012.15E − 140.378μx444.65694545σx41.598600μx536.898136.773634.8023σx50.64440.00e + 000.9057μOF−3.0655e + 04−3.0665e + 04−3.0531e + 04σOF43.01871.47E − 1165.8031Best OF−3.0666e + 04−3.0665e + 04−3.0660e + 04Worst OF−3.0375e + 04−3.0665e + 04−3.0378e + 04


## 5. Structural Example 1: Simply Supported Beam

In Figure5, it is considered as an ideal simply supported beam of length L with a constant cross-section A=b⋅h loaded with a distributed constant load q which is supposed to be much greater than the self-weight for the sake of simplicity. The aim is to minimize the weight of this structure respecting the tensional constraints and maximum deflection constraint due to only the q load. The self-weight is proportional to the volume V of the material, as stated by [47]. The objective function is fd=ρV=ρAL, where d is the design vector and ρ is the material density which is supposed to be constant. In this case, only stress constraints on normal stress σ, tangential stress τ, and maximum deflection vz are considered. This is a typical sizing optimization problem. The design vector d=d1,d2T contains the design variables which are changed during the optimization process, i.e., in this case, d1=b and d2=h. Since the cross-section as well as the cross-sectional dimensions must be greater than zero, it implies the presence of a new constraint to satisfy. Performing an elastic analysis, the maximum allowable stress is the yielding stress σy, and it is possible to use the Navier Formula and the Jourawsky Formula for the normal and tangential stress, respectively. The maximum moment is in the middle span, whereas the maximum shear force is in correspondence of the supports z=0 and z=L:(5)MEdz=L2=qL28,VEdz=0,L=qL2.Figure 5
Problem formulation: simply supported beam with constant cross-section.Recalling the elastic resistance modulus for a rectangular section and using the Navier formula, it is possible to write the maximum normal stress in the middle span as(6)Wel=bh26=d1d226,σz=L2=MEdWel=34qL2d1d22.Using the Jourawsky formula, it is possible to write the maximum tangential stress in the middle of cross-sectiony=0 (parabolic tangential stress diagram on a rectangular section) as(7)τz=0,z=L=VEdSx∗y=0Ixb=32VEdbh=34qLd1d2.In order to take into account both normal and tangential stresses, it is necessary to refer to a yield criterion. In this case, the Von Mises yield criterion is adopted:(8)σ2z+3τ2z≤σid.Substituting (6) and (7) into (8), respectively, it is possible to obtain the two expressions of the stress constraints in middle span z=L/2 (pure moment) and z=0,L (pure shear):(9)34qL2d1d22≤σid,34qLd1d2≤σid3.The maximum deflectionvz=L/2 can be calculated using the virtual work principle obtaining(10)vz=L2=5384qL4EIx=532qL4Ed1d23.The complete statement of the optimization problem is the following:(11)minfd1,d2=d1d2,s.t.d1>0,d2>0,34qL2d1d22−σid≤0,34qLd1d2−σid3≤0,532qL4Ed1d23−vmax≤0,where the constant ρL is dropped by the objective function as stated in similar problems analyzed in [48], σid is the ideal Von Mises normal stress, and vmax is the maximum deflection admissible by reference design codes, i.e., in this case, it is fixed to vmax=L/250. It is possible to define the fixed variables’ vector b=b1,b2,b3,b4,b5T which contains problem data which does not change during the optimization process, i.e., in this case, b1=q, b2=L, b3=σid, b4=E, and b5=vmax. Because of the amount of fixed parameters, it is more convenient to work with a dimensionless form of the same problem. Posing d˜1=b/L and d˜2=h/L, the new dimensionless objective function become f˜d˜1,d˜2=fd1,d2/L2=d˜1d˜2, whilst the normal stress constraint (9), the shear constraint (9), and the deflection constraint (10) become, respectively, as(12)34qLσid1d˜1d˜22≤1,334qLσid1d˜1d˜2≤1,5⋅12516qEL1d˜1d˜23≤1.It is useful to define two dimensionless nonnegative parameters, collected inb˜, which completely characterize the fixed variables of the problem:(13)ψσ=qLσid,ψE=qEL.Finally, the dimensionless version of problem (11) is given by(14)minf˜d˜1,d˜2=d˜1d˜2,s.t.d˜1>0,d˜2>0,34ψσ1d˜1d˜22≤1,334ψσ1d˜1d˜2≤1,5⋅12516ψE1d˜1d˜23≤1.To solve this problem, the PSO-SVM is adopted with the piecewise linear decreasing relax constraint function with a user coefficientλ for standard deviation of the unfeasible points fixed to λ=0.05. For academic purposes, in order to graphically analyze the behavior of the constraint handling, the dimensionless parameters are fixed as ψσ=ψE=0.2, looking for an optimal solution in the design domain for the dimensionless design variables as 0≤d˜1≤1.5 and 0≤d˜2≤1. In this way, with this particular choice of the dimensionless parameters, the constraints intersect each other creating discontinuous nonlinear boundary of the feasible region. The population size is always 50 individuals and the maximum iterations are 100. After 50 runs, the results show a quite great variability of the design variables but always at the same objective function value. This fact enlightens the presence of a front of possible optimal solutions. In fact, as shown in the graphical representations in Figure 6, for this specific choice of b˜, only two constraints are active, and there exists a region of optimal solutions on the τ constraint boundary line. The comparison in Table 4 among PSO-SVM, PSO-Penalty, and GA shows the mean value and the standard deviation of the best objective function value obtained after 50 runs.Figure 6
Structural example 1: simply supported beam, case piecewise linear decreasing relax function. (a) 3D graph of simply supported beam problem design space. (b) Generation 1. (c) Generation 2. (d) Generation 50. (e) Generation 100. (f) Objective function history.
(a)(b)(c)(d)(e)(f)Table 4
Structural Example 1: results from PSO-SVM piecewise linear decreasing relax constraints, PSO-Penalty, and GA.
PSO-SVMPSO-PenaltyGAμOF0.25980760.25980760.2610035σOF2.04e-095.61e-172.58e-03Best OF0.25980760.2598080.2598076Worst OF0.25980760.2598080.2664507In order to find the entire front of the all optimal possible solutions for this specific problem, it is necessary to find all the pairsd˜1,opt,d˜2,opt posing the objective function as fd˜1,d˜2=d˜1⋅d˜2≈0.2598. As one can see in Figure 6, this optimal front corresponds to a part of the τ constraint (12) posed as an equality. Referring to the d˜1 optimal possible values, the optimum front is upper bounded from the σ constraint (12) posed as an equality and lower bounded from the box search space limits. Then, to calculate the optimal upper bound d˜1,opt,UB, it is necessary to calculate the intersection between the two aforementioned constraints:(15)34ψσ1d˜1d˜22=1,334ψσ1d˜1d˜2=1,⟶34ψσ1d˜1=334ψσ1d˜1⟶d˜1,opt,UB=94ψσ=0.45.As before, to calculate the optimal lower boundd˜1,opt,LB, it is necessary to calculate the intersection between the equality τ constraint and the horizontal line d˜2=1, obtaining(16)334ψσ1d˜1d˜2=1,d˜2=1,⟶d˜1,opt,LB=334ψσ=0.2598.Finally, considering (15) and (16), it is possible to obtain all the optimal pairs d˜1,opt,d˜2,opt using the following equation:(17)d˜2,opt=0.2598d˜1,opt,with0.2598≤d˜1,opt≤0.45.Since the algorithm works with dimensionless parameters, in order to find the physical dimensions of the optimized cross-section, it is sufficient to multiply the obtained valuesd˜1,opt,d˜2,opt by L.To show a technical possible application coming from this simple example, only a concrete beam with span lengthL=3 m is now considered. Disregarding for the moment the self-weight load, the q load set to 15 kN/m represents only a live load. The concrete modulus is set to E=25 GPa and the Von Mises ideal stress is related to the tensile stress of concrete set to σid=3 MPa. Considering the box search space as 0≤b≤40 cm and 0≤h≤45 cm, the algorithm found the minimum weight respecting the constraints with b=16.67 cm and h=42 cm. Rounding-off these values, the self-weight associated to a concrete beam with b=18 cm and and h=45 cm is equal to(18)G=γconcrete⋅b⋅h=24⋅0.18⋅0.45=1.944kN/m.For sake of simplicity, addingG to q a new load equal to 16.944 kN/m which takes into account also the self-weight is defined. Launching again the algorithm, new optimal exact dimensions are now obtained: b=18.83 cm and h=45 cm. Rounding-up the exact solution, a new self-weight equal to G=2.16 kN/m is coming from a section with b=20 cm and h=45 cm. Now, the convergence is reached because the new optimal exact solution is b=19.07 cm and h=45 cm. Finally, the optimal cross-section for this concrete beam which minimizes the self-weight is given by b=20 cm and h=45 cm.


## 6. Structural Example 2: Optimization of a Warren Truss Beam

The second structural example comes from [49]. In that work, the weight optimization of an in-plane Warren truss simply supported beam, depicted in Figure 7, is performed with Differential Evolutionary Algorithm (DEA). The steel profile used for truss members is a square hollow core section, as shown in Figure 7. This kind of profile ensures good stability against buckling, and it represents a good solution for this type of structure because of its high strength-to-weight ratio [49]. On the contrary, joint connections are usually welded, and so, in order to reduce the total cost, it is important to reduce the size of sections to be welded. Regarding to the size optimization problem, as shown in Figure 7, this kind of sections are completely described by only two independent design variables: the outer dimension of the cross-section B and the thickness of the webs s. In this problem, the truss has m members belonging to four different types of cross-sections, as reported in Figure 7: lower chord (B1,s1), upper chord (B2,s2), internal webs (B3,s3), and external webs (B4,s4). To perform the shape optimization further, two design variables are considered: vertical height of the external webs Hmin and the maximum height Hmax. The design vector is therefore defined as(19)x=B1,B2,B3,B4,s1,s2,s3,s4,Hmin,Hmax,and the box search space domain Ω is defined by 60≤Bi≤360mm,4≤si≤30mm, and 50≤H≤L mm, where i=1,2,3,and4 and L is the total span length. Considering the maximum value of thickness smax=30mm, the minimum value of dimension B cannot be assumed less than 2smax due to geometric limits. The objective function is represented by the total weight of the structure [49]:(20)Wx=∑i=1mρiliAi,where ρi=7.85t/m3 is the steel density supposed equal for all members, li is the length, and Ai is the cross-section of the ith member. The structural steel used in this example is a S275, and the modulus of elasticity of the steel is 210 GPa. Regarding topology optimization, in [49], for fixed length L=20m, the optimal number of bays which divides the lower chord is 20. The external load is as a uniformly distributed load q=100kN/m applied on the lower chord acting as point loads in the nodes of the truss. The constraints are represented by the strength verifications about tensile stress (without any holes), compression stress, and buckling instability according to Eurocode 3 (EN 1993-1 2005 and EN 1993-2 2006). Despite the Eurocode γM0=1 and γM1=1.1 for bridges being recommended, to be more safe, the partial safety factors are set both equal to γM0=γM1=1.1. Another constraint to satisfy is the maximum deflection which is usually set to ulim=L/500 for bridges like that. For the current case, model uncertainty was not considered; however, the eventual adoption of specific model uncertainty coefficients (e.g., Castaldo et al., [50]) can be easily included in resistance or deformability checks. In order to make a comparison with the results of [49], for academic reasons, there is no distinction of the load combination for the strength verifications and for the deformability checks. The verification equations do not need deeper examination because this is beyond the scope of the present document. Therefore, the optimization problem statement is the following [49]: Findx∈Ωsuch that(21)minfx=Wx,s.t.NEdNt,Rd≤1,NEdNc,Rd≤1,NEdNb,Rd≤1,umax≤ulim.Figure 7
On the left, problem formulation of simply supported truss Warren beam. On the right, square hollow core tubular section with the indication of the design variables.In order to solve (21), the FEM structural analysis was performed in the Matlab® CALFEM and PSO-SVM was adopted for the optimization process. In the PSO-SVM, a population size of 100 individuals is set with kmax=200 iterations and a constant relax function with user parameter λ=3 is applied to standard deviation of unfeasible points. It is performed 50 time runs, and the best-obtained solutions are collected in Figure 8. As one can see in Figure 8, due to the complexity of the problem, sometimes the algorithm does not reach the optimum and stack in a local minimum. As shown in the graph, the optimum solution is around 3.1 t; then, it is possible to cut the graph considering only the 21 runs over the total 50 which are characterized by a best OF solution lower than 3.1 t (dashed line). In this way, the possible outliers are excluded and now it is possible to perform the postprocessing searching for the real best solution.Figure 8
Structural example 2: Warren truss. Results from 50 times run PSO-SVM.Considering the abovementioned solutions, the obtained results shown in Table5 have a quite large standard deviation in terms of design variables but very low in terms of objective function. This means that it is possible to find many combinations of design variables which are giving always the almost same objective function. The best design value of the 50 runs, also reported in Table 5, is taken into consideration. It is possible to compare this latter objective function (3.074t) with the optimal exact solution given by the original DEA code output [49]. The DEA optimal solution was characterized by the weight of 2.95t, so this is in the right order of magnitude. The comparison of the design variables is based only on general observation in accordance with the literature, as affirmed in [49]. In fact, for instance, it is expected that, mainly due to instability problems, upper chord and external diagonals would be bigger than the lower chord and internal diagonals. Finally, the best exact solution may be trivially rounded-up to get an industrialized more realistic design. The new design variables are reported in 5, and as one can see, this solution is more conservative, and it leads to an increased total weight (3.189t).Table 5
Structural example 2: mean valuesμ and standard deviations σ of best results from 21 solution over 50 runs of PSO-SVM with OF less than 3.1 t. Last three columns: best exact solution, trivial rounded-up solution, and refined industrial solution.
(mm)μ exact sol.σ exact sol.Best sol.Trivial industrial sol.Best industrial sol.B172.515.794.29595B2190.270.2128.1130105B3128.91.4128.8130130B4211.7107.9132.1135110s16.01.5444s214.77.418.82026s34.00.03444s414.08.414.71519Hmin399.421.5410.5410410Hmax4064.7113.54145.041454145OF (t)3.08980.00713.0743.1893.092If one wants to find a more accurate industrial solution, it is necessary to perform a more accurate analysis of the obtained results. Since the topology optimization was already taken into account in [49], one has to remember that, in the design variables, the algorithm is performing the size optimization and the shape optimization. This latter is regarding to the definition of the Hmin and Hmax values. Considering only the abovementioned 21 solutions and the standard deviations of Hmin and Hmax, it is possible to assume that the rounded-up values of the best exact solution can represent a good result for shape optimal parameters values: Hmin≈410mm and Hmax≈4145mm.Once the shape optimization is solved, regarding the size optimization, firstly the best exact solution as the optimal one, one has to remember thatBi and si were chosen as design parameters because of their independence, but in the optimization process, they are connected. In fact, in both objective function evaluation and constraints evaluation, these two parameters are combined into the resisting cross-section value. It is possible to obtain almost the same value of cross-section with different combinations of the design parameters. In particular, one can refer to the optimal exact solution in terms of resisting cross-sections which represent the best solution in terms of both strength verification and minimization of the weight. As one can check, for the best exact solution, the section class of all members is 1, but other optimal solutions within the 21 considered are characterized by class 4 profile. In this case, to get the strength verification satisfied, it is necessary to refer to the resisting effective area. Usually in the design, when it is possible, it is preferred to avoid class 4 profiles, and the best condition is to find an optimal solution with class 1 profile. Therefore, the best industrial solution which respects all the constraint is given by all the pairs Bi,si with i=1,2,3,and4 which gives class 1 profiles and the minimum value of area greater or equal to the effective areas requested by the best solution in Table 6. This procedure allows us to find the best solution which respects the strength verification only. For the instability verification, it is necessary to take into account also the second moment of inertia which conditions Euler’s critical load Ncr and consequently the dimensionless slenderness λ¯ which influences the reduction factor χ. Fixing the thickness si to discrete values (rounded with 1 mm of precision), starting from the best-found solution and making an iterative discrete research to find Bi (rounded with 5 mm of precision) respectful of our abovementioned design rule, the best optimal industrial solution is found and reported in the last column of Table 5. As one can check, the minimum cross area of internal webs could be assured by the pair B3,s3=105;5mm, but, due to instability problems, it is necessary to take into account the inertia and choose a profile that ensures both strength and instability verifications. The best industrial structure is verified by all strength and instability constraints. Making a comparison between the two last columns of Table 5, the trivial rounded-up solution represents a good optimal solution in terms of objective function. In fact, minimizing the weight is important, but the total cost is also affected by other aspects, e.g., welding and detailing and labour cost. So, the solution obtained by trivial rounded-up the exact one; it can be considered as an acceptable optimal result. In Figure 9(a), the undeformed and the deformed shape are depicted. It is possible to appreciate that the node 7 and, due to symmetry, the node 15 are the nodes that undergo the most deflection umax=39.8mm; however, it is respectfully of the service limit L/500=40 mm. In order to assess the validity of the results of the optimization process, the Warren truss beam is modelled with FEM professional software MIDAS Gen®. The assessment is made not only in terms of axial force for each member, but in particular, in terms of the performance ratio. This latter represents an efficiency percentage of the usage of the steel and is given by the strength ratio between the demand and the capacity. The simply supported Warren truss is modelled through truss elements in order to guarantee pure axial behaviour of each member. The section properties of the trivial best industrial results from Table 5 are assigned to each element. The adopted steel is always a S275 and the uniformly distributed load q=100 kN/m acting on the lower chord is reconducted as a hanged load directly applied at lower chord nodes acting as concentrated forces, as depicted in Figure 9(b). As already remarked, in order to get results which are directly comparable with the Matlab code and with DEA code from [49] and for the sake of simplicity, no load combination is considered. The aim is to demonstrate that the proposed algorithm provides comparable results with DEA code which is used as a benchmark and not to make a perfect design completely respectful of the current codes. In order to take into account this latter issue, it is sufficient to consider the correct the load combination at ULS for the strength verification and SLS for the maximum allowable displacement.Table 6
Structural example 2: comparison between best exact solution and best industrial solution and relative requested and provided areas.
(mm,mm2)B;sexactRequested areaB;sindustrialbestProvided areaLower chord(94.2; 4)1442.46(95; 4)1456Upper chord(128.1; 18.8)8206.18(105; 26)8216Internal webs(128.8; 4)1996.81(130; 4)2016External webs(132.1; 14.7)6907.21(110; 19)6916Figure 9
Structural example 2. (a) Optimal Warren truss, solid Line: undeformed shape; dashed line: deformed shape. (b) View of the model of the Warren truss beam on Midas Gen®. (c) Planar view of the Warren truss beam on Midas Gen® with the axial force values.
(a)(b)(c)As reported in Figure9(c), the results obtained from the Matlab code are equal to the FEM software. The overall behaviour shows that the upper chord is entirely compressed, whereas the lower chord is entirely tensed. The internal webs are alternatively compressed and tensed as usual for truss beam of this typology. Calculating the strength ratio, it is worth noting that a performance ratio between 75% and 98% is obtained for the members. These remarkable results show the importance of the optimization process during the design phase which can strongly support the decision process of the designer.


## 7. Conclusions

After an introductory part on PSO and the state of the art of constraint handling techniques, the paper presented a new valid alternative not-penalty method to solve constrained optimization problems. The main advantage with respect to the most-used nowadays penalty approach is represented by the generality of the machine learning SVM algorithm. Since it depends intrinsically on the inner product of the data, it is more adaptive even with discontinuous and nonlinear boundary of the feasible region in the design space. In order to improve the behavior of the proposed algorithm to deal with very sharp and narrow feasible regions, a relax constraint function is also implemented. From the computational point of view, the trade-off has been found adopting a population size not too large and using an incremental boundary update. It would be also possible sampling a huge initial random data and leaving the boundary fixed during the generations; however, this procedure does not lead to a good result in terms of objective function. Finally, the two numerical benchmark examples demonstrated the convergence of the new method in comparison with another penalty approach and with a GA. The last two examples highlighted the adaptability of this new method even into the structural optimization field. In particular, in the Warren truss beam problem, the optimization algorithm provided a numerical exact solution which can be easily industrialized by the designer with a trivial rounding-off without jeopardizing the optimization process. Although the Warren truss beam example is performed under simplified assumption in order to make comparisons with DEA code from [49], from the technical point of view, the new optimization algorithm becomes a really useful and powerful support for the designer during the design and the decision process. It is important to stress that working with metaheuristic algorithms always involves the definition of the value of many arbitrary parameters. Starting from literature suggestions for these values, it is always strongly suggested to perform a fine-tuning of some of these parameters also with a trial-and-error approach for each specific problem in order to find the best optimal results in terms of objective function convergence, computational effort, and elaboration time.Future applications of the proposed PSO-SVM framework in structural engineering are different. Among the different possibilities, an effective exploitation could address reliability-based calibration of partial safety factors or semiempirical coefficients included in design expressions (e.g., Mancini et al. [51]). Other applications could include seismic reliability or expected annual loss optimization for structural problems which could benefit (e.g., Di Trapani et al. [52] and Cavaleri et al. [53, 54]).



---

*Source: 6617750-2021-02-16.xml*