Suboptimal Impact-Angle-Constrained Guidance Law Using Linear Pseudospectral Model Predictive Spread Control

This paper aims at proposing a nonlinear suboptimal impact-angle-constrained guidance law for an interceptor engaging against an incoming high-speed ballistic missile in three dimensions. The method is based on an improved Linear Pseudospectral Model Predictive Control (LPMPC), in which the control variable is parameterized as a sum of general standard basis functions with few undetermined coefficients. Although the final solution is suboptimal, a series of analytical improved formulas can be also successively derived to eliminate the final predictive errors, which remains the character that few points will achieve high accuracy. Furthermore, the number of improved variables will significantly decrease, which will further improve the computational efficiency and provide smoother control history. Several numerical simulations are carried out to evaluate the performance of the proposed method. Results show that it performs well in achieving desired impact angles and the near-zero miss distance with high accuracy. In comparison with other typical methods, the proposed method has superior performance in computational efficiency. Therefore, it has the potential to be applied in the framework of guidance.


I. INTRODUCTION
Impact-angle-constrained guidance aims at achieving engagement with a desired impact angle, which will significantly increase the destructiveness of direct hit or warhead carried by missile. Therefore, it has drawn a great attention in modern warfare in recent decades. Different battle scenarios have different requirements on the impact angle. For the anti-tank or anti-ship missile it is desirable to achieve a top attack because the top of a tank or ship is usually vulnerable. For the anti-ballistic missile, the relative velocity between the ballistic target and the interceptor is so great that a head-on engagement is preferable to maximize the impact velocity and the lethality of the interceptor. It should be noted that terminal guidance design for intercepting a ballistic missile is one of the most challenging problems because it engages with high relative speed, high rate of line of sight and thin atmosphere. It has high requirements on miss distance and impact angle. Additionally, the lateral acceleration commands should be as The associate editor coordinating the review of this manuscript and approving it for publication was Rosario Pecora . small and smooth as possible, which is easy for the inner-loop autopilot to execute.
The classical proportional navigation (PN) guidance [1] is one of the most famous and widely used guidance all over the world since it is very easy to be implemented. It should be noted that original PN guidance does not consider the impact angle constraints. And following researchers have made great efforts to overcome this drawback. Kim, Lee and Ham firstly proposed a biased PNG law, in which a time-varying bias term is used to shape the trajectory so as to achieve the desired terminal impact [2]. Ratnoo and Ghose introduced a special PN guidance with switched gain, which is capable of satisfying the constraint on impact angle and has great performance [3]. Jeong et al derived the analytical bias term for a specified impact angle on the assumption of the linearized engagement, on which a new biased PN guidance is proposed [4]. Akhil and Ghose further developed a time-variant biased PNG law in which nonlinear engagement model is fully involved [5].
Note that, because above guidance laws are under the assumptions of the PN philosophy, they cannot provide optimal guidance commands and suffer from many inherent drawbacks such as sharp lateral accelerations at the end of engagement, no trajectory shaping capability, and etc., [1], [6], [7]. Therefore, many researchers pay attention to developing the impact-angle-constrained guidance under the framework of optimal control theory [8]- [15] which can fully exert the maneuverability so as to increase the operational effectiveness for missile weapon system. Bryson is one of the most notable pioneers pioneering in applying the optimal control theory into the missile guidance [8]. Bryson regarded the nonrestrictive and truly optimal impact-angle-constrained guidance problem as a typical optimal control problem, in which nonlinear engagement dynamics were fully formulated, zero miss distance and terminal angle constraints were enforced strictly, and the lateral acceleration demand throughout the engagement were regarded as performance. Undoubtedly, it is very hard to solve this complicated optimal control problem for previous computational technology. Therefore, simplified and linearized modeling assumptions were employed to derive analytical optimal solution. Kim and Grider firstly proposed an optimal terminal angle constraint guidance law on the foundation of the linearized modeling [9]. Ratnoo and Ghose have proposed a state-dependent Riccati equation technique to satisfy the impact angle constraint [10]. Song et al. developed an optimal impact angle guidance law considering maneuvering target [11].
With the development of numerical computational technology, many numerical algorithms capable of solving the nonlinear optimal control problem have been emerged. Padhi has proposed a nonlinear optimal control algorithm named model predictive static programing (MPSP) to solve the nonlinear optimal control problem with hard terminal constraints and quadratic performance [16]. This method combines the philosophy of approximate dynamic programming [17] and model predictive control [18]. The control solution can be analytically calculated in a recursive and iterative manner. Therefore, it has high computational efficiency in comparison with the conventional numerical methods and dynamic programming approach. Also, this method has been successfully applied to formulate several optimal guidance laws with terminal angle constraint [6], [19]- [21]. However, to ensure the Euler integration within a desired precision limitation, a large number of discrete points should be selected for MPSP. In order to overcome this drawback, Yang has proposed an optimal control method named linear pseudospectral model predictive control (LPMPC) to solve the same problem as the MPSP. This method combines the idea of linear quadrature optimal control, Gauss Pseudospectral method [22] and model predictive control. The nonlinear optimal control problem is transferred into a set of linear algebraic equations by linearization and Gauss Pseudospectral method. Then, the current control is obtained via iteratively solving a set of linear algebraic equations. Therefore, this method has attractive performance of high computational efficiency and accuracy with fewer collocation points. It is suitable for online applications. This method has been used to develop several challenging optimal guidance laws successfully [23]- [27].
Inspired by the idea of parameterized representation of control variable and the character of LPMPC that the control solution can be expressed as a smooth function at collocation points, a further extended LPMPC method is proposed for further improving the computational efficiency in this paper. This method considers a parameterized version of the control variable as a sum of general standard basis functions with undetermined coefficients. Those standard basis functions can be selected as Legendre polynomial, Power series polynomial, Chebyshev polynomial according to the specified engineering problem. It spreads the control requirement through the whole control history, which will provide smoother control history. Therefore, this extended method is called linear pseudospectral model predictive spread control (LPMPSC). Although the final solution is suboptimal, a series of analytical improved formulas can be also successively derived to eliminate the final predictive errors, which remains the character that few points will achieve high accuracy. Additionally, the most attractive point for this method is that the number of updated variables will significantly decrease, which will further improve the computational efficiency. Then, the LPMPSC is applied to solve the impact angle constrained guidance problem so as to intercept an incoming high-speed ballistic missile target. Several simulations have been carried out to evaluate the performance of the proposed method. Results show that it performs well in achieving desired impact angles and the near-zero miss distance with high accuracy. In comparison with other typical methods, the proposed method has superior performance in computational efficiency and providing more smooth control history. Therefore, it has the potential to be applied in the framework of guidance.
The rest of paper is structured as follows: Section 2 presents the mathematical formulation of LPMPSC; Section 3 develops the suboptimal impact-angle-constrained guidance law in detail; Section 4 provides the simulation results to demonstrate the performance of the method. Finally, conclusion is given in Section 5.

II. MATHEMATICAL FORMULATION OF LPMPSC
In this section, the mathematical formulation of LPMPSC is described in detail. The main improvement of LPMPSC here is a priori parameterizations of control, which leads to some further advantages over the previous methods such as decreasing the number of improved variables, improving the computational efficiency, and guaranteeing smoother control.

A. LINEARIZATION OF NONLINEAR DYNAMIACS AND PARAMETERIZATION OF CONTROL COMMAND
Without loss of generality, consider the general nonlinear dynamic system with hard terminal constraints as follows: where, x ∈ R n , u ∈ R m and t ∈ R denote the state vector, the control vector and the time variable, x f is the desired terminal constraints. δx f denotes the vector of terminal deviations from the desired state.
In the LPMPSC design, the control command is expressed as parameterized version using a general standard basis functions with undetermined coefficients. And then those coefficients are determined by LPMPSC optimization process. It should be noted that Legendre polynomial, Power series polynomial, Chebyshev polynomial and etc. can all be used as standard basis functions. Without loss of generality, quadratic formulation is considered as the standard basis function to express the control command in this paper.
Parameterized expression of each component u i of control vector u ∈ R m is given as where φ j (t) denotes the preselected basis functions, c ij is the weight of the function φ j (t) of component of control vector u i and needs to be determined by LPMPSC. N p gives the number of such basis functions in (3). This paper assumes that each component of the control vector can be represented as a sum of the same basis functions with different coefficients. An important advantage of this technique is that only a set of coefficients, c ij , instead of the control variables at every collection point should be updated. Thus, the number of updated variables will significantly decrease from N * m to N p * m, which will further improve the computational efficiency. Here, N is the number of Legendre Gauss (LG) points which will be introduced in subsection 2.2. Furthermore, N is much bigger than N p .
For example, if the optimal control is close to a linear function, the control can be considered to be a linear function of time in Fig. 1. The number of updated variables is only 2 * m.
If the optimal control is close to a quadratic function, the control can also be considered to be a quadratic function Next, expand (1) in Taylor series with high-order differentiation around the nominal trajectory. Then, neglecting the higher order terms and using the deviations as the independent variables, we can obtain a set of linear dynamics equations as where, the Jacobean matrices A andB are given by where x p , u p and c p are the nominal state vector, nominal control vector and nominal weight coefficient vector (which is the coefficients in (3)). Note that the real state vector is defined as x = x p − δx. Similarly, the real control vector and the real weight coefficient vector are given by u = u p − δu, c = c p − δc.

B. LINEAR PSEUDOSPECTRAL MODEL PREDICTIVE SPREAD CONTROL
In this section, the linear pseudospectral model predictive control is extended by considering the parameterized control command. In the linear pseudospectral method, the state and control variables are approximated using a basis of Lagrange interpolating polynomials, and the differential-algebraic equations are approximated via orthogonal collocation. Thus, the linear optimal control problem is transferred into the problem of solving a set of linear algebraic equations. And the computational time of solving this problem is reduced to only a fraction of a second. It is noted that Legendre Gauss (LG) [22], [28]- [30] collocation points are selected in this paper which are roots of a Legendre polynomial and has superior characteristics than other collocation points.
Firstly, because the LG points exit in the domain [−1, 1], we should convert the time domain [t 0 , t f ] of this linear optimal control problem (6) to [-1, 1] by the following transformation.
Then, the original linear optimal control problem is transferred into Then, we define that L N (τ ) is a Lagrange interpolating polynomial of degree N. And τ i are named LG points which are the roots of Legendre polynomial of degree N. It is noted that there is no analytical solution for these particular roots. Hence, we can obtain the LG points by using a numerical algorithm. According to Gauss Pseudospectral Method, the state, control and costate variables are approximated by a basis functions of Lagrange interpolating polynomials of degree N as According to the characteristics of the Lagrange interpolation polynomial, L N (τ ) has the following properties and At the LG points, the derivative of δx can be formulated using a differential approximation matrix as following (14).
where the differential approximation matrix D ∈ R N ×(N +1) can be obtained by the derivative of each Lagrange interpolating polynomial at the LG points, andD k ∈ R N is the first column of the differential approximation matrix D. The elements of D are determined as It is supposed that the state deviations can be given as a vector By substituting (14) into (8), the linear differential equations are not only transferred into a set of algebraic equations, but also represented by the deviations of the state vector at the LG points.
In order to derive a set of general analytical correction formulas for final errors correction, the approximating matrix is decomposed into two parts. And then, we rearrange the deviations of the state vector in (17), the algebraic equation is simplified as where, D 1 is one part of D associated with the initial variation of the state, D 2:n is the other part of D associated with the variation of the state at the LG points. The elements of D 1 and D 2:n are Note that superscript s gives the number of state variables, and the partial derivatives matrices A andB are represented as Next, rearranging (18) again, the deviations of the state vector at all LG points can be analytically represented as Then, because the LG points do not cover the boundary points, δx does not include the variation of final state. The variation of final state can be calculated via Gauss Quadrature rule. Then, integration of (10) is expressed as Hence, the final state deviation is represented as where, ω i are the integration weighting coefficients of Gaussian Quadrature rule. Rearranging (23), the variation of final state is represented as where, L is a column vector, W is a matrix which consists of the integration weighting coefficients of Gaussian Quadrature rule. The elements of those matrix and vector are derived as Adding (21) into (24) and rearranging it, it can be found that the terminal state variation can be expressed as a linear function of the deviation of initial state δx 0 and improvement of control parameter δc.
where, K x and K c are partial derivatives of the deviation of final state with regard to the deviation of initial state and the improvement of the control parameter, respectively. And they denote the effects of the initial state deviation and the improvement of control parameter on the final state deviation. The elements of K x and K c are given as where, the functions F x , F c are formulated in detail as It should be noted that the initial state deviation, δx 0 , is a zero vector. For some special cases, part of terminal state is free. Hence, matrix Y is used to eliminate the corresponding unknown terminal state deviation as following Here, Y is chosen as the unit matrix with corresponding rows deleted.
If the number of unknowns, δc, is the same as the number of equations in (29), then the solution of the improvement of control parameter, which is used to eliminate the predictive terminal state deviation, can be obtained as If the number of unknowns is greater than the number of equations in (29), however, the solution can be obtained by formulating additional cost functions that can be minimized or maximized. It is noted that the upper bound of control vector in (3) can be represented as Because φ j (t k ) are fixed, we can obtain the minimization of u (t k ) via minimizing c j . Then, such a cost function of the form can be chosen as where, R is the weighting matrix which needs to be chosen carefully by the control designer. Thus, (29) and (32) formulate an appropriate constrained static optimization problem. According to the static optimization theory [8], the augmented cost function is represented as Using the KKT conditions, it leads to Solving for δc from (35) leads to Adding (36) into (34), Lagrange multiplier λ can be expressed as Substituting λ into (36), the updated control parameter, c, is obtained as Finally, the updated control is represented as It can be observed that the update control parameter of LPMPSC is a closed form expression, which is one of the key reasons responsible for the computational efficiency of the proposed algorithm. In comparison with previous LPMPC, the number of updated variables is significantly reduced. Thus, the computational efficiency of LPMPSC is further improved over LPMPC. And the control function is certified to be a smooth function because it is a linear combination of smooth basis functions φ j (t).

