CN103995539A - Unmanned aerial vehicle autonomous formation evaluation index and MPC formation control method - Google Patents

Abstract

The invention discloses an unmanned aerial vehicle autonomous formation evaluation index and model prediction control (MPC) formation control method and belongs to the technical field of flight control. According to the method, a definition relevant to unmanned aerial vehicle autonomous formation, a computing method of the collision prevention probability of unmanned aerial vehicle autonomous formation and an unmanned aerial vehicle autonomous formation evaluation index are provided, and the MPC formation control method applicable to unmanned aerial vehicles is provided. The method solves the problem how to give out the collision prevention quantitative index which should be followed by design of an autonomous formation flight control system when it is known that the unmanned aerial vehicle autonomous formation scale is M, the average expectation distance between unmanned aerial vehicles in formation is dr, and the collision probability of the whole formation is required to be smaller than Pc. The method provides a quantitative index basis for engineering design of the unmanned aerial vehicle autonomous formation control system, and the collision problem the unmanned aerial vehicles in dense formation is solved.

Description

A kind of unmanned plane autonomous formation evaluation index and MPC formation control method

Technical field

The invention belongs to flight control technology field, relate to a kind of foundation of unmanned plane autonomous formation evaluation index and the method for designing of MPC (Model Predictive Control, Model Predictive Control) formation control device.

Background technology

Carry out intensive autonomous formation flight at multiple UAVs, distance between unmanned plane and unmanned plane approaches minimum safe distance, and for a certain specific unmanned plane in forming into columns, because it is distributed with other unmanned plane around, lack free movement space, this make the collision prevention of this unmanned plane motor-driven (for example evading the unexpected motion that contiguous unmanned plane produces due to various random disturbance or some emergency case or fault) than more difficult and complicated in loose formation situation.How to judge the dense degree of autonomous formation, and the performance that how to design and evaluate formation control device is that the Control System Design of unmanned plane autonomous formation needs one of problem that emphasis considers.

Summary of the invention

The present invention, in order to solve problems of the prior art, provides a kind of unmanned plane autonomous formation evaluation index, and a kind of MPC formation control method is provided simultaneously.

The invention provides a kind of unmanned plane autonomous formation evaluation index and MPC formation control method, wherein, the unmanned plane autonomous formation evaluation index providing is:

$\mathrm{\σ}\left(t\right)=\sqrt{E[\mathrm{\Δ}{d}_f\left(t\right){]}^2-\{E\left[\mathrm{\Δ}{d}_f\left(t\right)\right]{\}}^2}=\sqrt{E[\mathrm{\Δ}{d}_f\left(t\right){]}^2-[\mathrm{\Δ}{d}_{\mathrm{Ef}}\left(t\right){]}^2}$

Wherein, Δ d _f(t) represent the average connection safe distance surplus that the t moment forms into columns; σ (t) is t moment unmanned plane autonomous formation evaluation index, is Δ d _f(t) mean square deviation.

Described MPC formation control method comprises the steps:

The first step, sets up unmanned plane formation motion model;

Second step, realizes MPC formation control, comprising: forecast model, rolling optimization and feedback compensation.

The described first step is set up unmanned plane formation motion model, specifically:

If two unmanned plane W and L, speed is respectively V _wand V _l, flight path drift angle is respectively with unmanned plane W is d to the distance of unmanned plane L, x ₀, y ₀for apart from d at unmanned plane W flight path axis system x, the component on y axle; and V _wcbe respectively the instruction of flight path drift angle and the speed command of unmanned plane W, and V _lcbe respectively the instruction of flight path drift angle and the speed command of unmanned plane L, be respectively the time constant of unmanned plane L and the response of unmanned plane W flight path drift angle, the time constant that is respectively unmanned plane L and unmanned plane W speed responsive, the motion model of unmanned plane L and unmanned plane W is as follows:

Forecast model in described second step, specifically:

Utilize the unmanned plane formation motion model of the first step as forecast model, get state controlled being input as the flight path drift angle of node L and speed are regarded as can measurements interference output wherein Δ V=V _l– V _w, while ignoring D, the discrete form of unmanned plane formation motion model is:

