S. Huang, X. Lu, W. Zuo, X. Zhang, C. Liang 2019. "Model-based optimal operationof heating tower heat pump systems." Building and Environment, 160, 106199.DOI: el-based optimal operation of heating tower heat pumpsystemsShifang Huanga, b, Xing Lub, Wangda Zuob, Xiaosong Zhanga, *, Caihua LiangaabSchool of Energy and Environment, Southeast University, Nanjing, 210096, PR ChinaDepartment of Civil, Environmental and Architectural Engineering, University of Colorado Boulder, Boulder, CO 80309, U.S.A.* Corresponding author. E-mail address: [email protected] (Xiaosong Zhang)Abstract: In current applications of heating tower heat pumps (HTHPs), the systems tend to run with constant speed or fixed setpoints, which can be inefficient under varying weather data and building loads. To address this issue, this study proposes a model-basedoptimal operation of the HTHPs to achieve energy savings in both cooling and heating modes. Firstly, a physics-based model for anexisting HTHP system was developed. Then, artificial neural network (ANN) models were developed and trained with vast amount ofoperational data generated by the physics-based model. The ANN models were found to be highly accurate (average relative error lessthan 1%) and computationally efficient (about 300 times faster than the physics-based model). After that, three optimal approacheswere proposed to minimize the total energy consumption of the HTHP system. Approach 1 optimizes the load distribution betweendifferent heat pump units. Approach 2 optimizes the speed of fans and pumps by fixed approach and range of the condenser water (orevaporator solution). Approach 3 optimizes both the load distribution and the speed of fans and pumps. The optimization isimplemented by using the ANN models, proposed approaches, and a genetic algorithm via a case study. The results show that theenergy savings in the cooling season are 2.7%, 11.4%, and 14.8% by the three approaches, respectively. In the heating season, theenergy savings of the three approaches are 1.6%, -1.4%, and 4.7%, respectively. Moreover, the thermodynamic performance in typicaldays was analyzed to investigate how energy savings could be achieved.Key words: heating tower heat pump; model-based optimization; ANN model; genetic algorithm; energy saving1. IntroductionChillers with cooling towers are the most widely used systems for cooling supply in large commercial buildingsdue to their high energy efficiency [1]. Inspired by the heat rejection process in a standard cooling tower, previousresearch proposed a reversibly-used cooling tower, named as heating tower [2,3]. This tower can be used to transferheat from ambient air to water or solution with low temperature. Due to the high demand of building heating supplyin cool or cold regions, researchers proposed to replace the chillers with heat pumps and coupled them with theheating towers to fulfill the heating demand. The mean coefficient of performance (COP) of the heating tower heatpumps (HTHPs) in winter varied from 2.24 to 4.12 for different climate zones [4,5] And in summer, the HTHPcould function as a water-cooled chiller with high efficiency [6]. The HTHP has shown significant advantages overthe conventional heat pumps. For example, compared with the air-source heat pump (ASHP), the HTHP shows1

higher efficiency and completely addresses the frosting problem. Comparing to the ground-source heat pump(GSHP), the HTHP has lower geographic limitation and lower initial cost. Therefore, the HTHP can be served asan alternative to the conventional heat pump for building cooling and heating in practical applications [7].However, in previous research and applications of the HTHP [4,7], all components of the HTHP were runningwith constant speeds or fixed set points, which can be inefficient under varying weather conditions and buildingloads. Currently, few research has addressed the optimal operation of the HTHP system. Researches in optimaloperation of water-cooled chiller systems can be used as a reference because these two systems have similarities insummer. For the optimization of the chiller plants, the near-optimal control and model-based optimal control weremainly investigated. The near-optimal control developed by Yu et al. [8], adopted simple polynomials to realizedirect or indirect control of fan speed and pump speed. The wet-bulb temperature and part load ratio (PLR) ofchillers were taken as the input variables in the polynomials. Those methods are simple and fast, which can be easilyapplied in the chiller plants with limited computational resources. However, achieving the energy savings cannot beguaranteed due to the simplification in system modeling. To improve the accuracy, Huang et al. [9] employed aBayesian network model trained by optimal control setpoints computed by optimization algorithm with a highfidelity physic model. The results showed that the energy saving ratio was improved by about 5.1% compared withthe polynomial prediction, but still 0.2% lower than the model-based optimization method. Another problem in thenear optimization is that the PLR of the chiller is calculated by a fixed maximum cooling capacity, instead of theactual capacity which is changing with the operating conditions. As a result, the optimal number of operating chillersdetermined by this method may be inappropriate.To achieve the highest energy saving ratio, many researchers performed model-based optimizations for chillerplants. Three key steps of the model-based optimization are model development, control strategy selection, andoptimization algorithm implementation. The following part will discuss each step in detail.For the model development, regression models, physics-based models, and data-driven models are mostly used.The regression models are the fastest in computing speed, which then require much fewer computing efforts inmodel-based optimization [10,11]. However, they are reliable only for operating conditions within the range of theregression data, and extrapolation outside this range may lead to significant errors. In addition, all the existingregression models for chiller plants cannot be used to predict the performance of the HTHP system because theyhave different physical processes in winter. The physics-based models have no such limitations in analyzing theperformance of the HTHP system, since the components models, system configuration, and operation mode can bechanged by employing physical laws [12,13]. Although the physics-based models are effective, they always requirehigh computational costs. When combined with global optimization techniques to seek the optimal solutions, theslow computing speed can prevent their online applications. To address this problem, data-driven models areproposed and adopted, which can achieve both fast computing speed and high accuracy [14,15]. However, it requiresa vast amount of operational data to train the data-driven models. This is even more difficult for the newly sized or2

