基于MATLAB+YALMIP+Gurobi的源荷不确定性机组调度优化实现

电力系统机组调度碰上源荷不确定性，用matlab+yalmip+gurobi这条技术栈来啃，算是目前学术界和工程界都比较主流的一条路子。我前后也算趟过不少类似的优化问题，从最初的确定性经济调度，一路做到考虑风光出力波动和负荷预测误差的随机优化，踩过的坑比写过的约束条件还多。这篇文章就围绕这个标题，把这个事儿从头到尾捋一遍，包括问题怎么建模、不确定性怎么处理、yalmip里怎么声明变量、gurobi求解器怎么调参，以及实际跑算例时最容易出错的地方。

1. 内容整体设计与思路拆解

先说清楚这个题目到底在解决什么问题。传统的机组调度（Unit Commitment, UC）是在给定负荷曲线和机组参数的条件下，决定未来一段时间内（通常是24小时或更长）哪些机组开机、哪些停机、每台机组带多少出力，目标是让总发电成本最低，同时满足系统功率平衡、机组出力上下限、爬坡速率、最小启停时间等一系列约束。

但现实情况远比这个理想模型复杂得多。新能源大规模并网之后，风电和光伏的出力本身就是波动的、间歇性的，加上负荷预测不可能完全准确，源荷双侧都带有不确定性。如果你还按照确定性模型去调度，把预测值当成真实值来用，结果往往在实时运行中出现偏差——要么某些时段功率平衡被打破，要么机组爬坡跟不上实际需求变化，严重的甚至需要切负荷或者弃风弃光。这就是为什么“考虑源荷不确定性”成了这类问题里最核心的加分项。

用matlab+yalmip+gurobi这套组合，本质上就是干三件事：

• 把机组调度问题写成数学优化模型（目标函数+约束条件）
• 用yalmip作为建模语言，把数学模型翻译成求解器能识别的形式
• 调用gurobi求解混合整数线性规划（MILP）问题

这三件套装在一起，相当于你既不用手写复杂的求解算法，也不用费劲去处理大规模稀疏矩阵，只要把模型写清楚，剩下的交给求解器去算。yalmip的建模语法非常接近数学表达式的写法，这一点对于做电力系统研究的人来说特别友好，因为你可以把论文里那套数学公式几乎逐字翻译成matlab代码，不容易出现“模型写对了但代码实现变了味”的情况。

那为什么要用gurobi而不是matlab自带的linprog或者intlinprog？坦白说，机组调度问题一旦扩展到数十台机组、数百个时段，再加上不确定性场景集，模型规模会变成几千上万个变量和约束，其中还包含大量的0-1整数变量。matlab内置求解器在这种规模下往往表现吃力，而gurobi在求解MILP问题上的性能是商业求解器里数一数二的，尤其是它的分支切割算法和大规模稀疏矩阵处理能力，实测下来求解速度和稳定性都明显高一个档次，许可证对学术用户也免费，所以它在电力系统优化领域几乎是标配。

还有一个容易被忽略的点：选题里特意提到的“源荷不确定性”，决定了模型不能停留在单场景确定性阶段，你需要引入多场景随机规划、鲁棒优化或者分布鲁棒这类建模框架。不同的不确定性建模方式，带来的问题规模和求解难度完全不一样。比如简单的等概率场景集，可以保持MILP结构不变，只是把约束复制多份；但如果你采用鲁棒优化里的盒式不确定集，模型里会出现min-max结构，转成对偶后依然是线性规划或MILP，可理解和调试的难度却上升了一个台阶。所以做这个问题的第一步，不是急着写代码，而是先把不确定性的数学刻画方式定下来。

我把整体设计思路用一个简单类比来说明。确定性调度相当于你开车走一条完全已知的路，路况、车速、红绿灯时间全都预先确定，你只需要按最优路径开就行。考虑不确定性之后，相当于你开到一半可能遇到修路、堵车、天气变化，你不能再死板地按原定计划走，得预留出应对各种突发情况的调整空间。机组调度里的“预留空间”就体现在旋转备用容量、爬坡裕度、以及多场景下机组出力的可行性上。

