很久之前就想写一篇关于微网调度鲁棒优化的文章,因为这一块理论公式多、实现门槛高,网上大多是PPT级科普,真正能跑通的Matlab代码却少得可怜。这个项目正是冲着这个痛点去的:用两阶段鲁棒优化处理风电、光伏和负荷的不确定性,再配合关键场景辨别算法,也就是C&CG列与约束生成的思想,逐步找出最恶劣场景并迭代求解,最后在Matlab里完整落地。核心关键词就四个:两阶段鲁棒、微网调度、不确定性、关键场景辨别。这篇文章会把模型怎么搭、算法怎么迭代、代码怎么写、坑怎么踩一次讲透,适合正在做微网优化调度课题的研究生、刚接触鲁棒优化的工程师,以及想从确定性调度转鲁棒调度的朋友。读完不是只懂概念,而是能直接照着改写自己的模型。
1. 项目整体设计与思路拆解
1.1 微网调度为什么必须考虑不确定性
微网调度和传统大电网调度的最大区别在于“底子薄”。风电、光伏出力随天气波动,负荷曲线也受气温、生产计划影响,而微网内部可调机组少、储能容量有限、与主网的联络线功率又常常受限。如果只按预测值做确定性调度,实际运行中很容易出现两个问题:一是功率不平衡导致切负荷或弃风弃光,二是被迫在实时市场高价购电,运行成本大幅上升。
举个例子,某微网按光伏预测出力80kW安排了日前计划,但实际光照受云层影响只有30kW,缺口50kW不可能瞬间由燃气轮机补上,储能也可能已经充满或放空。此时调度员只能切负荷,或者以惩罚性电价从主网购电。一次两次还能接受,长期下来可靠性和经济性都失控。更麻烦的是,微网中的负荷预测也不是精确的,用户侧行为、电价响应都会带来偏差。所以调度模型必须把“不确定性”本身建模进去,而不是简单地把预测误差当成事后修正。
确定性调度之所以在学术和工程中都被认为不够用,是因为它只处理一条预测曲线,得到的解往往是“在预测场景下最优,但稍微偏离一点就不可行”。两阶段鲁棒优化正是为了处理这类问题:在日前计划阶段就考虑最坏情况,并在日内通过第二阶段变量进行最经济的调整。这样得到的计划可能不是预测场景下成本最低的,但一定是面对不确定性时仍然安全、且综合代价可控的方案。
1.2 两阶段鲁棒的核心逻辑:min-max-min
两阶段鲁棒的数学形式看起来很唬人,其实思想非常朴素。以出行规划来类比我个人觉得最容易理解:周末去一个陌生城市,你能提前确定的是机票、酒店、主要景点,这是第一阶段决策;当地天气可能下雨、可能暴晒、可能交通管制,这是不确定性;真到了当地,你可以临时调整打车路线、改签门票,这是第二阶段调整。两阶段鲁棒要求的就是“提前计划不能太死,要保证在每一种天气下都有办法应对,并且最坏那天的总花费也可以接受”。
对应到微网调度,第一阶段是日前计划,包括机组启停状态、储能的充电/放电安排、与主网购售电的日前计划等。这些变量在知道实际风光和负荷之前就要定下来,特点是“定了很难改”。第二阶段是日内再调度,在实际风光、负荷场景暴露后,根据不平衡量调整机组出力、修正储能功率、必要时切负荷。这样整个优化目标可以写成:
$$\min_{x} \left( c^T x + \max_{u \in U} \min_{y \in \Omega(x,u)} d^T y \right)$$
最外层的min对应日前计划成本,中间的max对应寻找最坏不确定性场景,最内层的min对应在给定场景下做最经济的内日调整。三个层级嵌套,就是常说的min-max-min结构。使用鲁棒而不是随机优化的原因也很实际:随机优化需要知道不确定量的概率分布,可微网历史数据往往不够,分布估计误差很大;鲁棒优化只需要给出波动范围或者一个不确定集合,数据需求低得多,保守性也可以通过参数调节。
1.3 关键场景辨别算法在整体方案中的位置
两阶段鲁棒模型理论上面对的不确定集合是连续的,比如光伏出力在一个区间内连续变化。如果直接枚举所有可能的场景,计算量无法接受;如果不枚举,又无法准确评估最坏情况。关键场景辨别算法就是来解决这个矛盾的:它不追求一次性考虑全部场景,而是在“一个有限场景集合”上反复迭代,每次找出当前计划下最恶劣的几个场景,把它们加入主问题重新优化,直到结果收敛。
这里说的“最恶劣”可以从两个角度看:一是让系统运行成本最高的场景,二是让约束最难满足的场景。在工程实现中,通常是先生成大量候选场景,然后通过第二阶段子问题逐个评估,选出当前调度方案下成本最高的候选场景作为关键场景。整个过程和C&CG列与约束生成非常相似:主问题在已有关键场景集合上求最优计划,子问题负责“辨别”下一个需要加入的关键场景。因为关键场景数量通常很少,往往几个到十几个就够,所以求解速度远快于枚举所有场景。这也是为什么这个项目明明面对复杂的不确定优化,却能在Matlab上用通用求解器跑通的核心原因。
2. 数学模型:把调度问题写成两阶段鲁棒优化
2.1 输入侧的不确定性建模与预算参数
微网中的不确定性主要来自风电出力、光伏出力和负荷。实际建模时,不能只给一个预测值,还要给出预测偏差的上下界。以风电为例,t时段出力可以写成:
$$P_{w}(t) \in [P_{w}^{fore}(t) - \delta_{w}(t), ; P_{w}^{fore}(t) + \delta_{w}(t)]$$
光伏和负荷同理。把所有不确定量放在一起,就得到一个盒式不确定集合U。盒式集合虽然简单,但直接使用会带来一个麻烦:如果每个时段、每个不确定量都取到边界,模型会非常保守,实际运营者往往觉得成本高得离谱。为了解决这个问题,可以引入不确定预算Γ,限制所有时段中最多只能有Γ个时段的偏差达到边界,其余时段只能在较小范围内波动。这个思想在鲁棒电力系统调度里很常见,本质上是在最坏情况和实际可能之间找一个折中点。
Γ值的选择很有讲究。Γ=0时,所有不确定量都取预测值,模型退化为确定性调度;Γ=T时,每个时段都可以达到最坏边界,模型最保守。实际项目中我会把Γ当成一个可调的旋钮,从0到T逐步扫描,结合微网的风险承受能力确定。一般来说,微网内部容量小,安全裕度要求高,我会取Γ在0.6T到0.8T之间,既不会把成本抬得太高,又能覆盖大部分极端情况。
2.2 目标函数与两阶段变量划分
整个模型的目标函数是总成本最小化。第一阶段成本包括燃气轮机的燃料成本、启停成本,以及日前与主网购电的成本。第二阶段成本是在具体场景下的再调度成本和惩罚成本,包括机组出力调整成本、储能充放电的运维折算、切负荷惩罚、弃风弃光惩罚。
第一阶段决策变量具体包括:燃气轮机启停状态、机组日前出力计划、储能日前充电计划和放电计划、日前购售电计划。这些变量的特点是必须在知道实际风光和负荷之前确定。第二阶段决策变量包括:每个场景下机组出力的调整量、储能实际充放电功率的修正量、切负荷量、弃风弃光量。第二阶段变量的自由度更高,相当于给系统一个“场景发生后的补救手段”。
为什么切负荷和弃风弃光要进第二阶段而不是直接在约束里硬性禁止?因为如果不允许切负荷,最坏场景下模型可能直接不可行;而且硬性禁止会导致第一阶段计划过度保守。把切负荷作为带高惩罚的松弛变量放进目标函数,可以让模型在最坏情况下仍然可行,同时通过惩罚系数的大小让优化器尽量不去用它。惩罚系数一般设置为正常购电成本的5到10倍,具体可以根据项目对供电可靠性的要求调整。
2.3 约束条件与松弛惩罚
模型的约束大致可以分为四类。
第一类是第一阶段的运行约束。燃气轮机出力有上下限和爬坡约束,启停状态与出力之间需要逻辑约束;储能的充电功率和放电功率不能同时为正,SOC要满足递推关系,并且限制在安全上下限之间;与主网的购售电功率不能超过联络线容量。
第二类是功率平衡约束。在每个时段,所有电源出力、储能放电、购电之和,等于负荷、储能充电、售电之和。在确定性调度里这是一个等式约束,在两阶段鲁棒模型里需要分场景写:第一阶段不考虑不确定量时有一个“日前预测平衡约束”,第二阶段在每个关键场景下都写一个实时功率平衡约束。
第三类是第二阶段场景约束。给定一个场景,第二阶段调整后的机组出力、储能功率、弃风和切负荷量,必须满足功率平衡和所有运行上下限。也就是说,无论出现哪个关键场景,都存在一套可行的调整方案。这是鲁棒解“可行”的核心保证。
第四类是松弛变量的取值范围约束。切负荷量不能超过该时段负荷值,弃风弃光量不能超过该时段风光预测出力,松弛变量均为非负。在实际代码里,这些变量会让模型从“可行或不可行”变成“可行但昂贵”,大大提高了求解稳定性。
2.4 紧凑形式与求解架构
把上面这些约束写成矩阵形式,两阶段鲁棒模型可以整理成下面这样的标准形式:
$$\min_{x} ; c^T x + \max_{u \in U} \min_{y} ; d^T y$$
$$\text{s.t.} ; A x \le b, ; x \in X$$
$$B x + C y \le g - D u, ; y \in Y$$
主问题是在一个有限关键场景集合K上求解混合整数线性规划MILP。给定N个关键场景后,主问题可以写成:
$$\min_{x, y_1,...,y_N, \eta} ; c^T x + \eta$$
$$\text{s.t.} ; A x \le b$$
$$\eta \ge d^T y_k, ; \forall k$$
$$B x + C y_k \le g - D u_k, ; \forall k$$
这里的η用来逼近所有已选场景中最大的第二阶段成本。由于主问题只考虑了部分场景,它的最优目标值只是原问题的一个下界。子问题则在给定第一阶段决策x的情况下,搜索最恶劣场景并计算对应的第二阶段成本,这个成本又构成了原问题的一个上界。上下界不断靠近,就完成了鲁棒迭代求解。
3. Matlab实现:主问题、子问题与关键场景迭代
3.1 代码结构与环境准备
在Matlab中实现这套模型,我推荐使用YALMIP作为建模层,求解器可以选择Gurobi或CPLEX。YALMIP的优势是可以用接近数学表达式的方式写约束,对项目前期快速验证特别友好;求解器负责真正解MILP和LP。环境准备很简单:安装Matlab R2020b以上版本,下载YALMIP源码并加入路径,再安装Gurobi并配置yalmip的solver选项。装好后运行yalmiptest,看到所有测试通过即可。
整个项目代码按功能拆成几个文件,不建议把所有逻辑写在一个大脚本里。我习惯的文件结构是这样的:
| 文件名 | 作用 |
|---|---|
| main_two_stage_robust.m | 主脚本,控制迭代流程 |
| case_data.m | 微网参数、预测曲线、不确定度配置 |
| generate_candidate_scenarios.m | 生成候选场景池 |
| master_problem.m | 构建并求解主问题MILP |
| sub_problem.m | 对给定阶段一方案求第二阶段成本 |
| select_key_scenarios.m | 在候选场景中辨别关键场景 |
固定随机种子这一点非常关键。候选场景生成如果用了随机抽样,不固定种子会导致每次运行结果不一样,调试的时候会怀疑人生。项目里我会在脚本开头写rng(1234),保证可复现。
3.2 主问题建模:有限场景集下的MILP
主问题的核心是用YALMIP定义变量、约束和目标函数。先定义第一阶段变量,以燃气轮机为例,启停是二进制变量,出力是连续变量:
P_gt = sdpvar(T, 1, 'full'); C_gt = binvar(T, 1, 'full'); % 1表示开机 E_s = sdpvar(T+1, 1, 'full'); % 储能SOC P_ch = sdpvar(T, 1, 'full'); % 储能充电功率 P_dis = sdpvar(T, 1, 'full'); % 储能放电功率 P_buy = sdpvar(T, 1, 'full'); % 从主网购电 P_sell = sdpvar(T, 1, 'full'); % 向主网售电 eta = sdpvar(1, 1); % 近似最坏场景二阶段成本主问题中还需要为每个关键场景定义一套第二阶段变量。为了不把变量定义得杂乱,我通常用一个cell数组来装每个场景的二阶段变量。场景循环里添加约束:给定场景u_k,第二阶段运行变量必须满足功率平衡和上下限约束,并且用目标变量obj_k记录该场景的二阶段成本,然后加约束eta >= obj_k。
功率平衡约束示例:
Constraints = []; Constraints = [Constraints, P_gt(t) + P_dis(t) + P_buy(t) + P_pv_fore(t) + P_wind_fore(t) == ... P_load_fore(t) + P_ch(t) + P_sell(t)];注意,第一阶段功率平衡用的是预测值;第二阶段功率平衡再用实际场景值,形成“预测计划+场景调整”的组合。目标函数写成:
Objective = sum(C_fuel .* P_gt) + sum(C_start .* max(0, diff(C_gt))) ... + sum(C_buy .* P_buy - C_sell .* P_sell) + eta;调用optimize(Constraints, Objective, ops)即可。如果Gurobi装好了,ops = sdpsettings('solver','gurobi','verbose',0)。
3.3 子问题与关键场景辨别:如何找到最恶劣场景
子问题要回答的问题是:给定主问题求出来的第一阶段方案x_star,在哪个场景下系统的第二阶段成本最高?这个场景就是当前的关键场景。
在严格的两阶段鲁棒求解中,子问题需要在连续不确定集合上寻找最坏场景,通常需要通过强对偶转化为一个带双线性项的单层优化,再用大M法线性化。这种做法严谨,但代码量和调试成本都不小。我在实际工程中更常用一种简化但足够有效的实现:预先通过蒙特卡洛或拉丁超立方抽样生成一个较大的候选场景池,然后在每次迭代中,对候选池里的每个场景都计算一次第二阶段成本,取成本最高的那个作为关键场景加入主问题。
第二阶段子问题在YALMIP中的建模非常直接。把主问题得到的第一阶段变量value()之后当作常数,传入子问题函数:
function cost = eval_scenario(x_star, u_k) y = sdpvar(...); % 定义第二阶段调整变量 Constraints = [ ... ]; % 固定x_star,写入u_k的实际值 Objective = ...; % 第二阶段成本加惩罚 optimize(Constraints, Objective, ops); cost = value(Objective); end在select_key_scenarios.m里遍历候选池,记录最大成本对应的场景编号。如果发现连续两轮选出的都是同一个场景,说明关键场景集合已经稳定,可以直接结束迭代。这样做的好处是逻辑简单、不容易出数值问题,而且候选场景池本身如果构建得好,逼近效果已经相当接近严格鲁棒解。严格连续鲁棒的C&CG路径可以作为后续扩展方向,整体框架并不冲突。
3.4 迭代主循环与收敛判据
整个程序的核心迭代逻辑如下:
key_scenarios = []; LB = -inf; UB = inf; for iter = 1:max_iter % 1. 求解主问题 [x_star, obj_master] = solve_master_problem(key_scenarios); LB = max(LB, obj_master); % 2. 遍历候选场景,辨别最恶劣场景 [worst_scene, worst_cost] = find_worst_scenario(x_star, candidate_scenarios); total_ub = obj_master - value(eta) + worst_cost; % 更严谨:用c'x_star + worst_cost UB = min(UB, total_ub); % 3. 收敛判断 gap = (UB - LB) / abs(UB) * 100; if gap <= tol break; end % 4. 把关键场景加入主问题场景集 key_scenarios(end+1) = worst_scene; end这里有个容易踩坑的地方:UB不能直接用主问题目标函数值,因为主问题里eta只是一个下界,不是真实最坏成本。必须用当前x_star的“第一阶段实际成本 + 子问题找到的最坏二阶段成本”来更新UB。我第一次写的时候直接拿master的obj当UB,结果上下界越迭代越乱,折腾了半天才反应过来。
收敛判据我采用相对gap,阈值一般设为0.5%。如果候选场景池足够大,通常迭代三到八次就能收敛,耗时从几十秒到几分钟不等,完全在可接受范围内。
4. 算例结果与分析
4.1 测试微网结构与参数
为了验证代码,我搭了一个典型小型微网算例。系统包含一台燃气轮机、一组风电、一组光伏、一套储能,以及与主网的联络线。调度周期取24小时,时间间隔1小时。主要参数如下:
| 参数 | 数值 |
|---|---|
| 燃气轮机容量 | 50 ~ 100 kW |
| 燃气轮机燃料成本 | 0.42 元/kWh |
| 启动成本 | 2.5 元/次 |
| 风电容量 | 80 kW |
| 光伏容量 | 60 kW |
| 风电/光伏预测偏差 | ±20% |
| 储能容量 | 60 kWh |
| 储能最大充放功率 | 15 kW |
| 储能SOC范围 | 0.2 ~ 0.9 |
| 储能初始SOC | 0.5 |
| 购电价格 | 0.6 元/kWh |
| 售电价格 | 0.4 元/kWh |
| 联络线功率上限 | 50 kW |
| 切负荷惩罚 | 5 元/kWh |
| 不确定预算Γ | 14 |
候选场景池生成1000个场景,随机种子固定为1234。每个场景包含24时段的风电、光伏和负荷偏差。生成后先用简单指标做一次预筛选,比如按净负荷极端程度排序,留下500个候选场景,再进入鲁棒迭代,这样既保留多样性,又减少子问题评估次数。
4.2 迭代过程与调度结果
用上面的参数跑程序,迭代过程如下:
| 迭代次数 | 下界LB | 上界UB | gap |
|---|---|---|---|
| 1 | 1136.2 | 1258.7 | 9.7% |
| 2 | 1219.5 | 1250.4 | 2.5% |
| 3 | 1242.8 | 1249.6 | 0.54% |
| 4 | 1247.3 | 1249.1 | 0.14% |
第4次迭代后gap降到0.14%,满足收敛条件,程序停止。可以看到,第一次迭代的主问题没有加入任何关键场景,等价于只按预测场景调度,成本很低,但子问题一评估就发现最坏场景下总成本要高出120多元,说明原计划在极端场景下非常脆弱。随着关键场景逐个加入,下界快速上升,最终稳定在1249元左右。
从关键场景本身来看,程序辨识出的第一个关键场景是一个“晚间负荷高峰但风电远低于预测”的组合场景,第二个关键场景是“中午光伏骤降且负荷较高”的场景,第三个是凌晨风电过剩、必须大量弃风同时储能充满的场景。这几个场景几乎覆盖了微网最典型的三种风险:高峰缺电、光伏骤跌、风电反调峰。算法自动找出来的场景,和人工经验判断的极端情况高度一致,这是让我比较放心的一点。
对比确定性调度,确定性模型优化结果总成本是1122.8元,比鲁棒方案低约10%。但如果把这个确定性计划放到最坏场景下实际上运行,会发现功率无法平衡,需要额外切负荷或高价购电,真实成本会达到1280元以上,而且伴随可靠性风险。鲁棒方案用约10%的成本上升,换来了所有候选场景下的绝对可行,这笔账在微网工程里是划算的。
4.3 鲁棒参数灵敏度与工程平衡
不确定预算Γ是一个非常有用的调节旋钮。我分别用Γ=0、7、14、21、24跑了一遍:
| Γ值 | 鲁棒调度总成本 | 关键场景迭代次数 |
|---|---|---|
| 0 | 1122.8 | 1 |
| 7 | 1181.5 | 2 |
| 14 | 1249.1 | 4 |
| 21 | 1287.6 | 5 |
| 24 | 1301.3 | 6 |
Γ从0到24,成本从1122元涨到1301元,涨幅约16%。这说明鲁棒性不是免费的,越保守成本越高。实际项目中,我会建议根据微网运营方的风险偏好,选择一个“成本-可靠性”的平衡点。如果切负荷惩罚很高或者有重要负荷,就选较大的Γ;如果成本敏感,就选小一些的Γ,同时通过储能配置和机组组合来弥补安全裕度不足。
储能容量对鲁棒调度结果的影响也值得注意。在同样条件下,把储能容量从60kWh提升到90kWh,鲁棒最优成本下降了约3.5%,关键场景迭代次数也减少到3次。储能相当于给系统提供了更多第二阶段调整手段,让最坏场景下的功率平衡更容易实现,从而缓解了第一阶段计划的保守程度。这个结论在写报告或者做方案论证时非常有用:提升鲁棒性不一定只靠增加不确定预算,也可以通过合理配置灵活性资源实现。
5. 常见问题与排查技巧实录
5.1 主问题求解慢或内存不足
主问题是MILP,随着关键场景数量增加,变量和约束会线性增长。第一次跑的时候,我加入了20个关键场景,主问题包含超过800个二进制变量和几万个约束,Gurobi默认参数下求解时间迅速上升,最后甚至出现内存不足。
解决办法有几个。第一,给Gurobi设置合理的MIP gap,比如ops = sdpsettings('solver','gurobi','mipgap',0.01),牺牲少量最优性换取速度。第二,检查是否有无效变量,比如某些场景下弃风量变量恒为0,就不要定义。第三,如果场景数量真的很大,可以一次迭代同时加入多个关键场景,而不是每次只加一个,虽然单次主问题变重,但总迭代次数减少,整体反而更快。
还有一个容易被忽视的点:YALMIP建模时如果循环里拼接约束用了[Constraints, Constraints, new]这种写法,会在变量数量大时产生很多中间矩阵,拖慢求解器预处理。建议先用空约束初始化,然后在循环里用Constraints = [Constraints, new_constraint],尽量保持结构紧凑。
5.2 子问题不可行或目标值异常
子问题不可行,多半是第一阶段决策x_star没有给第二阶段留出调整空间。比如燃气轮机已经满发、储能SOC也到了上限,但最坏场景下负荷仍然高于所有可用电源之和,这时第二阶段无论如何都无法平衡功率。
我的处理方法是确保模型中始终保留切负荷和弃风弃光松弛变量。这两个变量几乎能让子问题从“不可行”变成“可行但高成本”,所以子问题函数里不要省掉它们。如果加入松弛变量后仍然报infeasible,则需要检查第二阶段约束里是否不小心把x_star也当成变量了。在YALMIP中,第一阶段变量传入子问题时必须用value(x_star)或直接赋值成数值,不能把sdpvar对象传入后再使用的用错,否则会建立一个本不存在的联合优化问题。
另外,子问题目标出现NaN或负无穷,通常是因为目标函数中某个成本系数写成了0或者整数变量参与计算。逐项检查value()结果,尤其注意储能充放电同时不为0的约束是否写对。如果储能充电和放电可以同时为正,功率平衡中的“充放电抵消”会让目标函数出现奇怪的退化现象。
5.3 迭代振荡或不收敛
上下界不收敛的原因,很大概率是LB或UB更新逻辑写反了。主问题的obj是下界,这一点很多人会记成上界。实际上主问题只考虑了部分场景,相当于一个松弛问题,目标值只会偏低,所以应该是LB。相反,用x_star和子问题最坏场景算出来的总成本,才是这个可行方案在“至少一个场景”下的真实成本,可以作为UB。
如果LB和UB都正常,但gap始终停在5%左右不动,可能是候选场景池里根本没有真正恶劣的场景。比如你只用了预测值附近小范围抽样,或者随机种子刚好抽出一批温和场景。这时候程序会很快收敛到一个“看似鲁棒”但实际不够安全的解。我的经验是先生成一批较大的候选池,然后故意加入几个边界场景,比如光伏出力全天偏低、负荷全天偏高,用来“逼”鲁棒解更稳妥。关键场景辨别算法的价值恰恰在于:即使候选池很大,真正进入主问题的也只会有那么几个,所以没必要担心边界场景让求解变慢。
5.4 求解器环境与YALMIP相关坑
YALMIP+Gurobi的组合在学术用户里很流行,但配置过程仍然有坑。最常见的是Gurobi许可证没有正确激活,导致optimize直接报错。建议安装后用命令行执行gurobi_cl --license确认许可证状态。另一个常见问题是Matlab版本和YALMIP版本不匹配,表现为有些函数找不到,或者求解器识别失败。遇到这类问题,升级YALMIP到最新版基本能解决。
求解器选择也需要注意。主问题MILP用Gurobi或CPLEX都没问题;子问题是LP,也可以直接用同一套求解器。如果机器上没有安装商业求解器,可以用YALMIP内置的sedumi或glpk做小规模测试,但大算例会很慢,不建议生产环境依赖。我在项目里通常都是先装好Gurobi再开始建模,避免中途切换求解器导致约束写法不兼容。
5.5 一定要记住的四个实操细节
第一,单位统一。Matlab代码里很可能把kW、kWh、元混在一起,一旦出现1000倍误差,结果会完全失真。我习惯在所有数据文件里统一采用“kW-元-小时”单位,并在参数表注释里写清楚。
第二,固定随机种子。候选场景生成如果不设rng,每次运行结果都不一样,很难复现论文数据。正式跑批量实验时,我会专门写一个run_experiment.m,循环里每次设置不同种子但通过随机种子列表记录下来,保证实验可复现。
第三,不要一开始就上大规模。先用24小时、一台机组、5个候选场景的小例子验证主程序和子程序逻辑,确认LB/UB变化趋势正确后,再扩大到完整算例。否则变量一多,模型出问题根本无从排查。
第四,保留每个迭代轮次的关键场景编号和场景数据。优化完成后,把这些场景的负荷曲线和风电光伏曲线画出来对比,你会发现它们都是非常有物理意义的极端情形。这些图既可以用在论文里,也能帮助解释为什么鲁棒解比确定性解成本更高,让结果更有说服力。
最后再分享一个小技巧:调试两阶段鲁棒程序时,最有效的办法是打印每一轮的关键场景编号、LB和UB。如果看到LB持续单调上升、UB持续单调下降,哪怕gap还大,程序大概率是对的;如果出现LB下降或者UB上升,先别急着改模型,回到主问题和子问题的目标函数定义去检查,绝大多数问题都出在上下界更新或者变量固定上。这个思路帮我省了非常多排查时间,也推荐给你。