# Predictive Functional Control-Based Zenith Pass Controller Design for Roll-Pitch Seeker

Authors: Xinchun Wang; Bo Mo; Xuan Li; Shangen Su
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

A roll-pitch seeker has a wide field of view but suffers from a singularity as the sightline coincides with the outer gimbal (OG) axis. In the vicinity of the singularity, the tracking effectiveness is often degraded or even lost due to the high actuation demand on OG, which is known as the zenith pass problem. To solve this problem, this paper first proposes a novel motion model of sightline to predict the singularity in a receding horizon, where the model parameters are identified using a modified recursive least square estimator. And with the singularity predictions as set points, a predictive functional controller is then designed for the OG position control to minimize the tracking error. This novel combination control scheme is validated in MATLAB/Simulink. Simulation results have confirmed that the proposed scheme can significantly mitigate the zenith pass problem and be applied to the real-time tracking process.


---

## Body

## 1. Introduction

Taking into account the cost savings and size reduction, a two-axis gimbaled platform has been widely used in optics seekers of precision-guided missiles. It can be intelligibly used to track the target and to provide the inertial stabilization for the sensor’s pointing vector. There are usually two kinds of structural seeker: elevate-azimuth seeker and roll-pitch seeker, classified by the rotation axis. Between them, the roll-pitch seeker and the one this paper concerned can be designed smarter and of a larger gimbal angle, which attracted much more attention in the past two decades [1].A sketch of the roll-pitch seeker is shown in Figure1. In this configuration, the target detector and inertial sensors are mounted on the inner gimbal (IG), which is constrained to rotate with the pitch axis. And the outer gimbal (OG) rotates with the roll axis. Each gimbal is actuated by a brushless DC motor and angle information is measured by an encoder, represented by γs,ϑs for OG and IG, respectively. With combinations of γs,ϑs, the boresight of the seeker can point to any orientation in space. In practical applications, however, there is a fundamental limitation to this configuration. That is, the occasional situation where the target passes directly the zenith (defined as ϑs=0 degree), the roll axis has to rotate nearly 180 degrees instantly. It is hard to achieve such a large (mathematically infinite) rate and acceleration demand for a physical servo system. As a result, the roll-pitch seeker will be unable to remain pointed at the target and the tracking error is correspondingly magnified through the zenith pass event. Even the target will be lost if the error is too large to exceed the instantaneous field of view of the seeker, as shown in Figure 2. That is the so-called problem of zenith pass for the roll-pitch gimbal tracking system.Figure 1
Roll-pitch gimbaled seeker configuration.Figure 2
Schematic diagram of zenith pass problem (γsc and ϑsc denote the desired frame angles while γs and ϑs are the actual angles achieved).To avoid or mitigate this problem, several approaches have been proposed in the past [2–5]. They can be roughly summarized into three solutions, namely, three-axis design, tilt outer gimbal axis design, and program guidance method. Among them, this first solution and the one most straightforward is to add an extra axis to the stabilization and tracking system. With the three-axis design, the zenith pass problem can be completely overcome by altering the two rotation axes of three in different tracking scenarios. However, this solution can not only increase the size and cost but also bring some other complex issues like coupling of axes’ inertia. The second solution, the tilt outer gimbal design, can be subdivided into two: fixed tilt angle and adjustable tilt angle. The fixed tilt angle design is just to orient the boresight of the seeker so that the zenith singularity lies away from the working area of the system. It is certainly unsuitable as the area is unknown and the seeker demands a field of hemispherical coverage, while the adjustable tilt angle design can well handle this trouble by switching angles in different scenarios, it brings other issues like increasing size, cost, and complexity as the three-axis design. Unlike the previous two solutions to change the hardware structure, the third one, program guidance method, is to tackle the problem by preplanning the boresight trajectory within acceptable error range. In this method, the OG angle response of zenith pass is finally optimized under the actuator constraints, and the singularity is significantly mitigated. However, this method is basically only used in the design of satellite tracking antennas for the geometry trajectory of sightline has to be known a priori in the process. To achieve this application on the roll-pitch seeker tracking system, the trajectory must be predicted optimally online.Inspired by the last solution, this paper proposed an online receding prediction and optimization scheme for the OG position control around the zenith. The scheme is mainly implemented with two components, the zenith reference predictor and the predictive functional controller (PFC). The predictor uses a low-order sightline motion model to predict the zenith trajectory and to generate the tracking angle reference of OG in a short receding horizon, where the model parameters are updated with the newest measurements using the recursive least square method. While the PFC, a fast processing control technique [6], uses a generalized model of inner loop to predict the angle response over the same horizon, combining this with the zenith reference, a quadratic performance index optimum process can be constructed to minimize the tracking error. Emphatically, the optimization process is repeated with the newest reference for online real-time control and the constraints of actuator can be taken into account in optimizing the performance index. The proposed PFC with the zenith reference predictor is validated in MATLAB/Simulink. Simulation results have proved the effectiveness of the proposed control strategy, by tracking the zenith pass trajectory within a smaller error range than linear PI.The total paper is divided into five major sections starting from first, introduction. Section II elaborates the strategy of zenith reference generation. Section III gives a brief overview of PFC scheme and employs it in the gimbal position controller design. Section IV gives MATLAB simulation results and some observations around zenith, and section VI concludes the contribution of the paper.


## 2. Zenith Reference Predictor

### 2.1. Gimbal Motion Equations

