CN102890506B - Small body approaching section guidance control method based on constraint programming - Google Patents

Abstract

The invention discloses a small body approaching section guidance control method based on constraint programming. A target body gravity model is led to an orbit planning process, and simultaneously state and control are constrained in the approaching process; and a non-linear dynamic path planning problem with control constraint and orbit constraint is converted into a second-order cone programming problem (SOCP) taking optimizing burn-up as a performance index, so that the convex planning problem can be efficiently resolved through an interior point method. A path point generated by programming is not merely a target tracking state in the traditional significance, and the state meets a feasible solution of a dynamic system with state and control constraint under the situation of considering the effect of a target body gravity field. In addition, the obtained feedforward control input meets the requirements of closing time of a thruster required for an onboard device and filter estimation.

Description

A kind of small feature loss Approach phase Guidance and control method based on constraint planning

Technical field

The present invention relates to a kind of Guidance and control method to small celestial body exploration device, particularly a kind of Guidance and control method meeting state and control constraints in close to small feature loss process.

Background technology

To small feature loss close in process, due to gravitation, geometry landform and the probabilistic existence of correlative factor thereof, detector Guidance and control algorithm needs to have certain processing power to X factor.Usually for the detection mission of planet and the moon, due to the long-time observation to target celestial body, obtain comparatively complete celestial body characteristic information in advance, therefore become the normal mode it carried out in Proximity operation process based on the Guidance and control strategy of land station.And for small celestial body exploration task, do not have the condition of it being carried out to long-term observation at present, this just needs detector self to possess the ability of unknown situation being carried out to fast processing, in addition the existence of longer communication delay, makes utilization greatly reduce based on the possibility of the Guidance and control scheme of land station's pattern.Self-contained guidance control algolithm is conducive to properly settling this problem, and current spaceborne computing power and available calculation technology make this method become possibility.

Prior art is see R.R.Sostaric, J.R.Rea.Powered descent guidance methods for themoon and mars.San Francisco, USA:American Institute of Aeronautics and AstronauticsInc, 2005., traditional method for planning track utilizes polynomial expression carry out matching to detector current state and expectation state thus form track path point, using it as the analytical form of approaching system optimization solution.Method is simple and operand is low due to this, therefore become in the past detector close to the selection of task section Guidance and control strategy.But the solution of nonlinear detector dynamical system is difficult to replace by simple polynomial form, and the path point formed by traditional method for planning track is not the optimum solution of system, simultaneously, what require with detection along with the complicated of detection mission improves constantly, tradition method for planning track cannot meet state constraint strict in task process and control constraints, therefore the track Guidance and control method finding optimization becomes the problem of recent domestic focus of attention, see S.R.Ploen, A.B.Acikmese, A.Wolf.A comparison of powered descentguidance laws for Mars pinpoint landing.Reston, VA, USA:American Institute ofAeronautics and Astronautics Inc, 2006 pairs of multiple planet method of guidances have carried out detailed comparisons's analysis.

Summary of the invention

The present invention is directed to the defect of small feature loss Approach phase task tradition fitting of a polynomial Guidance and control algorithm, Dynamic Constraints, state constraint and control constraints are introduced in the process of trajectory planning simultaneously, produce by solving convex programming problem the feasible solution meeting and there is the Kind of Nonlinear Dynamical System of state and control constraints, thus improve small feature loss Approach phase Guidance and control performance.

This kind, based on the small feature loss Approach phase Guidance and control method of constraint planning, specifically comprises the following steps:

The first step: set up detector close to small feature loss kinetic model;

Second step: system model is carried out discretize, and gained model is used for the guidance of follow-up convex programming;

3rd step: according to detector current state, end expectation state and expectation time kept in reserve, obtain an initial reference track by linear or fitting of a polynomial;

4th step: by this initial reference track, introduces state constraint and control constraints, solves to have the kinetic pathways that control constraints and track state retrain and plan by iterating, and obtains close to guided paths feasible in target celestial body process.

Wherein first step medium power model is by following the Representation Equation

$\stackrel{..}{r}+\stackrel{.}{\mathrm{\ω}}\×r+2\mathrm{\ω}\×\stackrel{.}{r}+\mathrm{\ω}\×(\mathrm{\ω}\×r)=u+d+g\left(r\right)---\left(1\right)$

