Reverse Design of Solid Propellant Grain for a Performance-Matching Goal: Shape Optimization via Evolutionary Neural Network

: The reverse design of solid propellant grain for a performance-matching goal, one of the most challenging directions of the solid rocket motor designing work, is limited by the traditional semi-empirical parameter-driven optimization methods based on some predeﬁned grain conﬁgurations. Grain designers call for a new method that can automatically provide brand-new grain shapes beyond the traditional ones. In this work, a shape optimization method based on the evolutionary neural network is proposed to achieve the reverse design of two-dimensional (2D) grains. Firstly, the modiﬁed ellipse-form eikonal equation is solved by using the ﬁnite element method to realize the burn-back analysis of 2D grains in any shape on a ﬁxed unstructured mesh. Then, the neural network is introduced to determine the spatial distribution of the propellant to deﬁne the grain shape. The hyperparameters of the network are continuously evolved with the aid of the genetic algorithm. Finally, the optimal grain shape that matches the performance goal most is obtained. The method is veriﬁed in different scenarios. The result shows that the design can precisely match the given pressure-time curve of star grains and slotted-tube grains. Furthermore, the method can automatically evolve a new dendritic-shaped grain that matches the given dual-thrust pressure-time curve. Since the reverse design uses the concept of shape optimization, it does not require any pre-selection of the grain shape, and the designers shall be free from deﬁning different kinds of geometric parameters for speciﬁc grain conﬁgurations. Consequently, the method has the potential to apply in the reconstruction of an actual grain and the conceptual design of innovative grain conﬁgurations.