Consider a two-axis, roll-pitch gimbal system as depicted in Figure1. In the figure, the roll and pitch gimbals are indicated. The sensor is placed on the pitch gimbal. It detects and evaluates the sightline angle errors εy,εz, which are defined by the angles between sightline and antenna boresight in pitch and yaw axis plane [7]. The sightline angle errors are fed back to move the gimbals with desired angles, and they make the antenna track the target.To guarantee the track is error free, the precise motion equations of the gimbal angels and rates are derived in this part. The definitions of the coordinate systems used in derivation of the angles are summarized in Table1. In the table, Rxγs is the coordinate transform matrix representing the rotation of angle γs about x-axis. The relations among the coordinate systems are illustrated in Figures 1, 3, and 4.Table 1
Definitions of coordinate system.
Frame
Definition
B
(i) Missile body frame (regarded stationary for simplicity)(ii) XB: missile nose direction(iii) YB: right direction, on the yaw plane
S
(i) Seeker antenna heading frame(ii) XS: antenna boresight direction(iii) ZS: coincides with pitch gimbal axis(iv) TSB=RZϑSRXγS,γS,ϑS = seeker roll, pitch gimbal angles
P
(i) Target pointing frame(ii) XP: sightline direction(iii) TPS=RzεzRyεy
DS
(i) Target pointing frame(ii) XDS: sightline direction(iii) ZDS: coincides with pitch gimbal angles(iv) TDSB=RZϑDSRxγDS,ϑDS,γDS = desired roll, pitch gimbal angles of seeker
LOS
(i) LOS frame(ii) XLOS: sightline direction(iii) ZLOS: coincides with the z-axis of body frame(iv) TLOSB=RzλRyq,q,λ = azimuth and elevation angles of LOS based on missile body frameFigure 3
Sensor heading and target pointing frame.Figure 4
Relation between missile body and LOS frames.It should be emphasized that the frame of LOS, defined in the last row of the table, is introduced as an intermediate coordinate system to calculate the desired gimbal anglesγDS,ϑDSfrom the current angles γs,ϑsand the sightline angle errors εy,εz. This may seem redundant and conversely increase the complexity of derivation. But, as will be shown in the next part, the azimuth and elevation angles of the LOS frame are characteristic parameters to predict the sightline trajectory and used to generate the OG angle references around the zenith, such that the less additional complexity could be ignored.In most situations, the azimuth and elevation angles are found by calculating the unit vector representing line of sight, then using the information to recover the angles from the rotation matrix [8]. Here, the sightline vector can be first calculated with angles γs,ϑs,εy,εz, and since the orientation of the vector coincides with that of XP and XLOS axes, the sightline angles q,λ in terms of γs,ϑs,εy,εz can be easily derived:
(1)q=tan−1cosεzcosεysinϑssinγssinεzsinϑs−cosεzcosεycosϑs+sinεzcosϑssinγs−cosεzsinεycosγssinεzsinϑs−cosεzcosεycosϑs,λ=sin−1cosεzcosεysinϑscosγs+sinεzcosϑscosγs+cosεzsinεysinγs.With a similar derivation, the transformation fromq,λ to γDS,ϑDS is given by
(2)γDS=tan−1−sinqtanλ,ϑDS=cos−1cosqcosλ.Take the derivative of Equation (2), with respect to time, we have
(3)γ̇DS=sinqsin2λ+cos2λsin2qλ̇−cosqtanλtan2λ+sin2qq̇,ϑ̇DS=tanλcosqtan2λ+sin2qλ̇+sinqtan2λ+sin2qq̇,where γ̇DS,ϑ̇DS are the desired gimbal rates for tracking. q̇,λ̇˙ are the azimuth and elevation rates of sightline, respectively. When the sightline moves close to the zenith, i.e., q⟶0,λ⟶0, Equation (3) can be simplified with small angle approximation (sinq≈q,tanq≈q,cosq≈1, the same for λ), as shown
(4)γ̇DS=qλ̇−λq̇λ2+q2,ϑ̇DS=qq̇+λλ̇λ2+q2.It can be easily found from Equation (4) that the limit of γ̇DS tends to infinity as q⟶0,λ⟶0, which means when the sightline approaches around the zenith, the outer gimbal has to follow a very large rate demand. It is hard to achieve for a physical servo system. So loss of target is likely to happen in the circumstance, and this is the called zenith pass problem.


### 2.2. Zenith Reference Generation

Predictive functional control requires a set-point trajectory that reflects estimates of the demands for a period of time into the future (the prediction horizon). For a tracking task, this is not generally possible, but the seeker tracking scenario is such that because the crossing rate of sightline as seen from the missile is quite low [9]. Thus, it can be assumed, the azimuth/elevation sequence of sightline follows a stationary time series model in a relatively short period of time. And it intends to adopt a low-order stationary time series model to fit the motion of sightline (here, three-order is selected for use). By identifying the model parameters with past measurements, the required future azimuth/elevation can be predicted through extrapolating the motion model into the prediction horizon.Assuming that the time series model about azimuth/elevation of sightline (based on frame B) could be written as [10]
(5)yk=ϕTk−1θ,with
(6)ϕTk−1=yk−1,yk−2,⋯,yk−n,θ=a1,a2,⋯,anT.where the azimuth/elevation is modeled by a linear combination of its past values. θ is the unknown combination coefficients. ϕTk−1 represents the most recent n samples of sightline, of which the value is usually unmeasurable but calculated with the measurements of γs,ϑs,εy,εz through Equation (1).A recursive least-squares (RLS) estimator usually with some modification is the standard recursive estimation technique for determining the parameters in such a system [11]. Here, for tracking the time-varying azimuths, we employ two windowing techniques in RLS: one is the exponentially weighted window to emphasize recent data and the other is the sliding window to discard old data [12]. In this way, the most recent changes in the relative motion of the missile and target are sufficiently taken into account, and the past unfit motion is ignored. The corresponding cost function is defined as follows:
(7)Ek=∑i=k−N+1kwk−iyi−ϕTi−1θ^k2,where θ^k is the estimated vector of θ. N and w are a window size and a weighting factor, respectively. Differentiating Ek with respect to θ^k and setting the result to zero renders
(8)θ^k=Φ−1kΨk,where
(9)Φk=∑i=k−N+1kwk−iϕi−1ϕTi−1,Ψk=∑i=k−N+1kwk−iyiϕi−1.For iterative calculation of Equation (8), Φk+1 is expressed as
(10)Φk+1=∑i=k−N+2k+1wk+1−iϕi−1ϕTi−1=wΦk+ϕkϕTk−wNϕk−NϕTk−N.Similarly,(11)Ψk+1=wΨk+yk+1ϕk−wNyk−N+1ϕk−N,where Ψk can be replaced with Φk·θ^k.Substituting Equations (10) and (11) into Equation (8), the estimate θ^k+1 is then recursively computed as
(12)θ^k+1=θ^k+Φ−1k+1·ϕke0k+1−wNϕk−Nelk+1,where
(13)e0k+1=yk+1−ϕTk+1θ^k,elk+1=yk−N+1−ϕTk−Nθ^k.As for the calculation ofΦ−1k+1, we define an intermediate matrix:
(14)Mk=wΦk+ϕkϕTk.Then,(15)Φk+1=Mk−wNϕk−NϕTk−N.Using the matrix inversion lemma, Equations (14) and (15) become
(16)M−1k=1wΦ−1k−Φ−1kϕkϕTkΦ−1kw+ϕTkΦ−1kϕk,Φ−1k+1=M−1k+M−1kϕk−NϕTk−NM−1kw−N−ϕTk−NM−1kϕk−N.By means of Equations (12) and (16), the estimated combination parameters are recursively calculated. And the rolling azimuth/elevation predictions over the short prediction horizon are then obtained as
(17)q^k+i=ϕ^qTk+i−1θ^k,λ^k+i=ϕ^λTk+i−1θ^k,i=1,2,⋯,P,where
(18)ϕ^qTk+i−1=q^k+i−1,⋯,qk+i−N,ϕ^λTk+i−1=λ^k+i−1,⋯,λk+i−N.Substitute the predictions into Equation (2), the receding OG angle reference at time instant k is finally generated.


## 2.1. Gimbal Motion Equations

Consider a two-axis, roll-pitch gimbal system as depicted in Figure1. In the figure, the roll and pitch gimbals are indicated. The sensor is placed on the pitch gimbal. It detects and evaluates the sightline angle errors εy,εz, which are defined by the angles between sightline and antenna boresight in pitch and yaw axis plane [7]. The sightline angle errors are fed back to move the gimbals with desired angles, and they make the antenna track the target.To guarantee the track is error free, the precise motion equations of the gimbal angels and rates are derived in this part. The definitions of the coordinate systems used in derivation of the angles are summarized in Table1. In the table, Rxγs is the coordinate transform matrix representing the rotation of angle γs about x-axis. The relations among the coordinate systems are illustrated in Figures 1, 3, and 4.Table 1
Definitions of coordinate system.
Frame
Definition
B
(i) Missile body frame (regarded stationary for simplicity)(ii) XB: missile nose direction(iii) YB: right direction, on the yaw plane
S
(i) Seeker antenna heading frame(ii) XS: antenna boresight direction(iii) ZS: coincides with pitch gimbal axis(iv) TSB=RZϑSRXγS,γS,ϑS = seeker roll, pitch gimbal angles
P
(i) Target pointing frame(ii) XP: sightline direction(iii) TPS=RzεzRyεy
DS
(i) Target pointing frame(ii) XDS: sightline direction(iii) ZDS: coincides with pitch gimbal angles(iv) TDSB=RZϑDSRxγDS,ϑDS,γDS = desired roll, pitch gimbal angles of seeker
LOS
(i) LOS frame(ii) XLOS: sightline direction(iii) ZLOS: coincides with the z-axis of body frame(iv) TLOSB=RzλRyq,q,λ = azimuth and elevation angles of LOS based on missile body frameFigure 3
Sensor heading and target pointing frame.Figure 4
Relation between missile body and LOS frames.It should be emphasized that the frame of LOS, defined in the last row of the table, is introduced as an intermediate coordinate system to calculate the desired gimbal anglesγDS,ϑDSfrom the current angles γs,ϑsand the sightline angle errors εy,εz. This may seem redundant and conversely increase the complexity of derivation. But, as will be shown in the next part, the azimuth and elevation angles of the LOS frame are characteristic parameters to predict the sightline trajectory and used to generate the OG angle references around the zenith, such that the less additional complexity could be ignored.In most situations, the azimuth and elevation angles are found by calculating the unit vector representing line of sight, then using the information to recover the angles from the rotation matrix [8]. Here, the sightline vector can be first calculated with angles γs,ϑs,εy,εz, and since the orientation of the vector coincides with that of XP and XLOS axes, the sightline angles q,λ in terms of γs,ϑs,εy,εz can be easily derived:
(1)q=tan−1cosεzcosεysinϑssinγssinεzsinϑs−cosεzcosεycosϑs+sinεzcosϑssinγs−cosεzsinεycosγssinεzsinϑs−cosεzcosεycosϑs,λ=sin−1cosεzcosεysinϑscosγs+sinεzcosϑscosγs+cosεzsinεysinγs.With a similar derivation, the transformation fromq,λ to γDS,ϑDS is given by
(2)γDS=tan−1−sinqtanλ,ϑDS=cos−1cosqcosλ.Take the derivative of Equation (2), with respect to time, we have
(3)γ̇DS=sinqsin2λ+cos2λsin2qλ̇−cosqtanλtan2λ+sin2qq̇,ϑ̇DS=tanλcosqtan2λ+sin2qλ̇+sinqtan2λ+sin2qq̇,where γ̇DS,ϑ̇DS are the desired gimbal rates for tracking. q̇,λ̇˙ are the azimuth and elevation rates of sightline, respectively. When the sightline moves close to the zenith, i.e., q⟶0,λ⟶0, Equation (3) can be simplified with small angle approximation (sinq≈q,tanq≈q,cosq≈1, the same for λ), as shown
(4)γ̇DS=qλ̇−λq̇λ2+q2,ϑ̇DS=qq̇+λλ̇λ2+q2.It can be easily found from Equation (4) that the limit of γ̇DS tends to infinity as q⟶0,λ⟶0, which means when the sightline approaches around the zenith, the outer gimbal has to follow a very large rate demand. It is hard to achieve for a physical servo system. So loss of target is likely to happen in the circumstance, and this is the called zenith pass problem.


