A heat pump unit variable temperature difference control method and system based on load prediction

Abstract

The present application belongs to the technical field of heat pump system control, and particularly relates to a heat pump unit variable-temperature difference control method and system based on load prediction, which comprises the following steps: obtaining historical data, extracting environmental features by using B-spline function and covariance matrix characteristic decomposition, and generating periodic features by using trigonometric functions; calculating Euclidean distance to cluster mode clusters and constructing a search space to retrieve neighboring samples in combination with adjacent clusters; constructing a decay factor to correct feature Euclidean distance according to topological distance, inputting the neighboring sample with the minimum weighted distance into a Gaussian process regression model to output predicted load and posterior variance; generating relative variance by the ratio of posterior variance to predicted load square, dynamically updating PID parameters and previous weights, fusing weighted feedforward items and feedback items to output control instructions. The present application quantitatively measures and calculates prediction uncertainty, and improves unit variable load response stability and energy saving effect.

Description

Technical Field

[0001] This invention relates to the field of heat pump system control technology. More specifically, this invention relates to a variable temperature difference control method and system for heat pump units based on load prediction. Background Technology

[0002] The operating status of a heat pump system directly affects the energy consumption level of a building or production line. Its heating and cooling loads are influenced by factors such as outdoor temperature and humidity, unit operating parameters, personnel work and rest schedules, and production rhythm. To achieve refined and low-energy operation of the unit, it is necessary to accurately predict future loads and adjust the supply and return water temperature difference and circulation flow accordingly. Existing heat pump control systems typically employ fixed-parameter feedback control, which can meet basic temperature regulation requirements and has made initial attempts at predictive compensation in some systems, but the overall adjustment accuracy and adaptability still need improvement.

[0003] However, existing load forecasting methods often focus on single-point or shallow feature modeling, failing to adequately represent continuous environmental change curves, unit operating trends, and daily, weekly, and seasonal cyclical patterns. When dealing with large historical data sets and frequent changes in operating conditions, it is difficult to quickly identify truly similar and valuable operating segments, resulting in insufficient adaptability and stability of the forecast results for complex operating conditions. Furthermore, existing heat pump systems often employ fixed-parameter feedback control, which lags in response to sudden changes in the external environment or rapid load variations, easily leading to supply-demand discrepancies, increased energy consumption, and temperature fluctuations. While some systems introduce load forecasting feedforward control, they lack quantitative assessment of forecast uncertainty, making it impossible to dynamically allocate feedforward and feedback weights based on forecast reliability, and difficult to adjust PID control intensity in real time. When forecast deviations are large, control commands may deviate from the reasonable operating point, even inducing system oscillations.

[0004] Referring to Chinese patent application CN108800431A, an energy-saving control method and system for an air source heat pump system are disclosed. This method involves constructing a standard building cooling and heating load prediction model, performing refined building zoning, and generating dynamic curves of cooling and heating energy demand for each area. Based on this, the system performs phased control of the heat pump and its associated chilled water pumps and valve units according to the dynamic energy demand curves. Particularly during the steady-state operation phase, the system judges the difference between the actual temperature difference and the set value, and uses PID temperature difference control to adjust the chilled water pump frequency and the opening degree of related valves, thereby achieving a certain degree of regulation.

[0005] Although the aforementioned application documents proposed and implemented a technical solution combining load forecasting demand curves with PID temperature difference control, its load forecasting mechanism still falls short of the limitations of shallow modeling, failing to fully extract multi-dimensional spatiotemporal features and quickly match similar operating conditions. More importantly, the PID control merely remains at the level of mechanical feedback with fixed parameters, failing to take into account the posterior error of load forecasting. Consequently, it cannot adaptively coordinate feedforward and feedback control in conjunction with forecast uncertainty, and is still highly susceptible to response lag and oscillations when dealing with sudden disturbances. Summary of the Invention

[0006] To address the problems of insufficient load forecast feature representation, low efficiency of similar operating condition retrieval, difficulty in incorporating forecast uncertainty into control, and lag in fixed-proportion integral-derivative feedback response in existing technologies, this invention provides solutions in the following aspects.

[0007] In a first aspect, the present invention provides a variable temperature difference control method for heat pump units based on load forecasting, comprising: S1, acquiring historical datasets composed of heat pump unit operating parameters, environmental parameters, and actual load values, as well as current PID gain parameters and feedforward control weights; dimensionlessly converting continuous time-series data within a predetermined time window and extracting principal component scores through functional principal component analysis to form an environmental feature vector; generating periodic time feature vectors from discrete-time data through sine and cosine encoding; S2, using a self-organizing map network to cluster the periodic time feature vectors into multiple operating mode clusters on a two-dimensional rectangular topological grid network, with each cluster and its adjacent clusters forming a local search space, and controlling the temperature difference of the heat pump unit to be predicted. The periodic time feature vector of time is input into the network to determine the best matching cluster, and k nearest neighbor historical samples are retrieved in the local search space; S3, the topological distance is determined based on the coordinates of each cluster on the grid, the distance of adjacent cluster samples is divided by the attenuation factor negatively correlated with the topological distance, and the predicted load value and the variance of the posterior probability distribution are calculated by using the environmental feature vector of the nearest neighbor historical samples and the actual load value through a Gaussian process regression model; S4, the predicted load value is used as the feedforward signal to calculate the reference temperature difference setpoint, the current PID gain parameter and feedforward control weight are adjusted according to the variance of the posterior probability distribution, the reference temperature difference setpoint is weighted by the feedforward control weight to obtain the weighted feedforward term, and then superimposed with the updated PID feedback compensation term to generate an instruction and issue the instruction to drive the variable frequency motor to adjust the flow rate and perform flow limiting control.

[0008] This invention effectively extracts the environmental change trends in continuous time series and the periodic patterns in discrete time series by integrating functional principal component analysis and sine and cosine phase angle encoding. It employs a self-organizing mapping network combined with a topological distance decay mechanism to accurately screen similar operating conditions within the dimensionality-reduced local search space, improving retrieval efficiency and matching accuracy. By utilizing a Gaussian process regression model, it not only outputs the predicted load but also specifically calculates the variance of the predicted posterior probability distribution, providing a mathematical tolerance basis for subsequent regulation. By dynamically adjusting the feedforward control weights and proportional-integral-derivative feedback parameters, it overcomes the problems of lag in response and deviation of commands from reasonable operating points in traditional fixed-parameter feedback, thereby improving the control stability and energy-saving effect of heat pump units under varying operating conditions.

[0009] Preferably, the step of extracting principal component scores to form an environmental feature vector after dimensionless transformation of continuous time-series data within a predetermined time window using functional principal component analysis includes: standardizing the continuous time-series data of the operating parameters and the environmental parameters respectively; smoothing the standardized continuous time-series data using B-spline basis functions to generate continuous functional data curves; calculating the covariance matrix between all functional data curves and performing eigenvalue decomposition, extracting the eigenvectors corresponding to the first three largest non-zero eigenvalues ​​as the learning basis functions; and calculating the principal component scores of each functional data curve on each basis function through integration, wherein the first... The principal component score of the basis function is the sum of the functional data curve and the first... The time integral of the product of the basis functions; the scores of the first three principal components are concatenated in order to form the environmental feature vector.

[0010] This invention uses B-spline basis functions for smoothing and denoising, combined with eigenvalue decomposition and integral operations, to accurately extract the core evolution laws of continuous time series data, achieving effective dimensionality reduction of high-dimensional time series data and improving the noise resistance and computational efficiency of feature representation.

[0011] Preferably, the step of generating a periodic time feature vector from discrete-time data using sine and cosine encoding includes: extracting month values, weekday values, and hour values ​​from the historical dataset; normalizing the month values, weekday values, and hour values ​​using angular frequency parameters to obtain three normalized angular frequency values; substituting the three angular frequency values ​​into sine and cosine functions respectively to calculate six function result values, and concatenating the six function result values ​​to form the periodic time feature vector.

[0012] This invention utilizes the natural cyclical properties of sine and cosine trigonometric functions to transform increasing discrete-time variables into smooth and continuous two-dimensional correlated coordinates, eliminating abrupt changes in features at time boundaries and accurately reflecting the periodic fluctuation patterns of the data.

[0013] Preferably, the step of clustering periodic time feature vectors into multiple working mode clusters on a two-dimensional rectangular topological grid network using a self-organizing map network includes: initializing the weight vectors of the two-dimensional rectangular topological grid nodes with consistent length and width; calculating the Euclidean distance between each periodic time feature vector and the weight vectors of all grid nodes, and selecting the grid node weight vector with the smallest Euclidean distance as the winning neuron node; and synchronously updating the grid weight vectors of the winning neuron node and its neighboring nodes with a Chebyshev distance of 1 according to the Gaussian neighborhood topological update function and a set learning rate. The Gaussian neighborhood topological update function is... ,in, The neighborhood effect value. The square of the Chebyshev distance between the two nodes within the grid. The neighborhood radius decreases with each iteration; after a set number of iterations, historical samples mapped to the same winning neuron node are incorporated into the same working mode cluster.

[0014] This invention uses a Gaussian neighborhood topology update mechanism to automatically aggregate spatiotemporally related samples into a structured, high-purity working pattern cluster, thereby narrowing the search range for subsequent matching and improving the efficiency of optimization retrieval.

[0015] Preferably, dividing the distance between adjacent cluster samples by a decay factor negatively correlated with the topological distance includes: calculating the Euclidean distance between the center coordinates of the adjacent cluster and the center coordinates of the best-matching cluster within the two-dimensional rectangular topological grid network as the topological distance; the decay factor is determined using a negative exponential function with the natural constant e as the base, the exponent of which is the product of the topological distance value and a preset negative coefficient; calculating the feature space Euclidean distance between the periodic time feature vector of the time to be predicted and each nearest historical sample of the adjacent cluster, dividing the feature space Euclidean distance by the decay factor to obtain the weighted comprehensive distance; and using the weighted comprehensive distance instead of the original distance to determine the proximity of the query match when performing nearest neighbor query matching.

[0016] In the search phase, this invention actively applies a distance penalty amplification mechanism to adjacent cluster samples that cross cluster boundaries, and prioritizes core data with closer topology and more similar features as the prediction benchmark, thus avoiding the mixing of heterogeneous mutation conditions.

[0017] Preferably, the step of calculating the predicted load value and posterior probability distribution variance using the environmental feature vectors of nearest-neighbor historical samples and the actual load value via a Gaussian process regression model includes: establishing a prior distribution with the zero vector as the mean function, and setting a kernel function containing an exponential quadratic kernel function as the covariance kernel function; extracting... The environmental feature vectors of the nearest historical samples form a training input feature matrix. The corresponding actual load values ​​are extracted and standardized using the standard deviation of the training set load to form the observation target vector. Based on the covariance kernel function, the cross-covariance matrix, the autocovariance matrix, and the kernel covariance matrix are calculated. The standardized expected value of the prediction is obtained according to the regression formula and then de-standardized to the predicted load value. The standardized posterior variance is obtained according to the posterior covariance formula and multiplied by the square of the standard deviation of the training set load to obtain the posterior probability distribution variance with the dimension of the load square.

[0018] This invention uses a nonparametric Bayesian framework to deeply fit the nonlinear relationship between features and load, which not only accurately outputs the predicted load value, but also uses posterior covariance to calculate and output a prediction uncertainty index with clear physical dimensions.

[0019] Preferably, adjusting the current PID gain parameter and feedforward control weight based on the posterior probability distribution variance includes: calculating the ratio of the posterior probability distribution variance to the sum of the square of the predicted load value plus a preset small positive number to generate a relative variance; when the relative variance is greater than a set upper limit judgment threshold, reducing the feedforward control weight, and increasing the global proportional gain coefficient and integral increment gain coefficient in the PID gain parameter by a corresponding preset step size and limiting them to a preset upper limit; when the relative variance is less than a set lower limit judgment threshold, increasing the feedforward control weight, and decreasing the global proportional gain coefficient and integral increment gain coefficient by a corresponding preset step size and limiting them to a preset lower limit greater than 0; when the relative variance is between the lower limit judgment threshold and the upper limit judgment threshold, keeping the PID gain parameter and the feedforward control weight unchanged.

[0020] This invention achieves dual-channel dynamic coupling and optimal allocation of feedforward weights and feedback rigid parameters, which strengthens correction and ensures safety when the prediction deviation is large, and improves forward-looking tuning and reduces overshoot risk when the prediction is accurate, effectively compressing the system overshoot and speeding up the stabilization time.

[0021] Preferably, the acquisition of heat pump unit operating parameters, environmental parameters, and actual load values ​​to form a historical dataset includes: acquiring compressor power, inlet and outlet water temperatures, outdoor dry-bulb temperature, outdoor relative humidity, and actual cooling and heating load values ​​at a set sampling frequency to form the historical dataset; performing standard deviation normalization on the operating parameters and environmental parameters as input features, and performing standard deviation normalization on the actual load values ​​before saving the training set load mean and load standard deviation.

