CN104578157B - 一种分布式电源接入电网的潮流计算方法 - Google Patents

Abstract

一种分布式电源接入电网的潮流计算方法，包括以下步骤：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。本发明可以较好的估计输出随机变量的概率分布，能有效地处理电力市场中的不确定性问题，本发明节约了调试的人力和物力，降低了的生产成本。

Description

一种分布式电源接入电网的潮流计算方法

技术领域

本发明涉及电力***配电网的应用领域，尤其涉及一种分布式电源接入电网的潮流计算方法。

背景技术

风能、太阳能是绿色清洁能源，大力发展风电、光伏有利于减少化石燃料消耗、降低碳排放水平。但因其具有间歇性和随机性的特点，对电力***运行控制提出更高要求。近年来，我国风电、光伏发电发展迅速，装机容量迅速增加，消纳出现困难，在电力***规划设计和运行控制中更全面地考虑风电场、光伏电站的特性，掌握其波动规律，对提高***的安全性和经济性有重要意义。

风电、光伏出力受自然天气条件的影响很大，当***中有大规模风能及光伏接入时，其出力的波动性会相对以往的火电、水电的出力调度有所不同。如何在满足***功率供需平衡的条件下，优先调度新能源，在考虑新能源波动性的情况下，让火电机组承担基本负荷，不同时段出力变化较小；让水电调节峰值负荷，不同时段间可以有较大的波动；同时考虑有功出力优化及无功出力优化且使***总发电费用最低，这些都对最优潮流的建模提出了更高的要求。电网最优潮流具有很高的实用价值，它第一次将经济性与安全性、有功和无功优化近乎完美的结合在一起，满足了大***互连、电网规模扩大后***规划设计人员、运行调度人员的要求。

由于新能源出力具有不确定性，概率最优潮流计算也变得更复杂。目前，已有相关文献对含新能源的最优潮流进行研究。文献《考虑注入功率分布的随机最优潮流方法》考虑风机出力的不确定性，建立了机会约束的最优潮流模型。文献《基于概率最优潮流的风电接入能力分析》运用随机技术的粒子群优化算法求解概率最优潮流模型，对风电接入能力的可行性和有效性进行了评估。但上述研究一般只考虑风电场，而很少研究风电场和光伏电站同时接入***对最优潮流的影响。风电场和光伏电站出力均具有随机性且出力特性不同，增加了电力市场中的不确定因素。

目前，考虑随机性影响的最优潮流计算方法主要包括蒙特卡洛法、累积量法、点估计法、蚁群算法等。蒙特卡洛法可以很好的研究随机性因素对***最优潮流的影响，但该方法需要成千上万次模拟***不同运行状态才能得到合理的结果，计算时间长、占用内存大。在输入随机变量相互独立或满足线性关系的前提下，累积量法用Gram-Charlier展开级数、Cornish-Fisher展开级数等进行拟合，从而得到输出随机变量的概率密度函数，提高了计算效率。点估计法虽然具有较快的计算速度，但其输出随机变量的高阶矩误差较大。蚁群算法运算量很大直接影响计算速度。

发明内容

为了解决上述问题，本发明提供一种分布式电源接入电网的潮流计算方法，包括以下步骤：

S1：读取电力***的初始数据；

S2：确定采样次数N和输入随机变量的维数s；

S3：按照以下3步，生成s×N阶采样矩阵，形成点列中第个点(j＝1,…,s；n＝1,…)的步骤如下：

S3-1：把第N-1个整数用2进制数表示，即式(1)

N-1＝a_R-1a_R-2…a₂a₁ (1)

其中a_n∈Z_b，Z_b＝{0,1,…，b-1}，R为满足b^r≤N的r的最大值；

S3-2：对N-1＝a_R-1a_R-2…a₂a₁进行排序，得到排序后的序列[d₁d₂…d_n…d_R]^T为式(2)

其中，为生成矩阵，0≤d_n≤b-1；引入生成矩阵是为了重置a₁a₂···a_n···a_R-1中各个数字的位置；数字的位置经过重置后，每一维和其它维的数字大小相同，但排列顺序不同，从而保证了结果的均匀性；

S3-3：经过第S3-2步的计算，可以表示为式(3)的2进制形式：

最后，将2进制表示的根据式(2)转化为10进制数即可；

S4：将采样次数初始化：令n＝1；

S5：判断n和采样次数N的大小，若n＞N，直接输出变量的概率统计结果；若n≤N，转S6；

S6：确定风电和光伏发电出力模型，确定负荷随机模型

S6-1：风速服从韦布尔分布，风电场有功功率P_w的概率密度函数可表示为式(4)：

