Machine learning combined online raman for real-time monitoring of hmx polymorphic transformation

Abstract

This invention provides a real-time monitoring method for HMX crystal transformation using machine learning combined with online Raman spectroscopy, belonging to the field of chemical substance analysis technology. The method includes the following steps: Crystal transformation is performed under controlled conditions using a stirred tank reactor; spectral data is acquired in real time using an online Raman spectrometer and combined with X-ray diffraction results to construct a dataset. A partial least squares quantitative model is further established by preprocessing the dataset using Savitzky-Golay first derivatives and optimizing variables using an adaptive gray wolf optimization algorithm. This model is established using a five-fold cross-validation method and used to predict the performance of validation set samples, with the coefficient of determination R0... 2 Using mean absolute error (MAE) as the evaluation metric, the accuracy of the quantitative model is effectively improved. This invention employs a machine learning-based real-time monitoring method combined with online Raman spectroscopy for HMX crystal transformation, which has significant application value for real-time monitoring of crystal transformation processes.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of chemical substance analysis, in particular to a machine learning combined with online Raman HMX crystallization real-time monitoring method. BACKGROUND

[0002] HMX has excellent explosive performance and is widely used in national defense and civil fields. The performance of HMX is largely dependent on the crystal structure (alpha, beta, gamma, delta). At present, the main product obtained by acetic anhydride method is alpha crystal HMX crude product, which has relatively poor stability and needs to be converted into beta-HMX to obtain beta-HMX with high density, good stability and more practical value.

[0003] In industrial production, the conversion of alpha-HMX crude product into beta-HMX is a key step to improve its performance. The traditional crystallization process monitoring mainly relies on offline sampling analysis by X-ray diffraction technology, which can provide accurate analysis results, but has the following obvious disadvantages:

[0004] 1) Periodic sampling from the crystallization process may cause sample discontinuity and representativeness problems;

[0005] 2) Offline analysis is time-consuming and may become a bottleneck in production environments requiring rapid response;

[0006] 3) Frequent sampling and analysis not only increase the work intensity, but also increase the overall production cost.

[0007] Therefore, it is of great significance to develop a method for real-time monitoring of HMX crystallization process.

[0008] Online Raman spectroscopy technology has become an ideal real-time monitoring tool with its non-invasive, high sensitivity and ability to provide molecular structure information. This method does not require complex sample processing and can be directly used in the production process, avoiding most steps in traditional methods. However, so far there has been no report on using machine learning combined with online Raman to monitor HMX crystallization process in real time. SUMMARY

[0009] The purpose of the present application is to provide a machine learning combined with online Raman HMX crystallization real-time monitoring method, which has important application value for real-time monitoring of crystal transformation process.

[0010] To achieve the above purpose, the present application provides a machine learning combined with online Raman HMX crystallization real-time monitoring method, comprising the following steps:

[0011] S1, taking alpha-HMX solid and completely dissolving it in a mixed solvent, adding it to a stirred tank reactor and starting stirring after heating to a constant temperature;

[0012] The β-HMX is taken as an inducing crystal seed, added into a stirred tank reactor to form a mixed solution, and kept warm to convert the crystal;

[0013] During the crystal conversion process, the spectrum is collected every 3 minutes by using an online Raman spectrometer, and 1% of the mixed solution is taken, filtered, dried, and then the content of β-HMX is determined by using an X-ray diffractometer;

[0014] S2, the collected Raman spectrum and the corresponding chemical value matrix of the content of β-HMX determined by the X-ray diffractometer are composed of a data set, and a calibration set sample and a validation set sample are obtained according to a random division method;

[0015] S3, the spectrum matrix of the calibration set sample and the spectrum matrix of the validation set sample are pretreated by Savitzky-Golay first-order derivative;

[0016] S4, the spectrum matrix of the calibration set sample is selected by using an adaptive grey wolf optimization algorithm, and the extracted key spectral variables are used to screen data from the spectrum matrix of the calibration set sample and the spectrum matrix of the validation set sample, to obtain an optimized calibration set sample spectrum matrix subset and a validation set sample spectrum matrix subset;

