CN108664749B - 一种发电机调速器频率反馈系数的优化方法 - Google Patents

Abstract

本发明公开了一种发电机调速器频率反馈系数的优化方法，步骤包括：1)建立含有调速器的电网节点动力学模型；2)初始化多目标粒子群算法的参数；3)计算粒子群中每个粒子的性能目标函数和稀疏性目标函数；4)更新一般精英集和可行解精英集；5)更新个体向导粒子的位置向量

和全局向导粒子的位置向量X^gbest；6)更新每个粒子的位置和速度；7)在可行解精英集中寻找稀疏性目标函数J₂(X_m(k_max))最小对应的粒子，设该粒子的位置向量为X_m(k_max)，得到反馈系数矩阵K_m，即为最终优化结果。本发明的方法，在降低成本代价的同时能够提升电网频率响应性能。

Description

一种发电机调速器频率反馈系数的优化方法

技术领域

本发明属于电网频率控制技术领域，涉及一种发电机调速器频率反馈系数的优化方法。

背景技术

调速器是电网中发电节点一次调频的主要控制器，是根据发电节点的频率对输入机械功率进行反馈调节的手段，对电网频率稳定性具有重要作用，频率反馈控制系数对调速器性能和成本都有重要影响。随着电力系统同步相角测量单元(Phasor MeasurementUnit，PMU)和广域监测系统(Wide Area Measurement System，WAMS)等技术的发展，使节点间(区域间)可以实时通信，以实现利用其它节点信息的区域协同控制。

根据通信结构，控制方式分为分散式控制、区域式控制和集中式控制共3种。传统的分散式控制是根据自身节点信息进行功率控制，无需通信线路，成本代价小，但频率调节性能较低；集中式控制利用全网所有发电节点信息进行反馈控制，需要电网中的发电节点两两之间均有通信线路，频率调节性能高，但成本代价较大；区域式控制利用部分发电节点信息进行反馈控制，只需要电网中部分发电节点之间有通信线路，兼顾了分散式控制和集中式控制的优点。然而，确定利用哪个节点信息进行反馈控制和反馈系数大小以平衡性能和代价这对矛盾仍是需要解决的问题。

发明内容

本发明提供了一种发电机调速器频率反馈系数的优化方法，解决现有技术难以确定调速器的频率反馈系数矩阵，难以平衡控制性能和通信代价的问题。

本发明采用的技术方案是，一种发电机调速器频率反馈系数的优化方法，按照以下步骤实施：

步骤1，建立含有调速器的电网节点动力学模型，

含有调速器的电网节点动力学模型利用传统同步发电机转子的摇摆方程来描述，表达式为下式(1)：

其中，下标i,j＝1,2,…,n为节点编号,n是电网中发电节点个数；J_i,D_i分别是节点i的惯性常量和阻尼系数，θ_i,ω_i分别是在同步旋转参考坐标系下节点i的相位和角频率，

表示节点i相位的一阶导数，即节点i的角频率ω_i；θ_j是节点j在同步旋转参考坐标系下的相位；当电力系统额定频率为f^*＝50Hz，则额定角频率为ω^*＝2πf^*；

为标称机械输入功率，

表示取导纳虚部，即不考虑引起线路损耗的电阻部分，L_ij是电网的拉普拉斯矩阵L第i行第j列的元素；E_i,E_j是节点i,j的端电压的幅值，Y_ij∈Y，矩阵Y是电网经过Kron化简后仅保留发电节点的导纳矩阵；u_i是节点i的控制量，β_ij为节点i与节点j之间的频率反馈系数；

将式(1)改写为矩阵结构的表达式如下式(2)：

其中，n×1维的列向量θ＝[θ₁,θ₂,…,θ_n]^T为相位向量，n×1维的列向量ω＝[ω₁,ω₂,…,ω_n]^T为角频率向量；n×1维的列向量u＝[u₁,u₂,…,u_n]^T为控制向量；n×n维矩阵J为对角矩阵，第i行对角线元素为节点i的惯性常量J_i；n×n维矩阵D为对角矩阵，第i行对角线元素为节点i的阻尼系数D_i；n×n维矩阵K为频率反馈系数矩阵，第i行第j列元素为β_ij；

将实际电网的参数代入式(2)，得到含有调速器的电网节点动力学模型；

步骤2，初始化多目标粒子群算法的参数，

