- Email: [email protected]
Optimal control of indoor-air cooling in buildings using a reduced order model
Energy 116 (2016) 1191e1204
Contents lists available at ScienceDirect
Energy journal homepage: www.elsevier.com/locate/energy
Optimal control of indoor-air cooling in buildings using a reduced order model S. Ben Ayed a, *, D. Kim b, J.T. Borggaard a, E.M. Cliff a a b
Interdisciplinary Center for Applied Mathematics, Virginia Polytechnic Institute and State University, Blacksburg, VA, USA School of Mechanical Engineering, Purdue University, West Lafayette, IN, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 30 November 2015 Received in revised form 15 August 2016 Accepted 7 October 2016 Available online 27 October 2016
In this work, we use Computational Fluid Dynamics (CFD) to generate the distributed dynamic responses of temperature and moisture in a restaurant that correspond to perturbations of input variables. The CFD model is validated with experimental data of the roof top units return temperatures. A Reduced Order Model (ROM) is developed by approximating the responses to these perturbations using a linear timeinvariant model. The resulting indoor-air model is coupled to a dynamic envelope model with longer time scales. Assuming the outside temperature distribution and the occupants' loads are known throughout the day, optimal control is applied to our coupled model to minimize the cooling power subject to local and global comfort constraints. The results show that the applied method is more efficient and comfortable than the results given by the experimentally measured response of a conventional set-point control strategy. Published by Elsevier Ltd.
Keywords: Optimal control Building energy Reduced order model Berkeley library for optimization modeling Envelope model Computational fluid dynamics
1. Introduction According to The U.S. Energy Information administration, in 2013, 40% of total energy consumption in the U.S. was consumed in residential and commercial buildings. This generally consists of heating-ventilation-air-conditioning (HVAC), lighting, and plug loads. Thus, even a modest reduction in building energy use represents significant contributions toward the relief of a host of pressing environmental issues. Here we focus on cooling requirements and specifically on the development of models suitable for study of predictive control of the indoor environment. Optimal control of the HVAC energy in buildings has been the object of various research works during the last decades. Naidu and Rieger wrote an extensive review about advanced control strategies in HVAC systems ([1,2]). Sane et al. [3] presented a survey of problems in controls and optimization related to commercial chilled water cooling plants. They included in multivariable optimal control, dynamic resource allocation and control over networks. Several works have explicitly reported the savings due to the use of optimal control for HVAC systems. House and Smith [4] used optimal control of the supply air temperature to minimize the cost
* Corresponding author. E-mail address: [email protected] (S. Ben Ayed). http://dx.doi.org/10.1016/j.energy.2016.10.022 0360-5442/Published by Elsevier Ltd.
for energy consumption. Comparing this method to the conventional control, they found a cost savings of 11%. Yu and Chan ([5]) used an optimal control strategy to control the cooling towers and condenser water pumps in a water-cooled chiller system. They concluded that their procedure contributed to annual energy savings of 8.6% and cost savings of 9.9%. For a good estimate of the optimal energy cost, thermal load modeling is very important. Various models and computational techniques have been adopted. Oldewurtel et al. [6] claims that computer-aided modeling tools like TRNSYS [7] and EnergyPlus [8] are mainly designed to estimate the energy use of building and cannot be readily used for control. They suggested the use of a thermal Resistance-Capacitance (RC) network for control purposes. Many others [9e11] used the same simple RC model to estimate the thermal load. However, Ma et al. [12] used a black box model structure instead of an RC model. They used EnergyPlus to model a commercial building, generated outputs for a nominal configuration, used the input/output data to identify mathematical models for power and temperature, then generated the optimal control. Zlatanovi [13] used experimental data to generate the energy saving estimation model. This was based on monitoring the energy consumption of each HVAC system appliance in a large trade center on a daily basis for one year. Modeling the indoor environment is very important in small to
1192
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
medium commercial buildings that utilize rooftop units (RTUs) to condition large open spaces. Examples include retail stores, restaurants and factories. In fact, for these buildings, the well-mixed assumption is not practical due to non-uniformly distributed diffuser systems, local cooling/heating actions corresponding to each RTU, and non-uniform heat gains. In the present work, we develop a coupled CFD-envelope model having two distinct time scales for a specific restaurant. The CFD model captures important dynamics of the spatial distributions of temperature and moisture inside the building; while the envelope model studies possible inside/outside heat transfer mechanisms, namely convection, conduction and radiation. Various approaches to couple prediction codes for fluid dynamics and solid conduction in a problem of conjugate heat transfer have been studied [14,15]. We use optimal control with this coupled model to minimize the cooling energy in a given restaurant. Since the optimization procedure is iterative, a reduced order model that captures the essential dynamics of the problem and has a low computational cost [16,18,19] is developed for performance analysis of optimal control for open-spaces by utilizing ROM. For this study, the reduced order model is created from CFD simulations such that the steady-state responses are adequately captured and such that the resulting model easily couples to the envelope model. 2. CFD of the restaurant indoor air and reduced order model Our approach to developing a reduced-order model is to construct a linear time-invariant (LTI) model that approximates the input-output behavior observed in computational fluid dynamic (CFD) simulations. Accordingly, a CFD model for the restaurant was developed using FLUENT. The coupling between the indoor-air and the envelope models is naturally achieved by treating the bounding surface temperatures as inputs to the CFD and the corresponding surface heat fluxes as outputs; that is, we use CFD to approximate the Dirichlet-to-Neumann (dynamic) response. Other CFD inputs include the mass-flow-rate, temperature, and humidity supplied by each of the roof-top-units (RTU), the sensible and latent loads provided by occupants, and external loads generated by the kitchen activity. Additional outputs of interest include air temperature/ humidity at selected zones/locations, and return air properties. The mesh is formed of 4.70 m tetrahedral cells with 624 k vertices. We use a refinement of about 1/3 around the diffuser and returns and approximately 1/7 in the wall regions where the thermostat temperatures are averaged (around 250 tetrahedral faces). The moisture is treated by integrating the species transport model in FLUENT. Hence the fluid is a mixture of air and water vapor. This FLUENT feature models the mixing and transport of species by solving conservation equations describing convection and diffusion for the species. The turbulence model used in the work is a standard kEpsilon model with standard wall functions. In fact, it is a simple two-equation model in which the solution of two separate transport equations allows the turbulent velocity and length scales to be independently determined. Robustness, economy, and reasonable accuracy for a wide range of turbulent flows explain this model's popularity in industrial flow and heat transfer simulations [17].
North Wall: 27.00 C East Wall: 25.58 C West Wall: 28.43 C Floor: 24.92 C Ceiling: 25.48 C
2.1.2. Kitchen interface The south wall boundary of our computational domain (see Fig. 1) separates the dining zone from the kitchen, which is not modeled. For the solid part of this boundary we impose no-slip flow and adiabatic thermal conditions. For the opening(s) between the kitchen and the dining zones (red) we have prescribed the flow from the kitchen. mass-flow-rate: 0.591 kg s temperature: 25.00 C water mass-fraction: 16 g/kg
2.1.3. Supply air Each supply air vent is a 2 ft by 2 ft ceiling mounted square. The air supply from each RTU is uniformly distributed among the associated vents, and is supplied at 16.60 C. The velocity at each vent is normal (i.e. downward), and is spatially distributed using a flattened parabolic profile. RTU#1- mass-flow-rate: 2.9320 kg/s; water mass-fraction kg; distributed over 19 supply vents (blue) RTU#2- mass-flow-rate: 0.5320 kg/s; water mass-fraction kg; distributed over 04 supply vents (red) RTU#3- mass-flow-rate: 0.3990 kg/s; water mass-fraction kg; distributed over 03 supply vents (purple) RTU#4- mass-flow-rate: 0.5320 kg/s; water mass-fraction kg; distributed over 04 supply vents (green)
7 g/ 9 g/ 9 g/ 9 g/
Note that RTU#3 supplies four vents but only three of these are in the computational domain. 2.1.4. Return air Each of the return vents is ceiling mounted. These are specified as pressure outlet with the flow normal to the boundary, and the pressure at 40 Pa below that on the floor. The pressure difference approximates the static head loss under vertical equilibrium (i.e. DP¼r g Dh) at a ceiling height of 3.4 m.
2.1. Nominal boundary/source conditions The indoor-air response is decomposed into a nominally steady part and a perturbation. For the steady part we apply fixed boundary and source conditions enumerated here: 2.1.1. Bounding surfaces For five of the bounding surfaces we prescribe no-slip flow and fixed temperature (Dirichlet) data
Fig. 1. Grid with zones - the left (green) wall is North. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
2.1.5. Dining area source The dining area is the roughly 69 m3 volume shown in yellow in Fig. 1. The intent is to simulate the sensible/latent loads produced by the occupants. In this volume we uniformly distribute an energy source of approximately 8,750 W, and a moisture source of approximately 3.08 g/s. 2.1.6. Wine bar area source The winebar area is the roughly 23 m3 volume that surrounds the winebar in the north end of the domain. The intent is to simulate the sensible/latent loads produced by the occupants. In this volume we uniformly distribute an energy source of approximately 3,310 W, and a moisture source of approximately 1.10 g/s. 2.2. Steady state response The system was subjected to fixed boundary and source conditions chosen to represent a mean condition in the space (see Section 2.1). FLUENT was run in both steady and unsteady modes; the latter consisted of 7 h of simulated time with external (boundary & source) conditions fixed. Some responses for the last hour are shown in Fig. 2. The black horizontal lines are averages over the 1-h window. The meaning of these various temperatures can be understood in terms of Fig. 1. In particular dining zone is the yellow rectangular volume in the center of the figure. wine_bar zone is the yellow volume in the shape of a reflected D in the left of the figure. fireplace is the surface of the red box in the left foreground (to the left of the entrance bump-in).
1193
that a steady flow exists, moreover the computed results depend on the grid which is likely somewhat too coarse.
2.3. FLUENT model validation In the early morning of 21 August 2013, a test to validate the FLUENT results was conducted at the Harvest Grill site. At approximately 03:30 the supply air temperature from RTU1 was reduced from about 16.5 C to about 9.5 C as shown in Fig. 3. A FLUENT run was made wherein a piecewise linear approximation of this temperature profile was prescribed for the RTU1 supply air while (most) other source/boundary conditions were maintained at nominal values. The only exception was that the temperature of the (supply) kitchen air was changed from the nominal 25 Ce21 C. The change was made because at the time of the experiment there was no cooking activity in the kitchen. The chosen value (21 C) was the ambient temperature in the RTN3 return air which is closely coupled to the kitchen air. Initial data for the simulation was a result of previous steady state runs. Fig. 4 displays the comparisons of the FLUENT and measured data for the return air temperatures in RTN1 and RTN2. Note that these start at rather different values; that is, the initial condition in the FLUENT simulation does not well-represent the overnight conditions of the test. To compare the trends we shifted the FLUENT predictions so that the initial values match. Fig. 5 displays the shifted comparisons for all four return air locations. In the time interval [03:30, 05:00] the FLUENT predictions are very close to the data. After 05:00 the supply air temperatures for the other RTU's began to change; a feature not captured in the FLUENT simulation.
2.4. Input variables
Note that the room average is the least oscillatory: it is the massaveraged temperature over the entire room. The dining and wine_bar values are mass-averaged over the cited volumes. The fireplace value is averaged over the external surfaces of the designated red-box, and exhibits the largest oscillations. Also, note that the values at the right end (t ¼ 0) do not generally match the average, and this is particularly noticeable in the fireplace data. One might expect that the oscillations would decay and the internal flow in the room would become steady. In fact there is no guarantee
Twenty-two of the boundary conditions discussed in Section 2.1 were identified as useful for injecting control or disturbance to the system. These are enumerated in Table 1. For each input we ran a FLUENT simulation beginning from the final state of the steady simulation and proceeding for 7200 s. In each run, one of the inputs in Table 1 was subjected to a ramp input beginning at t ¼ 8 s, ramping to the appropriate perturbation value over a 60 s span, then remaining constant at the perturbed value.
Fig. 2. Steady response; Zone Average Temperatures.
Fig. 3. RTU1 supply air temperature.
1194
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
Fig. 4. RTN1 (left), RTN2 (right).
Fig. 5. Shifted comparisons: RTN1 (upper left), RTN2 (upper right) RTN3 (lower left), RTN4 (lower right).
2.5. LTI model Our approach is to construct an LTI model that approximates the perturbed CFD responses generated by FLUENT simulations. Whereas there are a number of approaches to extract reduced-order LTI models from response data (for example [18,19]), we are motivated to find a low order model with the appropriate asymptotic behavior in a computationally efficient way. Nominally Steady Outputs In order to describe perturbed responses in a meaningful way, we first deal with defining steady values of the output variables. The output variables are listed in Table 2. Outputs 1 to 5 represent the average thermal flux through the walls which will be coupled to a building envelope model as shown in Fig. 10. Outputs 6 to 9, 10 to 13 and 14 to 17 represent the return air conditions, i.e. the average of the mass flow rates,
temperatures and humidities respectively. These are used to calculate HVAC power consumption. Outputs 18 to 21 indicates local temperatures at thermostat locations. In order to measure comfort, air conditions at some local area of interests and zonal volume averaged temperatures and humidities were chosen as outputs. Outputs 22 to 33 are the conditions at the local points while the rests are the zonal conditions. The final two outputs 42 and 43 are left for the average temperature and moisture of the whole space. It is clear from Fig. 2 that notable oscillations persist in the output variables. There are several potential explanations for this observed behavior: 1. There is no guarantee that the Reynolds-Averaged Navier-Stokes PDE model has a steady solution;
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
1195
samples). For each input, the perturbed output values were computed as:
Table 1 Input variables: Nominal and perturbation values. Number
Name
Nominal Value
Increment
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
north_temp east_temp west_temp floor_temp ceiling_temp rtu_1_temp rtu_2_temp rtu_3_temp rtu_4_temp rtu_1_h2o rtu_2_h2o rtu_3_h2o rtu_4_h2o rtu_1_mass rtu_2_mass rtu_3_mass rtu_4_mass load_dining load_winebar kitchen_inflow_temp kitchen_inflow_h20 kitchen_inflow_mass
27.00 C 25.58 C 28.43 C 24.92 C 25.46 C 16.60 C 16.60 C 16.60 C 16.60 C 7.00 g/kg 9.00 g/kg 9.00 g/kg 9.00 g/kg 2.9320 kg/s 0.5320 kg/s 0.3990 kg/s 0.5320 kg/s 100% 100% 25.00 C 16.00 g/kg 0.5910 kg/s
5.00 C 5.00 C 5.00 C 5.00 C 5.00 C 5.00 C 5.00 C 5.00 C 5.00 C 1.40 g/kg 1.80 g/kg 1.80 g/kg 1.80 g/kg 0.5864 kg/s 0.1064 kg/s 0.0798 kg/s 0.1064 kg/s 20% 20% 5.00 C 3.20 g/kg 0.1180 kg/s
ypert tJ ¼ Ydata tJ Ynom
J ¼ 1; …; 1800:
Best Fit There are several elegant approaches to extract a model from time-response data ([19,20]). Our approach decomposes the problem into sub-problems by treating each input response separately. That is, we hypothesize a structure for a low-order model and formulate an optimization problem to find appropriate parameter values. Specifically, for each input we propose a firstorder model of the form
1 _ ¼ xðtÞ þ b uðtÞ xðtÞ
(1)
yðtÞ ¼ C xðtÞ þ D uðtÞ;
(2)
t
where t is a time constant, b is a scalar, and C; D2IR43 . We take b ¼ 0.01 (scaling). Most of the 2243 input/output pairs are represented as strictly-proper transfer functions (D ¼ 0). However, for the first five inputs (surface temperatures) the heat flux output for
Table 2 Selected nominal steady outputs. Number
Name
Value
Number
Name
Value
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
q_north q_east q_west q_floor q_ceiling rtn_1_mass rtn_2_mass rtn_3_mass rtn_4_mass rtn_1_temp rtn_2_temp rtn_3_temp rtn_4_temp rtn_1_h2o rtn_2_h2o rtn_ 3_h2o rtn_4_h2o thermo_1_temp thermo_2_temp thermo_3_temp thermo_4_temp dk_2_temp
11.5095 w/m2 10.4278 w/m2 19.3538 w/m2 15.8151 w/m2 4.9463 w/m2 2.4959 kg/s 0.8248 kg/s 0.8248 kg/s 0.8392 kg/s 22.1082 C 22.3562 C 21.7175 C 22.6413 C 9.3741 g/kg 9.2538 g/kg 9.2147 g/kg 10.5528 g/kg 21.5048 C 21.80948 C 21.7650 C 21.8446 C 21.7961 C
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
dk_4_temp dk_2_h2o dk_4_h2o rtu1_z1_temp rtu2_z1_temp rtu3_z1_temp rtu4_z1_temp rtu1_z1_h2o rtu2_z1_h2o rtu3_z1_h2o rtu4_z1_h2o dining_room_temp dining_room_h2o winebar_temp winebar_room_h2o winebar_skirt_temp winebar_skirt_h2o fireplace_temp fireplace_h2o room_avg_temp room_avg_h2o
21.7974 C 9.9237 g/kg 9.8813 g/kg 21.1414 C 22.2903 C 21.3036 C 22.0641 C 8.8190 g/kg 9.903 g/kg 9.4225 g/kg 9.9617 g/kg 21.9751 C 9.2293 g/kg 22.3145 C 10.0578 g/kg 21.65525 C 9.8574 g/kg 21.7940 C 9.6071 g/kg 22.0661 C 9.4163 g/kg
2. The grid may not be sufficiently resolved to represent the PDE behavior; and/or, 3. The iterative procedure for the nonlinear algebraic system may not be sufficiently converged. Whatever the cause(s) we simply define steady values for the outputs by averaging over a time-window. In particular, steady output values are computed by averaging the observed outputs over the final 60 min of the nominally steady solution. Steady values for the four output variables in Fig. 2 are indicated by solid lines. Steady values for all of the output variables are listed in Table 2. CFD Perturbations As noted above, in each run, one of the inputs in Table 1 was subjected to a ramp input beginning at t ¼ 8 s, ramping linearly to the corresponding perturbation value over 60 s, then remaining constant at the perturbed value. In each case the output variables (see Table 2 were recorded at 4 s intervals (1800
that surface has a single non-zero entry in D, since the surface temperature is directly coupled to the heat flux on that surface. Thus, we expect that for each surface heat flux output (the first five) we have DJ ¼ hdı J , whereas for all other outputs we have D ¼ 0. The time-scales of the building envelope are typically much longer than those of the indoor-air. Therefore, during most of the response the indoor-air system will be near equilibrium. For this reason we require that the ROM accurately represent the indoor-air dynamics for responses near equilibrium. This requirement is not naturally imposed in other model reduction methods and leads to our ad-hoc approach. Since the input perturbation function u goes to a constant value (say, uss), then with t > 0 the dynamics (1) will go to the stable equilibrium
xss ¼ tbuss 0yss ¼ ðC t b þ DÞuss : Solving for C we find
1196
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
Cðt; DÞ ¼
yss =uss D ; tb
(3)
where the vector D2IR43 has (at most) a single non-zero entry. To implement (3) we must specify a method for computing a steadystate value of the perturbed output, i.e. yss 2IR43 . Because the output data includes residual oscillations it is prudent to use some form of averaging. For the results presented here we averaged y(tı) over the last 120 samples (out of 1800 total), that is ı2f1681; 1682; …1800g. This amounts to an 8 min wide sampling window. With these preliminaries we define a scalar-valued error function
Jðt; DÞ ¼
1800 Xh
i2 ypert ðtı Þ yðtı ; t; DÞ ;
(4)
ı¼1
where yðtı ; t; DÞ is computed from (1)e(2) with the state-to-output map given by (3), and the input (u) is the step-like input described in Section 2.5. For each input we compute optimal values of t and (possibly) D by minimizing the error functional J in (4). These calculations were done using MATLAB's fminsearch procedure. The five non-zero values for the D arrays are; 1. 2. 3. 4. 5.
q_north: 2.912 W/m2/ C q_east: 2.038 W/m2/ C q_west: 3.514 W/m2/ C q_floor: 4.281 W/m2/ C q_ceiling: 1.423 W/m2/ C
Note that these can be interpreted as (time and spatial) averages of the film cooling coefficient on the respective surfaces.
2.6. ROM-FLUENT comparisons Fig. 6 provides comparisons of the FLUENT and the reduced-ordermodel (ROM) responses to a change in the east wall temperature. The input perturbation ramps to þ5 C over 60 s and the corresponding ramp in the East wall heat flux is seen in the upper graph. In the lower graph the ceiling heat flux is seen to decrease by 0.6 W/ m2, presumably the result of an increase in the air temperature engendered by the increased wall temperature. Note that the graphs display only the first 15 min of the response; the data fits are based on 2 h of data. Similarly, Fig. 7 provides comparisons of the FLUENT and the reduced-order-model (ROM) responses to a 5 C decrease in the RTU#1 supply air temperature. The upper graph displays the change in return 1 air temperature, whereas the lower graph characterizes the temperature in the dining zone. Fig. 8 provides also comparisons of the FLUENT and the reducedorder-model (ROM) responses to a 3.2 g/kg increase in the kitchen inflow moisture. The upper graph displays the change in return 1 water vapor mass fraction, whereas the lower graph shows the response of the dining zone water vapor mass fraction.
3. Development of a thermal building envelope model In order to connect a reduced-order indoor air model to a building envelope model, a general formula derived in Ref. [21] is adopted. The majority of the derivations in this section is the same as [21], but it is slightly modified for the coupling.
Fig. 6. Comparison of Fluent and ROM T_east responses.
3.1. Conduction through walls A finite volume formulation is used to describe the heat conduction through walls and is depicted in Fig. 9. For any jth node in a wall except the first and last nodes, an energy balance leads to
rij Cji wij
dTji dt
i i i i i i i ¼ hLcd Tj1 hLcd þ hRcd Tji þ hRcd Tjþ1 þ qgen ; j
j
j
j
j
(5) where rij is the density at jth node in ith wall [kg/m3], Cji is the specific heat capacity of jth node in ith wall [J/kgK] and wij is the i i kL j width of control volume of jth node in ith wall [m]. hC DL ¼ ji j wL j j i where kL is the thermal conductivity at left surface of the jth node j i in ith wall [W/m2K] and wL is the distance from the (j1)th node j i to the jth node in ith wall [m].qgen is an energy source [W/m2] j
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
Fig. 8. Comparison of Fluent and ROM h_kitchen responses. Fig. 7. Comparison of Fluent and ROM T_rtu_1 responses.
inside the jth finite control volume that belongs to the ith wall. 2 q 1 3 gen j 2 7 6 7 ! 6 Using the following matrix notation, q j ≡6 qgen 7, j 4 :: 5 Nw qgen j
2
Tj1
3
6 2 7 ! 6 Tj 7 7; T j ≡6 6 7 4 :: 5 Nw Tj
2 L=R 1 hcd j 6 6 ~ L=R ≡6 H 6 j 6 4
3
0 L=R 2 hcd j
::
0 1 2 rCw j 6 6 ~ ≡6 C j 6 4 0
L=R Nw hcd
7 7 7 7; and 7 5
j
3
0 2 rCw j
::
Nw
rCw
7 7 7; 7 5
j
Eq. (5) can be expressed in the following matrix form
Fig. 9. Notation for conduction through walls.
1197
1198
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
3 3 i i qi1 ¼ hicv;ex Ta þ 4sεi1 Fsky T sky Tsky þ Fgrd T grd Tgrd þ ai1 qiSWR ; T sky ¼
Tsky þ T1 ; 2
T grd ¼
Tgrd þ T1 : 2
The mean temperatures for long-wave exchange between the surface and sky and surface and ground are assumed to be the same for all outside surfaces. Eq. (8) can be generalized and written in compact matrix form as
! ! R ! ~ dT 1 ¼ H ~ ~ ~R ! ~ C 1 cd;1 þ H rad;ex þ H cv;ex T 1 þ H cd;1 T 2 þ q 1 ; dt
(9)
where
3 3 ! i i ð q 1 Þi ¼ hicv;ex Ta þ 4sεi1 Fsky T sky Tsky þ Fgrd T grd Tgrd
Fig. 10. Simplified scheme of ROM-envelope coupling.
þ ai1 qiSWR : ! L ! ~ dT j ¼ H ~L ! ~R ! ~R ! ~ C j cd;j T j1 H cd;j þ H cd;j T j þ H cd;j T jþ1 þ q j : dt
(6)
! Note that T j is a group of all of the wall temperature nodes belonging to an iTindividual zone. More precisely, h ! T j ¼ Tj1 Tj2 :: TjNw whereNw is the number of walls in a zone and the superscript T represents the transpose.
3.3. Heat balance at the interior surfaces For the nth zone and ith wall, the energy balance equation for the inside surface is
rin Cni win 3.2. Heat balance at the exterior surfaces For any ith outside wall (connected to the external ambient) belonging to an individual zone, the heat balance equation at the surface is
dT i ri1 C1i wi1 1 dt
¼
hicv;ex
Ta
T1i
þ
i hRcd T2i 1
T1i
þ
(7)
where T1i represents the wall temperature of the first node which is set to be an outside surface of the wall. With the assumptions that the outside surface is gray and diffuse and the air is a nonparticipating radiation media, net long wavelength interactions with the environment can be expressed as i qiLWR ¼ sεi1 Fsky
Tsky
4
4 4 4 i þ sεi1 Fgrd Tgrd T1i : T1i
Using a linear approximation of the long-wave heat exchange term gives
ri1 C1i wi1
i i dT1i ¼ hicv;ex þ hRcd þ hirad;ex T1i þ hRcd T2i þ qi1 ; 1 1 dt
where
hirad;ex
3 3 i i ; ¼ 4sεi1 Fsky T sky þ Fgrd T grd
(8)
(10)
where qinet;rad is net radiative flux out of the inside wall and qicv;in the net heat flux from the zone air and allows for the CFD coupling. The radiosity method is utilized to express the net flux under the assumption that the walls are gray, diffuse and opaque. The same linearization method used in Eq. (8) is employed leading to
h 0! i ! ~ 1 B ~ Tn! q net;rad ¼ A ho
ai1 qiSWR
þ qiLWR ;
i dTni i ¼ hLcd Tn1 Tni qinet;rad þ qicv;in ; dt n
(11)
where
~ ≡dij rj F ; A ij εj εj ij and 0 ~ ¼ 4s d F T 3 : B ij ij ij
Radiosity does not appear explicitly in the above expression, which is convenient for building simulations. Since hi0 represents an external radiative source acting on the ith surface, the effects of internal sources and transmitted solar energy though windows are treated in a consistent manner. For any shaped room, the net radiative flux can be explicitly calculated as a function of surface temperatures if the view factors and the external radiative sources ! ~ 1 ! ~ 1 ~ 0 ~ are known. By letting H h o, rad;in ≡A B and q n ≡A
! ! L ! ! ~nd T n ¼ H ~L ! ~ ~ C cd;n T n1 H cd;n þ H rad;in T n þ q n þ q n;cv;in : dt (12)
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
3.4. State-space representation of the thermal envelope network Gathering the system of equations that represent heat balance equations for external to internal wall elements, i.e. Eqs. (6), (9) and (12), _ ! ! ~w! ~ ww ! Tw ¼ H T w þ q w þ q cv;in ; C
(13)
where
2 ~w C
6 ¼6 4 2
~ C 1
3 ~ C 2
~ H 1 6 6 ~L H ¼6 6 cd;2 4
::
~n C
7 7; 5
~ ww H
~R H L cd;1 R ~ ~ H cd;2 þ H cd;2 ::
Table 3 Model variables. Variable Tc i ði ¼ 1 : 4Þ m_ s i ði ¼ 1 : 4Þ fs i ði ¼ 1 : 4Þ DL WBL KT m_ r i ði ¼ 1 : 4Þ Tr i ði ¼ 1 : 4Þ Tz i ði ¼ 1 : 4Þ fr ði ¼ 1 : 4Þ
Significance
ROM Input/Output
RTU control temperature
IN(6)to IN(9)
Ta di
i
3 ~R H cd;2 :: ~L H cd;n
li
7 7 7; 7 :: 5 ~ H n
~R þ H ~ cv;ex þ H ~ ~ n ¼ ðH ~L þ H ~ ~ ¼ ðH H H 1 rad;ex Þ, rad;in Þ, cd;1 cd;n h h i iT ! T T T T ! T T T ! ! ! ! ! Tw¼ T and q w ¼ ! q 1 q 2 :: q n . As T 2 :: T n 1 expected, a tri-diagonal block matrix is formed with parameters that characterize heat transfer due to conduction in the walls (from node number 2 to node number n-1) and radiative/convective heat ! transfer at the boundaries (first and last nodes only). The terms q j where j ranges from 2 to n1 vanish if there are no heat flux sources inside the wall such as embedded radiant heating or ! cooling systems. The term, q cv;in , represents the set of net convective heat fluxes to the zone air and allows CFD coupling. All entries in the vector are zeros except for inner walls. 3.5. Coupling with the restaurant CFD model The one dimensional heat balance equations, Eq. (13), are coupled to a reduced-order indoor air model. Although it is possible to subdivide wall surfaces sufficiently to model temperature variations across the walls and solve three dimension heat equations, we chose the 1-D form for computational simplicity. However Zhai and Chen [22] showed that the isothermal wall dynamics coupled to a CFD model leads to reasonable accuracy in predicting spatial temperature distribution of zone air and cooling/heating load. The state-space modeling (Eq. (13)) is applied to the restaurant walls. It results in the following system of equations:
X_ e ¼ Ae Xe þ Be u Ye ¼ Ce Xe þ De u;
(14)
1199
RTU supply mass flow rate
IN(14)to IN(17)
Supply water mass fraction
IN(10)to IN(13)
Dining load Wine bar load Kitchen temperature RTU return mass flow rate
IN(18) IN(19) IN(20) OUT(6)to OUT(9)
Return temperature
OUT(10) to OUT(13)
Local zone temperature
OUT(26)to OUT(29)
Return water mass fraction
OUT(14)to OUT(17)
Average zone temperature Fraction of outside air Fraction of recirculated air
OUT(42) NA NA
values as specified in Section 2. The supply air temperatures, occupancy and cooking loads, walls, ceiling and floor temperatures coming from the envelope model are considered as perturbations to the nominal configurations. Thus, due to the linearity of the model, their effects on the ROM are additive. The control variables are the four RTU supply air temperatures and the four RTU outside air fractions. All the variables and problem parameters are summarized in Tables 3 and 4. The details of loads and occupancy are summarized as follows: Kitchen temperature (KT): 24 C rising to 28 C between t ¼ 10 h and t¼22 h. Dining load (DL): 0.2 rising to 1.2 between t ¼ 12 h and t¼22 h. Wine bar load (WBL): 0.2 rising to 1.2 between t ¼ 12 h and t¼22 h. Occupancy: the occupancy function is 1 between hmin and hmax, and 0 otherwise
4.2. Cost function (equipment model) rates, namely, In this analysis, we distinguish several mass flow m_ s i as the supply mass flow rate of each RTU, m_ oa i as the mass flow rate supplied from outside for each RTU and m_ r i as the return mass flow rate to each RTU, where the subscript i designates quantities related to the ith RTU. Continuity of these rates gives the following balance:
m_ s i ¼ m_ oa i þ li m_ r i :
(15)
Where li is the fraction of the return air supplied in each RTU. We denote by di the fraction of the supply air from outside. This fraction satisfies the following equation:
! where Xe is the state vector ( T w in Eq. (13)), u is the input vector of 22 elements comprising the outdoor weather variables and the fluxes through the walls computed from the ROM. The finite difference scheme with 2e4 nodes per material layer of a wall results in 35 elements in Xe. The number was selected after checking the grid independence of a simulation result. The output Ye represents a vector of the computed wall temperatures that are inputs to the ROM. This coupling is illustrated in a simplified scheme in Fig. 10.
Combining Eqs. (15) and (16), we get the following relation between di and li, where di is a control variable and li is computed as:
4. Optimal control of the cooling energy
li ¼ ð1 di Þ
4.1. Problem formulation
Since li is supposed to be less than or equal to unity, we distinguish two cases:
This work studies the restaurant indoor-air model for use in performing optimal control calculations. Unless modified, the default values of all the input variables are fixed at the nominal
m_ oa i ¼ di m_ s i :
m_ s : m_ r i
_ r if di 1 m , then li ¼ 1 m_ s i m_ r if di > 1 m_ s i , then li is computed using Eq. (17)
(16)
(17)
1200
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204 Table 4 Model parameters. Parameter Ta
Significance
Value
Upper bound for comfort
23 C
hmin,max Cp dmin,max Tmin,max copi(i ¼ 1:4) capmax i ði ¼ 1 : 4Þ foa dt L
Occupancy bounds Heat capacity Air fraction bounds Bounds for control temp Total coefficient of performance Total cooling capacity of each RTU
11 h, 22 h 1 kJ/(kg K) 0.1, 1.0 10 C, 30 C 3.5, 2.6, 2.6, 2.6 53, 14, 10.5, 14 kW
Outside water mass fraction Time step Cooling power for removal of 1 g of water
17 g/kg 7.5 min 2.26 kW
max
The total cooling load at each RTU is divided into sensible Qs i and latent Ql i loads. The sensible load related to the temperature change is expressed as:
Qs i ¼ Cp m_ s i Tm i Tc i
(18)
0:05 Tc i Troof ; Tc i ¼ Tc i 1:004
(19)
where Cp is the heat capacity of the air at constant temperature, Tm i is the RTU inlet air temperature, Tc i is the supply RTU tem perature, and Tc i is the modified supply RTU temperature obtained after the circulation of the supply air in the attic. The RTU inlet air temperature or the mixing temperature can be written as:
m_ r i li Tr i þ m_ oa i Toa Tm i ¼ ; m_ s i
(20)
where Toa is the outside air temperature and Tr i is the temperature at each return. The moisture formulation included in the CFD model allows us to easily compute the latent load related to the change of state, which can be expressed as:
Ql i ¼ Lðm_ r i li fr i þ m_ oa i foa fs m_ s i Þ;
(21)
where fr i , foa and fs i are the mass fraction of water respectively for the return, outside air and supply flows, and L is the amount of cooling required to remove 1 g of water. Thus, the total power expression is:
Qs i þ Ql i Pt i ¼ ; copi
(22)
2. Set the upper comfort bound for the local temperature of the different zones at the occupancy time: for hmin t hmax, Tz i Ta max . 3. Set the upper comfort bound for the average zone temperature at the occupancy time. For hmin t hmax, Ta Tajmax . 4. Cooling coil can only decrease the temperature. Tc i Tm i . 5. Set the upper and lower bounds for the outside air fraction. 0.1 di 1. 6. Set the upper and lower bounds for the recirculated air fraction. 0li 1. 7. Set the upper and lower bounds for the cooling power of each RTU.
0 Qs i þ Ql i capmax i 8. Latent heat is only removed during the occupancy period.
4.4. Numerical approximation In optimal control, one is concerned with minimizing the functional (23) from a given initial state over a forward planning horizon ([0,T]). Here we briefly study a discretized version of the problem using the MATLAB/SIMULINK based tool: Berkeley Library for Optimization Modeling. At this time BLOM 2.0 is restricted to discrete-time dynamics, so it is necessary to convert from a continuous to a discrete time model. Thus, for example, the envelope model (15) is replaced by
Xe kþ1 ¼ Fe Xe k þ Ge ;
(24)
where where copi is the coefficient of performance of RTU #i. We note that the RTUs are constant fan speed regardless of RTU stages. Therefore, fan powers are assumed to be constant and do not appear in the objective function. The final cost expression can be written as:
C¼
ZT X 4 0
Pt i dt
Fe ¼ expðAe dtÞ Ge ¼
(23)
i¼1
(25)
Zdt expðAe ðdt sÞÞBe ds
(26)
0
and dt is the time step.
where T is a 24-hour period. 5. Results 4.3. Constraints Most of the constraints are applied for each RTU i (i ¼ 1, 2, 3, 4): 1. Specify the upper and lower bounds for the control temperature for different RTUs. 10 C Tcji 30 C.
5.1. Optimal control results From the figures plotting the average zone temperature (Fig. 11), the local zone temperatures (Fig. 12), the control temperatures (Fig. 13), the fraction of outside air (Fig. 14), the fraction of recirculated air (Fig. 15), the sensible load (Fig. 16) and the total power
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
1201
Fig. 14. Fraction of outside air in each RTU. Fig. 11. Average zone temperature profile compared to the outside temperature C.
Fig. 12. Local zone temperatures ( C).
Fig. 13. RTU control temperatures ( C).
used in each RTU (Fig. 17), we conclude the following. The control process can be divided into three regions: the first is from 0 h to
Fig. 15. Fraction of recirculated air in each RTU.
Fig. 16. Sensible load in each RTU (kW).
11 h where the comfort constraints are not active, the second is from 11 h to 22 h where the comfort constraints are active, and the
1202
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
Fig. 18. Zone temperatures under conventional control for August 28th. Fig. 17. Cooling coil power used in each RTU (kW).
third is from 22 h to 24 h where the comfort constraints are not active again. Concerning the first portion, all the units are off, and the supply temperatures are delivered in a way that guarantees a zero sensible load. There is no need for equipment pre-cooling here because the damper is fully open and the supplied mass flow rate comes completely from outside. Therefore, the temperature of the indoor air is already cool. Near the beginning of the occupancy period, the sensible loads at most of the zones (1, 2 and 4) become positive and the RTUs start to cool the area. This influences almost all time histories and shows an obvious discontinuity. In fact, the fraction of outside air jumps from 1 to 0.1 because the outside air is hot and it is more efficient to rely on the recirculated air. During this period, RTU 1 delivers greater levels of power than the other RTUs since it is the most efficient, while the sensible load at RTU 3 remains zero. Concerning the last period, the comfort constraints are relaxed, which drives the sensible load to return to zero at about 22 h.The total cost for this standard configuration is 689822 kJ. 5.2. Comparison with experimental conventional control Experimental measurements of the temperature at four data loggers in different zones of the restaurant were taken upon application of a classic conventional control. RTUs 2, 3 and 4 operate in only one compressor stage, whereas RTU 1 operates in two equivalent compressor stages. We present the cooling results for few days of August 2013. The data loggers sampling time is 5 min. The data represent the temperature at each zone and the corresponding power used to reach the given temperature. To find the power used, the data is integrated over a time interval of 24 h, and compared to the power used for the case of conventional control. We compare the two controllers for five consecutive days starting from August 25th, 2013 in Table 5, assuming that the
Fig. 19. Zone temperatures under optimal control for August 28th.
maximum temperature for the optimal control is 24.5 C. These comparison values are very sensitive to the maximum allowed indoor air temperature. The power reduction may reach 45% for specific cases. We plot the temperature response to both conventional (Fig. 18) and optimal controls (Figs. 19 and 20) on August 28th. The figures show that the optimal control gives more regular temperatures, whereas the conventional control shows a considerable difference between the zones and significant temperature fluctuations. Therefore, the optimal control is both more comfortable and energy efficient.
Table 5 Conventional/Optimal control cost. Date (August 2013)
25
26
27
28
29
30
Conventional cost(105 kJ) Optimal cost(105 kJ) Gain (%)
7.0083 5.7264 18.29
6.9046 5.4820 20.60
6.8156 5.3574 21.39
11.279 6.1908 45.11
7.2406 6.1410 15.18
9.4188 5.6805 39.69
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
1203
5.3. Parametric study of the optimal results
Fig. 20. Average temperature under optimal control for August 28th compared to the outside air temperature.
Table 6 Time-step sensitivity. Time-step dt (min)
7.5
3.75
1.875
Cost (kJ) Percentage
689822 0%
686909 0.4%
685600 0.6%
In this section, we study the impact of varying problem loads and variables on the final cost function. We start with the time convergence. All the previous simulations are done for dt ¼ 7.5 min. Table 6 shows the effects of the time step dt on the cost. It is shown that when the time step decreases, the cost also decreases but appears to converge. Fig. 21 plots the control supply air temperatures and shows the same trend as with the time step change. Next, we study the effects of some load variations on the final cost. Table 7 shows the effects of the dining load on the cost. For the reference configuration, the load rises from 0.2 to 1.2 during the occupancy period. We reduce the step size by 20%, then by 40%. The corresponding changes in cost are respectively a reduction of 6% and 12%. This result is predictable since reducing the load, requires less cooling power to reach the given constraint. Moreover, the ROM is linear. So when input (a) is reduced twice as much as input (b), the reduction in output (a) is twice as much as that of output (b). The same findings are revealed in Table 8 which shows a study of variations in the kitchen load. In fact, the kitchen load rises from 24 C to 28 C during occupancy time in the reference configuration. If we reduce the step size by 20%, then 40%, the change in the final cost is repectively 0.73% and 1.46%. Then, we check the effects of the occupancy period on the cost by simulating three starting hours: 11 (reference case), 12 and 13 as shown in Table 9. This variation is highly dependent on the external weather data and therefore is not expected to be linear.
Fig. 21. Supply control temperatures using different time steps.
1204
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
State University. The authors gratefully acknowledge constructive discussions exchanged with BLOM developers.
Table 7 Dining load sensitivity. Dining load
1.2
1.0
0.8
Cost (kJ) Percentage
689822 0%
648175 6%
606561 12%
Table 8 Kitchen load sensitivity. Kitchen load (C)
28
27.2
26.4
Cost (kJ) Percentage
689822 0%
684773 0.73%
679723 1.46%
Table 9 Occupancy period sensitivity. Control starting time (h)
11
12
13
Cost (kJ) Percentage
689822 0%
656182 4.9%
605763 12.2%
6. Conclusion A CFD model has been developed to simulate the indoor-air temperature and moisture of a restaurant. Different loads such as human occupancy and cooking emissions, and boundary conditions were modeled based on on-site and conventional data. After bringing the simulation to a steady state, single step input changes were applied to construct a linear time invariant representation of the coupled system. The reduced order model serves as a basis for the optimal control of the power required to operate the cooling rooftop units. The results show that our optimal control is more efficient and comfortable than a classic set-point conventional control. However, the real implementation of such controllers is not straightforward because it includes the actuation of the damper position and it assumes that the RTU stages in a continuous fashion, which requires the modification of the current cooling settings. This analysis can be used to determine if equipment upgrade would be worth the investment to reduce cooling power consumption. Acknowledgment This research was supported by the Air Force Office of Scientific Research under grant FA9550-10-1-0201 and by the DOE contract DE-EE0004261 under subcontract # 4345-VT-DOE-4261 from Penn
the
References [1] Naidu D, Rieger CG. Advanced control strategies for HVAC&R systemseAn overview: Part I: hard control. HVAC&R Res 2011:11216. [2] Naidu D, Rieger CG. Advanced control strategies for HVAC&R systemseAn overview: Part II: soft and fusion control. HVAC&R Res 2011;17(2):144. [3] Sane H, Haugstetter C, Bortoff S. Building HVAC control systemsrole of controls and optimization. In: Proceedings of the 2006 American control conference, Minneapolis, MN, vol. 17(1); 2006. p. 2e21. [4] House J, Smith T. Optimal control of building and HVAC systems. Proceedings of the American control conference. vol. 6 Seattle, WA, 432630. [5] Yu F, Chan K. Economic benefits of optimal control for water-cooled chiller systems serving hotels in a subtropical climate. Energy Build 2014;42(2):2039. [6] Oldewurtel F, Parisio A, Jones C, Gyalistras D, Gwerder M, Stauch V, et al. Use of model predictive control and weather forecasts for energy efficient building climate control. Energy Build 2012;45:15e27. [7] University of WisconsineMadison. Solar energy laboratory and Klein, Sanford A. TRNSYS, a transient system simulation program. Solar Energy Laborataory, University of WisconsineMadison; 1979. [8] Crawley DB, Lawrie LK, Winkelmann FC, Buhl WF, Huang WJ, Pedersen CO, et al. EnergyPlus: creating a new-generation building energy simulation program. Energy Build 2001;33(4):319e31. [9] Ma Y, Kelman A, Daly Y, Borrelli F. Predictive control for energy efficient buildings with thermal storage. IEEE Control Syst Mag 2012;32(1):44e64. [10] Ma Y, Borrelli F, Hencey B, Coffey B, Bengea S, Haves P. Model predictive control for the operation of building cooling systems. IEEE Trans Control Syst Technol 2012;20(3):796e803. [11] Siroky J, Oldewurtel F, Cigler J, Prvara S. Experimental analysis of model predictive control for an energy efficient building heating system 2011. Appl Energy 2011;88:30793087. [12] Ma J, Qin J, Salsbury T, Xu P. Demand reduction in building energy systems based on economic model predictive control. Chem Eng Sci 2012;67:92e100. [13] Zlatanovi I, Gligorevi K, Ivanovi S, Rudonja N. Energy-saving estimation model for hypermarket HVAC systems applications. Energy Build 2011;43:33533359. [14] Srebric J, Chen Q, Glicksman LG. A coupled airflow-and-energy simulation program for indoor thermal environmental studies. ASHRAE Trans 2000;106(1):465e76. [15] Zhai Z, Chen Q, Hayes P, Klems J. On approaches to couple energy simulations and computational fluid dynamics programs. Build Environ 2000;37:857e64. [16] Borggaard J, Cliff E, Gugercin S. Model reduction for indoor-air behavior in control design for energy-efficient buildings. In: Proceedings of 2012 ACC; 2012. [17] Fluent user's manual. 2006. [18] Gugercin S, Antoulas AC, Beattie CA. H 2 model reduction for large-scale linear dynamical systems. SIAM J Matrix Analysis 2008;30:609e38. [19] Kung SY. A new identification and model reduction algorithm via singular value decomposition. In: Proc 12th Asilomar Conf Circuits, Syst Comput; 1978. [20] Cliff EM, Borggaard JT, Braun JE, Gugercin S, Kim D. Coupled CFD/building envelope model for the purdue living lab. In: 2012 international high performance buildings conference at purdue. West Lafayette, IN: Purdue University; 2012. [21] Kim D, Braun JE. A general approach for generating reduced-order models for large multi-zone buildings. J Build Perform Simul 2015;8(6):435e48. [22] Zhai ZJ, Chen QY. Performance of coupled building energy and CFD simulations. Energy Build 2005;37:333e44.
Contents lists available at ScienceDirect
Energy journal homepage: www.elsevier.com/locate/energy
Optimal control of indoor-air cooling in buildings using a reduced order model S. Ben Ayed a, *, D. Kim b, J.T. Borggaard a, E.M. Cliff a a b
Interdisciplinary Center for Applied Mathematics, Virginia Polytechnic Institute and State University, Blacksburg, VA, USA School of Mechanical Engineering, Purdue University, West Lafayette, IN, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 30 November 2015 Received in revised form 15 August 2016 Accepted 7 October 2016 Available online 27 October 2016
In this work, we use Computational Fluid Dynamics (CFD) to generate the distributed dynamic responses of temperature and moisture in a restaurant that correspond to perturbations of input variables. The CFD model is validated with experimental data of the roof top units return temperatures. A Reduced Order Model (ROM) is developed by approximating the responses to these perturbations using a linear timeinvariant model. The resulting indoor-air model is coupled to a dynamic envelope model with longer time scales. Assuming the outside temperature distribution and the occupants' loads are known throughout the day, optimal control is applied to our coupled model to minimize the cooling power subject to local and global comfort constraints. The results show that the applied method is more efficient and comfortable than the results given by the experimentally measured response of a conventional set-point control strategy. Published by Elsevier Ltd.
Keywords: Optimal control Building energy Reduced order model Berkeley library for optimization modeling Envelope model Computational fluid dynamics
1. Introduction According to The U.S. Energy Information administration, in 2013, 40% of total energy consumption in the U.S. was consumed in residential and commercial buildings. This generally consists of heating-ventilation-air-conditioning (HVAC), lighting, and plug loads. Thus, even a modest reduction in building energy use represents significant contributions toward the relief of a host of pressing environmental issues. Here we focus on cooling requirements and specifically on the development of models suitable for study of predictive control of the indoor environment. Optimal control of the HVAC energy in buildings has been the object of various research works during the last decades. Naidu and Rieger wrote an extensive review about advanced control strategies in HVAC systems ([1,2]). Sane et al. [3] presented a survey of problems in controls and optimization related to commercial chilled water cooling plants. They included in multivariable optimal control, dynamic resource allocation and control over networks. Several works have explicitly reported the savings due to the use of optimal control for HVAC systems. House and Smith [4] used optimal control of the supply air temperature to minimize the cost
* Corresponding author. E-mail address: [email protected] (S. Ben Ayed). http://dx.doi.org/10.1016/j.energy.2016.10.022 0360-5442/Published by Elsevier Ltd.
for energy consumption. Comparing this method to the conventional control, they found a cost savings of 11%. Yu and Chan ([5]) used an optimal control strategy to control the cooling towers and condenser water pumps in a water-cooled chiller system. They concluded that their procedure contributed to annual energy savings of 8.6% and cost savings of 9.9%. For a good estimate of the optimal energy cost, thermal load modeling is very important. Various models and computational techniques have been adopted. Oldewurtel et al. [6] claims that computer-aided modeling tools like TRNSYS [7] and EnergyPlus [8] are mainly designed to estimate the energy use of building and cannot be readily used for control. They suggested the use of a thermal Resistance-Capacitance (RC) network for control purposes. Many others [9e11] used the same simple RC model to estimate the thermal load. However, Ma et al. [12] used a black box model structure instead of an RC model. They used EnergyPlus to model a commercial building, generated outputs for a nominal configuration, used the input/output data to identify mathematical models for power and temperature, then generated the optimal control. Zlatanovi [13] used experimental data to generate the energy saving estimation model. This was based on monitoring the energy consumption of each HVAC system appliance in a large trade center on a daily basis for one year. Modeling the indoor environment is very important in small to
1192
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
medium commercial buildings that utilize rooftop units (RTUs) to condition large open spaces. Examples include retail stores, restaurants and factories. In fact, for these buildings, the well-mixed assumption is not practical due to non-uniformly distributed diffuser systems, local cooling/heating actions corresponding to each RTU, and non-uniform heat gains. In the present work, we develop a coupled CFD-envelope model having two distinct time scales for a specific restaurant. The CFD model captures important dynamics of the spatial distributions of temperature and moisture inside the building; while the envelope model studies possible inside/outside heat transfer mechanisms, namely convection, conduction and radiation. Various approaches to couple prediction codes for fluid dynamics and solid conduction in a problem of conjugate heat transfer have been studied [14,15]. We use optimal control with this coupled model to minimize the cooling energy in a given restaurant. Since the optimization procedure is iterative, a reduced order model that captures the essential dynamics of the problem and has a low computational cost [16,18,19] is developed for performance analysis of optimal control for open-spaces by utilizing ROM. For this study, the reduced order model is created from CFD simulations such that the steady-state responses are adequately captured and such that the resulting model easily couples to the envelope model. 2. CFD of the restaurant indoor air and reduced order model Our approach to developing a reduced-order model is to construct a linear time-invariant (LTI) model that approximates the input-output behavior observed in computational fluid dynamic (CFD) simulations. Accordingly, a CFD model for the restaurant was developed using FLUENT. The coupling between the indoor-air and the envelope models is naturally achieved by treating the bounding surface temperatures as inputs to the CFD and the corresponding surface heat fluxes as outputs; that is, we use CFD to approximate the Dirichlet-to-Neumann (dynamic) response. Other CFD inputs include the mass-flow-rate, temperature, and humidity supplied by each of the roof-top-units (RTU), the sensible and latent loads provided by occupants, and external loads generated by the kitchen activity. Additional outputs of interest include air temperature/ humidity at selected zones/locations, and return air properties. The mesh is formed of 4.70 m tetrahedral cells with 624 k vertices. We use a refinement of about 1/3 around the diffuser and returns and approximately 1/7 in the wall regions where the thermostat temperatures are averaged (around 250 tetrahedral faces). The moisture is treated by integrating the species transport model in FLUENT. Hence the fluid is a mixture of air and water vapor. This FLUENT feature models the mixing and transport of species by solving conservation equations describing convection and diffusion for the species. The turbulence model used in the work is a standard kEpsilon model with standard wall functions. In fact, it is a simple two-equation model in which the solution of two separate transport equations allows the turbulent velocity and length scales to be independently determined. Robustness, economy, and reasonable accuracy for a wide range of turbulent flows explain this model's popularity in industrial flow and heat transfer simulations [17].
North Wall: 27.00 C East Wall: 25.58 C West Wall: 28.43 C Floor: 24.92 C Ceiling: 25.48 C
2.1.2. Kitchen interface The south wall boundary of our computational domain (see Fig. 1) separates the dining zone from the kitchen, which is not modeled. For the solid part of this boundary we impose no-slip flow and adiabatic thermal conditions. For the opening(s) between the kitchen and the dining zones (red) we have prescribed the flow from the kitchen. mass-flow-rate: 0.591 kg s temperature: 25.00 C water mass-fraction: 16 g/kg
2.1.3. Supply air Each supply air vent is a 2 ft by 2 ft ceiling mounted square. The air supply from each RTU is uniformly distributed among the associated vents, and is supplied at 16.60 C. The velocity at each vent is normal (i.e. downward), and is spatially distributed using a flattened parabolic profile. RTU#1- mass-flow-rate: 2.9320 kg/s; water mass-fraction kg; distributed over 19 supply vents (blue) RTU#2- mass-flow-rate: 0.5320 kg/s; water mass-fraction kg; distributed over 04 supply vents (red) RTU#3- mass-flow-rate: 0.3990 kg/s; water mass-fraction kg; distributed over 03 supply vents (purple) RTU#4- mass-flow-rate: 0.5320 kg/s; water mass-fraction kg; distributed over 04 supply vents (green)
7 g/ 9 g/ 9 g/ 9 g/
Note that RTU#3 supplies four vents but only three of these are in the computational domain. 2.1.4. Return air Each of the return vents is ceiling mounted. These are specified as pressure outlet with the flow normal to the boundary, and the pressure at 40 Pa below that on the floor. The pressure difference approximates the static head loss under vertical equilibrium (i.e. DP¼r g Dh) at a ceiling height of 3.4 m.
2.1. Nominal boundary/source conditions The indoor-air response is decomposed into a nominally steady part and a perturbation. For the steady part we apply fixed boundary and source conditions enumerated here: 2.1.1. Bounding surfaces For five of the bounding surfaces we prescribe no-slip flow and fixed temperature (Dirichlet) data
Fig. 1. Grid with zones - the left (green) wall is North. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
2.1.5. Dining area source The dining area is the roughly 69 m3 volume shown in yellow in Fig. 1. The intent is to simulate the sensible/latent loads produced by the occupants. In this volume we uniformly distribute an energy source of approximately 8,750 W, and a moisture source of approximately 3.08 g/s. 2.1.6. Wine bar area source The winebar area is the roughly 23 m3 volume that surrounds the winebar in the north end of the domain. The intent is to simulate the sensible/latent loads produced by the occupants. In this volume we uniformly distribute an energy source of approximately 3,310 W, and a moisture source of approximately 1.10 g/s. 2.2. Steady state response The system was subjected to fixed boundary and source conditions chosen to represent a mean condition in the space (see Section 2.1). FLUENT was run in both steady and unsteady modes; the latter consisted of 7 h of simulated time with external (boundary & source) conditions fixed. Some responses for the last hour are shown in Fig. 2. The black horizontal lines are averages over the 1-h window. The meaning of these various temperatures can be understood in terms of Fig. 1. In particular dining zone is the yellow rectangular volume in the center of the figure. wine_bar zone is the yellow volume in the shape of a reflected D in the left of the figure. fireplace is the surface of the red box in the left foreground (to the left of the entrance bump-in).
1193
that a steady flow exists, moreover the computed results depend on the grid which is likely somewhat too coarse.
2.3. FLUENT model validation In the early morning of 21 August 2013, a test to validate the FLUENT results was conducted at the Harvest Grill site. At approximately 03:30 the supply air temperature from RTU1 was reduced from about 16.5 C to about 9.5 C as shown in Fig. 3. A FLUENT run was made wherein a piecewise linear approximation of this temperature profile was prescribed for the RTU1 supply air while (most) other source/boundary conditions were maintained at nominal values. The only exception was that the temperature of the (supply) kitchen air was changed from the nominal 25 Ce21 C. The change was made because at the time of the experiment there was no cooking activity in the kitchen. The chosen value (21 C) was the ambient temperature in the RTN3 return air which is closely coupled to the kitchen air. Initial data for the simulation was a result of previous steady state runs. Fig. 4 displays the comparisons of the FLUENT and measured data for the return air temperatures in RTN1 and RTN2. Note that these start at rather different values; that is, the initial condition in the FLUENT simulation does not well-represent the overnight conditions of the test. To compare the trends we shifted the FLUENT predictions so that the initial values match. Fig. 5 displays the shifted comparisons for all four return air locations. In the time interval [03:30, 05:00] the FLUENT predictions are very close to the data. After 05:00 the supply air temperatures for the other RTU's began to change; a feature not captured in the FLUENT simulation.
2.4. Input variables
Note that the room average is the least oscillatory: it is the massaveraged temperature over the entire room. The dining and wine_bar values are mass-averaged over the cited volumes. The fireplace value is averaged over the external surfaces of the designated red-box, and exhibits the largest oscillations. Also, note that the values at the right end (t ¼ 0) do not generally match the average, and this is particularly noticeable in the fireplace data. One might expect that the oscillations would decay and the internal flow in the room would become steady. In fact there is no guarantee
Twenty-two of the boundary conditions discussed in Section 2.1 were identified as useful for injecting control or disturbance to the system. These are enumerated in Table 1. For each input we ran a FLUENT simulation beginning from the final state of the steady simulation and proceeding for 7200 s. In each run, one of the inputs in Table 1 was subjected to a ramp input beginning at t ¼ 8 s, ramping to the appropriate perturbation value over a 60 s span, then remaining constant at the perturbed value.
Fig. 2. Steady response; Zone Average Temperatures.
Fig. 3. RTU1 supply air temperature.
1194
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
Fig. 4. RTN1 (left), RTN2 (right).
Fig. 5. Shifted comparisons: RTN1 (upper left), RTN2 (upper right) RTN3 (lower left), RTN4 (lower right).
2.5. LTI model Our approach is to construct an LTI model that approximates the perturbed CFD responses generated by FLUENT simulations. Whereas there are a number of approaches to extract reduced-order LTI models from response data (for example [18,19]), we are motivated to find a low order model with the appropriate asymptotic behavior in a computationally efficient way. Nominally Steady Outputs In order to describe perturbed responses in a meaningful way, we first deal with defining steady values of the output variables. The output variables are listed in Table 2. Outputs 1 to 5 represent the average thermal flux through the walls which will be coupled to a building envelope model as shown in Fig. 10. Outputs 6 to 9, 10 to 13 and 14 to 17 represent the return air conditions, i.e. the average of the mass flow rates,
temperatures and humidities respectively. These are used to calculate HVAC power consumption. Outputs 18 to 21 indicates local temperatures at thermostat locations. In order to measure comfort, air conditions at some local area of interests and zonal volume averaged temperatures and humidities were chosen as outputs. Outputs 22 to 33 are the conditions at the local points while the rests are the zonal conditions. The final two outputs 42 and 43 are left for the average temperature and moisture of the whole space. It is clear from Fig. 2 that notable oscillations persist in the output variables. There are several potential explanations for this observed behavior: 1. There is no guarantee that the Reynolds-Averaged Navier-Stokes PDE model has a steady solution;
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
1195
samples). For each input, the perturbed output values were computed as:
Table 1 Input variables: Nominal and perturbation values. Number
Name
Nominal Value
Increment
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
north_temp east_temp west_temp floor_temp ceiling_temp rtu_1_temp rtu_2_temp rtu_3_temp rtu_4_temp rtu_1_h2o rtu_2_h2o rtu_3_h2o rtu_4_h2o rtu_1_mass rtu_2_mass rtu_3_mass rtu_4_mass load_dining load_winebar kitchen_inflow_temp kitchen_inflow_h20 kitchen_inflow_mass
27.00 C 25.58 C 28.43 C 24.92 C 25.46 C 16.60 C 16.60 C 16.60 C 16.60 C 7.00 g/kg 9.00 g/kg 9.00 g/kg 9.00 g/kg 2.9320 kg/s 0.5320 kg/s 0.3990 kg/s 0.5320 kg/s 100% 100% 25.00 C 16.00 g/kg 0.5910 kg/s
5.00 C 5.00 C 5.00 C 5.00 C 5.00 C 5.00 C 5.00 C 5.00 C 5.00 C 1.40 g/kg 1.80 g/kg 1.80 g/kg 1.80 g/kg 0.5864 kg/s 0.1064 kg/s 0.0798 kg/s 0.1064 kg/s 20% 20% 5.00 C 3.20 g/kg 0.1180 kg/s
ypert tJ ¼ Ydata tJ Ynom
J ¼ 1; …; 1800:
Best Fit There are several elegant approaches to extract a model from time-response data ([19,20]). Our approach decomposes the problem into sub-problems by treating each input response separately. That is, we hypothesize a structure for a low-order model and formulate an optimization problem to find appropriate parameter values. Specifically, for each input we propose a firstorder model of the form
1 _ ¼ xðtÞ þ b uðtÞ xðtÞ
(1)
yðtÞ ¼ C xðtÞ þ D uðtÞ;
(2)
t
where t is a time constant, b is a scalar, and C; D2IR43 . We take b ¼ 0.01 (scaling). Most of the 2243 input/output pairs are represented as strictly-proper transfer functions (D ¼ 0). However, for the first five inputs (surface temperatures) the heat flux output for
Table 2 Selected nominal steady outputs. Number
Name
Value
Number
Name
Value
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
q_north q_east q_west q_floor q_ceiling rtn_1_mass rtn_2_mass rtn_3_mass rtn_4_mass rtn_1_temp rtn_2_temp rtn_3_temp rtn_4_temp rtn_1_h2o rtn_2_h2o rtn_ 3_h2o rtn_4_h2o thermo_1_temp thermo_2_temp thermo_3_temp thermo_4_temp dk_2_temp
11.5095 w/m2 10.4278 w/m2 19.3538 w/m2 15.8151 w/m2 4.9463 w/m2 2.4959 kg/s 0.8248 kg/s 0.8248 kg/s 0.8392 kg/s 22.1082 C 22.3562 C 21.7175 C 22.6413 C 9.3741 g/kg 9.2538 g/kg 9.2147 g/kg 10.5528 g/kg 21.5048 C 21.80948 C 21.7650 C 21.8446 C 21.7961 C
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43
dk_4_temp dk_2_h2o dk_4_h2o rtu1_z1_temp rtu2_z1_temp rtu3_z1_temp rtu4_z1_temp rtu1_z1_h2o rtu2_z1_h2o rtu3_z1_h2o rtu4_z1_h2o dining_room_temp dining_room_h2o winebar_temp winebar_room_h2o winebar_skirt_temp winebar_skirt_h2o fireplace_temp fireplace_h2o room_avg_temp room_avg_h2o
21.7974 C 9.9237 g/kg 9.8813 g/kg 21.1414 C 22.2903 C 21.3036 C 22.0641 C 8.8190 g/kg 9.903 g/kg 9.4225 g/kg 9.9617 g/kg 21.9751 C 9.2293 g/kg 22.3145 C 10.0578 g/kg 21.65525 C 9.8574 g/kg 21.7940 C 9.6071 g/kg 22.0661 C 9.4163 g/kg
2. The grid may not be sufficiently resolved to represent the PDE behavior; and/or, 3. The iterative procedure for the nonlinear algebraic system may not be sufficiently converged. Whatever the cause(s) we simply define steady values for the outputs by averaging over a time-window. In particular, steady output values are computed by averaging the observed outputs over the final 60 min of the nominally steady solution. Steady values for the four output variables in Fig. 2 are indicated by solid lines. Steady values for all of the output variables are listed in Table 2. CFD Perturbations As noted above, in each run, one of the inputs in Table 1 was subjected to a ramp input beginning at t ¼ 8 s, ramping linearly to the corresponding perturbation value over 60 s, then remaining constant at the perturbed value. In each case the output variables (see Table 2 were recorded at 4 s intervals (1800
that surface has a single non-zero entry in D, since the surface temperature is directly coupled to the heat flux on that surface. Thus, we expect that for each surface heat flux output (the first five) we have DJ ¼ hdı J , whereas for all other outputs we have D ¼ 0. The time-scales of the building envelope are typically much longer than those of the indoor-air. Therefore, during most of the response the indoor-air system will be near equilibrium. For this reason we require that the ROM accurately represent the indoor-air dynamics for responses near equilibrium. This requirement is not naturally imposed in other model reduction methods and leads to our ad-hoc approach. Since the input perturbation function u goes to a constant value (say, uss), then with t > 0 the dynamics (1) will go to the stable equilibrium
xss ¼ tbuss 0yss ¼ ðC t b þ DÞuss : Solving for C we find
1196
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
Cðt; DÞ ¼
yss =uss D ; tb
(3)
where the vector D2IR43 has (at most) a single non-zero entry. To implement (3) we must specify a method for computing a steadystate value of the perturbed output, i.e. yss 2IR43 . Because the output data includes residual oscillations it is prudent to use some form of averaging. For the results presented here we averaged y(tı) over the last 120 samples (out of 1800 total), that is ı2f1681; 1682; …1800g. This amounts to an 8 min wide sampling window. With these preliminaries we define a scalar-valued error function
Jðt; DÞ ¼
1800 Xh
i2 ypert ðtı Þ yðtı ; t; DÞ ;
(4)
ı¼1
where yðtı ; t; DÞ is computed from (1)e(2) with the state-to-output map given by (3), and the input (u) is the step-like input described in Section 2.5. For each input we compute optimal values of t and (possibly) D by minimizing the error functional J in (4). These calculations were done using MATLAB's fminsearch procedure. The five non-zero values for the D arrays are; 1. 2. 3. 4. 5.
q_north: 2.912 W/m2/ C q_east: 2.038 W/m2/ C q_west: 3.514 W/m2/ C q_floor: 4.281 W/m2/ C q_ceiling: 1.423 W/m2/ C
Note that these can be interpreted as (time and spatial) averages of the film cooling coefficient on the respective surfaces.
2.6. ROM-FLUENT comparisons Fig. 6 provides comparisons of the FLUENT and the reduced-ordermodel (ROM) responses to a change in the east wall temperature. The input perturbation ramps to þ5 C over 60 s and the corresponding ramp in the East wall heat flux is seen in the upper graph. In the lower graph the ceiling heat flux is seen to decrease by 0.6 W/ m2, presumably the result of an increase in the air temperature engendered by the increased wall temperature. Note that the graphs display only the first 15 min of the response; the data fits are based on 2 h of data. Similarly, Fig. 7 provides comparisons of the FLUENT and the reduced-order-model (ROM) responses to a 5 C decrease in the RTU#1 supply air temperature. The upper graph displays the change in return 1 air temperature, whereas the lower graph characterizes the temperature in the dining zone. Fig. 8 provides also comparisons of the FLUENT and the reducedorder-model (ROM) responses to a 3.2 g/kg increase in the kitchen inflow moisture. The upper graph displays the change in return 1 water vapor mass fraction, whereas the lower graph shows the response of the dining zone water vapor mass fraction.
3. Development of a thermal building envelope model In order to connect a reduced-order indoor air model to a building envelope model, a general formula derived in Ref. [21] is adopted. The majority of the derivations in this section is the same as [21], but it is slightly modified for the coupling.
Fig. 6. Comparison of Fluent and ROM T_east responses.
3.1. Conduction through walls A finite volume formulation is used to describe the heat conduction through walls and is depicted in Fig. 9. For any jth node in a wall except the first and last nodes, an energy balance leads to
rij Cji wij
dTji dt
i i i i i i i ¼ hLcd Tj1 hLcd þ hRcd Tji þ hRcd Tjþ1 þ qgen ; j
j
j
j
j
(5) where rij is the density at jth node in ith wall [kg/m3], Cji is the specific heat capacity of jth node in ith wall [J/kgK] and wij is the i i kL j width of control volume of jth node in ith wall [m]. hC DL ¼ ji j wL j j i where kL is the thermal conductivity at left surface of the jth node j i in ith wall [W/m2K] and wL is the distance from the (j1)th node j i to the jth node in ith wall [m].qgen is an energy source [W/m2] j
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
Fig. 8. Comparison of Fluent and ROM h_kitchen responses. Fig. 7. Comparison of Fluent and ROM T_rtu_1 responses.
inside the jth finite control volume that belongs to the ith wall. 2 q 1 3 gen j 2 7 6 7 ! 6 Using the following matrix notation, q j ≡6 qgen 7, j 4 :: 5 Nw qgen j
2
Tj1
3
6 2 7 ! 6 Tj 7 7; T j ≡6 6 7 4 :: 5 Nw Tj
2 L=R 1 hcd j 6 6 ~ L=R ≡6 H 6 j 6 4
3
0 L=R 2 hcd j
::
0 1 2 rCw j 6 6 ~ ≡6 C j 6 4 0
L=R Nw hcd
7 7 7 7; and 7 5
j
3
0 2 rCw j
::
Nw
rCw
7 7 7; 7 5
j
Eq. (5) can be expressed in the following matrix form
Fig. 9. Notation for conduction through walls.
1197
1198
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
3 3 i i qi1 ¼ hicv;ex Ta þ 4sεi1 Fsky T sky Tsky þ Fgrd T grd Tgrd þ ai1 qiSWR ; T sky ¼
Tsky þ T1 ; 2
T grd ¼
Tgrd þ T1 : 2
The mean temperatures for long-wave exchange between the surface and sky and surface and ground are assumed to be the same for all outside surfaces. Eq. (8) can be generalized and written in compact matrix form as
! ! R ! ~ dT 1 ¼ H ~ ~ ~R ! ~ C 1 cd;1 þ H rad;ex þ H cv;ex T 1 þ H cd;1 T 2 þ q 1 ; dt
(9)
where
3 3 ! i i ð q 1 Þi ¼ hicv;ex Ta þ 4sεi1 Fsky T sky Tsky þ Fgrd T grd Tgrd
Fig. 10. Simplified scheme of ROM-envelope coupling.
þ ai1 qiSWR : ! L ! ~ dT j ¼ H ~L ! ~R ! ~R ! ~ C j cd;j T j1 H cd;j þ H cd;j T j þ H cd;j T jþ1 þ q j : dt
(6)
! Note that T j is a group of all of the wall temperature nodes belonging to an iTindividual zone. More precisely, h ! T j ¼ Tj1 Tj2 :: TjNw whereNw is the number of walls in a zone and the superscript T represents the transpose.
3.3. Heat balance at the interior surfaces For the nth zone and ith wall, the energy balance equation for the inside surface is
rin Cni win 3.2. Heat balance at the exterior surfaces For any ith outside wall (connected to the external ambient) belonging to an individual zone, the heat balance equation at the surface is
dT i ri1 C1i wi1 1 dt
¼
hicv;ex
Ta
T1i
þ
i hRcd T2i 1
T1i
þ
(7)
where T1i represents the wall temperature of the first node which is set to be an outside surface of the wall. With the assumptions that the outside surface is gray and diffuse and the air is a nonparticipating radiation media, net long wavelength interactions with the environment can be expressed as i qiLWR ¼ sεi1 Fsky
Tsky
4
4 4 4 i þ sεi1 Fgrd Tgrd T1i : T1i
Using a linear approximation of the long-wave heat exchange term gives
ri1 C1i wi1
i i dT1i ¼ hicv;ex þ hRcd þ hirad;ex T1i þ hRcd T2i þ qi1 ; 1 1 dt
where
hirad;ex
3 3 i i ; ¼ 4sεi1 Fsky T sky þ Fgrd T grd
(8)
(10)
where qinet;rad is net radiative flux out of the inside wall and qicv;in the net heat flux from the zone air and allows for the CFD coupling. The radiosity method is utilized to express the net flux under the assumption that the walls are gray, diffuse and opaque. The same linearization method used in Eq. (8) is employed leading to
h 0! i ! ~ 1 B ~ Tn! q net;rad ¼ A ho
ai1 qiSWR
þ qiLWR ;
i dTni i ¼ hLcd Tn1 Tni qinet;rad þ qicv;in ; dt n
(11)
where
~ ≡dij rj F ; A ij εj εj ij and 0 ~ ¼ 4s d F T 3 : B ij ij ij
Radiosity does not appear explicitly in the above expression, which is convenient for building simulations. Since hi0 represents an external radiative source acting on the ith surface, the effects of internal sources and transmitted solar energy though windows are treated in a consistent manner. For any shaped room, the net radiative flux can be explicitly calculated as a function of surface temperatures if the view factors and the external radiative sources ! ~ 1 ! ~ 1 ~ 0 ~ are known. By letting H h o, rad;in ≡A B and q n ≡A
! ! L ! ! ~nd T n ¼ H ~L ! ~ ~ C cd;n T n1 H cd;n þ H rad;in T n þ q n þ q n;cv;in : dt (12)
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
3.4. State-space representation of the thermal envelope network Gathering the system of equations that represent heat balance equations for external to internal wall elements, i.e. Eqs. (6), (9) and (12), _ ! ! ~w! ~ ww ! Tw ¼ H T w þ q w þ q cv;in ; C
(13)
where
2 ~w C
6 ¼6 4 2
~ C 1
3 ~ C 2
~ H 1 6 6 ~L H ¼6 6 cd;2 4
::
~n C
7 7; 5
~ ww H
~R H L cd;1 R ~ ~ H cd;2 þ H cd;2 ::
Table 3 Model variables. Variable Tc i ði ¼ 1 : 4Þ m_ s i ði ¼ 1 : 4Þ fs i ði ¼ 1 : 4Þ DL WBL KT m_ r i ði ¼ 1 : 4Þ Tr i ði ¼ 1 : 4Þ Tz i ði ¼ 1 : 4Þ fr ði ¼ 1 : 4Þ
Significance
ROM Input/Output
RTU control temperature
IN(6)to IN(9)
Ta di
i
3 ~R H cd;2 :: ~L H cd;n
li
7 7 7; 7 :: 5 ~ H n
~R þ H ~ cv;ex þ H ~ ~ n ¼ ðH ~L þ H ~ ~ ¼ ðH H H 1 rad;ex Þ, rad;in Þ, cd;1 cd;n h h i iT ! T T T T ! T T T ! ! ! ! ! Tw¼ T and q w ¼ ! q 1 q 2 :: q n . As T 2 :: T n 1 expected, a tri-diagonal block matrix is formed with parameters that characterize heat transfer due to conduction in the walls (from node number 2 to node number n-1) and radiative/convective heat ! transfer at the boundaries (first and last nodes only). The terms q j where j ranges from 2 to n1 vanish if there are no heat flux sources inside the wall such as embedded radiant heating or ! cooling systems. The term, q cv;in , represents the set of net convective heat fluxes to the zone air and allows CFD coupling. All entries in the vector are zeros except for inner walls. 3.5. Coupling with the restaurant CFD model The one dimensional heat balance equations, Eq. (13), are coupled to a reduced-order indoor air model. Although it is possible to subdivide wall surfaces sufficiently to model temperature variations across the walls and solve three dimension heat equations, we chose the 1-D form for computational simplicity. However Zhai and Chen [22] showed that the isothermal wall dynamics coupled to a CFD model leads to reasonable accuracy in predicting spatial temperature distribution of zone air and cooling/heating load. The state-space modeling (Eq. (13)) is applied to the restaurant walls. It results in the following system of equations:
X_ e ¼ Ae Xe þ Be u Ye ¼ Ce Xe þ De u;
(14)
1199
RTU supply mass flow rate
IN(14)to IN(17)
Supply water mass fraction
IN(10)to IN(13)
Dining load Wine bar load Kitchen temperature RTU return mass flow rate
IN(18) IN(19) IN(20) OUT(6)to OUT(9)
Return temperature
OUT(10) to OUT(13)
Local zone temperature
OUT(26)to OUT(29)
Return water mass fraction
OUT(14)to OUT(17)
Average zone temperature Fraction of outside air Fraction of recirculated air
OUT(42) NA NA
values as specified in Section 2. The supply air temperatures, occupancy and cooking loads, walls, ceiling and floor temperatures coming from the envelope model are considered as perturbations to the nominal configurations. Thus, due to the linearity of the model, their effects on the ROM are additive. The control variables are the four RTU supply air temperatures and the four RTU outside air fractions. All the variables and problem parameters are summarized in Tables 3 and 4. The details of loads and occupancy are summarized as follows: Kitchen temperature (KT): 24 C rising to 28 C between t ¼ 10 h and t¼22 h. Dining load (DL): 0.2 rising to 1.2 between t ¼ 12 h and t¼22 h. Wine bar load (WBL): 0.2 rising to 1.2 between t ¼ 12 h and t¼22 h. Occupancy: the occupancy function is 1 between hmin and hmax, and 0 otherwise
4.2. Cost function (equipment model) rates, namely, In this analysis, we distinguish several mass flow m_ s i as the supply mass flow rate of each RTU, m_ oa i as the mass flow rate supplied from outside for each RTU and m_ r i as the return mass flow rate to each RTU, where the subscript i designates quantities related to the ith RTU. Continuity of these rates gives the following balance:
m_ s i ¼ m_ oa i þ li m_ r i :
(15)
Where li is the fraction of the return air supplied in each RTU. We denote by di the fraction of the supply air from outside. This fraction satisfies the following equation:
! where Xe is the state vector ( T w in Eq. (13)), u is the input vector of 22 elements comprising the outdoor weather variables and the fluxes through the walls computed from the ROM. The finite difference scheme with 2e4 nodes per material layer of a wall results in 35 elements in Xe. The number was selected after checking the grid independence of a simulation result. The output Ye represents a vector of the computed wall temperatures that are inputs to the ROM. This coupling is illustrated in a simplified scheme in Fig. 10.
Combining Eqs. (15) and (16), we get the following relation between di and li, where di is a control variable and li is computed as:
4. Optimal control of the cooling energy
li ¼ ð1 di Þ
4.1. Problem formulation
Since li is supposed to be less than or equal to unity, we distinguish two cases:
This work studies the restaurant indoor-air model for use in performing optimal control calculations. Unless modified, the default values of all the input variables are fixed at the nominal
m_ oa i ¼ di m_ s i :
m_ s : m_ r i
_ r if di 1 m , then li ¼ 1 m_ s i m_ r if di > 1 m_ s i , then li is computed using Eq. (17)
(16)
(17)
1200
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204 Table 4 Model parameters. Parameter Ta
Significance
Value
Upper bound for comfort
23 C
hmin,max Cp dmin,max Tmin,max copi(i ¼ 1:4) capmax i ði ¼ 1 : 4Þ foa dt L
Occupancy bounds Heat capacity Air fraction bounds Bounds for control temp Total coefficient of performance Total cooling capacity of each RTU
11 h, 22 h 1 kJ/(kg K) 0.1, 1.0 10 C, 30 C 3.5, 2.6, 2.6, 2.6 53, 14, 10.5, 14 kW
Outside water mass fraction Time step Cooling power for removal of 1 g of water
17 g/kg 7.5 min 2.26 kW
max
The total cooling load at each RTU is divided into sensible Qs i and latent Ql i loads. The sensible load related to the temperature change is expressed as:
Qs i ¼ Cp m_ s i Tm i Tc i
(18)
0:05 Tc i Troof ; Tc i ¼ Tc i 1:004
(19)
where Cp is the heat capacity of the air at constant temperature, Tm i is the RTU inlet air temperature, Tc i is the supply RTU tem perature, and Tc i is the modified supply RTU temperature obtained after the circulation of the supply air in the attic. The RTU inlet air temperature or the mixing temperature can be written as:
m_ r i li Tr i þ m_ oa i Toa Tm i ¼ ; m_ s i
(20)
where Toa is the outside air temperature and Tr i is the temperature at each return. The moisture formulation included in the CFD model allows us to easily compute the latent load related to the change of state, which can be expressed as:
Ql i ¼ Lðm_ r i li fr i þ m_ oa i foa fs m_ s i Þ;
(21)
where fr i , foa and fs i are the mass fraction of water respectively for the return, outside air and supply flows, and L is the amount of cooling required to remove 1 g of water. Thus, the total power expression is:
Qs i þ Ql i Pt i ¼ ; copi
(22)
2. Set the upper comfort bound for the local temperature of the different zones at the occupancy time: for hmin t hmax, Tz i Ta max . 3. Set the upper comfort bound for the average zone temperature at the occupancy time. For hmin t hmax, Ta Tajmax . 4. Cooling coil can only decrease the temperature. Tc i Tm i . 5. Set the upper and lower bounds for the outside air fraction. 0.1 di 1. 6. Set the upper and lower bounds for the recirculated air fraction. 0li 1. 7. Set the upper and lower bounds for the cooling power of each RTU.
0 Qs i þ Ql i capmax i 8. Latent heat is only removed during the occupancy period.
4.4. Numerical approximation In optimal control, one is concerned with minimizing the functional (23) from a given initial state over a forward planning horizon ([0,T]). Here we briefly study a discretized version of the problem using the MATLAB/SIMULINK based tool: Berkeley Library for Optimization Modeling. At this time BLOM 2.0 is restricted to discrete-time dynamics, so it is necessary to convert from a continuous to a discrete time model. Thus, for example, the envelope model (15) is replaced by
Xe kþ1 ¼ Fe Xe k þ Ge ;
(24)
where where copi is the coefficient of performance of RTU #i. We note that the RTUs are constant fan speed regardless of RTU stages. Therefore, fan powers are assumed to be constant and do not appear in the objective function. The final cost expression can be written as:
C¼
ZT X 4 0
Pt i dt
Fe ¼ expðAe dtÞ Ge ¼
(23)
i¼1
(25)
Zdt expðAe ðdt sÞÞBe ds
(26)
0
and dt is the time step.
where T is a 24-hour period. 5. Results 4.3. Constraints Most of the constraints are applied for each RTU i (i ¼ 1, 2, 3, 4): 1. Specify the upper and lower bounds for the control temperature for different RTUs. 10 C Tcji 30 C.
5.1. Optimal control results From the figures plotting the average zone temperature (Fig. 11), the local zone temperatures (Fig. 12), the control temperatures (Fig. 13), the fraction of outside air (Fig. 14), the fraction of recirculated air (Fig. 15), the sensible load (Fig. 16) and the total power
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
1201
Fig. 14. Fraction of outside air in each RTU. Fig. 11. Average zone temperature profile compared to the outside temperature C.
Fig. 12. Local zone temperatures ( C).
Fig. 13. RTU control temperatures ( C).
used in each RTU (Fig. 17), we conclude the following. The control process can be divided into three regions: the first is from 0 h to
Fig. 15. Fraction of recirculated air in each RTU.
Fig. 16. Sensible load in each RTU (kW).
11 h where the comfort constraints are not active, the second is from 11 h to 22 h where the comfort constraints are active, and the
1202
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
Fig. 18. Zone temperatures under conventional control for August 28th. Fig. 17. Cooling coil power used in each RTU (kW).
third is from 22 h to 24 h where the comfort constraints are not active again. Concerning the first portion, all the units are off, and the supply temperatures are delivered in a way that guarantees a zero sensible load. There is no need for equipment pre-cooling here because the damper is fully open and the supplied mass flow rate comes completely from outside. Therefore, the temperature of the indoor air is already cool. Near the beginning of the occupancy period, the sensible loads at most of the zones (1, 2 and 4) become positive and the RTUs start to cool the area. This influences almost all time histories and shows an obvious discontinuity. In fact, the fraction of outside air jumps from 1 to 0.1 because the outside air is hot and it is more efficient to rely on the recirculated air. During this period, RTU 1 delivers greater levels of power than the other RTUs since it is the most efficient, while the sensible load at RTU 3 remains zero. Concerning the last period, the comfort constraints are relaxed, which drives the sensible load to return to zero at about 22 h.The total cost for this standard configuration is 689822 kJ. 5.2. Comparison with experimental conventional control Experimental measurements of the temperature at four data loggers in different zones of the restaurant were taken upon application of a classic conventional control. RTUs 2, 3 and 4 operate in only one compressor stage, whereas RTU 1 operates in two equivalent compressor stages. We present the cooling results for few days of August 2013. The data loggers sampling time is 5 min. The data represent the temperature at each zone and the corresponding power used to reach the given temperature. To find the power used, the data is integrated over a time interval of 24 h, and compared to the power used for the case of conventional control. We compare the two controllers for five consecutive days starting from August 25th, 2013 in Table 5, assuming that the
Fig. 19. Zone temperatures under optimal control for August 28th.
maximum temperature for the optimal control is 24.5 C. These comparison values are very sensitive to the maximum allowed indoor air temperature. The power reduction may reach 45% for specific cases. We plot the temperature response to both conventional (Fig. 18) and optimal controls (Figs. 19 and 20) on August 28th. The figures show that the optimal control gives more regular temperatures, whereas the conventional control shows a considerable difference between the zones and significant temperature fluctuations. Therefore, the optimal control is both more comfortable and energy efficient.
Table 5 Conventional/Optimal control cost. Date (August 2013)
25
26
27
28
29
30
Conventional cost(105 kJ) Optimal cost(105 kJ) Gain (%)
7.0083 5.7264 18.29
6.9046 5.4820 20.60
6.8156 5.3574 21.39
11.279 6.1908 45.11
7.2406 6.1410 15.18
9.4188 5.6805 39.69
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
1203
5.3. Parametric study of the optimal results
Fig. 20. Average temperature under optimal control for August 28th compared to the outside air temperature.
Table 6 Time-step sensitivity. Time-step dt (min)
7.5
3.75
1.875
Cost (kJ) Percentage
689822 0%
686909 0.4%
685600 0.6%
In this section, we study the impact of varying problem loads and variables on the final cost function. We start with the time convergence. All the previous simulations are done for dt ¼ 7.5 min. Table 6 shows the effects of the time step dt on the cost. It is shown that when the time step decreases, the cost also decreases but appears to converge. Fig. 21 plots the control supply air temperatures and shows the same trend as with the time step change. Next, we study the effects of some load variations on the final cost. Table 7 shows the effects of the dining load on the cost. For the reference configuration, the load rises from 0.2 to 1.2 during the occupancy period. We reduce the step size by 20%, then by 40%. The corresponding changes in cost are respectively a reduction of 6% and 12%. This result is predictable since reducing the load, requires less cooling power to reach the given constraint. Moreover, the ROM is linear. So when input (a) is reduced twice as much as input (b), the reduction in output (a) is twice as much as that of output (b). The same findings are revealed in Table 8 which shows a study of variations in the kitchen load. In fact, the kitchen load rises from 24 C to 28 C during occupancy time in the reference configuration. If we reduce the step size by 20%, then 40%, the change in the final cost is repectively 0.73% and 1.46%. Then, we check the effects of the occupancy period on the cost by simulating three starting hours: 11 (reference case), 12 and 13 as shown in Table 9. This variation is highly dependent on the external weather data and therefore is not expected to be linear.
Fig. 21. Supply control temperatures using different time steps.
1204
S. Ben Ayed et al. / Energy 116 (2016) 1191e1204
State University. The authors gratefully acknowledge constructive discussions exchanged with BLOM developers.
Table 7 Dining load sensitivity. Dining load
1.2
1.0
0.8
Cost (kJ) Percentage
689822 0%
648175 6%
606561 12%
Table 8 Kitchen load sensitivity. Kitchen load (C)
28
27.2
26.4
Cost (kJ) Percentage
689822 0%
684773 0.73%
679723 1.46%
Table 9 Occupancy period sensitivity. Control starting time (h)
11
12
13
Cost (kJ) Percentage
689822 0%
656182 4.9%
605763 12.2%
6. Conclusion A CFD model has been developed to simulate the indoor-air temperature and moisture of a restaurant. Different loads such as human occupancy and cooking emissions, and boundary conditions were modeled based on on-site and conventional data. After bringing the simulation to a steady state, single step input changes were applied to construct a linear time invariant representation of the coupled system. The reduced order model serves as a basis for the optimal control of the power required to operate the cooling rooftop units. The results show that our optimal control is more efficient and comfortable than a classic set-point conventional control. However, the real implementation of such controllers is not straightforward because it includes the actuation of the damper position and it assumes that the RTU stages in a continuous fashion, which requires the modification of the current cooling settings. This analysis can be used to determine if equipment upgrade would be worth the investment to reduce cooling power consumption. Acknowledgment This research was supported by the Air Force Office of Scientific Research under grant FA9550-10-1-0201 and by the DOE contract DE-EE0004261 under subcontract # 4345-VT-DOE-4261 from Penn
the
References [1] Naidu D, Rieger CG. Advanced control strategies for HVAC&R systemseAn overview: Part I: hard control. HVAC&R Res 2011:11216. [2] Naidu D, Rieger CG. Advanced control strategies for HVAC&R systemseAn overview: Part II: soft and fusion control. HVAC&R Res 2011;17(2):144. [3] Sane H, Haugstetter C, Bortoff S. Building HVAC control systemsrole of controls and optimization. In: Proceedings of the 2006 American control conference, Minneapolis, MN, vol. 17(1); 2006. p. 2e21. [4] House J, Smith T. Optimal control of building and HVAC systems. Proceedings of the American control conference. vol. 6 Seattle, WA, 432630. [5] Yu F, Chan K. Economic benefits of optimal control for water-cooled chiller systems serving hotels in a subtropical climate. Energy Build 2014;42(2):2039. [6] Oldewurtel F, Parisio A, Jones C, Gyalistras D, Gwerder M, Stauch V, et al. Use of model predictive control and weather forecasts for energy efficient building climate control. Energy Build 2012;45:15e27. [7] University of WisconsineMadison. Solar energy laboratory and Klein, Sanford A. TRNSYS, a transient system simulation program. Solar Energy Laborataory, University of WisconsineMadison; 1979. [8] Crawley DB, Lawrie LK, Winkelmann FC, Buhl WF, Huang WJ, Pedersen CO, et al. EnergyPlus: creating a new-generation building energy simulation program. Energy Build 2001;33(4):319e31. [9] Ma Y, Kelman A, Daly Y, Borrelli F. Predictive control for energy efficient buildings with thermal storage. IEEE Control Syst Mag 2012;32(1):44e64. [10] Ma Y, Borrelli F, Hencey B, Coffey B, Bengea S, Haves P. Model predictive control for the operation of building cooling systems. IEEE Trans Control Syst Technol 2012;20(3):796e803. [11] Siroky J, Oldewurtel F, Cigler J, Prvara S. Experimental analysis of model predictive control for an energy efficient building heating system 2011. Appl Energy 2011;88:30793087. [12] Ma J, Qin J, Salsbury T, Xu P. Demand reduction in building energy systems based on economic model predictive control. Chem Eng Sci 2012;67:92e100. [13] Zlatanovi I, Gligorevi K, Ivanovi S, Rudonja N. Energy-saving estimation model for hypermarket HVAC systems applications. Energy Build 2011;43:33533359. [14] Srebric J, Chen Q, Glicksman LG. A coupled airflow-and-energy simulation program for indoor thermal environmental studies. ASHRAE Trans 2000;106(1):465e76. [15] Zhai Z, Chen Q, Hayes P, Klems J. On approaches to couple energy simulations and computational fluid dynamics programs. Build Environ 2000;37:857e64. [16] Borggaard J, Cliff E, Gugercin S. Model reduction for indoor-air behavior in control design for energy-efficient buildings. In: Proceedings of 2012 ACC; 2012. [17] Fluent user's manual. 2006. [18] Gugercin S, Antoulas AC, Beattie CA. H 2 model reduction for large-scale linear dynamical systems. SIAM J Matrix Analysis 2008;30:609e38. [19] Kung SY. A new identification and model reduction algorithm via singular value decomposition. In: Proc 12th Asilomar Conf Circuits, Syst Comput; 1978. [20] Cliff EM, Borggaard JT, Braun JE, Gugercin S, Kim D. Coupled CFD/building envelope model for the purdue living lab. In: 2012 international high performance buildings conference at purdue. West Lafayette, IN: Purdue University; 2012. [21] Kim D, Braun JE. A general approach for generating reduced-order models for large multi-zone buildings. J Build Perform Simul 2015;8(6):435e48. [22] Zhai ZJ, Chen QY. Performance of coupled building energy and CFD simulations. Energy Build 2005;37:333e44.