C. IMPLEMENTION OF LPMPSC
The LPMPSC aims to solve the nonlinear optimal control problem with hard terminal constraints by parameterizing the control variables and solving the corresponding transformed linear optimal control problem. The LPMPSC relies on an initial control guess to start the algorithm and iterates control until it converges. The necessary steps for accomplishing the LPMPSC is included in the flow diagram is shown in Fig. 3 and described in detail as follows.
1) Initialization: set initial simulation parameters and select a proper initial control guess. The initial control parameter, c p , is obtained via fitting the initial control guess history.
2) Use the current control iterations to propagate state dynamics, and obtain the trajectory information and the deviation of terminal state, δx f .
3) Judge on the satisfaction of the deviation of terminal state, δx f , tolerance. If δx f is within the specified tolerance value, then go to next step, if not, go to step 6.
4) Update the current control U 0 . 5) Apply the current control U 0 to the real plant and go to step 2. 6) Compute matrix K c . Update the control parameter via (38) and obtain the updated control, u, by (39). Then, go to step 2.
These steps are repeated until the solution for nonlinear optimal control problem meets the terminal constraints. This improved method LPMPSC is easy to implement and remains the character that few points will achieve high accuracy. In addition, it offers further advantages over the previous methods such as 1) significantly decreasing the number of updated variables, 2) further computational efficiency, and 3) guaranteeing smoother control variable by enforcing the spread of control through the whole available time to go.

III. APPLICATION TO TERMINAL GUIDANCE OF INTERCEPTOR
In this section, the proposed method in the preceding section is applied to a three-dimensional interception guidance law with terminal angle constraints. The impact-angle-constraint engagement against a non-maneuvering incoming ballistic missile target by a surface to air interceptor is considered. The objective is to minimize the terminal constraints including the miss distance and terminal angles. Meanwhile, it is better to produce small and smooth lateral acceleration history. It is noted that the biased PN guidance law is used to obtain the initial control guess history for the interceptor.

A. THREE-DIMENSIONAL DYNAMICS AND CONTROL PARAMETER
The standard three-dimensional dynamics of interceptor in plane geodetic coordinate system can be given aṡ It is noted that a first-order autopilot lag is considered to simulate the flight control system. Where, a zc and a yc are lateral acceleration commands in z and y directions respectively, τ is the time constant of the first-order autopilot lag. The realized lateral accelerations in z and y directions are a z and a y respectively. It should be noted that it is not necessary to consider the first-order autopilot lag in the process of calculating the updated lateral acceleration commands. The acceleration vector is perpendicular to the velocity vector V m . (x m , y m , z m ) gives the instantaneous position of the interceptor in threedimensional space. γ m and ψ m are the flight-path angle and the heading angle respectively. Simulation parameters of the interceptor refer to the reference [31].
The standard point-mass dynamics of target over a flat earth is written asẋ where, V t , γ t and ψ t give the velocity vector, flight-path angle and heading angle of the target respectively. The three-dimensional coordinates of the target are given by x t , y t , z t . Note that the target is a coasting flighting vehicle without power maneuvering. The value of the ballistic cofficient of the target β is given as 900kg/m 2 . Next, the state and control variables of the interceptor are normalized to keep the value of those normalized variables in a similar numerical range. The main purpose of the normalization is to ensure numerical stability for solving linear algebraic equations. The normalized state and control variables of the interceptor are represented as where, the state and control variables with a second subscript n are normalized values, second subscript n shows the normalized values, and the normalizing quantities are given with superscript ( * ). The normalizing values and some key simulation parameters are shown in Table 1. It should be noted that the control variable is considered to be parameterized as quadratic function of time, t. Expressions of the basis function, φ (t), for quadratic parameterization are given as φ 1 (t) = t 2 , φ 2 (t) = t, and φ 3 (t) = 1. It should be noted that the selected quadratic function can fit the control curve of interceptor well. Additionally, there are only 3 coefficients of quadratic function that needed to be improved, which is convenient to improve the computational efficiency.

