Inverse Distance Weighting model An improved version of Inverse Distance Weighting metamodel assisted Harmony Search algorithm for truss design optimization 1

This paper focuses on a metamodel-based design optimization algorithm. The intention is to improve its computational cost and convergence rate. Metamodel-based optimization method introduced here, provides the necessary means to reduce the computational cost and convergence rate of the optimization through a surrogate. This algorithm is a combination of a high quality approximation technique called Inverse Distance Weighting and a meta-heuristic algorithm called Harmony Search. The outcome is then polished by a semi-tabu search algorithm. This algorithm adopts a filtering system and determines solution vectors where exact simulation should be applied. The performance of the algorithm is evaluated by standard truss design problems and there has been a significant decrease in the computational effort and improvement of convergence rate.

An improved version of Inverse Distance Weighting metamodel assisted Harmony Search algorithm for truss design optimization 1

INTRODUCTION
Employing exact simulations are very common in optimizing engineering design.The problem with some of the most common and robust optimization algorithms such as Genetic Algorithm, Ant Colony, and Harmony Search, is that they entail a large number of iterations for reaching the optimum solution [1][2][3][4].This is a significant barrier when applying exact simulations to real life engineering optimization problems.It has long been recognized that approximations or metamodeling techniques are most effective tools for reducing the computational effort in complex problems.The basic approach is to replace the computationally expensive simulation with a compatible approximate, which is then used in optimization runs.This inexact model is often referred to as metamodel (model of a model).A survey on the use of metamodels in structural optimization has been carried out by Barthelemy and Haftka [5] and recently, Jin et al. [6] and Martin and Simpson [7] have reviewed metamodeling techniques for some engineering optimizations.

Observation point and optimization variable
The new value of

