Optimized Erosion Prediction with MAGA Algorithm Based on BP Neural Network for Submerged Low-Pressure Water Jet

Featured Application: Underwater cleaning and rock breaking operations of ocean engineering. Abstract: In order to accurately predict the erosion e ﬀ ect of underwater cleaning with an angle nozzle under di ﬀ erent working conditions, this paper uses refractory bricks to simulate marine fouling as the erosion target, and studies the optimized erosion prediction model by erosion test based on the submerged low-pressure water jet. The erosion test is conducted by orthogonal experimental design, and experimental data are used for the prediction model. By combining with statistical range and variance analysis methods, the jet pressure, impact time and jet angle are determined as three inputs of the prediction model, and erosion depth is the output index of the prediction model. A virtual data generation method is used to increase the amount of input data for the prediction model. This paper also proposes a Mind-evolved Advanced Genetic Algorithm (MAGA), which has a reliable optimization e ﬀ ect in the veriﬁcation of four stand test functions. Then, the improved back-propagating (BP) neural network prediction models are established by respectively using Genetic Algorithm (GA) and MAGA optimization algorithms to optimize the initial thresholds and weights of the BP neural network. Compared to the prediction results of the BP and GA-BP models, the R 2 of the MAGA-BP model is the highest, reaching 0.9954; the total error is reduced by 47.31% and 35.01%; the root mean square error decreases by 51.05% and 31.80%; and the maximum absolute percentage error decreases by 65.79% and 64.01%, respectively. The average prediction accuracy of the MAGA-BP model is controlled within 3%, which has been signiﬁcantly improved. The results show that the prediction accuracy of the MAGA-BP prediction model is higher and more reliable, and the MAGA algorithm has a good optimization e ﬀ ect. This optimized erosion prediction method is feasible.


Introduction
Marine fouling is a serious problem in ship operations. The marine fouling adhering to the hull can increase navigation resistance and fuel consumption, up to 40% [1]. Submerged cavitation water jet technology is one of the cleaning methods that can be applied in the underwater cleaning of the ship hull. Submerged water jet refers to the water jet generated where the working and the surrounding media are the same. The high-speed water jet formed by the nozzle under water interacts with the low-speed water of the surrounding medium, generating strong disturbances, and forming periodic vortex rings on the periphery of the jet. Both the low-pressure region of the vortex center The prediction method of backpropagation (BP) neural network is similar to a black box that does not need detailed material constitutive parameters. The prediction model is established on the basis of the representative experimental data and BP neural network. Li et al. [23] used BP neural network to estimate the energy consumption of liquefied natural gas buses. Wei et al. [24] predicted fuel consumption with BP neural network to reduce unpredictable fuel for flights. Ning et al. [25] used the BP neural network to predict the temperature of equipment out of satellite. Wang et al. [26] established a model for predicting soil bulk density based on the BP Neural Network. It can be seen that the BP neural network has a wide range of applications in prediction [27]. However, the output of the network depends on the network input and weights matrix, and the final training result is obviously affected by the initial weights. Therefore, in order to improve the prediction performance of BP neural network, on one hand, when the number of tests is limited, a hybrid method of orthogonal experimental design and virtual data are used to obtain representative test data, which are sufficiently accurate and effective, as the input of prediction model. On the other hand, its initial weights and thresholds are also optimized. Han et al. [28] used the GA-BP algorithm to improve the accuracy of thermal fatigue failure prediction; Zhu et al. [29] relied on GA-BP neural network to forecast the short-term traffic flow at the intersection, and the prediction accuracy increased by 77%. Wang et al. [30] used the GA-BP neural network to predict the ablation wear value of artillery barrel, which greatly reduced the relative error. Tang et al. [31] applied the GA-BP algorithm in the prediction of catalyst volume in the SCR denitrification system, and the prediction accuracy and data fitting ability were improved. It can be seen that the optimized GA-BP model can indeed effectively improve the prediction accuracy of the BP neural network, but the GA optimization algorithm is liable to fall into local points. Therefore, an advanced GA algorithm based on mind evolution (MAGA) is proposed to improve the prediction performance of the BP neural network, making the prediction more reliable.
Most of the current cleaning methods are manual. Workers rely on their working experience to adjust the working conditions at will, which is not reliable and causes a lot of undesired damage. This manual method has a high cost and low efficiency, which will cause a considerable waste. In addition, divers are very difficult to work underwater, and there are great safety risks. Therefore, the ship hull underwater cleaning robot technology has greater advantages [1]. This paper is a part of the research on cleaning technology in the ship hull underwater cleaning robot project based on our patent. In the project, submerged cavitation water jet technology is adopted as the cleaning method. The cleaning equipment is mounted on the robot, and a hybrid adsorption method based on thrust adsorption is used to keep the robot staying on the hull surface for underwater operations [32]. Because this process depends on hydrodynamic performance, adsorption force and hovering attitude of an underwater robot, higher jet pressure will increase the complexity of hydrodynamic performance, make the adsorption of the robot on ship hull more difficult, cause the robot's unstable attitude by excessive reaction forces, increase energy consumption and improve the difficulty of control. Therefore, low-pressure conditions are selected for this research.
The submerged water jet is a cavitation water jet. Although the cavitation water jet is at low pressure, the impact pressure of the cavitation water jet is 8.6-124 times that of the continuous water jet [33]. So even when the pressure is lower, the damaging effect of higher pressure can be achieved. In underwater cleaning, higher pressure water jets may not get good erosion results. Due to the large damage ability of cavitation and lack of effective erosion control, this cleaning method is difficult to control cleaning effect, and can cause certain damage to hull surface and anticorrosion coating during the process of removing marine fouling [1,7,8]. This damage can cause the steel layer of ship hull to be directly exposed to seawater, which will accelerate the corrosion process of steel, promote hull aging, increase maintenance costs, and reduce operating life of the ship. In addition, when the anticorrosion coating is damaged, toxic copper compounds will be released [7] and pollute the marine environment, which is prohibited in many countries. Therefore, the cleaning operation must be accurately controlled by the optimized erosion prediction to ensure that while removing marine fouling, the hull and anticorrosion coating are not damaged. However, there is almost no comprehensive consideration Appl. Sci. 2020, 10, 2926 4 of 25 of erosion control of underwater cleaning in the existing research works. Only a few studies [7,8,34] have analyzed the erosion effect of underwater cleaning for a particular situation under one or several working conditions, and selected a constant working condition for cleaning. Considering that the situations in underwater cleaning areas are variable, the hull surface is not flat and the thickness of marine fouling is different, as shown in Figure 1, a single constant working condition in these existing studies cannot always be appropriate. The unsuitable erosion will waste power consumption of robots, affect the cleaning effect, and even damage the anticorrosion coating and hull surface. Therefore, the erosion prediction under different working conditions in this paper can be better used to achieve real-time erosion control of underwater cleaning robots. One of the essential factors to the successful application of underwater cleaning robot technology in actual cleaning operation is to control erosion to avoid the damage to anticorrosion coating and hull surface. Although this erosion control will slightly increase the development cost of robots, when the underwater cleaning robot is successfully applied in the ship underwater cleaning, the benefits will far outweigh the increased development costs. The underwater cleaning of the ship hull has been approved by several classification societies to extend the dock repair interval to two cycles [1]. Dock repair requires significant economic and time cost, while underwater cleaning can save 60% of dock repair cost [8]. In addition, it is estimated that a single underwater cleaning can improve fuel efficiency by 6%. For ships with an annual fuel consumption of more than 10 million tons, the annual energy-saving benefit of underwater cleaning by robots can reach 200 million CNY [35]. In terms of global shipping, underwater cleaning by robots can save 15 billion USD per year and reduce one billion tons of greenhouse gas emissions from fleets [8].
In order to improve the prediction accuracy and overcome these problems that the difficulty in obtaining fouling samples, testing all working conditions, showing obvious damage at low pressures, and defining target material constitutive models, this paper uses refractory bricks to simulate marine fouling as erosion target, and erosion depth as the output index of the prediction model. Erosion tests are carried out as the test scheme designed by the orthogonal table. Through statistical analysis of range and variance, the jet pressure, impact time and jet angle are three main factors and are used as the input variables of the prediction model. Based on the experimental data, an artificial virtual data generation method is used to increase the input data for the prediction model. Then, an improved optimization algorithm MAGA is proposed. After verifying the optimization effect of MAGA with four stand test functions, the initial thresholds and weights of the BP neural network are optimized by MAGA. Thereby, the erosion effect of the submerged low-pressure water jet under different working conditions is predicted. Finally, the optimized erosion prediction method can be applied to the formulation and control of the hull underwater cleaning scheme, which is an important part of the ship hull underwater cleaning robot project. This paper is organized as follows: Section 2 shows a test method based on orthogonal experimental design to obtain the data needed in the study, and uses a hybrid analysis based on range and variance to determine main factors affecting erosion results, so that the main factors are used as input variables of the prediction model. Section 3 presents an improved optimization algorithm and validates the improved optimization effect of this algorithm with the certain set parameters. Section 4 describes the prediction model established by the input variables in Section 2, and uses the optimization algorithm in Section 3 to optimize the initial thresholds and weights of the network model. Section 5 discusses the prediction accuracy and reliability of prediction models, proving that the optimized erosion prediction method is feasible. Finally, Section 6 shows the results and applications of this paper with recommendations for future research.