[0017] S5, the optimized calibration set sample spectrum matrix subset and the corresponding concentration matrix are used to establish a PLS quantitative model by using a five-fold cross-validation method of partial least squares, and the prediction performance of the PLS quantitative model is verified by using the validation set sample spectrum matrix subset and the corresponding concentration matrix;

[0018] S6, the newly collected Raman spectrum data are screened by the Savitzky-Golay first-order derivative of step S3 and the key spectral variables obtained by the adaptive grey wolf optimization algorithm in step S4 to obtain a prediction set sample spectrum matrix subset, and the PLS quantitative model established in step S5 is used to monitor the HMX crystal conversion process in real time.

[0019] Preferably, in step S1, the ratio of α-HMX solid to mixed solvent is 1g:5-20ml.

[0020] Preferably, in step S1, the mixed solvent is composed of acetone, ethyl acetate and water, and the ratio of acetone, ethyl acetate and water is 10-80:10-80:10-80, and the total of the three is 100%.

[0021] Preferably, in step S1, the stirring rate is 100-300 rpm / min; the constant temperature is 30-70℃; and the mass of β-HMX added is 5-20% of the mass of α-HMX solid.

[0022] Preferably, in step S3, the specific steps of Savitzky-Golay first derivative are as follows:

[0023] S31, select appropriate window size w, polynomial degree p and difference term d, where d = 1;

[0024] S32, construct a design matrix V, where each row corresponds to the Raman spectrum data of a given sample x k (x1, x2, x3……x n ) and each column corresponds to different power terms and derivative terms:

[0025]

[0026] S33, calculate the filter coefficient E by least squares method:

[0027] E = (V T V) -1 V T Y

[0028] Where Y is the matrix to be solved:

[0029] S34, use the filter coefficient E to calculate the output matrix Y' of the filtered given sample x k :

[0030] Y' = VE

[0031] S35, repeat steps S31 to S34, that is, obtain the corrected set sample spectrum matrix X1 and the verification set sample spectrum matrix X2 after Savitzky-Golay first derivative processing.

[0032] Preferably, in step S4, the specific steps of adaptive grey wolf optimization algorithm are as follows:

[0033] S41, initialize parameters: determine the dimension dim of the optimization problem, the number of grey wolf population N in the algorithm, the maximum number of iterations T, and set the upper bound ub and lower bound lb for each dimension;

[0034] The position S i of each individual in the population is randomly initialized in the search space, i = 1, 2, 3, …, N, that is:

[0035] S i = lb + (ub-lb)·rand(dim)

[0036] S42, define the fitness function: for each individual in the population, calculate its fitness value Fitness i under the fitness function f:

[0037]

[0038] where y j is the true value of the β-HMX content of the jth sample in the calibration set; y j is the predicted value of the β-HMX content of the jth sample in the calibration set, and n1 is the number of samples in the calibration set.

[0039] S43, coefficient dynamic updating mechanism: the updating formulae of the aggregation coefficient a and the step size g are respectively:

[0040]

[0041] where a new is the new aggregation coefficient used in the current iteration; b is the attenuation coefficient; t is the current iteration number; g new is the new step size used in the current iteration; mean Fitness is the average fitness of the population; std Fitness is the standard deviation of the fitness of the population.

[0042] S44, position updating: through the following position updating formula, the current optimal individual position S best is combined with the variable A, the variable C, the variable D, and the position S i of each individual to calculate the new position of each individual

[0043] A = 2 · a new · rand(dim) - a new

[0044] C = 2 · rand(dim)

[0045] D = |C · S best '- S i |

[0046]

[0047] where N(0, 1) represents a standard normal distribution.

[0048] S45, boundary check and fitness updating: the boundary check is performed on the new position, and if it exceeds the search space, it is corrected; then, the fitness value is recalculated based on the new position.

[0049] S46, repeating steps S41-S45 until the maximum iteration number T is reached.

[0050] S47, outputting the global optimal fitness value and its optimal position S best , wherein S best is the key spectral variable extracted.

[0051] Therefore, the application adopts the above-mentioned machine learning combined with online Raman HMX crystallization real-time monitoring method, and the technical effects are as follows:

[0052] (1) The machine learning combined with online Raman realizes real-time monitoring of the HMX crystallization process, the method is simple to operate, cost-effective, reduces the technical requirements for the operator, reduces the use of expensive instruments, and is easy to be widely used in the field of chemical material production.

