A method and system for multi-target collaborative control of powder grinding

Abstract

The application discloses a kind of powder grinding multi-objective collaborative control method and system, belong to industrial process automation control field, it includes extracting powder concentrator load and mill main motor power time series data, based on variation mode decomposition separates internal mode component and generates working condition quality mapping matrix in combination with residence time distribution;Least square polynomial regression and numerical domination comparison are executed to working condition quality mapping matrix to construct candidate control parameter;The instantaneous frequency variance of internal mode component is calculated to construct dynamic weight, select the candidate vector of minimum weighted sum and execute amplitude truncation to generate control instruction;New data is collected to calculate prediction residual vector, and least square iteration update is executed to generate revised regression equation;The abnormal row of working condition quality mapping matrix is removed using the upper quantile boundary of prediction residual sequence, new data is added and secondary calibration is executed.The application uses data-driven closed-loop control, and the adaptability and control precision of model to nonlinear working condition can be improved.

Description

Technical Field

[0001] This invention relates to the field of industrial process automation control, and in particular to a multi-objective collaborative control method and system for powder grinding. Background Technology

[0002] The powder grinding process typically includes grinding, classification, and material recycling. The production site generally includes the main mill motor, classifier, circulating fan, feeding device, and regulating valves. Under a distributed control architecture, data is collected during operation. The collected data may include classifier load time-series data, main mill motor power time-series data, classifier speed, circulating fan speed, and total feed rate, among other process variables. Product quality indicators can be characterized using offline particle size analysis data, and energy consumption indicators can be characterized using power per unit output.

[0003] In related technologies, Chinese patent CN121050232A discloses a method for operating and controlling a pulverizing unit in conjunction with a coal-fired boiler during peak shaving in thermal power plants. Specifically, it includes: determining a first typical particle size distribution range when the coal-fired boiler is in normal operation and a second typical particle size distribution range when it is in peak shaving operation; establishing a mathematical intelligent model between coal particle size and the performance of the thermal power generation system; configuring the number and operating parameters of the pulverizing mills; and, when the thermal power generation system is in normal or peak shaving operation, determining and starting the corresponding number of pulverizing mills based on the calculated coal particle size requirements, combined with real-time monitoring data and the optimization results of the mathematical intelligent model, and dynamically adjusting the operating parameters of the pulverizing mills.

[0004] However, when using offline quality indicators and online process variables for the coordinated control of powder grinding processes, constraints exist, including inconsistencies in the time bases of offline detection and online process variables, noise interference from process variables, and the difficulty in uniformly expressing multi-objective trade-offs. Offline particle size detection obtains results through sampling, conveying, sample preparation, and detection, resulting in a time lag between the detected value and the sampling time. Simultaneously, the residence time of materials in the grinding loop exhibits a distributed pattern, with different batches of samples arriving at the classification loop at discrete times, making the direct correspondence between offline particle size detection values ​​and online process variables unstable. Time-series signals such as classifier load and mill main motor power contain trend and fluctuation terms, accompanied by measurement noise or transient disturbances. When the raw signals are directly used for modeling, the regression parameters are sensitive to noise. Unit output power and particle size deviation are objectives with different dimensions, and the emphasis of the trade-off changes under different operating conditions. Fixed weights or fixed setpoints are insufficient to cover changes in operating conditions. Furthermore, changes in raw material characteristics, equipment wear, sensor drift, or intermittent failures can cause shifts in the mapping relationship between process variables and quality and energy consumption indicators. When there is a lack of model updates based on the latest data and abnormal sample cleaning, regression bias accumulates over multiple control cycles, affecting the selection of control parameters and the instruction generation process. Summary of the Invention

[0005] To address the aforementioned issues, this invention provides a multi-objective collaborative control method and system for powder grinding, employing data-driven closed-loop control, which can improve the model's adaptability to nonlinear operating conditions and control accuracy.

[0006] The above objectives can be achieved through the following approach:

[0007] A multi-objective collaborative control method for powder grinding includes extracting the time-series data of the classifier load and the mill main motor power from a distributed control system; performing signal denoising processing based on variational mode decomposition to separate the internal model component; acquiring offline particle size detection data from an automated laboratory; performing convolution mapping on the internal model component and the offline particle size detection data using the residence time distribution of materials in the grinding system to generate a working condition quality mapping matrix; performing least squares polynomial regression on the working condition quality mapping matrix, using the classifier speed, circulating fan speed, and total feed rate as independent variables, and the unit output power value and the sum of squares of particle size difference as dependent variables to obtain a regression equation; sampling within the domain of the independent variables to generate a test vector set and substituting it into the regression equation to calculate predicted values; performing numerical dominance comparison to remove dominated vectors and extracting undominated vectors to construct candidate control parameters; calculating the instantaneous frequency variance of the internal model component, and constructing dynamic weights for the unit output power value and the sum of squares of particle size difference. The sum of squared differences is weighted and summed to select the candidate control vector with the smallest weighted sum. The differential increment is calculated using the candidate control vector and the current setpoint, and amplitude truncation is performed to generate inverter frequency and valve opening commands for execution. The inverter frequency and valve opening commands are executed to collect new unit output power values ​​and new particle size detection values. The predicted residual vectors of the new unit output power values ​​and the new particle size detection values ​​relative to the regression equation are calculated, and the regression equation is updated using least squares iteration to generate a modified regression equation. The predicted residual sequence of historical data is extracted using the modified regression equation, and the upper quantile boundary of the predicted residual sequence within the sliding window is calculated. The operating condition quality mapping matrix is ​​traversed, and abnormal data rows with residual magnitudes exceeding the upper quantile boundary are removed. The new unit output power values ​​and the new particle size detection values ​​are appended to the operating condition quality mapping matrix. Least squares regression is performed to perform secondary calibration of the modified regression equation and archive it.

[0008] Optionally, the generation of the operating condition quality mapping matrix includes: performing variational mode decomposition on the time-series data of the classifier load and the time-series data of the mill main motor power, outputting a set of modal components, calculating the cross-correlation coefficient between each component in the set of modal components and the original time-series data, and locking the component with the highest cross-correlation coefficient as the inner mode component; extracting the sampling time from the offline particle size detection data, calculating the cross-correlation function between the inner mode component and the offline particle size detection data, and taking the time shift of the maximum value to generate the average lag time, constructing a residence time distribution function with the average lag time as the expectation and performing discretization processing to obtain a time lag weight sequence; using the sampling time minus the average lag time as an index, extracting historical time window data with the same length as the time lag weight sequence from the inner mode component, performing a weighted summation operation on the historical time window data and the time lag weight sequence to obtain a weighted operating condition feature value, and concatenating the weighted operating condition feature value with the offline particle size detection data to generate the operating condition quality mapping matrix.

[0009] Optionally, the step of extracting undominated vectors to construct candidate control parameters includes: extracting data columns based on the classifier speed, circulating fan speed, and total feed amount in the working condition quality mapping matrix; constructing a polynomial extended design matrix; and performing matrix operations in conjunction with the unit output power value and the sum of squared particle size differences to calculate a regression coefficient vector to establish a regression equation; extracting the domain boundary values ​​of the independent variables; performing equal-step grid scanning within the domain boundary values ​​to generate a test vector set; substituting the test vector set into the regression equation to generate a set of predicted values; traversing the set of predicted values ​​to perform pairwise numerical comparisons; if there exists a first vector whose unit output power value and the sum of squared particle size differences are both less than the second vector, then the second vector is marked as a dominated vector and removed, and the unmarked vectors are retained as candidate control parameters.

[0010] Optionally, the calculation of the regression coefficient vector includes: calculating the squared vector and the product vector of the data column; performing column vector concatenation on the data column, the squared vector, and the product vector to generate a polynomial extended design matrix; performing operations on the polynomial extended design matrix to obtain the inverse matrix; and performing chain matrix multiplication on the inverse matrix, the transpose of the polynomial extended design matrix, the unit output power value, and the sum of squares of the particle size difference to generate the regression coefficient vector.

[0011] Optionally, the generation and execution of the inverter frequency and valve opening command includes: counting the number of zero-crossings of the internal model component within the sampling time; dividing the sampling time by the number of zero-crossings to obtain the average period; calculating the variance of the amplitude of each point of the internal model component relative to the average period as the instantaneous frequency variance; mapping the instantaneous frequency variance to a normalized value using the hyperbolic tangent function as a dynamic weight for the sum of squares of the granularity difference; calculating the mean of the unit output power value and the mean of the sum of squares of the granularity difference using the candidate control parameters; and mapping the unit output power value to a normalized value using the hyperbolic tangent function as a dynamic weight for the sum of squares of the granularity difference. The unit output power value and the sum of squares of the particle size difference are respectively divided by the mean of the corresponding unit output power value and the mean of the sum of squares of the particle size difference to generate a dimensionless ratio. The dimensionless ratio is weighted and summed using the dynamic weight, and the vector with the smallest summation result is taken as the preferred vector. The standard deviation of the classifier speed is calculated using the working condition quality mapping matrix as a safety fluctuation threshold. The difference between the preferred vector and the current set value is calculated. If the absolute value of the difference is greater than the safety fluctuation threshold, the frequency converter frequency and valve opening command are generated.

[0012] Optionally, the method further includes: performing matrix multiplication with the polynomial extended design matrix and the regression coefficient vector to generate a theoretical recurrence value vector; calculating the difference vector between the theoretical recurrence value vector and the unit output power value, and defining the difference vector as a static fitting deviation vector.

[0013] Optionally, generating the modified regression equation includes: substituting the inverter frequency and valve opening command into the regression equation to perform forward calculation to obtain a theoretical predicted value; calculating the difference between the new unit output power value and the new particle size detection value relative to the theoretical predicted value to obtain a prediction residual vector; constructing a vector autocorrelation matrix using the inverter frequency and valve opening command; calculating the inverse matrix of the vector autocorrelation matrix to generate an orthogonal projection gain matrix; performing chain multiplication on the inverter frequency and valve opening command, the prediction residual vector, and the orthogonal projection gain matrix to generate a coefficient correction matrix; extracting the current coefficient matrix for the regression equation; performing matrix addition on the current coefficient matrix and the coefficient correction matrix to reconstruct the regression equation to generate the modified regression equation.

[0014] Optionally, the secondary calibration and archiving of the modified regression equation includes: calculating the coefficient difference between the modified regression equation and the original regression equation; multiplying the coefficient difference with the polynomial extended design matrix to obtain the model drift vector; calculating the difference between the static fit deviation vector and the model drift vector to generate a prediction residual sequence; extracting the upper quantile of the absolute value of the prediction residual sequence as a cleaning boundary; removing row vectors with residual moduli greater than the cleaning boundary from the working condition quality mapping matrix; defining the remaining matrix data as a reliable sample matrix; appending the new unit output power value and the new granularity detection value to the reliable sample matrix; solving the inverse matrix for the reliable sample matrix; generating secondary calibrated regression coefficients; performing parameter coverage on the modified regression equation; and archiving the data.

[0015] Optionally, the archiving includes: counting the number of rows in the trusted sample matrix as the current sample size, extracting the window length of the sliding window as the capacity limit; calculating the current sample size minus the capacity limit as the overflow quantity; if the overflow quantity is greater than zero, then starting from the starting row index of the trusted sample matrix, sequentially removing row vectors equal to the overflow quantity.

[0016] Based on the same inventive concept, this invention also provides a multi-objective collaborative control system for powder grinding. The system includes: a working condition quality mapping construction module, used to extract the classifier load time-series data and mill main motor power time-series data from the distributed control system; perform signal noise reduction processing based on variational mode decomposition to separate the internal model component; obtain offline particle size detection data from an automated laboratory; and perform convolution mapping on the internal model component and the offline particle size detection data using the material residence time distribution within the grinding system to generate a working condition quality mapping matrix; a regression prediction and undominated screening module, used to perform least-squares polynomial regression on the working condition quality mapping matrix, using the classifier speed, circulating fan speed, and total feed rate as independent variables, and the unit output power value and the sum of squares of particle size differences as dependent variables to obtain a regression equation; sampling within the domain of the independent variables to generate a test vector set and substituting it into the regression equation to calculate predicted values; performing numerical dominance comparison to remove dominated vectors; and extracting undominated vectors to construct candidate control parameters; and a candidate vector sorting and instruction generation module, used to calculate the instantaneous frequency variance of the internal model component and construct dynamic... The state weights perform a weighted summation of the unit output power value and the sum of squares of the particle size difference, select the candidate control vector with the smallest weighted sum, calculate the differential increment using the candidate control vector and the current setpoint, and perform amplitude truncation to generate inverter frequency and valve opening commands for execution; the prediction residual calculation and regression update module is used to execute the inverter frequency and valve opening commands to collect new unit output power values ​​and new particle size detection values, calculate the prediction residual vector of the new unit output power value and the new particle size detection value relative to the regression equation, and for the given... The regression equation is updated using least squares iteration to generate a modified regression equation. The residual threshold calculation and outlier removal module is used to extract the predicted residual sequence of historical data using the modified regression equation, calculate the upper quantile boundary of the predicted residual sequence within the sliding window, traverse the working condition quality mapping matrix and remove outlier data rows whose residual modulus exceeds the upper quantile boundary, append the new unit output power value and the new granularity detection value to the working condition quality mapping matrix, perform least squares regression to perform secondary calibration of the modified regression equation and archive it.

[0017] Compared with the prior art, the present invention has the following advantages:

[0018] 1. Variational mode decomposition and residence time distribution convolution mapping are adopted to correlate the time-series signals of classifier load and mill power with offline particle size detection data under the same time reference to form a working condition quality mapping matrix, thereby reducing the impact of offline particle size data sampling lag on control calculation;

[0019] 2. A test vector set is generated based on least squares polynomial regression and domain sampling. Candidate control parameters are screened through numerical dominance comparison. Then, dynamic weights are generated from instantaneous frequency variance for weighted selection, realizing closed-loop linkage between control command generation and regression equation iterative update.

[0020] Other features and advantages of the invention will be set forth in the description which follows, and will be apparent in part from the description, or may be learned by practicing the invention. The objects and other advantages of the invention may be realized and obtained by means of the structures pointed out in the description, claims and drawings. Attached Figure Description

[0021] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0022] Figure 1 This is a schematic flowchart of a multi-objective collaborative control method for powder grinding according to an embodiment of the present invention.

[0023] Figure 2 This is a graph showing the consistency between the regression equation prediction and the actual measurement in an embodiment of the present invention.

[0024] Figure 3 This is a Taylor plot comparing the prediction performance of the regression equation before and after secondary calibration in an embodiment of the present invention.

[0025] Figure 4 This is a schematic diagram of the structure of a multi-objective collaborative control system for powder grinding according to an embodiment of the present invention. Detailed Implementation

[0026] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0027] Reference Figure 1 One embodiment of the present invention proposes a multi-objective collaborative control method for powder grinding, which adopts data-driven closed-loop control, and can improve the model's adaptability to nonlinear working conditions and control accuracy.

[0028] The method described in this embodiment specifically includes:

[0029] Extract the time-series data of the classifier load and the time-series data of the mill main motor power from the distributed control system, perform signal noise reduction processing based on variational mode decomposition, separate the internal model component, obtain the offline particle size detection data of the automated laboratory, and perform convolution mapping on the internal model component and the offline particle size detection data using the residence time distribution of the material in the grinding system to generate the working condition quality mapping matrix;

[0030] The least squares polynomial regression is performed on the working condition quality mapping matrix. The classifier speed, circulating fan speed and total feed amount are used as independent variables, and the unit output power value and the sum of squares of particle size difference are used as dependent variables to obtain the regression equation. Test vector sets are generated by sampling within the domain of the independent variables and substituted into the regression equation to calculate the predicted value. Numerical dominance comparison is performed to remove dominated vectors and extract non-dominated vectors to construct candidate control parameters.

[0031] Calculate the instantaneous frequency variance of the internal model component, construct dynamic weights to perform a weighted summation of the unit output power value and the sum of squares of the particle size difference, select the candidate control vector with the smallest weighted sum, calculate the differential increment through the candidate control vector and the current set value and perform amplitude truncation, generate inverter frequency and valve opening commands and issue them for execution;

[0032] The inverter frequency and valve opening commands are executed to collect new unit output power value and new particle size detection value. The prediction residual vector of the new unit output power value and the new particle size detection value relative to the regression equation is calculated. The least squares iterative update is performed on the regression equation to generate a modified regression equation.

[0033] The predicted residual sequence of historical data is extracted using the modified regression equation. The upper quantile boundary of the predicted residual sequence within the sliding window is calculated. The working condition quality mapping matrix is ​​traversed and abnormal data rows with residual modulus exceeding the upper quantile boundary are removed. The new unit output power value and the new granularity detection value are appended to the working condition quality mapping matrix. Least square regression is performed to perform secondary calibration of the modified regression equation and archive it.

[0034] Optionally, the generated working condition quality mapping matrix includes:

[0035] Variational mode decomposition is performed on the time series data of the classifier load and the time series data of the mill main motor power to output a set of modal components. The cross-correlation coefficient between each component in the set of modal components and the original time series data is calculated, and the component with the first cross-correlation coefficient is locked as the inner mode component.

[0036] Variational mode decomposition (VM) was performed on the classifier load time-series data and the mill main motor power time-series data. The number of decomposition modes and the penalty factor were set to decompose the non-stationary time-series signals into multiple intrinsic mode functions (IMFs) with different center frequencies. Using the Pearson correlation coefficient formula, the correlation coefficient between each IMF and the original time-series data was calculated. All calculated correlation coefficients were sorted in descending order of absolute value. The IMF corresponding to the correlation coefficient with the largest absolute value was selected and defined as the internal mode component.

[0037] For example, the decomposition level is set to 5, and the penalty factor is 2000. Mill main motor power data is collected for 10 minutes and decomposed into 5 component signals. The correlation coefficients between these 5 components and the original data are calculated, and the results are 0.85, 0.42, 0.15, 0.08, and 0.02, respectively. Since 0.85 is the maximum value, the first component corresponding to 0.85 is locked as the internal mold component.

[0038] The sampling time is extracted from the offline granularity detection data. The cross-correlation function between the internal model component and the offline granularity detection data is calculated, and the time shift of the maximum value is used to generate the average lag time. The dwell time distribution function is constructed with the average lag time as the expectation and discretization is performed to obtain the time lag weight sequence.

[0039] The time stamps corresponding to the offline granularity detection data are read as sampling times. Cross-correlation is performed on the internal model component data sequence and the offline granularity detection data sequence to obtain the cross-correlation function curve. The point with the largest peak in the cross-correlation function curve is identified, and the corresponding time displacement value is read. This time displacement value is defined as the average lag time. Using this average lag time as the expected value, a dwell time distribution function is constructed. By employing a Gaussian distribution function, the dwell time distribution function is discretized and sampled within a fixed time range to generate a vector sequence composed of numerical values, i.e., the time lag weight sequence.

[0040] Using the sampling time minus the average lag time as an index, historical time window data with the same length as the time lag weight sequence is extracted from the internal model component. A weighted summation operation is performed on the historical time window data and the time lag weight sequence to obtain a weighted operating condition feature value. The weighted operating condition feature value is then concatenated with the offline granularity detection data to generate an operating condition quality mapping matrix.

[0041] A traceability index time point is obtained by subtracting the average lag time from the sampling time. Using this traceability index time point as a reference, a segment of data is extracted from the internal model component data stream, ensuring that the length of the extracted data is exactly equal to the number of elements in the time lag weight sequence. This segment of data constitutes the historical time window data. Each value in the historical time window data is multiplied by the corresponding value in the time lag weight sequence. All multiplication results are summed to obtain a unique scalar value as the weighted operating condition feature value. This weighted operating condition feature value is used as the feature variable, and the corresponding offline granular detection data is used as the target variable. These two values ​​are combined into a single data record and stored in the operating condition quality mapping matrix.

[0042] For example, the sampling time for the granular data is 2:00 PM. The calculated average lag time is 20 minutes. The constructed dwell time distribution function is discretized to generate a weighted sequence of 5 minutes, with values ​​of 0.1, 0.2, 0.4, 0.2, and 0.1 respectively. Subtracting 20 minutes from 2:00 PM positions the data to 1:40 PM. Data from 1:38 PM to 1:42 PM is extracted from the internal model component, with values ​​of 100, 102, 105, 103, and 101 respectively. A weighted summation operation is performed: 100 multiplied by 0.1, plus 102 multiplied by 0.2, plus 105 multiplied by 0.4, plus 103 multiplied by 0.2, plus 101 multiplied by 0.1, resulting in 103.1. The value 103.1 is concatenated with the granular data and stored.

[0043] Optionally, the extraction of non-dominated vectors to construct candidate control parameters includes:

[0044] Based on the classifier speed, circulating fan speed and total feed amount extracted from the working condition quality mapping matrix, a polynomial extended design matrix is ​​constructed, and matrix operations are performed in combination with the unit output power value and the sum of squared particle size differences to calculate the regression coefficient vector to establish the regression equation.

[0045] The classifier speed, circulating fan speed, and total feed rate are read from the working condition quality mapping matrix as independent variables. To capture nonlinear characteristics, feature expansion processing is performed on these three columns. The square of each column is calculated as a quadratic term feature, and the product of each pair of different columns is calculated as an interaction term feature. The original three columns of data, the calculated quadratic term feature data, and the interaction term feature data are concatenated column by column to construct a polynomial extended design matrix. Simultaneously, the unit output power value column and the sum of squares of particle size difference column are read from the working condition quality mapping matrix as the dependent variable matrix. The polynomial extended design matrix is ​​transposed, and the transposed matrix is ​​multiplied by the original matrix to obtain the normal equation matrix. The inverse matrix of this normal equation matrix is ​​calculated. This inverse matrix is ​​multiplied by the transposed matrix, and then multiplied by the dependent variable matrix. A regression coefficient vector is obtained through chain matrix operations. Each coefficient in the regression coefficient vector is multiplied by its corresponding independent variable feature term and summed to establish the regression equation used for prediction, such as... Figure 2 The chart shows the comparison between measured and predicted consistency of the regression equation. The horizontal axis represents measured values, and the vertical axis represents predicted values. The gray diagonal line indicates ideal consistency, and the light gray band represents the ±10% error range. Different symbols represent different operating condition zones, demonstrating the model's fitting stability across different load ranges.

[0046] Extract the domain boundary values ​​of the independent variable, perform equal-step grid scanning within the domain boundary values ​​to generate a test vector set, and substitute the test vector set into the regression equation to generate a set of predicted values.

[0047] The system reads the equipment parameter configuration file to obtain the minimum and maximum allowable operating values ​​for the classifier speed, circulating fan speed, and total feed rate, and defines these extreme values ​​as the domain boundary values. A scan step size is set for each independent variable. Using the domain boundary values ​​as the range limit, the system iterates incrementally through each independent variable according to the step size, generating all possible numerical combinations. Each numerical combination is encapsulated into a vector, forming a test vector set. Each vector in the test vector set is read sequentially and substituted into the previously established regression equation for forward calculation, outputting the corresponding predicted unit output power value and the predicted sum of squared particle size differences. These two prediction results are combined to form a set of predicted values.

[0048] For example, the defined boundary value of the classifier speed is 600 rpm to 800 rpm, with a scan step size of 10 rpm. The defined boundary value of the circulating fan speed is 1000 rpm to 1200 rpm, with a scan step size of 20 rpm. The defined boundary value of the total feed rate is 80 tons per hour to 100 tons per hour, with a scan step size of 5 tons per hour. A grid scan is performed within these ranges, generating 400 test vectors. Substituting the first test vector into the regression equation, the calculated unit output power is 45 kWh per ton, and the sum of squares of particle size differences is 0.05. Substituting the second test vector, the calculated unit output power is 48 kWh per ton, and the sum of squares of particle size differences is 0.08.

[0049] The predicted value set is traversed and pairwise numerical comparisons are performed. If the unit output power value and the sum of squared particle size differences of the first vector are both less than those of the second vector, then the second vector is marked as a dominated vector and removed, and the unmarked vectors are retained as candidate control parameters.

[0050] A double-loop iterative approach is used to perform a full traversal of the predicted value set. In the outer loop, a vector is selected sequentially from the predicted value set as a benchmark reference object, defined as the first vector. In the inner loop, another vector besides the benchmark reference object is selected from the predicted value set as the object to be judged, defined as the second vector. Numerical comparison logic is executed: the unit output power value of the first vector is compared with the unit output power value of the second vector, and the sum of squared granularity differences of the first vector is also compared with the sum of squared granularity differences of the second vector. If the unit output power value of the first vector is less than that of the second vector, and the sum of squared granularity differences of the first vector is also less than that of the second vector, then the second vector is determined to be a suboptimal solution. The second vector is marked as a dominated state. After completing the comparison of all possible combinations, all vectors marked as dominated states are identified and permanently deleted from the predicted value set. The remaining vectors in the set that have never been marked are retained; these vectors are the non-dominated vectors and are stored as candidate control parameters.

[0051] For example, the predicted value set includes vectors A, B, and C. In one comparison cycle, vector A is selected as the first vector, with a value of 42 kWh per ton for unit output power and a sum of squares of particle size difference of 0.03. Vector B is selected as the second vector, with a value of 45 kWh per ton for unit output power and a sum of squares of particle size difference of 0.06. The comparison shows that 42 is less than 45, and 0.03 is less than 0.06. Therefore, vector A is superior to vector B in both metrics, and vector B is determined to be dominated by vector A, so its label is removed. In another comparison, vector A (42, 0.03) is selected as the first vector, and vector C (40, 0.08) is selected as the second vector. The comparison shows that 42 is greater than 40, but 0.03 is less than 0.08. Since they do not dominate each other, no labeling operation is performed, and vector C is retained. Finally, vectors A and C are selected as candidate control parameters for the next stage.

[0052] Optionally, the calculation of the regression coefficient vector includes:

[0053] Calculate the squared vector and the product vector of the data column, and concatenate the data column, the squared vector and the product vector into column vectors to generate a polynomial extended design matrix;

[0054] The process involves reading the three columns of raw data—classifier speed, circulating fan speed, and total feed rate—stored in the working condition quality mapping matrix. For each column, a self-multiplication operation is performed to obtain three corresponding squared vectors. Then, pairwise multiplication operations are performed on the three columns: the product of classifier speed and circulating fan speed, the product of classifier speed and total feed rate, and the product of circulating fan speed and total feed rate, resulting in three cross-product vectors. An empty matrix is ​​constructed, with the three raw data columns as the base columns, the three squared vectors as the extension columns, and the three cross-product vectors as the interaction columns. These are then horizontally concatenated into the empty matrix in column order. Additionally, to handle the intercept term, a constant vector with all values ​​of one is appended before the first column of the matrix. The final matrix, containing columns of constant terms, raw data, squared terms, and cross-product terms, is defined as the polynomial extended design matrix.

[0055] For example, assume the original data contains 100 rows of samples. The extracted classifier speed column is vector X1, the circulating fan speed column is vector X2, and the total feed rate column is vector X3. Multiplying X1 by X1 yields the squared vector S1, and similarly yields S2 and S3. Multiplying X1 by X2 yields the product vector M1, X1 by X3 yields M2, and X2 by X3 yields M3. This generates a 100-row, 1-column vector C consisting entirely of 1s. Concatenating vectors C, X1, X2, X3, S1, S2, S3, M1, M2, and M3 horizontally in sequence generates a 100-row, 10-column two-dimensional array, which is the polynomial extended design matrix.

[0056] The inverse matrix is ​​obtained by performing operations on the polynomial extended design matrix, and a chain matrix multiplication is performed on the inverse matrix, the transpose of the polynomial extended design matrix, the unit output power value, and the sum of squared granularity differences to generate a regression coefficient vector.

[0057] Perform a matrix transpose operation on the polynomial extended design matrix, interchanging rows and columns to obtain the transpose matrix. Multiply this transpose matrix with the original polynomial extended design matrix to generate a square matrix, i.e., the normal equation matrix. Use Gaussian elimination to obtain the inverse matrix of this normal equation matrix. Read the unit output power value column and the particle size difference sum of squares column from the working condition quality mapping matrix to form the dependent variable matrix. Perform chained multiplication according to a specific multiplication order: first calculate the product of the inverse matrix and the transpose matrix, then multiply the product result by the dependent variable matrix. The final matrix obtained is the regression coefficient vector, where each row corresponds to a set of regression coefficients for one dependent variable.

[0058] For example, the polynomial extended design matrix is ​​A, with a dimension of 100 rows and 10 columns. The dependent variable matrix is ​​Y, containing two columns of data: unit output power value and sum of squared granularity differences, with a dimension of 100 rows and 2 columns. First, calculate the transpose matrix A of A. T The dimension is 10 rows and 100 columns. Calculate A. T The product of A and B yields the normal equation matrix B, with dimensions 10 rows and 10 columns. Calculate the inverse matrix Binverse. inv Calculate B inv With A T The product of P and Y yields the projection matrix P, which has dimensions of 10 rows and 100 columns. Finally, the product of P and Y is calculated to obtain a regression coefficient vector with dimensions of 10 rows and 2 columns. The first column contains 10 regression coefficients for the power value per unit output, and the second column contains 10 regression coefficients for the sum of squares of granularity differences.

[0059] Optionally, the generation and execution of the inverter frequency and valve opening command includes:

[0060] The number of zero-crossing points of the internal model component within the sampling time is counted. The sampling time is divided by the number of zero-crossing points to obtain the average period. The variance of the amplitude of each point of the internal model component relative to the average period is calculated as the instantaneous frequency variance. The instantaneous frequency variance is mapped to a normalized value using the hyperbolic tangent function as a dynamic weight for the sum of squares of the granularity difference.

[0061] Read the internal model component data sequence and identify the number of times the numerical sign changes in the sequence; define this number as the number of zero-crossings. Read the total sampling duration of the internal model components and divide the total sampling duration by the number of zero-crossings to calculate the average period. Calculate the amplitude of each data point in the internal model component sequence, calculate the variance of these amplitudes, and divide this variance by the average period to obtain the instantaneous frequency variance. Use the hyperbolic tangent function as the activation function, with the instantaneous frequency variance as the input variable, and calculate the output value. For positive inputs, since the output range of the hyperbolic tangent function is between zero and one, this output value is directly defined as the dynamic weight. This dynamic weight is specifically used to characterize the degree of importance attached to the product quality indicator, namely the sum of squared granularity differences.