Erosion Test Setup of Submerged Water Jet
In the submerged water jet erosion test, the working water medium in the storage water tank (Figure 2a) enters the high-pressure pump (Figure 2b) for pressurization. Then working water is delivered to the nozzle (Figure 2d), which is completely immersed in the operating water tank (Figure 2c) below the free surface of the ambient water. After working water passes through the nozzle, it forms a high-speed water jet and impacts the surface of the refractory brick. The damage is shown in Figure 3. The jet pressure is controlled by a pressure regulating valve and adjusted according to different working conditions. The bottom of the refractory brick is fixed on the adjusting base, which is convenient for adjusting the jet angle and standoff distance. Experimental equipment are as follows: (a) High-pressure pump: the rated power is 15 kW, and the rated flow is 41 L/min. (b) Pressure gauge: the measurement range is 0-40 MPa, and the accuracy level is 2.5. (c) Operating water tank: acrylic glass material, size is 60 × 60 × 50 cm. (d) Nozzle: stainless steel material, angle nozzle is used for better erosion performance. (e) Refractory brick: size is 230 × 55 × 65 mm, temperature resistance is 1350 • C, composition content is 40% alumina, 55% silica.

Orthogonal Experimental Design
The orthogonal experiment design method was proposed by Genichi Taguchi [36], and it is widely used [37][38][39]. The core of this method uses the orthogonal table to uniformly select representative conditions from all factors and levels instead of considering all the combinations. Without comprehensive tests, the effects of various factors can be analyzed to find the optimal level combination, which greatly reduces the test cost.
According to the existing research of Soyama et al. [10] and Kang et al. [12], factors such as jet pressure, jet angle, standoff distance, and impact time will have certain effects on the experimental results. They are also the setting parameters of the working conditions that can be controlled by the test equipment during the test, and have a major influence on the control of underwater robots. In order to further understand their effect on test results, a hybrid orthogonal test design L 16 (4 3 × 2 1 ) is selected for range and variance analysis. According to the method of orthogonal test design, three factors, namely jet pressure (p), jet angle (α), and impact time (t), have four levels, respectively. The four levels of jet pressure at low-pressure conditions are 2 MPa, 4 MPa, 6 MPa, and 8 MPa. The four levels of jet angle are 30 • , 45 • , 60 • , and 90 • . The four levels of impact time are 20 s, 40 s, 60 s, and 80 s, which can be enough to achieve a certain amount of damage. The factor of standoff distance (s) has two levels, 10 mm and 15 mm. These parameters are determined by considering the practical low-pressure operating conditions and the limitations of the experimental device. The jet pressure is shown on the pressure gauge. The jet angle is a complementary one between the central axis of the jet and the normal of the target surface, and the range is from 0 • to 90 • . The standoff distance is the distance from the nozzle outlet to the target surface. The erosion depth (H) is the maximum distance from the surface of refractory brick to the bottom of erosion pit caused by water jet. The experimental results of erosion depth are used as evaluation indexes to reflect the extent of erosion damage. The deeper the erosion, the stronger the erosion ability of the water jet. The erosion parameters of the water jet test are shown in Figure 4. In Figure 4b), the various variables of working conditions are shown when the water jet impinges on the target surface at a specific angle α. The factors and levels are shown in Table 1, and the orthogonal test results are shown in Table 2.   Table 3 shows the range analysis of the erosion depth of the submerged water jet. Range analysis can determine the influence order of each factor on test results, and obtain the possible optimal level of each factor to design a better test condition. R is the range, calculated by Equation (1), which reflects the fluctuation of test results and the significance of the corresponding factor's effect on the evaluation index. A larger R indicates that a change in this factor will cause a larger change in test results. This means that the influence of this factor on the evaluation index is greater. The actual definitions are: where k i is the average value of the indexes, which represents the average value of test results of the i level of the factor, and reflects the comprehensive effect of the i level on the evaluation index. k max and k min correspond to the maximum and minimum values of k i , respectively. K i is the sum of the indexes, which represents the sum of the results of the i level of the factor. y ij corresponds to the index value H of the j test in the i level. i represents the number of the level. r is the number of repeated tests at the i level of the factor. j is the number of the test from 1 to r. n is the total number of tests. q is the total number of levels of the factor. In the hybrid orthogonal table, the different number of levels of each factor will lead to a greater R for the factor with more levels, although this factor may have the same impact as other factors. Therefore, R is used to correct the range R for the factors with different levels [40], as shown in Equation (2), and the correction coefficients d are shown in Table 4. As shown in Table 3, R p > R t > R α > R s . Therefore, jet pressure has the most significant effect on erosion depth. Impact time and jet angle follow closely, and their correction ranges are relatively small, in which the effect of impact time is greater than the jet angle. The correction range of standoff distance is the smallest, so its change has the least influence on the test index. The factor and index chart can further analyze the influence of various factors in the test ranges on the index of erosion depth, as shown in Figure 5. The trend chart of the factors and indexes are based on the levels of each factor as the abscissa, and the average value k i of the test indexes in Table 3 is used as the ordinate. It can more intuitively show the change trend of the test index with the change of the factor level. For the jet pressure factor, the average value of the test index corresponding to the fourth level k 4 is the largest, so 8 MPa is the optimal level for this factor. Similarly, the optimal levels of impact time, jet angle and standoff distance are 80 s, 60 • , and 10 mm, respectively, so the optimal level combination in this orthogonal test design is (p,α,t,s)ˆ* = (p_4,α_3,t_4,s_1). Since this combination is not in the test list, a supplementary verification test is required. Under this condition, the erosion depth is 2.238 cm, which is indeed larger than the maximum index value of 2.108 cm on the fourteenth test in Table 2. This shows that the optimal working condition by range analysis is effective. In addition, Figure 5a reflects that when jet pressure is in the range of 2-8 MPa, erosion depth increases significantly with the increase of jet pressure. In Figure 5b, with the increase of jet angle, the optimal jet angle is around 60 • , which is consistent with the experimental result of Sun et al. [41]. Figure 5c describes that erosion depth is increasing with the increase of impact time in the range below 80 s. In Figure 5d, standoff distance is within the test range of 10-15 mm, and the erosion depth decreases as the standoff distance increases. This is because the distance for energy exchange between water jet and the surrounding medium has become longer, and the entrainment effect of the boundary layer causes more low-speed surrounding media to be mixed into water jet, which can reduce the erosion ability of water jet. This is in accordance with the results of Liu et al. [42].

Possible Errors a.
Errors from experimental equipment design.
Due to frictional drag in the pipeline, there is a slight pressure loss of the jet. When designing the test device, the length of the pipeline should be minimized, and the pressure gauge should be installed near the nozzle. In order to avoid the influence of operating water tank size on the flow field, a suitable size is adopted. According to the research of [2,17], this size will not affect the cavitation effect of the nozzle.

b.
Errors from experimental materials.
To minimize the effect of working water medium on the cavitation effect, after the storage water tank is filled with water, allow it to settle for two hours before starting the experiment. The same batch of refractory bricks are used as impact targets to ensure that the surface of each refractory brick has the same characteristic, and soaking in water for one hour after cleaning the surface impurities to ensure that refractory bricks contain the same water content. The operating water tank is always kept full of water to ensure that the erosion effect is not affected by the free surface during the test. The depth of the nozzle underwater is also fixed to maintain a constant pressure around the nozzle outlet. c.
Errors from experimental operation.
At the beginning of the test, pump pressure is unstable. Therefore, after the pump attains stability, the jet pressure is quickly adjusted to reduce its effect on test results during pressure changes. When the pressure reaches the setting conditions and the work is stable, test timing is started to reduce the error caused by pressure fluctuations. As the test progresses, more and more debris is washed away, which may affect the effect of cavitation. Therefore, after a large amount of debris appears, the ambient water in the working water tank is replaced to reduce the effect of turbidity on the test results. To improve the accuracy of the test results, multiple tests are performed for each working condition, and the average result is taken. d. Errors from reading data.
Due to the limitation of the accuracy level of the pressure gauge, there is a slight precision error. After the test, the morphology of erosion damage on the refractory brick is irregular. For the maximum erosion depth, multiple measurements with a vernier caliper are required to reduce measurement errors.