Introduction
Propellant grain design is a significant stage of solid rocket motor designing work and its major task is to match the given performance curves. The main internal ballistic performance curves include the variations of combustion chamber pressure with time (p c~t curve) and motor thrust with time (F~t curve). These curves are essential to the overall design of the missile/rocket and even the guidance, navigation, and control (GNC) system. Performance prediction, commonly used in grain design, calculates the performance curve according to given grain configurations. However, this design method often falls into a trialand-error approach and fails to obtain the best solution matching the given performance curve. In contrast, the reverse design, a concept rarely mentioned by previous scholars, is an inversed problem of the performance prediction. It seeks the optimal grain shape according to the given performance curve [1]. The relation between the performance prediction and grain reverse design is shown in Figure 1. The reverse design is an optimization problem of grain configuration essentially. It is always a nonconvex nonlinear problem, and the objective function does not appear to have a clear analytical formula. Therefore, it is hard to find the optimal global solution via traditional line search methods or trust-region methods. Sforzini [2,3] is the first to study reverse design. The developed solid rocket motor design and optimization (SRMDO) code uses the pattern search method to optimize the 14 parameters of a 3D grain. The objective function is to minimize the square sum of the gap between the designed thrust curve and the given one. Since the 21st century, with the emergence and broad application of heuristic optimization algorithms such as the particle swarm optimization (PSO) and the genetic algorithm (GA), it has become possible to solve the reverse design practically. Albarado [4] uses a hybrid optimizer combining pattern search and particle swarm to minimize the square sum of the gap between the designed pressure curve and the given one. Yücel [5] uses the GA to optimize a finocyl grain to ensure that the error of the thrust curve is minimum. Recently, the reverse design effort is also denoted as the performance matching by Wu [6] or predefined performance criteria by Hashish [7]. So far, the PSO and GA has become a common method for grain optimization design [8][9][10].
Although researchers have developed all kinds of procedures to achieve reverse design, these procedures all share the same pattern. Firstly, the geometry of the grain is strictly determined by several parameters. Then, the selection of the grain shape, such as the star, wagon wheel, and finocyl shape, is supposed to be consistent with historical experience. After that, the global optimization algorithms are conducted to optimize the geometric parameters so that the designed performance curve is as coherent as possible with the goal. If it fails to find the optimal solution, reselection of the grain shape and redesign are necessary. Consequently, these procedures are semi-empirical and parameter-driven.
The traditional semi-empirical parameter-driven optimization methods have significant shortcomings. Firstly, the design results are limited by the existing grain shapes and historical experience. Although there have been attempts to use expert systems to select grain configurations [11], it is still difficult to design new grain shapes with higher loading fraction and higher performance. Secondly, to prevent interference between the geometric features of the grain, there are numerous nonlinear constraints between the geometric parameters. It is challenging to fully understand these constraints and add them to the optimization algorithm. In addition, the number of optimized variables is equal to the number of geometric parameters, so the degree of freedom of the design is relatively low, making it frequently fail to meet the given performance curve. Finally, the geometric parameters corresponding to different grain shapes are different, so it is not easy to develop a unified optimization framework for multiple grain shapes.
To overcome the shortcomings of the traditional methods of the reverse design, it is necessary to introduce the shape optimization technique. Shape optimization, also considered as a special topology optimization, has been widely used in the mechanical structure design. It is dedicated to finding the optimal distribution of materials instead of optimal geometric parameters [12][13][14]. The materials distribution, that is, the propellant distribution, can be defined by the phase field. Conventionally, the phase field function can be expressed by basis functions, such as finite element basis functions, radial basis functions, and spectral parameterization [13]. With the help of machine learning, a simple feedforward neural network (FNN) is also capable of expressing the phase field. Once FNN parameterizes the phase field, one can easily define and manipulate an irregular grain.
Unlike mechanical structure optimization, the challenging part of grain shape optimization lies in the burn-back analysis of the irregular grain shape defined by the phase field. As the proceeding of combustion, the burning surface of the grain changes, which is known as surface regression. The motor's internal ballistic characteristics can only be obtained by accurately calculating the burning surface area. There are two types of burnback analysis methods: the non-meshing and meshing methods. Non-meshing methods include the analytical method [15,16], general coordinate method [17,18], and CAD modeling method [19,20]. Meshing methods include the minimum distance method [21,22], moving grid method [23,24], level-set method [9,25,26], and fast marching method [27]. On the one hand, the non-meshing methods are accurate, computationally efficient, and reliable. However, they can hardly adapt to the grains with irregular geometry or complex distribution of burning rate. Thus, the non-meshing methods are unfit to apply in shape optimization tasks. On the other hand, the meshing methods possess great generality, but they need a lot of pre-processing and post-processing work. Therefore, they are not yet widely used in engineering. To achieve shape optimization, designers urgently need a new meshing method that is easy to implement on any irregular grain shape-this will be one of the focuses in this work.
This work proposes the shape optimization method based on the evolutionary neural network to achieve the reverse design of two-dimensional (2D) grain. Given the pressuretime curve, the GA is used to optimize the hyperparameters of the neural network which can fully define the grain shape to minimize the gap between the designed results and the goal. As the method uses the concept of shape optimization, it does not need to pre-select any grain shape, and the designers can be free from defining different kinds of geometric parameters for specific grain shapes. Furthermore, the degree of freedom in the evolutionary neural network is controllable, so it will obtain much more optimal and suboptimal solutions than the traditional methods. The method has the potential to apply in many new application scenarios, such as the conceptual design of creative grain shapes, the reconstruction of an actual grain, and the shape design of the additive manufactured grain.
The structure of this work is as follows. The second chapter begins with the steps of grain reverse design. It will continue with a modified PEF method to realize the gas-solid unified burn-back analysis based on the phase field. Then, it will discuss the definition of the phase field by introducing the neural network. Finally, the chapter ends with the genetic algorithm framework of reverse design. The third chapter will provide some test results of the method in detail, and the fourth chapter will give the conclusions.

