随着生产力水平的快速提高, 人们对电力的需求不断增加, 传统火力发电的成本快速上升, 也带来了环境污染问题^[1]。中国在减排方面需要发挥重要作用^[2], 在“碳达峰”和“碳中和”的背景下^[3-4], 科学决策机组组合问题, 充分利用资源, 既可以减少环境污染问题, 又可以实现经济效益的最大化。

机组组合问题被认为是一个NP-hard问题, 具有非凸性、非线性、高维等特点^[5]。目前解决机组组合问题的方法很多, 经典的方法有优先顺序法、动态规划法^[6-7]、拉格朗日松弛(Lagrangian Relaxation, LR)算法^[8]等。优先顺序法较为简单, 但得到的解的质量通常不是很好; 动态规划法易于保持解的可行性, 但容易陷入维度灾难, 优化求解时间过长; 智能算法包括粒子群(Particle Swarm Optimization, PSO)算法^[9-11]、遗传算法(Genetic Algorithm, GA)^[12-14]等。PSO算法具有较好的鲁棒性, 但搜索速度较慢, 容易陷入局部最优值, 有学者将PSO算法与其他进化算法相结合来解决优化问题^[15-18]。改进遗传算法使GA在机组组合问题上得到了进一步的发展, 但计算效率取决于参数的选择; 禁忌搜索不使用随机数, 不受随机数质量的影响, 但容易陷入局部最优。总体来说, 经典的方法比较简单, 但是从数学的角度看不够严谨, 或者容易造成维数灾难; 智能算法通常需要大量计算, 耗时较长。因此, 许多学者使用混合方法来解决机组组合问题, 但这些算法在具体应用中存在局限性, 场景不同, 算法效果会受到影响。大量文献表明, 拉格朗日松弛算法是解决NP-hard问题的有效算法, 在实际应用中, 会出现对偶间隙不能收敛的现象, 通过不断调整拉格朗日乘数, 可以在合理的时间内得到更好的近似最优解, 但乘子的更新路径会影响算法的有效性, 因此需要采取一些措施来缩小对偶间隙^[19]。

本文提出了一种改进的拉格朗日松弛算法LR-CMSCA。首先建立大规模的机组组合模型, 然后采用LR-CMSCA对模型进行优化求解, 最后进行数值仿真和结果分析。

本文考虑一个单目标优化问题。我们每时每刻都需要不同的电量, 为保证电力供需平衡, 需要提前规划好开机机组。为了解决停机问题, 需要将启动和停机成本包含在总发电成本T_c中, 即

$ T_{\mathrm{c}}=\min \sum\limits_{t=1}^T \sum\limits_{i=1}^N\left[O_{\mathrm{c} i}\left(P_i^t\right) u_i^t+R_{\mathrm{c} i}^t\left(1-u_i^{t-1}\right) u_i^t\right] $

(1)

式中: T——总时段数;

t——第t个时段, t=1, 2, 3, …, T;

O_ci(P_i^t)——机组i(i=1, 2, 3, …, N)在时间t的发电成本, 即消耗化石能源的燃料成本;

u_i^t——机组i在时间t的状态;

R_ci^t——机组i的启动成本。

当u_i^t =1时, 表示机组i在时间t运行; 当u_i^t =0时, 表示t时刻机组i处于停机状态。机组i的运行成本可以表示为以下二次函数

$ O_{\mathrm{c} i}\left(P_i^t\right)=a_i\left(P_i^t\right)^2+b_i P_i^t+c_i $

(2)

式中: a_i、b_i、c_i——机组i的煤炭消耗成本系数;

P_i^t——机组i在时间t的出力大小。

当机组设备关闭后, 将开始冷却。根据机组停机时间的长短, 启动可分为热启动和冷启动, 即

