RealTimeSimulator_LoadForecasting.py 29 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655
  1. # # Day-ahead load forecasting
  2. #
  3. # DB : MS SQL
  4. #
  5. # Program Language : Python
  6. #
  7. # kgpark@hdc-icontrols.com
  8. # April 10, 2020
  9. # ### BEMS 데이터 수집 메카니즘
  10. # #### 데이터 별로 수집 타입에 따라 다르지만, Raw 테이블에 적산 값으로 저장이 되고 15min 테이블에서 해당 시간대와 그 전 시간대의 차이 값을 입력한다.
  11. # #### DGW 혹은 시스템에 이상이 생겼을 때, 데이터가 들어오지 않거나 0으로 입력된다.
  12. # #### 1시간 테이블은 15분 테이블에서 각 15분, 30분, 45분, 60분의 데이터 합산 값이 나왔다.
  13. # #### 합산 값으로 저장되다보니 4개 포인트 중 적어도 하나만 있어도 1시간 데이터로 저장이 된다.
  14. # #### 따라서, 15분 데이터를 전처리하는 것이 주효하고 데이터가 없거나 0값을 검출하여 비정상 데이터로 추정하는 것을 추천한다.
  15. # #### 또한, 1시간 단위로 데이터 주기를 변환한다면 15분 테이블의 4개 포인트 중 하나라도 값을 모른다면 그 시간의 데이터가 비정상이라고 가정하는 것을 추천한다.
  16. import matplotlib.pyplot as plt
  17. import pymssql
  18. import datetime
  19. import numpy as np
  20. import math
  21. from korean_lunar_calendar import KoreanLunarCalendar
  22. import calendar
  23. import configparser
  24. import sys
  25. import time
  26. # ## Define functions
  27. ### Define day-type
  28. def getDayName(year, month, day):
  29. return ['MON','TUE','WED','THU','FRI','SAT','SUN'][datetime.date(year, month, day).weekday()]
  30. def getDayType(DateinDay, Period, SpecialHoliday):
  31. DoW=[]; # Day of Week
  32. for i in range(Period):
  33. if DateinDay[i].year==2019 and DateinDay[i].month==5 and DateinDay[i].day==18:
  34. DoW.append([5, DateinDay[i]])
  35. elif getDayName(DateinDay[i].year,DateinDay[i].month,DateinDay[i].day) == 'MON':
  36. DoW.append([1, DateinDay[i]])
  37. elif getDayName(DateinDay[i].year,DateinDay[i].month,DateinDay[i].day) == 'TUE':
  38. DoW.append([2, DateinDay[i]])
  39. elif getDayName(DateinDay[i].year,DateinDay[i].month,DateinDay[i].day) == 'WED':
  40. DoW.append([3, DateinDay[i]])
  41. elif getDayName(DateinDay[i].year,DateinDay[i].month,DateinDay[i].day) == 'THU':
  42. DoW.append([4, DateinDay[i]])
  43. elif getDayName(DateinDay[i].year,DateinDay[i].month,DateinDay[i].day) == 'FRI':
  44. DoW.append([5, DateinDay[i]])
  45. elif getDayName(DateinDay[i].year,DateinDay[i].month,DateinDay[i].day) == 'SAT':
  46. DoW.append([6, DateinDay[i]])
  47. elif getDayName(DateinDay[i].year,DateinDay[i].month,DateinDay[i].day) == 'SUN':
  48. DoW.append([7, DateinDay[i]])
  49. for j in range(len(SpecialHoliday)):
  50. if SpecialHoliday[j] == datetime.date(DateinDay[i].year,DateinDay[i].month,DateinDay[i].day):
  51. DoW[-1][0] = 8
  52. break
  53. ### W:1, N:2, ### W: Workday, N: Non-workday
  54. DayType=[]
  55. for i in range(Period):
  56. if DoW[i][0] <= 5:
  57. DayType.append([1, DateinDay[i]])
  58. elif DoW[i][0] > 5:
  59. DayType.append([2, DateinDay[i]])
  60. return DoW, DayType
  61. def Reconstruction(DayType, DatainHour, mark, DataRes, isRecent):
  62. ReconstructedData=[]
  63. DayType1h=[]
  64. Day_len = len(DayType)
  65. # Rearrange data in hour unit
  66. for i in range(Day_len):
  67. if i == Day_len-1 and isRecent:
  68. Time_len = len(DatainHour) - i*DataRes
  69. else:
  70. Time_len=DataRes
  71. for j in range(Time_len):
  72. DayType1h.append([DatainHour[i*DataRes + j], DayType[i][0], datetime.datetime(DayType[i][1].year, DayType[i][1].month, DayType[i][1].day, j, 0)]) ## data, daytype, time
  73. # 비정상 데이터보다 앞선 시간의 데이터 중 DayType이 같고 시간이 같은 5개 날 데이터의 평균으로 복원함
  74. for i in reversed(range(len(DayType1h))):
  75. AccData=[]
  76. cnt=0
  77. if math.isnan(DayType1h[i][0]):
  78. for j in range(len(DayType1h)):
  79. if cnt > 5:
  80. break
  81. if i < j and DayType1h[j][1] == DayType1h[i][1] and DayType1h[j][2].hour == DayType1h[i][2].hour and (not math.isnan(DayType1h[j][0])):
  82. AccData.append(DayType1h[j][0])
  83. cnt += 1
  84. DayType1h[i][0] = np.mean(AccData)
  85. ReconstructedData.append(DayType1h[i][0])
  86. ReconstructedData.reverse()
  87. ### Double-checking for the data which is not reconstructed, especially in front
  88. for i in range(len(DayType1h)):
  89. AccData=[]
  90. cnt=0
  91. if math.isnan(DayType1h[i][0]):
  92. #print('Here is NaN!!',ReconstructedData[i],i,DayType1h[i][2].hour, DayType1h[i][1])
  93. for j in reversed(range(len(DayType1h))):
  94. if cnt > 5:
  95. break
  96. if i > j and DayType1h[j][1] == DayType1h[i][1] and DayType1h[j][2].hour == DayType1h[i][2].hour and (not math.isnan(DayType1h[j][0])):
  97. AccData.append(DayType1h[j][0])
  98. cnt += 1
  99. ReconstructedData[i] = np.mean(AccData)
  100. return ReconstructedData, DayType1h
  101. ## For day-ahead linear prediction
  102. def lpc_pred_DayAhead(Data_trn, DayType_trn, cov_lth, DayType_tst, DataRes):
  103. # Calculating the filter bank for each hour and day-type using traing set
  104. for c_w in range(1,3):
  105. DayType_trn[0,0]=0
  106. CP_pred_fb=np.zeros(Data_trn.shape)
  107. lpc_fb=np.zeros([cov_lth[c_w-1],DataRes])
  108. Prv_A=[]
  109. Prv_A=np.transpose(Data_trn[DataRes-cov_lth[c_w-1]:DataRes,np.where(DayType_trn == c_w)[0]-1])
  110. for hr_i in range(24):
  111. lpc_fb[:,hr_i]=np.dot(np.linalg.pinv(Prv_A), np.transpose(Data_trn[hr_i,np.where(DayType_trn == c_w)[0]]))
  112. if c_w == 1:
  113. lpc_fb1=lpc_fb
  114. elif c_w == 2:
  115. lpc_fb2=lpc_fb
  116. ## For testing
  117. if DayType_tst[0,0] == 1:
  118. lpc_t=lpc_fb1
  119. elif DayType_tst[0,0] == 2:
  120. lpc_t=lpc_fb2
  121. Data_tt=Data_trn[:,-1]
  122. # Load prediction for test day based on the filter bank
  123. CP_pred=np.transpose(np.dot(np.transpose(Data_tt[DataRes-cov_lth[DayType_tst[0,0]-1]:DataRes+1]),lpc_t))
  124. return CP_pred
  125. ## For step-ahead linear prediction
  126. def lpc_pred_OneStepAhead(Data_trn, DayType_trn, cov_lth, DayType_tst, DataRes):
  127. for c_w in range(1,3):
  128. DayType_trn[0,0]=0
  129. lpc_fb=np.zeros([cov_lth[c_w-1],DataRes])
  130. Prv_A=[]
  131. Prv_A=np.transpose(Data_trn[DataRes-cov_lth[c_w-1]:DataRes,np.where(DayType_trn == c_w)[0]-1])
  132. lpc_fb=np.dot(np.linalg.pinv(Prv_A), np.transpose(Data_trn[0,np.where(DayType_trn == c_w)[0]]))
  133. if c_w == 1:
  134. lpc_fb1=lpc_fb
  135. elif c_w == 2:
  136. lpc_fb2=lpc_fb
  137. ## Testing
  138. if DayType_tst[0,0] == 1:
  139. lpc_t=lpc_fb1
  140. elif DayType_tst[0,0] == 2:
  141. lpc_t=lpc_fb2
  142. Data_tt=Data_trn[:,-1]
  143. CP_pred=np.transpose(np.dot(np.transpose(Data_tt[DataRes-cov_lth[DayType_tst[0,0]-1]:DataRes+1]),lpc_t))
  144. return CP_pred
  145. ## Measure
  146. def MAPE(y_observed, y_pred):
  147. return np.mean(np.abs((y_observed - y_pred) / y_observed)) * 100
  148. def MAE(y_observed, y_pred):
  149. return np.mean(np.abs(y_observed - y_pred))
  150. def MBE(y_observed, y_pred):
  151. return (np.sum((y_observed - y_pred))/(len(y_observed)*np.mean(y_observed)))*100
  152. def CVRMSE(y_observed, y_pred):
  153. return (np.sqrt(np.mean((y_observed - y_pred)*(y_observed - y_pred)))/np.mean(y_observed))*100
  154. ## Check for normal time stamp
  155. def Check_AlivedTimeStamp(RawData, ComparedData, idx_raw, idx_comp):
  156. if datetime.date(RawData[idx_raw][4].year,RawData[idx_raw][4].month,RawData[idx_raw][4].day) == datetime.date(ComparedData[idx_comp].year, ComparedData[idx_comp].month, ComparedData[idx_comp].day) and datetime.time(RawData[idx_raw][4].hour,RawData[idx_raw][4].minute) == datetime.time(ComparedData[idx_comp].hour, ComparedData[idx_comp].minute):
  157. isAlived = True
  158. else:
  159. isAlived = False
  160. return isAlived
  161. if __name__ == "__main__" :
  162. ## Check every hour on the hour operating infinite loop
  163. while True:
  164. now = datetime.datetime.now().now()
  165. ## distinguish real time update and specific day
  166. ## 자정에 생기는 인덱싱 문제로 0시에는 16분에 업데이트
  167. if (now.hour != 0 and now.minute == 1) or (now.hour == 0 and now.minute == 16):
  168. PredctionActive = True
  169. else:
  170. PredctionActive = False
  171. if now.second > 55:
  172. print("[ Current Time -", now.hour,":", now.minute,":", now.second,"], " "Sleeping for 30 seconds... Prediction starts every hour")
  173. time.sleep(30)
  174. else:
  175. print("[ Current Time -", now.hour,":", now.minute,":", now.second,"], " "Sleeping for 60 seconds... Prediction starts every hour")
  176. time.sleep(60)
  177. if PredctionActive:
  178. ## Loading .ini file
  179. myINI = configparser.ConfigParser()
  180. myINI.read("Config.ini", "utf-8" )
  181. # MSSQL Access
  182. conn = pymssql.connect(host=myINI.get('LocalDB_Info','ip_address'), user=myINI.get('LocalDB_Info','user_id'), password=myINI.get('LocalDB_Info','user_password'), database=myINI.get('LocalDB_Info','db_name'), autocommit=True)
  183. # Create Cursor from Connection
  184. cursor = conn.cursor()
  185. # Execute SQL (Electric consumption)
  186. cursor.execute('SELECT * FROM BemsConfigData where SiteId = 1')
  187. rowDB_info = cursor.fetchone()
  188. conn.close()
  189. loadDBIP = rowDB_info[1]
  190. loadDBUserID = rowDB_info[2]
  191. loadDBUserPW = rowDB_info[3]
  192. loadDBName = rowDB_info[4]
  193. targetDBIP = rowDB_info[5]
  194. targetDBUserID = rowDB_info[6]
  195. targetDBUserPW = rowDB_info[7]
  196. targetDBName = rowDB_info[8]
  197. linearFilterLength = '24,24'
  198. print("=================== Prediction start! ===================")
  199. startday = datetime.date(2019,1,1)
  200. # ## Data accumulation
  201. isRecent = True
  202. lastday = datetime.date(now.year, now.month, now.day)
  203. if startday < datetime.date(2017,1,1):
  204. print('[ERROR] 데이터 최소 시작 시점은 2017.01.01 입니다')
  205. elif startday > lastday:
  206. print('[ERROR] 예측 타깃 시작시점이 데이터 시작 시점보다 작을 수 없습니다')
  207. now_ = datetime.date(now.year, now.month, now.day)
  208. # 학습데이터의 기간은 최대 2년으로 한정
  209. if (startday-now_).days > 730:
  210. Ago_2year = now_ + timedelta(days=-730)
  211. startday = datetime.date(Ago_2year.year, Ago_2year.month, Ago_2year.day)
  212. # MSSQL Access
  213. conn = pymssql.connect(host = loadDBIP, user = loadDBUserID, password = loadDBUserPW, database = loadDBName, autocommit=True)
  214. # Create Cursor from Connection
  215. cursor = conn.cursor()
  216. # Execute SQL (Electric consumption)
  217. cursor.execute('SELECT * FROM BemsMonitoringPointHistory15min where SiteId = 1 and FacilityTypeId = 99 and FacilityCode = 4863 and PropertyId = 1 order by CreatedDateTime desc')
  218. # 데이타 하나씩 Fetch하여 출력
  219. row = cursor.fetchone()
  220. DataRes_org=96
  221. DataRes_24=24
  222. rawData=[]
  223. while row:
  224. row = cursor.fetchone()
  225. if datetime.date(row[4].year,row[4].month,row[4].day) < startday:
  226. break
  227. rawData.append(row)
  228. rawData.reverse() # 오름차순 정렬
  229. # 연결 끊기
  230. conn.close()
  231. # 현장 데이터가 없을 경우 예외처리
  232. if now.hour == 0:
  233. hour_calib = 0
  234. else:
  235. hour_calib = 1
  236. if datetime.datetime(now.year, now.month, now.day, now.hour, 0, 0) - datetime.timedelta(hours=hour_calib) == datetime.datetime(rawData[-1][4].year, rawData[-1][4].month, rawData[-1][4].day, rawData[-1][4].hour, 0, 0):
  237. # MSSQL Access
  238. conn = pymssql.connect(host = loadDBIP, user = loadDBUserID, password = loadDBUserPW, database = loadDBName, autocommit=True)
  239. # Create Cursor from Connection
  240. cursor = conn.cursor()
  241. # SQL문 실행 (정기휴일)
  242. cursor.execute('SELECT * FROM CmHoliday where SiteId = 1 and IsUse = 1')
  243. # 데이타 하나씩 Fetch하여 출력
  244. row = cursor.fetchone()
  245. regularHolidayData = [row]
  246. while row:
  247. row = cursor.fetchone()
  248. regularHolidayData.append(row)
  249. regularHolidayData = regularHolidayData[0:-1]
  250. # SQL문 실행 (비정기휴일)
  251. cursor.execute('SELECT * FROM CmHolidayCustom where SiteId = 1 and IsUse = 1')
  252. # 데이타 하나씩 Fetch하여 출력
  253. row = cursor.fetchone()
  254. suddenHolidayData = [row]
  255. while row:
  256. row = cursor.fetchone()
  257. suddenHolidayData.append(row)
  258. suddenHolidayData = suddenHolidayData[0:-1]
  259. # 연결 끊기
  260. conn.close()
  261. # 공휴일의 음력 계산
  262. calendar_convert = KoreanLunarCalendar()
  263. SpecialHoliday = []
  264. for i in range(lastday.year-startday.year+1):
  265. for j in range(len(regularHolidayData)):
  266. if regularHolidayData[j][3] == 1:
  267. if regularHolidayData[j][1] == 12 and regularHolidayData[j][2] == 30: ## 설 하루 전 연휴 계산을 위함
  268. calendar_convert.setLunarDate(startday.year+i-1, regularHolidayData[j][1], regularHolidayData[j][2], False)
  269. SpecialHoliday.append(datetime.date(int(calendar_convert.SolarIsoFormat().split(' ')[0].split('-')[0]), int(calendar_convert.SolarIsoFormat().split(' ')[0].split('-')[1]), int(calendar_convert.SolarIsoFormat().split(' ')[0].split('-')[2])))
  270. else:
  271. calendar_convert.setLunarDate(startday.year+i, regularHolidayData[j][1], regularHolidayData[j][2], False)
  272. SpecialHoliday.append(datetime.date(int(calendar_convert.SolarIsoFormat().split(' ')[0].split('-')[0]), int(calendar_convert.SolarIsoFormat().split(' ')[0].split('-')[1]), int(calendar_convert.SolarIsoFormat().split(' ')[0].split('-')[2])))
  273. else:
  274. SpecialHoliday.append(datetime.date(startday.year+i,regularHolidayData[j][1],regularHolidayData[j][2]))
  275. for i in range(len(suddenHolidayData)):
  276. if suddenHolidayData[i][1].year >= startday.year:
  277. SpecialHoliday.append(datetime.date(suddenHolidayData[i][1].year, suddenHolidayData[i][1].month, suddenHolidayData[i][1].day))
  278. SpecialHoliday=list(set(SpecialHoliday))
  279. DayPeriod = (lastday - startday).days + 1
  280. print('First day:',startday,',', 'Last Day:', lastday,',','Current Time:', now)
  281. print('Day period :', DayPeriod)
  282. # ## Find unkown/zero data (Bad data)
  283. StartTime = datetime.datetime(int(startday.strftime('%Y')), int(startday.strftime('%m')), int(startday.strftime('%d')), 0, 0, 0)
  284. TimeStamp_DayUnit = []
  285. StandardTimeStamp = []
  286. # Create normal time stamp
  287. for idx_day in range(DayPeriod):
  288. TimeStamp_DayUnit.append(startday + datetime.timedelta(days=idx_day))
  289. if isRecent and idx_day == DayPeriod-1:
  290. if now.hour == 0: # 예외처리용 (자정에 Day count가 안되는 현상)
  291. tmp_len = 1
  292. else:
  293. tmp_len = now.hour*4 + int(now.minute/15)
  294. for idx_time in range(tmp_len):
  295. StandardTimeStamp.append(StartTime)
  296. StartTime += datetime.timedelta(minutes = 15)
  297. else:
  298. for idx_time in range(DataRes_org):
  299. StandardTimeStamp.append(StartTime)
  300. StartTime += datetime.timedelta(minutes = 15)
  301. RawDate=[] # raw data (date)
  302. RawElectricLoad=[] # raw data (electric load)
  303. for i in range(len(rawData)):
  304. if datetime.date(rawData[i][4].year,rawData[i][4].month,rawData[i][4].day) >= startday:
  305. if datetime.date(rawData[i][4].year,rawData[i][4].month,rawData[i][4].day) <= lastday:
  306. RawDate.append(rawData[i][4])
  307. RawElectricLoad.append(rawData[i][5])
  308. if datetime.date(rawData[i][4].year,rawData[i][4].month,rawData[i][4].day) > lastday:
  309. break
  310. Data_len=len(RawDate)
  311. if isRecent:
  312. DataAct_len = (DayPeriod-1)*DataRes_org + now.hour*4 + int(now.minute/15)
  313. else:
  314. DataAct_len = DayPeriod*DataRes_org
  315. ### Unknown/zero data counts
  316. DataCount=[]
  317. for i in range(len(TimeStamp_DayUnit)):
  318. cnt_unk=0 # For Unknown data count
  319. cnt_zero=0 # zero data count
  320. for j in range(Data_len):
  321. if TimeStamp_DayUnit[i] == datetime.date(RawDate[j].year,RawDate[j].month,RawDate[j].day):
  322. cnt_unk += 1
  323. if RawElectricLoad[j] == 0:
  324. cnt_zero += 1
  325. if isRecent and i==len(TimeStamp_DayUnit)-1:
  326. DataCount.append([TimeStamp_DayUnit[i], now.hour*4 + int(now.minute/15) - cnt_unk, cnt_zero])
  327. else:
  328. DataCount.append([TimeStamp_DayUnit[i], DataRes_org-cnt_unk, cnt_zero])
  329. ## Visualization
  330. ## 월 인덱스 설정 ##
  331. idxCal=[]
  332. idxCalName=[]
  333. idxCal.append(0)
  334. for y_idx in range(lastday.year - startday.year + 1):
  335. if startday.year == lastday.year:
  336. for m_idx in range(lastday.month - startday.month + 1):
  337. month = startday.month + m_idx
  338. idxCal.append(idxCal[-1] + calendar.monthrange(startday.year, month)[1])
  339. idxCalName.append(calendar.month_name[month])
  340. else:
  341. if y_idx == 0: ## 첫번째 해
  342. for m_idx in range(13-startday.month):
  343. month = startday.month + m_idx
  344. idxCal.append(idxCal[-1] + calendar.monthrange(startday.year, month)[1])
  345. idxCalName.append(calendar.month_name[month])
  346. elif y_idx !=0 and y_idx == lastday.year - startday.year: ## 마지막 해
  347. for m_idx in range(lastday.month):
  348. month = m_idx + 1
  349. idxCal.append(idxCal[-1] + calendar.monthrange(lastday.year, month)[1])
  350. idxCalName.append(calendar.month_name[month])
  351. else:
  352. for m_idx in range(12):
  353. month = m_idx + 1
  354. idxCal.append(idxCal[-1] + calendar.monthrange(startday.year+y_idx, month)[1])
  355. idxCalName.append(calendar.month_name[month])
  356. DataCountMat=np.matrix(DataCount)
  357. print("The number of unknown data:",sum(DataCountMat[:,1]), ", The number of zero data:", sum(DataCountMat[:,2]))
  358. plt.figure(figsize=(16,9))
  359. plt.subplot(311)
  360. plt.plot(DataCountMat[:,1],label='Unknown data', linewidth = 2)
  361. plt.plot(DataCountMat[:,2],label='Zero data', linewidth = 2)
  362. # plt.xlabel('Months', fontsize = 16)
  363. plt.ylabel('Data counts', fontsize = 14)
  364. plt.legend(loc='upper left', fontsize = 14)
  365. plt.title("Unknown/zero electric load data per 15min. unit ("+str(startday.year)+"."+str(startday.month)+"."+str(startday.day)+" - "+str(lastday.year)+"."+str(lastday.month)+"."+str(lastday.day)+")", fontsize = 14)
  366. plt.xlim(idxCal[0], idxCal[-1])
  367. plt.xticks(idxCal, idxCalName, fontsize=6.5)
  368. plt.yticks(fontsize=14)
  369. print("Bad data detection complete!")
  370. ### NaN-padding after finding unknown data
  371. ######## 현재 DB 특성상 값이 0으로 찍히거나 시간테이블의 행 자체가 없는 경우가 있고, 이 데이터 1 step 앞뒤로 데이터가 비정상일 확률이 높으므로 비정상데이터 뿐만 아니라 앞뒤 1 step까지 nan으로 처리함
  372. ElectricLoad_Un_ZP=[]
  373. RawDate=[]
  374. idx=0
  375. idx2=0
  376. isBadData = False
  377. for i in range(DataAct_len):
  378. if datetime.date(rawData[idx][4].year,rawData[idx][4].month,rawData[idx][4].day) >= startday and datetime.date(rawData[idx][4].year,rawData[idx][4].month,rawData[idx][4].day) <= lastday:
  379. RawDate.append(StandardTimeStamp[idx2])
  380. if isBadData == True:
  381. ElectricLoad_Un_ZP.append(np.nan)
  382. isBadData=False
  383. elif rawData[idx][5]==0:
  384. ElectricLoad_Un_ZP[-1]=np.nan
  385. ElectricLoad_Un_ZP.append(np.nan)
  386. if rawData[idx+1][5] > 0 and Check_AlivedTimeStamp(rawData, StandardTimeStamp, idx+1, idx2+1):
  387. isBadData = True
  388. elif Check_AlivedTimeStamp(rawData, StandardTimeStamp, idx, idx2):
  389. ElectricLoad_Un_ZP.append(rawData[idx][5])
  390. else:
  391. ElectricLoad_Un_ZP[-1]=np.nan
  392. ElectricLoad_Un_ZP.append(np.nan)
  393. if rawData[idx+1][5] > 0 and Check_AlivedTimeStamp(rawData, StandardTimeStamp, idx+1, idx2+1):
  394. isBadData = True
  395. idx -= 1
  396. idx2 += 1
  397. idx += 1
  398. print('NaN-padding complete!')
  399. # ## Decimation to 1-hour period
  400. ElectricLoad_1h = []
  401. for i in range(DayPeriod):
  402. if i == DayPeriod-1 and isRecent:
  403. Time_len = DataAct_len - i*DataRes_org + 1
  404. else:
  405. Time_len = DataRes_org
  406. isNaN=False
  407. for j in range(Time_len):
  408. if ElectricLoad_Un_ZP[i*4 + j] == np.nan:
  409. isNaN=True
  410. if j%4==3:
  411. if isNaN:
  412. ElectricLoad_1h.append(np.nan)
  413. else:
  414. ElectricLoad_1h.append(sum(ElectricLoad_Un_ZP[i*DataRes_org + j-3:i*DataRes_org + j+1]))
  415. print('Decimation to 1hour complete!')
  416. # ## Data reconstruction using similar-day approach
  417. DateinDay=[]
  418. for k in range(DayPeriod):
  419. DateinDay.append(RawDate[k*DataRes_org])
  420. DoW, DayType = getDayType(DateinDay, DayPeriod, SpecialHoliday)
  421. # Find the similar-day and reconstructed data
  422. marking=np.nan
  423. ReconstructedData, DayType1h = Reconstruction(DayType, ElectricLoad_1h, marking, DataRes_24, isRecent)
  424. plt.subplot(312)
  425. plt.plot(ReconstructedData, '*-', label='Reconstructed data',linewidth=3)
  426. plt.plot(ElectricLoad_1h, '--', label='Raw data',linewidth=3)
  427. plt.legend(loc='upper right', fontsize = 14)
  428. plt.ylabel('Power [kW]', fontsize = 14)
  429. plt.yticks(fontsize=14)
  430. plt.xticks([0],fontsize=14)
  431. plt.xlim((DayPeriod-10)*24, DayPeriod*24)
  432. plt.title('Raw & reconstructed data in the latest 10 days',fontsize=14)
  433. print('Reconstruct complete!')
  434. # ## Day-ahead load forecasting
  435. ####### Convert to matrix
  436. ReconstructedData_Arr=np.zeros((DataRes_24, DayPeriod))
  437. for i in range(DayPeriod):
  438. if isRecent and i==DayPeriod-1:
  439. for j in range(len(ReconstructedData)%DataRes_24):
  440. ReconstructedData_Arr[j,i]=ReconstructedData[i*DataRes_24+j]
  441. else:
  442. for j in range(DataRes_24):
  443. ReconstructedData_Arr[j,i]=ReconstructedData[i*DataRes_24+j]
  444. trn_period=DayPeriod - 1
  445. DayType_m=np.matrix(DayType)
  446. Data_trn=ReconstructedData_Arr[:,0:trn_period]
  447. Data_tst=ReconstructedData_Arr[:,trn_period]
  448. DayType_trn=DayType_m[0:trn_period,:]
  449. DayType_tst=DayType_m[trn_period,:]
  450. cov_lth=np.array([int(linearFilterLength.split(',')[0]),int(linearFilterLength.split(',')[1])])
  451. y_pred_dayAhead = lpc_pred_DayAhead(Data_trn, DayType_trn, cov_lth, DayType_tst, DataRes_24)
  452. print('-------------------------Day-ahead prediction result-------------------------')
  453. if isRecent:
  454. if now.hour == 0:
  455. print('MAPE :', MAPE(Data_tst[0],y_pred_dayAhead[0]), 'MAE :', MAE(Data_tst[0],y_pred_dayAhead[0]))
  456. else:
  457. print('MAPE :', MAPE(Data_tst[0:now.hour],y_pred_dayAhead[0:now.hour]), 'MAE :', MAE(Data_tst[0:now.hour],y_pred_dayAhead[0:now.hour]))
  458. else:
  459. print('MAPE :', MAPE(Data_tst,y_pred_dayAhead),'MAE :', MAE(Data_tst,y_pred_dayAhead))
  460. print('MBE :', MBE(Data_tst,y_pred_dayAhead), 'CVRMSE :', CVRMSE(Data_tst,y_pred_dayAhead))
  461. print('-------------------------------------------------------------------------------')
  462. # ## One-step-ahead load forecasting
  463. y_pred_oneStep=[]
  464. Data_tst_oneStep=[]
  465. if isRecent:
  466. dayHour = now.hour + 1
  467. else:
  468. dayHour = DataRes_24
  469. for i in range(dayHour):
  470. ####### Convert to matrix
  471. ReconstructedData_tmp=ReconstructedData[i:]
  472. if isRecent:
  473. for ii in range(DataRes_24-i):
  474. ReconstructedData_tmp.append(np.nan)
  475. for ii in range(i):
  476. ReconstructedData_tmp.append(np.nan)
  477. ReconstructedData_Arr_oneStep=np.zeros((DataRes_24, DayPeriod))
  478. for j in range(DayPeriod):
  479. for k in range(DataRes_24):
  480. ReconstructedData_Arr_oneStep[k,j]=ReconstructedData_tmp[j*DataRes_24+k]
  481. Data_trn=ReconstructedData_Arr_oneStep[:,0:trn_period]
  482. if isRecent:
  483. Data_tst_oneStep.append(ReconstructedData_Arr_oneStep[i,trn_period])
  484. else:
  485. Data_tst_oneStep=ReconstructedData_Arr[:,trn_period]
  486. y_pred_oneStep.append(lpc_pred_OneStepAhead(Data_trn, DayType_trn, cov_lth, DayType_tst, DataRes_24))
  487. print('-------------------------OneStep-ahead prediction result-------------------------')
  488. if isRecent:
  489. if now.hour == 0:
  490. print('MAPE :', MAPE(Data_tst[0],y_pred_oneStep[0]), 'MAE :', MAE(Data_tst[0],y_pred_oneStep[0]))
  491. else:
  492. print('MAPE :', MAPE(Data_tst[0:now.hour],y_pred_oneStep[0:now.hour]), 'MAE :', MAE(Data_tst[0:now.hour],y_pred_oneStep[0:now.hour]))
  493. else:
  494. print('MAPE :', MAPE(Data_tst_oneStep,y_pred_oneStep),'MAE :', MAE(Data_tst_oneStep,y_pred_oneStep))
  495. print('-------------------------------------------------------------------------------')
  496. plt.subplot(313)
  497. plt.grid(b=True, which='both',axis='y')
  498. if isRecent:
  499. plt.plot(ReconstructedData_Arr[0:now.hour,trn_period], label='Observed data', linewidth=3)
  500. else:
  501. plt.plot(ReconstructedData_Arr[:,trn_period], label='Observed data', linewidth=3)
  502. plt.plot(y_pred_dayAhead, '--', label='Day-ahead Prediction', linewidth=3)
  503. plt.plot(y_pred_oneStep, '*-.', label='OneStep-ahead Prediction', MarkerSize=10, linewidth=3)
  504. plt.xlabel('Time [hour]', fontsize = 14)
  505. plt.ylabel('Power [kW]', fontsize = 14)
  506. plt.legend(loc='upper right', fontsize = 14)
  507. plt.xticks([6,12,18,24],['6','12','18','24'], fontsize = 14)
  508. plt.yticks(fontsize = 14)
  509. plt.ylim(min(ReconstructedData)*0.9,max(ReconstructedData)*1.1)
  510. if isRecent:
  511. plt.title("Electric load forecasting on "+str(now.year)+"."+str(now.month)+"."+str(now.day)+" (Updated every hour) - DGB 2nd branch", fontsize=14)
  512. else:
  513. plt.title("Electric load forecasting on "+str(lastday.year)+"."+str(lastday.month)+"."+str(lastday.day)+" (Updated every hour) - DGB 2nd branch", fontsize=14)
  514. #plt.show()
  515. print("=================== Prediction was successfully finished! ===================")
  516. fig = plt.gcf()
  517. if isRecent:
  518. # Save the figure file of result
  519. # fig.savefig("Result of electric load forecasting on "+str(now.year)+"."+str(now.month)+"."+str(now.day)+" "+str(now.hour)+"h"+str(now.minute)+"m - DGB 2nd branch.png", dpi=fig.dpi)
  520. ### One-hour-ahead load forecasting updated every 1 minute
  521. # MSSQL Access
  522. conn = pymssql.connect(host = targetDBIP, user = targetDBUserID, password = targetDBUserPW, database = targetDBName, autocommit=True)
  523. # Create Cursor from Connection
  524. cursor = conn.cursor()
  525. #########################################################################
  526. #try:
  527. # cursor.execute("INSERT INTO " + targetDBName + ".dbo.BemsMonitoringPointForecastingHourAhead (SiteId,FacilityTypeId,FacilityCode,PropertyId,CreatedDateTime,TargetDateTime,ForecastedValue) VALUES(1,99,4863,1,'" + datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S') + "','"+ datetime.datetime(now.year,now.month,now.day,now.hour,0,0).strftime('%Y-%m-%d %H:00:00') + "',"+str(y_pred_oneStep[-1])+")")
  528. #except:
  529. # print('Hour-ahead forecasted data already exists! (' + (datetime.datetime(now.year,now.month,now.day,now.hour,0,0) - datetime.timedelta(hours=1)).strftime('%Y-%m-%d %H:00:00') + ')')
  530. #########################################################################
  531. ### Day-ahead load forecasting updated every midnight
  532. #if now.hour == 0:
  533. # # Create Cursor from Connection
  534. # cursor = conn.cursor()
  535. # for i in range(len(y_pred_dayAhead)):
  536. # try:
  537. # cursor.execute("INSERT INTO " + targetDBName + ".dbo.BemsMonitoringPointForecastingDayAhead (SiteId,FacilityTypeId,FacilityCode,PropertyId,CreatedDateTime,TargetDateTime,ForecastedValue) VALUES(1,99,4863,1,'" + datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S') + "','" + (datetime.datetime(now.year,now.month,now.day,0,0,0) + datetime.timedelta(hours=i)).strftime('%Y-%m-%d %H:00:00') + "'," + str(y_pred_dayAhead[i]) + ")")
  538. # except:
  539. # print('Day-ahead forecasted data already exists! ('+(datetime.datetime(now.year,now.month,now.day,0,0,0) + datetime.timedelta(hours=i)).strftime('%Y-%m-%d %H:00:00')+')')
  540. #########################################################################
  541. #########################################################################
  542. # Restore the previous data and save data
  543. # Hour-ahead (지우고 다시 쓰기)
  544. currentTime = datetime.datetime(now.year, now.month, now.day, now.hour, 0, 0)
  545. tmpTime = [datetime.datetime(now.year, now.month, now.day, 0, 0, 0)]
  546. while tmpTime[-1] < currentTime:
  547. tmpTime.append(tmpTime[-1] + datetime.timedelta(hours=1))
  548. cursor.execute("DELETE " + targetDBName + ".dbo.BemsMonitoringPointForecastingHourAhead where SiteId = 1 and FacilityTypeId = 99 and FacilityCode = 4863 and PropertyId = 1 and TargetDateTime >= '" + datetime.datetime.now().strftime('%Y-%m-%d 00:00:00') + "' and TargetDateTime < '" + (datetime.datetime.now() + datetime.timedelta(days=1)).strftime('%Y-%m-%d 00:00:00') + "'")
  549. for i in range(len(tmpTime)):
  550. cursor.execute("INSERT INTO " + targetDBName + ".dbo.BemsMonitoringPointForecastingHourAhead (SiteId,FacilityTypeId,FacilityCode,PropertyId,CreatedDateTime,TargetDateTime,ForecastedValue) VALUES(1,99,4863,1,'" + datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S') + "','"+ datetime.datetime(tmpTime[i].year,tmpTime[i].month,tmpTime[i].day,tmpTime[i].hour,0,0).strftime('%Y-%m-%d %H:00:00') + "',"+str(y_pred_oneStep[i])+")")
  551. # Day-ahead (지우고 다시 쓰기)
  552. FinalTime = datetime.datetime(now.year, now.month, now.day, 23, 0, 0)
  553. tmpTime = [datetime.datetime(now.year, now.month, now.day, 0, 0, 0)]
  554. while tmpTime[-1] < FinalTime:
  555. tmpTime.append(tmpTime[-1] + datetime.timedelta(hours=1))
  556. cursor.execute("DELETE " + targetDBName + ".dbo.BemsMonitoringPointForecastingDayAhead where SiteId = 1 and FacilityTypeId = 99 and FacilityCode = 4863 and PropertyId = 1 and TargetDateTime >= '" + datetime.datetime.now().strftime('%Y-%m-%d 00:00:00') + "' and TargetDateTime < '" + (datetime.datetime.now() + datetime.timedelta(days=1)).strftime('%Y-%m-%d 00:00:00') + "'")
  557. for i in range(len(tmpTime)):
  558. cursor.execute("INSERT INTO " + targetDBName + ".dbo.BemsMonitoringPointForecastingDayAhead (SiteId,FacilityTypeId,FacilityCode,PropertyId,CreatedDateTime,TargetDateTime,ForecastedValue) VALUES(1,99,4863,1,'" + datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S') + "','"+ datetime.datetime(tmpTime[i].year,tmpTime[i].month,tmpTime[i].day,tmpTime[i].hour,0,0).strftime('%Y-%m-%d %H:00:00') + "',"+str(y_pred_dayAhead[i])+")")
  559. #########################################################################
  560. conn.close()
  561. print("The result was saved!")
  562. else:
  563. fig.savefig("Result of electric load forecasting on "+str(lastday.year)+"."+str(lastday.month)+"."+str(lastday.day)+ "- DGB 2nd branch.png", dpi=fig.dpi)
  564. plt.show()
  565. print("Sleeping for 60 seconds ...")
  566. else:
  567. print("No data ... Sleeping for 60 seconds ...")
  568. time.sleep(60)