Evaluation of novel-objective functions in the design optimization of a transonic rotor by using deep learning

Design optimization of transonic airfoils for rotary blades is a challenging subject that remarkably affects the stage and overall performance of axial-flow compressors. This paper describes a surrogate-based multi-objective optimization process over a transonic rotary blade. This blade works in the first high-pressure stage of a pre-designed industrial axial compressor. It experiences a massive separation behind an impinging shock wave over its suction side, resulting in very low efficiency of the whole stage. The key components of the current approach involve the application of novel-objective functions over the pressure distribution of airfoils, called the location of the shock wave and a flat-roof-top factor, to design supercritical airfoils. Moreover, to ensure the advantages of having an attached boundary layer and a high efficient blade, the area of separated boundary layer is also defined alongside other well-known objective functions related to the polar loss diagram. Notably, a sequential feed-forward multi-layer perceptron is designed to construct a mapping between airfoil geometrical variables and the objective functions. A numerical simulation of the whole compressor has shown an efficiency improvement of about 10% and 0.17% for the first stage and the whole compressor, respectively, and an attached boundary layer with a supercritical pressure distribution when employing the optimized rotor blade at the design stage.


Introduction
Deep learning methods introduce an intriguing, promising application within the community of design optimization of transonic airfoils in achieving a highpressure, highly efficient stage for the rotary machines and are now common practice in several research works (Sekar et al., 2019;Singh et al., 2017). In every optimization process, especially in the transonic flow regime, the tendency towards flow divergence by a numerical solver might increase as geometrical parameter variation increases due to high correlation of airfoil performance with minor variations of geometrical parameters as well as highly complex flow by the formation of Shock Wave Boundary Layer Interaction (SWBLI) on the suction side of airfoils. Such unstable structure of the SWBLI has been experimentally investigated over a supercritical airfoil (Lee, 2001;Masdari et al., 2020). Under such an unstable flow condition, CFD solvers are prone to create an unstable solution through each optimization process. Handling of this issue relating to the mentioned instabilities is performed by the introduction of Artificial Neural Networks (ANNs), attempting to construct a mapping between inputs and outputs to discover correlations between patterns or feature of designed data sets and airfoil aerodynamic performance. In this regard, Sessarego et al. (2020) described the application of a neural network with a gradient-based optimizer for the design optimization of a curved wind turbine blade using an aero-elastic simulator with synthetic inflow turbulence. Moreover, Yilmaz and German (2017) employed a convolutional neural networks (CNNs) approach for training the performance of an airfoil over a wide range of geometries. Various applications of deep learning methods in the prediction of fluids properties (Ardabili et al., 2018;Baghban et al., 2019;Ghalandari et al., 2019;Maleki et al., 2020;Shamshirband et al., 2020) have been applied gradually and still remain a challenging topic. Theoretically, airfoil design optimization is an infinitedimensional problem in which an airfoil is in practice parametrized with a finite set of geometrical design variables. The distribution of the variables involves making a decision on how they have to move and are usually determined based upon the experience of designers or sensitivity analysis. Airfoil parametrization methods have been developed through several research works, e.g. the Parametric Section (PARSEC) method (Sobieczky, 1999), the Hicks-Henne function (Hicks & Henne, 1978), Class-Shape Function Transformation (CST) (Kulfan, 2008) and Free-Form Deformation (FFD) (Lassila & Rozza, 2010). A comprehensive review of airfoil shape parametrization methods can be found in Masters et al. (2017). In the present article, each airfoil will be represented by the Bezier parameters of its camber and thickness lines. Apart from the method of parametrization, sensitivity analysis might also be used for dimensionality reduction. In this regard, many research works used a large collection of airfoils and singular value decompositions to produce a universal set of modes representing the deformation of airfoils shapes (Ghoman et al., 2012;Poole et al., 2015;Toal et al., 2010).
Apart from the flow solver and the type of optimization algorithm, the expressive capability of defined objective functions remains challenging for identifying the physics of transonic flow over airfoils and their aerodynamic performance in each optimization algorithm. Koller (Willer et al., 1999) described the performance of an airfoil by defining a corresponding function over the total pressure loss to design a 'flat' loss curve over 80% of its operating range. Schnoes and Nicke (2017) employed a single objective function optimized approach that considers the average total pressure loss of all operating points, and described methods of classifying a large set of optimized airfoils according to inlet Mach number, blade stagger angle, etc. in a database.
In this article, the first rotor blade of a designed axial compressor that experiences a massive separation behind a spanwise-distributed SWBLI over the suction side of the blade at its designed operating point is selected to be optimized through a surrogate-based Multi-objective Genetic Algorithm (MOGA). This paper aims, firstly, to predict transonic airfoil performance by describing new objective functions defined over the pressure distribution of the airfoil. These objectives are the location of the shock wave and a flat-roof-top factor defined over the supersonic bubble. Moreover, to ensure the advantage of having an attached boundary layer and a highly efficient blade, the area of the separated boundary layer is also defined alongside other well-known objective functions relating to the polar loss diagram.
Secondly, these objectives are employed in the MOGA to evaluate the ability of these novel-functions in design optimization. The airfoils are parametrized at seven different locations of a transonic rotor blade, three of which, in the vicinity of the hub, mid-span and near the shroud, are selected for representing the optimization results. A correlation-based sensitivity analysis is employed for gaining a better understanding of the physics of the flow and dimensionality reduction. The approach applied and the results obtained are provided in the following sections.

