# Nonlinear Predictive Control of Mass Moment Aerospace Vehicles Based on Ant Colony Genetic Algorithm Optimization

Authors: Xiaoyu Zhang; Peng Li; Dexin Xu; Ben Mao; Kunpeng He
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

Based on the mathematical model of the mass moment aerospace vehicles (MMAV), a coupled nonlinear dynamical system is established by rational simplification. The flight control system of MMAV is designed via utilizing nonlinear predictive control (NPC) approach. Aiming at the parameters of NPC is generally used the trial-and-error method to optimize and design, a novel kind of NPC parameters optimization strategy based on ant colony genetic algorithm (ACGA) is proposed in this paper. The method for setting NPC parameters with ACA in which the routes of ants are optimized by the genetic algorithm (GA) is derived. And then, a detailed realized process of this method is also presented. Furthermore, this optimization algorithm of the NPC parameters is applied to the flight control system of MMAV. The simulation results show that the system not only meets the demands of time-response specifications but also has excellent robustness.


---

## Body

## 1. Introduction

The aerodynamic rudder control is generally used for the traditional flight missile in the atmosphere. But in upper air, especially in the exoatmospheric space above 30 kilometers, because of lower air density, the dynamic pressure becomes very low. In the far field of missile fighting, because the speed of missile reduces, the dynamic pressure is relatively low. This causes the dynamic characteristics of the missile to be decreased and the maneuverability to be weak. The weak maneuverability of the missile makes the missile miss the targets frequently and cannot meet the demands of precise guidance and control. This requests the new technical approach. In order to enhance the maneuverability, agility, and guidance precision of the antiaircraft missile, a new method of mass moment control is proposed in 2000 by USA. From the related literatures, the moving-mass-center control was studied by both Russia and USA. In the 1990s of the 20th century. But, because the maneuver control technology of the reentry warhead was kept secret extremely in the world, the public literatures and the related research are very few. The related research which can be seen in the literature has almost focused on reentry warhead at high Mach numbers. Several studies have suggested that mass moment control system (MMCS) appears to offer the greater design and cost advantages [1–8]. There are some advantages of mass moment control as follows. (1) All the mechanism of MMCS is in aerospace vehicles, which will not affect the aerodynamic configuration and is better to accuracy of terminal attack. (2) The actuators of MMCS are internal moving masses, which decrease the thermal load of aerodynamic configuration and avoid the gap on the surface of the vehicle and ablation steering surface. (3) By using aerodynamic forces generated by high-speed flight of vehicle, we can decrease the energy consumption and get effective control avoiding conflict between fuel consumption and control moment generated by lateral jet engine.The stability control mode of the three channels was used for MMAV in this paper. The two orthogonal moving masses were arranged on radial of MMAV, and another moving mass was arranged on axis. Each of the channels has strong coupled aerodynamic, control, inertia, and dynamics. So MMAV becomes one multivariable system model of the nonlinearity and uncertainty. A coupled nonlinear dynamic model is obtained by the rational simplification for this system. The nonlinear predictive control (NPC) of MMAV is investigated in the paper. By considering the insufficiency of NPC adjustment parameters, one parameter optimization method is proposed based on ant colony genetic algorithm (ACGA). With the merits of GA, this method has the global search ability, latent parallel computing, and simple computation process and can solve the contradiction of the optimal solution and the searching space. The contradiction of searching space and time performance is also overcome, which easy to convergent in nonglobal optimal solution and computing time too long when using ant colony algorithm (ACA) to solute optimization question. This method is applied to system design of flight control on MMAV.During the last decade, much effort has been made in GA to investigate aerospace vehicles systems [9–18]. For instance, [9] demonstrates the advantages of using a real coded genetic algorithm (GA) for aerospace engineering design applications. In [10], in order to realize the high precision attitude determination of the space crafts by star sensors, an optimization method of star map distorted model based on improved genetic algorithm is presented. Reference [12] presented a new method in hierarchical genetic algorithms (HGAs) to speed up the optimization of aerodynamic shapes. In [13], an improved GA for airfoil shape parameterization is presented which takes into consideration the characteristics of viscous transonic flow particularly around the trailing edge. In [14], a novel parsimonious genetic programming (PGP) algorithm together with a novel aeroengine optimum data-driven dynamic start process model based on PGP is proposed. In order to obtain the global solution efficiently, [15] proposed a hybrid optimizer combining the advantages of a floating-coded genetic algorithm and simplex method. In [16], an optimization process, based on a genetic algorithm, is used to meet the frequency domain handling qualities requirements for the longitudinal plane. It is well known that GA cannot make good use of system feedback information. It often results in a lot of redundancy and reduces the computing efficiency. In contrast, ant colony algorithm accumulates and renews the optimal path information continuously. It has the ability of parallel processing and global searching. Many scholars have researched the hybrid algorithm including GA and ACA and got some good results in applications. For example, [19] gave an application of ant colony, genetic algorithm, and data mining-based techniques for scheduling. [20] used the hybrid algorithm in the contrast of images. Reference [21] addressed the optimization of cloud database route scheduling based on combination of GA and ACA. More results can be found from [22–28]. However, to the best of our knowledge, there are few works appeared to investigate the MMCS by the GA combined with ACA. This motivates the research in this paper.This paper is organized as follows. The mathematical model of MMAV is derived in Section2. In Section 3, nonlinear predictive control and its robustness analysis are employed to design MMCS. Ant colony genetic algorithm is applied to parameter optimization of MMCS in Section 4. In Section 5, simulation demonstrates the ability of the optimized controller to effectively control the MMAV’s motion. Conclusion is given in Section 6.


## 2.  The Mathematical Model of MMAV

### 2.1. The General Dynamics Model of the MMAV

