Testing heuristic optimisation methods for vibration-based detection of damage

Considerable efforts have been made in recent years for the development of nondestructive techniques for detection of structural damage based on changes of dynamic parameters of structures. The tests are made on a model of a simple beam with cantilevers. Two heuristic optimisation techniques are used in the proposed procedure for detecting the location and level of damage: taboo search and simulated annealing. The results show that both proposes procedures, with appropriate target functions and weight factors, are highly efficient for detecting the location and level of damage.


Introduction
The presence of damage in the structure may reduce its performance by decreasing its service life, or the damage may progress and eventually cause catastrophic effects [1].Nondestructive testing techniques are generally used to investigate critical changes in structural parameters, thus preventing any unexpected failure.These techniques are included in the so-called Structural Health Monitoring -(SHM) system, which is the art of instructing information about the existence, location, and extent of damage in a structure using non-destructive methods, [2][3][4][5].One possible approach is to monitor and interpret changes in structural dynamic properties extracted from measurements during experimental modal analyses, or ambient vibration measurements and appropriate signal-processing techniques [6].The occurrence of damage in a structure produces changes in its natural frequencies and mode shapes.However, it is quite difficult to locate the damage just by using natural frequencies, because modal frequencies are a global property of the structure, and cannot provide spatial information on the change in structural modal properties.Therefore, the mode shape information is required to identify location of damage.Different methods have been used for damage detection in literature, such as measurement of changes in natural frequencies, evaluation of changes in the structure flexibility matrix, and definition of change in the dynamically measured stiffness matrix and flexural stiffness of a structure [7,8].Changes in modal properties, structural natural frequencies, and mode shapes, are compared with the same modal properties extracted through experimental analysis due to reduction in flexural stiffness.Most damage detection methods include an optimization procedure for exploring the existence and severity of structural damage, which is an essential part of the entire process.The search (design) space may be so large that the global optimum cannot be found within a reasonable time.The existing direct linear or nonlinear optimization methods may not be efficient or could be computationally expensive for solving such problems [9].Heuristic optimization that can be used to find global or near-global minimum could be suitable for such complex and difficult optimization problems.Several heuristic optimization methods are used: Genetic algorithms (GA); Particle Swarm Optimization (PSO); Artificial Neural Network (ANN); Tabu Search (TS); and Simulated Annealing (SA).The TS method is a meta-heuristic iterative procedure that starts with some initial feasible solution and attempts to determine a better solution.The TS procedure initially starts from several neighbourhood hyper-points (position), and selects the new point (move), producing the best solution among all candidate points (moves) for the current iteration.The process of selecting the best move (which may or may not improve the best solution) is based on the assumption that good moves are more likely to reach an optimal or a near-optimum solution [10].The SA is mimicry of the metal annealing principle from materials science.It relies on the fact that certain material alloys have multiple stable states, with different molecular arrangements and energy levels.During the annealing process, an alloy is initially heated up to a certain temperature, with all particles randomly distributed in the liquid phase.Afterwards, the temperature slowly decreases, until the material solidifies.If the temperature decreases sufficiently slowly, the annealing process always guarantees that the alloy will reach its global minimum energy state [4,11].The TS and SA methods are tested in this study because of their quick convergence performances [10,11] and proven ability to avoid any existing local minimums in searching space with variation of objective function as the main part of each optimization method.Two different objective function types are implemented in the scope of the tests for application of TS and SA methods.

Testing heuristic optimisation methods for vibration-based detection of damage
The research methods used are presented in a block diagram, as shown in Figure 1.It shows the procedures that are required for detecting damage within structures.

Adopted heuristic damage detection methods
The aim of any optimization method is to find an optimum or near-optimum solution using a reasonable number of iterations, in the range involving 5 % to 20 %out of the total number of trials (points) in the searching space.Therefore, the adopted optimization method has to manage and find an optimum or near optimum solution, while exploring just a part of the available domain (searching space), within the lowest possible number of iterations.Heuristic methods, as one of the most powerful hard non-linear optimization tools, exhibit a very good performance in solving this type of problem.