## 2.2. Zenith Reference Generation

Predictive functional control requires a set-point trajectory that reflects estimates of the demands for a period of time into the future (the prediction horizon). For a tracking task, this is not generally possible, but the seeker tracking scenario is such that because the crossing rate of sightline as seen from the missile is quite low [9]. Thus, it can be assumed, the azimuth/elevation sequence of sightline follows a stationary time series model in a relatively short period of time. And it intends to adopt a low-order stationary time series model to fit the motion of sightline (here, three-order is selected for use). By identifying the model parameters with past measurements, the required future azimuth/elevation can be predicted through extrapolating the motion model into the prediction horizon.Assuming that the time series model about azimuth/elevation of sightline (based on frame B) could be written as [10]
(5)yk=ϕTk−1θ,with
(6)ϕTk−1=yk−1,yk−2,⋯,yk−n,θ=a1,a2,⋯,anT.where the azimuth/elevation is modeled by a linear combination of its past values. θ is the unknown combination coefficients. ϕTk−1 represents the most recent n samples of sightline, of which the value is usually unmeasurable but calculated with the measurements of γs,ϑs,εy,εz through Equation (1).A recursive least-squares (RLS) estimator usually with some modification is the standard recursive estimation technique for determining the parameters in such a system [11]. Here, for tracking the time-varying azimuths, we employ two windowing techniques in RLS: one is the exponentially weighted window to emphasize recent data and the other is the sliding window to discard old data [12]. In this way, the most recent changes in the relative motion of the missile and target are sufficiently taken into account, and the past unfit motion is ignored. The corresponding cost function is defined as follows:
(7)Ek=∑i=k−N+1kwk−iyi−ϕTi−1θ^k2,where θ^k is the estimated vector of θ. N and w are a window size and a weighting factor, respectively. Differentiating Ek with respect to θ^k and setting the result to zero renders
(8)θ^k=Φ−1kΨk,where
(9)Φk=∑i=k−N+1kwk−iϕi−1ϕTi−1,Ψk=∑i=k−N+1kwk−iyiϕi−1.For iterative calculation of Equation (8), Φk+1 is expressed as
(10)Φk+1=∑i=k−N+2k+1wk+1−iϕi−1ϕTi−1=wΦk+ϕkϕTk−wNϕk−NϕTk−N.Similarly,(11)Ψk+1=wΨk+yk+1ϕk−wNyk−N+1ϕk−N,where Ψk can be replaced with Φk·θ^k.Substituting Equations (10) and (11) into Equation (8), the estimate θ^k+1 is then recursively computed as
(12)θ^k+1=θ^k+Φ−1k+1·ϕke0k+1−wNϕk−Nelk+1,where
(13)e0k+1=yk+1−ϕTk+1θ^k,elk+1=yk−N+1−ϕTk−Nθ^k.As for the calculation ofΦ−1k+1, we define an intermediate matrix:
(14)Mk=wΦk+ϕkϕTk.Then,(15)Φk+1=Mk−wNϕk−NϕTk−N.Using the matrix inversion lemma, Equations (14) and (15) become
(16)M−1k=1wΦ−1k−Φ−1kϕkϕTkΦ−1kw+ϕTkΦ−1kϕk,Φ−1k+1=M−1k+M−1kϕk−NϕTk−NM−1kw−N−ϕTk−NM−1kϕk−N.By means of Equations (12) and (16), the estimated combination parameters are recursively calculated. And the rolling azimuth/elevation predictions over the short prediction horizon are then obtained as
(17)q^k+i=ϕ^qTk+i−1θ^k,λ^k+i=ϕ^λTk+i−1θ^k,i=1,2,⋯,P,where
(18)ϕ^qTk+i−1=q^k+i−1,⋯,qk+i−N,ϕ^λTk+i−1=λ^k+i−1,⋯,λk+i−N.Substitute the predictions into Equation (2), the receding OG angle reference at time instant k is finally generated.


## 3. PFC Position Controller

Most high-performance servo systems adopt a three-loop control strategy: an inner current loop for an accurate torque, a speed loop for platform stability, and a position loop for angle tracking [13]. In view of this, we proposed a cascade scheme for combining PFC with PID in the servomechanism control. That is, the inner current loop and speed loop adopt in the classical PID control method to guarantee a good interference rejection, while the position loop employs a novel PFC controller to minimize the tracking error with respect to the zenith reference. The minimization is specifically achieved through the controller to optimize the rate demand for the speed loop. To illustrate this process, this section first outlines the design procedures of PFC. And a mathematical model of the internal loop including actuator is then derived as the generalized control plant for application. The actuator constraint is also taken into account by transferring to the generalized model in this process. Finally, completing appropriate parameters, a sightline position controller is designed on the proposed technique.


### 3.1. Principle of Predictive Functional Control

In contrast to other predictive control techniques, PFC can significantly reduce the online computation time, which makes it possible to apply the method in real-time control of fast response system. Within the scope of model predictive control (MPC), it is mainly composed of the following parts [14]:


#### 3.1.1. Base Function

The future controlled variable is structured as a linear combination of several prior known functions:(19)uk+i=∑j=1nBμjfji,i=0,1,2,⋯,P−1,where nB is the number of base functions, P is the length of prediction horizon, and μj is the weight coefficient corresponding to the base function fj. The selection of base functions is determined by the set-point trajectory and the nature of the process. Generally, canonical functions are used, e.g., step, ramp, or exponential function [15]. Depending on the selected base functions, the output response of the object can be calculated offline a priori. And only the weighting coefficients need to be found online. For the gimbal servo control system, we will use the polynomial type as
(20)fji=ij−1j=1,2,⋯,nB,j=1,2,⋯,nB.


#### 3.1.2. Prediction Model