[0022] This invention covers key underlying data sources that affect unit heat transfer and energy efficiency, and eliminates scaling differences between heterogeneous parameters through unified dimensionless processing, thereby accelerating the optimization and convergence speed of the underlying algorithm model.

[0023] Preferably, the covariance kernel function is a composite covariance kernel function formed by superimposing the exponential quadratic kernel function and the white noise kernel function, including: optimizing the hyperparameters of the composite covariance kernel function by maximizing the marginal log-likelihood function and using the finite-memory BFGS algorithm; and adding a white noise variance term to the kernel covariance matrix to improve the stability of matrix inversion.

[0024] This invention eliminates the risk of kernel matrix divergence under massive nodes by superimposing white noise to cohesively absorb system measurement and control errors and using the limited memory BFGS optimization algorithm to efficiently find hyperparameters, thus ensuring the numerical stability of the algorithm.

[0025] Secondly, the present invention provides a load-predictive-based heat pump unit temperature difference control system, including a processor and a memory, wherein the memory stores computer program instructions, and when the computer program instructions are executed by the processor, the above-mentioned load-predictive-based heat pump unit temperature difference control method is implemented.

[0026] By adopting the above technical solution, a computer program is generated from the above-mentioned heat pump unit temperature difference control method based on load prediction and stored in the memory so that it can be loaded and executed by the processor. In this way, a terminal device can be made based on the memory and the processor for convenient use.

[0027] The beneficial effects of this invention are as follows:

[0028] This invention comprehensively captures the dynamic evolution and periodic patterns of unit parameters through functional principal component analysis and sine and cosine phase angle encoding. At the same time, it combines two-dimensional topological grid and distance decay retrieval mechanism to accurately screen similar operating conditions in the local search space, thereby strengthening the completeness of multi-dimensional feature expression from the bottom layer and effectively improving the clustering retrieval and matching efficiency in the context of massive historical data.

[0029] Furthermore, relying on the variance of the posterior probability distribution of the numerical output of the Gaussian process regression model, the uncertainty of the prediction is directly introduced into the underlying control loop, and a dual-line adaptive coordinated adjustment strategy of feedforward weight and feedback compensation is constructed. This effectively avoids the risk of lag and derailment caused by a single fixed parameter, and can dynamically weigh the correction force according to the sudden changes in the operating environment, significantly improving the response speed and energy-saving temperature control capability of heat pump equipment under complex operating conditions. Attached Figure Description

[0030] Figure 1 The flowchart shows the variable temperature difference control method for heat pump units based on load forecasting. Figure 2 To smooth the fitted curve for continuous time series data; Figure 3 This is a graph showing the regression prediction and uncertainty curve for a Gaussian process. Detailed Implementation

[0031] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0032] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0033] This invention discloses a variable temperature difference control method for heat pump units based on load forecasting, referring to... Figure 1 This includes steps S1-S4: S1. Acquire multi-source data and extract spatiotemporal features.

[0034] The processor, through a temperature sensor, relative humidity sensor, flow meter, and data acquisition module, reads compressor power, inlet and outlet water temperature, outdoor dry-bulb temperature, outdoor relative humidity, and actual heating and cooling load values ​​at a set sampling frequency via the Modbus communication protocol, forming an initial multivariate historical dataset. Subsequently, linear interpolation of missing values ​​is performed on the historical dataset to restore the time series to an equally spaced complete sequence. Further, the operating parameters and environmental parameters used as input features are normalized to standard deviation, i.e., converted into dimensionless variables based on their mean and standard deviation. Simultaneously, standard deviation normalization is performed on the actual heating and cooling load values, and the load mean and standard deviation of the training set are saved as the basis for subsequent anti-standardization and variance dimension recovery of the prediction results. The dimensionless variables obtained through standard deviation normalization are floating-point values, while the saved load mean and standard deviation both have kW dimensions.

[0035] Specifically, when a sensor loses connection during a sampling period or its collected data exceeds the aforementioned operating range, the system executes the following protection logic: for short-term missing or out-of-range sampling points, linear interpolation of the preceding and following valid sampling values ​​is used to complete the connection; for channels that are continuously out of contact, the system switches to the redundant backup sensor reading for that measurement point; if there is no redundant channel, the system maintains the current measurement point's output with the most recent valid value and marks the data as low confidence. During the execution of the protection logic, the data of the affected measurement points remains in a safe placeholder state without interrupting the subsequent process until the channel is restored and normal acquisition automatically resumes. In this way, different physical quantities are unified into dimensionless variables, missing and abnormal sampling are repaired, and a complete and equally spaced time sequence is provided for subsequent continuous curve fitting.

[0036] Standard deviation normalization is only one implementation method to achieve dimensional uniformity. Without departing from the concept of this invention, any dimensionless means that can map different physical quantities to comparable numerical scales, such as range normalization or stability normalization, can be applied to this step. Similarly, the Modbus communication protocol can be replaced with CAN bus or industrial Ethernet protocol without affecting the implementation of data access function.

[0037] Furthermore, low-dimensional and complete spatiotemporal features are extracted from the continuous changing trends of operating and environmental parameters. First, the sliding time window is set to the past 24 hours. With a sampling frequency of once every 15 minutes, a single time window contains 96 discrete time nodes; the time window length is adjustable to adapt to different operating condition switching rhythms. A B-spline basis function system of order 4 (cubic polynomial) with 20 internal nodes is selected, with the internal nodes distributed at equal intervals within the predetermined time window. The discrete sequences of each variable within the window are smoothed and fitted using the least squares method with a roughness penalty. The roughness penalty term is the square integral of the second derivative of the functional data curve. The smoothing parameter λ is adaptively determined using generalized cross-validation to achieve a balance between fidelity and smoothness in the fitted curve. (Refer to...) Figure 2 As shown, a continuous functional data curve is generated that can filter out some high-frequency measurement noise while retaining the main trend of change. Then, the standardized functional data curves of compressor power, outdoor dry-bulb temperature, and outdoor relative humidity are used as multivariate functional samples to construct a joint covariance matrix and perform eigenvalue decomposition. The eigenvectors corresponding to the three largest non-zero eigenvalues ​​are extracted as the basis functions for learning.