installed the HTHP systems. In this paper, to address the above mentioned issues, we developed a physics-basedmodel and validated it with a small amount of experimental data. Then, this new physics-based model was run underdifferent combinations of weather and operation conditions to generate a vast amount of operational data for thedata-driven models (artificial neural network models).For the control strategy selection, the optimal load distribution and optimal fan and pump speeds were mainlydiscussed in the previous studies. Those two approaches were combined in some studies to achieve the highestpossible energy saving ratio. The optimal load distribution method works because the chiller’s COPs varies with itsPLRs [16]. On one side, chillers, whose PLRs are adjusted by inlet guide vanes or modulating sliding valves, havethe highest COP at full load due to the lower mechanical loss [17]. On the other side, chillers, whose PLR areadjusted by variable frequency drive (VFD), perform better at low PLRs than high PLRs due to the higher thermalheat exchange efficiency at low PLR [17]. The optimal load distribution method can still work when the chillers arereplaced by heat pumps in the HTHP systems since the heat pump’s COP is also a function of its PLR. In the optimalfan and pump speeds method, various indirect control variables (e.g. condenser water range, cooling tower approach,and condenser water supply temperature) are adopted to adjust the fan and pump speeds. Reducing the fan and pumpspeed at low PLRs can significantly reduce their energy consumption with little negative effects on the performanceof chillers since the heat rejection capacity of the cooling tower is redundant at low PLRs. In the HTHP system, theheat rejection capacity of the tower in summer was about 3-4 times of its heat absorption capacity in winter [18].Thus, the heating tower are sized to satisfy the heat absorption in winter. As a result, it has more redundant heatrejection capacity than the conventional cooling towers, which indicates that the heating tower has higher energysaving potential than conventional cooling towers. In winter, the energy use of the solution regeneration is anadditional factor, which can significantly impact the optimization results.For the optimization algorithm, the gradient-based algorithm was often adopted, including reduced gradientsearch, Hooke–Jeeves search, and Lagrangian method [10,11]. These methods are fast in convergence, but they canbe trapped into local optima and require the gradient of objective function. Due to these limitations, they are oftenused in simple regression and need to couple with starting point selecting methods. The heuristic algorithms, suchas genetic algorithm (GA), simulated annealing, and particle swarm optimization, were more popular in theoptimization of chiller plants because they do not require the gradient of objective function [19,20]. Since theoptimization of HTHP is a multi-variable, non-differentiable optimization problem, this study selects the GA whichis a popular heuristic algorithm successfully used for the optimization of chiller plants.In this paper, we described an existing HTHP and developed the physics-based model for it. Then, the artificialneural network (ANN) models were developed and trained with vast amount of operational data generated by thephysics-based model. After that, the optimization problem, optimal approaches, and optimization algorithm wereintroduced. Using the GA and ANN models, the model-based optimization was conducted in a case study. Finally,3

the annual energy savings of the proposed approaches were presented. Moreover, the thermodynamic performanceof typical days was analyzed to investigate how energy savings could be achieved.2. Modeling of the HTHP system2.1. System descriptionFig. 1 shows the schematics of a typical HTHP system in summer and winter conditions, respectively. ThisHTHP system consists of three identical HTHPs, which are practical HTHPs. Each HTHP has one heat pump, oneheating tower, one tower side pump, and one user side pump. All equipment is dedicated to that HTHP. To investigatethe energy-saving potential of different optimal approaches, this study adopted variable frequency drives for thecompressors, tower fans, and tower side pumps.In the summer condition, the HTHP runs as a water chiller with cooling tower. The condenser of the heat pumpis connected with the heating tower which is used as cooling tower in summer with water as the working fluid. Heatrejection of the heat pump is done by the water evaporation in the tower, and the mass balance in this process isachieved by adding make-up water to the loop. The evaporator is connected with the building system to providechilled water for the buildings.In the winter condition, the connection between components is switched by external valves. Opposite to thesummer condition, the evaporator is connected with the heating tower in winter, and water is replaced by solutionwith low freezing point (e.g. glycol aqueous) to avoid system freeze. The solution absorbs both sensible and latentheat from the atmosphere, and provides heat to the heat pump. Since the solution can be diluted in this process, asolution regeneration system based on vacuum boiling and condensation is equipped to achieve mass balance in thisstudy. The details of the regeneration system can be found in our previous studies. This regeneration device isindependent from the HTHP system, and the power used by it is to satisfy the gasification latent heat of water.Therefore, the regeneration device power has no influence on the heating capacity of the HTHP system[21,22].(a) Summer condition(b) Winter conditionFig. 1. The schematic of the HTHP system (HT, heating tower; HP, heat pump; TSP, tower side pump; USP, userside pump; CWS/CWR, supply/return condenser water; CHWS/CHWR, supply/return chilled water; SS/SR,supply/return solution; HWS/HWR, supply/return hot water; VFD, variable frequency drive.)4

