Multi-Objective Optimization of Cascade Blade Profile Based on Reinforcement Learning

Abstract

The multi-objective optimization of compressor cascade rotor blade is important for aero engine design. Many conventional approaches are thus proposed; however, they lack a methodology for utilizing existing design data/experiences to guide actual design. Therefore, the conventional methods require and consume large computational resources due to their need for large numbers of stochastic cases for determining optimization direction in the design space of problem. This paper proposed a Reinforcement Learning method as a new approach for compressor blade multi-objective optimization. By using Deep Deterministic Policy Gradient (DDPG), the approach modifies the blade profile as an intelligent designer according to the design policy: it learns the design experience of cascade blade as accumulated knowledge from interaction with the computation-based environment; the design policy can thus be updated. The accumulated computational data is therefore transformed into design experience and policies, which are directly applied to the cascade optimization, and the good-performance profiles can be thus approached. In a case study provided in this paper, the proposed approach is applied on a blade profile, which is thus optimized in terms of total pressure loss and laminar flow area. Compared with the initial profile, the total pressure loss coefficient is reduced by $3.59\%$, and the relative laminar flow area at the suction surface is improved by $25.4\%$.

1. Introduction

For complex aerodynamic design cases, such as turbomachinery blades, single optimization objective is not enough, according to the complicated flow phenomena. Multi-objective optimization is usually required and becomes more and more essential in the aerodynamic shape optimization. The compressor cascade design is an optimization with multiple objectives, including, e.g., total pressure loss, outlet flow angles (which need to be constrained during total pressure loss opitmizatiton), laminar flow area increasing (for skin friction drag reduction [1,2]), etc. [3]. The increased laminar flow area not only benefits the total pressure loss but also results in the decrease of the turbulence of the wake, which improves the quality of the incoming flow of downstream blades. Furthermore, the optimization of the compressor cascade is not only at one design point, but at multiple operating conditions [3], which makes the multi-objective optimization of cascade more complicated and a big challenge to the optimization methods. As the number of optimization objectives increases, the computation required for optimization increases, and there will be difficulty to obtain the globally optimized result [4]. Since aerodynamic shape design requires good understanding of aerodynamics and sufficient experience, one solution to the above problem is to provide the guidance of experience to the optimization process for a better use of accumulated data, thus saving computation resource. Therefore, a more global, efficient, and intelligent multi-objective optimization method based on accumulated design experience is proposed in this paper.

The optimization methods used in aerodynamic shape design are conventionally gradient-based methods and stochastic methods. Gradient-based methods rely on the gradient of the objection function, usually converge rapidly, and have low requirements for the computational cost. Adjoint methods are typical gradient-based methods and have been used in many aerodynamic shape optimizations [5,6,7], such as turbomachinery optimization [8,9,10]. However, gradient-based methods easily get trapped in local optima and perform poorly for strongly non-linear or discontinuous objective functions [11]. With the development of computing power, stochastic methods are more and more widely used in aerodynamic shape optimization, e.g., Particle Swarm Optimization (PSO) [12], Simulated Annealing (SA) [13], Genetic Algorithm (GA) [14], and Evolutionary Algorithm (EA) [15], which are good at finding global optima, in many multi-objective optimizations [16,17,18,19]. Multi-objective stochastic methods have been used in the optimization of turbomachinery (Multi-objective Genetic Algorithm (MOGA) [20,21,22], Multi-objective Evolutionary Algorithm (MOEA) [23,24,25]). However, stochastic methods need large numbers of stochastic cases to decide evolutionary trends via probability-based selection, whereas the existing data is not capable of being utilized in the above process to determine optimization direction. It results in the drawbacks of stochastic methods, like its slow convergence and high computational cost (especially for high-fidelity computational fluid dynamics (CFD)). These limitations are overwhelmingly disturbing in global multi-objective optimization with a larger number of design parameters.

In recent years, many artificial intelligence techniques have caught the attention of researchers. Reinforcement Learning (RL) is one such approach that interacts with the environment according to policies, learns from experience, and improves policies. For example, the RL algorithm, Q-learning [26], computes the Q value of every state and action pair and stores it in one table. This value-based algorithm completes each iterative calculation by looking up the table, so it is suitable for the conditions wherein the space of state and action are discrete and the dimensions are not high. The experience-based learning feature of RL makes it advantageous for optimization application [27]. In its application on aerodynamic shape design, the optimizer modifies the shape according to the design policy as an intelligent designer, interacts with the environment (obtaining aerodynamic data through CFD simulation or surrogate model), evaluates the policy according to the aerodynamic data, learns design experience, and updates the design policy. Some RL algorithms have been developed and applied in aerospace optimization. The application of RL on the direct optimization of aerodynamic shape was presented by Viquerat et al. [28] and led to new generic optimization strategies. Li et al. [29] exploited an RL algorithm to learn the policy for drag reduction of supercritical airfoils, where the mean drag reduction of 200 airfoils had been effectively accomplished.

