A method for calibrating parameters of a temperature field model of a grain bin

Abstract

The application provides a kind of calibration method of granary temperature field model parameters, belongs to the technical field of grain condition measurement and control and numerical simulation.The method comprises the following steps: collecting and preprocessing granary environment and state data to obtain granary temperature field data, and constructing a three-dimensional temperature field mechanism model of the granary based on a three-dimensional unsteady heat conduction equation.A phased elite filing calibration method is used to obtain the optimal parameter combination, and the combination is subjected to multidimensional accuracy verification and physical consistency test to confirm the effectiveness of the simulation result.The application takes the three-dimensional unsteady heat conduction equation as the core, constructs a mechanism model containing boundary heat exchange and layered heat source, improves the parameter identification efficiency by using a phased search strategy, combines the search efficiency of machine learning and the explanation ability of mechanism model, and avoids the blindness of traditional calibration.Provides a general solution for the calibration of mechanism model in sparse observation scenarios, has physical interpretability and engineering usability, and provides a scientific basis for the safety management of grain storage.

Description

Technical Field

[0001] This invention belongs to the field of grain condition monitoring and numerical simulation technology, specifically a calibration method for the parameters of a grain storage temperature field model. Background Technology

[0002] The temperature field in grain storage is a complex three-dimensional spatiotemporal process closely coupled with grain respiration, moisture migration, mold growth, and quality deterioration. Traditional single-point or average temperature monitoring is insufficient to reveal the formation and evolution of localized hot spots, which are often precursors to grain quality deterioration. Therefore, constructing a high-precision three-dimensional temperature field simulation model is of great significance for ensuring the safety of grain storage. Internationally, research on grain storage temperature fields has evolved from early simplified one-dimensional or two-dimensional models to a three-dimensional multiphysics framework capable of coupling heat, humidity, airflow, and gas component migration. The advantage of mechanistic models lies in their clear physical meaning and strong interpretability, but they require dense and stable observational data to constrain the parameters.

[0003] With the development of grain condition monitoring systems, machine learning methods have been widely applied in the field of grain temperature prediction. Researchers have achieved excellent fitting results in short-term grain temperature prediction by utilizing advanced algorithms such as Support Vector Regression (SVR), ensemble learning, Long Short-Term Memory (LSTM) networks, Temporal Convolutional Networks (TCN), and even Transformers. These methods represent different technical directions such as traditional machine learning, ensemble learning, temporal deep learning, and multi-scale feature extraction, and their theoretical foundations can be traced back to relevant classic research.

[0004] However, the "black box" nature of machine learning models makes it difficult to meet the interpretive needs of warehouse management. Especially in real-world scenarios with sparse observations, pure machine learning models are prone to learning sensor noise, equipment malfunctions, or random events as "patterns," leading to incorrect decision signals in practical applications. Related research emphasizes outlier diagnosis, interpolation completion, and signal complexity analysis in the data preprocessing stage, such as using permutation entropy to identify sequence mutations. This work reminds us that model performance is not solely determined by the algorithm; the quality of the observation chain also profoundly affects the final results.

[0005] In real-world grain storage scenarios, observational data often suffers from difficulties in acquisition, sparsity, and limited quality, leading to unconstrained traditional mechanism calibration and a lack of physical interpretability in pure machine learning. In existing technologies, the elite mixing in stage two generates coarse-grained elite interference, and the elite crossover in stage three results in ineffective crossover operations. This causes the RMSE curve to exhibit abnormal and repeated fluctuations in the later stages, failing to stabilize at the minimum point. High-cost operations (crossover) are used to replace low-cost operations (perturbation) without any additional benefits, resulting in low computational efficiency. Summary of the Invention

[0006] This invention provides a calibration method for the parameters of a grain storage temperature field model, in order to overcome the deficiencies in the prior art.

[0007] This invention provides a method for calibrating parameters of a grain storage temperature field model, comprising: The environmental and state data of the grain warehouse were collected and preprocessed to obtain the temperature field data of the grain warehouse. Based on the three-dimensional unsteady heat conduction equation, a three-dimensional temperature field mechanism model of the grain warehouse was constructed.

[0008] Implement a phased elite archiving rate determination method, including: Phase 1: Using random sampling and extreme perturbations, a large-scale search is conducted to escape local optima and establish an initial elite pool.

[0009] Phase 2: Employ a hybrid elite approach with small perturbations to perform a medium-range search, followed by a fine-grained search around the elite solution to quickly converge to the region of minimum error and reduce RMSE.

[0010] Phase 3: Employ focused perturbation combined with elite cross-referencing to conduct a small-scale search, optimizing behavioral constraints while maintaining accuracy.

[0011] After phased elite archiving and calibration of the three-dimensional temperature field mechanism model of the grain warehouse, the optimal parameter combination was obtained. The combination was then subjected to multi-dimensional accuracy verification and physical consistency test to confirm the effectiveness of the simulation results.

[0012] According to the present invention, a method for calibrating parameters of a grain storage temperature field model is provided, wherein the elite mixing process includes: Two elite samples are selected from the elite sample library as parents. The difference vector between the two parents is calculated, and a new combination of candidate parameters is generated along the direction of the difference vector.

[0013] The elite sample library obtained in Phase 1 is updated to generate new candidate parameter combinations. The average sliding window value of the elite mixed entry rate is analyzed to determine the coarse-grained elite interference execution strategy.

[0014] According to the calibration method for grain warehouse temperature field model parameters provided by the present invention, the process of analyzing the sliding window average value of the elite mixed storage rate includes: Before the start of Phase Two, set a starting value for the elite mixed-input round.

[0015] In Phase Two, the elite mixing operation is performed. When the number of elite mixing entry rounds reaches the start value of the elite mixing entry round, the sliding window length is calculated, and the elite mixing entry judgment value is calculated.

[0016] Based on the sliding window length and the elite mixed entry judgment value, the average sliding window value of the elite mixed entry rate is obtained.

[0017] According to the calibration method for grain warehouse temperature field model parameters provided by the present invention, the process of determining the coarse-grained elite interference execution strategy includes: Obtain the average threshold of the sliding window for the pre-defined elite mixed entry rate in the database.

[0018] If the average sliding window value of the elite mixed entry rate is less than the average threshold of the elite mixed entry rate, then the coarse-grained elite interference execution strategy is recorded as continuing local refinement.

[0019] If the average value of the sliding window for the elite mixed entry rate is greater than or equal to the average threshold of the sliding window for the elite mixed entry rate, then the coarse-grained elite interference execution strategy is recorded as performing elite interference cleaning and adjustment.

[0020] When making a judgment, intermittent detection is performed.

[0021] Perform elite interference cleaning and adjustment.

[0022] According to the calibration method for grain warehouse temperature field model parameters provided by the present invention, the process of performing intermittent detection includes: The system status is used to determine the number of intermittent detection rounds of samples.

[0023] The sample intermittent detection round value is obtained by analyzing the grain warehouse temperature field data acquired by the system at this time.

[0024] Based on the intermittent detection round value, an intermittent detection strategy is executed to determine coarse-grained elite interference.

