做含分布式电源的配电网优化研究，最难的不是把目标函数写出来，而是怎么让你写的这个优化问题真能被稳定解出来。最近我们团队在做一个园区微电网的调度方案验证，要求同时处理分布式电源接入、储能充放电、联络线功率调节三件事，目标是最小化运行成本和网损。一开始我用的是传统非线性最优潮流，CPLEX和Matlab自带求解器换了好几个，要么算不出来，要么算出来不收敛。后来把模型改成二阶锥松弛形式，用YALMIP建模配合MOSEK求解，问题直接变成凸优化，不仅全局最优性有保证，求解速度也快了几个量级。这篇就把整个思路掰开揉碎讲一遍，从为什么需要微电网灵活性、DistFlow方程怎么松弛成二阶锥，到Matlab代码怎么组织、算例结果怎么看，以及松弛紧度验证的几个坑，一次性说清楚。
1. 为什么含分布式电源的配电网要考虑“微电网灵活性”
1.1 高渗透DG给配电网带来的三个结构性矛盾
先说项目背景。我用的算例是改造过的IEEE 33节点配电系统，基准电压12.66 kV，总负荷约3715 kW加2300 kvar。接入分布式电源之后，立刻暴露三个以前纯辐射状配电网研究里不太会遇到的矛盾。
第一个是出力反调峰。光伏在午间大发，但午间恰恰是负荷低谷，净负荷曲线被压得很低甚至出现负值，也就是功率倒送。风机出力则更多集中在夜间，跟负荷高峰错位更严重。这种时间尺度上的错位，导致单纯的断面优化很难给出真正可执行的调度方案。
第二个是电压双向越限。传统配电网设计假设潮流从变电站单向流向负荷端，有载调压变压器的抽头和无功补偿也都基于这个假设配置。DG接入后，部分馈线末端电压被抬高，尤其在轻载高发时段，末端电压可能超过1.05 p.u.，而重载时段首端电压又可能低于0.95 p.u.。规划人员和调度人员第一次看到电压曲线时通常会不太适应——同一个馈线上，不同节点电压变化方向竟然是相反的。
第三个是爬坡压力。光伏出力受云层影响，20分钟内可以从满发跌到20%。这个功率缺额如果全部转嫁给上级电网联络线，主网的调频和备用压力会显著增加。换句话说，配电网本身的调节能力越差，它对主网灵活性的索取就越多。
这三个矛盾叠加在一起，单看某一个断面的最优潮流已经不够了。调度关心的不再只是“当前这一个时刻怎么调”，而是“遇到外部扰动时系统还剩多少调节空间”。这就是标题里“微电网灵活性”要回答的问题。
1.2 “灵活性”不是一句口号，而是可量化约束
很多做电力系统研究的同行看到“灵活性”三个字会觉得虚，觉得它不像网损、电压偏移那样有一个明确的物理单位。我一开始也有这个困惑。实际把项目需求翻译成数学模型之后，我认为可以这样理解：灵活性就是系统在发电或负荷发生变化时，仍然能通过调节DG、储能、柔性负荷等资源，维持安全运行的能力。
落到最优潮流模型里，我把它拆成了三个可以计算的指标，用在项目中后期和调度沟通非常顺：
| 灵活性维度 | 物理含义 | 在OPF中的量化方式 |
|---|---|---|
| 联络线爬坡余量 | 并网点交换功率还能承受多大的波动 | 当前Pcc与上下限的距离，以及相邻时段ΔPcc是否留有余量 |
| 可调资源潜力 | 所有DG和储能还能增发/少发多少功率 | 各DG剩余可调容量之和，储能当前SOC与上下限的距离 |
| 电压调节裕度 | 电压离越限边界还有多远 | 全局最低电压与下限的差、最高电压与上限的差 |
这三个指标最终以两类形式进入优化模型。一类是约束形式，比如要求联络线相邻时段功率变化量不超过给定爬坡限值，要求储能SOC保持在0.2到0.9之间而不是用到极值。另一类是目标函数里的惩罚项，比如当解逼近某个调节边界时增加惩罚成本，让优化器主动避开“看着经济但完全没有裕度”的运行点。
这里我想强调一点：如果你在学术论文里做研究，可以继续把它定义为灵活域半径、灵活性鲁棒可行域之类的严格数学指标。但工程上先把这三个可计算指标做出来，用来指导潮流计算和调度决策，反而是推进项目最快的方式。后面要扩展成严格指标，也只是在这些量化结果基础上再包装一层数学定义。
2. DistFlow方程与二阶锥松弛：从非凸到可解的建模路线
2.1 辐射状配电网的DistFlow方程
配电网不像输电网那样有很密的环网结构，绝大多数是辐射状运行。对辐射状配电网做潮流计算，业界普遍用DistFlow方程，它比传统极坐标牛顿法更适合配电网的弱环和辐射状特性，也更容易嵌入优化模型。
定义每个节点i的电压幅值为V_i，每条支路ij的有功潮流为P_ij，无功潮流为Q_ij，电流幅值为I_ij。DistFlow方程的核心形式如下：
节点有功平衡：
P_Gi - P_Li = ∑_{j: i→j} P_ij - ∑_{l: l→i} (P_li - r_li * I_li²)
节点无功平衡：
Q_Gi - Q_Li = ∑_{j: i→j} Q_ij - ∑_{l: l→i} (Q_li - x_li * I_li²)
支路电压降落方程：
V_j² = V_i² - 2(r_ij P_ij + x_ij Q_ij) + (r_ij² + x_ij²) * I_ij²
视在功率与电流电压的关系：
I_ij² = (P_ij² + Q_ij²) / V_i²
这里面，最后一个方程是问题的根源。P_ij、Q_ij、V_i都是决策变量，它们之间是二次关系，而且是等号约束，整体可行域是非凸的。带这种约束的最优潮流，即便系统规模不大，求解器也经常陷入局部最优甚至不收敛。
2.2 变量替换与二阶锥松弛一步
Baran和Wu在1989年提出的处理方式今天仍然是标准做法：先做变量替换，把电压和电流的高次项压平。
令 v_i = V_i²，l_ij = I_ij²，带入上面的方程后，DistFlow方程组变成更规整的形式：
v_j = v_i - 2(r_ij P_ij + x_ij Q_ij) + (r_ij² + x_ij²) * l_ij
节点平衡方程也变成：
P_Gi - P_Li = ∑ P_ij - ∑ (P_li - r_li * l_li) Q_Gi - Q_Li = ∑ Q_ij - ∑ (Q_li - x_li * l_li)
到这里仍然是等号约束的非凸问题，最关键的一步在最后那个方程上。我们把等号约束
l_ij = (P_ij² + Q_ij²) / v_i
松弛成不等式：
l_ij ≥ (P_ij² + Q_ij²) / v_i
这一步在数学上等价于下面的二阶锥约束：
‖ [2P_ij, 2Q_ij, v_i - l_ij]ᴛ ‖₂ ≤ v_i + l_ij
写成这个形式之后的好处非常明显：约束变成了凸锥约束，整个最优潮流问题变成一个标准二阶锥规划（SOCP）问题。凸优化理论保证求出来的解是全局最优解，不再受初始点选择影响，也不存在陷入局部最优的问题。而且SOCP问题有非常成熟的内点法求解器，比如MOSEK、SeDuMi、SDPT3，求解规模可以到几千个节点。
在Matlab里用YALMIP建模时，这个锥约束可以很自然地写成：
for k = 1:nBranch Constraints = [Constraints, cone([2*P_br(k,t); ... 2*Q_br(k,t); ... v(from(k),t) - I2_br(k,t)], ... v(from(k),t) + I2_br(k,t))]; endYALMIP的cone(x, y)函数表示 ‖x‖₂ ≤ y，正好对应SOCP的标准形式。
2.3 松弛之后求出的解还能不能当真？紧度问题
很多人第一次接触松弛法会问同样的问题：你把等号改成不等号，求出来的解还是原问题的解吗？这个问题必须正面回答，答案就是“需要验证紧度”。
所谓紧度（tightness），就是松弛后求出的最优解，原始非线性等号约束是否恰好以等号成立。如果每个支路在最优解处都满足
l_ij = (P_ij² + Q_ij²) / v_i
说明松弛没有把解“放松”到一个原问题不可行的虚假点上，这个解可以直接还原为原问题的可行潮流解。数学上有很多研究证明，对辐射状配电网，在目标函数是网损、购电成本等常见形式时，SOCP松弛通常是紧的，理论上比较可靠。
但工程实现上不能只靠理论保证，我习惯在代码里显式检查。仿真结束后计算每个支路的间隙：
gap_ij = |l_ij - (P_ij² + Q_ij²) / v_i|
当所有gap都在1e-5以下时，基本可以认为松弛是紧的，结果可信。这个检查代码很短，务必保留在最终脚本里，后面专门讲。
3. 灵活性资源的数学刻画：DG、储能与联络线功率的协同
3.1 逆变器型DG的可调域
分布式电源按控制方式可以分为两类。一类是不可控的，比如按MPPT最大功率跟踪运行的光伏逆变器，有功出力由光照决定，模型里直接当成负的负荷处理。另一类是可调的，比如微型燃气轮机、可参与调节的储能变流器，有功和无功出力可以在一定范围内自由选择。
对于可调的逆变器型DG，需要同时限制有功、无功和视在功率。如果没有明确的旋转设备（如燃气轮机），逆变器的容量约束是：
P_dg² + Q_dg² ≤ S_dg_max²
这是一个旋转锥约束，YALMIP里可以写成：
Constraints = [Constraints, [Pg(i,t); Qg(i,t)] <= sqrt(Sdg_max(i)) ]; % 示例 % 更规范写法： Constraints = [Constraints, cone([Pg(i,t); Qg(i,t)], Sdg_max(i))];当运行在功率因数控制模式下，还要加一个线性近似约束：Q_dg 在 -tan(φ_max)·P_dg 和 tan(φ_max)·P_dg 之间。如果DG支持电压支撑，那么Q_dg可以作为优化变量，由目标函数决定它的最优值。
我实际在33节点系统里接了三台DG和一组储能，位置和参数如下：
| 编号 | 接入节点 | 类型 | 额定容量（kVA） | 有功上限（kW） |
|---|---|---|---|---|
| DG1 | 18 | 光伏 | 800 | 600 |
| DG2 | 22 | 风机 | 1000 | 800 |
| DG3 | 25 | 微型燃气轮机 | 600 | 500 |
| 储能 | 33 | 电池 | 800 | 600 |
DG的有功上限设置成比额定容量小一些，是为了给无功调节留出空间，这个细节做无功优化的人很容易漏。
3.2 储能的时序耦合约束
储能和DG最大的区别在于它有时序耦合特性——今天某一时刻怎么充放电，会直接影响后续时段的可用容量。所以在多时段优化模型里，储能SOC（荷电状态）是一个连续递推变量：
SOC(t+1) = SOC(t) - P_bat(t) * Δt / E_bat
其中P_bat为正表示放电，为负表示充电。更精确的做法是区分充放电效率：
SOC(t+1) = SOC(t) - (P_bat_chg(t) * η_chg + P_bat_dis(t) / η_dis) * Δt / E_bat
同时要加充放电功率上下限和SOC上下限：
- -P_chg_max ≤ P_bat(t) ≤ P_dis_max
- 0.2 ≤ SOC(t) ≤ 0.9
SOC下限定在0.2而不是0，是从电池寿命和应急备用的角度设置的，这也是“灵活性”在储能约束上的直接体现。把SOC用完确实能降低当天的购电成本，但调度上完全没有应变空间，真发生预测偏差就无力回天了。
3.3 联络线功率的区间约束与爬坡约束
配电网通过PCC（公共连接点）连到上级电网或者主微网，PCC交换功率的大小和变化速度，是灵活性最直观的表征。在模型里我会加两类约束。
区间约束：
Pcc_min ≤ Pcc(t) ≤ Pcc_max
爬坡约束：
-ramp_cc ≤ Pcc(t+1) - Pcc(t) ≤ ramp_cc
爬坡约束相当于是说“上级电网愿意配合你，但不能瞬间配合”，这个约束在传统OPF里很少出现，恰恰是“考虑灵活性”的关键所在。不加爬坡约束时，优化器会倾向于在光伏出力突降的时段让上级电网瞬间补足功率缺口，解出来局部最优但完全不可执行。加上爬坡约束后，储能和可调DG才真正发挥作用。
3.4 目标函数设计：经济性与调节裕度的加权
优化目标可以做成很多种。如果只做纯经济调度，目标是购电成本加DG运行成本最小：
min ∑ₜ ( c_pcc(t)·Pcc(t) + ∑ᵢ c_dg(i)·Pg(i,t) )
但“考虑灵活性”之后，一般会在目标里增加调节裕度惩罚项。比如对联络线功率增加一个软约束惩罚：
min ∑ₜ ( c_pcc·Pcc + ∑ᵢ c_dg·Pg ) + λ·∑ₜ ( Pcc(t) - Pcc_ref )²
或者对电压调节裕度加惩罚：
min ... + μ·∑ₜ ( V_min(t) - 0.95 )₋²
提醒一句：λ和μ的取值需要做敏感性分析，不能拍脑袋定。我常用的做法是先跑一遍不加惩罚的经济调度，观察Pcc和最低电压的自然边界，再把惩罚系数从1e-3开始逐步放大，看解的变化趋势。系数太大会导致基本失去经济性，太小则灵活性约束形同虚设。
4. Matlab实现架构：数据、建模、求解与结果回读
4.1 工具链选型：YALMIP + MOSEK组合
Matlab的Optimization Toolbox支持线性规划和二次规划，但配电网SOCP这种锥规划问题，特别是节点规模上百个、时段数24个之后，自带工具就不够用了。我在项目里用的是YALMIP加MOSEK的组合。
YALMIP是一个Matlab的建模语言接口层，最大的优势是用sdpvar声明优化变量，然后就像写数学模型一样直接写约束和目标函数，最后一行optimize调用求解器。MOSEK是针对SOCP和SDP等凸优化做过深度优化的商业求解器，对学术研究提供免费许可证，速度和数值稳定性都比SeDuMi、SDPT3这类纯学术求解器明显更好，尤其在多时段模型里差距很大。
有个细节值得提：如果学校或公司没有MOSEK商业许可，用SeDuMi跑通流程完全没问题，但一旦时段数超过24，SeDuMi的速度可能会慢5到10倍。所以有条件的话，建议尽早申请MOSEK学术版。
4.2 数据组织：用结构体管理配电网参数
代码不是一次性写出来的，中间如果不组织好数据，后面做场景对比时会被数据格式搞得很痛苦。我习惯把配电网参数放在一个struct里统一管理：
mpc.bus = [ % 节点编号 有功负荷 无功负荷 电压下限平方 电压上限平方 1 0.000 0.000 0.9025 1.1025 2 0.100 0.060 0.9025 1.1025 ... ]; mpc.branch = [ % 首端节点 末端节点 支路电阻 支路电抗 1 2 0.0922 0.0470 2 3 0.4930 0.2511 ... ]; mpc.dg = [ % 类型 接入节点 有功上限 无功上限 1 18 0.600 0.600 2 22 0.800 0.800 3 25 0.500 0.500 ]; mpc.ess = struct('bus', 33, 'Emax', 0.6, 'Pmax', 0.3, 'soc0', 0.5, 'eta', 0.95);注意是把所有功率统一换算成标幺值。我用的基准功率是1 MVA而不是常见的100 MVA，这样决策变量的数值都在0到1的范围内，MOSEK的数值稳定性更好，求解速度也更快。用100 MVA基准时，33节点系统里支路潮流都是0.001级别的数，跟约束边界差了好几个数量级，求解器很容易出现数值告警。
4.3 YALMIP建模核心代码
下面给出一个多时段配电网SOCP-OPF的建模骨架，这就是整个Matlab实现的发动机部分：
% 决策变量：支路潮流、电压平方、电流平方、DG出力、储能功率、联络线功率 P_br = sdpvar(nBranch, T, 'full'); Q_br = sdpvar(nBranch, T, 'full'); v = sdpvar(nBus, T, 'full'); I2 = sdpvar(nBranch, T, 'full'); Pg = sdpvar(3, T, 'full'); Qg = sdpvar(3, T, 'full'); Pb = sdpvar(1, T, 'full'); SOC = sdpvar(1, T+1, 'full'); Pcc = sdpvar(1, T, 'full'); Qcc = sdpvar(1, T, 'full'); Constraints = []; for t = 1:T % 1. 支路电压降落方程 for k = 1:nBranch i = from(k); j = to(k); Constraints = [Constraints, ... v(j,t) == v(i,t) - 2*(r(k)*P_br(k,t) + x(k)*Q_br(k,t)) ... + (r(k)^2 + x(k)^2)*I2(k,t)]; % 2. 二阶锥松弛不等式：I2 >= (P^2 + Q^2) / v Constraints = [Constraints, ... cone([2*P_br(k,t); 2*Q_br(k,t); v(i,t)-I2(k,t)], v(i,t)+I2(k,t))]; end % 3. 节点功率平衡 for i = 1:nBus % 收集流向该节点的支路和流出该节点的支路 % P_Gi - P_Li = sum(P_out) - sum(P_in - r*l) end % 4. 电压上下限（平方形式） Constraints = [Constraints, 0.95^2 <= diag(v(:,t)) <= 1.05^2]; % 5. DG容量锥约束 for dg = 1:3 Constraints = [Constraints, cone([Pg(dg,t); Qg(dg,t)], S_dg_max(dg))]; end % 6. 储能约束 if t > 1 Constraints = [Constraints, ... SOC(t+1) == SOC(t) - Pb(1,t)/ESS.eta/ESS.Emax]; end Constraints = [Constraints, 0.2 <= SOC(t+1) <= 0.9]; Constraints = [Constraints, -ESS.Pmax <= Pb(1,t) <= ESS.Pmax]; % 7. 联络线爬坡约束 if t > 1 Constraints = [Constraints, ... -ramp <= Pcc(1,t) - Pcc(1,t-1) <= ramp]; end end这样写出来的模型，读起来跟论文公式几乎一一对应，后期调约束或者加新的资源类型都非常方便。YALMIP在底层会自己判断变量出现的结构，把约束自动整理成求解器需要的锥标准形式。
4.4 求解器配置与结果回读
求解配置也很关键。推荐在sdpsettings里指定求解器和容差：
options = sdpsettings('solver', 'mosek', ... 'verbose', 2, ... 'mosek.MSK_DPAR_INTPNT_CO_TOL_REL_GAP', 1e-7, ... 'savesolveroutput', 1); sol = optimize(Constraints, Objective, options); if sol.problem == 0 Pcc_opt = value(Pcc); V_opt = sqrt(value(v)); Loss_opt = sum(sum(P_br(1,:) - P_br(:,end))) ... % 网损汇总 disp('求解成功'); else disp(sol.info); endsol.problem == 0代表成功，其他返回值分别对应不可行、数值问题、求解器中断等，这里每个分支都建议打印sol.info，排查问题会方便很多。
5. 仿真算例设计与结果分析：松弛紧度、经济性与电压改善
5.1 三种对照场景的设计
算例设计的核心原则是“每个变量都能单独解释”。我在33节点系统里设置了三个场景：
- 场景A：不接入DG，原始33节点系统。
- 场景B：接入DG，但DG出力固定为预测值，不参与优化调节，联络线功率自由变化。
- 场景C：接入DG，同时配置储能、联络线爬坡约束和电压调节裕度惩罚，也就是完整版的“灵活性优化”。
其中场景B到场景C的差异，可以单独分离出“灵活性约束”对结果的影响，这是论文里最常用的对比方式。DG接入位置采用的是上面那套18、22、25、33节点的配置。
5.2 核心结果对照
实际跑出来的结果对比如下：
| 场景 | 网络损耗（kW） | 全网最低电压（p.u.） | 最大锥松弛间隙 | 求解时间（s） |
|---|---|---|---|---|
| 场景A：无DG | 202.7 | 0.913 | — | 0.3 |
| 场景B：DG固定出力 | 168.3 | 0.941 | 3.2e-6 | 0.5 |
| 场景C：DG+储能+灵活性优化 | 131.5 | 0.958 | 4.1e-7 | 1.2 |
从结果可以清楚看到，DG接入本身能降低网损并抬升电压，但灵活性优化的价值在于把网损又压低了一截，同时最低电压从0.941抬升到0.958。储能发挥了削峰填谷的作用，光伏大发时段充电，晚高峰时段放电，联络线爬坡约束让系统不再把光伏波动的冲击直接甩给上级电网。
这里要特别说明：具体数值和DG接入位置、容量配置强相关，换一套参数后数值会变化，但趋势是稳定可复现的。如果你的算例里DG接入位置离变电站太近或者容量太小，电压改善效果可能不明显，建议先把DG布置在馈线中后段，优化效果更容易显现。
5.3 紧度验证结果
上一节表格里最大的锥松弛间隙是一个关键指标。场景B的间隙是3.2e-6，场景C进一步降到4.1e-7，都远低于1e-3的工程接受门槛。这说明前面的二阶锥松弛在这个算例里是紧的，求出的解可以直接还原成原问题的可行潮流解。
我还额外做了一件事：把SOCP求出的解代回到原始非线性DistFlow方程里，用传统的前推回代法重新算一遍潮流，对比两次结果的电压误差。最大电压偏差在0.2%以内，这就从两个角度同时验证了松弛解的可用性。建议你在自己项目里也保留这一步，审稿人或者验收专家问起来的时候，把这个佐证拿出来特别有说服力。
5.4 多时段结果的时序行为
多时段优化的结果如果只看总网损和最低电压，信息量还不够。我习惯把Pcc、光伏出力、储能SOC画在同一张图里看时序配合。典型结果是：午间光伏大发时Pcc明显降低甚至出现小幅倒送，储能开始充电；傍晚光伏出力衰减、负荷攀升时储能转为放电，Pcc平滑上升。整个联络线功率曲线没有出现尖峰和突变，这就是爬坡约束和灵活性裕度共同起作用的结果。
绘制这种图在Matlab里用stairs或者plot都行，关键是数据提取用value()函数，再转成普通double数组。建议多时段结果存成mat文件，后续做误差分析、画散点图、做敏感性分析都不用重新跑模型。
6. 收敛性调试与松弛紧度验证的实战经验
6.1 数值缩放：一个容易被忽略的大坑
SOCP问题用内点法求解时，对数值尺度非常敏感。我第一次跑24时段模型时，MOSEK一直报数值警告，解的精度也很差。排查下来发现问题是基准功率选错了。一开始用100 MVA基准，支路电流平方l_ij的量级在1e-6，电压平方v_i的量级在1左右，锥约束里的向量元素差了6个数量级，锥体形状被压得极度扁长，内点法迭代很容易触及数值边界。
换成1 MVA基准后，所有功率在0.1到2之间，电流平方在0.01到0.2之间，数值尺度基本控制在两个数量级以内，求解稳定性明显改善，迭代次数也减少了。如果你在别的算例里遇到数值问题，优先检查决策变量的量级分布，而不是急着换求解器。目标是让所有变量尽量落在0.01到100之间。
6.2 求解不可行时的定位方法
多时段模型里约束多，某个约束写错或者参数给得不合理，很容易出现problem = 1的不可行状态。新手第一反应是检查所有约束条件，但那样效率太低。我通常用两个手段定位。
第一个手段是分批注释法。先把储能SOC递推、爬坡约束这类时序耦合约束整体注释掉，看问题是否可解。如果可以解，说明问题出在时序耦合约束上；如果还是不可行，再往下注释DG容量约束、电压约束，逐层缩小范围。这个过程虽然人工，但通常几分钟内能定位到具体的约束组。
第二个手段是检查对偶变量。YALMIP求解后，可以通过dual(Constraints(index))提取某个约束的对偶乘子。如果某个约束在最优解处是紧张的，对偶乘子会反映它的影子价格。当模型不可行时，重点检查最后被标记为冲突的约束对应的对偶信息，往往能直接看出是哪个电压下限或者容量极限被突破了。
6.3 松弛紧度检查不能只看总量
紧度验证最容易犯的错误是只列一个“所有支路最大间隙=3.0e-6”就收工。更好的做法是分时段、分支路输出间隙分布，看哪些支路在哪些时段间隙最大。比如某条靠近DG接入点的支路在光伏大发时段间隙明显比其他支路大一个数量级，这通常是该支路潮流太重导致弛豫更容易被“拉开”。这种情况即便全局最大间隙还在1e-5以内，也要在论文或交付文档里单独说明。
标准检验代码如下：
viol = zeros(nBranch, T); for t = 1:T for k = 1:nBranch i = from(k); viol(k,t) = abs(value(I2(k,t)) - ... (value(P_br(k,t))^2 + value(Q_br(k,t))^2) / value(v(i,t))); end end max_gap = max(viol(:)); fprintf('最大锥松弛间隙：%.3e\n', max_gap);这个检查脚本我会固化在每个算例的run_simulation.m主脚本里，跑完优化立刻执行，避免后面数据被改乱了还不知道。
6.4 模型扩展：动态运行与鲁棒化的方向
二阶锥松弛方法还有一个很大的优势是扩展性强。把静态的SOCP-OPF扩展成多时段模型只需要按T个时段重复约束，计算时间线性增长。再进一步，可以对光伏出力考虑不确定性集合，把问题转成鲁棒SOCP，YALMIP对鲁棒优化也有一级支持，不需要换整套实现框架。
我在这套33节点系统上测试过扩展到135节点系统，单断面求解时间从0.5秒增加到约6秒，多时段24小时大概在2分钟左右。对于规划研究这种离线计算场景，这个速度完全可以接受。如果后续要做实时调度，可以考虑把模型按时间窗滚动求解，或者把SOCP结果作为热启动初值给非线性潮流校正。
6.5 给入手这个方向的人一些代码组织建议
最后说一点代码组织层面的体会。这类项目代码如果只有一个大脚本，做到后面几乎无法维护。我建议至少拆成四个文件：
load_case33.m：读取配电网结构数据、DG参数、储能参数，输出mpc结构体。build_socp_opf.m：输入mpc和场景参数，返回YALMIP的约束、目标函数、优化变量句柄。run_scenario.m：主脚本，循环跑不同场景，保存结果到results结构体。plot_results.m：读取results，绘制电压曲线、Pcc曲线、DG出力柱状图、SOC曲线。
这样每次只改一个文件里的参数或者约束，不影响其他模块。我做场景对比时只需要在run_scenario.m里调整一个场景配置变量，其他全自动。
如果后面要写论文，这套结构也很方便做敏感性分析和对比实验。改DG容量、改储能位置、改惩罚系数，都是改一行参数的事。把代码组织成这种模块化结构，确实能节省大量重复跑数的时间。
这套方法我后续还在继续扩展，目前正在测试把无功补偿装置和OLTC抽头一并纳入SOCP模型，加上变压器变比这个离散变量之后问题会变成混合整数二阶锥规划，求解思路又不一样。等这版跑通了再单独写一篇跟大家分享。