[0038] Next, using Simpson's integral rule or the trapezoidal numerical integral formula, the multivariate smoothed functional data curve of a certain day is projected onto the three learned basis functions mentioned above. The score of each principal component is the time integral of the product of the continuous functional data curve and the corresponding learned basis function, thus obtaining the scores of the first to third joint principal components. Taking a certain time window as an example, if the scores of the three principal components are 12.5, -3.2 and 1.1 respectively, they are sequentially concatenated to form a real-valued environmental feature vector with a dimension of 3 and components of 12.5, -3.2 and 1.1 respectively, which is used to characterize the comprehensive spatiotemporal variation characteristics of the operating parameters and environmental parameters within that time window.

[0039] In this way, the multidimensional time series originally consisting of 96 discrete sampling points is compressed into a continuous variation feature containing only three components. While significantly reducing the data dimensionality, it retains the dominant trend of unit operation and environmental evolution, providing a denoised and compact input for subsequent load forecasting.

[0040] In another embodiment, the B-spline basis functions used for smooth fitting can also be replaced by Fourier basis functions or wavelet basis functions to directly obtain the corresponding functional data curves. This substitution method has a more compact expression advantage compared to the B-spline system used in this embodiment when the time series data exhibits obvious periodic oscillation characteristics. Further optionally, the number of extracted learning basis functions can be fixed at three, or it can be adaptively determined according to the standard that the cumulative variance contribution rate reaches a set proportion.

[0041] To achieve the sub-objective of expressing periodic time patterns in the data access and feature construction module, discrete calendar time attributes are transformed into continuous features at the boundaries of the period. Specifically, the processor extracts three-dimensional integer values ​​from the timestamps of the historical dataset: month values ​​ranging from 1 to 12, day values ​​ranging from 1 to 7, and hour values ​​ranging from 0 to 23 within a 24-hour period. Then, angular frequency parameters are used to normalize these three values. The month angular frequency is determined by subtracting 1 from the month value, dividing by 12, and then multiplying by 2π; the week angular frequency is determined by subtracting 1 from the day value, dividing by 7, and then multiplying by 2π; and the hour angular frequency is determined by dividing the hour value by 24 and then multiplying by 2π, resulting in three normalized angular frequency values ​​distributed between 0 and 2π.

[0042] Next, the three normalized angular frequency values ​​are independently substituted into the sine and cosine functions to obtain six function results. These results are then concatenated in the order of their corresponding sine and cosine values ​​to form a periodic time feature vector of length 6. Taking the predicted time as August, Thursday, and 2 PM as an example, the sine and cosine values ​​corresponding to the month, day of the week, and hour are -0.500, -0.866, 0.434, -0.901, -0.500, and -0.866, respectively. This forms a 6-dimensional real number array whose elements are these six floating-point values, serving as a seamless periodic time feature vector for that moment.

[0043] In this way, the originally one-dimensional increasing discrete time variable is mapped to a smooth continuous coordinate on a two-dimensional unit circle, and the numerical abrupt changes at the periodic switching points such as the end of the year to the beginning of the year and the end of the week to the beginning of the week are eliminated, so that the time characteristics have a continuous expression at the periodic boundaries.

[0044] In another embodiment, the sine and cosine encoding used can also be replaced by radial basis-based periodic embedding or learnable periodic embedding. The alternative means can also achieve the technical effect of continuous expression of time variables at the periodic boundary. Furthermore, periodic dimensions other than months, days of the week, and hours, such as quarters or holiday attributes, can also be extended and added to the periodic time feature vector in the same way.

[0045] S2, Topological network clustering and nearest neighbor sample retrieval.

[0046] In one embodiment, a two-dimensional rectangular topological grid network with uniform length and width is established. In this embodiment, its size is set to 10 rows and 10 columns, with 100 competing layer neuron nodes. By randomly sampling within a continuous uniform distribution interval from -1 to 1, a 6-dimensional initial weight vector consistent with the dimension of the periodic time feature vector is assigned to each grid node. Subsequently, all periodic time feature vectors in the historical dataset are input sequentially, and the Euclidean distance between the input vector and the weight vector of each grid node is calculated iteratively. The node with the smallest Euclidean distance is selected as the winning neuron node for this iteration.

[0047] In the weight update stage, a learning rate with an initial amplitude of 0.1 that decays exponentially with the number of training steps to 0.01 is used. This learning rate is then combined with a Gaussian neighborhood topology update function centered on the winning neuron to adjust the weight space. This update function is introduced to numerically represent the relationship between the influence of the winning neuron on its surrounding nodes and the decay of influence as the grid topology distance increases, thus making adjacent nodes in the topology more likely to respond with similar time-period patterns. The Gaussian neighborhood topology update function is as follows:

[0048] in, To win the neuron node For the first The neighborhood influence value of each neuron node is a dimensionless real number. For nodes With the The Chebyshev distance of each node within a two-dimensional rectangular topological mesh network is taken from the integer planar coordinates of the mesh nodes. This represents the neighborhood radius, which decreases with the number of iterations. This function indicates that the smaller the topological distance between a node and the winning neuron, the closer its neighborhood influence value is to 1; conversely, the larger the distance, the closer the influence value is to 0. Furthermore, the neighborhood radius... As the iteration process decreases, the network undergoes large-scale coarse adjustments in the early stages of training and small-scale fine adjustments in the later stages of training. This corresponds to the physical process of the working mode cluster gradually converging from a coarse division to a clear boundary.

[0049] After obtaining the neighborhood action value, the weight vectors of the winning neuron node and its neighboring nodes with a Chebyshev distance of 1 are updated synchronously according to the following relationship:

[0050] in, and These are the first and second pages after the update, respectively. The node Each weight component The learning rate for the current iteration. For the input vector of the th Each component represents a component. This formula indicates that the node weights move along the input vector direction by an amount determined by both the neighborhood influence value and the learning rate. The larger the neighborhood influence value, the greater the magnitude by which the node moves closer to the input sample, thus making it easier for adjacent grids to respond together to similar time-period patterns. It should be noted that the specific forms in the above two formulas are only one implementation method for evaluating the relationship between neighborhood decay and weight movement in this embodiment. Other equivalent mathematical expressions based on the same physical relationship, such as using the Mexican hat function to characterize the neighborhood influence or adjusting the weights in a batch update manner, can also be applied to this step.

[0051] After 500 complete iterations of the full data loop, the entire network converges. All historical samples mapped to the same winning neuron are grouped into the same working pattern cluster, ultimately establishing 100 working pattern clusters with specific temporal behavioral boundaries. In this way, a large number of historical samples are organized into adjacent working pattern clusters according to time periodic patterns, laying the topological foundation for subsequent local retrieval based on working conditions.

