NEURAL NETWORK-BASED NARMA-L2 MULTIVARIABLE CONTROL METHOD

Abstract

The present invention discloses a neural network-based NARMA-L2 multivariable control method, which includes: deriving a NARMA-L2 multivariable control law of a nonlinear discrete system; using neural networks to offline identify nonlinear functions in the control law; designing the controller with the control law, designing a closed-loop control system for a certain type of turbofan engine, and adding the function of making the online update of neural network errors on the basis of this to correct the controller parameters, and studying the tracking performance of the control system. Aimed to minimize the quadratic performance index for the errors of the engine's high-pressure rotating speed and pressure ratio, by using neural networks in combination with the deduced NARMA-L2 multivariable control law, the present invention has designed a dual-variable controller, calculated the engine's control variables and controlled the fuel flow of the engine and the critical interface area of the exhaust nozzle.

Claims

1. A neural network-based NARMA-L2 multivariable control method, comprising following steps: Step A) deducing the NARMA-L2 equations and multivariable control law of a multi-input multi-output (MIMO) nonlinear discrete system; Step B) using neural networks to offline identify nonlinear functions in the muitivariable control law; using the identified neural networks to design a dual-variable controller, establishing a dual-variable closed-loop control system of a turbofan engine, and verifying a stable and dynamic state performance of the dual-variable controller, wherein specific sub-steps of the Step A) are as follows: Step A1) obtaining NARMA equations of the MIMO system through recursion according to a state space description of a nonlinear discrete system, and apply it to the MIMO system; Step A2) according to the NARMA equations of the MIMO system, performing multivariate Taylor expansion near an equilibrium point and ignore Taylor high-order remainders to obtain the NARMA-L2 equations of the NEMO nonlinear system; Step A3) performing matrix operations by using the NARMA-L2 equations of the MIMO nonlinear system to obtain the multivariable control law of the MIMO system, wherein the state space description of the MIMO nonlinear discrete system in the Step A1) is as follows: $x\left[k+1\right]=f\left(x\left[k\right],u\left[k\right]\right)y_1\left[k\right]=h_1\left(x\left[k\right]\right)y_2\left[k\right]=h_2\left(x\left[k\right]\right).Math.y_m\left[k\right]=h_m\left(x\left[k\right]\right)$ wherein, u[k]=[u.sub.1[k]u.sub.2[k]. . . . u.sub.m[k]]ϵR.sup.m is a control variable input of the system, y.sub.1[k]ϵR, y.sub.2[k]ϵR, . . . , y.sub.m[k]ϵR is a controlled variable output of the system, m represents a number of the selected control variables and controlled variables, a state vector x[k]ϵR.sup.n, n represents a dimension of the state vector, functions f(.), h.sub.1(.), . . . , h.sub.m(.) ϵC.sup.∞, and an original point is an equilibrium state, and k represents time; the NARMA equations of the MIMO system are: $\begin{array}{c}{y}_1\left[k+1\right]=F\left(y_1\left[k-n+1\right],y_1\left[k-n+2\right],.Math.,y_1\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k\right]\right)\\ y_2\left[k+1\right]=G\left(y_2\left[k-n+1\right],y_2\left[k-n+2\right],.Math.,y_2\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k\right]\right)\\ .Math.\\ y_m\left[k+1\right]=H\left(y_m\left[k-n+1\right],y_m\left[k-n+2\right],.Math.,y_m\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k\right]\right)\end{array}$ wherein, an input variable a u[k]=[u.sub.1[k]u.sub.2[k]. . . . u.sub.m[k]]ϵR.sup.m, output values y.sub.1[k]ϵR, y.sub.2[k]ϵR, . . . , y.sub.m[k]ϵR, functions F(.), G(.), H(.)ϵC.sup.∞ and k represents time, wherein the NARMA-L2 equations of the MIMO system in the Step A2) are: $\begin{array}{c}{y}_1\left[k+1\right]=F\left(M_1,u_1\left[k\right],u_2\left[k\right],.Math.,u_m\left[k\right]\right)\\ =F\left(M_1,0,0,.Math.,0\right)\\ {{{+\frac{\partial F}{\partial u_1\left[k\right]}.Math.}_{\left(M_1,0,0,.Math.,0\right)}{u}_1\left[k\right]+\frac{\partial F}{\partial u_2\left[k\right]}.Math.}_{\left(M_1,0,0,.Math.,0\right)}{u}_2\left[k\right]+.Math.+\frac{\partial F}{\partial u_m\left[k\right]}.Math.}_{\left(M_1,0,0,.Math.,0\right)}{u}_m\left[k\right]\end{array}$ $\begin{array}{c}{y}_2\left[k+1\right]=G\left(M_2,u_1\left[k\right],u_2\left[k\right],.Math.,u_m\left[k\right]\right)\\ =G\left(M_2,0,0,.Math.,0\right)\\ {{{+\frac{\partial G}{\partial u_1\left[k\right]}.Math.}_{\left(M_2,0,0,.Math.,0\right)}{u}_1\left[k\right]+\frac{\partial G}{\partial u_2\left[k\right]}.Math.}_{\left(M_2,0,0,.Math.,0\right)}{u}_2\left[k\right]+.Math.+\frac{\partial G}{\partial u_m\left[k\right]}.Math.}_{\left(M_2,0,0,.Math.,0\right)}{u}_m\left[k\right]\\ .Math.\end{array}$ $\begin{array}{c}{y}_m\left[k+1\right]=H\left(M_2,u_1\left[k\right],u_2\left[k\right],.Math.,u_m\left[k\right]\right)\\ =H\left(M_m,0,0,.Math.,0\right)\\ {{{+\frac{\partial H}{\partial u_1\left[k\right]}.Math.}_{\left(M_m,0,0,.Math.,0\right)}{u}_1\left[k\right]+\frac{\partial H}{\partial u_2\left[k\right]}.Math.}_{\left(M_m,0,0,.Math.,0\right)}{u}_2\left[k\right]+.Math.+\frac{\partial H}{\partial u_m\left[k\right]}.Math.}_{\left(M_m,0,0,.Math.,0\right)}{u}_m\left[k\right]\end{array}$ wherein, $M_1=\left(y_1\left[k-n+1\right],y_1\left[k-n+2\right],.Math.,y_1\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k-1\right]\right)M_2=\left(y_2\left[k-n+1\right],y_2\left[k-n+2\right],.Math.,y_2\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k-1\right]\right).Math.M_m=\left(y_m\left[k-n+1\right],y_m\left[k-n+2\right],.Math.,y_m\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k-1\right]\right).$