$\left\{\begin{array}{c}X(k+1)=\mathrm{AX}\left(k\right)+\mathrm{Bu}\left(k\right)\\ Y\left(k\right)=\mathrm{CX}\left(k\right)\end{array}\right.$

Wherein, X (k) represents the state of moment k, and u (k) represents the controlled input of moment k, and Y (k) represents the output of moment k, and A, B and C represent matrix of coefficients; If u (k+i|k) is illustrated in the controlled input value u in the following k+i moment of moment k supposition; X (k+i|k), Y (k+i|k) is illustrated respectively in moment k, utilizes u (k+j|k), (j=0,1 ..., i-1) to X, the predicted value that Y makes; The predictive equation of the state X obtaining is:

Wherein, P is prediction time domain, and N is for controlling time domain, and u (k-1) represents the controlled input in k-1 moment, difference DELTA u (k+j|k)=u (k+j|k)-u (k+j-1|k);

The predictive equation of output Y is:

Y(k+j|k)＝CX(k+j|k),j＝1,2,...,P

Rolling optimization in described second step, specifically:

If formation flight at the cost function J of moment k is:

$J=\underset{l=1}{\overset{P}{\mathrm{\Σ}}}\left(\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\right||{\mathrm{\ω}}_y(l,m)\·(Y_m(k+l|k)-{\mathrm{ref}}_m(k+l))|{|}^2+\underset{m=1}{\overset{N}{\mathrm{\Σ}}}\left|\right|{\mathrm{\ω}}_{\mathrm{\Δu}}(l,m)\·\mathrm{\Δ}{u}_m(k+l|k)\left|{|}^2\right)$

Wherein, ω _yfor output punishment weight matrix, ω _y(l, m) represents the output punishment weight that time domain is corresponding to (l, m); Y _m(k+l|k) be under control time m, the predicted value in k+l moment of calculating in k moment; Ref is spacing hold instruction, ref _m(k+l) the spacing hold instruction in the k+l moment of be under control time m, the k moment calculating; ω _{Δ u}for input slew rate punishment weight matrix, ω _{Δ u}(l, m) represents the input slew rate punishment weight that time domain is corresponding to (l, m); Δ u _m(k+l|k) be the k+l moment input slew rate calculating in the k moment under control time m; Make each time domain equate the weight of (l, m), ω _yand ω _{Δ u}deteriorate to respectively a row vector;

Constrained optimization problem below the controlled quentity controlled variable u of moment k (k) is:

Δu _opt＝arg?minJ

Wherein, with be respectively flight path drift angle and speed command restriction; Δ u _optfor first controlled quentity controlled variable of optimal sequence, u ₁, u ₂be respectively and V _c; I _{2 × 2}represent 2 × 2 unit matrix, O _{(N-1) × (N-1)}represent 0 matrix of (N-1) × (N-1);

Feedback compensation in described second step, comprising: (1) is the original state using X (k) as system directly; (2) the flight path drift angle of unmanned plane L and speed being regarded as can measurements interference by network delivery to unmanned plane W and append to the predicted value of forecast model.

Advantage of the present invention and good effect are:

(1) the autonomous formation evaluation index proposing in the present invention is for evaluating the density of formation, and the present invention has put forward the collision prevention quantizating index that autonomous formation Design of Flight Control should be followed;

(2) the MPC formation control method proposing in the present invention has solved the collision prevention problem of unmanned plane close/intra;

(3) this collision prevention method for designing for autonomous formation control system proposed by the invention, for the practical implementation of unmanned plane autonomous formation Control System Design provides the index reference quantizing.

Brief description of the drawings

Fig. 1 is the safe distance schematic diagram of unmanned plane member in the present invention;

Fig. 2 is the possible position graph of a relation of unmanned plane i and neighbour's unmanned plane j in the present invention;

Fig. 3 is the internodal relative motion relation of forming into columns;

Fig. 4 is Model Predictive Control block diagram;

Fig. 5 is Model Predictive Control principle;

Fig. 6 is the simulation result of PID formation control device in the present invention, and a and b are respectively unmanned plane relative distance d at certain unmanned aerial vehicle flight path coordinate system x, the component on y axle, and c is the schematic diagram of unmanned plane relative distance d;

Fig. 7 is the simulation result figure that adopts MPC formation control method of the present invention, and A and B are respectively unmanned plane relative distance d at certain unmanned aerial vehicle flight path coordinate system x, the component on y axle, and C is the schematic diagram of unmanned plane relative distance d.

Embodiment

Below in conjunction with drawings and Examples, the present invention is described in detail.

Model of the present invention unmanned plane autonomous formation evaluation index.

The probability bumping when hypotelorism between the unmanned plane of tight formation flight is higher, and unmanned plane all can be determined a safe distance d according to the actual conditions of self conventionally _s.D _sbe so a kind of distance, in the time that the distance of certain unmanned plane and another unmanned plane or threat is less than this distance, this unmanned plane must be taked corresponding measure.Distance when between two unmanned planes is greater than d _stime, the collision opportunity between unmanned plane is less, at this moment forms into columns that it is loose to think, i.e. less demanding to tail clearance retentive control device; Distance when between unmanned plane approaches d _stime, unmanned plane needs to evade possible collision threat at any time, and at this moment, to the having relatively high expectations of tail clearance retentive control device, formation can be thought closely.

The present invention will address the problem: be M in the scale that requires autonomous formation to have, the average expectation spacing between the unmanned plane in formation is d _r, the collision probability of whole formation is less than P _cin situation, how to provide the collision prevention quantizating index that autonomous formation Design of Flight Control should be followed.

Utilize topography and geomorphology shielding, in long and narrow escape way, carry out penetration cruise with flat close/intra form, saturation attack mode with highly dense intensity surprise attack formula is the main patterns of warfare that unmanned plane is formed into columns, now in this flat close/intra between unmanned plane, unmanned plane and obstacle or threaten between the probability that bumps increase, pursuing under high penetraton probability principle, for forming into columns, collaborative guidance controller performance has proposed very high requirement.Flat close/intra form is equivalent to whole formation does not almost have the motor-driven free space of collision prevention in the vertical direction, it is motor-driven that whole formation can only be carried out collision prevention in the horizontal direction, this be equivalent to be longitudinally unlimited intensive, there is no a motor-driven free space of collision prevention, be can not by adjusting the method for flying height, to carry out collision prevention motor-driven, therefore in two-dimensional level face, the collision prevention problem of close/intra be the most complicated Basic Problems.To expand to easily three-dimensional formation problem to close/intra way to solve the problem and conclusion in two-dimensional level face.

The density degree (abbreviation density) that unmanned plane is formed into columns be one describe form into columns in the key concept of distance size characteristic between unmanned plane member, it has material impact for the collaborative guidance of unmanned plane autonomous formation Control System Design, will the density of unmanned plane autonomous formation be analyzed below.

In order to describe the density of unmanned plane autonomous formation, given first form into columns in unmanned plane member safe distance concept and definition.Conventionally be suitable for unmanned plane formation at some and carry out in saturation attack task, need to utilize topography and geomorphology shielding, to form into columns form penetration cruise in narrow escape way, geography, the weather environment of formation flight are complicated and changeable.If the hypotelorism in forming into columns between unmanned plane, between unmanned plane and obstacle, the probability bumping so will increase, thereby the flight safety that threatens whole unmanned plane to form into columns reduces the fighting efficiency that unmanned plane is formed into columns.

How to be definite suitable safe distance between the unmanned plane in formation flight, be one of the basic link of collaborative guidance Control System Design of forming into columns.Conventionally determine a safe distance d according to unmanned plane actual conditions separately _si, and utilize it to control the distance between unmanned plane member in formation (be called for short machine spacing) d (t) _ij(t).Here d, _ij(t)=d _ji(t), and i ≠ j, the distance between unmanned plane i and unmanned plane j represented, d _si(t) safe distance of expression unmanned plane i.D _ij(t) number wherein i, j=1,2 ..., M (t), is the numbering of unmanned plane member in forming into columns, i namely unmanned plane i of unmanned plane member, and M (t) is current t moment unmanned plane sum.

Definition 1: unmanned plane member's safe distance d _si(t) be so a kind of distance, as the machine spacing d between other any unmanned plane j in i frame unmanned plane and formation _ij(t) when being less than it, i.e. d _ij(t) < d _si(t), i frame unmanned plane must be taked corresponding collision prevention measure; Work as d _ij(t)=d _si(t), time, i frame unmanned plane is in taking the standby condition of corresponding collision prevention measure.

Definition 2: unmanned plane member's safe distance surplus Δ d _ij(t), refer to the machine spacing d of unmanned plane member in formation _ij(t) with its safe distance d _si(t) poor.Be Δ d _ij(t)=d _ij(t)-d _si(t).

As shown in Figure 1, unmanned plane member's safe distance d in impact formation _si(t) factor mainly contains three: sensor measurement error component, mobility component and the network is guided component.Unmanned plane member's sensor measurement error component is mainly by the sensor technology Determines of unmanned plane self, if the distance between two unmanned planes is less than sensor measurement error component, thinks that collision has occurred for they.Mobility component has reflected that unmanned plane member has the ability of evading the collision threat that happens suddenly.In the time of single frame unmanned plane solo hop, generally in the time that safe distance is less than sensor measurement error component and mobility component sum, unmanned plane will be taked collision prevention measure.The network is guided component be due to formation supporting network to information transmit the factor such as time delay, packet loss to form into columns in unmanned plane member safe distance newly-increased restriction.

Except above-mentioned employing machine spacing d _ij(t) with safe distance d _si(t) the safe distance surplus Δ d of magnitude relationship _ij(t) describe outside the density of formation, also should investigate unmanned plane member free space situation around in formation, investigate the dense degree (abbreviation closeness) of forming into columns.Because the machine spacing d between i frame unmanned plane member and other unmanned planes _ij(t) be less than its safe distance d _si(t), i.e. safe distance surplus Δ d _ij(t) when < 0, it is motor-driven that this unmanned plane need to carry out collision prevention for the collision threat facing, and this evasion manoeuvre is usually towards its free space around.If when now i frame unmanned plane does not around have free space again, so it this how to carry out collision threat evasion manoeuvre? as can be seen here, except needs use safe distance d _si(t) concept is weighed outside the density of formation, also needs to introduce a variable of describing unmanned plane member free space around in formation, is called the variable n of connection degree _ci(t), reflecting that formation may have many " close " on earth, how is the closeness of forming into columns?

Connection degree n _ci(t) refer to that on two dimensional surface, distance unmanned plane i approaches d _si(t) number of adjacent unmanned plane.Find by analysis, in two dimensional surface, the i frame unmanned plane member of close/intra around distance approaches d _si(t) mean value of adjacent number of members is about 6 to 7 left and right, the connection degree n of thick in forming into columns _ci(t) be about 6 to 7 left and right.In order to describe the average closeness level of whole formation, introduce average connection degree therefore, according to the average connection degree n forming into columns _c(t), can carry out to judge approx the dense degree of formation.

Collision prevention problem between unmanned plane member is one of the basic problem that the formation flight of highly dense intensity needs emphasis to solve, and it is directly connected to independence that unmanned plane forms into columns and the fighting efficiency of concertedness level and entirety thereof.Provide one below and take the relational expression between the motor-driven probability of collision prevention and density and closeness about unmanned plane member and formation.

Unmanned plane member's safe distance surplus Δ d in formation _ij(t)=d _ij(t)-d _si(t) be ergodic stochastic process, below by with the variable of time correlation in will save the time and write a Chinese character in simplified form, for example Δ d _ij(t) be abbreviated as Δ d _ij, other is also done this and writes a Chinese character in simplified form.According to central limit theorem, suppose under the combined action of various mutual inapparent enchancement factors unmanned plane member's safe distance surplus Δ d _ijthe approximate safe distance surplus Δ d obeying taking formation as its setting _rijfor mathematical expectation, with σ _ij ²for the normal distribution of variance, in forming into columns, unmanned plane member's machine spacing is not more than its safe distance, i.e. Δ d _ij≤ 0 o'clock, unmanned plane member just needed or prepares to take the motor-driven measure of collision prevention, as shown in Figure 2.Therefore, between unmanned plane i and unmanned plane j, need or prepare the Probability p of taking collision prevention motor-driven _b(i, j) is:

p _b(i,j)＝P{Δdi _j≤0}

$\begin{array}{c}={\&Integral;}_{-\&infin;}^{\frac{-\mathrm{\&Delta;}{d}_{\mathrm{rij}}}{{\mathrm{\&sigma;}}_{\mathrm{ij}}}}\frac{1}{\sqrt{2\mathrm{\&pi;}}}\&CenterDot;e^{-\frac{{\mathrm{\&tau;}}^2}{2}\mathrm{d\&tau;}}\\ =\mathrm{\&Phi;}\left(\frac{-\mathrm{\&Delta;}{d}_{\mathrm{rij}}}{{\mathrm{\&sigma;}}_{\mathrm{ij}}}\right)\\ =\mathrm{\&Phi;}(-{\mathrm{\&lambda;}}_{\mathrm{rij}})\end{array}---\left(1\right)$

Between unmanned plane i and unmanned plane j without taking the motor-driven probability q of collision prevention _b(i, j) is:

q _b(i,j)＝P{Δd _ij＞0}

$\begin{array}{c}=1-{\&Integral;}_{-\&infin;}^{\frac{-\mathrm{\&Delta;}{d}_{\mathrm{rij}}}{{\mathrm{\&sigma;}}_{\mathrm{ij}}}}\frac{1}{\sqrt{2\mathrm{\&pi;}}}\&CenterDot;e^{-\frac{{\mathrm{\&tau;}}^2}{2}\mathrm{d\&tau;}}\\ =1-\mathrm{\&Phi;}\left(\frac{-\mathrm{\&Delta;}{d}_{\mathrm{rij}}}{{\mathrm{\&sigma;}}_{\mathrm{ij}}}\right)\\ =1-\mathrm{\&Phi;}(-{\mathrm{\&lambda;}}_{\mathrm{rij}})\end{array}---\left(2\right)$

In formula, safe distance surplus ratio, for normal probability integral.

Above formula is the motor-driven new probability formula of collision prevention of only having considered density factor, because all likely bump danger in any direction of unmanned plane i, so also closeness should be reflected the impact of the motor-driven probability of collision prevention, so the motor-driven Probability p of the collision prevention of unmanned plane i oneself _b(i) be:

$\begin{array}{c}{p}_b\left(i\right)\underset{j=1}{\overset{n_{\mathrm{ci}}}{\mathrm{\&Sigma;}}}{p}_b(i,j)[q_b(i,j){]}^{n_{\mathrm{ci}}-1}+\underset{j=1}{\overset{C_{n_{\mathrm{ci}}}^2}{\mathrm{\&Sigma;}}}[p_b(i,j){]}^2[q_b(i,j){]}^{n_{\mathrm{ci}}-2}+\&CenterDot;\&CenterDot;\&CenterDot;\\ +\underset{j=1}{\overset{C_{n_{\mathrm{ci}}}^k}{\mathrm{\&Sigma;}}}\left[p_b(i,j){]}^k\right[q_b(i,j){]}^{n_{\mathrm{ci}}-k}+\&CenterDot;\&CenterDot;\&CenterDot;+\underset{j=1}{\overset{n_{\mathrm{ci}}}{\mathrm{\&Pi;}}}{p}_b(i,j)\\ =\underset{k=1}{\overset{n_{\mathrm{ci}}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{C_{n_{\mathrm{ci}}}^k}{\mathrm{\&Sigma;}}}\left[p_b(i,j){]}^k\right[q_b(i,j){]}^{n_{\mathrm{ci}}-k}\\ =\underset{k=1}{\overset{n_{\mathrm{ci}}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{C_{n_{\mathrm{ci}}}^k}{\mathrm{\&Sigma;}}}\left[\mathrm{\&Phi;}(-{\mathrm{\&lambda;}}_{\mathrm{rij}}){]}^k\right[1-\mathrm{\&Phi;}(-{\mathrm{\&lambda;}}_{\mathrm{rij}}){]}^{n_{\mathrm{ci}}-k},(n_{\mathrm{ci}}\&GreaterEqual;2)\end{array}---\left(3\right)$

Wherein, n _cifor unmanned plane i neighbour's number around, i.e. degree of connection, for number of combinations.For connection degree n _ci=1 o'clock, the Probability p that collision prevention is motor-driven _b(i) adopt formula (1) to calculate.

Especially, compare λ in the safe distance surplus of setting of forming into columns _rij=0, by the safe distance of unmanned plane member i as under the expectation machine spacer conditions of setting, as connection degree n _ci, utilize formula (3) can calculate the motor-driven Probability p of collision prevention of unmanned plane member i at=2 o'clock _b(i)=0.75; But as connection degree n _ci=6 o'clock, p _b(i)=0.984375.That is to say, if the surplus of setting compares λ _rij=0, so for having compared with for the formation of high density and closeness, unmanned plane member i almost must take the motor-driven measure of collision prevention.

The motor-driven Probability p of collision prevention is taked in whole formation _bfor:

$\begin{array}{c}{p}_b=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{p}_b\left(i\right)[1-p_b\left(i\right){]}^{M-1}+\underset{i=1}{\overset{C_M^2}{\mathrm{\&Sigma;}}}[p_b\left(i\right){]}^2[1-p_b\left(i\right){]}^{M-2}+\&CenterDot;\&CenterDot;\&CenterDot;\\ +\underset{i=1}{\overset{C_M^k}{\mathrm{\&Sigma;}}}\left[p_b\left(i\right){]}^k\right[1-p_b\left(i\right){]}^{M-k}+\&CenterDot;\&CenterDot;\&CenterDot;\underset{i=1}{\overset{M}{\mathrm{\&Pi;}}}{p}_b\left(i\right)\\ =\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{C_M^k}{\mathrm{\&Sigma;}}}\left[p_b\left(i\right){]}^k\right[1-p_b\left(i\right){]}^{M-k}\\ =\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{C_M^k}{\mathrm{\&Sigma;}}}\left[p_b\left(i\right){]}^k\right[q_b\left(i\right){]}^{M-k}\end{array}---\left(4\right)$

In formula, q _b(i)=1-p _b(i) for unmanned plane member i is without taking the motor-driven probability of collision prevention, the nodes that M is whole formation, k is the number of nodes that needs or prepare to take the motor-driven measure of collision prevention in whole formation.

From formula (1)-Shi (4), safe distance surplus compare λ _rijunder stable condition, the motor-driven Probability p of collision prevention is taked in whole formation _balong with connection degree n _ciincrease with the increase of formation scale M.

Especially, work as p _bwhen (i, j)=α (0≤α≤1) is normal value, formula (3) and formula (4) become respectively:

$\begin{array}{c}{p}_b\left(i\right)=p=n_{\mathrm{ci}}\&CenterDot;\mathrm{\&alpha;}{\mathrm{\&beta;}}^{n_{\mathrm{ci}}-1}+C_{n_{\mathrm{ci}}}^2{\mathrm{\&alpha;}}^2\&CenterDot;{\mathrm{\&beta;}}^{n_{\mathrm{ci}}-2}+\&CenterDot;\&CenterDot;\&CenterDot;+C_{n_{\mathrm{ci}}}^k\&CenterDot;{\mathrm{\&alpha;}}^k{\mathrm{\&beta;}}^{n_{\mathrm{ci}}-k}+\&CenterDot;\&CenterDot;\&CenterDot;+{\mathrm{\&alpha;}}^{n_{\mathrm{ci}}}\\ =\underset{k=1}{\overset{n_{\mathrm{ci}}}{\mathrm{\&Sigma;}}}{C}_{n_{\mathrm{ci}}}^k\&CenterDot;{\mathrm{\&alpha;}}^k{\mathrm{\&beta;}}^{n_{\mathrm{vi}}-k}\\ =1-{\mathrm{\&beta;}}^{n_{\mathrm{ci}}}\end{array}---\left(5\right)$

$\begin{array}{c}{P}_b=M\&CenterDot;p{q}^{M-1}+C_M^2\&CenterDot;p^2{q}^{M-2}+\&CenterDot;\&CenterDot;\&CenterDot;+C_M^k\&CenterDot;p^k{q}^{M-k}+\&CenterDot;\&CenterDot;\&CenterDot;+p^M\\ =\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_M^k\&CenterDot;p^k{q}^{M-k}\\ =1-q^M\end{array}---\left(6\right)$

Wherein, β=1-α, q=1-p.

If by unmanned plane member's safe distance d _sibe set as its sensor measurement error component, the safe distance surplus that order is set compares λ _rij=0, the probability that what calculated formula (1)-(6) so is exactly probability that corresponding unmanned plane member collides with each other and formation bump.

Just utilize the concept that formation dynamic stability characteristic is had to the motor-driven probability of collision prevention of crucial influence below, from safe distance surplus and around free space two aspects unmanned plane autonomous formation is carried out to class definition.According to density and closeness concept, unmanned plane autonomous formation can be divided into three types: loosen and form into columns, closely form into columns and close/intra.

Definition 3: the safe distance surplus that loose formation refers to all unmanned plane members than be all not less than 3 or all between 2 and 3 and on average degree of connection be not more than 3 unmanned plane formation.Be λ _rij>=3 or 2≤λ _rij< 3 & n _c≤ 3, that is to say, all unmanned plane members' safe distance surplus be all not less than its mean square deviation of 3 times or all between the mean square deviation of 2 times and 3 times and on average degree of connection be not more than 3, i.e. Δ d _rij>=3 σ _ijor 2 σ _ij≤ Δ d _rij< 3 σ _ijaMP.AMp.Amp n _c≤ 3.In the time that unmanned plane is formed into columns the loose formation flight of employing, in formation, unmanned plane member need to take the motor-driven probability of collision prevention lower, between unmanned plane member, there is hardly the chance bumping, the density of forming into columns and closeness are all in very low state, and the performance requirement of the collaborative guidance controller of forming into columns for unmanned plane is in the case not high.

Definition 4: the safe distance surplus that refers to all unmanned plane members of closely forming into columns is than all between between 2 and 3 and the on average unmanned plane formation of degree of connection between 3 and 6.I.e. 2≤λ _rij< 3 and 3 < n _c< 6.That is to say, all unmanned plane members' safe distance surplus all between the mean square deviation of its 2 times and 3 times and on average degree of connection between 3 and 6, i.e. 2 σ _ij≤ Δ d _rij< 3 σ _ijand 3 < n _c< 6.From the density of forming into columns, different from loose formation, in the time adopting tight formation flight, in formation, unmanned plane member need to take the motor-driven probability of collision prevention larger, exists the probability bumping higher between unmanned plane member; And from form into columns closeness, due to form into columns average connection degree be less than 6, these unmanned planes member towards periphery needed free space of evasion manoeuvre exists, these free spaces for unmanned plane member towards periphery evasion manoeuvre be moderate.Compared with loose formation situation, the performance requirement of the collaborative guidance controller of at this moment unmanned plane being formed into columns is higher.

Definition 5: the safe distance surplus that close/intra refers to all unmanned plane members than be all less than 2 or all between 2 and 3 and on average degree of connection be not less than 6 unmanned plane formation.Be λ _rij< 2 or 2≤λ _rij< 3 & n _c>=6.That is to say, all unmanned plane members' safe distance surplus be all less than the mean square deviation of its twice or on average degree of connection be not less than 6, i.e. Δ d _rij< 2 σ _ijor 2 σ _ij≤ Δ d _rij< 3 σ _ijaMP.AMp.Amp n _c>=6.From the density of forming into columns, close/intra and the not too large difference of closely forming into columns, just make machine spacing d _ijfurther close to its safe distance d _si, make unmanned plane member need to take the tight formation of the motor-driven likelihood ratio of collision prevention further to increase.Meanwhile, because the closeness of forming into columns improves, the average connection degree of forming into columns is not less than 6, causes unmanned plane member almost there is no the needed free space of evasion manoeuvre towards periphery, exists the probability bumping very high between unmanned plane member.Therefore,, compared with tight formation situation, the performance requirement of the collaborative guidance controller that tight formation flight is formed into columns to unmanned plane is higher.

Definition 6: the average safe distance d of formation _sf(t), refer to all unmanned plane member safe distance d in formation _si(t) mean value.? $d_{\mathrm{sf}}\left(t\right)=\frac{1}{M\left(t\right)}\underset{i=1}{\overset{M\left(t\right)}{\mathrm{\&Sigma;}}}{d}_{\mathrm{si}}\left(t\right).$

Definition 7: the average junctor spacing d of formation _f(t), refer to the machine spacing d of all unmanned plane members in formation and its contiguous unmanned plane _ij(t) mean value.? $d_f\left(t\right)=\frac{1}{M\left(t\right)}\underset{i=1}{\overset{M\left(t\right)}{\mathrm{\&Sigma;}}}\left\{\frac{1}{n_{\mathrm{ci}}\left(t\right)}\underset{j=1,i\&NotEqual;j}{\overset{n_{\mathrm{ci}}\left(t\right)}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}\left(t\right)\right\}.$

Definition 8: the average connection safe distance surplus Δ d of formation _f(t), refer to the average junctor spacing d of formation _f(t) with average safe distance d _sf(t) poor.That is:

Δd _f(t)=d _f(t)-d _sf(t)

$\begin{array}{c}=\frac{1}{M\left(t\right)}\underset{i=1}{\overset{M\left(t\right)}{\mathrm{\&Sigma;}}}\left\{\frac{1}{n_{\mathrm{ci}}\left(t\right)}\underset{j=1,i\&NotEqual;j}{\overset{n_{\mathrm{ci}}\left(t\right)}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}\left(t\right)\right\}-\frac{1}{M\left(t\right)}\underset{i=1}{\overset{M\left(t\right)}{\mathrm{\&Sigma;}}}{d}_{\mathrm{si}}\left(t\right)\\ =\frac{1}{M\left(t\right)}\underset{i=1}{\overset{M\left(t\right)}{\mathrm{\&Sigma;}}}\left\{\right[\frac{1}{n_{\mathrm{ci}}\left(t\right)}\underset{j=1,i\&NotEqual;1}{\overset{n_{\mathrm{ci}}\left(t\right)}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}\left(t\right)]-d_{\mathrm{si}}\left(t\right)\}\end{array}$

Formula (1)-(6) also can be used as the judgement schematics of the collaborative guidance controller performance of forming into columns.For the formation with same regular formation, as the safe distance d of unmanned plane _siwith the formation setpoint distance d forming into columns _rijall identical, safe distance surplus Δ d _ijone timing, the size of the mean square deviation of available average machine spacing is as evaluating the collaborative guidance controller performance index of forming into columns, i.e. index σ:

$\mathrm{\&sigma;}\left(t\right)=\sqrt{E[\mathrm{\&Delta;}{d}_f\left(t\right){]}^2-\{E\left[\mathrm{\&Delta;}{d}_f\left(t\right)\right]{\}}^2}=\sqrt{E[\mathrm{\&Delta;}{d}_f\left(t\right){]}^2-[\mathrm{\&Delta;}{d}_{\mathrm{Ef}}\left(t\right){]}^2}---\left(7\right)$

Wherein, ex represents to ask for the mean value of x, and the less expression formation of the value performance of σ is better.

The MPC formation control method of unmanned plane provided by the invention, comprises the steps:

The first step, sets up unmanned plane formation motion model.

The flight formation control system of present stage, is all based on lead aircraft-wing plane pattern, adopts wing plane to follow the formation mode of lead aircraft, and lead aircraft can be that true lead aircraft can be also virtual lead aircraft.As Fig. 3, each node represents a unmanned plane, V _wand V _lbe respectively the speed of node W and node L, with be respectively the flight path drift angle of node W and node L, d is the distance of node W to node L, the difference of the flight path drift angle of node W and node L x and y are the Orthogonal Decomposition of d on node W velocity reversal.(in Fig. 3, be numbered N with node W _w) flight path axis system (velocity reversal of node W is x axle, with the direction of speed vertical-right be y axle) for relative coordinate system, ground coordinate is fixed coordinate system, utilize the absolute velocity=relative velocity+convected velocity relation in theoretical mechanics, the kinematical equation of setting up node L and node W is as follows:

Wherein, represent respectively the velocity vector of node L, node W, represent the distance vector of node W to node L.

Above formula is decomposed in the flight path axis system of node W:

Wherein x, y is that the distance d of the relative node L of node W is at node W flight path axis system x, the component on y axle.Linearization obtains to above formula (9) to adopt low-angle and microvariations hypothesis:

X ₀, y ₀for distance d under a certain disturbance is at unmanned plane W flight path axis system x, the component on y axle.

The motion model of node L and node W:

Wherein, and V _wcbe respectively the instruction of flight path drift angle and the speed command of node W, and V _lcbe respectively the instruction of flight path drift angle and the speed command of node L, and V _lc, V _wcfor flight path drift angle and the speed controlled quentity controlled variable of outer shroud, act on interior ring flight control system as instruction, be respectively the time constant of node L and the response of node W flight path drift angle, be respectively the time constant of node L and node W speed responsive.

Second step, realizes MPC formation control, and design is based on MPC formation control device.

For formation flight control, can be found out by formula (11), movement mechanism between two unmanned planes is clear, model simple, therefore be well suited for application model PREDICTIVE CONTROL, and along with the fast development of robot calculator software and hardware technology, calculate and carrying cost more and more lower, MPC method is applied in the such rapid system of flight control system more and more widely, and therefore the present invention will design the formation control device based on state space MPC according to the feature of formation flight.

MPC control method based on model has three major parts conventionally: forecast model, rolling optimization and feedback compensation, as shown in Figure 4.When in network inducement delay situation, forecast model obtains as shown in step 2.4.

Step 2.1: forecast model.

Forecast model is only focused on the function of model, and do not focus on the form of model, as long as have the model of forecast function, as state equation, this quasi-tradition model of transport function, and for example step response, this class nonparametric model of impulse response etc. all can be used as forecast model use.In formation control of the present invention, utilize formula (11) as forecast model.

Choose state according to formula (11) wherein x, y is the projection in node W flight path axis system of the relative distance of node L and node W; for the flight path drift angle of node W; V _wfor the speed of node W; Controlled being input as wherein and V _wcbe respectively the instruction of flight path drift angle and the speed command of node W.For alleviating computation burden, here the flight path drift angle of node L and speed being regarded as can measurements interference do not use first order modeling to predict it, directly transmit flight path drift angle and the speed of node L by supporting network; Output wherein while ignoring D (D can be added in forecast model subsequently, is considered as the correction of model), the discrete form of formation motion model is:

$\left.\begin{array}{c}X(k+1)=\mathrm{AX}\left(k\right)+\mathrm{Bu}\left(k\right)\\ Y\left(k\right)=\mathrm{CX}\left(k\right)\end{array}\right\}---\left(12\right)$

