Building-Level Energy Transactive [Part I]

Dr.Kongrith Komasatid

Published in

Super AI Engineer

6 min readMar 27, 2021

Part I: Optimal Operation of Battery Energy Storage under Arbitrage Mode

รูปที่ 1 Pilot project สำหรับการพิสูจน์แนวคิด National Energy Trading Platform ของ กฟผ.

1. เกริ่นนำ

นับตั้งแต่ การสาธิตการซื้อขายไฟแบบ Peer-to-Peer (P2P) ในอุตสาหกรรมไฟฟ้าครั้งแรกเมื่อ เม.ย. 2559 ณ เมืองบรูคลิน นิวยอร์ค ประเทศสหรัฐอเมริกา ซึ่งบ้านที่ติดตั้งแผงผลิตไฟฟ้าจากโซลาร์เซลล์ได้ขายพลังงานส่วนเกินไปยังเพื่อนบ้านผ่านสัญญาแบบฉลาด (Smart Contact) ด้วยบล๊อคเชนอีเธอร์เลียม (Ethereum) ซึ่งริเริ่มโดยบริษัทขนาดเล็ก (LO3 Energy) ทำให้มองเห็นภาพในอนาคตที่จะกำลังจะเกิดขึ้นในอุตสาหกรรมไฟฟ้าต่อจากการแปรรูปตลาดซื้อขายไฟฟ้าในอเมริกาและในยุโรป

รูปที่ 2 ระบบไฟฟ้าจำหน่ายปัจจุบันซึ่งมีแนวโน้มการเติบโตของรถไฟฟ้า (EV) และการผลิตไฟฟ้าจากโซลาร์บนหลังคาบ้าน (rooftop solar) source: IRENA

อะไรคือ Energy Arbitrage (การเก็งมูลค่าพลังงาน) ?

การกักเก็บพลังงานในขณะที่มีราคาถูก และนำพลังงานที่เก็บไว้ออกมาขายในขณะที่มีราคาแพง

รูปที่ 3เปรียบเทียบไอซครีมกับการ Enerygy Arbitrage [ref]

แตกต่างจาก Load Leveling (การปรับระดับของโหลด) ?

คล้ายกัน เนื่องจากเป็นการบริหารระดับโหลดช่วงที่มีความต้องการต่ำ (Valley Load) ซึ่งมีต้นทุนไฟฟ้าถูก นำมาตัดยอดการใช้พลังงานไฟฟ้าที่มีความต้องการไฟฟ้าสูง (Peak Load) ซึ่งมีต้นทุนไฟฟ้าราคาแพง

รูปที่ 4การทำงานของระบบกักเก็บพลังงานภายใต้ Typical Load Profile

ใช้หลักการเศรษฐศาสตร์ด้านไหน ?

ใช้หลักการของ Energy Marginal Analysis โดยปัจจัยที่เข้ามาเกี่ยว คือ round-trip efficiency, energy capacity, self-loss และ system marginal cost เป็นต้น

รูปที่ 5 การทำ Load Leveling ตามหลักการของ Energy Marginal Analysis

สามารถติดตาม VDO ปูพื้นฐาน “ Energy-as-a-service (EaaS) และกรณีศึกษาการใช้ Battery Energy Storage System เพื่อทำ Energy Arbitrage [ตอนที่1] ” ทาง Youtube ได้ที่นี่

2. Programming Tutorial

ในหัวข้อนี้ จะเป็นตัวอย่างการนำ Grid-Scale LiFePO4 ขนาด 110kW/110kWh มาทำหน้าที่ Energy Arbitrage สำหรับตึก Eureka ของ กฟผ.

2.1 Installation

ลงโปรแกรมในส่วนที่จำเป็น

!pip install pyunpack
!pip install patool
!pip install PuLP

Import Library ที่จำเป็น

import json
import pandas as pd
import numpy as np
import time
import matplotlib.pyplot as plt
from matplotlib import dates as mpl_dates
from pyunpack import Archive
import seaborn as sns
import datetime as dt
import pulp

