Correlation Study of Operational Data and System Performance of District Cooling System with Ice Storage

Hui Cao¹, Nan Li^1,2*, Jiajing Lin¹

¹School of Civil Engineering, Chongqing University, Chongqing, China.
²National Centre for International Research of Low-Carbon and Green Buildings, Ministry of Science & Technology, Chongqing University, Chongqing, China.

DOI: 10.4236/jpee.2024.123006                

Abstract

The district cooling system (DCS) with ice storage can reduce the peak electricity demand of the business district buildings it serves, improve system efficiency, and lower operational costs. This study utilizes a monitoring and control platform for DCS with ice storage to analyze historical parameter values related to system operation and executed operations. We assess the distribution of cooling loads among various devices within the DCS, identify operational characteristics of the system through correlation analysis and principal component analysis (PCA), and subsequently determine key parameters affecting changes in cooling loads. Accurate forecasting of cooling loads is crucial for determining optimal control strategies. The research process can be summarized briefly as follows: data preprocessing, parameter analysis, parameter selection, and validation of load forecasting performance. The study reveals that while individual devices in the system perform well, there is considerable room for improving overall system efficiency. Six principal components have been identified as input parameters for the cold load forecasting model, with each of these components having eigenvalues greater than 1 and contributing to an accumulated variance of 87.26%, and during the dimensionality reduction process, we obtained a confidence ellipse with a 95% confidence interval. Regarding cooling load forecasting, the Relative Absolute Error (RAE) value of the light gradient boosting machine (lightGBM) algorithm is 3.62%, Relative Root Mean Square Error (RRMSE) is 42.75%, and R-squared value (R²) is 92.96%, indicating superior forecasting performance compared to other commonly used cooling load forecasting algorithms. This research provides valuable insights and auxiliary guidance for data analysis and optimizing operations in practical engineering applications.

Keywords

DCS, Correlation Coefficient, PCA, Hourly Cooling Load, System Performance

1. Introduction

Since the Industrial Revolution, human civilization has mainly relied on and consumed fossil energy, allowing global industrialization and urbanization to develop rapidly, with civilization and prosperity growing by the day, while at the same time, the rapid development has also brought about environmental pollution, energy shortages, climate change and many other problems [1]. The Paris Agreement’s central aim is to strengthen the global response to the threat of climate change by keeping global temperature rise this century well below 2˚C above pre-industrial levels and to pursue efforts to limit the temperature increase even further to 1.5˚C (United Nations, 2015). Based on climate model simulations, experts suggest that we will reach 1.5˚C of global warming between 2030 and 2052, as show in Figure 1 [2]. Approximately one-third of the world’s total energy consumption is used in buildings, heating ventilation and air conditioning energy consumption in developed countries accounts for half of the energy use in buildings [3] [4], as show in Figure 2. At a global level, space cooling makes up only 6% of energy demand in buildings. This is because of low ownership of air-conditioners in many of the world’s warmest regions due to low incomes. As incomes increase over the Outlook in these regions, the number of households with space cooling rises rapidly, adding to the growth in electricity demand [4].

The electricity consumption for building operations accounts for nearly one quarters of the total electricity consumption in china, of which 60% can be attributed to heating, ventilation, and air conditioning equipment [5]. Cooling and heating loads are significant portions that can be shifted throughout time, as power grids are increasingly becoming stressed, commercial companies and institutions are hoping that incentives and electric cost plans to be encouraged on demand side load reduction [5] [6]. DCS techniques provide an alternative solution to enhance the energy efficiency of a central cooling system during part load conditions, it has been recognized as one of the effective methods to enhance the energy efficiency in buildings, especially for places where different tariff is adopted for top/peak and flat/valley section of energy used [7], At the same time, DCS in conjunction with the use of renewable energy sources will also help to reduce the rise in global temperature.

To ensure the efficiency of DCS, the control of the system needs to adapt to the aging of the system, require the least human intervention, and adapt to the changes of electricity prices in real time [8]. With the emergence of smart grid and the more penetration of unstable and intermittent renewable energy, DCS

Figure 1. A simple model for estimating the time when global warming is likely to reach 1.5˚C.

Figure 2. Energy demand by service in key countries and regions.

will be highly needed to balance the supply and demand sides of the power grid [9] [10]. Various factors influence the DCS loads such as building type, occupancy, meteorological conditions, equipment efficiency, thermostat setpoints, season and building controls [11] [12] [13]. Many artificial intelligence (AI) algorithms are used for data-driven load forecasting models. Such as artificial neural networks (ANN), random forests (RF), generalized regression neural network (GRNN), eXtreme gradient boosting (XGBoost), suppport vector machine regression (SVMR), long and short term memory(LSTM), deep neural networks(DNN) and light gradient boosting machine (LightGBM) algorithm, etc. [14]-[22]. For AI load forecasting models, the quality of the training samples largely determines the forecasting accuracy. There are always some algorithms applies better, that can produce better results for lower dimensional data inputs, with greater interpretability and smaller computational resources [23]. Based on a large number of literature studies, we chose Algorithm LSTM to predict the hourly cooling loads and to compare the study with some other load forecasting algorithms that are commonly used.

In the evaluation of the system, the decision variables of load forecasting should be determined first. Given the new trends of data collection and data analytics, data-driven solution can more accurately solve the optimization problem for the DCS. Our research evaluated the system’s current operating status and contributed to the subsequent improvement of performance and to the full realization of the savings potential of the system. Also, the field measurements are collected and will be used to develop a mathematical model for the DCS and to evaluate the proposed optimization control strategy under a demand response program. Laying the groundwork for the optimize control research process in the future which runs in real time to properly update the ice storage charge and adapt to weather forecasts and tariff changes.