Basic design layout and information
A multi-stage transonic axial compressor was designed based on the mean-line method under ISO inlet conditions. In this machine, the first rotor blade, which experiences a massive separation behind a spanwisedistributed shock wave over the suction side of the blade, is selected to be optimized through a surrogate-based multi-objective optimization process. The blade requires as high a turning angle as possible to provide the required high pressure ratio at the design point. Figure 1 illustrates a conceptual layout of the assembeled rotors.

MISES solver
The MISES solver is a combination of FORTRAN codes provided and integrated by Drela and Youngren (1998). MISES has several advantages for design applications, including accuracy, speed, and a robust inverse design capability. Primarily, its speed makes this approach attractive for design optimization. A typical case is solved in 3-10 Newton cycles, indeed only a matter of minutes on a fast workstation. The accuracy also has a strong influence on the design process -the drag predictions from ISES are sufficiently good that they can be used reliably by a designer to optimize an airfoil in both on-design and off-design conditions.
Owing to some simplifications implemented in full Navier-Stokes schemes, MISES is able to simulate subsonic and low inflow angles more precisely with respect to their higher inflow angles. It incorporates a modified version of the bypass transition model, namely the Abu-Ghannam-Shaw (AGS) model, which can be activated by users in the code. In this article, a recommended grid topology for transonic airfoils near hubs, mid-spans, and shrouds is selected to meet the requirement of near orthogonality of the mesh cells. The cascade mesh consists of 30 streamlines with 80 and 60 grid points at the suction and pressure sides, respectively. For precise numerical calculations near the leading edge, the aspect ratio of the thick airfoils is considered to be 0.1, and for the thin ones it is set as 0.05 in the MISES grid setup. Figure 2 shows the grid lines for the airfoil geometry located at the mid-span position. During the optimization, the solver mesh will serve as the representation of the cascade geometry. In this condition, it is necessary to tackle away mesh deformation for different new airfoil geometries. To accomplish this, there are two available parameters in MISES named the mesh inlet and outlet angles, which are set according to the airfoil leading edge and trailing edge metal angles.
In the case of the validation of MISES solver results with some experimental tests conducted over a turbine cascade, Andrew and Kahveci (2007) carried out an experimental analysis over two different geometries of the turbine blades of an air breathing engine in General Electric and made a comparison between experimental and MISES results over their selected airfoil. Since the location of the shock wave in this paper as well as the flat-roof-top factor are used as objective functions, a further validation of the pressure distribution calculated by MISES by means of the experimental results conducted in Fuchs et al. (1993) is performed. According to that article, an inviscid-viscous interaction technique, derived from Denton and Walz methods, was utilized for the calculation of the blade-to-blade flow field over a designed transonic compressor cascade. As shown in Figure 3, MISES has an acceptable accuracy for prediction of the shock wave mean location. As the precise flow solver, moreover, a combination of the k − ω turbulence model and the Abu-Ghannam-Shaw (AGS) transition model (Abu-Ghannam & Shaw, 1980) based on Drela's modification (Drela & Aero-Astro, 1998), which has been successfully used in the MISES code, are implemented in the quasi-three-dimensional version of HSTAR (Arima et al., 1997). Furthermore, the reliability of the design process and the application of MISES as a flow solver in the design process of new airfoils have been investigated by Schreiber et al. (2004).  Since there was more concern about a reliable numerical simulation of the polar diagrams over compressor airfoils, a comparison between FLUENT R and MISES simulations is also conducted for an arbitrary airfoil of the axial compressor. Figure 4 represents the polar diagrams of an arbitrary thin airfoil simulated by MISES and FLUENT R at M ∞ = 0.5 and 0.8. From the figure, there is a remarkable discrepancy in the calculated loss and outflow angles at higher inflow angles and Mach numbers, respectively, which is due to the Kutta condition on the trailing edge of the airfoil which prevents the flow being completely solved in this region at high inflow angles and Mach numbers. In short, it can be inferred that MISES is reliable in predicting the aerodynamic performance of airfoils, and consequently a numerical dataset of the design optimization process is built.