A linear numerical model called internal model is required for PFC to predict the future outputs of plant. It can be expressed in state space as(21)Xmk+1=AmXmk+Bmuk,ymk=CmXmk,where Am,Bm, and Cm are, respectively, coefficient vectors or matrixes of the state equation. Xm is the state vector, u is the input scalar, and ym is the model output.Through recursion, the model predictive output at sampling timek+i is divided into two parts:
(22)ymk+i=yufk+i+yfk+ii=1,2,⋯,P,with
(23)yufk+i=CmAmiXmk,yfk+i=∑j=0i−1CmAmi−1−jBmuk+j,where yuf is the free output (the past output) response to u=0 and yfis the forced output (the future output) response to the control variable given by Equation (19). It can be abbreviated as
(24)yfk+i=μTgi,where μ=μ1,⋯,μnBT,gi=g1i,⋯gnBiTand gji=∑l=0i−1CmAmi−1−lBmfjl.


#### 3.1.3. Error Compensation

Due to the model mismatch and the impact of noises, there is always an error between the model and the practical outputs. In PFC, the errors over the prediction horizon are predicted and then counted as the feedforward into the reference trajectory for compensation. For a stable process, the future error can be expressed as(25)ek+i=ek=yk−ymk,where ymk,yk are the model and practical outputs at the sampling tim k, respectively. Hence, the practical output prediction ypk+i can be expressed as
(26)ypk+i=ymk+i+ek+i


#### 3.1.4. Reference Trajectory

For preventing drastic changes and overshoot, the reference trajectory is constructed to provide a smooth transition towards the future set points within a certain prediction horizon. The new expected path defines the behavior of the closed-loop system. In its construction, an exponential function is often employed for some excellent mathematical characteristics [6]. The reference trajectory is given by
(27)yrk+i=ck+1−αick−yk,where yr represents the reference trajectory, c is the set value, and α=e−Ts/Tr is the attenuation coefficient, determined by the sampling period Ts and the reference trajectory time constant Tr.


#### 3.1.5. Performance Index

A sum of the errors between the predicted actual output and the reference trajectory is selected as the performance index. It is defined as follows:(28)Jk=∑i=1Sypk+hi−yrk+hi2,where S represents the number of coincident points, hi is the coincident point on the prediction horizon. In the absence of constraints, the minimization of Jk with respect to the coefficient leads to the optimal control sequence and only the first term is effectively applied to control, that is,
(29)uk=k0ck−yk+kmXmk+VTCk,with
(30)Ck=ck+h1−ck,⋯,ck+hs−ckT.Apparently, the control variableuk at the time instant k is comprised of three parts, namely, the position tracking error term, the model compensation term, and the future planning term. All gains of the three terms, k0,km,V can be calculated off-line as follows:
(31)k0=VT1−αh11−αh2⋮1−αh3,km=−VTCmTAmh1−ICmTAmh2−I⋮CmTAmhs−I,V=RTf10f20⋮fnB0,where
(32)R=∑j=1SghjghjT−1g1ig2i⋯gnBi


#### 3.1.6. Constraints

Constrained optimization problems are usually solved by quadratic programming, which is time consuming and unsuitable for online tracking processes. While in PFC, the input constraints are directly taken into account by passing the projections of the input produced by the regulator through appropriate limiters. That is,(33)uk=maxufk,umax,ufk>0,minufk,umin,ufk<0,where ufk is obtained through Equation (3) without considering constraints, umax,umin are the constraints imposed on the input. By considering the future scenario for the current sample point only and to repeat at each future sampling point, this strategy turns out to be quite sufficient for most applications.


### 3.2. Application to the Position Controller

Each gimbal is driven by a brushless dc motor via a pulse-width modulation. Beyond the elastic gears, the rigid-body gimbal dynamics can be modeled in terms of the relevant angular position variableθs in the following manner:
(34)Jeqθ¨s+Bθ̇s=Tm,where Jeq,B are the moment of inertia of the gimbal about its rotation axis and viscous friction coefficient, respectively. Tm stands for the control force actuated by the torque motor, which is modeled with armature resistance Ra, inductance La, and back-electromotive ke through simple circuit equivalence expressions as follows [13]:
(35)Raiat+Ladiatdt=uat−keθ̇s,Tm=ktiat,ke=kt.Applying the model to the proposed control strategy, an OG position controller block diagram is constructed, as shown in Figure5. In the figure, the PFC controller works after the latest zenith reference, and behind the controller is the internal loop, which is taken as a whole object Gs for PFC control. Neglecting the non-linear links, the transfer function of G(s) is expressed by:
(36)Gs=Gωs·1s,Gωs=b0s2+b1s+b2a0s4+a1s3+a2s2+a3s+a4,with
(37)a0=JeqLa,a1=JeqRa+KPi+LaB,a2=JeqRa+KPi+LaB,a3=KIiB+ktKIiKPω+KIωKPi,a4=ktKIiKIω,b0=ktKPiKPω,b1=ktKIiKPω+KIωKPi,b2=ktKIiKIω,where KP,KI are proportional constant and integral constant of PI controller and superscripts i and ω denote the current loop and speed loop, respectively. Their numerical values as well as other model parameters of an actual system are summarized in Table 2. Gωs stands for the transfer function of speed loop. For such an asymptotically stable high-order object, Gωs can be equivalent to a second-order model through fitting and simplifying. Thus, the generalized controlled object is simplified to a three-order form Gms as
(38)Gms=3.665×106s3+1.58×104s2+3.569×106sFigure 5
OG position controller block diagram.Table 2
Model parameters.
Model parameters
Numerical values
Jeq,B
0.024kg·m2,0Ns/rad
La,Ra
0.6mH,0.8Ω
kt,ke
0.26N/A,0.26Ns/rad
KPω,KIω
0.4,0.1
KPi,KIi
8.4,400
Ts,Tsω,Tsi
0.02s,0.001s,0.0001sDiscretizing Equation (38) at each sampling instant Ts, the internal model for prediction Am,Bm,Cm is obtained.Another procedure that must be specially designed is the constraints imposed on the control input of the PFC controller. The control input now is the desired motor speed. For its constrained optimization, the physical actuator constraints including speed and torque limits should be transferred to the position controller. Specifically, the speed limit of the actuator is directly imposed on the input through an appropriate limiter as mentioned before (the maximum of the limiter here is selected as 8 rad/s for infrared image processing limit), while the torque limit (i.e., the current limit) is respected by converting into instantaneous speed limit with the discretization of transfer function fromi to ω, as shown:
(39)θ̇mn=1JeqktimaxTs−Bθ̇sn−1+θ̇sn−1,where the current is assumed to reach the constraint over the sampling time. θ̇s is the actual rate response, and θ̇m is the instantaneous constraint, which is obviously variable because of θ̇s. By calculating and preloading θ̇mon the PFC controller at each sample point, it is ensured that the speed demand provided to the internal loop is appropriate and feasible to improve the regulation performance. And the internal model derived above is also guaranteed valid in this way.The other parameters for implementing the PFC controller are chosen after experiments asTr=1ms,nB=2,P=20,S=20. With these, the tracking law of OG around zenith is deduced based on the principle discussed above.