Steps of the Grain Reverse Design
The performance prediction of solid rocket motor can be divided into two steps in practice, as shown on the left side of Figure 2. The first step is the burn-back analysis and the second step is the internal ballistic calculation. Similarly, the grain reverse design, essentially an inversed problem of the performance prediction, can also be divided into two steps, as shown on the right side of Figure 2. The first step, denoted as reverse internal ballistic calculation, converts the performance curve (p c~t curve) into variation of burning surface area with burned web thickness (A b~w curve), and the given variables are listed in Table 1. The second step, denoted as grain reconstruction, is an optimization process that finds the grain shape matches the A b~w curve most. This division strategy decomposes the reverse design, which can effectively reduce the complexity.  To solve the reverse internal ballistic calculation, the following three assumptions are proposed: Assumption 1: The nozzle throat to grain port area ratio of the motor is small enough so that the gas flow speed is low, and there is no erosive combustion phenomenon.
Assumption 2: The combustion chamber pressure does not change drastically during the whole working process, which makes the equilibrium pressure formula reasonable.
Assumption 3: The characteristic velocity and density of the propellant do not vary with the combustion chamber pressure.
For brevity, the concept of similar grains is proposed. As shown in Figure 3, two 2D grains are admitted being similar if all their sides are proportional. The dimensionless burn perimeterl b is defined as Equation (1) where R is the outer radius of the grain. For 2D grains, the burn perimeter l b is defined as Equation (2): where A b is the burning surface area of the grain, and L is the length of the grain. It is worth noting that the grain shape can often be divided into the several minor geometry units that can repeat themselves by rotational and mirror transformations. Therefore, concentrating on one single unit provides much more efficiency than the whole grain. The number of units (half an angle) is denoted as N. For example, in Figure 3, one geometric unit has half a slot, yielding N = 12. The dimensionless burned web thicknessŵ is also defined as Equation (3): For 2D grains, the reverse internal ballistic calculation needs to convert the known p c~t curve into the unknownl b ∼ŵ curve. As shown in Figure 3, although the grain in the physical space may be much smaller (or greater) than the one in the cyber space, they are still similar grains and sharing the samel b ∼ŵ curve.
The advantage of nondimensionalization is that it can significantly simplify the second step of reverse design, denoted as grain reconstruction. Any grain in physical space can be transformed into a standard grain (outer radius R 0 is 1 m) in the cyber space through scaling transformation, so the burn-back analysis can be performed within a unified standard framework. The designers only need to focus on optimizing the internal profile without considering a mutable outer radius. As soon as the grain in the cyber space meets the requirements, it can be directly returned to the physical space through simple scaling transformation.
The calculation formula of the burn rate r b can usually be expressed as Equation (4) [28]: where a and n are the burn rate coefficient and pressure index, respectively. Derived from the steady-state continuity equation, the equilibrium pressure of the combustion chamber is Equation (5): where ρ p , c * and A t are the propellant density, characteristic velocity, and nozzle throat area, respectively. According to Equations (4) and (5), the expression forl b andŵ over time is as Equation (6): where t a is the action time, and the subscript T in the formula represents the design target, which is the requirement of designers. The C T is a dimensionless variable expressed as Equation (7): Since the designers are uncertain about the L and A t from the beginning, the C T is also a quantity that needs to be designed.l b,T (t)C T is the product ofl b,T (t) and C T . It is an independent variable and can also be regarded as a dimensionless burn perimeter. It is important to emphasize that N is unknown in practice, but to reduce the computation cost, N is considered to be a given integer in this work. Since N is a finite integer, the optimal N can be found by enumeration.
Given the variables in Table 1, the l b,T C T ∼ŵ T curve will be obtained via Equation (6). So far, the reverse internal ballistic calculation is completed.

Step2: Grain Reconstruction
Given the l b,T C T ∼ŵ T curve, grain reconstruction refers to finding a 2D grain that satisfies this curve. The following three requirements are proposed: Requirement 1: The burning starts from the internal port surface of the grain. This kind of grain is known as internal burning grain.

Requirement 2:
The solid region for one geometric unit shall be connected without holes.

Requirement 3:
The gas-solid interface can be an irregular curve.
Using the subscript D to represent the design result, the grain reconstruction is to optimize the internal profile in the cyber space so that thel b,D ∼ŵ D curve obtained by the burn-back analysis matches thel b,T ∼ŵ T curve. In practice, a method based on linear regression can evaluate the matching degree quantitatively according to thel b,D ∼ŵ D curve and l b,T C T ∼ŵ T curve calculated from Equation (6)    So far, the procedure of the reverse design is clear. The specific methods will be described further in the next sections.

