CN105222648A - A kind of linear pseudo-spectrum broad sense mark control miss distance Guidance and control method - Google Patents

Abstract

The invention discloses a kind of linear pseudo-spectrum broad sense mark control miss distance Guidance and control method, it is the concept based on broad sense mark control miss distance, in conjunction with None-linear approximation Model Predictive Control, linear quadratic optimum control and the pseudo-spectrometry of Gauss, compose discrete by linearisation and Gauss's puppet, original nonlinear optimal control problem is converted into the problem of one group of continuous solving linear algebraic equation systems.Advantage of the present invention is to have very high computational efficiency to solving of optimal control problem, only need several node just can obtain good computational accuracy, and final solution can be expressed as the smooth function about controlling discrete nodes, is highly suitable for line computation.The method being applied to the latter end with terminal point constraint attacks in guidance, simulation result shows, relative to MPSP method, the present invention not only has higher computational efficiency and precision, and can be applicable in the guidance framework of terminal guidance completely, relative to self adaptation latter end proportional guidance, the present invention will produce less required overload instruction.

Description

A kind of linear pseudo-spectrum broad sense mark control miss distance Guidance and control method

Technical field

The invention belongs to the Guidance and control technical field of aircraft, be specifically related to a kind of linear pseudo-spectrum broad sense mark control miss distance Guidance and control method.

Background technology

Broad sense mark control miss distance (GNEM) refers to from t ₀moment starts aircraft according to standard control u _pflight, at t _fthe state of flight of moment aircraft and the deviation of given last current state.Broad sense mark control miss distance not only comprises traditional miss distance, also can comprise terminal impingement angle deviation, even comprises terminal impact velocity deviation.The acquisition of broad sense mark control miss distance, can based on the thinking of Model Predictive Control, from off-line to solve optimal control problem (optimal control problem consisted of the full nonlinear equation of direct solution) different, Model Predictive Control normally solves passes through linearisation, and one group that obtains is closed on optimal control problem, often through solving a two-point boundary value problem, and obtain current controlled quentity controlled variable.

In order to increase computational efficiency, Ohtsuka and Fujii has expanded stable continuous control method, thus obtain one can the optimization method of line solver Control of Nonlinear Systems.LuPing proposes a kind of path tracking method of closed loop, and former problem, by multistep Taylor expansion and Euler-Simpson integration method, is converted into a unconfined quadratic programming problem, obtains the analytical relation of current controlled quentity controlled variable afterwards by the method.Yan utilizes Legendre pseudo-spectrometry Discrete Linear Quadratic optimal control, initial control just can be obtained afterwards by resolving the one group of linear algebraic equation be transformed by former problem, and the satellite attitude stabilization under being successfully applied to magnetic fields, the method hypothesis state error is an a small amount of relative to former reference state.In order to overcome this drawback, PaulWilliams adopts Quasilinearization Method that former problem is converted into one group of linear optimal control problem, then the pseudo-spectrometry of Jacobi this problem discrete is adopted, it is converted into the problem that one group of Algebraic Equation set solves, rolling time horizon integration is for overcoming larger state deviation.

Method mentioned above is all based on roll stablized loop, this control method can only solve the local optimum control problem of a Short Time Domain, and therefore, method can obtain local solution fast, but can not ensure Global Optimality, these class methods also depend on the good base of trajectory of an off-line.

In order to solve the nonlinear optimal control problem with strong end conswtraint and quadratic performance index, India scholar Padhi proposes the NONLINEAR OPTIMAL CONTROL method that one is called model prediction static programming (MPSP), the method combines the theoretical thought with approaching Dynamic Programming of non-linear mould predictive, just the control sequence meeting global optimum can be obtained by solving static programming problem continuously, but the method needs to select a large amount of discrete nodes just can have enough Euler's integral precision.

Summary of the invention

The object of the present invention is to provide a kind of linear pseudo-spectrum broad sense mark control miss distance Guidance and control method (LPGNEMG & C), state variable and association's state variable are expressed as the linear combination that one group take Lagrange interpolation polynomial as basic function by the method, can be just one group of Algebraic Constraint by differential dynamics equation constraints conversion by positive mating point, therefore, the linear optimal control problem with quadratic performance index is just converted into the problem solving one group of linear algebraic equation.The method is relative to numerical solution nonlinear programming problem, and the time loss of solving linear algebric equation is very little, has very high computational efficiency in this way.

For achieving the above object, the present invention by the following technical solutions.

A kind of linear pseudo-spectrum broad sense mark control miss distance Guidance and control method, the method comprises the steps:

(1) initialize: initial computer sim-ulation parameter is set, by trajectory optimisation or the rational primary standard control sequence of guidance emulation acquisition of off-line;

(2) trajectory integration is predicted: use the quantity of state of current time as initial value for integral, described primary standard control sequence is as control inputs, and the controlled quentity controlled variable wherein corresponding to current time is designated as U ₀, carry out prediction trajectory integration, obtain broad sense mark control miss distance (SOT state of termination constraint deviation) d _ψand overall trajectory message X _k, U _k;

(3) precision of described broad sense mark control miss distance is judged: if described d _ψmeet the required precision that described step (1) is arranged, enter step 4; If do not meet required precision, enter step 6;

(4) initially control U is upgraded ₀: judge described U ₀whether meet and control border restriction, as described U ₀exceed boundary Control U _max, described U ₀equal U _max, otherwise, described U ₀equal U ₀, enter step 5;

(5) control instruction U is performed ₀: by described U ₀be applied in actual system and go, return step 2, and by current control sequence U _kas the control inputs of next step trajectory integration;

(6) control sequence upgrades: around prediction integration trajectory, carry out linearization process, in conjunction with optimum control first order necessary condition, and use pseudo-spectrometry carry out discrete after, resolve linear algebraic equation and obtain control sequence and upgrade.

Control method as above, preferably, described step (6) specifically comprises:

First, around prediction integration trajectory, carry out linearization process, in conjunction with the first order necessary condition of optimum control, linear optimal control problem be converted into the two-point boundary value problem (TPBVP) meeting dynamics constraint condition:

$\left[\begin{array}{c}\mathrm{\δ}\stackrel{\·}{x}\\ \stackrel{\·}{\mathrm{\λ}}\end{array}\right]=\left[\begin{array}{cc}A& -{\mathrm{BR}}^{-1}{B}^T\\ -Q& -A^T\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}x\\ \mathrm{\λ}\end{array}\right]+\left[\begin{array}{c}B{u}_p\\ {\mathrm{Qx}}_p\end{array}\right]---\left(1\right)$

Wherein, δ u=u _p-R ^-1b ^tλ;