设置惯性权重ω、加速常数c₁,c₂、最大迭代次数k_max，k为粒子群算法的迭代次数，令k＝1，粒子的总数为M个，每个粒子的位置和速度均是长度为n²的行向量，n为电网中发电节点个数，将第m个粒子在第k次迭代时的位置向量和速度向量分别记为X_m(k)和V_m(k)，初始化第1代时每个粒子的位置向量和速度向量中的每个元素为0-1均匀分布的随机数，设置可行解精英集最大容量N_{fea max}，即最多存储N_{fea max}个粒子的信息，设置初始惯性权重ω₀，以及最后一次迭代时的惯性权重ω_end；

初始化时，即迭代次数k＝1时，第m个粒子的位置向量和速度向量分别为X_m(1)和V_m(1)；X_m(1),V_m(1)为0-1均匀分布的随机向量；每个体向导粒子的位置向量

为该粒子的初始粒子的位置向量；

步骤3，计算粒子群中每个粒子的性能目标函数和稀疏性目标函数，

每个粒子的位置向量X_m(k)都代表一个解，即一个频率反馈系数矩阵K_m，位置向量的第1到第n个元素对应反馈系数矩阵K_m的第1行，第n+1到第2n个元素对应系数矩阵K_m的第2行，依次类推；

将每个粒子的位置向量对应的频率反馈系数矩阵K_m代入式(2)中进行时域仿真，采用式(3)计算第k代每个粒子的性能目标函数J₁(x)：

式(3)中，ω_m(t)代表采用反馈系数矩阵K_m时第t秒的角频率向量；u_m(t)代表采用反馈系数矩阵K_m时第t秒的控制向量；矩阵Q和矩阵R为单位对角矩阵；

采用式(4)计算第k代每个粒子的稀疏性目标函数J₂(x)：

J₂(X_m(k))＝||K_m||₀ (4)

步骤4，更新一般精英集和可行解精英集，

若对于粒子i在第k代时的两个目标函数J₁(X_i(k))和J₂(X_i(k))，如果满足J₁(X_i(k))＜J₁(X_j(k))且J₂(X_i(k))＜J₂(X_j(k))，那么称粒子i的解X_i(k)为非支配解；所有非支配解构成的集合称为一般精英集，也称为Pareto最优解集；

称Pareto最优解集中满足稳定性约束条件的解为可行Pareto最优解，反之称为不可行Pareto最优解；定义一个阈值J_threshold，当J₁≤J_threshold时，J₁所对应的Pareto最优解是可行Pareto最优解；反之，当J₁＞J_threshold时，J₁所对应的Pareto最优解是不可行Pareto最优解；

设第k次迭代产生了N_fea(k)可行Pareto最优解，如果此时可行解精英集中已经存储了N(k)个历史可行解，若N_fea(k)+N(k)＜N_{fea max}，则将N_fea(k)个可行Pareto最优解对应的粒子的信息全部存储到可行解精英集中；否则，在N_fea(k)+N(k)个可行Pareto最优解采用欧式距离法选出距离全局向导粒子距离最小的前N_feamax个可行解存入可行解精英集；

步骤5，更新个体向导粒子的位置向量

和全局向导粒子的位置向量X^gbest，

根据粒子m第k代的两个目标函数更新粒子m的个体向导粒子的目标函数为J_m1min和J_m2min，

具体过程是：

若J₁(X_m(k))＜J_m1min且J₂(X_m(k))＜J_m2min，则令J_m1min＝J₁(X_m(k))，J_m2min＝J₂(X_m(k))且

否则，J_m1min、J_m2min和

保持不变；

采用欧式距离法在一般精英集中确定全局向导粒子X^gbest；

步骤6，更新每个粒子的位置和速度，

速度向量和位置向量的更新表达式为式(5)和式(6)：

其中，X_m(k)为第k代时粒子m的位置向量，V_m(k)为第k代时粒子m的速度向量，m＝1,2,…,M，M是粒子总数，

粒子m的位置向量的第p个元素，

是粒子m的速度向量的第p个元素，个体向导粒子的位置向量

和全局向导粒子的位置向量X^gbest由步骤5产生，

是个体向导粒子的位置向量

的第p个元素，X^p,gbest是全局向导粒子的位置向量X^gbest的第p个元素；c₁,c₂是加速常数，r_m,1,r_m,2是对应粒子m的两个随机数向量，