Tabu search optimization method
Tabu Search (TS) is a meta-strategy for combinatorial optimization problems, Glover [17], which allows random selection of possible solutions during iterative procedure.
It is classified as a practical tool for solving hard nonlinear optimization problems and could be considered as an improvement to the hill-climbing optimization method.In order to avoid repetition of solutions in the search process, recent solutions are stored in the so-called short-term memory.The algorithm could be improved by adding the longterm memory for the same purpose, which enhances the performance simultaneously within iterations [10].Both memory types contain lists of forbidden solutions memorized in the so-called Tabu-lists.The adopted TS procedure, with an advanced long-term memory performance, is summarized in Figure 2.
The algorithm starts with the definition of the iterations termination criterion.It has to be defined based on the nature of the problem, maximum number of iterations according to allowable computing, and satisfactory level of objective function improvement over a certain number of subsequent iterations.It is assumed that the termination (stopping) criterion is a part of the total number of possible solutions in the searching space, ranging from 5 % to 20 %.The maximum number of iterations for TS application was N T TS = 120, which represents 14.3 % of the total number of solutions In order to investigate the entire searching space (set of all possible solutions), it is useful to initially divide it into a certain number of subsets, R, which was adopted as 10 in the study.The long-term memory, for storage of all solutions checked during the iteration process, has to be initialised prior to the process.Also, the initial value of best approximation of the objective function optimal value, f BEST , has to be assumed, as a large number for the case of objective function minimization, as well as the solutions which correspond to the optimal value of the objective function as a trivial solution, S*.The iteration process is performed in three loops: the main loop, loop over subsets, and loop over chosen solutions.The main (outer) loop over iterations corresponds to the total allowable moves, N MOVE , in each of the created subsets.The short-term memory memorizes (stores) just checked (visited) points in the actual move, and has to be initialised as zero vector at each move.The searching process in actual move begins randomly by choosing the initial solution, over all subsets, (S 0 ) j , j = 1, 2, 3,...R.In addition, a certain number of K neighbouring solutions (trials), N(S i ) j , i = 1, 2, 3,...K, has also been randomly selected around the initial solution of each set.Three neighbouring solutions (points) have been adopted, and they form the socalled candidate lists, N(S i ) j of each set j = 1, 2, 3,...R.For all randomly-generated solutions of each set, either initial or neighbouring, it is necessary to determine their presence in the short and long-term memory.If better approximation of the objective function value appears at some point, the best approximation, f BEST , as well as the optimal solution S, have to be updated.The iteration process terminates when the maximum number of iterations is reached, which is manifested Ahmed Alalikhan, Saad Al-Wazni, Zoran Mišković, Ratko Salatić, Ljiljana Mišković by the number of checked solutions, N T TS = R • K • N MOVE , in the adopted TS method.

Simulated annealing optimization method
The Simulated Annealing (SA) heuristic optimization procedure starts by initializing the value of the objective function F(x initial ), which corresponds to a randomly selected initial point in the searching space, x initial .The initial solution depends on one or more initial values of variable(s) or parameter(s), related to the structure of problem, which also have to be chosen at random.The total number of iterations N T SA , has to be proposed with respect to all possible solutions in the searching space.Hence, the initial solution x initial has to be assumed as the current solutionx current , which is generated randomly.The process of searching for a better solution begins in the neighbourhood of the current solution, during the so-called cooling schedule.A set of new variable(s) is generated randomly from searching space in the iteration, and the new solution is defined as x new = x j .The new value of the objective function F j (x new ) comes from the vector of new parameter values, x new in the j th iteration.The new solution x new will or will not be accepted, depending on the acceptance probability p, defined by Equation ( 1), [11]: The acceptance probability depends on the internal energy level difference, defined by the objective function values, ∆F = (x new ) -F(x current ) = ∆E.If F(x new ) results ∆E ≤ 0, the new state, which corresponds to x new vector of variable(s), will always be accepted, and the current solution will be updated, x current = x new .However, if F(x new ) results in ∆E ≤ 0, the acceptance of the x new state depends on the acceptance probability defined by p = e -∆E/T at a particular annealing temperature T, compared to the randomly generated parameter r in the range of [0,1].If p > r the x current = x new will be accepted, otherwise the new state, defined by x new , will be rejected.According to acceptance probability, there is a high probability that the change of variable(s) will be accepted, if the annealing temperature T is high, and vice versa.Thus, the SA algorithm accepts some changes of state without improvement of the solution, in order to avoid being trapped in a local minima.In case the new state is accepted, x new becomes x current in a new iteration, and F(x current ) = F(x new ).In the next iteration, the annealing temperature T is going to be decreased, as a starting parameter of the proposed so-called cooling schedule during the iterations, which reduces the acceptance probability p.An important issue of SA algorithm is the initial temperature T 0 , as an initial parameter of the algorithm and cooling schedule during iterations, which defines decrease in temperature and reduction in acceptance probability.Some authors [12] proposed various approaches for determining initial temperature and cooling-schedule.Those approaches could be defined either by linear, exponential or logarithmic relation.In the given study, linear relation is used with the reduction factor of 0.8 between two subsequent iterations, and the initial temperature parameterT 0 = 20.This terminates when the objective function reaches a certain level of convergence (adopted as 0.1) between two subsequent iterations.Termination criteria have to be viewed as a part of the total number of possible solutions in the searching space.In case of the SA procedure, it is assumed that the maximum number of iterations is N T SA = 200.