[0053] (2) The application of the self-adaptive grey wolf optimization algorithm optimizes the variable selection process of the model, improves the prediction accuracy and generalization ability of the model, and reduces the complexity of the model. BRIEF DESCRIPTION OF DRAWINGS

[0054] Figure 1 is a flow chart of the machine learning combined with online Raman HMX crystallization real-time monitoring method of the application;

[0055] Figure 2 is a comparison chart of the real-time monitoring of the content of β-HMX in the HMX crystallization process and the offline analysis results of the X-ray diffractometer. DETAILED DESCRIPTION

[0056] The technical solutions of the application are further described below through the drawings and examples.

[0057] Unless otherwise defined, the technical terms or scientific terms used in the application should be understood as the usual meanings understood by those skilled in the art in the field to which the application belongs.

[0058] Example 1

[0059] As shown in Figure 1 is a flow chart of the machine learning combined with online Raman HMX crystallization real-time monitoring method of the application, which specifically includes the following steps:

[0060] S1, in a stirred tank reactor with a jacket, 60 g of α-HMX crude product is dissolved in 300 mL of mixed solvent (acetone, ethyl acetate and water in a volume ratio of 5:4:1) to prepare a saturated solution; under a stirring speed of 115 rpm / min, the temperature is raised to 55℃;

[0061] After the temperature in the kettle is constant, 3.6 g of β-HMX is taken as an induction crystal seed and added to the above-mentioned saturated solution to form a mixed suspension, and under the elapse of time, the transformation of α-HMX to β-HMX is realized.

[0062] During the whole crystallization process, the spectra were collected by blaze online Raman spectrometer every 3 minutes (a total of 220 Raman spectra), and a small part of the above suspension was taken, filtered, dried, and then the β-HMX content was determined by X-ray diffractometer;

[0063] S2, the spectrum matrix of the 220 collected Raman spectra and the corresponding β-HMX content chemical value matrix determined by the X-ray diffractometer constitute a data set, and are randomly divided according to a ratio of 7:3 to obtain 154 calibration set samples and 66 validation set samples;

[0064] S3, the spectrum matrix of the 154 obtained calibration set samples and 66 validation set samples is pretreated by Savitzky-Golay first derivative;

[0065] The specific steps of Savitzky-Golay first derivative are as follows:

[0066] a) selecting a suitable window size w, polynomial degree p and difference term d (d = 1);

[0067] b) constructing a design matrix V, wherein each row corresponds to the Raman spectrum data of a given sample x k (x1,x2,x3……x n ), and each column corresponds to different power terms and derivative terms:

[0068]

[0069] c) calculating the filter coefficient E by least squares method, wherein Y is the matrix to be solved:

[0070] E=(V T V) -1 V T Y

[0071] d) for a given sample x k , using the filter coefficient E to calculate the output matrix Y' of the filtered given sample:

[0072] Y′=VE

[0073] e) repeating steps a to d, that is, obtaining the calibration set sample spectrum matrix X1 and the validation set sample spectrum matrix X2 after Savitzky-Golay first derivative processing.

[0074] In specific use, the window size w and the polynomial degree p are 50 and 2 respectively.

[0075] S4, using the adaptive grey wolf optimization algorithm to select variables from the spectrum matrix of the calibration set samples, and using the extracted key spectral variables to screen data from the spectrum matrix of the calibration set samples and the spectrum matrix of the validation set samples, to obtain an optimized calibration set sample spectrum matrix subset and a validation set sample spectrum matrix subset;

[0076] The specific steps of the adaptive grey wolf optimization algorithm are as follows:

[0077] a) Initialization parameters: determine the dimension dim of the optimization problem, the number of grey wolf populations N in the algorithm, the maximum number of iterations T, and set the upper bound ub and the lower bound lb for each dimension.

[0078] The position S of each individual in the population i (i = 1, 2, 3, …, N) is randomly initialized in the search space, that is:

[0079] S i = lb + (ub-lb)·rand(dim)

[0080] b) Define the fitness function: for each individual in the population, calculate its fitness value Fitness under the fitness function f i :

[0081]

