清楚明了的凸松弛最优潮流！基于混合整数二阶锥规划的主动配电网最优潮流研究程序代码！

前言

最优潮流(optimal power flow，OPF)问题，是电力系统中最常见、最基础的一类优化问题。在满足基尔霍夫定律、线路容量约束以及运行安全约束等电力网络物理约束的前提下，OPF问题旨在寻找一个最优的潮流稳态工作点，使得在该工作点上系统总发电成本、总网损等目标函数达到最优。OPF问题的相关研究在电力系统的经济调度、机组组合、状态估计、稳定性/可靠性评估、电压/无功控制、需求侧响应等关键领域得到了极为广泛的应用。由于模型中二次潮流等式约束的非凸性，OPF问题是一个难以求解的非凸优化问题。其在求解过程中容易陷入局部最优，并被证明是一个NP-hard问题。因此，研究能够同时保证求解效率与质量的OPF求解方法一直是电力系统领域乃至整个优化领域最为关注的热点之一。

OPF求解方法分类

目前，OPF问题的求解方法主要可分为以下3类：

1）寻求局部最优解。

早期的经典解法，如简化梯度法、牛顿法、序列二次规划法、内点法以及近年来发展迅速的启发式算法，都属于该类求解方法。然而，由于OPF问题本身是非凸的，这些方法无法保证解的质量，在得到局部最优解后无法衡量其与全局最优解之间的差距，并且即使得到全局最优解也无法证明其全局最优性。

2）对潮流等式约束作近似处理。

该类方法主要针对OPF模型非凸性的来源，对潮流等式约束进行线性化近似处理。其中具有代表性的方法是将交流潮流约束近似为线性的直流潮流约束，然后求解近似后的直流最优潮流(DC optimal power flow，DC-OPF)问题。其存在二个明显缺陷：一是较难应用于电压、无功相关的研究领域以及馈线R/X较高的配电网中；二是DC-OPF问题的最优解并不一定是原OPF问题的可行解，导致实际优化过程中需要不断调整DC-OPF部分约束的松紧程度并重新求解。

3）对潮流等式约束进行凸松弛。

这类方法有以下几个优势：一是在松弛精确的情况下，可保证在多项式求解时间内获得原问题的全局最优解；二是在松弛不够精确的情况下，可以为原优化问题提供最优值的一个下界；三是可以通过松弛后的凸问题无可行解来判定原非凸问题无解。然而，对于凸松弛技术来说，一个至关重要的问题是如何保证松弛的精确性。在松弛不精确的情况下，所得到的解没有实际的物理意义。

凸松弛技术概述

2006年，Jabr首次运用SOCP松弛技术将辐射状电网的潮流模型转化为一个二阶锥模型，并随后将其推广至环网，标志着现代凸松弛理论在OPF应用中的开端。Bai等将经典OPF问题转化为一个二次约束二次规划问题，并首次将其近似为SDP问题进行求解，让SDP规划进入电力领域优化者的视野之内。随后，Lavaei和Low首次提出了使用SDP技术求解OPF过程中“精确松弛”的概念，并在IEEE测试系统中完成了对所提精确松弛充分条件的检验，表明了凸松弛技术在电力系统优化领域的应用潜力。Low对截至2014年的OPF凸松弛理论首次做了系统性的梳理，给出了OPF模型中SOCP、SDP松弛的完整数学形式，比较了两类松弛各自的特点和适用范围，总结了精确松弛的充分条件，使得这两篇文献成为了引用甚广的OPF凸松弛基础教程。

在之后的研究中，出现了基于凸包络(convex hull/envelope)的凸松弛技术。在凸优化理论中，一个非凸区域的凸包络，被定义为该区域最紧的凸松弛。这种松弛技术与SDP松弛相比，需要的计算量更小，通常通过线性松弛或SOCP松弛来求解原OPF问题。针对极坐标系下OPF模型中含三角函数的约束。

尽管凸松弛技术的理论基础严格，研究热度也逐年上升，但其在实际OPF应用中仍存在较大的困难。首先是来自计算量方面的挑战。对于节点数成百上千的实际电力系统来说，直接运用内点法求解松弛后SDP问题的计算量过大，计算效率较低。

SDP与SOCP等凸松弛技术在数学领域中可被概括为经典的“升维松弛–返回映射”(lift-and-project)过程，即首先引入新变量(如SDP中的矩阵变量)完成对原问题的升维，然后在升维后的问题中松弛掉非凸约束并求解，最后从中通过返回映射来恢复原问题的解。由此过程可以看出，凸松弛在应用过程中的另一障碍是能否由松弛后的解恢复出原OPF解，这与松弛是否精确是等价的。若松弛不够精确，求解OPF仅得原问题(目标函数最小化)最优值的一个下界，松弛所得的高维矩阵解没有明确的物理含义，无法从中恢复原OPF解。因此，在使用凸松弛技术时，如何确保松弛的精确程度是亟需解决的核心问题。对于这类问题，目前有2种研究方向：1）寻求精确凸松弛的充分条件，即不对松弛过程本身作处理，而是寻求外生的、被证明可以保证松弛精确的应用场景；2）对松弛过程本身进行“加强”，构造尽可能紧的凸松弛，使其精确度内生化地提高。

