Coffee roasting process precision control method and system based on digital twinning

Abstract

This invention discloses a method and system for precise control of the coffee roasting process based on digital twins, specifically relating to the fields of industrial process control and digital twin technology. It addresses the problem of inaccurate control decisions caused by model prediction errors during nonlinear dynamic abrupt changes in the process, as seen in existing digital twin-based control methods. This invention identifies in real-time whether the coffee roasting process has entered a nonlinear dynamic abrupt change stage. Within this stage, it simultaneously quantifies the prediction reliability of the digital twin model and analyzes the strength of path dependence effects generated by its historical predictions. Based on the quantification results, it dynamically adjusts the model prediction weights. Finally, based on the adjusted weights, it fuses the prediction results of the digital twin model to generate precise control commands, thereby achieving stable and precise control of the coffee roasting process.

Description

Technical Field

[0001] This invention relates to the fields of industrial process control and digital twin technology, and more specifically, to a method and system for precise control of a coffee roasting process based on digital twins. Background Technology

[0002] In the field of industrial process control, to improve the control quality and intelligence level of complex physical processes, digital twin technology can be used to construct virtual models that operate synchronously with the physical entities. Especially in processes involving complex thermodynamic changes and chemical reactions, such as coffee roasting, existing methods use integrated real-time sensor data to drive digital twin models for state simulation and trend prediction, and generate control commands based on this, aiming to achieve precise tracking and control of the roasting curve and the quality of the final product.

[0003] However, such physical processes often include stages with highly significant nonlinear dynamic characteristics. During these stages, intense thermochemical reactions and physical structural transformations occur within the material, causing rapid drift in the parameters of the dynamic model of its state evolution. This leads to significant fluctuations and uncertainties in the prediction accuracy of digital twin models trained on steady-state or normal-state data. Existing control schemes typically treat the model's prediction outputs at different stages indiscriminately as the direct and sole basis for decision-making. Consequently, when the process reaches the aforementioned nonlinear abrupt change stage, the control system still generates control commands based on a model prediction whose reliability has been substantially reduced. This results in inaccurate control decisions and makes it difficult to reliably guarantee the final quality indicators of the process, such as the consistency of product flavor and texture. Summary of the Invention

[0004] In order to overcome the above-mentioned defects of the prior art, the present invention provides a method and system for precise control of coffee roasting process based on digital twin to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, the present invention provides the following technical solution: Precise control methods for the coffee roasting process based on digital twins include: S1: Real-time acquisition of multi-dimensional sensor data of the physical entity during coffee roasting, and simultaneous acquisition of prediction data generated by the digital twin model; S2: Map multi-dimensional sensor data to a low-dimensional feature space, and identify whether the coffee roasting process has entered a nonlinear dynamic mutation stage based on the geometric characteristics of the historical motion trajectory of the process state points in the low-dimensional feature space. S3: If the nonlinear dynamic mutation stage is entered, the deviation sequence between the predicted data and the multi-dimensional sensor data is calculated for multiple key process variables. The unexpected coupling relationship characteristics that emerge between the deviation sequences in the nonlinear dynamic mutation stage are analyzed to quantify the prediction credibility of the digital twin model. S4: Based on the identification results of the nonlinear dynamic mutation stage, combined with the historical control command sequence and the current process state, the strength of the path dependence effect of the digital twin model due to the influence of historical prediction is analyzed. S5: Dynamically adjust the model prediction weights based on prediction confidence and the strength of path dependence effect; S6: Based on the model prediction weights and the prediction results of the digital twin model, control instructions are generated to control the coffee roasting process.

[0006] Furthermore, S1 includes: Real-time acquisition of multi-dimensional sensor data reflecting the thermodynamic and chemical states of coffee roasting process; Input multi-dimensional sensor data into the digital twin model; Obtain predictive data that is temporally aligned with the multidimensional sensor data, generated by forward simulation calculations based on multidimensional sensor data and the internal state of the model, from the digital twin model.

[0007] Furthermore, S2 includes: Dimensionality reduction processing is performed on multi-dimensional sensor data to obtain process state points represented in a low-dimensional feature space; By connecting continuous process state points in time sequence, the historical motion trajectory of process state points in a low-dimensional feature space is formed. Calculate the changes in geometric characteristics of historical motion trajectories within the current and adjacent time windows; Determine whether the change in geometric characteristics continuously exceeds the preset abrupt change threshold and reaches the preset duration; If so, the coffee roasting process is determined to have entered a non-linear dynamic mutation stage.

[0008] Furthermore, S3 includes: For each key process variable, within the time window of the nonlinear dynamic mutation stage, calculate the time-ordered sequence of deviations between the predicted data and the multi-dimensional sensor data; Calculate the coupling strength measure between each pair of deviation sequences to obtain the coupling relationship matrix in the nonlinear dynamic mutation stage; Obtain the baseline coupling relationship matrix established in advance during the steady-state phase of the process; Calculate the difference measure between the coupling relationship matrix and the baseline coupling relationship matrix in the nonlinear dynamic mutation stage, and use the difference measure to characterize the features of unexpected coupling relationships; Based on the saliency of the unexpected coupling relationship features, a predictive reliability quantification of the digital twin model is generated by mapping.

[0009] Furthermore, the coupling strength measure between each pair of deviation sequences is calculated, including: for each pair of deviation sequences corresponding to key process variables, within the time window of the nonlinear dynamic mutation stage, the maximum information coefficient between the two deviation sequences is calculated, and the value of the maximum information coefficient is used as the coupling strength measure between the deviation sequences corresponding to this pair of key process variables.

[0010] Furthermore, S4 includes: Obtain the historical control command sequence and current process status within the time period preceding the entry into the nonlinear dynamic mutation stage; By combining the historical control command sequence with the reference model of the process's steady phase, a reference process state sequence that is expected to be reached if there is no model prediction bias is generated. Calculate the deviation between the current process state and the expected state of the reference process state sequence at the same time. The path dependence effect strength is obtained by weighting and correcting the state deviation based on the radicalness of the command changes in the historical control command sequence.

[0011] Furthermore, S5 includes: Obtain the quantified prediction credibility and the analyzed path dependence effect strength; Based on the preset weight adjustment mapping relationship, the combined input of prediction confidence and path dependence effect strength is mapped to the corresponding model prediction weight adjustment amount; The current model prediction weights are updated based on the model prediction weight adjustment amount to generate dynamically adjusted model prediction weights.

[0012] Furthermore, based on the preset weight adjustment mapping relationship, the combined input of prediction confidence and path dependence effect strength is mapped to the corresponding model prediction weight adjustment amount, including: The prediction confidence and the strength of the path dependence effect are input into a pre-defined fuzzy inference system; Through the inference calculation of the fuzzy inference system, the model prediction weight adjustment amount corresponding to the combination of prediction credibility and path dependence effect strength is output.

[0013] Furthermore, S6 includes: Based on the dynamically adjusted model prediction weights, the prediction results of the digital twin model are weighted and fused to generate a weighted process state prediction. Based on the weighted process state prediction and the preset baking process target state, the target control command is calculated through the control algorithm; The target control command is output to the actuator of the coffee roasting process to drive the actuator to control the coffee roasting process.

[0014] On the other hand, the present invention provides a precise control system for a coffee roasting process based on digital twins, comprising: The data acquisition module is used to collect multi-dimensional sensor data of the physical entity during the coffee roasting process in real time, and simultaneously acquire the prediction data generated by the digital twin model; The mutation identification module is used to map multi-dimensional sensor data to a low-dimensional feature space and identify whether the coffee roasting process has entered a nonlinear dynamic mutation stage based on the geometric characteristics of the historical motion trajectory of the process state points in the low-dimensional feature space. The credibility assessment module is used to calculate the deviation sequence between the predicted data and the multi-dimensional sensor data for multiple key process variables if the nonlinear dynamic mutation stage is entered. It analyzes the unexpected coupling relationship characteristics that emerge between the deviation sequences in the nonlinear dynamic mutation stage to quantify the prediction credibility of the digital twin model. The effect analysis module is used to analyze the strength of the path dependence effect of the digital twin model caused by the influence of historical predictions, based on the identification results of the nonlinear dynamic mutation stage combined with the historical control command sequence and the current process state. The weight adjustment module is used to dynamically adjust the model prediction weights based on prediction confidence and the strength of path dependence effect; The instruction generation module is used to generate control instructions based on the model prediction weights and the prediction results of the digital twin model to control the coffee roasting process.

