简介面向电力系统调度与优化研究的一份MATLAB资源包聚焦安全约束机组组合SCUC问题同时覆盖基于交流潮流方程与直流潮流方程两种建模方式契合电力系统分析中对精度与计算效率的双重需求。资源包共有9个文件其中7个m脚本构成核心代码涵盖系统模型定义、发电机出力与线路热容量约束、遗传算法/粒子群等优化求解、负荷数据预处理及结果可视化等完整流程另附1个txt说明和1个zip补充包整体仅263KB便于快速部署。该项目已有86人学习浏览适合电力工程、控制理论或计算科学方向的学生和研究人员也适合希望将人工智能技术引入电力调度优化的开发者。包内代码可帮助理解交流/直流潮流方程在SCUC中的建模差异掌握优化算法设计与MATLAB实现并可探索机器学习预测负荷、改进机组启停决策的AI融合路径对经济调度研究与智能算法实践有直接参考价值。1. 安全约束机组承诺模型交流潮流与直流潮流两种建模范式的分水岭拿到《电力系统安全约束单位承诺模型包括基于交流潮流方程和直流潮流方程的.zip》这个压缩包先别急着解压。标题里的“单位承诺”是英文 Unit Commitment 的直译国内工程与学术界通常叫“机组组合”。机组组合解决的是未来 24 小时里每台机组何时启动、何时停机、各发多少电加上“安全约束”之后问题从单纯的经济调度升级为在机组运行约束的基础上计及输电网给潮流限制保证调度结果送得出去、线路不越限。而这个标题最值得拆解的点是把交流潮流方程和直流潮流方程放在同一个模型里对照——几乎所有做电力系统优化的人都会在这里遇到分岔路口交流模型精确但非线性非凸直流模型线性可解但牺牲了电压和无功信息。这篇文章把这两条路的理论、建模代码、求解技巧和验证方法讲透适合刚接触 SCUC 的研究生、做电网调度的工程师以及想把这类开源包改造到自己项目里的算法开发人员。2. 潮流方程如何进入优化模型从交流完整形式到直流三步近似2.1 交流潮流方程的极坐标形式与其数学本质交流潮流方程描述的是节点注入功率与节点电压幅值、相角之间的非线性关系。极坐标形式下节点 i 的有功注入方程写成P_i V_i Σ V_j ( G_ij cosδ_ij B_ij sinδ_ij )无功注入方程对应Q_i V_i Σ V_j ( G_ij sinδ_ij - B_ij cosδ_ij )其中 V_i 是节点电压幅值δ_ij δ_i - δ_j 是节点间相角差G_ij 和 B_ij 是导纳矩阵的实部与虚部。这股非线性来自两处V_i 与 V_j 的乘积项以及三角函数项。当把这样的约束嵌入机组组合模型时变量集从有功功率 P 扩展为 P、Q、V、δ 四个维度机组的无功出力和节点电压上限也一并进入约束集合问题性质变成了混合整数非线性规划MINLP。这个“变身”带来的求解代价是数量级的提升。机组组合本身已包含 0-1 整数变量处理启停状态交流潮流把约束变成非凸可行域整数变量和非线性约束交织在一起全局最优解基本无望工程上只能接受局部最优或近似解。所以当你看到某个 SCUC 项目包同时给了交流和直流两个版本通常的用意是直流版跑日前计划交流版做安全校核或作为精确对照。2.2 直流潮流方程的三步近似从三角函数到线性方程直流潮流DC Power Flow不是“直流电”的潮流而是交流潮流的线性化近似。从交流支路有功潮流公式出发忽略支路电阻和并联电容支路导纳只剩电抗分量于是 B_ij -1/x_ijG_ij 0支路有功潮流简化为P_ij V_i V_j sinδ_ij / x_ij接下来做三步近似。第一步假设节点电压幅值都在额定值附近取 V_i V_j 1.0 pu消去电压乘积项。第二步假设线路两端相角差很小sinδ_ij ≈ δ_ij。第三步忽略支路电阻和所有并联支路。最终得到P_ij (δ_i - δ_j) / x_ij线路有功潮流只与两端相角差和线路电抗有关是线性关系。节点层面的功率平衡方程为P_i Σ_j (δ_i - δ_j) / x_ij写成矩阵形式就是经典的 P Bδ其中 B 是节点电纳矩阵对角线元素是连接该节点所有支路电抗倒数之和非对角线元素是 -1/x_ij。需要指定一个平衡节点作为角度参考通常取其相角为 0。这套近似在多数输电网规划与调度中精度可接受。对以有功潮流为主的机组组合问题直流潮流把非线性约束全部线性化模型退化为混合整数线性规划MILP商业求解器可以直接求出全局最优。代价是丢失了电压幅值、无功功率和网损信息当系统运行在重负荷或电压敏感区域时直流模型的调度结果可能与真实交流潮流有较大出入。2.3 两种建模范式放入 SCUC 后的差异对照对比维度交流潮流 SCUC直流潮流 SCUC变量集合有功、无功、电压幅值、相角有功、相角约束形态非线性非凸等式约束线性等式约束问题类型MINLPMILP求解难度极高一般只能得到局部最优较低可求得全局最优典型求解器IPOPT、KNITRO、CONOPTCPLEX、Gurobi、CBC、GLPK是否含电压约束是可校核电压越限否电压信息被屏蔽网损处理隐含在潮流方程中忽略产生较小误差选型逻辑很简单若要研究网络阻塞对日前市场出清的影响直流模型已经够用速度快且结果稳定若要研究电压支撑、无功优化或高新能源渗透下的电压越限问题必须上交流模型。很多生产系统采用“直流优化 交流校核”的两级架构而不是直接求交流 SCUC。这个思路会在后面代码实战与验证章节中具体展开。3. 用 Pyomo 搭建可复现的 SCUC 模型从直流版本到交流版本3.1 数据集设计三节点小系统与机组参数表教学和验证场景不需要 IEEE 300 节点那么大的规模一个能看出“安全约束起作用”的三节点系统足够。节点 1 接便宜大机组节点 2 接贵一些的中型机组节点 3 接负荷节点 1 到节点 3 的线路容量设置偏小。这样在没有网络约束时便宜的机组 1 会尽量多发电导致线路 1-3 过载加上直流潮流约束后机组 2 必须多出力才能把线路潮流压回限额以内。负荷按 24 时段设置三段式曲线夜间 150MW白天高峰 220MW晚高峰 180MW。三台机组参数如下机组节点Pmin/MWPmax/MW固定成本/元边际成本/元每MWh最小开机/h最小停机/h爬坡/MW每hG1110300400203290G2210200500352260G325100300501140线路参数1-2 电抗 0.25 pu容量 150MW1-3 电抗 0.20 pu容量 120MW2-3 电抗 0.3 pu容量 130MW。基准容量取 100MVA机组出力和负荷的标幺值由程序折算。这些参数直接写进 Python 字典便于后续扩展成读取 Excel 或 CSV 的版本。T 24 lines [ # (起点, 终点, 电抗x, 容量上限MW) (1, 2, 0.25, 150), (1, 3, 0.20, 120), (2, 3, 0.30, 130), ] # 机组用两级成本函数近似固定成本 线性边际成本 units { 1: {bus: 1, pmin: 10, pmax: 300, fixed: 400, slope: 20, UT: 3, DT: 2, RU: 90}, 2: {bus: 2, pmin: 10, pmax: 200, fixed: 500, slope: 35, UT: 2, DT: 2, RU: 60}, 3: {bus: 2, pmin: 5, pmax: 100, fixed: 300, slope: 50, UT: 1, DT: 1, RU: 40}, } demand { 3: [150, 150, 150, 160, 170, 180, 190, 200, 210, 220, 220, 210, 200, 190, 180, 180, 190, 200, 210, 220, 200, 180, 160, 150] }说明固定成本对应机组空载热备用成本边际成本对应燃煤或燃气可变成本。负荷放在节点 3节点 1 和 2 无本地负荷。线路潮流上限单位为 MW程序内部使用标幺值时需要统一除以基准功率。3.2 直流 SCUC 建模整数变量、爬坡约束与线性潮流约束下面这段 Pyomo 代码实现了直流 SCUC 的主体框架。集合 G 和 T 分别表示机组和时段变量 u、start、stop 是 0-1 变量p 是有功出力连续变量theta 是节点相角。目标函数是所有时段的固定成本和边际成本之和。import pyomo.environ as pyo model pyo.ConcreteModel() model.G pyo.Set(initializeunits.keys()) model.T pyo.Set(initializerange(1, T 1)) model.B pyo.Set(initialize[1, 2, 3]) model.L pyo.Set(initialize[(i, j) for (i, j, x, cap) in lines]) model.u pyo.Var(model.G, model.T, domainpyo.Binary) model.start pyo.Var(model.G, model.T, domainpyo.Binary) model.stop pyo.Var(model.G, model.T, domainpyo.Binary) model.p pyo.Var(model.G, model.T, domainpyo.NonNegativeReals) model.theta pyo.Var(model.B, model.T, domainpyo.Reals) def objective_rule(model): expr 0 for g in model.G: for t in model.T: expr units[g][fixed] * model.u[g, t] expr units[g][slope] * model.p[g, t] return expr model.objective pyo.Objective(ruleobjective_rule, sensepyo.minimize)目标函数里用固定成本乘以机组状态变量 u边际成本乘以有功出力 p。当固定成本被设为 0 时模型退化为纯经济负荷分配但机组频繁启停会频繁发生所以实际中不能省略固定成本项。单台机组在一个时段内的启停损耗若要更精确可以再加启动成本项这里用 start 变量乘以单位启动费用即可扩展到目标函数中。出力上下限和启停状态绑定def p_max_rule(model, g, t): return model.p[g, t] units[g][pmax] * model.u[g, t] model.pmax_con pyo.Constraint(model.G, model.T, rulep_max_rule) def p_min_rule(model, g, t): return model.p[g, t] units[g][pmin] * model.u[g, t] model.pmin_con pyo.Constraint(model.G, model.T, rulep_min_rule)启停机逻辑约束 start[t] - stop[t] u[t] - u[t-1]保证“启动”和“停机”两个二进制变量与原状态变量的切换一致。首时段用初始状态替代这里假设所有机组从零状态开始即 u[g, 0] 0。def u_transition_rule(model, g, t): if t 1: return model.start[g, t] - model.stop[g, t] model.u[g, t] return (model.start[g, t] - model.stop[g, t] model.u[g, t] - model.u[g, t - 1]) model.transition pyo.Constraint(model.G, model.T, ruleu_transition_rule)最小开机与最小停机的处理方式是 SCUC 建模中最容易写错的部分。这里采用标准展开式写法若机组在时段 t 启动则从 t 时刻起连续 UT 个时段都必须处于开机状态若在 t 时段停机则从 t 时刻起连续 DT 个时段都必须处于停机状态。def min_uptime_rule(model, g, t, k): if t k T: return pyo.Constraint.Skip return model.u[g, t k] model.start[g, t] def min_downtime_rule(model, g, t, k): if t k T: return pyo.Constraint.Skip return model.u[g, t k] 1 - model.stop[g, t] max_ut max(units[g][UT] for g in model.G) max_dt max(units[g][DT] for g in model.G) model.mut pyo.Constraint(model.G, model.T, range(max_ut 1), rulemin_uptime_rule) model.mdt pyo.Constraint(model.G, model.T, range(max_dt 1), rulemin_downtime_rule)最小开机/停机约束中的索引 k 从 0 到最大时长但只有在 t k 不超过 T 时限才生效。这种写法直观但会产生较多冗余约束在面对大规模系统时工程技术上通常改成“启动后累计状态”的动态表达减少约束数量。本案例规模小冗余可以忽略。爬坡约束需要考虑启停瞬间的特殊处理。若机组在 t-1 时段停机、t 时段开机出力从 0 直接跳到某个水平不能用正常爬坡速率限制。工程上的常见处理是给公式加一个 Big-M 项允许停机状态切换到开机状态时跳过爬坡限制def ramp_up_rule(model, g, t): if t 1: return pyo.Constraint.Skip relax 2 * units[g][pmax] * (1 - model.u[g, t - 1]) return (model.p[g, t] - model.p[g, t - 1] units[g][RU] relax) model.ramp_up pyo.Constraint(model.G, model.T, ruleramp_up_rule) def ramp_down_rule(model, g, t): if t 1: return pyo.Constraint.Skip relax 2 * units[g][pmax] * (1 - model.u[g, t]) return (model.p[g, t - 1] - model.p[g, t] units[g][RU] relax) model.ramp_down pyo.Constraint(model.G, model.T, ruleramp_down_rule)爬坡约束中的 relax 项是常见教学处理当启动或停机发生时给爬坡速率上限放大一个足够大的量使其不实际起到约束作用。系数 2 乘以 pmax 可以保证大于任何可能的出力跳变。这个处理牺牲了一定的严格性但换来模型线性稳定性。节点功率平衡是直流模型的核心也是把潮流方程嵌入优化的关键。每个节点每个时段满足注入功率等于流出功率。节点 n 的注入包括连接在 n 上的机组出力和负荷流出则通过所有连接线路的直流潮流公式计算def balance_rule(model, t, n): total_gen sum(model.p[g, t] for g in model.G if units[g][bus] n) total_demand sum(demand.get(n, {}).get(t, 0) for t in model.T) if isinstance(demand.get(n), list) else 0 total_demand 0 if n in demand: total_demand demand[n][t - 1] flow_out 0 for (i, j, x, cap) in lines: if i n: flow_out (model.theta[t, n] - model.theta[t, j]) / x elif j n: flow_out (model.theta[t, n] - model.theta[t, i]) / x return total_gen - total_demand - flow_out 0 model.balance pyo.Constraint(model.B, model.T, rulebalance_rule)代码里 flow_out 的符号统一为“从节点 n 流出为正”所以表达式中机组发电减去负荷再减去流出功率等于 0。直流潮流忽略网损所有节点的有功注入之和自然等于负荷之和不需要额外的系统功率平衡约束。注意参考节点约束 theta0 必须显式添加否则相角有无穷多组解def ref_rule(model, t): return model.theta[t, 1] 0 model.ref pyo.Constraint(model.T, ruleref_rule)线路容量约束直接对每条线路的直流潮流值加上下限def flow_limit_rule(model, t, i, j): f (model.theta[t, i] - model.theta[t, j]) / next(x for (a, b, x, c) in lines if (a, b) (i, j)) return [-c for (a, b, x, c) in lines if (a, b) (i, j)][0] f \ next(c for (a, b, x, c) in lines if (a, b) (i, j))这里的写法为了演示简化实际上会变得冗长。生产代码通常把线路参数放入字典直接按 (i,j) 索引。潮流方向需要双向约束因为相角差可能是负值潮流可能反向流动。求解用 CBC 或 GLPK 都能完成小规模 MILP 几秒内收敛。solver pyo.SolverFactory(cbc) result solver.solve(model, teeTrue)3.3 交流 SCUC 的关键替换电压与无功进入约束集从直流版切到交流版不是简单换一个潮流函数而是要把整个网络约束集替换为非线性的交流潮流方程并新增电压幅值变量和节点无功平衡。机组侧需要补充无功出力上限数据这里取有功上限的 60% 作为简化。代码以 Pyomo 的 ConstraintList 方式给出差异部分model.V pyo.Var(model.B, model.T, bounds(0.95, 1.05), initialize1.0) model.Q pyo.Var(model.G, model.T, boundslambda m, g, t: (-0.6 * units[g][pmax], 0.6 * units[g][pmax])) model.ac_balance pyo.ConstraintList() for t in model.T: for n in model.B: p_inj sum(model.p[g, t] for g in model.G if units[g][bus] n) q_inj sum(model.Q[g, t] for g in model.G if units[g][bus] n) p_load demand[n][t - 1] if n in demand else 0 # 节点i注入功率的交流方程构造求和项 def p_flow_expr(n, t): expr 0 for (a, b, x, cap) in lines: if a n or b n: # G0, B-1/x 的简化交流支路模型 expr model.V[n, t]**2 * 0 # 实际应为完整支路参数这里保留展开写法便于扩展 return expr model.ac_balance.add( p_inj - p_load p_flow_expr(n, t) )这段代码故意留了扩展位置因为完整的交流支路潮流需要导纳矩阵实虚部参数。实际建模中建议用 Matpower 的 case 数据结构导出导纳矩阵再按索引注入约束。求解器换成 IPOPT并且不能再用整数变量——IPOPT 只能处理连续问题。常见的做法是先固定整数解把交流 SCUC 降为交流最优潮流问题。这也说明了为什么工程上很少直接求解完整的交流 SCUC而是把整数决策与交流潮流分开迭代。4. 求解器参数与收敛加速让 SCUC 从教学样例走向日尺度4.1 四类求解器在 DC-SCUC 与 AC-SCUC 上的边界求解器DC-SCUCAC-SCUC关键参数设置CBCMILP 可求解速度一般不支持非线性mipgap、timelimitGLPKMILP 可求解速度慢不支持非线性intopt、tmlimCPLEXMILP 专业级不支持非线性mip_tolerances_mipgap、time_limitGurobiMILP 专业级不支持非线性MIPGap、TimeLimitIPOPT不可处理整数连续 NLP 可求解tol、max_iter、bound_relax_factor直流模型在 CPLEX 或 Gurobi 上三节点系统秒级收敛IEEE 118 节点加 24 时段大约在几分钟内收敛这取决于整数变量的松弛质量和初始界。交流 SCUC 若用 IPOPT 直接处理 MINLP需要外部框架如 Pyomo 的 MindtPy做分支定界收敛速度对初值极其敏感。先解直流模型得到机组启停方案再把状态固定后交给 IPOPT 解交流潮流是工业界最常见的组合方式。求解器参数设置直接影响收敛质量和速度。Gurobi 中常见设置如下solver pyo.SolverFactory(gurobi) solver.options[MIPGap] 0.01 # 相对间隙目标1% 足够日前计划使用 solver.options[TimeLimit] 300 # 秒 solver.options[MIPFocus] 2 # 优先改进下界适合长时间尺度问题MIPGap 设成 0.01 意味着允许与最优解的偏差不超过 1%。日前市场出清对这个精度已经足够盲目追求 gap 到 0.001% 会显著拉长求解时间。TimeLimit 必须与求解周期匹配滚动调度中 5 分钟窗口内必须出结果这时候反而要提高 gap 到 3% 到 5% 换取快速响应。4.2 热启动把上一日解作为本日初始可行点求解 MILP 时一个好的整数变量初始值能把求解时间缩短一个数量级以上。常见做法是取前一日的调度计划作为当日初始点即 warm start。Pyomo 中设置变量的 start 属性并让求解器读取# 用上一日机组状态字典 initial_status 初始化 for g in model.G: for t in model.T: model.u[g, t].set_value(initial_status[g][t]) model.p[g, t].set_value(initial_power[g][t]) solver.options[MIPStart] 1 # Gurobi 读取初始整数解设置 MIPStart 后Gurobi 会先尝试修补这个不完整或不可行的解作为启发式起点。初始解不需要完全可行但要保证机组启停状态大致延续昨日形态。日尺度滚动计划中凌晨时段的负荷与昨日同时段接近热启动的有效性非常高。另一个手法是分两阶段求解先忽略网络约束只做普通机组组合拿到启停方案后再补上潮流约束。这个方法得到的整数解通常已经接近 SCUC 最优解剩下的只是调整机组出力和部分启停顺序。4.3 收敛失败诊断从参考节点到爬坡约束DC-SCUC 求解失败时首先检查平衡节点相角约束是否添加。漏掉参考节点会让潮流方程有无穷多解求解器报出“numerical issues”或“infeasible”。其次检查首时段的机组初始状态若模型假设所有机组从停机开始而实际系统中首时段已有部分机组在线会导致后续时段连锁不可行。最后是爬坡约束与最小启停时间的交叉作用机组刚启动时被最小开机时间强制保持开机同时爬坡约束又限制了出力变化幅度可能造成无解。AC-SCUC 的失败更多来自数值病态。交流潮流方程中的电压幅值在 0.95 到 1.05 之间而相角在 -0.5 到 0.5 弧度之间量级差异大容易导致求解器在违约数判断上失真。把变量按标幺值归一化、为电压设定更紧的初始界、关闭某些 IPOPT 的输出选项都是常用手段。遇到“Restoration phase failed”之类的报错优先怀疑初值质量而不是立刻调求解器参数。# 若使用 IPOPT可用以下参数提高鲁棒性 # --mu_strategy adaptive # --bound_relax_factor 1e-8 # --linear_solver mumps # 在 Pyomo 中对应 # solver.options[mu_strategy] adaptiveMUMPS 作为 IPOPT 默认线性求解器在多数场合可用但遇到大规模稀疏矩阵换用 MA57 或 MA27如许可证允许能明显改善收敛速度。交流潮流方程是等式非线性约束IPOPT 默认用的二阶导数近似在高曲率区域可能失真适度放宽 tol 到 1e-5 有助于收敛但会降低节点功率平衡的精度校核时要同步观察最大不平衡量。5. 快速验证技巧先用直流算完再用交流潮流逐点校核拿到一个 SCUC 模型包验证正确性的最快路径不是拿它跑一个大系统而是用“直流定计划、交流校核”两阶段法。第一步用直流 SCUC 求出机组的启停状态和出力计划第二步固定这些状态把交流最优潮流AC-OPF当作校核工具求解连续变量模型第三步逐条线路对比交流潮流结果与线路容量上限找出越限点。若存在越限把该线路在直流模型中的容量上限按交流结果的越限比例收紧重新求解。迭代两到三轮通常能得到一个在交流潮流下安全可行的调度计划。这个迭代流程的代码结构如下# pseudo-codeDC solve - AC check - adjust limit for it in range(3): dc_results solve_scuc_dc() fix_commitment_to(dc_results) ac_results solve_acopf_fixed_commitment() violation find_worst_line_overload(ac_results) if not violation: break tighten_line_limits(violation)校核结果中有一个非常灵敏的诊断指标——节点边际电价LMP。直流模型的 LMP 只反映有功阻塞价格交流模型则包含电压和无功的边际成本信息。对比两种模型下同一节点的 LMP 差值差值超过 20% 的节点往往是电压支撑薄弱的节点或者无功补偿不足导致线路损耗异常升高。这类节点在纯直流模型中看不到任何问题却可能在真实调度中引发电压越限。最后分享一个工程上常用的小技巧不要等交流校核之后才看电压而是在直流 SCUC 求解时就把线路潮流容量上限预留 5% 的安全裕度。这样既保留了 MILP 的求解速度又能吸收交流与直流之间的大部分偏差。如果校核后仍有少量线路越限针对个别线路单独收紧容量而不是整体加裕度收敛更快经济性损失也更小。这个策略的实际效果可以在 IEEE 30 节点系统上通过对比迭代前后的机组总成本验证。本文还有配套的精品资源点击获取 交流与直流潮流双范式下安全约束机组组合建模与Pyomo实现
简介面向电力系统调度与优化研究的一份MATLAB资源包聚焦安全约束机组组合SCUC问题同时覆盖基于交流潮流方程与直流潮流方程两种建模方式契合电力系统分析中对精度与计算效率的双重需求。资源包共有9个文件其中7个m脚本构成核心代码涵盖系统模型定义、发电机出力与线路热容量约束、遗传算法/粒子群等优化求解、负荷数据预处理及结果可视化等完整流程另附1个txt说明和1个zip补充包整体仅263KB便于快速部署。该项目已有86人学习浏览适合电力工程、控制理论或计算科学方向的学生和研究人员也适合希望将人工智能技术引入电力调度优化的开发者。包内代码可帮助理解交流/直流潮流方程在SCUC中的建模差异掌握优化算法设计与MATLAB实现并可探索机器学习预测负荷、改进机组启停决策的AI融合路径对经济调度研究与智能算法实践有直接参考价值。1. 安全约束机组承诺模型交流潮流与直流潮流两种建模范式的分水岭拿到《电力系统安全约束单位承诺模型包括基于交流潮流方程和直流潮流方程的.zip》这个压缩包先别急着解压。标题里的“单位承诺”是英文 Unit Commitment 的直译国内工程与学术界通常叫“机组组合”。机组组合解决的是未来 24 小时里每台机组何时启动、何时停机、各发多少电加上“安全约束”之后问题从单纯的经济调度升级为在机组运行约束的基础上计及输电网给潮流限制保证调度结果送得出去、线路不越限。而这个标题最值得拆解的点是把交流潮流方程和直流潮流方程放在同一个模型里对照——几乎所有做电力系统优化的人都会在这里遇到分岔路口交流模型精确但非线性非凸直流模型线性可解但牺牲了电压和无功信息。这篇文章把这两条路的理论、建模代码、求解技巧和验证方法讲透适合刚接触 SCUC 的研究生、做电网调度的工程师以及想把这类开源包改造到自己项目里的算法开发人员。2. 潮流方程如何进入优化模型从交流完整形式到直流三步近似2.1 交流潮流方程的极坐标形式与其数学本质交流潮流方程描述的是节点注入功率与节点电压幅值、相角之间的非线性关系。极坐标形式下节点 i 的有功注入方程写成P_i V_i Σ V_j ( G_ij cosδ_ij B_ij sinδ_ij )无功注入方程对应Q_i V_i Σ V_j ( G_ij sinδ_ij - B_ij cosδ_ij )其中 V_i 是节点电压幅值δ_ij δ_i - δ_j 是节点间相角差G_ij 和 B_ij 是导纳矩阵的实部与虚部。这股非线性来自两处V_i 与 V_j 的乘积项以及三角函数项。当把这样的约束嵌入机组组合模型时变量集从有功功率 P 扩展为 P、Q、V、δ 四个维度机组的无功出力和节点电压上限也一并进入约束集合问题性质变成了混合整数非线性规划MINLP。这个“变身”带来的求解代价是数量级的提升。机组组合本身已包含 0-1 整数变量处理启停状态交流潮流把约束变成非凸可行域整数变量和非线性约束交织在一起全局最优解基本无望工程上只能接受局部最优或近似解。所以当你看到某个 SCUC 项目包同时给了交流和直流两个版本通常的用意是直流版跑日前计划交流版做安全校核或作为精确对照。2.2 直流潮流方程的三步近似从三角函数到线性方程直流潮流DC Power Flow不是“直流电”的潮流而是交流潮流的线性化近似。从交流支路有功潮流公式出发忽略支路电阻和并联电容支路导纳只剩电抗分量于是 B_ij -1/x_ijG_ij 0支路有功潮流简化为P_ij V_i V_j sinδ_ij / x_ij接下来做三步近似。第一步假设节点电压幅值都在额定值附近取 V_i V_j 1.0 pu消去电压乘积项。第二步假设线路两端相角差很小sinδ_ij ≈ δ_ij。第三步忽略支路电阻和所有并联支路。最终得到P_ij (δ_i - δ_j) / x_ij线路有功潮流只与两端相角差和线路电抗有关是线性关系。节点层面的功率平衡方程为P_i Σ_j (δ_i - δ_j) / x_ij写成矩阵形式就是经典的 P Bδ其中 B 是节点电纳矩阵对角线元素是连接该节点所有支路电抗倒数之和非对角线元素是 -1/x_ij。需要指定一个平衡节点作为角度参考通常取其相角为 0。这套近似在多数输电网规划与调度中精度可接受。对以有功潮流为主的机组组合问题直流潮流把非线性约束全部线性化模型退化为混合整数线性规划MILP商业求解器可以直接求出全局最优。代价是丢失了电压幅值、无功功率和网损信息当系统运行在重负荷或电压敏感区域时直流模型的调度结果可能与真实交流潮流有较大出入。2.3 两种建模范式放入 SCUC 后的差异对照对比维度交流潮流 SCUC直流潮流 SCUC变量集合有功、无功、电压幅值、相角有功、相角约束形态非线性非凸等式约束线性等式约束问题类型MINLPMILP求解难度极高一般只能得到局部最优较低可求得全局最优典型求解器IPOPT、KNITRO、CONOPTCPLEX、Gurobi、CBC、GLPK是否含电压约束是可校核电压越限否电压信息被屏蔽网损处理隐含在潮流方程中忽略产生较小误差选型逻辑很简单若要研究网络阻塞对日前市场出清的影响直流模型已经够用速度快且结果稳定若要研究电压支撑、无功优化或高新能源渗透下的电压越限问题必须上交流模型。很多生产系统采用“直流优化 交流校核”的两级架构而不是直接求交流 SCUC。这个思路会在后面代码实战与验证章节中具体展开。3. 用 Pyomo 搭建可复现的 SCUC 模型从直流版本到交流版本3.1 数据集设计三节点小系统与机组参数表教学和验证场景不需要 IEEE 300 节点那么大的规模一个能看出“安全约束起作用”的三节点系统足够。节点 1 接便宜大机组节点 2 接贵一些的中型机组节点 3 接负荷节点 1 到节点 3 的线路容量设置偏小。这样在没有网络约束时便宜的机组 1 会尽量多发电导致线路 1-3 过载加上直流潮流约束后机组 2 必须多出力才能把线路潮流压回限额以内。负荷按 24 时段设置三段式曲线夜间 150MW白天高峰 220MW晚高峰 180MW。三台机组参数如下机组节点Pmin/MWPmax/MW固定成本/元边际成本/元每MWh最小开机/h最小停机/h爬坡/MW每hG1110300400203290G2210200500352260G325100300501140线路参数1-2 电抗 0.25 pu容量 150MW1-3 电抗 0.20 pu容量 120MW2-3 电抗 0.3 pu容量 130MW。基准容量取 100MVA机组出力和负荷的标幺值由程序折算。这些参数直接写进 Python 字典便于后续扩展成读取 Excel 或 CSV 的版本。T 24 lines [ # (起点, 终点, 电抗x, 容量上限MW) (1, 2, 0.25, 150), (1, 3, 0.20, 120), (2, 3, 0.30, 130), ] # 机组用两级成本函数近似固定成本 线性边际成本 units { 1: {bus: 1, pmin: 10, pmax: 300, fixed: 400, slope: 20, UT: 3, DT: 2, RU: 90}, 2: {bus: 2, pmin: 10, pmax: 200, fixed: 500, slope: 35, UT: 2, DT: 2, RU: 60}, 3: {bus: 2, pmin: 5, pmax: 100, fixed: 300, slope: 50, UT: 1, DT: 1, RU: 40}, } demand { 3: [150, 150, 150, 160, 170, 180, 190, 200, 210, 220, 220, 210, 200, 190, 180, 180, 190, 200, 210, 220, 200, 180, 160, 150] }说明固定成本对应机组空载热备用成本边际成本对应燃煤或燃气可变成本。负荷放在节点 3节点 1 和 2 无本地负荷。线路潮流上限单位为 MW程序内部使用标幺值时需要统一除以基准功率。3.2 直流 SCUC 建模整数变量、爬坡约束与线性潮流约束下面这段 Pyomo 代码实现了直流 SCUC 的主体框架。集合 G 和 T 分别表示机组和时段变量 u、start、stop 是 0-1 变量p 是有功出力连续变量theta 是节点相角。目标函数是所有时段的固定成本和边际成本之和。import pyomo.environ as pyo model pyo.ConcreteModel() model.G pyo.Set(initializeunits.keys()) model.T pyo.Set(initializerange(1, T 1)) model.B pyo.Set(initialize[1, 2, 3]) model.L pyo.Set(initialize[(i, j) for (i, j, x, cap) in lines]) model.u pyo.Var(model.G, model.T, domainpyo.Binary) model.start pyo.Var(model.G, model.T, domainpyo.Binary) model.stop pyo.Var(model.G, model.T, domainpyo.Binary) model.p pyo.Var(model.G, model.T, domainpyo.NonNegativeReals) model.theta pyo.Var(model.B, model.T, domainpyo.Reals) def objective_rule(model): expr 0 for g in model.G: for t in model.T: expr units[g][fixed] * model.u[g, t] expr units[g][slope] * model.p[g, t] return expr model.objective pyo.Objective(ruleobjective_rule, sensepyo.minimize)目标函数里用固定成本乘以机组状态变量 u边际成本乘以有功出力 p。当固定成本被设为 0 时模型退化为纯经济负荷分配但机组频繁启停会频繁发生所以实际中不能省略固定成本项。单台机组在一个时段内的启停损耗若要更精确可以再加启动成本项这里用 start 变量乘以单位启动费用即可扩展到目标函数中。出力上下限和启停状态绑定def p_max_rule(model, g, t): return model.p[g, t] units[g][pmax] * model.u[g, t] model.pmax_con pyo.Constraint(model.G, model.T, rulep_max_rule) def p_min_rule(model, g, t): return model.p[g, t] units[g][pmin] * model.u[g, t] model.pmin_con pyo.Constraint(model.G, model.T, rulep_min_rule)启停机逻辑约束 start[t] - stop[t] u[t] - u[t-1]保证“启动”和“停机”两个二进制变量与原状态变量的切换一致。首时段用初始状态替代这里假设所有机组从零状态开始即 u[g, 0] 0。def u_transition_rule(model, g, t): if t 1: return model.start[g, t] - model.stop[g, t] model.u[g, t] return (model.start[g, t] - model.stop[g, t] model.u[g, t] - model.u[g, t - 1]) model.transition pyo.Constraint(model.G, model.T, ruleu_transition_rule)最小开机与最小停机的处理方式是 SCUC 建模中最容易写错的部分。这里采用标准展开式写法若机组在时段 t 启动则从 t 时刻起连续 UT 个时段都必须处于开机状态若在 t 时段停机则从 t 时刻起连续 DT 个时段都必须处于停机状态。def min_uptime_rule(model, g, t, k): if t k T: return pyo.Constraint.Skip return model.u[g, t k] model.start[g, t] def min_downtime_rule(model, g, t, k): if t k T: return pyo.Constraint.Skip return model.u[g, t k] 1 - model.stop[g, t] max_ut max(units[g][UT] for g in model.G) max_dt max(units[g][DT] for g in model.G) model.mut pyo.Constraint(model.G, model.T, range(max_ut 1), rulemin_uptime_rule) model.mdt pyo.Constraint(model.G, model.T, range(max_dt 1), rulemin_downtime_rule)最小开机/停机约束中的索引 k 从 0 到最大时长但只有在 t k 不超过 T 时限才生效。这种写法直观但会产生较多冗余约束在面对大规模系统时工程技术上通常改成“启动后累计状态”的动态表达减少约束数量。本案例规模小冗余可以忽略。爬坡约束需要考虑启停瞬间的特殊处理。若机组在 t-1 时段停机、t 时段开机出力从 0 直接跳到某个水平不能用正常爬坡速率限制。工程上的常见处理是给公式加一个 Big-M 项允许停机状态切换到开机状态时跳过爬坡限制def ramp_up_rule(model, g, t): if t 1: return pyo.Constraint.Skip relax 2 * units[g][pmax] * (1 - model.u[g, t - 1]) return (model.p[g, t] - model.p[g, t - 1] units[g][RU] relax) model.ramp_up pyo.Constraint(model.G, model.T, ruleramp_up_rule) def ramp_down_rule(model, g, t): if t 1: return pyo.Constraint.Skip relax 2 * units[g][pmax] * (1 - model.u[g, t]) return (model.p[g, t - 1] - model.p[g, t] units[g][RU] relax) model.ramp_down pyo.Constraint(model.G, model.T, ruleramp_down_rule)爬坡约束中的 relax 项是常见教学处理当启动或停机发生时给爬坡速率上限放大一个足够大的量使其不实际起到约束作用。系数 2 乘以 pmax 可以保证大于任何可能的出力跳变。这个处理牺牲了一定的严格性但换来模型线性稳定性。节点功率平衡是直流模型的核心也是把潮流方程嵌入优化的关键。每个节点每个时段满足注入功率等于流出功率。节点 n 的注入包括连接在 n 上的机组出力和负荷流出则通过所有连接线路的直流潮流公式计算def balance_rule(model, t, n): total_gen sum(model.p[g, t] for g in model.G if units[g][bus] n) total_demand sum(demand.get(n, {}).get(t, 0) for t in model.T) if isinstance(demand.get(n), list) else 0 total_demand 0 if n in demand: total_demand demand[n][t - 1] flow_out 0 for (i, j, x, cap) in lines: if i n: flow_out (model.theta[t, n] - model.theta[t, j]) / x elif j n: flow_out (model.theta[t, n] - model.theta[t, i]) / x return total_gen - total_demand - flow_out 0 model.balance pyo.Constraint(model.B, model.T, rulebalance_rule)代码里 flow_out 的符号统一为“从节点 n 流出为正”所以表达式中机组发电减去负荷再减去流出功率等于 0。直流潮流忽略网损所有节点的有功注入之和自然等于负荷之和不需要额外的系统功率平衡约束。注意参考节点约束 theta0 必须显式添加否则相角有无穷多组解def ref_rule(model, t): return model.theta[t, 1] 0 model.ref pyo.Constraint(model.T, ruleref_rule)线路容量约束直接对每条线路的直流潮流值加上下限def flow_limit_rule(model, t, i, j): f (model.theta[t, i] - model.theta[t, j]) / next(x for (a, b, x, c) in lines if (a, b) (i, j)) return [-c for (a, b, x, c) in lines if (a, b) (i, j)][0] f \ next(c for (a, b, x, c) in lines if (a, b) (i, j))这里的写法为了演示简化实际上会变得冗长。生产代码通常把线路参数放入字典直接按 (i,j) 索引。潮流方向需要双向约束因为相角差可能是负值潮流可能反向流动。求解用 CBC 或 GLPK 都能完成小规模 MILP 几秒内收敛。solver pyo.SolverFactory(cbc) result solver.solve(model, teeTrue)3.3 交流 SCUC 的关键替换电压与无功进入约束集从直流版切到交流版不是简单换一个潮流函数而是要把整个网络约束集替换为非线性的交流潮流方程并新增电压幅值变量和节点无功平衡。机组侧需要补充无功出力上限数据这里取有功上限的 60% 作为简化。代码以 Pyomo 的 ConstraintList 方式给出差异部分model.V pyo.Var(model.B, model.T, bounds(0.95, 1.05), initialize1.0) model.Q pyo.Var(model.G, model.T, boundslambda m, g, t: (-0.6 * units[g][pmax], 0.6 * units[g][pmax])) model.ac_balance pyo.ConstraintList() for t in model.T: for n in model.B: p_inj sum(model.p[g, t] for g in model.G if units[g][bus] n) q_inj sum(model.Q[g, t] for g in model.G if units[g][bus] n) p_load demand[n][t - 1] if n in demand else 0 # 节点i注入功率的交流方程构造求和项 def p_flow_expr(n, t): expr 0 for (a, b, x, cap) in lines: if a n or b n: # G0, B-1/x 的简化交流支路模型 expr model.V[n, t]**2 * 0 # 实际应为完整支路参数这里保留展开写法便于扩展 return expr model.ac_balance.add( p_inj - p_load p_flow_expr(n, t) )这段代码故意留了扩展位置因为完整的交流支路潮流需要导纳矩阵实虚部参数。实际建模中建议用 Matpower 的 case 数据结构导出导纳矩阵再按索引注入约束。求解器换成 IPOPT并且不能再用整数变量——IPOPT 只能处理连续问题。常见的做法是先固定整数解把交流 SCUC 降为交流最优潮流问题。这也说明了为什么工程上很少直接求解完整的交流 SCUC而是把整数决策与交流潮流分开迭代。4. 求解器参数与收敛加速让 SCUC 从教学样例走向日尺度4.1 四类求解器在 DC-SCUC 与 AC-SCUC 上的边界求解器DC-SCUCAC-SCUC关键参数设置CBCMILP 可求解速度一般不支持非线性mipgap、timelimitGLPKMILP 可求解速度慢不支持非线性intopt、tmlimCPLEXMILP 专业级不支持非线性mip_tolerances_mipgap、time_limitGurobiMILP 专业级不支持非线性MIPGap、TimeLimitIPOPT不可处理整数连续 NLP 可求解tol、max_iter、bound_relax_factor直流模型在 CPLEX 或 Gurobi 上三节点系统秒级收敛IEEE 118 节点加 24 时段大约在几分钟内收敛这取决于整数变量的松弛质量和初始界。交流 SCUC 若用 IPOPT 直接处理 MINLP需要外部框架如 Pyomo 的 MindtPy做分支定界收敛速度对初值极其敏感。先解直流模型得到机组启停方案再把状态固定后交给 IPOPT 解交流潮流是工业界最常见的组合方式。求解器参数设置直接影响收敛质量和速度。Gurobi 中常见设置如下solver pyo.SolverFactory(gurobi) solver.options[MIPGap] 0.01 # 相对间隙目标1% 足够日前计划使用 solver.options[TimeLimit] 300 # 秒 solver.options[MIPFocus] 2 # 优先改进下界适合长时间尺度问题MIPGap 设成 0.01 意味着允许与最优解的偏差不超过 1%。日前市场出清对这个精度已经足够盲目追求 gap 到 0.001% 会显著拉长求解时间。TimeLimit 必须与求解周期匹配滚动调度中 5 分钟窗口内必须出结果这时候反而要提高 gap 到 3% 到 5% 换取快速响应。4.2 热启动把上一日解作为本日初始可行点求解 MILP 时一个好的整数变量初始值能把求解时间缩短一个数量级以上。常见做法是取前一日的调度计划作为当日初始点即 warm start。Pyomo 中设置变量的 start 属性并让求解器读取# 用上一日机组状态字典 initial_status 初始化 for g in model.G: for t in model.T: model.u[g, t].set_value(initial_status[g][t]) model.p[g, t].set_value(initial_power[g][t]) solver.options[MIPStart] 1 # Gurobi 读取初始整数解设置 MIPStart 后Gurobi 会先尝试修补这个不完整或不可行的解作为启发式起点。初始解不需要完全可行但要保证机组启停状态大致延续昨日形态。日尺度滚动计划中凌晨时段的负荷与昨日同时段接近热启动的有效性非常高。另一个手法是分两阶段求解先忽略网络约束只做普通机组组合拿到启停方案后再补上潮流约束。这个方法得到的整数解通常已经接近 SCUC 最优解剩下的只是调整机组出力和部分启停顺序。4.3 收敛失败诊断从参考节点到爬坡约束DC-SCUC 求解失败时首先检查平衡节点相角约束是否添加。漏掉参考节点会让潮流方程有无穷多解求解器报出“numerical issues”或“infeasible”。其次检查首时段的机组初始状态若模型假设所有机组从停机开始而实际系统中首时段已有部分机组在线会导致后续时段连锁不可行。最后是爬坡约束与最小启停时间的交叉作用机组刚启动时被最小开机时间强制保持开机同时爬坡约束又限制了出力变化幅度可能造成无解。AC-SCUC 的失败更多来自数值病态。交流潮流方程中的电压幅值在 0.95 到 1.05 之间而相角在 -0.5 到 0.5 弧度之间量级差异大容易导致求解器在违约数判断上失真。把变量按标幺值归一化、为电压设定更紧的初始界、关闭某些 IPOPT 的输出选项都是常用手段。遇到“Restoration phase failed”之类的报错优先怀疑初值质量而不是立刻调求解器参数。# 若使用 IPOPT可用以下参数提高鲁棒性 # --mu_strategy adaptive # --bound_relax_factor 1e-8 # --linear_solver mumps # 在 Pyomo 中对应 # solver.options[mu_strategy] adaptiveMUMPS 作为 IPOPT 默认线性求解器在多数场合可用但遇到大规模稀疏矩阵换用 MA57 或 MA27如许可证允许能明显改善收敛速度。交流潮流方程是等式非线性约束IPOPT 默认用的二阶导数近似在高曲率区域可能失真适度放宽 tol 到 1e-5 有助于收敛但会降低节点功率平衡的精度校核时要同步观察最大不平衡量。5. 快速验证技巧先用直流算完再用交流潮流逐点校核拿到一个 SCUC 模型包验证正确性的最快路径不是拿它跑一个大系统而是用“直流定计划、交流校核”两阶段法。第一步用直流 SCUC 求出机组的启停状态和出力计划第二步固定这些状态把交流最优潮流AC-OPF当作校核工具求解连续变量模型第三步逐条线路对比交流潮流结果与线路容量上限找出越限点。若存在越限把该线路在直流模型中的容量上限按交流结果的越限比例收紧重新求解。迭代两到三轮通常能得到一个在交流潮流下安全可行的调度计划。这个迭代流程的代码结构如下# pseudo-codeDC solve - AC check - adjust limit for it in range(3): dc_results solve_scuc_dc() fix_commitment_to(dc_results) ac_results solve_acopf_fixed_commitment() violation find_worst_line_overload(ac_results) if not violation: break tighten_line_limits(violation)校核结果中有一个非常灵敏的诊断指标——节点边际电价LMP。直流模型的 LMP 只反映有功阻塞价格交流模型则包含电压和无功的边际成本信息。对比两种模型下同一节点的 LMP 差值差值超过 20% 的节点往往是电压支撑薄弱的节点或者无功补偿不足导致线路损耗异常升高。这类节点在纯直流模型中看不到任何问题却可能在真实调度中引发电压越限。最后分享一个工程上常用的小技巧不要等交流校核之后才看电压而是在直流 SCUC 求解时就把线路潮流容量上限预留 5% 的安全裕度。这样既保留了 MILP 的求解速度又能吸收交流与直流之间的大部分偏差。如果校核后仍有少量线路越限针对个别线路单独收紧容量而不是整体加裕度收敛更快经济性损失也更小。这个策略的实际效果可以在 IEEE 30 节点系统上通过对比迭代前后的机组总成本验证。本文还有配套的精品资源点击获取