凸松弛技术对比

根据已有的研究成果定性地刻画了3类凸松弛模型解集之间的关系，如下图所示。从图中可以看出，SOCP松弛模型的解集包含SDP松弛和QC松弛模型的解集，而SDP松弛和QC松弛模型的解集互不包含。换言之，一般情况下SDP松弛和QC松弛对潮流方程进行的松弛均比SOCP松弛更紧，但是SDP松弛和QC松弛两者的松紧程度需要视情况而定。

值得说明的是，不同凸松弛技术的精确程度在2种情况下会呈现完全相同的关系：1）当相角差的边界条件被逐渐放开时，三角函数项的凸包络构造效果与该边界条件有关，QC松弛的精确程度会越来越低，直至最终与SOCP松弛相同；2）当电网呈辐射状结构时，SDP松弛与SOCP松弛的精确程度相同。SOCP松弛没有SDP松弛精确的原因是其额外松弛了“相角循环条件”(cycle condition)。该条件直观的物理意义是对于电网中的每一个回路，按一定方向(如顺时针)循环累加的电压相角差之和对2取模后的结果等于0。由于辐射状网络中不含环，无需考虑“相角循环条件”，因此SOCP与SDP松弛的精确程度相同。

对于一个求最小值的OPF问题，内点法所得最优值≥全局最优值≥凸松弛所得最优值，因此可考虑使用凸松弛与内点法所得最优值的差异来衡量松弛解与原全局最优解的差异，即松弛的精确程度。其中优化间隙(optimalitygap)的定义为

凸松弛的精确性

由最优化理论可知，凸松弛的精确性与原OPF问题的目标函数及边界条件有着密切的关系。然而，在实际应用场景中，由于目标函数各不相同，且出于应用需要在模型中设定的边界条件多种多样，期望通过一个比较普适化的方法来保证凸松弛的精确性是极为困难的。针对凸松弛的精确性，目前主要有2种研究方向：一种是在理论上寻求可以保证精确凸松弛的充分条件，即找到目标函数和边界条件满足特定要求的OPF应用场景；还有一种是在SDP、SOCP等松弛技术的基础上构造更紧更精确的松弛。

保证精确松弛的充分条件

目前关于精确松弛充分条件的研究多选择在拓扑简单的辐射状电网中进行，因为从理论上推导复杂电网的精确松弛充分条件并完成相应的数学证明有较大的难度。选择的凸松弛技术对象主要为SOCP松弛，因为SOCP松弛在辐射状电网中与SDP松弛的精确度相同但计算效率远高于SDP。

研究结果表明，辐射状电网中SOCP精确松弛的充分条件对OPF模型的目标函数和边界条件有着较为苛刻甚至不符实际的要求。早期文献的研究成果按照充分条件对OPF模型边界条件的要求总结为3类：

1）对节点净注入功率(发电功率减去负荷功率)约束的要求。

对于电网中任一线路两端节点，共计8条净注入功率约束(2个节点、有功/无功、上界/下界)，且每条净注入功率约束与复平面单位圆上的点一一对应。如果每对节点的8条约束对应的点在复平面单位圆周上均具有“线性可分”(linear separability)的性质，即处于复平面内某条经过原点的直线的同一侧，那么整个OPF模型的SOCP松弛就是精确的。这一类条件的几个特例为：允许所有节点的负荷功率达到无穷大 (松弛所有节点的净注入功率下界约束)；在OPF模型中仅保留部分节点注入的有功功率下界约束。可以看出此类充分条件往往不符实际，难以应用。

2）对节点电压幅值约束的要求。

若原OPF问题没有电压幅值上界约束，或是松弛模型解的任一节点电压幅值都未达到预设的上界，且电网中不存在反向潮流，则SOCP松弛是精确的。与忽略电压幅值约束相比，检验“电压幅值是否达到上界”这一条件具有更强的可操作性，但仅为“后验式”条件。

3）对线路两端相角差约束的要求 。

通过研究不同边界条件下OPF可行域的几何性质提出，当线路两端的电压相角差足够接近时，原可行域与松弛后凸包络区域的帕累托最优集(最优解集)可能完全相同。涉及相角差约束的充分条件通常也包含部分对节点幅值或者节点净注入功率的要求。

更精确松弛的构造方法

该领域近年来的另一种代表性方法是通过矩(moment)信息来构造更紧的SDP松弛。该方法提出的多项式优化理论，使用了高维的矩信息，以增加计算量的代价提高了SDP松弛的精确程度。然而，该方法在使用中有2点不足：一是通过该方法获得全局最优解所要求的矩信息阶数非常高，会将SDP问题原本就比较大的计算量提升到一个更难以处理的程度；二是得到的SD矩阵解不一定是秩为1的矩阵，难以从中恢复出原OPF问题的最优解。