[0052] Optionally, the 10x10 grid size of the two-dimensional rectangular topological grid used in this step is only a specific value for this embodiment. The grid size can be enlarged or reduced while maintaining the same length and width to balance the cluster granularity and clustering stability. At the same time, in addition to Euclidean distance, the similarity measure between the input vector and the node weight can also be Manhattan distance or cosine similarity. It can be understood that the above adjustments can all achieve the function of working condition clustering under the same competitive learning framework.

[0053] Furthermore, the processor is responsible for locating the working mode cluster to which the time to be predicted belongs during the prediction task, and delineating a targeted local search space accordingly. Specifically, the processor extracts the 6-dimensional periodic temporal feature vector of the time to be predicted and inputs it into the pre-trained self-organizing map network. It calculates the Euclidean distance between the vector and the weight vectors of each grid node, and selects the working mode cluster corresponding to the smallest distance as the best matching cluster. Then, the best matching cluster is jointly constructed with all historical samples in the eight neighboring clusters with a Chebyshev distance of 1 on the two-dimensional rectangular topological grid to form the local search space of the best matching cluster.

[0054] It should be noted that when the number of available historical samples within the best matching cluster and its neighboring clusters is insufficient to support the sample size required for subsequent nearest neighbor retrieval, the system executes the following degradation protection logic: automatically expanding the neighborhood range from a Chebyshev distance of 1 to a Chebyshev distance of 2, progressively incorporating samples from outermost neighboring clusters until the number of samples in the local search space meets the retrieval requirements. During the execution of the protection logic, the local search space expands in a controllable manner without degenerating into a global retrieval, until the sample sufficiency condition is met, at which point it reverts to the default eight-neighbor range. In this way, the retrieval originally targeting all historical samples is shrunk to a local subspace composed of similar conditions, significantly reducing the candidate size for subsequent nearest neighbor retrieval.

[0055] As another implementation, the adjacency range of the local search space can also be pre-configured by the user as a fixed larger or smaller neighborhood to directly obtain local subspaces with different coverage. The larger neighborhood has a smoother cross-condition coverage advantage compared to the eight adjacency range used in this embodiment when the operating conditions change frequently, while the smaller neighborhood has higher retrieval efficiency when the operating condition boundaries are clear.

[0056] Furthermore, a decay correction negatively correlated with topological distance is introduced for the retrieval distance of adjacent cluster samples. The processor first determines the affiliation of candidate historical samples. Specifically, if a candidate sample belongs to the best-matching cluster, its distance decay factor is set to 1; if a candidate sample belongs to an adjacent cluster, the Euclidean distance between the center coordinates of the adjacent cluster and the center coordinates of the best-matching cluster within the two-dimensional rectangular topological grid is calculated as the topological distance value. For example, the topological distance value between two clusters in adjacent rows and columns is 1.414. The decay factor is determined using a negative exponential function with the natural constant e as the base, where the exponent is the product of the topological distance value and a preset negative coefficient. The specific relationship is as follows:

[0057] in, The attenuation factor is a dimensionless real number greater than 0 and not greater than 1. The aforementioned topological distance value is a dimensionless real number, and the coefficient 0.5 is a preset negative coefficient value used in this embodiment to control the rate attenuation that increases with the increase of topological distance. In this way, the confidence reduction of cross-cluster samples relative to samples in the same cluster is evaluated in a continuously adjustable manner: the larger the topological distance, the smaller the attenuation factor, and the lower the similarity between the working condition to which the corresponding sample belongs and the working condition to be predicted.

[0058] For example, with a topological distance of 1.414 and a decay factor of 0.493, after obtaining the decay factor, the feature space Euclidean distance between the periodic time feature vector of the time to be predicted and each nearest historical sample of the adjacent cluster is calculated. This feature space Euclidean distance is then divided by the decay factor to obtain the weighted comprehensive distance. For example, when the feature space Euclidean distance of a sample in an adjacent cluster is 2.5, dividing it by the decay factor of 0.493 yields a weighted comprehensive distance of 5.07. When performing nearest neighbor query matching, the weighted comprehensive distance is used instead of the original distance to determine the proximity of the query match. In this embodiment, based on the KD tree spatial partitioning algorithm, the 15 nearest historical samples with the smallest weighted comprehensive distance are retrieved in the local search space. The number of nearest neighbors can be adjusted according to the data density.

[0059] In this way, the far-end samples across clusters are subjected to a distance penalty amplification corresponding to the topological distance. While ensuring the stability of retrieval across cluster boundaries, the retrieval results tend to retrieve historical samples with closer topological distance and smaller feature space distance, thereby improving the targeting of the nearest neighbor candidate set for the prediction case.

[0060] In another embodiment, the negative exponential decay function with the natural constant as the base can also be replaced by any other decay function that monotonically decreases with topological distance and has bounded values, such as inverse proportional decay or truncated linear decay. The replacement function can also achieve the technical effect of reducing the confidence of cross-cluster samples.

[0061] S3. Use Gaussian process regression to predict loads and calculate variance.

[0062] In one embodiment, the processor declares the prior mean function of the Gaussian process regression model as a zero function, that is, the prior expected value of the objective function is set to 0 before the introduction of observation data, and sets the exponential quadratic kernel function, that is, the composite covariance kernel function formed by the superposition of the squared exponential kernel function and the white noise kernel function, as the covariance kernel function to characterize the correlation between different environmental feature vectors.

[0063] In this embodiment, the length scale parameter of the squared exponential kernel is set to 1.5 to control the rate of covariance decay as the environmental feature distance increases; the signal variance parameter is set to 2 to characterize the magnitude of change in the target load function; and the white noise variance parameter is set to 0.1 to absorb the influence of sensor measurement errors and random disturbances on the observed load value. The initial values ​​of the above hyperparameters can be adjusted according to the data characteristics and optimized by maximizing the marginal log-likelihood function and using the finite-memory BFGS algorithm. Subsequently, 15 nearest-neighbor historical sample environmental feature vectors of length 3 are extracted and stacked row by row to construct a training input feature matrix with a dimension of 15 rows and 3 columns. The corresponding actual load values ​​are extracted simultaneously and standardized according to the load mean and load standard deviation of the training set to form an observation target vector of 15 rows and 1 column.

[0064] It should be noted that under normal operation, the kernel covariance matrix already includes white noise variance terms on the diagonal according to the composite covariance kernel function to ensure the numerical stability of matrix inversion. When the distribution of nearest neighbor samples is abnormally dense, causing the kernel covariance matrix to still be close to singular and the inversion to exhibit numerical ill-conditioned results, the system executes the following fault-tolerant protection logic: based on the already included default white noise variance terms, the value of the white noise variance parameter is temporarily increased for enhanced regularization processing to restore the matrix to a benign and invertible state; if it is still close to ill-conditioned after the increase, the white noise variance parameter is further increased and the current prediction is marked as low confidence. During the execution of the protection logic, the model inference maintains a computable stable state without interruption until the matrix benign condition is met, after which the default hyperparameters are restored. In this way, the model obtains a local covariance structure and tuned hyperparameters that match the operating conditions to be predicted, providing a stable numerical basis for subsequent prediction and inference.

[0065] Furthermore, the processor performs a cross-kernel operation on the 3D environmental feature vector at the time to be predicted and the training input feature matrix to obtain a 1x15 cross-covariance matrix. It then performs a kernel operation on the feature vector to be predicted and itself to obtain a 1x1 autocovariance matrix. Finally, it calculates the variance of the samples within the training input feature matrix to obtain a 15x15 full-rank symmetric kernel covariance matrix, and includes a white noise variance term in this kernel covariance matrix. The kernel covariance matrix is ​​then inverted, multiplied by the cross-covariance matrix, and then by the observed target vector. The resulting single-dimensional value is the standardized expected value. After performing destandardization, a predicted load value with kW dimensions is obtained, for example, 120.5kW. Subtracting the product of the cross-covariance matrix, the inverse of the kernel covariance matrix, and the transpose of the cross-covariance matrix from the self-covariance matrix yields the posterior variance in the standardized space, for example, 3.2. Multiplying this value by the square of the global training set load standard deviation saved above restores the posterior probability distribution variance with the dimension of the load square, i.e., kW², as shown below. Figure 3 As shown. The reason for using the standard deviation of the global training set load instead of the standard deviation of the local nearest neighbor samples for destandardization is that the standardization of the load values ​​has already been completed uniformly based on the global training set in step S1, and the destandardization of the predicted expected value also uses this global statistic; to maintain the consistency between the predicted value and its uncertainty dimension, the destandardization of the posterior variance uses the same benchmark. The predicted load value is used as the load basis for feedforward control, while the variance of the posterior probability distribution serves as the prediction uncertainty index in subsequent control and scheduling.

[0066] After this step, the system provides the expected load at the time to be predicted, along with the predicted dispersion value, so that the downstream control link can obtain both the feedforward quantity and the prediction reliability.

[0067] S4. Dynamically update parameters and output control commands.

[0068] In one embodiment, the processor substitutes the predicted load value into the thermodynamic energy balance relationship of the heat pump unit. That is, under the relationship that the predicted load value equals the product of the specific heat capacity of water, the density of water, the preset reference volumetric flow rate, and the feedforward reference temperature difference setpoint, the processor inversely solves for the feedforward reference temperature difference setpoint that satisfies the predicted load demand. Here, the specific heat capacity and density of water are known physical property constants, the preset reference volumetric flow rate is the system setpoint, and the feedforward reference temperature difference setpoint is the variable to be determined, which has the dimension of temperature difference (°C). In this way, the abstract load prediction is transformed into a feedforward control reference directly facing the supply and return water temperature difference, providing a clear setting target for composite control.

[0069] As another implementation method, the energy balance relationship of back-calculating the feedforward reference temperature difference from the predicted load can also be replaced by a pre-calibrated load-temperature difference lookup table or fitting curve to directly obtain the feedforward reference temperature difference setpoint.

[0070] Because the predicted load value output is uncertain due to variations in operating conditions, if feedforward and feedback are fused with fixed weights, the control command may deviate from the reasonable operating point when the prediction deviation is large. Therefore, the control parameters are scheduled online based on the variance of the posterior probability distribution. Specifically, the processor first calculates the ratio of the variance of the posterior probability distribution to the sum of the square of the predicted load value and a preset small positive number as the relative variance. The preset small positive number is 1kW² to avoid the denominator being 0, and the relative variance is a dimensionless quantity. Taking a posterior probability distribution variance of 200kW² and a predicted load value of 100kW as an example, the relative variance is 0.02. The processor pre-stores a dead zone comparison interval that reflects the prediction reliability. The upper and lower limit thresholds are determined as follows: Historical operational data is collected, and the ratio of the variance of the posterior probability distribution of the Gaussian process regression output at each historical moment to the sum of the square of the predicted load value and a preset small positive number is calculated to obtain a historical distribution sequence of relative variance. The 85th percentile of this historical distribution sequence is taken as the upper limit threshold, and the 15th percentile as the lower limit threshold. When historical data is lacking, the upper limit threshold is set to 0.03 by default, and the lower limit threshold is set to 0.01 by default, and these values ​​are continuously updated as the system operates. In this embodiment, the upper limit threshold determined by the above method is 0.03, and the lower limit threshold is 0.01.

[0071] In each real-time control cycle, such as every 250ms scan interval, the relative variance is compared with the upper and lower thresholds for judgment. Specifically, when the relative variance is greater than the upper threshold, it indicates that the predicted dispersion range is wide. The feedforward control weight corresponding to the predicted load is reduced, for example from 0.7 to 0.55, and the global proportional gain is increased by 0.1, for example from 1.5 to 1.6, and the integral increment gain is increased by 0.05, for example to 0.45. The increased gain factors are then limited to preset upper limits of no more than 3 and 1, respectively. When the relative variance is less than the lower limit judgment threshold, it indicates that the model prediction is relatively accurate. The feedforward control weight is increased, for example, from 0.7 to 0.85, and the global proportional gain is reduced by 0.1 and the integral increment gain is reduced by 0.05. The two reduced gain factors are then limited to preset lower limits greater than 0, which are not lower than 0.5 and 0.1, respectively. When the relative variance is between the lower limit judgment threshold and the upper limit judgment threshold, the PID gain parameter and the feedforward control weight remain unchanged, where the feedforward control weight is a real number between 0 and 1.

[0072] It should be noted that when the relative variance falls exactly within the insensitive dead zone between the lower and upper judgment thresholds, the system executes the following protection logic: the current state of the relevant PID gain parameters and feedforward control weights remains unchanged. During the execution of this protection logic, the actuator maintains the safe state of the current control parameters to avoid invalid actions and communication bus burden caused by frequent recalculation near the critical point, until the relative variance crosses any judgment threshold, at which point parameter updates resume. In this way, the control system adaptively allocates the feedforward dominance and feedback correction strength according to the prediction reliability; it enhances feedback control when prediction uncertainty is high and enhances feedforward following when prediction reliability is high.

