CN104881394B - Single-mode system harmonic balance subtraction unit - Google Patents
Single-mode system harmonic balance subtraction unit Download PDFInfo
- Publication number
- CN104881394B CN104881394B CN201510295236.7A CN201510295236A CN104881394B CN 104881394 B CN104881394 B CN 104881394B CN 201510295236 A CN201510295236 A CN 201510295236A CN 104881394 B CN104881394 B CN 104881394B
- Authority
- CN
- China
- Prior art keywords
- harmonic
- amplitude
- phase
- average
- subharmonic
- 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
Landscapes
- Complex Calculations (AREA)
Abstract
The invention discloses a kind of single-mode system harmonic balance subtraction unit, belong to computer engineering application field.The device includes nonlinear system arrangement components, symbolic variable definition component, calculating unit, output block.An objective function unit, the object function needed during for dynamic construction Optimization Solution average, the amplitude of each harmonic and phase are provided with calculating unit.The harmonic wave equilibrium method that the device is used, only considers the limited subharmonic in response, is classified as low-order harmonic part and high-frequency harmonic part.Optimization, iteration improve solving precision since low frequency part only includes average with first harmonic, it is continuously increased the number of times for the highest subharmonic that low frequency part is included, reduce the minimum subharmonic of high-frequency harmonic part simultaneously, obtained result is calculated as initial value using low-order harmonic part, average, the amplitude and phase of limited subharmonic in response is finally solved.The device that the present invention is provided has convergence rate and very high precision quickly.
Description
Technical field
The invention belongs to computer engineering application field, in particular to a kind of single-mode system harmonic wave equilibrium method
Device.
Background technology
Harmonic wave equilibrium method is a kind of very effective method for handling nonlinear system problem.Many actual control systems
All meet should be in this way condition, therefore obtain a wide range of applications, such as gear drive, shake in field of engineering technology
The analysis and design of the machines such as dynamic conveyer, vibration cooling machine, vibrating mill, engine, motor turning, electronic circuit and system
In.Wherein single-degree-of-freedom nonlinear system is most commonly seen, is also Nonlinear System with Multi degree of Freedom research and the basis of application, obtains
Domestic and international experts and scholars, the extensive concern of engineers and technicians, and put into substantial amounts of research.Basis and use that it is applied
Processing method be：Output/response is expressed as to the form of Fourier space, limited subharmonic is only considered, passes through mathematical method
The amplitude and phase of harmonic wave are solved, so as to obtain the frequency domain characteristic of system.In nonlinear system equation substitute into input/excitation with
The Fourier space of output/response, compares both sides coefficient to set up equation group, and then solve the amplitude and phase of harmonic wave；But
When nonlinear system becomes complexity, it is impossible to obtain explicit equation group, the method will be faced with many applications above it is difficult.Directly
Scoop out with Fast Fourier Transform (FFT) (FFT), using the extremely complex nonlinear system of numerical method analysis, then with very strong logical
The property used.But when the harmonic wave item number to be analyzed is many in the harmonic wave item number or response that system input/excitation is included, and in view of being
The subharmonic and ultraharmonics of system output, then be related to harmonic amplitude and a large amount of unknown numbers to be solved such as phase, have impact on the effect of analysis
Rate and precision, implement extremely complex, low precision, and the calculating time sharply increases, or even can not effectively analyze.
The content of the invention
Realize that complicated, analysis harmonic wave item number is few, time-consuming and precision is low during for harmonic Balance Analysis nonlinear system
Problem, it is an object of the invention to provide a kind of single-degree-of-freedom nonlinear system harmonic balance subtraction unit, to solve the above problems.
The object of the present invention is achieved like this：A kind of single-mode system harmonic balance subtraction unit, using a kind of harmonic wave
Balancing method, it considers the limited subharmonic in response when solving the response of single-degree-of-freedom nonlinear system, only, will be described limited
Subharmonic is divided into low-frequency harmonics part and high-frequency harmonic part, and only include average from the low-frequency harmonics part opens with first harmonic
Begin, by Optimization Method average, the harmonic amplitude of the low-frequency harmonics part and phase, estimate the high-frequency harmonic part
Harmonic amplitude and phase, and precision is improved by iteration, is continuously increased the most high order that the low-frequency harmonics part includes humorous
The number of times of ripple, increase by one every time, while removing the minimum subharmonic of high-frequency harmonic part, in terms of the low-frequency harmonics part
Obtained result finally solves average, the amplitude and phase of the limited subharmonic in response as initial value,
The single-mode system harmonic balance subtraction unit is included with lower component：
Nonlinear system arrangement components, expression formula and excitation for defining nonlinear system,
Symbolic variable definition component, for define average in response of nonlinear system, the amplitude of the limited subharmonic with
The symbolic variable of phase,
Calculating unit, for by the harmonic wave equilibrium method calculate response in average, the amplitude of the limited subharmonic with
Phase, output block, for storing result that the calculating unit obtains, result being shown in the way of image conversion,
The calculating unit includes：Computing module preparation module, computing module main body module, computing module post processing mould
Block, wherein, the computing module preparation module, for set the initial value of iteration, in a primitive period time dispersion number
Mesh, the computing module main body module, for calculating average, the amplitude of the limited subharmonic in response by the method for iteration
With phase, and evaluating, the computing module main body module is made up of two layers of circulation, and outer loop controls the low frequency humorous
The number of times of ripple part highest subharmonic, interior loop by Optimization Method respond described in average in low-frequency harmonics part,
The amplitude and phase of each harmonic；
An objective function unit is provided with the computing module main body module, is solved for constitution optimization described
Object function in low-frequency harmonics part when average, the amplitude of each harmonic and phase, in an iterative process, object function
The number of average, the amplitude of each harmonic and phase is in dynamic change in expression formula form, the low-frequency harmonics part to be optimized
Change, the objective function unit is exactly to be used for constructing this object function.
Further, applicable single-degree-of-freedom nonlinear system is expressed as equation
In formula, t is the time；Y (t) is the response of system, that is, is exported；Single order, the second order respectively responded
Derivative；For onWithFunction；F (t) is excitation, that is, is inputted；
In the calculating unit, the method bag that average, the amplitude of the limited subharmonic and phase are used in response is solved
Include following steps：
S1：The cyclic variable k of the highest subharmonic number of times of the low-frequency harmonics part, interior loop variable i is represented to put just
Value：Put initial value in k=1, i=1, the high-frequency harmonic part：yHF(t)=0, the initial value X0 of Optimization Solution is set；
S2：Optimization method calculates average, the amplitude of each harmonic and phase in the low-frequency harmonics part, including：
S21：Calculate the low-frequency harmonics part yLF(t) expression formula
In formula, A0/ 2 be average, AnWithThe respectively amplitude and phase of n-th harmonic, k is the low-frequency harmonics portion
The number of times of highest subharmonic in point；Ω1For the angular frequency of the angular frequency of fundamental wave, i.e. first harmonic,
S22：Calculate the response y (t) of nonlinear system expression formula (6)
Y (t)=yLF(t)+yHF(t) (6)
S23：Seek error in equation sequence e (tm) character expression, seek its FFT character expression
In formula, tmFor discrete time point, m=1,2,3 ..., N, N is the points of the error in equation series of discrete；
S24：Seek the expression formula of object function
S25：With formula (9) minimum optimization aim, obtain in the low-frequency harmonics part average, the amplitude of each harmonic with
Phase A0, A1 A2,..., Ak,Value；
S26：By the amplitude A of each harmonic in the low-frequency harmonics part1, A2..., AkIt is converted into nonnegative value, phaseBe transformed into it is interval [0,2 π) on, method is：
S3：Iterative initial value, high-frequency harmonic part are updated, including：
S31：Calculate the error in equation sequence e (tm) numerical expression, accounting equation error mean square root erms；
S32：Judge whether interior loop meets exit criteria：1. the error in equation root-mean-square value ermsIt is reduced to setting
Precision, 2. the iterations i of interior loop reach the upper limit of setting,
If meeting any of which condition, turning S4 increases the number of times k of the low-frequency harmonics part highest subharmonic, continues
Iteration,
If be all unsatisfactory for, turn S33；
S33：The high-frequency harmonic part is updated, the initial value X0 of nonlinear optimization is set, cyclic variable increases by 1, i=i+1,
Turn S22 and continue iteration, wherein, use equation below and step：
S331：The concrete numerical value of calculating formula (11),
In formula, yLFFor yLF(t) shorthand, yHFFor yHF(t) shorthand,For yLF(t) first derivative,For yHF(t) first derivative；
S332：The amplitude and phase of the harmonic wave included with FFT calculating formulas (11), are translated into cosine function and sinusoidal letter
The form of number sum, obtains amplitude anWith bn, seek the expression formula (12) of the estimate of the high-frequency harmonic part
S4：Judge whether the cyclic variable k for representing the highest subharmonic number of times of the low-frequency harmonics part is not more than to set
Upper limit number of times, or not up to precision conditions,
If it is, the initial value X0 of the nonlinear optimization is set, cyclic variable k increases by 1：K=k+1, goes to step S2；
Otherwise, calculating terminates.
Further, in step S4, upper limit number of times is not more than the points N of error in equation series of discrete described in step S2
Half N/2.
Further, in step S4, precision conditions are that outer loop previous cycle, i.e. kth wheel are circulated, with last round of circulation
The obtained error in equation root-mean-square value it is equal or difference be less than setting the limits of error.
In some embodiments, in step S3, the initial value X0 of nonlinear optimization method to set up is：Interior loop
In this i wheel circulations from 1 to i, when ith is circulated, by optimizing obtained average, harmonic amplitude and phase in step S2,
In step S4, the initial value X0 of the nonlinear optimization includes 2k+3 element, is this k wheel of the outer loop from 1 to k
In circulation, under kth time circulation, in corresponding interior loop, when ith is circulated, the result obtained in step S2 and S3, wherein 2k
+ 1 element is by optimizing obtained average, harmonic amplitude and phase in step S2；Other 2 elements are FFT in step S3
Calculate the amplitude and phase of obtained corresponding harmonic wave.
In some embodiments, in step S3, the initial value X0 of nonlinear optimization method to set up is：Interior loop
In this i wheel circulations from 1 to i, when ith is circulated, by optimizing obtained average, harmonic amplitude and phase in step S2,
In step S4, the initial value X0 of the nonlinear optimization includes 2k+3 element, is this k wheel of the outer loop from 1 to k
In circulation, under kth time circulation, in this i wheel circulations of the corresponding interior loop from 1 to i, minimum error mean square root ermsIt is right
The result answered, wherein 2k+1 element is by optimizing obtained average, harmonic amplitude and phase in step S2；Other 2 members
Element calculates the amplitude and phase of obtained corresponding harmonic wave for FFT in step S3.
In some embodiments, step S33 detailed step is：Judge the error in equation root-mean-square value ermsWhether
The value obtained less than last round of iteration, if it is not, then the initial value X0 of nonlinear optimization is added at random on the basis of original
Disturbance, cyclic variable increase by 1：I=i+1, turns S2 and continues iteration；Otherwise, S331 and S332 is performed, iterative initial value X0 is set to step
By optimizing obtained average, harmonic amplitude and phase in rapid S2, if the error in equation root-mean-square value ermsChanged with last round of
The difference for the error mean square root that generation obtains is less than the limits of error of setting, then in X0 along with a random perturbation,
In step S4, the initial value X0 of the nonlinear optimization includes 2k+3 element, is this k wheel of the outer loop from 1 to k
In circulation, under kth time circulation, in corresponding interior loop, when ith is circulated, the result obtained in step S2 and S3, wherein 2k
+ 1 element is by optimizing obtained average, harmonic amplitude and phase in step S2；Other 2 elements are FFT in step S3
Calculate the amplitude and phase of obtained corresponding harmonic wave.
Further, step S21~S24 symbolization computings, step S3 uses numerical operation.
Further, in the symbolic variable definition component, when being realized using perceptive construction on mathematics, represent equal
The symbolic variable of value, amplitude and phase is made up of prefix with digital two parts, in Part I before average, amplitude symbolic variable
Sew identical, the prefix of phase symbol variable is identical, Part II numeral is positive integer numeral, the determination method of digit is：
If the highest overtone order of limited subharmonic described in response is M, then the minimum value of the digit of Part II numeral isIf digit not enough, is above supplied with numeral 0.
Yet further, the method that the objective function unit is used comprises the following steps：
a：Expression formula (4) is converted into character string myObjFunStr1 with char functions；
b：Then character string myObjFunStr1 average is replaced with into character string with character string replacement function strrep：x
(1), amplitude replaces with character string：X (2), x (3) ..., x (k), x (k+1), phase replaces with character string：x(k+2)、x(k+
3) ..., x (2k), x (2k+1), obtain representing the character string type expression formula myObjFunStr2 of object function；
c：Before character string myObjFunStr2 plus character string '@(x) ', obtain new character string
myObjFunStr3；
d：Character string myObjFunStr3 is converted into object function with eval functions.
The present invention has significant advantage and beneficial effect compared with prior art：
(1) in the present invention, symbolization computing is combined with numerical operation, and computational accuracy is high, efficiency high.Symbolic operation
Precision is ensure that, numerical operation improves operational efficiency, reduces the calculating time.
(2) the limited subharmonic considered in the response of nonlinear system is divided into low-frequency harmonics portion by the method that the present invention is used
Point with high-frequency harmonic part, since low-frequency harmonics part only comprising average and first harmonic, only three are to be optimized unknown
Amount, nonlinear optimization is simple and easy to apply；Gradually increase low-frequency harmonics part includes the item number of harmonic wave, and one high order of increase is humorous every time
Ripple, introduces two unknown quantitys, and nonlinear optimization is from simple to complexity, it is ensured that optimization precision and efficiency.
(3) an objective function unit is provided with the device that the present invention is provided, makes non-linear objective function all the time
Unknown quantity to be optimized is only included, i.e., there was only 2k+1 unknown quantity per suboptimization, the unknown quantity without redundancy improves optimization effect
Rate, substantially reduce run time.
(4) in step S3, random perturbation is added when initial value X0 is set, it is to avoid repeat what is undergone during follow-up Optimized Iterative
Calculating process, and add with randomization the speed of iteration convergence, can significantly improve the operational precision of device, reduce it is whole
The number of times of body iteration.
(5) in step S4, when initial value X0 is set, FFT in step S3 is used for two newly-increased unknown quantitys to be optimized
The amplitude and phase of obtained corresponding subharmonic are calculated, in the case of no prior information, the reasons why having abundant thinks it most
Close to actual value, therefore the efficiency of next round nonlinear optimization can be improved, improve computational accuracy.
(6) symbolic variable Part II numeral positive integer digit groups as defined in digit of average, amplitude and phase are represented
Into, and digit associates with the finite harmonic number of times to be considered, it is ensured that and it is readable, succinct that objective function unit is realized
Property and convenience, using the present invention provide device when, it is only necessary to change highest overtone order be M value, conveniently make
With.
Brief description of the drawings
The structure composition figure for the device that Fig. 1 provides for the present invention；
Fig. 2 is the schematic diagram for the harmonic wave equilibrium method that apparatus of the present invention are used；
Fig. 3 is the flow chart for the harmonic wave equilibrium method that apparatus of the present invention are used；
Fig. 4 is objective function element method figure；
Fig. 5 is the average and harmonic amplitude figure of the embodiment of the present invention 1；
Fig. 6 is the average and harmonic amplitude figure of the embodiment of the present invention 2；
Fig. 7 is the error in equation root-mean-square value of the embodiment of the present invention 3 with iterations variation diagram；
Fig. 8 is the average and harmonic amplitude figure of the embodiment of the present invention 3.
Embodiment
Below in conjunction with accompanying drawing and preferred embodiment, to the architectural feature provided according to the present invention, embodiment and its
Effect, is described in detail as follows.
The invention provides a kind of single-mode system harmonic balance subtraction unit, structure is as shown in Figure 1.What this device was applicable
Single-degree-of-freedom nonlinear system is represented by equation below
In formula, t is the time；Y (t) is the response of system, that is, is exported, convenient following for statement, and some ground places are omitted
(t)；The single order that respectively responds, second dervative；For onWithLetter
Number；F (t) is excitation, that is, is inputted.
When solving the response of nonlinear system with harmonic wave equilibrium method, output/response is expressed as to the form of Fourier space,
Only consider limited subharmonic, the number of times of highest subharmonic in response is set to a fixed value M.
Then output is expressed as
In formula, A0/ 2 be average, A1、A2、A3、…、AM-1、AMWithRespectively M humorous
The amplitude and phase of ripple, Ω1For the angular frequency of fundamental wave.
If y is not accurate solution, just there is an error e (t) on equation (1) the equal sign left side and the right
Harmonic balance subtraction unit needs average, M harmonic amplitude, M harmonic phase in solution formula (2), altogether 2M+1
Individual unknown quantity.The single-mode system harmonic balance subtraction unit that the present invention is provided uses method as shown in Figure 2, will respond
All harmonic waves in formula (2) are divided into the low-frequency harmonics part and the high-frequency harmonic part.The low-frequency harmonics part is included
Harmonic wave item number, that is, highest subharmonic number of times, controlled with cyclic variable k.
The low-frequency harmonics are partially shown as
In formula, A0/ 2 be average；AnWithThe respectively amplitude and phase of n-th harmonic；K is the low-frequency harmonics part
The number of times of middle highest subharmonic；Ω1For the angular frequency of the angular frequency of fundamental wave, i.e. first harmonic.
The low-frequency harmonics part described by formula (4) correspond to 1 average, k harmonic amplitude, k phase, so institute
Low-frequency harmonics part is stated altogether comprising 2k+1 unknown quantity.
Then the high-frequency harmonic is partially shown as
The number of times for the minimum subharmonic that the higher hamonic wave part described by formula (5) is included is k+1, then the harmonic high frequency
The harmonic wave item number that ripple part is included is M-k.This M-k harmonic wave corresponds to M-k amplitude and M-k phase, so the high frequency
Harmonic includes altogether 2 (M-k) individual unknown quantitys.
During beginning, k=1, the low-frequency harmonics part only includes average and first harmonic, now there is average, first harmonic
Amplitude and 3 unknown quantitys of phase, this 3 unknown quantitys are asked by nonlinear optimization, the initial value of optimization includes three elements, can be with
Random generation is set according to the priori of nonlinear system.With the concrete numerical value of formula (3) accounting equation error, pass through FFT
Calculating obtains the amplitude of the high-frequency harmonic part and the estimate of phase, altogether 2 (M-1) individual amount.According to the above results, repeatedly
Iteration, improves average, the precision of 3 unknown quantitys of first harmonic amplitude and phase.
It is continuously increased the number of times for the highest subharmonic that the low-frequency harmonics part is included, the increase by 1 every time of k values.That is
Low-frequency harmonics portions per increases by one, while removing the minimum subharmonic of the high-frequency harmonic part, i.e. k+1 subharmonic
.Using outer loop as 1,2,3 ..., iteration result during k iterates, finally asked as the initial value of+1 iteration of kth
Solve average in response, the amplitude and phase of M harmonic waves.
As shown in figure 1, the single-mode system harmonic balance subtraction unit is included with lower component：
Nonlinear system arrangement components, expression formula and excitation for defining nonlinear system,
Symbolic variable definition component, for define average in response of nonlinear system, the amplitude of the limited subharmonic with
The symbolic variable of phase,
Calculating unit, for by the harmonic wave equilibrium method calculate response in average, the amplitude of the limited subharmonic with
Phase, output block, for storing result that the calculating unit obtains, result being shown in the way of image conversion,
The calculating unit includes：Computing module preparation module, computing module main body module, computing module post processing mould
Block,
Wherein, the computing module preparation module, for set the initial value of iteration, in a primitive period time from
Dissipate number, the computing module main body module, for calculating average in response, the limited subharmonic by the method for iteration
Amplitude and phase, and evaluating, the computing module main body module are made up of two layers of circulation, and outer loop control is described low
The number of times of frequency harmonic highest subharmonic, interior loop by Optimization Method respond described in low-frequency harmonics part
Value, the amplitude and phase of each subharmonic；
An objective function unit is provided with the computing module main body module, is solved for constitution optimization described
Object function in low-frequency harmonics part when average, the amplitude of each harmonic and phase, in an iterative process, object function
The number 2k+1 of average, the amplitude of each harmonic and phase is dynamic in expression formula form, the low-frequency harmonics part to be optimized
State changes, and the objective function unit is exactly to be used for constructing this object function.
Specifically, the single-mode system harmonic balance subtraction unit that provides of the present invention mainly by two layers of circulation come
Realize.The item number for the harmonic wave that the outer loop variable k controls low-frequency harmonics part is included；Internal memory cyclic variable i is controlled in k
In the case of certain, the precision of average, the phase of harmonic wave and amplitude that the low-frequency harmonics part is included is improved by iteration, together
When also improve the amplitude and the precision of phase estimation value of the harmonic wave that the high-frequency harmonic part is included.Specifically come institute, mainly include
Following steps.
S1：The cyclic variable k of the highest subharmonic number of times of the low-frequency harmonics part, interior loop variable i is represented to put just
Value：K=1, i=1, the high-frequency harmonic part yHF(t) initial value is put：yHF(t)=0, the initial value X0 of Optimization Solution is set, initially
X0 can be generated at random, can also be determined according to the priori of non-thread system.
S2：Average, the amplitude of each harmonic and phase in the low-frequency harmonics part are calculated by optimization method, including：
S21：Calculate the low-frequency harmonics part yLF(t) expression formula (4)；
S22：Calculate the response y (t) of nonlinear system expression formula (6)；
Y (t)=yLF(t)+yHF(t) (6)
S23：And asked with formula (3) discrete type error in equation sequence e (tm)
e(tm)=yLF(tm)+yHF(tm) (7)
In formula, tmFor discrete time point, m=1,2,3 ..., N, N is the points of the error in equation series of discrete,
And seek error in equation sequence e (tm) Fast Fourier Transform (FFT) expression formula, and seek FFT expression formula
S24：Seek the expression formula of non-linear objective function
S25：With formula (9) minimum optimization aim, obtain in the low-frequency harmonics part average, the amplitude of each harmonic with
Phase A0, A1,A2,..., Ak,Value；
S26：By the harmonic amplitude A in the low-frequency harmonics part1, A2..., AkIt is converted into nonnegative value, phaseBe transformed into it is interval [0,2 π) on, method is：
S3：Iterative initial value, high-frequency harmonic part are updated, including：
S31：Calculate the error in equation sequence e (tm), calculate its error mean square root erms；
S32：Judge whether interior loop meets exit criteria：1. the error mean square root ermsIt is reduced under setting
Limit, 2. iterations i reaches the upper limit of setting,
If meeting one of condition, turn the number of times k of S4 increase subharmonic, continue iteration,
If be all unsatisfactory for, turn S33；
S33：Update the high-frequency harmonic part, S35：The initial value X0 of nonlinear optimization, cyclic variable increase by 1 are set：I=
I+1, turns S22 and continues iteration, wherein, use equation below and step：
S331：The value of calculating formula (11),
In formula, yLFFor yLF(t) shorthand, yHFFor yHF(t) shorthand,For yLF(t) first derivative,For yHF(t) first derivative；
S332：The amplitude and phase of the harmonic wave then included with FFT calculating formulas (11), are translated into cosine function and sine
The form of function sum, obtains amplitude anWith bn, seek the expression formula of the estimate of high-frequency harmonic part
S4：Judge whether the cyclic variable k of the highest subharmonic number of times of low-frequency harmonics part described in table is not more than setting
Upper limit number of times, or not up to precision conditions,
If it is, the initial value X0 of the nonlinear optimization is set,
Cyclic variable k increases by 1：K=k+1, goes to step S2；Otherwise, calculating terminates.
In step S4, upper limit number of times is not more than the half N/2 of the points N of error in equation series of discrete described in step S2.
With the increase of k values, the result precision more and more higher that the system with one degree of freedom harmonic balance subtraction unit is obtained.Work as k
When increasing to a certain degree, the error mean square root erms keeps constant or reduction very slow, at this moment, it is believed that iteration is whole
Only.
Step S3 and S4 are respectively that the next round of internal layer and outer loop is circularly set iterative initial value.As shown in Fig. 2 step
X0 has 2k+1 element in S3, corresponding to average, the amplitude and phase of k low-order harmonics, is that the next round circulation of internal layer prepares
Optimize initial value；X0 has 2k+3 element in step S4, is the next of outer layer corresponding to average, the amplitude and phase of k+1 harmonic waves
Wheel circulation prepares optimization initial value.Outer loop increase by 1, it is meant that add a harmonic wave in the low-frequency harmonics part.Step S4
In 2 elements having more than the X0 in step S3 of X0 be the corresponding amplitude of this harmonic wave that the low-frequency harmonics part is increased newly
With phase.For the setting of iterative initial value X0 in step S3 and S4, there are following 3 kinds technical sides that are all practical, working well
Case.
Technical scheme 1：In step S3, the initial value X0 of nonlinear optimization method to set up is：Interior loop is from 1 to i
This i wheel circulations in, when ith is circulated, by optimizing obtained average, harmonic amplitude and phase in step S2.
In step S4, the initial value X0 of the nonlinear optimization includes 2k+3 element, is this k wheel of the outer loop from 1 to k
In circulation, under kth time circulation, in corresponding interior loop, when ith is circulated, the result obtained in step S2 and S3, wherein 2k
+ 1 element is by optimizing obtained average, harmonic amplitude and phase in step S2；Other 2 elements are FFT in step S3
Calculate the amplitude and phase of obtained corresponding harmonic wave.
Technical scheme 2：In step S3, the initial value X0 of nonlinear optimization method to set up is：Interior loop is from 1 to i
This i wheel circulations in, when ith is circulated, by optimizing obtained average, harmonic amplitude and phase in step S2.
In step S4, the initial value X0 of the nonlinear optimization includes 2k+3 element, is this k wheel of the outer loop from 1 to k
In circulation, under kth time circulation, in this i wheel circulations of the corresponding interior loop from 1 to i, minimum error mean square root ermsIt is right
The result answered, wherein 2k+1 element is by optimizing obtained average, harmonic amplitude and phase in step S2；Other 2 members
Element calculates the amplitude and phase of obtained corresponding harmonic wave for FFT in step S3.
Technical scheme 3：Step S33 detailed step is：Judge the error in equation root-mean-square value ermsWhether upper one is less than
The value that wheel iteration is obtained, if it is not, then the initial value X0 of nonlinear optimization adds random perturbation, circulation on the basis of original
Variable increase by 1：I=i+1, turns S2 and continues iteration；Otherwise, S331 and S332 is performed, iterative initial value X0 is set to lead in step S2
Average, harmonic amplitude and phase that optimization is obtained are crossed, if the error in equation root-mean-square value ermsObtained with last round of iteration
The difference of error mean square root is less than the limits of error of setting, then in X0 along with a random perturbation, available random function is given birth at random
Into.
In step S4, the initial value X0 of the nonlinear optimization includes 2k+3 element, is this k wheel of the outer loop from 1 to k
In circulation, under kth time circulation, in corresponding interior loop, when ith is circulated, the result obtained in step S2 and S3, wherein 2k
+ 1 element is by optimizing obtained average, harmonic amplitude and phase in step S2；Other 2 elements are FFT in step S3
Calculate the amplitude and phase of obtained corresponding harmonic wave.
In step s 2, formula (8) is with all containing 2k+1 for unknown quantity to be asked in (9), and the form of these formulas in itself
Change with nonlinear system equation (1), M and k change, and difference is huge；Therefore, in step S21~S24 using symbol
Number computing, can make full use of the mathematical operation functions, pole such as powerful symbolic operation function, the derivation of engineering software or mathematical software
The earth reduces the complexity of programming.Step S3 is to estimate the amplitude and phase of each harmonic of the high-frequency harmonic part
Position, arithmetic speed can be greatly enhanced using numerical operation.
The single-mode system harmonic balance subtraction unit that the present invention is provided can be soft in many high-level programming languages, engineering
Implement in part, mathematical software, such as MATLAB.During specific implementation, average, harmonic amplitude and phase in formula (2) have altogether
2M+1 amounts to be solved are, it is necessary to define corresponding variable to handle.In order to construct the convenience of object function, average, width are represented
The symbolic variable of value and phase is made up of prefix with digital two parts.Average, the prefix of amplitude symbolic variable are identical, phase symbol
The prefix of variable is identical；Part II numeral is positive integer numeral, and the determination method of digit is：If needing the response analyzed
The number of times of middle highest harmonic wave is M, then the minimum value of the digit of Part II numeral isIf digit is not enough,
It is then high-order to be supplied by numeral 0.
For example, average, the prefix of amplitude symbolic variable are An, the prefix of phase symbol variable is Phin.Work as M<When 100,Then the method for the symbolic variable defined in MATLAB is as follows：
If M=12, above method define 13 symbolic variable An01 for representing average and harmonic amplitude, An02,
An03 ..., An11, An12, An13, also define 13 symbolic variable Phin01 for representing corresponding with above-mentioned amplitude phase,
Phin02, Phin03 ..., Phin11, Phin12, Phin13, and also define two group numbers array An and Phin, all include
13 elements, the respectively symbolic variable of amplitude and phase.Because An01 is corresponding with average, symbolic variable Phin01 is to be
The convenience of matching program, makes the amplitude of correspondence subharmonic with phase variant subscript consistent and definition, herein and without using.
Based on notation defined above variable, the method that the objective function unit is used comprises the following steps：
a：Expression formula (4) is converted into character string myObjFunStr1 with char functions；
b：Then character string myObjFunStr1 average is replaced with into character string with character string replacement function strrep：x
(1), amplitude replaces with character string：X (2), x (3) ..., x (k), x (k+1), phase replaces with character string：x(k+2)、x(k+
3) ..., x (2k), x (2k+1), obtain representing the character string type expression formula myObjFunStr2 of object function；
c：Before character string myObjFunStr2 plus character string '@(x) ', obtain new character string
myObjFunStr3；
d：Character string myObjFunStr3 is converted into object function with eval functions.
Step b mainly uses strrep functions, and specific method is：
Represent average, amplitude and phase Part II numeral digit minimum value asA high position mend 0 be for
Ensure to occur without wrong replacement during strrep functions are replaced.If not fixing digit, mistake occurs.Example
Such as：If by variable-definition be An1, An2, An3 ..., An11, An12, An13, function strrep when A1 replaces with x (1),
An10, An11, An12, An13 mistake can be substituted for x (1) 0, x (1) 1, x (1) 2, x (1) 3.
Implement the present invention in perceptive construction on mathematics.Illustrate beneficial effects of the present invention below with 3 embodiments.
Embodiment 1：The equation of single-degree-of-freedom nonlinear system is
Assuming that the angular frequency of first harmonic is 1 in response, iterative initial value X0=[0,0,0] is set, only considered less than 32 times
Harmonic wave, i.e. M=32, the maximum times of interior loop are 5, sampling number N=64, and the limits of error is set to 1 × 10 in step S4-16,
The scope of random perturbation is [- 5 × 10-5, 5 × 10-5].Implement technical scheme 3.With the increase of k values, the error in equation is equal
Root value reduces rapidly, and error in equation mean-square value during preceding 6 circulations is followed successively by：
5.656854249492381e+000,2.496859179098087e-003,8.156122475628342e-008,
4.360844843951736e-012,4.360844843951736e-012,4.360844843951736e-012.
The amplitude of each harmonic is with phase as shown in figure 5, wherein 0 subharmonic corresponds to average, and 1~6 subharmonic is by step
S2 nonlinear optimization is obtained, and 7~32 subharmonic are obtained by the numerical value FFT in step S3.
From above-mentioned data, it is apparent that working as k>After 4, iteration result tends to be constant, and precision grade is 10e-12,
I.e. 10-12。
Embodiment 2：Excitation in embodiment 1 is revised as cos2t, it is assumed that the angular frequency of first harmonic is 2 in response, if
It is random generation, other specification be the same as Example 1 to put iterative initial value X0.Implement technical scheme 3.The amplitude and phase of each harmonic are such as
Shown in Fig. 6, wherein 0 subharmonic correspondence average, 1~6 subharmonic obtains by step S2 nonlinear optimization, 7~32 subharmonic by
Numerical value FFT in step S3 is obtained.Work as k>After 2, iteration result tends to be constant, and precision grade is 10-10。
Embodiment 3：Excitation in embodiment 1 is revised as cos3t, it is assumed that the angular frequency of first harmonic is 3 in response, if
Put iterative initial value X0 to generate to be random, only the harmonic wave of consideration less than 32 times, i.e. M=32, the maximum times of interior loop are 4, are adopted
Number of samples N=64,.Implement technical scheme 3.With the increase of k values, the error in equation root-mean-square value reduces rapidly, such as Fig. 7 institutes
Show, iteration result precision grade is 10-12.The amplitude of each harmonic is with phase as shown in figure 8, the corresponding average of wherein 0 subharmonic, 1
~6 subharmonic are obtained by step S2 nonlinear optimization, and 7~32 subharmonic are obtained by the numerical value FFT in step S3.
The device that the present invention is provided only needs to change the number of times of the M i.e. controllable harmonic wave of value, using convenient.Above-mentioned 3
The result of embodiment, illustrates that the device of the invention provided has convergence rate and very high precision quickly.
If analyzing subharmonic, it is only necessary to set the frequency of first harmonic according to the angular frequency of subharmonic is wanted.
Object function can also be changed to the harmonic wave root-mean-square value of error in equation, good effect is equally can reach.
It should be noted that what above-mentioned specific embodiment was merely exemplary, under the above-mentioned teaching of the present invention, this area
Technical staff can carry out various improvement and deformation on the basis of above-described embodiment, and these are improved or deformation all falls within this
In the protection domain of invention.It will be understood by those skilled in the art that specific descriptions above are intended merely to explain the mesh of the present invention
, it is not intended to limit the present invention.Protection scope of the present invention is limited by claim and its equivalent.
Claims (10)
1. a kind of single-mode system harmonic balance subtraction unit, it is characterised in that use a kind of harmonic wave equilibrium method, it is solving list
During the response of free degree nonlinear system, only consider the limited subharmonic in response, the limited subharmonic is divided into low frequency humorous
Ripple part and high-frequency harmonic part,
Only including average with first harmonic since the low-frequency harmonics part, pass through Optimization Method average, the low frequency
The harmonic amplitude and phase of harmonic, estimate the harmonic amplitude and phase of the high-frequency harmonic part, and are carried by iteration
High accuracy, is continuously increased the number of times for the highest subharmonic that the low-frequency harmonics part is included, every time increase by one, while removing height
The minimum subharmonic of frequency harmonic, obtained result is calculated as initial value using the low-frequency harmonics part, final to solve
Go out average, the amplitude and phase of the limited subharmonic in response,
The single-mode system harmonic balance subtraction unit is included with lower component：
Nonlinear system arrangement components, expression formula and excitation for defining nonlinear system,
Symbolic variable definition component, for defining average in response of nonlinear system, the amplitude and phase of the limited subharmonic
Symbolic variable,
Calculating unit, for calculating average, the amplitude and phase of the limited subharmonic in response by the harmonic wave equilibrium method,
Output block, for storing result that the calculating unit obtains, result being shown in the way of image conversion,
The calculating unit includes：Computing module preparation module, computing module main body module, computing module post-processing module,
Wherein, the computing module preparation module, for set the initial value of iteration, in a primitive period time dispersion number
Mesh,
The computing module main body module, average in responding, the width of the limited subharmonic are calculated for the method by iteration
Value and phase, and evaluating, the computing module main body module are made up of two layers of circulation, and outer loop controls the low frequency
The number of times of harmonic highest subharmonic, interior loop by Optimization Method respond described in low-frequency harmonics part
Value, the amplitude and phase of each harmonic；
An objective function unit is provided with the computing module main body module, the low frequency is solved for constitution optimization
Object function in harmonic when average, the amplitude of each harmonic and phase, in an iterative process, the expression of object function
The number of average, the amplitude of each harmonic and phase is in dynamic change, institute in formula form, the low-frequency harmonics part to be optimized
It is exactly to be used for constructing this object function to state objective function unit.
2. single-mode system harmonic balance subtraction unit according to claim 1, it is characterised in that applicable single-degree-of-freedom
Nonlinear system is expressed as equation
In formula, t is the time；Y (t) is the response of system, that is, is exported；The single order that respectively responds, second dervative；For onWithFunction；F (t) is excitation, that is, is inputted；
Characterized in that, in the calculating unit, solving what average, the amplitude of the limited subharmonic and phase in response were used
Method comprises the following steps：
S1：Represent the cyclic variable k of the highest subharmonic number of times of the low-frequency harmonics part, interior loop variable i and put initial value：k
Put initial value in=1, i=1, the high-frequency harmonic part：yHF(t)=0, the initial value X0 of Optimization Solution is set；
S2：Optimization method calculates average, the amplitude of each harmonic and phase in the low-frequency harmonics part, including：
S21：Calculate the low-frequency harmonics part yLF(t) expression formula
In formula, A0/ 2 be average, AnWithThe respectively amplitude and phase of n-th harmonic, k be in the low-frequency harmonics part most
The number of times of higher hamonic wave；Ω1For the angular frequency of the angular frequency of fundamental wave, i.e. first harmonic；
S22：Calculate the response y (t) of nonlinear system expression formula (6)
Y (t)=yLF(t)+yHF(t) (6)
S23：Seek error in equation sequence e (tm) character expression, seek its FFT character expression
In formula, tmFor discrete time point, m=1,2,3 ..., N, N is the points of the error in equation series of discrete；
S24：Seek the expression formula of object function
S25：With formula (9) minimum optimization aim, average, the amplitude of each harmonic and phase in the low-frequency harmonics part are obtainedValue；
S26：By the amplitude A of each harmonic in the low-frequency harmonics part1, A2..., AkIt is converted into nonnegative value, phaseBe transformed into it is interval [0,2 π) on, method is：
S3：Iterative initial value, the high-frequency harmonic part are updated, including：
S31：Calculate the error in equation sequence e (tm) numerical expression, accounting equation error mean square root erms；
S32：Judge whether interior loop meets exit criteria：1. the error in equation root-mean-square value ermsIt is reduced to the essence of setting
Degree, 2. the iterations i of interior loop reaches the upper limit of setting,
If meeting any of which condition, turning S4 increases the number of times k of the low-frequency harmonics part highest subharmonic, continues to change
In generation, if be all unsatisfactory for, turn S33；
S33：The high-frequency harmonic part is updated, the initial value X0 of nonlinear optimization, cyclic variable increase by 1 are set：I=i+1, turns
S22 continues iteration, wherein, use equation below and step：
S331：The concrete numerical value of calculating formula (11),
In formula, yLFFor yLF(t) shorthand, yHFFor yHF(t) shorthand,For yLF(t) first derivative,
For yHF(t) first derivative；
S332：The amplitude and phase of the harmonic wave included with FFT calculating formulas (11), be translated into cosine function and SIN function it
The form of sum, obtains amplitude anWith bn, seek the expression formula (12) of the estimate of the high-frequency harmonic part
S4：Judge whether the cyclic variable k for representing the highest subharmonic number of times of the low-frequency harmonics part is not more than the upper of setting
Number of times, or not up to precision conditions are limited,
If it is, the initial value X0 of the nonlinear optimization is set, cyclic variable k increases by 1：K=k+1, goes to step S2；
Otherwise, calculating terminates.
3. single-mode system harmonic balance subtraction unit according to claim 2, it is characterised in that in step S4, the upper limit
Number of times is not more than the half N/2 of the points N of error in equation series of discrete described in step S2.
4. single-mode system harmonic balance subtraction unit according to claim 3, it is characterised in that in step S4, precision
Condition is that outer loop previous cycle, i.e. kth wheel are circulated, equal with the error in equation root-mean-square value that last round of circulation is obtained
Or difference be less than setting the limits of error.
5. single-mode system harmonic balance subtraction unit according to claim 4, it is characterised in that
In step S3, the initial value X0 of nonlinear optimization method to set up is：In this i wheel circulations of the interior loop from 1 to i,
When ith is circulated, by optimizing obtained average, harmonic amplitude and phase in step S2,
In step S4, the initial value X0 of the nonlinear optimization includes 2k+3 element, is this k repeating query of outer loop from 1 to k
In ring, under kth time circulation, in corresponding interior loop, when ith is circulated, the result obtained in step S2 and S3, wherein 2k+1
Individual element is by optimizing obtained average, harmonic amplitude and phase in step S2；Other 2 elements are FFT meters in step S3
The amplitude and phase of obtained corresponding harmonic wave.
6. single-mode system harmonic balance subtraction unit according to claim 4, it is characterised in that
In step S3, the initial value X0 of nonlinear optimization method to set up is：In this i wheel circulations of the interior loop from 1 to i,
When ith is circulated, by optimizing obtained average, harmonic amplitude and phase in step S2,
In step S4, the initial value X0 of the nonlinear optimization includes 2k+3 element, is this k wheel circulation of the outer loop from 1 to k
In, under kth time circulation, in this i wheel circulations of the corresponding interior loop from 1 to i, minimum error mean square root ermsIt is corresponding
As a result, wherein 2k+1 element is by optimizing obtained average, harmonic amplitude and phase in step S2；Other 2 elements are
FFT calculates the amplitude and phase of obtained corresponding harmonic wave in step S3.
7. single-mode system harmonic balance subtraction unit according to claim 4, it is characterised in that
Step S33 detailed step is：Judge the error in equation root-mean-square value ermsWhether it is less than what last round of iteration was obtained
Value, if it is not, then the initial value X0 of nonlinear optimization adds random perturbation, cyclic variable increase by 1 on the basis of original：i
=i+1, turns S2 and continues iteration；Otherwise, S331 and S332 is performed, iterative initial value X0 is set to what is obtained in step S2 by optimization
Average, harmonic amplitude and phase, if the error in equation root-mean-square value ermsThe error mean square root obtained with last round of iteration
Difference be less than the limits of error of setting, then in X0 along with random perturbation,
In step S4, the initial value X0 of the nonlinear optimization includes 2k+3 element, is this k wheel circulation of the outer loop from 1 to k
In, under kth time circulation, in corresponding interior loop, when ith is circulated, the result obtained in step S2 and S3, wherein 2k+1
Element is by optimizing obtained average, harmonic amplitude and phase in step S2；Other 2 elements are FFT calculating in step S3
The amplitude and phase of obtained corresponding harmonic wave.
8. the single-mode system harmonic balance subtraction unit according to claim 5,6 or 7, it is characterised in that step S21~
S24 symbolization computings, step S3 uses numerical operation.
9. single-mode system harmonic balance subtraction unit according to claim 8, it is characterised in that the symbolic variable is determined
In adopted part, when being realized using perceptive construction on mathematics, the symbolic variable of average, amplitude and phase is represented by prefix and numeral
Two parts are constituted, and average, the prefix of amplitude symbolic variable are identical in Part I, and the prefix of phase symbol variable is identical, and second
Part number is positive integer numeral, and the determination method of digit is：If the highest harmonic wave of limited subharmonic described in response
Number is M, then the minimum value of the digit of Part II numeral isIf digit is not enough, above with numeral 0
Supply.
10. single-mode system harmonic balance subtraction unit according to claim 9, it is characterised in that the object function
The method that structural unit is used comprises the following steps：
a：Expression formula (4) is converted into character string myObjFunStr1 with char functions；
b：Then character string myObjFunStr1 average is replaced with into character string with character string replacement function strrep：X (1),
Amplitude replaces with character string：X (2), x (3) ..., x (k), x (k+1), phase replaces with character string：x(k+2)、x(k+3)、…、x
(2k), x (2k+1), obtains representing the character string type expression formula myObjFunStr2 of object function；
c：Before character string myObjFunStr2 plus character string '@(x) ', obtain new character string myObjFunStr3；
d：Character string myObjFunStr3 is converted into object function with eval functions.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510295236.7A CN104881394B (en) | 2015-06-03 | 2015-06-03 | Single-mode system harmonic balance subtraction unit |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510295236.7A CN104881394B (en) | 2015-06-03 | 2015-06-03 | Single-mode system harmonic balance subtraction unit |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104881394A CN104881394A (en) | 2015-09-02 |
| CN104881394B true CN104881394B (en) | 2017-08-18 |
Family
ID=53948889
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510295236.7A Active CN104881394B (en) | 2015-06-03 | 2015-06-03 | Single-mode system harmonic balance subtraction unit |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104881394B (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112924798B (en) * | 2021-02-08 | 2022-11-08 | 北京中电普华信息技术有限公司 | Power quality monitoring method and device and electronic equipment |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102636693A (en) * | 2012-05-04 | 2012-08-15 | 重庆大学 | Harmonic analysis algorithm combining fast Fourier transform (FFT) and nonlinear least square |
| CN103399203A (en) * | 2013-08-09 | 2013-11-20 | 重庆大学 | High-precision harmonic parameter estimation method based on composite iterative algorithm |
-
2015
- 2015-06-03 CN CN201510295236.7A patent/CN104881394B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102636693A (en) * | 2012-05-04 | 2012-08-15 | 重庆大学 | Harmonic analysis algorithm combining fast Fourier transform (FFT) and nonlinear least square |
| CN103399203A (en) * | 2013-08-09 | 2013-11-20 | 重庆大学 | High-precision harmonic parameter estimation method based on composite iterative algorithm |
Non-Patent Citations (4)
| Title |
|---|
| 单自由度迟滞阻尼系统的谐波平衡分析;王轲;《第八届全国振动理论及应用学术会议论文集》;20031130;第2003年卷;全文 * |
| 基于傅里叶和小波变换的电力系统谐波分析;龚黎明 等;《测试测量技术》;20101031;第2010年卷(第10期);全文 * |
| 幅值波动电力信号谐波分析的HHT-FFT方法;梅永 等;《计量技术》;20130430;第2013年卷(第4期);全文 * |
| 非线性振动系统的主共振的增量谐波平衡法;熊铁华 等;《武汉大学学报》;20020430;第35卷(第2期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104881394A (en) | 2015-09-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| He | Homotopy perturbation method with an auxiliary term | |
| Garcke et al. | Suboptimal feedback control of PDEs by solving HJB equations on adaptive sparse grids | |
| LaBryer et al. | A harmonic balance approach for large-scale problems in nonlinear structural dynamics | |
| Benner et al. | Generalised tangential interpolation for model reduction of discrete-time MIMO bilinear systems | |
| CN105574809B (en) | Electromagnetic transient simulation graphics processor parallel calculating method based on matrix exponetial | |
| Adams et al. | Base-2 expansions for linearizing products of functions of discrete variables | |
| JP5316433B2 (en) | Optimization processing program, method and apparatus | |
| Zhang et al. | Machine learning building blocks for real-time emulation of advanced transport power systems | |
| CN104881394B (en) | Single-mode system harmonic balance subtraction unit | |
| Ou et al. | Unsteady adjoint method for the optimal control of advection and Burger's equations using high-order spectral difference method | |
| Truhar et al. | Damping optimization in mechanical systems with external force | |
| Gou et al. | The method of lower and upper solutions for impulsive fractional evolution equations in Banach spaces | |
| Wang et al. | Two-sided projection methods for model reduction of MIMO bilinear systems | |
| Xu et al. | Disturbance-observer-based LQR control of singularly perturbed systems via recursive decoupling methods | |
| Jiang et al. | Transformation matrix for time discretization based on Tustin’s method | |
| Liu et al. | A novel linear algebra method for the determination of periodic steady states of nonlinear oscillators | |
| Zhang et al. | Monotone iterative method for retarded evolution equations involving nonlocal and impulsive conditions | |
| Cen et al. | Corrected L-type Method for Multi-singularity Problems Arising from Delay Fractional Equations | |
| Khan et al. | Parameters approach applied on nonlinear oscillators | |
| CN113158447A (en) | Large-step-length frequency-shift electromagnetic transient simulation method and system | |
| Zhang et al. | Fractional Rayleigh–Duffing-like system and its synchronization | |
| CN106227971A (en) | Dynamic derivative fast prediction technology based on harmonic wave equilibrium method | |
| JP2017194811A (en) | Frequency characteristic analysis device, method, and program | |
| Yao et al. | A comparative study of time-marching schemes for fluid-structure interactions | |
| Akgün et al. | Frequency response investigations of multi-input multi-output nonlinear systems using automated symbolic harmonic balance method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| EXSB | Decision made by sipo to initiate substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |