Energy 55 (2013) 392e402
Contents lists available at SciVerse ScienceDirect
Energy journal homepage: www.elsevier.com/locate/energy
Model-based optimization for vapor compression refrigeration cycle Lei Zhao, Wenjian Cai*, Xudong Ding, Weichung Chang School of Electrical and Electronic Engineering, Nanyang Technological University, Singapore 639798, Singapore
a r t i c l e i n f o
a b s t r a c t
Article history: Received 26 September 2012 Received in revised form 25 February 2013 Accepted 26 February 2013 Available online 26 April 2013
This paper presents a model-based optimization strategy for vapor compression refrigeration cycle. Through analyzing each component characteristics and interactions within the cycle, the optimization problem is formulated as minimizing the total operating cost of the energy consuming devices subject to the constraints of mechanical limitations, component interactions, environment conditions and cooling load demands. A MGA (modified genetic algorithm) together with a solution strategy for a group of nonlinear equations is proposed to obtain optimal set point under different operating conditions. Simulation studies are conducted to compare the proposed method with traditional oneoff control strategy to evaluate its performance. Experiment results of a real practical system are also presented to demonstrate its feasibility. Ó 2013 Elsevier Ltd. All rights reserved.
Keywords: Vapor compression refrigeration cycle Hybrid components models Global optimization Modified genetic algorithm System simulation and testing
1. Introduction Air conditioning and refrigeration systems, as common residential and industry devices, are widely used to transfer heat between different locations [1]. They control space temperatures for comfort of human life, food storage and transportation, etc., and consume a lot of energy for both developing and developed countries. For instance, statistical data shows air conditioner and refrigerator account for 28% of home energy consumption in US [2]. For hot and humid tropical country such as Singapore, this ratio can even rise to over 50% [3]. Among all types of cooling systems, electricity-based vapor compression cooling systems are still dominant in the current market. Since the vapor compression refrigeration cycle is the heart of and consumes most of the energy in any cooling system, the effort to reduce the energy consumption through system control and optimization in vapor compression cooling system is of practical significance due to both energy shortage and global warming concerns. All optimization techniques of vapor compression refrigeration cycle involve components and system modeling. Even though there are some discussions on the modeling of compressors, expansion valves and overall systems, the main works in this area is still focused on the development of heat exchanger models. So far, several models have developed including: distributed model, NTU
* Corresponding author. Tel.: þ65 6790 6862; fax: þ65 6793 3318. E-mail address:
[email protected] (W. Cai). 0360-5442/$ e see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.energy.2013.02.071
(Number of Transfer Units) model, black box model and hybrid model. In the distributed model, the heat exchanger is divided into consecutive segments, the inlet of each segment equals to the outlet of its former adjacent segment. The final static state is determined by recursive calculation. Jia et al. built a distributed model for dry expansion evaporators to analyze their steady and dynamic characteristics [4]. To accurately design plate-fin heat exchanger, Zhang et al. developed a 3D distributed parameter model [5]. In 3-NTU model, the static states of heat exchanger are determined by a dimensional parameter called the heat transfer effectiveness 3 . Ren presented an accurate model of counterflow cooling tower with modified 3-NTU method [6]. Arsenyeva et al. used 3-NTU model of plate and frame heat exchanger in design procedure to improve heat recovery [7]. The black box model uses a predetermined function format as its heat transfer function, the coefficients are determined by curve fitting or identification techniques according to experiment data. Ertunc and Hosoz proposed an adaptive neuro-fuzzy inference system to predict the heat rejection rate and the temperatures of refrigerant and leaving air of evaporative condensers. Simulation results showed that the prediction errors are within 5% [8]. Kim et al. derived a polynomial equation with seven variables to predict the effectiveness of a pilot dry coil indirect evaporative cooler based on experiment data [9]. The hybrid model lumps the physical information as constants, calculate the model coefficients through least-square methods based on experiment data, therefore it takes advantages of both physical and empirical models, which makes it suitable for online optimization with high accuracy. Wang et al. derived the hybrid
L. Zhao et al. / Energy 55 (2013) 392e402
Nomenclature Av c f F K H Hc,fg Hc,r,i Hc,r,o He,fg He,r,i He,r,o _ c;air m _ c;air;nom m _ e;air m _ e;air;nom m _r m _ r;max m _ r;min m n Pc Pc,max Pc,min Pe Pe,max Pe,min Q_ c Q_ e Q_ com Q_ req
Tc,air,i Tc,max
opening percentage of electronic expansion valve constants of hybrid model function frequency of fans (Hz) penalty enthalpy (kJ/kg) enthalpy difference of gas and liquid saturated refrigerant in condenser (kJ/kg) condenser inlet refrigerant enthalpy (kJ/kg) condenser outlet refrigerant enthalpy (kJ/kg) enthalpy difference of gas and liquid saturated refrigerant in evaporator (kJ/kg) evaporator inlet refrigerant enthalpy (kJ/kg) evaporator outlet refrigerant enthalpy (kJ/kg) air flow rate of condenser (kg/s) nominal air flow rate of condenser (kg/s) air flow rate of evaporator (kg/s) nominal air flow rate of evaporator (kg/s) refrigerant mass flow rate (kg/s) maximal refrigerant mass flow rate (kg/s) minimal refrigerant mass flow rate (kg/s) polytropic exponent condenser refrigerant saturated pressure (bar) maximal condenser saturated pressure allowed (bar) minimal condenser saturated pressure allowed (bar) evaporator refrigerant saturated pressure (bar) maximal evaporator saturated pressure allowed (bar) minimal evaporator saturated pressure allowed (bar) condenser energy exchange rate (kJ/kg) evaporator energy exchanger rate (kJ/kg) mechanical work of compressor (kJ/kg) required cooling load (kJ/kg) condenser inlet air temperature ( C) maximal refrigerant saturation temperature in condenser ( C)
model of cooling coil based on the flow characteristics and nonlinear least-square method [10]. Considering the refrigerant phase change in heat exchanger, Ding et al. developed the hybrid model for evaporator with three parameters [11], they further applied the hybrid method to analyze condenser that contains gas, liquid and two phase refrigerant, derived the condenser hybrid model including four parameters [12]. Aiming at ES (energy saving) vapor compression refrigeration cycles, many research works have been done on improving the system COP (coefficient of performance). Stoecker claimed to achieve higher COP, both superheat and subcool should be decreased to near zero, this setting will increase the efficiency of heat exchangers [13]. Jensen and Skogestad proposed to add the refrigerant volume (active charge) in the system as an additional control variable, they also discussed the effect of subcooling on system efficiency [14]. To realize self-optimizing control, they discussed the selection of the controlled variable to improve the system efficiency and performance of control method [15]. In Ref. [16], they focused on the available degrees of freedom since it is the basis of optimal operation and plant wide control. Larsen et al. proposed a gradient method to find the suboptimal solution for condensing pressure, while keeping the superheat and evaporating pressure constant [17]. Selbas et al. optimized the components of subcooled
393
Tc,min
minimum refrigerant saturation temperature in condenser ( C) Tc,r,i condenser inlet refrigerant temperature ( C) Tc,sc condenser subcool temperature ( C) Tc,r,sat condenser refrigerant saturated temperature ( C) Tc,r,i condenser refrigerant inlet temperature ( C) Tc,r,o condenser refrigerant outlet temperature ( C) Te,r,sat refrigerant saturated temperature in evaporator ( C) Te,max maximal refrigerant saturation temperature in evaporator ( C) Te,min minimal refrigerant saturation temperature in evaporator ( C) Te,sh evaporator superheat temperature ( C) Va volume at bottom dead center (m3) Vd volume when suction valve open (m3) _ W c;fan condenser fan power (kW) _ W c;fan;nom condenser fan power when air flow rate is mc,air,nom (kW) _ com electricity power consumption of compressor (kW) W _ W evaporator fan power (kW) c;fan _ W e;fan;nom evaporator fan power when air flow rate is me,air,nom (kW) _ total power (kW) W total h enthalpy delivery efficiency u compressor rotation speed (r/s) r inlet refrigerant density (kg/m3) Subscripts air feature of air c condenser com compressor e evaporator ev expansion valve i inlet o outlet m mass flow rate q heat exchange h enthalpy delivery efficiency
and superheated vapor compression refrigeration cycle separately based on exergy analysis [18]. The entransy dissipation-based optimization technique was introduced to find the optimal design and to minimize the operation cost of evaporative cooling system by Yuan and Chen [19]. Since system states of vapor compression refrigeration cycle are severely interacted, it is highly possible that solutions of classical methods are trapped in local minimum instead of global minimum. On the other hand, the newly developed artificial intelligence computation techniques; have been successfully used in optimization of air conditioning and HVAC (heating, ventilation, and air conditioning) systems. Fong et al. employed differential evolution to optimize solar thermal refrigeration system [20]. Kusiak et al. proposed to optimize the commercial HVAC system with evolutionary computation algorithm [21]. They also used multi-objective particle swarm algorithm to find the optimal solution that can balance energy consumption and thermal comfort [22]. Kusiak and Xu then used dynamic neural network to construct black box models for system performance prediction with controllable and uncontrollable inputs and outputs [23]. Simulation studies based on experiment conducted at the ERS (Energy Resource Station) showed that system level optimization that they proposed can improve overall system operating performance significantly.
394
L. Zhao et al. / Energy 55 (2013) 392e402
Furthermore, Hovgaard et al. developed a novel economic model predictive control scheme to reduce the operating cost of supermarket refrigeration system [24]. In this paper, we propose a modified genetic optimization technique to tackle the optimization problem of vapor compression refrigeration cycle. The objective function is formulated based on mathematical models of the major components. The thermodynamic characteristics of the heat exchangers are predicted with the hybrid models, which are simple yet accurate with some physical significance. The power consumptions of the power consuming equipments are predicted also using the hybrid models. A MGA (modified genetic algorithm) is developed to search the optimal settings. Simulation and experimental results on a lab scale pilot plant demonstrate that a significant operating cost can be saved by the proposed method. 2. Working principle and component models The vapor compression refrigeration cycle is used to remove heat from a space of lower temperature (cold reservoir) to an environment of high temperature (hot reservoir). As shown schematically in Fig. 1, it basically consists of four components e Evaporator, Compressor, Condenser, and Expansion Valve which are connected in a closed loop so that the refrigerant is continuously circulated. A peh chart reflecting the variation of the refrigerant states is shown as in Fig. 2 which is briefly explained as below. 1. Refrigerant starts from state 1, inlet of compressor, in gas phase at low pressure and low temperature. As the refrigerant exits
the compressor at state 2, it is still in gas phase but with high pressure (saturation pressure) and high temperature. The superheated vapor refrigerant enters the condenser at state 2, where the heat from the working fluid is transferred to secondary fluid as the secondary fluid has relatively lower temperature. Inside the condenser, the refrigerant is in two phases (liquid mixed with gas). When the refrigerant exits the condenser after transferring the heat to the secondary fluid, the temperature is below the condensation temperature to ensure the fully liquidated refrigerant. 2. The refrigerant in liquid phase enters the expansion valve at the state 3 where the pressure of the fluid is decreased through the throttling effect of the expansion valve. The sudden decrease in the pressure of the refrigerant causes it to partially evaporate which in turn reduces the saturation temperature. Thus, when the refrigerant exits the valve it is in low pressure, low temperature and two phases (liquid mixed with gas). 3. At state 4, the refrigerant enters the evaporator at low temperature and low pressure condition. Since the refrigerant temperature is lower than the temperature of cold reservoir, the refrigerant completely evaporates at a constant temperature with the help of evaporator fan to increase the heat transfer rate. At the exit of the evaporator, the temperature is slightly above the evaporating temperature for the safe operation of compressor. The mathematical models of the four components which realize the transition of the refrigerant between the states and are instrumental in the model-based optimization are given below.
Fig. 1. Vapor compression refrigeration cycle.
L. Zhao et al. / Energy 55 (2013) 392e402
395
2.3. Compressor In a compressor, the mass flow rate and the refrigerant energy change during compression stroke can be expressed [25,26]:
ccom;m;3 Pc u Pe
(5)
n1 Pc n _r 1 um Pe
(6)
_r ¼ m
ccom;m;1 ccom;m;2
and
Q_ com ¼
Fig. 2. peh chart of vapor compression cycle.
2.1. Evaporator Based on mass and energy balance, a simple hybrid model to describe the heat transfer properties in evaporator is given as [11]
_ r þ ce;1 m _ cr e;3 Te;air;i Te;r;sat He;g He;r;i m Q_ e ¼ !ce;3 _r m 1 þ ce;2 _ e;air m
(1)
where ccom,m,1, ccom,m,2 and ccom,m,3 are constants determined by curve fitting, Pc, Pe, u, n, Va, Vd, and Q_ com are condensing temperature, evaporating pressure, compressor rotation speed, polytropic exponent, volume at bottom dead center, volume when suction valve open and mechanical work input to refrigerant by compressor respectively. Since the parameters n, Va and Vd are constants for a given compressor and at a specified working environment, Eq. (6) can be further simplified to a hybrid model form:
_ r Pe Q_ com ¼ ccom;q;1 um
ccom;q;2 Pc 1 Pe
(7)
where ccom;q;1 ¼ n=n 1ðVa Vd Þ and ccom;q;2 ¼ n 1=n.
where ce,1, ce,2 and ce,3 are constants obtained by fitting experiment _ r, m _ e;air , Te,air,i, Te,r,sat and Q_ e are the enthalpy of data, He,g, He,r,i, m saturated gas phase refrigerant in evaporator, refrigerant enthalpy of evaporator inlet, mass flow rates of refrigerant, air outside evaporator, temperature of inlet air and saturated refrigerant of evaporator, heat exchanging rate of evaporator, respectively (see Appendix B for the calculations of He,g, He,r,i and Te,r,sat). In addition, an energy balance equation, i.e., energy absorbed by the refrigerant is equal to energy reduced in the cold reservoir, i.e.
_ r He;r;o He;r;i Q_ e ¼ m
n ðVa Vd ÞPe n1
(2)
will also be used in later development, where He,r,o is refrigerant enthalpy at evaporator outlet (see Appendix B for the calculation of He,r,o).
2.4. Expansion valve The mass flow rate of expansion valve is determined by valve opening percentage, pressure difference and inlet refrigerant density. Its mass flow rate is given by [25]
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi _ r ¼ cev;1 þ cev;2 Av rðPc Pe Þ m
(8)
where cev,1 and cev,2 are constants, Av and r are opening percentage of electronic expansion valve and density of inlet refrigerant respectively (see Appendix B for the calculation of r). Since the expansion process is isenthalpic, it is assumed that refrigerant enthalpy is constant which implies Qev ¼ 0. 2.5. Energy balance of overall cycle
2.2. Condenser Similar to the evaporator, condenser can also be represented by a hybrid model [12]:
_ cr c;4 Tc;r;sat Tc;air;i þ cc;2 m _ r Tc;r;i Tc;r;sat þ Hc;fg m _r cc;1 m _ Qc ¼ !cc;4 _r m 1 þ cc;3 _ c;air m (3) where cc,1, cc,2, cc,3 and cc,4 are constants calculated by fitting experiment data, Hc,fg is enthalpy difference between saturated _ c;air is the mass flow liquid and gas phase refrigerant in condenser, m rate of air outside condenser. Tc,r,sat, Tc,r,i, Tc,air,i are temperature of saturated refrigerant, inlet refrigerant and inlet air of condenser, Q_ c is the heat transferring rate of condenser (see Appendix B for the calculations of Hc,fg and Tc,r,sat). The corresponding energy balance equation is
_ r Hc;r;i Hc;r;o Q_ c ¼ m
(4)
where Hc,r,i and Hc,r,o are enthalpy of inlet and outlet refrigerant (see Appendix B for the calculations of Hc,r,i, and Hc,r,o).
Finally, the corresponding energy balance equation for the overall system is given as
Q_ c ¼ Q_ com þ Q_ e
(9)
3. Optimization problem formulation The objective of global optimization for vapor compression refrigeration cycle is to satisfy cooling load requirement of cold reservoir with minimal energy consumption. Mathematically, it can be formulated as
_ _ _ _ Min W total ¼ W com þ W c;fan þ W e;fan _ _ Subject to Q e ¼ Q req
(10)
_ _ _ _ _ where W total ; W com ; W c;fan ; W e;fan and Q req are total power consumption, power consumption of compressor, condenser fan power consumption, evaporator fan power consumption and required cooling load respectively. The power consumption models of the compressor, condenser fan and evaporator fan are formulated according to their working principles.
396
L. Zhao et al. / Energy 55 (2013) 392e402
3. Condensing and evaporating pressure
3.1. Power consumption of compressor Using a hybrid model to describe the delivery coefficient of compressor, hcom [26]
Pc;min Pc Pc;max
(16a)
hcom ¼ ccom;h;1 þ ccom;h;2 ðPc =Pe Þccom;h;3
Pe;min Pe Pe;max
(16b)
(11)
the power consumption of compressor is given as
_ com ¼ W
where Pc,min, Pc,max, Pe,min, and Pe,max are minimal and maximal condensing pressure, minimal and maximal evaporating pressure.
Q_ com
(12)
hcom
4. Superheat
where ccom;h;1 ; ccom;h;2 and ccom;h;3 are constants calculated by catalog or experiment data. 3.2. Power consumptions of condenser fan and evaporator fan The power consumptions of fans are influenced by two parameters: mass flow rates of fluids and the pressure difference between the inlets and outlets and can be described [27]:
_ _ W c;fan ¼ W c;fan;nom cc;fan;0 þ cc;fan;1 þ cc;fan;2
_ c;air m _ c;air;nom m
_ c;air m _ c;air;nom m
!2 þ cc;fan;3
7 + C Te;sh 25 + C
(17)
Superheat temperature too low will cause oscillation in the cycle [28], while the COP of system will decrease steeply if the superheat temperature is too high [13]. The upper and lower limits depend on the system configurations. 5. Subcool
! Tc;sc 0
_ c;air m _ c;air;nom m
!3 ! (13a)
(18)
Subcool only has lower bound to ensure the refrigerant at expansion valve inlet is at liquid phase for flashing. 6. Compressor frequency (for the testing system)
_ _ W e;fan ¼ W e;fan;nom ce;fan;0 þ ce;fan;1 þ ce;fan;2
_ e;air m _ me;air;nom
_ e;air m _ e;air;nom m
!2 þ ce;fan;3
!
_ e;air m _ me;air;nom
30 Hz F 50 Hz !3 ! (13b)
(19)
This physical constraint is directly determined by compressor working range. 7. Mass flow rate of refrigerant
where cc,fan,0 to cc,fan,3 and ce,fan,0 to ce,fan,3 are coefficients deter_ _ _ mined by catalog or experiment data, W c;fan ; W e;fan ; W c;fan;nom ; _ _ _ W e;fan;nom ; mc;air;nom and me;air;nom are measured power consumptions of condenser fan, measured power consumption of evaporator fan, nominal power consumption of condenser fan, nominal power consumption of evaporator fan, nominal mass flow rate of condenser fan and mass flow rate of evaporator fan respectively. In the optimization problem, there are two types of constraints, i.e., interaction between the components represented by Eqs. (1)e (5) and (7)e(9), and physical limitations of the components and fluids, including: 1. Air mass flow rates of condenser and evaporator fans
_r m _ r;max _ r;min m m
(20)
_ r;min and m _ r;max are minimal and maximal refrigerant where m mass flow rate. The upper bound of mass flow rate of refrigerant is determined by expansion valve and compressor capacities, the lower bound is determined by flow meter working range. 8. EEV opening percentage
0 < Av 1
(21)
To reduce the variables and simplify the solution procedures, we notice that
_ c;air m _ c;air;max _ c;air;min m m
(14a)
_ e;air m _ e;air;max _ e;air;min m m
(14b)
_ c;air;min ; m _ c;air;max ; m _ e;air;min and m _ e;air;max are minimal and where m maximal condenser air mass flow rates, minimal and maximal evaporator air mass flow rates respectively. 2. Condensing and evaporating temperatures
Tc;air;i Tc;r;o Tc;min Tc;r;sat Tc;max Tc;r;i
(15a)
Te;r;i Te;min Te;r;sat Te;max Te;air;i
(15b)
where Tc,min, Tc,max, Te,min and Te,max are minimal and maximal condensing temperatures, minimal and maximal evaporating temperatures.
Pc, Hc,fg and Tc,r,sat Pe, He,fg and Te,r,sat Hc,r,o and He,r,i have same physical properties, therefore, we can use Pc; Pe; and Hc,r,o to replace Tc,r,sat and Hc,fg; Te,r,sat and He,fg; He,r,i, in the system objective, constraints and power consumption equations. All variables used in the optimization problem, according to their rules, can be further classified into three groups: Uncontrollable variables (Q_ e , Tc,air,i, Te,air,i): They are determined by user demand or environmental conditions, these variables are fixed during each calculation step. Among these three variables, the outdoor air temperature, Tc,air,i, is the most critical data, the other two uncontrollable variables, Te,air,i and Q_ e , are also affected by Tc,air,i.
L. Zhao et al. / Energy 55 (2013) 392e402
Independent variables (Pc, Pe, u, Tc,r,o): Independent variables can be selected randomly within the constraint limits to meet the optimization objective. The obvious choices for indepen_ c;air ; m _ e;air as they are directly dent variables are u; Av ; m controlled by compressor, EEV (electronic expansion valve) and fan frequencies. However, Pe, Pc, u and Tc,r,o are assigned to be independent variables in this work in order to reduce the complexity of solving process as well as to increase solution accuracy. _ r; T ; m _ _ Dependent variables (m ; Te;r;o ; m ; Q_ ; Q_ c;r;i
c;air
e;air
c
com
and Av): These values can be resolved through the interaction and physical constraint equations for the given independent variables. These variables can then be used to calculate the power consumption equations.
397
a constant which can be specified by user, the rule of thumb of the selection is 3 times of the rated compressor power. By calculating Eq. (21), the fitness value of each individual chromosome is converted into predetermined interval to control the evolution speed in an acceptable range. 4.3. Evolution In evolution, four major evolutionary operations: selection, crossover, mutation and reinsertion are performed. The elitist roulette wheel selection [29], single point crossover and single bit mutation operators are employed. The crossover and mutation points are selected randomly in each generation. 4.4. Termination
4. Solution procedures To solve the constrained nonlinear optimization problem, a modified genetic algorithm (MGA) combined with a solution strategy for a group of nonlinear equations is devised for higher efficiency and less computation complexity. This algorithm is divided into four parts: encoding, construction of fitness function, evolution and termination [29], which are briefly discussed below.
The genetic algorithm terminates if one of the following criterion is satisfied: 1. Predetermined maximal generation number is reached.
4.1. Encoding
start
Each independent variable in constraint Eqs. (15a), (16a), (16b) and (19) is converted into a binary string. All of these strings are combined together to form a chromosome. Gray coding is preferred for overcoming the hidden representative bias in binary representation. The lengths of the binary strings are determined by the control precision of the corresponding variables: the more precise set point control requires the longer binary string. For instance, if the compressor controller’s accuracy is 102, then 10 bits is enough for u to describe the compressor setting. After encoding, a random of population of chromosomes, is generated such that, when decoded to real values, the variable values lie within their limitations. The whole generation should possibly try to include the entire range of values.
determine component model parameters
measure air temperature and calculate cooling load
reduce variables by practical consideration initialize parameters of modified genetic algorithm
gen>max gen? Y
4.2. Constructing fitness function
N
After setting the random variables, eight dependent variables are solved according to interaction constraint equations based on the coded independent variables (see Appendix B). Since the _ objective of this optimization problem is to minimize W total , the fitness function should be constructed to give the maximum value _ to the smallest W total . In addition, the fitness function should reflect whether the chromosome violates the constraints. Based on these considerations, the fitness function is constructed as follows:
fitness ¼
1 _ _ _ com þ W W e;fan þ W c;fan þ Ktotal
gen=gen+1
calculate dependent variables
calcuate fitness
evolution and store the best one
(21) termination criterion satisfied?
where the definition of total system penalty function Ktotal is as following:
Ktotal ¼
8 > > > <
0
8 > P > > Ki : cp þ i¼1
if if
8 P i¼1 8 P
N
Y find indiviual with largest fitness among all generations
Ki ¼ 0
i¼1
(22)
calculate power consumption of optimal set point
Ki s0 end
where Ki, i ¼ 1,.,8 are the individual state penalty functions (see Appendices A and B for the detailed definition and calculation), cp is
Fig. 3. Procedure flowchart.
398
L. Zhao et al. / Energy 55 (2013) 392e402
Table 1 Physical limits of variables.
1.5
Variables
Lower bound
Upper bound
Fc,fan Fe,fan Tc,r,sat Te,r,sat Pc Pe Te,sh Tc,sc Fcom _r m
15 Hz 15 Hz Tc,air,i 12 C 8 bar 2 bar 5 C 0 C 30 Hz 0.008 kg/s
35 Hz 30 Hz Tc,r,i Te,air,i 15 bar 5 bar 25 C NA 50 Hz 0.03 kg/s
MGA on-off control
1.45
Energy consumption(kW)
1.4 1.35 1.3 1.25 1.2 1.15
33
1.1 1.05
32
8
Temperature(°C)
12
14 Time(h)
16
18
20
Fig. 6. Energy consumptions for working hours.
31
30
4.5. Remark
29
28
27
10
0
4
8
12
16
20
24 28 Time(h)
32
36
40
44
48
Fig. 4. Measured temperatures for two consecutive days.
2. The optimal fitness of group does not change much over several steps, the bound for determining whether its value changes too much is predetermined according to experience. 3. All chromosomes converge to nearly the same fitness groups. Compared with the traditional genetic algorithm, there is a difference of the proposed algorithm, particularly modified for such kind of problems, which is discussed in the following remark.
A main drawback of traditional genetic algorithm when applied to vapor compression refrigeration cycle is that when Q_ e is approaching to system cooling capacity, it is highly possible that all individuals in last generation do not satisfy the physical constraints, even though several individuals in previous generation do, the good chromosome might lost after evolution. The proposed algorithm is to store the individuals that satisfy the physical constraints with maximum fitness value in each generation step when Q_ e is large. As a result, only chromosome with maximal fitness is stored for each generation and to compare with the final generation, this modification increases the possibility of finding a feasible solution. The flowchart of MGA for optimization of vapor compressor cycle is illustrated in Fig. 3 which can be summarized as follows: Step 1: Determine the parameters of the component models by catalog or experimental data. Step 2: Measure Te,air,i and Tc,air,i, and calculate Q_ e .
25 4.6 4.4
20 Percentage of time interval(%)
4.2
Cooling load(kW)
4 3.8 3.6 3.4 3.2
15
10
5
3 2.8 2.6
0
4
8
12
16
20
24 28 Time(h)
32
36
40
Fig. 5. Assumed cooling load profiles in simulation.
44
48
0 29.5
30
30.5
31 31.5 Temperature(°C)
32
Fig. 7. Bar chart of one day temperature profile.
32.5
33
L. Zhao et al. / Energy 55 (2013) 392e402
399
3.5 Energy consumptions:MGA Energy savings:MGA
Energy consumption(kWh)
3
2.5
2
1.5
1
0.5
0 29.5
30
30.5
31 31.5 Temperature(°C)
32
32.5
33
Fig. 8. System energy consumptions and savings.
Step 3: Initialize the parameters of MGA: population size, maximum number of generations, precisions of variables, generation gap and probabilities of crossover and mutation. Step 4: Perform the evolutionary operations. Step 5: Repeat Step 5 and store chromosome according to Remark until termination criterion is satisfied. Step 6: Find the chromosome with largest fitness value among all generations and calculate corresponding power consumption. After calculating the solution, control system will set the frequencies of compressor, evaporator fan, condenser fan and expansion valve opening according to the solution, resulting in all operating states equal or close to their corresponding optimal value.
Fig. 10. Picture of the lab scale vapor compression refrigeration system.
Considering the large influences of evaporating and condensing temperatures on operating states and system energy consumption, the fans will be used to slightly adjust them to decrease the error caused by model inaccuracy. In practice, the precisions of the state variables in optimization algorithm are usually higher than that of fans’ and compressor’s controllers, the closest values that controllers can achieve will be chosen as the set point. This approach can minimize the difference between desired and practical output.
Fig. 9. Schematic of the vapor compression refrigeration system.
400
L. Zhao et al. / Energy 55 (2013) 392e402
schemes have the same trends, but MGA is more sensitive to the change of cooling load. The largest and smallest energy saving of MGA compared with oneoff control scheme is 1.44% and 14.78% respectively. To analyze the energy consumptions at different temperature conditions, the measured temperatures are regrouped according to the temperature range. The percentage of temperature in the measured day from 8am to 19pm is shown as in Fig. 7. To estimate the energy saving effect of MGA, the total energy consumptions for both MGA and oneoff control schemes of the temperature profile in Fig. 7 are calculated and presented in Fig. 8, which shows that the average energy saving of 8.52% is achieved.
1.5 MGA test on-off test MGA sim on-off sim
1.45
Energy consumption(kW)
1.4 1.35 1.3 1.25 1.2 1.15
6. Experiment
1.1 1.05
8
10
12
14 Time(h)
16
18
20
Fig. 11. Energy consumptions.
5. Simulation The simulation studies are based on the models of a lab scale pilot plant. The coefficients of models are determined through catalog fitting 200 groups of experiment data which are obtained by different compressor, condenser fan and evaporator fan frequencies. All constraints also follow the system physical limits listed as in Table 1. To simulate the energy saving effect in practical situations, environmental temperatures of two consecutive days are measured as shown in Fig. 4. Assuming the set point of indoor temperature is 23 C and the highest cooling load is equal to the maximum cooling capacity of the pilot plant, the measured temperature data can then be proportionally converted to the cooling load shown as in Fig. 5. The simulation results for energy consumption of the proposed MGA and traditional oneoff control scheme from 8am to 19pm are illustrated in Fig. 6. It can be seen that energy consumptions of MGA are lower than that of oneoff control scheme for all cooling loads, but the difference is smaller when the cooling load is near its maximum cooling capacity. The energy consumptions of the two
To verify the effect of the proposed MGA, experimental tests are conducted on a lab scale multi-evaporator vapor compression cycle. The schematic and photo of the experimental system are illustrated in Figs. 9 and 10, respectively. The test bench consists of a single-stage vapor compression formed by basic components: variable speed compressor, fin-tube condenser, thermostatic expansion valve and isolated fin-tube evaporator. The heat exchangers work with secondary fluids loops, which represents and simulates environmental and thermal load. The experiment system includes a semi-hermetic reciprocating compressor, an air-cooled finned-tube condenser, three electronic expansion valves and three evaporators (one air-cooled finnedtube evaporator and two electronic evaporators). One air duct heater controls the inlet air temperature of condenser for simulating outdoor condition, and the inlet air temperature of evaporator is constantly kept as 25 C by HVAC system. The working fluid used for the system is R134a. The compressor, the condenser fan and the evaporator fan are equipped with inverters to adjust their corresponding frequencies. An air duct heater is installed in front of the condenser to control the temperature of condenser inlet air. In addition, several temperature and pressure sensors are installed for detecting state variables. To measure the mass flow rate, a flow meter with maximum 4% full scale error is used. The measurement range of the pressure transducers and the temperature transmitters are 0e16 bar and 40 Ce200 C, with their maximum full scale
Fig. 12. System energy consumptions and savings.
L. Zhao et al. / Energy 55 (2013) 392e402
error are within 0.3% and 0.3 C respectively. The positions of these sensors are as following (refer to Fig. 9): refrigerant temperature and pressure at compressor outlet (T1 and P1); temperature and pressure at condenser inlet (T2 and P2); temperature and pressure at condenser outlet (T3 and P3); temperature and pressure at EEV inlet (T5 and P5); temperature and pressure at evaporator outlet (T6 and P6); _ r ). mass flow rate after receiver (m humidities at evaporator inlet and outlet (R.H.1, R.H.2). To simplify the test and calculation, the average temperatures and cooling load demands of each hour from 8am to 19pm are used during testing. The energy consumptions of both MGA and oneoff control schemes for results of both simulation and testing are compared in Fig. 11. It shows that EC (energy consumptions) are higher in the experiment results than that of simulation for nearly all testing groups, the average and largest difference between experiment and simulation results are 3.43% and 10.47% respectively, which might be caused by the modeling errors. However, the results are accurate enough for the MGA to be implemented in practical applications. Moreover, the difference between experimental results of MGA and oneoff control scheme verifies the energy saving effect of proposed method. The smallest and largest differences between these two groups of data are 0.03% and 13.14% respectively with the average overall energy saving (ES) of 8.45% as shown in Fig. 12. 7. Conclusion In this paper, a model-based optimization strategy for vapor compression refrigeration cycle was developed, which was formulated as minimizing the total operating cost of the energy consuming devices subject to the constraints of mechanical limitations, component interactions, environment conditions and cooling load demands. A modified genetic algorithm together with a solution strategy for a group of nonlinear equations was proposed to effectively obtain the optimal set point under different operating conditions. The experimental results showed that the set points calculated by MGA can reduce energy consumption compared with traditional oneoff control, the overall energy saving for a typical day during 8ame19pm is 8.45%. Furthermore, the average difference between simulated and experimental results for energy consumptions is 3.43%, which demonstrated the effectiveness of the proposed method. An interesting topic along this direction is to study the quantitative effect of each independent variable to the overall system performance, therefore to reduce the search space, which is currently under study and the findings will be reported later. Acknowledgements The work was funded by National Research Foundation of Singapore: NRF2008EWT-CERP002-010. The other project partners are also acknowledged.
401
Evaporator air mass flow rate penalty:
8 _ _ _ 0 < if me;air;min me;air me;air;max _ e;air;min m _ e;air;min > m _ e;air m _ e;air _ e;air if m K2 ¼ m : _ e;air;max m _ e;air > m _ e;air m _ e;air;max _ e;air m if m (A2) Condensing saturated temperature penalty:
K3 ¼
0 Tc;r;i Tc;r;sat Tc;r;sat
if Tc;r;o Tc;r;sat Tc;r;i if Tc;r;i < Tc;r;sat (A3)
Evaporating saturated temperature penalty:
K4 ¼
if Te;r;sat Te;air;i if Te;r;sat > Te;air;i
0 Te;r;sat Te;air;i Te;r;sat
Superheat penalty:
8 0 < if Te;r;sat þ 7 Te;r;o Te;r;sat þ 25 if Te;r;o < Te;r;sat þ 7 K5 ¼ Te;r;o Te;r;sat Te;r;o : Te;r;o Te;r;sat Te;r;sat if Te;r;o > Te;r;sat þ 25 (A5) Subcool penalty:
K6 ¼
0 Tc;r;o Tc;r;sat Tc;r;sat
if Tc;r;o Tc;r;sat if Tc;r;o > Tc;r;sat
(A6)
Refrigerant mass flow rate penalty:
8 0 < m _ r;min m _r m _ K7 ¼ r : _ r;max m _r m _ r m
_rm _ r;max _ r;min m if m _ r;min > m _r if m _r>m _ r;max if m
(A7)
EEV opening percentage penalty:
8 < K8
0 1 : jðAv 1Þ=Av j
0 Av 1 Av < 0 Av > 1
(A8)
Appendix B. Calculation of state penalty functions Following procedures are used to calculate penalty for individual system state variable (dependent variable): _ r and Step 1: for the assigned independent variables Pc, Pe and u, m Av are determined by Eqs. (5) and (8), where r is a linear function of Pc and Tc,r,o
r ¼ fr Pc ; Tc;r;o ¼ ar Pc þ br Tc;r;o þ cr
(B1)
and the coefficients ar, br, cr can be obtained by curve fitting for give refrigerant. Step 2: Q_ com and Q_ c are determined by Eqs. (7) and (9) Step 3: for He,r,o and Hc,r,o determined through
He;r;o ¼ fHe;r;o Pe ; Te;r;o ¼ aHe;r;o Pe þ bHe;r;o Te;r;o þ cHe;r;o
Appendix A. Definition of state penalty functions
(A4)
(B2)
and
Condenser air mass flow rate penalty:
8 _ _ _ 0 < if mc;air;min mc;air mc;air;max _ c;air;min m _ c;air;min > m _ c;air m _ c;air _ c;air if m K1 ¼ m : _ c;air;max m _ c;air > m _ c;air m _ c;air;max _ c;air m if m (A1)
Hc;r;o ¼ fHc;r;o Pc ; Tc;r;o ¼ aHc;r;o Pc þ bHc;r;o Tc;r;o þ cHc;r;o
(B3)
where the coefficients aHe;r;o , bHe;r;o , cHe;r;o , aHc;r;o , bHc;r;o and cHc;r;o can be obtained by curve fitting for the given refrigerant. Te,r,o can then be solved by
402
Te;r;o ¼
L. Zhao et al. / Energy 55 (2013) 392e402
He;r;o cHe;r;o
. bHe;r;o aHe;r;o Pe
(B4)
_ r þ He;r;i and He,r,i ¼ Hc,r,o where He;r;o ¼ Q_ e =m Step 4: for the Hc,r,i determined through
Hc;r;i ¼ fHc;r;i Pc ; Tc;r;i ¼ aHc;r;i Pc þ bHc;r;i Tc;r;i þ cHc;r;i
(B5)
where the coefficients aHc;r;i , bHc;r;i and cHc;r;i can be obtained by curve fitting for the given refrigerant. Tc,r,i can then be determined by
Tc;r;i ¼
. Hc;r;i aHc;r;i Hc;r;i cHc;r;i bHc;r;i
(B6)
_ r þ Hc;r;o where Hc;r;i ¼ Q_ c =m Step 5: for He,g, Te,r,sat, Hc,fg and Tc,r,sat expressed by
He;g ¼ fHe;g ðPe Þ ¼ aHe;g Pe2 þ bHe;g Pe þ cHe;g
(B7)
Te;r;sat ¼ fTe;r;sat ðPe Þ ¼ aTe;r;sat Pe2 þ bTe;r;sat Pe þ cTe;r;sat
(B8)
Hc;fg ¼ fHc;fg ðPc Þ ¼ aHc;fg Pc2 þ bHc;fg Pc þ cHc;fg
(B9)
and
Tc;r;sat ¼ fTc;r;sat ðPc Þ ¼ aTc;r;sat Pc2 þ bTc;r;sat Pc þ cTc;r;sat
(B10)
where the coefficients aHe;g , bHe;g , cHe;g , aTe;r;sat , bTe;r;sat , cTe;r;sat , aHc;fg , bHc;fg , cHc;fg , aTc;r;sat , bTc;r;sat , cTc;r;sat can be obtained by curve fitting for the _ e;air and m _ c;air given refrigerant and operating conditions. Then m are determined through Eqs. (1) and (3). References [1] Dossat RJ, Horan TJ. Principles of refrigeration. 5th ed. Upper Saddle River, N.J., U.S.A: Prentice Hall; 2002. [2] US household electricity report. Washington D.C, U.S.A: Energy Information Administration; 2005. [3] Green building design guide. Singapore: Building & Construction Authority; 2007. [4] Jia X, Tso CP, Jolly P, Wong YW. Distributed steady and dynamic modelling of dry-expansion evaporators. International Journal of Refrigeration 1999;22: 126e36. [5] Zhang L, Yang C, Zhou J. A distributed parameter model and its application in optimizing the plate-fin heat exchanger based on the minimum entropy generation. International Journal of Thermal Sciences 2010;49:1427e36. [6] Ren C-Q. Corrections to the simple effectiveness-NTU method for counterflow cooling towers and packed bed liquid desiccanteair contact systems. International Journal of Heat and Mass Transfer 2008;51:237e45.
[7] Arsenyeva OP, Tovazhnyansky LL, Kapustenko PO, Khavin GL. Optimal design of plate-and-frame heat exchangers for efficient heat recovery in process industries. Energy 2011;36:4588e98. [8] Metin Ertunc H, Hosoz M. Comparative analysis of an evaporative condenser using artificial neural network and adaptive neuro-fuzzy inference system. International Journal of Refrigeration 2008;31:1426e36. [9] Kim M-H, Kim J-H, Choi A-S, Jeong J-W. Experimental study on the heat exchange effectiveness of a dry coil indirect evaporation cooler under various operating conditions. Energy 2011;36:6479e89. [10] Wang Y-W, Cai W-J, Soh Y-C. A simplified modeling of cooling coils for control and optimization of HVAC systems. Energy Conversion and Management 2004;45:2915e30. [11] Xudong D, Wenjian C, Lei J. Evaporator modeling e a hybrid approach. Applied Energy 2009;86:81e8. [12] Xudong D, Wenjian C, Lei J. A hybrid condenser model for real-time applications in performance monitoring, control and optimization. Energy Conversion and Management 2009;50:1513e21. [13] Stoecker WF. Industrial refrigeration handbook. New York, U.S.A: McGrawHill; 1998. [14] Jensen JB, Skogestad S. Optimal operation of simple refrigeration cycles: part I: degrees of freedom and optimality of sub-cooling. Computers & Chemical Engineering 2007;31:712e21. [15] Jensen JB, Skogestad S. Optimal operation of simple refrigeration cycles: part II: selection of controlled variables. Computers & Chemical Engineering 2007;31:1590e601. [16] Jensen JB, Skogestad S. Steady-state operational degrees of freedom with application to refrigeration cycles. Industrial and Engineering Chemistry Research 2009;48:6652e9. [17] Larsen LFS, Thybo C, Stoustrup J, Rasmussen H, A method for online steady state energy minimization, with application to refrigeration systems. In: 2004 43rd IEEE Conference on Decision and Control (CDC), December 14, 2004e December 17, 2004, Nassau, Bahamas, 2004, pp. 4708e4713. [18] Selbas¸ R, Kızılkan Ö, S¸encan A. Thermoeconomic optimization of subcooled and superheated vapor compression refrigeration cycle. Energy 2006;31: 2108e28. [19] Yuan F, Chen Q. A global optimization method for evaporative cooling systems based on the entransy theory. Energy 2012;42:181e91. [20] Fong KF, Lee CK, Chow CK, Yuen SY. Simulationeoptimization of solarethermal refrigeration systems for office use in subtropical Hong Kong. Energy 2011;36:6298e307. [21] Kusiak A, Tang F, Xu G. Multi-objective optimization of HVAC system with an evolutionary computation algorithm. Energy 2011;36:2440e9. [22] Kusiak A, Xu G, Tang F. Optimization of an HVAC system with a strength multi-objective particle-swarm algorithm. Energy 2011;36:5935e43. [23] Kusiak A, Xu G. Modeling and optimization of HVAC systems using a dynamic neural network. Energy 2012;42:241e50. [24] Hovgaard TG, Larsen LFS, Edlund K, Jørgensen JB. Model predictive control technologies for efficient and flexible power consumption in refrigeration systems. Energy 2012;44:105e16. [25] Arora CP. Refrigeration and air conditioning. 2nd ed. New Delhi, India: Tata McGraw-Hill; 2000. [26] Ding X, Jia L, Cai W, A hybrid modeling for the real-time control and optimization of compressors. In: 2009 4th IEEE Conference on Industrial Electronics and Applications, ICIEA 2009, May 25, 2009eMay 27, 2009, Xi’an, China, 2009, pp. 3256e3261. [27] Dixon SL. Fluid mechanics and thermodynamics of turbomachinery. 4th ed. Boston, U.S.A: Butterworth-Heinemann; 1998. [28] Mithraratne P, Wijeysundera NE. An experimental and numerical study of hunting in thermostatic-expansion-valve-controlled evaporators. International Journal of Refrigeration 2002;25:992e8. [29] Goldberg DE. Genetic algorithms in search, optimization, and machine learning. Reading, Mass, U.S.A: Addison-Wesley Pub. Co.; 1989.