[0062] For example, the sampling time is 10 seconds. The internal model component sequence changes from positive to negative or vice versa a total of 50 times within these 10 seconds, i.e., the number of zero crossings is 50. The average period is calculated to be 0.2 seconds. The statistical variance of the internal model component amplitude is 2.0. Dividing 2.0 by 0.2 seconds yields an instantaneous frequency variance of 10. Substituting the value 10 into the hyperbolic tangent function, we obtain 0.99. 0.99 is determined as the dynamic weight for the sum of squared granularity differences. This indicates that the current operating condition fluctuates drastically, and the control strategy will allocate 99% of its weight to ensuring granularity quality.

[0063] The average value of unit output power and the average value of the sum of squares of particle size difference are calculated using the candidate control parameters. The average value of unit output power and the average value of the sum of squares of particle size difference are divided by the corresponding average value of unit output power and the average value of the sum of squares of particle size difference to generate a dimensionless ratio. The dimensionless ratio is weighted and summed using the dynamic weights. The vector with the smallest summation result is selected as the preferred vector.

[0064] The candidate control parameter set is traversed. The unit output power value of all vectors is summed and divided by the total number of vectors to obtain the average unit output power value. The sum of squared granularity differences of all vectors is also summed and divided by the total number of vectors to obtain the average sum of squared granularity differences. For each candidate vector in the set, its unit output power value is divided by the average unit output power value to obtain the dimensionless power ratio; its sum of squared granularity differences is divided by the average sum of squared granularity differences to obtain the dimensionless granularity ratio. Using dynamic weights, a weighted evaluation formula is used to calculate the comprehensive evaluation value of each candidate vector. The formula for the weighted summation is as follows:

[0065] ,

[0066] in, Indicates the first The weighted sum of the candidate vectors is the comprehensive evaluation value. The smaller the value, the better the comprehensive performance of the control vector under the current working conditions. This represents the dynamic weight, with a value ranging from 0 to 1. This parameter originates from the fluctuation of the internal model component and is significant in adjusting the emphasis on quality and energy consumption. When the fluctuation is large, When the value approaches 1, the formula is mainly dominated by the mass term; when it is stable, As the value approaches zero, the formula is mainly dominated by the energy consumption term. Indicates the first The sum of squared differences in granularity of each candidate vector is derived from the predicted output of the regression equation and represents the deviation between product quality and the target. This represents the mean of the sum of squared granular differences among all candidate vectors. This denominator is introduced to eliminate dimensions and normalize physical quantities of different magnitudes to the same scale. Indicates the first The unit output power value of each candidate vector is derived from the predicted output of the regression equation and represents the energy consumption level. This represents the mean power per unit output of all candidate vectors; This represents the complementary weights for power values ​​per unit output. The formula uses mean normalization to eliminate the order-of-magnitude difference between power units and particle size units, ensuring that they can be directly added. This is achieved through dynamic weighting. The adjustment achieves the adaptive control objective of sacrificing energy consumption to ensure quality when the operating conditions are unstable, and sacrificing quality deviation to save energy when the operating conditions are stable.

[0067] For example, assume dynamic weights The power is 0.8. The value is 45, indicating a difference in particle size. The power is 0.02. The value is 40, and the particle size difference is significant. The mean power of all vectors is 0.08. The average particle size difference is 42.5. The score is 0.05. Therefore, the score for vector A is: The score for vector B is: Since 0.5316 is less than 1.4682, the weighted sum of vector A is smaller, therefore vector A is selected as the preferred vector.

[0068] The standard deviation of the classifier speed is calculated using the working condition quality mapping matrix as a safety fluctuation threshold. The difference between the preferred vector and the current set value is calculated. If the absolute value of the difference is greater than the safety fluctuation threshold, the frequency converter frequency and valve opening command are generated.

[0069] The historical classifier speed data column stored in the working condition quality mapping matrix is ​​extracted. The dispersion of this data column is statistically analyzed using the standard deviation calculation formula, and the calculated standard deviation value is defined as the safety fluctuation threshold. The currently executing inverter frequency and valve opening setpoints are read. The target parameter value contained in the optimized vector is subtracted from the current setpoint to obtain the difference. The absolute value of this difference is calculated and compared with the safety fluctuation threshold. If the absolute value of the difference is greater than the safety fluctuation threshold, it indicates that the adjustment range required by the optimized vector is too large, which may cause oscillations. In this case, the adjustment range is limited using the safety fluctuation threshold. The safety fluctuation threshold is multiplied by the sign of the difference, and the result is added to the current setpoint to generate the final inverter frequency and valve opening command. If the absolute value of the difference is less than or equal to the safety fluctuation threshold, the optimized vector is directly issued as the command.

[0070] For example, the standard deviation of the air classifier's rotational speed in historical data is 15 revolutions per minute (rpm), meaning the safety fluctuation threshold is 15. The current air classifier speed is 800 rpm. The preferred vector suggests adjusting the speed to 850 rpm. The calculated difference is 50 rpm. Since 50 is greater than 15, the safety limit logic is triggered. The generated new instruction is 800 plus 15, i.e., 815 rpm. This 815 rpm instruction is then sent to the frequency converter for execution.

[0071] Optionally, the method further includes:

[0072] The theoretical recurrence value vector is generated by performing matrix multiplication operations between the polynomial extended design matrix and the regression coefficient vector.

[0073] The pre-constructed polynomial extended design matrix is ​​read from memory. Each row of this matrix represents a polynomial feature combination of a historical sample, with the number of rows representing the total number of samples and the number of columns representing the total number of features. Simultaneously, the regression coefficient vector calculated using the least squares method is read; the dimension of this vector is consistent with the number of columns in the polynomial extended design matrix. Following the rules of linear algebra matrix multiplication, the polynomial extended design matrix and the regression coefficient vector are multiplied. Specifically, the first row vector of the polynomial extended design matrix is ​​taken, and each element is multiplied by the corresponding element of the regression coefficient vector. All products are summed to obtain the theoretical recurrence value of the first sample. This inner product operation is performed sequentially on each row of the matrix, ultimately generating a column vector, which is defined as the theoretical recurrence value vector.

[0074] Calculate the difference vector between the theoretical reproduction value vector and the unit output power value, and define the difference vector as the static fitting deviation vector.

[0075] Extract the unit output power value column from the working condition quality mapping matrix and define it as the actual power vector. Check if the number of rows in this actual power vector is strictly equal to the number of rows in the theoretical reproduction value vector. Perform vector subtraction, using the theoretical reproduction value vector as the minuend and the actual power vector as the subtrahend. That is, read the first value in the theoretical reproduction value vector and subtract the first value in the actual power vector to obtain the first deviation value. Perform this subtraction operation on the values ​​at each position in the sequence. Combine all the calculated differences into a new column vector in their original order, and define this column vector as the static fitting deviation vector. This vector physically represents the distribution of the fitting error of the current regression equation on historical data.

[0076] For example, the multinomial extended design matrix contains 100 historical samples. After feature expansion, each sample contains 10 feature terms, resulting in a matrix dimension of 100 rows and 10 columns. The regression coefficient vector has a dimension of 10 rows and 1 column. After performing matrix multiplication, a 100-row, 1-column theoretical recurrence value vector is generated. The theoretical recurrence value calculated for the first row of samples is 45.5 kWh per ton. The corresponding actual unit output power is 45.0 kWh per ton. Subtraction is performed: 45.5 minus 45.0 equals 0.5. The theoretical recurrence value calculated for the second row of samples is 48.0, and the actual value is 48.2. Subtraction yields -0.2. The final static fit bias vector contains 100 values, including 0.5 and -0.2.

