简介:面向电力系统调度与优化研究人员的SCUC建模与求解资料包,聚焦安全约束机组组合问题,基于交流潮流方程精确计算与直流潮流方程快速估算相结合的方式,模拟电网运行状态并实现发电单元最优启停规划。压缩包内有9个文件,以7个MATLAB脚本(.m)为核心,涵盖模型定义、约束条件、优化求解与结果分析等环节,附带1个说明文档和1个GitHub原始压缩包,整体仅263KB,便于快速下载和查阅。目前已有86人学习浏览,内容适合电力工程、控制理论或计算科学方向的本科高年级学生和研究者作为入门参考。通过学习项目文件,读者可掌握基于MATLAB的SCUC建模流程,理解遗传算法、粒子群或线性规划等优化手段在实际调度中的应用,同时可结合人工智能方法改进负荷预测与决策策略,是一份兼具教学与科研价值的实用代码资源。
1. 做 SCUC 模型时,比“要不要启动机组”更容易踩坑的是潮流方程选型
这类安全约束机组组合(SCUC)项目最费时间的往往不是启停逻辑,而是潮流方程的接口。直流潮流方程只需要解一组线性等式,速度是够了,但没法回答“电压会不会越限”这类安全约束;交流潮流方程用电压幅值和相角把网络还原得更真实,却让整个优化变成混合整数非线性问题。这套 MATLAB 工程把两种方程放在同一个框架里:SCUC-GitHub DC目录下是可直接线性化的 DC 版本,function和AC目录下是交流校验与迭代求解的片段,example脚本适合直接套入 IEEE 节点数据。我一般会先用 DC 版本拿到可用的启停计划,再带着这个计划去交流验证,省掉的排错时间比想象的多。
2. SCUC 的数学结构和 AC/DC 潮流方程的三层耦合
2.1 目标函数和机组级约束的线性化写法
SCUC 不是单纯的经济调度,它在“发多少电”之上多了一层“这台机组开不开”的 0/1 决策。常见的目标函数写法是:
[ \min \sum_{t\in T}\sum_{i\in G}\left(c0_i u_{i,t} + c1_i P_{i,t} + SU_i y_{i,t}\right) ]
其中 (u_{i,t}) 是 0/1 状态变量,(y_{i,t}) 表示启动动作,(SU_i) 是启动成本。(P_{i,t}) 是对应时段的出力。除了负荷平衡,还要处理最小开机时间、停机时间、爬坡以及网络潮流安全约束。约束一多,靠经验判断基本失效,必须把问题组织成求解器能够接受的线性或混合整数线性形式。
很多新手第一步就错在把 (u) 和 (P) 关系写成 if 判断。实际上,intlinprog只接受矩阵形式的线性约束,所以要把“机组关闭时出力必须为 0”转换成两个不等式:
function c = build_unit_bounds(u, P, Pmin, Pmax) % u: 0/1 状态; P: 出力向量 % 约束1: P <= u * Pmax,限制出力上限 % 约束2: u * Pmin <= P,防止运行时机组低于最小技术出力 c = [P - u .* Pmax; u .* Pmin - P]; end这段代码是后面所有 SCUC 求解脚本的公共基础,它把非线性逻辑拆成了线性不等式。注意 (u_i=0) 时两个约束会强制 (P_i=0),(u_i=1) 时则退化成正常的最小/最大出力边界。实际系统里有成百上千台机组时,这个线性化会带来大量类似的稀疏行,MATLAB 里最好直接拼接稀疏矩阵,而不要用for逐行构建,否则大规模算例会慢得无法接受。
2.2 交流潮流方程如何把电压问题带进启停优化
交流潮流方程描述的是节点功率平衡:
[ P_i = V_i\sum_j V_j\left(G_{ij}\cos\theta_{ij}+B_{ij}\sin\theta_{ij}\right) ]
[ Q_i = V_i\sum_j V_j\left(G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij}\right) ]
这里 (V_i) 是节点电压幅值,(\theta_{ij}=\theta_i-\theta_j) 是相角差,(G_{ij}, B_{ij}) 来自节点导纳矩阵。这个方程组本身是非线性的,如果直接作为 SCUC 的等约束参与整数优化,会得到一个混合整数非线性规划(MINLP),求解器往往要跑很久还不一定收敛。
工程上比较实用的做法是先把交流潮流函数单独封装成校验工具,用直流 SCUC 根据 0/1 状态解出初始出力,再把这组出力带进交流潮流迭代。交流潮流最关键的一步是组装节点导纳矩阵,这部分在很多AC目录实现里都有一份类似的代码:
% branch: [from, to, r, x, 对地电纳/2] % nbus: 节点总数 Y = zeros(nbus, nbus); for k = 1:size(branch, 1) y = 1 / (branch(k,3) + 1j*branch(k,4)); f = branch(k,1); t = branch(k,2); Y(f,f) = Y(f,f) + y + 1j*branch(k,5)/2; Y(t,t) = Y(t,t) + y + 1j*branch(k,5)/2; Y(f,t) = Y(f,t) - y; Y(t,f) = Y(t,f) - y; end以上代码把每条支路的阻抗转换成导纳,并累加到对角元和非对角元,最后形成复数稀疏矩阵 (Y)。branch(k,5)通常保存线路对地电容的一半,SCUC 做稳态校验时一般影响不大,但在考虑充电功率比较明显的长线路时不能省略。得到 (Y) 之后再进行牛顿-拉夫逊迭代,每次迭代都需要解一个线性方程组,所以效率瓶颈在矩阵分解,而不是潮流公式本身。
2.3 直流潮流这条捷径到底牺牲了什么
直流潮流是交流潮流的线性近似,假设所有电压幅值近似为 1、相角差很小、支路电阻远小于电抗。最后得到:
[ P_i = \sum_j \frac{\theta_i-\theta_j}{x_{ij}} ]
写成矩阵形式就是 (P=B'\theta),其中 (B') 是只保留支路电抗倒数的节点导纳矩阵。好处非常明显:约束变成线性,和整数变量一起交给intlinprog或分支定界就能直接求解。但代价是它把无功功率和电压幅值完全剔除了,线路潮流只看有功,对高比例受端电网而言,系统可能在直流 SCUC 下“看着安全”,实际交流校验时电压已经越限。
| 对比维度 | 交流潮流方程 | 直流潮流方程 |
|---|---|---|
| 基本变量 | 电压幅值、相角、有功、无功 | 主要是有功、相角 |
| 等式类型 | 非线性 | 线性 |
| 无功与电压表达能力 | 能 | 不能 |
| 求解速度 | 慢,需迭代 | 快,一次线性求解 |
| 适合场景 | 区域网、受端网、电压敏感分析 | 主网快速分析、长周期机组组合 |
这套资源把两种模式都保留,正是为了让你在训练数据或算例规模变化时来回切换。先跑直流版本找到整数解空间的大致结构,再补交流校验,是目前很多工业级软件的做法。
3. MATLAB 下 SCUC 的工程实现:从数据矩阵到 intlinprog 与负荷预测
3.1 原始数据如何组织成求解器认识的矩阵
拿到压缩包之后,第一件事不是改目标函数,而是认清数据格式。IEEE 算例最常见的组织方式是把数据拆成bus、gen、branch、load四个部分,每行对应一条记录,每列对应一个属性。下面这个对应关系在example脚本里经常出现:
| 数据名 | 关键字段 | 在 SCUC 中的作用 |
|---|---|---|
bus | 节点编号、节点类型、基准电压 | 定义网络拓扑和电压约束 |
gen | 所属节点、Pmin/Pmax、启动成本 | 构建机组的 0/1 约束和目标函数 |
branch | 起点、终点、r/x/电纳、容量 | 组装潮流方程和线路限额 |
load | 节点负荷或时段负荷序列 | 作为常数项进入功率平衡方程 |
我一般会把四部分读入 struct 或直接读成表格,再用sparse构建节点关联矩阵。不要把这些参数硬编码,因为后期要做 N-1 安全约束或引入人工智能预测时,需要频繁替换负荷和线路参数。function目录下的公共函数往往会接收mpc结构体,里面包含mpc.bus、mpc.gen、mpc.branch,这样修改算例时主程序不用动。
3.2 两节点 DC-SCUC 的intlinprog最小工作示例
要理解整套代码,最直接的办法是用一个两节点系统把整数变量和直流潮流约束拼起来。下面这个单时段例子可实际运行,变量顺序是 [u1, u2, Pg1, Pg2, theta1],假设节点 2 为平衡节点,相角为 0。
% 两节点系统:节点1带100MW负荷,G1在节点1,G2在节点2 % 发电机数据:空载成本, 边际成本, Pmin, Pmax gen = [100, 30, 50, 200; 80, 35, 30, 150]; D1 = 100; X12 = 0.2; Fmax = 30; nvar = 5; f = [gen(1,1); gen(2,1); gen(1,2); gen(2,2); 0]; % 目标:空载+出力+相角系数0 intcon = [1 2]; % 只有状态变量是整数 Aineq = [ -gen(1,4), 0, 1, 0, 0; % Pg1 - u1*Pmax <= 0 gen(1,3), 0, -1, 0, 0; % u1*Pmin - Pg1 <= 0 0, -gen(2,4), 0, 1, 0; 0, gen(2,3), 0, -1, 0; 0, 0, 0, 0, 1/X12; % 线路有功 <= Fmax 0, 0, 0, 0, -1/X12 ]; % -线路有功 <= Fmax bineq = [0; 0; 0; 0; Fmax; Fmax]; Aeq = [0, 0, 1, 0, -1/X12; % 节点1功率平衡:Pg1 - F = D1 0, 0, 0, 1, +1/X12 ]; % 节点2功率平衡:Pg2 + F = 0 beq = [D1; 0]; lb = [0; 0; 0; 0; -inf]; ub = [1; 1; inf; inf; inf]; x = intlinprog(f, intcon, Aineq, bineq, Aeq, beq, lb, ub); if ~isempty(x) fprintf('u=[%d %d], Pg=[%.2f %.2f], theta1=%.4f\n', ... x(1), x(2), x(3), x(4), x(5)); end这里intcon专门指定前两个变量为整数 0/1;theta1虽然是连续变量,但它可以取负值,所以下界设为-inf。两行等式约束对应两个节点的直流潮流平衡,线路潮流由(theta1 - 0)/X12表示,线路容量约束被拆成上下两个不等式。这个例子没有包含启动成本和爬坡约束,但已经把安全约束、整数变量和潮流矩阵三者结合在一个问题里,后续扩展到多时段时只需要把时间维度的变量拼接在一起。
扩展多时段时,常见做法是额外引入启动变量 (y_{i,t}),并在约束中加入 (y_{i,t} \ge u_{i,t} - u_{i,t-1}),再把 (SU_i y_{i,t}) 放进目标函数。因为 (y_{i,t}) 也是 0/1 整数,所以仍然适合intlinprog。
3.3 给 SCUC 增加人工智能负荷预测的接入点
MATLAB 里的人工智能模块很适合做 SCUC 前置数据预测。调度程序通常在每天整点前拿到未来 24 小时负荷曲线,模型越准,机组组合结果越不容易因为预测偏差触发安全约束。下面是一个极简的负荷预测函数,用前几小时负荷滚动预测下一小时总负荷:
function D_next = forecast_next_hour(D_hist) % D_hist: 过去 T 小时的负荷列向量,T >= 5 D_hist = D_hist(:); T = length(D_hist); X = D_hist(1:T-1)'; % 输入:t-1 时刻负荷 Y = D_hist(2:T)'; % 目标:t 时刻负荷 net = feedforwardnet([8 4]); net.trainParam.showWindow = false; net.trainParam.epochs = 100; net = train(net, X, Y); D_next = net(D_hist(end)); end这段代码把历史负荷构造成监督学习样本:输入是前一个点的负荷,输出是下一个点的负荷,用前馈神经网络训练。feedforwardnet([8 4])表示隐藏层分别为 8 个和 4 个神经元,适合做平滑的时间序列映射。训练完成后,net(D_hist(end))就是用最近时刻预测下一时刻的总负荷。工程中不会直接用原始预测值进入 SCUC,一般会乘以一个 1%~2% 的安全系数,再把结果写回load向量,作为功率平衡约束的常数项。
需要说明的是,负荷预测误差对 SCUC 的影响往往体现在峰值时段:预测偏低会导致机组启动数量不足,交流约束越限。所以在example脚本里接入预测结果时,我会对高峰时段额外加一个保守的备用容量约束,避免优化结果过于逼近线路极限。
4. 从 DC 转到 AC 的调试清单:约束越限与求解器参数联动
4.1 为什么必须先进 DC 再切 AC
直流 SCUC 的解通常可以在几十秒内得到,但交流校验可能要反复迭代。常见流程是:先用直流模型生成一张机组启停计划,然后把发电机状态固定,让出力作为交流潮流中的给定注入功率,判断线路是否过载、节点电压是否越限。这样做的好处是避免了整数变量和非线性方程直接耦合,排错时更容易定位问题。
如果交流校验发现某条线路越限,不要立刻去改求解器的目标函数,而应该先看直流模型的线路功率和交流模型下的视在功率差在哪。直流模型忽略无功,所以交流越限通常是线路的无功流动造成,表现为末端电压偏低和视在功率超过热稳定限额。
4.2 交流越限扫描函数的复用方式
下面这段代码可以当作通用工具,放在function目录里复用:
function viol = ac_scan_violation(yij, Vf, Vt, Fmax, thresh) % yij: 线路串联导纳; Vf: 送端复电压; Vt: 受端复电压 % thresh: 过载阈值,默认1.0 S_ij = Vf * conj((Vf - Vt) * yij); overload = abs(S_ij) / Fmax; viol = overload > thresh; end这里S_ij是送端视在功率,工程上线路热稳定容量通常也是按视在功率给定,所以用abs(S_ij)与Fmax比较是合适的。thresh可以设为 0.95,表示预留 5% 的交流安全裕度。如果返回viol为逻辑数组,就可以在循环中定位到对应的线路编号,再反查这条线路两侧的发电机和母线,决定下一轮是否要让它们改变启停状态。
4.3 求解器参数和容差的工程建议
intlinprog的收敛质量取决于分支定界策略,工程上建议显式设置容差和时间上限,避免默认值把服务器拖死。下面是我常用的配置:
opt = optimoptions('intlinprog'); opt.RelativeGapTolerance = 0.01; % 相对MIP gap opt.MaxTime = 600; % 秒 opt.Display = 'iter'; [x, fval, exitflag] = intlinprog(f, intcon, Aineq, bineq, Aeq, beq, lb, ub, opt);exitflag为 1 代表找到最优解,为 -2 时通常是问题不可行,为 -3 时是时间超限。遇到不可行,先不要怀疑求解器,优先检查负荷平衡等式是否和线路容量约束冲突。交流转直流时也有关键参数:直流线路容量需要降低 5%~10%,才能留出交流校验中无功和电压下降带来的额外负担。
| 参数名 | 常见取值 | 调整目的 |
|---|---|---|
RelativeGapTolerance | 0.005~0.02 | 控制解的质量,越小越精确,但耗时增长 |
MaxTime | 300~1200 | 限制单次 SCUC 最长时间 |
dc_line_limit_factor | 0.9~0.98 | 直流模型中线路容量折扣系数 |
penalty_viol | 1000~10000 | 目标函数里的越限惩罚项系数 |
这里的penalty_viol适合在遗传算法、粒子群等启发式 SCUC 实现里使用,当线路过载时给目标函数加惩罚值,逼迫个体往安全域内移动。intlinprog不适合非线性惩罚,所以在线性版本里更多是收紧容量系数而不是加惩罚。
5. 用交流校验反哺 DC-SCUC:一个不改求解器也能提升安全性的循环技巧
既然完整 AC-SCUC 的整数规划难解,可以退一步:把交流校验当作 DC-SCUC 的“二次筛选器”,每次生成新的线路容量系数后重跑直流版本,循环两三轮就能找到兼顾计算速度和安全裕度的计划。
具体做法是:第 1 轮用原始直流容量算出机组启停;将这组出力注入交流潮流,统计越限线路;第 2 轮把这些线路的容量乘以 0.9,再跑一次 DC-SCUC。如果交流越限仍然严重,就继续收紧,最多迭代 5 次。这个过程有一个很直观的收敛信号:相邻两轮得到的启停计划不再变化。
line_limit = mpc.branch(:, 6) * 0.95; % 初始直流容量按95%取 for iter = 1:5 [Pg, status] = solve_scuc_dc(mpc, line_limit); [V, Sij, ok] = solve_ac_pf(mpc, Pg); if ok && max(abs(Sij) ./ mpc.branch(:,6)) <= 1.0 break; end for k = 1:length(line_limit) if abs(Sij(k)) > mpc.branch(k,6) line_limit(k) = 0.9 * line_limit(k); end end end这段代码中的solve_scuc_dc和solve_ac_pf可以分别映射到压缩包里 DC 目录和 AC 目录下的两个函数。循环里每次更新的是线路容量,而不是整数变量的边界,所以求解器没必要从头设计,工程改动很小。值得注意的一个细节是:不要急着把所有越限线路一次性收紧到 0.9,因为交流潮流中的线路功率受电压影响很大,可能会从一条线路转移到另一条;更稳妥的做法是每次只收紧越限最严重的 3~5 条线路,然后重新评估。如果连续三轮都收敛不了,就要回去检查直流模型的阻抗参数是否与交流分支参数不一致,尤其是串联电抗 (X) 的单位和基准功率是否匹配。这个循环技巧本质上是用交流潮流给直流模型加了一个自适应安全边际,对 IEEE 节点算例和高比例可再生能源场景都适用。
本文还有配套的精品资源,点击获取