Variance Analysis
In order to distinguish test errors from test index variation caused by changes in test conditions, an analysis of variance is used to estimate the test error that must exist during the test and the result measurement. The sum of the square of test error, SSe, was 0.041 by Equation (3). SS T is the total sum of square, SS p,α,t,s is the sum of the square of jet pressure, jet angle, impact time and standoff distance, respectively. With the sum of the square of each factor, the F value and P-value corresponding to each factor can be calculated to determine whether the influence of the factor on erosion depth is significant and reliable. The smaller the P-value, the more significant the effect. In statistics [43], according to P-value is less than 0.05, jet pressure, jet angle and impact time can be considered to have significant and reliable effects on erosion depth, as shown in Table 5. The P-value corresponding to jet pressure is less than 0.01, indicating that the influence of jet pressure is very significant, and much greater than impact time and jet angle. The effect of impact time is also greater than jet angle, which is consistent with the results of the range analysis above. However, the effect of standoff distance is the smallest. Since the P-value of standoff distance is greater than 0.05, this factor can be considered to have no statistically significant effect on erosion depth within the range of 10-15 mm. Therefore, in this paper, jet pressure, jet angle and impact time are three main factors and are determined as the input variables of the prediction model.

Genetic Algorithm
The genetic algorithm (GA) [44] is a single population algorithm. Individuals in the population evolve their excellent genes through selection, crossing, and mutation operations, and enrich the gene pool of each generation to form a new generation of populations. Better individuals emerge with improved individual fitness. When the immature convergence occurs randomly, which drastically reduces the diversity of individuals in the population, and the viable new individuals are difficult to generate by the constant cross and mutation operators to break the limits, the result is easy to fall into a local solution. Therefore, both retaining good individuals and maintaining the diversity of the group should be considered to avoid falling into a local solution prematurely. The basic flow of GA is shown in Figure 6. The GA algorithm steps are as follows: Step 1: The population is initialized to generate individuals carrying specific genetic information; Step 2: The first generation individuals are real-coded to keep them within the specified range; Step 3: The fitness function is determined and the fitness value of each individual is calculated; Step 4: According to the select operator P g , individuals with better fitness values are selected for the next generation; Step 5: According to the cross operator P c , the genes a kz and a lz at the z position of the k and l two individuals in the population are selected for cross operation by Equation (4); where b is a random number between 0 and 1.
Step 6: According to the mutation operator P m , the gene a uz at the z position of the u individual in the population is selected for mutation operation by Equations (5) and (6); where a max and a min are the upper and lower limits of the gene a uz respectively; v 2 is a random number; num is the current iteration number; G max is the maximum iteration number; v is a random number between 0 and 1.
Step 7: Update fitness values for all individuals in the population; Step 8: Calculate the individual with the best fitness value in the current generation; Step 9: If the termination condition "maximum iteration number" is satisfied, the calculation proceeds to the next step, otherwise, it returns to Step 4; Step 10: Output the best individual and its fitness value as the optimal solution.