2.2. System performance indexesFor the heating tower, the heat and mass transfer capacities are adopted:𝑄𝑠 (πΆπ‘π‘Ž πœ”π‘Ž 𝐢𝑝𝑣 ) π‘€π‘Ž (π‘‡π‘Ž,π‘œ π‘‡π‘Ž,𝑖 ) ,(1)𝑄𝑙 π‘Ÿπ‘£ π‘€π‘Ž (πœ”π‘Ž,π‘œ πœ”π‘Ž,𝑖 ) ,(2)where 𝑄𝑠 is the sensible heat transfer capacity, and 𝑄𝑙 is the latent heat transfer capacity. Here, the positive valuesof 𝑄𝑠 and 𝑄𝑙 mean that the heat and mass transfer directions are from air to condenser water or solution. When thevalues are negative, the directions are the opposite. The 𝐢𝑝, 𝑀, 𝑇, πœ”, and π‘Ÿ represent the specific heat, mass flowrate, temperature, humidity ratio, and vaporization latent heat, respectively. The subscripts π‘Ž, 𝑣, 𝑖, and π‘œ representthe air, water vapor, tower inlet, and tower outlet, respectively.Besides, the approach, 𝜏, and range, π›₯𝑇, which are commonly used in the conventional cooling towers, areemployed. For the winter condition, the wet-bulb temperature, 𝑇𝑀𝑏 , is replaced by the equivalent wet-bulbtemperature, 𝑇𝑀𝑏,π‘’π‘ž , at which the equivalent enthalpy of solution is equal to the enthalpy of the ambient air [4].𝑇𝑐𝑀𝑠 𝑇𝑀𝑏 , summer condition𝜏 {,𝑇𝑀𝑏,π‘’π‘ž 𝑇𝑠𝑠 , winter condition(3)𝑇 𝑇𝑐𝑀𝑠 , summer conditionπ›₯𝑇 { π‘π‘€π‘Ÿ.𝑇𝑠𝑠 π‘‡π‘ π‘Ÿ ,winter condition(4)The energy performance of the heat pump is measured by the coefficient of performance, 𝐢𝑂𝑃:𝐢𝑝𝑀 𝑀𝑀 (π‘‡π‘β„Žπ‘€π‘Ÿ π‘‡π‘β„Žπ‘€π‘  )𝐢𝑂𝑃 {π‘Šπ»π‘ƒπΆπ‘π‘€ 𝑀𝑀 (π‘‡β„Žπ‘€π‘  π‘‡β„Žπ‘€π‘Ÿ )π‘Šπ»π‘ƒ, summer condition,,winter condition(5)where subscript 𝑀 represents the water. The π‘Šπ»π‘ƒ is the power consumption of the heat pump.The energy performance of the HTHP system is indicated by energy efficiency ratio, 𝐸𝐸𝑅:𝐢𝑝𝑀 𝑀𝑀 (π‘‡π‘β„Žπ‘€π‘Ÿ π‘‡π‘β„Žπ‘€π‘  )𝐸𝐸𝑅 {π‘Šπ»π‘ƒ π‘Šπ‘‡π‘†π‘ƒ π‘Šπ»π‘‡πΆπ‘π‘€ 𝑀𝑀 (π‘‡β„Žπ‘€π‘  π‘‡β„Žπ‘€π‘Ÿ ),π‘Šπ»π‘ƒ π‘Šπ‘‡π‘†π‘ƒ π‘Šπ»π‘‡ π‘Šπ‘…π·summer condition,, winter condition(6)where π‘Šπ‘‡π‘†π‘ƒ , π‘Šπ»π‘‡ , and π‘Šπ‘…π· are the power consumption of the tower side pump, fan, and regeneration device,respectively.2.3. Physics-based modelThe physics-based models of all equipment of the studied HTHP system were developed separately, and thencoupled together using energy and mass balances between different equipment. The following sections describe themodels for heat pump, heating tower, fan and pumps.2.3.1. Heat pump modelThe heat pump consists of four main components, including the scroll compressors, shell-tube evaporator,shell-tube condenser, and thermostatic expansion valve. Models of the components are developed as follows.5

CompressorThe refrigerant mass flow rate, 𝑀𝑅 , and power consumption, π‘Šπ‘π‘œπ‘šπ‘ , for fixed speed scroll compressors canbe expressed by a function of the evaporating temperature, 𝑇𝑒 , and condensing temperature, 𝑇𝑐 [23].𝑀𝑅,π‘Ÿπ‘Žπ‘‘π‘’π‘‘ πœ‹1 πœ‹2 𝑇𝑒 πœ‹3 𝑇𝑐 πœ‹4 𝑇𝑒2 πœ‹5 𝑇𝑒 𝑇𝑐 πœ‹6 𝑇𝑐2 πœ‹7 𝑇𝑒3 πœ‹8 𝑇𝑒2 𝑇𝑐 πœ‹9 𝑇𝑒 𝑇𝑐2 πœ‹10 𝑇𝑐3 ,(7)π‘Šπ‘π‘œπ‘šπ‘,π‘Ÿπ‘Žπ‘‘π‘’π‘‘ πœ‹1 πœ‹2 𝑇𝑒 πœ‹3 𝑇𝑐 πœ‹4 𝑇𝑒2 πœ‹5 𝑇𝑒 𝑇𝑐 πœ‹6 𝑇𝑐2 πœ‹7 𝑇𝑒3 πœ‹8 𝑇𝑒2 𝑇𝑐 πœ‹9 𝑇𝑒 𝑇𝑐2 πœ‹10 𝑇𝑐3 ,(8)where the subscript π‘Ÿπ‘Žπ‘‘π‘’π‘‘ represents the performance under rated speed. The regression coefficients are listed inTable 1.Table 1. Regression coefficients of Eqs. (1) and (2)Regression coefficients of Eq. (1)Regression coefficients of Eq. (2)Ο€12.696 102Ο€14.505Ο€28.461Ο€23.547 10-2Ο€3-1.881Ο€31.107 10-1Ο€41.261 10-1Ο€45.832 10-6Ο€5-3.693 10-2Ο€58.214 10-5Ο€65.741 10-2Ο€6-2.706 10-4Ο€71.154 10-3Ο€7-8.487 10-7Ο€8-1.006 10-3Ο€81.705 10-6Ο€99.011 10-4Ο€9-9.203 10-6Ο€10-7.020 10-4Ο€102.340 10-5For variable speed scroll compressors, 𝑀𝑅 and π‘Šπ‘π‘œπ‘šπ‘ were found proportional to the rotation speed, π‘π‘π‘œπ‘šπ‘[24]. The above models with speed ratio correction can be used to express the performance of the variable speedscroll compressors [25]:𝑀𝑅 𝑀𝑅,π‘Ÿπ‘Žπ‘‘π‘’π‘‘ οΏ½οΏ½οΏ½π‘’π‘‘π‘Šπ‘π‘œπ‘šπ‘ π‘Šπ‘π‘œπ‘šπ‘,π‘Ÿπ‘Žπ‘‘π‘’π‘‘ ���𝑑𝑒𝑑(10).Evaporator and condenserThe classical logarithmic mean temperature difference method is adopted in both evaporator and condensermodels. The cooling capacity of the evaporator, 𝑄𝑒 , and the heating capacity of the condenser, 𝑄𝑐 , can be expressedas follows:𝑄𝑒 𝐾𝑒 𝐴𝑒 𝐿𝑀𝑇𝐷𝑒 ,(11)𝑄𝑐 𝐾𝑐 𝐴𝑐 𝐿𝑀𝑇𝐷𝑐 ,(12)where 𝐾 , 𝐴, and 𝐿𝑀𝑇𝐷 are the heat transfer coefficient, heat transfer area, and logarithmic mean temperaturedifference between refrigerant and water/ solution. The subscript 𝑒 represents the evaporator, and 𝑐 represents thecondenser.6

For the shell-tube evaporator, the refrigerant evaporates inside the tube, and the water/solution flows acrossthe tube. The heat transfer coefficient for the evaporation of R22 inside the tube, 𝐾𝑅 , can be described as follows:πœ†π‘…,𝑙0.80.4𝐾𝑅 51 π‘₯ 0.8 πœŒπ‘…,𝑔𝑐0 (π‘₯()πœŒπ‘…,𝑙)25𝐺 2𝑐5π‘žπ‘4) [𝑐1 (𝑐0 )𝑐2 (9.8𝜌2 𝑅𝑑 ) 2.2𝑐3 (𝐺 𝑅,𝑖) ] ,π‘Ÿπ‘…,𝑙 𝑖(13)𝑅 𝑅(14),where 𝑅𝑒 is the Reynolds number, and π‘ƒπ‘Ÿ is the Prandtl number. The πœ† , 𝐺 , 𝜌 , π‘ž , and π‘₯ are the thermalconductivity coefficient, mass flow flux, density, inner heat flux, and dryness, respectively. The subscript 𝑅 is therefrigerant. The subscript 𝑙 and 𝑔 represent the liquid phase and gas phase, respectively. The 𝑑𝑖 is the innerdiameter of the tube. The 𝑐0 is the characteristic number of convection heat transfer, and 𝑐1 -𝑐5 are constantsdepending on 𝑐0 . The heat transfer coefficient of the water/solution across the tube, 𝐾𝑀 (or 𝐾𝑠 ), can be expressedas follow:1/3 πœ†π‘€π‘‘0𝐾𝑀 0.22𝑅𝑒𝑀0.6 π‘ƒπ‘Ÿπ‘€,1/3 πœ†π‘ π‘‘0𝐾𝑠 0.22𝑅𝑒𝑠0.6 π‘ƒπ‘Ÿπ‘ (15),(16)where the subscript 𝑠 represents the solution. The π‘‘π‘œ is the external diameter of the tube.For the shell-tube condenser, the refrigerant condenses outside the tube, and the water flows inside the tube.The heat transfer coefficient for the condensation of R22 outside the tube, 𝐾𝑅 , and the heat transfer coefficient ofthe water flowing inside the tube, 𝐾𝑀 , are shown as follows:29.8πœ†3𝑅,𝑙 πœŒπ‘…,π‘™π‘Ÿπ‘…πΎπ‘… 0.725 (𝑒𝑅,𝑙𝐾𝑀 0.023𝑅𝑒𝑀0.8 π‘ƒπ‘Ÿπ‘€0.40.25)πœ†π‘€π‘‘π‘–π‘‘π‘œ 0.25 (𝑇𝑐 π‘‡π‘€π‘Žπ‘™π‘™ ) 0.25 𝛹1 πœ€1 ,,(17)(18)where 𝑒𝑅,𝑙 is liquid phase dynamic viscosity of the refrigerant, and π‘‡π‘€π‘Žπ‘™π‘™ is the temperature of the tube wall. The𝛹1 and πœ€1 are correction factors depending on the size and arrangement of tubes.Based on the coefficients of refrigerant, water and solution for different processes, the overall heat transfercoefficient for the evaporator or condenser, 𝐾𝑒 (𝐾𝑐 ), can be expressed as a function of the heat transfer coefficientinside the tube, 𝐾𝑖 , and outside the tube, πΎπ‘œ :𝐾𝑒 (𝐾𝑐 ) 11πΎπ‘–π΄π‘œπ›Ώπ΄π‘œ1 π‘…π‘œ 𝐴𝑖 πœ†π‘€π‘Žπ‘™π‘™ π΄π‘šπΎπ‘œ( 𝑅𝑖 ),(19)where 𝑅 is the heat transfer resistance. The 𝛿 is the thickness of the wall. The subscripts 𝑖 and π‘œ represent theproperty inside and outside the tube, respectively. And the subscripts π‘š is the mean value of the inside and outsideproperty.The heat transfer capacities of the evaporator and condenser can be expressed by the energy variations ofrefrigerant and water/solution.𝑄𝑒 𝑀𝑅 (β„Ž1 β„Ž4 ) ,(20)7

𝑄𝑒 𝑀𝑀 𝐢𝑝𝑀 (π‘‡π‘β„Žπ‘€π‘Ÿ π‘‡π‘β„Žπ‘€π‘  ) or 𝑄𝑒 𝑀𝑠 𝐢𝑝𝑠 (𝑇𝑠𝑠 π‘‡π‘ π‘Ÿ ) ,(21)𝑄𝑐 𝑀𝑅 (β„Ž2 β„Ž3 ) ,(22)𝑄𝑐 𝑀𝑀 𝐢𝑝𝑀 (π‘‡π‘π‘€π‘Ÿ 𝑇𝑐𝑀𝑠 ) π‘œπ‘Ÿ 𝑄𝑐 𝑀𝑀 𝐢𝑝𝑀 (π‘‡β„Žπ‘€π‘  π‘‡β„Žπ‘€π‘Ÿ ) ,(23)where β„Ž1 and β„Ž4 are the enthalpy of the refrigerant in the inlet and outlet of the evaporator, β„Ž2 and β„Ž3 are theenthalpy of the refrigerant in the inlet and outlet of the condenser.Expansion valveThe expansion process in the expansion valve is taken as an isenthalpic process as shown in Eq.(24). The massflow rate of the refrigerant can be calculated by Eq. (25) [26].β„Ž 3 β„Ž4 ,(24)𝑀𝑅 𝐢𝐷 π΄π‘‘β„Ž πœŒπ‘…,𝑙 (𝑃𝑐 𝑃𝑒 ),(25)where, 𝐢𝐷 is the constant mass flow coefficient. The π΄π‘‘β„Ž represents the geometric throat area of the thermostaticexpansion, which is adjustable and controlled by the superheat. The 𝑃𝑐 and 𝑃𝑒 are the pressure of condenser andevaporator, respectively.2.3.2. Heating tower modelThe model of the heating tower in winter is developed using a finite difference method [3]. Eqs.(26) and (27)express the energy and mass balances between air and solution. Eq.(28) describes the solute balance of the solution.π‘šπ‘Ž π‘‘β„Žπ‘Ž 𝐢𝑝𝑠 π‘šπ‘  𝑑𝑇𝑠 𝐢𝑝𝑠 𝑇𝑠 π‘‘π‘šπ‘  ,(26)π‘‘π‘šπ‘  π‘šπ‘Ž π‘‘πœ”π‘Ž ,(27)𝑋𝑠 π‘šπ‘  (𝑋𝑠 𝑑𝑋𝑠 )(π‘šπ‘  π‘‘π‘šπ‘  ) ,(28)where π‘‘β„Žπ‘Ž is the enthalpy variation of the air through an element. The 𝑑𝑇𝑠 , π‘‘π‘šπ‘  and 𝑑𝑋𝑠 represent the variationsof the solution in temperature, mass flow rate and concentration through an element, respectively. The 𝑋𝑠 is themass concentration of the solution.The convective heat and mass transfer are also applied as shown in Eqs.β„Žπ‘ 𝐿 𝑑π‘₯ 𝑑𝑦 𝛼𝑀 (𝑇𝑠 π‘‡π‘Ž ) π‘šπ‘Ž (πΆπ‘π‘Ž πœ”π‘Ž 𝐢𝑝𝑣 )π‘‘π‘‡π‘Ž ,(29)β„Žπ‘‘ 𝐿 𝑑π‘₯ 𝑑𝑦 𝛼𝑀 (πœ”π‘  πœ”π‘Ž ) π‘šπ‘Ž π‘‘πœ”π‘Ž ,(30)where πœ”π‘  is the equivalent humidity ratio of the solution, π‘‘π‘‡π‘Ž and π‘‘πœ”π‘Ž are the temperature and humidity ratiovariation of air through an element. The 𝑑π‘₯ 𝑑𝑦 represents the size of each element, 𝐿 is the length of the packing,and 𝛼𝑀 is the specific area of the packing. The β„Žπ‘ is the heat transfer coefficient, β„Žπ‘‘ is the mass transfercoefficient. These two coefficients are expressed as functions of the solution mass flow flux, 𝐺𝑠 , and air mass flowflux, πΊπ‘Ž , using iteration method and data from a practical heating tower [3].πœ‹πœ‹(31)πœ‹πœ‹(32)β„Žπ‘ πœ‹1 𝐺𝑠 2 πΊπ‘Ž 3 ,β„Žπ‘‘ πœ‹1 𝐺𝑠 2 πΊπ‘Ž 3 .8

By replacing the subscript 𝑠 with 𝑀 and setting 𝑋𝑠 to zero, the model listed above can also be used tosimulate the performance in summer.2.3.3. Pump and fan modelThe characteristic curves of pump and fan are often provided by manufacturers, while the characteristic curvesof pipeline are difficult to know in the practice. Therefore, it is hard to determine the power consumption of pumpand fan, which is calculated by the equipment characteristic curve and the pipeline characteristic curve. In this study,regression models were adopted to calculate the performance of fan and pump.π‘Šπ‘‡π‘†π‘ƒ πœ‹1 πœ‹2 (π‘Šπ»π‘‡ πœ‹1 πœ‹2 𝑑𝑁𝐻𝑇𝑁𝐻𝑇,π‘Ÿπ‘Žπ‘‘π‘’π‘‘) πœ‹3 �𝑑) πœ‹3 () πœ‹4 οΏ½) πœ‹4 οΏ½))(33),(34).And the data in our previous study was adopted to fit the coefficients of the above equations [4]. The results arelisted in Table 2.Table 2. Coefficients of Eqs. (33) and (34)Regression coefficientπœ‹1πœ‹2πœ‹3πœ‹4Eq. (33)1.761 10-48.264 10-1-7.857 10-13.707Eq. (34)3.088 10-4-1.558 10-13.8522.483 10-12.4 Data-driven modelThe physics-based model adopting physical laws has high accuracy. However, it is time consuming and needs alot of computing resources, since it contains several iterations and a great number of grids for the finite differencescheme applied on the heating tower model. This study used a Lenovo T450 computer with a four-core processor(I5-5200U, 2.2 GHz). The physics-based model developed in this study takes about 41.4 s to simulate one workingcondition in summer (cooling mode), and 101.7 s for one working condition in winter (heating mode). This isunacceptable in the model-based optimization since we need to run hundreds of simulations to identify the controlsettings for one working condition, and it takes several hours. Therefore, it is helpful to replace the complex physicalmodels for the heat pump and heating tower by fast computing models. We still use the models for pump and fan asthey are simple and easy to compute. The ANN models were selected due to its high accuracy and fast computing.2.4.1. ANN model developmentTo develop the ANN models, we first define the dependent and independent variables for each equipment. Forthe heat pump model in summer, as shown in Fig. 2 (a), five independent variables are identified, including thespeed of the compressor (π‘π‘π‘œπ‘šπ‘ ), speed of tower side pump (π‘π‘π‘’π‘šπ‘ π‘œπ‘Ÿ 𝑀𝑐𝑀 ), speed of the user side pump(π‘π‘π‘’π‘šπ‘ π‘œπ‘Ÿ π‘€π‘β„Žπ‘€ ), condenser water supply temperature (𝑇𝑐𝑀𝑠 ), and chilled water supply temperature (π‘‡π‘β„Žπ‘€π‘  ). Thereare two dependent variables, such as the coefficient of performance (𝐢𝑂𝑃) and cooling capacity of the heat pump9