B. INITIAL CONTROL GUESS HISTORY FOR LPMPSC
The initial control guess history for LPMPSC in this note is obtained by applying the biased PN guidance (BPNG) [5] to intercept a non-maneuvering ballistic missile. Since BPNG can accomplish the terminal angle-constrained engagement. It is necessary to introduce the conventional expression of this guidance law. The lateral acceleration command is given as where, N e gives the proportional navigation constant, V c denotes the closing velocity,λ is the LOS angle rate. Note thatλ b is the bias value of the LOS angle rate that aims to eliminate the LOS angle error. The expression ofλ b is given asλ where, λ f is the final LOS angle corresponing to the impact angle constrained, and t go is the time to go until intercepting. λ 0 andλ 0 are the instantaneous LOS angle and LOS angle rate respectively. The expression of the final LOS angle in pitch and yaw are given by where, γ tf and ψ tf is the terminal flight-path angle and heading angle of target. Hence, we can obtain the final interceptor flight-path angle, γ mf , and heading angle, ψ mf , corresponding to desired impact angle. Then, the final LOS angle in pitch and yaw are obtained via (45). σ denotes the three-dimensional line-of-sight angle between the interceptor and target in the inertial frame, and the line-of-sight rateσ is expressed aṡ where, the relative range between the interceptor and target is presented as r = r 2 x + r 2 y + r 2 z , and r x , r y and r z are its components. Next, line-of-sight rate is converted into the LOS frame as followinġ λ z = sin(ψ)σ x − cos(ψ)σ ẏ λ y = − sin(γ )(cos(ψ)σ x + sin(ψ)σ y ) + cos(γ )σ z (47) The closing velocity is given by The BPNG demands the lateral acceleration command as