[0077] Optionally, generating the modified regression equation includes:

[0078] The inverter frequency and valve opening command are substituted into the regression equation to perform forward calculation to obtain the theoretical prediction value. The difference between the new unit output power value and the new particle size detection value and the theoretical prediction value is calculated to obtain the prediction residual vector.

[0079] The system reads the actual inverter frequency and valve opening values ​​issued to the equipment in the current cycle and combines these two values ​​into a command vector. This command vector is then used as the independent variable and substituted into the regression equation with a defined control cycle. Following the mathematical structure of the regression equation, the system calculates the sum of the products of each coefficient and the independent variable, outputting a theoretical estimate of the unit output power value and the sum of squared particle size differences. This result is defined as the theoretical prediction value. Simultaneously, real-time unit output power values ​​are collected via sensors, and new particle size detection data is obtained through laboratory testing. Subtraction operations are performed: the power component of the theoretical prediction value is subtracted from the newly collected unit output power value, and the particle size component of the theoretical prediction value is subtracted from the sum of squared particle size differences corresponding to the newly acquired particle size detection data. These two calculated differences are combined into a column vector, defined as the prediction residual vector.

[0080] For example, the instruction vector includes an inverter frequency of 50 Hz and a valve opening of 80 degrees. Substituting these into the regression equation, the theoretical prediction value is: power 45 kWh per ton, and the sum of squares of particle size differences is 0.05. The actual collected new unit output power value is 46 kWh per ton, and the sum of squares of particle size differences calculated from the new particle size detection value is 0.04. Calculating the difference: the power difference is 46 minus 45 equals 1; the particle size difference is 0.04 minus 0.05 equals -0.01. The generated prediction residual vector contains the value 1 and the value -0.01.

[0081] A vector autocorrelation matrix is ​​constructed using the inverter frequency and valve opening command. The inverse matrix of the vector autocorrelation matrix is ​​calculated to generate an orthogonal projection gain matrix. The inverter frequency and valve opening command, the prediction residual vector, and the orthogonal projection gain matrix are multiplied by a chain to generate a coefficient correction matrix.

[0082] The inverter frequency and valve opening commands are extracted to form an input feature vector. A vector self-external product operation is performed on this input feature vector, i.e., the product of the vector and its transpose is calculated, generating a square matrix, which is defined as the vector autocorrelation matrix. The inverse matrix of this vector autocorrelation matrix is ​​calculated using a matrix inversion algorithm, and the result is defined as the orthogonal projection gain matrix. Chained matrix multiplication is performed: first, the orthogonal projection gain matrix is ​​multiplied by the input feature vector to obtain an intermediate projection vector; then, this intermediate projection vector is multiplied by the prediction residual vector, generating a matrix with the same dimensions as the regression coefficient matrix, which is defined as the coefficient correction matrix.

[0083] Extract the current coefficient matrix from the regression equation, perform matrix addition on the current coefficient matrix and the coefficient correction matrix, and reconstruct the regression equation to generate a modified regression equation.

[0084] Read all the regression coefficients currently in use in the regression equation and arrange them into the current coefficient matrix according to the order of the polynomial terms. Read the coefficient correction matrix calculated in the previous step. Perform matrix addition: add the element value at each position in the current coefficient matrix to the corresponding element value in the coefficient correction matrix. Replace the original coefficient values ​​with the new sum. Reconstruct the mathematical expression of the regression equation using the updated coefficient matrix, define the updated regression equation as the corrected regression equation, and use the corrected regression equation for prediction in the next control period.

[0085] Optionally, the step of performing a second calibration and archiving on the modified regression equation includes:

[0086] Calculate the coefficient difference between the modified regression equation and the regression equation, multiply the coefficient difference by the polynomial extended design matrix to obtain the model drift vector, and calculate the difference between the static fit bias vector and the model drift vector to generate a prediction residual sequence.

[0087] Read the new coefficient vector composed of all regression coefficients in the revised regression equation, and simultaneously read the old coefficient vector from the original regression equation. Perform vector subtraction, subtracting the old coefficient vector from the new coefficient vector to obtain the coefficient difference vector. Read the polynomial extended design matrix from memory. Perform matrix multiplication between the polynomial extended design matrix and the coefficient difference vector to calculate the change in predicted values ​​on historical data due to coefficient adjustments; define this change as the model drift vector. Read the static fit bias vector calculated and stored in the previous steps. Perform vector subtraction, subtracting the model drift vector from the static fit bias vector. This operation uses the difference principle to quickly calculate the residuals of historical data relative to the revised model; the result is defined as the predicted residual sequence.

[0088] For example, a term in the old coefficient vector has a value of 2.0, and the corresponding term in the corrected new coefficient vector has a value of 2.1. The calculated coefficient difference is 0.1. The eigenvalue of a row in the polynomial extended design matrix is ​​10. The calculated model drift is 10 multiplied by 0.1, which equals 1.0. The static fit bias corresponding to this row of data is 1.5. Subtraction is performed: 1.5 minus 1.0 equals 0.5. This value of 0.5 is the prediction residual of this historical sample under the new model.

[0089] The absolute quantile of the predicted residual sequence is extracted as the cleaning boundary. Row vectors with residual moduli greater than the cleaning boundary are removed from the working condition quality mapping matrix. The remaining matrix data is defined as a reliable sample matrix.

[0090] For each value in the predicted residual sequence, perform an absolute value operation and sort the resulting sequence in ascending order. Based on a preset confidence level, such as 95%, locate the value at that level and extract it as the upper quantile. This upper quantile is directly defined as the cleaning boundary. Iterate through each row of data in the working condition quality mapping matrix, indexing the corresponding absolute value of the predicted residual. Perform a logical judgment: if the absolute value of the predicted residual for a row is greater than the cleaning boundary, the row is determined to be an abnormal sample subject to interference, and a deletion operation is performed, removing the row from the working condition quality mapping matrix. After completing the traversal and removal operations, all remaining data rows in the matrix are combined into a new matrix, defined as the reliable sample matrix.

[0091] The new unit output power value and the new particle size detection value are appended to the trusted sample matrix. The inverse matrix of the trusted sample matrix is ​​solved to generate regression coefficients after secondary calibration. The parameters of the modified regression equation are then covered and archived.

[0092] The newly collected unit output power value, new particle size detection value, and corresponding operating condition characteristics are combined into a new row of sample data. This new sample data is appended to the end of the reliable sample matrix. Based on the updated reliable sample matrix, the polynomial extended design matrix and dependent variable matrix are reconstructed. The new polynomial extended design matrix is ​​transposed and multiplied by itself to construct a new normal equation matrix. The inverse matrix of this normal equation matrix is ​​solved using Gaussian elimination. This inverse matrix is ​​multiplied sequentially by the transposed matrix and the dependent variable matrix to calculate a new regression coefficient vector. This regression coefficient vector is defined as the regression coefficients after secondary calibration. The existing coefficients in the corrected regression equation are replaced with these secondary calibration regression coefficients to complete parameter coverage. The final determined equation parameters are saved to the historical database for archiving, for use in the next control cycle, such as... Figure 3 As shown, using the observed sequence as a reference, the "predicted values ​​of the corrected regression equation" and the "predicted values ​​after secondary calibration" are simultaneously projected onto the Taylor plot coordinate system. The radial distance from the point to the origin represents the standard deviation of the sequence, corresponding to the scales on the lower and left coordinate axes. The angle between the point and the positive direction of the horizontal axis corresponds to the correlation coefficient. The scale on the outer arc ranges from 0 to 1, with closer values ​​to 1 indicating higher trend consistency. The dotted-dash arc centered on the observation point represents the centered root mean square error contour lines, with smaller values ​​indicating smaller overall error. Comparing the positions of the two predicted points relative to the observation point reveals that if the predicted point after secondary calibration is closer to the observation point and radially closer to the observed standard deviation, falling on a smaller contour line, it indicates that secondary calibration, without changing the evaluation benchmark, simultaneously improves the predicted sequence in three aspects: "reproduction of fluctuation amplitude," "trend consistency," and "overall bias after mean removal," demonstrating the effectiveness of secondary calibration in improving model prediction performance.