(𝑄𝑒 ). The other parameters (e.g. the condenser water return temperature, chilled water return temperature, andcapacity of the condenser) can be calculated by the independent and dependent variables. Therefore, the input matrixcan be expressed as [π‘π‘π‘œπ‘šπ‘ , 𝑀𝑐𝑀 , π‘€π‘β„Žπ‘€ , 𝑇𝑐𝑀𝑠 , π‘‡π‘β„Žπ‘€π‘  ], and the output matrix is [𝐢𝑂𝑃, 𝑄𝑒 ]. In winter, besides the fiveindependent variables mentioned above, the concentration of the solution ( 𝑋𝑠 ) should also be added as anindependent variable. Thus, the inputs and outputs can be expressed as [π‘π‘π‘œπ‘šπ‘ , π‘€β„Žπ‘€ , 𝑀𝑠 , π‘‡β„Žπ‘€π‘  , 𝑇𝑠𝑠 , 𝑋𝑠 ] and[𝐢𝑂𝑃, 𝑄𝑐 ].For the heating tower model in summer, independent variables are identified as the outdoor air dry-bulb (𝑇𝑑𝑏 ),wet-bulb temperature (𝑇𝑀𝑏 ), speed of the tower fan (π‘π‘“π‘Žπ‘› π‘œπ‘Ÿ π‘€π‘Ž ), mass flow rate of the condenser water (𝑀𝑐𝑀 ), andcondenser water return temperature (π‘‡π‘π‘€π‘Ÿ ). The sensible heat transfer capacity (𝑄𝑠 ) and latent heat transfer capacity(𝑄𝑙 ) of the tower are the dependent variables. Therefore, the inputs and outputs of the tower model can be expressedas [𝑇𝑑𝑏 , 𝑇𝑀𝑏 , π‘€π‘Ž , 𝑀𝑐𝑀 , π‘‡π‘π‘€π‘Ÿ ] and [𝑄𝑠 , 𝑄𝑙 ], respectively. Similarly, the inputs and outputs of the tower model inwinter are [𝑇𝑑𝑏 , 𝑇𝑀𝑏 , π‘€π‘Ž , 𝑀𝑠 , π‘‡π‘ π‘Ÿ , 𝑋𝑠 ] and [𝑄𝑠 , 𝑄𝑙 ]. Different from the conventional cooling tower model which onlytakes the 𝑇𝑀𝑏 as input weather data, the 𝑇𝑑𝑏 is also considered here, since both heat and mass transfer cannot beignored in the heating tower.According to recommended numbers of the nodes and hidden layers [27], a back-propagation (BP) networkwith two hidden layers and ten nodes in each hidden layer was employed. The BP network was implemented underthe MATLAB environment. We tried several commonly used activation functions (e.g. Sigmold, tanh, ReLu), andselected the functions by assessing the calculation accuracy and time of the ANNs. The Tansig function was adoptedas the transfer function between the first three layers, and Purelin function served the last two layers, respectively.The structures of the BP networks adopted for the four models are presented in Fig. 2.(a) Heat pump(b) Heating towerFig. 2. The structure of the BP network for modeling2.4.2. Model training and validationFor each ANN model, 18,000 data sets were used for training, the rest 2,000 data sets were adopted forvalidation. In order to improve the prediction agreement, all the data were normalized into the range of -1 to 1. The10

