Pig feed production quality tracing method based on cloud computing

Abstract

The invention provides a pig feed production quality tracing method based on cloud computing. Comprising the steps of whole-process data acquisition and structuring of multi-source heterogeneous process parameters and quality detection indexes, sliding window segmentation normalization and time alignment processing, causal relationship identification based on Granger causal test and information transfer entropy, and dynamic optimization of a causal atlas by an improved time-aware structure equation model. According to the method, the accuracy and the analysis efficiency of feed production abnormity traceability are remarkably improved, key influence procedures can be conveniently and rapidly positioned, and the accuracy and the analysis efficiency of feed production abnormity traceability are improved, so that the accuracy and the analysis efficiency of the feed production abnormity traceability are improved, and the accuracy and the analysis efficiency of the feed production abnormity traceability are improved. And quality control and process optimization are facilitated.

Description

Technical Field

[0001] This invention relates to the field of quality traceability and intelligent diagnostic technology in feed production processes, and in particular to a cloud computing-based method for traceability of pig feed production quality. Background Technology

[0002] With the continuous improvement of livestock farming scale and automation, quality management and traceability of the entire feed production process has become a core technology for ensuring food safety and product quality in the farming process. Currently, the industry generally adopts methods based on data collection and statistical analysis, covering processes such as raw material pretreatment, mixing and formulation, high-temperature pelleting, and finished product testing. Mainstream technical solutions include using IoT sensors and automated testing equipment to collect process parameters (such as moisture, temperature, and mixing uniformity) and finished product quality indicators in real time. Data integration and standardization are achieved through cloud platforms or local area network systems. Further, correlation analysis, statistical regression modeling, simple rule diagnosis, or graph-structured traceability methods are used for anomaly identification and source analysis. Some companies or research teams have attempted to introduce modeling techniques such as Bayesian networks and graph neural networks to support correlation analysis and propagation path mining between multi-level processes; however, their core logic still focuses on variable correlation, lacking in-depth expression of causal relationships and interpretable inferences. Existing feed production quality traceability systems typically emphasize the complete collection of information flow and data-driven anomaly detection. For example, they employ methods such as quality information coding, process identification, and batch traceability to track the production process, facilitating the later location of problematic batches and the allocation of quality management responsibilities. However, in the specific anomaly propagation path identification and modeling stage, they largely rely on correlation analysis and one-way causal attribution, failing to reveal the true impact path of abnormal parameters across processes and levels, and struggling to effectively distinguish between direct causal drivers and indirect correlations between variables under complex operating conditions. Some solutions attempt to model the transition probabilities between process variables using Bayesian networks, but face technical bottlenecks such as high model structure rigidity, difficulty in real-time dynamic adaptation of parameter updates, and insufficient topological inference capabilities for complex multi-process data. Furthermore, anomaly traceability schemes based on graph neural networks suffer from a "black box effect"—while they can output propagation paths, the model's interpretability is poor, making it difficult to effectively correlate specific causal nodes with the actual state of the production process, thus weakening the practical diagnostic value of anomaly traceability. Summary of the Invention

[0003] In order to solve the above-mentioned technical problems, the present invention provides a cloud computing-based method for tracing the quality of pig feed production.

[0004] The technical solution of this invention is implemented as follows: A cloud computing-based method for tracing the quality of pig feed production, comprising: S1: Collect process parameters and quality test indicators for each stage of the feed production process, including raw material pretreatment, mixing and formulation, and high-temperature pelleting. The process parameters include time-series data of moisture content, temperature gradient curve and dispersion value of mixing uniformity. The quality test indicators include the deviation rate of finished product protein content and the dispersion coefficient of pellet hardness. S2: Perform adaptive filtering and time alignment processing on the collected multi-source heterogeneous data, and use a sliding window strategy to segment and normalize the process parameters to generate a standardized time series dataset containing process feature labels. The window length is dynamically adjusted according to the feature time constant of each stage. S3: Based on the joint criterion of Granger causality test and information transfer entropy, the causal strength of variable pairs in the standardized time series dataset is calculated, and an initial causal topology is constructed, where nodes represent process quality characteristics, edge weights represent cross-process anomaly propagation strength, and directionality reflects the propagation time sequence relationship. S4: An improved time-aware structural equation model (T-SEM) is used to dynamically optimize the initial causal topology. The causal edge weights are updated in real time using a variational Bayesian inference method to generate a causal graph sequence that changes with production conditions. The model parameters are iteratively updated based on the latest quality inspection data within the sliding window. S5: When the finished product quality deviation rate is detected to exceed the preset threshold, the counterfactual reasoning mechanism is activated. A virtual intervention matrix is ​​constructed based on the causal graph sequence. By calculating the hypothetical intervention effect of 'if no parameter shift occurred in a specific upstream process', the attribution contribution of each process to the final anomaly is quantified. S6: Input the causal path with an attribution contribution exceeding the preset significance level into the expert rule knowledge base, perform semantic annotation and path compression operations, and generate a visual tracing report containing the abnormal propagation time sequence characteristics, key nodes and propagation intensity. The annotation basis includes process specification constraints and historical failure mode library.

[0005] The present invention provides a cloud computing-based method for tracing the quality of pig feed production, which has the following beneficial effects: (1) This invention significantly improves the accuracy and mechanism transparency of anomaly source identification by constructing a modeling framework based on dynamic causal graphs and counterfactual reasoning. It adopts an improved time-aware structural equation model (T-SEM) and combines a sliding window strategy to segment and model historical working conditions. This can adaptively capture the dynamic causal dependencies between key parameters at different production stages and generate a causal graph sequence that evolves over time, thus overcoming the limitations of static modeling in adapting to non-stationary processes. At the same time, by introducing the joint criterion of Granger causality test and information transfer entropy to construct the initial topology, the robustness of causal direction judgment is enhanced, avoiding the problems of false or missed connections caused by a single criterion, and providing reliable basic topological support for subsequent accurate attribution. (2) This invention innovatively integrates counterfactual reasoning mechanism to construct an attribution calculation path oriented towards "what-if" hypothesis intervention, effectively overcoming the technical bottleneck of existing black box models lacking decision-making basis output. When a batch of finished products has quality deviations, the counterfactual scenario of "if a certain upstream process did not have parameter shifts" is simulated through virtual intervention, thereby assessing the degree of impact of the intervention on the final output result, and calculating the attribution score of each potential node accordingly. This mechanism not only realizes the ability to trace back from the result to the source, but also quantifies the contribution weight of different process links to abnormal events, making the diagnostic conclusions have clear logical support and engineering understandability. In addition, the variational Bayesian inference method is combined to update the causal edge weights online, enabling the model to respond to causal relationship drift caused by factors such as equipment aging and raw material fluctuations in actual production, significantly enhancing the dynamic adaptability and long-term operational stability of the method. The final generated visual traceability report marks high contribution paths by integrating semantic rules in the expert knowledge base, further improving the professionalism and practicality of the diagnostic results, making it easier for technicians to quickly locate key control points and take corrective measures; (3) This invention abandons the high dependence on large-scale training data and computing resources, and constructs a lightweight, efficient, and interpretable modeling paradigm that does not require complex parameter tuning, thus possessing stronger industrial feasibility. The entire modeling process does not rely on an end-to-end black-box learning architecture, avoiding the problems of trust loss and debugging difficulties caused by uncontrollable models, and significantly reducing the deployment threshold while ensuring accuracy. The application of adaptive filtering technology ensures that time-series data from different sensor sources are accurately aligned on the time axis, solving the risk of causal misjudgment caused by asynchronous data acquisition across processes; while the modular processing flow design makes the system have good scalability and can flexibly adapt to changes in production processes of different formulas and production lines. This method is not only applicable to the whole-process quality monitoring of pig feed production, but can also be extended to quality traceability scenarios in other process manufacturing industries, with broad industry application prospects and technology transfer value. Attached Figure Description

[0006] Figure 1 This is a flowchart of a cloud computing-based pig feed production quality traceability method according to the present invention; Figure 2 This is a sub-flowchart of a cloud computing-based pig feed production quality traceability method according to the present invention; Figure 3 This is another sub-flowchart of a cloud computing-based pig feed production quality traceability method according to the present invention. Detailed Implementation

[0007] Embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention.

[0008] The following disclosure provides many different embodiments or examples for implementing different structures of the invention. To simplify the disclosure, specific examples of components and arrangements are described below. Of course, these are merely examples and are not intended to limit the invention. Furthermore, reference numerals and / or letters may be repeated in different examples; such repetition is for simplification and clarity and does not in itself indicate a relationship between the various embodiments and / or arrangements discussed.