## 3.1. Principle of Predictive Functional Control

In contrast to other predictive control techniques, PFC can significantly reduce the online computation time, which makes it possible to apply the method in real-time control of fast response system. Within the scope of model predictive control (MPC), it is mainly composed of the following parts [14]:


### 3.1.1. Base Function

The future controlled variable is structured as a linear combination of several prior known functions:(19)uk+i=∑j=1nBμjfji,i=0,1,2,⋯,P−1,where nB is the number of base functions, P is the length of prediction horizon, and μj is the weight coefficient corresponding to the base function fj. The selection of base functions is determined by the set-point trajectory and the nature of the process. Generally, canonical functions are used, e.g., step, ramp, or exponential function [15]. Depending on the selected base functions, the output response of the object can be calculated offline a priori. And only the weighting coefficients need to be found online. For the gimbal servo control system, we will use the polynomial type as
(20)fji=ij−1j=1,2,⋯,nB,j=1,2,⋯,nB.


### 3.1.2. Prediction Model

A linear numerical model called internal model is required for PFC to predict the future outputs of plant. It can be expressed in state space as(21)Xmk+1=AmXmk+Bmuk,ymk=CmXmk,where Am,Bm, and Cm are, respectively, coefficient vectors or matrixes of the state equation. Xm is the state vector, u is the input scalar, and ym is the model output.Through recursion, the model predictive output at sampling timek+i is divided into two parts:
(22)ymk+i=yufk+i+yfk+ii=1,2,⋯,P,with
(23)yufk+i=CmAmiXmk,yfk+i=∑j=0i−1CmAmi−1−jBmuk+j,where yuf is the free output (the past output) response to u=0 and yfis the forced output (the future output) response to the control variable given by Equation (19). It can be abbreviated as
(24)yfk+i=μTgi,where μ=μ1,⋯,μnBT,gi=g1i,⋯gnBiTand gji=∑l=0i−1CmAmi−1−lBmfjl.


### 3.1.3. Error Compensation

Due to the model mismatch and the impact of noises, there is always an error between the model and the practical outputs. In PFC, the errors over the prediction horizon are predicted and then counted as the feedforward into the reference trajectory for compensation. For a stable process, the future error can be expressed as(25)ek+i=ek=yk−ymk,where ymk,yk are the model and practical outputs at the sampling tim k, respectively. Hence, the practical output prediction ypk+i can be expressed as
(26)ypk+i=ymk+i+ek+i


### 3.1.4. Reference Trajectory

For preventing drastic changes and overshoot, the reference trajectory is constructed to provide a smooth transition towards the future set points within a certain prediction horizon. The new expected path defines the behavior of the closed-loop system. In its construction, an exponential function is often employed for some excellent mathematical characteristics [6]. The reference trajectory is given by
(27)yrk+i=ck+1−αick−yk,where yr represents the reference trajectory, c is the set value, and α=e−Ts/Tr is the attenuation coefficient, determined by the sampling period Ts and the reference trajectory time constant Tr.


### 3.1.5. Performance Index

A sum of the errors between the predicted actual output and the reference trajectory is selected as the performance index. It is defined as follows:(28)Jk=∑i=1Sypk+hi−yrk+hi2,where S represents the number of coincident points, hi is the coincident point on the prediction horizon. In the absence of constraints, the minimization of Jk with respect to the coefficient leads to the optimal control sequence and only the first term is effectively applied to control, that is,
(29)uk=k0ck−yk+kmXmk+VTCk,with
(30)Ck=ck+h1−ck,⋯,ck+hs−ckT.Apparently, the control variableuk at the time instant k is comprised of three parts, namely, the position tracking error term, the model compensation term, and the future planning term. All gains of the three terms, k0,km,V can be calculated off-line as follows:
(31)k0=VT1−αh11−αh2⋮1−αh3,km=−VTCmTAmh1−ICmTAmh2−I⋮CmTAmhs−I,V=RTf10f20⋮fnB0,where
(32)R=∑j=1SghjghjT−1g1ig2i⋯gnBi


### 3.1.6. Constraints

Constrained optimization problems are usually solved by quadratic programming, which is time consuming and unsuitable for online tracking processes. While in PFC, the input constraints are directly taken into account by passing the projections of the input produced by the regulator through appropriate limiters. That is,(33)uk=maxufk,umax,ufk>0,minufk,umin,ufk<0,where ufk is obtained through Equation (3) without considering constraints, umax,umin are the constraints imposed on the input. By considering the future scenario for the current sample point only and to repeat at each future sampling point, this strategy turns out to be quite sufficient for most applications.


## 3.1.1. Base Function

The future controlled variable is structured as a linear combination of several prior known functions:(19)uk+i=∑j=1nBμjfji,i=0,1,2,⋯,P−1,where nB is the number of base functions, P is the length of prediction horizon, and μj is the weight coefficient corresponding to the base function fj. The selection of base functions is determined by the set-point trajectory and the nature of the process. Generally, canonical functions are used, e.g., step, ramp, or exponential function [15]. Depending on the selected base functions, the output response of the object can be calculated offline a priori. And only the weighting coefficients need to be found online. For the gimbal servo control system, we will use the polynomial type as
(20)fji=ij−1j=1,2,⋯,nB,j=1,2,⋯,nB.


