はじめに
電力系統において、負荷が変動した場合は発電機の慣性力によって瞬間的な変動は抑制される。一方で、それよりも期間が長い中程度の擾乱はガバナフリー制御(GF)によって火力入力に対して負のフィードバックをかけることで制御する。また、それよりもさらに長い数十分程度の負荷変動などは負荷周波数制御（LFC）によって調整される。そこで、今回は以下の記事のモデルを参考にし、現代制御の状態方程式を用いることで負荷変動に対する周波数応答についてシミュレーションを行う。具体的にはGFとLFCをそれぞれ使用した場合において周波数がどのように応答するのかを調査する。
GF制御におけるP動作
以下の記事より、GF制御は電力変動に対してP（比例）動作で周波数を制御する。したがって、P動作の宿命上、定常偏差が残ってしまう。（詳細は2番目の記事を参照されたい）
GF制御について
PID制御について
そこで、LFC制御によるI(積分)動作などによって目標周波数と現状周波数の偏差を少なくさせる制御を用いる。
ブロック図
以下のようなブロック図を用いて負荷変動に対する周波数の応答を制御する。
GFとAFCによる周波数の安定化
ただし、このブロック図は以下の記事を参考にした。GFはP動作でLFCはI動作で負のフィードバックを発電式の動揺方程式にかけている。
発電機の動揺方程式
発電機の動揺方程式の簡略系を以下に示す。
M\frac{df}{dt}+Df=P_m-P_e
ただし、$M,D,P_m,P_e$をそれぞれ、慣性定数、ダンピングトルク定数、発電機に係る機械的入力、発電機の電気出力とする。今回は、この簡略式を用いる。
ここで、系統の影響を無視し、負荷の消費電力$P_L=P_e$とする。
なので、系統の周波数$f$を監視し、急な負荷の変動$\Delta P_L$が発生したとき、発電機出力$P_m$を変化させることで過度の系統周波数の変動を抑制する。
プログラム
さて、上記のブロック図を状態方程式に直して数値計算することを考える。
以下のPID制御に関する記事のプログラムを改造する。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
#発電機のパラメータ
#慣性定数
M=1
#ダンピングトルク
D=2
t_max=30
#発電機の状態方程式
A=np.array([[1,1],[1,1]])*(D/M)
b=np.array([1,1])*(1/M)
c=np.array([1,0])
delta_f_opt=0
z=0
delta_t=1.0*10**-5
def control_main(K_p,K_i):
global x,v,z
#初期条件
t=-1
x_1=0
x_2=0
v_1=0
v_2=0
v=np.array([v_1,v_2])
x=np.array([x_1,x_2])
#変数セット
x_1_ary=[]
x_2_ary=[]
v_ary=[]
t_ary=[]
z_ary=[]
delta_f_ary=[]
delta_P_L_ary=[]
while t<t_max:
if t>0:
delta_P_L=-0.1
else:
delta_P_L=0
delta_f=c@x
u=-delta_P_L+K_i*z+K_p*(delta_f_opt-delta_f)
v= A@x+b*u
x=x+v*delta_t
z=z+(delta_f_opt-delta_f)*delta_t
v_ary.append(v)
x_1_ary.append(x[0])
x_2_ary.append(x[1])
t_ary.append(t)
delta_f_ary.append(delta_f)
delta_P_L_ary.append(delta_P_L)
z_ary.append(z)
t=t+delta_t
plt.title("gf="+str(K_p)+",LFC="+str(K_i))
plt.plot(t_ary,delta_f_ary,color="blue",label="周波数変動Δｆ分")
plt.plot(t_ary,delta_P_L_ary,color="red",label="負荷電力変動分")
plt.legend()
#GFの定数
K_p=10
#LFCの定数
K_i=0
control_main(K_p,K_i)
plt.savefig("gf="+str(K_p)+",LFC="+str(K_i)+"_3.png")
plt.show()
結果
負荷が重くなると周波数が低下する。一方で、負荷が軽くなると周波数は増加する。また、GFのみだと定状偏差（オフセット）が残ってしまう。それを防ぐために、I動作を保持しているLFC制御もかける。
負荷が軽くなった場合
GFのみ
周波数が増加し、そのまま新たなる平衡点に行ってしまった。
GF+LFC
周波数が一瞬増加したが、もとの周波数に回復し安定した。
負荷が重くなった場合
GFのみ
GF+LFC
まとめ
今回は、GF,LFC制御による周波数制御について状態方程式を用いることでシミュレーションを行った。結果、PID制御と同等に周波数変動を抑制させることができた。しかし、GF制御のみだとP動作しか働かないので、オフセットが残ってしまった。そこで、LFC制御によるI動作をかけることでオフセットを無くす方向に制御をかけることができた。ただし、今回の扱ったモデルはかなり単純化したモデルである。
参考
さて、負荷変動$\Delta P_L$が発生した場合は、系統周波数変動$\Delta f$を抑制するために、発電機出力$P_m$を変化させる。（負荷側を切り離しという手もあるが、今回は考えないものとする。）
ここで、動揺方程式の状態方程式の入力部$u$が、発電機に係る機械的入力$\Delta P_m$となる。したがって、以下のようなプログラムを作成した。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
#発電機のパラメータ
#慣性定数
M=1
#ダンピングトルク
D=2
t_max=30
#発電機の状態方程式
A=np.array([[1,1],[1,1]])*(D/M)
b=np.array([1,1])*(1/M)
c=np.array([1,0])
delta_f_opt=0
z=0
delta_t=1.0*10**-5
def control_main(K_p,K_i):
global x,v,z
#初期条件
t=-1
x_1=0
x_2=0
v_1=0
v_2=0
v=np.array([v_1,v_2])
x=np.array([x_1,x_2])
#変数セット
x_1_ary=[]
x_2_ary=[]
v_ary=[]
t_ary=[]
z_ary=[]
delta_f_ary=[]
delta_P_L_ary=[]
delta_P_m_ary=[]
while t<t_max:
if t>0:
delta_P_L=-0.1
else:
delta_P_L=0
delta_f=c@x
u=-delta_P_L+K_i*z+K_p*(delta_f_opt-delta_f)
v= A@x+b*u
x=x+v*delta_t
z=z+(delta_f_opt-delta_f)*delta_t
delta_P_m_ary.append(u)
v_ary.append(v)
x_1_ary.append(x[0])
x_2_ary.append(x[1])
t_ary.append(t)
delta_f_ary.append(delta_f)
delta_P_L_ary.append(delta_P_L)
z_ary.append(z)
t=t+delta_t
plt.title("gf="+str(K_p)+",LFC="+str(K_i))
plt.plot(t_ary,delta_f_ary,color="blue",label="周波数変動Δｆ分")
plt.plot(t_ary,delta_P_L_ary,color="red",label="負荷電力変動分")
plt.plot(t_ary,delta_P_m_ary,color="green",label="発電機出力変動分")
plt.legend()
#GFの定数
K_p=10
#LFCの定数
K_i=1
control_main(K_p,K_i)
plt.savefig("gf="+str(K_p)+",LFC="+str(K_i)+"_3.png")
plt.show()
これを実行すると以下のようになった。
ガバナフリー運転のみだと、負荷の脱落変動が起きた瞬間に発電機の出力を急激に絞る。そうすることで、周波数の上昇自体は抑えることができるが、最初の周波数の安定状態よりも高い周波数で平衡してしまった。
一方でガバナフリー運転とLFC運転を合わせた場合、以下のようになった。
負荷の脱落変動とともに、周波数は上昇するが、発電機出力を緩やかに0へ絞ることで、周波数を定常状態に戻している。