[发明专利]一种分布式电源接入电网的潮流计算方法

摘要

一种分布式电源接入电网的潮流计算方法，包括以下步骤S1读取电力系统的初始数据；S2确定采样次数N和输入随机变量的维数s；S3生成s×N阶采样矩阵；S4将采样次数初始化令n＝1；S5判断n和采样次数N的大小，若n＞N，直接输出变量的概率统计结果；若n≤N，转S6；S6确定风电和光伏发电出力模型，确定负荷随机模型；S7确定潮流计算模型；S8确定最优经济模型；S9潮流计算S10记录第n组节点电压、支路功率及发电成本等数据；S11进行下一轮潮流计算，t＝t+1，转S5。本发明可以较好的估计输出随机变量的概率分布，能有效地处理电力市场中的不确定性问题，本发明节约了调试的人力和物力，降低了的生产成本。

主权项

一种分布式电源接入电网的潮流计算方法，其特征在于，包括以下步骤：S1：读取电力系统的初始数据；S2：确定采样次数N和输入随机变量的维数s；S3：按照以下3步，生成s×N阶采样矩阵，形成点列中第个点(j＝1,…,s；n＝1,…)的步骤如下：S3‑1：把第N‑1个整数用2进制数表示，即式(1)N‑1＝aR‑1aR‑2…a2a1 (1)其中an∈Zb，Zb＝{0,1,…,b‑1}，R为满足br≤N的r的最大值；S3‑2：对N‑1＝aR‑1aR‑2…a2a1进行排序，得到排序后的序列[d1d2…dn…dR]T为式(2)[d1d2...dn...dR]T=CN×Ni[a1a2...an...aR-1]T---(2)]]>其中，为生成矩阵，0≤dn≤b‑1；引入生成矩阵是为了重置a1a2…an…aR‑1中各个数字的位置；数字的位置经过重置后，每一维和其它维的数字大小相同，但排列顺序不同，从而保证了结果的均匀性；S3‑3：经过第S3‑2步的计算，可以表示为式(3)的2进制形式：xnj=0.d1d2...dR---(3)]]>最后，将2进制表示的根据式(2)转化为10进制数即可；S4：将采样次数初始化：令n＝1；S5：判断n和采样次数N的大小，若n＞N，直接输出变量的概率统计结果；若n≤N，转S6；S6：确定风电和光伏发电出力模型，确定负荷随机模型；S6‑1：风速服从韦布尔分布，风电场有功功率Pw的概率密度函数可表示为式(4)：f(Pw)=kk1c(Pw-k2k1c)exp[-(Pw-k2k1c)k]---(4)]]>式中：k,c分别为韦布尔分布的形状参数和尺度参数，k2＝‑k1vci，Pr为风机额定功率，vr,vci分别为额定风速和切入风速；风电处理为PQ节点，令潮流计算中风机功率因数恒定不变，则无功功率按下式(5)计算：式中：为功率因数角，对并网风机而言，位于第四象限，为负值；S6‑2：光伏出力随机模型一定时间段内，太阳光照强度可认为服从贝塔分布，则光伏电站输出功率Ppv的概率密度函数表示为式(6)：f(Ppv)=Γ(α+β)Γ(α)+Γ(β)(PpvRpv)α-1(1-PpvRpv)β-1---(6)]]>式中：Rpv＝Aηγmax为仿真最大输出功率，A为太阳能电池仿真总面积，η为仿真总的光电转换效率，γmax为一段时间内的最大光照强度，Γ为Gamma函数，α,β均为贝塔分布的形状参数；与风电相同，潮流计算中将光伏电站也作为PQ节点；S6‑3：负荷随机模型作为中长期的负荷预测结果，负荷的概率分布规律基本符合于正态分布；其均值和方差均可以由大量的历史统计数据得到；这样，负荷的有功和无功功率的概率密度函数分别为式(7)和(8)：f(P)=12πδPexp(-(P-μP)22δP2)---(7)]]>f(P)=12πδQexp(-(Q-μQ)22δQ2)---(8)]]>式中：μp为有功功率的均值，δp2为有功功率的方差，μQ为无功功率的均值，δQ2为无功功率的方差；S7：确定潮流计算模型S7‑1：目标函数构建的发电优化模型如下：minF=Σi∈Ng[CGpi(PGi(t))+CGqi(Qi(t))]+Σj∈NqCgqj·(Qgj(t))---(9)]]>式(9)中CGpi、CGqi为机组i的有功和无功发电费用函数，Cgqj为无功补偿装置j的无功发电费用函数，PGi(t)、QGi(t)为第i台发电机组在时段t的有功出力和无功出力，Qgj(t)为第j台无功补偿装置在时段t的无功出力；Ng、Nq为发电机节点个数与无功补偿设备个数；目标函数使各个时段系统的发电费用最小；S7‑2：等式约束等式约束为各个时段的节点潮流平衡约束：PGi-PDi-ViΣj=1NVj(Gijcosθij+Bijsinθij)=0---(10)]]>QGi+Qgi-QDi-ViΣj=1NVj(Gijsinθij+Bijcosθij)=0---(11)]]>式(10)、式(11)中：Vi、θi为节点电压与相角，θij＝θi‑θj；PDi、QDi为有功负荷与无功负荷；Gij、Bij为节点导纳矩阵的电导和电纳；S7‑3：不等式约束式(12)PGi‾≤PGi≤PGi‾QGi‾≤QGi≤QGi‾Qgi‾≤Qgi≤Qgi‾Vi‾≤Vi≤Vi‾Si≤Si‾,i∈Ngi∈Ngi∈Nqi∈Ni∈Nb---(12)]]>式中，PGi为发电机i有功出力上下限；QGi为发电机i无功出力上下限；Qgi为无功补偿设备i无功出力上下限；Vi为节点电压幅值上下限；为线路i持续输送容量极限(MVA)；N、Nb为节点集、支路集合；PGT,i(t+1)‑PGT,i≤Ri,up (13)PGT,i(t)‑PGT,i(t+1)≤Ri,down (14)式(13)中，Ri，up为第i台火电机组的向上爬坡速率；式(14)中，Ri，down为第i台火电机组的向下爬坡速率；S8：确定最优经济模型S8‑1：火电厂的发电费用火电燃煤机组的有功出力是以煤耗量为标准进行计费的，机组i有功出力费用函数CGpi以式(15)进行计算；式中ai、bi、ci为第i台火电机组的煤耗费用系数；CGpi=aiPGi2+biPGi+ci---(15)]]>无功机会成本是该发电机因输出无功功率而损失的有功功率发电容量所对应的利润；如果忽略原动机的出力极限，并假设该无功机会成本Cop(QGi)可表示如式(16)；CGqi(QGi)=Cop(QGi)=k[CGpi(SGi,max)-CGpi(SGi,max2-QGi2)]---(16)]]>将式(15)代入到式(16)中，并进行泰勒展开，保留到项，忽略高次项后并整理得到式(17)CGqi(QGi)=dQGi2+e---(17)]]>CGqi(QGi)为发电机组i的无功出力费用函数，SGi,max为发电机组i的额定视在功率，QGi为发电机组i的无功出力值，k为发电厂生产有功功率的利润率，一般为5％‑10％；S8‑2：水电厂的发电费用采用水电有功发电成本式(15)的形式进行计费，而其具体参数的取值近似与火电中ai，bi，ci相差m倍，m为火电运行成本与水电运行成本电价的比值，ai，bi，ci取值有微调变化，以区分同类电站的发电费用；水电厂的无功发电费用也采用类似火电厂无功出力费用，并按照式(16)的计费方式，其中的CGpi取相应水电厂的有功发电费用函数；S8‑3：光伏电站及风电场发电费用以最大限度优先调用新能源为准则，令补贴后光伏发电费用及风力发电费用低于火电及水电的上网发电价格，其有功费用函数的选取与水电的有功费用选取方式相同；S8‑4：无功补偿设备的发电费用以电容器、电抗器、同步调相机、SVC无功费用为固定成本表达式(18)：Cgqj(Qgj)=CfY×365×24×pQgj=fqQgj---(18)]]>其中Y为并联电容器的使用寿命，取15年；p为平均使用率，近似取为2/3，Cf为电容器单位容量的固定成本，平均可取为62500元/MVar，以此数据计算得出fq＝1.97；S9：潮流计算利用拉格朗日函数法来处理优化问题中的等式约束，从而将具有等式约束的优化问题转化为无约束的优化问题；利用对数障碍函数法的罚函数方法处理不等式约束，最后用牛顿法来求解无约束优化问题最优解；将非线性问题用以下数学公式表示：obj min.f(x)s.t. h(x)＝0 (19)g‾≤g(x)≤g‾]]>其中：min.f(x)为目标函数，是一个非线性函数；h(x)＝[h1(x),...,hm(x)]T为非线性等式约束条件，g(x)＝[g1(x),...,gr(x)]T为非线性不等式约束；假设在以上模型中共有k个变量，m个等式约束，r个不等式约束；用内点法求解问题(19)时，先将不等式约束转化为等式约束，同时构造障碍函数；为此先引入松弛变量l＞0，u＞0，l∈Rr，u∈Rr，将式(19)的不等式约束转化为等式约束，并把目标函数改造成障碍函数，可以得到以下优化问题A：objmin.f(x)-μ(Σj=1rlnlj+Σj=1rlnuj)]]>s.t. h(x)＝0 (20)g(x)+u-g‾=0]]>g(x)‑l‑g＝0其中扰动因子u＞0；当li或ui靠近边界时，以上函数趋于无穷大，因此满足以上障碍目标函数的极小解不可能在边界上找到，只能在满足l＞0，u>0时才可能得到最优解；这样，就通过目标函数的变换把含有不等式限制的优化问题变成了只含等式约束限制的优化问题A，因此可以直接用拉格朗日乘子法来求解；优化模型A的拉格朗日函数为：L=f(x)-yTh(x)-zT[g(x)-l-g‾]-wT[g(x)+u-g‾]-μ(Σj=1rlnlj+Σj=1rlnuj)---(21)]]>式中：y＝[y1,...,ym]T，z＝[z1,...,zr]T，w＝[w1,...,wr]T均为拉格朗日乘子；该问题极小值存在的必要条件是拉格朗日函数对所有变量及乘子的偏导数为0，从而将有约束优化转化为无约束优化，接下来可以使用现有技术中的牛顿法求解；S10：记录第n组节点电压、支路功率及发电成本等数据；S11：进行下一轮潮流计算，t＝t+1，转S5。