2. Methodology

2.1. Information of the system

The operational structural configurations for the primary and secondary sides of the system are depicted in Figure 3. The system is designed with a specified cooling total load of 187,517 kW and a heating total load of 85,465 kW. Its structural arrangement features a conventional setup where refrigeration units are positioned upstream, followed by a series arrangement of open-type discharge ice systems with connected ice tanks. The designed capacity of these ice tanks is 330,785 kWh.

The primary pumps mainly address the resistance encountered from equipment such as refrigeration units and filters. In contrast, the secondary pumps are designed to counteract the resistance from the plate heat exchangers at the terminal user and the associated distribution network. Given the configuration of the system, where refrigeration units are positioned upstream, followed by a series-connected open-type discharge ice system with ice tanks, when the units collaborate with the ice storage tanks for cooling supply, the network’s return water undergoes primary cooling through plate heat exchangers corresponding to both the heat pump unit and the duplex unit. Subsequently, it undergoes secondary cooling and heat exchange via the discharge ice supply plate heat exchanger, ensuring the maintenance of a relatively stable supply water temperature.

The system mainly has 18 heat pump units, including 10 heat pump units and 8 duplex chiller units. The system is also equipped with 4 river water intake pumps, 5 primary chilled water pumps, 8 external network circulation pumps, 8 sets of cooling plate heat exchange, 6 sets of ice discharge plate heat exchange, and 8 ice discharge chilled water pumps. There is also a comprehensive energy consumption monitoring and control platform, which is the basis for efficient energy saving and safe operation and maintenance of the energy station. The platform can monitor and record the operating status parameters of DCS online 24 hours.

2.2. Data Preprocessing

The cold load calculation formula is

$Load=Q\rho C(T_{in}-T_{out})/3600$ (1)

Figure 3. Sketch of the operating mode structure on the primary and secondary side of the system.

where Q is the chilled water flow (m³/h); $\rho$ is the chilled water density (kg/m³); C is the specific heat capacity (kJ/kg ˚C); $T_{out}$ is the supply temperature (˚C); $T_{in}$ is the return temperature (˚C).

The preprocessing of load data primarily is data smoothing techniques: Determine whether the data in the load sequence for day d at time t exhibits any anomalies. Here, $t=0,1,2,\mathrm{...},23$ , $d=1,2,3,\mathrm{...},n$ . The specific steps for smoothing treatment are as follows:

(1) Calculate the mean $E(t)$ and standard deviation $\sigma (t)$ of the cooling load at time t over n days.

$E(t)=\frac{1}{n}{\displaystyle {\sum }_{d=1}^{n}x(d,t)}$ (2)

$\sigma (t)=\sqrt{\frac{1}{n}{\displaystyle {\sum }_{d=1}^n{\left[x\left(d,t\right)-E\left(t\right)\right]}^2}}$ (3)

(2) Calculate the deviation rate $\eta (d,t)$ of the cooling load at time t on day d and the maximum deviation rate ${\eta }_{\max }(d,t)$ at time t over n days.

$\eta (d,t)=\frac{\left|x(d,t)-E(t)\right|}{\sigma (t)}$ (4)

(3) Identify potential outliers. Set the maximum allowable deviation rate as C. If ${\eta }_{\max }(d,t)>C$ , it is deemed that an outlier exists at time t. Subsequently, verify if $\eta (d,t)>C$ holds true, thereby completing the assessment of the deviation rates at time t across all historical days to pinpoint all outliers. In this study, the empirical value set for C is 3.

(4) Correct the identified outliers. Replace the outlier $x(d,t)$ with the mean $\bar{x}(d,t)$ of time t from one day prior and one day subsequent (or multiple days) to the outlier.

