# Neural Network Control for the Probe Landing Based on Proportional Integral Observer

Authors: Yuanchun Li; Tianhao Ma; Bo Zhao
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

For the probe descending and landing safely, a neural network control method based on proportional integral observer (PIO) is proposed. First, the dynamics equation of the probe under the landing site coordinate system is deduced and the nominal trajectory meeting the constraints in advance on three axes is preplanned. Then the PIO designed by using LMI technique is employed in the control law to compensate the effect of the disturbance. At last, the neural network control algorithm is used to guarantee the double zero control of the probe and ensure the probe can land safely. An illustrative design example is employed to demonstrate the effectiveness of the proposed control approach.


---

## Body

## 1. Introduction

The exploration mission to near-earth asteroids (NEAs) would be one of the most complex tasks in the future deep space exploration [1, 2]. There is surge in NEAs mission activities, for which various space agencies around the world (e.g., NASA, European Space Agency, Japan Aerospace Exploration Agency, etc.) were commissioning researches about NEAs to determine the feasible exploration missions, including the (1) NEAR probe launched by NASA which can realize the fly-around to 433 Eros whose shape is like a potato with a size of 34.4 km × 11.2 km × 11.2 km and which verified the gravitational field model of 433 Eros and the stability of frozen orbit around the asteroid [3]; (2) the Hayabusa probe from JAEA which had successfully achieved to be attached and sample to the 25143 Itokawa (due to the smaller size and quality of the Itokawa asteroid, this mission realized the detection to the asteroid by hovering way [4]); (3) ROSETTA implemented by ESA which will arrive in the Churyumov-Gerasimenko comet in 2014 after a decade of interstellar flight and will make the comprehensive observation of the comet for a long time [5].In view of the complex environment around the small body, together with the long distance between the probe and the surface of the earth [6], a variety of accurate physical parameters and motion information of small bodies cannot be obtained through optical telescopes on the ground or radio telescopes. In addition, the complex process uncertainty, large time delay, nonlinearity, and multivariable coupling always exist in the probe dynamic model, so ground control for deep space exploration mission has become no longer appropriate; as a consequence, it puts forward a new challenge to autonomous navigation, guidance, and control (GNC) technology of landing softly on a small body. To cope with these problems, both domestic and foreign scholars have paid a great deal of attention to the GNC problem of landing small objects. As is well known, the accurate physical parameters and motion information of small bodies are the important premises of the probe softly landing. Misu et al. [6] proposed an autonomous optical navigation and guidance method, which extracted visual small features from the images taken by the navigation camera and tracked them robustly and accurately. Kawaguchi et al. [7] discussed an autonomous optical guidance and navigation strategy to approaching small bodies. Horneman and Kluever [8] presented a terminal area energy management (TAEM) guidance methodology which employed a trajectory planning algorithm to compute a feasible path from the current state to the desired approach and landing target state rather than relying on a precalculated one, stored database of neighboring TAEM trajectories. However, even if the accurate physical parameters and motion information of small bodies are gained, the controller is difficult to be designed to make the probe system meet the key performance indicators of the probe softly landing. In order to solve this problem, Furfaro et al. [9] presented a high order sliding mode variable structure control method to make the probe reach the sliding surface in finite time and overcame the chattering effect, generally existing in the common sliding mode control. Crassidis et al. [10] introduced a variable-structure controller based on a Gibbs vector parameterization, a modified-Rodrigues parameterization, and a quaternion parameterization. Blackmore [11] studied the robust path and feedback control under the condition of existing uncertainty; through this control method, the stability of the system is ensured. Meissinger and Greenstadt [12] proposed a soft landing scheme, which used a feedback control with a radar altimeter and a three-beam Doppler radar system to achieve landing spacecraft at Eros’ north polar region with a low-impact velocity. In [13], a novel robust stability condition was obtained for sliding mode dynamics by using Lyapunov theory in delta domain. Some other approaches for analysis and design of sliding mode control were presented in [14–16]. Apart from the position and the velocity of the probe, the attitude dynamics analyses also play an important role in the probe softly landing. Kumar and Shah [17] set up the general formulation of the spacecraft equations of motion in an equatorial eccentric orbit using Lagrangian method and did some analysis on the stability. Then the control laws for three-axis attitude control of spacecrafts had been developed and a closed-form solution of the system had been derived. Liang and Li [18] designed a robust adaptive backstepping sliding mode control law to make the attitude of the probe stabilized and respond accurately to the expectation in the presence of disturbances and parametric uncertainties. Nonetheless, these methods of dealing with the interference all made the inhibition of bounded disturbances implicit in the above autonomous GNC rather than using the interference information effectively, so the designed controller cannot meet the control requirement of the system when there exists a larger interference in the system.As a result of the complex environment in the deep space around the small bodies and the coupling effect of the detector itself, it leads to a great deal of uncertainties in the dynamic model and makes the system include the complex external disturbance. At present, the main approaches to process the external disturbances include disturbance decoupling, disturbance compensation, and robust control, especially disturbance compensation. There are many scholars proposing a variety of stability control strategies based on the observer aimed at different objects. Chadli and Karimi [19] dealt with the observer design for Takagi-Sugeno (T-S) fuzzy models subject to unknown inputs and disturbance affecting both states and outputs of the system. Chong et al. [20] designed a robust circle criterion observer and applied it to neural mass models; Sun et al. [21] proposed a novel speed observation scheme using artificial neural network (ANN) inverse method to effectively reject the influence of speed detection on system stability and precision for a bearingless induction motor. Above all, the observer can inhibit the effect of disturbance in the system by accurately measuring the unknown disturbance.The main advantages of the presented approach are generalized into two aspects: one is that it combines the characteristics of the probe dynamic model and the good estimation performance of observer, to eliminate the effect of the unknown disturbance and to avoid the chattering of the control signal caused by the large disturbance. This paper designs PIO using LMI technique, which can estimate the system states and unknown input disturbance simultaneously. The other is that PID neural network control algorithm is introduced in the design of the controller. It combines the advantages of traditional PID controller and learning memory function of neural networks. So it improves the convergence rate close to the ideal position on the condition that the convergence of the system can be ensured and simultaneously avoids the effect of nonlinear and strong coupling features of the system in a wide range, compared with the sliding mode control strategy.This paper proceeds as follows. In Section2, the dynamics equation of the probe is deduced under the landing site coordinate system and the interference outside system is treated as the known bounded function. In Section 3, firstly, the nominal trajectories based on the theory of suboptimal fuel are planned. Then PIO is designed by using LMI technique to estimate the unknown disturbance. Finally, PID neural network control algorithm is used to design the controller to ensure the stability and control performance of the system. In Section 4, Eros 433 is employed to demonstrate the effectiveness of the proposed control approach. Conclusions are presented in Section 5.