2. (canceled)

3. The neural network-based NARMA-L2 multivariable control method according to claim 1, wherein specific sub-steps of the Step B) are as follows: Step B1) first, taking the fuel flow of the turbofan engine and a critical cross-sectional area of an exhaust nozzle as control variables, and taking a high-pressure rotating speed and a pressure ratio of the turbofan engine as controlled variables, using an engine model to generate input and output data, and using neural networks to offline identify the nonlinear functions in the control law; Step B2) according to the control law, constructing a matrix with outputs of the nonlinear functions already identified by the neural networks, establishing a closed-loop controller for the high-pressure rotating speed and the pressure ratio of the turbofan engine, and obtain values of the control variables online, establishing quadratic performance index according to an error between actual output and a command output, and optimizing neural network parameters online; Step B3) after the dual-variable controller is designed, verifying a tracking performance of the dual-variable controller at design points and non-design points and within a full flight envelope.

4-5. (canceled)

6. The neural network-based NARMA-L2 multivariable control method according to claim 1, wherein the multivariable control law of the system in the Step A3) is: $\left[\begin{array}{c}{u}_1\left(k\right)\\ u_2\left(k\right)\\ .Math.\\ u_m\left(k\right)\end{array}\right]={\left[\begin{array}{cccc}{F}_1& F_2& .Math.& F_m\\ G_1& G_2& .Math.& G_m\\ .Math.& .Math.& & .Math.\\ H_1& H_2& .Math.& H_m\end{array}\right]}^{-1}\left(\left[\begin{array}{c}{y}_1\left(k+1\right)\\ y_2\left(k+1\right)\\ .Math.\\ y_m\left(k+1\right)\end{array}\right]-\left[\begin{array}{c}{F}_0\\ G_0\\ .Math.\\ H_0\end{array}\right]\right)$ ${{{{{{{F_m=\frac{\partial F}{\partial u_m\left(k\right)}.Math.}_{\left(M_1,0,0,.Math.,0\right)},G_0=G,G_1=\frac{\partial G}{\partial u_1\left(k\right)}.Math.}_{\left(M_2,0,0,.Math.,0\right)},G_2=\frac{\partial G}{\partial u_2\left(k\right)}.Math.}_{\left(M_2,0,0,.Math.,0\right)},G_m=\frac{\partial G}{\partial u_m\left(k\right)}.Math.}_{\left(M_2,0,0,.Math.,0\right)},H_0=H,H_1=\frac{\partial H}{\partial u_1\left(k\right)}.Math.}_{\left(M_m,0,0,.Math.,0\right)},H_2=\frac{\partial H}{\partial u_2\left(k\right)}.Math.}_{\left(M_m,0,0,.Math.,0\right)},H_m=\frac{\partial H}{\partial u_m\left(k\right)}.Math.}_{\left(M_m,0,0,.Math.,0\right)}.$

