Energy Spatio-Temporal Pattern Prediction for Electric Vehicle Networks

1 Ener gy Spatio-T emporal P attern Prediction for Electric V ehicle Networks Qinglong W ang Abstract Information about the spatio-temporal pattern of electricity energy carried by EVs, instead of EVs themselv es, is crucial for EVs to establish more effecti ve and intelligent interactions with the smart grid. In this paper, we propose a framework for predicting the amount of the electricity energy stored by a large number of EVs aggregated within different city-scale regions, based on spatio-temporal pattern of the electricity energy . The spatial pattern is modeled via using a neural network based spatial predictor , while the temporal pattern is captured via using a linear-chain conditional random field (CRF) based temporal predictor . T wo predictors are fed with spatial and temporal features respectiv ely , which are extracted based on real trajectories data recorded in Beijing. Furthermore, we combine both predictors to b uild the spatio-temporal predictor , by using an optimal combination coefficient which minimizes the normalized mean square error (NMSE) of the predictions. The prediction performance is ev aluated based on extensiv e experiments cov ering both spatial and temporal predictions, and the improvement achieved by the combined spatio- temporal predictor . The experiment results sho w that the NMSE of the spatio-temporal predictor is maintained below 0 . 1 for all in vestig ate regions of Beijing. W e further visualize the prediction and discuss the potential benefits can be brought to smart grid scheduling and EV charging by utilizing the proposed framework. I . I N T RO D U C T I O N Electric v ehicles (EVs) [1] are increasingly regarded as not only a promising alternati ve to fossil fuel engine vehicles for en vironmental and efficienc y concerns, but also flexible distributed ener gy facilities which are capable of both consuming and pro viding electricity energy via interacting with the smart grid [2 – 4]. Recent thri ving dev elopment of battery techniques has greatly accelerated this tendency [5]. While batteries of EVs hav e been advanced in term of both capacity and power density , this in turn indicates that EVs will consume more energy with lar ger po wer . F or example, according to the U.S. Energy Information Administration [6], in 2013, a U.S. residential utility customer consumes an a verage of 909 kWh of electricity per month, while an EV with a 24 kWh battery pack consumes around 720 kWh if charged once daily . Meanwhile, the char ging power can be as high as 19 . 2 kW when charged by le vel 2 charging standard February 15, 2018 DRAFT 2 ( 240 V and 80 A) [7]. Obviously , electricity load from EV charging would pose significant impacts on the smart grid if not handled properly [8, 9]. Meanwhile, EVs also bring new opportunities to the smart grid. Utilizing EVs to provide regulation services and de veloping V ehicle-to-Grid (V2G) technique have been more and more intensiv ely discussed [2, 3, 5]. Indeed, it is actually the electricity energy carried by EVs that can influence and also benefit the smart grid. As EVs moving across streets, districts and e ven cities, the electricity ener gy is carried along. Hence, the spatial pattern of EVs’ mobility is inherited by the electricity energy stored in battery packs. Meanwhile, the stored energy keeps decreasing as an EV keeps driv en by electric po wer . The changes of the stored energy also possess temporal pattern. Therefore, to mitigate the impact and also improv e the efficiency of using EVs’ batteries for pro viding regulation and V2G services, it is crucial to discover the spatio-temporal pattern of the electricity energy carried by EVs. Spatially related energy patterns can be explored by observing the mobility pattern of EVs. There have been many research works aiming at studying traffic pattern [10 – 12] and driving pattern of indi vidual vehicle [13 – 15]. Howe ver , spatial pattern of the energy carried by EVs are rarely explored but is critical for lar ge-scale energy planning. Meanwhile, temporal changes of energy are due to the usage of electricity energy for propulsion. Research w orks studying the char ging and discharging scheduling of EVs either focus on modeling energy changes of an indi vidual EV [15, 16], or of aggregated EVs restricted in static scenarios, where EVs are plugged-in and the scheduling is only for the plugged-in period [17 – 19]. [20] has considered the influence of EVs’ mobility on the stored ener gy and proposed a coordinated charging strategy . Our work is different in that we view EVs’ energy pattern through both spatial and temporal lens. W e aim at building an energy spatio-temporal pattern prediction framew ork, upon which crucial questions including when, where and how much is the av ailable energy for effecti ve scheduling can be answered. T o realize this framew ork, v arious types of sensing and monitoring facilities are needed. Although there hav e been man y e xisting techniques a vailable to capture and transmit EVs’ mobility information, e.g. GPS, vehicular ad-hoc networks (V ANETs) [14, 21, 22], dedicated short range communication (DSRC) [23], cellular networks [24, 25] and W i-Fi [26, 27], etc. Howe ver , when the focus is on energy pattern, these techniques are not suf ficient. W e need additional functionality of providing information of energy storage status, e.g. state-of-charge (SOC). This can be realized through transmitting SOC information from an EV to road-side communicating facilities. e.g. road-side-units (RSUs) [20, 28], and further submitting SOC information from RSUs to central control facilities. Our work shows that, with communicating techniques February 15, 2018 DRAFT 3 and facilities that are capable of collecting SOC information, EVs’ energy pattern can be discov ered, thus contributes to establishing more intelligent interaction with the smart grid [29]. In summary , the contributions of this paper are as follows: 1) W e propose a spatio-temporal prediction framew ork to predict the aggregated energy carried by a large number of EVs within city-scale regions. Both spatial and temporal patterns are utilized for building a combined spatio-temporal predictor . 2) W e extract both spatial and temporal features for exploring the hidden correlation between them and the changing pattern of the aggregated energy . 3) W e utilize a real city-scale trajectories dataset for feature extraction and prediction performance e valuation. Lo w prediction errors demonstrate the ef fecti veness of the proposed spatio-temporal frame work. The rest of this paper is or ganized as follo ws: Section II introduce the data set used and spatial and temporal features extraction procedures. In Section III we provide an ov erview of the energy spatio-temporal pattern prediction framew ork. Section IV presents the design of the spatial predictor and the temporal predictor in detail. W e ev aluate the prediction performance of the proposed predictor in Section. V , including the demonstration of the influence on the prediction performance by selecting different features, and the comparison of the prediction performance of the spatial predictor , the temporal predictor and the spatio- temporal predictor . Finally , Section VI concludes the paper and address our future work. I I . D A TA S E T A N D F E A T U R E E X T R A C T I O N The dataset used in this paper is based on a GPS trajectories dataset of 10 , 357 taxis from Feb . 3 to Feb . 8, 2008 in Beijing [30]. As most of taxis recorded tra vel within the 4 th Ring Road of Beijing, we tar get on this area and further divide it into sixteen disjointed regions, as shown in Fig. 1. W ithin each re gion, we assume there e xist at least one road-side facility responsible for collecting the SOC information reported by each EV , and transmitting it to the central controller . Assuming all the trajectories data is generated by EVs, we further introduce the extraction of both spatial and temporal features based on the trajectories data. In term of spatial features, for each di vided region, information related to the quantity , velocities and moving directions of EVs within such a region is considered. In term of temporal features, the hourly aggregated energy of all EVs within each region and the associated recording time is included. T able. I sho ws the extracted spatial and temporal features for the k th region, k = 1 , · · · , 16 . February 15, 2018 DRAFT 4 4 3 2 1 8 7 6 5 12 11 10 9 16 15 14 13 Fig. 1. Spatial distribution of all recorded taxis over sixteen divided regions within 4th Ring Road of Beijing. T ABLE I A L L E X T R AC T E D F E A T U R E S F O R R E G I O N k = 1 , · · · , 16 Features Components Meaning Spatial features F k N Number of EVs F k V A verage velocity: V k ave V elocity variance: V k var F k D Number of EVs in direction 1: D k 1 Number of EVs in direction 2: D k 2 Number of EVs in direction 3: D k 3 Number of EVs in direction 4: D k 4 V ariance of the number of EVs in four directions : D k var F k E Aggregated energy: E k sum Aggregated energy variance: E k var T emporal features F k E Aggregated energy: E k sum F k H Recording hour: H A. Spatial F eature Extraction There ha ve been many models proposed for studying traffic patterns. For example, [31, 32] draw an analogy between circuit and traffic networks, and use Kirchhof f ’ s law to constrain traf fic network models. Building models of this type can be ef fectiv e for a single road se gment or an intersect, as the traf fic flow is constrained by explicit and simple conditions. Howe ver , if vie wed at a much larger scale (for e xample, the February 15, 2018 DRAFT 5 Fig. 2. Demonstration of energy temporal patterns of a selected region. From left to right are aggregated energy , recorded number of EVs, normalized aggregated energy and discrete normalized aggregated energy of 24 hours of five days. Especially , in second graph on the left, the recorded numbers of EVs of all observed days share a similar varying pattern. This implies that most taxis maintain similar daily driving pattern, which can be applicable for priv ate EVs. whole city as studied in this paper), models of this type [33] can bring much complexity when considering the countless roads, lanes and intersects, etc. Moreover , as we aim to explore the hidden correlation between the aggregated energy , instead of the traffic, among adjacent areas, simply employing any specific models may lead to sev ere bias problems. Up to the present, there is no clear form for modeling the flow patterns of energy when having moving vehicles as the carriers. Howe ver , luckily there hav e been many research w orks emplo ying neural networks to handle the complex modeling of traf fic networks [34, 35]. Therefore, as sho wn in T able. I, we extract dynamic spatial features representing general aspects of a selected area (compared with specific ones representing a single road or intersect). Meanwhile, v arious static features including terrain, real estates along road segments, road number and capacity have been proposed in [36]. Note that although these static features also represent a specific re gion, the y remain unchanged and can hardly represent the changing pattern of energy . As sho wn in T able. I, for each region k = 1 , · · · , 16 , spatial features contain follo wing four feature sets: • Firstly , we use F k N ( t ) to denote the number of EVs within region k during t th recording hour . F k N ( t ) is calculated based on the GPS location and timing records contained in the original trajectories dataset. • Secondly , we calculate the driving distance S n m − 1 ,m between two adjacent GPS locations P n m − 1 and P n m for EV n . W e also calculate the time dif ference ∆ t n m − 1 ,m between corresponding time records. Hence, the velocity of EV n at recording time t can be obtained as V n m − 1 ,m ( t + 1) = S n m − 1 ,m / ∆ t n m − 1 ,m . Furthermore, we get additional features including the a verage and variance of the velocity as V k av e ( t ) and V k v ar ( t ) for region k . February 15, 2018 DRAFT 6 • Thirdly , we utilize GPS location records to calculate the moving direction d n m of EV n , as d n m =  P n m , P n m − 1  T . d n m is further determined to be within one of four quadrants D 1 , D 2 , D 3 and D 4 . W ith d m of all EVs, we calculate the v ariance of the moving directions as D k v ar for region k . • Fourthly , we obtain the energy related feature set F k E and recording time F k H as introduced in Section. II-B. After calculating the spatial features of each re gion, these features are fed to a neural network based spatial predictor to capture the hidden correlation of the aggregated ener gy among adjacent regions. The input for the spatial predictor consists of the four spatial feature sets, F N , F V , F D and F E of all neighbours of region k (For e xample, for region 6 , its neighbours include region 1 , 2 , 3 , 5 , 7 , 9 , 10 and 11 . While for regions residing at the margins of Fig. 1, take region 1 as an example, we pick regions 2 , 3 , 5 , 6 , 7 , 9 , 10 and 11 as its neighbours). B. T emporal F eatur e Extraction T emporal features consists of recording hour H h , where h = 1 , ..., 24 and energy related feature set F E . For the former , it can be easily obtained by checking the GPS timing records. The latter is extracted based on previously calculated distance S n m − 1 ,m of EV n . W e first assume all the EVs have the same battery capacity C = 24 k W h . Then we extract the SOC of each EV by setting it linearly dependent on the distance this EV is driv en, i.e. the more distance an EV travels, the lower is its SOC. This simplified model has also been used in [20, 37]. Since the trajectories data is generated by taxis, we notice that most taxis recorded keep moving e ven at midnight. This implies that there is no uni versal shift time for these taxis. Therefore, we need to assign an uni versal recharging time for all these taxis, giv en we hav e assumed them to be EVs. W e assume they are all recharged to full capacity (i.e. SOC is initiated to 1 ) at 0 a.m. of e very observed day from Feb . 3 to Feb . 8, 2008. This assumption is practical when trajectories data of common pri vate EVs is av ailable, since people usually have their EVs recharged to full capacity for ne xt day’ s driv e. In this case, SOC of pri vate EVs are most likely to be initiated to 1 before an univ ersal work time. Similar to spatial feature extraction in Section. II-A, we calculate the v ariance of remaining electricity E k v ar ( t ) (obtained via using all EVs’ reported SOC and capacity C ) and the aggreg ated remaining energy E k sum ( t ) of all EV within k th region at time t . The hourly updated E k sum ( t ) of six observed days of an example region is sho wn in the first graph on the left in Fig. 2. From Fig. 2, we notice that all observed days of this region share a similar varying pattern of aggregated energy , although with dif ferent scales. This February 15, 2018 DRAFT 7 can be explained by observing the fluctuation of the number of EVs in this region (shown in the second graph on the left in Fig. 2). F or all observed days, EVs crowd into and depart from this region at similar time. T ake 12 p.m. of e very day as an example, when the aggregated energy starts to drop, the number of EVs remains at the maximal value for hours. This is because when the number of EVs first reaches the highest value, EVs have more remaining ener gy , compared with hours later when the number of EVs starts to drop. The peak of aggregated energy lasts shorter than that of the number of aggregated EVs. As it is the v arying pattern of aggregated energy that we aim to capture, it is necessary to mitigate the impact of the scale differences of the aggre gated ener gy among each observed days. W e normalize the aggregated energy E k sum ( t ) at time t by the summation P 24 t =1 E k sum ( t ) of the day it is recorded. The normalized result is sho wn in the third graph on the left in Fig. 2. Clearly , the varying pattern of aggregated energy is well preserved. Furthermore, in general, continuous variable CRF [38, 39] often leads to exponential computation complexity [40, 41]. In this paper , we adopt discrete variable CRF . W e discretize the normalized aggregated energy into 10 le vels, which is sufficient to represent the varying patterns. The discretized aggregated energy is shown in the fourth graph on the left in Fig. 2. I I I . F R A M E W O R K O V E RV I E W The ener gy spatio-temporal prediction framew ork consists of three phases: (1) offline training of predictors, (2) optimal combining of both predictors and (3) online predicting, as sho wn in Fig. 3. W e split both temporal and spatial features into a training set, a v alidation set and a testing set respecti vely . During of fline training, the spatial and temporal predictors are trained using the training set independently . In the second phase, two predictors are combined via using a combination coef ficient λ . The optimal λ ∗ is obtained by minimizing the normalized mean square error (NMSE) of the v alidation set. Finally , both trained predictors and the optimal combination coef ficient λ ∗ are tested by the testing set. A. Offline T raining During offline training, we first select temporal and spatial features from Feb . 3 to Feb . 6, 2008 as the training set for the temporal predictor and the spatial predictor respecti vely . 1) T raining T emporal Predictor: The linear-chain CRF based temporal predictor is fed with temporal features (i.e. recording hour F k H and energy related features F k E ) of each region k = 1 , ..., 16 , and then trained independently . February 15, 2018 DRAFT 8 Spatial predictor Temporal predictor Offline training Combination coefficient optimization Data from Feb. 2 to Feb 6, 2008 Data on Feb. 7, 2008 Optimal combining Combination coefficient Data on Feb. 8, 2008 Spatio-Temporal predictor Online predicting Offline training Optimal combining Online predicting Fig. 3. Overview of energy spatio-temporal pattern prediction framew ork. 2) T raining Spatial Pr edictor: Spatial features (i.e. F k N , F k V and F k D of all neighbour regions of region k ) are gathered independently for each region. The training for the spatial predictor is hourly based, meaning that the spa tial predictor takes all spatial features at current time t as input, while taking the actual aggreg ated energy v alue at time t + 1 as target. B. Optimal Combining W e take temporal and spatial features on Feb . 7, 2008 as the validation set for obtaining the optimal combination coef ficient λ ∗ . Both temporal and spatial features are fed into the trained temporal and spatial predictor respectiv ely . Since the SOC of all EVs are initiated to 1 at the be ginning hour of each day (see Section. II for detail), we need to assign an initial v alue for the temporal predictor to start. W e av erage the aggregated ener gy values of the starting hour of the entire training set as the initial value. The predictions of both temporal predictor and spatial predictor are conducted hourly and independently . The temporal prediction results are the hourly probability distributions of different lev els of aggre gated energy for each region, while the spatial prediction results are the aggregated energy v alues of all regions. After obtaining the prediction results of all 24 hours, the temporal and spatial prediction results for all regions are combined via λ , for which the optimal λ ∗ minimizes the NMSE of predictions for the validation set. February 15, 2018 DRAFT 9 C. Online Pr edicting Online predicting is similar to the procedure of optimizing combination coef ficient, except that the optimal λ ∗ has been set. Both spatial predictor and temporal predictor are retrained with features extracted from Feb . 2 to Feb . 7, 2008, and tested by features extracted on Feb . 8, 2008. Additionally , during online predicting, spatial predictors are retained repeatedly e very hour when ne w observations are a vailable. For temporal predictor , it is retrained when the observations of an entire ne w day are obtained. I V . S PA T I A L & T E M P O R A L L E A R N I N G In this section, we introduce the design of the spatial predictor and the temporal predictor in detail. Then the combining method for building the spatio-temporal predictor is presented. The optimal combining coef ficient is obtained through minimizing the NMSE of the prediction results generated by the spatio- temporal predictor . A. Design of the Spatial Pr edictor W e utilize neural netw ork to capture the complex spatially related energy pattern, which inherits the changing pattern of mo ving traffic. For the latter , neural network has been widely employed for performing traf fic predictions [34, 35]. W e utilize a two-layer feedforward neural network, of which the neurons in the hidden layer are associated with a sigmoid transfer function, and a linear function for the output layer . For each region k , its spatial features are I k =  F k N , F k V , F k D , F k E  T , of which each component denotes a sequence of hourly recorded spatial features of region k . For example, F k N =  F k N (1) , · · · , F k N ( T train )  T , where T train is the total number of time slots used for training. The input to the neural netw ork for re gion k is I k =  I k 1 , · · · , I k N  , which is the set of spatial features of all its neighbor regions k 1 , · · · , k N ∈ R k , where R k denotes the neighbour set of region k . Meanwhile, the target is O k = E k sum , where E k sum is a sequence of aggreg ated energy within the k th region, i.e. E k sum =  E k sum (1 + ∆ T ) , · · · , E k sum ( T train + ∆ T )  T . ∆ T is the prediction time interval, which can be set from 1 to 24 . The testing of the spatial predictor is conducted hourly via using ne wly measured features at time t as input to predict the aggregated ener gy value in the following time slot t + 1 . When the measured features of time t + 1 are av ailable, these newly obtained features are appended to the training set, with which the spatial predictor is retrained. February 15, 2018 DRAFT 10 B. Design of the T emporal Pr edictor W e employ a linear-chain CRF to b uild our temporal predictor , which utilizes temporal features of each region to predict the future aggreg ated energy within that region. A linear-chain CRF has the advantages of relaxing independence assumption o ver hidden Mark ov models and a voiding the fundamental label bias problem over maximum entropy Markov models and other discriminativ e Markov models based on directed graphical models [40, 42]. A linear -chain CRF is for b uilding undirected probabilistic models for segmenting sequence data, as the graph G = ( V , E ) shown in Fig. 4. X is the observation set, consisting of a sequence of observations X t , i.e. X = { X 1 , X 2 , · · · X T } . In our case, we set each observation X t as the recording hour F k H ( t ) of k th region at time t . Random state v ariable Y is index ed by the vertices of G , and consists of components Y t , which varies ov er the label set comprising 10 discrete levels of normalized aggregated energy . Giv en X and Y , ( X , Y ) is a CRF with a graph structure as sho wn in Fig. 4. Therefore, according to the Markov property , the conditional distribution of random state variables Y t follo ws: p ( Y t | X , Y τ , τ 6 = t ) = p ( Y t | X , Y τ , τ ∼ t ) , (1) where τ ∼ t means that τ and t are neighbors in graph G . Note that the linear -chain CRF is globally conditioned on observation set X . Denote an observation sequence as x and a label sequence as y respecti vely , and ha ve y | e and y | v represent the set of components of y associated with the edges and v ertices of graph G = ( V , E ) . The conditional distribution of y gi ven x is as follows: p θ ( y | x ) ∝ exp  X e ∈ E , t γ t f t ( e, y | e , x ) + X v ∈ V , t µ t g t ( v , y | v , x )  , (2) where f t and g t are the edge and vertex potential functions which represent the features associated with edges and vertices respectiv ely . θ = ( γ 1 , γ 2 , · · · ; µ 1 , µ 2 , · · · ) are numerical weights assigned to all features. 1 t X  t X 1 + t X 1 t Y  t Y 1 t Y  Fig. 4. Graphic structure of a linear-chain CRF . February 15, 2018 DRAFT 11 In our case, we set both f t and g t as Boolean variables to indicates recording hour F k H ( t ) . Giv en training data D = { x n , y n } N n =1 , where each x n is a sequence of observ ation, and each y t is a sequence of state v ariables, parameters θ can be estimated via maximizing the following conditional log-likelihood function: L ( θ ) = N X t = n log p ( y n | x n ) . (3) T o combat o ver -fitting, we further use L 2 regularization, which penalizes v ectors with too large norms. The maximization of (3) over θ is solved via gradient descent. Additionally , the computation of gradient requires marginal distribution for each edge transition, which is solved by forward-backward recursions. C. Optimal Combining Method After training the temporal predictor and spatial predictor , these tw o predictors are fed with a validation set that does not overlap with the training set already used. For the temporal predictor , its predicted result is the probability distribution of the normalized aggregated energy over 10 discrete lev els. For example, we use p k 1 , · · · , p k 10 to denote the probabilities for le vel 1 , · · · , 10 of region k respectiv ely . Meanwhile, each class 1 , · · · , 10 is associated with Y k 1 , · · · , Y k 10 respecti vely , which is reco vered via a veraging pre vious observ ations (see Section. III for detail). For the spatial predictor , the predicted result of time t is ˆ E k S P ( t ) , which is further associated with a combination coef ficient λ . In this way , we combine the prediction results from both predictors at time t as follows: 10 X i =1 p k i ( t ) P 10 m =1 p k m ( t ) + λ Y k i ( t ) + λ P 10 m =1 p k m ( t ) + λ ˆ E k S P ( t ) , = 10 X i =1 p k i ( t ) 1 + λ Y k i ( t ) + λ 1 + λ ˆ E k S P ( t ) , = 1 1 + λ Y k T P ( t ) + λ 1 + λ ˆ E k S P ( t ) , (4) where Y k T P ( t ) = P 10 i =1 p k i ( t ) Y k i ( t ) . While (4) is only for region k , we further combine the prediction results of all regions. An uni versal optimal λ ∗ is obtained through minimizing the NMSE of predictions, as shown in (5): min λ K X k =1    E k − ( 1 1+ λ Y k T P + λ 1+ λ ˆ E k S P )    2 2 k E k k 2 2 , (5) where E k =  E k (1) , · · · , E k ( T )  T , Y k T P =  Y k T P (1) , · · · , Y k T P ( T )  T , ˆ E k S P = h ˆ E k S P (1) , · · · , ˆ E k S P ( T ) i T . February 15, 2018 DRAFT 12 V . E X P E R I M E N T S In this section, we first conduct experiments to sho w the influence on the prediction performance of using different features. Then we fix the feature sets to which generate the best prediction performance, and demonstrate the performance of the temporal and spatial predictor respectiv ely , as well as the performance gain obtained by combining both predictors with the optimal combination coefficient. W e further sho w the prediction performance changes when the prediction interv al increases to more than one hour , to explore the ability of the spatio-temporal predictor for performing longer term prediction. At last, we provide a direct visualization of the prediction. A. Pr ediction P erformance of Using Differ ent F eatur e Sets W e ha ve extracted four cate gories of features (i.e. F V , F D , F N and F E ) as input for the spatial predictor , as shown in T able. I in Section II. Among these, F V , F D and F N are more closely related to EVs’ mobility pattern. T o elaborate their influence on the prediction performance, we first di vide the whole feature sets into se veral combinations, and sho w the prediction performance of the spatial predictor and the spatio-temporal predictor based on these combinations respectiv ely (the prediction performance of the temporal predictor is not sho wn since it utilizes only temporal features, which do not change in this experiment). Due to space limit, we use abbre viations SP , TP and STP for spatial predictor , temporal predictor and spatio-temporal predictor in following tables and figures. In T able. II, different feature combinations are permutated in the descending order according to the prediction performance (NMSE) achie ved based on them. As shown in T able. II, not all features contrib ute to improving prediction performance equally . Simply using all features does not lead to the best performance (but the second worst performance in this e xperiment). The worst performance is obtained by using all feature sets except F D , while the best performance is obtained when the whole feature sets are used except F V . This result indicates that mo ving direction related features are more important for achie ving accurate prediction compared with velocity related features. This may because F D features can represent the interaction between each re gion and their neighbours, as they show directly the number of EVs mo ving between adjacent regions. Therefore, F D features are indirectly related to the energy flo w among these regions. As for F V , since the SOC of an EV is assumed to be linearly dependent on its trav eling distance, velocity features F V hav e little correlation with energy varying pattern. Furthermore, it is also noticed that using either F N or F E solely cannot lead to the best prediction performance. This indicates that traffic February 15, 2018 DRAFT 13 related spatial features indeed help improve the prediction performance. In following experiments, we select the feature combination of F D , F N , F E due to their lowest prediction NMSE. B. Hourly Pr ediction P erformance In this e xperiment, we demonstrate the prediction performance of the temporal, spatial and spatio-temporal predictor for 24 hours of a testing day in Fig. 5. W e select four regions as examples due to space constraint. Fig. 5 shows that the combination of spatial predictor and temporal predictor neutralizes the prediction errors from both. Although the spatio-temporal predictor does not generate better prediction results in all in vestigated hours, it does achie ve a better prediction performance in general. Additionally , it is also noticeable that the superiority of the spatio-temporal predictor remains consistent in all regions. Since the same combination coefficient λ is utilized by all regions, this implies that an univ ersal λ is applicable for all regions. C. A verag e Pr ediction P erformance for All Re gions In this experiment, we compare the av erage prediction performance of the temporal, spatial and spatio- temporal predictor over 24 hours of all sixteen regions. As shown in Fig. 6, the temporal predictor maintains an a verage prediction NMSE below 0 . 1 in all sixteen regions. This indicates that temporal predictor is robust to spatial changes, as it only focuses on temporal changes within each region. Meanwhile, the 5 10 15 20 0 5000 10000 15000 24 Hours Predicted energy (kWh) Real TP SP STP 5 10 15 20 0 1 2 3 x 10 4 24 Hours Predicted energy (kWh) Real TP SP STP 5 10 15 20 0 5000 10000 15000 24 Hours Predicted energy (kWh) Real TP SP STP 5 10 15 20 0 0.5 1 1.5 2 x 10 4 24 Hours Predicted energy (kWh) Real TP SP STP Fig. 5. Prediction performance of the temporal, spatial and spatio-temporal predictor , compared with the real aggregated energy value of four example regions on a testing day . The prediction results of the temporal predictor and the spatial predictor deviate from the real value in different directions. In most inv estigated hours, the spatial predictor is more intended to overestimate, while the temporal predictor is more intended to do the opposite. This performance difference between both predictors is further neutralized by the spatio-temporal predictor . February 15, 2018 DRAFT 14 R1 R2 R3 R4 R5 R6 R7 R8 0 0.05 0.1 0.15 0.2 Region 1−8 NMSE TP SP STP R9 R10 R11 R12 R13 R14 R15 R16 0 0.05 0.1 0.15 0.2 Region 9−16 NMSE TP SP STP Fig. 6. The average prediction performance over 24 hours for all regions. This results shows that the spatial predictor gets worse prediction results in most regions. Howe ver , it is also noticeable that the performance of the spatio-temporal predictor is improved in most regions. This is due to the cases when the spatial predictor generates more accurate results compared with the temporal predictor . spatial predictor shows a general inferiority compared with the temporal predictor . This causes the spatio- temporal predictor to hav e larger prediction NMSE compared with the temporal predictor for some regions. Ho wev er , as shown in Fig. 5, this inferiority of the spatial predictor does not hold constantly . Since the spatio-temporal predictor is a weighted combination of both temporal predictor and spatial predictor , it inherits the robustness property from the temporal predictor , and has both temporal predictor and spatial predictor compensate each other’ s prediction errors. Therefore, the spatio-temporal predictor also maintains an av erage prediction NMSE below 0 . 1 in all sixteen regions. W e further show the performance gain of the spatio-temporal predictor over temporal predictor and spatial predictor respecti vely in Fig. 7. R2 R4 R6 R8 R10 R12 R14 R16 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 Region 1−16 Performance gain over TP and SP TP SP Fig. 7. The performance gain of the spatio-temporal predictor . In most regions, the spatio-temporal predictor shows an improv ement of performance, especially when compared with the spatial predictor . For regions where the spatio-temporal predictor generates worse prediction results compared with the temporal predictor , larger combination weights should be assigned to temporal predictor . February 15, 2018 DRAFT 15 STP TP SP 1 hour ahead 2 hours ahead 3 hours ahead 0 0.1 0.2 0.3 0.4 NMSE Fig. 8. Comparison of the NMSE of spatial, temporal and spatio-temporal prediction over three different intervals. The largest NMSE of the spatial predictor ( 0 . 25 ) and the spatio-temporal predictor ( 0 . 09 ) are generated by 3 hours ahead predictions. The temporal predictor has the largest NMSE (0.06) for 1 hour ahead prediction. D. Long T erm Pr ediction P erformance Although hourly prediction has been employed by the ener gy market decades ago, energy prediction which can provide a forecast with longer interval also has traits when the scheduling requires long term information. In this experiment, we e xplore the ability of the spatio-temporal predictor to perform two longer term (i.e. two-hour ahead and three-hour ahead) predictions. As sho wn in Fig. 8, as the prediction interval increases, the prediction performance of three predictors changes dif ferently . Both the spatial predictor and the spatio-temporal predictor show performance degradation (the NMSE results obtained by both predictors increase), especially for the spatial predictor . In term of spatial features, the correlation between current observ ation and future results fades as the time interv al increases. Since the spatial predictor only utilizes spatial features, its prediction performance is influenced more. This influence is further inherited by the spatio-temporal predictor . Meanwhile, the prediction performance of the temporal predictor even shows a slight improvement. This indicates that the linear-chain CRF based temporal predictor is capable of performing longer term prediction. This property pre vents the performance of spatio-temporal predictor from degrading more sev erely . February 15, 2018 DRAFT 16 Fig. 9. V isualization of predicted aggregated energy of each region within Beijing 4 th Ring Road. E. V isualization of Pr edicted Ener gy W e provide a realtime visualization of the aggregated energy for practical use. A realtime visualization can assist various aspects of the interactions between the smart gird and aggregated EVs. W e utilize the spatio-temporal predictor to predict the aggreg ated energy of all di vided regions within Beijing 4 th Ring Road at 12 p.m. Feb . 8, 2008, as shown in Fig. 9. The color mosaic cov ering the Beijing map indicates dif ferent amount of aggregated energy of each region. According to Fig. 9, the regions within 3 rd Ring Road of Beijing ha ve more aggregated energy , compared with the rest regions. This is reasonable in that the traffic within 3 rd Ring Road is commonly more cro wded, and most EVs still hav e suf ficient stored energy by the time of 12 p.m. V I . C O N C L U S I O N A N D F U T U R E W O R K In this paper , we focus on exploring spatial and temporal pattern of the aggregated energy carried by a large number of EVs. A spatio-temporal combined prediction framew ork has been dev eloped to provide predictions of future aggregated energy of multiple city-scale regions. W e hav e extracted both temporal and spatial features from real recorded trajectories data of taxis within Beijing, and studied their influences on the prediction performance. The prediction performance has been demonstrated from both a temporal and spatial perspectiv e. W e ha ve shown that the proposed spatio-temporal predictor can maintain a prediction NMSE below 0 . 1 in all in vestigated regions, and has an improv ed performance compared with both the temporal predictor and the spatial predictor . Furthermore, we hav e demonstrated that the spatio-temporal predictor is capable of providing longer term predictions with low NMSE. February 15, 2018 DRAFT 17 Our future works will utilize the proposed frame work in follo wing aspects, to improv e the interaction between EVs and the smart grid: 1) Utilizing energy pattern prediction for establishing better interactions between EVs and the smart grid. W e plan to utilize the predicted energy pattern to impro ve the operating ef ficiency of the smart grid. Since the aggregated energy can be predicted, the smart grid can perform more flexible scheduling. 2) Utilizing ener gy pattern prediction to benefit EV char ging service providers. More economic ef ficiency can be brought to EV char ging service providers, giv en hourly prediction of the aggregated energy of their residing areas is av ailable. 3) Utilizing energy pattern prediction to improve the quality of EV char ging service. Predicted energy pattern can help EV char ging service providers to schedule and purchase ener gy in adv ance, thus providing better service to EV customers. R E F E R E N C E S [1] Q. W ang, X. Liu, J. Du, and F . K ong, “Smart charging for electric vehicles: A surve y from the algorithmic perspectiv e, ” IEEE Communications Surve ys T utorials , vol. 18, no. 2, pp. 1500–1517, 2016. [2] C. Liu, K. Chau, D. W u, and S. Gao, “Opportunities and challenges of v ehicle-to-home, vehicle-to- vehicle, and vehicle-to-grid technologies, ” Pr oc. IEEE , v ol. 101, no. 11, pp. 2409–2427, Nov 2013. [3] W . Su, H. Eichi, W . Zeng, and M.-Y . Chow , “ A survey on the electrification of transportation in a smart grid en vironment, ” IEEE T rans. Ind. Informat , vol. 8, no. 1, pp. 1–10, Feb 2012. [4] J. Du, S. Ma, Y .-C. W u, and H. V . Poor, “Distrib uted hybrid power state estimation under PMU sampling phase errors, ” IEEE T ransactions on Signal Pr ocessing , vol. 62, no. 16, pp. 4052–4063, 2014. [5] D. T uttle and R. Baldick, “The ev olution of plug-in electric vehicle-grid interactions, ” IEEE T rans. Smart Grid , vol. 3, no. 1, pp. 500–505, March 2012. [6] U.S. Energy Information Administration, “How much electricity does an American home use?” http: //www .eia.gov/tools/f aqs/faq.cfm?id=97&t=3, 2015, [Online; accessed 30-July-2015]. [7] M. Y ilmaz and P . Krein, “Re view of battery charger topologies, charging po wer lev els, and infras- tructure for plug-in electric and hybrid vehicles, ” IEEE T rans. P ower Electr on , v ol. 28, no. 5, pp. 2151–2169, May 2013. [8] R. C. G. II, L. W ang, and M. Alam, “The impact of plug-in hybrid electric vehicles on distribution February 15, 2018 DRAFT 18 networks: A revie w and outlook, ” Rene wable and Sustainable Ener gy Re views , vol. 15, no. 1, pp. 544 – 553, 2011. [9] K. Clement-Nyns, E. Haesen, and J. Driesen, “The impact of charging plug-in h ybrid electric v ehicles on a residential distribution grid, ” IEEE T rans. P ower Syst , vol. 25, no. 1, pp. 371–380, Feb 2010. [10] C. de Fabritiis, R. Ragona, and G. V alenti, “T raf fic estimation and prediction based on real time floating car data, ” in Intelligent T ransportation Systems, 2008. ITSC 2008. 11th International IEEE Confer ence on , Oct 2008, pp. 197–203. [11] J. Y uan, Y . Zheng, C. Zhang, W . Xie, X. Xie, G. Sun, and Y . Huang, “T -driv e: Driving directions based on taxi trajectories, ” in Pr oceedings of the 18th SIGSP ATIAL International Confer ence on Advances in Geographic Information Systems , ser . GIS ’10. Ne w Y ork, NY , USA: A CM, 2010, pp. 99–108. [12] J. Y uan, Y . Zheng, X. Xie, and G. Sun, “Dri ving with kno wledge from the physical world, ” in Pr oceedings of the 17th A CM SIGKDD International Confer ence on Knowledge Discovery and Data Mining , ser . KDD ’11. New Y ork, NY , USA: A CM, 2011, pp. 316–324. [13] Q. W u, A. Nielsen, J. stergaard, S. T . Cha, F . Marra, Y . Chen, and C. Trholt, “Driving pattern analysis for electric vehicle (ev) grid inte gration study , ” in Innovative Smart Grid T echnologies Confer ence Eur ope (ISGT Eur ope), 2010 IEEE PES , Oct 2010, pp. –. [14] L. Zhang, B. Y u, and J. Pan, “Geomob: A mobility-aware geocast scheme in metropolitans via taxicabs and buses, ” in INFOCOM, 2014 Pr oceedings IEEE , April 2014, pp. 1279–1787. [15] T .-K. Lee, Z. Bareket, T . Gordon, and Z. Filipi, “Stochastic modeling for studies of real-world phev usage: Driving schedule and daily temporal distributions, ” IEEE T rans. V ehicular T echnology , vol. 61, no. 4, pp. 1493–1502, May 2012. [16] J. Donadee and M. Ilic, “Stochastic optimization of grid to vehicle frequency regulation capacity bids, ” IEEE T rans. Smart Grid , vol. 5, no. 2, pp. 1061–1069, March 2014. [17] Y . He, B. V enkatesh, and L. Guan, “Optimal scheduling for charging and discharging of electric vehicles, ” IEEE T rans. Smart Grid , vol. 3, no. 3, pp. 1095–1105, Sept 2012. [18] O. Ardakanian, S. K esha v , and C. Rosenber g, “Real-time distributed control for smart electric vehicle chargers: From a static to a dynamic study , ” IEEE T rans. Smart Grid , vol. 5, no. 5, pp. 2295–2305, Sept 2014. [19] E. Sortomme, M. Hindi, S. MacPherson, and S. V enkata, “Coordinated char ging of plug-in h ybrid electric vehicles to minimize distrib ution system losses, ” IEEE T rans. Smart Grid , vol. 2, no. 1, pp. 198–205, March 2011. [20] M. W ang, H. Liang, R. Zhang, R. Deng, and X. Shen, “Mobility-aware coordinated charging for February 15, 2018 DRAFT 19 electric vehicles in vanet-enhanced smart grid, ” IEEE J. Sel. Ar eas in Commun , vol. 32, no. 7, pp. 1344–1360, July 2014. [21] Y . Li and W . W ang, “Horizon on the mo ve: Geocast in intermittently connected v ehicular ad hoc networks, ” in INFOCOM, 2013 Pr oceedings IEEE , April 2013, pp. 2553–2561. [22] J. Du, X. Liu, and L. Rao, “Proactiv e doppler shift compensation in vehicular cyber -physical systems, ” IEEE/A CM T ransactions on Networking, arXiv pr eprint arXiv:1710.00778 , v ol. PP , no. 99, pp. 1–12, 2018. [23] J. K enney , “Dedicated short-range communications (dsrc) standards in the united states, ” Pr oc. IEEE , vol. 99, no. 7, pp. 1162–1182, July 2011. [24] A. Eltahir , R. Saeed, and R. Mokhtar , “V ehicular communication and cellular network inte gration: Gate way selection perspecti ve, ” in Computer and Communication Engineering (ICCCE), 2014 International Conference on , Sept 2014, pp. 64–67. [25] J. Du and Y .-C. W u, “Network-wide distributed carrier frequency of fsets estimation and compensation via belief propagation, ” IEEE T ransactions on Signal Pr ocessing , vol. 61, no. 23, pp. 5868–5877, 2013. [26] S. Hu, H. Liu, L. Su, H. W ang, T . Abdelzaher , P . Hui, W . Zheng, Z. Xie, and J. Stanko vic, “T ow ards automatic phone-to-phone communication for vehicular networking applications, ” in INFOCOM, 2014 Pr oceedings IEEE , April 2014, pp. 1752–1760. [27] J. Du and Y .-C. W u, “Fully distrib uted clock ske w and of fset estimation in wireless sensor networks, ” in IEEE International Confer ence on Acoustics, Speech and Signal Pr ocessing (ICASSP) , 2013, pp. 4499–4503. [28] H. Liang and W . Zhuang, “Double-loop receiv er-initiated mac for cooperativ e data dissemination via roadside wlans, ” IEEE T rans. Communications , vol. 60, no. 9, pp. 2644–2656, September 2012. [29] J. Du, S. Ma, Y . C. W u, and H. V . Poor , “Distributed bayesian hybrid power state estimation with pmu synchronization errors, ” in 2014 IEEE Global Communications Confer ence , 2014, pp. 3174–3179. [30] Y . Zheng, “T -drive trajectory data sample, ” August 2011. [Online]. A vailable: http://research. microsoft.com/apps/pubs/default.aspx?id=152883 [31] H. Holden and N. H. Risebro, “ A mathematical model of traffic flow on a network of unidirectional roads, ” SIAM Journal on Mathematical Analysis , vol. 26, no. 4, pp. 999–1017, 1995. [32] S. L ¨ ammer and D. Helbing, “Self-control of traf fic lights and vehicle flo ws in urban road networks, ” J ournal of Statistical Mechanics: Theory and Experiment , vol. 2008, no. 04, p. P04019, 2008. [33] J. Du and Y .-C. W u, “Distributed clock ske w and of fset estimation in wireless sensor networks: February 15, 2018 DRAFT 20 Asynchronous algorithm and con ver gence analysis, ” IEEE T rans. W ireless Commun. , vol. 12, no. 11, 2013. [34] H. Dia, “ An object-oriented neural network approach to short-term traffic forecasting, ” Eur opean J ournal of Operational Resear ch , vol. 131, no. 2, pp. 253 – 261, 2001, artificial Intelligence on T ransportation Systems and Science. [35] W . Zheng, D.-H. Lee, and Q. Shi, “Short-term free way traffic flow prediction: Bayesian combined neural network approach, ” J ournal of transportation engineering , vol. 132, no. 2, pp. 114–121, 2006. [36] Y . Zheng, F . Liu, and H.-P . Hsieh, “U-air: When urban air quality inference meets big data, ” in Pr oceedings of the 19th A CM SIGKDD International Confer ence on Knowledge Discovery and Data Mining , ser . KDD ’13. New Y ork, NY , USA: A CM, 2013, pp. 1436–1444. [37] H. Liang, B. J. Choi, W . Zhuang, and X. Shen, “Optimizing the energy deliv ery via v2g systems based on stochastic in ventory theory , ” IEEE T rans. Smart Grid , vol. 4, no. 4, pp. 2230–2243, Dec 2013. [38] J. Du, S. Ma, Y .-C. W u, S. Kar , and J. M. Moura, “Conv ergence analysis of the information matrix in Gaussian belief propagation, ” 2017 IEEE International Confer ence on Acoustics, Speech and Signal Pr ocessing (ICASSP), arXiv pr eprint arXiv:1704.03969 , 2017. [39] J. Du, S. Kar , and J. M. Moura, “Distributed con vergence verification for Gaussian belief prop- agation, ” Asilomar Confer ence on Signals, Systems and Computers (ASILOMAR), arXiv pr eprint arXiv:1711.09888 , 2017. [40] J. Lafferty , A. McCallum, and F . C. Pereira, “Conditional random fields: Probabilistic models for segmenting and labeling sequence data, ” June. 2001. [41] J. Du, S. Ma, Y .-C. W u, S. Kar , and J. M. Moura, “Con ver gence analysis of belief propagation for pairwise linear Gaussian models, ” 2017 IEEE Global Confer ence on Signal and Information Pr ocessing (GlobalSIP), arXiv pr eprint arXiv:1706.04074 , 2017. [42] J. Du, S. Ma, Y .-C. W u, S. Kar , and J. M. Moura, “Con ver gence analysis of distributed inference with vector -v alued Gaussian belief propagation, ” arXiv pr eprint arXiv:1611.02010 , 2016. February 15, 2018 DRAFT 21 T ABLE II P R E D I C T I O N P E R F O R M A N C E B A S E D O N D I FF E R E N T F E AT U R E S E T S . Feature sets (SP) NMSE ave SP NMSE min SP NMSE max SP F V , F N , F E 0 . 2311 0 . 9322 0 . 0305 F D , F V , F N , F E 0 . 1585 0 . 3541 0 . 0367 F E , F N 0 . 1570 0 . 3215 0 . 0653 F N 0 . 1450 0 . 3588 0 . 0615 F E 0 . 1247 0 . 2593 0 . 0359 F D , F N , F E 0 . 1185 0 . 1851 0 . 0649 Feature sets (STP) NMSE ave STP NMSE min STP NMSE max STP F V , F N , F E 0 . 0800 0 . 3265 0 . 0124 F D , F V , F N , F E 0 . 0558 0 . 1310 0 . 0211 F E , F N 0 . 0675 0 . 1274 0 . 0369 F N 0 . 0446 0 . 0825 0 . 0186 F E 0 . 0459 0 . 0897 0 . 0171 F D , F N , F E 0 . 0485 0 . 0703 0 . 0254 February 15, 2018 DRAFT