## 2. Small Body and Probe Dynamic Model

In this section, the body-fixed coordinate system of small body is set up, which is shown in Figure1. Let the oa-xayaza coordinate system be fixed on small body with the origin coinciding with the mass center of small body, xa-axis coinciding with the minimum inertia axis of small body, za-axis coinciding with the spin axis of small body, and ya-axis meeting the condition that the xa, ya, and za axes compose the right-handed coordinate system. The oc-xcyczc coordinate system is fixed on optical navigation camera (ONC), and the image plane of ONC is defined as ocxcyc, and zc axis is parallel to the optical axis of ONC and is directed to the surface.Figure 1
Geometrical relationship of coordinate systems.The dynamic equations of the probe in the fixed-body coordinate system are given as [22]
(1)R¨+2ω×R˙+ω×ω×R+ω˙×R=a+UR+fd,
where R, R˙, R¨, ω, a, UR, and fd are the position vector from the target small body mass center of the spacecraft, the first and second time derivatives with respect to the body-fixed rotating frame, the instantaneous rotation vector of the small body, the control acceleration, the gradient of the gravitational potential U, and the components of unmodeled perturbation accelerations mainly from the solar radiation pressure and the solar gravitation.Considering the origin ofΣl is the vector ρ in the Σa, which is the vector from the target small body mass center to the landing site, the vector R has the following satisfaction relation in the body-fixed coordinate system Σa:
(2)R=Clar+ρ,
where r and Cla are the vector from the landing site to the probe in the Σl and coordinate transform matrix from the Σl to the Σa. The transform matrix is given as follows:
(3)Cla=cos⁡φsinθ-sinφcos⁡φcos⁡θsinφsinθcos⁡φsinφcos⁡θ-cos⁡θ0sinθ.Suppose that the small body rotates around thez-axis and rotation velocity ω is a constant; we can get the final expression of dynamic models as
(4)x¨=ω2x+2ωy˙+Ux+acx+fdx,y¨=ω2y-2ωx˙+Uy+acy+fdy,z¨=Uz+acz+fdz,
where fdx, fdy, and fdz are the components of unmodelled perturbation accelerations mainly from the solar radiation pressure and the solar gravitation.Generally, given the small size, irregular shape, and variable surface properties of small bodies, orbital dynamics became complicated; thus it is difficult to obtain the gravitational field of the small bodies accurately. Considering that the gravitational potential is related to the distance, the latitude, and the longitude, it can be expanded into a series of spherical harmonics and can be expressed as(5)U=μr1+Rar212C203sin2φ-1Rar2Rar212C203sin2φ-1Rar2+3C22cos⁡2φcos⁡2θ,
where μ, Ra, φ, θ, and r are the product of the gravitational constant and the mass of the target small body, the referenced radius which is similar to the large equatorial radius, the latitude and longitude in the same coordinate system whose origins are at the center of body mass, and the distance from the mass center of small body to the probe, respectively.According to the relationship between the rectangular coordinate and polar coordinate, one obtains(6)sinφ=zr,cos⁡2φ=x2+y2r2,cos⁡2θ=1-2sin2θ=x2-y2x2+y2.Introduce (6) into (5), and one can obtain
(7)U=μr1+Rar212C202z2-x2-y2r2x4-y4r2(x2+y2)Rar212C202z2-x2-y2r2+3C22x4-y4r2x2+y2.Furthermore, the derivatives ofU can be computed explicitly with respect to x, y, and z, respectively, as
(8)Ux=-μxr31+32C20Rar25z2r2-1+3C22Rar25x2-y2r2-2,Uy=-μyr31+32C20Rar25z2r2-1+3C22Rar25x2-y2r2-2,Uz=-μzr31+32C20Rar25z2r2-1+3C22Rar25x2-y2r2-2.


## 3. Guidance Law and Control Law Design

Considering the probe achieves the vertical soft landing within the expected timeτ, this paper presents the nominal trajectory guidance law based on the theory of suboptimal fuel, and the nominal trajectories of three-axis direction are preplanned. Then use neural network control method based on PIO to track the planned ideal nominal trajectories.


### 3.1. The Nominal Trajectory Planning