Objective functions
Generally, Figure 5(A) shows a schematic of the objective functions defined over the polar loss of the airfoils. It is vitally important that the blade sections should operate with a low loss (ω) represented in Equation (1), wide operating range ( β) represented in Equation (2), and acceptable chock (ξ c ) and stall (ξ s ) margins shown in Equations (3) and (4), respectively, also shown in Figure 5(A). In this figure, the Design Point (DP) implies the minimum loss and the Operating Point (OP) is a condition where an airfoil is operating based on CFX simulation results. The dotted horizontal line in Figure 5(B) indicates the critical pressure coefficient over a transonic airfoil (i.e. the pressure coefficient where an isentropic flow would reach the sonic condition defined in Equation (5). This line can faithfully imply the shock wave location at a certain Mach number, and can be considered as an objective function to be maximized for the airfoils at higher inflow Mach number at mid-span as well as shroud positions. To achieve an attached boundary layer, it is necessary to minimize the area of separated flow. Such an objective function, which also implies that the separated boundary thickness itself is calculated based on Equation (6), is illustrated in Figure 5(C) at each iteration of the process. For a transonic rotor blade, the pressure distribution requirements over its airfoil at different spanwise locations alter from hub to shroud position. Thanks to the introduction of a ski-slope pressure distribution (Sieverding et al., 2003), which removes separation near the hub, the pressure distribution should have a peak less than 5% of the chord and allow diffusion to start from a region where the boundary layer thickness is still thin. In order to obtain better high-speed performance at the operating point for the transonic airfoils standing at midspan as well as shroud positions, it is desirable to shift the position of maximum thickness and camber towards the trailing edge of the airfoils to accomplish weakening of the shock wave on their suction surfaces. Furthermore, positioning of the maximum camber close to the maximum thickness location ensures that the region behind the shock wave, individually, will properly diffuse the flow. Such transonic airfoils, additionally, should utilize a flat-roof-top pressure distribution to avoid sudden acceleration of the flow upstream of the maximum thickness location. In this paper, a flat-roof-top property is defined over the pressure distribution of these airfoils by using Equation (7). Minimizing the front-part pressure distribution points' mean deviation with respect to a referencē C p is considered as an objective function.

Airfoil parametrization method
The three airfoils are parametrized using a five-ordered Bezier representation of their camber and thickness distribution lines as described in Pakatchian et al. (2019) and illustrated in Figure 6. Note that, if an airfoil has a common semi-circular trailing edge, a cut-off is required before importing into the MISES solver. Therefore, a blunt trailing edge is specified by the blade 'open', and the first and last coordinate points do not coincide. The design space consists of 37 geometrical variables (x i , y i ), seven objective functions and five constraints (V1, V2, V14, V16, V18). Figure 7 shows a boxplot of nondimensional range of the involved variables from hub to shroud of selected airfoils. In this figure, Vi stands for x i where (i = 2n + 1, n ∈ R), and for y i where (i = 2n, n ∈ R). For instance, (V1, V2) implies the (x, y) related to the first control point of the thickness distribution line. There are specified variables that are held constant, as determined by zero values in Figure 7, due to initial design considerations and tangency requirements of the Bezier representation; for example, (V7, V8) and (V5, V6), and (V15, V16) and (V17, V18) conform to this requirement. V14, V16, and V18 are related to the maximum thickness of the airfoil and are considered to be constant during optimization. The constants V25 and V27 coincide on the leading edge and the trailing edge, respectively, and determine the chord of the airfoil. Note that the meridional  chord is held constant in order to meet the requirements of spacing between the rotor and stator of the first stage. This spacing varies from hub to shroud respectively in the basic design. The fill factor and outflow angle are allowed to vary in a determined bound, in case thickening of the blade at the hub or even at the shroud are required.

Neural network design
In this study, sequential feed-forward neural networks with two hidden layers are designed to construct a mapping between airfoil geometrical variables and an objective function. Detailed information relating to the setup of neural networks is listed in Table 1.
Each neural network is trained over a design space, and is also responsible for predicting an objective function. In this regard, the Latin hypercube method is utilized for near-random sampling of the geometrical variables. The initial dataset including 500 members is divided to training, validation, and test data, which are normalized with 'MinMaxScaler' from the scikit-learn library in Python TM . All members are solved by MISES and then related objective functions are calculated. The training data are used during the learning process to fit the data-driven model. Since the validation dataset is also used to fit the model, it is also regarded as part of the training set. The test data are completely independent of the training data but follow the same probability distribution and are used to obtain the performance characteristics of the final model. In this study, 5 to 10% of data are used for the validation, 10 to 15% for testing, and the remaining dataset is applied for the training of each neural network.  Figure 9. Robustness of surrogate-based (A) and traditional (B) methods. Figure 8 shows the procedure of the approach applied, which is shown by numbers (1)-(9). In the first steps (1) and (2) the initial dataset is built over the design space using the Latin hypercube method. Then the objective functions relating to each member in this dataset are calculated by the MISES solver and used for initial training of the neural networks (3). This step takes approximately 13 h12 min. Note that, to avoid overfitting, each generated dataset in the current approach must not be generalized for different airfoil shapes at other spanwise locations due to different boundary conditions. Furthermore, the best weights during the training process are found using the checkpoint module in the KERAS library from Python. To avoid overtraining, a dropout layer is also considered as a regularization method in which the fraction of the input units to drop is 0.02. In the next step (4), the trained neural networks rather than the MISES solver predict the involved objective functions for the genetic algorithm in each optimization iteration. When the optimization terminates or reaches a certain tolerance, a Pareto front will be constructed containing a reduced-elite of optimum airfoils. This elite is imported into the MISES solver (6) to calculate their performance (7) and validate the predicted performance (5) by means of the neural networks. Then the error of each objective function can be calculated by Equation (8):

Design optimization algorithm
If one objective error exceeds a specified limit (10e-3), this elite will be imported into the initial dataset and the neural networks will be retrained until the errors between the calculated and estimated objective functions satisfy the limit. In this study, the convergence time speed as well as the robustness of the surrogate base and the traditional method are discussed. The main advantage of using this surrogate-base optimization is that the convergence time speed dramatically reduces. In this regard, the total computational time of the surrogate-based optimization is presented in Equation (9), where t DOE stands for the time required to solve the initial data by MISES. Training each neural network takes t train . Moreover, t SGA and t val. refer to the computational time related to running the genetic algorithm and the validation steps, respectively.
In the case of the traditional optimization approach, the MISES solver would directly calculate the objective functions for all generation, N GEN , each of which consists of N case new cases. One can, therefore, estimate the computational time from Equation (13): According the information in Table 2, Equations (14) and (15) can be derived. In these equations, N I stands for the number of process iterations, i.e. steps (1)-(8) in Figure 8: In this approach, an acceptable error limit is usually attained at N I ≤ 4 for the optimization of each airfoil. So the surrogate-based optimization improves the convergence speed by about 70% compared with the traditional method.
The robustness of the traditional as well as surrogatebased optimization methods is studied for 10 individual optimization runs. The values of the selected objective functions are illustrated in Figure 9 using box plots. From this figure it can be inferred that both methods are able to converge the optimization algorithm to acceptable local minima with very similar improvement intervals.

Results
The shock wave strength in the blade-to-blade domain strongly depends on for a certain inflow Mach number. Figure 10 shows the Mach number contour for an airfoil placed in the vicinity of the hub. It is seen that the lowest back-pressure forms a fully pitchwise-extended shock wave plus a massive separated flow behind it. As the backpressure rises, the shock wave location moves towards the leading edge of the airfoil. In this flow regime, the cascade is totally in stall mode. It should be noted that, in this research, the three selected base airfoils experience a similar flow condition, about p2/p01 = 0.74 and p2/p01 = 0.81 as shown in Figure 10, at their design points. Figures 11(A)-11(D) show scatter plots of the fitness functions related to the initial dataset grouped by the location of the shock wave, X/C s , and separation factor, λ. Figure 11(B) shows that there are some airfoils that operate in a high design loss condition while their separation factor, λ, is low with respect to other airfoils in the population. These special airfoils are shown by the letters O1 and O2 in Figure 11(B). It is inferred that the design loss might be independent of the separated boundary layer factor λ and, indeed, minimizing the design loss might not ensure having an attached boundary layer over the suction side for such airfoils. Note that this phenomenon is not observed in the design space of other airfoils at the mid-span and shroud with higher inlet Mach number. To cope with this, it is vital to consider λ as an objective function during optimization of the airfoil near the hub. Figures 12(A)-12(D) show the same scatter plots for airfoils at the mid-span. It is seen that higher operating ranges are achieved at the expense of higher design loss, as shown in Figure 12(A). Moreover, the flatroof-top pressure distribution, i.e. low values of κ, is often formed for those cases with a shock wave location near 0.4 of the chord. Indeed, those airfoils with the most aft shock-wave position would be desirable if they had low values of κ. The trend shows that the rate of change of the factor κ is proportional to the factorλ, as shown by smaller markers in Figure 12(B). In fact, those airfoils that do not meet the flat-roof-top requirement tend to have massive separated boundary layers over their suction side. From Figure 12(C) it can be inferred that those airfoils with large chock margins tend to reproduce lower operating loss in this design space. Figures 13(A)-13(D) show the same scatter plots for airfoils near the shroud position. It is seen from Figure 13(B) that, in contrast to the presented trend in the design space at the hub (as shown in Figure 11), there is a positive correlation between the separated boundary layer area, λ, and the design loss, ω D . The objective functions are influenced by some major geometrical variables. Likewise, there are some minor ones that can be neglected. The calculation of the proposed correlation will lead to a reduction in the size of the search space and to an acceleration of the process. In this article, the correlation of the design variables and the objective functions is used to reduce the dimensions of the design space by eliminating some minor geometrical variables in the optimization algorithm.
The robustness of the optimization depends on the number of design variables and the process speed. Once the number of design variables increases, the optimizer is prone to explore airfoils with more curvature variations, contributing to difficulties for the flow solver. In this study, a correlation-based sensitivity analysis is employed for efficient computation and dimensionality reduction. Figures 14(A)-14(D) show the heat maps of all variables and objective functions at the neighbourhood of the shroud, all of which are symmetric matrices including arrays representing the correlation coefficients. For instance, as shown in Figure 14(A), variable no. 9 (V9) could not make a considerable change in ξ c , since its correlation coefficient is approximately −0.13. Figure 14 shows that the rate of change of V13, which stands for the maximum thickness location, could be inversely proportional to the design loss. As the maximum thickness location increases and reaches towards the leading edge, both design and mean loss values could decrease. V29, V30, and V32, which are related to the front cambering of an airfoil, have considerable correlation with both mean and design loss values, as shown in Figure 14(C). Moreover, V33, V34, and V36, which stand for aft-cambering, could increase the loss values as they increase. In this study, some variables for which their correlation with respect to the specified objective functions is within the bounds of ±0.1 are removed from the space of involved variables. The same process is applied for other airfoils at mid-span and hub locations. Figures 15(A)-15(G) show the loss-monitor of the involved neural networks inside the scatter plots of the network's predicted performance during the initial training for each epoch. The loss function is considered as mean square error when setting up sequential neural networks. As may be seen in Figures 15(A)-15(G), all of the obtained values of R 2 , related to the trained and the test data, are above 0.95 and 0.88, respectively. Note that the initial training is prepared based on the initial dataset population and this dataset utilizes the merit of being increased by means of the Pareto front members in each iteration of the optimization process shown in Figure 8. This merit and considering a 'Mod-elCheckpoint' for saving the best weights as well as a dropout layer during 400 epochs can positively affect the accuracy of the final neural networks with two layers in predicting aerodynamic performance. For instance, comparing the prediction performance related to the flatroof-top network (G in Figure 15) with the results of pressure distribution related to the optimum airfoil in Figure 17(B) or Figure 18(B), implies the potential of this retrained network for accurately predicting the flat-rooftop function. Figure 16 shows the final result of the optimization in comparison with the base airfoil at the neighbourhood of the hub. It is seen that the selected airfoil from the Pareto front with a relatively higher operating loss compared with the base airfoil, as shown in Figure 16(A), is able to produce an attached boundary layer. As may be seen in Figure 16(B), according the contour of the Mach number in the pitch-wise plane, there is a relatively thick separated boundary layer over the suction side of the base airfoil due to an intense adverse pressure gradient in the pressure recovery region. By deliberately selecting a vast range of displacement for the variables relating to the aft-thickness variation (V19, V21, and V26), the optimization process is permitted to search the space relating to the thicker airfoils. Such permission enables the blade to tolerate more centrifugal tension. Furthermore, the maximum thickness and camber locations of the airfoil have been shifted towards the leading edge, which results in generating a ski-slope airfoil with a peak pressure close to its leading edge. Note that the optimization algorithm is allowed to make a reduction in the maximum camber value to enable it to remove separation at the expense of decreasing the pressure ratio just for this airfoil near the hub position.   Figures 17 and 18 show the same results for the airfoils at the mid-span and at the neighbouring shroud, where the optimization goals are selected as an attached boundary layer behind a weak shock wave as well as the formation of a flat-roof-top pressure distribution to avoid reaccelerating flow. As shown in Figure 17(A), a more than 80% reduction in the loss of operating points in comparison with the base airfoil is obtained in the ω versus β 1 diagram. As can be seen from the β 2 versus β 1 diagram, the increment of outflow angle of the operating points with respect to the base airfoil implies a reduction in the turning angle of the flow. A significant change in the airfoil wake can be inferred from the stream tube number versus the M out diagram, which implies an attached boundary layer over the suction side of the airfoil. It can be also seen from Figure 16(B) that there is no thick separated boundary layer over the suction side of the optimized airfoil due to the above-mentioned refinements. Moreover, the strong shock wave as shown in the centre of this figure is weakened thanks to the formation of a flat-roof-top pressure distribution over the optimized airfoil.
For airfoils in the neighbourhood of the shroud, as shown in Figure 18, the design loss, flat-roof-top factor, and shock wave location are selected as the objectives of the genetic algorithm. These objectives ensure that a supercritical airfoil at higher Mach number than the midspan airfoil, with lower design loss with respect to the base, as well as with an aft-position shock wave as far as possible will be achieved. This condition allows the aft-part of the airfoil to diffuse the flow near the shroud position. From the ω versus β 1 diagram of Figure 18(A), it is seen that this airfoil experiences a thick separated boundary layer behind a strong shock wave and high loss due to high inflow Mach number near the tip of the blade.
As shown in the β 2 versus β 1 diagram, the outflow angles at different operating points are increased with respect to the base airfoil after the optimization process. Moreover, p2/p01 versus Q shows that such an optimum passage is able to pass more mass flow per unit area than that at the base in all off-design points. A similar smoothness in the wake of the optimum airfoil can be seen from the streamlines versus M out diagram, which implies an attached boundary layer over the suction side of this airfoil. Figure 19(A) shows the local and global separation over the base compressor blades from IGV to stator no. 3. These regions are classified as turbulent, laminar, and corner stall separation regions. Additionally, the contour of velocity at those different spanwise locations, Y/C = 0, Y/C = 25, Y/C = 0.50, Y/C = 0.95, of the first rotor blade is shown. The turbulent separation along the spanwise direction arises from the boundary layer separation behind the shock wave at the pressure recovery region of the blade. The pattern of corner stall explained by Gbadebo et al. (2005) is also seen near the hub. Such patterns are formed as a consequence of interaction between the shroud and blade surface turbulent boundary layers. Such patterns, and moreover the existence of tip clearance, might lead to induced flow from the pressure side of a rotary blade toward its suction side, thereby forming a tip vortex. Figure 19(B) shows CFX simulation results of a configuration in which rotor no. 1 has been optimized. As shown in the contour of Mach number over its airfoils at different spanwise locations, the massive separation behind the shock waves of each airfoil is almost removed due to the smearing of the impinging shock wave formed on the suction side of the blade.     Figure 20(A) shows that, at a constantly maintained inlet Mach number distribution, a smooth outlet Mach number is achieved after optimizing the rotor blade. In fact, a smooth outflow Mach number guarantees a smooth wake profile of the blade. Moreover, Figure 19(B) shows inlet and outlet flow angles for both rotors. It can be inferred that the outflow angle increases over a half of the span after optimizing the rotor blade, which means that the optimization algorithm has reached to a point where the optimized blade will remove separated flow, reducing the total loss as well as the outlet temperature -as shown in Figures 20(C) and 20(D), respectively -at the expense of a small decrease in the pressure ratio of the rotor for the same inlet flow conditions. In total, Table 3 contains information on the numerical simulation of the whole compressor after replacing the optimized rotor with the base one, from three different viewpoints. As shown, the calculated pressure ratios for the three steps have remained fairly constant. In contrast, the effect of optimization is considerable in all three steps. Firstly, improvement in the first row efficiency is about 1.8%, which shows a remarkable effect on the first stage efficiency with a more than 10% improvement at the expense of a slight decrease in the first stage pressure ratio and a small increase in inlet mass flow by 0.74 and 0.25%, respectively. Note that full compressor results show a 0.17% enhancement in total efficiency. Although the surge margin is highly dependent on the tendency to flow divergence relating to the last stages at higher pressure ratios, a notable increase of about 0.51% in the surge margin has been achieved after optimization of this rotor. This implies an improvement in the velocity triangles for the last stages at off-design operating points.

Conclusions
In the present work, three different novel-objective functions, namely the area of the separated boundary layer, the location of the shock wave, and the flat-roof-top pressure distribution factor, were introduced in a surrogatebased optimization process. These objectives alongside other well-known objective functions were purposely  selected for the optimization of each airfoil across the span of a transonic rotor. The space of airfoil variables, designed by the Latin hypercube method and described by the scatter plots of the objective functions, demonstrated that the area of separated boundary layer correlated directly with the design loss for the airfoil near the shroud, where the presented airfoil experienced the highest inflow Mach number, in contrast to the airfoil near the hub position. The surrogate-based optimization results showed that these objectives are efficiently able to provide a supercritical pressure distribution for airfoils at mid-span and in positions neighbouring the shroud, and interestingly are able to remove the massive separation behind the impinging shock wave, possibly increasing the efficiency of the transonic rotor. The simulation of the whole compressor with the optimum rotor blades showed that supercritical transonic airfoils could dramatically increase the efficiency of the rotor blades by about 1.8% with the desired attached boundary layer at the expense of a slight decrease in the first stage pressure ratio and a small increase in inlet mass flow by −0.74 and 0.25%, respectively. As the main disadvantage of surrogate-based design optimization, it should be noted that surrogate-based design optimization might limit the maximum attainable pressure ratio of blades. Moreover, for upcoming research in this field it is recommended that a combination of experimental and validated numerical data be gathered in a dataset, the classification of which can also be done by employing relevant classification algorithms such as the support vector machine method. This could eventually lead to the derivation of more accurate as well as general data-driven models.