CN110514608B - Unbiased estimation method of reaction kinetic rate constant based on spectrum - Google Patents

Unbiased estimation method of reaction kinetic rate constant based on spectrum Download PDF

Info

Publication number
CN110514608B
CN110514608B CN201910802085.8A CN201910802085A CN110514608B CN 110514608 B CN110514608 B CN 110514608B CN 201910802085 A CN201910802085 A CN 201910802085A CN 110514608 B CN110514608 B CN 110514608B
Authority
CN
China
Prior art keywords
matrix
rate constant
reaction
spectrum
reaction kinetic
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
Application number
CN201910802085.8A
Other languages
Chinese (zh)
Other versions
CN110514608A (en
Inventor
陈伟锋
何腾
张贵军
胡俊
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University of Technology ZJUT
Original Assignee
Zhejiang University of Technology ZJUT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University of Technology ZJUT filed Critical Zhejiang University of Technology ZJUT
Priority to CN201910802085.8A priority Critical patent/CN110514608B/en
Publication of CN110514608A publication Critical patent/CN110514608A/en
Application granted granted Critical
Publication of CN110514608B publication Critical patent/CN110514608B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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

Abstract

The invention provides an unbiased estimation method of a reaction kinetics rate constant based on a spectrum, belonging to the field of chemical reaction kinetics. The method comprises the steps of firstly, constructing a spectrum data structure, calculating a concentration distribution curve according to initial concentration of reactants and an initial estimation value of a reaction kinetic rate constant, obtaining a spectrum matrix, and comparing the spectrum matrix with the spectrum matrix obtained by measurement to obtain an error matrix; then, adjusting the estimated value of the reaction dynamics rate constant according to the error matrix, recalculating the error matrix, and repeating iteration until the error reaches the minimum value; and finally, estimating the noise variance by using the final error matrix, adding noise in the measured spectrum, recalculating the reaction kinetic rate constant estimation value, and obtaining a result closer to a true value by using an unbiased estimation method. Due to the fact that an unbiased estimation method is used, other information such as a pure spectrum of a reaction substance is not relied on, the influence caused by noise is effectively reduced, and the obtained reaction kinetic rate constant is more accurate.

Description

Unbiased estimation method of reaction kinetic rate constant based on spectrum
Technical Field
The invention belongs to the field of chemical reaction kinetics, and particularly relates to an unbiased estimation method of a reaction kinetics rate constant based on a spectrum.
Background
The method has the advantages that an accurate dynamic model of the chemical reaction process is established so as to reasonably predict and effectively control the chemical reaction process, the safety of the whole production process can be guaranteed, the chemical production efficiency can be improved, the cost is reduced, the economic benefit is maximized, and the method has great significance for the modern chemical industry. The accurate reaction kinetics model includes model structure and model parameters. Establishing a model structure requires not only a good understanding of the reaction mechanism of the chemical reaction process, but also a reasonable assumption according to the type of the reactor and the characteristics of the problem under study. It is also very important how to optimize the model parameters so that the deviation of the data obtained by model simulation from the actual values is minimum, which is a key point of the current chemical reaction kinetics research.
Through the development of many years, a large number of scholars develop researches on model parameter estimation from the aspects of optimal setting, optimal control and the like. The method for fitting the kinetic model parameters based on the spectrum is also continuously developed, and a corresponding strategy is provided for factors influencing the result precision. The classical curve analysis method is a nonlinear least square fitting method, which is a gradient method and iterates to an optimal value according to gradient information from an initial value. The reaction rate constant can be estimated from the spectrum relatively quickly by using the method, but experimental errors and measurement noise are not further researched. To supplement this drawback, a constrained least squares method using spectra to estimate the reaction rate constant has also been proposed, but this method has been found to be a major compromise when pure spectra of the participating reactants are difficult to obtain. Therefore, the existing reaction kinetic rate constant estimation method still has defects and needs to be improved.
Disclosure of Invention
In order to overcome the defects of the existing reaction kinetic parameter estimation method, the invention provides an unbiased estimation method of a reaction kinetic rate parameter based on a spectrum.
The technical scheme adopted by the invention for solving the technical problems is as follows: an unbiased estimation method of a reaction kinetic rate constant based on a spectrum specifically comprises the following steps:
the method comprises the following steps: the method comprises the steps of constructing a structure of spectral data based on the Beer-Lambert law, establishing a reaction kinetic model by using an ordinary differential equation system, calculating concentration distribution curves of all substances related to chemical reaction along with time according to initial concentrations of reactants and initial estimation values of reaction kinetic rate constants, and calculating light absorption of each substance by using least squaresCalculating the coefficient, calculating the spectrum data matrix, comparing the calculated spectrum data matrix with the experimentally measured spectrum to obtain the error matrix E (k)0)。
Step two: said error matrix E (k) according to step one0) Adjusting the initial estimation value of the reaction kinetics rate constant, calculating to obtain a new concentration distribution curve, calculating the absorption coefficient of each substance by using least square, calculating a new spectral data matrix, comparing the calculated new spectral data matrix with the experimentally measured spectrum, and recalculating to obtain a new error matrix E (k)1). Repeating the iteration by the method to obtain a final error matrix E (k)n-1) Until the final error matrix E (k)n-1) The sum of squares of each element in the error matrix E (k) of the last timen-2) The difference between the square sums of each element in the series is less than 10-6. n is iteration times, and n is more than or equal to 1.
Step three: obtaining a final error matrix E (k) according to the step twon-1) The square sum of each element in the spectrum is used for estimating the variance of the noise, the noise with the variance is added in the measured spectrum to obtain the spectrum data with the superimposed noise, the estimation value of the reaction kinetic rate constant is recalculated, and the estimation of the reaction kinetic rate constant is completed by using an unbiased estimation method.
Further, the first step is specifically as follows:
the spectral data matrix D is: d ═ CST+E(kn-1) Wherein C represents a concentration matrix of a substance having a light absorption property with time, S is an absorption coefficient matrix of each substance at different measurement wavelengths, D represents a change in the degree of absorption of the reaction solution with time for each measurement wavelength, and E (k)n-1) The same as the D dimension.
Establishing a reaction kinetic model by using a time-independent ordinary differential equation set, giving a spectral data matrix D and initial concentrations of all reactants, and according to an initial estimation value k of a reaction kinetic rate constant0C, S, E (k) is calculatedn-1) Wherein, C+Is the pseudo-inverse of C:
E(kn-1)＝D-CST
ST＝C+D
C+＝(CTC)-1CT
further, the second step comprises the following substeps:
(2.1)E(k0) Each element E ofi,jThe solution for the sum of squares ssq is:
Figure BDA0002182608480000021
nt is the number of sampling points and n λ is the number of measurement wavelengths.
(2.2) calculating the shift Δ k of the reaction kinetic rate constant k, expressed as k, from ssq moving towards its minimum value0After Δ k ═ Δ k1,Δk2,…,Δknk) After a small offset, E (k) is estimated using Taylor expansion0) Wherein E (k)0) The Taylor expansion of (1) contains a partial derivative matrix.
(2.3) mixing E (k)0) And vectorizing the partial derivative matrix, wherein the vectorized partial derivative matrix is called a Jacobian matrix J, and solving delta k, e (k) by the following formula0) Is E (k)0) Vectorized vector:
e(k0)＝-JΔk
J+＝(JTJ)-1JT
Δk＝-J+e(k0)
(2.4) calculating the obtained delta k and k0Adding to obtain new reaction kinetic rate constant estimated value, repeating steps (2.1) - (2.3) until E (k)n-1) Ssq and E (k)n-2) Ssq is less than 10-6The cycle is terminated and the reaction kinetic rate constant k (D) is obtained at this point.
Further, the third step includes the following substeps:
(3.1) calculating E (k) corresponding to the rate constant k (D) obtained in the step (2.4)n-1) Variance σ of all elements contained2And the variance of the noise is taken as experimental measurement.
(3.2) adding the mean value of 0 and the variance of sigma to the measured spectral data matrix D2The white gaussian noise (k) is recorded as Z, and the estimated value k (Z) of the reaction kinetics rate constant is recalculated to obtain an unbiased estimation result k ═ 2k (d) -k (Z).
(3.3) repeating the step (3.2), and averaging the results to obtain the final unbiased estimation result of the reaction kinetic rate constant.
The invention has the beneficial effects that: the method starts from the spectrum, and estimates the reaction kinetic rate constant by using a classical curve analysis method. Because measurement noise exists in the experiment and the obtained result contains certain errors, an unbiased estimation method is used, other information such as a pure spectrum of a reaction substance is not relied on, the influence caused by the noise is effectively reduced, and the obtained reaction kinetic rate constant is more accurate.
Drawings
FIG. 1 is a schematic flow diagram of the process of the present invention;
fig. 2 is a schematic flow chart of the finite difference adaptive step size selection method in step (2.2).
Detailed Description
In order to make the purpose and technical solution of the present invention more apparent, the present invention will be further described in detail with reference to the accompanying drawings and examples. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.
Referring to fig. 1, a method for unbiased estimation of a reaction kinetic rate constant based on spectra, comprising the steps of:
the method comprises the following steps: the method comprises the steps of constructing a structure of spectral data based on the Beer-Lambert law, establishing a reaction kinetic model by using an ordinary differential equation system, calculating concentration distribution curves of all substances related to chemical reaction along with time according to initial concentrations of reactants and initial estimation values of reaction kinetic rate constants, calculating concentration distribution curves of all substances related to chemical reaction along with time, calculating light absorption coefficients of all substances by using least squares, and obtaining a spectral data matrixComparing the calculated spectrum matrix with the experimentally measured spectrum to obtain an error matrix E (k)0)。
Step two: adjusting the initial estimation value of the reaction kinetics rate constant according to the error matrix in the step one, calculating to obtain a new concentration distribution curve, calculating the absorption coefficient of each substance by using least square to obtain a new spectral data matrix, comparing the new spectral data matrix obtained by calculation with the spectrum obtained by measurement, and recalculating to obtain a new error matrix E (k)1). Repeating the iteration by the method to obtain a final error matrix E (k)n-1) Up to the final error matrix E (k)n-1) The sum of squares of each element in the error matrix E (k) of the last timen-2) The difference between the square sums of each element in the series is less than 10-6. n is iteration times, and n is more than or equal to 1.
Step three: e (k) obtained according to step twon-1) The square sum of each element in the spectrum is used for estimating the variance of the noise, the noise with the variance is added in the spectrum measured by the experiment, the spectrum data with the superimposed noise is obtained, the estimation value of the reaction kinetic rate constant is recalculated, and the estimation of the reaction kinetic rate constant is completed by using an unbiased estimation method.
The method of the first step specifically comprises the following steps: at time t, according to Beer-Lambert's law, the absorbance of a solution at a particular wavelength λ is the sum of the absorptions from all the species that absorb that wavelength of light. The contribution of each substance to absorbance scales linearly with concentration and path length (the distance traveled by the incident beam in solution). If the path length and the proportionality constant consist of a single constant si(λ), then the Beer-Lambert law is represented by the following equation:
d(t,λ)＝c1(t)s1(λ)+c2(t)s2(λ)+…+cnc(t)snc(λ)
nc is the number of light-absorbing substances at the wavelength λ, ci(t) is the concentration of the substance i at time t, and d (t, λ) is the sum of the absorptions of the solution at time t for light of wavelength λ. The change of the chemical reaction is tracked, the obtained spectrum is a function of the reaction time, and the measurement result isTwo-dimensional spectral data matrix: d ═ CST+E(kn-1) Wherein C represents the time-varying concentration of each substance and is a matrix of nt × nc; s is the absorption coefficient of each substance under different measurement wavelengths and is a matrix of n lambda multiplied by nc; d is a spectrum data matrix of nt multiplied by n lambda and represents the change of the light absorption degree of the reaction solution for each measurement wavelength along with the time; due to measurement errors, E (k) is introducedn-1) Same as the D dimension.
The chemical reaction of the invention is characterized in that various materials in a reactor are assumed to be uniformly mixed, a model of a study object under the assumption is used for establishing a reaction kinetic model by using an Ordinary Differential Equations (ODES) with time as an independent variable, a spectral data matrix D and the initial concentration of all substances participating in the reaction are given, and the initial estimation value k of the reaction kinetic rate constant is used0C, E (k) was calculated in conjunction with ODESn-1) Wherein, C+Is the pseudo-inverse of C:
E(kn-1)＝D-CST
ST＝C+D
C+＝(CTC)-1CT
the second step also comprises the following substeps:
(2.1) calculation of E (k) Using classical Curve analysis method0) Each element E ofi,jThe solving process for the sum of squares ssq (sum of squares, ssq) is:
Figure BDA0002182608480000051
nt is the number of sampling points and n λ is the number of measurement wavelengths.
(2.2) the shift Δ k of the reaction kinetic rate constant k is calculated in such a way that ssq moves towards its minimum, for which purpose it is emphasized that E (k)0) And subsequently ssq is a function of only the reaction kinetic rate constant k, the initial value of which is expressed as k0After Δ k ═ Δ k1,Δk2,…,Δknk) After a slight shift by TaylorExpansion to estimate E (k)0) Wherein E (k)0) The Taylor expansion of (1) contains a partial derivative matrix, and only uses a first derivative to obtain the following linear expression:
Figure BDA0002182608480000052
although the above formula is an approximate expression method, it is easy to handle since it is a linear expression. The goal is to make E (k)0The value of + Δ k) tends to 0, so E (k) is replaced by 00+ Δ k), the rearrangement result of the above equation is as follows:
Figure BDA0002182608480000053
error matrix E (k)0) Can be expressed by the above equation, and the partial derivative can be calculated according to the following method:
Figure BDA0002182608480000061
here, forward differencing is used, here selecting Δ k according to the algorithm proposed by GilliThe algorithm flow refers to fig. 2, and the steps are detailed as follows:
Figure BDA0002182608480000062
Figure BDA0002182608480000063
Figure BDA0002182608480000064
Figure BDA0002182608480000065
Figure BDA0002182608480000066
Figure BDA0002182608480000067
εA＝|f(x)|εR
wherein h is the finite difference self-adaptive step length,
Figure BDA0002182608480000068
for the forward difference of the function f (x),
Figure BDA0002182608480000069
for the backward difference of the function f (x),
Figure BDA00021826084800000610
is the central difference of the function f (x), phi (f, h) is the second derivative of the function f (x), epsilonRTo machine accuracy,. epsilonAIn order to be the minimum step size,
Figure BDA00021826084800000611
is used as a first judgment variable and is used as a second judgment variable,
Figure BDA00021826084800000612
is a second decision variable that is a function of the first decision variable,
Figure BDA00021826084800000613
is the third discriminant variable.
Step 1: identifying the function f (x) and the argument x, selecting the parameter M such that M is 0, hs＝-1，
Figure BDA00021826084800000614
Figure BDA00021826084800000615
The following formula:
Figure BDA00021826084800000616
calculating each parameter in the formula (9); if it is not
Figure BDA00021826084800000617
Then h iss←h0。
Step 2: m is m +1, hm＝10hm-1Then, each parameter in the equation (9) is recalculated. If h iss<0 and
Figure BDA00021826084800000618
then h iss←hm. If it is satisfied with
Figure BDA00021826084800000619
Then h isΦ←hmStopping step 2, and entering the next step; otherwise, repeating the step 2 until M is equal to M, and entering the next step;
and step 3: m is m +1, hm＝hm-1And/10, recalculating the parameters in equation (9). If it is not
Figure BDA00021826084800000620
Then h iss←hm. If it is not
Figure BDA00021826084800000621
Then h isΦ←hmStopping step 3, and entering the next step; otherwise, repeating the step 3 until M is equal to M, stopping the step 3, and entering the next step;
and 4, step 4: h isF＝hΦOutput hFThe algorithm terminates;
and 5: if h iss<0, then get
Figure BDA0002182608480000071
If it is not
Figure BDA0002182608480000072
Then take hF←hs(ii) a If it is not
Figure BDA0002182608480000073
Then take hF←hM. Output hFThe algorithm terminates.
(2.3) mixing E (k)0) And partial derivative matrix vectorization, namely unfolding long column vectors, and folding the long column vectors into a product of a matrix and the vectors; the partial derivative matrix after vectorization is called a Jacobian matrix J, and the following formula is used for solving delta k, e (k)0) Is E (k)0) Vectorized vector:
e(k0)＝-JΔk
J+＝(JTJ)-1JT
Δk＝-J+e(k0)
(2.4) calculating the obtained delta k and k0Adding to obtain new reaction kinetic rate constant estimated value, repeating steps (2.1) - (2.3) until E (k)n-1) Ssq and E (k)n-2) Ssq is less than 10-6The cycle is terminated and the reaction kinetic rate constant k (D) is obtained at this point.
The third step also comprises the following substeps:
(3.1) calculating E (k) corresponding to the rate constant k (D) obtained in the step (2.4)n-1) Variance σ of all elements contained2And the variance of the noise is taken as experimental measurement.
(3.2) adding the mean value of 0 and the variance of sigma to the measured spectral data matrix D2The white gaussian noise (k) is recorded as Z, and the estimated value k (Z) of the reaction kinetics rate constant is recalculated to obtain an unbiased estimation result k ═ 2k (d) -k (Z).
(3.3) it is emphasized that the practical meaning of unbiased means that there is no systematic deviation. No matter how many times the same chemical reaction process is measured, the reaction rate constants obtained through parameter estimation are higher and lower than the true values. And unbiased means that these errors are averaged and have a value of zero, i.e., only random errors exist in the unbiased estimate. Therefore, the step (3.2) is repeated 150 times, and the results are averaged to obtain the final unbiased estimation result of the reaction kinetic rate constant.
Examples
(1) At time t, according to Beer-Lambert's law, the absorbance of a solution at a particular wavelength λ is the sum of the absorptions from all the species that absorb that wavelength of light. The contribution of each substance to absorbance scales linearly with concentration and path length (the distance traveled by the incident beam in solution). Tracking the change of chemical reaction, obtaining a spectrum which is a function of reaction time, wherein the measurement result is a two-dimensional matrix which is a spectrum data matrix: d ═ CST+E(kn-1) Wherein C represents the concentration of each substance with time, S is the absorption coefficient of each substance at different measurement wavelengths, D is a spectral data matrix representing the change with time of the absorption degree of the reaction solution for each measurement wavelength, and E (k)n-1) Representing the error matrix, the same as the D dimension.
Establishing a reaction kinetic model by using ODES with time as an independent variable, wherein an example demonstration in an MCR-ALS toolkit is selected, and ODES is as follows:
Figure BDA0002182608480000081
Figure BDA0002182608480000082
true value of reaction rate constant kr＝[2,0.2]Initial concentration of substance [ A ]]0＝0.001mol·L-1And the balance 0. Sampling time t is 10s, t0The number of sampling points nt is 300. The initial estimate of the reaction kinetic parameter is set to k, given that the spectral matrix D is measured0＝[1,0.5]The matrix C of the distribution of the concentrations of the three substances over time can be calculated in combination with ODES, and the error matrix E can be calculated as follows, where C is+Is the pseudo-inverse of C:
E＝D-CST
ST＝C+D
C+＝(CTC)-1CT
(2.1) calculating the error matrix E (k) using classical curve analysis methods0) Each element E ofi,jThe solution for the sum of squares ssq is:
Figure BDA0002182608480000083
nt is the number of sampling points and n λ is the number of measurement wavelengths.
(2.2) once ssq is determined, the shift in reaction rate constant is calculated in such a way that ssq moves towards its minimum. For this reason, it is emphasized that E (k)0) And subsequently ssq is simply a function of the rate constant k if the initial value of the rate constant is denoted as k0The parameter is then changed by (Δ k) to (Δ k)1,Δk2,…Δknk) After a slight shift, the taylor expansion can be used to estimate E (k)0) If only the first derivative is used, the following linear expression is obtained:
Figure BDA0002182608480000084
wherein
Figure BDA0002182608480000085
Is E (k)0) For knkPartial derivative matrix of (a). The goal is to make E (k)0The value of + Δ k) tends to 0, so E (k) is replaced by 00+ Δ k), the rearrangement result of the above equation is as follows:
Figure BDA0002182608480000086
error matrix E (k)0) Can be expressed by the above equation, and the partial derivative can be calculated according to the following method:
Figure BDA0002182608480000087
(2.3) mixing E (k)0) And partial derivative matrix vectorization, namely unfolding long column vectors, and folding the long column vectors into a product of a matrix and the vectors; the partial derivative matrix after vectorization is called a Jacobian matrix J, and the following formula is used for solving delta k, e (k)0) Is E (k)0) Vectorized vector:
e(k0)＝-JΔk
J+＝(JTJ)-1JT
Δk＝-J+e(k0)
(2.4) calculating the obtained delta k and k0Adding to obtain new reaction kinetic rate constant estimated value, repeating steps (2.1) - (2.3) until E (k)n-1) Ssq and E (k)n-2) Ssq is less than 10-6The cycle is terminated and the reaction kinetic rate constant at this point is obtained.
Generally, the method converges rapidly around a minimum value as described above. However, if the initial estimate is poor, the functional approximation of the taylor series expansion and linearization of the problem may result in divergence of ssq and failure of the algorithm. The modification proposed by Marquardt based on the concept of Levenberg is to add a certain number of Marquardt parameters mp to the diagonal elements of the Hessian matrix during the parameter calculation. The new expression of the Hessian matrix and Δ k is as follows:
H＝JTJ
Δk＝-(H+mp*I)-1JTe(k0)
when the value of mp is significantly larger than the diagonal elements of the Hessian matrix, Δ k is expressed as follows, with the property of linear decline:
Figure BDA0002182608480000091
when mp is small, the Δ k expression is unchanged. The management parameter mp policy is as follows: initially set to 1, if parameter ssq begins to diverge, mp increases (multiply by 10 per iteration) and ssq converges and decreases (divide by v 10 per iteration), resetting to 1 after the iteration ends.
(3.1) the rate constant obtained in step (2.4) was recorded as k (d) ═ 1.9766,0.19476]Calculate the E (k) corresponding to itn-1) Variance σ of all elements contained24.3, the variance of the experimental measurement noise can be regarded as;
(3.2) adding the mean value 0 and the variance σ to the measured spectrum matrix D2The white gaussian noise (k) is recorded as Z, and the estimated value k (Z) of the reaction kinetics rate constant is recalculated to obtain an unbiased estimation result k ═ 2k (d) -k (Z).
And (3.3) repeating the step (3.2)150 times, averaging the results, and calculating as an unbiased estimation calculation. By true value kr＝[2,0.2]As a standard, the average error rate can be used as an index for evaluating whether a result is excellent.
The result obtained without taking into account the measurement noise was k (d) ═ 1.9766,0.19476, with an average error rate of 1.90%. After the unbiased estimation method was used, in 100 replicates, there were 69 times where k (z) was closer to true than k (d), and some of the results are shown in table 1. It is shown that in this example, the unbiased estimation can effectively reduce the error caused by the noise.
TABLE 1 partial results of unbiased estimation of MCR-ALS example simulation
Figure BDA0002182608480000101
While there has been shown and described what are at present considered to be the fundamental principles of the invention and its essential features and advantages, the invention is not limited to the details of the description and the embodiments, rather it is capable of modification in various other respects, all within the scope of the invention, and various modifications may be readily made by those skilled in the art without departing from the general concept defined by the claims and their equivalents.