这个类比在后面的建模细节里会反复印证。严格来说，不确定性模型并不追求“预测准”，它追求的是“无论真实情况偏离预测多少，系统都能安全经济运行”。所以你会看到，场景法需要在目标函数里对所有场景成本求期望，鲁棒法则需要保证最坏情况下依然可行。

2. 核心细节解析与实操要点

2.1 决策变量的声明与分类

用yalmip做机组调度，第一步是理清变量类型。这个看起来简单，但实际建模时很多人栽跟头就栽在变量声明不清晰上。机组调度里有两类核心决策变量：连续变量和二进制变量。

连续变量包括每台机组在每个时段的出力功率P(i,t)，以及某些模型里会用到的启动功率、停机功率等辅助量。二进制变量包括机组启停状态u(i,t)和启动操作变量start(i,t)、停机操作变量stop(i,t)。此外，如果模型里考虑了储能装置，还需要额外的充电/放电功率变量和充电状态变量。

在yalmip中，这几类变量的声明方式分别是：

% 假设有N台机组，T个时段 P = sdpvar(N, T, 'full'); % 连续变量：机组出力 u = binvar(N, T, 'full'); % 0-1变量：启停状态 start = binvar(N, T, 'full'); % 0-1变量：启动操作 stop = binvar(N, T, 'full'); % 0-1变量：停机操作

这里有一个关键细节：yalmip默认变量是对称矩阵还是全矩阵，取决于第三个参数。我以前吃过亏，sdpvar(N, T)默认生成的是对称方阵，N不等于T的时候会报错，N等于T的时候会莫名其妙引入多余的对称约束，结果模型约束数量翻倍，求解时间暴涨。正确做法是明确加上'full'参数。

还有一点，binvar和sdpvar不要混用。有些新手图省事，直接用sdpvar声明0-1变量然后靠约束u<=1、u>=0、u取整来限制，这在数值上非常不推荐。gurobi对整数变量有专门的处理策略，你强行用连续变量加取整约束，不仅模型规模变大，而且求解器很难有效剪枝，容易卡在整数性上。

2.2 目标函数：从确定性成本到随机期望成本

机组调度的目标函数一般是总成本最小化，包含燃料成本、启动成本和停机成本。最常见的燃料成本是机组出力的二次函数：

$$F_i(P_{i,t}) = a_i P_{i,t}^2 + b_i P_{i,t} + c_i$$

gurobi本身是求解线性规划和混合整数线性规划的，二次目标函数如果是凸的，也可以用gurobi的二次规划（QP或MIQP）功能直接求解。不过实际中更常见的做法是分段线性化（piecewise linearization），把二次成本曲线分成若干段，用一组连续变量加上SOS2约束或者直接线性插值来逼近。

用yalmip实现分段线性化最方便的方式是用sos2约束，或者干脆提前算好折线段的断点和斜率。我个人的习惯是用分段线性化，因为MIQP虽然gurobi能解，但实际工程和学术论文里MILP的通用性更强，而且分段线性化之后问题变纯线性，数值稳定性更好。具体做法是：

% 将每台机组的成本曲线分为K段，每段区间长度为L_k，斜率为s_k P_seg = sdpvar(N, T, K, 'full'); % 每段出力 for k = 1:K Constraints = [Constraints, 0 <= P_seg(:,:,k) <= L_k]; end % 总出力 = sum(P_seg, 3)，总成本 = sum(s_k .* P_seg, 3) + 固定成本

注意每一段的出力上限L_k需要根据机组的出力区间和分段方式提前算好，分段数K一般取3到5段就足够逼近二次曲线的精度要求。分段太多会显著增加变量和约束数量，对求解器形成不必要的负担。

当模型引入不确定性场景后，目标函数从确定性总成本变成所有场景下成本的期望值。假设有S个场景，每个场景s对应的负荷和新能源出力不同，那目标函数写成：

$$f_{total} = \sum_{s=1}^{S} \rho_s \sum_{t=1}^{T} \sum_{i=1}^{N} \left[ F_i(P_{i,t,s}) + SU_i \cdot start_{i,t} + SD_i \cdot stop_{i,t} \right]$$

这里的rho_s是场景s的概率权重，满足所有场景概率之和为1。要注意的是，机组启停变量u、start、stop是第一阶段的决策变量，在所有场景下必须保持一致，不能随场景变化而改变。而机组出力P(i,t,s)是第二阶段的决策变量，可以针对每个场景进行调整。这就是经典的“非预期性约束”（non-anticipativity constraint），在数学上体现为u和start、stop变量不带有场景下标。

在yalmip里实现这一点，只需要让u和start、stop的维度带上NT而不带S，而P的维度带上NTS。这就自然保证了启停决策在所有场景中保持一致。如果代码里不小心把u也定义成了NT*S，那模型就会变成每个场景各调各的启停计划，这在实际中毫无意义，模型也会因为约束过松得到过于乐观的结果。

2.3 约束条件：平衡约束与机组运行约束的细节

目标函数理解了，剩下的核心约束包括以下五类，每类在yalmip里都有对应的写法，但细节上有不少坑。

第一是功率平衡约束。不考虑不确定性时，每个时段系统总出力等于总负荷加上网损。网损如果不考虑，就是简单的等式约束：

Constraints = [Constraints, sum(P, 1) == Load];

其中Load是1*T的行向量。如果Load是列向量，需要转置对齐维度。yalmip的维度匹配不像matlab矩阵运算那样自动广播到所有维度，尤其在三维变量和二维变量混在一起运算时，经常出现“尺寸不符”报错。我的建议是提前把所有向量统一成行向量或者列向量，并且在每个约束添加后检查一下维度。

引入多场景后，功率平衡约束变成：

for s = 1:S Constraints = [Constraints, sum(P(:,:,s), 1) == Load(s,:)]; end

这里的Load(s,:)是场景s下的负荷曲线，注意必须把场景下标放在第一个维度上。

第二是机组出力上下限约束。这个约束要和机组状态变量挂钩：

$$u_{i,t} \cdot P_{i,min} \le P_{i,t} \le u_{i,t} \cdot P_{i,max}$$

在yalmip里这样写：

Constraints = [Constraints, P >= Pmin .* u, P <= Pmax .* u];

Pmin和Pmax是N1的向量，u是NT的矩阵，乘法要利用matlab的广播机制，把Pmin扩展为N*T的矩阵。好在yalmip底层支持隐式广播，但前提是维度要能对上。如果u定义的是矩阵而Pmin是列向量，yalmip会自动扩展，但出于代码清晰考虑，我通常直接repmat成矩阵再乘。

第三是爬坡约束。机组从一个时段的出力调整到下一个时段，出力变化率不能超过爬坡速率限制。这里有个常见的错误是只把爬坡约束写成：

$$P_{i,t+1} - P_{i,t} \le R_i^{up}$$

这忽略了一个重要事实：当机组在相邻时段从停机状态转为开机状态或者正好相反时，功率变化可能跨越整个出力范围，而这个范围远大于爬坡速率允许的增量。所以爬坡约束需要写成带状态变量的形式：

Constraints = [Constraints, P(:,2:T) - P(:,1:T-1) <= R_up .* u(:,1:T-1) + Pmax .* (1 - u(:,1:T-1) + start(:,2:T))]; Constraints = [Constraints, P(:,1:T-1) - P(:,2:T) <= R_down .* u(:,2:T) + Pmax .* (1 - u(:,2:T) + stop(:,1:T-1))];