A, B and C represent matrix of coefficients, namely k the sampling period of moment k herein, if X (k|k)=X (k) is illustrated in moment k, according to the state value of Y (k) and u (k-1) acquisition, because u (k) does not now also calculate; U (k+i|k) is illustrated in the input value u in the following k+i moment of moment k supposition; X (k+i|k), Y (k+i|k) is illustrated in moment k, utilizes u (k+j|k), (j=0,1 ..., i-1) to X, the predicted value that Y makes; P is prediction time domain; N is for controlling time domain, i.e. u (k+i|k)=u (k+N-1|k) N < i < P-1, can obtain by iterative (12):

$\left.\begin{array}{c}X(k+1|k)=\mathrm{AX}\left(k\right)+\mathrm{Bu}\left(k\right|k)\\ X(k+2|k)=\mathrm{AX}(k+1|k)+\mathrm{Bu}(k+1|k)\\ =A^{2}X\left(k\right)+\mathrm{ABu}\left(k\right|k)+\mathrm{Bu}(k+1|k)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ X(k+P|k)=A^{P}X\left(k\right)+A^{P-1}\mathrm{Bu}\left(k\right|k)+...+\mathrm{Bu}(k+P-1|k)\end{array}\right\}---\left(13\right)$

That is:

$X(k+j|k)=A^{j}X\left(k\right)+\left[\begin{array}{cccc}{A}^{j-1}& A^{j-2}& ...& I\end{array}\right]B\left[\begin{array}{c}u\left(k\right|k)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ u(k+j-1|k)\end{array}\right]j=\mathrm{1,2},...,P---\left(14\right)$

