风电接入电网之后调度问题就不再是“给定负荷曲线安排机组启停”这么简单了。风电场站的出力天然带有随机性和波动性你说它是电源吧它又不能像火电机组那样随时听调度指令你说它是负荷吧它又在往系统里注入能量。这种“半可控”属性让传统的确定型经济调度模型越来越吃力——风电大发时可能高估成本、低估备用需求风电骤减时又可能面临切负荷风险。我最近刚好完成了一个基于Matlab的风电随机性动态经济调度模型从建模、场景生成到求解、结果分析完整走了一遍今天把思路和代码细节整理出来给正在做相关课题或者工程落地参考的朋友一些实在的经验。这个模型解决的核心问题非常明确在风电出力不确定的前提下如何安排常规机组火电为主的动态出力计划使得未来若干个调度时段内的总运行成本最低同时保证系统安全稳定运行。动态经济调度的“动态”二字强调的是机组爬坡约束和跨时段耦合而不是简单的静态功率分配而“随机性”则意味着风电不能被当作一个固定的已知量必须在模型里显式处理它的不确定性区间或概率分布。这套模型比较适合电力系统方向的研究生、从事调度自动化系统开发、以及做新能源并网影响分析的工程师参考代码层面基于Matlab YALMIP 商用求解器实现既有教学意义上的可读性也具备工程落地的可行性。1. 模型设计的核心思路与关键点1.1 为什么动态经济调度必须显式处理风电随机性先明确一下电力系统经济调度经历了几个阶段静态经济调度只考虑单一时间断面的负荷分配、动态经济调度考虑多时段爬坡约束、含新能源的经济调度引入风电光伏、随机经济调度把不确定性变量纳入优化模型。我见过很多初学者直接把风电预测值当成确定值塞进约束里这样算出来的结果在预测误差小时没问题一旦风功率实际出力与预测值偏差超过预期就可能出现两种情况一是实际需要的备用容量不足频率波动甚至越限二是调度计划与实际出力不匹配被迫在实时阶段大幅调整机组出力经济性大打折扣。风电随机性对调度的影响不是线性的。举个我在算例里常碰到的现象某时段风功率预测为600 MW但95%置信区间下实际出力可能落在400~800 MW之间。如果把600 MW当作确定值处理系统只要安排600 MW对应的常规机组出力就行但如果实际风功率只有450 MW缺口150 MW需要其他机组马上顶上而备用机组必须有足够的爬坡速率和容量裕度。反过来如果实际风功率到了800 MW常规机组又要快速压出力否则弃风。这就意味着调度计划必须预留“双向”的调节能力而不只是简单按预测值匹配。从数学建模角度风电随机性进入优化模型通常有三条路径场景法Scenario-based、区间法Interval-based和机会约束规划Chance-constrained。场景法把风电出力的概率分布离散成若干个典型场景每个场景对应一组风电出力时序目标函数变成所有场景下成本的期望值区间法只考虑风功率上下界鲁棒性最强但经济性偏保守机会约束则允许以一定置信水平违反约束折中性较好。我这次选用的是场景法因为它在工程实践里最直观、最容易和现有的混合整数线性规划框架结合而且可以通过场景削减控制计算量。1.2 目标函数与约束条件的数学刻画模型的目标函数是让调度周期内所有常规机组的总运行成本期望值最小。这里的成本构成要分三块第一块是煤耗成本一般表示成出力的二次函数但为了用线性规划求解需要分段线性化第二块是机组启停成本启动成本和停机成本分别计费这对应着二进制状态变量也是模型变成混合整数规划的原因第三块是惩罚成本包括弃风惩罚和失负荷惩罚用来在目标函数里“拦住”那些虽然满足约束但会造成严重弃风或切负荷的方案。目标函数写成数学形式大致是[ \min \sum_{s\in S} \rho_s \sum_{t\in T} \sum_{i\in G} \left[ C_i(P_{i,t,s}) SU_i \cdot u_{i,t}^{su} SD_i \cdot u_{i,t}^{sd} c_w \cdot P_{t,s}^{curtail} c_{loss} \cdot P_{t,s}^{loss} \right] ]其中 (\rho_s) 是场景 (s) 的概率权重(S) 是场景集合(T) 是调度时段集合(G) 是常规机组集合。(P_{i,t,s}) 是机组 (i) 在时段 (t) 场景 (s) 下的出力(u_{i,t}^{su}) 和 (u_{i,t}^{sd}) 分别表示启动和停机动作弃风量和失负荷量也被显式建模为决策变量并赋予惩罚系数。这样设计的好处是当模型在极端场景下无法完全平衡功率时它会自动选择代价最小的“软约束”松弛方式而不是直接给出无可行解。约束条件是这个模型真正的灵魂所在我逐条列一下每一条在实际编程时都有对应的实现技巧功率平衡约束(\sum_i P_{i,t,s} P_{t,s}^{wind} - P_{t,s}^{curtail} P_{t,s}^{loss} D_t)对每个时段、每个场景成立。这里风电实际出力是场景给定的弃风和失负荷是松弛变量。机组出力上下限(P_i^{min} \cdot u_{i,t} \le P_{i,t,s} \le P_i^{max} \cdot u_{i,t})注意只有当机组处于开机状态时出力才有意义。爬坡约束(-RD_i \le P_{i,t,s} - P_{i,t-1,s} \le RU_i)这是动态调度区别于静态调度的关键约束跨时段耦合由此产生。最小启停时间约束机组一旦启动必须持续运行至少 (T_{i,on}^{min}) 个时段一旦停机必须持续停机至少 (T_{i,off}^{min}) 个时段。这个约束相当难以处理我后面会给出具体的线性化公式。旋转备用约束(\sum_i \min(RU_i, P_i^{max} - P_{i,t,s}) \ge R_t)要求系统在风电出力突然下降时有足够的向上调节能力。这些约束全部都要对每一个场景、每一个时段建模所以场景数量直接决定模型规模。我最早试过直接生成500个场景加上10台机组、24个时段整个模型大概是几十万行约束YALMIP构建模型花了近一分钟求解器跑起来就更慢了。后来用场景削减把场景数压到20~30个计算时间降了两个数量级而调度结果的期望成本变化不到1%这说明场景削减是非常必要的一步。2. 风电场景生成与削减的实现细节2.1 基于预测误差分布的场景生成方法场景生成的第一步是确定风电出力的概率分布。工程上最常见的做法是假设风电预测误差服从正态分布或者t分布均值取预测值标准差取预测值的某个百分比比如10%~20%。这个假设虽然简化但实际偏差数据统计下来大多数情况下也基本符合至少做研究和初步工程分析是够用的。更精细的做法是用历史数据拟合比如非参数核密度估计或者把误差分成不同风况区间分别建模。我这次用的是拉丁超立方抽样LHS加Cholesky分解处理时段间相关性。为什么不用普通的蒙特卡洛抽样因为蒙特卡洛抽样在样本数不大时容易出现“聚团”现象就是所有场景的风电功率曲线在某个时段都偏低或偏高导致场景代表性差。LHS先把每个时段的风电误差分布均匀分成N个等概率区间再从每个区间里抽取一个样本然后随机配对形成场景这样能保证每个时段的采样覆盖全部分布区间。但LHS有个问题它只保证了单变量单时段的均匀覆盖时段之间仍然是独立抽样的实际风电出力有很强的时间相关性——上一时段风大下一时段风大概率也大。所以还需要用Cholesky分解把独立均匀采样变成符合目标相关系数矩阵的多元正态样本。具体做法是在Matlab里写这么几步用lhsdesign或者手动实现LHS生成维度为 (24 \times N_{scenario}) 的独立均匀分布矩阵 (U)。通过逆累积分布函数转换成标准正态分布矩阵 (Z)Z norminv(U, 0, 1)。设定目标时间相关系数矩阵 (\Sigma)可以用历史风电出力数据计算也可以人工设定做Cholesky分解 (\Sigma L L^T)。计算相关后的正态样本Z_corr L * Z。把 (Z_corr) 反变换回风电误差空间err mu sigma .* Z_corr这里mu取0sigma是预测误差标准差。最终场景功率为P_wind max(0, P_forecast err)注意下限截断到0上限截断到装机容量。这里有个容易踩的坑Cholesky分解要求相关系数矩阵必须是正定矩阵如果直接用样本协方差矩阵有时候会奇异或非正定尤其是时段数多、样本量不够时。解决办法是对相关系数矩阵做特征值修正把所有小于某个小阈值的特征值设为该阈值再重构或者用nearestSPD这类函数找最近的正定矩阵。我在代码里加了这一段处理后续再详细说。2.2 场景削减同步回代消除法的Matlab实现场景生成之后如果直接把全部场景放进模型计算量大到难以接受。场景削减的目的就是从原始 (N) 个场景中挑出 (K) 个最有代表性的场景并且重新分配它们的概率权重使得削减后的场景集合在概率分布意义上最接近原场景集合。我实现的是经典的同步回代消除法Simultaneous Backward Reduction它的核心思想是每次迭代中找出一个场景删除它导致的概率距离增量最小然后把它的概率转移给离它最近的那个场景。重复这个过程直到剩下目标数量的场景。步骤用代码写出来是% 输入scenarios为 N_scen x 24 的风电功率矩阵各场景等概率 1/N % 目标削减到 K 个场景 while size(scenarios, 1) K N size(scenarios, 1); % 计算场景间欧氏距离矩阵 dist pdist2(scenarios, scenarios); % 对每个场景i找它与最近场景j的距离以及对应的j for i 1:N dist(i, i) Inf; % 自身距离设为无穷大 [min_dist(i), nearest(i)] min(dist(i, :)); end % 选择被删除的场景min_dist * 概率 最小的那个 probs ones(N, 1) / N; % 初始等概率 [~, del_idx] min(min_dist .* probs); % 把被删场景的概率加到最近场景上 probs(nearest(del_idx)) probs(nearest(del_idx)) probs(del_idx); % 删除场景 scenarios(del_idx, :) []; probs(del_idx) []; end这个算法虽然简单但我实际用下来有几个细节需要注意距离度量用欧氏距离是最直接的但它把所有时段同等对待这不太合理——峰荷时段的功率平衡误差比低谷时段影响更大如果做精细化场景削减建议用加权距离权重跟该时段的负荷水平或者系统调节能力挂钩。另外削减结果对初始场景数量敏感如果初始场景只有200个而你要削减到5个那代表性会很差我一般是生成1000个场景削减到20个效果比较稳定。还有一点每次迭代重新计算距离矩阵是 (O(N^2)) 的复杂度初始场景多的时候会有点慢好在Matlab的pdist2是编译过的1000x24的矩阵迭代下来大概也就几十秒完全可以接受。2.3 场景质量评估削减前后分布对比场景削减做完之后务必要做一个质量评估不然你不知道削减后的场景集合是不是把关键信息丢掉了。我的标准做法是计算两个误差指标一是场景削减前后的风电出力期望值曲线偏差二是方差偏差。期望值偏差控制在1%以内基本就没问题方差偏差稍微大一点没关系但超过5%就要警惕削减过头了。exp_before mean(scenarios_origin, 1); exp_after sum(scenarios_reduced .* probs_reduced, 1); var_before var(scenarios_origin, 1); var_after sum((scenarios_reduced - exp_after).^2 .* probs_reduced, 1); fprintf(期望值最大偏差: %.4f MW\n, max(abs(exp_before - exp_after))); fprintf(方差最大偏差: %.4f (MW)^2\n, max(abs(var_before - var_after)));我做过一组测试1000个原始场景削减到20个期望值偏差不到0.3 MW方差偏差约2.8%而调度目标函数值与用500个场景直接求解的结果相比只差了0.5%左右这个精度对工程决策来说是完全可以接受的。所以强烈建议朋友们在做场景法调度时不要嫌麻烦跳过场景削减这一步也不要凭感觉定场景数用数据说话。3. Matlab代码实现的框架与核心函数拆解3.1 整体代码架构设计整个项目我按模块化思路组织每个功能块独立成函数主脚本只负责读取参数、调用模块、输出结果。这样做的好处是后续要换测试系统、换求解器或者改场景生成方式只需要动对应模块就行。目录结构大概是wind_stochastic_eed/ ├── main_script.m % 主脚本入口 ├── data/ │ ├── system_data.m % 10机24节点系统参数 │ └── wind_data.m % 风功率预测值与装机容量 ├── scenario/ │ ├── gen_scenarios.m % 场景生成LHS Cholesky │ ├── reduce_scenarios.m % 场景削减同步回代消除 │ └── eval_scenarios.m % 场景质量评估 ├── model/ │ ├── build_eed_model.m % 构建优化模型YALMIP │ └── linearize_cost.m % 煤耗成本分段线性化 ├── solve/ │ └── solve_model.m % 调用求解器并处理结果 └── plot/ └── plot_results.m % 绘制调度结果对比图主脚本的逻辑很直接先载入系统参数然后调用场景生成和削减模块接着构建优化模型求解最后画图分析。我写代码的习惯是在主脚本里用tic/toc记录每段耗时这样能很清楚看到瓶颈在哪一步。3.2 YALMIP建模的核心代码逻辑模型构建是整套代码最核心的部分。我用的建模工具是YALMIP它最大的优势是把复杂的混合整数规划约束用接近数学表达式的语法写出来可读性极强而且底层可以无缝切换CPLEX、Gurobi、MOSEK等求解器。如果你还在用intlinprog手工拼矩阵我强烈建议试试YALMIP尤其是做多场景、多时段耦合的模型手拼约束矩阵的维护成本实在太高了。核心建模代码框架大概是这样的function model build_eed_model(sys, wind_scen, weights, params) T params.T; % 时段数 G sys.num_units; % 机组数 S length(weights); % 场景数 % 决策变量出力 P开停机状态 u启动/停机动作 v_start/v_stop P sdpvar(G, T, S, full); % 机组出力 u binvar(G, T, full); % 开机状态0/1 v_start binvar(G, T, full); % 启动动作 v_stop binvar(G, T, full); % 停机动作 % 软约束松弛变量 curt sdpvar(1, T, S, full); % 弃风量 loss sdpvar(1, T, S, full); % 失负荷量 % 目标函数初始化 objective 0; constraints []; for s 1:S for t 1:T for i 1:G % 煤耗成本分段线性化后用罚函数累加 objective objective weights(s) * ... (linear_fuel_cost(P(i,t,s), sys.fuel(i,:), params.C_penalty)); end % 弃风/失负荷惩罚 objective objective weights(s) * ... (params.c_curtail * curt(1,t,s) params.c_loss * loss(1,t,s)); % 功率平衡约束 constraints [constraints, ... sum(P(:,t,s)) wind_scen(s,t) - curt(1,t,s) loss(1,t,s) sys.load(t)]; end end % 机组出力上下限与爬坡约束 for i 1:G for t 1:T for s 1:S constraints [constraints, ... sys.Pmin(i) * u(i,t) P(i,t,s) sys.Pmax(i) * u(i,t)]; end if t 1 for s 1:S constraints [constraints, ... P(i,t,s) - P(i,t-1,s) sys.RU(i)]; constraints [constraints, ... P(i,t-1,s) - P(i,t,s) sys.RD(i)]; end end end end % 最小启停时间约束的线性化表达 constraints [constraints, min_up_time_constraints(u, v_start, v_stop, sys)]; model.objective objective; model.constraints constraints; model.P P; model.u u; model.v_start v_start; model.v_stop v_stop; model.curt curt; model.loss loss; end这里需要重点解释两个地方。第一(P) 变量的维度是 (G \times T \times S)也就是说每个场景下每台机组每个时段都有一个独立的出力变量场景间的耦合完全通过目标函数里概率加权来实现约束条件本身不强耦合场景。这种“场景解耦”的结构是场景法能用分解算法提速的理论基础。第二为什么 (u) 变量不随场景变化因为机组的开停机计划是由调度员在日前阶段就制定好的不能随风电场景变化而变化这是非预期性non-anticipativity约束的要求——调度决策只能基于预测信息不能在看到实际风电后再去做决定。这一点初学者特别容易搞混如果把 (u) 也定义成 (G \times T \times S)相当于模型获得了“预知未来”的能力算出来的结果在现实中根本没法执行。3.3 最小启停时间约束怎么线性化最小启停时间约束是让很多人头疼的地方。这个约束的本质是如果机组在时段 (t) 启动则它在 (t, t1, ..., tT_{on}^{min}-1) 时段内都必须保持开机如果机组在时段 (t) 停机则它在 (t, t1, ..., tT_{off}^{min}-1) 时段内都必须保持停机。线性化的经典做法是引入启动动作变量 (v_{i,t}^{su} 1) 表示机组 (i) 在时段 (t) 开始启动停机动作变量 (v_{i,t}^{sd} 1) 表示开始停机然后建立状态与动作之间的逻辑关系[ u_{i,t} - u_{i,t-1} v_{i,t}^{su} - v_{i,t}^{sd} ]这条约束把“状态切换”和“动作发生”联系起来开机状态下启动动作取1停机状态下动作取0。然后最小开机时间约束写成[ \sum_{kt}^{tT_{on,i}^{min}-1} u_{i,k} \ge T_{on,i}^{min} \cdot v_{i,t}^{su}, \quad \forall t ]最小停机时间约束同理[ \sum_{kt}^{tT_{off,i}^{min}-1} (1 - u_{i,k}) \ge T_{off,i}^{min} \cdot v_{i,t}^{sd}, \quad \forall t ]这两条约束的直观理解是如果机组在 (t) 时段启动了(v^{su}1)那么未来 (T_{on}^{min}) 个时段里必须全部开机否则左边求和会小于右边约束被违反。注意约束里 (k) 的索引超过调度周期末端时要截断处理否则下标越界或者引入虚假约束这个问题我记得第一次跑的时候就把T弄错了导致边界时段无解。3.4 煤耗成本的分段线性化处理火电机组的煤耗成本通常是出力的二次函数比如 (C_i(P) a_i P^2 b_i P c_i)。YALMIP本身支持二次目标函数CPLEX和Gurobi也都能直接求解MIQP混合整数二次规划理论上可以直接写。但我建议还是做分段线性化原因有三一是MIQP的求解速度比MILP慢不少场景数多时差距显著二是分段线性化在工程中更通用有些求解器不带MIQP功能三是分段线性化的成本曲线对结果精度的影响完全可以控制在0.5%以内没必要为这个牺牲计算性能。分段线性化的做法是把出力范围 ([P^{min}, P^{max}]) 分成若干段每段首尾相连用线性插值逼近原二次函数。Matlab代码实现如下function cost linear_fuel_cost(P, fuel_coeff, seg_points) % fuel_coeff [a, b, c] 二次项、一次项、常数项 % seg_points: 分段断点如 [Pmin, P1, P2, ..., Pmax] a fuel_coeff(1); b fuel_coeff(2); c fuel_coeff(3); n_seg length(seg_points) - 1; % 计算每个断点的真实成本 cost_pts a * seg_points.^2 b * seg_points c; % 生成分段线性插值YALMIP的pwlp函数或者手动实现 % 这里推荐用YALMIP内置的pwlppiecewise linear programming cost pwlp(P, seg_points, cost_pts); endpwlp是YALMIP里处理分段线性函数的内置方法使用起来很简洁它内部会自动引入二进制变量和辅助变量来建模分段区间选择。如果不用YALMIP手写分段线性化也有标准做法把出力决策变量拆成各分段区间上的子变量之和每个子变量限幅落在该分段区间内再引入二进制变量保证相邻分段按顺序填充。这个手写过程容易出错我在这里就不展开了有兴趣的朋友可以找线性规划的建模资料看看。4. 求解策略与测试算例分析4.1 求解器选择与配置经验这套模型最终是混合整数线性规划所以求解器必须支持MILP。我配置的是YALMIP Gurobi的组合如果没装GurobiCPLEX效果也差不多。为什么要用商用求解器说实话Matlab自带的intlinprog也能跑但求解这种多场景多时段的混合整数规划问题中小规模还行规模一大就明显慢而且数值稳定性不够好有时候同一个模型多跑几次结果都可能不一样虽然这更多是容差设置的问题。商用求解器的预求解presolve、割平面加速、启发式分支策略都要成熟很多同一个模型用intlinprog可能跑10分钟还在挣扎Gurobi 30秒内给出接近最优的可行解差距不是一点半点。配置YALMIP Gurobi要注意版本匹配问题特别是Matlab版本和Gurobi的接口版本要对应否则yalmip(clear)之后solvesdp大概率报找不到求解器的错误。我常用的配置方式是% 在Matlab中设置Gurobi路径 addpath(C:\gurobi1100\win64\matlab); gurobi_setup; % Gurobi自带的Matlab接口初始化脚本 % 检查YALMIP是否能找到Gurobi yalmiptest;如果yalmiptest显示Gurobi状态为available那就没问题了。求解时记得设置求解参数我常用的几个参数是options sdpsettings(solver, gurobi, ... verbose, 2, ... % 输出求解日志 gurobi.MIPGap, 0.01, ... % 最优性间隙设为1% gurobi.TimeLimit, 300, ... % 时间限制5分钟 gurobi.Threads, 8); % 并行线程数MIPGap的设置很关键。如果你设成0求解器可能需要非常久才能证明最优性设成1%或2%通常能在几十秒内给出一个非常接近最优的可行解而且目标函数值的差距通常小于0.3%完全够用。4.2 测试系统与参数设置我用的测试系统是经典的10机组系统24个调度时段负荷曲线采用典型夏季日负荷数据峰值负荷在时段19约2100 MW谷值在时段4约1200 MW。风电装机容量设为600 MW占系统峰荷比例约29%这个渗透率在当前很多实际电网中已经相当常见了算出来的结果有参考价值。风光预测曲线我默认给定一个典型形状夜间风大、白天风小、傍晚略有回升——这是我根据实际风电出力的“反调峰”特性设计的恰好和负荷曲线形成鲜明对比。预测误差标准差设为预测值的15%用LHS方法生成1000个原始场景然后削减到20个代表场景。置信度约束那块我暂时用的是软惩罚方式先不收窄场景集合而是让失负荷惩罚系数足够大使模型自动避免失负荷场景。4.3 确定性调度与随机性调度的结果对比为了直观展示随机性建模的价值我做了三组对比实验方案A确定性只考虑风电预测值按一个固定确定性场景建模。方案B随机性风电出力用20个削减场景描述目标函数为期望成本最小。方案C鲁棒性只考虑风电出力波动区间的最坏情况目标函数为最坏场景下成本最小。三组实验的模型参数完全一致只改风电出力的建模方式。结果如下指标方案A确定性方案B随机性方案C鲁棒期望总成本万元126.4118.2137.8失负荷概率4.7%0.3%0%弃风率0.2%0.1%8.5%计算时间秒8.242.615.3方案A的期望成本比方案B高了不少原因在于确定性模型因为“看到”的风电出力偏高安排了较少的常规机组开机导致某些场景下需要调用昂贵的高爬坡率机组来补缺口或者失负荷产生巨额惩罚。方案C虽然彻底消除了失负荷风险但代价是大量弃风——因为方案C认为风功率可能随时跌到最低区间为了保险让常规机组踩满出力结果在实际风功率正常偏高的场景下只能弃风。方案B在这两个极端之间做了平衡它不追求绝对安全而是追求“期望成本最优”这正是随机规划的哲学。再单独看一下方案B的机组出力曲线我发现一个有意思的现象风电预测出力较高的时段比如时段3~5模型仍然保留了一台小容量快速机组待命没有因为它“大概率不需要”就把它停机。这是因为场景集中存在风电偏低的情形保留这台快速机组的额外煤耗成本远小于万一需要它而它没法在短时间内启动造成的失负荷惩罚。这就是随机调度区别于确定性调度的本质特征——它不是用预测值做决策而是用整个概率分布做决策。4.4 场景数量与置信水平对结果的影响最后说一下场景数对求解质量和速度的影响。我扫了一组场景数从5到100的测试场景数期望成本万元求解时间秒相对偏差%5119.83.21.3510118.98.50.5920118.242.60.0850118.1187.40.00100118.1超过10分钟—可以看到场景数从5增加到20期望成本下降了约1.3%但这个差距在工程上也不算特别大从20增加到50成本变化只有0.08%基本可以忽略。考虑到计算时间的爆炸式增长20个场景在这套测试系统下是性价比最好的选择。当然这个结论跟系统规模、风电渗透率、机组数量都有关系具体项目里还是要做一次敏感性扫描来定场景数不能盲目照抄。5. 常见问题与排查实录5.1 求解器报“无可行解”怎么排查这个错误是最常见的我几乎每次改动模型参数都会遇到几次。排查思路我整理成一个标准流程第一先检查松约束。把硬约束改成软约束看看目标函数里的惩罚项是不是有值——如果有值说明问题出在约束过紧与目标函数惩罚的权衡之间如果还是无可行解那就是约束本身逻辑错误。第二检查功率平衡系数和负荷数据量纲。我踩过一次很蠢的坑负荷数据单位是MW风电数据单位是kW二者差了三个数量级平衡约束永远不成立模型当然无解。第三检查最小启停时间约束的初始状态匹配问题。如果系统数据里给定的是机组在调度周期前的初始状态比如开机时长不够导致强制继续开机而你的模型没有正确建模初始状态就可能出现矛盾约束。解决办法是专门为前几个时段补上初始状态约束比如 (u_{i,1} 1) 或 (\sum_{k1}^{T_i^{min}-T_0} u_{i,k} \ge ...)。第四最直接的方法把模型里的二进制变量全部松弛为连续变量在YALMIP里把binvar改成sdpvar如果松弛后仍无可行解那问题出在线性约束本身如果松弛后有解但整数模型无解那问题出在整数的组合约束上大概率是最小启停时间或者分段线性化的辅助变量约束写错了。5.2 场景生成中的正态性假设失效怎么办我之前说过风电预测误差很多时候用正态分布近似是够用的但在极端天气情况下比如台风前、辐合带云系过境误差分布会出现明显的重尾或者双峰特征此时单纯的正态假设会低估极端场景的发生概率导致调度计划在大风极端场景下失衡。如果数据足够我建议用历史误差数据直接做经验分布采样具体操作是把实测误差按分位数排布然后用datasample加权重直接抽样绕开分布拟合这一步。如果数据不足退而求其次用t分布自由度5~8替代正态分布t分布的尾巴更厚对极端事件的覆盖更保守。我做敏感性分析时还会额外做一组“压力测试”人为把某个风电出力极低的场景概率权重从1%调到10%观察目标函数和开机计划是否发生显著改变如果变化太大说明模型的鲁棒性不足这时候要稍微调大旋转备用约束或者调低置信度。5.3 求解速度太慢的优化套路最后送几个让MILP求解速度提升的经验这些都是我在实际调参过程中验证过的一是初始化可行解warm start。先用确定性场景预测值快速求解得到一组机组开停机计划然后把这组 (u) 作为随机模型的初始解传给求解器。商用求解器非常擅长在好的人工解基础上做改进这往往能把求解时间压缩一半以上。YALMIP里可以通过assign和sdpsettings(gurobi.StartNodeLimit, ...)来设置初始解。二是善用对称性破缺symmetry breaking。10机组测试系统里有两三台类型相同的机组场景模型对相同机组做对称的调度决策会让分支定界过程大量重复搜索等价节点。通过添加机组优先顺序约束比如 (P_1 \ge P_2 \ge P_3)对于同型机组可以显著收紧搜索空间。三是调整求解器的容差参数。MIPGap从默认的万分之一放宽到百分之一求解时间往往能降低一个数量级而经济效益损失不到千分之三。在工程决策场景中这个交换是非常划算的。记住这个原则先跑一版严格最优解确认模型正确后续批量做敏感性分析时就放宽容差加速效率和安全兼顾。四是考虑用Benders分解或拉格朗日松弛做大规模场景加速但这就是另一个层面的算法优化了。对于大多数教学场景和中型算例单机MILP已经足够不用一开始就上分解算法避免把简单问题复杂化。6. 这套模型的扩展方向与实际应用价值最后再说几句扩展方向。这套随机动态经济调度模型框架本身是“骨架”风电只是第一个实例化的不确定源。分布式光伏、电动汽车负荷、储能系统调度它们的随机性和波动性本质上是同一类问题——无非是改掉不确定变量的概率分布和参数边界再把储能充放电约束加进去。我在实际项目里就把这套框架扩展成了风光储联合调度模型做法是在目标函数里新增储能充放电的损耗成本和充放电功率约束风电和光伏各生成一组场景用联合场景矩阵代入模型求解。扩展过程很顺滑这得益于底层YALMIP建模的高可读性和模块化代码结构。另一个值得尝试的方向是反馈调度rolling horizon。日前调度做了24小时计划实际执行时风电预测每4小时刷新一次与其固定执行已过时的日前计划不如每4小时重跑一次滚动调度只执行未来8~12小时的部分计划。这种模型预测控制MPC式的调度方式和动态经济调度在数学框架上是完全兼容的只多了一步“每次求解后只取前几个时段的解”这也是从学术模型走向工程落地的关键一步。我个人在写完这套模型后最大的感触是随机调度的重点不在算法本身有多复杂而在于你是否准确理解了“决策必须在不确定性被揭示之前做出”这一非预期性约束。很多跑出错误结果的项目翻来覆去检查也找不到问题最后发现就是某个决策变量被错误地定义成了场景相关导致模型拥有了现实世界不可能有的“预知能力”。希望这篇文章能帮大家避开这个最大的坑也欢迎有同行在类似问题上一起来交流心得。
基于Matlab的风电随机性动态经济调度模型设计与实现
风电接入电网之后调度问题就不再是“给定负荷曲线安排机组启停”这么简单了。风电场站的出力天然带有随机性和波动性你说它是电源吧它又不能像火电机组那样随时听调度指令你说它是负荷吧它又在往系统里注入能量。这种“半可控”属性让传统的确定型经济调度模型越来越吃力——风电大发时可能高估成本、低估备用需求风电骤减时又可能面临切负荷风险。我最近刚好完成了一个基于Matlab的风电随机性动态经济调度模型从建模、场景生成到求解、结果分析完整走了一遍今天把思路和代码细节整理出来给正在做相关课题或者工程落地参考的朋友一些实在的经验。这个模型解决的核心问题非常明确在风电出力不确定的前提下如何安排常规机组火电为主的动态出力计划使得未来若干个调度时段内的总运行成本最低同时保证系统安全稳定运行。动态经济调度的“动态”二字强调的是机组爬坡约束和跨时段耦合而不是简单的静态功率分配而“随机性”则意味着风电不能被当作一个固定的已知量必须在模型里显式处理它的不确定性区间或概率分布。这套模型比较适合电力系统方向的研究生、从事调度自动化系统开发、以及做新能源并网影响分析的工程师参考代码层面基于Matlab YALMIP 商用求解器实现既有教学意义上的可读性也具备工程落地的可行性。1. 模型设计的核心思路与关键点1.1 为什么动态经济调度必须显式处理风电随机性先明确一下电力系统经济调度经历了几个阶段静态经济调度只考虑单一时间断面的负荷分配、动态经济调度考虑多时段爬坡约束、含新能源的经济调度引入风电光伏、随机经济调度把不确定性变量纳入优化模型。我见过很多初学者直接把风电预测值当成确定值塞进约束里这样算出来的结果在预测误差小时没问题一旦风功率实际出力与预测值偏差超过预期就可能出现两种情况一是实际需要的备用容量不足频率波动甚至越限二是调度计划与实际出力不匹配被迫在实时阶段大幅调整机组出力经济性大打折扣。风电随机性对调度的影响不是线性的。举个我在算例里常碰到的现象某时段风功率预测为600 MW但95%置信区间下实际出力可能落在400~800 MW之间。如果把600 MW当作确定值处理系统只要安排600 MW对应的常规机组出力就行但如果实际风功率只有450 MW缺口150 MW需要其他机组马上顶上而备用机组必须有足够的爬坡速率和容量裕度。反过来如果实际风功率到了800 MW常规机组又要快速压出力否则弃风。这就意味着调度计划必须预留“双向”的调节能力而不只是简单按预测值匹配。从数学建模角度风电随机性进入优化模型通常有三条路径场景法Scenario-based、区间法Interval-based和机会约束规划Chance-constrained。场景法把风电出力的概率分布离散成若干个典型场景每个场景对应一组风电出力时序目标函数变成所有场景下成本的期望值区间法只考虑风功率上下界鲁棒性最强但经济性偏保守机会约束则允许以一定置信水平违反约束折中性较好。我这次选用的是场景法因为它在工程实践里最直观、最容易和现有的混合整数线性规划框架结合而且可以通过场景削减控制计算量。1.2 目标函数与约束条件的数学刻画模型的目标函数是让调度周期内所有常规机组的总运行成本期望值最小。这里的成本构成要分三块第一块是煤耗成本一般表示成出力的二次函数但为了用线性规划求解需要分段线性化第二块是机组启停成本启动成本和停机成本分别计费这对应着二进制状态变量也是模型变成混合整数规划的原因第三块是惩罚成本包括弃风惩罚和失负荷惩罚用来在目标函数里“拦住”那些虽然满足约束但会造成严重弃风或切负荷的方案。目标函数写成数学形式大致是[ \min \sum_{s\in S} \rho_s \sum_{t\in T} \sum_{i\in G} \left[ C_i(P_{i,t,s}) SU_i \cdot u_{i,t}^{su} SD_i \cdot u_{i,t}^{sd} c_w \cdot P_{t,s}^{curtail} c_{loss} \cdot P_{t,s}^{loss} \right] ]其中 (\rho_s) 是场景 (s) 的概率权重(S) 是场景集合(T) 是调度时段集合(G) 是常规机组集合。(P_{i,t,s}) 是机组 (i) 在时段 (t) 场景 (s) 下的出力(u_{i,t}^{su}) 和 (u_{i,t}^{sd}) 分别表示启动和停机动作弃风量和失负荷量也被显式建模为决策变量并赋予惩罚系数。这样设计的好处是当模型在极端场景下无法完全平衡功率时它会自动选择代价最小的“软约束”松弛方式而不是直接给出无可行解。约束条件是这个模型真正的灵魂所在我逐条列一下每一条在实际编程时都有对应的实现技巧功率平衡约束(\sum_i P_{i,t,s} P_{t,s}^{wind} - P_{t,s}^{curtail} P_{t,s}^{loss} D_t)对每个时段、每个场景成立。这里风电实际出力是场景给定的弃风和失负荷是松弛变量。机组出力上下限(P_i^{min} \cdot u_{i,t} \le P_{i,t,s} \le P_i^{max} \cdot u_{i,t})注意只有当机组处于开机状态时出力才有意义。爬坡约束(-RD_i \le P_{i,t,s} - P_{i,t-1,s} \le RU_i)这是动态调度区别于静态调度的关键约束跨时段耦合由此产生。最小启停时间约束机组一旦启动必须持续运行至少 (T_{i,on}^{min}) 个时段一旦停机必须持续停机至少 (T_{i,off}^{min}) 个时段。这个约束相当难以处理我后面会给出具体的线性化公式。旋转备用约束(\sum_i \min(RU_i, P_i^{max} - P_{i,t,s}) \ge R_t)要求系统在风电出力突然下降时有足够的向上调节能力。这些约束全部都要对每一个场景、每一个时段建模所以场景数量直接决定模型规模。我最早试过直接生成500个场景加上10台机组、24个时段整个模型大概是几十万行约束YALMIP构建模型花了近一分钟求解器跑起来就更慢了。后来用场景削减把场景数压到20~30个计算时间降了两个数量级而调度结果的期望成本变化不到1%这说明场景削减是非常必要的一步。2. 风电场景生成与削减的实现细节2.1 基于预测误差分布的场景生成方法场景生成的第一步是确定风电出力的概率分布。工程上最常见的做法是假设风电预测误差服从正态分布或者t分布均值取预测值标准差取预测值的某个百分比比如10%~20%。这个假设虽然简化但实际偏差数据统计下来大多数情况下也基本符合至少做研究和初步工程分析是够用的。更精细的做法是用历史数据拟合比如非参数核密度估计或者把误差分成不同风况区间分别建模。我这次用的是拉丁超立方抽样LHS加Cholesky分解处理时段间相关性。为什么不用普通的蒙特卡洛抽样因为蒙特卡洛抽样在样本数不大时容易出现“聚团”现象就是所有场景的风电功率曲线在某个时段都偏低或偏高导致场景代表性差。LHS先把每个时段的风电误差分布均匀分成N个等概率区间再从每个区间里抽取一个样本然后随机配对形成场景这样能保证每个时段的采样覆盖全部分布区间。但LHS有个问题它只保证了单变量单时段的均匀覆盖时段之间仍然是独立抽样的实际风电出力有很强的时间相关性——上一时段风大下一时段风大概率也大。所以还需要用Cholesky分解把独立均匀采样变成符合目标相关系数矩阵的多元正态样本。具体做法是在Matlab里写这么几步用lhsdesign或者手动实现LHS生成维度为 (24 \times N_{scenario}) 的独立均匀分布矩阵 (U)。通过逆累积分布函数转换成标准正态分布矩阵 (Z)Z norminv(U, 0, 1)。设定目标时间相关系数矩阵 (\Sigma)可以用历史风电出力数据计算也可以人工设定做Cholesky分解 (\Sigma L L^T)。计算相关后的正态样本Z_corr L * Z。把 (Z_corr) 反变换回风电误差空间err mu sigma .* Z_corr这里mu取0sigma是预测误差标准差。最终场景功率为P_wind max(0, P_forecast err)注意下限截断到0上限截断到装机容量。这里有个容易踩的坑Cholesky分解要求相关系数矩阵必须是正定矩阵如果直接用样本协方差矩阵有时候会奇异或非正定尤其是时段数多、样本量不够时。解决办法是对相关系数矩阵做特征值修正把所有小于某个小阈值的特征值设为该阈值再重构或者用nearestSPD这类函数找最近的正定矩阵。我在代码里加了这一段处理后续再详细说。2.2 场景削减同步回代消除法的Matlab实现场景生成之后如果直接把全部场景放进模型计算量大到难以接受。场景削减的目的就是从原始 (N) 个场景中挑出 (K) 个最有代表性的场景并且重新分配它们的概率权重使得削减后的场景集合在概率分布意义上最接近原场景集合。我实现的是经典的同步回代消除法Simultaneous Backward Reduction它的核心思想是每次迭代中找出一个场景删除它导致的概率距离增量最小然后把它的概率转移给离它最近的那个场景。重复这个过程直到剩下目标数量的场景。步骤用代码写出来是% 输入scenarios为 N_scen x 24 的风电功率矩阵各场景等概率 1/N % 目标削减到 K 个场景 while size(scenarios, 1) K N size(scenarios, 1); % 计算场景间欧氏距离矩阵 dist pdist2(scenarios, scenarios); % 对每个场景i找它与最近场景j的距离以及对应的j for i 1:N dist(i, i) Inf; % 自身距离设为无穷大 [min_dist(i), nearest(i)] min(dist(i, :)); end % 选择被删除的场景min_dist * 概率 最小的那个 probs ones(N, 1) / N; % 初始等概率 [~, del_idx] min(min_dist .* probs); % 把被删场景的概率加到最近场景上 probs(nearest(del_idx)) probs(nearest(del_idx)) probs(del_idx); % 删除场景 scenarios(del_idx, :) []; probs(del_idx) []; end这个算法虽然简单但我实际用下来有几个细节需要注意距离度量用欧氏距离是最直接的但它把所有时段同等对待这不太合理——峰荷时段的功率平衡误差比低谷时段影响更大如果做精细化场景削减建议用加权距离权重跟该时段的负荷水平或者系统调节能力挂钩。另外削减结果对初始场景数量敏感如果初始场景只有200个而你要削减到5个那代表性会很差我一般是生成1000个场景削减到20个效果比较稳定。还有一点每次迭代重新计算距离矩阵是 (O(N^2)) 的复杂度初始场景多的时候会有点慢好在Matlab的pdist2是编译过的1000x24的矩阵迭代下来大概也就几十秒完全可以接受。2.3 场景质量评估削减前后分布对比场景削减做完之后务必要做一个质量评估不然你不知道削减后的场景集合是不是把关键信息丢掉了。我的标准做法是计算两个误差指标一是场景削减前后的风电出力期望值曲线偏差二是方差偏差。期望值偏差控制在1%以内基本就没问题方差偏差稍微大一点没关系但超过5%就要警惕削减过头了。exp_before mean(scenarios_origin, 1); exp_after sum(scenarios_reduced .* probs_reduced, 1); var_before var(scenarios_origin, 1); var_after sum((scenarios_reduced - exp_after).^2 .* probs_reduced, 1); fprintf(期望值最大偏差: %.4f MW\n, max(abs(exp_before - exp_after))); fprintf(方差最大偏差: %.4f (MW)^2\n, max(abs(var_before - var_after)));我做过一组测试1000个原始场景削减到20个期望值偏差不到0.3 MW方差偏差约2.8%而调度目标函数值与用500个场景直接求解的结果相比只差了0.5%左右这个精度对工程决策来说是完全可以接受的。所以强烈建议朋友们在做场景法调度时不要嫌麻烦跳过场景削减这一步也不要凭感觉定场景数用数据说话。3. Matlab代码实现的框架与核心函数拆解3.1 整体代码架构设计整个项目我按模块化思路组织每个功能块独立成函数主脚本只负责读取参数、调用模块、输出结果。这样做的好处是后续要换测试系统、换求解器或者改场景生成方式只需要动对应模块就行。目录结构大概是wind_stochastic_eed/ ├── main_script.m % 主脚本入口 ├── data/ │ ├── system_data.m % 10机24节点系统参数 │ └── wind_data.m % 风功率预测值与装机容量 ├── scenario/ │ ├── gen_scenarios.m % 场景生成LHS Cholesky │ ├── reduce_scenarios.m % 场景削减同步回代消除 │ └── eval_scenarios.m % 场景质量评估 ├── model/ │ ├── build_eed_model.m % 构建优化模型YALMIP │ └── linearize_cost.m % 煤耗成本分段线性化 ├── solve/ │ └── solve_model.m % 调用求解器并处理结果 └── plot/ └── plot_results.m % 绘制调度结果对比图主脚本的逻辑很直接先载入系统参数然后调用场景生成和削减模块接着构建优化模型求解最后画图分析。我写代码的习惯是在主脚本里用tic/toc记录每段耗时这样能很清楚看到瓶颈在哪一步。3.2 YALMIP建模的核心代码逻辑模型构建是整套代码最核心的部分。我用的建模工具是YALMIP它最大的优势是把复杂的混合整数规划约束用接近数学表达式的语法写出来可读性极强而且底层可以无缝切换CPLEX、Gurobi、MOSEK等求解器。如果你还在用intlinprog手工拼矩阵我强烈建议试试YALMIP尤其是做多场景、多时段耦合的模型手拼约束矩阵的维护成本实在太高了。核心建模代码框架大概是这样的function model build_eed_model(sys, wind_scen, weights, params) T params.T; % 时段数 G sys.num_units; % 机组数 S length(weights); % 场景数 % 决策变量出力 P开停机状态 u启动/停机动作 v_start/v_stop P sdpvar(G, T, S, full); % 机组出力 u binvar(G, T, full); % 开机状态0/1 v_start binvar(G, T, full); % 启动动作 v_stop binvar(G, T, full); % 停机动作 % 软约束松弛变量 curt sdpvar(1, T, S, full); % 弃风量 loss sdpvar(1, T, S, full); % 失负荷量 % 目标函数初始化 objective 0; constraints []; for s 1:S for t 1:T for i 1:G % 煤耗成本分段线性化后用罚函数累加 objective objective weights(s) * ... (linear_fuel_cost(P(i,t,s), sys.fuel(i,:), params.C_penalty)); end % 弃风/失负荷惩罚 objective objective weights(s) * ... (params.c_curtail * curt(1,t,s) params.c_loss * loss(1,t,s)); % 功率平衡约束 constraints [constraints, ... sum(P(:,t,s)) wind_scen(s,t) - curt(1,t,s) loss(1,t,s) sys.load(t)]; end end % 机组出力上下限与爬坡约束 for i 1:G for t 1:T for s 1:S constraints [constraints, ... sys.Pmin(i) * u(i,t) P(i,t,s) sys.Pmax(i) * u(i,t)]; end if t 1 for s 1:S constraints [constraints, ... P(i,t,s) - P(i,t-1,s) sys.RU(i)]; constraints [constraints, ... P(i,t-1,s) - P(i,t,s) sys.RD(i)]; end end end end % 最小启停时间约束的线性化表达 constraints [constraints, min_up_time_constraints(u, v_start, v_stop, sys)]; model.objective objective; model.constraints constraints; model.P P; model.u u; model.v_start v_start; model.v_stop v_stop; model.curt curt; model.loss loss; end这里需要重点解释两个地方。第一(P) 变量的维度是 (G \times T \times S)也就是说每个场景下每台机组每个时段都有一个独立的出力变量场景间的耦合完全通过目标函数里概率加权来实现约束条件本身不强耦合场景。这种“场景解耦”的结构是场景法能用分解算法提速的理论基础。第二为什么 (u) 变量不随场景变化因为机组的开停机计划是由调度员在日前阶段就制定好的不能随风电场景变化而变化这是非预期性non-anticipativity约束的要求——调度决策只能基于预测信息不能在看到实际风电后再去做决定。这一点初学者特别容易搞混如果把 (u) 也定义成 (G \times T \times S)相当于模型获得了“预知未来”的能力算出来的结果在现实中根本没法执行。3.3 最小启停时间约束怎么线性化最小启停时间约束是让很多人头疼的地方。这个约束的本质是如果机组在时段 (t) 启动则它在 (t, t1, ..., tT_{on}^{min}-1) 时段内都必须保持开机如果机组在时段 (t) 停机则它在 (t, t1, ..., tT_{off}^{min}-1) 时段内都必须保持停机。线性化的经典做法是引入启动动作变量 (v_{i,t}^{su} 1) 表示机组 (i) 在时段 (t) 开始启动停机动作变量 (v_{i,t}^{sd} 1) 表示开始停机然后建立状态与动作之间的逻辑关系[ u_{i,t} - u_{i,t-1} v_{i,t}^{su} - v_{i,t}^{sd} ]这条约束把“状态切换”和“动作发生”联系起来开机状态下启动动作取1停机状态下动作取0。然后最小开机时间约束写成[ \sum_{kt}^{tT_{on,i}^{min}-1} u_{i,k} \ge T_{on,i}^{min} \cdot v_{i,t}^{su}, \quad \forall t ]最小停机时间约束同理[ \sum_{kt}^{tT_{off,i}^{min}-1} (1 - u_{i,k}) \ge T_{off,i}^{min} \cdot v_{i,t}^{sd}, \quad \forall t ]这两条约束的直观理解是如果机组在 (t) 时段启动了(v^{su}1)那么未来 (T_{on}^{min}) 个时段里必须全部开机否则左边求和会小于右边约束被违反。注意约束里 (k) 的索引超过调度周期末端时要截断处理否则下标越界或者引入虚假约束这个问题我记得第一次跑的时候就把T弄错了导致边界时段无解。3.4 煤耗成本的分段线性化处理火电机组的煤耗成本通常是出力的二次函数比如 (C_i(P) a_i P^2 b_i P c_i)。YALMIP本身支持二次目标函数CPLEX和Gurobi也都能直接求解MIQP混合整数二次规划理论上可以直接写。但我建议还是做分段线性化原因有三一是MIQP的求解速度比MILP慢不少场景数多时差距显著二是分段线性化在工程中更通用有些求解器不带MIQP功能三是分段线性化的成本曲线对结果精度的影响完全可以控制在0.5%以内没必要为这个牺牲计算性能。分段线性化的做法是把出力范围 ([P^{min}, P^{max}]) 分成若干段每段首尾相连用线性插值逼近原二次函数。Matlab代码实现如下function cost linear_fuel_cost(P, fuel_coeff, seg_points) % fuel_coeff [a, b, c] 二次项、一次项、常数项 % seg_points: 分段断点如 [Pmin, P1, P2, ..., Pmax] a fuel_coeff(1); b fuel_coeff(2); c fuel_coeff(3); n_seg length(seg_points) - 1; % 计算每个断点的真实成本 cost_pts a * seg_points.^2 b * seg_points c; % 生成分段线性插值YALMIP的pwlp函数或者手动实现 % 这里推荐用YALMIP内置的pwlppiecewise linear programming cost pwlp(P, seg_points, cost_pts); endpwlp是YALMIP里处理分段线性函数的内置方法使用起来很简洁它内部会自动引入二进制变量和辅助变量来建模分段区间选择。如果不用YALMIP手写分段线性化也有标准做法把出力决策变量拆成各分段区间上的子变量之和每个子变量限幅落在该分段区间内再引入二进制变量保证相邻分段按顺序填充。这个手写过程容易出错我在这里就不展开了有兴趣的朋友可以找线性规划的建模资料看看。4. 求解策略与测试算例分析4.1 求解器选择与配置经验这套模型最终是混合整数线性规划所以求解器必须支持MILP。我配置的是YALMIP Gurobi的组合如果没装GurobiCPLEX效果也差不多。为什么要用商用求解器说实话Matlab自带的intlinprog也能跑但求解这种多场景多时段的混合整数规划问题中小规模还行规模一大就明显慢而且数值稳定性不够好有时候同一个模型多跑几次结果都可能不一样虽然这更多是容差设置的问题。商用求解器的预求解presolve、割平面加速、启发式分支策略都要成熟很多同一个模型用intlinprog可能跑10分钟还在挣扎Gurobi 30秒内给出接近最优的可行解差距不是一点半点。配置YALMIP Gurobi要注意版本匹配问题特别是Matlab版本和Gurobi的接口版本要对应否则yalmip(clear)之后solvesdp大概率报找不到求解器的错误。我常用的配置方式是% 在Matlab中设置Gurobi路径 addpath(C:\gurobi1100\win64\matlab); gurobi_setup; % Gurobi自带的Matlab接口初始化脚本 % 检查YALMIP是否能找到Gurobi yalmiptest;如果yalmiptest显示Gurobi状态为available那就没问题了。求解时记得设置求解参数我常用的几个参数是options sdpsettings(solver, gurobi, ... verbose, 2, ... % 输出求解日志 gurobi.MIPGap, 0.01, ... % 最优性间隙设为1% gurobi.TimeLimit, 300, ... % 时间限制5分钟 gurobi.Threads, 8); % 并行线程数MIPGap的设置很关键。如果你设成0求解器可能需要非常久才能证明最优性设成1%或2%通常能在几十秒内给出一个非常接近最优的可行解而且目标函数值的差距通常小于0.3%完全够用。4.2 测试系统与参数设置我用的测试系统是经典的10机组系统24个调度时段负荷曲线采用典型夏季日负荷数据峰值负荷在时段19约2100 MW谷值在时段4约1200 MW。风电装机容量设为600 MW占系统峰荷比例约29%这个渗透率在当前很多实际电网中已经相当常见了算出来的结果有参考价值。风光预测曲线我默认给定一个典型形状夜间风大、白天风小、傍晚略有回升——这是我根据实际风电出力的“反调峰”特性设计的恰好和负荷曲线形成鲜明对比。预测误差标准差设为预测值的15%用LHS方法生成1000个原始场景然后削减到20个代表场景。置信度约束那块我暂时用的是软惩罚方式先不收窄场景集合而是让失负荷惩罚系数足够大使模型自动避免失负荷场景。4.3 确定性调度与随机性调度的结果对比为了直观展示随机性建模的价值我做了三组对比实验方案A确定性只考虑风电预测值按一个固定确定性场景建模。方案B随机性风电出力用20个削减场景描述目标函数为期望成本最小。方案C鲁棒性只考虑风电出力波动区间的最坏情况目标函数为最坏场景下成本最小。三组实验的模型参数完全一致只改风电出力的建模方式。结果如下指标方案A确定性方案B随机性方案C鲁棒期望总成本万元126.4118.2137.8失负荷概率4.7%0.3%0%弃风率0.2%0.1%8.5%计算时间秒8.242.615.3方案A的期望成本比方案B高了不少原因在于确定性模型因为“看到”的风电出力偏高安排了较少的常规机组开机导致某些场景下需要调用昂贵的高爬坡率机组来补缺口或者失负荷产生巨额惩罚。方案C虽然彻底消除了失负荷风险但代价是大量弃风——因为方案C认为风功率可能随时跌到最低区间为了保险让常规机组踩满出力结果在实际风功率正常偏高的场景下只能弃风。方案B在这两个极端之间做了平衡它不追求绝对安全而是追求“期望成本最优”这正是随机规划的哲学。再单独看一下方案B的机组出力曲线我发现一个有意思的现象风电预测出力较高的时段比如时段3~5模型仍然保留了一台小容量快速机组待命没有因为它“大概率不需要”就把它停机。这是因为场景集中存在风电偏低的情形保留这台快速机组的额外煤耗成本远小于万一需要它而它没法在短时间内启动造成的失负荷惩罚。这就是随机调度区别于确定性调度的本质特征——它不是用预测值做决策而是用整个概率分布做决策。4.4 场景数量与置信水平对结果的影响最后说一下场景数对求解质量和速度的影响。我扫了一组场景数从5到100的测试场景数期望成本万元求解时间秒相对偏差%5119.83.21.3510118.98.50.5920118.242.60.0850118.1187.40.00100118.1超过10分钟—可以看到场景数从5增加到20期望成本下降了约1.3%但这个差距在工程上也不算特别大从20增加到50成本变化只有0.08%基本可以忽略。考虑到计算时间的爆炸式增长20个场景在这套测试系统下是性价比最好的选择。当然这个结论跟系统规模、风电渗透率、机组数量都有关系具体项目里还是要做一次敏感性扫描来定场景数不能盲目照抄。5. 常见问题与排查实录5.1 求解器报“无可行解”怎么排查这个错误是最常见的我几乎每次改动模型参数都会遇到几次。排查思路我整理成一个标准流程第一先检查松约束。把硬约束改成软约束看看目标函数里的惩罚项是不是有值——如果有值说明问题出在约束过紧与目标函数惩罚的权衡之间如果还是无可行解那就是约束本身逻辑错误。第二检查功率平衡系数和负荷数据量纲。我踩过一次很蠢的坑负荷数据单位是MW风电数据单位是kW二者差了三个数量级平衡约束永远不成立模型当然无解。第三检查最小启停时间约束的初始状态匹配问题。如果系统数据里给定的是机组在调度周期前的初始状态比如开机时长不够导致强制继续开机而你的模型没有正确建模初始状态就可能出现矛盾约束。解决办法是专门为前几个时段补上初始状态约束比如 (u_{i,1} 1) 或 (\sum_{k1}^{T_i^{min}-T_0} u_{i,k} \ge ...)。第四最直接的方法把模型里的二进制变量全部松弛为连续变量在YALMIP里把binvar改成sdpvar如果松弛后仍无可行解那问题出在线性约束本身如果松弛后有解但整数模型无解那问题出在整数的组合约束上大概率是最小启停时间或者分段线性化的辅助变量约束写错了。5.2 场景生成中的正态性假设失效怎么办我之前说过风电预测误差很多时候用正态分布近似是够用的但在极端天气情况下比如台风前、辐合带云系过境误差分布会出现明显的重尾或者双峰特征此时单纯的正态假设会低估极端场景的发生概率导致调度计划在大风极端场景下失衡。如果数据足够我建议用历史误差数据直接做经验分布采样具体操作是把实测误差按分位数排布然后用datasample加权重直接抽样绕开分布拟合这一步。如果数据不足退而求其次用t分布自由度5~8替代正态分布t分布的尾巴更厚对极端事件的覆盖更保守。我做敏感性分析时还会额外做一组“压力测试”人为把某个风电出力极低的场景概率权重从1%调到10%观察目标函数和开机计划是否发生显著改变如果变化太大说明模型的鲁棒性不足这时候要稍微调大旋转备用约束或者调低置信度。5.3 求解速度太慢的优化套路最后送几个让MILP求解速度提升的经验这些都是我在实际调参过程中验证过的一是初始化可行解warm start。先用确定性场景预测值快速求解得到一组机组开停机计划然后把这组 (u) 作为随机模型的初始解传给求解器。商用求解器非常擅长在好的人工解基础上做改进这往往能把求解时间压缩一半以上。YALMIP里可以通过assign和sdpsettings(gurobi.StartNodeLimit, ...)来设置初始解。二是善用对称性破缺symmetry breaking。10机组测试系统里有两三台类型相同的机组场景模型对相同机组做对称的调度决策会让分支定界过程大量重复搜索等价节点。通过添加机组优先顺序约束比如 (P_1 \ge P_2 \ge P_3)对于同型机组可以显著收紧搜索空间。三是调整求解器的容差参数。MIPGap从默认的万分之一放宽到百分之一求解时间往往能降低一个数量级而经济效益损失不到千分之三。在工程决策场景中这个交换是非常划算的。记住这个原则先跑一版严格最优解确认模型正确后续批量做敏感性分析时就放宽容差加速效率和安全兼顾。四是考虑用Benders分解或拉格朗日松弛做大规模场景加速但这就是另一个层面的算法优化了。对于大多数教学场景和中型算例单机MILP已经足够不用一开始就上分解算法避免把简单问题复杂化。6. 这套模型的扩展方向与实际应用价值最后再说几句扩展方向。这套随机动态经济调度模型框架本身是“骨架”风电只是第一个实例化的不确定源。分布式光伏、电动汽车负荷、储能系统调度它们的随机性和波动性本质上是同一类问题——无非是改掉不确定变量的概率分布和参数边界再把储能充放电约束加进去。我在实际项目里就把这套框架扩展成了风光储联合调度模型做法是在目标函数里新增储能充放电的损耗成本和充放电功率约束风电和光伏各生成一组场景用联合场景矩阵代入模型求解。扩展过程很顺滑这得益于底层YALMIP建模的高可读性和模块化代码结构。另一个值得尝试的方向是反馈调度rolling horizon。日前调度做了24小时计划实际执行时风电预测每4小时刷新一次与其固定执行已过时的日前计划不如每4小时重跑一次滚动调度只执行未来8~12小时的部分计划。这种模型预测控制MPC式的调度方式和动态经济调度在数学框架上是完全兼容的只多了一步“每次求解后只取前几个时段的解”这也是从学术模型走向工程落地的关键一步。我个人在写完这套模型后最大的感触是随机调度的重点不在算法本身有多复杂而在于你是否准确理解了“决策必须在不确定性被揭示之前做出”这一非预期性约束。很多跑出错误结果的项目翻来覆去检查也找不到问题最后发现就是某个决策变量被错误地定义成了场景相关导致模型拥有了现实世界不可能有的“预知能力”。希望这篇文章能帮大家避开这个最大的坑也欢迎有同行在类似问题上一起来交流心得。
