CN113359484B - Aero-engine model prediction control method and device based on semi-alternative optimization - Google Patents
Aero-engine model prediction control method and device based on semi-alternative optimization Download PDFInfo
- Publication number
- CN113359484B CN113359484B CN202110841852.3A CN202110841852A CN113359484B CN 113359484 B CN113359484 B CN 113359484B CN 202110841852 A CN202110841852 A CN 202110841852A CN 113359484 B CN113359484 B CN 113359484B
- Authority
- CN
- China
- Prior art keywords
- control
- ctrl
- engine
- vector
- con
- Prior art date
- Legal status (The legal status 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 status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 62
- 238000005457 optimization Methods 0.000 title claims abstract description 52
- 230000007246 mechanism Effects 0.000 claims abstract description 13
- 238000004088 simulation Methods 0.000 claims description 18
- 230000008859 change Effects 0.000 claims description 15
- 239000011159 matrix material Substances 0.000 claims description 12
- 238000005192 partition Methods 0.000 claims description 11
- 230000003044 adaptive effect Effects 0.000 claims description 7
- 230000009471 action Effects 0.000 claims description 6
- 238000000638 solvent extraction Methods 0.000 claims description 5
- 230000017105 transposition Effects 0.000 claims description 5
- 238000010276 construction Methods 0.000 claims description 3
- 230000008569 process Effects 0.000 description 10
- 230000005477 standard model Effects 0.000 description 8
- 230000000694 effects Effects 0.000 description 7
- 238000004422 calculation algorithm Methods 0.000 description 5
- 238000013461 design Methods 0.000 description 4
- 239000000446 fuel Substances 0.000 description 4
- 238000005070 sampling Methods 0.000 description 4
- MZWGYEJOZNRLQE-KXQOOQHDSA-N 1-stearoyl-2-myristoyl-sn-glycero-3-phosphocholine Chemical compound CCCCCCCCCCCCCCCCCC(=O)OC[C@H](COP([O-])(=O)OCC[N+](C)(C)C)OC(=O)CCCCCCCCCCCCC MZWGYEJOZNRLQE-KXQOOQHDSA-N 0.000 description 3
- 230000008901 benefit Effects 0.000 description 3
- 150000001875 compounds Chemical class 0.000 description 3
- 230000007547 defect Effects 0.000 description 3
- 230000004044 response Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 2
- 238000002485 combustion reaction Methods 0.000 description 2
- 238000010586 diagram Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000001105 regulatory effect Effects 0.000 description 2
- RZVHIXYEVGDQDX-UHFFFAOYSA-N 9,10-anthraquinone Chemical compound C1=CC=C2C(=O)C3=CC=CC=C3C(=O)C2=C1 RZVHIXYEVGDQDX-UHFFFAOYSA-N 0.000 description 1
- 230000009194 climbing Effects 0.000 description 1
- 230000001276 controlling effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000002372 labelling Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
Images
Classifications
-
- 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
- G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
- G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
- G05B13/04—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators
- G05B13/042—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric involving the use of models or simulators in which a parameter or coefficient is automatically adjusted to optimise the performance
Landscapes
- Engineering & Computer Science (AREA)
- Health & Medical Sciences (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Computation (AREA)
- Medical Informatics (AREA)
- Software Systems (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Automation & Control Theory (AREA)
- Feedback Control In General (AREA)
Abstract
The invention discloses a prediction control method of an aero-engine model based on semi-alternating optimization, which comprises the steps of constructing a corresponding secondary planning problem with constraint by utilizing a discrete state space model and combining the limitation quantity of an engine, constraint conditions of an execution mechanism, instruction input and an objective function, and solving the problem; applying the obtained control sequence to the engine control of the current control period; tracking the real engine working state in real time through the self-adaptive airborne model, and carrying out real-time online linearization to obtain a discrete state space model for prediction; the method comprises the steps of introducing a semi-alternating control sequence to be optimized into a model predictive control structure to replace an existing control sequence to be optimized, constructing a new parameter predictive equation form, reconstructing a new constrained quadratic programming problem to solve, and further integrally realizing multivariable control, limited protection control and optimized control on the aircraft engine, so that the scale of the optimized problem is greatly reduced, the real-time performance is improved, and the control precision is improved.
Description
Technical Field
The invention relates to an aero-engine control method, in particular to an aero-engine model prediction control method, and belongs to the technical field of aero-engine control.
Background
An aircraft engine control system is a necessary guarantee for safe and efficient operation of an aircraft engine. With the continuous improvement of the design and manufacturing technology of the aircraft engine, the variable which can be regulated by the engine is increasing, so that the control system of the engine is greatly advancing towards the direction of multi-variable control. In addition, the work of the aircraft engine needs to meet the requirements of flight tasks and corresponding safety boundary requirements, and the tasks to be processed by the control system are increasingly complicated and diversified.
The traditional aeroengine control system usually adopts a mode of combining a plurality of control loops to realize multivariable control, or relies on single-variable Min-Max logic to realize the integration of a main control task and a limiting protection control task. In addition, due to the strong non-linear characteristic of the aircraft engine, multiple sets of controllers are required to be designed for parameter scheduling, so that the control effect is poor at an operating point without controller design.
In contrast, novel model predictive control is becoming increasingly favored over highly integrated control architectures. The original complex multi-control loop design problem is converted into a constrained optimization problem, real-time optimal control on the controlled object is realized by solving the optimization problem in real time on line, and repeated design work caused by nonlinearity of the controlled object can be effectively avoided. This control architecture is very attractive for aircraft engine control. Because the requirements in terms of the aircraft engine limit protection control can be effectively abstracted into a series of constraints, and the flight control tasks that are desired to be achieved can be effectively abstracted into corresponding objective functions. In other words, the aircraft engine control problem can be effectively abstracted into an optimization problem required to be solved by model predictive control, so that the model predictive control is used for replacing a traditional engine control structure to form a feasible development direction.
However, model predictive control is also becoming a computationally intensive method because it requires on-line real-time solution of constrained optimization problems, which also greatly limits the application of this control method in aircraft engine control systems. This is because the computing power of the aircraft engine control system is very limited, and the scale of the optimization problem constructed by the model predictive control increases significantly as the control time domain and the prediction time domain increase. Particularly, when facing an aircraft engine which is an object with a plurality of regulating variables, the size of the optimization problem of the required solution can be multiplied by the requirement of the multivariable control, and the optimization problem is far beyond the computing capability of an engine control system and can not meet the requirement of real-time control. However, the solution accuracy is sacrificed in order to meet the real-time control requirement, so that the output of the obtained controller is difficult to reach the optimal solution, the control effect is deteriorated, and the advantages of the model predictive control method cannot be exerted.
Disclosure of Invention
The invention aims to overcome the defects of the prior art and provides a model predictive control method of an aero-engine based on semi-alternative optimization, which can effectively reduce the calculation scale of the optimization solution problem required by the model predictive control method, greatly improve the real-time performance, further improve the control effect, realize multi-mode and multi-variable control of the aero-engine, and simultaneously ensure that the model predictive control method meets the flight task requirement and the safe working requirement of the engine.
The invention specifically adopts the following technical scheme to solve the technical problems:
a prediction control method of an aircraft engine model based on semi-alternating optimization is characterized in that a discrete state space model is utilized and combined with a limiting quantity constraint condition, an execution mechanism constraint condition, instruction input and an objective function of an engine to construct a corresponding secondary planning problem with constraint and solve the problem; applying the obtained control sequence to the engine control of the current control period; tracking the real engine working state in real time through the self-adaptive airborne model, and carrying out real-time online linearization to obtain a discrete state space model for prediction to be used for the optimization problem construction at the next moment; the constrained quadratic programming problem is constructed according to the following method:
step A, constructing a semi-alternate control sequence to be optimized delta U specially used for the current time k:
where the subscripts k +0, k +1, …, k + r-1, k + r denote the control vector blocks for engine control at the current time k, future times k +1, k + r-1, and k + r, and Δ u1、Δu2、Δui、Δui+1、ΔurInputting a partitioning vector of a vector u for the discrete state space model, wherein a superscript T represents transposition, a represents deviation, p is the dimension of a control vector of the discrete state space model, r is the number of partitions of all elements in the control vector of the discrete state space model, and the dimension of each component partition is z1、…、zr，
Representing a set of rational numbers; the subscript j indicates the value for the future time k + nu-1 control vector component index, which takes on values of r and nuTo changeIt is determined as follows:
the component indices are repeated sequentially in the order of i, i +1, i +2, …, r,1,2, …, i-1, i, i +1, i +2, …, r,1,2, …, i-1 …, the nth from the first iuThe number is the value of j;
step B, according to the control time domain nuAnd predicting the time domain nyConstructing a prediction equation of the controlled quantity and the limited quantity of the engine with the following form:
where the subscript ctrl represents the controlled quantity, con represents the limiting quantity, ε represents the additional state increment due to linearization when the discrete-state model is obtained by linearization, Δ YctrlAnd Δ YconRespectively representing the predicted values of the controlled quantity and the limited quantity in the prediction time domain, Pctrl、Hctrl、Lctrl、Pcon、HconAnd LconCoefficient matrices representing controlled quantity and limiting quantity prediction equations, respectively, and HctrlAnd HconAt each simulation instant k changes with the change in Δ U, and
in the formula, x and y represent a state vector and an output vector of the discrete state space model, respectively, and y representsctrlAnd yconSubscripts "ctrl" and "con" of (c) are written uniformly in the lower right-hand corner of the square brackets to simplify the expression;
step C, obtaining a secondary objective function only containing delta U as the quantity to be optimized according to the objective function and by using the prediction equation; and constructing a linear constraint condition by using a definitional formula of a control sequence delta U to be optimized and the prediction equation according to the restriction of the engine restriction, the restriction of the upper limit and the lower limit of the actuating mechanism and the restriction of the maximum action value of the actuating mechanism at each step, thereby forming a quadratic programming problem with restriction.
Preferably, the objective function is selected from a plurality of different objective functions according to a control target.
Further preferably, the plurality of different objective functions includes at least a conventional control objective function and a thrust optimization control objective function;
the conventional control objective function is specifically as follows:
Jobj＝eTW1e+ΔUTW2ΔU
wherein e represents a control error vector, W1And W2Is a weighting coefficient diagonal matrix with proper dimension; wherein e is defined as follows:
e＝rctrl-Yctrl
in the formula, rctrlTo predict time domain nyInstruction vector of upper controlled quantity, YctrlTo predict time domain nyA predicted value of the upper controlled quantity;
the thrust optimization control objective function is specifically as follows:
in the formula, W2And W3Is an adaptive diagonal array of weighting coefficients, ΓctrlTo predict time domain nyAnd (4) predicting the upward thrust.
Preferably, when solving the quadratic programming problem with constraints, an external penalty function method is used for converting the constraint problem into an unconstrained problem to be solved.
Based on the same inventive concept, the following technical scheme can be obtained:
an aircraft engine model predictive control device based on semi-alternating optimization comprises:
the control sequence generation module is used for constructing a corresponding secondary planning problem with constraint by using a discrete state space model and combining a limiting quantity constraint condition, instruction input and an objective function of the engine, solving the problem, and applying the obtained control sequence to the engine control of the current control period;
the self-adaptive airborne model is used for tracking the real engine working state in real time and carrying out real-time online linearization to obtain a discrete state space model for prediction;
the constrained quadratic programming problem is constructed according to the following method:
step A, constructing a semi-alternate control sequence to be optimized delta U specially used for the current time k:
where the subscripts k +0, k +1, …, k + r-1, k + r denote the control vector blocks for engine control at the current time k, future times k +1, k + r-1, and k + r, and Δ u1、Δu2、Δui、Δui+1、ΔurInputting a partitioning vector of a vector u for the discrete state space model, wherein a superscript T represents transposition, a represents deviation, p is the dimension of a control vector of the discrete state space model, r is the number of partitions of all elements in the control vector of the discrete state space model, and the dimension of each component partition is z1、…、zr，
Representing a set of rational numbers; the subscript j indicates the value for the future time k + nu-1 control vector component index, which takes on values of r and nuBut instead, is determined as follows:
the component indices are repeated sequentially in the order of i, i +1, i +2, …, r,1,2, …, i-1, i, i +1, i +2, …, r,1,2, …, i-1 …, the nth from the first iuThe number is the value of j;
step B, according to the control time domain nuAnd predicting the time domain nyConstructing a prediction equation of the controlled quantity and the limited quantity of the engine with the following form:
where the subscript ctrl represents the controlled quantity, con represents the limiting quantity, ε represents the additional state increment due to linearization when the discrete-state model is obtained by linearization, Δ YctrlAnd Δ YconRespectively representing the predicted values of the controlled quantity and the limited quantity in the prediction time domain, Pctrl、Hctrl、Lctrl、Pcon、HconAnd LconCoefficient matrices representing controlled quantity and limiting quantity prediction equations, respectively, and HctrlAnd HconAt each simulation instant k changes with the change in Δ U, and
in the formula, x and y represent a state vector and an output vector of the discrete state space model, respectively, and y representsctrlAnd yconThe subscripts "ctrl" and "con" of (a) are written uniformly in the lower right corner of the square bracket to simplify the expression;
step C, obtaining a secondary objective function only containing delta U as the quantity to be optimized according to the objective function and by using the prediction equation; and constructing a linear constraint condition by using a definitional formula of a control sequence delta U to be optimized and the prediction equation according to the restriction of the engine restriction, the restriction of the upper limit and the lower limit of the actuating mechanism and the restriction of the maximum action value of the actuating mechanism at each step, thereby forming a quadratic programming problem with restriction.
Preferably, the objective function is selected from a plurality of different objective functions according to a control target.
Further preferably, the plurality of different objective functions includes at least a conventional control objective function and a thrust optimization control objective function;
the conventional control objective function is specifically as follows:
Jobj＝eTW1e+ΔUTW2ΔU
wherein e represents a control error vector, W1And W2Is a weighting coefficient diagonal matrix with proper dimension; it is composed ofWherein e is defined as follows:
e＝rctrl-Yctrl
in the formula, rctrlTo predict time domain nyInstruction vector of upper controlled variable, YctrlTo predict time domain nyA predicted value of the upper controlled quantity;
the thrust optimization control objective function is specifically as follows:
Jobj＝ΔUTW2ΔU-ΓcTtrlW3Γctrl
in the formula, W2And W3Is an adaptive diagonal array of weighting coefficients, ΓctrlTo predict time domain nyAnd (4) predicting the upward thrust.
Preferably, when solving the quadratic programming problem with constraints, an external penalty function method is used for converting the constraint problem into an unconstrained problem to be solved.
Compared with the prior art, the technical scheme of the invention has the following beneficial effects:
(1) capability of multi-mode control: the invention fully utilizes the advantage of the optimization solution of the control method, sets the optimization objective functions under different control targets, and can realize different control tasks according to different optimization objective functions;
(2) high real-time performance: by introducing the semi-alternating control sequence to be optimized, the invention effectively reduces the scale of the quadratic programming problem to be solved at each sampling moment, obviously reduces the demand of on-line calculated amount and can effectively improve the real-time performance of control;
(3) high control precision: according to the invention, by introducing the semi-alternating control sequence to be optimized, the scale of the quadratic programming problem to be solved at each sampling moment is effectively reduced, so that the solving difficulty of the whole optimization problem is reduced, the optimal solution can be more effectively found, and the high-precision engine control is realized;
(4) the method has the following self-adaptive control capacity: according to the method, the available margin of the actually available limiting protection variable can be estimated by means of an airborne model according to the current working point of the engine, the thrust of the engine is increased in a self-adaptive manner, and the working potential of the engine is fully utilized;
(5) universality and portability: the method is implemented based on an airborne model and an optimization algorithm, and is suitable for various aero-engines capable of generating thrust.
Drawings
FIG. 1 is a block diagram of a typical two-shaft mixed-row turbofan engine;
FIG. 2 is a schematic diagram of the overall structure of a model predictive control apparatus for an aircraft engine according to the present invention;
FIG. 3 illustrates flight conditions used in simulation experiments in an embodiment of the present invention;
FIG. 4 is a graph comparing the thrust under the conventional control objective and the standard model predictive control according to the present invention;
FIG. 5 is a comparison of total inlet temperature of a high pressure turbine under conventional control objectives and under standard model predictive control in accordance with the present invention;
FIG. 6 is a graph comparing fan surge margin under the conventional control objective of the present invention with that under the standard model predictive control;
FIG. 7 is a graph comparing the change in main fuel flow under the conventional control objective of the present invention with that under the standard model predictive control;
FIG. 8 is a curve of the change in throat area of the nozzle under the conventional control objective and the standard model predictive control in accordance with the present invention;
FIG. 9 is a force comparison graph under the conventional control objective of the whole process of the present invention and under the multi-mode control of the present invention;
FIG. 10 is a comparison of total inlet temperature of a high pressure turbine under conventional control objectives throughout the present invention versus multi-mode control in accordance with the present invention;
FIG. 11 is a graph comparing fan surge margin under the global conventional control objective of the present invention with the multi-mode control of the present invention;
FIG. 12 is a graph comparing the main fuel flow variation under the conventional control objective of the present invention and under the multi-mode control of the present invention during the whole course;
FIG. 13 is a curve of the change of throat area of the lower jet nozzle under the conventional control target in the whole process and the multi-mode control of the invention.
Detailed Description
Aiming at the defects of the existing prediction control method of the model of the aero-engine, the invention is based on the idea of semi-alternating optimization, a control sequence to be optimized in a semi-alternating form is defined and introduced into the model prediction control structure to replace the existing control sequence to be optimized, a new parameter prediction equation form is constructed, a new quadratic programming problem with constraint is reconstructed to solve, and then the multivariable control, the limited protection control and the optimization control of the aero-engine are integrally realized, so that the scale of the optimization problem is greatly reduced, the real-time performance is improved, and the control precision is improved. The method is applicable to various power machines which can establish a complex analytical model, have a plurality of adjustable variables and can generate thrust, and the method comprises but is not limited to a turbojet engine, a turbofan engine, a turboprop engine, a variable cycle engine, a turbine-based ramjet combined engine and the like.
The invention provides a prediction control method of an aero-engine model based on semi-alternating optimization, which utilizes a discrete state space model and combines a limiting quantity constraint condition, an execution mechanism constraint condition, instruction input and an objective function of an engine to construct a corresponding quadratic programming problem with constraint and solve the problem; applying the obtained control sequence to the engine control of the current control period; tracking the real engine working state in real time through the self-adaptive airborne model, and carrying out real-time online linearization to obtain a discrete state space model for prediction to be used for the optimization problem construction at the next moment;
wherein the constrained quadratic programming problem is constructed according to the following method:
step A, constructing a semi-alternate control sequence to be optimized delta U specially used for the current time k:
where the subscripts k +0, k +1, …, k + r-1, k + r denote the control vector blocks for engine control at the current time k, future times k +1, k + r-1, and k + r, and Δ u1、Δu2、Δui、Δui+1、ΔurInputting a partitioning vector of a vector u for the discrete state space model, wherein a superscript T represents transposition, a represents deviation, p is the dimension of a control vector of the discrete state space model, r is the number of partitions of all elements in the control vector of the discrete state space model, and the dimension of each component partition is z1、…、zr，
Representing a set of rational numbers; the subscript j indicates the value for the future time k + nu-1 control vector component index, which takes on values of r and nuBut instead, is determined as follows:
the component indices are repeated sequentially in the order of i, i +1, i +2, …, r,1,2, …, i-1, i, i +1, i +2, …, r,1,2, …, i-1 …, the nth from the first iuThe number is the value of j;
step B, according to the control time domain nuAnd predicting the time domain nyConstructing a prediction equation of the controlled quantity and the limited quantity of the engine with the following form:
where the subscript ctrl represents the controlled quantity, con represents the limiting quantity, ε represents the additional state increment due to linearization when the discrete-state model is obtained by linearization, Δ YctrlAnd Δ YconRespectively representing the predicted values of the controlled quantity and the limited quantity in the prediction time domain, Pctrl、Hctrl、Lctrl、Pcon、HconAnd LconCoefficient matrices representing controlled quantity and limiting quantity prediction equations, respectively, and HctrlAnd HconAt each simulation instant k changes with the change in Δ U, and
in the formula, x,y represents the state vector and the output vector of the discrete state space model, respectively, yctrlAnd yconThe subscripts "ctrl" and "con" of (a) are written uniformly in the lower right corner of the square bracket to simplify the expression;
step C, obtaining a secondary objective function only containing delta U as the quantity to be optimized according to the objective function and by using the prediction equation; and constructing a linear constraint condition by using a definitional formula of a control sequence delta U to be optimized and the prediction equation according to the restriction of the engine restriction, the restriction of the upper limit and the lower limit of the actuating mechanism and the restriction of the maximum action value of the actuating mechanism at each step, thereby forming a quadratic programming problem with restriction.
Preferably, the objective function is selected from a plurality of different objective functions according to a control target.
Further preferably, the plurality of different objective functions includes at least a conventional control objective function and a thrust optimization control objective function;
the conventional control objective function is specifically as follows:
Jobj＝eTW1e+ΔUTW2ΔU
wherein e represents a control error vector, W1And W2Is a weighting coefficient diagonal matrix with proper dimension; wherein e is defined as follows:
e＝rctrl-Yctrl
in the formula, rctrlTo predict time domain nyInstruction vector of upper controlled variable, YctrlTo predict time domain nyA predicted value of the upper controlled quantity;
the thrust optimization control objective function is specifically as follows:
in the formula, W2And W3Is an adaptive diagonal array of weighting coefficients, ΓctrlTo predict time domain nyAnd (4) predicting the upward thrust.
Preferably, when solving the constrained quadratic programming problem, the constrained problem is converted into an unconstrained problem for solving by using an external penalty function method, so as to ensure that a reasonable approximate solution can be found even if the constructed constrained quadratic programming problem has no feasible solution.
After a quadratic programming problem with constraints is formed by the method, solving the quadratic programming problem to obtain an optimal control sequence delta U at the current moment; the value of the control vector index i in the control sequence is then updated: if i +1> r, i equals 1, otherwise i equals i + 1.
For the public understanding, the following takes the process of applying the model predictive control method of the semi-alternating optimization-based aircraft engine described in the invention to the double-shaft mixed-row turbofan engine shown in fig. 1 as an example to describe the technical scheme of the invention in detail.
Section 1 in fig. 1 is an inlet of an air inlet; 2, the cross section is an air inlet outlet and a fan inlet; the section of the 22 is a fan outlet; 25 and 15 sections are inlets and outlets of inner and outer culverts; the section 3 is a compressor outlet and a combustion chamber inlet; 4 the section is a combustion chamber outlet; the section of the high-pressure turbine is 41; the section 42 is a high-pressure turbine outlet; the 45 section is a low-pressure turbine inlet; the section of 46 is a low-pressure turbine outlet; 16 and 6 sections are respectively an outer culvert outlet and an inner culvert outlet; 8, the section is a tail nozzle throat; the section 9 is a tail nozzle outlet.
FIG. 2 shows the overall structure of the model predictive control device for an aircraft engine based on semi-alternating optimization, which tracks the real engine working state in real time through an adaptive airborne model and obtains a discrete state space model for prediction through real-time online linearization; constructing a corresponding secondary planning problem with constraint by using a discrete state space model and combining a limiting quantity constraint condition, instruction input and an objective function of the engine, and solving the problem; and applying the obtained control sequence to the engine control of the current control period. Y in FIG. 2obThe representation sensor measures engine parameters, the delta u is a control variable which is obtained by optimizing and solving the controller and is used for controlling at the current moment, and the u is an engine control variable output by the controller.
In this embodiment, the illustrated engine has two adjustable variables, namely the main fuel flow mfbAnd the throat area A of the tail nozzle8. The engine control system realizes the multi-variable control of the corresponding engine by adjusting the two variables. Thus, the discrete state space model control vector dimension p of the engine is 2, and in this embodiment, the vector is divided into two groups, i.e., r is 2, which is denoted as Δ u1＝Δmfb，Δu2＝ΔA8So that the dimension of each group is z1＝z21. Therefore, in the present embodiment, each partition vector is reduced to a scalar form, which is beneficial for public understanding. However, it should be noted that the implementation process using the block vector form is the same as that of the present embodiment, and only the dimension of the matrix specifically involved may be different, so that the present embodiment is still representative.
According to the actual control requirement, the controlled quantity is the thrust force F in the embodiment, and the limited quantity is the low-pressure rotor speed nfHigh-pressure rotor speed ncHigh pressure turbine inlet temperature T41Surge margin S of fanmlSurge margin S of compressormhThe relevant limits are shown in the following table:
TABLE 1 limiting parameters used for model predictive control simulation
In the invention, the real engine working state is tracked in real time through the self-adaptive airborne model, and the discrete state space model for prediction is obtained through real-time online linearization. Various self-adaptive airborne models which can effectively reflect the current working state of the engine in real time, have self-updating capability and can generate a discrete state space model can be used in the invention. Without loss of generality, in the embodiment, steps A, B and C in a chinese patent "online prediction method of air path parameters of an aircraft engine" (patent No. ZL201910788314.5) are adopted to realize real-time updating of a discrete state space model for prediction required in the present invention.
For convenience of description of the embodiment, it is assumed without loss of generality that the current time is time k, this time i is 1, and the discrete state space model for prediction has been obtained at the last time:
where x, y, u represent the engine state vector, output vector and input vector, respectively, A and B are the system matrices, C and D are the output matrices, subscript m represents the time of the discrete state space model itself, xk-1,uk-1,yk-1The operating point of the engine at the time k-1 is shown, i.e. the discrete model is linearized at this operating point.
Note that, if the engine control is started for the first time, i is initialized to 1, and the discrete state space model is initialized.
A control sequence Δ U of a semi-alternating form specific to the current instant k is defined as follows:
if the time domain n is controlleduIs odd, then
If the time domain n is controlleduIs an even number, then
Wherein the subscripts k +0, k +1, k +2, …, k + nu-1 denotes for the current time k, future times k +1, k +2 and k + nu-1 control quantity block of engine control.
For the purposes of understanding, the control sequence Δ U at the next time k +1 is given in particular as a comparison, noting that this time i will be equal to 2, so i +1 will be greater than r, so the value of i +1 needs to be taken to be 1:
if the time domain n is controlleduIs odd, then
If the control time domain nuIs an even number, then
Wherein the subscripts (k +1) +0, (k +1) +1, (k +1) +2, …, (k +1) + nu-1 denotes for the current time k +1 (for the next time k +1), the future time (k +1) +1, (k +1) +2 and (k +1) + nu-1 control quantity block of engine control.
According to the discrete state space model obtained by the k-1 time linearization, rewriting the following discrete state space model:
wherein the subscript ctrl represents the controlled amount, con represents the limiting amount, and
and in this embodiment, yctrl＝F，ycon＝[nf,nc,T41,Sml,Smh]T。
Using equation (4), a corresponding prediction equation is constructed:
in the formula,. DELTA.YctrlAnd Δ YconRespectively representing the predicted values of the controlled quantity and the limited quantity in the prediction time domain, Pctrl、Hctrl、 Lctrl、Pcon、HconAnd LconCoefficient matrices representing controlled quantity and limiting quantity prediction equations, respectively, and HctrlAnd HconAt each simulation timeThe moment k varies with the change of Δ U, and
note that, in the equations (8) and (9), for the sake of convenience of labeling, subscripts ctrl of C and D matrices are written in the lower right corner of the square bracket, and the corresponding C and D matrices each represent a matrix related to a controlled quantity.
Limited space due to Pcon、HconAnd LconAnd Pctrl、HctrlAnd LctrlHave exactly the same formal structure, and differ only in the matrix values and dimensions, and in the subscript "con", so nothing is omitted here.
Here too, H is given at the next moment k +1ctrlTo facilitate a public-based understanding:
the desired objective function is selected based on the engine control objective. It should be noted that different control target selections may affect the control effect of the model predictive control method of the present invention on the controlled object, but may not affect the implementation process of the technical solution of the present invention. Therefore, the engine control target is selected as the normal control target without loss of generality, and the following objective function J is adoptedobjDefining:
Jobj＝eTW1e+ΔUTW2ΔU (13)
wherein e represents a control error vector, W1And W2Is a diagonal matrix of adaptive weighting coefficients. Wherein e is defined as follows:
e＝rctrl-Yctrl (14)
in the formula, rctrlTo predict time domain nyAnd an instruction vector of the upper controlled quantity.
In the formula, rk+0、rk+1、rk+nyThe thrust control command at each time in the time domain is predicted.
It should be noted that other optimization control targets such as the thrust optimization control target may be selected as needed, and the implementation process of the subsequent step is not affected. The upper and lower limits of the engine limit are respectively marked as yubAnd ylb:：
In the formula, the subscript ub represents an upper limit and lb represents a lower limit.
Then the derived achievable limit constraint is:
in the formula (I), the compound is shown in the specification,
considering the actuator constraints, there are the upper and lower limits of the actuator and the maximum value of each step of motion:
in the formula (I), the compound is shown in the specification,
in the formula u2,ub,u2,lb,u1,ub,u1,lb,u2,a,u1,aAre each u2And u1Corresponding upper and lower limits of the actuating mechanism and the maximum value of each step of action.
Note that Uub、UlbAnd UaDynamically changing according to different time deltau.
From this, the following quadratic programming problem with constraints is constructed:
in the formula (I), the compound is shown in the specification,
in this embodiment, an external penalty function method is used to convert the constraint problem into an unconstrained problem. Thus, a new equivalent objective function J is redefined by using the penalty function methodnew：
In which σ is a penalty factor and has
To solve equation (23), a number of existing, well-established iterative optimization algorithms may be selected. In this embodiment, the DFP algorithm is used to solve equation (23) to obtain an optimal control sequence Δ U.
Updating the value of the control vector index i in the control sequence: since i +1 is 2, r, i is i + 1.
After the steps, according to the solved optimal control sequence delta U, outputting the following control signals to control the engine at the current time:
then, updating the self-adaptive model and obtaining a new discrete state space model; and judging whether the control is finished or not, if so, finishing the whole process, and if not, entering the next control period if k is k + 1.
In order to verify the effectiveness of the invention, a simulation experiment within the flight envelope is developed. The simulated sampling period is 20 milliseconds and the simulated flight conditions are shown in figure 3.
FIGS. 4 to 8 show thrust F and total high-pressure turbine inlet temperature T under control of the method of the invention when the control target is selected as the conventional control target41And fan surge margin SmlThe remaining engine variables are no longer given due to the untriggered limit. And in order to facilitate comparison of results, simulation results of the standard model predictive control method are also provided, so that control effects of different control algorithms can be compared. In the figure, "Command" indicates a thrust Command, "Limit" indicates a Limit line, "SMPC" indicates standard model predictive control, and "SAOMPC" indicates the present invention. SMPC and SAOMPC adopted by simulation both adopt the same parameter setting, including predicting time domain nyControl time domain nuWeight coefficient W1And W2. Predicting the time domain nyAnd control time domain nuSet to 6 and 3, respectively, that is, predict the change of the output parameter at the future simulation time k +0 to k +5, optimize the change of the control quantity at the future simulation time k +0 to k +2And (4) transforming. In addition, the weight matrix W1Arranged as a unit array, W2Arranged as a diagonal matrix and in which m is related to the fuel quantityfbThe relevant element is set as 1 and the area A of the throat of the tail nozzle8The associated element is set to 5. It can be seen from fig. 4 to 8 that the present invention can effectively control the engine.
To further evaluate the control effect over a large envelope, the root mean square error was calculated as follows:
wherein N is the number of samples, refIs the command value and y is the controlled quantity. For a simulation duration of 60 seconds, the sampling period is 20 milliseconds, so the number of samples N is 3000.
Table 2 shows the root mean square error for different combinations of control and prediction time domains. Note in particular that the root mean square error in the table is x 10-2To display more decimal place numbers.
TABLE 2 comparison of control errors for two model predictive control algorithms
As can be seen from table 2, under different combinations of the control time domain and the prediction time domain, the method of the present invention can achieve a smaller control error than the existing standard model prediction control method, which indicates that the method of the present invention has advantages in control accuracy compared with the existing method.
The calculation time consumption test is carried out on the STM32F 767-based microcontroller unit, and the maximum time consumption t under different control time domain and prediction time domain combinations is shown in Table 3mavg. The maximum time consumption refers to the average time consumption of 10% of samples with the maximum time consumption in the time consumption for solving the optimization problem in the whole simulation process.
TABLE 3 comparison of time consumption for optimization solutions under different combinations of control time domain and prediction time domain
As can be seen from tables 2 and 3, compared with the conventional SMPC control method, the method of the present invention has a better control effect and can greatly improve the real-time performance.
To verify the multimode control capability of the present invention, the following simulations were performed: the change in the flight condition was the first 20s of the change in the flight condition shown in fig. 3, and thereafter the flight condition was kept constant. The thrust optimization control target is selected at the beginning of the climbing phase, and the normal control target is switched back after the flat flight stability is achieved (16 th s of simulation time). Control time domain is taken to be 3, prediction time domain is taken to be 6, and weight W is divided3The parameter settings other than the above are also selected as the parameter settings in the conventional control target. Weight W3The diagonal matrix is taken, the values of the nonzero elements are the same and are all taken as 0.02. I.e. the simulation has selected different control targets at different time periods.
FIGS. 9 to 13 show thrust F and total inlet temperature T of the high-pressure turbine in the case of multi-mode control under the control of the method of the present invention41And fan surge margin SmlThe response curve and the control variable variation curve of (a), the remaining engine variables are not given because of the triggerless limit. In the figure, "Command" indicates a thrust Command change, "Normal throttle control" indicates a parameter response in the full-stroke Normal control target mode, and "w ═ 0.02" indicates a parameter response in the multi-mode control. As can be seen from the figure, the invention can effectively improve the thrust of the engine, and can effectively change the control target according to the requirement in the flight process, thereby realizing the multi-mode control.
In conclusion, the method can effectively make up for the defects of the prior control technology, effectively improve the control precision, greatly reduce the time consumption of the optimization solving process, improve the real-time performance, effectively realize the multi-mode control and effectively lay a foundation for engineering application.
Claims (8)
1. A prediction control method of an aircraft engine model based on semi-alternating optimization is characterized in that a discrete state space model is utilized and combined with a limiting quantity constraint condition, an execution mechanism constraint condition, instruction input and an objective function of an engine to construct a corresponding secondary planning problem with constraint and solve the problem; applying the obtained control sequence to the engine control of the current control period; tracking the real engine working state in real time through the self-adaptive airborne model, and carrying out real-time online linearization to obtain a discrete state space model for prediction to be used for the optimization problem construction at the next moment; the method is characterized in that the secondary planning problem with constraint is constructed according to the following method:
step A, constructing a semi-alternate control sequence to be optimized delta U specially used for the current time k:
where the subscripts k +0, k +1, …, k + r-1, k + r denote the control vector blocks for engine control at the current time k, future times k +1, …, k + r-1, and k + r, Δ u1、Δu2、Δui、Δui+1、ΔurInputting a partitioning vector of a vector u for the discrete state space model, wherein a superscript T represents transposition, a represents deviation, p is the dimension of a control vector of the discrete state space model, r is the number of partitions of all elements in the control vector of the discrete state space model, and the dimension of each component partition is z1、…、zr，
Representing a set of rational numbers; the subscript j indicates the value for the future time k + nu-1 control vector component index, which takes on values of r and nuBut instead, is determined as follows:
the component indices are repeated sequentially in the order of i, i +1, i +2, …, r,1,2, …, i-1, i, i +1, i +2, …, r,1,2, …, i-1 …, the nth from the first iuThe number is the value of j;
step B, according to the control time domain nuAnd predicting the time domain nyConstructing a prediction equation of the controlled quantity and the limited quantity of the engine with the following form:
where the subscript ctrl represents the controlled quantity, con represents the limiting quantity, ε represents the additional state increment due to linearization when the discrete-state model is obtained by linearization, Δ YctrlAnd Δ YconRespectively representing the predicted values of the controlled quantity and the limited quantity in the prediction time domain, Pctrl、Hctrl、Lctrl、Pcon、HconAnd LconCoefficient matrices representing controlled quantity and limiting quantity prediction equations, respectively, and HctrlAnd HconAt each simulation instant k changes with the change in Δ U, and
in the formula, x and y represent a state vector and an output vector of the discrete state space model, respectively, and y representsctrlAnd yconThe subscripts "ctrl" and "con" of (a) are written uniformly in the lower right corner of the square bracket to simplify the expression;
step C, obtaining a secondary objective function only containing delta U as the quantity to be optimized according to the objective function and by using the prediction equation; and constructing a linear constraint condition by utilizing a definition formula of a control sequence delta U to be optimized and the prediction equation according to engine limit quantity constraint, actuator upper and lower limit constraint and actuator maximum value constraint of each step of action, thereby forming a quadratic programming problem with constraint.
2. The method of claim 1, wherein the objective function is selected from a plurality of different objective functions based on a control objective.
3. The semi-alternating optimization-based aircraft engine model predictive control method of claim 2, wherein the plurality of different objective functions includes at least a conventional control objective function and a thrust optimization control objective function;
the conventional control objective function is specifically as follows:
Jobj＝eTW1e+ΔUTW2ΔU
wherein e represents a control error vector, W1And W2Is a weighting coefficient diagonal matrix with proper dimension; wherein e is defined as follows:
e＝rctrl-Yctrl
in the formula, rctrlTo predict time domain nyInstruction vector of upper controlled variable, YctrlTo predict time domain nyA predicted value of the upper controlled quantity;
the thrust optimization control objective function is specifically as follows:
in the formula, W2And W3Is an adaptive diagonal array of weighting coefficients, ΓctrlTo predict time domain nyAnd (4) predicting the upward thrust.
4. The method for predictive control of an aircraft engine model based on semi-alternating optimization according to claim 1, wherein in solving the constrained quadratic programming problem, an external penalty function method is used to convert the constrained problem into an unconstrained problem for solving.
5. An aircraft engine model predictive control apparatus based on semi-alternating optimization, comprising:
the control sequence generation module is used for constructing a corresponding secondary planning problem with constraint by using a discrete state space model and combining a limiting quantity constraint condition, instruction input and an objective function of the engine, solving the problem, and applying the obtained control sequence to the engine control of the current control period;
the self-adaptive airborne model is used for tracking the real engine working state in real time and carrying out real-time online linearization to obtain a discrete state space model for prediction;
the constrained quadratic programming problem is constructed according to the following method:
step A, constructing a semi-alternate control sequence to be optimized delta U specially used for the current time k:
where the subscripts k +0, k +1, …, k + r-1, k + r denote the control vector blocks for engine control at the current time k, future times k +1, k + r-1, and k + r, and Δ u1、Δu2、Δui、Δui+1、ΔurInputting a partitioning vector of a vector u for the discrete state space model, wherein a superscript T represents transposition, a represents deviation, p is the dimension of a control vector of the discrete state space model, r is the number of partitions of all elements in the control vector of the discrete state space model, and the dimension of each component partition is z1、…、zr，
Representing a set of rational numbers; the subscript j indicates the value for the future time k + nu-1 control vector component index, which takes on values of r and nuBut instead, is determined as follows:
the component indices are repeated sequentially in the order of i, i +1, i +2, …, r,1,2, …, i-1, i, i +1, i +2, …, r,1,2, …, i-1 …, the nth from the first iuThe number is the value of j;
step B, according to the control time domain nuAnd predicting the time domain nyConstructing a prediction equation of the controlled quantity and the limited quantity of the engine with the following form:
where the subscript ctrl represents the controlled quantity, con represents the limiting quantity, ε represents the additional state increment due to linearization when the discrete-state model is obtained by linearization, Δ YctrlAnd Δ YconRespectively representing the predicted values of the controlled quantity and the limited quantity in the prediction time domain, Pctrl、Hctrl、Lctrl、Pcon、HconAnd LconCoefficient matrices representing controlled quantity and limiting quantity prediction equations, respectively, and HctrlAnd HconAt each simulation instant k changes with the change in Δ U, and
in the formula, x and y represent a state vector and an output vector of the discrete state space model, respectively, and y representsctrlAnd yconThe subscripts "ctrl" and "con" of (a) are written uniformly in the lower right corner of the square bracket to simplify the expression;
step C, obtaining a secondary objective function only containing delta U as the quantity to be optimized according to the objective function and by using the prediction equation; and constructing a linear constraint condition by using a definitional formula of a control sequence delta U to be optimized and the prediction equation according to the restriction of the engine restriction, the restriction of the upper limit and the lower limit of the actuating mechanism and the restriction of the maximum action value of the actuating mechanism at each step, thereby forming a quadratic programming problem with restriction.
6. An aircraft engine model predictive control apparatus based on semi-alternating optimization according to claim 5, wherein the objective function is selected from a plurality of different objective functions based on a control objective.
7. The semi-alternating optimization-based aircraft engine model predictive control of claim 6, wherein said plurality of different objective functions includes at least a conventional control objective function and a thrust optimization control objective function;
the conventional control objective function is specifically as follows:
Jobj＝eTW1e+ΔUTW2ΔU
wherein e represents a control error vector, W1And W2Is a weighted coefficient diagonal matrix with proper dimension; wherein e is defined as follows:
e＝rctrl-Yctrl
in the formula, rctrlTo predict time domain nyInstruction vector of upper controlled variable, YctrlTo predict time domain nyA predicted value of the upper controlled quantity;
the thrust optimization control objective function is specifically as follows:
in the formula, W2And W3Is an adaptive diagonal array of weighting coefficients, ΓctrlTo predict time domain nyAnd (4) predicting the upward thrust.
8. The aircraft engine model predictive control apparatus based on semi-alternating optimization of claim 5, wherein in solving the constrained quadratic programming problem, an external penalty function method is used to convert the constrained problem into an unconstrained problem for solution.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110841852.3A CN113359484B (en) | 2021-07-26 | 2021-07-26 | Aero-engine model prediction control method and device based on semi-alternative optimization |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110841852.3A CN113359484B (en) | 2021-07-26 | 2021-07-26 | Aero-engine model prediction control method and device based on semi-alternative optimization |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113359484A CN113359484A (en) | 2021-09-07 |
| CN113359484B true CN113359484B (en) | 2022-06-10 |
Family
ID=77540400
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110841852.3A Active CN113359484B (en) | 2021-07-26 | 2021-07-26 | Aero-engine model prediction control method and device based on semi-alternative optimization |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113359484B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113759727B (en) * | 2021-09-30 | 2023-08-29 | 中国航发控制系统研究所 | Comprehensive optimization design method for multi-variable controller of aero-engine |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102855349A (en) * | 2012-08-06 | 2013-01-02 | 南京航空航天大学 | Quick prototype design method and platform for gas path fault diagnosis for aeroengine |
| CN105114189A (en) * | 2015-06-09 | 2015-12-02 | 吉林大学 | Model prediction control system for electronic throttle valve based on FPGA |
| CN106886151A (en) * | 2017-04-17 | 2017-06-23 | 大连理工大学 | The design and dispatching method of constrained forecast controller under a kind of aero-engine multi-state |
| CN107193212A (en) * | 2017-06-26 | 2017-09-22 | 南京航空航天大学 | Aero-engine non-linear predication control method based on new grey wolf optimized algorithm |
| CN112483261A (en) * | 2020-11-16 | 2021-03-12 | 大连理工大学 | Method for resisting stress application disturbance of aircraft engine |
-
2021
- 2021-07-26 CN CN202110841852.3A patent/CN113359484B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102855349A (en) * | 2012-08-06 | 2013-01-02 | 南京航空航天大学 | Quick prototype design method and platform for gas path fault diagnosis for aeroengine |
| CN105114189A (en) * | 2015-06-09 | 2015-12-02 | 吉林大学 | Model prediction control system for electronic throttle valve based on FPGA |
| CN106886151A (en) * | 2017-04-17 | 2017-06-23 | 大连理工大学 | The design and dispatching method of constrained forecast controller under a kind of aero-engine multi-state |
| CN107193212A (en) * | 2017-06-26 | 2017-09-22 | 南京航空航天大学 | Aero-engine non-linear predication control method based on new grey wolf optimized algorithm |
| CN112483261A (en) * | 2020-11-16 | 2021-03-12 | 大连理工大学 | Method for resisting stress application disturbance of aircraft engine |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113359484A (en) | 2021-09-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108762089B (en) | Model prediction-based aero-engine online optimization and multivariable control design method | |
| Frederick et al. | Turbofan engine control design using robust multivariable control technologies | |
| Zheng et al. | Equilibrium running principle analysis on an adaptive cycle engine | |
| CN108828947B (en) | Modeling method for time-lag-containing uncertain fuzzy dynamic model of aircraft engine | |
| CN109162813B (en) | One kind being based on the modified Aeroengine Smart method for controlling number of revolution of iterative learning | |
| CN110502840B (en) | Online prediction method for gas circuit parameters of aero-engine | |
| US20140005909A1 (en) | Real time linearization of a component-level gas turbine engine model for model-based control | |
| Pang et al. | Improved nonlinear MPC for aircraft gas turbine engine based on semi-alternative optimization strategy | |
| CN113359484B (en) | Aero-engine model prediction control method and device based on semi-alternative optimization | |
| CN111679576B (en) | Variable cycle engine controller design method based on improved deterministic strategy gradient algorithm | |
| Bing et al. | Wide-range model predictive control for aero-engine transient state | |
| CN105676640A (en) | Turbofan engine acceleration process control law design method based on Bezier curve | |
| Zhou et al. | HNN-based generalized predictive control for turbofan engine direct performance optimization | |
| CN111475889A (en) | Physical and data-driven aeroengine digital twin modeling method | |
| Hao et al. | Acceleration technique for global optimization of a variable cycle engine | |
| Cai et al. | Full-envelope acceleration control method of turbofan engine based on variable geometry compound adjustment | |
| CN114237029A (en) | Turbofan engine acceleration process control method and device based on active switching logic | |
| CN110985216A (en) | Intelligent multivariable control method for aero-engine with online correction | |
| CN111734533B (en) | Turbofan engine-based model prediction method and system | |
| CN114637210A (en) | Multi-loop comprehensive control method and system for aircraft engine | |
| Juan et al. | Maneuver control at high angle of attack based on real-time optimization of integrated aero-propulsion | |
| CN113759727A (en) | Comprehensive optimization design method for multiple variable controllers of aircraft engine | |
| Cai et al. | Predictive control method for mode transition process of multi-mode turbine engine based on onboard adaptive composite model | |
| CN112523874B (en) | Multivariable limit protection control method for aircraft engine | |
| Nie et al. | Control law design of variable cycle engine based on DQN |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |