最近有不少做配电网方向的同学在问最优潮流(OPF)怎么用二阶锥松弛来求解,特别是怎么在Matlab里把代码跑通。这确实是个很实用的课题,不管是做毕业设计、发论文,还是工程前期分析,二阶锥松弛(SOCP)都是处理配电网最优潮流问题的主流方法之一。这篇文章我就把这个课题从头到尾讲透,包括数学原理、建模过程、代码实现和常见避坑经验,希望能帮到你。
1. 整体设计思路:为什么不直接用传统最优潮流算法
1.1 配电网最优潮流到底难在哪
最优潮流的本质是在满足潮流方程、电压限值、线路容量等一系列约束的前提下,通过调节发电机出力、无功补偿等控制手段,让某个目标函数达到最优——比如网损最小、运行成本最低。传统的输电网最优潮流用线性规划、牛顿法这类方法已经比较成熟,但到了配电网层面,问题就变得麻烦起来。
麻烦的根源很简单:配电网的潮流方程是非线性的,而且是非凸的。节点电压方程里面有电压幅值的平方项、相角差三角函数项,这些东西混在一起,导致整个可行域不是凸集。非凸意味着什么?意味着你用传统的梯度下降或者牛顿法去求解,得到的很可能只是一个局部最优解,而不是全局最优解。对辐射状配电网来说,线路电阻较大、网损占比高、无功流动明显,局部最优和全局最优之间的差距可能相当可观。
我打个比方,非凸优化就像你在连绵起伏的山脉里找最低点,如果只靠“沿着下坡走”的规则,很容易掉进一个山坳里就以为自己到了最低处,实际上旁边还有更深的谷底。而凸优化就像在一个大碗里找最低点,不管从哪个方向走,到达的一定是碗底——全局最优解。
这就是二阶锥松弛价值所在:它把原来非凸的潮流方程进行巧妙的变换和松弛,把整个问题转化为一个凸优化问题。凸优化有个特别好的性质:局部最优解就是全局最优解,而且有非常成熟的求解算法和商业求解器,计算效率有保障。
1.2 二阶锥松弛的基本思路:把非凸问题“掰弯”成凸问题
二阶锥松弛的核心思想并不复杂,一句话概括就是:通过变量替换,把非凸的约束改写成一个旋转二阶锥约束。
具体来说,配电网支路潮流模型(Branch Flow Model)中,每条支路从节点i流向节点j的功率可以表示为:
- P_ij 表示支路首端有功功率
- Q_ij 表示支路首端无功功率
- V_i 表示节点i的电压幅值
- I_ij 表示支路i-j的电流幅值
在支路潮流方程中,存在一个核心约束:
V_i² = V_j² + 2(R_ij * P_ij + X_ij * Q_ij) + (R_ij² + X_ij²) * I_ij²
这个式子里同时有电压平方、电流平方和功率乘积项,直接放在一起是非凸的。为了把它变凸,我们做两个变量替换:
- 令 u_i = V_i²,表示节点电压幅值的平方
- 令 l_ij = I_ij²,表示支路电流幅值的平方
替换之后,再经过Kron简化推导,会得到一个关键的不等式关系:u_i * u_j ≥ P_ij² + Q_ij²。这个不等式经过整理,恰好就是一个标准的旋转二阶锥约束:
|| [2P_ij, 2Q_ij, u_i - u_j] || ≤ u_i + u_j
这个约束用二阶锥形式表达之后,整个问题就变成了一个二阶锥规划问题(SOCP)。二阶锥规划属于凸优化,所以可以直接用成熟的商业求解器高效求解。
这里必须强调一个关键点:为什么能做这个松弛,而且松弛完了还不失真?这涉及到二阶锥松弛的“精确性”问题。对于辐射状配电网,在目标函数是关于支路电流的增函数(比如网损最小)、且节点电压没有上限约束或者上限约束不起作用时,松弛后的解会自动满足等式约束,也就是说松弛是精确的,解出来的结果等价于原非凸问题的全局最优解。这个性质在文献中已经被严格证明了,实际工程中大部分辐射状配电网的算例都满足这个条件。
2. 配电网最优潮流的数学模型构建
2.1 目标函数的选择与差异
配电网最优潮流的目标函数有很多种选择方式,不同目标对应不同的应用场景。
最常见的两种是网损最小和运行成本最小。网损最小的表达式是:
min Σ (R_ij * l_ij)
其中R_ij是支路电阻,l_ij是支路电流平方。这个目标函数的好处是它直接是二阶锥变量l_ij的线性函数,而且它天然满足松弛精确性条件——网损越小,电流越小,松弛间隙会自动趋于零。再加上网损直接反映了配电网的经济性,所以无论是学术研究还是工程分析,网损最小都是最经典的选择。
运行成本最小则通常用在含分布式电源(DG)的场景,目标函数是:
min Σ (a_i * P_gi² + b_i * P_gi + c_i)
其中P_gi表示节点i处分布式电源的有功出力,a_i、b_i、c_i是成本系数。这里P_gi是优化变量,如果成本函数是二次的,在YALMIP中处理为二次目标,整体仍然可以保持SOCP形式。
我个人的经验是:如果你是第一次做这个课题,先把网损最小的模型跑通,再去扩展含DG的多目标场景。把最基础的版本搞清楚,后面的扩展其实就是加变量、加约束的问题,模型框架不会变。
2.2 约束条件的完整梳理
配电网最优潮流的约束条件由以下几部分组成。
第一是潮流平衡约束。对每个节点i,注入功率之和等于负荷功率加上流出到下游支路的功率之和。有功和无功分别写:
Σ P_ki - Σ P_ij = P_load_i - P_gen_i
Σ Q_ki - Σ Q_ij = Q_load_i - Q_gen_i
其中k是节点i的父节点,j是节点i的子节点。这组约束保证了功率在整个网络中的流动平衡,是潮流计算的基本要求。
第二是支路电压降落方程。这是把原来节点电压方程通过变量替换改写而成的等式约束:
u_j = u_i - 2(R_ij * P_ij + X_ij * Q_ij) + (R_ij² + X_ij²) * l_ij
注意这里的正负号。在支路潮流模型中,电流从i流向j,经过线路阻抗后,电压会有压降,所以节点j的电压平方等于节点i的电压平方减去压降项。这个细节如果搞反了,整个潮流结果都会出错。
第三是二阶锥松弛约束本身:
|| [2P_ij, 2Q_ij, u_i - u_j] || ≤ u_i + u_j
这是整个模型最核心的约束,它把原问题中隐含的电流与功率关系进行了松弛处理。在YALMIP中,这个约束可以直接用旋转锥命令cone来表示。
第四是节点电压安全约束:
U_min² ≤ u_i ≤ U_max²
电压约束在配电网中特别重要,因为配电网的供电半径长、负荷点多,电压偏差问题比输电网严重得多。标准规定10kV配电网电压偏差一般不超过额定电压的±7%,在建模时可以取U_min = 0.93,U_max = 1.07(标幺值)。注意我们用的是电压幅值平方,所以边界要做平方运算。
第五是支路容量约束:
l_ij ≤ l_ij_max
线路的最大载流量限制了支路电流的平方值,这个约束既考虑了线路热稳定极限,也防止模型出现不切实际的过载结果。
2.3 为什么用辐射状网络结构
二阶锥松弛在配电网最优潮流中的一个重要前提是网络呈辐射状结构。辐射状结构意味着网络中没有环路,每个节点有且只有一个父节点,拓扑关系是一棵树。
为什么辐射状结构这么重要?因为支路潮流模型的推导本身就基于功率的单向流动关系。在辐射状网络中,从根节点(通常是变电站母线)出发,功率沿着唯一的路径流向各个负荷节点,父节点和子节点关系清晰,每条支路的潮流量可以明确写出。如果网络存在环网,功率流动方向不唯一,支路潮流方程需要额外增加环路约束才能成立,问题复杂性会大幅上升。
好消息是,实际配电网绝大多数是辐射状或近似辐射状运行的——即使有些地方有联络开关形成备用环路,正常运行时开关也是断开的。所以二阶锥松弛的方法在实际配电网中应用非常自然。如果你手头遇到的是弱环网,可以先用“线路开断”处理把环解开,再应用SOCP;或者使用扩展版的SOCP方法处理环网约束,但代码复杂度会高不少。
3. Matlab代码实现全流程解析
3.1 算例选取与数据处理
做完数学建模,接下来就是最关键的代码实现环节。我以电气工程领域最经典的IEEE 33节点配电网系统作为算例。这个系统是1991年提出的经典测试系统,包含33个节点、32条支路、根节点电压为12.66kV,总负荷约3.72MW + 2.3Mvar。系统数据网上很容易找到,MATPOWER或者各类论文附录里都有,这里我直接用一段代码把数据和拓扑结构读进来。
IEEE 33节点系统非常适合入门,因为它规模适中——节点数足够展示方法的有效性,又不会因为规模太大导致调试困难。等你把代码跑通了,切换到IEEE 69节点、119节点系统只需要替换数据文件就行,程序主体根本不用改。
首先在Matlab中定义基本参数和网络拓扑。节点的父节点编号、线路阻抗数据和节点负荷数据是三个核心数组:
%% IEEE 33节点系统数据定义 % bus, 父节点, 子节点, 电阻(ohm), 电抗(ohm), 有功负荷(kW), 无功负荷(kvar) bus_data = [ 1, 0, 1, 0.0922, 0.0470, 100, 60; 2, 1, 2, 0.4930, 0.2511, 90, 40; 3, 2, 3, 0.3660, 0.1864, 120, 80; ]; % 这里需要补充完整33条支路的数据,网络上有很多现成版本 % 节点电压基准值 V_base = 12.66; % kV % 功率基准值 S_base = 10; % MVA数据的组织方式决定后续代码的简洁程度。我的习惯是用一个矩阵存全部数据,每行代表一条支路,列依次是支路编号、父节点、子节点、电阻、电抗、节点有功负荷、节点无功负荷。这样后续在建约束时,只需要循环遍历支路编号即可。
在数据处理阶段有一个非常容易踩的坑:计算中的标幺化。IEEE 33节点系统的原始数据都是有名值(欧姆、千瓦、千乏),而在优化模型中,潮流约束中的电压平方项和功率项混合在一起,数量级差异巨大。如果直接拿有名值算,电压平方(数值上百)和线路阻抗(数值零点几)乘在一起,会导致数值条件数极差,求解器容易出现数值问题。
我的习惯是统一使用标幺值。取电压基准值V_base = 12.66kV,功率基准值S_base = 10MVA,那么阻抗基准值Z_base = V_base² / S_base = 16.02Ω,负荷标幺值等于有名值除以S_base。把所有数据都标幺化之后再建模,数值范围基本都在0到2之间,求解器跑起来很舒服。很多同学代码跑不出来,最后发现问题就出在单位的混乱上。
3.2 YALMIP建模:变量定义与约束铺设
在Matlab中求解二阶锥规划问题,我的首选工具箱是YALMIP。YALMIP是一个免费的建模工具箱,提供了统一的建模语法,能把你写的优化问题自动转换为求解器需要的标准形式。有了它,你再也不用关心怎么把锥约束转成Cplex或者Gurobi的内置格式,省去了大量琐碎的格式转换工作。
第一步是定义优化变量。这里我们定义三类决策变量:
%% 定义决策变量 n_bus = 33; % 节点数 n_branch = 32; % 支路数 P = sdpvar(n_branch, 1); % 支路有功 Q = sdpvar(n_branch, 1); % 支路无功 u = sdpvar(n_bus, 1); % 节点电压幅值平方 L = sdpvar(n_branch, 1); % 支路电流幅值平方变量的维度设计要和数据结构对齐:P、Q、L是支路相关的变量,维度是支路数;u是节点相关的变量,维度是节点数。这里要在心里建立一个映射:支路和它的首端节点、末端节点如何对应。
第二步是约束集合的定义。用Constraints = []新建一个空约束集,然后逐条添加。先加潮流方程约束和电压降落约束:
%% 约束条件 Constraints = []; % 潮流平衡约束 for k = 1:n_bus % 找到节点k的子节点集合 child_idx = find(parent_node == k); % 找到节点k的父节点 parent_idx = parent_node(k); % 构建功率平衡方程 if parent_idx ~= 0 % 非根节点,入流功率 = 流出功率 + 负荷 - 发电 Constraints = [Constraints, P(parent_idx) == ... sum(P(child_idx)) + P_load(k) - P_gen(k)]; Constraints = [Constraints, Q(parent_idx) == ... sum(Q(child_idx)) + Q_load(k) - Q_gen(k)]; else % 根节点(变电站),出流功率 = 总负荷 - 总发电 Constraints = [Constraints, sum(P(child_idx)) == ... sum(P_load) - sum(P_gen)]; Constraints = [Constraints, sum(Q(child_idx)) == ... sum(Q_load) - sum(Q_gen)]; end end这段代码的逻辑需要仔细揣摩。对非根节点,流入该节点的功率(来自父支路)应该等于流出该节点的功率(流向所有子支路)加上该节点自身的净负荷。对根节点,它没有父支路,所以潮流量直接等于总负荷减去总发电。
第三步是把支路电压方程和二阶锥松弛约束加进去。这里要用到YALMIP的cone命令。YALMIP的cone(x, y)表示约束||x||_2 ≤ y,其中x是向量,y是标量。对于旋转二阶锥约束,我们需要构造一个三维向量和一个标量:
% 支路电压降落方程和二阶锥松弛 for j = 1:n_branch i = branch_from(j); % 支路首端节点 k = branch_to(j); % 支路末端节点 % 电压降落方程 Constraints = [Constraints, u(k) == u(i) - 2*(R(j)*P(j) + X(j)*Q(j)) + ... (R(j)^2 + X(j)^2)*L(j)]; % 二阶锥松弛约束:||2P, 2Q, u_i - u_j|| <= u_i + u_j Constraints = [Constraints, cone([2*P(j); 2*Q(j); u(i) - u(k)], u(i) + u(k))]; end这里的cone([2*P(j); 2*Q(j); u(i)-u(k)], u(i)+u(k))就对应前面讲的旋转二阶锥约束。YALMIP在传递给求解器时会把旋转锥自动处理成标准二阶锥,不需要你手动做矩阵变换。但理解这个约束的来龙去脉非常重要,因为后续要检验松弛是否精确、如何分析松弛间隙,都离不开这个式子。
第四步是添加电压和容量上下限约束:
% 电压上下限约束(标幺值),配电网允许±7%偏差 U_min = 0.93^2; % 电压下限的平方 U_max = 1.07^2; % 电压上限的平方 Constraints = [Constraints, U_min <= u <= U_max]; % 支路电流上限约束(标幺值) I_max = 0.1; % 这里根据具体线路参数设定 Constraints = [Constraints, 0 <= L <= I_max^2];到这里,模型的全部约束就铺设完成了。我把约束分为四组:功率平衡、电压降落、二阶锥松弛、运行安全限值,分别对应模型中的物理潮流关系、电气关系和安全要求,逻辑非常清晰。
3.3 目标函数设置与求解器调用
目标函数选择网损最小,用支路电流平方和电阻乘积求和:
%% 目标函数:网损最小 objective = sum(R .* L);这里R是标幺化后的支路电阻向量,L是支路电流平方变量,R.*L是逐项乘积,求和就是全网的有功网损。
然后是关键的求解环节。我建议用optimize函数,并指定求解器为Cplex或者Gurobi:
%% 求解 ops = sdpsettings('solver', 'cplex', 'verbose', 2); result = optimize(Constraints, objective, ops);如果你的电脑没有安装Cplex或Gurobi这些商业求解器,也可以用开源的Sedumi或SDPT3。但我的经验是,对于这种规模的SOCP问题,Sedumi和SDPT3虽然能解,但速度比Cplex、Gurobi慢不少。33节点系统还好,如果网络规模上百节点,求解效率的差距就很明显了。
对学术研究来说,Gurobi和Cplex通常有学术免费许可,学生凭学校邮箱就可以申请,不妨下载一个用上。
verbose参数控制求解器输出的详细程度,设为2可以看到迭代过程,方便排查问题;调试完成后可以改为0,避免运行结果被大量日志刷屏。
求解完成后,检查状态标志:
%% 检查求解状态 if result.problem == 0 disp('求解成功!'); elseif result.problem == 1 disp('求解器返回不可行,请检查约束设置'); else disp('求解出错,查看result.info获取详细信息'); endresult.problem == 0表示求解成功。如果返回不可行,问题通常出在约束冲突上——例如电压下限设得太高、负荷和发电不平衡等,需要逐条检查。
3.4 结果提取与松弛间隙验证
求解完成后,用value函数提取各个优化变量的取值,然后做后处理分析:
%% 提取结果 P_opt = value(P); Q_opt = value(Q); u_opt = value(u); L_opt = value(L); % 计算网损 loss_total = sum(R .* L_opt) * S_base; % 换算到有名值 fprintf('全网有功网损:%.4f MW\n', loss_total); % 节点电压分布 V_opt = sqrt(u_opt); disp('节点电压标幺值:'); disp(V_opt);结果提取出来后,最重要的一步验证是检查二阶锥松弛是否精确。前面提到过,松弛把原问题的等式约束变成了不等式约束,只有当松弛间隙足够小、最优解自动满足等式时,得到的才是原问题的真实最优解。验证方法是对每条支路计算松弛间隙,然后取最大值:
%% 计算每条支路的松弛间隙 gap = zeros(n_branch, 1); for j = 1:n_branch i = branch_from(j); k = branch_to(j); % 松弛间隙 = (u_i + u_j)^2 - (4*P_j^2 + 4*Q_j^2 + (u_i-u_j)^2) % 简化计算:||2P, 2Q, u_i-u_j|| 应该等于 u_i+u_j lhs = norm([2*P_opt(j); 2*Q_opt(j); u_opt(i) - u_opt(k)]); rhs = u_opt(i) + u_opt(k); gap(j) = rhs - lhs; end max_gap = max(gap); fprintf('最大松弛间隙:%.6e\n', max_gap); if max_gap < 1e-5 disp('松弛精确,结果可信'); else disp('松弛间隙较大,需要检查是否有特殊约束或网络结构问题'); end松弛间隙的阈值通常取1e-5或更小。如果这个数值过大,说明最优解没有落在原问题的可行域内,你的结果可能只是松弛问题的解,而原问题在这个工况下可能有更差的目标值。此时需要检查目标函数是否严格递增、网络是否存在环网、电压上限约束是否激活等问题。
4. 拓展方向与进阶玩法
4.1 含分布式电源的扩展模型
基础代码跑通之后,最简单的进阶方向就是加入分布式电源。配电网最优潮流在“双碳”背景下,分布式光伏、风电接入的比例越来越高,配电网从无源网络变成了有源网络,潮流方向、电压分布都发生了根本性变化。
在代码中加入分布式电源只需要做两件事:第一,在原来P_gen = zeros(n_bus, 1)的位置改成一个变量,表示各节点DG的有功出力;第二,在目标函数中加入运行成本项。修改后的核心代码如下:
%% 含分布式电源的最优潮流 P_gen = sdpvar(n_bus, 1); % DG有功出力 Q_gen = sdpvar(n_bus, 1); % DG无功出力 % 约束DG出力上下限 P_gen_min = 0; P_gen_max = 2.0; % 标幺值 Q_gen_min = -1; Q_gen_max = 1; Constraints = [Constraints, P_gen_min <= P_gen <= P_gen_max]; Constraints = [Constraints, Q_gen_min <= Q_gen <= Q_gen_max]; % 目标函数:网损 + DG成本 alpha = 30; % 成本系数,权重可以调整 objective = sum(R .* L) + alpha * sum(P_gen); % 注意:潮流平衡约束中要把P_gen当作变量参与计算这里alpha是成本权重系数,它控制网损和DG成本之间的权衡。alpha取大值,模型会尽可能减少DG出力甚至不出力;alpha取小值,则更倾向于压低网损、增加DG出力。实际应用中,这个权重系数要结合上网电价、网损电价等经济参数来标定。
加入DG后,一个重要现象是电压分布会被抬高——特别是DG集中接入在馈线末端时,末端电压可能从原来的偏低变成偏高,甚至越上限。二阶锥松弛在这个场景下的一个优势是,它能同时处理电压约束和DG出力约束,自动寻找最优的DG有功无功组合。我在实际算例中测过,当DG出力足够大时,系统的网损甚至会出现负值——因为DG发出的功率就地消纳,大大减少了长距离输送的损耗。
4.2 和牛拉法、遗传算法的定位差异
做一个技术对比是很有意义的,能帮你更清楚二阶锥松弛在整个工具谱系中的位置。我在下表里列出了几种常用方法的定位:
| 方法 | 问题处理能力 | 计算速度 | 全局最优保证 | 适用场景 |
|---|---|---|---|---|
| 牛拉法(Newton-Raphson) | 非线性潮流计算 | 快 | 无 | 给定工况下的潮流分布计算 |
| 遗传算法(GA) | 各种复杂约束、混合整数 | 慢(几十秒到几分钟) | 无(启发式) | 网络重构、DG选址定容 |
| 二阶锥松弛(SOCP) | 连续变量凸优化 | 快(秒级) | 有(凸优化理论保证) | 配电网最优潮流、无功优化 |
| 混合整数二阶锥(MISOCP) | SOCP + 离散变量 | 中等 | 有(分支定界框架) | 含分段开关、电容器组的优化 |
牛拉法解决的是“给定运行状态,算出潮流分布”的问题;遗传算法这类启发式方法虽然能处理各种复杂约束,但本质上没有最优性保证,运行一次的结果可能是这次好、下次差;而二阶锥松弛站在两者中间,既能处理优化目标,又有严格的数学基础保证解的质量,而且求解速度远快于启发式方法。
对于配电网最优潮流这类问题,我的建议很直接:如果模型是连续变量的,优先试SOCP;如果需要同时处理离散决策,比如有载调压变压器、电容器投切、开关状态,那就要升级到混合整数二阶锥规划(MISOCP),在YALMIP中只需要把离散变量定义为binvar或intvar,然后指定支持MISOCP的求解器(比如Gurobi的MIP+SOCP模式),代码框架几乎可以复用。
4.3 从33节点扩展到大规模系统的性能调优
把IEEE 33节点系统的代码迁移到更大规模时,有几个细节值得提前注意。
第一个是对稀疏性的利用。当网络节点数达到几百甚至上千时,直接把所有约束用矩阵大块拼接会让YALMIP的建模速度变慢,求解器的预处理时间也会增加。此时可以考虑把循环构建的约束改为向量化表达,或者利用YALMIP的repmat、kron操作生成系数矩阵。不过对大多数学生项目来说,几百个节点的系统直接循环就行,问题不大。
第二个是求解器参数调优。Gurobi和Cplex对SOCP问题都有专门的参数控制项。比如Gurobi的PreDual、Crossover等参数会影响求解路径和收敛速度。我的经验是,处理较大规模问题时,把TimeLimit设置一个上限,比如300秒,避免某些病态算例无限期运行;同时开启MIPGap(如果是MISOCP)让求解器在可接受的次优性范围内提前终止,换取计算效率。
第三个是校验发电机出力的可行性。大规模配电网在仿真时,有时会出现DG总出力超过总负荷加网损的情况,这会导致潮流反向流动、根节点母线出现功率倒送,电压控制策略需要相应调整。在代码中加一个功率校验环节,确保求解结果满足物理常识,可以避免后续结果分析时的诸多困惑。
第四个是我特别想强调的:用规划结果做交流潮流校验。二阶锥松弛在数学上虽然大多数情况下精确,但工程上我们仍需要在结果出来后,把SOCP求得的发电机出力、无功补偿量代入回交流潮流计算程序(比如MATPOWER的runpf),逐一校验节点电压和支路潮流是否真的满足交流潮流方程。如果校验结果与SOCP结果一致,说明你的模型和代码完全可靠,这个习惯在论文审稿阶段非常加分。
5. 常见问题与排查技巧实录
5.1 YALMIP和求解器版本适配问题
很多同学问过我怎么装YALMIP和求解器,其实步骤很清晰:YALMIP从GitHub上下载压缩包,解压后把文件夹加入Matlab路径即可;Gurobi类似,安装后在Matlab里用gurobi_setup脚本配置路径。但有个细节容易出问题:版本兼容性。
YALMIP对求解器的接口通常比较稳定,但如果你用很老版本的YALMIP配最新版本的Gurobi,有时会出现“无法识别求解器”或者“符号变量类型错误”的警告。我建议是尽量从YALMIP官方GitHub仓库获取最新版本,不要用一些第三方分享的所谓“破解版”或者“绿色版”。求解器方面也一样,先确认你的Matlab版本是否支持该求解器版本——比如Gurobi 10.x对Matlab R2022a之前的版本支持可能已经丢弃,安装前先在Gurobi官网上查一下兼容性矩阵。
5.2 “Infeasible problem”——建模中最常见的坑
求解器返回不可行,是所有初学者最头疼的问题。结合我自己调试的经验,不可行的原因通常集中在以下几类:
- 负荷和发电的量纲不统一,导致功率平衡无法满足
- 电压上下限设得太紧,比如把U_min设成0.95甚至0.98,配电网在重载工况下根本达不到
- 节点数据里的父节点编号有误,导致拓扑关系错乱,出现了“孤儿节点”或“环”
- 忘记了DG出力变量需要设置上下限,默认是0,导致系统功率失衡
排查的时候有一个非常有效的操作方法:先把网络约束全部放开,比如把电压约束设为0到2,把支路容量上限设成一个很大的数,看能不能解出来。如果放开后可以解,那就说明可行域非空,问题出在你加的某些收紧约束上,然后逐步加回来,用二分法定位是哪条约束导致的不可行。这个排查思路在做优化建模时屡试不爽。
5.3 松弛不精确怎么办
如果执行完代码后,最大松弛间隙大于1e-5,说明二阶锥松弛的精确性出了问题。最常见的原因有三个:目标函数不是电流平方的严格增函数、网络结构不是辐射状、电压上限约束被激活在运行边界上。
处理思路也分三步走:先检查网络拓扑,确认是不是有环网;再检查目标函数,确认目标函数中网的变量确实是网损、且没有加一些奇怪的非线性项;最后检查结果里哪些节点的电压顶到了上限,如果在边界上,可以考虑稍微放宽电压上限,看一下松弛间隙是否会下降。
如果这些检查都没有问题,那说明这个特定工况下SOCP确实无法给出精确解,可以考虑改用SDP(半定规划)松弛,或者使用迭代求解SOCP的方法逐步逼近原问题的最优解。这些内容比较深,但对发论文的同学来说,是一个值得展开的研究点。
5.4 求解速度慢和数值异常的处理
对于中等规模配电网,SOCP的求解时间通常在秒级。如果你发现求解时间异常长,除了网络规模大以外,还有一个容易被忽略的原因:目标函数或约束中存在比较大的系数差异。比如电压平方是0.93,而网损项是0.0001,两者差了好几个数量级,求解器在内部处理时会出现数值困难。
解决办法是对模型进行归一化处理,把所有变量和约束都控制在同一个量级内。我在前面强调过标幺化,这里再次强调一次:标幺化的意义不只是物理量纲的统一,更是数值稳定性的保证。如果标幺化后仍然出现数值警告,可以在sdpsettings中设置求解器的数值稳健性选项,比如Gurobi的NumericFocus设为1或2,让求解器花更多精力在数值处理上。
6. 实操总结与个人经验
二阶锥松弛在配电网最优潮流中的应用,整体来说是一个“性价比”非常高的课题。原理层面,它只需要理解变量替换和凸松弛两个核心操作;实现层面,有YALMIP这种成熟的建模工具箱帮我们省掉了大量底层工作;效果层面,它能保证求到全局最优解,比传统启发式算法有质的优势。
我在自己跑这个课题的时候,最大的体会是:数学推导和代码实现之间是有过渡地带的。仅仅理解了二阶锥的数学原理,不代表你能顺利跑出结果——中间还隔着数据标幺化、拓扑映射、YALMIP语法、求解器参数配置这一大堆工程细节。反过来,如果只是把代码跑通,却不理解锥约束的来龙去脉,出了问题也无从下手。所以我的建议是:先手推一遍33节点的简单算例,再打开代码一句一句对照,把模型从数学形式到代码形式完整映射一遍,之后不管换系统、换目标、加DG,都很顺畅了。
最后再分享一个小技巧。在你调试SOCP代码的时候,可以用一个非常小的测试网络——3节点或5节点的链式配电网——来做正确性验证。小网络的解析解你可以手算出来,拿SOCP结果去对比,一旦代码在小网络上跑不出正确结果,就不要贸然拿到33节点上调试。小网络上验证通过后,再切换到标准算例集。我见过太多人一上来就跑IEEE 33节点系统,遇到问题根本不知道是数学模型错了还是数据读错了还是代码写错了,排错成本极高。从小处着手,逐步放大,反而是最快的路径。
希望这篇文章对正在做这个课题的你有帮助。如果还有具体问题,欢迎在评论区交流,我可以针对性地展开讲。