做含分布式电源的配电网日前两阶段优化调度绕不开三个坎不确定性怎么建模、两阶段结构怎么设计、Matlab代码怎么落地。我最近完整跑通了一个典型方案——以IEEE 33节点配电网为算例用YalmipGurobi实现“日前决策实时调整”的两阶段鲁棒优化调度光伏、风电、储能、微型燃气轮机全部纳入从数学建模到迭代求解再到结果分析走了一遍。这篇就把整个模型的思路拆解、关键约束形式、CCG求解流程以及我实际调试中踩过的坑完整记录下来供正在做配电网优化调度课题的研究生、需要算法对比写论文的同学以及做园区微网或台区能量管理的工程师参考。1. 项目背景与问题本质为什么配电网调度要“两阶段”1.1 分布式电源并网后调度前提彻底变了传统配电网的调度逻辑很单纯负荷是预测出来的电源在变电站侧功率从变电站单向流向用户调度中心只需要按负荷曲线倒推变电站出力再把无功补偿和变压器分接头调一调就行。分布式电源大规模接入之后这个前提被打破了。光伏和风电的出力由天气决定既不连续也不完全可预测正午阳光好时光伏大发负荷却不一定处在高峰夜间风电出力猛增负荷反而在低谷。配电网从一个“被动接受功率的终端网络”变成了“功率可能双向流动的有源网络”。这就带来三个直接后果第一DG出力超过就地负荷时潮流反向变电站母线以下的线路可能反向过载第二分布式电源接入点附近的节点电压会被抬高极端情况下末端电压越上限第三调度员手里可以调的开关和机组变多了但不确定性也更大了用单一确定性的负荷和光伏预测值去排计划到实际运行时往往完全跑不动。我在算例里特意构造了一个高渗透率场景33节点系统在末端接入2组光伏和1组风电光伏渗透率约30%。中午12点附近光伏出力接近1.2 MW而5号节点到18号节点这一带的负荷总和只有0.9 MW末端节点电压直接从0.98 pu被抬到1.06 pu如果不做优化调度而是按传统方式直接并网电压越限是必然的。1.2 “日前”和“两阶段”到底在解决什么先看“日前”。日前调度指的是提前24小时、以1小时或15分钟为步长典型是24个时段或96个时段制定第二天的运行计划。它的本质是“预先排班”哪些机组明天要开机、储能什么时候充什么时候放、可中断负荷切哪几户、联络开关要不要倒闸这些都需要提前定下来。再看“两阶段”。为什么一天之内的调度还要分两个阶段因为日前做计划时光伏和风电的真实出力还不知道只能依据预测值。预测值和实际值之间的偏差就是不确定性的来源。两阶段结构把决策分成两类第一阶段now-here决策在不确定量实现之前就要拍板的事比如机组启停、储能充放电状态、联络开关状态。这类决策涉及0-1整数变量一旦定了就不能随便改改起来代价很大。第二阶段wait-and-see决策等在某个具体场景下看到实际DG出力后再做功率调整比如常规机组出力、储能实际功率、向上级电网购电功率、必要时弃光弃风甚至切负荷。这类决策是连续的可以随场景变化。对配电网来说两阶段的合理性用一句话就能说清日前计划定的是骨架日内调整补的是血肉。如果只有一个阶段那就只能赌预测完全准确赌错了要么电压越限要么失负荷。如果每个时段都做多场景随机优化计算量又爆炸。两阶段结构是计算精度和求解复杂度之间的一个平衡点。1.3 与输电网经济调度的本质差异很多初次做配电网优化的人会想当然地把输电网经济调度的模型搬过来改一改。我在项目早期也这么干过结果发现水土不服。核心差异有三个对比项输电网经济调度配电网两阶段优化调度网络结构环网为主潮流方向清晰辐射状为主潮流可能双向潮流模型直流潮流或交流潮流R/X小P-Q解耦可用R/X大P-Q耦合强需用DistFlow前推回代模型电压问题主要通过无功调节有功与电压弱相关DG有功出力直接抬升/拉低节点电压有功电压强耦合不确定性负荷预测误差为主DG出力强波动 负荷双重不确定性决策对象大机组开停机、联络线功率小型MT启停、储能SOC、联络开关、弃光弃风、可中断负荷一句话总结输电网调度主要管“发多少电、谁发”配电网调度要管“功率从哪来、往哪去、电压在哪儿越限”问题是强耦合、强不确定、强约束的。这也是为什么必须用两阶段鲁棒/随机优化框架来处理。2. 数学模型拆解目标函数与约束条件的工程解读2.1 第一阶段模型定骨架的“上层决策”第一阶段的目标函数我在代码里写成了这样常规机组微型燃气轮机的启停费用和空载费用储能充放电切换的惩罚项避免频繁切换状态联络开关动作惩罚配电网开关设备寿命有限倒闸次数不能太多上级电网购电的日前计划成本按预测场景估算。决策变量是MT机组的开停机状态 u_t ∈ {0,1}MT机组的启动/停机动作变量 v_t, w_t ∈ {0,1}储能充电/放电状态变量如果需要建模状态约束联络开关状态当两阶段中嵌套了网络重构时。第一阶段的约束主要包括最小开停机时间约束机组不能频繁启停功率平衡的粗略形式用预测值校验备用容量约束保证最坏情况下系统也能保住负荷。这里有一个容易被忽略的细节第一阶段其实不需要把全部潮流约束都铺上去。因为第一阶段用的是预测场景你在这个阶段做精细潮流计算没有意义反而会让MIP问题规模膨胀。正确的做法是第一阶段的计划只保证“基础可行”把精确的潮流校验、电压越限判断全部留给第二阶段子问题去做。这也是两阶段分解算法的一个基本原则。2.2 第二阶段模型应对具体场景的“下层调整”第二阶段是在第一阶段计划 x 已确定的情况下对每个实际可能出现的DG出力场景 ξ 做再调度目标是最小化在该场景下的调整总成本。如果采用鲁棒优化表述完整的目标是$$ \min_{x} \left( c_1^T x \max_{\xi \in U} \min_{y \in F(x, \xi)} c_2^T y \right) $$其中 x 是第一阶段的启停/状态决策U 是光伏和风电出力的不确定集合y 是第二阶段的连续运行决策包括MT机组实际有功/无功出力储能充放电功率向上级电网的购电功率弃光、弃风功率切负荷功率DG的无功出力如果DG具备无功调节能力。第二阶段约束必须包含完整的DistFlow潮流方程、节点电压上下限、支路电流容量约束、机组爬坡约束、储能SOC连续性约束、弃风和弃光功率上下限等。注意切负荷费用和弃风弃光惩罚系数一般都设置得很高这样优化模型会尽量避免这些不友好的运行方式只在确实无法消纳时才动用。2.3 分布式电源出力模型与不确定集合工程上怎么处理“说不准”光伏和风电的出力随机性建模方式大致有三条路线场景法用蒙特卡洛抽样生成大量出力场景再做场景削减典型用同步回代法削减到10~20个场景适合随机规划不确定集法不枚举具体场景而是用一个集合描述出力可能波动的范围适合鲁棒优化机会约束用概率约束描述“电压越限概率不超过5%”需要对约束做确定性转化。我在这个项目里用的是鲁棒优化路线这也是“两阶段”最常搭配的路线。不确定集采用“盒式预算约束”的标准形式$$ U \left{ P_i^{DG} \in \left[ P_i^{pre} - \Delta P_i, P_i^{pre} \Delta P_i \right], \quad \sum_i \frac{|P_i^{DG} - P_i^{pre}|}{\Delta P_i} \leq \Gamma \right} $$这里面有两个关键参数ΔP_i 是单个DG出力的最大偏差由历史预测误差统计得到一般取预测值的15%~25%Γ 是“预算系数”它的物理含义是同时偏离预测最坏值的DG数量上限。Γ0时退化为确定性模型Γ越大鲁棒水平越高但运行成本也越保守。Γ 的选取直接影响调度结果的激进程度。实际做仿真时我建议做一组Γ从0到最大值的灵敏度曲线一方面是论文需要另一方面能直观看到“花多少钱买多少鲁棒性”这种权衡曲线放到论文里非常有说服力。光伏出力本身我更推荐用Beta分布描述风电用Weibull分布。如果图省事直接用正态分布抽样很容易抽到负的出力值这在物理上是荒谬的还会让后面的功率平衡约束直接不可行——我早期就踩过这个坑后面会专门讲。2.4 DistFlow支路潮流模型配电网优化不可回避的“正确打开方式”配电网R/X比大输电网那套P-Q解耦的直流潮流在这里误差很大。配电网优化调度里最经典的潮流模型是Baran和Wu提出的DistFlow支路潮流模型。对任一条支路 i-ji为父节点j为子节点DistFlow写成$$ P_{ij} - \sum_{k:j\to k} P_{jk} - r_{ij} l_{ij} -P_j^{inj} $$$$ Q_{ij} - \sum_{k:j\to k} Q_{jk} - x_{ij} l_{ij} -Q_j^{inj} $$$$ V_j V_i - 2(r_{ij} P_{ij} x_{ij} Q_{ij}) (r_{ij}^2 x_{ij}^2) l_{ij} $$其中 l_{ij} |I_{ij}|^2即支路电流平方。这里出现了非线性项直接放进优化模型会变成非凸问题求解非常困难。工程上的标准做法是“二阶锥松弛”引入中间变量把非线性约束等价转化为一个凸锥约束$$ \left| \begin{bmatrix} 2P_{ij} \ 2Q_{ij} \ l_{ij} - V_i \end{bmatrix} \right| \leq l_{ij} V_i $$这个松弛在辐射状配电网加上合适的凸目标函数下是精确的也就是说松弛解就是原问题的可行解。这也是为什么现在配电网优化调度几乎都走SOCP路线——模型既能被求解器高效处理结果又可靠。3. 求解策略与Matlab实现架构3.1 为什么要用CCG而不是硬解一个单层大模型把前面两层模型直接“摊平”成一个大规模混合整数二阶锥规划MISOCP丢给求解器理论上是可以的。但实际跑起来你会发现两个问题规模爆炸。33节点、96时段、3台DG、1台MT、1套储能0-1变量和二阶锥约束叠加在一起模型维度非常可观直接求解经常几小时出不来还要占用大量内存鲁棒优化里那个 max-min 嵌套结构没法直接交给求解器必须手动把子问题的内层最小化转换成对偶最大化再交给求解器处理。所以实际工程中都采用分解算法。最常用的两个是Benders分解和列与约束生成CCG。两者的区别我直接说结论Benders分解通过割平面逐步逼近第二阶段的代价函数次梯度收敛在某些场景下迭代次数多、收敛慢CCG的思路更直接——它不逼近代价函数而是直接把第二阶段变量的“副本”和对应场景约束添加到主问题里让主问题逐步“长出”第二阶段的结构。对线性问题和混合整数问题CCG的迭代次数明显少于Benders而且越到后期优势越大。我在同一套算例上对比过CCG大约迭代6~9次达到0.1%的gapBenders则需要20次以上而且后期收敛明显变慢。所以这个项目里我最终选的是CCG。3.2 CCG迭代流程四步循环别搞乱顺序CCG的完整流程是这样的初始化给定第一阶段决策的初始可行解 x0可以用确定性模型快速算一个设置下界 LB-∞上界 UB∞迭代次数 k0收敛精度 ε0.001。求解主问题MP在已有场景集合中最小化第一阶段成本已有场景下的第二阶段成本得到 x*更新 LB 当前MP最优值。固定 x求解子问题SP*子问题对不确定量求 max对第二阶段变量求 min即寻找当前最优第一阶段决策下“最恶劣的场景”和对应的最小调整成本。子问题求解结果得到最恶劣场景 ξ* 和目标值 R(x*)。更新 UB min(UB, min_x c1 x R(x*))。收敛判断如果 UB - LB ≤ ε停止否则把 ξ* 对应的第二阶段变量副本和约束加入主问题kk1回到步骤2。这里有个关键点子问题内部是 max-min 嵌套实际求解时要把内层 min 线性规划问题写成 KKT 条件或者对偶形式把整个子问题变形成一个单层 max 问题。对线性子问题直接取对偶然后做一个 max 化的线性规划就行对二阶锥子问题需要用到二阶锥对偶理论Yalmip可以自动处理一部分但理解原理对调试很重要。3.3 Matlab工程文件架构代码怎么组织才不变成“一坨”我强烈建议工程化组织代码不要一个main脚本从头写到尾。我这个项目的文件结构是这样的project/ ├── main_两阶段调度.m # 主入口定义系统参数调CCG循环 ├── case33_parameters.m # IEEE 33节点所有数据和负荷/DG曲线 ├── build_MP.m # 构建主问题Yalmip模型 ├── build_SP.m # 构建子问题Yalmip模型 ├── solve_subproblem_dual.m # 子问题对偶化处理 ├── plot_results.m # 输出调度曲线、电压剖面、收敛曲线 └── data/ ├── load_curve.csv # 典型日负荷曲线96点 ├── pv_curve.csv # 光伏出力预测曲线 ├── wind_curve.csv # 风电出力预测曲线 └── scenarios.mat # 削减后的场景集主程序只负责“数据读取→初始化→CCG迭代→结果保存”模型构建和求解全部封装在函数里。这样后期换算例比如从33节点换到123节点、换求解器、换不确定性参数都只需要改对应的函数不会牵一发动全身。求解器方面我的选择是Yalmip Gurobi。Yalmip负责建模翻译Gurobi负责真正求解MISOCP。为什么不用Matlab内置求解器因为内置的 linprog/quadprog 不支持二阶锥约束intlinprog不支持锥约束而这个问题本质上是混合整数二阶锥规划必须要Gurobi或CPLEX这类商业求解器。学生可以用免费的学术license也可以用SCS等开源求解器顶替不过混合整数部分开源求解器效率会差不少。4. 核心代码实现细节与关键操作要点4.1 数据准备IEEE 33节点算例与负荷/DG曲线IEEE 33节点系统是配电网优化调度的“标配”算例母线上手必须要会用。它的标准参数我整理了一份速查参数数值基准电压12.66 kV节点数33支路数37含5条联络开关支路系统总负荷约3.72 MW 2.30 Mvar根节点电压1.00 pu变电站母线电压允许范围0.95~1.05 pu我按国标习惯取的原始算例的数据在网上很好找33个节点的支路阻抗、节点负荷表都是公开的。拿到之后要做的第一件事是画一下拓扑图把DG接入位置标注上去确认根节点的位置和支路方向——DistFlow对支路方向很敏感方向搞反了潮流约束就全错了。负荷曲线和DG出力曲线我建议统一用“有名值”而不是标幺值。理由很简单后续要对照论文结果、要验证潮流结果时有名值kW、kVar、kV一目了然标幺值虽然方便计算但错一个基准值就全线崩溃调试成本太高。光伏预测曲线可以用典型晴天曲线加一个正弦包络近似风电曲线则是“夜间大、白天小”的典型反调峰形状。这些曲线不用太精细但一定要保证在某个时段构造出“光伏大发负荷低估”的极端情况否则鲁棒优化最恶劣场景的求解结果会显得平淡看不出两阶段模型的优势。4.2 Yalmip建模关键代码示例我在build_MP.m里用Yalmip定义变量的方式如下% 定义第一阶段决策变量以小时为周期24时段 u binvar(ng, T, full); % MT机组开停机状态 v binvar(ng, T, full); % 启动动作 w binvar(ng, T, full); % 停机动作 s_ch binvar(n_ess, T, full); % 储能充电状态 s_dis binvar(n_ess, T, full); % 储能放电状态 % 第二阶段场景相关变量以场景s为例 p_dg sdpvar(ndg, T, full); % DG有功出力 p_mt sdpvar(ng, T, full); % MT机组出力 p_ess sdpvar(n_ess, T, full); % 储能功率正为放电 soc sdpvar(n_ess, T1, full); % 储能SOC p_curtail sdpvar(ndg, T, full); % 弃风弃光功率 p_cutload sdpvar(n_load, T, full); % 切负荷功率DistFlow锥约束我直接用Yalmip的norm命令写成constraints [constraints, norm([2*p_ij; 2*q_ij; I2_ij - V_j]) I2_ij V_j];这是二阶锥约束的标准写法Yalmip会识别成锥约束并传给Gurobi。注意 I2_ij 是电流平方的连续变量V_i、V_j 是节点电压平方的连续变量DistFlow里电压用平方形式。别直接把 I_ij 当电流变量写会把约束变成非凸的。4.3 CCG迭代的收敛判据与热启动技巧CCG迭代的收敛判据我写在主函数里UB inf; LB -inf; gap 1; k 0; while gap 1e-3 k 30 % 求解主问题 [x_opt, LB_new] solve_MP(); LB max(LB, LB_new); % 固定x求解子问题 [worst_scenario, R_value] solve_SP(x_opt); UB min(UB, c1*x_opt R_value); gap abs(UB - LB) / abs(UB); % 把最恶劣场景加入主问题 MP add_scenario_to_MP(MP, worst_scenario); k k 1; end这里有一个实操细节CCG的主问题每迭代一轮就要重新求解一次模型规模会逐渐增大。如果每次都用冷启动求解时间会指数级上升。我用的是Gurobi的warm start机制——把上一轮的x_opt作为初始解传给下一次求解很多约束和变量的解在相邻两轮之间变化不大热启动能省掉大量分支定界时间。Yalmip里设置方法如下ops sdpsettings(solver, gurobi, verbose, 2, gurobi.WarmStart, 1);实测下来带热启动的CCG在33节点系统上单轮MP求解时间从最初的十几秒降到后来的几秒整个循环加起来不到两分钟这个效率写论文做批量仿真完全够用。4.4 结果输出画图要画哪几张才“出活”调度结果出来之后绘图部分我建议至少输出三张图这也是论文里标配调度曲线图横轴是24小时或96时段纵轴是功率画出DG出力、MT出力、储能充放电、上级电网购电、弃风弃光功率能直观看出调度策略。电压剖面图选几个关键时段光伏高峰、负荷晚高峰画出33个节点的电压分布电压曲线不能有越限这是两阶段模型有效性的直接证据。迭代收敛曲线图画出UB和LB随迭代次数的变化两条线逐渐收敛到同一个值。审稿人看到这张图就会认可你的算法实现是可靠的。5. 常见问题与排查技巧实录5.1 求解报错速查表我把这个项目里遇到过的典型问题和排查方法整理成一张速查表照着查可以省很多时间报错/现象可能原因处理办法求解结果返回Infeasible problem功率平衡约束与DG/负荷数据冲突储能SOC初值不闭合用diagnostics optimize(...)查看diagnostics.problem具体错误先删掉部分约束做二分定位结果全是NaNYalmip变量维度不匹配某变量未在约束中被限定检查size()打印每个变量的维度检查约束拼接是否有[]维度冲突锥约束报错“nonconvex”笔误把I2_ij写成I_ij^2导致约束变成二次非凸将电流平方设为独立变量I2统一用norm写锥约束子问题无界第二阶段变量切负荷、弃光没有加上下限或上下限太宽松给所有第二阶段变量手动加上物理边界比如弃光功率 ∈ [0, P_dg_max]CCG迭代不收敛主问题在下一次迭代时没有正确保留历史场景约束检查add_scenario_to_MP是否把新场景的新变量副本完整加入而不是覆盖了之前的约束求解极慢1小时0-1变量过多MP没有热启动先用热启动考虑将联络开关变量固定为常开先跑通基本流程再做扩展5.2 我踩过的几个“教科书不写”的坑第一个坑是光伏出力抽样出现负值。用正态分布拟合光伏出力误差抽出来一堆负数功率平衡直接崩掉。后来查资料才发现光伏出力更适合Beta分布如果你一定要用正态分布也要做截断处理——小于0的样本直接置0。第二个坑是电压基准选错。我一开始图方便把整个系统电压标幺化结果忘记根节点电压基准是12.66kV而不是0.4kV电压越限检查形同虚设。后来全部改成有名值用V单位kV写死约束V_min0.9512.66V_max1.0512.66问题立刻清晰。第三个坑是储能SOC约束的索引错误。我写第t1时段的SOC等于第t时段SOC加充电功率减去放电功率结果在96时段循环里有一个地方把t写成了t1导致SOC“凭空增加”储能变成永动机优化结果里储能一个劲放电赚钱。后来对着能量守恒检查才发现。这里我建议用向量化写法避免循环里的索引错误SOC(:, 2:end) SOC(:, 1:end-1) eta_ch * p_ch - p_dis / eta_dis;5.3 怎么判断调度结果“物理上是对的”求解器说收敛了不代表结果就是对的。我每次跑完都会做三层合理性检查功率是否守恒所有时段的发电购电 负荷网损充电误差超过1%一定有问题电压是否在限值内任意节点的电压幅值曲线不能越0.95~1.05pu如果某个节点电压贴着上限走说明DG出力确实把它抬起来了这是符合物理预期的策略是否符合常识光伏大发的正午时段如果结果依然是“拼命购电而不充电、不弃光”那模型大概率有问题合理的策略应该是正午储能充电、光伏满足本地负荷傍晚光伏衰退后储能放电、补齐晚高峰缺口。我举个典型结果供参考在Γ8的鲁棒配置下与确定性模型相比购电成本上升约6.4%但最恶劣场景下的电压越限概率从18%降到0弃光率从12%降到2%以内。这就是“花一点经济代价换供电安全”的典型权衡。6. 扩展思路从33节点模型走向更实用的形态这个项目做完之后我还有几个方向想继续尝试也给正在做类似课题的朋友一些参考。第一个扩展方向是把单目标扩展为多目标。当前模型只做了经济成本最小但配电网调度里网损、开关动作次数、DG利用率、碳排放都是重要指标。把碳排放或网损加进目标函数就成了多目标优化问题可以用加权法或者epsilon约束法处理得到的Pareto前沿在论文里很有说服力。第二个方向是把日前两阶段改成“日前日内滚动”的模型预测控制结构。日前做一整天的骨架计划日内每隔15分钟拿最新预测数据滚动修正一次这样能更精细地应对光伏/风电的短时波动。代价是计算量更大但工程落地时更有价值。第三个方向是考虑网络重构和两阶段调度的联合优化。把联络开关的开合状态也变成第一阶段的0-1变量问题会更有挑战性但网络结构的灵活性会给DG消纳带来额外空间。不过要注意配电网重构会让DistFlow的父-子节点关系在运行时改变建模复杂度会上一个台阶建议先把固定拓扑调顺了再扩展。真要说这个项目带给我最大的体会反而不是算法本身而是“建模思路比代码重要”这件事。两阶段鲁棒优化的代码实现其实就那么几步拆层、对偶、迭代、加场景。但每个环节的选择——不确定集形式、预算系数Γ、锥松弛写法——背后都有清晰的物理逻辑。把逻辑想清楚了代码只是一个翻译过程逻辑没想清楚就急着跑仿真大概率是在一堆报错和NaN里反复打转。如果你也在做类似的配电网调度课题建议按这个路线走先把确定性单阶段模型跑通验证潮流和成本结果符合物理直觉再加入不确定集合和子问题跑CCG迭代最后再调Γ、换场景、加扩展目标。每一步都踩实了再往前走比一次到位顺畅得多。