Claims (2)

1. An unbiased estimation method of reaction kinetic rate constants based on spectra is characterized in that: the method specifically comprises the following steps:
the method comprises the following steps: the method comprises the steps of constructing a structure of spectral data based on the Beer-Lambert law, establishing a reaction kinetic model by using an ordinary differential equation system, calculating concentration distribution curves of all substances related to chemical reaction along with time according to initial concentrations of reactants and initial estimation values of reaction kinetic rate constants, calculating light absorption coefficients of all the substances by using least squares, calculating a spectral data matrix, and comparing the calculated spectral data matrix with experimentally measured spectra to obtain an error matrix
Figure DEST_PATH_IMAGE001
；
Step two: the error matrix according to step one
Figure 760067DEST_PATH_IMAGE001
Adjusting the initial estimation value of the reaction kinetics rate constant, calculating to obtain a new concentration distribution curve, calculating the absorption coefficient of each substance by using least square, calculating a new spectral data matrix, comparing the calculated new spectral data matrix with the experimentally measured spectrum, and recalculating to obtain a new error matrix
Figure DEST_PATH_IMAGE002
(ii) a Repeating the iteration by the method to obtain the final error matrix
Figure DEST_PATH_IMAGE003
Until the final error matrix
Figure 631202DEST_PATH_IMAGE003
The sum of squares of each element in the matrix and the last error matrix
Figure DEST_PATH_IMAGE004
The difference between the square sums of each element in the series is less than 10-6(ii) a n is iteration times, and n is more than or equal to 1;
the second step comprises the following substeps:
（2.1）
Figure 269994DEST_PATH_IMAGE001
each element of
Figure DEST_PATH_IMAGE005
The solution for the sum of squares ssq is:
Figure DEST_PATH_IMAGE006
，
Figure DEST_PATH_IMAGE007
is the number of sampling points that are to be sampled,
Figure DEST_PATH_IMAGE008
to measure the number of wavelengths;
(2.2) calculation of the reaction kinetic Rate constant in such a way that ssq is shifted towards its minimum
Figure DEST_PATH_IMAGE009
Is offset from
Figure DEST_PATH_IMAGE010
The initial value of the rate constant of the reaction kinetics is expressed as
Figure DEST_PATH_IMAGE011
Through which is passed
Figure DEST_PATH_IMAGE012
After a small offset, estimate by Taylor expansion
Figure 163738DEST_PATH_IMAGE001
Wherein
Figure 509269DEST_PATH_IMAGE001
The Taylor expansion contains partial derivative matrixes;
(2.3) mixingThe above-mentioned
Figure 66152DEST_PATH_IMAGE001
And vectorizing the partial derivative matrix, wherein the vectorized partial derivative matrix is called a Jacobian matrix
Figure DEST_PATH_IMAGE013
Solving by the following formula
Figure 696635DEST_PATH_IMAGE010
，
Figure DEST_PATH_IMAGE014
Is composed of
Figure 364377DEST_PATH_IMAGE001
Vectorized vector:
Figure DEST_PATH_IMAGE015
(2.4) will calculate
Figure 220206DEST_PATH_IMAGE010
And
Figure 885674DEST_PATH_IMAGE011
adding to obtain new reaction kinetic rate constant estimated value, repeating steps (2.1) - (2.3) until
Figure 797260DEST_PATH_IMAGE003
Ssp and
Figure 65430DEST_PATH_IMAGE004
ssp of less than 10-6At the end of the cycle, the rate constant of the reaction kinetics at that time is obtained
Figure DEST_PATH_IMAGE016
；
Step three: obtaining a final error matrix according to the step two
Figure 979029DEST_PATH_IMAGE003
The square sum of each element in the spectrum is used for estimating the variance of the noise, the noise with the variance is added in the measured spectrum to obtain the spectrum data with the superimposed noise, the estimation value of the reaction kinetics rate constant is recalculated, and the estimation of the reaction kinetics rate constant is completed by using an unbiased estimation method;
the third step comprises the following substeps:
(3.1) calculating the Rate constant obtained in step (2.4)
Figure 815398DEST_PATH_IMAGE016
Corresponding to
Figure 588182DEST_PATH_IMAGE003
Variance of all elements contained
Figure DEST_PATH_IMAGE017
The variance of the experimental measurement noise is regarded as;
(3.2) adding the mean value of 0 and the variance of 0 to the measured spectrum data matrix
Figure 548791DEST_PATH_IMAGE017
Is recorded as Z, recalculating the estimated value of the reaction kinetics rate constant
Figure DEST_PATH_IMAGE018
Obtaining an unbiased estimation result of
Figure DEST_PATH_IMAGE019
；
(3.3) repeating the step (3.2), and averaging the results to obtain the final unbiased estimation result of the reaction kinetic rate constant.
2. The unbiased estimation method of reaction kinetic rate constants according to claim 1, characterized in that step one is specifically:
the spectral data matrixDComprises the following steps:
Figure DEST_PATH_IMAGE020
wherein, in the step (A),Crepresenting a time-varying concentration matrix of a substance having light-absorbing properties,Sis a matrix of the absorption coefficients of each substance at different measurement wavelengths,Dshowing the change in the degree of absorption of the reaction solution for each measurement wavelength with time,
Figure 457841DEST_PATH_IMAGE003
andDthe dimensions are the same;
establishing a reaction kinetic model by using an ordinary differential equation set with time as an independent variable, and giving a spectrum data matrixDAnd initial concentrations of all reactants, based on initial estimates of reaction kinetic rate constants
Figure 747002DEST_PATH_IMAGE011
Calculate outC、S、
Figure 210345DEST_PATH_IMAGE003
Wherein, in the step (A),C +is composed ofCPseudo-inverse of (2):
Figure DEST_PATH_IMAGE021
。
CN201910802085.8A 2019-08-28 2019-08-28 Unbiased estimation method of reaction kinetic rate constant based on spectrum Active CN110514608B (en)

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 CN110514608A (en) 2019-11-29
CN110514608B true 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)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
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 (9)