[0073] To achieve the instruction output sub-objective of the composite variable temperature difference control module, the feedforward and feedback control quantities are merged and sent to the execution layer. The processor subtracts the obtained feedforward reference temperature difference setpoint from the actual supply and return water temperature difference collected by the real-time temperature sensor to obtain the temperature difference deviation, which is then input to the gain-scheduled PID controller that has completed parameter updates. The controller performs proportional, integral, and derivative operations on the temperature difference deviation to generate a PID feedback compensation term. Simultaneously, the feedforward control weight is multiplied by the feedforward reference temperature difference setpoint to obtain a weighted feedforward term. The weighted feedforward term and the feedback compensation term are then superimposed to generate the composite variable temperature difference control instruction. This composite variable temperature difference control instruction is converted into a target flow control quantity by the programmable logic controller according to a preset temperature difference-flow mapping relationship, and then converted into an analog voltage signal from 0 to 10V and transmitted to the circulating water pump frequency converter. The frequency converter adjusts the power supply frequency accordingly to drive the variable frequency asynchronous motor to change its speed to regulate the flow rate of the circulating water in the pipeline. At the same time, the electronic electric regulating valve adjusts its valve opening according to the control instruction to perform flow limiting regulation, thereby realizing the composite control of the system-side flow rate and the supply and return water temperature difference.

[0074] It should be noted that when the inverter communication link times out or the analog voltage output reaches the saturation limit, the system executes the following degraded protection logic: during the communication timeout period, the actuator continues to operate according to the most recent valid control command and triggers reconnection; when the output is saturated, the control quantity is limited to the safe output range and the integral accumulation is frozen to prevent integral saturation. During the execution of the protection logic, the unit maintains a safe and controllable operating state until communication is restored or the output exits the saturation zone, after which normal control operations are automatically resumed. In this way, predictive feedforward and temperature difference feedback are integrated into a single executable command and stably issued to the inverter motor and regulating valve, completing the real-time control of the system-side flow rate and supply / return water temperature difference.

[0075] Optionally, the 0 to 10V analog voltage signal transmitted to the execution layer can be replaced with a 4 to 20mA analog current signal or a digital control quantity transmitted via fieldbus.

[0076] This embodiment uses continuous operating records of a large commercial heat pump system over the past 12 months, with a sampling interval of 15 minutes, as the dataset. The records include outdoor dry-bulb temperature, relative humidity, compressor operating power, and actual heating and cooling loads. The first 10 months, totaling 28,800 records, are used as the training set, and the last 2 months' records are used as the test set. Load prediction performance is evaluated using root mean square error and mean absolute error, while control tracking performance is evaluated using mean steady-state response time and maximum temperature overshoot. The control setup includes a basic control group, an ablation group one with the grid topology distance attenuation factor removed, an ablation group two with the posterior variance scheduling control mechanism disabled, and an experimental group applying the complete scheme. The basic control group uses traditional fixed-parameter PID feedback control without introducing load prediction feedforward compensation; the PID gain parameter and feedforward control weight remain unchanged from their initial settings throughout the test.

[0077] Under this test dataset, the root mean square error (RMSE) of load prediction for the basic control group was 12.4 kW, the mean absolute error (MAE) was 8.6 kW, the mean steady-state response time was 320 s, and the maximum temperature overshoot was 2.5 ℃. In ablation group one, the nearest neighbor retrieval accuracy decreased after removing the topological distance decay mechanism, with a RMSE of 9.8 kW and an MAE of 7.1 kW. Ablation group two retained feature extraction and weighted nearest neighbor retrieval but disabled control parameter scheduling based on posterior variance; its prediction error was close to that of the complete scheme, but the MAE remained at 265 s and the maximum temperature overshoot was 1.8 ℃. The complete scheme had a RMSE of 5.8 kW, an MAE of 3.9 kW, a MAE of 185 s, and a maximum temperature overshoot of 0.4 ℃.

[0078] Based on the above results, under this test dataset, the complete scheme exhibits lower load prediction error and better control tracking performance compared to the baseline control group. The combination of periodic sine and cosine characteristics with functional principal component environmental characteristics enhances the expressive power of historically similar operating conditions; the topological distance decay mechanism improves the screening effect of local nearest neighbor samples; and the control scheduling mechanism based on the posterior variance of the Gaussian process can dynamically adjust the feedforward weights and PID feedback intensity according to prediction uncertainty. When prediction uncertainty is high, the system reduces the feedforward control weights and enhances feedback correction; when prediction reliability is high, the system increases the feedforward control weights and appropriately smooths the feedback adjustment, thus balancing response speed and operational stability to a certain extent.

[0079] This invention also discloses a load-predictive-based variable temperature difference control system for heat pump units, including a processor and a memory. The memory stores computer program instructions, which, when executed by the processor, implement a load-predictive-based variable temperature difference control method for heat pump units according to this invention.

[0080] The aforementioned load-prediction-based heat pump unit temperature difference control system also includes other components well known to those skilled in the art, such as communication buses and communication interfaces. Their settings and functions are known in the art and will not be described in detail here.

[0081] In the description of this specification, "multiple" or "several" means at least two, such as two, three or more, unless otherwise expressly and specifically defined.

Claims

1. A method for controlling the temperature difference of a heat pump unit based on load forecasting, characterized in that, include: S1. Obtain the operating parameters, environmental parameters, and actual load values ​​of the heat pump unit to form a historical dataset, as well as the current PID gain parameters and feedforward control weights; After dimensionless transformation of continuous time-series data within a predetermined time window, principal component scores are extracted using functional principal component analysis to construct environmental feature vectors. Discrete-time data are encoded using sine and cosine to generate periodic time feature vectors. Heat pump unit operating parameters include compressor power and inlet / outlet water temperature; environmental parameters include outdoor dry-bulb temperature and outdoor relative humidity; and actual load values ​​include actual cooling and heating load values. S2. A self-organizing mapping network is used to cluster the periodic time feature vectors into multiple operating mode clusters on a two-dimensional rectangular topological grid network. Local searches are then conducted using historical samples from each cluster and its neighboring clusters. In the spatial search, the periodic time feature vector of the time to be predicted is input into the network to determine the best matching cluster, and k nearest neighbor historical samples are retrieved in the local search space; S3, the topological distance is determined based on the coordinates of each cluster on the grid, and the distance between adjacent cluster samples is divided by a decay factor that is negatively correlated with the topological distance. The predicted load value and the variance of the posterior probability distribution are calculated using the environmental feature vector of the nearest neighbor historical samples and the actual load value through a Gaussian process regression model; the decay factor is determined by a negative exponential function with the natural constant e as the base, and its exponent term is the product of the topological distance value and a preset negative coefficient; S4. The predicted load value is used as the feedforward signal to calculate the reference temperature difference setpoint. The current PID gain parameter and feedforward control weight are adjusted according to the variance of the posterior probability distribution. The reference temperature difference setpoint is weighted by the feedforward control weight to obtain the weighted feedforward term. Then, it is superimposed with the updated PID feedback compensation term to generate an instruction and issue the instruction to drive the variable frequency motor to adjust the flow rate and perform current limiting control. Adjusting the current PID gain parameters and feedforward control weights based on the posterior probability distribution variance includes: calculating the ratio of the posterior probability distribution variance to the sum of the square of the predicted load value and a preset small positive number to generate a relative variance; When the relative variance is greater than the set upper limit judgment threshold, the feedforward control weight is reduced, and the global proportional gain coefficient and integral increment gain coefficient in the PID gain parameter are increased by a corresponding preset step size and limited to a preset upper limit; when the relative variance is less than the set lower limit judgment threshold, the feedforward control weight is increased, and the global proportional gain coefficient and integral increment gain coefficient are decreased by a corresponding preset step size and limited to a preset lower limit greater than 0; when the relative variance is between the lower limit judgment threshold and the upper limit judgment threshold, the PID gain parameter and the feedforward control weight remain unchanged.

