本記事では、数値解析法の1つである「ニュートン・ラフソン法」を用いた潮流計算の手法について、計算例を通じて解説する。
本記事では、数値解析法の1つである「ニュートン・ラフソン法」を用いた潮流計算の手法について解説する。ニュートン・ラフソン法の概要ニュートン・ラフソン法（Newton-Raphson method）は数値解析法の1つで、非線形[…]
計算準備
適用系統と計算条件
ニュートン・ラフソン法を潮流計算に用いるにあたり、適用例として電験一種二次試験「電力・管理」令和3年度問3にて出題された、図1に示すような3母線の電力系統を扱う。
図1 3母線の電力系統
図1の系統における計算条件は、下記とする（いずれの値も単位は$[\mathrm{p.u.}]$とする）。
- 母線1は基準ノード（スラックノード）で、電圧$\dot{V}_1=1.0+j0=1.0\angle0$である。
- 母線2は$P-V$指定ノードで、有効電力$P_2=0.7$，電圧の大きさ$V_2=\left|\dot{V}_2\right|=1.05$である。
- 母線3は$P-Q$指定ノードで、有効電力$P_3=-1.0$，無効電力$Q_3=-0.2$である。
- ノードアドミタンス行列$\boldsymbol{\dot{Y}}$は、下記である。
\boldsymbol{\dot{Y}}&=\displaystyle\left(\begin{array}{ccc} \dot{Y}_{11} & \dot{Y}_{12} & \dot{Y}_{13} \\\dot{Y}_{21} & \dot{Y}_{22} & \dot{Y}_{23}\\ \dot{Y}_{31} & \dot{Y}_{32} & \dot{Y}_{33}\end{array}\right)\\\\
&=\displaystyle\left(\begin{array}{ccc} -j18 & j10 & j10 \\ j10 & -j18 & j10 \\ j10 & j10 & -j18\end{array}\right)
\end{align*}$$
今回の計算の目的として、上記の条件を満たすような母線電圧$\dot{V}_2,\ \dot{V}_3$の値（大きさおよび位相）を求めていくことになる。
ニュートン・ラフソン法の計算フロー
ニュートン・ラフソン法では、後に示す修正方程式について、得られた近似解が希望の値になるまで繰り返し解いていくことになる。
この一連の計算の流れを表したフローチャートを図2に示す。
図2 ニュートン・ラフソン法の計算フロー
なお今回、収束条件としては誤差の値が$\varepsilon=1\times10^{-7}$以下となった場合に計算を終了するものとする。
図2のカッコ付き添字「$^\left(i\right)$」については文献によって若干異なるものの、本記事では収束計算の「番号」を示すこととする（単純に回数ではない）。
例えば本記事では1回目の収束計算を示すが、この場合は計算フロー上の最初の計算であり、$i=0$となる。
ニュートン・ラフソン法の計算（直交座標表示）
⓪初期値の設定
まず、母線電圧について、母線1は前項の計算条件より、
$$\dot{V}_1=e_1+jf_1=1.0+j0$$
となる。
次に、母線2および3の電圧$\dot{V}_2,\ \dot{V}_3$を、変数$e_2,\ f_2,\ e_3,\ f_3$を用いて、
$$\begin{cases}
\dot{V}_2=e_2+jf_2\\\\
\dot{V}_3=e_3+jf_3
\end{cases}$$
これらの変数の初期値として、今回は次のようにおく。
①指定値に対する誤差の計算
ここで、与えられたノードアドミタンス行列は、次のように書き換えることができる。
\boldsymbol{\dot{Y}}&=\displaystyle\left(\begin{array}{ccc} G_{11}+jB_{11} & G_{12}+jB_{12} & G_{13}+jB_{13} \\ G_{21}+jB_{21} & G_{22}+jB_{22} & G_{23}+jB_{23}\\ G_{31}+jB_{31} & G_{32}+jB_{32} & G_{33}+jB_{33}\end{array}\right)\\\\
&=\displaystyle\left(\begin{array}{ccc} 0-j18 & 0+j10 & 0+j10 \\ 0+j10 & 0-j18 & 0+j10 \\ 0+j10 & 0+j10 & 0-j18\end{array}\right)
\end{align*}$$
上記の値を用いると、母線2における電力方程式は、
P_2&=\displaystyle \sum_{l=1}^3\left\{G_{2l}\left(e_2e_l+f_2f_l\right)-B_{2l}\left(e_2f_l-e_lf_2\right)\right\}\\\\
&=G_{21}\left(e_2e_1+f_2f_1\right)-B_{21}\left(e_2f_1-e_1f_2\right)+G_{22}\left(e^2_2+f^2_2\right)-B_{22}\left(e_2f_2-e_2f_2\right)\\\\
&\qquad+G_{23}\left(e_2e_3+f_2f_3\right)-B_{23}\left(e_2f_3-e_3f_2\right)\\\\
&=10f_2-10e_2f_3+10e_3f_2 ・・・(1)\\\\\\
V_2&=e^2_2+f^2_2 ・・・(2)
\end{align*}$$
同様に、母線3における電力方程式は、
P_3&=\displaystyle \sum_{l=1}^3\left\{G_{3l}\left(e_3e_l+f_3f_l\right)-B_{3l}\left(e_3f_l-e_lf_3\right)\right\}\\\\
&=G_{31}\left(e_3e_1+f_3f_1\right)-B_{31}\left(e_3f_1-e_1f_3\right)+G_{32}\left(e_3e_2+f_2f_3\right)-B_{32}\left(e_3f_2-e_2f_3\right)\\\\
&\qquad+G_{33}\left(e^2_3+f^2_3\right)-B_{33}\left(e_3f_3-e_3f_3\right)\\\\
&=10f_3-10e_3f_2+10e_2f_3 ・・・(3)\\\\\\
Q_3&=\displaystyle \sum_{l=1}^3\left\{G_{3l}\left(e_lf_3-e_3f_l\right)-B_{3l}\left(e_3e_l+f_3f_l\right)\right\}\\\\
&=G_{31}\left(e_1f_3-e_3f_1\right)-B_{31}\left(e_3e_1+f_3f_1\right)+G_{32}\left(e_2f_3-e_3f_2\right)-B_{32}\left(e_3e_2+f_3f_2\right)\\\\
&\qquad+G_{33}\left(e_3f_3-e_3f_3\right)-B_{33}\left(e^2_3+f^2_3\right)\\\\
&=-10e_3-10e_2e_3-10f_2f_3+18e^2_3+18f^2_3 ・・・(4)
\end{align*}$$
本記事では、潮流計算を行うために必要となる電力方程式について解説する。潮流計算を行うための要素潮流計算（Load flow（またはPower flow）analysis）は、電力系統における発電機の発電出力や負荷の消費電力が[…]
このとき、先程設定した初期値$e^{\left(0\right)}_2,\ f^{\left(0\right)}_2,\ e^{\left(0\right)}_3,\ f^{\left(0\right)}_3$における計算値$P^{\left(0\right)}_2,\ V^{\left(0\right)2}_2,\ P^{\left(0\right)}_3,\ Q^{\left(0\right)}_3$は、$(1)\sim(4)$式より、
P^{\left(0\right)}_2&=10f^{\left(0\right)}_2-10e^{\left(0\right)}_2f^{\left(0\right)}_3+10e^{\left(0\right)}_3f^{\left(0\right)}_2=0\\\\
V^{\left(0\right)2}_2&=e^{\left(0\right)2}_2+f^{\left(0\right)2}_2=1.0\\\\
P^{\left(0\right)}_3&=10f^{\left(0\right)}_3-10e^{\left(0\right)}_3f^{\left(0\right)}_2+10e^{\left(0\right)}_2f^{\left(0\right)}_3=0\\\\
Q^{\left(0\right)}_3&=-10e^{\left(0\right)}_3-10e^{\left(0\right)}_2e^{\left(0\right)}_3-10f^{\left(0\right)}_2f^{\left(0\right)}_3+18e^{\left(0\right)2}_3+18f^{\left(0\right)2}_3\\\\
&=-10-10+18=-2.0
\end{align*}$$
また、指定値については、計算条件より、
$$\begin{cases}
P_{2s}&=0.7\\\\
V^2_{2s}&=1.05^2\\\\
P_{3s}&=-1.0\\\\
Q_{3s}&=-0.2
\end{cases}$$
したがって、指定値と計算値との誤差$\Delta P^{\left(0\right)}_2,\ \Delta V^{\left(0\right)2}_2,\ \Delta P^{\left(0\right)}_3,\ \Delta Q^{\left(0\right)}_3$は、
\Delta P^{\left(0\right)}_2&=P_{2s}-P^{\left(0\right)}_2=0.7-0=0.7\\\\
\Delta V^{\left(0\right)2}_2&=V^2_{2s}-V^{\left(0\right)2}_2=1.05^2-1.0=0.1025\\\\
\Delta P^{\left(0\right)}_3&=P_{3s}-P^{\left(0\right)}_3=-1.0-0=-1.0\\\\
\Delta Q^{\left(0\right)}_3&=Q_{3s}-Q^{\left(0\right)}_3=-0.2-\left(-2.0\right)=1.8
\end{align*}$$
②ヤコビ行列の計算
ニュートン・ラフソン法において、1回目（$i=0$）の収束計算における修正方程式は次のようになる。
\left(\begin{array}{c} \Delta P^\left(0\right)_2 \\ \Delta V^{\left(0\right)2}_2 \\ \Delta P^\left(0\right)_3 \\ \Delta Q^\left(0\right)_3 \end{array}\right)&=\left(\begin{array}{cccc}\displaystyle\frac{\partial P_2}{\partial e_2} & \displaystyle\frac{\partial P_2}{\partial f_2} & \displaystyle\frac{\partial P_2}{\partial e_3} & \displaystyle\frac{\partial P_2}{\partial f_3} \\ \displaystyle\frac{\partial V^2_2}{\partial e_2} & \displaystyle\frac{\partial V^2_2}{\partial f_2} & \displaystyle\frac{\partial V^2_2}{\partial e_3} & \displaystyle\frac{\partial V^2_2}{\partial f_3} \\ \displaystyle\frac{\partial P_3}{\partial e_2} & \displaystyle\frac{\partial P_3}{\partial f_2} & \displaystyle\frac{\partial P_3}{\partial e_3} & \displaystyle\frac{\partial P_3}{\partial f_3} \\ \displaystyle\frac{\partial Q_3}{\partial e_2} & \displaystyle\frac{\partial Q_3}{\partial f_2} & \displaystyle\frac{\partial Q_3}{\partial e_3} & \displaystyle\frac{\partial Q_3}{\partial f_3} \end{array}\right)\left(\begin{array}{c} \Delta e^\left(1\right)_2 \\ \Delta f^\left(1\right)_2 \\ \Delta e^\left(1\right)_3 \\ \Delta f^\left(1\right)_3\end{array}\right)\\\\
&=\boldsymbol{J}^{\left(0\right)}\left(\begin{array}{c} \Delta e^\left(1\right)_2 \\ \Delta f^\left(1\right)_2 \\ \Delta e^\left(1\right)_3 \\ \Delta f^\left(1\right)_3\end{array}\right) ・・・(5)\\\\
\left(\begin{array}{c} e^\left(1\right)_2 \\ f^\left(1\right)_2 \\ e^\left(1\right)_3 \\ f^\left(1\right)_3\end{array}\right)&=\left(\begin{array}{c} e^\left(0\right)_2 \\ f^\left(0\right)_2 \\ e^\left(0\right)_3 \\ f^\left(0\right)_3\end{array}\right)+\left(\begin{array}{c} \Delta e^\left(1\right)_2 \\ \Delta f^\left(1\right)_2 \\ \Delta e^\left(1\right)_3 \\ \Delta f^\left(1\right)_3\end{array}\right) ・・・(6)
\end{align*}$$
なお、$(5)$式のヤコビ行列$\boldsymbol{J}^{\left(0\right)}$の係数はすべて$e^{\left(0\right)}_2,\ f^{\left(0\right)}_2,\ e^{\left(0\right)}_3,\ f^{\left(0\right)}_3$における値とする。
ここで、修正方程式の係数を計算する式は、$(1)\sim(4)$式をそれぞれの変数で偏微分して、
\frac{\partial P_2}{\partial e_2}&=-10f_3,\ \frac{\partial P_2}{\partial f_2}=10+10e_3\\\\
\frac{\partial P_2}{\partial e_3}&=10f_2,\ \frac{\partial P_2}{\partial f_3}=-10e_3\\\\
\frac{\partial V^2_2}{\partial e_2}&=2e_2,\ \frac{\partial V^2_2}{\partial f_2}=2f_2\\\\
\frac{\partial V^2_2}{\partial e_3}&=0,\ \frac{\partial V^2_2}{\partial f_3}=0\\\\
\frac{\partial P_3}{\partial e_2}&=10f_3,\ \frac{\partial P_3}{\partial f_2}=-10e_3\\\\
\frac{\partial P_3}{\partial e_3}&=-10f_2,\ \frac{\partial P_3}{\partial f_3}=10+10e_2\\\\
\frac{\partial Q_3}{\partial e_2}&=-10e_3,\ \frac{\partial Q_3}{\partial f_2}=-10f_3\\\\
\frac{\partial Q_3}{\partial e_3}&=-10-10e_2+36e_3,\ \frac{\partial Q_3}{\partial f_3}=-10+36f_3
\end{align} ・・・(7)$$
$(7)$式より、初期値$e^{\left(0\right)}_2,\ f^{\left(0\right)}_2,\ e^{\left(0\right)}_3,\ f^{\left(0\right)}_3$におけるヤコビ行列$\boldsymbol{J}^{\left(0\right)}$は、
上記の結果より、逆行列$\left[\boldsymbol{J}^{\left(0\right)}\right]^{-1}$を計算すると、
今回、逆行列の計算は、後に示すExcelの計算表によって行った。
逆行列を求めるのはこのニュートン・ラフソン法の計算の中では特に煩雑になり、かつ収束するまで何度も計算を行わなければならない。
手計算で行うにしても（そのような事態があるかは分からないが）、逆行列を含む行列の計算はこちらのようなサイトを利用することをお勧めする。
③収束判定と近似解の設定
求めた逆行列$\left[\boldsymbol{J}^{\left(0\right)}\right]^{-1}$を用いて、$(5)$式の修正方程式により誤差$\Delta e^\left(1\right)_2,\ \Delta f^\left(1\right)_2,\ \Delta e^\left(1\right)_3,\ \Delta f^\left(1\right)_3$を計算すると、
\left(\begin{array}{c} \Delta e^\left(1\right)_2 \\ \Delta f^\left(1\right)_2 \\ \Delta e^\left(1\right)_3 \\ \Delta f^\left(1\right)_3\end{array}\right)&=\left[\boldsymbol{J}^{\left(0\right)}\right]^{-1}\left(\begin{array}{c} \Delta P^\left(1\right)_2 \\ \Delta V^{\left(1\right)2}_2 \\ \Delta P^\left(1\right)_3 \\ \Delta Q^\left(1\right)_3 \end{array}\right)\\\\
&=\left(\begin{array}{cccc}0 & 0.500000 & 0 & 0 \\ 0.066667 & 0 & 0.033333 & 0 \\ 0.020833 & 0.312500 & 0.041667 & 0.062500 \\ 0.033333 & 0 & 0.066667 & 0 \end{array}\right)\left(\begin{array}{c} 0.7 \\ 0.1025 \\ -1.0 \\ 1.8 \end{array}\right)\\\\
&=\left(\begin{array}{c} 0.512500 \\ 0.013333 \\ 0.117448 \\ -0.043333 \end{array}\right)
\end{align*}$$
上記より、誤差の大きさは最大で$\left|\Delta e^\left(1\right)_2\right|=0.512500>1\times10^{-7}$となることから、計算フロー上の収束判定はNGとなり、また①のプロセスに戻って計算を続行する。
また、1回目の収束計算における近似解$e^\left(1\right)_2,\ f^\left(1\right)_2,\ e^\left(1\right)_3,\ f^\left(1\right)_3$は、$(6)$式より、
\left(\begin{array}{c} e^\left(1\right)_2 \\ f^\left(1\right)_2 \\ e^\left(1\right)_3 \\ f^\left(1\right)_3\end{array}\right)&=\left(\begin{array}{c} e^\left(0\right)_2 \\ f^\left(0\right)_2 \\ e^\left(0\right)_3 \\ f^\left(0\right)_3\end{array}\right)+\left(\begin{array}{c} \Delta e^\left(1\right)_2 \\ \Delta f^\left(1\right)_2 \\ \Delta e^\left(1\right)_3 \\ \Delta f^\left(1\right)_3\end{array}\right)\\\\
&=\left(\begin{array}{c} 1.0 \\ 0 \\ 1.0 \\ 0 \end{array}\right)+\left(\begin{array}{c} 0.512500 \\ 0.013333 \\ 0.117448 \\ -0.043333 \end{array}\right)\\\\
&=\left(\begin{array}{c} 1.512500 \\ 0.013333 \\ 1.117448 \\ -0.043333 \end{array}\right)
\end{align*}$$
計算結果
上記の①〜③の計算プロセスを、誤差の値が$\varepsilon=1\times10^{-7}$以下になるまで繰り返し、次の結果が得られた。
今回は5回目の計算（$i=4$）で収束する（誤差の最大値は$\left|\Delta e_3\right|=4.242\times10^{-9}<\varepsilon$）という結果となった。
先述のように、今回の収束計算はExcelで行った。
Excelでは行列の計算も行うことができ、例えば1回目の収束計算において、下記のような入力をおこなっている。
- 逆行列：セルに「=MINVERSE(範囲)」を入力
- 行列の積：セルに「=MMULT(行列1の範囲,行列2の範囲)」を入力
上記より、得られた解$e_2,\ f_2,\ e_3,\ f_3$を元に算出した電圧$\dot{V}_2,\ \dot{V}_3$の値は、次のようになる。
$$\begin{cases}
\dot{V}_2&=1.04993+j0.012119\\\\
\dot{V}_3&=1.12717-j0.042119
\end{cases}$$
なお、上記の解および電力方程式により、母線1における電力$P_1,\ Q_1$を計算すると、
P_1&=\displaystyle \sum_{l=1}^3\left\{G_{1l}\left(e_1e_l+f_1f_l\right)-B_{1l}\left(e_1f_l-e_lf_1\right)\right\}\\\\
&=G_{11}\left(e^2_1+f^2_1\right)-B_{11}\left(e_1f_1-e_1f_1\right)+G_{12}\left(e_1e_2+f_1f_2\right)-B_{12}\left(e_1f_2-e_2f_1\right)\\\\
&\qquad+G_{13}\left(e_1e_3+f_1f_3\right)-B_{13}\left(e_1f_3-e_3f_1\right)\\\\
&=-10f_2-10f_3\\\\
&=-10\times0.012119-10\times\left(-0.042119\right)\\\\
&=0.3\\\\\\
Q_1&=\displaystyle \sum_{l=1}^3\left\{G_{1l}\left(e_lf_1-e_1f_l\right)-B_{1l}\left(e_1e_l+f_1f_l\right)\right\}\\\\
&=G_{11}\left(e_1f_1-e_1f_1\right)-B_{11}\left(e^2_1+f^2_1\right)+G_{12}\left(e_2f_1-e_1f_2\right)-B_{12}\left(e_1e_2+f_1f_2\right)\\\\
&\qquad+G_{13}\left(e_3f_1-e_1f_3\right)-B_{13}\left(e_1e_3+f_1f_3\right)\\\\
&=18-10e_2-10e_3\\\\
&=18-10\times1.04993-10\times1.12717\\\\
&=-3.77097
\end{align*}$$
上記の結果より、今回は図1の系統にコンダクタンス成分がなく、電力損失がないため$P_1+P_2+P_3=0.3+0.7-1.0=0$となることがわかる。
また、今回は母線1の電圧が$1.0\mathrm{p.u.}$である一方で母線2および3の電圧ともに高い値であり、このような状況にするためには母線1から大きな進相無効電力（∵符号が負）を注入しなければならないことがわかる。
ニュートン・ラフソン法の計算（極座標表示）
⓪初期値の設定
次に極座標表示においても、直交座標表示と同様のプロセスにおいて計算を行っていく。
まず、母線電圧について、母線1は前項の計算条件より、
$$\dot{V}_1=V_1\angle\delta_1=1.0\angle0$$
となる。
次に、母線2は$P-V$指定ノードであり、電圧$\dot{V}_2$の大きさ$V_2=\left|\dot{V}_2\right|$は与えられている。
したがって、今回の計算における変数は母線2の電圧位相$\delta_2$，および母線3の電圧の大きさ$V_3=\left|\dot{V}_3\right|$および位相$\delta_3$となる。
これらの変数の初期値として、今回は次のようにおく。
①指定値に対する誤差の計算
ここで、与えられたノードアドミタンス行列は、次のように書き換えることができる。
\boldsymbol{\dot{Y}}&=\displaystyle\left(\begin{array}{ccc} Y_{11}\angle\theta_{11} & Y_{12}\angle\theta_{12} & Y_{13}\angle\theta_{13} \\ Y_{21}\angle\theta_{21} & Y_{22}\angle\theta_{22} & Y_{23}\angle\theta_{23}\\ Y_{31}\angle\theta_{31} & Y_{32}\angle\theta_{32} & Y_{33}\angle\theta_{33}\end{array}\right)\\\\
&=\displaystyle\left(\begin{array}{ccc} 18\angle-\displaystyle\frac{\pi}{2} & 10\angle\displaystyle\frac{\pi}{2} & 10\angle\displaystyle\frac{\pi}{2} \\ 10\angle\displaystyle\frac{\pi}{2} & 18\angle-\displaystyle\frac{\pi}{2} & 10\angle\displaystyle\frac{\pi}{2} \\ 10\angle\displaystyle\frac{\pi}{2} & 10\angle\displaystyle\frac{\pi}{2} & 18\angle-\displaystyle\frac{\pi}{2} \end{array}\right)
\end{align*}$$
上記の値を用いるとともに、$\delta_{kl}=\delta_{k}-\delta_{l}$であることに注意すると、母線2における電力方程式は、
P_2&=\displaystyle \sum_{l=1}^3Y_{2l}V_lV_2\cos\left(\delta_{2l}-\theta_{2l}\right)\\\\
&=Y_{21}V_1V_2\cos\left(\delta_{21}-\theta_{21}\right)+Y_{22}V^2_2\cos\left(\delta_{22}-\theta_{22}\right)\\\\
&\qquad+Y_{23}V_1V_3\cos\left(\delta_{23}-\theta_{23}\right)\\\\
&=Y_{21}V_1V_2\cos\left(\delta_{21}-\frac{\pi}{2}\right)+Y_{22}V^2_2\cos\left(\delta_{22}+\frac{\pi}{2}\right)\\\\
&\qquad+Y_{23}V_1V_3\cos\left(\delta_{23}-\frac{\pi}{2}\right)\\\\
&=Y_{21}V_1V_2\sin\delta_{21}-Y_{22}V^2_2\sin\delta_{22}+Y_{23}V_1V_3\sin\delta_{23}\\\\
&=10\cdot1.0\cdot1.05\sin\left(\delta_2-\delta_1\right)-18\cdot1.05^2\sin\left(\delta_2-\delta_2\right)\\\\
&\qquad+10\cdot1.05V_3\sin\left(\delta_2-\delta_3\right)\\\\
&=10.5\sin\delta_2+10.5V_3\sin\left(\delta_2-\delta_3\right) ・・・(8)
\end{align*}$$
同様に、母線3における電力方程式は、
P_3&=\displaystyle \sum_{l=1}^3Y_{3l}V_lV_3\cos\left(\delta_{3l}-\theta_{3l}\right)\\\\
&=Y_{31}V_1V_3\cos\left(\delta_{31}-\theta_{31}\right)+Y_{32}V_2V_3\cos\left(\delta_{32}-\theta_{32}\right)\\\\
&\qquad+Y_{33}V^2_3\cos\left(\delta_{33}-\theta_{33}\right)\\\\
&=Y_{31}V_1V_3\cos\left(\delta_{31}-\frac{\pi}{2}\right)+Y_{32}V_2V_3\cos\left(\delta_{32}-\frac{\pi}{2}\right)\\\\
&\qquad+Y_{33}V^2_3\cos\left(\delta_{33}+\frac{\pi}{2}\right)\\\\
&=Y_{31}V_1V_3\sin\delta_{31}+Y_{32}V_2V_3\sin\delta_{32}-Y_{33}V^2_3\sin\delta_{33}\\\\
&=10\cdot1.0\cdot V_3\sin\left(\delta_3-\delta_1\right)+10\cdot1.05\cdot V_3\sin\left(\delta_3-\delta_2\right)\\\\
&\qquad-18V^2_3\sin\left(\delta_3-\delta_3\right)\\\\
&=10V_3\sin\delta_3+10.5V_3\sin\left(\delta_3-\delta_2\right) ・・・(9)
\end{align*}$$
Q_3&=\displaystyle \sum_{l=1}^3Y_{3l}V_lV_3\sin\left(\delta_{3l}-\theta_{3l}\right)\\\\
&=Y_{31}V_1V_3\sin\left(\delta_{31}-\theta_{31}\right)+Y_{32}V_2V_3\sin\left(\delta_{32}-\theta_{32}\right)\\\\
&\qquad+Y_{33}V^2_3\sin\left(\delta_{33}-\theta_{33}\right)\\\\
&=Y_{31}V_1V_3\sin\left(\delta_{31}-\frac{\pi}{2}\right)+Y_{32}V_2V_3\sin\left(\delta_{32}-\frac{\pi}{2}\right)\\\\
&\qquad+Y_{33}V^2_3\sin\left(\delta_{33}+\frac{\pi}{2}\right)\\\\
&=-Y_{31}V_1V_3\cos\delta_{31}-Y_{32}V_2V_3\cos\delta_{32}+Y_{33}V^2_3\cos\delta_{33}\\\\
&=-10\cdot1.0\cdot V_3\cos\left(\delta_3-\delta_1\right)-10\cdot1.05\cdot V_3\cos\left(\delta_3-\delta_2\right)\\\\
&\qquad+18V^2_3\cos\left(\delta_3-\delta_3\right)\\\\
&=-10V_3\cos\delta_3-10.5V_3\cos\left(\delta_3-\delta_2\right)+18V^2_3 ・・・(10)
\end{align*}$$
このとき、先程設定した初期値$\delta^{\left(0\right)}_2,\ \delta^{\left(0\right)}_3,\ V^{\left(0\right)}_3$における計算値$P^{\left(0\right)}_2,\ P^{\left(0\right)}_3,\ Q^{\left(0\right)}_3$は、$(8)\sim(10)$式より、
P^{\left(0\right)}_2&=10.5\sin\delta^{\left(0\right)}_2+10.5V^{\left(0\right)}_3\sin\left(\delta^{\left(0\right)}_2-\delta^{\left(0\right)}_3\right)=0\\\\
P^{\left(0\right)}_3&=10V_3\sin\delta^{\left(0\right)}_3+10.5V^{\left(0\right)}_3\sin\left(\delta^{\left(0\right)}_3-\delta^{\left(0\right)}_2\right)=0\\\\
Q^{\left(0\right)}_3&=-10V^{\left(0\right)}_3\cos\delta^{\left(0\right)}_3-10.5V^{\left(0\right)}_3\cos\left(\delta^{\left(0\right)}_3-\delta^{\left(0\right)}_2\right)+18V^{\left(0\right)2}_3\\\\
&=-10-10.5+18=-2.5
\end{align*}$$
また、指定値については、計算条件より、
$$\begin{cases}
P_{2s}&=0.7\\\\
P_{3s}&=-1.0\\\\
Q_{3s}&=-0.2
\end{cases}$$
したがって、指定値と計算値との誤差$\Delta P^{\left(0\right)}_2,\ \Delta P^{\left(0\right)}_3,\ \Delta Q^{\left(0\right)}_3$は、
\Delta P^{\left(0\right)}_2&=P_{2s}-P^{\left(0\right)}_2=0.7-0=0.7\\\\
\Delta P^{\left(0\right)}_3&=P_{3s}-P^{\left(0\right)}_3=-1.0-0=-1.0\\\\
\Delta Q^{\left(0\right)}_3&=Q_{3s}-Q^{\left(0\right)}_3=-0.2-\left(-2.5\right)=2.3
\end{align*}$$
②ヤコビ行列の計算
ニュートン・ラフソン法において、1回目（$i=0$）の収束計算における修正方程式は次のようになる。
\left(\begin{array}{c} \Delta P^\left(0\right)_2 \\ \Delta P^\left(0\right)_3 \\ \Delta Q^\left(0\right)_3 \end{array}\right)&=\left(\begin{array}{ccc}\displaystyle\frac{\partial P_2}{\partial \delta_2} & \displaystyle\frac{\partial P_2}{\partial \delta_3} & \displaystyle\frac{\partial P_2}{\partial V_3} \\ \displaystyle\frac{\partial P_3}{\partial \delta_2} & \displaystyle\frac{\partial P_3}{\partial \delta_3} & \displaystyle\frac{\partial P_3}{\partial V_3} \\ \displaystyle\frac{\partial Q_3}{\partial \delta_2} & \displaystyle\frac{\partial Q_3}{\partial \delta_3} & \displaystyle\frac{\partial Q_3}{\partial V_3} \end{array}\right)\left(\begin{array}{c} \Delta\delta^\left(1\right)_2 \\ \Delta\delta^\left(1\right)_3 \\ \Delta V^\left(1\right)_3\end{array}\right)\\\\
&=\boldsymbol{J}^{\left(0\right)}\left(\begin{array}{c} \Delta\delta^\left(1\right)_2 \\ \Delta\delta^\left(1\right)_3 \\ \Delta V^\left(1\right)_3\end{array}\right) ・・・(11)\\\\
\left(\begin{array}{c} \delta^\left(1\right)_2 \\ \delta^\left(1\right)_3 \\ V^\left(1\right)_3 \end{array}\right)&=\left(\begin{array}{c} \delta^\left(0\right)_2 \\ \delta^\left(0\right)_3 \\ V^\left(0\right)_3\ \end{array}\right)+\left(\begin{array}{c} \Delta\delta^\left(1\right)_2 \\ \Delta\delta^\left(1\right)_3 \\ \Delta V^\left(1\right)_3\end{array}\right) ・・・(12)
\end{align*}$$
なお、$(12)$式のヤコビ行列$\boldsymbol{J}^{\left(0\right)}$の係数はすべて$\delta^{\left(0\right)}_2,\ \delta^{\left(0\right)}_2,\ V^{\left(0\right)}_3$における値とする。
ここで、修正方程式の係数を計算する式は、$(8)\sim(10)$式をそれぞれの変数で偏微分して、
\frac{\partial P_2}{\partial \delta_2}&=10.5\cos\delta_2+10.5V_3\cos\left(\delta_2-\delta_3\right)\\\\
\frac{\partial P_2}{\partial \delta_3}&=-10.5V_3\cos\left(\delta_2-\delta_3\right)\\\\
\frac{\partial P_2}{\partial V_3}&=10.5\sin\left(\delta_2-\delta_3\right)\\\\
\frac{\partial P_3}{\partial \delta_2}&=-10.5V_3\cos\left(\delta_3-\delta_2\right)\\\\
\frac{\partial P_3}{\partial \delta_3}&=10.5V_3\cos\delta_3+10.5V_3\cos\left(\delta_3-\delta_2\right)\\\\
\frac{\partial P_3}{\partial V_3}&=10\sin\delta_3+10.5\sin\left(\delta_3-\delta_2\right)\\\\
\frac{\partial Q_3}{\partial \delta_2}&=-10.5V_3\sin\left(\delta_3-\delta_2\right)\\\\
\frac{\partial Q_3}{\partial \delta_3}&=10V_3\sin\delta_3+10.5V_3\sin\left(\delta_3-\delta_2\right)\\\\
\frac{\partial Q_3}{\partial V_3}&=-10\cos\delta_3+10.5\cos\left(\delta_3-\delta_2\right)+36V_3
\end{align} ・・・(13)$$
$(13)$式より、初期値$\delta^{\left(0\right)}_2,\ \delta^{\left(0\right)}_2,\ V^{\left(0\right)}_3$におけるヤコビ行列$\boldsymbol{J}^{\left(0\right)}$は、
上記の結果より、逆行列$\left[\boldsymbol{J}^{\left(0\right)}\right]^{-1}$を計算すると、
③収束判定と近似解の設定
求めた逆行列$\left[\boldsymbol{J}^{\left(0\right)}\right]^{-1}$を用いて、$(11)$式の修正方程式により誤差$\Delta\delta^\left(1\right)_2,\ \Delta\delta^\left(1\right)_3,\ \Delta V^\left(1\right)_3$を計算すると、
\left(\begin{array}{c} \Delta\delta^\left(1\right)_2 \\ \Delta\delta^\left(1\right)_3 \\ \Delta V^\left(1\right)_3\end{array}\right)&=\left[\boldsymbol{J}^{\left(0\right)}\right]^{-1}\left(\begin{array}{c} \Delta P^\left(1\right)_2 \\ \Delta P^\left(1\right)_3 \\ \Delta Q^\left(1\right)_3 \end{array}\right)\\\\
&=\left(\begin{array}{ccc}0.064012 & 0.032787 & 0 \\ 0.032787 & 0.065574 & 0 \\ 0 & 0 & 0.064516 \end{array}\right)\left(\begin{array}{c} 0.7 \\ -1.0 \\ 2.3 \end{array}\right)\\\\
&=\left(\begin{array}{c} 0.012022 \\ -0.042623 \\ 0.148387 \end{array}\right)
\end{align*}$$
上記より、誤差の大きさは最大で$\left|\Delta V^\left(1\right)_3\right|=0.148387>1\times10^{-7}$となることから、計算フロー上の収束判定はNGとなり、また①のプロセスに戻って計算を続行する。
また、1回目の収束計算における近似解$\delta^\left(1\right)_2,\ \delta^\left(1\right)_3,\ V^\left(1\right)_3$は、$(12)$式より、
\left(\begin{array}{c}\delta^\left(1\right)_2 \\ \delta^\left(1\right)_3 \\ V^\left(1\right)_3\end{array}\right)&=\left(\begin{array}{c} \delta^\left(0\right)_2 \\ \delta^\left(0\right)_3 \\ V^\left(0\right)_3\end{array}\right)+\left(\begin{array}{c} \Delta\delta^\left(1\right)_2 \\ \Delta\delta^\left(1\right)_3 \\ \Delta V^\left(1\right)_3\end{array}\right)\\\\
&=\left(\begin{array}{c} 0 \\ 0 \\ 1.0 \end{array}\right)+\left(\begin{array}{c} 0.012022 \\ -0.042623 \\ 0.148387 \end{array}\right)\\\\
&=\left(\begin{array}{c} 0.012022 \\ -0.042623 \\ 1.148387 \end{array}\right)
\end{align*}$$
計算結果
上記の①〜③の計算プロセスを、誤差の値が$\varepsilon=1\times10^{-7}$以下になるまで繰り返し、次の結果が得られた。
今回は5回目の計算（$i=4$）で収束する（誤差の最大値は$\left|\Delta\delta_3\right|=2.1812\times10^{-8}<\varepsilon$）という結果となった。
上記より、得られた解$\delta_2,\ \delta_3,\ V_3$を元に算出した電圧$\dot{V}_2,\ \dot{V}_3$の値は、次のようになる。
$$\begin{cases}
\dot{V}_2&=1.05\angle0.011542\\\\
\dot{V}_3&=1.12795\angle-0.037349
\end{cases}$$
なお、上記を直交座標に直すと、
\dot{V}_2&=1.05\left(\cos0.011542+j\sin0.011542\right)\\\\
&=1.04993+j0.012186\\\\
\dot{V}_3&=1.12795\left\{\cos\left(-0.037349\right)+j\sin\left(-0.037349\right)\right\}\\\\
&=1.12717-j0.042119
\end{cases}$$
となり、直交座標表示における計算結果と一致していることがわかる。
本記事では、潮流計算の手法の1つである直流法について解説する。直流法の考え方直流法とは直流法は、電力系統の潮流分布（有効電力）について、収束計算を用いずに1回の計算のみで近似値を求める手法である。この直流法に対し[…]
関連する例題（「電験王」へのリンク）
電験一種
参考文献
- 電気学会『電気工学ハンドブック 第7版』オーム社，2013
- 関根泰次『電力系統過渡解析論』オーム社，1984
- 新田目倖造『電力系統技術計算の基礎』電気書院，1981
- 加藤政一『詳解電力系統工学』東京電機大学出版局，2017
著書・製品のご紹介
『書籍×動画』が織り成す、未だかつてない最高の学習体験があなたを待っている！
※本ページはプロモーションが含まれています。―『書籍×動画』が織り成す、未だかつてない最高の学習体験があなたを待っている― 当サイト「電気の神髄」をいつもご利用ありがとうございます。管理人の摺り足の加藤です。[…]
この講座との出会いは、数学が苦手なあなたを救う！
電験アカデミアにテキストを書き下ろしてもらい、電験どうでしょうの川尻将先生により動画解説を行ない、電験3種受験予定者が電…
すべての電験二種受験生の方に向けて「最強の対策教材」作りました！
※本ページはプロモーションが含まれています。すべての電験二種受験生の方に向けて「最強の対策教材」作りました！ 当サイト「電気の神髄」をいつもご愛読ありがとうございます。管理人の摺り足の加藤です。 […]
初学者が躓きがちなギモンを、電験アカデミアがスッキリ解決します！
※本ページはプロモーションが含まれています。 当サイト「電気の神髄」をいつもご利用ありがとうございます。管理人の摺り足の加藤です。 2022年5月18日、オーム社より「電験カフェへようこそ[…]