☰
数据中心微网两阶段鲁棒优化:CCG算法与Matlab复现全记录
2026/10/11 3:07:46 网站建设 项目流程

最近在复现一篇EI论文里的数据中心微网规划方法,核心是两阶段鲁棒优化,模型里还把数据中心IT负载、储能、制冷的灵活性资源一起考虑进去了。前前后后折腾了两周,从数学建模到CC&G迭代求解,最后把Matlab代码完整跑通,这期间踩了不少坑,也把关键参数和调试经验都沉淀了下来。这篇博文就是这次复现的完整记录,从数学模型到列与约束生成算法,再到算例设计和排错心得,我会尽量把每一步讲透。如果你正在做微电网规划、综合能源系统优化或者鲁棒优化相关课题,又或者你只是想知道两阶段鲁棒规划在Matlab里到底怎么落地,这篇文章应该能帮你省下不少摸索时间。

1. 先搞懂模型在做什么:灵活性加鲁棒规划的核心逻辑

1.1 数据中心微网规划:为什么不能套普通微网的思路

数据中心和普通园区微网最大的区别在于它的负荷结构。一个数据中心的用电量里,IT设备占四到六成,制冷系统占三到四成,这两块不是简单并联关系,而是互相耦合的:IT负载越高,设备发热越大,制冷需要消耗的能量也越多。传统微网规划里把负荷当成一个固定曲线或者仅按季节性变化来处理,数据中心这种强耦合、可变性强的场景根本没法反应真实运行情况。

另一个特点是可靠性与经济性之间的张力。数据中心通常要求供电可用率极高,所以规划时不能只看平均场景下能不能经济运行,还得考虑光伏出力偏差较大、负荷突增这些恶劣情形会不会导致供电不足。这就意味着规划模型必须在投资阶段就给系统留有足够的"余地",而这个余地的量化方法,恰恰是鲁棒优化的看家本领。

我在复现时的第一个体会是:不要直接拿论文公式就开始敲代码,先画一张系统的能量流图,把光伏、储能、燃气轮机、主网交互、数据中心负载、制冷负载之间的关系理清楚。这张图画清楚了,后面的目标函数和约束项基本就水到渠成了。

1.2 "灵活性"在模型里究竟怎么表达

这个词在论文里出现频率很高,落到数学上其实是一组可调节的决策变量和相应的时段耦合约束。数据中心微网里的灵活性来源主要有三个:

第一个是IT负载的可转移性。数据中心里有一部分计算任务并不是必须即时完成的,比如夜间的大规模数据备份、非交互式的训练任务、数据分析批作业,它们可以在一个时间窗口内挪动。反映到模型里,就是每个时段的数据中心总功耗约束里,增加一个可平移的负载量,同时保证整个调度周期内被转移的总工作量不变。

第二个是储能和UPS的协同。数据中心的UPS电池本身是一个稀缺的灵活性资源,规划阶段往往只考虑它的容量需求,但实际上在鲁棒运行模型中,UPS和储能电池一起参与功率平衡,可以在电价高峰或光伏不足时放电,给主网交互功率"削峰"。这个协同机制如果不在模型里体现,储能容量往往会配置偏大。

第三个是制冷系统的热惯性。冷冻水系统本身有一定的蓄冷能力,室温在一定范围内可以波动,所以制冷功率不必时刻匹配IT散热负荷,而是可以在时间维度上有一定的弹性。论文模型一般把这个热惯性简化成一个制冷功率的调节范围加能量平衡约束,不需要非常复杂的建筑热动态模型。

我在代码实现中,把这三个灵活性都反映在第二阶段运行约束里,后面对比算例时也单独测试了“去掉可转移负载约束”的情况,储能规划容量直接上升了约两成。这说明灵活性资源不是摆设,它对投资决策的影响是实打实的。

1.3 为什么是两阶段鲁棒而不是随机规划

做规划问题,最直觉的思路是用确定性优化:给一组典型的负荷曲线、光伏出力曲线,优化出设备容量。但光伏出力的随机性很强,一个夏天午后突然飘来的云就能让出力在几分钟内掉四成,确定性方案的容量配置在这种场景下就是裸奔。

