A grinding process crushing efficiency prediction method based on PINN and drift correction

Abstract

The application discloses a kind of based on PINN and drift correction's grinding process crushing efficiency prediction method, belong to mineral processing equipment technical field.The method includes: obtaining DEM-BPM simulation, industrial sensing and ore dressing granularity test etc. multimodal data, unified representation is realized by space-time alignment and feature level fusion;Physical information neural network (PINN) is constructed embedding data fitting, mass conservation and Bond crushing energy balance constraint, and based on residual gradient norm dynamic distribution weight;Quantitative prediction uncertainty is carried out using deep integration method, and the crushing efficiency prediction value is output;Design parameter drift correction mechanism based on cumulative sum control chart test and Kalman filtering parameter tracking, realize the model parameter update when the change of ore property.The application solves the problems that the existing prediction method lacks physical law constraint, leading to poor output consistency and model parameters cannot adapt to the long-term prediction accuracy attenuation caused by the time-varying drift of working condition, and can realize long-term, stable prediction of grinding process crushing efficiency.

Description

Technical Field

[0001] This invention relates to the field of mineral processing equipment technology, and specifically to a method for predicting crushing efficiency in grinding processes based on PINN and drift correction. Background Technology

[0002] Ball mills, as core equipment in the mineral processing industry, primarily undertake the tasks of ore crushing and grinding. Their working principle involves the rotation of the cylinder driving the movement of the internal ore and steel balls, utilizing the collision, impact, and friction between particles to achieve ore crushing and grinding. However, the complexity of ore properties, the uneven particle size distribution, and the dynamic fluctuations of operating parameters can lead to unstable crushing efficiency in ball mills. This is often accompanied by problems such as over-grinding, under-grinding, excessive energy consumption, and high steel consumption, severely impacting production efficiency and economic benefits.

[0003] Currently, the main shortcomings of ball mill crushing efficiency prediction methods are as follows: First, there are data-driven machine learning methods. These methods utilize operational data collected by field sensors to establish a mapping relationship between operating parameters and crushing efficiency through regression models. However, these methods have the following problems: (1) they ignore the physical constraints of the crushing process, which may lead to biased predictions that violate basic physical laws such as mass conservation and energy balance; (2) they only output deterministic single-point values ​​and cannot provide a reliable quantitative evaluation of the prediction results. Second, there are mechanism prediction methods based on pure discrete element simulation. These methods simulate the entire ore crushing process by establishing a DEM-BPM numerical model, which has clear physical interpretability, but suffers from bottlenecks such as large computational load and long simulation time, and cannot meet the real-time requirements of industrial sites. Thirdly, there is the hybrid method of "simulation + artificial intelligence algorithm". Its essence is still to use artificial intelligence algorithm as an accelerator. It does not solve the following defects: (1) Physical constraints are embedded in a fixed hyperparameter manner. During the training process, the gradient magnitude of each loss term is significantly different, which leads to the physical constraint terms being "submerged" or "over-constrained" too early, affecting the model convergence quality; (2) In terms of uncertainty quantification, existing methods mostly use MC Dropout to approximate Bayesian inference, but MC Dropout cannot separate cognitive uncertainty from accidental uncertainty, and has limited support for diagnosing the cause of prediction failure; (3) In the face of time-varying drift of ore properties, existing methods can only detect the occurrence of drift, but cannot quickly and dynamically track and correct key model parameters under sparse sample conditions. The parameter drift response has obvious lag.

[0004] In summary, developing a ball mill crushing efficiency prediction method that can adjust physical constraint weights, provide cognitive / accidental uncertainty decomposition output, and has the ability to quickly track and correct parameters using Kalman filtering is a technical challenge that the industry urgently needs to overcome. Summary of the Invention

[0005] To overcome the above-mentioned technical problems, the present invention provides a method for predicting crushing efficiency in grinding processes based on PINN and drift correction.

[0006] To achieve the above technology, the following steps are included: S1. Construct a multimodal data system that includes discrete element simulation of the ore crushing process, real-time acquisition of the grinding process by sensor array, and manual ore discharge particle size analysis. S2. Perform unified alignment and fusion operations on the multimodal data system obtained in S1 in time and feature space to generate fused feature vectors; S3. Input the fused feature vector obtained in S2 into the physical information neural network PINN with a physical constraint weight adjustment mechanism to predict the breaking efficiency, and obtain the trained network. S4. Construct a preset number of network structures that are identical to the network in S3 but satisfy different random initializations. Estimate the predicted distribution of the crushing efficiency using the mean and variance, and output the expected value and probability prediction interval. S5. Continuously monitor the prediction error sequence of S4 based on the CUSUM control chart method; S6. When the CUSUM statistic exceeds the control threshold, it is determined that the ore properties have drifted, triggering the two-layer drift correction process. S7. After training is completed, read the real-time multi-source data of the grinding process and output the final expected value of crushing efficiency and its probability prediction range. Based on S1 to S7, a method for predicting crushing efficiency in the grinding process was constructed.

[0007] Furthermore, in step S1, the data acquisition step for discrete element simulation of the ore crushing process includes: First, a three-dimensional geometric model including the cylinder and liners uniformly distributed along the internal axis of the cylinder is created using SolidWorks. The model is exported as a .STL file and imported into EDEM. The BPM model is used to construct aggregates formed by bonding bonds at a preset size to represent the ore. Then, material properties are set for the cylinder, liners, ore, and steel balls generated by the simulation software. Secondly, the normal and tangential contact forces between particles were calculated using the Hertz-Mindlin model for the established geometric model; the rolling friction torque was calculated using the Standard Rolling Friction model; the bond breakage mechanism in the ore was simulated using the particle bonding model BPM; and the physical rules for a single simulation were defined by randomly stacking ore and steel balls in the model through initialization settings. Next, based on the physical rules of a single simulation, multiple sets of working conditions were obtained through orthogonal design to provide experimental schemes for the simulation process. Among them, the orthogonal experiment was set as a three-factor, three-level orthogonal experiment. The three factors included: cylinder rotation speed, total filling rate, and steel ball filling rate. The three levels included: cylinder rotation speed set to 14, 17, and 20 r / min; total filling rate set to 25%, 35%, and 45%; and steel ball filling rate set to 8%, 11.5%, and 15%. Finally, by performing DEM-BPM simulation, the number of bond fractures was recorded, the number of bond fractures was converted into fracture efficiency, and the simulation dataset was stored.

[0008] Furthermore, in S1, the data acquisition step of the sensor array in real time collecting data of the grinding process includes: using preset type sensors to collect five types of high-frequency signals: feed rate, water supply rate, return sand rate, discharge slurry flow rate, and discharge slurry concentration.

[0009] Furthermore, in step S1, the data acquisition steps for the manual ore discharge particle size analysis process include: Samples of discharged ore are collected at a fixed sampling period. The Bond work index is calculated based on the sieve aperture diameter of particles passing through a preset percentage in the discharged ore and the sieve aperture diameter of particles passing through a preset percentage in the raw ore. The calculated Bond work index is used as the numerator, and the unit mass power consumption measured on-site is used as the denominator to calculate the true crushing efficiency. All test data are stored as a calibration dataset according to the timestamp.

[0010] Furthermore, step S2 includes: First, an alignment window is defined with the test time as the main anchor point. Sensor statistical features are aggregated within this window. The most recent simulation conditions are found and the simulation conditions closest to the test time are extracted. The spliced ​​features are obtained by splicing the data. Secondly, the splicing features are standardized using Z-Score. Finally, a multi-head attention mechanism is used to output the fused feature vector after projection of the spliced ​​features.

[0011] Furthermore, step S3 includes: First, a PINN network derived from a fully connected feedforward neural network is constructed, with the following settings: the input layer dimension is the same as the fused feature vector dimension; 6 feature extraction layers, each containing a linear transformation, BatchNorm, and ReLU activation function; a physical constraint embedding branch is set between the 3rd and 4th layers to calculate the mass conservation constraint residual and the energy balance constraint residual; the output layer outputs the predicted value. Secondly, the computational network includes: a total loss function comprising the data fitting loss function, the mass conservation constraint loss function, and the Bond breakage energy balance constraint loss function, expressed as follows: In the formula, L data This represents the loss function for data fitting. L mass Represents the loss function due to the mass conservation constraint; L energy This represents the loss function for the energy balance constraint of Bond breakage. α (t) and β (t) represent the first and second dynamic weights, respectively; Next, using the GradNorm algorithm, by calculating the gradient norm and relative training rate for each loss task, the optimization objective is set to make the gradient norm of each loss task equal to the average gradient norm multiplied by the relative training rate, as shown in the following expression: In the formula, Represents the gradient norm of the quality task; This represents the mean of the gradient norms for each task; Indicates the relative training rate of the quality task; Represents the energy task gradient norm; This represents the relative training rate of the energy task; γ is an asymmetric hyperparameter with a value of 1.5. Then, after each gradient step, the constraint weights are updated according to the following rules: In the formula, or w The learning rate is updated for the weights, with a value of 0.025. C w This is the normalization constant for the sum of the constraint weights, and its value is 1. Finally, the Adam optimizer is used for training.

[0012] Furthermore, step S4 includes: Five PINN networks with the same structure but different initialization parameters were constructed, and the training data were randomly shuffled to ensure diversity. For each PINN network, the predicted mean and predicted variance are simultaneously output, and it is trained with negative log-likelihood loss. The ensemble prediction distribution is calculated by fusing the prediction mean and prediction variance of each PINN network using an equal-weight Gaussian mixture approximation. This distribution includes the ensemble prediction mean and ensemble prediction variance. The ensemble prediction variance comprises cognitive uncertainty and random uncertainty, expressed as: In the formula, This indicates cognitive uncertainty. Indicates randomness and uncertainty; The integrated prediction mean is used as the final predicted value for the crushing efficiency point; under the assumption of normal distribution, the mean is used as the final predicted value for the crushing efficiency point. As a probability prediction interval with a 95% confidence level, where, Indicates the standard deviation of the integrated forecast; or ens Indicates the ensemble forecast mean; The indicator for calculating the proportion of cognitive uncertainty is: in, Indicates the variance of the integrated predictions; when R epi When the accuracy is >60%, it is determined that the current prediction is mainly driven by insufficient model knowledge, and simulation or test data for the corresponding working condition area should be supplemented. when s ens / or ens When the accuracy rate is greater than 15%, a warning about the reliability of the prediction will be issued.

[0013] Furthermore, step S5 includes: Calculate the predicted relative error at each test calibration time. : The error sequence is standardized using the mean of the historical steady-state window. m 0 and standard deviation s 0 as the baseline: In the formula, Indicates the first The relative error of the prediction time for the next test; Indicates the first Model predictions at the time of the next test; This represents the standardization error.

[0014] Calculate the two-sided CUSUM statistic to detect the positive and negative drift of the error mean, as shown in the following expression: In the formula, Indicates the first Positive CUSUM statistic at the time of the test; Indicates the first The negative CUSUM statistic at the time of the test; Indicates the first Positive CUSUM statistic at the time of the test; Indicates the first The negative CUSUM statistic at the time of the test; As a reference value, 0.5 is used to detect mean shifts greater than 1σ; initial value of CUSUM. = =0.

[0015] When the CUSUM statistic exceeds the control limit At that time, it was determined that the properties of the ore had undergone a significant shift: like > or > If this occurs, a drift warning will be triggered.

[0016] To control the limit, we take 4, which corresponds to the average running length under normal operating conditions without ore property drift. After the drift alarm is triggered, the drift start time is recorded, the CUSUM statistic is reset, and the two-layer drift correction process in step S6 is started.

[0017] Furthermore, the two-layer drift correction process of S6 is as follows: First layer: Use Kalman filtering to perform state estimation and dynamic tracking update of the bias parameters of the PINN output layer; The second layer involves constructing an incremental training dataset and updating all parameters of PINN in a limited number of rounds using a forgetting factor-weighted loss function.

[0018] Furthermore, step S6 includes: The bias parameter vector of the PINN output layer is modeled as the hidden state of a random walk process, and Kalman filtering is used for online dynamic estimation of the optimal value, including: In the prediction phase, the Kalman state transition equation is expressed as: In the formula, This represents the prior bias estimate; Indicates the first The posterior bias estimate after correction; Represents the prior error covariance matrix; Indicates the first The posterior error covariance matrix; Represents the process noise covariance matrix; During the correction phase, the Kalman update equations include: Calculate the observation residuals based on the newly collected laboratory calibration data: In the formula, Indicates the first The observation residuals of the step; Indicates the first The true crushing efficiency of the test; No. Predicted crushing efficiency from the first test; Calculate Kalman gain : In the formula, Indicates the first Kalman gain of the step; This represents the transpose of the observation matrix; The observation matrix; This represents the observation noise covariance, taken as 0.01; Posterior estimation of bias parameters and update of error covariance: In the formula, Indicates the first Posterior bias estimation of the step; Represents the identity matrix; Update the bias parameters By directly replacing the PINN output layer bias, fast parameter correction can be achieved without retraining. Forgetting factor weighted incremental fine-tuning: After a drift alarm is triggered, newly added drift region samples are continuously accumulated. When the number of new samples reaches the batch threshold of 10, a full parameter fine-tuning is performed, as follows: In the formula, This represents the total loss function for incremental fine-tuning; This represents the forgetting factor, with an initial value of 0.8, which decays to 0.9 times the previous value after each trigger. Represents the loss on the core historical dataset, where, This represents the core historical dataset; This represents the loss in the newly added drift region data, where, This indicates a newly added drift region dataset; This represents the EWC regularization coefficient, set to 0.01; This represents all parameters of the model; This represents the model parameters before fine-tuning; This indicates that the flexible weights are fixed regularization terms to prevent catastrophic forgetting; Using the current model parameters as initial values, set the learning rate to 10.-4 Perform training updates in a limited number of rounds; The dual-layer calibration mechanism is set up as follows: the first layer takes effect in the first test cycle after drift; the second layer provides depth adaptation to the new working condition distribution after accumulating 10 samples. If the CUSUM statistic continues to exceed the control limit after three consecutive drift triggers, a major change alarm for ore properties will be sent and the full DEM simulation case extension will be re-executed to supplement the simulation knowledge base.

[0019] Compared with the prior art, the beneficial effects of the present invention are reflected in the following aspects: (1) The physical constraint weight adjustment mechanism based on the GradNorm algorithm solves the problem of poor convergence quality caused by gradient magnitude imbalance in fixed-weight PINN. By dynamically balancing the gradient contribution of each loss term, the model can effectively learn the data fitting rules and physical consistency constraints at the same time, thereby improving the model's prediction accuracy and physical rationality.

[0020] (2) The model prediction uncertainty was quantified in detail using deep learning methods, and the overall uncertainty was decomposed into cognitive uncertainty and accidental uncertainty, providing operators with richer diagnostic information than the traditional Monte Carlo dropout algorithm.

[0021] (3) The CUSUM control chart is used to continuously monitor the operating condition drift of the forecast results. Compared with the traditional sliding window mean test method, the drift trend can be detected 2 to 5 calibration cycles in advance. At the same time, the Kalman filter parameter correction mechanism can complete the optimal parameter estimate update with only a single new data point, which effectively shortens the parameter correction cycle. Attached Figure Description

[0022] Figure 1 This is a schematic diagram of the steps of the present invention; Figure 2 This is a schematic diagram of the discrete element simulation model of the ball mill constructed in this invention; Figure 3 The BPM particle bonding model constructed for this invention, wherein, Figure 3 (a) is the overall structure diagram; Figure 3 (b) is a local particle bonding model; Figure 4 This is a visualization example of the expected value curve of crushing efficiency prediction and the prediction result of probability prediction interval in an embodiment of the present invention. Figure 5 This is a visualization example of the time series prediction results of cognitive / accidental uncertainty decomposition in an embodiment of the present invention. Detailed Implementation

[0023] The present invention will be further described in detail below with reference to specific embodiments.

[0024] See Figure 1 A method for predicting crushing efficiency in grinding processes based on PINN and drift correction includes the following steps: S1. Construct a multimodal data acquisition system, including three data sources. The construction steps are as follows: S1.1, Simulation data acquisition; A discrete element simulation model (DEM-BPM) for the grinding process was established. The Hertz-Mindlin contact model and the BPM particle bonding model were used to simulate the ore crushing process. The number of bond breakages was used as the simulation characterization index of crushing efficiency. Orthogonal working condition batch offline simulations were performed to obtain a low-frequency, high-precision simulation dataset. The specific steps included: S1.1.1 Establish the three-dimensional geometric model and BPM model of the ball mill, and set the actual parameters of the ball mill; The three-dimensional geometric model of the ball mill includes: a cylinder and liners evenly distributed along the internal axis of the cylinder; the liners are further divided into feed end liners, cylinder liners, and discharge end liners, which are used to bear particle collisions and record changes in bonding bonds; The BPM model parameters include: number of subparticles, size, bond strength, bond stiffness, and bond radius; The actual parameters of the ball mill include: the material properties of the cylinder and the liner, among which: The material properties of the cylinder and liner include: shear modulus, Young's modulus, density, and Poisson's ratio; The geometry of a steel ball is spherical; the material properties of a steel ball include: shear modulus, Young's modulus, density, and Poisson's ratio. The geometry of the ore is determined by the geometry of the BPM model; the material properties of the ore include: shear modulus, Young's modulus, density, and Poisson's ratio; A ball mill model was built using SolidWorks software. After mesh reconstruction, the model was exported as an .STL file and then imported into EDEM software to complete the construction of the discrete element simulation model. like Figure 2 As shown, a discrete element simulation model of a ball mill is established, and the material property parameters of each component are set as follows, where the mass of the steel ball is 2kg and the radius is 40mm. The specific settings are shown in Table 1. Table 1: Material parameter settings for the discrete element simulation model of the ball mill Furthermore, such as Figure 3 As shown, the geometry of the ore is determined by the BPM model, with a radius of 20 mm, and consists of 127 sub-particles with a diameter of 6 mm connected by bonding bonds to form an agglomerate.