[0015] Compared with the prior art, the present invention has the following beneficial effects: The beneficial effects provided by this invention are mainly reflected in the following two aspects.

[0016] 1. By introducing a refined real-time evaluation mechanism for the reliability of digital twin model outputs and a dynamic adjustment mechanism for decision weights, the traditional control scheme's reliance on model predictions is changed. By identifying nonlinear dynamic mutation stages and deeply analyzing the unexpected coupling relationships emerging between model prediction errors and the path dependence effect caused by historical predictions during these stages, the true reliability of the model under the current complex dynamics can be quantitatively and multi-dimensionally evaluated. Based on this evaluation result, the model prediction weights are dynamically adjusted, so that the decision basis of the control system is no longer a single and fixed model output, but an adaptive fusion information based on the real-time changes in model performance. This effectively overcomes the problem of inaccurate control decisions caused by blindly relying on inaccurate models during drastic process changes, thereby significantly improving the control accuracy and stability of process states in key process stages.

[0017] 2. By generating the final control command through a weighted fusion of model predictions and predictions based on classical feedback principles, the system can fully utilize its predictive advantage when the model's reliability is high, while automatically enhancing the robust control component based on real-time feedback when the model's reliability is low. This smooth decision transition ensures the consistency and stability of control behavior, avoiding drastic changes in control commands caused by model performance fluctuations. By ensuring that control commands are always generated based on the most reliable information source throughout the baking process, especially in the nonlinear dynamic mutation stage, the system guarantees high consistency of the quality indicators of the end product. Attached Figure Description

[0018] Figure 1 This is a flowchart of the precise control method for the coffee roasting process based on digital twins according to the present invention; Figure 2 This is a schematic diagram of the structure of the precision control system for the coffee roasting process based on digital twins according to the present invention. Detailed Implementation

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

[0020] Example 1: Figure 1 The present invention provides a method for precise control of the coffee roasting process based on digital twins, comprising: S1: Real-time acquisition of multi-dimensional sensor data of the physical entity during coffee roasting, and simultaneous acquisition of prediction data generated by the digital twin model; S2: Map multi-dimensional sensor data to a low-dimensional feature space, and identify whether the coffee roasting process has entered a nonlinear dynamic mutation stage based on the geometric characteristics of the historical motion trajectory of the process state points in the low-dimensional feature space. S3: If the nonlinear dynamic mutation stage is entered, the deviation sequence between the predicted data and the multi-dimensional sensor data is calculated for multiple key process variables. The unexpected coupling relationship characteristics that emerge between the deviation sequences in the nonlinear dynamic mutation stage are analyzed to quantify the prediction credibility of the digital twin model. S4: Based on the identification results of the nonlinear dynamic mutation stage, combined with the historical control command sequence and the current process state, the strength of the path dependence effect of the digital twin model due to the influence of historical prediction is analyzed. S5: Dynamically adjust the model prediction weights based on prediction confidence and the strength of path dependence effect; S6: Based on the model prediction weights and the prediction results of the digital twin model, control instructions are generated to control the coffee roasting process.

[0021] Step S1 is implemented in the following manner: Multiple sensors are deployed at key locations within the physical structure of the coffee roaster to collect multi-dimensional sensor data reflecting the thermodynamic and chemical states of the process in real time. Thermodynamic state sensors specifically include K-type thermocouple temperature sensors installed in the coffee bean bed area, measuring the coffee bean bed temperature over a range of 50 to 250 degrees Celsius; platinum resistance temperature sensors installed in the hot air duct, measuring the hot air temperature; and surface-mount thermocouples installed on the drum wall, measuring the drum wall temperature. Chemical state sensors include a metal oxide semiconductor gas sensor array installed in the exhaust duct, containing sensor elements sensitive to carbon dioxide and volatile organic compounds; and a near-infrared spectral probe installed at the viewing window, illuminating the coffee bean surface with a specific wavelength range and receiving the reflected spectrum. All sensors operate at a fixed sampling frequency, determined based on the fastest timescale of key physical changes during coffee roasting, for example, setting the sampling frequency to collect data 10 times per second. Analog measurement signals from all sensors are synchronously sampled and converted from analog to digital via a multi-channel data acquisition card, forming a time-strictly aligned multi-dimensional sensor data stream. Each data point in the multi-dimensional sensor data stream contains the measurement values ​​of all sensors at the same moment, and each data point is marked with a timestamp generated by a unified clock source. The multi-dimensional sensor data is transmitted in real time to a computing unit deployed with a digital twin model via an industrial Ethernet network based on the TCP / IP protocol.

[0022] The computing unit runs a pre-established digital twin model of the coffee roasting process. The digital twin model of the coffee roasting process is a computational model that integrates the physical mechanisms of heat and mass transfer with the kinetics of key chemical reactions. The internal state variables of the digital twin model of the coffee roasting process include coffee bean core temperature, coffee bean surface temperature, overall coffee bean moisture content, Maillard reaction progress, and caramelization reaction progress.

[0023] The method for establishing a digital twin model of the coffee roasting process is as follows: Based on the physical process of heat exchange between coffee beans and hot air and the drum wall inside the roasting drum, a heat transfer equation describing the temperature changes of the coffee bean core and surface is established; based on the diffusion and evaporation of moisture inside the coffee beans, a mass transfer equation describing the changes in the moisture content of the coffee beans is established; for the Maillard reaction and caramelization reaction, a reaction kinetic equation describing the changes in the reaction progress is established; the heat transfer equation, mass transfer equation, and reaction kinetic equation are combined to form a set of differential equations describing the complete physicochemical behavior of the coffee roasting process.

[0024] The rate of change of coffee bean core temperature is calculated as: dTcore / dt=(hcore×Acore×(Tsurface−Tcore)) / (m×Cpcore); where dTcore / dt represents the rate of change of coffee bean core temperature, hcore represents the equivalent heat transfer coefficient between the bean core and the bean surface, Acore represents the equivalent heat transfer area of ​​the bean core, Tsurface represents the surface temperature of the coffee bean, Tcore represents the temperature of the coffee bean core, m represents the mass of the coffee bean, and Cpcore represents the specific heat capacity of the coffee bean core.

[0025] The rate of change of coffee bean surface temperature is calculated as: dTsurface / dt=(hair×Abean×(Tair−Tsurface)+hwall×Awall×(Twall−Tsurface)−hcore×Acore×(Tsurface−Tcore)) / (msurface×Cpsurface); where dTsurface / dt represents the rate of change of coffee bean surface temperature, hair represents the convective heat transfer coefficient between hot air and coffee bean surface, Abean represents the coffee bean surface area, Tair represents the hot air temperature, hwall represents the equivalent heat transfer coefficient between the drum wall and coffee bean, Awall represents the contact area between coffee bean and drum wall, Twall represents the drum wall temperature, msurface represents the coffee bean surface mass, and Cpsurface represents the coffee bean surface specific heat capacity.

[0026] The rate of change of coffee bean moisture content is calculated as: dM / dt = −kevap × (M − Meq(Tair)); where dM / dt represents the rate of change of coffee bean moisture content, kevap represents the moisture evaporation rate constant, M represents the current moisture content of coffee beans, and Meq(Tair) represents the equilibrium moisture content at the hot air temperature Tair, which is a function of the hot air temperature.

[0027] The Maillard reaction progress rate is calculated as: dRmaillard / dt = kmaillard × exp(−Eamaillard / (Rgas × Tcore)) × fmaillard(M); where dRmaillard / dt represents the Maillard reaction progress rate, kmaillard represents the Maillard reaction pre-exponential factor, Eamaillard represents the Maillard reaction activation energy, Rgas represents the universal gas constant, Tcore represents the coffee bean core temperature, and fmaillard(M) represents the correction function of moisture content on the Maillard reaction rate. The specific form of fmaillard(M) is determined by fitting experimental data.

[0028] The rate of change of the caramelization reaction progress is calculated as: dRcaramel / dt=kcaramel×exp(−Eacaramel / (Rgas×Tcore))×fcaramel(M); where dRcaramel / dt represents the rate of change of the caramelization reaction progress, kcaramel represents the pre-exponential factor of the caramelization reaction, Eacaramel represents the activation energy of the caramelization reaction, Rgas represents the universal gas constant, Tcore represents the coffee bean core temperature, and fcaramel(M) represents the correction function of moisture content on the caramelization reaction rate. The specific form of fcaramel(M) is determined by fitting experimental data.

[0029] The model parameters in the aforementioned differential equation system include equivalent heat transfer coefficients, convective heat transfer coefficients, moisture evaporation rate constants, pre-exponential factors, activation energies, and empirical coefficients in the correction function. These model parameters are determined as follows: For a specific coffee bean variety, roasting calibration experiments are designed under different temperature conditions. During the experiments, the changes in hot air temperature, drum wall temperature, coffee bean bed temperature, and coffee bean moisture content are recorded simultaneously. At the same time, the content of Maillard reaction products and caramelization reaction products in the coffee beans at different times is measured offline using chemical analysis. The experimental data is used as a training set, and the least squares method is employed to optimize and fit the model parameters in the differential equation system, minimizing the error between the model simulation output and the experimental measurement data. The calibrated model parameters are then stored and used for online prediction of that specific coffee bean variety. For different varieties of green coffee beans, the above parameter calibration process must be performed separately to obtain the corresponding model parameter sets.

[0030] Real-time multi-dimensional sensor data is input into the digital twin model of the coffee roasting process. Specifically, the hot air temperature measurement and drum wall temperature measurement from the multi-dimensional sensor data are used as the thermal boundary conditions input to the digital twin model. The coffee bean bed temperature measurement, exhaust gas sensor response value, and near-infrared spectral characteristic value from the multi-dimensional sensor data are used as the state observation inputs to the digital twin model. Based on the input multi-dimensional sensor data and the internal state variable values ​​updated and saved in the previous calculation cycle, the digital twin model performs forward simulation calculations.

[0031] The forward simulation calculation process is achieved by numerically solving a set of differential equations, which describe the conservation of energy, mass, and chemical reaction kinetics. The digital twin model of the coffee roasting process uses the input hot air temperature and drum wall temperature, combined with convection and conduction heat transfer coefficients, to calculate the heat transferred to the coffee beans per unit time, and updates the values ​​of the coffee bean core temperature variable and the coffee bean surface temperature variable accordingly. Using the updated temperature and coffee bean moisture content variables, the digital twin model calculates the instantaneous reaction rates of the Maillard reaction and caramelization reaction using reaction kinetic formulas. The reaction kinetic formulas adopt the form of the Arrhenius equation, and the activation energy parameter and pre-exponential factor parameter in the Arrhenius equation are obtained through calibration using roasting experimental data for specific coffee bean varieties. The digital twin model updates the reaction progress variable based on the instantaneous reaction rate and simultaneously updates the moisture evaporation rate. The change in reaction progress variable drives the color evolution sub-model and volatile matter generation sub-model within the digital twin model of the coffee roasting process. The color evolution sub-model outputs a predicted change in the surface color of the coffee beans based on the reaction progress, while the volatile matter generation sub-model outputs a predicted change in the rate of volatile gas generation based on the reaction progress and temperature. The coupled differential equation system is solved numerically, using methods such as the Euler method or the Runge-Kutta method. The digital twin model of the coffee roasting process extrapolates forward by a predetermined time step. The time step must satisfy a numerical stability condition, requiring that the time step be less than a certain proportion of the system's fastest dynamic time constant; for example, a time step of 0.1 seconds can be used. The predicted values ​​of all internal state variables of the digital twin model of the coffee roasting process at this future moment are obtained through this extrapolation. These predicted internal state variables are transformed into physical quantities that are exactly the same type of measurements as those taken by the physical sensors through output mapping functions. These output mapping functions include functions that map the predicted coffee bean surface temperature to the predicted bean bed region temperature, functions that map the predicted volatile gas concentration to the predicted metal-oxide-semiconductor sensor response, and functions that map the predicted color state to the predicted near-infrared spectral characteristic values. The data generated by these output mapping functions is the prediction data generated by the digital twin model.

[0032] To ensure temporal comparability between predicted data and real-time acquired multi-dimensional sensor data, the computation cycle and data sampling cycle of the digital twin model for the coffee roasting process are coordinated. Each forward simulation calculation begins at the time corresponding to the latest received multi-dimensional sensor data, and timestamp management ensures that the future time point corresponding to the output predicted data corresponds to the time point of subsequently received multi-dimensional sensor data. When computation time results in a fixed delay in the output of predicted data, a caching and pairing mechanism is employed to pair the predicted data with subsequent multi-dimensional sensor data that matches the timestamp, ensuring that the data pairs used for subsequent comparative analysis point to the same physical time. The predicted data generated by the acquired digital twin model and the real-time acquired multi-dimensional sensor data are used together as input for subsequent steps.

[0033] Step S2 is implemented in the following manner: Dimensionality reduction is performed on the multi-dimensional sensor data from step S1. Multi-dimensional sensor data is a high-dimensional dataset containing time series of multiple different physical quantities. Examples of multi-dimensional sensor data include coffee bean bed temperature series, hot air temperature series, drum wall temperature series, response value series of multiple gas sensors, and absorbance series of multiple wavelengths in the near-infrared spectrum. Dimensionality reduction is achieved using principal component analysis. Specifically, a multi-dimensional sensor data sample within a historical time window is extracted from the real-time data stream. The length of the historical time window should cover a relatively stable stage in the coffee roasting process; for example, the length of the historical time window can be set to 120 seconds. These multi-dimensional sensor data samples undergo standardization preprocessing, which involves subtracting the mean of each sensor data sequence within the historical time window and dividing by its standard deviation within the historical time window. The covariance matrix of the standardized preprocessed data is calculated. Eigenvalue decomposition is performed on the covariance matrix to calculate the eigenvalues ​​and corresponding eigenvectors. The eigenvectors corresponding to the two largest eigenvalues ​​are selected as projection basis vectors. The two-dimensional plane spanned by the eigenvectors corresponding to the two largest eigenvalues ​​constitutes the required low-dimensional feature space. The two coordinate axes of the low-dimensional feature space have physical meaning: the first principal component direction typically represents the total energy input of the process, and the second principal component direction typically represents the intensity of the chemical reaction. The multi-dimensional sensor data vectors arriving in real time undergo the same standardized preprocessing, and then a dot product operation is performed with each of the two projected basis vectors. The result of the dot product operation is the coordinate of that time point in the low-dimensional feature space; this coordinate point is called the process state point. Through continuous calculation, a series of process state points arranged in chronological order are obtained.

[0034] By connecting consecutive process state points in time sequence, a historical trajectory of these process state points in a low-dimensional feature space is formed. The connection method involves sequentially connecting process state points with adjacent time stamps using straight line segments, forming a polygonal trajectory. This polygonal trajectory records the path of the process state's evolution over time in a two-dimensional feature plane defined by total energy input and the intensity of chemical reactions. The length of the historical trajectory is determined by the number of process state points within the considered time window; for example, continuously tracking all process state points generated within the last 60 seconds to form the current historical trajectory.

[0035] Calculate the changes in geometric properties of the historical trajectory within the current and adjacent time windows. The geometric property chosen is the curvature of the historical trajectory. Curvature is a geometric quantity that describes the degree of bending of a curve. For a historical trajectory composed of discrete process state points, the instantaneous curvature approximation of each intermediate point at the current moment is calculated using the adjacent three-point method. Specifically, for three consecutive process state points, denoted as points A, B, and C, with the time sequence corresponding to points A, B, and C increasing sequentially; calculate the vector AB pointing from point A to point B and the vector BC pointing from point B to point C; calculate the change in the angle between vectors AB and BC; divide the change in the angle between vectors AB and BC by the average time interval from point B to points A and C to obtain the rate of change of the angular velocity of the trajectory near point B; divide the rate of change of the angular velocity of the trajectory near point B by the approximate velocity modulus of the trajectory near point B, which is obtained by calculating the average of the modulus of vectors AB and BC and then dividing by the time interval; finally, a dimensionless quantity reflecting the instantaneous curvature of the trajectory at point B is obtained, namely the instantaneous curvature value. The change in geometric characteristics of a historical trajectory within the current and adjacent time windows is measured by calculating the absolute value of the difference between the instantaneous curvature value at the current moment and the instantaneous curvature value a short time ago. For example, the absolute value of the difference between the current instantaneous curvature and the instantaneous curvature 5 seconds ago is defined as the curvature change. The curvature change characterizes whether the drastic change in the direction of motion of the process state per unit time is intensifying.