Combined with Deep Learning, RL algorithm can extract features from data in high-dimensional state spaces. Recently, Google Deepmind developed Deep Reinforcement Learning (DRL). Deep Q-learning (DQN) [30] is one such example in that it is the improvement of Q-learning with three innovations: (1) DQN uses deep neural networks to approximate the value function; (2) DQN uses a replay buffer to reduce the correlation between samples; and (3) DQN sets up an independent target network to handle the error in the temporal difference. However, DQN is restricted to cases of the discrete and low-dimensional state and action. For high-dimensional continuous action problems, Deterministic Policy Gradient (DPG) was proposed. Different from the value-based learning in DQN, DPG learning is policy-based in that directly optimizes deterministic policies in form of parameters of the neural networks via gradient-based methods [31]. Combined with Deep Learning, Deep Deterministic Policy Gradient (DDPG) can thus be utilized in cases of more complicated action [32]. In its applications on aerodynamic optimizations, Yan et al. [33] proposed a DDPG-based aerodynamic shape optimization. They optimized a preliminary missile shape design involving lift coefficient and drag coefficient, where the DDPG network was applied to guide the shape parameters variation by interacting with the environment based on CFD and semi-empirical method DATCOM.

In this paper, DDPG is applied in a multi-objective optimization of compressor cascade blade profile. The geometric parameters of blade profile are taken as State, and the modification of parameters are taken as Action. From interaction with the simulation-based environment, the aerodynamic parameters corresponding to a certain State (profile geometry) can be achieved. The DDPG network can learn the design experience of cascade blade with the feedback according to the aerodynamic parameters, which can guide the modification of the profile optimization. The existing computational data is now utilized to build the policy to guide the optimization direction in the design space, which saves many calculations. With the learnt policy, the DDPG network can optimize the blade profile as an experienced designer to achieve the good-performance profile.

In order to reduce the computation and shorten the optimization period, an ANN (artificial neural network)-based surrogate model is used in this paper to approximate high-fidelity CFD simulations [34]. The model is trained with profile geometric parameters and aerodynamic parameters by numerical calculation. The surrogate model is utilized to feed back the aerodynamic information to the DDPG network as the environment in the optimization. To achieve the multiple objectives interconnection, as well as an appropriate objective function, analysis on the objective space, composed of optimization aerodynamic parameters, is done with Principal Component Analysis (PCA) [35]. The flowchart of the optimization is illustrated in Figure 1.

This paper is organized as follows: In Section 2, the parameterization of the blade profile is introduced. In Section 3, the cascade flow field calculation is validated. Section 4 introduces the surrogate model based on ANN. Section 5 introduces the DDPG algorithm, and the optimization application is presented. Section 6 shows the optimization results with some discussions, and the paper is concluded in Section 7.

2. Parameterization of the Blade Profile

The size of the database for training networks grows fast with the increased number of the input and output variables of the networks. For computation reduction, the number of parameters is required to be as few as possible. One approach is to parameterize the perturbation instead of complete blade profile. It is necessary to have an adequate number of variables to ensure geometric accuracy if parameterization is directly applied to complete description of the blade profile. However, the required number is less if the parameterization is applied on the perturbation of the blade profile.

The parameterization is used to obtain a library of modified profiles with related parameters. The parameterization of blade profiles in this paper is achieved by using the Hicks-Henne bump functions to perturb the original airfoil with assigned coefficients. For situations without fitting the leading edge and trailing edge, the shape functions are defined as:

$$b_i\left(x\right)=si{n}^w\left(\pi x^{ln0.5/ln{x}_i}\right),$$

(1)

where $x_i$ represents the normalized position of the bump, and w controls the width of the bump (w equals 3 in this paper). The deformed profile curve can thus be described by:

$$y\left(x\right)=y_{basic}\left(x\right)+\sum _{i=1}^N{\alpha }_i{b}_i\left(x\right),$$

(2)

where ${\alpha }_i$ presents the magnitude of bumps.

According to previous numerical simulation of the internal flow field of the initial cascade done by the authors, the front middle of the blade suction surface is where the transition and a strong passage shock wave takes place. In order to control the transition position and the passage shock, for both suction surface and pressure surface, five bump functions are used to perturb the profile, and the $x_i$ equals $0.2,0.3,0.4,0.5,0.6$$(i=1,2,3,4,5)$ in Equation (1). The shape functions of profile perturbation are shown in Figure 2.

The geometric variable consists of ten coefficients of bump functions ${\alpha }_i$. With Latin Hypercube Sampling (LHS), the geometric parameters database is generated, thereby generating the blade profiles library. Ranges of variation of each coefficient are given in

Table 1. The range of bump function coefficients ${\alpha }_i$.

${\mathit{x}}_{\mathit{i}}$	${\mathit{\mu }}_{{\mathit{\alpha }}_{\mathit{i}}}$	${\mathit{\sigma }}_{{\mathit{\alpha }}_{\mathit{i}}}$
0.2 (SS)	−0.12	0.046
0.3 (SS)	0	0.046
0.4 (SS)	0.06	0.046
0.5 (SS)	0.04	0.046
0.6 (SS)	0	0.046
0.2 (PS)	−0.04	0.046
0.3 (PS)	−0.04	0.046
0.4 (PS)	−0.04	0.046
0.5 (PS)	−0.04	0.046
0.6 (PS)	−0.04	0.046

, where ${\mu }_{{\alpha }_i}$ and ${\sigma }_{{\alpha }_i}$ represent expectation and standard deviation of the normal distribution of bump coefficient ${\alpha }_i$. The distribution of the profiles is shown in Figure 3.

3. Numerical Method

3.1. Governing Equations, Turbulence, and Transition Model

The governing equations of general Reynolds-averaged Navier–-Stokes (RANS) equations are as follows:

$$\frac{\partial }{\partial t}{\int }_{\mathsf{\Omega }}\mathbf{U}d\mathsf{\Omega }+{\int }_S{\overrightarrow{\mathbf{F}}}_I·d\overrightarrow{S}+{\int }_S{\overrightarrow{\mathbf{F}}}_V·d\overrightarrow{S}={\int }_{\mathsf{\Omega }}\mathbf{S}d\mathsf{\Omega },$$