PEF Method (Poisson Equation-Eikonal Equation-Finite Element Method)
Since the method for burn-back analysis in this work uses the Poisson equation, the eikonal equation, and the finite element method, it is called the PEF method. Similar to the propagation of light [29], burning surface regression of complex grains can be described by the eikonal equation. In the grain with uniform burning rate, it can be written as Equation (8): where W, in meters, is the field of burned web thickness w, which gives the thickness of the burned web corresponding to any spatial coordinate. Any isosurface of the W field represents the burning surface. The meaning of Equation (8) is that the gradient of the W field is 1 everywhere, which is consistent with the assumption of parallel-layer combustion widely adopted in solid rocket motor design. There are two kinds of boundary conditions, that is, the Dirichlet boundary condition W = 0 at the initial burning surface and the zero-flux boundary condition ∂W/∂n = 0 at the flame-retardant surface or symmetry surface. However, for the grain with the non-convex gas cavity, as the burning surface moves, the W field is no longer derivable at the intersection of burning surfaces (see Figure 5). One possible way to ensure that the burning surface keeps smooth numerically is to introduce a burning rate that varies with the curvature. The essence is to add a diffusion term based on Equation (8). Then, Equation (9) is obtained [30]: where α is the diffusion term coefficient, a positive number in meters. If the diffusion term coefficient is too large, the calculation accuracy will decrease. On the other hand, if it is too small, the convergence of the finite element method will deteriorate. The recommended value range is: where δl is the local element size. Figure 6 shows the results of burning surface regression for a typical star grain (R = 0.5 m, N = 12) with different meshes. Figure 7 shows that with the refinement of the element size, compared with the analytic solution, the accuracy of PEF method will continue improving. Compared to the traditional method of burn-back analysis, the advantages of the PEF method are significant. It can perform the burn-back analysis of 3D grain with irregular geometry robustly. In addition, since the PEF method converts the burn-back analysis into an analogous steady heat conduction problem, it can be directly modeled and solved using COMSOL, ANSYS, ABAQUS, and other commercial finite element software, which are relatively easy to implement.

Modified PEF Method (ϕ-PEF)
It is worth noting that the solution domain of the PEF method is only the solid-phase region, excluding the gas-phase region. Therefore, the gas-solid unified burn-back analysis cannot be accomplished. Since the exact shape of the solid-phase region is unsure in the current reverse design, the PEF method must be modified.
The gas-solid unified burn-back analysis can be established by introducing the virtual burning rate and the source term factor. Rewrite the governing equation as Equation (11): where r is the dimensionless virtual burning rate, and s is the dimensionless source term factor. Inspired by the level-set method [31] and phase field method [32] of structural topology optimization, the concept of phase field is then introduced to clearly describe the gas-phase region, solid-phase region, and their interfaces. A dimensionless scalar ϕ(x, y) is defined as Equation (12): The distribution of ϕ(x, y) in the whole domain constitutes the phase field. As long as this distribution is obtained, it is equivalent to obtain a 2D grain shape. The initial burning surface is given by the zero level ϕ(x, y) = 0.
Based on the phase field, the virtual burning rate field r can be constructed, subsequently, which is defined as Equation (13): where K is a constant much larger than 1. The physical meaning of the virtual burning rate r is that the burning rate in the gas region is so large that the web thickness field W is approximately equal at any position. Unfortunately, when the virtual burning rate r presents in the gas-phase region, the source term near the gas-solid interface has a great jump. The iteration of finite element method (FEM) is difficult to converge due to its nonlinearity. Therefore, the source term factor s is then introduced, which is defined as Equation (14): The physical meaning of the source term factor is as follows. In the solid-phase region, the original governing equation is still maintained. While in the gas-phase region, the source term of the governing equation is set to 0, and the governing equation degenerates into a steady-state heat conduction equation without any internal heat source, also known as the Laplace equation. Solving the Laplace equation in the gas-phase region will significantly improve the convergence of the FEM.
Since the method uses phase field, it is called the ϕ-PEF method in this work. It has many advantages. Firstly, the burn-back analysis can be completed by using a fixed mesh without any extra work to track the moving interface. It avoids remeshing after changing the grain shape, which is quite useful in the shape optimization. Secondly, since the mesh is fixed, setting up and defining Dirichlet boundary conditions are simplified. Finally, inheriting the advantages of the PEF method, it can be directly modeled and solved in the commercial finite element software, which is relatively easy to implement.
To compare the difference between the PEF method and ϕ-PEF method, the calculation results for the same grain are illustrated in Figure 8. The solid-region of the grain is formed by the intersection of a square with a side length of 2 m and a circle with a radius of √ 5/2 m. The workflow of the PEF method is as follows. Firstly, the geometry is drawn using CAD software. The right-hand-side boundary condition is set to the flame-retardant layer and the rest of the boundary are set to the Dirichlet boundary condition. Then with the help of unstructured mesh, Equation (9) is discretized into an algebraic equation to solve the W field, while, the ϕ-PEF method uses the phase field to describe the shape of the solid region of the grain in a fixed square region, and the phase field for this problem has an analytical expression as Equation (15): The W field is obtained by solving the algebraic equation derived from Equation (11). It can be seen from Figure 8 that the results obtained by both methods are identical.
Once the W field is obtained, the burning surfaces can be obtained by extracting the isosurfaces of the W field, and the burn perimeter l b can be calculated from these burning surfaces. Since the ϕ-PEF method is much more flexible, we chose it to achieve the burn-back analysis of the irregular grain shape. In Section 2.2, the grain design problem is transformed into the design work of a continuous two-dimensional function ϕ(x, y), and this function will become the optimization variable. The corresponding optimal solution can be obtained by using indirect methods (IM) such as optimal control theory [33]. Although the IM can effectively guarantee optimality and high accuracy, there are still great difficulties in deriving the first-order optimal conditions and solving two-point boundary value problems [34]. Therefore, this work will focus on the more practical method denoted as the direct method (DM) to solve the optimization problem.
The DM, approximating the function to be solved as a nonlinear combination of finite parameters, transforms the design problem of the function into a multi-parameter optimization problem [35]. Finally, the infinite-dimensional optimization problem is transformed into a finite-dimensional one, as shown in the Equation (16): where p is the independent vector of the optimization problem,φ is a predetermined function and is often selected as a linear combination of basis functions. However, its accuracy is generally insufficient especially for nonlinear problems, which cannot meet the demands for an engineering application, so the neural network is proposed next.