The basic principle by which MMCS is able to control the vehicle’s motion is to produce the control torque by using the aerodynamic forces and moving the masses within the MMAV to offset the c. m. of system.Suppose that the MMAV includesn moving masses and the mass of MMAV’s shell is mB. The mass of the ith moving mass is mi. So, the total mass of MMAV is mC=mB+∑i=1nmi. The mass ratio of the ith moving mass is μi=mi/mC. The coordinates in the body fixed frame are δ1=(x1m,y1m,z1m)T, i=1,2,…,n. In the ground frame, the velocity of the center of MMAV is VAB=δ˙AB, and the acceleration is V˙AB=δ¨AB. Let the coordinates in the ground with the body fixed frame of the ith moving mass be δdi and the coordinates in the ground frame δAi. There is the relationship δdi=LA1δ1i, δAi=δdi+δAB.The coordinates of MMAV’s c. m. in the ground frame are given by(1)δAC=1mC(mBδAB+∑i=1nmiδAi)=δAB+∑i=1nμiLA1δ1i=δAB+LA1∑i=1nμiδ1i.After derivation, the translational equation of the MMAV in the ground frame can be presented below:(2)V˙AB=RmC+g-((ω1)A×*+(ω1)A×2)×LA1∑i=1nμiδ1i-2(ω1)A×LA1×∑i=1nμiδ˙1i-LA1∑i=1nμiδ¨1i.Then, the translational equation of the MMAV in the body fixed frame is as follows:(3)V˙1B+(ω1)1×V1B=R1mC+g1-((ω1)1×*+(ω1)1×2)∑i=1nμiδ1i-2(ω1)1×∑i=1nμiδ˙1i-∑i=1nμiδ¨1i,
where LA1 is the transformation matrix from the body fixed frame to the ground frame and (ω1)A× is the antisymmetric matrix of the angular velocity of the MMAV in the ground frame. (ω1)1× is the antisymmetric matrix of the angular velocity of the MMAV in the body fixed frame.Correspondingly, the force equation of theith moving mass in the body fixed frame is given by
(4)F1i=μiR1+mi[((ω1)1×*+(ω1)1×2)δ1iF1i=μiR1+mij+2(ω1)1×δ˙1i+δ¨1i]F1i=-mi[∑i=1nμiδ˙1i((ω1)1×*+(ω1)1×2)F1i=μiR1-×∑i=1nμiδ1i+2(ω1)1×F1i=-mi[×∑i=1nμiδ˙1i+∑i=1nμiδ¨1i].According to D’lembert principle, the rotational equation in body coordinates is obtained as follows:(5)(dHdt)1+M1F=dH1dt+(ω1)1×H1+M1F=M1R,
where M1F=-∑i=1nS1iF1i, S1i is the antisymmetric matrix of the ith moving mass in the body fixed frame representing the position coordinates δ1i.


### 2.2. The Dynamics Model of the MMAV with Three Moving Masses

The structure diagram of the MMAV with three moving masses is shown in Figure1. To quickly adjust the flying attitude and decrease coupling, one mass is fixed at x-axis in the body fixed frame. Other two masses are fixed at y-axis and z-axis in radial direction through the MMAV’s axis.Figure 1
. MMAV with three moving masses.The mass of the MMAV’s shell ismB. The mass of axial moving mass A1 is m1, and the coordinate in the body fixed frame is δ11. The mass of radial moving mass A2 is m2, and the coordinate in the body fixed frame is δ12. The mass of radial moving mass A3 is m3, and the coordinate in the body fixed frame is δ13. So, the total mass is mC=mB+m1+m2+m3, and the mass ratios are μ1=m1/mC, μ2=m2/mC, and μ3=m3/mC, respectively.This section derives the equations of motion fully accounted for the dynamic coupling between the four bodies. The moving masses are allowed to translate with respect to the MMAV’s shell but are not allowed to rotate with respect to the MMAV’s shell. Both the MMAV and the moving masses are assumed to be rigid bodies.In the body fixed frame, the interaction between axial moving massA1 and MMAV is F11, the interaction between radial moving mass A2 and MMAV is F12, and the interaction between radial moving mass A3 and MMAV is F13. Equation (4) can be presented below:
(6)F11=m1[V˙1B+(ω1)1×V1B+((ω1)1×*+(ω1)1×2)δ11+2(ω1)1×δ˙11+δ¨11-g1],F12=m1[V˙1B+(ω1)1×V1B+((ω1)1×*+(ω1)1×2)δ12+2(ω1)1×δ˙12+δ¨12-g1],F13=m3[V˙1B+(ω1)1×V1B+((ω1)1×*+(ω1)1×2)δ13+2(ω1)1×δ˙13+δ¨13-g1].The vector translational dynamics of MMAV can be obtained by (3) as follows:
(7)V˙1B+(ω1)1×V1B+(ω1)1×*(μ1δ11+μ2δ12+μ3δ13)=R1mC+g1-[(ω1)1×2(μ1δ11+μ2δ12+μ3δ13)=R1mC+g1-L+2(ω1)1×(μ1δ˙11+μ2δ˙12+μ3δ˙13)=R1mC+g1-L+(μ1δ¨11+μ2δ¨12+μ3δ¨13)].The rotational dynamics of MMAV obtained by (5) are given by
(8)[J1+(1-μ1)m1S11TS11+(1-μ2)m2S12TS12+(1-μ3)m3S13TS13-μ2m1S11S12T-μ3m1S11S13T-μ1m2S12S11T-μ3m2S12S13T-μ1m3S13S11T-μ2m3S13S12T](ω1)1*=M1R-(μ1S11+μ2S12+μ3S13)R1-(ω1)1×J1(ω1)1+[μ1m2S12+μ1m3S13-(1-μ1)m1S11]×((ω1)1×2δ11+2(ω1)1×δ˙11+δ¨11)+[μ2m1S11+μ2m3S13-(1-μ2)m2S12]×((ω1)1×2δ12+2(ω1)1×δ˙12+δ¨12)+[μ3m1S11+μ3m2S12-(1-μ3)m3S13]×((ω1)1×2δ13+2(ω1)1×δ˙13+δ¨13),
where the antisymmetric matrixes of position coordinates of the three moving masses in the body fixed frame are  S11, S12, and S13, respectively.Furthermore, the equations of motion of MMAV system also include some relative movement functions and nonlinear aerodynamic functions. The equations of motion clearly indicate that the MMCS is a complex nonlinear system which has the variable coefficients and large disturbances caused by the accelerations and velocities of masses.


## 2.1. The General Dynamics Model of the MMAV

The basic principle by which MMCS is able to control the vehicle’s motion is to produce the control torque by using the aerodynamic forces and moving the masses within the MMAV to offset the c. m. of system.Suppose that the MMAV includesn moving masses and the mass of MMAV’s shell is mB. The mass of the ith moving mass is mi. So, the total mass of MMAV is mC=mB+∑i=1nmi. The mass ratio of the ith moving mass is μi=mi/mC. The coordinates in the body fixed frame are δ1=(x1m,y1m,z1m)T, i=1,2,…,n. In the ground frame, the velocity of the center of MMAV is VAB=δ˙AB, and the acceleration is V˙AB=δ¨AB. Let the coordinates in the ground with the body fixed frame of the ith moving mass be δdi and the coordinates in the ground frame δAi. There is the relationship δdi=LA1δ1i, δAi=δdi+δAB.The coordinates of MMAV’s c. m. in the ground frame are given by(1)δAC=1mC(mBδAB+∑i=1nmiδAi)=δAB+∑i=1nμiLA1δ1i=δAB+LA1∑i=1nμiδ1i.After derivation, the translational equation of the MMAV in the ground frame can be presented below:(2)V˙AB=RmC+g-((ω1)A×*+(ω1)A×2)×LA1∑i=1nμiδ1i-2(ω1)A×LA1×∑i=1nμiδ˙1i-LA1∑i=1nμiδ¨1i.Then, the translational equation of the MMAV in the body fixed frame is as follows:(3)V˙1B+(ω1)1×V1B=R1mC+g1-((ω1)1×*+(ω1)1×2)∑i=1nμiδ1i-2(ω1)1×∑i=1nμiδ˙1i-∑i=1nμiδ¨1i,
where LA1 is the transformation matrix from the body fixed frame to the ground frame and (ω1)A× is the antisymmetric matrix of the angular velocity of the MMAV in the ground frame. (ω1)1× is the antisymmetric matrix of the angular velocity of the MMAV in the body fixed frame.Correspondingly, the force equation of theith moving mass in the body fixed frame is given by
(4)F1i=μiR1+mi[((ω1)1×*+(ω1)1×2)δ1iF1i=μiR1+mij+2(ω1)1×δ˙1i+δ¨1i]F1i=-mi[∑i=1nμiδ˙1i((ω1)1×*+(ω1)1×2)F1i=μiR1-×∑i=1nμiδ1i+2(ω1)1×F1i=-mi[×∑i=1nμiδ˙1i+∑i=1nμiδ¨1i].According to D’lembert principle, the rotational equation in body coordinates is obtained as follows:(5)(dHdt)1+M1F=dH1dt+(ω1)1×H1+M1F=M1R,
where M1F=-∑i=1nS1iF1i, S1i is the antisymmetric matrix of the ith moving mass in the body fixed frame representing the position coordinates δ1i.


## 2.2. The Dynamics Model of the MMAV with Three Moving Masses

The structure diagram of the MMAV with three moving masses is shown in Figure1. To quickly adjust the flying attitude and decrease coupling, one mass is fixed at x-axis in the body fixed frame. Other two masses are fixed at y-axis and z-axis in radial direction through the MMAV’s axis.Figure 1
. MMAV with three moving masses.The mass of the MMAV’s shell ismB. The mass of axial moving mass A1 is m1, and the coordinate in the body fixed frame is δ11. The mass of radial moving mass A2 is m2, and the coordinate in the body fixed frame is δ12. The mass of radial moving mass A3 is m3, and the coordinate in the body fixed frame is δ13. So, the total mass is mC=mB+m1+m2+m3, and the mass ratios are μ1=m1/mC, μ2=m2/mC, and μ3=m3/mC, respectively.This section derives the equations of motion fully accounted for the dynamic coupling between the four bodies. The moving masses are allowed to translate with respect to the MMAV’s shell but are not allowed to rotate with respect to the MMAV’s shell. Both the MMAV and the moving masses are assumed to be rigid bodies.In the body fixed frame, the interaction between axial moving massA1 and MMAV is F11, the interaction between radial moving mass A2 and MMAV is F12, and the interaction between radial moving mass A3 and MMAV is F13. Equation (4) can be presented below:
(6)F11=m1[V˙1B+(ω1)1×V1B+((ω1)1×*+(ω1)1×2)δ11+2(ω1)1×δ˙11+δ¨11-g1],F12=m1[V˙1B+(ω1)1×V1B+((ω1)1×*+(ω1)1×2)δ12+2(ω1)1×δ˙12+δ¨12-g1],F13=m3[V˙1B+(ω1)1×V1B+((ω1)1×*+(ω1)1×2)δ13+2(ω1)1×δ˙13+δ¨13-g1].The vector translational dynamics of MMAV can be obtained by (3) as follows:
(7)V˙1B+(ω1)1×V1B+(ω1)1×*(μ1δ11+μ2δ12+μ3δ13)=R1mC+g1-[(ω1)1×2(μ1δ11+μ2δ12+μ3δ13)=R1mC+g1-L+2(ω1)1×(μ1δ˙11+μ2δ˙12+μ3δ˙13)=R1mC+g1-L+(μ1δ¨11+μ2δ¨12+μ3δ¨13)].The rotational dynamics of MMAV obtained by (5) are given by
(8)[J1+(1-μ1)m1S11TS11+(1-μ2)m2S12TS12+(1-μ3)m3S13TS13-μ2m1S11S12T-μ3m1S11S13T-μ1m2S12S11T-μ3m2S12S13T-μ1m3S13S11T-μ2m3S13S12T](ω1)1*=M1R-(μ1S11+μ2S12+μ3S13)R1-(ω1)1×J1(ω1)1+[μ1m2S12+μ1m3S13-(1-μ1)m1S11]×((ω1)1×2δ11+2(ω1)1×δ˙11+δ¨11)+[μ2m1S11+μ2m3S13-(1-μ2)m2S12]×((ω1)1×2δ12+2(ω1)1×δ˙12+δ¨12)+[μ3m1S11+μ3m2S12-(1-μ3)m3S13]×((ω1)1×2δ13+2(ω1)1×δ˙13+δ¨13),
where the antisymmetric matrixes of position coordinates of the three moving masses in the body fixed frame are  S11, S12, and S13, respectively.Furthermore, the equations of motion of MMAV system also include some relative movement functions and nonlinear aerodynamic functions. The equations of motion clearly indicate that the MMCS is a complex nonlinear system which has the variable coefficients and large disturbances caused by the accelerations and velocities of masses.


## 3. Nonlinear Predictive Control and Robustness Analysis

### 3.1. Nonlinear Predictive Control

Consider the following MMAV nonlinear system [29]:
(9)x˙=f(x)+g(x)uy=h(x),
where [ωx1ωy1ωz1γnynz]T=[x1x2x3x4x5x6]T is the state variable, u=[δxδyδz]T∈R3 is the control input, and y=[x4x5x6]T is the output to be controlled by the control input u.Reference output signals of the system are(10)yr(t)=[y1r(t)y2r(t)y3r(t)]T.Generally,y1r(t), y2r(t), and y3r(t) are chosen to be step signal.Then the output tracking errors of the system are(11)[e1e2e3]=[y1-y1r(t)y2-y2r(t)y3-y3r(t)].According to the mathematical model forms of MMAV, select functionS(x) is as follows:
(12)S(x)=[e˙1+c1e1+c2xS1e˙2+c3e2+c4xS2e˙3+c5e3+c6xS3],
where x˙S1=e1, x˙S2=e2, and x˙S3=e3.So,S(x) can be as the linear combination of the output tracking error, its differential term, and its integral term. Furthermore the choice from c1 to c6 should meet the following conditions:
(13)λ2+c1λ+c2=(λ+λ1)(λ+λ2)λ2+c3λ+c4=(λ+λ3)(λ+λ4)λ2+c5λ+c6=(λ+λ5)(λ+λ6),
where λi>0, i=1,2,3,4,5,6.Derivation of functionS(x) is
(14)S˙(x)=A(x)+D(x)u.In order to obtain the predictive control law, the quadratic objective function is chosen as follows:(15)J(u)=12[ST(t+h)QS(t+h)+uT(t)Ru(t)].Q is a positive definite symmetric matrix, R is a positive semidefinite symmetric matrix, h is a small positive number, and S(t+h) is the predictive value of S at time t+h. S(t+h) is expanded approximately at t according to the first-order Taylor’s series:
(16)S(t+h)≈S(t)+hS˙(t)=S(t)+h(A+Du).In order to makeJ(u) least, set ∂J/∂u=0 in (15). The control law is obtained below:
(17)u=-h(R+h2DTQD)-1DTQ(S+hA),
when R=0, and D-1 exist, substitute formula (17) into formula (16) and obtain formula (18) below:
(18)S˙(x)+S(x)h=0.Obviously, whent→∞, S(x) tends to zero. When R≠0 and D-1 exist, S(x) also tends to be boundary by simulation. Because S(x) includes integral term of tracking error, the output tracking error tends to zero inevitably.


### 3.2. Robustness Analysis of the System

Consider the uncertain nonlinear system(19)x˙=f^(x)+Δf(x)+[g^(x)+Δg(x)]uy=h(x),
where f^(x) and g^(x) are the parts of the nominal system and Δf(x) and Δg(x) are the uncertainty parts of the system.Suppose thatf^(x) and g^(x) are the known bounded continuous functions and Δf(x) and Δg(x) satisfy the following conditions:
(20)∥Δf(x)∥<∥f^(x)∥∥Δg(x)∥<∥g^(x)∥.Then, each term in (17) can be written as the following form:
(21)A(x)=A^(x)+ΔA(x)D(x)=D^(x)+ΔD(x)S(x)=S^(x)+ΔS(x),
where A^(x), D^(x), and S^(x) correspond to the known nominal system and ΔA(x), ΔD(x), and ΔS(x) correspond to the uncertain parts of system.Firstly, consider the situation isR=0. According to the parameter of nominal system by formula (17), the control law is obtained by
(22)u=-h-1D^-1(S^+hA^).Substitute formula (21) and formula (22) into formula (14):
(23)S˙=-h-1(1+ΔDD^-1)S-h-1(1+ΔDD^-1)ΔS+ΔA=-h-1(1+ΔDD^-1)(S-ΔS)+ΔA.Because of∥ΔDD^-1∥<1 and ∥ΔS∥<∥S∥, S also is boundary which should be known.SetS(x) as input and e as output. The transfer functions are obtained by (12):
(24)G(s)=[ss2+c1s+c2ss2+c3s+c4ss2+c5s+c6].From formula (24), G(s) is strict stability convergence, S(x) is boundary input, and e is boundary output inevitably.  Because S(x) includes the integral term of e,e(∞) should tend to zero ultimately.To sum up, whenR is equal to 0, the system has strong robustness. Similarly, when R≠0, the system still has strong robustness by analysis above.


## 3.1. Nonlinear Predictive Control

Consider the following MMAV nonlinear system [29]:
(9)x˙=f(x)+g(x)uy=h(x),
where [ωx1ωy1ωz1γnynz]T=[x1x2x3x4x5x6]T is the state variable, u=[δxδyδz]T∈R3 is the control input, and y=[x4x5x6]T is the output to be controlled by the control input u.Reference output signals of the system are(10)yr(t)=[y1r(t)y2r(t)y3r(t)]T.Generally,y1r(t), y2r(t), and y3r(t) are chosen to be step signal.Then the output tracking errors of the system are(11)[e1e2e3]=[y1-y1r(t)y2-y2r(t)y3-y3r(t)].According to the mathematical model forms of MMAV, select functionS(x) is as follows:
(12)S(x)=[e˙1+c1e1+c2xS1e˙2+c3e2+c4xS2e˙3+c5e3+c6xS3],
where x˙S1=e1, x˙S2=e2, and x˙S3=e3.So,S(x) can be as the linear combination of the output tracking error, its differential term, and its integral term. Furthermore the choice from c1 to c6 should meet the following conditions:
(13)λ2+c1λ+c2=(λ+λ1)(λ+λ2)λ2+c3λ+c4=(λ+λ3)(λ+λ4)λ2+c5λ+c6=(λ+λ5)(λ+λ6),
where λi>0, i=1,2,3,4,5,6.Derivation of functionS(x) is
(14)S˙(x)=A(x)+D(x)u.In order to obtain the predictive control law, the quadratic objective function is chosen as follows:(15)J(u)=12[ST(t+h)QS(t+h)+uT(t)Ru(t)].Q is a positive definite symmetric matrix, R is a positive semidefinite symmetric matrix, h is a small positive number, and S(t+h) is the predictive value of S at time t+h. S(t+h) is expanded approximately at t according to the first-order Taylor’s series:
(16)S(t+h)≈S(t)+hS˙(t)=S(t)+h(A+Du).In order to makeJ(u) least, set ∂J/∂u=0 in (15). The control law is obtained below:
(17)u=-h(R+h2DTQD)-1DTQ(S+hA),
when R=0, and D-1 exist, substitute formula (17) into formula (16) and obtain formula (18) below:
(18)S˙(x)+S(x)h=0.Obviously, whent→∞, S(x) tends to zero. When R≠0 and D-1 exist, S(x) also tends to be boundary by simulation. Because S(x) includes integral term of tracking error, the output tracking error tends to zero inevitably.


## 3.2. Robustness Analysis of the System

Consider the uncertain nonlinear system(19)x˙=f^(x)+Δf(x)+[g^(x)+Δg(x)]uy=h(x),
where f^(x) and g^(x) are the parts of the nominal system and Δf(x) and Δg(x) are the uncertainty parts of the system.Suppose thatf^(x) and g^(x) are the known bounded continuous functions and Δf(x) and Δg(x) satisfy the following conditions:
(20)∥Δf(x)∥<∥f^(x)∥∥Δg(x)∥<∥g^(x)∥.Then, each term in (17) can be written as the following form:
(21)A(x)=A^(x)+ΔA(x)D(x)=D^(x)+ΔD(x)S(x)=S^(x)+ΔS(x),
where A^(x), D^(x), and S^(x) correspond to the known nominal system and ΔA(x), ΔD(x), and ΔS(x) correspond to the uncertain parts of system.Firstly, consider the situation isR=0. According to the parameter of nominal system by formula (17), the control law is obtained by
(22)u=-h-1D^-1(S^+hA^).Substitute formula (21) and formula (22) into formula (14):
(23)S˙=-h-1(1+ΔDD^-1)S-h-1(1+ΔDD^-1)ΔS+ΔA=-h-1(1+ΔDD^-1)(S-ΔS)+ΔA.Because of∥ΔDD^-1∥<1 and ∥ΔS∥<∥S∥, S also is boundary which should be known.SetS(x) as input and e as output. The transfer functions are obtained by (12):
(24)G(s)=[ss2+c1s+c2ss2+c3s+c4ss2+c5s+c6].From formula (24), G(s) is strict stability convergence, S(x) is boundary input, and e is boundary output inevitably.  Because S(x) includes the integral term of e,e(∞) should tend to zero ultimately.To sum up, whenR is equal to 0, the system has strong robustness. Similarly, when R≠0, the system still has strong robustness by analysis above.


## 4. Ant Colony Genetic Algorithm

### 4.1. Optimization Principle of Ant Colony Algorithm

When ants in the nature seek foods or meet obstacles, they always can find one optimal route from their nests to food source or round the obstacles. The reason is that the ants can release one special secretion called pheromone on the roads which they pass through (pheromone, this material can volatilize gradually along with time); then the following ants may choose next route according to the pheromone which the fore ants leaved. The pheromone in one route is higher, and the number of times which should be selected is more. That is, the performance of this route will be superior, and the probability that the following ants choose this route will be larger. So, a positive feedback process of learning information is established by this way. At last, the optimal solution can be obtained [30–32]. The ant colony algorithm is based on this principle.For formula (13), when control system is designed by nonlinear predictive control, the parameters to be optimized are ci>0, i=1,2,3,4,5,6. In the generality, set c1=c3=c5 and c2=c4=c6. So, only c1 and c2 need to be optimized.


### 4.2. Generation of Nodes and Routes

By the inspiration of [33], the two parameters to be optimized are expressed in a plane for using the ant colony algorithm expediently. The graph is shown in Figure 2.Figure 2
Schematic plot on the generation of nodes and routes.L
1~L5 and L6~L9 express c1 and c2, respectively. According to the test results in [7], the initial value c1 is 200.00 and the initial value c2 is 10.00. c1 and c2 are calculated, respectively, by the formulas below:
(25)c1=y1,j×102+y2,j×101+y3,j×100+y4,j×10-1+y5,j×10-2(26)c2=y6,j×101+y7,j×100+y8,j×10-1+y9,j×10-2.Take the fold line in Figure2 as example. The expressed parameters are c1=218.54 and c2=15.62.


### 4.3. Establishment of the Objective Function

When the objective function is established, the dynamic performances of the system on establishing overload instructions should be made as basis according to system requests for design, which mainly consider system overshoot, the rise time, and the settling time. The objective function is defined as follows:(27)F=λδ(δδ0)+λtr(trtr0)+λts(tsts0),
where δ0, tr0, and ts0 are the overshoot, the rise time, and the settling time which guidance control system requested. λδ, λtr, and λts are the weight coefficients, respectively. Their values are 0.4, 0.2, and 0.4.


### 4.4. The Choice of the Route Nodes

Suppose that it takes each ant the equal time to crawl from any node on lineLi to the next any node on line Li+1 and has nothing to do with the distance between the nodes. Therefore, if all ants start from the origin of coordinates, they will simultaneously arrive at each line Li (i = 1~9). Finally, they will also simultaneously arrive at their end points on the line L9. Then, one circulation is completed. Define τ(xi,yi,j,t) to express the remaining pheromone on node (xi,yi,j) at t. The pheromone is equal on each node in the initial time. In order to make the algorithm search more feasible solutions in the preliminary stage, τ(xi,yi,j)=τmax and Δτ(xi,yi,j,0)=0 are set at initialization. Suppose that Pk(xi,yi,j,t) expresses the probability in which the ant k crawls from any node on Li-1 to node (xi,yi,j) at t; then
(28)Pk(xi,yi,j,t)=τα(xi,yi,j,t)ηβ(xi,yi,j,t)∑j=09τα(xi,yi,j,t)ηβ(xi,yi,j,t),in whichη(xi,yi,j,t) is the visibility on node (xi,yi,j),
(29)(xi,yi,j,t)=10-|yi,j-yi,j*|10,
where the value yi,j* is obtained as the following method. In the first circulation of ant colony algorithm, yi,j* is given by the y-coordinate value of each node which c1 and c2 map in Figure 2. In each circulation later, yi,j* is given by the y-coordinate value of each node in which c1 and c2 correspond to the optimal route produced in the previous circulation (the optimal performances in previous circulation) map in Figure 2.


### 4.5. Renewal and Limitation of the Pheromone

Suppose thatt=0 at initial time; all ants are located at the origin of coordinates O. After 9 units of time, all ants crawl from the initial point to the end point. The pheromone in the route may be adjusted by the following formula:
(30)τ(xi,yi,j,t+9)=ρ(xi,yi,j,t)+Δτ(xi,yi,j).For using more optimal solution information during the running period of algorithm, the renewal and limitary mechanism of the pheromone can use the maximal minimal ant system (MMAS), in which only the one optimal ant can increase pheromone after each iterative. And there isτ(xi,yi,j)∈[τmin,τmax]:
(31)Δτ(xi,yi,j)best=qFbest(t).Fbest(t) represents the optimal objective function value in this circulation. It can be calculated by (27).


### 4.6. The Steps of Optimizing MMCS Parameter Using ACGA

The hybridization of AC and GA has been studied in many works. The methods have been successfully applied to diverse combinatorial optimization problems including traveling salesman, quadratic assignment, vehicle routing, telecommunication networks, image processing, constraint satisfaction, and scheduling. Recently, researchers have been dealing with the relation of ACGA to the fields of optimal control and reinforcement learning, in which, with the merits of GA, this method has the global search ability, latent parallel computing, and simple computation process and can solve the contradiction of the optimal solution and the searching space. It also can overcome the contradiction of searching space and time performance, in contrast with easing to converge to nonglobal optimal solution and computing time too long when using AC to solute optimization question. So, the control parameters optimization method of MMAV is proposed based on ACGA. The goals are to meet the demands of time-response specifications and improve the robustness for MMCS. Several studies have suggested that MMCS appears to offer the greater design and cost advantages in flight vehicles control. We present the mathematical modeling process of MMAV and describe the method of nonlinear predictive control for MMCS. But it is very difficult to design the MMCS because there are many parameters to adjust manually. In this paper, a new method is proposed using ACGA, by which automatic parameter optimization can be realized instead of manually adjusting.When the MMCS parameters are optimized by the ant colony algorithm, the variation operation of GA can be used for the second optimization of the obtained result [19, 34]. That is, to avoid to fall into partial optimal by using ant colony algorithm to obtain a group partial solution that makes the objective function expressed by formula (27) is minimal. Then, the coding is carried on (the binary system is used in this paper), and the variation of certain digits is caused stochastically. If the variation solution surpasses the original solution in genetic and variation process, then the variation solution is retained. Otherwise, the original solution is used. After this variation, seeking in the partial minimal region can be jumped out, and the quality of solutions can be enhanced; the convergence velocity of the optimal solution is speeded up. So, ACGA can realize the highly effective, fast, and complete seeking to the search space.The basic steps to optimize MMCS parameters with ACGA are described as follows:Step 1.
Initialize the MMCS parameters and the correlation ACGA parameters. The initial values ofc1 and c2 can be seen above.Step 2.
Set ant number asm. Define the one-dimensional array Routek including 9 elements for each ant. The y-coordinate values of 9 nodes which the ant will pass through are saved in this array, that is, the route of the ant crawling.Step 3.
Make time countert=0 and cycle number Nc=0. Define the biggest cycle number and the pheromone values on each node at initial time; set Δτ(xi,yi,j)=0. Put all the ants at the initial station O.Step 4.
Seti=1.Step 5.
Compute the transfer probabilities that the ant moves to each node on the lineLi by formula (28). According to these probabilities, a node for each ant on the line Li is chosen by roulette principle, and the ant also will be moved to this node. Simultaneously, the y-coordinate value of this node is saved in ith element of Routek.Step 6.
Seti=i+1. If there is i≤9, jump to Step 4. Otherwise, jump to Step 7.Step 7.
c
1
k and c2k corresponding to this route are calculated by formula (25) and formula (26) according to Routek, that is, the routes passed by the ants. The dynamic performances δk,trk, and tsk can be calculated by carrying on computer simulation applying these parameters to the MMCS system of MMAV. Then, the objective function Fk corresponding to the ant k is calculated by formula (27), and Fbest(t) is obtained. The optimal solution is saved in this circulation. The MMCS parameters corresponding to these also are saved into c1* and c2*. Then, the second optimization is carried on by GA. The optimal solution from the original solution and the mutation solution is chosen.Step 8.
Sett←t+9 and Nc←Nc+1. According to formula (30) and formula (31), the pheromones on each node are renewed, and Routek is reset to zero.Step 9.
IfNc<Ncmax and the whole ant colony is not convergent to one route, all ants will be put at the initial station O again, and then jump to Step 4. If Nc<Ncmax, but the whole ant colony converges to one same route, the algorithm will be ended and the corresponding optimized MMCS parameters c1* and c2* are achieved.


## 4.1. Optimization Principle of Ant Colony Algorithm

When ants in the nature seek foods or meet obstacles, they always can find one optimal route from their nests to food source or round the obstacles. The reason is that the ants can release one special secretion called pheromone on the roads which they pass through (pheromone, this material can volatilize gradually along with time); then the following ants may choose next route according to the pheromone which the fore ants leaved. The pheromone in one route is higher, and the number of times which should be selected is more. That is, the performance of this route will be superior, and the probability that the following ants choose this route will be larger. So, a positive feedback process of learning information is established by this way. At last, the optimal solution can be obtained [30–32]. The ant colony algorithm is based on this principle.For formula (13), when control system is designed by nonlinear predictive control, the parameters to be optimized are ci>0, i=1,2,3,4,5,6. In the generality, set c1=c3=c5 and c2=c4=c6. So, only c1 and c2 need to be optimized.


## 4.2. Generation of Nodes and Routes

By the inspiration of [33], the two parameters to be optimized are expressed in a plane for using the ant colony algorithm expediently. The graph is shown in Figure 2.Figure 2
Schematic plot on the generation of nodes and routes.L
1~L5 and L6~L9 express c1 and c2, respectively. According to the test results in [7], the initial value c1 is 200.00 and the initial value c2 is 10.00. c1 and c2 are calculated, respectively, by the formulas below:
(25)c1=y1,j×102+y2,j×101+y3,j×100+y4,j×10-1+y5,j×10-2(26)c2=y6,j×101+y7,j×100+y8,j×10-1+y9,j×10-2.Take the fold line in Figure2 as example. The expressed parameters are c1=218.54 and c2=15.62.


## 4.3. Establishment of the Objective Function

When the objective function is established, the dynamic performances of the system on establishing overload instructions should be made as basis according to system requests for design, which mainly consider system overshoot, the rise time, and the settling time. The objective function is defined as follows:(27)F=λδ(δδ0)+λtr(trtr0)+λts(tsts0),
where δ0, tr0, and ts0 are the overshoot, the rise time, and the settling time which guidance control system requested. λδ, λtr, and λts are the weight coefficients, respectively. Their values are 0.4, 0.2, and 0.4.


## 4.4. The Choice of the Route Nodes

Suppose that it takes each ant the equal time to crawl from any node on lineLi to the next any node on line Li+1 and has nothing to do with the distance between the nodes. Therefore, if all ants start from the origin of coordinates, they will simultaneously arrive at each line Li (i = 1~9). Finally, they will also simultaneously arrive at their end points on the line L9. Then, one circulation is completed. Define τ(xi,yi,j,t) to express the remaining pheromone on node (xi,yi,j) at t. The pheromone is equal on each node in the initial time. In order to make the algorithm search more feasible solutions in the preliminary stage, τ(xi,yi,j)=τmax and Δτ(xi,yi,j,0)=0 are set at initialization. Suppose that Pk(xi,yi,j,t) expresses the probability in which the ant k crawls from any node on Li-1 to node (xi,yi,j) at t; then
(28)Pk(xi,yi,j,t)=τα(xi,yi,j,t)ηβ(xi,yi,j,t)∑j=09τα(xi,yi,j,t)ηβ(xi,yi,j,t),in whichη(xi,yi,j,t) is the visibility on node (xi,yi,j),
(29)(xi,yi,j,t)=10-|yi,j-yi,j*|10,
where the value yi,j* is obtained as the following method. In the first circulation of ant colony algorithm, yi,j* is given by the y-coordinate value of each node which c1 and c2 map in Figure 2. In each circulation later, yi,j* is given by the y-coordinate value of each node in which c1 and c2 correspond to the optimal route produced in the previous circulation (the optimal performances in previous circulation) map in Figure 2.


## 4.5. Renewal and Limitation of the Pheromone

Suppose thatt=0 at initial time; all ants are located at the origin of coordinates O. After 9 units of time, all ants crawl from the initial point to the end point. The pheromone in the route may be adjusted by the following formula:
(30)τ(xi,yi,j,t+9)=ρ(xi,yi,j,t)+Δτ(xi,yi,j).For using more optimal solution information during the running period of algorithm, the renewal and limitary mechanism of the pheromone can use the maximal minimal ant system (MMAS), in which only the one optimal ant can increase pheromone after each iterative. And there isτ(xi,yi,j)∈[τmin,τmax]:
(31)Δτ(xi,yi,j)best=qFbest(t).Fbest(t) represents the optimal objective function value in this circulation. It can be calculated by (27).


## 4.6. The Steps of Optimizing MMCS Parameter Using ACGA

The hybridization of AC and GA has been studied in many works. The methods have been successfully applied to diverse combinatorial optimization problems including traveling salesman, quadratic assignment, vehicle routing, telecommunication networks, image processing, constraint satisfaction, and scheduling. Recently, researchers have been dealing with the relation of ACGA to the fields of optimal control and reinforcement learning, in which, with the merits of GA, this method has the global search ability, latent parallel computing, and simple computation process and can solve the contradiction of the optimal solution and the searching space. It also can overcome the contradiction of searching space and time performance, in contrast with easing to converge to nonglobal optimal solution and computing time too long when using AC to solute optimization question. So, the control parameters optimization method of MMAV is proposed based on ACGA. The goals are to meet the demands of time-response specifications and improve the robustness for MMCS. Several studies have suggested that MMCS appears to offer the greater design and cost advantages in flight vehicles control. We present the mathematical modeling process of MMAV and describe the method of nonlinear predictive control for MMCS. But it is very difficult to design the MMCS because there are many parameters to adjust manually. In this paper, a new method is proposed using ACGA, by which automatic parameter optimization can be realized instead of manually adjusting.When the MMCS parameters are optimized by the ant colony algorithm, the variation operation of GA can be used for the second optimization of the obtained result [19, 34]. That is, to avoid to fall into partial optimal by using ant colony algorithm to obtain a group partial solution that makes the objective function expressed by formula (27) is minimal. Then, the coding is carried on (the binary system is used in this paper), and the variation of certain digits is caused stochastically. If the variation solution surpasses the original solution in genetic and variation process, then the variation solution is retained. Otherwise, the original solution is used. After this variation, seeking in the partial minimal region can be jumped out, and the quality of solutions can be enhanced; the convergence velocity of the optimal solution is speeded up. So, ACGA can realize the highly effective, fast, and complete seeking to the search space.The basic steps to optimize MMCS parameters with ACGA are described as follows:Step 1.
Initialize the MMCS parameters and the correlation ACGA parameters. The initial values ofc1 and c2 can be seen above.Step 2.
Set ant number asm. Define the one-dimensional array Routek including 9 elements for each ant. The y-coordinate values of 9 nodes which the ant will pass through are saved in this array, that is, the route of the ant crawling.Step 3.
Make time countert=0 and cycle number Nc=0. Define the biggest cycle number and the pheromone values on each node at initial time; set Δτ(xi,yi,j)=0. Put all the ants at the initial station O.Step 4.
Seti=1.Step 5.
Compute the transfer probabilities that the ant moves to each node on the lineLi by formula (28). According to these probabilities, a node for each ant on the line Li is chosen by roulette principle, and the ant also will be moved to this node. Simultaneously, the y-coordinate value of this node is saved in ith element of Routek.Step 6.
Seti=i+1. If there is i≤9, jump to Step 4. Otherwise, jump to Step 7.Step 7.
c
1
k and c2k corresponding to this route are calculated by formula (25) and formula (26) according to Routek, that is, the routes passed by the ants. The dynamic performances δk,trk, and tsk can be calculated by carrying on computer simulation applying these parameters to the MMCS system of MMAV. Then, the objective function Fk corresponding to the ant k is calculated by formula (27), and Fbest(t) is obtained. The optimal solution is saved in this circulation. The MMCS parameters corresponding to these also are saved into c1* and c2*. Then, the second optimization is carried on by GA. The optimal solution from the original solution and the mutation solution is chosen.Step 8.
Sett←t+9 and Nc←Nc+1. According to formula (30) and formula (31), the pheromones on each node are renewed, and Routek is reset to zero.Step 9.
IfNc<Ncmax and the whole ant colony is not convergent to one route, all ants will be put at the initial station O again, and then jump to Step 4. If Nc<Ncmax, but the whole ant colony converges to one same route, the algorithm will be ended and the corresponding optimized MMCS parameters c1* and c2* are achieved.


## 5. Simulation Result and Discussion

Because the control performances of MMAV need to be realized through an axial moving mass and two-radial-moving-mass to and fro movement on the guide rail, the control effects of the system directly are affected by the characteristics of servo system. So, in the design of MMCS, the dynamic characteristics of servo system have to be considered.Considering the band width of the measure elements for system (accelerometer and rate gyro) is relatively large, the simplification model of the servo system is obtained by repeated test:(32)GR(s)=1TR2s2+2ζRTRs+1,
where TR=0.01 and ζR=0.8 to the axial moving mass, TR=0.005 and ζR=0.8 to the radial moving mass. Simultaneously, the saturated constraints of the displacements should be considered for the whole design assigns when the system is tested.According to the ant colony genetic algorithm proposed in this paper to carry on parameter optimization, the final parameters of MMAV in nonlinear predictive control system are shown in Figure2. Simultaneously, 6DOF mathematical simulation of MMAV was carried on. The signals y˙ir and y¨ir(i=1,2,3) in the flight control system are obtained by the second-order filter. The simulation results in which the aerodynamic parameters have 50% uncertainty are shown at length of article limited and the system design results in nominal state as the foundation. The simulation graphs are shown in Figures 3 and 4, which are state responses of MMAV using NPC with the original GA and ACGA, respectively. The time response and system stability of ACGA are better than those of the original GA. The simulation curves of the normal overload and each control of MMAV are only shown in this paper (because the radial moving mass is symmetry, the displacement curves of oy1 axis and oz1 axis for the moving mass are the same). (Only the displacement curve of oy1 axis for the moving mass is given). From the graphs, the displacement of the axial moving mass is limited in the scope ±0.16m, and the displacements of the two radial moving masses are limited in the scope ±0.08m. That flight control of MMAV may be realized by only very short displacement of the moving mass. The dynamic performances of control system for MMAV still satisfy the design requests although aerodynamic parameters have been changed 50% with their nominal values.Figure 3
State response of MMAV using NPC based on GA optimization.Figure 4
State response of MMAV using NPC based on ACGA optimization.As shown in simulation results, the nonlinear predictive control is proposed to design MMCS, in which controller parameters are optimized with the original GA and ACGA. The simulation results show that MMCS with ACGA has better dynamic behavior and strong robustness than that with the original GA. This approach can improve the robustness of MMCS and make it work effectively. The theoretical analysis and the simulation results of the flight control system for MMAV show that the strategy of the NPC parameter optimization based on the ant colony genetic algorithm proposed in this paper is effective and feasible. The algorithm of the flight control system for MMAV designed by using NPC based on ACGA optimization is simple. Moreover, the dynamic performances of the system satisfy the design requests completely and also have strong robustness. Simultaneously, the time is saved. The optimization algorithm proposed in this paper does not rely on the precise mathematical model of MMAV and makes processing question more flexible, adaptable, and robust. This method can not only improve the quality of control system design but also reduce the design difficulty and may solve optimization problems of the control system parameters which are designed by other control methods. So, it has important reference value in the control engineering domain.


## 6. Conclusion

An autopilot for a nonlinear six-degree-of-freedom MMAV is introduced in this paper. Because the stability control system used by MMAV with three moving masses is three-channel, and the deduced mathematical model is too complicated, the dynamical model of MMAV is still a nonlinear system after reasonable simplification. Therefore, this control system is hard to design. The nonlinear predictive control method is proposed to design MMCS, in which controller parameters are optimized based on ACGA. This approach can improve the robustness of MMCS and make it work effectively. The simulation results show that MMCS has good dynamic behavior and strong robustness. As for the mechanism of mass moment control, it is very complicated and the related research is only in the beginning. The public references are insufficient; we only give some useful discussion about the control system design method of MMAV in this paper.



---

*Source: 858190-2013-10-22.xml*