2.2 Data Preparation

max_discharge_power_capaacity กับ max_charge_power_capacuty จะเป็นพิกัดกำลังของแบตเตอร์รีในหน่วย kWh efficiency จะเป็น round-trip efficiency ซึ่งไม่มีหน่วย ส่วน discharge_energy_capacity จะเกี่ยวข้องกับความจุของแบตเตอร์รี

max_discharge_power_capacity = 110
max_charge_power_capacity = 110
efficiency = 0.85
discharge_energy_capacity = 200
max_daily_discharged_throughput = 200

ข้อมูลตัวอย่างครั้งนี้เป็นไฟล์ JSON ของเดือน ก.พ. 2564 โดยเราสามารถใช้คำสั่ง pd_read_json ในการอ่านข้อมูลได้ โดยจะมีประเด็นที่น่าสนใจตรง timestamp ซึ่งเป็นเวลามาตรฐาน โดยเราสามารถแปลงมาเป็นเวลาท้องถิ่นได้โดยคำสั่ง df['Date']= df['timestamp'].dt.tz_localize(None).astype('datetime64[s]') จากนั้นก็จัด date format ให้อยู่ในรูปแบบ "%Y-%m-%d %H:%M:%S" เพื่อแสดง ปี-เดือน-วัน ชั่วโมง:นาที:วินาที

รูปที่ 4 แสดงข้อมูลใน JSON ที่เก็บ log ของ API ที่วัดค่าต่างๆ

จากนั้นให้ปัดหน่วย micro second ของ timestamp โดยใช้คำสั่ง round ซึ่งในขั้นตอนนี้เอง อาจจะมีข้อมูลที่เกิดซ้ำซ้อนเนื่องจากการปัดหน่วยเวลาได้ จึงต้องทำการลบข้อมูลซ้ำโดยใช้ drop_duplicatesจากนั้นใช้ timestamp สำหรับเป็น index ของ dataframe

df['Date']=df['Date'].dt.round('1min')
df.drop_duplicates(subset=['Date', 'name'], keep='last', inplace=True)
df.set_index('Date',inplace = True)
df.drop(['timestamp','id','point_ref'], axis=1, inplace=True)
df

จากนั้นให้สกัดข้อมูลที่จำเป็นต้องใช้ในแต่ละ row ออกจาก dataframe โดยจะมีสมการที่สำคัญ ดังนี้

ความต้องการไฟฟ้าสุทธิ = ความต้องการไฟฟ้าที่แท้จริง - (การผลิตไฟฟ้าของโซลาร์เซล + การชาร์จ หรือ ดิสชาร์จของแบตเตอร์รี)

ความต้องการไฟฟ้าที่แท้จริง= ความต้องการไฟฟ้าที่แท้จริง การผลิตไฟฟ้าของโซลาร์เซลล์ (pv_out) ความต้องการไฟฟ้าสุทธิ (P_f1)

รูปที่ 6 แสดงข้อมูลของเดือน ก.พ. 2564 ที่ผ่านการปัดค่า timestamp

จากนั้นให้ทำการ reindex จากข้อมูลรายนาที มาเป็น ข้อมูลรายชั่วโมง

minutely_index = pd.date_range('2021-02-01 00:00:00', '2021-02-28 23:59:00', freq='1min')data1 = P_f1.value.reindex(minutely_index)
data2 = PV1_out.value.reindex(minutely_index)
data3 = Output_HomeBatt1_FLR1.value.reindex(minutely_index)
data4 = Output_HomeBatt2_FLR1.value.reindex(minutely_index)
data5 = Output_HomeBatt3_FLR1.value.reindex(minutely_index)
Load_FLR1 = pd.DataFrame(data1 + data2 + data3 + data4 + data5, index= minutely_index)plt.figure(figsize = (1,5))
sns.heatmap(Load_FLR1.isnull(), cbar=False)

