机器学习在电能需求预测中的应用

问题：电力公司会提前一天或者几天跟发电厂购买电量，买多了势必会造成不必要的浪费，而买少了临时向发电厂买电代价是比平时高得多的，如何根据以往的数据对未来的电量需求做出精准的预测，是电力公司非常重视的一个问题。

解决方案：通过历年来该地区的实际电量需求数据和该地区的天气数据进行建模，建立一个回归模型，预测未来几天的用电需求。

数据来源：美国纽约州电力NYISO公司公布了历年来纽约州各个地区的实际用电数据，该公司还公布了历年来提前一天的电量预测数据，因此可以从http://mis.nyiso.com/public/P-58Blist.htm（需要VPN）下载历年的实际用电量loads，从http://www.nyiso.com/public/markets_operations/market_data/custom_report/index.jsp?report=load_forecast下载NYISO公司提前一天的电能预测数据。这里我们只是获取了电量的使用情况，按建模的话说也就是只知道标签Y，特征X就需要从天气数据中获取，天气数据从https://www.wunderground.com/history/airport/获取。

1. 数据抓取：首先是抓取电量数据loads，由于NYISO网站上的数据是zip，因此可以通过python的urllib库进行抓取下载到本地：

import urllibimport pandas as pd#首先通过pandas的date_range产出日期格式，再通过遍历改变url的形式逐个抓取数据，通过urllib.urlretrieve下载zip数据dates=pd.date_range(pd.to_datetime('2005-01-01'),pd.to_datetime('2016-12-31'),freq='M')for date in dates:    url='http://mis.nyiso.com/public/csv/pal/{0}{1}01pal_csv.zip'.format(date.year,  str(date.month).zfill(2))     urllib.urlretrieve(url, "{0}".format(url.split('/')[-1]))    print date

通过上述操作，把2005年1月1日到16年12月31日的所有数据下载到本地，下载下来的数据为zip格式，因此通过zipfile库来进行解压：

import osimport zipfileimport pandas as pd#通过os.listdir打开指定目录所有文件，在通过if筛选所有的zip文件，把文件名（不带路径）都append到一个list中zips=[]for file in os.listdir("../data/zip/"):    if file.endswith(".zip"):       zips.append(file)

建立一个解压函数

def unzip(fileName,fileDir):    with zipfile.ZipFile(fileName) as zf:         try:             zf.extractall(fileDir)         except:             return

遍历zips文件名list，逐个解压到fileDir(“../data/zip/csvData/”)文件夹下

for z in zips:    try:        unzip("../data/zip/"+z, "../data/zip/csvData/")    except:        print "../data/zip/csvData/"        continue#把所有的csv文件合并到同一个csv文件中csvs=[] for file in os.listdir("../data/zip/csvData/"):    if file.endswith("pal.csv"):#合并到同一个csv，由于每一个csv文件都有头部信息，因此先读取第一个csv的头部信息作为总的csv的头部信息，剩下的直接跳过fout=open("../data/combine.csv","a")for line in open("../data/zip/csvData/"+csvs[0]):    fout.write(line)for file in csvs[1:]:    f=open("../data/zip/csvData/"+file)    f.next()#跳过头部信息,从第二行开始    for line in f:        fout.write(line)    f.close()fout.close()

按地区保存数据，由于下载下来的数据分为好几个地区，因此需要进行地区分组

df=pd.read_csv("../data/combine.csv")cols=df.columnsdf.columns=[col.lower().replace(' ','') for col in cols] #转小写，去空格df=df[['timestamp','name','ptid','load']]#时间戳，地区名，ID，用电量#由于天气数据地区（台站名）和电量数据地区（地区名）不一致，因此需要建立一个字典，把它们联系起来（地区包含台站名）regionsName=list(df.name.unique()) #地区名weatherName= ['kalb', 'ksyr', 'klga', 'kroc', 'kpou', 'kjfk', 'krme', 'klga', 'kjfk', 'kpbg', 'kbuf']***#本文为了简便，就不在处理所有地区，这里只选择 地区：north  台站：kjfk***#因此电量数据只保存north地区的数据subset=df[df.name=='north'].copy()filename='north.csv'subset.to_csv("../data"+filename,index=False)#想要完整保存所有地区数据的操作如下："""df.name.unique() #每个地区，每个地区又有不同的城市##建立地区和城市字典regions = list(df.name.unique())region_names = ['Capital', 'Central', 'Dunwoodie', 'Genese', 'Hudson Valley', 'Long Island', 'Mohawk Valley', 'Millwood', 'NYC', 'North', 'West']cities = ['Albany', 'Syracuse', 'Yonkers', 'Rochester', 'Poughkeepsie', 'NYC', 'Utica', 'Yonkers', 'NYC', 'Plattsburgh', 'Buffalo']weather_stations = ['kalb', 'ksyr', 'klga', 'kroc', 'kpou', 'kjfk', 'krme', 'klga', 'kjfk', 'kpbg', 'kbuf']weather_dict = dict(zip(regions, zip(weather_stations, region_names, cities)))import cPickle as picklepickle.dump(weather_dict, open('weather_dict.pkl','wb'))for region in weather_dict.keys():    subset = df[df.name == region].copy()    filename = weather_dict[region][1].lower().replace(' ', '') + '.csv'    subset.to_csv("../data" + filename, index=False)"""           

2.至此，电量的数据已经下载保存完毕，以csv格式、分地区保存，下面来抓取天气数据：
天气的数据不能像电量数据一样直接下载，它是直接显示在html网页上的，但可以通过?format=1把html网页转为表格数据（只显示html中的表格数据），再通过pandas的read_csv来读取到本地

import pandas pdimport os#同一产出时间dates=pd.date_range(pd.to_datetime('2005-01-01'),pd.to_datetime('2016-12-31'),freq='D')#定义一个抓取weather的函数，传入地区和时间，返回weather数据def getWeather(WeatherStation,dates):    for d in dates:        try:  df=pd.read_csv('https://www.wunderground.com/history/airport/{0}/{1}/{2}/{3}/DailyHistory.html?format=1'.format(WeatherStation, d.year, d.month, d.day)) #自行按照Python格式对其了，这里网页显示问题            cols=df.columns            df.columns=[clo.lower().replace(' ','').replace('<br/>','')for col in clos]            df.dateutc=df.dateutc.apply(lambda x: pd.to_datetime(x.erplace('<br />','')))            #不同地区的格式有点不一致            try:                df.gustspeedmph = df.gustspeedmph.replace('-', 0)                df.windspeedmph = df.windspeedmph.replace('Calm', 0)                df.precipitationin = df.precipitationin.replace('NaN', 0)                df.events = df.events.replace('NaN', 0)                print '1111'            except:                df.gustspeedkmh = df.gustspeedkmh.replace('-', 0)                df.windspeedkmh = df.windspeedkmh.replace('Calm', 0)                df.precipitationmm= df.precipitationmm.replace('NaN', 0)                df.events = df.events.replace('NaN', 0)                print '2222'           filepath="../data/weather/"+WeatherStation+str(d.date()).replace('-','') + '.csv'           df.to_csv(filepath,index=False)           if type(df.dateutc[0])==pd.tslib.Timestamp:              continue           else:              print 'wrong!!!'              break     except:            print "date ",d ," can't be downloaded!"            continue print "Files for %s have been written" % WeatherStation return#调用函数，获取weather数据  ***#本文为了简便，就不在处理所有地区，这里只选择 地区：north  台站：kjfk***getWeather('kjfk',dates) #只抓取一个地区的数据，如果想抓所有地区，可写一个for遍历weather_dict

合并数据并按地区存储

WeatherStation='kjfk'csvs=[]for file in os.listdir("../data/weather"):    if file.endswith(".csv"):       csvs.append(file)fout=open("../data/weather/"+str(WeatherStation)+"_all.csv","a")#同上操作，先第一个文件读头部信息，后续文件跳过for line in open("../data/weather/"+csvs[0]):    fout.write(line)for file in csvs[1:]:    f=open("../data/weather/"+file)    f.next()    for line in f:        fout.write(line)    f.close()fout.close()print "have been done %s"%WeatherStation            

3.至此，天气数据已经下载合并分区完毕，下面进行天气数据电量数据的合并，以及天气数据特征提取
这里需要注意的是，两者的时间格式是不同的，先进行格式统一，最后再通过时间来进行合并数据。电量数据是每隔5分钟一组数据，而天气是每小时一组，因此通过sklearn的KNN来进行差值。

首先先对时间格式进行处理：

import pandas as pdfrom sklearn.neighbors import NearestNeighbors#台站kjfk对应地区northloads = pd.read_csv("../data/north.csv") #电能需求Yweather = pd.read_csv("../data/weather/kjfk_all.csv")#天气Xweather['date']=weather.dateutc.apply(lambda x: pd.to_datetime(x).date())weather['timeest']=weather.timeest.apply(lambda x:pd.to_datetime(x).time())f=weather[['date','timeest']].astype(str)weather['timestamp']=pd.to_datetime(f['date']+' '+f['timeest'])loads['timestamp']=loads.timestamo.apply(lambda x: pd.to_datetime(x))

用knn进行差值

nbrs=NearestNeighbors(1).fit(weather['timestamp']).values[:,None]#这里1表示参数n_neighbors个数dist,ind=nbrs.kneighbors(loads['timestamp'].values[:,None])#ind和dist都是一个一维向量loads['nearesttime']=weather['timestamp'].values[ind.ravel()]full=loads.merge(weather,left_on='nearesttime',right_on='timestamp')#合并#去掉无关列，重命名full = full[['timestamp_x', 'load', 'nearesttime', 'temperaturec',                 'dewpointc', 'humidity', 'sealevelpressurehpa', 'winddirection', 'windspeedkmh',                'precipitationmm']].rename(columns={'timestamp_x': 'timestamp', 'nearesttime':'weathertime'})

构建weather特征，根据这篇文章来提取天气特征的时间特征https://arxiv.org/pdf/1506.06972.pdf

def get_prev_days(x, n_days):    '''先前推n天的电能需求数据'''    try:        lo = full[full.timestamp == x - n_days*pday].load.values[0]    except:        lo = full[full.timestamp == x].load.values[0]    return lodef get_prev_hours(x, n_hours):    '''n小时前的电能需求'''    try:        lo = full[full.timestamp == x - n_hours*phour].load.values[0]    except:        lo = full[full.timestamp == x].load.values[0]    return lofull = full.copy()full['dow'] = full.timestamp.apply(lambda x: x.dayofweek) #一周中的一天full['doy'] = full.timestamp.apply(lambda x: x.dayofyear) #一年中的一天full['day'] = full.timestamp.apply(lambda x: x.day) #每天full['month'] = full.timestamp.apply(lambda x: x.month) #每月full['year'] = full.timestamp.apply(lambda x: x.year) #每年full['hour'] = full.timestamp.apply(lambda x: x.hour) #每小时full['minute'] = full.timestamp.apply(lambda x: x.hour*60 + x.minute)#每分full['t_m24'] = full.timestamp.apply(get_prev_days, args=(1,))#向前推1天的电能需求数据full['t_m48'] = full.timestamp.apply(get_prev_days, args=(2,))# #两天full['tdif'] = full['load'] - full['t_m24']#今天和昨天的区别full['t_m1'] = full.timestamp.apply(get_prev_hours, args=(1,))#先前推1小时的电能需求（这里计算非常耗时，大概15个小时）full.to_csv('full_{0}_features.csv'.format('north'), index=False)#数据保存

4.终于把数据采集、清洗、特征提取给搞完了，接下来建模。嗯，建模！！！