Method and system for optimizing self-adaptive operation strategy of refrigerating system

Abstract

The invention provides a self-adaptive operation strategy optimization method and system for a refrigeration system, and relates to the technical field of automatic control, and the method comprises the steps: obtaining temperature and humidity data and fluid boundary parameters, reconstructing a temperature and humidity distribution field through a physical constraint neural network, extracting enthalpy difference features, constructing a topological coupling network, and calculating and correcting a refrigeration load. And on the basis of energy consumption and enthalpy value uniformity, refrigeration resource distribution parameters are optimized, accurate regulation and control of the dynamic load area are achieved, the energy efficiency of the system is improved, the enthalpy value uniformity requirement of each area is met, and the environment change adapting capacity of the system is enhanced.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of automatic control, in particular to a self-adaptive operation strategy optimization method and system for a refrigeration system. BACKGROUND

[0002] A refrigeration system is an air conditioning system widely used in data centers, clean workshops, precision machining workshops and other places. This kind of system maintains the enthalpy value of the environment in an appropriate range by regulating the temperature and humidity in the space, ensuring the normal operation of the equipment and the stability of the process. At present, the operation strategy of this kind of refrigeration system is mainly based on the traditional temperature and humidity control method, which adjusts by setting fixed temperature and humidity thresholds. However, there are the following defects and deficiencies:

[0003] This kind of refrigeration system operation strategy lacks accurate perception ability of the enthalpy value distribution in the space, and only relies on limited temperature and humidity sampling points for control, which cannot accurately reflect the dynamic change characteristics of the enthalpy value of the whole space, resulting in unreasonable allocation of refrigeration resources and serious energy waste.

[0004] The prior art does not fully consider the thermal-hygroscopic coupling relationship between different areas in the refrigeration space, the mutual influence between areas is ignored, and the overall collaborative optimization control cannot be realized, resulting in the possibility of overcooling or overheating in local areas, affecting the stability and uniformity of the system.

[0005] The control strategy lacks self-adaptive adjustment mechanism and cannot flexibly adjust the control parameters according to the dynamic load change. The system response is lagging, and it is difficult to respond to sudden load changes, especially in the scene of high-density equipment distribution, the balance between refrigeration effect and energy consumption is difficult to optimize. SUMMARY

[0006] The embodiments of the present application provide a self-adaptive operation strategy optimization method and system for a refrigeration system, which can solve the problems in the prior art.

[0007] In a first aspect of the embodiments of the present application, a self-adaptive operation strategy optimization method for a refrigeration system is provided, comprising:

[0008] Obtaining temperature and humidity data and fluid boundary parameters in the space acted on by the refrigeration system;

[0009] Reconstructing a continuous temperature and humidity distribution field based on a physical constraint neural network, calculating an enthalpy distribution field and extracting a gradient change trajectory, dividing dynamic load areas according to gradient direction turning points and gradient periodic fluctuation points, and extracting load fluctuation intensity of each area as a local enthalpy difference feature;

[0010] The mutual correlation coefficient and phase difference of the time series of the load fluctuation intensity of adjacent areas are calculated to construct a topological coupling network; the enthalpy value decay rate of air supply arriving at each area is calculated according to the fluid boundary parameter; for each area, the local enthalpy difference feature, the coupling weight of the topological coupling network and the enthalpy value decay rate are calculated to obtain a modified refrigeration load through energy balance calculation;

[0011] The system energy consumption and the enthalpy value uniformity of each area are calculated based on the modified refrigeration load as an optimization target, the target weight coefficient is dynamically adjusted according to the load fluctuation intensity, and the refrigeration resource allocation parameter is solved;

[0012] Actual temperature and humidity data after the refrigeration resource allocation parameter is executed are collected, and the dynamic load area boundary is adjusted according to the deviation type of the actual enthalpy value distribution and the expected enthalpy value distribution.

[0013] The continuous temperature and humidity distribution field is obtained by field reconstruction of the temperature and humidity data based on the physical constraint neural network, and the steps of calculating the enthalpy value distribution field and extracting the gradient change trajectory include:

[0014] The physical constraint neural network includes a temperature prediction branch and a humidity prediction branch, the temperature prediction branch receives spatial position and time information as input, and the humidity prediction branch receives spatial position, time information and hidden layer features of the temperature prediction branch as input;

[0015] The loss function of the physical constraint neural network includes a data fitting term and a physical constraint term, the physical constraint term includes a thermodynamic conservation constraint and a time causality constraint; the time causality constraint calculates the enthalpy value according to the output of the temperature prediction branch and the humidity prediction branch, and constrains the time derivative sign of the enthalpy value to be consistent with the direction of energy inflow and outflow at this position; the weight proportion of the data fitting term and the physical constraint term is dynamically adjusted based on the density of temperature and humidity measurement points;

[0016] The continuous temperature and humidity distribution field is obtained by field reconstruction of the temperature and humidity data based on the physical constraint neural network, the enthalpy value distribution field is calculated, and the time gradient and spatial gradient of the enthalpy value distribution field are calculated, and the evolution process of the time gradient and spatial gradient in the time-space domain is taken as the gradient change trajectory.

[0017] The dynamic load area is divided according to the gradient direction turning point and the gradient periodic fluctuation point, and the load fluctuation intensity of each area is extracted as the local enthalpy difference feature, and the steps include:

[0018] The direction of the enthalpy value spatial gradient in the gradient change trajectory is tracked, and the spatial position where the gradient direction vector angle change rate exceeds a first set threshold is identified as the gradient direction turning point;

[0019] performing time domain analysis on the enthalpy value time gradient in the gradient change trajectory, and identifying a spatial position where a gradient amplitude change rate at adjacent time points exceeds a second set threshold and a change direction is reversed as a gradient periodic fluctuation point;

[0020] performing clustering analysis on the gradient direction turning points and the gradient periodic fluctuation points in spatial distribution to obtain a plurality of spatial clustering clusters, and taking a boundary of each spatial clustering cluster as a boundary of a dynamic load area;

[0021] calculating fluctuation characteristics of the enthalpy value time gradient in each dynamic load area as a load fluctuation intensity, and combining regional position information to form a local enthalpy difference feature.

[0022] The steps of calculating a cross-correlation coefficient and a phase difference of a load fluctuation intensity time sequence of adjacent areas to construct a topological coupling network include:

[0023] forming a load fluctuation intensity time sequence of each dynamic load area based on the load fluctuation intensity of each dynamic load area in the local enthalpy difference feature;

[0024] determining a spatial adjacent relationship between dynamic load areas, for two spatially adjacent dynamic load areas, calculating a cross-correlation function of their load fluctuation intensity time sequences, extracting a maximum value of the cross-correlation function as a cross-correlation coefficient, and extracting a time offset corresponding to the maximum value as a phase difference;

[0025] taking each dynamic load area as a network node, taking a cross-correlation coefficient between spatially adjacent dynamic load areas as a weight of a connection edge, and taking the phase difference as a direction attribute of the connection edge to construct a topological coupling network.

[0026] The steps of performing energy balance calculation on the local enthalpy difference feature, a coupling weight of the topological coupling network, and an enthalpy value decay rate to obtain a corrected refrigeration load include:

[0027] For each dynamic load area, extracting the load fluctuation intensity and regional boundary information of the area from the local enthalpy difference feature, and extracting the coupling weight and phase difference between the area and adjacent areas from the topological coupling network;

[0028] calculating a basic refrigeration load based on the load fluctuation intensity and regional boundary information of the area;

[0029] According to the coupling weight and phase difference between the dynamic load area and each adjacent area, calculating a heat and moisture transfer influence amount of the adjacent area on the area, which is equal to a product of the load fluctuation intensity and the coupling weight of the adjacent area and a time offset corrected cumulative value according to the phase difference;

[0030] According to the enthalpy attenuation rate of the dynamic load area, the effective refrigeration capacity of the supply air reaching the area is calculated, which is equal to the energy value corresponding to the difference between the nominal refrigeration capacity of the supply air and the enthalpy attenuation rate;

[0031] The base refrigeration load of the dynamic load area is added to the heat and moisture transfer impact quantity to obtain the total heat and moisture load of the area, and energy balance calculation is performed on the total heat and moisture load of the area and the effective refrigeration capacity to obtain the corrected refrigeration load of the dynamic load area.

[0032] Based on the corrected refrigeration load, system energy consumption and enthalpy uniformity of each area are calculated as optimization targets, target weight coefficients are dynamically adjusted according to load fluctuation intensity, and the steps of solving refrigeration resource allocation parameters include:

[0033] A state vector containing load fluctuation intensity and area coupling characteristics is constructed, and a target weight coefficient of each dynamic load area is generated based on a reinforcement learning strategy network according to the state vector, and the generation of the target weight coefficient is based on a collaborative adjustment strategy adopted by the topological coupling network;

[0034] Based on the coupling weight and phase difference in the topological coupling network, an influence matrix representing the degree of influence and time delay of heat and moisture transfer between areas is constructed;

[0035] Taking the supply air parameters of each dynamic load area as decision variables, taking the sum of the system energy consumption and the enthalpy uniformity index of each area weighted by the target weight coefficient as the optimization target, and taking the heat and moisture transfer constraint between areas defined by the influence matrix and the physical capacity constraint of the refrigeration equipment as the constraint condition, a multi-objective optimization function is constructed, and the refrigeration resource allocation parameters of each dynamic load area are obtained by solving the multi-objective optimization function;

[0036] The refrigeration resource allocation parameters are executed, the actual running effect is monitored, the parameters of the reinforcement learning strategy network are updated according to the monitoring results, and online learning and continuous optimization of the target weight coefficient adjustment strategy are realized.

[0037] The steps of generating the target weight coefficient based on the collaborative adjustment strategy adopted by the topological coupling network include:

[0038] Based on the distance relationship and connection mode between areas in the topological coupling network, a collaborative adjustment matrix is constructed combined with the enthalpy gradient change trajectory of each area, and the elements in the collaborative adjustment matrix represent the influence relationship of the adjustment of the target weight coefficient between areas;

[0039] According to the element values of the collaborative adjustment matrix, the areas are divided into groups, the areas with an influence relationship exceeding a preset influence threshold are determined as a collaborative area group, the master-slave relationship of the areas in the group is determined based on the size of the influence relationship, and an initial adjustment scheme of the target weight coefficient of the collaborative areas is generated;

[0040] According to the cooperative adjustment matrix, a change trend of the enthalpy gradient in the cooperative region group is calculated, when the change trend of the enthalpy gradient exceeds a preset range, a target weight coefficient of a corresponding region in the initial adjustment scheme is adjusted until the change trend of the enthalpy gradient meets a preset range requirement;

[0041] An adjustment scheme meeting the preset range requirement is taken as a final target weight coefficient, and is used for solving the multi-objective optimization function.

[0042] A second aspect of the embodiment of the present application provides a refrigeration system adaptive operation strategy optimization system, comprising:

[0043] A first unit is configured to acquire temperature and humidity data and fluid boundary parameters in a space acted on by a refrigeration system;

[0044] A second unit is configured to perform field reconstruction on the temperature and humidity data based on a physical constraint neural network to obtain a continuous temperature and humidity distribution field, calculate an enthalpy distribution field and extract a gradient change trajectory, divide dynamic load regions according to gradient direction turning points and gradient periodic fluctuation points, and extract load fluctuation intensities of the regions as local enthalpy difference features;

[0045] A third unit is configured to calculate mutual correlation coefficients and phase differences of time series of the load fluctuation intensities of adjacent regions to construct a topological coupling network, calculate enthalpy attenuation rates of air supply when reaching the regions according to the fluid boundary parameters, and perform energy balance calculation on the local enthalpy difference features, coupling weights of the topological coupling network, and the enthalpy attenuation rates to obtain a corrected refrigeration load for each region;

[0046] A fourth unit is configured to calculate system energy consumption and enthalpy uniformity of the regions as optimization objectives based on the corrected refrigeration load, dynamically adjust target weight coefficients according to the load fluctuation intensities, and solve to obtain refrigeration resource allocation parameters;

[0047] A fifth unit is configured to collect actual temperature and humidity data after the refrigeration resource allocation parameters are executed, and adjust boundaries of the dynamic load regions according to a deviation type of actual enthalpy distribution and expected enthalpy distribution.

[0048] A third aspect of the embodiment of the present application,

[0049] An electronic device is provided, comprising:

[0050] a processor;

[0051] a memory for storing processor-executable instructions;

[0052] The processor is configured to invoke the instructions stored in the memory to execute the method described above.

[0053] A fourth aspect of the embodiment of the present application,

[0054] A computer readable storage medium is provided, and computer program instructions are stored on the computer readable storage medium, and the computer program instructions are executed by a processor to implement the method.

[0055] The present application realizes accurate reconstruction of the temperature and humidity field by physically constraining the neural network, scientifically divides the dynamic load area based on the gradient change trajectory, accurately captures the complex heat flow dynamic characteristics, and greatly improves the operation state perception accuracy.

[0056] The topology coupling network constructed by the present application can represent the transmission relationship of load fluctuations between regions, and combined with enthalpy decay rate for energy balance calculation, a modified refrigeration load more consistent with the actual physical process is obtained, avoiding the load estimation deviation caused by ignoring the heat exchange between regions in the traditional method, and making the refrigeration resource distribution more reasonable.

[0057] The present application dynamically adjusts the optimization target weight according to the load fluctuation intensity, realizes adaptive balance between system energy consumption and enthalpy uniformity, and adjusts the dynamic load area boundary through the deviation feedback of actual enthalpy distribution and expected enthalpy distribution, forming a closed-loop adaptive control mechanism, improving the adaptability of the system to environmental changes, and significantly reducing the system energy consumption while ensuring the comfort of each region. BRIEF DESCRIPTION OF DRAWINGS

[0058] Figure 1 A flowchart of the refrigeration system adaptive operation strategy optimization method of the embodiment of the present application is shown in

[0059] Figure 2 A flowchart of the energy balance modified refrigeration load calculation based on local features and topology coupling is shown in DETAILED DESCRIPTION

[0060] To make the purpose, technical scheme and advantages of the embodiments of the present application clearer, the technical scheme in the embodiments of the present application will be described clearly and completely below in combination with the drawings of the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0061] The technical scheme of the present application will be described in detail below with specific embodiments. The following specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described in some embodiments.

[0062] Figure 1 A flowchart of the refrigeration system adaptive operation strategy optimization method of the embodiment of the present application is shown in Figure 1 As shown in the figure, the method comprises:

[0063] acquire temperature and humidity data and fluid boundary parameters in a space of a refrigeration system;

[0064] reconstruct a continuous temperature and humidity distribution field based on a physically constrained neural network, calculate an enthalpy distribution field and extract a gradient change trajectory, divide dynamic load areas according to gradient direction turning points and gradient periodic fluctuation points, and extract load fluctuation intensity of each area as a local enthalpy difference feature;

[0065] calculate a cross-correlation coefficient and a phase difference of a time series of load fluctuation intensity of adjacent areas to construct a topological coupling network, calculate an enthalpy decay rate of air supply when reaching each area according to the fluid boundary parameters, and perform energy balance calculation on the local enthalpy difference feature, coupling weight of the topological coupling network, and the enthalpy decay rate to obtain a corrected refrigeration load for each area.

[0066] calculate system energy consumption and enthalpy uniformity of each area based on the corrected refrigeration load as an optimization target, dynamically adjust a target weight coefficient according to the load fluctuation intensity, and solve to obtain refrigeration resource allocation parameters;

[0067] collect actual temperature and humidity data after the refrigeration resource allocation parameters are executed, and adjust the boundary of the dynamic load area according to the deviation type of the actual enthalpy distribution and the expected enthalpy distribution.

[0068] In an optional implementation, the step of reconstructing a continuous temperature and humidity distribution field based on a physically constrained neural network, calculating an enthalpy distribution field, and extracting a gradient change trajectory includes:

[0069] The physically constrained neural network includes a temperature prediction branch and a humidity prediction branch, the temperature prediction branch receives spatial position and time information as input, and the humidity prediction branch receives spatial position, time information, and hidden layer features of the temperature prediction branch as input.

[0070] The loss function of the physically constrained neural network includes a data fitting term and a physical constraint term, the physical constraint term includes a thermodynamic conservation constraint and a time causality constraint, the time causality constraint calculates an enthalpy value according to the output of the temperature prediction branch and the humidity prediction branch, and constrains the time derivative sign of the enthalpy value to be consistent with the direction of energy inflow and outflow at the position, and dynamically adjusts the weight proportion of the data fitting term and the physical constraint term based on the density of temperature and humidity measurement points.

[0071] reconstruct a continuous temperature and humidity distribution field based on the physically constrained neural network, calculate an enthalpy distribution field, and calculate a time gradient and a spatial gradient of the enthalpy distribution field, and take the evolution process of the time gradient and the spatial gradient in the time-space domain as a gradient change trajectory.

[0072] Exemplarily, the physical constraint neural network adopts a double-branch architecture to realize the joint reconstruction of temperature and humidity fields. The temperature prediction branch receives three-dimensional spatial coordinates and a timestamp as input, and the spatial coordinates are normalized to the interval of 0 to 1 in meters, and the timestamp is normalized to the interval of 0 to 1 of the current control period in seconds. The temperature prediction branch includes four fully connected layers, and the number of neurons of each layer is 128, 256, 256, and 128 in turn. The activation function adopts the hyperbolic tangent function, and the last layer outputs a single scalar representing the predicted temperature of the space-time position. The humidity prediction branch receives spatial coordinates and a timestamp, as well as a 256-dimensional hidden layer feature vector of the third fully connected layer of the temperature prediction branch. The feature vector is concatenated with the input of the humidity prediction branch and sent to the fully connected layer of the humidity prediction branch. The humidity prediction branch also includes four fully connected layers, and the number of neurons is configured as 384, 512, 256, and 128. The first layer has an input dimension of 387, including 3-dimensional spatial coordinates, 1-dimensional time, 256-dimensional temperature hidden layer features, and 127-dimensional self-hidden layer supplement. The activation function is consistent with the temperature branch, and the last layer outputs a single scalar representing the predicted humidity. The inter-layer connection weights of the temperature branch and the humidity branch are initialized using the Xavier method, and the bias term is initialized to 0.01.

[0073] The loss function of the physical constraint neural network is composed of a data fitting term and a physical constraint term. The data fitting term calculates the mean square error between the predicted value of the neural network and the actual measured point data. For temperature measurement points, the square of the difference between the predicted temperature and the measured temperature is calculated, and for humidity measurement points, the square of the difference between the predicted humidity and the measured humidity is calculated. The sum of the error squares of all measurement points is divided by the total number of measurement points to obtain the data fitting term value. The physical constraint term includes two parts: the thermodynamic conservation constraint and the time causality constraint. The thermodynamic conservation constraint requires that the temperature field and the humidity field satisfy the energy conservation equation. In the specific implementation, 1000 validation points are uniformly sampled in the field domain. For each validation point, the partial derivative of temperature with respect to time, the second-order partial derivative of temperature with respect to three spatial directions, the partial derivative of humidity with respect to time, and the second-order partial derivative of humidity with respect to three spatial directions are calculated. The automatic differentiation mechanism is used to calculate the derivatives of the neural network output. The temperature time derivative and the temperature spatial second-order derivative are linearly combined according to the thermal diffusion coefficient 0.0023 m² / s to calculate the residual square. The humidity time derivative and the humidity spatial second-order derivative are linearly combined according to the humidity diffusion coefficient 0.0019 m² / s to calculate the residual square. The sum of the thermodynamic conservation residual squares of all validation points is divided by the number of validation points to obtain the thermodynamic conservation constraint term value.

[0074] The time causality constraint requires the time derivative of enthalpy to be positive if the energy flux is positive, and negative if the energy flux is negative. For each validation point, the enthalpy is calculated from the outputs of the temperature prediction branch and the humidity prediction branch. The enthalpy calculation includes two terms, the first term is the temperature value multiplied by the air constant pressure specific heat 1.005 kJ / (kg·℃), the second term is the humidity value multiplied by the water vapor latent heat of vaporization 2501 kJ / kg, the sum of the two terms is the enthalpy of the point. The first-order partial derivative of enthalpy with respect to time is obtained, and the first-order partial derivatives of temperature with respect to three spatial directions are obtained to obtain the temperature gradient vector. The first-order partial derivatives of humidity with respect to three spatial directions are obtained to obtain the humidity gradient vector. The energy flux density vector is calculated from the temperature gradient vector and the humidity gradient vector. Each component of the energy flux density vector is calculated in two parts, the first part is the product of the negative thermal diffusion coefficient and the corresponding component of the temperature gradient, the second part is the product of the negative humidity diffusion coefficient and the corresponding component of the humidity gradient multiplied by the latent heat of vaporization, and the first part minus the second part is equal to the component of the energy flux density vector. The divergence of the energy flux density vector is calculated, and the divergence value is the sum of the partial derivatives of the three components of the energy flux density vector with respect to the corresponding spatial direction. The time causality constraint requires the product of the time derivative of enthalpy and the divergence of the energy flux density to be greater than or equal to zero. When the product is less than zero, a penalty term is generated, and the value of the penalty term is equal to the absolute value of the product squared. The time causality constraint term value is obtained by dividing the sum of the time causality penalty terms of all validation points by the number of validation points.

[0075] The measurement point density is defined as the number of measurement points divided by the field volume, and the field volume is in cubic meters. When the measurement point density is less than 0.1 / m³, the data fitting term weight is set to 0.3, the thermodynamic conservation constraint weight in the physical constraint term is set to 0.5, and the time causality constraint weight is set to 0.2. When the measurement point density is between 0.1 / m³ and 0.5 / m³, the data fitting term weight linearly increases to 0.6, the thermodynamic conservation constraint weight in the physical constraint term linearly decreases to 0.3, and the time causality constraint weight remains 0.1. When the measurement point density exceeds 0.5 / m³, the data fitting term weight is fixed at 0.8, the thermodynamic conservation constraint weight in the physical constraint term is fixed at 0.15, and the time causality constraint weight is fixed at 0.05. The total value of the loss function is calculated by adding three parts, the first part is the data fitting term multiplied by its weight, the second part is the thermodynamic conservation constraint term multiplied by its weight, and the third part is the time causality constraint term multiplied by its weight.

[0076] The neural network training adopts the Adam optimizer, and the initial learning rate is set to 0.001, which is attenuated by 10% every 500 iterations. The batch size is set to 256, including randomly sampled measurement point data and validation point data. The training process lasts for 5000 iterations, and the loss function value is evaluated on the validation set every 100 iterations. When the validation set loss does not decrease for three consecutive evaluations, the early stopping mechanism is triggered to terminate the training. After training, the neural network model saves the weight parameters and network structure configuration for subsequent field reconstruction reasoning.

[0077] In the field reconstruction reasoning stage, a three-dimensional space grid is established in the refrigeration system action space, divided into 50 nodes along the length direction, 30 nodes along the width direction, and 20 nodes along the height direction, totaling 30000 space grid points. For a given time, the coordinates of the 30000 space grid points and the time stamp of that time are input into the trained physical constraint neural network, and the temperature prediction branch and the humidity prediction branch output 30000 temperature prediction values and 30000 humidity prediction values respectively, forming a continuous temperature and humidity distribution field. According to the temperature distribution field and the humidity distribution field, the enthalpy distribution field is calculated. The enthalpy calculation of each grid point includes two additions, the first is the temperature multiplied by the air constant pressure specific heat, and the second is the humidity multiplied by the water vapor latent heat of vaporization. The enthalpy distribution field is stored in a three-dimensional array, with dimensions of 50x30x20.

[0078] The time gradient and spatial gradient of the enthalpy distribution field are calculated. The time gradient calculation uses the central difference method. For the current time, the enthalpy distribution field of the previous time and the next time is obtained. The enthalpy time gradient of each grid point at the current time is calculated as the difference between the enthalpy of the next time and the enthalpy of the previous time, and the difference is divided by the time interval. The time interval is 1 / 10 of the control period, typically 6s. The spatial gradient calculation also uses the central difference method. For each grid point, the enthalpy spatial gradient along the length direction is calculated as the difference between the enthalpy of the right adjacent point and the enthalpy of the left adjacent point, and the difference is divided by the length direction grid spacing. The grid spacing is obtained by dividing the length of the refrigeration system action space by 49. The enthalpy spatial gradient along the width direction and height direction is calculated in a similar manner, respectively divided by the grid spacing in the corresponding direction. Each grid point has a time gradient scalar and three spatial gradient components, forming a four-dimensional gradient vector.

[0079] The gradient change trajectory records the evolution process of the enthalpy time gradient and the spatial gradient in the time-space domain. For a control period, the enthalpy distribution field snapshot contains 10 time points, each time point has 30000 grid points, and each grid point has a four-dimensional gradient vector, which totals 300000 gradient vectors to form the gradient change trajectory of the control period. The gradient change trajectory is stored in a time sequence data structure, supporting multi-dimensional queries by time index, spatial location index, and gradient component index. For scenarios that require extraction of gradient evolution in a specific spatial region, the spatial coordinate range can be specified to extract the gradient vector sequence of all time points in the range. For scenarios that require analysis of gradient mutation, a gradient change rate threshold can be set, and when the difference between the gradient vectors of the same grid point at adjacent time points exceeds the threshold, it is marked as a mutation point, and the time-space location and gradient value of the mutation point are recorded to the mutation event list.

[0080] The present application realizes high-precision enthalpy distribution field construction. Through the collaborative design of temperature prediction branch and humidity prediction branch, the physical principles and data fitting are fused to ensure that the field reconstruction results meet the data consistency and physical rationality. In particular, the introduction of time causality constraint ensures that the enthalpy change conforms to the energy flow law, solving the limitations of traditional interpolation methods in physical consistency. Through dynamic weight adjustment based on the density of measuring points, the method can maintain good reconstruction accuracy in data sparse areas.

[0081] In an optional embodiment, the step of dividing the dynamic load area according to the gradient direction turning point and the gradient periodic fluctuation point, and extracting the load fluctuation intensity of each area as the local enthalpy difference feature includes:

[0082] Direction tracking is performed on the enthalpy spatial gradient in the gradient change trajectory, and a spatial position where the gradient direction vector angle change rate exceeds a first set threshold is identified as a gradient direction turning point;

[0083] Time domain analysis is performed on the enthalpy time gradient in the gradient change trajectory, and a spatial position where the gradient amplitude change rate at adjacent time points exceeds a second set threshold and the change direction is reversed is identified as a gradient periodic fluctuation point;

[0084] The gradient direction turning point and the gradient periodic fluctuation point are subjected to cluster analysis in spatial distribution to obtain a plurality of spatial cluster clusters, and the boundaries of each spatial cluster cluster are taken as the boundaries of the dynamic load area;

[0085] The fluctuation characteristics of the enthalpy time gradient in each dynamic load area are calculated as the load fluctuation intensity, and the local enthalpy difference feature is formed in combination with the regional position information.

[0086] Exemplarily, the gradient change trajectory contains evolution data of the time gradient and the space gradient of the enthalpy value in the entire control period, and the direction tracking is performed on the space gradient to identify the gradient direction turning point. The space gradient of the enthalpy value has three components at each grid point, which correspond to the gradient values of the length, width and height of the three spatial directions respectively. The method for calculating the gradient direction vector is to form a three-dimensional vector with the three components, and the direction angle of the vector is determined by the ratio of the three components. For the space gradient vectors of a grid point at two adjacent time instants, the included angle between the two vectors is calculated. The included angle is calculated by using the ratio of the dot product of the two gradient vectors to the product of the two vector modules, that is, the sum of the products of the corresponding components of the two gradient vectors is divided by the product of the two vector modules to obtain the cosine value of the included angle, and the radian value of the included angle is obtained by the inverse cosine function. The gradient direction vector angle change rate is defined as the included angle of the gradient vectors at adjacent time instants divided by the time interval, and the time interval is the time difference between adjacent snapshots in the control period. When the gradient direction vector angle change rate of a grid point exceeds a first set threshold, the grid point is marked as a gradient direction turning point. The first set threshold is determined according to the spatial scale of the space acted on by the refrigeration system and the load change rate, and the typical value range is 0.5 to 2.0 rad / s, and the default value is set to 1.0 rad / s. The distribution of the gradient direction turning points in space is recorded as a coordinate list, and the list elements contain three-dimensional space coordinates and the time when the turning is detected.

[0087] The time domain analysis of the enthalpy time gradient is used to identify the gradient periodic fluctuation point. Each grid point has a plurality of time gradient values of the enthalpy at time instants in the control period, forming a time sequence. The time sequence is analyzed point by point, and the change rate of the enthalpy time gradient between adjacent time points is calculated. The change rate is defined as the difference between the gradient value at the next time instant and the gradient value at the previous time instant divided by the time interval. The method for judging whether the change direction is reversed is to check the signs of the gradient change rates of the three consecutive time points. When the first change rate is positive and the second change rate is negative, or the first change rate is negative and the second change rate is positive, it is considered that the change direction is reversed. At the same time, it is required that the absolute value of the change rate exceeds a second set threshold, which is determined according to the fluctuation amplitude range of the enthalpy time gradient, and the typical value range is 0.01 to 0.1 kJ / (kg·s²), and the default value is set to 0.05 kJ / (kg·s²). The grid point that meets the two conditions of the change rate exceeding the threshold and the change direction being reversed is marked as a gradient periodic fluctuation point. The gradient periodic fluctuation points are also recorded as a coordinate list, and the list elements contain three-dimensional space coordinates and the time when the fluctuation is detected.

[0088] The gradient direction turning points and the gradient periodic fluctuation points are merged into a unified feature point set, and each element of the feature point set is a three-dimensional space coordinate. The clustering analysis adopts a density-based spatial clustering algorithm, which realizes automatic clustering by defining a neighborhood radius and a minimum point threshold. The neighborhood radius is defined as the maximum spatial distance of a feature point within which other feature points are considered neighbors. The typical value range is 0.5 to 2.0 m, and the default value is set to 1.0 m. The minimum point threshold is defined as the minimum number of feature points required to form a cluster. The typical value range is 5 to 20, and the default value is set to 10. The clustering algorithm traverses all feature points. For each feature point, the number of other feature points within its neighborhood radius range is calculated. When the number of feature points within the neighborhood reaches or exceeds the minimum point threshold, the feature point is marked as a core point. From any core point, all points within its neighborhood are included in the same cluster, and the neighborhood of the newly included points is recursively checked. Points with reachable density are continuously added to the cluster. Feature points not belonging to any cluster are marked as noise points. After clustering, multiple spatial clusters are obtained, and each cluster corresponds to a dynamic load area.

[0089] The boundary of the dynamic load area is determined by the peripheral points of the cluster. For each cluster, the minimum and maximum values of all feature points in the three spatial directions are calculated to form a cuboid bounding box that encloses the cluster. The six faces of the bounding box are determined by the length minimum value, length maximum value, width minimum value, width maximum value, height minimum value, and height maximum value, respectively. Each face of the bounding box is extended outward by a boundary buffer distance, which is set to half of the neighborhood radius to ensure that the bounding box completely covers the influence range of load changes in the region. The extended bounding box is used as the final boundary of the dynamic load area. When the bounding boxes of multiple clusters overlap, the overlapping region is assigned to the cluster with more feature points. The number of feature points of each cluster in the overlapping region is calculated, and the overlapping region is divided into the cluster with more feature points. The number of dynamic load areas is equal to the number of clusters, and each area has a unique area number that starts from 1 and increments. The region boundary is stored in the form of six boundary values, including the length minimum value, length maximum value, width minimum value, width maximum value, height minimum value, and height maximum value.

[0090] The load fluctuation intensity of each dynamic load region is obtained by calculating the fluctuation characteristics of the enthalpy time gradient in the region. For a certain dynamic load region, the enthalpy time gradient values of all grid points within the boundary range of the region in the entire control period are extracted. The standard deviation of these gradient values is calculated as the basic indicator of the fluctuation characteristics. The standard deviation is calculated as the square sum of the difference between all gradient values and the average gradient value divided by the total number of gradient values, and then the square root is taken. The standard deviation reflects the dispersion degree of the enthalpy time gradient, and the larger the standard deviation, the more intense the load fluctuation. In addition to the standard deviation, the peak-to-peak value of the gradient value is also calculated, which is defined as the maximum gradient value minus the minimum gradient value. The peak-to-peak value reflects the range of extreme values of the enthalpy time gradient, and the larger the peak-to-peak value, the greater the load fluctuation amplitude. The standard deviation and the peak-to-peak value are combined by weighting to obtain the load fluctuation intensity, and the weighting coefficients are 0.6 and 0.4 respectively. The load fluctuation intensity is equal to the standard deviation multiplied by 0.6 plus the peak-to-peak value multiplied by 0.4. The unit of load fluctuation intensity is kJ / (kg·s), and the numerical range is determined according to the load characteristics of the space acted on by the refrigeration system. The typical range is 0.05 to 0.5 kJ / (kg·s).

[0091] The local enthalpy difference characteristics are formed by combining the load fluctuation intensity and the region position information. The region position information includes the geometric center coordinates of the region boundary box, and the three components of the geometric center coordinates are the average of the minimum and maximum values of the length direction, the average of the minimum and maximum values of the width direction, and the average of the minimum and maximum values of the height direction. The region position information also includes the region volume, which is equal to the product of the length direction boundary difference, the width direction boundary difference and the height direction boundary difference. The local enthalpy difference characteristics are stored in a structured data form, and each dynamic load region corresponds to a record. The record fields include region number, geometric center length coordinate, geometric center width coordinate, geometric center height coordinate, region volume, load fluctuation intensity, length direction minimum boundary, length direction maximum boundary, width direction minimum boundary, width direction maximum boundary, height direction minimum boundary, and height direction maximum boundary. The local enthalpy difference characteristic set contains the records of all dynamic load regions, which are used for subsequent control strategy optimization.

[0092] The application innovatively expands the enthalpy difference principle from a test evaluation tool to a basic theory of refrigeration system operation control, realizes adaptive optimization of the refrigeration system operation strategy through accurate perception of the enthalpy value distribution field in the space, dynamic extraction of the regional enthalpy difference characteristics, and enthalpy difference-based energy balance calculation, and changes the enthalpy difference principle from a static test evaluation tool to a dynamic operation control method. The application breaks through the limitation of traditional fixed region division, can adaptively capture the dynamic change characteristics of the actual load distribution by identifying the direction turning points and periodic fluctuation points of the enthalpy value gradient, and uniformly considers the space gradient and the time gradient to realize automatic identification and boundary division of the load region. The dynamic division method can timely respond to the load distribution change and avoid the problem of uneven load in the region caused by fixed region division.

[0093] In an alternative embodiment, the step of calculating the cross-correlation coefficient and the phase difference of the load fluctuation intensity time series of adjacent regions to construct a topological coupling network comprises:

[0094] Based on the load fluctuation intensity of each dynamic load region in the local enthalpy difference characteristics, a load fluctuation intensity time series of each dynamic load region is formed;

[0095] The spatial adjacent relationship between the dynamic load regions is determined, for two spatially adjacent dynamic load regions, the cross-correlation function of the load fluctuation intensity time series thereof is calculated, the maximum value of the cross-correlation function is extracted as the cross-correlation coefficient, and the time offset corresponding to the maximum value of the cross-correlation function is extracted as the phase difference;

[0096] Each dynamic load region is taken as a network node, the cross-correlation coefficient between the spatially adjacent dynamic load regions is taken as the weight of the connecting edge, and the phase difference is taken as the direction attribute of the connecting edge to construct a topological coupling network.

[0097] For example, the local enthalpy difference characteristic of each dynamic load area contains a load fluctuation intensity value, which is the load fluctuation characteristic of the area at a certain time. Forming a load fluctuation intensity time series requires repeatedly calculating the load fluctuation intensity of each dynamic load area over multiple consecutive control periods. The time interval of the control period is set to 60s, and the continuous collection time window is set to 30 control periods, corresponding to an observation time length of 1800s. For a certain dynamic load area, at the end of each control period, the enthalpy time gradient values of all grid points within the boundary range of the area are extracted, the standard deviation and peak-to-peak value are calculated, and the load fluctuation intensity of the control period is obtained by weighting the standard deviation multiplied by 0.6 plus the peak-to-peak value multiplied by 0.4. The load fluctuation intensities of the 30 control periods are arranged in chronological order to form a time series with a length of 30. Each element in the time series corresponds to a load fluctuation intensity value at a certain time, and the element index ranges from 0 to 29, with index 0 corresponding to the first control period and index 29 corresponding to the thirtieth control period. The load fluctuation intensity time series is stored in the form of a one-dimensional array, with array elements being floating-point numbers, and the value range being determined according to the load characteristics of the area. The typical range is 0.05 to 0.5 kJ / (kg·s). Each dynamic load area generates a load fluctuation intensity time series, and if there are 6 dynamic load areas, 6 time series are generated, each with a length of 30.

[0098] The spatial adjacency relationship between dynamic load areas is determined by the geometric position of the area boundary. For any two dynamic load areas, the six boundary values of their bounding boxes are obtained, including the minimum length, maximum length, minimum width, maximum width, minimum height, and maximum height. The method of determining whether two areas are spatially adjacent is to check the shortest distance between the two bounding boxes. The coordinates of the eight vertices of the first area bounding box are calculated, each vertex coordinate being composed of the boundary values of length, width, and height, and there are eight combinations corresponding to eight vertices. Similarly, the coordinates of the eight vertices of the second area bounding box are calculated. Traverse each vertex of the first area and each face of the second area to calculate the shortest distance from the vertex to the face. Similarly, traverse each vertex of the second area and each face of the first area to calculate the shortest distance from the vertex to the face. The minimum value of all vertex-to-face distances is defined as the shortest distance between the two bounding boxes. When the shortest distance is less than or equal to the adjacency threshold, the two areas are considered spatially adjacent. The adjacency threshold is determined according to the spatial scale of the refrigeration system space and the area granularity, and the typical value range is 0.5 to 2.0 m, with a default value of 1.5 m. The spatial adjacency relationship is stored in the form of an adjacency matrix, with the row and column indices of the matrix corresponding to the dynamic load area numbers, and the matrix elements being Boolean values. When two areas are spatially adjacent, the corresponding element is true, otherwise it is false. The adjacency matrix is a symmetric matrix, and the diagonal elements are false, indicating that the area itself is not adjacent to itself.

[0099] For two dynamically loaded regions that are spatially adjacent, the cross-correlation function of their load fluctuation intensity time series is calculated. The cross-correlation function describes the similarity of two time series at different time offsets. Let the time series of the first region be sequence A, the time series of the second region be sequence B, and the time offset be parameter tau, tau ranges from negative sequence length plus one to positive sequence length minus one. For a given time offset tau, the cross-correlation function value is calculated by multiplying the elements of sequence A and the shifted sequence B point by point and summing them up. When tau is positive, sequence B is shifted backward by tau time steps, and the product of the first several elements of sequence A and the corresponding elements of the shifted sequence B is calculated. When tau is negative, sequence B is shifted forward by absolute value tau time steps, and the product of the last several elements of sequence A and the corresponding elements of the shifted sequence B is calculated. The effective calculation range is the overlapping part of the two sequences, and the length of the overlapping part is equal to the sequence length minus the absolute value of the offset. The point-by-point product of the overlapping part is summed up to obtain the cross-correlation original value at this offset. In order to eliminate the influence of sequence length and amplitude on the cross-correlation function, the original value is normalized. The normalization method is to divide the original value by the product of the standard deviations of the two sequences in the overlapping part and then divide by the overlapping length. The normalized cross-correlation function value ranges from negative one to positive one, positive values indicate positive correlation, negative values indicate negative correlation, and absolute values indicate correlation strength. By traversing all time offsets tau, the cross-correlation function value corresponding to each tau is calculated to form a cross-correlation function curve. The cross-correlation function curve is stored in the form of a one-dimensional array, and the array index corresponds to the time offset, and the array element is the normalized cross-correlation function value.

[0100] The maximum value of the cross-correlation function is taken as the cross-correlation coefficient, which reflects the strongest correlation degree of the load fluctuation intensity time series of the two regions. By traversing all elements of the cross-correlation function curve, the element with the largest absolute value is found, and the value of this element is the cross-correlation coefficient. The cross-correlation coefficient ranges from negative one to positive one, close to positive one indicates that the load fluctuations of the two regions are highly synchronized, close to negative one indicates that the load fluctuations of the two regions are highly reversed, and close to zero indicates that there is no obvious correlation between the load fluctuations of the two regions. The time offset corresponding to the maximum value of the cross-correlation function is taken as the phase difference, which reflects the time sequence of the load fluctuations of the two regions. The array index corresponding to the cross-correlation coefficient is extracted, and the index value is converted to the time offset, and the unit of the time offset is the number of control periods. A positive phase difference indicates that the load fluctuation of the second region lags behind the first region, and a negative phase difference indicates that the load fluctuation of the second region leads the first region. The absolute value of the phase difference reflects the degree of time delay, and the actual time delay is obtained by multiplying the phase difference by the control period time interval, with the unit of seconds. The cross-correlation coefficient and the phase difference are stored in the form of key-value pairs, with the key being the combination of the two region numbers and the value being a structure containing the cross-correlation coefficient and the phase difference.

[0101] The topology coupling network is represented by a data structure of graph, which contains a set of nodes and a set of edges. Each dynamic load region is a node of the network, and the number of nodes is equal to the number of dynamic load regions. Each node has a unique node identifier, which is consistent with the region number. The node attributes include the geometric center coordinates of the region, the volume of the region, and the boundary range of the region, which are extracted from the local enthalpy difference characteristics. The connection edges are established between the spatially adjacent dynamic load regions. The existence of the edge is determined by the adjacency matrix. For the elements with true values in the adjacency matrix, an edge is created between the corresponding two nodes. The weight of the edge is the cross-correlation coefficient, and the weight value ranges from -1 to 1. The greater the absolute value of the weight, the higher the coupling strength between the two regions. The direction attribute of the edge is the phase difference, which is an integer with the unit of control period number and the value range from negative sequence length plus one to positive sequence length minus one. The data structure of the edge contains four fields: starting node identifier, ending node identifier, weight, and phase difference. Since the cross-correlation relationship has symmetry, the weight of the edge from node i to node j is the same as that from node j to node i, but the sign of the phase difference is opposite. When the topology coupling network is represented by an undirected weighted graph, the edge is recorded only once, and the sign of the phase difference is agreed to be the phase difference when pointing from the node with smaller number to the node with larger number. The topology coupling network is stored in the form of adjacency list, and each node maintains a neighbor list. The list elements contain neighbor node identifier, edge weight, and phase difference. The global attributes of the network include the total number of nodes, the total number of edges, the average weight, and the standard deviation of the weight, which are used for network topology analysis and control strategy optimization.

[0102] The implementation method for calculating the enthalpy decay rate of air supply reaching each area according to the fluid boundary parameters is as follows: the fluid boundary parameters refer to a set of key parameters for describing the physical characteristics of the air supply system, including the temperature, humidity, air speed of the air supply outlet, the length of the air supply pipeline, the heat conduction coefficient of the pipeline material, the cross-sectional size of the pipeline, the number of elbows, the branch structure and the like. These parameters directly affect the energy loss of air in the transmission process. When calculating the enthalpy decay rate, the path from the air supply outlet to each dynamic load area is decomposed into three basic units of straight pipe sections, elbows and branches. For each straight pipe section, the heat loss coefficient a is calculated according to the pipe length L (m), the pipe diameter D (m), the wall heat conduction coefficient K (W / m·K) and the pipe wall thickness δ (m); for each elbow, the local resistance coefficient β is calculated according to the bending angle and the cross-section change; for each branch, the shunt influence coefficient γ is calculated according to the flow distribution ratio. These coefficients are combined with the air supply temperature difference ΔT (the difference between the air supply temperature and the ambient temperature) and the humidity difference ΔW to obtain the total energy loss by using the cumulative calculation method, and the energy loss is divided by the initial enthalpy to obtain the enthalpy decay rate. For example, a certain path contains a 15-meter straight pipe (a = 0.02 / m), two 90-degree elbows (β = 0.15 / each) and one three-way branch (γ = 0.1), the air supply temperature is 15°C, and the ambient temperature is 25°C. The calculated enthalpy decay rate is 15 x 0.02 x 10 + 2 x 0.15 x 10 + 0.1 x 10 = 5.0%. In this way, the complex heat fluid transmission process is simplified into an engineering calculation model, which takes into account accuracy and practicality.

[0103] The present application reveals the time sequence correlation and influence propagation characteristics of the load changes between regions by calculating the cross-correlation of the time series of the load fluctuation intensity of adjacent regions. The cross-correlation coefficient is introduced as the connection edge weight and the phase difference is introduced as the direction attribute, so that the network can represent the spatial adjacency relationship and the time sequence influence relationship at the same time. This representation method breaks through the limitation of only considering spatial adjacency in the traditional method and can capture the dynamic coupling characteristics between regions.

[0104] In an optional implementation, the step of performing energy balance calculation on the local enthalpy difference characteristics, the coupling weight of the topological coupling network and the enthalpy decay rate to obtain the modified refrigeration load comprises:

[0105] For each dynamic load area, the load fluctuation intensity and the region boundary information of the area are extracted from the local enthalpy difference characteristics, and the coupling weight and the phase difference between the area and the adjacent areas are extracted from the topological coupling network;

[0106] The basic refrigeration load is calculated based on the load fluctuation intensity and the region boundary information of the region;

[0107] According to the coupling weight and the phase difference between the dynamic load area and each adjacent area, the heat and moisture transfer influence of the adjacent area on the area is calculated, which is equal to the product of the load fluctuation intensity of the adjacent area and the coupling weight and the time offset correction value according to the phase difference;

[0108] According to the enthalpy decay rate of the dynamic load area, the effective cooling capacity when the supply air reaches the area is calculated, which is equal to the energy value corresponding to the difference between the nominal cooling capacity of the supply air and the enthalpy decay rate;

[0109] The total heat and moisture load of the dynamic load area is obtained by adding the basic cooling load of the dynamic load area and the heat and moisture transfer influence, and the energy balance calculation is performed between the total heat and moisture load of the area and the effective cooling capacity, to obtain the modified cooling load of the dynamic load area.

[0110] In combination Figure 2 Based on the energy balance modified cooling load calculation flowchart of local feature and topological coupling, for example, when the modified cooling load of each dynamic load area is calculated, the relevant parameters need to be extracted from the local enthalpy difference feature and topological coupling network. The local enthalpy difference feature is stored in a structured table form, and each row corresponds to a dynamic load area. For a dynamic load area, the load fluctuation intensity field and the area boundary information field of the area are read from the table. The load fluctuation intensity is a floating point number, with the unit of kJ / (kg·s), reflecting the fluctuation degree of the enthalpy time gradient of the area. The area boundary information includes six boundary values, which are the minimum boundary in the length direction, the maximum boundary in the length direction, the minimum boundary in the width direction, the maximum boundary in the width direction, the minimum boundary in the height direction, and the maximum boundary in the height direction, with the unit of meters. The topological coupling network is stored in the form of an adjacency table, and for a node corresponding to a dynamic load area, the node identifier, coupling weight and phase difference of all adjacent areas are extracted from its neighbor list. The coupling weight is a floating point number, with the value range of negative one to positive one, reflecting the correlation intensity of the load fluctuation of the area and the adjacent area. The phase difference is an integer, with the unit of control period number, reflecting the time offset amount of the load fluctuation of the adjacent area relative to the area. The extraction process is completed through the region number index, and the region number is used as the table row index and network node identifier to ensure consistent access of data.

[0111] The calculation of the basic cooling load is based on the load fluctuation intensity of the region and the region boundary information. The region boundary information is used to calculate the region volume, which is equal to the length direction boundary difference multiplied by the width direction boundary difference multiplied by the height direction boundary difference. The length direction boundary difference is the length direction maximum boundary minus the length direction minimum boundary, and the width direction and height direction boundary differences are calculated in a similar manner. The region volume is in units of cubic meters, reflecting the spatial scale of the region. The basic cooling load is defined as the load fluctuation intensity multiplied by the region volume multiplied by the air density. The air density is determined according to the temperature and humidity conditions of the space affected by the cooling system, and the typical value is 1.2 kg / m³, which is applicable under standard atmospheric pressure and room temperature conditions. When the experimental cabin environment deviates from the standard conditions, it needs to be corrected. The basic cooling load is in units of kilowatts, and the calculation process is to multiply the load fluctuation intensity value by the region volume value by the air density value. Since the unit of load fluctuation intensity is kJ / (kg·s), and the units of region volume and air density are kg / m³ and kg / m³ respectively, the unit is converted to kJ / s, i.e. kilowatts. The basic cooling load reflects the heat and humidity load generated by the region itself, without considering the influence of adjacent regions.

[0112] The calculation of the heat and moisture transfer influence quantity considers the coupling relationship between the dynamic load region and each adjacent region. For each adjacent region of the region, the load fluctuation intensity of the adjacent region, the coupling weight between the region and the adjacent region, and the phase difference are obtained. The load fluctuation intensity of the adjacent region is extracted from the local enthalpy difference characteristic table according to the load fluctuation intensity field of the corresponding row indexed by the adjacent region number. The instantaneous influence quantity of the adjacent region on the region is calculated as the load fluctuation intensity of the adjacent region multiplied by the coupling weight. The absolute value of the coupling weight reflects the influence intensity, and the sign reflects the influence direction. A positive value indicates that an increase in the load fluctuation of the adjacent region will lead to an increase in the load of the region, and a negative value indicates the opposite effect. The phase difference is used for time offset correction. When the phase difference is positive, it means that the load fluctuation of the adjacent region lags behind the region, and when calculating the influence quantity of the region at the current time, the load fluctuation intensity of the adjacent region at the past time should be used. The specific time offset is the current time index minus the phase difference value, which gives the index of the adjacent region at a certain past time, and the value at that time is extracted from the load fluctuation intensity time series of the adjacent region. When the phase difference is negative, it means that the load fluctuation of the adjacent region leads the region, and the load fluctuation intensity of the adjacent region at the future time should be used. If the index after time offset exceeds the time series range, the influence quantity of the adjacent region on the region at the current time is set to zero. The influence quantities of all adjacent regions after time offset correction are added up to get the total value of the heat and moisture transfer influence quantity. The accumulation process is to traverse all adjacent regions of the region, calculate the influence quantity of each adjacent region, and sum them up. The unit of heat and moisture transfer influence quantity is kJ / (kg·s), consistent with the unit of load fluctuation intensity.

[0113] The enthalpy decay rate describes the enthalpy change characteristics of the supply air during the process of reaching the dynamic load area from the supply air outlet. The enthalpy decay rate is calculated by analyzing the temperature and humidity field distribution on the supply air path. For a certain dynamic load area, the spatial distance between the area and the nearest supply air outlet is determined. The spatial distance is calculated as the Euclidean distance between the coordinates of the geometric center of the area and the coordinates of the supply air outlet. The distance formula is the square root of the sum of the squares of the differences in the three coordinate axes. The coordinates of the supply air outlet are determined by the layout of the space affected by the refrigeration system, including the position coordinates and the supply air direction. According to the spatial distance and the enthalpy distribution field, the enthalpy sequence on the supply air path is extracted, which includes the enthalpy at the supply air outlet position and the enthalpy at the geometric center of the area. The enthalpy decay rate is defined as the difference between the enthalpy at the supply air outlet and the enthalpy at the geometric center of the area, divided by the spatial distance, with the unit of kJ / (kg·m). When the enthalpy decay rate is positive, it indicates that the enthalpy of the supply air decreases during the transmission process, and the refrigeration capacity weakens. When the enthalpy decay rate is negative, it indicates that the enthalpy of the supply air increases during the transmission process, and it is affected by the heat source. The typical value range of the enthalpy decay rate is 0 to 2.0 kJ / (kg·m), and the default value is set to 0.8 kJ / (kg·m).

[0114] The nominal refrigeration capacity of the supply air is defined as the difference between the enthalpy at the supply air outlet position and the indoor target enthalpy multiplied by the mass flow rate of the supply air. The enthalpy at the supply air outlet position is extracted from the enthalpy distribution field, and the indoor target enthalpy is set according to the control target, with a typical value of 48 kJ / kg. The mass flow rate of the supply air is determined according to the operating parameters of the supply fan, with a typical value range of 0.5 to 5.0 kg / s. The nominal refrigeration capacity of the supply air is in kilowatts, calculated as the difference in enthalpy multiplied by the mass flow rate, with the unit conversion process being kJ / kg multiplied by kg / s to obtain kJ / s, which is kilowatts. The calculation of the effective refrigeration capacity needs to deduct the capacity loss caused by the enthalpy decay, and the loss amount is equal to the enthalpy decay rate multiplied by the spatial distance multiplied by the mass flow rate of the supply air. The unit of the enthalpy decay rate is kJ / (kg·m), multiplied by the spatial distance in meters to obtain kJ / kg, and then multiplied by the mass flow rate in kg / s to obtain kilowatts. The effective refrigeration capacity is equal to the nominal refrigeration capacity of the supply air minus the loss amount, with the unit of kilowatts. When the enthalpy decay rate is large or the spatial distance is far, the effective refrigeration capacity is significantly lower than the nominal value, which needs to be improved by increasing the supply air volume or optimizing the supply air layout.

[0115] The calculation method of the total heat and moisture load of the area is to add the basic refrigeration load and the heat and moisture transfer impact. The basic refrigeration load is in kilowatts, and the heat and moisture transfer impact is in kJ / (kg·s). To ensure unit consistency, the heat and moisture transfer impact needs to be multiplied by the area volume and air density to convert it to kilowatts. The conversion formula is the heat and moisture transfer impact multiplied by the area volume multiplied by the air density, resulting in a transfer impact power with a unit of kilowatts. The total heat and moisture load of the area is equal to the basic refrigeration load plus the transfer impact power, with a unit of kilowatts, reflecting the actual heat and moisture load level of the area considering the coupling effect of adjacent areas.

[0116] The energy balance calculation compares the total heat and humidity load of the area with the effective refrigeration capacity to obtain a modified refrigeration load. The modified refrigeration load is defined as the difference between the total heat and humidity load of the area and the effective refrigeration capacity. When the difference is positive, it means that the heat and humidity load of the area exceeds the effective refrigeration capacity of the existing supply air, and the refrigeration capacity needs to be increased or the supply air parameters need to be adjusted to achieve energy balance. The positive value of the modified refrigeration load is the refrigeration capacity that needs to be supplemented, with the unit being kilowatts. When the difference is negative, it means that the effective refrigeration capacity of the area exceeds the heat and humidity load, and there is a risk of supercooling. The negative value of the modified refrigeration load represents the refrigeration capacity that can be reduced. When the modified refrigeration load is zero, it means that the area has reached an energy balance state and does not need to be adjusted. The absolute value of the modified refrigeration load reflects the degree of deviation of the energy balance of the area. The larger the value, the more serious the deviation, and the higher the control priority.

[0117] The modified refrigeration load calculation method based on energy balance of the present application comprehensively considers the basic refrigeration load, the influence of heat and humidity transfer between areas, and the supply air attenuation factor, realizes accurate estimation of actual refrigeration demand, makes up for the deficiency of traditional load calculation ignoring transmission loss, ensures that the allocated refrigeration resources can meet the actual demand of the area, and avoids the problems of local supercooling or insufficient refrigeration.

[0118] In an optional embodiment, the system energy consumption and enthalpy uniformity of each area are calculated based on the modified refrigeration load as optimization targets, the target weight coefficient is dynamically adjusted according to the load fluctuation intensity, and the steps of solving the refrigeration resource allocation parameters include:

[0119] A state vector containing the load fluctuation intensity and the coupling characteristics of the area is constructed, the target weight coefficient of each dynamic load area is generated based on the reinforcement learning strategy network according to the state vector, and the generation of the target weight coefficient is based on the collaborative adjustment strategy of the topological coupling network;

[0120] Based on the coupling weight and phase difference in the topological coupling network, an influence matrix representing the influence degree and time delay of heat and humidity transfer between areas is constructed;

[0121] Taking the supply air parameters of each dynamic load area as decision variables, taking the sum of the system energy consumption and the enthalpy uniformity index of each area weighted by the target weight coefficient as the optimization target, and taking the heat and humidity transfer constraint between areas defined by the influence matrix and the physical capacity constraint of the refrigeration equipment as the constraint condition, a multi-objective optimization function is constructed, and the refrigeration resource allocation parameters of each dynamic load area are obtained by solving the multi-objective optimization function;

[0122] The refrigeration resource allocation parameters are executed, the actual running effect is monitored, the parameters of the reinforcement learning strategy network are updated according to the monitoring results, and the online learning and continuous optimization of the target weight coefficient adjustment strategy are realized.

[0123] Exemplarily, the construction of the state vector includes two parts: the load fluctuation intensity feature and the region coupling feature. The load fluctuation intensity feature is extracted from the local enthalpy difference feature. For the six dynamic load regions, the load fluctuation intensity values of the regions are extracted to form a load fluctuation intensity vector with a length of 6. The region coupling feature is extracted from the topological coupling network, including the neighbor number, the average coupling weight, the maximum coupling weight, and the average phase difference absolute value of each region. The neighbor number is defined as the total number of adjacent regions connected in the topological coupling network from the neighbor list length of the region node. The average coupling weight is calculated as the arithmetic average of the coupling weights between the region and all adjacent regions, and is set to zero when the region has no adjacent region. The maximum coupling weight is the maximum value of the coupling weight absolute values between the region and all adjacent regions. The average phase difference absolute value is calculated as the arithmetic average of the phase difference absolute values between the region and all adjacent regions. For the six dynamic load regions, the neighbor number, the average coupling weight, the maximum coupling weight, and the average phase difference absolute value of each region form vectors with a length of 6, respectively. The state vector splices the load fluctuation intensity vector and the four region coupling feature vectors to obtain a one-dimensional vector with a length of 30. The elements of the vector are floating-point numbers, which are normalized to the interval of zero to one. The normalization method is to divide the load fluctuation intensity by the maximum value of the load fluctuation intensity, which is 0.5 kJ / (kg·s), to divide the coupling weight by 1.0, and to divide the phase difference absolute value by the time series length 30. The state vector is stored in the form of a one-dimensional array. The array indexes 0 to 5 correspond to the load fluctuation intensity of each region, the indexes 6 to 11 correspond to the neighbor number, the indexes 12 to 17 correspond to the average coupling weight, the indexes 18 to 23 correspond to the maximum coupling weight, and the indexes 24 to 29 correspond to the average phase difference absolute value.

[0124] The reinforcement learning policy network is implemented with a deep neural network, which receives a state vector as input and outputs target weight coefficients for each dynamic load zone. The policy network contains three fully connected layers. The first layer has an input dimension of 30 and 64 neurons with a rectified linear unit (ReLU) activation function. The second layer has 32 neurons with a ReLU activation function. The third layer has an output dimension of 12, which includes the system energy weight and enthalpy uniformity weight for each of the six zones. The third layer uses a Softmax function as the activation function to ensure that the sum of the two weight coefficients for each zone is 1.0 and that all weight coefficients are positive. The Softmax function groups and normalizes the 12 outputs of the third layer. For zone i, the system energy weight coefficient and the enthalpy uniformity weight coefficient are calculated by taking the exponential of the corresponding output values and dividing by the sum of the two exponential values. The weight parameters of the policy network are initialized using the Xavier method, and the bias terms are initialized to 0.01. The collaborative adjustment policy is implemented based on a topological coupling network. For pairs of adjacent zones with a coupling weight absolute value greater than the collaboration threshold, the target weight coefficients are adjusted to be close to each other. The collaboration threshold is set to 0.7. When the absolute value of the coupling weight between two zones exceeds 0.7, the average of the system energy weight coefficient and the enthalpy uniformity weight coefficient for each zone is calculated, and the average is used as the adjusted weight coefficient, which replaces the original output of the policy network. The collaborative adjustment ensures that strongly coupled zones use similar control strategies and avoids oscillations caused by local optimization.

[0125] The influence matrix represents the degree of influence and time delay of heat and moisture transfer between zones. The influence matrix is a square matrix with the same number of rows and columns as the number of dynamic load zones. For six zones, the influence matrix has a dimension of 6x6. The element (i,j) of the matrix represents the influence coefficient of zone j on zone i. The influence coefficient is calculated based on the coupling weight and phase difference in the topological coupling network. When zone j is an adjacent zone of zone i, the influence coefficient is equal to the absolute value of the coupling weight multiplied by the time delay attenuation factor. The time delay attenuation factor is determined based on the absolute value of the phase difference and is calculated as 1.0 minus the absolute value of the phase difference divided by the time series length 30, multiplied by the attenuation coefficient 0.5. The attenuation coefficient reflects the weakening of the influence strength due to time delay, with a typical value range of 0.3 to 0.7 and a default value of 0.5. When zone j is not an adjacent zone of zone i, the influence coefficient is set to zero. The diagonal elements of the influence matrix are set to zero, indicating that a zone is not influenced by itself. The influence matrix is stored as a two-dimensional array with floating-point elements ranging from zero to one. The influence matrix is used to construct the constraint conditions to ensure that the control strategy considers the interaction between zones.

[0126] The decision variables of the multi-objective optimization function are the supply parameters of each dynamic load zone, including supply temperature, supply humidity, and supply mass flow rate. For six zones, the total number of decision variables is 18, with three supply parameters for each zone. The supply temperature ranges from 12 to 18 degrees Celsius, the supply humidity ranges from 0.006 to 0.010 kg / kg of dry air, and the supply mass flow rate ranges from 0.5 to 5.0 kg / s. The decision variables are represented in a one-dimensional vector form, with a vector length of 18. Indices 0 to 2 correspond to the supply temperature, supply humidity, and supply mass flow rate of Zone 1, indices 3 to 5 correspond to Zone 2, and so on. The optimization objectives include system energy consumption and enthalpy uniformity index. System energy consumption is defined as the sum of the energy consumption of all zone supply fans and the energy consumption of the chiller unit. The energy consumption of the supply fan is directly proportional to the cube of the supply mass flow rate, with a proportionality coefficient of 0.0015 kW / (kg / s) cubed. The energy consumption of the chiller unit is directly proportional to the total cooling capacity, with a proportionality coefficient of the inverse of the energy efficiency ratio, which is typically 3.5, resulting in a proportionality coefficient of 0.286. The total cooling capacity is equal to the sum of the modified cooling loads of all zones, which are recalculated based on the current supply parameters. The actual enthalpy of each zone is calculated based on the temperature and humidity distribution field of that zone, taking the arithmetic mean of the enthalpy values of all grid points within the boundary range of the zone. The target enthalpy is set to 48.0 kJ / kg. The enthalpy uniformity index is defined as the square root of the sum of the squared deviations of the actual enthalpy of each zone from the target enthalpy. The unit of the enthalpy uniformity index is kJ / kg, and a smaller value indicates that the enthalpy of each zone is closer to the target value and the distribution is more uniform. The optimization objective function is the product of the system energy consumption and its target weight coefficient plus the product of the enthalpy uniformity index and its target weight coefficient. The total objective function value is obtained by summing the weighted objectives of the six zones. The unit of the objective function value is kilowatts, and since the enthalpy uniformity index and system energy consumption have different dimensions, they need to be dimensionally normalized. The enthalpy uniformity index is multiplied by a normalization coefficient of 10.0 to make its numerical range comparable to that of the system energy consumption.

[0127] The constraints include inter-zone heat and moisture transfer constraints and physical capability constraints of the refrigeration equipment. The inter-zone heat and moisture transfer constraints are defined based on an influence matrix, which requires that the actual cooling effect of a zone is equal to the modified cooling load of the zone minus the heat and moisture transfer effects of all adjacent zones on the zone. The heat and moisture transfer effects are calculated as the sum of the elements of the corresponding row of the influence matrix multiplied by the modified cooling load of each zone. The constraints ensure that the control strategy considers the coupling relationship between zones when allocating the cooling resources, and avoids a situation where the cooling amount of a zone is sufficient, but the actual effect is insufficient due to the heat and moisture transfer of adjacent zones. The physical capability constraints of the refrigeration equipment include upper and lower limits of the supply air temperature, upper and lower limits of the supply air humidity, upper and lower limits of the supply air mass flow rate, an upper limit of the total supply air flow rate, and an upper limit of the refrigeration unit capacity. The upper and lower limits of the supply air temperature are 12 to 18 degrees Celsius, the upper and lower limits of the supply air humidity are 0.006 to 0.010 kg / kg of dry air, and the upper and lower limits of the supply air mass flow rate are 0.5 to 5.0 kg / s. The upper limit of the total supply air flow rate is that the sum of the supply air mass flow rates of all zones does not exceed 20.0 kg / s, which is derived from the total power constraint of the supply air fan. The upper limit of the refrigeration unit capacity is that the total cooling capacity does not exceed 50.0 kW, which is derived from the rated power of the refrigeration unit. The constraints are expressed in the form of inequalities, and the optimization solver checks whether the decision variables satisfy all the constraints. Solutions that do not satisfy the constraints are rejected.

[0128] The solution of the multi-objective optimization function uses a particle swarm optimization algorithm. The particle swarm contains 50 particles, each representing a set of supply air parameter configurations. The position vector of a particle has a dimension of 18, corresponding to the 18 decision variables. The velocity vector of a particle also has a dimension of 18, representing the rate of change of the position. When the particle swarm is initialized, the positions of each particle are randomly sampled within the value range of the decision variables, and the velocity is initialized to zero. During the iteration process of the particle swarm, each particle updates its velocity and position based on its historical optimal position and the global optimal position. The historical optimal position of a particle is the position that achieves the minimum objective function value in the previous iterations, and the global optimal position is the position that achieves the minimum objective function value among all particles in the previous iterations. The velocity update formula includes three parts: the inertia term, the cognitive term, and the social term. The inertia term is the current velocity multiplied by the inertia weight 0.7, the cognitive term is the difference between the historical optimal position and the current position multiplied by the cognitive coefficient 1.5 and then multiplied by a random number between zero and one, and the social term is the difference between the global optimal position and the current position multiplied by the social coefficient 1.5 and then multiplied by a random number between zero and one. The sum of the three terms is the updated velocity. The position update is the current position plus the updated velocity. The velocity is limited between the negative maximum velocity and the positive maximum velocity, and the maximum velocity is set to 20% of the value range of the decision variables. After the position update, the constraints need to be checked. If the constraints are violated, the position is adjusted to the constraint boundary. The particle swarm iterates 100 times, and the objective function value of all particles is calculated in each iteration to update the historical optimal position and the global optimal position. After the iteration is terminated, the global optimal position is taken as the optimization solution, i.e., the allocation parameters of the cooling resources of each dynamic load zone.

[0129] The refrigeration resource allocation parameters include the supply air temperature, supply air humidity, and supply air mass flow rate of each zone, and are output in vector form. When the refrigeration resource allocation parameters are executed, the supply air parameters are sent to the air conditioning controller through the communication interface. The communication interface uses the Modbus protocol, the controller address and register mapping are pre-configured, the supply air temperature is written to the temperature set point register, the supply air humidity is written to the humidity set point register, the supply air mass flow rate is achieved by adjusting the air valve opening, and the mapping relationship between the air valve opening and the flow rate is determined through calibration. After the parameter sending is completed, the air conditioning system adjusts the supply fan speed, the refrigeration unit load, and the air valve opening according to the new set point to achieve the target supply air parameters. The parameter execution delay is the control cycle time interval of 60s, and the actual running effect is monitored after execution.

[0130] The actual running effect is monitored by collecting the temperature and humidity sensor data of each zone. The monitoring period is 60s, corresponding to one control cycle. The monitoring indicators include the average temperature, average humidity, average enthalpy of each zone, and the total energy consumption of the system. The average temperature and average humidity of each zone are calculated from the arithmetic mean of all sensor readings in the zone, and the average enthalpy is calculated based on the average temperature and average humidity. The total energy consumption of the system is obtained by accumulating the readings of the supply fan power meter and the refrigeration unit power meter. The monitoring indicators are compared with the optimization target, the deviation between the actual system energy consumption and the predicted system energy consumption is calculated, and the deviation between the actual enthalpy uniformity index and the predicted enthalpy uniformity index is calculated. The deviation is used as the reward signal of reinforcement learning, and the reward value is defined as the weighted sum of negative deviation. The energy consumption deviation is multiplied by the energy consumption penalty coefficient 0.5, the enthalpy uniformity deviation is multiplied by the uniformity penalty coefficient 5.0, and the sum of the two is taken as the reward. The smaller the deviation, the greater the reward. When the actual running effect is better than the prediction, the reward is positive, and when the actual running effect is worse than the prediction, the reward is negative.

[0131] The policy gradient method is adopted for the reinforcement learning strategy network parameter update, and the advantage function is combined to improve the training stability. The system stores the state vector of the current control period, the target weight coefficient, the distribution parameter, and the reward value as experience samples into a replay buffer with a capacity of 1000. When 32 samples are accumulated, random sampling is performed for batch training. The loss function includes two parts, policy loss and value loss. The policy network and the value network share the first two network structures. The policy loss is calculated as the logarithmic probability of the weight coefficient multiplied by the negative mean value of the advantage value, and the value loss is the mean square error of the value estimate and the actual reward. The advantage value is equal to the actual reward minus the value estimate, and a positive value indicates better performance than the average level. The total loss function combines the two parts, and the value loss weight coefficient is set to 0.5. The network parameter update uses the Adam optimizer, and the learning rate is 0.0001. The optimizer calculates the gradient through backpropagation, maintains the first and second moment estimates of each parameter, and updates the parameter step based on the learning rate and the moment estimate. The system saves a checkpoint containing network weights and optimizer state every 10 batches of training, and triggers a rollback mechanism when a new version performs poorly on the validation set. The online learning process forms a complete closed loop: the policy network continuously receives the state vector, generates the weight coefficient, executes the optimization solution, monitors the effect, and updates the parameters, realizing the continuous optimization of the weight adjustment strategy.

[0132] The reinforcement learning-based target weight coefficient dynamic adjustment method of the present application realizes the adaptive balance between system energy consumption and enthalpy uniformity. The influence matrix introduced incorporates the inter-regional heat and moisture transfer constraints into the optimization framework, ensuring the physical feasibility of the optimization results. The closed-loop feedback mechanism and online learning capability enable the system to continuously accumulate operating experience and continuously optimize the control strategy, achieving adaptive optimization and long-term performance improvement of refrigeration resource distribution.

[0133] In an alternative embodiment, the step of generating the target weight coefficient based on the topology coupling network using a collaborative adjustment strategy comprises:

[0134] Based on the distance relationship and connection mode between regions in the topology coupling network, and combined with the enthalpy gradient change trajectory of each region, a collaborative adjustment matrix is constructed, and the elements in the collaborative adjustment matrix represent the influence relationship of the target weight coefficient adjustment between regions;

[0135] According to the element value of the collaborative adjustment matrix, the region groups are divided, the regions with an influence relationship exceeding a preset influence threshold are determined as a collaborative region group, the master-slave relationship of the regions within the group is determined based on the size of the influence relationship, and an initial adjustment scheme of the target weight coefficient of the collaborative region is generated;

[0136] According to the synergic adjustment matrix, a change trend of enthalpy gradient in the synergic region group is calculated, and when the change trend of the enthalpy gradient exceeds a preset range, a target weight coefficient of a corresponding region in the initial adjustment scheme is adjusted until the change trend of the enthalpy gradient meets a preset range requirement.

[0137] The adjustment scheme meeting the preset range requirement is taken as a final target weight coefficient, and is used for solving the multi-objective optimization function.

[0138] For example, the construction of the synergic adjustment matrix is based on the distance relationship between regions, the connection mode and the enthalpy gradient change trajectory of each region in the topological coupling network. The distance relationship between regions is obtained by calculating the Euclidean distance between the geometric center coordinates of each dynamic load region. For region i and region j, the distance is calculated as the square root of the sum of the square of the coordinate difference in length, width and height of the geometric center of the two regions. The distance matrix is a 6x6 symmetric matrix, and the matrix element (i, j) represents the spatial distance between region i and region j, with the unit of meters, and the diagonal element is zero. The connection mode is extracted from the adjacency matrix of the topological coupling network, and the adjacency matrix element is a Boolean value. When two regions are spatially adjacent, the corresponding element is true, otherwise it is false. The enthalpy gradient change trajectory includes the evolution data of the enthalpy time gradient and the spatial gradient of each region in the control period. For a certain region, the average value of the enthalpy time gradient of the region at 30 time points is extracted to form a time sequence with a length of 30, which reflects the change trend of the enthalpy of the region with time. The difference between the enthalpy time gradient of two adjacent time points is calculated to obtain the enthalpy gradient change rate sequence with a length of 29. The standard deviation of the enthalpy gradient change rate is taken as the activity index of the enthalpy gradient change of the region. The larger the standard deviation is, the more intense the enthalpy gradient change of the region is.

[0139] The calculation of the cooperative adjustment matrix element (i, j) integrates the distance relationship, connection mode, and enthalpy gradient change activity. When region i and region j are adjacent in the adjacency matrix, the basic influence coefficient is set as the absolute value of the corresponding coupling weight in the topological coupling network. When region i and region j are not adjacent, the basic influence coefficient is set to zero. The distance correction factor is calculated according to the spatial distance between the two regions. The calculation method is 1.0 divided by 1.0 plus the distance value divided by the distance attenuation constant. The distance attenuation constant is set to 5.0 m, which reflects the weakening degree of the distance on the influence relationship. The distance correction factor ranges from zero to one, and the greater the distance, the smaller the correction factor. The activity correction factor is calculated according to the geometric mean of the enthalpy gradient change activity of the two regions. The geometric mean is the square root of the product of the standard deviations of the two regions. The activity correction factor is obtained by normalizing the geometric mean to the interval of zero to one, and the normalization divisor is the maximum value of the standard deviation of all regions. The cooperative adjustment matrix element is equal to the basic influence coefficient multiplied by the distance correction factor and then multiplied by the activity correction factor. The cooperative adjustment matrix is a 6x6 asymmetric matrix, and the element value ranges from zero to one. The greater the element value, the stronger the influence of region j on the target weight coefficient adjustment of region i. The cooperative adjustment matrix is stored in the form of a two-dimensional array, and the array elements are floating-point numbers.

[0140] The division of the region group is based on the comparison of the cooperative adjustment matrix element value and the preset influence threshold. The preset influence threshold is set to 0.5, which is determined according to the region coupling strength distribution of the refrigeration system action space. The typical value range is 0.4 to 0.7. Traverse all elements of the cooperative adjustment matrix. When the value of element (i, j) exceeds the preset influence threshold, region i and region j are determined as a region pair with cooperative relationship. The union-find set data structure is used to manage the cooperative relationship. When the union-find set is initialized, each region is an independent set. When region i and region j are found to have a cooperative relationship, the sets in which the two regions are located are merged. After traversal is completed, each connected component in the union-find set corresponds to a cooperative region group. The number of cooperative region groups is equal to the number of root nodes in the union-find set. Each cooperative region group contains one or more dynamic load regions. For a cooperative region group containing a single region, the region independently adjusts the target weight coefficient and is not affected by other regions. For a cooperative region group containing multiple regions, the regions in the group need to be adjusted cooperatively.

[0141] The master-slave relationship of the regions in the collaborative region group is determined based on the size of the influence relationship. For a collaborative region group, the relevant rows and columns of all regions in the group in the collaborative adjustment matrix are extracted to form a sub-matrix. The sum of the elements of each row of the sub-matrix is calculated, and the row element sum represents the total strength of the region affected by other regions in the group. The sum of the elements of each column of the sub-matrix is calculated, and the column element sum represents the total strength of the region affecting other regions in the group. The dominance degree is defined as the column element sum minus the row element sum. The positive value of the dominance degree indicates that the influence of the region on other regions is greater than the influence of other regions on the region, and the region is in a dominant position in the collaborative region group. The negative value of the dominance degree indicates that the region is in a subordinate position in the collaborative region group. All regions in the collaborative region group are sorted in descending order of the dominance degree, and the region with the largest dominance degree is determined as the master region, and the other regions are determined as the slave regions. The target weight coefficient adjustment priority of the master region is the highest, and the target weight coefficient adjustment of the slave region needs to refer to the adjustment direction of the master region.

[0142] The initial adjustment scheme of the target weight coefficient of the collaborative region is generated based on the output of the reinforcement learning policy network. For a collaborative region group, the initial target weight coefficients of each region in the group are obtained, and the initial target weight coefficients are the system energy weight coefficient and the enthalpy uniformity weight coefficient directly output by the policy network. For the master region, the initial adjustment scheme keeps the target weight coefficients output by the policy network unchanged. For the slave region, the initial adjustment scheme moves the target weight coefficient of the slave region towards the target weight coefficient of the master region. The convergence coefficient is determined according to the influence coefficient of the slave region to the master region in the collaborative adjustment matrix, and the convergence coefficient is equal to the influence coefficient multiplied by the collaborative strength parameter. The collaborative strength parameter is set to 0.6, and the typical value range is 0.4 to 0.8. The adjusted system energy weight coefficient of the slave region is calculated as the initial system energy weight coefficient of the slave region multiplied by 1.0 minus the convergence coefficient, plus the system energy weight coefficient of the master region multiplied by the convergence coefficient. The adjusted enthalpy uniformity weight coefficient of the slave region is calculated as 1.0 minus the adjusted system energy weight coefficient, to ensure that the sum of the two weight coefficients is 1.0. The initial adjustment scheme is stored in the form of a dictionary, with the region number as the key and the tuple containing the system energy weight coefficient and the enthalpy uniformity weight coefficient as the value.

[0143] For a collaborative region group, the average of the enthalpy time gradient of each region in the group in the last 10 control cycles is extracted to form a time series with a length of 10. The linear fitting slope of the time series is calculated as the enthalpy gradient change trend. The linear fitting uses the least squares method, and the slope is calculated as the covariance of the time series value and the time index divided by the variance of the time index. The time index is an integer sequence from 0 to 9, the covariance is calculated as the difference between the time series value and its average multiplied by the difference between the time index and its average, summed over 10 times and divided by 10. The variance of the time index is the sum of the squares of the difference between the index value and the average index value divided by 10. The enthalpy gradient change trend slope unit is kJ / (kg·s²), a positive value indicates that the enthalpy time gradient is rising, the load fluctuation intensity is increasing, and a negative value indicates that the enthalpy time gradient is falling, the load fluctuation intensity is decreasing. The weighted average of the enthalpy gradient change trend slope of all regions in the collaborative region group is calculated, the weight is the region volume of each region, and the weighted average slope is used as the overall enthalpy gradient change trend of the collaborative region group.

[0144] The preset range is defined as the allowed interval of the enthalpy gradient change trend slope, the lower limit of the interval is set to -0.005 kJ / (kg·s²), and the upper limit is set to 0.005 kJ / (kg·s²). The range reflects the desired load fluctuation stability state of the system, when the enthalpy gradient change trend slope is within the range, it is considered that the load fluctuation is in a controllable state. When the overall enthalpy gradient change trend slope of the collaborative region group exceeds the preset range, the target weight coefficient of the corresponding region in the initial adjustment scheme needs to be adjusted. The adjustment strategy is determined according to the positive and negative directions of the enthalpy gradient change trend slope, when the slope is positive and exceeds the upper limit, it indicates that the load fluctuation intensity increases too fast, the enthalpy uniformity weight coefficient needs to be increased to strengthen the uniformity control, and the system energy consumption weight coefficient needs to be reduced. When the slope is negative and below the lower limit, it indicates that the load fluctuation intensity decreases too fast leading to excessive refrigeration, the system energy consumption weight coefficient needs to be increased to reduce energy consumption, and the enthalpy uniformity weight coefficient needs to be reduced.

[0145] The adjustment amount of the target weight coefficient is calculated according to the deviation degree of the enthalpy gradient change trend slope from the preset range boundary. The deviation degree is defined as the absolute value of the difference between the actual slope and the preset range boundary. The adjustment step is equal to the deviation degree multiplied by the adjustment sensitivity coefficient, which is set to 20.0, and the typical value range is 10.0 to 30.0. This coefficient reflects the response speed of the weight coefficient to the enthalpy gradient change trend. The value range of the adjustment step is limited to 0.01 to 0.15 to avoid oscillation caused by too large adjustment amplitude in a single adjustment. When it is necessary to increase the enthalpy uniformity weight coefficient, the enthalpy uniformity weight coefficients of all regions in the collaborative region group in the initial adjustment scheme are increased by the adjustment step, and the system energy consumption weight coefficient is decreased by the same adjustment step. When it is necessary to increase the system energy consumption weight coefficient, the adjustment direction is opposite. The adjusted target weight coefficient needs to meet the non-negative constraint and the constraint that the sum is 1.0. If the adjusted weight coefficient is less than zero, it is set to zero, and the other weight coefficient is set to 1.0. If the adjusted weight coefficient is greater than 1.0, it is set to 1.0, and the other weight coefficient is set to zero.

[0146] After the adjustment is completed, the overall enthalpy gradient change trend of the collaborative region group is recalculated. Based on the adjusted target weight coefficient, the multi-objective optimization function is simulated to obtain new refrigeration resource allocation parameters. The enthalpy time gradient of each region is predicted based on the new allocation parameters, and the predicted enthalpy gradient change trend slope is calculated. The prediction method uses a physically constrained neural network inference, which inputs the new supply air parameters and the current time, and outputs the temperature and humidity field distribution in the next 10 control periods. The enthalpy time gradient of each region is calculated, and the change trend slope is extracted. It is checked whether the predicted enthalpy gradient change trend slope meets the preset range requirement. If it meets the requirement, the current adjustment scheme is taken as the final target weight coefficient. If it does not meet the requirement, the adjustment continues. The maximum number of iterations is set to 5. If the preset range requirement is still not met after 5 iterations, the adjustment scheme with the enthalpy gradient change trend slope closest to the preset range center value of zero in the 5 iterations is selected as the final target weight coefficient.

[0147] The final target weight coefficient is used to replace the target weight coefficient directly output by the reinforcement learning strategy network when solving the multi-objective optimization function. For regions belonging to the collaborative region group, the target weight coefficient adjusted by cooperation is used. For independent regions not belonging to any collaborative region group, the target weight coefficient directly output by the strategy network is used. The target weight coefficient is passed to the multi-objective optimization function in the form of a dictionary, with the region number as the key and the tuple containing the system energy consumption weight coefficient and the enthalpy uniformity weight coefficient as the value. The multi-objective optimization function constructs a weighted objective function according to the target weight coefficient and solves the refrigeration resource allocation parameters by the particle swarm optimization algorithm.

[0148] The present application can quantitatively characterize the mutual influence degree of inter-regional weight adjustment by constructing a synergistic adjustment matrix. Based on the synergistic region group and master-slave relationship determination mechanism divided according to the influence relationship, the weight adjustment can be carried out on the premise of considering the global influence. The prediction and verification mechanism of enthalpy gradient change trend is introduced to ensure that the weight adjustment will not lead to the imbalance of enthalpy distribution between regions. This synergistic adjustment method avoids the optimization conflict problem caused by traditional independent adjustment, and improves the overall coordination and operation stability of the multi-region refrigeration system.

[0149] The actual temperature and humidity data after the refrigeration resource allocation parameter is executed is collected, and the implementation of adjusting the dynamic load region boundary according to the deviation type of the actual enthalpy distribution and the expected enthalpy distribution is as follows: the actual temperature and humidity data after the refrigeration resource allocation parameter is executed is collected by the temperature and humidity sensor network distributed in the refrigeration space with a period of 60 seconds. According to the collected temperature and humidity data, the actual enthalpy of each point is calculated using the enthalpy calculation formula, and the enthalpy is equal to the sum of the dry air enthalpy and the water vapor enthalpy, the dry air enthalpy is the temperature multiplied by the specific heat capacity coefficient 1.005, and the water vapor enthalpy is the humidity multiplied by the latent heat of vaporization 2501kJ / kg. The actual enthalpy calculated is compared with the expected enthalpy predicted by the physical constraint neural network to generate an enthalpy deviation field. The system identifies three typical deviation types in the deviation field: boundary drift type, which is characterized by steep deviation gradient at the region boundary; area change type, which is characterized by overall region shrinkage or expansion; and splitting and merging type, which is characterized by single region splitting or multiple region fusion. For the boundary drift type deviation, the boundary coordinates are corrected by moving, and the moving amount is proportional to the gradient direction; for the area change type deviation, the six boundary values of the region bounding box are adjusted in proportion to realize; for the splitting and merging type deviation, the region division algorithm based on clustering is re-executed. After the boundary adjustment is completed, the region boundary information library is updated and used as the initial boundary of the next control period. The boundary adjustment frequency is dynamically adjusted according to the system stability, and the adjustment is performed once every 10 control periods in the stable state and once every 2-3 periods in the unstable state, so that the boundary reflects the actual load distribution.

[0150] In a second aspect of the embodiments of the present application, a refrigeration system adaptive operation strategy optimization system is provided, comprising:

[0151] A first unit is configured to acquire temperature and humidity data and fluid boundary parameters in a refrigeration system action space;

[0152] A second unit is configured to reconstruct a continuous temperature and humidity distribution field based on a physical constraint neural network, calculate an enthalpy distribution field and extract a gradient change trajectory, divide dynamic load regions according to gradient direction turning points and gradient periodic fluctuation points, and extract load fluctuation intensity of each region as a local enthalpy difference feature;

[0153] The third unit is configured to calculate a mutual correlation coefficient and a phase difference of a load fluctuation intensity time sequence of adjacent areas to construct a topological coupling network, calculate an enthalpy value decay rate when air supply reaches each area according to the fluid boundary parameter, and perform energy balance calculation on the local enthalpy difference feature, the coupling weight of the topological coupling network, and the enthalpy value decay rate to obtain a corrected refrigeration load for each area.

[0154] The fourth unit is configured to calculate system energy consumption and enthalpy value uniformity of each area as an optimization target based on the corrected refrigeration load, dynamically adjust a target weight coefficient according to the load fluctuation intensity, and solve to obtain refrigeration resource allocation parameters.

[0155] The fifth unit is configured to collect actual temperature and humidity data after the refrigeration resource allocation parameters are executed, and adjust a dynamic load area boundary according to a deviation type of an actual enthalpy value distribution and an expected enthalpy value distribution.

[0156] The third aspect of the embodiment of the present application provides an electronic device, comprising:

[0157] a processor;

[0158] a memory for storing processor-executable instructions;

[0159] The processor is configured to call the instructions stored in the memory to execute the method described above.

[0160] The fourth aspect of the embodiment of the present application provides a computer-readable storage medium, which stores computer program instructions, and the computer program instructions are executed by a processor to implement the method described above.

[0161] The present application can be a method, device, system and / or computer program product. The computer program product can include a computer readable storage medium having computer readable program instructions stored therein, which are used to perform various aspects of the present application.

Claims

1. A method for optimizing the adaptive operation strategy of a refrigeration system, characterized in that, include: Acquire temperature and humidity data, and fluid boundary parameters within the operating space of the refrigeration system; The temperature and humidity data are reconstructed using a physical constraint neural network to obtain a continuous temperature and humidity distribution field. The enthalpy distribution field is calculated and the gradient change trajectory is extracted. Dynamic load regions are divided according to the gradient direction inflection points and gradient periodic fluctuation points. The load fluctuation intensity of each region is extracted as a local enthalpy difference feature. Calculate the cross-correlation coefficient and phase difference of the load fluctuation intensity time series of adjacent areas to construct a topologically coupled network; calculate the enthalpy decay rate of the supplied air when it reaches each area based on the fluid boundary parameters; For each region, the corrected cooling load is obtained by energy balance calculation based on local enthalpy difference characteristics, coupling weights of the topological coupling network, and enthalpy decay rate. Based on the energy consumption of the modified cooling load calculation system and the uniformity of enthalpy values ​​in each region as optimization objectives, the target weight coefficients are dynamically adjusted according to the intensity of load fluctuations to obtain the cooling resource allocation parameters. Collect actual temperature and humidity data after executing cooling resource allocation parameters, and adjust the dynamic load zone boundary according to the deviation type between the actual enthalpy distribution and the expected enthalpy distribution.

2. The method according to claim 1, characterized in that, The steps of reconstructing the temperature and humidity data based on a physically constrained neural network to obtain a continuous temperature and humidity distribution field, calculating the enthalpy distribution field, and extracting the gradient change trajectory include: The physical constraint neural network includes a temperature prediction branch and a humidity prediction branch. The temperature prediction branch receives spatial location and time information as input, and the humidity prediction branch receives spatial location, time information, and hidden layer features of the temperature prediction branch as input. The loss function of the physical constraint neural network includes a data fitting term and a physical constraint term. The physical constraint term includes thermodynamic conservation constraints and time causality constraints. The time causality constraint calculates the enthalpy value based on the outputs of the temperature prediction branch and the humidity prediction branch, and constrains the sign of the time derivative of the enthalpy value to be consistent with the direction of energy inflow and outflow. The weight ratio of the data fitting term and the physical constraint term is dynamically adjusted based on the density of temperature and humidity measuring points. Based on the physical constraint neural network, the temperature and humidity data are reconstructed to obtain a continuous temperature and humidity distribution field. The enthalpy distribution field is calculated, and the temporal gradient and spatial gradient of the enthalpy distribution field are obtained. The evolution process of the temporal gradient and spatial gradient in the spatiotemporal domain is taken as the gradient change trajectory.

3. The method according to claim 1, characterized in that, The steps for dividing dynamic load regions based on gradient direction inflection points and gradient periodic fluctuation points, and extracting the load fluctuation intensity of each region as a local enthalpy difference characteristic, include: The enthalpy spatial gradient in the gradient change trajectory is oriented, and the spatial position where the rate of change of the gradient direction vector angle exceeds a first set threshold is identified as the gradient direction inflection point. A time-domain analysis is performed on the enthalpy time gradient in the gradient change trajectory to identify spatial locations where the rate of change of the gradient amplitude at adjacent time points exceeds a second set threshold and the direction of change reverses as gradient periodic fluctuation points. Cluster analysis is performed on the spatial distribution of the gradient direction inflection points and the gradient periodic fluctuation points to obtain multiple spatial clusters, and the boundary of each spatial cluster is used as the boundary of the dynamic load region. The fluctuation characteristics of the enthalpy time gradient within each dynamic load region are calculated as the load fluctuation intensity, and combined with regional location information to form local enthalpy difference characteristics.

4. The method according to claim 1, characterized in that, The steps for calculating the cross-correlation coefficient and phase difference of the load fluctuation intensity time series of adjacent regions and constructing the topologically coupled network include: Based on the load fluctuation intensity of each dynamic load region in the local enthalpy difference characteristics, a time series of load fluctuation intensity for each dynamic load region is formed; To determine the spatial adjacency between dynamic load areas, for two spatially adjacent dynamic load areas, calculate the cross-correlation function of their load fluctuation intensity time series, extract the maximum value of the cross-correlation function as the cross-correlation coefficient, and extract the time offset corresponding to the maximum value of the cross-correlation function as the phase difference; A topologically coupled network is constructed by using each dynamic load region as a network node, the cross-correlation coefficient between spatially adjacent dynamic load regions as the weight of the connecting edge, and the phase difference as the directional attribute of the connecting edge.

5. The method according to claim 1, characterized in that, The steps for calculating the corrected cooling load by using local enthalpy difference characteristics, coupling weights of the topologically coupled network, and enthalpy decay rate to perform energy balance calculations include: For each dynamic load region, the load fluctuation intensity and region boundary information of the region are extracted from the local enthalpy difference characteristics, and the coupling weight and phase difference between the region and its neighboring regions are extracted from the topological coupling network. The base cooling load is calculated based on the intensity of load fluctuations in the region and the region boundary information. Based on the coupling weight and phase difference between the dynamic load area and each adjacent area, the heat and moisture transfer influence of the adjacent area on the area is calculated. The heat and moisture transfer influence is equal to the product of the load fluctuation intensity of the adjacent area and the coupling weight, and the accumulated value after time offset correction according to the phase difference. The effective cooling capacity when the supply air reaches the area is calculated based on the enthalpy decay rate of the dynamic load area. The effective cooling capacity is equal to the energy value corresponding to the difference between the nominal cooling capacity of the supply air and the enthalpy decay rate. The total heat and humidity load of the region is obtained by adding the basic cooling load of the dynamic load region to the heat and humidity transfer influence. The energy balance calculation is performed on the total heat and humidity load of the region and the effective cooling capacity to obtain the corrected cooling load of the dynamic load region.

6. The method according to claim 1, characterized in that, The steps for obtaining the cooling resource allocation parameters, based on the energy consumption of the corrected cooling load calculation system and the enthalpy uniformity of each region as optimization objectives, and dynamically adjusting the objective weight coefficients according to the intensity of load fluctuations, include: A state vector containing load fluctuation intensity and regional coupling characteristics is constructed. Based on the reinforcement learning policy network, target weight coefficients for each dynamic load region are generated according to the state vector. The generation of the target weight coefficients is based on the topological coupling network using a cooperative adjustment strategy. Based on the coupling weights and phase differences in the topologically coupled network, an influence matrix is ​​constructed to characterize the degree of influence and time delay of heat and moisture transfer between regions. Using the air supply parameters of each dynamic load region as decision variables, the sum of the system energy consumption after weighting by the target weight coefficient and the enthalpy uniformity index of each region as the optimization objective, and the inter-regional heat and moisture transfer constraints and the physical capacity constraints of the refrigeration equipment defined by the influence matrix as constraints, a multi-objective optimization function is constructed. Solving the multi-objective optimization function yields the refrigeration resource allocation parameters of each dynamic load region. The cooling resource allocation parameters are executed, the actual operating effect is monitored, and the parameters of the reinforcement learning strategy network are updated based on the monitoring results, so as to realize the online learning and continuous optimization of the target weight coefficient adjustment strategy.

7. The method according to claim 6, characterized in that, The steps for generating the target weight coefficients based on the cooperative adjustment strategy of the topologically coupled network include: Based on the distance relationships and connection methods between regions in the topologically coupled network, and combined with the enthalpy gradient change trajectory of each region, a collaborative adjustment matrix is ​​constructed. The elements in the collaborative adjustment matrix represent the influence relationship of the adjustment of the target weight coefficients between regions. The region groups are divided according to the element values ​​of the collaborative adjustment matrix. Regions whose influence relationship exceeds the preset influence threshold are identified as collaborative region groups. The master-slave relationship of the regions within the group is determined based on the magnitude of the influence relationship, and an initial adjustment scheme for the target weight coefficient of the collaborative region is generated. The enthalpy gradient change trend within the collaborative region group is calculated based on the collaborative adjustment matrix. When the enthalpy gradient change trend exceeds the preset range, the target weight coefficient of the corresponding region in the initial adjustment scheme is adjusted until the enthalpy gradient change trend meets the preset range requirement. The adjustment scheme that meets the preset range requirements is used as the final target weight coefficient for solving the multi-objective optimization function.

8. A refrigeration system adaptive operation strategy optimization system, used to implement the method of any one of claims 1-7, characterized in that, include: The first unit is used to acquire temperature and humidity data and fluid boundary parameters within the operating space of the refrigeration system. The second unit is used to reconstruct the temperature and humidity data based on the physical constraint neural network to obtain a continuous temperature and humidity distribution field, calculate the enthalpy distribution field and extract the gradient change trajectory, divide the dynamic load region according to the gradient direction inflection point and the gradient periodic fluctuation point, and extract the load fluctuation intensity of each region as the local enthalpy difference feature. The third unit is used to calculate the cross-correlation coefficient and phase difference of the load fluctuation intensity time series of adjacent areas to construct a topological coupling network; and to calculate the enthalpy decay rate when the air supply reaches each area based on the fluid boundary parameters. For each region, the corrected cooling load is obtained by energy balance calculation based on local enthalpy difference characteristics, coupling weights of the topological coupling network, and enthalpy decay rate. The fourth unit is used to calculate the system energy consumption and the enthalpy uniformity of each region based on the corrected cooling load as optimization objectives. The objective weight coefficients are dynamically adjusted according to the intensity of load fluctuations to obtain the cooling resource allocation parameters. The fifth unit is used to collect actual temperature and humidity data after executing the refrigeration resource allocation parameters, and adjust the dynamic load zone boundary according to the deviation type between the actual enthalpy distribution and the expected enthalpy distribution.

9. An electronic device, characterized in that, include: processor; Memory used to store processor-executable instructions; The processor is configured to invoke instructions stored in the memory to execute the method according to any one of claims 1 to 7.

10. A computer-readable storage medium having computer program instructions stored thereon, characterized in that, When the computer program instructions are executed by the processor, they implement the method described in any one of claims 1 to 7.

Images