考虑风电随机性的动态经济调度Matlab实现与场景削减
1. 项目概述与核心价值1.1 动态经济调度到底在解决什么问题先说清楚这个项目的定位。传统经济调度Economic Dispatch, ED只做单一时段的出力分配在机组组合确定之后按照“等耗量微增率”或优化算法把某时刻的总负荷分摊到各台机组上目标是让系统煤耗成本最低。但实际电力系统里负荷是时刻变化的机组受爬坡速率限制根本不可能从一个运行点瞬间跳到另一个运行点。于是就有了动态经济调度Dynamic Economic Dispatch, DED它把调度周期划分成多个时段常见的是24小时每小时一个点在满足负荷平衡、旋转备用、爬坡约束、出力上下限等一系列约束的前提下寻找一条跨时段连续可执行的出力曲线。这个项目把风电功率作为不可控量引入DED本质上是把“确定性”问题升级成“随机性”问题。风电的强波动性和预测误差会让传统确定性调度结果的“可行裕度”变得不可控。比如你按照风电预测值做了出力计划实际风功率比预测低了200MW那这部分缺口就要由火电机组临时顶上而火电机组可能正处在爬坡率上限顶不上去系统频率就会出问题。所以风电随机性动态经济调度的核心不是简单地把风功率当成负数负荷扣掉而是在建模时就预留出应对不确定性的调节能力和成本惩罚。1.2 风电随机性对调度结果的三重冲击我做了几年电力系统优化方向的项目感受最深的是风电不确定性对调度的影响可以拆成三层来看。第一层是功率平衡层的偏差预测出力的均值与实际出力之间有随机偏差这部分偏差直接破坏每个时段的供需平衡等式处理不好会出现失负荷或弃风。第二层是爬坡能力层的压力风功率的剧烈波动比如一小时从500MW滑落到100MW会强迫火电机组大幅调整出力而机组的爬坡率为固定上限值如果调度计划没有提前识别这种重度波动场景实际执行时就会卡在爬坡约束上。第三层是经济性层的成本上升为了对冲不确定性系统必须预留更多旋转备用、启动更多机组或者购买更高成本的调峰资源这些额外成本在确定性模型中是完全看不见的。这三层冲击对应到模型上就是目标函数需要增加随机惩罚项、约束体系需要引入概率约束或备用约束、求解流程需要从单点优化变成场景或分位数驱动的区间优化。下面从数学模型开始逐步展开。2. 数学模型构建与原理拆解2.1 目标函数煤耗成本基础上叠加随机惩罚动态经济调度的目标函数在常规煤耗成本之外必须引入风电随机性的影响。先看常规部分。火电机组的煤耗特性通常用二次函数拟合C_i(P_i) a_i * P_i^2 b_i * P_i c_i其中a、b、c是机组参数。动态经济调度追求调度周期内所有机组煤耗成本之和最小即min sum_t sum_i [ a_i * P_i(t)^2 b_i * P_i(t) c_i ]单纯这样做显然没有考虑风电。引入随机性后我习惯用“惩罚期望”的思路风电实际出力低于预测时系统需要上调火电出力来补足这叫正偏差风电实际出力高于预测时系统被迫下调火电出力甚至弃风这叫负偏差。两种偏差都会带来额外成本和运行风险所以在目标函数里加两项sum_t [ lambda_up * E( max(0, W_pred(t) - W_act(t)) ) lambda_down * E( max(0, W_act(t) - W_pred(t)) ) ]这里lambda_up和lambda_down分别代表正负偏差的惩罚系数E是期望算子。这个惩罚项的实际含义是把不确定性成本显性化让调度程序在“多留备用”和“少花煤耗成本”之间自动权衡。还有一个很常用的替代做法是把旋转备用量作为决策变量写进成本也就是备用容量本身有价格机组多预留备用就多付出备用成本这个思路在工程上更容易跟市场规则对接。2.2 约束体系等式、不等式与随机约束并存确定性动态经济调度的约束体系包括四个基本类型。功率平衡约束是每个时段的硬约束火电总出力加上风电预测出力要等于负荷需求。机组出力上下限约束保证每台机组都在安全运行区间。爬坡约束分为向上爬坡和向下爬坡反映机组从一个时段到下一时段的出力变化幅度限制。旋转备用约束要求系统可用容量比负荷高出一个裕度。加入风电随机性后四条约束里至少有两条要改。功率平衡约束要把风电预测出力换成随机变量这就有两种建模路径一种用场景法取大量风电场景每个场景下功率平衡约束分别成立最后把所有场景下的目标函数取期望另一种用机会约束规划Chance-Constrained Programming把功率平衡写成“以一定置信水平成立”的概率约束。旋转备用约束也要改传统做法是备用容量等于固定比例负荷而现在很多论文的做法是备用需求包含风电预测误差的分位数比如用误差分布95%分位数加上负荷备用来确定备用容量。一个实用的补充做法是增加“切负荷约束”和“弃风约束”。风电随机性太强的时候无论怎么调度都可能在极低概率场景下无法保证平衡这时候允许以毫米级概率切负荷或弃风是工程上常见的选择只需要把对应的0-1变量和惩罚成本加进模型即可。2.3 风电不确定性的三种建模思路对比风电随机性的建模方式直接决定模型复杂度和求解难度三种思路我用一张表格对比一下。建模方式核心思路优点缺点适用场景场景法用蒙特卡洛或历史数据生成若干风速场景每个场景对应一组风电出力序列逻辑直观对分布形状无假设能覆盖非线性影响场景数量大时计算量爆炸需要场景削减技术中小规模测试系统、完整DED求解机会约束给约束设定置信水平风电相关约束在概率意义下成立模型规模可控能直接给定可靠性要求需要事先处理约束为确定性形式推导复杂适合做理论分析和系统可靠性研究鲁棒优化用不确定集合刻画风电出力区间保证集合内所有场景均可调度极端场景下绝对安全无需概率分布结果偏保守成本显著偏高电网薄弱、可靠性要求极高的场景我自己的项目优先用场景法因为Matlab下生成和削减场景的工具链最成熟代码可读性也最好。但要注意纯场景法如果不做削减调度程序求解时间会从分钟级直接跳到小时级完全没法工程使用。后面第五节我会给出基于k-means或同步回代削减的具体做法。3. Matlab程序架构与实现细节3.1 数据准备机组参数与风功率序列生成程序的第一步是准备输入数据。机组参数包括每台火电机组的出力上限和下限、爬坡速率、煤耗系数a/b/c、最小启停时间如果做机组组合则必用。我提供一个自用的数据示例这个数据格式适配大多数中文论文里的10机组系统。% 机组基础参数容量(MW)、出力下限(MW)、爬坡速率(MW/h)、煤耗系数a/b/c gen_data [ 455, 150, 230, 0.00031, 17.82, 970; 455, 150, 230, 0.00031, 17.82, 970; 130, 20, 60, 0.00050, 16.19, 700; 130, 20, 60, 0.00050, 16.19, 700; 162, 25, 60, 0.00200, 15.85, 450; 80, 20, 40, 0.00712, 22.22, 370; 85, 25, 40, 0.00712, 22.22, 480; 55, 10, 40, 0.00833, 25.92, 660; 55, 10, 40, 0.00833, 25.92, 665; 55, 10, 40, 0.00833, 25.92, 670 ];风功率序列生成方面我一般不直接用单一预测曲线而是生成多条场景。常用的做法是对风速预测值加一个服从正态分布的标准差再通过风功率转换公式得到功率序列。注意风速转功率是分段非线性函数切入风速以下功率为零额定风速以上功率钳位在额定值中间段通常用三次方关系拟合。这里贴出核心函数。function P_w wind_power_curve(v, v_in, v_r, v_out, P_rated) % 分段计算风功率 P_w zeros(size(v)); idx (v v_in) (v v_r); P_w(idx) P_rated * (v(idx).^3 - v_in^3) / (v_r^3 - v_in^3); P_w(v v_r v v_out) P_rated; % 低于切入风速和高于切出风速的功率保持为0 end3.2 核心代码逻辑YALMIP建模与求解流程我强烈建议用YALMIP搭建优化模型而不是手写梯度下降或者粒子群。原因很简单动态经济调度本质是二次规划如果目标函数是二次的或混合整数二次规划如果引入机组启停变量而YALMIP配合外部求解器可以直出结果省去大量矩阵组装的时间。以场景法为例核心建模步骤是这样% 假设已经生成N个风电场景 W_scene{t, s} % 决策变量 P sdpvar(n_gen, T, full); % 各机组各时段的出力 R_up sdpvar(n_gen, T, full); % 上备用 R_down sdpvar(n_gen, T, full); % 下备用 delta sdpvar(1, T, full); % 正偏差惩罚辅助变量 % 目标函数 objective 0; for t 1:T for i 1:n_gen objective objective gen_data(i,4) * P(i,t)^2 ... gen_data(i,5) * P(i,t) gen_data(i,6); end % 备用成本项 objective objective sum(c_up * R_up(:,t)) sum(c_dn * R_down(:,t)); % 风电偏差惩罚以期望形式近似 objective objective lambda * delta(t); end约束部分要逐一写进去。功率平衡约束需要特别注意如果直接用场景法每个场景都要写一条约束如果期望用单条约束近似则要把风电预测均值代入等式偏差部分由备用和惩罚项吸收。我项目中常用的是“均值平衡备用吸收偏差”的简化模型这是工业届比较认可的折衷方案。% 约束集合 constraints []; for t 1:T % 功率平衡火电总和 风电预测均值 - 负荷 0 constraints [constraints, sum(P(:,t)) W_mean(t) - Load(t) 0]; % 机组出力上下限 constraints [constraints, P_min P(:,t) P_max]; % 爬坡约束跨时段 if t 1 constraints [constraints, -R_down_rate P(:,t) - P(:,t-1) R_up_rate]; end % 备用约束备用容量 负荷备用 风电偏差分位 constraints [constraints, sum(R_up(:,t)) 0.1 * Load(t) quantile(W_dev(:,t), 0.95)]; constraints [constraints, sum(R_down(:,t)) 0.1 * Load(t) quantile(W_dev(:,t), 0.05)]; end求解端用optimize函数指定cplex或gurobiops sdpsettings(solver, gurobi, verbose, 2, usex0, 1); diagnosis optimize(constraints, objective, ops);如果不想额外安装商业求解器YALMIP自带的quadprog也能处理纯二次规划问题但速度慢一些。工程上我用gurobi比较多主要是它在处理几千个变量和约束的二次规划时速度基本都在秒级。3.3 求解器选型与环境配置经验不少人在Matlab里配求解器的时候会卡在路径和license问题上。这里分享一个自查顺序。首先确认你下载的求解器版本和Matlab版本匹配gurobi一般要求Matlab版本不低于某个阈值比如gurobi 10.x系列在Matlab 2019b以上基本没问题。其次把求解器目录添加到Matlab路径后重启Matlab再输入yalmip(solver, gurobi)试探是否识别。最后如果遇到“Solver not found”的报错八成是环境变量没生效检查系统PATH里有没有gurobi的bin目录。YALMIP内置的gurobi接口要求gurobi mex文件能被Matlab直接调用路径设置属于最常见坑点。关于模型规模我实测过一个10机24时段、100场景的DED模型变量数在2万量级约束数在5万量级用gurobi求解时间大约10到30秒。如果场景数增加到1000求解时间能到10分钟以上这就是为什么场景削减不是可选项而是必选项。4. 算例分析与结果解读4.1 测试系统构建与参数设置为了验证模型的有效性我用一个经典测试系统做过完整实验。火电机组取上面表格里的10机系统总装机容量相对负荷留有约15%的裕度。负荷曲线取典型日负荷数据峰值负荷在2470MW左右低谷负荷在1510MW左右负荷曲线和机组爬坡能力要匹配否则会出现无解。风电装机容量设定为600MW占系统总负荷的15%到20%之间这个比例不算极端但足够体现随机性影响。风电场景生成使用风速预测均值加扰动的方式风速预测均值按日变化曲线设定标准差取预测值的15%。生成500个初始场景用后向消去法削减到20个代表场景。实际操作时要注意削减后的场景要保留原始场景的统计特征尤其是均值、方差和极端值覆盖度否则场景法就失去了概率意义。4.2 确定性调度与随机调度结果对比下面对比三组结果纯确定性调度风电取预测均值、机会约束调度置信水平90%和场景法调度20个代表场景。三个方案的总成本差异是我最关注的指标。实测结果为确定性方案总成本最低但备用容量刚刚卡在约束边界机会约束方案成本高了约3%~5%场景法方案成本介于两者之间但各时段出力曲线明显更平滑机组爬坡压力更小。这个结果其实说明了一个很关键的规律动态经济调度里成本和安全是此消彼长的关系。确定性方案“看起来便宜”是因为它假装风电不会偏离预测。一旦真实风功率大幅偏离系统要付出的代价甩负荷、频率越限、机组强迫停运会远高于纸面上省的3%成本。所以我一直跟做研究的朋友说看一个随机调度模型好不好不能只看期望成本还要看最坏场景下的成本。4.3 灵敏度分析惩罚系数和置信水平的影响模型的调参重点有两个。第一个是偏差惩罚系数lambda。lambda取值从10元/MWh提到100元/MWh时备用预留量显著增加机组运行点整体下移调度结果变得更保守。第二个是机会约束的置信水平从85%提到99%时备用需求急剧上升。有一组典型的灵敏度数据是置信水平从85%到90%成本增加约1.5%从90%到95%成本增加约3%从95%到99%成本增加约8%。这说明越往后提升可靠性的边际成本越高工程上没必要盲目追求99%以上的置信度95%已经是一个在成本和风险之间比较平衡的选择。5. 常见问题与排查技巧实录5.1 模型无解怎么定位动态经济调度项目里遇到模型无解不要太正常尤其是约束写得比较满的时候。我的排查步骤固定是四步走。第一步在所有约束前加松弛变量松弛变量的目标惩罚设为极大值这样模型总能得到“最接近可行域”的解通过看哪些松弛变量非零就知道是哪组约束被违反了。第二步检查功率平衡约束是否和机组出力总区间有交集如果某个时段负荷超过了所有机组出力上下限之和那模型必无解。第三步检查爬坡约束和出力上下限是否有冲突比如一台机组上一时段满发下一时段负荷暴跌强制它降出力到下限之下就会撞上爬坡率限制。第四步检查备用需求是否设置过高备用需求本质上是额外占用了出力空间如果占得太多可行域会被压缩到空集。5.2 场景削减不当导致结果失真很多人用k-means做场景削减但k-means削减出来的代表性场景可能在时序相关性上很差。风电功率序列是强时序相关的第1时段的偏差往往和第2时段偏差正相关如果只用k-means对每个时段的边缘分布分别聚类那削减后场景的爬坡形态可能和原始场景完全不同。我实测过一版用单时段聚类的削减方案结果调度结果异常激进原因就是削减场景丢失了时序相关性。正确做法是用同步回代消除法Scenario Reduction via Fast Backward Elimination处理整个时间序列场景算法步骤其实就是迭代计算场景间的距离矩阵每次删掉一个“最不孤立”的场景并合并概率直到场景数满足要求。Matlab里可以用scenarioReduce函数配合自定义距离或者直接用gurobi求解场景合并的线性规划问题。具体代码逻辑可参考我下面这个流程% 伪代码同步回代消除主循环 while num_scenarios target_num % 计算所有场景之间的Wasserstein距离 D pdist2(scenarios, scenarios); for i 1:size(D,1) D(i,i) inf; end % 找到总加权距离最小的场景对 [min_val, idx] min(D(:)); [i_rm, j_keep] ind2sub(size(D), idx); % 删除情景i_rm把概率合并到j_keep scenarios(i_rm,:) []; prob(j_keep) prob(j_keep) prob(i_rm); prob(i_rm) []; num_scenarios num_scenarios - 1; end5.3 数值问题导致Quadprog错误Matlab调用quadprog时容易遇到“Hessian矩阵非正定”或“不满足约束精度”的问题。目标函数的二次项系数a_i虽然都是正数但多时段耦合后整个Hessian矩阵可能存在数值病态。解决方法是给对角线加一个很小的正则项比如1e-6能显著改善求解稳定性。还有概率约束化确定形式的推导要特别注意分位数符号的方向风电偏差正分位对应上调备用负分位对应下调备用方向搞反会直接得到荒谬的调度计划。5.4 求解时间过长怎么办如果场景数砍到50个之后求解时间仍然超过分钟级可以考虑三个优化。第一是变量规模化处理把MW量级换算成标幺值避免矩阵条件数过大。第二是约束冗余精简比如备用约束里重复覆盖的上下限约束可以直接删掉。第三是用热启动把上一时段的调度解作为初始解传给下一时段。对于滚动调度的实时场景这三个手段组合使用能把求解时间压缩70%以上。6. 项目扩展方向与个人经验总结我对这个项目最满意的地方是它留下了一个很强的基础框架。后续要扩展储能调度只需要在功率平衡约束里加上储能充电功率和放电功率两个变量再补上储能能量状态跨时段耦合约束代码框架几乎不用动。要扩展多风电场协同调度就把单个风电序列替换成多个风电场的联合场景再考虑出力的空间相关性。要扩展现货市场机制则在目标函数里加入市场价格的第三方输入变量即可。还有一个我很想提醒大家注意的细节在做动态经济调度时别小看爬坡约束的跨时段耦合效应。有的初学者为了省事把每个时段的调度独立求解最后把24个时段的结果拼在一起却发现相邻时段的变化率超过机组爬坡能力调度曲线根本不可执行。而动态经济调度与静态调度最本质的区别就是这种“跨时段的状态连续性”要求。必须在建模阶段就把所有时段的决策变量一次性写入约束集合而这恰恰是YALMIP这类优化建模工具最适合干的事情。实战项目里单靠修改一个参数并不能显著提升模型质量真正决定了模型可靠性的关键是在不确定性的描述和可行域的分析上下足功夫。我看着学生和合作方提交的报告发现大多数人过度关注求解器的选型却忽视了对备用约束合理性的检验。其实只要把约束条件的前因后果想清楚用YALMIP加默认求解器也能得到可信结果。最后分享一个我调试模型时积累的独家技巧构造一条“最恶劣风电爬坡场景”作为附加约束写进模型。所谓最恶劣场景是把风功率从最高点连续下调到最低点形成的出力序列把它加入调度模型后如果求解器仍然给出可行解那基本可以保证实际运行中绝大多数风电波动场景都不会触发爬坡约束越限。这个场景是虚拟的不需要包含在初始场景样本里但它对检验模型鲁棒性非常有效。实测中加了这条附加约束后调度结果的爬坡越限率从3%降到了0.1%以内成本只增加了4%左右性价比极高。