Mind-Evolved Advanced Genetic Algorithm
Because of the above-mentioned shortcomings of GA algorithm, this paper proposes MAGA. This MAGA algorithm has a competitive relationship between multiple populations. Its main idea is to transform the single population of GA algorithm into several competing subpopulations based on the mind evolution method [45], and enhance the communication between groups, which can reflect the mind evolutionary characteristics to complicate and deepen a simple mind. Following the principle of survival of the fittest, the superior groups and the best individual are retained, while the substandard sub-populations are eliminated. New sub-populations are generated near the best individual to replace the eliminated ones, so that the individuals participating in GA algorithm can be continuously updated and the search is moving towards the optimal solution. Within the sub-population, the select operator, cross operator and mutation operator in the genetic process are improved. When a new superior individual appears in the sub-population, there is a greater opportunity to break the limitations of the current gene and the balance of certain comparative advantage, and then evolving a more adaptive individual to improve the competitiveness of this sub-population. This strategy not only preserves the individuals with excellent genes, but also ensures individual diversity in the population. In addition, the different search scope near the optimal solution can improve the local search ability in different regions of the global space.
This MAGA algorithm uses a rank-base fitness assignment to design the selection operator P g , which is only related to the ranking of individual fitness values, and is not directly affected by the specific value of fitness values. The individual with the best fitness has the highest probability to live to the next generation, and the second-best individual has the second-highest probability, and so on. The fitness value is only used to determine the ranking of individuals in the population, which can overcome the scale problem of proportional fitness. The select probability of individuals in the sorting space P g is calculated by Equation (7).
It should be ensured that the select probability of the last individual is the smallest, c ∈ (0.5, 1), and d g=1 P g = 1. g is the individual ranking number, d is the total number of individuals.
Cross and mutation operators determine the evolutionary ability of individual genes. This MAGA algorithm presents the improved cross and mutation operators, which are adjusted with the change of fitness value. When the individual fitness value is smaller than the average fitness value, the probability of this individual being changed is appropriately reduced to protect this better individual for the next generation. When the individual fitness value is larger than the average fitness value, the probability of this worse individual being crossed or mutated is increased to keep the diversity. Because the best individual may not be the global optimal solution, a certain minimum cross and mutation probability is also assigned to the best individual to increase the possibility of jumping out of a local optimum. The improved cross operator P c and mutation operator P m are updated by Equations (8) and (9). (9) where P c1 and P c2 are the initial set values of the cross operator, P m1 and P m2 are the initial set values of mutation operator, f avg is the average fitness value of the current generation population, f ' is the better fitness value of the two individuals to be crossed, f is the fitness value of the individual to be mutated. The basic flow of MAGA is shown in Figure 7. The MAGA algorithm steps are as follows: Step 1: The population is initialized to generate individuals carrying specific genetic information; Step 2: The first generation individuals are real-coded to keep them within the specified range; Step 3: The fitness function is determined and the fitness values of all individual is calculated; Step 4: Select some winning individuals and generate corresponding sub-populations near these individuals. The subsequent Step 5-14 are performed simultaneously in each sub-population; Step 5: Update the fitness values of all individuals in the sub-population. If a new winning individual appears, the calculation proceeds to the next step; otherwise, skip to Step 15; Step 6: Update the best individual in the current sub-population; Step 7: Update select operator P g according to Equation (7), and use the roulette method to perform the selection operation; Step 8: Update cross operator P c according to Equation (8), and perform the cross operation by Equation (4); Step 9: Update mutation operator P m according to Equation (9), and perform the mutation operation by Equations (5) and (6); Step 10: Update fitness values for all individuals in the sub-population, and calculate the individual with the best fitness value; Step 11: If the termination condition "maximum iteration number" is satisfied, the calculation proceeds to the next step, otherwise, it returns to Step 7; Step 12: Output the individual with the best fitness value. If this individual is better than the best individual in the current sub-population, update the best individual in the current sub-population; otherwise, skip this step; Step 13: Generate new sub-population near the best individual in the current sub-population, and update the fitness values of all individuals in the new sub-population; Step 14: If a new winning individual appears, the calculation returns to Step 6; otherwise, proceed to the next step; Step 15: Compare the fitness value of the best individuals of each sub-population, and update the optimal individual of the total population; Step 16: If the fitness value of the best individual in the sub-population meets the fitness value standard, the sub-population is retained; otherwise, all individuals in the sub-population are disbanded, and a new sub-population is randomly added near the optimal individual of the total population to keep the total number of sub-populations constant; Step 17: If the termination condition "maximum mind evolutionary number" is satisfied, the calculation proceeds to the next step, otherwise, all sub-populations return to Step 5; Step 18: Output the optimal individual of the total population and its fitness value as the optimal solution.

Verification of Improved Optimization Algorithms
In order to verify the effect of the advanced algorithm, four standard functions are used to test the global optimization capacity of the GA algorithm and MAGA algorithm. f 1 (x) is the Rosenbrock function, as shown in Equation (S1). It is a single-peak function, but it is pathological. f 2 (x), f 3 (x), and f 4 (x) are the Griewank function, the Rastrigrin function, and the Shaffer function, respectively, as shown in Equation (S2), Equation (S3), and Equation (S4). They are multi-peak and typical nonlinear multi-modal functions with a wide search space. All functions have many local minimum points but only one global minimum point. The setting parameters of different algorithms are shown in Table 6. These setting parameters are also applied in the optimization of neural network. The number of individuals in the population is set to 40. All the optimization processes are done on a 64-bit computer with a 3.80 GHz, AMD Ryzen 9 3900X 12-Core Processor CPU and 16 GB RAM.  Table S1 lists the optimal solutions of different optimization algorithms for different test functions under dimensions D = 2 and D = 20, respectively. Each value is the best solution obtained in the 100 independent runs of each algorithm. The theoretical minimum value of each test function is zero, which represents the global optimal solution. It can be seen from the comparison of the results that the optimal solution of the MAGA algorithm is obviously better than GA algorithm and closer to the theoretical minimum value in the definition domain. When D = 2, for f 1 (x), f 2 (x), f 3 (x), and f 4 (x), the values of the MAGA algorithm are very close to the minimum value of zero, where the value of f 1 (x) is slightly larger. In contrast, the values of GA algorithm are greater than those of the MAGA algorithm, and it is prematurely trapped in the local solution. As shown in Figures S1-S4, the position of the best individual can be directly observed, where the yellow dot with red circle f4(x) are the Griewank function, the Rastrigrin function, and the Shaffer function, respectively, as shown in Equation (S2), Equation (S3), and Equation (S4). They are multi-peak and typical nonlinear multi-modal functions with a wide search space. All functions have many local minimum points but only one global minimum point.
The setting parameters of different algorithms are shown in Table 6. These setting parameters are also applied in the optimization of neural network. The number of individuals in the population is set to 40. All the optimization processes are done on a 64-bit computer with a 3.80 GHz, AMD Ryzen 9 3900X 12-Core Processor CPU and 16 GB RAM.

PC Pm PC1 PC2 Pm1 Pm2
Mindevolution iteration  Table S1 lists the optimal solutions of different optimization algorithms for different test functions under dimensions D = 2 and D = 20, respectively. Each value is the best solution obtained in the 100 independent runs of each algorithm. The theoretical minimum value of each test function is zero, which represents the global optimal solution. It can be seen from the comparison of the results that the optimal solution of the MAGA algorithm is obviously better than GA algorithm and closer to the theoretical minimum value in the definition domain. When D = 2, for f1(x), f2(x), f3(x), and f4(x), the values of the MAGA algorithm are very close to the minimum value of zero, where the value of f1(x) is slightly larger. In contrast, the values of GA algorithm are greater than those of the MAGA algorithm, and it is prematurely trapped in the local solution. As shown in Figures S1-S4, the position of the best individual can be directly observed, where the yellow dot with red circle marks the global minimum value, and the red dot with green circle marks the optimal solution of the algorithm. It can be seen that the best individuals of the GA algorithm are all trapped in the local position, and its distance to the global minimum point is much larger than that of the MAGA algorithm. This shows that MAGA has a stronger ability to jump out of the local area and find the global optimal solution. When D = 20, the optimal values of GA and MAGA algorithms on all test functions are larger than those at D = 2. Howerver, comparing the two algorithms at D=20, the optimal solutions of the MAGA algorithm are still much smaller than the GA algorithm, where f2(x) has the best effect. Figures S5 and S6 show the iterative history of the two algorithms in different test functions, where MAGA convergence performance is better than the GA algorithm. It can be seen that the MAGA algorithm is suitable for functions of different dimensions with single or multi peaks, and its optimal performance is superior to the GA algorithm. Therefore, the improved optimization algorithm MAGA is reliable and can be applied to the optimization of the prediction model with the same setting parameters in Section 4. By optimizing the initial thresholds and weights of the neural network, the best initial values can be obtained. The fitness value is the test error of the neural network to obtain the smallest prediction error and improve the model accuracy.

Preprocessing of Sample Data
marks the global minimum value, and the red dot with green circle f4(x) are the Griewank function, the Rastrigrin function, and the Shaffer function, respectively, as shown in Equation (S2), Equation (S3), and Equation (S4). They are multi-peak and typical nonlinear multi-modal functions with a wide search space. All functions have many local minimum points but only one global minimum point.
The setting parameters of different algorithms are shown in Table 6. These setting parameters are also applied in the optimization of neural network. The number of individuals in the population is set to 40. All the optimization processes are done on a 64-bit computer with a 3.80 GHz, AMD Ryzen 9 3900X 12-Core Processor CPU and 16 GB RAM.  Table S1 lists the optimal solutions of different optimization algorithms for different test functions under dimensions D = 2 and D = 20, respectively. Each value is the best solution obtained in the 100 independent runs of each algorithm. The theoretical minimum value of each test function is zero, which represents the global optimal solution. It can be seen from the comparison of the results that the optimal solution of the MAGA algorithm is obviously better than GA algorithm and closer to the theoretical minimum value in the definition domain. When D = 2, for f1(x), f2(x), f3(x), and f4(x), the values of the MAGA algorithm are very close to the minimum value of zero, where the value of f1(x) is slightly larger. In contrast, the values of GA algorithm are greater than those of the MAGA algorithm, and it is prematurely trapped in the local solution. As shown in Figures S1-S4, the position of the best individual can be directly observed, where the yellow dot with red circle marks the global minimum value, and the red dot with green circle marks the optimal solution of the algorithm. It can be seen that the best individuals of the GA algorithm are all trapped in the local position, and its distance to the global minimum point is much larger than that of the MAGA algorithm. This shows that MAGA has a stronger ability to jump out of the local area and find the global optimal solution. When D = 20, the optimal values of GA and MAGA algorithms on all test functions are larger than those at D = 2. Howerver, comparing the two algorithms at D=20, the optimal solutions of the MAGA algorithm are still much smaller than the GA algorithm, where f2(x) has the best effect. Figures S5 and S6 show the iterative history of the two algorithms in different test functions, where MAGA convergence performance is better than the GA algorithm. It can be seen that the MAGA algorithm is suitable for functions of different dimensions with single or multi peaks, and its optimal performance is superior to the GA algorithm. Therefore, the improved optimization algorithm MAGA is reliable and can be applied to the optimization of the prediction model with the same setting parameters in Section 4. By optimizing the initial thresholds and weights of the neural network, the best initial values can be obtained. The fitness value is the test error of the neural network to obtain the smallest prediction error and improve the model accuracy.

Preprocessing of Sample Data
marks the optimal solution of the algorithm. It can be seen that the best individuals of the GA algorithm are all trapped in the local position, and its distance to the global minimum point is much larger than that of the MAGA algorithm. This shows that MAGA has a stronger ability to jump out of the local area and find the global optimal solution. When D = 20, the optimal values of GA and MAGA algorithms on all test functions are larger than those at D = 2. Howerver, comparing the two algorithms at D = 20, the optimal solutions of the MAGA algorithm are still much smaller than the GA algorithm, where f 2 (x) has the best effect. Figures  S5 and S6 show the iterative history of the two algorithms in different test functions, where MAGA convergence performance is better than the GA algorithm. It can be seen that the MAGA algorithm is suitable for functions of different dimensions with single or multi peaks, and its optimal performance is superior to the GA algorithm. Therefore, the improved optimization algorithm MAGA is reliable and can be applied to the optimization of the prediction model with the same setting parameters in Section 4. By optimizing the initial thresholds and weights of the neural network, the best initial values can be obtained. The fitness value is the test error of the neural network to obtain the smallest prediction error and improve the model accuracy.

