CN105977961A - 一种基于自动微分技术的温度状态估计方法 - Google Patents

Abstract

本发明提出了一种基于自动微分(automatic differentiation,AD)技术的温度状态估计方法。在传统的状态估计计算过程中始终保持线路电阻不变，但因为输电线路电阻是会随着外界温度的变化而变化的，本发明将温度作为新的状态量引入到状态估计过程中，建立了计及温度变化的状态估计模型。另外，状态估计中需要计算量测函数的雅克比矩阵，由于考虑了支路温度变量，因而手工计算数量巨大的微分函数和编写微分代码，工作过于繁琐且容易出错，所以本发明还通过使用AD技术替代传统的手工计算雅克比矩阵，减少了人工编写代码的繁琐和出错的机率，有效避免了截断误差，提高了程序的开发效率，而且本发明易于与已有的状态估计软件相结合。

Description

一种基于自动微分技术的温度状态估计方法

技术领域

发明涉及一种基于自动微分技术的温度状态估计方法，属于电力系统运行和控制技术领域。

背景技术

作为能量管理系统(Energy Management System，EMS)的核心，电力系统状态估计通过对生数据的处理，获得状态量的最佳估计值。传统的加权最小二乘法(Weighted LeastSquares，WLS)状态估计算法估计质量和收敛性能很好，是状态估计的经典解法和理论基础，适应各种类型的量测系统。

传统的状态估计计算过程中假定输电线路电阻参数始终不变，以恒定的电网节点导纳矩阵进行计算。但是根据电热耦合的相关研究，输电线路的温度及电阻随着环境温度的不同而变化，所以需要考虑温度对电阻的影响。在常规的状态估计中都忽略了电阻的变化，状态量中也不包含温度这一变化量，导致了计算结果与实际情况间的较大误差。

基本加权最小二乘法是电力系统状态估计的最基本方法，该方法模型简单，收敛性能好，估计质量高，是目前应用最为广泛的方法之一。该方法的不足就是计算雅克比矩阵时需要手工推导微分公式和编写微分代码，工作繁琐且容易出错，效率低下。当支路温度作为新增的状态量后，雅克比元素会更加复杂，该方法的缺点也更明显。

自动微分(automatic differentiation,AD)技术的出现克服了这一缺点，它是计算机数值计算和分析领域内的一项完全崭新的技术。AD将微分定义为代数运算，与其他微分方法(如数值微分、符号微分)相比，它可以自动计算函数的任意阶导数，减少了人工编写代码的繁琐和出错的机率，从而避免了截断误差，提高了计算精度。

发明内容

发明目的：现有状态估计过程中并未考虑温度与电阻的关系，缺少相应的处理模型，并且考虑了支路温度这一状态量后，量测函数的雅克比矩阵元素会相当复杂，求导工作量巨大，而本发明提出的一种基于自动微分技术的温度状态估计方法能够解决现有技术存在的缺陷。

技术方案：本发明提出一种基于自动微分技术的温度状态估计方法，其特征在于包括以下步骤：

获取电力系统的网络参数，包括：母线编号、名称、补偿电容，输电线路的支路号、首端节点和末端节点编号、串联电阻、串联电抗、并联电导、并联电纳、 变压器变比和阻抗，系统支路的当前温度；

程序初始化，包括：对状态量设置初值、节点次序优化、形成节点导纳矩阵、设置收敛精度λ、分配内存、声明活跃变量、一个规定的参考温度t_Ref以及相关参数的基准值；

输入遥测数据z，包括：电压幅值、发电机有功功率、发电机无功功率、负荷有功功率、负荷无功功率、线路首端有功功率、线路首端无功功率、线路末端有功功率、线路末端无功功率以及支路上的温度量测；

根据电热耦合思想建立支路温度与电阻的模型，用公式表示为：

其中，R是导体的电阻，R_Ref是导体在参考温度下的电阻，t_Ref是参考温度，t_F是固定的温度。普通线路上导体的温度t是外界环境温度t_Amb与导体上升的温度t_Rise之和，可表示如下：