(3)

where $\mathsf{\Omega }$ is the volume, and $\mathbf{S}$ is the surface of the control volume. The conservative variables $\mathbf{U}$, the inviscid flux ${\overrightarrow{\mathbf{F}}}_I$, and the viscid flux ${\overrightarrow{\mathbf{F}}}_V$ in the time-averaged form can be written as:

$$\mathbf{U}=\left[\begin{array}{c}\overline{\rho }\\ \overline{\rho }{\tilde{v}}_1\\ \overline{\rho }{\tilde{v}}_2\\ \overline{\rho }{\tilde{v}}_3\\ \overline{\rho }\tilde{E}\end{array}\right]{\overrightarrow{\mathbf{F}}}_I=\left[\begin{array}{c}\overline{\rho }\widetilde{v_i}\\ \overline{\rho }\widetilde{v_1}{\tilde{v}}_i+\overline{p^{\ast }}{\delta }_{1i}\\ \overline{\rho }\widetilde{v_2}{\tilde{v}}_i+\overline{p^{\ast }}{\delta }_{2i}\\ \overline{\rho }\widetilde{v_3}{\tilde{v}}_i+\overline{p^{\ast }}{\delta }_{3i}\\ (\overline{\rho }\tilde{E}+\overline{p^{\ast }})\tilde{v}\end{array}\right]-{\overrightarrow{\mathbf{F}}}_V=\left[\begin{array}{c}0\\ {\overline{\tau }}_{i1}\\ {\overline{\tau }}_{i2}\\ {\overline{\tau }}_{i3}\\ {\overline{q}}_i+{\tilde{v}}_i{\overline{\tau }}_{ij}\end{array}\right].$$

(4)

Based on Boussinesq’s assumption, the Reynolds stress can be calculated with turbulent eddy viscosity ${\mu }_t$, which is obtained by solving turbulence models.

For turbulence calculation, Spalart-Allmaras turbulence model (SA) coupled with Abu-Ghannam and Shaw transition model (AGS) is used in this paper. SA turbulence model is a one-equation turbulence model proposed and developed by Spalart and Allmaras in 1992 [36]. SA model is often used in the calculation of the internal flow field of turbomachinery ([37,38,39,40], [etc.]) because of its robustness and its ability to treat complex flows [41,42].

The principle of SA model is aimed at solving an additional transport equation for kinematic eddy viscosity. In a non-conservative manner, the equation contains an advective, a diffusive, and a source term. The transport equation is as follows:

$$\frac{\partial \tilde{\nu }}{\partial t}+\overrightarrow{v}·\nabla \tilde{\nu }=\frac{1}{\sigma }\{\nabla ·[(\nu +(1+c_{b2})\tilde{\nu })\nabla \tilde{\nu }]-c_{b2}\tilde{\nu }∆\tilde{\nu }\}+S_T.$$

(5)

The turbulent viscosity is given by

$${\nu }_t=\tilde{\nu }{f}_{v1}$$

(6)

and