Levenberg-Marquardt acted as the training function, and the mean squared error was used as the performance index,during the training process. This section introduces the generation of training and validation data, and the resultsare shown in Section 4.2.20,000 input data sets for heat pump model in summer were generated randomly in the following ranges:π‘π‘π‘œπ‘šπ‘ [580, 2900] Rev. min-1, 𝑀𝑐𝑀 [18.5, 37] m3 h-1, π‘€π‘β„Žπ‘€ [11, 22] m3 h-1, 𝑇𝑐𝑀𝑠 [15, 35] C, π‘‡π‘β„Žπ‘€π‘  [7,9] C. Similarly, the ranges of the 20,000 input data sets for heat pump model in winter are π‘π‘π‘œπ‘šπ‘ [580,2900]Rev. min-1, π‘€β„Žπ‘€ [11,22] m3 h-1, 𝑀𝑠 [18.5,37] m3 h-1, π‘‡β„Žπ‘€π‘  [44,46] C, 𝑇𝑠𝑠 [-12,15] C, and 𝑋𝑠 [0,0.4].After the generation of each 𝑋𝑠 , the freezing point of the solution will be verified to make sure it is lower than thesolution temperature. For the heating tower model in summer, the ranges of the 20,000 data sets are 𝑇𝑑𝑏 [15,38] C, π‘€π‘Ž [4300,43000] m3 h-1, 𝑀𝑐𝑀 [18.5,37] m3 h-1, π‘‡π‘π‘€π‘Ÿ [20,40] C. For the heating tower model inwinter, the ranges of the 20,000 data sets are 𝑇𝑑𝑏 [-6,20] C, π‘€π‘Ž [4300,43000] m3 h-1, 𝑀𝑠 [18.5,37] m3 h-1,π‘‡π‘ π‘Ÿ [-13,15] C, 𝑋𝑠 [0,0.4]. Particularly, the 𝑇𝑀𝑏 is calculated by 𝑇𝑑𝑏 and the relative humidity (πœ‘π‘Ž ) from 15%to 100%, to avoid physical violations. With the input data sets mentioned above, the corresponding outputs werecomputed by the physics-based model in both summer and winter conditions.3. Model-based optimization of the HTHP system3.1. Optimization formulationThe aim of the optimization work is to minimize the power consumption of the HTHP system with timevarying building load and weather data. For a typical HTHP system described in Section 2.1, the major energyconsumers in summer are heat pumps, tower side pumps, and tower fans. In winter, besides the above threeequipment, the energy consumption of the regeneration drive should be also counted. Thus, the objective functionscan be expressed as:π½π‘ π‘’π‘šπ‘šπ‘’π‘Ÿ min( 3𝑖 1 π‘Šπ»π‘ƒ,𝑖 3𝑖 1 π‘Šπ‘‡π‘†π‘ƒ,𝑖 3𝑖 1 π‘Šπ»π‘‡,𝑖 ) ,(35)π½π‘€π‘–π‘›π‘‘π‘’π‘Ÿ min( 3𝑖 1 π‘Šπ»π‘‡,𝑖 3𝑖 1 π‘Šπ‘‡π‘†π‘ƒ,𝑖 3𝑖 1 π‘Šπ»π‘‡,𝑖 π‘Šπ‘…π· ) .(36)When the heat pump model, heating tower model and building system are coupled together, some interactionconstraints need to be satisfied. The actual cooling/heating load provided by the system should be the same as thebuilding load. In addition, the properties of the condenser water (or solution) entering and exiting the heat pumpshould be the same as those of the heating tower. Besides, the mechanical constraints provided by the manufacturersshould be 𝑑 π‘π‘π‘œπ‘šπ‘ 100%π‘π‘π‘œπ‘šπ‘,π‘Ÿπ‘Žπ‘‘π‘’π‘‘ ,(37)50%𝑀𝑐𝑀,π‘Ÿπ‘Žπ‘‘π‘’π‘‘ 𝑀𝑐𝑀 100%𝑀𝑐𝑀,π‘Ÿπ‘Žπ‘‘π‘’π‘‘ ,(38)50%𝑀𝑠,π‘Ÿπ‘Žπ‘‘π‘’π‘‘ 𝑀𝑠 100%𝑀𝑠,π‘Ÿπ‘Žπ‘‘π‘’π‘‘ ,(39)10%π‘€π‘Ž,π‘Ÿπ‘Žπ‘‘π‘’π‘‘ π‘€π‘Ž 100%π‘€π‘Ž,π‘Ÿπ‘Žπ‘‘π‘’π‘‘ .(40)11

In addition, some constraints of the operation parameters should also be satisfied. In the summer condition, the𝑇𝑐𝑀𝑠 should be higher than 15 C as required by the manufacturers to avoid the small pressure difference of thecompressor. In the winter condition, the freezing point of the solution should be lower than the π‘‡π‘ π‘Ÿ , which is thelowest temperature in the system, to prevent system from freeze. Here, the penalty function is adopted to addressthe constrained optimization problems. By introducing the pe

compressors, tower fans, and tower side pumps. In the summer condition, the HTHP runs as a water chiller with cooling tower. The condenser of the heat pump is connected with the heating tower which is used as cooling tower