* Cited by examiner, † Cited by third party
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
CN109270484A (en) * 2018-07-24 2019-01-25 南京航空航天大学 A kind of multiple source DOA estimation method based on movement integrated array

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015066052A1 (en) * 2013-10-28 2015-05-07 New York University Methods, computer-accessible medium and systems to model disease progression using biomedical data from multiple patients

Patent Citations (9)

* Cited by examiner, † Cited by third party
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
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)

* Cited by examiner, † Cited by third party
Title
New Approach for Investigating Reaction Dynamics and Rates with Ab Initio Calculations;Kelly L. Fleming et al.;《The Journal of Physical Chemistry A》;20151221;第120卷;第299-305页 *
基于连串反应在线光谱监测数据的动力学曲线拟合及速率常数的估计;陈迪钊 等;《计算机与应用化学》;20020728;第19卷(第4期);第399-408页 *
气体No／s02／N2／02等离子体反应机理及动力学分析;李婷婷 等;《工程热物理学报》;20070531;第28卷(第3期);第532-533页 *

Also Published As

Publication number Publication date
CN110514608A (en) 2019-11-29

Similar Documents

Publication Publication Date Title
Kelleher et al. Atomic transition probabilities of sodium and magnesium. A critical compilation
Wu et al. The QCD renormalization group equation and the elimination of fixed-order scheme-and-scale ambiguities using the principle of maximum conformality
Noor et al. Consistent estimation of Gibbs energy using component contributions
Cyburt et al. New BBN limits on physics beyond the standard model from 4He
Bhatt et al. Incremental identification of reaction systems—A comparison between rate-based and extent-based approaches
Bardow et al. Incremental and simultaneous identification of reaction kinetics: methods and comparison
Zong et al. Convergence and stability of the semi-tamed Euler scheme for stochastic differential equations with non-Lipschitz continuous coefficients
Zaika et al. Parametric identification of hydrogen permeability model by delay times and conjugate equations
Yurchenko et al. Ab initio dipole moment and theoretical rovibrational intensities in the electronic ground state of PH3
Donskoi et al. Optimization of coal pyrolysis modeling
CN110514608B (en) Unbiased estimation method of reaction kinetic rate constant based on spectrum
Wentworth et al. Rigorous least-squares estimation of molecular complex equilibriums. I. Single intermolecular complex utilizing spectrophotometric data
Jaumot et al. Noise propagation and error estimations in multivariate curve resolution alternating least squares using resampling methods
Opoku et al. Electron propagator self-energies versus improved GW100 vertical ionization energies
Lazopoulos et al. NLO QCD corrections to the production of t t¯ Z in gluon fusion
Kuksa et al. Factorization in the model of unstable particles with continuous masses
Ratcliffe Matrix approach to a numerical solution of the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi evolution equations
Gracey Four loop MS mass anomalous dimension in the Gross–Neveu model
Knoll et al. The thermodynamic scaling function of a polymer solution
CN114334030A (en) Method for evaluating high molecular polymerization reaction product based on quantum support vector machine
Messick et al. Selecting factors for partial least squares
Dijoux Construction of the tetration distribution based on the continuous iteration of the exponential‐minus‐one function
Liu et al. Parameter and uncertainty estimation in stochastic differential equation models with multi-rate data and nonstationary disturbances
Kharroubi et al. Posterior simulation via the signed root log-likelihood ratio
Pumputis et al. Estimation of parameters of finite population L-statistics

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