[0093] For example, the original reliable sample matrix had 98 rows of data. After adding one new row, it became 99 rows. Using these 99 rows of data, a matrix was constructed and the inverse matrix operation was performed to obtain a new set of regression coefficients, where the coefficient of the first-order term was calculated to be 2.15. This coefficient in the original modified regression equation was 2.1. 2.15 was used to overwrite 2.1, and 2.15 was stored as the final determined model parameter in the hard disk database.

[0094] Optionally, the archive includes:

[0095] The number of rows in the reliable sample matrix is ​​counted to determine the current sample size, and the window length of the sliding window is extracted to determine the upper limit of the capacity.

[0096] Access the trusted sample matrix in memory, traverse the matrix using a row counting algorithm, count the total number of data rows contained within the matrix, and define this statistical result as the current sample size. Read the sliding window parameter value set in the configuration file, which represents the maximum time span or number of samples allowed when the model is trained or calibrated on historical data, and define this value as the capacity limit.

[0097] The current sample size minus the capacity limit is calculated as the overflow quantity. If the overflow quantity is greater than zero, then starting from the starting row index of the trusted sample matrix, row vectors equal to the overflow quantity are removed sequentially.

[0098] Perform a subtraction operation, subtracting the capacity limit from the current sample size, to obtain a difference value, which is defined as the overflow quantity. Perform a logical check on this overflow quantity: check if the value is greater than zero. If the result is yes, it indicates that the currently stored data volume has exceeded the allowed range. Locate the starting row index of the trusted sample matrix. Starting from this starting row index, select row vectors consecutively in the direction of increasing rows, with the number equal to the overflow quantity. Perform a physical deletion operation on these selected row vectors, releasing the storage space they occupy, thereby achieving first-in-first-out (FIFO) management of the sample database. If the result is no, no deletion operation is performed.

[0099] For example, the sliding window length is set to 5000 rows. After data appending and cleaning, the number of rows in the statistically reliable sample matrix is ​​5005. The current sample size of 5005 is subtracted from the capacity limit of 5000, resulting in an overflow of 5. Since 5 is greater than zero, rows 1 to 5 of the matrix are located. A deletion command is executed to permanently remove these 5 rows of data. The remaining 5000 rows of data in the matrix are retained for the next round of calculation.

[0100] Based on the same inventive concept, the present invention also provides a multi-objective collaborative control system for powder grinding, such as... Figure 4 As shown, the system includes:

[0101] The working condition quality mapping construction module is used to extract the time-series data of the classifier load and the time-series data of the mill main motor power in the distributed control system, perform signal noise reduction processing based on variational mode decomposition, separate the internal model component, obtain the offline particle size detection data of the automated laboratory, and perform convolution mapping on the internal model component and the offline particle size detection data using the residence time distribution of the material in the grinding system to generate the working condition quality mapping matrix.

[0102] The regression prediction and undominated screening module is used to perform least squares polynomial regression on the working condition quality mapping matrix, with the classifier speed, circulating fan speed and total feed amount as independent variables, and the unit output power value and the sum of squares of particle size difference as dependent variables to obtain the regression equation. Within the domain of the independent variables, a test vector set is generated and substituted into the regression equation to calculate the predicted value. Numerical dominance comparison is performed to remove dominated vectors and extract undominated vectors to construct candidate control parameters.

[0103] The candidate vector sorting and instruction generation module is used to calculate the instantaneous frequency variance of the internal model component, construct dynamic weights to perform weighted summation on the unit output power value and the sum of squares of the particle size difference, select the candidate control vector with the smallest weighted sum, calculate the differential increment through the candidate control vector and the current set value and perform amplitude truncation, and generate frequency converter frequency and valve opening instructions for issuance and execution.

[0104] The prediction residual calculation and regression update module is used to execute the frequency converter frequency and valve opening command to collect new unit output power value and new particle size detection value, calculate the prediction residual vector of the new unit output power value and the new particle size detection value relative to the regression equation, and perform least squares iterative update for the regression equation to generate a modified regression equation.

[0105] The residual threshold calculation and outlier removal module is used to extract the predicted residual sequence of historical data using the modified regression equation, calculate the upper quantile boundary of the predicted residual sequence within the sliding window, traverse the working condition quality mapping matrix and remove outlier data rows whose residual modulus exceeds the upper quantile boundary, append the new unit output power value and the new granularity detection value to the working condition quality mapping matrix, perform least squares regression to perform secondary calibration of the modified regression equation and archive it.

[0106] It should be noted that the functional division and information interaction between the various modules described above are logical, but in terms of physical implementation, they can be integrated on the same software platform or deployed in a distributed manner. The connections between them represent data flow and control flow, aiming to collaboratively achieve the objectives of this invention. The above descriptions are merely exemplary embodiments of this invention and should not be construed as limiting the scope of protection of this invention.

Claims

1. A multi-objective collaborative control method for powder grinding, characterized in that, The method includes: Extract the time-series data of the classifier load and the time-series data of the mill main motor power from the distributed control system, perform signal noise reduction processing based on variational mode decomposition, separate the internal model component, obtain the offline particle size detection data of the automated laboratory, and perform convolution mapping on the internal model component and the offline particle size detection data using the residence time distribution of the material in the grinding system to generate the working condition quality mapping matrix; The least squares polynomial regression is performed on the working condition quality mapping matrix. The classifier speed, circulating fan speed and total feed amount are used as independent variables, and the unit output power value and the sum of squares of particle size difference are used as dependent variables to obtain the regression equation. Test vector sets are generated by sampling within the domain of the independent variables and substituted into the regression equation to calculate the predicted value. Numerical dominance comparison is performed to remove dominated vectors and extract non-dominated vectors to construct candidate control parameters. Calculate the instantaneous frequency variance of the internal model component, construct dynamic weights to perform a weighted summation of the unit output power value and the sum of squares of the particle size difference, select the candidate control vector with the smallest weighted sum, calculate the differential increment through the candidate control vector and the current set value and perform amplitude truncation, generate inverter frequency and valve opening commands and issue them for execution; The inverter frequency and valve opening commands are executed to collect new unit output power value and new particle size detection value. The prediction residual vector of the new unit output power value and the new particle size detection value relative to the regression equation is calculated. The least squares iterative update is performed on the regression equation to generate a modified regression equation. The predicted residual sequence of historical data is extracted using the modified regression equation. The upper quantile boundary of the predicted residual sequence within the sliding window is calculated. The working condition quality mapping matrix is ​​traversed and abnormal data rows with residual modulus exceeding the upper quantile boundary are removed. The new unit output power value and the new granularity detection value are appended to the working condition quality mapping matrix. Least square regression is performed to perform secondary calibration of the modified regression equation and archive it.