Design of Feedforward Neural Network
Artificial neural network (ANN) refers to a series of mathematical models inspired by biology and neuroscience. These models can simulate the biological neural network by mathematically abstracting the biological nervous system. After decades of development, ANN, the backbone of artificial intelligence technology, has been widely used in mathematical and engineering fields such as fitting, classification, and clustering.
Feedforward neural network (FNN) is the first ANN invented. In each layer, neuron units can receive signals from units in the previous layer and generate signals to transmit to the next layer. According to the Universal Approximation Theorem, for an FNN with a linear output layer and at least one hidden layer using a nonlinear activation function, as long as the number of neurons in the hidden layer is large enough, the FNN can approximate any bounded closed set function defined in the real number space with any accuracy [36,37].
After considering the computational cost and fitting accuracy, we designed an FNN of one single hidden layer to construct the functionφ(p, x, y). The input layer is denoted as layer 0, the hidden layer is denoted as layer 1, and the output layer is denoted as layer 2. Table 2 gives the relevant notation of the FNN used in this work.  If a (0) = (x, y) T , the FNN can propagate information by iterating Equation (17): First, the net activation of layer l is calculated according to the activation of layer l-1 and the weight matrix and bias vector of layer l. Then, the activation of layer l is calculated according to the nonlinear activation function. Equation (17) can be written in a compact form: According to several tests, setting M 1 = 20 can provide enough design freedom without significantly increasing the computational burden of the optimization. The parameters of the FNN structure used in this study are shown in Table 3. So far, the specific structure of the FNN has been determined, as shown in Figure 9. Table 3. Parameters of the FNN structure.

Items
Value The total number of tunable hyperparameters of the neural network can be found in Equation (19): Substituting the data in Table 3 yields N hp = 81. After comparing with Equation (16), it can be seen that the vector p is an 81-dimensional vector, and its elements are composed of all the elements in the matrices W (1) , W (2) and the vectors b (1) , b (2) .

Pre-Training Based on Supervised Learning
To obtain faster optimization speed, an appropriate initial guess of vector p 0 needs to be selected according to specific prior knowledge to transform the problem into optimizing correction ∆p, as shown in Equation (20): where the range of ∆p is chosen to be [−1, 1] to avoid holes in the grain. This technique is denoted as pre-training in this work. In the rest of this section, p 0 will be obtained via supervised learning from a typical tube grain shape. The phase field of a typical tube grain shape can be written in analytical form. Then, the position coordinates and phase field of the grain are taken as training samples, which can be expressed as Equation (21): where M is the total number of the samples and ϕ (n) 0 is the phase field of the typical tube grain shape.
The loss function in the form of mean square error (MSE) is defined as follows: whereφ p, x (n) , y (n) represents the estimated value of the phase field. This work obtains the corresponding parameter vectors using the Bayesian regularization training method [38]. This method effectively avoids overfitting and has an excellent fitting effect for complex problems. The variation curve of the loss function with the training epoch is shown in Figure 10. The output error of the neural network is shown in Figure 11. The mapped phase field of the tube grain is shown in Figure 12.   As shown in Figure 10, the MSE decreases monotonically with the training epoch. When the training epoch reached 1000, the MSE reached the minimum value 5.8844 × 10 −8 . In addition, the variation of MSE between the training set and the test set is relatively small, indicating that no overfitting has occurred in the training process. It can be seen from Figure 11 that the output error of the network is at a small level. In summary, the pre-training can obtain relatively accurate initial parameter values p 0 .