[0082] where y j is the chemical true value of the β-HMX content of the jth sample in the calibration set; y j ' is the predicted value of the β-HMX content of the jth sample in the calibration set, and n1 is the number of calibration set samples.

[0083] c) Coefficient dynamic updating mechanism: the improvement of the adaptive grey wolf optimization algorithm compared with the traditional grey wolf optimization algorithm lies in the dynamic updating of the aggregation coefficient a and the step size g.

[0084] The purpose of dynamic updating of the aggregation coefficient a is to adaptively adjust the search behavior according to the progress of the algorithm, so as to better balance the exploration and development stages. The dynamic updating of the step size g is to control the search granularity by adjusting the step size, and then realize different search behaviors in different stages of the algorithm.

[0085] The updating formulas of the aggregation coefficient a and the step size g are as follows:

[0086]

[0087]

[0088] In the formula, a new is the new aggregation coefficient used in the current iteration; b is the attenuation coefficient; t is the current iteration number; gnew new step size used in the current iteration; mean Fitness mean fitness of the population; std Fitness standard deviation of the population fitness.

[0089] d) Position Update: The new position of each individual S best ' is calculated by the following position update equation: i A, C, D and the position of each individual S new new

[0090] A = 2 · a best · rand(dim) - a i

[0091] C = 2 · rand(dim)

[0092] D = | C · S best ' - S best |

[0093]

[0094] where N(0, 1) denotes the standard normal distribution;

[0095] e) Boundary Check and Fitness Update: Boundary check is performed for the new position, if it is out of the search space, then it is corrected. Then, the fitness value is recalculated based on the new position.

[0096] f) Repeat the above steps until the maximum number of iterations T is reached;

[0097] g) Output the global optimal fitness value and its optimal position S best , which is the extracted key spectral variable.

[0098] In a specific use, the dimension dim is 886, the population number N is 30, the maximum iteration number T is 100, the upper bound ub and the lower bound lb are 1 and 0 respectively, the aggregation coefficient a is 1.202, the attenuation coefficient b is 1.147, the step length g is 1.298, and the extracted key spectral variables are 16, 19, 20, 25, 27, 36, 50, 63, 69, 79, 94, 101, 103, 132, 135, 140, 147, 159, 167, 171, 174, 211, 215, 216, 221, 226, 228, 243, 246, 248, 275, 283, 297, 298, 301, 305, 306, 318, 326, 333, 355, 356, 358, 360, 368, 379, 387, 394, 396, 398, 400, 403, 409, 410, 411, 412, 432, 460, 461, 476, 485, 497, 498, 502, 507, 526, 530, 532, 534, 539, 552, 558, 564, 566, 569, 584, 586, 589, 595, 597, 605, 610, 612, 619, 622, 625, 642, 647, 687, 708, 716, 720, 733, 737, 752, 753, 757, 762, 765, 791, 799, 815, 824, 835, 873, 878, 879, 880, 883, 885.

[0099] S5, using the optimized corrected sample spectrum matrix subset and the corresponding concentration matrix, a PLS quantitative model is established by using a five-fold cross-validation method of partial least squares (PLS). The variable-optimized sample spectrum matrix subset of the validation set is predicted, and the prediction performance of the PLS quantitative model on the validation set samples is evaluated according to the determination coefficient R 2 and the mean absolute error MAE;

[0100] In step S5, the determination coefficient R 2 and the mean absolute error MAE are:

[0101] ① The determination coefficient R 2 : R 2 closer to 1 indicates that the prediction effect of the model is better, and the calculation formula is:

[0102]

[0103] In the formula, y is the predicted value of the lth sample in the validation set; y l is the chemical true value of the lth sample in the validation set; and y is the average value of the chemical true values of the validation set samples.

[0104] ②Mean Absolute Error (MAE): The smaller the MAE value, the better the prediction accuracy of the model. The calculation formula is:

[0105]

[0106] In the formula, n2 is the number of samples in the verification set.

[0107] Table 1 shows the results of variable optimization of Raman spectrum data by different variable optimization algorithms (full spectrum, traditional grey wolf optimization algorithm and adaptive grey wolf optimization algorithm of the present application), and the results after variable optimization are used to establish PLS quantitative model by five-fold cross-validation method of partial least squares (PLS).