สุดท้ายสามารถหาข้อมูล missing data ได้จากคำสั่ง sns.heatmap และข้อมูลที่ป้อนเข้าให้ใส่ property isnull()

รูปที่ 7 แสดง data missing และ ข้อมูลที่อยู่ระหว่าง sampling period

สำคัญมาก! สามารถทำการ interpolate ข้อมูลที่ null โดยใช้ linear interpolation และหากข้อมูลความต้องการไฟฟ้าเป็น behind-the-meter ให้ทำการสร้างข้อมูล ความต้องการไฟฟ้าที่แท้จริง กลับมาตามสมการที่ให้ข้างบน

ข้อมูลตัวอย่างนี้เป็นตึก 3 ชั้น (Eureka) ซึ่งมี Solar PV ทั้งสามชั้น และมีแบตเตอร์รีแบบ auto smoothing ที่ชั้น 1 และ ชั้น 3 โดยแบตเตอร์รีที่ชั้น 2 จะเป็น Arbitage Mode

## Floor 1 ##
Load_FLR1.interpolate(method='linear',limit_direction='both',inplace=True)
Load_FLR1['value'].clip(lower =0).plot(figsize=(12,6))## Floor 2 ##
data1 = P_f2.value.reindex(minutely_index)
data2 = PV2_out.value.reindex(minutely_index)
data3 = p_out.value.reindex(minutely_index)
Load_FLR2 = pd.DataFrame(data1 + data2 + data3, index= minutely_index)
Load_FLR2.interpolate(method='linear',limit_direction='both',inplace=True)
Load_FLR2['value'].clip(lower =0).plot(figsize=(12,6))## Floor 3 ##
data1 = P_f3.value.reindex(minutely_index)
data2 = PV3_out.value.reindex(minutely_index)
data3 = Output_HomeBatt1_FLR3.value.reindex(minutely_index)
data4 = Output_HomeBatt2_FLR3.value.reindex(minutely_index)
data5 = Output_HomeBatt3_FLR3.value.reindex(minutely_index)TotalLoad = Load_FLR1 + Load_FLR2 + Load_FLR3

2.3 Data Visualization

สามารถ plot กราฟเพื่อดูความต้องการไฟฟ้าสุทธิแยก 3 ชั้นได้ ดังนี้

รูป 8 ข้อมูลความต้องการไฟฟ้าสุทธิ (Net Load) แยกรายชั้น

เขียนโปรแรมแสดงกราฟข้อมูลความต้องการไฟฟ้าที่แท้จริงของทั้งตึก ได้ดังนี้

n = 1
xrange=(1440*7*n, 1440*7*n+1440*7)
yrange=(0,150)
ticks = np.arange(14)*720+1440*7*nlabels = ['','Monday','', 'Tuesday','','Wednesday','','Thursday','','Friday','','Saturday','','Sunday']plt.style.use('ggplot')
plt.figure(figsize=(13,3.5), tight_layout=True)
plt.xticks(ticks, labels)ax = TotalLoad['value'].plot(color='blue', xlim=xrange, ylim=yrange)
ax.set_title('Eureka')
for i in np.arange(7)*1440+1440*7*n:
 ax.axvline(x=i, color='grey',linestyle='--'

สามารถ plot กราฟเพื่อดูข้อมูลความต้องการไฟฟ้าที่แท้จริงของทั้งตึกได้ ดังนี้

รูป 9 ข้อมูลความต้องการไฟฟ้าที่แท้จริงของทั้งตึก (Actual Load)

เขียนโปรแรมแสดง Heatmap ของความต้องการไฟฟ้าสูงสุดของเดือน ก.พ. 2564 ได้ ดังนี้

df_map['Month'] = df_map.index.month
heatmap = pd.pivot_table(data = df_map, index='Month', values='value', columns='day', aggfunc=np.max)
sns.heatmap(heatmap, cmap="YlOrBr", center = (df_map.value.max() - df_map.value.min()) / 2, vmin=0, vmax=140)

จะได้ Heatmap ของความต้องการไฟฟ้าสูงสุดของเดือน ก.พ. 2564 ได้ ดังนี้

รูปที่ 10 Heatmap ของความต้องการไฟฟ้าสูงสุดในแต่ละวันของเดือน ก.พ. 2564

2.4 Model

ทำการสร้าง class ของ Battery โดยมีตัวแปรตัดสินใจเป็น charge power และ discharge power ซึ่งเป็น non-negative real number และกำหนดเป็นปัญหา maximization problem ดังนี้

เพิ่มฟังก์ชันเงื่อนไขใน class ของ Battery :

เงื่อนไข 1) charging constraint ซึ่งเป็นเงื่อนไขที่การชาร์จพลังงานของแบตเตอร์รีที่เวลาใดๆ ต้องมีค่ามากกว่าศูนย์

