CN104009494B - 一种环境经济发电调度方法 - Google Patents

Abstract

本发明提供一种环境经济发电调度方法，包括建立考虑火电机组和水电机组出力平稳性的多目标调度模型，所述多目标调度模型包括目标函数和约束条件；输入预设的系统实时电力不足期望值；用多目标粒子群算法求解环境经济调度的帕累托最优前沿，获取历史帕累托最优解集；计算各帕累托最优解的满意度，取满意度最高的最优解作为本发明所需的环境经济发电调度结果。本发明可提供含运行风险水平约束的风-火-水协调优化调度方案，合理可行。

Description

一种环境经济发电调度方法

技术领域

本发明属于电力系统多目标优化调度技术领域，特别涉及基于多目标粒子群算法的环境经济发电调度方法。

背景技术

风电并网对于电网节能减排具有重要意义，但风电出力具有随机波动性，且现有出力预测方法存在较大的预测误差，这给电力系统调度带来严峻考验。为适应风电出力的波动性，常规火电机组或水电机组需要频繁调整机组出力，这会缩短机组使用寿命、增加调度员的工作量，影响系统安全。因此，制定合理的考虑风-火-水各自特性的调度方案显得尤为重要。

目前，国内外针对电力系统调度问题主要集中于传统经济调度和含风电场调度的风火联合调度，针对风-火-水协调调度的文献还较少。对于风电并网的处理也主要集中于风电出力不定性处理，但其预测误差固有不可消除，研究通过增减火电机组出力以适应风电出力的波动性和预测误差会不可避免的缩减常规机组使用寿命，增大系统运行风险。有文献指出，由于水电机组出力可以快速调整出力，提供系统备用，纳入水电可以提高系统纳入风电并网的能力。因此，有必要在制定调度方案时将水电与风电、火电一起考虑。电网现行水电调度方案是根据以往水库运行情况计划未来一日水电站出力，在制定调度方案时令未来一日的调度时段内水电一日总出力等于水电站计划出力。这种方法操作简单，且不会浪费水电资源。但是，随着风电并入电网，导致所需系统备用容量提高，而上述调度水电的方法会制约水电提供系统备用的能力，从而会降低系统可纳入风电的容量。因此，有必要制定更加节能的风-火-水协调调度方案。

发明内容

本发明的目的在于建立考虑火电、水电机组出力平稳性的风、火、水协调多目标动态调度模型，该模型在风电并网后的响应风险水平约束的前提下，不仅考虑火电、水电机组出力的平稳性，同时还考虑使得系统火电机组发电成本和排放，同时使水电机组出力尽量靠近计划出力，减少浪费。为求解该模型，通过引入教与学优化算法中的学习环节对传统多目标粒子群优化算法进行改进，提高算法的全局搜索能力。

本发明的技术方案提供一种环境经济发电调度方法，包括以下步骤：

步骤1、建立考虑火电机组和水电机组出力平稳性的多目标调度模型，所述多目标调度模型包括目标函数和约束条件；

目标函数如下，

$\mathrm{\ϵ}=\left(\underset{t=1}{\overset{\mathrm{NT}-1}{\mathrm{\Σ}}}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\right|P_{\mathrm{Gi}}(t+1)-P_{\mathrm{Gi}}\left(t\right)|+\underset{j=1}{\overset{\mathrm{Nh}}{\mathrm{\Σ}}}|P_{\mathrm{Hj}}(t+1)-P_{\mathrm{Hj}}\left(t\right)\left|\right)\right)/E$

F_P＝F_u+F_m+F_h

式中，ε为机组出力平稳性指标，E为某一日内参与调度的火电、水电机组的总装机容量；NT为某一调度日内所含调度时段，N为系统含火电机组数，P_Gi(t)为时段t火电机组i出力，i的取值为1,2,…,N；Nh为系统含水电机组数，P_Hj(t)为时段t水电机组j出力，j的取值为1,2,…,Nh；F_p为经济性指标目标函数值，F_u为能耗成本，F_m为废弃物造成的环境成本，F_h为偏离计划出力的惩罚项；

约束条件如下，

$\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{Gi}}\left(t\right)+\underset{j=1}{\overset{\mathrm{Nh}}{\mathrm{\Σ}}}{P}_{\mathrm{Hj}}\left(t\right)+\underset{k=1}{\overset{\mathrm{Nw}}{\mathrm{\Σ}}}{P}_{\mathrm{wk}}\left(t\right)=P_D\left(t\right)$

式中，Nw为系统含风电场个数，P_wk(t)为时段t风电场k出力，k的取值为1,2,…,Nw，P_D(t)为时段t系统负荷；

P_Gimin≤P_Gi(t)≤P_Gimax

式中，P_Gimin为火电机组i最小出力，P_Gimax为火电机组i最大出力；

P_Hjmin≤P_Hj(t)≤P_Hjmax

式中，P_Hjmin为水电机组j最小出力，P_Himax为水电机组j最大出力；

DR_Gi≤P_Gi(t+1)-P_Gi(t)≤UR_Gi

式中，DR_Gi为火电机组i的向下爬坡速率极值，UR_Gi为火电机组i向上爬坡速率极值；

$R_{\mathrm{EDNS}}^{\max }\≥\mathrm{REDNS}\left(t\right),t=1,...,\mathrm{NT}$

式中，为预设的调度时段t内系统实时电力不足期望值；

步骤2、输入预设的值；

步骤3、用多目标粒子群算法求解环境经济调度的帕累托最优前沿，获取历史帕累托最优解集；

步骤4、在步骤3得到历史帕累托最优解集基础上，计算各帕累托最优解的满意度，取满意度最高的最优解作为本发明所需的环境经济发电调度结果。

而且，步骤3包含以下子步骤，

步骤3.1、初始化及约束出力，包括初始相关参数和种群，令当前迭代次数为0，然后初始化全局帕累托最优解集及历史帕累托最优解集、独立备用选择集；

所述初始相关参数和种群，包括设置最大迭代次数K_max，初始化种群的速度和位置，初始化的种群中包括多个粒子，每个粒子代表一种出力方案，所述出力方案包括一个系统中所有机组分别在某一调度日内各调度时段的出力；对初始化种群结果进行约束处理，直至所产生的粒子满足约束；初始化粒子的个体最优解为该粒子初始化所得位置对应的目标函数值，初始化种群的全局最优解通过随机选取一组初始化粒子得到；

所述全局帕累托最优解集是用于存储本次迭代产生的非支配解的解集，历史帕累托最优解集是用于存储在整个迭代过程中所产生的非支配解的解集，独立备用选择集用于存储当前种群中的可行解；

步骤3.2、迭代次数q加1；