The desired descent altitude and velocity are planned in order to satisfy the requirements of soft landing on the surface of small bodies. The constraint condition is defined as [23]
(9)z˙0=z˙0,z(0)=z0,zτ=ρ,z˙τ=0,
where z0 and z˙0 denote the initial altitude and altitude change rate, z˙0 and z0 are the planned altitude and altitude change rate, and τ is the descent time. The cubic curve to satisfy the boundary condition is given by
(10)znt=z0+z1t+z2t2+z3t3,
where z0, z1, z2, and z3 are the cubic function coefficients.Using (9), the coefficients are determined. The descent curve is given by
(11)znt=z0+z˙0t-3z0+2z˙0τ-3ρτ2t2+2z0+z˙0τ-2ρτ3t3,
where ρ is the altitude of the landing site.Next, the time derivatives of (12) are given by
(12)z˙nt=z˙0-6z0+4z˙0τ-6ρτ2t+6z0+3z˙0τ-6ρτ3t2,(13)z¨nt=-6z0+4z˙0τ-6ρτ2+12z0+6z˙0τ-12ρτ3t.Similarly, the ideal nominal trajectories can be planned on the other two axis directions.


### 3.2. Control Law Design

#### 3.2.1. Proportional Integral Observer Design

Let(14)xt=x,y,z,x˙,y˙,z˙T,ut=000uxt-Uxuyt-Uyuzt-UzT.Convert (4) into nonlinear time-invariant system as
(15)x˙t=Axt+But+ft,yt=Cxt,
where
(16)A=000100000010000001ω20002ω00ω20-2ω00000000,B=000b1b2b3,ft=000fdxfdyfdz,C=111111,
where ω is the instantaneous rotation vector of the small body and fdx, fdy, and fdz are the components of unmodelled perturbation accelerations mainly from the solar radiation pressure and the solar gravitation.Next, the PIO is designed as follows [24]:
(17)x^˙t=Ax^t+But+f^t+Kpyt-y^t,y^t=Cx^t,f^˙t=KIyt-y^t,
where Kp and KI are the observer gain matrix and the integral coefficient of estimated unknown disturbance, respectively.Note the state error and unknown disturbance error as(18)et=x^t-xt,eft=f^t-ft,
where x^ and f^ are the estimations of the state vector x and the unknown disturbance f, respectively.Using (17) and (18) the error dynamics are as follows:
(19)e˙t=A-KPCet+ft-f^t,e˙f=-KIC+KvCA-KpCet-KvCef-f˙t.The augmented estimator system could be rewritten as(20)e˙te˙ft=Aefeef+Beff˙t,
where
(21)Aef=A-KPC1-KIC-KvCA+KvCKpC0,Bef=01.Lemma 1 (see [25]).
Givenγ>0 and (20), if there exist symmetric matrixes P, Q and two matrices KP, KI of appropriate dimension as well as LMI such that
(22)A11A120A12TA22-QS0-QSTγI1<0,
where
(23)A11=PA-KpC+A-KpCTP+Ie,A12=P+QSKiC+KvA-KpCT,A22=-QSKvC-KiFS+KvCQST+Is
hold, then system (20) is stable and satisfies corresponding performance index.Theorem 2.
For a given positive constantγ>0 and (20), if there exist symmetric matrixes P, Q and two matrices KP, KI of appropriate dimension such that (22) LMI holds, then system (20) is robust and asymptotically stable and satisfies the performance index as follows:
(24)eef≤γf^2+V0,
where
(25)eef2=∫0t1eefeefTdt,f^2=∫0t1f^f^Tdt.Ie, Is, and I1 are unit matrixes with appropriate dimension.Proof.
Choose the Lyapunov candidate as(26)V=eefP1eefT;
then taking the derivative ofV with respect to time along the trajectories of (19), one obtains
(27)V˙=eefAefTP1+P1AefeefT+2eefTP1Beff^.
Now define performance indicators as follows:(28)J=∫0t1eefeefT-γf^f^Tdt.
Then(29)J=∫0t1eefeefT-γf^f^T+V˙dt-∫0t1V˙dt≤∫0t1eefeefT-γf^f^T+V˙dt+V0=∫0t1eefAefTP1+P1Aef+IeefT-eefTP1Beff˙-γf^f^Tdt+V0=∫0t1eeff^AefTP1+P1Aef+IP1BefBefTP1-γI1eeff^Tdt+V0.
If there exists(30)AefTP1+P1Aef+IP1BefBefTP1-γI1<0,
then theH∞ tracking performance can be satisfied.
Next, note symmetric positive-definite matrix(31)P1=P00Q.
Thus(32)P1Bef=0-QS,AefTP1+P1Aef+I=A11A12A12TA22,
where A11, A12, and A22 are defined as (23); then
(33)∫0t1eefeefTdt≤∫0t1γf^f^Tdt+V0.
Then (24) can be obtained by simplifying (33); thus it is verified that unknown disturbance observer error ef of PIO can converge to zero in finite time, as well as the estimated interference f^t converging to actual interference ft.


#### 3.2.2. PID Neural Network Structure and Calculation Method

