Optimal Price Based Demand Response of HVAC Systems in Commercial Buildings Considering Peak Load Reduction

: Electric utility companies (EUCs) play an intermediary role of retailers between wholesale market and end-users, maximizing their proﬁts. Retail pricing can be well deployed with the support of EUCs to promote demand response (DR) programs for heating, ventilating, and air-conditioning (HVAC) systems in commercial buildings. This paper proposes a pricing strategy to help EUCs and building operators achieve an optimal DR of price-elastic HVAC systems, considering peak load reduction. The proposed strategy is implemented by adopting a bi-level decision model. The nonlinear thermal response of an experimental building room is modeled using piecewise linear equations, which helps convert the bi-level model to the single-level model. The pricing strategy is implemented considering a time-of-use (TOU) pricing scheme, leading to low price volatility. Case studies are conducted for two types of load curves and the results demonstrate that the proposed strategy helps EUC promote the price-based DR of the commercial buildings for conventional load curves. However, EUC cannot reduce the peak load on duck curve caused by the large introduction of photovoltaic generators, even with price-sensitive HVAC systems in commercial building. This will be addressed in future studies by inducing DR participation of HVAC systems in residential buildings.


Introduction
The power consumption in the distribution network is predicted to significantly increase owing to the electrification of transport and heating [1]. In addition, the peak load is expected to double by 2050 [2]. The peak load in the conventional load curve usually occurs during the daylight hours. Therefore, the issue of peak load may seem to be insignificant in the new load curve changed by large introduction of photovoltaics (PV). However, the large introduction of PV has caused the peak load to occur after sunset, but the importance of peak load reduction remains unchanged [3]. Increasing the peak load requires to upgrade conventional substations and lines or adding new substations and lines. As this increases the operating costs, the distribution system operators (DSOs) or electric utility companies (EUCs) try to reduce the peak load in the distribution network [4].
In this environment, a EUC should maintain the operating stability of the distribution network while maximizing its profit. With the development of information technology, the demand side has become increasingly important to solve the problem such as peak load reduction in a power system. To accomplish one's role, the EUC uses the demand response (DR), which is the behavior of consumers adjusting their normal power consumption patterns through voluntary participation owing sunset, when PV is no longer available. However, few studies have attempted to solve the duck curve problem using DR resources.
In this paper, we propose a bi-level optimal pricing model for the EUC to reduce the peak load of two types of load curves using the experimental data-driven model of an HVAC system in commercial buildings. Specifically, in the upper level, the EUC maximize its profit based on advantages of locational marginal pricing (LMP) and time-of-use (TOU) rates. While in the lower level, building end-users schedule the optimal power inputs of HVAC systems to minimize the electricity bills according to the optimal retail rates. This allows the EUC to obtain more accurate and useful insights into the load shifting or curtailment capacities of HVAC systems, when the proposed pricing strategy is applied to DR programs in practice. Moreover, this enables building managers to ensure the thermal comfort of occupants and, consequently, participate more fully in DR programs due to improved comprehension of the inherent thermal energy storage capacity in their building structures.
In addition, the proposed pricing strategy is applied to the peak load reduction of two types of load curve patterns: a conventional load curve and duck curve. The proposed pricing method induces DR of HVAC system in commercial buildings, and through this, we analyzed the effect of the proposed pricing strategy on the peak load reduction. Further, the proposed pricing model was demonstrated through explicit analysis of various case studies such as change of DR purpose and constraints' parameter.
Section 2 explains the pricing strategy framework using the thermal response of the HVAC system in a commercial building. Section 3 presents the optimization problem formulation for the proposed method. Section 4 discusses the simulation case study results for the proposed pricing strategy. Section 5 provides our conclusions. Figure 1 shows the framework of the proposed optimal retail pricing strategy. The EUC buys electricity from the day-ahead wholesale market at LMP and sells it to consumers at two retail prices, namely TOU and real-time pricing (RTP). The LMP represents the price to buy and sell power at different locations within wholesale electricity markets. In this paper, LMP can be considered simply as a wholesale electricity price. The TOU rates typically applies to usage over broad blocks of hours where the electricity price for each period is predetermined and constant. In this paper, we consider the TOU as simply the retail price paid by the end-users. The conventional building operator operates the HVAC systems to maintain a set temperature regardless of the electricity price. The conventional HVAC systems do not participate in the DR service and pay for electricity at the conventional price scheme, TOU (i.e., common price scheme at the distribution level) [24]. However, the proposed retail price on an hourly basis is applied to HVAC systems participating price-based DR service. We have assumed that the end-users are rational to adjust the power inputs of HVAC systems, so that they minimize the electricity bills as long as the thermal comforts are satisfied; i.e., the indoor temperatures are maintained within the acceptable ranges. In addition, smart grid technologies such as advanced metering infrastructures (AMIs) and behind-the-meter energy storage systems have been developing rapidly and continuously during the last decade. Moreover, significant attention has been given to developing energy-efficient devices, homes, or buildings and integrating a large number of renewable energy resources into power grid, to mitigate global warming and climate change. The development of such technologies and applications makes end-users become more aware of the rational use of electricity. For the same reason, the assumption on minimizing the end-users' electricity bills can be commonly made in many previous studies on price-based DR scheduling [15] and [25,26]. The EUC ensures that the proposed price is always less than the TOU to induce DR participation of HVAC systems in commercial buildings. Therefore, the 24-h operating costs of the price-sensitive HVAC units  As shown in Figure 1, the EUC determines the optimal profile of the retail price for the scheduling time (i.e., 24 h) to maximize its net profit, which is equivalent to the difference between the total revenue of selling electricity to HVAC systems and the total cost of purchasing electricity from the wholesale market. The EUC seeks to reduce peak loads to bypass costly additional facility (i.e., substations or lines) expansion. Therefore, the EUC changes the retail price to exploit the demand-side flexibility of the price-sensitive HVAC loads for peak load reduction. For the optimal profile of retail price, the building operators schedule the optimal power consumptions of HVAC systems so as to minimize their electricity bills while ensuring the thermal comforts of occupants.

