CN102890506A - 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 little celestial body of Constraint-based planning approaches section guidance control method

Technical field

The present invention relates to a kind of guidance control method to the small celestial body exploration device, particularly a kind of guidance control method that near little celestial body process, satisfies state and control constraint.

Background technology

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

Prior art is referring to R.R.Sostaric, J.R.Rea.Powered descent guidance methods for the moon and mars.San Francisco, USA:American Institute of Aeronautics and Astronautics Inc, 2005., thereby traditional method for planning track is to utilize polynomial expression that detector current state and expectation state are carried out match to form the track path point, with it as the analytical form of approaching system optimization solution.Because this method is simple and operand is low, therefore become in the past detector near the selection of task section guidance control strategy.Yet the solution of nonlinear detector dynamical system is difficult to replace with simple polynomial form, is not the optimum solution of system by the formed path of traditional method for planning track point; Simultaneously, along with the complicated of detection mission and improving constantly that detection requires, the tradition method for planning track can't satisfy state constraint strict in the task process and control constraint, therefore seeking the track guidance control method of optimizing becomes the problem of recent domestic focus of attention, referring to S.R.Ploen, A.B.Acikmese, A.Wolf.A comparison of powered descent guidance laws for Mars pinpoint landing.Reston, VA, USA:American Institute of Aeronautics and Astronautics Inc, 2006 pairs of multiple planet method of guidances have carried out detailed comparative analysis.

Summary of the invention

The present invention is directed to the defective that little celestial body approaches section task tradition fitting of a polynomial guidance control algolithm, Dynamic Constraints, state constraint and control constraint are introduced in the process of trajectory planning simultaneously, satisfy the feasible solution of the Kind of Nonlinear Dynamical System with state and control constraint by finding the solution the convex programming problem generation, guide a control performance thereby improve the approaching section of little celestial body.

The little celestial body of this kind Constraint-based planning approaches section guidance control method, specifically may further comprise the steps:

The first step: set up detector near little astrodynamics model;

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

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

The 4th step: begun by this initial reference track, introduce state constraint and control constraint, find the solution the dynamics path planning with control constraint and track state constraint by iterating, obtain near guided paths feasible in the target celestial body process.

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

In the formula, r ∈ R ³Be the radius vector of detector under little celestial body barycenter is connected coordinate system; ω ∈ R ³Be the celestial body angular spin rate; U ∈ R ³Control acceleration for effect; D is disturbing acceleration; G is gravitational acceleration; Suppose that little celestial body angular spin rate is stable and be a normal value, namely

Then the kinetic model in (1) formula is by following state space description form

In the 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] ^TC _r=[0 I] and

Matrix representation for vectorial multiplication cross ω * ().

Wherein discretize adopts following methods in the second step:

At Fixed Time Interval Δ t, wherein Δ t 〉=δ _f+ δ _sIn, control acceleration input u (t) is 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 the formula, t _K+1-t _k=Δ t and u _kBe limited control inputs, in the thruster opening process, the solution of formula (3) institute descriptive system is at t=t _k+ δ _fConstantly be

$x(t_k+{\mathrm{\δ}}_f)=e^{{\mathrm{A\δ}}_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\τ}\·g_k---\left(5\right)$

In the formula, g _kFor the celestial body gravitation potential function at state x _kOn gradient, namely with the numerical value the Gradient of Gravitation as gravitational acceleration input item in the discrete model; In thruster closing process subsequently, the solution of formula (3) institute descriptive system is at t=t _K+1Constantly be

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

With formula (5) substitution 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\τ}\·g_k---\left(7\right)$

With x (t _k) brief note is for x _k, at each time interval Δ t following discrete equation form is arranged

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

In the 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 match obtains the initial reference track in the 4th step

To determine

K=0, L, N-1; This match track satisfies at initial time t=0 (k=0)

x _SBe current state, terminal constantly t=T, k=N satisfies

x _FBe terminal 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, expectation SOT state of termination x _FAnd initial reference track

Carry out j=0, L, M-1 goes on foot iteration:

(1) guided paths of utilizing the j time iteration to produce

Calculate and upgrade discrete time-varying model parameter (k=0, L, N-1);

(2) find the solution following SOCP problem 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{subjectto}\\ 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,L,M-1)\end{array}\right\}---\left(9\right)$

In the formula, α in the cost function, β are respectively the weight of burnup and energy consumption, and (α, β)=(1,0) be used for to optimize burnup, and (α, β)=(0,1) is used for optimizing energy consumption; Relax terminal constraint, the terminal velocity constraint is added cost function, get ε＞0, E _r=[I 0], E _v=[0 I]; State constraint territory X and controlled quentity controlled variable constrained domain U are the convex set of utilizing second order circular cone constraint specification; Get constrained parameters 0＜κ＜1, to guarantee being a Cauchy sequence at each status switch constantly, so far resulting convex programming problem can utilize interior point method to carry out Efficient Solution.