这个形式里的非线性部分（两个0-1变量相乘）需要做线性化，常见的做法是引入辅助变量替换乘积项。不过初学者可以先简化处理，忽略启动/停机状态下的特殊爬坡约束，只考虑正常运行时段的爬坡限制，这样模型更干净，也更容易通过求解调试。

第四是最小启停时间约束。这是UC问题里最麻烦的约束之一，因为它天然带有“记忆性”，需要用状态变量的历史来限制当前决策。经典写法是：

% 最小开机时间 for i = 1:N for t = 2:T % 如果机组在t-1时段停机、t时段启动，那么未来MUT个时段必须保持开机 Constraints = [Constraints, u(i,t) - u(i,t-1) <= u(i,min(t+MUT(i)-1,T))]; end end

这段代码的意思是：启动动作发生意味着u(i,t)从0变到1，那么u(i,t)-u(i,t-1)=1，这个启动信号必须被未来若干时段的开机状态覆盖。这类约束的维度遍历会让模型约束数量暴涨，如果T=24、N=10，这个循环会生成几百条约束，还不算大问题；但如果T=96（15分钟一个点），N=100，这个循环会生成上万条约束，求解会明显变慢。所以实际中可以采用所谓的“紧密约束形式”（tight formulation），用累积开机时数的变量来替代这种逐条遍历的写法，效果会好很多。

2.4 源荷不确定性的三种主流建模方式

标题里的“源荷不确定性”是核心，实际工程和学术里处理它有三条常见路线，适用场景差异很大。

第一种是场景法（Scenario-based Stochastic Programming）。把风电、光伏出力和负荷预测误差用一组离散场景来表示，每个场景对应一条可能的源荷曲线，同时给每个场景分配一个概率。场景怎么来？通常用历史预测误差数据做统计拟合，然后用蒙特卡洛抽样或者拉丁超立方抽样生成大量样本，再用场景削减技术（比如快速前向选择法）把上千个场景缩减为几十个有代表性的场景。yalmip建模里，场景法相当于把每个场景下的变量都复制一份，目标函数变成所有场景成本的加权平均。好处是模型结构直观，gurobi直接能解；坏处是场景数量一多，变量规模线性增长，求解时间会快速上升。

第二种是鲁棒优化（Robust Optimization）。不确定参数用一个集合来描述，比如盒式集合或者预算椭球集合，目标是即使在最坏情况下也能满足所有约束。在机组调度里，最典型的做法是把功率平衡约束中的负荷和风电出力替换成含不确定参数的表达式，然后利用强对偶理论把内层的max-min问题转化成等价的线性约束。yalmip处理这种问题的难点在于对偶转化的过程没法直接写成一两条约束，需要你先在纸面上推导出对偶问题的紧凑形式，再手工把对偶变量和约束翻译成yalmip代码。好处是计算代价相对固定，不会因为不确定性刻画得太细导致组合爆炸；坏处是结果往往偏保守，成本会比场景法高不少。

第三种是分布鲁棒优化（Distributionally Robust Optimization）。它介于随机规划和鲁棒优化之间，假设不确定参数的真实分布落在某个以经验分布为中心的球内，优化目标是在这个分布集合内最坏情况下的期望成本。这类问题的建模通常需要引入矩信息或者Wasserstein距离的概念，最终形式往往转化为半定规划（SDP）或者带二阶锥约束的MILP。gurobi 9.0以后的版本支持二阶锥约束，所以用yalmip也能做，但建模复杂度明显比前两种高出一大截，适合有一定优化理论功底的读者。

对于刚上手的人，我建议从场景法开始，因为它的数学直观性和代码可实现性最好。先把场景法跑通，理解两阶段随机规划的变量结构，再逐渐引入鲁棒模型去对比保守性和经济性，这条学习路径实测下来最顺畅。

2.5 场景削减的实现方式

讲场景法就绕不开场景削减。很多人忽略这一步，直接抽500个场景丢给求解器，然后发现求解时间爆炸，模型根本算不动。实际上，好的场景削减可以在保证精度损失很小的情况下，把场景数量从上千降到几十，计算时间能缩短一到两个数量级。