Framework of the Proposed Pricing Strategy
For the proposed framework, a bi-level optimization problem is formulated to reflect the hierarchical decision models of the EUC and the consumers (i.e., leader and followers) in the upperand lower-decision levels, respectively. In this sense, the proposed pricing framework is related to the Stackelberg games in economics [27]. The optimization problem is transformed to an equivalent single-level problem using the Karush-Kuhn-Tucker (KKT) condition and the penalty method [9]. The EUC solves the single-level problem to maximize its profit, considering that the consumers minimize their electricity bills while satisfying the occupants' thermal comfort.
In addition, an accurate but simple thermal dynamic model (i.e., an experimental, data-driven model of the actual HVAC system) needs to be developed to estimate the indoor temperature variations for the power input changes of HVAC system, which requires a comprehensive understanding of the thermal dynamics of the target buildings under time-varying indoor and outdoor environmental conditions. We also need to pay careful attention when integrating the thermal dynamic model into the optimization problem, so that common solvers in the off-the-shelf software programs can be readily applied to find the optimal solution within a reasonable computational time. Comprehensive models of building thermal dynamics can be developed using a range of commercial software, for example, TRNSYS, Dymola/Modelica, and EnergyPlus. However, these comprehensive models are often too sophisticated and slow to be applied to the 24-h optimal scheduling problem where the EUC communicates with a number of buildings located throughout the distribution network. To resolve this challenge, the experimental data-driven models of a building room and HVAC system were adopted and integrated into the optimal scheduling problem in the proposed retail pricing strategy. As shown in Figure 1, the EUC determines the optimal profile of the retail price for the scheduling time (i.e., 24 h) to maximize its net profit, which is equivalent to the difference between the total revenue of selling electricity to HVAC systems and the total cost of purchasing electricity from the wholesale market. The EUC seeks to reduce peak loads to bypass costly additional facility (i.e., substations or lines) expansion. Therefore, the EUC changes the retail price to exploit the demand-side flexibility of the price-sensitive HVAC loads for peak load reduction. For the optimal profile of retail price, the building operators schedule the optimal power consumptions of HVAC systems so as to minimize their electricity bills while ensuring the thermal comforts of occupants.
For the proposed framework, a bi-level optimization problem is formulated to reflect the hierarchical decision models of the EUC and the consumers (i.e., leader and followers) in the upperand lower-decision levels, respectively. In this sense, the proposed pricing framework is related to the Stackelberg games in economics [27]. The optimization problem is transformed to an equivalent single-level problem using the Karush-Kuhn-Tucker (KKT) condition and the penalty method [9]. The EUC solves the single-level problem to maximize its profit, considering that the consumers minimize their electricity bills while satisfying the occupants' thermal comfort.
In addition, an accurate but simple thermal dynamic model (i.e., an experimental, data-driven model of the actual HVAC system) needs to be developed to estimate the indoor temperature variations for the power input changes of HVAC system, which requires a comprehensive understanding of the thermal dynamics of the target buildings under time-varying indoor and outdoor environmental conditions. We also need to pay careful attention when integrating the thermal dynamic model into the optimization problem, so that common solvers in the off-the-shelf software programs can be readily applied to find the optimal solution within a reasonable computational time. Comprehensive models of building thermal dynamics can be developed using a range of commercial software, for example, TRNSYS, Dymola/Modelica, and EnergyPlus. However, these comprehensive models are often too sophisticated and slow to be applied to the 24-h optimal scheduling problem where the EUC communicates with a number of buildings located throughout the distribution network. To resolve this Energies 2020, 13, 862 5 of 20 challenge, the experimental data-driven models of a building room and HVAC system were adopted and integrated into the optimal scheduling problem in the proposed retail pricing strategy.

Modeling the Thermal Response of HVAC Systems in Commercial Buildings
In this paper, the thermal dynamic response of an HVAC system in a small office building of United States was modeled based on the experimental building room thoroughly described in [13]. As shown in Figure 2, it is divided into a climate room and a test room, both in a larger laboratory room with a constant indoor temperature. The test room includes lights and general heat sources to simulate internal heat gains for a common office room. In [13], the thermal dynamic model of the test room was implemented using a regression-based inverse transfer function to characterize the linear dependencies on the external and internal environmental parameters during the past and current times as

