# A New Method for Initial Parameters Optimization of Guided Projectiles

Authors: Feng Bi-Ming; Nie Wan-Sheng
Journal: International Journal of Aerospace Engineering
 (2013)

---

## Abstract

A new algorithm was developed for the initial parameters optimization of guided projectiles with multiple constraints. Due to the relationship between the analytic guidance logic and state variables of guided projectiles, the Radau pseudospectral method was used to discretize the differential equations including control variables and state variables with multiple constraints into series algebraic equations that were expressed only by state variables. The initial parameter optimization problem was transformed to a nonlinear programming problem, and the sequential quadratic programming algorithm was used to obtain the optimal combinations of initial height and range to target for the final velocity of guided projectiles maximum with constraints. Comparing with the appropriate initial conditions solved by Monte Carlo method and the flight characteristics solved by integrating the original differential equations in the optimal initial parameters computed by the new algorithm, the feasibility of new algorithm was validated.


---

## Body

## 1. Introduction

The initial parameter optimization problem of guided projectiles with multiple constraints is considered. In addition to the constraints on its terminal position coordinates, the projectile must also impact the target from a specified direction with the required circular error probable (CEP) and make sure the impact velocity maximum.Presently, several studies are conducted for developing the guidance logic with multiple constraints [1, 2] and searching the launch acceptable region (LAR) of a guided bomb, while the initial conditions are known [3, 4]. However, a few studies pay attention to the initial parameters optimization for the guided projectiles with multiple constraints. The appropriate parameters are calculated based on different combinations of initial conditions (height above sea level, velocity, and path angles) by searching method to find out the largest LAR, but they are not optimal conditions [5]; the indirect method and direct method both are called calculative method which transforms the parameters optimization to the optimal control problem. Indirect method based on the Pontryagin maximum principle works out the analytical optimal control logic to obtain the optimal index for the guided projectiles in known initial conditions, for instance, maximum impact velocity or least energy [6, 7]; direct method based on pseudospectral method discretizes the equations of motion to obtain the numerical optimal control law for making parameters optimal with constraints [8, 9]. Now, when the analytical guidance logic exists, how to find out the optimal initial parameters of the guided projectiles with multiple constraints to make the final velocity maximum becomes a new problem.Due to the relationship between the analytical guidance logic and state variables, the 2-dimensional equations of motion with both control variables and state variables are described by the differential equations only with the state variables. Then the new differential equations are discretized to be series algebraic equations expressed by state variables, while the pseudospectral method is used [10]. So, the initial parameters optimization is transformed to a nonlinear programming problem which can be solved using a sequential quadratic programming (SQP) algorithm to improve computation efficiency and reduce the complexity of initial condition searching.


## 2. Theoretical Development

### 2.1. Guidance Logic with Impact Direction Constraint in Longitudinal Plane

For development purposes, a simplified flat earth coordinate system is defined to be shown in Figure1. The target is at the origin of the coordinate system. We will focus on a fixed target in the following discussion. The x-axis is pointed to the projectile, and the y-axis is vertical to the horizontal plane. The line-of-sight (LOS) from the origin to the projectile is defined by the elevation angle λD, where 0≤λD≤π/2 is measured from positive x-axis in a clockwise direction. The flight path angle θ is the angle between the relative velocity vector and the local horizontal plane, and -π/2≤θ≤0 is measured from local horizon in a clockwise direction as defined in Figure 1. γD (-π/2≤γD≤0) is the azimuth angle of relative velocity vector in vertical plane as defined in Figure 1.Figure 1
Coordinate system and geometry.By the definition of the coordinate system shown in Figure1, the standard 2D equations of motion of the projectile in the longitudinal plane can be described as follows:
(1)V˙=-ρV2SCx(α)2M+gxh,θ˙=1V(ρV2SCy(α)2M+gyh),x˙=-Vcosθ,y˙=Vsinθ.As described in (1), the position coordinates are x and y. The earth-relative velocity is V, the mass of projectile is M, the atmospheric density model is ρ=ρ0·e(-y/H), the angle of attack is α, and gxh and gyh can be defined as follows:
(2)gxh=-ghsinθ,gyh=-ghcosθ.The additional parameters are defined as follows:gh=9.81m/s2: gravity acceleration,ρ0=1.226kg/m3: standard sea-level density,H=7254.24m: atmospheric scale height,S=0.102m2: aerodynamic reference area.In this vertical plane the guidance logic of the projectile with impact direction constraint may be represented by [6]
(3)γ˙D=KGDλ˙D+KLD(λD+γDF)Tg.HereγDF is defined as the desired impact direction.


### 2.2. Differential Equations Expressed by State Variables