7. The neural network-based NARMA-L2 multivariable control method according to claim 3, wherein a function of the quadratic performance index established in the Step B2) is: $$\begin{gathered} J=\frac{1}{2}\left(e_{c1}^2\left[k\right]+e_{c2}^2\left[k\right]\right) \\ {e}_{c1}=y_1\left[k\right]-r_1\left[k\right] \\ {e}_{c2}=y_2\left[k\right]-r_2\left[k\right] \end{gathered}$$ wherein, r.sub.1 and r.sub.2 represent high-pressure rotating speed and commands of the pressure ratio of the turbofan engine, respectively; and y.sub.1 and y.sub.2 represent actual high-pressure rotating speed and actual pressure ratio of the turbofan engine, respectively; after the dual-variable controller obtaining a command of controlled variable at each moment, the dual-variable controller calculates the control variable corresponding to the command according to the control law and inputs the control variable into the engine, the engine then outputs a real value of the controlled variable, in combination with the quadratic performance index, the command of the controlled variable and an actual output value are used for an online correction of the neural network that approximates the nonlinear function, so as to achieve the adaptive control of the engine in the full flight envelope.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0024] FIG. 1 is the offline identification structure diagram of NARMA-L2 neural network;

[0025] FIG. 2 is the off-line identification error curve diagram of the neural network;

[0026] FIG. 3 is the structure diagram of NARMA-L2 dual-variable closed-loop control system;

[0027] FIG. 4 is the tracking response diagram of the rotating speed outputted at the design point (H=0 m and Ma=0);

[0028] FIG. 5 is the tracking response diagram of the pressure ratio outputted at the design point (H=0 m and Ma=0);

[0029] FIG. 6 is the tracking response diagram of the rotating speed outputted at the non-design point (H=4,000 m and Ma=0.8);

[0030] FIG. 7 is the tracking response diagram of the pressure ratio outputted at the non-design point (H=4,000 m and Ma=0.8);

[0031] FIG. 8 is the change diagram of altitude and Mach number during the take-off, climb, descent, acceleration and deceleration processes of the aircraft;

[0032] FIG. 9 is the tracking response diagram of the rotating speed outputted during the take-off, climb, descent, acceleration and deceleration processes of the aircraft;

[0033] FIG. 10 is the tracking response diagram of the pressure ratio outputted during the take-off, climb, descent, acceleration and deceleration processes of the aircraft.

DESCRIPTION OF THE EMBODIMENTS

[0034] The technical solution of the present invention will be further illustrated in detail below in conjunction with the accompanying drawings.

[0035] The idea of the present invention is first to obtain the NARMA equations of the MIMO nonlinear discrete system through recursion and to perform the multivariate Taylor expansion at the equilibrium point to obtain the NARMA-L2 equations of the system, and then to get the control law from the equations through simple derivation. In response to the multi-variable control requirements of the turbofan engine, the dual-variable controller is designed by using the derived NARMA-L2 multi-variable control law to approximate the nonlinear functions in the control law based on the neural network on the basis of the component-level model of the turbofan engine. Compared with the decentralized loop dual-variable control system, the controller eliminates the impact of coupling and has better tracking effects and a little higher accuracy. Compared with other multivariable control methods, this control method has a simple structure and can easily design the control law, and also has good stable and dynamic state performance.

[0036] The embodiment of the present invention takes the design of a dual-variable controller involving the high-pressure rotating speed and pressure ratio of a certain type of turbofan engine as an example. The NARMA-L2 multivariable control method includes the following steps:

[0037] Step A) Deduce the NARMA-L2 equations and multivariable control law of the MIMO nonlinear discrete system;

[0038] Step B) Design a NARMA-L2 dual-variable controller and perform performance verification at design points and non-design points and within the full flight envelope.

[0039] Wherein, the specific sub-steps of Step A) are as follows:

[0040] Step A1) The state space of the MIMO nonlinear discrete system is described as:

[00001]$\begin{array}{c}x\left[k+1\right]=f\left(x\left[k\right],u\left[k\right]\right)\\ y_1\left[k\right]=h_1\left(x\left[k\right]\right)\\ y_2\left[k\right]=h_2\left(x\left[k\right]\right)\\ .Math.\\ y_m\left[k\right]=h_m\left(x\left[k\right]\right)\end{array}$

[0041] Wherein, the input variable u[k]=[u.sub.1[k]u.sub.2[k]. . . u.sub.m[k]]ϵR.sup.m, the state vector x[k]ϵR.sup.m, the output values, y.sub.1[k]ϵR, y.sub.2[k]ϵR, . . . , y.sub.m[k]ϵR, the functions f(.),h.sub.1(.), . . . , h.sub.m(.)ϵC.sup.∞, and the original point is the equilibrium state, and k is the time.

[0042] In practical applications, because usually only input and output signals are obtained, it is needed to establish the input and output description of the system. Get the following results through recursion:

[00002]$\begin{array}{c}{y}_1\left[k\right]=h_1\left(x\left[k\right]\right)={\Psi }_1\left(x\left[k\right]\right)\\ y_1\left[k+1\right]=h_1\left(x\left[k+1\right]\right)=h_1\left(f\left(x\left[k\right],u\left[k\right]\right)\right)={\Psi }_2\left(x\left[k\right],u\left[k\right]\right)\\ y_1\left[k+2\right]=h_1\left(x\left[k+2\right]\right)=h_1\left(f\left(x\left[k+1\right],u\left[k+1\right]\right)\right)\\ =h_1\left(f\left(f\left(x\left[k\right]\right),u\left[k\right]\right),u\left[k+1\right]\right))={\Psi }_3\left(x\left[k\right],u\left[k\right],u\left[k+1\right]\right)\\ .Math.\\ y_1\left[k+n-1\right]=h_1\circ f^{n-1}\left(x\left[k\right],u\left[k\right]\right)={\Psi }_n\left(x\left[k\right],u\left[k\right],.Math.,u\left[k+n-2\right]\right)\end{array}$

[0043] Define:

y.sub.1n[k]=[y.sub.1[k], y.sub.1[k+1], . . . , y.sub.1[k+n−1]]=Ψ(x[k], u.sub.n−1[k]) u.sub.n−1[k]=[u[k],u[k+1], . . . , u[k+n−2]]

[0044] Then the following equation can be deduced:

Ψ(x[k], u.sub.n−1[k])=Y.sub.1n[k]

[0045] Supposed that ∂Ψ/∂x is nonsingular, due to the existence theorem of the implicit function, it can be deduced that the state vector x[k] can be expressed with Y.sub.1n[k] and u.sub.n−1[k].

[0046] Also:

[00003]$\begin{array}{c}x\left[k+n\right]=f\left(x\left[k+n-1\right],u\left[k+n-1\right]\right)\\ =f\left(f\left(x\left[k+n-2\right],u\left[k+n-2\right]\right),u\left[k+n-1\right]\right)\\ =.Math.\\ =f\circ f^{n-1}\left(x\left[k\right],u\left[k\right]\right)\\ =\Phi \left(x\left[k\right],u\left[k\right],u\left[k+1\right],.Math.,u\left[k+n-1\right]\right)\\ =\Omega \left(Y_{1n}\left[k\right],u_{n-1}\left[k\right],u\left[k\right],u\left[k+1\right],.Math.,u\left[k+n-1\right]\right)\\ =\Omega \left(Y_{1n}\left[k\right],u\left[k\right],u\left[k+1\right],.Math.,u\left[k+n-1\right]\right)\end{array}$

[0047] Wherein, Φ(.) and Ω(.) are smoothing functions.

[0048] According to y.sub.1[k+n]=h.sub.1(x[k+n]), the following can be obtained:

[00004]$\begin{array}{c}{y}_1\left[k+n\right]=h_1\left(x\left[k+n\right]\right)\\ =h_1\left(\Omega \left(y_1\left[k\right],y_1\left[k+1\right],.Math.,y_1\left[k+n-1\right],u\left[k\right],u\left[k+1\right],.Math.,u\left[k+n-1\right]\right)\right)\\ =F\left(y_1\left[k\right],y_1\left[k+1\right],.Math.,y_1\left[k+n-1\right],u\left[k\right],u\left[k+1\right],.Math.,u\left[k+n-1\right]\right)\end{array}$

[0049] Wherein, F(.), Φ(.) and Ω(.) are smoothing functions.

[0050] The following can be obtained by the same way:

[00005]$\begin{array}{c}{y}_2\left[k+n\right]=h_2\left(x\left[k+n\right]\right)\\ =h_2\left(\psi \left(y_2\left[k\right],y_2\left[k+1\right],.Math.,y_2\left[k+n-1\right],u\left[k\right],u\left[k+1\right],.Math.,u\left[k+n-1\right]\right)\right)\\ =G\left(y_2\left[k\right],y_2\left[k+1\right],.Math.,y_2\left[k+n-1\right],u\left[k\right],u\left[k+1\right],.Math.,u\left[k+n-1\right]\right)\end{array}$

[00006]$\begin{array}{c}{y}_m\left[k+n\right]=h_m\left(x\left[k+n\right]\right)\\ =h_m\left(P\left(y_m\left[k\right],y_m\left[k+1\right],.Math.,y_m\left[k+n-1\right],u\left[k\right],u\left[k+1\right],.Math.,u\left[k+n-1\right]\right)\right)\\ =H\left(y_m\left[k\right],y_m\left[k+1\right],.Math.,y_m\left[k+n-1\right],u\left[k\right],u\left[k+1\right],.Math.,u\left[k+n-1\right]\right)\end{array}$

[0051] Wherein, ψ(.), G(.), P(.) and H(.) are smoothing functions.

[0052] Therefore, the NARMA equations of the MIMO system are recursively obtained as:

[00007]$\begin{array}{c}{y}_1\left[k+1\right]=F\left(y_1\left[k-n+1\right],y_1\left[k-n+2\right],.Math.,y_1\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k\right]\right)\\ y_2\left[k+1\right]=G\left(y_2\left[k-n+1\right],y_2\left[k-n+2\right],.Math.,y_2\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k\right]\right)\\ .Math.\\ y_m\left[k+1\right]=H\left(y_m\left[k-n+1\right],y_m\left[k-n+2\right],.Math.,y_m\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k\right]\right)\end{array}$

[0053] Wherein, the input variable u[k]=[u.sub.1[k]u.sub.2[k]. . . u.sub.m[k]]ϵR.sup.m, the output values y.sub.1[k]ϵR, y.sub.2[k]ϵR, . . . , y.sub.m[k]ϵR. the functions, F(.), G(.), H(.)ϵC.sup.∞ and k is the time.

[0054] Step A2) The NARMA model is actually an accurate mathematical expression of the nonlinear system in the equilibrium state, but it is not suitable for the actual system due to a large amount of calculation, so the NARMA model needs to be approximated and simplified.

[0055] Make:

[00008]$M_1=\left(y_1\left[k-n+1\right],y_1\left[k-n+2\right],.Math.,y_1\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k-1\right]\right)M_2=\left(y_2\left[k-n+1\right],y_2\left[k-n+2\right],.Math.,y_2\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k-1\right]\right).Math.M_m=\left(y_m\left[k-n+1\right],y_m\left[k-n+2\right],.Math.,y_m\left[k\right],u\left[k-n+1\right],u\left[k-n+2\right],.Math.,u\left[k-1\right]\right)$

[0056] By performing multivariate Taylor expansion for the NARMA equations near the equilibrium point and ignoring Taylor high-order remainders, the NARMA-L2 equations of the system are obtained as:

[00009]$\begin{array}{c}{y}_1\left[k+1\right]=F\left(M_1,u_1\left[k\right],u_2\left[k\right],.Math.,u_m\left[k\right]\right)\\ =F\left(M_1,0,0,.Math.,0\right)\\ {{{+\frac{\partial F}{\partial u_1\left[k\right]}.Math.}_{\left(M_1,0,0,.Math.,0\right)}{u}_1\left[k\right]+\frac{\partial F}{\partial u_2\left[k\right]}.Math.}_{\left(M_1,0,0,.Math.,0\right)}{u}_2\left[k\right]+.Math.+\frac{\partial F}{\partial u_m\left[k\right]}.Math.}_{\left(M_1,0,0,.Math.,0\right)}{u}_m\left[k\right]\end{array}$$\begin{array}{c}{y}_2\left[k+1\right]=G\left(M_2,u_1\left[k\right],u_2\left[k\right],.Math.,u_m\left[k\right]\right)\\ =G\left(M_2,0,0,.Math.,0\right)\\ {{{+\frac{\partial G}{\partial u_1\left[k\right]}.Math.}_{\left(M_2,0,0,.Math.,0\right)}{u}_1\left[k\right]+\frac{\partial G}{\partial u_2\left[k\right]}.Math.}_{\left(M_2,0,0,.Math.,0\right)}{u}_2\left[k\right]+.Math.+\frac{\partial G}{\partial u_m\left[k\right]}.Math.}_{\left(M_2,0,0,.Math.,0\right)}{u}_m\left[k\right]\\ .Math.\end{array}$$\begin{array}{c}{y}_m\left[k+1\right]=G\left(M_m,u_1\left[k\right],u_2\left[k\right],.Math.,u_m\left[k\right]\right)\\ =G\left(M_m,0,0,.Math.,0\right)\\ {{{+\frac{\partial G}{\partial u_1\left[k\right]}.Math.}_{\left(M_m,0,0,.Math.,0\right)}{u}_1\left[k\right]+\frac{\partial G}{\partial u_2\left[k\right]}.Math.}_{\left(M_m,0,0,.Math.,0\right)}{u}_2\left[k\right]+.Math.+\frac{\partial G}{\partial u_m\left[k\right]}.Math.}_{\left(M_m,0,0,.Math.,0\right)}{u}_m\left[k\right]\end{array}$

[0057] Step A3) Convert the NARMA-L2 equations of the MIMO system into a matrix form:

[00010]$\left[\begin{array}{c}{y}_1\left(k+1\right)\\ y_2\left(k+1\right)\\ .Math.\\ y_m\left(k+1\right)\end{array}\right]=\left[\begin{array}{c}{F}_0\\ G_0\\ .Math.\\ H_0\end{array}\right]+\left[\begin{array}{cccc}{F}_1& F_2& .Math.& F_m\\ G_1& G_2& .Math.& G_m\\ .Math.& .Math.& & .Math.\\ H_1& H_2& .Math.& H_m\end{array}\right]\left[\begin{array}{c}{u}_1\left(k\right)\\ u_2\left(k\right)\\ .Math.\\ u_m\left(k\right)\end{array}\right]$

[0058] Wherein,

[00011]${{{{{{{{F_0=F,F_1=\frac{\partial F}{\partial u_1\left(k\right)}.Math.}_{\left(M_1,0,0,.Math.,0\right)},{F_2=\frac{\partial F}{\partial u_2\left(k\right)}.Math.}_{\left(M_1,0,0,.Math.,0\right)},F_m=\frac{\partial F}{\partial u_m\left(k\right)}.Math.}_{\left(M_1,0,0,.Math.,0\right)},G_0=G,G_1=\frac{\partial G}{\partial u_1\left(k\right)}.Math.}_{\left(M_2,0,0,.Math.,0\right)},G_2=\frac{\partial G}{\partial u_2\left(k\right)}.Math.}_{\left(M_2,0,0,.Math.,0\right)},G_m=\frac{\partial G}{\partial u_m\left(k\right)}.Math.}_{\left(M_2,0,0,.Math.,0\right)},H_0=H,H_1=\frac{\partial H}{\partial u_1\left(k\right)}.Math.}_{\left(M_m,0,0,.Math.,0\right)},H_2=\frac{\partial H}{\partial u_2\left(k\right)}.Math.}_{\left(M_m,0,0,.Math.,0\right)},H_m=\frac{\partial H}{\partial u_m\left(k\right)}.Math.}_{\left(M_m,0,0,.Math.,0\right)}$

