同期発電機のV字曲線

# はじめに
出力の有効電力を一定とした場合、同期発電機の励磁電流と電機子電流をプロットした場合、力率を遅れから進みにするにつれて、V字の曲線を描くことが知られている。これは、数学的なモデルを用いることで証明することができるのだが、今回はPythonを用いて以下のようなV字曲線(遅れを正)を描写することができるかどうかを調査することを目的とする。

![同期発電機のV字曲線_力率角.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/7462e8cf-a54c-7b9d-86c1-8cd988c6ba24.png)

# 導入

V曲線の本質は、同期発電機の起電力（電圧）が励磁電流と比例するということである。

このことを念頭におくと以下のようなベクトル図を考えることができる。（分かりやすいTairaさんのツイッターを引用した）

https://x.com/Taira25619745/status/1584302909713702918

# 同期発電機の諸式

同期インピーダンスのうち抵抗成分は無視して、以下のような等価回路モデルで考える。（１相分）

![同期発電機等価回路_1相分.JPG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/2400118d-918b-8de5-e10e-30eee04416b1.jpeg)

したがって、このような同期発電機モデルには、以下の3式が成立することに気をつける。

**1.受電端電力Pについて**
``` math
P=Real(E_0\bar{I})
```

**2.キルヒホッフの法則について**

``` math
\dot{I}=\frac{\dot{E}-E_0}{jX}
```

**3.励磁電流と電機子起電力の関係について**

励磁電流と電機子起電力には比例関係があり、$k$を比例定数とした場合、以下の式が成立する。

``` math
I_f=\frac{E}{k}
```

**補足**

ただし、前提条件として、発電機側の電圧$\dot{E}$は負荷側（系統側）の電圧と比較して$\delta$分だけ進んでいるとする。

``` math
\dot{E}=E(cos \delta+j sin \delta)
```

という関係があるとする。

また、負荷側の力率に関しては、遅れを正とするので、電流$I$は$E_0$を基準として以下のように表すことができる。

``` math
\dot{I}=I \{cos (-\theta)+j sin (-\theta) \}=I(cos \theta -j sin \theta)
```


## プログラム
### 有効電力分に注目した場合
これによって、Pythonプログラムを作成すると以下のようになる。


``` python doukiki_v.py
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math 

k=1.0
X=1
E_0=1.0
n=1000
#送電端電圧
E=np.linspace(0,2,n)
delta=np.linspace(0,math.pi/2,n)
E,delta=np.meshgrid(E,delta)
E_dot=E*(np.cos(delta)+1j*np.sin(delta))

#励磁電流・電機子電流・受電端電力
I_dot=(E_dot-E_0)/(1j*X)
I_f_dot=E_dot/k
P=(E_0*(I_dot.conjugate())).real
I=abs(I_dot)
I_f=abs(I_f_dot)
plt.xlabel("励磁電流")
plt.ylabel("電機子電流")
plt.contourf(I_f,I,P)
plt.colorbar(label="受電端電力")
plt.savefig("同期発電機のV字曲線_受電端電力.png")
plt.show()
```
これを実行すると以下のようになる。

![同期発電機のV字曲線_受電端電力.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/3bb772b4-6593-2f4c-c65a-fc0ee91dedcc.png)

### 力率に注目した場合

一方で、力率に注目すると以下のようなプログラムを書くことができる。

``` python doukiki_v_t.py
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math 

k=1.0
X=1
E_0=1.0
n=1000
#送電端電圧
E=np.linspace(0,2,n)
delta=np.linspace(0,math.pi/2,n)
E,delta=np.meshgrid(E,delta)
E_dot=E*(np.cos(delta)+1j*np.sin(delta))

#励磁電流・電機子電流・受電端電力
I_dot=(E_dot-E_0)/(1j*X)
I_f_dot=E_dot/k
P=(E_0*(I_dot.conjugate())).real
I=abs(I_dot)
I_f=abs(I_f_dot)

theta=np.angle(I_dot)

pf=np.cos(theta)

plt.xlabel("励磁電流")
plt.ylabel("電機子電流")
plt.contourf(I_f,I,pf)
plt.colorbar(label="力率")
plt.savefig("同期発電機のV字曲線_力率.png")
plt.show()
```

これを実行すると以下のようなV字曲線が描けている事がわかる。
![同期発電機のV字曲線_力率.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/c866b15f-6448-e991-5512-738e4e7b1125.png)

このように、力率が1付近のときに、真ん中に来る。

ただし、この場合だと遅れか進みかを判断できなかったので、力率角に注目してみる。

### 力率角に注目した場合

遅れを正として力率角に注目した場合、以下のようなプログラムを書くことができる。

``` python doukiki_v_t.py
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math 

k=1.0
X=1
E_0=1.0
n=1000
#送電端電圧
E=np.linspace(0,2,n)
delta=np.linspace(0,math.pi/2,n)
E,delta=np.meshgrid(E,delta)
E_dot=E*(np.cos(delta)+1j*np.sin(delta))

#励磁電流・電機子電流・受電端電力
I_dot=(E_dot-E_0)/(1j*X)
I_f_dot=E_dot/k
P=(E_0*(I_dot.conjugate())).real
I=abs(I_dot)
I_f=abs(I_f_dot)


#力率の定義
theta=np.angle(E_0/I_dot)

#弧度法から度数法へ変換
theta = theta*180/math.pi

#pf=np.cos(theta)

plt.xlabel("励磁電流")
plt.ylabel("電機子電流")
#-90<theta<90の範囲を表示する
plt.contourf(I_f,I,theta,np.linspace(-90, 90, 10))
plt.colorbar(label="力率角[度]")
plt.savefig("同期発電機のV字曲線_力率角.png")
plt.show()
``` 

![同期発電機のV字曲線_力率角.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/7462e8cf-a54c-7b9d-86c1-8cd988c6ba24.png)

このように、励磁電流を増やすほど力率角は遅れ方向に上昇する。なので、同期発電機は一般に遅相運転（高励磁運転）を行う。

理由としては、起電力が大きくなるので、電圧の高い運転が行えるため、安定度が向上するといったことなどが挙げられる。


# 代数解析

## 目的
1,2,3式を変形することにより、$I,I_f$は変数$\delta$の関数で表すことを目標とする。これができれば、V曲線を描写することができるはずである。（以下の考察はベクトル図と余弦定理を用いても行うことができる。もう一度Tairaさんの画像を引用するので興味のある方は考えてみて欲しい）

https://x.com/Taira25619745/status/1584302909713702918

## 解析

``` math
P=Real(E_0\bar{I})=Real(E_0 \overline{\frac{E(cos\delta+jsin\delta)-E_0}{jX}})
```

より、発電機出力は

``` math
P=\frac{EE_0}{X}sin\delta 
```

と表されるので、起電力と励磁電流は比例関係にあることから、

``` math
P=\frac{k I_f E_0}{X}sin\delta 
```

したがって、

``` math
I_f=\frac{PX}{kE_0 sin\delta} =f_1(\delta)
```

とすることができる。

一方で、

``` math
\dot{I}=\frac{E(cos\delta + j sin\delta )-E_0}{jX} 
```

より、

``` math
I=\frac{\sqrt{(E cos\delta -E_0)^2 + E^2sin^2 \delta }}{X}=\frac{\sqrt{E^2-2EE_0 cos\delta +E_0^2}}{X}=\frac{\sqrt{(PX)^2-E_0^2 PXsin 2\delta + E_0^4sin^2\delta}}{E_0Xsin\delta}　=f_2(\delta)
``` 

ゆえに、$I,I_f$は変数$\delta$のみに依存する関数であるといえる。言い換えれば、$I,I_f$を$\delta$は媒介変数表示で表すことができる。プログラムでは、この演算を行っている。

# 電験1種との関係
## 同期電動機のV字曲線
最後に、電験一種の二次試験の機械制御で出た問題について少しだけ述べる。（後半は電験についての雑談なので、興味ない方は飛ばして欲しい）

令和1年の問2である。

https://denken-ou.com/c1/kikaiseigyor1-2/

これは、同期電動機のベクトル図に余弦定理を用いることで$I_f,I$の関係を考察している。ただし、この問題では$\delta$を$sin^2 \delta + cos^2 \delta =1$を用いて消去している。複雑な計算に見えるが、原理としては、$E=kI_f$を用いている。プログラムを書くときに、どうして複素数演算を用いたのかというと、このように数式が煩雑になってしまうからである。

このように、電験の問題では計算量が多く本質を突いている問題（重要事項の証明問題）が多い。例を上げれば、電力円線図、電圧降下の近似式、ノーズカーブなどきりがない。一方で、暗記していないと時間内に解けないような問題も存在する。例えば、同期機の過渡応答やニュートン法を用いた潮流解析が挙げられる。これは、重要な背景知識を含んではいるが知っていないとまず解くことができない。

話を戻すと、個人的には令和1年度の二次試験において合否を分けたのは本問ではないかと思う。たしかに、類題は有名な参考書にはあまり載っていないように見える。しかし、本質である起電力と励磁電流の関係性を理解している受験生にとっては格闘することができたのではないだろうか。

電験の計算の世界はある意味で暗記とも言える。便利な公式や重要事項を丸暗記してしまえば、問題を効率良く解くことはできる。また、そうしないと行けない場合も多い。しかし、時には焦らずじっくりと原理原則に立ち返ってみると面白いのではないだろうか。

以下、同期電動機のV字曲線を描写するプログラムを作成した。

``` python doukiki_v_t_moter.py
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math 

k=2.0
X=1
E_0=2.0
n=1000
#送電端電圧
E=np.linspace(0,4,n)
delta=np.linspace(-math.pi/2,0,n)
E,delta=np.meshgrid(E,delta)
E_dot=E*(np.cos(delta)+1j*np.sin(delta))

#励磁電流・電機子電流・受電端電力
I_dot=(E_0-E_dot)/(1j*X)
I_f_dot=E_dot/k
P=(E_dot*(I_dot.conjugate())).real
I=abs(I_dot)
I_f=abs(I_f_dot)


#力率の定義
theta=np.angle(E_0/I_dot)

theta = theta


#弧度法から度数法へ変換
theta = theta*180/math.pi

#pf=np.cos(theta)

plt.xlabel("励磁電流")
plt.ylabel("電機子電流")
#-90<theta<90の範囲を表示する
plt.contourf(I_f,I,theta,np.linspace(-90, 90, 10))
plt.colorbar(label="力率角[度]")
plt.savefig("同期電動機のV字曲線_力率角.png")
plt.show()
``` 

![同期電動機のV字曲線_力率角.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/ad79df98-7b7f-083b-bb1e-4682dcc48da8.png)

ここで、気を付けて欲しいのが、等価回路では発電機と電動機で電流の向きが異なるということである。したがって、相差角の正負も変化するため、力率角の正負が発電機とは逆になる。


# まとめ
今回は同期発電機のV字曲線について等価回路モデルを用いてPythonで描写することができるかを調査した。その結果、見事にV字曲線を描写できた。また、同じV字曲線でも、出力有効電力と力率といった目線で見方を変えると異なった景色が見えた。V字曲線は、電験の試験でも頻出なことから、ぜひとも理解しておこう。