[0009] like Figure 1 As shown, this invention provides a cloud computing-based method for tracing the quality of pig feed production, specifically including: S1: Collect process parameters and quality test indicators for each stage of the feed production process, including raw material pretreatment, mixing and formulation, and high-temperature pelleting. The process parameters include time-series data of moisture content, temperature gradient curve and dispersion value of mixing uniformity. The quality test indicators include the deviation rate of finished product protein content and the dispersion coefficient of pellet hardness. S2: Perform adaptive filtering and time alignment processing on the collected multi-source heterogeneous data, and use a sliding window strategy to segment and normalize the process parameters to generate a standardized time series dataset containing process feature labels. The window length is dynamically adjusted according to the feature time constant of each stage. S3: Based on the joint criterion of Granger causality test and information transfer entropy, the causal strength of variable pairs in the standardized time series dataset is calculated, and an initial causal topology is constructed, where nodes represent process quality characteristics, edge weights represent cross-process anomaly propagation strength, and directionality reflects the propagation time sequence relationship. S4: An improved time-aware structural equation model (T-SEM) is used to dynamically optimize the initial causal topology. The causal edge weights are updated in real time using a variational Bayesian inference method to generate a causal graph sequence that changes with production conditions. The model parameters are iteratively updated based on the latest quality inspection data within the sliding window. S5: When the finished product quality deviation rate is detected to exceed the preset threshold, the counterfactual reasoning mechanism is activated. A virtual intervention matrix is ​​constructed based on the causal graph sequence. By calculating the hypothetical intervention effect of 'if no parameter shift occurred in a specific upstream process', the attribution contribution of each process to the final anomaly is quantified. S6: Input the causal path with an attribution contribution exceeding the preset significance level into the expert rule knowledge base, perform semantic annotation and path compression operations, and generate a visual tracing report containing the abnormal propagation time sequence characteristics, key nodes and propagation intensity. The annotation basis includes process specification constraints and historical failure mode library.

[0010] Step S1: Collect process parameters and quality inspection indicators for each stage of the feed production process, including raw material pretreatment, mixing and formulation, and high-temperature pelleting. The process parameters include time-series data on moisture content, temperature gradient curves, and dispersion values ​​of mixing uniformity. The quality inspection indicators include the deviation rate of finished product protein content and the dispersion coefficient of pellet hardness. Specifically, this includes: S1.1: Based on the moisture and temperature sensors in the raw material pretreatment process, collect time-series data of raw material moisture content and temperature gradient curves to reflect the stability of the initial state of the raw material. Based on the output signals of moisture and temperature sensors deployed in the raw material pretreatment stage, an analog-to-digital conversion acquisition module (parameters: sampling frequency 2Hz, resolution 16bit) is used to achieve accurate digital sampling of the moisture content and temperature of the raw materials. Furthermore, by using a signal integrity verification algorithm (parameter: Cyclic Redundancy Check CRC-16), the validity of the sampled data packets is determined, and sensor pulse sequences that fail the verification are removed, resulting in a continuous and error-free original sensor data stream. Furthermore, a moving average filtering method (parameter: window length of 5 points) is used to perform noise suppression processing on the moisture content time series data and obtain a smoothed moisture trend curve to reduce the impact of random fluctuations on the raw material stability assessment. Furthermore, a polynomial fitting algorithm (order: second) is used to fit the instantaneous sampled values ​​of the temperature sensor to generate a temperature gradient curve, and the rate of change of the temperature gradient is calculated to reflect the uniformity of heat transfer during the raw material pretreatment process. Furthermore, the process stability criterion function is used. The stability of the moisture content sequence was quantitatively assessed, among which... The standard deviation of moisture content, This is the average moisture content; the smaller the value, the more stable the initial moisture state of the raw material. Through the above-mentioned sensor data acquisition, filtering, fitting and stability calculation processing methods, the simulation signal results of the previous step are transformed into structured moisture content time series data and temperature gradient curves, so as to realize the quantitative characterization of the initial state stability of the raw materials. For example, a feed mill configures a capacitive moisture sensor (model HC-WMS100, sampling frequency 2Hz, measurement range 0%~50%, accuracy ±0.2%) and a silicon temperature sensor (model KTY84-130, sampling frequency 2Hz, measurement range 0℃~200℃, accuracy ±0.5℃) in the raw material pretreatment stage. The sensor output signals are sampled by a 24-bit ADC module and subjected to CRC-16 verification, eliminating all invalid frame signals appearing within 2 seconds. A moving average filter with a window length of 5 points is applied to the moisture data, smoothing it and reducing the maximum fluctuation amplitude to 0.3%. A second-order polynomial fitting is performed on the temperature data to obtain the temperature gradient curve, and the gradient change rate is calculated. ℃ / s indicates high stability of heat transfer. Combined with the standard deviation of moisture content... with the mean Calculate the process stability criterion value =0.004, which greatly improves the accuracy of the initial stability assessment of raw materials and outputs a structured moisture time series and temperature gradient curve dataset for subsequent process correlation analysis of mixing and granulation. S1.2: Based on the weighing system and mixing uniformity detection device in the mixing and batching process, collect discrete value data of mixing uniformity to quantify the uniformity control level of the batching process; Based on the real-time output signals of the weighing system and the mixing uniformity detection device in the mixing and batching process, a parallel sampling synchronous interface (parameters: sampling frequency 1Hz, interface delay ≤5ms) is used to realize the parallel acquisition and time calibration of weighing data and raw mixing uniformity signals. Furthermore, through a discretization and quantization algorithm (parameters: the output of the detection device is an analog signal, the quantization bit is 12 bits, and the sampling period is 1 second), the discrete value conversion of the continuous signal of mixing uniformity is realized, and a uniformity sequence bound to the batch identifier of the ingredients is obtained. Furthermore, a sliding window statistical method is adopted (parameters: window length is 1 / 5 of the batching cycle, step size is half of the window length) to calculate the mean and variance of the discrete values ​​of mixing uniformity within each time window, and generate a set of statistical indicators reflecting the level of local uniformity. Furthermore, by using a weighted mean square error normalization algorithm (with weights set based on the reliability scores of each weighing channel), the uniformity data of different weighing channels are fused, and a global mixed uniformity discrete value is generated as a quantitative value for the uniformity control level of the batching process. Furthermore, an anomaly detection method (joint criterion: the comprehensive discrete value of uniformity exceeds the lower limit of uniformity in the process specification and the variance fluctuation exceeds the set threshold) is adopted to identify abnormal events in the mixing uniformity and provide anomaly labeling data for subsequent quality causal modeling; Through the above-mentioned algorithm chain processing method, the collected signal in the previous step is transformed into mixed uniformity feature data containing batch identifier, timestamp, comprehensive discrete value and anomaly mark, so as to realize the precise quantification of the uniformity control level of the batching process. For example, in a pig feed production batch, the weighing system in the mixing and batching stage is configured as a four-channel electronic scale, with a single channel capacity of 50kg and a sampling frequency of 1Hz. The uniformity detection device uses an image recognition particle distribution analyzer, outputting an analog voltage signal with a quantization bit depth of 12 bits and a sampling period of 1s. During data acquisition, a parallel sampling synchronization interface enables the synchronous acquisition of four-channel weighing data and uniformity analog signals, with a delay controlled within 5ms. After processing by a discretization and quantization algorithm, a uniformity discrete value is obtained per second, with a value range between 0 and 1024. The sliding window statistical method sets the window length to 1 / 5 of the batching cycle (e.g., 200s), i.e., 40s, and the step size to 20s. The mean and variance of the uniformity discrete value for each window are calculated to obtain a local uniformity index, such as a mean value between 512 and 530 and a variance between 6 and 8. The weighted mean squared error normalization algorithm integrates the uniformity data of each weighing channel based on the credibility score of the weighing channel (e.g., channel 1 weight 0.3, channel 2 weight 0.2, channel 3 weight 0.3, channel 4 weight 0.2) to generate a comprehensive discrete value, which is concentrated between 520 and 525. When the comprehensive discrete value is lower than the lower limit of the process specification (e.g., 515) and the variance fluctuation exceeds the threshold (e.g., 10), the anomaly feature detection method automatically marks the time segment as a uniformity anomaly event and outputs the batch identifier, anomaly timestamp, anomaly comprehensive value and corresponding variance, providing highly reliable uniformity anomaly source data for cloud-based causal modeling. S1.3: Based on the temperature controller and pressure sensor in the high-temperature granulation process, collect the temperature change curve and pressure fluctuation data during the granulation process to reflect the stability of the granulation conditions; The data collected consisted of the output signal of the temperature controller and the sampling data of the pressure sensor in the high-temperature granulation process. The goal was to obtain the temperature change curve and pressure fluctuation characteristics that reflect the stability of the granulation conditions. A high-precision thermocouple temperature sensor (parameters: sampling frequency 1Hz, range 0~200℃) is used to obtain the instantaneous temperature values ​​at different locations in the granulation chamber, thereby realizing the time-series recording of the temperature field during the production process. Furthermore, by using a multi-point temperature data fusion algorithm (parameters: weighting coefficients are determined based on the sensor location and the heat flow distribution characteristics of the granulation chamber), the spatial mean curve of the temperature field is generated, and the temperature gradient curve data results over time are obtained. Furthermore, a strain gauge pressure sensor (parameters: range 0~5MPa, sampling frequency 10Hz) was used to collect the working pressure between the pelletizing roller and the die to obtain the raw time-series data of pressure fluctuations; Furthermore, the pressure time series data is filtered by a bandpass filter method (parameters: frequency range 0.1~5Hz, matching the vibration characteristics of the granulation machinery) to remove high-frequency noise and low-frequency drift, and generate a purified pressure fluctuation curve. Furthermore, the instantaneous rate of change of the temperature curve and the peak-valley difference of the pressure fluctuation curve are obtained through statistical calculation methods, and the two are combined into a joint index to characterize the stability of the granulation process. By simultaneously acquiring and processing temperature gradient curves and pressure fluctuation curves, the raw signals from the previous step are transformed into structural chemical condition stability data, providing key inputs for subsequent causal modeling of quality anomalies. For example, a feed mill configures three thermocouple temperature sensors in its high-temperature pelleting process, located at the feed inlet, the middle of the pelleting chamber, and the discharge outlet, respectively. The sampling frequency is 1Hz, and the measurement range is 0~200℃. The weighting coefficients of the temperature fusion algorithm are set to 0.2, 0.5, and 0.3, respectively, to obtain a spatially averaged temperature curve and calculate the temperature gradient. The gradient calculation formula is as follows:

[0011] in, The average temperature. Sampling time, The sampling interval is [value]. The pressure acquisition device has a range of 0~5MPa and a sampling frequency of 10Hz. After bandpass filtering (0.1~5Hz), the peak-to-trough difference in pressure fluctuation is calculated using the following formula:

[0012] in, This represents the maximum pressure value within the time window. This represents the minimum pressure value within the time window. Under stable production conditions during granulation, the temperature gradient curve fluctuation is less than 5℃ / s, and the pressure fluctuation ΔP remains within 0.8MPa. The output stability indicators are significantly improved, effectively supporting the accurate identification of subsequent abnormal propagation paths. S1.4: Based on the finished product quality testing equipment, collect the deviation rate of finished product protein content and the dispersion coefficient of particle hardness to characterize the quality consistency and physical properties of the final feed product; The input conditions for the finished product quality inspection process include the protein content detection value sequence obtained by the infrared protein analyzer and the hardness value sequence corresponding to the crushing force of the particles collected by the particle hardness tester. A content deviation rate calculation method based on standard deviation (parameters: detection value sequence, target content value) is adopted to quantify the dispersion of protein content in each batch of finished products and obtain time-series data of content deviation rate reflecting the degree of deviation from the target value; Furthermore, by using the dispersion coefficient calculation method for particle hardness (parameters: hardness value sequence, batch average hardness), the uniformity of physical strength of finished particles can be evaluated, and the corresponding dispersion coefficient index can be generated. Furthermore, an outlier removal algorithm (parameters: content deviation rate time series data, hardness dispersion coefficient data) is applied to optimize the stability of the test indicators and retain test data that represent the actual working conditions. By combining the protein content deviation rate and particle hardness dispersion coefficient with the batch number and detection timestamp, the data is transformed into structured finished product quality detection data, thereby achieving a complete characterization of the final feed product's quality consistency and physical properties. For example, in a practical application scenario, the infrared protein analyzer is set to detect once per minute, the target protein content is 21.5%, and the crushing force (in N) of particles with a diameter range of 3.5 mm to 3.8 mm is recorded by a hardness tester every 5 seconds. The average batch hardness is assumed to be 45.0 N. For calculating the deviation rate of protein content, assuming the detection sequence values ​​are {21.1, 21.6, 21.4, 21.9} and the target value is 21.5%, the following formula is used:

[0013] in, To detect the standard deviation of the sequence, This represents the target protein content. It was obtained through calculation. =0.29, =21.5, with a deviation rate of approximately 0.0134, achieving a quantitative assessment of protein content consistency. Regarding the calculation of the particle hardness dispersion coefficient, assuming the detection sequence values ​​are {44.8, 45.2, 46.0, 44.0}, the batch average hardness... =45.0N, standard deviation =0.74, the formula for the coefficient of variation is:

[0014] The calculated coefficient of variation is approximately 0.0164, reflecting the stable state of the hardness distribution of this batch of products. After outlier removal, the above test results are bound with batch number BK20240608 and sampling time series to form standardized finished product quality test data, which is then input into S1.5. This provides accurate and interpretable terminal quality indicator data for subsequent adaptive filtering and causal modeling, significantly improving the diagnostic efficiency and reliability of results in the process of tracing the source of anomalies. S1.5: Based on the collected multi-source heterogeneous process parameters and quality inspection indicators, generate an original dataset containing process feature labels to provide structured input for subsequent adaptive filtering and time alignment processing.

[0015] Step S2: Adaptive filtering and time alignment processing are performed on the collected multi-source heterogeneous data. A sliding window strategy is used to segment and normalize the process parameters, generating a standardized time-series dataset containing process feature labels. The window length is dynamically adjusted according to the characteristic time constant of each stage. Specifically, this includes: S2.1: Perform missing value imputation and outlier removal operations on the collected multi-source heterogeneous raw process data. Based on the sliding window mid-value filtering algorithm, perform preliminary cleaning on the time series data of moisture content, temperature gradient curve and discrete value of mixing uniformity of each process to obtain the cleaned process parameter data sequence. For the collected multi-source heterogeneous raw process data, based on the time series data of moisture content, temperature gradient curve and discrete value of mixing uniformity output by the sensor, a missing value interpolation algorithm (parameters: interpolation type is linear, interpolation range threshold is 5 sampling intervals) is used to realize the continuous completion of the time series, so as to eliminate the data gaps caused by instantaneous sensor failure or communication delay during the acquisition process. Furthermore, by using an outlier removal method (parameters: 3σ principle, detection window length is the number of sampling points corresponding to the characteristic time constant of each process), outlier samples that exceed the normal process fluctuation range in the time series data of moisture content, temperature gradient curve and discrete value of mixing uniformity are removed, and valid data sequences that conform to the process rules are obtained. Furthermore, a sliding window midpoint filtering algorithm (parameters: window length dynamically set according to the process characteristic time constant, step size of 1 sampling point) is adopted to smooth the remaining valid data and generate a smooth sequence with reduced instantaneous noise interference. The filtering formula can be expressed as:

[0016] in, Half the length of the window. The original sequence, This is the filtered output sequence; Furthermore, a scale consistency check is performed on the filtered data (parameter: maximum allowable standard deviation fluctuation rate is 0.1) to achieve statistical scale uniformity of data from different processes, thereby ensuring unbiased processing in the subsequent normalization stage; By imputing missing values, removing outliers, filtering the middle value in the sliding window, and checking the scale consistency, the multi-source heterogeneous raw process data from the previous step is transformed into a preliminary cleaned data sequence with good time continuity, low noise level, and uniform scale, thus achieving the expected technical effect of subsequent sliding window length calculation and segmented normalization. For example, in the moisture content data collection during the raw material pretreatment stage, there was an approximately 5-second data acquisition interruption. Linear interpolation was selected as the interpolation algorithm, with an interpolation threshold corresponding to 5 sampling points. After interpolation, the temporal continuity of the data was significantly improved. During outlier removal, a standard deviation multiple threshold of 3 was set, and the detection window length was converted from the 10-minute time constant corresponding to this stage's characteristic time constant to a calculation of sampling points, eliminating several points with instantaneous fluctuations exceeding the normal range. The sliding window median filter had a window length of 21 sampling points and a step size of 1 sampling point, effectively smoothing small-amplitude spike noise. The above filtering was performed according to the formula... The standard deviation volatility of the output smoothed sequence is controlled at 0.08, which meets the threshold requirement of scale consistency check, and finally a stable and continuous preliminary cleaned dataset is obtained for the calculation of sliding window length. S2.2: Calculate the sliding window length parameter based on the characteristic time constant of each process. The raw material pretreatment process uses a 10-minute window, the mixing and batching process uses a 5-minute window, and the high-temperature granulation process uses a 3-minute window to match the dynamic response characteristics of each process and generate a set of window configuration parameters bound to the process type. S2.3: The process parameters after cleaning are segmented and normalized using the window configuration parameter set. The moisture content, temperature gradient and mixing uniformity of each process are transformed dimensionlessly based on the Z-score standardization method to obtain a standardized process parameter subsequence with a mean of 0 and a standard deviation of 1 within the time window. S2.4: Perform time alignment operation on the standardized process parameter subsequence, map multi-process data under different acquisition frequencies to a unified time base based on the timestamp matching algorithm, and use linear interpolation to compensate asynchronous sampling points for time, so as to obtain a process feature data set with consistent time axis. S2.5: Perform feature labeling on the time-aligned process feature data set. Based on the process specification database, add process type labels, time window numbers and feature dimension identifiers to each process data segment to generate a structured and standardized time series dataset to support the subsequent causal modeling module in identifying and processing process features. For the time-aligned process feature data set, a label mapping algorithm (parameters: process specification database, time window number index) is used to achieve automatic matching of process type labels. This algorithm establishes a one-to-one mapping relationship between the data fragment's acquisition source and the corresponding process type by retrieving a predefined association table between process steps and acquisition device identifiers in the process specification database, and obtains a process type label matrix. Furthermore, by using a time window encoding algorithm (parameters: time base mapping table, window length configuration set), the time window number of each data segment is generated, and the number is indexed and bound to the process type label to form a label-window lookup table to ensure the distinguishability of different process data within the same time period; Furthermore, by using the feature dimension identifier generation method (parameters: feature field definition set, structured data pattern), unique identifiers for feature data dimensions of each process are assigned, and a feature dimension index table is generated so that feature data can be quickly located when the causal modeling module is called later. Furthermore, through a label fusion processing algorithm (parameters: process type label matrix, time window number table, feature dimension index table), the generated process type labels, time window numbers, and feature dimension identifiers are merged into a unified structured meta information table, and field-level binding is performed with the standardized time series dataset to ensure that each column of feature data in the dataset is associated with complete label information. By using the label consistency verification method (parameters: process specification consistency rule set, data label meta information table), the results of the previous step are transformed into structured and standardized time-series data with a complete label system, so as to achieve accurate calling effect of the subsequent causal modeling module in process feature identification and path analysis; For example, in a pig feed production plant, the time window length for the raw material pretreatment stage is set to 10 minutes, for the mixing and batching stage to 5 minutes, and for the high-temperature pelleting stage to 3 minutes. The time-aligned moisture content, temperature gradient curve, mixing uniformity, finished product protein content deviation rate, and pellet hardness dispersion coefficient are used as five feature dimensions, and assigned dimension identifiers F1 to F5 respectively. The process specification database defines the sensor ID range for the raw material pretreatment stage as 1000-1999, the mixing and batching stage as 2000-2999, and the high-temperature pelleting stage as 3000-3999. The label mapping algorithm matches the acquisition device ID and labels the normalized data segments as three types of process labels: PRE_10, MIX_5, and PEL_3. The time window encoding algorithm calculates the window number based on the acquisition timestamp and window length; for example, the first acquisition window in the raw material pretreatment stage is encoded as W_PRE_001. The feature dimension identifier generation method identifies moisture content as F1, temperature gradient as F2, mixing uniformity as F3, protein content deviation rate as F4, and particle hardness dispersion coefficient as F5. After label fusion processing, a metadata table containing process type labels, time window numbers, and feature dimension identifiers is generated and bound to a standardized time-series dataset to form structured entry data. In subsequent causal modeling calls, the moisture content data of the first window in the raw material pretreatment stage can be directly located using label combinations such as (W_PRE_001, F1), achieving efficient and accurate path analysis.

[0017] like Figure 2 As shown, step S3 involves calculating the causal strength of variable pairs in the standardized time-series dataset based on the Granger causality test and the information transmission entropy joint criterion, and constructing an initial causal topology. Nodes represent process quality characteristics, edge weights represent the cross-process anomaly propagation strength, and directionality reflects the propagation time-series relationship. Specifically, this includes: S3.1: Perform Granger causality tests on the process quality characteristic variable pairs in the standardized time series dataset, model the lagged correlation of the time series, calculate the linear causal strength between variables, and identify process node pairs with significant linear causal relationships; The Granger causality test method (parameter: lag order p is adaptively determined according to the Akaike Information Criterion AIC) is used to model linear causal relationships based on time series for process quality characteristic variables in the standardized time series dataset. Furthermore, by constructing a bivariate vector autoregressive (VAR) model (parameters: including the current value and p-order lag terms), the time-series dependency fitting between variables is achieved, and multi-step prediction residual vector data for causality testing is obtained; Furthermore, the F-statistic calculation method (parameters: significance level α is 0.05, and the degree of freedom for testing is determined by the sample length and lag order) is used to compare the sum of squared predicted residuals of the restricted and unrestricted models, and to generate a linear causal strength index for the variable pairs.

[0018] in, For the sum of squared residuals of the constrained model, For the sum of squared residuals of the unconstrained model, To constrain the quantity, The total length of the sample. This represents the number of variables in the model. Furthermore, the F-statistic results are subjected to significance determination processing (parameter: the critical value is determined by the F-distribution table) to identify the significance of the linear causal relationship for each pair of process nodes, and to select the set of node pairs that meet the significance level requirements. By using the Granger causality test algorithm, the standardized time-series dataset from the previous step is transformed into process node pairs with significant linear causal strength and their corresponding strength indices, thereby achieving preliminary identification of cross-process linear causal propagation relationships. For example, in a standardized time-series dataset of a feed production batch, two process quality characteristics—raw material moisture content and mixing uniformity—are selected for Granger causality testing. A lag order of p=3 (chosen based on the minimum AIC) is set, and a bivariate VAR model is constructed. The model has degrees of freedom of sample length T=500, number of variables n=2, and number of constraints m=3. The sum of squared residuals of the constrained model is calculated. Unconstrained model residual sum of squares Substituting into the formula, we get the value of F. The critical value for the F-distribution at a significance level of α=0.05 is found in the F-distribution table. Since the F-value is greater than the critical value, it is determined that there is a significant linear causal relationship between the raw material moisture content and the mixing uniformity. The output is a process node pair (raw material moisture content → mixing uniformity), with a linear causal strength of 3.75, which serves as input data for subsequent nonlinear causal identification and causal topology construction. S3.2: Based on the Granger causality test results, the information transfer entropy method is used to calculate the nonlinear causal strength of process node pairs. The difference between mutual information and conditional mutual information is used to evaluate the nonlinear dependency between variables, so as to supplement the complex causal structure that linear methods cannot capture. S3.3: The Granger causality strength and information transmission entropy results are weighted and fused, and a comprehensive causality strength index is generated based on a dynamic weight allocation strategy. The weights are adaptively adjusted according to the nonlinearity of the quality characteristics of each process to improve the comprehensiveness and accuracy of causal identification. The input consists of the linear causal intensity index matrix and the nonlinear causal intensity index matrix obtained by S3.1 and S3.2 respectively. The two are organized in units of process quality characteristic variable pairs and numbered according to time windows. A dynamic weight allocation method (parameters: nonlinearity threshold set, weight adjustment coefficient set) is adopted to achieve adaptive fusion of linear causal intensity and nonlinear causal intensity; Furthermore, the degree of nonlinearity for each variable pair is quantified by using a nonlinearity evaluation algorithm (parameters: mutual information statistic of variable pairs, Pearson correlation coefficient), and a nonlinearity score matrix is ​​obtained. Furthermore, by using a weight mapping function (parameters: nonlinearity score, preset mapping curve parameter set), the nonlinearity score is converted into a fusion weight of linear and nonlinear components, generating a weight matrix to control the fusion ratio. Furthermore, a weighted fusion calculation formula is used to fuse the two types of causal strength element by element. The calculation formula is as follows:

[0019] in, To form a comprehensive causal strength index matrix, To fuse the weight matrix, It is a linear Granger causality strength matrix. This is the nonlinear information transmission entropy intensity matrix; Furthermore, a normalization algorithm (parameters: the maximum and minimum values ​​of the fused index matrix) is used to make the comprehensive causal strength index dimensionless, so that it can be used for subsequent construction of a directed weighted causal graph. Through the above dynamic weight allocation and weighted fusion processing, the results of the previous step are transformed into a causal strength index that reflects the combined effect of linear and nonlinear dependencies, thereby achieving the expected technical effect of high coverage and accuracy in causal identification. For example, in a standardized time-series dataset of a certain pig feed production batch, the Granger causality strength of the deviation rate between raw material moisture content and finished product protein content is 0.68, the nonlinear information transfer entropy strength is 0.45, and the score obtained after nonlinearity evaluation is 0.40. This score is converted into a fusion weight w=0.35 through a weight mapping function. The comprehensive causality strength is calculated using the weighted fusion formula: The result was 0.53, which was normalized to 0.76. After multivariate pair fusion and normalization, the final output comprehensive causal intensity matrix showed significant differences among the variable pairs of different processes. The high intensity index was concentrated between the raw material moisture content and mixing uniformity, and between mixing uniformity and particle hardness. This verifies that the fusion strategy can significantly improve the comprehensiveness and accuracy of anomaly propagation path identification, and in this scenario, it achieves effective quantification of cross-process quality anomaly correlation. S3.4: Construct an initial causal topology based on the comprehensive causal strength index, take process quality characteristics as graph nodes, causal strength as edge weights, and determine the edge direction according to the time lag relationship between variables to form a directed weighted causal graph. Using the comprehensive causal strength index obtained by weighted fusion as input, node mapping is performed on the process quality feature variables in the standardized time series dataset. A node construction method (parameters: process feature label, quality detection index set) is adopted to assign each different process quality feature as an independent node in the initial causal topology. Furthermore, by using the edge weight assignment method (parameter: comprehensive causal strength index matrix), the comprehensive causal strength values ​​of each variable pair are mapped to the edge weights between corresponding nodes, and weighted directed connection relationship data is obtained. Furthermore, a time lag analysis method (parameter: the significant lag duration of each variable pair) is used to determine the directionality of each weighted edge and generate a set of edges containing directional attributes, where the direction points from the variable with smaller lag to the variable with larger lag. Furthermore, a preliminary causal topological structure is constructed and a topological structure dataset is generated through a graph structure generation algorithm (parameters: node set, weighted directed edge set); Furthermore, by using a causal graph construction tool (parameters: topology dataset, node temporal attributes), the topology is visualized as a directed weighted causal graph, and the attribute information of nodes and edges is labeled. Through the above graph structure generation and attribute annotation processing methods, the comprehensive causal strength result of the previous step is transformed into a directed weighted causal graph that can be used for subsequent dynamic optimization, thereby realizing the preliminary identification and structured representation of cross-process anomaly propagation paths. For example, in a batch analysis of feed production, the standardized time-series dataset includes variables such as raw material moisture content, high-temperature pelleting temperature gradient, mixing uniformity, and finished product protein content deviation rate. The node construction method maps these variables to four nodes. The comprehensive causal strength index matrix is ​​a 4×4 real number matrix with elements ranging from 0.15 to 0.92. The edge weight assignment method directly uses these values ​​as the edge weights. The time lag analysis is based on the lag time calculated within the sliding window. The lag between raw material moisture content and mixing uniformity is 5 minutes, so the directionality points to the mixing uniformity node. The lag between high-temperature pelleting temperature gradient and finished product protein content deviation rate is 3 minutes, so the directionality points to the finished product protein content deviation rate node. The graph structure generation algorithm combines the node set with the directed weighted edge set to obtain a causal topology structure containing six directed weighted edges. In the directed weighted causal graph generated by the causal graph construction tool, each edge and node is labeled with a weight or lag time. The final output graph clearly represents the key anomaly propagation link from raw material pretreatment to finished product, achieving high-precision identification of the initial causal structure. S3.5: The initial causal graph is sparsified by using a threshold pruning strategy to remove weak connections with causal strength below a set threshold, thereby optimizing the graph structure, improving model interpretability, and reducing the computational complexity of subsequent dynamic updates. The initial causal graph, constructed based on a comprehensive causal strength index, takes as input the quality feature nodes of each process and their directed weighted edges. A threshold pruning algorithm is then employed (parameter: pruning threshold). The weights of all edges in the graph are compared one by one, and weights less than 1 are removed. The edges are used to eliminate weak causal relationships; Furthermore, by using the edge weight normalization method (parameter: normalization modulus range [0,1]), the retained edge weights are scaled according to the maximum value ratio to ensure the comparability of causal strength between different processes and obtain a normalized weighted directed graph data structure. Furthermore, based on the node degree distribution equalization algorithm (parameter: threshold for the difference between node in-degree and out-degree), Adjusting the edge retention ratio of highly connected nodes reduces the interference of high-density regions on model interpretability, resulting in a structurally balanced sparse graph. Furthermore, a connectivity-preserving constraint strategy is adopted (parameter: minimum connected subgraph size). During the pruning process, the set of edges that meet the minimum connectivity requirement is preserved to avoid the breakage of causal propagation paths due to pruning, thereby maintaining the integrity of key causal paths. Furthermore, the sparsity index of the graph is calculated using a complexity evaluation method for the sparsified graph (parameters: number of nodes N, number of edges E). This ensures that the sparsity is within a preset range, thereby reducing the computational load of subsequent dynamic update processes; By combining the above-mentioned threshold pruning, weight normalization, node degree equalization and connectivity preservation processing methods, the initial causal graph of the previous step is transformed into sparse causal graph data with optimized structure and reduced computational complexity, thereby improving the interpretability of the model and improving the computational efficiency. For example, in an initial causal graph containing 50 process quality feature nodes and 200 weighted edges, a pruning threshold is set. The weight is set to 0.15. Eighty edges with a weight less than 0.15 are removed, leaving 120 edges. The retained edges are then normalized, with the maximum weight scaled from 0.92 to 1.0, and the remaining edges adjusted proportionally. A threshold for the difference in node in-degree is set. Set the value to 5, reduce the number of edges on nodes with an in-degree greater than 20 by 10%, thus reducing the total number of in-edges and out-edges on high-density nodes. Set the minimum connected subgraph size. Set the value to 3 to ensure that each process node maintains a causal connection with at least two other nodes after pruning. Calculate the sparsity index, derived from the sparsity of the original graph. =0.163 decreased to sparseness after pruning =0.098, realizing the transformation of the graph structure from dense to balanced sparse, significantly improving the clarity of causal path presentation and reducing the computational cost during dynamic updates.

[0020] like Figure 3 As shown, step S4 involves dynamically optimizing the initial causal topology using an improved time-aware structural equation model (T-SEM), updating the causal edge weights in real time using a variational Bayesian inference method, and generating a causal graph sequence that changes with production conditions. The model parameters are iteratively updated based on the latest quality inspection data within the sliding window. Specifically, this includes: S4.1: Based on the latest quality inspection data and process parameter time series data within the sliding window, initialize the node variables and structural paths of the time-aware structural equation model (T-SEM). The node variables include raw material moisture content, mixing uniformity, high-temperature granulation temperature gradient, finished product protein content deviation rate, and particle hardness dispersion coefficient. The structural paths are initialized based on the initial causal topology constructed in S3 to establish the initial state of the model. Based on the latest quality inspection data and process parameter time series data within the sliding window, a node variable definition method (parameters: moisture content, mixing uniformity, high-temperature granulation temperature gradient, protein content deviation rate, particle hardness dispersion coefficient) is adopted to initialize the node set of the time-aware structural equation model (T-SEM). Furthermore, the model structure path matrix is ​​initialized through the structure path mapping method (parameter source: the initial causal topology constructed in step S3), and the consistency of node and path indices is maintained to ensure the traceability of the causal relationship graph in the time dimension. Furthermore, a process feature binding strategy is adopted (based on process type labels and feature dimension identifiers) to realize the mapping and binding of node variables with specific process steps, and generate a node-process mapping table to support the differentiated processing of different processes in the subsequent parameter update process; Furthermore, by employing a numerical normalization method (using Z-score standardization), the initial node values ​​are made dimensionless, and the mean value is obtained. Standard deviation is The node state vectors provide scale-consistent data input for variational Bayesian inference; Furthermore, a parameter initialization function is used. (in (For the historical observation count of the nodes), the initial values ​​of the structural path coefficients are set, and a model state data package containing node variables, structural paths and initial coefficients is generated; Through the above initialization process, the latest quality inspection data and process parameter time series data within the sliding window are transformed into a T-SEM initial state with time traceability, structural integrity and scale consistency, thus realizing the foundation for interpretability of the model in dynamic optimization. For example, in a production batch where the moisture content of the raw materials is 8.2%, the mixing uniformity dispersion value is 0.12, the high-temperature granulation temperature gradient range is 85℃ to 95℃, the finished product protein content deviation rate is 0.03, and the particle hardness dispersion coefficient is 0.05, a sliding window length of 10 minutes is selected, and the node sets are initialized as follows: , , , , The path matrix is ​​constructed using the causal topology output from step S3, and then processed according to the formula. When assigning values ​​to the path coefficients, the initial values ​​of the coefficients are as follows: when the number of historical observations of a node is 10, 8, 12, 15, and 9, respectively. , , , , After the initialization of this batch, the model state data package contains complete node variable vectors, path matrices, and coefficient sets. It has been verified that in the subsequent variational Bayesian inference stage, the model can converge stably and significantly improve the accuracy and interpretability of anomaly propagation path identification. S4.2: Perform variational Bayesian inference on the structural path coefficients in the T-SEM model. Use the standardized time series dataset within the sliding window to estimate the approximate posterior distribution of the model parameters to obtain the optimal causal edge weights under the current working condition. The inference process adopts the mean field approximation strategy to decompose the joint distribution into independent approximate distributions of each path coefficient to reduce computational complexity. For the initialized time-aware structural equation model (T-SEM) node variables and structural paths, the variational Bayesian inference method (parameters: sliding window width, standardized time series dataset) is used to approximate the posterior distribution of the structural path coefficients, thereby achieving dynamic optimization of the causal edge weights under the current production conditions. Furthermore, by utilizing a mean-field approximation strategy (parameters: set of path coefficients, initial causal topology), the joint probability distribution is... Decomposed into multiple independent univariate approximate distributions This reduces the computational complexity of joint inference; Furthermore, by minimizing the variational free energy The optimization objective is to calculate the difference between the expected log-likelihood and KL divergence of the standardized time series dataset under the model structure, and to perform gradient updates on the approximate posterior distribution parameters corresponding to the coefficients of each structural path. Furthermore, the observed values ​​of each node variable within the sliding window The expected value matrix of causal edge weights is generated by performing a weighted convolution with the approximate distribution parameters of the path coefficients. Through the above inference process, the current working condition data of the sliding window is transformed into the optimal edge weight distribution that reflects the causal strength between variables, thereby realizing the dynamic interpretability update of the causal graph. For example, in the third sliding window of a certain batch of feed production, the window length is 3 minutes, and the corresponding standardized time-series dataset includes three node variables: raw material moisture content, high-temperature pelleting temperature gradient, and mixing uniformity. In the mean-field approximation strategy, the initial distribution of each path coefficient is set to a normal distribution. The learning rate is 30 steps through gradient updates. Free energy is generated during the iteration process. Descending to The model achieves stable convergence. In the expected value matrix of edge weights generated by this window, the edge weights pointing from the mixing uniformity to the deviation rate of the finished product protein content are significantly increased, reflecting the enhanced impact of this process on quality anomalies under the current operating conditions. The model clearly reveals the dynamic changes of this path in terms of interpretability. S4.3: The T-SEM model is updated with parameters iteratively based on the variational free energy minimization criterion. The coordinate ascent method is used to optimize the coefficients of each path one by one until they converge to a local optimum, thereby obtaining a dynamic causal map that reflects the current production status. The convergence condition is set as the change in free energy between two adjacent iterations being less than a preset threshold. S4.4: Bind the updated causal edge weights to timestamps to generate a time-labeled causal graph sequence, where each graph corresponds to a time segment of a sliding window, to support the tracing and comparative analysis of historical causal relationships in subsequent counterfactual reasoning; S4.5: Perform graph structure stability assessment on the causal graph sequence. Calculate the graph evolution trend index based on the graph node degree distribution and edge weight change rate. When the evolution trend index exceeds the preset fluctuation threshold, trigger the causal structure relearning mechanism, re-execute the Granger causality test and the information transfer entropy joint criterion to update the initial causal topology and feed it back to the T-SEM model, thereby achieving closed-loop optimization of the causal graph. Based on the node set and edge weight time series data of the causal graph sequence, a node degree distribution evaluation algorithm (parameters: causal graph node set, edge set, timestamp) is used to realize the quantitative analysis of the redundancy and centralization of the graph structure in the current time segment. Furthermore, by using the edge weight change rate calculation method (parameters: edge weight sequence within the sliding window, time interval), the fluctuation of the causal relationship strength over time is measured, and the change rate vector result for each edge is obtained; Furthermore, using a comprehensive evaluation formula for graph evolution trends, the standard deviation of node degree distribution and the mean of edge weight change rates are normalized and weighted to generate an evolution trend index characterizing the overall structural stability. The formula is as follows:

[0021] in, The standard deviation of the node degree distribution. The mean of the rate of change of edge weights. These are the weighting coefficients. As an indicator of evolutionary trends; Furthermore, by setting a fluctuation threshold determination function (parameter: preset threshold) Evolutionary trend indicators This enables the detection of out-of-limit states in the evolution stability of causal graphs and outputs an out-of-limit flag signal. Furthermore, based on the out-of-limit flag signal, the causal structure relearning module is invoked to re-execute the joint criterion calculation of Granger causality test and information transmission entropy, and update the initial causal topology data; By feeding the relearned topological edge set and node set back to the T-SEM model, the original structural path is covered, and the dynamic optimization process of the next iteration cycle is triggered, thereby achieving closed-loop optimization of the causal graph and adaptive maintenance of structural stability. For example, in the cloud-based analysis of a pig feed production batch, weighting coefficients are set. for Preset threshold for In a sliding window with a length of Standard deviation of node degree distribution in a causal graph sequence of minutes Calculated Mean rate of change of edge weights for Substitute the values ​​into the formula to calculate the evolution trend index: The result is Because this value exceeds the preset threshold. This triggers a causal structure relearning mechanism. The system recalculates Granger causality and information transfer entropy, updates the topology, significantly improves the stability of the graph in subsequent analysis cycles and the accuracy of abnormal path identification, and ensures that the causal graph input to the counterfactual reasoning stage maintains high confidence.

[0022] Step S5: When the finished product quality deviation rate is detected to exceed a preset threshold, a counterfactual reasoning mechanism is activated. A virtual intervention matrix is ​​constructed based on the causal graph sequence. By calculating the hypothetical intervention effect of 'if no parameter shift occurred in a specific upstream process,' the attribution contribution of each process to the final anomaly is quantified. Specifically, this includes: S5.1: Based on the dynamic causal edge weights between nodes in the causal graph sequence, a virtual intervention matrix is ​​constructed, where the matrix rows represent upstream process nodes, the columns represent downstream process nodes, and the matrix elements are the normalized propagation strengths of the corresponding causal edges, in order to support the quantitative calculation of subsequent counterfactual intervention effects. The causal graph sequence output from step S4 is used as input conditions, and the dynamic causal edge weights between nodes in the causal graph are selected as the basic data for constructing the matrix. A node index mapping method is used (parameter: mapping upstream process node identifiers to matrix row indices, and mapping downstream process node identifiers to matrix column indices) to achieve a bidirectional mapping relationship between graph nodes and matrix coordinates; Furthermore, a normalization algorithm (parameters: edge weight dataset, normalization range set to 0 to 1) is used to normalize the propagation strength of each causal edge, resulting in a normalized weight set. This calculation process uses the following formula:

[0023] in, To normalize the propagation intensity, The original causal edge weights, The minimum value of the weight set. The maximum value of the weight set; Furthermore, by using a matrix filling method (parameter: matrix size = number of upstream nodes × number of downstream nodes), the normalized weight set is filled into the corresponding matrix positions, and a full virtual intervention matrix is ​​generated. Furthermore, by using a matrix sparsity control algorithm (parameter: threshold is 0.05), matrix elements with propagation strength below the threshold are removed and assigned zero, thereby reducing subsequent computational complexity and improving the significance of the intervention effect. Through the matrix construction and sparsification methods described above, the dynamic causal edge weights extracted in the previous step are transformed into a structured virtual intervention matrix, thus providing a quantitative input basis for the calculation of the counterfactual intervention effect in S5.2. For example, in a pig feed production quality anomaly tracing application scenario, the causal graph sequence includes five types of nodes: raw material moisture content, mixing uniformity, high-temperature pelleting temperature gradient, finished product protein content deviation rate, and pellet hardness dispersion coefficient. The dynamic causal edge weights range from 0.12 to 0.87. A 5×5 matrix coordinate system is established using a node index mapping method, and the weights between node pairs are filled into the corresponding matrix coordinates. During normalization processing... =0.12, =0.87, for example, if the original weights of the upstream node "raw material moisture content" and the downstream node "mixing uniformity" are 0.56, then the calculated normalized propagation intensity is 0.87. =0.5867. All normalized results were filled into the matrix and sparsified, removing elements less than 0.05 to form the final virtual intervention matrix. In this matrix, non-zero cells only retain the intensity values ​​corresponding to significant propagation paths, ensuring the controllability of computational costs in subsequent counterfactual intervention analysis and improving the sensitivity of path attribution. S5.2: Perform a hypothetical intervention operation on each row of the virtual intervention matrix, that is, set the upstream process variable value corresponding to the row to the baseline state (no parameter shift occurs), and perform forward propagation simulation based on the causal graph structure to obtain the expected quality output of the downstream node under the intervention condition; S5.3: Calculate the difference between the expected quality output under the intervention conditions and the actual detected quality deviation rate to obtain the abnormal mitigation effect of each intervention path, which serves as the potential contribution index of the upstream process to the final abnormal result. S5.4: Based on the anomaly mitigation effect size, combined with the path length and propagation delay information in the causal graph, perform weighted attribution calculations on each intervention path to generate an attribution contribution vector containing the propagation path weight and the intensity of causal influence, so as to reflect the comprehensive influence of different upstream processes on the finished product anomalies. Based on the abnormal mitigation effect size data obtained in step S5.3 and the structural information in the causal graph sequence, a propagation path weighted attribution algorithm (parameters: abnormal mitigation effect size matrix, path length vector, propagation delay vector) is used to realize the comprehensive attribution calculation for each intervention path; Furthermore, by using the path length normalization method (parameter: number of nodes in each path), the path length is used to correct the attribution contribution, appropriately amplifying the effect size of shorter paths and appropriately attenuating the effect size of longer paths, and obtaining a length-weighted factor matrix. Furthermore, by using a propagation delay weighting calculation method (parameter: the sum of delays of each path in the causal graph), the contribution of delay to attribution is adjusted, giving higher weights to paths with shorter delays and lower weights to paths with longer delays, and generating a delay weighting factor matrix. Furthermore, a comprehensive weighted fusion method (parameters: length-weighted factor matrix, time-delay-weighted factor matrix, and anomaly mitigation effect matrix) is employed to achieve the synergistic effect of multiple indicator weights. The anomaly mitigation effect is multiplied according to the corresponding length and time-delay weighting factors to obtain the final weighted anomaly mitigation effect matrix. The calculation formula is as follows:

[0024] in, This is a weighted abnormal mitigation effect matrix. This is the abnormal mitigation effect size matrix. This is the path length weighting factor matrix. The propagation delay weighting factor matrix; Furthermore, an attribution vector generation algorithm (parameter: weighted anomaly mitigation effect matrix) is adopted to aggregate the matrix according to the path, sum the weighted anomaly mitigation effects of the same upstream process, and generate an attribution contribution vector containing the propagation path weight and the causal influence intensity. By using the above algorithm, the anomaly mitigation effect of the previous step is transformed into a structured attribution contribution vector, thereby quantifying the degree of influence of different upstream processes on finished product anomalies under the condition of comprehensively considering path length and propagation delay. For example, in a certain pig feed production batch, the anomaly mitigation effect for deviations in raw material pretreatment and moisture content is 0.35, the anomaly mitigation effect for deviations in mixing uniformity is 0.28, and the anomaly mitigation effect for deviations in high-temperature pelleting temperature gradient is 0.42. The path lengths are 3, 4, and 2 nodes, and the propagation delays are 12 minutes, 18 minutes, and 8 minutes, respectively. The length weighting factor is set to... The time delay weighting factor is set to The weighted calculation process is as follows: The calculation results show that the attribution contribution of the first path is 0.441, the attribution contribution of the second path is 0.292, and the attribution contribution of the third path is 0.579. The aggregated attribution contribution vector reflects that the third path has a significant impact, higher than the other paths. The model can significantly improve the accuracy of anomaly source localization in multiple batch validations. S5.5: Perform a significance test on the attribution contribution vector, and select process nodes and their propagation paths whose contribution exceeds the preset significance threshold to form a set of high-impact anomaly propagation paths, which will serve as the core input data for subsequent semantic annotation and visual tracing.

[0025] Step S6: Input causal paths with attribution contribution exceeding a preset significance level into the expert rule knowledge base, perform semantic annotation and path compression operations, and generate a visualized tracing report containing anomaly propagation timing features, key nodes, and propagation intensity. The annotation is based on process specification constraints and a historical failure mode library. Specifically, this includes: S6.1: Perform path filtering processing on causal paths whose attribution contribution exceeds the preset significance level. Based on the temporal relationship of nodes in the causal graph and the propagation intensity threshold, extract a set of key propagation paths with significant influence to remove redundant causal paths and improve the efficiency of subsequent semantic annotation. S6.2: Based on process specification constraints and historical failure mode library, construct a semantic mapping rule set in expert rule knowledge base. The semantic mapping rules include typical failure modes, anomaly propagation timing feature labels and process parameter deviation types corresponding to each process node, so as to realize semantic annotation of causal path. Based on the set of key propagation paths whose attribution contribution exceeds a preset significance level, the input conditions include the selected causal path sequences and their nodes' process parameter deviation types, propagation time sequence information, and historical quality inspection data index; A rule extraction method (parameters: process specification database, failure mode library index, node process link code) is used to parse the link identifier of the process node and the corresponding process specification constraint, so as to build a node-to-process link mapping table. Furthermore, by using a pattern matching algorithm (parameters: anomaly type encoding rules, historical fault mode template set), the similarity between the process parameter deviation features of each node and the typical patterns in the historical fault mode library is calculated, and the typical fault mode label set corresponding to the node is obtained. Furthermore, a temporal feature extraction method (parameters: propagation path timestamp sequence, stage division rules) is adopted to analyze the propagation interval and sequence relationship between nodes in the causal path, and generate abnormal propagation temporal feature labels, including three types of labels: initiation stage, intermediate diffusion stage and terminal result stage. Furthermore, by applying a parameter classification method (parameters: deviation quantification index, deviation type definition table), the process parameter offset values ​​of nodes are mapped to process parameter deviation type labels, thereby realizing the classification of parameters with different dimensions in a unified label space. Through the above rule construction process, the set of key propagation paths in the previous step is transformed into a semantic mapping rule set with process link names, typical failure modes, abnormal propagation timing feature labels and process parameter deviation types, so as to realize the generation of basic data for semantic annotation of causal paths. For example, in the high-temperature pelleting stage of feed production, the node parameter deviation is the increase of the pelleting temperature gradient. ℃, this step is coded as HDG-03 in the process specification database, and the corresponding specification constraint is that the temperature gradient should be maintained at ℃. The temperature range is specified. A pattern matching algorithm retrieves a historical fault mode database, identifies a typical pattern encoding FM-HDG-TEMP, tagged as "temperature overshoot leading to protein denaturation." Similarity is calculated based on the cosine similarity between the deviation amplitude and the pattern feature vector. And select the label. Time-series feature extraction and analysis of the propagation path's timestamp sequence revealed that this node is located at the [number]th [timest point] after the anomaly occurred. The first hour, belonging to the intermediate diffusion stage, is labeled STAGE-MID. The parameter classification method is based on the temperature offset definition table. The ℃ offset is categorized as a high-amplitude temperature deviation type: TYPE-TEMP-HIGH. Ultimately, this results in the following entries in the semantic mapping rule set: process node "high-temperature granulation stage"—fault mode label "temperature overshoot leading to protein denaturation"—temporal feature label "intermediate diffusion stage"—process parameter deviation type "high-amplitude temperature deviation". This rule set provides a precise label data foundation for subsequent semantic annotation and visualization report generation, significantly improving the interpretability and credibility of the traceability results. S6.3: Perform semantic annotation on each path in the selected set of critical propagation paths. Based on the semantic mapping rule set in the expert rule knowledge base, annotate each node and edge in the path with the corresponding process link name, anomaly type label and propagation time sequence stage to generate a structured set of semantically enhanced causal paths. S6.4: Perform path compression processing on the semantically enhanced causal path set. Based on the consistency of the topological structure and semantic labels between nodes in the path, merge continuous path segments with the same process links and abnormality type labels to reduce the number of paths and improve the readability of the visual traceability report. S6.5: Based on the compressed semantically enhanced causal path set, generate a visual source tracing report. The visual source tracing report includes abnormal propagation time sequence features, key nodes and their semantic labels, propagation intensity heatmap and attribution contribution distribution map to support users in quickly identifying and analyzing the root causes of feed quality abnormalities. Based on the compressed semantically enhanced causal path set, a multi-layer temporal mapping construction method (parameters: timestamp of path node, process step identifier, and anomaly type label) is used to generate the anomaly propagation temporal feature matrix. Furthermore, by using a node semantic aggregation algorithm (parameters: key node set, label consistency coefficient), the binding mapping between key nodes and their semantic labels is realized, and a node semantic description table is obtained; Furthermore, by using a propagation intensity calculation method (parameters: causal edge weights, time delay coefficients), the path propagation intensity matrix is ​​generated, forming a two-dimensional numerical array that can be used for heatmap rendering; Furthermore, an attribution contribution distribution calculation method (parameters: attribution contribution vector, path weight coefficient) is adopted to construct the dataset of the attribution contribution distribution map; By using visualization rendering, the temporal feature matrix, node semantic description table, propagation intensity matrix and attribution contribution distribution data generated in the previous step are transformed into a visualization traceability report, enabling rapid identification and analysis of the root causes of abnormal feed quality. For example, in a scenario containing a compressed semantically enhanced causal path set, the path set includes 5 significant anomaly propagation paths. The timestamps of each node are accurate to the second level. The process steps are identified as raw material pretreatment, mixing and batching, high-temperature granulation, and finished product inspection. Anomaly type labels include high moisture content, low temperature, uneven mixing, and low protein content. A multi-layer temporal mapping construction method is used to normalize the timestamps of the path nodes to a unified time base of 0-3600 seconds, generating a 5×4 temporal feature matrix. In the node semantic aggregation algorithm, the label consistency coefficient is set to 0.85, and consecutive nodes with the same anomaly type label are merged to obtain a compressed node semantic description table. A propagation intensity calculation method is used, based on the causal edge weight range of 0.2-0.9 and the time delay coefficient range of 1-5 seconds, to calculate the propagation intensity value of each path, forming a 5×3 intensity matrix. For example, the intensity value of a certain path in the intensity matrix is... This indicates that the path is at a high intensity level in the anomaly propagation chain. The attribution contribution distribution calculation method is based on the attribution contribution vector [ , , , The path weight coefficients range from 0.5 to 1.0, and the attribution score for each node is calculated. For example, in path 3, the node attribution score is... This indicates that it has a strong impact on finished product anomalies. Inputting the above matrix and dataset into a visualization rendering engine generates a visual source tracing report containing a line graph of anomaly propagation time-series characteristics, a key node label mapping table, a propagation intensity heatmap, and an attribution contribution distribution histogram. This significantly improves the efficiency of anomaly root cause localization for inspection personnel under multi-process conditions.

[0026] The technical solution of the present invention has been described above with reference to the preferred embodiments shown in the accompanying drawings. However, it will be readily understood by those skilled in the art that the scope of protection of the present invention is obviously not limited to these specific embodiments. Without departing from the principles of the present invention, those skilled in the art can make equivalent changes or substitutions to the relevant technical features, and the technical solutions after these changes or substitutions will all fall within the scope of protection of the present invention.

[0027] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and rules of the present invention should be included within the scope of protection of the present invention.

Claims

1. A cloud computing-based method for tracing the quality of pig feed production, characterized in that, Includes the following steps: S1: Collect process parameters and quality inspection indicators for each stage of the feed production process, including raw material pretreatment, mixing and batching, and high-temperature pelleting. The process parameters include time-series data of moisture content, temperature gradient curves, and dispersion values ​​of mixing uniformity. The quality inspection indicators include the deviation rate of finished product protein content and the dispersion coefficient of pellet hardness. S2: Perform adaptive filtering and time alignment processing on the collected multi-source heterogeneous data, and use a sliding window strategy to segment and normalize the process parameters to generate a standardized time series dataset containing process feature labels. S3: Based on the joint criterion of Granger causality test and information transfer entropy, calculate the causal strength of the variable pairs in the standardized time series dataset and construct the initial causal topology. S4: The initial causal topology is dynamically optimized using an improved time-aware structural equation model. The causal edge weights are updated in real time using a variational Bayesian inference method to generate a causal graph sequence that changes with production conditions. S5: When the finished product quality deviation rate is detected to exceed the preset threshold, the counterfactual reasoning mechanism is activated, a virtual intervention matrix is ​​constructed based on the causal graph sequence, the hypothetical intervention effect is calculated, and the attribution contribution of each process to the final abnormality is quantified. S6: Input the causal path whose attribution contribution exceeds the preset significance level into the expert rule knowledge base, perform semantic annotation and path compression operations, and generate a visual source tracing report containing abnormal propagation time sequence features, key nodes and propagation intensity.

2. The method for traceability of pig feed production quality based on cloud computing according to claim 1, characterized in that, Step S1 specifically includes: Based on the moisture and temperature sensors in the raw material pretreatment process, time-series data of raw material moisture content and temperature gradient curves are collected. Based on the weighing system and mixing uniformity detection device in the mixing and batching process, the discrete value data of mixing uniformity are collected. Based on the temperature controller and pressure sensor in the high-temperature granulation process, temperature change curves and pressure fluctuation data are collected during the granulation process. Based on the finished product quality testing equipment, the deviation rate of protein content and the dispersion coefficient of particle hardness of the finished product were collected. Based on the multi-source heterogeneous process parameters containing the time-series data of the raw material moisture content, the temperature gradient curve, and the discrete value of the mixing uniformity, and the quality inspection indicators containing the deviation rate of the finished product protein content and the discrete coefficient of particle hardness, an original dataset containing process feature labels is generated.

3. The method for traceability of pig feed production quality based on cloud computing according to claim 2, characterized in that, The time-series data of raw material moisture content and temperature gradient curve were acquired using moisture and temperature sensors, and were collected by an analog-to-digital conversion acquisition module and verified by cyclic redundancy. Moisture data was filtered for noise reduction using a moving average, and temperature data was fitted with a second-order polynomial to generate a temperature gradient curve. The state of the raw materials was then quantitatively evaluated using process stability criteria.

4. The method for traceability of pig feed production quality based on cloud computing according to claim 1, characterized in that, Step S2 specifically includes: Missing value imputation and outlier removal operations are performed on the collected multi-source heterogeneous process parameter data. Based on the sliding window mid-value filtering algorithm, the time series data of moisture content, temperature gradient curve and discrete value of mixing uniformity of each process are initially cleaned to obtain the cleaned process parameter data sequence. The sliding window length parameter is calculated based on the characteristic time constant of each process, and a set of window configuration parameters bound to the process type is generated. The window configuration parameter set is used to perform segmented normalization processing on the process parameter data sequence after cleaning. Based on the standardization method, the moisture content, temperature gradient and mixing uniformity discrete values ​​of each process are transformed into dimensionless values ​​to obtain a standardized process parameter subsequence. Time alignment is performed on the standardized process parameter subsequence. Based on the timestamp matching algorithm, multi-process data at different acquisition frequencies are mapped to a unified time base. Linear interpolation is used to compensate asynchronous sampling points for time, thereby obtaining a set of process feature data with a consistent time axis. The feature labeling operation is performed on the process feature data set. Based on the process specification database, process type labels, time window numbers and feature dimension identifiers are added to each process data segment to generate a structured and standardized time series dataset.

5. The method for traceability of pig feed production quality based on cloud computing according to claim 1, characterized in that, Step S3 specifically includes: Granger causality tests were performed on the process quality characteristic variable pairs in the standardized time series dataset. Based on the lagged correlation model of the time series, the linear causal strength between variables was calculated, and process node pairs with significant linear causal relationships were identified. Based on the Granger causality test results, the information transfer entropy method is used to calculate the nonlinear causal strength of the process node pairs, and the difference between mutual information and conditional mutual information is used to evaluate the nonlinear dependency between variables. The Granger causality test results and the information transmission entropy results are weighted and fused to generate a comprehensive causality strength index based on a dynamic weight allocation strategy. An initial causal topology is constructed based on the comprehensive causal strength index, forming a directed weighted causal graph; The directed weighted causal graph is sparsified, and a threshold pruning strategy is used to remove weak connections with causal strength below a set threshold.

6. The method for traceability of pig feed production quality based on cloud computing according to claim 5, characterized in that, The initial causal topology structure includes using process quality characteristics as graph nodes, causal strength as edge weights, and determining edge directions based on the time lag relationship between variables.

7. The method for traceability of pig feed production quality based on cloud computing according to claim 1, characterized in that, Step S4 specifically includes: Based on the latest quality inspection data and process parameter time series data within the sliding window, initialize the node variables and structural paths of the time-aware structural equation model and establish the initial state of the model; The variational Bayesian inference process is performed on the structural path coefficients in the time-aware structural equation model. The model parameters are approximate posterior distribution estimation is performed using the standardized time series dataset within the sliding window to obtain the optimal causal edge weights under the current working condition. The time-aware structural equation model is iteratively updated based on the variational free energy minimization criterion. The coefficients of each path are successively optimized using the coordinate ascent method until they converge to a local optimum, thus obtaining a dynamic causal graph reflecting the current production status. The updated causal edge weights are bound to timestamps to generate a causal graph sequence with time tags; The graph structure stability assessment is performed on the causal graph sequence. The graph evolution trend index is calculated based on the graph node degree distribution and edge weight change rate. When the evolution trend index exceeds the preset fluctuation threshold, the causal structure relearning mechanism is triggered, the Granger causality test and the information transmission entropy joint criterion are re-executed, the initial causal topology is updated and fed back to the time-aware structural equation model.

8. The method for traceability of pig feed production quality based on cloud computing according to claim 7, characterized in that, In the variational Bayesian inference process, the inference process adopts a mean-field approximation strategy, which decomposes the joint distribution into independent approximate distributions of the path coefficients.

9. A cloud computing-based method for tracing the quality of pig feed production according to claim 1, characterized in that, Step S5 specifically includes: A virtual intervention matrix is ​​constructed based on the dynamic causal edge weights between nodes in the causal graph sequence; Perform a hypothetical intervention operation on each row of the virtual intervention matrix, set the upstream process variable value corresponding to the row as the baseline state, and perform forward propagation simulation based on the causal graph structure to obtain the expected quality output of the downstream node under the intervention condition. The difference between the expected quality output and the actual detected quality deviation rate is calculated to obtain the abnormal mitigation effect size under each intervention path. Based on the aforementioned abnormal mitigation effect size, and combined with the path length and propagation delay information in the causal graph, a weighted attribution calculation is performed on each intervention path to generate an attribution contribution vector containing the propagation path weight and the causal influence intensity. A significance test is performed on the attribution contribution vector to screen out process nodes and their propagation paths whose contribution exceeds a preset significance threshold, forming a set of high-impact anomaly propagation paths.

10. A cloud computing-based method for tracing the quality of pig feed production according to claim 9, characterized in that, In the virtual intervention matrix, the rows represent upstream process nodes, the columns represent downstream process nodes, and the matrix elements are the normalized propagation strengths of the corresponding causal edges.

Images