In formula, r ∈ R ³for the radius vector of detector under small feature loss barycenter is connected coordinate system; ω ∈ R ³for celestial body angular spin rate; U ∈ R ³for the control acceleration of effect; D is disturbing acceleration; G is gravitational acceleration; Suppose that small feature loss angular spin rate is stable and be a constant value, namely kinetic model then in (1) formula is by following state space description form

$x={\left[\begin{array}{cc}r& \stackrel{.}{r}\end{array}\right]}^T---\left(2\right)$

$\stackrel{.}{x}=\mathrm{Ax}+\mathrm{Bu}+\mathrm{Bg}\left(C_{r}x\right)---\left(3\right)$

In formula, $A=\left[\begin{array}{cc}{\left[\begin{array}{cc}0& -{\widehat{\mathrm{\ω}}}^2\end{array}\right]}^T& {\left[\begin{array}{cc}I& -2\widehat{\mathrm{\ω}}\end{array}\right]}^T\end{array}\right];$ B=[0 I] ^t; C _r=[0 I] and for the matrix representation of vectorial multiplication cross ω × ().

Wherein in second step, discretize adopts following methods:

At Fixed Time Interval Δ t, wherein Δ t>=δ _f+ δ _sin, control acceleration input u (t) and be described as

$u\left(t\right)=\left\{\begin{array}{cc}{u}_k& t\∈[t_k,t_k+{\mathrm{\δ}}_f]\\ 0& t\∈(t_k+{\mathrm{\δ}}_f,t_{k+1})\end{array}\right.---\left(4\right)$

In formula, t _k+1-t _k=Δ t and u _kfor limited control inputs, in thruster opening process, the solution of the described system of formula (3) is at t=t _k+ δ _fmoment is

$x(t_k+{\mathrm{\δ}}_f)=e^{A{\mathrm{\δ}}_f}x\left(t_k\right)+{\∫}_0^{{\mathrm{\δ}}_f}{e}^{A({\mathrm{\δ}}_f-\mathrm{\τ})}\mathrm{Bd\τ}\·u_k+{\∫}_0^{{\mathrm{\δ}}_f}{e}^{A({\mathrm{\δ}}_f-\mathrm{\τ})}\mathrm{Bd}{\mathrm{\τ}\·g}_k---\left(5\right)$

In formula, g _kfor celestial body gravitation potential function is at state x _kon gradient, namely using numerical value the Gradient of Gravitation as gravitational acceleration input item in discrete model; In thruster closing process subsequently, the solution of the described system of formula (3) is at t=t _k+1moment is

$x\left(t_{k+1}\right)=e^{A(\mathrm{\Δt}-{\mathrm{\δ}}_f)}x(t_k+{\mathrm{\δ}}_f)+{\∫}_{{\mathrm{\δ}}_f}^{\mathrm{\Δt}}{e}^{A(\mathrm{\Δt}-\mathrm{\τ})}\mathrm{Bd\τ}\·g_k---\left(6\right)$

Formula (5) is substituted into formula (6),

$x\left(t_{k+1}\right)=e^{\mathrm{A\Δt}}x\left(t_k\right)+e^{A(\mathrm{\Δt}-{\mathrm{\δ}}_f)}{\∫}_0^{{\mathrm{\δ}}_f}{e}^{A({\mathrm{\δ}}_f-\mathrm{\τ})}\mathrm{Bd\τ}\·u_k+{\∫}_0^{\mathrm{\Δt}}{e}^{A(\mathrm{\Δt}-\mathrm{\τ})}\mathrm{Bd}{\mathrm{\τ}\·g}_k---\left(7\right)$

By x (t _k) referred to as x _k, each time interval Δ t has following discrete equation form

x _k+1＝A _kx _k+B _ku _k+E _kg _k(8)

In formula, A _k=e ^{a Δ t}; $B_k=e^{A(\mathrm{\Δt}-{\mathrm{\δ}}_f)}{\∫}_0^{{\mathrm{\δ}}_f}{e}^{A({\mathrm{\δ}}_f-\mathrm{\τ})}\mathrm{Bd\τ};$ $E_k={\∫}_0^{\mathrm{\Δt}}{e}^{A(\mathrm{\Δt}-\mathrm{\τ})}\mathrm{Bd\τ}.$

Wherein in the 4th step, matching obtains initial reference track to determine k=0 ..., N-1; This matching track meets at initial time t=0 (k=0) x _sfor current state, end moment t=T, k=N, meet x _ffor end expectation state; This solution procedure is described as the process of following iterative convex programming problem:

Given current state x _s: measured value or estimated value, expect SOT state of termination x _fand initial reference track perform j=0 ..., M-1 walks iteration:

(1) guided paths that jth time iteration produces is utilized calculate and upgrade discrete time-varying model parameter $g_k^j(k=0,\·\·\·,N-1);$

(2) following SOCP problem is solved to produce guided paths with corresponding control sequence

$\left.\begin{array}{c}\mathrm{Minimize}\underset{k=0}{\overset{N-1}{\mathrm{\Σ}}}\left(\mathrm{\α}\right|\left|u_k^{j+1}\right||+\mathrm{\β}{\left|\right|u_k^{j+1}\left|\right|}^2)+\mathrm{\ϵ}\left|\right|E_v(x_N^{j+1}-x_F)\left|\right|\\ \mathrm{subject\; to}\\ x_{k+1}^{j+1}=A_k{x}_k^{j+1}+B_k{u}_k^{j+1}+E_k{g}_k^j\\ x_0^{j+1}=x_S\\ E_r(x_N^{j+1}-x_F)=0\\ x_k^{j+1}\∈X\\ u_k^{j+1}\∈U\\ \left|\right|x_k^{j+1}-x_k^j\left|\right|\≤\mathrm{\κ}\left|\right|x_k^j-x_k^{j-1}\left|\right|(j=1,\·\·\·,M-1)\end{array}\right\}---\left(9\right)$