$\bar{x}(d,t)=\frac{1}{2}\left[x(d-1(n),t+x(d+1(n),t)\right]$ (5)

Repeat the process from steps 2 to 4 to identify all potential outliers, thereby completing the data smoothing procedure. It is also worth noting that in this study, data preprocessing is not solely limited to the cooling load data; it is also applicable to the data concerning factors affecting the cooling load. When necessary, preprocessing techniques can be employed for these datasets as well.

2.3. Correlation Study

Given our data collection involves a preprocessing step, this study opts for the Spearman rank correlation coefficient. The Spearman rank correlation coefficient is a non-parametric measure, denoted as $r_{\text{s}}$ , which is independent of the distribution [24]. It quantifies the linear correlation between ordinal variables. Assuming that the original sample data $(x_i,y_i)$ is arranged in ascending order, denoting $(x_i{}^{\text{'}},y_i{}^{\text{'}})$ as the positions of the sorted data from $(x_i,y_i)$ , $(x_i{}^{\text{'}},y_i{}^{\text{'}})$ is referred to as the ranks of $(x_i,y_i)$ . $d_i=x_i{}^{\text{'}}-y_i{}^{\text{'}}$ represents the differences in the ranks of $(x_i,y_i)$ .

If there are no tied ranks, $r_{\text{s}}$ is given by:

$r_{\text{s}}=\frac{n{\displaystyle \sum x_i^{\prime }}{y}_i^{\prime }-{\displaystyle \sum x_i^{\prime }}{\displaystyle \sum y_i^{\prime }}}{\sqrt{n{\displaystyle \sum x_i^{\prime 2}}-{\left({\displaystyle \sum x_i^{\prime }}\right)}^2}\sqrt{n{\displaystyle \sum y_i^{\prime 2}}-{\left({\displaystyle \sum y_i^{\prime }}\right)}^2}}=1-\frac{6{\displaystyle \underset{i=1}{\overset{n}{\sum }}{d}_i^2}}{n\left(n^2-1\right)}$ (6)

The Spearman rank correlation coefficient is often regarded as the Pearson correlation coefficient between the variables after ranking. If tied ranks are present, it becomes necessary to compute the Pearson correlation coefficient between the ranks. In such a case, $r_{\text{s}}$ is given by:

$r_{\text{s}}=\frac{{\displaystyle \underset{i}{\sum }\left(x_i^{\prime }-\overline{x^{\text{'}}}\right)}\left(y_i^{\prime }-\overline{y^{\text{'}}}\right)}{\sqrt{{\displaystyle \underset{i}{\sum }{\left(x_i^{\prime }-\overline{x^{\text{'}}}\right)}^2}{\displaystyle \underset{i}{\sum }{\left(y_i^{\prime }-\overline{y^{\text{'}}}\right)}^2}}}$ (7)

For the aforementioned correlation coefficients, n ideally should not be less than 10. When n is around 10, a correlation coefficient of at least approximately 0.7 indicates a close relationship. For $n\ge \text{10}$ , a correlation coefficient greater than 0.3 suggests that the two variables have met the threshold for a close relationship. In this study, n corresponds to the volume of data over a two-month period, with the order of magnitude being 10³. Therefore, concerning the degree of correlation between parameters studied in this paper, our findings are summarized in Table 1.

2.4. Principal Component Study

Regarding the input parameters for the cooling load, PCA is employed to perform linearly independent transformations and dimensionality reduction on the input parameters, ultimately determining the input parameters for the cooling load forecasting model.

It should be clarified that principal components can be derived either from the covariance matrix $\Sigma$ (where variables are not standardized) or from the correlation matrix R (where variables are standardized). In this study, the option of standardized variables is chosen.

Suppose there are n samples, with each sample observing m variables. This yields the original sample data matrix:

$X={\left(\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1m}\\ x_{21}& x_{22}& \cdots & x_{2m}\\ ⋮& ⋮& & ⋮\\ x_{n1}& x_{n2}& \cdots & x_{nm}\end{array}\right)}_{n\times m}$ (8)

x is an m-dimensional vector composed of m random variables $x_1,x_2,\mathrm{...},x_m$ , denoted as $x={(x_1,x_2,\cdot \cdot \cdot ,x_m)}^{{\rm T}}$ . Its mean vector is represented as $\mu$ , the covariance is $\Sigma$ , and the linear combination is:

$\{\begin{array}{l}{F}_1=a_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1m}{x}_m\hfill \\ F_2=a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2m}{x}_m\hfill \\ \cdots \hfill \\ F_m=a_{m1}{x}_1+a_{m2}{x}_2+\cdots +a_{mm}{x}_m\hfill \end{array}$ (9)

Table 1. The values of $\left|r_{\text{s}}\right|$ and their corresponding degrees of correlation.

To standardize the columns of X, we have:

$x_{ij}{}^{\ast }=\frac{x_{ij}-{\bar{x}}_j}{{\sigma }_j}$ (10)

where $i=1,2,\cdots ,n$ , $j=1,2,\cdots ,m$ , ${\bar{x}}_j=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_{ij}}$ , ${\sigma }_j{}^2=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(x_{ij}-{\bar{x}}_j\right)}^2}$ .

Clearly:

${\bar{x}}_j{}^{\ast }=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(x_{ij}{}^{\ast }\right)}=0$ (11)

${\sigma }^{\ast }{{}_j}^2=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(x_{ij}-{\bar{x}}_j\right)}^2}=1$ (12)

$\Vert x_j{}^{\ast }\Vert =\sqrt{x_j{{}^{\ast }}^{{\rm T}}{x}_j{}^{\ast }}=\sqrt{n}$ (13)

where $x_j{}^{\ast }={(x_{1j}{}^{\ast },x_{2j}{}^{\ast },\cdots ,x_{nj}{}^{\ast })}^{{\rm T}}$ .

At this point, for the sake of convenience in notation, the standardization of the variable $x_j{}^{\ast }$ is represented without the subscript $\ast$ . Consequently, in the sample covariance matrix S, the variances of variables $x_j$ and $x_k$ is:

${\sigma }_{jk}=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(x_{ij}-{\bar{x}}_j\right)}\left(x_{ik}-{\bar{x}}_k\right)=\frac{1}{n}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{x}_{ij}}{x}_{ik}=\frac{1}{n}{x}_j{}^{{\rm T}}{x}_k$ (14)

simultaneously, the correlation coefficient between the variables $x_j$ and $x_k$ is:

$r_{jk}=\frac{{\sigma }_{jk}}{\sqrt{{\sigma }_{jj}}\sqrt{{\sigma }_{kk}}}={\sigma }_{jk}=\frac{1}{n}{x}_j{}^{{\rm T}}{x}_k$ (15)

From this, it can be deduced that after standardization, the sample covariance matrix S and the correlation matrix R of the variables are identical. Moreover:

$R=S=\frac{1}{n}{X}^{{\rm T}}X$ (16)

where X represents the data after standardization, and this property can be succinctly denoted as

$\mathrm{cov}(X)={(V^{1/2})}^{-1}\Sigma {(V^{1/2})}^{-1}\triangleq R$ (17)

here $V^{1/2}=diag(\sqrt{{\sigma }_{11}},\sqrt{{\sigma }_{22}},\cdots ,\sqrt{{\sigma }_{mm}})$ . Then the problem of PCA transforms into the task of determining the principal components starting from the correlation matrix R.

2.5. Cooling Load Forecasting Algorithm

The XGBoost algorithm has been recognized as one of the most accurate methods for cooling load forecasting. Prior to the introduction of LightGBM, XGBoost was the most renowned gradient boosting decision tree (GBDT) tool available [25]. However, when dealing with large datasets with multiple features, training with XGBoost can be time-consuming and memory-intensive. LightGBM offers optimizations over traditional GBDT algorithms, addressing the aforementioned limitations of XGBoost [14]. Without compromising accuracy, LightGBM accelerates the training speed of GBDT models. Hence, this study adopts the LightGBM algorithm. LightGBM is a Gradient Boosting Machine (GBM) algorithm utilized for both classification and regression tasks. It employs an ensemble learning method based on trees and leverages gradient boosting to combine multiple weak learners, typically decision trees, into a robust model.

3. Results and Discussion

3.1. The Correlation between System Data and System Performance

3.1.1. Cooling Load Data with System Performance

Define After data preprocessing, the hourly cooling load from June 1st to August 31st during the cooling season of 2022 was determined. Within this dataset, we selected the hourly cooling load for a typical day. Subsequently, a comparative analysis was conducted with the hourly cooling load of the design day, as illustrated in Figure 4.

Figure 4. Hourly cooling loads and ratio for design day and typical day.

Based on the two figures above, it is evident that the typical daily load of the system closely aligns with the design day load, exhibiting a trend of increase followed by a decline. This pattern aligns with the work and life routines of individuals primarily associated with office buildings served by this system. On a daily basis, the peak load predominantly occurs between 10:00 AM and 5:00 PM, with a significant drop in load demand after 7:00 PM. Between 11:00 PM and 6:00 AM, there is minimal demand for load. The design day load is substantially higher than the typical daily load, with an average load rate difference of 0.4. This discrepancy arises because the district is an emerging economic zone where the DCS may potentially cater to a larger number of terminal-users in the future. This study, grounded in historical and current cooling load data of the DCS, takes into account the anticipated increased cooling demand due to regional growth. By integrating ice storage technology, it offers crucial parameter data and feasibility analysis for optimizing the operation of the DCS.

3.1.2. Refrigeration Unit Data with System Performance

Through the energy management platform, hourly data for the evaporator and condenser sides, including inlet and outlet flow rates, inlet and outlet water temperatures, unit power, current ratios, etc., of a typical cooling month (August 2022) for the refrigeration units were obtained. This study collects data from the typical month, focusing on the duplex chiller units during the nighttime (00:00- 7:00) for ice charge and the heat pump units during the day (8:00-20:00) for cooling.

The hourly values and distribution of the COP for the unit throughout this month can be determined, as illustrated in Figure 5.

It can be observed that within this typical month, the hourly COP values for the duplex chiller units predominantly range between 3.5 and 4.4, with a monthly average COP value of 3.96. The rated COP value for the duplex chiller

Figure 5. Visualization of hourly COP statistics for a typical month.

units is 4.37. The hourly and monthly average COP values for the duplex chiller units are quite close to their rated values, indicating their satisfactory operational performance. This is attributed to the fact that during nighttime ice charge, the duplex chiller units often operate at maximum capacity. Apart from start-up and shut-down phases, the load ratio of the duplex chiller units remains relatively high most of the time.

Similarly, for the heat pump units, the actual hourly COP values predominantly range between 4.3 and 5.4, with a monthly average COP value of 4.73. However, the rated COP value for the heat pump units under standard conditions is 5.40. There’s a considerable deviation between the hourly and monthly average COP values and their rated values. One contributing factor could be that during the summer, the sand content in the river water exceeds the filtration capacity of the intake system, leading to reduced water supply on the condensing side of the units. Additionally, it’s evident that the current DCS experiences relatively low-end user loads. During the daytime, cooling predominantly relies on discharging ice. The heat pump units are used to complement the loads that exceed the capacity of the ice storage tanks, resulting in a lower load ratio for these units. The primary reason behind this phenomenon is the poor self-regulation capability of the DCS. There’s a frequent need to manually start or stop the heat pump units in response to fluctuations in the system’s terminal demand, causing the units to operate at prolonged low load ratios. To enhance the operational efficiency of the heat pump units, improvements in the control system are essential. It’s crucial to allocate the load responsibilities of the system equipment reasonably and plan the loads of the refrigeration units based on load forecasts. Optimizing the start-stop operations of the units will ensure their operation at higher load ratios.

3.1.3. Ice Tank Data with System Performance

Figure 6 depicts the cooling details of the DCS within this typical month for the district. The ice discharge for cooling is divided into two sections based on tariff: valley/flat section and top/peak section. Within this typical month, the cooling load for most days is predominantly managed by ice discharge cooling. However, there are only 9 days where the ice cooling ratio is at 1, and these days predominantly fall on weekends. On these specific days, the cooling units within the system remain inactive.

Furthermore, the cumulative ice cooling volume has never reached the designed storage capacity of the ice storage tanks. Regardless of whether the terminal-user cooling demand surpasses the designed storage capacity of the ice storage tanks, except for those 9 days, the refrigeration units are activated. This implies that the potential of the ice storage tanks is not fully utilized. This phenomenon often occurs during top/peak tariff periods when the system activates the refrigeration units to meet the terminal-user demands. Initiating the refrigeration units during top/peak tariff periods is typically uneconomical. The average daily ice cooling volume for this typical month stands at 212,908.87 kWh,

Figure 6. Details of cooling in a typical day and month for DCS.

representing only 64.36% of the ice storage tank’s designed capacity. Given this scenario, the average ice cooling ratio for the ice storage tanks is 0.85, indicating excellent defrosting performance. However, the potential of the DCS to optimize peak shaving and load balancing remains underutilized, underscoring the need for an hourly load optimization for the system’s components.

3.1.4. The current Overall Performance of the System

Figure 7 illustrates the daily fluctuations of the Energy Efficiency Ratio (EER) and the Cost Per Unit of Cooling Capacity (CPOC) during the typical summer month. EER represents the ratio of the system’s cooling capacity to its total power consumption, while CPOC denotes the ratio of the system’s operating cost to its cooling capacity. The EER predominantly falls within the range of 2.23 to 3.11, with an average of 2.79, slightly below the empirical average (3.02) derived from extensive studies and statistical analyses of relevant engineering data in the country. Meanwhile, the CPOC is mainly concentrated between 0.13 and 0.19 yuan/kWh, with an average of 0.16 yuan/kWh, also lower than the empirical average (0.2 yuan/kWh) based on certain engineering benchmarks. This indicates that the system’s operational cost-effectiveness is relatively favorable.

Compared to the current control strategies, if this study can further enhance the energy-saving rate and cost-saving rate of the system while fully leveraging the advantages of peak-shifting and valley-filling through the ice storage system, it implies that by improving the system’s EER, the CPOC can be further reduced. Therefore, there is considerable room for improvement in both the energy efficiency and economic aspects of the system.

Figure 7. CPOC and EER for the DCS during the typical summer month.

3.2. Preliminary Determination of input Parameters

Regarding the DCS’s cooling load and its forecasting models, the influencing factors on the cooling load primarily encompass meteorological factors and intrinsic factors of the DCS itself.

For meteorological influencing factors, this study primarily investigates the factors that relatively have a significant impact on the cooling load: dry bulb temperature T, relative humidity φ, solar radiation intensity R, precipitation amount W, atmospheric pressure P, and atmospheric wind speed V.

For the intrinsic factors of the DCS itself, factors such as indoor set parameters, internal disturbance parameters, occupancy rate of individuals indoors, and usage rate of indoor equipment cannot be overlooked. Given the broad supply range of the DCS and the multitude of terminal-users, obtaining comprehensive and precise data is impractical. However, the factors affecting the cooling load by the DCS itself are highly correlated with time. Therefore, we have comprehensively introduced the time point count factor t and the historical cooling load factor $L\text{-}n,n\in (1,2,\cdots ,24)$ .

In this study, data on cooling load and factors influencing the cooling load from July 1st to August 31st, 2022, were collected. After preprocessing this data, the Spearman rank correlation coefficient was utilized to perform a 24h correlation analysis on these load influencing factors, assessing the correlation level between each factor and the predicted cooling load L at the specific time point.

To avoid redundant analysis, this study focuses solely on the correlation analysis between T and L, as well as the selection of input parameters for the cooling load forecasting model related to T, as a case for investigation. Figure 8 depicts the Spearman rank correlation coefficient between $T\text{-}n,n\in (0,1,2,\cdots ,24)$ and L, while Figure 9 illustrates the correlation trend between $T\text{-}n,n\in (0,1,2,\cdots ,24)$ and L over a 24-hour period.

The correlation changes between non-adjacent factors exhibit corresponding characteristics. As indicated by Figure 9, with increasing time intervals, the variation trend of the correlation coefficient between T and L resembles a sine curve. Beyond a time interval of 16, the correlation reverts to positive, peaking at a correlation coefficient value of 0.57 when the time interval reaches 22.

Using an absolute value of the correlation coefficient between T and L set at 0.4 as the threshold, the values for n within set $T\text{-}n,n\in (0,1,2,\cdots ,24)$ should be 0, 1, 9, 10, 11, 19, 20, 21, 22, 23, and 24. From this, we deduce the input parameters for the cooling load forecasting model associated with T.

Figure 8. Plot of calculated Spearman rank correlation coefficients for T and L.

Figure 9. Trend curve of correlation between T-n and L over 24 hours.

Subsequently, we analyze the correlation between relative humidity φ-n, solar radiation intensity R-n, hourly precipitation W-n, atmospheric pressure P-n, atmospheric wind speed V-n, time point count t-n, and the historical cooling load L-n with the predicted cooling load L. This will allow us to preliminarily determine the input parameters for the cooling load forecasting model. The initially identified input parameters for the cooling load forecasting model total 44 items, as shown in Table 2.

3.3. Dimensionality Reduction of input Parameters

For machine learning-based cooling load forecasting model, having too many input parameters can increase computational complexity. In this study, after carefully considering various influencing factors, our focus was on selecting the parameters that have the most significant impact on the cooling load forecasting model to enhance both its accuracy and generalization capability. Among the 44 preliminary cooling load input parameters identified in Table 2, PCA was employed to perform a linearly independent transformation and dimensionality reduction of the inputs, ultimately determining the input parameters for the cooling load forecasting model.

Sequentially compute the mean and standard deviation for each input parameter. Then standardize the data to derive the correlation matrix R for the 44 input parameters. Based on R, determine its eigenvalues, resulting in a scree plot that illustrates the relationship between each parameter (principal component) and its corresponding eigenvalue. In the scree plot, the x-axis indicates the order of the eigenvalues (or principal components), while the y-axis represents the values of the eigenvalues, as illustrated in Figure 10. Alongside the eigenvalues, measures such as variance contribution rates and cumulative variance contribution rates are also obtained, detailed in Table 3.

Table 2. Preliminary determination of input parameters for the cooling load forecasting model.

Figure 10. Scree plot for PCA of each parameter.

Table 3. The eigenvalues and variance contribution rates corresponding to the correlation matrix.

For all parameters (principal components), one can obtain an eigenvector corresponding to each eigenvalue. The parameter loading matrix is a matrix composed of eigenvectors corresponding to the eigenvalues. In simpler terms, the loading vector is synonymous with the eigenvector and simultaneously represents the direction of the principal component. Mathematically speaking, the score of a principal component can be understood as the projection of a point onto the loading vector, representing the weights of each parameter point with respect to that principal component. The principal component is essentially a combination of these weights.

Firstly, corresponding to the 44 eigenvalues presented in Table 3, a biplot is generated as illustrated in Figure 11. The distances between points representing parameters (principal components) in the graph approximately signify the similarity between those parameters. The cosine value of the angle between the loading vectors approximately indicates the correlation between variables of the parameters. The projection of a point onto a vector approximately represents the interaction, or weight, between the parameter (principal component) and its variables. We obtain a 95% confidence ellipse, implying that with a probability as high as 95%, this is a relatively small three-dimensional ellipse for variables PC1(50.6%), PC2(18.1%), and PC3(8.8%). This indicates that the values of the parameters we have chosen closely approximate the true values, suggesting a high precision in estimating the sample parameters from the population parameters.

Secondly, as observed back from Figure 10, the scree plot starts to plateau after

Figure 11. Biplot for PCA.

the sixth point, suggesting it can be theoretically disregarded. The criteria for selecting principal components are based on: identifying inflection points of eigenvalues, considering the magnitude of these eigenvalues, and evaluating the contribution rate of variance. It is evident that the sixth point in the plot signifies the inflection point of the eigenvalues. Referring to Table 3, this sixth point corresponds to T-19. Furthermore, if we focus solely on parameters with eigenvalues exceeding 1, it’s evident that all eigenvalues are greater than 1 before the inflection point; post the inflection point, the eigenvalues corresponding to the eigenroots are all less than 1. Considering the cumulative variance contribution rate, when this rate surpasses 80%, it implies that the principal components essentially encapsulate all the information inherent in the parameters. The cumulative variance contribution rate corresponding to the sixth point is 87.26%. Consequently, in descending order of eigenvalues, this study selects the first six principal components from Table 3 as input parameters for the cooling load forecasting model.

Lastly, let’s denote the eigenvalues of the first six principal components as ${\lambda }_{\text{1}}\ge {\lambda }_2\ge {\lambda }_{\text{3}}\ge {\lambda }_{\text{4}}\ge {\lambda }_{\text{5}}\ge {\lambda }_{\text{6}}$ corresponding to eigenvectors $e_{\text{1}},e_2,e_{\text{3}},e_{\text{4}},e_{\text{5}},e_{\text{6}}$ , respectively. With this, dimensionality reduction can be carried out. Additionally, the first six principal components are denoted as $F_{\text{1}},F_2,F_{\text{3}},F_{\text{4}},F_{\text{5}},F_{\text{6}}$ . Listing the 44 parameters from Table 3 sequentially as $x_i(i=1,2,\cdots ,44)$ , we conclusively determine the relationship between these 6 principal components and the selected 44 parameters. Essentially, they encapsulate the entirety of the information possessed by all parameters. These 6 principal components will serve as the input parameters for the cooling load prediction model, as illustrated by Equation (18).

$F_n=\left(x_1{x}_2\cdots x_{44}\right)e_n,n\in \left(1,2,\cdots ,6\right)$ (18)

3.4. The Performance of the Cooling Load Forecasting Model

Using the LightGBM algorithm and based on operational data from June 1st to August 31st, 2022, as the sample set, the model was trained to forecast the hourly cooling load for the entire month of August 2022. Figure 12 presents the hourly cooling load forecasting results over 31 days. For the majority of the time periods, the forecasted values of the hourly cooling load align closely with the actual measured values. Given that the end-users of the DCS are predominantly office buildings, there is a distinct difference in curve variations between weekdays and weekends.

We compared the forecasting performance of the LightGBM algorithm for cooling load with other algorithms mentioned in current relevant literature, specifically the GRNN, SVMR, and XGBoost. Using data from the entire month, we evaluated their forecasting capabilities using three commonly used metrics, the metrics are RAE(Relative Absolute Error), RRMSE(Relative Root Mean Square Error), and R-squared value (R²), they are defined as [8] [26]:

Figure 12. Visual comparison of predicted and measured hourly cooling loads.

$\text{RAE}=\frac{\frac{1}{n}\underset{m=1}{\overset{n}{{\displaystyle \sum }}}\left|y_m-y_{pre,m}\right|}{\frac{1}{n}\underset{m=1}{\overset{n}{{\displaystyle \sum }}}{y}_{pre,m}}\times 100\%$ (19)

$\text{RRMSE}=100\times \frac{\sqrt{\frac{{\displaystyle {\sum }_{m=1}^n{(y_m-y_{pre,m})}^2}}{n}}}{\frac{1}{n}\underset{m=1}{\overset{n}{{\displaystyle \sum }}}{y}_{pre,m}}\times 100\%$ (20)