常用的方法是快速前向选择（Fast Forward Selection），思路很直观：从原始大场景集中逐个挑选最有代表性的场景，每次选一个能使剩余场景到已选场景的最小距离最大化的场景加入备选集，然后用概率分配的方式给每个选中场景重新赋权，让它在某种意义上逼近原始场景分布。

matlab里实现这个算法需要自己写，代码量不大，大概几十行。核心步骤是：

% 假设原始场景矩阵Scen原始，维度 S_orig * (T_load + T_wind) % 计算场景之间的欧氏距离矩阵D D = squareform(pdist(Scen原始)); % 初始化选中场景集合为空，剩余场景全量 % 迭代选择，直到选中场景数量达到目标值

如果你不想自己写，也可以找现成的工具箱，比如SCENRED（场景缩减工具），但那个工具接口比较老，和现代matlab版本的兼容性一般。后来自学者用python比较多，但那又是另一套技术栈，这里就不展开。

我自己的习惯是将削减后的场景概率重新归一化，然后把概率向量作为参数传入yalmip模型。注意此时场景的数量已经大幅下降，对应的变量维度也变小，求解效率显著提升。

3. 实操过程与核心环节实现

这一节直接给出一套可以完整跑通的最小案例。案例规模这样设定：10台火电机组、24个调度时段、风速和负荷不确定性各生成50个场景然后削减到10个场景。这个规模用gurobi求解，在普通笔记本上通常几十秒到几分钟内就能算完，非常适合用来验证模型正确性和理解代码结构。

3.1 数据准备与单位归一化

做电力系统优化，单位问题是最容易出错的隐形地雷。我的建议是全部采用标幺值（p.u.），基准容量取系统最大负荷或者100 MW。为什么？因为机组参数比如Pmax、爬坡速率、成本系数原始数值大小差异很大（成本可能上千，出力可能小到几十），求解器在读入这些数据时数值尺度不一致，会导致数值稳定性下降，特别是在处理整数变量和线性约束时容易出现小数值被大数值掩盖的精度问题。

把数据归一到p.u.以后，单位成本系数也要相应折算。举个例子，如果基准容量是100 MW，某台机组Pmax是200 MW，归一化后Pmax就是2.0 p.u.，对应的成本系数a、b、c也要乘上相应的倍数，具体换算方式跟成本函数的形式有关，以二次函数F(P)=aP^2+bP+c为例，原始P单位是MW，转换为p.u.后需要重新拟合成本曲线，不能直接拿原始系数套用。

我在实际案例里直接造了一套简化数据，不追求严格物理意义，但保证模型能启动、约束能找到可行解。机组台数N=10、时段T=24，系统负荷曲线用一条典型的双峰形状，早晚两个峰值。风力场景在负荷低谷时段出力高、高峰时段出力不稳定，这样才能体现源荷不确定性的真实博弈。

3.2 yalmip模型主框架代码

下面这段代码是我实际调试后能直接跑的框架，关键部分都加了注释。