步骤3.3、调整惯性权重和粒子速度及位置；

步骤3.4、对本次迭代执行步骤3.3更新后的粒子进行约束处理；

步骤3.5、通过计算目标函数值获取各个粒子的适应度值；

步骤3.6、寻找非支配解，然后更新全局帕累托最优解集，根据全局帕累托最优解集画出帕累托最优前沿，并将非支配解存入历史帕累托最优解集中；

步骤3.7、使非支配解进行学习，包括对本次迭代执行步骤3.6所得帕累托最优前沿中的帕累托前沿解进行变异，对全局帕累托最优解集进行修整；

步骤3.8、对当前的全局帕累托最优解集进行约束处理，然后使满足约束的粒子与其他帕累托最优解中的粒子进行支配关系判断，得到新的全局帕累托最优解集，构建新的帕累托前沿；

再一次更新全局帕累托最优解集，并将非支配解存入历史帕累托最优解集中，将本次迭代执行步骤3.3～3.8中产生的非支配解全部存入备用选择集中；

步骤3.9、基于占优定义和K最邻近距离和更新全局帕累托最优解集和历史帕累托最优解集，包括在独立备用选择集和历史帕累托最优解集中，分别计算各个粒子的K最邻近距离和，并判断是否超出相应个数限制，当超出限制时根据K最邻近距离和由小到大进行修剪，直至满足相应个数限制；

步骤3.10、更新全局最优解和个体最优解；

步骤3.11、判断是否达到最大迭代次数K_max，若否则转步骤3.2，若是则根据全局帕累托最优解集输出帕累托最优前沿，保存历史帕累托最优解集。

而且，步骤3.7中，对帕累托前沿解进行变异实现如下，

设某一帕累托前沿解相应粒子为X＝[x₁,x₂,...,x_n]，则对其通过学习操作得到新的粒子，用新的粒子取代原来的粒子，学习操作见下式，

$x_{i^{\′}}=x_{i^{\′}}+\frac{1}{3}\×p(\mathrm{rand}1\×x_{i^{\′}-1}+\mathrm{rand}2\×x_{i^{\′}}+\mathrm{rand}3\×x_{i^{\′}+1})(K_{\max }-k)/K_{\max }$

式中，p为初始学习因子，K_max为最大迭代次数，k为当前迭代次数，rand1、rand2、rand3为（0,1）均匀分布随机数，分别为该帕累托最优解在帕累托前沿上左右相邻的粒子。

而且，步骤3.11所述根据全局帕累托最优解集输出帕累托最优前沿，包括计算全局帕累托最优解集各个粒子的K最邻近距离和，并判断是否超出预设的全局帕累托最优解集个数限制，当超出限制时根据K最邻近距离和由小到大进行修剪，直至满足全局帕累托最优解集个数限制。

而且，步骤3.9和步骤3.11中，计算K最邻近距离和时，求取两个粒子的距离d_i,j如下式实现，

$d_{i^{\′},j^{\′}}=\sqrt{\underset{k^{\′}=1}{\overset{N_{\mathrm{obj}}}{\mathrm{\Σ}}}{\left(\frac{f_k^{i^{\′}}-f_k^{j^{\′}}}{f_{k^{\′}}^{\max }-f_{k^{\′}}^{\min }}\right)}^2}$

式中，N_obj为目标函数的个数，和分别为第i'个粒子和第j'个粒子的第k'维目标函数值，和分别为第k'维目标函数值的最大值和最小值。

而且，步骤3.1、3.4和3.8的约束处理中，

对于等式约束，采用减少维数法将等式约束转变为不等式约束；

对于不等式约束，采用倒退机制进行处理，包括若某一时段粒子不满足不等式约束，则倒退两个时段重新寻找。

而且，步骤4采用折中解方案时，采用模糊隶属函数法来评估各帕累托最优解的满意度如下，

第i'个最优解在第k'维目标函数下的满意度如下，

${\mathrm{\μ}}_{i^{\′}}^{k^{\′}}=\left\{\begin{array}{cc}1,& f_{i^{\′}}^{k^{\′}}\≤f^{k^{\′},\min }\\ \frac{f^{k^{\′},\max }-f_i^{k^{\′}}}{f^{k^{\′},\max }-f^{k^{\′},\min }},& f^{k^{\′},\min }\≤f_{i^{\′}}^{k^{\′}}\≤f^{k^{\′},\max }\\ 0,& f_{i^{\′}}^{k^{\′}}\≥f^{k^{\′},\max }\end{array}\right.$

各个粒子的满意度如下

${\mathrm{\μ}}_{i^{\′}}=\underset{k^{\′}=1}{\overset{N_{\mathrm{obj}}}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i^{\′}}^{k^{\′}}$

式中，N_obj为目标函数的个数，为第i'个最优解的第k'维目标函数值，分别为第k'维目标函数值的最小值和最大值。

与现有技术相比，本发明具有以下有益效果：

1、本发明通过定义系统实时响应风险水平约束，相比较与传统确定性确定系统备用容量约束的方法，可以使系统在整个调度周期内的不同时段响应风险确定在同一水平。此外，仿真结果表明，采用该约束取代确定性方法确定系统备用，可以提供指标更优的调度方案。

2、本发明提在模型中在目标函数1中通过定义机组出力平稳性指标以衡量系统所包含的火电、水电机组出力平稳性，从而可以有效降低火电水电机组的出力波动性，减少调度员的工作量，有效防止因为频繁调整机组出力导致机组使用寿命缩减。在目标函数2中，不仅包含火电机组发电成本和排放物所导致的环境污染，还包含水电偏离计划出力所带来的浪费惩罚量，从而可以使得系统中火电、水电和风电协调优化调度。

3、为求解所提出的模型，本发明进一步提出了一种改进多目标粒子群优化算法（NMPSO），与传统多目标粒子群算法相比较能够有效地跳出局部最优，找到更好的帕累托前沿；所提出的约束处理方法能够减少不断尝试罚函数中的惩罚因子所带来的工作量；所提出的构建独立最优解备选择集的做法能够有效避免因迭代次数的不断增加导致算法速度下降。

4、本发明提出的实时响应风险约束水平能够满足运营者不同响应风险水平要求的调度策略。

附图说明

图1为本发明实施例求解模型的NMPSO算法流程图。

图2为本发明实施例系统某一日24时段负荷值和风电场C出力预测曲线图。

图3为本发明实施例情形1时MPSO得到的帕累托前沿。

图4为本发明实施例情形1时NMPSO得到的帕累托前沿。

图5为本发明实施例情形2时MPSO得到的帕累托前沿。

图6为本发明实施例情形2时NMPSO得到的帕累托前沿。

图7为本发明实施例不同REDNS水平约束下的帕累托前沿和折中解。

具体实施方式

本发明提出的调度模型包含了两个目标函数：1.在风电全额上网的前提下，通过将水电偏离水电一日预测出力值作为惩罚项，结合火电机组发电成本和废气物造成的环境污染一起构成系统发电环境成本目标函数。2.以火电和水电机组出力的波动量构建机组出力平稳性目标函数，建立考虑机组出力平稳性的风-火-水协调调度模型。环境成本目标函数的构建旨在使火电发电成本最小、废弃物排放最低、同时使水电在未来一日调度方案中的发电量尽量接近机组的安排出力，废弃物和水电偏离机组的安排出力量通过价格罚因子和火电发电成本以统一量纲。机组出力平稳性目标函数旨在使火电、水电机组在运行时相邻调度时段出力波动尽量小。

考虑机组出力平稳性的风火水协调优化调度模型是一个多维、多约束、多目标、非凸的非线性规划问题，传统求解算法失效。智能算法具有隐含的并行搜索特性，可以用于求解该问题。近年来，多目标群搜索、差分进化、粒子群等智能算法纷纷用于求解多目标电力系统调度问题。多目标粒子算法由于算法简单、搜索能力强等特点，在求解高维、多约束、多目标非线性规划问题中得到了广泛的应用。但是传统多目标粒子群算法有易陷入局部最优解陷阱中，容易错过全局最优解，如何采取合适的方法避免陷入局部最优是提高多目标粒子群搜索能力的关键。

本发明提出的技术方案可采用计算机软件技术实现自动运行流程。以下结合附图和实施例详细说明本发明技术方案。

本发明实施例包括以下步骤：

步骤1、建立考虑火电、水电机组出力平稳性的风、火、水协调多目标动态调度模型

目标函数1：机组出力平稳性指标最小

本发明采用机组平稳运行指标来衡量机组出力平稳性，该指标越小则机组出力平稳性越好。

火电、水电机组的出力频繁调整不仅会影响机组使用寿命，同时还会增加调度员的工作量。因此，在制定调度方案时有必要使火电、水电机组出力变化尽量减小。本发明通过定义机组出力平稳性指标来衡量火电、水电出力在某一调度日内的调整程度。

系统某一调度日内火电、水电机组出力波动量ΔQ为：

$\mathrm{\ΔQ}=\underset{t=1}{\overset{\mathrm{NT}-1}{\mathrm{\Σ}}}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\right|P_{\mathrm{Gi}}(t+1)-P_{\mathrm{Gi}}\left(t\right)|+\underset{j=1}{\overset{\mathrm{Nh}}{\mathrm{\Σ}}}|P_{\mathrm{Hj}}(t+1)-P_{\mathrm{Hj}}\left(t\right)\left|\right)---\left(1\right)$