随机规划是另一种常用思路,它需要给不确定参数指定概率分布,然后抽样生成大量场景。但在工程实际里,光伏出力预测的误差分布很难准确刻画,场景数量一多计算量也大,场景数量少了又容易低估风险。

鲁棒优化走的是一条更"保守"的路:你不需要知道不确定参数的具体概率,只需要给出它的变化范围,也就是不确定集,然后优化在最坏情况下的表现。这样得到的方案天然具备抗风险能力,而且当不确定集参数设计合理时,不会过度牺牲经济性。

两阶段结构则对应了规划与运行的自然层级。第一阶段是投资决策——装多少光伏、配多大储能,属于"现在就要定下来"的变量;第二阶段是运行调度——每个时段的购售电功率、储能充放电、可转移负载安排,属于"等实际场景发生后再调整"的变量。两阶段鲁棒优化把这两层决策耦合起来,第一阶段决策要能应对第二阶段可能出现的所有恶劣场景,这个逻辑天然匹配微网规划问题,所以我看到论文标题时就觉得建模框架选得很准确。

2. 数学模型:目标函数、约束与不确定集的逐条拆解

2.1 第一阶段规划决策变量与投资成本

第一阶段的决策变量通常包含光伏装机容量 (C_{pv})、储能额定容量 (C_{ess})、储能额定功率 (P_{ess}^{rate}),有些模型还会有燃气轮机容量 (C_{gt})。这些是规划层面的选择,一旦确定在整个调度周期内不变。

目标函数的第一项是年化投资成本,形式一般是:

[ C_{inv} = \alpha_{pv} C_{pv} + \alpha_{ess,c} C_{ess} + \alpha_{ess,p} P_{ess}^{rate} + \alpha_{gt} C_{gt} ]

这里的 (\alpha) 是各设备的年化投资系数,由单位投资成本乘以年折现因子得到。我建议注意一下论文里给的投资成本单位是元/kW还是万元/MW,复现时最容易在这里差几个数量级。

第一阶段的约束相对简单,主要是各设备的规划容量上限:

[ 0 \le C_{pv} \le C_{pv}^{max}, \quad 0 \le C_{ess} \le C_{ess}^{max}, \quad 0 \le P_{ess}^{rate} \le P_{ess}^{rate,max} ]

有些论文还会加入储能功率与容量的匹配约束,比如容量不能小于功率乘以最低持续充放电时长,这个约束也是投资决策里很常见的。

2.2 第二阶段运行调度约束:数据中心特征的核心体现

第二阶段是在给定设备容量和某个不确定场景后,求解每个时段的运行调度问题。决策变量包括主网购电功率 (P_{grid,t})、售电功率 (P_{sell,t})、燃气轮机出力 (P_{gt,t})、储能充电功率 (P_{ch,t})、放电功率 (P_{dis,t}) 以及可转移负载量 (\Delta_t)。

数据中心微网和普通微网的最大区别在这里以约束的形式体现出来。首先是数据中心负载平衡约束,可以写成以下形式:

[ P_{pv,t} + P_{gt,t} + P_{grid,t} - P_{sell,t} + P_{dis,t} - P_{ch,t} = P_{base,t} + P_{cool,t} + \Delta_t ]

其中 (P_{base,t}) 是IT设备的基础功耗,(P_{cool,t}) 是制冷功耗。注意这里我采用了单位时段能量平衡的形式,实际代码里每个时段都是一条这样的平衡方程。

制冷功率的约束可以写成范围加时段间爬坡的形式:

[ P_{cool,min} \le P_{cool,t} \le P_{cool,max} ]

[ -P_{cool}^{ramp} \le P_{cool,t} - P_{cool,t-1} \le P_{cool}^{ramp} ]

爬坡约束体现的就是热惯性——制冷功率不能瞬间大幅跳变,因为温度变化有滞后。

储能约束则是经典的SOC递推:

[ SOC_{t+1} = SOC_t + \eta_{ch} P_{ch,t} - \frac{P_{dis,t}}{\eta_{dis}} ]

[ SOC_{min} \le SOC_t \le SOC_{max}, \quad SOC_1 = SOC_T ]

可转移负载的约束要重点说。论文里通常的处理方式是设定转移比例上限和周期转移量守恒:

[ |\Delta_t| \le \delta_{max} \cdot P_{base,t} ]

[ \sum_{t=1}^{T} \Delta_t = 0 ]

第一条是说每个时段允许转移的负载量有上限,防止模型把所有负载都堆到某个时段;第二条是周期守恒,即数据中心一段时期内实际执行的计算总量不变,只是时间上挪动了。这个约束一旦写错,结果可能出现"模型通过无限转移负载来逃避购电成本"的荒谬现象,我在调试早期就撞上过。

2.3 不确定集设计:保守度是怎么被旋钮控制的

鲁棒优化的核心是不确定集 (U)。我复现的论文里采用了盒式加预算的经典组合。光伏出力和数据中心基础负荷分别引入不确定性。

光伏出力的不确定集可以写成:

[ P_{pv,t} \in [P_{pv,t}^{f} - \hat{P}{pv,t}, P{pv,t}^{f} + \hat{P}_{pv,t}] ]

其中 (P_{pv,t}^{f}) 是预测值,(\hat{P}_{pv,t}) 是允许的偏差范围,一般用预测值乘以一个百分比。负荷不确定性同理。

单独用盒式不确定集的问题在于它允许"所有时段同时达到最坏情况",这在实际中几乎不会发生,导致结果非常保守。所以引入了预算约束:

[ \sum_{t=1}^{T} \frac{|P_{pv,t} - P_{pv,t}^{f}|}{\hat{P}_{pv,t}} \le \Gamma ]

这个 (\Gamma) 就是控制保守度的关键参数。它的含义是:在调度周期内,最多有多少个时段的不确定参数同时取到边界值。如果 (\Gamma=0),模型退化为确定性优化;如果 (\Gamma=T),就变成了完全盒式鲁棒。实际使用中可以通过扫描 (\Gamma) 来画出"成本-保守度"的帕累托曲线,帮决策者选择合适的配置方案。

这里的实现细节是,我将不确定集以显式的 (\pm) 偏差形式直接写进第二阶段约束,没有用场景枚举。这样在CC&G框架下,每个迭代步只需要处理一个被不断收紧的"不确定参数到约束右侧"的线性问题,代码会简洁很多。

3. CC&G求解算法与Matlab代码实现:从原理到可运行

3.1 列与约束生成(CC&G)的迭代逻辑

两阶段鲁棒问题的一般形式可以写成:

[ \min_{x} \left[ c^T x + \max_{u \in U} \min_{y \in F(x,u)} d^T y \right] ]

CC&G算法的核心思想是把外层min和内层max-min问题交替求解。具体来说,主问题是在当前已知的若干"最坏场景"下做联合优化,子问题则是给定第一阶段决策后,找出当前最坏的场景和对应的运行成本,并据此生成新的约束(割)反馈给主问题。

迭代过程大体如下:

  • 初始化一个最坏场景 (u_1)(通常取不确定集的一个极值点)
  • 迭代开始,求解主问题,得到第一阶段解和当前下界 (LB)
  • 将第一阶段解代入子问题,求解最坏场景和对应运行成本,得到当前上界 (UB)
  • 把子问题求解出来的最坏场景加入主问题的场景集合,生成割约束
  • 重复直到上下界之差满足收敛精度

这个算法在工程中的表现比Benders分解稳定不少,因为割是直接针对场景生成的,而不是针对对偶约束生成的,所以迭代次数通常更少,收敛性也更好。我在实际复现中,24时段、3类不确定参数的算例大概迭代8到12次就能收敛,相比Benders动辄几十次已经快很多了。

3.2 主问题MP的YALMIP表述

主问题里除了投资决策变量,还有每个已有场景对应的运行变量。随着迭代进行,场景数越来越多,运行变量组也会越来越多,这是YALMIP建模时需要注意的。

主问题的结构可以抽象为:

[ \min \quad c^T x + \eta ]

[ s.t. \quad Ax \le b ]

[ \eta \ge d^T y_k, \quad \forall k \in {1,...,K} ]

[ By_k + Cx \le g + D u_k, \quad \forall k ]

这里的 (k) 是已生成的场景索引,(u_k) 是第k次子问题算出的最坏场景。(\eta) 是一个辅助变量,用来近似表示最坏运行成本。每一次迭代多一组 (y_k) 变量和对应的约束,这是CC&G"列"与"约束"同时生成的含义。

3.3 子问题SP的对偶变换与求解

子问题是给定第一阶段 (x = \bar{x}) 后求解:

[ \max_{u \in U} \min_{y} d^T y ]

[ s.t. \quad By \le g + D u - C \bar{x} ]

内层min是一个线性规划,所以可以用强对偶定理转换成max问题。这样整体的max-min问题就变成了一个单层的max问题。关键来了,对偶变换后,约束右边出现了 (u) 与对偶变量的乘积项,这就产生了双线性项。

处理这个双线性项,我在复现中采用了两种方案,取决于 (u) 的取值特性。如果不确定集是离散的(比如简化为若干种典型偏差场景),直接枚举即可。如果是连续的,那就需要用大M法引入辅助变量和0-1变量来线性化。我实际采用的是后者,因为论文里的不确定集是连续盒式加预算约束。

经过对偶变换后,子问题变成了一个含有0-1变量的混合整数线性规划MILP,交给Cplex求解器计算。这里要特别注意对偶变量与大M参数的设置,后面在踩坑部分我再详细说。

3.4 Matlab代码框架:CC&G主循环的骨架

下面是我整理的CC&G主循环Matlab代码框架,基于YALMIP建模、Cplex求解:

%% 参数初始化 T = 24; % 调度周期时段数 max_iter = 30; % 最大迭代次数 tol = 1e-4; % 收敛精度 LB = -1e8; UB = 1e8; iter = 0; U_scenarios = {}; % 存最坏场景 U_scenarios{1} = initial_scenario; % 初始场景 %% 定义主问题变量(一次定义,迭代中复用) x = sdpvar(1, 3); % 第一阶段:光伏、储能容量、储能功率 eta = sdpvar(1, 1); % 辅助变量 % 注意:y_k 变量在循环内动态添加 %% CC&G主循环 while ((UB - LB) / abs(UB) > tol) && (iter < max_iter) iter = iter + 1; fprintf('---- 迭代 %d ----\n', iter); % ========== 步骤1:求解主问题 ========== mp_constr = []; mp_cost = alpha * x' + eta; % 第一阶段基本约束 mp_constr = mp_constr + [lb_x <= x <= ub_x]; % 对每个已生成场景,加入运行约束 for k = 1:length(U_scenarios) u_k = U_scenarios{k}; % 定义该场景对应的第二阶段变量 y_k = sdpvar(ny, T); % 添加运行约束 By + C*x <= g + D*u_k mp_constr = mp_constr + [runtime_constr(y_k, x, u_k)]; % 添加成本约束 eta >= d'*y_k mp_constr = mp_constr + [eta >= cost_coef * y_k(:)]; end ops = sdpsettings('solver', 'cplex', 'verbose', 0); optimize(mp_constr, mp_cost, ops); LB = value(mp_cost); x_best = value(x); % ========== 步骤2:求解子问题 ========== % 固定 x,求最坏场景下的运行成本 [worst_cost, u_worst] = solve_subproblem(x_best); UB = min(UB, alpha * x_best' + worst_cost); % ========== 步骤3:收敛判断与场景添加 ========== fprintf('LB = %.2f, UB = %.2f, gap = %.4f\n', LB, UB, (UB-LB)/abs(UB)); if (UB - LB) / abs(UB) <= tol break; end % 把新场景加入集合 U_scenarios{end+1} = u_worst; end

solve_subproblem函数内部实现的是对偶转换后的MILP求解。这里有一个很重要的实现细节:每轮迭代的主问题里,新加入的场景对应的变量组是全新的,但已有的场景变量组不需要重新优化吗?答案是需要的,因为第一阶段决策 (x) 变了,所有场景下的运行变量都在同一轮优化里重新求解。所以我把mp_constr和mp_cost放在同一个针对 (x) 和所有 (y_k) 的优化问题里,而不是只对 (\eta) 做增量更新。

4. 算例设计与结果解读:参数整定与数据可视化实战

4.1 系统参数与基准场景设置

算例参数的设置直接决定复现结果能否和论文对上。我采用的系统参数是参考目标论文算例并做了适当简化调整,主要数据如下:

参数取值说明
数据中心IT基础负荷峰值800 kW含可转移负载潜力
可转移负载比例上限20%相对基础负荷
光伏规划容量上限1200 kW可用屋顶与场地限制
储能规划容量上限600 kWh—
储能充放电功率上限200 kW—
主网交互功率上限500 kW变压器容量限制
光伏预测偏差±15%相对预测出力
负荷预测偏差±10%相对基线负荷
预算不确定集参数 Γ624时段中最多6个偏差

负荷曲线我采用了典型的数据中心日负荷形态,凌晨三到五点有一个低谷,白天逐渐升高,晚上八点左右到达峰值。光伏出力曲线则按晴朗夏季典型日生成,正午达到峰值。预测曲线作为名义值输入,偏差在一定范围内波动。

4.2 不确定集预算值对配置结果的影响

预算值 (\Gamma) 是最值得扫描的参数。我分别计算了 (\Gamma = 0)(确定性)、(\Gamma = 6)(中等保守)、(\Gamma = 12)(高度保守)三种情况下的规划结果,数据如下:

(\Gamma)光伏容量(kW)储能容量(kWh)年化总成本(万元)迭代次数
0720280286.51
6850430312.89
12920560344.211

这个趋势非常符合直觉:不确定集越大,系统需要配置更多的光伏和储能来应对最坏情况,投资成本和预期运行成本自然水涨船高。但注意光伏容量的增长幅度明显小于储能容量的增长,这说明在应对短期功率波动时,储能是更有效的灵活性手段,光伏本身解决不了"夜里没太阳"的时段匹配问题。

从工程决策的角度看,(\Gamma = 12) 意味着假设一天中有一半天数都在最坏情况下运行,这显然过于保守。而 (\Gamma = 6) 对应的方案在成本上只比确定性方案高约9%,但能在光伏波动时保证供电可靠性,算是性价比比较高的选择。

4.3 灵活性参数变化对容量配置的敏感性

这个部分的仿真结果是我个人认为整个项目最有价值的信息。我固定 (\Gamma = 6),调整数据中心可转移负载比例上限,从0逐步提高到30%,观察储能配置的变化:

可转移负载比例储能容量(kWh)总成本(万元)储能利用率
0%510326.768%
10%465318.972%
20%430312.878%
30%405308.681%

可见,只要把10%的IT负载变成可调度的灵活性资源,储能需求就下降了约9%;到30%可转移比例时,储能容量相比无灵活性方案减少了约20%,总成本下降了约5.5%。这组数据说明,数据中心微网的规划不能只盯着储能和光伏的容量选择,"需求侧的灵活性"同样是影响投资决策的重要杠杆。

这也让我理解了论文为什么要花大篇幅把数据中心负载特征建模进去——因为数据中心的IT负载天然具有可转移性,这是它区别于普通工业园区的独特优势,把这种灵活性资产纳入规划,能直接反映为投资成本的下降。

5. 复现避坑指南:双线性项、数值不稳定与结果对齐

5.1 子问题双线性项处理的三个陷阱

子问题对偶变换后出现的 (u \cdot \lambda) 双线性项,是整个复现过程中我最头疼的部分。这里有三个具体问题值得展开。

第一个陷阱是把大M值设得过大。线性化双线性项时,通常引入辅助变量 (v = u \cdot \lambda),并添加约束:

[ v \le M \cdot \lambda, \quad v \ge -M \cdot \lambda ]

[ v \le u + M(1 - \lambda), \quad v \ge u - M(1 - \lambda) ]

如果M取得太大,比如一上来就设1e6,Cplex求解时的数值稳定性会严重恶化,会出现主问题可行但子问题退化解的情况,上下界振荡。我的做法是先把每个约束的量级估算出来,然后取最优预算范围内最大可能值的1.2倍作为M值,并测试一下结果对M值的敏感性。

第二个陷阱是忽略了 (u) 与对偶变量 (\lambda) 的取值范围差异。光伏偏差 (u) 是归一化在[-1,1]区间的,但对偶变量的量级取决于运行成本系数和约束右侧的量级,可能到10的3次方以上。如果对两个变量用同一个M值,M就会被取到一个很大值来覆盖对偶变量范围,又回到了第一个陷阱。所以分开设定M值效果要好得多。

第三个陷阱是不确定集预算约束在对偶后如何保留。预算约束本身是一堆绝对值不等式,对偶化之后绝对值仍然会带来非线性的额外项。我在复现中的做法是先把不确定集写成显式的参数化形式,用一组非负的辅助变量来表示偏移量,避免绝对值运算出现在对偶问题里。

5.2 求解慢与不收敛的调试思路

如果你遇到了CC&G迭代不收敛,不要先去翻算法文献,先检查以下三件事。

第一,检查子问题是否真的求到了全局最优。Cplex求解MILP本身能保证全局最优,但如果您把内层min问题手动对偶转换时写错了对偶变量的正负号,那子问题的"最优"就是错的。这个检查办法很简单,构造一个固定场景,直接把内层min用LP求解一次,再把对偶后的max求解一次,对比两个目标函数值是否一致。

第二,检查主问题的割约束是否在不断重复已经存在的场景。如果新生成的场景与之前的场景非常接近,主问题会退化为同一个结构,导致下界提升缓慢。这时候可以在子问题中加入一个小扰动项,或者设置一个"场景去重"机制,保证每次加入的场景在欧几里得距离上显著区别于已有场景。

第三,检查收敛判据。我推荐用相对缺口 ((UB - LB) / \max(1, |UB|)) 而不是绝对差值,因为投资成本和运行成本的量级差异很大,绝对差值容易导致过早终止或永不终止。

在调试一个典型的24时段问题时,我遇到的另一个情况是子问题求解时间随迭代次数线性增长,因为主问题的变量数量在增加。这里可以做的优化是只让每个新增场景对应的变量组参与求解,并以热启动方式把上一轮主问题解作为当前轮初始解,实测能减少约30%的求解时间。

5.3 与论文结果对不上的排查清单

复现论文时最尴尬的事情莫过于代码跑通了,结果和论文对不上。我把排查顺序整理成了一份清单,按优先级排列。

先查量纲和归一化。论文里可能出现的单位组合非常多,功率用kW、能量用kWh、投资成本用万元、运行成本用元。我遇到过光伏容量被系统自动换算成MW,导致投资成本比论文高三个数量级的情况。最佳实践是在代码里统一使用一套单位并加注释,关键时刻多打印几个关键数值来验证。

再查时间粒度。有些论文的优化周期是24小时,但运行约束是按4小时一个时段聚合的,也就是6个时段。如果代码直接用了24个时段但是参数还是按聚合时段给的,相当于每个时段的允许充放电能量被缩小了。这类问题通常表现为储能利用率异常低或可转移负载数值不合理。

最后查目标函数里的常数项。年化投资成本的计算涉及折旧年限和折现率,不同论文假设差别不小。如果你的投资成本项和论文差异在20%以内,大概率是折现参数选取不同,这个属于合理差异。但如果差异超过50%,就要回看公式里的年化系数是不是算错了。

最后的实操体会

复现做完之后,我自己最大的收获其实不是跑通了CC&G算法,而是理解了灵活性的价值不能只靠定性描述,它需要通过精确的约束建模才能在规划决策中释放出来。数据中心的IT负载可转移比例每提高10%,储能需求能降低约一成,这个量化结论对实际项目的前期设计很有指导意义。

如果你准备动手做类似的复现,我建议先把简化算例跑通,比如只考虑光伏不确定性和储能两个灵活性资源,把CC&G的迭代逻辑和子问题对偶彻底吃透,再逐步加入制冷耦合、可转移负载、燃气轮机这些复杂元素。另外,Matlab和YALMIP的组合在快速原型验证上效率确实高,但如果后续要把模型做到更大规模,可以考虑转到Python的Pyomo或Julia的JuMP框架,它们在处理大规模线性规划和复杂割结构时更灵活。

最后补一个小的代码技巧:在所有约束和成本函数里统一使用列向量,避免YALMIP自动广播带来的维度混乱;每次迭代开始前把上一次求解的x_best作为变量初值传给下一次求解,能显著减少MILP的求解时间。这些小习惯积累起来,复现效率会有明显提升。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询