Principle of work of the present invention: the target celestial body gravity model is introduced in the trajectory planning process, add simultaneously near in the process to the constraint of state and control, it is that the second order circular cone of performance index is planned (SOCP) problem that the nonlinear kinetics path planning problem that will have control constraint and track state constraint is converted into to optimize burnup.And then this convex programming problem can carry out Efficient Solution to it by interior point method.The path point of planning generation is not only traditional target following state thus, and this state is to consider that the target celestial body gravitational field affects satisfied feasible solution with state and control constrained dynamics system in the situation; In addition, on-board equipment and the required thruster shut-in time requirement of filtering estimation are satisfied in the input of gained feedforward control.

Beneficial effect of the present invention:

The convex programming guidance algorithm has provided one and has satisfied simultaneously state constraint and the nominal trajectory of controlling constraint.Compare with traditional fitting of a polynomial guidance algorithm, because Dynamic Constraints, state constraint and control constraint are introduced in the process of trajectory planning simultaneously, therefore the path point that produces based on the guidance algorithm of convex programming is to satisfy the feasible solution with state and Kind of Nonlinear Dynamical System of control constraint, thereby has improved the guidance control performance.Simultaneously, because the high efficiency that convex programming problem is found the solution, this guidance control method has the potentiality of online planning, can satisfy small-sized spacecraft and cheap spacecraft because computational resource is limited and need online situation about processing in real time.

Embodiment

For making the purpose, technical solutions and advantages of the present invention clearer; the below elaborates to embodiments of the invention: present embodiment is implemented under take technical scheme of the present invention as prerequisite; provided detailed embodiment and concrete operating process, but protection scope of the present invention is not limited to following embodiment.

Under being connected coordinate system, little celestial body barycenter sets up detector near little astrodynamics model, and can be by following the Representation Equation

In the formula, r ∈ R ³Be the radius vector of detector under little celestial body barycenter is connected coordinate system; ω ∈ R ³Be the celestial body angular spin rate; U ∈ R ³Control acceleration for effect; D is disturbing acceleration (mainly launching generation by extraneous gravitation interference, sun optical pressure, comet dust); G is gravitational acceleration.

Suppose that little celestial body angular spin rate is stable and be that one often value is (namely

), then the kinetic model in (1) formula can have following state space description form

In the 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] ^TC _r=[0 I] and

Matrix representation for vectorial multiplication cross ω * ().It should be noted that here not to g (C _rX) carry out linearization process, in subsequent design, will utilize the numerical value the Gradient of Gravitation to replace; And will being used as the uncertain factor of nominal system dynamics model (3), processes by FEEDBACK CONTROL disturbing acceleration d.

Consider the limited opening time δ of thruster _fAnd shut-in time δ _sRequirement, system (3) is carried out discretize, the 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) can 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 the formula, t _K+1-t _k=Δ t and u _kBe limited control inputs.In the thruster opening process, the solution of formula (3) institute descriptive system is at t=t _k+ δ _fConstantly can be approximately

$x(t_k+{\mathrm{\δ}}_f)=e^{{\mathrm{A\δ}}_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\τ}\·g_k---\left(5\right)$

In the formula, g _kFor the celestial body gravitation potential function at state x _kOn gradient (namely with the numerical value the Gradient of Gravitation as gravitational acceleration input item in the discrete model).In thruster closing process subsequently, the solution of formula (3) institute descriptive system is at t=t _K+1Constantly can be approximately

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

With formula (5) substitution formula (6), can get

$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\τ}\·g_k---\left(7\right)$

With x (t _k) brief note is for x _k, at each time interval Δ t following discrete equation form is arranged

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

In the 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\τ}.$

In the discrete system model that formula (8) is described, because g _kDepend on guided paths x _k, therefore can obtain an initial reference track by linearity or fitting of a polynomial first

To determine

(k=0, L, N-1).This match track satisfies at initial time t=0 (k=0)

(x _SBe current state), terminal constantly t=T (k=N) satisfies

(x _FBe terminal expectation state), begun by this initial guess again, introduce state constraint and control constraint, find the solution the dynamics path planning problem with control constraint and track state constraint by iterating, with searching near guided paths { u feasible in the 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), expectation SOT state of termination x _FAnd initial reference track

Carry out j=0, L, M-1 goes on foot iteration:

(1) guided paths of utilizing the j time iteration to produce

Calculate and upgrade discrete time-varying model parameter

(k=0, L, N-1);

(2) find the solution following SOCP problem 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{subjectto}\\ 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,L,M-1)\end{array}\right\}---\left(9\right)$

In the formula, α in the cost function, β are respectively the weight of burnup and energy consumption, and (α, β)=(1,0) be used for to optimize burnup, and (α, β)=(0,1) is used for optimizing energy consumption; Calculate for the ease of processing, relax terminal constraint, the terminal velocity constraint is added cost function, get ε＞0, E _r=[I 0], E _v=[0 I]; State constraint territory X and controlled quentity controlled variable constrained domain U are the convex set of utilizing second order circular cone constraint specification; For guaranteeing that the convex programming guidance algorithm produces the convergence of track, get constrained parameters 0＜κ＜1, to guarantee being a Cauchy sequence at each status switch constantly.So far resulting convex programming problem can utilize interior point method to carry out Efficient Solution.

Approach the section task as example take asteroid Eros, the detector model parameter is as shown in table 1.Target celestial body quadravalence gravitational potential function model is suc as formula shown in (10).Wherein, μ _A=GM is the product of universal gravitational constant and little day weight; A is the nominal radius of little celestial body; R, θ,

Being respectively that little celestial body barycenter is connected is detector radial distance, right ascension and declination under the coordinate system; C _MnBe humorous coefficient of ball.

Table 1 detector simulation parameters

Target celestial body Proximity operation guidance is in expectation time kept in reserve T=300s with the control target, make detector under control strategy u+ δ u effect by current state x _S=[8,950 20-50 1.5 2 0] ^TMove to expectation end state x along path planning _F=[8,450 0000 0] ^T, and in this process, satisfy following state constraint and control constraint:

Near the I hypothesis approximated position the about 8445m of catalog of celestial bodies identity distance celestial body barycenter bumps against for avoiding detector and celestial body, sets nominal state constraint set X={x|c ^TX 〉=1, c=[1/8445 0 _{1 * 5}] ^T;

Be limited on the II control inputs amplitude $U_{\max }=\frac{125}{m}m/s^2;$

The III thruster opening time is δ _f=10s, the thruster shut-in time is δ _s=20s is because guidance control time interval of delta t 〉=δ _f+ δ _s, therefore get Δ t=30s.

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

Claims (4)

1. the little celestial body of a Constraint-based planning approaches section guidance control method, it is characterized in that, specifically may further comprise the steps:

The first step: set up detector near little astrodynamics model;

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

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

The 4th step: begun by this initial reference track, introduce state constraint and control constraint, find the solution the dynamics path planning with control constraint and track state constraint by iterating, obtain near guided paths feasible in the target celestial body process.

2. the little celestial body of a kind of Constraint-based planning as claimed in claim 1 approaches section guidance control method, it is characterized in that, wherein first step medium power is learned model by following the Representation Equation

In the formula, r ∈ R ³Be the radius vector of detector under little celestial body barycenter is connected coordinate system; ω ∈ R ³Be the celestial body angular spin rate; U ∈ R ³Control acceleration for effect; D is disturbing acceleration; G is gravitational acceleration; Suppose that little celestial body angular spin rate is stable and be a normal value, namely

Then the kinetic model in (1) formula is by following state space description form

In the 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] ^TC _r=[0 I] and

Matrix representation for vectorial multiplication cross ω * ().

3. the little celestial body of a kind of Constraint-based planning as claimed in claim 1 approaches section guidance control method, it is characterized in that wherein discretize adopts following methods in the second step:

At Fixed Time Interval Δ t, wherein Δ t 〉=δ _f+ δ _sIn, control acceleration input u (t) is 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 the formula, t _K+1-t _k=Δ t and u _kBe limited control inputs, in the thruster opening process, the solution of formula (3) institute descriptive system is at t=t _k+ δ _fConstantly be

$x(t_k+{\mathrm{\δ}}_f)=e^{{\mathrm{A\δ}}_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\τ}\·g_k---\left(5\right)$

In the formula, g _kFor the celestial body gravitation potential function at state x _kOn gradient, namely with the numerical value the Gradient of Gravitation as gravitational acceleration input item in the discrete model; In thruster closing process subsequently, the solution of formula (3) institute descriptive system is at t=t _K+1Constantly be

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

With formula (5) substitution 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\τ}\·g_k---\left(7\right)$

With x (t _k) brief note is for x _k, at each time interval Δ t following discrete equation form is arranged

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

In the 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\τ}.$

4. the little celestial body such as claim 1 or 2 or 3 described a kind of Constraint-based planning approaches section guidance control method, it is characterized in that, wherein match obtains the initial reference track in the 4th step

To determine

K=0, L, N-1; This match track satisfies at initial time t=0 (k=0) x _SBe current state, terminal constantly t=T, k=N satisfies x _FBe terminal 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, expectation SOT state of termination x _FAnd initial reference track

Carry out j=0, L, M-1 goes on foot iteration:

(1) guided paths of utilizing the j time iteration to produce Calculate and upgrade discrete time-varying model parameter (k=0, L, N-1);

(2) find the solution following SOCP problem 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{subjectto}\\ 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,L,M-1)\end{array}\right\}---\left(9\right)$

In the formula, α in the cost function, β are respectively the weight of burnup and energy consumption, and (α, β)=(1,0) be used for to optimize burnup, and (α, β)=(0,1) is used for optimizing energy consumption; Relax terminal constraint, the terminal velocity constraint is added cost function, get ε＞0, E _r=[I 0], E _v=[0 I]; State constraint territory X and controlled quentity controlled variable constrained domain U are the convex set of utilizing second order circular cone constraint specification; Get constrained parameters 0＜κ＜1, to guarantee being a Cauchy sequence at each status switch constantly, so far resulting convex programming problem can utilize interior point method to carry out Efficient Solution.