Optimization of a Screw Centrifugal Blood Pump Based on Random Forest and Multi-Objective Gray Wolf Optimization Algorithm

The centrifugal blood pump is a commonly used ventricular assist device. It can replace part of the heart function, pumping blood throughout the body in order to maintain normal function. However, the high shear stress caused by the impeller rotating at high speeds can lead to hemolysis and, as a consequence, to stroke and other syndromes. Therefore, reducing the hemolysis level while ensuring adequate pressure generation is key to the optimization of centrifugal blood pumps. In this study, a screw centrifugal blood pump was used as the research object. In addition, pressure generation and the hemolysis level were optimized simultaneously using a coupled algorithm composed of random forest (RF) and multi-objective gray wolf optimization (MOGWO). After verifying the prediction accuracy of the algorithm, three optimized models were selected and compared with the baseline model in terms of pressure cloud, 2D streamline, SSS distribution, HI distribution, and vortex distribution. Finally, via a comprehensive evaluation, the optimized model was selected as the final optimization design, in which the pressure generation increased by 24% and the hemolysis value decreased by 48%.


Introduction
Heart failure (HF) is a cardiovascular disease with high morbidity and mortality [1]. According to relevant studies [2,3], there are more than 64.3 million heart failure patients worldwide, and the incidence of heart failure is increasing with global population growth. Thus, it represents a serious issue for the global healthcare economy. For end-stage heart failure, heart transplantation is an effective treatment. However, the severe shortage of donors is an issue, and the long waiting times mean that a large number of patients miss the optimal time for heart transplantation. In addition, the increasing gap between the number of heart donors and heart failure patients is gradually increasing the mortality rate [4,5]. For patients with end-stage heart failure who are unable to undergo heart transplantation and are not suitable, a ventricular assist device (VAD) becomes the last option to significantly improve their chances of survival and their quality of life [6][7][8][9]. Centrifugal blood pumps are a common life support device [10] and are widely used in left heart assist devices and extracorporeal membrane pulmonary oxygenation (ECMO) [11,12]. These operations are characterized by the need for a high rotational speed in order to assist the left ventricle to return to a normal pumping volume [13]. Under high rotational speed, high shear stresses are generated near the impeller and at various gaps in the pump body [14]. This causes a large number of erythrocytes to rupture, which initiates the release of hemoglobin from the cells into the plasma, provoking hemolysis [15]. This also induces acquired vascular hemophilia and platelet activation [16,17], which can be life-threatening in severe cases. In addition, prolonged contact between the blood and foreign surfaces can cause platelets to activate and aggregate, forming clots [18]. Current methods of reducing blood damage are method, and a 3D model is built and numerically simulated to obtain a database consisting of impeller parameters and pressure generation, and HI values. Thereafter, the random forest prediction model is trained. The second part is multi-objective optimization, using the multi-objective gray wolf optimizer (MOGWO) to calculate the Pareto front. This involves selecting three different optimized models in the Pareto front for modeling and performing the numerical simulation. In the third part, the internal flow differences between the baseline model and the optimized model are analyzed, and the model with the best-combined situation is selected. Each step is described carefully in the following section.

Optimization Process
The overall optimization process for the blood pump is divided into three parts (Figure 1). The first step is the creation of a random forest prediction model. First, the parameters and boundary values controlling the impeller shape are determined. Then, 240 sets of impeller parameters are randomly obtained using the Latin hypercube sampling method, and a 3D model is built and numerically simulated to obtain a database consisting of impeller parameters and pressure generation, and HI values. Thereafter, the random forest prediction model is trained. The second part is multi-objective optimization, using the multi-objective gray wolf optimizer (MOGWO) to calculate the Pareto front. This involves selecting three different optimized models in the Pareto front for modeling and performing the numerical simulation. In the third part, the internal flow differences between the baseline model and the optimized model are analyzed, and the model with the best-combined situation is selected. Each step is described carefully in the following section.

Model Introduction and Parameterization
For this study, we used a self-developed screw centrifugal pump as the optimized object, because a screw centrifugal pump combines the advantages of volumetric pumps and centrifugal pumps [27]. The impeller was designed with a unique three-dimensional screw structure, which has superior nondestructive performance, a wider high-efficiency range, and better non-clogging characteristics than ordinary centrifugal pump impellers [28,29]. These advantages of the screw centrifugal pump reduce the amount of damage to the blood. A three-dimensional diagram and a cross section of the baseline model pump are shown in Figure 2.