是向量r_m,1第p个元素，为0-1均匀分布的随机数，

是向量r_m,2第p个元素，也为0-1均匀分布的随机数，p＝1,2,…,n²，n是频率反馈系数矩阵K的维数；ω(k)为第k代的权重系数，按照式(7)所示的线性递减权值确定ω(k)的大小：

ω(k)＝(ω₀-ω_end)(k_max-k)/k_max+ω_end (7)

步骤7，判断迭代次数是否到达初始化所设置的最大迭代次数k_max，若没有达到最大迭代次数k_max，则返回步骤3；否则，结束循环，以可行解精英集为最终的优化解集，在可行解精英集中寻找稀疏性目标函数J₂(X_m(k_max))最小对应的粒子，设该粒子的位置向量为X_m(k_max)，则位置向量X_m(k_max)的第1到第n个元素对应反馈系数矩阵K_m的第1行，第n+1到第2n个元素对应反馈系数矩阵K_m的第2行，依次类推，得到反馈系数矩阵K_m，即为最终优化结果。

本发明的有益效果是，该方法基于多目标粒子群的稀疏提升算法，在降低成本代价的同时能够提升电网频率响应性能。

附图说明

图1是本发明方法采用欧式距离法确定全局向导粒子的示意图；

图2是本发明方法实施例采用的IEEE-39标准测试系统示意图；

图3是对比方法针对IEEE-39标准测试系统得到频率反馈系数矩阵的非零元素示意图；

图4是本发明方法针对的IEEE-39标准测试系统得到频率反馈系数矩阵的非零元素示意图。

具体实施方式

下面基于IEEE-39标准测试系统，对本发明调速器频率反馈系数矩阵的优化多对目标粒子群的稀疏提升优化过程进行具体说明。

本发明方法，按照以下步骤实施：

步骤1，建立含有调速器的电网节点动力学模型，

含有调速器的电网节点动力学模型利用传统同步发电机转子的摇摆方程来描述，表达式为下式(1)：

其中，下标i,j＝1,2,…,n为节点编号，n是电网中发电节点个数；J_i,D_i分别是节点i的惯性常量和阻尼系数，θ_i,ω_i分别是在同步旋转参考坐标系下节点i的相位和角频率，

表示节点i相位的一阶导数，即节点i的角频率ω_i；θ_j是节点j在同步旋转参考坐标系下的相位；一般情况下，电力系统额定频率为f^*＝50Hz，则额定角频率为ω^*＝2πf^*；

为标称机械输入功率，

表示取导纳虚部，即不考虑引起线路损耗的电阻部分，L_ij是电网的拉普拉斯矩阵L第i行第j列的元素；E_i,E_j是节点i,j的端电压的幅值，Y_ij∈Y，矩阵Y是电网经过Kron化简后仅保留发电节点的导纳矩阵；u_i是节点i的控制量，β_ij为节点i与节点j之间的频率反馈系数；

此处假设：1)传输线路损耗忽略不计；2)节点的端电压E_i,E_j工作在其额定工作点，为常数；

将式(1)改写为矩阵结构的表达式如下式(2)：

其中，n×1维的列向量θ＝[θ₁,θ₂,…,θ_n]^T为相位向量，n×1维的列向量ω＝[ω₁,ω₂,…,ω_n]^T为角频率向量；n×1维的列向量u＝[u₁,u₂,…,u_n]^T为控制向量；n×n维矩阵J为对角矩阵，第i行对角线元素为节点i的惯性常量J_i；n×n维矩阵D为对角矩阵，第i行对角线元素为节点i的阻尼系数D_i；n×n维矩阵K为频率反馈系数矩阵，第i行第j列元素为β_ij；

参照图2，实施例采用IEEE-39标准测试系统，其中共包含10个发电机，分别位于母线30至母线39上，令这10个发电机所在的母线依次为发电节点1至发电节点10，则节点数为n＝10；在电网节点动力学模型中，节点i的惯性常量J_i、阻尼系数D_i、标称机械输入功率

初始相位θ_i和初始角频率ω_i具体数值，如表1所示。

表1、IEEE-39标准测试系统含调速器的电网节点动力学模型参数

将实际电网的参数代入式(2)，得到含有调速器的电网节点动力学模型。

步骤2，初始化多目标粒子群算法的参数，

