本发明涉及一种对量测噪声鲁棒的数据驱动电网潮流方程线性化计算方法，属于电网潮流计算领域以及数据驱动技术领域。
背景技术：
电网潮流分析与优化是电力系统规划、运行、控制的基础，然而电网潮流方程的非线性给电网分析与优化算法带来了高计算负担以及较差收敛性的挑战。尤其是对于涉及较多不确定性分析的电网风险评估计算，以及涉及高实时性计算的电力市场节点电价计算，都对算法的计算时间与计算收敛性有较高要求。以往基于模型的电网潮流方程线性化计算方法难以保证计算精度。
近期文献Liu Y,Zhang N,Wang Y,et al.Data-Driven Power Flow Linearization:A Regression Approach[J].IEEE Transactions on Smart Grid,2018使用数据驱动方法对潮流方程进行线性化，这种方法利用电网的历史运行数据量测数据线性化潮流。由于是针对特定电网的情况“定制”的线性化潮流方程，这种方法可以大幅度提升线性化潮流方程的计算精度。
但是电网量测数据含有噪声，文献Liu Y,Zhang N,Wang Y,et al.Data-Driven Power Flow Linearization:A Regression Approach[J].IEEE Transactions on Smart Grid,2018中的方法在量测数据含有噪声时精度大幅降低。目前尚没有研究提出对量测噪声鲁棒的数据驱动电网潮流方程线性化计算方法。
除此之外，本专利所涉及的步骤中包括调用商用最优化求解软件包Mosek。Mosek是求解大规模线性规划、混合整数线性规划、二次规划、混合整数二次规划、锥优化等凸优化的商用软件包。Mosek与C、C#、Java、Python、Matlab等语言都有较好的交互。Mosek使用内点法为主要的最优化求解算法，具体原理见文献Andersen E D,Andersen K D.The MOSEK interior point optimizer for linear programming:an implementation of the homogeneous algorithm[M]//High performance optimization.Springer,Boston,MA,2000:197-232。
技术实现要素：
本发明的目的是提出一种对量测噪声鲁棒的数据驱动电网潮流方程线性化计算方法。利用电网历史量测数据对潮流方程进行线性化，对量测数据中的噪声鲁棒，以提高电网潮流方程计算的精度、减少电网运行过程中的计算量，有助于降低电网运行成本，在涉及较多不确定性分析的电网风险评估计算等领域有广阔应用场景。
本发明提出的对量测噪声鲁棒的数据驱动电网潮流方程线性化计算方法，包括以下步骤：
(1)对电网历史运行数据中的量测数据进行预处理，量测数据包括节点有功功率注入、节点无功功率注入、节点电压相角和节点电压幅值；具体步骤如下：
(1-1)将电力系统的所有N个节点划分为PQ、PV、Vθ节点，其中PQ节点代表节点有功功率注入和无功功率注入为已知量的节点，PV节点代表节点有功功率注入和电压幅值为已知量的节点，Vθ节点代表节点电压幅值和电压相角为已知量的节点，根据节点类型，将电网运行的有功功率注入、无功功率注入、电压幅值和电压相角数据分别按照PQ、PV、Vθ节点的顺序排列：
其中，P表示电网各个节点有功功率注入向量，P是N×1维向量，表示PQ节点有功功率注入向量PL的转置，表示PV节点有功功率注入向量PS的转置，表示Vθ节点有功功率注入向量PR的转置，Q表示电网各个节点无功功率注入向量，Q是N×1维向量，表示PQ节点无功功率注入向量QL的转置，表示PV节点无功功率注入向量QS的转置，表示Vθ节点无功功率注入向量的转置；V表示电网各个节点电压幅值向量，V是N×1维向量，表示PQ节点电压幅值向量VL的转置，表示PV节点电压幅值向量VS的转置，表示Vθ节点电压幅值向量VR的转置；θ表示电网各个节点电压相角向量，θ是N×1维向量，表示PQ节点电压相角向量θL的转置，表示PV节点电压相角向量θS的转置，表示Vθ节点电压相角向量θR的转置；
(1-2)将电力系统不同时间点的历史量测数据以步骤1-1中的格式为单位，整理成X矩阵与Y矩阵的形式：
其中，M为历史量测数据的时间点个数，上标1...m...M表示历史量测数据的时间点，其中xm和ym都为2N×1维向量，X和Y都是2N×M维矩阵；
(2)求解线性化的潮流方程参数，具体步骤如下：
(2-1)采用内点法求解以下(3)中的最优化问题，得到潮流方程参数初值：
在上述最优化问题中，(3a)是目标函数，(3b)(3c)是约束条件，其中，||·||F表示矩阵的Frobenius范数，即对给定矩阵D，tr(·)表示矩阵的迹，迹是矩阵的对角元素加和；
在上述最优化问题中，E、F、和C是待求未知变量，其中，E是N×N对称矩阵，即矩阵中元素与转置矩阵ET中相同位置的元素相等：E＝ET，F是N×N对称矩阵，即矩阵中元素与转置矩阵FT中相同位置的元素相等：F＝FT，是N×N对角矩阵：是N×N对角矩阵：C是N×1维向量，代表以为对角元素构建的对角矩阵，代表以为对角元素构建的对角矩阵，代表取矩阵对角元构建的向量，代表取矩阵对角元构建的向量；
在上述最优化问题中，X、Y、Pmin、Pmax、Qmin和Qmax是已知量，其中，Pmin和Pmax分别是历史数据中每个节点从时刻1到T的最小和最大有功功率注入，Qmin和Qmax分别是历史数据中每个节点从时刻1到T的最小和最大无功功率注入，Pmin、Pmax、Qmin和Qmax都是N×1维向量；
求解上述最优化问题，并利用公式(3b)计算得到2N×2N维矩阵AJ'；
(2-2)对上述潮流方程参数进行粗调，对上述步骤(2-1)的2N×2N维矩阵AJ'进行更新，包括以下步骤：
(2-2-1)采用内点法求解以下(4)中的最优化问题：
Y-εY＝AJ'X+dAJ'X-AJ'εX+C1T (4b)
在最优化问题(4)中，(4a)是目标函数，(4b)-(4d)是约束条件，其中，表示关于矩阵Σ-1的范数，即对给定矩阵D，Σ-1是Σ的逆矩阵，Σ为设定的对角阵，对角元为量测误差的方差：和分别表示εX和εY的方差；
在最优化问题(4)中，εX、εY、E、F、和C是待求未知变量，其中，εX和εY都是2N×M维矩阵，εX和εY分别表示量测矩阵X和Y的量测噪声，dAJ'是2N×2N维矩阵，表示矩阵AJ'的变化量；
在最优化问题(4)中，AJ'、Σ、X、Y、Pmin、Pmax、Qmin和Qmax是已知量，其中，AJ'由上述步骤(2-1)计算得到，Σ为设定的对角阵，对角元为量测误差的方差：和分别表示εX和εY的方差，由电网量测设备的说明上获取；
求解最优化问题(4)，利用公式(4c)计算dAJ'，以该计算得到的dAJ'更新矩阵AJ'：AJ'←AJ'+dAJ'；
(2-2-2)对不等式进行判断，若则执行步骤(2-3)；若则执行步骤(2-2-1)，其中，W为权重矩阵：σX和σY分别表示矩阵X和矩阵Y的标准差；表示Hadamard算子，计算两个矩阵对应位置元素的乘积，max|·|表示求矩阵中元素绝对值的最大值，α表示停止迭代的阈值参数，取α＝0.1；
(2-3)对上述潮流方程参数进行细调，对上述步骤(2-2)的2N×2N维矩阵AJ'进行更新，包括以下步骤：
(2-3-1)采用内点法求解(5)中的最优化问题：
Y-εY＝AX+dAX-AεX+C1T (5b)
在最优化问题(5)中，(5a)是目标函数，(5b)是约束条件，其中，||·||1表示矩阵的1-范数，即对给定矩阵D，dij表示矩阵D第i行j列的元素。
在最优化问题(5)中，εX、εY、dA和C是待求未知变量，其中，dA是2N×2N维矩阵，表示矩阵A的变化量；
在最优化问题(5)中，A、Σ、X和Y是已知量，其中，A是上述步骤(2-2)中的AJ'；
求解最优化问题(5)后，更新矩阵A：A←A+dA，并记录变量C；
(2-3-2)对不等式进行判断：若则执行步骤(3-1)，若执行步骤(2-3-1)；
(3)根据上述求解得到的潮流方程参数，得到电网的线性化潮流方程，并进行电网潮流计算，具体步骤如下：
(3-1)根据上述步骤(2)得到的参数A与C，以及上述步骤(1)中划分的节点类型PQ、PV和Vθ，写出电网的线性化潮流方程：
其中，Ci表示上述步骤(2-3-1)向量C按照(6)排列后的子向量，Aij表示上述步骤(2-3-1)矩阵A按照(6)排列后的子矩阵，在计算潮流时，对于式(6)映射左侧的因变量，PR，QS和QR是未知量，PL，PS，和QL是已知量，对于式(6)映射右侧的自变量，θL，θS和VL是未知量，θR，VS和VR是已知量，进而将矩阵(6)写成以下分块矩阵的形式：
其中，x2＝[θR,VS,VR]T和y1＝[PL,PS,QL]T是已知量，x1＝[θL,θS,VL]T和y2＝[PR,QS,QR]T是未知量，和分别表示式(6)中Aij矩阵的对应部分：
(3-2)根据下式求解上述步骤(3-1)的线性化潮流方程，得到电网线性化潮流方程的实时解：
其中x1＝[θL,θS,VL]T和y2＝[PR,QS,QR]T即为线性化的潮流方程的解。
本发明提出的对量测噪声鲁棒的数据驱动电网潮流方程线性化计算方法，其优点是：
1、考虑了电力系统量测数据中含有噪声这一特点，本发明对电力系统量测噪声鲁棒，相比于之前研究提出的数据驱动潮流方程线性化计算方法，本发明充分考虑了电力系统中量测数据的实际情况。
2、本发明方法的训练数据是量测的历史数据，反映了特定电力系统的真实运行状态，另外本发明对电力系统中的量测数据噪声比较鲁棒，因此线性化后的方程具有更高的计算精度。
3、本发明方法得到的是线性化的电网潮流计算方程，在直接计算以及作为优化问题的约束进行计算的问题中，都能够有效减少计算量，尤其是对于涉及较多不确定性分析的电网风险评估计算，以及涉及高实时性计算的电力市场节点电价计算，具有较为明显的优势。
附图说明
图1为本发明实施例中IEEE-33系统300组测试结果电压相角与幅值示意图。
具体实施方式
本发明提出的对量测噪声鲁棒的数据驱动电网潮流方程线性化计算方法，包括以下步骤：
(1)对电网历史运行数据中的量测数据进行预处理，量测数据包括节点有功功率注入、节点无功功率注入、节点电压相角和节点电压幅值；具体步骤如下：
(1-1)将电力系统的所有N个节点划分为PQ、PV、Vθ节点，其中PQ节点代表节点有功功率注入和无功功率注入为已知量的节点，PV节点代表节点有功功率注入和电压幅值为已知量的节点，Vθ节点代表节点电压幅值和电压相角为已知量的节点，根据节点类型，将电网运行的有功功率注入、无功功率注入、电压幅值和电压相角数据分别按照PQ、PV、Vθ节点的顺序排列：
其中，P表示电网各个节点有功功率注入向量，P是N×1维向量，表示PQ节点有功功率注入向量PL的转置，表示PV节点有功功率注入向量PS的转置，表示Vθ节点有功功率注入向量PR的转置，Q表示电网各个节点无功功率注入向量，Q是N×1维向量，表示PQ节点无功功率注入向量QL的转置，表示PV节点无功功率注入向量QS的转置，表示Vθ节点无功功率注入向量的转置；V表示电网各个节点电压幅值向量，V是N×1维向量，表示PQ节点电压幅值向量VL的转置，表示PV节点电压幅值向量VS的转置，表示Vθ节点电压幅值向量VR的转置；θ表示电网各个节点电压相角向量，θ是N×1维向量，表示PQ节点电压相角向量θL的转置，表示PV节点电压相角向量θS的转置，表示Vθ节点电压相角向量θR的转置；
(1-2)将电力系统不同时间点的历史量测数据以步骤1-1中的格式为单位，整理成X矩阵与Y矩阵的形式：
其中，M为历史量测数据的时间点个数，上标1...m...M表示历史量测数据的时间点，其中xm和ym都为2N×1维向量，X和Y都是2N×M维矩阵；需要指出，步骤1中记录和整理的所有数据都是含有噪声的量测数据，而不是真实数据。
(2)求解线性化的潮流方程参数，具体步骤如下：
(2-1)采用内点法求解以下(3)中的最优化问题，得到潮流方程参数初值：
在上述最优化问题中，(3a)是目标函数，(3b)(3c)是约束条件，其中，||·||F表示矩阵的Frobenius范数，即对给定矩阵D，tr(·)表示矩阵的迹，迹是矩阵的对角元素加和；
在上述最优化问题中，E、F、和C是待求未知变量，其中，E是N×N对称矩阵，即矩阵中元素与转置矩阵ET中相同位置的元素相等：E＝ET，F是N×N对称矩阵，即矩阵中元素与转置矩阵FT中相同位置的元素相等：F＝FT，是N×N对角矩阵：是N×N对角矩阵：C是N×1维向量，代表以为对角元素构建的对角矩阵，代表以为对角元素构建的对角矩阵，代表取矩阵对角元构建的向量，代表取矩阵对角元构建的向量；
在上述最优化问题中，X、Y、Pmin、Pmax、Qmin和Qmax是已知量，其中，Pmin和Pmax分别是历史数据中每个节点从时刻1到T的最小和最大有功功率注入，Qmin和Qmax分别是历史数据中每个节点从时刻1到T的最小和最大无功功率注入，Pmin、Pmax、Qmin和Qmax都是N×1维向量；
求解上述最优化问题，并利用公式(3b)计算得到2N×2N维矩阵AJ'；
(2-2)对上述潮流方程参数进行粗调，对上述步骤(2-1)的2N×2N维矩阵AJ'进行更新，包括以下步骤：
(2-2-1)采用内点法求解以下(4)中的最优化问题：
Y-εY＝AJ'X+dAJ'X-AJ'εX+C1T (4b)
在最优化问题(4)中，(4a)是目标函数，(4b)-(4d)是约束条件，其中，表示关于矩阵Σ-1的范数，即对给定矩阵D，Σ-1是Σ的逆矩阵，Σ为设定的对角阵，对角元为量测误差的方差：和分别表示εX和εY的方差；
在最优化问题(4)中，εX、εY、E、F、和C是待求未知变量，其中，εX和εY都是2N×M维矩阵，εX和εY分别表示量测矩阵X和Y的量测噪声，dAJ'是2N×2N维矩阵，表示矩阵AJ'的变化量；
在最优化问题(4)中，AJ'、Σ、X、Y、Pmin、Pmax、Qmin和Qmax是已知量，其中，AJ'由上述步骤(2-1)计算得到，Σ为设定的对角阵，对角元为量测误差的方差：和分别表示εX和εY的方差，由电网量测设备的说明上获取；
求解最优化问题(4)，利用公式(4c)计算dAJ'，以该计算得到的dAJ'更新矩阵AJ'：AJ'←AJ'+dAJ'；
(2-2-2)对不等式进行判断，若则执行步骤(2-3)；若则执行步骤(2-2-1)，其中，W为权重矩阵：σX和σY分别表示矩阵X和矩阵Y的标准差；表示Hadamard算子，计算两个矩阵对应位置元素的乘积，max|·|表示求矩阵中元素绝对值的最大值，α表示停止迭代的阈值参数，取α＝0.1；
(2-3)对上述潮流方程参数进行细调，对上述步骤(2-2)的2N×2N维矩阵AJ'进行更新，包括以下步骤：
(2-3-1)采用内点法求解(5)中的最优化问题：
Y-εY＝AX+dAX-AεX+C1T (5b)
在最优化问题(5)中，(5a)是目标函数，(5b)是约束条件，其中，||·||1表示矩阵的1-范数，即对给定矩阵D，dij表示矩阵D第i行j列的元素。
在最优化问题(5)中，εX、εY、dA和C是待求未知变量，其中，dA是2N×2N维矩阵，表示矩阵A的变化量；
在最优化问题(5)中，A、Σ、X和Y是已知量，其中，A是上述步骤(2-2)中的AJ'；
求解最优化问题(5)后，更新矩阵A：A←A+dA，并记录变量C；
(2-3-2)对不等式进行判断：若则执行步骤(3-1)，若执行步骤(2-3-1)；
(3)根据上述求解得到的潮流方程参数，得到电网的线性化潮流方程，并进行电网潮流计算，具体步骤如下：
(3-1)根据上述步骤(2)得到的参数A与C，以及上述步骤(1)中划分的节点类型PQ、PV和Vθ，写出电网的线性化潮流方程：
其中，Ci表示上述步骤(2-3-1)向量C按照(6)排列后的子向量，Aij表示上述步骤(2-3-1)矩阵A按照(6)排列后的子矩阵，在计算潮流时，对于式(6)映射左侧的因变量，PR，QS和QR是未知量，PL，PS，和QL是已知量，对于式(6)映射右侧的自变量，θL，θS和VL是未知量，θR，VS和VR是已知量，进而将矩阵(6)写成以下分块矩阵的形式：
其中，x2＝[θR,VS,VR]T和y1＝[PL,PS,QL]T是已知量，x1＝[θL,θS,VL]T和y2＝[PR,QS,QR]T是未知量，和分别表示式(6)中Aij矩阵的对应部分：
(3-2)根据下式求解上述步骤(3-1)的线性化潮流方程，得到电网线性化潮流方程的实时解：
其中x1＝[θL,θS,VL]T和y2＝[PR,QS,QR]T即为线性化的潮流方程的解。需要指出，步骤(1)以及步骤(2)中所涉及的电网量测数据都是历史数据，步骤(3)中的已知数据是新的电网数据，新的电网数据可以是电网实时数据、设想情况下生成的电网数据等等，根据实际应用需求而定。
以下结合附图介绍本发明方法的实施例：
本发明以Baran M E,Wu F F.Network reconfiguration in distribution systems for loss reduction and load balancing[J].IEEE Transactions on Power delivery,1989,4(2):1401-1407.中的IEEE-33测试系统为例对本发明所提出的方法进行视力验证。其中量测数据由蒙特卡洛采样与量测噪声的加和构成。蒙特卡洛采样根据利用Matpower软件包在特定初始条件下计算得到。其中有功功率负荷数据根据测试系统中的基值与随机生成的乘子想成而得，乘子数值介于[0.8-1.2]之间；无功功率负荷数据根据有功功率负荷值与随机生成的乘子相乘而得，乘子数值介于[0.15-0.25]之间。量测噪声按照信噪比(Signal to noise ratio,SNR)40dB生成，噪声服从高斯分布。信噪比的定义服从下式：
其中asignal/anoise分别表示信号/噪声的均方根幅值。本实例共取了300组数据用于训练，300组数据用于测试。
根据本发明提出的方法得到的线性化电网潮流方程计算潮流，并与精确的交流潮流计算结果对比，其结果如图1所示。图1子图(a1)与(b1)描述了电压相角的计算误差，子图(a2)与(b2)描述了电压幅值的计算误差。子图(a1)与(a2)展示了300组测试结果中的1组测试结果，横轴表示节点号，纵轴右侧表示误差，左侧表示具体计算数值。子图(b1)与(b2)展示了300组结果平均绝对误差的直方图，其中横轴表示误差，纵轴表示频数。在所有子图中，计算误差都用对数坐标轴表示。为展示本发明的有益效果，实例中使用如下方法进行对比：
1、Noiseless：利用没有噪声的量测数据，使用文献Liu Y,Zhang N,Wang Y,et al.Data-Driven Power Flow Linearization:A Regression Approach[J].IEEE Transactions on Smart Grid,2018中的方法得到线性化的潮流方程。这种方法可以看作是数据驱动方法在含有噪声条件下的“理论极限”。
2、M1：利用有含有噪声的量测数据，使用本发明的方法，得到线性化的潮流方程。
3、M2：利用文献J.Yang,N.Zhang,C.Kang,and Q.Xia,"A State-Independent Linear Power Flow Model with Accurate Estimation of Voltage Magnitude,"IEEE Trans.Power Syst.,vol.22,pp.3607-3617,2017中的方法得到线性化的潮流方程，这种方法是传统的基于模型的潮流方程线性化计算方法的代表。
4、M3：利用含有噪声的量测数据，使用文献Liu Y,Zhang N,Wang Y,et al.Data-Driven Power Flow Linearization:A Regression Approach[J].IEEE Transactions on Smart Grid,2018中的方法得到线性化的潮流方程。
从图1中可以得出以下结论：
1、在40dB量测噪声情况下，方法M3的误差比方法Noiseless的误差大三个数量级，这说明传统的数据驱动方法如果不经过对噪声鲁棒的特殊设计，其线性化潮流计算结果会在一定噪声情况下产生很大误差。
2、方法M3不仅精度差，而且其不同节点之间的计算误差相比其他方法会有明显的波动，这意味着计算结果相邻结果的差值会产生更大误差，而相邻节点差值计算的精准对潮流分析非常重要。
3、方法M1的误差比方法M2和M3的误差小两个数量级，方法M1的误差仅比方法Noiseless大不到一个数量级。这说明本发明提出的方法在含有量测噪声的情况下，显著地提升了电网线性化潮流的计算精度，且线性化精度接近没有噪声的“理论极限”。
4、方法M1的误差比方法Noiseless和M3的误差具有更小的方差，这说明方法M3在不同组计算结果中表现非常稳定。