As it is difficult to acquire the physical parameters and motion information of small bodies accurately, there exists a highly nonlinear dynamic model of the small body. PID neural network control algorithm not only has the advantages of conventional PID controller, but also owns parallel structure and function of learning and memory of neural network and the ability of multilayer networks to approximate arbitrary functions. Therefore, the algorithm shows good superiority and stability performances in the control for the probe.The PID neural network is introduced as follows. The PID neural network is a three-forward neural network. Suppose that the controlled object has three inputs and three outputs, which is a nonlinear and strong coupling system with three variables. There exists a three-layer neural network comprising proportional neurons, integral neurons, and derivative neurons between the input layer and hidden layers. In addition, connected weights exist between the hidden layer and output layer. Figure2 shows a multivariable control structure based on a PID neural network of probe power decline period.Figure 2
Multivariable control structure based on a PID neural network of probe power decline period.(1) PID Neural Network Forward Algorithm. At any sampling time k, the forward calculation equations of the PID neural network are as follows.(a) The input-output function of input-layer neurons is(34)xsik=usik,
where xsi, ss=1,2,3, and ii=1,2 are input values of input-layer neurons, output values of input-layer neurons, and the number of the subnet input layers, respectively.Define the position error inz-axis orientation as e; then
(35)e=zt-znt,
where zt and znt are the actual position on the za-axis at time t and the nominal position on the za-axis at corresponding time t, respectively.Introduce a simple filters as new state variable, and the input of input layer is defined as follows:
(36)s=e˙+λe,
where λ is a positive scalar.(b) Hidden layer contains nine neurons (three proportional neurons, three integral neurons, and three derivative neurons); the input values of these neurons can be calculated as follows:(37)netsi(k)=∑i=12ωsijxsik,j=1,2,3,….For subnetworki, the formula of the output of hidden layer neurons is given by
(38)us1k=nets1k,us2k=nets2k+us2k-1,us3k=nets3k-nets3k-1,
where nets1(k), usj(k), ωsij, and j are input value of neurons in the hidden layer, the output value of neurons in the hidden layer, weight between input layer and hidden layer in each subnet, and the hidden layer neuron number in the subnet (j=1,2,3), respectively.(c) The input and output of output-layer neurons: the output of output-layer neurons is the sum of output weights of all hidden layer neurons as(39)yhk=∑s=1n∑j=13ωsjhusjk,
where yhk, ωsjh, and s are output value of output-layer neurons, connected weight between hidden layer and output layer, and sequence number of output-layer neurons (s=1,2,3), respectively.(2) PID Neural Network Learning Algorithm. In this subsection, a multivariable probe control system based on the PID neural network algorithm is regarded as a generalized network, using the backpropagation (BP) learning algorithm to minimize the criterion function within the scope of the requirements. Criterion function is given by [25]
(40)J=E=∑p=1nEp=12∑p=1nrpk-ypk2=12∑pe2k≤ε.The weight of the PID neural network can be adjusted by virtue of the gradient method, trained and learned throughk steps, and then determined depending on the following equation.(a) The iterative equation of weight between input layer and hidden layer is(41)ωijk+1=ωijk-η∂J∂ωij+η1ωijk-ωijk+1.(b) The iterative equation of weight between hidden layer and output layer is(42)ωjhk+1=ωjh(k)-η∂J∂ωjh+η2ωjhk-ωjhk+1.Proof.
Choose the Lyapunov candidate as(43)VJ,ω=αJ+12β∂J∂ω2,
where the symbol · signifies the quadratic norm of the vector and the parameters α, β are strict positive constants and they are utilized to determine the degree of proportion; then
(44)∂J∂ωl2=∂J∂ωl∂J∂ωT=∂J∂ω∂J∂ωT,
where ∂J/∂ω=∂J/∂ω1,…,∂J/∂ωl(45)∂J∂ωT=∂J∂ω1,…,∂J∂ωl.
Clearly,J and ∂J/∂ω2 are equal to the minimal neighborhood; in addition, the parameters α>0, β>0, so the function VJ,ω is positive definite.
Taking the derivative ofVJ,ω with respect to time, one can obtain
(46)V˙=α∂J∂ωω˙+∂J∂xx˙+∂J∂t+β∂J∂ω∂2J∂t∂ωT+∂2J∂ω∂ωTω˙+∂2J∂x∂ωTx˙=∂J∂ωαIl+β∂2J∂ω∂ωTω˙+α∂J∂x+β∂J∂ω∂2J∂x∂ωTx˙=∂J∂ωαIl+β∂2J∂ω∂ωTω˙+α∂J∂t+β∂J∂ω∂2J∂t∂ωT=-u∂J∂ω2-VJ2,if=∂J∂ω≠0α∂J∂t,if=∂J∂ω=0.
Above all, considering the functionV is positive definite and another function V˙ the BP learning algorithm, it is certified that the BP learning algorithm has the internality of making the error converge to the minimum.


## 3.1. The Nominal Trajectory Planning

The desired descent altitude and velocity are planned in order to satisfy the requirements of soft landing on the surface of small bodies. The constraint condition is defined as [23]
(9)z˙0=z˙0,z(0)=z0,zτ=ρ,z˙τ=0,
where z0 and z˙0 denote the initial altitude and altitude change rate, z˙0 and z0 are the planned altitude and altitude change rate, and τ is the descent time. The cubic curve to satisfy the boundary condition is given by
(10)znt=z0+z1t+z2t2+z3t3,
where z0, z1, z2, and z3 are the cubic function coefficients.Using (9), the coefficients are determined. The descent curve is given by
(11)znt=z0+z˙0t-3z0+2z˙0τ-3ρτ2t2+2z0+z˙0τ-2ρτ3t3,
where ρ is the altitude of the landing site.Next, the time derivatives of (12) are given by
(12)z˙nt=z˙0-6z0+4z˙0τ-6ρτ2t+6z0+3z˙0τ-6ρτ3t2,(13)z¨nt=-6z0+4z˙0τ-6ρτ2+12z0+6z˙0τ-12ρτ3t.Similarly, the ideal nominal trajectories can be planned on the other two axis directions.


