CN110514608A - A kind of unbiased esti-mator method of the kinetics rate constant based on spectrum - Google Patents
A kind of unbiased esti-mator method of the kinetics rate constant based on spectrum Download PDFInfo
- Publication number
- CN110514608A CN110514608A CN201910802085.8A CN201910802085A CN110514608A CN 110514608 A CN110514608 A CN 110514608A CN 201910802085 A CN201910802085 A CN 201910802085A CN 110514608 A CN110514608 A CN 110514608A
- Authority
- CN
- China
- Prior art keywords
- matrix
- rate constant
- kinetics
- spectrum
- kinetics rate
- 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.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
- G01N21/27—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands using photo-electric detection ; circuits for computing concentration
- G01N21/274—Calibration, base line adjustment, drift correction
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/17—Systems in which incident light is modified in accordance with the properties of the material investigated
- G01N21/25—Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
- G01N21/31—Investigating relative effect of material at wavelengths characteristic of specific elements or molecules, e.g. atomic absorption spectrometry
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
- G01N21/00—Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
- G01N21/84—Systems specially adapted for particular applications
- G01N2021/8411—Application to online plant, process monitoring
Landscapes
- Physics & Mathematics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Analytical Chemistry (AREA)
- Biochemistry (AREA)
- Pathology (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Chemical & Material Sciences (AREA)
- Immunology (AREA)
- General Physics & Mathematics (AREA)
- General Health & Medical Sciences (AREA)
- Theoretical Computer Science (AREA)
- Engineering & Computer Science (AREA)
- Mathematical Physics (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
- Investigating Or Analysing Materials By The Use Of Chemical Reactions (AREA)
Abstract
The unbiased esti-mator method for the kinetics rate constant based on spectrum that the present invention provides a kind of, belongs to Chemical Kinetics field.This method constructs spectroscopic data structure first, according to the initial estimate of reactant initial concentration and kinetics rate constant, calculates concentration profile, obtains spectrum matrix, compared with the spectrum matrix that measurement obtains, obtains error matrix；Then the estimated value that kinetics rate constant is adjusted according to error matrix, recalculates error matrix, iteration, until error reaches minimum value；Final error Matrix Estimation noise variance is finally used, increases noise in the spectrum that measurement obtains, recalculates kinetics rate constant estimated value, using unbiased esti-mator method, obtained result is more nearly true value.Due to using unbiased esti-mator method, the other informations such as the pure spectrum of reactive material are not depended on, the influence of noise generation is effectively reduced, the kinetics rate constant made is more accurate.
Description
Technical field
The invention belongs to Chemical Kinetics fields, normal more particularly to a kind of kinetics rate based on spectrum
Several unbiased esti-mator methods.
Background technique
Establish the kinetic model of accurate chemical reaction process, so as to chemical reaction process carry out reasonable prediction and
Effective control, both can guarantee the safety of entire production process, and can also improve Chemical Manufacture efficiency, and reduced cost, maximize
Economic benefit is of great significance for modern chemical industry.Accurate reaction Kinetics Model includes model structure and mould
Shape parameter.Establishing model structure not only needs to be well understood by the reaction mechanism of chemical reaction process, it is also necessary to according to reactor
Type, it is studied a question the characteristics of reasonably assumed.Equally significantly how Optimized model parameter so that mould
The quasi- data obtained of pattern and actual value deviation are minimum, are an emphasis of current Chemical Kinetics research.
By the development of many years, a large amount of scholars expand model parameter estimation from optimal setting and optimum control etc.
Research.Method based on spectrum simulation kinetic parameters also continues to develop, and proposes for the factor for influencing result precision
Corresponding strategy.Classical tracing analysis method is a kind of nonlinear least square fitting method, and this method is a kind of gradient side
Method iterates to optimal value from initial value according to gradient information.Can relatively rapid it be estimated from spectrum using this method
Reaction rate constant out, but experimental error and measurement noise are not conducted further research.In order to supplement the defect, benefit
It is also suggested with the Constraint least square algorithm of light spectrum estimation reaction rate constant, but this method is in the pure light for participating in reactive material
When spectrum is difficult to obtain, effect is had a greatly reduced quality.Therefore existing kinetics rate constant estimation method is there are still defect,
It needs to improve.
Summary of the invention
In order to solve deficiency existing for existing reactive kinetics parameters estimation method, the present invention proposes a kind of based on spectrum
Kinetics rate parameter unbiased esti-mator method.
The technical solution adopted by the present invention to solve the technical problems is: a kind of kinetics rate based on spectrum is normal
Several unbiased esti-mator methods, specifically includes the following steps:
Step 1: the structure based on Beer-Lambert law building spectroscopic data is established anti-using ordinary differential system
Kinetic model is answered, according to the initial concentration of reactant and kinetics rate constant initial estimate, calculates chemical reaction
Involved in the concentration profile that changes over time of all substances, utilize the absorptivity of each substance of least-squares calculation, meter
Spectrum data matrix is calculated, by the spectrum data matrix being calculated compared with the spectrum of experiment measurement, obtains error matrix E
(k0)。
Step 2: the error matrix E (k according to step 10) adjustment kinetics rate constant initial estimate, meter
Calculation obtains new concentration profile, using the absorptivity of each substance of least-squares calculation, calculates new spectrum data matrix,
By the new spectrum data matrix being calculated compared with the spectrum of experiment measurement, recalculate to obtain new error matrix E (k1)。
Through above method iteration, final error matrix E (k is obtainedn-1), until final error matrix E (kn-1) in each element
Quadratic sum and last time error matrix E (kn-2) in each element quadratic sum difference less than 10-6.N is the number of iterations, n >=1.
Step 3: the final error matrix E (k obtained according to step 2n-1) in each element quadratic sum estimation noise
Variance increases the noise with the variance in the spectrum for measuring and obtaining, obtains the spectroscopic data with superimposed noise,
Kinetics rate constant estimated value is recalculated, using unbiased esti-mator method, is completed to kinetics rate constant
Estimation.
Further, step 1 specifically:
The spectrum data matrix D are as follows: D=CST+E(kn-1), wherein C indicates the substance with optical absorption characteristics at any time
Between the concentration matrix that changes, S is absorptivity matrixes of the every kind of substance under different measurement wavelength, D indicate reaction solution for
The extinction degree of every kind of measurement wavelength changes with time, E (kn-1) identical as D dimension.
Reaction Kinetics Model is established by the ordinary differential system of independent variable of the time, gives spectrum data matrix D and institute
There is the initial concentration of reactant, according to the initial estimate k of kinetics rate constant0, calculate C, S, E (kn-1),
In, C+For the pseudoinverse of C:
E(kn-1)=D-CST
ST=C+D
C+=(CTC)-1CT。
Further, step 2 includes following sub-step:
(2.1)E(k0) in each element Ei,jQuadratic sum ssq solution procedure are as follows:
Nt is the quantity of sampled point, and n λ is the quantity for measuring wavelength.
(2.2) the shifted by delta k of kinetics rate constants k is calculated in a manner of mobile to its minimum value by ssq, reaction is dynamic
Mechanics rate constant initial-value table is shown as k0, by Δ k=(Δ k1,Δk2,…,Δknk) minor shifts after, use Taylor expansion
To estimate E (k0), wherein E (k0) Taylor expansion in contain partial derivative matrix.
(2.3) by E (k0) and the partial derivative matrix vector quantization, wherein the partial derivative matrix after vector quantization be known as it is refined can
Than matrix J, Δ k, e (k are solved by following equation0) it is E (k0) vector after vector quantization:
e(k0)=- J Δ k
J+=(JTJ)-1JT
Δ k=-J+e(k0)
(2.4) resulting Δ k and k will be calculated0It is added, obtains new kinetics rate constant estimated value, repeat to walk
Suddenly (2.1)-(2.3), until E (kn-1) ssq and with E (kn-2) ssq difference less than 10-6, loop termination obtains at this time anti-
Answer Kinetics Rate Constants By Using k (D).
Further, the step 3 includes following sub-step:
(3.1) E (k corresponding to rate constants k (D) obtained in step (2.4) is calculatedn-1) all elements that contain
Variances sigma2, it is considered as the variance of experiment measurement noise.
It (3.2) is 0 plus mean value on measure spectrum data matrix D, variance σ2White Gaussian noise, be denoted as Z, again
The estimated value k (Z) for calculating kinetics rate constant, obtaining unbiased esti-mator result is k=2k (D)-k (Z).
(3.3) step (3.2) are repeated, result is averaged, end reaction Kinetics Rate Constants By Using unbiased esti-mator knot is obtained
Fruit.
The beneficial effects of the present invention are: the present invention is realized dynamic to reaction from spectrum using classical tracing analysis method
The estimation of mechanics rate constant.Due to there is measurement noise in experiment, obtained result contains certain error, therefore uses nothing
Bias estimation does not depend on the other informations such as the pure spectrum of reactive material, is effectively reduced the influence of noise generation, makes to obtain
Kinetics rate constant it is more accurate.
Detailed description of the invention
Fig. 1 is the flow diagram of the method for the present invention；
Fig. 2 is the flow diagram of finite difference adaptive step selection method in step (2.2).
Specific embodiment
In order to which the purpose of the present invention and technical solution is more clearly understood, below in conjunction with drawings and examples, to this hair
It is bright to be further elaborated.It should be appreciated that the specific embodiments described herein are merely illustrative of the present invention, and do not have to
It is of the invention in limiting.
Referring to Fig.1, a kind of unbiased esti-mator method of the kinetics rate constant based on spectrum, comprising the following steps:
Step 1: the structure based on Beer-Lambert law building spectroscopic data is established anti-using ordinary differential system
Kinetic model is answered, according to the initial concentration of reactant and kinetics rate constant initial estimate, calculates chemical reaction
Involved in the concentration profile that changes over time of all substances, calculate all substances involved in chemical reaction and become at any time
The concentration profile of change obtains spectrum data matrix using the absorptivity of each substance of least-squares calculation, will be calculated
Spectrum matrix with experiment measurement spectrum compared with, obtain error matrix E (k0)。
Step 2: the error matrix according to step 1 adjusts kinetics rate constant initial estimate, calculates
New spectrum data matrix is obtained using the absorptivity of each substance of least-squares calculation to new concentration profile, will be calculated
Obtained new spectrum data matrix recalculates to obtain new error matrix E (k compared with the spectrum for measuring and obtaining1).Through
Above method iteration obtains final error matrix E (kn-1), until final error matrix E (kn-1) in each element
Quadratic sum and last time error matrix E (kn-2) in each element quadratic sum difference less than 10-6.N is the number of iterations, n >=1.
Step 3: the E (k obtained according to step 2n-1) in each element quadratic sum estimation noise variance, described
The noise for increasing in the spectrum of measurement and there is the variance is tested, the spectroscopic data with superimposed noise is obtained, is recalculated anti-
Kinetics Rate Constants By Using estimated value is answered, using unbiased esti-mator method, completes the estimation to kinetics rate constant.
Wherein, the method for step 1 specifically: according to Beer-Lambert law, when the time is t, in specific wavelength X
Place, the absorbance of solution are the extinction summations from all substances for absorbing the wavelength light.Tribute of the every kind of substance to absorbance
It offers with concentration and optical path length (incident beam advance in the solution at a distance from) linearly.If path length and ratio are normal
Number is by single constant si(λ) indicates that then Beer-Lambert law is indicated by following formula:
D (t, λ)=c1(t)s1(λ)+c2(t)s2(λ)+…+cnc(t)snc(λ)
Nc is the quantity in af at wavelength lambda extinction material, ci(t) be the time be t when i substance concentration, d (t, λ) indicate when
Between be t when, solution to wavelength be λ light absorption summation.The variation of chemical reaction is tracked, the spectrum of acquisition is the reaction time
Function, measurement result be two-dimensional spectrum data matrix: D=CST+E(kn-1), wherein C indicates that every kind of substance changes over time
Concentration, be nt × nc matrix；S is absorptivity of the every kind of substance under different measurement wavelength, is the matrix of n λ × nc；D
For nt × n λ spectrum data matrix, the extinction degree for indicating that reaction solution measures wavelength for every kind changes with time；Due to
There are errors for measurement, introduce E (kn-1), it is identical as D dimension.
The chemical reaction that the present invention studies, it is assumed that various materials uniformly mix in reactor, research object under the hypothesis
Model is established to the ordinary differential system (ordinary differential equations, ODES) that the time is independent variable
Reaction Kinetics Model gives spectrum data matrix D and all initial concentrations for participating in reactive material, according to kinetics speed
The initial estimate k of rate constant0, in conjunction with ODES, calculate C, E (kn-1), wherein C+For the pseudoinverse of C:
E(kn-1)=D-CST
ST=C+D
C+=(CTC)-1CT
Step 2 further includes following sub-step:
(2.1) using classical tracing analysis method, E (k is calculated0) in each element Ei,jQuadratic sum ssq (sum of
Squares, ssq) solution procedure are as follows:
Nt is the quantity of sampled point, and n λ is the quantity for measuring wavelength.
(2.2) the shifted by delta k of kinetics rate constants k is calculated in a manner of mobile to its minimum value by ssq, for this purpose,
It is emphasized that E (k0) and subsequent ssq be only the function of kinetics rate constants k, kinetics rate constant
Initial-value table is shown as k0, by Δ k=(Δ k1,Δk2,…,Δknk) minor shifts after, E (k is estimated with Taylor expansion0),
Wherein E (k0) Taylor expansion in contain partial derivative matrix, be used only first derivative, obtain following linear representation:
Although above-mentioned formula is approximate representation method, since it is linear representation, this is easily handled it.Target is
Make E (k0The value of+Δ k) tends to 0, therefore replaces E (k with 00+ Δ k), it is as follows that above formula rearranges result:
Error matrix E (k0) can be expressed by above formula, and partial derivative method can calculate according to the following formula:
Used herein is forward difference, carries out selection Δ k herein according to the Gill algorithm proposedi, algorithm flow reference
Fig. 2, step detailed annotation are as follows:
εA=| f (x) | εR
Wherein, h is finite difference adaptive step,For the forward difference of function f (x),For function f
(x) backward difference,For the centered difference of function f (x), Φ (f, h) is the second dervative of function f (x), εRFor machine
Device precision, εAFor minimum step,For the first discrimination variable,For the second discrimination variable,For third differentiation
Variable.
Step 1: confirmation function f (x) and independent variable x, selection parameter M enable m=0, hs=-1, Such as following formula:
Parameters in calculating formula (9)；IfThen hs←h0。
Step 2:m=m+1, hm=10hm-1, recalculate the parameters in formula (9).If hs< 0 andThen hs←hm.If metThen hΦ←hm, step 2 stops, into next step；
Otherwise step 2 is repeated until m=M, into next step；
Step 3:m=m+1, hm=hm-1/ 10, recalculate the parameters in formula (9).If
Then hs←hm.IfThen hΦ←hm, step 3 stops, into next step；Otherwise step 3 is repeated until m=M, step
Rapid 3 stop, into next step；
Step 4:hF=hΦ, export hF, algorithm termination；
Step 5: if hs< 0, then it takesIfThen take hF←hs；IfThen take
hF←hM.Export hF, algorithm termination.
(2.3) by E (k0) and the partial derivative matrix vector quantization, that is, it is launched into long column vector, the long column vector is rolled over
Build up the product of matrix and vector；Partial derivative matrix after vector quantization is known as Jacobian matrix J, solves Δ k by following equation,
e(k0) it is E (k0) vector after vector quantization:
e(k0)=- J Δ k
J+=(JTJ)-1JT
Δ k=-J+e(k0)
(2.4) resulting Δ k and k will be calculated0It is added, obtains new kinetics rate constant estimated value, repeat to walk
Suddenly (2.1)-(2.3), until E (kn-1) ssq and with E (kn-2) ssq difference less than 10-6, loop termination obtains at this time anti-
Answer Kinetics Rate Constants By Using k (D).
Step 3 further includes following sub-step:
(3.1) E (k corresponding to rate constants k (D) obtained in step (2.4) is calculatedn-1) all elements that contain
Variances sigma2, it is considered as the variance of experiment measurement noise.
It (3.2) is 0 plus mean value on measure spectrum data matrix D, variance σ2White Gaussian noise, be denoted as Z, again
The estimated value k (Z) for calculating kinetics rate constant, obtaining unbiased esti-mator result is k=2k (D)-k (Z).
(3.3) it is emphasized that the practical significance of unbiasedness refers to no systematic deviation.No matter measure same several times
One chemical reaction process, the reaction rate constant obtained by parameter Estimation for true value, always have it is some higher, one
It is a little relatively low.And unbiasedness indicates, on an average these errors, value zero, i.e. unbiased estimator only exist random error.Cause
This is repeated step (3.2) 150 times, and result is averaged, end reaction Kinetics Rate Constants By Using unbiased esti-mator result is obtained.
Embodiment
(1) according to Beer-Lambert law, when the time is t, in specific af at wavelength lambda, the absorbance of solution is to come from
The extinction summation of all substances for absorbing the wavelength light.Every kind of substance (enters the contribution of absorbance with concentration and optical path length
The distance that irradiating light beam is advanced in the solution) linearly.The variation of chemical reaction is tracked, the spectrum of acquisition is the reaction time
Function, measurement result is a two-dimensional matrix, as spectrum data matrix: D=CST+E(kn-1), wherein C indicates every kind of object
The concentration that matter changes over time, S are absorptivity of the every kind of substance under different measurement wavelength, and D is spectrum data matrix, is indicated
The extinction degree that reaction solution measures wavelength for every kind changes with time, E (kn-1) indicate error matrix, with D dimension phase
Together.
Reaction Kinetics Model is established by the ODES of independent variable of the time, selects the demonstration in MCR-ALS kit here
Example, ODES are as follows:
Reaction rate constant true value kr=[2,0.2], substance initial concentration [A]0=0.001molL-1, remaining substance
It is 0.Sampling time t=10s, t0=0, number of sampling points nt=300.Measure spectrum matrix D is it is known that reactive kinetics parameters
Initial estimate is set as k0=[1,0.5] can calculate the Matrix C that three kinds of material concentrations are distributed at any time in conjunction with ODES, accidentally
The calculation method of poor matrix E is as follows, wherein C+For the pseudoinverse of C:
E=D-CST
ST=C+D
C+=(CTC)-1CT
(2.1) using classical tracing analysis method, the error matrix E (k is calculated0) in each element Ei,jQuadratic sum
The solution procedure of ssq are as follows:
Nt is the quantity of sampled point, and n λ is the quantity for measuring wavelength.
(2.2) once it is determined that ssq, just calculates the offset of reaction rate constant by ssq to its minimum value in a manner of mobile.
For this reason, it may be necessary to it is emphasised that E (k0) and subsequent ssq be only the function of rate constants k, if rate constant initial-value table is shown as
k0, then parameter passes through Δ k=(Δ k1,Δk2,…Δknk) minor shifts after, E (k can be estimated with Taylor expansion0),
If first derivative is used only, following linear representation is obtained:
WhereinFor E (k0) for knkPartial derivative matrix.Target is to make E (k0The value of+Δ k) tends to 0, therefore with 0
Instead of E (k0+ Δ k), it is as follows that above formula rearranges result:
Error matrix E (k0) can be expressed by above formula, and partial derivative method can calculate according to the following formula:
(2.3) by E (k0) and the partial derivative matrix vector quantization, that is, it is launched into long column vector, the long column vector is rolled over
Build up the product of matrix and vector；Partial derivative matrix after vector quantization is known as Jacobian matrix J, solves Δ k by following equation,
e(k0) it is E (k0) vector after vector quantization:
e(k0)=- J Δ k
J+=(JTJ)-1JT
Δ k=-J+e(k0)
(2.4) resulting Δ k and k will be calculated0It is added, obtains new kinetics rate constant estimated value, repeat to walk
Suddenly (2.1)-(2.3), until E (kn-1) ssq and with E (kn-2) ssq difference less than 10-6, loop termination obtains at this time anti-
Answer Kinetics Rate Constants By Using.
In general, method fast convergence near minimum value as described above.However, if initial estimation is poor, Taylor's grade
The linearisation of the function approximation and problem of number expansion may cause the diverging of ssq and the failure of algorithm.Marquardt is based on
The modification that the thought of Levenberg proposes is to be added to a certain number of Marquardt parameter mp in parameter calculation procedure
On the diagonal element of Hessian matrix.Hessian matrix and Δ k it is new expression formula it is as follows:
H=JTJ
Δ k=- (H+mp*I)-1JTe(k0)
When the value of mp is noticeably greater than the diagonal element of Hessian matrix, Δ k expression formula is as follows, with linear decline
Property:
When mp very little, Δ k expression formula is constant.Management parameters mp strategy is as follows: 1 is initially set to, if parameter ssq is opened
Originate scattered, then mp is continuously increased (each iteration is multiplied by 10), and ssq convergence then constantly reduces (each iteration is divided by √ 10), repeatedly
1 is reset to after generation.
(3.1) rate constant obtained in step (2.4) is denoted as k (D)=[1.9766,0.19476], and it is right to calculate its institute
E (the k answeredn-1) variances sigma of all elements that contains2=4.3, the variance of experiment measurement noise can be considered as；
It (3.2) is 0 plus mean value in measure spectrum matrix D, variance σ2White Gaussian noise, be denoted as Z, recalculate
The estimated value k (Z) of kinetics rate constant, obtaining unbiased esti-mator result is k=2k (D)-k (Z).
(3.3) it repeats step (3.2) 150 times, is as a result averaged, can be regarded as a unbiased esti-mator and calculate.With true value kr=
It [2,0.2] is standard, average error rate can be used as the whether outstanding index of one result of evaluation.
Do not consider that measuring the resulting result of noise is k (D)=[1.9766,0.19476], average error rate 1.90%.
After having used unbiased esti-mator method, in 100 repetition experiments, there are 69 k (Z) true value more closer than k (D), part
As a result such as table 1.Illustrate in this example, unbiased esti-mator can more efficiently reduce noise bring error.
1 MCR-ALS Simulation Example unbiased esti-mator partial results of table
The above shows and describes the basic principles and main features of the present invention and the advantages of the present invention, but the present invention not
Be only limited to that specification and embodiments are listed to be used, it can be fully applied to various fields suitable for the present invention, for
For those skilled in the art, other modifications may be easily implemented, therefore without departing substantially from claim and equivalency range institute
Under the universal of limitation, the present invention is not limited to specific details and embodiment shown here.
Claims (4)
1. a kind of unbiased esti-mator method of the kinetics rate constant based on spectrum, it is characterised in that: specifically include following
Step:
Step 1: it is dynamic to establish reaction using ordinary differential system for the structure based on Beer-Lambert law building spectroscopic data
Mechanical model is calculated and is related in chemical reaction according to the initial concentration of reactant and kinetics rate constant initial estimate
And the concentration profile that changes over time of all substances, utilize the absorptivity of each substance of least-squares calculation, calculate light
Modal data matrix obtains error matrix E (k by the spectrum data matrix being calculated compared with the spectrum of experiment measurement0)。
Step 2: the error matrix E (k according to step 10) adjustment kinetics rate constant initial estimate, it calculates
New spectrum data matrix is calculated using the absorptivity of each substance of least-squares calculation to new concentration profile, will be counted
Obtained new spectrum data matrix recalculates to obtain new error matrix E (k compared with the spectrum of experiment measurement1).Through upper
Method iteration is stated, final error matrix E (k is obtainedn-1), until final error matrix E (kn-1) in each element square
With with last error matrix E (kn-2) in each element quadratic sum difference less than 10-6.N is the number of iterations, n >=1.
Step 3: the final error matrix E (k obtained according to step 2n-1) in each element quadratic sum estimation noise side
Difference increases the noise with the variance in the spectrum for measuring and obtaining, and obtains the spectroscopic data with superimposed noise, weight
New to calculate kinetics rate constant estimated value, using unbiased esti-mator method, kinetics rate constant is estimated in completion
Meter.
2. the unbiased esti-mator method of kinetics rate constant according to claim 1, which is characterized in that step 1 is specific
Are as follows:
The spectrum data matrix D are as follows: D=CST+E(kn-1), wherein C indicates that the substance with optical absorption characteristics changes over time
Concentration matrix, the absorptivity matrixes that S is every kind of substance under different measurement wavelength, D indicates that reaction solution surveys every kind
The extinction degree of amount wavelength changes with time, E (kn-1) identical as D dimension.
Reaction Kinetics Model is established by the ordinary differential system of independent variable of the time, gives spectrum data matrix D and all anti-
The initial concentration for answering object, according to the initial estimate k of kinetics rate constant0, calculate C, S, E (kn-1), wherein C+For
The pseudoinverse of C:
E(kn-1)=D-CST
ST=C+D
C+=(CTC)-1CT。
3. the unbiased esti-mator method of kinetics rate constant according to claim 1, which is characterized in that step 2 includes
Following sub-step:
(2.1)E(k0) in each element Ei,jQuadratic sum ssq solution procedure are as follows:
Nt is the quantity of sampled point, and n λ is the quantity for measuring wavelength.
(2.2) the shifted by delta k of kinetics rate constants k, kinetics are calculated in a manner of mobile to its minimum value by ssq
Rate constant initial-value table is shown as k0, by Δ k=(Δ k1,Δk2,…,Δknk) minor shifts after, estimated with Taylor expansion
Count E (k0), wherein E (k0) Taylor expansion in contain partial derivative matrix.
(2.3) by the E (k0) and the partial derivative matrix vector quantization, wherein the partial derivative matrix after vector quantization is known as Jacobi
Matrix J solves Δ k, e (k by following equation0) it is E (k0) vector after vector quantization:
e(k0)=- J Δ k
J+=(JTJ)-1JT
Δ k=-J+e(k0)
(2.4) resulting Δ k and k will be calculated0It is added, obtains new kinetics rate constant estimated value, repeat step
(2.1)-(2.3), until E (kn-1) ssq and with E (kn-2) ssq difference less than 10-6, loop termination obtains reaction at this time
Kinetics Rate Constants By Using k (D).
4. the unbiased esti-mator method of kinetics rate constant according to claim 1, which is characterized in that the step 3
Including following sub-step:
(3.1) E (k corresponding to rate constants k (D) obtained in step (2.4) is calculatedn-1) variance of all elements that contains
σ2, it is considered as the variance of experiment measurement noise.
It (3.2) is 0 plus mean value on measure spectrum data matrix D, variance σ2White Gaussian noise, be denoted as Z, recalculate
The estimated value k (Z) of kinetics rate constant, obtaining unbiased esti-mator result is k=2k (D)-k (Z).
(3.3) step (3.2) are repeated, result is averaged, end reaction Kinetics Rate Constants By Using unbiased esti-mator result is obtained.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910802085.8A CN110514608B (en) | 2019-08-28 | 2019-08-28 | Unbiased estimation method of reaction kinetic rate constant based on spectrum |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910802085.8A CN110514608B (en) | 2019-08-28 | 2019-08-28 | Unbiased estimation method of reaction kinetic rate constant based on spectrum |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110514608A true CN110514608A (en) | 2019-11-29 |
| CN110514608B CN110514608B (en) | 2021-08-24 |
Family
ID=68628676
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910802085.8A Active CN110514608B (en) | 2019-08-28 | 2019-08-28 | Unbiased estimation method of reaction kinetic rate constant based on spectrum |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110514608B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112818558A (en) * | 2021-02-22 | 2021-05-18 | 浙江工业大学 | Method for selecting estimable parameter subset aiming at mechanism modeling of reaction process |
Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7181323B1 (en) * | 2004-10-25 | 2007-02-20 | Lockheed Martin Corporation | Computerized method for generating low-bias estimates of position of a vehicle from sensor data |
| CN1916605A (en) * | 2006-08-31 | 2007-02-21 | 山东省科学院海洋仪器仪表研究所 | Method for analyzing kinetic parameters of ozone oxidation reaction through dynamics curve of chemiluminescence |
| CN101001219A (en) * | 2006-01-09 | 2007-07-18 | 电子科技大学中山学院 | Fast convergence rate adaptive blind estimation method for characteristic parameter |
| CN101218242A (en) * | 2005-04-01 | 2008-07-09 | 孟山都技术公司 | Control of N-(phosphonomethyl) iminodiacetic acid conversion in manufacture of glyphosate |
| CN101260037A (en) * | 2008-03-17 | 2008-09-10 | 华东理工大学 | Modeling method for alkylarene liquid phase oxidation dynamics mechanism model |
| CN101508768A (en) * | 2009-03-20 | 2009-08-19 | 华东理工大学 | Intelligent modeling method in industrial polyester production process |
| CN102744379A (en) * | 2012-03-07 | 2012-10-24 | 中冶南方工程技术有限公司 | Crystallizer control system state estimation method based on Kalman filtering |
| CN103105778A (en) * | 2013-02-05 | 2013-05-15 | 中原工学院 | Estimation method for industrial process simulation mathematical model parameters |
| US20160300036A1 (en) * | 2013-10-28 | 2016-10-13 | New York University | Methods, computer-accessible medium and systems to model disease progression using biomedical data from multiple patients |
| CN109270484A (en) * | 2018-07-24 | 2019-01-25 | 南京航空航天大学 | A kind of multiple source DOA estimation method based on movement integrated array |
-
2019
- 2019-08-28 CN CN201910802085.8A patent/CN110514608B/en active Active
Patent Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7181323B1 (en) * | 2004-10-25 | 2007-02-20 | Lockheed Martin Corporation | Computerized method for generating low-bias estimates of position of a vehicle from sensor data |
| CN101218242A (en) * | 2005-04-01 | 2008-07-09 | 孟山都技术公司 | Control of N-(phosphonomethyl) iminodiacetic acid conversion in manufacture of glyphosate |
| CN101001219A (en) * | 2006-01-09 | 2007-07-18 | 电子科技大学中山学院 | Fast convergence rate adaptive blind estimation method for characteristic parameter |
| CN1916605A (en) * | 2006-08-31 | 2007-02-21 | 山东省科学院海洋仪器仪表研究所 | Method for analyzing kinetic parameters of ozone oxidation reaction through dynamics curve of chemiluminescence |
| CN101260037A (en) * | 2008-03-17 | 2008-09-10 | 华东理工大学 | Modeling method for alkylarene liquid phase oxidation dynamics mechanism model |
| CN101508768A (en) * | 2009-03-20 | 2009-08-19 | 华东理工大学 | Intelligent modeling method in industrial polyester production process |
| CN102744379A (en) * | 2012-03-07 | 2012-10-24 | 中冶南方工程技术有限公司 | Crystallizer control system state estimation method based on Kalman filtering |
| CN103105778A (en) * | 2013-02-05 | 2013-05-15 | 中原工学院 | Estimation method for industrial process simulation mathematical model parameters |
| US20160300036A1 (en) * | 2013-10-28 | 2016-10-13 | New York University | Methods, computer-accessible medium and systems to model disease progression using biomedical data from multiple patients |
| CN109270484A (en) * | 2018-07-24 | 2019-01-25 | 南京航空航天大学 | A kind of multiple source DOA estimation method based on movement integrated array |
Non-Patent Citations (3)
| Title |
|---|
| KELLY L. FLEMING ET AL.: "New Approach for Investigating Reaction Dynamics and Rates with Ab Initio Calculations", 《THE JOURNAL OF PHYSICAL CHEMISTRY A》 * |
| 李婷婷 等: "气体No／s02／N2／02等离子体反应机理及动力学分析", 《工程热物理学报》 * |
| 陈迪钊 等: "基于连串反应在线光谱监测数据的动力学曲线拟合及速率常数的估计", 《计算机与应用化学》 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112818558A (en) * | 2021-02-22 | 2021-05-18 | 浙江工业大学 | Method for selecting estimable parameter subset aiming at mechanism modeling of reaction process |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110514608B (en) | 2021-08-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106443285B (en) | Multiple-harmonic-source harmonic responsibility quantitative analysis method based on total least square method | |
| CN108520325B (en) | Integrated life prediction method based on accelerated degradation data in variable environment | |
| CN103605860B (en) | A kind of imperfect light source analogy method based on equivalent phase screen method | |
| CN109088407B (en) | Power distribution network state estimation method based on deep belief network pseudo-measurement modeling | |
| CN110543618A (en) | roundness uncertainty evaluation method based on probability density function estimation | |
| CN103617816A (en) | Reactor core power distribution measuring method | |
| CN110826021B (en) | Robust identification and output estimation method for nonlinear industrial process | |
| CN114235111B (en) | Ultrasonic water meter flow calibration method based on model optimization | |
| CN104504475A (en) | AR*-SVM (support vector machine) hybrid modeling based haze time series prediction method | |
| CN108595892A (en) | Soft-measuring modeling method based on time difference model | |
| CN106610587A (en) | Temperature multi-model prediction function control method and device | |
| Salas et al. | Framework design for weight-average molecular weight control in semi-batch polymerization | |
| CN110514608A (en) | A kind of unbiased esti-mator method of the kinetics rate constant based on spectrum | |
| CN109698505B (en) | Regulation and control quantitative mapping calculation method for large power grid static voltage stability online prevention and control | |
| Kramer et al. | Computation of the posterior entropy in a Bayesian framework for parameter estimation in biological networks | |
| CN104865228A (en) | Quantitative laser-induced breakdown spectroscopy (LIBS) detecting method based on fusion entropy optimization | |
| CN115165770B (en) | Water COD and turbidity simultaneous detection method based on broad spectrum and BPNN | |
| CN108120694A (en) | For the multivariate calibration methods and system of Dark sun-cured chemical composition analysis | |
| KR20150010103A (en) | Method for quantitative and comparative analysis of distributions of molecular orbitals between different charge states and system using the same | |
| CN113702896B (en) | System and method for measuring direct-current electric energy standard meter error based on voltage reference | |
| Paturi et al. | Prediction of weld-line width and sink-mark depth of plastic injection moulded parts using neural networks | |
| CN107944552A (en) | A kind of industrial Internet of Things parameter prediction method based on Elman neutral nets | |
| CN110455902B (en) | Method for rapidly calibrating multiple standard samples in environment detection | |
| CN114487976A (en) | Method and system for evaluating traceability uncertainty of MCM electronic transformer calibrator | |
| Fikri et al. | Predicting Moroccan Real Network's Power Flow Employing the Artificial Neural Networks |
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 |