Preprocessing of Sample Data
According to the analysis in Section 2.3, this paper selects three main factors, which are jet pressure, jet angle and impact time, as the input parameters of the prediction model, and erosion depth as the predicted output. The sample data are obtained by the orthogonal experimental design of L 25 (5 3 ). The orthogonal experimental design has the characteristic of "uniformity and decentralization, orderliness and comparable". The 125 data obtained by this hybrid orthogonal experimental design method are sufficient and reliable as the whole input data set for neural network, and satisfactory prediction results can be obtained [30,40,46,47]. Factors and levels of the tests are shown in Table 7. Table 8 records the results of the submerged jet erosion test on refractory bricks by an orthogonal table of L 25 (5 3 ). The standoff distance is fixed at 10 mm. BP neural network requires a large amount of effective sample data for training, but due to the constraints of the experimental conditions, it is difficult to obtain a sufficient number of experimental samples. Therefore, the virtual sample method [46] can be used to increase the number of samples data. In the experiment, due to the inevitable errors caused by the instrumentation, operation and measurement, there is a slight deviation between the input variable value and the experimental value, and this deviation is normal. Existing research [47] considers that all the input values within this deviation range correspond to the same output value, and usually the fluctuation range is ±∆ = 0.2%. Then, each experimental sample can generate 2 3 virtual samples, although not all these virtual samples are needed for training. With the orthogonal experimental design of L 4 (2 3 ), four virtual samples can be used to represent all these virtual samples corresponding to each experimental sample. As shown in Table 9, four virtual samples are generated based on the first experimental sample of test No.1. In this way, a total of 100 virtual samples can be generated. In order to avoid the difference in the magnitude of the data affecting the recognition accuracy of the network, the sample is normalized by the proportional compression method, as the Equation (10), so that the data mapping range is from −1 to 1.
where x is the original input data, x max and x min are the maximum and minimum values of the original data, respectively, and x' is the normalized input data.

Prediction Model of Optimized BP Neural Network
BP neural network [48] is a multi-layer feedforward network with a strong prediction ability. The topology of the single hidden layer is shown in Figure 8a. Figure   Variables x 1 , x 2 , . . . x n are the inputs from each neuron in the previous layer. ωj 1 , ωj 2 , . . . ωj n are the connection weights from each neuron in the previous layer to the j neuron. Sj is the input of the j neuron. θ j is the threshold. y i is the output of the j neuron. Activation function is Equation (11). The activation function of the hidden layer uses the tan-sigmoid function with Equation (12). The Levenberg-Marquardt algorithm is used to train the network. According to the principle of minimum root mean square error, the neuron number in a single hidden layer was set to 5 by the empirical formula as Equation (13) [40].
tan sig(S j ) = 2 where n i is the number of input layer neurons, n h is the number of hidden layer neurons, n o is the number of output layer neurons and a is an empirical constant from 0 to 10. Because the final prediction result is significantly affected by the initial thresholds and weights, GA and MAGA algorithms are used to optimize the initial values, respectively. According to the topology structure of the BP neural network 3-5-1, the optimized dimension is 26. The network test error is used as the fitness value to ensure that the prediction error is minimized. The initial thresholds and weights of the neural network are used to generate individuals in the population by genetic coding. When the optimization is completed, the best individual is decoded to obtain the corresponding optimal weights and thresholds, which are updated to initialize the BP neural network. Then the preprocessed data set is used to train, test, and regulate the network to obtain satisfactory model accuracy. Finally, erosion prediction is performed. The flow chart of the optimized BP neural network for prediction is shown in Figure 9.

Performance Evaluation
This paper uses the following evaluation functions as the indicators for the performance of the prediction model. The experimental value of erosion depth is denoted as y i , the predicted value of erosion depth is denoted asŷ i , and the number of prediction data is denoted as n.
(1) Errorsum is the sum of the deviation of all predicted points from the corresponding experimental value by Equation (14).
(2) APE max is the maximum absolute percentage error, which reflects the deviation of the worst point by Equation (15). APE reflects the deviation of each predicted point from the experimental value by Equation (16).
(3) MAPE is the mean absolute percentage error, which reflects the accuracy of the prediction results by Equation (17).
(4) RMSE is the root mean square error and is the most important indicator by Equation (18). This indicator reflects the accuracy of the prediction model. The smaller the value, the closer the predicted result is to the experimental value. It means that the prediction accuracy of this model is higher.
(5) R 2 is the coefficient of determination, which reflects the reliability of the prediction model, and is an important indicator of the goodness of fit between the predicted and experimental values. It is calculated by Equation (19). When 0.8 ≤ R 2 ≤ 1, it means a very strong correlation; When 0.6 ≤ R 2 ≤ 0.8, it means a strong correlation; 0.2 ≤ R 2 ≤ 0.4 means a weak correlation; 0 ≤ R 2 ≤ 0.2 means no correlation. The closer R 2 is to 1, the higher the prediction accuracy and goodness of fit.