## 3.2. Control Law Design

### 3.2.1. Proportional Integral Observer Design

Let(14)xt=x,y,z,x˙,y˙,z˙T,ut=000uxt-Uxuyt-Uyuzt-UzT.Convert (4) into nonlinear time-invariant system as
(15)x˙t=Axt+But+ft,yt=Cxt,
where
(16)A=000100000010000001ω20002ω00ω20-2ω00000000,B=000b1b2b3,ft=000fdxfdyfdz,C=111111,
where ω is the instantaneous rotation vector of the small body and fdx, fdy, and fdz are the components of unmodelled perturbation accelerations mainly from the solar radiation pressure and the solar gravitation.Next, the PIO is designed as follows [24]:
(17)x^˙t=Ax^t+But+f^t+Kpyt-y^t,y^t=Cx^t,f^˙t=KIyt-y^t,
where Kp and KI are the observer gain matrix and the integral coefficient of estimated unknown disturbance, respectively.Note the state error and unknown disturbance error as(18)et=x^t-xt,eft=f^t-ft,
where x^ and f^ are the estimations of the state vector x and the unknown disturbance f, respectively.Using (17) and (18) the error dynamics are as follows:
(19)e˙t=A-KPCet+ft-f^t,e˙f=-KIC+KvCA-KpCet-KvCef-f˙t.The augmented estimator system could be rewritten as(20)e˙te˙ft=Aefeef+Beff˙t,
where
(21)Aef=A-KPC1-KIC-KvCA+KvCKpC0,Bef=01.Lemma 1 (see [25]).
Givenγ>0 and (20), if there exist symmetric matrixes P, Q and two matrices KP, KI of appropriate dimension as well as LMI such that
(22)A11A120A12TA22-QS0-QSTγI1<0,
where
(23)A11=PA-KpC+A-KpCTP+Ie,A12=P+QSKiC+KvA-KpCT,A22=-QSKvC-KiFS+KvCQST+Is
hold, then system (20) is stable and satisfies corresponding performance index.Theorem 2.
For a given positive constantγ>0 and (20), if there exist symmetric matrixes P, Q and two matrices KP, KI of appropriate dimension such that (22) LMI holds, then system (20) is robust and asymptotically stable and satisfies the performance index as follows:
(24)eef≤γf^2+V0,
where
(25)eef2=∫0t1eefeefTdt,f^2=∫0t1f^f^Tdt.Ie, Is, and I1 are unit matrixes with appropriate dimension.Proof.
Choose the Lyapunov candidate as(26)V=eefP1eefT;
then taking the derivative ofV with respect to time along the trajectories of (19), one obtains
(27)V˙=eefAefTP1+P1AefeefT+2eefTP1Beff^.
Now define performance indicators as follows:(28)J=∫0t1eefeefT-γf^f^Tdt.
Then(29)J=∫0t1eefeefT-γf^f^T+V˙dt-∫0t1V˙dt≤∫0t1eefeefT-γf^f^T+V˙dt+V0=∫0t1eefAefTP1+P1Aef+IeefT-eefTP1Beff˙-γf^f^Tdt+V0=∫0t1eeff^AefTP1+P1Aef+IP1BefBefTP1-γI1eeff^Tdt+V0.
If there exists(30)AefTP1+P1Aef+IP1BefBefTP1-γI1<0,
then theH∞ tracking performance can be satisfied.
Next, note symmetric positive-definite matrix(31)P1=P00Q.
Thus(32)P1Bef=0-QS,AefTP1+P1Aef+I=A11A12A12TA22,
where A11, A12, and A22 are defined as (23); then
(33)∫0t1eefeefTdt≤∫0t1γf^f^Tdt+V0.
Then (24) can be obtained by simplifying (33); thus it is verified that unknown disturbance observer error ef of PIO can converge to zero in finite time, as well as the estimated interference f^t converging to actual interference ft.


### 3.2.2. PID Neural Network Structure and Calculation Method