Because u (k) is unknown in the k moment, and u (k-1) is known, establishes difference

Δu(k+j|k)＝u(k+j|k)-u(k+j-1|k) (15)

Have:

$\left.\begin{array}{c}u\left(k\right|k)=\mathrm{\&Delta;u}\left(k\right|k)+u(k-1)\\ u(k+1|k)=\mathrm{\&Delta;u}(k+1|k)+\mathrm{\&Delta;u}\left(k\right|k)+u(k-1)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ u(k+N-1|k)=\mathrm{\&Delta;u}(k+N-1|k)+\&CenterDot;\&CenterDot;\&CenterDot;\mathrm{\&Delta;u}\left(k\right|k)+u(k-1)\end{array}\right\}---\left(16\right)$

Have for 1≤j≤N:

$\left.\begin{array}{c}X(k+1|k)=\mathrm{AX}\left(k\right)+B[\mathrm{\&Delta;u}\left(k\right|k)+u(k-1)]\\ X(k+2|k)=A^{2}X\left(k\right)+\mathrm{AB}[\mathrm{\&Delta;u}\left(k\right|k)+u(k-1|k)]\\ +B[\mathrm{\&Delta;u}(k+1|k)+\mathrm{\&Delta;u}\left(k\right|k)+u(k-1)]\\ =A^{2}X\left(k\right)+(A+I)\mathrm{B\&Delta;u}\left(k\right|k)+\mathrm{B\&Delta;u}(k+1|k)\\ +(A+I)\mathrm{Bu}(k-1)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ X(k+N|k)=A^{N}X\left(k\right)+(A^{N-1}+\&CenterDot;\&CenterDot;\&CenterDot;+A+I)\mathrm{Bu\&Delta;}\left(k\right|k)+...\\ +\mathrm{B\&Delta;u}(k+N-1|k)+(A^{N-1}+\&CenterDot;\&CenterDot;\&CenterDot;+A+I)\mathrm{Bu}(k-1)\end{array}\right\}---\left(17\right)$