$$\left\{\begin{array}{c}{f}_{v1}=\frac{{\chi }^3}{{\chi }^3+c_{v1}^3},\chi =\frac{\tilde{\nu }}{\nu }\\ S_T=\tilde{\nu }P\left(\tilde{\nu }\right)-\tilde{\nu }D\left(\tilde{\nu }\right),\tilde{\nu }P\left(\tilde{\nu }\right)=c_{b1}\tilde{S}\tilde{\nu },\tilde{\nu }D\left(\tilde{\nu }\right)=c_{w1}{f}_w{\left(\frac{\tilde{\nu }}{d}\right)}^2\\ \tilde{S}=S{f}_{v3}+\frac{\tilde{\nu }}{{\kappa }^2{d}^2}{f}_{v2},f_{v2}=\frac{1}{{(1+\chi /c_{v2})}^3},f_{v3}=\frac{(1+\chi f_{v1})(1-f_{v2})}{\chi }\\ f_w=g{\left(\frac{1+c_{w3}^6}{g^6+c_{w3}^6}\right)}^{\frac{1}{6}},g=r+c_{w2}(r^6-r),r=\frac{\tilde{\nu }}{\tilde{S}{\kappa }^2{d}^2}\end{array},\right.$$

(7)

where $\tilde{\nu }$ is the turbulent working variable, $\nu$ is the molecular viscosity, d is the distance to the closest wall, and S denotes the magnitude of vorticity. The constants arising in the model are referred to Reference [42].

Abu-Ghannam and Shaw [43] investigated the experimental data of transition on a flat plate with pressure gradients and proposed the model to compute the location of transition from the flow solution by using empirical relations related to external parameters. The AGS model dictates that transition starts at a momentum thickness Reynolds number:

$$R{e}_{\mathsf{\Theta }s}=163+e^{F\left({\lambda }_{\mathsf{\Theta }}\right)-\frac{F\left({\lambda }_{\mathsf{\Theta }}\right)}{6.91}Tu}.$$

(8)

The dimensionless pressure gradient is defined as:

$${\lambda }_{\mathsf{\Theta }}=\frac{{\mathsf{\Theta }}^2}{\nu }\frac{d{U}_e}{ds},$$

(9)

where $U_e$ is the velocity at the edge of boundary layer, S is the streamwise distance from the leading edge, and $\mathsf{\Theta }$ denotes the momentum thickness in the laminar region.

For the SA model, the free stream turbulence level $Tu$ is computed based on the upstream eddy viscosity and assuming a turbulent length scale of about $10\%$ of the reference length scale:

$$Tu=86\sqrt{\frac{{\nu }_t}{U{L}_{ref}}}.$$

(10)

The function $F\left({\lambda }_{\mathsf{\Theta }}\right)$ depends on the pressure gradient:

$$F\left({\lambda }_{\mathsf{\Theta }}\right)=\left\{\begin{array}{c}6.91+12.75{\lambda }_{\mathsf{\Theta }}+63.64{\lambda }_{\mathsf{\Theta }}^2,{\lambda }_{\mathsf{\Theta }}<0\\ 6.91+2.48{\lambda }_{\mathsf{\Theta }}-12.27{\lambda }_{\mathsf{\Theta }}^2,{\lambda }_{\mathsf{\Theta }}\ge 0\end{array}.\right.$$

(11)

Therefore, according to these relations, transition is promoted in the adverse pressure gradient, whereas it is retarded in the favorable pressure gradient. The AGS model proves sufficiently accurate [44,45,46] and is used in predictions of transition in turbomachinery [47,48,49,50].

3.2. Cascade Performance Calculation

The three-dimensional compressible Navier–Stokes equation is solved by using a finite-volume based RANS solver NUMECA FineTurbo. Multigrid technique is used during the solving process in order to accelerate the convergence.

To test the accuracy of the transition position prediction, the VKI turbine blade is chosen as the validation case. The VKI turbine blade, as is shown in Figure 4, has been investigated a lot, both numerically and experimentally, and the experimental data on its transition position can be referred to in Ubaldi [51]. The blade-to-blade grid used to validate the calculation is shown in Figure 5.

Figure 6 shows the experimental and the calculated data of the transition position of VKI turbine. The experiment was carried out by Ubaldi [51]. Details about experimental measurement can be referred to in Appendix A.1. With Laser Doppler Velocimetry (LDV), the boundary layer integral parameters distribution is measured, and the start of transition (ST), transition position (TP), and end of transition (ET), marked in Figure 6, were concluded by Ubaldi, according to the boundary layer integral parameter distributions from the experiment [51]. The intermittency factor is used to trigger the transition locally, which is computed on each point of the blade surface and at each iteration via the solver. The intermittency factor is 0 before transition and 1 after the transition onset. According to the calculated intermittency factor distribution in Figure 6, the computational transition position fits well with the experimental transition position (magenta dashed line in Figure 6), which validates the accuracy of the calculation of transition position.

The geometry of the initial cascade blade profile to be optimized in this paper is shown in Figure 7, with the basic parameters listed in

Table 2. Geometric parameters of cascade.

Parameter	Value
Chord length c (mm)	109.5
Solidity	$1.792$
Blade geometric inlet angle ${\beta }_{1k}$ (${}^{\circ }$)	50.0
Blade geometric outlet angle ${\beta }_{2k}$ (${}^{\circ }$)	11.4
Blade camber angle (${}^{\circ }$)	$38.6$
Blade stagger angle (${}^{\circ }$)	$32.8$

. The blade-to-blade grid is generated as a structured grid, with the topology shown in Figure 8. The length of the inlet domain is as long as the chord line, and the outlet domain is twice as long as the chord line.

A three-dimensional grid is used as the cascade grid, with five points across the span of the cascade, generated via NUMECA AutoGrid. The boundary condition at the hub and shroud is set to the mirror boundary condition, to eliminate the three dimensional effect. The grid points distribution is listed in

Table 3. Grid parameters.

Block Name	Type	Points $(\mathit{x},\mathit{\theta },\mathit{r})$	Total Points
Inlet	H	$193\times 41\times 5$	39,565
Up blade region	H	$721\times 29\times 5$	104,545
Blade passage	O	$673\times 33\times 5$	111,045
Down blade region	H	$721\times 29\times 5$	104,545
Outlet	H	$233\times 41\times 5$	47,765
Total			407,465

. The total number of grid points is $407,465$. The height of the cells at the walls is set to $3\times {10}^{-6}$ m to ensure that $y^{+}$ is below 10. The 2D cross-sectional view of the computational grid is shown in Figure 9.

At the inlet boundary, the total pressure, the total temperature and the flow angle are specified. The total pressure of 119,600 Pa and the total temperature of 293 K are given at the inlet boundary. At the outlet boundary, the static pressure at standard sea level was specified. On the solid wall boundary, the no-slip and adiabatic conditions were used. In these conditions, the cascade flow is examined with an inlet Mach number of $0.8$ at the operating point (${\beta }_1={45.66}^{\circ }$). From the simulation, the relative laminar region at the operating point is $34.56\%$ at the suction surface (SS) and $29.26\%$ at the pressure surface (PS).

In this paper, total pressure loss coefficient $\overline{\omega }$ is defined as: $\overline{\omega }=\frac{p_0^{in}-p_0^{out}}{p_0^{in}-p^{in}}$, where $p_0^{in}$ and $p_0^{out}$ are the total pressure of the inflow and outflow. $p^{in}$ is the static pressure of the inflow. From multiple computations with the inlet flow of the same Mach number and different ${\beta }_1$, the aerodynamic performance curve can be obtained, as in Figure 10. The operating point is at the inlet flow angel ${\beta }_1$ of ${45.66}^{\circ }$, and the operating range is from ${45.16}^{\circ }$ to ${48.66}^{\circ }$ for ${\beta }_1$. The result of grid independence test is shown in Figure 11.

4. Surrogate Model Based on ANN

The optimization, especially with multiple objectives, are computationally intensive and time-costly. Using a surrogate model to approximate and replace high-fidelity CFD simulations is necessary for the complicated optimization of aerodynamic shapes. Surrogate models can provide aerodynamic information rapidly, which reduces the computation during optimization and shortens optimization period [52,53]. It is attractive for complex multi-objective optimizations for turbomachineries and is used in blade shape optimizations [54,55]. ANN has excellent nonlinear mapping and learning ability, which makes it advantageous in approximating the strong nonlinear relationship bewteen the geometric shape and multiple simulated aerodynamic parameters [34,56]. Using ANN as the surrogate model, with prediction instead of computation, eliminates the need for CFD simulations during optimization, which improves optimization efficiency and automation [56].

In this paper, a surrogate model based on ANN is used to predict the total pressure loss coefficients, outlet flow angles, and laminar flow areas of the deformed blade profiles. One separate neural network is used for prediction of each aerodynamic parameter. During the training process of the DDPG, the surrogate model acts as the role of environment, giving the predicted aerodynamic parameters and calculating the reward of the given action and state.

The structure of the neural network used in this paper is shown in Figure 12. The input variables are the geometrical parameters, and the output variable is the aerodynamic parameter. The network has a six-layer structure with four hidden layers. The numbers of nodes of the four hidden layers are $60,100,80,40$, respectively. The selection of the ANN structure, including the number of hidden layers and the number of nodes in each hidden layer, is experience-based [57]. In this paper, referring to Reference [34], the structure is defined by adjusting and trying many times. The activation function of each layer is ReLU (Rectified Linear Unit) function $f\left(x\right)$, in which the expression can be written as: $f\left(x\right)=\max (0,x)$. The data is divided into a training set $\left(1000\right)$ and test set $\left(200\right)$. The optimizer of the network is Adam [58]. The learning rate is set as $0.03$ to make sure the network converges fast.

Since the transition location and the total pressure loss vary consistently in a small range near the design point when the modification is reasonable and small, while the outlet flow angles vary less consistently over the operating range, the total pressure loss coefficient and the laminar flow areas at the operating point and multiple on-design and off-design outlet flow angles are considered in the optimization. One network with the above structure and hyperparameters is built and trained for each output variable in aerodynamic parameters ($\overline{\omega }$, ${\tilde{A}}_{SS}$, ${\tilde{A}}_{PS}$, ${\beta }_2$ of ${\beta }_1$ in (${44.66}^{\circ },{45.66}^{\circ },{46.66}^{\circ },{47.66}^{\circ },{48.66}^{\circ },{49.66}^{\circ }$)). The trained networks can represent the correspondence between geometric data and aerodynamic data. The aerodynamic parameters of one specific deformed rotor can be predicted easily by feeding the geometric parameters to the trained networks.

5. DDPG-Based Optimization for Blade Profile

5.1. Deep Deterministic Policy Gradient (DDPG)

In this paper, DDPG [32] is used as the intelligent designer to find the optimial policy for the multi-objective design of cascade blade. DDPG is one reinforcement learning algorithm which can choose the appropriate action according to specific state. Different from the value-based learning methods, e.g., deep Q-learning (DQN), which is restricted to discrete and low-dimensional state and action space, DDPG can handle the high-dimensional complicated action problems with the continuity of state and action space.

DDPG is developed from DPG, by using deterministic policy instead of random policy [31]. Since the policy is deterministic, the output action is determined according to the input state and does not need to be sampled from a probability distribution. Therefore, the calculation of internal expectation for policies is avoided, and external expectation only needs to be calculated according to the environment. The computational resources are saved, and the variance becomes small, improving the convergence.

DDPG is a policy gradient method with an actor-critic structure [59]. The actor network approximates to the deterministic policy function $a=\mu \left(x\right)$, which learns a policy that yields the highest possible reward so that a definite action from the current state can be given. The critic network approximates the action-state value function $Q^{\mu }(s,a)$ to criticize the given action, according to the action $a_t$ in state $s_t$ and following policy $\mu$. With the network parameter ${\mathsf{\Theta }}^Q$, the critic network can be written as $Q(s,a|{\mathsf{\Theta }}^Q)$. The actor network can be written as $\mu \left(s\right|{\mathsf{\Theta }}^{\mu })$ with parameter ${\mathsf{\Theta }}^{\mu }$. With the value function, the actor network can be updated by the gradient of the policy $J\left({\mathsf{\Theta }}^{\mu }\right)$, according to the deterministic policy gradient theorem:

$${\nabla }^{2}J\left({\mathsf{\Theta }}^{\mu }\right)=E_{s\sim {\rho }^{\mu }}\left[{\nabla }_{\mathsf{\Theta }\mu }\mu \left(s\right|{\mathsf{\Theta }}^{\mu }){\nabla }_{a}Q(s,a|{\mathsf{\Theta }}^Q){|}_{a={\mu }_{\mathsf{\Theta }}\left(s\right)}\right].$$

(12)

Drawing on the idea of DQN, DDPG has two frameworks of actor-critic structure: main framework and target framework. Each framework contains one actor network and one critic network. The main actor network gives the action $a_t$. The main critic network outputs the value function $Q(s_t,a_t)$ and updates the actor network. The action $a_t$ is performed in the environment, to achieve the next moment state $s_{t+1}$ and the calculated reward $r_t$. The target framework stores the fixed target networks. The target actor network gives $a_{t+1}$ according to the state $s_{t+1}$ and the policy ${\mu }^{\prime }\left(s_{t+1}\right|{\mathsf{\Theta }}^{{\mu }^{\prime }})$. The target critic network outputs the value $Q^{\prime }(s_{t+1},a_{t+1})$. When the target networks are updated, soft update is used: ${\mathsf{\Theta }}^{Q^{\prime }}\leftarrow \tau {\mathsf{\Theta }}^Q+(1-\tau ){\mathsf{\Theta }}^{Q^{\prime }}$,${\mathsf{\Theta }}^{{\mu }^{\prime }}\leftarrow \tau {\mathsf{\Theta }}^{\mu }+(1-\tau ){\mathsf{\Theta }}^{{\mu }^{\prime }}$, where $\tau$ is the soft-update parameter. Therefore, the target network parameters can be updated slowly, with improved learning stability. The main critic network is updated to minimize the loss function:

$$L=E\left[{(r_t+\gamma Q(s_{t+1},a_{t+1})-Q(s_t,a_t))}^2\right].$$

(13)

DDPG uses an experience replay memory to store transitions $(s_t,a_t,r_t,s_{t+1})$ from interacting with the environment. The experience transitions are stored in the memory and then fed to update the networks in batches with a random sample strategy. Random sampling can eliminate the correlation between transitions, improving the stability of network training. The structure of DDPG is shown in Figure 13.

Based on DPG, DDPG uses deep neural networks. A multilayer feed-forward network with enough neurons in hidden layers can approximate any continuous function with arbitrary accuracy [60]. With the deeper network structure, the actor and the critic network can approximate the policy and the value function more accurately. The deeper network structure also increases the dimension of action that can be handled. But, for problems with complex action space, increasing the number of hidden layers and the number of nodes of each layer may increases the training time of the network, and even lead to the tendency of over-fitting.

5.2. Application of DDPG and Optimization

DDPG is used in this paper to approach optimization policies. Since the optimization is to achieve the blade profile with the optimal aerodynamic parameters by modifying the geometry, it is suitable to take the geometric parameters space as the state of DDPG network, and to choose the variations of geometric parameters as the action.

The hyperparameters of the DDPG network are as follows: The target smoothing coefficient is 0.005. The discounted factor is 0.99. The batch size equals 64, and the replay buffer size equals 50,000. The learning rate is selected as 0.01, to make sure the network can converge faster. The other hyperparameters are selected by referring to Reference [61,62] with little adjustment, which is common and widely-used in the applications of DDPG. The hyperparameters are selected in order to ensure the actor and critic networks can converge and the DDPG network can achieve high rewards.

For multi-objective optimization, the analysis of objective space can help us better understand the relationship between optimization variables, which is directly associated with a more problem-oriented formulation of optimization objective. The aerodynamic parameter database generated to train the network can provide a large number of samples for objective space analysis. Principal Component Analysis (PCA) [63] is an effective tool to extract the main modes from objective space considered to be the essential information that may help establish a more reasonable optimization objective.

For blade profiles in the profile library sampled corresponding to given n, the aerodynamic parameters of n profiles achieved through CFD simulation constitute a database represented in Equation (14), where the eight columns are the total pressure loss, relative laminar flow areas, outlet flow angles of different inlet flow angles, respectively. To compare the parameters according to the principal component, the aerodynamic parameters in $\mathbf{A}$ are normalized. To obtain the relative distribution matrix $\tilde{\mathbf{A}}$, Each column of the data matrix $\mathbf{A}$ is subtracted from the mean value of the column. According to Equation (15), the data from objective space systematically forms the covariance matrix $\mathbf{P}$. The eigendecomposition is shown in Equation (16), where $\mathsf{\Lambda }$ represents diagonal matrix of eigenvalues and $\mathbf{U}$ represents agglomeration of eigenvectors.

$$\mathbf{A}=\left[\begin{array}{ccccccccc}{\overline{\omega }}^{\left(1\right)}& {\tilde{A}}_{SS}^{\left(1\right)}& {\tilde{A}}_{PS}^{\left(1\right)}& {\beta }_{21}^{\left(1\right)}& {\beta }_{22}^{\left(1\right)}& {\beta }_{23}^{\left(1\right)}& {\beta }_{24}^{\left(1\right)}& {\beta }_{25}^{\left(1\right)}& {\beta }_{26}^{\left(1\right)}\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ {\overline{\omega }}^{\left(n\right)}& {\tilde{A}}_{SS}^{\left(n\right)}& {\tilde{A}}_{PS}^{\left(n\right)}& {\beta }_{21}^{\left(n\right)}& {\beta }_{22}^{\left(n\right)}& {\beta }_{23}^{\left(n\right)}& {\beta }_{24}^{\left(n\right)}& {\beta }_{25}^{\left(n\right)}& {\beta }_{26}^{\left(n\right)}\end{array}\right],$$

(14)

$${\mathbf{P}}_{n\times n}=\frac{1}{n}{\tilde{\mathbf{A}}}^T\tilde{\mathbf{A}},$$

(15)

$${\mathbf{P}}_{n\times n}=\mathbf{U}\mathsf{\Lambda }{\mathbf{U}}^T.$$

(16)

Through the analysis of eigenvectors, the relationship between the parameters in the multi-objective optimization of the cascade blade can be known and help to determine the reasonable reward function of reinforcement learning method and appropriate objective function of optimization. The first and second principal component of the aerodynamic data are listed in

Table 4. The 1st and 2nd principal component of objective space analysis.

	$\overline{\mathit{\omega }}$	${\tilde{\mathit{A}}}_{\mathit{SS}}$	${\tilde{\mathit{A}}}_{\mathit{PS}}$	${\mathit{\beta }}_{21}$	${\mathit{\beta }}_{22}$	${\mathit{\beta }}_{23}$	${\mathit{\beta }}_{24}$	${\mathit{\beta }}_{25}$	${\mathit{\beta }}_{26}$
1st PC	−3.4 $\times {10}^{-16}$	7.6 $\times {10}^{-2}$	−1.3 $\times {10}^{-1}$	5.6 $\times {10}^{-1}$	−5.3 $\times {10}^{-2}$	1.5 $\times {10}^{-1}$	2.8 $\times {10}^{-1}$	3.8 $\times {10}^{-1}$	−1.6 $\times {10}^{-1}$
2nd PC	−4.1 $\times {10}^{-1}$	5.7 $\times {10}^{-1}$	3.9 $\times {10}^{-2}$	−6.9 $\times {10}^{-2}$	−1.2 $\times {10}^{-1}$	−2.4 $\times {10}^{-1}$	−3.0 $\times {10}^{-1}$	−3.0 $\times {10}^{-2}$	6.8 $\times {10}^{-2}$

. The first principal component represents the internal relation of outlet flow angle variations. The second principal component contains information: the decrease of total pressure loss coefficient has a strong dependence on the increase of laminar flow area. It may be inferred that the increase of laminar flow area is beneficial to the reduction of total pressure loss, and the influence of the laminar flow area of the suction surface is much greater than that of the pressure surface.

The environment of the DDPG network in this paper is the ANN-based surrogate model. With the given state $s_t$ and action $a_t$ (the profile geometry and geometric change), the next step state $s_{t+1}$ (updated geometry) can be obtained. The aerodynamic parameters according to the state and action can be predicted with the surrogate model, to calculate the reward function $r_t$. The reward function represents our actual evaluation of an action. It is used to adjust the critic network to output value function closer to reality. Since a blade profile is expected with lower total pressure loss coefficient, larger laminar flow area, and as small as possible change in outlet flow angle, the reward function can be expressed as:

$$\begin{array}{cc}\hfill r=& ({\overline{\omega }}^0-\overline{\omega })\ast w_0-\max (0,({\tilde{A}}_{SS}^0\ast 1.15-{\tilde{A}}_{SS}))\ast w_1\hfill \\ \hfill & -\max (0,({\tilde{A}}_{PS}^0\ast 1.15-{\tilde{A}}_{PS}))\ast w_2-\sum _i{\varphi }_i\ast w_3,\hfill \end{array}$$

(17)

where $\varphi =\max (0,{\beta }_2-{\beta }_2^{+})+\max (0,{\beta }_2^{-}-{\beta }_2)$, and ${\beta }_2^{+}$ and ${\beta }_2^{-}$ are boundaries to constrain ${\beta }_2$ of different ${\beta }_1$. ${\overline{\omega }}^0$, ${\tilde{A}}_{SS}^0$, and ${\tilde{A}}_{PS}^0$ represent the initial profile’s $\overline{\omega }$, ${\tilde{A}}_{SS}$, and ${\tilde{A}}_{PS}$. According to the PCA, the weights in Equation (17) are: $w_0=100,w_1=3,w_2=1.5,w_3=1/3$.

The target of the optimization is to achieve a blade profile with reduced total pressure loss and improved laminar flow area. During the training process of the DDPG network, huge numbers of states (profiles) interact with the environment, many of which have good performance. Each time the reward value of one state is calculated, the aerodynamic parameters of the blade profile of the state are predicted by the surrogate models. The state (the geometric parameters) of the blade profile with the smaller objective function is stored in the training process of the DDPG network to optimize the blade profile. Compared with the initial blade profile, the optimized blade profile should have as small total pressure loss coefficient as possible, enough large laminar flow area, and enough small change in outlet flow angle, so the objective function of the optimization can be expressed as:

$$\begin{array}{cc}\hfill loss=& \overline{\omega }\ast w_0+\max (0,({\tilde{A}}_{SS}^0\ast 1.15-{\tilde{A}}_{SS}))\ast w_1\hfill \\ \hfill & +\max (0,({\tilde{A}}_{PS}^0\ast 1.15-{\tilde{A}}_{PS}))\ast w_2,\hfill \end{array}$$

(18)

with a constraint on ${\beta }_2$. The weights are: $w_0=1,w_1=0.15,w_2=0.075$.

With the optimization policy from DDPG, the optimizer can approach the good-performance profile (high reward state) with not many iterations of geometric modification (action).

6. Results and Discussion

The convergence curves of the neural networks are shown in Figure 14. The aerodynamic parameters, except ${\tilde{A}}_{SS}$, have a good convergence result with MSEs (mean square errors) by less than ${10}^{-5}$. ${\tilde{A}}_{SS}$ varies much for profiles in the profile library; thus, it does not converge as well as other parameters, but it is still within small MSE of ${10}^{-4}$. The small convergence errors show that the ANNs have good performance of prediction on the test sets. The independence study of division of train-set and test-set is shown in Appendix A.2. Figure 15 shows the reward curve of the DDPG algorithm and the objective function curve of the optimization. The DDPG network can be trained well and obtain a high reward in 100 episodes. With the trained DDPG network, the optimizer can find and record a profile with very good performance in no more than 200 episodes, which, as it turns out, means that this optimization method based on DDPG can rapidly achieve a fairly good result. The better-performance profile updates slow and varies slightly after 200 episodes, as well as the objective function.

With the ANN-based surrogate model and the DDPG based optimization, a new rotor blade profile with lower total pressure loss coefficient and larger laminar flow area is achieved. By CFD validation, the performance characteristics and the flow field characteristics of the cascade with the optimized blade profile can be known and compared with the original cascade.

Figure 16 shows the original and the optimized rotor blade profiles. Compared to the initial profile, the optimized profile has following characteristics: (1) The thickness of the front of the profile is slightly reduced; (2) the thickness of the middle part of the profile is increased, where the thickness at the pressure surface is increased slightly resulting in a flatter shape, and the thickness and curvature at the suction surface is obviously increased.

The aerodynamic performance curves of the cascades of the original blade profile and the optimized profile are shown in Figure 17, respectively. In the operating region, with ${\beta }_1$ being ${45.16}^{\circ }$ to ${48.66}^{\circ }$, $\overline{\omega }$ of the cascade with the optimized blade profile decreases considerably, with a relative reduction up to $6.07\%$. The variation of ${\beta }_2$ is small, and the relative change in the operating region is no more than $0.31\%$.

The performance characteristics of the cascades at the operating point are listed in

Table 5. Performance characteristics of the cascades at the operating point.

	$\overline{\mathit{\omega }}$	${\tilde{\mathit{A}}}_{\mathit{SS}}$	${\tilde{\mathit{A}}}_{\mathit{PS}}$
initial	$0.01642$	$34.56\%$	$29.26\%$
optimized	$0.01583$	$43.35\%$	$30.05\%$

. $\overline{\omega }$ of the cascade goes from $0.01642$ to $0.01583$, with a relative reduction of $3.59\%$. The relative laminar flow area at the suction surface ${\tilde{A}}_{SS}$ rises from $34.56\%$ to $43.35\%$, with massive increase. The relative laminar flow area at the pressure surface is basically unchanged. It can be concluded from the results that the DDPG-based optimizer successfully optimized the cascade blade profile with $\overline{\omega }$$3.59\%$ relatively being reduced, and ${\tilde{A}}_{SS}$$25.4\%$ higher than the initial profile. The large improvement of ${\tilde{A}}_{SS}$ is in line with the information from the PCA that ${\tilde{A}}_{SS}$ contributes more on reducing $\overline{\omega }$ than ${\tilde{A}}_{PS}$ (cf. Table 4). CFD simulation and ANN prediction on performance characteristics of the initial and optimized cascades are compared in

Table 6. Comparison of CFD simulation and ANN prediction on performance of the initial and optimized cascades.

	${\overline{\mathit{\omega }}}^{\left(\mathit{i}\right)}$	${\tilde{\mathit{A}}}_{\mathit{SS}}^{\left(\mathit{i}\right)}$	${\tilde{\mathit{A}}}_{\mathit{PS}}^{\left(\mathit{i}\right)}$	${\overline{\mathit{\omega }}}^{\left(\mathit{o}\right)}$	${\tilde{\mathit{A}}}_{\mathit{SS}}^{\left(\mathit{o}\right)}$	${\tilde{\mathit{A}}}_{\mathit{PS}}^{\left(\mathit{o}\right)}$
CFD simulation	$0.01642$	$34.56\%$	$29.26\%$	$0.01583$	$43.35\%$	$30.05\%$
ANN prediction	$0.01641$	$35.40\%$	$29.00\%$	$0.01567$	$41.84\%$	$29.58\%$

. Compared to the CFD results, the ANN prediction has an error of about $0.8\%$ on the total pressure loss coefficients and an error of about $3.5\%$ on the laminar flow areas, which proves that the ANN-based model can accurately predict the aerodynamic performance

Figure 18 shows $C_f$ distributions on suction and pressure surfaces. Compared with the initial profile, the position of the jump of $C_f$ is moved aft at the suction surface, which indicates an obvious delay of transition. The $C_f$ distributions on pressure surfaces are basically the same, consistent with the laminar flow areas obtained with the intermittency factor distributions, as in Table 5.

The absolute Mach number contours of the original and the optimized blades are plotted in Figure 19. Compared with the initial profile, the shock wave at the suction surface is weakened and moved aft. This results in a reduction in shock loss, which is also a contribution to the decrease of the total pressure loss.

7. Conclusions

In this paper, a multi-objective optimization method based on reinforcement learning technique is proposed and applied in a hybrid optimization of the compressor cascade blade profile on the total pressure loss and the laminar flow area. Firstly, an ANN-based surrogate model is used as the environment in the DDPG network, feeding back the aerodynamic parameters corresponding to the deformed profiles with surrogate model trained with parameterized profiles and aerodynamic data calculated via CFD solver. Secondly, based on the reinforcement learning technique of DDPG, a multi-objective optimization method is proposed, where the profile geometry is taken as the State of the DDPG network, the modification of geometry is taken as the Action, and the DDPG network acts as an intelligent optimizer learning the design experience from the interaction with the environment. The optimizer can approach the good-performance profiles rapidly according to the learnt design experience.

With the optimization on total pressure loss and laminar flow area, the cascade with modified blade profile has an improved performance and increased laminar flow area. The total pressure loss coefficient of the optimized profile is reduced by $3.59\%$, and the relative laminar flow area at the suction surface is $25.4\%$ higher than the initial profile. The result shows that the DDPG-based method is one successful multi-objective optimization of cascade blade profile, with an improvement on using the accumulated calculation data to extract design experience and making design policies for guiding optimization. The cascade blade profile optimized in this case is the section profile at $40\%$ span of one compressor rotor blade. This optimization method can be applied to profiles at different spans of the rotor. With multiple profiles at various spans considered, this method can be extended to 3D shape optimization method and be applied to 3D rotor blades.