2. The multi-objective collaborative control method for powder grinding according to claim 1, characterized in that, The generated working condition quality mapping matrix includes: Variational mode decomposition is performed on the time series data of the classifier load and the time series data of the mill main motor power to output a set of modal components. The cross-correlation coefficient between each component in the set of modal components and the original time series data is calculated, and the component with the first cross-correlation coefficient is locked as the inner mode component. The sampling time is extracted from the offline granularity detection data. The cross-correlation function between the internal model component and the offline granularity detection data is calculated, and the time shift of the maximum value is used to generate the average lag time. The dwell time distribution function is constructed with the average lag time as the expectation and discretization is performed to obtain the time lag weight sequence. Using the sampling time minus the average lag time as an index, historical time window data with the same length as the time lag weight sequence is extracted from the internal model component. A weighted summation operation is performed on the historical time window data and the time lag weight sequence to obtain a weighted operating condition feature value. The weighted operating condition feature value is then concatenated with the offline granularity detection data to generate an operating condition quality mapping matrix.

3. The multi-objective collaborative control method for powder grinding according to claim 1, characterized in that, The extraction of undominated vectors to construct candidate control parameters includes: Based on the classifier speed, circulating fan speed and total feed amount extracted from the working condition quality mapping matrix, a polynomial extended design matrix is ​​constructed, and matrix operations are performed in combination with the unit output power value and the sum of squared particle size differences to calculate the regression coefficient vector to establish the regression equation. Extract the domain boundary values ​​of the independent variable, perform equal-step grid scanning within the domain boundary values ​​to generate a test vector set, and substitute the test vector set into the regression equation to generate a set of predicted values. The predicted value set is traversed and pairwise numerical comparisons are performed. If the unit output power value and the sum of squared particle size differences of the first vector are both less than those of the second vector, then the second vector is marked as a dominated vector and removed, and the unmarked vectors are retained as candidate control parameters.

4. The multi-objective collaborative control method for powder grinding according to claim 3, characterized in that, The vector for calculating regression coefficients includes: Calculate the squared vector and the product vector of the data column, and concatenate the data column, the squared vector and the product vector into column vectors to generate a polynomial extended design matrix; The inverse matrix is ​​obtained by performing operations on the polynomial extended design matrix, and a chain matrix multiplication is performed on the inverse matrix, the transpose of the polynomial extended design matrix, the unit output power value, and the sum of squared granularity differences to generate a regression coefficient vector.

5. The multi-objective collaborative control method for powder grinding according to claim 1, characterized in that, The generation and execution of the inverter frequency and valve opening command include: The number of zero-crossing points of the internal model component within the sampling time is counted. The sampling time is divided by the number of zero-crossing points to obtain the average period. The variance of the amplitude of each point of the internal model component relative to the average period is calculated as the instantaneous frequency variance. The instantaneous frequency variance is mapped to a normalized value using the hyperbolic tangent function as a dynamic weight for the sum of squares of the granularity difference. The average value of unit output power and the average value of the sum of squares of particle size difference are calculated using the candidate control parameters. The average value of unit output power and the average value of the sum of squares of particle size difference are divided by the corresponding average value of unit output power and the average value of the sum of squares of particle size difference to generate a dimensionless ratio. The dimensionless ratio is weighted and summed using the dynamic weights. The vector with the smallest summation result is selected as the preferred vector. The standard deviation of the classifier speed is calculated using the working condition quality mapping matrix as a safety fluctuation threshold. The difference between the preferred vector and the current set value is calculated. If the absolute value of the difference is greater than the safety fluctuation threshold, the frequency converter frequency and valve opening command are generated.

6. The multi-objective collaborative control method for powder grinding according to claim 4, characterized in that, The method further includes: The theoretical recurrence value vector is generated by performing matrix multiplication operations between the polynomial extended design matrix and the regression coefficient vector. Calculate the difference vector between the theoretical reproduction value vector and the unit output power value, and define the difference vector as the static fitting deviation vector.

7. The multi-objective collaborative control method for powder grinding according to claim 1, characterized in that, The generation of the modified regression equation includes: The inverter frequency and valve opening command are substituted into the regression equation to perform forward calculation to obtain the theoretical prediction value. The difference between the new unit output power value and the new particle size detection value and the theoretical prediction value is calculated to obtain the prediction residual vector. A vector autocorrelation matrix is ​​constructed using the inverter frequency and valve opening command. The inverse matrix of the vector autocorrelation matrix is ​​calculated to generate an orthogonal projection gain matrix. The inverter frequency and valve opening command, the prediction residual vector, and the orthogonal projection gain matrix are multiplied by a chain to generate a coefficient correction matrix. Extract the current coefficient matrix from the regression equation, perform matrix addition on the current coefficient matrix and the coefficient correction matrix, and reconstruct the regression equation to generate a modified regression equation.

8. The multi-objective collaborative control method for powder grinding according to claim 6, characterized in that, The step of performing secondary calibration and archiving of the modified regression equation includes: Calculate the coefficient difference between the modified regression equation and the regression equation, multiply the coefficient difference by the polynomial extended design matrix to obtain the model drift vector, and calculate the difference between the static fit bias vector and the model drift vector to generate a prediction residual sequence. The absolute quantile of the predicted residual sequence is extracted as the cleaning boundary. Row vectors with residual moduli greater than the cleaning boundary are removed from the working condition quality mapping matrix. The remaining matrix data is defined as a reliable sample matrix. The new unit output power value and the new particle size detection value are appended to the trusted sample matrix. The inverse matrix of the trusted sample matrix is ​​solved to generate regression coefficients after secondary calibration. The parameters of the modified regression equation are then covered and archived.

9. The multi-objective collaborative control method for powder grinding according to claim 8, characterized in that, The archive includes: The number of rows in the reliable sample matrix is ​​counted to determine the current sample size, and the window length of the sliding window is extracted to determine the upper limit of the capacity. The current sample size minus the capacity limit is calculated as the overflow quantity. If the overflow quantity is greater than zero, then starting from the starting row index of the trusted sample matrix, row vectors equal to the overflow quantity are removed sequentially.

10. A multi-objective collaborative control system for powder grinding, applied to the multi-objective collaborative control method for powder grinding as described in any one of claims 1-9, characterized in that, The system includes: The working condition quality mapping construction module is used to extract the time-series data of the classifier load and the time-series data of the mill main motor power in the distributed control system, perform signal noise reduction processing based on variational mode decomposition, separate the internal model component, obtain the offline particle size detection data of the automated laboratory, and perform convolution mapping on the internal model component and the offline particle size detection data using the residence time distribution of the material in the grinding system to generate the working condition quality mapping matrix. The regression prediction and undominated screening module is used to perform least squares polynomial regression on the working condition quality mapping matrix, with the classifier speed, circulating fan speed and total feed amount as independent variables, and the unit output power value and the sum of squares of particle size difference as dependent variables to obtain the regression equation. Within the domain of the independent variables, a test vector set is generated and substituted into the regression equation to calculate the predicted value. Numerical dominance comparison is performed to remove dominated vectors and extract undominated vectors to construct candidate control parameters. The candidate vector sorting and instruction generation module is used to calculate the instantaneous frequency variance of the internal model component, construct dynamic weights to perform weighted summation on the unit output power value and the sum of squares of the particle size difference, select the candidate control vector with the smallest weighted sum, calculate the differential increment through the candidate control vector and the current set value and perform amplitude truncation, and generate frequency converter frequency and valve opening instructions for issuance and execution. The prediction residual calculation and regression update module is used to execute the frequency converter frequency and valve opening command to collect new unit output power value and new particle size detection value, calculate the prediction residual vector of the new unit output power value and the new particle size detection value relative to the regression equation, and perform least squares iterative update for the regression equation to generate a modified regression equation. The residual threshold calculation and outlier removal module is used to extract the predicted residual sequence of historical data using the modified regression equation, calculate the upper quantile boundary of the predicted residual sequence within the sliding window, traverse the working condition quality mapping matrix and remove outlier data rows whose residual modulus exceeds the upper quantile boundary, append the new unit output power value and the new granularity detection value to the working condition quality mapping matrix, perform least squares regression to perform secondary calibration of the modified regression equation and archive it.

Images