式中：NT为某一调度日内所含调度时段，实施例以小时为单位划分时段，NT为24；N为系统含火电机组数；P_Gi(t)为时段t火电机组i出力，i的取值为1,2,…,N；Nh为系统含水电机组数；P_Hj(t)为时段t水电机组j出力，j的取值为1,2,…,Nh。

机组出力波动量能够反映参与调度的火电、水电机组在调度周期内出力变化情况。但若两种调度方案中参与调度的机组数不一样，机组出力波动量则不能准确反映这两种调度方案中机组出力调整程度的不同。因此，为比较这两种调度方案中火电、水电机组调整程度，定义机组出力平稳性指标ε为：

ε＝ΔQ/E（2）

式中：E为某一日内参与调度的火电、水电机组的总装机容量；ΔQ为某一调度日内火电、水电机组出力波动量。

目标函数2：经济性指标最小

电网根据水库的运行情况计划未来一日用于发电的电量，现行的调度方式是使未来一日参与调度每个水电站出力严格等于日计划出力，这种方式简单易行，且不会造成水电浪费。但是，这种方式会限制水电提供系统备用，降低系统纳入风电的能力。因此，本发明对于水电的出力原则是尽量使水电满足日计划出力，对于浪费的水资源采用价格罚因子纳入环境成本目标函数中去。

风-火-水协调优化调度模型的经济性目标函数不仅应该包括机组能耗成本和污染气体排放所造成的环境成本，还应该反映水电机组偏离计划出力所造成浪费量。

经济性指标目标函数表达式为：

F_P＝F_u+F_m+F_h（3）

式中：F_p为经济性指标目标函数值，F_u为能耗成本，F_m为废弃物造成的环境成本，F_h为偏离计划出力的惩罚项。

1）火电机组能耗成本可表示为：

$F_u=\underset{t=1}{\overset{\mathrm{NT}}{\mathrm{\Σ}}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}(a_i{P}_{\mathrm{Gi}}^2\left(t\right)+b_j{P}_{\mathrm{Gi}}\left(t\right)+c_i)---\left(4\right)$

式中：a_i、b_i、c_i为机组能耗特性系数。

2）火电机组废弃物的环境成本表达式为：

$F_m=\underset{t=1}{\overset{\mathrm{NT}}{\mathrm{\Σ}}}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}C\left(t\right)E_i\left(t\right)---\left(5\right)$

$E_i\left(t\right)={\mathrm{\α}}_i{P}_{\mathrm{Gi}}^2\left(t\right)+{\mathrm{\β}}_i{P}_{\mathrm{Gi}}\left(t\right)+{\mathrm{\ξ}}_i\exp ({\mathrm{\η}}_i\·P_{\mathrm{Gi}}\left(t\right))---\left(6\right)$

式中：C(t)为时段t的价格罚因子，具体实施时可由本领域技术人员预先设置；E_i(t)为火电机组i在时段t的废弃物排放量；α_i、β_i、ξ_i、η_i分别为火电机组i废弃物排放特性系数。

3）水电偏离计划出力惩罚项表达式为：

$F_h=\underset{j=1}{\overset{\mathrm{Nh}}{\mathrm{\Σ}}}{C}_h\left(j\right)\·(\mathrm{NT}\·P_{\mathrm{Hj}}-\underset{t=1}{\overset{\mathrm{NT}}{\mathrm{\Σ}}}{P}_{\mathrm{Hj}}\left(t\right))---\left(7\right)$

式中：C_h(j)为水电机组j惩罚因子，实施时可采用该水电机组的上网电价；P_Hj(t)为时段t水电机组j出力；P_Hj为水电机组j的调度日每个时段计划出力。

约束条件

1）系统功率平衡