%% 参数设置 N = 10; T = 24; S = 10; % 负荷基准曲线（1*T，单位p.u.） Load = ... % 场景概率（1*S） Rho = ... % 由场景削减得到 % 场景下的负荷矩阵（S*T），场景下的风电出力矩阵（S*T） LoadScen = ... WindScen = ... %% 决策变量 P = sdpvar(N, T, S, 'full'); % 出力 u = binvar(N, T, 'full'); % 启停状态 start = binvar(N, T, 'full'); % 启动操作 stop = binvar(N, T, 'full'); % 停机操作 %% 目标函数 Cost = 0; for s = 1:S for t = 1:T for i = 1:N % 分段线性化成本，这里简化为二次直接相加，实际应用请用线性化 Cost = Cost + Rho(s) * (a(i)*P(i,t,s)^2 + b(i)*P(i,t,s) + c(i)*u(i,t)); end end end % 启停成本不随场景变化 for t = 1:T for i = 1:N Cost = Cost + SU(i)*start(i,t) + SD(i)*stop(i,t); end end %% 约束集合 Constraints = []; % 功率平衡（考虑源荷不确定性，场景s下净负荷 = 负荷 - 风电） for s = 1:S Constraints = [Constraints, sum(P(:,:,s), 1) - LoadScen(s,:) + WindScen(s,:) == 0]; end % 出力上下限 for i = 1:N for s = 1:S Constraints = [Constraints, P(i,:,s) >= Pmin(i)*u(i,:)]; Constraints = [Constraints, P(i,:,s) <= Pmax(i)*u(i,:)]; end end % 启停逻辑约束 for t = 2:T Constraints = [Constraints, u(:,t) - u(:,t-1) == start(:,t) - stop(:,t)]; end % 自变量的0-1约束 Constraints = [Constraints, start <= 1, stop <= 1]; %% 求解 ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'gurobi.MIPGap', 0.001); sol = optimize(Constraints, Cost, ops);

这里有几个细节需要注意。

第一，我故意简化了最小启停时间约束和爬坡约束的写法，目的是让主框架清晰。实际套用的时候把这些约束加回去需要额外编写循环，代码会膨胀不少。但主框架完整保证了模型是一个正确的两阶段随机UC雏形。

第二，功率平衡约束里，我把风电出力当作负的负荷来处理，这相当于假设风电全部被消纳，不考虑弃风限电。如果你想考虑弃风，需要把风电出力也设成决策变量，并且乘以一个0到1之间的消纳系数。这就又增加了一组约束和变量。初学者可以先从全消纳模型入手，跑通之后再扩展。

第三，成本函数里直接用了二次项P^2。gurobi能求解二次目标，但有个重要前提：目标函数的Hessian矩阵必须是半正定的（凸的）。机组成本函数系数a是正数，所以整个二次项是凸的，没问题。但如果你把模型扩展到含网损的潮流约束，非凸项会很麻烦，gurobi就未必能直接处理。这也是为什么很多电力系统优化文章宁可用分段线性化也不用二次目标——MILP的通用性永远是第一位的。

3.3 求解器配置与参数调整

调用gurobi的代码看着简单，但求解效果往往取决于几个关键参数怎么设置。

第一个是MIPGap。默认值是1e-4，但对于大规模场景问题，这个精度要求偏高，可能导致求解时间过长。实际工程里把MIPGap放宽到1e-3或者5e-3，成本误差几乎可以忽略，但求解时间可能缩短好几倍。

ops = sdpsettings('solver', 'gurobi', 'gurobi.MIPGap', 0.005);

第二个是TimeLimit。对于大算例，设定一个求解时间上限非常有必要。gurobi到时间后会返回当前最优可行解和gap，这样你可以得到一个“能用的解”而不是等到天荒地老。

ops = sdpsettings('solver', 'gurobi', 'gurobi.TimeLimit', 600);

第三个是线程数。gurobi默认利用所有物理核心，但在共享计算资源的服务器上，你可能希望限制线程数，避免影响其他人。

ops = sdpsettings('solver', 'gurobi', 'gurobi.Threads', 8);

另外有一个很容易被忽略的点：yalmip在调用外部求解器时，默认做了一层预处理和变量转换。有时候你明明觉得模型规模不大，但求解器报出来的变量数量比预期多很多，这可能是yalmip把非线性约束做了额外线性化。如果你想看清楚求解器到底在读什么模型，可以用export命令把模型导出成lp或mps文件，然后用文本编辑器打开检查。这一步在排查模型错误时特别有效。

export(Constraints, Cost, 'model.mps');

我第一次调试一个带场景的UC模型时，一直觉得约束逻辑没问题，但求解结果里某些时段的机组出力莫名其妙为零，查了很久没头绪。后来把模型导出成mps文件，用gurobi自带的工具读取，才发现是yalmip在做维度对齐时把一个变量矩阵转置了，导致部分场景下的约束条件张冠李戴。从那以后我养成了习惯：凡是模型行为异常，第一件事就是导出中间模型检查，而不是对着matlab代码死磕。

3.4 结果可视化与合理性检查

求解完成之后，代码还没结束。从sol里取出变量值，第一件事是画图而不是急着分析数据。看曲线能最快暴露模型异常。

% 取出机组出力 P_value = value(P); u_value = value(u); % 画典型场景下的机组出力堆叠图 for s = 1:S figure; bar(squeeze(P_value(:,:,s))', 'stacked'); hold on; plot(LoadScen(s,:), 'r-', 'LineWidth', 2); title(['场景', num2str(s), '下机组出力及负荷曲线']); end

从图里主要看三点：

• 各时段总出力（堆叠图高度）是否和该场景净负荷曲线一致，如果不一致说明功率平衡约束可能没生效
• 机组出力是否在上下限范围内，有没有越过Pmax的时段
• 机组启停状态是否合理，比如高负荷时段是不是应该有更多机组在线，低负荷时段是不是明显少开几台

另外一个很像“坑”的细节：风电出力高的时段，净负荷曲线可能低于单台最小技术出力机组的出力下限，这时候模型就找不到可行解。处理办法一般是允许弃风，或者把部分机组的最小出力下限调低，或者增加储能装置。理解了这一点，你就会明白为什么实际电力系统里需要灵活性资源——不是电网不聪明，而是物理约束摆在那里，必须有东西去“填坑”。

4. 常见问题与排查技巧实录

这一节整理我在实际跑这类模型时遇到的典型问题，每个都附上排查思路。做优化模型和写普通程序有个很大的区别：程序报错往往有明确的语法提示，但优化模型报“infeasible problem”的时候，它不会告诉你是哪条约束导致的，需要你自己去定位。这个定位过程占整个调试时间的大头。

4.1 模型无可行解（Infeasible Problem）

这是最常见的问题，尤其是第一次把不确定性场景加进去之后。原因多种多样，但按照下面的顺序排查，命中率很高。

第一步，先跑确定性模型（S=1，且不确定参数全取预测值），确认基础模型本身没问题。如果确定性模型都无解，那大概率是基础约束自相矛盾，比如负荷曲线超过了所有机组的总出力上限，或者某个时段的净负荷低于所有开机机组的最小出力之和。

第二步，确认场景削减后的场景数据没有越界。比如负荷场景里可能出现负值或者特别极端的值，影响功率平衡约束。把每个场景的负荷曲线画出来看一遍，肉眼检查比写代码检查更快。

第三步，检查机组启停变量和出力变量之间的关联约束。最常见的问题是最小启停时间约束写得过紧，导致没有机组能在某些时段满足启停序列要求。

第四步，利用yalmip的diagnostic工具查看模型信息，定位不可行约束集合。

diagnostic = optimize(Constraints, Cost, ops); if diagnostic.problem == 1 disp('模型不可行'); end

第五步，审查约束维度匹配。维度错误经常导致部分约束根本没覆盖所有时段或所有场景，造成局部时段无界或不可行。这个检查其实很烦琐，但一旦出现，往往是最耗时间的问题来源。

4.2 求解速度极慢

模型明明规模不大，但gurobi迟迟不收敛。这种情况通常有三个原因。

第一个原因是整数变量数量太多。UC模型里每个时段每台机组都有一组binvar，N=10、T=24就有240个整数变量，场景法还不会增加整数变量，但如果模型扩展了启动/停机序列变量，整数规模会更大。MILP的复杂度随整数变量数量指数增长，所以尽量保持整数变量只用于机组启停状态，不要因为建模方便给其他决策变量也套上整数属性。

第二个原因是约束里出现了大M约束但M取值过大。大M约束是线性化非线性的常用技巧，但M值过大会显著恶化线性松弛，导致分支定界效率低下。gurobi对这类约束提供的割平面效果有限，最好根据物理含义给M取一个尽量紧的上界。比如描述启动状态的约束，M取机组的Pmax就足够了，不要夸张地取成几千上万。

第三个原因是目标函数数值尺度相差太大。如果成本系数里c的量级是几千、a的量级是0.01，混合在一起数值尺度跨越五个数量级，求解器的数值容差会非常吃力。解决办法是调整单位，最好把所有成本值控制在0.1到1000这个区间内。

4.3 结果不合理：风电出力高的时段反而开更多机组

这个问题我见过不止一次，原因几乎都是忽略了旋转备用约束。实际系统里，风电出力预测高不代表真实出力一定高，为了应对风电突然降低，电网必须预留足够的上调备用，也就是多开机组保持低负荷运行，让系统具备向上调节的能力。如果模型里没有备用约束，求解器会贪图经济性，把机组几乎全停，只靠风电供负荷。一旦真实风电下降，系统根本没有调整空间。

解决方式是加旋转备用约束：

% 上调备用容量约束：每个场景每个时段开机机组的总可用上调空间 >= 备用需求 for s = 1:S for t = 1:T Constraints = [Constraints, sum(Pmax .* u(:,t), 1) - sum(P(:,t,s), 1) >= Reserve_up(t)]; end end

这一步是源荷不确定性调度和确定性调度之间最重要的区别之一，很多初学的人卡在这个地方，逻辑上说不通，本质上就是忘了“备用”这个物理需求在数学上如何表达。

4.4 场景数量太多导致计算爆炸

场景数量直接决定问题规模。S=10的模型跑得动，S=100就未必了。如果场景削减做得不够，或者不确定参数维度太高导致削减效果不好，就会面临这个困境。

有一种常见的进阶手段是Benders分解或者拉格朗日松弛，把两阶段问题拆成主问题和子问题迭代求解。但这个方向实现复杂度高，对初学者不友好。更简单的替代方案是换用鲁棒优化模型，把大量场景换成几个不确定集参数，求解速度急剧提升，代价是结果保守性增加。

我实际的做法是：先跑一个S=5的小规模场景模型，确认结果合理，再逐步增加场景到S=10、S=20，观察目标函数值的变化趋势。如果目标值趋于稳定，就说明场景数量已经足够反映不确定性特征，继续加场景只是浪费算力。这个收敛性检验在写论文时也是评委关心的重点。

5. 常见问题速查表

6. 写在最后的实操心得

回头把整个流程串一遍，这个问题的核心其实不在于matlab代码怎么写，而在于你头脑里对不确定性物理图景的把握有多清晰。机组调度不是一道静态算术题，它要回答的是“在面临一系列可能发生的未来时，我现在该做什么准备”。yalmip帮我们把数学翻译成代码，gurobi帮我们把模型变成最优解，但判断哪些不确定性需要刻画、哪些约束必须保留、结果是否违背物理常识，这些判断只能来自对电力系统运行机制本身的深入理解。

我个人的经验是，拿到一个新算例不要急着写完整代码，先用手头的简化数据跑通一个最小模型，哪怕只有3台机组、6个时段、1个场景，先把所有约束和变量全部验证正确，再逐渐扩展到完整规模。跳步的结果往往是模型bug被大算例的复杂度淹没，排查难度成倍上升。

还建议你把功率平衡、备用约束和机组出力三个核心结果的曲线放到同一张图里对照着看。做任何优化模型，只有最终结果能直观解释物理过程，才能确定这个模型是真的对了，否则gap收敛到0.1%也不能给你任何信心。

这个方向后续能扩展的空间相当大。比如加入储能系统联合调度，用价格型需求响应提高负荷灵活性，或者在鲁棒模型里引入分布鲁棒框架减少保守性。但万变不离其宗，核心还是两阶段决策框架加不确定性建模，先把这一步打扎实，后面什么花样都接得住。