[0036] The system determines whether the change in geometric characteristics continuously exceeds a preset mutation threshold and reaches a preset duration. The preset mutation threshold is a key parameter, determined through analysis of extensive historical coffee roasting data. Specifically, it involves collecting data from multiple normal roasting batches and calculating the statistical distribution of curvature changes throughout the roasting process, such as calculating the mean and standard deviation. Simultaneously, it collects data from multiple roasting batches empirically considered to have experienced significant nonlinear dynamic mutations, identifies the mutation periods, and calculates typical values ​​of curvature changes within those periods. By comparing the normal fluctuation range with the typical mutation value, a value between the two is determined as the preset mutation threshold; for example, the preset threshold could be set as the normal fluctuation mean plus three times the standard deviation. The preset duration is another parameter, used to avoid misjudgments due to instantaneous spikes caused by noise. It is determined based on the typical shortest duration of the nonlinear dynamic mutation phase; for example, through analysis of historical mutation events, the preset duration could be set to 8 seconds. The judgment logic is as follows: starting from the current moment, backtracking forward, check whether the curvature change at each calculation moment is greater than the preset mutation judgment threshold within a continuous period of time, such as 8 consecutive seconds; if the curvature change at all moments exceeds the preset mutation judgment threshold, and the duration of this continuous state exceeding the preset mutation judgment threshold reaches the preset duration, then the judgment condition is met.

[0037] If the judgment result is yes, meaning the change in geometric characteristics continuously exceeds the preset mutation judgment threshold and reaches the preset duration, then the coffee roasting process is determined to have entered the nonlinear dynamic mutation stage. This judgment signifies a fundamental change in the dynamic characteristics of the coffee roasting process. This judgment result is output as a Boolean type flag signal and used to trigger the specific analysis processes in subsequent steps S3 and S4. If the judgment result is no, the judgment that the coffee roasting process has not entered the nonlinear dynamic mutation stage is maintained, and the specific analysis processes for the nonlinear dynamic mutation stage in subsequent steps S3 and S4 will not be activated or will run using the default smooth mode parameters.

[0038] Step S3 is implemented in the following manner: When step S2 determines that the coffee roasting process has entered a nonlinear dynamic mutation phase, the specific analysis process for this step is initiated. Several key process variables are selected from the real-time multi-dimensional sensor data obtained in step S1 and the predicted data generated by the digital twin model. Key process variables refer to physicochemical quantities that have a decisive impact on coffee roasting quality and can be directly or indirectly measured by sensors. Examples of key process variables include coffee bean bed temperature, hot air temperature, and characteristic values ​​of carbon monoxide concentration in the exhaust gas. For each selected key process variable, specific processing is performed within the duration of the nonlinear dynamic mutation phase. The starting point of the nonlinear dynamic mutation phase time window is determined by the determination time in step S2, and the ending point is determined by the continuous real-time determination results of the nonlinear dynamic mutation phase. For example, the time window of the nonlinear dynamic mutation phase may cover a period of 20 seconds from the determination time. Within this time window, for each sampling time, the difference between the predicted data value and the multi-dimensional sensor data value corresponding to the key process variable is calculated. The difference is calculated by subtracting the multi-dimensional sensor data value from the predicted data value. The differences of all sampling times are arranged in chronological order to form a deviation sequence for the key process variable. A bias sequence is a time series in which each data point represents the error in the digital twin model's prediction of that variable at a specific time.

[0039] The coupling strength metric between each pair of deviation sequences is calculated to obtain the coupling relationship matrix during the nonlinear dynamic mutation stage. The maximum information coefficient method is used to calculate the coupling strength metric. For each pair of deviation sequences corresponding to key process variables, such as the coffee bean bed temperature deviation sequence and the hot air temperature deviation sequence forming a pair, the following calculations are performed. Data points at the same time stamp of the two deviation sequences are extracted to form a set of two-dimensional data points. The numerical range of each deviation sequence is discretized and binned. Binning is the process of dividing a continuous numerical interval into a finite number of discrete intervals. The number of bins needs to be determined experimentally to ensure statistical significance; for example, the number of bins can be set to an integer between 10 and 30. Based on the binning results, the joint probability distribution of the two deviation sequences is calculated. The joint probability distribution is obtained by counting the frequency of each two-dimensional data point falling into different bin combinations and dividing by the total number of data points. At the same time, the marginal probability distribution of each deviation sequence is calculated. The mutual information between two biased sequences is calculated based on their joint and marginal probability distributions. The method involves multiplying the joint probability by the base-2 logarithm of the ratio of the joint probability to the product of the two marginal probabilities for all possible binning combinations, and then summing the results for all combinations. Calculating the maximum information coefficient requires searching for different binning schemes to maximize mutual information. These schemes include changing the number of bins or adjusting the binning boundaries. The final maximum information coefficient is the standardized mutual information value, which is obtained by dividing the found maximum mutual information value by the minimum entropy of the two sequences under a given binning. The calculated maximum information coefficient is used as a measure of the coupling strength between the pair of biased sequences. Repeat the above calculation process for all key process variables pairwise, and fill the calculated coupling strength measure into a symmetric matrix. The rows and columns of the symmetric matrix correspond to each key process variable. The element in the i-th row and j-th column of the symmetric matrix is ​​the coupling strength measure between the deviation sequences of the i-th variable and the j-th variable. This symmetric matrix is ​​the coupling relationship matrix in the nonlinear dynamic mutation stage, where i is the row number in the symmetric matrix and j is the column number in the symmetric matrix.

[0040] Obtain a pre-established baseline coupling matrix during the stable phase of the process. A stable phase refers to a stage in the coffee roasting process where nonlinear dynamic characteristics are not significant. Such a phase may include the bean preheating stage in the early stages of roasting or the stable temperature rise stage before the abrupt change. The baseline coupling matrix is ​​established by identifying multiple time periods of stable phases in historical roasting data. Within these time periods, the deviation sequences of each key process variable are calculated using the same method as for calculating the coupling matrix, and the maximum information coefficients between each pair of variables are calculated. The calculation results from multiple stable phase time periods are averaged using an arithmetic mean to obtain a set of baseline coupling strength measures representing typical coupling relationships during the stable phase. These baseline coupling strength measures are arranged into a symmetric matrix according to the same variable order; this symmetric matrix is ​​the baseline coupling matrix. The baseline coupling matrix is ​​pre-calculated and stored during system initialization. The baseline coupling matrix characterizes the normal statistical correlation strength between the prediction errors of each variable when the digital twin model prediction is relatively accurate and the process dynamics are stable.

[0041] A difference measure is calculated between the coupling matrix in the nonlinear dynamic mutation phase and the baseline coupling matrix to characterize unexpected coupling features. The difference measure quantifies the overall deviation between the two matrices. One method for calculating the difference measure is to calculate the Frobenius norm difference between the two matrices. Specifically, the coupling matrix in the nonlinear dynamic mutation phase is subtracted element-wise from the baseline coupling matrix to obtain the difference matrix. Then, the Frobenius norm of the difference matrix is ​​calculated by squaring each element in the difference matrix, summing all the squares, and then taking the square root of the sum. The resulting norm value is the difference measure. The difference measure is a non-negative scalar value. A larger difference measure value indicates a more significant change in the coupling pattern between the prediction errors of various variables in the nonlinear dynamic mutation phase compared to the normal pattern in the stationary phase; this change is defined as an unexpected coupling feature. The significance of the unexpected coupling feature is directly reflected by the magnitude of the difference measure value.

[0042] Based on the saliency of the unexpected coupling relationship characteristics, a predictive credibility quantification value for the digital twin model is generated through mapping. This mapping relationship is implemented using a pre-defined function. This function takes the difference metric as input and outputs a predictive credibility quantification value between 0 and 1, where 1 represents complete confidence and 0 represents complete unconfidence. The function is established based on the analysis of historical data. Specifically, historical data containing different baking scenarios are collected, and the difference metric for each scenario's nonlinear dynamic mutation stage is calculated. Experts or the final product quality is used to inversely evaluate the predictive credibility of the digital twin model under that scenario. The correspondence between the difference metric value and the predictive credibility quantification value is determined through regression analysis or piecewise linear interpolation. For example, a piecewise linear mapping function can be established. When the difference metric is below a set low threshold, the predictive credibility quantification value remains at a high level, for example, above 0.9. When the difference metric is between a set low threshold and a set high threshold, the predictive credibility quantification value linearly decreases from 0.9 to 0.2 as the difference metric increases. When the difference metric exceeds a set high threshold, the predictive credibility quantification value remains at a low level below 0.2. The determination of the set low and high thresholds is based on statistical analysis of historical data. For example, the low threshold can be set as the upper limit of the fluctuation range of the baseline coupling relationship matrix itself, and the high threshold can be set as twice the low threshold. Through this mapping function, the calculated difference metric is converted into a specific, standardized predictive credibility metric of the digital twin model. This predictive credibility metric of the digital twin model is used for decision fusion in the subsequent step S5.

[0043] Step S4 is implemented in the following manner: When step S2 determines that the coffee roasting process has entered a nonlinear dynamic mutation phase, the specific analysis process for this step is initiated. The historical control command sequence for the period preceding the entry into the nonlinear dynamic mutation phase is obtained. The historical control command sequence refers to a series of control commands actually sent to actuators such as heaters or fans within a certain period prior to the moment the nonlinear dynamic mutation phase was identified in step S2, extracted from the control command execution log of the coffee roasting control system. The determination of the length of the period preceding the entry into the nonlinear dynamic mutation phase needs to consider the response delay of the control system and the inertia of the process dynamics. The length of the period preceding the entry into the nonlinear dynamic mutation phase must be set to ensure that it covers past control actions that significantly affect the current state. For example, the set value for the length of the period preceding the entry into the nonlinear dynamic mutation phase can be the most recent 60 seconds before the moment the nonlinear dynamic mutation phase begins. The historical control command sequence is an array arranged in chronological order. Each element in the historical control command sequence contains a timestamp and a command value. The command value is, for example, the power percentage setting value of the hot air heater, with a range of 0 to 100.

[0044] Simultaneously, the current process state is acquired. The current process state refers to the process state vector corresponding to the moment in step S2 when the nonlinear dynamic mutation stage is determined. The current process state is not directly derived from a single sensor reading, but rather is a comprehensive state representation obtained by feature extraction and fusion of the multi-dimensional sensor data acquired in step S1. One method for acquiring the current process state is to project the multi-dimensional sensor data vector at the moment of determining the entry into the nonlinear dynamic mutation stage onto the same low-dimensional feature space as in step S2, using the same projection basis vector as in the dimensionality reduction process. The resulting low-dimensional coordinate vector is then used as the current process state. The current process state is a vector containing multiple dimensions, with the dimension of the current process state being the same as the dimension of the low-dimensional feature space. For example, the current process state could be a two-dimensional vector containing coordinate values ​​representing the total energy input and coordinate values ​​representing the intensity of the chemical reaction.

[0045] By combining historical control command sequences with a reference model of the process's steady-state phase, a reference process state sequence expected to be achieved without model prediction bias is generated. The reference model for the steady-state phase is a simplified but reliable dynamic model used to describe the ideal situation where the process state should evolve in response to control commands during the coffee roasting process's steady-state phase, assuming the digital twin model's predictions are completely accurate. The method for establishing the reference model for the steady-state phase involves selecting a segment of historical data from the steady-state phase during system initialization or calibration. This historical data includes the control command sequence and the corresponding verified accurate process state sequence. A dynamic model is fitted using a system identification method, such as a least-squares parameter estimation algorithm. The fitted dynamic model can predict the evolution of the historical process state sequence with a small error based on the historical control command sequence; this fitted dynamic model is the reference model for the steady-state phase. Generating the reference process state sequence is an open-loop simulation process. The known process state corresponding to the start time of the time period before entering the nonlinear dynamic mutation phase is used as the initial simulation state. The acquired historical control command sequences are then progressively input into the reference model for the steady-state phase in chronological order. During the steady-state phase, the reference model calculates the predicted process state for the next step based on the input commands and the model's current internal state, with a fixed simulation step size. The simulation step size must meet numerical stability requirements; for example, it can be set to 1 second. This calculation is performed iteratively, ultimately yielding a predicted process state sequence from the start time to the time of determination of the nonlinear dynamic mutation phase. This sequence is the reference process state sequence. The reference process state sequence describes the expected trajectory that the coffee roasting process should have followed under the same historical control commands, provided it consistently adhered to the ideal dynamic model of the steady-state phase and was not affected by prediction errors from the digital twin model.

[0046] Calculate the state deviation between the current process state and the expected state of the reference process state sequence at the same time. The same time refers to the moment in step S2 when the process is determined to enter the nonlinear dynamic mutation stage. Find the expected state corresponding to the moment the nonlinear dynamic mutation stage is determined from the generated reference process state sequence; the expected state is the last state vector in the reference process state sequence. The state deviation is a scalar value quantifying the overall difference between the current process state vector and the expected state vector. The state deviation can be calculated using the Euclidean distance formula. Specifically, calculate the difference between the current process state vector and the expected state vector in each corresponding dimension. Square the difference in each dimension. Summate the squared values ​​in all dimensions. Take the square root of the summation result; the result is the state deviation. The state deviation is a non-negative value; the larger the state deviation, the greater the deviation between the actual process state and the ideal reference trajectory at the critical point of entering the nonlinear dynamic mutation stage. This deviation is initially attributed to the influence of inaccurate model predictions on historical control commands, leading the process to be guided into an unexpected state region.

[0047] The path dependence effect strength is obtained by weighting and correcting the state deviation based on the aggression of command changes in the historical control command sequence. The aggression of command changes measures whether historical control commands were frequently and significantly adjusted before entering the nonlinear dynamic mutation phase. The aggression of command changes is quantified by calculating the sum of the absolute values ​​of the differences between adjacent commands in the historical control command sequence. Specifically, for each pair of adjacent commands in the historical control command sequence, the absolute value of the difference between the subsequent command value and the preceding command value is calculated. The absolute values ​​of all adjacent command differences are summed to obtain the total absolute change. The total absolute change is divided by the total time length covered by the historical control command sequence, which is equal to the length of the time period before entering the nonlinear dynamic mutation phase, to obtain the average command change intensity per unit time. This average command change intensity is the quantified value of the aggression of command changes. The purpose of weighting correction is to correlate the state deviation with the aggression of command changes to reflect the contribution weight of inappropriate control commands to the state deviation. One method of weighting correction is multiplicative weighting. Multiplicative weighting multiplies the state deviation by a weighting factor determined by the aggressiveness of the command change. Determining the weighting factor requires a baseline value. This baseline value is obtained by analyzing the aggressiveness of command changes during normal control in historical stationary phases. Specifically, this is done by collecting control command sequences from multiple historical stationary phases, calculating the quantified aggressiveness of command changes for each sequence, and then taking the average of these quantified values ​​as the baseline value. The weighting factor is designed as the ratio of the quantified aggressiveness of command changes to the baseline value. When the quantified aggressiveness of command changes equals the baseline value, the weighting factor is 1, and the state deviation is neither amplified nor reduced. When the quantified aggressiveness of command changes is greater than the baseline value, the weighting factor is greater than 1, and the state deviation is amplified, indicating that the aggressive command amplifies the path dependence effect. When the quantified aggressiveness of command changes is less than the baseline value, the weighting factor is less than 1, and the state deviation is reduced. Multiplying the state deviation by this weighting factor yields the strength of the path dependence effect. The strength of the path dependence effect is a comprehensive indicator. The path dependence effect strength reflects both the magnitude of state deviation caused by model prediction bias and the amplifying effect of the aggressiveness of historical control commands on this deviation. A larger path dependence effect strength indicates a stronger cumulative adverse impact of the digital twin model's historical predictions on the current process state, which is difficult to correct through immediate control. This path dependence effect strength value is output and used for decision fusion in subsequent step S5.

[0048] Step S5 is implemented in the following manner: Obtain the quantized prediction confidence and the analyzed path dependence effect strength. The quantized prediction confidence is derived from the calculation output of step S3, and is a value between 0 and 1. The analyzed path dependence effect strength is derived from the calculation output of step S4, and is a non-negative scalar value.

[0049] Based on a pre-defined weight adjustment mapping relationship, the combined input of prediction confidence and path dependence effect strength is mapped to the corresponding model prediction weight adjustment amount. This pre-defined weight adjustment mapping relationship is implemented through a pre-defined fuzzy inference system. The design process of the pre-defined fuzzy inference system is as follows: Two input linguistic variables are defined. The first input linguistic variable is prediction confidence. The universe of discourse for prediction confidence is a continuous interval from 0 to 1. Three fuzzy sets are defined for prediction confidence. The first fuzzy set is low confidence. The membership function of the low confidence fuzzy set uses a trapezoidal function. The parameters of the membership function of the low confidence fuzzy set are determined by analyzing the distribution of prediction confidence quantification values ​​corresponding to severe model prediction inaccuracies in historical data. For example, the membership function of the low confidence fuzzy set takes a value of 1 when the prediction confidence is less than 0.3, takes a value of 0 when the prediction confidence is greater than 0.5, and decreases linearly between 0.3 and 0.5. The second fuzzy set is medium confidence. The membership function of the medium confidence fuzzy set uses a triangular function. The membership function parameters of medium-confidence fuzzy sets are determined by analyzing the quantified values ​​of predicted confidence when the model prediction has general biases. For example, the membership function of a medium-confidence fuzzy set takes a value of 1 when the predicted confidence is 0.4, 0 when the predicted confidence is 0.2, and 0 when the predicted confidence is 0.6, increasing linearly between 0.2 and 0.4, and decreasing linearly between 0.4 and 0.6. The third type of fuzzy set is high-confidence. The membership function of a high-confidence fuzzy set uses a trapezoidal function. The membership function parameters of a high-confidence fuzzy set are determined by analyzing the distribution of quantified values ​​of predicted confidence when the model prediction is relatively accurate. For example, the membership function of a high-confidence fuzzy set takes a value of 1 when the predicted confidence is greater than 0.7, and 0 when the predicted confidence is less than 0.5, increasing linearly between 0.5 and 0.7.

[0050] The second input linguistic variable is the strength of the path dependence effect. The universe of discourse for the strength of the path dependence effect is determined based on historical data analysis, typically within a continuous interval of 0 to 10. Three fuzzy sets are defined for the strength of the path dependence effect. The first fuzzy set is for weak path dependence. The membership function of the weak path dependence fuzzy set is set by analyzing the distribution of quantified values ​​when the path dependence effect is insignificant in historical data. For example, the membership function of the weak path dependence fuzzy set is 1 when the path dependence effect strength is less than 2, 0 when the path dependence effect strength is greater than 4, and linearly decreases between 2 and 4. The second fuzzy set is for moderate path dependence. The membership function of the moderate path dependence fuzzy set is set by analyzing the quantified values ​​under moderate cumulative effect conditions. For example, the membership function of the moderate path dependence fuzzy set is 1 when the path dependence effect strength is 5, 0 when the path dependence effect strength is less than 3, 0 when the path dependence effect strength is greater than 7, linearly increasing between 3 and 5, and linearly decreasing between 5 and 7. The third fuzzy set is for strong path dependence. The membership function of a strong path-dependent fuzzy set is set by analyzing the distribution of quantized values ​​when the path-dependent effect is significant in historical data. For example, the membership function of a strong path-dependent fuzzy set takes a value of 1 when the path-dependent effect strength is greater than 8, takes a value of 0 when the path-dependent effect strength is less than 6, and increases linearly between 6 and 8.

[0051] The output linguistic variable is defined as the model prediction weight adjustment. The universe of discourse for the model prediction weight adjustment is set to a continuous interval from -0.5 to +0.5. Five fuzzy sets are defined for the model prediction weight adjustment: significantly reduced, slightly reduced, unchanged, slightly increased, and significantly increased. The membership function of each fuzzy set is a triangular function. The center position and width of the membership function of each fuzzy set are determined according to the control strategy design objective.

[0052] A fuzzy rule base is established. The fuzzy rule base contains a series of if-then rules. The construction of the fuzzy rule base is based on empirical knowledge of the control logic of the coffee roasting process and analysis of the impact on the error of the digital twin model. The fuzzy rule base contains nine rules. The first rule is: if the prediction confidence is high and the path dependence effect is weak, then the model prediction weight adjustment is slightly increased; the second rule is: if the prediction confidence is high and the path dependence effect is moderate, then the model prediction weight adjustment remains unchanged; the third rule is: if the prediction confidence is high and the path dependence effect is strong, then the model prediction weight adjustment is slightly decreased; the fourth rule is: if the prediction confidence is moderate and the path dependence effect is weak, then the model prediction weight adjustment remains unchanged; the fifth rule is: if the prediction confidence is moderate and the path dependence effect is strong, then the model prediction weight adjustment is slightly decreased; the sixth rule is: if the prediction confidence is moderate and the path dependence effect is weak, then the model prediction weight adjustment is: if the prediction confidence is moderate and the path dependence effect is weak, then the model prediction weight adjustment is: if the prediction confidence is moderate and the path dependence effect is weak, then the model prediction weight adjustment is: if the prediction confidence is moderate and the path dependence effect is strong ... The sixth rule states that if the prediction confidence is medium and the path dependence effect is strong, the model prediction weight adjustment will be significantly reduced. The seventh rule states that if the prediction confidence is low and the path dependence effect is weak, the model prediction weight adjustment will be slightly reduced. The eighth rule states that if the prediction confidence is low and the path dependence effect is medium, the model prediction weight adjustment will be significantly reduced. The ninth rule states that if the prediction confidence is low and the path dependence effect is strong, the model prediction weight adjustment will be significantly reduced.

[0053] The fuzzy inference calculation process is as follows: The prediction confidence and path dependence effect strength of the acquired specific numerical values ​​are used as inputs to the preset fuzzy inference system. Fuzzification calculation is performed, calculating the membership degree of each input value with respect to all its respective fuzzy sets. The fuzzy rule base is activated. For each fuzzy rule, the minimum value of the membership degree of the corresponding input fuzzy set in the antecedent of the rule is taken as the activation strength of the rule. Based on the activation strength of each rule, the Mamdani inference method is used to segment the output fuzzy set specified in the consequent of the rule. The output results of all activated rules are aggregated using the union method. Defuzzification calculation is performed, converting the aggregated output fuzzy set into a precise numerical value as the model prediction weight adjustment amount. The defuzzification calculation uses the centroid method. The centroid method calculates the abscissa value corresponding to the centroid of the area under the membership function curve of the entire output fuzzy set. This abscissa value is the precise model prediction weight adjustment amount obtained after final defuzzification. The entire fuzzy inference calculation process is executed by the preset fuzzy inference system software library.

[0054] The current model prediction weights are updated based on the model prediction weight adjustment amount to generate dynamically adjusted model prediction weights. The current model prediction weight is a stored value, and its initial value is set at system startup. For example, the initial value of the current model prediction weight can be set to 0.8. The update operation adds the model prediction weight adjustment amount to the current model prediction weight. After the update operation, an out-of-bounds check and handling are performed. If the value after the update operation is greater than 1, the dynamically adjusted model prediction weight is set to 1. If the value after the update operation is less than 0, the dynamically adjusted model prediction weight is set to 0. If the value after the update operation is between 0 and 1, this value is directly used as the dynamically adjusted model prediction weight. The dynamically adjusted model prediction weight is output and stored to replace the current model prediction weight for use in step S6 when generating control commands.

[0055] Step S6 is implemented in the following manner: Based on the dynamically adjusted model prediction weights, the prediction results of the digital twin model are weighted and fused to generate a weighted process state prediction. The dynamically adjusted model prediction weights are derived from the output of step S5, and are values ​​between 0 and 1, representing the confidence level of the digital twin model's prediction results in the final decision during the current control cycle. The prediction results of the digital twin model are derived from the prediction data generated by the digital twin model in step S1, which is time-aligned with the multi-dimensional sensor data. Weighted fusion also requires another prediction source as a benchmark. This other prediction source employs a simple predictor based on real-time feedback; this other prediction source is called the feedback controller prediction. The method for generating feedback controller predictions involves first acquiring multi-dimensional sensor data at the current moment and extracting key process state variables from this data, such as coffee bean bed temperature. Then, the target value corresponding to a preset roasting process target state is obtained, for example, the target bean temperature at the current moment. The deviation between the current state and the target state is calculated. Based on this deviation, a preset proportional-integral-derivative (PID) control law is used to calculate the predicted value of the control command applied to the actuator to bring the process state toward the target. The proportional, integral, and derivative coefficients of the PID control law are predetermined using engineering tuning methods. Finally, a simplified process response model is used to deduce the predicted value of the process state at the next moment based on the predicted control command; this predicted value is the feedback controller prediction. Feedback controller prediction aims to reflect the expected state achievable solely through classical feedback control based on the current error, without relying on complex digital twin models. The specific calculation for weighted fusion is as follows. The weighted process state prediction is calculated as: Pweighted = W × Pdt + (1 − W) × Pfb; where Pweighted represents the weighted process state prediction, W represents the dynamically adjusted model prediction weights, Pdt represents the prediction result of the digital twin model, and Pfb represents the feedback controller prediction. The value of W ranges from 0 to 1; for example, W can be set to 0.8. This calculation is performed separately for each key process state variable that needs to be predicted. The weighted process state prediction is a vector that includes a comprehensive prediction of the future process state, integrating the intelligence of the digital twin model with the robustness of classical feedback control.

[0056] Based on the weighted process state prediction and the preset roasting process target state, the target control command is calculated through a control algorithm. The preset roasting process target state is a state trajectory that changes over time. This preset target state is pre-set and stored before roasting begins, based on the coffee bean variety and the desired roasting level. Examples of preset roasting process target states include bean temperature target curves, exhaust temperature target curves, and color development target curves. The control algorithm employs a proportional-integral-derivative (PID) control algorithm. The calculation process of the PID control algorithm is as follows: First, for each controlled key process variable, such as hot air temperature, the predicted value of that variable at the next prediction time is extracted from the weighted process state prediction. The target value of that variable at the same time is obtained from the preset roasting process target state. The deviation between the predicted value and the target value is calculated; this deviation is called the prediction deviation. The PID control algorithm calculates the target control command based on the prediction deviation. The output of the PID control algorithm consists of the sum of three parts: a proportional term, an integral term, and a derivative term. The proportional term is calculated as: proportional term = Kp × prediction deviation; where Kp represents the proportional coefficient. The integral term is calculated as: Integral term = Ki × Integral sum; where Ki represents the integral coefficient, the integral sum is the accumulated value of historical prediction deviations, and an anti-saturation limit is set for the integral sum, which stops accumulating when the integral sum exceeds a preset upper limit. The derivative term is calculated as: Derivative term = Kd × (Current prediction deviation − Previous prediction deviation); where Kd represents the derivative coefficient. The values ​​of the proportional, integral, and derivative terms are added together to obtain the initial control command value. The proportional coefficient Kp, integral coefficient Ki, and derivative coefficient Kd are tuned using engineering tuning methods, such as the Ziegler-Nichols tuning method. The specific tuning process of the Ziegler-Nichols method is as follows: First, set the integral coefficient Ki and derivative coefficient Kd to 0, then gradually increase the proportional coefficient Kp until the control system exhibits constant-amplitude oscillations. Record the critical value of the proportional coefficient Ku and the critical value of the oscillation period Tu at this point. Then, set the parameters according to the Ziegler-Nichols formula. The formula for calculating the proportional coefficient Kp is: Kp = 0.6 × K The integral coefficient Ki is calculated using the formula: Ki = Kp / (0.5 × Tu); the derivative coefficient Kd is calculated using the formula: Kd = Kp × 0.125 × Tu. The empirical coefficients 0.6, 0.5, and 0.125 are recognized empirical values ​​set for standard proportional-integral-derivative controllers in the Ziegler-Nichols method. These values ​​are derived from classic experimental summaries in control engineering and are used to calculate controller parameters that provide good control performance and stability when the system is in a critical oscillating state. In practical applications, these values ​​are then fine-tuned based on the control effect. The calculated preliminary control command values ​​also need to undergo output limiting processing. Output limiting processing sets the upper and lower limits of the command based on the physical operating range of the actuator.For example, for a hot air heater, the physical range of its power control command is 0 to 100, where 0 to 100 represents shutting off to full power. If the initial control command value is greater than the upper limit of 100, the target control command is set to 100; if the initial control command value is less than the lower limit of 0, the target control command is set to 0; if the initial control command value is between 0 and 100, it is directly used as the target control command. For multiple actuators requiring coordinated control, such as the hot air heater and the drum speed motor, the above proportional-integral-derivative (PID) control algorithm is run independently for each actuator to generate its own target control command. The final target control command is a set of commands containing the set values ​​for multiple actuators.

[0057] The target control commands are output to the actuators in the coffee roasting process to drive them and control the roasting process. The actuators include a power controller for the hot air heater, a speed controller for the drum drive motor, and a cooling valve on / off controller. The output process is implemented through a standard industrial control communication protocol, for example, sending the setpoint of the target control command to the local controller of the corresponding actuator via a 4-20 mA analog current signal. Upon receiving the new power percentage setpoint, the power controller for the hot air heater adjusts the conduction angle of the thyristor to change the input power of the heating element, thereby bringing the hot air temperature closer to the value desired by the control algorithm. Upon receiving the speed setpoint, the speed controller for the drum drive motor adjusts the power supply frequency and voltage of the motor via a frequency converter, thereby changing the drum's rotation speed. The cooling valve on / off controller opens or closes the cooling air valve according to the command. The actions of the actuators collectively alter the thermodynamic and hydrodynamic environment within the coffee roaster, thus affecting the coffee bean roasting process. By continuously and periodically executing steps S1 to S6, the entire system achieves dynamic evaluation of model credibility and path dependence effects based on the identification of nonlinear dynamic mutation stages, and intelligent adjustment of control decision weights. Ultimately, it generates precise control commands that integrate the advantages of model prediction and feedback control, thereby achieving closed-loop precise control of the coffee roasting process and improving the stability of the quality of the final roasted product.

[0058] Example 2: Figure 2 A schematic diagram of the digital twin-based precision control system for coffee roasting is provided. The digital twin-based precision control system for coffee roasting includes: The data acquisition module is used to collect multi-dimensional sensor data of the physical entity during the coffee roasting process in real time, and simultaneously acquire the prediction data generated by the digital twin model; The mutation identification module is used to map multi-dimensional sensor data to a low-dimensional feature space and identify whether the coffee roasting process has entered a nonlinear dynamic mutation stage based on the geometric characteristics of the historical motion trajectory of the process state points in the low-dimensional feature space. The credibility assessment module is used to calculate the deviation sequence between the predicted data and the multi-dimensional sensor data for multiple key process variables if the nonlinear dynamic mutation stage is entered. It analyzes the unexpected coupling relationship characteristics that emerge between the deviation sequences in the nonlinear dynamic mutation stage to quantify the prediction credibility of the digital twin model. The effect analysis module is used to analyze the strength of the path dependence effect of the digital twin model caused by the influence of historical predictions, based on the identification results of the nonlinear dynamic mutation stage combined with the historical control command sequence and the current process state. The weight adjustment module is used to dynamically adjust the model prediction weights based on prediction confidence and the strength of path dependence effect; The instruction generation module is used to generate control instructions based on the model prediction weights and the prediction results of the digital twin model to control the coffee roasting process.

[0059] All calculations involved in the embodiments are dimensionless numerical calculations, and the preset parameters and thresholds in the calculations are set by those skilled in the art according to the actual situation.

[0060] It should be noted that this invention can be deployed on the device itself to realize embedded applications, or it can run on a PC or other terminal with a user interface, thereby meeting various hardware environments and usage requirements.

[0061] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, as a computer program product. A computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer programs are loaded or executed on a computer, all or part of the processes or functions according to the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. Computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wireless or wired transmission; wired transmission methods include optical fiber, twisted pair, coaxial cable, etc.; wireless transmission includes infrared, microwave, etc. Computer-readable storage media can be any available medium that a computer can access or a data storage device such as a server or data center that contains one or more sets of available media. Available media can be magnetic media (e.g., floppy disks, hard disks, magnetic tapes), optical media (e.g., DVDs), or semiconductor media. Semiconductor media can be solid-state drives.

[0062] Those skilled in the art will understand that, for the sake of convenience and brevity, the specific working processes of the systems, devices, and modules described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0063] In the several embodiments provided in this application, it should be understood that the disclosed systems, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of modules is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple modules or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between apparatuses or modules may be electrical, mechanical, or other forms.

[0064] The modules described as separate components may or may not be physically separate. The components shown as modules may or may not be physical modules; they may be located in one place or distributed across multiple network modules. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs.

[0065] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0066] If a function is implemented as a software module and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods in the various embodiments of this application. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0067] The above are merely specific embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0068] In conclusion, the above are merely preferred embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for precise control of the coffee roasting process based on digital twins, characterized in that, include: S1: Real-time acquisition of multi-dimensional sensor data of the physical entity during coffee roasting, and simultaneous acquisition of prediction data generated by the digital twin model; S2: Map multi-dimensional sensor data to a low-dimensional feature space, and identify whether the coffee roasting process has entered a nonlinear dynamic mutation stage based on the geometric characteristics of the historical motion trajectory of the process state points in the low-dimensional feature space. S3: If the nonlinear dynamic mutation stage is entered, the deviation sequence between the predicted data and the multi-dimensional sensor data is calculated for multiple key process variables. The unexpected coupling relationship characteristics that emerge between the deviation sequences in the nonlinear dynamic mutation stage are analyzed to quantify the prediction credibility of the digital twin model. S4: Based on the identification results of the nonlinear dynamic mutation stage, combined with the historical control command sequence and the current process state, the strength of the path dependence effect of the digital twin model due to the influence of historical prediction is analyzed. S5: Dynamically adjust the model prediction weights based on prediction confidence and the strength of path dependence effect; S6: Based on the model prediction weights and the prediction results of the digital twin model, control instructions are generated to control the coffee roasting process.

2. The method for precise control of coffee roasting process based on digital twin as described in claim 1, characterized in that, S1 includes: Real-time acquisition of multi-dimensional sensor data reflecting the thermodynamic and chemical states of coffee roasting process; Input multi-dimensional sensor data into the digital twin model; Obtain predictive data that is temporally aligned with the multidimensional sensor data, generated by forward simulation calculations based on multidimensional sensor data and the internal state of the model, from the digital twin model.

3. The method for precise control of coffee roasting process based on digital twin according to claim 1, characterized in that, S2 include: Dimensionality reduction processing is performed on multi-dimensional sensor data to obtain process state points represented in a low-dimensional feature space; By connecting continuous process state points in time sequence, the historical motion trajectory of process state points in a low-dimensional feature space is formed. Calculate the changes in geometric characteristics of historical motion trajectories within the current and adjacent time windows; Determine whether the change in geometric characteristics continuously exceeds the preset abrupt change threshold and reaches the preset duration; If so, the coffee roasting process is determined to have entered a non-linear dynamic mutation stage.

4. The method for precise control of coffee roasting process based on digital twin according to claim 1, characterized in that, S3 include: For each key process variable, within the time window of the nonlinear dynamic mutation stage, calculate the time-ordered sequence of deviations between the predicted data and the multi-dimensional sensor data; Calculate the coupling strength measure between each pair of deviation sequences to obtain the coupling relationship matrix in the nonlinear dynamic mutation stage; Obtain the baseline coupling relationship matrix established in advance during the steady-state phase of the process; Calculate the difference measure between the coupling relationship matrix and the baseline coupling relationship matrix in the nonlinear dynamic mutation stage, and use the difference measure to characterize the features of unexpected coupling relationships; Based on the saliency of the unexpected coupling relationship features, a predictive reliability quantification of the digital twin model is generated by mapping.

5. The method for precise control of coffee roasting process based on digital twin according to claim 4, characterized in that, The coupling strength measure between each pair of deviation sequences is calculated, including: for each pair of deviation sequences corresponding to key process variables, within the time window of the nonlinear dynamic mutation stage, the maximum information coefficient between the two deviation sequences is calculated, and the value of the maximum information coefficient is used as the coupling strength measure between the deviation sequences corresponding to this pair of key process variables.

6. The method for precise control of coffee roasting process based on digital twin according to claim 1, characterized in that, S4 include: Obtain the historical control command sequence and current process status within the time period preceding the entry into the nonlinear dynamic mutation stage; By combining the historical control command sequence with the reference model of the process's steady phase, a reference process state sequence that is expected to be reached if there is no model prediction bias is generated. Calculate the deviation between the current process state and the expected state of the reference process state sequence at the same time. The path dependence effect strength is obtained by weighting and correcting the state deviation based on the radicalness of the command changes in the historical control command sequence.

7. The method for precise control of coffee roasting process based on digital twin according to claim 1, characterized in that, S5 include: Obtain the quantified prediction credibility and the analyzed path dependence effect strength; Based on the preset weight adjustment mapping relationship, the combined input of prediction confidence and path dependence effect strength is mapped to the corresponding model prediction weight adjustment amount; The current model prediction weights are updated based on the model prediction weight adjustment amount to generate dynamically adjusted model prediction weights.

8. The method for precise control of coffee roasting process based on digital twin according to claim 7, characterized in that, Based on the preset weight adjustment mapping relationship, the combined input of prediction confidence and path dependence effect strength is mapped to the corresponding model prediction weight adjustment amount, including: The prediction confidence and the strength of the path dependence effect are input into a pre-defined fuzzy inference system; Through the inference calculation of the fuzzy inference system, the model prediction weight adjustment amount corresponding to the combination of prediction credibility and path dependence effect strength is output.

9. The method for precise control of coffee roasting process based on digital twin according to claim 1, characterized in that, S6 include: Based on the dynamically adjusted model prediction weights, the prediction results of the digital twin model are weighted and fused to generate a weighted process state prediction. Based on the weighted process state prediction and the preset baking process target state, the target control command is calculated through the control algorithm; The target control command is output to the actuator of the coffee roasting process to drive the actuator to control the coffee roasting process.

10. A precise control system for a coffee roasting process based on digital twins, used to implement the precise control method for a coffee roasting process based on digital twins as described in any one of claims 1-9, characterized in that, include: The data acquisition module is used to collect multi-dimensional sensor data of the physical entity during the coffee roasting process in real time, and simultaneously acquire the prediction data generated by the digital twin model; The mutation identification module is used to map multi-dimensional sensor data to a low-dimensional feature space and identify whether the coffee roasting process has entered a nonlinear dynamic mutation stage based on the geometric characteristics of the historical motion trajectory of the process state points in the low-dimensional feature space. The credibility assessment module is used to calculate the deviation sequence between the predicted data and the multi-dimensional sensor data for multiple key process variables if the nonlinear dynamic mutation stage is entered. It analyzes the unexpected coupling relationship characteristics that emerge between the deviation sequences in the nonlinear dynamic mutation stage to quantify the prediction credibility of the digital twin model. The effect analysis module is used to analyze the strength of the path dependence effect of the digital twin model caused by the influence of historical predictions, based on the identification results of the nonlinear dynamic mutation stage combined with the historical control command sequence and the current process state. The weight adjustment module is used to dynamically adjust the model prediction weights based on prediction confidence and the strength of path dependence effect; The instruction generation module is used to generate control instructions based on the model prediction weights and the prediction results of the digital twin model to control the coffee roasting process.

Images