As it is difficult to acquire the physical parameters and motion information of small bodies accurately, there exists a highly nonlinear dynamic model of the small body. PID neural network control algorithm not only has the advantages of conventional PID controller, but also owns parallel structure and function of learning and memory of neural network and the ability of multilayer networks to approximate arbitrary functions. Therefore, the algorithm shows good superiority and stability performances in the control for the probe.The PID neural network is introduced as follows. The PID neural network is a three-forward neural network. Suppose that the controlled object has three inputs and three outputs, which is a nonlinear and strong coupling system with three variables. There exists a three-layer neural network comprising proportional neurons, integral neurons, and derivative neurons between the input layer and hidden layers. In addition, connected weights exist between the hidden layer and output layer. Figure2 shows a multivariable control structure based on a PID neural network of probe power decline period.Figure 2
Multivariable control structure based on a PID neural network of probe power decline period.(1) PID Neural Network Forward Algorithm. At any sampling time k, the forward calculation equations of the PID neural network are as follows.(a) The input-output function of input-layer neurons is(34)xsik=usik,
where xsi, ss=1,2,3, and ii=1,2 are input values of input-layer neurons, output values of input-layer neurons, and the number of the subnet input layers, respectively.Define the position error inz-axis orientation as e; then
(35)e=zt-znt,
where zt and znt are the actual position on the za-axis at time t and the nominal position on the za-axis at corresponding time t, respectively.Introduce a simple filters as new state variable, and the input of input layer is defined as follows:
(36)s=e˙+λe,
where λ is a positive scalar.(b) Hidden layer contains nine neurons (three proportional neurons, three integral neurons, and three derivative neurons); the input values of these neurons can be calculated as follows:(37)netsi(k)=∑i=12ωsijxsik,j=1,2,3,….For subnetworki, the formula of the output of hidden layer neurons is given by
(38)us1k=nets1k,us2k=nets2k+us2k-1,us3k=nets3k-nets3k-1,
where nets1(k), usj(k), ωsij, and j are input value of neurons in the hidden layer, the output value of neurons in the hidden layer, weight between input layer and hidden layer in each subnet, and the hidden layer neuron number in the subnet (j=1,2,3), respectively.(c) The input and output of output-layer neurons: the output of output-layer neurons is the sum of output weights of all hidden layer neurons as(39)yhk=∑s=1n∑j=13ωsjhusjk,
where yhk, ωsjh, and s are output value of output-layer neurons, connected weight between hidden layer and output layer, and sequence number of output-layer neurons (s=1,2,3), respectively.(2) PID Neural Network Learning Algorithm. In this subsection, a multivariable probe control system based on the PID neural network algorithm is regarded as a generalized network, using the backpropagation (BP) learning algorithm to minimize the criterion function within the scope of the requirements. Criterion function is given by [25]
(40)J=E=∑p=1nEp=12∑p=1nrpk-ypk2=12∑pe2k≤ε.The weight of the PID neural network can be adjusted by virtue of the gradient method, trained and learned throughk steps, and then determined depending on the following equation.(a) The iterative equation of weight between input layer and hidden layer is(41)ωijk+1=ωijk-η∂J∂ωij+η1ωijk-ωijk+1.(b) The iterative equation of weight between hidden layer and output layer is(42)ωjhk+1=ωjh(k)-η∂J∂ωjh+η2ωjhk-ωjhk+1.Proof.
Choose the Lyapunov candidate as(43)VJ,ω=αJ+12β∂J∂ω2,
where the symbol · signifies the quadratic norm of the vector and the parameters α, β are strict positive constants and they are utilized to determine the degree of proportion; then
(44)∂J∂ωl2=∂J∂ωl∂J∂ωT=∂J∂ω∂J∂ωT,
where ∂J/∂ω=∂J/∂ω1,…,∂J/∂ωl(45)∂J∂ωT=∂J∂ω1,…,∂J∂ωl.
Clearly,J and ∂J/∂ω2 are equal to the minimal neighborhood; in addition, the parameters α>0, β>0, so the function VJ,ω is positive definite.
Taking the derivative ofVJ,ω with respect to time, one can obtain
(46)V˙=α∂J∂ωω˙+∂J∂xx˙+∂J∂t+β∂J∂ω∂2J∂t∂ωT+∂2J∂ω∂ωTω˙+∂2J∂x∂ωTx˙=∂J∂ωαIl+β∂2J∂ω∂ωTω˙+α∂J∂x+β∂J∂ω∂2J∂x∂ωTx˙=∂J∂ωαIl+β∂2J∂ω∂ωTω˙+α∂J∂t+β∂J∂ω∂2J∂t∂ωT=-u∂J∂ω2-VJ2,if=∂J∂ω≠0α∂J∂t,if=∂J∂ω=0.
Above all, considering the functionV is positive definite and another function V˙ the BP learning algorithm, it is certified that the BP learning algorithm has the internality of making the error converge to the minimum.


## 3.2.1. Proportional Integral Observer Design