Modeling the Thermal Response of HVAC Systems in Commercial Buildings
In this paper, the thermal dynamic response of an HVAC system in a small office building of United States was modeled based on the experimental building room thoroughly described in [13]. As shown in Figure 2, it is divided into a climate room and a test room, both in a larger laboratory room with a constant indoor temperature. The test room includes lights and general heat sources to simulate internal heat gains for a common office room. In [13], the thermal dynamic model of the test room was implemented using a regression-based inverse transfer function to characterize the linear dependencies on the external and internal environmental parameters during the past and current times as (1) In (1), Tht is the indoor temperature adjusted by HVAC system h, Tat is the adjacent room temperature, and Txt is the ambient temperature. In addition, Qht is the cooling rate provided by HVAC unit h, and Qct and Qrt are the heat gains from convective and radiative internal loads, respectively. The coefficients aht-kht depend on building orientation, structure, and materials. The adjacent and ambient temperatures, as well as the internal heat gains, can be measured and forecast using historical data acquired during normal building operation [14]. In this study, we focus on the relations between Tht and Qht by simplifying (1) to (2) where lht represents the building environmental parameters. The coefficients aht, sht, dht, fht, ght, and kht in (1) were estimated as shown in Table 1. The coefficients were obtained by developing and validating simulation models in [28].  Table 1. Coefficients for the inverse transfer function (ITF) model (1) that were experimentally obtained in [28].
The thermal dynamic model is incorporated particularly with the model of a variable speed heat pump (VSHP), which is a common example of an HVAC system. This leads the cooling rate provided by the VSHP unit h at time t to be represented [13] as: Figure 2. Experimental room setup for estimating the thermal dynamic response to the input power of an HVAC system.
In (1), T ht is the indoor temperature adjusted by HVAC system h, T at is the adjacent room temperature, and T xt is the ambient temperature. In addition, Q ht is the cooling rate provided by HVAC unit h, and Q ct and Q rt are the heat gains from convective and radiative internal loads, respectively. The coefficients a ht -k ht depend on building orientation, structure, and materials. The adjacent and ambient temperatures, as well as the internal heat gains, can be measured and forecast using historical data acquired during normal building operation [14]. In this study, we focus on the relations between T ht and Q ht by simplifying (1) to (2) where l ht represents the building environmental parameters. The coefficients a ht , s ht , d ht , f ht , g ht , and k ht in (1) were estimated as shown in Table 1. The coefficients were obtained by developing and validating simulation models in [28]. Table 1. Coefficients for the inverse transfer function (ITF) model (1) that were experimentally obtained in [28].
Energies 2020, 13, 862 6 of 20 Although the small test building, shown in Figure 2, has been used here for simplicity, the thermal dynamic modeling and, hence, the proposed optimal pricing strategy can be easily applied to large-scale, multi-zone commercial buildings using different thermal coefficients in (1) and (2) [13,29].
The thermal dynamic model is incorporated particularly with the model of a variable speed heat pump (VSHP), which is a common example of an HVAC system. This leads the cooling rate provided by the VSHP unit h at time t to be represented [13] as: In (3) and (4), δ mht is the m th power segment of the VSHP unit h at time t and z 0,ht -z 3,ht are the coefficients dependent on the steady-state operating conditions (i.e., δ mht and T xt ) of the VSHP unit h at time t. As discussed in [13], Q ht is also affected by the mixed-air temperature in the air-handling unit, which in general is maintained at 26 • C. We assumed that the return air temperature is kept the same as the indoor temperature if a supply air lower than 25 • C (i.e., maximum indoor temperature) is introduced into the room. Therefore, the return air temperature is lower than 26 • C. Under these conditions, it is assumed that the return air is mixed with a small amount of outdoor air in order to reduce the loss of cooling energy in buildings and to maintain indoor air quality above a certain level [30]. In [30], the mixed air temperature was set to constant 26 • C, but it showed the same effect even if the mixed air temperature is higher than 26 • C. By combining (2) and (3), the indoor temperature for the input power of the VSHP can be simply expressed as Figure 3 shows an example of the nonlinear curves obtained using (4). In Figure 3, T z,ht is defined as the indoor temperature at time t for Q hk = 0 for 0 ≤ k ≤ t; i.e., the HVAC system is turned off during the period. Furthermore, F mhtj is the linear gradient of T ht at time t = j for the m th power segment of HVAC system h at time t = k (i.e., k ≤ j). N L represents the number of piecewise linear blocks. The nonlinear thermal dynamic response of the building room to the input power variation of HVAC unit can be approximated to (5) using a piecewise linearization approach where the operating range from 0 to P h,max of HVAC system is divided into N L segments.
Energies 2020, 13, x FOR PEER REVIEW 6 of 19 (3) In (3) and (4), δmht is the mth power segment of the VSHP unit h at time t and z0,ht-z3,ht are the coefficients dependent on the steady-state operating conditions (i.e., δmht and Txt) of the VSHP unit h at time t. As discussed in [13], Qht is also affected by the mixed-air temperature in the air-handling unit, which in general is maintained at 26 °C. We assumed that the return air temperature is kept the same as the indoor temperature if a supply air lower than 25 °C (i.e., maximum indoor temperature) is introduced into the room. Therefore, the return air temperature is lower than 26 °C. Under these conditions, it is assumed that the return air is mixed with a small amount of outdoor air in order to reduce the loss of cooling energy in buildings and to maintain indoor air quality above a certain level [30]. In [30], the mixed air temperature was set to constant 26 °C, but it showed the same effect even if the mixed air temperature is higher than 26 °C. By combining (2) and (3), the indoor temperature for the input power of the VSHP can be simply expressed as (4) Figure 3 shows an example of the nonlinear curves obtained using (4). In Figure 3, Tz,ht is defined as the indoor temperature at time t for Qhk = 0 for 0 ≤ k ≤ t; i.e., the HVAC system is turned off during the period. Furthermore, Fmhtj is the linear gradient of Tht at time t = j for the mth power segment of HVAC system h at time t = k (i.e., k ≤ j). NL represents the number of piecewise linear blocks. The nonlinear thermal dynamic response of the building room to the input power variation of HVAC unit can be approximated to (5) using a piecewise linearization approach where the operating range from 0 to Ph,max of HVAC system is divided into NL segments. (5) This approach aims to develop a set of linear constraints on the indoor temperature variation, as discussed in Section 3, which can be applied to the KKT conditions. Consequently, the bi-level decision model for the proposed optimal pricing strategy can be converted into the single-level decision model.