设置惯性权重ω、加速常数c₁,c₂、最大迭代次数k_max，k为粒子群算法的迭代次数，令k＝1，粒子的总数为M个，每个粒子的位置和速度均是长度为n²的行向量，n为电网中发电节点个数，将第m个粒子在第k次迭代时的位置向量和速度向量分别记为X_m(k)和V_m(k)，初始化第1代时每个粒子的位置向量和速度向量中的每个元素为0-1均匀分布的随机数，设置可行解精英集最大容量N_{fea max}，即最多存储N_{fea max}个粒子的信息，设置初始惯性权重ω₀，以及最后一次迭代时的惯性权重ω_end；

本实施例多目标粒子群算法的具体参数值见表2，以下按照3个粒子M＝3为例对多目标粒子群算法进行描述。

表2 多目标粒子群算法的具体参数

初始化时，即迭代次数k＝1时，第1个粒子的位置向量和速度向量分别为X₁(1)和V₁(1)；第2个粒子的位置向量和速度向量分别为X₂(1)和V₂(1)；第3个粒子的位置向量和速度向量分别为X₃(1)和V₃(1)；将3个粒子的位置向量X₁(1),X₂(1),X₃(1)设为0-1均匀分布的随机向量，将3个粒子的速度向量V₁(1),V₂(1),V₃(1)设为0-1均匀分布的随机向量；因为节点数n＝10，所以3个粒子的位置向量和速度向量均为n²×1维，即100×1维。

步骤3，计算粒子群中每个粒子的性能目标函数和稀疏性目标函数，

每个粒子的位置向量X_m(k)都代表一个解，即一个频率反馈系数矩阵K_m，位置向量的第1到第n个元素对应反馈系数矩阵K_m的第1行，第n+1到第2n个元素对应系数矩阵K_m的第2行，依次类推；

将每个粒子的位置向量对应的频率反馈系数矩阵K_m代入式(2)中进行时域仿真，采用式(3)计算第k代(本代)每个粒子的性能目标函数J₁(x)：

式(3)中，ω_m(t)代表采用反馈系数矩阵K_m时第t秒的角频率向量；u_m(t)代表采用反馈系数矩阵K_m时第t秒的控制向量；矩阵Q和矩阵R为单位对角矩阵；

采用式(4)计算第k代(本代)每个粒子对应的稀疏性目标函数J₂(x)：

J₂(X_m(k))＝||K_m||₀ (4)

本实施例中将每个粒子的位置向量转换为频率反馈系数矩阵K的形式，即第1个粒子的位置向量的第1个元素到第10个元素为矩阵K的第1行，第11个元素到第20个元素为矩阵K的第2行，以此类推；第2个粒子和第3个粒子采用同样的转换方式得到矩阵K；

将3个粒子所得的3个频率反馈系数矩阵K代入(2)进行仿真的结果分别代入式(4)，按照表1设置的参数值进行30s的时域仿真，用式(3)计算3个粒子的性能目标函数J₁分别为J₁(X₁(1))＝25.67，J₁(X₂(1))＝29.33，J₁(X₃(1))＝26.51。

将3个粒子所得的3个频率反馈系数矩阵K分别代入式(4)，计算3个粒子的稀疏性目标函数J₂分别为J₂(X₁(1))＝100，J₂(X₂(1))＝98，J₂(X₃(1))＝98。

步骤4，更新一般精英集和可行解精英集，

若对于粒子i在第k代时的两个目标函数J₁(X_i(k))和J₂(X_i(k))，如果满足J₁(X_i(k))＜J₁(X_j(k))且J₂(X_i(k))＜J₂(X_j(k))，那么称粒子i的解X_i(k)为非支配解；所有非支配解构成的集合称为一般精英集，也称为Pareto最优解集；

称Pareto最优解集中满足稳定性约束条件的解为可行Pareto最优解，反之称为不可行Pareto最优解；由于可行Pareto最优解所代表的频率反馈系数矩阵K所对应的含有调速器的电网节点动力学模型式(2)是稳定的，所以由式(3)所得的目标函数J₁较小；而不可行Pareto最优解所代表的频率反馈系数矩阵K所对应的含有调速器的电网节点动力学模型式(2)是不稳定的，所以由式(3)所得的目标函数J₁较大；为此，本发明步骤定义一个阈值J_threshold，当J₁≤J_threshold时，J₁所对应的Pareto最优解是可行Pareto最优解；反之，当J₁＞J_threshold时，J₁所对应的Pareto最优解是不可行Pareto最优解；