${\text{R}}^2=100\times \left(1-\frac{\frac{{\displaystyle {\sum }_{m=1}^n{(y_m-y_{pre,m})}^2}}{n}}{\frac{1}{n}\underset{m=1}{\overset{n}{{\displaystyle \sum }}}{(y_{pre,m}-\frac{1}{n}\underset{m=1}{\overset{n}{{\displaystyle \sum }}}{y}_{pre,m})}^2}\right)\times 100\%$ (21)

where $y_m$ is the actual value, $y_{pre,m}$ is the predicted value, and n is the number of sample, $m\in [1,n]$ .

RAE and RRMSE are used to measure the deviation between forecasted values and actual measurements, with lower values indicating better performance. R² quantifies the degree of correlation between forecasted values and measured values, with a higher value suggesting better correlation. Figure 13 illustrates the performance of various load forecasting algorithms. It is evident that the LightGBM algorithm outperforms other forecasting methods. Specifically, its RAE is 3.62%, RRMSE is 42.75%, and R² is 92.96%. Additionally, it is observed that the SVM algorithm exhibits the least accuracy in the cooling load forecasting process.

Figure 13. Forecasting performance of various algorithms.

4. Conclusions

Within a typical month, the hourly COP values for the duplex chiller units predominantly range between 3.5 and 4.4, with a monthly average COP value of 3.96. Both the hourly and monthly average COP values for the duplex chiller units closely align with their rated values. For the heat pump unit, the actual hourly COP values primarily lie between 4.3 and 5.4, resulting in a monthly average COP value of 4.73, indicating a significant deviation from its rated value. Throughout the typical month, the system’s daily EER is primarily concentrated within the range of 2.23 to 3.11, with an average value of 2.79. Similarly, the CPOC is mainly situated within the range of 0.13 to 0.19 yuan/kWh, averaging at 0.16 yuan/kWh. Both average values are below the empirical benchmarks. If there is an opportunity to further enhance the system’s EER, the CPOC for cooling capacity would consequently decrease.

Performing PCA and dimensionality reduction on the input parameters of the cooling load forecasting model, we derived principal components from the correlation matrix R (with variables standardized). Out of the 44 input parameters (principal components), 6 principal components were ultimately selected as the input parameters for the cooling load forecasting model. The eigenvalues of these six principal components are all greater than 1, and the cumulative variance contribution rate reaches 87.26%, essentially encapsulating the information contained in all parameters. During the dimensionality reduction process, we obtained a confidence ellipse with a 95% confidence interval. With a probability as high as 95%, this represents a smaller three-dimensional ellipse spanned by PC1 (50.6%), PC2 (18.1%), and PC3 (8.8%). This indicates that the parameter values we selected closely approximate the true values, suggesting a high precision in estimating population parameters based on sample parameters.

In this research, we initially delved into the correlation between the data and the system’s performance to discern their interrelationships. Upon acquiring a comprehensive understanding of the system, we proceeded to investigate hourly cooling load forecasting using a data-driven approach. Among the commonly employed algorithms for cooling load forecasting, the lightGBM algorithm demonstrated superior performance. Specifically, the RAE value for the lightGBM algorithm was 3.62%, the RRMSE value stood at 42.75%, and the R² value reached 92.96%. By analyzing, organizing, and categorizing the data, as well as refining and optimizing the input parameters of the load forecasting model, we lay the groundwork for further research and insights. This paves the way for enhancing the operational efficiency of regional cooling systems, achieving cost savings, and elevating the comfort levels for terminal-users through various optimization measures.

Acknowledgements