Model Introduction and Parameterization
For this study, we used a self-developed screw centrifugal pump as the optimized object, because a screw centrifugal pump combines the advantages of volumetric pumps and centrifugal pumps [27]. The impeller was designed with a unique three-dimensional screw structure, which has superior nondestructive performance, a wider high-efficiency range, and better non-clogging characteristics than ordinary centrifugal pump impellers [28,29]. These advantages of the screw centrifugal pump reduce the amount of damage to the blood. A three-dimensional diagram and a cross section of the baseline model pump are shown in Figure 2. In the design process of centrifugal pumps, the axial projection diagram of the impeller is crucial. This is because a large number of parameters in the diagram can control the shape of the impeller [30,31] and changing any of these parameters will deform the impeller. This provides very convenient conditions for impeller parameterization. Figure  3 is the axial projection of the baseline model impeller. The shaded area in the figure represents the shape of the blade axis surface. Moreover, d is the impeller's outer diameter. In the design process of centrifugal pumps, the axial projection diagram of the impeller is crucial. This is because a large number of parameters in the diagram can control the shape of the impeller [30,31] and changing any of these parameters will deform the impeller. This provides very convenient conditions for impeller parameterization. Figure 3 is the axial projection of the baseline model impeller. The shaded area in the figure represents the shape of the blade axis surface. Moreover, d 2 is the impeller's outer diameter. In this study, this parameter was fixed at 30 mm. In addition, dS denotes the impeller inlet diameter, b 2 denotes the blade exit width, and the two circular curves represent the front cover line of the impeller and the rear cover line of the impeller, respectively. The straight line running from A to B is the position of the blade inlet side and it is defined by u 1 and u 2 , which denote the radial length from the intersection of the blade inlet edge line with the front and rear cover profiles to the impeller inlet, respectively. It can be seen that the position of the blade inlet edge is controlled by u 1 and u 2 together. There are certain impeller parameters that cannot be described visually in the axial projection diagram, such as the blade wrap angle (φ), the blade exit angle (β 2 ), and the axis ratio of the leading edge of the blade (A r ). These are also part of impeller parameterization as they have an impact on the impeller shape, hydraulic performance, and hemolysis performance of the blood pump [32,33].
In the design process of centrifugal pumps, the axial projection diagram of the impeller is crucial. This is because a large number of parameters in the diagram can control the shape of the impeller [30,31] and changing any of these parameters will deform the impeller. This provides very convenient conditions for impeller parameterization. Figure  3 is the axial projection of the baseline model impeller. The shaded area in the figure represents the shape of the blade axis surface. Moreover, 2 d is the impeller's outer diameter. In this study, this parameter was fixed at 30 mm. In addition, dS denotes the impeller inlet diameter, 2 b denotes the blade exit width, and the two circular curves represent the front cover line of the impeller and the rear cover line of the impeller, respectively. The straight line running from A to B is the position of the blade inlet side and it is defined by 1 u and 2 u , which denote the radial length from the intersection of the blade inlet edge line with the front and rear cover profiles to the impeller inlet, respectively. It can be seen that the position of the blade inlet edge is controlled by 1 u and 2 u together. There are certain impeller parameters that cannot be described visually in the axial projection diagram, such as the blade wrap angle (φ), the blade exit angle ( 2 β ), and the axis ratio of the leading edge of the blade ( r A ). These are also part of impeller parameterization as they have an impact on the impeller shape, hydraulic performance, and hemolysis performance of the blood pump [32,33].  It is worth mentioning that, in this study, we used an elliptical leading-edge design scheme [34]. The top view of the blade's leading-edge shape is shown in Figure 4. The sharpness of the ellipse is expressed by Equation (1): where w L and w S represent the long and short axes of the ellipse, respectively. The blade thickness increases linearly from 1 mm at the leading edge (w S ) to 2 mm at the trailing edge, and the axial thickness remains constant. The inlet and outlet diameters of the volute are 10 mm. where L w and S w represent the long and short axes of the ellipse, respectively. The blade thickness increases linearly from 1 mm at the leading edge ( S w ) to 2 mm at the trailing edge, and the axial thickness remains constant. The inlet and outlet diameters of the volute are 10 mm. Table 1 summarizes all the parameter values of the baseline model and their variation ranges, where L x and U x represent the lower and upper bounds of the impeller parameters, respectively.  When building proxy models for expensive functions or experimental data, the preexperimental design is particularly important for the selection of the initial sampling points. Thus, how to achieve good space-filling effects with as few samples as possible is a hot research topic. In computer experiments, classical sampling methods include Monte Carlo (MCS), Latin hypercube sampling (LHS), Sobol sequence, etc. [35]. Among them, Latin hypercube sampling is popular. It is a random sampling method that was improved by Mckay et al. [36,37] based on the Monte Carlo sampling method. They enhanced a core sampling strategy with stratified sampling, which ensures full coverage of the multivariate distribution and accuracy, while significantly reducing the sample size. In this study, 240 sets of impeller parameters were randomly generated using the Latin hypercube sampling method in the range of impeller parameter variation. Thereafter, all of them were modeled and numerically simulated to obtain a database consisting of impeller parameters and head and HI values.  When building proxy models for expensive functions or experimental data, the preexperimental design is particularly important for the selection of the initial sampling points. Thus, how to achieve good space-filling effects with as few samples as possible is a hot research topic. In computer experiments, classical sampling methods include Monte Carlo (MCS), Latin hypercube sampling (LHS), Sobol sequence, etc. [35]. Among them, Latin hypercube sampling is popular. It is a random sampling method that was improved by Mckay et al. [36,37] based on the Monte Carlo sampling method. They enhanced a core sampling strategy with stratified sampling, which ensures full coverage of the multivariate distribution and accuracy, while significantly reducing the sample size. In this study, 240 sets of impeller parameters were randomly generated using the Latin hypercube sampling method in the range of impeller parameter variation. Thereafter, all of them were modeled and numerically simulated to obtain a database consisting of impeller parameters and head and HI values.