Optimization Problem Formulation
For the proposed framework discussed in Section 2, the optimal electricity price schedule and corresponding input power profile of HVAC units can be determined by solving the upper-level problem and lower-level problem, respectively.  This approach aims to develop a set of linear constraints on the indoor temperature variation, as discussed in Section 3, which can be applied to the KKT conditions. Consequently, the bi-level decision model for the proposed optimal pricing strategy can be converted into the single-level decision model.

Optimization Problem Formulation
For the proposed framework discussed in Section 2, the optimal electricity price schedule and corresponding input power profile of HVAC units can be determined by solving the upper-level problem and lower-level problem, respectively.

Upper-Level Problem
The EUC determines the optimal retail price C t to maximize its profit J DP by solving the upper-level problem (6)- (10). In (6), the EUC's profit J DP is calculated as the difference between the revenue of selling the electricity to the price-sensitive HVAC systems at C t and the cost of purchasing the electricity at M t from the wholesale market. N T is the number of scheduling time intervals in a day. The conventional HVAC systems were assumed to consume pre-determined input power P C,ht for 5 ≤ t ≤ 19 to maintain the indoor temperature at 22.5 • C during 8 ≤ t ≤ 19. The EUC's profit J DC = t (U t -M t ) × h, P C,ht in business with the conventional HVAC systems is constant, because the TOU price U t is also determined in advance. It can then be added to (6), resulting in the total profit of the EUC J D = J DP + J DC .
In (7), the ratio of the power consumption of HVAC systems in the commercial buildings to the total power consumption in the distribution network is set to R H and used to determine the total number of commercial HVAC systems in the distribution network. N H is the number of the types of HVAC systems and N HVAC is the total number of HVAC systems in the distribution system divided by N H .
The EUCs are regulated not to maximize their profit by manipulating the prices; therefore, we set (8) to prevent EUC from setting its price to maximize profit. In addition, (8) represents the upper and lower limits of C t at time t, which are equivalent to U t and M t , respectively. This allows consumers to reduce their electricity bills via the proposed retail pricing scheme, and, consequently, the EUC can attract a number of HVAC loads to its optimal price-based DR strategy, procuring substantial load flexibility. The constraints (8) can also reduce the exposure consumers have to highly volatile prices and alleviate the EUC's market power by preventing the EUC from increasing its profit excessively, which appears to be unfair to consumers. It was assumed in [31] that the EUC and consumers reach such an agreement to bring (8) into effect. In (9), P R,t is defined as the remainder load other than the load used in conventional HVAC systems. N HVAC · h m δ mht + P R,t is established to estimate the peak limit P peak because of price-sensitive HVAC loads in (10).

Lower-Level Problem
For the price-sensitive HVAC units, the optimal power inputs m δ mht for all h are scheduled to minimize their operating cost J BP , given C t , by solving the lower-level problem (11)- (17). The operations of the conventional HVAC units result in the constant operating cost J BC = t U t × h P C,ht , which can be simply added to (11) as J B = J BP + J BC without affecting the optimal schedules of δ mht . To complete the piecewise linear approximation, (12)- (14) are established to describe the boundary conditions on δ mht , as shown in Figure 3, using the binary variables b mht . For example, for b 2ht = 1, b 3ht = 0, and b 4ht = 0, δ 3ht can be scheduled to the value between 0 and δ 3h,max , whereas δ 4ht is forced to be zero. Moreover, (15) specifies that the scheduled power inputs m δ mht of HVAC systems need to be less than their rated power P h,max , and (16) specifies the limits D L and D U on the positive and negative ramp rates of m δ mht , respectively, during the unit time interval ∆t unit = 1 h. In general, D U needs to be carefully determined by considering practical operation of heat pumps (HPs) to prevent severe operational stress on HP compressors. Furthermore, the inequality constraints (17) specify that the indoor temperature T ht should be maintained within an acceptable range between T ht,min and T ht,max so as to satisfy the occupants' thermal comfort. Note that T ht is approximated using a piecewise linear equation (5) for the application of the KKT conditions.