## 3.1.2. Prediction Model

A linear numerical model called internal model is required for PFC to predict the future outputs of plant. It can be expressed in state space as(21)Xmk+1=AmXmk+Bmuk,ymk=CmXmk,where Am,Bm, and Cm are, respectively, coefficient vectors or matrixes of the state equation. Xm is the state vector, u is the input scalar, and ym is the model output.Through recursion, the model predictive output at sampling timek+i is divided into two parts:
(22)ymk+i=yufk+i+yfk+ii=1,2,⋯,P,with
(23)yufk+i=CmAmiXmk,yfk+i=∑j=0i−1CmAmi−1−jBmuk+j,where yuf is the free output (the past output) response to u=0 and yfis the forced output (the future output) response to the control variable given by Equation (19). It can be abbreviated as
(24)yfk+i=μTgi,where μ=μ1,⋯,μnBT,gi=g1i,⋯gnBiTand gji=∑l=0i−1CmAmi−1−lBmfjl.


## 3.1.3. Error Compensation

Due to the model mismatch and the impact of noises, there is always an error between the model and the practical outputs. In PFC, the errors over the prediction horizon are predicted and then counted as the feedforward into the reference trajectory for compensation. For a stable process, the future error can be expressed as(25)ek+i=ek=yk−ymk,where ymk,yk are the model and practical outputs at the sampling tim k, respectively. Hence, the practical output prediction ypk+i can be expressed as
(26)ypk+i=ymk+i+ek+i


## 3.1.4. Reference Trajectory

For preventing drastic changes and overshoot, the reference trajectory is constructed to provide a smooth transition towards the future set points within a certain prediction horizon. The new expected path defines the behavior of the closed-loop system. In its construction, an exponential function is often employed for some excellent mathematical characteristics [6]. The reference trajectory is given by
(27)yrk+i=ck+1−αick−yk,where yr represents the reference trajectory, c is the set value, and α=e−Ts/Tr is the attenuation coefficient, determined by the sampling period Ts and the reference trajectory time constant Tr.


## 3.1.5. Performance Index

A sum of the errors between the predicted actual output and the reference trajectory is selected as the performance index. It is defined as follows:(28)Jk=∑i=1Sypk+hi−yrk+hi2,where S represents the number of coincident points, hi is the coincident point on the prediction horizon. In the absence of constraints, the minimization of Jk with respect to the coefficient leads to the optimal control sequence and only the first term is effectively applied to control, that is,
(29)uk=k0ck−yk+kmXmk+VTCk,with
(30)Ck=ck+h1−ck,⋯,ck+hs−ckT.Apparently, the control variableuk at the time instant k is comprised of three parts, namely, the position tracking error term, the model compensation term, and the future planning term. All gains of the three terms, k0,km,V can be calculated off-line as follows:
(31)k0=VT1−αh11−αh2⋮1−αh3,km=−VTCmTAmh1−ICmTAmh2−I⋮CmTAmhs−I,V=RTf10f20⋮fnB0,where
(32)R=∑j=1SghjghjT−1g1ig2i⋯gnBi


## 3.1.6. Constraints

Constrained optimization problems are usually solved by quadratic programming, which is time consuming and unsuitable for online tracking processes. While in PFC, the input constraints are directly taken into account by passing the projections of the input produced by the regulator through appropriate limiters. That is,(33)uk=maxufk,umax,ufk>0,minufk,umin,ufk<0,where ufk is obtained through Equation (3) without considering constraints, umax,umin are the constraints imposed on the input. By considering the future scenario for the current sample point only and to repeat at each future sampling point, this strategy turns out to be quite sufficient for most applications.


## 3.2. Application to the Position Controller

Each gimbal is driven by a brushless dc motor via a pulse-width modulation. Beyond the elastic gears, the rigid-body gimbal dynamics can be modeled in terms of the relevant angular position variableθs in the following manner:
(34)Jeqθ¨s+Bθ̇s=Tm,where Jeq,B are the moment of inertia of the gimbal about its rotation axis and viscous friction coefficient, respectively. Tm stands for the control force actuated by the torque motor, which is modeled with armature resistance Ra, inductance La, and back-electromotive ke through simple circuit equivalence expressions as follows [13]:
(35)Raiat+Ladiatdt=uat−keθ̇s,Tm=ktiat,ke=kt.Applying the model to the proposed control strategy, an OG position controller block diagram is constructed, as shown in Figure5. In the figure, the PFC controller works after the latest zenith reference, and behind the controller is the internal loop, which is taken as a whole object Gs for PFC control. Neglecting the non-linear links, the transfer function of G(s) is expressed by:
(36)Gs=Gωs·1s,Gωs=b0s2+b1s+b2a0s4+a1s3+a2s2+a3s+a4,with
(37)a0=JeqLa,a1=JeqRa+KPi+LaB,a2=JeqRa+KPi+LaB,a3=KIiB+ktKIiKPω+KIωKPi,a4=ktKIiKIω,b0=ktKPiKPω,b1=ktKIiKPω+KIωKPi,b2=ktKIiKIω,where KP,KI are proportional constant and integral constant of PI controller and superscripts i and ω denote the current loop and speed loop, respectively. Their numerical values as well as other model parameters of an actual system are summarized in Table 2. Gωs stands for the transfer function of speed loop. For such an asymptotically stable high-order object, Gωs can be equivalent to a second-order model through fitting and simplifying. Thus, the generalized controlled object is simplified to a three-order form Gms as
(38)Gms=3.665×106s3+1.58×104s2+3.569×106sFigure 5
OG position controller block diagram.Table 2
Model parameters.
Model parameters
Numerical values
Jeq,B
0.024kg·m2,0Ns/rad
La,Ra
0.6mH,0.8Ω
kt,ke
0.26N/A,0.26Ns/rad
KPω,KIω
0.4,0.1
KPi,KIi
8.4,400
Ts,Tsω,Tsi
0.02s,0.001s,0.0001sDiscretizing Equation (38) at each sampling instant Ts, the internal model for prediction Am,Bm,Cm is obtained.Another procedure that must be specially designed is the constraints imposed on the control input of the PFC controller. The control input now is the desired motor speed. For its constrained optimization, the physical actuator constraints including speed and torque limits should be transferred to the position controller. Specifically, the speed limit of the actuator is directly imposed on the input through an appropriate limiter as mentioned before (the maximum of the limiter here is selected as 8 rad/s for infrared image processing limit), while the torque limit (i.e., the current limit) is respected by converting into instantaneous speed limit with the discretization of transfer function fromi to ω, as shown:
(39)θ̇mn=1JeqktimaxTs−Bθ̇sn−1+θ̇sn−1,where the current is assumed to reach the constraint over the sampling time. θ̇s is the actual rate response, and θ̇m is the instantaneous constraint, which is obviously variable because of θ̇s. By calculating and preloading θ̇mon the PFC controller at each sample point, it is ensured that the speed demand provided to the internal loop is appropriate and feasible to improve the regulation performance. And the internal model derived above is also guaranteed valid in this way.The other parameters for implementing the PFC controller are chosen after experiments asTr=1ms,nB=2,P=20,S=20. With these, the tracking law of OG around zenith is deduced based on the principle discussed above.