[0025] S1.1.2. An ore crushing model is constructed using the Hertz-Mindlin (no slip) model, Standard Rolling Friction model, and Bonded-particle model (BPM) in the discrete element method. This model is used for numerical simulation to calculate the number of bond breaks. The steps are as follows: S1.1.2.1 Constructing the Hertz-Mindlin model to calculate the normal contact force of the particles. and tangential contact force The expression is as follows: In the formula, The normal elastic modulus of the particle; The tangential elastic modulus of the particle; This represents the normal overlap between particles; This represents the tangential overlap between particles; , , , The calculation method is as follows: In the formula, Represents the equivalent elastic modulus; Indicates the equivalent radius of the particle; Indicates the equivalent shear modulus; Indicates particles The equivalent elastic modulus; Indicates particles The equivalent elastic modulus; Indicates particles Poisson's ratio; Indicates particles Poisson's ratio; Indicates particles The equivalent contact radius; Indicates particles The equivalent contact radius; Indicates particles The equivalent shear modulus; Indicates particles The equivalent shear modulus; S1.1.2.2 Constructing the Standard Rolling Friction model to calculate rolling friction. The expression is as follows: In the formula, The coefficient of rolling friction between particles; This is the distance from the particle contact point to its center of mass. The vector of the unit angular velocity of the particle at the point of contact; S1.1.2.3 Constructing the Bonded-Particle Model (BPM): The bond parameters that need to be determined in the simulation are as follows: normal stiffness, shear stiffness, tensile strength, shear strength, and bond radius, with the following expressions: In the formula: , These represent the increments of axial force and shear force experienced by the material per unit time, respectively. , These are the axial displacement increment and shear displacement increment per unit time, respectively; A The cross-sectional area of ​​the bonding bond; When the maximum tensile stress exceeds the tensile strength ( ) or the maximum shear stress exceeds the shear strength ( If the bond breaks, the adhesive bond will break. Where: In the formula: , These represent the normal force and shear force acting on the particle, respectively. , These are the normal torque and shear torque acting on the particle, respectively. I The moment of inertia of the bond; R The radius of the bond is 4 mm.

[0026] S1.1.2.4 Simulation initialization: Set up several ores and steel balls, which are randomly and naturally piled in the model cylinder. The calculation formula between the number of ores and steel balls and the filling rate is as follows: In the formula, This indicates that the total number of particles is 127, with a diameter of 6 mm; Fill rate; The effective radius of the cylinder; The length of the cylinder; This indicates the particle size, specifically the particle radius.

[0027] S1.1.3. An orthogonal experimental design was adopted, with the mill speed as the starting point. n Range 14~20 r / min, total fill rate f total Range 25%~45%, steel ball filling rate f ball The range of 8% to 15% is a variable, forming a group that includes... N group = Table 2 of orthogonal working conditions for 27 working conditions; Table 2: Orthogonal Operating Conditions S1.1.4 Perform DEM-BPM simulation for each working condition and record the number of bond breaks. N broken Its expression is as follows: In the formula, N broken is the total number of bond breaks per unit time; T is the total number of simulation time steps, which is set to 10s in this embodiment; M is the total number of broken model agglomerates in the model. For the t-th time step The number of bond breaks in an aggregate.

[0028] S1.1.5, Further calculate the simulated crushing efficiency: In the formula, The initial total number of bond bonds is expressed as follows: In the formula, C These are the fitting coefficients. This represents the total number of sub-particles that make up the BPM model.

[0029] S1.1.6. Store the data as a simulation knowledge dataset. ,in, Used to represent simulation knowledge datasets, Indicates the first The simulation input feature vector of the group working conditions, Indicates the first Simulated crushing efficiency under various working conditions; S1.2 Industrial field sensor data acquisition; Through an industrial field sensor array, high-frequency parameter data during the grinding process are collected in real time, including: feed rate, water rate, return sand rate, discharge slurry flow rate, and discharge slurry concentration, specifically including: feed rate. F feed ( t ), water supply F water ( t ), amount of returned sand F return ( t ), discharge slurry flow rate Q dis ( t ) and ore discharge slurry concentration Q sl ( t ).

[0030] After preprocessing (noise reduction, filtering), the raw signals from each sensor are used to extract time-domain statistical features (mean, root mean square, peak factor, kurtosis) to form a high-frequency sensor feature vector. X sens ( t In this invention, the sampling frequencies of each sensor are set as follows: 2000Hz for the accelerometer, 100Hz for the torque sensor, 500kHz for the acoustic emission sensor, and 1Hz each for the flow meter / concentration meter / belt scale.

[0031] S1.3, Collection of ore particle size analysis data; The particle size analysis data of the discharged ore is collected manually at regular intervals and converted into a true crushing efficiency calibration value using the Bond work index method. The specific steps include: S1.3.1, Sampling according to a fixed sampling period T sample =8h to collect ore discharge samples, sieve analysis and extraction P 80 and F 80 The parameters, using the Bond work index method, are used to calculate the theoretical specific crushing energy: In the formula, W i The Bond work index of the ore; P 80 The sieve aperture diameter at which 80% of the particles in the ore discharge pass ( μm ); F 80 The sieve aperture diameter through which 80% of the particles in the raw ore pass ( μm ); S1.3.2 Calculate the actual crushing efficiency: In the formula,W actual The power consumption per unit mass was measured on-site; all test data were stored as a calibration dataset according to timestamps. ,in, Used to represent laboratory calibration datasets, Indicates the first The time of the next test Indicates the first Operating parameters at the time of the next test; Indicates the first The actual crushing efficiency corresponding to the test result; S2, Spatiotemporal alignment and feature-level fusion; The design of the spatiotemporal alignment and feature-level fusion module unifies and fuses the three types of multimodal data (simulation data, sensor data, and laboratory data) in terms of time scale and feature space to generate a fused feature vector, including the following steps: S2.1 Time Scale Alignment: Based on the test time Use the main anchor point to define the alignment window Δ T align =4h, aggregated sensor statistical characteristics X sens,agg ( Find the most recent simulation conditions. extract Concatenate and merge candidate vectors: In the formula, Indicates the first The fusion candidate vector at each test time; Indicates the first Statistical characteristics of the sensors aggregated at each test time; Indicates the relationship with the first The simulation conditions most closely approximating the time of the next test; Indicates the first On-site operating parameters at the time of the next test.

[0032] S2.2 Feature Space Normalization: Z-Score normalization is performed on each feature dimension of the fused candidate vector. In the formula, X norm The features are normalized. X The original input features, m This is the mean of the feature. s This is the standard deviation of the feature.

[0033] S2.3. Employing Multi-Head Attention (MHA) to achieve deep fusion of heterogeneous features: In the formula, Q For query matrix; K The key matrix; V It is a value matrix; For the first The output of each attention head, ; W O To output the projection matrix; d k The dimension of the key vector is 64; The final output of the bulls' attention. H =8; , and They represent the first The projection matrix of the query, key, and value of each attention head; The final output is a fused feature vector. ,in, d fused This indicates the dimension of the fused feature vector, which is 512.

[0034] S3. Construct a PINN with embedded physical constraints; A Physics-Informed Neural Network (PINN) with an embedded physical constraint weight adjustment mechanism is constructed. The network input is a fused feature vector, and the output is a predicted value of the crushing efficiency. The loss function of the PINN consists of three weighted parts: a data fitting loss function, and a loss function. L data Mass conservation physical constraint loss function L mass and Bond breakage energy balance constraint loss function L energy The weights of each constraint term are dynamically assigned using the GradNorm algorithm based on the residual gradient norm. The specific steps include: S3.1 The physical information neural network is a fully connected feedforward neural network with the following layer configuration: Input layer: Receives the fused feature vector Z fused Dimension d fused =512; Feature extraction layer: composed ofN fe = It consists of 6 fully connected layers, each with 256 hidden units; each layer contains a linear transformation, BatchNorm and ReLU activation function; Physical constraint embedding branch: This branch is introduced between layers 3 and 4 of the feature extraction layer to calculate the mass conservation constraint residuals. r mass and energy balance constraint residuals r energy ; Output layer: Fully connected layer outputs predicted crushing efficiency. .

[0035] S3.2 The total loss function of the physical information neural network is: In the formula, α (t) and β (t) represents the first and second dynamic weights; data fitting loss function L data The mean square error form is used: In the formula, Represents the total number of training samples, where An index belonging to the total number of training samples; Indicates the first Predicted values ​​for each sample; Indicates the first The true value of each sample; Mass conservation constraint loss function L mass According to the material balance in the ball milling process, the total amount of material entering the mill equals the total amount of discharged ore plus the change in the amount of material retained in the mill. In the formula, Indicates the first The total amount of grinding material in each sample; express Ore feed rate at the time; express The amount of sand returned at that time; This indicates the average discharge slurry flow rate within the aligned window; This indicates the average discharge slurry concentration within the alignment window; Indicates the density of the discharged slurry; This indicates the change in the amount of material remaining inside the mill; This represents the relative residual due to the conservation of mass. Bond breakage energy balance constraint loss function L energy for: In the formula, This represents the actual fracture work quantity predicted by the model; Indicates the first Measured power consumption per sample; Indicates the amount of ore processed; This represents the relative residual of energy balance; This indicates a loss in energy balance.

[0036] S3.3 The adaptive physical constraint weight adjustment mechanism of the physical information neural network based on the GradNorm algorithm of residual gradient norm is as follows: To avoid the problem of unbalanced magnitudes among the loss terms, the GradNorm algorithm is introduced to dynamically adjust the weights of the physical constraints. α (t) and β (t); The GradNorm algorithm automatically adjusts the weights of each task to make the gradient norm of each task tend to be balanced by monitoring the gradient norm of each loss term with respect to the shared network parameters. The gradient norms for each task are: In the formula, Represents the gradient norm of the data task; Indicates network shared parameters Find the gradient; , and They represent Data fitting loss function, mass conservation physical constraint loss function, and Bond breakage energy balance constraint loss function; Represents the gradient norm of the quality task; Represents the energy task gradient norm; The relative training rates for each task are: In the formula, Indicates the relative training rate of the quality task; Indicates the relative training rate of the energy task; , and These represent the initial data fitting loss function, the mass conservation physical constraint loss function, and the Bond breakage energy balance constraint loss function, respectively. The goal of GradNorm weight updates is: In the formula, Let γ be the mean of the gradient norms for each task, and γ be the asymmetric hyperparameter with a value of 1.5. After each gradient step, the constraint weights are updated according to the following rules: In the formula, or w The learning rate is updated for the weights, with a value of 0.025. C w This is the normalization constant for the sum of the constraint weights, and its value is 1. Through the above dynamic adjustment, PINN can automatically reduce the physical constraint strength in the early stage of training when the data fitting loss is large, and automatically increase the physical constraint strength in the later stage of training when the data fitting tends to stabilize, thus achieving gradient balance throughout the training process.

[0037] S3.4 The training strategy for the physical information neural network is as follows: The Adam optimizer was used, with an initial learning rate set to 10. -3 , batch size =32; Introducing a cosine annealing learning rate scheduling strategy to train the network, with a maximum number of training rounds. max poch =500; During the process, an Early Stopping mechanism is used to monitor the loss function value on the validation set, and the patience value is... P patience The value is 50.

[0038] S4, Deep Integration Uncertainty Quantification; The uncertainty of the prediction results is quantified using a deep ensemble method: training Two identical PINN member networks, but with random initialization, are fused. The predicted output of each member network is used to estimate the predicted distribution of the crushing efficiency through the mean and variance, and the expected value and probability prediction interval are output. S4.1 Training A member PINN network (in this invention) =5), each member network has the same structure but uses different random seeds to initialize parameters, and the dataset is randomly shuffled in different orders during training to ensure the prediction diversity of the member networks; S4.2, For each member network ( =1,2,..., The output layer uses a mean-variance parameterization method and outputs the predicted mean of crushing efficiency. and prediction variance (The variance is ensured to be non-negative using the Softplus activation function), and each member network is trained using negative log-likelihood loss instead of mean squared error loss: In the formula, Indicates the first The negative log-likelihood loss of a network; Indicates the first The official value of each sample is obtained through testing to determine the true crushing efficiency. Indicates the first The network is for the first The predicted mean of each sample, where For the first The predicted mean of the nth sample, i.e. The fusion feature vector corresponding to each test; S4.3, Fusion The predicted outputs of each member network are approximated using an equal-weight Gaussian mixture method to calculate the ensemble prediction distribution, where the ensemble prediction mean is... or ens and integrated prediction variance They are respectively: In the formula, Indicates the first The average prediction of each network; Indicates the first The prediction variance of the network; S4.4, Integrated Prediction Variance It contains two parts of uncertainty, namely: Cognitive uncertainty, reflected by the dispersion of the predicted means among members, characterizes the uncertainty of the model's knowledge. In the formula, This indicates cognitive uncertainty; Random uncertainty is reflected by the mean of the variances of the predictions of each member, characterizing the inherent randomness of data noise, etc. The total forecast uncertainty is the sum of the two: .

[0039] S4.5, with or ens As the final predicted value for crushing efficiency; under the assumption of normal distribution, with As a probability prediction interval with a 95% confidence level, where, Indicates the standard deviation of the integrated forecast; S4.6, Define the indicator of the proportion of cognitive uncertainty as follows: when R epi When the accuracy is >60%, it is determined that the current prediction is mainly driven by insufficient model knowledge, and simulation or test data for the corresponding working condition area should be supplemented; when s ens / or ens When the accuracy rate is greater than 15%, a warning about the reliability of the prediction will be issued.

[0040] S5, CUSUM control chart drift detection; The model prediction error sequence is continuously monitored using the Cumulative Sum (CUSUM) control chart method. When the CUSUM statistic exceeds the control threshold, it is determined that the ore properties have drifted significantly, triggering a two-layer drift correction process. The specific steps are as follows: S5.1 Calculate the relative prediction error at each test calibration time. : The error sequence is standardized using the mean of the historical steady-state window. m 0 and standard deviation s 0 is the baseline (in this embodiment, the values ​​are 2.8% and 1.2% respectively): In the formula, Indicates the first The relative error of the prediction time for the next test; Indicates the first Model predictions at the time of the next test; This represents the standardization error.

[0041] S5.2 Calculate the two-sided CUSUM statistic to detect positive drift (systematic overestimation) and negative drift (systematic underestimation) of the mean error, respectively: In the formula, Indicates the first Positive CUSUM statistic at the time of the test; Indicates the first The negative CUSUM statistic at the time of the test; Indicates the first Positive CUSUM statistic at the time of the test; Indicates the first The negative CUSUM statistic at the time of the test; As a reference value, 0.5 is used to detect mean shifts greater than 1σ; initial value of CUSUM. = =0.

[0042] S5.3 When the CUSUM statistic exceeds the control limit At that time, it was determined that the properties of the ore had undergone a significant shift: like > or > If this occurs, a drift warning will be triggered.

[0043] In the formula, For control limits (in this invention) =4, corresponding to average running length ARL 0. Under normal operating conditions without ore property drift, the time interval between approximately 500 false alarms in laboratory calibrations.

[0044] S5.4 After the drift alarm is triggered, record the drift start time. Reset the CUSUM statistic and initiate the two-layer drift correction process in step S6.

[0045] S6, Double-layer first-order drift correction; Upon triggering drift detection, a two-layer online drift correction is executed. The first layer uses Kalman filtering to estimate the state and dynamically update the bias parameters of the PINN output layer. The second layer constructs an incremental training dataset and updates all PINN parameters in a limited number of rounds using a forgetting factor-weighted loss function. Specific steps include: S6.1, First Layer: Fast Tracking and Correction of Kalman Filter Output Layer Parameters The bias parameter vector of the PINN output layer b out The hidden states are modeled as a random walk process, and online dynamic estimation is performed using a linear Kalman filter: Kalman state transition equations (prediction phase): In the formula, This represents the prior bias estimate; Indicates the first The posterior bias estimate after correction; Represents the prior error covariance matrix; Indicates the first The posterior error covariance matrix; The process noise covariance matrix (reflecting the parameter drift rate, with values ​​of diag(10)) represents the process noise covariance matrix. -4 ); Kalman update equation (correction phase): Calculate the observation residuals based on the newly collected laboratory calibration data: In the formula, Indicates the first The observation residuals of the step; Indicates the first The true crushing efficiency of the test; No. Predicted crushing efficiency from the first test; Kalman gain : In the formula, Indicates the first Kalman gain of the step; This represents the transpose of the observation matrix; The observation matrix (the Jacobian matrix of the output layer bias); This represents the observation noise covariance, taken as 0.01; Posterior estimation of bias parameters and update of error covariance: In the formula, Indicates the first Posterior bias estimation of the step; Represents the identity matrix; Update the bias parameters By directly replacing the PINN output layer bias, rapid parameter correction without retraining is achieved. After Kalman filter correction, the system can use the updated parameters for prediction online. It does not rely on the batch accumulation of new samples. The optimal linear unbiased estimate of the output layer bias can be updated based on a single new test data, shortening the parameter correction response time to within a single test cycle (8 hours).

[0046] S6.2, Second Layer: Fine-tuning of the weighted increment of the forgetting factor Once a drift alarm is triggered, newly added drift region samples are continuously accumulated. When the number of new samples reaches the batch threshold... N batch When =10, perform a full parameter fine-tuning, which is represented as: In the formula, This represents the total loss function for incremental fine-tuning; This represents the forgetting factor (initial value is 0.8, decays to 0.9 times the previous value after each trigger). Represents the loss on the core historical dataset, where, This represents the core historical dataset; This represents the loss in the newly added drift region data, where, This indicates a newly added drift region dataset; This represents the EWC regularization coefficient, set to 0.01; This represents all parameters of the model; This represents the model parameters before fine-tuning; This indicates a regularization term for Elastic Weight Consolidation (EWC) to prevent catastrophic forgetting. Using the current model parameters as initial values, a relatively low learning rate is adopted. lr finetune =10 -4 Perform training updates in a limited number of rounds (20 rounds).

[0047] S6.3, Dual-layer calibration collaborative mechanism; The first layer of Kalman filter correction takes effect in the first test cycle after drift (approximately 8 hours), quickly correcting the systematic bias; the second layer of incremental fine-tuning provides deep adaptation to the new operating condition distribution after accumulating 10 samples; the two layers of correction work together to achieve parameter drift compensation on both fast and slow time scales.

[0048] S6.4 If the CUSUM statistic continues to exceed the control limit after three consecutive drift triggers, a major change alarm for ore properties will be sent, and the full DEM simulation case extension will be re-executed to supplement the simulation knowledge base.

[0049] S7. After training, read real-time multi-source data of the grinding process and output the final expected value of the crushing efficiency and its probability prediction interval. The results are as follows: Figure 4 The bandwidth of the probability prediction interval shown intuitively reflects the level of total prediction uncertainty; Figure 5 The exhibit U epi and U ale As time changes, when U epiWhen there is a sudden increase, the system will automatically prompt you to supplement the corresponding operating condition data.

[0050] Specific embodiments of the invention have now been described. Other embodiments are within the scope of the appended claims. In some cases, the described actions can be performed in a different order and still achieve the desired result. Furthermore, the processes depicted in the drawings do not necessarily require a specific or sequential order to achieve the desired result. In some embodiments, multitasking and parallel processing can be advantageous.

[0051] Although embodiments of this application have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting this application. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of this application.

Claims

1. A method for predicting crushing efficiency in grinding processes based on PINN and drift correction, characterized in that, Includes the following steps: S1. Construct a multimodal data system that includes discrete element simulation of the ore crushing process, real-time acquisition of the grinding process by sensor array, and manual ore discharge particle size analysis. S2. Perform unified alignment and fusion operations on the multimodal data system obtained in S1 in time and feature space to generate fused feature vectors; S3. Input the fused feature vector obtained in S2 into the physical information neural network PINN with a physical constraint weight adjustment mechanism to predict the breaking efficiency, and obtain the trained network. S4. Construct a preset number of network structures that are identical to the network in S3 but satisfy different random initializations. Estimate the predicted distribution of the crushing efficiency using the mean and variance, and output the expected value and probability prediction interval. S5. Continuously monitor the prediction error sequence of S4 based on the CUSUM control chart method; S6. When the CUSUM statistic exceeds the control threshold, it is determined that the ore properties have drifted, triggering the two-layer drift correction process. S7. After training is completed, read the real-time multi-source data of the grinding process and output the final expected value of crushing efficiency and its probability prediction range. Based on S1 to S7, a method for predicting crushing efficiency in the grinding process was constructed.

2. The method for predicting crushing efficiency in grinding processes based on PINN and drift correction according to claim 1, characterized in that: In step S1, the data acquisition steps for discrete element simulation of the ore crushing process include: First, a three-dimensional geometric model including the cylinder and liners uniformly distributed along the internal axis of the cylinder is created using SolidWorks. The model is exported as a .STL file and imported into EDEM. The BPM model is used to construct aggregates formed by bonding bonds at a preset size to represent the ore. Then, material properties are set for the cylinder, liners, ore, and steel balls generated by the simulation software. Secondly, the normal and tangential contact forces between particles were calculated using the Hertz-Mindlin model for the established geometric model; the rolling friction torque was calculated using the Standard Rolling Friction model; the bond breakage mechanism in the ore was simulated using the particle bonding model BPM; and the physical rules for a single simulation were defined by randomly stacking ore and steel balls in the model through initialization settings. Next, based on the physical rules of a single simulation, multiple sets of working conditions were obtained through orthogonal design to provide experimental schemes for the simulation process. Among them, the orthogonal experiment was set as a three-factor, three-level orthogonal experiment. The three factors included: cylinder rotation speed, total filling rate, and steel ball filling rate. The three levels included: cylinder rotation speed set to 14, 17, and 20 r / min; total filling rate set to 25%, 35%, and 45%; and steel ball filling rate set to 8%, 11.5%, and 15%. Finally, by performing DEM-BPM simulation, the number of bond fractures was recorded, the number of bond fractures was converted into fracture efficiency, and the simulation dataset was stored.

3. The method for predicting crushing efficiency in grinding processes based on PINN and drift correction according to claim 1, characterized in that: In step S1, the real-time data acquisition step of the sensor array for the grinding process includes: using preset type sensors to acquire five types of high-frequency signals: feed rate, water supply rate, return sand rate, discharge slurry flow rate, and discharge slurry concentration.

4. The method for predicting crushing efficiency in grinding processes based on PINN and drift correction according to claim 1, characterized in that: In step S1, the data acquisition steps for the manual ore discharge particle size analysis process include: Samples of discharged ore are collected at a fixed sampling period. The Bond work index is calculated based on the sieve aperture diameter of particles passing through a preset percentage in the discharged ore and the sieve aperture diameter of particles passing through a preset percentage in the raw ore. The calculated Bond work index is used as the numerator, and the unit mass power consumption measured on-site is used as the denominator to calculate the true crushing efficiency. All test data are stored as a calibration dataset according to the timestamp.

5. The method for predicting crushing efficiency in grinding processes based on PINN and drift correction according to claim 1, characterized in that: Step S2 includes: First, an alignment window is defined with the test time as the main anchor point. Sensor statistical features are aggregated within this window. The most recent simulation conditions are found and the simulation conditions closest to the test time are extracted. The spliced ​​features are obtained by splicing the data. Secondly, the splicing features are standardized using Z-Score. Finally, a multi-head attention mechanism is used to output the fused feature vector after projection of the spliced ​​features.

6. The method for predicting crushing efficiency in a grinding process based on PINN and drift correction as described in claim 5, characterized in that: The steps in S3 include: First, a PINN network derived from a fully connected feedforward neural network is constructed, with the following settings: the input layer dimension is the same as the fused feature vector dimension; 6 feature extraction layers, each containing a linear transformation, BatchNorm, and ReLU activation function; a physical constraint embedding branch is set between the 3rd and 4th layers to calculate the mass conservation constraint residual and the energy balance constraint residual; the output layer outputs the predicted value. Secondly, the computational network includes: a total loss function comprising the data fitting loss function, the mass conservation constraint loss function, and the Bond breakage energy balance constraint loss function, expressed as follows: In the formula, L data This represents the loss function for data fitting. L mass Represents the loss function under the mass conservation constraint; L energy This represents the loss function for the energy balance constraint of Bond breakage. α (t) and β (t) represent the first and second dynamic weights, respectively; Next, using the GradNorm algorithm, by calculating the gradient norm and relative training rate for each loss task, the optimization objective is set to make the gradient norm of each loss task equal to the average gradient norm multiplied by the relative training rate, as shown in the following expression: In the formula, Represents the gradient norm of the quality task; This represents the mean of the gradient norms for each task; Indicates the relative training rate of the quality task; Represents the energy task gradient norm; This represents the relative training rate of the energy task; γ is an asymmetric hyperparameter with a value of 1.

5. Then, after each gradient step, the constraint weights are updated according to the following rules: In the formula, η w The learning rate is updated for the weights, with a value of 0.

025. C w This is the normalization constant for the sum of the constraint weights, and its value is 1. Finally, the Adam optimizer is used for training.

7. The method for predicting crushing efficiency in a grinding process based on PINN and drift correction as described in claim 6, characterized in that: The steps in S4 include: Five PINN networks with the same structure but different initialization parameters were constructed, and the training data were randomly shuffled to ensure diversity. For each PINN network, the predicted mean and predicted variance are simultaneously output, and it is trained with negative log-likelihood loss. The ensemble prediction distribution is calculated by fusing the prediction mean and prediction variance of each PINN network using an equal-weight Gaussian mixture approximation. This distribution includes the ensemble prediction mean and ensemble prediction variance. The ensemble prediction variance comprises cognitive uncertainty and random uncertainty, expressed as: In the formula, This indicates cognitive uncertainty. Indicates randomness and uncertainty; The integrated prediction mean is used as the final predicted value for the crushing efficiency point; under the assumption of normal distribution, the mean is used as the final predicted value for the crushing efficiency point. As a probability prediction interval with a 95% confidence level, where, Indicates the standard deviation of the integrated forecast; η ens Indicates the ensemble forecast mean; The indicator for calculating the proportion of cognitive uncertainty is: in, Indicates the variance of the integrated predictions; when R epi When the accuracy is >60%, it is determined that the current prediction is mainly driven by insufficient model knowledge, and simulation or test data for the corresponding working condition area should be supplemented. when σ ens / η ens When the accuracy rate is greater than 15%, a warning about the reliability of the prediction will be issued.

8. The method for predicting crushing efficiency in grinding processes based on PINN and drift correction according to claim 7, characterized in that: The steps in S5 include: Calculate the predicted relative error at each test calibration time. : The error sequence is standardized using the mean of the historical steady-state window. μ 0 and standard deviation σ 0 as the baseline: In the formula, Indicates the first The relative error of the prediction time for the next test; Indicates the first Model predictions at the time of the next test; Indicates standardized error; Calculate the two-sided CUSUM statistic to detect the positive and negative drift of the error mean, as shown in the following expression: In the formula, Indicates the first Positive CUSUM statistic at the time of the test; Indicates the first The negative CUSUM statistic at the time of the test; Indicates the first Positive CUSUM statistic at the time of the test; Indicates the first The negative CUSUM statistic at the time of the test; As a reference value, 0.5 is used to detect mean shifts greater than 1σ; initial value of CUSUM. = =0; When the CUSUM statistic exceeds the control limit At that time, it was determined that the properties of the ore had undergone a significant shift: like > or > If so, a drift warning will be triggered; To control the limit, we take 4, which corresponds to the average running length under normal operating conditions without ore property drift. After the drift alarm is triggered, the drift start time is recorded, the CUSUM statistic is reset, and the two-layer drift correction process in step S6 is started.

9. The method for predicting crushing efficiency in a grinding process based on PINN and drift correction as described in claim 8, characterized in that: The two-layer drift correction process of S6 is as follows: First layer: Use Kalman filtering to perform state estimation and dynamic tracking update of the bias parameters of the PINN output layer; The second layer involves constructing an incremental training dataset and updating all parameters of PINN in a limited number of rounds using a forgetting factor-weighted loss function.

10. The method for predicting crushing efficiency in a grinding process based on PINN and drift correction according to claim 9, characterized in that: Step S6 includes: The bias parameter vector of the PINN output layer is modeled as the hidden state of a random walk process, and Kalman filtering is used for online dynamic estimation of the optimal value, including: In the prediction phase, the Kalman state transition equation is expressed as: In the formula, This represents the prior bias estimate; Indicates the first The posterior bias estimate after correction; Represents the prior error covariance matrix; Indicates the first The posterior error covariance matrix; Represents the process noise covariance matrix; During the correction phase, the Kalman update equations include: Calculate the observation residuals based on the newly collected laboratory calibration data: In the formula, Indicates the first The observation residuals of the step; Indicates the first The true crushing efficiency of the test; No. Predicted crushing efficiency from the first test; Calculate Kalman gain : In the formula, Indicates the first Kalman gain of the step; This represents the transpose of the observation matrix; The observation matrix; This represents the observation noise covariance, taken as 0.01; Posterior estimation of bias parameters and update of error covariance: In the formula, Indicates the first Posterior bias estimation of the step; Represents the identity matrix; Update the bias parameters By directly replacing the PINN output layer bias, fast parameter correction can be achieved without retraining. Forgetting factor weighted incremental fine-tuning: After a drift alarm is triggered, newly added drift region samples are continuously accumulated. When the number of new samples reaches the batch threshold of 10, a full parameter fine-tuning is performed, as follows: In the formula, This represents the total loss function for incremental fine-tuning; This represents the forgetting factor, with an initial value of 0.8, which decays to 0.9 times the previous value after each trigger. Represents the loss on the core historical dataset, where, This represents the core historical dataset; This represents the loss in the newly added drift region data, where, This indicates a newly added drift region dataset; This represents the EWC regularization coefficient, set to 0.01; This represents all parameters of the model; This represents the model parameters before fine-tuning; This indicates a regularization term for the fixed elastic weights; Using the current model parameters as initial values, set the learning rate to 10. -4 Perform training updates in a limited number of rounds; The dual-layer calibration mechanism is set up as follows: the first layer takes effect in the first test cycle after drift; the second layer provides depth adaptation to the new working condition distribution after accumulating 10 samples. If the CUSUM statistic continues to exceed the control limit after three consecutive drift triggers, a major change alarm for ore properties will be sent and the full DEM simulation case extension will be re-executed to supplement the simulation knowledge base.

Images