Secondly, by described equation (1) from time domain [t ₀, t _f] be transformed between [-1,1] by time-varying function, and discrete at LG Nodes:

$\left[\begin{array}{c}\mathrm{\δ}\stackrel{\·}{x}\\ \stackrel{\·}{\mathrm{\λ}}\end{array}\right]=\left[\begin{array}{cc}\frac{t_f-t_0}{2}A& -\frac{t_f-t_0}{2}{\mathrm{BR}}^{-1}{B}^T\\ -\frac{t_f-t_0}{2}Q& -\frac{t_f-t_0}{2}{A}^T\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}x\\ \mathrm{\λ}\end{array}\right]+\frac{t_f-t_0}{2}\left[\begin{array}{c}B{u}_p\\ {\mathrm{Qx}}_p\end{array}\right]---\left(2\right)$

$\left\{\begin{array}{c}\underset{l=0}{\overset{N}{\Σ}}{D}_{kl}{\mathrm{\δx}}_l-\frac{t_f-t_0}{2}\left(A_k{\mathrm{\δx}}_k-B_k{R}^{-1}{B}_k^T{\mathrm{\λ}}_k\right)=\frac{t_f-t_0}{2}{B}_k{u}_{pk}\\ \underset{l=1}{\overset{N+1}{\Σ}}{D}_{kl}^{*}{\mathrm{\λ}}_l+\frac{t_f-t_0}{2}\left(Q_k{\mathrm{\δx}}_k+A_k^T{\mathrm{\λ}}_k\right)=\frac{t_f-t_0}{2}{Q}_k{x}_{pk}\end{array}\right.---\left(3\right)$

$u_k=u_{pk}-{\mathrm{\δu}}_k=R_k^{-1}{B}_k{\mathrm{\λ}}_k---\left(4\right)$

Wherein, $u_k=u_{pk}-{\mathrm{\δu}}_k=R_k^{-1}{B}_k{\mathrm{\λ}}_k,k=1,2,...,N;$

Arrange described equation (3), described linear optimal control problem is converted to one group of linear algebraic equation, and the analytical expression of the state variable obtained on described LG node and association's state variable:

Sz＝K(5)

Wherein, described z is the column vector about described state variable and described association state variable, being specifically expressed as follows of element in described S and K:

$S=\left[\begin{array}{cc}{S}^{xx}& S^{x\mathrm{\λ}}\\ S^{\mathrm{\λ}x}& S^{\mathrm{\λ}\mathrm{\λ}}\end{array}\right];K=\left[\begin{array}{cc}{\left(K^x\right)}^T& {\left(K^{\mathrm{\λ}}\right)}^T\end{array}\right]---\left(6\right)$

Association state variable λ is solved by described equation (5) _k, and substituted in described equation (4), the control sequence u after renewal _kcan separate.

Below by matrix S and the solution throughway of K and the renewal process of controlled quentity controlled variable in detailed description the method.The linearisation of nonlinear dynamical equation is one of foremost method in commercial Application, although a large amount of documents is linearization technique described in more detail, still be necessary again to describe in detail before derivation, consider the Kind of Nonlinear Dynamical System as follows with strong end conswtraint

$\stackrel{\·}{x}=f(x,u,t)---\left(7\right)$

ψ(x(t _f))＝0(8)

Wherein, x ∈ R ⁿstate vector, u ∈ R ^mbe dominant vector, t ∈ R is time variable, ψ ∈ R ^sbe terminal function constraint vector, based on one group of reference trajectory, equation (7) carried out the Taylor expansion of high-order, and neglect higher order term, just can obtain one group of linear dynamics equation (9) using state deviation as quantity of state

$\mathrm{\δ}\stackrel{\·}{x}=A\mathrm{\δ}x+B\mathrm{\δ}u---\left(9\right)$

Wherein, A=f _x(x _p, u _p) and B=f _u(x _p, u _p), x _pand u _pfor standard state and control variables, δ x is the state deviation between reference trajectory, it should be noted that, in the linearization procedure of nonlinear equation, the relation between virtual condition and standard state is expressed as x=x _p-δ x, with the x=x described in other documents _p+ δ x is different.

The optimal control problem that this method can solve, its performance indications are necessary for the form of Linear-Quadratic Problem.Equation (10) describes the general type of such performance indications.

$\begin{array}{c}J={\mathrm{\δx}}^T\left(t_f\right)P_f\mathrm{\δ}x\left(t_f\right)+v^T\left(\frac{\∂\mathrm{\ψ}}{\∂x_f}\mathrm{\δ}x\left(t_f\right)-d\mathrm{\ψ}\right)+\\ \frac{1}{2}{\∫}_{t_0}^{t_f}{\left(x_p-\mathrm{\δ}x\right)}^{T}Q\left(x_p-\mathrm{\δ}x\right)+{\left(u_p-\mathrm{\δ}u\right)}^{T}R\left(u_p-\mathrm{\δ}u\right)dt\end{array}---\left(10\right)$

Two is terminal capabilities index above, P _ffor the SOT state of termination weight function matrix of positive semidefinite, state error between the SOT state of termination constraint that d ψ is the acquisition of employing ballistic prediction integration and the SOT state of termination needed for reality retrain, Q and R is respectively state and controls the weight matrix in performance indications integration, is all the matrix of positive semidefinite.It is pointed out that these performance indications are actual state and the weighted sum of squares of control, instead of the weighted sum of squares between state deviation, although adopt the weighted sum of squares of state deviation to bring calculating into do not affect derivation.

According to the first order necessary condition of optimum control, a set condition equation (11) can be obtained, (12), (13) and (14).

$H=\frac{1}{2}{(x_p-\mathrm{\δ}x)}^{T}Q(x_p-\mathrm{\δ}x)+\frac{1}{2}{(u_p-\mathrm{\δ}u)}^{T}R(u_p-\mathrm{\δ}u)+{\mathrm{\λ}}^T(A\mathrm{\δ}x+B\mathrm{\δ}u)---\left(11\right)$

$\frac{\∂H}{\∂\mathrm{\δ}u}=-R(u_p-\mathrm{\δ}u)+B^T\mathrm{\λ}=0\→\mathrm{\δ}u=u_p-R^{-1}{B}^T\mathrm{\λ}---\left(12\right)$

$-\stackrel{\·}{\mathrm{\λ}}=\frac{\∂H}{\∂\mathrm{\δ}x}=-Q(x_p-\mathrm{\δ}x)+A^T\mathrm{\λ}\→\stackrel{\·}{\mathrm{\λ}}=-Q\mathrm{\δ}x-A^T\mathrm{\λ}+{\mathrm{Qx}}_p---\left(13\right)$

$\mathrm{\δ}\stackrel{\·}{x}=A\mathrm{\δ}x-{\mathrm{BR}}^{-1}{B}^T\mathrm{\λ}+{\mathrm{Bu}}_p---\left(14\right)$

Wherein, H represents the hamilton's function of linear dynamics equation, and λ is association's state variable of current sequence.

For such problem, there are two kinds of typical truncation conditions, one is for work as x _fwhen being present in terminal capabilities target function, equation (15) is for calculating association's state function value on corresponding border.

$\mathrm{\λ}\left(t_f\right)=P_f\mathrm{\δ}x\left(t_f\right)+v^T\left(\frac{\∂\mathrm{\ψ}}{\∂x_f}\right)---\left(15\right)$

Another is exactly work as x _fwhen not being present among terminal capabilities target function, in order to eliminate the impact of corresponding terminal state, assuming that at weight function P _fin, that a part of component of this state is zero, and so, corresponding border association state function value is also zero.

λ(t _f)＝0(16)

If when the SOT state of termination is constrained to a particular value, be just equivalent to the special circumstances of the first typical transversality condition, SOT state of termination constraint just can be expressed as x (t _f)=x _f, corresponding border association state function value is just expressed as (17):

λ(t _f)＝v(17)

Wherein, v is Lagrange multiplier.

In brief, linear optimal control problem is just converted into the two-point boundary value problem (TPBVP) that meets Dynamic Constraints (9) condition, and the matrix form that this problem can be described as shown in equation (18), its transversality condition is determined by SOT state of termination constraint (8).

$\left[\begin{array}{c}\mathrm{\δ}\stackrel{\·}{x}\\ \stackrel{\·}{\mathrm{\λ}}\end{array}\right]=\left[\begin{array}{cc}A& -{\mathrm{BR}}^{-1}{B}^T\\ -Q& -A^T\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}x\\ \mathrm{\λ}\end{array}\right]+\left[\begin{array}{c}B{u}_p\\ {\mathrm{Qx}}_p\end{array}\right]---\left(18\right)$