Let(14)xt=x,y,z,x˙,y˙,z˙T,ut=000uxt-Uxuyt-Uyuzt-UzT.Convert (4) into nonlinear time-invariant system as
(15)x˙t=Axt+But+ft,yt=Cxt,
where
(16)A=000100000010000001ω20002ω00ω20-2ω00000000,B=000b1b2b3,ft=000fdxfdyfdz,C=111111,
where ω is the instantaneous rotation vector of the small body and fdx, fdy, and fdz are the components of unmodelled perturbation accelerations mainly from the solar radiation pressure and the solar gravitation.Next, the PIO is designed as follows [24]:
(17)x^˙t=Ax^t+But+f^t+Kpyt-y^t,y^t=Cx^t,f^˙t=KIyt-y^t,
where Kp and KI are the observer gain matrix and the integral coefficient of estimated unknown disturbance, respectively.Note the state error and unknown disturbance error as(18)et=x^t-xt,eft=f^t-ft,
where x^ and f^ are the estimations of the state vector x and the unknown disturbance f, respectively.Using (17) and (18) the error dynamics are as follows:
(19)e˙t=A-KPCet+ft-f^t,e˙f=-KIC+KvCA-KpCet-KvCef-f˙t.The augmented estimator system could be rewritten as(20)e˙te˙ft=Aefeef+Beff˙t,
where
(21)Aef=A-KPC1-KIC-KvCA+KvCKpC0,Bef=01.Lemma 1 (see [25]).
Givenγ>0 and (20), if there exist symmetric matrixes P, Q and two matrices KP, KI of appropriate dimension as well as LMI such that
(22)A11A120A12TA22-QS0-QSTγI1<0,
where
(23)A11=PA-KpC+A-KpCTP+Ie,A12=P+QSKiC+KvA-KpCT,A22=-QSKvC-KiFS+KvCQST+Is
hold, then system (20) is stable and satisfies corresponding performance index.Theorem 2.
For a given positive constantγ>0 and (20), if there exist symmetric matrixes P, Q and two matrices KP, KI of appropriate dimension such that (22) LMI holds, then system (20) is robust and asymptotically stable and satisfies the performance index as follows:
(24)eef≤γf^2+V0,
where
(25)eef2=∫0t1eefeefTdt,f^2=∫0t1f^f^Tdt.Ie, Is, and I1 are unit matrixes with appropriate dimension.Proof.
Choose the Lyapunov candidate as(26)V=eefP1eefT;
then taking the derivative ofV with respect to time along the trajectories of (19), one obtains
(27)V˙=eefAefTP1+P1AefeefT+2eefTP1Beff^.
Now define performance indicators as follows:(28)J=∫0t1eefeefT-γf^f^Tdt.
Then(29)J=∫0t1eefeefT-γf^f^T+V˙dt-∫0t1V˙dt≤∫0t1eefeefT-γf^f^T+V˙dt+V0=∫0t1eefAefTP1+P1Aef+IeefT-eefTP1Beff˙-γf^f^Tdt+V0=∫0t1eeff^AefTP1+P1Aef+IP1BefBefTP1-γI1eeff^Tdt+V0.
If there exists(30)AefTP1+P1Aef+IP1BefBefTP1-γI1<0,
then theH∞ tracking performance can be satisfied.
Next, note symmetric positive-definite matrix(31)P1=P00Q.
Thus(32)P1Bef=0-QS,AefTP1+P1Aef+I=A11A12A12TA22,
where A11, A12, and A22 are defined as (23); then
(33)∫0t1eefeefTdt≤∫0t1γf^f^Tdt+V0.
Then (24) can be obtained by simplifying (33); thus it is verified that unknown disturbance observer error ef of PIO can converge to zero in finite time, as well as the estimated interference f^t converging to actual interference ft.


## 3.2.2. PID Neural Network Structure and Calculation Method

As it is difficult to acquire the physical parameters and motion information of small bodies accurately, there exists a highly nonlinear dynamic model of the small body. PID neural network control algorithm not only has the advantages of conventional PID controller, but also owns parallel structure and function of learning and memory of neural network and the ability of multilayer networks to approximate arbitrary functions. Therefore, the algorithm shows good superiority and stability performances in the control for the probe.The PID neural network is introduced as follows. The PID neural network is a three-forward neural network. Suppose that the controlled object has three inputs and three outputs, which is a nonlinear and strong coupling system with three variables. There exists a three-layer neural network comprising proportional neurons, integral neurons, and derivative neurons between the input layer and hidden layers. In addition, connected weights exist between the hidden layer and output layer. Figure2 shows a multivariable control structure based on a PID neural network of probe power decline period.Figure 2
Multivariable control structure based on a PID neural network of probe power decline period.(1) PID Neural Network Forward Algorithm. At any sampling time k, the forward calculation equations of the PID neural network are as follows.(a) The input-output function of input-layer neurons is(34)xsik=usik,
where xsi, ss=1,2,3, and ii=1,2 are input values of input-layer neurons, output values of input-layer neurons, and the number of the subnet input layers, respectively.Define the position error inz-axis orientation as e; then
(35)e=zt-znt,
where zt and znt are the actual position on the za-axis at time t and the nominal position on the za-axis at corresponding time t, respectively.Introduce a simple filters as new state variable, and the input of input layer is defined as follows:
(36)s=e˙+λe,
where λ is a positive scalar.(b) Hidden layer contains nine neurons (three proportional neurons, three integral neurons, and three derivative neurons); the input values of these neurons can be calculated as follows:(37)netsi(k)=∑i=12ωsijxsik,j=1,2,3,….For subnetworki, the formula of the output of hidden layer neurons is given by
(38)us1k=nets1k,us2k=nets2k+us2k-1,us3k=nets3k-nets3k-1,
where nets1(k), usj(k), ωsij, and j are input value of neurons in the hidden layer, the output value of neurons in the hidden layer, weight between input layer and hidden layer in each subnet, and the hidden layer neuron number in the subnet (j=1,2,3), respectively.(c) The input and output of output-layer neurons: the output of output-layer neurons is the sum of output weights of all hidden layer neurons as(39)yhk=∑s=1n∑j=13ωsjhusjk,
where yhk, ωsjh, and s are output value of output-layer neurons, connected weight between hidden layer and output layer, and sequence number of output-layer neurons (s=1,2,3), respectively.(2) PID Neural Network Learning Algorithm. In this subsection, a multivariable probe control system based on the PID neural network algorithm is regarded as a generalized network, using the backpropagation (BP) learning algorithm to minimize the criterion function within the scope of the requirements. Criterion function is given by [25]
(40)J=E=∑p=1nEp=12∑p=1nrpk-ypk2=12∑pe2k≤ε.The weight of the PID neural network can be adjusted by virtue of the gradient method, trained and learned throughk steps, and then determined depending on the following equation.(a) The iterative equation of weight between input layer and hidden layer is(41)ωijk+1=ωijk-η∂J∂ωij+η1ωijk-ωijk+1.(b) The iterative equation of weight between hidden layer and output layer is(42)ωjhk+1=ωjh(k)-η∂J∂ωjh+η2ωjhk-ωjhk+1.Proof.
Choose the Lyapunov candidate as(43)VJ,ω=αJ+12β∂J∂ω2,
where the symbol · signifies the quadratic norm of the vector and the parameters α, β are strict positive constants and they are utilized to determine the degree of proportion; then
(44)∂J∂ωl2=∂J∂ωl∂J∂ωT=∂J∂ω∂J∂ωT,
where ∂J/∂ω=∂J/∂ω1,…,∂J/∂ωl(45)∂J∂ωT=∂J∂ω1,…,∂J∂ωl.
Clearly,J and ∂J/∂ω2 are equal to the minimal neighborhood; in addition, the parameters α>0, β>0, so the function VJ,ω is positive definite.
Taking the derivative ofVJ,ω with respect to time, one can obtain
(46)V˙=α∂J∂ωω˙+∂J∂xx˙+∂J∂t+β∂J∂ω∂2J∂t∂ωT+∂2J∂ω∂ωTω˙+∂2J∂x∂ωTx˙=∂J∂ωαIl+β∂2J∂ω∂ωTω˙+α∂J∂x+β∂J∂ω∂2J∂x∂ωTx˙=∂J∂ωαIl+β∂2J∂ω∂ωTω˙+α∂J∂t+β∂J∂ω∂2J∂t∂ωT=-u∂J∂ω2-VJ2,if=∂J∂ω≠0α∂J∂t,if=∂J∂ω=0.
Above all, considering the functionV is positive definite and another function V˙ the BP learning algorithm, it is certified that the BP learning algorithm has the internality of making the error converge to the minimum.