เงื่อนไข 2) discharging constraint ซึ่งเป็นเงื่อนไขที่การปล่อยพลังงานของแบตเตอร์รีที่เวลาใดๆ ไม่สามารถมากกว่าความจุแบตเตอร์รีได้

เงื่อนไข 3) add_throughput_constraint ซึ่งเป็นเงื่อนไขความจุแบตเตอร์รี

ทำการสร้าง class ของ Tariff โดยอ้างอิง อัตราตามช่วงเวลาของการใช้ (Time of Use Tariff: TOU Tariff) แบบ 2.1.1 ซึ่งหมายถึงกิจการขนาดเล็ก (แรงดัน 12–24 กิโลโวลต์) จากเวปของ กฟน. ดังรูป

รูปที่ 11 อัตราตามช่วงเวลาของการใช้ (Time of Use Tariff: TOU Tariff) ของกิจการขนาดเล็ก

โดยสามารถเขียน class ของ Tariff ได้ ดังนี้

class Tariff():def __init__(self, minutely_index, BusinessType):
 self.index = pd.date_range(start = minutely_index[0], end = minutely_index[-1], freq ='1H')
 self.BusinessType = BusinessType
 if self.BusinessType == 1:
  self.expensive_price = 5.1135
  self.cheap_price = 2.6037
 else:
  self.expensive_price = 5.7982
  self.cheap_price = 2.6369def create_price(self):
 price_df = pd.DataFrame(index = self.index)
 price_df['dow'] = price_df.index.dayofweek
 price_df['hour'] = price_df.index.hour
 price_df['WorkingDay'] = price_df['dow'].apply(lambda x: 1 if x <5 else 0) price_df['WorkingHour'] = price_df['hour'].apply(lambda x: 1 if ((x >= 9) & (x <= 22)) else 0) price_df['OnPeak'] = price_df['WorkingDay'] & price_df['WorkingHour'] price_df['price'] = price_df['OnPeak'].apply(lambda x: self.expensive_price if x == 1 else self.cheap_price)
 price_df.drop(['dow','hour','WorkingDay','WorkingHour','OnPeak'], axis = 1, inplace=True) return price_df

ลองทดสอบการทำงานของ class ของ Tariff

Eureka = Tariff(minutely_index, 1)
price_df = Eureka.create_price()
price_df['price'].plot()
plt.ylabel('Tariff (baht/kWh)')

จะได้กราฟของ Tariff เมื่อพิจารณาช่วง On-Peak และ Off-Peak ได้ ดังรูป

รูปที่ 12 กราฟ Tariff ของกิจการขนาดเล็กภายใต้ On-Peak และ Off-Peak

ทำการสร้าง function simulation สำหรับจำลองการทำ Energy Arbitrage

all_hourly_charges = np.empty(0)
all_hourly_discharges = np.empty(0)
all_hourly_state_of_energy = np.empty(0)
all_daily_discharge_throughput = np.empty(0)battery = Battery(
time_horizon = time_horizon,
max_discharge_power_capacity = max_discharge_power_capacity,
max_charge_power_capacity = max_charge_power_capacity)for day_count in range(28):
 print('Trying day {}'.format(day_count))
 start_time = start_day + pd.Timedelta(day_count, unit='days')
 end_time = start_time + pd.Timedelta(time_horizon-1, unit='hours')
 prices = price_data[start_time:end_time]['price'].valuesbattery.set_objective(prices)
   battery.add_storage_constraints(efficiency=efficiency,min_capacity=0,discharge_energy_capacity=discharge_energy_capacity,
initial_level=initial_level)battery.add_throughput_constraints(max_daily_discharged_throughput=
max_daily_discharged_throughput)battery.solve_model()
 hourly_charges, hourly_discharges = battery.collect_output()daily_discharge_throughput = sum(hourly_discharges)
 net_hourly_activity = (hourly_charges * efficiency) -  hourly_discharges
 cumulative_hourly_activity = np.cumsum(net_hourly_activity)
 state_of_energy_from_t2 = initial_level + cumulative_hourly_activityall_hourly_charges = np.append(all_hourly_charges, hourly_charges)
 all_hourly_discharges = np.append(all_hourly_discharges, hourly_discharges)
 all_hourly_state_of_energy = np.append(all_hourly_state_of_energy, state_of_energy_from_t2)
 all_daily_discharge_throughput=np.append(all_daily_discharge_throughput,daily_discharge_throughput)initial_level = state_of_energy_from_t2[-1]

2.5 Result

ดูการทำงานของ Battery โดยเขียน code ข้างล่าง

plt.rcParams["figure.figsize"] = [5,3]
plt.rcParams["figure.dpi"] = 100
plt.rcParams.update({"font.size":7})
fig1 = plt.figure()
plt.xlabel('kW')
plt.title('Hourly power output')
fig2 = plt.figure()
plt.xlabel('kWh')
plt.title('Hourly state of energy')
ax1 = fig1.add_subplot(111)
ax1.hist(all_hourly_discharges - all_hourly_charges)
ax2 = fig2.add_subplot(111)
ax2.hist(all_hourly_state_of_energy)
plt.show()

เราจะได้ histogram ซ้ายมือแสดง charge/discharge ของแบตเตอร์รีรายชั่วโมง ซึ่ง +100kW, -100kW หมายถึง discharge และ charge ตามลำดับ

histogram ขวามือแสดงพลังงานที่เก็บอยู่ในแบตเตอร์รีรายชั่วโมง โดยมีค่าตั้งแต่ 0 kWh จนถึว 200 kWh (ความจุของแบตเตอร์รี)

รูปที่ 13 กราฟ Histogram ของการทำงานของแบตเตอร์รี และ พลังงานที่สะสมรายชั่วโมง

สรุปการทำงานของ Energy Arbitrage ดังข้างล่าง โปรแกรมแสดงคำตอบ total discharging cost is 20454.00 baht ,total_charging_cost is 10725.85 baht และ total_profit is 9728.15baht

all_data_sim_time = price_df[
pd.Timestamp(year=2021, month=2, day=1, hour=0):
pd.Timestamp(year=2021, month=2, day=28, hour=23)].copy()all_data_sim_time['Charging power'] = all_hourly_charges
all_data_sim_time['Discharging power'] = all_hourly_discharges
all_data_sim_time['Power output'] = all_hourly_discharges - all_hourly_charges
all_data_sim_time['State of Energy'] = np.append(initial_level, all_hourly_state_of_energy[0:-1])
all_data_sim_time['Revenue generation'] = all_data_sim_time['Discharging power'] * all_data_sim_time['price']     all_data_sim_time['Charging cost'] = all_data_sim_time['Charging power'] * all_data_sim_time['price']all_data_sim_time['Profit'] = all_data_sim_time['Revenue generation'] - all_data_sim_time['Charging cost']print('Total discharging cost is %.2f baht' %(all_data_sim_time['Revenue generation'].sum()))print('Total charging cost is %.2f baht' (all_data_sim_time['Charging cost'].sum()))
print('Total profit is %.2f baht \n' %(all_data_sim_time['Profit'].sum()))

