Pythonを用いたノーズカーブの描写方法

# はじめに
ノーズカーブは、負荷側の電圧と有効電力の関係をグラフ化したものであり、人の鼻のような形をしていることからその名になった。ノーズカーブを描写するためには、電力円線図から導出される電圧の二乗からなる4次方程式を解く必要性があるが、今回はそれをPythonのSympyというライブラリを用いてプログラムに代数方程式を解かせることにより描写することができるかを試みた。また、その比較のために媒介変数表示を用いた方法と2つの電圧解を用いた方法について紹介する。

# 送電モデル
まず、送電モデルはいかに示すような中距離送電線モデルを考え、かつ単位法を用いるものとする。
![送電モデル_単位法.JPG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/e3e4b7d5-441c-2a85-e7a4-003faa5f8e73.jpeg)

送電電流は以下のように表すことができる。

```math
\dot{I}=I (cos \theta - j sin \theta)
```
ゆえに、送電端電圧は以下のように表すことができる。
```math
\dot{V_s}=V_r+(r+jx)I(cos\theta -jsin\theta)={V_r+(r cos\theta + x sin\theta)I}+jI(x cos\theta -r sin \theta)

```

ここで、

```math
\begin{equation}
\left\{ \,
    \begin{aligned}
    & P = V_rI cos\theta \\
    & Q = V_rI sin\theta \\

    \end{aligned}
\right.
\end{equation}


```

を用いて、
```math
\dot{V_s}={V_r+(r P/V_r + x Q/V_r)}+j(x P/V_r -r Q/V_r)
```

と表すことができるので、絶対値を取って整理すると以下のような電力円線図の式が得られる。
```math
(V_r^2+Pr+Qx)^2+(xP-rQ)^2=(V_s V_r)^2
```

これを$V_r$について得くことでノーズカーブを描写することができる。
# Sympyによる代数解析による描写方法

```python nose_curve_sympy.py
from sympy import Symbol, symbols, solve, simplify
from sympy.plotting import plot

V_r = Symbol('V_r')
P, Q, x,r,V_s = symbols('P Q x r V_s')
#y = 0となる方程式つまり電力方程式をV＿ｒを軸にして解く
y=(P*r+Q*x+V_r**2)**2+(Q*r-P*x)**2-(V_r*V_s)**2
solved = solve(y,V_r)
solved
#受電端電圧の解全ての表示
print(solved)

V_r_1=solved[1].subs([(Q,0),(r,0),(x,1),(V_s,2),])
V_r_2=solved[3].subs([(Q,0),(r,0),(x,1),(V_s,2),])
#グラフ化、2つのグラフを重ねる
p = plot(V_r_1, V_r_2,(P, 0, 2), legend = True, show = False)

p[0].line_color = 'b'#配色
p[1].line_color = 'r'

p.show()
```

![nose_curve_sympy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/ec547219-7349-8c13-dc6c-cb0518b23219.png)

この方法では、解を求めるために若干時間がかかってしまうのと、グラフのパラメータを微調整しにくいといったデメリットがある。

# 解の表示による描写方法

次に、あらかじめ以下の方程式の解を与えておいてグラフ化させるといった方法もある。

```math
V_r^4+(2(Pr+Qx)-V_s^2)V_r^2+(xP-rQ)=0
```

これの解をプロットする以下のプログラムを作るとノーズカーブをより高速で描写することができる。
```python nose_curve_4pr.py
import numpy as np
import matplotlib.pyplot as plt

P = np.arange(0.0, 1.0, 0.01)


Q=0.0

V_s=1.0

x=1.0
r=0

#y=ax**2+b*x+c(x=V,y=P)

a=1
b=2*(P*r+Q*x)-V_s**2
c=(x*P-r*Q)**2

#安定解
V_r_1=((-b+(b**2-4*a*c)**0.5)/2*a)**0.5
#不安定解
V_r_2=((-b-(b**2-4*a*c)**0.5)/2*a)**0.5




plt.plot(P,V_r_1, 'r-')
plt.plot(P,V_r_2,'b-')
plt.xlim(0,0.5)
plt.ylim(0,1)

plt.xlabel('有効電力P',fontname="MS Gothic")
plt.ylabel('受電端電圧V_r',fontname="MS Gothic")


plt.show()
```
![ノーズカーブ2本.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/3edb7504-2321-fe11-2522-ed2e32456602.png)

この描写方法だと前者と比較して高速に描写することができた。


# 媒介変数表示
最後に、媒介変数表示を用いた描写方法について示す。

受電端電力円線図の複素電力

```math
\dot{S_s}=\dot{E_s}\bar{I}=\frac{E_s e^{j\delta}(E_s e^{-j(\delta - \alpha)}-E_r e^{j \alpha})}{Z}
=\frac{{E_s}^2 e^{j \alpha}-E_s E_r e^{\alpha + \delta}}{Z}=(-\frac{E_r E_s cos(\delta + \alpha)}{Z}+\frac{E_s^2}{Z}cos\alpha)+j(-\frac{E_r E_s sin(\delta + \alpha)}{Z}+\frac{E_s^2}{Z} sin \alpha)

```

と$Q=P tan \theta$を使用することによって、Pと$E_r$を$\delta$の式で表すことを考えると、以下のようなプログラムを書くことができる。
``` python  nose_curve_par.py
import matplotlib.pyplot as plt
import math 
import numpy as np


#パラメータ
E_s=2.0
r=0.1
x=0.5

Z=(r**2+x**2)**0.5

#力率
theta = 20*math.pi/180


alpha = np.arctan(x/r)
#変数項

delta = np.linspace(0,math.pi,100)

E_r = E_s*(np.tan(theta)*np.cos(delta-alpha)-np.sin(delta-alpha))/(np.tan(theta)*np.cos(alpha)-np.sin(alpha))


P=(E_r*E_s/Z)*np.cos(delta-alpha)-(E_r**2/Z)*np.cos(alpha)



plt.xlabel("P")
plt.ylabel("V(V_R)")



plt.plot(P,E_r)
plt.xlim(0,2*E_s) #x軸の範囲
plt.ylim(0,1.01*E_s) #y軸の範囲

plt.show()
```

これを実行すると以下のような画像が出力される。

![ノーズカーブ.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1686019/7cc00fe0-1bf0-72e3-6fa5-fa39a300ddee.png)

媒介変数表示のメリットは、2つに場合分けすることなく描写することができる点にある。

# まとめ
今回はPythonを用いて電力円線図の式からノーズカーブが描写できるかを試みた。結果、代数方程式を解く方法と媒介変数表示を用いる方法を行うことによって上手く描写することができた。しかし、Sympyを用いた代数方程式を解く方法では若干計算に時間がかかることから、計算スピードを考えると、予めある程度代数計算を行っておいた方が良いだろう。

# 参考文献

[電圧安定性に関する解析例](https://www.iee.jp/assets/pes/pdf/ele_systems/ex_model/3.2.1.pdf "電圧安定性に関する解析例")

https://denki-no-shinzui.com/voltage-stability/