基于内点法的最优潮流计算.doc
《基于内点法的最优潮流计算.doc》由会员分享，可在线阅读，更多相关《基于内点法的最优潮流计算.doc（33页珍藏版）》请在冰点文库上搜索。
摘要
内点法是一种能在可行域内部寻优的方法，即从初始内点出发，沿着中心路径方向在可行域内部直接走向最优解的方法。
其中路径跟踪法是目前最具有发展潜力的一类内点算法，该方法鲁棒性强，对初值的选择不敏感，在目前电力系统优化问题中得到了广泛的应用。
本文采用路径跟踪法进行最优求解，首先介绍了路径跟踪法的基本模型，并且结合具体算例，用编写的Matlab程序进行仿真分析，验证了该方法在最优潮流计算中的优越性能。
关键词：
最优潮流、内点法、路径跟踪法、仿真
非物质文化遗产是指各族人民世代传承的，与群众生活密切相关的各种传统文化表现形式和文化空间，包括民俗活动、表演艺术、传统知识和技能以及与之相关的器具、实物、手工制品等
电力系统稳态分析课程文论第Ⅰ页共Ⅰ页
目次
0、引言 1
1、路径跟踪法的基本数学模型 2
2、路径跟踪法的最优潮流求解思路 3
3、具体算例及程序实现流程 6
3.1、算例描述 6
3.2、程序具体实现流程 7
4、运行结果及分析 11
4．1运行结果 11
4．2结果分析 16
5、结论 17
6、编程中遇到的问题 17
参考文献 19
附录 20
电力系统稳态分析课程论文第31页共30页
0、引言
电力系统最优潮流，简称OPF（OptimalPowerFlow）。
OPF问题是一个复杂的非线性规划问题，要求满足待定的电力系统运行和安全约束条件下，通过调整系统中可利用控制手段实现预定目标最优的系统稳定运行状态。
针对不同的应用，OPF模型课以选择不同的控制变量、状态变量集合，不同的目标函数，以及不同的约束条件，其数学模型可描述为确定一组最优控制变量u，以使目标函数取极小值，并且满足如下等式和不等式。
（0-1）
其中为优化的目标函数，可以表示系统运行成本最小、或者系统运行网损最小。
为等式约束，表示满足系统稳定运行的功率平衡。
为不等式约束，表示电源有功出力的上下界约束、节点电压上下线约束、线路传输功率上下线约束等等。
电力系统最优潮流算法大致可以分为两类：
经典算法和智能算法。
其中经典算法主要是指以简化梯度法、牛顿法、内点法和解耦法为代表的基于线性规划和非线性规划以及解耦原则的算法，是研究最多的最优潮流算法,这类算法的特点是以一阶或二阶梯度作为寻找最优解的主要信息。
智能算法主要是指遗传算法和模拟退火发等，这类算法的特点是不以梯度作为寻优信息，属于非导数的优化方法。
因此经典算法的优点是能按目标函数的导数信息确定搜索方向，计算速度快，算法比较成熟，结果可信度高。
缺点是对目标函数及约束条件有一定的限制，可能出现局部极小时难以收敛。
而智能算法的优点是计算与导数无关，灵活性高，随机性强，缺点是算法不稳定，结果不可信，并且控制参数需凭经验给出。
通过对这些常见算法的简单比较，内点法具有其优越的性能，特别是路径跟踪法，其算法收敛迅速，鲁棒性强，对初值的选择不敏感，其迭代次数与系统规模或控制变量的数目关系不大，因此本文采用该方法进行最优计算。
1、路径跟踪法的基本数学模型
内点法最初的基本思路是希望通过寻优迭代过程始终在可行域内进行，因此，初始点应在可行域内，并在可行域的边界设置‘障碍’使迭代点接近边界时其目标函数迅速增大，从而保证迭代点均在可行域的内点。
但是对于大规模实际问题而言，寻找初始点往往十分困难。
为此许多学者长期以来致力于内点算法初始‘内点’条件的改进。
以下介绍的路径跟踪法只要求在寻优过程中松弛变量和拉格朗日乘子满足简单的大于0或者小于0的条件，即可代替原来必须在可行域内求解的要求，使得计算过程大为简化。
一般可以将最优潮流模型简化为如下的非线性优化模型。
Obj.（1-1）
s.t.（1-2）
（1-3）
其中为优化的目标函数，为等式约束，为不等式约束，路径跟踪内点法的基本思路是：
首先将式（1-3）的不等约束变成等式约束：
（1-4）
（1-5）
其中松弛变量，，应满足
u>0，l>0
这样原问题就转化为问题A：
Obj.
S.t.
然后，把目标函数改造成障碍函数，该函数在可行域内应接近于原函数f（x），而在边界时变得很大。
一次可得带优化问题B：
obj.
s.t.
其中扰动因子或者障碍因子u>0。
当l或u接近边界时，以上函数将趋于无穷大，因此满足以上障碍目标函数的极小解不可能在边界上找到。
这样就通过目标函数的变化把含不等式限制的优化问题A变成只含等式限制优化的问题B了，因此可以直接用拉格朗日乘子法来求解。
优化问题B的拉格朗日目标函数为：
式中：
和均为拉格朗日乘子。
因此最后简化的求解问题就是求取上述表达式的极小解。
2、路径跟踪法的最优潮流求解思路
路径跟踪法的最优潮流求解过程就是对拉格朗日目标函数求极小值问题：
式中：
和均为拉格朗日乘子。
该问题极小值存在的必要条件是拉格朗日函数对所有变量及乘子的偏导数为0。
即：
（2-1）
通过上述表达式可以解出：
（2-2）
定义：
，称为互补间隙。
可得：
（2-3）
如果x*是优化问题A的最优解，当u固定时，x（u）是优化问题B的解，那么当Gap0，u0时，产生的序列{x（u）}收敛至x*。
建议采用：
。
式中称为中心参数，一般取0.1，在大多数场合可获得较好的收敛效果。
通过偏导数为0的表达式可以可得内点法的修正方程为：
（2-4）
求解方程可得到第k次迭代的修正量，于是最优解的一个新的近似解为：
（2-5）
式中，和为步长：
（2-6）
其潮流计算的流程图如下图1所示，其中初始化部分包括：
（1）、设置松弛变量l和u，保证[l，u]T>0
（2）、设置拉格朗日乘子w、y、z，满足[w<0,z>0，Y！
=0]T
（3）、设置优化问题的初值。
（4）、取中心参数，给定计算精度，迭代次数初值K=0。
输出最优解，停止计算
计算互补间隙Gap
计算扰动因子
计算和
更新原始变量及拉格朗日乘子
kGap<
求解修正方程，求出，，，，，
输出“计算不收敛”
初始化
是
是
否
否
图1.内点法潮流计算流程图
3、具体算例及程序实现流程
这部分主要有算例描述以及程序的实现流程两部分，其中算例描述主要是对系统参数以及优化问题进行说明。
而程序的实现流程主要描述的是最优潮流计算中所涉及的矩阵方程的描述。
3.1、算例描述
该算例为王锡凡编写的《现代电力系统分析》中的3-1的例题，是以系统燃料最省为最优潮流的目标函数。
选择该题目作为算例分析的原因是，该题目有比较详细的解题思路以及列写出了比较详细的迭代结果，方便对编写程序的运行结果进行比对。
求如下图所示简化系统的系统燃料最省的最优潮流计算：
1
22
4
32
5
1:
1.05
1.05:
1
j0.015
j0.03
j0.25
j0.25
0.08+j0.30
0.1+j0.35
0.04+j0.25
j0.25
j0.25
3.7+j1.3
2+j1
1.6+j0.8
除了由上图所提供的系统母线负荷功率数据、线路参数和变压器之路参数数据、变压器便比数据之外，以下顺序给出了线路传输功率边界（表3-1），发电机有功无功出力上下界和燃料耗费曲线参数（表3-2）。
若不作特殊说明，所有数据都是以标幺值形式给出，功率基准值为100MVA，母线电压上下界分别为1.1和0.9。
表3-1线路传输功率边界
支路号
首末端母线号
线路传输功率边界
1
1-2
2
2
1-3
0.65
3
2-3
2
4
2-4
6
5
3-5
5
表3-2发电机数据
发电机序号
母线号
出力上界
出力下界
燃料耗费曲线参数
有功
无功
有功
无功
二次系数
一次系数
常数
1
4
8
3
1
-3
50.4395
200.4355
1200.64
2
5
8
5
1
-2.1
200.55
500.745
1857.20
3.2、程序具体实现过程
针对上述系统，首先我们先列写出该算例的数学模型和有关计算公式。
在该算例中，共有5个节点，相应的状态量为：
系统中有2台发电机，没有其他无功源，因此控制变量为：
应该指出，此处发电机和无功源的编号与及诶单编号无关，是独立编号的。
这是因为系统中一个节点可能接有多台发电机的缘故。
因此系统中总变量共有14个：
最优潮流的数学模型为：
目标函数：
约束条件：
每个节点有两个潮流方程，共有10个等式约束条件，对非发电机而言：
（i=1,2,3）
对发电机节点：
（i=4,5）
式中：
表示第k台发电机接在节点i上。
不等式约束共有14个条件，分别是：
根据以上模型可以形成修正方程。
该方程包括形成等式左边的系数矩阵和等式右边的常数项两部分。
1、形成系数矩阵
1）、等式约束的雅克比矩阵
等式右端包括3个子矩阵：
其中：
其中：
式中：
i为发电机的序号；j为节点号；表示第i台发电机是在节点j上的。
（潮流计算中的雅克比矩阵）
2）、不等式约束的雅克比矩阵
式中：
、、和依次表示电源有功出力的上下界约束，无功电源出力的上下界约束，节点电压赋值的上下界约束和线路潮流约束。
，，，
，，，
，
式中：
第行列元素为1，其他元素均为0。
3）、对角矩阵
4）、海森伯矩阵
这是最复杂的部分，共包含四项。
有上述推导已经可以得到其中的第四项为：
其余三项是：
目标函数的海森伯矩阵、等式约束海森伯矩阵与拉格朗日乘子y的乘积和不等式约束海森伯矩阵与拉格朗日乘子的乘积。
2、形成常数项
均可根据定义直接求得。
可以表示为：
当知道目标函数梯度矢量
之后，再根据以上等式和不等式约束的雅克比矩阵公式就可以求得。
4、运行结果及分析
4．1运行结果
以下对该算例的寻优过程用数字加以说明，设4、5节点发电机均能有算法调节其出力。
在初始化过程中各变量初值根据实际问题自行设置的，我们给出所用变量的处置如下：
节点电压；平衡节点；发电机出力有功取其边界值；松弛因子，当收敛条件时，需要迭代进行23次（例题所给出的迭代次数为17次）。
表4-1迭代过程中各节点电压增量的变化情况
迭代次数
1
-0.
-0.
-0.26638
-0.
-0.
2
0.
0.
0.
0.
-0.
3
2.
1.
3.
1.
3.
4
0.
-0.
2.
-0.
2.
5
0.
0.
-0.18991
0.
-0.
6
-0.
-0.
-0.31384
-0.
-0.
7
-0.
-0.
-0.52191
-0.
-0.
8
0.50336
-0.
1.
-0.
1.
9
0.
-0.05228
0.
-0.
0.
10
-0.
-0.
-0.13409
-0.
-0.
11
-0.
-0.
-0.135
-0.02272
-0.
12
0.
0.
0.
-0.
0.
13
0.
-0.
0.0966
-0.
0.
14
0.
-0.
0.
-0.
0.
15
-0.
-0.
-0.01655
1.24E-05
-0.
16
-0.
-0.
-0.01235
-0.
-0.
17
-0.
-0.
-0.00473
-0.
-0.
18
0.
-0.00223
0.
-0.
-0.
19
0.
-0.
0.
-0.
-0.
20
0.
-0.
0.
-0.
-0.
21
0.
-0.
0.
0.
-0.
22
4.42E-08
-7.14E-08
8.44E-06
1.71E-05
-6.67E-08
23
-9.43E-09
7.65E-09
1.38E-07
2.40E-08
1.58E-08
表4-2迭代过程中各节点相角增量的变化情况
迭代次数
1
10.
10.
9.
10.
9.
2
0.
-0.
0.
-0.
0.
3
0.
-0.
0.
-0.
0.
4
0.
0.
0.
0.
-0.
5
-28.46696
-28.48845
-28.5226
-28.58592
-28.
6
32.
35.
29.1472
35.
29.
7
-0.
-0.
-0.30025
-0.
-0.
8
-0.
-0.
-0.3624
-0.07659
-0.
9
-0.
0.
-0.27117
0.
-0.
10
6054.9484
-6300.974
6058.914
3059.0885
-121.
11
-0.
0.
-0.13752
0.
-0.
12
-131.1933
-133.1411
-129.232
-131.1213
-131.21147
13
0.
0.
0.
0.
0.
14
0.
-0.
-0.11121
0.
-0.
15
0.
0.03296
-0.04828
6.43E-02
-0.
16
0.
0.
0.
-0.
0.
17
-0.
-0.
-0.18778
-0.
-0.
18
-6.68934
-6.67987
-6.6504
-6.62907
-6.
19
-0.06977
-0.
-0.00065
-0.
-0.
20
-0.00417
-0.
-0.00209
-0.00978
-0.
21
0.973E-06
0.966E-06
0.97E-06
0.9718E-06
0.9797E-06
22
-1.42E-08
-1.43E-08
-1.43E-08
-1.43E-08
-1.41E-08
23
1.01E-09
1.01E-09
1.1E-08
1.01E-09
1.01E-09
表4-3迭代过程中有功源有功、无功源无功增量的变化情况
迭代次数
有功源有功出力增量
无功源无功出力增量
1
6.
0.
0.
2.
2
1.
0.
3.
-3.
3
0.
-0.
-11.7288
-1.
4
-0.
-0.
-15.2923
8.
5
0.
0.
3.
-1.
6
1.
0.
3.
0.
7
0.
0.
4.
-0.
8
-0.
0.
-6.52471
4.
9
0.
0.
-3.01366
3.
10
0.
0.
1.
0.
11
0.
0.
1.47339
0.
12
0.
-0.
-3.81381
3.
13
-0.
0.
-0.45866
0.
14
0.
-0.
-0.47294
1.
15
0.
-0.
0.
3.58E-01
16
0.
0.
0.
0.
17
-0.081
0.
-0.24799
0.
18
-0.
0.
-0.22899
-0.
19
-0.3602
0.
-0.22777
-0.
20
-0.
0.
0.
-0.
21
-0.
0.
0.
-0.
22
-1.41E-04
1.34E-04
2.33E-04
-1.17E-04
23
-2.60E-06
2.54E-06
3.21E-08
-1.69E-07
将各次迭代过程中Gap变化情况绘制成曲线，可以显示出路劲跟踪法最优潮流计算的收敛特性，如图4-1所示：
图4-15节点系统最优潮流内点法收敛特性
图4-2为5节点系统最优潮流计算结果截图，其中包括迭代次数、燃料总费用、发电机有功无功出力、各节点电压幅值与相角、以及各支路有功功率。
（注：
结果中的值为标幺值，功率的基准值为100MVA）
4．2结果分析
将最优潮流计算的结果和普通潮流计算结果进行比较，其中PF表示为普通潮流计算。
普通潮流计算中，发电机不会调节其出力。
即4节点为PQ节点，5节点为平衡节点。
见表4-4和表4-5。
从表中可以看出，由于4机组比5机组的燃料曲线系数小，因此4机组有功出力增加，5机组有功出力减少。
同时系统的网损、无功功率都有所增加。
这是由于要将1节点电压抬高至其下界以满足不等式约束的要求而产生的副作用。
但是网损的增加并不影响目标函数的优化。
整个系统的燃料费用与不优化的潮流计算相比仍然减少了243.76$。
表4-4各有功源有功和无功源无功出力
发电机序号
母线
序号
有功出力
无功出力
燃料费用/$
OPF
PF
OPF
PF
OPF
PF
1
4
5.5056
5.00
1.7780
1.88311
3833.06
3463.80
2
5
2.1568
2.5794
2.6194
2.2994
3870.13
4483.15
总计
7.6624
7.5794
4.3974
4.1305
7703.20
7946.95
表4-5各节点电压向量
母线序号
电压幅值
电压相角/rad
OPF
PF
OPF
PF
1
0.9000
0.8822
-0.007
-0.08340
2
1.1000
1.0779
0.4049
0.31160
3
1.0818
1.0364
-0.0571
-0.7473
4
1.0697
1.0500
0.47867
0.31160
5
1.100
1.0500
0
0
5、结论
路径跟踪法在电力系统中应用与求解线性规划和非线性规划模型中，还是比较有优势的，具有算法收敛迅速，鲁棒性强，对初值的选择不敏感，其迭代次数与系统规模或控制变量的数目关系不大等特点。
即在该简化模型中迭代次数为23次，但是由于其迭代次数是与系统规模关系不大，对IEEE30、IEEE118节点系统的计算结果其迭代次数始终保持在21到27次之间。
另外将水火电最优潮流问题分解为火电最优潮流子问题和水电子问题，提供了有效的协调算法。
但是对于路径跟踪法的影响因素还是比较多的，比如初始点的选择、迭代步长的选取、壁垒参数的调整、离散变量的处理等等，如果选取不当可能会出现不恰当的结果，因此还需要研究者们做大量的工作。
6、编程中遇到的问题
上图为例题所给出的迭代次数与仿真结果，可以看出除了迭代次数与书中不一致以外，迭代的结果基本上完全一样。
但在程序的实现过程中，很难按照书上的流程编写出结果一致的程序，因为程序中有大量的矩阵计算，而且也不能将书上列写的矩阵直接翻译成MATLAB语言，需要进行一些不同的处理方式，所以需要在网上去寻找一些算法的实现方式。
另外在编程中发现，王锡凡的电力系统一书中，也有一些公式书写有误的现象，所以需要对公式进行一定的验证推导，因此很大程度上会出现问题。
另外对于路径跟踪法来说，其初始点的选择、迭代步长的选取、壁垒参数的调整等都对计算结果又一定的影响。
而在本程序中，我选取的初始点以及拉格朗日因子也都与例题所提供的有些许不同，如果选择和书中相同的参数，计算结果就会出现问题，这其中的原因，我的猜测是可能在用matalb语言实现过程中，细节方面可能与例题所展示的有所出入，由于时间关系，暂时还未找到原因。
参考文献
[1]张伯明.高等电力网络分析[M].清华大学出版社,2007.
[2]王锡凡.现代电力系统分析[M].北京科学出版社,2003
[3]张江红.最优潮流算法综述[J].华北电力，2010,07
[4]赫玉国.一种基于KarmarKar内点法的最优潮流算法[J].中国机电工程学报，1996,11
附录
clc
clear
%%%%%%%%%%%%%%%%%%%%%%%读取计算参数%%%%%%%%%%%%%%%%
Branch=load（'Branch.txt'）;%支路参数文档
Node=load（'Node.txt'）;%节点参数文档
Generator=load（'Generator.txt'）;%发电机参数文档
%支路数据提取
Nbr=Branch（:
1）;%支路号
Nl=Branch（:
2）;%支路首节点
Nr=Branch（:
3）;%支路末节点
%节点数据提取
N=Node（:
1）;%节点号
Type=Node（:
2）;%节点类型
Uamp=Node（:
3）;%节点电压幅值
Dlta=Node（:
4）;%节点电压相角
Pd=Node（:
5）;%节点负荷有功
Qd=Node（:
6）;%节点负荷无功
Pg=Node（:
7）;%节点出力有功
Qg=Node（:
8）;%节点出力无功
%发电机数据提取
Ng=G
升级会员 