Genetic Algorithm Framework
The genetic algorithm (GA) uses neither sensitivity derivatives nor reasonable starting points. As a direct search-based method, it can solve the global optimization problems with a non-analytical objective function, which is suitable for optimizing hyperparameters of neural networks. Therefore, in this work, GA is used to optimize the neural network to achieve a phase field matching the performance best. The neural network combined with GA is also known as the evolutionary neural network. Figure 13 illustrates the framework of the GA, where the optimized variable is the ∆p in Equation (20), containing a total number of 81 hyperparameters. The next section will focus on the computation of the objective function of GA.

Objective Function
Determining a suitable objective function is one of the most challenging tasks in the reverse design. Using different objective functions will lead to different results. The objective function shall compare the similarity betweenl b,D ∼ŵ D curve and l b,T C T ∼ŵ T curve, and output the performance-matching degree (see Figure 4).
To evaluate the similarity between curves, a method based on linear interpolation combined with linear regression is proposed. For the convenience of writing, the function obtained by linear interpolation ofl b,D ∼ŵ D is denoted as f 1 , and the function obtained by linear interpolation of l b,T C T ∼ŵ T is denoted as f 2 . Take n nodes at equal intervals within the value range of the independent variable to form a vector x: The vectors y 1 and y 2 is calculated using f 1 and f 2 , and the formula is as follows: If f 1 and f 2 are approximately proportional, then there is: Equation (25) contains one unknown variable C D and n equations. Since the number of equations is much larger than the number of unknowns, it is overdetermined and generally has no solution. However, the linear regression method can be used to find the solution in the presence of the smallest residual error, that is, the least squared solution. The unknown variable C D is approached using Equation (26): According to the discussion in Section 2.1.2, C D can be used as the design result of LR/A t . In order to quantitatively evaluate the similarity between f 1 and f 2 , the coefficient of determination (COD) can be calculated as Equation (27): where y 2 is an n-dimensional vector consisting of the mean value of y 2 . The COD is a number never greater than 1. The closer it is to 1, the more similar are the f 1 and f 2 . When COD is equal to 1, Equation (25) will be strictly satisfied. Since the GA algorithm always seeks the smallest objective function, we have Equation (28): where FV is the abbreviation of Fitness Value. In summary, the goal of the GA algorithm is to minimize Equation (28).

Consideration of Constraints
Constraints can be added to GA in the form of penalty functions. The specific implementation is to add a very large number to the objective function if the current solution violates the constraints, shown as Equation (29): where P is the penalty functions. There are two types of constraints that need to be emphasized.

1.
Rationality of the grain structure.
As shown in Figure 14, if the grain is completely separated from the chamber wall or there are holes in the grain, it is unacceptable. In practice, holes can be basically eliminated by choosing the range of ∆p to be [−1, 1]. If for every x ∈ wall there is ϕ < 0, then it is known that the grain is completely separated from the chamber wall, and the objective function should be penalized for these cases.
C D in Equation (26) is the design result of LR/A t . It can qualitatively describe the loading fraction. The loading fraction decreases with the increase of C D . In order to obtain a more reasonable loading fraction in engineering, C D range should be selected according to the actual needs. The objective function should be penalized for cases that violate the range.
So far, the method of grain reconstruction is established. The reverse design can be achieved via two separate steps, including reverse internal ballistic calculation and grain reconstruction.

Results and Discussion
This chapter will give three examples of shape optimization. Sections 3.1 and 3.2 are the reconstruction of given grains, excluding the reverse internal ballistic calculation.
Starting from the l b,T C T ∼ŵ T curves of given grains, they are dedicated to finding results most matching the given grains. Based on Sections 3.1 and 3.2, a complete reverse design, including two steps denoted as the reverse internal ballistic calculation and the grain reconstruction, is conducted in Section 3.3. The purpose is to find the grains most matching the given p c ∼ t curve without any additional prior knowledge.

Design Target
Star shape is one of the most fundamental grain shapes. In this section, starting from the l b,T C T ∼ŵ T curve of a provided star grain, the reconstruction of the grain is conducted with the evolutionary neural network. The geometry of the provided star grain is shown in Figure 15. The number of the geometry units is 12, and the throat area is 5 × 10 −2 m. The PEF method is used to obtain thel b,T ∼ŵ T curve. Considering Equation (7), there is Equation (30): Then the l b,T C T ∼ŵ T curve is calculated and plotted in Figure 16. Figure 16. Relation between the dimensionless burn perimeter and the dimensionless web thickness of the given star grain.