represent a nonlinear function respectively, and their input parameters are (y.sub.1[k−n+1], y.sub.1[k−n+2], . . . , y.sub.1[k], u[k−1], . . . , u[k−n+1]).

[0059] Therefore, according to the matrix operation, the multivariable control law is obtained as:

[00012]$\left[\begin{array}{c}{u}_1\left(k\right)\\ u_2\left(k\right)\\ .Math.\\ u_m\left(k\right)\end{array}\right]={\left[\begin{array}{cccc}{F}_1& F_2& .Math.& F_m\\ G_1& G_2& .Math.& G_m\\ .Math.& .Math.& & .Math.\\ H_1& H_2& .Math.& H_m\end{array}\right]}^{-1}\left(\left[\begin{array}{c}{y}_1\left(k+1\right)\\ y_2\left(k+1\right)\\ .Math.\\ y_m\left(k+1\right)\end{array}\right]-\left[\begin{array}{c}{F}_0\\ G_0\\ .Math.\\ H_0\end{array}\right]\right)$

[0060] The specific sub-steps of Step B) are as follows:

[0061] Step B1) Take the fuel flow of the turbofan engine and the critical cross-sectional area of the exhaust nozzle as control variables, and take the high-pressure rotating speed and pressure ratio of the engine as controlled variables to design a dual-variable controller. For the dual-variable control system, its control law is:

[00013]$\left[\begin{array}{c}{u}_1\left(k\right)\\ u_2\left(k\right)\end{array}\right]={\left[\begin{array}{cc}{F}_1& F_2\\ G_1& G_2\end{array}\right]}^{-1}\left(\left[\begin{array}{c}{y}_1\left(k+1\right)\\ y_2\left(k+1\right)\end{array}\right]-\left[\begin{array}{c}{F}_0\\ G_0\end{array}\right]\right)$

[0062] Wherein, non-linear functions F0, F1, F2, G0, G1, G2 can be fitted with a neural network, respectively. The six nonlinear functions in the dual-variable control law are identified offline through six neural networks, and the network model adopts the BP-NN network topology.