Set of possible range of each decision variable
The maximum value of kth decision variable The minimum value of kth decision variable random variable to be estimated The value at the observation point Optimization using metamodels entails [8]: 1.The exact model is replaced by a low order polynomial.This is to reduce the number of exact simulations and smooth out the numerical noise.
2. Enables separation of the analysis code from the optimization routine and eases the integration of codes from various disciplines.
Latin American Journal of Solids and Structures 10(2013) 283-300 3. Provides an overall view of the design search space.4. Provides an in depth knowledge of the problem domain and hence makes it easier to adjust important design parameters during the optimization process.The shortcoming of metamodel, especially those using polynomials for their approximations, is that the cost of providing precise data for fitting the global approximations increases rapidly with the number of design variables, which makes it difficult to construct a good global approximation with low order polynomials.Some well-known metamodel methods are: Inverse Distance Weighting (IDW) [9]; Polynomial Regression (PR) [6,8]; Moving Least Squares Method (MLSM) [8]; Kriging (KG) [7,[10][11][12][13][14]; Multivariate Adaptive Regression Splines (MARS) [15]; and Radial Basis Function (RBF) [6,8].
The algorithm presented in this paper is a combination of IDW and HS, which provides a more robust and global search capabilities.An immature version of IDW+HS algorithm has been already published by the authors [16].In this research, the prior algorithm is enhanced by normalizing the results and then polishing the final outcome using a semi-tabu search method.IDW is selected for its applicability and ease of use for multidimensional metamodeling problems and HS is selected for its ability to support continuous variable functions [3,[17][18][19][20][21].The numerical results shows that the enhanced IDW+HS algorithm (IDW+HS+Tabu), in comparison to the pure HS algorithm as well as the other conventional meta-heuristic algorithms (GA and ACO), leads to a lower computational cost and higher convergence rate.

APPROXIMATION ALGORITHM SELECTION
As mentioned above, the IDW model is employed in this research for approximation.The easier extendibility and competitive computational cost of this model compared to the other interpolation methods such as Kriging (KG), Polynomial Regression (PR), Multivariate Adaptive Regression Splines (MARS), Radial Basis Functions (RBF) caused the IDW to be chosen and employed in this work [6,22].The dominant problems of these methods are their high computational cost and dimension limit.
One serious problem of the KG, for example, is the large number of observation points needed to estimate the result [23].On the other hand, as in the IDW, the kriged estimate is a weighted average of the values at the observation points in which the weights are obtained by solving (1).
(1) where is the fitted semi-variance of and , y is a Lagrange multiplier and is the fitted semi-variance between and , is the distance from test point to observation point i x , and z is the random variable to be estimated [22].This weighting formulation signifi- cantly increases the computational effort, while having a small difference in results compared to IDW [24].The other models involve similar drawbacks.For instance, in spite of the quick convergence [25] and easy implementation of PR [6], some instabilities may arise when applying this technique to highly nonlinear problems with high-order polynomials [6,26] or it may be too difficult to take sufficient sample data to estimate all of the coefficients in the polynomial equation, particularly in large dimen-sions.MARS, also did not show a satisfying performance when a scarce set of samples is applied [6].As far as RBF is concerned, although it has performed well in metamodeling problems, when the sample points grow up, the performance of the method decreases [8].Although a Shepard like modification [9] seems to be applicable for enhancing the performance, it is not considered in this research because of its complexity and high computational effort.On the other hand, the overall accuracy essentially depends on the selected basis function for a given set of data samples in the modified versions [27].
In terms of cross-validation, the accuracy of the approximation is not highly important at the early stages.Instead, computational cost must be of a great importance in selecting the approximation method.On the other hand, while the algorithm closes to the final stages, the predicted values are becoming closer to the global optimum (Figure 1).The IDW is a multivariable model based on an interpolation, which is well suited to irregularly spaced data.In two-dimensional space, there are two general types of exact interpolation methods: Single global function and a collection of simple and local functions.The former is usually accompanied by an unmanageable complexity, while the latter is well defined and matches appropriately at its boundaries and is continuously differentiable even at the local junctions [9].
The value at any point in a plane is a weighted average of values at the data point in which weighting is a function of distance to those points.Shepard [9] has introduced (2) for interpolating values at .The function indicates that the points further away from will have lower weights. ( where is the value at the data point and is the Cartesian distance between and ( ).As approaches to a data point , approaches to zero.For , both left and right partial derivatives exist and for , there is no derivative.Empirical data shows that the higher exponents ( ) tend to make the surface relatively flat near all data points with very steep gradi- Latin American Journal of Solids and Structures 10(2013) 283-300 ent over a small interval between data points.Lower exponents produce a surface relatively flat with short blips.An exponent of not only gives seemingly satisfactory results for the general surface mapping and description purposes, but also presents the simplest form of calculation as Shepard [9] recommended.Figure 2 proves that leads to the best concurrence with the exact simulation.In Cartesian coordinates the distance between two points can be expressed as (3).).These extra modifications not only increase the computational load, but also are not easily applicable to the multidimensional space hence they are not considered in this research.

HARMONY SEARCH ALGORITHM
Harmony Search (HS) algorithm is a replication of a musical performance process.A musician's search to find a better state of musical harmony (a perfect state) [3] is similar to optimization process that seeks to find a global solution (a perfect state) as determined by an objective function.The pitch of each musical instrument determines the aesthetic quality; similarly, the set of values assigned to each decision variable determines the value of the objective function.The optimization procedure of HS algorithm consists of the following steps: (5) where = an objective function; = optimization variable; and = the set of possible range of each decision variable.The parameters of HS algorithm, that is, HMS, HMCR, and PAR are also specified at this step.HM is a memory space where all the solution vectors (sets of decision variables) are stored which is similar to the genetic pool in GA.This will be shared with IDW as its initial points.
2. Initializing HM: HM matrix shown in ( 6) is filled with as many solution vectors as the value of HMS.
These solutions are randomly generated and sorted according to the values of the objective function, .
3. Improvise a new harmony from HM: Memory considerations and pitch adjustments determine if the new harmony vector should be generated from HM.That is, (the value of variable for the new vector) is chosen from the values in HM or .The value of HMCR which varies between 0 and 1 will determine where to choose the possible new value as indicated by ( 7). ( Every component of the new harmony vector, is examined to determine whether it should be pitch-adjusted.This procedure uses the PAR parameter that sets the rate of adjustment for the pitch chosen from the HM as shown in (8). ( The HMCR and PAR parameters, introduced in Harmony Search, help the algorithm find globally and locally improved solutions, correspondingly. 4. Update the HM: The new harmony vector replaces the worst harmony if a better objective function value is obtained.HM is then resorted. 5. Check if the termination criterion is satisfied: When the termination criterion is not satisfied steps 3 and 4 are repeated.
For further information about the HS algorithm see Refs [3,[17][18][19][20][21].In this study, a combination of IDW and HS is employed to decrease the computational effort and improve the convergence rate.The presented method has been applied to a number of truss design problems, which are presented in the next section.

THE NEW MULTILEVEL OPTIMIZATION ALGORITHM
An optimal truss design is one in which the optimal cross-sectional area assigned to each member satisfies the given constraints.Design constraints typically consist of the maximum allowable compressive and/or tensile stress in any member of the truss and the maximum allowable deflection of any node [1].The truss optimization objective function can be expressed as (9).(9) where = weight of the truss; = number of members making up the truss; = material density of member ; = length of member ; and = cross-sectional area of member i, chosen from a set of areas between and , where = lower bound and = upper bound.The parameters and represent the stress and deflection of the member of the truss respectively.The algorithm consists of the following steps: 1. Generating the sample points and initializing HM: sample points are generated by selecting crosssections from an eligible set of cross-sections.At least one cross-section should be chosen near the upper bound and one near the lower bound (e.g. and respectively).The remaining selections should be distributed uniformly between these bounds.Then a sample grid matrix ( ) is generated from the chosen cross-sections, where is the number of truss members and should be greater than or equal to 2 (the minimum number of sample point entries).In other words, the entries of sample grid matrix are chosen from those sample points.For instance, (10) shows the formation of SG matrix for 25-bar truss from the sample points 0.5 and 3.0. ( HMS determines the number of cross-sections from SG matrix that HM should be filled with.An exact simulation is applied to each point in HM.Both IDW and HS get their initial points from the same HM. Note: In (9) the cross-sectional constraint is zero cost, while and impose computational efforts.
In the case where and are not satisfied, HM would not be updated.In this case, the algorithm recovers the unacceptable solution and by applying superposition principle, makes it more acceptable.That is, the cross-sectional area of each member is multiplied by (11) and then it would be added to HM. Applying such a factor to the cross-sectional area makes the stress and deflection ratios to be as close as possible to 1.The above operation will satisfy the first two constraints mentioned in ( 9) for all solution vectors, but the cross-sectional constraint may be violated.As HM is filled by order of priority from the lightest to the heaviest truss, which generally implies that HM will be filled from the smallest to the largest cross-sections, the repeated application of algorithm through its convergence process will force the satisfaction of cross-sectional area constraint.For a detailed explanation of the normalization process, assume the algorithm found the solution [0.236, 0.409, 2.094, 1.992, 0.426, 0.459, 1.106, 2.026] which is not acceptable according to the stress and deflection ratios (1.5 and 0.216 respectively).However, the solution could be recovered through the normalization process by multiplying the solution to the normalization factor, 1.5, which is obtained using (11).The recovered solution would be [0.354,0.614, 3.141, 2.988, 0.639, 0.688, 1.658, 3.039], with stress and deflection ratios of 1.000 and 0.144, which satisfies the criteria and could be added to HM.It should also be noted that the only parameter, which affects the performance of the algorithm, is the initial sample points. (11) 2. Generating a new solution vector: A new solution vector is generated according to (12) where the cross-sectional area of the member is determined using HMCR and PAR (see ( 7) and ( 8) respectively). ( where = new cross-sectional area of member; = total number of iterations; and = current step.
3. Estimating the approximate response: IDW is used to estimate the new solution response.If the new approximate response is better than the worst one stored in HM, the algorithm will proceed to the next step.Otherwise, it will return to the previous step. 4. Evaluating the exact response: An exact response evaluation process is applied to the new generated solution.If the accurate response is acceptable, the worst solution in HM is replaced.Otherwise, the algorithm will return to step 3.In further iterations, as the responses in HM get closer to the optimal response, the distance between the approximate and accurate responses is reduced.Numerical examples show that the most of unnecessary and expensive calculations are eliminated in the early iterations where the responses are far enough from the optimal response.
5. Checking the termination criteria: The algorithm terminates in a certain number of iterations ( i.e. ).Although the most of the meta-heuristic algorithms have some other criteria together with maximum iterations, we avoided to adopt such criteria here, because the IDW reduces the number of HS (exact analysis) iterations.On the other hand, the criteria excluding maximum itera- tions may lead to a premature solution.If the termination criterion is not met, the algorithm repeats steps 3 and 4.
6. Polishing the result: A semi-tabu algorithm is applied to the final solution to obtain a better solution.
The semi-tabu is to find a smaller cross section area for each member to reach a better solution, while satisfying the criteria.This is achieved by reducing each area by multiplying by 0.95 in a loop, while satisfying the criteria.It should be noted that only one member is reduced in each loop.
Figure 3 shows the optimization procedure of the IDW+HS algorithm.The efficiency of the algorithm was evaluated through some numerical examples.To investigate the computational cost and convergence rate, three examples (10, 25, and 72-bar truss) are selected from refs [1][2][3][4].The algorithm was developed using Python programming language and ran on a Centrino, 1.4 GHZ computer.The structural analysis also was done using OpenSees code.The 10-bar truss is discussed in detail to show the efficacy of the algorithm, and the other examples are given for further investigations.These problems have been widely used as benchmarks to test and verify the efficiency of many different engineering optimization methods.

Ten-bar truss
Latin American Journal of Solids and Structures 10(2013) 283-300 tion Figure 4 shows the geometry and support conditions for two-dimensional cantilever 10-bar truss with a single loading condition. .The upper and lower limits of cross-section areas of each truss member .The sample points by which the sample grid is generated are 1.0 and 30.0. Figure 5 compares the convergence rate of the proposed algorithm with that of some other algorithms (GA, ACO, HS, IDW+HS, etc.).The IDW+HS has been already compared with the other pure meta-heuristic algorithms by the authors [16].Here the IDW+HS is compared with the new modified algorithm.The comparison indicates that although the IDW+HS+Tabu fits the IDW+HS at the early stages, in the final steps, the semi-tabu algorithm leads to a lighter design.The exact design weight is shown in Table 1.It can be seen that the new algorithm designed a lighter structure compared to the IDW+HS algorithm (4893 vs 4963), in a fewer iterations (682 vs 1870).; sample points = 0.5 and 3.0.This paper describes a surrogate-based optimization algorithm.It is a combination of a metamodel and a meta-heuristic algorithm.Most of the meta-heuristic algorithms (e.g.GA, HS, and ACO) employ a large initial population size, which leads to a large and costly number of function evaluation and suffer from slow convergence.These prohibitive factors are more pronounced when applied to exact simulation.IDW+HS+Tabu algorithm implements the optimization process in two steps.First, the IDW is applied to eliminate the unnecessary calculations.Then the remaining solutions are processed for exact simulation.Then superposition principle is applied to normalize the result of exact simulation.Finally HM is updated by the results.The algorithm was applied to four truss design problems.
The obtained results showed that the proposed method not only reduces the computational effort tion but also significantly improves convergence rate.As shown above, the proposed algorithm saves the computational cost from 86 to 99 percent.
Notations Normalization factor Fitted semi-variance of and fitted semi-variance between ( ) Search domain of area of each member Lagrange multiplier The minimum allowable area of truss member The maximum allowable area of truss member Cross sectional area of each truss member Cross sectional area of each truss member generated at new steps Normalized cross-sectional area of each member The number of data points included in the circle with radius r harmony memory Current iteration number of the algorithm Total number of iterations Length of each truss member The number of truss members The number of decision variables Total number of data points Test point PAR Pitch adjusting rate The minimum number if sample point entries Radius of the circle containing the data points Sample grid matrix Structure weight Test point

Figure 1
Figure 1 Closeness of the IDW result to the global optimum 3 INVERSE DISTANCE WEIGHTING MODEL (IDW)

Figure 2
Figure 2 Influence of on Chi-Square test function Latin American Journal of Solids and Structures 10(2013) 283-300 tion 1. Initializing the optimization problem and algorithm parameters: The general form of optimization problem is specified as follows: of Solids and Structures 10(2013) 283-300 Latin American Journal of Solids and Structures 10(2013) 283-300

Figure 3
Figure 3 Optimization procedure of IDW+HS algorithm

Figure 4
Figure 4 Configuration of 10-bar truss

Figure 5
Figure 5 The convergence diagram for minimum weight design of 10-bar truss with various algorithms

Figure 8
Figure 8 The convergence diagram for minimum weight design of 25-bar truss under multiple load case with various algorithms

FFigure 10
Figure 10 The convergence diagram for minimum weight design of 72-bar truss with various algorithms

Table 1
The best solutions obtained from various algorithms for 10-bar truss

Table 2
Single load case for 25-bar truss

Table 3
Multiple load case for 25-bar truss

Table 4
The best solutions obtained from various algorithms for 25-bar truss under single load case Latin American Journal of Solids and Structures 10(2013) 283-300

Table 5
The best solutions obtained from various algorithms for 25-bar truss under multiple load case

Table 8
The best solutions obtained from various algorithms for 72-bar truss