$ R_{\mathrm{ci}}^t=\left\{\begin{array}{l} H_{\mathrm{s} i}, M_{\mathrm{T} i, \mathrm{~d}} \leqslant T_{i, \mathrm{off}}^t \leqslant M_{\mathrm{T} i, \mathrm{~d}}+C_{\mathrm{T} i} \\ C_{\mathrm{s} i}, T_{i, \mathrm{off}}^t>M_{\mathrm{T} i, \mathrm{~d}}+C_{\mathrm{T} i} \end{array}\right. $

(3)

式中: H_si——机组i的热启动费用, 热启动费用相对较低;

C_si——机组i的冷启动费用, 冷启动费用相对较高;

M_{Ti, d}——机组i的最小连续停机时间;

T_{i, off}^t——机组i在t时刻的停机时长;

C_Ti——机组i的冷启动时间。

为了使模型更贴合实际, 需要考虑以下约束条件。

(1) 系统功率平衡约束。在任何时间段开启的所有机组的输出之和等于所需负载功率。

$ \sum\limits_{i=1}^N u_i^t P_i^t=P_{\mathrm{D}}^t $

(4)

式中: P_D^t——t时段所需的总功率, 即负载功率。

(2) 最小开停机约束。机组一般不能频繁启停, 需要满足最小启停时间。设T_{i, on}^t代表第i个机组在t时段的连续运行时间, M_{Ti, up}表示第i个机组的最小连续运行时间, 如果连续运行时间T_{i, on}^t小于最小连续运行时间M_{Ti, up}, 则该机组继续保持开机; 如果连续停机时间T_{i, off}^t小于最小连续停机时间M_{Ti, d}, 则该机组不能开机。

$ T_{i, \mathrm{on}}^t \geqslant M_{\mathrm{Ti}, \mathrm{up}}, T_{i, \mathrm{on}}^{t} \geqslant M_{\mathrm{Ti}, \mathrm{d}} $

(5)

(3) 机组出力上下限约束。机组在运行时需要满足自身出力大小约束, 即

$ P_{i, \min } \leqslant P_i^t \leqslant P_{i, \max } $

(6)

式中: P_i、min、P_{i, max}——第i个机组的最小和最大发电量。

(4) 系统旋转备用约束。电力的突然中断会带来许多不便, 甚至造成严重后果。为了应对一些突发状况, 如机组故障或者负荷大幅波动, 需要引入旋转备用来保障用电的安全和稳定。

$ \sum\limits_{i=1}^N u_i^t P_{i, \max } \geqslant P_D^t+S_{\mathrm{R} t} $

(7)

式中: S_Rt——系统在t时刻的旋转备用约束量, 在实验中取S_Rt=0.1P_D^t。

在众多算法中, LR算法因其计算速度快、适用于大规模混合整数规划问题而得到广泛应用。该算法引入了拉格朗日乘子, 形成了不包含耦合约束的原优化问题的对偶问题。这个问题对应的最优值(最小化)只是原问题最优值的一个下界, 因此需要调整拉格朗日乘子多次求解。拉格朗日乘子的选择和修正策略对算法的效率有一定的影响^[20-21]。

将系统约束条件加入目标函数中, 引入拉格朗日乘子λ_t和μ_t, 可得到拉格朗日函数

$ \begin{gathered} L_1=T_{\mathrm{c}}+\sum\limits_{t=1}^T \lambda_t\left(P_{\mathrm{D}}^t-\sum\limits_{i=1}^N u_i^t P_i^t\right)+ \\ \sum\limits_{t=1}^T \mu_t\left(P_{\mathrm{D}}^t+S_{\mathrm{R} t}-\sum\limits_{i=1}^N u_i^t P_{i, \max }\right) \end{gathered} $

(8)

式中: λ_t——t时段负荷平衡约束的松弛因子;

μ_t——t时段备用约束的松弛因子。

将式(1)代入式(8)可得:

$ \begin{gathered} L_1=\sum\limits_{t=1}^T \sum\limits_{i=1}^N\left[O_{\mathrm{c} i}\left(P_i^t\right) u_i^t+R_{\mathrm{c} i}^t\left(1-u_i^{t-1}\right) u_i^t-\right. \\ \left.\lambda_t u_i^t P_i^t-\mu_t u_i^t P_{i, \max }\right]+ \\ \sum\limits_{t=1}^T\left[\lambda_t P_{\mathrm{D}}^t+\mu_t\left(P_{\mathrm{D}}^t+S_{\mathrm{R} t}\right)\right] \end{gathered} $

(9)

采用拉格朗日松弛算法求解, 原函数的最优值J^*与对偶问题的最优值q^*之间的差异称为“对偶间隙”。本文采用相对对偶间隙Q^*来表示, 即

$ Q^*=\frac{J^*-q^*}{q^*} $

(10)

随着求解过程的进行, 对偶间隙逐渐缩小, 并且对偶间隙随着机组数目的增加而减小, 可将对偶间隙的大小作为是否停止迭代的准则。

正弦余弦算法(Sine Cosine Algorithm, SCA)由澳大利亚学者MIRJALILI S在2016年提出^[22]。与其他元启发式智能优化算法不同, SCA利用三角函数Sine和Cosine的周期性和波动性特征来实现算法的搜索和优化, 最终得到满足适应度要求的满意解。本文利用SCA的随机搜索性及其较强的全局搜索能力和局部开发能力, 将其与LR算法相结合, 优化拉格朗日乘子的迭代路径, 缩小相对对偶间隙, 提高可行解的精确度。

到目前为止, SCA及其许多改进算法已成功应用于优化领域^[23-24]。SCA在标准测试函数中其优化精度和收敛速度均优于PSO和GA等。

设X_i为第i个个体的位置, i=1, 2, 3, …, N, N为种群规模, f(X_i)为第i个个体的目标函数值, X_i^*为种群中最优个体的空间位置。在搜索过程中, 第i个个体在下一次迭代过程中空间位置的更新公式为

$ \begin{aligned} & X_i(t+1)= \\ & \left\{\begin{array}{l} X_i(t)+r_1 \times \sin r_2 \times \\ \left|r_3 \times X_i^*-X_i(t)\right|, \quad r_4<0.5 \\ X_i(t)+r_1 \times \cos r_2 \times \\ \left|r_3 \times X_i^*-X_i(t)\right|, \quad r_4 \geqslant 0.5 \end{array}\right. \end{aligned} $

(11)

$ r_1=a-\theta \frac{a}{\theta_{\max }} $

(12)

式中: X_i(t+1)——第i个个体第t+1次迭代解的位置;

X_i(t)——第i个个体第t次迭代解的位置;

r₁——线性递减函数, 描述了迭代时的步长, 大小根据式(12)进行调整;

r₂——[0, 2π]的随机数, 表示解的更新方向;

r₃——[0, 2]的随机权重, 表示当前最优解的影响力;

r₄——[0, 1]之间的随机数, 用于选择位置更新方程;

a——常数, 取a=2;

θ_max——最大迭代次数。

在SCA中, r₃是当前最优解的随机权重。为了能最大程度地保留优势粒子, 淘汰劣势粒子, 本文采用自适应权重更新策略来更新r₃。具体表达式为

$ r_3=\left\{\begin{array}{l} r_{\min }+\left(r_{\max }-r_{\min }\right) \frac{f_i-f_{\min }}{f_{\mathrm{avg}}-f_{\min }}, f \leqslant f_{\mathrm{avg}} \\ r_{\max }, f>f_{\mathrm{avg}} \end{array}\right. $

(13)

式中: r_min、r_max ——r₃的最小值与最大值;

f_i——粒子i的适应度值;

f_min——当前所有粒子的最小适应度值;

f_avg——当前所有粒子的平均适应度值。

对于适应度值大于平均适应度值的粒子, 应该有较大的权重, 以最大程度保留该粒子; 反之, 取较小的权重。

早有学者证明相比传统的高斯变异算子, 柯西变异具有更强的搜索能力。一维柯西分布的概率密度函数为

$ d(x)=\frac{1}{\pi} \cdot \frac{t}{t^2+x^2}, \quad-\infty<x<+\infty $

(14)

当t=1时为标准柯西分布。柯西分布函数为

$ D(x)=\frac{1}{\pi} \arctan \left(\frac{x}{t}\right)+\frac{1}{2} $

(15)

在原点处柯西分布的峰值小于高斯分布, 两端接近于零的速度也小于高斯分布, 使得个体有更高的概率逃离局部最优, 产生的后代与父本差异较大, 因此本文引入柯西变异来增加种群的多样性, 逃离局部最优解。

对粒子的状态X_i (t) 进行自适应变异, 即

$ X_i^{\prime}(t)=X_i(t)+X_i(t) C(0, 1) $

(16)

式中: C(0, 1)——柯西分布产生的随机值。

当粒子与当前最优粒子的位置距离越近, 表明种群间的相似性越大, 将距离当前最优粒子最远的个体进行柯西变异, 如果变异后粒子的适应度值大于原适应度值, 则变异后的粒子代替原粒子, 否则保持不变。

机组组合问题约束条件较多, 本身具有一定的难度, 当机组规模较大时, 则问题更具有挑战性。本文提出的LR-CMSCA算法总体步骤如下:

(1) 导入所有电力系统数据, 如最大和最小发电功率、启动成本、最小启停时间、成本系数等;

(2) 初始化参数及拉格朗日乘子;

(3) 判断是否达到最大迭代次数, 如果是, 则转到步骤(10), 否则进入步骤(4);

(4) 建立T个时段的启停问题及原问题的对偶问题;

(5) 求解对偶函数的最优解及各机组的出力大小;

(6) 判断是否为可行解, 如果不是, 转到步骤(7), 如果是, 转到步骤(8);

(7) 用SCA对拉格朗日乘子λ_t和μ_t进行调整, 并对距离当前最优解较远的解进行柯西变异, 将变异后的解与原解进行比较, 保留较好的解, 迭代次数加1, 返回步骤(3);

(8) 计算原问题的目标函数值;

(9) 判断对偶间隙是否满足要求, 如果满足, 转到步骤(10), 否则转到步骤(7);

(10) 算法结束。

图 1为LR-CMSCA算法流程示意。

为了验证本文所提算法的有效性, 将LR算法与LR-CMSCA在10机组系统中的相对对偶间隙变化进行比较, 如图 2所示。由图 2可知, 在LR算法的求解过程中, 相对对偶间隙先快速减小, 然后振荡, 不能收敛, 最终的解决方案未必可行; LR-CMSCA可以很好地改善振荡现象, 最终稳定收敛, 这主要得益于SCA的搜索随机性。

表 1为不同机组系统中相对对偶间隙变化比较。由表 1可知, 相同迭代次数时, 机组规模越大, 相对对偶间隙越小, 如10机组系统第8次迭代时相对对偶间隙是0.023, 而80机组系统第8次迭代时相对对偶间隙下降到了0.010, 表明机组规模越大精度越高, 收敛速度越快。

为了进一步验证LR-CMSCA的有效性, 本文分别对不同规模的火电机组进行仿真, 仿真在MATLAB 2020a环境下进行, 并与其他文献所提算法进行比较。在本案例中, 首先在一个典型的10机组系统中进行仿真, 10机组系统的总负荷需求如表 2所示, 特性参数如表 3所示, LR-CMSCA最终得到的10机组调度结果如表 4所示。表 4中, “1”表示投运, “0”表示关闭。10机组出力总值等于负荷总需求, 各机组出力结果如表 5所示, 出力平衡约束如图 3所示。调度时, 在满足负荷需求的情况下, 尽可能优先考虑成本较低的机组。因此, 由表 4可知, 1^#和2^#机组都是24 h投入的; 由表 5可知, 总成本为564 000 ＄。

为了避免偶然性, 本文将LR-CMSCA与文献[8]的LR算法、文献[13]的GA算法、文献[25]的ELRPSO算法、文献[22]的SCA算法在10、20、40、60、80、100机组系统中进行了仿真分析比较, 每种算法分别运行20次, 记录最优值, 仿真结果如表 6所示。由表 6可知, 与其他4种算法相比, 本文算法提供了较好的解决方案, 且机组规模越大效果越明显。这表明SCA是更新拉格朗日乘子的一种有效方法。通过引入柯西变异算子和自适应权重更新策略, 可以增加种群的多样性, 使粒子更快地向最优解移动。

针对LR算法在解决大规模机组组合时对偶间隙无法收敛的现象, 利用SCA随机搜索能力强的优点, 提出了一种新的改进算法LR-CMSCA, 并在不同规模机组系统中采用仿真实验测试该算法。实验结果表明, LR-CMSCA具有优越性, 传统LR算法对偶间隙不能收敛的现象得到抑制, 运行效率高, 且机组规模越大, 效果越好, 具有较好的应用前景。