The authors would like to thank the Chongqing Construction Science and Technology Plan Project of China (2019 No. 2-2-3, No. 2-2-4) for their financial support.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	IEA. World Energy Outlook 2022. http://www.iea.org
[2]	B-Open/Aureliana Barghini. Copernicus Climate Change Service Application Documentation Global Temperature Trend Monitor 2021). http://www.copernicus.eu
[3]	Angelidis, O., Ioannou, A., Friedrich D., Thomson, A. and Falcone, G. (2023) District Heating and Cooling Networks with Decentralised Energy Substations: Opportunities and Barriers for Holistic Energy System Decarbonisation. Energy, 269. https://doi.org/10.1016/j.energy.2023.126740
[4]	bp. bp Energy Outlook 2023. http://www.bp.com
[5]	Zhao, J. and Liu, X. (2018) A Hybrid Method of Dynamic Cooling and Heating Load Forecasting for Office Buildings Based on Artificial Intelligence and Regression Analysis. Energ Buildings, 174, 293-308. https://doi.org/10.1016/j.enbuild.2018.06.050
[6]	Gelazanskas, L. and Gamage, K.A.A. (2014) Demand Side Management in Smart Grid: A Review and Proposals for Future Direction. Sustain Cities Soc, 11, 22-30. https://doi.org/10.1016/j.scs.2013.11.001
[7]	Lei, Y., Wang, D., Jia, H., Chen, J., Li, J., Song, Y., et al. (2020) Multi-Objective Stochastic Expansion Planning Based on Multi-Dimensional Correlation Scenario Generation Method for Regional Integrated Energy system Integrated Renewable Energy. Appl Energ, 276. https://doi.org/10.1016/j.apenergy.2020.115395
[8]	Luo, N., Hong, T.Z., Li, H., Jia, R.X. and Weng, W.G. (2017) Data Analytics and Optimization of an Ice-Based Energy Storage System for Commercial Buildings. Appl Energ, 204, 459-475. https://doi.org/10.1016/j.apenergy.2017.07.048
[9]	Inayat, A. and Raza, M. (2019) District Cooling System via Renewable Energy Sources: A Review. Renew Sust Energ Rev, 107, 360-373. https://doi.org/10.1016/j.rser.2019.03.023
[10]	Vandermeulen, A., van der Heijde, B. and Helsen, L. (2018) Controlling District Heating and Cooling Networks to Unlock Flexibility: A Review. Energy, 151, 103- 115. https://doi.org/10.1016/j.energy.2018.03.034
[11]	Sanaye, S. and Hekmatian, M. (2016) Ice Thermal Energy Storage (ITES) for Air- Conditioning Application in Full and Partial Load Operating Modes. International Journal of Refrigeration-Revue Internationale Du Froid, 66, 181-197. https://doi.org/10.1016/j.ijrefrig.2015.10.014
[12]	Sait, H.H. and Selim, A.M. (2014) Charging and Discharging Characteristics of Cool Thermal Energy Storage System with Horizontal Pipes Using Water as Phase Change Material. Energy Conversion and Management, 77, 755-762. https://doi.org/10.1016/j.enconman.2013.10.034
[13]	Ruan, Y.J., Liu, Q.R., Li, Z.W. and Wu, J.Z. (2016) Optimization and Analysis of Building Combined Cooling, Heating and Power (BCHP) Plants with Chilled Ice Thermal Storage System. Appl Energ, 179, 738-754. https://doi.org/10.1016/j.apenergy.2016.07.009
[14]	Kang, X., Wang, X., An, J. and Yan, D. (2022) A Novel Approach of Day-Ahead Cooling Load Prediction and Optimal Control for Ice-Based Thermal Energy Storage (TES) System in Commercial Buildings. Energ Buildings, 275. https://doi.org/10.1016/j.enbuild.2022.112478
[15]	Hu, J., Zheng, W., Zhang, S., Li, H., Liu, Z., Zhang, G., et al. (2021) Thermal Load Prediction and Operation Optimization of Office Building with a Zone-Level Artificial Neural Network and Rule-Based Control. Appl Energ, 300. https://doi.org/10.1016/j.apenergy.2021.117429
[16]	Wang, Z., Hong, T. and Piette, M.A. (2020) Building Thermal Load Prediction through Shallow Machine Learning and Deep Learning. Appl Energ, 263. https://doi.org/10.1016/j.apenergy.2020.114683
[17]	Xu, Y., Li, F. and Asgari, A. (2022) Prediction and Optimization of Heating and Cooling Loads in a Residential Building Based on Multi-Layer Perceptron Neural Network and Different Optimization Algorithms. Energy, 240. https://doi.org/10.1016/j.energy.2021.122692
[18]	Ben-Nakhi, A.E. and Mahmoud, M.A. (2004) Cooling Load Prediction for Buildings Using General Regression Neural Networks. Energ Convers Manage, 45, 2127-2141. https://doi.org/10.1016/j.enconman.2003.10.009
[19]	Cao, L., Li, Y., Zhang, J., Jiang, Y., Han, Y. and Wei, J. (2020) Electrical Load Prediction of Healthcare Buildings through Single and Ensemble Learning. Energy Rep, 6, 2751-2767. https://doi.org/10.1016/j.egyr.2020.10.005
[20]	Li, Q., Meng, Q., Cai, J., Yoshino, H. and Mochida, A. (2009) Predicting Hourly Cooling Load in the Building: A Comparison of Support Vector Machine and Different Artificial Neural Networks. Energ Convers Manage, 50, 90-96. https://doi.org/10.1016/j.enconman.2008.08.033
[21]	Cui, B., Gao, D.-C., Xiao, F. and Wang, S. (2017) Model-Based Optimal Design of Active Cool Thermal Energy Storage for Maximal Life-Cycle Cost Saving from Demand Management in Commercial Buildings. Appl Energ, 201, 382-396. https://doi.org/10.1016/j.apenergy.2016.12.035
[22]	Shan, K., Fan, C. and Wang, J.Y. (2019) Model Predictive Control for Thermal Energy Storage Assisted Large Central Cooling Systems. Energy, 179, 916-927. https://doi.org/10.1016/j.energy.2019.04.178
[23]	Tang, H., Yu, J., Geng Y., Liu, X. and Lin, B. (2023) Optimization of Operational Strategy for Ice Thermal Energy Storage in a District Cooling System Based on Model Predictive Control. J Energy Storage, 62. https://doi.org/10.1016/j.est.2023.106872
[24]	Cao, H., Lin, J. and Li, N. (2023) Optimal Control and Energy Efficiency Evaluation of District Ice Storage System. Energy, 276. https://doi.org/10.1016/j.energy.2023.127598
[25]	Ma, X., Sha, J., Wang, D., Yu, Y., Yang, Q. and Niu, X. (2018) Study on a Prediction of P2P Network Loan Default Based on the Machine Learning LightGBM and XGboost Algorithms According to Different High Dimensional Data Cleaning. Electronic Commerce Research and Applications, 31, 24-39. https://doi.org/10.1016/j.elerap.2018.08.002
[26]	Heine, K., Tabares-Velasco, P.C. and Deru, M. (2021) Design and Dispatch Optimization of Packaged Ice Storage Systems within a Connected Community. Appl Energ, 298. https://doi.org/10.1016/j.apenergy.2021.117147