做微网优化调度这块也有几年了Matlab里PSO粒子群算法解决过不少问题但“考虑需求响应的微网优化调度模型”这个组合我每次看都还是觉得挺有代表性的。它表面上是跑一个粒子群算法实际上把微网的源、荷、储、网都给串起来了还要把需求响应这种柔性的负荷侧资源加进去。今天就把这个模型的全部细节拆开讲一遍包括目标函数怎么写、约束怎么加、粒子怎么编码、代码怎么搭以及我实际跑仿真时踩过的坑。这篇内容主要适合正在做微网经济调度、研究需求响应策略、或者想用粒子群算法做能源优化方向的同学参考。1. 先搞清楚这个模型在干什么微网调度和需求响应到底优化了什么1.1 微网的源、荷、储、网结构微网本质上是一个小型的电力系统麻雀虽小五脏俱全。常见的组成包括分布式电源光伏PV、风电WT、微型燃气轮机MT、燃料电池FC、储能系统BESS、负荷以及跟大电网之间的联络线。优化调度的目标就是未来24小时通常以1小时为调度间隔内决定每个时刻各台机组出力多少、储能充放多少、跟电网交换功率多少使得总运行成本最低同时满足安全运行约束。这里边有一个容易忽略的点微网不是越大越好也不是发电越多越好。因为光伏、风电这类可再生电源虽然有发电成本低的好处但有间歇性和波动性不可控。微型燃气轮机虽说可调可控但有燃料成本。储能则是一个“时间搬运工”低电价时充、高电价时放。而需求响应又加了一层部分负荷是可以调节、平移、削减的这给了调度更多腾挪空间。1.2 需求响应是怎么把“负荷”变成可调变量的需求响应简单说就是让用户侧根据电价信号或激励政策主动调整用电行为。在调度模型里通常把负荷分成刚性负荷和柔性负荷两类刚性负荷比如医院设备、生产线核心工艺必须按时供电。柔性负荷比如空调、洗衣机、热水器、充电桩可以在一定时间范围内调整。柔性负荷又可以进一步细化成可削减负荷、可转移负荷和可平移负荷。可转移负荷的特点是“总量守恒、时段可变”比如工业车间的某道工序可以从高峰时段挪到低谷时段但总的用电量不能变。可削减负荷则是直接减少某个时段的用电用户会获得补偿但总用电量下降了。这个特性加到调度模型里边最核心的变化就是功率平衡等式里多了一批“负荷侧决策变量”。调度结果不再是单纯告诉机组“你什么时候发多少电”还要告诉负荷“你什么时候用多少电、用电的时段怎么重新分配”。1.3 目标函数不是只有一个成本要拆开算很多新手写调度模型的时候以为目标函数就是一段简单的总成本写成 min C 购电成本 燃料成本。一旦需求响应加入目标函数至少要拆成下面几块第一块是向大电网购电成本价格通常采用分时电价峰平谷三个时段电价不同。这个直接跟联络线功率和当前时刻电价相乘即可累加。第二块是微型燃气轮机的燃料成本。常见方法是把燃料成本拟合成一个二次函数或者是分段线性函数表达式一般是 C_fuel a * P^2 b * P c 这种形式。第三块是储能和分布式电源的运行维护成本。这部分一般按功率线性折算即 C_om k * Pk是每个单位发电量对应的维护费用。第四块就是需求响应成本。如果做了可削减负荷就要给用户补偿费用一般按削减电量乘一个补偿系数。如果做了可转移负荷模型里通常不直接算钱的转移但要让转移前后的总用电量保持不变。整体目标函数写下来就是这几个部分的分时累加然后乘上时间段长度得到全天24小时的总成本。需要注意电价单位、功率单位、时间单位要统一不然算出来结果差一个数量级还找不到原因。1.4 约束条件才是调度的“硬骨头”目标函数好说约束条件才是真正决定了模型能不能跑、结果合不合理的核心。我见过太多人只把目标函数建得很漂亮约束条件写得太宽松结果粒子群一搜调度方案完全违反物理规律。最基本的约束包括功率平衡约束任一时刻所有电源出力加上电网购电、储能放电等于负荷需求加上储能充电、负荷削减或转移后的净负荷。这是调度模型最核心的等式约束。机组出力上下限光伏、风电、燃气轮机的出力不能超过额定范围。爬坡约束燃气轮机在相邻时段功率变化量有限制不能突然从零跳到满发。储能SOC约束荷电状态要保持在安全区间比如0.1到0.9之间每个时段的充放电功率也有上限充放电效率要体现出来一天结束时的SOC最好跟初始时刻一致保证日循环的可持续性。联络线功率约束微网跟大电网交换功率不能超过线路容量。需求响应约束可削减负荷量有上限可转移负荷必须满足转移前后总电量守恒负荷调节量要在用户可接受的范围内。这些约束放到粒子群算法里处理起来不太轻松尤其是等式约束。我的习惯做法是用惩罚函数把约束不满足的量加到目标函数上让不满足约束的粒子在迭代中被淘汰掉。当然也有人偏好通过编码方式直接满足部分约束后面在代码部分会细说。2. 为什么用粒子群算法选型思路和算法原理2.1 调度问题为什么可以交给智能算法严格来说微网优化调度本质上是一个混合整数规划问题因为储能充放电状态、机组启停状态都是0-1整数变量数学上可以用CPLEX、Gurobi这类商业求解器去解。那为什么还有很多人用粒子群算法我的理解有三点原因第一模型一旦加入需求响应、非线性爬坡约束、非凸运行区域问题会变成非凸非线性混合整数规划商业求解器未必能在合理时间内给出全局最优解而智能算法可以在可接受的计算时间内给出一个质量不错的近似最优解。第二粒子群算法实现简单没有梯度信息要求对目标函数是否可导、是否连续不那么敏感。只要能把解编码成粒子位置并且能算出适应度值就能用。对做科研、做课题验证的同学来说这大大降低了门槛。第三粒子群算法的参数含义直观收敛过程可视往论文里放收敛曲线也方便。虽然理论上不保证全局最优但通过多策略改进惯性权重递减、变异操作等实际效果一般都够用。2.2 粒子群算法的核心公式粒子群算法模拟鸟群觅食行为。每个粒子代表一个候选解粒子在搜索空间中飞行通过不断修正自己的速度来逼近全局最优位置。核心就是速度和位置更新公式v_{i,d}^{(k1)} w * v_{i,d}^{(k)} c1 * r1 * (pbest_{i,d} - x_{i,d}^{(k)}) c2 * r2 * (gbest_d - x_{i,d}^{(k)})x_{i,d}^{(k1)} x_{i,d}^{(k)} v_{i,d}^{(k1)}其中w是惯性权重控制前一时刻速度对当前速度的影响w大有利于全局搜索w小有利于局部精细搜索。c1和c2是学习因子c1偏向个体认知c2偏向社会群体。r1和r2是[0,1]之间的随机数。这里有一个工程上很实用的改进惯性权重线性递减。比如初始w从0.9慢慢降到0.4让算法前期在整个空间里找大方向后期在最优解附近精细搜索。我在代码里基本都这么做收敛速度和精度都比固定w好不少。2.3 调度方案是怎么编码成粒子的粒子群算法处理微网调度最关键的一步就是定义粒子位置对应什么决策变量。我常用的编码方式是把每个时段的可控变量拼接成一个长向量。假设调度周期H24小时决策变量包括燃气轮机出力P_mt(1:24)、储能充放电功率P_bat(1:24)、联络线功率P_grid(1:24)、可削减负荷P_cl(1:24)、可转移负荷P_tr(1:24)。如果都用实数编码粒子位置维数就是 5×24120 维。每一维都有明确的物理意义。储能充放电这里有个处理技巧可以用一个带符号的变量表示净充放电功率正值表示放电负值表示充电但需要在约束里限制充电功率和放电功率分别不超过上限同时不允许同时充放。转速问题在于这种自然编码方式的可行域是“条带状”的。粒子更新之后的解很可能违反约束所以后面必须紧跟约束修正步骤或罚函数评价。3. Matlab代码实现从数据到结果的完整流程3.1 场景设置和数据准备我实际搭建模型的时候会先建一个场景结构体把仿真日期、调度时段、负荷曲线、风光出力曲线、价格参数全部放进去。这里给出一个可以参考的示例设置。一天分为24个时段光伏和风电数据采用典型日的归一化曲线然后乘以额定容量。负荷曲线用典型的居民/商业混合负荷形状峰时段出现在上午和晚间谷时段出现在凌晨。分时电价设置为峰时1.2元/kWh、平时0.8元/kWh、谷时0.4元/kWh。下面是一个简单的数据准备代码骨架% 基本参数设置 H 24; % 调度时段数 P_load_base [ ... ]; % 基础负荷曲线1x24向量 P_pv [ ... ]; % 光伏预测出力 P_wt [ ... ]; % 风电预测出力 % 分时电价 price zeros(1, H); price(10:12) 1.2; % 峰 price(19:21) 1.2; % 峰 price(1:7) 0.4; % 谷 price(find(price 0)) 0.8; % 平时 % 燃气轮机参数 mt struct(Pmin, 20, Pmax, 200, a, 0.001, b, 0.2, c, 10, ... ramp, 50, om, 0.05); % 储能参数 bess struct(cap, 500, soc_init, 0.5, soc_min, 0.1, ... soc_max, 0.9, eta_ch, 0.95, eta_dis, 0.95, ... Pmax_ch, 100, Pmax_dis, 100, om, 0.02); % 需求响应参数 dr struct(rate, 0.5, Pcl_max, 30, Pshift_max, 40, ... compensation, 0.6);这里需要注意所有数组都应该是1×24的行向量否则后面矩阵运算容易出错。我自己踩过这个坑在构建适应度函数时维度不一致导致矩阵相加失败浪费了不少时间。3.2 目标函数的代码化目标函数是整个程序的核心。为了便于粒子群调用我会把它封装成一个独立的函数文件输入是粒子位置向量x输出是总成本和约束违反量。用上面提到的编码方式粒子x的前24维是燃气轮机出力25到48维是储能功率带符号49到72维是联络线功率73到96维是可削减负荷97到120维是可转移负荷。下面是目标函数的核心计算逻辑function [cost, viol_sum] objective_func(x, params) H params.H; % 解码决策变量 P_mt x(1:H); P_bat x(H1:2*H); % 正放电负充电 P_grid x(2*H1:3*H); P_cl x(3*H1:4*H); P_tr x(4*H1:5*H); % 更新负荷基础负荷 - 可削减负荷 可转移负荷影响 P_load_net params.P_load - P_cl P_tr; % 成本计算 cost_fuel sum(params.mt.a .* P_mt.^2 params.mt.b .* P_mt params.mt.c); cost_grid sum(P_grid .* params.price); cost_om sum(params.mt.om .* P_mt) sum(params.bess.om .* abs(P_bat)); cost_dr sum(P_cl) * params.dr.compensation; cost cost_fuel cost_grid cost_om cost_dr; % 约束违反量先留空在下一步计算 viol_sum 0; end这里解释一下为什么需求响应成本里只计算可削减负荷的补偿。可转移负荷虽然改变了用电时段但总用电量没有变化在“用电量守恒”的约束下它不会改变总用电成本只影响时序上的成本分布所以目标函数里不需要单独给它计费。这个思路很多初看的同学会绕不过弯仔细想一下就明白了。3.3 约束处理罚函数与边界控制粒子群本身不知道“哪些位置是合法的”所以约束处理必须由我们自己实现。这里我采用组合策略一部分简单约束直接通过边界截断来处理另一部分复杂约束通过罚函数来惩罚。边界截断很好理解。比如燃气轮机出力如果超出最大/最小出力直接把它的值拉到边界上。再比如储能功率绝对值超过额定功率就直接截到额定值。这样能保证粒子位置的变量值都在“物理可行范围”内。但像功率平衡这种等式约束没有办法靠截断解决只能计算偏差并加入罚函数% 功率平衡等式约束 bal P_grid P_pv P_wt P_mt P_bat - P_load_net; % 储能SOC计算 SOC(1) params.bess.soc_init; for t 1:H if P_bat(t) 0 SOC(t1) SOC(t) - P_bat(t) / params.bess.eta_dis / params.bess.cap; else SOC(t1) SOC(t) - P_bat(t) * params.bess.eta_ch / params.bess.cap; end end % SOC必须在区间内 soc_min_viol sum(max(0, params.bess.soc_min - SOC(2:end))); soc_max_viol sum(max(0, SOC(2:end) - params.bess.soc_max)); % 爬坡约束 ramp_viol sum(max(0, abs(diff(P_mt)) - params.mt.ramp)); % 可转移负荷总量守恒 shift_viol abs(sum(P_tr)); % 其它约束按需加入最后加权求和 viol_sum sum(abs(bal)) soc_min_viol soc_max_viol ramp_viol ... shift_viol;罚函数的核心思想是把“不满足约束的程度”加到目标函数上让违反约束的粒子适应度变差。所以最终的适应度函数是fit cost penalty_factor * viol_sum;惩罚系数很有讲究调太大会导致算法过早收敛到不可行域边缘调太小则最优解可能绕过约束。我一般先用较小的值让算法能搜开再用较大的值做二次校验。3.4 PSO主循环典型写法粒子群主循环的写法有很多种我建议把它写成“参数结构体 主循环”的形式方便调参。首先初始化粒子群nPop 100; % 种群规模 maxIt 200; % 最大迭代次数 dim 5 * H; % 决策变量维数 X initialization(nPop, dim); % 初始化粒子位置 V zeros(nPop, dim); % 初始速度为零初始化粒子位置时要注意不能全部用随机均匀分布。比如联络线功率的初值我习惯根据负荷和可再生能源出力的差来设置这样初始种群中有一部分粒子已经离可行域比较近了能显著加快收敛速度。然后进入迭代主循环for it 1:maxIt w 0.9 - (0.9 - 0.4) * it / maxIt; % 惯性权重线性递减 for i 1:nPop % 计算适应度 [cost_i, viol_i] objective_func(X(i, :), params); fit_i cost_i penalty * viol_i; % 更新个体最优 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 % 更新粒子速度和位置 for i 1:nPop V(i, :) w * V(i, :) ... c1 * rand(1, dim) .* (pbest(i, :) - X(i, :)) ... c2 * rand(1, dim) .* (gbest - X(i, :)); X(i, :) X(i, :) V(i, :); % 边界约束处理 X(i, :) boundary_check(X(i, :), params); end % 记录每代最优适应度 history(it) fit_gbest; end这里还有一个小细节速度V如果不加限制粒子可能飞出去很远导致算法后期震荡。我习惯同时限制速度上下限是位置范围的10%~20%左右。在储能功率、机组出力等不同量纲的变量之间速度限幅要分别设置不能用一个统一的数值。3.5 结果可视化收敛曲线和调度方案图程序跑完最后一步是把结果画出来。一般至少需要两类图一类是适应度收敛曲线用来展示算法的迭代优化过程另一类是功率平衡调度方案图把24小时的电源出力、储能充放电、负荷需求画成堆叠面积图或条形图。功率平衡图是判断结果是否合理的直接手段。我常用的绘图逻辑是横轴是24小时纵轴是功率。正方向画光伏、风电、燃气轮机、储能放电、购电负方向画负荷、储能充电、售电如果有。这样堆叠起来每个时刻正方向总面积减去负方向总面积应该近似为零说明功率平衡。Matlab里画堆叠图用area函数比较方便figure; area(t, [P_pv; P_wt; P_mt; P_bat_dis; P_grid], LineWidth, 0.5); hold on; plot(t, P_load_net, k-, LineWidth, 1.5); legend(光伏, 风电, 燃气轮机, 储能放电, 电网购电, 净负荷); xlabel(时间/h); ylabel(功率/kW);如果堆叠图里某个时刻的气轮机出力负值、储能同时充放电这类不符合逻辑的情况出现基本可以断定约束没写对或者解码有bug这是快速自查的好办法。4. 仿真结果怎么看效益对比和参数调优4.1 需求响应前后对比加入需求响应前后的差别理论上应该是峰值负荷降低、低谷负荷抬高、总运行成本下降。为什么会有这样的效果因为在分时电价下负荷峰值时段往往也是电价最高的时段。如果此时微网自身出力不足需要高价从电网买电成本非常高。需求响应把一部分负荷转移到电价较低的时段相当于把购电成本的峰值抹平了一部分。我某次仿真的典型结果对比如下数据经过了归一化处理不代表真实系统但反映趋势指标无需求响应有需求响应变化峰值负荷/kW1000872-12.8%谷值负荷/kW45051815.1%峰谷差/kW550354-35.6%总运行成本/元1286011540-10.2%可以看到峰谷差明显被削峰填谷同时成本下降了约10%。这就是需求响应在微网调度里的核心价值也是论文里最常放的对比分析。4.2 粒子群参数怎么试出来粒子群算法最让人头疼的就是参数调优。不同问题的参数响应面不一样网上传的“标准参数”只能作为起点。我从大量实验里总结出来一套经验值种群规模nPop对于120维的决策变量100个粒子是最低配置150个更稳妥。太少容易早熟太多则单次迭代计算量过大。最大迭代次数maxIt一般200到500代足够。判断标准是收敛曲线后期是否平坦。惯性权重w线性递减从0.9到0.4是最常规的设置基本不用改。个体学习因子c1取1.5到2.0越大越强调个体经验。社会学习因子c2取1.5到2.0越大越强调群体共享。速度限幅系数取位置边界的10%~20%。调参的时候我建议不要一次调好几个参数很难判断是哪个起了作用。先固定c1c22和w递减把种群规模和迭代次数定下来再微调学习因子。每次调完看三样东西最终适应度下降了多少、收敛曲线是否平滑、全局最优解对应约束是否完全满足。4.3 初始种群和随机种子对结果的影响实际跑粒子群的人会发现同一个参数、同一个模型每次运行结果可能都不一样。这是因为PSO本质是随机搜索算法初始种群和随机数序列不同都会带来差异。如果做研究写论文最忌讳只跑一次就报结果。可靠的做法是对同一组参数运行10次甚至20次记录每次的成本和运行时间。计算平均值、标准差、最优值和最差值。用最优值展示调度方案用收敛曲线的平均趋势展示算法性能。标准差如果特别大说明算法稳定性差需要增大种群或者增加迭代次数。标准差很小但平均值不好说明算法整体陷入了局部最优要调整参数或者引入变异机制。在Matlab里可以通过设置随机种子让实验时可复现rng(42); % 固定随机种子这个细节在复现别人代码时尤其重要。很多网上的代码没有固定随机种子你跑出来的结果跟论文对不上不一定是代码错了而是随机数序列不一样。5. 常见问题和排错经验5.1 罚函数惩罚系数怎么定惩罚系数太小粒子会停留在不可行区域因为罚函数“力度不够”约束违反带来的额外成本比省下的成本还要少。惩罚系数太大算法会非常保守所有粒子只敢在可行域边缘活动不利于搜索到真正的最优解。我实际操作的时候用的是动态罚函数策略迭代前期惩罚系数稍微小一点比如100让粒子敢于探索更广的空间迭代后期逐渐加大到1000甚至更高迫使粒子最终回到可行域。这样做的好处是前期不会漏掉潜在优解后期又能保证解的可行性。另外不同约束的违反量量纲不一样。功率平衡偏差单位是kWSOC越界单位是%爬坡偏差单位也是kW。如果直接全部乘同一个系数会出现某个约束主导罚函数的情况。我建议每个约束单独设置系数先归一化再求和或者根据量级手动放大。5.2 SOC曲线出现锯齿或越界SOC曲线锯齿状波动通常原因是储能充放电功率频繁变号。也就是这个时段在充电下个时段放电再下个时段又充电。虽然每个时刻的功率上限满足约束但频繁切换充放电状态在实际系统中会严重降低储能寿命调度结果也不够真实。这类问题可以适当改动惩罚函数增加一个对相邻时段充放电状态切换的惩罚项让算法倾向于连续充或连续放。也可以设置一个充放电最小持续时间约束但用粒子群实现起来比较麻烦。SOC越界的原因一般不是约束没写而是SOC约束被放到了罚函数里但惩罚系数不够导致算法认为轻微越界比换一个更优方案更划算。解决方法是把SOC越界项的惩罚系数单独调大。5.3 算法陷入局部最优粒子群最容易出现的问题是早熟收敛整群粒子抱团在一个局部最优解附近全局搜索能力完全丧失。我常用的排查方法是打印出每一代全局最优的多样性指标比如全体粒子位置的标准差。如果在迭代初期标准差就趋近于零说明种群多样性已经崩溃无论迭代多少次结果都不会再变好。解决办法有几个方向。第一增大惯性权重的初值提高前期探索能力第二引入变异操作按一定概率随机重置部分粒子的部分维度第三引入扰动策略当全局最优连续若干代无改善时给当前粒子群加一个随机偏移。这些都算是粒子群的经典改进方向实现起来也简单效果立竿见影。5.4 仿真太慢怎么办如果调度时段增加比如15分钟一个时段24小时有96个时段决策变量维数会从120维飙升到480维粒子群运行时间会快速上涨。最直接的加速手段是向量化计算。早期版本我用for循环逐粒子逐一计算目标函数跑200代、100个粒子就要两三分分钟。后来改成把目标函数里的部分计算向量化、利用矩阵一次性处理多个粒子或者在主循环里用parfor并行计算粒子群的适应度速度能提升很多倍。Matlab并行计算池可以通过以下方式启动parpool;然后把粒子群主循环里的for i 1:nPop改成parfor i 1:nPop。但要注意目标函数里如果有全局变量或者随机数流并行环境下需要处理好随机种子的独立性问题否则每个worker产生的随机数可能关联影响解的质量。除了并行还可以降低不必要的约束计算频率。比如SOC计算必须逐时段递推这部分没法完全向量化但可以提前把效率矩阵算好减少循环里的重复运算。最后再分享一个我自己的体会模型不是越复杂越好。对一个刚入门的同学来说建议先把基础算例跑通——只考虑价格型需求响应储能和燃气轮机都用最简单的模型验证粒子群算法确实能收敛并满足功率平衡之后再逐步加入可转移负荷、爬坡约束、SOC循环约束这些复杂因素。我在实际调试中踩过几次坑之后发现大部分跑不出来或者结果诡异的情况根源都不是算法的问题而是模型和数据的物理逻辑没有对清楚。把这个顺序走对了后面的研究会顺畅很多。