Test System and Simulation Conditions
In the case studies, the effect of the proposed pricing strategy on the peak load reduction are explicitly investigated for the cases of conventional load demand and duck curve, and the data in Figure 4b-d are from United States (i.e., test site). A 24-h load demand curve during summertime was adopted from [34] and scaled down, so that the distribution network has the peak load demand of 7.04 MW at t = 17 h, as shown by the blue line in Figure 4a. The PV production time-series (i.e., red line in Figure 4a) was acquired from [35] and scaled up to 20% of the peak of the conventional load curve (i.e., the blue line in Figure 4a). The load curve that changed with the introduction of PV is represented by the green line shown in Figure 4a. The peak load of the green line is 6.61 MW at t = 20 h. The 24-h curve of the ambient air temperature [36], shown in Figure 4b, is similar to the conventional load demand curve. In particular, both the maximum temperature and the peak load took place at t = 17 h. In addition, Figure 4c indicates the curves of the wholesale price M t [37] and retail price U t Energies 2020, 13, 862 10 of 20 for the TOU scheme [38]. It was assumed that the EUC was informed of M t one day ahead to schedule the optimal C t considering the price-based DR of HVAC units. For the optimal C t , the input power of HVAC systems was then determined, given the N H profiles of the building thermal loads, which were developed as shown in Figure 4d based on the real data in [30] and [39]. Specifically, the VSHP load profiles in a real commercial building were measured and then scaled down with the different ratios to determine the rated power inputs of the VSHPs. For the rated inputs, the corresponding cooling rates Q ht were then estimated using (3), and the scaled fractions of Q ht were applied to the control of T ht . As shown in Figure 4d, the N H type profiles of thermal energy loads Q ct and Q rt were then determined for the test buildings. In this paper, we assumed that the internal heat gain could be forecasted as Figure 4d. Based on the real profiles of m δ mht , Q ht , and T xt , for simplicity, it was assumed that the group of buildings has the same thermal load profile. Figure 4d shows that the thermal loads were increased at t = 7 h (when people started coming to work) and were maintained at high levels until t = 19 h, when the workers left the buildings. Table 2 lists the rated values of m δ mht and Q ht of HVAC units and maximum value of (Q ct + Q rt ) per building room.   As mentioned in Section 3, NH types of HVAC systems were used and the total number of HVAC systems included in the distribution network was NH × NHVAC. In this paper, NH was set to 15 as shown in Table 2. As mentioned in Section 1, more than 13% of the total energy usage is consumed by commercial HVAC systems, so the case studies were simulated with RH = 0.13. In the case of the conventional load profiles, NHVAC is set to 6 by (7), and in the case of duck curve, NHVAC is set to 4.6 by (7). Figure 5a shows that the input power of the conventional HVAC systems varied over the scheduling time period to maintain Tht at Tset = 22.5 °C during the working hours (i.e., 8 h ≤ t ≤ 19 h), as shown in Figure 5b, given the time-varying thermal conditions Txt, Qct, and Qrt of the buildings. It   As mentioned in Section 3, N H types of HVAC systems were used and the total number of HVAC systems included in the distribution network was N H × N HVAC . In this paper, N H was set to 15 as shown in Table 2. As mentioned in Section 1, more than 13% of the total energy usage is consumed by commercial HVAC systems, so the case studies were simulated with R H = 0.13. In the case of the conventional load profiles, N HVAC is set to 6 by (7), and in the case of duck curve, N HVAC is set to 4.6 by (7). Figure 5a shows that the input power of the conventional HVAC systems varied over the scheduling time period to maintain T ht at T set = 22.5 • C during the working hours (i.e., 8 h ≤ t ≤ 19 h), as shown in Figure 5b, given the time-varying thermal conditions T xt , Q ct , and Q rt of the buildings. It was assumed that the conventional HVAC systems are operated during 5 h ≤ t ≤ 19 h. As mentioned in Section 3, NH types of HVAC systems were used and the total number of HVAC systems included in the distribution network was NH × NHVAC. In this paper, NH was set to 15 as shown in Table 2. As mentioned in Section 1, more than 13% of the total energy usage is consumed by commercial HVAC systems, so the case studies were simulated with RH = 0.13. In the case of the conventional load profiles, NHVAC is set to 6 by (7), and in the case of duck curve, NHVAC is set to 4.6 by (7). Figure 5a shows that the input power of the conventional HVAC systems varied over the scheduling time period to maintain Tht at Tset = 22.5 °C during the working hours (i.e., 8 h ≤ t ≤ 19 h), as shown in Figure 5b, given the time-varying thermal conditions Txt, Qct, and Qrt of the buildings. It was assumed that the conventional HVAC systems are operated during 5 h ≤ t ≤ 19 h. In addition, the load factor and load leveling index are introduced to analyze the effect of the proposed pricing strategy on the peak load reduction. In [40], the load factor can be defined as the ratio of the average demand to the peak load as: The load leveling index is the MSE with respect to the reference value as [41]: Figure 6a shows the optimal schedules of C t for the proposed pricing strategy considering peak load reduction applied to conventional load curve, and Figure 6b,d show the corresponding data on the total power inputs of the price-elastic HVAC units and the indoor temperatures in the building rooms. In Figure 6c, the red line shows the total load curve when HVAC systems are operated in the conventional manner, and the blue line shows the total load curve when HVAC systems are operated by the price sensitive method. Figure 6a shows the optimal schedules of Ct for the proposed pricing strategy considering peak load reduction applied to conventional load curve, and Figure 6b,d show the corresponding data on the total power inputs of the price-elastic HVAC units and the indoor temperatures in the building rooms. In Figure 6c, the red line shows the total load curve when HVAC systems are operated in the conventional manner, and the blue line shows the total load curve when HVAC systems are operated by the price sensitive method. In Figure 6a, the EUC maintained Ct as sufficiently low in the early morning (i.e., 1 h ≤ t ≤ 8 h), taking into account that the end-users aim to minimize the operating cost of their HVAC systems. This led to the pre-cooling of HVAC systems, reducing the indoor temperature before people come to work. Specifically, the large input power of the HVAC systems was scheduled during the time period of 3 h ≤ t ≤ 7 h, even though the ambient temperature and thermal energy loads were maintained high during 8 h ≤ t ≤ 19 h. Owing to the low retail price in the early morning, HVAC loads were shifted away from on-peak hours to off-peak hours (i.e., early morning), leading to the precooling operation and a reduction in the EUC's profit and building operators' operation cost. In Figure 6d, the pre-cooled indoor temperature increased gradually from 20 °C at t = 8 h to 25 °C at t = In Figure 6a, the EUC maintained C t as sufficiently low in the early morning (i.e., 1 h ≤ t ≤ 8 h), taking into account that the end-users aim to minimize the operating cost of their HVAC systems. This led to the pre-cooling of HVAC systems, reducing the indoor temperature before people come to work. Specifically, the large input power of the HVAC systems was scheduled during the time period of 3 h ≤ t ≤ 7 h, even though the ambient temperature and thermal energy loads were maintained high during 8 h ≤ t ≤ 19 h. Owing to the low retail price in the early morning, HVAC loads were shifted away from on-peak hours to off-peak hours (i.e., early morning), leading to the pre-cooling operation and a reduction in the EUC's profit and building operators' operation cost. In Figure 6d, the pre-cooled indoor temperature increased gradually from 20 • C at t = 8 h to 25 • C at t = 12 h, as the ambient temperature and internal heat gains increased. Because of the limited thermal capacity in the building structures, HVAC loads started increasing at t = 13 h for satisfying occupants' comforts (i.e., temperature constraints). During 11 h ≤ t ≤ 19 h, the EUC sets C t to the upper limit to compensate for the reduction in the EUC's profit, resulting from the low value of C t for the pre-cooling. Specially, the EUC sets C t to the upper limit to lower power consumption during on-peak hours (i.e., 16 h ≤ t ≤ 19 h). To prevent excessive compensation of the EUC's profit, the upper limit was set as the TOU price U t , as shown in (8), mitigating the retail price volatility and ensuring the financial benefit of the consumers participating in the proposed pricing scheme. In Figure 6c, the peak load of the total load containing the conventional HVAC system is 7.04 MW, while the peak load of the total load containing the proposed price-sensitive HVAC system is 6.75 MW and the reduction rate is approximately 4.18%. In the proposed method, the peak of the total load curve takes place at t = 20 h. The peak load reduction rate cannot be increased further because HVAC system does not need to operate for 20 h ≤ t ≤ 24 h owing to the mitigation of the indoor temperature constraints. However, the difference between the peak and minimum load values shows that the proposed method does not improve in terms of load leveling with 2.33 MW for the conventional HVAC system and 2.32 MW for the price-sensitive HVAC system. Table 3 summarizes the J BP and J BC for individual HVAC units for the proposed pricing scheme and the TOU pricing scheme, respectively, for peak load reduction of the conventional load curve. The 24-h operating costs J BP of the price-sensitive HVAC units are cheaper than those of the conventional units. This demonstrates that using the proposed pricing scheme, the EUC can effectively attract more HVAC units to its DR program. As shown in Figure 7, the optimal schedules of C t and m δ mht were determined by proposed pricing strategy considering peak load reduction applied to duck curve. The peak load of the duck curve occurs at t = 20 h, as shown in Figure 7c. The commercial HVAC system does not need to be operated for 20 h ≤ t ≤ 24 h owing to mitigation of the indoor temperature conditions. As a result, the commercial HVAC system cannot affect the load curve after 20 h. In the proposed method (i.e., blue line) in Figure 7c, the difference between the peak load and the minimum value of the load curve is larger than the conventional strategy (i.e., red line). In other words, even if the EUC utilizes the commercial HVAC system as a DR resource, the peak load of duck curve cannot be reduced.

