WO2019237320A1 - 一种基于模型预测的航空发动机在线优化及多变量控制设计方法 - Google Patents

一种基于模型预测的航空发动机在线优化及多变量控制设计方法 Download PDF

Info

Publication number
WO2019237320A1
WO2019237320A1 PCT/CN2018/091442 CN2018091442W WO2019237320A1 WO 2019237320 A1 WO2019237320 A1 WO 2019237320A1 CN 2018091442 W CN2018091442 W CN 2018091442W WO 2019237320 A1 WO2019237320 A1 WO 2019237320A1
Authority
WO
WIPO (PCT)
Prior art keywords
engine
control
state
model
aero
Prior art date
Application number
PCT/CN2018/091442
Other languages
English (en)
French (fr)
Inventor
杜宪
马艳华
孙希明
王智民
Original Assignee
大连理工大学
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 大连理工大学 filed Critical 大连理工大学
Priority to US16/462,519 priority Critical patent/US11454177B2/en
Priority to PCT/CN2018/091442 priority patent/WO2019237320A1/zh
Publication of WO2019237320A1 publication Critical patent/WO2019237320A1/zh

Links

Images

Classifications

    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F02—COMBUSTION ENGINES; HOT-GAS OR COMBUSTION-PRODUCT ENGINE PLANTS
    • F02C—GAS-TURBINE PLANTS; AIR INTAKES FOR JET-PROPULSION PLANTS; CONTROLLING FUEL SUPPLY IN AIR-BREATHING JET-PROPULSION PLANTS
    • F02C9/00—Controlling gas-turbine plants; Controlling fuel supply in air- breathing jet-propulsion plants
    • F02C9/26—Control of fuel supply
    • F02C9/28—Regulating systems responsive to plant or ambient parameters, e.g. temperature, pressure, rotor speed
    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B17/00—Systems involving the use of models or simulators of said systems
    • G05B17/02—Systems involving the use of models or simulators of said systems electric
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B64—AIRCRAFT; AVIATION; COSMONAUTICS
    • B64D—EQUIPMENT FOR FITTING IN OR TO AIRCRAFT; FLIGHT SUITS; PARACHUTES; ARRANGEMENT OR MOUNTING OF POWER PLANTS OR PROPULSION TRANSMISSIONS IN AIRCRAFT
    • B64D31/00—Power plant control systems; Arrangement of power plant control systems in aircraft
    • B64D31/02—Initiating means
    • B64D31/06—Initiating means actuated automatically
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F05—INDEXING SCHEMES RELATING TO ENGINES OR PUMPS IN VARIOUS SUBCLASSES OF CLASSES F01-F04
    • F05D—INDEXING SCHEME FOR ASPECTS RELATING TO NON-POSITIVE-DISPLACEMENT MACHINES OR ENGINES, GAS-TURBINES OR JET-PROPULSION PLANTS
    • F05D2220/00—Application
    • F05D2220/30—Application in turbines
    • F05D2220/32—Application in turbines in gas turbines
    • F05D2220/323—Application in turbines in gas turbines for aircraft propulsion, e.g. jet engines
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F05—INDEXING SCHEMES RELATING TO ENGINES OR PUMPS IN VARIOUS SUBCLASSES OF CLASSES F01-F04
    • F05D—INDEXING SCHEME FOR ASPECTS RELATING TO NON-POSITIVE-DISPLACEMENT MACHINES OR ENGINES, GAS-TURBINES OR JET-PROPULSION PLANTS
    • F05D2260/00—Function
    • F05D2260/80—Diagnostics
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F05—INDEXING SCHEMES RELATING TO ENGINES OR PUMPS IN VARIOUS SUBCLASSES OF CLASSES F01-F04
    • F05D—INDEXING SCHEME FOR ASPECTS RELATING TO NON-POSITIVE-DISPLACEMENT MACHINES OR ENGINES, GAS-TURBINES OR JET-PROPULSION PLANTS
    • F05D2260/00—Function
    • F05D2260/81—Modelling or simulation
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F05—INDEXING SCHEMES RELATING TO ENGINES OR PUMPS IN VARIOUS SUBCLASSES OF CLASSES F01-F04
    • F05D—INDEXING SCHEME FOR ASPECTS RELATING TO NON-POSITIVE-DISPLACEMENT MACHINES OR ENGINES, GAS-TURBINES OR JET-PROPULSION PLANTS
    • F05D2270/00—Control
    • F05D2270/01—Purpose of the control system
    • F05D2270/20—Purpose of the control system to optimize the performance of a machine
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F05—INDEXING SCHEMES RELATING TO ENGINES OR PUMPS IN VARIOUS SUBCLASSES OF CLASSES F01-F04
    • F05D—INDEXING SCHEME FOR ASPECTS RELATING TO NON-POSITIVE-DISPLACEMENT MACHINES OR ENGINES, GAS-TURBINES OR JET-PROPULSION PLANTS
    • F05D2270/00—Control
    • F05D2270/40—Type of control system
    • F05D2270/44—Type of control system active, predictive, or anticipative
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Definitions

  • the invention provides an aero-engine online optimization and multi-variable control design method based on model prediction, and belongs to the field of aerospace propulsion system control and simulation technology.
  • Model predictive control as a model-based multi-variable control algorithm, can not only achieve effective control of multiple variables in aero engine systems, but also deal with various constraints that need to be addressed in the process of aero engine control. It can fully realize the potential of aero engines. Improve the performance of the entire aviation propulsion system.
  • model predictive control has many excellent characteristics, it also has the disadvantage that it is difficult to balance controller performance and control system real-time performance.
  • the present invention continuously establishes a linearized model of the engine for each control cycle, and designs a multivariable model predictive controller based on the consideration of the control system actuator to achieve real-time online optimization and reduce the amount of model predictive control calculations. While improving the performance of the control system.
  • the present invention proposes a model-based aviation Engine online optimization and multivariable control design method.
  • the control system structure is mainly composed of two parts.
  • the first part is the prediction model acquisition layer.
  • the actual working state and external environmental parameters of each control cycle of the aeroengine are Based on the benchmark, continuously build a small deviation linear model of the engine near different steady-state points, and provide the model parameters to the controller;
  • the second part is the control law optimization decision-making layer, which consists of the model prediction controller and external output feedback. Constrained linear optimization problem and closed-loop feedback to determine the control output of the controller; the specific steps are as follows:
  • a small deviation linear model of the aero engine is established in real time and continuously based on the current actual working state of the aero engine and the external environment; the actual input of the aero engine includes the fuel flow in the main combustion chamber, the fuel flow in the afterburner, and the tail nozzle cutoff.
  • Engine controllable parameters such as area, fan stator vane angle, compressor stator vane angle, etc .
  • external environmental parameters include environmental parameters such as flying height, flying Mach number, etc .
  • small deviation linear models of aero-engines are obtained by the method of pumping power. Spatial expression
  • step S2 Normalize the small deviation linear model of the aero-engine obtained in step S1, that is, convert the absolute increments of the input variables, state variables, and output variables into relative increments about the steady-state points by linear transformation; normalize The subsequent model is discretized to make the parameters of the model convenient for digital control.
  • the sampling period for discretization is the same as the control period.
  • the model to predict the multivariable controller includes: the normalized and discretized linear model parameters of the aeroengine, the discretized model of the actuator, the constraints of the control system, and the control structure of the control system
  • the actuator and the aero engine are two subsystems connected in series. Therefore, each control cycle combines the aero engine and the actuator in series to form a state space expression of the overall system, and expands to an expanded state space expression.
  • the objective function of minimizing the index is constructed by the factors such as the difference between the predicted output and the expected output, and the magnitude of the change in the input amount, and then the actual input and output limit of the control system is used as the optimization constraint. Constrained linear optimization problem to obtain the next control output;
  • step S4 Apply the controller output obtained in step S3 to the controlled system, compare it with the actual output of the aero-engine obtained by the sensor, obtain the control error, and use it as a model predictive controller input to guide the controller's decision on the next controller output .
  • the corresponding steady state value of the amount of fuel supplied to the engine at the previous time by the S1.1 control system is used as the steady state point parameter of the small deviation linear model, and the corresponding dynamic equilibrium equation is established based on the engine's rotor as the energy storage component:
  • J is a diagonal matrix composed of the rotor inertia, and the dimension of the diagonal matrix depends on the number of rotors of the engine.
  • I the derivative vector of each rotor of the engine with respect to time, and ⁇ M is the residual torque vector of each rotor;
  • ⁇ M ⁇ M (n, u in , u out )
  • n is a vector composed of the rotational speed of each rotor of the engine
  • u in is a vector composed of a controlled variable
  • u out is a vector composed of external environmental variables to be considered
  • Jac is the Jacobian matrix of ⁇ M versus n and u in , ⁇ n and ⁇ u in are column vectors, the elements of ⁇ n are the absolute increase of the speed of each rotor of the engine relative to the steady-state speed, and the elements of ⁇ u in are the controls Absolute increment of the quantity relative to the steady-state control quantity;
  • ⁇ y is a column vector, and its elements are the absolute increments of the parameters of the engine except the speed relative to the steady state point;
  • the subscript s represents a dimensionless vector
  • the subscript m represents the state quantity, input quantity, and output quantity of the engine at a steady state point
  • a a is a diagonal matrix composed of the negative derivative of the time constant of the actuator
  • B a is a diagonal matrix composed of the derivative of the time constant of the actuator
  • a ad and B ad are the system matrix and input matrix after discretization of the state space expression of the actuator, respectively;
  • x t (k + 1) A t x t (k) + B t u t (k)
  • S3.3 further forms a new augmented state vector from the state vector and the input vector, and uses the increment of the input vector as the new input vector to obtain the augmented state space expression:
  • ⁇ u (k) is the increment vector of the control amount in the k-th control cycle
  • new state space expression is written in a compact form:
  • x t ′ (k + 1) A t ′ x t ′ (k) + B t ′ ⁇ u (k)
  • the online optimization objective can be described as a simplified form of the objective function as follows:
  • r is the controller input used as the reference sequence.
  • online rolling optimization can be achieved by solving the quadratic programming problem.
  • the sequence of incremental control variables obtained by the optimization the first element of the sequence is taken and summed with the control variables of the previous control cycle as the next one. The actual amount of control during the control cycle.
  • An aeroengine on-line optimization and multivariable control design method based on model prediction provided by the present invention For the prediction model acquisition, no traditional interpolation or fitting method is used. Instead, based on the existing aero-engine nonlinear component-level model, the steady-state model under current input conditions is obtained in real time as an important part of the current prediction model through the power extraction method. Compared with the traditional interpolation and fitting methods, it can be closer to the real characteristic curve of the engine. Compared with directly using the engine nonlinear model to solve the nonlinear optimization problem, it can greatly reduce the calculation time and improve the real-time performance of online optimization.
  • the model predictive multivariable controller designed by the present invention can directly introduce the saturation limit of the actuator and the safety boundaries and constraints of the aero-engine-related parameters to achieve online optimization without the need to plan in advance the acceleration and deceleration plans that ensure the safe operation of the system
  • the curve and limit protection greatly reduce the complexity of the control system design, fully realize the potential of the engine, and help improve the overall performance of the flight propulsion system.
  • the aeroengine multivariable control system structure based on model prediction designed by the present invention has a large closed loop of the actual output of the engine, which can reflect the true response of the engine to the output of the controller, so as to adjust the model to predict the expected command of the controller and compensate
  • the influence of model mismatch and external disturbance on the system improves the control accuracy of the control system.
  • Figure 1 is the curve of the rotation speed of the high-pressure rotor.
  • Figure 2 is the curve of turbine drop pressure.
  • Figure 3 is the temperature change curve before the turbine.
  • Figure 4 is the change curve of fuel flow.
  • FIG. 5 is a change curve of the cross-sectional area of the tail nozzle.
  • Figure 6 is a curve of fan surge margin.
  • Fig. 7 is a curve of the compressor surge margin.
  • FIG. 8 is a curve of engine thrust.
  • Figure 9 is a comparison curve of high-speed rotor speed under different operating conditions.
  • Figure 10 is a comparison curve of the turbine drop pressure ratio under different operating conditions.
  • Figure 11 is a comparison curve of the temperature before the turbine under different operating conditions.
  • Figure 12 is a comparison curve of fuel flow under different operating conditions.
  • Figure 13 is a comparison curve of the cross-sectional area of the tail nozzle under different working conditions.
  • Figure 14 is a comparison curve of fan surge margin under different operating conditions.
  • Figure 15 is the change curve of compressor surge margin under different operating conditions.
  • Figure 16 is a comparison of engine thrust under different operating conditions.
  • FIG. 17 is a structural diagram of an aero-engine online optimization and multivariable control system based on model prediction.
  • This embodiment is a design method of aero-engine online optimization and multi-variable control based on model prediction.
  • the detailed design steps are as follows:
  • Step 1 Establish a small deviation gain model of the aeroengine with the current aeroengine actual input and external environmental parameters as steady-state points; first, obtain the steady-state point for calculating the linearized small deviation model.
  • a small deviation model of the current time is calculated by using the input amount of the engine at the last sampling time and the environmental parameters as the steady state points.
  • the method of calculating the small deviation model of the engine near the steady state point is as follows.
  • the thermal inertia of the component is much smaller than that of the rotor. Therefore, the heat transfer process is ignored in the modeling, and the engine rotor is used as the energy storage component.
  • the engine in the embodiment is a typical two-rotor engine, so its dynamic characteristics depend on the torque balance equation of the two rotors, namely:
  • J H and J L are the moments of inertia of the high and low pressure rotors of the engine
  • ⁇ M H and ⁇ M L are the residual torques of the high and low pressure rotors, respectively.
  • the residual torque of the high and low pressure rotors can be expressed as:
  • M TH and M TL are high and low pressure turbine torques respectively
  • M CH and M CL are high and low pressure compressor torques, respectively.
  • M fr, H , M fr, L are the friction torques of the high and low pressure rotors , respectively. Since the values are much smaller than the other four terms, they are ignored during the calculation of the linearization model.
  • ⁇ M H and ⁇ M L are non-linear functions of high and low pressure rotor speed, main fuel flow, cross section of tail nozzle and pressure and temperature of working environment, so it can be expanded into the following linearized form according to Taylor's formula. :
  • ⁇ n H , ⁇ n L , ⁇ q mf , ⁇ A 8 , ⁇ p 2 , and ⁇ T 2 are absolute increments of each variable relative to the steady state point.
  • the pressure, temperature, and engine thrust along the cross section of the engine are functions of speed, input, and environmental interference, and can be expressed as:
  • x x (n H , n L , q mf , A 8 , p 2 , T 2 ) (10) x is an arbitrary parameter of the engine other than the speed. After ignoring the change of the working environment in the current control cycle, it can be further Can be linearized as:
  • the present invention constructs a linear model of the engine by using a linear state-space expression.
  • the state vector is composed of two elements.
  • the rotation speeds ⁇ n H and ⁇ n L of the two rotors of the dual-rotor engine are selected as state variables, and the state equation can be written as:
  • the coefficients in the model expressions can be obtained through the following work.
  • the method we use to obtain the partial derivative of the steady-state point is based on the engine's non-linear calculation program.
  • the balance calculation will include the following five states:
  • the column principal element elimination method can be used to program, and each pending partial derivative can be solved to complete the linearization of the model near the steady-state operating point of the aero engine.
  • Step 2 The aeroengine small deviation gain model obtained in step 1 is normalized and discretized. Because the units used for each element in the input vector, output vector, and even the state vector are not the same, and the magnitudes of different variables are quite different, the magnitude of the elements in the parameter matrix of the engine linearization model is greatly different, and it is easy to form a ill-conditioned matrix. It is very easy to affect the accuracy of the calculation, so the model needs to be normalized.
  • the elements in the input vector, state vector and output vector in equations (12) and (13) are absolute increments.
  • the magnitude difference between the absolute increments is large, which is easy to cause in the calculation process. Large error. Therefore, the present invention adopts an engine linear normalization model obtained by converting an absolute increment into a relative increment for a steady-state operating point.
  • the absolute increment needs to be dimensionlessly converted to a relative increment.
  • the expression of the relative increment is The formulas are as follows:
  • x i , u i , y i are the elements in the corresponding vector
  • x mi , u mi , y mi are the values of the parameters of the steady state point
  • x si , u si , y si are the relative increments of the elements in each vector.
  • the model predictive control algorithm is a computer control algorithm, it is necessary to discretize the model.
  • the sampling period of discretization should be consistent with the control period of the control system.
  • the general form of the discrete state space expression is as follows:
  • a ed , B ed , C ed , and D ed are parameter matrices after discretization.
  • Step 3 Design the model to predict the multivariable controller.
  • the main work of the model predictive controller includes constructing an overall predictive model and performing a limited time-domain linear optimization with constraints based on the expected instructions of the input controller to determine the output of the next controller.
  • a ad and B ad are discretized parameter matrices of the actuator.
  • the state space expressions of the actuator and the engine can be combined in series to form the prediction model of the model prediction controller as the following form:
  • x t (k + 1) A t x t (k) + B t u t (k) (31)
  • x t [x 1 , x 2 , ⁇ n Hs , ⁇ n Ls ] T
  • u t [q mfs , A 8 ] T
  • y t [ ⁇ n Hs , ⁇ n Ls , ⁇ T 4s , ⁇ Ts ] T ;
  • the state vector and the input vector are constituted as a new augmented state vector, and the increment of the input vector is used as the new input vector, so we obtain the augmented state space expression of the aero engine:
  • x t ′ (k + 1) A t ′ x t ′ (k) + B t ′ ⁇ u (k) (35)
  • r is the controller input used as the reference sequence. After expansion and merging and removing terms that have no effect on J, the objective function function is finally reduced to the form of the following quadratic form:
  • U are the input matrix upper limit and input matrix lower limit, respectively.
  • Step 4 Apply the controller output obtained in Step 3 to the controlled system, obtain the actual output of the aero engine through the sensor, and compare it with the expected value to obtain the control error.
  • the control error will be used as a new reference sequence in the optimization process of the model predictive controller, which will affect the form of the objective function of the model prediction. Objectively, it will compensate the effect of the mismatch of the prediction model and external disturbances on the control performance to ensure the control system. Stability and control accuracy.
  • the output of the engine related to the constraints will be directly fed back into the model predictive controller to make reasonable adjustments to the upper and lower limits of the output variable matrix. Ensure that the input and output of the engine do not exceed the limit, and ensure the stable and safe operation of the system.
  • the effectiveness of the method of the present invention is verified by a simulation experiment of a group of engines starting from acceleration, turning on afterburner, turning off afterburner, and entering a deceleration process.
  • the simulation running time is 70s, the engine high-pressure rotor speed and the turbine drop pressure ratio are controlled variables.
  • the front turbine temperature of 1650K and 1700K are used as engine safety constraints, and the application of the same is required.
  • the control signal to the control actuator meets the physical constraints of the amplitude and rate of change.
  • the pressure ratio is 11.
  • the fuel flow rate of the afterburner will be increased to 5000kg / h at the 20th moment, and the core machine control variables will be kept unchanged in the afterburner state.
  • the afterburner chamber will be closed. Deceleration begins at the 50s until a predetermined new steady-state point is reached.
  • the speed of the new steady-state high-pressure rotor is 12000 r / min, and the turbine drop pressure ratio is 6.5.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Engineering & Computer Science (AREA)
  • Combustion & Propulsion (AREA)
  • Chemical & Material Sciences (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Geometry (AREA)
  • Data Mining & Analysis (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Automation & Control Theory (AREA)
  • Mechanical Engineering (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Feedback Control In General (AREA)

Abstract

一种基于模型预测控制的航空发动机在线优化及多变量控制设计方法，控制系统包括两部分：第一部分是预测模型获取层，以航空发动机每个控制周期的实际工作状态和外部环境参数为基准，连续建立不同稳态点附近的发动机小偏离线性模型，并将模型参数提供给控制器；第二个部分是控制规律决策层，由模型预测控制器和外部输出反馈构成的闭环结构，模型预测控制器根据当前状态下的发动机模型，控制指令和相关约束限制，并确定控制器下一时刻输出，外部输出反馈则将航空发动机实际输出引入到控制器对未来控制量的决策中，以补偿模型失配和外部扰动的影响，实现航空发动机多个变量的控制与在线优化，可提高控制系统的控制精度。

Description

一种基于模型预测的航空发动机在线优化及多变量控制设计方法 技术领域
本发明提供了一种基于模型预测的航空发动机在线优化及多变量控制设计方法，属于航空航天推进系统控制与仿真技术领域。
背景技术
随着航空发动机技术的不断进步，现代航空发动机的结构日益复杂，工作范围也不断扩大，航空发动机控制系统的要求也越来越高，传统的单输入/单输出的控制系统已经难以满足控制需求。因此，选取更多的控制变量以实现对发动机多个参数进行控制的多变量控制系统成为提高发动机控制系统性能的重要手段。模型预测控制作为一种基于模型的多变量控制算法，不仅可以实现对航空发动机系统多个变量的有效控制，而且可以处理航空发动机控制过程需要处理各种约束问题，可以充分发挥航空发动机的潜力，提高整个航空推进系统性能。尽管模型预测控制具有许多优秀的特性，但也存在着控制器性能与控制系统实时性难以兼顾的缺点。本发明基于发动机的实际输入量，每个控制周期连续建立发动机线性化模型，并在考虑控制系统执行机构的基础上设计多变量模型预测控制器，实现实时在线优化，在减少模型预测控制计算量的同时提高了控制系统的性能。
发明内容
为了保证航空发动机在整个飞行包线内，考虑执行机构自身特性、物理约束以及工作限制条件下，实现多个控制量下的稳态控制和过渡态控制，本发明提出一种基于模型预测的航空发动机在线优化及多变量控制设计方法。
本发明的技术方案：
一种基于模型预测的航空发动机在线优化及多变量控制设计方法，控制系统结构主要由两部分组成，第一部分是预测模型获取层，以航空发动机每个控制周期的实际工作状态和外部环境参数为基准，连续建立不同稳态点附近的发动机小偏离线性模型，并将模型参数提供给控制器；第二个部分是控制规律优化决策层，由模型预测控制器和外部输出反馈构成，通过求解带约束的线性优化问题与闭环反馈作用，确定控制器的控制输出；具体步骤如下：
S1.实时连续地建立以航空发动机当前实际工作状态和外部环境为稳态点的航空发动机小偏离线性模型；航空发动机的实际输入包括主燃烧室燃油流量、加力燃烧室燃油流量、尾喷口截面积、风扇静子导流叶片角度、压气机静子导流叶片角度等发动机可控参数；外部环境参数包括飞行高度、飞行马赫数等环境参数；航空发动机小偏离线性模型通过抽功法得到，以状态空间形式表达；
S2.将步骤S1获得的航空发动机小偏离线性模型进行归一化处理，即将输入变量、状态变量、输出变量的绝对增量通过线性变换转变为关于稳态点的相对增量；将归一化之后的模型进行离散化处理，使得模型的参数便于数字控制使用，离散化的采样周期与控制周期相同；
S3.设计模型预测多变量控制器；控制器中存储的数据包括：航空发动机归一化和离散化的线性模型参数，执行机构的离散化模型，控制系统的约束条件；在控制系统的控制结构中，执行机构和航空发动机为串联方式连接的两个子系统，因此每个控制周期将航空发动机和执行机构通过串联方式合并为一个整体系统的状态空间表达式，并拓展为增广状态空间表达式；对于在线优化部分，以预测输出与期望输出的差值、输入量变化幅度等因素构建指标最小化的目标函数，再以控制系统的实际输入输出量限制为优化约束条件，通过在线求解上述带约束的线性优化问题获得下一步的控制输出；
S4.将步骤S3得到的控制器输出作用于被控系统，与传感器获得航空发动机的实际输出进行对比，得到控制误差，作为模型预测控制器的输入，指导控制器对下一步控制器输出的决策。
建立以航空发动机实际输入和外部环境为稳态点的航空发动机小偏离线性模型的步骤如下：
S1.1控制系统上一时刻供给发动机的燃油量相应的稳态值作为小偏离线性模型的稳态点参数，并根据发动机的转子作为储能部件建立相应动力学平衡方程：
Figure PCTCN2018091442-appb-000001
其中，J为由转子转动惯量组成的对角矩阵，对角矩阵的维数取决于发动机的转子个数，
Figure PCTCN2018091442-appb-000002
为发动机各转子对时间的导数向量，ΔM为各转子剩余力矩向量；
S1.2根据发动机原理，各转子的剩余力矩表示为如下形式：
ΔM＝ΔM(n,u in,u out)
其中，n为发动机各转子转速组成的向量，u in为控制量组成的向量，u out为需要考虑的外部环境变量组成的向量；
S1.3由于飞行器惯性较发动机惯性大得多，所以建立动态方程时不考虑环境参数的变化，将转子剩余力矩的表达式近似线性化为：
Figure PCTCN2018091442-appb-000003
其中，Jac为ΔM对n和u in的雅可比矩阵，Δn和Δu in均为列向量，Δn的元素为发动机各转子转速相对于稳态点转速的绝对增量，Δu in的元素为各控制量相对于稳态点控制量的绝对增量；
S1.4将
Figure PCTCN2018091442-appb-000004
J一项移动到等号右端，并将Jac拆分为矩阵A和矩阵B，得到如下形式的转子动力学方程：
Figure PCTCN2018091442-appb-000005
S1.5在只考虑转子惯性的条件下，发动机各截面的压力、温度以及发动机推力都是转速、输入量和环境参数的函数，表示为：y＝y(n,u in,u out)，y表示为除转速以外的发动机任意参数，在忽略当前控制周期内环境参数干扰后，航空发动机参数表达式线性化为：
Δy＝CΔn+DΔu in
其中Δy为列向量，其元素为发动机除转速外的各参数相对于稳态点的绝对增量；
S1.6矩阵A,B,C,D的元素均由相应偏导数确定，相应偏导数通过小扰动所带来的稳态参数的变化值获得；将步骤S1.4中的航空发动机转子动力学方程与步骤S1.5中的发动机参数表达式组合起来，得到航空发动机小偏离线性模型的状态空间表达式：
Figure PCTCN2018091442-appb-000006
对获得的航空发动机小偏离线性模型进行归一化与离散化处理的步骤如下：
S2.1将步骤S1中得到的状态空间表达式中输入向量、状态向量和输出向量中的绝对增量变换为相对增量，无量纲的相对增量的表达形式表示成如下形式：
Figure PCTCN2018091442-appb-000007
其中，下角标s表示无量纲化后的向量，下角标m表示发动机在稳态点的状态量、输入量和输出量；
S2.2将发动机在稳态点附近的状态空间表达式转换为如下形式：
Figure PCTCN2018091442-appb-000008
y＝C eΔn+D eΔu in
其中，A e＝N -1AN，B e＝N -1BU，C e＝Y -1CN，D e＝Y -1DU，N,Y,U分别为由稳态点状态量，输出量，输入量组成的对角矩阵；
S2.3由于模型预测控制算法是一种计算机控制算法，因此需要对模型进行离散化处理， 离散化处理的采样周期应与控制系统的控制周期保持一致，离散后的状态空间表达式的一般形式表示如下：
Figure PCTCN2018091442-appb-000009
设计模型预测多变量控制器及在线优化的步骤如下：
S3.1考虑执行机构的影响，并将执行机构近似为一阶惯性环节，执行机构的状态空间表达式写成如下形式：
Figure PCTCN2018091442-appb-000010
其中，A a为执行机构时间常数负导数组成的对角矩阵，B a为执行机构时间常数导数组成的对角矩阵，按照与控制周期相同的采样时间可以将其离散化为：
Figure PCTCN2018091442-appb-000011
其中，A ad和B ad分别为执行机构状态空间表达式离散化后的系统矩阵和输入矩阵；
S3.2将执行机构与发动机模型按照串联方式整合为统一的状态空间表达式为：
x t(k+1)＝A tx t(k)+B tu t(k)
y t(k)＝C tx t(k)
式中，x t＝[x,Δn] T，u t＝Δu in T，y t＝[Δn,Δy] T；
Figure PCTCN2018091442-appb-000012
S3.3进一步将状态向量和输入向量构成新的增广状态向量，并将输入向量的增量作为新的输入向量，得到增广状态空间表达式：
Figure PCTCN2018091442-appb-000013
y t(k)＝C tx t(k)
其中，Δu(k)为第k个控制周期内，控制量的增量向量，将新的状态空间表达式写成紧凑的形式：
x t′(k+1)＝A t′x t′(k)+B t′Δu(k)
y t(k)＝C tx t(k)
根据航空发动机的状态空间表达式，令预测时域为p和控制时域为q，得到状态序列与输出序列的预测公式的一般形式如下：
Figure PCTCN2018091442-appb-000014
Figure PCTCN2018091442-appb-000015
式中，
Figure PCTCN2018091442-appb-000016
Figure PCTCN2018091442-appb-000017
Figure PCTCN2018091442-appb-000018
Figure PCTCN2018091442-appb-000019
Figure PCTCN2018091442-appb-000020
根据系统的状态空间表达式、状态序列预测公式以及输出序列预测公式，在线优化目标可以描述成如下的目标函数的简化形式：
Figure PCTCN2018091442-appb-000021
式中r为作为参考序列的控制器输入，经过展开，合并剔除对J无影响的项之后，目标函数最终可化简为如下二次型的形式：
Figure PCTCN2018091442-appb-000022
对于输入约束，可表述为如下紧凑形式：
Figure PCTCN2018091442-appb-000023
式中，
Figure PCTCN2018091442-appb-000024
Figure PCTCN2018091442-appb-000025
和 U分别为输入量上限和输入矩阵下限；
对于输出约束，表述为如下形式：
Figure PCTCN2018091442-appb-000026
Figure PCTCN2018091442-appb-000027
式中，
Figure PCTCN2018091442-appb-000028
和 Y分别是输出量上限和下限；
根据目标函数和约束条件，可以通过求解二次规划问题实现在线滚动优化，根据优化获得的控制量增量序列，采取序列的首元素，与上一控制周期的控制量求和后，作为下一个控制周期内的实际控制量。
本发明的有益效果：
(1)本发明提出的一种基于模型预测的航空发动机在线优化及多变量控制设计方法。对于预测模型的获取，并没有采用传统的插值法或拟合法。而是根据现有的航空发动机非线性部件级模型，通过抽功法，实时获取当前状输入条件下的稳态模型作为当前预测模型的重要组成部分。相较于传统插值法和拟合法可以更贴近发动机的真实特性曲线，而与直接使用 发动机非线性模型求解非线性优化问题相比，可以大大减少计算时间，提高在线优化的实时性。
(2)本发明在构建预测模型的过程中，注意到了发动机线性化模型归一化的必要性，并考虑了执行机构动态特性对控制系统影响。在航空发动机模型中，输入向量，输出向量和状态向量的单位不同，且向量内的部分元素间数量级相差较大，因此会造成矩阵元素在数量级上相差较大，甚至产生病态矩阵影响进一步计算的精度。为此，本发明在生成线性化模型的过程中引入了归一化的环节以提高运算精度。同时，虽然执行机构的时间常数往往小于航空发动机的时间常数，但对于基于模型的控制算法，执行机构的引入依然会对控制系统的控制效果造成影响。在本发明中，将执行机构和发动机这两个串联的系统合并为一个整体，虽然提高了系统的阶次，但是取得了更好的控制效果。
(3)本发明设计的模型预测多变量控制器可以直接引入执行机构的饱和限制和航空发动机相关参数的安全边界和约束条件，实现在线优化，不需要提前离线规划保证系统安全运行的加减速计划曲线和限制保护，大大降低了控制系统设计的复杂性，充分发挥了发动机的潜力，有助于提高飞行推进系统的整体性能。
(4)本发明设计的基于模型预测的航空发动机多变量控制系统结构拥有一个发动机实际输出大闭环，可以反映出发动机对控制器输出的真实响应，以此调整模型预测控制器的期望指令，补偿模型失配与外部扰动对系统的影响，提高控制系统的控制精度。
附图说明
图1是高压转子转速变化曲线。
图2是涡轮落压变化曲线。
图3是涡轮前温度变化曲线。
图4是燃油流量变化曲线。
图5是尾喷管截面积变化曲线。
图6是风扇喘振裕度变化曲线。
图7是压气机喘振裕度变化曲线。
图8是发动机推力变化曲线。
图9是不同工况下高压转子转速对比曲线。
图10是不同工况下涡轮落压比对比曲线。
图11是不同工况下涡轮前温度对比曲线。
图12是不同工况下燃油流量对比曲线。
图13是不同工况下尾喷管截面积对比曲线。
图14是不同工况下风扇喘振裕度对比曲线。
图15是不同工况下压气机喘振裕度变化曲线。
图16是不同工况下发动机推力对比曲线。
图17是基于模型预测的航空发动机在线优化及多变量控制系统结构图。
具体实施方式
以下结合附图和技术方案，进一步说明本发明的具体实施方式。
本实施例是一种基于模型预测的航空发动机在线优化及多变量控制设计方法。具体详细设计步骤如下：
步骤1：建立以当前航空发动机实际输入和外部环境参数为稳态点的航空发动机小偏离增益模型；首先获得计算线性化小偏离模型的稳态点。本发明中以发动机上一采样时刻的输入量与环境参数近似作为稳态点来计算当前时刻的小偏离模型。计算发动机在稳态点附近的小偏离模型的方法如下所述。
在航空发动机工作的过程中，部件热惯性相对转子惯性小得多，因此在建模时忽略传热过程，将发动机转子作为储能部件。实施例中的发动机为一种典型的双转子发动机，因此其动态特性取决于两个转子的力矩平衡方程，即：
Figure PCTCN2018091442-appb-000029
式中，J H，J L分别为发动机高、低压转子的转动惯量，ΔM H，ΔM L分别为高低压转子的剩余力矩。
而根据发动机原理，高、低压转子的剩余力矩可以表示为：
Figure PCTCN2018091442-appb-000030
式中，M TH，M TL分别为高、低压涡轮转矩，M CH，M CL分别为高、低压压气机转矩。M fr,H，M fr,L分别为高、低压转子的摩擦力矩，由于数值远远小于其他四项，所以在计算线性化模型的过程中忽略不计。由式(2)可知，ΔM H，ΔM L分别是高、低压转子转速，主燃油流量，尾喷管截面积以及工作环境压力温度的非线性函数，因此可以根据泰勒公式展开成如下线性化形式：
Figure PCTCN2018091442-appb-000031
Figure PCTCN2018091442-appb-000032
式中，Δn H，Δn L，Δq mf，ΔA 8，Δp 2，ΔT 2分别为各变量相对于稳态点的绝对增量。
因为飞行器惯性比发动机的惯性大很多，所以在建立动态方程式时，可以假设在当前控制周期内飞行条件保持不变，也就是暂时不考虑p 2和T 2的影响。结合上面列出的式子，可以得到高、低压转子的运动方程如下：
Figure PCTCN2018091442-appb-000033
式中，
Figure PCTCN2018091442-appb-000034
Figure PCTCN2018091442-appb-000035
Figure PCTCN2018091442-appb-000036
Figure PCTCN2018091442-appb-000037
在只考虑转子惯性的条件下，发动机沿各截面的压力、温度以及发动机推力都是转速、输入量和环境干扰量的函数，可以表示为：
x＝x(n H,n L,q mf,A 8,p 2,T 2)       (10)x表示为除转速以外的发动机任意参数，在忽略当前控制周期内工作环境的变化后，可以进一步可以线性化为：
Figure PCTCN2018091442-appb-000038
于是可以得到包括涡轮前温度和涡轮落压比在内的发动机其他参数的线性化动态方程。
在得到航空发动机状态变量和输出变量的近似线性函数表达式之后，本发明通过线性状态空间表达式来构建发动机的线性模型。
在双轴发动机中，当只考虑两个独立转子部件为储能元件时，状态向量由两个元素组成。本发明选择双转子发动机的两个转子的转速Δn H和Δn L作为状态变量，状态方程可以写成：
Figure PCTCN2018091442-appb-000039
实施例中选取高、低压转子转速，涡轮前温度，涡轮落压比作为输出变量时，可以得到如下输出方程：
Figure PCTCN2018091442-appb-000040
式中，
Figure PCTCN2018091442-appb-000041
Figure PCTCN2018091442-appb-000042
Figure PCTCN2018091442-appb-000043
Figure PCTCN2018091442-appb-000044
在获得了双轴发动机用状态空间表达式描述的线性方程后，可以通过下面的工作求取模型表达式中的系数，也就是各项偏导数。我们采用的获取稳态点偏导数的方法是以发动机的非线性计算程序为基础，通过抽取发动机高压转子、低压转子上的功率，以及改变燃油流量和尾喷管截面面积，通过平衡计算获取高低压转子转速，涡轮前温度和压气机总压比的变化量，并根据变化量求出偏导数。其中平衡计算将包括以下5种状态的计算：
状态一：根据给定的稳态工作点，获取发动机在计算点的稳态参数值，记为n H0，n L0，T 40，π T0。
状态二：其他条件不变，在低压轴上抽取小量功率，计算低压轴抽取小量功率时发动机转速及其他参数值，记为n H1，n L1，T 41，π T1。
状态三：其他条件不变，在高压轴上抽取小量功率，计算高压轴抽取功率时发动机各参数值记为n H2，n L2，T 42，π T2。
状态四：其他条件不变，燃油流量相对于状态一有小量变化时，计算燃油流量偏离稳态值时发动机各参数值n H3，n L3，T 43，π T3。
状态五：其他条件不变，尾喷口截面面积相对于状态一有小量变化，计算尾喷口截面面积偏离稳态值时发动机各参数数值n H4，n L4，T 44，π T4。
在整理计算状态二、三、四、五相对于状态一的偏离值之后我们可以得到以下三组方程组：
Figure PCTCN2018091442-appb-000045
Figure PCTCN2018091442-appb-000046
Figure PCTCN2018091442-appb-000047
得到如上所示的三组线性方程组后，通过编程利用列主元消去法，求解各待定的偏导数即可完成对航空发动机稳态工作点附近模型的线性化。
步骤2：将步骤1获得的航空发动机小偏离增益模型进行归一化、离散化处理。由于在 输入向量，输出向量乃至状态向量中的各个元素所使用的单位并不一致，且不同变量的数量级相差较大，造成发动机线性化模型的参数矩阵内的元素数量级相差较大，容易形成病态矩阵，极易影响计算的精度，所以需要对模型进行归一化处理。
式子(12)，(13)中的输入向量，状态向量和输出向量中的元素均为绝对增量，在实际应用中，绝对增量之间的数量级相差较大，在计算过程中容易造成较大误差。因此本发明采用将绝对增量转化为对于稳态工作点的相对增量的方式来获取的发动机线性归一化模型，首先需要将绝对增量无量纲化为相对增量，相对增量的表达式分别如下：
Figure PCTCN2018091442-appb-000048
式中
Figure PCTCN2018091442-appb-000049
x i，u i，y i为对应向量中的元素，x mi，u mi，y mi分别为稳态点各参数的值，
Figure PCTCN2018091442-appb-000050
x si，u si，y si分别为各向量中元素的相对增量。
整理之后的状态空间表达式可以写成如下形式：
Figure PCTCN2018091442-appb-000051
Figure PCTCN2018091442-appb-000052
式中：
Figure PCTCN2018091442-appb-000053
Figure PCTCN2018091442-appb-000054
由于模型预测控制算法是一种借助于计算机控制算法，因此需要对模型进行离散化处理。离散化处理的采样周期应与控制系统的控制周期保持一致。离散后的状态空间表达式的一般形式表示如下：
Figure PCTCN2018091442-appb-000055
式中，A ed，B ed，C ed，D ed为离散化后的参数矩阵。
步骤3：设计模型预测多变量控制器。模型预测控制器的主要工作包括构建整体的预测模型以及根据输入控制器的期望指令进行有限时域的带约束线性优化，以确定下一步控制器的输出。
由于本发明设计的控制系统中考虑到了执行机构的影响，并将执行机构近似一阶惯性环节，所以需要将执行机构也考虑到预测模型之中。首先将执行机构的近似模型写成如下所示的状态空间表达式：
Figure PCTCN2018091442-appb-000056
Figure PCTCN2018091442-appb-000057
将式(27)，(28)离散化后得到：
Figure PCTCN2018091442-appb-000058
Figure PCTCN2018091442-appb-000059
式中，A ad和B ad为执行机构离散化后的参数矩阵。
由于执行机构与发动机在系统中是串联的，所以可以将执行机构和发动机的状态空间表达式按照串联方式合并作为模型预测控制器的预测模型写成如下形式：
x t(k+1)＝A tx t(k)+B tu t(k)   (31)
y t(k)＝C tx t(k)    (32)
式中，x t＝[x 1,x 2,Δn Hs,Δn Ls] T，u t＝[q mfs,A 8] T，y t＝[Δn Hs,Δn Ls,ΔT 4s,Δπ Ts] T；
Figure PCTCN2018091442-appb-000060
对于已经得到的预测模型，将状态向量和输入向量构成为新的增广状态向量，并将输入向量的增量作为新的输入向量，于是我们得到航空发动机的增广状态空间表达式：
Figure PCTCN2018091442-appb-000061
y t(k)＝C tx t(k)   (34)
其中Δu(k)为第k个控制周期内，控制量的增量向量，将新的状态空间表达式写成紧凑的形式：
x t′(k+1)＝A t′x t′(k)+B t′Δu(k)   (35)
y t(k)＝C tx t(k)   (36)
根据航空发动机的状态空间表达式，令预测时域为10和控制时域位3，可以得到状态序列与输出序列的预测公式的一般形式如下：
Figure PCTCN2018091442-appb-000062
Figure PCTCN2018091442-appb-000063
式中，
Figure PCTCN2018091442-appb-000064
Figure PCTCN2018091442-appb-000065
Figure PCTCN2018091442-appb-000066
Figure PCTCN2018091442-appb-000067
Figure PCTCN2018091442-appb-000068
根据系统的状态空间表达式，状态序列预测公式以及输出序列预测公式后。我们可以写出如下的目标函数的简化形式：
Figure PCTCN2018091442-appb-000069
式中r为作为参考序列的控制器输入，经过展开，合并剔除对J无影响的项之后，目标函数函数最终被化简为如下二次型的形式：
Figure PCTCN2018091442-appb-000070
对于输入约束，我们可以表述为如下紧凑形式：
Figure PCTCN2018091442-appb-000071
Figure PCTCN2018091442-appb-000072
式中，
Figure PCTCN2018091442-appb-000073
Figure PCTCN2018091442-appb-000074
和 U分别位输入矩阵上限和输入矩阵下限。
对于输出矩阵的约束，则可以表示为如下形式：
Figure PCTCN2018091442-appb-000075
Figure PCTCN2018091442-appb-000076
式中，
Figure PCTCN2018091442-appb-000077
和 Y分别是输出矩阵的上限和下限。
根据目标函数和约束条件，我们可以通过求解二次规划问题实现滚动优化。在取得优化后的控制输出增量序列之后，我们取序列的首元素作为下一个控制周期内的实际控制输出。
步骤4.将步骤3得到的控制器输出作用于被控系统，通过传感器获得航空发动机的实际输出，将其与期望值进行对比，得到控制误差。控制误差将作为模型预测控制器优化过程中新的参考序列，对模型预测的目标函数的形式产生影响，客观上对预测模型失配和外部扰动对控制性能的影响产生补偿作用，保证科控制系统的稳定性和控制精度。而发动机与约束条件相关的输出量将直接被反馈到模型预测控制器当中，以便对输出变量矩阵的上限和下限做出合理的调整。保证发动机的输入输出不超限，确保系统稳定安全运行。
为了进一步说明本实施例中模型预测多变量控制器的效果，通过一组发动机从开始加速，打开加力，关闭加力到到进入减速过程的仿真实验，来验证本发明中方法的有效性。
仿真运行时间为70s，以发动机高压转子转速和涡轮落压比为被控变量，在高度0km，飞行速度为0Ma的条件下，分别以涡轮前温度1650K和1700K作为发动机安全约束条件，同时要求施加到控制执行机构的控制信号满足变化幅值和变化速率的物理约束条件。在0-5s内保持高压转子转速为11622r/min，涡轮落压比为6；在5s后开始加速直至不开加力的最大 状态，该状态下，高压转子的转速为14500r/min，涡轮落压比为11。达到稳定后在第20s时刻开始增加加力燃烧室的燃油流量至5000kg/h，并要求核心机控制变量在加力状态下保持不变；第40s关闭加力燃烧室；待其稳定后，在第50s开始减速，直至达到预定的新稳态点，新稳态高压转子转速为12000r/min，涡轮落压比为6.5。
通过观察图1和图2，可以看到发动机的控制变量和控制目标符合的很好。在发动机打开加力后也可以很好的将各控制变量稳定在设定的水平；通过观察图3，可发现控制器可以很好的将涡轮前温度限制在安全的范围内；通过观察图4和图5，可以发现控制器能够满足执行机构的幅值限制和速率限制的要求；通过图6和图7可以发现，由于被控变量满足了执行机构的限制要求，使得发动机的风扇喘振裕度和压气机喘振裕度满足安全需求；通过观察图8，可以发现，虽然存在限制条件的差异，控制器对发动机实际推力控制依然十分有效。
为了验证控制系统在飞行包线内的适用性，在飞行包线内分别取(2km,0.2Ma)，(5km,0.5Ma)，(8km,0.8Ma)三个点，分别进行发动机从开始加速，打开加力，关闭加力到到进入减速过程的仿真实验。根据发动机的工作环境发生变化和部件实际物理特性，将涡轮前温度的限制条件调整为1800K。控制执行机构的控制信号满足变化幅值和变化速率的物理约束条件不变。70s的仿真过程安排如下：
在0-5s内保持高压转子转速为11622r/min，涡轮落压比为7；在5s后开始加速直至不开加力的最大状态，该状态下，高压转子的转速为14500r/min，涡轮落压比为11。达到稳定后在第20s时刻开始增加加力燃烧室的燃油流量至5000kg/h，并要求核心机控制变量在加力状态下保持不变；第40s关闭加力燃烧室；待其稳定后，在第50s开始减速，直至达到预定的新稳态点，新稳态高压转子转速为12000r/min，涡轮落压比为7.5。
通过观察图9和图10，可以看到发动机在不同工况下稳态的控制变量和控制目标符合的很好，动态过程略有差异，这与发动机不同工作状态下控制输入的变化有关。在发动机打开加力后不同工况下同样可以很好的将各控制变量稳定在设定的水平；通过观察图11，可发现控制器均可以很好的将涡轮前温度限制在安全的范围内；通过观察图12和图13，可以发现不同工作状态下发动机控制输入并不相同，因此造成了图9和图10中动态过程的差异；通过图14和图15可以发现，发动机的风扇喘振裕度和压气机喘振裕度满足安全需求，但是可以明显的发现不同工况下的安全裕度并不相同；通过观察图16，可以发现，虽然高压转子转速和涡轮落压比相同，发动机在不同的工况下，实际推力会有所不同。

Claims (5)

  1. 一种基于模型预测的航空发动机在线优化及多变量控制设计方法，控制系统结构主要由两部分组成，第一部分是预测模型获取层，以航空发动机每个控制周期的实际工作状态和外部环境参数为基准，连续建立不同稳态点附近的发动机小偏离线性模型，并将模型参数提供给控制器；第二个部分是控制规律优化决策层，由模型预测控制器和外部输出反馈构成，通过求解带约束的线性优化问题与闭环反馈作用，确定控制器的控制输出；其特征在于，具体步骤如下：
    S1.实时连续地建立以航空发动机当前实际工作状态和外部环境为稳态点的航空发动机小偏离线性模型；航空发动机的实际输入包括主燃烧室燃油流量、加力燃烧室燃油流量、尾喷口截面积、风扇静子导流叶片角度、压气机静子导流叶片角度发动机可控参数；外部环境参数包括飞行高度、飞行马赫数环境参数；航空发动机小偏离线性模型通过抽功法得到，以状态空间形式表达；
    S2.将步骤S1获得的航空发动机小偏离线性模型进行归一化处理，即将输入变量、状态变量、输出变量的绝对增量通过线性变换转变为关于稳态点的相对增量；将归一化之后的模型进行离散化处理，使得模型的参数便于数字控制使用，离散化的采样周期与控制周期相同；
    S3.设计模型预测多变量控制器；控制器中存储的数据包括：航空发动机归一化和离散化的线性模型参数，执行机构的离散化模型，控制系统的约束条件；在控制系统的控制结构中，执行机构和航空发动机为串联方式连接的两个子系统，每个控制周期将航空发动机和执行机构通过串联方式合并为一个整体系统的状态空间表达式，并拓展为增广状态空间表达式；对于在线优化部分，以预测输出与期望输出的差值、输入量变化幅度等因素构建指标最小化的目标函数，再以控制系统的实际输入输出量限制为优化约束条件，通过在线求解上述带约束的线性优化问题获得下一步的控制输出；
    S4.将步骤S3得到的控制器输出作用于被控系统，与传感器获得航空发动机的实际输出进行对比，得到控制误差，作为模型预测控制器的输入，指导控制器对下一步控制器输出的决策。
  2. 根据权利要求1所述的基于模型预测的航空发动机在线优化及多变量控制设计方法，其特征在于，建立以航空发动机实际输入和外部环境为稳态点的航空发动机小偏离线性模型的步骤如下：
    S1.1控制系统上一时刻供给发动机的燃油量相应的稳态值作为小偏离线性模型的稳态点参数，并根据发动机的转子作为储能部件建立相应动力学平衡方程：
    Figure PCTCN2018091442-appb-100001
    其中，J为由转子转动惯量组成的对角矩阵，对角矩阵的维数取决于发动机的转子个数，
    Figure PCTCN2018091442-appb-100002
    为发动机各转子对时间的导数向量，ΔM为各转子剩余力矩向量；
    S1.2根据发动机原理，各转子的剩余力矩表示为如下形式：
    ΔM＝ΔM(n,u in,u out)
    其中，n为发动机各转子转速组成的向量，u in为控制量组成的向量，u out为需要考虑的外部环境变量组成的向量；
    S1.3由于飞行器惯性较发动机惯性大得多，所以建立动态方程时不考虑环境参数的变化，将转子剩余力矩的表达式近似线性化为：
    Figure PCTCN2018091442-appb-100003
    其中，Jac为ΔM对n和u in的雅可比矩阵，Δn和Δu in均为列向量，Δn的元素为发动机各转子转速相对于稳态点转速的绝对增量，Δu in的元素为各控制量相对于稳态点控制量的绝对增量；
    S1.4将
    Figure PCTCN2018091442-appb-100004
    J一项移动到等号右端，并将Jac拆分为矩阵A和矩阵B，得到如下形式的转子动力学方程：
    Figure PCTCN2018091442-appb-100005
    S1.5在只考虑转子惯性的条件下，发动机各截面的压力、温度以及发动机推力都是转速、输入量和环境参数的函数，表示为：y＝y(n,u in,u out)，y表示为除转速以外的发动机任意参数，在忽略当前控制周期内环境参数干扰后，航空发动机参数表达式线性化为：
    Δy＝CΔn+DΔu in
    其中Δy为列向量，其元素为发动机除转速外的各参数相对于稳态点的绝对增量；
    S1.6矩阵A,B,C,D的元素均由相应偏导数确定，相应偏导数通过小扰动所带来的稳态参数的变化值获得；将步骤S1.4中的航空发动机转子动力学方程与步骤S1.5中的发动机参数表达式组合起来，得到航空发动机小偏离线性模型的状态空间表达式：
    Figure PCTCN2018091442-appb-100006
  3. 根据权利要求1或2所述的基于模型预测的航空发动机在线优化及多变量控制设计方法，其特征在于，对获得的航空发动机小偏离线性模型进行归一化与离散化处理的步骤如下：
    S2.1将步骤S1中得到的状态空间表达式中输入向量、状态向量和输出向量中的绝对增量变换为相对增量，无量纲的相对增量的表达形式表示成如下形式：
    Figure PCTCN2018091442-appb-100007
    其中，下角标s表示无量纲化后的向量，下角标m表示发动机在稳态点的状态量、输入量和输出量；
    S2.2将发动机在稳态点附近的状态空间表达式转换为如下形式：
    Figure PCTCN2018091442-appb-100008
    y＝C eΔn+D eΔu in
    其中，A e＝N -1AN，B e＝N -1BU，C e＝Y -1CN，D e＝Y -1DU，N,Y,U分别为由稳态点状态量，输出量，输入量组成的对角矩阵；
    S2.3对航空发动机小偏离线性模型进行离散化处理，离散化处理的采样周期应与控制系统的控制周期保持一致，离散后的状态空间表达式的一般形式表示如下：
    Figure PCTCN2018091442-appb-100009
  4. 根据权利要求1或2所述的基于模型预测的航空发动机在线优化及多变量控制设计方法，其特征在于，设计模型预测多变量控制器及在线优化的步骤如下：
    S3.1考虑执行机构的影响，并将执行机构近似为一阶惯性环节，执行机构的状态空间表达式写成如下形式：
    Figure PCTCN2018091442-appb-100010
    其中，A a为执行机构时间常数负导数组成的对角矩阵，B a为执行机构时间常数导数组成的对角矩阵，按照与控制周期相同的采样时间将其离散化为：
    Figure PCTCN2018091442-appb-100011
    其中，A ad和B ad分别为执行机构状态空间表达式离散化后的系统矩阵和输入矩阵；
    S3.2将执行机构与发动机模型按照串联方式整合为统一的状态空间表达式为：
    x t(k+1)＝A tx t(k)+B tu t(k)
    y t(k)＝C tx t(k)
    式中，x t＝[x,Δn] T，u t＝Δu in T，y t＝[Δn,Δy] T；
    Figure PCTCN2018091442-appb-100012
    C t＝[D edC ad,C ed]；
    S3.3进一步将状态向量和输入向量构成新的增广状态向量，并将输入向量的增量作为新的输入向量，得到增广状态空间表达式：
    Figure PCTCN2018091442-appb-100013
    y t(k)＝C tx t(k)
    其中，Δu(k)为第k个控制周期内，控制量的增量向量，将新的状态空间表达式写成紧凑的形式：
    x′ t(k+1)＝A′ tx′ t(k)+B′ tΔu(k)
    y t(k)＝C tx t(k)
    根据航空发动机的状态空间表达式，令预测时域为p和控制时域为q，得到状态序列与输出序列的预测公式的一般形式如下：
    Figure PCTCN2018091442-appb-100014
    Figure PCTCN2018091442-appb-100015
    式中，
    Figure PCTCN2018091442-appb-100016
    Figure PCTCN2018091442-appb-100017
    Figure PCTCN2018091442-appb-100018
    Figure PCTCN2018091442-appb-100019
    Figure PCTCN2018091442-appb-100020
    根据系统的状态空间表达式、状态序列预测公式以及输出序列预测公式，在线优化目标可以描述成如下的目标函数的简化形式：
    Figure PCTCN2018091442-appb-100021
    式中r为作为参考序列的控制器输入，经过展开，合并剔除对J无影响的项之后，目标函数最终可化简为如下二次型的形式：
    Figure PCTCN2018091442-appb-100022
    对于输入约束，表述为如下紧凑形式：
    Figure PCTCN2018091442-appb-100023
    Figure PCTCN2018091442-appb-100024
    式中，
    Figure PCTCN2018091442-appb-100025
    和
    Figure PCTCN2018091442-appb-100026
    分别为输入量上限和输入矩阵下限；
    对于输出约束，表述为如下形式：
    Figure PCTCN2018091442-appb-100027
    Figure PCTCN2018091442-appb-100028
    式中，
    Figure PCTCN2018091442-appb-100029
    和
    Figure PCTCN2018091442-appb-100030
    分别是输出量上限和下限；
    根据目标函数和约束条件，通过求解二次规划问题实现在线滚动优化，根据优化获得的控制量增量序列，采取序列的首元素，与上一控制周期的控制量求和后，作为下一个控制周期内的实际控制量。
  5. 根据权利要求3所述的基于模型预测的航空发动机在线优化及多变量控制设计方法，其特征在于，设计模型预测多变量控制器及在线优化的步骤如下：
    S3.1考虑执行机构的影响，并将执行机构近似为一阶惯性环节，执行机构的状态空间表达式写成如下形式：
    Figure PCTCN2018091442-appb-100031
    其中，A a为执行机构时间常数负导数组成的对角矩阵，B a为执行机构时间常数导数组成的对角矩阵，按照与控制周期相同的采样时间将其离散化为：
    Figure PCTCN2018091442-appb-100032
    其中，A ad和B ad分别为执行机构状态空间表达式离散化后的系统矩阵和输入矩阵；
    S3.2将执行机构与发动机模型按照串联方式整合为统一的状态空间表达式为：
    x t(k+1)＝A tx t(k)+B tu t(k)
    y t(k)＝C tx t(k)
    式中，x t＝[x,Δn] T，u t＝Δu in T，y t＝[Δn,Δy] T；
    Figure PCTCN2018091442-appb-100033
    C t＝[D edC ad,C ed]；
    S3.3进一步将状态向量和输入向量构成新的增广状态向量，并将输入向量的增量作为新的输入向量，得到增广状态空间表达式：
    Figure PCTCN2018091442-appb-100034
    y t(k)＝C tx t(k)
    其中，Δu(k)为第k个控制周期内，控制量的增量向量，将新的状态空间表达式写成紧凑的形式：
    x′ t(k+1)＝A′ tx′ t(k)+B′ tΔu(k)
    y t(k)＝C tx t(k)
    根据航空发动机的状态空间表达式，令预测时域为p和控制时域为q，得到状态序列与输出序列的预测公式的一般形式如下：
    Figure PCTCN2018091442-appb-100035
    Figure PCTCN2018091442-appb-100036
    式中，
    Figure PCTCN2018091442-appb-100037
    Figure PCTCN2018091442-appb-100038
    Figure PCTCN2018091442-appb-100039
    Figure PCTCN2018091442-appb-100040
    Figure PCTCN2018091442-appb-100041
    根据系统的状态空间表达式、状态序列预测公式以及输出序列预测公式，在线优化目标可以描述成如下的目标函数的简化形式：
    Figure PCTCN2018091442-appb-100042
    式中r为作为参考序列的控制器输入，经过展开，合并剔除对J无影响的项之后，目标函数最终可化简为如下二次型的形式：
    Figure PCTCN2018091442-appb-100043
    对于输入约束，表述为如下紧凑形式：
    Figure PCTCN2018091442-appb-100044
    Figure PCTCN2018091442-appb-100045
    式中，
    Figure PCTCN2018091442-appb-100046
    和
    Figure PCTCN2018091442-appb-100047
    分别为输入量上限和输入矩阵下限；
    对于输出约束，表述为如下形式：
    Figure PCTCN2018091442-appb-100048
    Figure PCTCN2018091442-appb-100049
    式中，
    Figure PCTCN2018091442-appb-100050
    和
    Figure PCTCN2018091442-appb-100051
    分别是输出量上限和下限；
    根据目标函数和约束条件，通过求解二次规划问题实现在线滚动优化，根据优化获得的控制量增量序列，采取序列的首元素，与上一控制周期的控制量求和后，作为下一个控制周期内的实际控制量。
PCT/CN2018/091442 2018-06-15 2018-06-15 一种基于模型预测的航空发动机在线优化及多变量控制设计方法 WO2019237320A1 (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
US16/462,519 US11454177B2 (en) 2018-06-15 2018-06-15 Method of aero-engine on-line optimization and multivariable control based on model prediction
PCT/CN2018/091442 WO2019237320A1 (zh) 2018-06-15 2018-06-15 一种基于模型预测的航空发动机在线优化及多变量控制设计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2018/091442 WO2019237320A1 (zh) 2018-06-15 2018-06-15 一种基于模型预测的航空发动机在线优化及多变量控制设计方法

Publications (1)

Publication Number Publication Date
WO2019237320A1 true WO2019237320A1 (zh) 2019-12-19

Family

ID=68839793

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2018/091442 WO2019237320A1 (zh) 2018-06-15 2018-06-15 一种基于模型预测的航空发动机在线优化及多变量控制设计方法

Country Status (2)

Country Link
US (1) US11454177B2 (zh)
WO (1) WO2019237320A1 (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111931366A (zh) * 2020-07-31 2020-11-13 中国航发贵阳发动机设计研究所 一种航空发动机可调喷管反馈钢索行程的计算方法
CN112282957A (zh) * 2020-11-11 2021-01-29 西华大学 一种二冲程航空活塞发动机性能优化的热管理系统与方法
CN112668162A (zh) * 2020-12-17 2021-04-16 江苏航空职业技术学院 一种基于惯性滑模的航空发动机建模方法
CN114781153A (zh) * 2022-04-18 2022-07-22 北京航空航天大学 一种整机变维度仿真性能仿真流程控制方法
CN114781153B (zh) * 2022-04-18 2024-06-04 北京航空航天大学 一种整机变维度仿真性能仿真流程控制方法

Families Citing this family (22)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2019144386A1 (zh) * 2018-01-26 2019-08-01 大连理工大学 一种航空发动机过渡态关键性能参数预测方法
CN111651913B (zh) * 2020-05-13 2023-03-24 上海交通大学 一种汽车发动机性能预测及辅助标定方法及系统
CN111608808A (zh) * 2020-06-15 2020-09-01 西北工业大学 输入受限航空发动机增益调度容错控制器
CN111852662A (zh) * 2020-06-15 2020-10-30 西北工业大学 航空发动机最大推力状态容错二自由度h∞控制器
CN113836645B (zh) * 2020-10-27 2024-02-02 深圳三零三防务科技有限公司 一种运载火箭的在线飞行程序重构及轨道高度保持控制方法
CN113447271A (zh) * 2020-11-30 2021-09-28 中国人民解放军火箭军工程大学 一种基于修正偏导数的航空发动机气路在线监测方法
CN112464478B (zh) * 2020-11-30 2023-06-30 中国长江电力股份有限公司 一种水轮机调速系统的控制规律优化方法及装置
WO2022126472A1 (zh) * 2020-12-17 2022-06-23 大连理工大学 一种多几何参数可调的进/排/发一体化航空推进系统建模方法
CN112836290B (zh) * 2021-01-08 2023-08-11 聊城大学 一种基于系统输出的自由活塞直线电机匹配优化方法
CN112861258A (zh) * 2021-01-14 2021-05-28 西北工业大学 一种航空发动机最大推力控制优化方法
CN112879165A (zh) * 2021-01-14 2021-06-01 西北工业大学 考虑气路部件故障的航空发动机加速过程最优控制方法
US11859539B2 (en) * 2021-02-01 2024-01-02 General Electric Company Aircraft propulsion system with inter-turbine burner
CN113268000B (zh) * 2021-05-20 2022-08-09 大连理工大学 一种航空发动机多模型预测控制的软切换方法
CN113279997B (zh) * 2021-06-04 2022-04-12 大连理工大学 一种基于控制器模糊切换的航空发动机喘振主动控制系统
CN113495486B (zh) * 2021-08-06 2023-11-24 南京工业大学 一种结构热试验基于扩展状态观测器的模型预测控制方法
CN113759727B (zh) * 2021-09-30 2023-08-29 中国航发控制系统研究所 航空发动机多变量控制器的综合优化设计方法
CN114036449B (zh) * 2021-10-15 2023-07-11 北京航空航天大学 热声稳定性预测方法和装置
CN114115276A (zh) * 2021-11-26 2022-03-01 江苏科技大学 一种基于在线分组优化模型预测的船舶动力定位控制方法
CN114415506B (zh) * 2022-01-07 2023-08-04 大连理工大学 航空发动机双模跟踪预测控制系统设计方法
CN115434802B (zh) * 2022-09-15 2024-05-07 西安交通大学 氨-氢双燃料航空转子发动机多目标优化控制策略及系统
CN116049977B (zh) * 2022-12-26 2024-04-12 西南科技大学 一种航空发动机燃烧室的参数多目标优化方法
CN116776654B (zh) * 2023-08-24 2023-11-03 北京航空航天大学 一种航空发动机流道方案合理性评估方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6459963B1 (en) * 2000-07-31 2002-10-01 General Electric Company Methods and apparatus for trimming engine control systems
CN103984242A (zh) * 2014-05-19 2014-08-13 上海交通大学 基于模型预测控制的分层预测控制系统及其控制方法
CN105114189A (zh) * 2015-06-09 2015-12-02 吉林大学 基于fpga实现的电子节气门模型预测控制系统
CN106647268A (zh) * 2016-12-21 2017-05-10 东南大学 基于模型预测控制的mgt‑cchp分层最优控制系统
CN107193212A (zh) * 2017-06-26 2017-09-22 南京航空航天大学 基于新型灰狼优化算法的航空发动机非线性预测控制方法
CN108762089A (zh) * 2018-06-15 2018-11-06 大连理工大学 一种基于模型预测的航空发动机在线优化及多变量控制设计方法

Family Cites Families (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7219040B2 (en) * 2002-11-05 2007-05-15 General Electric Company Method and system for model based control of heavy duty gas turbine
EP1938499A4 (en) * 2005-09-19 2011-06-29 Univ State Cleveland CONTROLS, OBSERVERS AND APPLICATIONS THEREOF
US8849542B2 (en) * 2012-06-29 2014-09-30 United Technologies Corporation Real time linearization of a component-level gas turbine engine model for model-based control
CN102855349B (zh) * 2012-08-06 2015-07-01 南京航空航天大学 航空发动机气路故障诊断的快速原型设计方法及平台
CN104102769B (zh) * 2014-06-27 2017-02-22 南京航空航天大学 基于人工智能的涡轴发动机实时部件级模型建立方法
US9732625B2 (en) * 2015-01-19 2017-08-15 United Technologies Corporation System and method for controlling a gas turbine engine
CN105652665B (zh) * 2016-03-03 2018-03-20 东南大学 一种微型燃气轮机冷热电三联供系统的协调控制方法
CN106886151B (zh) * 2017-04-17 2019-08-13 大连理工大学 一种航空发动机多工况下约束预测控制器的设计及调度方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6459963B1 (en) * 2000-07-31 2002-10-01 General Electric Company Methods and apparatus for trimming engine control systems
CN103984242A (zh) * 2014-05-19 2014-08-13 上海交通大学 基于模型预测控制的分层预测控制系统及其控制方法
CN105114189A (zh) * 2015-06-09 2015-12-02 吉林大学 基于fpga实现的电子节气门模型预测控制系统
CN106647268A (zh) * 2016-12-21 2017-05-10 东南大学 基于模型预测控制的mgt‑cchp分层最优控制系统
CN107193212A (zh) * 2017-06-26 2017-09-22 南京航空航天大学 基于新型灰狼优化算法的航空发动机非线性预测控制方法
CN108762089A (zh) * 2018-06-15 2018-11-06 大连理工大学 一种基于模型预测的航空发动机在线优化及多变量控制设计方法

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111931366A (zh) * 2020-07-31 2020-11-13 中国航发贵阳发动机设计研究所 一种航空发动机可调喷管反馈钢索行程的计算方法
CN111931366B (zh) * 2020-07-31 2024-05-24 中国航发贵阳发动机设计研究所 一种航空发动机可调喷管反馈钢索行程的计算方法
CN112282957A (zh) * 2020-11-11 2021-01-29 西华大学 一种二冲程航空活塞发动机性能优化的热管理系统与方法
CN112282957B (zh) * 2020-11-11 2022-08-19 西华大学 一种二冲程航空活塞发动机性能优化的热管理系统与方法
CN112668162A (zh) * 2020-12-17 2021-04-16 江苏航空职业技术学院 一种基于惯性滑模的航空发动机建模方法
CN114781153A (zh) * 2022-04-18 2022-07-22 北京航空航天大学 一种整机变维度仿真性能仿真流程控制方法
CN114781153B (zh) * 2022-04-18 2024-06-04 北京航空航天大学 一种整机变维度仿真性能仿真流程控制方法

Also Published As

Publication number Publication date
US11454177B2 (en) 2022-09-27
US20190383221A1 (en) 2019-12-19

Similar Documents

Publication Publication Date Title
WO2019237320A1 (zh) 一种基于模型预测的航空发动机在线优化及多变量控制设计方法
CN108762089B (zh) 一种基于模型预测的航空发动机在线优化及多变量控制设计方法
CN111006843B (zh) 一种暂冲式超声速风洞的连续变速压方法
CN112286047B (zh) 基于神经网络的narma-l2多变量控制方法
CN111731490B (zh) 变旋翼转速直升机综合控制方法、装置
CN109896026A (zh) 变旋翼转速直升机-涡轴发动机综合控制方法及装置
CN109446605B (zh) 涡轴发动机非线性动态逆控制方法及装置
CN112483261B (zh) 一种航空发动机抗加力扰动方法
CN108196443A (zh) 变循环发动机的非线性预测控制设计方法
Wang et al. Overshoot-free acceleration of aero-engines: An energy-based switching control method
Liu et al. Model reference adaptive control for aero-engine based on system equilibrium manifold expansion model
CN112594069B (zh) 基于功率需求预测的涡轴发动机抗扰控制方法及装置
Yu et al. A new method for the design of optimal control in the transient state of a gas turbine engine
Smith et al. Optimizing aircraft performance with adaptive, integrated flight/propulsion control
Lv et al. Extended-Kalman-filter-based equilibrium manifold expansion observer for ramjet nonlinear control
CN110985216A (zh) 一种含在线修正的航空发动机智能多变量控制方法
Jia et al. Multi-variable anti-disturbance controller with state-dependent switching law for adaptive cycle engine
CN114637210A (zh) 航空发动机多回路综合控制方法及系统
CN112668162A (zh) 一种基于惯性滑模的航空发动机建模方法
CN112363411A (zh) 一种航空发动机动态矩阵控制器的设计方法
Liu et al. Design and implementation for the state time-delay and input saturation compensator of gas turbine aero-engine control system
CN113669167B (zh) 涡扇发动机加力燃油执行机构自适应控制方法及装置
Yu et al. Direct Thrust Control of Aeroengines Based on Nonlinear Model Predictive Control
Wang et al. Research on Integrated Control Method for Helicopter/Turboshaft Engine with Variable Rotor Speed Based on the Error Between Engine Required and Real Output Torque
Chen et al. Analysis and simulation of thrust management system for large plane

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 18922736

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 18922736

Country of ref document: EP

Kind code of ref document: A1