In formula, in cost function, α, β are respectively the weight of burnup and energy consumption, and (α, β)=(1,0) is for optimizing burnup, and (α, β)=(0,1) is for optimizing energy consumption; Relax end-fixity, terminal velocity constraint is added cost function, gets ε > 0, E _r=[I 0], E _v=[0 I]; State constraint territory X and control quantity constraint territory U is the convex set utilizing second order circular cone constraint specification; Get constrained parameters 0 < κ < 1, to guarantee to be a Cauchy sequence at the status switch in each moment, so far obtained convex programming problem can utilize interior point method to carry out Efficient Solution.

Principle of work of the present invention: target celestial body gravity model is introduced in trajectory planning process, add simultaneously close in process to the constraint of state and control, be converted into optimize the second order circular cone that burnup is performance index plan (SOCP) problem by having nonlinear kinetics path planning problem that control constraints and track state retrain.And then this convex programming problem carries out Efficient Solution by interior point method to it.Plan that the path point of generation is not only traditional target following state thus, and this state is considering to meet the feasible solution with state and control constraints dynamical system under target celestial body gravitational field affects situation; In addition, the input of gained feedforward control meets on-board equipment and the thruster shut-in time requirement needed for filtering estimation.

Beneficial effect of the present invention:

Convex programming guidance algorithm gives the nominal trajectory that meets state constraint and control constraints simultaneously.Compared with traditional fitting of a polynomial guidance algorithm, due to Dynamic Constraints, state constraint and control constraints are introduced in the process of trajectory planning simultaneously, therefore the path point produced based on the guidance algorithm of convex programming meets the feasible solution with the Kind of Nonlinear Dynamical System of state and control constraints, thus improve Guidance and control performance.Meanwhile, due to the high efficiency that convex programming problem solves, this Guidance and control method has the potentiality of online planning, can meet small-sized spacecraft and cheap spacecraft due to computational resource limited and need the situation of online process in real time.

Embodiment

For making the object, technical solutions and advantages of the present invention clearly; below embodiments of the invention are elaborated: the present embodiment is implemented under premised on technical scheme of the present invention; give detailed embodiment and concrete operating process, but protection scope of the present invention is not limited to following embodiment.

Detector is set up close to small feature loss kinetic model under small feature loss barycenter is connected coordinate system, can by following the Representation Equation

$\stackrel{..}{r}+\stackrel{.}{\mathrm{\&omega;}}\&times;r+2\mathrm{\&omega;}\&times;\stackrel{.}{r}+\mathrm{\&omega;}\&times;(\mathrm{\&omega;}\&times;r)=u+d+g\left(r\right)---\left(1\right)$

In formula, r ∈ R ³for the radius vector of detector under small feature loss barycenter is connected coordinate system; ω ∈ R ³for celestial body angular spin rate; U ∈ R ³for the control acceleration of effect; D is disturbing acceleration (launching generation primarily of extraneous gravitation interference, solar light pressure, comet dust); G is gravitational acceleration.

Suppose that small feature loss angular spin rate is stable and be a constant value kinetic model then in (1) formula can have following state space description form

$x={\left[\begin{array}{cc}r& \stackrel{.}{r}\end{array}\right]}^T---\left(2\right)$

$\stackrel{.}{x}=\mathrm{Ax}+\mathrm{Bu}+\mathrm{Bg}\left(C_{r}x\right)---\left(3\right)$

In formula, $A=\left[\begin{array}{cc}{\left[\begin{array}{cc}0& -{\widehat{\mathrm{\&omega;}}}^2\end{array}\right]}^T& {\left[\begin{array}{cc}I& -2\widehat{\mathrm{\&omega;}}\end{array}\right]}^T\end{array}\right];$ B=[0 I] ^t; C _r=[0 I] and for the matrix representation of vectorial multiplication cross ω × ().It should be noted that here not to g (C _rx) carry out linearization process, numerical value the Gradient of Gravitation will be utilized in subsequent designs to replace; And disturbing acceleration d is processed by the uncertain factor as nominal system dynamics model (3) by FEEDBACK CONTROL.

Consider the limited opening time δ of thruster _fand shut-in time δ _srequirement, system (3) is carried out discretize, and gained model will be used to the design of convex programming guidance algorithm.At Fixed Time Interval Δ t (Δ t>=δ _f+ δ _s) in, control acceleration input u (t) and can be described as