设第k次迭代产生了N_fea(k)可行Pareto最优解，如果此时可行解精英集中已经存储了N(k)个历史可行解，若N_fea(k)+N(k)＜N_{fea max}，则将N_fea(k)个可行Pareto最优解对应的粒子的信息全部存储到可行解精英集中；否则，在N_fea(k)+N(k)个可行Pareto最优解采用欧式距离法选出距离全局向导粒子距离最小的前N_{fea max}个可行解存入可行解精英集；

一般精英集不设置最大容量，存储每一代的所有Pareto最优解；

实施例中，由于J₁(X₁(1))＜J₁(X₂(1))但J₂(X₁(1))＞J₂(X₂(1))，所以X₁(1)和X₂(1)是非支配关系。同理，判断X₁(1)、X₂(1)和X₃(1)两两之间的支配关系，可知X₁(1)、X₂(1)和X₃(1)都是非支配解，则X₁(1)、X₂(1)和X₃(1)构成Pareto最优解集，即一般精英集。

设置阈值J_threshold＝50，由于J₁(X₁(1))＝25.67＜J_threshold，

J₁(X₂(1))＝29.33＜J_threshold以及J₁(X₃(1))＝26.51＜J_threshold，所以Pareto最优解X₁(1)、X₂(1)和X₃(1)均为可行Pareto最优解，即N_fea(1)＝3。可行解精英集中之前未存储任何粒子的信息，即N(1)＝0，即可得到N_fea(1)+N(1)＝3＜N_{fea max}＝50，则将3个可行Pareto最优解X₁(1)、X₂(1)和X₃(1)存储到可行解精英集中。

步骤5，更新个体向导粒子的位置向量

和全局向导粒子的位置向量X^gbest，

若k＝1，则每个初始粒子的位置向量设置为该粒子的个体向导粒子的位置向量

m＝1...n；否则，根据粒子m第k代的两个目标函数更新粒子m的个体向导粒子的目标函数为J_m1min和J_m2min，

具体过程是：

若J₁(X_m(k))＜J_m1min且J₂(X_m(k))＜J_m2min，则令J_m1min＝J₁(X_m(k))，J_m2min＝J₂(X_m(k))且

否则，J_m1min、J_m2min和

保持不变；

采用欧式距离法在一般精英集中确定全局向导粒子X^gbest，具体过程是：

首先，将一般精英集中所有粒子以目标函数J₁为横坐标，目标函数J₂为纵坐标排列在二维笛卡尔坐标系中，找出一般精英集中粒子的目标函数J₁的最小值J_1min和目标函数J₂的最小值J_2min；

然后，过J_1min沿垂直横轴方向作延长线，过J_2min沿垂直纵轴方向作延长线，定义两延长线交点为期望全局向导粒子的性能指标所在位置(J_1min,J_2min)，如图1所示；

最后，比较笛卡尔坐标系中每个粒子的性能指标到(J_1min,J_2min)的距离，性能指标距(J_1min,J_2min)的欧式距离最小的粒子为全局向导粒子X^gbest；

对于本实施例，在k＝1时，则第1个粒子的个体向导粒子的位置向量

第2个粒子的个体向导粒子的位置向量

第2个粒子的个体向导粒子的位置向量

参照图1，采用欧式距离法在一般精英集中确定全局向导粒子。以性能目标函数J₁为横轴，以稀疏性目标函数J₂为纵轴，构建二维直角坐标系，则第1个粒子的坐标为(25.67,100)，第2个粒子的坐标为(29.33,98)，第3个粒子的坐标为(26.51,98)，那么期望全局向导粒子的性能指标所在位置为(25.67,98)。3个粒子的坐标与(25.67,98)的欧式距离分别为：

因此，D₃＜D₁＜D₂，所以全局向导粒子为第3个粒子，对应的全局向导粒子的位置向量X^gbest＝X₃(1)；若k＝1，也可进行类似操作。

步骤6，更新每个粒子的位置和速度，

速度向量和位置向量的更新表达式为式(5)和式(6)：

其中，X_m(k)为第k代时粒子m的位置向量，V_m(k)为第k代时粒子m的速度向量，m＝1,2,…,M，M是粒子总数，