式中：k,c分别为韦布尔分布的形状参数和尺度参数，P_r为风机额定功率，v_r,v_ci分别为额定风速和切入风速；

风电处理为PQ节点，令潮流计算中风机功率因数恒定不变，则无功功率按下式(5)计算：

式中：为功率因数角，对并网风机而言，一般位于第四象限，为负值。

S6-2：光伏出力随机模型

一定时间段内，太阳光照强度可认为服从贝塔分布，则光伏电站输出功率P_pv的概率密度函数表示为式(6)：

式中：R_pv＝Aηγ_max为仿真最大输出功率，A为太阳能电池仿真总面积，η为仿真总的光电转换效率，γ_max为一段时间内的最大光照强度，Γ为Gamma函数，α,β均为贝塔分布的形状参数；

与风电相同，潮流计算中将光伏电站也作为PQ节点；

S6-3：负荷随机模型

负荷具有时变性，很多有关文献都提出了对区域负荷进行预测的方法得到其概率分布；而作为中长期的负荷预测结果，负荷的概率分布规律基本符合于正态分布；其均值和方差均可以由大量的历史统计数据得到；这样，负荷的有功和无功功率的概率密度函数分别为式(7)和(8)：

式中：μ_p为有功功率的均值，δ_p ²为有功功率的方差，μ_Q为无功功率的均值，δ_Q ²为无功功率的方差；

S7：确定潮流计算模型

本发明建立各类能源有功及无功发电总费用最小为目标函数的最优潮流模型，尽可能调整发电机出力及无功源出力来满足负荷需要和***运行约束，在确保当前负荷需要以及满足各节点电压上下限和传输线路的传输功率极限情况下，寻找可行且总发电费用最小的发电机出力安排和电网潮流分布状态；

S7-1：目标函数

本发明构建的发电优化模型如下：

式(9)中C_Gpi、C_Gqi为机组i的有功和无功发电费用函数，C_gqj、C_gqj为无功补偿装置j的无功发电费用函数，P_Gi(t)、Q_Gi(t)为第i台发电机组在时段t的有功出力和无功出力，Q_gj(t)为第j台无功补偿装置在时段t的无功出力；N_g、N_q为发电机节点个数与无功补偿设备个数；目标函数是各个时段***的发电费用最小；

S7-2：等式约束

等式约束为各个时段的节点潮流平衡约束：

式(10)、式(11)中：V_i、θ_i为节点电压与相角，θ_ij＝θ_i-θ_j；P_Di、Q_Di为有功负荷与无功负荷；G_ij、B_ij为节点导纳矩阵的电导和电纳；

S7-3：不等式约束式(12)

式中，为发电机i有功出力上下限；为发电机i无功出力上下限；Q_gi为无功补偿设备i无功出力上下限；为节点电压幅值上下限；为线路ij持续输送容量极限(MVA)；N、N_b为节点集、支路集合；

P_GT,i(t+1)-P_GT,i≤R_i,up (13)

P_GT,i(t)-P_GT,i(t+1)≤R_i,down (14)

式(13)中，R_i，up为第i台火电机组的向上爬坡速率；式(14)中，R_i，down为第i台火电机组的向下爬坡速率；

S8：确定最优经济模型

S8-1：火电厂的发电费用

火电燃煤机组的有功出力是以煤耗量为标准进行计费的，机组i有功出力费用函数C_Gpi以式(15)进行计算。式中a_i、b_i、c_i为第i台火电机组的煤耗费用系数；

发电侧的无功电价分为两部分：无功容量电价和无功电量电价。无功电量电价主要涉及的是发电机的无功机会成本及有功损耗费用，本发明将无功机会成本作为发电机侧的总无功发电费用；

无功机会成本是该发电机因输出无功功率而损失的有功功率发电容量所对应的利润；如果忽略原动机的出力极限，并假设该无功机会成本C_op(Q_Gi)可表示如式(16)；

将式(15)代入到式(16)中，并进行泰勒展开，保留到项，忽略高次项后并整理得到式(17)；

C_Gqi(Q_Gi)为发电机组i的无功出力费用函数，S_Gi,max为发电机组i的额定视在功率，Q_Gi为发电机组i的无功出力值，k为发电厂生产有功功率的利润率，一般为5％-10％；

S8-2：水电厂的发电费用

目前我国水电运行成本一般是4～9分/千瓦时，而我国火电运行成本约为0.09-0.19元/千瓦时，本发明采用水电有功发电成本式(15)的形式进行计费，而其具体参数的取值近似与火电中a_i，b_i，c_i相差m倍，m为火电运行成本与水电运行成本电价的比值，a_i，b_i，c_i取值有微调变化，以区分同类电站的发电费用；水电厂的无功发电费用也采用类似火电厂无功出力费用，并按照式(16)的计费方式，其中的C_Gpi取相应水电厂的有功发电费用函数；