## 4. Simulation Results

(a) According to Theorem2, PIO parameters can be derived by using the LMI toolbox as follows: (47)Kp=1.8820-11.21214.9258-11.2727-14.75485.87730.55219.56234.92583.77460.2464-7.3346-12.34703.22469.275825.08644.908428.90355.2324-12.38699.2758-6.9024-3.7824-2.93573.419624.1648-4.973513.93405.983512.3925-5.20152.97387.17481.38699.0247-8.0924,KI=4.377812.9673-1.37829.0357-14.78931.5044.The initial values of proportional neurons, integral neurons, and three derivative neurons are implemented as follows:(48)ωs1j=0.1,ωs2j=-0.1,ωs3j=0.2.The initial values of connected weight between hidden layer and output layer are defined, respectively, as follows:(49)ωsj1=0.4,ωsj2=-0.5,ωsj3=0.8.(b) The asteroid Eros 433 is taken as the target small body for simulation to verify the feasibility of the presented control scheme. The parameters of the small body are gained from [26] and are shown in  Table 1.Table 1
Parameters
Real world
Simulation parameters
μ
4
.
749
E
-
04
4.800
E
-
04
Spin period (h)
10.54
10.55
R
0 (km)
1.150
1.148
Gravitational coefficients
C20
−0.113
−0.110
C22
0.0396
0.0397
Initial position (m)
[400, 400, 11000]
Initial velocity (m/s)
[−0.9, 1.2,  −1]
Terminal site
[100, 100, 8000]In this paper, compared with the perturbation uncertainties proposed in [27], the larger perturbation uncertainties are chosen as follows:
(50)fdx=150sin2t,fdy=160sin1.5t,fdz=140sin3t.From Figure3, it can be seen that the actual trajectory of probe exhibits evident chattering in the system suffering a larger disturbance. Inherent robustness of the sliding mode control algorithm is not sufficient to guarantee the actual trajectory to track the desired trajectory. On this occasion, the neural network control algorithm based on PIO in this paper is utilized to compensate the unknown disturbance and eliminate the chattering problem of trajectory. Meanwhile, it can track the required ideal location quickly.Figure 3
Landing trajectory curve of the probe.Figures4, 5, 6, 7, 8, and 9 show the error curves between ideal location and actual locations and the velocity curves as a function of time on the three axis directions. For the system exhibiting large initial error and perturbation uncertainties, on the condition that the convergence of the system can be ensured, compared with the sliding mode control algorithm [27] the neural network control algorithm based on PIO can improve the convergence rate of the actual position error and the actual velocity; namely, the actual trajectory can fast and accurately track the planned trajectory on the condition that there exist parameter, feedback state error, and external larger disturbance in the system. Thereby it can satisfy the probe to land smoothly on the surface of small body and avoid the occurrence of the probe crashing due to the excessive landing speed.Figure 4
Position error component (x-axis) as function of time.Figure 5
Position error component (y-axis) as function of time.Figure 6
Position error component (z-axis) as function of time.Figure 7
Velocity component (x-axis) as function of time.Figure 8
Velocity component (y-axis) as function of time.Figure 9
Velocity component (z-axis) as function of time.


## 5. Conclusion

This paper has presented a neural network control algorithm based on PIO. In view of the power descent section of soft landing on small bodies, the system dynamic models of the small bodies under the body-fixed coordinate system are given with ignoring the attitude control. The solar radiation pressure and the third-body’s gravity are treated as the perturbation, which is viewed as a bounded function. The nominal trajectories meeting the constraints on the threeaxes are preplanned. The simulation results show that the neural network control algorithm based on PIO can ensure fast and accurate response to parameter uncertainty, feedback state error, and external disturbances. Moreover, for the system exhibiting larger interference, it can overcome the inherent chattering problem of sliding mode control algorithm and make the position error and the velocity error converge to the small finite value, realizing the aim to softly land.



---

*Source: 102424-2015-02-22.xml*