Design Result
Firstly, the computational domain is restricted to a fan shape which is 1/N of a whole circle with an inner radius of 0.2R and an outer radius of R. The initial combustion surface is set to be a circular arc at the inner radius, and the rest of the boundaries are flame retardant. Then, the computational domain is divided into 1000 s-order quadrilateral elements, as shown in Figure 17. Secondly, a phase field is defined using the typical tube grain. According to this phase field, the pre-training of the neural network is completed, as shown in Figure 12.
Thirdly, taking Equation (29) as the objective function, setting the population to 200 and generations to 50, the neural network is evolved via the GA algorithm. If the reduction rate of the objective function is negligible, it can be considered that the optimal shape of the grain is obtained.
The evolutionary process in Figure 18 shows that the shape of the grain changes as the GA algorithm proceeds. The figure presents the best individuals in populations during generations. As the result of pre-training, G0 is a typical tube grain. After one generation of reproduction, G1 tends to be a star grain. However, the shape of G5 appears to be highly irregular. After G10, the irregular shape vanishes, and a continuous and smooth star grain is gradually formed. Using G50 as the final design result, linear regression is used to evaluate the gap between the design result and the goal, as shown in Figure 19. The linear regression is indicated in Equation (31): Comparing C D with C T , given in Equation (30), the relative error is calculated in Equation (32): Since R = 0.5 m, the designed grain must satisfy Equation (33): Equation (33) indicates that the length of the grain strictly corresponds to the throat area. Particularly, when A t,D = 5 × 10 −2 , the designed grain replicates the given one precisely. Furthermore, the dimensionless burn perimeter is plotted in Figure 20 as a function of dimensionless web thickness. Although the tail-off stage has few computational points, the linear interpolation in Equation (24) can ensure the accuracy of the objective function. After considering the combustion characteristics of the propellant in Table 4, the combustion chamber pressure is plotted in Figure 21 as a function of time. If the nozzle area ratio and ambient pressure are known, the thrust-time curve can also be obtained.  It can be seen from the figures that the design result has matched the target curve well. Consequently, the shape optimization method can precisely reconstruct a given star grain.

Design Target
Theoretically, the GA algorithm can easily find the global optimal solution of the hyperparameters in the neural network. However, designers often need to consider several constraints such as loading fraction, grain length, and throat to port area. To obtain more candidates, the designers always wish to find a series of suboptimal feasible solutions instead of a single optimal solution.
In this section, taking a given typical slotted-tube grain as the design target, several suboptimal feasible solutions satisfying the burning surface area requirements are found by setting suitable criteria for C D . On one hand, a larger C D criteria will return grains with a lower loading fraction. On the other hand, a smaller C D criteria will return grains with a higher loading fraction.
The geometry of the provided slotted-tube grain is shown in Figure 22. The number of the geometry units is 12, and the throat area is 5 × 10 −2 m. Then, the l b,T C T ∼ŵ T curve is calculated from this grain.

Design Result
The meshing technique and the pre-training of the neural network are not repeated here. By appropriately constraining the range of C D , four suboptimal solutions with different loading fractions can be obtained, as shown in Figure 23. The dimensionless burn perimeter of these schemes is shown in Figure 24.  The design result reveals that although the grain shapes of the schemes are quite different, they all share identical burning surface areas or identical performance. The significant difference between them is the loading fraction. The loading fraction of Scheme A is 66.26%, which means the grain length must be long enough to satisfy the performance of the given slotted-tube grain. Contrarily, the loading fraction of Scheme D is 86.5%, which means the grain length can be much shorter. Since the loading fraction of Scheme C and Scheme D are higher than the given one, the method successfully obtains several solutions even better than the original one. Consequently, these schemes could provide designers with more options.

Design Target
Dual-thrust motors are widely used in tactical missiles to increase the flight range significantly. The booster stage provides a high but short-term thrust, followed by the cruising stage with low and long-term thrust. Some 2D grain shapes, such as wagon wheel, dendrite, and dog bone [15], might produce dual-thrust internal ballistics. However, there are still questions about whether these grain shapes are the best options. Therefore, it is necessary to design a more suitable grain shape using shape optimization.
The dual-thrust performance curve is given with the aid of a piecewise function shown as Equation (34): The parameters of the dual-thrust grain are given in Table 5. According to Equations (6) and (34), and Table 5, the l b,T C T ∼ŵ T curve can be derived from the p c~t curve, as shown in Figure 25.