S8-3：光伏电站及风电场发电费用

目前光伏电站及风电场的上网电价仍高于传统能源，但是随着光伏设备和风电设备成本的降低，及国家针对新能源发电补贴政策的加强，光伏发电及风能发电的上网电价的进一步降低是可以预期的。本发明中以最大限度优先调用新能源为准则，令补贴后光伏发电费用及风力发电费用低于火电及水电的上网发电价格，其有功费用函数的选取与水电的有功费用选取方式相同；

S8-4：无功补偿设备的发电费用

以电容器、电抗器、同步调相机、SVC无功费用为固定成本表达式(18)：

其中Y为并联电容器的使用寿命，通常取15年；p为平均使用率，近似取为2/3，C_f为电容器单位容量的固定成本，平均可取为62500元/MVar，以此数据计算得出f_q＝1.97；

S9：潮流计算

利用拉格朗日函数法来处理优化问题中的等式约束，从而将具有等式约束的优化问题转化为无约束的优化问题；利用对数障碍函数法的罚函数方法处理不等式约束，最后用牛顿法来求解无约束优化问题最优解；

将非线性问题用以下数学公式表示：

obj min.f(x)

s.t.h(x)＝0 (19)

其中：min.f(x)为目标函数，是一个非线性函数；h(x)＝[h₁(x),...,h_m(x)]^T为非线性等式约束条件，g(x)＝[g₁(x),...,g_r(x)]^T为非线性不等式约束。假设在以上模型中共有k个变量，m个等式约束，r个不等式约束。用内点法求解问题(19)时，先将不等式约束转化为等式约束，同时构造障碍函数；为此先引入松弛变量l＞0，u＞0，l∈R^r，u∈R^r，将式(19)的不等式约束转化为等式约束，并把目标函数改造成障碍函数，可以得到以下优化问题A：

s.t.h(x)＝0 (11)

其中扰动因子u＞0；当l_i或u_i靠近边界时，以上函数趋于无穷大，因此满足以上障碍目标函数的极小解不可能在边界上找到，只能在满足l＞0，u>0时才可能得到最优解；这样，就通过目标函数的变换把含有不等式限制的优化问题变成了只含等式约束限制的优化问题A，因此可以直接用拉格朗日乘子法来求解。

优化模型A的拉格朗日函数为：

式中：y＝[y₁,...,y_m]^T，z＝[z₁,...,z_r]^T，w＝[w₁,...,w_r]^T均为拉格朗日乘子；该问题极小值存在的必要条件是拉格朗日函数对所有变量及乘子的偏导数为0，从而将有约束优化转化为无约束优化，接下来可以使用现有技术中的牛顿法求解；

S10：记录第n组节点电压、支路功率及发电成本等数据；

S11：进行下一轮潮流计算，t＝t+1，转S5；

本发明同现有技术相比，具有以下优点和有益效果：

1.本发明计算速度快、准确性高，得到的概率统计信息能全面地反应电力市场的运行状况，能有效地处理电力市场中的不确定性问题，具有较好的工程实用价值；

2.通过实施例仿真和验证，本发明可以提升分布式电源接入电网的安全经济运行，同时降低了网损，改善了节点电压水平，有效的提高了的经济性和实用性；

3.本发明采用的风电场和光伏电站混合***的节点电价、网损及支路功率波动情况比单独风电场***更小，接入***的光伏容量越大，节点电价越低，能够更加全面、有效、快捷地充分发挥潮流计算和控制的作用，本发明节约了调试的人力和物力，降低了的生产成本，有一定经济效益。

附图说明

图1是本发明的步骤流程图；

图2是实施例光伏接入电网时，节点电价的期望值；

图3是实施例光伏接入电网时，节点电价的标准差。

具体实施方式

一种分布式电源接入电网的潮流计算方法，包括以下步骤：

S1：读取电力***的初始数据；

S2：确定采样次数N和输入随机变量的维数s；

S3：按照以下3步，生成s×N阶采样矩阵，形成点列中第个点(j＝1,···,s；n＝1,···)的步骤如下：

S3-1：把第N-1个整数用2进制数表示，即式(1)

N-1＝a_R-1a_R-2···a₂a₁ (1)

其中a_n∈Z_b，Z_b＝{0,1,···,b-1}，R为满足b^r≤N的r的最大值；

S3-2：对N-1＝a_R-1a_R-2···a₂a₁进行排序，得到排序后的序列[d₁d₂···d_n···d_R]^T为式(2)