$u\left(t\right)=\left\{\begin{array}{cc}{u}_k& t\&Element;[t_k,t_k+{\mathrm{\&delta;}}_f]\\ 0& t\&Element;(t_k+{\mathrm{\&delta;}}_f,t_{k+1})\end{array}\right.---\left(4\right)$

In formula, t _k+1-t _k=Δ t and u _kfor limited control inputs.In thruster opening process, the solution of the described system of formula (3) is at t=t _k+ δ _fmoment can be approximately

$x(t_k+{\mathrm{\&delta;}}_f)=e^{A{\mathrm{\&delta;}}_f}x\left(t_k\right)+{\&Integral;}_0^{{\mathrm{\&delta;}}_f}{e}^{A({\mathrm{\&delta;}}_f-\mathrm{\&tau;})}\mathrm{Bd\&tau;}\&CenterDot;u_k+{\&Integral;}_0^{{\mathrm{\&delta;}}_f}{e}^{A({\mathrm{\&delta;}}_f-\mathrm{\&tau;})}\mathrm{Bd}{\mathrm{\&tau;}\&CenterDot;g}_k---\left(5\right)$

In formula, g _kfor celestial body gravitation potential function is at state x _kon gradient (namely using numerical value the Gradient of Gravitation as gravitational acceleration input item in discrete model).In thruster closing process subsequently, the solution of the described system of formula (3) is at t=t _k+1moment can be approximately

$x\left(t_{k+1}\right)=e^{A(\mathrm{\&Delta;t}-{\mathrm{\&delta;}}_f)}x(t_k+{\mathrm{\&delta;}}_f)+{\&Integral;}_{{\mathrm{\&delta;}}_f}^{\mathrm{\&Delta;t}}{e}^{A(\mathrm{\&Delta;t}-\mathrm{\&tau;})}\mathrm{Bd\&tau;}\&CenterDot;g_k---\left(6\right)$

Formula (5) is substituted into formula (6), can obtain

$x\left(t_{k+1}\right)=e^{\mathrm{A\&Delta;t}}x\left(t_k\right)+e^{A(\mathrm{\&Delta;t}-{\mathrm{\&delta;}}_f)}{\&Integral;}_0^{{\mathrm{\&delta;}}_f}{e}^{A({\mathrm{\&delta;}}_f-\mathrm{\&tau;})}\mathrm{Bd\&tau;}\&CenterDot;u_k+{\&Integral;}_0^{\mathrm{\&Delta;t}}{e}^{A(\mathrm{\&Delta;t}-\mathrm{\&tau;})}\mathrm{Bd}{\mathrm{\&tau;}\&CenterDot;g}_k---\left(7\right)$

By x (t _k) referred to as x _k, each time interval Δ t has following discrete equation form

x _k+1＝A _kx _k+B _ku _k+E _kg _k(8)

In formula, A _k=e ^{a Δ t}; $B_k=e^{A(\mathrm{\&Delta;t}-{\mathrm{\&delta;}}_f)}{\&Integral;}_0^{{\mathrm{\&delta;}}_f}{e}^{A({\mathrm{\&delta;}}_f-\mathrm{\&tau;})}\mathrm{Bd\&tau;};$ $E_k={\&Integral;}_0^{\mathrm{\&Delta;t}}{e}^{A(\mathrm{\&Delta;t}-\mathrm{\&tau;})}\mathrm{Bd\&tau;}.$

In the discrete system model that formula (8) describes, due to g _kdepend on guided paths x _k, therefore first can obtain an initial reference track by linear or fitting of a polynomial to determine this matching track meets at initial time t=0 (k=0) (x _sfor current state), end moment t=T (k=N) meets (x _ffor end expectation state), again by this initial guess, introducing state constraint and control constraints, solving the kinetic pathways planning problem that there is control constraints and track state and retrain by iterating, to find close to guided paths { u feasible in target celestial body process _k, { x _k.This solution procedure can be described as the process of following iterative convex programming problem:

Given current state x _s(measured value or estimated value), expects SOT state of termination x _fand initial reference track perform j=0 ..., M-1 walks iteration:

(1) guided paths that jth time iteration produces is utilized calculate and upgrade discrete time-varying model parameter $g_k^j(k=0,\&CenterDot;\&CenterDot;\&CenterDot;,N-1);$

(2) following SOCP problem is solved to produce guided paths with corresponding control sequence

$\left.\begin{array}{c}\mathrm{Minimize}\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\left(\mathrm{\&alpha;}\right|\left|u_k^{j+1}\right||+\mathrm{\&beta;}{\left|\right|u_k^{j+1}\left|\right|}^2)+\mathrm{\&epsiv;}\left|\right|E_v(x_N^{j+1}-x_F)\left|\right|\\ \mathrm{subject\; to}\\ x_{k+1}^{j+1}=A_k{x}_k^{j+1}+B_k{u}_k^{j+1}+E_k{g}_k^j\\ x_0^{j+1}=x_S\\ E_r(x_N^{j+1}-x_F)=0\\ x_k^{j+1}\&Element;X\\ u_k^{j+1}\&Element;U\\ \left|\right|x_k^{j+1}-x_k^j\left|\right|\&le;\mathrm{\&kappa;}\left|\right|x_k^j-x_k^{j-1}\left|\right|(j=1,\&CenterDot;\&CenterDot;\&CenterDot;,M-1)\end{array}\right\}---\left(9\right)$

In formula, in cost function, α, β are respectively the weight of burnup and energy consumption, and (α, β)=(1,0) is for optimizing burnup, and (α, β)=(0,1) is for optimizing energy consumption; Calculate for the ease of process, relax end-fixity, terminal velocity constraint is added cost function, gets ε > 0, E _r=[I 0], E _v=[0 I]; State constraint territory X and control quantity constraint territory U is the convex set utilizing second order circular cone constraint specification; For ensureing that convex programming guidance algorithm produces the convergence of track, get constrained parameters 0 < κ < 1, to guarantee to be a Cauchy sequence at the status switch in each moment.So far obtained convex programming problem can utilize interior point method to carry out Efficient Solution.

For asteroid Eros Approach phase task, detector model parameter is as shown in table 1.Target celestial body quadravalence gravitational potential function model is such as formula shown in (10).Wherein, μ _a=GM is the product of universal gravitational constant and small feature loss quality; A is the nominal radius of small feature loss; R, θ, being respectively that small feature loss barycenter is connected is detector radial distance, right ascension and declination under coordinate system; C _mnfor the humorous term coefficient of ball.

Table 1 detector simulation parameters

Target celestial body Proximity operation Control and Guild target is in expectation time kept in reserve T=300s, makes detector under control strategy u+ δ u effect by current state x _s=[8,950 20-50 1.5 2 0] ^tmove to along path planning and expect end state x _f=[845 00000 0] ^t, and meet following state constraint and control constraints in the process:

Catalog of celestial bodies identity distance celestial body barycenter near I hypothesis approximated position is about 8445m, for avoiding detector and celestial body to bump against, and setting nominal state constraint set $X=\left\{x\right|c^{T}x\&GreaterEqual;1,c={\left[\begin{array}{cc}1/8445& 0_{1\&times;5}\end{array}\right]}^T\};$

The II control inputs amplitude upper limit is

III thruster opening time was δ _f=10s, the thruster shut-in time is δ _s=20s, due to Guidance and control time interval Δ t>=δ _f+ δ _s, therefore get Δ t=30s.

Utilize convex programming guidance algorithm gained feedforward control amount detector can be sent within the scope of the 2m of expectation approximated position point, final speed can be controlled in about 0.02m/s, and the control that guidance algorithm provides meets the thrust amplitude upper limit constraint of thruster within the limited opening time.

Claims (3)

1., based on a small feature loss Approach phase Guidance and control method for constraint planning, it is characterized in that, specifically comprise the following steps:

The first step: set up detector close to small feature loss kinetic model;

Second step: system model is carried out discretize, and gained model is used for the guidance of follow-up convex programming;

3rd step: according to detector current state, end expectation state and expectation time kept in reserve, obtain an initial reference track by linear or fitting of a polynomial;

4th step: by this initial reference track, introduces state constraint and control constraints, solves to have the kinetic pathways that control constraints and track state retrain and plan by iterating, and obtains close to guided paths feasible in target celestial body process;

Matching obtains initial reference track to determine g _kfor celestial body gravitation potential function is at state x _kon gradient, k=0 ..., N-1; This matching track meets at initial time t=0 (k=0) x _sfor current state, end moment t=T, k=N, meet x _ffor end expectation state; This solution procedure is described as the process of following iterative convex programming problem:

Given current state x _s: measured value or estimated value, expect SOT state of termination x _fand initial reference track perform j=0 ..., M-1 walks iteration:

(1) guided paths that jth time iteration produces is utilized calculate and upgrade discrete time-varying model parameter (k=0 ..., N-1);

(2) following SOCP problem is solved to produce guided paths with corresponding control sequence

$\left.\begin{array}{c}\mathrm{Minimize}\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\left(\mathrm{\&alpha;}\right|\left|u_k^{j+1}\right||+\mathrm{\&beta;}{\left|\right|u_k^{j+1}\left|\right|}^2)+\mathrm{\&epsiv;}\left|\right|E_v(x_N^{j+1}-x_F)\left|\right|\\ \mathrm{subject\; to}\\ x_{k+1}^{j+1}=A_k{x}_k^{j+1}+B_k{u}_k^{j+1}+E_k{g}_k^j\\ x_0^{j+1}=x_S\\ E_r(x_N^{j+1}-x_F)=0\\ x_k^{j+1}\&Element;X\\ u_k^{j+1}\&Element;U\\ \left|\right|x_k^{j+1}-x_k^j\left|\right|\&le;\mathrm{\&kappa;}\left|\right|x_k^j-x_k^{j-1}\left|\right|,(j=1,...,M-1)\end{array}\right\}---\left(9\right)$

In formula, in cost function, α, β are respectively the weight of burnup and energy consumption, and (α, β)=(1,0) is for optimizing burnup, and (α, β)=(0,1) is for optimizing energy consumption; Relax end-fixity, terminal velocity constraint is added cost function, gets ε > 0, E _r=[I 0], E _v=[0 I]; State constraint territory X and control quantity constraint territory U is the convex set utilizing second order circular cone constraint specification; Get constrained parameters 0 < κ < 1, to guarantee to be a Cauchy sequence at the status switch in each moment, so far obtained convex programming problem can utilize interior point method to carry out Efficient Solution.

2. a kind of based on retraining the small feature loss Approach phase Guidance and control method planned as claimed in claim 1, it is characterized in that, wherein first step medium power model is by following the Representation Equation

$\stackrel{\&CenterDot;\&CenterDot;}{r}+\stackrel{\&CenterDot;}{\mathrm{\&omega;}}\&times;r+2\mathrm{\&omega;}\&times;\stackrel{\&CenterDot;}{r}+\mathrm{\&omega;}\&times;(\mathrm{\&omega;}\&times;r)=u+d+g\left(r\right)---\left(1\right)$

In formula, r ∈ R ³for the radius vector of detector under small feature loss barycenter is connected coordinate system; ω ∈ R ³for celestial body angular spin rate; U ∈ R ³for the control acceleration of effect; D is disturbing acceleration; G is gravitational acceleration; Suppose that small feature loss angular spin rate is stable and be a constant value, namely kinetic model then in (1) formula is by following state space description form

$x={\left[\begin{array}{cc}r& \stackrel{\&CenterDot;}{r}\end{array}\right]}^T---\left(2\right)$

$\stackrel{\&CenterDot;}{x}=\mathrm{Ax}+\mathrm{Bu}+\mathrm{Bg}\left(C_{r}x\right)---\left(3\right)$

In formula, $A=\left[\begin{array}{cc}{\left[\begin{array}{cc}0& -{\widehat{\mathrm{\&omega;}}}^2\end{array}\right]}^T& {\left[\begin{array}{cc}I& -2\widehat{\mathrm{\&omega;}}\end{array}\right]}^T\end{array}\right];$ B=[0 I] ^t; C _r=[0 I] and for the matrix representation of vectorial multiplication cross ω × ().

3. as claimed in claim 1 a kind of based on constraint planning small feature loss Approach phase Guidance and control method, it is characterized in that, wherein in second step discretize adopt following methods:

At Fixed Time Interval Δ t, wherein Δ t>=δ _f+ δ _sin, control acceleration input u (t) and be described as

$u\left(t\right)=\left\{\begin{array}{cc}{u}_k& t\&Element;[t_k,t_k+{\mathrm{\&delta;}}_f]\\ 0& t\&Element;(t_k+{\mathrm{\&delta;}}_f,t_{k+1})\end{array}\right.----\left(4\right)$

In formula, t _k+1-t _k=Δ t and u _kfor limited control inputs, in thruster opening process, the solution of the described system of formula (3) is at t=t _k+ δ _fmoment is

$x(t_k+{\mathrm{\&delta;}}_f)=e^{A{\mathrm{\&delta;}}_f}x\left(t_k\right)+{\&Integral;}_0^{{\mathrm{\&delta;}}_f}{e}^{A({\mathrm{\&delta;}}_f-\mathrm{\&tau;})}\mathrm{Bd\&tau;}\&CenterDot;u_k+{\text{\&Integral;}}_0^{{\mathrm{\&delta;}}_f}{e}^{A({\mathrm{\&delta;}}_f-\mathrm{\&tau;})}\mathrm{Bd\&tau;}\&CenterDot;g_k---\left(5\right)$

In formula, g _kfor celestial body gravitation potential function is at state x _kon gradient, namely using numerical value the Gradient of Gravitation as gravitational acceleration input item in discrete model; In thruster closing process subsequently, the solution of the described system of formula (3) is at t=t _k+1moment is

$x\left(t_{k+1}\right)=e^{A(\mathrm{\&Delta;t}-{\mathrm{\&delta;}}_f)}x(t_k+{\mathrm{\&delta;}}_f)+{\&Integral;}_{{\mathrm{\&delta;}}_f}^{\mathrm{\&Delta;t}}{e}^{A(\mathrm{\&Delta;t}-\mathrm{\&tau;})}\mathrm{Bd\&tau;}\&CenterDot;g_k---\left(6\right)$

Formula (5) is substituted into formula (6),

$x\left(t_{k+1}\right)=e^{\mathrm{A\&Delta;t}}x\left(t_k\right)+e^{A(\mathrm{\&Delta;t}-{\mathrm{\&delta;}}_f)}{{\&Integral;}_0^{{\mathrm{\&delta;}}_f}e}^{A({\mathrm{\&delta;}}_f-\mathrm{\&tau;})}\mathrm{Bd\&tau;}\&CenterDot;u_k+{\text{\&Integral;}}_0^{\mathrm{\&Delta;t}}{e}^{A(\mathrm{\&Delta;t}-\mathrm{\&tau;})}\mathrm{Bd\&tau;}\&CenterDot;g_k---\left(7\right)$

By x (t _k) referred to as x _k, each time interval Δ t has following discrete equation form

x _k+1＝A _kx _k+B _ku _k+E _kg _k(8)

In formula, A _k=e ^{a Δ t}; $B_k=e^{A(\mathrm{\&Delta;t}-{\mathrm{\&delta;}}_f)}{\&Integral;}_0^{{\mathrm{\&delta;}}_f}{e}^{A({\mathrm{\&delta;}}_f-\mathrm{\&tau;})}\mathrm{Bd\&tau;};E_k={\&Integral;}_0^{\mathrm{\&Delta;t}}{e}^{A(\mathrm{\&Delta;t}-\mathrm{\&tau;})}\mathrm{Bd\&tau;}.$