Numerical Simulation and Boundary Conditions
The model for the area calculation used for the numerical simulation is shown in Figure 5. The whole stream field area was divided into three areas: the inlet domain, the impeller domain, and the volute domain. The structure of the inlet domain and volute domain is relatively simple and can be drawn manually using the 3D modeling software NX12.0. The impeller domain was automatically generated using CFTurbo2021 R2.2 by inputting parameters.
In this study, Fluent Meshing, a fast and high-quality pre-processing meshing software, was used to discretize the computational domain. High-adaptability polyhedral mesh was used inside the computational domain. Six-layer structured mesh was used for the walls. In addition, the mesh of the blade wall was encrypted to ensure a more accurate calculation of flow details. Finally, the orthogonal qualities of the meshes were all greater than 0.3, and y + near the wall was less than 5, conforming to SST K − ω turbulence modeling standards. To reduce computational costs, the mesh independence of the baseline model was verified at a speed of 3900 rpm and a flow rate of 5 L/min. The mesh-independence verification carried out is presented in this paper. Table 2 shows that when the number of meshes was more than 3,825,107, the pressure generation value no longer changed. To reduce the simulation time, the number of meshes in option 3 was chosen.
where H represents the pressure generation, out p denotes the outlet pressure (Pa), in p denotes the inlet pressure (Pa), and ρ is the blood density.
The Lagrangian method and the power-law model based on scalar shear force and red blood cell exposure time proposed by Giersiepen et al. [48] were used to calculate the intra-pump hemolysis values. The mathematical formulation of the model is as follows: where Hb is total hemoglobin, dH b is total free hemoglobin, C α β 、 、 are the empirical constants ( In conjunction with the SSS calculation proposed by Bludszuweit [49], The final expression of SSS takes the following form:   The finite element analysis software Fluent (Ansys, Canonsburg, PA, USA) was used as the solver for the numerical simulations in this paper. Under high shear stress, blood can usually be considered an incompressible Newtonian fluid with a viscosity of 0.0035 kg/m.s and a density of 1050 kg/m 3 [38][39][40]. Flow in the pump is controlled by the Reynolds averaged Navistokes (RANS) equations [41]. In this study, steady-state simulations were performed using the SST K − ω turbulence model with dual controlling equations [42,43]. The control equation of the SST K − ω turbulence model is shown in Equation (2): where ω is the unknown specific dissipation rate; µ t is the eddy viscosity; F B is the mixing function; β * = 0.09; and σ ω out = 1.168. According to the flow rate of 5 L/min and the inlet cross-sectional area, the blood pump inlet was set to a velocity inlet boundary condition of 1.061 m/s, the blood pump outlet was set to the outflow outlet condition, the rotational speed was controlled at 3900 rpm, the impeller walls were set as rotating walls, and the remaining walls were assumed to be slip-free. The computational domain contains one rotational domain and two stationary domains. Based on the numerical simulation study of the blood pump [17,[44][45][46][47], the multiple reference coordinate system (MRF) was selected, in which the impeller domain was the rotational coordinate system, and the inlet domain and the worm-casing domain were the stationary coordinate systems. The surfaces overlapping between the three computational domains were set as interfaces. The semi-implicit SIMPLEC algorithm was used to solve the pressure-velocity-coupled equations, and the convergence criterion was set to 10 −5 . After spatial discretization, all equations were set as second-order windward equations.

Calculation of Pressure Generation, Scalar Shear Stress, and HI Value
The calculation of the pressure generation in this study is as follows: where H represents the pressure generation, p out denotes the outlet pressure (Pa), p in denotes the inlet pressure (Pa), and ρ is the blood density. The Lagrangian method and the power-law model based on scalar shear force and red blood cell exposure time proposed by Giersiepen et al. [48] were used to calculate the intra-pump hemolysis values. The mathematical formulation of the model is as follows: where Hb is total hemoglobin, dHb is total free hemoglobin, C, α, β are the empirical constants (C = 3.62 × 10 −7 , α = 0.785, β = 2.418), t is the erythrocyte exposure time, and τ is the scalar shear stress (later referred to as SSS). SSS is the sum of the viscous shear stress (σ ij ) and Reynolds shear stress (ρu i u j ), i.e., In conjunction with the SSS calculation proposed by Bludszuweit [49], The final expression of SSS takes the following form:

Random Forest (RF)
Random forest (RF) was proposed by Leo Breiman [50] in 2001. As shown in Figure 6, it is a prediction algorithm that combines bagging integrated learning theory with a random subspace. Its core lies in N CART (classification and regression trees) composed of g(D, θ n ). The final prediction result can be obtained by averaging the regression results of all regression trees, which means the random forest algorithm is more resistant to overfitting, has a lower error rate [51], and is more accurate when making predictions on small samples. Random forest (RF) was proposed by Leo Breiman [50] in 2001. As shown in Figure  6, it is a prediction algorithm that combines bagging integrated learning theory with a random subspace. Its core lies in N CART (classification and regression trees) composed of ( , ) n g D θ . The final prediction result can be obtained by averaging the regression results of all regression trees, which means the random forest algorithm is more resistant to overfitting, has a lower error rate [51], and is more accurate when making predictions on small samples.
In this study, we randomly selected 200 sets of data from the database mentioned above and used them as the training set. In addition, we trained and built two independent random forest proxy models to predict pressure generation and HI. The input values of the proxy models were the impeller parameters, the output values were HI and pressure generation, and the rest of the data were used as the prediction set to evaluate the prediction accuracy of the models. To quantitatively evaluate the prediction accuracy of the random forest prediction models, in this study, mean absolute percentage error (MAPE) [52] and mean square error (MSE) [53] were used as model accuracy evaluation metrics for the HI prediction set and the pressure generation prediction set. MAPE and MSE are calculated as follows: In this study, we randomly selected 200 sets of data from the database mentioned above and used them as the training set. In addition, we trained and built two independent random forest proxy models to predict pressure generation and HI. The input values of the proxy models were the impeller parameters, the output values were HI and pressure generation, and the rest of the data were used as the prediction set to evaluate the prediction accuracy of the models.
To quantitatively evaluate the prediction accuracy of the random forest prediction models, in this study, mean absolute percentage error (MAPE) [52] and mean square error (MSE) [53] were used as model accuracy evaluation metrics for the HI prediction set and the pressure generation prediction set. MAPE and MSE are calculated as follows: where y i andŷ i are the true and predicted values of the i-th sample, y is the average of the actual values of the samples, and n is the capacity of the corresponding sample. Figure 7 shows the predicted and true values of the random forest algorithm, where the blue curve represents the predicted value and the red curve represents the true value. The curve trends of the prediction set and the training set are consistent with a high degree of similarity, with some errors existing at extreme values only. The quantitative model accuracy evaluation is shown in Table 3.    In Table 3, the MSE and MAPE of both the pressure generation training set and the HI training set were maintained below 0.5, which demonstrated that the random forest model has a high prediction accuracy and can accurately regress the mapping relationship between the design parameters and the pressure generation and HI levels.

Multi-Objective Gray Wolf Optimization Algorithm (MOGWO)
The gray wolf optimization algorithm (GWO) is a pack intelligence optimization algorithm designed by Mirjalili [54]. It was inspired by the social stratification characteristics and hunting and trapping behavior of wolves. It has the advantages of strong convergence performance, a simple structure, and easy implementation. There exist convergence factors with self-adaptive adjustment and information feedback mechanisms. Thus, a balance between local optimization and global optimization is achieved. For this reason, it has good problem-solving accuracy. Gray wolves are divided into four social classes: α, β, δ, ω. After sorting according to the adaptation function, α wolves represent the optimal solutions. β wolves and δ wolves represent the suboptimal solutions, with their role to assist α wolves in making decisions. The remaining candidate solutions are defined as ω wolves. Classes of α wolves, β wolves, and δ wolves command the hunting behavior, and ω wolves follow the above-mentioned higher-level wolves in hunting. Since the position of the prey (the optimal solution) in the solution space is not known, it is assumed that the positions of wolf α, wolf, and wolf δ are closest to the prey.
As shown in Figure 8 [54], the hunting process of the gray wolf surrounding prey is represented by the following mathematical model: where → C and → A are coefficient vectors determined by random values, → X p (t) and → X(t) represent the position vector of the prey and the position vector of the gray wolf at iteration up to the t-th generation, respectively [54]. Since the location of the optimal solution in the solution space is not known, it is assumed that the locations of wolf α, wolf β, and wolf δ are closest to the optimal solution. After recording the positions of these three wolves, the ω wolf is ordered to approach the α wolf, and the β wolf is ordered to approach the δ wolf. During each iteration, the positions of the α wolf, the β wolf, and the δ wolf are updated via the formulae shown in Equation (12) [54]: Averaging the positions of the α wolf, the β wolf, the δ wolf, and the result obtained is regarded as the final position after this iteration is updated. When | → A| < 1 in Equation (12), the next generation of gray wolves can be located anywhere near the prey. The constant repetition of this behavior is to hunt the prey. It is worth mentioning that in order to prevent falling into local optimal solutions during the optimization process, the gray wolf chooses to move away from the prey when | → A| > 1.   Figure 10. The results show that the predicted values of the three optimized models are in good agreement with the simulation values, so the prediction using the random forest algorithm can be considered a reliable method.  To enable the gray wolf optimization algorithm to perform multi-objective optimization, Mirjalili provided two components [55]. One is a storage component responsible for storing non-dominated Pareto optimal solutions, which implements the storage function of several optimal solutions. The other is a leader-selection strategy component, which can obtain dominance relations with the help of the Pareto global optimum concept but cannot compare solutions that are not dominated by each other. This component selects a new leader in the uncrowded region of non-dominated solutions according to the roulette wheel method, marked as α, β, δ, which is then archived.

Results
In this study, two objective functions were chosen as constraints and the expressions are as shown in Equation (14): where P and HI represent the random forest HI prediction model and the random forest indentation pressure generation prediction model, respectively, x i (i = 1, 2, . . . , 6) represents the i-th design parameter, and x L and x U are the upper and lower bounds of each design parameter, respectively. All codes in this paper were written in the MATLAB environment, the number of populations in each iteration was set to 100, and the maximum number of external archives was chosen to be 100, for a total of 400 iterations, running on a PC with Inlet(R) Core(TM) i7-12700 k CPU 5.00 GHz 32 GB RAM.  Figure 10. The results show that the predicted values of the three optimized models are in good agreement with the simulation values, so the prediction using the random forest algorithm can be considered a reliable method. function optimal solution point (the point of the highest-pressure generation) and point B is the second objective function optimal solution point (the point of the lowest HI value). Point C is the intermediate point near point B. To verify the prediction accuracy of the random forest and MOGWO models, a numerical simulation was performed for models A, B, and C. The simulated and predicted values of each objective function for points A, B, and C are shown in Figure 10. The results show that the predicted values of the three optimized models are in good agreement with the simulation values, so the prediction using the random forest algorithm can be considered a reliable method.

Impeller Parameter Analysis
The impeller parameters of the four models before and after optimization are shown in Table 4. The pressure generation at point A increased by 41% as compared with the baseline model, but the HI increased by 105%. The axis ratio ( r A ) did not change significantly as compared with the benchmark model, but the wrap angle decreased by 27 degrees and the blade outlet width increased by 1.1 mm. In addition, the A-point model exhibited the largest blade outlet width and exit blade angle. Therefore, the area of blood worked by the blades was larger than that of all the other model pumps. In terms of the blade inlet edge position, the A-point model 2 u value was greatly reduced, which means the blade inlet edge along the rear cover plate profile produced a large forward extension, the blade shape was distorted, and the hub side shrank significantly toward the center. B-point model pressure generation was improved by 19% as compared with the base model, and HI was reduced by 50%. At this time, the B-point model wrap angle was reduced by 60° and the 1 u and 2 u values were lower than the base model. This allowed

Impeller Parameter Analysis
The impeller parameters of the four models before and after optimization are shown in Table 4. The pressure generation at point A increased by 41% as compared with the baseline model, but the HI increased by 105%. The axis ratio (A r ) did not change significantly as compared with the benchmark model, but the wrap angle decreased by 27 degrees and the blade outlet width increased by 1.1 mm. In addition, the A-point model exhibited the largest blade outlet width and exit blade angle. Therefore, the area of blood worked by the blades was larger than that of all the other model pumps. In terms of the blade inlet edge position, the A-point model u 2 value was greatly reduced, which means the blade inlet edge along the rear cover plate profile produced a large forward extension, the blade shape was distorted, and the hub side shrank significantly toward the center. B-point model pressure generation was improved by 19% as compared with the base model, and HI was reduced by 50%. At this time, the B-point model wrap angle was reduced by 60 • and the u 1 and u 2 values were lower than the base model. This allowed for the blade inlet side to be moved forward and positioned more axially. The blade exit width of the B-point model increased by 0.3 mm, which located the area where the blade worked on the blood and pressure generation between the baseline model pump and the A-point model pump. In terms of the wrap angle, the B-point model exhibited the most significant change, decreasing by 60 • as compared with the baseline model, while the shaft ratio and blade exit did not exhibit a significant change.
The pressure generation of the C-point model was improved by 24%, and the HI value was reduced by 48%. The parameters of the C-point model and the B-point model were similar, with the main differences being related to the blade outlet width. Specifically, the blade exit width of the C-point model was slightly larger than that of the B-point model. The larger work area and surface area resulted in a simultaneous increase in pressure generation and the HI value; however, as can be seen from Figure 10, the increase in the HI value for the C-point model was almost negligible, while there was a significant increase in the pressure generation value. Figure 11 depicts the pressure distribution clouds of the baseline model and the optimized model in the XY direction. All models had a low-pressure zone at the impeller inlet, with the impeller diameter gradually increasing the pressure and finally reaching its maximum at the trailing edge of the blade and the volute. The optimized pressure values at the outlet of all models were improved to different degrees, with the largest pressure of 25,028 Pa being observed at the outlet of the A-point model, although there was a region of uneven pressure distribution at the volute. The pressures at the outlet of the B-point model and C-point model were 19,535 Pa and 19,679 Pa, respectively. Both exhibited steady pressure growth, with no obvious excessive pressure gradient area, but the outlet pressure value of the C-point model was slightly higher than that of the B-point model. Therefore, the C-point model can provide higher pressure generation than the B-point model under the same operating conditions. sure of 25,028 Pa being observed at the outlet of the A-point model, although there was a region of uneven pressure distribution at the volute. The pressures at the outlet of the Bpoint model and C-point model were 19,535 Pa and 19,679 Pa, respectively. Both exhibited steady pressure growth, with no obvious excessive pressure gradient area, but the outlet pressure value of the C-point model was slightly higher than that of the B-point model. Therefore, the C-point model can provide higher pressure generation than the B-point model under the same operating conditions.

Stream Field Analysis
In this study, the three-dimensional streamlines were discretized into two-dimensional streamlines on the monitoring surface to analyze the stream velocity and stream pattern in the impeller calculation domain. Figure 12 illustrates the selection of the monitoring surfaces. A total of eight monitoring surfaces were selected from , , , , [1,4] i

i i i a b c h i ⋅⋅⋅⋅ ∈
(the subscripts i from 1 to 4 represent the base model, the A-point model, the B-point model, and the C-point model, respectively) from the axis upward, and the angle between adjacent monitoring surfaces was 45 degrees. Figure 13 shows the stream field comparison of the same monitoring surfaces. In terms of the stream pattern, all models exhibited a stable stream near the impeller inlet with a gentle velocity. This

Stream Field Analysis
In this study, the three-dimensional streamlines were discretized into two-dimensional streamlines on the monitoring surface to analyze the stream velocity and stream pattern in the impeller calculation domain. Figure 12 illustrates the selection of the monitoring surfaces. A total of eight monitoring surfaces were selected from a i , b i , c i , · · ··, h i i ∈ [1,4] (the subscripts i from 1 to 4 represent the base model, the A-point model, the B-point model, and the C-point model, respectively) from the axis upward, and the angle between adjacent monitoring surfaces was 45 degrees. Figure 13 shows the stream field comparison of the same monitoring surfaces. In terms of the stream pattern, all models exhibited a stable stream near the impeller inlet with a gentle velocity. This stream pattern changed rapidly to a turbulent pattern when the blood flowed through the impeller. The overall stream chaos was the highest in the A-point model, and an obvious low-speed vortex area could be observed in the monitoring surface of Figure 13 b 2 , c 2 , d 2 . As a result of the large contraction of the rear cover side of the blade in the A-point model, the difference between the front cover side and the rear cover side of the blade in the axial control direction of blood was obvious. Thus, we observed a low-speed stream separation area touching the suction side of the blade in all monitoring surfaces of the A-point model. The overall stream pattern in the B-point and C-point models was stable, and there was no obvious stream separation in the stream channel. During the rotation of monitoring surface a to monitoring surface e, the streamline of the impeller outlet part was gradually distorted, and a local vortex formed at monitoring surfaces e and f (Figure 13 e 3 , e 4 , f 3 , f 4 ). As for the rate of stream velocity change, in the stream channel between adjacent blades of the baseline model (Figure 13 d 1 , e 1 , f 1 , g 1 , h 1 ), there was a large velocity gradient region, which resulted from the narrowing of the stream channel. This, in turn, was caused by the large blade wrap angle of the baseline model and the crowding-out effect of the blades on the fluid in the stream channel, which increased the stream velocity. The reduction in the wrap angle of the optimized model widened the stream channel between the blades and reduced the bloodstream velocity in the stream channel. As a result, no high-velocity zone was generated in the stream channel of any model after optimization, and only a low-velocity return zone was observed in the axial direction on the b 2 , d 2 monitoring surface of the A-point model. The velocity variation in the impeller exit region of the A-point model was the most dramatic, and there were obvious high-velocity zones from monitoring surface e to monitoring surface h. The coexistence of the maximum and minimum velocity could even be observed at the monitoring surface of f 2 , g 2 . The overall stream velocity was stable from monitoring surface a to monitoring surface c in the B-point model and C-point model (Figure 13 a 3 , a 4 , b 3 , b 4 , c 3 , c 4 ), and only a high-speed vortex area appeared on surface f. Therefore, the blood exhibited a good stream within the B-point model and C-point model.
The stream condition in the C-point model was similar to that of the B-point model, but the larger blade area in the C-point model caused a slight increase in the HI value. 2 2 exit region of the A-point model was the most dramatic, and there were obvious highvelocity zones from monitoring surface e to monitoring surface h. The coexistence of the maximum and minimum velocity could even be observed at the monitoring surface of 2 f , 2 g . The overall stream velocity was stable from monitoring surface a to monitoring surface c in the B-point model and C-point model (Figure 13 3  The blade-to-blade sections were intercepted in the radial direction for spans of 0 0.5, and 0.9, respectively. There were more vortices in the A-point model as compar with the other models. In the cross section with a span = 0.2 (Figure 14b), there wa significant backflow in the stream channels of both blade-leading edges for the A-po model. As the impeller diameter increased, the stream velocity in the middle section The blade-to-blade sections were intercepted in the radial direction for spans of 0.2, 0.5, and 0.9, respectively. There were more vortices in the A-point model as compared with the other models. In the cross section with a span = 0.2 (Figure 14b), there was a significant backflow in the stream channels of both blade-leading edges for the A-point model. As the impeller diameter increased, the stream velocity in the middle section of the pressure side of the blade decreased, and the streamline spread toward the impeller outlet. The velocity vector on the suction side of the trailing edge of the blade suddenly changed and a stagnation zone was formed. In addition, there were two high stream velocity zones at the gap between the impeller and the volute, which were the same as the monitoring surface e-h, so this was not repeated. In the span = 0.5 cross section (Figure 15b) and the span = 0.9 cross section (Figure 16b), it can be clearly seen that a relatively large low-speed backflow area covered almost all of the pressure surface of the blade, indicating that the most serious area of backflow occurred on the front cover side of the blade working surface. From the three cross sections of the B-point model (Figures 14c, 15c and 16c), it can be seen that the streamlines within the impeller of the B-point model were uniformly distributed and close to the blade profile. Moreover, they fit closely to the blade and there was no stream separation region in the blade suction surface. In addition, the stream velocity in the stream channel was controlled in the range of 1 m/s~4 m/s, and the stable stream pattern minimized the destruction of red blood cells in the stream channel. Three radial sections of the C-point model exhibited similar stream conditions as those of the B-point model, but in the span = 0.5 and 0.9 sections (Figures 15d and 16d), the left blade suction surface had fewer low-velocity areas and reflux areas (Figures 15c and 16c), which meant the C-point model had the lowest risk of thrombosis.   (Figures 15d and 16d), the l blade suction surface had fewer low-velocity areas and reflux areas (Figures 15c and 16 which meant the C-point model had the lowest risk of thrombosis.

SSS and HI Index Analysis
In the HI index formula, the SSS value dominates, so the distribution of SSS correlat with the HI value distribution of the blood pump. Figure 17 compares the SSS distributio

SSS and HI Index Analysis
In the HI index formula, the SSS value dominates, so the distribution of SSS correlates with the HI value distribution of the blood pump. Figure 17 compares the SSS distribution clouds on the blade surfaces of the four models and the HI index distribution clouds. In all models, the SSS on the pressure surface of the blade was generally higher than that on the suction surface. The majority of the HI values of the baseline model blade (Figure 17a) were between 1.894 × 10 −3 and 2.502 × 10 −3 , and the high SSS region (100 Pa-140 Pa) was mainly concentrated in the flow path near the leading edge of the blade. Thus, there was a local high HI region at the same location on the HI distribution map. The high SSS region and the high HI region of the A-point model were the most widely distributed. The larger blade area and the chaotic flow pattern on the pressure side of the blade (Figure 15b) led to a sudden increase in SSS at the leading and trailing edges of the front cover side of the blade (Figure 17b), so a large high SSS region (>100 Pa) and a large HI region covered the entire trailing edge of the blade. The B-point model (Figure 17c) exhibited the smallest high SSS region, except for the high shear region at the leading edge of the blade, with SSS in the rest of the parts remaining below 100 Pa. From the HI distribution, it can be seen that the middle section of the blade caused the least damage to the blood. This was essentially below 1.742 × 10 −3 . In the leading and trailing edges of the blade, the HI value rose to 2.654 × 10 −3 . The C-point model maintained a similar HI distribution and SSS distribution as the B-point model, but the degree of blood damage at the impeller inlet was slightly higher than that of the B-point model due to the difference in the shape of the leading edge of the blade. Finally, in this study, the vortex situation in all model fluid domains was analyzed using a vortex identification method based on the Q-criterion. This is, in turn, based on the second matrix invariant of the velocity gradient tensor [56] and Equation (15): Finally, in this study, the vortex situation in all model fluid domains was analyzed using a vortex identification method based on the Q-criterion. This is, in turn, based on the second matrix invariant of the velocity gradient tensor [56] and Equation (15): where F denotes the Frobenius parametrization of the matrix, and A and B represent the deformation and rotation of a point in the flow field, respectively. In this study, Q = 0.01 was used to identify the vortex in the impeller domain, as shown in Figure 18.

Discussion
Related studies have confirmed that hemolysis is an important cause of blood damage in blood pumps [57][58][59]. In the published literature, most scholars reduce the HI value in blood pumps and thrombus by manually adjusting the impeller and volute parameters [60][61][62]. However, trial-and-error analysis and numerical simulations can be time-consuming and ignore the co-influence of the impeller parameters on the hydraulic performance and HI performance, which may result in bypassing the optimal combined solution. The multi-objective optimization algorithm can obtain a set of optimal solutions in a short time and select the final optimization result independently according to the requirements, which greatly improves optimization efficiency and speed. In this study, we developed an optimization process for a screw centrifugal blood pump using a multi-algorithm machine-learning approach. The whole process is divided into three parts: random forest algorithm prediction, MOGWO optimization, and internal flow field analysis.
During the random forest algorithm prediction stage, 240 sets of impeller parameters were randomly generated using the Latin hypercube sampling method. Then, these were modeled and numerically simulated to obtain a dataset consisting of impeller parameters, pressure generation, and HI values. Thereafter, two independent random forest prediction models were established with impeller parameters as input values, and pressure generation and HI values as output values. These were used as the objective functions in the multi-objective optimization algorithm. During the multi-objective optimization stage, the multi-objective gray wolf optimization algorithm was selected as the optimizer for this study. The model that exhibited the highest-pressure generation, the model that exhibited the lowest HI value, and the middle point model were selected as the optimized models in the Pareto front with iterations up to 400 generations.
Finally, the hydraulic performance and HI performance of all blood pumps were compared using pressure cloud, velocity flow line, SSS distribution, HI distribution, and vortex distribution plots. In terms of hydraulic performance, all optimized models exhibited different degrees of pressure generation increase. The A-point model had the largest pressure generation increase because the maximum vane outlet angle and vane outlet width provided the largest effective work area. The B-point model and C-point model demonstrated a lesser pressure generation increase than the A-point model, but the C-

Discussion
Related studies have confirmed that hemolysis is an important cause of blood damage in blood pumps [57][58][59]. In the published literature, most scholars reduce the HI value in blood pumps and thrombus by manually adjusting the impeller and volute parameters [60][61][62]. However, trial-and-error analysis and numerical simulations can be time-consuming and ignore the co-influence of the impeller parameters on the hydraulic performance and HI performance, which may result in bypassing the optimal combined solution. The multi-objective optimization algorithm can obtain a set of optimal solutions in a short time and select the final optimization result independently according to the requirements, which greatly improves optimization efficiency and speed. In this study, we developed an optimization process for a screw centrifugal blood pump using a multi-algorithm machine-learning approach. The whole process is divided into three parts: random forest algorithm prediction, MOGWO optimization, and internal flow field analysis.
During the random forest algorithm prediction stage, 240 sets of impeller parameters were randomly generated using the Latin hypercube sampling method. Then, these were modeled and numerically simulated to obtain a dataset consisting of impeller parameters, pressure generation, and HI values. Thereafter, two independent random forest prediction models were established with impeller parameters as input values, and pressure generation and HI values as output values. These were used as the objective functions in the multiobjective optimization algorithm. During the multi-objective optimization stage, the multiobjective gray wolf optimization algorithm was selected as the optimizer for this study. The model that exhibited the highest-pressure generation, the model that exhibited the lowest HI value, and the middle point model were selected as the optimized models in the Pareto front with iterations up to 400 generations.
Finally, the hydraulic performance and HI performance of all blood pumps were compared using pressure cloud, velocity flow line, SSS distribution, HI distribution, and vortex distribution plots. In terms of hydraulic performance, all optimized models exhibited different degrees of pressure generation increase. The A-point model had the largest pressure generation increase because the maximum vane outlet angle and vane outlet width provided the largest effective work area. The B-point model and C-point model demonstrated a lesser pressure generation increase than the A-point model, but the C-point model exhibited a 5% higher pump pressure generation than the B-point model. Thus, the C-point model had a higher hydraulic efficiency than the B-point model. The hydraulic efficiency of the C-point model was higher than that of the B-point model. With regards to HI performance, a power law equation based on SSS and exposure time was used to assess the blood compatibility of the blood pump. In the baseline model, the overall flow state of the blood was smooth, but a high shear zone caused by a high-flow velocity existed in the flow channel near the leading edge of the impeller, resulting in a high HI value at this location. The inlet side of the blade in the A-point model extended significantly forward along the rear cover plate profile, the hub side of the blade shrank toward the center, and the irregular impeller shape led to a large flow separation and backflow near the working surface of the blade, which elevated the thrombus formation risk. In the radial flow diagram, it can be observed that the flow separation area with the large velocity gradient was touching the middle and trailing edge of the impeller, so the blade of the A-point model had the largest high SSS area and high HI area, and the blood damage degree of the A-point model was the most serious. The B-point model blade on the inlet side, along the front cover plate profile forward extension, caused the blade to be biased to the axial placement. Each height of the blade was similar in terms of the direction of action to the blood, so the flow pattern was stable and close to the blade surface. In addition, the spacious flow channel caused a small change in the flow velocity and no flow separation in the flow field. The majority of the blade SSS was maintained below 100 Pa. In the HI distribution, the HI values of all parts of the blade were below 2.654 × 10 −3 . The SSS and HI distributions of the C-point model were similar to those of the B-point model, but the HI values were 2% higher than those of the B-point model. This may indicate that the increase in the HI value for the C-point model originated from the increase in the impeller surface area rather than from the difference in the impeller shape, and if the pressure generation of the C-point model was reduced to the same value as that of the B-point model by reducing the speed, its HI value would decrease accordingly. In summary, the C-point model exhibited the largest increase in pressure generation while reducing the HI value under the same operating conditions. Thus, it is the best solution among all the optimized models.
There are some limitations to this study. First, because the impeller diameter and the base circle diameter of the pressurized water chamber were fixed during the design process, the width of the gap between the blade and the volute sidewall was not included as a design variable. In fact, the width of the gap also had an impact on the HI level [22]. In addition, the friction between the rotor and the stator during rotation in the blood pump generated thermal effects, which may have also had an impact on the red blood cells. In this study, the thermal effect was not considered during the simulation.

Conclusions
This study provides a complete multi-objective optimization process for a spiral centrifugal blood pump. The process contains model parameterization, Latin hypercube sampling, random forest prediction, a numerical simulation, and a multi-objective gray wolf optimization algorithm. In the first step, six impeller parameters are selected to control the shape of the impeller. Then, 240 sets of impeller parameters are randomly generated using the Latin hypercube sampling method, and these are numerically simulated to obtain a database for training the random forest prediction model. In the second step, the Pareto front is calculated using the multi-objective gray wolf optimization algorithm, and the highest point of the lift model (the A-point model), the lowest point of the HI model (the B-point model), and the middle point model (the C-point model) are selected as the optimized models. In the last step, the pressure cloud, the velocity flow line, the SSS distribution, the HI distribution, and the vortex distribution of all the optimized models are compared with the baseline model. The comparison results show that the pressure generation in the A-point model has increased by 41% and the HI value has increased by 105%. The pressure generation in the B-point model has increased by 19% and the HI value has decreased by 50%. The pressure generation in the C-point model has increased by 24% and the HI value has decreased by 48%, and the increase in the HI value is almost negligible as compared with the increase in pressure generation. Thus, the C-point model is the best-optimized model. The coupled optimization algorithm proposed in this paper exhibited a high prediction accuracy and fast optimization speed and can be used as a highly efficient optimization method for screw centrifugal blood pumps. In future research, the method proposed in this paper can be used to optimize other types of left ventricular assist devices such as axial flow pumps.

Data Availability Statement:
The data presented in this article are available upon request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.