2. The method for controlling the temperature difference of a heat pump unit based on load forecasting according to claim 1, characterized in that, The step of dimensionlessly converting continuous time-series data within a predetermined time window and extracting principal component scores through functional principal component analysis to construct an environmental feature vector includes: The continuous time-series data of the operating parameters and environmental parameters are standardized respectively; the standardized continuous time-series data are smoothed using B-spline basis functions to generate continuous functional data curves; the covariance matrix among all functional data curves is calculated and eigenvalue decomposition is performed, extracting the eigenvectors corresponding to the first three largest non-zero eigenvalues ​​as the learning basis functions; the principal component scores of each functional data curve on each basis function are calculated by integration, where the... The principal component score of the basis function is the sum of the functional data curve and the... The time integral of the product of the basis functions; the scores of the first three principal components are concatenated in order to form the environmental feature vector.

3. The method for controlling the temperature difference of a heat pump unit based on load forecasting according to claim 1, characterized in that, The step of generating periodic time feature vectors from discrete-time data using sine and cosine encoding includes: Extract the month, week, and hour values ​​from the historical dataset; normalize the month, week, and hour values ​​using angular frequency parameters to obtain three normalized angular frequency values; substitute the three angular frequency values ​​into the sine and cosine functions respectively to calculate six function result values, and concatenate the six function result values ​​to form the periodic time feature vector.

4. The method for controlling the temperature difference of a heat pump unit based on load forecasting according to claim 1, characterized in that, The method of clustering periodic time feature vectors into multiple working mode clusters on a two-dimensional rectangular topological grid network using a self-organizing map network includes: initializing the weight vectors of the two-dimensional rectangular topological grid nodes with consistent length and width; calculating the Euclidean distance between each periodic time feature vector and the weight vectors of all grid nodes, and selecting the grid node weight vector with the smallest Euclidean distance as the winning neuron node; and synchronously updating the grid weight vectors of the winning neuron node and its neighboring nodes with a Chebyshev distance of 1 according to the Gaussian neighborhood topological update function and a set learning rate. The Gaussian neighborhood topological update function is... ,in, The neighborhood effect value. The square of the Chebyshev distance between the two nodes within the grid. The neighborhood radius decreases with each iteration; after a set number of iterations, historical samples mapped to the same winning neuron node are incorporated into the same working mode cluster.

5. The method for controlling the temperature difference of a heat pump unit based on load forecasting according to claim 1, characterized in that, The division of the distance between adjacent cluster samples by a decay factor negatively correlated with the topological distance includes: The Euclidean distance between the center coordinates of the adjacent cluster and the center coordinates of the best-matching cluster within the two-dimensional rectangular topological grid network is calculated as the topological distance. The feature space Euclidean distance between the periodic time feature vector of the time to be predicted and each nearest historical sample of the adjacent cluster is calculated, and the feature space Euclidean distance is divided by the decay factor to obtain the weighted comprehensive distance. When performing nearest neighbor query matching, the weighted comprehensive distance is used instead of the original distance to determine the proximity of the query match.

6. The method for controlling the temperature difference of a heat pump unit based on load forecasting according to claim 1, characterized in that, The method of calculating the predicted load value and the variance of the posterior probability distribution using the environmental feature vectors of nearby historical samples and the actual load value through a Gaussian process regression model includes: Establish a prior distribution with the zero vector as the mean function, and define a kernel function containing an exponential quadratic kernel function as the covariance kernel function; extract... The environmental feature vectors of the nearest historical samples form a training input feature matrix. The corresponding actual load values ​​are extracted and standardized using the standard deviation of the training set load to form the observation target vector. Based on the covariance kernel function, the cross-covariance matrix, the autocovariance matrix, and the kernel covariance matrix are calculated. The standardized expected value of the prediction is obtained according to the regression formula and then de-standardized to the predicted load value. The standardized posterior variance is obtained according to the posterior covariance formula and multiplied by the square of the standard deviation of the training set load to obtain the posterior probability distribution variance with the dimension of the load square.

7. The method for controlling the temperature difference of a heat pump unit based on load forecasting according to claim 1, characterized in that, The acquired heat pump unit operating parameters, environmental parameters, and actual load values ​​constitute a historical dataset, including: The historical dataset is constructed by acquiring compressor power, inlet and outlet water temperature, outdoor dry bulb temperature, outdoor relative humidity, and actual cooling and heating load values ​​at a set sampling frequency. The operating parameters and environmental parameters, which are used as input features, are normalized according to their standard deviations. After normalizing the standard deviations of the actual load values, the training set load mean and load standard deviation are saved.

8. The method for controlling the temperature difference of a heat pump unit based on load forecasting according to claim 6, characterized in that, The covariance kernel function is a composite covariance kernel function formed by superimposing the exponential quadratic kernel function and the white noise kernel function, including: The hyperparameters of the composite covariance kernel function are optimized by maximizing the marginal log-likelihood function and using the finite-memory BFGS algorithm; a white noise variance term is added to the kernel covariance matrix to improve the stability of matrix inversion.

9. A variable temperature difference control system for a heat pump unit based on load forecasting, characterized in that, include: A processor and a memory, the memory storing computer program instructions that, when executed by the processor, implement a variable temperature difference control method for a heat pump unit based on load forecasting according to any one of claims 1-8.

Images