[0108] As shown in Table 1, the PLS quantitative model established by full spectrum has determination coefficients R 2 and mean absolute error MAE of 0.99879, 0.84904, 0.99673 and 1.41297 in the calibration set and the verification set, respectively. After variable optimization by traditional grey wolf optimization algorithm and adaptive grey wolf optimization algorithm, the number of variables is greatly reduced, and the determination coefficients R 2 and mean absolute error MAE in the calibration set and the verification set are improved, indicating that after variable optimization, the information of redundant variables is eliminated, and the adaptability of the model to unknown data is enhanced, i.e. the generalization ability is improved. Compared with the traditional grey wolf optimization algorithm, the adaptive grey wolf optimization algorithm improves the global search ability and local search efficiency due to the dynamic updating of the aggregation coefficient α and the step size γ, thereby improving the search performance in the variable space and more conducive to reducing the number of variables while improving the overall prediction performance of the model. After Savitzky-Golay first derivative preprocessing and variable optimization by adaptive grey wolf optimization algorithm, the determination coefficients R 2 and mean absolute error of the model in the calibration set are 0.99940 and 0.58389, respectively, and the determination coefficients R 2 and mean absolute error of the model in the verification set are 0.99880 and 0.83965, respectively, indicating the feasibility of the method of the present application, which can be used for real-time monitoring of the HMX crystallization process.

[0109] Table 1 Comparison results of different variable optimization algorithms

[0110]

[0111] S6, at room temperature, 48 g of crude a-HMX was completely dissolved in 240 mL of mixed solvent (the mixed solvent was composed of acetone, ethyl acetate and water, the volume ratio was 5:4:1), added to the stirred tank reactor with stirring rate of 115 rpm / min, heated to 55°C and kept constant; 2.88 g of β-HMX was taken as the induction seed and added to the stirred tank reactor to form a mixed solution, and the crystal transformation was carried out; the spectra were collected by blaze online Raman spectrometer in a certain time (sampling 25 times), and after the Savitzky-Golay first-order derivative pretreatment and the screening of the key spectral variables obtained by the adaptive grey wolf optimization algorithm in step S4, the prediction set sample spectrum matrix subset was obtained, and the β-HMX content at this time was calculated by the PLS quantitative model established in step S5, and compared with the mixed solution taken synchronously, which was filtered, dried and treated, and then the β-HMX content was measured by X-ray diffractometer (see Figure 2 ). It was found that the results obtained by real-time monitoring by the established PLS quantitative model were basically consistent with the results measured by X-ray diffractometer, indicating that the method constructed could be used for real-time monitoring of the HMX crystal transformation process.

[0112] Therefore, the HMX crystal transformation real-time monitoring method of the present application combined with online Raman has important application value for real-time monitoring of the crystal transformation process.

[0113] Finally, it should be noted that: the above examples are only used to illustrate the technical solutions of the present application but not to limit it, although the present application has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that: the technical solutions of the present application can still be modified or replaced by equivalents, and these modifications or equivalent replacements cannot make the modified technical solutions deviate from the spirit and scope of the technical solutions of the present application.

Claims