Case Study A: Two Types of Load Profiles
improve in terms of load leveling with 2.33 MW for the conventional HVAC system and 2.32 MW for the price-sensitive HVAC system. Table 3 summarizes the JBP and JBC for individual HVAC units for the proposed pricing scheme and the TOU pricing scheme, respectively, for peak load reduction of the conventional load curve. The 24-h operating costs JBP of the price-sensitive HVAC units are cheaper than those of the conventional units. This demonstrates that using the proposed pricing scheme, the EUC can effectively attract more HVAC units to its DR program. As shown in Figure 7, the optimal schedules of Ct and ∑m δmht were determined by proposed pricing strategy considering peak load reduction applied to duck curve. The peak load of the duck curve occurs at t = 20 h, as shown in Figure 7c. The commercial HVAC system does not need to be operated for 20 h ≤ t ≤ 24 h owing to mitigation of the indoor temperature conditions. As a result, the commercial HVAC system cannot affect the load curve after 20 h. In the proposed method (i.e., blue line) in Figure 7c, the difference between the peak load and the minimum value of the load curve is larger than the conventional strategy (i.e., red line). In other words, even if the EUC utilizes the commercial HVAC system as a DR resource, the peak load of duck curve cannot be reduced.

Case Study B: Purpose of Demand Response
We demonstrated that the proposed pricing strategy could be applied to other purpose such as load leveling besides peak load reduction. The simple constraint (i.e., P min < N HVAC · m δ mht + P R,t ) was added to (10) to show that price-sensitive HVAC systems can have positive effects even in terms of load leveling, and, the simulation results of applying the proposed pricing strategy to the conventional load curve are shown in Figure 8. To satisfy the added constraint, the EUC lowered C t for 6 h ≤ t ≤ 13 h to induce the power consumption of HVAC system. In particular, the blue line in Figure 8a had a minimum value at t = 6 h, and the EUC set C t close to the lower limit to increase power consumption during this time period. For HVAC system, the input power of the previous time zone affects the indoor temperature of the subsequent time zone, as shown in (4). Therefore, an increase in power consumption for 6 h ≤ t ≤ 8 h could violate the indoor temperature conditions at t = 8 h; this problem could be solved by increasing C t for 1 h ≤ t ≤ 3 h to reduce the power consumption for 1 h ≤ t ≤ 3 h. The difference between the peak load and minimum values is 1.51 MW, a significant reduction from 2.33 MW for conventional HVAC systems. power consumption during this time period. For HVAC system, the input power of the previous time zone affects the indoor temperature of the subsequent time zone, as shown in (4). Therefore, an increase in power consumption for 6 h ≤ t ≤ 8 h could violate the indoor temperature conditions at t = 8 h; this problem could be solved by increasing Ct for 1 h ≤ t ≤ 3 h to reduce the power consumption for 1 h ≤ t ≤ 3 h. The difference between the peak load and minimum values is 1.51 MW, a significant reduction from 2.33 MW for conventional HVAC systems.  Figure 9 shows the simulation results of applying the proposed pricing strategy considering load leveling to the duck curve. The blue line in Figure 9c shows that the difference between the peak load and the minimum load is similar to that of the conventional method (i.e., red line); however, the blue line also shows a more rapid change in power consumption for 16 h ≤ t ≤ 20 h. In other words, the total load profile (i.e., blue line) of the price-sensitive HVAC systems requires larger ramp rate of generators for 16 h ≤ t ≤ 20 h. As a result, even if the EUC utilizes the commercial HVAC system as a DR resource, the duck curve cannot be improved in terms of load leveling, as shown in Figure 9c.  Figure 9 shows the simulation results of applying the proposed pricing strategy considering load leveling to the duck curve. The blue line in Figure 9c shows that the difference between the peak load and the minimum load is similar to that of the conventional method (i.e., red line); however, the blue line also shows a more rapid change in power consumption for 16 h ≤ t ≤ 20 h. In other words, the total load profile (i.e., blue line) of the price-sensitive HVAC systems requires larger ramp rate of generators for 16 h ≤ t ≤ 20 h. As a result, even if the EUC utilizes the commercial HVAC system as a DR resource, the duck curve cannot be improved in terms of load leveling, as shown in Figure 9c.   Figure 9 shows the simulation results of applying the proposed pricing strategy considering the load leveling to the duck curve. The blue line in Figure 9c shows that the difference between the peak load and the minimum load is similar to that of the conventional method (i.e., red line); however, the blue line also shows a more rapid change in power consumption for 16 h ≤ t ≤ 20 h. In other words, the total load profile (i.e., blue line) of the price-sensitive HVAC systems requires larger ramp rate of generators for 16 h ≤ t ≤ 20 h. As a result, even if the EUC utilizes the commercial HVAC system as a DR resource, the duck curve cannot be improved in terms of load leveling, as shown in Figure 9c.