So, just can obtain optimum state variable and association's state variable value by solving the two-point boundary value problem with initial and end conswtraint.Corresponding optimum control also obtains by equation (12).

To linear puppet spectrum broad sense mark control miss distance Guidance and control method (LPGNEMG & C) proposed by the invention be introduced, for resolving such optimal control problem fast below.

In the method, state variable and association's state variable are expressed as the linear combination that one group take Lagrange interpolation polynomial as basic function, can be just one group of Algebraic Constraint by differential dynamics equation constraints conversion by positive mating point, therefore, the linear optimal control problem with quadratic performance index is just converted into the problem solving one group of linear algebraic equation.Relative to numerical solution nonlinear programming problem, the time loss of solving linear algebric equation is very little, will have very high computational efficiency.Linear pseudo-spectrum broad sense mark control miss distance Guidance and control algorithm is adopted to carry out discrete concrete steps as follows.

Without loss of generality, consider the as above linear time varying system mentioned in a joint, its Dynamic Constraints and performance indications are respectively (9) and (10).The strong point due to Lagrange interpolation polynomial is all orthogonal points, and all between [-1,1], the discrete first step is by time zone [t ₀, t _f] be transformed between [-1,1] by time-varying function.

$t=\frac{t_f-t_0}{2}\mathrm{\τ}+\frac{t_f+t_0}{2}---\left(19\right)$

Original performance indications and kinetics equation just can be converted to equation (20) and (21).

$J={\mathrm{\δx}}^T\left(1\right)P_f\mathrm{\δ}x\left(1\right)+v^T\left(\frac{\∂\mathrm{\ψ}}{\∂x_f}\mathrm{\δ}x\right(1)-d\mathrm{\ψ})+\frac{t_f-t_0}{4}{\∫}_{-1}^1\begin{array}{c}{(x_p-\mathrm{\δ}x)}^{T}Q(x_p-\mathrm{\δ}x)+\\ {(u_p-\mathrm{\δ}u)}^{T}R(u_p-\mathrm{\δ}u)d\mathrm{\τ}\end{array}---\left(20\right)$

$\mathrm{\δ}\stackrel{\·}{x}=\frac{t_f-t_0}{2}A\mathrm{\δ}x+\frac{t_f-t_0}{2}B\mathrm{\δ}u---\left(21\right)$

Corresponding two-point boundary value problem is just converted into equation (22), (23) and (24).

$\left[\begin{array}{c}\mathrm{\δ}\stackrel{\·}{x}\\ \stackrel{\·}{\mathrm{\λ}}\end{array}\right]=\left[\begin{array}{cc}\frac{t_f-t_0}{2}A& -\frac{t_f-t_0}{2}{\mathrm{BR}}^{-1}{B}^T\\ -\frac{t_f-t_0}{2}Q& -\frac{t_f-t_0}{2}{A}^T\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}x\\ \mathrm{\λ}\end{array}\right]+\frac{t_f-t_0}{2}\left[\begin{array}{c}B{u}_p\\ {\mathrm{Qx}}_p\end{array}\right]---\left(22\right)$

δx(-1)＝δx ₀(23)

$\mathrm{\λ}\left(1\right)=P_f\left(\mathrm{\δ}x\right(1\left)\right)+v^T\left(\frac{\∂\mathrm{\ψ}}{\∂x_f}\right)---\left(24\right)$

δ x ₀for initial state deviation vector, by the linearisation of prediction trajectory anomalous integral, continuous print solves the solution that above linear optimal control problem just can obtain primary nonlinear optimal control problem.

Because trajectory integration is using current state amount as initial value, so, initial state deviation is just zero, therefore, and δ x ₀zero is always in computational process.

L _n(τ) and L ^* _n(τ) Lagrange interpolation polynomial on N rank is respectively, τ _l(l=1,2 ..., N) and be the strong point of N rank Lagrange interpolation polynomial, it is the zero point of Legnedre polynomial, is also referred to as LG node.Asking in process the solution formula acquisition not having to resolve, the LG node needed for obtaining can only calculated by numerical algorithm.Afterwards, state variable, association's state variable and control variables just can be expressed as one group of Lagrange interpolation polynomial is the linear combination of basic function.

${\mathrm{\δx}}^N\left(\mathrm{\τ}\right)=\underset{l=0}{\overset{N}{\Σ}}\mathrm{\δ}x\left({\mathrm{\τ}}_l\right)L_l\left(\mathrm{\τ}\right);{\mathrm{\δu}}^N\left(\mathrm{\τ}\right)=\underset{l=0}{\overset{N}{\Σ}}\mathrm{\δ}u\left({\mathrm{\τ}}_l\right)L_l\left(\mathrm{\τ}\right);{\mathrm{\λ}}^N\left(\mathrm{\τ}\right)=\underset{l=1}{\overset{N+1}{\Σ}}\mathrm{\λ}\left({\mathrm{\τ}}_l\right)L_l^{*}\left(\mathrm{\τ}\right)---\left(25\right)$

According to the basic theories of Lagrange interpolation polynomial, L _n(τ) and L ^* _n(τ) following character is met.

$L_l\left({\mathrm{\τ}}_k\right)=\left\{\begin{array}{c}1,l=k\\ 0,l\≠k\end{array}\right.;L_l^{*}\left({\mathrm{\τ}}_k\right)=\left\{\begin{array}{c}1,l=k\\ 0,l\≠k\end{array}\right.---\left(26\right)$

δx ^N(τ _l)＝δx(τ _l)；δu ^N(τ _l)＝δu(τ _l)；λ ^N(τ _l)＝λ(τ _l)(27)

At each x of LG Nodes ⁿjust can be expressed as differential to the derivative of time and approach the quantic that matrix is multiplied with the state of start node with LG node.Differential approaches the matrix of a matrix D formula N × (N+1), and by obtaining each part differentiate of Lagrange interpolation polynomial, each element in D is:

$\mathrm{\δ}{\stackrel{\·}{x}}^N\left({\mathrm{\τ}}_k\right)=\underset{l=0}{\overset{N}{\Σ}}{D}_{kl}\mathrm{\δ}x\left({\mathrm{\τ}}_l\right);{\stackrel{\·}{\mathrm{\λ}}}^N\left({\mathrm{\τ}}_k\right)=\underset{l=1}{\overset{N+1}{\Σ}}{D}_{kl}^{*}\mathrm{\λ}\left({\mathrm{\τ}}_l\right)---\left(28\right)$

$D_{kl}={\stackrel{\·}{L}}_i\left({\mathrm{\τ}}_k\right)=\underset{l=0}{\overset{N}{\Σ}}\frac{\underset{j=0,j\≠i,l}{\overset{N}{\Π}}({\mathrm{\τ}}_k-{\mathrm{\τ}}_j)}{\underset{j=0,j\≠i}{\overset{N}{\Π}}({\mathrm{\τ}}_k-{\mathrm{\τ}}_j)}---\left(29\right)$

Differential approaches matrix D and its adjoint differential and approaches matrix D * and have following relation:

$D_{ik}^{*}=-\frac{{\mathrm{\ω}}_k}{{\mathrm{\ω}}_i}{D}_{kl}---\left(30\right)$

And approaching in matrix D at differential, there is following relation in first, last column element and other column elements:

${\stackrel{\‾}{D}}_i=-\underset{k=1}{\overset{N}{\Σ}}{D}_{ik}---\left(31\right)$

Wherein, for differential approaches the first, last column element of matrix, according to above two formulas, approach matrix D * with differential and just directly can approach matrix D by differential and determined.

Assuming that state variable and association's state variable have following form.

${\mathrm{\δx}}_l=\left(\begin{array}{ccc}{\mathrm{\δx}}_1^T& ...& {\mathrm{\δx}}_{N+1}^T\end{array}\right);{\mathrm{\λ}}_l=\left(\begin{array}{ccc}{\mathrm{\λ}}_0^T& {\mathrm{\λ}}_1^T& {\mathrm{\λ}}_{N+1}^T\end{array}\right)---\left(32\right)$

Be brought into by equation (28) in equation (22), the kinetics equation of state and association's state is not only converted into one group of algebraic equation, and they also can be expressed as the function of state variable and association's state variable on LG node.

$\left\{\begin{array}{c}\underset{l=0}{\overset{N}{\Σ}}{D}_{kl}{\mathrm{\δx}}_l-\frac{t_f-t_0}{2}\left(A_k{\mathrm{\δx}}_k-B_k{R}^{-1}{B}_k^T{\mathrm{\λ}}_k\right)=\frac{t_f-t_0}{2}{B}_k{u}_{pk}\\ {\mathrm{\δx}}_{N+1}={\mathrm{\δx}}_0+\frac{t_f-t_0}{2}\underset{k=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_k\left(A_k{\mathrm{\δx}}_k-B_k{R}^{-1}{B^T}_k{\mathrm{\λ}}_k+B_k{u}_{pk}\right)\\ \underset{l=1}{\overset{N+1}{\Σ}}{D}_{kl}^{*}{\mathrm{\λ}}_l+\frac{t_f-t_0}{2}\left(Q_k{\mathrm{\δx}}_k+A_k^T{\mathrm{\λ}}_k\right)=\frac{t_f-t_0}{2}{Q}_k{x}_{pk}\\ {\mathrm{\λ}}_0={\mathrm{\λ}}_{N+1}+\frac{t_f-t_0}{2}\underset{k=1}{\overset{N}{\Σ}}{\mathrm{\ω}}_k\left(Q_k{\mathrm{\δx}}_k+A_k^T{\mathrm{\λ}}_k-Q_k{x}_{pk}\right)\end{array}\right.---\left(33\right)$

Wherein, k=1,2 ..., N, it is to be noted, Dynamic Constraints is only configured on LG node, and does not comprise first and last two end points, therefore, utilize Gaussian integral formula and Associative algcbra constraint equation (28), two additional constraints is in order to ensure that state variable on first, last end points and association's state variable meet Dynamic Constraints.

Arrange equation (33), linear optimal control problem is just converted to one group of linear algebraic equation, and the analytical expression of the state variable that can obtain on LG node and association's state variable:

Sz＝K(34)

Wherein, z is the column vector about state variable and association's state variable, being specifically expressed as follows of element in S and K:

$S=\left[\begin{array}{cc}{S}^{xx}& S^{x\mathrm{\λ}}\\ S^{\mathrm{\λ}x}& S^{\mathrm{\λ}\mathrm{\λ}}\end{array}\right];K=\left[\begin{array}{cc}{\left(K^x\right)}^T& {\left(K^{\mathrm{\λ}}\right)}^T\end{array}\right]---\left(35\right)$

Wherein,

In linear algebraic equation, K value is made up of constant value element, and is made up of two parts, and a part is the part about SOT state of termination constraint, and another part is the integration trajectory state about prediction, and its expression is as follows:

K＝-K ₁+K ₂(40)

Wherein,

$K_1=\left[\begin{array}{c}{S}_{1+s(N+1):m+s(N+1)}^{xx}\\ S_{1+s(N+1):m+s(N+1)}^{\mathrm{\λ}x}\end{array}\right]{\mathrm{\δx}}_f---\left(41\right)$

$\begin{array}{c}{K}_2=\left[\begin{array}{cccc}\frac{t_f-t_0}{2}{B}_1{u}_{p1}& ...& \frac{t_f-t_0}{2}{B}_N{u}_{pN}& {\Σ}_1^N{\mathrm{\ω}}_k\frac{t_f-t_0}{2}{B}_k{u}_{pk}\end{array}\right.\\ {\left.\begin{array}{cccc}-{\Σ}_1^N{\mathrm{\ω}}_k\frac{t_f-t_0}{2}{Q}_k{x}_{pk}& \frac{t_f-t_0}{2}{Q}_1{x}_{p1}& ...& \frac{t_f-t_0}{2}{Q}_N{x}_{pN}\end{array}\right]}^T\end{array}---\left(42\right)$

Wherein, s is the state variable number of system, and m is in the about intrafascicular state number of the SOT state of termination.

It should be noted that the end conswtraint function in linear optimal control problem needs the functional form concrete according to it to represent.Therefore, be difficult to carry out unified expression in superincumbent matrix, but, when end conswtraint is the particular value about quantity of state, only need to assist state variable to consider as independent variable accordingly.Afterwards, the number due to constraint equation equals unknown state number, by the solution of existence anduniquess.

When the number of the state variable not in terminal capabilities index is k, just can determine to assist state variable accordingly according to transversality condition, the number of known variables also will be reduced to (2N+3) s-m-k, and equal the number of constraint equation.State variable on LG node and association's state variable just can be obtained by the system of linear equations above solving, and equally, the control variables on LG node can be calculated by following formula.

$u_k=u_{pk}-{\mathrm{\δu}}_k=R_k^{-1}{B}_k{\mathrm{\λ}}_k---\left(43\right)$

The step repeated above just can obtain the solution of nonlinear optimal control problem, but, along with real system is close to target, the number of LG node and the step number of prediction trajectory integration all will reduce, this method is very easy to realize, and is combined with linearization of nonlinear system and Model Predictive Control by pseudo-for Gauss spectrometry, therefore, not only remain the solving precision of indirect method, also due to need not nonlinear programming problem be solved, there is good computational efficiency.

Linear pseudo-spectrum broad sense mark control miss distance Guidance and control method (LPGNEMG & C) of one of the present invention, it is the concept based on broad sense mark control miss distance, in conjunction with None-linear approximation Model Predictive Control, linear quadratic optimum control and the pseudo-spectrometry of Gauss, by the discrete and linearisation of Gauss's puppet spectrum, original nonlinear optimal control problem is converted into the problem of one group of continuous solving linear algebraic equation systems.

Advantage of the present invention is to have very high computational efficiency to solving of optimal control problem, only need several node just can obtain good computational accuracy, and final solution can be expressed as the smooth function about controlling discrete nodes, is highly suitable for line computation.

The area-wide optimal control of satisfied strong end conswtraint that what the present invention and Paul method maximum different were to obtain is therefore, be not the locally optimal solution based on rolling time horizon, and linearisation neither based on a fixing reference trajectory.Another one is not both, and relative to the discrete approximation differential matrix of the pseudo-spectrometry of Legendre and the pseudo-spectrometry of Jacobi, the discrete approximation differential matrix of the pseudo-spectrometry of Gauss used in the present invention has completeness, can not cause numerical solution difficulty due to discrete method.The present invention by testing in the end game guidance simulation of particular orientation, simulation result shows, relative to MPSP method and self adaptation end proportional guidance, the present invention not only has higher computational efficiency and precision, and can be applicable to completely in the guidance framework of terminal guidance.

Accompanying drawing explanation

Fig. 1 is the implementing procedure figure of control method of the present invention.

Fig. 2 is the motion schematic diagram of guided missile in large ground level.

Fig. 3 is three-dimensional flight path schematic diagram.

Fig. 4 is that in different example situation, change curve is in time transshipped in longitudinal instruction.

Fig. 5 is that in different example situation, change curve is in time transshipped in horizontal instruction.Fig. 6 is the trajectory tilt angle change curve in time in different example situation.

Fig. 7 is the trajectory course angle change curve in time in different example situation.Fig. 8 is the end error convergency factor of different example.

Fig. 9 is that self adaptation proportional guidance compares with the trajectory tilt angle of LPGNEMG & C.

Figure 10 is that self adaptation proportional guidance compares with the course angle of LPGNEMG & C.

Figure 11 is that self adaptation proportional guidance compares with the actual overload of LPGNEMG & C.

Figure 12 is that self adaptation proportional guidance compares with total overload of LPGNEMG & C.

Detailed description of the invention

Below in conjunction with accompanying drawing and example, the present invention is described in further detail.

The inventive method is applied in the terminal guidance of the air-to-ground guided missile ensureing direction of attack, and adopts three dimensional non-linear engagement model.Use line straightening machine moving-target Example Verification the method for different azimuth for the applicability in different angle of attack direction, the simulation example of Different L G node illustrates the impact of LG node on guidance precision and computational efficiency.The guidance effect of the method also contrasts with other Guidance Law, and the simulation result of all equal usage ratio guidings of example is brought guidance into as initial value and is calculated, whole simulated programs is complete in the personal computer of 3.3GHz and under MATLAB2008b simulated environment at CPU frequency, and more efficient translation and compiling environment will improve computational efficiency.

1. air-to-ground guided missile terminal guidance problem describes

Under plane earth supposed situation, the three-dimensional particle movement kinetics equation of standard is expressed as:

$\begin{array}{c}{\stackrel{\·}{x}}_m=V_m{\mathrm{cos\γ}}_m{\mathrm{cos\ψ}}_m;{\stackrel{\·}{y}}_m=V_m{\mathrm{cos\γ}}_m{\mathrm{sin\ψ}}_m;{\stackrel{\·}{z}}_m=V_m{\mathrm{sin\γ}}_m;\\ {\stackrel{\·}{V}}_m=-\frac{D_m}{m_m}-g{\mathrm{sin\γ}}_m;{\stackrel{\·}{\mathrm{\γ}}}_m=\frac{-a_z-g{\mathrm{cos\γ}}_m}{V_m};{\stackrel{\·}{\mathrm{\ψ}}}_m=\frac{-a_y}{V_m{\mathrm{cos\γ}}_m};\\ {\stackrel{\·}{a}}_y=\frac{\left(a_{yc}-a_y\right)}{t_{\mathrm{\τ}}};{\stackrel{\·}{a}}_z=\frac{\left(a_{zc}-a_z\right)}{t_{\mathrm{\τ}}};\end{array}---\left(44\right)$

It should be noted that the carryover effects that a joint link produces for analog auto pilot link, wherein, a _zcand a _ycbe respectively z to y to command acceleration, a _zand a _ybe respectively guided missile z to y to suffered actual acceleration, and the acceleration of two is all perpendicular to velocity direction.

What need proposition is in guidance solution process, and in order to reduce computing time, increase computational efficiency, consideration one does not save the impact of link.X _m, y _m, z _mbe respectively three coordinate positions of aircraft at large ground level, γ _mfor trajectory tilt angle, be the angle of velocity and XY plane projection, ψ _mfor course angle curve, it is the angle that speed projection planar and x-axis are formed.In order to better understanding, Fig. 2 will pass through the motion of view shows guided missile intuitively.

The simulation parameter relevant to aircraft is (as quality, area of reference etc.) be all content in existing document, and all quantity of states are all carried out normalized by emulation, and the object of normalized resolves environment to obtain better numerical value exactly.

$\begin{array}{c}{x}_n=\frac{x_m}{x_{*}};y_n=\frac{y_m}{y_{*}};z_n=\frac{z_m}{z_{*}};a_{yn}=\frac{a_{ym}}{a_{y*}};\\ a_{zn}=\frac{a_{zm}}{a_{z*}};V_n=\frac{V_m}{V_{*}};{\mathrm{\γ}}_n=\frac{{\mathrm{\γ}}_m}{{\mathrm{\γ}}_{*}};{\mathrm{\ψ}}_n=\frac{{\mathrm{\ψ}}_m}{{\mathrm{\ψ}}_{*}};\end{array}---\left(45\right)$

Wherein, subscript " n " represents the value after corresponding state normalization, subscript " * " represents the normalization size of corresponding state, and what these needs were careful chooses to make corresponding state and other states all on the same order of magnitude, and its value and some special simulation parameter are described in Table 1.

Table 1 simulation parameter is arranged

Because higher height atmospheric density is thinner, thus causing the steerage of guided missile lower, the acceleration of generation will be very limited, thus be necessary the acceleration process constraints of increase by about height, assuming that the peak acceleration of guided missile in two planes is in length and breadth limited by following formula.

$a_{max}=\left\{\begin{array}{cc}2,& ifz\≥4000\\ -\frac{z}{2000}+4,& ifz\≤4000\end{array}\right.---\left(46\right)$

Increase this process constraints guidance environment that not only more real reflection is actual, and can verify that linear pseudo-spectrum broad sense mark control miss distance Guidance and control algorithm is controlling the guidance effect under saturation conditions.

2. the Selection Strategy of initial value and object module

In this section, use and expand the primary standard control sequence that proportional guidance carries out guiding the linear pseudo-spectrum broad sense mark control miss distance Guidance and control algorithm of emulation acquisition, controlled quentity controlled variable is an array, from current time until deadline, not in the same time, corresponding different controlled quentity controlled variables, so controlled quentity controlled variable is string number, so be called control sequence, the controlled quentity controlled variable that wherein current time is corresponding is exactly the aforementioned U mentioned ₀, because proportional guidance is foremost Guidance Law, and have a large amount of documents to be described in detail it, will not tire out at this and state its Guidance Law principle, equation (47) and (48) be respectively z to y to acceleration instruction.

$a_{zc}=\left\{\begin{array}{c}{N}_e{V}_c{\stackrel{\·}{\mathrm{\σ}}}_{pitch}+\frac{1}{2}{a}_{tpitch}+g{\mathrm{cos\γ}}_m,if\left|a_{zc}\right|\≤a_{cmax}\\ a_{cmax}sign\left(a_{zc}\right),if\left|a_{zc}\right|>a_{c\max }\end{array}\right.---\left(47\right)$

$a_{yc}=\left\{\begin{array}{c}{N}_e{V}_c{\stackrel{\·}{\mathrm{\σ}}}_{yaw}+\frac{1}{2}{a}_{tyaw},if\left|a_{yc}\right|\≤a_{cmax}\\ a_{cmax}sign\left(a_{yc}\right),if\left|a_{yc}\right|>a_{cmax}\end{array}\right.---\left(48\right)$

Wherein, with for the rate of change of the angle of sight is in the projection of two planes, V _cfor guided missile closing speed, N _efor expanding the guidance coefficient of proportional guidance, a _cmaxfor the acceleration of aircraft limits, a _tpitchand a _tyawbe respectively the acceleration function of target, and be expressed as follows:

$\begin{array}{c}{a}_{tpitch}=-{\mathrm{sin\γ}}_m({\mathrm{cos\ψ}}_m{\stackrel{\·\·}{x}}_t+{\mathrm{sin\ψ}}_m{\stackrel{\·\·}{y}}_t)\\ a_{tyaw}=-{\mathrm{sin\γ}}_m{\stackrel{\·\·}{x}}_t+{\mathrm{cos\γ}}_m{\stackrel{\·\·}{y}}_t\end{array}---\left(49\right)$

The speed that is assumed to be about target is constant value, and planar only have one perpendicular to the acceleration of speed, the kinetics equation of target is expressed as:

${\stackrel{\·}{x}}_t=V_t{\mathrm{cos\ψ}}_t;{\stackrel{\·}{y}}_t=V_t{\mathrm{sin\ψ}}_t;{\stackrel{\·}{\mathrm{\ψ}}}_t=\frac{a_{yt}}{V_t};---\left(50\right)$

The simulation example of 3 line straightening machine moving-targets

The motor-driven simulation example of different target straight line, to be used for verifying linear pseudo-spectrum broad sense mark control miss distance Guidance and control algorithm, applicability in the terminal guidance with the constraint of angle of attack orientation angles, except constraint terminal position, aircraft must hit the mark with specific angle of attack direction, the original state of aircraft and target, required end angle of attack orientation angles and LG Node configuration as shown in table 2

Initial and the SOT state of termination of the different simulation example of table 2

Example	ψ t(deg)	γ mf deg	ψ mf deg	LG node number
					Example-1	15	-65	-40	10
Example-2	0	-75	-15	10
					Example-3	-10	-30	20	10

The three-dimensional aerial vehicle trajectory of different example as shown in Figure 3, as can be seen from the figure, linear puppet spectrum broad sense mark control miss distance Guidance and control algorithm can hit the target of movement by directing aircraft smoothly, and the angle, target initial heading of different example is respectively 15 degree, 0 degree and-10 degree.

Curve is over time transshipped in the longitudinal direction that Fig. 4 and Fig. 5 is respectively different example situation and horizontal instruction, aircraft alters a great deal at flight initial period command acceleration, this is owing to there is first-order lag link, initial command acceleration is zero, be because the SOT state of termination error at initial time is by the maximum moment for whole flight on the other hand, need larger acceleration correction could meet required the end game orientation angles requirement.

Also very large in the change of attacking time instructions acceleration, this is also because whole emulation all considers the first-order lag of acceleration instruction, when command acceleration is non-vanishing, effect caused by the deviation between actual acceleration and instruction acceleration will sharply increase at the eleventh hour, thus make aircraft need larger command acceleration could eliminate the impact of error generation, therefore, the miss distance that the first-order lag link selecting the rational guidance cycle that time constant can be made to be 0.2s produces is limited in very little scope, in miss distance statistics afterwards, such inference result will be reflected.Fig. 4, Fig. 5 disclose the later stage that the energy management with terminal attack angle constraint terminal guidance mainly concentrates on flight too.

For different examples, command acceleration will produce very large difference after 20 seconds.Another it should be noted that the command acceleration in example-3 has contacted the boundary value of command acceleration in flight early stage, approximately experienced by the saturated flight of command acceleration of nearly 6 seconds, the command acceleration that guidance algorithm calculates is less than the restriction of maximum command acceleration, finally, all examples can hit the mark according to the orientation of setting separately by directing aircraft smoothly, even if simulation result also illustrate that there is the saturated restriction of instruction, linear pseudo-spectrum broad sense mark control miss distance Guidance and control algorithm also can well guide online.

Fig. 6, Fig. 7 are trajectory tilt angle and the course angle curve over time of different example situation, can find out easily, and trajectory tilt angle and course angle all meet required requirement at end, γ _mf=-65deg, ψ _mf=-40deg; γ _mf=-75deg, ψ _mf=-15deg; γ _mf=-30deg, ψ _mf=20deg.

Equally, the Changing Pattern of and instruction acceleration is compared, and two angles occur that larger difference time will after difference appears in acceleration, and this will mean that the energy management process ensureing terminal angle of attack orientation angles is a macrocyclic process,

In order to study the convergence efficiency of linear pseudo-spectrum broad sense mark control miss distance Guidance and control algorithmic method, another group numerical simulation will be carried out, the simulation parameter of three class examples arrange as table 1 and table 2 introduce the same, but the single order link of instruction overload will not be considered in simulations, Fig. 8 is the logarithm change curve of end error convergence efficiency with update times of different example, it should be noted that dYN ²for the weighted sum of squares of different terminals error, simulation result shows, only needs little a few step iteration, just can obtain very high computational accuracy, and therefore, linear pseudo-spectrum broad sense mark control miss distance Guidance and control algorithm has very high computational efficiency.

Table 3 is the deviation statistics of the terminal location of different example, terminal trajectory tilt angle and terminal course angle, in order to ensure that aircraft can to hit the mark with specifically attacking orientation, all quantity of states except speed are all considered in guidance computational process as end conswtraint, can see, the position deviation of terminal can not be greater than 0.02m, and the deviation of two terminal point also can not be greater than 0.05 degree.These three examples can obtain following inference: terminal point constraint (especially terminal trajectory tilt angle) is less, by the less position deviation of generation and terminal point deviation.

The middle lonely deviation statistics of the different example of table 3

Another needs the problem of research to be exactly about the impact of LG interstitial content on guidance effect, typically, the number of LG node and the number of state equation determine the scale of linear algebraic equation systems, and when state variable smooth enough, state variable just can high-precisionly by the interpolation polynomial that exponent number is less be approached, use the discrete solution that just can obtain degree of precision of the Gauss of less node puppet spectrum, this uses less solution that just can obtain degree of precision computing time by meaning, therefore, LG node number will be a very important simulation parameter.In order to the time efficiency of linear pseudo-spectrum broad sense mark control miss distance Guidance and control algorithm is described, by comparing of the averaging instruction rise time of LPGNEMG & C and MPSP method, and result is displayed by table 4.

The instruction rise time contrast of table 4 distinct methods

Initial and the SOT state of termination of simulation example and the same of example-1, statistics shows, along with the minimizing of LG node, the averaging instruction rise time will significantly reduce, when LG node number is 6, averaging instruction computing time is 0.01131 second, the scale that this result discloses linear algebraic equation too along with LG node number minimizing exponentially form reduce.It is pointed out that and adopt MPSP and LPGNEMG & C to guide, can obtain less SOT state of termination error, therefore, compared with MPSP, LPGNEMG & C will have higher computational efficiency and precision.

4 with the comparing of self adaptation proportional guidance

Self adaptation proportional guidance is a kind of proportional guidance of new model, can directing aircraft with specific angle target of attack, command acceleration is obtained by resize ratio steering coefficient according to current guided missile, dbjective state and required angle-of-attack, and the proportional guidance in three-dimensional planar has following form:

${\stackrel{\·}{\mathrm{\ψ}}}_{com}=-{\mathrm{\λ}}_1\stackrel{\·}{\mathrm{\θ}}---\left(51\right)$

${\stackrel{\·}{\mathrm{\γ}}}_{com}=-{\mathrm{\λ}}_2\stackrel{\·}{\mathrm{\φ}}---\left(52\right)$

Wherein, λ ₁and λ ₂be respectively the self adaptation proportional guidance coefficient in longitudinal and transverse guidance plane, θ and φ is respectively the angle that sight line is formed in XY plane and XZ plane, and they can be calculated by the current state of guided missile and target the derivative of time.It is worthy of note to only have and work as λ ₁and λ ₂when being greater than 2 respectively, the self adaptation proportional guidance of place plane could hit the mark with specific angle by effective directing aircraft.These two parameters are also determine according to the current location of aircraft and velocity direction, therefore, before entering self adaptation proportional guidance, need to arrange two guidance logics, enter the guidanuce condition meeting self adaptation proportional guidance for directing aircraft at different guidance plane.Once λ ₁or λ ₂be greater than 2, the self adaptation proportional guidance of respective planes will activate.

The need command acceleration that relatively more linear pseudo-spectrum broad sense mark control miss distance Guidance and control algorithm and self adaptation proportional guidance generate, simulation parameter is as shown in table 1 too, and target is that terminal point is constrained to γ without any motor-driven static object _mf=-65deg, ψ _mf=-40deg.The constant value command acceleration of a 2g is using the guidance logic as transverse plane, and it is 40 degree that aircraft completes the course angle preset value laterally aimed at, and at horizontal alignment stage, aircraft longitudinally will fly according to the overload planning of following (53).

$n_{com}=2-2e^{-t/T_{trans}}---\left(53\right)$

Wherein, t is the flight time, T _transfor the parameter selected in advance smoothly transits in order to ensure that instruction crosses to be loaded between 0 and 2g.

The terminal location deviation that two method of guidances are formed all is less than 0.2m, therefore, meets guidance requirement completely.Fig. 9, Figure 10 are trajectory tilt angle and the course angle curve over time of LPGNEMG & C and self adaptation proportional guidance two kinds of methods, as can be seen from Figure 7, two kinds of methods can meet terminal point requirement smoothly, equally also clearly disclose self adaptation proportional guidance respectively at 3 seconds with within 12 seconds, enter into self adaptation proportional guidance stage of corresponding guidance plane separately, this figure also illustrates that LPGNEMG & C algorithm will produce more smooth trajectory tilt angle and course angle curve.

Figure 11, Figure 12 generates by two kinds of methods the command acceleration of Different Plane and the total overload change curve with horizontal journey, for self adaptation proportional guidance, to reflect more obvious on acceleration diagram by the transformer effect guidance logical stage to terminal guidance, laterally guiding switch instant, acceleration instruction has jumped to 2g from-2g, at longitudinal guidance switch instant, acceleration instruction has jumped to-6g from 2g, but, for LPGNEMG & C, acceleration instruction change unusual light, in whole guidance process, do not have to experience larger acceleration change.

From the change curve of total overload, can be easy to find out, total mistake that the method produces in guidance process is loaded in the most of the time the total overload being less than self adaptation proportional guidance and producing, therefore, linear puppet spectrum broad sense mark control miss distance Guidance and control algorithm is relative to self adaptation proportional guidance, to more smooth overload instruction and need less overload demand be produced, thus be convenient to automatic pilot and follow the tracks of.

Claims (2)

1. a linear pseudo-spectrum broad sense mark control miss distance Guidance and control method, it is characterized in that, the method comprises the steps:

(1) initialize: initial computer sim-ulation parameter is set, by trajectory optimisation or the rational primary standard control sequence of guidance emulation acquisition of off-line;

(2) trajectory integration is predicted: use the quantity of state of current time as initial value for integral, described primary standard control sequence is as control inputs, and the controlled quentity controlled variable wherein corresponding to current time is designated as U ₀, carry out prediction trajectory integration, obtain broad sense mark control miss distance d ψ and overall trajectory message X _k, U _k;

(3) precision of described broad sense mark control miss distance is judged: if described d ψ meets the required precision that described step (1) is arranged, enter step 4; If do not meet required precision, enter step 6;

(4) initially control U is upgraded ₀: judge described U ₀whether meet and control border restriction, as described U ₀exceed boundary Control U _max, described U ₀equal U _max, otherwise, described U ₀equal U ₀, enter step 5;

(5) control instruction U is performed ₀: by described U ₀be applied in actual system and go, return described step 2, and by current control sequence U _kas the control inputs of next step trajectory integration;

(6) control sequence upgrades: around prediction integration trajectory, carry out linearization process, in conjunction with optimum control first order necessary condition, and use pseudo-spectrometry carry out discrete after, resolve linear algebraic equation and obtain control sequence and upgrade.

2. control method as claimed in claim 1, it is characterized in that, described step (6) specifically comprises:

First, around prediction integration trajectory, carry out linearization process, in conjunction with the first order necessary condition of optimum control, linear optimal control problem be converted into the two-point boundary value problem meeting dynamics constraint condition:

$\left[\begin{array}{c}\mathrm{\δ}\stackrel{\·}{x}\\ \stackrel{\·}{\mathrm{\λ}}\end{array}\right]=\left[\begin{array}{cc}A& -{\mathrm{BR}}^{-1}{B}^T\\ -Q& -A^T\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}x\\ \mathrm{\λ}\end{array}\right]+\left[\begin{array}{c}{\mathrm{Bu}}_p\\ {\mathrm{Qx}}_p\end{array}\right]---\left(1\right)$

Wherein, δ u=u _p-R ^-1b ^tλ;

Secondly, by described equation (1) from time domain [t ₀, t _f] be transformed between [-1,1] by time-varying function, and discrete at LG Nodes:

$\left[\begin{array}{c}\mathrm{\δ}\stackrel{\·}{x}\\ \stackrel{\·}{\mathrm{\λ}}\end{array}\right]=\left[\begin{array}{cc}\frac{t_f-t_0}{2}A& -\frac{t_f-t_0}{2}{\mathrm{BR}}^{-1}{B}^T\\ -\frac{t_f-t_0}{2}Q& -\frac{t_f-t_0}{2}{A}^T\end{array}\right]\left[\begin{array}{c}\mathrm{\δ}x\\ \mathrm{\λ}\end{array}\right]+\frac{t_f-t_0}{2}\left[\begin{array}{c}{\mathrm{Bu}}_p\\ {\mathrm{Qx}}_p\end{array}\right]---\left(2\right)$

$\{\begin{array}{c}\underset{l=0}{\overset{N}{\Σ}}{D}_{kl}{\mathrm{\δx}}_l-\frac{t_f-t_0}{2}\left(A_k{\mathrm{\δx}}_k-B_k{R}^{-1}{B}_k^T{\mathrm{\λ}}_k\right)=\frac{t_f-t_0}{2}{B}_k{u}_{pk}\\ \underset{l=1}{\overset{N+1}{\Σ}}{D}_{kl}^{*}{\mathrm{\λ}}_l+\frac{t_f-t_0}{2}\left(Q_k{\mathrm{\δx}}_k+A_k^T+{\mathrm{\λ}}_k\right)=\frac{t_f-t_0}{2}{Q}_k{x}_{pk}\end{array}---\left(3\right)$

$u_k=u_{pk}-{\mathrm{\δu}}_k=R_k^{-1}{B}_k{\mathrm{\λ}}_k---\left(4\right)$

Wherein, $u_k=u_{pk}-{\mathrm{\δu}}_k=R_k^{-1}{B}_k{\mathrm{\λ}}_k,k=1,2,...,N;$

Arrange described equation (3), described linear optimal control problem is converted to one group of linear algebraic equation, and the analytical expression of the state variable obtained on described LG node and association's state variable:

Sz＝K(5)

Wherein, described z is the column vector about described state variable and described association state variable, being specifically expressed as follows of element in described S and K:

$S=\left[\begin{array}{cc}{S}^{xx}& S^{x\mathrm{\λ}}\\ S^{\mathrm{\λ}x}& S^{\mathrm{\λ}\mathrm{\λ}}\end{array}\right];K=\[\begin{array}{cc}{\left(K^x\right)}^T& {\left(K^{\mathrm{\λ}}\right)}^T\end{array}\]---\left(6\right)$

Association state variable λ is solved by described equation (5) _k, and substituted in described equation (4), the control sequence u after renewal _kcan separate.