Objective functions in proposed damage detection procedures
The objective function is an extremely important segment of every optimization procedure because it has to enable explorations aimed at finding global solution to the problem.It has to take into consideration the most important parameters, according to the structure of the problem.Generally, in case of a vibration-based damage detection task, an objective function has to represent the difference between the experimentally estimated and numerically simulated modal properties.In the damage detection process, the objective function usually takes into the account different parameters as a measure of compatibility of results, such as modal frequencies, mode shape vectors, MAC (Modal Assurance Criterion) values, strain energy residuals, and flexibility matrices [14].Researchers use different combinations of these fundamental parameters.Two forms of the objective function have been adopted in the proposed damage detection procedures.

Objective function forms adopted in this study
The first objective function form adopted in the part of the study for the application of TS optimization for damage detection, includes three parameters: relative differences in modal frequencies, difference in normalized mode shape vectors, and modal assurance criteria (MAC) values.These parameters have been selected due to their high relevance to the subject under study, [5].
Participation of error in the frequencies between experimentally estimated and numerically simulated values is taken into the account according to Equation ( 2). ( The symbols f i E and f i N represent natural frequencies estimated based on the experimental investigation of the structure and numerical computations, respectively, while n represents the number of modes.

Testing heuristic optimisation methods for vibration-based detection of damage
The differences between mode shapes are taken into account according to two parameters.The first one represents the difference between the absolute values of displacement components of all modal vectors examined, Equation (3), assuming that such differences are minor. (3) , in Equation ( 3), represent the normalized reduced mode shape vector, with components solely in vertical direction at measurement points, obtained from experimental and numerically computed data, respectively.The measure of the mode shape vector difference D is based on the assumption that there are slight differences between the same sign of the mode shape vector components.Symbol m denotes the total number of measuring points (transducers) located at particular nodes of the FE structural model.The orthogonality of mode shapes is represented with the parameter MAC [14].The participation of diagonal MAC ii values, the significant indication of mode shapes compatibility, could be taken into the account according to Equation ( 4), which represents the second parameter of the mode shape comparison. (4) The first form of the objective function adopted for application of TS method is defined by Equation ( 5). (5) The factors W F , W D and W M n Equation ( 5) are the so-called weighting factors of each of the previously defined three parts of the objective function, Equations ( 2), ( 3) and ( 4).The second proposed form of objective function for the application of SA method, includes slight modifications of the objective function form used in [5], Equation ( 6).(6) According to the targets in damage detection procedures, detection extent and location of damage, an accurate mathematically defined standard is adopted as a measure of difference between mode shape vectors.
The standard N of differences in reduced mode shape vectors is defined by Equation (7), with components as previously explained in Equation ( 3).An advantage of this type of measure of displacement differences is that it covers a possibility of different signs between displacement components. (7)

Tested structural model
The study was carried out based on a simply supported steel overhang beam as a case for testing damage detection procedures.The total length of the beam is 1500 mm and it measures 5 x 50 mm in cross section.The 20 mm overhang was made on the left-hand side to provide support, while the overhang of 380 mm was made on the right-hand side to make the structure asymmetrical.The mid-span length between the supports was 1100 mm, Figure 3   The element length is assumed to enable identification of the precise damage location for the adopted cases of damage.Numerical analysis was conducted by means of ANSYS package for structural analysis using the Beam4 element type and a linear elastic isotropic material.Due to the low mass of the model, masses of instruments and supporting cubes were included in the FE model.Masses of seven transducers used during the experimental investigations were added to the FEM using the mass element type, with an average mass of 34.525x10 -3 kg and an appropriate mass moment of inertia of 5.3945 x10 -6 kgm 2 .The first four modes in the intact state were computed using the initial FE model for assumed initial values of material Ahmed Alalikhan, Saad Al-Wazni, Zoran Mišković, Ratko Salatić, Ljiljana Mišković parameters, Table 1.The corresponding modal frequencies and mode shapes, extracted from experimental data, are also given in the same table.

Initial experimental analysis using ambient vibration
The experimental study was carried out in order to extract modal properties of the overhang beam model excited by simulated ambient vibrations.The model of the structure, with geometric properties defined in Figure 3, is adopted for all tests in the Laboratory for Structures of the Faculty of Civil Engineering, University of Belgrade.Seven accelerometers were installed on the model and placed on cubes glued onto the beam.Ambient vibration was simulated using a music device, bass shaker, which produces low vibrations in the range from 5 Hz to 200 Hz.It was installed on the independent wood beam connected to both supports of the beam model.Low noise, high sensitivity and low frequency accelerometers, Silicon Designs Model 2400 with the range of DC-600 Hz, were used for vibration measurements.The data were acquired via the 24-bits 8 channel HBM -Hottinger Baldwin Messtechnik, QuantumX measuring amplifier.Modal frequencies and mode shapes were extracted using the ARTeMIS -extractor software [15], which is a state-of-the-art program for ambient vibration analysis.Among several currently available extraction procedures, the authors selected the frequency domain decomposition (FDD) technique with required settings of parameters for the analysis (sampling frequency: 600 Hz, number of frequency lines: 1024, number of channels: 7).Several trial measurements were made to provide a greater confidence (5 trials for intact case), with the corresponding estimation of modal properties according to spectral density matrices, using the recorded ambient vibration data.In all cases, the structural modes were estimated using the peak-picking method, according to FDD estimation procedure [13], Figure 6.Frequency

Calibration of FE model
In order to match closer numerical and experimental modal frequencies, FE model was calibrated by tuning various model properties: mass, modulus of elasticity, geometry, etc., [16], which affect modal properties.Also, the differences between experimental and numerical analyses are the consequence of the noise that is constantly present in the recorded data, as well as imperfections of the real structure or structural model.The calibration can be conducted either manually by a parametric study, or automatically, using the SA optimisation procedure.Both applied approaches, manual model tuning by parametric study, and model updating using the SA optimization procedure, resulted in practically the same results of the varied parameters.

Calibration of initial FE model by parametric study
It is assumed in parametric study that just first three modes could exhibit enough indication for the sake of damage detection.So, the first three modal frequencies were selected according to FDD modal estimation technique, as peak values should be adopted as the requested mode.Numerically computed initial modal frequencies listed in Table 1 could be adjusted by parametric study to get closer values with respect to experimental results.According to high sensitivity of the modal frequencies with respect to the change of elasticity modulus, a variation of the parameter was observed, in the range of 5 % of the initial value.Table 2 shows the first three computed modal frequencies for modulus of elasticity values varied in the proposed range.Also, the difference between the computed and experimentally estimated frequencies for the first two modes is shown in the same table.The best agreement between the computed and experimentally estimated values is reached for the modulus of elasticity of E optTS = 2.09 • 10 5 N/mm 2 , according to the sum of absolute differences of the first two best tuned modal frequencies (Equation (8). Error The third modal frequency was not taken into consideration because, unlike the first two, it exhibits divergent behaviour.
Figure 7 shows the dependence between modal frequencies and modulus of elasticity.By using the optimal value of the modulus of elasticity, the correlated FE model was implemented for further numerical computations and simulations, with the application of the TS optimization procedure in damage detection.

Calibration of initial FE model using SA optimization procedure
As can be seen in Figure 4, the initial FE model correlation was conducted by means of the SA optimization algorithm for application of the damage detection procedure via SA optimization.The optimisation was made by varying two parameters, modulus of elasticity and thickness of beam section, in the range of 5 %.The choice of varied parameters

Mode
No. been significantly reduced for the first two modal frequencies, to less than 0.5 %.

Implemented damage
Two cases of damage, differing by damage extent, were implemented to investigate the structural model behaviour and the efficiency of the proposed procedure for detecting small and large crack/notch occurrences.The location of damage was fixed after the study of few positions on the calibrated FE beam model by numerical analyses, in order to find the location with the most significant impact on modal frequency changes.The damage was located at the distance of about 400 mm from the left edge of the beam, close to one of the measuring points, while the studied damage cases exhibit different levels of damage.Two damage cases were created by cutter machine from the bottom side of the cross section.The first one was a small damage, which afterwards extended, thus generating the second damage case.

Damage cases
The first damage case, DC-1, consists of a crack 1.9 mm in depth and 25 mm in width on the structural model, situated 387.5 mm away from the left support, Figure 9.The position of damage coincides with the element number 17 of FE model, with element numbering from the left end of the beam, Figure 3.
GRAĐEVINAR 68 (2016) 7, 543-557 Ahmed Alalikhan, Saad Al-Wazni, Zoran Mišković, Ratko Salatić, Ljiljana Mišković was made based on the investigation of their uncertainty and high sensitivity of natural frequencies of those two parameters.The total mass of the beam, as well as additional masses, were accurately measured, and no variations were registered.In addition, the width of the beam section did not significantly affect the natural frequency values for minor changes.
The objective function for FE model calibration is represented by the sum of absolute differences between experimental and numerical natural frequency values, Equation ( 8), including the first two modes.The correlation improvement during SA iterations is shown in Figure 8.The best correlated results correspond to the optimum modulus of elasticity of E opt-SA = 2.09 • 10 5 N/mm 2 , and optimum cross-sectional thickness of t opt-SA = 4.99 mm.For these values, the first four natural frequencies of the correlated FE model increased by about 2 % with respect to the initial FE model, Table 3.
Better agreement between the calibrated FE model results and experimental results is evident.The differences have

Table 3. Computed natural frequencies by correlated FE model after application of SA optimization for parameters calibration and comparison with experimentally estimated values
Testing heuristic optimisation methods for vibration-based detection of damage Further extension of the damage produced the second case of damage, DC-2, by increasing the crack depth from 1.9 mm to 2.75 mm at the entire damage width of 25 mm, in the same position, Figure 9.

Experimental analysis of damaged beam based on ambient vibration measurements
Three trial ambient vibration measurements were carried out for each of the two damage cases in order to improve reliability of the estimated modal properties.Estimated modal frequencies for both damage cases and the first four modes are listed in Table 4.
Compared to the intact case, different frequency values were extracted in the case of DC-1 for the second, third and fourth modes from experimentally recorded data, while in case of DC-2 different frequency values were extracted for all four first modes with more significant differences.In the first case, it was more difficult to detect the position and extent of damage than in the second case (DC-2), which is mostly due to insignificant differences, Table 4.

Tests of objective function sensitivity on studied damage cases
Objective function efficiency tests were conducted prior to the application of the proposed damage detection procedures.Both proposed objective function forms, defined by Equations ( 5) and ( 6), were analysed in order to find the most effective weightingfactor values for different severities of damage.
In the present study, the severity of damage is represented by the ratio of crack_depth/overall_depth, R D -damage ratio parameter, as defined by Equation ( 9). ( Symbol d in Equation ( 9), represents crack depth from the bottom side of the cross section, while t represents the overall depth of the cross section.The stiffness of the entire element has to be reduced during numerical simulation of damage.Hence, in order to get an adequate stiffness, because the crack width could be smaller than the element length, the parameter a needs to be introduced to weigh an adequate stiffness of the element and its real stiffness with narrow crack, Equation ( 9).This parameter takes value in the range of 0 < a ≤ 1.
It is assumed that a = 1 , because the width of the actual crack was the same as the element length (25 mm).The actual damage ratios were computed according to Equation ( 9), Table 5.  11.Participation of weighting factors for the studied damage case DC-1 is shown in Figure 12.Significant influence of mode shape vectors can be noted.The same study was carried out for a more serious damage case -DC-2.The damage ratio varied between 0.4 and 0.7, while the actual damage ratio of DC-2 amounted to 0.55.In this case, the objective function, with the weighting factor values assumed as in the previous case, exhibits similar behaviour, Figure 13.For the actual damage ratio of 0.55 and its position at element No.17, the objective minimum value was reached at the position of the element number 17, exactly at the actual position.In addition, minimums achieved for other examined damage ratio values were very close to the actual position, Figure 13.These values were used for computation of effective stiffness of the damage element during numerical simulations.

Tests of objective function for the use of TS optimization in damage detection
In order to ensure that the adopted objective function defined by Equation ( 5) will be sensitive enough to any changes in modal properties due to existence of damage, a parametric study was conducted to examine the influence of weighting factor values.The study was carried out using the FE structural model, Figure 4, with calibrated properties.The first three modes were taken into account based on the assumption that they contain enough information on the change of modal structural response due to damage.Table 6 shows variation of weighting factors for the first damage case, DC-1.The most suitable values in terms of efficiency of the objective function are W F = 1, W D = 1 and W M = 1, for differences in frequencies, normalized mode shapes, and main MAC values, respectively.Testing heuristic optimisation methods for vibration-based detection of damage Participation of weighting factors in the observed damage case DC-2 is shown in Figure 14, with significant influence of mode shape vectors, as in the previous case.

Objective function tests for SA optimization in damage detection
As in the previous case, the parametric study of the weighting factor influence on objective function efficiency for SA optimization, Equation ( 6), was conducted to determine the most suitable values.The studied values of weighting factors W F and W D varied in the range shown in Table 7.According to the results given in Table 4, the first four mode shapes were taken into consideration in the parametric study.In the observed damage case DC-1, with small differences in frequencies in the range of 3 %, the differences in displacements (mode shapes) with adopted weighting factors participate more significantly in the total value of the objective function, compared to natural frequency differences.The most suitable weighting factor values are practically explored to amplify the participation of the small frequency differences.Also, the displacements have to be included to get the minimum value of the objective function at the exact position of damage, at element No. 17, as shown in Figure 15.
As can be seen in Table 5, the same analysis was conducted for the adopted DC-2 damage scenario with a more significant damage ratio of R D = 0.55, using the same proposed weighting factor values, W F = 10 and W D = 1.The behaviour of the adopted modified objective function shows that the proposed weighting factors are adequate and that they enable attainment of the global minimum at the exact position, i.e. at element No.17, Figure 16.Further investigation was conducted to examine the efficiency of objective function defined by Equation ( 6), with the previously defined most suitable weighting factor values and different Value of the weighting factor 0 / 1 / 10 / 100 0 / 1 the actual damage location and severity with its global minimum.

Proposed damage detection procedures
Using the TS and SA optimization, the proposed procedures were implemented by routines developed in MATLAB, which is a fundamental tool for scientific computing.They include automatic numerical simulation of damage by calling the ANSYS FE package for computation of modal properties.During the optimization process, the main routine modifies the bending stiffness of the model by selecting the damaged element and damage ratio at the current stage [8].Furthermore, it computes the objective function value based on updated modal properties computed by FE package.The proposed damage-detection procedures were applied for both damage cases, DC-1 and DC-2.

Damage detection by TS optimization
The proposed objective function and weighting factor values, Equation ( 5), are used in damage detection by TS optimization for the observed damage cases.The adopted searching range of damage ratio was R D = [0.55,0.7] with an increment of ∆R D = 0.05, i.e. with 14 different discrete values.Hence, for the total number of 60 elements in the FE model, the number of possible solutions in the searching space was 840.The maximum number of iterations was limited to 120, which is within a reasonable range of 20 % of total points in the searching space.The proposed procedure included the first three modes, frequency and transverse displacement vector of mode shapes.The improvement of the objective function global minimum value approximation is shown in Figure 19 for the adopted DC-1, and the optimum value was reached after 19 iterations.The exact position of the damaged element No. 17 was explored with the damage ratio of R D * = 0.4, which is the best approximation of the actual damage ratio implemented in DC-1 for the adopted accuracy of ∆R D = 0.05.In DC-2 damage case, an optimum solution was reached after 120 iterations, Figure 20.In this case, the exact optimum solution explored the damage element No.17 with the damage ratio of R D * = 0.55, because the assumed accuracy of ∆R D = 0.05 was appropriate to locate the exact position of the objective function minimum.
The procedure was repeated several times, and it demonstrated the same results in all cases under the limited number of iterations.

Damage detection by SA optimization
The SA optimization, with the proposed modified objective function form and weighting factor values, Equation ( 6   The study shows that the modified objective function, Equation ( 6), with the adopted weighing factors, reflects Testing heuristic optimisation methods for vibration-based detection of damage applied for damage detection in both damage cases adopted in this study.The damage ratio was considered to be a discrete variable with the accuracy of 0.01 in the range of R D = [0,01, 0,80].
While the number of elements in the FE model was 60, Figure 4, the total number of all possible solutions in the searching space was 4800.Unlike the TS optimization procedure, the searching space was significantly extended in order to improve the damage ratio determination accuracy.Due to the width of the searching space, the maximum number of iterations was limited to 200, i.e. to 4 % of the total number of points in the searching space.The proposed procedure included the first four modes, frequency and transverse displacement vector of mode shapes.The improvement of solution during SA optimization, with evident quick convergence, is shown in Figure 21.The optimum solution was reached after 184 iterations, which fully corresponds to the actual DC-1, damaged element no.17, with the damage ratio of R D * = 0,38.The points generated in the searching space, Figure 21.b, represent their adequate random distribution.It covers the whole searching space that is needed to explore the global optimal solution and avoid local minimums.
In case of DC-2, the proposed damage detection procedure based on SA optimization exhibits quick convergence, Figure 22.a, and an optimum solution is reached after 179 iterations.In this case, the so-called near-optimal solution of R D * = 0,56, was reached, which is very close to the actual damage ratio with the exact damage position, element No.17.Just like in previous case, the checked solutions are randomly well distributed in the searching space, Figure 22.b.Very good agreement between the implemented and estimated damage location and ratio points to the robustness of the proposed procedure and quick convergence during the limited number of iterations.

Conclusion
Two heuristic optimization methods, TS and SA, are used in this study for detecting structural damage based on the changes of modal properties of structures.Ambient vibration measurements and FDD method of extraction of modal properties from a simple structural model are used in accordance with a state-of-the-art modal identification procedure for large civil structures.
Using a parametric study and automatic updating by an optimization procedure, the FE model is calibrated according to the modal properties estimated from experimental analysis in the intact (undamaged) state of the structure.It resulted in modal frequencies that were close to experimental values, with differences ranging from 2 to 3 %, which is the level of differences that could be achieved in case of real civil structures.According to the results, it would seem advisable to use more parameters in the calibration process by applying an automatic updating based on heuristic optimization.It could be particularly useful in case of more complex structures, where variation of a greater number of parameters could be made.The proposed weighting factor values for both adopted objective function types reflect the damage state of a simple structural model, which was experimentally investigated in this study.
The damage detection procedure is based on TS optimization, according to the ability to explore the position and severity of damage under a limited number of iterations, which enhances good performance of this technique for the applied accuracy.In both cases, close or similar severity of damage was explored for the adopted accuracy and at exact position.The proposed damage detection procedure based on the SA optimization technique exhibits high accuracy and quick convergence in a very large searching space.In addition, damage was explored exactly and very closely for both cases of damage, in accordance with the adopted accuracy.
A slightly better convergence of SA compared to TS optimization was observed during the research within a smaller percentage of checking solutions with respect to the total solutions in the searching space.Both proposed procedures exhibit excellent agreement with experimental data.These first findings suggest that the adopted procedures could be applied on real structures.Nevertheless, more tests should be carried out on complex structural models to prove high potential of the adopted heuristic optimization techniques in damage detection.In addition, as there are no restrictions with respect to the use of multiple damage scenarios, the performance of the adopted techniques should also be studied for more cases of damage, starting with simple structures.

Figure 1 .
Figure 1.Chart of overall research procedure: adopted heuristic optimization, tuning of objective weighting factors, finite element model correlation, and damage detection procedure

Figure 2 .
Figure 2. Flow chart of Tabu search optimization method adopted for damage detection . Initial material properties are: modulus of elasticity E initial = 2x10 5 N/mm 2 , Poisson's ratio n = 0.3 and mass density m initial = 7.860 x10 3 kg/m 3 .

Figure 4 .
Figure 4. Overhang FE model for damage detection procedure

Figure 7 .
Figure 7. Differences in first three natural frequencies between experimentally estimated and numerically computed frequencies for different values of modulus of elasticity E

Figure 9 .
Figure 9. Tested overhang steel beam model with damage and transducer locations

Figure 8 .
Figure 8. Improvement of FE model correlation with respect to differences in first four modal frequencies by SA optimization procedure

Figure 10 .
Figure 10.Size of damage in DC-2

FigureFigure 13 .
Figure Influence of weighting factors in objective function for TS application, Equation (5), for studied damage scenario DC-1: R D = 0,38, and variation of damage position

Figure 11 . 1
Figure 11.Objective function graph with damage location along the beam for different damage ratios -simulation by FE correlated structural model and DC-1

Figure 14 .
Figure 14.Influence of weighting factors in objective function for TS application, Equation (5), for studied damage scenario DC-2: R D = 0,55 and variation of the position of damage

For
the first case of damage, DC-1, with the damage ratio of R D = 0,38, the most suitable values of weighing factors W F = 10 and W D = 1 were found in the modified objective function, Equation (6), which results in the exact minimum at the actual position of damage, Figure 15.

Figure 15 .
Figure 15.Influence of weighting factors in objective function for SA application, Equation (6), for observed damage scenario DC-1: R D = 0,38, and variation of positions of damage

Figure 16 .
Figure 16.Influence of weighting factors in objective function for SA application, Equation (6), for studied damage scenario DC-2: R D = 0,55, and variation of the position of damage ), was also GRAĐEVINAR 68 (2016) 7, 543-557 Ahmed Alalikhan, Saad Al-Wazni, Zoran Mišković, Ratko Salatić, Ljiljana Mišković severity of damage.First, the study was carried out for the adopted scenario of damage DC-1, by varying the damage ratio R D in the range of [0.1, 0.6] and the position of damage along the beam.Figure17shows the change of the objective function during numerical simulation of damage for the actual DC-1 results.The global minimum of the objective function achieved for the damage ratio of R D = 0,38 at the exact damage position, element No. 17, with more local minimums, is shown in Figure17.

Figure 17 .
Figure 17.Objective function for SA application, numerically simulated for DC-1, with global and local minimums

Figure 18 .
Figure 18.Objective function for SA application, numerically simulated for DC-2, with global and local minimums

Figure 21 .Figure 20 .
Figure 21.Improvement of objective function global minimum approximation during a) SA based damage detection procedure for implemented damage case DC-1; B) with distribution of solutions checked in searching space

Figure 22 .
Figure 22.Improvement of objective function global minimum approximation during: a) proposed SA based damage detection procedure for damage case DC-2; b) with distribution of checked solutions in searching space

Table 1 . Experimental and numerical modes in intact state with initial values of material parameters
GRAĐEVINAR 68 (2016) 7, 543-557 Testing heuristic optimisation methods for vibration-based detection of damage

Table 4 . Experimentally estimated modal frequencies of studied model in intact state and with implemented damages
value of the objective function for the proposed weighting factor values coincides with the damage position, element No.17.For other examined damage ratio values, minimum global positions are in the neighbourhood of the exact position, Figure minimum