$\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{Gi}}\left(t\right)+\underset{j=1}{\overset{\mathrm{Nh}}{\mathrm{\Σ}}}{P}_{\mathrm{Hj}}\left(t\right)+\underset{k=1}{\overset{\mathrm{Nw}}{\mathrm{\Σ}}}{P}_{\mathrm{wk}}\left(t\right)=P_D\left(t\right)---\left(8\right)$

式中：Nw为风电场个数，P_wk(t)为时段t风电场k出力，k的取值为1,2,…,Nw，P_D(t)为时段t系统负荷。

2）火电机组出力约束

P_Gimin≤P_Gi(t)≤P_Gimax（9）

式中：P_Gimin为火电机组i最小出力，P_Gimax为火电机组i最大出力。

3）水电机组出力约束

P_Hjmin≤P_Hj(t)≤P_Hjmax（10）

式中：P_Hjmin为水电机组j最小出力，P_Himax为水电机组j最大出力。

4）火电机组爬坡速率约束

DR_Gi≤P_Gi(t+1)-P_Gi(t)≤UR_Gi（11）

式中，DR_Gi为火电机组i的向下爬坡速率极值，UR_Gi为火电机组i向上爬坡速率极值，P_Gi(t+1)为时段t+1火电机组i出力。

5）实时响应风险上限

$R_{\mathrm{EDNS}}^{\max }\≥\mathrm{REDNS}\left(t\right),t=1,...,\mathrm{NT}---\left(12\right)$

式中，为预设的调度时段t内系统REDNS（Real-timeExpectedDemandNotSupplied，实时电力不足期望）值。某调度时段t内实际的系统实时电力不足期望（REDNS）值REDNS(t)求解方式如下：

$\mathrm{EDNS}\left(t\right)=\underset{K\∈S}{\mathrm{\Σ}}{C}_K^t{P}_K---\left(13\right)$

REDNS(t)＝EDNS(t)/P_D(t)（14）

$C_K^t=\max (0,L),K\∈S---\left(15\right)$

$L=\underset{m\∈U}{\mathrm{\Σ}}{P}_{g,m}\left(t\right)+P_D\left(t\right)\×\mathrm{Ld}\%+P_{\mathrm{Wk}}\left(t\right)\×W\%-\underset{n\∈A}{\mathrm{\Σ}}{\mathrm{SR}}_n\left(t\right)---\left(16\right)$

式中：REDNS(t)为t时段的实时系统响应风险水平，S为投运机组状态集，设其中任一状态记为状态K；为t时段系统处于状态K时需要消减的负荷功率，P_K为t时段系统处于状态K的概率，P_D(t)为t时段系统负荷；L表示系统在状态K时系统备用不足量；P_g,m(t)为组合A内某受迫停运火电机组m在t时段出力；U为状态K时受迫停运机组组合，A为状态K时可用机组组合；Ld%为负荷预测误差，P_wk(t)为时段t风电场k出力，W%为风电场出力预测误差，本发明中Ld%、W%分别取30%、3%；SR_n(t)为在状态K时组合A内某机组n在响应时间段t可以释放的备用容量，由于水电机组一般不计及爬坡速率约束，因此，若是火电机组则按式（18）计算，若是水电机组，则按式（19）计算；

传统的系统备用容量约束为：

$\mathrm{SSR}\left(t\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{SR}}_i\left(t\right)+\underset{j=1}{\overset{\mathrm{Nh}}{\mathrm{\Σ}}}{\mathrm{SR}}_j\left(t\right)\≥L\left(t\right)\×\mathrm{Ld}\%+W\left(t\right)\×W\%---\left(17\right)$

SR_i(t)＝min(P_Gimax-P_Gi(t),UR_Gi)（18）

SR_j(t)＝(P_Hjmax-P_j(t))（19）

式中：SSR(t)为系统时段t备用需求；SR_i(t)为火电机组i时段t能提供备用；SR_j(t)为水电机组j时段t能提供备用；L(t)为系统t时段预测负荷；W(t)为风电场t时段总的预测出力。

最终，实施例建立考虑火电、水电机组出力平稳性的多目标调度模型如下：

$\left\{\begin{array}{c}\mathrm{\ϵ}=\left(\underset{t=1}{\overset{\mathrm{NT}-1}{\mathrm{\Σ}}}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\right|P_{\mathrm{Gi}}(t+1)-P_{\mathrm{Gi}}\left(t\right)|+\underset{j=1}{\overset{\mathrm{Nh}}{\mathrm{\Σ}}}|P_{\mathrm{Hj}}(t+1)-P_{\mathrm{Hj}}\left(t\right)\left|\right)\right)/E\\ F_P=F_u+F_m+F_h\\ s.t.\\ \underset{i=1}{\overset{N}{\mathrm{\Σ}}}{P}_{\mathrm{Gi}}\left(t\right)+\underset{j=1}{\overset{\mathrm{Nh}}{\mathrm{\Σ}}}{P}_{\mathrm{Hj}}\left(t\right)+\underset{k=1}{\overset{\mathrm{Nw}}{\mathrm{\Σ}}}{P}_{\mathrm{wk}}\left(t\right)=P_D\left(t\right)\\ P_{\mathrm{Gi}\min }\≤P_{\mathrm{Gi}}\left(t\right)\≤P_{\mathrm{Gi}\max }\\ P_{\mathrm{Hj}\min }\≤P_{\mathrm{Hj}}\left(t\right)\≤P_{\mathrm{Hj}\max }\\ {\mathrm{DR}}_{\mathrm{Gi}}\≤P_{\mathrm{Gi}}(t+1)-P_{\mathrm{Gi}}\left(t\right)\≤{\mathrm{UR}}_{\mathrm{Gi}}\\ R_{\mathrm{EDNS}}^{\max }\≥\mathrm{REDNS}\left(t\right),t=1,...,\mathrm{NT}\end{array}\right.$

步骤2、输入值，以作为约束。

具体实施时，可由调度员根据实际情况和相关规定预先给定调度值。

步骤3、用多目标粒子群算法求解环境经济调度的帕累托最优前沿，获取历史帕累托最优解集。具体实施时，可以采用传统的多目标粒子群算法。为进一步提高效果起见，实施例提出改进多目标粒子群算法，算法流程如附图1所示，具体步骤如下：

步骤3.1、初始化机组出力。

首先，初始化相关参数和种群，令当前迭代次数q为0，然后初始化独立备用选择集合和帕累托最优解集。包括初始化种群的速度和位置，设置最大迭代次数K_max，给定学习因子C₁、C₂，初始化的种群中包括多个粒子，实施例中用每个粒子代表一种出力方案，即一个帕累托可行解，包括一个系统中所有机组分别在某一调度日内各调度时段的出力。

可根据本发明以下方法对对初始化种群结果处理约束，直至初始化种群所产生的粒子满足约束：