IV. SIMULATION RESULTS
In this section, several simulations have been carried out to evaluate the effectiveness of the suboptimal impact angle constrained guidance law in intercepting a high-speed nonmaneuvering ballistic target. In addition, a comparison with other typical methods such as BPNG, MPSP and LPMPC is also provided to demonstrate its superior performance in computational accuracy and efficiency. It is assumed that the velocity of the target is higher than that of interceptor. The initial states such as positions, velocities, flight-path angle, and heading angle of the target and interceptor are listed in Table 2. In order to maximize the impact velocity and penetration of warhead, the impact angle for intercepting a high-speed ballistic target is desired to be 180 deg, which means a head-on engagement. For increasing the observability and lethality of directional warhead, the impact angle is desired to be 165 to 175 deg. Hence, the desired terminal flight-path angle, γ df , and heading angle, ψ df , of the interceptor are chosen to be 53 deg and 40 deg respectively. The three-dimensional flight paths of the interceptor and target during the terminal phase engagement using BPNG, MPSP, LPMPC and LPMPSC guidance are depicted in Fig. 4. It can be observed that the flight paths of the interceptor using LPMPSC, MPSP and LPMPC guidance are substantially coincident. And the curvatures of those trajectories are less than that generated by BPNG. The curvature of a flight path is indicative of the achieved lateral acceleration that the smaller the lateral acceleration, the straighter the flight path, and the smaller the aerodynamic drag. Thus, a flight path with small curvature is advantageous in the context of aerodynamic drag minimization. The curvatures of different flight paths are shaped by achieved lateral acceleration generated by BPNG, MPSP, LPMPC and LPMPSC.   Fig. 6 show the achieved vertical and horizontal acceleration histories generated by those four different guidance laws respectively. It is obvious that the lateral  accelerations in vertical and horizontal directions achieved by BPNG change sharply at the end of the interception, which will lead to a great increase in miss distance. It is apparent that LPMPSC, MPSP and LPMPC demand less lateral acceleration in comparison with BPNG. Moreover, the lateral accelerations achieved by LPMPSC, MPSP and LPMPC change gently throughout the whole flight, which leads to smoother trajectory and higher accuracy for terminal constraints including miss distance, and flight-path angle and heading angle. Furthermore, in comparison with LPMPC and MPSP, the lateral acceleration histories obtained by LPMPSC are much smoother. It should be noted that the lateral acceleration histories obtained by LPMPC and MPSP guidance are both optimal, and therefore they coincide well. The suboptimal solution of lateral acceleration achieved by LPMPSC is very close to the optimal solution. Fig. 7 and Fig. 8 show the flight-path angle and heading angle histories of the interceptor during the engagement respectively. It is clear that the terminal angle constraints are all satisfied by those four guidance laws (BPNG, MPSP, LPMPC and LPMPSC). Obviously, the terminal heading  angles for MPSP, LPMPC and LPMPSC are much smoother than that of BPNG. The reason is that BPN will generate the sharp lateral acceleration at the end of interception.
In order to study the computational efficiency and accuracy of LPMPSC guidance law, a comparison among those four guidance laws is provided in Table 3. Table 3 includes the terminal flight-path angle and heading angle, miss distance, the sum of square terminal errors, dYN 2 , the mean computing time, and iterations. It is obvious that LPMPSC, MPSP and LPMPC all produce more precise results in comparison with BPNG, and LPMPSC achieves less miss distance and more accurate terminal angles in comparison with the other three guidance laws. It is noted that the sum of square terminal errors for LPMPSC stands at a lower level in comparison with the other guidance laws.
It should be noted that the time step of prediction simulation for MPSP, LPMPC and LPMPSC is 50 ms. The number of collocation points for MPSP is 121, and the collocation points for LPMPC and LPMPSC are 6 LG points. In comparison with LPMPC, the number of updated variables for LPMPSC is decreased from 12 to only 6. The simulation results reveal that it only takes 3 steps for the LPMPSC to arrive at the desired terminal precision. It takes the same few steps as MPSP and LPMPC to reach a high accuracy. Consequently, the proposed method has a fast convergence rate. The mean computing time for updating the control variable in each iteration is only 0.0015 sec for LPMPSC, which is much less than 0.012 sec and 0.004 sec for MPSP and LPMPC respectively. Therefore, it is easy to conclude that the proposed method has superior performance in computational efficiency in comparison with MPSP and LPMPC methods and is more suitable for on-line application.
Finally, some engagement cases with different terminal conditions are carried out to test the applicability of LPMPSC. The initial states of the target and interceptor for those engagement cases are the same as that listed in Table 2. In case 1, the desired terminal flight-path angles, γ df , are selected as 59 deg, 50 deg, 41deg, and 30deg respectively. The engagement with varied terminal flight-path angles are shown in Fig. 9. Fig. 10. shows the flight-path angle histories with legend on desired terminal flight-path angles and achieved corresponding errors. All cases greatly meet the desired terminal contaminants with high accuracy. In case 2, the desired FIGURE 9. Engagement with fixed γ m0 , ψ m0 and varied γ df . terminal heading angles, ψ df , are selected as 61 deg, 54 deg, 45deg, and 37deg respectively. The desired terminal heading angles and obtained corresponding deviations are shown in Fig. 11. Fig. 12. shows the corresponding trajectories. As can FIGURE 11. Heading angle with fixed γ m0 , ψ m0 and varied ψ df . be seen from these figures, the deviations of terminal angles achieved by the proposed method are 0.0209 deg, 0.0403 deg, 0.0766 deg, and 0.0035 deg respectively. It is easy to conclude that LPMPSC performs well in achieving the desired terminal angles and providing steady and smooth flight curves, which has the great potential in being applied on the framework of guidance.

V. CONCLUSION
In this paper, the linear pseudospectral model predictive control method has been further extended by considering a priori parameterization of the control variable as polynomial functions of time or time to go. The presented method is named as the linear pseudospectral model predictive spread control. The advantage of parameterized expression of the control variable is that it only needs to optimize a few coefficients of the polynomial expression instead of optimizing the control variable at all LG points, which means the number of updated variables is much smaller. Hence, the computational efficiency is improved obviously over the LPMPC, which is better for online applications. And the control variable is guaranteed to be smooth, which leads to great input for the inner-loop autopilot. The proposed method is successfully applied to intercept incoming high-speed ballistic missiles with impact angle constraint. Simulation results demonstrate that the proposed guidance law has good performance in meeting hard terminal constraints (both in miss distance and angle constraint) and has higher computational efficiency in comparison with the MPSP and LPMPC.