A Kullback-Leibler Divergence-based Distributionally Robust Optimization Model for Heat Pump Day-ahead Operational Schedule in Distribution Networks
1 A bstract—For its high coefficient of performance and zero local emissions, the heat pump (HP) has recently become popular in North Europe and China. However, the integration of HPs may aggravate the daily peak-valley gap in distribution networks significantly. In this paper, we describe a distributionally robust optimization (DRO)-based heat pump day-ahead operational schedule model (HP-DOSM) to shave the peak and reduce residents’ costs under time-of-use. The ambiguity set of the DRO is constructed using Kullback–Leibler divergence with a nominal distribution. This model can well capture the uncertainties of weather, photovoltaic, and load prediction errors. Moreover, this DRO based HP-DOSM can be transformed into a tractable deterministic model. Compared with robust optimization (RO) models, our model is less conservative since more statistical information of the uncertainties is utilized. Numerical tests were conducted to demonstrate its performance, compared with the RO model via Monte Carlo simulations. Index Terms —Distribution network, distributionally robust optimization, day-ahead operational schedule, Kullback–Leibler divergence N OMENCLATURE N Total number of houses equipped with heat pump and water tank H Time horizon in 1 day , k t T Indoor air temperature of the k th house at time t , , w k t T Tank water temperature of the k th house at time t , out t T Outdoor temperature at time t k T Vector ,0 ,1 ,
Tk k k H
T T T , w k T Vector , ,0 , ,1 , ,
Tw k w k w k H
T T T out T Vector ,0 ,1 ,
Tout out out H
T T T k T Upper temperature boundary of the k th house k T Lower temperature boundary of the k th house k R Thermal resistance of the k th house Manuscript received XXX, 2017. This work was supported in part by the National Science Foundation of China (51477083), in part by Key Science and Technology Project of State Grid Corporation of China (Grant. 5202011600U4). The authors are with the State Key Laboratory of Power Systems, Department of Electrical Engineering, Tsinghua University, Beijing 100084, China (e-mail: [email protected]). k C Thermal capacitor of the k th house , w k R Thermal resistance of the k th water tank , w k C Thermal capacitor of the k th water tank , k t x ON/OFF (1/0) state of the k th heat pump k x Vector ,0 ,1 ,
Tk k k H x x x , HP k P Rated power of the k th heat pump , COP k Coefficient of performance of the k th HP w h k Water tank to house heating efficiency of the k th house Z Zone of local transformer and feeder lines , ,
Load k t P Load of the k th house except heat pump at time t , , PV k t P Photovoltaic power of the k th house at time t , Ztrans t P Transformer power of zone Z at time t , P t Forecast errors of power at time t , T t Forecast errors of outdoor temperature at time t t P Distribution of power forecast error at time t t Ambiguity set of t P t P Center/nominal distribution of t t Q Distribution of outdoor temperature forecast error at time t Q t Ambiguity set of t Q t Q Center/nominal distribution of Q t t Farthest KL divergence in the ambiguity set of power forecast error at time t t Farthest KL divergence in the ambiguity set of outdoor temperature forecast error at time t I. I NTRODUCTION eat pumps (HPs) have seen a revival in recent years, in which their potential role with renewable energy sources is being investigated. HP technologies have experienced great improvements in operating stability and have subzero coefficient of performance [1], providing confidence for consumers in their reliability. Unsurprisingly, more than 750,000 units were sold in the EU-20, reducing greenhouse gas emissions by 6.8 Mt in 2010 [2]. However, HPs have placed a considerable load on power systems, although recognized as a promising resource for demand-side management in using renewable energy [3]. In recent years, a policy named ‘Coal to Electricity’ has been promulgated in Beijing, to replace coal stoves in the nearby countryside with HPs for heating, which may be beneficial in reducing the gray smog sky over Beijing [4]. By the end of November 2016,
A Kullback–Leibler Divergence-based Distributionally Robust Optimization Model f or Heat Pump Day-ahead Operational Schedule in Distribution Networks
Zihao Li, Wenchuan Wu , Senior Member, IEEE,
Boming Zhang,
Fellow, IEEE H D ETERMINISTIC HEAT PUMP D AY -A HEAD O PERATIONAL S CHEDULE MODEL
Figure 1 shows a local distribution network integrated with photovoltaic panels and consumers with air-source HPs. For comfort considerations, the HP first heats the water in a tank to supply heat for the house. The tank is installed indoors and links to a finned tube with constant water flux.
Distribution network
Figure 1. Local distribution network
We apply a one-order equivalent thermal RC model to simulate the house and the water tank. The RC model is formulated as follows, ,, , , ,, k w kk COP k HP k w k w kw k
T Tx P C TR (1) ,2 , , w k kout k w h k k kk w k
T TT T C TR R (2) where , k k R C and , , , w k w k R C are the 1 st -order thermal parameters of the house and water tank, respectively, of the th k house, and they can be identified from historical measurements for real application. k T , , w k T , and out T are the indoor air, water tank, and outdoor air temperatures, respectively. , HP k P , , COP k , w h k , and k x are the rated power, constant HP to water coefficient of performance (COP), constant water-to-house efficiency, and the ON-OFF 3 state of the th k house’s HP. The difference schemes of (1) and (2) are ( ) ( ){ [ ]} k k w h k k w k t k tk t out tw k tw h k k w k t k t R Ck t out t w k R T TT TR R T TT T eR (3) , , , , 1 , , , , ,, , , , , , , ( ( )) w k w k w k t w k COP k HP k k t k t tR Cw k t w k COP k HP k k t k t
T R P x TT R P x T e (4) where t is the time interval. The HP operates in two alternative modes: OFF and ON. Thus, the power consumed by a HP is a function of its operational state: ,, , , , , , k tHP k t k t HP k k t k t for state ONP x P here x for state OFF , (5) The objective function of the HP day-ahead operational schedule we propose includes shaving the peak load at the local transformer and minimizing the electricity cost: minimize max , , Z t k t HP kZ Trans k Z t
P x P (6) where max Z P represents one-day’s maximum power load on the transformer Z , is a penalty weight for peak shaving, called the peak cost factor [22], and t represents the time of use (TOU) electricity price. A DC power flow model is used for the distribution network and the power balance constraint for the transformer of zone Z is: , , , , , , , , ZTrans t k t HP k t Load k t PV k tk Z k Z k Z
P x P P P (7) where , ,
Load k t P and , , PV k t P are the residential loads (except the HPs) and the PV output. The capacity constraints of the transformer are max , , =1,2,..., Z ZTrans t
P P t H (8) max _ Z Ztrans capacity
P P (9) where _ Ztrans capacity P is the rated capacity of the local transformer Z and H is the total schedule time periods of 1 day. Thermal constraints include (3), (4), indoor temperature comfort constraints: , k k t k T T T (10) The initial temperature settings are , 0 , 0, k t k t set
T T (11) , , 0 , , 0, w k t w k t set T T (12) The daily water temperature hold is , , , , 0 w k t H w k t T T (13) and the HP ON-OFF switching time-delay limitations are , , 1 , 1 k t k t k t x x x (14) where index t represents the schedule time, , 0, k t set T and , , 0, w k t set T are the initial temperatures of the indoor air and the tank water, and , k k T T are the indoor air temperature comfort range bounds. The switching time-delay constraints indicate that the HP needs to keep the same mode for at least two time intervals. Then, the deterministic HP-DOSM can be formulated as:
Model 1: Deterministic heat pump day-ahead operational schedule model min (6) . . s t (3)(4),(7)-(14) , {0,1} k t x ,
1, , , 1,2, , , t H k N Z
With the given prediction data, this MILP model can solved efficiently. However, the schedule decision based on bias prediction may lead to infeasibility in real operation. Accordingly, we introduce a KL divergence-based DRO in the next section, considering the uncertainties of forecast errors in outdoor temperature, load prediction, and PV prediction. III. KL DIVERGENCE - BASED
DRO M ODEL FOR HP D AY -A HEAD O PERATIONAL S CHEDULE
In this section, we first present the KL divergence-based DRO (
KL-DRO ) formulation to minimize total cost given the expectation constraints. Then, this
KL-DRO model is transformed to a deterministic MILP. Finally, we discuss how to use historical data to build an ambiguity set for the
KL-DRO model as well as the risk meaning of the radius distance of the ambiguity set. There are several statistical distance functions, such as the Prokhorov metric, the KL divergence, and the Wasserstein metric, which show common properties and can co-transfer in some situations [17]. Here, the KL divergence is adopted and is defined as ( )( ( ) || ( )) ( )ln ( ) KL pD p q p dq (15) where , p q are distribution functions in measure space . A. KL-DRO model
To capture the uncertainties in prediction errors, we reform the deterministic constraints, (8), (10), and (13), as expected constraints under the worst distribution in the ambiguity sets: i.e., constraints (7) and (8) are reformed as ,, , max max [ ] 0, 1,2, ,
Zsum tP k t HP kP k Z
E x P P P t H (16) Here, , , , , , Zsum t Load k t PV k tk Z k Z
P P P is the summation of load demand and PV output in zone Z. , Zsum t P is further expressed as , , , Z Zsum t sum t P t
P P (17) where , Zsum t P is the prediction value and , P t represents the prediction error. The distribution function t P of , P t belongs to the ambiguity set t , defined as { | ( || ) } t t t KL t t t P D P P (18) where t is the farthest KL divergence with the nominal distribution t P , also referred to as the radius. 4 In a similar way, we can also reformulate constraints (10) and (13) as Q max [ ] 0 outQ k k k k kQ E J x K T l T (19) Q max [ ] 0 outQ k k k k kQ E J x K T l T (20) max ( ) 0 outQ w k k w k w k w k t setQ E Z M x N T p T (21) where k x , k T , w k T , out T are the vectors explained in the nomenclature, Z , and the equalities (3) and (4) are merged in these formulae, , , , out t out t T t T T , for dividing the prediction value and the uncertain part. The distribution function t Q of , T t belongs to the ambiguity set Q t , which is defined as Q { Q | ( || ) } t t t KL t t t
Q D Q Q and t is the farthest KL divergence (radius) with the nominal distribution t Q . Details of this reformulation are provided in Appendix A. Thus, the KL-DRO model for the HP day-ahead operational schedule is described as follows.
Model 2: KL-DRO model min (6) . . s t (9), (11)(12)(14), (16)-(21) , {0,1} k t x ,
1, , , 1,2, , , t H k N Z
Selection of the nominal distribution and the radius will be further discussed in subsections B and D . B. Reformulation
The nominal distribution P of the ambiguity set contains all of the information driven from historical data. The historical data are the deviation records of predicted versus real values. Referring to [18], a DRO expectation constrained program under KL divergence ambiguity such as min ( ). .max [ ( , )] 0 x PP h xs t E H x (22) { | ( || ) } KL P D P P can be transformed into the following form on the basis of strong duality theory: ( , )/0 min ( ). .min{ ln [ ] } 0 x H xP h xs t E e (23) if for every x X , ( , ) { | 0, [ ] } sH xP S s s E e is not empty. Based on the theorem, the expectation constraints, (16)-(21), in Model 2 can be transformed into , max0 [ ]/0 min{ ln [ ] } 0
Z Zk HP k sum Pk Z x P P PP
E e (24) [ ]/0 min{ ln [ ] } 0 k k k out k k k T J x K T l T KQ
E e (25) [ ]/0 min{ ln [ ] } 0 k k k out k k k T J x K T l T KQ
E e (26) [ ( ) ]/0 min{ ln [ ] } 0 w k k w k out w k w k T w k t set
Z M x N T p N TQ
E e (27) where ,1 ,2 , TZ Z Z Zsum sum sum sum H
P P P P and has the same definition format for , , ,
P T . If the expectation expressions exist, then the minimization problem over can be solved and these constraints are transformed to deterministic ones. Thus, the nominal distributions , t t P Q need to take an exponential form. Two ways are proposed here to obtain the nominal distribution according to the different methods of using the historical data. One is the Gaussian assumption (GA), in which it is assumed that the nominal distribution obeys a Gaussian distribution, the mean and variance of which are calculated from the first and second moments of the historical data, respectively. The other way is the KDE approximation (KDEA) [23], in which the uncertain variable has the following distribution function:
1( ) iNN iN N f HNh h (28) Here, N is the amount of valid historical predicted-actual errors data i , N h is a positive constant number, and ( ) H is a smooth function satisfying ( ) 0 H , ( ) 1 H d , ( ) 0 H d , and ( ) 0 H d [15]. We choose ( ) H as a standard normal distribution, (0,1) N . Then, the ( ) N f takes the sum of N normal distributions, which have the same variance N h but a different mean i : ( )21 iN N hN i N f eN h (29) It was shown in [24] that the KDE function converges to the real distribution in a 1-norm sense, | ( ) ( ) | 0 as K N f f d N (30) Thus, we reach the following formulations for constraints (24)–(27) and show an explanation for (24) here. (1) Gaussian Assumption (GA) Suppose
2, , , ( , )
P t P t P t N , and for t H , (24) equals
2, ,, , , max 2, , ( )[ ]/ / 2 ,0 ,
P t P tZ Zk t HP k sum t P t P tk Z x P P P P t tP t e e e d Determining the minimum of , we have , , , max , , Z Zk t HP k sum t P t P t tk Z x P P P (31), which is a tightened deterministic constraint of (1). The derivation procedure is listed in Appendix B1. (2) KDE Approximation (KDEA) Under KDEA, for t H , the constraint (24) is transformed into ,, iP tP t NN PZ Zk t HP k sum t tk Z iP t hx P P P eN (32) where , P t N is the amount of power forecast errors for time t (details are shown in Appendix B2). We show that with given specific t and , N P h and historical power forecast error data, , iP t , the function , ( ) P t g ,
1( ) ln , 02 iP t
NN PP t t i hg eN (33) is a convex function (Appendix C). Then, we denote the minimum value of , ( ) P t g as , ,min ( ) P t g and take ,min ,1,min ,2,min , ,min ( ) ( ) ( ) ( ) TP P P P H g g g g . Thus, constraint (24) is equivalent to , max ,min ( ) 0
Z Zk HP k sum Pk Z x P P P g (34) which is also a deterministic constraint. Similarly, (25)–(27) also can be transformed to deterministic forms: ,min ( ) 0 k k k out k T k k J x K T K g l T (35) ,min ( ) 0 k k k out k T k k
J x K T K g l T (36) ( ( ) ) 0 w k k w k out w k T w k w k t set
Z M x N T N g p T (37) where ,min ,1,min ,2,min , ,min ( ) ( ) ( ) ( ) TT T T T H g g g g , ,,
1( ) ln , 02 iT tT t
NN TT t t iT t hg eN (38) where , T t N represents the temperature forecast error. The KL-DRO with KDEA methods can be reformulated as Model 3, shown below. Model 3: KL-DRO-KDEA min (6) . . s t (9), (11)(12)(14), (34)-(37) , {0,1} k t x ,
1, , , 1,2, , , t H k N Z C. Historical data utilization with information loss
It is clear that the knowledge of uncertainties, concentrated from the available data, dominates the quality of the DRO program. Based on information theory, the uncertain variables and the historical data could be regarded as an information source and many independent observations. We apply the minimum information loss theorem as Bb bloss b b
PI f P B f (39) where ( ) b P and ( ) b f are the integral probability and the dropped-into frequency of bin b , respectively, which means using frequency to estimate probability by depicting a histogram of all the b . The measured space is divided equally into B bars. When B , we have: B b KLB b b
P PP d D P fB f f (40), which indicates that if we use the histogram method to estimate the real distribution of the uncertain vectors, we would, at the same time, show the least information loss using the KL divergence. The KDE function is an approximation of the histogram with an appropriate N h and converges to the histogram when , N B [23], further demonstrating the superiority of this method. D. Risk meaning of the radius of an ambiguity set
The radius of the ambiguity set dominates the possible deviated degree of the nominal distribution and is obviously linked to the conservatism of the DRO program. In [18], it was shown that a EC-DRO is equivalent to a chance-constrained program, where radius reflects a decision maker’s risk level, , with equation e . For example, when risk level , the radius, . This provides a reference for selecting a suitable radius of our DRO model. IV. N UMERICAL TESTS
In this section, we compared performance between the conventional RO, the KL-DRO with Gaussian assumption (GA-DRO), and the KL-DRO with KDEA assumption (KDEA-DRO) via Monte Carlo simulations. A. Test settings
The test system includes 10 heterogeneous houses with HPs and water tanks; their parameters are listed in Table 1. The COP of the HPs and the water-to-house efficiency are set as constants: COP and w h . The comfort range for indoor temperature is set to be [18,24] C , and S kW . A finite difference model of the house [25] is used for the thermal simulations with 5 s for each step. The period number is 288 per day with 5 min for each time interval. Table 1 . Parameters of 10 houses, tanks, and heat pumps
House 1 2 3 4 5 6 7 8 9 10
R ( °C/kW) kJ/°C) (°C)
19 20 21 20 19 19 20 21 20 19 R w (°C/kW) w (kJ/°C) w0 (°C)
42 45 47 45 42 42 45 47 45 42 P HP (kW)
5 4.7 4.3 4.7 4.8 5 4.7 4.3 4.7 4.8 -2 -1 0 1 2024 Prediction errors F r equen cy F r equen cy
12 h-2 -1 0 1 200.511.5 Prediction errors F r equen cy
18 h -2 -1 0 1 200.51 Prediction errors F r equen cy
24 h
Figure 2. Histogram of 92 days’ outdoor temperature prediction errors for four typical hours from day-ahead 23:30 to the coming day.
As shown in Figure 2, the day-ahead predicted errors of outdoor temperature were recorded for 92 days from November 1, 2016, to January 31, 2017, in Beijing, based on data in [26]. It can be seen that the early time nodes show relatively less forecast deviation. The prediction errors of the zonal power, the summation of PV, and load are generated randomly with an asymmetric χ distribution, P rand for 5000 sample data for every moment. We tested the performance of different , N P h of the KDE function with power predicted errors as shown in 6 Figure 3. This shows that a larger , N P h may lead to a closer KDE function with more burr. We finally chose , N P h and , N T h for an appropriately smooth KDE function. -4 -2 0 2 400.511.5 Prediction errors F r equen cy hn = 0.01 -4 -2 0 2 400.511.5 Prediction errors F r equen cy hn = 0.1-4 -2 0 2 400.511.5 Prediction errors F r equen cy hn = 1 -4 -2 0 2 400.511.5 Prediction errors F r equen cy hn = 10 Figure 3. Histogram and KDE function of 5000 zonal power prediction errors
The predictions of PV power and the outdoor temperature curve are from real measurements on February 1, 2017, Tsinghua University, Beijing (Fig. 4). The TOU electricity price is shown in Figure 5; the mean value is kWh . The penalty weight for transformer capacity is
10 $ / kW for the test cases. P o w e r f o r e c a s t ( k W ) LoadPhotovoltaic
Figure 4. Day-ahead PV and load prediction T i m e o f u s e e l e c t r i c i t y p r i c e ( $ ) Figure 5. Curve of TOU electricity price
For the radius distance of the ambiguity set, we let t t t for a 0.9 confidence level. Numerical experiments were conducted in MATLAB 2014b with a core-i7 laptop. The minimization of ( ) g was calculated by an inner-point method and the run time increased slightly when the historical data set increased. The mixed integer program was solved using CPLEX 12.5 with a 1% tolerance, and the average solving time for a KDEA-DRO problem was about 10 s. B. Comparison of KDEA-DRO with unscheduled conditions
For unscheduled HP operational conditions, the ON/OFF state of the HP obeys a simple hysteresis rule: ,, , , ,, 1 , , ,,
1, 0 & 0, 1 & , w kk t w k t w kk t k t w k tk t if x T Tx if x T Tx other conditions (41) The water temperature bounds ,, [ , ] w kw k T T were set at [40,45] °C here. The results for KDEA-DRO and unscheduled operation are shown in Figure 6 and Table 2. For the unscheduled scenario, Figure 6 shows that the transformer is overloaded from 17:30 to 19:30, when the electricity price is high. In contrast, the load profile of the transformer is flat and kept within its security limits with KDEA-DRO. Beyond that, HPs in the scheduled model turn off at high price moments to save money, and turn on at night for thermal energy storage, indicating good performance in load shifting. The large numerical differences of these models are shown in Table 3. The scheduled model cuts around 34% of the maximum daily transformer power and reduces the electricity cost for residents by about 18%. P o w e r ( k W ) P o w e r ( k W ) TransformerHeat pumps
Figure 6. Operational power of transformer and heat pumps of the unscheduled model and KDEA-DRO model with y = 10 kW/$. Table 2 . Main results of unscheduled operational model and KDEA-DRO model P max ( kW) Peak cost ($) Elec. cost ($) Total ($) Unscheduled 73.204 732.0 962.7 1694.7 KDEA-DRO 48.612 486.1 787.6 1273.7 C. Comparison of deterministic, KDEA-DRO, GA-DRO, and RO models
Operational strategies were generated using the deterministic, KDEA-DRO, GA-DRO, and RO models with the same setting and confirmed with 1000 Monte Carlo experiments. The intervals for RO were chosen with a 95% falling rate of 7 historical data and are symmetrical about zero (the concrete formulation of RO is shown in Appendix D). The forecast errors of outdoor temperature were generated by normally distributed random numbers, the means and variances of which were calculated from historical data, and the errors in power prediction were generated by the manipulated χ distribution mentioned above. The comfort rate (CR) is defined as the time between the comfort bounds divided by the total time. , ( )min ( ) k k t kday k Num T T TCR Num t (42)
Table 3 . Results of Monte Carlo experiments for the operational strategies generated with different models
Method Deterministic KDEA-DRO GA-DRO RO P max /kW Best day 44.38 45.82 48.09 48.76 Worst day
Mean
Elec cost /$ Best day 734.4 702.8 720.5 725.7 Worst day
Mean
Comfort rate Best day 100.0% 100.0% 100.0% 100.0% Worst day
Mean T e m pe r a t u r e ( ℃ ) Deterministic 0 6 12 18 2402040 KDEA-DRO0 6 12 18 2402040 Time (h) T e m pe r a t u r e ( ℃ ) GA-DRO 0 6 12 18 2402040 Time (h)ROComfort bound indoor airwateroutdoor air
Figure 7. Temperature simulations of deterministic, KDEA-DRO, GA-DRO, and RO model for 10 houses
From Table 3, we can see that the DRO methods have relatively high comfort satisfaction with less energy cost than the RO method. Because the uncertainties are not considered, the deterministic method has a poorer comfort rate and a higher cost. Among the two DRO methods, max P and the energy cost with KDEA-DRO are less than with GA-DRO, because the KDE function is more accurate than the Gaussian function in representing historical information, while the comfort satisfaction of KDEA-DRO is acceptable. Regarding the indoor comfort aspect, outage of indoor temperature of one house can be observed in the deterministic curves in Figure 7, while the robustness is better in the three robust methods. D. Discussion: Influence of the divergences t and t The ambiguity sets of the probability distribution of the uncertainties in our KDEA-DRO model are controlled by the divergences t and t . As discussed in Section III, for an ambiguity set, the risk level and the set’s divergence , it holds that e . When the risk level increases to one, the divergence decreases to zero, and the DRO problem reduces to a stochastic optimization problem. We first analyze the ,min ( ) P g and four typical hours’ , ,min ( ) T t g , which represent a robust part in the KDEA-DRO model with different divergence settings. As shown in Figure 8, with the increase in divergence and the decrease in risk level, the ,min ( ) P g and , ,min ( ) T t g increase with a gradually decreasing rate due to the robustness of the relative constraints. Due to the higher forecast accuracy, the ,6 ,min ( ) T h g is smaller than at the three other times (Fig. 8). g m i n ( α ) g P,min (α)g
T,6h,min (α)g
T,12h,min (α)g
T,18h,min (α)g
T,24h,min (α) β=0.05 β=0.01 β=0.001β=0.1
Figure 8. g min ( ) of power and four typical hours prediction error data with different divergence η of the ambiguity set Next, we discuss the influence of divergences t and t on the results of KDEA-DRO via 1000 time Monte Carlo simulations; the other parameters are the same as those in subsection C . We select four risk levels. β = 0.1, 0.05, 0.01, and 0.001, for the power and temperature ambiguity sets for comparison; the outcomes are shown in Table 4. With the increase in power divergence, t , the power peak increases, the electricity cost decreases, and the comfort rate rises. Thus, a large power divergence could provide an advantage in reducing the electricity cost and improving the users’ thermal satisfaction. Moreover, with an increase in outdoor temperature divergence, t , the power peak increases slightly, the electricity cost increases, and the comfort rate rises. Based on these results, we can control the robustness and conservatism of the KDEA-DRO by varying the divergence settings. Table 4 . Results of the KDEA-DRO under different divergences of t and t for different risk levels β (a) P max (kW) t (risk level) t (risk level) (β = 0.05) 4.6052 (β = 0.01) (b) Electricity cost ($) t (risk level) t (risk level) (c) Comfort rate t (risk level) t (risk level) V. C ONCLUSIONS AND FUTURE WORK
In this paper, we developed a KL distance-based DRO model of HP-DOSM with residential temperature constraints , which can both decrease the peak-valley gap and the cost to residents under a TOU electricity price. This distance-based DRO model can well capture the uncertainties of weather prediction, photovoltaic prediction, and load prediction errors, while it is tractable. Numerical tests showed that our distance-based DRO was robust with less conservatism than the conventional RO model. Moreover, the robustness of this model can be adjusted by tuning the risk level, which has an explicit meaning in the optimization problem. In future work, a distributed algorithm is needed to solve HP-DOSM for large-scale distribution networks integrated with massive heat pumps. A
PPENDIX A. Detailed structural transformation of the heat pump day- ahead operational schedule model
We first apply several simplifications. Transform the (3) and (4) equalities to a simplified form: , 1 , , , , k t k k t k w k t k out t
T T T T (43) , , 1 , , , , , , , w k t w k k t w k w k t w k k t T T T x (44) With the definition of the vectors in the nomenclature, we have: , , 0, k kk w kk k k t setk outk T T TT . (45) This can be rewritten as , k k k w k k out k A T B T C T d (46) In the same way, we can state that , , , , , w k w k w k k w k k w k
E T F T G x h (47) It is easy to show that k A and , w k E are non-singular; thus, we have ( ) k k w k w k k k w k w k k k out k w k w k k A B E F T B E G x C T B E h d (48) ( ) w k w k k k w k w k k k out w k k w k w k k k
E F A B T F A C T G x h F A d (49) Due to the Schur complement, the matrices
1, , ( ) k k w k w k
A B E F and
1, , ( ) w k w k k k
E F A B are non-singular, so we have: ( ) ( ) k k k w k w k k w k w k k k out k w k w k k T A B E F B E G x C T B E h d (50) ( ) ( ) w k w k w k k k w k k k out w k k w k w k k k
T E F A B F A C T G x h F A d (51) represented by k k k k out k
T J x K T l (52) , , , , w k w k k w k out w k
T M x N T p (53) where k T and , w k T are decoupled. B. Detailed reformulation of Gaussian Assumption and KDE Approximation in constraint (24)
B.1 Gaussian Assumption
The inner formula of (24) equals
2, ,, , , max 2, , ( )[ ]/ / 2 ,,
P t P tZ Zk t HP k sum t P t P tk Z x P P P P t tP t e e e d P tZ Zk t HP k sum t P t tk Z x P P P (54) Due to
2, ,0 min{ } 22
P tt P t t , we take (24) as , , , max , , Z Zk t HP k sum t P t P t tk Z x P P P (55) B.2 KDE Approximation
The inner formula of (24) equals
2, ,,, , , max 2, , ( )[ ]/ / 2 ,1, , iP t P tZ Z P tk t HP k sum t P t N Pk Z
Nx P P P h P t tiP t N P e e e dN h ,, iP tP t NN PZ Zk t HP k sum t tk Z iP t hx P P P eN (56) C. Convexity proof of ( ) g The convexity of ( ) g can be calculated by a second order differential function:
1( ) ln , 02 i NN i hg eN (57) 9 eeln e'( ) 2 ii ii NN iN i Ni h Ng (58)
212 123 31 21 ( )''( ) ee 0ee ii ii
N i Ni iN NN ii i hg (59) Thus, ( ) g is convex. D. Robust optimization model min (6) . . s t , max[ , ] max 0
P P P
Z Zk HP k sum Pk Z x P P P (60) [ , ] max 0 T T T k k k out k k k T
J x K T l T K (61) [ , ] max 0 T T T k k k out k k k T
J x K T l T K (62) max ( ) 0 T T T w k k w k out w k w k T w k t set
Z M x N T p N T (63) (9), (11)(12)(14) , {0,1} k t x ,
1, , , 1,2, , , t H k N Z R EFERENCES [1] C. R. Inc, "Heat pump characterization study," 2010. [2] EHPA, "European Heat Pump Association,"
Outlook European Heat Pump Statistics,
Applied Thermal Engineering, vol. 51, no. 1-2, pp. 155-165, 2013. [4] G. Z. C. D.-m. Y. J.-h. J. B. F. X.-m. Z. Liang, "Study on Rural Low-Voltage Distribution Network Planning with Coal-to-Electricity Project,"
Electrical & Electronics, http://news.xinhuanet.com/local/2017-01/25/c_1120381959.htm . [6] H. O. R. Howlader, H. Matayoshi, and T. Senjyu, "Distributed generation incorporated with the thermal generation for optimum operation of a smart grid considering forecast error,"
Energy Conversion and Management, vol. 96, pp. 303-314, 2015. [7] D. Vanhoudt, D. Geysen, B. Claessens, F. Leemans, L. Jespers, and J. Van Bael, "An actively controlled residential heat pump: Potential on peak shaving and maximization of self-consumption of renewable energy,"
Renewable Energy, vol. 63, pp. 531-543, 2014. [8] D. K. Wolfram Wiesemann, and Melvyn Sim, "Distributionally Robust Convex Optimization," 2013. [9] Y. Zhang, S. Shen, and J. Mathieu, "Distributionally Robust Chance-Constrained Optimal Power Flow with Uncertain Renewables and Uncertain Reserves Provided by Loads,"
IEEE Transactions on Power Systems, pp. 1-1, 2016. [10] S. Zymler, D. Kuhn, and B. Rustem, "Distributionally robust joint chance constraints with second-order moment information,"
Mathematical Programming, vol. 137, no. 1-2, pp. 167-198, 2011. [11] P. Xiong, P. Jirutitijaroen, and C. Singh, "A Distributionally Robust Optimization Model for Unit Commitment Considering Uncertain Wind Power Generation,"
IEEE Transactions on Power Systems, vol. 32, no. 1, pp. 39-49, 2017. [12] E. D. a. Y. Ye, "Distributionally Robust Optimization under Moment Uncertainty with Application to Data-Driven Problems,"
Operation Research,
Optimization and Control,
Mathematical Programming, journal article vol. 158, no. 1, pp. 291-327, 2016. [16] X. Chen and Y. Zhang, "Uncertain Linear Programs: Extended Affinely Adjustable Robust Counterparts,"
Operations Research, vol. 57, no. 6, pp. 1469-1482, 2009. [17] D. d. H. By Aharon Ben-Tal, Anja De Waegenaere, Bertrand Melenberg, Gijs Rennen, "Robust solutions of optimization problems affected by uncertain probabilities," 2011. [18] Z. Hu, Hong, L. Jeff, "Kullback-Leibler Divergence Constrained Distributionally Robust Optimization,"
Conference Proceedings,
European Journal of Operational Research, vol. 233, no. 3, pp. 459-473, 2014. [21] W. Wu, Xin Chen, B. Zhang, C. Lin, "Data-driven DG Capacity Assessment Method for Active Distribution Networks,"
IEEE Transactions on Power Systems,
IEEE Transactions on Sustainable Energy, vol. 4, no. 1, pp. 182-191, 2013. [23] J. Pintér, "Deterministic approximations of probability inequalities," 1989. [24] L. G. L. Devroye, "Nonparametric Density Estimation: The ℓ1 View," 1985. [25] Q. Z. Qisen Yan,
Building thermal process . 1986. [26] http://data.cma.cn/site/index.html.
China Meteorological Administration ..