Results and Discussion
All data for the prediction model include two different data sets, the input data set and the prediction data set, in order to ensure the generalization ability of the network. The input data set is a combination of the orthogonal test results and artificial virtual data in Section 4. The prediction data set comes from the orthogonal test results in Section 2. The Holdout-validation method is used to randomly divide the input data set into three data sets: training data set, validation data set, and test data set according to the proportion of 70%, 15%, and 15%, so as to separate the test data set and the training data set. In addition, the Early stopping strategy is adopted to prevent overfitting. The best accuracy model is selected by Holdout-validation with a certain number of times, and then used to predict the erosion depth. This method is effective and feasible [49]. Table 10 shows the prediction results of BP neural network optimized by different algorithms under different working conditions with the prediction data set. The deviation between the predicted value and the experimental value can be seen in Figure 10. In general, the prediction results of the MAGA-BP model are closer to the points of experimental values. In Figure 11, the fitting results between the prediction results and experimental values of the BP, GA-BP, and MAGA-BP models are shown, and their R 2 are 0.9848, 0.9904, and 0.9954, respectively. This shows that all three models have good reliability, and the fitting effect of MAGA-BP is the best.    Table 11 shows the five performance indicators of the three prediction models with Errorsum, RMSE, APE max , MAPE, and R 2 . ∆ BP-GABP represents the percentage change amount of GA-BP model relative to BP model; ∆ BP-MAGABP represents the percentage change amount of MAGA-BP model relative to BP model; ∆ GABP-MAGABP represents the percentage change amount of MAGA-BP model relative to GA-BP model. It can be seen that the sum of the errors of the BP, GA-BP and MAGA-BP models are 0.3531, 0.2862 and 0.1860, respectively. The Errorsum of the models optimized by GA and MAGA are reduced by 18.93% and 47.31% respectively compared to the BP model. The RMSE optimized by GA and MAGA algorithms reduces the error values from 0.0738 to 0.0530 and 0.0361, respectively, by -28.22% and -51.05%, which has the same order of magnitude in [30,46,47,50]. For APE max , the optimization effect of GA is not obvious and it still exceeds 10%, while the value of the model optimized by MAGA is 65.79% lower than it of the GA-BP model, and drops to 4%. The accuracy is significantly improved by MAGA. The MAPE of the BP, GA-BP, and MAGA-BP models are 5.08%, 4.57%, and 2.23%, respectively, of which the value of the MAGA-BP model is the smallest with a decrease of 56%. It can be seen that the results of all indicators of MAGA-BP are better than the other two models, so the optimization effect of the MAGA algorithm is obvious. In Figure 12, the APE max of the MAGA-BP prediction model is only about 4%, while the APE of other prediction points is less than 2%, and the MAPE is 2.23%, which is also less than 3%. According to existing studies [29,40,[51][52][53][54], this prediction error is acceptable. Therefore, the MAGA-BP prediction model is reliable. Compared with the BP prediction model, only the sixth prediction point has a slightly higher deviation, and the deviations at other points are lower than those of the BP prediction model. Compared with the GA-BP model, only the deviations of the fourth and fifth points are slightly larger. This individual deviation phenomenon is normal, which is mentioned by Wang et al. [30]. As a whole, the prediction model optimized by MAGA is more accurate.

Conclusions
This paper mainly proposes an optimized erosion prediction method to predict the erosion effect of the cavitation nozzle under different working conditions. On one hand, a prediction model based on erosion depth is established. On the other hand, an improved optimization algorithm is proposed. Then the improved optimization algorithm MAGA is applied to the prediction model to optimize the BP neural network and improve the prediction performance. Finally, the optimized erosion prediction method is intended to be applied in the formulation and control of the hull underwater cleaning solution, which is an important part of the ship hull underwater cleaning robot project. By this prediction method, the erosion effect of the nozzle under different working conditions can be obtained. Then, in actual operation, for different curved areas on the hull surface and different thicknesses of marine fouling, different schemes are formulated to timely modify the working condition parameters of the nozzle, so as to avoid the water-jet directly damaging the hull surface and anti-corrosion coating. In addition, feedback on the adjustment of the operating condition parameters to the cleaning robot can save energy consumption and optimize the control of the robot. For example, the change in jet pressure affects the regulation of the robot's adsorption and propulsion forces. The change in jet angle affects the control angle of the nozzle and the attitude of the underwater robot. The erosion time determines the robot's stay time and moving speed. With the establishment of a large database, the prediction method can also be used to select the appropriate nozzle. In addition, the optimized erosion prediction method of this paper can be summarized as follows: (1) Refractory brick is a simulated marine fouling, and its erosion depth is used as a test index, which can well reflect the erosion performance of submerged low-pressure water jets produced by angle nozzles under different working conditions. By predicting the erosion effect of water jets to select appropriate nozzles and corresponding optimal working conditions, the application of water jet technology in marine equipment can be promoted, and the efficiency of underwater operations can also be improved. The erosion depth can be controlled by adjusting the operating conditions, to ensure that while removing marine fouling, the hull surface and the anti-corrosion coating cannot be damaged. (2) Since the average error of the prediction result is less than 3%, the prediction model using the three factors of jet pressure, impact time, and jet angle as inputs and erosion depth as output is reliable. This also shows that it is feasible to analyze the primary and secondary factors of the model by using statistical analysis of range and variance based on orthogonal experimental design. The erosion depth increases with the increase of jet pressure and impact time, and decreases with the increase of standoff distance. With the increase of jet angle, the optimal angle appears around 60 • . The influence of different experimental factors on the index is jet pressure, impact time, jet angle and target distance in descending order. The jet pressure, impact time and jet angle are the main factors, and the standoff distance is the secondary factor in this experimental study. (3) In the complex test functions and BP neural network prediction, the performance index of the MAGA optimization algorithm has been greatly improved, and the optimization effect is obvious. It is feasible to apply the improved optimization algorithms to neural networks for improving prediction performance. The MAGA-BP neural network prediction model is also suitable for material research without an accurate constitutive model and can obtain higher prediction accuracy. Therefore, the optimized erosion prediction model plays a key role in ship hull underwater cleaning robot to achieve accurate cleaning control for different regions and different marine fouling.
Due to the limitation of experimental conditions, there are still many factors affecting erosion effect. Different nozzle structures, different working and environmental water quality, different underwater depths, and different variable value ranges may all affect the erosion ability. In addition, the evaluation index of erosion ability is not only the erosion depth, but also can be analyzed in combination with indicators such as the pit diameter and the erosion volume of the erosion damage to describe the damage caused by water jets more accurately. Therefore, the prediction model can be improved. In addition, different types of marine fouling can be simulated by changing the surface hardness of refractory bricks by burning them brittle before the test, or leaving them under-hardened when produced, etc. Finally, the optimized erosion prediction method can be applied in more practical fields.  Table S1: Comparison of optimization results on four test functions, Figure S1. The optimal individual location of f1(x) with different algorithms (D = 2), Figure S2: The optimal individual location of f2(x) with different algorithms (D = 2), Figure S3: The optimal individual location of f3(x) with different algorithms (D = 2), Figure S4