粒子m的位置向量的第p个元素，

是粒子m的速度向量的第p个元素，个体向导粒子的位置向量

和全局向导粒子的位置向量X^gbest由步骤5产生，

是个体向导粒子的位置向量

的第p个元素，X^p,gbest是全局向导粒子的位置向量X^gbest的第p个元素；c₁,c₂是加速常数，r_m,1,r_m,2是对应粒子m的两个随机数向量，

是向量r_m,1第p个元素，为0-1均匀分布的随机数，

是向量r_m,2第p个元素，也为0-1均匀分布的随机数，p＝1,2,…,n²，n是频率反馈系数矩阵K的维数；ω(k)为第k代的权重系数，按照式(7)所示的线性递减权值确定ω(k)的大小：

ω(k)＝(ω₀-ω_end)(k_max-k)/k_max+ω_end (7)

对于本实施例，由表2可知，初始权重ω₀＝0.9，最终惯性权重ω_end＝0.4，最大迭代次数k_max＝200，则根据线性递减权重的式(7)，第1代的权重ω(1)为：ω(1)＝(0.9-0.4)(200-1)/200+0.4＝0.8975。

将步骤3产生的个体向导粒子

及

和全局向导粒子X^gbest＝X₃(1)的对应元素以及第1代权重系数ω(1)＝0.8975，加速常数c₁＝2.8，c₂＝1.3，以及随机数向量r_1,1，r_2,1，r_2,1，r_2,2，r_3,1，r_3,2的对应元素代入式(5)和式(6)进行粒子群算法的迭代，3个粒子第2代的位置向量X₁(2)、X₂(2)及X₃(2)和速度向量V₁(2)、V₂(2)及V₃(2)的第p个元素的具体迭代表达式为：

由上面3组方程即可得到3个粒子第2代的位置向量X₁(2)、X₂(2)及X₃(2)和速度向量V₁(2)、V₂(2)及V₃(2)。

步骤7，判断迭代次数是否到达初始化所设置的最大迭代次数k_max，若没有达到最大迭代次数k_max，则返回步骤3开始重新迭代；否则，结束循环，以可行解精英集为最终的优化解集，在可行解精英集中寻找稀疏性目标函数J₂(X_m(k_max))最小对应的粒子，设该粒子的位置向量为X_m(k_max)，则位置向量X_m(k_max)的第1到第n个元素对应反馈系数矩阵K_m的第1行，第n+1到第2n个元素对应反馈系数矩阵K_m的第2行，依次类推，得到反馈系数矩阵K_m，即为最终优化结果。

对于上面实施例，因为此时k＝2≤k_max＝200，所以返回步骤3。

当运行到第k_max＝200代时，可行解精英集中共有50个可行Pareto最优解，其中稀疏性目标函数最小对应的粒子的位置向量为X_m(k_max)，该粒子对应的性能目标函数为J₁(X_m(k_max))＝15.966，对应的稀疏性目标函数为J₂(X_m(k_max))＝2。将该粒子位置向量为X_m(k_max)变换为频率系数矩阵K_m，具体方式为：将该粒子的位置向量X_m(k_max)的第1个元素到第10个元素变为矩阵K_m的第1行，第11个元素到第20个元素为矩阵K的第2行，以此类推，频率系数矩阵K_m即为最终优化结果。

对比仿真验证

对比方法为文献[Fu L,Fardad M,

M R.Design of Optimal SparseFeedback Gains via the Alternating Direction Method ofMultipliers[J].IEEETransactions on Automatic Control,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.696	45
分散式控制	18.865	0
			区域式控制(对比方法)	16.321	3
区域式控制(本发明方法)	15.966	2

接下来，通过算法运行时间衡量本发明方法和对比方法的效率。由于粒子群算法的运行时间具有一定随机性，取30次实验的平均运行时间作为参考，每次实验均基于相同的计算机环境(Intel Corei5-3470,3.20GHz CPU,4GB RAM)。具体运行时间如表4所示。由表4可以看出本发明方法求解频率反馈系数矩阵的运行时间比对比方法少，因此本发明方法具有更高的效率。

表4、求解IEEE-39标准测试系统频率反馈系数矩阵的运行时间

算法	运行时间(单位：秒/s)
		对比方法	22.35
本发明方法	19.21

综上所述，从对比仿真结果可以看出，本发明方法所得的优化结果不论是在频率响应性能还是成本代价方面都具有优越性，同时具有更高的运算效率。

Claims (3)

1.一种发电机调速器频率反馈系数的优化方法，其特征在于，按照以下步骤实施：

步骤1，建立含有调速器的电网节点动力学模型，

含有调速器的电网节点动力学模型利用传统同步发电机转子的摇摆方程来描述，表达式为下式(1)：

其中，下标i,j＝1,2,…,n为节点编号,n是电网中发电节点个数；J_i,D_i分别是节点i的惯性常量和阻尼系数，θ_i,ω_i分别是在同步旋转参考坐标系下节点i的相位和角频率，

表示节点i相位的一阶导数，即节点i的角频率ω_i；θ_j是节点j在同步旋转参考坐标系下的相位；当电力系统额定频率为f^*＝50Hz，则额定角频率为ω^*＝2πf^*；

为标称机械输入功率，

表示取导纳虚部，即不考虑引起线路损耗的电阻部分，L_ij是电网的拉普拉斯矩阵L第i行第j列的元素；E_i,E_j是节点i,j的端电压的幅值，Y_ij∈Y，矩阵Y是电网经过Kron化简后仅保留发电节点的导纳矩阵；u_i是节点i的控制量，β_ij为节点i与节点j之间的频率反馈系数；

将式(1)改写为矩阵结构的表达式如下式(2)：

其中，n×1维的列向量θ＝[θ₁,θ₂,…,θ_n]^T为相位向量，n×1维的列向量ω＝[ω₁,ω₂,…,ω_n]^T为角频率向量；n×1维的列向量u＝[u₁,u₂,…,u_n]^T为控制向量；n×n维矩阵J为对角矩阵，第i行对角线元素为节点i的惯性常量J_i；n×n维矩阵D为对角矩阵，第i行对角线元素为节点i的阻尼系数D_i；n×n维矩阵K为频率反馈系数矩阵，第i行第j列元素为β_ij；

将实际电网的参数代入式(2)，得到含有调速器的电网节点动力学模型；

步骤2，初始化多目标粒子群算法的参数，

设置惯性权重ω、加速常数c₁,c₂、最大迭代次数k_max，k为粒子群算法的迭代次数，令k＝1，粒子的总数为M个，每个粒子的位置和速度均是长度为n²的行向量，n为电网中发电节点个数，将第m个粒子在第k次迭代时的位置向量和速度向量分别记为X_m(k)和V_m(k)，初始化第1代时每个粒子的位置向量和速度向量中的每个元素为0-1均匀分布的随机数，设置可行解精英集最大容量N_feamax，即最多存储N_feamax个粒子的信息，设置初始惯性权重ω₀，以及最后一次迭代时的惯性权重ω_end；

初始化时，即迭代次数k＝1时，第m个粒子的位置向量和速度向量分别为X_m(1)和V_m(1)；X_m(1),V_m(1)为0-1均匀分布的随机向量；每个体向导粒子的位置向量

为该粒子的初始粒子的位置向量；

步骤3，计算粒子群中每个粒子的性能目标函数和稀疏性目标函数，

每个粒子的位置向量X_m(k)都代表一个解，即一个频率反馈系数矩阵K_m，位置向量的第1到第n个元素对应反馈系数矩阵K_m的第1行，第n+1到第2n个元素对应系数矩阵K_m的第2行，依次类推；

将每个粒子的位置向量对应的频率反馈系数矩阵K_m代入式(2)中进行时域仿真，采用式(3)计算第k代每个粒子的性能目标函数J₁(x)：

式(3)中，ω_m(t)代表采用反馈系数矩阵K_m时第t秒的角频率向量；u_m(t)代表采用反馈系数矩阵K_m时第t秒的控制向量；矩阵Q和矩阵R为单位对角矩阵；

采用式(4)计算第k代每个粒子的稀疏性目标函数J₂(x)：

J₂(X_m(k))＝||K_m||₀ (4)

步骤4，更新一般精英集和可行解精英集，

若对于粒子i在第k代时的两个目标函数J₁(X_i(k))和J₂(X_i(k))，如果满足J₁(X_i(k))＜J₁(X_j(k))且J₂(X_i(k))＜J₂(X_j(k))，那么称粒子i的解X_i(k)为非支配解；所有非支配解构成的集合称为一般精英集，也称为Pareto最优解集；

称Pareto最优解集中满足稳定性约束条件的解为可行Pareto最优解，反之称为不可行Pareto最优解；定义一个阈值J_threshold，当J₁≤J_threshold时，J₁所对应的Pareto最优解是可行Pareto最优解；反之，当J₁＞J_threshold时，J₁所对应的Pareto最优解是不可行Pareto最优解；

设第k次迭代产生了N_fea(k)可行Pareto最优解，如果此时可行解精英集中已经存储了N(k)个历史可行解，若N_fea(k)+N(k)＜N_feamax，则将N_fea(k)个可行Pareto最优解对应的粒子的信息全部存储到可行解精英集中；否则，在N_fea(k)+N(k)个可行Pareto最优解采用欧式距离法选出距离全局向导粒子距离最小的前N_feamax个可行解存入可行解精英集；

步骤5，更新个体向导粒子的位置向量

和全局向导粒子的位置向量X^gbest，

根据粒子m第k代的两个目标函数更新粒子m的个体向导粒子的目标函数为J_m1min和J_m2min，

具体过程是：

若J₁(X_m(k))＜J_m1min且J₂(X_m(k))＜J_m2min，则令J_m1min＝J₁(X_m(k))，J_m2min＝J₂(X_m(k))且

否则，J_m1min、J_m2min和

保持不变；

采用欧式距离法在一般精英集中确定全局向导粒子X^gbest；

步骤6，更新每个粒子的位置和速度，

速度向量和位置向量的更新表达式为式(5)和式(6)：

其中，X_m(k)为第k代时粒子m的位置向量，V_m(k)为第k代时粒子m的速度向量，m＝1,2,…,M，M是粒子总数，

粒子m的位置向量的第p个元素，

是粒子m的速度向量的第p个元素，个体向导粒子的位置向量

和全局向导粒子的位置向量X^gbest由步骤5产生，

是个体向导粒子的位置向量

的第p个元素，X^p,gbest是全局向导粒子的位置向量X^gbest的第p个元素；c₁,c₂是加速常数，r_m,1,r_m,2是对应粒子m的两个随机数向量，

是向量r_m,1第p个元素，为0-1均匀分布的随机数，

是向量r_m,2第p个元素，也为0-1均匀分布的随机数，p＝1,2,…,n²，n是频率反馈系数矩阵K的维数；ω(k)为第k代的权重系数，按照式(7)所示的线性递减权值确定ω(k)的大小：

ω(k)＝(ω₀-ω_end)(k_max-k)/k_max+ω_end (7)

步骤7，判断迭代次数是否到达初始化所设置的最大迭代次数k_max，若没有达到最大迭代次数k_max，则返回步骤3；否则，结束循环，以可行解精英集为最终的优化解集，在可行解精英集中寻找稀疏性目标函数J₂(X_m(k_max))最小对应的粒子，设该粒子的位置向量为X_m(k_max)，则位置向量X_m(k_max)的第1到第n个元素对应反馈系数矩阵K_m的第1行，第n+1到第2n个元素对应反馈系数矩阵K_m的第2行，依次类推，得到反馈系数矩阵K_m，即为最终优化结果。

2.根据权利要求1所述的发电机调速器频率反馈系数的优化方法，其特征在于：所述的步骤4中，一般精英集不设置最大容量，存储每一代的所有Pareto最优解。

3.根据权利要求1所述的发电机调速器频率反馈系数的优化方法，其特征在于：所述的步骤5中，采用欧式距离法在一般精英集中确定全局向导粒子X^gbest，具体过程是：

首先，将一般精英集中所有粒子以目标函数J₁为横坐标，目标函数J₂为纵坐标排列在二维笛卡尔坐标系中，找出一般精英集中粒子的目标函数J₁的最小值J_1min和目标函数J₂的最小值J_2min；

然后，过J_1min沿垂直横轴方向作延长线，过J_2min沿垂直纵轴方向作延长线，定义两延长线交点为期望全局向导粒子的性能指标所在位置(J_1min,J_2min)；

最后，比较笛卡尔坐标系中每个粒子的性能指标到(J_1min,J_2min)的距离，性能指标距(J_1min,J_2min)的欧式距离最小的粒子为全局向导粒子X^gbest。

Images