最近在复现一篇EI论文题目是“基于粒子群算法的风电-水电抽水蓄能联合优化调度”。论文不算新但涉及风电、常规水电和抽水蓄能三种电源加上粒子群算法的实现内容非常典型。Matlab代码断断续续写了两周第一阶段算是完整跑通了。这篇博文把整个复现过程的思路、模型、代码结构和容易踩的坑全部记录下来给后面做同类研究的同学一个可以直接参考的路线图。我默认你手里的EI论文是英文或者中文核心问题建模基本围绕“如何让风电与抽水蓄能协同使系统运行成本最小或弃风率最低”展开。因为输入信息里没有给具体的目标函数表达式我采用的是一套在EI文献里非常常见的组合系统运行成本最小化加上弃风惩罚项。后面会解释为什么这样设计以及怎么替换成论文里的其他目标。1. 为什么风电要拉上抽水蓄能一起调度问题本质与EI复现的逻辑起点1.1 风电出力的“反调峰”困境风电出力最大的特点就是波动性和随机性。很多时候夜间负荷低风机却因为风速大而满发白天负荷高风速反而降下来。这种“反调峰”特性让电网调度员非常头疼。负荷高峰期需要风电顶上时它发不出来负荷低谷期又容易发得太多导致电网无法完全消纳只能弃风。在纯火电或纯水电系统里调度问题相对简单因为火电和水电机组的出力都可以主动调节。风电加入后问题变成了一个不确定环境下的决策问题在满足负荷需求、保证电网安全的前提下尽可能多地消纳风电同时让其他电源承担调节压力。这就是风电-水电联合优化调度的第一个逻辑起点风电出力必须先给定或预测出来然后再计算其他机组怎么配合。1.2 抽水蓄能在联合调度中扮演的角色抽水蓄能电站本质上是储能电站。它有两个工作状态抽水时消耗电能把下水库的水抽到上水库相当于给系统增加负荷发电时把上水库的水放下来推动机组发电相当于给系统增加电源。这个“削峰填谷”能力正好弥补风电的波动性。比如夜间风电大发时负荷又低抽水蓄能可以开抽水模式把多余的风电“吞掉”存成水的势能白天负荷高峰或者风电出力下降时再放水发电把能量“吐”出来。这样一来风电的弃风率能降下来系统的峰谷差也能被削平。这里要注意抽水蓄能电站不是常规水电站。常规水电站靠天然来水发电受河流来水量的约束抽水蓄能电站主要靠“水循环”工作抽上去多少水发下来多少电中间还要乘一个综合效率一般在0.75到0.8左右。这个效率因子特别容易在复现代码时写错后面我会专门讲。1.3 EI复现第一步先把论文的“系统拓扑”画出来很多同学拿到论文就直接看公式、抄代码结果跑了半天不知道每个变量代表什么。我的建议是不管EI论文里画没画系统图你都要自己先画一张“电源-负荷”框图。典型的风电-水电抽水蓄能联合系统长这样风电场出力曲线给定可弃风是一个可调减的电源。抽水蓄能电站有抽水工况和发电工况不能同时进行。常规水电站有入库径流有水库库容限制发电功率与发电流量相关。电网负荷每个时段的负荷需求给定是功率平衡约束的右边。可能还有一个火电机组作为基荷或调节电源目的是让系统供电可靠但会增加成本。把这几个部分画出来后你才知道决策变量有哪些约束条件有哪些。这一步看起来简单但能帮你避开“抄错变量含义”这种最蠢的错误。2. 联合优化调度的数学模型目标函数与约束条件的逐条拆解2.1 目标函数运行成本最小还是弃风量最小EI论文里目标函数可以是多种形式最常见的有系统总运行成本最小火电煤耗成本 启停成本 弃风惩罚等弃风量最小系统发电收益最大多目标加权经济性和环保性。我复现时采用的是典型的“运行成本最小化 弃风惩罚”形式。因为单纯弃风量最小没法体现经济性单纯成本最小又可能让系统把风电全弃掉——风电边际成本接近零所以必须加一个弃风惩罚项让算法尽量多用风电。可以写成min F sum_t ( C_f * P_f(t) ) lambda * sum_t ( P_wind_max(t) - P_wind(t) )其中P_f(t)是火电或外购电功率C_f是单位发电成本系数P_wind_max(t)是t时段风电最大可发功率P_wind(t)是t时段实际并网功率lambda是弃风惩罚系数一般取一个相对较大的数比如火电成本的5到10倍。如果你的论文没有火电只有风电、水电和抽蓄那目标函数可能会变成“抽水蓄能网损最小”或“水电出力最大化”但惩罚项思路是一样的。EI复现时先确认目标函数的形式再写代码。2.2 功率平衡与机组出力上下限约束功率平衡是所有电力调度模型必须满足的核心等式约束。每个时段系统总发电功率必须等于负荷功率。P_wind(t) P_hydro(t) P_ps_g(t) - P_ps_p(t) P_f(t) P_load(t)这里P_hydro(t)是常规水电出力P_ps_g(t)是抽蓄发电功率P_ps_p(t)是抽水功率P_load(t)是负荷。注意抽水功率是消耗电能的所以公式左边要减去它。很多初学者把抽水功率从等式右边减去最后功率平衡怎么都调不对。等式约束之外就是各类机组的出力上下限约束0 P_wind(t) P_wind_max(t) P_hydro_min P_hydro(t) P_hydro_max 0 P_ps_g(t) P_ps_g_max * u_g(t) 0 P_ps_p(t) P_ps_p_max * u_p(t) u_g(t) u_p(t) 1最后这个互斥约束特别重要抽水蓄能机组在同一个时段要么抽水要么发电绝对不能同时进行。这个约束如果写漏了算法很容易给出一个“既抽水又发电”的荒谬解从能量角度看等于让水在上水库和下水库之间白做功效率还不为0结果会非常失真。2.3 抽水蓄能电站的核心约束抽发循环的能量守恒抽水蓄能本质上是在“搬运”能量所以必须满足一个跨时段的水量/能量守恒约束。如果不加这个约束电站可以凭空发很多电系统功率平衡也能满足但那是不合理的。最常见的模型是用上水库的蓄能状态来描述S(t1) S(t) eta_p * P_ps_p(t) * delta_t - P_ps_g(t) * delta_t / eta_g这里S(t)是上水库t时段的蓄能单位MWh或m3水当量eta_p是抽水效率通常0.85左右eta_g是发电效率通常0.9左右delta_t是每个时段的时长单位小时24时段模型里通常为1。这个公式的意思是抽水时蓄能增加但因为有损耗所以乘以效率发电时蓄能减少因为同样有损耗所以除以效率。如果觉得两个效率相乘麻烦也可以用一个综合效率eta_total eta_p * eta_g公式变成S(t1) S(t) eta_total * P_ps_p(t) * delta_t - P_ps_g(t) * delta_t但它只能用于粗略估算严谨一点还是分开写。另外还需要上下水库容量约束S_min S(t) S_max S(1) S(24) S_init这个“调度周期始末蓄能相等”的约束非常关键否则算法会把上水库的水一次性放完得到一个不可持续的解。我在最初复现时漏了这个结果收敛曲线很漂亮但调度结果完全不合理——抽蓄电站最后库容见底这个解放到实际系统中根本没法执行。2.4 常规水电及水库水量平衡约束常规水电和抽水蓄能不一样它有一个天然来水过程。典型模型里水电出力与发电流量和水库水头相关。简单做法是用“水库蓄水量”作为状态变量V(t1) V(t) Q_in(t) - Q_hydro(t) - Q_spill(t)其中V(t)是水库库容Q_in(t)是天然入库流量Q_hydro(t)是发电流量Q_spill(t)是弃水流量。水电出力可以近似为P_hydro(t) k * Q_hydro(t) * H_net(t)如果假设水头恒定那么出力就和发电流量近似成正比变成一个线性关系。EI复现中很多论文都做了这个简化你不需要一上来就搞水头非线性模型先复现线性版本再考虑扩展。库容约束和初始库容约束同样要写V_min V(t) V_max V(1) V_init V(24) V_end如果论文里没给水库初始库容你可以设置一个合理值比如库容上限的50%并在代码里作为参数保留方便后面测试不同场景。3. 粒子群算法如何落到调度场景编码方式、适应度函数与参数整定3.1 为什么选粒子群而不是遗传算法或线性规划这个问题我复现前也纠结过。联合优化调度其实可以建模成混合整数线性规划用cplex或gurobi这类求解器也能解。但EI论文明确用了粒子群算法复现的目的就是要还原论文的方法。而且粒子群在非线性、非凸目标函数上有天然优势不需要求导不依赖初始点实现代码也短适合快速验证模型。和遗传算法相比粒子群收敛速度更快、参数更少但容易早熟。不过对于24时段、决策变量不算太多的调度问题粒子群的早熟问题可以通过增加种群规模、调节惯性权重来缓解。所以它成了很多电力调度论文的首选智能算法。3.2 粒子编码把24小时的调度序列排成一个向量粒子群算法里每个粒子代表一个候选解。在调度问题中候选解就是一整天的调度计划。怎么把调度计划编码成一个“位置向量”是写代码前必须想清楚的事情。我采用的编码方式非常直接位置向量 [所有时段的抽蓄发电功率P_ps_g(1..24), 所有时段的抽蓄抽水功率P_ps_p(1..24), 所有时段的常规水电出力P_hydro(1..24), 所有时段的并网风电功率P_wind(1..24)]假设调度周期有T个时段那么每个粒子的维度就是D 4 * T。如果T24D96。你可能要问为什么不把火电功率也作为决策变量因为火电功率可以由功率平衡约束直接算出来P_f(t) P_load(t) - P_wind(t) - P_hydro(t) - P_ps_g(t) P_ps_p(t)这是处理等式约束常用的技巧把等式约束变成一个“计算式”只要其余变量定了火电功率就唯一确定等式约束自动满足。这样粒子编码维度少了一个变量算法搜索空间也更小。位置向量的每个元素都有物理边界。比如风电并网功率不能超过该时段最大可发功率抽蓄功率不能超过额定功率。这些边界在初始化位置时就要设好后面迭代更新位置时也要做边界处理。3.3 适应度函数与惩罚函数设计粒子群算法需要用一个适应度值来衡量每个粒子的好坏。在最小化问题里适应度值越小越好。我直接把目标函数值作为基础再加上约束违反惩罚。为什么需要惩罚因为粒子更新过程中位置向量是随机变化的很可能不满足所有约束。比如抽水蓄能水库的容量约束可能被打破常规水电的库容可能超出上限。如果这些不满足约束的粒子直接参与比较搜索会被带偏。常用方法是“外点罚函数法”fitness F sigma1 * sum(约束违反量^2) sigma2 * ...具体到我的代码里约束包括抽蓄上水库蓄能上下限约束常规水电库容上下限约束抽蓄始末蓄能相等约束水电站始末库容相等约束抽水与发电状态互斥约束火电功率上下限约束用功率平衡等式算出的火电可能超限。每个约束违反的绝对值会乘一个惩罚系数然后平方累加。惩罚系数不能太小否则约束形同虚设也不能太大否则会主导适应度函数让算法只顾满足约束而忽略经济性。我的经验是先跑一次不惩罚的版本看目标函数量级然后把惩罚系数设为目标函数量级的10到100倍再逐步调整。3.4 参数整定的几个经验值粒子群算法的标准参数有四个种群规模N、最大迭代次数iter_max、惯性权重w、加速常数c1和c2。我复现时用的一组默认参数如下参数取值说明种群规模N100对于96维问题100个粒子不算大但已经能收敛最大迭代次数300后期看收敛曲线调整如果没平就加惯性权重w0.9到0.4线性递减前期全局搜索后期局部搜索学习因子c12.0个体认知权重学习因子c22.0社会经验权重速度上限Vmax每个变量范围的20%防止粒子飞得太远惯性权重递减是非常经典的策略公式是w w_max - (w_max - w_min) * iter / iter_max前期的w比较大粒子飞行速度快能覆盖整个搜索空间后期w变小粒子在局部精细搜索。这个策略对调度模型特别管用因为初始阶段很多粒子可能落在不可行域需要大范围探索才能找到可行区域。4. Matlab代码架构详解主程序、约束处理与循环迭代的骨架4.1 主程序流程从参数初始化到结果出图写Matlab代码时我习惯把整个复现过程分成逻辑清晰的三层主脚本、函数文件和绘图脚本。主脚本负责参数设置、初始化、调用优化循环、输出结果。函数文件包括目标函数、约束校验、粒子更新等。这样代码好读也好改。主程序的核心流程可以描述为设置基础数据负荷曲线、风电最大出力曲线、水电参数、抽蓄参数、算法参数。初始化粒子群随机生成N个位置向量和N个速度向量。计算初始适应度对每个粒子调用适应度函数更新个体最优和全局最优。迭代循环更新速度、更新位置、边界处理、计算适应度、更新pbest和gbest。输出最优解把gbest位置解码成调度计划绘制功率平衡图、库容曲线、收敛曲线。我给的代码骨架是经过简化的但核心逻辑和完整版一致。实际复现时你可以在每个循环中间加一个调试打印观察适应度值的变化趋势。主脚本开头通常是数据区我用结构体把系统参数存起来避免函数参数列表太长%% 参数设置 T 24; % 时段数 N 100; % 粒子群规模 iter_max 300; % 最大迭代次数 c1 2.0; c2 2.0; % 学习因子 w_max 0.9; w_min 0.4; % 负荷与风电预测数据示例 load_curve [ /* 24个负荷值 */ ]; wind_max [ /* 24个风电最大出力值 */ ]; % 抽水蓄能参数 P_ps_g_max 80; % 发电最大功率 P_ps_p_max 80; % 抽水最大功率 S_max 500; % 上水库最大蓄能 S_min 50; % 上水库最小蓄能 S_init 250; % 初始蓄能 eta_p 0.85; % 抽水效率 eta_g 0.90; % 发电效率 % 常规水电参数 P_hydro_max 120; P_hydro_min 20; V_max 1000; V_min 200; V_init 600; Q_in ones(1, T) * 20; % 入库流量 k_h 0.8; % 水电出力系数这些参数只代表一个测试算例实际复现时请根据论文中的测试系统填写。4.2 核心代码片段PSO速度与位置更新PSO更新是算法的核心。标准的速度更新公式v w * v c1 * rand(size(v)) .* (pbest - x) c2 * rand(size(v)) .* (gbest - x);位置更新x x v;速度限制v max(v, -Vmax); v min(v, Vmax);位置边界处理可以用“边界吸收”x(x lb) lb(x lb); x(x ub) ub(x ub);其中lb和ub是所有决策变量的上下界向量。以抽蓄发电功率为例下界是0上界是P_ps_g_max抽水功率同理。风电并网功率上界是wind_max中的对应值。完整的更新循环大概是这样的for iter 1:iter_max w w_max - (w_max - w_min) * (iter / iter_max); for i 1:N v(i,:) w * v(i,:) c1*rand(1,D).*(pbest(i,:) - x(i,:)) ... c2*rand(1,D).*(gbest - x(i,:)); v(i,:) max(v(i,:), -Vmax); v(i,:) min(v(i,:), Vmax); x(i,:) x(i,:) v(i,:); x(i,:) max(x(i,:), lb); x(i,:) min(x(i,:), ub); fit(i) fitness_func(x(i,:), sys_param); if fit(i) fit_pbest(i) pbest(i,:) x(i,:); fit_pbest(i) fit(i); end if fit(i) fit_gbest gbest x(i,:); fit_gbest fit(i); end end best_history(iter) fit_gbest; end这段代码看起来简单但有几个细节需要注意。rand(1,D)每次在速度更新公式里应该生成不同的随机数所以必须在表达式里直接调用不能提前保存一个固定向量否则算法会失去随机性。另一个细节是粒子的历史最优。在电力调度问题中即使当前迭代的适应度值不是很好历史最优里可能已经存了一个可行且经济的解。所以pbest的更新必须严格用“小于”而不是“小于等于”否则可能出现历史最优被一个“更差但数值上等于”的解覆盖的问题。4.3 约束处理的实现细节边界吸收与状态互斥边界吸收是处理上下界约束最简单的方法。位置更新后如果某个决策变量超过上界就直接拉回上界低于下界就拉回下界。这个方法的好处是实现简单不容易产生意料之外的解。但缺点是会让靠近边界的粒子失去速度信息可能降低搜索效率。另一种方法是“随机重置到边界附近”实现稍微复杂一点但效果更好。抽水与发电状态互斥约束不能靠边界处理因为它是两个变量之间的耦合关系。我写了一个专门的处理函数如果一个粒子里同一时段的P_ps_g(t)和P_ps_p(t)都大于0就按一个随机偏向保留其中一个另一个置0。比如生成一个随机数如果大于0.5就保留发电、清掉抽水否则保留抽水、清掉发电。这样能保证互斥约束成立。不过要注意这种“修复策略”在每一代都做等于人为改变了粒子的位置。如果你的惩罚函数设计得好也可以不强制修复而是让算法通过惩罚自动学会“不要同时抽水和发电”。我个人的经验是先强制修复保证每个粒子都满足互斥约束这样跑出来的结果物理上一定合理对排查bug也有帮助。4.4 代码组织建议脚本、函数与测试用例完整代码里我建议至少包含这几个文件main_PSO_optimization.m主脚本fitness_func.m适应度函数目标 惩罚decode_solution.m把位置向量解码成调度计划constraint_fix.m互斥约束修复、边界处理plot_results.m出图脚本。函数模块化的好处是出了问题可以单点调试。我最开始把所有逻辑都写在一个脚本里出了bug定位很慢。后来拆成函数通过打印中间变量很快就发现是抽蓄能量守恒公式里的效率除反了。5. 复现结果怎么判断对不对收敛曲线、调度出力与几个容易翻车的细节5.1 收敛曲线的判读标准不管用什么优化算法第一步都要画收敛曲线。横坐标是迭代次数纵坐标是全局最优适应度值。一个正常的PSO收敛过程应该分三个阶段快速下降期前50代左右缓慢下降期100代以后平缓收敛期200代以后几乎不再变化。如果你的收敛曲线在100代以后还在剧烈波动一般是惯性权重设置太大或者速度上限太高。如果一开始下降很慢可能是惩罚系数太大算法在约束边界附近来回试探迟迟找不到好解。我这里有一个实用技巧画收敛曲线时把种群最优适应度和种群平均适应度同时画出来。如果平均适应度稳步下降且接近最优适应度说明种群多样性保持得不错如果最优适应度不动但平均适应度波动很大说明粒子还在乱飞需要降低速度上限或增大惯性权重下降速度。5.2 调度出力结果验证功率平衡图要逐时段闭合迭代结束后把gbest向量解码成调度计划然后画一个堆叠面积图直观展示每个时段的电源构成。把负荷曲线画在同一个图上逐时段检查P_wind(t) P_hydro(t) P_ps_g(t) P_f(t) - P_ps_p(t) P_load(t)由于我在编码时用功率平衡式计算了火电功率所以这个等式理论上自动满足。但实际运行时可能因为浮点数精度问题两边有微小误差。误差在1e-6量级可以接受如果出现几十千瓦的偏差说明解码函数里某个功率的符号写反了。另外还要单独画一张抽水蓄能的蓄能量S(t)曲线。它的形状应该像波浪一样负荷低谷时上升负荷高峰时下降并且调度周期末回到初始值。如果S(24)和S(1)相差很大说明能量守恒约束没被满足需要回到公式检查。5.3 我在复现中踩过的三个大坑第一个坑是抽蓄综合效率写反。最开始我用的是综合效率eta_total 0.75但公式里写成了S(t1) S(t) eta_total * P_ps_p(t) - (1/eta_total) * P_ps_g(t)这本身没问题。问题出在我把发电效率设置成了1.0等于发电不损耗抽水损耗25%最后算出来的库容曲线没有周期末回到初始值反而越用越多。后来检查才发现综合效率里必须同时考虑抽水和发电损耗不能只在一个环节打折。第二个坑是初始种群全部落在不可行域。96维的决策变量如果不加处理地随机生成很容易生成一组决策变量算出来的火电功率超出限值或者抽蓄库容约束大量违反。这时候适应度函数里惩罚项会非常大粒子需要很多代才能爬进可行域。我的解决办法是初始化解时先让抽蓄功率和水库出力清零然后在范围内随机分配最基础的出力保证初始解接近可行域。这个操作对收敛速度的提升非常明显。第三个坑是Matlab版本兼容性。我在用旧版Matlab重跑代码时发现rand的种子设置函数已经变了还有一些自带的优化工具箱函数名被改过。为了不依赖特定工具箱我把算法核心部分全部用最基础的原生函数实现只用了矩阵运算、循环和plot。这样换个版本也能直接跑。5.4 与论文结果对比的合理误差范围EI复现最终要面对的问题是我复现出的结果和论文里的结果差多少算正常很多EI论文给出的调度结果是基于某一组特定参数和负荷数据。你手中如果没有完全相同的原始数据想1:1复现是做不到的。合理的做法是使用论文的算例拓扑和参数范围自己设一套合理的负荷和风电曲线复现出“相同趋势”的结果。比如论文里说抽蓄能在夜间填谷、白天削峰你的结果也应该能看到这个规律论文里说系统弃风率降低了10个百分点你的模型在相同定性条件下也应该能看出来弃风率下降具体数值可以有一定浮动。如果完全复现不出来先不要怀疑算法先去检查你的模型约束是不是比论文里更严格。我之前复现时发现结果总比论文弃风率高一大截后来发现是论文把风电最大出力曲线平滑处理过而我用的是原始波动数据。把输入数据对齐后结果就靠近了。最后再分享一个小经验粒子群算法是随机算法每次运行结果都会略有不同。所以不要用一次跑出来的结果直接写进论文。我的做法是连续运行30次取最优值、平均值和标准差。如果标准差很小说明算法已经稳定收敛如果标准差很大说明参数还需要调整。这个习惯能帮你避免很多不必要的学术争议。如果你也在复现类似的电力调度模型我希望上面这些内容能帮你少走弯路。先画系统拓扑再列模型公式然后再写代码顺序千万不要反。代码出问题的时候利用Matlab的断点一步步看比瞎猜快得多。祝你的EI复现顺利。
最近在复现一篇EI论文题目是“基于粒子群算法的风电-水电抽水蓄能联合优化调度”。论文不算新但涉及风电、常规水电和抽水蓄能三种电源加上粒子群算法的实现内容非常典型。Matlab代码断断续续写了两周第一阶段算是完整跑通了。这篇博文把整个复现过程的思路、模型、代码结构和容易踩的坑全部记录下来给后面做同类研究的同学一个可以直接参考的路线图。我默认你手里的EI论文是英文或者中文核心问题建模基本围绕“如何让风电与抽水蓄能协同使系统运行成本最小或弃风率最低”展开。因为输入信息里没有给具体的目标函数表达式我采用的是一套在EI文献里非常常见的组合系统运行成本最小化加上弃风惩罚项。后面会解释为什么这样设计以及怎么替换成论文里的其他目标。1. 为什么风电要拉上抽水蓄能一起调度问题本质与EI复现的逻辑起点1.1 风电出力的“反调峰”困境风电出力最大的特点就是波动性和随机性。很多时候夜间负荷低风机却因为风速大而满发白天负荷高风速反而降下来。这种“反调峰”特性让电网调度员非常头疼。负荷高峰期需要风电顶上时它发不出来负荷低谷期又容易发得太多导致电网无法完全消纳只能弃风。在纯火电或纯水电系统里调度问题相对简单因为火电和水电机组的出力都可以主动调节。风电加入后问题变成了一个不确定环境下的决策问题在满足负荷需求、保证电网安全的前提下尽可能多地消纳风电同时让其他电源承担调节压力。这就是风电-水电联合优化调度的第一个逻辑起点风电出力必须先给定或预测出来然后再计算其他机组怎么配合。1.2 抽水蓄能在联合调度中扮演的角色抽水蓄能电站本质上是储能电站。它有两个工作状态抽水时消耗电能把下水库的水抽到上水库相当于给系统增加负荷发电时把上水库的水放下来推动机组发电相当于给系统增加电源。这个“削峰填谷”能力正好弥补风电的波动性。比如夜间风电大发时负荷又低抽水蓄能可以开抽水模式把多余的风电“吞掉”存成水的势能白天负荷高峰或者风电出力下降时再放水发电把能量“吐”出来。这样一来风电的弃风率能降下来系统的峰谷差也能被削平。这里要注意抽水蓄能电站不是常规水电站。常规水电站靠天然来水发电受河流来水量的约束抽水蓄能电站主要靠“水循环”工作抽上去多少水发下来多少电中间还要乘一个综合效率一般在0.75到0.8左右。这个效率因子特别容易在复现代码时写错后面我会专门讲。1.3 EI复现第一步先把论文的“系统拓扑”画出来很多同学拿到论文就直接看公式、抄代码结果跑了半天不知道每个变量代表什么。我的建议是不管EI论文里画没画系统图你都要自己先画一张“电源-负荷”框图。典型的风电-水电抽水蓄能联合系统长这样风电场出力曲线给定可弃风是一个可调减的电源。抽水蓄能电站有抽水工况和发电工况不能同时进行。常规水电站有入库径流有水库库容限制发电功率与发电流量相关。电网负荷每个时段的负荷需求给定是功率平衡约束的右边。可能还有一个火电机组作为基荷或调节电源目的是让系统供电可靠但会增加成本。把这几个部分画出来后你才知道决策变量有哪些约束条件有哪些。这一步看起来简单但能帮你避开“抄错变量含义”这种最蠢的错误。2. 联合优化调度的数学模型目标函数与约束条件的逐条拆解2.1 目标函数运行成本最小还是弃风量最小EI论文里目标函数可以是多种形式最常见的有系统总运行成本最小火电煤耗成本 启停成本 弃风惩罚等弃风量最小系统发电收益最大多目标加权经济性和环保性。我复现时采用的是典型的“运行成本最小化 弃风惩罚”形式。因为单纯弃风量最小没法体现经济性单纯成本最小又可能让系统把风电全弃掉——风电边际成本接近零所以必须加一个弃风惩罚项让算法尽量多用风电。可以写成min F sum_t ( C_f * P_f(t) ) lambda * sum_t ( P_wind_max(t) - P_wind(t) )其中P_f(t)是火电或外购电功率C_f是单位发电成本系数P_wind_max(t)是t时段风电最大可发功率P_wind(t)是t时段实际并网功率lambda是弃风惩罚系数一般取一个相对较大的数比如火电成本的5到10倍。如果你的论文没有火电只有风电、水电和抽蓄那目标函数可能会变成“抽水蓄能网损最小”或“水电出力最大化”但惩罚项思路是一样的。EI复现时先确认目标函数的形式再写代码。2.2 功率平衡与机组出力上下限约束功率平衡是所有电力调度模型必须满足的核心等式约束。每个时段系统总发电功率必须等于负荷功率。P_wind(t) P_hydro(t) P_ps_g(t) - P_ps_p(t) P_f(t) P_load(t)这里P_hydro(t)是常规水电出力P_ps_g(t)是抽蓄发电功率P_ps_p(t)是抽水功率P_load(t)是负荷。注意抽水功率是消耗电能的所以公式左边要减去它。很多初学者把抽水功率从等式右边减去最后功率平衡怎么都调不对。等式约束之外就是各类机组的出力上下限约束0 P_wind(t) P_wind_max(t) P_hydro_min P_hydro(t) P_hydro_max 0 P_ps_g(t) P_ps_g_max * u_g(t) 0 P_ps_p(t) P_ps_p_max * u_p(t) u_g(t) u_p(t) 1最后这个互斥约束特别重要抽水蓄能机组在同一个时段要么抽水要么发电绝对不能同时进行。这个约束如果写漏了算法很容易给出一个“既抽水又发电”的荒谬解从能量角度看等于让水在上水库和下水库之间白做功效率还不为0结果会非常失真。2.3 抽水蓄能电站的核心约束抽发循环的能量守恒抽水蓄能本质上是在“搬运”能量所以必须满足一个跨时段的水量/能量守恒约束。如果不加这个约束电站可以凭空发很多电系统功率平衡也能满足但那是不合理的。最常见的模型是用上水库的蓄能状态来描述S(t1) S(t) eta_p * P_ps_p(t) * delta_t - P_ps_g(t) * delta_t / eta_g这里S(t)是上水库t时段的蓄能单位MWh或m3水当量eta_p是抽水效率通常0.85左右eta_g是发电效率通常0.9左右delta_t是每个时段的时长单位小时24时段模型里通常为1。这个公式的意思是抽水时蓄能增加但因为有损耗所以乘以效率发电时蓄能减少因为同样有损耗所以除以效率。如果觉得两个效率相乘麻烦也可以用一个综合效率eta_total eta_p * eta_g公式变成S(t1) S(t) eta_total * P_ps_p(t) * delta_t - P_ps_g(t) * delta_t但它只能用于粗略估算严谨一点还是分开写。另外还需要上下水库容量约束S_min S(t) S_max S(1) S(24) S_init这个“调度周期始末蓄能相等”的约束非常关键否则算法会把上水库的水一次性放完得到一个不可持续的解。我在最初复现时漏了这个结果收敛曲线很漂亮但调度结果完全不合理——抽蓄电站最后库容见底这个解放到实际系统中根本没法执行。2.4 常规水电及水库水量平衡约束常规水电和抽水蓄能不一样它有一个天然来水过程。典型模型里水电出力与发电流量和水库水头相关。简单做法是用“水库蓄水量”作为状态变量V(t1) V(t) Q_in(t) - Q_hydro(t) - Q_spill(t)其中V(t)是水库库容Q_in(t)是天然入库流量Q_hydro(t)是发电流量Q_spill(t)是弃水流量。水电出力可以近似为P_hydro(t) k * Q_hydro(t) * H_net(t)如果假设水头恒定那么出力就和发电流量近似成正比变成一个线性关系。EI复现中很多论文都做了这个简化你不需要一上来就搞水头非线性模型先复现线性版本再考虑扩展。库容约束和初始库容约束同样要写V_min V(t) V_max V(1) V_init V(24) V_end如果论文里没给水库初始库容你可以设置一个合理值比如库容上限的50%并在代码里作为参数保留方便后面测试不同场景。3. 粒子群算法如何落到调度场景编码方式、适应度函数与参数整定3.1 为什么选粒子群而不是遗传算法或线性规划这个问题我复现前也纠结过。联合优化调度其实可以建模成混合整数线性规划用cplex或gurobi这类求解器也能解。但EI论文明确用了粒子群算法复现的目的就是要还原论文的方法。而且粒子群在非线性、非凸目标函数上有天然优势不需要求导不依赖初始点实现代码也短适合快速验证模型。和遗传算法相比粒子群收敛速度更快、参数更少但容易早熟。不过对于24时段、决策变量不算太多的调度问题粒子群的早熟问题可以通过增加种群规模、调节惯性权重来缓解。所以它成了很多电力调度论文的首选智能算法。3.2 粒子编码把24小时的调度序列排成一个向量粒子群算法里每个粒子代表一个候选解。在调度问题中候选解就是一整天的调度计划。怎么把调度计划编码成一个“位置向量”是写代码前必须想清楚的事情。我采用的编码方式非常直接位置向量 [所有时段的抽蓄发电功率P_ps_g(1..24), 所有时段的抽蓄抽水功率P_ps_p(1..24), 所有时段的常规水电出力P_hydro(1..24), 所有时段的并网风电功率P_wind(1..24)]假设调度周期有T个时段那么每个粒子的维度就是D 4 * T。如果T24D96。你可能要问为什么不把火电功率也作为决策变量因为火电功率可以由功率平衡约束直接算出来P_f(t) P_load(t) - P_wind(t) - P_hydro(t) - P_ps_g(t) P_ps_p(t)这是处理等式约束常用的技巧把等式约束变成一个“计算式”只要其余变量定了火电功率就唯一确定等式约束自动满足。这样粒子编码维度少了一个变量算法搜索空间也更小。位置向量的每个元素都有物理边界。比如风电并网功率不能超过该时段最大可发功率抽蓄功率不能超过额定功率。这些边界在初始化位置时就要设好后面迭代更新位置时也要做边界处理。3.3 适应度函数与惩罚函数设计粒子群算法需要用一个适应度值来衡量每个粒子的好坏。在最小化问题里适应度值越小越好。我直接把目标函数值作为基础再加上约束违反惩罚。为什么需要惩罚因为粒子更新过程中位置向量是随机变化的很可能不满足所有约束。比如抽水蓄能水库的容量约束可能被打破常规水电的库容可能超出上限。如果这些不满足约束的粒子直接参与比较搜索会被带偏。常用方法是“外点罚函数法”fitness F sigma1 * sum(约束违反量^2) sigma2 * ...具体到我的代码里约束包括抽蓄上水库蓄能上下限约束常规水电库容上下限约束抽蓄始末蓄能相等约束水电站始末库容相等约束抽水与发电状态互斥约束火电功率上下限约束用功率平衡等式算出的火电可能超限。每个约束违反的绝对值会乘一个惩罚系数然后平方累加。惩罚系数不能太小否则约束形同虚设也不能太大否则会主导适应度函数让算法只顾满足约束而忽略经济性。我的经验是先跑一次不惩罚的版本看目标函数量级然后把惩罚系数设为目标函数量级的10到100倍再逐步调整。3.4 参数整定的几个经验值粒子群算法的标准参数有四个种群规模N、最大迭代次数iter_max、惯性权重w、加速常数c1和c2。我复现时用的一组默认参数如下参数取值说明种群规模N100对于96维问题100个粒子不算大但已经能收敛最大迭代次数300后期看收敛曲线调整如果没平就加惯性权重w0.9到0.4线性递减前期全局搜索后期局部搜索学习因子c12.0个体认知权重学习因子c22.0社会经验权重速度上限Vmax每个变量范围的20%防止粒子飞得太远惯性权重递减是非常经典的策略公式是w w_max - (w_max - w_min) * iter / iter_max前期的w比较大粒子飞行速度快能覆盖整个搜索空间后期w变小粒子在局部精细搜索。这个策略对调度模型特别管用因为初始阶段很多粒子可能落在不可行域需要大范围探索才能找到可行区域。4. Matlab代码架构详解主程序、约束处理与循环迭代的骨架4.1 主程序流程从参数初始化到结果出图写Matlab代码时我习惯把整个复现过程分成逻辑清晰的三层主脚本、函数文件和绘图脚本。主脚本负责参数设置、初始化、调用优化循环、输出结果。函数文件包括目标函数、约束校验、粒子更新等。这样代码好读也好改。主程序的核心流程可以描述为设置基础数据负荷曲线、风电最大出力曲线、水电参数、抽蓄参数、算法参数。初始化粒子群随机生成N个位置向量和N个速度向量。计算初始适应度对每个粒子调用适应度函数更新个体最优和全局最优。迭代循环更新速度、更新位置、边界处理、计算适应度、更新pbest和gbest。输出最优解把gbest位置解码成调度计划绘制功率平衡图、库容曲线、收敛曲线。我给的代码骨架是经过简化的但核心逻辑和完整版一致。实际复现时你可以在每个循环中间加一个调试打印观察适应度值的变化趋势。主脚本开头通常是数据区我用结构体把系统参数存起来避免函数参数列表太长%% 参数设置 T 24; % 时段数 N 100; % 粒子群规模 iter_max 300; % 最大迭代次数 c1 2.0; c2 2.0; % 学习因子 w_max 0.9; w_min 0.4; % 负荷与风电预测数据示例 load_curve [ /* 24个负荷值 */ ]; wind_max [ /* 24个风电最大出力值 */ ]; % 抽水蓄能参数 P_ps_g_max 80; % 发电最大功率 P_ps_p_max 80; % 抽水最大功率 S_max 500; % 上水库最大蓄能 S_min 50; % 上水库最小蓄能 S_init 250; % 初始蓄能 eta_p 0.85; % 抽水效率 eta_g 0.90; % 发电效率 % 常规水电参数 P_hydro_max 120; P_hydro_min 20; V_max 1000; V_min 200; V_init 600; Q_in ones(1, T) * 20; % 入库流量 k_h 0.8; % 水电出力系数这些参数只代表一个测试算例实际复现时请根据论文中的测试系统填写。4.2 核心代码片段PSO速度与位置更新PSO更新是算法的核心。标准的速度更新公式v w * v c1 * rand(size(v)) .* (pbest - x) c2 * rand(size(v)) .* (gbest - x);位置更新x x v;速度限制v max(v, -Vmax); v min(v, Vmax);位置边界处理可以用“边界吸收”x(x lb) lb(x lb); x(x ub) ub(x ub);其中lb和ub是所有决策变量的上下界向量。以抽蓄发电功率为例下界是0上界是P_ps_g_max抽水功率同理。风电并网功率上界是wind_max中的对应值。完整的更新循环大概是这样的for iter 1:iter_max w w_max - (w_max - w_min) * (iter / iter_max); for i 1:N v(i,:) w * v(i,:) c1*rand(1,D).*(pbest(i,:) - x(i,:)) ... c2*rand(1,D).*(gbest - x(i,:)); v(i,:) max(v(i,:), -Vmax); v(i,:) min(v(i,:), Vmax); x(i,:) x(i,:) v(i,:); x(i,:) max(x(i,:), lb); x(i,:) min(x(i,:), ub); fit(i) fitness_func(x(i,:), sys_param); if fit(i) fit_pbest(i) pbest(i,:) x(i,:); fit_pbest(i) fit(i); end if fit(i) fit_gbest gbest x(i,:); fit_gbest fit(i); end end best_history(iter) fit_gbest; end这段代码看起来简单但有几个细节需要注意。rand(1,D)每次在速度更新公式里应该生成不同的随机数所以必须在表达式里直接调用不能提前保存一个固定向量否则算法会失去随机性。另一个细节是粒子的历史最优。在电力调度问题中即使当前迭代的适应度值不是很好历史最优里可能已经存了一个可行且经济的解。所以pbest的更新必须严格用“小于”而不是“小于等于”否则可能出现历史最优被一个“更差但数值上等于”的解覆盖的问题。4.3 约束处理的实现细节边界吸收与状态互斥边界吸收是处理上下界约束最简单的方法。位置更新后如果某个决策变量超过上界就直接拉回上界低于下界就拉回下界。这个方法的好处是实现简单不容易产生意料之外的解。但缺点是会让靠近边界的粒子失去速度信息可能降低搜索效率。另一种方法是“随机重置到边界附近”实现稍微复杂一点但效果更好。抽水与发电状态互斥约束不能靠边界处理因为它是两个变量之间的耦合关系。我写了一个专门的处理函数如果一个粒子里同一时段的P_ps_g(t)和P_ps_p(t)都大于0就按一个随机偏向保留其中一个另一个置0。比如生成一个随机数如果大于0.5就保留发电、清掉抽水否则保留抽水、清掉发电。这样能保证互斥约束成立。不过要注意这种“修复策略”在每一代都做等于人为改变了粒子的位置。如果你的惩罚函数设计得好也可以不强制修复而是让算法通过惩罚自动学会“不要同时抽水和发电”。我个人的经验是先强制修复保证每个粒子都满足互斥约束这样跑出来的结果物理上一定合理对排查bug也有帮助。4.4 代码组织建议脚本、函数与测试用例完整代码里我建议至少包含这几个文件main_PSO_optimization.m主脚本fitness_func.m适应度函数目标 惩罚decode_solution.m把位置向量解码成调度计划constraint_fix.m互斥约束修复、边界处理plot_results.m出图脚本。函数模块化的好处是出了问题可以单点调试。我最开始把所有逻辑都写在一个脚本里出了bug定位很慢。后来拆成函数通过打印中间变量很快就发现是抽蓄能量守恒公式里的效率除反了。5. 复现结果怎么判断对不对收敛曲线、调度出力与几个容易翻车的细节5.1 收敛曲线的判读标准不管用什么优化算法第一步都要画收敛曲线。横坐标是迭代次数纵坐标是全局最优适应度值。一个正常的PSO收敛过程应该分三个阶段快速下降期前50代左右缓慢下降期100代以后平缓收敛期200代以后几乎不再变化。如果你的收敛曲线在100代以后还在剧烈波动一般是惯性权重设置太大或者速度上限太高。如果一开始下降很慢可能是惩罚系数太大算法在约束边界附近来回试探迟迟找不到好解。我这里有一个实用技巧画收敛曲线时把种群最优适应度和种群平均适应度同时画出来。如果平均适应度稳步下降且接近最优适应度说明种群多样性保持得不错如果最优适应度不动但平均适应度波动很大说明粒子还在乱飞需要降低速度上限或增大惯性权重下降速度。5.2 调度出力结果验证功率平衡图要逐时段闭合迭代结束后把gbest向量解码成调度计划然后画一个堆叠面积图直观展示每个时段的电源构成。把负荷曲线画在同一个图上逐时段检查P_wind(t) P_hydro(t) P_ps_g(t) P_f(t) - P_ps_p(t) P_load(t)由于我在编码时用功率平衡式计算了火电功率所以这个等式理论上自动满足。但实际运行时可能因为浮点数精度问题两边有微小误差。误差在1e-6量级可以接受如果出现几十千瓦的偏差说明解码函数里某个功率的符号写反了。另外还要单独画一张抽水蓄能的蓄能量S(t)曲线。它的形状应该像波浪一样负荷低谷时上升负荷高峰时下降并且调度周期末回到初始值。如果S(24)和S(1)相差很大说明能量守恒约束没被满足需要回到公式检查。5.3 我在复现中踩过的三个大坑第一个坑是抽蓄综合效率写反。最开始我用的是综合效率eta_total 0.75但公式里写成了S(t1) S(t) eta_total * P_ps_p(t) - (1/eta_total) * P_ps_g(t)这本身没问题。问题出在我把发电效率设置成了1.0等于发电不损耗抽水损耗25%最后算出来的库容曲线没有周期末回到初始值反而越用越多。后来检查才发现综合效率里必须同时考虑抽水和发电损耗不能只在一个环节打折。第二个坑是初始种群全部落在不可行域。96维的决策变量如果不加处理地随机生成很容易生成一组决策变量算出来的火电功率超出限值或者抽蓄库容约束大量违反。这时候适应度函数里惩罚项会非常大粒子需要很多代才能爬进可行域。我的解决办法是初始化解时先让抽蓄功率和水库出力清零然后在范围内随机分配最基础的出力保证初始解接近可行域。这个操作对收敛速度的提升非常明显。第三个坑是Matlab版本兼容性。我在用旧版Matlab重跑代码时发现rand的种子设置函数已经变了还有一些自带的优化工具箱函数名被改过。为了不依赖特定工具箱我把算法核心部分全部用最基础的原生函数实现只用了矩阵运算、循环和plot。这样换个版本也能直接跑。5.4 与论文结果对比的合理误差范围EI复现最终要面对的问题是我复现出的结果和论文里的结果差多少算正常很多EI论文给出的调度结果是基于某一组特定参数和负荷数据。你手中如果没有完全相同的原始数据想1:1复现是做不到的。合理的做法是使用论文的算例拓扑和参数范围自己设一套合理的负荷和风电曲线复现出“相同趋势”的结果。比如论文里说抽蓄能在夜间填谷、白天削峰你的结果也应该能看到这个规律论文里说系统弃风率降低了10个百分点你的模型在相同定性条件下也应该能看出来弃风率下降具体数值可以有一定浮动。如果完全复现不出来先不要怀疑算法先去检查你的模型约束是不是比论文里更严格。我之前复现时发现结果总比论文弃风率高一大截后来发现是论文把风电最大出力曲线平滑处理过而我用的是原始波动数据。把输入数据对齐后结果就靠近了。最后再分享一个小经验粒子群算法是随机算法每次运行结果都会略有不同。所以不要用一次跑出来的结果直接写进论文。我的做法是连续运行30次取最优值、平均值和标准差。如果标准差很小说明算法已经稳定收敛如果标准差很大说明参数还需要调整。这个习惯能帮你避免很多不必要的学术争议。如果你也在复现类似的电力调度模型我希望上面这些内容能帮你少走弯路。先画系统拓扑再列模型公式然后再写代码顺序千万不要反。代码出问题的时候利用Matlab的断点一步步看比瞎猜快得多。祝你的EI复现顺利。