对于等式约束（8）：采用减少维数法进行处理，具体操作见步骤3.4；

对于约束（9）-（10），在产生初始粒子时，由于rand在（0,1）之间，这就能够保证粒子初始值满足约束；

对于约束（11），由于当前为初始迭代次数，不存在爬坡约束问题，在后续迭代过程中；

对于约束（12），重复产生初始粒子初始值的过程，直到满足该约束。

需要说明的是：尽管本发明提出的出力约束（12）的方法需要不断尝试，在该步骤增加了程序运行时间，但是本发明相比较于普遍采用罚函数法处理约束省去了寻找合适罚函数的时间，所以总的运行时间要短于罚函数法。

再根据目标函数评价初始种群，初始化各粒子的个体最优解、种群的全局最优解、全局帕累托最优解集及历史帕累托最优解集，令当前迭代次数q＝0。

本发明的初始化种群速度、位置、各粒子的个体最优解、种群全局最优解、全局帕累托最优解集及历史帕累托最优解集、独立备用选择集详细操作阐述如下：

初始化粒子位置：其中表示粒子i'初始化值（即第0次迭代）； 为该粒子上下限值，分别对应机组出力上下限值，rand表示随机数值，范围为（0,1）；可得到粒子i'在第z维空间中的初始飞行速度和位置

初始化种群速度： $x_{i^{\′}}^0=\mathrm{rand}\×0.005\×({\mathrm{ub}}^{i^{\′}}-{\mathrm{lb}}^{i^{\′}});$

初始化粒子的个体最优解：本发明取该粒子初始化值对应的根据（2）和（3）计算的目标函数值；可得粒子i'的个体最优值在第z维空间中的初始位置分量

初始化种群的全局最优解：随机选取一组初始化粒子作为全局最优解；

初始化全局帕累托最优解集及历史帕累托最优解集：全局帕累托最优解集是本发明设置的一个用于存储帕累托最优前沿可行解（即本次迭代产生的非支配解）的解集，每次一迭代结束时就是根据全局帕累托最优解集画出帕累托最优前沿；历史帕累托最优解集是用于存储在整个迭代过程中所产生的所有帕累托最优可行解集（即整个迭代过程中产生的非支配解）；本发明中全局帕累托最优解集和历史帕累托最优解集初始值全部随机取初始化产生的粒子。

初始化独立备用选择集：本发明构建独立备用选择集及修剪策略，单独构建备用选择集来存储当前种群中的可行解。独立备用选择集初始化为空，随着迭代，独立备用选择集中粒子个数会增加，计算速度会变慢，有必要对选择集的个体数加以限定（本专利在具体实施中取60），对超过限制值情况下的粒子采用基于K最邻近距离和进行修剪。个体的K最邻近距离和能够反映选择集的密度信息。

修剪原则为：若选择集的个体数大于设定值，则删除K最邻近距离和最小的若干个体，从而限定选择集的个体数，并保证选择集分布的均匀性。

K最邻近距离和计算：本发明提出对计算距离粒子最近的K个个体距离之和的具体实现进行改进，包括如式（20）所示求取两个粒子的距离对某粒子的距离最近的K个粒子求和，使其既能反应该粒子在种群中的相对位置，也能反应最邻近K个个体的分布情况。

$d_{i,j}=\sqrt{\underset{k^{\′}=1}{\overset{N_{\mathrm{obj}}}{\mathrm{\Σ}}}{\left(\frac{f_k^{i^{\′}}-f_k^{j^{\′}}}{f_{k^{\′}}^{\max }-f_{k^{\′}}^{\min }}\right)}^2}---\left(20\right)$

式中，N_obj为目标函数的个数，本发明值为2，k'=1,2；和分别为第i'个粒子和第j'个粒子的第k'维目标函数值，和分别为第k'维目标函数值的最大值和最小值。K的取值可预先设定，实施例取5。

步骤3.2、迭代次数q加1。

步骤3.3、调整惯性权重和粒子速度及位置，包括根据式（21）计算本次迭代的惯性权重w，采用式（22）更新粒子速度和位置。

$w=w_{\max }-\frac{w_{\max }-w_{\min }}{K_{\max }}q---\left(21\right)$

式中，K_max为最大迭代次数，q为当前迭代次数，惯性权重最大值w_max＝0.9，惯性权重最小值w_min＝0.4。

$\left\{\begin{array}{c}{v}_{i^{\′},z}^q=w\×v_{i,z}^{q-1}+c_1\×r_1\×(p_{i^{\′},z}^{q-1}-x_{i^{\′},z}^{q-1})+c_2\×r_2\×(g_z^{q-1}-x_{i^{\′},z}^{q-1})\\ x_{i^{\′},z}^q=x_{i^{\′},z}^{q-1}+v_{i^{\′},z}^q\end{array}---\left(22\right)\right.$

式中，r₁、r₂为[0,1]区间上服从均匀分布的随机数，分别为粒子i'在第q次迭代中在第z维空间中的飞行速度和位置，分别为粒子i'在第q-1次迭代中在第z维空间中的飞行速度和位置，为粒子i'的个体最优值在第q-1次迭代中在第z维空间中的位置分量，为种群的全局最优值在第q-1次迭代中在第z维空间中的位置分量；

第一次执行到步骤3.3时，q＝1，根据步骤3.1的初始化结果得到的采用式（22）更新得到q=2,3,…K_max时，同理根据第q-1次迭代所得相应数据求取。

步骤3.4、对本次迭代执行步骤3.3更新后的粒子约束处理。

本发明提出的新的约束处理方法：

对于等式约束（8），本发明采用减少维数法将等式约束转变为不等式约束，从而省去罚函数法不断尝试惩罚因子的工作，减少约束方法说明如下：

若x₁,x₂,...,x_n分别为某调度时段第1,2,…,n台机组出力，且要求满足等式约束：x₁+x₂+...+x_n＝P_D，P_D为该时段的系统负荷，则只需要初始化前n-1台机组出力，第n台机组出力为：x_n＝P_D-(x₁+x₂+...+x_n-1)，这样问题就由n维减少到n-1维；

对于不等式约束（10）～（12），传统的方法是采用罚函数法，这而确定罚函数因子需要大量的工作，本发明具体实施时可采用倒退机制进行处理。例如若某一时段粒子不满足不等式约束，则倒退两个时段，即当粒子在时段t的出力方案不满足不等式约束，则返回到时段t-2重新寻找。这样能够提高寻找到符合约束的粒子的概率，重新更新粒子位置和粒子群惯性权重，直至满足约束。

步骤3.5、计算各个粒子的适应度值：即根据式（2）、（3）计算目标函数值。

步骤3.6、寻找非支配解：

寻找支配解方法：对于所有粒子，按目标函数1、2求解目标函数值。根据目标函数值判断所有例子的支配关系。以粒子i为例，为其两个目标函数值，为粒子j的两个目标函数值，其中i≠j，若且或者且则称i支配j，j为支配解，若粒子i不被任何粒子支配，则将i称为非支配解。

根据上述方法得到帕累托非支配解后，更新全局帕累托最优解集，根据全局帕累托最优解集画出帕累托最优前沿，并将非支配解存入历史帕累托最优解集中。

步骤3.7、使非支配解进行学习，即根据式（23）对全局帕累托最优解集进行修整。

通过引入教与学算法中的学习环节，使独立备用选择集中的每一个粒子与其左右相邻的两个粒子进行学习，从而得到新的粒子。

变异的具体操作如下：

设本次迭代在步骤3.6根据全局帕累托最优解集画出的帕累托最优前沿中，某一帕累托前沿解包含的粒子为X＝[x₁,x₂,...,x_n]，n为该帕累托前沿解的粒子总数，则对每个粒子（i'＝1,2,...,n）通过学习操作得到新的粒子，用该粒子取代原来的粒子，操作见表达式（23）。

$x_{i^{\′}}=x_{i^{\′}}+\frac{1}{3}\×p(\mathrm{rand}1\×x_{i^{\′}-1}+\mathrm{rand}2\×x_{i^{\′}}+\mathrm{rand}3\×x_{i^{\′}+1})(K_{\max }-q)/K_{\max }---\left(23\right)$

式中：p为初始学习因子，本领域技术人员可自行预先设定，实施例根据实验结果预设为0.002；K_max为最大迭代次数；q为当前迭代次数；rand1、rand2、rand3为（0,1）均匀分布随机数；分别为该独立备用选择集中左右相邻的粒子。帕累托前沿当两端的粒子没有左或右边的相邻粒子时，可无需学习。

该操作粒子的学习环节不仅针对当前非支配粒子进行自学，还针对独立备用选择集中的粒子使其与其左右相邻的两个粒子进行学习。自学可以是粒子跳出局部最优陷阱，与左右粒子学习可以使新得到全局帕累托最优解集更加丰富。

步骤3.8、约束处理，采用步骤3.4中的方法处理约束。步骤3.1、3.4、3.8的约束处理方式一致。

对当前的全局帕累托最优解集进行约束处理，然后使满足约束的粒子与当前的全局帕累托最优解集中其他帕累托最优解中的粒子进行支配关系判断，当该粒子不是非支配解时，从当前的全局帕累托最优解集中删除，得到新的全局帕累托最优解集，从而构建新的帕累托前沿。进行支配关系判断可参见步骤3.6寻找非支配解的具体实现。

再一次更新全局帕累托最优解集，并将非支配解存入历史帕累托最优解集中。将本次迭代执行步骤3.3～3.8中产生的非支配解全部存入备用选择集中。

步骤3.9、基于占优定义和K最邻近距离和更新独立备用选择集及历史帕累托最优解集：在独立备用选择集和历史帕累托最优解集中，分别根据式（20）计算各个粒子的K最邻近距离和，并判断是否超出预设的独立备用选择集或历史帕累托最优解集个数限制，当超出限制时根据K最邻近距离和由小到大进行修剪，直至满足相应个数限制。

步骤3.10、更新全局最优解和个体最优解：将原来的全局最优解（个体最优解）和本次得到的全局最优解（个体最优解）进行比较，若小于本次全局最优解（个体最优解），则用本次全局最优解（个体最优解）替换原来的全局最优解（个体最优解）；若大于本次全局最优解（个体最优解），则不变。

步骤3.11、判断是否达到最大迭代次数K_max，若否，转步骤3.2；若是，根据全局帕累托最优解集输出帕累托最优前沿，保存历史帕累托最优解集。

此时全局帕累托最优解集中的粒子会很多，为了使得根据其输出的帕累托前沿更为光滑、均匀，具体实施时，可进一步根据K最邻近距离和对全局帕累托最优解集进行修改，根据修剪后的全局帕累托最优解集输出帕累托最优前沿。即根据式（20）计算全局帕累托最优解集各个粒子的K最邻近距离和，并判断是否超出预设的全局帕累托最优解集个数限制，当超出限制时根据K最邻近距离和由小到大进行修剪，直至满足全局帕累托最优解集个数限制。

步骤4、在步骤3得到的历史帕累托最优解集基础上，采用模糊数学法计算历史帕累托最优解集中各粒子的满意度，得到折中解方案。然后，取满意度最高的最优解作为本发明所需的环境经济发电调度结果。

实施例的步骤4所用折中解求解方案：采用模糊隶属函数法来评估历史帕累托最优解集中各帕累托最优解的满意度，第i'个最优解在第k'维目标函数下的满意度为：

${\mathrm{\μ}}_{i^{\′}}^{k^{\′}}=\left\{\begin{array}{cc}1,& f_{i^{\′}}^{k^{\′}}\≤f^{k^{\′},\min }\\ \frac{f^{k^{\′},\max }-f_i^{k^{\′}}}{f^{k^{\′},\max }-f^{k^{\′},\min }},& f^{k^{\′},\min }\≤f_{i^{\′}}^{k^{\′}}\≤f^{k^{\′},\max }\\ 0,& f_{i^{\′}}^{k^{\′}}\≥f^{k^{\′},\max }\end{array}\right.$

从而各个粒子的满意度为：

${\mathrm{\μ}}_{i^{\′}}=\underset{k^{\′}=1}{\overset{N_{\mathrm{obj}}}{\mathrm{\Σ}}}{\mathrm{\μ}}_{i^{\′}}^{k^{\′}}---\left(25\right)$

式中：k'∈{1,2,...,N_obj}，为第i'个最优解的第k'维目标函数值，分别为第k'维目标函数值的最小值和最大值。

为便于理解本发明技术效果起见，提供本实施例的应用例子如下：

1.考虑系统备用约束的算例设计与求解

在标准的10机测试系统中并入两个水电站A、B及一个风电场C，A、B的装机容量分别为150MW、80MW，风电场C的装机容量为150MW。水电机组相关参数如表1所示。火电机组废弃物价格罚因子如表2所示。系统负荷和风电场C预测出力值如表3和图2所示。

表1水电站A组和B组的参数

表2各时段的价格罚因子

表3系统某一日24时段负荷值和风电场C出力预测值

为了全面对比是否纳入水电的影响，设置了下面两个情形：

1）为计及系统备用约束的风-火联合调度，总装机容量为2631MW；

2）为计及系统备用约束的风-火-水联合调度模型，总装机容量为2861MW。

分别用多目标粒子群算法（Multi-objectiveParticleSwarmOptimization，MPSO和本发明提出的改进多目标粒子群算法（NewMulti-objectiveParticleSwarmOptimization，NMPSO两种算法独立计算25次并求均值。

情形1、2下得到的计算结果分别如表4和表5所示，帕累托（Pareto）前沿分别如图3-图6所示。情形2的折中解各台机组详细出力见附表A。

表4情形1时MPSO算法与NMPSO算法计算结果比较

表5情形2时MPSO算法与NMPSO算法计算结果比较

3实时等响应风险约束的风-火-水联合调度算例分析

各台火电机组的强迫停运概率见表6，按式（6）-式（9）求解的情形3中折中解各个时段内的实时响应风险水平如表7所示。

表6常规发电机组强迫停运的概率

表7情形2中NMPSO算法所得折中解各个时段对应的REDNS值

时段	1	2	3	4	5	6	7	8
									EDNSR(t)	0	0	0	0	0	0	4.62	18.3
时段	9	10	11	12	13	14	15	16
									EDNSR(t)	57.8	83.7	101.6	101.8	95.6	55.7	18.4	0
时段	17	18	19	20	21	22	23	24
									EDNSR(t)	0	7.04	27.5	67.9	48.5	2.42	0	0

分析表7可以发现，当采用风电预测出力比例确定系统备用容量时，时段10、11、12和时段13的运行风险较大，而且不同时间段的REDNS不相等。设系统在某一日的实时响应风险水平维持在某一相等水平并以此作为约束进行求解。

将REDNS约束分别设置为0.001018（情形2中NMPSO算法所得折中解中实时响应风险水平最高值）、0.0008、0.0002，采用NMPSO算法求解，结果如表8所示，帕累托前沿如图7所示。

表8不同REDNS约束下的解

由不同REDNS约束下的求解结果以及与情形3的对比分析，可以得出以下几点结论：

1）当REDNS设置为情形2中折中解对于的最大响应风险水平0.00108时，以环境成本最小为单目标时，环境成本为2608152.14$，比情形2中的2667196.13$减少了59043.99$。以机组出力平稳性指标最小为单目标时，平稳性指标为1.1794，比情形2减少了0.53%，对应的机组出力波动量为3374.36MW，比情形2中的3378.44MW减少了4.08MW。所求折中解两目标分别为2624450.73$、1.2145，比情形3的2682695.61$、1.2153分别减少了58244.88$、0.066%，对应的机组出力波动量减少了2.82MW。从而验证了以等实时响应为约束时能够获得比按比例型设定系统备用约束更小的环境成本和更好的机组出力波动量的调度方案。

2）比较不同REDNS约束下求解的结果可以得出结论，REDNS约束越低，环境成本越高和机组出力平稳性指标越大；当REDNS为0.00015时，无可行解。

对上述两种情形求解结果进行对比分析，可以得出以下几点结论：

1）不含水电时，以NMPSO算法求解时，以环境成本最小为单目标时为3045652.35$，含水电时为2621475.02$；以机组出力平稳性指标最小为单目标时为1.2854，对应的机组出力波动量为3381.83MW，含水电时为1.1857，对应的波动量为3392.41MW。不含水电时以传统方案得到的折中解两目标值分别为3047998.54$、1.2153，对应的波动量为3471.51MW；含水电时两目标值分别为2635780.97$、1.2153，对应的机组出力波动量为3477.11MW。

对比上述数据可以看出，纳入水电可以降低发电环境成本，能够提高机组出力平稳性，但由于水电机组也承担出力波动，因此总的机组出力波动量会增大。

2）比较表4、表5中MPSO算法和NMPSO算法得到的各自解，加入学习环节的NMPSO算法在情形1、2中均能找到比MPSO算法更优的解。

对比图3、图4、图5、图6，可以发现构建独立的备用选择集后的NMPSO算法所求得的帕累托更为均匀、光滑，从而验证了本发明所提的NMPSO算法能够有效地提高算法性能。

本发明中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (7)

1.一种环境经济发电调度方法，其特征在于，包括以下步骤：

步骤1、建立考虑火电机组和水电机组出力平稳性的多目标调度模型，所述多目标调度模型包括目标函数和约束条件；

目标函数如下，

$\mathrm{\ϵ}=\left(\underset{t=1}{\overset{NT-1}{\Σ}}\right(\underset{i=1}{\overset{N}{\Σ}}\left|P_{Gi}\right(t+1)-P_{Gi}(t\left)\right|+\underset{j=1}{\overset{Nh}{\Σ}}\left|P_{Hj}\right(t+1)-P_{Hj}(t\left)\right|))/E$

F_P＝F_u+F_m+F_h

式中，ε为机组出力平稳性指标，E为某一日内参与调度的火电、水电机组的总装机容量；NT为某一调度日内所含调度时段，N为系统含火电机组数，P_Gi(t)为时段t火电机组i出力，i的取值为1,2,…,N；Nh为系统含水电机组数，P_Hj(t)为时段t水电机组j出力，j的取值为1,2,…,Nh；F_p为经济性指标目标函数值，F_u为能耗成本，F_m为废弃物造成的环境成本，F_h为偏离计划出力的惩罚项；

约束条件如下，

$\underset{i=1}{\overset{N}{\Σ}}{P}_{Gi}\left(t\right)+\underset{j=1}{\overset{Nh}{\Σ}}{P}_{Hj}\left(t\right)+\underset{k=1}{\overset{Nw}{\Σ}}{P}_{wk}\left(t\right)=P_D\left(t\right)$

式中，Nw为系统含风电场个数，P_wk(t)为时段t风电场k出力，k的取值为1,2,…,Nw，P_D(t)为时段t系统负荷；

P_Gimin≤P_Gi(t)≤P_Gimax

式中，P_Gimin为火电机组i最小出力，P_Gimax为火电机组i最大出力；

P_Hjmin≤P_Hj(t)≤P_Hjmax

式中，P_Hjmin为水电机组j最小出力，P_Himax为水电机组j最大出力；

DR_Gi≤P_Gi(t+1)-P_Gi(t)≤UR_Gi

式中，DR_Gi为火电机组i的向下爬坡速率极值，UR_Gi为火电机组i向上爬坡速率极值；

$R_{EDNS}^{\max }\≥REDNS\left(t\right),t=1,...,NT$

式中，为预设的调度时段t内系统实时电力不足期望值，REDNS(t)为t时段的实时系统响应风险水平；

步骤2、输入预设的值；

步骤3、用多目标粒子群算法求解环境经济调度的帕累托最优前沿，获取历史帕累托最优解集；

步骤4、在步骤3得到历史帕累托最优解集基础上，计算各帕累托最优解的满意度，取满意度最高的最优解作为本发明所需的环境经济发电调度结果。

2.根据权利要求1所述环境经济发电调度方法，其特征在于：步骤3包含以下子步骤，

步骤3.1、初始化及约束出力，包括初始相关参数和种群，令当前迭代次数为0，然后初始化全局帕累托最优解集及历史帕累托最优解集、独立备用选择集，

所述初始相关参数和种群，包括设置最大迭代次数K_max，初始化种群的速度和位置，初始化的种群中包括多个粒子，每个粒子代表一种出力方案，所述出力方案包括一个系统中所有机组分别在某一调度日内各调度时段的出力；对初始化种群结果进行约束处理，直至所产生的粒子满足约束；初始化粒子的个体最优解为该粒子初始化所得位置对应的目标函数值，初始化种群的全局最优解通过随机选取一组初始化粒子得到；

所述全局帕累托最优解集是用于存储本次迭代产生的非支配解的解集，历史帕累托最优解集是用于存储在整个迭代过程中所产生的非支配解的解集，独立备用选择集用于存储当前种群中的可行解；

步骤3.2、迭代次数q加1；

步骤3.3、调整惯性权重和粒子速度及位置；

步骤3.4、对本次迭代执行步骤3.3更新后的粒子进行约束处理；

步骤3.5、通过计算目标函数值获取各个粒子的适应度值；

步骤3.6、寻找非支配解，然后更新全局帕累托最优解集，根据全局帕累托最优解集画出帕累托最优前沿，并将非支配解存入历史帕累托最优解集中；

步骤3.7、使非支配解进行学习，包括对本次迭代执行步骤3.6所得帕累托最优前沿中的帕累托前沿解进行变异，对全局帕累托最优解集进行修整；

步骤3.8、对当前的全局帕累托最优解集进行约束处理，然后使满足约束的粒子与其他帕累托最优解中的粒子进行支配关系判断，得到新的全局帕累托最优解集，构建新的帕累托前沿；

再一次更新全局帕累托最优解集，并将非支配解存入历史帕累托最优解集中，将本次迭代执行步骤3.3～3.8中产生的非支配解全部存入备用选择集中；

步骤3.9、基于占优定义和K最邻近距离和更新全局帕累托最优解集和历史帕累托最优解集，包括在独立备用选择集和历史帕累托最优解集中，分别计算各个粒子的K最邻近距离和，并判断是否超出相应个数限制，当超出限制时根据K最邻近距离和由小到大进行修剪，直至满足相应个数限制，K最邻近距离和为距离粒子最近的K个个体距离之和，K的取值预先设定；

步骤3.10、更新全局最优解和个体最优解；

步骤3.11、判断是否达到最大迭代次数K_max，若否则转步骤3.2，若是则根据全局帕累托最优解集输出帕累托最优前沿，保存历史帕累托最优解集。

3.根据权利要求2所述环境经济发电调度方法，其特征在于：步骤3.7中，对帕累托前沿解进行变异实现如下，

设某一帕累托前沿解相应粒子为X＝[x₁,x₂,...,x_n]，n为该帕累托前沿解的粒子总数，则对每个粒子x_i'通过学习操作得到新的粒子，i'＝1,2,...,n，用新的粒子取代原来的粒子，学习操作见下式，

$x_{i^{\′}}=x_{i^{\′}}+\frac{1}{3}\×p(rand1\×x_{i^{\′}-1}+rand2\×x_{l^{\′}}+rand3\×x_{i^{\′}+1})(K_{max}-k)/K_{max}$

式中，p为初始学习因子，K_max为最大迭代次数，k为当前迭代次数，rand1、rand2、rand3为(0,1)均匀分布随机数，x_i'-1、x_i'+1分别为该帕累托最优解在帕累托前沿上左右相邻的粒子。

4.根据权利要求2所述环境经济发电调度方法，其特征在于：步骤3.11所述根据全局帕累托最优解集输出帕累托最优前沿，包括计算全局帕累托最优解集各个粒子的K最邻近距离和，并判断是否超出预设的全局帕累托最优解集个数限制，当超出限制时根据K最邻近距离和由小到大进行修剪，直至满足全局帕累托最优解集个数限制。

5.根据权利要求4所述环境经济发电调度方法，其特征在于：步骤3.9和步骤3.11中，计算K最邻近距离和时，求取两个粒子的距离d_i',j'如下式实现，

$d_{i^{\′},j^{\′}}=\sqrt{\underset{k^{\′}=1}{\overset{N_{obj}}{\mathrm{\Σ}}}{\left(\frac{f_{k^{\′}}^{i^{\′}}-f_{k^{\′}}^{j^{\′}}}{f_{k^{\′}}^{\max }-f_{k^{\′}}^{min}}\right)}^2}$

式中，N_obj为目标函数的个数，和分别为第i'个粒子和第j'个粒子的第k'维目标函数值，和分别为第k'维目标函数值的最大值和最小值。

6.根据权利要求2或3或4或5所述环境经济发电调度方法，其特征在于：步骤3.1、3.4和3.8的约束处理中，

对于等式约束，采用减少维数法将等式约束转变为不等式约束；

对于不等式约束，采用倒退机制进行处理，包括若某一时段粒子不满足不等式约束，则倒退两个时段重新寻找。

7.根据权利要求2或3或4或5所述环境经济发电调度方法，其特征在于：步骤4采用折中解方案时，采用模糊隶属函数法来评估各帕累托最优解的满意度如下，

第i'个最优解在第k'维目标函数下的满意度如下，

${\mathrm{\μ}}_{i^{\′}}^k=\left\{\begin{array}{cc}1,& f_{k^{\′}}^{i^{\′}}\≤f_{k^{\′}}^{\min }\\ \frac{f_{k^{\′}}^{max}-f_{k^{\′}}^{i^{\′}}|}{f_{k^{\′}}^{\max }-f_{k^{\′}}^{\min }},& f_{k^{\′}}^{min}\≤f_{k^{\′}}^{i^{\′}}\≤f_{k^{\′}}^{max}\\ 0,& f_{k^{\′}}^{i^{\′}}\≥f_{k^{\′}}^{\max }\end{array}\right.$

各个粒子的满意度如下

${\mathrm{\μ}}_{i^{\′}}=\underset{k^{\′}=1}{\overset{N_{obj}}{\Σ}}{\mathrm{\μ}}_{i^{\′}}^{k^{\′}}$

式中，N_obj为目标函数的个数，为i'个粒子的第k'维目标函数值，表示第i'个最优解的第k'维目标函数值，和分别为第k'维目标函数值的最大值和最小值。