In formula, I is unit matrix, further obtains:

$X(k+j|k)=A^{j}X\left(k\right)+\left[\begin{array}{ccc}{\mathrm{\&Sigma;}}_{i=0}^{j-1}{A}^{i}B& \&CenterDot;\&CenterDot;\&CenterDot;& B\end{array}\right]\left[\begin{array}{c}\mathrm{\&Delta;u}\left(k\right|k)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \mathrm{\&Delta;u}(k+j-1|k)\end{array}\right]+{\mathrm{\&Sigma;}}_{i=0}^{j-1}{A}^i\mathrm{Bu}(k-1)---\left(18\right)$

Have for N≤j≤P:

$\left.\begin{array}{c}X(k+N1|k)=A^{N+1}X\left(k\right)+(A^N+\&CenterDot;\&CenterDot;\&CenterDot;+A+I)\mathrm{Bu\&Delta;}\left(k\right|k)+...\\ +(A+I)\mathrm{B\&Delta;u}(k+N-1|k)\\ +(A^N+\&CenterDot;\&CenterDot;\&CenterDot;+A+I)\mathrm{Bu}(k-1)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ X(k+P|k)=A^{P}X\left(k\right)+(A^{P-1}+\&CenterDot;\&CenterDot;\&CenterDot;A+I)\mathrm{Bu\&Delta;u}\left(k\right|k)+...\\ +(A^{P-N}+\&CenterDot;\&CenterDot;\&CenterDot;+A+I)\mathrm{B\&Delta;u}(k+N-1|k)\\ +(A^{P-1}+\&CenterDot;\&CenterDot;\&CenterDot;+A+I)\mathrm{Bu}(k+1)\end{array}\right\}---\left(19\right)$

That is:

$\begin{array}{c}X(k+j|k)=A^{j}X\left(k\right)+\left[\begin{array}{ccc}{\mathrm{\&Sigma;}}_{i=0}^{j-1}{A}^{i}B& \&CenterDot;\&CenterDot;\&CenterDot;& {\mathrm{\&Sigma;}}_{i=0}^{j-N}\end{array},A^{i}B\right]\left[\begin{array}{c}\mathrm{\&Delta;u}\left(k\right|k)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ \mathrm{\&Delta;u}(k+\mathrm{aN}-1|k)\end{array}\right]---\left(20\right)\\ +{\mathrm{\&Sigma;}}_{i=0}^{j-1}{A}^i\mathrm{Bu}(k-1)\end{array}$

The predictive equation that to sum up can obtain state X is always:

Finally, the predictive equation of output Y is:

Y(k+j|k)＝CX(k+j|k),j＝1,2,...,P (22)

Step 2.2: rolling optimization.

Rolling optimization is the main feature of MPC, the cost function that is mainly manifested in rolling optimization with the maximum difference of common optimal control algorithm (LQR) is not unalterable overall situation function, but in each sampling instant, cost function only relates to the limited period (prediction time domain and control time domain) from present to future, and to next sampling instant, prediction time domain is passed forward as shown in Figure 5.Therefore MPC each time be carved with a cost function with respect to this moment, the relative form of cost function is not identical in the same time, but its absolute form (is comprised time zone, weight etc.) be different, optimizing is not that an off-line completes, but repeatedly carry out online, the implication of Here it is rolling optimization, is also the essential characteristics that MPC is different from traditional optimum control.MPC adopt this limited period to optimize to take into account due to model mismatch, time the uncertainty that causes such as change, interference, make up in time the impact that these factors cause, and all the time new optimization is based upon on the basis of real process, for leading machine formation control, the equilibrium point of model, the interference of lead aircraft is in continuous variation, therefore, is based upon the feature that Rolling optimal strategy on the limited period meets formation flight more.

If formation flight at the cost function J of moment k is:

$J=\underset{l=1}{\overset{P}{\mathrm{\&Sigma;}}}\left(\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}\right||{\mathrm{\&omega;}}_y(l,m)\&CenterDot;(Y_m(k+l|k)-{\mathrm{ref}}_m(k+l))|{|}^2+\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|\right|{\mathrm{\&omega;}}_{\mathrm{\&Delta;u}}(l,m)\&CenterDot;\mathrm{\&Delta;}{u}_m(k+l|k)\left|{|}^2\right)---\left(23\right)$

Wherein, P is prediction time domain; N is for controlling time domain; Y _m(k+l|k) be under control time m, the predicted value in k+l moment of calculating in k moment; Ref is spacing hold instruction, ref _m(k+l) the spacing hold instruction in the k+l moment of be under control time m, the k moment calculating; ω _yfor output punishment weight matrix, reflect the dynamics of mating of the maintenance dynamics of node W and L spacing and flight path drift angle, speed, be the of paramount importance weight of maintenance of adjusting the distance, ω _y(l, m) represents the output punishment weight that time domain is corresponding to (l, m); ω _{Δ u}for input slew rate punishment weight matrix, reflect the size of instruction rate of change, there is the character of integration, make system realize indifference control, ω _{Δ u}(l, m) represents the input slew rate punishment weight that time domain is corresponding to (l, m); Δ u _m(k+l|k) be the k+l moment input slew rate calculating in the k moment under control time m.For simplifying solving of quadratic programming (Quadratic programming, QP), make each time domain equate the weight of (l, m) here, two weight matrix deteriorate to respectively row vector with ω _x, ω _yrepresent the maintenance dynamics of node W and L spacing, ω _{Δ V}represent the flight path drift angle of node W, the coupling dynamics of speed, ω _{Δ Vc}represent the size of controlled input slew rate.Visible in the situation that not adding end conswtraint, the physical significance of weight is very clear and definite, but in order to ensure the stable of system, has taked larger prediction time domain here, gets P=50 through design of Simulation; N=40; ω _y=(200,250,200,50); ω _{Δ u}=(5,0.1).

Finally, the constrained optimization problem of MPC below the controlled quentity controlled variable u of moment k (k) is:

Wherein, for the restriction of flight path drift angle, represent respectively minimum, maximum yaw angle, for speed command restriction, V _min, V _maxrepresent respectively minimum speed and maximal rate.Δ u _optfor first controlled quentity controlled variable of optimal sequence, u1, u ₂be respectively and V _c.I _{2 × 2}represent 2 × 2 unit matrix, O _{(N-1) × (N-1)}represent 0 matrix of (N-1) × (N-1).

Shape can turn to QP typical problem (quadratic programming problem) suc as formula the cost function of (24), when weight is positive, is Strict Convex, and detailed process adopts more ripe KWIK Algorithm for Solving, is not repeating herein.KWIK algorithm reference paper: Schmid, C., L.T.Biegler.Quadratic programming methods for reduced hessian SQP[J] .Computers & Chemical Engineering, 18 (9): 817832,1994.

Step 2.3: feedback compensation.

In real system, exist non-linear, time the factor such as change, model mismatch, interference, can not conform to completely with actual conditions based on the prediction of varying model not, therefore, feedback policy is indispensable.Rolling optimization is only based upon on the basis of feedback compensation, just can embody its superiority.Therefore, MPC is determining by optimization after the control action in a series of futures, control departing from perfect condition in order to prevent that model mismatch or environmental interference from causing, be not that these control actions are all implemented one by one, and just apply the control action of current time, first controlled quentity controlled variable of the optimal control sequence of trying to achieve, to next sampling instant, first the actual output of monitoring target, and by various feedback policy correction forecast models or compensated, and then carry out new optimization.

In sum, MPC comprehensive utilization historical information and model information, constantly carry out rolling optimization to cost function, and revise or compensation prediction model according to the actual object output recording.

Obtain the current state X (k) of system by supporting network at each computation period, may not be inconsistent with the predicted value X of state (k|k-1), but state space MPC method is the original state using X (k) as system directly, and this has played the effect of similar feedback compensation; In addition as previously mentioned, the flight path drift angle of node L and speed being regarded as can measurements interference also append in the predicted value of forecast model by formula (11) to node W by network delivery, also played the effect of proofreading and correct, and typical PID controller does not process to the interference input of node L.

Step 2.4: have the forecast model in network inducement delay situation.

In the time that network inducement delay be can not ignore, the discretize of linear movement model need to be considered the impact of time delay, investigates the general type of linear time invariant continuous system:

$\left.\begin{array}{c}\frac{\mathrm{dX}\left(t\right)}{\mathrm{dt}}\mathrm{AX}\left(t\right)+\mathrm{Bu}\left(t\right)\\ Y=\mathrm{CX}\left(t\right)\end{array}\right\}---\left(25\right)$

Its time solution is:

$X\left(t\right)=e^{A(t-t_0)}X\left(t_0\right)+{\&Integral;}_{t_0}^t{e}^{A(t-\mathrm{\&mu;})}\mathrm{Bu}\left(\mathrm{\&mu;}\right)\mathrm{d\&mu;}---\left(26\right)$

T ₀the initial time of expression system, in the time there is network inducement delay τ, because the input u of MPC is discrete, the sampling period of establishing u is T _s, in the time that new state does not arrive, u will equal the input in a upper moment, so have