[0025] According to the calibration method for grain storage temperature field model parameters provided by the present invention, the process of elite disturbance cleaning and adjustment includes: Obtain the timestamp of the samples added to the elite library.

[0026] Get the current iteration round number and the total number of iteration rounds in Phase 2 for the sample.

[0027] The difference in the number of iteration rounds is obtained based on the current iteration round number and the total number of iteration rounds in Phase 2.

[0028] Obtain the preset elimination ratio coefficient.

[0029] The dynamic elimination value is obtained by using the difference in the number of iteration rounds, the total number of iteration rounds, and the elimination ratio coefficient.

[0030] Elite interference cleaning and adjustment are performed based on dynamic elimination values.

[0031] According to the calibration method for grain warehouse temperature field model parameters provided by the present invention, the elite cross-operation process includes: Two elite samples are selected from the elite sample library as parents. Mixing coefficients are independently assigned to each dimension in the parameter space. The parent parameter values ​​of each dimension are weighted and combined dimension by dimension to generate new candidate parameter combinations.

[0032] While maintaining the sample accuracy of Phase 2, perform final parameter fitting, identify invalid crossover strategies, and adjust invalid crossover accordingly.

[0033] According to the calibration method for grain warehouse temperature field model parameters provided by the present invention, the process of determining invalid crossover strategies includes: Obtain the current perturbation amplitude of the sample.

[0034] The average distance between two samples in the elite database in the normalized parameter space is denoted as the sample average distance.

[0035] If the average distance between samples is greater than or equal to the current perturbation magnitude, the invalid crossover strategy is recorded as continuing behavior correction.

[0036] If the average distance between samples is less than the current perturbation amplitude, then the invalid crossover strategy is recorded as invalid crossover adjustment.

[0037] According to the calibration method for grain warehouse temperature field model parameters provided by the present invention, the process of performing ineffective cross-adjustment includes: The ratio of the average distance of the samples to the current perturbation amplitude is denoted as the elite cross-dynamic ratio.

[0038] Obtain the minimum and maximum execution probabilities of elite crossover.

[0039] Obtain the judgment coefficient for invalid crossover.

[0040] Ineffective cross-adjustment is performed based on the data obtained above.

[0041] This invention provides a calibration method for grain storage temperature field model parameters. It employs a Bayesian-inspired phased elite archiving calibration method, drawing on its 'prior-evidence-posterior' update logic. A phased elite search strategy replaces the expensive surrogate model, and innovatively introduces a 'behavioral penalty term' into the objective function, using the physical laws governing grain storage temperature changes (such as heat transfer direction and respiration heat variation trends) as constraints. This method aims to solve the under-constrained parameter problem under sparse data with lower computational cost, achieving high-precision and physically consistent three-dimensional temperature field simulations. The phased search strategy improves parameter identification efficiency and avoids the blindness of traditional calibration. This method integrates machine learning search efficiency with mechanistic model interpretation capabilities, providing a general solution for calibrating mechanistic models in sparse observation scenarios, and has significant practical application value for ensuring the safety of grain storage.

[0042] The elite hybrid optimization strategy in Stage 2 eliminates the impact of differences in the physical parameters of the grain storage on the search direction, solves the diversity depletion problem caused by inbreeding in the elite pool, and achieves a smooth transition from coarse-grained elites to fine-grained elites, effectively improving the convergence speed and final accuracy of the local refining stage. The elite crossover optimization strategy in Stage 3 adaptively matches the crossover execution probability with the clustering degree of the elite pool, aligns the parent selection criteria with the behavior correction target, and correlates the perturbation direction with the penalty gradient direction, effectively solving problems such as excessively high invalid crossover ratios, deviation in correction direction, and mismatched perturbation step size. Attached Figure Description

[0043] The invention will now be further described with reference to the accompanying drawings.

[0044] Figure 1 This is a flowchart illustrating a method for calibrating parameters of a grain storage temperature field model in this invention. Figure 2 This is a flowchart illustrating the phased elite archiving rate determination method in this invention. Figure 3 This is a flowchart illustrating the process of determining the coarse-grained elite interference execution strategy in this invention. Figure 4 This is a schematic diagram of the coordinate arrangement of the grain temperature sensing network in this invention; Figure 5 This is a schematic diagram of the time sequence comparison of the stratified temperature in the grain warehouse in this invention; Figure 6 This is a schematic diagram comparing the planar distribution of the temperature field in the grain warehouse in this invention; Figure 7 This is a schematic diagram of the vertical temperature profile of the grain silo in this invention; Figure 8 This is a schematic diagram of the boundary driving and overall temperature equalization timing in this invention; Figure 9 This is a schematic diagram of the parameter calibration optimization trajectory in this invention. Detailed Implementation

[0045] To make the technical means, creative features, objectives and effects of this invention easier to understand, the invention will be further described below in conjunction with specific embodiments.

[0046] like Figures 1 to 4 As shown in the figure, an embodiment of the present invention provides a method for calibrating the parameters of a grain storage temperature field model, comprising: refer to Figure 4This invention uses a specific grain warehouse as the research object. The warehouse measures 24.0m in length, 18.0m in width, and 8.0m in height. The observation period spans from July 18th to November 7th of a given year, covering the critical transition from summer heat to autumn cooling. Four layers of temperature-measuring cables are installed inside the warehouse, with 24 measuring points in each layer, forming a 4×6×4 three-dimensional observation network. A total of 1152 temperature observations were collected at 12 time intervals. The coordinates of the grain temperature sensor network layout are shown in [reference needed]. Figure 2 It is displayed in the middle.

[0047] The environmental and state data of the grain warehouse were collected and preprocessed to obtain the temperature field data of the grain warehouse. Based on the three-dimensional unsteady heat conduction equation, a three-dimensional temperature field mechanism model of the grain warehouse was constructed.

[0048] Data preprocessing followed the principle of "preserving physical events and removing equipment anomalies." Suspicious cold spots significantly below the layer mean and outliers with abnormal spatial variance within the layer were identified and removed. The removal rule was that observations causing the variance of the air temperature in that layer to be greater than 2.5℃ were used. Short-term missing air temperature data within the warehouse were completed using time interpolation. After preprocessing, 1064 valid observation data points were obtained, providing a reliable data foundation for the parameter calibration of the three-dimensional temperature mechanism model of the grain warehouse.

[0049] The process of constructing a three-dimensional temperature field mechanism model for grain storage includes: S1, Construct the governing equations: The grain pile is considered an anisotropic porous medium, and its temperature field satisfies the three-dimensional unsteady-state heat conduction equation: , In the formula, T represents the grain pile temperature (°C), which is the three-dimensional temperature field to be solved. t represents time (s). x, y, and z represent spatial coordinates (m), corresponding to the length, width, and height of the grain pile, respectively. αh represents the horizontal thermal diffusivity (m² / s), which controls the temperature diffusivity within the same layer. αv represents the vertical thermal diffusivity (m² / s), which controls the interlayer penetration rate of cold and heat. Q(z,t,T) represents the internal heat source term (W / m³), characterizing the internal thermal support at different depths and stages. ρc represents the volumetric heat capacity (J / (m³·K), converting the heat source intensity into the rate of temperature change.

[0050] The model employs an explicit finite difference method for time integration, with a basic output step size set to 1 hour (3600 seconds). To ensure numerical stability, each output step size is automatically decomposed into multiple internal computational substeps based on the Courant-Friedrichs-Lewy (CFL) conditions of the heat conduction equation.

[0051] The maximum stable time step is determined by both the mesh size and the thermal diffusivity, and is calculated using the following formula: , In the formula, The maximum stable time step (s). , These are the horizontal and vertical thermal diffusivity (m² / s), respectively. , , They are respectively , , Grid spacing in the direction (m). , , They are respectively Directional sidewalls Equivalent heat transfer coefficients (dimensionless) for the directional sidewalls and bottom boundary.

[0052] In the actual simulation, each 1-hour output step is decomposed into several internal computational substeps, the number of which is... The calculation formula is: , In the formula, Output step size (3600 seconds). This indicates rounding up. This adaptive time step strategy ensures both numerical stability and the ability to accurately capture transient changes in the temperature field.

[0053] In the simulation cycle, boundary conditions (such as outdoor temperature and indoor temperature) are linearly interpolated in each internal substep to ensure the continuity of boundary-driven operations. The observation data time interval is 10–22 days, and the model integrates continuously with a step size of 1 hour. Simulation results are extracted at the observation time for comparison and verification.

[0054] S2, Determine boundary conditions: Boundary conditions act as a bridge connecting the internal and external environments of the grain pile, determining how external heat enters and how internal heat dissipates. This study establishes explicit flux equations for the top, sidewalls, and floor boundaries.

[0055] The heat flux at the top boundary consists of two parts: convective heat transfer between the air inside the silo and the grain temperature at the top, and shortwave solar radiation absorbed by the roof. The calculation formula is as follows: , In the formula, The heat flux at the top boundary (W / m²) represents the intensity of heat exchange between the top boundary and the grain pile. The equivalent heat transfer coefficient at the top (W / (m²·K)) represents the strength of the heat transfer from the top boundary into the grain pile. The temperature inside the chamber (°C) represents the ambient temperature inside the chamber controlled by the top boundary. The temperature of the top layer grain (°C) represents the temperature obtained by the numerical model on the top boundary element. The equivalent shortwave input at the top (W / m²) represents the contribution of roof radiative heat storage to the temperature of the upper grain layer.

[0056] The sidewall boundary considers the convective heat transfer between the outside air temperature conducted through the wall and the edge grain layer. The calculation formula is as follows: , In the formula, The heat flux at the sidewall boundary (W / m²) represents the intensity of heat exchange between the sidewall boundary and the grain pile. The equivalent heat transfer coefficient of the sidewall (W / (m²·K)) represents the strength of the heat transfer from the sidewall boundary into the grain pile. The equivalent temperature of the wall (°C) represents the equivalent temperature of the external air temperature driven by the control sidewall boundary after it is conducted through the wall. The edge grain temperature (°C) represents the temperature obtained by the numerical model on the edge boundary element.

[0057] The floor boundary considers the combined effects of soil thermal inertia and the bottom insulation layer; the calculation formula is as follows: , In the formula, The heat flux at the floor boundary (W / m²) represents the intensity of heat exchange between the floor boundary and the grain pile. The equivalent heat transfer coefficient of the floor (W / (m²·K)) represents the strength of the heat transfer from the floor boundary into the grain pile. Ground temperature (°C) represents the soil temperature driven by the bottom boundary. The bottom grain temperature (°C) represents the temperature obtained by the numerical model on the bottom boundary unit.

[0058] The annual cycle fluctuation of geothermal temperature is calculated as follows: , In the formula, The soil temperature (°C) varies over time, representing the soil temperature boundary condition considering annual cycle fluctuations. The base value for geothermal temperature (°C) represents the annual average temperature of the deep soil layers. The amplitude of the geothermal temperature (°C) represents the annual fluctuation range of the geothermal temperature. Let d be time, representing the number of days in a year. The geothermal phase lag (d) represents the number of days that the geothermal peak lags behind the air temperature peak.

[0059] The formula for calculating the equivalent heat transfer coefficient after considering the thermal resistance of the insulation layer and the soil is: , In the formula, The floor's equivalent heat transfer coefficient (W / (m²·K)) represents the overall heat transfer capacity considering the thermal resistance of the insulation layer and the soil. The floor foundation heat transfer coefficient (W / (m²·K)) represents the convective heat transfer capacity between the earth's surface and the air. The thickness (m) of the bottom insulation layer represents the geometric dimensions of the insulation layer that control the bottom insulation performance. The thermal conductivity of the bottom insulation layer (W / (m·K)) represents the thermal conductivity of the insulation material. The soil representative depth (m) represents the characteristic depth that determines the magnitude of soil thermal inertia. Soil thermal conductivity (W / (m·K)) represents the soil's ability to conduct heat.

[0060] S3, Construct a layered heat source model: To accurately depict the complex thermal support structure inside the grain pile, the total heat source is decomposed into three parts: a basic uniform heat source, a breathing heat source, and a stratified Gaussian heat source. , In the formula, The total heat source term (W / m³) represents the sum of internal volume heat sources entering the governing equations. The background heat source is represented by a uniform heat source (W / m³), which is controlled by a constant term and a temperature coefficient. The heat source for respiration (W / m³) represents the internal heat source generated by the respiration of grains. The term "layered heat source" (W / m³) indicates the layered thermal support that exists inside the grain pile after the membrane is removed.

[0061] (1) Basic uniform heat source: The formula for calculating a basic uniform heat source is: , In the formula, The background heat source is represented by a uniform heat source (W / m³) and the term represents the background heat source term as it varies with temperature. The basic heat source constant (W / m³) represents the basic heat production intensity that does not change with temperature. The heat source temperature coefficient (W / (m³·K)) represents the sensitivity of the basic heat source to temperature changes. The temperature of the grain pile (°C) represents the actual temperature of the grain pile at present. The reference temperature (°C) represents the base temperature used to calculate the basic heat source.

[0062] In this study, both Q0 and b_Q are fixed at 0 to reduce unnecessary degrees of freedom.

[0063] (2) Layered Gaussian heat source: The grain pile exhibits distinct layered thermal support, represented by a superposition of three Gaussian distributions: , In the formula, The heat source for the stratification (W / m³) represents the strength of the stratified thermal support inside the grain pile after the membrane is removed. Establish a time modulation function for the heat source, representing a dimensionless factor that controls the speed at which the heat source builds up after the film is removed. Let be the heat source decay time modulation function, representing the dimensionless factor controlling the later decay process of the middle layer heat source. The values ​​represent the heat source intensities (W / m³) for the shallow, upper-middle, and lower-middle layers, respectively, indicating the peak parameters controlling the thermal support intensity at different layers. These are Gaussian profiles of the shallow, upper-middle, and lower-middle layers, representing the spatial distribution of each heat source in the vertical direction.

[0064] The formula for calculating the Gaussian profile of the upper heat source space is: , In the formula, This is a shallow spatial Gaussian profile, representing the spatial distribution of shallow heat sources in the vertical direction. The coordinates are spatial height (m), representing the vertical position within the grain pile. The warehouse height (m) represents the geometric parameters used to normalize the height coordinates. The location of the center of the shallow heat source indicates the center position of the shallow heat source at the normalized height. The thickness of the shallow heat source represents the spatial distribution width of the shallow heat source at the normalized height.

[0065] The formula for calculating the Gaussian profile of the upper and middle heat source space is: , In the formula, This is a Gaussian profile of the upper and middle layers, representing the spatial distribution of the main supporting heat source in the upper and middle layers in the vertical direction. The coordinates are spatial height (m), representing the vertical position within the grain pile. The warehouse height (m) represents the geometric parameters used to normalize the height coordinates. The center position of the upper-middle layer heat source indicates the center position of the upper-middle layer heat source at the normalized height. The thickness of the upper-middle layer heat source represents the spatial distribution width of the upper-middle layer heat source at the normalized height.

[0066] The formula for calculating the Gaussian profile of the middle and lower heat source space is: , In the formula, The Gaussian profile of the middle and lower layers represents the spatial distribution of the supporting heat sources in the vertical direction. The coordinates are spatial height (m), representing the vertical position within the grain pile. The warehouse height (m) represents the geometric parameters used to normalize the height coordinates. The location of the center of the middle and lower layer heat source indicates the center position of the middle and lower layer heat source at the normalized height. The thickness of the lower-middle layer heat source represents the spatial distribution width of the lower-middle layer heat source at the normalized height.

[0067] (3) Time modulation function: The heat source settling time modulation function is expressed as: , In the formula, Establish a time modulation function for the heat source, representing a dimensionless factor that controls the speed at which the heat source builds up after the film is removed. Time (d) represents the number of days after the film is removed. The establishment timescale (d) for the heat source represents the characteristic time it takes for the heat source to establish itself from zero to a steady state.

[0068] The heat source decay time modulation function is expressed as: , In the formula, Let f(x) be the heat source decay time modulation function, which represents the dimensionless factor controlling the later decay process of the middle layer heat source. Time (d) represents the number of days after the film is removed. The decay start time (d) represents the time point at which the mid-layer heat source begins to de-heat. The attenuation floor value (—) represents the proportion of residual intensity after the heat source has attenuated. The decay timescale (d) represents the characteristic time that determines the rate of heat reduction in the later stages.

[0069] Implement a phased elite archiving rate determination method.

[0070] This method framework includes: Drawing inspiration from Bayesian calibration, parameter identification is achieved through a "prior-evidence-posterior" update process. However, instead of employing a computationally expensive Gaussian process surrogate model, it utilizes a phased elite archiving strategy to achieve efficient search. The core idea is: using a mechanistic model as the framework, observation errors as the supervisory signal, and phased search and elite sample memorization to drive parameter optimization. The core of Bayesian calibration can be described as follows: , In the formula, This is the posterior distribution, representing the updated understanding of the parameters based on observational evidence; This is the model parameter vector, representing the set of physical parameters to be identified. The observed data represents evidence used to update the confidence level of the parameters. Let be the likelihood function, representing the likelihood given parameters. The following observation data The probability of occurrence. The prior distribution of the parameters represents the understanding of the reasonable range and central location of the parameters before calibration begins.

[0071] The calibration method proposed in this paper logically follows the core idea of ​​Bayesian inference, but improves the implementation path for engineering applications: 1) Prior distribution: By setting the physical feasible region and initial distribution of the parameters, the prior knowledge of the parameters is reflected; 2) Evidence update: Unlike traditional methods that construct complex likelihood functions, this paper constructs a comprehensive scoring function that includes a precision term and a behavior penalty term. The behavior penalty term acts as a 'physical regularization' term, effectively eliminating parameter combinations that, although having small fitting errors, violate thermodynamic laws; 3) Posterior search: Abandoning the computationally inefficient traditional sampling method, a phased elite archiving strategy is adopted. This "Bayesian-inspired" strategy essentially deeply integrates the search efficiency of machine learning with the physical interpretability of the mechanistic model.

[0072] The design of the objective function of this method includes: Traditional calibration methods only use RMSE as the objective function, which can easily lead to physical distortion under conditions of sparse observations. This invention constructs a comprehensive scoring function that includes both a precision term and a behavior penalty term: , In the formula, Root mean square error (°C) is the basic metric representing the accuracy of model simulation. The total number of observation samples represents the number of observation data points involved in error calculation. For the first The simulated temperature (°C) at each point represents the predicted temperature value output by the numerical model. For the first The observed temperature (°C) at each point represents the actual temperature value measured by the sensor.

[0073] The overall score function equation is as follows: , In the formula, Let be the comprehensive scoring function (°C), representing the comprehensive optimization objective that includes both accuracy and behavioral penalty terms. Root mean square error (°C) represents the fundamental simulation error of the model. The upper-level cold penalty item (°C) represents a behavioral constraint to prevent the upper level from remaining cold for an extended period. The penalty (°C) for excessively rapid cooling of the upper layer represents a behavioral constraint to prevent abnormal cooling rates in the upper layer. The penalty (°C) for insufficient cooling of the mid-layer indicates a behavioral constraint to prevent insufficient cooling of the mid-layer. The penalty (°C) for warmer lower layers represents a behavioral constraint to prevent prolonged warmer lower layers.

[0074] By introducing behavioral penalties, the model is forced to conform to the physical laws governing temperature changes in the grain silo, avoiding unreasonable phenomena such as excessively cold upper layers, insufficient cooling in the middle layers, and excessively warm lower layers. This method is essentially a practical application of the objective function shaping concept from machine learning in the parameter tuning of mechanistic models.

[0075] A phased search strategy includes: To address the challenge of searching high-dimensional parameter spaces, this invention employs a "phased elite archiving" strategy. The entire calibration process is divided into three stages: S1, Phase 1, namely global exploration, uses random sampling and extreme perturbation to conduct a large-scale search, escape local optima, and establish an initial elite pool.

[0076] Global exploration refers to performing random sampling and extreme perturbation operations within a preset feasible region of parameters to conduct a large-scale global search, traversing various regions of the feasible region of parameters, and establishing an initial elite sample library. The initial elite sample library contains multiple sets of candidate parameter combinations that perform well in the evaluation of the objective function.

[0077] S2, the second stage, is local refinement. It uses a mixture of elite solutions with small perturbations to perform a medium-range search, and then performs a fine search around the elite solution to quickly converge to the region of minimum error and reduce RMSE.

[0078] Local refinement refers to: taking the high-quality solutions in the elite sample library as the center, performing elite hybrid operation and small perturbation operation, and conducting a moderate-range local fine search within the parameter space spanned by the elite sample library, so that the root mean square error of the model output quickly converges to the minimum region.

[0079] Two elite samples are selected from the elite sample library as parents. The difference vector between the two parents is calculated, and a new combination of candidate parameters is generated along the direction of the difference vector.

[0080] The elite sample library obtained in Phase 1 is updated to generate new candidate parameter combinations. The average sliding window value of the elite mixed entry rate is analyzed to determine the coarse-grained elite interference execution strategy.

[0081] Before the start of Phase Two, set a starting value for the elite mixed-input round.

[0082] The "elite blending and warehousing round start value" refers to the value of N rounds of elite blending iterations that begin monitoring the elite blending and warehousing rate at the start of Phase Two. This value is obtained from the grain warehouse temperature field data acquired by the system.

[0083] In Phase Two, the elite mixing operation is performed. When the number of elite mixing entry rounds reaches the start value of the elite mixing entry round, the sliding window length is calculated, and the elite mixing entry judgment value is calculated.

[0084] The sliding window length is a preset positive integer constant representing the "number of historical iterations" used for statistics. That is, it averages the elite mixed-inclusion judgment values ​​from the most recent t iterations to reflect the "recent trend of the proportion of elite mixed-inclusion." Using a sliding window avoids interference from single outliers, reflects recent trends rather than global history, and allows statistics to track real-time changes in the algorithm's state.

[0085] The elite pooling judgment value is a binary indicator variable (i.e., 0 or 1, denoted by T), used to record which operation generated the candidate solutions newly entering the elite pool in the t-th iteration. If T=1, it means that the elite samples newly entered in the t-th iteration were generated by the elite pooling operation; if T=0, it means that the elite samples newly entered in the t-th iteration were generated by the small-step perturbation operation. This value is automatically recorded in the "elite pool update" phase of each iteration.

[0086] The elite pool update record refers to the following: each new sample is appended with a label field, indicating whether the sample was generated by "elite mixing" or directly by "small-step perturbation". (Note: This label is only used for statistics and belongs to the algorithm's internal log record variable).

[0087] Based on the sliding window length and the elite mixed entry judgment value, the average sliding window value of the elite mixed entry rate is obtained.

[0088] The sliding window average of the elite pooling inclusion rate is the proportion of elite pooling inclusions calculated based on the inclusion records of the most recent W rounds at the k-th iteration. Its physical meaning is: in the most recent W iterations, what proportion of all newly included elite samples were generated by elite pooling operations? A higher value indicates that large-scale macro-level steps (elite pooling) are dominating the elite pool update, with coarse-grained samples continuously crowding out the inclusion opportunities of fine-grained samples; a lower value indicates that small-step micro-level search (small-step perturbation) is dominating the elite pool update.

[0089] , Where P refers to the average value of the elite mixed entry rate within a sliding window, L refers to the length of the sliding window, T refers to the elite mixed entry judgment value, and k and t represent the meanings explained above.

[0090] Obtain the average threshold of the sliding window for the pre-defined elite mixed entry rate in the database.

[0091] If the average sliding window value of the elite mixed entry rate is less than the average threshold of the elite mixed entry rate, then the coarse-grained elite interference execution strategy is recorded as continuing local refinement.

[0092] It should be noted that if the average sliding window value of the elite mixed entry rate is less than the average threshold of the sliding window value of the elite mixed entry rate, it means that the current search state of Phase 2 is normal and no intervention is required. At this time, the algorithm is converging towards the minimum value region as expected. If forced cleaning or changing the strategy is performed, it will destroy the current benign convergence trend. Therefore, the coarse-grained elite interference execution strategy is called elite interference cleaning adjustment.

[0093] If the average value of the sliding window for the elite mixed entry rate is greater than or equal to the average threshold of the sliding window for the elite mixed entry rate, then the coarse-grained elite interference execution strategy is recorded as performing elite interference cleaning and adjustment.

[0094] It should be noted that if the average sliding window value of the elite pool inclusion rate is greater than or equal to the average threshold of the elite pool inclusion rate, it indicates that the macroscopic large-step search in Stage 2 is dominating the update of the elite pool. This will cause coarse-grained samples to continuously crowd out fine-grained samples. Once the elite pool is dominated by coarse-grained samples, a vicious cycle of "self-reinforcement" is formed, and fine convergence will be continuously interrupted. When coarse-grained candidate solutions generated by elite pooling continue to be included in the pool, each round may have a solution that "leaps" into the elite pool and may become the new current optimal solution. This will cause the search center of the algorithm to be unable to stabilize within a small region. The current optimal solution will jump continuously in the parameter space (instead of gradually converging), causing the RMSE curve to exhibit sawtooth fluctuations and failing to smoothly decrease to a minimum value. Ultimately, the obtained parameter accuracy is lower than the normal convergence level, wasting the computational cost of model operation.

[0095] When making a judgment, intermittent detection is performed.

[0096] Perform elite interference cleaning and adjustment.

[0097] The system status is used to determine the number of intermittent detection rounds of samples.

[0098] The intermittent detection round value refers to the number of iteration rounds between two cleaning detections. The cleaning detection itself requires calculating statistics such as the average distance to the elite database and the proportion of sources from which data enters the database, incurring computational overhead. Detecting in every round would add unnecessary computational burden. In most cases, the state of the elite database changes slowly and gradually, making each round of detection unnecessary. Using intermittent detection can reduce computational costs while ensuring timely response.

[0099] The sample intermittent detection round value is obtained by analyzing the grain warehouse temperature field data acquired by the system at this time, which is beneficial for elite mixing in stage two.

[0100] Based on the intermittent detection round value, an intermittent detection strategy is executed to determine coarse-grained elite interference.

[0101] Assuming that elite interference detection begins in the 10th round after the start of Phase 2, the average value of the sliding window for the elite mixed entry rate is obtained by acquiring the sliding window length and elite mixed entry judgment value of the next 10 rounds. Based on this, it is determined whether elite interference cleaning and adjustment should be performed. Elite interference detection is not performed in the next 5 rounds, and the interval is rested. After the rest, the next round of elite interference detection and judgment is performed.

[0102] Obtain the timestamp of the samples added to the elite library.

[0103] The entry timestamp of an elite sample refers to the iteration round number at which the elite sample was added to the elite pool. It is a positive integer that records the time when the sample was "born". In each iteration, when a candidate solution is successfully added to the pool, the algorithm assigns the current iteration round number to the entry timestamp of that sample. Therefore, it is an existing variable that is automatically recorded when a sample is added to the pool.

[0104] Get the current iteration round number and the total number of iteration rounds in Phase 2 for the sample.

[0105] The current iteration count represents the number of iterations the algorithm has executed so far. It is a monotonically increasing integer variable, starting from 0 or 1 and incrementing by 1 after each iteration. This is the most fundamental state variable during the algorithm's operation and does not require external intervention. The total number of iterations in Phase 2 is the preset maximum number of iterations for Phase 2 (local refinement phase). This determines the maximum number of iterations Phase 2 will execute before automatically terminating, or before entering Phase 3. It is a hyperparameter preset during algorithm initialization and is determined before Phase 2 begins.

[0106] The difference in the number of iteration rounds is obtained based on the current iteration round number and the total number of iteration rounds in Phase 2.

[0107] The difference in the number of iterations is obtained by subtracting the total number of iterations in Phase 2 from the current iteration number.

[0108] Obtain the preset elimination ratio coefficient.

[0109] The elimination ratio coefficient is a ratio coefficient between 0 and 1, used to calculate the elimination threshold and determine which "old enough" samples should be eliminated. It is a hyperparameter preset during algorithm initialization.

[0110] The dynamic elimination value is obtained by using the difference in the number of iteration rounds, the total number of iteration rounds, and the elimination ratio coefficient.

[0111] , Where B refers to the dynamic elimination value, and b refers to the total number of iteration rounds. This refers to the difference in the number of iteration rounds. This refers to the elimination ratio coefficient.

[0112] The dynamic elimination value is a dynamically changing value that changes dynamically with different iteration rounds and different elite mixing results. This is intended to interfere with the removal of elites from the pool, thereby ensuring the purity of the elite pool and providing accurate data for subsequent judgments.

[0113] Elite interference cleaning and adjustment are performed based on dynamic elimination values.

[0114] Samples with elimination values ​​less than the dynamic elimination value in the current elite pool are eliminated. After adjustment, "outdated coarse-grained elites" are removed in batches, forcing the elite pool to consist only of recently added fine-grained samples. Since most of the samples added in the later stages of Phase Two come from small perturbations, the "purity" of the elite pool is guaranteed after cleaning. Even if subsequent elite mixing is done randomly, it can only be performed among fine-grained samples, and the difference vector naturally shrinks to a microscale.

[0115] S3, the third stage, namely behavior correction, uses focused perturbation plus elite cross-referencing to conduct a small-scale search, optimizing behavior constraints while maintaining accuracy, to ensure that the final parameters are not only "well-fitted" but also "physically reasonable".

[0116] Behavioral correction refers to: taking the optimal solution obtained from the second stage convergence as the initial starting point, performing focused perturbation operation and elite cross operation, and correcting the physical behavior constraints of the model parameters while keeping the root mean square error within the preset tolerance range, so that the final calibration parameters simultaneously meet the requirements of mathematical accuracy and physical rationality.

[0117] Phase three enters the "behavioral correction" phase, where the goal shifts from "finding the mathematical optimum" to "finding the physically reasonable mathematical optimum." At this point, the samples in the elite pool are very close to the region of the final solution, and the distance between samples is very small.

[0118] The elite crossover operation is as follows: select two elite samples from the elite sample library as parents, independently assign mixing coefficients to each dimension in the parameter space, and perform a weighted combination of the parent parameter values ​​in each dimension to generate new candidate parameter combinations.

[0119] While maintaining the sample accuracy of Phase 2, perform final parameter fitting, identify invalid crossover strategies, and adjust invalid crossover accordingly.

[0120] Obtain the current perturbation amplitude of the sample.

[0121] Here, the current perturbation amplitude represents the magnitude of the random displacement applied to the parameters by the focused perturbation operation in a certain iteration. It determines the "step size" of the small perturbation; the larger the current perturbation amplitude, the greater the change in parameter values ​​caused by the perturbation; the smaller the current perturbation amplitude, the more refined the perturbation. The current perturbation amplitude is a variable calculated in real time by the decay formula during algorithm execution, and its initial value and decay method are preset during algorithm initialization.

[0122] The typical acquisition method is exponential decay: , in, This refers to the current perturbation amplitude of the sample. This refers to the initial disturbance amplitude (a hyperparameter preset by the system). This refers to the attenuation rate (a system-preset hyperparameter). This refers to the current iteration round (which is an existing variable in the system).

[0123] The average distance between two samples in the elite database in the normalized parameter space is denoted as the sample average distance.

[0124] The average sample distance refers to the average Euclidean distance between all pairs of samples in the elite pool in the normalized parameter space. It reflects the dispersion of the samples in the current elite pool. The larger the average sample distance, the more dispersed the samples are in the elite pool; the smaller the average sample distance, the more concentrated (highly clustered) the samples are in the elite pool. The average sample distance is calculated in real time each round from the coordinate values ​​of all samples in the elite pool.

[0125] If the average distance between samples is greater than or equal to the current perturbation magnitude, the invalid crossover strategy is recorded as continuing behavior correction.

[0126] It should be noted that if the average distance between samples is greater than or equal to the current perturbation amplitude, it means that the distance between the two parents is large enough, and the displacement of the offspring generated by elite crossover relative to the parents can reach or exceed the displacement generated by the perturbation. In this case, elite crossover can still generate effective new directions, and its search contribution is no less than that of the perturbation operation. Therefore, the crossover operation is still valuable and can allow the algorithm to continue to perform behavior correction. Forcibly reducing the crossover probability will result in the loss of effective search directions. Therefore, invalid crossover strategies are recorded as continuing to perform behavior correction.

[0127] If the average distance between samples is less than the current perturbation amplitude, then the invalid crossover strategy is recorded as invalid crossover adjustment.

[0128] If the average distance between samples is less than the current perturbation magnitude, it means that the average distance between samples in the elite pool is less than the current perturbation step size, the distance between the two parents is very small, and the displacement of the offspring generated by elite crossover is less than the displacement generated by perturbation. This will result in the offspring generated by crossover being almost identical to the parents. The high cost and low benefit of crossover are mismatched; that is, a higher-cost operation produces a smaller search displacement, and computational resources are wasted. Crossover operations crowd out computational opportunities for perturbation operations, and rounds that should have been used for perturbation (the only effective search method at this time) are wasted, which will significantly reduce the fine convergence speed in the later stages of stage three, or even cause it to stagnate. To adjust invalid crossovers to change the current problem, invalid crossover strategies are denoted as invalid crossover adjustment.

[0129] The ratio of the average distance of the samples to the current perturbation amplitude is denoted as the elite cross-dynamic ratio.

[0130] , Where X refers to the elite cross-dynamic ratio. This refers to the average distance between samples. This refers to the current disturbance magnitude.

[0131] Obtain the minimum and maximum execution probabilities of elite crossover.

[0132] Among them, the minimum and maximum execution probabilities of elite crossover are preset hyperparameters, that is, the optimal boundary values ​​of elite crossover execution probabilities set by the system, which are obtained directly from the system when used.

[0133] Obtain the coefficients for determining invalid crossovers; The judgment coefficient refers to a preset threshold parameter used in Phase 3 to determine whether elite crossover is becoming ineffective. It is a dimensionless positive number used to compare with the ratio of the average distance of the elite pool to the current perturbation amplitude to determine whether the execution probability of elite crossover needs to be reduced.

[0134] Ineffective cross-adjustment is performed based on the data obtained above.

[0135] , in, This refers to the probability of elite cross-validation (i.e., the object to be adjusted for ineffective cross-validation). This refers to the minimum execution probability of elite crossover, which is a preset hyperparameter set by the system. This refers to the highest execution probability of elite crossover, which is a preset hyperparameter set by the system. sigmoid(x) is a sigmoid function that maps the input to the range of 0 to 1. This refers to the judgment coefficient for invalid crossovers.

[0136] The elite crossover execution probability is a core parameter controlling the algorithm's behavioral tendency in Phase 3 (behavior correction phase). It represents the probability that the algorithm will execute the elite crossover operation in each iteration. The elite crossover execution probability is set to a fixed constant during algorithm initialization.

[0137] at this time This is the adjusted elite crossover execution probability, which is a control parameter already present in Phase 3 (required in each iteration to determine which operation to execute). In traditional schemes, it is a fixed constant; this method only changes its assignment method—from "fixed preset" to "dynamically calculated by system state." When the elite pool is highly clustered and crossover becomes ineffective, the crossover probability is automatically reduced to the lower limit, freeing up computational resources from ineffective crossover operations and shifting them to the only effective focus perturbation operation at this point. The core task in the later stages of Phase 3 is fine-tuning, and focus perturbation is the most effective search method at this time. After reducing the ineffective crossover probability, focus perturbation gets more execution opportunities, accelerating the decrease in RMSE and Score, avoiding the misjudgment of "false convergence," and improving the accuracy of the final parameters.

[0138] After the above three stages are completed, the process of generating candidate parameters in each round can be formalized as follows: , In the formula, For the first Round candidate parameters represent new parameter combinations generated by the search strategy. The feasible region pruning operator means pruning the generated candidate parameters to a physically reasonable feasible region. Inside. For the first The candidate generation operator for the round represents the operator based on the current stage. Functions for adjusting search scale and strategy. For the first The elite sample library of the round represents a set that preserves historically optimal parameters and provides search directions. For the first The optimal parameters of the wheel represent the best combination of parameters for the current search center point. The random perturbation term represents a random variable used to avoid the search from getting trapped in a local optimum.

[0139] The advantage of a phased search strategy is that it first broadly searches for directions and then narrows down the search to specific areas. This strategy not only improves search efficiency but also ensures the robustness of the final results.

[0140] After phased elite archiving and calibration of the three-dimensional temperature field mechanism model of the grain warehouse, the optimal parameter combination was obtained. The combination was then subjected to multi-dimensional accuracy verification and physical consistency test to confirm the effectiveness of the simulation results.

[0141] By integrating the obtained grain warehouse temperature field observation data, analyzing the three-dimensional temperature field model, and determining the validity of the simulation results from multiple perspectives.

[0142] 1. Overall accuracy: After multiple rounds of iterative optimization, the model achieved a historical low RMSE of 0.898°C and a comprehensive score of 18.224. Table 1 presents the statistical results of the stratified accuracy.

[0143]

[0144] 2. Time series comparison: Figure 5 The time series comparison between simulated and observed temperatures for each layer is presented. It can be seen that the model can capture the temperature change trends of each layer well, especially the upper and lower layers, where the fitting effect is the best (RMSE = 0.816°C and 0.785°C, respectively). The upper layer, due to direct boundary forcing, has a relatively large error (RMSE = 1.074°C).

[0145] 3. Spatial distribution of temperature field: Figure 6 The diagram shows the planar distribution of the temperature field at times with larger and smaller errors, where actual observations are represented numerically and simulated values ​​are rendered in color. At times with larger errors (late September), the hot spot at the top is obvious, and the cold band on the sidewalls is clear; at times with smaller errors (late October), the temperature distribution tends to be uniform. Figure 7 The vertical profile shows the average temperature distribution along the length of the chamber. The simulation results agree well with the measured values, verifying the model's accurate characterization of the interlayer heat transfer process.

[0146] 4. Boundary-driven analysis: Figure 8 The study demonstrates the trends in outdoor temperature, indoor temperature, and overall average temperature within the grain pile. Outdoor temperature exhibits the largest fluctuations (range 10–30°C), followed by indoor temperature. The overall average temperature of the grain pile changes most gradually, reflecting the thermal inertia effect of the grain pile. The simulated average temperature effectively follows the observed average temperature changes, indicating that the model's response to boundary-driven conditions is reasonable.

[0147] 5. Optimize trajectory: Figure 9The optimization trajectory of parameter calibration is shown. Looking at the 452 iterations, the RMSE decreased rapidly in the early stage (first 100 iterations), from approximately 1.8°C to around 1.1°C; in the middle stage (100-350 iterations), it continued to decrease slowly and entered a plateau; in the later stage (after 350 iterations), it fluctuated slightly around 1.0°C, eventually reaching the historical low of RMSE = 0.898°C. The score continued to optimize in the later stage, indicating that the behavioral penalty term effectively guided the model to evolve in a direction more consistent with physical laws.

[0148] It's important to note that this step ensures the effectiveness of the calibration parameters, rather than resulting in "overfitting pseudooptimality." Simply reducing RMSE might only mean "memorizing" the training data, not "learning" the true physical laws. Through multi-dimensional validation, it was confirmed that the model performed well at spatial points (different layers) and time points (different time intervals) where calibration was not involved, proving that the parameters have generalization ability, rather than overfitting. This ensures that the calibration parameters have reliable predictive power in practical applications, rather than just performing well on training data.

[0149] To provide systematic evidence for the credibility of a model, a single indicator can only indicate that it "fits well," not that it is "physically reasonable." By cross-validating from multiple perspectives, such as stratified statistics, time series comparisons, spatial distribution comparisons, and boundary response analysis, comprehensive evidence is provided to support the model's credibility. Users have confidence in the model's prediction results and are willing to make decisions based on the model's output.

[0150] To support the final confirmation of "physical rationality," Phase Three introduces a behavioral penalty term to constrain physical laws. However, this penalty term is only a soft constraint and cannot guarantee that the final parameters are physically rational. Through physical consistency checks (such as whether the temperature gradient direction is correct, whether the temperature change rate of each layer is reasonable, and whether the thermal inertia conforms to the characteristics of the grain pile), the simulation results are rigorously examined to confirm the true effect of the behavioral correction. The final parameters not only meet mathematical accuracy standards but also pass independent physical tests, proving their thermodynamic reliability and completely avoiding the problem of "mathematical pseudo-optimality."

[0151] Identify the sources of model bias to guide future improvements.

[0152] A reproducibility standard for the calibration method was established. By defining clear verification indicators, a quantifiable evaluation standard was created for the application of this method in different grain warehouses and different grain varieties. After using this method, any implementer can use the same verification process to determine whether the calibration results are qualified. This method has universality and reproducibility across scenarios and varieties, rather than being a special case solution applicable only to a single grain warehouse.

[0153] In summary, this embodiment provides a calibration method for grain storage temperature field model parameters. It proposes a Bayesian-inspired, phased elite archiving calibration method, using a three-dimensional unsteady-state heat conduction equation as its core to construct a mechanistic model incorporating boundary heat transfer and stratified heat sources. Machine learning-based search is then used to efficiently optimize the parameters in a high-dimensional parameter space. This method effectively solves the problem of under-constrained parameters in traditional calibration methods under sparse data scenarios, possessing both physical interpretability and engineering usability, thus providing a scientific basis for grain storage safety management.

[0154] By judging and adjusting the coarse elite interference that occurs in the elite mixing of Phase 2, the search bias caused by the difference in dimensions is eliminated, synchronous convergence is achieved, invalid mixing is avoided, and computational efficiency is improved. When the elite library is highly clustered, the mixing operation is automatically skipped, and only small perturbations are performed. This allows invalid mixing to be accurately identified and skipped, reduces the number of model evaluations, and focuses computational resources on truly effective search directions, thus reducing the overall time consumption of Phase 2. The continuous interference of coarse-grained elites on the later fine search is eliminated, and the "average age" of the elite library is automatically reduced as the iteration progresses. The later elite library is dominated by fine-grained samples, and the direction generated by mixing always fits the current fine search scale. The RMSE curve declines smoothly and is no longer frequently interrupted by the early large step-size directions. Adaptive matching between the mixing amplitude and the search progress is achieved. Phase 2 smoothly transitions from "large step crossing" to "small step fine-tuning", eliminating the "scale conflict" between the mixing step size and the perturbation step size. The convergence curve is stable and smooth, without sawtooth oscillations.

[0155] By adjusting the invalid crossovers that occur in the elite crossover process of Stage 3, precise allocation of computational resources is achieved, eliminating invalid crossovers and reducing the proportion of invalid crossovers from a fixed probability. This allows computational resources to be precisely allocated to the focused perturbation operation (the only effective search method at this point), significantly improving the fine convergence speed in the later stages of Stage 3. Parent selection is aligned with the optimization objective, ensuring the correct correction direction. The offspring combinations generated by elite crossovers always point towards physically reasonable and acceptable accuracy directions, avoiding the "direction misalignment" problem that occurs when using the RMSE standard to select parents but still results in physically unreasonable regions. Perturbations are automatically amplified at penalty boundaries to ensure boundary reachability. The behavioral penalty term is only triggered occasionally when parameters violate physical constraints. When active penalty is detected, the exponential decay of the perturbation amplitude is automatically paused or even temporarily amplified. This makes candidate solutions longer near the physical constraint boundaries, capable of crossing the boundaries into physically reasonable regions, and shorter when far from the boundaries, ensuring normal fine convergence. This avoids the problem of "not being able to reach" the boundaries due to excessively small step sizes.

[0156] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods described in the various embodiments or some parts of the embodiments.

[0157] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for calibrating the parameters of a grain storage temperature field model, characterized in that, include: Collect and preprocess environmental and status data of grain warehouses to obtain temperature field data of grain warehouses, and construct a three-dimensional temperature field mechanism model of grain warehouses based on three-dimensional unsteady heat conduction equations; Implement a phased elite archiving rate determination method, including: Phase 1: Using random sampling and extreme perturbations, a large-scale search is conducted to escape local optima and establish an initial elite pool; Phase 2: Employ a hybrid elite solution with small perturbations to perform a medium-range search, followed by a fine-grained search around the elite solution to quickly converge to the region of minimum error and reduce RMSE; Phase 3: Employ focused perturbation combined with elite cross-validation to conduct a small-scale search and optimize behavioral constraints while maintaining accuracy; After phased elite archiving and calibration of the three-dimensional temperature field mechanism model of the grain warehouse, the optimal parameter combination was obtained. The combination was then subjected to multi-dimensional accuracy verification and physical consistency test to confirm the effectiveness of the simulation results.

2. The calibration method for the parameters of a grain storage temperature field model according to claim 1, characterized in that, The process of elite mixing includes: Two elite samples are selected from the elite sample library as parents, the difference vector between the two parents is calculated, and a new combination of candidate parameters is generated along the direction of the difference vector. The elite sample library obtained in Phase 1 is updated to generate new candidate parameter combinations. The average sliding window value of the elite mixed entry rate is analyzed to determine the coarse-grained elite interference execution strategy.

3. The calibration method for the parameters of a grain storage temperature field model according to claim 2, characterized in that, The process of analyzing the moving window average of the elite mixed inbound rate includes: Before the start of Phase Two, set a starting value for the elite mixed-input round; Perform the elite mixing operation in Phase 2. When the elite mixing entry round reaches the elite mixing entry round start value, calculate the sliding window length and the elite mixing entry judgment value. Based on the sliding window length and the elite mixed entry judgment value, the average sliding window value of the elite mixed entry rate is obtained.

4. The calibration method for the parameters of a grain storage temperature field model according to claim 2, characterized in that, The process of determining the coarse-grained elite interference execution strategy includes: Obtain the average threshold of the sliding window for the pre-defined elite mixed entry rate in the database; If the average sliding window value of the elite mixed entry rate is less than the average threshold of the elite mixed entry rate, then the coarse-grained elite interference execution strategy is recorded as continuing local refinement. If the average value of the sliding window of the elite mixed inbound rate is greater than or equal to the average threshold of the sliding window of the elite mixed inbound rate, then the coarse-grained elite interference execution strategy is recorded as performing elite interference cleaning and adjustment. When making a judgment, intermittent detection is performed; Perform elite interference cleaning and adjustment.

5. The calibration method for the parameters of a grain storage temperature field model according to claim 4, characterized in that, The process of performing intermittent testing includes: Based on the system status, determine the number of intermittent detection rounds for the samples; The sample intermittent detection round value is obtained by analyzing the grain warehouse temperature field data acquired by the system at this time; Based on the intermittent detection round value, a strategy for judging coarse-grained elite interference is implemented to perform intermittent detection.

6. The calibration method for the parameters of a grain storage temperature field model according to claim 4, characterized in that, The process of elite interference cleaning and adjustment includes: Obtain the timestamp of the samples added to the elite library; Obtain the current iteration round number and the total number of iteration rounds in Phase 2 for the sample; The difference in the number of iteration rounds is obtained based on the current iteration round number and the total number of iteration rounds in stage two. Obtain the preset elimination ratio coefficient; The dynamic elimination value is obtained by using the difference in the number of iteration rounds, the total number of iteration rounds, and the elimination ratio coefficient. Elite interference cleaning and adjustment are performed based on dynamic elimination values.

7. The calibration method for the parameters of a grain storage temperature field model according to claim 2, characterized in that, The process of elite cross-operation includes: Two elite samples are selected from the elite sample library as parent samples. Mixing coefficients are independently assigned to each dimension in the parameter space. The parent parameter values ​​of each dimension are weighted and combined dimension by dimension to generate new candidate parameter combinations. While maintaining the sample accuracy of Phase 2, perform final parameter fitting, identify invalid crossover strategies, and adjust invalid crossover accordingly.

8. The method for calibrating the parameters of a grain storage temperature field model according to claim 7, characterized in that, The process of determining an invalid crossover strategy includes: Obtain the current perturbation amplitude of the sample; The average distance between two samples in the elite pool in the normalized parameter space is denoted as the sample average distance. If the average distance of the samples is greater than or equal to the current perturbation amplitude, the invalid crossover strategy is marked as continuing the behavior correction; If the average distance between samples is less than the current perturbation amplitude, then the invalid crossover strategy is recorded as invalid crossover adjustment.

9. The method for calibrating the parameters of a grain storage temperature field model according to claim 8, characterized in that, The process of performing ineffective cross-adjustment includes: The ratio of the average distance of the samples to the current perturbation amplitude is calculated and denoted as the elite crossover dynamic ratio. Obtain the minimum and maximum execution probabilities of elite crossover; Obtain the coefficients for determining invalid crossovers; Ineffective cross-adjustment is performed based on the data obtained above.

Images