赶在交稿前整理一下这段时间做的一个项目基于两阶段鲁棒优化算法的微网电源容量优化配置仿真整体在MATLAB环境里用YALMIP建模求解器走CPLEX。核心关键词是微网容量配置、两阶段鲁棒、YALMIP、CPLEX算得上电力系统优化方向里非常经典也很有代表性的一套组合。如果你正在做微网规划、储能容量优化、或者刚接触鲁棒优化想找个实际算例练手这篇应该能给你省下不少踩坑时间。先说结论两阶段鲁棒优化做微网容量配置难度不在于数学推导本身而在于怎么把min-max-min结构落地成可迭代的求解流程再在MATLAB里稳定跑通。我这次踩过的坑、试过的方法、验证过可行的参数设置下面全部展开。1. 项目整体设计与解题思路1.1 容量配置问题的本质是什么微网电源容量优化配置说白了就是回答一个问题在一个给定的微网系统里风机装多少、光伏装多少、储能配多少才能让全生命周期的投资成本和运行成本最低同时保证各种场景下供电都不掉链子。这不是简单的“多装点就完事”。设备装多了初始投资高可能利用率很低装少了遇到风光出力差、负荷尖峰的时候就只能切负荷可靠性不达标。所以这是个典型的多目标协调优化问题而且因为涉及不确定性风电、光伏出力说不准负荷也在波动直接拍脑袋定容量根本不靠谱。我做这个项目时的目标设定是微网含风力发电、光伏发电、储能系统可选配柴油发电机作为备用通过两阶段鲁棒优化确定各电源的安装容量目标函数是最小化年化投资成本与运行成本之和约束覆盖功率平衡、储能SOC连续性、机组出力上下限、联络线功率限制等。1.2 为什么选两阶段鲁棒而不是普通确定性优化确定性优化最简单给定一组典型日数据直接求最优容量。但问题在于典型日只是“平均情况”真实的风光出力可能远差于典型日你按典型日配出来的容量在极端天气下大概率不够用储能可能提前耗尽最终切负荷。随机规划也能考虑不确定性它给不确定性参数预设概率分布目标变成期望成本最小化。问题是概率分布很难准确估计而且需要大量场景才能收敛计算成本偏高。更关键的是随机规划关注的是“平均表现”不能保证极端场景下系统可行。两阶段鲁棒优化的思路不一样。它不依赖概率分布而是构造一个不确定集合uncertainty set把可能出现的风光出力、负荷波动都圈在这个集合里然后优化“最坏情况下的成本”。这就是经典的min-max-min结构第一阶段在不确定性实现前做投资决策容量配置第二阶段在不确定性实现后做运行调度经济调度目标是最小化“最坏场景下的总成本”。用大白话说鲁棒优化相当于你买保险时按最坏理赔情况来评估保费而确定性优化相当于只看平均理赔情况。前者偏保守但能保证系统在所有预设的不确定场景下都可行这对电力系统这种可靠性要求极高的场合很合适。1.3 两阶段结构具体怎么切分在容量配置这个问题里阶段划分非常自然第一阶段决策变量here-and-now风电机组装机容量、光伏装机容量、储能额定容量和额定功率这些是长期投资决策一旦定了就改不了。约束包括投资预算上限、安装场地面积限制、各设备容量上下限。第二阶段决策变量wait-and-see每个调度时段的机组出力、储能充放电功率、切负荷量、柴油机油耗等这些是在看到实际风光出力和负荷之后才能做的运行决策。目标函数相应分成两部分投资成本第一阶段加上运行成本期望中的最坏值第二阶段。运行成本用最坏情况代表正是鲁棒优化的核心特征。这个划分的逻辑在于投资决策必须提前锁定运行调度可以事后调整两阶段结构恰好模拟了这个时间先后关系。这也是为什么用两阶段鲁棒而不是单阶段鲁棒——单阶段鲁棒把容量和出力同时放在一个阶段决策过于保守而且不符合实际决策流程算出来的结果很可能是浪费投资。1.4 为什么用YALMIPCPLEX这个组合MATLABYALMIPCPLEX是电力系统优化领域事实上的标准组合我选择它有三个原因第一YALMIP建模快。YALMIP是MATLAB里的一个建模层它不需要你把优化问题手动转成标准矩阵形式直接用符号变量写目标函数和约束内部自动转换成求解器需要的格式。对于微网容量配置这种约束多、变量类型杂连续变量0-1变量的问题这个抽象层能省大量时间。第二CPLEX求解能力强。CPLEX是IBM的商业求解器对线性规划、混合整数线性规划、二次规划都有极强的大规模求解能力。两阶段鲁棒优化的子问题经过对偶变换后是个max问题主问题是带二进制变量的MILPCPLEX处理这两类问题都很快几百个节点几十个时段的算例基本秒解。第三调试和可视化方便。用YALMIP建模后可以直接输出约束数量、变量数量、求解状态等诊断信息配合MATLAB的绘图功能可以快速可视化不确定集合、迭代收敛曲线、典型日运行结果。用C语言直接调CPLEX的API调试体验会差很多。2. 两阶段鲁棒优化核心细节解析2.1 不确定集合的构建——不是越大越好两阶段鲁棒优化的关键第一步是构建不确定集合。我这次考虑三类不确定性风电出力、光伏出力、负荷需求。最常用的不确定集合是盒式集合加上预算约束。具体来说对每个时段t风电出力预测值为实际出力的波动范围为其中是波动幅度上限。盒式集合把所有时段的波动都限制在最大范围但这会过于保守——它允许所有时段同时达到最差值实际情况中不太可能所有时段都这么恶劣。所以我引入了预算约束budget of uncertainty控制最多只有个时段的波动能同时达到极端值。这个预算值是一个超参数它的取值范围和决策者的风险偏好直接相关越大越保守越小越乐观。工程实践中我建议先跑几组不同预算值的对比实验观察目标函数值和最坏场景运行成本的变化趋势再决定取值。还有个细节负荷的不确定性方向要单独处理。风光出力不确定是“实际值小于预测值”时对系统最不利负荷不确定是“实际值大于预测值”时对系统最不利所以在子问题里要让这两类不确定性的极值方向相反否则可能把不确定集合的方向搞反导致子问题算出的最坏场景其实不是最优解上的最坏场景。2.2 CCG列与约束生成算法流程两阶段鲁棒优化求解的核心是CCG算法Column-and-Constraint Generation列与约束生成。它的基本思想是把原问题拆成一个主问题和一个子问题通过不断迭代让主问题逐步逼近原问题的最优解。主问题MP只包含第一阶段变量和一个不确定场景集合初始时场景为空或只有一个预测场景是标准的MILP。子问题SP是“给定第一阶段决策变量后找到最坏不确定场景下的最小运行成本”本质上是max-min问题需要经过对偶变换转成单层max问题。CCG的迭代流程初始化设定下界LB-inf、上界UBinf迭代计数k1求解主问题MP得到第一阶段最优解和最优目标值更新LB主问题最优值将第一阶段解代入子问题SP求解最坏不确定场景子问题目标值加上第一阶段成本更新UB若UB-LB相对偏差小于阈值或达到最大迭代次数停止否则把子问题找到的最坏场景作为新场景加入主问题引入与该场景对应的第二阶段变量和约束kk1回到步骤2CCG的名字里“列生成”体现在新增场景约束“约束生成”体现在每个新场景对应的第二阶段变量和约束被添加进主问题。CCG和经典的Benders分解不同点在于Benders是往主问题加对偶割平面CCG是直接加变量和约束。实践下来CCG在两阶段鲁棒优化上收敛速度显著优于Benders尤其在第二阶段是LP的情况下往往几轮迭代就能收敛。2.3 子问题的对偶变换——max-min怎么解子问题是个max-min问题外层是找最坏不确定场景内层是给定场景下的最小运行成本。大多数求解器不能直接处理max-min结构需要把它转化为可求解的形式。我用的方法是内层LP对偶。因为内层min问题是个线性规划满足强对偶条件我可以把内层min问题写成对偶max问题然后和外层max合并成一个单层max问题。这样整个子问题就变成一个带双线性项不确定变量乘以对偶变量的优化问题。双线性项的处理有讲究。当不确定集合是预算约束盒式集合时可以通过KKT条件或线性化技巧处理双线性项。这里可以用一个不太严谨但实际很好用的说法因为在最坏场景下不确定参数倾向于取边界值所以双线性项可以通过一系列线性约束等价转换。YALMIP内置的鲁棒优化模块在某些情况下可以自动处理这个过程但手动实现CCG时建议自己构建对偶问题并线性化这样更可控。我在实际代码里用的做法是把子问题的内层min问题用YALMIP建模求解时先用dual()获取对偶变量再手工构建外层max问题。在CCG迭代初期这个对偶求解过程需要严格验证否则后面几轮迭代容易发散。2.4 为什么主问题不需要重复建模完整的运行约束刚开始接触CCG时容易犯一个错误在主问题里把所有可用时段的运行约束一次性全加进去。这样做主问题规模会爆炸而且违背了CCG“按需添加场景”的初衷。正确做法是主问题最开始只放预测场景或一个初始场景对应的运行约束随着迭代不断从子问题中提取新的最坏场景再把新场景对应的运行约束加入主问题。每个场景在主问题里都对应一套自己的第二阶段变量互不干扰这样才能保证主问题是原问题的松弛版本下界估计才有效。我在代码里用一个incremental数组存储已经生成的场景集合每轮CCG迭代后把新场景push进去主问题里用循环遍历这个场景集合生成约束。这种方式代码结构清晰也方便调试时查看每个迭代轮次新增了什么约束。3. MATLABYALMIPCPLEX实操过程3.1 环境配置与版本选型先说版本因为我在这上面栽过跟头。MATLAB我只推荐R2020a之后的版本越新越好主要原因是YALMIP对R2020a之后的接口适配更完善。CPLEX我用的IBM ILOG CPLEX Optimization Studio 12.10或20.1两者都能跑但要注意32位和64位的匹配问题——MATLAB、YALMIP、CPLEX三者必须同为64位否则会报solver not found或者直接崩溃。YALMIP的安装最推荐从GitHub仓库yalmip/YALMIP直接下载最新版然后把整个文件夹加入MATLAB路径。用老版本的话对uncertain这类函数或者optimize的求解器参数支持不好容易出怪问题。CPLEX的配置相对麻烦一点。安装完CPLEX后需要把cplex/binary/x64_win64等路径下的文件加入MATLAB路径。配置完成后在MATLAB里运行yalmiptest命令检查CPLEX是否被正确识别。这一步强烈建议做可以看到YALMIP支持的求解器列表和对应状态省得老以为问题出在自己代码上结果其实是求解器没接好。3.2 YALMIP建模的核心变量定义在YALMIP里定义变量时要区分三类第一阶段变量用sdpvar定义连续变量风机容量、光伏容量、储能容量用binvar或integer可以定义整数变量比如柴油机台数。这些变量在主问题中始终存在。第二阶段变量每个场景都有自己独立的一套变量比如该场景下每个时段的火电出力、储能充放电功率、切负荷量等。我记得YALMIP里多维变量可以直接sdpvar(nHorizon, nScenario)这样每个场景对应一列循环写约束时非常方便。不确定变量在鲁棒优化模块中可以用uncertain()声明但在手动CCG实现中不确定变量本质上也是sdpvar只是取值范围被约束在不确定集合内。手动实现时不需要额外声明只要把它们当普通连续变量处理并在约束中添加不确定集合的约束即可。这里引出一个关键选择用YALMIP内置的鲁棒优化模块还是手动实现CCG我的经验YALMIP内置的optimize(Constraints, Objective, sdpsettings(robust.lplp, enumeration))可以自动处理鲁棒问题但只适用于有一定结构的问题一旦模型复杂或者第二阶段是非线性比如有储能SOC的连续性约束、有0-1变量内置模块可能处理不了或者效率很低。手动实现CCG虽然代码量大但可控性强、调试方便、计算效率也高最终我选了手动实现。3.3 主问题MP的代码骨架主问题的建模核心是遍历已有场景集合为每个场景添加对应的运行约束。这里有个重要技巧每个场景的变量要单独定义不要复用同一组变量名。我最初图省事把所有场景共用一套变量结果主问题约束数量没增加迭代完全不收敛最后发现是场景之间变量互相覆盖了。正确做法如下% 场景集合 scenarios 是一个 cell每个元素是一个场景下所有不确定参数的值 % x_fixed 是第一阶段变量容量决策 Constraints []; Objective InvestmentCost(x_fixed); % 第一阶段投资成本 for k 1:length(scenarios) % 为场景 k 定义第二阶段变量 y sdpvar(ny, nHorizon); % 该场景下的运行调度变量 % 该场景下的运行约束 Constraints [Constraints, OperationalConstraints(y, scenarios{k})]; % 该场景下的运行成本 Objective Objective OperationCost(y); % 这里通常乘一个权重或直接求和 end % 求解主问题 ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.mip.tolerances.mipgap 0.01; % 设置MIP间隙容忍度 sol optimize(Constraints, Objective, ops);注意投资成本只在第一阶段计算一次运行成本对每个场景分别计算并累加。如果场景数量多了主问题会很大所以要靠CCG迭代控制场景数量而不是一开始就把所有潜在场景全部列进去。3.4 子问题SP的建模与对偶技巧子问题给定了第一阶段决策变量x_fixed后需要在不确定集合内寻找最坏场景并计算对应的最小运行成本。子问题原形式是max-min我通过对偶转成单层max问题。这里的关键步骤先写出内层min问题确认它是LP满足强对偶然后用YALMIP求解一次得到对偶变量再构建外层max问题。外层max的目标函数里会出现不确定变量乘以对偶变量的乘积项。举个例子如果内层约束里包含了负荷平衡约束负荷不确定参数会出现在等号右侧对偶后会在目标函数里产生“负荷值乘对偶变量”的双线性项。处理这种双线性项在盒式不确定集合下可以转换% y_dual 是内层问题某个约束的对偶变量 % d_uncertain 是不确定负荷负荷值 % 双线性项 y_dual * d_uncertain 的线性化 % 当 d_uncertain 在有界范围内变化时若 y_dual 非负 % 则最坏情况取上界若 y_dual 非正则取下界。具体处理要看对偶变量的符号方向可以用-inf和inf的约束限制辅助引入二进制变量把双线性项转为一系列线性约束或者如果对偶变量符号已知直接代入边界值这部分需要特别小心一个符号搞错最坏场景就会算错。我实际做法是先验证符号规律比如负荷平衡约束的对偶变量在正常运行时通常是正的相当于节点电价那么负荷最坏场景就会取上界这种情况直接代入上界不需要二进制变量。但如果不确定变量同时出现在多个约束中不能用这种简化必须走完整的大M法线性化。建议初学者先做单类不确定性的简化版跑通后再加复杂度。3.5 关键参数设置与计算细节CPLEX参数设置对求解性能影响非常大。我最常用的几个设置ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.mip.tolerances.mipgap 0.005; % 主问题MIP的间隙容忍度 ops.cplex.mip.tolerances.absmipgap 1e-3; % 绝对间隙 ops.cplex.timelimit 600; % 单次求解时间上限防止死循环 ops.cplex.mip.strategy.startalgorithm 1; % 初始求解算法mipgap的设置要平衡精度和速度。我之前为了追求精确解把mipgap设成0结果因为大规模MILP本身存在数值误差导致主问题的下界估计不准确反而影响了收敛性。后来放宽到0.5%到1%迭代速度明显提升最终结果跟精确解相差不到1%完全可以接受。不确定集合的预算值Γ也要设置。我做了几个实验当Γ从0增加到总时段数的一半时最坏场景运行成本会显著上升但系统总成本也会相应增加。选择一个折中的Γ值我例子中是总时段的30%左右能体现鲁棒优化的价值又能避免过度保守。数据预处理也得注意。微网优化中的数量级差异很大投资成本可能是百万级运行成本是千元/小时级储能SOC是0到1之间。建议统一到标幺值或按比例缩放否则求解器可能出现数值问题报infeasible或unbounded。我实际是把功率统一到kW或MW成本统一到千元效果稳定。4. 常见问题与排查技巧实录4.1 CPLEX识别不了或YALMIP报solver not found95%的情况是路径没配对或者版本位数不一致。先检查四件事CPLEX的bin目录是否加入MATLAB路径MATLAB和CPLEX是否同为64位是否运行过yalmiptest确认YALMIP能发现CPLEX是否存在多个版本的YALMIP被同时加入搜索路径有个坑装了Gurobi之后再去装CPLEX某些情况下YALMIP可能默认调用Gurobi而不是CPLEX。这不是错误但如果你特意想用CPLEX对比结果记得在sdpsettings里显式指定solvercplex。4.2 子问题对偶后求解报infeasible或unbounded这个我调试了比较久。原因通常是内层min问题不是严格的LP可能我建模时不小心引入了整数变量或者某个约束写反了方向导致对偶不可行。排查方法是先把内层min问题单独拿出来固定不确定参数为预测场景看看内层LP是否可行。如果内层LP本身不可行那说明运行约束建模有硬伤比如功率平衡约束里发电机出力范围没覆盖负荷水平这不怪对偶得先修正约束。内层可行但外层对偶后的max问题不可行那就要检查对偶变量符号处理是否正确双线性项线性化的符号是否弄反。另外提醒一下dual()函数获取对偶变量时必须保证上次求解使用的求解器确实返回了对偶值。CPLEX返回的dual一般是可靠的但如果内部自动用了presolve或barrier求解有些版本可能不返回对偶值这时候可以在sdpsettings里关闭某些预求解选项或者换simplex算法。4.3 CCG迭代不收敛或上下界不对齐上下界不收敛主要看两点主问题是否真的是原问题的松弛版本子问题是否真的找到了最坏场景。如果主问题没有包含所有第一阶段的运营约束那LB可能高于UB这一定是主问题漏约束了回去检查主问题场景集合更新逻辑。如果子问题在固定第一阶段变量后找不到可行解那UB会变成NaN或inf这时候要把该子问题视为不可行场景需要在主问题中添加可行性割feasibility cut这就是CCG里的另一种割平面代码里要留出这个分支。我调试时的做法每轮迭代把主问题解、子问题解、最坏场景的具体数值都打印出来手动检查几轮看是不是真实反映了物理规律。比如最坏场景下风光出力应该是下限负荷应该是上限如果算出的场景不是这样子问题建模肯定有问题。4.4 求解时间太长怎么加速两阶段鲁棒容量配置的算例规模一上来CPLEX单次求解时间可能很长。加速手段按优先级排序设置合理的mipgap0.5%到1%足够工程使用给CPLEX提供好的初始解用预测场景下的确定性优化结果作为热启动对主问题MILP开启CPLEX的剪枝策略和启发式算法如果场景多不确定集合可以只保留关键时段或者用聚类降维把24个时段减少到代表性时段储能SOC约束如果导致时间耦合可以尝试用滚动优化分解时段但容量配置问题一般不需要这么激进我最终把mipgap设为0.5%迭代次数限制在20轮以内单案例跑下来大概5到8分钟工程实践完全够用。4.5 数值问题导致的结果不合理有一次我算出来的储能容量是0.02kWh显然不合理。排查发现是储能SOC的初值没设置好导致连续性约束形同虚设储能可以白嫖能量。后来在每个场景里显式约束SOC初始值和终值相等结果就正常了。还有一个容易忽略的是可再生能源“同时出力为0”的场景如果风光全为0、负荷又是最大值柴油机又因为容量限制顶不上去系统连基准场景都可能不可行。这种情况下鲁棒优化会把容量配得偏大所以要在成本目标和可靠性之间做权衡。建议可以引入一个很小的切负荷惩罚系数让系统在最极端情况下允许少量切负荷而不是强制所有场景严格满足功率平衡。这个改进能让目标函数值更平滑迭代收敛也更快。5. 实操总结与扩展方向跑完整个项目我最大的感受是两阶段鲁棒优化做容量配置核心价值在于它把“不确定性”从一个风险因素变成了一个可量化的约束条件让你在设计阶段就能清楚地知道“如果明天风电完全不出力我现在的配置能不能撑住”。这种底气是确定性优化给不了的。对刚开始做这个方向的朋友建议从一个小规模算例入手先不急着上完整模型。我推荐先做单时段容量配置确认CCG流程清晰、对偶正确再扩展到24时段或全年8760小时。中间遇到任何调试问题优先怀疑建模而不是算法毕竟YALMIP和CPLEX这两个工具经过大量项目验证出问题的概率远小于自己代码的概率。如果要在这个方向上继续深入有几个扩展值得考虑考虑储能系统寿命损耗的容量衰减模型、多微网间的功率交互与共享储能配置、需求响应资源对容量配置的影响以及把碳交易成本纳入目标函数。这些方向在两阶段鲁棒框架下都能自然扩展只需要修改第二阶段运行约束和目标函数CCG求解流程基本不变。最后分享一个小技巧运行结果出来后画一张三维图展示不确定集合中不同预算值Γ对应的总成本与最坏场景成本变化曲线或者画一张各设备容量占比的堆叠条形图在论文和方案汇报中非常有说服力。参数敏感性分析时批量跑程序建议把每组参数的结果保存为.mat文件最后统一汇总成表格比每次重新跑要高效得多。
基于两阶段鲁棒优化的微网电源容量配置仿真
赶在交稿前整理一下这段时间做的一个项目基于两阶段鲁棒优化算法的微网电源容量优化配置仿真整体在MATLAB环境里用YALMIP建模求解器走CPLEX。核心关键词是微网容量配置、两阶段鲁棒、YALMIP、CPLEX算得上电力系统优化方向里非常经典也很有代表性的一套组合。如果你正在做微网规划、储能容量优化、或者刚接触鲁棒优化想找个实际算例练手这篇应该能给你省下不少踩坑时间。先说结论两阶段鲁棒优化做微网容量配置难度不在于数学推导本身而在于怎么把min-max-min结构落地成可迭代的求解流程再在MATLAB里稳定跑通。我这次踩过的坑、试过的方法、验证过可行的参数设置下面全部展开。1. 项目整体设计与解题思路1.1 容量配置问题的本质是什么微网电源容量优化配置说白了就是回答一个问题在一个给定的微网系统里风机装多少、光伏装多少、储能配多少才能让全生命周期的投资成本和运行成本最低同时保证各种场景下供电都不掉链子。这不是简单的“多装点就完事”。设备装多了初始投资高可能利用率很低装少了遇到风光出力差、负荷尖峰的时候就只能切负荷可靠性不达标。所以这是个典型的多目标协调优化问题而且因为涉及不确定性风电、光伏出力说不准负荷也在波动直接拍脑袋定容量根本不靠谱。我做这个项目时的目标设定是微网含风力发电、光伏发电、储能系统可选配柴油发电机作为备用通过两阶段鲁棒优化确定各电源的安装容量目标函数是最小化年化投资成本与运行成本之和约束覆盖功率平衡、储能SOC连续性、机组出力上下限、联络线功率限制等。1.2 为什么选两阶段鲁棒而不是普通确定性优化确定性优化最简单给定一组典型日数据直接求最优容量。但问题在于典型日只是“平均情况”真实的风光出力可能远差于典型日你按典型日配出来的容量在极端天气下大概率不够用储能可能提前耗尽最终切负荷。随机规划也能考虑不确定性它给不确定性参数预设概率分布目标变成期望成本最小化。问题是概率分布很难准确估计而且需要大量场景才能收敛计算成本偏高。更关键的是随机规划关注的是“平均表现”不能保证极端场景下系统可行。两阶段鲁棒优化的思路不一样。它不依赖概率分布而是构造一个不确定集合uncertainty set把可能出现的风光出力、负荷波动都圈在这个集合里然后优化“最坏情况下的成本”。这就是经典的min-max-min结构第一阶段在不确定性实现前做投资决策容量配置第二阶段在不确定性实现后做运行调度经济调度目标是最小化“最坏场景下的总成本”。用大白话说鲁棒优化相当于你买保险时按最坏理赔情况来评估保费而确定性优化相当于只看平均理赔情况。前者偏保守但能保证系统在所有预设的不确定场景下都可行这对电力系统这种可靠性要求极高的场合很合适。1.3 两阶段结构具体怎么切分在容量配置这个问题里阶段划分非常自然第一阶段决策变量here-and-now风电机组装机容量、光伏装机容量、储能额定容量和额定功率这些是长期投资决策一旦定了就改不了。约束包括投资预算上限、安装场地面积限制、各设备容量上下限。第二阶段决策变量wait-and-see每个调度时段的机组出力、储能充放电功率、切负荷量、柴油机油耗等这些是在看到实际风光出力和负荷之后才能做的运行决策。目标函数相应分成两部分投资成本第一阶段加上运行成本期望中的最坏值第二阶段。运行成本用最坏情况代表正是鲁棒优化的核心特征。这个划分的逻辑在于投资决策必须提前锁定运行调度可以事后调整两阶段结构恰好模拟了这个时间先后关系。这也是为什么用两阶段鲁棒而不是单阶段鲁棒——单阶段鲁棒把容量和出力同时放在一个阶段决策过于保守而且不符合实际决策流程算出来的结果很可能是浪费投资。1.4 为什么用YALMIPCPLEX这个组合MATLABYALMIPCPLEX是电力系统优化领域事实上的标准组合我选择它有三个原因第一YALMIP建模快。YALMIP是MATLAB里的一个建模层它不需要你把优化问题手动转成标准矩阵形式直接用符号变量写目标函数和约束内部自动转换成求解器需要的格式。对于微网容量配置这种约束多、变量类型杂连续变量0-1变量的问题这个抽象层能省大量时间。第二CPLEX求解能力强。CPLEX是IBM的商业求解器对线性规划、混合整数线性规划、二次规划都有极强的大规模求解能力。两阶段鲁棒优化的子问题经过对偶变换后是个max问题主问题是带二进制变量的MILPCPLEX处理这两类问题都很快几百个节点几十个时段的算例基本秒解。第三调试和可视化方便。用YALMIP建模后可以直接输出约束数量、变量数量、求解状态等诊断信息配合MATLAB的绘图功能可以快速可视化不确定集合、迭代收敛曲线、典型日运行结果。用C语言直接调CPLEX的API调试体验会差很多。2. 两阶段鲁棒优化核心细节解析2.1 不确定集合的构建——不是越大越好两阶段鲁棒优化的关键第一步是构建不确定集合。我这次考虑三类不确定性风电出力、光伏出力、负荷需求。最常用的不确定集合是盒式集合加上预算约束。具体来说对每个时段t风电出力预测值为实际出力的波动范围为其中是波动幅度上限。盒式集合把所有时段的波动都限制在最大范围但这会过于保守——它允许所有时段同时达到最差值实际情况中不太可能所有时段都这么恶劣。所以我引入了预算约束budget of uncertainty控制最多只有个时段的波动能同时达到极端值。这个预算值是一个超参数它的取值范围和决策者的风险偏好直接相关越大越保守越小越乐观。工程实践中我建议先跑几组不同预算值的对比实验观察目标函数值和最坏场景运行成本的变化趋势再决定取值。还有个细节负荷的不确定性方向要单独处理。风光出力不确定是“实际值小于预测值”时对系统最不利负荷不确定是“实际值大于预测值”时对系统最不利所以在子问题里要让这两类不确定性的极值方向相反否则可能把不确定集合的方向搞反导致子问题算出的最坏场景其实不是最优解上的最坏场景。2.2 CCG列与约束生成算法流程两阶段鲁棒优化求解的核心是CCG算法Column-and-Constraint Generation列与约束生成。它的基本思想是把原问题拆成一个主问题和一个子问题通过不断迭代让主问题逐步逼近原问题的最优解。主问题MP只包含第一阶段变量和一个不确定场景集合初始时场景为空或只有一个预测场景是标准的MILP。子问题SP是“给定第一阶段决策变量后找到最坏不确定场景下的最小运行成本”本质上是max-min问题需要经过对偶变换转成单层max问题。CCG的迭代流程初始化设定下界LB-inf、上界UBinf迭代计数k1求解主问题MP得到第一阶段最优解和最优目标值更新LB主问题最优值将第一阶段解代入子问题SP求解最坏不确定场景子问题目标值加上第一阶段成本更新UB若UB-LB相对偏差小于阈值或达到最大迭代次数停止否则把子问题找到的最坏场景作为新场景加入主问题引入与该场景对应的第二阶段变量和约束kk1回到步骤2CCG的名字里“列生成”体现在新增场景约束“约束生成”体现在每个新场景对应的第二阶段变量和约束被添加进主问题。CCG和经典的Benders分解不同点在于Benders是往主问题加对偶割平面CCG是直接加变量和约束。实践下来CCG在两阶段鲁棒优化上收敛速度显著优于Benders尤其在第二阶段是LP的情况下往往几轮迭代就能收敛。2.3 子问题的对偶变换——max-min怎么解子问题是个max-min问题外层是找最坏不确定场景内层是给定场景下的最小运行成本。大多数求解器不能直接处理max-min结构需要把它转化为可求解的形式。我用的方法是内层LP对偶。因为内层min问题是个线性规划满足强对偶条件我可以把内层min问题写成对偶max问题然后和外层max合并成一个单层max问题。这样整个子问题就变成一个带双线性项不确定变量乘以对偶变量的优化问题。双线性项的处理有讲究。当不确定集合是预算约束盒式集合时可以通过KKT条件或线性化技巧处理双线性项。这里可以用一个不太严谨但实际很好用的说法因为在最坏场景下不确定参数倾向于取边界值所以双线性项可以通过一系列线性约束等价转换。YALMIP内置的鲁棒优化模块在某些情况下可以自动处理这个过程但手动实现CCG时建议自己构建对偶问题并线性化这样更可控。我在实际代码里用的做法是把子问题的内层min问题用YALMIP建模求解时先用dual()获取对偶变量再手工构建外层max问题。在CCG迭代初期这个对偶求解过程需要严格验证否则后面几轮迭代容易发散。2.4 为什么主问题不需要重复建模完整的运行约束刚开始接触CCG时容易犯一个错误在主问题里把所有可用时段的运行约束一次性全加进去。这样做主问题规模会爆炸而且违背了CCG“按需添加场景”的初衷。正确做法是主问题最开始只放预测场景或一个初始场景对应的运行约束随着迭代不断从子问题中提取新的最坏场景再把新场景对应的运行约束加入主问题。每个场景在主问题里都对应一套自己的第二阶段变量互不干扰这样才能保证主问题是原问题的松弛版本下界估计才有效。我在代码里用一个incremental数组存储已经生成的场景集合每轮CCG迭代后把新场景push进去主问题里用循环遍历这个场景集合生成约束。这种方式代码结构清晰也方便调试时查看每个迭代轮次新增了什么约束。3. MATLABYALMIPCPLEX实操过程3.1 环境配置与版本选型先说版本因为我在这上面栽过跟头。MATLAB我只推荐R2020a之后的版本越新越好主要原因是YALMIP对R2020a之后的接口适配更完善。CPLEX我用的IBM ILOG CPLEX Optimization Studio 12.10或20.1两者都能跑但要注意32位和64位的匹配问题——MATLAB、YALMIP、CPLEX三者必须同为64位否则会报solver not found或者直接崩溃。YALMIP的安装最推荐从GitHub仓库yalmip/YALMIP直接下载最新版然后把整个文件夹加入MATLAB路径。用老版本的话对uncertain这类函数或者optimize的求解器参数支持不好容易出怪问题。CPLEX的配置相对麻烦一点。安装完CPLEX后需要把cplex/binary/x64_win64等路径下的文件加入MATLAB路径。配置完成后在MATLAB里运行yalmiptest命令检查CPLEX是否被正确识别。这一步强烈建议做可以看到YALMIP支持的求解器列表和对应状态省得老以为问题出在自己代码上结果其实是求解器没接好。3.2 YALMIP建模的核心变量定义在YALMIP里定义变量时要区分三类第一阶段变量用sdpvar定义连续变量风机容量、光伏容量、储能容量用binvar或integer可以定义整数变量比如柴油机台数。这些变量在主问题中始终存在。第二阶段变量每个场景都有自己独立的一套变量比如该场景下每个时段的火电出力、储能充放电功率、切负荷量等。我记得YALMIP里多维变量可以直接sdpvar(nHorizon, nScenario)这样每个场景对应一列循环写约束时非常方便。不确定变量在鲁棒优化模块中可以用uncertain()声明但在手动CCG实现中不确定变量本质上也是sdpvar只是取值范围被约束在不确定集合内。手动实现时不需要额外声明只要把它们当普通连续变量处理并在约束中添加不确定集合的约束即可。这里引出一个关键选择用YALMIP内置的鲁棒优化模块还是手动实现CCG我的经验YALMIP内置的optimize(Constraints, Objective, sdpsettings(robust.lplp, enumeration))可以自动处理鲁棒问题但只适用于有一定结构的问题一旦模型复杂或者第二阶段是非线性比如有储能SOC的连续性约束、有0-1变量内置模块可能处理不了或者效率很低。手动实现CCG虽然代码量大但可控性强、调试方便、计算效率也高最终我选了手动实现。3.3 主问题MP的代码骨架主问题的建模核心是遍历已有场景集合为每个场景添加对应的运行约束。这里有个重要技巧每个场景的变量要单独定义不要复用同一组变量名。我最初图省事把所有场景共用一套变量结果主问题约束数量没增加迭代完全不收敛最后发现是场景之间变量互相覆盖了。正确做法如下% 场景集合 scenarios 是一个 cell每个元素是一个场景下所有不确定参数的值 % x_fixed 是第一阶段变量容量决策 Constraints []; Objective InvestmentCost(x_fixed); % 第一阶段投资成本 for k 1:length(scenarios) % 为场景 k 定义第二阶段变量 y sdpvar(ny, nHorizon); % 该场景下的运行调度变量 % 该场景下的运行约束 Constraints [Constraints, OperationalConstraints(y, scenarios{k})]; % 该场景下的运行成本 Objective Objective OperationCost(y); % 这里通常乘一个权重或直接求和 end % 求解主问题 ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.mip.tolerances.mipgap 0.01; % 设置MIP间隙容忍度 sol optimize(Constraints, Objective, ops);注意投资成本只在第一阶段计算一次运行成本对每个场景分别计算并累加。如果场景数量多了主问题会很大所以要靠CCG迭代控制场景数量而不是一开始就把所有潜在场景全部列进去。3.4 子问题SP的建模与对偶技巧子问题给定了第一阶段决策变量x_fixed后需要在不确定集合内寻找最坏场景并计算对应的最小运行成本。子问题原形式是max-min我通过对偶转成单层max问题。这里的关键步骤先写出内层min问题确认它是LP满足强对偶然后用YALMIP求解一次得到对偶变量再构建外层max问题。外层max的目标函数里会出现不确定变量乘以对偶变量的乘积项。举个例子如果内层约束里包含了负荷平衡约束负荷不确定参数会出现在等号右侧对偶后会在目标函数里产生“负荷值乘对偶变量”的双线性项。处理这种双线性项在盒式不确定集合下可以转换% y_dual 是内层问题某个约束的对偶变量 % d_uncertain 是不确定负荷负荷值 % 双线性项 y_dual * d_uncertain 的线性化 % 当 d_uncertain 在有界范围内变化时若 y_dual 非负 % 则最坏情况取上界若 y_dual 非正则取下界。具体处理要看对偶变量的符号方向可以用-inf和inf的约束限制辅助引入二进制变量把双线性项转为一系列线性约束或者如果对偶变量符号已知直接代入边界值这部分需要特别小心一个符号搞错最坏场景就会算错。我实际做法是先验证符号规律比如负荷平衡约束的对偶变量在正常运行时通常是正的相当于节点电价那么负荷最坏场景就会取上界这种情况直接代入上界不需要二进制变量。但如果不确定变量同时出现在多个约束中不能用这种简化必须走完整的大M法线性化。建议初学者先做单类不确定性的简化版跑通后再加复杂度。3.5 关键参数设置与计算细节CPLEX参数设置对求解性能影响非常大。我最常用的几个设置ops sdpsettings(solver, cplex, verbose, 2); ops.cplex.mip.tolerances.mipgap 0.005; % 主问题MIP的间隙容忍度 ops.cplex.mip.tolerances.absmipgap 1e-3; % 绝对间隙 ops.cplex.timelimit 600; % 单次求解时间上限防止死循环 ops.cplex.mip.strategy.startalgorithm 1; % 初始求解算法mipgap的设置要平衡精度和速度。我之前为了追求精确解把mipgap设成0结果因为大规模MILP本身存在数值误差导致主问题的下界估计不准确反而影响了收敛性。后来放宽到0.5%到1%迭代速度明显提升最终结果跟精确解相差不到1%完全可以接受。不确定集合的预算值Γ也要设置。我做了几个实验当Γ从0增加到总时段数的一半时最坏场景运行成本会显著上升但系统总成本也会相应增加。选择一个折中的Γ值我例子中是总时段的30%左右能体现鲁棒优化的价值又能避免过度保守。数据预处理也得注意。微网优化中的数量级差异很大投资成本可能是百万级运行成本是千元/小时级储能SOC是0到1之间。建议统一到标幺值或按比例缩放否则求解器可能出现数值问题报infeasible或unbounded。我实际是把功率统一到kW或MW成本统一到千元效果稳定。4. 常见问题与排查技巧实录4.1 CPLEX识别不了或YALMIP报solver not found95%的情况是路径没配对或者版本位数不一致。先检查四件事CPLEX的bin目录是否加入MATLAB路径MATLAB和CPLEX是否同为64位是否运行过yalmiptest确认YALMIP能发现CPLEX是否存在多个版本的YALMIP被同时加入搜索路径有个坑装了Gurobi之后再去装CPLEX某些情况下YALMIP可能默认调用Gurobi而不是CPLEX。这不是错误但如果你特意想用CPLEX对比结果记得在sdpsettings里显式指定solvercplex。4.2 子问题对偶后求解报infeasible或unbounded这个我调试了比较久。原因通常是内层min问题不是严格的LP可能我建模时不小心引入了整数变量或者某个约束写反了方向导致对偶不可行。排查方法是先把内层min问题单独拿出来固定不确定参数为预测场景看看内层LP是否可行。如果内层LP本身不可行那说明运行约束建模有硬伤比如功率平衡约束里发电机出力范围没覆盖负荷水平这不怪对偶得先修正约束。内层可行但外层对偶后的max问题不可行那就要检查对偶变量符号处理是否正确双线性项线性化的符号是否弄反。另外提醒一下dual()函数获取对偶变量时必须保证上次求解使用的求解器确实返回了对偶值。CPLEX返回的dual一般是可靠的但如果内部自动用了presolve或barrier求解有些版本可能不返回对偶值这时候可以在sdpsettings里关闭某些预求解选项或者换simplex算法。4.3 CCG迭代不收敛或上下界不对齐上下界不收敛主要看两点主问题是否真的是原问题的松弛版本子问题是否真的找到了最坏场景。如果主问题没有包含所有第一阶段的运营约束那LB可能高于UB这一定是主问题漏约束了回去检查主问题场景集合更新逻辑。如果子问题在固定第一阶段变量后找不到可行解那UB会变成NaN或inf这时候要把该子问题视为不可行场景需要在主问题中添加可行性割feasibility cut这就是CCG里的另一种割平面代码里要留出这个分支。我调试时的做法每轮迭代把主问题解、子问题解、最坏场景的具体数值都打印出来手动检查几轮看是不是真实反映了物理规律。比如最坏场景下风光出力应该是下限负荷应该是上限如果算出的场景不是这样子问题建模肯定有问题。4.4 求解时间太长怎么加速两阶段鲁棒容量配置的算例规模一上来CPLEX单次求解时间可能很长。加速手段按优先级排序设置合理的mipgap0.5%到1%足够工程使用给CPLEX提供好的初始解用预测场景下的确定性优化结果作为热启动对主问题MILP开启CPLEX的剪枝策略和启发式算法如果场景多不确定集合可以只保留关键时段或者用聚类降维把24个时段减少到代表性时段储能SOC约束如果导致时间耦合可以尝试用滚动优化分解时段但容量配置问题一般不需要这么激进我最终把mipgap设为0.5%迭代次数限制在20轮以内单案例跑下来大概5到8分钟工程实践完全够用。4.5 数值问题导致的结果不合理有一次我算出来的储能容量是0.02kWh显然不合理。排查发现是储能SOC的初值没设置好导致连续性约束形同虚设储能可以白嫖能量。后来在每个场景里显式约束SOC初始值和终值相等结果就正常了。还有一个容易忽略的是可再生能源“同时出力为0”的场景如果风光全为0、负荷又是最大值柴油机又因为容量限制顶不上去系统连基准场景都可能不可行。这种情况下鲁棒优化会把容量配得偏大所以要在成本目标和可靠性之间做权衡。建议可以引入一个很小的切负荷惩罚系数让系统在最极端情况下允许少量切负荷而不是强制所有场景严格满足功率平衡。这个改进能让目标函数值更平滑迭代收敛也更快。5. 实操总结与扩展方向跑完整个项目我最大的感受是两阶段鲁棒优化做容量配置核心价值在于它把“不确定性”从一个风险因素变成了一个可量化的约束条件让你在设计阶段就能清楚地知道“如果明天风电完全不出力我现在的配置能不能撑住”。这种底气是确定性优化给不了的。对刚开始做这个方向的朋友建议从一个小规模算例入手先不急着上完整模型。我推荐先做单时段容量配置确认CCG流程清晰、对偶正确再扩展到24时段或全年8760小时。中间遇到任何调试问题优先怀疑建模而不是算法毕竟YALMIP和CPLEX这两个工具经过大量项目验证出问题的概率远小于自己代码的概率。如果要在这个方向上继续深入有几个扩展值得考虑考虑储能系统寿命损耗的容量衰减模型、多微网间的功率交互与共享储能配置、需求响应资源对容量配置的影响以及把碳交易成本纳入目标函数。这些方向在两阶段鲁棒框架下都能自然扩展只需要修改第二阶段运行约束和目标函数CCG求解流程基本不变。最后分享一个小技巧运行结果出来后画一张三维图展示不确定集合中不同预算值Γ对应的总成本与最坏场景成本变化曲线或者画一张各设备容量占比的堆叠条形图在论文和方案汇报中非常有说服力。参数敏感性分析时批量跑程序建议把每组参数的结果保存为.mat文件最后统一汇总成表格比每次重新跑要高效得多。