其中，P_Loss是设备内部所有的损耗，t_RatedRise是额定的升高温度，P_RatedLoss是对应额定的损耗；

构造计及温度变量的量测方程；

将雅克比矩阵中不变元素的位置和数值存到一个链表中；

恢复迭代计数器：k＝1；

应用AD技术计算雅克比矩阵中的可变元素，同时读取链表中相应矩阵的不变元素，获得所需的雅克比矩阵；其中含有状态量电压相角θ、电压幅值U和支路温度t的分块扩展雅克比矩阵为：

其中，H为雅克比矩阵；P、Q和L分别表示有功、无功和温度相关的自由矢量；

求解下述方程得到状态修正量Δx^(k)，选取并修正状态量得到x^(k+1)：

Δx^(k)＝[H^T(x^(k))WH(x^(k))]^-1H^T(x^(k))W[z-h(x^(k))]

x^(k+1)＝x^(k)+Δx^(k)

其中，x为状态量，z为量测量，h为非线性量测函数，H(x)为h(x)的雅克比矩阵，T表示矩阵的转置，W为对角权重矩阵，W_ii＝1/σ_i ²，σ_i为标准差；

当小于设定的收敛精度λ时，则计算结束，输出状态估计结果，否则修正状态量，计算新温度下的导纳矩阵，进行第k+1次状态估计。

有益效果：本发明提供一种基于AD技术的温度状态估计方法，通过建立计及支路温度变化的状态估计模型，不仅能够提高状态估计结果的精度，体现电阻随着温度的变化而变化，而且通过AD技术自动计算获得雅克比矩阵，有效避免了截断误差，提高了算法的计算效率和程序的开发效率，具有工程应用前景。

附图说明

图1：本发明方法流程图；

图2：本发明电力系统支路部分简化的热阻模型。

具体实施方式

下面结合附图对发明的技术流程进行详细说明：

通常金属导体的电阻与温度有如下关系：

其中，R是导体的电阻，t是导体的温度，R_Ref是导体在参考温度下的电阻，t_Ref是参考温度，t_F是固定的温度(一般铜线是234.5℃，硬铝线是228.1℃，铝变压器绕组是225.0℃)。

著名的热阻模型如图2所示，描绘了设备的热耦合现象。在热阻模型中，设备温度的上升量与设备的损耗近似成线性关系，系数为R_θ，用公式表示如下：

其中，t_Rise是设备高出外界环境的温度，P_Loss是设备内部所有的损耗，t_RatedRise是额定的升高温度，P_RatedLoss是对应额定的损耗。

普通线路上导体的温度t是外界环境温度t_Amb与导体上升的温度t_Rise之和，表示为：

其中，某条支路ij的P_Loss,ij可通过下式计算得到：

式中：g_ij表示节点i和节点j之间的电导，U_i和U_j分别表示节点i和j的电 压幅值，θ_i和θ_j分别表示节点i和j的电压相角。

对于变压器支路，其温度模型可表示为：

其中，I和I_Rated分别表示变压器支路上的电流和额定电流，n取决于变压器冷却方式：密封式变压器取0.7，自然冷式变压器取0.8，强制风冷变压器取1.0。

电力系统的实时运行和控制需要了解系统的真实运行工况，由于测量和传输等方面的原因，得到的生数据难免存在误差，状态估计能在一定程度上提高数据的精度。自1969年美国麻省理工学院的许怀丕(F.C.Schweppe)等人提出了电力系统状态估计的最基本算法——基本加权最小二乘(weighted least squares，WLS)状态估计算法以来，加权最小二乘法成为电力系统状态估计中应用最多的算法。其基本思想是以量测量和量测估计值之差的平方和最小为目标准则的估计方法。本发明提出的一种基于AD技术的温度状态估计方法，该方法是建立在非线性WLS状态估计的模型基础上，考虑了支路电阻受温度的影响，运用AD技术计算非线性WLS状态估计方法中的雅克比矩阵。