其中，为生成矩阵，0≤d_n≤b-1；引入生成矩阵是为了重置a₁a₂···a_n···a_R-1中各个数字的位置；数字的位置经过重置后，每一维和其它维的数字大小相同，但排列顺序不同，从而保证了结果的均匀性；

S3-3：经过第S3-2步的计算，可以表示为式(3)的2进制形式：

最后，将2进制表示的根据式(2)转化为10进制数即可；

S4：将采样次数初始化：令n＝1；

S5：判断n和采样次数N的大小，若n＞N，直接输出变量的概率统计结果；若n≤N，转S6；

S6：确定风电和光伏发电出力模型，确定负荷随机模型

S6-1：风速服从韦布尔分布，风电场有功功率P_w的概率密度函数可表示为式(4)：

式中：k,c分别为韦布尔分布的形状参数和尺度参数，P_r为风机额定功率，v_r,v_ci分别为额定风速和切入风速；

风电处理为PQ节点，令潮流计算中风机功率因数恒定不变，则无功功率按下式(5)计算：

式中：为功率因数角，对并网风机而言，一般位于第四象限，为负值；

S6-2：光伏出力随机模型

一定时间段内，太阳光照强度可认为服从贝塔分布，则光伏电站输出功率P_pv的概率密度函数表示为式(6)：

式中：R_pv＝Aηγ_max为仿真最大输出功率，A为太阳能电池仿真总面积，η为仿真总的光电转换效率，γ_max为一段时间内的最大光照强度，Γ为Gamma函数，α,β均为贝塔分布的形状参数；

与风电相同，潮流计算中将光伏电站也作为PQ节点；

S6-3：负荷随机模型

负荷具有时变性，很多有关文献都提出了对区域负荷进行预测的方法得到其概率分布；而作为中长期的负荷预测结果，负荷的概率分布规律基本符合于正态分布。其均值和方差均可以由大量的历史统计数据得到；这样，负荷的有功和无功功率的概率密度函数分别为式(7)和(8)：

式中：μ_p为有功功率的均值，δ_p ²为有功功率的方差，μ_Q为无功功率的均值，δ_Q ²为无功功率的方差；

S7：确定潮流计算模型

本发明建立各类能源有功及无功发电总费用最小为目标函数的最优潮流模型，尽可能调整发电机出力及无功源出力来满足负荷需要和***运行约束，在确保当前负荷需要以及满足各节点电压上下限和传输线路的传输功率极限情况下，寻找可行且总发电费用最小的发电机出力安排和电网潮流分布状态；

S7-1：目标函数

本发明构建的发电优化模型如下：

式(9)中C_Gpi、C_Gqi为机组i的有功和无功发电费用函数，C_gqj、C_gqj为无功补偿装置j的无功发电费用函数，P_Gi(t)、Q_Gi(t)为第i台发电机组在时段t的有功出力和无功出力，Q_gj(t)为第j台无功补偿装置在时段t的无功出力；N_g、N_q为发电机节点个数与无功补偿设备个数；目标函数是各个时段***的发电费用最小；

S7-2：等式约束

等式约束为各个时段的节点潮流平衡约束：

式(10)、式(11)中：V_i、θ_i为节点电压与相角，θ_ij＝θ_i-θ_j；P_Di、Q_Di为有功负荷与无功负荷；G_ij、B_ij为节点导纳矩阵的电导和电纳；

S7-3：不等式约束式(12)

式中，为发电机i有功出力上下限；为发电机i无功出力上下限；Q_gi为无功补偿设备i无功出力上下限；为节点电压幅值上下限；为线路ij持续输送容量极限(MVA)；N、N_b为节点集、支路集合；

P_GT,i(t+1)-P_GT,i≤R_i,up (13)

P_GT,i(t)-P_GT,i(t+1)≤R_i,down (14)

式(13)中，R_i，up为第i台火电机组的向上爬坡速率；式(14)中，R_i，down为第i台火电机组的向下爬坡速率；

S8：确定最优经济模型

S8-1：火电厂的发电费用

火电燃煤机组的有功出力是以煤耗量为标准进行计费的，机组i有功出力费用函数C_Gpi以式(15)进行计算。式中a_i、b_i、c_i为第i台火电机组的煤耗费用系数；

发电侧的无功电价分为两部分：无功容量电价和无功电量电价。无功电量电价主要涉及的是发电机的无功机会成本及有功损耗费用，本发明将无功机会成本作为发电机侧的总无功发电费用；

无功机会成本是该发电机因输出无功功率而损失的有功功率发电容量所对应的利润；如果忽略原动机的出力极限，并假设该无功机会成本C_op(Q_Gi)可表示如式(16)；

将式(15)代入到式(16)中，并进行泰勒展开，保留到项，忽略高次项后并整理得到式(17)；