[0063] The three neural networks approximating the nonlinear functions F0, F1 , F2 all take M.sub.1=(y.sub.1[k−n+1], y.sub.1[k−n+2], . . . , y.sub.1[k], u[k−n+1], u[k−n+2], . . . , u[k−1]) as the input, and their outputs are the approximate values of F.sub.0(y.sub.1[k−n+1], y.sub.1[k−n+2], . . . , y.sub.1[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), F.sub.1(y.sub.1[k−n+1], y.sub.1[k−n+2], . . . , y.sub.1[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), F.sub.2(y.sub.1[k−n+1], y.sub.1[k−n+2], . . . , y.sub.1[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), respectively. The turbofan engine is generally recognized as a second-order system, so n is equal to 3, the input neuron of the neural network is 9, the hidden layer's neuron is 25, the output layer's neuron is 1 and the hidden layer's activation function is the tanh function.

[0064] The three neural networks approximating the nonlinear functions G0, G1, G2 all take M.sub.2=(y.sub.2[k−n+1], y.sub.2[k−n+2], . . . , y.sub.2[k], u[k−n+1], u[k−n+2], . . . , u[k−1]) as the input, and their outputs are the approximate values of G.sub.0(y.sub.2[k−n+1], y.sub.2[k−n+2], . . . , y.sub.2[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), G.sub.1(y.sub.2[k−n+1], y.sub.2[k−n+2], . . . , y.sub.2[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), G.sub.2=(y.sub.2[k−n+1], y.sub.2[k−n+2], . . . , y.sub.2[k], u[k−n+1], u[k−n+2], . . . , u[k−1], u[k]), respectively. Also, n is equal to 3, the input neuron of the neural network is 9, the hidden layer's neuron is 25, the output layer's neuron is 1 and the hidden layer's activation function is the tanh function.

[0065] The engine model is used to generate a large amount of ground point input and output data. Such data are used as the input of each neural network. The model output of the NARMA-L2 equation is calculated through the output of each neural network, and the error is backpropagated to train each neural network.

[0066] The performance index of the model is:

[00014]$$\begin{gathered} J=\frac{1}{2}\left(e_1^2\left[k\right]+e_2^2\left[k\right]\right) \\ {e}_1=y_1\left[k\right]-y{m}_1\left[k\right] \\ {e}_2=y_2\left[k\right]-y{m}_2\left[k\right] \end{gathered}$$

[0067] Wherein, ym.sub.1 and ym.sub.2 respectively represent the model outputs calculated according to the NARMA-L2 equations, and y.sub.1 and y.sub.2 respectively represent the actual high-pressure rotating speed and pressure ratio of the engine.

[0068] The offline identification structure of the NARMA-L2 model is shown in FIG. 1, where the bolded variables are all matrices, namely u [k]=[u.sub.1[k]u.sub.2[k]].sup.T, y[k]=[y.sub.1[k]y.sub.2[k]].sup.T, ym[k]=[ym[k]ym.sub.2[k]].sup.T, e=[e.sub.1 e.sub.2 ].sup.T. The NARMA-L2 model identification error curve is shown in FIG. 2. As the number of training increases, the neural network identification is basically completed. The theoretical output values calculated with the NARMA-L2 equations are basically consistent with the actual output values.

[0069] Step B2) According to the NARMA-L2 dual-variable control law, construct a matrix with the outputs of the nonlinear functions already identified by the neural networks, establish a closed-loop controller for the high-pressure rotating speed and pressure ratio of the engine, and get the values of the control variables online. Moreover, establish a quadratic performance index according to the error between the actual output and the command output, and use the backpropagation algorithm to optimize the neural network parameters online.

[0070] The performance index is:

[00015]$$\begin{gathered} J=\frac{1}{2}\left(e_{c1}^2\left[k\right]+e_{c2}^2\left[k\right]\right) \\ {e}_{c1}=y_1\left[k\right]-r_1\left[k\right] \\ {e}_{c2}=y_2\left[k\right]-r_2\left[k\right] \end{gathered}$$

[0071] Wherein, r.sub.1 and r.sub.2 represent the engine's high-pressure rotating speed and pressure ratio commands, respectively; and y.sub.1 and y.sub.2 represent the engine's actual high-pressure rotating speed and pressure ratio, respectively.

[0072] The built control structure is shown in FIG. 3, where the bolded variables are all matrices, namely u[k]=[u.sub.1[k]u.sub.2[k]].sup.T, y[k]=[y.sub.1[k]y.sub.2[k]].sup.T, r[k]=[k]r.sub.1[k]r.sub.2[k]].sup.T, e.sub.c=[e .sub.c1 e.sub.c2].sup.T.

[0073] Step B3) After the dual-variable controller is designed, the dual-variable closed-loop control system of the turbofan engine is established and its tracking performance is verified at design points and non-design points and within the full flight envelope.

[0074] The simulation results under the state of the design point (H=0 m and Ma=0) are shown in FIG. 4 and FIG. 5. It can be seen from the figures that the high-pressure rotating speed and pressure ratio under this state can track the command step signal stably without steady-state error; the overshoot of the high-pressure rotating speed is 6.32% and the adjustment time is less than 2 s; the overshoot of the pressure ratio is 5.53% and the adjustment time is less than 1.4 s, so it has good stable and dynamic state performance.

[0075] The simulation results at the non-design point (H=4,000 m and Ma=0.8) are shown in FIG. 6 and FIG. 7. It can be seen from the figures that the high-pressure rotating speed and pressure ratio under this state can track the change of the command signal stably, basically without steady-state error; the overshoot of the high-pressure rotating speed is 7.95% and the adjustment time is less than 2.5 s; and the overshoot of the pressure ratio is 6.34% and the adjustment time is less than 2 s, so it has good dynamic characteristics.

[0076] The changes of the altitude and Mach number during the take-off, climb, descent, acceleration and deceleration process of the aircraft are shown in FIG. 8, and the output tracking effects of the high-pressure rotating speed and pressure ratio are shown in FIG. 9 and FIG. 10. It can be seen from the figures that in various processes, the high-pressure rotating speed and pressure ratio can basically track the change of the command stably, and the system responds quickly, so the model can well meet the needs of the dual-variable control, has good control quality and can adapt to the changes of various working conditions of the engine.

[0077] It should be pointed out that what mentioned above is only the embodiment of the present invention, but the protection scope of the present invention is not limited to this. Any changes and replacements that any technician familiar with this technical field can easily make within the technical scope disclosed by the present invention should fall into the protection scope of the present invention. Therefore, the protection scope of the present invention shall be the protection scope of the following claims.