## 4. Simulation and Results

Ji et al. [16] have analyzed the factors influencing the zenith pass problem and concluded that the tracking performance near the zenith is significantly affected by the approaching degree and the angular rate of LOS. In view of this, a set of straight close zenith passes at different degrees and rates are selected to validate the proposed control strategy. The specified reference trajectories of qref,λrefis described as
(40)qref=tan−1tanσcosωt−π/3,λref=sin−1cosσsinωt−π/3,where for the sake of generality, the elevation and azimuth references to be predicted are expressed in a form of a nonlinear function. ω is the sightline angular rate. σ is the zenith pass angle, of which the sine function determines the proximity to the zenith. By performing MATLAB/Simulink simulations of the trajectory tracking on the computer (Intel(R) Core(TM) i5-4200 H, 2.80 GHz, RAM, 8.00 GB), a comparison of the PFC controller with the standard linear PI in position loop is shown from Figures 6–11.Figure 6
OG angle response of PFC and PI (1 deg, 60 deg/s).Figure 7
OG angle response of PFC and PI (1 deg, 60 deg/s).Figure 8
Tracking errors of PFC and PI (5 deg, 30 deg/s).Figure 9
Tracking errors of PFC and PI (1 deg, 30 deg/s).Figure 10
Tracking errors of PFC and PI (5 deg, 60 deg/s).Figure 11
Tracking errors of PFC and PI (1 deg, 60 deg/s).Figures6 and 7 show the OG angle and rate responses, respectively, for an 1 deg zenith pass at 60 deg/s. The curve shown by a broken line is the reference trajectory while the solid line and the dot dash line indicate the tracking performance of PFC and PI, respectively. It can be seen that as t⟶1s, a drastic change occurs in the OG angle reference trajectory. And the angular rate demand at this time tends to be a large value, about 54 rad/s. For tracking this scenario, the preemptive control effect of the PFC control strategy is evident that the PFC controller returns to the accurate track in a short setting time, while the linear PI exhibits a significant hysteresis. Also evident is that the OG angle response error of PFC occurs earlier than that of PI, which seems to indicate PFC has a poor performance. But since the pass is much close to the zenith at this point, the elevation of IG is quite so small that the PFC response error will not lead to a large pointing deviation. The tracking effect is fully demonstrated in Figures 8–11.Figures8–11 show a range of trajectories of tracking errors at different zenith pass angles and rates. In each figure, a 4°×4° square is shown to simulate the field of view of the seeker, and a simple table of lost time to quantify the tracking performance improvement. It is evident that the receding horizon control effectively minimizes the global tracking error for all scenarios. Especially in Figures 9 and 10, the PFC controller kept the target in the field of view of the seeker while the PI controller lost track for 0.24 s, 0.22 s, respectively. Even though the target was impossible to track exactly for both the methods in Figure 11, the PFC controller, compared with PI, still minimized the range of tracking error and provided an extra time of 0.18 s pointing to the target. These indicate PFC a certain promotion of the tracking performance. On the other hand, through comparison among the figures, it can also be seen that the critical pass angle and rate of the track loss is much laxer for PFC than PI: from tracking of (5 deg, 30 deg/s) to losing of (1 deg, 30 deg/s), (5 deg, 60 deg/s) for PI and from tracking of (1 deg, 30 deg/s), (5 deg, 60 deg/s) to losing of (1 deg, 60 deg/s) for PFC.As for the real-time problem of the algorithm, Table3 summarized the single-step time spent in different prediction horizons. It can be easily found that as the prediction horizon expends, the computation time request becomes larger. This is because the increase of predicted points and coincident points makes the dimension of matrices lager and the calculation burden becomes heavier. However, the amount of calculation is so small that even the prediction horizon extends to P=30, the computation time is just 1.76×10−3s, which is far less than 0.02s, the control period of angle tracking.Table 3
Computation time in different prediction horizons.
Prediction horizon
Computation time
P=10
0.98×10−3s
P=20
1.41×10−3s
P=30
1.76×10−3s


## 5. Conclusions

The paper has developed a novel scheme of combining zenith trajectory generation and predictive functional control to deal with the zenith pass problem of the roll-pitch seeker. The generation is achieved over a finite receding horizon under a reasonable assumption, and the PFC controller is implemented in the position loop to minimize the tracking error. Simulation tests around zenith confirm that this scheme is superior to the linear PI control method on all performance indicators and can handle more tracking scenarios of close zenith pass. These suggest the novel control scheme is a practical design technique for zenith controller of the roll-pitch seeker.



---

*Source: 9709341-2020-02-20.xml*