基于内点法的最优潮流计算

基于内点法的最优潮流计算Prepared on 22 November 2020摘要内点法是一种能在可行域内部寻优的方法，即从初始内点出发，沿着中心路径方向在可行域内部直接走向最优解的方法。

其中路径跟踪法是目前最具有发展潜力的一类内点算法，该方法鲁棒性强，对初值的选择不敏感，在目前电力系统优化问题中得到了广泛的应用。

本文采用路径跟踪法进行最优求解，首先介绍了路径跟踪法的基本模型，并且结合具体算例，用编写的Matlab程序进行仿真分析，验证了该方法在最优潮流计算中的优越性能。

关键词：最优潮流、内点法、路径跟踪法、仿真目次0、引言电力系统最优潮流，简称OPF（Optimal Power Flow）。

OPF问题是一个复杂的非线性规划问题，要求满足待定的电力系统运行和安全约束条件下，通过调整系统中可利用控制手段实现预定目标最优的系统稳定运行状态。

针对不同的应用，OPF模型课以选择不同的控制变量、状态变量集合，不同的目标函数，以及不同的约束条件，其数学模型可描述为确定一组最优控制变量u，以使目标函数取极小值，并且满足如下等式和不等式。

{min u f(x,u)S.t.ℎ(x,u)=0g(x,u)≤0（0-1）其中min u f(x,u)为优化的目标函数，可以表示系统运行成本最小、或者系统运行网损最小。

S.t.ℎ(x,u)=0为等式约束，表示满足系统稳定运行的功率平衡。

g(x,u)≤0为不等式约束，表示电源有功出力的上下界约束、节点电压上下线约束、线路传输功率上下线约束等等。

电力系统最优潮流算法大致可以分为两类：经典算法和智能算法。

其中经典算法主要是指以简化梯度法、牛顿法、内点法和解耦法为代表的基于线性规划和非线性规划以及解耦原则的算法，是研究最多的最优潮流算法, 这类算法的特点是以一阶或二阶梯度作为寻找最优解的主要信息。

智能算法主要是指遗传算法和模拟退火发等，这类算法的特点是不以梯度作为寻优信息，属于非导数的优化方法。

因此经典算法的优点是能按目标函数的导数信息确定搜索方向，计算速度快，算法比较成熟，结果可信度高。

缺点是对目标函数及约束条件有一定的限制，可能出现局部极小时难以收敛。

而智能算法的优点是计算与导数无关，灵活性高，随机性强，缺点是算法不稳定，结果不可信，并且控制参数需凭经验给出。

通过对这些常见算法的简单比较，内点法具有其优越的性能，特别是路径跟踪法，其算法收敛迅速，鲁棒性强，对初值的选择不敏感，其迭代次数与系统规模或控制变量的数目关系不大，因此本文采用该方法进行最优计算。

1、路径跟踪法的基本数学模型内点法最初的基本思路是希望通过寻优迭代过程始终在可行域内进行，因此，初始点应在可行域内，并在可行域的边界设置‘障碍’使迭代点接近边界时其目标函数迅速增大，从而保证迭代点均在可行域的内点。

但是对于大规模实际问题而言，寻找初始点往往十分困难。

为此许多学者长期以来致力于内点算法初始‘内点’条件的改进。

以下介绍的路径跟踪法只要求在寻优过程中松弛变量和拉格朗日乘子满足简单的大于0或者小于0的条件，即可代替原来必须在可行域内求解的要求，使得计算过程大为简化。

一般可以将最优潮流模型简化为如下的非线性优化模型。

Obj. min u f(x,u)（1-1）. S.t.ℎ(x,u)=0（1-2）g−≤g(x,u)≤g−（1-3）其中min u f(x,u)为优化的目标函数， S.t.ℎ(x,u)=0为等式约束， g(x,u)为不等式约束，路径跟踪内点法的基本思路是：首先将式（1-3）的不等约束变成等式约束：g (x,u )+u =g − （1-4） g (x,u )−l =g − （1-5）其中松弛变量 l =[l 1,…,l r ]T， u =[u 1,…,u r ]T，应满足u>0，l>0这样原问题就转化为问题A ：Obj. min u f (x,u ). ()0()()h x g x u g g x l g=+=-=然后，把目标函数改造成障碍函数，该函数在可行域内应接近于原函数f(x)，而在边界时变得很大。

一次可得带优化问题B ：obj. 11min.()log()log()rrr r j j f x u l u u ==--∑∑. ()0()()h x g x u g g x l g=+=-=其中扰动因子或者障碍因子u>0。

当l 或u 接近边界时，以上函数将趋于无穷大，因此满足以上障碍目标函数的极小解不可能在边界上找到。

这样就通过目标函数的变化把含不等式限制的优化问题A 变成只含等式限制优化的问题B 了，因此可以直接用拉格朗日乘子法来求解。

优化问题B 的拉格朗日目标函数为：11()()[()][()]log()log()r rTTTr r j j L f x y h x z g x l g w g x u g u l u u ===-----+---∑∑ 式中：y ,z 和w 均为拉格朗日乘子。

因此最后简化的求解问题就是求取上述表达式的极小解。

2、 路径跟踪法的最优潮流求解思路路径跟踪法的最优潮流求解过程就是对拉格朗日目标函数求极小值问题： 式中：y ,z 和w 均为拉格朗日乘子。

该问题极小值存在的必要条件是拉格朗日函数对所有变量及乘子的偏导数为0。

即：11()()()(+w)0()0()0()+u 00U 0μμμμμμ--∂==∇-∇-∇=∂∂===∂∂==--=∂∂==-=∂∂==-⇒=-=∂∂==--⇒=+=∂x x x x y z w l lu uLL f x h x y g x z x L L h x y L L g x l g z L L g x g w L L z L e L LZe e l L L w e L UWe e u （2-1）通过上述表达式可以解出：2T T l z u w rμ-= （2-2）定义：T T Gap l z u w =-，称为互补间隙。

可得：2Gaprμ=（2-3）如果x *是优化问题A 的最优解，当u 固定时，x(u)是优化问题B 的解，那么当Gap →0，u →0时，产生的序列{x(u)}收敛至x *。

建议采用：2Gaprμσ=。

式中(0,1)σ∈称为中心参数，一般取，在大多数场合可获得较好的收敛效果。

通过偏导数为0的表达式可以可得内点法的修正方程为：111'0000000()000000000()00000'()0000()0l T z x u Tw x x x T y x L L z I L ZL l I g x U L w I L u I g x L x H h x L y h x μμ----∆⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥∆-∇⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-∆⎢⎥=⎢⎥⎢⎥⎢⎥-∆∇⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥∆∇⎢⎥⎢⎥⎢⎥-∆∇⎢⎥⎣⎦⎣⎦⎣⎦（2-4） 求解方程可得到第k 次迭代的修正量，于是最优解的一个新的近似解为：(1)()(1)()(1)()(1)()(1)()(1)()k k p k k p k k p k k d k k d k k d x x a x l l a l u u a u yya yz z a z w w a w++++++=+∆=+∆=+∆=+∆=+∆=+∆ （2-5）式中，p a 和d a 为步长：0.9995min{min(,0;,0),1}i i p i i i il ua l u l u --=∆<∆<∆∆ 0.9995min{min(,0;,0),1}i i d i i i iz wa z w z w --=∆<∆>∆∆ （2-6） 其潮流计算的流程图如下图1所示，其中初始化部分包括： (1)、设置松弛变量l 和u ，保证[l ，u]T >0(2)、设置拉格朗日乘子w 、y 、z ，满足[w<0,z>0，Y ！=0]T (3)、设置优化问题的初值。

(4)、取中心参数(0,1)σ∈，给定计算精度，迭代次数初值K=0。

3、 具体算例及程序实现流程这部分主要有算例描述以及程序的实现流程两部分，其中算例描述主要是对系统参数以及优化问题进行说明。

而程序的实现流程主要描述的是最优潮流计算中所涉及的矩阵方程的描述。

、算例描述该算例为王锡凡编写的《现代电力系统分析》中的3-1的例题，是以系统燃料最省为最优潮流的目标函数。

选择该题目作为算例分析的原因是，该题目有比较详细的解题思路以及列写出了比较详细的迭代结果，方便对编写程序的运行结果进行比对。

求如下图所示简化系统的系统燃料最省的最优潮流计算：除了由上图所提供的系统母线负荷功率数据、线路参数和变压器之路参数数据、变压器便比数据之外，以下顺序给出了线路传输功率边界（表3-1），发电机有功无功出力上下界和燃料耗费曲线 参数（表3-2）。

若不作特殊说明，2435所有数据都是以标幺值形式给出，功率基准值为100MVA，母线电压上下界分别为和。

表3-1线路传输功率边界表3-2 发电机数据、程序具体实现过程针对上述系统，首先我们先列写出该算例的数学模型和有关计算公式。

在该算例中，共有5个节点，相应的状态量为：系统中有2台发电机，没有其他无功源，因此控制变量为：应该指出，此处发电机和无功源的编号与及诶单编号无关，是独立编号的。

这是因为系统中一个节点可能接有多台发电机的缘故。

因此系统中总变量共有14个：最优潮流的数学模型为：目标函数：约束条件：每个节点有两个潮流方程，共有10个等式约束条件，对非发电机而言：5151(cos sin )0(sin cos )0θθθθ===--+==-+-=∑∑i Di i j ij ij ij ij j i Di i j ij ij ij ij j P P V V G B Q Q V V G B （i=1,2,3）对发电机节点：5151(cos sin )0(sin cos )0θθθθ∈=∈==--+==-+-=∑∑∑∑i Gk Di i j ij ij ij ij k ij i Gk Di i j ij ij ij ij k ij P P P V V G B Q Q Q V V G B （i=4,5）式中：∈k i 表示第k 台发电机接在节点i 上。

不等式约束共有14个条件，分别是：根据以上模型可以形成修正方程。

该方程包括形成等式左边的系数矩阵和等式右边的常数项两部分。

1、形成系数矩阵1）、等式约束的雅克比矩阵 等式右端包括3个子矩阵： 其中： 其中：式中：i 为发电机的序号；j 为节点号；()∈i j 表示第i 台发电机是在节点j 上的。

111111111111111155551111555511111010.................θθθθθθθθ⨯∂∆∂∆∂∆∂∆⎡⎤⎢⎥∂∂∂∂⎢⎥⎢⎥∂∆∂∆∂∆∂∆⎢⎥∂∂∂∂⎢⎥∂⎢⎥=∂⎢⎥∂∆∂∆∂∆∂∆⎢⎥⎢⎥∂∂∂∂⎢⎥∂∆∂∆∂∆∂∆⎢⎥⎢⎥∂∂∂∂⎣⎦P Q P Q P P P P V V V V h x P Q P Q P Q P Q V V V V （潮流计算中的雅克比矩阵）2）、不等式约束的雅克比矩阵式中：1g 、2g 、3g 和4g 依次表示电源有功出力的上下界约束，无功电源出力的上下界约束，节点电压赋值的上下界约束和线路潮流约束。