Before the parameter optimization, the equations of motion with control variable are transformed to the differential equations only including state variables. By the discretization of new differential equations, the equation of motion can be described by some algebraic equations. So, the initial parameter optimization can be transformed to the nonlinear programming problem.As shown in Figure1, the λD is also easily obtained:
(4)λD=arctan(yx).So, the derivative ofλD can be described as
(5)λ˙D=-VcosθsinλD+VsinθcosλDx2+y2,Tg=x2+y2VcosθcosλD+VsinθsinλD.Using (5), the guidance logic equation (3) can be represented as follows:
(6)γ˙D=KGD-VcosθsinλD+VsinθcosλDx2+y2+KLD(λD+γDF)VcosθcosλD+VsinθsinλDx2+y2.While the equationθ˙=γ˙D satisfies in the vertical plane, the lift coefficient Cy(α) yields from (1):
(7)Cy(α)=2(γ˙DV-gyh)·Mρ0·e(-y/H)V2S.The aerodynamic coefficientsCyα, Cx0, and Cxα come from [11], and the lift coefficient can be defined as Cy(α)=Cyαα, so the α is obtained:
(8)α=Cy(α)Cyα=1Cyα·2(γ˙DV-gyh)·Mρ0·e(-y/H)V2S.Because the drag coefficient is defined asCx(α)=Cx0+Cxαα2, the drag force D can be readily achieved:
(9)D=ρ0·e(-y/H)V2S2{Cx0+Cxα[Cy(α)Cyα]2}.By (9) and (6), (1) with control variable α can be transformed to differential equations expressed by state variables [x,y,V,θ]; it is described as follows:
(10)dVdt=-ρ0·e(-y/H)V2S2M{Cx0+Cxα[2(γ˙DV-gyh)·MCyαρ0·e(-y/H)V2S]2}+gxh,dθdt=KGD-VcosθsinλD+VsinθcosλDx2+y2+KLD(λD+γDF)VcosθcosλD+VsinθsinλDx2+y2,dxdt=-Vcosθ,dydt=Vsinθ.If the initial statex0 is known, the only final state xf is readily obtained by integrating (10). So, the initial parameter optimization problem can be described as follows.Cost function:(11)minφ(x0,xf).Subject to,the dynamic constraints:(12)x˙=f(x),the inequality constraints:(13)|α|=|1Cyα·2(γ˙DV-gyh)·Mρ0·e(-y/H)V2S|≤αLIMIT,|ny|=|γ˙DV-gyh|gh≤nyLIMIT,|xf|≤CEP,x0>0,y0>0,Vf>0,the boundary conditions:(14)θ0=θinitial,V0=Vinitial,yf=yterminal,θf=γDF.


### 2.3. Discretization of Differential Equations

For instance, the continuous differential equations are discretized into the algebraic equations by the Radau pseudospectral method. The parameter optimization problem is described as follows.Cost functional:(15)minJ≈Φ(X1(1),t0,XNK+1(K),tK).Equation constraint:equation of motion:(16)∑j=1Nk+1Dij(k)Xj(k)-tk-tk-12f(Xi(k),τi(k);tk-1,tk)=0(i=1,…,Nk),equation of boundary condition:(17)ϕ(X1(1),t0,XNK+1(K),tK)=0.Inequation constraint:angle of attack and overloading constraint:(18)C(k)(Xi(k),τi;tk-1,tk)≤0(i=1,…,Nk),inequation of boundary condition:(19)ϕ(X1(1),t0,XNK+1(K),tK)≤0.Then the initial parameter optimization problem of the equations of motion that satisfied the guidance logic is transformed to the nonlinear programming problem of algebraic equations expressed by state variables. The optimal initial parameter can be computed by SQP.


## 2.1. Guidance Logic with Impact Direction Constraint in Longitudinal Plane