Design Result
The evolution process of the dual-thrust grain is illustrated in Figure 26. The figure presents the best individuals in populations during generations. G0 is a typical tube grain whose performance fails to match the given one. After the random population initialization, one individual denoted as G1 among the whole population, appears to match the performance most. Interestingly, the grain shape of G1 seems to be a combination of two different star grains. After 50 generations, the grain shape gradually stabilizes and becomes a star-wagon-wheel combination. After several generations, the results are almost identical, and it is reasonable to believe that a globally optimal solution is obtained. Coincidentally, the design result is similar to a new type of combined dendrite grain proposed by Jiang [1,39] through a traditional approach. The coefficient of determination (COD) is 0.95162, proving the proportional relationship between the two curves. As R = 0.5 m, the designed grain must satisfy Equation (36): If the outer radius of the grain is fixed to 0.5 m, the length of the grain will strictly correspond to the throat area.
The dimensionless burn perimeter is plotted in Figure 28 as a function of the dimensionless web thickness. Similarly, the combustion chamber pressure is plotted in Figure 29 as a function of time. It can be found that the designed pressure of the booster stage and the cruising stage is stable, which are 10 MPa and 5 MPa, respectively, so the design results have well-matched the expected dual-thrust performance. However, due to the slow burning rate of the residual propellant, the tail-off stage lasts up to 5 s. In order to eliminate the tail-off impulse, it is necessary to replace the residual propellant with a non-flammable filler in further consideration.  It should be noticed that the current results are obtained under the condition of N = 12. However, N is not a given integer in practice, so it is worth investigating the design result with different N, as shown in Figure 30. The relation between the chamber pressure and the burning time with different N is shown in Figure 31.  As can be seen from the above figure, the reverse design result is related to N. When N is small (8 and 10), the results are similar to the dog bone grain. When N is 12 and 16, the results turn into the combined dendrite grain as mentioned before. When N is large (greater than 20), the results turn into the wagon wheel grain. As the outer radius of the grain limits the length of the dendritic wheel arm, and if N is sufficiently small, the combined dendrite grain will no longer maintain the performance curve, and the design result will turn into a dog bone. It can also be found that with the increase of N, the residual propellant becomes smaller, and results with a higher performance-matching degree can be obtained.
In conclusion, the shape optimization method based on the evolutionary neural network can realize the reverse design through an example of a dual-thrust propellant grain. The brand-new grain shape will better match the performance curve and provide designers with numerous innovative ideas.

Conclusions
The reverse design of solid propellant grain for a performance-matching goal is investigated, and the conclusions are as follows: 1.
The major contribution of this work is to propose a shape optimization method based on the evolutionary neural network to achieve grain reverse design. It does not require any pre-selection of the grain shape and the designers shall be free from defining different kinds of geometric parameters for specific grain configurations.

2.
Burn-back analysis method (ϕ-PEF) of irregular grain shape on a fixed unstructured mesh is proposed. The method can be directly solved by the steady-state heat conduction modules in commercial finite element platforms.

3.
The neural network is designed to determine and manipulate the phase field of the grain. Since the degree of freedom is controllable (much higher), it will form much more grain shapes and gives designers more options.

4.
The genetic algorithm framework is established to optimize the hyperparameters in the neural network. The objective function that evaluates the performance-matching degree is the coefficient of determination (COD) derived from the linear regression. 5.
The computation results show that the method can precisely match the performance of the given star grain and slotted-tube grain. By adding constraints, it can achieve other performance requirements, such as higher loading fraction. 6.
The method can automatically evolve a new dendritic-shaped grain that matches the given dual-thrust pressure-time curve. Consequently, the method has the potential to apply to the conceptual design of innovative grain shapes through the reverse design procedure, which can hardly be predefined and achieved by the traditional approaches.
Future work will focus on extending the method to three-dimensional (3D) grains. Due to the heavy computational burden of 3D problems, the difficulties will lie in improving the time efficiency of the burn-back analysis, introducing a well-designed deep neural network, and adopting a parallel genetic algorithm. In addition, manufacturability, mechanical performance, and new 'pre-set' parameters of grains should also be carefully considered.