智能配电网络建模与分析
配电网状态估计
2020年
本次课的内容
（1）了解配电网状态估计的背景、现状、作用；
（2）用一个简单例子说明状态估计的基本方法；
（3）极大似然加权最小二乘法；
（4）一种以合格率最大为目标的状态估计新方法；
背景
用户对电能质量和可靠性的要求越来越高，配电SCADA系统的应用越来越广；
监控电力系统的功率潮流和电压，对于维持系统的安全性、稳定性非常重要；
出于经济性的考虑，量测设备安装的数目是有限的；
实时数据不足；
传送到控制中心的数据不准确、不正常或者具有时延。
配电网状态估计是在获知全网网络结构的条件下，结合从馈线FTU和母线RTU得到的实时功率和电压信息，补充负荷预测数据以及抄表数据，运用新型的数学和计算机手段估计配电网用户实时负荷，由此获得全网当前时刻各部分的运行状态和参数的过程，可以为其他配电系统高级应用软件提供可靠的实时数据信息。
配电网状态估计特点
配电网状态估计具有如下不同于输电网状态估计的特点：
（1）配电网中的量测量包括实时量测、伪量测和零注入量测等，由于配电馈线分支数量庞大，不可能对所有馈线分支配置实时量测，往往利用用户数据库中的数据得到的预测值作为伪量测使用。显然这种量测的误差较大，且随时间段的不同而发生变化；而实时量测的标准差则较小，这就需要每次状态估计之后，对于量测的权重重新赋值，以改善预测可信度。
（2）关于第一次状态估计伪量测权重值，可以根据各类用电的预测可信度以及量测值来源对不同类型的伪量测赋予相应权值。
（3）在配电网中可以对电流幅值与电压实时测量，这些量测数量较大，而且较容易获得，可以有效提高配电网状态估计的精度。
（4）在辐射状配电网状态估计中变电站出口电压一般视为精确值而不参与状态估计，并且馈线之间除根节点以外没有电气联系，所以各条馈线可以分别进行状态估计，即以馈线作为状态估计基本单元。
配电网量测与调度系统
配电网量测与调度系统如图所示：
每个环节都可能对数据产生误差，导致测量误差。对一个经过良好的校对的量测系统来说，其误差具有正态分布的性质，即对于每一个量测量，量测误差标准差σ为正常量测范围的0.5%至2%。配电网调度中心收到的不良数据主要来自于测量与传送系统受到的较大随机干扰或偶然故障、配电网快速变化过程中各测量点间的非同时测量以及系统过渡过程。
配电网状态估计的主要作用
（1）提高数据精度。根据量测量的精度和基尔霍夫定律按最佳估计准则对生数据进行计算，得到最接近于系统真实状态的最佳估计值。
（2）提高数据系统的可靠性。对生数据进行不良数据的检测与辨识，对不良数据进行删除或修正。
（3）推算出完整而精确的各种电气量。例如根据周围相邻的变电站的量测量推算出某一装设有远动装置的变电站的各种电气量:或者根据现有类型的量测量推算其他难以量测的电气量，如根据有功功率量测值推算各节点电压的相角。
（4）网络结线辨识或开关状态辨识。根据遥测量估计电网的实际开关状态，纠正偶然出现的错误的开关状态信息，以保证数据库中电网结线方式的正确性。
（5）数据预测。可以应用状态估计算法以现有的数据预测未来的趋势和可能出现的状态。丰富数据库的内容，为安全分析与运行计划等提供必要的计算条件。
（6）参数估计。如果把某些可疑或未知的参数作为状态量处理时，也可以用状态估计的方法估计出这些参数的值。例如带负荷自动改变分接头位置的变压器， 如果分接头位置信号没有传送到中调，可以把它作为参数估计求出。
（7）确定合适的测点数量及其合理分布。通过状态估计的离线模拟试验，可以确定配电网合理的数据收集与传送系统，用以改进现有的远动系统或规划未来的远动系统，使软件与硬件联合以发挥更大的效益，既保证了数据的质量，又降低了整个系统的投资。
实际应用：估计计算和结线分析,简单的不良数据检测与辨识、结线辨识以及在线应用的功能。
综上所述，配电网状态估计是远动装置与数据库之间的重要一环。通过配电网状态估计使从远动装置接收的低精度、不完整的数据，转变成最终数据库中高精度的、完整且可靠的数据。状态估计不仅提高了数据精度，滤掉了不良数据，还补充了一些测点，能够得到某些难以直接量测的物理量（如节点电压的相角）。
电力系统状态量测一个简单的例子
一个简单的直流潮流模型如图所示：
假设三母线直流系统在左图中所示的负荷和发电机出力条件下运行，系统的量测信息由3个功率潮流量测提供，量测放置位置如右图所示。在这些量测信息中，只需要两个就可以计算得到所有母线的相角以及所有负载和发电机的功率大小，也就是系统的当前状态。
量测结果无误差的情况
若使用\(M_{13}\)和\(M_{32}\)这两个量测值进行计算，并进一步假设它们与真实值的误差为0，即能够对潮流进行准确测量，它们的量测值如下：
\[M_{13}=5\;{\rm{MW}}, M_{32}=40\;{\rm{MW}}\]
由直流潮流公式可知，线路1-3传输的功率大小为：
\[f_{13}=\frac{1}{x_{13}}(\theta_1-\theta_3)=M_{13}=5\;{\rm{MW}}\]
线路3-2传输的功率大小为：
\[f_{32}=\frac{1}{x_{23}}(\theta_3-\theta_2)=M_{32}=40\;{\rm{MW}}\]
选择母线3为参考点，则\(\theta_3 = 0\)，代入上述方程解得：
\[\theta_1=0.02\;{\rm{rad}},\theta_2=-0.10\;{\rm{rad}}\]
量测结果有微小误差的情况
假设量测值分别为：
\[M_{13}=6\;{\rm{MW}},\ M_{32}=37\;{\rm{MW}},\ M_{12}=62\;{\rm{MW}}\]
如果按照上一小节的方法，仍然选择量测仪表\(M_{13}\)和\(M_{32}\)的两个量测值进行计算，可以求得相角为：
\[\theta_1=0.024\;{\rm{rad}},\theta_2=-0.0925\;{\rm{rad}},\theta_3=0\;{\rm{rad}}\]
左图：使用量测仪表\(M_{13}\)和\(M_{32}\)计算得到的潮流结果；线路1-2上的功率潮流计算结果与量测结果\(M_{12}\)不符
右图：忽略量测仪表\(M_{13}\)的测量结果，而采用量测仪表\(M_{12}\)和\(M_{32}\)的两个量测值进行计算得到的潮流结果。线路1-2上的功率潮流计算结果与量测结果\(M_{12}\)相符，但线路1-3上的功率潮流计算结果与量测结果\(M_{13}\)不相符。
=》需要采用一种方法，能够充分利用所有三个量测仪表的量测信息，得到系统的最佳状态估计值，包括相角、线路潮流、母线负荷和发电机输出功率等。
统计原理
统计估计是用量测值（采样值）估计电力系统中一个或多个状态量（未知参数）的过程。如何利用已知的量测值对状态量进行最佳的估计，研究人员总结出了最常用的三条准则：
（1）极大似然准则。
其目的是使得状态量的估计值\(\hat{x}\)是其真实值\(x\)的可能性最大。
（2）加权最小二乘准则。
其目的是最小化估计量测值\(\hat{z}\)与真实量测值\(z\)的加权偏差的平方和。
（3）最小方差准则。
其目的是使状态变量向量的估计值与真实值的偏差平方和的数学期望最小。
极大似然加权最小二乘状态估计法
假定测量仪器的误差服从正态无偏分布，则根据上述准则中的任一条都能得出相同的估计量。我们将采用极大似然法进行后续分析，该方法直接引入了测量误差的权重矩阵\(R\)。
极大似然法主要讨论对于已经通过测量仪器得到的量测值，在不同的状态变量描述的系统状态下出现这样量测值的概率是多少，而这个概率取决于测量仪器的随机误差和待估计的未知参数。
在应用极大似然估计法的过程中，选择使得上述概率最大的估计值作为估计的真实值，就完成了一次估计的过程，
极大似然估计法的前提是假设测量随机误差的概率密度函数已知。
除了极大似然估计法，其他的估计方法也可以使用，其中，最小二乘估计法不需要知道样本或测量误差的概率密度函数，若认为样本或测量误差的概率密度函数呈正态分布，则可以得到与极大似然估计法相同的估计公式。
随机测量误差是从测量设备得到的某个参数量测值和其真实值之间的偏差，量测值与真实值之间的关系如下式所示：
\[\begin{align}
z_{\text c}=z_{\text t}+ \eta
\end{align}\]
式中：\(z_{\text c}\)为量测值；\(z_{\text t}\)为真实值；\(\eta\)为随机测量误差，用来表征测量中的不确定性。
如果测量误差是无偏的，那么\(\eta\)的概率密度函数服从均值为0的正态分布。需要注意的是，其他的概率密度函数也可用极大似然法进行分析。\(\eta\)的概率密度函数表达式如下：
\[\begin{align}
f\left ( \eta \right )=\frac{1}{\sigma \sqrt{2\pi }}e^{-\frac{\eta ^{2}}{2\sigma ^{2}}}
\end{align}\]
式中：\(\sigma\)为标准差；\(\sigma ^{2}\)为方差；正态分布\(f\left ( \eta \right )\)的曲线如图所示。标准差\(\sigma\)反映了随机测量误差的大小，\(\sigma\)越大，测量就越不准确。这里采用正态分布来描述测量误差分布的原因是当多重因素作用在测量误差上时，整体的测量误差就趋向于呈现正态分布。
极大似然估计的一个例子
带有电流量测的直流电路如图所示。本例将使用带有误差的电流表来估算电压源的值\(V\)，其中电流表误差的标准差已知。电流表显示的值记为\(I_1\)，等于电路中电流的真实值\(I_{\text{1t}}\)加上电流表测量的随机误差\(\eta _1\)，即：
\[\begin{align}
I_1=I_{\text{1t}}+ \eta_1
\end{align}\]
由于\(\eta _1\)的均值为0，所以\(I_1\)与\(I_{\text{1t}}\)的均值相等。因此，\(I_1\)的概率密度函数为：
\[\begin{align}
f\left ( I_1 \right )=\frac{1}{\sigma_1 \sqrt{2\pi }}e^{-\frac{\left ( I_1-I_{\text {1t}} \right )^{2}}{2\sigma_1 ^{2}}}
\end{align}\]
式中：\(\sigma_1\)是随机误差\(\eta_1\)的标准差。
假设电路中电阻的值已知为\(r_1\)，那么上式可以改写成：
\[\begin{align}
f\left ( I_1 \right )=\frac{1}{\sigma_1 \sqrt{2\pi }}e^{-\frac{\left ( I_1-\frac{V}{r_1} \right )^{2}}{2\sigma_1 ^{2}}}
\end{align}\]
根据极大似然估计法的要求，现在要做的是找到\(V\)的估计值\(V_{\rm{est}}\)，使得出现电流表示数为\(I_1\)情况的概率最大。由于\(I_1\)的概率密度函数已知，则：
\[\begin{align}
p\left ( I_1 \right )=\int_{I_1}^{I_1+\text{d}I_1}f\left ( I_1 \right )\text{d}I_1=f\left ( I_1 \right )\text{d}I_1
\end{align}\]
上式中是将\(V\)作为变量，则\(f\left ( I_1 \right )\)为\(V\)的函数，而且当\(f\left ( I_1 \right )\)取最大值时\(p\left ( I_1 \right )\)最大。为了运算方便，可以取\(f\left ( I_1 \right )\)的自然对数，即
\[\begin{align}
\ln \left ( f \left ( I_1 \right ) \right )=-\left ( \ln\left ( \sigma _1\sqrt{2\pi } \right )+\frac{\left ( I_1-\frac{V}{r_1} \right )^{2}}{2\sigma _1^{2}} \right )
\end{align}\]
上式括号中第一项为常数，因此问题转化为求括号中第二项的最大值，对第二项求导得
\[\begin{align}
\label{qiudaogongshi1}
\frac{\text{d}}{\text{d}V}\left [ \frac{\left ( I_1-\frac{V}{r_1} \right )^{2}}{2\sigma _1^{2}} \right ]=-\frac{I_1-\frac{V}{r_1}}{r_1\sigma _1^{2}}=0
\end{align}\]
求解上式得到\(V\)的估计值\(V_{\rm{est}}\)：
\[V_{\rm{est}}=I_1 \times r_1\]
在原电路上再并联一条支路，得到的并联直流电路如下图所示。
假设两条支路的电阻均为已知，分别记为\(r_1\)、\(r_2\)。两个电流表的量测值分别记为\(I_1\)、\(I_2\)。由上例可得
\[\begin{align}
\begin{split}
I_1=I_{\text{1t}}+ \eta _1 \\
I_2=I_{\text{2 t}}+ \eta _2
\end{split}
\end{align}\]
\(I_{\text{1t}}\)、\(I_{\text{2t}}\)分别为对应支路电流的真实值，\(\eta _1\)、\(\eta _2\)分别为对应电流表的随机误差。两个误差的概率密度函数相互独立，均值为0，呈正态分布，因此：
\[\begin{align}
\begin{split}
f\left ( \eta _1 \right )=\frac{1}{\sigma_1 \sqrt{2\pi }}e^{-\frac{\eta_1 ^{2}}{2\sigma_1 ^{2}}}
f\left ( \eta _2 \right )=\frac{1}{\sigma_2 \sqrt{2\pi }}e^{-\frac{\eta_2 ^{2}}{2\sigma_2 ^{2}}}
\end{split}
\end{align}\]
由上例同理可得\(I_1\)、\(I_2\)的概率密度函数为
\[\begin{align}
\begin{split}
f\left ( I_1 \right )=\frac{1}{\sigma_1 \sqrt{2\pi }}e^{-\frac{\left ( I_1-\frac{V}{r_1} \right )^{2}}{2\sigma_1 ^{2}}}
f\left ( I_2 \right )=\frac{1}{\sigma_2 \sqrt{2\pi }}e^{-\frac{\left ( I_2-\frac{V}{r_2} \right )^{2}}{2\sigma_2 ^{2}}}
\end{split}
\end{align}\]
这里的似然函数是\(I_1\)、\(I_2\)量测值同时出现的概率，由于之前已经假定它们的随机误差\(\eta _1\)、\(\eta _2\)相互独立，所以：
\[\begin{align}
p\left ( I_1 I_2 \right )=p\left ( I_1 \right )\times p\left ( I_2 \right )=f\left ( I_1 \right ) f\left ( I_2 \right )\text{d}I_1\text{d}I_2\\
=\frac{1}{\sigma_1 \sqrt{2\pi }}e^{-\frac{\left ( I_1-\frac{V}{r_1} \right )^{2}}{2\sigma_1 ^{2}}}\times \frac{1}{\sigma_2 \sqrt{2\pi }}e^{-\frac{\left ( I_2-\frac{V}{r_2} \right )^{2}}{2\sigma_2 ^{2}}}\text{d}I_1\text{d}I_2
\end{align}\]
对上式取对数再求导，由极值的必要条件可得
\[\begin{align}
\label{qiudaogongshi2}
\frac{\text{d}}{\text{d}V}\left [ \frac{\left ( I_1-\frac{V}{r_1} \right )^{2}}{2\sigma _1^{2}}+\frac{\left ( I_2-\frac{V}{r_2} \right )^{2}}{2\sigma _2^{2}} \right ]=-\frac{I_1-\frac{V}{r_1}}{r_1\sigma _1^{2}}-\frac{I_2-\frac{V}{r_2}}{r_2\sigma _2^{2}}=0
\end{align}\]
解得
\[\begin{align}
V_{\rm{est}}=\frac{\frac{I_1}{r_1\sigma _1^{2}}+\frac{I_2}{r_2\sigma _2^{2}}}{\frac{1}{r_1^{2}\sigma _1^{2}}+\frac{1}{r_2^{2}\sigma _2^{2}}}
\end{align}\]
当其中一个电流表的方差远小于另一个，如\(\sigma _2^{2} \ll \sigma _1^{2}\)时，电压估计值为
\[V_{\rm{est}} \simeq I_2\times r_2\]
因此，通过极大似然估计法来估计未知参数，能够根据测量结果的好坏来给它们赋予不同的权重值，以期得到更准确的估计结果。
通过上面的例子中两个求导公式可以发现，极大似然估计法的目标函数是量测值与真实值误差平方和，求此目标函数的最小值。其中，真实值是未知参数的函数。因此，用\(n\)个量测值来估计参数\(V\)时，得到表达式如下,记为方程A：
\[\begin{eqnarray} \label{7s15}
\min_{V}J\left ( V \right )=\sum_{i=1}^{n}\frac{\left [ I_i-f_i\left ( V \right ) \right ]^{2}}{\sigma _i^{2}}
\end{eqnarray}\]
如果需要用\(n\)个量测值来估计\(m\)个参数，则表达式为方程B
\[\begin{eqnarray} \label{7s16}
\min_{V_1,V_2,...,V_m}J\left ( V_1,V_2,...,V_m \right )=\sum_{i=1}^{n}\frac{\left [ I_i-f_i\left ( V_1,V_2,...,V_m \right ) \right ]^{2}}{\sigma _i^{2}}
\end{eqnarray}\]
上面的算法称为加权最小二乘估计法，如前面所说，如果测量误差服从正态分布，则加权最小二乘估计法等价于极大似然估计法。
加权最小二乘法的矩阵表达形式
如果函数\(f_i(x_1,x_2,\cdots ,x_{N_{\text{s}}})\)是线性函数，则方程B有一个封闭解。将函数\(f_i(x_1,x_2,\cdots ,x_{N_{\text{s}}})\)写成：
\[\begin{eqnarray} \label{7s18}
f_i(x_1,x_2,\cdots ,x_{N_{\text{s}}}) = f_i(\boldsymbol x) = h_{i1}x_1 + h_{i2}x_2 + \cdots +h_{iN_{\text{s}}}x_{N_{\text{s}}}
\end{eqnarray}\]
将所有\(f_i\)函数写成向量形式，则
\[\begin{eqnarray} \label{7s19}
\boldsymbol f(\boldsymbol x) = \begin{bmatrix}
f_1(\boldsymbol x)\\
f_2(\boldsymbol x)\\
\vdots \\
f_{N_{\text{m}}}(\boldsymbol x)
\end{bmatrix} = \boldsymbol H \boldsymbol x
\end{eqnarray}\]
式中：\(\boldsymbol H\)是包含线性函数\(f_i(\boldsymbol x)\)系数的\(N_{\text{m}} \times N_{\text{s}}\)矩阵；\(N_{\text{m}}\) 是量测值的个数；\(N_{\text{s}}\)是待估计的未知参数的个数。
将量测值写成向量形式：
\[\begin{eqnarray} \label{7s20}
\boldsymbol z^{\rm{meas}} = \begin{bmatrix}
z_1^{\rm{meas}}\\
z_2^{\rm{meas}}\\
\vdots \\
z_{N_{\text{m}}}^{\rm{meas}}
\end{bmatrix}
\end{eqnarray}\]
则方程B可以写成，记为B*：
\[\begin{eqnarray} \label{7s21}
\underset{\boldsymbol x}{\min J(\boldsymbol x)} = [\boldsymbol z^{\rm{meas}} - \boldsymbol f(\boldsymbol x)]^{\text{T}} \boldsymbol R^{-1} [\boldsymbol z^{\rm{meas}} - \boldsymbol f(\boldsymbol x)]
\end{eqnarray}\]
式中：矩阵\(R\)是测量误差的协方差矩阵。
\[
R= \begin{bmatrix}
\sigma _1^2 & & & \\
& \sigma _2^2 & & \\
& & \ddots & \\
& & & \sigma _{N_{\text{m}}}^2
\end{bmatrix}
\]
为获得式B最小值的一般表达式，展开并用\(Hx\)替换\(f(x)\)得
\[\begin{multline}
\underset{\boldsymbol x}{\min J(\boldsymbol x)} = ({{\boldsymbol z^{\rm{meas}}})^{\text{T}} \boldsymbol R^{-1} \boldsymbol z^{\rm{meas}} - \boldsymbol x^{\text{T}} \boldsymbol H^{\text{T}} \boldsymbol R^{-1} \boldsymbol z^{\rm{meas}}}-({\boldsymbol z^{\rm{meas}}})^{\text{T}}\boldsymbol R^{-1}\boldsymbol H \boldsymbol x + \boldsymbol x^{\text{T}}\boldsymbol H^{\text{T}} \boldsymbol R^{-1} \boldsymbol H \boldsymbol x
\end{multline}\]
当\(\partial J(\boldsymbol x)/\partial x_i = 0,i = 1,\cdots ,N_{\text{s}}\) 时，\(J( \boldsymbol x)\)最小，即\(\nabla J(x)\)为零。而\(J(\boldsymbol x)\)的梯度为:
\[
\nabla J(\boldsymbol x) = -2\boldsymbol H^{\text{T}}\boldsymbol R^{-1} \boldsymbol z^{\rm{meas}} + 2\boldsymbol H^{\text{T}}\boldsymbol R^{-1}\boldsymbol H\boldsymbol x
\]
由\(\nabla J(\boldsymbol x) = 0\)得
\[\begin{eqnarray} \label{7s23}
\boldsymbol x^{\rm{est}} = (\boldsymbol H^{\text{T}} \boldsymbol R^{-1} \boldsymbol H)^{-1}\boldsymbol H^{\text{T}} \boldsymbol R^{-1} \boldsymbol z^{\rm{meas}}
\end{eqnarray}\]
注意当\(N_{\text{s}} < N_{\text{m}}\)，即当需要估计的参数的数量小于量测值的数量时，上式成立。 当\(N_{\text{s}} = N_{\text{m}}\)时，上述估计问题简化为
\[\begin{eqnarray} \label{7s24}
\boldsymbol x^{\rm{est}} = \boldsymbol H^{-1} \boldsymbol z^{\rm{meas}}
\end{eqnarray}\]
当\(N_{\text{s}} > N_{\text{m}}\)时，这个问题依然有一个封闭解，但是，在这种情况下不再估计使似然函数最大的\(\boldsymbol x\)值，因为\(N_{\text{s}} > N_{\text{m}}\)通常意味着能找到很多不同的\(\boldsymbol x^{\rm{est}}\)，使得对于所有的\(i = 1,\cdots ,N_{\text{m}}\)，\(f_i( \boldsymbol x^{\rm{est}})\)等于\(z_i^{\rm{meas}}\)。取而代之的目标是找到使\(x_i^{\rm{est}}\)的平方和最小的\(\boldsymbol x^{\rm{est}}\)，即
\[\begin{eqnarray} \label{7s25}
\underset{\boldsymbol x}{\min}\sum_{i=1}^{N_{\text{s}}} x_i^2 = \boldsymbol x^{\text{T}} \boldsymbol x
\end{eqnarray}\]
受\(\boldsymbol z^{\rm{meas}} = \boldsymbol H \boldsymbol x\)的限制，这种情况下的封闭解为
\[\begin{eqnarray} \label{7s26}
\boldsymbol x^{\rm{est}} = \boldsymbol H^{\text{T}}(\boldsymbol H\boldsymbol H^{\text{T}} )^{-1} \boldsymbol z^{\rm{meas}}
\end{eqnarray}\]
在电力系统状态估计中，欠确定问题（即\(N_{\text{s}}>N_{\text{m}}\)的情况）无法解决，如上式所示。因此，需要向量测值集合中添加伪量测量，以使得它变为完全确定或过确定的问题。
算例一
在如下图所示的三母线直流潮流模型中，有三个量测值来确定母线1、2的相位角\(\theta _1\)、 \(\theta _2\)。从上面内动的展开式中，可以知道状态\(\theta _1\) 、 \(\theta _2\)可以通过最小化残差\(J(\theta _1, \theta _2)\)来估计，其中\(J(\theta _1, \theta_2)\)是每个量测值的残差平方除以量测值分布的方差的总和。假设算例中的三个仪表的满量程值为100 MW，仪表精度为\(\pm 3\) MW（表示在大约99％的时间内，仪表测量误差不超过\(\pm 3\) MW）。
在数学上，认为测量误差服从标准差为\(\sigma\)的正态分布，如图所示。随着误差绝对值的增大，出现此误差的概率会降低。对概率密度函数在-3\(\sigma\)和+3\(\sigma\)之间积分，得到的结果约为0.99。假设仪表的精度（在本算例的情况下为\(\pm 3\) MW）等于概率密度函数的3\(\sigma\)点，则\(\sigma\)=1 MW。
【解答】
加权最小二乘估计的公式为
\[
\boldsymbol x^{\rm{est}} = (\boldsymbol H^{\text{T}} \boldsymbol R^{-1}\boldsymbol H)^{-1}\boldsymbol H^{\text{T}} \boldsymbol R^{-1} \boldsymbol z^{\rm{meas}}
\]
式中：\(\boldsymbol x^{\rm{est}}\)为状态变量的估计值；\(\boldsymbol H\)为量测系数矩阵；\(\boldsymbol R\)为权重系数，也是量测值的协方差矩阵；\(\boldsymbol z^{\rm{meas}}\)是量测值向量。
对于三母线问题，有
\[\begin{equation}
\boldsymbol x^{\rm{est}} = \begin{bmatrix}
\theta_1^{\rm{est}}\\
\theta_2^{\rm{est}}
\end{bmatrix} \label{5.27}
\end{equation}\]
为了推导\(\boldsymbol H\)矩阵，需要将量测量\(M_{12}\)、\(M_{13}\)、\(M_{32}\)用状态变量\(\theta_1\)和\(\theta_2\)表示：
\[\begin{align}
M_{12} & = f_{12} = \frac{1}{0.2}(\theta_1 - \theta_2) = 5\theta_1 - 5\theta_2 \nonumber \\
M_{13} & = f_{13} = \frac{1}{0.4}(\theta_1 - \theta_3) = 2.5\theta_1 \label{5.28} \\
M_{32} & = f_{32} = \frac{1}{0.25}(\theta_3 - \theta_2) = -4\theta_2 \nonumber
\end{align}\]
参考母线相位角\(\theta_3\)假定为零，则：
\[
\boldsymbol H = \begin{bmatrix}
5 & -5\\
2.5 & 0\\
0 & -4
\end{bmatrix}
\]
权重矩阵\(\boldsymbol R\)为
\[
\boldsymbol R = \begin{bmatrix}
\sigma_{M12}^2 & & \\
&
\sigma_{M13}^2 & \\
& &
\sigma_{M32}^2
\end{bmatrix} =\begin{bmatrix}
0.0001 & & \\
& 0.0001 & \\
& & 0.0001
\end{bmatrix}
\]
注意，由于\(\boldsymbol H\)的系数是经过归一化的，因此需要用归一化量表示\(\boldsymbol R\)和\(\boldsymbol z^{\rm{meas}}\)。
通过最小二乘法计算得到\(\theta_1\)和\(\theta_2\)的最佳估计值为
\[\begin{align}
\begin{bmatrix}
\theta_1^{\rm{est}}\\
\theta_2^{\rm{est}}
\end{bmatrix} = & \left[ \begin{bmatrix}
5 & 2.5 & 0\\
-5 & 0 & -4
\end{bmatrix} \begin{bmatrix}
0.0001 & & \\
& 0.0001 & \\
& & 0.0001
\end{bmatrix}^{-1} \begin{bmatrix}
5 & -5\\
2.5 & 0\\
0 & -4
\end{bmatrix}\right]^{-1} \nonumber \\
& \times \begin{bmatrix}
5 & 2.5 & 0\\
-5 & 0 & -4
\end{bmatrix} \begin{bmatrix}
0.0001 & & \\
& 0.0001 & \\
& & 0.0001
\end{bmatrix}^{-1} \begin{bmatrix}
0.62\\
0.06\\
0.37
\end{bmatrix} \nonumber \\
= & \begin{bmatrix}
312500 & -250000\\
-250000 & 410000
\end{bmatrix}^{-1}\begin{bmatrix}
32500\\
-45800
\end{bmatrix} \nonumber \\
= & \begin{bmatrix}
0.028571\\
-0.094286
\end{bmatrix} \nonumber
\end{align}\]
其中：
\[
\boldsymbol z^{\rm{meas}} = \begin{bmatrix}
0.62\\
0.06\\
0.37
\end{bmatrix}
\]
根据\(\theta_1\)和\(\theta_2\)的最佳估计值得到的潮流结果如图所示。计算\(J(\theta_1,\theta_2)\)的值如下：
\[\begin{align}
J(\theta_1,\theta_2) & = \frac{[z_{12}-f_{12}(\theta_1,\theta_2)]^2}{\sigma_{12}^2} + \frac{[z_{13}-f_{13}(\theta_1,\theta_2)]^2}{\sigma_{13}^2} + \frac{[z_{32}-f_{32}(\theta_1,\theta_2)]^2}{\sigma_{32}^2} \nonumber \\
& = \frac{[0.62-(5\theta_1-5\theta_2)]^2}{0.0001} + \frac{[0.06-(2.5\theta_1)]^2}{0.0001} + \frac{[0.37+(4\theta_2)]^2}{0.0001} \nonumber \\
& = 2.14 \nonumber
\end{align}\]
假设\(M_{13}\)输电线上的仪表在测量精度上比\(M_{12}\)和\(M_{32}\)上的仪表高，可以推断从\(M_{13}\)得到的量测值将比\(M_{12}\)和\(M_{32}\)的更接近线路的真实功率值。如果在仪表参数上体现了\(M_{13}\)仪表精度更高，那么也期望在状态估计计算结果中利用这一优势，以期得到更加精确的估计结果。因此，采用如下的仪表参数数据进行分析：
则权重值矩阵\(\boldsymbol R\)为
\[
\boldsymbol R = \begin{bmatrix}\sigma_{M12}^2 & & \\ & \sigma_{M13}^2 & \\ & & \sigma_{M32}^2\end{bmatrix} =\begin{bmatrix}1\times 10^{-4} & & \\ & 1\times 10^{-6} & \\ & & 1\times 10^{-4}\end{bmatrix}
\]
代入
\[\begin{eqnarray} \label{zzz}
\boldsymbol x^{\rm{est}} = (\boldsymbol H^{\text{T}} \boldsymbol R^{-1} \boldsymbol H)^{-1}\boldsymbol H^{\text{T}} \boldsymbol R^{-1} \boldsymbol z^{\rm{meas}}
\end{eqnarray}\]
解得：
\[\begin{align}
\begin{bmatrix}
\theta_1^{\rm{est}}\\
\theta_2^{\rm{est}}
\end{bmatrix} = & \left[ \begin{bmatrix}
5 & 2.5 & 0\\
-5 & 0 & -4
\end{bmatrix} \begin{bmatrix}
1\times 10^{-4} & & \\
& 1\times 10^{-6} & \\
& & 1\times 10^{-4}
\end{bmatrix}^{-1} \begin{bmatrix}
5 & -5\\
2.5 & 0\\
0 & -4
\end{bmatrix}\right]^{-1} \nonumber \\
& \times \begin{bmatrix}
5 & 2.5 & 0\\
-5 & 0 & -4
\end{bmatrix} \begin{bmatrix}
1\times 10^{-4} & & \\
& 1\times 10^{-6} & \\
& & 1\times 10^{-4}
\end{bmatrix}^{-1} \begin{bmatrix}
0.62\\
0.06\\
0.37
\end{bmatrix} \nonumber \\
= & \begin{bmatrix}
6.5 \times 10^6 & -2.5 \times 10^5\\
-2.5 \times 10^5 & 4.1 \times 10^5
\end{bmatrix}^{-1}\begin{bmatrix}
1.81 \times 10^5 \\
-0.458\times 10^5
\end{bmatrix} \nonumber \\
= & \begin{bmatrix}
0.024115\\
-0.097003
\end{bmatrix} \nonumber
\end{align}\]
由上述结果可以得到当测量仪表\(M_{13}\)精度更高时的潮流结果，如下图所示。比较1-3线路上的潮流估计值可知，将\(\sigma_{M13}\)设置为0.1 MW使得1-3线路的潮流估计值更接近6.0 MW的仪表读数（量测值）。此外，线路1-2和3-2上的潮流估计值在本例下分别比\(M_{12}\)和\(M_{32}\)的仪表读数误差更大。
三相配电网状态估计方法
传统配电网中的功率流动是单向的，相对容易预测和管理。如今以小型分布式发电、智能家居和储能系统等为代表的分布式能源技术发展迅速，配电网变得越来越复杂，需要对其进行更加有效的状态估计。接下来将用两种模型对负载和线路参数具有显著不对称性的三相配电线路进行状态估计，第一个模型考虑相间耦合，第二个模型不考虑相间耦合。
模型建立
IEEE 13节点算例
典型的小型配电网，包含不对称馈线如单相馈线、两相馈线等，同时网络中含有不对称负荷，如三相不对称负荷、两相负荷和单相负荷。
网络及量测配置如下图所示，其中馈线632-645和645-646只包含\(B\)相，馈线671-684包含\(A\)相和\(C\)相，馈线684-611只包含\(A\)相，馈线684-652只包含\(C\)相。图中配电网的三相具有明显的不对称性。
为简化起见，模型中忽略了电压调节器、并联电容器组和变压器等元件。用于分析的第一个模型考虑相间耦合，称为耦合模型，第二个忽略相间相互耦合，称为解耦模型。下表列出了该模型初始情况下的潮流计算结果，功率和电压基准值分别为1.4 MVA和2.4 kV。
算法的数学原理
配电网络的状态可以由系统中每个节点的电压幅值和相角的矢量\(\boldsymbol x\)表示，状态变量的集合\(\boldsymbol x\)和量测值的关系为
\[\boldsymbol z= \boldsymbol h(\boldsymbol x)- \boldsymbol e\]
式中：\(\boldsymbol z=\left [ z_{1},z_{2}\cdots z_{m}\right ]^{\text{T}}\)为量测值矩阵；\(\boldsymbol x=\left [ x_{1},x_{2}\cdots x_{n}\right ]^{\text{T}}\)为状态变量矩阵；\(\boldsymbol e=\left [ e_{1},e_{2}\cdots e_{m}\right ]^{\text{T}}\)为量测误差；\(\boldsymbol h=\left [ h_{1}(\boldsymbol x),h_{2}(\boldsymbol x)\cdots h_{m}(\boldsymbol x)\right ]^{\text{T}}\)，其中\(h_{i}(\boldsymbol x)\) 是非线性函数，代入状态向量\(\boldsymbol x\)可以计算得到量测值；\(m\)是独立量测值的数量；\(n\)是状态变量的数量。
将加权最小二乘状态估计法应用于目标函数：
\[J(\boldsymbol x)=\sum_{i=1}^{m}W_{ii}\cdot (z_{j}-h_{i}(\boldsymbol x))^{2}\]
式中：\(\boldsymbol W\)是测量权重矩阵，为对角阵。
三相线路和单相线路状态估计之间的主要区别在于函数集\(h_{i}\)。在三相状态估计中，它由潮流方程确定，同时需要考虑网络相位间的关系。三相状态估计中节点的有功功率和无功功率（注入）为
线路有功功率和无功功率的量测值为
式中：\(G_{ij}^{k,l}\)是导纳矩阵中节点\(i\)第\(k\)相与节点\(k\)第\(l\)相之间的电导；\(B_{ij}^{k,l}\)是导纳矩阵中节点\(i\)第\(k\)相与节点\(k\)第\(l\)相之间的电纳；\(g_{ij}^{k,l}\)是连接节点\(i\)和\(j\)的线路的第\(k\)相和第\(l\)相之间的串联电导；\(b_{ij}^{k,l}\)是连接节点\(i\)和\(j\)的线路的第\(k\)相和第\(l\)相之间的串联电纳；NB为节点数；\(i,j\)为节点编号；\(k,l\)为相编号，取决于网络拓扑结构；\(P_{i}^{k},Q_{i}^{k}\)是节点\(i\)第\(k\)相的有功功率和无功功率注入；\(P_{ij}^{k},Q_{ij}^{k}\)是从节点\(i\)传输到节点\(j\)的第\(k\)相有功功率和无功功率；\(V_{i}^{l}\)是节点\(i\)上第\(l\)相的电压幅值；\(\delta _{i}^{l}\)是节点\(i\)上第\(l\)相的电压相角。
算例测试结果
馈线模型采用分布式参数。在基础潮流计算结果上，叠加服从不同方差标准正态分布的随机误差，模拟测量仪器的误差。根据测量仪表的估计精度确定权重矩阵的对角元素，算法收敛条件设置为\(10^{-4}\)。结果表明，耦合模型的电压相角更接近真实值，而解耦模型的电压相角与真实值存在较大的误差。
以合格率最大为目标的状态估计新方法
以合格率最大为目标，考虑潮流约束和运行中上下界约束的电力系统状态估计方法。
与以往方法相比，该方法可用现代内点法求解。
优点：合格率高，估计结果有较强的抗差性，计算中无须做坏数据检验、可观性校验，无须主观设定测点权重，大幅减少了调试和维护工作量。
投票表决的评价策略
如前所述，电力系统的遥测量方程可表示为
\(Z_{i}=h_{i}(\boldsymbol x)+e_{i},\ \ i=1,2,3,\cdots ,N\)
在实际应用中，对测点\(i\)若：
\[\begin{equation}
\left | e_{i} \right |=\left | Z_{i}-h_{i}(\boldsymbol x) \right |\leq \alpha _{i}\label{12.28}
\end{equation}\]
则称测点\(i\)合格，反之则不合格。\(\alpha _{i}\)为给定常量，可按电力公司已有标准选取。
经典最小二乘状态估计方法的基本思想：在已知量测情况下，找到某个状态\(x\)使得残差平方和最小。
缺点：当某个测点量测数据为坏数据时，其残差较大，对最小二乘结果影响也较大，造成所谓“残差污染”。实际应用中，为克服这一问题，往往要通过某些方法判断出该量测数据为坏数据，计算时给予较小的权重，来减少坏数据对估计结果的影响（即提高“抗差性”）。这种方法需要主观确定加权因子，调试和维护较为困难。
采用一种全新思路：
(1)任给系统某个状态\(x\)。
(2)在这一状态下，对任一测点\(i\)，若其是合格测点，则认为该测点“投票同意”这一状态；反之，认为该测点“投票反对”这一状态。
(3)认为最多测点赞同的状态即为待求的系统状态。显然，在该状态下系统合格点数最多（即合格率最大）。
测点评价函数的数学模型
可建立如下测点评价函数：
\[\begin{equation}
g_{i}(e_{i})=\left\{\begin{matrix}
0 & \left | e_{i} \right |\leq\alpha _{i}\\
1 & \left | e_{i} \right |> \alpha _{i}
\end{matrix}\right.
\end{equation}\]
上式意义为当测点\(i\)合格（\(\left | e_{i} \right |\leq \alpha _{i}\)）时，\(g_{i}(e_{i})\)为0；反之，\(g_{i}(e_{i})\)为1。
\(g_{i}(e_{i})\)形状见下图，其缺点是在\(\pm \alpha_i\)处有间断点，并且不是处处连续可导，不便于实际应用。
为此，建立\(g_{i}(e_{i})\)的近似函数，令
\[\begin{equation}
f_{i}(e_{i})=\frac{1}{1+e^{-\frac{(e_{i}-\alpha _{i})c}{\alpha _{i}}}}+\frac{1}{1+e^{\frac{(e_{i}+\alpha _{i})c}{\alpha _{i}}}}
\end{equation}\]
式中：\(c\)一般为大于3的常数。
设系统中不合格测点数目为\(k\)，则\(\sum_{i=1}^{N}f_{i}(e_{i})\approx k\)。而要寻求某一系统状态，使得合格测点数最多，即相当于寻求某一系统状态，使得在该状态下\(\sum_{i=1}^{N}f_{i}(e_{i})\)最小。得到新的状态估计模型如下：
\[\begin{equation}
\left\{\begin{matrix}
\min_{\boldsymbol x}\sum_{i=1}^{N}f_{i}(e_{i}) & \\
s.t.\ Z_{i}=h_{i}(\boldsymbol x)+e_{i}& i=1,2,\cdots ,N
\end{matrix}\right.
\end{equation}\]
由于实际运行中，真实系统状态必然满足潮流约束和运行中上下界约束，考虑这些约束后，得到改进的状态估计模型：
\[\begin{equation}
\left\{\begin{matrix}
\min_{\boldsymbol x}\sum_{i=1}^{N}f_{i}(e_{i}) & \\
s.t.\ Z_{i}=h_{i}(\boldsymbol x)+e_{i} &i=1,2,\cdots ,N \\
g(\boldsymbol x)=0 & \\
\underline{\boldsymbol x}\leq \boldsymbol x \leq \bar{\boldsymbol x} &i=1,2,\cdots ,N
\end{matrix}\right.
\end{equation}\]
式中：\(g(\boldsymbol x)=0\)代表潮流约束；\(\underline{\boldsymbol x}\leq \boldsymbol x \leq \bar{\boldsymbol x}\)代表各测点运行中上下界约束。
将\(e_{i}\)代入上式中的目标函数即得
\[\begin{equation}
\left\{\begin{matrix}
\min_{\boldsymbol x}\sum_{i=1}^{N}f_{i}(Z_{i}-h_{i}(\boldsymbol x)) & \\
s.t.\ g(\boldsymbol x)=0 & \\
\underline{\boldsymbol x}\leq \boldsymbol x \leq \bar{\boldsymbol x} &i=1,2,\cdots ,N
\end{matrix}\right.\label{12.29}
\end{equation}\]
以合格率最高为目标的配电网状态估计
设某一三相不平衡配电系统的支路数目为\(m\)，节点数目为\(n+1\)（其中包含一个参考节点）。系统中\(\varphi\)相第\(l\)条支路电流\({\dot I_{l,\varphi }}\)及\(\varphi\)相第\(k\)个节点电压如下：
\[\begin{align}
\begin{split}
{\dot I_{l,\varphi }} = {I_{{\rm{r}}l,\varphi }} + {\rm{j}}{I_{{\rm{x}}l,\varphi }}\ \ \ \ l = 1,2 \cdots ,m; \varphi = a,b,c \\
{\dot U_{k,\varphi }} = {U_{{\rm{r}}k,\varphi }} + {\rm{j}}{U_{{\rm{x}}k,\varphi }}\ \ \ \ k = 1,2 \cdots ,n;\varphi = a,b,c
\end{split}
\end{align}\]
以三相支路电流相量的实部和虚部、节点电压相量的实部和虚部作为状态变量，则该配电系统的状态变量共\(6m+6n\)个。
配电网状态估计的目标函数，即为前面的式子：
\[\begin{equation}
\left\{\begin{matrix}
\min_{\boldsymbol x}\sum_{i=1}^{N}f_{i}(Z_{i}-h_{i}(\boldsymbol x)) & \\
s.t.\ g(\boldsymbol x)=0 & \\
\underline{\boldsymbol x}\leq \boldsymbol x \leq \bar{\boldsymbol x} &i=1,2,\cdots ,N
\end{matrix}\right.\label{mmm}
\end{equation}\]
不同量测的测量方程
实际配电系统中，需要考虑的量测包括节点电压幅值、支路电流幅值、节点注入功率（负荷功率）、支路功率。下面分别就不同量测给出相应的量测方程。
（1）支路电流幅值量测
支路电流幅值量测可由下式描述：
\[\begin{equation}
|{\dot I_{l,\varphi }}| = \sqrt {{I_{{\rm{r}}l,\varphi }}^2 + {I_{{\rm{x}}l,\varphi }}^2} {\rm{ }}\ \ \ \ l = 1,2 \cdots ,m;{\rm{ }}\varphi = a,b,c \label{7-5-31}
\end{equation}\]
由于上式中存在根号，雅可比矩阵对应电流幅值部分的表示形式过于复杂，求解困难。为了简化雅可比矩阵，将电流幅值量测转变为电流幅值平方的量测，如式()所示：
\[\begin{equation}
|{\dot I_{l,\varphi }}{|^2} = {I_{{\rm{r}}l,\varphi }}^2 + {I_{{\rm{x}}l,\varphi }}^2{\rm{ }}\ \ \ \ l = 1,2 \cdots ,m;\varphi = a,b,c \label{7-5-32}
\end{equation}\]
（2）节点电压幅值量测
与上述电流幅值量测的处理方式相同，电压幅值量测由如下量测方程进行描述：
\[\begin{equation}
|{\dot U_{k,\varphi }}{|^2} = {U_{{\rm{r}}k,\varphi }}^2 + {U_{{\rm{x}}k,\varphi }}^2{\rm{ }}\ \ \ \ k = 1,2 \cdots ,n;\varphi = a,b,c \label{7-5-33}
\end{equation}\]
（3）节点注入功率量测
规定电流的正方向是从节点编号较小的节点流向节点编号大的节点。对于编号为\(k\)的节点，假设有\(p+q\)个节点与之相连，其中有\(p\)个节点编号小于\(k\)，\(q\)个节点编号大于\(k\)。则有\(p\)条支路电流流向节点\(k\)，\(q\)条支路电流流出节点\(k\)，如图所示。
根据上图可知节点\(k\)的\(\varphi\)相注入功率为
\[\begin{equation}
{P_{k,\varphi }} + {\rm{j}}{Q_{k,\varphi }} = {\dot U_{k,\varphi }}{(\sum\limits_{l = p + 1}^{p + q} {{{\dot I}_{l,\varphi }} - } \sum\limits_{i = 1}^p {{{\dot I}_{l,\varphi }}} )^ * } \label{7-5-34}
\end{equation}\]
用状态量表示如下：
\[\begin{gather}
{P_{k,\varphi }} = \sum\limits_{l = p + 1}^{p + q} {({U_{{\rm{r}}k,\varphi }}{I_{{\rm{r}}l,\varphi }} + {U_{{\rm{x}}k,\varphi }}{I_{{\rm{x}}l,\varphi }})} - \sum\limits_{l = 1}^p {({U_{{\rm{r}}k,\varphi }}{I_{{\rm{r}}l,\varphi }} + {U_{{\rm{x}}k,\varphi }}{I_{{\rm{x}}l,\varphi }})} \label{7-5-35}\\
{Q_{k,\varphi }} = \sum\limits_{l = p + 1}^{p + q} {({U_{{\rm{x}}k,\varphi }}{I_{{\rm{r}}l,\varphi }} - {U_{{\rm{r}}k,\varphi }}{I_{{\rm{x}}l,\varphi }})} - \sum\limits_{l = 1}^p {({U_{{\rm{x}}k,\varphi }}{I_{{\rm{r}}l,\varphi }} - {U_{{\rm{r}}k,\varphi }}{I_{{\rm{x}}l,\varphi }})} \label{7-5-36}
\end{gather}\]
（4）支路功率量测
支路\(l\)（首端和末端节点分别为\(i\)和\(j\)）的首末端功率可表示如下：
\[\begin{equation}
\left\{ \begin{array}{l}
{P_{ij,\varphi }} + {\rm{j}}{Q_{ij,\varphi }} = {{\dot U}_{i,\varphi }}{{\dot I}_{l,\varphi }}^*\\
{P_{ji,\varphi }} + {\rm{j}}{Q_{ji,\varphi }} = {{\dot U}_{j,\varphi }}{( - {{\dot I}_{l,\varphi }})^*}
\end{array} \right. \label{7-5-37}
\end{equation}\]
用状态量表示如下：
\[\begin{gather}
\left\{ \begin{array}{l}
{P_{ij,\varphi }} = {U_{{\rm{r}}i,\varphi }}{I_{{\rm{r}}l,\varphi }} + {U_{{\rm{x}}i,\varphi }}{I_{{\rm{x}}l,\varphi }}\\
{Q_{ij,\varphi }} = - {U_{{\rm{r}}i,\varphi }}{I_{{\rm{x}}l,\varphi }} + {U_{{\rm{x}}i,\varphi }}{I_{{\rm{r}}l,\varphi }}
\end{array} \right. \label{7-5-38}\\
\left\{ \begin{array}{l}
{P_{ji,\varphi }} = - {U_{{\rm{r}}j,\varphi }}{I_{{\rm{r}}l,\varphi }} - {U_{{\rm{x}}j,\varphi }}{I_{{\rm{x}}l,\varphi }}\\
{Q_{ji,\varphi }} = {U_{{\rm{r}}j,\varphi }}{I_{{\rm{x}}l,\varphi }} - {U_{{\rm{x}}j,\varphi }}{I_{{\rm{r}}l,\varphi }}
\end{array} \right. \label{7-5-39}
\end{gather}\]
不同量测的约束条件
（1）电压约束
除首节点外，每个子节点\(t\)都存在父节点\(s\)，如图所示。它们的节点电压存在下式所示关系。
\[\begin{equation}
\left[ \begin{array}{l}
{{\dot U}_{s,a}}\\
{{\dot U}_{s,b}}\\
{{\dot U}_{s,c}}
\end{array} \right] = \left[ \begin{array}{l}
{{\dot U}_{t,a}}\\
{{\dot U}_{t,b}}\\
{{\dot U}_{t,c}}
\end{array} \right] + \left[ \begin{array}{l}
{Z_{aa}}{\quad}{Z_{ab}}{\quad}{Z_{ac}}\\
{Z_{ba}}{\quad}{Z_{bb}}{\quad}{Z_{bc}}\\
{Z_{ca}}{\quad}{Z_{cb}}{\quad}{Z_{cc}}
\end{array} \right]\left[ \begin{array}{l}
{{\dot I}_{l,a}}\\
{{\dot I}_{l,b}}\\
{{\dot I}_{l,c}}
\end{array} \right] \label{7-5-40}
\end{equation}\]
对\(\varphi\)相电压，用状态变量表示式()如下：
\[\begin{equation}
\left\{ \begin{array}{l}
{U_{{\rm{r}}t,\varphi }} - {U_{{\rm{r}}s,\varphi }} + \sum\limits_{\phi = a}^c {({R_{\varphi ,\phi }}{I_{{\rm{r}}l,\phi }}} - {X_{\varphi ,\phi }}{I_{{\rm{x}}l,\phi }}) = 0\\
{U_{{\rm{x}}t,\varphi }} - {U_{{\rm{x}}s,\varphi }} + \sum\limits_{\phi = a}^c {({R_{\varphi ,\phi }}{I_{{\rm{x}}l,\phi }}} + {X_{\varphi ,\phi }}{I_{{\rm{r}}l,\phi }}) = 0
\end{array} \right. \label{7-5-41}
\end{equation}\]
式中：\({{\rm{Z}}_{\varphi ,\phi }} = {R_{\varphi ,\phi }} + {\rm{j}}{X_{\varphi ,\phi }}\)。
（2）电流约束
对于配电网内部的连接节点（无功率注入），如上上图中的节点\(k\)，有\(p\)条支路流向节点\(k\)，\(q\)条支路电流流出节点\(k\)。根据基尔霍夫电流定律，流入节点\(k\)的电流等于流出该节点的电流的总和，即
\[\begin{equation}
\sum\limits_{l = 1}^p {{{\dot I}_{l,\varphi }}} - \sum\limits_{l = p + 1}^{p + q} {{{\dot I}_{l,\varphi }}} = 0,{\rm{ }}\varphi = a,b,c \label{7-5-42}
\end{equation}\]
用状态变量的形式表示如下：
\[\begin{equation}
\left\{ \begin{array}{l}
\sum\limits_{l = 1}^p {{I_{{\rm{r}}l,\varphi }}} - \sum\limits_{l = p + 1}^{p + q} {{I_{{\rm{r}}l,\varphi }}} = 0,{\rm{ }}\varphi = a,b,c\\
\sum\limits_{l = 1}^p {{I_{{\rm{x}}l,\varphi }}} - \sum\limits_{l = p + 1}^{p + q} {{I_{{\rm{x}}l,\varphi }}} = 0,{\rm{ }}\varphi = a,b,c
\end{array} \right. \label{7-5-43}
\end{equation}\]
求解算法
\[\begin{equation}
\left\{\begin{matrix}
\min_{\boldsymbol x}\sum_{i=1}^{N}f_{i}(Z_{i}-h_{i}(\boldsymbol x)) & \\
s.t.\ g(\boldsymbol x)=0 & \\
\underline{\boldsymbol x}\leq \boldsymbol x \leq \bar{\boldsymbol x} &i=1,2,\cdots ,N
\end{matrix}\right.\label{ttt}
\end{equation}\]
上式是一最优潮流问题，可用求解最优潮流的算法求解。由于现代内点法在最优潮流求解中的成功应用，最优潮流的收敛性和计算速度均得到保证，因而上式所示模型也可用现代内点法顺利求解。
引入松弛变量\({s}_{{xl}}\)、\({s}_{{xu}}\)，将不等式约束变成等式约束：
\[\begin{equation}
\left\{ \begin{array}{l}
{x} - {s}_{{xl}} - \underline{{x}} = {0}\\
{x} + {s}_{{xu}} - \bar{{x}} = {0}
\end{array} \right. ({s}_{{xl}},{s}_{{xu}} \ge {0}) \label{7-5-44}
\end{equation}\]
进一步，将目标函数改造为障碍函数，可得，记其为方程A
\[\begin{equation}
\left\{ \begin{array}{l}
\mathop {\min }\limits_x F({x}) - \mu [\sum\limits_{j = 1}^{6m + 6n} {\ln (} {s_{xl(j)}}) + \sum\limits_{j = 1}^{6m + 6n} {\ln (} {s_{xu(j)}})]\\
s.t.{\rm{ }}g({x}) = 0\\
{\rm{ }}{{x}} - {s}_{{xl}} - \underline{{x}} = 0\\
{\rm{ }}{\bar{x}} - {x} - {s}_{{xu}} = 0
\end{array} \right. \label{7-5-45}
\end{equation}\]
上式为只含等式约束的优化问题，可用拉格朗日乘子法求解，其增广拉格朗日函数为
\({y}\)、\({y}_{{xl}}\)、\({y}_{{xu}}\)为拉格朗日乘子。进一步，通过KKT条件得到：
\[\begin{equation}
\left\{ \begin{array}{l}
{\nabla _{{x}}}L = \nabla F({{x}}) - \nabla g({{x}})^{\rm{T}}{{y}} - {{{y}}_{{{xl}}}} + {{{y}}_{{{xu}}}}\\
{\nabla _{{y}}}L = - g({{x}}) = 0\\
{\nabla _{{{{y}}_{{{xl}}}}}}L = {{{s}}_{{{xl}}}} + \underline{{x}} - {{x}} = 0\\
{\nabla _{{{{y}}_{{{xu}}}}}}L = {x} + {{{s}}_{{{xu}}}} - \bar{{x}} = 0\\
{\nabla _{{{{s}}_{{{xl}}}}}}L = {{{S}}_{{{xl}}}}{{{Y}}_{{{xl}}}}{{e}} - {{\mu }} = 0\\
{\nabla _{{{{s}}_{{{xu}}}}}}L = {{{S}}_{{{xu}}}}{{{Y}}_{{{xu}}}}{{e}} - {{\mu }} = 0
\end{array} \right. \label{7-5-47}
\end{equation}\]
方程A的最优解可以通过求解上面的方程组求得，而上面的方程组则可用牛顿法求解。当松弛因子\(\mu\)趋近于零时，方程A和第一个式子有相同的最优解。
算例三
采用IEEE 13节点系统对以合格率最大为目标的配电网状态估计方法进行算例分析。量测数据用如下方法来模拟。首先用梯形迭代法计算网络潮流，然后将潮流计算结果加上2\(\%\)的高斯噪声，得到试验用的生数据。进一步通过将生数据改变符号、置零或加减量测值20\(\%\)以上等方法得到试验用的不良数据。用加权最小二乘法（WLS法）计算时，测点\(i\)的权重取\(1/\sigma _{i}^{2}\)，\(\sigma _{i}\) 为正态分布的标准差。本算例中取\(\alpha_{i}=3σ_{i}\)，\(c=3\)。
评价配电网状态估计算法主要从算法性能（收敛性和计算时间）和估计结果合理性2个方面衡量。采用文献[4]中真值已知情况下的结果评价指标\(\xi_{3\sigma}\)、\(\xi_{2\sigma}\)、\(\xi_{\sigma}\)对该方法和WLS法的估计结果合理性进行评价分析。
\(\xi_{3\sigma}\)、\(\xi_{2\sigma}\)、\(\xi_{\sigma}\)分别表示系统中估计值与真值之差绝对值在3倍标准差、2倍标准差、1倍标准差之内的测点数目与总测点数目的比值，表征不同标准差下系统中测点估计值与真值的靠近程度。指标越大，方法的估计结果合理性越高。
计算结果
下表给出了IEEE 13节点系统在坏数据比例为0、3\(\%\)、6\(\%\)时的量测配置
下表给出了本节方法和WLS法的收敛性、计算时间和结果合理性指标的比较结果。
计算结果分析
对算例结果进行分析，可以得出如下结论：
（1）当系统中不存在坏数据时，该方法和WLS方法的迭代次数、计算时间以及3个指标基本相同。
（2）当系统中存在坏数据时，该方法的迭代次数和计算时间低于WLS方法，保持很好的收敛性。随着坏数据的增多，该方法的3个指标变化较小，测点正常率保持在98\(\%\)以上，估计结果的合理性较高。WLS法的3个指标急剧减小，估计结果的合理性急剧降低。
从测试结果可以看出，该方法在系统中存在坏数据时，具有明显的优势。该方法不易受坏数据影响，具有很强的抗差性。随着坏数据增多，合理性指标仍保持较高水平，估计结果具有稳定性和合理性。同时，随着系统规模的增大和测点数目的增多，该方法的迭代次数和计算时间只是略有增加，仍保持较好的收敛性。
算法特点总结
与以往状态估计方法相比，以合格率最大为目标的状态估计方法具有以下特点：
（1）估计准确性不易受不良数据影响，具有很强的抗差性。以往状态估计方法易受不良数据影响，且估计值偏离量测值越远，受影响就越大。该方法中，当测点不合格时，无论估计值偏离量测值多远, 其反映在目标函数中大小都是1，因而估计准确性不易受不良数据影响，抗差性强。
（2）可自动对坏数据点进行识别。该方法所估计状态为绝大多数测点所赞同的状态，由此推论：若绝大多数点赞同某一状态，则投票反对该状态的测点必为坏量测点。因而该方法在进行状态估计时，自动对坏数据点进行了识别。
（3）其他特点。所求得的状态估计结果为潮流解，且满足各种物理约束，更加接近实际；计算中无须做坏数据校验、可观性校验、权重因子设置，调试和维护极为简单；求解算法为现代内点法，收敛性好、计算速度快。
参考文献
- 何正友. 配电网分析及应用[M]. 北京: 科学出版社, 2017.
- Wood A J, Wollenberg B F. Power generation, operation and control[M]. Beijing: Tsinghua University press, 2003.
- Chusovitin P, Polyakov I, Pazderin A. Three-phase state estimation model for distribution grids[C]//. 2016 International Conference on the Science of Electrical Engineering, 2016.
- 董树锋, 何光宇, 孙英云, 等. 以合格率最大为目标的电力系统状态估计新方法[J]. 电力系统自动化, 2009, 33(16): 40-43.
- 何光宇, 董树锋. 基于测量不确定度的电力系统状态估计（一）结果评价[J]. 电力系统自动化, 2009, 33(19): 21-35.
- 何光宇, 董树锋. 基于测量不确定度的电力系统状态估计（二）方法研究[J]. 电力系统自动化, 2009, 33(20): 32-36.
- 何光宇, 董树锋. 基于测量不确定度的电力系统状态估计（三）算法比较[J]. 电力系统自动化, 2009, 33(21): 28-71.
- 王雅婷, 何光宇, 董树锋. 基于测量不确定度的配电网状态估计新方法[J]. 电力系统自动化, 2010, 34(7): 40-44.