For development purposes, a simplified flat earth coordinate system is defined to be shown in Figure1. The target is at the origin of the coordinate system. We will focus on a fixed target in the following discussion. The x-axis is pointed to the projectile, and the y-axis is vertical to the horizontal plane. The line-of-sight (LOS) from the origin to the projectile is defined by the elevation angle λD, where 0≤λD≤π/2 is measured from positive x-axis in a clockwise direction. The flight path angle θ is the angle between the relative velocity vector and the local horizontal plane, and -π/2≤θ≤0 is measured from local horizon in a clockwise direction as defined in Figure 1. γD (-π/2≤γD≤0) is the azimuth angle of relative velocity vector in vertical plane as defined in Figure 1.Figure 1
Coordinate system and geometry.By the definition of the coordinate system shown in Figure1, the standard 2D equations of motion of the projectile in the longitudinal plane can be described as follows:
(1)V˙=-ρV2SCx(α)2M+gxh,θ˙=1V(ρV2SCy(α)2M+gyh),x˙=-Vcosθ,y˙=Vsinθ.As described in (1), the position coordinates are x and y. The earth-relative velocity is V, the mass of projectile is M, the atmospheric density model is ρ=ρ0·e(-y/H), the angle of attack is α, and gxh and gyh can be defined as follows:
(2)gxh=-ghsinθ,gyh=-ghcosθ.The additional parameters are defined as follows:gh=9.81m/s2: gravity acceleration,ρ0=1.226kg/m3: standard sea-level density,H=7254.24m: atmospheric scale height,S=0.102m2: aerodynamic reference area.In this vertical plane the guidance logic of the projectile with impact direction constraint may be represented by [6]
(3)γ˙D=KGDλ˙D+KLD(λD+γDF)Tg.HereγDF is defined as the desired impact direction.


## 2.2. Differential Equations Expressed by State Variables

Before the parameter optimization, the equations of motion with control variable are transformed to the differential equations only including state variables. By the discretization of new differential equations, the equation of motion can be described by some algebraic equations. So, the initial parameter optimization can be transformed to the nonlinear programming problem.As shown in Figure1, the λD is also easily obtained:
(4)λD=arctan(yx).So, the derivative ofλD can be described as
(5)λ˙D=-VcosθsinλD+VsinθcosλDx2+y2,Tg=x2+y2VcosθcosλD+VsinθsinλD.Using (5), the guidance logic equation (3) can be represented as follows:
(6)γ˙D=KGD-VcosθsinλD+VsinθcosλDx2+y2+KLD(λD+γDF)VcosθcosλD+VsinθsinλDx2+y2.While the equationθ˙=γ˙D satisfies in the vertical plane, the lift coefficient Cy(α) yields from (1):
(7)Cy(α)=2(γ˙DV-gyh)·Mρ0·e(-y/H)V2S.The aerodynamic coefficientsCyα, Cx0, and Cxα come from [11], and the lift coefficient can be defined as Cy(α)=Cyαα, so the α is obtained:
(8)α=Cy(α)Cyα=1Cyα·2(γ˙DV-gyh)·Mρ0·e(-y/H)V2S.Because the drag coefficient is defined asCx(α)=Cx0+Cxαα2, the drag force D can be readily achieved:
(9)D=ρ0·e(-y/H)V2S2{Cx0+Cxα[Cy(α)Cyα]2}.By (9) and (6), (1) with control variable α can be transformed to differential equations expressed by state variables [x,y,V,θ]; it is described as follows:
(10)dVdt=-ρ0·e(-y/H)V2S2M{Cx0+Cxα[2(γ˙DV-gyh)·MCyαρ0·e(-y/H)V2S]2}+gxh,dθdt=KGD-VcosθsinλD+VsinθcosλDx2+y2+KLD(λD+γDF)VcosθcosλD+VsinθsinλDx2+y2,dxdt=-Vcosθ,dydt=Vsinθ.If the initial statex0 is known, the only final state xf is readily obtained by integrating (10). So, the initial parameter optimization problem can be described as follows.Cost function:(11)minφ(x0,xf).Subject to,the dynamic constraints:(12)x˙=f(x),the inequality constraints:(13)|α|=|1Cyα·2(γ˙DV-gyh)·Mρ0·e(-y/H)V2S|≤αLIMIT,|ny|=|γ˙DV-gyh|gh≤nyLIMIT,|xf|≤CEP,x0>0,y0>0,Vf>0,the boundary conditions:(14)θ0=θinitial,V0=Vinitial,yf=yterminal,θf=γDF.


## 2.3. Discretization of Differential Equations

For instance, the continuous differential equations are discretized into the algebraic equations by the Radau pseudospectral method. The parameter optimization problem is described as follows.Cost functional:(15)minJ≈Φ(X1(1),t0,XNK+1(K),tK).Equation constraint:equation of motion:(16)∑j=1Nk+1Dij(k)Xj(k)-tk-tk-12f(Xi(k),τi(k);tk-1,tk)=0(i=1,…,Nk),equation of boundary condition:(17)ϕ(X1(1),t0,XNK+1(K),tK)=0.Inequation constraint:angle of attack and overloading constraint:(18)C(k)(Xi(k),τi;tk-1,tk)≤0(i=1,…,Nk),inequation of boundary condition:(19)ϕ(X1(1),t0,XNK+1(K),tK)≤0.Then the initial parameter optimization problem of the equations of motion that satisfied the guidance logic is transformed to the nonlinear programming problem of algebraic equations expressed by state variables. The optimal initial parameter can be computed by SQP.