Case Study C: Parameter Sensitivity
In Section 4.4, we performed the simulation with several value of TOU pricing to demonstrate the applicability of the proposed pricing strategy, because the TOU pricing varies depending on the   Figure 9 shows the simulation results of applying the proposed pricing strategy considering the load leveling to the duck curve. The blue line in Figure 9c shows that the difference between the peak load and the minimum load is similar to that of the conventional method (i.e., red line); however, the blue line also shows a more rapid change in power consumption for 16 h ≤ t ≤ 20 h. In other words, the total load profile (i.e., blue line) of the price-sensitive HVAC systems requires larger ramp rate of generators for 16 h ≤ t ≤ 20 h. As a result, even if the EUC utilizes the commercial HVAC system as a DR resource, the duck curve cannot be improved in terms of load leveling, as shown in Figure 9c.

Case Study C: Parameter Sensitivity
In Section 4.4, we performed the simulation with several value of TOU pricing to demonstrate the applicability of the proposed pricing strategy, because the TOU pricing varies depending on the country or month. Figure 10 shows the effect of the price cap U t on the optimal profile of the C t . Figure 10a,b show that as U t max increased, C t increased accordingly; however, the overall profiles of C t and, consequently, the overall profiles of the EUC's profit remained similar. Note that in Figure 10c, the difference value between U t max and U t min were maintained at the same values. Figure 10b,d show the optimal profiles of C t for the different profiles of U t shown in Figure 10a,c, respectively. The C t profiles, apart from the magnitudes, look similar for all cases particularly during the time period from 1 h to 19 h; note that HVAC systems were turned off in the period 20 h ≤ t ≤ 24 h because of low occupancy.