$\frac{\mathrm{dX}\left(t\right)}{\mathrm{dt}}=\begin{array}{cc}\mathrm{AX}\left(t\right)+\mathrm{Bu}\left((k-1)T_s\right),& {\mathrm{kT}}_s\&le;t\&le;{\mathrm{kT}}_s+\mathrm{\&tau;}\\ \mathrm{AX}\left(t\right)+\mathrm{Bu}\left(k{T}_s\right),& k{T}_s+\mathrm{\&tau;}\&le;t<(k+1)T_s\end{array}---\left(27\right)$

Work as kT _s≤ t < kT _swhen+τ, consider that state X is from kT _stransfer to kT _s+ τ, has:

$\begin{array}{c}X(k{T}_s+\mathrm{\&tau;})=e^{\mathrm{A\&tau;}}X\left(k{T}_s\right)+{\&Integral;}_{k{T}_s}^{{\mathrm{kT}}_s+\mathrm{\&tau;}}{e}^{A(k{T}_s+\mathrm{\&tau;}-\mathrm{\&mu;})}\mathrm{Bu}\left((k-1)T_s\right)\mathrm{d\&mu;}\\ =e^{\mathrm{A\&tau;}}X\left(k{T}_s\right)+{\&Integral;}_0^{\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;}\&CenterDot;\mathrm{Bu}\left((k-1)T_s\right)\end{array}---\left(28\right)$

Work as kT _s+ τ≤t < (k+1) T _stime, consider that state X is from kT _s+ τ transfers to (k+1) T _s, have:

$\begin{array}{c}X\left((k+1)T_s\right)=e^{A(T_s-\mathrm{\&tau;})}X({\mathrm{kT}}_s+\mathrm{\&tau;})+{\&Integral;}_{{\mathrm{kT}}_s+\mathrm{\&tau;}}^{(k+1)T_s}{e}^{A((k+1)T_s-\mathrm{\&mu;})}\mathrm{Bu}\left({\mathrm{kT}}_s\right)\mathrm{d\&mu;}\\ =e^{A(T_s-\mathrm{\&tau;})}[e^{\mathrm{A\&tau;}}X\left({\mathrm{kT}}_s\right)+{\&Integral;}_0^{\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;}\&CenterDot;\mathrm{Bu}\left((k-1)T_s\right)]+{\&Integral;}_0^{T_s-\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;}\&CenterDot;\mathrm{Bu}\left(k{T}_s\right)\\ =e^{A{T}_s}X\left(k{T}_s\right)+e^{A(T_s-\mathrm{\&tau;})}{\&Integral;}_0^{\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;B}\&CenterDot;u\left((k-1)T_s\right)+{\&Integral;}_0^{T_s-\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;B}\&CenterDot;u\left(k{T}_s\right)\\ =A_{1}X\left(k{T}_s\right)+B_{1}u\left((k-1)T_s\right)+B_{2}u\left(k{T}_s\right)\end{array}---\left(29\right)$

Wherein:

$\begin{array}{ccc}{A}_1=e^{A{T}_s},& B_1=e^{A(T_s-\mathrm{\&tau;})}{\&Integral;}_0^{\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;B},& B_2=\end{array},{\&Integral;}_0^{T_s-\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;B}---\left(30\right)$

In the situation that not there is not ambiguity, k represents k sampling period, if any Z (k)=Z (kT _s), establish

$Z\left(k\right)=\left(\begin{array}{c}X\left(k\right)\\ u(k-1)\end{array}\right)---\left(31\right)$

Exist in the situation of network inducement delay, the discretization model of MPC is:

$\left.\begin{array}{c}Z(k+1)=A_{z}Z\left(k\right)+B_{z}u\left(k\right)\\ Y\left(k\right)=C_{z}Z\left(k\right)\end{array}\right\}---\left(32\right)$

Wherein

$A_z=\left(\begin{array}{cc}{A}_1& B_1\\ 0& 0\end{array}\right),B_z=\left(\begin{array}{c}{B}_2\\ I\end{array}\right),C_z=\left(\begin{array}{cc}C& 0\end{array}\right)---\left(33\right)$

Then can obtain forecast model.

Below MPC formation control method of the present invention is verified and evaluated.

First, by safe distance surplus compare λ _ritake the motor-driven probability of collision prevention with calculating to form into columns.Suppose the safe distance d of unmanned plane _si=200m, meansquaredeviationσ _ij=10m, the formation setpoint distance d of formation _rij=230m, wherein i, j=1,2 ..., 7.Safe distance surplus ratio according to definition 3, this formation belongs to loose and forms into columns.Between unmanned plane member i and j, need the Probability p of taking collision prevention motor-driven _b(i, j)=α=Φ (3)=0.00135.If forming regular formation by 7 unmanned planes is polygonal formation, i.e. N=7, n _c=6, calculate unmanned plane member by formula (5) and formula (6) and take the motor-driven Probability p of collision prevention _b(i)=0.008073 < 1%, forms into columns and takes the motor-driven Probability p of collision prevention _b=0.05516 < 6%.If other assumed condition is constant, only by the formation setpoint distance of forming into columns closer to safe distance d _si(t), change d into _rij(t)=210m, λ _rij=1, according to definition 5, this formation belongs to close/intra, p _b(i, j)=α=Φ (1)=0.1587, in the case, unmanned plane member takes the motor-driven Probability p of collision prevention _b(i)=0.64543 > 60%, forms into columns and takes the motor-driven Probability p of collision prevention _b=0.999295 ≈ 100%.That is to say, by safe distance surplus than by 3 being reduced to 1 (be safe distance surplus be reduced to 1 σ by 3 σ), form into columns and need to take the motor-driven probability of collision prevention just to increase to close to 100% from being less than 6%, formation type also becomes close/intra from loose formation.Formation need to take the motor-driven probability of collision prevention close to 1, forms into columns and almost must take the motor-driven measure of collision prevention.Visible, safe distance surplus compare λ _ribeing not only the significant parameter of dividing formation type, and being the variable that affects the dynamic stability characteristic of forming into columns, is to judge the density of formation and the important indicator of closeness, is also the important parameter of the collaborative guidance of examination formation controller performance.

Suppose to have the unmanned plane of N=32 to form into columns, on average degree of connection n _c(t)=7, the safe distance d of unmanned plane _si=200m, the formation setpoint distance d of formation _rij=230m, wherein i, j=1,2 ..., 32.

Unmanned plane member's safe distance surplus ratio if require the collision probability of forming into columns to be not more than 30% in statistical significance, i.e. p _b≤ 30%.Can calculate according to formula (1)-Shi (6) according to definition 5, this formation belongs to close/intra, proposes the collaborative meansquaredeviationσ of controller to 32 unmanned aerial vehicle (UAV) control precision that guide of formation of meter _ijshould be not more than 10.17m.

Then,, for the directive function of unmanned plane autonomous formation evaluation index for formation control device design is described, select typical PID formation control device and MPC formation control device proposed by the invention to carry out contrast simulation checking.

Two-shipper is formed into columns and is tested flight path scene: minimum speed 50m/s; Maximal rate 200m/s; Initial velocity 100m/s; Formation keeping is apart from x=100m; Y=-173.2m; D=200m; Negative 60 degree of flight path deflection when t=5s; When t=30s, deflection positive 60 is spent again.

(1) PID formation control device simulation result, as shown in a of Fig. 6, b and c, c is the schematic diagram of unmanned plane relative distance d, a and b be respectively unmanned plane apart from d at certain unmanned aerial vehicle flight path coordinate system x, the component on y axle;

(2) MPC formation control device simulation result, as shown in the A of Fig. 7, B and C, C is the schematic diagram of unmanned plane relative distance d, A and B are respectively unmanned plane relative distance d at certain unmanned aerial vehicle flight path coordinate system x, the component on y axle.

In the simulation time of test flight path 70s, adopt the poor σ=32.7089m of separation criteria of PID formation control device; Adopt the poor σ=0.8393m of separation criteria of MPC formation control device of the present invention.Simulation result shows, PID formation control device differs at most and has exceeded 60m with desired spacing in the process of maintenance of forming into columns, and within MPC formation control device and desired spacing maintain 2.75m, verify the advantage of MPC formation control method provided by the invention in formation flight control.

Claims (5)

1. unmanned plane autonomous formation evaluation index and Model Predictive Control (MPC) formation control method, is characterized in that:

Described unmanned plane autonomous formation evaluation index is as follows:

$\mathrm{\&sigma;}\left(t\right)=\sqrt{E[\mathrm{\&Delta;}{d}_f\left(t\right){]}^2-\{E\left[\mathrm{\&Delta;}{d}_f\left(t\right)\right]{\}}^2}=\sqrt{E[\mathrm{\&Delta;}{d}_f\left(t\right){]}^2-[\mathrm{\&Delta;}{d}_{\mathrm{Ef}}\left(t\right){]}^2}---\left(1\right)$

Wherein, Δ d _f(t) represent the average connection safe distance surplus that the t moment forms into columns; σ (t) is t moment unmanned plane autonomous formation evaluation index, is Δ d _f(t) mean square deviation; Parameter Δ d _ef(t)=E[Δ d _f(t)];

Described MPC formation control method, specific implementation step is as follows:

The first step, set up unmanned plane formation motion model:

If two unmanned plane W and L, speed is respectively V _wand V _l, flight path drift angle is respectively with unmanned plane W is d to the distance of unmanned plane L, x ₀, y ₀for apart from d at unmanned plane W flight path axis system x, the component on y axle; The motion model of unmanned plane L and unmanned plane W is as follows:

and V _wcbe respectively the instruction of flight path drift angle and the speed command of unmanned plane W, and V _lcbe respectively the instruction of flight path drift angle and the speed command of unmanned plane L, be respectively the time constant of unmanned plane L and the response of unmanned plane W flight path drift angle, be respectively the time constant of unmanned plane L and unmanned plane W speed responsive;

Second step, realizes MPC formation control, comprising: forecast model, rolling optimization and feedback compensation;

Step 2.1: forecast model;

Utilize formula (2) as forecast model, get state controlled being input as the flight path drift angle of node L and speed are regarded as can measurements interference output wherein while ignoring D, the discrete form of unmanned plane formation motion model is:

$\left\{\begin{array}{c}X(k+1)=\mathrm{AX}\left(k\right)+\mathrm{Bu}\left(k\right)\\ Y\left(k\right)=\mathrm{CX}\left(k\right)\end{array}\right.---\left(3\right)$

Wherein, X (k) represents the state of moment k, and u (k) represents the controlled input of moment k, and Y (k) represents the output of moment k, and A, B and C represent matrix of coefficients; If u (k+i|k) is illustrated in the controlled input value u in the following k+i moment of moment k supposition, X (k+i|k), Y (k+i|k) is illustrated respectively in moment k, utilize u (k+j|k), (j=0,1, ..., i-1) to X, the predicted value that Y makes; The predictive equation of the state X obtaining is:

Wherein, P is prediction time domain, and N is for controlling time domain, and u (k-1) represents the controlled input in k-1 moment, difference DELTA u (k+j|k)=u (k+j|k)-u (k+j-1|k);

The predictive equation of output Y is: Y (k+j|k)=CX (k+j|k), and j=1,2 ..., P;

Step 2.2: rolling optimization;

If formation flight at the cost function J of moment k is:

$J=\underset{l=1}{\overset{P}{\mathrm{\&Sigma;}}}\left(\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}\right||{\mathrm{\&omega;}}_y(l,m)\&CenterDot;(Y_m(k+l|k)-{\mathrm{ref}}_m(k+l))|{|}^2+\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}\left|\right|{\mathrm{\&omega;}}_{\mathrm{\&Delta;u}}(l,m)\&CenterDot;\mathrm{\&Delta;}{u}_m(k+l|k)\left|{|}^2\right)---\left(4\right)$

Wherein, ω _yfor output punishment weight matrix, ω _y(l, m) represents the output punishment weight that time domain is corresponding to (l, m); Y _m(k+l|k) be under control time m, the predicted value in k+l moment of calculating in k moment; Ref is spacing hold instruction, ref _m(k+l) the spacing hold instruction in the k+l moment of be under control time m, the k moment calculating; ω _{Δ u}for input slew rate punishment weight matrix, ω _{Δ u}(l, m) represents the input slew rate punishment weight that time domain is corresponding to (l, m); Δ u _m(k+l|k) be the k+l moment input slew rate calculating in the k moment under control time m; Make each time domain equate the weight of (l, m), ω _yand ω _{Δ u}deteriorate to respectively a row vector;

Constrained optimization problem below the controlled quentity controlled variable u of moment k (k) is:

Δu _opt＝arg?minJ

Wherein, with be respectively flight path drift angle and speed command restriction; Δ u _optfor first controlled quentity controlled variable of optimal sequence, u ₁u ₂be respectively and V _c; I _{2 × 2}represent 2 × 2 unit matrix, O _{(N-1) × (N-1)}represent 0 matrix of (N-1) × (N-1);

Step 2.3: feedback compensation, comprising: (1) is the original state using X (k) as system directly; (2) the flight path drift angle of unmanned plane L and speed being regarded as can measurements interference by network delivery to unmanned plane W and append to the predicted value of forecast model.

2. unmanned plane autonomous formation evaluation index according to claim 1 and MPC formation control method, is characterized in that: the average connection safe distance surplus Δ d that the described t moment forms into columns _f(t) acquisition methods is:

Δd _f(t)＝d _f(t)-d _sf(t)

$\begin{array}{c}=\frac{1}{M\left(t\right)}\underset{i=1}{\overset{M\left(t\right)}{\mathrm{\&Sigma;}}}\left\{\frac{1}{n_{\mathrm{ci}}\left(t\right)}\underset{j=1,i\&NotEqual;j}{\overset{n_{\mathrm{ci}}\left(t\right)}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}\left(t\right)\right\}-\frac{1}{M\left(t\right)}\underset{i=1}{\overset{M\left(t\right)}{\mathrm{\&Sigma;}}}{d}_{\mathrm{si}}\left(t\right)\\ =\frac{1}{M\left(t\right)}\underset{i=1}{\overset{M\left(t\right)}{\mathrm{\&Sigma;}}}\left\{\right[\frac{1}{n_{\mathrm{ci}}\left(t\right)}\underset{j=1,i\&NotEqual;1}{\overset{n_{\mathrm{ci}}\left(t\right)}{\mathrm{\&Sigma;}}}{d}_{\mathrm{ij}}\left(t\right)]-d_{\mathrm{si}}\left(t\right)\}\end{array}---\left(5\right)$

Wherein, d _f(t) represent the average junctor spacing that the t moment forms into columns, d _sf(t) represent the average safe distance that the t moment forms into columns, M (t) is t moment unmanned plane sum, d _ij(t) represent the distance between unmanned plane i and unmanned plane j, d _si(t) safe distance of expression unmanned plane i, n _ci(t) the connection degree of expression t moment unmanned plane i.

3. unmanned plane autonomous formation evaluation index according to claim 1 and 2 and MPC formation control method, is characterized in that: described unmanned plane autonomous formation is evaluated, by formula (6)-(11) below as judgement schematics, specifically:

If d _ij(t) represent the distance between unmanned plane i and unmanned plane j, d _si(t) safe distance of expression unmanned plane i, the safe distance surplus Δ d of unmanned plane i and unmanned plane j _ij(t)=d _ij(t)-d _si(t), below the time t in the variable relevant to time t is saved and write a Chinese character in simplified form;

If unmanned plane member's safe distance surplus Δ d _ijnormal Distribution Δ d _rijfor the safe distance surplus of setting, σ _ij ²for variance; As Δ d _ij≤ 0 o'clock, unmanned plane member just needed or prepares to take the motor-driven measure of collision prevention, needs or prepare the Probability p of taking collision prevention motor-driven between unmanned plane i and unmanned plane j _b(i, j) is:

p _b(i,j)＝P{Δd _ij≤0}

$\begin{array}{c}={\&Integral;}_{-\&infin;}^{\frac{-\mathrm{\&Delta;}{d}_{\mathrm{rij}}}{{\mathrm{\&sigma;}}_{\mathrm{ij}}}}\frac{1}{\sqrt{2\mathrm{\&pi;}}}\&CenterDot;e^{-\frac{{\mathrm{\&tau;}}^2}{2}\mathrm{d\&tau;}}\\ =\mathrm{\&Phi;}\left(\frac{-\mathrm{\&Delta;}{d}_{\mathrm{rij}}}{{\mathrm{\&sigma;}}_{\mathrm{ij}}}\right)\\ =\mathrm{\&Phi;}(-{\mathrm{\&lambda;}}_{\mathrm{rij}})\end{array}---\left(6\right)$

Between unmanned plane i and unmanned plane j without taking the motor-driven probability q of collision prevention _b(i, j) is:

q _b(i,j)＝P{Δd _ij＞0}

$\begin{array}{c}=1-{\&Integral;}_{-\&infin;}^{\frac{-\mathrm{\&Delta;}{d}_{\mathrm{rij}}}{{\mathrm{\&sigma;}}_{\mathrm{ij}}}}\frac{1}{\sqrt{2\mathrm{\&pi;}}}\&CenterDot;e^{-\frac{{\mathrm{\&tau;}}^2}{2}\mathrm{d\&tau;}}\\ =1-\mathrm{\&Phi;}\left(\frac{-\mathrm{\&Delta;}{d}_{\mathrm{rij}}}{{\mathrm{\&sigma;}}_{\mathrm{ij}}}\right)\\ =1-\mathrm{\&Phi;}(-{\mathrm{\&lambda;}}_{\mathrm{rij}})\end{array}---\left(7\right)$

In formula, safe distance surplus ratio, for normal probability integral;

Closeness is reflected to the motor-driven Probability p of collision prevention of unmanned plane i oneself to the impact of the motor-driven probability of collision prevention _b(i) be:

$\begin{array}{c}{p}_b\left(i\right)\underset{j=1}{\overset{n_{\mathrm{ci}}}{\mathrm{\&Sigma;}}}{p}_b(i,j)[q_b(i,j){]}^{n_{\mathrm{ci}}-1}+\underset{j=1}{\overset{C_{n_{\mathrm{ci}}}^2}{\mathrm{\&Sigma;}}}[p_b(i,j){]}^2[q_b(i,j){]}^{n_{\mathrm{ci}}-2}+\&CenterDot;\&CenterDot;\&CenterDot;\\ +\underset{j=1}{\overset{C_{n_{\mathrm{ci}}}^k}{\mathrm{\&Sigma;}}}\left[p_b(i,j){]}^k\right[q_b(i,j){]}^{n_{\mathrm{ci}}-k}+\&CenterDot;\&CenterDot;\&CenterDot;+\underset{j=1}{\overset{n_{\mathrm{ci}}}{\mathrm{\&Pi;}}}{p}_b(i,j)\\ =\underset{k=1}{\overset{n_{\mathrm{ci}}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{C_{n_{\mathrm{ci}}}^k}{\mathrm{\&Sigma;}}}\left[p_b(i,j){]}^k\right[q_b(i,j){]}^{n_{\mathrm{ci}}-k}\\ =\underset{k=1}{\overset{n_{\mathrm{ci}}}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{C_{n_{\mathrm{ci}}}^k}{\mathrm{\&Sigma;}}}\left[\mathrm{\&Phi;}(-{\mathrm{\&lambda;}}_{\mathrm{rij}}){]}^k\right[1-\mathrm{\&Phi;}(-{\mathrm{\&lambda;}}_{\mathrm{rij}}){]}^{n_{\mathrm{ci}}-k},(n_{\mathrm{ci}}\&GreaterEqual;2)\end{array}---\left(8\right)$

Wherein, n _cifor the connection degree of unmanned plane i, for number of combinations; For connection degree n _ci=1 o'clock, the Probability p that collision prevention is motor-driven _b(i) employing formula (6) is calculated;

The motor-driven Probability p of collision prevention is taked in whole formation _bfor:

$\begin{array}{c}{p}_b=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{p}_b\left(i\right)[1-p_b\left(i\right){]}^{M-1}+\underset{i=1}{\overset{C_M^2}{\mathrm{\&Sigma;}}}[p_b\left(i\right){]}^2[1-p_b\left(i\right){]}^{M-2}+\&CenterDot;\&CenterDot;\&CenterDot;\\ +\underset{i=1}{\overset{C_M^k}{\mathrm{\&Sigma;}}}\left[p_b\left(i\right){]}^k\right[1-p_b\left(i\right){]}^{M-k}+\&CenterDot;\&CenterDot;\&CenterDot;\underset{i=1}{\overset{M}{\mathrm{\&Pi;}}}{p}_b\left(i\right)\\ =\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{C_M^k}{\mathrm{\&Sigma;}}}\left[p_b\left(i\right){]}^k\right[1-p_b\left(i\right){]}^{M-k}\\ =\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{C_M^k}{\mathrm{\&Sigma;}}}\left[p_b\left(i\right){]}^k\right[q_b\left(i\right){]}^{M-k}\end{array}---\left(9\right)$

In formula, q _b(i)=1-p _b(i) for unmanned plane member i is without taking the motor-driven probability of collision prevention, the unmanned plane number that M is whole formation;

Especially, work as p _bwhen (i, j)=α (0≤α≤1) is normal value, formula (8) and formula (9) become respectively formula (10) and formula (11):

$\begin{array}{c}{p}_b\left(i\right)=p=n_{\mathrm{ci}}\&CenterDot;\mathrm{\&alpha;}{\mathrm{\&beta;}}^{n_{\mathrm{ci}}-1}+C_{n_{\mathrm{ci}}}^2{\mathrm{\&alpha;}}^2\&CenterDot;{\mathrm{\&beta;}}^{n_{\mathrm{ci}}-2}+\&CenterDot;\&CenterDot;\&CenterDot;+C_{n_{\mathrm{ci}}}^k\&CenterDot;{\mathrm{\&alpha;}}^k{\mathrm{\&beta;}}^{n_{\mathrm{ci}}-k}+\&CenterDot;\&CenterDot;\&CenterDot;+{\mathrm{\&alpha;}}^{n_{\mathrm{ci}}}\\ =\underset{k=1}{\overset{n_{\mathrm{ci}}}{\mathrm{\&Sigma;}}}{C}_{n_{\mathrm{ci}}}^k\&CenterDot;{\mathrm{\&alpha;}}^k{\mathrm{\&beta;}}^{n_{\mathrm{vi}}-k}\\ =1-{\mathrm{\&beta;}}^{n_{\mathrm{ci}}}\end{array}---\left(10\right)$

$\begin{array}{c}{P}_b=M\&CenterDot;p{q}^{M-1}+C_M^2\&CenterDot;p^2{q}^{M-2}+\&CenterDot;\&CenterDot;\&CenterDot;+C_M^k\&CenterDot;p^k{q}^{M-k}+\&CenterDot;\&CenterDot;\&CenterDot;+p^M\\ =\underset{k=1}{\overset{M}{\mathrm{\&Sigma;}}}{C}_M^k\&CenterDot;p^k{q}^{M-k}\\ =1-q^M\end{array}---\left(11\right)$

Wherein, β=1-α, q=1-p.

4. unmanned plane autonomous formation evaluation index according to claim 1 and MPC formation control method, is characterized in that: the forecast model in described second step, exist in network inducement delay situation, and concrete acquisition methods is:

In the time that network inducement delay can not be ignored, the discretize of linear movement model need to be considered the impact of time delay, and the general type of investigating linear time invariant continuous system is:

$\left.\begin{array}{c}\frac{\mathrm{dX}\left(t\right)}{\mathrm{dt}}\mathrm{AX}\left(t\right)+\mathrm{Bu}\left(t\right)\\ Y=\mathrm{CX}\left(t\right)\end{array}\right\}---\left(12\right)$

Its time solution is:

$X\left(t\right)=e^{A(t-t_0)}X\left(t_0\right)+{\&Integral;}_{t_0}^t{e}^{A(t-\mathrm{\&mu;})}\mathrm{Bu}\left(\mathrm{\&mu;}\right)\mathrm{d\&mu;}---\left(13\right)$

T ₀the initial time of expression system, in the time there is network inducement delay τ, because the input u of MPC is discrete, the sampling period of establishing u is T _s, in the time that new state does not arrive, u will equal the input in a upper moment, so have:

$\frac{\mathrm{dX}\left(t\right)}{\mathrm{dt}}=\begin{array}{cc}\mathrm{AX}\left(t\right)+\mathrm{Bu}\left((k-1)T_s\right),& {\mathrm{kT}}_s\&le;t\&le;{\mathrm{kT}}_s+\mathrm{\&tau;}\\ \mathrm{AX}\left(t\right)+\mathrm{Bu}\left(k{T}_s\right),& k{T}_s+\mathrm{\&tau;}\&le;t<(k+1)T_s\end{array}---\left(14\right)$

Work as kT _s≤ t < kT _swhen+τ, consider that state X is from kT _stransfer to kT _s+ τ, has:

$\begin{array}{c}X(k{T}_s+\mathrm{\&tau;})=e^{\mathrm{A\&tau;}}X\left(k{T}_s\right)+{\&Integral;}_{k{T}_s}^{{\mathrm{kT}}_s+\mathrm{\&tau;}}{e}^{A(k{T}_s+\mathrm{\&tau;}-\mathrm{\&mu;})}\mathrm{Bu}\left((k-1)T_s\right)\mathrm{d\&mu;}\\ =e^{\mathrm{A\&tau;}}X\left(k{T}_s\right)+{\&Integral;}_0^{\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;}\&CenterDot;\mathrm{Bu}\left((k-1)T_s\right)\end{array}---\left(15\right)$

Work as kT _s+ τ≤t < (k+1) T _stime, consider that state X is from kT _s+ τ transfers to (k+1) T _s, have

$\begin{array}{c}X\left((k+1)T_s\right)=e^{A(T_s-\mathrm{\&tau;})}X({\mathrm{kT}}_s+\mathrm{\&tau;})+{\&Integral;}_{{\mathrm{kT}}_s+\mathrm{\&tau;}}^{(k+1)T_s}{e}^{A((k+1)T_s-\mathrm{\&mu;})}\mathrm{Bu}\left({\mathrm{kT}}_s\right)\mathrm{d\&mu;}\\ =e^{A(T_s-\mathrm{\&tau;})}[e^{\mathrm{A\&tau;}}X\left({\mathrm{kT}}_s\right)+{\&Integral;}_0^{\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;}\&CenterDot;\mathrm{Bu}\left((k-1)T_s\right)]+{\&Integral;}_0^{T_s-\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;}\&CenterDot;\mathrm{Bu}\left(k{T}_s\right)\\ =e^{A{T}_s}X\left(k{T}_s\right)+e^{A(T_s-\mathrm{\&tau;})}{\&Integral;}_0^{\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;B}\&CenterDot;u\left((k-1)T_s\right)+{\&Integral;}_0^{T_s-\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;B}\&CenterDot;u\left(k{T}_s\right)\\ =A_{1}X\left(k{T}_s\right)+B_{1}u\left((k-1)T_s\right)+B_{2}u\left(k{T}_s\right)\end{array}---\left(16\right)$

Wherein parameter:

$\begin{array}{ccc}{A}_1=e^{A{T}_s},& B_1=e^{A(T_s-\mathrm{\&tau;})}{\&Integral;}_0^{\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;B},& B_2=\end{array},{\&Integral;}_0^{T_s-\mathrm{\&tau;}}{e}^{\mathrm{A\&mu;}}\mathrm{d\&mu;B}---\left(17\right)$

K represents k sampling period, Z (k)=Z (kT _s), establish

$Z\left(k\right)=\left(\begin{array}{c}X\left(k\right)\\ u(k-1)\end{array}\right)---\left(18\right)$

Exist in the situation of network inducement delay, the discretization model of MPC is:

$\left.\begin{array}{c}Z(k+1)=A_{z}Z\left(k\right)+B_{z}u\left(k\right)\\ Y\left(k\right)=C_{z}Z\left(k\right)\end{array}\right\}---\left(19\right)$

Wherein,

$A_z=\left(\begin{array}{cc}{A}_1& B_1\\ 0& 0\end{array}\right),B_z=\left(\begin{array}{c}{B}_2\\ I\end{array}\right),C_z=\left(\begin{array}{cc}C& 0\end{array}\right)---\left(20\right)$

Finally obtain forecast model according to the discretization model of formula (19) and (20).

5. according to the unmanned plane autonomous formation evaluation index described in claim 1 or 4 and MPC formation control method, it is characterized in that: described prediction time domain P=50; Described control time domain N=40; Described output punishment weights omega _y=(200,250,200,50), described input slew rate punishment weights omega _{Δ u}=(5,0.1).