## 3. Results and Discussions

Assuming that the initial path angle and velocity exist, the optimal initial height and range to the target are solved to get the maximum impact velocity. The problem can be described as follows.Cost function:(20)min-Vf.Subject to,the dynamic constraints:(21)x˙=f(x),the inequality constraints:(22)[-30∘-6000-10m]≤[αnyy0x0xf]≤[30∘6020km50km10m],the boundary conditions:(23)[θ0V0yfθf]=[-15∘Vinitial0-90∘].Suppose that the mass of guided projectile is 450 kg; the optimal parameters for maximum final velocity in different initial velocities are shown in Table1. With the increase of initial velocity, the optimal initial height becomes higher, the optimal range to target becomes farther, and the maximum impact velocity becomes faster. Comparing with the results as shown in Figure 2, the optimal heights and ranges to target for maximum impact velocities in different initial velocities are almost similar to the results computed by Monte Carlo method.Table 1
The optimal initial position in different velocities for maximum impact velocity.
θ
0 (°)
V
0 (m/s)
y
0 (km)
x
0 (km)
V
f (m/s)
−15
3200
12.66
10.58
1945
−15
3400
13.80
11.48
2045
−15
3600
15.28
13.19
2134
−15
3800
16.76
14.69
2205Distribution of the impact velocity and range to target computed by Monte Carlo method.
(a)
Distribution for initial velocity 3200 m/s
(b)
Distribution for initial velocity 3400 m/s
(c)
Distribution for initial velocity 3600 m/s
(d)
Distribution for initial velocity 3800 m/sFor instance, when the initial velocity and path angle are, respectively, 3400 m/s and −15°, the optimal height and range to target are, respectively, 13.8 km and 11.48 km (shown in Table1), and the maximum terminal velocity is 2045 m/s. Compared to the results shown in Figure 2(b), the results find that while the initial height becomes higher or lower than 13.8 km, the maximum final velocity becomes slower. As shown in Figure 2(b), when the initial height is 13.8 km, the range to target for the maximum final velocity is close to 11.5 km, and the impact velocity is close to 2050 m/s.Using the initial parameters as shown in Table1, the equation of motion (1) with guidance logic (3) is integrated to get the flight characteristics of projectile in different initial velocities. Analyzing the results shown in Figure 3(a) to Figure 3(d), the result finds that the trajectories of projectile for maximum impact velocity are close to the parabolas, and the final impact direction, AOA, and overloading satisfy the constraints. So, the feasibility of new initial parameter optimization method used in this work is validated by comparison to the results shown in Figures 2 and 3. The new method can obtain the optimal initial parameters of guided projectile for different cost functions with the multiple constraints.Comparison of flight characteristics of projectile in different initial condition.
(a)
Comparison of height and range to target
(b)
Comparison of velocity and path angle
(c)
Comparison of AOA
(d)
Comparison of overloadingAs shown in Figure3(d), the maximum overloading of projectile is saturated for long time during flight in different initial conditions; it goes against the adaptability of projectile to the error and disturbance. In fact, the maximum overloading should be smaller than the limit. With the constraint ny<60, the optimal initial parameters for maximum impact velocity are solved by the new method as shown in Table 2.Table 2
Optimal initial position with more severe constraint.
θ
0 (°)
V
0 (m/s)
y
0 (km)
x
0 (km)
V
f (m/s)
−15
3400
15.781
14.098
2021
−15
3600
16.562
14.149
2114
−15
3800
18.232
16.094
2175Contrasting the results shown in Table1, in the same initial velocities, the initial heights are higher, and the ranges to target are farther, but the maximum final velocities with constraint ny<60 are slower than the results with constraint ny≤60.The comparison of flight characteristics in different initial conditions set out in Table2 is shown in Figure 4. The maximum overloading in different initial conditions as shown in Figure 4(d) is smaller than 60; it is propitious to the adaptability of projectile to the error and disturbance.Comparison of flight characteristics with severe overloading constraint.
(a)
Comparison of height and range to target
(b)
Comparison of velocity and path angle
(c)
Comparison of AOA
(d)
Comparison of overloading


## 4. Conclusions

A new algorithm is developed for initial parameters optimization of guided projectiles by discretizing equations of motion with guidance logic based on pseudospectral method. The method can solve the initial parameters optimization problem with multiple constraints using SQP. The initial parameters solved by the new method for maximum impact velocity are similar to the results solved by Monte Carlo method, and the flight characteristics solved by integrating the equations of motion in the optimal initial conditions satisfy the all constraints.



---

*Source: 486020-2013-12-11.xml*