除此之外，在SDP模型的目标函数中加入了与矩阵解的秩相关的罚项来尽可能获得一个秩为1的半定矩阵解，从而提高松弛的精确度。

程序介绍

程序建立考虑DG、ESS、SVC和CB等连续、离散控制变量，以降低配电网运行成本、提高DG并网能力、消除过电压为目的的主动配电网有功无功协调的三相多时段优化模型；然后采用二阶锥松弛技术(second-order cone relaxation， SOCR)将其中的三相潮流方程作凸化松弛处理，将原问题转化为一个可被有效求解的混合整数二阶锥优化问题；最后，采用扩展的 IEEE 33 节点三相测试系统仿真计算，采用Cplex等算法包求得原问题的最优解；程序中算例丰富，注释清晰，干货满满，创新性和可扩展性很高，足以撑起一篇高水平论文！下面对程序做简要介绍！

程序适用平台：Matlab+Yalmip+Cplex

参考文献：《基于混合整数二阶锥规划的主动配电网有功–无功协调多时段优化运行》-中国电机工程学报

程序结果

部分程序

%% 1.设参
mpc = IEEE33BW;wind = mpc.wind;    pload = mpc.pload;    
pload_prim = mpc.pload_prim/1000;  %化为标幺值
pload = pload/a;%得到各个时段与单时段容量的比例系数
qload = pload/b;%假设有功负荷曲线与无功负荷变化曲线相同
pload = pload_prim*pload;   %得到33*24的负荷值，每一个时间段每个节点的负荷     
branch(:,3) = branch(:,3)*1/(12.66^2);%求阻抗标幺值              
T = 24;%时段数为24小时              
0.94*0.94*ones(1,T)];%加入变压器后，根节点前移，因此不是恒定值1.06
QCB_step = 100/1000;  %单组CB无功,100Kvar 转标幺值     
%% 2.设变量
V = sdpvar(nb,T);%电压的平方I = sdpvar(nl,T);%支路电流的平方
P = sdpvar(nl,T);%线路有功（是不是平方我就不清楚了，应该不是）
Q = sdpvar(nl,T);%线路无功Pg = sdpvar(nb,T);%发电机有功Qg = sdpvar(nb,T);%发电机无功
theta_CB = binvar(ncb,T,5); %CB档位选择，最大档为5
theta1_IN = binvar(noltc,T);%OLTC档位增大标识位theta1_DE = binvar(noltc,T);%OLTC档位减小标识位
%% 3.设约束   %% 储能装置（ESS）约束     %充放电状态约束        
C = [C, u_dch + u_ch <= 1];%表示充电，放电，不充不放三种状态
%功率约束
C = [C, 0 <= p_dch(1,:) <= u_dch(1,:)*0.3];C = [C, 0 <= p_dch(2,:) <= u_dch(2,:)*0.2];​
%容量约束
C = [C, E_ess(:,t+1) == E_ess(:,t) + 0.9*p_ch(:,t) - 1.11*p_dch(:,t)];   %效率   
%投入节点选择（两电池充放电状态）
P_dch = [zeros(14,T);p_dch(1,:);zeros(16,T);p_dch(2,:);zeros(1,T)];   %电池放在第15节点和第32节点
%% 风机（光伏)约束             ​     
P_wt = [zeros(16,24);p_wt(1,:);zeros(14,24);p_wt(2,:);zeros(1,24)];     %风机放在17和32节点
%% 有载调压变压器（OLTC）约束        
rjs = zeros(1,12);%相邻2个抽头的变比  平方之差     
C = [C, r1(1,t) == 0.94^2+ sum(rjs.*theta_OLTC(1,t,:))];  %%各个档位变比dita^2 * 开关档位状态
C = [C, theta_OLTC(:,:,i) >= theta_OLTC(:,:,i+1)];   %0下面不能有1
% theta_OLTC = value(theta_OLTC);
C = [C, V(33,:) == r1];   %%最大值是1.06^2,放在  主网到33节点   
k = sum(theta_OLTC,3);     %有载调压变压器投切状态个数求和 （1*24）   
C = [C, k(:,t+1) - k(:,t) <= theta1_IN(:,t)*12 - theta1_DE(:,t) ];  %升压不能超过12档，减压不能小于1档
C = [C, sum(theta1_IN + theta1_DE,2) <= 5 ];  %限制有载调压变压器日调节次数为5次
%% 连续无功补偿装置（SVC）约束      
Q_SVC = [zeros(4,T);q_SVC(1,:);zeros(9,T);q_SVC(2,:);zeros(15,T);q_SVC(3,:);zeros(2,T)];%SVC投入节点选择5、15、31
%% 离散无功补偿装置（CB）约束     
Q_cb = sum(theta_CB,3).*QCB_step;     

部分内容源自网络，侵权联系删除！

欢迎感兴趣的小伙伴关注并私信获取完整版代码，小编会不定期更新高质量的学习资料、文章和程序代码，为您的科研加油助力！