ทำการสร้าง function สำหรับการผลิตพลังานของแบตเตอร์รี ได้ ดังนี้

max_profit_week = (all_data_sim_time['Profit'].resample('W').sum() == all_data_sim_time['Profit'].resample('W').sum().max()).values
print(all_data_sim_time['Profit'].resample('W').sum()[max_profit_week])
print('\n')
plt.rcParams["figure.figsize"] = [10,4]
plt.rcParams["figure.dpi"] = 150
plt.rcParams.update({"font.size":7})most_profit_week_start = pd.Timestamp(year=2021, month=2, day=1)
ax = all_data_sim_time[most_profit_week_start:most_profit_week_start+pd.Timedelta(weeks=1)][['State of Energy', 'price']].plot(secondary_y='price', mark_right=False)ax.set_ylabel('State of energy (kWh)')
ax.right_ax.set_ylabel('Tariff (baht/kWh)')
ax.get_legend().set_bbox_to_anchor((0.3, 1))

จะได้กราฟแสดงการ charge/discharge ของวันที่ 1 ก.พ. 2564 ได้ดังนี้

รูปที่ 14 การ charge/discharge ของแบตเตอร์รีช่วงวันที่ 1–8 ก.พ. 2564

ทำการสร้าง function สำหรับแสดง profit รายวันได้ ดังนี้

plt.rcParams["figure.figsize"] = [10,4]
plt.rcParams["figure.dpi"] = 100
plt.rcParams.update({"font.size":8})
all_data_sim_time['Profit'].resample('D').sum().plot()
plt.ylabel('Total Dayly profit (baht)')

จะได้กราฟแสดง profit รายวันของเดือน ก.พ. 2564 ได้ดังนี้

รูปที่ 15 กำไร หรือ ผลประหยัดรายวันของการทำ Energy Arbitrage

3. สรุป

ดังนั้นจึงสามารถสรุปสั้นๆ ได้ว่า

“Energy Arbitrage” เป็นฟังก์ชันพื้นฐานที่ agent ต้องสามารถทำได้สำหรับ HEM (Home Energy Managmenet System)และ BEM (Building Energy Managmenet System) เพราะทำให้สามารถลดค่าใช้จ่ายด้านไฟฟ้าลงได้ จึงเป็นสิ่งจำเป็นสำหรับ EaaS (Energy-as-a-Service) สำหรับตอนต่อไปจะนำเสนอ Deep Q Network สำหรับการทำ Energy Arbitrage ครับ :)

ทิ้งท้าย :

กระผม ดร.คงฤทธิ์ โกมาสถิตย์ เป็นพนักงานรัฐวิสาหกิจที่การไฟฟ้าฝ่ายผลิตแห่งประเทศไทย ตำแหน่งวิศวกรไฟฟ้า โดยบทความนี้เป็นส่วนที่สองของผมสำหรับ โครงการ Super AI Engineer รหัสผู้สมัคร 22p23c0652 หากมีข้อบกพร่องประการผู้เขียนขอน้อมรับและจะนำไปปรับปรุงในอนาคตครับ :)

สุดท้ายนี้ขอขอบคุณแหล่งข้อมูลอ้างอิงที่สำคัญ ดังนี้

http://pierrepinson.com แนะนำให้อ่านสำหรับทำความเข้าใจพื้นฐานเรื่อง Energy Trading
Load Leveling Application of Energy Storage System for Generation Expansion Planning เป็น conference paper ซึ่งตีพิมพ์ที่เกาหลีใต้ ที่แนะนำเกี่ยวกับ Load Leveling
Peer-to-Peer Electricity Trading เป็น technical report ของ IRENA สำหรับปูพื้นฐานด้าน P2P Trading

Building-Level Energy Transactive [Part I]

1. เกริ่นนำ

2. Programming Tutorial

ทิ้งท้าย :

Written by Dr.Kongrith Komasatid