简介本资源是一套面向电力系统优化方向研究生与科研工程师的原创微电网两阶段鲁棒优化实现代码聚焦解决含高比例可再生能源的微电网在不确定性下的经济调度难题。代码完整复现《中国电机工程学报》刘一欣文献模型并修正原文笔误基于MATLABYALMIPCPLEX框架采用列与约束生成CCG算法求解混合整数线性化的主-子问题结构支持任意随机生成的光伏与负荷场景收敛。压缩包共8个文件6个核心M脚本、1个说明文档、1个文本说明总大小272KB涵盖主程序、主/子问题建模、场景生成、KKT条件处理、结果可视化及参数配置等模块注释详尽、逻辑分层清晰、图形输出直观。已有7736人学习下载读者可直接运行获得最恶劣不确定性场景下的最优调度策略并掌握鲁棒优化建模转化、CCG迭代实现及微电网多源协同控制建模的关键技术路径。1. 项目缘起为什么微电网优化需要“两阶段”和“鲁棒性”如果你正在研究微电网的能量管理大概率已经见过不少确定性优化模型。这些模型通常假设光伏出力、负荷需求等参数是精确已知的然后求解一个“最优”调度方案。但现实情况是这些参数充满了不确定性——天气预报不准、用户用电行为随机、设备故障偶发。一个在“理想数据”下算出来的最优方案一旦遇到实际波动轻则经济性变差重则导致系统越限甚至崩溃。这就是确定性优化的致命短板。为了解决这个问题“鲁棒优化”应运而生。它的核心思想不是追求在“最理想”情况下的最优而是追求在“最恶劣”但“可能发生”的情况下的最优可行解。换句话说它要找到一个方案即使不确定性参数在其波动范围内“捣乱”我的系统依然能安全、经济地运行。这就像给调度方案穿上了一层“防弹衣”。而“两阶段”则是鲁棒优化中一个非常经典且实用的建模框架。它把决策变量分成了两类第一阶段决策和第二阶段决策。在微电网场景下这通常对应着第一阶段决策“这里-现在”决策在不确定性实际发生之前就必须做出的、不可更改的决策。比如在一天开始前与主网签订的日购电计划、储能系统的充放电计划基线、可控分布式发电机如柴油发电机的启停计划等。这些决策一旦制定在当天内调整成本很高或无法调整。第二阶段决策“等待-观望”决策在不确定性实际显现后可以快速调整的决策。比如光伏出力比预测低了我可以实时调整储能的实际放电功率、启动备用柴油机、或者向主网购买额外的实时电来弥补缺口。两阶段鲁棒优化的目标就是找到一个第一阶段决策使得无论不确定性如何在其预设的“不确定集”内变化我都能通过调整第二阶段决策来应对并且使得“第一阶段成本 最坏情况下的第二阶段成本”的总和最小。我这次分享的就是基于MATLAB平台使用YALMIP建模语言和CPLEX求解器完整实现一个微电网两阶段鲁棒优化模型的过程。代码是我自己一行行敲出来、调试通过的不仅复现了理论更解决了一系列工程实现中的实际问题。你会发现从理论公式到可运行的代码中间有不少“坑”需要填平。2. 工具箱选型为什么是MATLABYALMIPCPLEX这个“黄金三角”在动手之前工具链的选择至关重要。我选择MATLAB YALMIP CPLEX这个组合是基于多年科研和工程实践的权衡绝非随意搭配。2.1 MATLAB快速原型与算法验证的基石MATLAB在科学计算领域的地位无需多言。对于微电网优化这类涉及矩阵运算、优化算法和结果可视化的任务MATLAB提供了无与伦比的便利性。矩阵操作原生支持优化问题的约束和目标函数本质上都是矩阵和向量的运算MATLAB处理起来得心应手。丰富的内置函数从基本的数学计算到高级的绘图功能都能找到现成的、高效的函数极大加速了开发过程。调试环境友好其工作区Workspace和调试器Debugger可以让你清晰地观察每一步的变量状态对于理解复杂的优化问题迭代过程至关重要。2.2 YALMIP建模语言连接问题与求解器的桥梁这是整个项目的灵魂。YALMIP不是一个求解器而是一个建模语言和建模工具箱。它的价值在于让你可以用近乎数学公式的方式描述优化问题而无需关心底层求解器的具体调用语法。声明式建模你只需要声明变量sdpvar、定义目标函数和约束用,,等符号连接YALMIP会自动将其转化为求解器能识别的标准形式。这比直接调用求解器API写矩阵要直观无数倍。求解器无关性同一套YALMIP模型代码只需更改一个参数就可以无缝切换使用CPLEX、Gurobi、MOSEK甚至开源的GLPK等求解器。这为算法对比和方案选型提供了极大便利。高级功能支持YALMIP对鲁棒优化、双线性规划、整数规划等复杂问题有很好的内置支持特别是其robustify框架或uncertain变量声明为实现两阶段鲁棒优化提供了清晰的路径。2.3 CPLEX工业级求解器可靠性与性能的保障CPLEX是IBM旗下顶尖的商业数学规划求解器尤其在线性规划LP、混合整数线性规划MILP和二次规划QP上性能卓越。我们的两阶段鲁棒优化问题最终会通过CCG或Benders分解等算法转化为一系列主问题和子问题这些问题大多是MILP或LP正是CPLEX的强项。求解速度与稳定性对于中等及以上规模的微电网优化问题CPLEX的求解速度和找到最优解的成功率远高于许多开源求解器。强大的诊断功能当问题不可行或无界时CPLEX能提供详细的诊断信息帮助你快速定位模型中的错误比如矛盾的约束条件。学术许可对于高校和研究机构IBM通常提供免费的学术版许可使得学习和研究成本大大降低。这个组合形成了一个完美的工作流在MATLAB的舒适环境中用YALMIP像写公式一样构建模型然后调用强大的CPLEX引擎进行求解最后再用MATLAB分析和可视化结果。它平衡了开发效率、代码可读性和计算性能。3. 模型核心微电网两阶段鲁棒优化的数学骨架在写代码之前我们必须把数学模型梳理清楚。这是整个项目的“图纸”任何模糊都会导致代码错误。下面我以一个典型的并网型微电网为例拆解其两阶段鲁棒优化模型。3.1 系统结构与不确定性定义假设我们的微电网包含光伏PV、风力发电机WT、柴油发电机DG、蓄电池ESS以及从主电网购电/售电的联络线。负荷分为固定负荷和可中断负荷IL。不确定性主要来自可再生能源出力光伏P_{pv, t}^u和风电P_{wt, t}^u。其预测值为P_{pv, t}^f,P_{wt, t}^f实际值在预测值上下波动波动范围由不确定集定义。负荷需求P_{load, t}^u。同样围绕预测值P_{load, t}^f波动。我们定义一个最常用的“盒式不确定集”Box Uncertainty SetU { u_t | u_t [ΔP_{pv,t}, ΔP_{wt,t}, ΔP_{load,t}]^T, |Δ| ≤ Γ * Δ_max }其中Γ是鲁棒调节参数保守度Γ0退化为确定性优化Γ越大考虑的不确定性范围越广解越保守成本也可能越高。3.2 决策变量划分这是两阶段建模的关键第一阶段变量 (x)在不确定性揭示前决定。例如u_{dg,t}柴油发电机t时段的启停状态0/1变量。P_{dg,t}柴油发电机t时段的计划出力。P_{grid,t}^{day-ahead}与主网签订的t时段日前购电计划功率。P_{ess,t}^{sch}蓄电池t时段的计划充放电功率正为放电。SOC_t^{sch}蓄电池t时段的计划荷电状态。第二阶段变量 (y)在不确定性u_t实际发生后实时调整的变量。例如P_{pv,t}^{spill}光伏弃光量当出力过高时。P_{curt,t}可中断负荷的中断量。P_{grid,t}^{real-time}实时与主网的功率交换可能不同于日前计划。r_{dg,t}^,r_{dg,t}^-柴油发电机的向上/向下备用调用。P_{ess,t}^{adj}蓄电池的实时功率调整。3.3 目标函数与约束目标函数最小化总成本包括第一阶段成本和最坏情况下的第二阶段成本。min_x { C_day-ahead(x) max_{u∈U} min_y C_real-time(y) }其中C_day-ahead(x) 柴油发电机燃料成本 日前购电成本 蓄电池折旧成本。C_real-time(y) 实时购电成本通常更贵 负荷中断惩罚成本 弃光惩罚成本 备用调用成本。约束条件需要分阶段写第一阶段约束仅涉及第一阶段变量x。柴油发电机爬坡约束|P_{dg,t} - P_{dg,t-1}| ≤ Ramp_{dg}蓄电池能量平衡与容量约束SOC_{t1} SOC_t (η_ch * P_ch,t - P_dis,t/η_dis) * ΔtSOC_min ≤ SOC_t ≤ SOC_max日前购电功率上下限P_grid_day-ahead_min ≤ P_{grid,t}^{day-ahead} ≤ P_grid_day-ahead_max第二阶段约束涉及第一阶段变量x、第二阶段变量y和不确定性u。这是鲁棒性的体现要求对于不确定集U内的所有可能情况都必须成立。实时功率平衡约束核心(P_{pv,t}^f ΔP_{pv,t} - P_{pv,t}^{spill}) (P_{wt,t}^f ΔP_{wt,t}) P_{dg,t} r_{dg,t}^ - r_{dg,t}^- P_{grid,t}^{day-ahead} P_{grid,t}^{real-time} (P_{ess,t}^{sch} P_{ess,t}^{adj}) (P_{load,t}^f ΔP_{load,t} - P_{curt,t})这个等式必须对所有ΔP在不确定集内成立。柴油发电机实时出力约束P_{dg,t} r_{dg,t}^ ≤ P_{dg,max} * u_{dg,t},P_{dg,t} - r_{dg,t}^- ≥ P_{dg,min} * u_{dg,t}联络线功率总约束|P_{grid,t}^{day-ahead} P_{grid,t}^{real-time}| ≤ P_grid_max第二阶段变量非负约束等。这个max-min问题就是典型的两阶段鲁棒优化模型。直接求解非常困难我们需要借助算法。4. 算法实现CCG算法在YALMIP中的工程化落地理论上的两阶段鲁棒优化是一个NP难问题。在实践中我们采用列与约束生成算法来求解。CCG算法将原问题分解为一个主问题和一个子问题通过迭代逼近最优解。4.1 CCG算法流程详解初始化设定下界LB -inf上界UB inf迭代次数k0收敛容差ε。生成一个初始的不确定性场景比如所有波动都为0即标称场景将其对应的约束加入到主问题中。求解主问题主问题是在当前已知的“最恶劣场景”集合下寻找最优的第一阶段决策x。随着迭代进行主问题会不断增加来自子问题的“最恶劣场景”及其对应的第二阶段决策约束因此解的质量会越来越好。主问题的解给出当前最优的x_k和一个目标值这个值作为原问题的下界。求解子问题给定主问题求出的x_k子问题要在不确定集U内寻找一个使第二阶段成本最大化的不确定性场景u。即max_{u∈U} min_y C_real-time(y)。子问题是一个max-min问题通常通过对偶理论或KKT条件将其转化为一个单层的最大化问题如果第二阶段问题是线性规划。子问题的解给出在当前x_k下最坏场景u_k及其对应的第二阶段成本η_k。C_day-ahead(x_k) η_k构成了原问题的一个上界。更新与判断更新上下界LB max(LB, 主问题目标值)UB min(UB, C_day-ahead(x_k) η_k)。如果(UB - LB) / LB ≤ ε则算法收敛输出当前解。否则将子问题找到的最坏场景u_k及其对应的第二阶段最优反应即子问题中min_y部分的解y_k作为一个新的场景以约束的形式添加到主问题中。令k k1返回步骤2。4.2 YALMIP中的关键实现技巧在YALMIP中实现CCG核心在于如何动态地构建主问题和子问题。主问题的构建主问题是一个混合整数线性规划。我们需要定义两类变量第一阶段变量x和一系列辅助的第二阶段变量y_l每个l对应一个已发现的最坏场景。每次迭代我们不是修改旧的约束而是新建一个y_l变量并添加与之对应的约束% 在第k次迭代添加新场景约束 y_new sdpvar(num_time, num_second_stage_vars, full); % 定义新的第二阶段变量 Constraints_MP [Constraints_MP, ...]; % 添加针对新场景u_k的功率平衡约束、变量上下限约束等 Constraints_MP [Constraints_MP, PowerBalance(x, y_new, u_k) 0]; Constraints_MP [Constraints_MP, 0 y_new y_max]; % 更新目标函数包含所有场景中最坏的那个成本通过一个辅助变量eta和约束实现 Objective_MP C_day_ahead(x) eta; Constraints_MP [Constraints_MP, eta C_real_time(y_new, u_k)];这里的关键是eta是一个标量变量约束eta C_real_time(y_l, u_l)对于所有已发现的场景l都成立。这样主问题目标C_day_ahead(x) eta自然就会去最小化所有已知场景中最坏的那个总成本。子问题的构建与求解子问题是在x_k固定后关于不确定性u和第二阶段变量y的max-min问题。如果第二阶段问题是线性规划我们可以利用强对偶定理将对偶变量引入将min_y问题转化从而将子问题变成一个单层的最大化问题通常是双线性规划但由于不确定集是盒式的可以进一步处理。在YALMIP中我们可以利用uncertain声明不确定性变量并用robustify命令但更直观的方式是手动实现。 一个更稳定、更通用的方法是将子问题本身也看作一个两阶段问题并用迭代法求解。即内层min_y用线性规划求解给定u外层max_u用启发式算法或者将其转化为混合整数线性规划。在我的实现中由于不确定集是盒式且目标函数对u是线性的子问题的最优解必然在不确定集的顶点上取得。因此我可以将子问题等价地转化为一个混合整数线性规划通过引入辅助整数变量来表征u取正边界还是负边界然后直接用CPLEX求解。这比通用的robustify更高效、更可控。% 假设不确定集为盒式-Δ_max Δ Δ_max % 引入二进制变量 z_t 来表示波动方向 z binvar(3, T, full); % 3种不确定性T个时段 % 将不确定性Δ表示为Δ Δ_max * (2*z - 1) 这样z1时ΔΔ_max, z0时Δ-Δ_max Delta diag(Delta_max) * (2*z - 1); % 构建子问题的目标函数最大化第二阶段成本此时x_k是已知常数 Objective_SP C_real_time(y, x_k, Delta); % y是子问题中的第二阶段变量 % 约束包括第二阶段约束含Delta以及Delta与z的关系 Constraints_SP [SecondStageConstraints(y, x_k, Delta), ...]; % 求解这个MILP得到最坏场景z_opt进而得到u_opt optimize(Constraints_SP, -Objective_SP, ops); % 注意这里是最大化所以加负号4.3 迭代循环与收敛判断将主问题和子问题的求解放入一个while循环中。每次迭代记录主问题的x_k和目标值obj_MP记录子问题的u_k和obj_SP。计算UB min(UB, C_day_ahead(x_k) obj_SP)LB max(LB, obj_MP)。判断相对间隙(UB-LB)/UB是否小于阈值如1e-3。一个重要的工程细节是为主问题和子问题设置不同的求解器参数。主问题是MILP需要设置合适的时间限制和容差避免在早期迭代中花费过多时间求精确解。子问题如果转化为MILP也需要合理设置。同时要保存每次迭代的模型和结果便于调试和绘制收敛曲线。5. 代码实战从零构建MATLAB程序的完整步骤纸上得来终觉浅绝知此事要躬行。下面我结合代码片段讲解关键部分的实现。请注意为了清晰我简化了部分细节但保留了核心逻辑。5.1 环境准备与数据加载首先确保你的MATLAB已安装YALMIP并能成功调用CPLEX求解器yalmiptest命令可以测试。% 清空环境 clear; close all; clc; % 添加路径如果你的CPLEX安装路径特殊 % addpath(genpath(your_cplex_path)); % 测试求解器 yalmiptest(cplex) % 加载微电网数据时间尺度、负荷、光伏/风电预测、电价、设备参数等 load(microgrid_data.mat); T 24; % 24小时调度周期 dt 1; % 1小时为间隔 % 定义设备参数 P_dg_max 500; % kW P_dg_min 100; % kW ramp_dg 200; % kW/h cost_dg 0.5; % $/kWh P_grid_max 1000; % kW price_day_ahead [0.05, 0.05, ...]; % 24小时日前电价 $/kWh price_real_time [0.08, 0.08, ...]; % 24小时实时电价 $/kWh ESS_capacity 1000; % kWh ESS_P_max 250; % kW SOC_min 0.2; SOC_max 0.9; eta_ch 0.95; eta_dis 0.95; % 不确定性参数 P_pv_forecast data.P_pv; % 光伏预测 P_wt_forecast data.P_wt; % 风电预测 P_load_forecast data.P_load; % 负荷预测 delta_max_pv 0.3 * P_pv_forecast; % 最大波动为预测值的30% delta_max_wt 0.4 * P_wt_forecast; delta_max_load 0.1 * P_load_forecast; Gamma 3; % 鲁棒调节参数控制不确定性的保守度5.2 定义优化变量使用YALMIP的sdpvar定义变量。注意区分第一阶段和第二阶段变量。% 第一阶段变量 % 柴油发电机 u_dg binvar(T, 1, full); % 启停状态 P_dg sdpvar(T, 1, full); % 计划出力 % 日前购电 P_grid_da sdpvar(T, 1, full); % 蓄电池计划 P_ess_sch sdpvar(T, 1, full); % 计划充放电功率正为放电 SOC sdpvar(T1, 1, full); % 荷电状态多一个初始时刻 % 将第一阶段变量收集到一个结构体中方便管理 x.P_dg P_dg; x.u_dg u_dg; x.P_grid_da P_grid_da; x.P_ess_sch P_ess_sch; x.SOC SOC; % 用于主问题的第二阶段变量每个场景一份 % 我们将动态添加这里先创建空细胞数组 y_mp_list {}; % 存储每个场景对应的第二阶段变量 u_scenario_list {}; % 存储每个场景对应的不确定性实现值 % 子问题中的变量 % 不确定性变量用于子问题盒式不确定集顶点 z_pv binvar(T, 1, full); z_wt binvar(T, 1, full); z_load binvar(T, 1, full); % 子问题第二阶段变量 P_pv_spill_sp sdpvar(T, 1, full); P_curt_sp sdpvar(T, 1, full); P_grid_rt_sp sdpvar(T, 1, full); r_dg_up_sp sdpvar(T, 1, full); r_dg_down_sp sdpvar(T, 1, full); P_ess_adj_sp sdpvar(T, 1, full); % 将子问题变量也收集到结构体 y_sp.P_pv_spill P_pv_spill_sp; y_sp.P_curt P_curt_sp; y_sp.P_grid_rt P_grid_rt_sp; y_sp.r_dg_up r_dg_up_sp; y_sp.r_dg_down r_dg_down_sp; y_sp.P_ess_adj P_ess_adj_sp; u_sp.z_pv z_pv; u_sp.z_wt z_wt; u_sp.z_load z_load;5.3 构建主问题函数我们将主问题的构建写成一个函数buildMasterProblem它接收当前已发现的场景列表返回主问题的约束、目标函数和变量。这样在迭代中调用更清晰。function [Constraints_MP, Objective_MP, x, eta] buildMasterProblem(scenario_list, system_data) % scenario_list: 细胞数组每个元素是一个结构体包含 u_scenario 和对应的第二阶段成本系数 % system_data: 包含所有设备参数和预测数据的结构体 T system_data.T; % 定义第一阶段变量 (同上略) x defineFirstStageVariables(T); % 定义主问题中的辅助变量 eta (代表最坏情况下的第二阶段成本) eta sdpvar(1,1); Constraints_MP []; Objective_MP system_data.cost_dg * sum(x.P_dg) ... % 柴油机成本 system_data.price_day_ahead * x.P_grid_da * system_data.dt ... % 日前购电成本 0.01 * sum(abs(x.P_ess_sch)) * system_data.dt ... % 蓄电池折旧简化模型 eta; % 最坏情况下的实时调整成本 % 第一阶段约束 % 柴油机约束 Constraints_MP [Constraints_MP, ... system_data.P_dg_min * x.u_dg x.P_dg system_data.P_dg_max * x.u_dg, ... -system_data.ramp_dg diff(x.P_dg) system_data.ramp_dg]; % 蓄电池约束 Constraints_MP [Constraints_MP, ... system_data.SOC_min * system_data.ESS_capacity x.SOC system_data.SOC_max * system_data.ESS_capacity]; for t 1:T Constraints_MP [Constraints_MP, ... x.SOC(t1) x.SOC(t) (system_data.eta_ch * max(0, -x.P_ess_sch(t)) - max(0, x.P_ess_sch(t))/system_data.eta_dis) * system_data.dt]; end % 日前购电约束 Constraints_MP [Constraints_MP, 0 x.P_grid_da system_data.P_grid_max]; % 为每一个已发现的场景添加对应的第二阶段变量和约束 for s 1:length(scenario_list) scenario scenario_list{s}; u_s scenario.u; % 该场景下的不确定性具体值 % 为该场景创建一份独立的第二阶段变量 y_s defineSecondStageVariables_MP(T); % 添加该场景下的功率平衡约束此时x和u_s是已知或变量y_s是变量 % 注意这里的约束是“对于这个特定场景u_s必须存在可行的y_s” [power_balance_con, other_con] buildSecondStageConstraints(x, y_s, u_s, system_data); Constraints_MP [Constraints_MP, power_balance_con, other_con]; % 计算该场景下的第二阶段成本 cost_s calculateSecondStageCost(y_s, u_s, system_data); % 关键约束eta必须大于等于所有已知场景的第二阶段成本 Constraints_MP [Constraints_MP, eta cost_s]; end end5.4 构建子问题函数子问题函数solveSubProblem接收固定的第一阶段决策x_fixed返回最坏场景u_worst和对应的第二阶段成本obj_sp。function [u_worst, obj_sp, status] solveSubProblem(x_fixed, system_data) % x_fixed: 从主问题得到的第一阶段决策数值不是sdpvar T system_data.T; % 定义子问题变量如5.2节所示 [y_sp, u_sp] defineSubProblemVariables(T); % 根据二进制变量z计算实际的不确定性Delta Delta_pv system_data.delta_max_pv .* (2*u_sp.z_pv - 1); Delta_wt system_data.delta_max_wt .* (2*u_sp.z_wt - 1); Delta_load system_data.delta_max_load .* (2*u_sp.z_load - 1); % 构建子问题目标最大化第二阶段成本 Objective_SP calculateSecondStageCost(y_sp, Delta_pv, Delta_wt, Delta_load, system_data); % 构建子问题约束 Constraints_SP []; % 1. 第二阶段运行约束给定x_fixed和不确定性Delta [power_balance_con_sp, other_con_sp] buildSecondStageConstraints_fixedX(x_fixed, y_sp, Delta_pv, Delta_wt, Delta_load, system_data); Constraints_SP [Constraints_SP, power_balance_con_sp, other_con_sp]; % 2. 鲁棒调节参数Gamma约束限制总的不确定性“预算” % 例如限制所有时段所有不确定性源的波动“激活”次数总和不超过Gamma*T % 这是一种常见的处理方式比简单的盒式集更精细。 Constraints_SP [Constraints_SP, ... sum(u_sp.z_pv) sum(u_sp.z_wt) sum(u_sp.z_load) system_data.Gamma * 3 * T]; % 3种不确定性 % 求解子问题这是一个MILP最大化问题 ops_sd sdpsettings(solver, cplex, verbose, 0, cplex.timelimit, 60); diagnostics optimize(Constraints_SP, -Objective_SP, ops_sd); % 注意负号 status diagnostics.problem; if status 0 u_worst.Delta_pv value(Delta_pv); u_worst.Delta_wt value(Delta_wt); u_worst.Delta_load value(Delta_load); obj_sp value(Objective_SP); else warning(子问题求解失败); u_worst []; obj_sp inf; end end5.5 CCG主循环这是整个算法的驱动引擎。% 算法参数 max_iter 20; tol 1e-3; UB inf; LB -inf; iter 0; scenario_list {}; % 存储发现的恶劣场景 history_LB []; history_UB []; history_x {}; % 初始场景标称场景所有波动为0 initial_scenario.u.Delta_pv zeros(T,1); initial_scenario.u.Delta_wt zeros(T,1); initial_scenario.u.Delta_load zeros(T,1); scenario_list{1} initial_scenario; while iter max_iter iter iter 1; fprintf( 第 %d 次迭代 \n, iter); % --- 步骤1求解主问题 --- fprintf(求解主问题...); [Constraints_MP, Objective_MP, x_mp, eta] buildMasterProblem(scenario_list, system_data); ops_mp sdpsettings(solver, cplex, verbose, 0, cplex.timelimit, 100); diagnostics_mp optimize(Constraints_MP, Objective_MP, ops_mp); if diagnostics_mp.problem 0 x_opt getFirstStageSolution(x_mp); % 从x_mp结构体中提取数值解 LB_current value(Objective_MP); LB max(LB, LB_current); fprintf(完成。LB %.2f, UB %.2f, Gap %.2f%%\n, LB, UB, 100*(UB-LB)/UB); else error(主问题求解失败); end % --- 步骤2求解子问题 --- fprintf(求解子问题...); [u_worst, obj_sp, status_sp] solveSubProblem(x_opt, system_data); if status_sp 0 UB_current calculateFirstStageCost(x_opt, system_data) obj_sp; UB min(UB, UB_current); fprintf(完成。最坏场景成本 %.2f, UB %.2f\n, obj_sp, UB); else error(子问题求解失败); end % 记录历史 history_LB(end1) LB; history_UB(end1) UB; history_x{end1} x_opt; % --- 步骤3收敛判断 --- gap (UB - LB) / UB; if gap tol gap 0 fprintf(\n 算法收敛 \n); fprintf(最终解: LB %.2f, UB %.2f, Gap %.4f%%\n, LB, UB, 100*gap); break; end % --- 步骤4添加新场景到主问题 --- new_scenario.u u_worst; scenario_list{end1} new_scenario; fprintf(添加第%d个恶劣场景到主问题。\n\n, length(scenario_list)); end if iter max_iter warning(达到最大迭代次数未完全收敛。); end % 输出最终的第一阶段调度方案 final_solution history_x{end}; disp(第一阶段调度方案柴油机出力:); disp(final_solution.P_dg); disp(第一阶段调度方案日前购电:); disp(final_solution.P_grid_da);5.6 结果可视化与分析算法收敛后绘制收敛曲线和各设备出力图是必不可少的。% 绘制上下界收敛曲线 figure; plot(1:length(history_LB), history_LB, b-o, LineWidth, 1.5, DisplayName, 下界 (LB)); hold on; plot(1:length(history_UB), history_UB, r-s, LineWidth, 1.5, DisplayName, 上界 (UB)); xlabel(迭代次数); ylabel(成本 ($)); title(CCG算法收敛过程); legend(show); grid on; % 绘制最终调度方案 figure; t 1:24; subplot(3,1,1); plot(t, final_solution.P_dg, m-^, LineWidth, 1.5); hold on; plot(t, final_solution.P_grid_da, g--, LineWidth, 1.5); ylabel(功率 (kW)); legend(柴油发电机, 日前购电); title(第一阶段决策柴油机与日前购电计划); grid on; subplot(3,1,2); plot(t, P_load_forecast, k-, LineWidth, 2); hold on; plot(t, P_pv_forecast, y-, LineWidth, 1.5); plot(t, P_wt_forecast, c-, LineWidth, 1.5); ylabel(功率 (kW)); legend(负荷预测, 光伏预测, 风电预测); title(预测值与负荷); grid on; subplot(3,1,3); SOC_plot final_solution.SOC(1:end-1); % 去掉最后一个冗余值 plot(t, SOC_plot, b-, LineWidth, 1.5); ylabel(SOC); xlabel(时间 (h)); title(蓄电池计划荷电状态 (SOC)); grid on;6. 避坑指南与性能优化心得实现过程中我踩过不少坑这里分享几个关键点能帮你节省大量调试时间。6.1 模型不可行如何快速定位问题这是最常见也最令人头疼的问题。当CPLEX返回“infeasible”时不要慌。先检查确定性模型将鲁棒调节参数Gamma设为0运行确定性优化。如果确定性模型就不可行那问题出在基础约束上比如设备容量太小根本无法满足负荷。你需要检查功率平衡等式左右是否匹配单位是否统一上下限设置是否合理例如放电功率上限是否设成了负数。使用CPLEX的冲突refiner在sdpsettings中设置cplex.conflict选项。求解失败后YALMIP/CPLEX可以给出导致不可行的最小约束集这是定位问题的神器。分阶段调试单独构建并求解主问题用初始场景和子问题给定一个合理的x。看是哪个问题不可行。子问题不可行往往意味着你给定的x太“激进”在某些极端场景下无法通过第二阶段调整来弥补这说明你的鲁棒模型可能过于严格或者第二阶段调整能力如备用、负荷中断不足。放松约束尝试暂时注释掉一些非核心约束如爬坡约束、SOC约束看模型是否变得可行。逐步添加约束找到导致不可行的那个。6.2 求解速度慢算法迭代次数多每次求解耗时久两阶段鲁棒优化的计算复杂度很高性能优化至关重要。削减约束在将子问题场景加入主问题时不是所有约束都需要严格对偶化。对于某些简单的上下界约束可以直接用原问题的约束这能显著减少主问题的规模。设置合适的求解器参数对于主问题MILP设置cplex.mip.tolerances.mipgap为一个稍大的值如1e-2在迭代前期快速得到一个可行解即可不必追求绝对最优。在迭代后期再缩小这个间隙。利用热启动CPLEX支持热启动。在迭代求解主问题时可以将上一次迭代的解作为初始解传入能加速求解。不确定集的选择盒式不确定集最直观但可能过于保守。Budget不确定集用Gamma控制总波动量更符合实际也能减少极端场景的数量从而可能加快收敛。我的代码中已经体现了这一点。并行计算如果子问题可以独立求解例如不同时段耦合不强可以考虑并行求解多个场景的子问题。但在标准的CCG中子问题是串行的。6.3 结果不鲁棒调度方案在实际波动下表现不佳这可能是因为你的不确定集没有很好地反映真实的不确定性。校准不确定集delta_max和Gamma不是随便设的。应该基于历史预测误差数据进行分析。例如delta_max可以取历史误差的某个分位数如95%分位数。Gamma可以根据你对风险的态度来调整。进行蒙特卡洛仿真验证得到鲁棒优化方案后不要只看最坏场景。应该用大量随机生成的可能场景符合历史误差分布去测试这个方案统计其约束违反概率和经济成本分布。这才是评估方案鲁棒性的金标准。考虑相关性光伏和风电出力可能有相关性负荷波动也可能有模式。简单的盒式或Budget不确定集忽略了相关性可能导致结果保守。更高级的集合如基于数据驱动的多面体集合可以考虑相关性但建模和求解会更复杂。6.4 YALMIP使用技巧避免在循环中重复定义变量像我在主问题函数中那样将变量定义放在函数外或初始化部分。在循环内只添加约束。使用value()函数提取解求解后用value(P_dg)来获取变量的数值。对于结构体中的变量需要逐层提取。清理旧模型每次迭代构建新主问题时旧的约束和变量如果不再需要确保MATLAB工作空间不会积累太多对象以免内存溢出。可以通过在函数内定义变量、函数结束时自动清除来实现。保存和加载模型对于复杂模型使用save和load命令保存中间结果和模型便于断点调试。实现一个完整的微电网两阶段鲁棒优化模型就像完成一个精密的机械组装。你需要深刻理解数学原理熟练运用建模工具并耐心处理工程实现中的各种细节。这个过程充满挑战但当看到算法收敛并生成一个能抵御风险的经济调度方案时那种成就感是无与伦比的。希望我的这份“踩坑实录”和完整代码思路能为你点亮前行的路。本文还有配套的精品资源点击获取