C_Gqi(Q_Gi)为发电机组i的无功出力费用函数，S_Gi,max为发电机组i的额定视在功率，Q_Gi为发电机组i的无功出力值，k为发电厂生产有功功率的利润率，一般为5％-10％；

S8-2：水电厂的发电费用

目前我国水电运行成本一般是4～9分/千瓦时，而我国火电运行成本约为0.09-0.19元/千瓦时，本发明采用水电有功发电成本式(15)的形式进行计费，而其具体参数的取值近似与火电中a_i，b_i，c_i相差m倍，m为火电运行成本与水电运行成本电价的比值，a_i，b_i，c_i取值有微调变化，以区分同类电站的发电费用。水电厂的无功发电费用也采用类似火电厂无功出力费用，并按照式(16)的计费方式，其中的C_Gpi取相应水电厂的有功发电费用函数；

S8-3：光伏电站及风电场发电费用

目前光伏电站及风电场的上网电价仍高于传统能源，但是随着光伏设备和风电设备成本的降低，及国家针对新能源发电补贴政策的加强，光伏发电及风能发电的上网电价的进一步降低是可以预期的。本发明中以最大限度优先调用新能源为准则，令补贴后光伏发电费用及风力发电费用低于火电及水电的上网发电价格，其有功费用函数的选取与水电的有功费用选取方式相同；

S8-4：无功补偿设备的发电费用

以电容器、电抗器、同步调相机、SVC无功费用为固定成本表达式(18)：

其中Y为并联电容器的使用寿命，通常取15年；p为平均使用率，近似取为2/3，C_f为电容器单位容量的固定成本，平均可取为62500元/MVar，以此数据计算得出f_q＝1.97；

S9：潮流计算

利用拉格朗日函数法来处理优化问题中的等式约束，从而将具有等式约束的优化问题转化为无约束的优化问题；利用对数障碍函数法的罚函数方法处理不等式约束，最后用牛顿法来求解无约束优化问题最优解；

将非线性问题用以下数学公式表示：

obj min.f(x)

s.t.h(x)＝0 (19)

其中：min.f(x)为目标函数，是一个非线性函数；h(x)＝[h₁(x),...,h_m(x)]^T为非线性等式约束条件，g(x)＝[g₁(x),...,g_r(x)]^T为非线性不等式约束。假设在以上模型中共有k个变量，m个等式约束，r个不等式约束。用内点法求解问题(19)时，先将不等式约束转化为等式约束，同时构造障碍函数。为此先引入松弛变量l＞0，u＞0，l∈R^r，u∈R^r，将式(19)的不等式约束转化为等式约束，并把目标函数改造成障碍函数，可以得到以下优化问题A：

s.t.h(x)＝0 (11)

其中扰动因子u＞0。当l_i或u_i靠近边界时，以上函数趋于无穷大，因此满足以上障碍目标函数的极小解不可能在边界上找到，只能在满足l＞0，u>0时才可能得到最优解；这样，就通过目标函数的变换把含有不等式限制的优化问题变成了只含等式约束限制的优化问题A，因此可以直接用拉格朗日乘子法来求解。

优化模型A的拉格朗日函数为：

式中：y＝[y₁,...,y_m]^T，z＝[z₁,...,z_r]^T，w＝[w₁,...,w_r]^T均为拉格朗日乘子；该问题极小值存在的必要条件是拉格朗日函数对所有变量及乘子的偏导数为0，从而将有约束优化转化为无约束优化，接下来可以使用现有技术中的牛顿法求解；

S10：记录第n组节点电压、支路功率及发电成本等数据；

S11：进行下一轮潮流计算，t＝t+1，转S5。

本发明可以采用IEEE30节点实例进行仿真和验证。下面将结合实施例附图，对本发明的技术方案进行清楚、完整地描述。

本算例直接运用采样规模为512次的DN法分析***接入风电场和光伏电站后的最优潮流概率统计特性。

表1 节点电价期望值比较

为了分布式能源出力的随机性对节电电价的影响，下面分为2种情况进行讨论:案例1，对于8个节点20、61、104、123、138、171、198和207,每个节点均接入一个装机容量为60MW的风电场；案例2，对于8个节点20、61、104、123、138、171、198和207,每个节点均接入一个装机容量为30MW的风电场和一个装机容量为30MW的光伏电站。

运用本文介绍的算法进行概率最优潮流计算，得到分布式能源接入节点的节点电价，如表1所示。

通过比较两种情况下的计算结果可以看出，风电场和光伏混合***的节点电价相对于仅有风电场***要低。

表2为两种情况下***的网损期望值，风电场和光伏混合***的网损小于仅有风电场***的网损，可见风电场和光伏混合***更有利于***经济运行。

表2 网损期望值比较

分布式能源接入方式	网损/(MW)
		案例1	232.622
案例2	230.495

表3为两种情况下支路13～20和支路181～138的期望值和标准差。由表中结果可以看出，***单独接入风电场时的支路功率均值和标准差比风电场和光伏混合***要大，支路功率波动更大，出现重载和越限的概率也更大，不利于线路安全校核。

表3 支路功率比较

为了衡量不同容量的光伏接入***对概率最优潮流的影响，在节点20、61、104、123，138、171、198和207接入8个装机容量相等的光伏电站，8个光伏电站总装机从50MW依次增加到500MW(每次步进增加50MW)，每种容量下分别进行概率最优潮流计算。

图2和图3所示为不同容量光伏电站接入***时节点电价的期望值和标准差。图2结果表明,随着接入***光伏容量的不断增加，光伏节点的节点电价呈现降低的趋势。这是因为光伏出力可以代替一部分传统火电机组出力，从而使节点电价降低。图3结果表明，光伏出力的随机性及不确定性会带来节点电价的波动。

Claims (1)

1.一种分布式电源接入电网的潮流计算方法，其特征在于，包括以下步骤：

S1：读取电力***的初始数据；

S2：确定采样次数N和输入随机变量的维数s；

S3：按照以下3步，生成s×N阶采样矩阵，形成点列中第个点(j＝1,…,s；n＝1,…)的步骤如下：

S3-1：把第N-1个整数用2进制数表示，即式(1)

N-1＝a_R-1a_R-2…a₂a₁ (1)

其中a_n∈Z_b，Z_b＝{0,1,…,b-1}，R为满足b^r≤N的r的最大值；

S3-2：对N-1＝a_R-1a_R-2…a₂a₁进行排序，得到排序后的序列[d₁d₂…d_n…d_R]^T为式(2)

${\[d_1{d}_2...d_n...d_R\]}^T=C_{N\×N}^i{\[a_1{a}_2...a_n...a_{R-1}\]}^T---\left(2\right)$

其中，为生成矩阵，0≤d_n≤b-1；引入生成矩阵是为了重置a₁a₂…a_n…a_R-1中各个数字的位置；数字的位置经过重置后，每一维和其它维的数字大小相同，但排列顺序不同，从而保证了结果的均匀性；

S3-3：经过第S3-2步的计算，可以表示为式(3)的2进制形式：

$x_n^j=0.d_1{d}_2...d_R---\left(3\right)$

最后，将2进制表示的根据式(2)转化为10进制数即可；

S4：将采样次数初始化：令n＝1；

S5：判断n和采样次数N的大小，若n＞N，直接输出变量的概率统计结果；若n≤N，转S6；

S6：确定风电和光伏发电出力模型，确定负荷随机模型；

S6-1：风速服从韦布尔分布，风电场有功功率P_w的概率密度函数可表示为式(4)：

$f\left(P_w\right)=\frac{k}{k_{1}c}\left(\frac{P_w-k_2}{k_{1}c}\right)\exp \[-{\left(\frac{P_w-k_2}{k_{1}c}\right)}^k\]---\left(4\right)$

式中：k,c分别为韦布尔分布的形状参数和尺度参数，k₂＝-k₁v_ci，P_r为风机额定功率，v_r,v_ci分别为额定风速和切入风速；

风电处理为PQ节点，令潮流计算中风机功率因数恒定不变，则无功功率按下式(5)计算：

式中：为功率因数角，对并网风机而言，位于第四象限，为负值；

S6-2：光伏出力随机模型

一定时间段内，太阳光照强度可认为服从贝塔分布，则光伏电站输出功率P_pv的概率密度函数表示为式(6)：

$f\left(P_{pv}\right)=\frac{\mathrm{\Γ}(\mathrm{\α}+\mathrm{\β})}{\mathrm{\Γ}\left(\mathrm{\α}\right)+\mathrm{\Γ}\left(\mathrm{\β}\right)}{\left(\frac{P_{pv}}{R_{pv}}\right)}^{\mathrm{\α}-1}{(1-\frac{P_{pv}}{R_{pv}})}^{\mathrm{\β}-1}---\left(6\right)$

式中：R_pv＝Aηγ_max为仿真最大输出功率，A为太阳能电池仿真总面积，η为仿真总的光电转换效率，γ_max为一段时间内的最大光照强度，Γ为Gamma函数，α,β均为贝塔分布的形状参数；

与风电相同，潮流计算中将光伏电站也作为PQ节点；

S6-3：负荷随机模型

作为中长期的负荷预测结果，负荷的概率分布规律基本符合于正态分布；其均值和方差均可以由大量的历史统计数据得到；这样，负荷的有功和无功功率的概率密度函数分别为式(7)和(8)：

$f\left(P\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\δ}}_P}\exp (-\frac{{(P-{\mathrm{\μ}}_P)}^2}{2{{\mathrm{\δ}}_P}^2})---\left(7\right)$

$f\left(P\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\δ}}_Q}\exp (-\frac{{(Q-{\mathrm{\μ}}_Q)}^2}{2{{\mathrm{\δ}}_Q}^2})---\left(8\right)$

式中：μ_p为有功功率的均值，δ_p ²为有功功率的方差，μ_Q为无功功率的均值，δ_Q ²为无功功率的方差；

S7：确定潮流计算模型

S7-1：目标函数

构建的发电优化模型如下：

$\min F=\underset{i\∈Ng}{\Σ}\[C_{Gpi}\left(P_{Gi}\right(t\left)\right)+C_{Gqi}\left(Q_i\right(t\left)\right)\]+\underset{j\∈Nq}{\Σ}{C}_{gqj}\·\left(Q_{gj}\right(t\left)\right)---\left(9\right)$

式(9)中C_Gpi、C_Gqi为机组i的有功和无功发电费用函数，C_gqj为无功补偿装置j的无功发电费用函数，P_Gi(t)、Q_Gi(t)为第i台发电机组在时段t的有功出力和无功出力，Q_gj(t)为第j台无功补偿装置在时段t的无功出力；N_g、N_q为发电机节点个数与无功补偿设备个数；目标函数使各个时段***的发电费用最小；

S7-2：等式约束

等式约束为各个时段的节点潮流平衡约束：

$P_{Gi}-P_{Di}-V_i\underset{j=1}{\overset{N}{\Σ}}{V}_j(G_{ij}{\mathrm{cos\θ}}_{ij}+B_{ij}{\mathrm{sin\θ}}_{ij})=0---\left(10\right)$

$Q_{Gi}+Q_{gi}-Q_{Di}-V_i\underset{j=1}{\overset{N}{\Σ}}{V}_j(G_{ij}{\mathrm{sin\θ}}_{ij}+B_{ij}{\mathrm{cos\θ}}_{ij})=0---\left(11\right)$

式(10)、式(11)中：V_i、θ_i为节点电压与相角，θ_ij＝θ_i-θ_j；P_Di、Q_Di为有功负荷与无功负荷；G_ij、B_ij为节点导纳矩阵的电导和电纳；

S7-3：不等式约束式(12)

$\left\{\begin{array}{c}\underset{\‾}{P_{Gi}}\≤P_{Gi}\≤\stackrel{\‾}{P_{Gi}}\\ \underset{\‾}{Q_{Gi}}\≤Q_{Gi}\≤\stackrel{\‾}{Q_{Gi}}\\ \underset{\‾}{Q_{gi}}\≤Q_{gi}\≤\stackrel{\‾}{Q_{gi}}\\ \underset{\‾}{V_i}\≤V_i\≤\stackrel{\‾}{V_i}\\ S_i\≤\stackrel{\‾}{S_i}\end{array}\right.,\left\{\begin{array}{c}i\∈N_g\\ i\∈N_g\\ i\∈N_q\\ i\∈N\\ i\∈N_b\end{array}\right.---\left(12\right)$

式中， P _Gi 为发电机i有功出力上下限； Q _Gi 为发电机i无功出力上下限；Q_gi为无功补偿设备i无功出力上下限； V _i 为节点电压幅值上下限；为线路i持续输送容量极限(MVA)；N、N_b为节点集、支路集合；

P_GT,i(t+1)-P_GT,i≤R_i,up (13)

P_GT,i(t)-P_GT,i(t+1)≤R_i,down (14)

式(13)中，R_i，up为第i台火电机组的向上爬坡速率；式(14)中，R_i，down为第i台火电机组的向下爬坡速率；

S8：确定最优经济模型

S8-1：火电厂的发电费用

火电燃煤机组的有功出力是以煤耗量为标准进行计费的，机组i有功出力费用函数C_Gpi以式(15)进行计算；式中a_i、b_i、c_i为第i台火电机组的煤耗费用系数；

$C_{Gpi}=a_i{P}_{Gi}^2+b_i{P}_{Gi}+c_i---\left(15\right)$

无功机会成本是该发电机因输出无功功率而损失的有功功率发电容量所对应的利润；如果忽略原动机的出力极限，并假设该无功机会成本C_op(Q_Gi)可表示如式(16)；

$C_{Gqi}\left(Q_{Gi}\right)=C_{op}\left(Q_{Gi}\right)=k\[C_{Gpi}\left(S_{Gi,\max }\right)-C_{Gpi}\left(\sqrt{S_{Gi,\max }^2-Q_{Gi}^2}\right)\]---\left(16\right)$

将式(15)代入到式(16)中，并进行泰勒展开，保留到项，忽略高次项后并整理得到式(17)

$C_{Gqi}\left(Q_{Gi}\right)={\mathrm{dQ}}_{Gi}^2+e---\left(17\right)$

C_Gqi(Q_Gi)为发电机组i的无功出力费用函数，S_Gi,max为发电机组i的额定视在功率，Q_Gi为发电机组i的无功出力值，k为发电厂生产有功功率的利润率，一般为5％-10％；

S8-2：水电厂的发电费用

采用水电有功发电成本式(15)的形式进行计费，而其具体参数的取值近似与火电中a_i，b_i，c_i相差m倍，m为火电运行成本与水电运行成本电价的比值，a_i，b_i，c_i取值有微调变化，以区分同类电站的发电费用；水电厂的无功发电费用也采用类似火电厂无功出力费用，并按照式(16)的计费方式，其中的C_Gpi取相应水电厂的有功发电费用函数；

S8-3：光伏电站及风电场发电费用

以最大限度优先调用新能源为准则，令补贴后光伏发电费用及风力发电费用低于火电及水电的上网发电价格，其有功费用函数的选取与水电的有功费用选取方式相同；

S8-4：无功补偿设备的发电费用

以电容器、电抗器、同步调相机、SVC无功费用为固定成本表达式(18)：

$C_{gqj}\left(Q_{gj}\right)=\frac{C_f}{Y\×365\×24\×p}{Q}_{gj}=f_q{Q}_{gj}---\left(18\right)$

其中Y为并联电容器的使用寿命，取15年；p为平均使用率，近似取为2/3，C_f为电容器单位容量的固定成本，平均可取为62500元/MVar，以此数据计算得出f_q＝1.97；

S9：潮流计算

利用拉格朗日函数法来处理优化问题中的等式约束，从而将具有等式约束的优化问题转化为无约束的优化问题；利用对数障碍函数法的罚函数方法处理不等式约束，最后用牛顿法来求解无约束优化问题最优解；

将非线性问题用以下数学公式表示：

obj min.f(x)

s.t. h(x)＝0 (19)

$\underset{\‾}{g}\≤g\left(x\right)\≤\stackrel{\‾}{g}$

其中：min.f(x)为目标函数，是一个非线性函数；h(x)＝[h₁(x),...,h_m(x)]^T为非线性等式约束条件，g(x)＝[g₁(x),...,g_r(x)]^T为非线性不等式约束；假设在以上模型中共有k个变量，m个等式约束，r个不等式约束；用内点法求解问题(19)时，先将不等式约束转化为等式约束，同时构造障碍函数；为此先引入松弛变量l＞0，u＞0，l∈R^r，u∈R^r，将式(19)的不等式约束转化为等式约束，并把目标函数改造成障碍函数，可以得到以下优化问题A：

$\begin{array}{cc}obj& min.f\left(x\right)-\mathrm{\μ}(\underset{j=1}{\overset{r}{\Σ}}\ln l_j+\underset{j=1}{\overset{r}{\Σ}}\ln u_j)\end{array}$

s.t. h(x)＝0 (20)

$g\left(x\right)+u-\stackrel{\‾}{g}=0$

g(x)-l-g＝0

其中扰动因子u＞0；当l_i或u_i靠近边界时，以上函数趋于无穷大，因此满足以上障碍目标函数的极小解不可能在边界上找到，只能在满足l＞0，u>0时才可能得到最优解；这样，就通过目标函数的变换把含有不等式限制的优化问题变成了只含等式约束限制的优化问题A，因此可以直接用拉格朗日乘子法来求解；

优化模型A的拉格朗日函数为：

$L=f\left(x\right)-y^{T}h\left(x\right)-z^T\[g\left(x\right)-l-\underset{\‾}{g}\]-w^T\[g\left(x\right)+u-\stackrel{\‾}{g}\]-\mathrm{\μ}(\underset{j=1}{\overset{r}{\Σ}}\ln l_j+\underset{j=1}{\overset{r}{\Σ}}\ln u_j)---\left(21\right)$

式中：y＝[y₁,...,y_m]^T，z＝[z₁,...,z_r]^T，w＝[w₁,...,w_r]^T均为拉格朗日乘子；该问题极小值存在的必要条件是拉格朗日函数对所有变量及乘子的偏导数为0，从而将有约束优化转化为无约束优化，接下来可以使用现有技术中的牛顿法求解；

S10：记录第n组节点电压、支路功率及发电成本等数据；

S11：进行下一轮潮流计算，t＝t+1，转S5。