1. A method for real-time monitoring of HMX crystallization by machine learning combined with online Raman, characterized in that, Comprise the following steps: S1, take α-HMX solid and completely dissolved in mixed solvent, add to the stirred tank reactor after starting stirring, heating to constant temperature; Take β-HMX as inducing crystal seed and add it to the stirred tank reactor to form a mixed solution, and carry out heat preservation and crystal transformation; During the crystal transformation process, collect spectra every 3 minutes using an online Raman spectrometer, and at the same time, take 1% of the mixed solution, filter, dry, and then use an X-ray diffractometer to determine the content of β-HMX therein; S2, the collected Raman spectrum spectrum matrix and the corresponding β-HMX content chemical value matrix determined by the X-ray diffractometer constitute a data set, and the calibration set sample and the validation set sample are obtained according to the random division method; S3, the calibration set sample spectrum matrix and the validation set sample spectrum matrix are pretreated by Savitzky-Golay first-order derivative; S4, the spectrum matrix of the calibration set sample is selected by using the adaptive grey wolf optimization algorithm, and the extracted key spectral variables are used to screen data from the spectrum matrix of the calibration set sample and the spectrum matrix of the validation set sample, to obtain the optimized calibration set sample spectrum matrix subset and the validation set sample spectrum matrix subset; S5, using the optimized calibration set sample spectrum matrix subset and the corresponding concentration matrix, a PLS quantitative model is established by using the five-fold cross-validation method of partial least squares; and the prediction performance of the PLS quantitative model is verified by using the validation set sample spectrum matrix subset and the corresponding concentration matrix; S6, the newly collected Raman spectrum data are screened by the Savitzky-Golay first-order derivative of step S3 and the key spectral variables obtained by the adaptive grey wolf optimization algorithm in step S4 to obtain a prediction set sample spectrum matrix subset, and the PLS quantitative model established in step S5 is used to monitor the HMX crystal transformation process in real time. 2.The HMX crystallization real-time monitoring method of machine learning combined with online Raman according to claim 1, wherein, In step S1, the ratio of α-HMX solid to mixed solvent is 1g:5-20ml. 3.The HMX crystallization real-time monitoring method of machine learning combined with online Raman according to claim 1, wherein, In step S1, the mixed solvent is composed of acetone, ethyl acetate and water, and the ratio of acetone, ethyl acetate and water is 10-80:10-80:10-80, and the total of the three is 100%. 4.The HMX crystallization real-time monitoring method of machine learning combined with online Raman according to claim 1, wherein, In step S1, the stirring rate is 100-300 rpm / min; the constant temperature is 30-70℃; and the mass of β-HMX added is 5-20% of the mass of α-HMX solid. 5.The HMX crystallization real-time monitoring method of machine learning combined with online Raman according to claim 1, wherein, In step S3, the specific steps of Savitzky-Golay first-order derivative are as follows: S31, select appropriate window size w, polynomial degree p and difference item d, wherein d=1; S32, construct a design matrix V, where each row corresponds to a given sample x k (x1,x2,x3……x n ) of Raman spectral data, each column corresponds to a different power term and derivative term: S33, calculate the filter coefficient E by least squares method: E = (V T V) -1 V T Y Wherein, Y is the matrix to be solved: S34, calculating a filtered given sample x using the filter coefficient E k the output matrix Y' : Y'=VE S35, repeat steps S31 to S34, that is, obtain the calibration set sample spectrum matrix X1 and the validation set sample spectrum matrix X2 after Savitzky-Golay first-order derivative processing.

6. The machine learning combined with online Raman HMX crystallization real-time monitoring method according to claim 5, characterized in that, In step S4, the specific steps of the adaptive grey wolf optimization algorithm are as follows: S41, initialize parameters: determine the dimension dim of the optimization problem, the number of grey wolf populations N in the algorithm, the maximum number of iterations T, and set the upper bound ub and the lower bound lb for each dimension; Position S of each individual in the population i Randomly initialize within the search space, i = 1, 2, 3, …, N, i.e.: S i = lb + (ub - lb) * rand(dim) S42, define fitness function: for each individual of the population, calculate its fitness value under the fitness function f i : wherein y j is the chemical true value of the β-HMX content of the jth sample in the calibration set sample; y j is the predicted value of the β-HMX content of the jth sample in the calibration set sample, and n1 is the number of calibration set samples. S43, coefficient dynamic updating mechanism: the updating formula of the aggregation coefficient a and the step size g is respectively: where a new is the new aggregation coefficient used in the current iteration; b is the decay coefficient; t is the current iteration number; g new is the new step size used in the current iteration; mean Fitness is the mean fitness of the population; std Fitness is the standard deviation of the population fitness; S44, position update: by the following position update formula, combined with the current optimal individual position S best , through the variable A, the variable C, the variable D, and the position S of each individual i Calculate the new position of each individual A = 2 · a new • rand(dim) - a new C = 2 · rand(dim) D = | C • S best - S i | Wherein, N(0, 1) indicates a standard normal distribution; S45, boundary check and fitness updating: boundary check is performed on the new position, if it is out of the search space, it is modified; then, the fitness value is recalculated based on the new position; S46, repeat steps S41-S45 until the maximum iteration number T is reached; S47, output the global optimal fitness value and its optimal position S best where S best That is, the extracted key spectral variables.

Images