在给定网络接线、支路参数和量测系统的条件下，非线性量测方程可表示为：

z＝h(x)+ε

式中：x为状态量；z为量测值矢量即遥测数据；h为由基尔霍夫等基本电路定律所建立的量测函数；ε为量测误差。

在计及支路温度变化的电力系统状态估计中，量测量配置的类型要比常规潮流多，不仅包括各节点的注入功率量测P_i、Q_i，各节点的电压幅值量测U_i，还包括支路的功率量测P_ij、Q_ij、P_ji、Q_ji以及支路的温度量测L_ij、L_ji。考虑温度影响后，它们的量测函数分别如下所示：

节点i的电压幅值：

U_i＝U_i

节点注入的有功功率和无功功率：

式中：P_i和Q_i分别表示节点i注入的有功和无功；U_i和U_j分别表示节点i和j的电压幅值；θ_ij是节点i到节点j的电压相角差；G_ij(t)和B_ij(t)则表示温度为t时节点导纳矩阵对应节点i和j之间的实部和虚部；n是系统节点总数。

非变压器支路i-j上始端有功功率和无功功率：

式中：P_ij和Q_ij分别为线路i-j上始端有功功率和无功功率，其方向规定：由i流向j为正，由j流向i为负；g_ij(t)和b_ij(t)表示t温度时支路ij上节点i和节点j之间的电导和电纳，y_c为线路对地导纳值。

非变压器支路i-j上终端有功功率和无功功率：

式中：P_ji和Q_ji分别为线路i-j上终端有功功率和无功功率，其方向规定：由j流向i为正，由i流向j为负。

对于变压器支路，量测方程如下式所示：

变压器支路i-j上始端有功功率和无功功率：

式中：K为变压器非标准变比；j为标准侧，变比为1，i为非标准侧，变比为K；b_T(t)为变压器标准侧(j侧)的电纳。其方向规定：由i流向j为正，由j流向i为负。

变压器支路i-j上终端有功功率和无功功率：

其方向规定：由j流向i为正，由i流向j为负。

本发明基于电热耦合思想，考虑到电阻随温度的变化，需要对支路上的温度构造新的温度量测，可表示为：

支路i-j上的温度量测：

式中：L_ij表示支路ij上的温度变化量，t_ij表示支路ij上的温度，t_Amb表示外界的环境温度，R_θ,ij表示支路ij对应的热阻系数。

给定量测矢量z以后，状态估计问题就是求使目标函数J(x)＝[z-h(x)]^TW[z-h(x)]达到最小时的x的值。其中J是目标函数，T表示矩阵的转置，W为对角权重矩阵，W_ii＝1/σ_i ²，σ_i为标准差。

根据目标函数J(x)＝[z-h(x)]^TW[z-h(x)]，由极值函数可以得出：

由于h(x)和H(x)是x的非线性函数，所以无法直接计算得出故采用迭代的方法求解。令x₀是x的某一近似值，可以在x₀附近对h(x)进行泰勒展开，保留一次项，并忽略二次以上的非线性项，得到：

h(x)≈h(x₀)+H(x₀)Δx

式中Δx＝x-x₀，H(x)为h(x)的雅克比矩阵。将此式代入目标函数中，可得到：

J(x)＝[Δz-H(x₀)Δx]^TW[Δz-H(x₀)Δx]

式中Δz＝z-h(x₀)，欲使目标函数最小，将上式展开配方得到：

J(x)＝Δz^T[W-WH(x₀)∑(x₀)H^T(x₀)W]Δz

+[Δx-∑(x₀)H^T(x₀)WΔz]^TΣ^-1(x₀)[Δx-∑(x₀)H^T(x₀)WΔz]

式中∑(x₀)＝[H^T(x₀)WH(x₀)]^-1。

上式中右边第一项与Δx无关。因此，欲使J(x)极小，第二项应为0，即Δx＝∑(x₀)H^T(x₀)WΔz，因此