Discussion
The results for various case studies are summarized in Tables 4 and 5. In Table 4, the proposed method reduced the peak load of the conventional load curve by 4.12% compared to the conventional method. In terms of load leveling, the proposed method, which only considers the peak load reduction, reduced the difference between the peak load and the minimum load by 0.43% compared to the conventional method, and reduced by 35.19% when the constraint for the load leveling is added. The load factor has similar values in the cases of (a)-(c). The proposed method did not improve the load curve in terms of load factor. However, the proposed method reduced the load leveling index by 39.59% in the case of (b) and 62.80% in the case of (c), compared to the conventional method. The load leveling index is further reduced in the case of (c) where load leveling constraints are added compared to the case of (b) with only constraints of peak load reduction. The proposed method effectively achieved peak load reduction and load leveling. These allow for the postponement of investments in grid upgrades or in new substation capacity. The EUC's profit for commercial HVAC systems (i.e., JDP) with the constraint for the peak load reduction is 56.51% less than that of the conventional method. The JDP is reduced by 67.38% with the addition of load leveling constraint. In other words, by load leveling and reducing the peak load, the profit for commercial HVAC systems is greatly reduced. While the JDP is greatly reduced, the total profit of the distribution network decreased by 6.81% (peak load reduction) and 7.94% (load leveling), respectively. Although this may initially limit the profit of the EUC, the reduction in profit can be sufficiently compensated for by the savings in the network operating costs that result from the improved load flexibility, e.g., a decrease in the installation costs of new substations and lines.

Discussion
The results for various case studies are summarized in Tables 4 and 5. In Table 4, the proposed method reduced the peak load of the conventional load curve by 4.12% compared to the conventional method. In terms of load leveling, the proposed method, which only considers the peak load reduction, reduced the difference between the peak load and the minimum load by 0.43% compared to the conventional method, and reduced by 35.19% when the constraint for the load leveling is added. The load factor has similar values in the cases of (a)-(c). The proposed method did not improve the load curve in terms of load factor. However, the proposed method reduced the load leveling index by 39.59% in the case of (b) and 62.80% in the case of (c), compared to the conventional method. The load leveling index is further reduced in the case of (c) where load leveling constraints are added compared to the case of (b) with only constraints of peak load reduction. The proposed method effectively achieved peak load reduction and load leveling. These allow for the postponement of investments in grid upgrades or in new substation capacity. The EUC's profit for commercial HVAC systems (i.e., J DP ) with the constraint for the peak load reduction is 56.51% less than that of the conventional method. The J DP is reduced by 67.38% with the addition of load leveling constraint. In other words, by load leveling and reducing the peak load, the profit for commercial HVAC systems is greatly reduced. While the J DP is greatly reduced, the total profit of the distribution network decreased by 6.81% (peak load reduction) and 7.94% (load leveling), respectively. Although this may initially limit the profit of the EUC, the reduction in profit can be sufficiently compensated for by the savings in the network operating costs that result from the improved load flexibility, e.g., a decrease in the installation costs of new substations and lines.   Table 5 reveals the proposed method did not reduce the peak load of the duck curve. The proposed method for peak load reduction increases the difference between peak and minimum value by −19.05%. Applying the constraint for the load leveling increases the difference between the peak and minimum value by −2.65%. The proposed strategy reduced the load factor in the cases (b) and (c), and has an adverse effect on the load curve in terms of load factor. In addition, the proposed method increased the load leveling index by −29.96% in the case of (b) and −12.12% in the case of (c), compared to the conventional method. However, in the proposed method, the retail price is lower than the TOU price in the conventional method, which reduces the J DP . As a result, the peak load reduction or load leveling of the duck curve cannot be achieved by using the DR of the commercial HVAC systems. This will be addressed in future studies by inducing DR participation of DR resources such as HVAC systems in residential buildings.

Conclusions
This paper proposed an optimal retail pricing scheme to support the EUC and building operators in determining the day-ahead optimal retail price and the corresponding input power of price-elastic HVAC units, respectively, reducing the peak load in the distribution network. A bi-level decision model was adopted to develop the proposed scheme, i.e., in the upper level, the EUC determines the optimal price to maximize its profit and in the lower level, the building operators schedule the optimal input power of HVAC units to minimize their operating costs. The proposed bi-level optimal pricing model was applied to two types of load curves using the experimental data-driven model of an HVAC system in commercial buildings. The optimal retail price was determined to reduce peak load demand, avoiding the cost for reinforcing power equipment and distribution lines. The upper bound of the retail price was set using the TOU pricing rates to achieve the low price volatility in the DR program. The nonlinear thermal dynamic response of an experimental building room to the input power of the VSHP was modeled using a set of piecewise linear equations, which enabled the application of the KKT conditions and the transformation of the bi-level decision model to the single-level model. The results of the case studies demonstrated the effectiveness of the proposed strategy with the following main features:

•
The EUC induced the pre-cooling of HVAC systems, considering that the consumers minimize the operating cost of HVAC units. To compensate for the reduction of the EUC's profit because of the pre-cooling, the EUC determined high retail price during on-peak times. In addition, the building operators scheduled small input power of HVAC system during the on-peak hours to minimize the electricity bills.

•
For the conventional load curve, the peak load could be successfully reduced using the pre-cooling of HVAC units. The corresponding reduction of the EUC's profit can be justified by the opportunity cost of installing an additional line or substation and the reduction in consumers' electricity bills.

•
For the duck curve, as participation in DR of the commercial HVAC systems does not result in peak load reduction or load leveling, residential HVAC systems will also need to be addressed by inducing participation in DR. In addition, we will research additional case studies on the application of the proposed pricing strategy to large-scale, multi-zone buildings.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: A. Acronyms: