Method for the production of dairy products for lactose intolerant people based on enzymatic technology

Abstract

The present application relates to a method for producing dairy products for lactose-intolerant people based on enzymatic technology, and in particular to the field of intelligent industrial control of food biological manufacturing. By collecting multi-source historical data and extracting time sequence characteristics, combining deep long short-term memory network and multi-head attention mechanism, precise prediction of lactose residue and flavor is realized. At the same time, the lactose concentration is detected in real time by using near-infrared spectrum, and the detection delay is compressed to a very short time. The compound control law fuses feedback, feedforward, disturbance compensation and flavor constraint, dynamically adjusts the enzyme dosage, and automatically switches the enzymolysis tank temperature mode according to the filling downtime, effectively inhibits excessive hydrolysis, significantly improves the flavor stability under the premise of ensuring the lactose standard, reduces the waste of raw materials and energy consumption, and realizes the intelligent and high-quality continuous production of low-lactose dairy products.

Description

Technical Field

[0001] This invention relates to the field of intelligent industrial control in food biomanufacturing, and more specifically, to a method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology. Background Technology

[0002] Low-lactose dairy products are functional dairy products developed for people with lactose intolerance. Their core production process involves adding exogenous lactase to hydrolyze lactose in milk into glucose and galactose, resulting in a final product with lactose residue below the national standard of 0.5 grams per 100 milliliters. Currently, large-scale manufacturers generally use continuous production lines for the industrial production of low-lactose dairy products. These lines typically include a raw milk receiving and standardized pretreatment unit, enzymatic hydrolysis reactors, an online sterilization unit, and a downstream aseptic filling unit. In actual production, the enzymatic hydrolysis section faces continuous influence from multiple sources of dynamic disturbances. These internal disturbances manifest as seasonal or batch-specific fluctuations in the natural lactose content of each batch of raw milk, influenced by factors such as the dairy cow's feeding cycle and feed formulation. These fluctuations can reach significant levels. ±2%; Due to differences in strain culture conditions and purification processes, the enzyme activity of commercial lactase preparations can fluctuate by ±15% between different batches; Downstream filling units may experience random fluctuations in material residence time within the upstream enzymatic hydrolysis tank due to repacking operations, equipment disinfection, or temporary shutdowns; External disturbances manifest as seasonal changes in workshop ambient temperature and deviations in enzymatic hydrolysis reaction temperature caused by fluctuations in steam pipeline pressure, with a deviation range of approximately ±1 degree Celsius; The combined effect of the above disturbance factors, especially under the combined conditions of high raw milk lactose content, temporary filling shutdowns, and low enzyme activity, can easily lead to the residual lactose at the end of enzymatic hydrolysis exceeding the control limit, directly affecting product safety and its suitability for lactose-intolerant individuals.

[0003] In existing technologies, the production control of low-lactose dairy products generally adopts a single-variable proportional-integral-derivative (PID) control strategy based on a fixed formula. This control strategy relies on a static enzymatic hydrolysis kinetic model established under laboratory conditions, pre-setting a constant enzyme dosage, a fixed reaction temperature of 42 degrees Celsius, and a fixed reaction time of 30 minutes. The reaction temperature is only adjusted in a closed loop using a PID controller. This approach has two significant drawbacks: First, the control system is limited to the equipment-level control of a single enzymatic hydrolysis tank, failing to incorporate rapid detection data of raw milk lactose content, batch calibration data of enzyme activity, and real-time status signals of downstream filling cycles into a unified control loop. This results in the system being unable to adjust based on the characteristics of upstream raw materials. Firstly, feedforward compensation based on changes and downstream process demands is essentially open-loop control. Secondly, online lactose detection relies on near-infrared spectroscopy, but the measurement model for this technology in emulsion systems is not dynamically calibrated. The scattering effects of bubbles, fat globules, and protein particles cause significant lag in the measured values, with lag times exceeding fifteen minutes, making it impossible to provide effective feedback signals for real-time control. Therefore, existing technologies struggle to accurately control the residual lactose at the end of enzymatic hydrolysis under complex conditions of multi-source dynamic disturbances. This often results in lactose residue exceeding the national standard limit or excessive enzyme addition leading to glucose exceeding the standard and flavor deterioration, thus hindering the consistency of low-lactose dairy product quality and production stability. Summary of the Invention

[0004] This invention addresses the technical problems existing in the prior art by providing a method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology, thereby resolving the issues raised in the background section.

[0005] The technical solution of this invention to solve the above-mentioned technical problems is as follows: specifically, it includes the following steps:

[0006] Step S1: Collect full historical data from three consecutive years of large-scale low-lactose dairy production lines. The full historical data includes batch rapid test values ​​of raw milk lactose and protein content, batch activity calibration values ​​of commercial lactase, temperature and pH time-series curves of enzymatic hydrolysis tanks, workshop ambient temperature records, material residence time records in tanks, raw spectral matrix output by online near-infrared spectrometer, and offline test values ​​and flavor scores of finished product lactose residue. The data is then aligned with timestamps by production batch. The coefficient of variation of raw milk lactose content for each batch is calculated, the exponential decay curve of enzyme activity with storage time is fitted, and a sliding window standard deviation analysis is performed on the residence time series composed of residence time records to obtain a structured time-series feature vector set.

[0007] Step S2: Using the temporal feature vector set as input, construct a deep long short-term memory network containing an input gate, a forget gate, an output gate, and cell states. Introduce a multi-head attention mechanism after the hidden layer of the deep long short-term memory network and train it with a joint loss function containing a mean squared error term, a cross-entropy loss term, and an L2 regularization term. Output the predicted lactose residue and flavor score at the target time. The forget gate weight matrix incorporates a modulation factor for the ambient temperature deviation, which is calculated based on the difference between the workshop ambient temperature record and the standard temperature of 42 degrees Celsius to achieve adaptivity.

[0008] Step S3: Perform discrete wavelet decomposition on the original spectral matrix acquired online, filter out high-frequency noise and scattering components, and reconstruct the denoised spectrum. Use the absorbance value of the denoised spectrum as the independent variable and the offline test value of the lactose residue in the finished product as the dependent variable to establish a partial least squares regression model, thereby obtaining the real-time lactose detection value and compressing the lactose concentration detection delay to within 30 seconds.

[0009] Step S4: The real-time lactose detection value obtained in step S3 is used as a feedback signal, the lactose residue prediction value output in step S2 is used as a feedforward compensation benchmark, the start-stop cycle signal of the downstream filling unit is used as a feedforward disturbance input, and the flavor score prediction value output in step S2 is used as a flavor constraint index. All of these are connected to the program controller. The program controller executes a composite control law to calculate the real-time adjustment increment of the enzyme dosage, and automatically triggers the cooling standby mode and heating recovery mode of the enzymatic hydrolysis tank according to the cumulative downtime of filling.

[0010] In a preferred embodiment, in step S1, before aligning the timestamps according to the production batch, missing value filling is performed: missing values ​​in the sensor channel caused by communication interruption or sensor saturation are filled using cubic spline interpolation, with the interpolation node spacing being twice the original sampling interval;

[0011] After completing the missing value filling, the specific operation of aligning the timestamps by production batch is as follows: using the temperature data in the temperature and pH time series curve of the enzymatic hydrolysis reactor as a reference, the data of all other sensor channels are shifted and corrected respectively. For each sensor channel, the cross-correlation function between the data sequence of the sensor channel and the temperature data sequence is calculated, the time offset that maximizes the cross-correlation function value is found, and the data of each sensor channel is shifted as a whole by a negative time offset sampling point.

[0012] After alignment, for each production batch, a multidimensional time series is constructed using the temperature data from the temperature-pH time series curve of the enzymatic hydrolysis reactor, the pH data from the temperature-pH time series curve of the enzymatic hydrolysis reactor, the first principal component score obtained by principal component analysis of the original spectral matrix output by the online near-infrared spectrometer, and the residence time record of the material in the enzymatic hydrolysis reactor. The local anomaly factor algorithm is used to calculate the anomaly score at each time point. Time points with anomaly scores greater than 1.5 are considered as anomaly segments and are removed. The length of the continuous valid interval formed after removal must not be less than 80% of the original length; otherwise, the entire batch of data is marked as invalid and excluded.

[0013] In a preferred embodiment, the specific process of calculating the coefficient of variation of lactose content in each batch of raw milk, fitting the exponential decay curve of enzyme activity with storage time, and performing sliding window standard deviation analysis on the residence time series composed of residence time records is as follows:

[0014] Calculate the coefficient of variation of raw milk lactose content in each batch. The coefficient of variation is equal to the weighted standard deviation divided by the weighted mean. The weighted standard deviation is the square root of the weighted average of the squared deviations of raw milk lactose content from the weighted mean at each sampling point. The weighted mean is the weighted average of raw milk lactose content at each sampling point.

[0015] The Levenberg-Marquardt algorithm was used to fit the exponential decay curve of the batch activity calibration value of commercial lactase as the storage days. After fitting, the initial enzyme activity and decay constant were output as two independent features.

[0016] For the residence time series consisting of records of the residence time of materials in the tank, a sliding window is used for standard deviation analysis. Specifically, the scaling value of the median absolute deviation within the window is calculated. The window length is fixed at five consecutive sampling points, and the window sliding step size is one sampling point. The maximum value of the scaling value of the median absolute deviation of all windows is taken as the result of the sliding window standard deviation analysis. The scaling value of the median absolute deviation is equal to 1.4826 multiplied by the median of the absolute deviation of the residence time within the window from the median of the median of the window. The coefficient 1.4826 makes the scaling value of the median absolute deviation a consistent estimator of the standard deviation. When calculating the median of the absolute deviation of the residence time within the window from the median of the window, the median of the residence time within the window is first calculated, then the absolute value of the difference between each residence time within the window and the median is calculated, and finally these absolute differences are sorted and the median is taken.

[0017] In a preferred embodiment, the process of obtaining a structured temporal feature vector set specifically includes:

[0018] After feature extraction, the coefficient of variation, initial enzyme activity, decay constant, sliding window standard deviation analysis results, arithmetic mean of protein content in the batch rapid test values ​​of raw milk lactose and protein content, and arithmetic mean of workshop ambient temperature records for each batch are combined into a structured feature vector. The structured feature vectors of all batches are arranged in chronological order, and each feature dimension is standardized by subtracting the global mean and dividing by the global standard deviation to obtain the final time-series feature vector set.

[0019] In a preferred embodiment, in step S2, the temporal feature vector set obtained in step S1 is first organized into a three-dimensional tensor in chronological order, wherein the first dimension is the batch number, the second dimension is the time step number within each batch, and the third dimension is the feature dimension with a size of six. The temporal feature vector of each time step is the feature vector in the temporal feature vector set. The three-dimensional tensor is subjected to layer normalization processing, and the processed three-dimensional tensor is input into a deep long short-term memory network.

[0020] The deep long short-term memory network is constructed as a two-layer stacked structure, with each layer containing 128 hidden units. At each time step, the deep long short-term memory network outputs hidden states and cell states, where the hidden states are 128-dimensional vectors. The hidden states output by the deep long short-term memory network at all time steps constitute a hidden state sequence. The network output layer contains two independent branch structures. One branch is a fully connected layer used to output the predicted lactose residue value at the target time step. The other branch connects the fully connected layer and the Softmax activation function in sequence to output the predicted flavor score value.

[0021] During network training, a joint loss function is used, which consists of the sum of the mean squared error term, the cross-entropy loss term, and the L2 regularization term.

[0022] In deep long short-term memory networks, the calculation process of the forget gate activation value introduces a modulation factor of the ambient temperature deviation. The calculation process is as follows: the difference between the workshop ambient temperature record at the current time step and the standard temperature of 42 degrees Celsius is obtained as the ambient temperature deviation. The ambient temperature deviation is compressed by the hyperbolic tangent function, multiplied by the learnable modulation factor, and then multiplied element-wise with the learnable projection vector and superimposed on the linear transformation result of the forget gate.

[0023] In a preferred embodiment, the specific operation of introducing a multi-head attention mechanism after the hidden layer of the deep long short-term memory network is as follows:

[0024] First, the hidden state sequence output by the last deep long short-term memory network is used as the input of the multi-head attention mechanism, and the number of attention heads is set to eight. Each attention head independently performs a linear transformation on the hidden state sequence to generate a query matrix, a key matrix, and a value matrix, respectively. The feature dimensions of the query matrix and the key matrix are set to sixteen, and the feature dimension of the value matrix is ​​set to sixteen.

[0025] Subsequently, each attention head calculates the product of the query matrix and the transpose of the key matrix, divides the product by the square root of the feature dimension of the key matrix for scaling, and then performs Softmax normalization on the scaled result to obtain the attention weight matrix. The attention weight matrix is ​​then multiplied by the value matrix to calculate the weighted output representation of each attention head.

[0026] Finally, the outputs of all eight attention heads are concatenated along the feature dimension, and the total dimension of the concatenation is mapped back to the hidden layer dimension through a linear projection layer, where the hidden layer dimension is 128. The mapping result is further processed by a non-linear activation function to obtain the final attention representation. After the attention representation is calculated, it is input into the network output layer to obtain the lactose residue prediction value and flavor score prediction value.

[0027] In a preferred embodiment, the specific operation of performing discrete wavelet decomposition on the online acquired original spectral matrix in step S3 is as follows:

[0028] Discrete wavelet decomposition uses the Dobessi type 5 mother wavelet with four decomposition levels. Each level generates a set of approximation coefficients and a set of detail coefficients. Adaptive soft thresholding is applied to the detail coefficients to filter out high-frequency noise and scattering components. The calculation process for the adaptive soft thresholding is as follows:

[0029] First, obtain the detail coefficients of the current layer. Then, take the absolute value of the detail coefficients and calculate the median. Divide the median by 0.6745 to obtain the estimated noise standard deviation of the current layer. Then, multiply the estimated noise standard deviation by the square root of the number of wavelengths under the natural logarithm to calculate the threshold of the current layer.

[0030] The specific operation of soft thresholding is as follows: iterate through each detail coefficient of the current layer. If the absolute value of the detail coefficient is less than or equal to the calculated threshold of the current layer, then set the value of the detail coefficient to zero. If the absolute value of the detail coefficient is greater than the threshold of the current layer, then subtract the threshold of the current layer from the detail coefficient and keep its original positive or negative sign unchanged.

[0031] The inverse discrete wavelet transform is performed on all the layer detail coefficients after thresholding and the approximation coefficients of the fourth layer to reconstruct a denoised spectrum. The same discrete wavelet decomposition and adaptive soft thresholding are then performed on each column of the original spectral matrix at all time points. The denoised spectra at all time points together form the denoised spectral matrix.

[0032] In a preferred embodiment, the specific process of establishing a partial least squares regression model with the absorbance value of the denoised spectrum as the independent variable and the offline test value of the lactose residue in the finished product as the dependent variable is as follows:

[0033] Using the absorbance value of the denoised spectral matrix as the independent variable and the offline test value of the lactose residue in the finished product accumulated in step S1 as the dependent variable, the potential components of the independent and dependent variables are extracted by nonlinear iterative partial least squares algorithm, and the number of potential components is determined by cross-validation.

[0034] After extracting the latent components, the weight matrix, loading matrix, and dependent variable loading vector are calculated sequentially. Then, the regression coefficient vector is solved through matrix operations. The solution process is as follows:

[0035] First, calculate the inverse matrix of the product of the weight matrix and the loading matrix. Then, perform matrix multiplication of the inverse matrix with the dependent variable loading vector. The result is the regression coefficient vector, whose dimension is consistent with the number of spectral wavelength points. For each new original spectral vector acquired online, first perform the same discrete wavelet decomposition and adaptive soft thresholding as in step S3 to obtain the corresponding denoised spectrum. Then, perform an inner product operation on the denoised spectrum and the regression coefficient vector. The inner product result is the real-time lactose detection value.

[0036] In a preferred embodiment, in step S4, the composite control law is composed of a superposition of a feedback term, a feedforward compensation term, a feedforward disturbance term, and a flavor constraint term.

[0037] The calculation process of the feedback term is as follows: First, calculate the difference between the real-time lactose detection value and the preset lactose residue target value as the feedback error. Then, take the absolute value of the feedback error and perform a power nonlinear operation. Multiply the result of the power nonlinear operation by the positive or negative sign of the feedback error, and then multiply it by the feedback gain coefficient to obtain the value of the feedback term.

[0038] The calculation process of the feedforward compensation term is as follows: First, calculate the difference between the predicted lactose residue value and the preset lactose residue target value as the feedforward error. Then, perform an integral operation on the feedforward error from time zero to the current time. Multiply the integral result by the feedforward integral gain coefficient to obtain the value of the feedforward compensation term.

[0039] The calculation process for the feedforward disturbance term is as follows: multiply the start-stop cycle signal of the downstream filling unit by the disturbance feedforward gain coefficient to obtain the value of the feedforward disturbance term;

[0040] The calculation process of the flavor constraint term is as follows: First, calculate the difference between the preset flavor threshold and the predicted flavor score. If the difference is positive, retain the difference; if the difference is negative, take zero. Multiply the retained difference by the flavor constraint weight coefficient to obtain the value of the flavor constraint term.

[0041] The real-time adjustment increment of enzyme dosage is obtained by adding the feedback term, feedforward compensation term, feedforward perturbation term, and flavor constraint term.

[0042] In a preferred embodiment, the process of automatically triggering the cooling standby mode and the heating recovery mode of the enzymatic hydrolysis tank based on the cumulative downtime of filling specifically includes:

[0043] The program controller monitors the cumulative downtime of filling in real time. When the cumulative downtime exceeds the preset downtime threshold, the program controller automatically triggers the cooling standby mode and starts the cooling system to reduce the temperature inside the enzymatic hydrolysis tank to ten degrees Celsius to inhibit enzyme activity.

[0044] When filling resumes and the cumulative downtime of filling is less than the preset restart temperature threshold, the program controller automatically triggers the temperature recovery mode, starts the heating system to raise the temperature inside the enzymatic hydrolysis tank to 42 degrees Celsius, and maintains it until the enzymatic hydrolysis is completed. A hysteresis interval is set between the cooling trigger point and the heating trigger point. When the cumulative downtime of filling is between the preset downtime threshold and the preset restart temperature threshold, the program controller maintains the current temperature mode unchanged.

[0045] The beneficial effects of this invention are as follows: by collecting multi-source historical data and extracting temporal features, and combining deep long short-term memory networks with multi-head attention mechanisms, accurate prediction of lactose residue and flavor is achieved; at the same time, near-infrared spectroscopy is used to detect lactose concentration in real time, compressing the detection delay to an extremely short time; the composite control law integrates feedback, feedforward, disturbance compensation and flavor constraints, dynamically adjusts the enzyme dosage, and automatically switches the temperature mode of the enzymatic hydrolysis tank according to the filling downtime, effectively inhibiting excessive hydrolysis, significantly improving flavor stability while ensuring lactose compliance, reducing raw material and energy waste, and realizing intelligent, high-quality continuous production of low-lactose dairy products. Attached Figure Description

[0046] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0047] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0048] In the description of this application, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include one or more features. In the description of this application, "multiple" means two or more, unless otherwise explicitly specified.

[0049] In the description of this application, the term "for example" is used to mean "used as an example, illustration, or description." Any embodiment described as "for example" in this application is not necessarily to be construed as being more preferred or advantageous than other embodiments. The following description is provided to enable any person skilled in the art to make and use the invention. Details are set forth in the following description for purposes of explanation. It should be understood that those skilled in the art will recognize that the invention can be made without using these specific details. In other instances, well-known structures and processes will not be described in detail to avoid obscuring the description of the invention with unnecessary detail. Therefore, the invention is not intended to be limited to the embodiments shown, but is consistent with the broadest scope of the principles and features disclosed in this application.

[0050] Example 1: This example provides the following... Figure 1 The method shown is a production control method for dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology, which specifically includes the following steps:

[0051] Step S1: Collect full historical data from three consecutive years of large-scale low-lactose dairy production lines. The full historical data includes batch rapid test values ​​of raw milk lactose and protein content, batch activity calibration values ​​of commercial lactase, temperature and pH time-series curves of enzymatic hydrolysis tanks, workshop ambient temperature records, material residence time records in tanks, raw spectral matrix output by online near-infrared spectrometer, and offline test values ​​and flavor scores of finished product lactose residue. The data is then time-stamped by production batch to remove abnormal segments. The coefficient of variation of raw milk lactose content for each batch is calculated, the exponential decay curve of enzyme activity with storage time is fitted, and the standard deviation analysis of residence time series composed of residence time records is performed to obtain a structured time-series feature vector set.

[0052] Step S2: Using the temporal feature vector set as input, construct a deep long short-term memory network containing an input gate, a forget gate, an output gate, and cell states. Introduce a multi-head attention mechanism after the hidden layer of the deep long short-term memory network and train it with a joint loss function containing a mean squared error term, a cross-entropy loss term, and an L2 regularization term. Output the predicted lactose residue and flavor score at the target time. The forget gate weight matrix incorporates a modulation factor for the ambient temperature deviation, which is calculated based on the difference between the workshop ambient temperature record and the standard temperature of 42 degrees Celsius to achieve adaptivity.

[0053] Step S3: Perform discrete wavelet decomposition on the original spectral matrix acquired online, filter out high-frequency noise and scattering components, and reconstruct the denoised spectrum. Use the absorbance value of the denoised spectrum as the independent variable and the offline test value of the lactose residue in the finished product as the dependent variable to establish a partial least squares regression model, thereby obtaining the real-time lactose detection value and compressing the lactose concentration detection delay to within 30 seconds.

[0054] Step S4: The real-time lactose detection value obtained in step S3 is used as a feedback signal, the lactose residue prediction value output in step S2 is used as a feedforward compensation benchmark, the start / stop cycle signal of the downstream filling unit is used as a feedforward disturbance input, and the flavor score prediction value output in step S2 is used as a flavor constraint index. These are all connected to the program controller. The program controller executes a composite control law to calculate the real-time adjustment increment of the enzyme dosage, and automatically triggers the cooling standby mode and heating recovery mode of the enzymatic hydrolysis tank according to the cumulative downtime of filling. When the flavor score prediction value is lower than the preset flavor threshold, the program controller increases the weight coefficient of the feedforward compensation term in the composite control law to reduce the enzymatic hydrolysis intensity and prevent flavor deterioration.

[0055] In this embodiment, it is specifically necessary to explain that in step S1, before aligning the timestamps according to the production batch, missing value filling is performed: for missing values ​​in the sensor channel caused by communication interruption or sensor saturation, cubic spline interpolation is used to fill them, and the interpolation node spacing is twice the original sampling interval to avoid overfitting; when using cubic spline interpolation to fill, the interpolation node spacing is set to twice the original sampling interval, the purpose of which is to avoid overfitting caused by too dense nodes while ensuring the smoothness of interpolation. The cubic spline interpolation method obtains piecewise cubic polynomial coefficients by solving a tridiagonal system of equations, ensuring that the interpolation curve has continuous second derivatives at the nodes, thereby approximating the changing trend of the real signal;

[0056] After completing the missing value filling, the specific operation for timestamp alignment by production batch is as follows: Using the temperature data in the temperature and pH time series curves of the enzymatic hydrolysis reactor as a reference, the data of all other sensor channels are shifted and corrected. For each sensor channel, the cross-correlation function between the sensor channel's data sequence and the temperature data sequence is calculated. The time offset that maximizes the cross-correlation function value is found, and the data of each sensor channel is shifted as a whole by a negative time offset of sampling points, thereby eliminating the phase deviation caused by sensor clock asynchrony. When calculating the cross-correlation function, the corresponding elements of the sensor channel's data sequence and the temperature data sequence are multiplied one by one and summed. The length of the two sequences minus one is the maximum offset. The process involves iterating through all possible time offsets and calculating the corresponding sum of products. The offset corresponding to the maximum sum of products is the time offset. The cross-correlation function measures the similarity between two sequences at different time offsets. During calculation, the corresponding elements of the data sequence from the sensor channel and the temperature data sequence are multiplied one by one and summed. The resulting sum of products is the cross-correlation function value. The process involves iterating through all possible time offsets (ranging from negative sequence length plus one to positive sequence length minus one) to find the offset that maximizes the sum of products. The absolute value of the offset reflects the phase difference between the two sequences, and its sign indicates which sequence is ahead or behind. By shifting the entire sensor channel data by a negative offset by a number of sampling points, precise alignment with the temperature sequence can be achieved.

[0057] After alignment, for each production batch, a multidimensional time series is constructed using the temperature data from the temperature-pH time series curve of the enzymatic hydrolysis reactor, the pH data from the temperature-pH time series curve of the enzymatic hydrolysis reactor, the first principal component score obtained by principal component analysis of the original spectral matrix output from the online near-infrared spectrometer, and the residence time record of the material in the enzymatic hydrolysis reactor. The local anomaly factor algorithm is used to calculate the anomaly score at each time point. The local anomaly factor is based on the ratio of the local reachability density of a sample point to the local reachability density of other points in its neighborhood. The process of calculating the local reachability density is as follows: first, determine the neighborhood radius; then, calculate the reachability distance of each sample point to all points in its neighborhood, and thus obtain the local reachability density. Time points with anomaly scores greater than 1.5 are considered anomaly segments and are removed. The length of the continuous valid interval formed after removal must not be less than 80% of the original length; otherwise, the entire batch of data is marked as invalid and excluded. The k-distance-based method is adopted, whereby for each sample point, its Euclidean distance to all other points is calculated and sorted in ascending order, and the k-th nearest distance is taken as the k-distance for that point, where k is taken as an empirical value of 5. The reachability distance is defined as the larger of the actual Euclidean distance between the sample point and another point in its neighborhood and the k-distance of that other point. The local reachability density is the reciprocal of the sum of the reachability distances between the sample point and all points in its neighborhood, multiplied by the number of points in the neighborhood. The local anomaly factor score is the average ratio of the local reachability density of the sample point to the local reachability density of all points in its neighborhood. Time points with a score greater than 1.5 are considered as anomaly segments, and the threshold is determined based on statistical analysis of historical normal production data, which can keep the false alarm rate within 5%. The length of the continuous valid interval formed after elimination must not be less than 80% of the original length. If it is less than 80%, the entire batch of data is marked as invalid and excluded to ensure the integrity of the data used for subsequent modeling.

[0058] The specific process for calculating the coefficient of variation of lactose content in each batch of raw milk, fitting the exponential decay curve of enzyme activity with storage time, and performing sliding window standard deviation analysis on the residence time series composed of residence time records is as follows:

[0059] The coefficient of variation (COP) of raw milk lactose content within each batch was calculated, with weights determined by the time interval between adjacent sampling points, assigning higher weights to more recent data. The COP equals the weighted standard deviation divided by the weighted mean. The weighted standard deviation is the square root of the weighted average of the squared deviations of the raw milk lactose content from the weighted mean at each sampling point, and the weighted mean is the weighted average of the raw milk lactose content at each sampling point, used to quantify the short-term fluctuation intensity of raw milk lactose. The formula for calculating the weighted standard deviation is as follows: first, calculate the squared deviation of the raw milk lactose content from the weighted mean at each sampling point, then multiply by the corresponding weight, sum the results, divide by the sum of the weights, and finally take the square root. The weighted mean is the sum of the raw milk lactose content multiplied by the corresponding weight at each sampling point, then divided by the sum of the weights. The weights are determined by the reciprocal of the time interval between adjacent sampling points; the shorter the time interval, the greater the weight, thus highlighting the contribution of recent data. The weighted COP can sensitively reflect the short-term fluctuation intensity of raw milk lactose, providing a preliminary reference for subsequent enzyme dosage adjustment.

[0060] To investigate the variation of commercial lactase batch activity calibration values ​​with storage days, the Levenberg-Marquardt algorithm was used to fit the exponential decay curve. The initial enzyme activity and decay constant were output as two independent features, representing the initial potency and decay rate of the enzyme, respectively. The Levenberg-Marquardt algorithm was used for nonlinear least-squares fitting, with the initial values ​​set as follows: the initial enzyme activity was the maximum value from the calibration data, and the decay constant was 0.01 day⁻¹. The algorithm iteration terminated when the parameter change between two consecutive iterations was less than 1e-6 or the number of iterations reached 100. The initial enzyme activity and decay constant output after fitting represent the enzyme preparation's activity level at the time of manufacture and the rate of activity decrease with storage time, respectively. These two parameters directly determine the calculation basis for enzyme dosage.

[0061] For the residence time series composed of records of material residence time in the tank, a sliding window standard deviation analysis is performed. Specifically, the scaled value of the median absolute deviation within the window is calculated. The window length is fixed at five consecutive sampling points, and the window sliding step size is one sampling point. The maximum scaled value of the median absolute deviation across all windows is taken as the standard deviation analysis result of the sliding window, used to capture extreme residence time abrupt changes caused by filling, repackaging, or disinfection operations. The scaled value of the median absolute deviation is equal to 1.4826 multiplied by the median of the absolute deviations between the residence time within the window and the median within the window. The coefficient 1.4826 makes the scaled value of the median absolute deviation a consistent estimator of the standard deviation. When calculating the median of the absolute deviations between the residence time within the window and the median within the window, the median residence time within the window is first calculated, then the absolute value of the difference between each residence time within the window and the median is calculated separately. Finally, these absolute differences are sorted and the result is taken. The median is used as a robust estimate of the standard deviation by employing a scaled value of the median absolute deviation. This is because outliers may exist in the dwell time series due to repacking operations, and the median absolute deviation is insensitive to outliers. The coefficient 1.4826 is derived from the asymptotic relationship between the median absolute deviation and the standard deviation under a normal distribution. That is, when the data follows a normal distribution, the expected value of the median absolute deviation is equal to the standard deviation multiplied by 0.6745. Taking its reciprocal, 1.4826, makes the scaled value an unbiased estimate of the standard deviation. The window length is fixed at five consecutive sampling points, balancing sensitivity to local fluctuations with computational stability. The sliding step size is one sampling point to ensure that no possible abrupt changes are missed. The maximum value of the scaled value of the median absolute deviation of all windows is used as the result of the sliding window standard deviation analysis. This can capture extreme dwell time abrupt changes caused by filling, repacking, or disinfection operations, and the feature has an important warning effect on subsequent control decisions.

[0062] The specific process of obtaining a structured time-series feature vector set is as follows:

[0063] After feature extraction, the coefficient of variation, initial enzyme activity, decay constant, sliding window standard deviation analysis results, arithmetic mean of protein content in the batch rapid test values ​​of raw milk lactose and protein content, and arithmetic mean of workshop ambient temperature records for each batch are combined into a structured feature vector. The structured feature vectors of all batches are arranged in chronological order, and each feature dimension is standardized by subtracting the global mean and dividing by the global standard deviation to obtain the final time-series feature vector set. The six features extracted from each batch (coefficient of variation, initial enzyme activity, decay constant, sliding window standard deviation analysis results, arithmetic mean of protein content, and arithmetic mean of workshop ambient temperature) are arranged in a fixed order. The features are combined into a one-dimensional vector. The feature vectors of all batches are arranged in chronological order of production time to form a two-dimensional matrix, where rows represent batches and columns represent feature dimensions. During standardization, for each column of features, the arithmetic mean (global mean) and sample standard deviation (global standard deviation) of all batches in the column are calculated first. Then, for each element of the column, the global mean is subtracted and the global standard deviation is divided. The standardized features have a mean of 0 and a variance of 1, which eliminates the adverse effects of differences in feature dimensions and scales on the subsequent training of deep learning models, while retaining the relative differences between batches. The resulting set of temporal feature vectors is the input data for the deep long short-term memory network in step S2.

[0064] In this embodiment, it is specifically necessary to explain step S2, where the temporal feature vector set obtained in step S1 is first organized into a three-dimensional tensor in chronological order. The first dimension is the batch number, the second dimension is the number of time steps within each batch, and the third dimension is the feature dimension with a size of six. The temporal feature vector of each time step is the feature vector in the temporal feature vector set. The division of time steps is based on the temporal continuity of production batches. Several consecutively produced batches are divided into a time window, and the number of batches within each time window is the size of the second dimension. If the number of batches contained in a time window is insufficient, the dimension is maintained by truncation or zero padding. The default length of the time window is set to ten batches. When the number of consecutive batches is less than ten... At any given time step, a virtual batch vector with a value of zero is added to the end of the sequence to complete the dimension. When the number of consecutive batches exceeds ten, only the ten most recent batches are retained, and the data of earlier batches are discarded to ensure that the dimension of all time windows is consistent. To ensure gradient stability, layer normalization is performed on the three-dimensional tensor. That is, the mean and variance of the temporal feature vector at each time step are calculated independently and then standardized. The processed three-dimensional tensor is then input into the deep long short-term memory network. The specific calculation of layer normalization is as follows: for the six-dimensional feature vector at each time step, the arithmetic mean and sample standard deviation of the six elements of the vector are calculated. Then, the mean is subtracted from each element and the result is divided by the standard deviation to obtain the normalized vector. This operation does not share statistics across time steps to preserve the differences between time steps.

[0065] The Deep Long Short-Term Memory (LSTM) network is constructed as a two-layer stacked structure, with each layer containing 128 hidden units. The two-layer stacking means that the hidden state sequence output by the first LSM layer is used as the input to the second LSM layer, which outputs the final hidden state sequence. No additional skip connections or residual structures are introduced between the two layers. At each time step, the Deep LSM network outputs both hidden states and cell states, where each hidden state is a 128-dimensional vector. The hidden states output by the Deep LSM network at all time steps constitute the hidden state sequence. The network output layer contains... The network employs two independent branch structures. One branch is a fully connected layer that outputs the predicted lactose residue value at the target time. The other branch connects the fully connected layer and a Softmax activation function to output the predicted flavor score, which is divided into five levels. The output dimension of the fully connected layer is either one-dimensional (predicted lactose residue value) or five-dimensional (unnormalized scores for each flavor score level). The Softmax function converts the five-dimensional scores into a probability distribution, with the sum of the probabilities being one. A joint loss function is used during network training, consisting of a mean squared error term, a cross-entropy loss term, and L2 regularization. The joint loss function is composed of several terms: a mean squared error term quantifies the deviation between the predicted lactose residue value and the offline test value of the finished product; a cross-entropy loss term quantifies the difference between the predicted probability distribution of the flavor score and the offline test value of the flavor score; and an L2 regularization term constrains all weight parameters of the network to prevent overfitting. The calculation process of the joint loss function is as follows: first, the arithmetic mean of the squares of the differences between the offline test value of the finished product's lactose residue and the predicted lactose residue value in all training samples is calculated to obtain the mean squared error term; then, the offline test value of the flavor score after probability distribution transformation is compared with the predicted probability... The arithmetic mean of the cross-entropy between distributions yields the cross-entropy loss term. This cross-entropy loss term is multiplied by a balance coefficient, added to the mean squared error term, and then added to the sum of squares of all network weight parameters multiplied by a regularization coefficient to obtain the final joint loss function. The balance coefficient adjusts the contribution weight of the flavor prediction task in the total loss, while the regularization coefficient suppresses model overfitting. The optimal value for the balance coefficient is 0.3, ensuring that the flavor prediction task contributes approximately 30% to the lactose prediction task. The optimal value for the regularization coefficient is 10 to the power of negative 5, at which the model exhibits the lowest generalization error on the validation set.

[0066] In deep long short-term memory networks, the calculation of the forget gate activation value incorporates a modulation factor based on the ambient temperature deviation. The calculation process is as follows: First, the hidden state from the previous time step and the temporal feature vector from the current time step are obtained. The hidden state from the previous time step is a 128-dimensional vector, and the temporal feature vector from the current time step is a 6-dimensional vector. Then, the linear transformation result of the forget gate is calculated. This involves using the forget gate weight matrix to perform a linear transformation on the concatenated vector of the hidden state from the previous time step and the temporal feature vector from the current time step, and adding a forget gate bias vector. The forget gate weight matrix has a dimension of 128 x 134, and the forget gate bias vector has a dimension of 128. Next, the difference between the current ambient temperature record and the standard temperature of 42 degrees Celsius is obtained, and this difference is taken as the ambient temperature deviation. Then, a hyperbolic tangent function compression operation is performed on the ambient temperature deviation, and the compressed result is multiplied by a learnable modulation factor. This learnable modulation factor is a scalar with a value ranging from zero to one, and an initial value of zero. First, the product is multiplied element-wise with a learnable projection vector, which is a 128-dimensional vector initially initialized with a normal distribution with a mean of zero and a standard deviation of 0.01. Then, the multiplication result is superimposed on the linear transformation result of the forget gate. Finally, the activation value of the forget gate is calculated, i.e., the superimposed linear transformation result is input into the Sigmoid function for mapping. The hyperbolic tangent function compresses the ambient temperature deviation to between -1 and +1, ensuring the modulation amplitude does not exceed the absolute value of the learnable modulation factor, thus avoiding gradient explosion caused by extreme temperature values. The learnable modulation factor is automatically updated through backpropagation during training, and its final value reflects the strength of the influence of ambient temperature on the forget gate. The design allows the ambient temperature deviation to dynamically adjust the degree of retention of historical information by the forget gate. When the ambient temperature deviation is positive, the forget gate tends to retain more historical information; when the ambient temperature deviation is negative, the forget gate tends to discard more historical information, thereby simulating the Arrhenius effect of temperature on reaction rate in enzymatic reactions.

[0067] The specific operation of introducing a multi-head attention mechanism after the hidden layers of a deep long short-term memory network is as follows:

[0068] First, the hidden state sequence output from the last deep long short-term memory network is used as the input to the multi-head attention mechanism, with eight attention heads. Each attention head independently performs a linear transformation on the hidden state sequence, generating a query matrix, a key matrix, and a value matrix. The feature dimensions of the query matrix and key matrix are set to sixteen, and the feature dimension of the value matrix is ​​also set to sixteen. The linear transformation process is as follows: each attention head holds an independent set of weight parameters, including a query weight matrix, a key weight matrix, and a value weight matrix, all of which have a dimension of 128 x 16. Each attention head performs matrix multiplication on the input hidden state sequence and the corresponding weight matrix to generate the query matrix, key matrix, and value matrix, respectively.

[0069] Subsequently, each attention head calculates the product of the query matrix and the transpose of the key matrix. The product is then scaled by dividing by the square root of the feature dimension of the key matrix. Softmax normalization is then performed on the scaled result to obtain the attention weight matrix. The attention weight matrix is ​​then multiplied by the value matrix to calculate the weighted output representation of each attention head. The purpose of dividing by the square root of the feature dimension is to prevent the softmax gradient from vanishing due to an excessively large inner product. This scaling factor originates from the standard design of the Transformer architecture.

[0070] Finally, the outputs of all eight attention heads are concatenated along the feature dimension, and the total dimension of the concatenation is mapped back to the hidden layer dimension through a linear projection layer, where the hidden layer dimension is 128. The mapping result is further processed by a non-linear activation function to obtain the final attention representation. The non-linear activation function used is the ReLU function, which takes the maximum value between the input and zero to introduce sparsity and non-linearity. The weight matrix dimension of the linear projection layer is 128 by 128, and the bias vector dimension is 128. After calculating the attention representation, the attention representation is input into the network output layer to obtain the lactose residue prediction value and flavor score prediction value. The attention representation is a vector sequence with the same length as the hidden state sequence, but the vector at each time step incorporates global context information. The network output layer only takes the attention representation of the last time step as the input of the fully connected layer to output the prediction value at the target time. The design enables the network to dynamically focus on the information most important to the current prediction in the historical time steps. When the raw milk lactose fluctuates drastically, the attention weights automatically focus on recent fluctuation features and ignore earlier stable periods, thereby improving the accuracy of the lactose residue prediction value and flavor score prediction value.

[0071] In this embodiment, it is specifically necessary to explain the following steps in step S3: the discrete wavelet decomposition of the original online spectral matrix is ​​performed as follows:

[0072] Discrete wavelet decomposition employs the Dobercy type 5 mother wavelet with four decomposition levels. Each level generates a set of approximation coefficients and a set of detail coefficients. The Dobercy type 5 mother wavelet possesses orthogonality and compact support properties, making it suitable for processing sharp absorption peaks in near-infrared spectra. The four-level decomposition effectively separates high-frequency noise and scattering components while preserving low-frequency chemical information. Adaptive soft thresholding is used on the detail coefficients to filter out high-frequency noise and scattering components. The calculation process for the adaptive soft thresholding is as follows:

[0073] First, obtain the detail coefficients of the current layer. Then, calculate the median after taking the absolute value of the detail coefficients. Divide the median by 0.6745 to obtain the estimated standard deviation of the noise in the current layer. 0.6745 is the reciprocal of the ratio of the absolute deviation of the median to the standard deviation under a standard normal distribution. This estimation method is not sensitive to outliers and is suitable for potential spike interference in spectral data. Subsequently, multiply the estimated standard deviation of the noise by the square root of the number of wavelengths under two natural logarithms to calculate the threshold of the current layer, where the number of wavelengths is the wavelength dimension of the original spectral matrix. This threshold formula is derived from the minimax criterion and minimizes the mean square error after denoising under the assumption that the noise follows a Gaussian distribution.

[0074] The specific operation of soft thresholding is as follows: iterate through each detail coefficient of the current layer. If the absolute value of the detail coefficient is less than or equal to the calculated threshold of the current layer, then set the value of the detail coefficient to zero. If the absolute value of the detail coefficient is greater than the threshold of the current layer, then subtract the threshold of the current layer from the detail coefficient and keep its original positive or negative sign unchanged. Compared with hard thresholding, soft thresholding can avoid discontinuous jumps at the threshold, making the reconstructed spectrum smoother and reducing pseudo-Gibbs oscillations.

[0075] All layer detail coefficients after thresholding and the approximation coefficients of the fourth layer are subjected to inverse discrete wavelet transform to reconstruct a denoised spectrum. The same discrete wavelet decomposition and adaptive soft thresholding are then performed on each column of the original spectral matrix at all time points. The denoised spectra at all time points together form the denoised spectral matrix. The denoised spectral matrix retains the chemical absorption characteristics related to lactose content, while suppressing interference caused by particle scattering, baseline drift and electronic noise, providing a clean spectral signal for subsequent quantitative modeling.

[0076] The specific process of establishing a partial least squares regression model with the absorbance value of the denoised spectrum as the independent variable and the offline test value of the lactose residue in the finished product as the dependent variable is as follows:

[0077] Using the absorbance value of the denoised spectral matrix as the independent variable and the offline test value of the finished product lactose residue accumulated in step S1 as the dependent variable, a nonlinear iterative partial least squares algorithm is used to extract the latent components of the independent and dependent variables. The number of latent components is determined by cross-validation, and the value of this number is in the range of five to ten. Partial least squares regression is suitable for situations where the independent variables are highly correlated (such as adjacent wavelength points in the near-infrared spectrum) and the number of samples is small. By extracting latent components, it maximizes the covariance of the independent and dependent variables at the same time, thus overcoming the multicollinearity problem of ordinary least squares. The number of latent components is determined by leave-one-out cross-validation, and the number that minimizes the sum of squares of the predicted residuals is selected, which is usually between five and ten. This range covers most industrial near-infrared application scenarios.

[0078] After extracting the latent components, the weight matrix, loading matrix, and dependent variable loading vector are calculated sequentially. Then, the regression coefficient vector is solved through matrix operations. The solution process is as follows:

[0079] First, the inverse matrix of the product of the weight matrix and the loading matrix is ​​calculated. Then, the inverse matrix is ​​multiplied by the dependent variable loading vector. The result is the regression coefficient vector, whose dimension is consistent with the number of spectral wavelengths. This solution process is equivalent to mapping the original spectral space to the latent component space and then regressing. Each element in the regression coefficient vector represents the contribution weight of the absorbance at the corresponding wavelength point to the lactose content. Positive values ​​indicate positive correlation, and negative values ​​indicate interference or cancellation. For each new original spectral vector acquired online, the same discrete wavelet decomposition and adaptive soft thresholding processing as in step S3 are first performed to obtain the corresponding denoised spectrum. Then, the denoised spectrum and the regression coefficient vector are multiplied by an inner product operation. The inner product result is the real-time lactose detection value, ensuring that the entire calculation process is completed within 30 seconds after a single spectral acquisition. The inner product operation involves only one vector multiplication, with extremely low computational complexity. Combined with the pre-stored regression coefficient vector, prediction can be completed in microseconds. The entire process, from spectral acquisition to output detection value, including wavelet decomposition, thresholding, reconstruction, and multiplication, takes no more than 30 seconds, meeting the requirements of online real-time control.

[0080] In this embodiment, it is specifically noted that in step S4, the composite control law is composed of a feedback term, a feedforward compensation term, a feedforward disturbance term, and a flavor constraint term superimposed on each other; the preset lactose residue target value is the upper limit of lactose residue concentration set according to the product specifications; the preset lactose residue target value is usually 0.1 g per liter, which corresponds to the requirement in the low lactose dairy product industry standard that the lactose content is less than 0.5 g per 100 ml, while reserving a safety margin;

[0081] The calculation process of the feedback term is as follows: First, the difference between the real-time lactose detection value and the preset lactose residue target value is calculated as the feedback error. After taking the absolute value of the feedback error, a power-law nonlinear operation is performed. The exponent of the power-law nonlinear operation ranges from 0.1 to 0.9. The result of the power-law nonlinear operation is multiplied by the sign of the feedback error, and then multiplied by the feedback gain coefficient to obtain the value of the feedback term. The feedback gain coefficient is a positive real number used to adjust the response strength of the feedback loop, and its value ranges from 0.5 to 2.0. When the exponent of the power-law nonlinear operation is 0.5, the feedback term has a high gain for small errors and saturation characteristics for large errors, which can effectively suppress overshoot and accelerate the convergence speed. When the feedback gain coefficient is 1.0, the feedback term outputs about 1 unit adjustment when the error is 1 gram per liter, which can be adjusted online according to the actual response effect.

[0082] The calculation process of the feedforward compensation term is as follows: First, calculate the difference between the predicted lactose residue value and the preset lactose residue target value as the feedforward error. Then, perform an integral operation on the feedforward error from time zero to the current time. Multiply the integral result by the feedforward integral gain coefficient to obtain the value of the feedforward compensation term. The feedforward integral gain coefficient is a positive real number used to adjust the strength of the feedforward compensation, and its value ranges from 0.01 to 0.1. When the feedforward integral gain coefficient is 0.05, the cumulative error per minute can generate a compensation of about 5%, which can effectively eliminate the steady-state deviation caused by enzyme activity decay or raw material fluctuations.

[0083] The calculation process of the feedforward disturbance term is as follows: multiply the start and stop cycle signal of the downstream filling unit by the disturbance feedforward gain coefficient to obtain the value of the feedforward disturbance term; the disturbance feedforward gain coefficient is a positive real number used to adjust the intensity of disturbance compensation, and the value range is 0.5 to 2.0; when the disturbance feedforward gain coefficient is 1.0, the disturbance compensation generated at the moment of filling stop just offsets the accelerated effect of enzymatic reaction caused by flow interruption, thus avoiding excessive hydrolysis of lactose;

[0084] The calculation process for the flavor constraint term is as follows: First, calculate the difference between the preset flavor threshold and the predicted flavor score. If the difference is positive, retain the difference; if the difference is negative, set it to zero. Multiply the retained difference by the flavor constraint weight coefficient to obtain the value of the flavor constraint term. The flavor constraint weight coefficient is a non-negative real number with an initial value of 1. When the predicted flavor score is lower than the preset flavor threshold, the program controller increases the feedforward integral gain coefficient in the feedforward compensation term, and the flavor constraint weight coefficient automatically increases. The two work together to reduce the enzymatic hydrolysis intensity and prevent flavor degradation. The typical maximum value of the flavor constraint weight coefficient does not exceed 5. The increase of the flavor constraint weight coefficient adopts an exponential growth method, that is, the flavor constraint weight coefficient is equal to the initial value multiplied by the natural constant α times the flavor deviation raised to a power of 2, so that for every 0.5 points decrease in the flavor score, the constraint intensity doubles, quickly curbing excessive enzymatic hydrolysis.

[0085] The real-time adjustment increment of enzyme dosage is obtained by adding the feedback term, feedforward compensation term, feedforward perturbation term, and flavor constraint term.

[0086] The process of automatically triggering the cooling and warming recovery modes of the enzymatic hydrolysis tank based on the cumulative downtime during filling is as follows:

[0087] The program controller monitors the cumulative downtime of filling in real time. When the cumulative downtime exceeds the preset downtime threshold, the program controller automatically triggers the cooling standby mode and starts the cooling system to lower the temperature inside the enzymatic hydrolysis tank to 10 degrees Celsius to inhibit enzyme activity. In the cooling standby mode, the cooling system operates at maximum power, and it takes about 15 minutes for the temperature inside the tank to drop from 42 degrees Celsius to 10 degrees Celsius. After the temperature drops to 10 degrees Celsius, the enzyme activity drops to less than 5% of that at room temperature, which can effectively prevent lactose from continuing to hydrolyze and produce bitter peptides during long-term downtime.

[0088] When filling resumes and the cumulative downtime is less than the preset restart temperature threshold, the program controller automatically triggers the temperature recovery mode, starting the heating system to raise the temperature inside the enzymatic hydrolysis tank to 42 degrees Celsius and maintain it until the enzymatic hydrolysis is complete. The preset downtime threshold is 30 minutes, and the preset restart temperature threshold is 5 minutes. The 30-minute downtime threshold is set because even with low enzyme activity, lactose residue may still exceed the limit after 30 minutes, thus requiring forced cooling. The 5-minute restart temperature threshold is set because short-term downtime (such as repackaging operations) does not require... Cooling can be achieved by simply raising the temperature back up. A hysteresis interval is set between the cooling trigger point and the heating trigger point. When the cumulative downtime for filling is between the preset downtime threshold and the preset restart temperature threshold, the program controller maintains the current temperature mode to avoid frequent switching. The hysteresis interval is 30 minutes to 5 minutes. That is, when the downtime exceeds 30 minutes, the system enters the cooling mode. Even if the downtime is shortened to less than 30 minutes thereafter, as long as it is not less than 5 minutes, the cooling mode is still maintained to prevent repeated switching between cooling and heating due to fluctuations in downtime, thereby extending the equipment life and saving energy.

[0089] It should be noted that the descriptions of each embodiment in the above embodiments have different focuses. For parts that are not described in detail in a certain embodiment, please refer to the relevant descriptions in other embodiments.

[0090] Those skilled in the art will understand that embodiments of the present invention can be provided as methods, systems, or computer program products. Therefore, the present invention can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention can take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0091] This invention is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart illustrations. Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0092] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0093] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0094] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention.

[0095] Obviously, those skilled in the art can make various modifications and variations to this invention without departing from its spirit and scope. Therefore, if these modifications and variations fall within the scope of the claims of this invention and their equivalents, this invention also intends to include these modifications and variations.

Claims

1. A method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology, characterized in that, Specifically, the steps include the following: Step S1: Collect full historical data from three consecutive years of large-scale low-lactose dairy production lines. The full historical data includes batch rapid test values ​​of raw milk lactose and protein content, batch activity calibration values ​​of commercial lactase, temperature and pH time-series curves of enzymatic hydrolysis tanks, workshop ambient temperature records, material residence time records in tanks, raw spectral matrix output by online near-infrared spectrometer, and offline test values ​​and flavor scores of finished product lactose residue. The data is then aligned with timestamps by production batch. The coefficient of variation of raw milk lactose content for each batch is calculated, the exponential decay curve of enzyme activity with storage time is fitted, and a sliding window standard deviation analysis is performed on the residence time series composed of residence time records to obtain a structured time-series feature vector set. Step S2: Using the temporal feature vector set as input, construct a deep long short-term memory network containing an input gate, a forget gate, an output gate, and cell states. Introduce a multi-head attention mechanism after the hidden layer of the deep long short-term memory network and train it with a joint loss function containing a mean squared error term, a cross-entropy loss term, and an L2 regularization term. Output the predicted lactose residue and flavor score at the target time. The forget gate weight matrix incorporates a modulation factor for the ambient temperature deviation, which is calculated based on the difference between the workshop ambient temperature record and the standard temperature of 42 degrees Celsius to achieve adaptivity. Step S3: Perform discrete wavelet decomposition on the original spectral matrix acquired online, filter out high-frequency noise and scattering components, and reconstruct the denoised spectrum. Use the absorbance value of the denoised spectrum as the independent variable and the offline test value of the lactose residue in the finished product as the dependent variable to establish a partial least squares regression model, thereby obtaining the real-time lactose detection value and compressing the lactose concentration detection delay to within 30 seconds. Step S4: The real-time lactose detection value obtained in step S3 is used as a feedback signal, the lactose residue prediction value output in step S2 is used as a feedforward compensation benchmark, the start-stop cycle signal of the downstream filling unit is used as a feedforward disturbance input, and the flavor score prediction value output in step S2 is used as a flavor constraint index. All of these are connected to the program controller. The program controller executes a composite control law to calculate the real-time adjustment increment of the enzyme dosage, and automatically triggers the cooling standby mode and heating recovery mode of the enzymatic hydrolysis tank according to the cumulative downtime of filling.

2. The method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology according to claim 1, characterized in that: In step S1, before aligning the timestamps according to the production batch, missing values ​​are filled: missing values ​​in the sensor channel caused by communication interruption or sensor saturation are filled using cubic spline interpolation, with the interpolation node spacing being twice the original sampling interval; After completing the missing value filling, the specific operation of aligning the timestamps by production batch is as follows: using the temperature data in the temperature and pH time series curve of the enzymatic hydrolysis reactor as a reference, the data of all other sensor channels are shifted and corrected respectively. For each sensor channel, the cross-correlation function between the data sequence of the sensor channel and the temperature data sequence is calculated, the time offset that maximizes the cross-correlation function value is found, and the data of each sensor channel is shifted as a whole by a negative time offset sampling point. After alignment, for each production batch, the temperature data in the temperature and pH time series curve of the enzymatic hydrolysis reactor, the pH data in the temperature and pH time series curve of the enzymatic hydrolysis reactor, the first principal component score obtained by principal component analysis of the original spectral matrix output by the online near-infrared spectrometer, and the residence time record of the material in the enzymatic hydrolysis reactor are used to form a multidimensional time series. The local anomaly factor algorithm is used to calculate the anomaly score at each time point. Time points with anomaly scores greater than 1.5 are considered as anomaly segments and are removed. The length of the continuous valid interval formed after removal must not be less than 80% of the original length; otherwise, the entire batch of data is marked as invalid and excluded.

3. The method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology according to claim 2, characterized in that: The specific process for calculating the coefficient of variation of lactose content in each batch of raw milk, fitting the exponential decay curve of enzyme activity with storage time, and performing sliding window standard deviation analysis on the residence time series composed of residence time records is as follows: Calculate the coefficient of variation of raw milk lactose content in each batch. The coefficient of variation is equal to the weighted standard deviation divided by the weighted mean. The weighted standard deviation is the square root of the weighted average of the squared deviations of raw milk lactose content from the weighted mean at each sampling point. The weighted mean is the weighted average of raw milk lactose content at each sampling point. The Levenberg-Marquardt algorithm was used to fit the exponential decay curve of the batch activity calibration value of commercial lactase as the storage days. After fitting, the initial enzyme activity and decay constant were output as two independent features. For the residence time series consisting of records of the residence time of materials in the tank, a sliding window is used for standard deviation analysis. Specifically, the scaling value of the median absolute deviation within the window is calculated. The window length is fixed at five consecutive sampling points, and the window sliding step size is one sampling point. The maximum value of the scaling value of the median absolute deviation of all windows is taken as the result of the sliding window standard deviation analysis. The scaling value of the median absolute deviation is equal to 1.4826 multiplied by the median of the absolute deviation of the residence time within the window from the median of the median of the window. The coefficient 1.4826 makes the scaling value of the median absolute deviation a consistent estimator of the standard deviation. When calculating the median of the absolute deviation of the residence time within the window from the median of the window, the median of the residence time within the window is first calculated, then the absolute value of the difference between each residence time within the window and the median is calculated, and finally these absolute differences are sorted and the median is taken.

4. The method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology according to claim 3, characterized in that: The process of obtaining the structured temporal feature vector set is as follows: After feature extraction, the coefficient of variation, initial enzyme activity, decay constant, sliding window standard deviation analysis results, arithmetic mean of protein content in the batch rapid test values ​​of raw milk lactose and protein content, and arithmetic mean of workshop ambient temperature records for each batch are combined into a structured feature vector. The structured feature vectors of all batches are arranged in chronological order, and each feature dimension is standardized by subtracting the global mean and dividing by the global standard deviation to obtain the final time-series feature vector set.

5. The method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology according to claim 4, characterized in that: In step S2, the temporal feature vector set obtained in step S1 is first organized into a three-dimensional tensor in chronological order. The first dimension is the batch number, the second dimension is the time step number in each batch, and the third dimension is the feature dimension with a size of six. The temporal feature vector of each time step is the feature vector in the temporal feature vector set. The three-dimensional tensor is subjected to layer normalization processing, and the processed three-dimensional tensor is input into a deep long short-term memory network. The deep long short-term memory network is constructed as a two-layer stacked structure, with each layer containing 128 hidden units. At each time step, the deep long short-term memory network outputs hidden states and cell states, where the hidden states are 128-dimensional vectors. The hidden states output by the deep long short-term memory network at all time steps constitute a hidden state sequence. The network output layer contains two independent branch structures. One branch is a fully connected layer used to output the predicted lactose residue value at the target time step. The other branch connects the fully connected layer and the Softmax activation function in sequence to output the predicted flavor score value. During network training, a joint loss function is used, which consists of the sum of the mean squared error term, the cross-entropy loss term, and the L2 regularization term. In deep long short-term memory networks, the calculation process of the forget gate activation value introduces a modulation factor of the ambient temperature deviation. The calculation process is as follows: the difference between the workshop ambient temperature record at the current time step and the standard temperature of 42 degrees Celsius is obtained as the ambient temperature deviation. The ambient temperature deviation is compressed by the hyperbolic tangent function, multiplied by the learnable modulation factor, and then multiplied element-wise with the learnable projection vector and superimposed on the linear transformation result of the forget gate.

6. The method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology according to claim 5, characterized in that: The specific operation of introducing a multi-head attention mechanism after the hidden layers of a deep long short-term memory network is as follows: First, the hidden state sequence output by the last deep long short-term memory network is used as the input of the multi-head attention mechanism, and the number of attention heads is set to eight. Each attention head independently performs a linear transformation on the hidden state sequence to generate a query matrix, a key matrix, and a value matrix, respectively. The feature dimensions of the query matrix and the key matrix are set to sixteen, and the feature dimension of the value matrix is ​​set to sixteen. Subsequently, each attention head calculates the product of the query matrix and the transpose of the key matrix, divides the product by the square root of the feature dimension of the key matrix for scaling, and then performs Softmax normalization on the scaled result to obtain the attention weight matrix. The attention weight matrix is ​​then multiplied by the value matrix to calculate the weighted output representation of each attention head. Finally, the outputs of all eight attention heads are concatenated along the feature dimension, and the total dimension of the concatenation is mapped back to the hidden layer dimension through a linear projection layer, where the hidden layer dimension is 128. The mapping result is further processed by a non-linear activation function to obtain the final attention representation. After the attention representation is calculated, it is input into the network output layer to obtain the lactose residue prediction value and flavor score prediction value.

7. The method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology according to claim 6, characterized in that: In step S3, the specific operation of performing discrete wavelet decomposition on the original spectral matrix acquired online is as follows: Discrete wavelet decomposition uses the Dobessi type 5 mother wavelet with four decomposition levels. Each level generates a set of approximation coefficients and a set of detail coefficients. Adaptive soft thresholding is applied to the detail coefficients to filter out high-frequency noise and scattering components. The calculation process for the adaptive soft thresholding is as follows: First, obtain the detail coefficients of the current layer. Then, take the absolute value of the detail coefficients and calculate the median. Divide the median by 0.6745 to obtain the estimated noise standard deviation of the current layer. Then, multiply the estimated noise standard deviation by the square root of the number of wavelengths under the natural logarithm to calculate the threshold of the current layer. The specific operation of soft thresholding is as follows: iterate through each detail coefficient of the current layer. If the absolute value of the detail coefficient is less than or equal to the calculated threshold of the current layer, then set the value of the detail coefficient to zero. If the absolute value of the detail coefficient is greater than the threshold of the current layer, then subtract the threshold of the current layer from the detail coefficient and keep its original positive or negative sign unchanged. The inverse discrete wavelet transform is performed on all the layer detail coefficients after thresholding and the approximation coefficients of the fourth layer to reconstruct a denoised spectrum. The same discrete wavelet decomposition and adaptive soft thresholding are then performed on each column of the original spectral matrix at all time points. The denoised spectra at all time points together form the denoised spectral matrix.

8. The method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology according to claim 7, characterized in that: The specific process of establishing a partial least squares regression model with the absorbance value of the denoised spectrum as the independent variable and the offline test value of the lactose residue in the finished product as the dependent variable is as follows: Using the absorbance value of the denoised spectral matrix as the independent variable and the offline test value of the lactose residue in the finished product accumulated in step S1 as the dependent variable, the potential components of the independent and dependent variables are extracted by nonlinear iterative partial least squares algorithm, and the number of potential components is determined by cross-validation. After extracting the latent components, the weight matrix, loading matrix, and dependent variable loading vector are calculated sequentially. Then, the regression coefficient vector is solved through matrix operations. The solution process is as follows: First, calculate the inverse matrix of the product of the weight matrix and the loading matrix. Then, perform matrix multiplication of the inverse matrix with the dependent variable loading vector. The result is the regression coefficient vector, whose dimension is consistent with the number of spectral wavelength points. For each new original spectral vector acquired online, first perform the same discrete wavelet decomposition and adaptive soft thresholding as in step S3 to obtain the corresponding denoised spectrum. Then, perform an inner product operation on the denoised spectrum and the regression coefficient vector. The inner product result is the real-time lactose detection value.

9. The method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology according to claim 8, characterized in that: In step S4, the composite control law is composed of a superposition of a feedback term, a feedforward compensation term, a feedforward disturbance term, and a flavor constraint term. The calculation process of the feedback term is as follows: First, calculate the difference between the real-time lactose detection value and the preset lactose residue target value as the feedback error. Then, take the absolute value of the feedback error and perform a power nonlinear operation. Multiply the result of the power nonlinear operation by the positive or negative sign of the feedback error, and then multiply it by the feedback gain coefficient to obtain the value of the feedback term. The calculation process of the feedforward compensation term is as follows: First, calculate the difference between the predicted lactose residue value and the preset lactose residue target value as the feedforward error. Then, perform an integral operation on the feedforward error from time zero to the current time. Multiply the integral result by the feedforward integral gain coefficient to obtain the value of the feedforward compensation term. The calculation process for the feedforward disturbance term is as follows: multiply the start-stop cycle signal of the downstream filling unit by the disturbance feedforward gain coefficient to obtain the value of the feedforward disturbance term; The calculation process of the flavor constraint term is as follows: First, calculate the difference between the preset flavor threshold and the predicted flavor score. If the difference is positive, retain the difference; if the difference is negative, take zero. Multiply the retained difference by the flavor constraint weight coefficient to obtain the value of the flavor constraint term. The real-time adjustment increment of enzyme dosage is obtained by adding the feedback term, feedforward compensation term, feedforward perturbation term, and flavor constraint term.

10. The method for controlling the production of dairy products for lactose-intolerant individuals based on enzymatic hydrolysis technology according to claim 9, characterized in that: The process of automatically triggering the cooling standby mode and heating recovery mode of the enzymatic hydrolysis tank based on the cumulative downtime of filling is as follows: The program controller monitors the cumulative downtime of filling in real time. When the cumulative downtime exceeds the preset downtime threshold, the program controller automatically triggers the cooling standby mode and starts the cooling system to reduce the temperature inside the enzymatic hydrolysis tank to ten degrees Celsius to inhibit enzyme activity. When filling resumes and the cumulative downtime of filling is less than the preset restart temperature threshold, the program controller automatically triggers the temperature recovery mode, starts the heating system to raise the temperature inside the enzymatic hydrolysis tank to 42 degrees Celsius, and maintains it until the enzymatic hydrolysis is completed. A hysteresis interval is set between the cooling trigger point and the heating trigger point. When the cumulative downtime of filling is between the preset downtime threshold and the preset restart temperature threshold, the program controller maintains the current temperature mode unchanged.

Images