這篇文章會從克西荷夫電流定律、歐姆定律和電功率公式中推導出電力潮流分析的基本公式，並說明常用的求解工具。一些電力系統的基本原理和名詞請參考電力系統運作基本原理。

Youtube | Power Flow Analysis (1): The Theory

電力潮流分析的基本公式

我們考慮一個電力網絡，其中有$ m $個線路和$ n $個節點。我們定義一$m \times n$的事件矩陣（incidence matrix）$[N]$：當線路#$i$的起點為節點#$j$時，該矩陣中第$i$列第$j$行之值為1；當線路#$i$的終點為節點#$j$時，該矩陣中第$i$列第$j$行之值為-1；其他情況下該矩陣之各單項之值為0。起點和終點係由線路的方向所決定，每條線路都有兩種方向可以任意選定。

根據克西荷夫電流定律，我們可以寫出下列關係式：

$$ \begin{equation} \{I_0\}^\mathsf{T}+\{I\}^\mathsf{T}[N]=0 \end{equation} $$

其中${I}$為各線路上的電流、${I_0}$則為各節點上流入或流出所考慮電力網絡外的電流（又稱為漏電流）。

同時，如果我們定義一$m \times m$傳導矩陣（admittance matrix）$[Y]$（這是阻抗矩陣impedence matrix的反矩陣），則根據歐姆定律可以寫出下列關係式：

$$ \begin{equation} \{I\}=[Y][N]\{U\} \end{equation} $$

其中${U}$為各節點上的節點電壓。當所考慮之線路中有電池或其他電壓源（voltage source）時，我們會以Thevenin-Norton理論求出等效併聯之電流源（current source），置入${I_0}$當中。

將克西荷夫電流定律與歐姆定律兩者衍生的關係式合併，我們有底下的節點壓流關係式：

$$ \begin{equation} [N]^\mathsf{T}[Y][N]\{U\}=-\{I_0\} \end{equation} $$

或者

$$ \begin{equation} [Y_{eq}]\{U\}=\{I_{eq}\} \end{equation} $$

因此，各節點之等效傳導矩陣和等效電流構成的線性系統的解，便是節點電壓。

補充說明：$[Y_{eq}]$ 的性質

$[Y_{eq}]$為$n \times n$之方矩陣，但其rank數最多僅有$n-1$。證明方式如下：

將$[N]$的第i列記為${r}_i$，則

$$ \begin{equation} [Y_{eq}] = \sum^m_{i=1} y_{i} \{r\}_i \otimes \{r\}_i \end{equation} $$$$ \begin{equation} \left( \{r\}_i \otimes \{r\}_{i} \right)_{jk} = \left\{\begin{array}{rcll} 1 &,& \: \left( \{r\}_i \right)_{j} \neq 0 &\&& \:j = k \\ -1 &,& \: \left( \{r\}_i \right)_{j} \left( \{r\}_i \right)_{k} \neq 0 &\&& \:j \neq k \\ 0 &,& \: \left( \{r\}_i \right)_{j} \left( \{r\}_i \right)_{k} = 0 \end{array}\right. \end{equation} $$

這代表$\{r\}_i \otimes \{r\}_{i}$的每一行（或每一列）各項相加的值必為0。由於$[Y_{eq}]$是各個$\{r\}_i \otimes \{r\}_{i}$的線性疊加，這代表$[Y_{eq}]$的任一行（或任一列）可透過高斯消去法完全消除。因此，$[Y_{eq}]$之rank數最多為$n-1$。事實上，在各節點彼此間至少有一條路徑相連接的系統中，$[Y_{eq}]$之rank數恰為$n-1$。證明方式如下：

$[Y_{eq}]$的rank必須和$[N]$相同：如果存在一向量$\{x\}$使$[Y_{eq}] \{x\} = 0$，則$\{x\}^\mathsf{T}([Y_{eq}]\{x\}) = \{x\}^\mathsf{T}[N]^\mathsf{T}[Y][N]\{x\} = \{[\sqrt Y][N]\{x\}\}^\mathsf{T}\{[\sqrt Y][N]\{x\}\}=0$，此意味$[N]\{x\}=0$；另一方面，如果$[N]\{x\}=0$，$[Y_{eq}] \{x\} = [N]^\mathsf{T}[Y]([N]\{x\}) = 0$；亦即，$\mathrm{null}([Y_{eq}]) = \mathrm{null}(N)$ 。而根據Rank–nullity theorem，$\mathrm{null}(N) = n - \mathrm{rank}([N])$，而$\mathrm{rank}([Y_{eq}]) = n- \mathrm{null}([Y_{eq}]) = \mathrm{rank}([N])$。
由於系統中各節點至少都會接上一條線段，這代表系統中至少有$n-1$條線段。當$m=n-1$時，$[N]_{ij}$可以整理成$\delta(i-j) - \delta(i - j + 1)$，此時$[N]$的rank數為$n-1$，因此$[Y_{eq}]$的rank數亦為$n-1$。另一方面，當$m \geq n$時，考慮底下的線路調整：假設線路$\# \alpha$連接節點$\#a$至$\#b$，線路$\# \beta$連接節點$\#a$至$\#c$，線路$\# \gamma$連接節點$\#b$至$\#c$，若我們將線路$\# \beta$的起始點改成節點$\#b$，此操作等同於對$[N]$做以下高斯消去法：將$\# \beta$行減去$\# \alpha$行，以使$\# \beta$行和$\# \gamma$相等。利用同樣的程序，我們可以將整個$[N]$漸進化簡成只有$n-1$種互相線性獨立的行向量，因此$[N]$的rank數為$n-1$，$[Y_{eq}]$的rank數亦為$n-1$。

Rank of Laplacian matrix

從以上的討論可知，在各節點彼此間至少有一條路徑相連接的系統中，必須先給定一節點之參考電壓，線性系統$[Y_{eq}]\{U\}=\{I_{eq}\}$才有唯一解。此數學定理的物理意義為：線路上電流和電流的移動，僅和各節點的相對電壓差有關，因此某節點上的參考電壓可任意設定，不影響最後的求解。同理也可以知道，當有些節點彼此之間沒有路徑相連時，$[Y_{eq}]$的rank會下降（因為彼此獨立的子系統都需要先設定一參考電壓）；網絡理論靠訴我們$\mathrm{Rank}([Y_{eq}]) = \mathrm{Rank}(N) = n - c(G)$(注)，其中$c(G)$為整個系統中所具有的元件（Component，在我們的應用中即獨立運作的子系統）數，此和我們的物理直覺相符合。
（注：其實利用我們上述證明就可以得到此結果：由於將每個獨立子系統對應在$Y_{eq}$的子方陣的rank求出並相加，就是整個系統的$Y_{eq}$的rank。）

而各節點上等效電流造成的電力潮流，則可利用電功率公式求出：

$$ \begin{equation} \{S\}=\{U\} \odot \{I_{eq}\}^* = \{U\} \odot \left([Y_{eq}]^*\{U\}^*\right) \end{equation} $$

其中$\{S\}$為各節點上等效電流造成的電力潮流。此即電力潮流公式（power flow equation）。為避免冗贅，本文底下不再寫出$[Y_{eq}]$的下標${}_{eq}$。

電力潮流分析的方法

電力潮流公式中，每個節點有兩個未知複數（電力潮流和電壓），但只有一個等式。如果將電力潮流分解成有效電力和無效電力、電壓分解成絕對值與相位角，則每節點有4個未知數和兩個等式需求解，故單只有電力潮流公式時，系統為未定（underdertermined）狀態。實際進行電力潮流分析時，我們需對每個節點的4個未知數中任意兩個事先做出假設，所以電力潮流分析中會有$C^4_2=6$種基本節點，雖然底下3種最為常見。

$P-Q$控制：節點輸出／入的有效電力和無效電力為事前已知。通常在模擬負載和沒有電網建構或支援控制功能的分散式再生能源電廠時使用。
$P-U$控制：節點電壓絕對值和輸出／入有效電力為事前已知。通常在模擬傳統電廠和進行電網建構或支援控制的分散式再生能源電廠時使用。
$U-\theta$控制：節點電壓絕對值和相位角為事前已知。因為電力系統中線路的電流僅受電壓差影響，因此有必要選出至少一個節點的電壓絕對值和相位角，當作其他節點的參考基準。理論上，電力網絡中所有能進行電網建構控制的電廠，其節點電壓絕對值和相位角應等同於此參考基準，因此這些電廠也都能採取$U-\theta$控制；實際模擬時這可能會導致系統沒有可行解，因此通常僅會選一節點做參考基準。

牛頓求解法

牛頓求解法即透過牛頓法（Newton’s method）找出同時滿足電力潮流公式和事先假設的節點電壓或電力潮流。假設電路上參考節點外每個節點都設成$P-Q$控制，則在一解附近可以對系統做線性化處理：

$$ \begin{equation} \left\{\begin{array}{c} \Delta P \\ \Delta Q \end{array}\right\}= [J]\left\{\begin{array}{c} \Delta ||U|| \\ \Delta \theta \end{array}\right\} \end{equation} $$

其中$\Delta P$和$\Delta Q$為特定節點電壓下，對應到的有效與無效電力和原先設定的差值：

$$ \begin{equation} \Delta P_k = -P_k + ||U_k|| \sum^n_{i=1} ||U_i|| \:||Y_{ki}|| \:\cos(\theta_{U_k}-\theta_{U_i}-\theta_{Y_{ki}}) \end{equation} $$$$ \begin{equation} \Delta Q_k = -Q_k + ||U_k|| \sum^n_{i=1} ||U_i|| \:||Y_{ki}|| \:\sin(\theta_{U_k}-\theta_{U_i}-\theta_{Y_{ki}}) \end{equation} $$

$[J]$則為Jacobian矩陣，即

$$ \begin{equation} [J]=\left[\begin{array}{cc} \frac{\partial \Delta P}{\partial ||U||} & \frac{\partial \Delta P}{\partial \theta}\\ \\ \frac{\partial \Delta Q}{\partial ||U||} & \frac{\partial \Delta Q}{\partial \theta} \end{array}\right] \end{equation} $$

這裡面的4種偏微分關係如下：

$$ \begin{equation} \begin{array} {ccl} \frac{\partial \Delta P_k}{\partial ||U_i||} &=& (1 - \delta_{ki}) \:||U_k|| \:||Y_{ki}|| \:\cos(\theta_{U_k}-\theta_{U_i}-\theta_{Y_{ki}}) \\ &+& \delta_{ki} \sum^n_{l = 1} (1 + \delta_{kl}) \:||U_l|| \:||Y_{kl}|| \:\cos(\theta_{U_k}-\theta_{U_l}-\theta_{Y_{kl}}) \end{array} \end{equation} $$$$ \begin{equation} \begin{array} {ccl} \frac{\partial \Delta Q_k}{\partial ||U_i||} &=& (1 - \delta_{ki}) \:||U_k|| \:||Y_{ki}|| \:\sin(\theta_{U_k}-\theta_{U_i}-\theta_{Y_{ki}}) \\ &+& \delta_{ki} \sum^n_{l = 1} (1 + \delta_{kl}) \:||U_l|| \:||Y_{kl}|| \:\sin(\theta_{U_k}-\theta_{U_l}-\theta_{Y_{kl}}) \end{array} \end{equation} $$$$ \begin{equation} \begin{array} {ccl} \frac{\partial \Delta P_k}{\partial \theta_{U_i}} &=& (1 - \delta_{ki}) \:||U_k|| \:||U_i|| \:||Y_{ki}|| \:\sin(\theta_{U_k}-\theta_{U_i}-\theta_{Y_{ki}}) \\ &-& \delta_{ki} \sum^n_{l = 1} (1 - \delta_{kl}) \:||U_k|| \:||U_l|| \:||Y_{kl}|| \:\sin(\theta_{U_k}-\theta_{U_l}-\theta_{Y_{kl}}) \end{array} \end{equation} $$$$ \begin{equation} \begin{array} {ccl} \frac{\partial \Delta Q_k}{\partial \theta_{U_i}} &=& -(1 - \delta_{ki}) \:||U_k|| \:||U_i|| \:||Y_{ki}|| \:\cos(\theta_{U_k}-\theta_{U_i}-\theta_{Y_{ki}}) \\ &+& \delta_{ki} \sum^n_{l = 1} (1 - \delta_{kl}) \:||U_k|| \:||U_l|| \:||Y_{kl}|| \:\cos(\theta_{U_k}-\theta_{U_l}-\theta_{Y_{kl}}) \end{array} \end{equation} $$

因此，給定各節點電壓絕對值與相位角初始猜測$\{||U_{old}||\}$和$\{\theta_{old}\}$後，我們可以計算此時的Jacobian矩陣，並且更新電壓絕對值與相位角的猜測：

$$ \begin{equation} \left\{\begin{array}{c} ||U_{new}|| \\ \theta_{new} \end{array}\right\} = \left\{\begin{array}{c} ||U_{old}|| \\ \theta_{old} \end{array}\right\} - [J]^{-1}\left\{\begin{array}{c} \Delta P \\ \Delta Q \end{array}\right\} \end{equation} $$

這樣持續迭代數次，順利的話便可收斂至有物理意義的可行解上。

牛頓求解法簡單易懂，但有底下幾個缺陷：首先，此法為局部性求解法，亦即每次初始條件不同時，必須重新求出電壓／相位對有效／無效電力的敏感度矩陣，影響求解效率；再來，因為其為局部性求解法，無法保證迭代過程一定收斂至有應用意義的高壓解上，而即使迭代結果發散也無法因而確定可行解不存在。凡此種種，皆讓此法在實際應用上，有許多隱憂。

HELM求解法

Youtube | Power Flow Analysis (3): Holomorphic embedding load flow method (HELM)

為了解決牛頓求解法分散性、非單一解的問題，在2012年，一種基於複數分析技巧的新潮流分析方法－全形嵌入潮流分析法（Holomorphic embedding load flow method，HELM）被發展了。這個潮流分析方法透過設置與電壓相對的鏡像變數，將原來的物理流問題嵌入一全形函數中，保證了最後求出的解具有明確性（univocality）：分析後無解時，必然代表具應用價值的解不存在；而若分析後有解，其亦則必定是具有實際應用意義的高電壓解（可參考Fundamentals of the Holomorphic Embedding Load-Flow Method，Antonio Trias，2015年）。

Holomorphic Embedding Load-flow method - Wikipedia

HELM的基本精神如下：假設和前節相同，電路上參考節點外每個節點都設成$P-Q$控制，則我們寫出一組嵌入系統。

$$ \begin{equation} \sum^{n}_{k=1} Y_{ik} U_k(s) = s \frac{S^*_i}{\hat U_{i}(s)} \end{equation} $$$$ \begin{equation} \sum^{n}_{k=1} Y^*_{ik} \hat U_{k}(s) = s \frac{S_i}{U_{i}(s)} \end{equation} $$

這會確保$U(s)$和$\hat{U}(s)$皆為複數平面上的全形函數（holomorphic function），使求解結果具有前述的明確性。而為了確保解有物理意義，在$s=1$時，$U(s)$和$\hat{U}(s)$必須滿足底下的反射條件（reflection condition）：

$$ \begin{equation} \{\hat{U}\}(s) = \{U^*\}(s^*) \end{equation} $$

補充說明：$U(s)$和$\hat{U}(s)$滿足全形函數的定義且可還原原始電力潮流問題

假設$U(s)$滿足全形函數之定義，我們會有

$$ \begin{equation} \frac{\partial \Re(U)}{\partial \Re(s)} = \frac{\partial \Im(U)}{\partial \Im(s)} \end{equation} $$

以及

$$ \begin{equation} \frac{\partial \Re(U)}{\partial \Im(s)} = -\frac{\partial \Im(U)}{\partial \Re(s)} \end{equation} $$

假設我們直接取$U(s)$的共軛複數$U^*(s)= \Re(U(s)) - j\Im(U(s))$，則

$$ \begin{equation} \frac{\partial \Re(U^*)}{\partial \Re(s)} = \frac{\partial \Re(U)}{\partial \Re(s)} = \frac{\partial \Im(U)}{\partial \Im(s)} = -\frac{\partial \Im(U^*)}{\partial \Im(s)} \end{equation} $$$$ \begin{equation} \frac{\partial \Re(U^*)}{\partial \Im(s)} = \frac{\partial \Re(U)}{\partial \Im(s)} = -\frac{\partial \Im(U)}{\partial \Re(s)} = \frac{\partial \Im(U^*)}{\partial \Re(s)} \end{equation} $$

很顯然地，$U^*(s)$並非全形函數。另一方面

$$ \begin{equation} \frac{\partial \Re(U^*)}{\partial \Re(s^*)} = \frac{\partial \Re(U)}{\partial \Re(s)} = \frac{\partial \Im(U)}{\partial \Im(s)} = \frac{\partial \Im(U^*)}{\partial \Im(s^*)} \end{equation} $$$$ \begin{equation} \frac{\partial \Re(U^*)}{\partial \Im(s^*)} = -\frac{\partial \Re(U)}{\partial \Im(s)} = \frac{\partial \Im(U)}{\partial \Re(s)} = -\frac{\partial \Im(U^*)}{\partial \Re(s^*)} \end{equation} $$

因此$U^*(s^*)$為全形函數。這給了我們設定$\hat{U}(s)$以及$U(s)$和$\hat{U}(s)$的反射條件的理由。而根據反射條件，在$s=1$時

$$ \begin{equation} \sum^{n}_{k=1} Y_{ik} U_k(1) = \frac{S^*_i}{\hat U_{i}(1)} = \frac{S^*_i}{U^*_{i}(1)} \end{equation} $$$$ \begin{equation} \sum^{n}_{k=1} Y^*_{ik} \hat U_{k}(s) = \sum^{n}_{k=1} Y^*_{ik} U^*_{i}(1) = \frac{S_i}{U_{i}(1)} \end{equation} $$

這兩個式子皆和原本的電力潮流問題一樣。

現在，我們把4個變數寫成$s$的多項式：$\{U\}(s)= \sum^{\infty}_{j=0}\{a_j\}s^j$、$\{\frac{1}{\hat{U}(s)}\}= \sum^{\infty}_{j=0}\{b_j\}s^j$、$\{\hat{U}\}(s)= \sum^{\infty}_{j=0}\{c_j\}s^j$、$\{\frac{1}{U(s)}\}= \sum^{\infty}_{j=0}\{d_j\}s^j$。

我們從$s=0$出發，這個點代表的是電路中完全沒有電流、各節點電壓皆相同的非零顯明高壓解（non-zero trivial high voltage solution；由於$[Y_{eq}]$的rank為n-1，故必存有一非零向量$\{\eta\}$使$[Y_{eq}]\{\eta\} = 0$，此$\{\eta\}$即為非零顯明高壓解）。故我們可知$\{a\}_0 = \{U_1\}$、$\{c\}_0 = \{U^*_1\}$、$\{a\}_0 \odot \{d\}_0 = \{b\}_0 \odot \{c\}_0 = 1$（假設節點$\#1$便是參考節點）。

對於一次項係數而言，我們有

$$ \begin{equation} \sum^{n}_{k=1} Y_{ik} a_{1,\:k} = S^*_i \:b_{0, \:i} \end{equation} $$$$ \begin{equation} \sum^{n}_{k=1} Y^*_{ik} c_{1, \:k} = S_i d_{0, \:i} \end{equation} $$

這組線性系統可以解出$\{a\}_1$和$\{c\}_1$，同時從倒數關係我們也有

$$ \begin{equation} \{a\}_1 \odot \{d\}_0 + \{a\}_0 \odot \{d\}_1 = \{b\}_0 \odot \{c\}_1 + \{b\}_1 \odot \{c\}_0 = 0 \end{equation} $$

如此一來也能求出$\{b\}_1$和$\{d\}_1$。接下來二式項、三次項的係數也都可以用類似的方法求出，直到任意多項係數為止。

在求出任意多項係數後，我們便有節點電壓在顯明解附近的多項式近似。如果取$s=1$而此多項式仍然收斂，則我們便有節點電壓之數值近似解。這代表在$s=0$和$s=1$之間全形函數$U(s)$必須是解析連續（analytic continuous）的；確保這件事情的充必條件為$|s|=1$的圓盤內不存在奇異點（singularity point）。使用近對角的Padé近似法（near-diagonal Padé approximant）將多項式近似替代成代數分式，除了可以加速收斂的速度（若解存在），也可以透過分母多項式的根，找出函數的奇異點¹。

Youtube | Padé Approximants

HELM雖然有諸多優勢，但也有缺點：就數值分析來說，當Padé近似法求出的分母含有許多根的時候，求解器可能會出現數值不穩定的問題，造成判讀解析連續性值的困難。

電力系統潮流分析的方法與應用

P-U控制下的HELM求解方法

前述的HELM求解法，在面臨P-U控制的節點時，必須做出調整。首先是在這些節點的電力潮流公式中增加一道等式

$$ \begin{equation} U_k(s) \hat{U}_k(s) = U_1^2 + s(U_k^2 - U_1^2) \end{equation} $$

這道等式會確保在$s=0$的時候求出的解為前節所提到之顯明解，但當$s=1$時則滿足該節點事先設定的電壓絕對值限制。同一時間，由於無效電力在P-U控制下的節點屬於自由變數，因此額外將其寫成$s$之多項式$\{Q\}(s)= \sum^{\infty}_{j=0}\{e_j\}s^j$。此時的嵌入系統被修正為

$$ \begin{equation} \sum^{n}_{k=1} Y_{ik} U_k(s) = \frac{sP_i -jQ_i(s)}{\hat U_{i}(s)} \end{equation} $$$$ \begin{equation} \sum^{n}_{k=1} Y^*_{ik} \hat U_{k}(s) = \frac{s P_i + jQ_i(s)}{U_{i}(s)} \end{equation} $$

從顯明解可以推得$\{e\}_0 = \{0\}$。

此系統中的一次項對應的等式為

$$ \begin{equation} \sum^{n}_{k=1} Y_{ik} a_{1,\:k} + j \:b_{0, \:i} \:e_{1, \:i}= P_i \:b_{0, \:i} \end{equation} $$$$ \begin{equation} \sum^{n}_{k=1} Y^*_{ik} c_{1, \:k} - j \:b_{0, \:i} \:e_{1, \:i} = P_i d_{0, \:i} \end{equation} $$$$ \begin{equation} \{a\}_1 \odot \{d\}_0 + \{a\}_0 \odot \{d\}_1 = \{b\}_0 \odot \{c\}_1 + \{b\}_1 \odot \{c\}_0 = 0 \end{equation} $$$$ \begin{equation} \{a\}_1 \odot \{c\}_0 + \{a\}_0 \odot \{c\}_1 = U_k^2 - U_1^2 \end{equation} $$

此時便必須同時解$\{a\}_1$、$\{b\}_1$、$\{c\}_1$、$\{d\}_1$、和$\{e\}_1$所構成的線性系統，會比純P-Q控制的案例中更費時，但方法學上並無太大差異－一次項解完之後解二次項，如此漸次操作直到任意多項時，再以Padé近似法求出代數分式並確認節點電壓是否收斂。

底下為能同時存在P-Q和P-U控制節點的HELM電力潮流求解器程式碼。

Tony Yen

Other Platforms

Blogs

電力潮流分析原理與方法

電力潮流分析的基本公式

電力潮流分析的方法

牛頓求解法

HELM求解法

P-U控制下的HELM求解方法

延伸閱讀