接下来，进行逐次迭代计算就可以求出状态量的估计值若以k表示迭代序号，则迭代计算公式如下：

Δx^(k)＝[H^T(x^(k))WH(x^(k))]^-1H^T(x^(k))W[z-h(x^(k))]

x^(k+1)＝x^(k)+Δx^(k)

因为考虑了支路温度的变化，本发明的状态估计模型在基本加权最小二乘法的基础上，引入温度t作为状态量，所以状态估计的修正量扩展到Δx＝[Δθ ΔU Δt]^T，得到新的含有t的分块扩展雅可比矩阵为：

其中，H为雅克比矩阵；P、Q和L分别表示有功、无功和温度相关的自由矢量。因为电阻会随着温度的变化而变化，所以导纳矩阵也要随着各状态量的修正而不断修正，才能获得准确的状态估计值。

在算法程序中雅克比矩阵是运用AD技术自动求得的，不需要手工推导微分公式和编写微分代码。雅克比矩阵中有些元素的值随着状态量变化而变化，但也有很大一部分元素的值均为常数零，其在迭代过程中不随状态量的变化而变化。而AD工具在生成雅克比矩阵时，处理可变元素和不变元素的过程是相同的，因此AD对雅克比矩阵中的不变元素应用链式法则求导是在重复运算，降低了程序 的计算效率。针对这个问题，本发明提出的处理方法是：对雅克比矩阵中的不变元素，在迭代前先将其各自的位置和数值存储在程序中开辟的一个链表中，在生成雅克比矩阵时直接将其从链表中读出放入对应的位置，这样就可以避免对雅克比矩阵中不变元素的重复计算；对雅克比矩阵中的可变元素，在迭代中运用AD技术自动求出，最终获得完整的雅克比矩阵，提高程序的计算效率。

上述方法具体步骤如下：

获取电力系统的网络参数，包括：母线编号、名称、补偿电容，输电线路的支路号、首端节点和末端节点编号、串联电阻、串联电抗、并联电导、并联电纳、变压器变比和阻抗，系统支路的当前温度；

程序初始化，包括：对状态量设置初值、节点次序优化、形成节点导纳矩阵、设置收敛精度λ、分配内存、声明活跃变量、一个规定的参考温度t_Ref以及相关参数的基准值；

输入遥测数据z，包括：电压幅值、发电机有功功率、发电机无功功率、负荷有功功率、负荷无功功率、线路首端有功功率、线路首端无功功率、线路末端有功功率、线路末端无功功率以及支路上的温度量测；

根据电热耦合思想建立支路温度与电阻的模型，用公式表示为：

其中，R是导体的电阻，R_Ref是导体在参考温度下的电阻，t_Ref是参考温度，t_F是固定的温度。普通线路上导体的温度t是外界环境温度t_Amb与导体上升的温度t_Rise之和，可表示如下：

其中，P_Loss是设备内部所有的损耗，t_RatedRise是额定的升高温度，P_RatedLoss是对应额定的损耗；

构造计及温度变量的量测方程；

将雅克比矩阵中不变元素的位置和数值存到一个链表中；

恢复迭代计数器：k＝1；

应用AD技术计算雅克比矩阵中的可变元素，同时读取链表中相应矩阵的不变元素，获得所需的雅克比矩阵；其中含有状态量电压相角θ、电压幅值U和支路温度t的分块扩展雅克比矩阵为：

其中，H为雅克比矩阵；P、Q和L分别表示有功、无功和温度相关的自由矢量；

求解下述方程得到状态修正量Δx^(k)，选取并修正状态量得到x^(k+1)：

Δx^(k)＝[H^T(x^(k))WH(x^(k))]^-1H^T(x^(k))W[z-h(x^(k))]

x^(k+1)＝x^(k)+Δx^(k)

其中，x为状态量，z为量测量，h为非线性量测函数，H(x)为h(x)的雅克比矩阵，T表示矩阵的转置，W为对角权重矩阵，W_ii＝1/σ_i ²，σ_i为标准差；

当小于设定的收敛精度λ时，则计算结束，输出状态估计结果，否则修正状态量，计算新温度下的导纳矩阵，进行第k+1次状态估计。

本发明提出的一种基于自动微分技术的温度状态估计方法，将温度引入到状态变量中，并为支路温度构造了新的零注入量测，最终建立了计及温度变化的状态估计模型。在求解雅克比矩阵过程中，本发明利用AD技术代替传统的手工计算，减少了人工编写代码的繁琐和出错的机率，有效避免了截断误差，提高了程序的开发效率，而且本发明易于与已有的状态估计软件相结合，具有工程应用价值。

Claims (1)

1.一种基于自动微分技术的温度状态估计方法，其特征在于包括以下步骤：

获取电力系统的网络参数，包括：母线编号、名称、补偿电容，输电线路的支路号、首端节点和末端节点编号、串联电阻、串联电抗、并联电导、并联电纳、变压器变比和阻抗，系统支路的当前温度；

程序初始化，包括：对状态量设置初值、节点次序优化、形成节点导纳矩阵、设置收敛精度λ、分配内存、声明活跃变量、一个规定的参考温度t_Ref以及相关参数的基准值；

输入遥测数据z，包括：电压幅值、发电机有功功率、发电机无功功率、负荷有功功率、负荷无功功率、线路首端有功功率、线路首端无功功率、线路末端有功功率、线路末端无功功率以及支路上的温度量测；

根据电热耦合思想建立支路温度与电阻的模型，用公式表示为：

$R=R_{\mathrm{Re}f}.\frac{t+t_F}{t_{\mathrm{Re}f}+t_F}$

其中，R是导体的电阻，R_Ref是导体在参考温度下的电阻，t_Ref是参考温度，t_F是固定的温度。普通线路上导体的温度t是外界环境温度t_Amb与导体上升的温度t_Rise之和，可表示如下：

$t=t_{Amb}+t_{Rise}=t_{Amb}+\left(\frac{P_{Loss}}{P_{RatedLoss}}\right)t_{RatedRise}$

其中，P_Loss是设备内部所有的损耗，t_RatedRise是额定的升高温度，P_RatedLoss是对应额定的损耗；

构造计及温度变量的量测方程；

将雅克比矩阵中不变元素的位置和数值存到一个链表中；

恢复迭代计数器：k＝1；

应用AD技术计算雅克比矩阵中的可变元素，同时读取链表中相应矩阵的不变元素，获得所需的雅克比矩阵；其中含有状态量电压相角θ、电压幅值U和支路温度t的分块扩展雅克比矩阵为：

$H=\left[\begin{array}{ccc}\frac{\∂P}{\∂\mathrm{\θ}}& \frac{\∂P}{\∂U}& \frac{\∂P}{\∂t}\\ \frac{\∂Q}{\∂\mathrm{\θ}}& \frac{\∂Q}{\∂U}& \frac{\∂Q}{\∂t}\\ \frac{\∂L}{\∂\mathrm{\θ}}& \frac{\∂L}{\∂U}& \frac{\∂L}{\∂t}\end{array}\right]$

其中，H为雅克比矩阵；P、Q和L分别表示有功、无功和温度相关的自由矢量；

求解下述方程得到状态修正量Δx^(k)，选取并修正状态量得到x^(k+1)：

Δx^(k)＝[H^T(x^(k))WH(x^(k))]^-1H^T(x^(k))W[z-h(x^(k))]

x^(k+1)＝x^(k)+Δx^(k)

其中，x为状态量，z为量测量，h为非线性量测函数，H(x)为h(x)的雅克比矩阵，T表示矩阵的转置，W为对角权重矩阵，σ_i为标准差；

当小于设定的收敛精度λ时，则计算结束，输出状态估计结果，否则修正状态量，计算新温度下的导纳矩阵，进行第k+1次状态估计。