本发明属于电网频率控制
技术领域：
，涉及一种发电机调速器频率反馈系数的优化方法。
背景技术：
：调速器是电网中发电节点一次调频的主要控制器，是根据发电节点的频率对输入机械功率进行反馈调节的手段，对电网频率稳定性具有重要作用，频率反馈控制系数对调速器性能和成本都有重要影响。随着电力系统同步相角测量单元(phasormeasurementunit，pmu)和广域监测系统(wideareameasurementsystem，wams)等技术的发展，使节点间(区域间)可以实时通信，以实现利用其它节点信息的区域协同控制。根据通信结构，控制方式分为分散式控制、区域式控制和集中式控制共3种。传统的分散式控制是根据自身节点信息进行功率控制，无需通信线路，成本代价小，但频率调节性能较低；集中式控制利用全网所有发电节点信息进行反馈控制，需要电网中的发电节点两两之间均有通信线路，频率调节性能高，但成本代价较大；区域式控制利用部分发电节点信息进行反馈控制，只需要电网中部分发电节点之间有通信线路，兼顾了分散式控制和集中式控制的优点。然而，确定利用哪个节点信息进行反馈控制和反馈系数大小以平衡性能和代价这对矛盾仍是需要解决的问题。技术实现要素：本发明提供了一种发电机调速器频率反馈系数的优化方法，解决现有技术难以确定调速器的频率反馈系数矩阵，难以平衡控制性能和通信代价的问题。本发明采用的技术方案是，一种发电机调速器频率反馈系数的优化方法，按照以下步骤实施：步骤1，建立含有调速器的电网节点动力学模型，含有调速器的电网节点动力学模型利用传统同步发电机转子的摇摆方程来描述，表达式为下式(1)：其中，下标i,j＝1,2,…,n为节点编号,n是电网中发电节点个数；ji,di分别是节点i的惯性常量和阻尼系数，θi,ωi分别是在同步旋转参考坐标系下节点i的相位和角频率，表示节点i相位的一阶导数，即节点i的角频率ωi；θj是节点j在同步旋转参考坐标系下的相位；当电力系统额定频率为f*＝50hz，则额定角频率为ω*＝2πf*；为标称机械输入功率，表示取导纳虚部，即不考虑引起线路损耗的电阻部分，lij是电网的拉普拉斯矩阵l第i行第j列的元素；ei,ej是节点i,j的端电压的幅值，yij∈y，矩阵y是电网经过kron化简后仅保留发电节点的导纳矩阵；ui是节点i的控制量，βij为节点i与节点j之间的频率反馈系数；将式(1)改写为矩阵结构的表达式如下式(2)：其中，n×1维的列向量θ＝[θ1,θ2,…,θn]t为相位向量，n×1维的列向量ω＝[ω1,ω2,…,ωn]t为角频率向量；n×1维的列向量u＝[u1,u2,…,un]t为控制向量；n×n维矩阵j为对角矩阵，第i行对角线元素为节点i的惯性常量ji；n×n维矩阵d为对角矩阵，第i行对角线元素为节点i的阻尼系数di；n×n维矩阵k为频率反馈系数矩阵，第i行第j列元素为βij；将实际电网的参数代入式(2)，得到含有调速器的电网节点动力学模型；步骤2，初始化多目标粒子群算法的参数，设置惯性权重ω、加速常数c1,c2、最大迭代次数kmax，k为粒子群算法的迭代次数，令k＝1，粒子的总数为m个，每个粒子的位置和速度均是长度为n2的行向量，n为电网中发电节点个数，将第m个粒子在第k次迭代时的位置向量和速度向量分别记为xm(k)和vm(k)，初始化第1代时每个粒子的位置向量和速度向量中的每个元素为0-1均匀分布的随机数，设置可行解精英集最大容量nfeamax，即最多存储nfeamax个粒子的信息，设置初始惯性权重ω0，以及最后一次迭代时的惯性权重ωend；初始化时，即迭代次数k＝1时，第m个粒子的位置向量和速度向量分别为xm(1)和vm(1)；xm(1),vm(1)为0-1均匀分布的随机向量；每个体向导粒子的位置向量为该粒子的初始粒子的位置向量；步骤3，计算粒子群中每个粒子的性能目标函数和稀疏性目标函数，每个粒子的位置向量xm(k)都代表一个解，即一个频率反馈系数矩阵km，位置向量的第1到第n个元素对应反馈系数矩阵km的第1行，第n+1到第2n个元素对应系数矩阵km的第2行，依次类推；将每个粒子的位置向量对应的频率反馈系数矩阵km代入式(2)中进行时域仿真，采用式(3)计算第k代每个粒子的性能目标函数j1(x)：式(3)中，ωm(t)代表采用反馈系数矩阵km时第t秒的角频率向量；um(t)代表采用反馈系数矩阵km时第t秒的控制向量；矩阵q和矩阵r为单位对角矩阵；采用式(4)计算第k代每个粒子的稀疏性目标函数j2(x)：j2(xm(k))＝||km||0(4)步骤4，更新一般精英集和可行解精英集，若对于粒子i在第k代时的两个目标函数j1(xi(k))和j2(xi(k))，如果满足j1(xi(k))＜j1(xj(k))且j2(xi(k))＜j2(xj(k))，那么称粒子i的解xi(k)为非支配解；所有非支配解构成的集合称为一般精英集，也称为pareto最优解集；称pareto最优解集中满足稳定性约束条件的解为可行pareto最优解，反之称为不可行pareto最优解；定义一个阈值jthreshold，当j1≤jthreshold时，j1所对应的pareto最优解是可行pareto最优解；反之，当j1＞jthreshold时，j1所对应的pareto最优解是不可行pareto最优解；设第k次迭代产生了nfea(k)可行pareto最优解，如果此时可行解精英集中已经存储了n(k)个历史可行解，若nfea(k)+n(k)＜nfeamax，则将nfea(k)个可行pareto最优解对应的粒子的信息全部存储到可行解精英集中；否则，在nfea(k)+n(k)个可行pareto最优解采用欧式距离法选出距离全局向导粒子距离最小的前nfeamax个可行解存入可行解精英集；步骤5，更新个体向导粒子的位置向量和全局向导粒子的位置向量xgbest，根据粒子m第k代的两个目标函数更新粒子m的个体向导粒子的目标函数为jm1min和jm2min，具体过程是：若j1(xm(k))＜jm1min且j2(xm(k))＜jm2min，则令jm1min＝j1(xm(k))，jm2min＝j2(xm(k))且否则，jm1min、jm2min和保持不变；采用欧式距离法在一般精英集中确定全局向导粒子xgbest；步骤6，更新每个粒子的位置和速度，速度向量和位置向量的更新表达式为式(5)和式(6)：其中，xm(k)为第k代时粒子m的位置向量，vm(k)为第k代时粒子m的速度向量，m＝1,2,…,m，m是粒子总数，粒子m的位置向量的第p个元素，是粒子m的速度向量的第p个元素，个体向导粒子的位置向量和全局向导粒子的位置向量xgbest由步骤5产生，是个体向导粒子的位置向量的第p个元素，xp,gbest是全局向导粒子的位置向量xgbest的第p个元素；c1,c2是加速常数，rm,1,rm,2是对应粒子m的两个随机数向量，是向量rm,1第p个元素，为0-1均匀分布的随机数，是向量rm,2第p个元素，也为0-1均匀分布的随机数，p＝1,2,…,n2，n是频率反馈系数矩阵k的维数；ω(k)为第k代的权重系数，按照式(7)所示的线性递减权值确定ω(k)的大小：ω(k)＝(ω0-ωend)(kmax-k)/kmax+ωend(7)步骤7，判断迭代次数是否到达初始化所设置的最大迭代次数kmax，若没有达到最大迭代次数kmax，则返回步骤3；否则，结束循环，以可行解精英集为最终的优化解集，在可行解精英集中寻找稀疏性目标函数j2(xm(kmax))最小对应的粒子，设该粒子的位置向量为xm(kmax)，则位置向量xm(kmax)的第1到第n个元素对应反馈系数矩阵km的第1行，第n+1到第2n个元素对应反馈系数矩阵km的第2行，依次类推，得到反馈系数矩阵km，即为最终优化结果。本发明的有益效果是，该方法基于多目标粒子群的稀疏提升算法，在降低成本代价的同时能够提升电网频率响应性能。附图说明图1是本发明方法采用欧式距离法确定全局向导粒子的示意图；图2是本发明方法实施例采用的ieee-39标准测试系统示意图；图3是对比方法针对ieee-39标准测试系统得到频率反馈系数矩阵的非零元素示意图；图4是本发明方法针对的ieee-39标准测试系统得到频率反馈系数矩阵的非零元素示意图。具体实施方式下面基于ieee-39标准测试系统，对本发明调速器频率反馈系数矩阵的优化多对目标粒子群的稀疏提升优化过程进行具体说明。本发明方法，按照以下步骤实施：步骤1，建立含有调速器的电网节点动力学模型，含有调速器的电网节点动力学模型利用传统同步发电机转子的摇摆方程来描述，表达式为下式(1)：其中，下标i,j＝1,2,…,n为节点编号，n是电网中发电节点个数；ji,di分别是节点i的惯性常量和阻尼系数，θi,ωi分别是在同步旋转参考坐标系下节点i的相位和角频率，表示节点i相位的一阶导数，即节点i的角频率ωi；θj是节点j在同步旋转参考坐标系下的相位；一般情况下，电力系统额定频率为f*＝50hz，则额定角频率为ω*＝2πf*；为标称机械输入功率，表示取导纳虚部，即不考虑引起线路损耗的电阻部分，lij是电网的拉普拉斯矩阵l第i行第j列的元素；ei,ej是节点i,j的端电压的幅值，yij∈y，矩阵y是电网经过kron化简后仅保留发电节点的导纳矩阵；ui是节点i的控制量，βij为节点i与节点j之间的频率反馈系数；此处假设：1)传输线路损耗忽略不计；2)节点的端电压ei,ej工作在其额定工作点，为常数；将式(1)改写为矩阵结构的表达式如下式(2)：其中，n×1维的列向量θ＝[θ1,θ2,…,θn]t为相位向量，n×1维的列向量ω＝[ω1,ω2,…,ωn]t为角频率向量；n×1维的列向量u＝[u1,u2,…,un]t为控制向量；n×n维矩阵j为对角矩阵，第i行对角线元素为节点i的惯性常量ji；n×n维矩阵d为对角矩阵，第i行对角线元素为节点i的阻尼系数di；n×n维矩阵k为频率反馈系数矩阵，第i行第j列元素为βij；参照图2，实施例采用ieee-39标准测试系统，其中共包含10个发电机，分别位于母线30至母线39上，令这10个发电机所在的母线依次为发电节点1至发电节点10，则节点数为n＝10；在电网节点动力学模型中，节点i的惯性常量ji、阻尼系数di、标称机械输入功率初始相位θi和初始角频率ωi具体数值，如表1所示。表1、ieee-39标准测试系统含调速器的电网节点动力学模型参数将实际电网的参数代入式(2)，得到含有调速器的电网节点动力学模型。步骤2，初始化多目标粒子群算法的参数，设置惯性权重ω、加速常数c1,c2、最大迭代次数kmax，k为粒子群算法的迭代次数，令k＝1，粒子的总数为m个，每个粒子的位置和速度均是长度为n2的行向量，n为电网中发电节点个数，将第m个粒子在第k次迭代时的位置向量和速度向量分别记为xm(k)和vm(k)，初始化第1代时每个粒子的位置向量和速度向量中的每个元素为0-1均匀分布的随机数，设置可行解精英集最大容量nfeamax，即最多存储nfeamax个粒子的信息，设置初始惯性权重ω0，以及最后一次迭代时的惯性权重ωend；本实施例多目标粒子群算法的具体参数值见表2，以下按照3个粒子m＝3为例对多目标粒子群算法进行描述。表2多目标粒子群算法的具体参数初始化时，即迭代次数k＝1时，第1个粒子的位置向量和速度向量分别为x1(1)和v1(1)；第2个粒子的位置向量和速度向量分别为x2(1)和v2(1)；第3个粒子的位置向量和速度向量分别为x3(1)和v3(1)；将3个粒子的位置向量x1(1),x2(1),x3(1)设为0-1均匀分布的随机向量，将3个粒子的速度向量v1(1),v2(1),v3(1)设为0-1均匀分布的随机向量；因为节点数n＝10，所以3个粒子的位置向量和速度向量均为n2×1维，即100×1维。步骤3，计算粒子群中每个粒子的性能目标函数和稀疏性目标函数，每个粒子的位置向量xm(k)都代表一个解，即一个频率反馈系数矩阵km，位置向量的第1到第n个元素对应反馈系数矩阵km的第1行，第n+1到第2n个元素对应系数矩阵km的第2行，依次类推；将每个粒子的位置向量对应的频率反馈系数矩阵km代入式(2)中进行时域仿真，采用式(3)计算第k代(本代)每个粒子的性能目标函数j1(x)：式(3)中，ωm(t)代表采用反馈系数矩阵km时第t秒的角频率向量；um(t)代表采用反馈系数矩阵km时第t秒的控制向量；矩阵q和矩阵r为单位对角矩阵；采用式(4)计算第k代(本代)每个粒子对应的稀疏性目标函数j2(x)：j2(xm(k))＝||km||0(4)本实施例中将每个粒子的位置向量转换为频率反馈系数矩阵k的形式，即第1个粒子的位置向量的第1个元素到第10个元素为矩阵k的第1行，第11个元素到第20个元素为矩阵k的第2行，以此类推；第2个粒子和第3个粒子采用同样的转换方式得到矩阵k；将3个粒子所得的3个频率反馈系数矩阵k代入(2)进行仿真的结果分别代入式(4)，按照表1设置的参数值进行30s的时域仿真，用式(3)计算3个粒子的性能目标函数j1分别为j1(x1(1))＝25.67，j1(x2(1))＝29.33，j1(x3(1))＝26.51。将3个粒子所得的3个频率反馈系数矩阵k分别代入式(4)，计算3个粒子的稀疏性目标函数j2分别为j2(x1(1))＝100，j2(x2(1))＝98，j2(x3(1))＝98。步骤4，更新一般精英集和可行解精英集，若对于粒子i在第k代时的两个目标函数j1(xi(k))和j2(xi(k))，如果满足j1(xi(k))＜j1(xj(k))且j2(xi(k))＜j2(xj(k))，那么称粒子i的解xi(k)为非支配解；所有非支配解构成的集合称为一般精英集，也称为pareto最优解集；称pareto最优解集中满足稳定性约束条件的解为可行pareto最优解，反之称为不可行pareto最优解；由于可行pareto最优解所代表的频率反馈系数矩阵k所对应的含有调速器的电网节点动力学模型式(2)是稳定的，所以由式(3)所得的目标函数j1较小；而不可行pareto最优解所代表的频率反馈系数矩阵k所对应的含有调速器的电网节点动力学模型式(2)是不稳定的，所以由式(3)所得的目标函数j1较大；为此，本发明步骤定义一个阈值jthreshold，当j1≤jthreshold时，j1所对应的pareto最优解是可行pareto最优解；反之，当j1＞jthreshold时，j1所对应的pareto最优解是不可行pareto最优解；设第k次迭代产生了nfea(k)可行pareto最优解，如果此时可行解精英集中已经存储了n(k)个历史可行解，若nfea(k)+n(k)＜nfeamax，则将nfea(k)个可行pareto最优解对应的粒子的信息全部存储到可行解精英集中；否则，在nfea(k)+n(k)个可行pareto最优解采用欧式距离法选出距离全局向导粒子距离最小的前nfeamax个可行解存入可行解精英集；一般精英集不设置最大容量，存储每一代的所有pareto最优解；实施例中，由于j1(x1(1))＜j1(x2(1))但j2(x1(1))＞j2(x2(1))，所以x1(1)和x2(1)是非支配关系。同理，判断x1(1)、x2(1)和x3(1)两两之间的支配关系，可知x1(1)、x2(1)和x3(1)都是非支配解，则x1(1)、x2(1)和x3(1)构成pareto最优解集，即一般精英集。设置阈值jthreshold＝50，由于j1(x1(1))＝25.67＜jthreshold，j1(x2(1))＝29.33＜jthreshold以及j1(x3(1))＝26.51＜jthreshold，所以pareto最优解x1(1)、x2(1)和x3(1)均为可行pareto最优解，即nfea(1)＝3。可行解精英集中之前未存储任何粒子的信息，即n(1)＝0，即可得到nfea(1)+n(1)＝3＜nfeamax＝50，则将3个可行pareto最优解x1(1)、x2(1)和x3(1)存储到可行解精英集中。步骤5，更新个体向导粒子的位置向量和全局向导粒子的位置向量xgbest，若k＝1，则每个初始粒子的位置向量设置为该粒子的个体向导粒子的位置向量m＝1...n；否则，根据粒子m第k代的两个目标函数更新粒子m的个体向导粒子的目标函数为jm1min和jm2min，具体过程是：若j1(xm(k))＜jm1min且j2(xm(k))＜jm2min，则令jm1min＝j1(xm(k))，jm2min＝j2(xm(k))且否则，jm1min、jm2min和保持不变；采用欧式距离法在一般精英集中确定全局向导粒子xgbest，具体过程是：首先，将一般精英集中所有粒子以目标函数j1为横坐标，目标函数j2为纵坐标排列在二维笛卡尔坐标系中，找出一般精英集中粒子的目标函数j1的最小值j1min和目标函数j2的最小值j2min；然后，过j1min沿垂直横轴方向作延长线，过j2min沿垂直纵轴方向作延长线，定义两延长线交点为期望全局向导粒子的性能指标所在位置(j1min,j2min)，如图1所示；最后，比较笛卡尔坐标系中每个粒子的性能指标到(j1min,j2min)的距离，性能指标距(j1min,j2min)的欧式距离最小的粒子为全局向导粒子xgbest；对于本实施例，在k＝1时，则第1个粒子的个体向导粒子的位置向量第2个粒子的个体向导粒子的位置向量第2个粒子的个体向导粒子的位置向量参照图1，采用欧式距离法在一般精英集中确定全局向导粒子。以性能目标函数j1为横轴，以稀疏性目标函数j2为纵轴，构建二维直角坐标系，则第1个粒子的坐标为(25.67,100)，第2个粒子的坐标为(29.33,98)，第3个粒子的坐标为(26.51,98)，那么期望全局向导粒子的性能指标所在位置为(25.67,98)。3个粒子的坐标与(25.67,98)的欧式距离分别为：因此，d3＜d1＜d2，所以全局向导粒子为第3个粒子，对应的全局向导粒子的位置向量xgbest＝x3(1)；若k＝1，也可进行类似操作。步骤6，更新每个粒子的位置和速度，速度向量和位置向量的更新表达式为式(5)和式(6)：其中，xm(k)为第k代时粒子m的位置向量，vm(k)为第k代时粒子m的速度向量，m＝1,2,…,m，m是粒子总数，粒子m的位置向量的第p个元素，是粒子m的速度向量的第p个元素，个体向导粒子的位置向量和全局向导粒子的位置向量xgbest由步骤5产生，是个体向导粒子的位置向量的第p个元素，xp,gbest是全局向导粒子的位置向量xgbest的第p个元素；c1,c2是加速常数，rm,1,rm,2是对应粒子m的两个随机数向量，是向量rm,1第p个元素，为0-1均匀分布的随机数，是向量rm,2第p个元素，也为0-1均匀分布的随机数，p＝1,2,…,n2，n是频率反馈系数矩阵k的维数；ω(k)为第k代的权重系数，按照式(7)所示的线性递减权值确定ω(k)的大小：ω(k)＝(ω0-ωend)(kmax-k)/kmax+ωend(7)对于本实施例，由表2可知，初始权重ω0＝0.9，最终惯性权重ωend＝0.4，最大迭代次数kmax＝200，则根据线性递减权重的式(7)，第1代的权重ω(1)为：ω(1)＝(0.9-0.4)(200-1)/200+0.4＝0.8975。将步骤3产生的个体向导粒子及和全局向导粒子xgbest＝x3(1)的对应元素以及第1代权重系数ω(1)＝0.8975，加速常数c1＝2.8，c2＝1.3，以及随机数向量r1,1，r2,1，r2,1，r2,2，r3,1，r3,2的对应元素代入式(5)和式(6)进行粒子群算法的迭代，3个粒子第2代的位置向量x1(2)、x2(2)及x3(2)和速度向量v1(2)、v2(2)及v3(2)的第p个元素的具体迭代表达式为：由上面3组方程即可得到3个粒子第2代的位置向量x1(2)、x2(2)及x3(2)和速度向量v1(2)、v2(2)及v3(2)。步骤7，判断迭代次数是否到达初始化所设置的最大迭代次数kmax，若没有达到最大迭代次数kmax，则返回步骤3开始重新迭代；否则，结束循环，以可行解精英集为最终的优化解集，在可行解精英集中寻找稀疏性目标函数j2(xm(kmax))最小对应的粒子，设该粒子的位置向量为xm(kmax)，则位置向量xm(kmax)的第1到第n个元素对应反馈系数矩阵km的第1行，第n+1到第2n个元素对应反馈系数矩阵km的第2行，依次类推，得到反馈系数矩阵km，即为最终优化结果。对于上面实施例，因为此时k＝2≤kmax＝200，所以返回步骤3。当运行到第kmax＝200代时，可行解精英集中共有50个可行pareto最优解，其中稀疏性目标函数最小对应的粒子的位置向量为xm(kmax)，该粒子对应的性能目标函数为j1(xm(kmax))＝15.966，对应的稀疏性目标函数为j2(xm(kmax))＝2。将该粒子位置向量为xm(kmax)变换为频率系数矩阵km，具体方式为：将该粒子的位置向量xm(kmax)的第1个元素到第10个元素变为矩阵km的第1行，第11个元素到第20个元素为矩阵k的第2行，以此类推，频率系数矩阵km即为最终优化结果。对比仿真验证对比方法为文献[ful,fardadm,mr.designofoptimalsparsefeedbackgainsviathealternatingdirectionmethodofmultipliers[j].ieeetransactionsonautomaticcontrol,2013,58(9):2426-2431]中提出的基于交叉乘子法的稀疏提升算法，运用加权系数的遍历和梯度信息，加权遍历使算法整体效率较低，且容易陷入局部极值。参照图3，是不同方法频率反馈系数矩阵k非零元素的示意图，其中对角线的实心圆点代表发电节点自身的频率反馈，无需通信线路；非对角线的空心圆点代表发电节点之间的频率反馈，需要通信线路。图3是对比方法的结果，可见发电节点1与发电节点3之间，发电节点1与发电节点6之间以及发电节点1与发电节点9之间存在频率反馈，因此需要3条通信线路。图4是本发明方法的结果，可见发电节点5与发电节点6之间以及发电节点7与发电节点8之间存在频率反馈，因此需要2条通信线路。可见本发明的方法比对比方法减少了1条通信线路。如表3，是ieee-39标准测试系统频率反馈系数矩阵采用不同方法时的性能指标和对应的通信线路数。可见本发明方法在保证频率响应性能的同时，明显降低了成本代价。表3、ieee-39标准测试系统频率反馈系数矩阵主要优化结果算法性能指标通信线路条数集中式控制15.69645分散式控制18.8650区域式控制(对比方法)16.3213区域式控制(本发明方法)15.9662接下来，通过算法运行时间衡量本发明方法和对比方法的效率。由于粒子群算法的运行时间具有一定随机性，取30次实验的平均运行时间作为参考，每次实验均基于相同的计算机环境(intelcorei5-3470,3.20ghzcpu,4gbram)。具体运行时间如表4所示。由表4可以看出本发明方法求解频率反馈系数矩阵的运行时间比对比方法少，因此本发明方法具有更高的效率。表4、求解ieee-39标准测试系统频率反馈系数矩阵的运行时间算法运行时间(单位：秒/s)对比方法22.35本发明方法19.21综上所述，从对比仿真结果可以看出，本发明方法所得的优化结果不论是在频率响应性能还是成本代价方面都具有优越性，同时具有更高的运算效率。当前第1页12