数学规划-基于Matlab符号计算工具箱的内点法最优潮流研究
的不敏感性。
#! 0#&（$）)# $ 0 ’(1
#" 0 3（#&（$）)# $ 0 ’)1）
#
)
0
3’3&
,#&（$）)#
$3
’3&（,’
(1
"
2’ *1
）
#
*
0
.
3&/#&（$）)#
$
2.
3&（/’
)
1
"
3’ +1
）
! "! " ! " 1（2）3（$）) #$ 03 $ （2，"）
3（$） 1 # (
"
’ * "’,-3" - 0 1
"
’ + "./- 2" - 0 1
（&:） （&4）
并在迭代过程中利用扰动变量和互补间隙的关系让
其逐步动态调整到零。
这样，上述扰动的 (() 条件就可用牛顿法解开，
修正方程式如下：
’&10 （#%%（$）( 2#%&（$）（)2*）3#% # （$）#$2
#%（$）# ( 2#&（$）（# )2#*）
%#9
$
（K）
+J" # +I" G ［!"（(% &"% F!% *"%）F ("（!% &"%/(% *"%）］
%#9
" E9，*，…，$
不等式约束 $（"）是系统的运行约束，分别为线
路约束、电压约束、发电机有功出力约束及无功电源
出力约束，即
00
)00 "%!)"%!) "%， （"，% ）$,"L
通过对式（&6）;（%&）进行适当的化简，可以得到
如 下 形 式 的 简 约 修 正 方 程 ，这 种 处 理 的 结 果 是 简 约
修 正 方 程 的 维 数 仅 取 决 于 等 式 约 束 的 维 数 ，与 不 等
式 约 束 的 维 数 无 关 ，有 效 地 保 证 了 算 法 对 系 统 规 模
& ’ &41
562 27
!"#$ 8：打印“迭代不收敛”，停机。
尽 管 路 径 跟 踪 内 点 算 法 逻 辑 清 晰 ，主 程 序 步 骤
简单，不需用户设定起始参数，且算法具有局部的二
次 收 敛 特 性 。 但 代 价 是 在 形 成 修 正 方 程 （9% ）时
（!"#$ :），需严格计算所有等式、不等式约束和目标
+,-.［,-.
（/!" !!"
：! !" 0 (；
/ #" !#"
：!#" 0 (）， 1］
&"#$2 ’
()***
+,-.［,-.
（/ $" ! $"
：! $" 0 (；
/ %" !%"
：! %" 3 (），
1］
!"#$ +：更新原始和对偶变量
!! $!! $ !! ! $!" $!" $ !! " $ ""# %%’""# %%4 &"#$$ ""! # %%，""% %%’""% %%4 &"#$2 ""! % %% #$ &#$ & #! $ &#& &#& & #! & &
制 次 数 /*?@ 0 41，中 心 参 数 %&（1，&），允 许 误差 & 0 &1A36。分别设置原始和对偶变量
的初值［!，"］)B 1，［) B 1，* C 1，( 0 1］)
DE+F=（/ C /*?@）GH </=> &：计算互补间隙 0#?>
0#?> " -)（’,3./）+I（0#?>C &），输出最优解，停机。否则 </=> %：由 0#?> 和 % 计算扰动参数 "
（&6）
式中
%’’1 ,#%（$）)# $
（&7）
%’(1,#&（$）)# $%# !
（&8）
ห้องสมุดไป่ตู้
%’)1,#&（$）)# $-# "
（&9）
"
%’ *1 , ,# ! 0 ’# )
（%1）
"
%’ +1 ,/# " - .# *
（%&）
#% %（$），#% &（$）分别是 %（$）和 &（$）海森矩
阵的线性组合。
! 引言
自法国学者 "#$%&’()&$ 在 *! 世纪 +! 年代初定 义了最优潮流 ,-.（,%()/#0 -12&$ .012）以来，作为 电 力 系 统 运 行 和 分 析 的 一 种 强 有 力 的 工 具 ，最 优 潮 流吸引了众多研究者的注意。在过去的四十年里， 各种不同的最优化方法，诸如：线性规划，二次规划， 非 线 性 规 划 以 及 解 耦 法 和 牛 顿 法 等 ，都 被 相 继 应 用 到这一领域。
（6）
式中 ’（5）0’（$，!，"；(，)，*，) ，*）；()，))，*) 为拉格 朗日乘子；); )，*; ) 为库恩乘子。
注意到有如下关系的成立：
!’（5）0 ) 3 ;)0 1， !’（5）0 3 *3*; 0 1
!!
!"
于是可以推导出式（4）的 (() 条件为
’&"##（$）+#%（$）( %#&（$）（)2*）0 1 （7）
简约修正方程（%:）的形成和求解上。
另外，为获得最佳的收敛特性，在扰动变量的 调
整 计 算 中 引 入 一 中 心 参 数 %，合 理 地 设 置 该 参 数将
使寻优过程在最优性和可行性之间找到一最佳的平
衡点，使算法具有良好的超线性收敛特性［’］。
路径跟踪法主程序流程为
</=> 1：初始化；设迭代标记单元 / 0 1，最大迭代限
%（$）
（%%） （%’） （%:）
其中
1（2）0［#% %（$）(2#% &（$）（)-*）%#%#（$）］2
#&（$）$#&（$）)， 4". 3&/3’3&,
3（$）0#%（$）)
$（5，"）0#%（$）( 3##（$）-#&（$）.
（. 3&/’)1 % ’3&,’(1 %"（. 3&%’3&）-） （%4） 可 见 ，路 径 跟 踪 法 的 主 要 计 算 负 荷 集 中 在 每 次
4#(0#= 作为一种强大的矩阵计算工具，内 建 丰 富高效的矩阵运算函数库。这使得在 4#(0#= 平台上 的 ,-. 程序开发极具效率。本文利用 4#(0#= 中的符 号计算工具箱 >4?（>@/=10)A 4#(5 ?110=1B）以符号 为变量完成了整个内点最优潮流的计算。可以获得 系统状态变量的显式符号结果，使得复杂的 ,-. 计 算仅限于在每次迭代中进行一次简单的代数替换， 大 大 简 化 了 ,-. 问 题 的 建 模 和 程 序 编 制 的 复 杂 程 度，提高 了 代 码 的 通 用 性 和 易 维 护 性 ，有 效 地 缩 短 了 ,-. 程序的开发周期。
)"%
#（( * "
/
!
* "
’(" (% ’
!"
!%）&"% /（("
!% ’(%
!"）*"%
1
* "%
!（(
*/
"
!
* "
）!1
* "
，
" # 9，*，…，$
（<）
00
)00 H"!)H"!) H"， "$,H
00
+00 J"!+J"!+ J"， "$,J
式中 -!"，-9"，-*" 为火电厂 " 的发电费用系数；)H"，+J"
’’" %（$）0 1 ’(" &（$）3! 35& 0 1 ’)" &（$）2" 3& 0 1 ’*" ’,- 0 1 ’+" ./- 0 1 （!，"；)）!1，*$1，(%1
（8） （9） （&1） （&&） （&%） （&’）
式中 ’，.，,，/ 是对角矩阵。
理论上可用牛顿 法 直 接 求 解 上 述 (() 条 件 而
（9）
%$!$（"）!$
其 中 ，目 标 函 数 !（"）分 别 设 为 有 功 网 损 最 小 ，无 功
网 损 最 小 、燃 料 费 用 最 小 和 固 定 负 荷 下 的 发 电 量 最
小 ，即
!#F"
$ #9
%
$ #9
潮流计算(matlab)实例计算
潮流例题：根据给定的参数或工程具体要求（如图），收集和查阅资料；学习相关软件（软件自选：本设计选择Matlab进行设计）。
2.在给定的电力网络上画出等值电路图。
3.运用计算机进行潮流计算。
4.编写设计说明书。
一、设计原理1．牛顿-拉夫逊原理牛顿迭代法是取x0 之后，在这个基础上，找到比x0 更接近的方程的跟，一步一步迭代，从而找到更接近方程根的近似跟。
牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0 的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复根。
电力系统潮流计算，一般来说，各个母线所供负荷的功率是已知的，各个节点电压是未知的（平衡节点外）可以根据网络结构形成节点导纳矩阵，然后由节点导纳矩阵列写功率方程，由于功率方程里功率是已知的，电压的幅值和相角是未知的，这样潮流计算的问题就转化为求解非线性方程组的问题了。
为了便于用迭代法解方程组，需要将上述功率方程改写成功率平衡方程，并对功率平衡方程求偏导，得出对应的雅可比矩阵，给未知节点赋电压初值，一般为额定电压，将初值带入功率平衡方程，得到功率不平衡量，这样由功率不平衡量、雅可比矩阵、节点电压不平衡量（未知的）构成了误差方程，解误差方程，得到节点电压不平衡量，节点电压加上节点电压不平衡量构成新的节点电压初值，将新的初值带入原来的功率平衡方程，并重新形成雅可比矩阵，然后计算新的电压不平衡量，这样不断迭代，不断修正，一般迭代三到五次就能收敛。
牛顿—拉夫逊迭代法的一般步骤：（1）形成各节点导纳矩阵Y。
（2）设个节点电压的初始值U和相角初始值e 还有迭代次数初值为0。
（3）计算各个节点的功率不平衡量。
（4）根据收敛条件判断是否满足，若不满足则向下进行。
（5）计算雅可比矩阵中的各元素。
（6）修正方程式个节点电压（7）利用新值自第（3）步开始进入下一次迭代，直至达到精度退出循环。
（8）计算平衡节点输出功率和各线路功率2．网络节点的优化１）静态地按最少出线支路数编号这种方法由称为静态优化法。
基于内点法最优潮流计算
1.1246e-001
4
-1.2326e-002
-1.8264e-001
1.9823e-001
-2.0804e-002
-1.9440e-002
5.0985e-002
5
-3.8403e-004
-7.6535e-002
7.7332e-002
-5.7025e-002
-2.2982e-003
5.3726e-002
基于内点法最优潮流计算
1
主要内容
1、 课题研究旳意义和现状 2、 最优潮流旳原对偶内点算法 3、 最优潮流旳预测校正内点算法 4、 结论
2
一、课题研究旳意义和现状
概念：
最优潮流问题（OPF）就是在系统构造参数及负荷 给定旳情况下，经过优选控制变量，拟定能满足全 部旳指定约束条件，并使系统旳某个性能指标到达 最优时旳潮流分布。
5 0.7754 0.9323
11 0.9995 0.9995
6 0.9716 0.9838
12 0.9995 0.9995
18
收敛特征分析
下表为计算过程中5节点系统旳迭代步长：
迭代次
数
1
2
3
4
5
6
7
8
ap
0.5222 0.0016 0.4027 0.3444 0.0016 0.9365 0.3377 0.6978
三个系统旳迭代次数分别为16、11、12次，迭代次数较少，计算时
间短，收敛特征好。
系统规模扩大时，迭代次数不会明显增长，阐明算法对系统规模不
敏感。
初始点为非内点时，算法也能够收敛至最优解，阐明算法对初始点
不敏感。
17
收敛特征分析
matlab潮流计算
前言电力系统潮流计算是研究电力系统稳态运行情况的一种计算，它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态：各母线的电压，各元件中流过的功率，系统的功率损耗等等。
在电力系统规划的设计和现有电力系统运行方式的研究中，都需要利用潮流计算来定量地分析比较供电方案或运行方式的合理性。
可靠性和经济性。
此外，电力系统潮流计算也是计算系统动态稳定和静态稳定的基础。
所以潮流计算是研究电力系统的一种很重要和很基础的计算。
随着科学技术的发展，电力系统变得越来越复杂，电气工程师掌握一种好的能对电力系统进行仿真的软件是学习和研究的需要。
文章简要介绍了MATLAB发展历史、组成和强大的功能，并用简单例子分别就编程和仿真两方面分析了MATIAB软件在电力系统研究中的具体应用。
采取等效电路法，能对特殊、复杂地电力系统进行高效仿真研究，因此，掌握编程和仿真是学好MATLAB的基础。
与众多专门的电力系统仿真软件相比，MATLAB软件具有易学、功能强大和开放性好，是电力系统仿真研究的有力工具。
1电力系统的基本概念电力系统：发电机把机械能转化为电能，电能经变压器和电力线路输送并分配到用户，在那里经电动机、电炉和电灯等设备又将电能转化为机械能、热能和光能等。
这些生产、变换、输送、分配、消费电能的发电机、变压器、变换器、电力线路及各种用电设备等联系在一起组成的统一整体称为电力系统。
电力网：电力系统中除发电机和用电设备外的部分。
动力系统：电力系统和“动力部分”的总和。
2潮流计算2.1潮流计算概述与发展电力系统潮流计算也分为离线计算和在线计算两种，前者主要用于系统规划设计和安排系统的运行方式，后者则用于正在运行系统的经常监视及实时控制。
利用电子数字计算机进行电力系统潮流计算从50年代中期就已经开始。
在这20年内，潮流计算曾采用了各种不同的方法，这些方法的发展主要围绕着对潮流计算的一些基本要求进行的。
对潮流计算的要求可以归纳为下面几点：1）计算方法的可靠性或收敛性；2）对计算机内存量的要求；3）计算速度；4）计算的方便性和灵活性。
matlab数学工具箱介绍
Matlab符号数学工具箱应用简介Matlab符号运算是通过集成在Matlab中的符号数学工具箱（Symbolic Math Toolbox）来实现的。
和别的工具箱有所不同，该工具箱不是基于矩阵的数值分析，而是使用字符串来进行符号分析与运算。
实际上，Matlab中的符号数学工具箱是建立在Maple基础上的，当进行Matlab符号运算时，它就请求Maple软件去计算并将结果返回给Matlab。
Matlab的符号数学工具箱可以完成几乎所有得符号运算功能。
这些功能主要包括：符号表达式的运算，符号表达式的复合、化简，符号矩阵的运算，符号微积分、符号函数画图，符号代数方程求解，符号微分方程求解等。
此外，工具箱还支持可变精度运算，既支持符号运算并以指定的精度返回结果。
在一般的Matlab书籍中都会对Matlab的符号运算做一些介绍，本文将略去这些简单的部分，主要对比较复杂的部分做一些介绍，另外，限于篇幅，和前面几篇一样，在此也仅仅列出函数的名称和功能，至于其参数设置，可借助Matlab的帮助系统一、符号表达式的运算[n,d]=numden(a) 提取符号表达式a的分子和分母，并将其存放在n和d中n=numden(a) 提取符号表达式a的分子和分母，只将分子存放在n中symadd(a,b) 返回符号表达式a和b的和，也可直接用a+bsymsub(a,b) 返回符号表达式a和b的差，也可直接用a-bsymmul(a,b) 返回符号表达式a和b的积，也可直接用a*bsymdiv(a,b) 返回符号表达式a和b的商，也可直接用a/bsympow(a,b) 返回符号表达式a的b次幂，也可直接用a^bcompose(f,g) 返回复合函数f(g(y))compose(f,g,z) 返回自变量为z的复合函数f(g(z))compose(f,g,x,z) 返回复合函数f(g(z))，并使x成为f函数的独立变量。
即，如果f=cos(x/t)，则compose(f,g,x,z)返回复合函数cos(g(z)/t)，而compose(f,g,t,z)返回cos(x/g(z))compose(f,g,x,y,z) 返回复合函数f(g(z))，并且使x与y分别成为f与g函数的独立变量。
基于MATLAB的电力系统潮流计算_毕业设计论文
基于MATLAB的电⼒系统潮流计算_毕业设计论⽂基于MATLAB的电⼒系统潮流计算摘要潮流计算是电⼒系统最基本、最常⽤的计算。
根据系统给定的运⾏条件、⽹络接线及元件参数，通过潮流计算可以确定各母线的电压（幅值及相⾓），各元件中流过的功率、整个系统的功率损耗等。
潮流计算是实现电⼒系统安全经济发供电的必要⼿段和重要⼯作环节。
因此潮流计算在电⼒系统的规划设计、⽣产运⾏、调度管理及科学研究中都有着⼴泛的应⽤。
本次设计的主要⽬的就是⾯向⼀般的电⼒⽹络，形成节点导纳矩阵，确定合适的算法，编写通⽤的计算程序，得到计算结果。
设计中主要介绍了⽜顿拉夫逊和PQ分解两种算法，PQ分解法虽然在结构上⽐⽜顿法更加简化，但是针对⼀般⽹络现代计算机在存储空间及计算速度上已经⼗分强⼤，鉴于对⽜顿法的熟悉与其算法的直观性，本次设计在编程时采⽤了⽜顿拉夫逊法的直⾓坐标形式。
解⽅程的过程利⽤Matlab的强⼤计算功能，编写M语⾔，合理设置变量，实现通⽤计算功能。
关键词: 电⼒系统，潮流计算，⽜顿—拉夫逊法，Matlab。
AbstractPower system load flow calculation is the most basic and commonly used calculations. Given according to the system operating conditions, the network connection and device parameters can be determined by power flow calculation of the bus voltage (magnitude and phase angle), the power flowing through the components, overall system power consumption and so on. Flow calculation is to achieve economic development of power system supply the necessary means and important part of the work. Therefore flow calculation in power system planning and design, production and operation, scheduling management, and scientific research have a wide range of applications.The main purpose of this design is for the general electricity network, the formation of the node admittance matrix, determine the appropriate method, the preparation of general-purpose computer program to get results. Introduces the design and the PQ decomposition Newton Raphson two algorithms, PQ decomposition although the structure is more streamlined than the Newton method, but for the general network of modern computer storage space and computing speed has been very strong, in view of the Newton Familiar with its intuitive algorithm, this design in programming using Newton Raphson polar form. The process of solving equations using matlab powerful computing capabilities, the preparation of M language, a reasonable set variables, to achieve general-purpose computing functions.Keywords: power system, power flow calculation, Newton - Raphson method, Matlab.⽬录摘要 (I)Abstract ................................................................................................................................................ II ⽬录.................................................................................................................................................... I II 1 引⾔ .. (1)1.1 潮流计算⽬的 (1)1.2 潮流计算意义 (1)1.3 潮流计算发展史 (1)1.4基于MATLAB 的电⼒系统潮流计算发展前景 (2)2简单电⼒系统潮流计算 (4)2.1简单辐射⽹络的潮流计算 (4)2.1.1简单⽀路的潮流分布和电压降落 (4)2.1.2 辐射型⽹络的⼿⼯潮流计算⽅法 (6)2.2 简单环⽹的潮流计算 (7)2.2.1两端电压相等 (7)2.2.2两端电压不相等 (8)3 复杂电⼒系统潮流计算的计算机算法 (10)3.1电⼒⽹络⽅程及等值电路 (10)3.2节点导纳矩阵形成及修改 (11)3.3节点的分类 (14)3.3.1 PQ节点 (14)3.3.2 PV节点 (14)3.3.3 平衡节点 (14)3.4潮流计算的约束条件 (15)3.5⽜顿-拉夫逊法（直⾓坐标） (15)3.5.1⽜顿-拉夫逊法的推导过程 (15)3.5.2潮流计算时的修正⽅程（直⾓坐标） (17) 3.5.3雅可⽐矩阵的特点: (19)3.5.4⽜顿-拉夫逊法计算步骤 (19)3.6 P-Q分解法潮流计算 (20)3.6.1 P-Q分解法潮流计算概述 (20)3.6.2 P-Q分解法的潮流计算步骤 (20)3.6.3 P-Q分解法的特点 (21)4 Matlab概述 (22)4.1Matlab简介 (22)4.2 Matlab中的变量 (22)4.3 Matlab编程 (23)4.3.1矩阵的输⼊ (23)4.3.2矩阵的运算 (24)4.3.3 MatLab的控制流 (24)5 ⽜顿法潮流计算程序设计及实例 (26)5.1⼿算 (26)5.2计算机算法的数据输⼊ (29)5.3潮流计算程序 (30)5.3 计算结果分析 (36)结论 (37)参考⽂献 (38)附录A 程序流程图 (39)附录B Matlab仿真 (40)致谢 (1)1 引⾔1.1 潮流计算⽬的电⼒系统潮流计算是研究电⼒系统稳态运⾏情况的⼀种基本电⽓计算。
matlab的dace工具箱加点准则
matlab的dace工具箱加点准则MATLAB是广泛应用于各个领域的数学软件，它拥有着丰富的工具箱，其中DACE工具箱是其中的一个。
DACE（Design and Analysis of Computer Experiments）工具箱是用于建立计算机实验的数据分析和建模的工具箱。
DACE工具箱主要用于优化和建模工作，而加准则则是在DACE工具箱中常用的方法之一。
在DACE工具箱中，加点准则是一种优化方法，它通过为模型选择适当的样本点来最小化模型误差，并提高模型预测的准确性。
加点准则主要涉及到以下几个步骤：1. 定义矩阵定义矩阵是建立DACE模型之前的第一步。
该矩阵包括了实验设计的输入和输出值，可以通过给定的随机设计或者Latin-hypercube采样方法进行生成。
2. 设计DACE模型设计DACE模型的目的是建立起原始数据之间的映射关系，以便将任意新的输入变量映射到相应的输出变量。
DACE建立的是一种高斯过程模型（gussian process model），并具有许多优点，例如可以利用小样本数据进行建模，映射函数的预测误差较小等。
3. 计算加点准则在DACE模型基础上，通过加点准则，选取出具有代表性的样本点，并利用这些样本点来重新构建模型从而提高模型的预测精度。
加点准则主要有几种常见的计算方式：A. 最小边际方差准则（Minimum Variance Criterion，MVC），它能够选择对预测输出变量贡献最大的样本点。
B. 最大自相关准则（Maximum Correlation Criterion，MXC），它能够选择对模型贡献最大的样本点，从而协同计算挑选哪些样本点最能提高模型质量。
C. 最大熵准则（Maximum Entropy Criterion，MEC），适用于大数据样本，充分利用样本点之间的信息，计算出样本点之间的相对独立度，从而确定选取哪些样本点有可能提高模型的精度。
4. 计算预测值在得到加点准则选出的样本点后，计算预测值是最后一个步骤。
内点法最优潮流MATLAB算法
内点法最优潮流MATLAB算法clear;%clc;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%数据加载n=input('请输入要计算的节点系统(5):')load Node5.txt;%节点数据load Branch5.txt;%支路数据load Generator5.txt;%发电机数据Node=Node5;Branch=Branch5;Generator=Generator5;%节点数据处理N=Node(:,1);%节点号Type=Node(:,2);%节点类型Uamp=Node(:,3);%节点电压幅值Dlta=Node(:,4);%节点电压相角Pd=Node(:,5);%节点负荷有功Qd=Node(:,6);%节点负荷无功Pg=Node(:,7);%节点出力有功Qg=Node(:,8);%节点出力无功Umax=Node(:,9);%节点电压幅值上限 Umin=Node(:,10);%节点电压幅值下限Bc=Node(:,11);%节点补偿电容电纳值 %支路数据处理Nbr=Branch(:,1);%支路号Nl=Branch(:,2);%支路首节点Nr=Branch(:,3);%支路末节点R=Branch(:,4);%支路电阻X=Branch(:,5);%支路电抗Z=R+1i*X;%支路阻抗=支路电阻+支路电抗 Bn=Branch(:,6);%支路对地电纳K=Branch(:,7);%支路变压器变比，0表示无变压器 Ptmax=Branch(:,8);%线路传输功率上限%发电机数据处理Ng=Generator(:,1);%发电机序号Nbus=Generator(:,2);%所在母线号Pumax=Generator(:,3);%发电机有功出力上界 Qumax=Generator(:,4);%发电机无功出力上界 Pumin=Generator(:,5);%发电机有功出力下界Qumin=Generator(:,6);%发电机无功出力下界a2=Generator(:,7);%燃料耗费曲线二次系数a1=Generator(:,8);%燃料耗费曲线一次系数a0=Generator(:,9);%燃料耗费曲线常数项%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%n=length(N);%节点个数ng=length(Ng);%发电机台数nbr=length(Nbr);%支路个数x=zeros(2*(ng+n),1);%控制变量+状态变量x(1:ng)=Pg(Nbus);x(ng+1:2*ng)=Qg(Nbus);x((2*ng+2):2:2*(ng+n))=Uamp; x((2*ng+1):2:2*(ng+n)-1)=Dlta; l=0.8*ones(2*ng+n+nbr,1);%松弛变量u=1.1*ones(2*ng+n+nbr,1);%松弛变量w=-1.5*ones(2*ng+n+nbr,1);%拉格朗日乘子z=ones(2*ng+n+nbr,1);%拉格朗日乘子y=zeros(2*n,1);%拉格朗日乘子y(1:2:2*n-1)=1e-3;y(2:2:2*n)=-1e-3;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算不等式约束的上下限%%%%%%%%%%%%%%%%%%%%%%%%%gmingmin=zeros(2*ng+n+nbr,1);gmin(1:ng)=Pumin;gmin(ng+1:2*ng)=Qumin;gmin(2*ng+1:2*ng+n)=Umin;gmin(2*ng+n+1:2*ng+n+nbr)=-Ptmax; %gmaxgmax=zeros(2*ng+n+nbr,1);gmax(1:ng)=Pumax;gmax(ng+1:2*ng)=Qumax;gmax(2*ng+1:2*ng+n)=Umax;gmax(2*ng+n+1:2*ng+n+nbr)=Ptmax;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成导纳矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Y=zeros(n,n);%%%%%%%%%%%%%%%%%%%%计算非对角元素%%%%%%%%%%%%%%%%%%%%% for ii=1:nbr if K(ii)==0%非变压器支路Y(Nl(ii),Nr(ii))=-1/Z(ii);Y(Nr(ii),Nl(ii))=Y(Nl(ii),Nr(ii));else%变压器支路Y(Nl(ii),Nr(ii))=-1/Z(ii)/K(ii);Y(Nr(ii),Nl(ii))= Y(Nl(ii),Nr(ii));endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%计算对角元素%%%%%%%%%%%%%%%%%%%%%%for ii=1:n%将支路导纳加入到对角元素中for jj=1:nbrif K(jj)==0&&(Nl(jj)==ii||Nr(jj)==ii)%非变压器支路Y(ii,ii)=Y(ii,ii)+1/Z(jj);else if K(jj)~=0&&(Nl(jj)==ii||Nr(jj)==ii)%变压器支路Y(ii,ii)=Y(ii,ii)+1/Z(jj)/K(jj);endendendendfor ii=1:nbr%将对地电纳加入到对角元素中if K(ii)==0%非变压器支路Y(Nl(ii),Nl(ii))=Y(Nl(ii),Nl(ii))+1i*Bn(ii);Y(Nr(ii),Nr(ii))=Y(Nr(ii),Nr(ii))+1i*Bn(ii);else%变压器支路Y(Nr(ii),Nr(ii))=Y(Nr(ii),Nr(ii))+(K(ii)-1)/K(ii)/Z(ii);Y(Nl(ii),Nl(ii))=Y(Nl(ii),Nl(ii))+(1-K(ii))/K(ii)/K(ii)/Z(ii);endendfor ii=1:nY(ii,ii)=Y(ii,ii)+i*Bc(ii);endG=real(Y);%电导B=imag(Y);%电纳%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%k=0;%迭代次数Kmax=150;%最大迭代次数iteration=1e-4;%误差精度delta=0.08;Gap=(l'*z-u'*w)*ones(Kmax,1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%主程序%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% while k<50%计算互补间隙GapGap(k+1)=l'*z-u'*w;if Gap>iterationmiu=delta*Gap(k+1)/(2*(2*ng+n+nbr)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%形成系数矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%相角差计算%%%%%%%%%%%%%%%%%%%%%%theta=zeros(n,n);for ii=1:nfor jj=1:ntheta(ii,jj)=Dlta(ii)-Dlta(jj);endend %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%1、等式约束雅克比矩阵%%%%%%%%%%%%%%%%pxh=zeros(2*(ng+n),2*n); %%%%%%%%%%%%%%%%%%%%%%%ah/aP%%%%%%%%%%%%%%%%%%%%%%%for ii=1:ngpxh(Ng(ii),2*Nbus(ii)-1)=1;end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%ah/aQ%%%%%%%%%%%%%%%%%%%%%%%for ii=1:ngpxh(Ng(ii)+ng,2*Nbus(ii))=1;end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%ah/ax%%%%%%%%%%%%%%%%%%%%%%%HH=zeros(n,n);JJ=zeros(n,n);NN=zeros(n,n);LL=zeros(n,n);for ii=1:nfor jj=1:nif ii~=jj%i!=j时的情况%非对角元素HH(ii,jj)=-Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));JJ(ii,jj)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin (theta(ii,jj)));NN(ii,jj)=-Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));LL(ii,jj)=-Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角元素HH(ii,ii)=HH(ii,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));JJ(ii,ii)=JJ(ii,ii)-Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)) );NN(ii,ii)=NN(ii,ii)-Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));LL(ii,ii)=LL(ii,ii)-Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));endendNN(ii,ii)=NN(ii,ii)-2*Uamp(ii)*G(ii,ii);LL(ii,ii)=LL(ii,ii)+2*Uamp(ii)*B(ii,ii);endpxh(1+2*ng:2:2*(n+ng)-1,1:2:2*n-1)=HH';pxh(1+2*ng:2:2*(n+ng)-1,2:2:2*n)=JJ';pxh(2+2*ng:2:2*(n+ng),1:2:2*n-1)=NN';pxh(2+2*ng:2:2*(n+ng),2:2:2*n)=LL';%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2、不等式约束的雅克比矩阵%%%%%%%%%%%%%%%%%%%% %g1:电源有功出力上下限约束ag1aP=eye(ng,ng);ag1aQ=zeros(ng,ng);ag1ax=zeros(2*n,ng);%g2:电源无功出力上下限约束ag2aP=zeros(ng,ng);ag2aQ=eye(ng,ng);ag2ax=zeros(2*n,ng);%g3:节点电压幅值上下限约束ag3aP=zeros(ng,n);ag3aQ=zeros(ng,n);ag3ax=zeros(2*n,n);for ii=1:nag3ax(2*ii,ii)=1;end%g4:线路潮流上下限约束ag4aP=zeros(ng,nbr);ag4aQ=zeros(ng,nbr);ag4ax=zeros(2*n,nbr);for ii=1:nfor jj=1:nbrif Nl(jj)==iiag4ax(2*ii-1,jj)=-Uamp(Nl(jj))*Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*sin(theta(Nl(jj),N r(jj)))-B(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj))));ag4ax(2*ii,jj)=Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj )))+B(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj))))-2*Uamp(Nl(jj))*G(Nl(jj),Nr(jj));endif Nr(jj)==iiag4ax(2*ii-1,jj)=Uamp(Nl(jj))*Uamp(Nr(jj))*(G(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr (jj)))-B(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj))));ag4ax(2*ii,jj)=Uamp(Nl(jj))*(G(Nl(jj),Nr(jj))*cos(theta(Nl(jj),Nr(jj )))+B(Nl(jj),Nr(jj))*sin(theta(Nl(jj),Nr(jj))));endendendpxg=[ag1aP ag2aP ag3aP ag4aP;ag1aQ ag2aQ ag3aQ ag4aQ;ag1ax ag2ax ag3ax ag4ax];%此即为不等式约束的雅克比矩阵%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%3、对角矩阵%%%%%%%%%%%%%%%%%%%%%%%% L_1Z=zeros(2*ng+n+nbr,2*ng+n+nbr);U_1W=zeros(2*ng+n+nbr,2*ng+n+nbr);for ii=1:2*ng+n+nbrL_1Z(ii,ii)=z(ii)/l(ii);U_1W(ii,ii)=w(ii)/u(ii);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%海森伯矩阵%%%%%%%%%%%%%%%%%%%%%%%%%% %将海森伯矩阵分为4块:H1,H2,H3,H4%%%%%%%%%%%%%%%%%%%%%H1%%%%%%%%%%%%%%%%%%%%%%A2=diag(a2);H1=zeros(2*(ng+n),2*(ng+n));H1(1:ng,1:ng)=2*A2;%%%%%%%%%%%%%%%%%%%%H2%%%%%%%%%%%%%%%%%%%%%%H2=zeros(2*(ng+n),2*(ng+n));A=zeros(2*n,2*n);Apb=zeros(2*n,2*n,n);Aqb=zeros(2*n,2*n,n);for ii=1:nfor jj=1:n %元素位置为:1 2if ii~=jj % 3 4%对角线上与ii对应的元素%ApApb(2*ii-1,2*ii-1,ii)=Apb(2*ii-1,2*ii-1,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(i i,jj)));%对角线处1号元素Apb(2*ii-1,2*ii,ii)=Apb(2*ii-1,2*ii,ii)+Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角线处2号元素%%3号元素与之相等%AqAqb(2*ii-1,2*ii-1,ii)=Aqb(2*ii-1,2*ii-1,ii)+Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角线处1号元素Aqb(2*ii-1,2*ii,ii)=Aqb(2*ii-1,2*ii,ii)-Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%对角线处2号元素%%3号元素与之相等%对角线上与jj对应的元素%ApApb(2*jj-1,2*jj-1,ii)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(i i,jj)));%对角线处1号元素Apb(2*jj-1,2*jj,ii)=-Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj))); %对角线处2号元素Apb(2*jj,2*jj-1,ii)=Apb(2*jj-1,2*jj,ii);%3号元素与2号元素相等%AqAqb(2*jj-1,2*jj-1,ii)=Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%对角线处1号元素Aqb(2*jj-1,2*jj,ii)=Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj ))); %对角线处2号元素Aqb(2*jj,2*jj-1,ii)=Aqb(2*jj-1,2*jj,ii);%3号元素与2号元素相等%4号元素为0%非对角线行元素%ApApb(2*ii-1,2*jj-1,ii)=-Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)) );%非对角线行处1号元素Apb(2*ii-1,2*jj,ii)=Uamp(ii)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处2号元素Apb(2*ii,2*jj-1,ii)=-Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处3号元素Apb(2*ii,2*jj,ii)=-(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处4号元素%AqAqb(2*ii-1,2*jj-1,ii)=-Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处1号元素Aqb(2*ii-1,2*jj,ii)=-Uamp(ii)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处2号元素Aqb(2*ii,2*jj-1,ii)=Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(theta(ii,jj)));%非对角线行处3号元素Aqb(2*ii,2*jj,ii)=-(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));%非对角线行处4号元素%非对角线列元素%ApApb(2*jj-1,2*ii-1,ii)=Apb(2*ii-1,2*jj-1,ii);%非对角线列处1号元素Apb(2*jj-1,2*ii,ii)=Apb(2*ii,2*jj-1,ii);%非对角线列处2号元素Apb(2*jj,2*ii-1,ii)=Apb(2*ii-1,2*jj,ii);%非对角线列处3号元素Apb(2*jj,2*ii,ii)=Apb(2*ii,2*jj,ii);%%非对角线列处4号元素%AqAqb(2*jj-1,2*ii-1,ii)=Aqb(2*ii-1,2*jj-1,ii);%非对角线列处1号元素Aqb(2*jj-1,2*ii,ii)=Aqb(2*ii,2*jj-1,ii);%非对角线列处2号元素Aqb(2*jj,2*ii-1,ii)=Aqb(2*ii-1,2*jj,ii);%非对角线列处3号元素Aqb(2*jj,2*ii,ii)=Aqb(2*ii,2*jj,ii);%%非对角线列处4号元素endend%对角线上与ii对应的元素%ApApb(2*ii,2*ii-1,ii)=Apb(2*ii-1,2*ii,ii);%对角线处3号元素与2号元素相等Apb(2*ii,2*ii,ii)=-2*G(ii,ii);%对角线处4号元素%AqAqb(2*ii,2*ii-1,ii)=Aqb(2*ii-1,2*ii,ii);%对角线处3号元素与2号元素相等Aqb(2*ii,2*ii,ii)=2*B(ii,ii);%对角线处4号元素endfor ii=1:nA=A+Apb(:,:,ii)*y(2*ii-1)+Aqb(:,:,ii)*y(2*ii);endH2(2*ng+1:2*(ng+n),2*ng+1:2*(ng+n))=A;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%H3%%%%%%%%%%%%%%%%%%%%%%H3=zeros(2*(ng+n),2*(ng+n));A3=zeros(2*n,2*n);Apc=zeros(2*n,2*n,nbr);for ii=1:nbr%对角线上iiApc(2*Nl(ii)-1,2*Nl(ii)-1,ii)=-Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))+B( Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii))));Apc(2*Nl(ii)-1,2*Nl(ii),ii)=-Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii))));Apc(2*Nl(ii),2*Nl(ii)-1,ii)=Apc(2*Nl(ii)-1,2*Nl(ii),ii);Apc(2*Nl(ii),2*Nl(ii),ii)=-2*G(Nl(ii),Nr(ii));%对角线上jjApc(2*Nr(ii)-1,2*Nr(ii)-1,ii)=-Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii)))+B( Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii))));Apc(2*Nr(ii)-1,2*Nr(ii),ii)=Uamp(Nl(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii))));Apc(2*Nr(ii),2*Nr(ii)-1,ii)=Apc(2*Nr(ii)-1,2*Nr(ii),ii);Apc(2*Nr(ii),2*Nr(ii),ii)=0;%非对角线ijApc(2*Nl(ii)-1,2*Nr(ii)-1,ii)=Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii )))+B(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii))));Apc(2*Nl(ii)-1,2*Nr(ii),ii)=-Uamp(Nl(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii))));Apc(2*Nl(ii),2*Nr(ii)-1,ii)=Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)))-B(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii))));Apc(2*Nl(ii),2*Nr(ii),ii)=G(Nl(ii),Nr(ii))*cos(theta(Nl(ii),Nr(ii))) +B(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii)));%非对角线jiApc(2*Nr(ii)-1,2*Nl(ii)-1,ii)=Apc(2*Nl(ii)-1,2*Nr(ii)-1,ii);Apc(2*Nr(ii)-1,2*Nl(ii),ii)=Apc(2*Nl(ii),2*Nr(ii)-1,ii);Apc(2*Nr(ii),2*Nl(ii)-1,ii)=Apc(2*Nl(ii)-1,2*Nr(ii),ii);Apc(2*Nr(ii),2*Nl(ii),ii)=Apc(2*Nl(ii),2*Nr(ii),ii);%求和c=z+w;A3=A3+Apc(:,:,ii)*c(2*ng+n+ii);endH3(2*ng+1:2*(ng+n),2*ng+1:2*(ng+n))=A3;%%%%%%%%%%%%%%%%%%%%%%%H4%%%%%%%%%%%%%%%%%%%%%%%%%H4=pxg*(L_1Z-U_1W)*pxg';%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%H=-H1+H2+H3-H4;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%形成常数项%%%%%%%%%%%%%%%%%%%%%%%%% %Lyh=zeros(2*n,1);for ii=1:nh(2*ii-1)=Pg(ii)-Pd(ii);h(2*ii)=Qg(ii)-Qd(ii);for jj=1:nh(2*ii-1)=h(2*ii-1)-Uamp(ii)*Uamp(jj)*(G(ii,jj)*cos(theta(ii,jj))+B(ii,jj)*sin(t heta(ii,jj)));h(2*ii)=h(2*ii)-Uamp(ii)*Uamp(jj)*(G(ii,jj)*sin(theta(ii,jj))-B(ii,jj)*cos(theta(ii,jj)));endendLy=h;%Lz%g(x)gx=zeros(2*ng+n+nbr,1);gx(1:ng)=x(1:ng);gx(ng+1:2*ng)=x(ng+1:2*ng);gx(2*ng+1:2*ng+n)=x(2*ng+2:2:2*(ng+n));for ii=1:nbrgx(2*ng+n+ii)=Uamp(Nl(ii))*Uamp(Nr(ii))*(G(Nl(ii),Nr(ii))*cos(theta( Nl(ii),Nr(ii)))+B(Nl(ii),Nr(ii))*sin(theta(Nl(ii),Nr(ii))))-Uamp(Nl(ii))*Uamp(Nl(ii))*G(Nl(ii),Nr(ii));endLz=gx-l-gmin;%LwLw=gx+u-gmax;%Lle=ones(2*ng+n+nbr,1);LZ=zeros(2*ng+n+nbr,2*ng+n+nbr);for ii=1:2*ng+n+nbr;LZ(ii,ii)=l(ii)*z(ii);endLl=LZ*e-miu*e;%LuUW=zeros(2*ng+n+nbr,2*ng+n+nbr);for ii=1:2*ng+n+nbrUW(ii,ii)=u(ii)*w(ii);endLu=UW*e+miu*e;%Lx'Lx1=zeros(2*(ng+n),1);Lx1(1:ng)=2*a2.*x(1:ng)+a1;Lx2=pxh*y;Lx3=pxg*c;Lx41=zeros(2*(ng+n),1);Lx42=zeros(2*(ng+n),1);for ii=1:2*ng+n+nbrLx41(ii)=(Ll(ii)+z(ii)*Lz(ii))/l(ii);Lx42(ii)=(Lu(ii)-w(ii)*Lw(ii))/u(ii);endLx4=pxg*(Lx41+Lx42);Lx=Lx1-Lx2-Lx3;Lxx=Lx+Lx4; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%求出修正量%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %dx,dyHxy=[H pxh;pxh' zeros(2*n,2*n)];LxLy=[Lxx;-Ly];I=eye(2*(ng+n)+2*n);dxdy=I/Hxy*LxLy;dx=dxdy(1:2*(ng+n));dy=dxdy(2*(ng+n)+1:2*(ng+n)+2*n);%dldl=pxg'*dx+Lz;%dudu=-pxg'*dx-Lw;%dzdz=zeros(2*ng+n+nbr,1);for ii=1:2*ng+n+nbrdz(ii)=(-Ll(ii)-z(ii)*dl(ii))/l(ii);end%dwdw=zeros(2*ng+n+nbr,1);for ii=1:2*ng+n+nbrdw(ii)=(-Lu(ii)-w(ii)*du(ii))/u(ii);end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%计算alfap和alfad%%%%%%%%%%%%%%%%%%%%%%%% alfap=1;alfad=1;for ii=1:2*ng+n+nbrif dl(ii)<0&&-l(ii)/dl(ii)<alfapalfap=-l(ii)/dl(ii);endif du(ii)<0&&-u(ii)/du(ii)<alfapalfap=-u(ii)/du(ii);endif dz(ii)<0&&-z(ii)/dz(ii)<alfadalfad=-z(ii)/dz(ii);endif dw(ii)>0&&-w(ii)/dw(ii)<alfadalfad=-w(ii)/dw(ii);endendalfap=0.9995*alfap;alfad=0.9995*alfad;x=x+alfap*dx;l=l+alfap*dl;u=u+alfap*du;y=y+alfad*dy;z=z+alfad*dz;w=w+alfad*dw;%迭代功率、电压幅值和相角for ii=1:ngPg(Nbus(ii))=x(ii);Qg(Nbus(ii))=x(ng+ii);endfor ii=1:nUamp(ii)=x(2*(ng+ii));Dlta(ii)=x(2*(ng+ii)-1);end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% k=k+1;elsebreak;endendfcost=0;for ii=1:ngfcost=fcost+a2(ii)*Pg(Nbus(ii))*Pg(Nbus(ii))+a1(ii)*Pg(Nbus(ii))+a0( ii);endfcostkplot(0:k,Gap(1:k+1),':*');PgQgUampfor ii=1:nif Type(ii)==3Dlta=Dlta-Dlta(ii)*ones(n,1);endendDlta。
最优潮流现代内点算法
n
j 1 n
i 1,..., n
P Gi PGi P Gi
Q
Ri
i SG
QRi Q Ri
2
i SR i SB
i j
i ST
2 Vi (ei2 f i2 ) V i
P ij Pij P ij
T i Ti T i
二、最优潮流计算方法现状
• 形成拉格朗日函数：
~u L f ( x) y T h( x) z T [ g ( x) l g ] wT [ g ( x) u g ] ~ zl w
三、现代内点Biblioteka 法• 导出KKT一阶最优性条件：
L x0 x f ( x) x h( x) y x g ( x)( z w) 0
五、仿真结果
算法性能
Ô Å ¶ ¼ ¼ ä Ï ¶ 1.E+3 1.E+2 1.E+1 1.E+0 1.E-1 1.E-2 1.E-3 1.E-4 1.E-5 1.E-6 1 6 11 16 21 ü ´ µ ú ´ Î Ê ý
Â Õ ú ó ¡ Ê Ï è
六、结论
a.
采用节点电压直角坐标模型，使其Hessian矩阵元素为常数， 不需每次迭代形成，方便编程的同时，加快计算速度； 新颖的数据结构定义了一个 4 4 的块矩阵，使待分解的系数 矩阵某部分具有与节点导纳矩阵相同的稀疏结构，方便使用稀 疏编程技巧，减少算法的计算时间； 新颖的数据结构减少注入元的产生，大大节约计算机内存，提 高算法的计算速度； 现代内点算法的超线性收敛性保证了算法的速度，其多项式时 间特性使算法具有良好的鲁棒性，更适合于大型电力系统的应 用。
matlab内点法
matlab内点法
Matlab内点法是一种数值优化方法，用于求解约束优化问题。
它通过将原始优化问题转化为一个无约束问题，并引入一系列罚函数或惩罚项来逼近原始问题的约束条件。
内点法的基本思想是从可行域内部（即“内点”）开始迭代，逐步逼近可行域的边界。
在Matlab中，可以使用内置的优化工具箱函数或编写自定义代码来实现内点法。
内点法的主要步骤包括初始化、构造目标函数和约束条件、选择合适的初始点、迭代更新内点、计算搜索方向和步长等。
Matlab提供了一些常用的内点法函数，例如fmincon，可以用于求解约束优化问题。
使用fmincon时，需要定义目标函数、约束条件、初始点等，并指定所需的优化算法和停止准则。
同时，还可以通过设置选项来调整算法的参数以获取更好的性能。
总之，Matlab内点法提供了一种有效的数值优化方法，可用于求解约束优化问题。
通过合理选择参数和算法，结合问题的特点，可以得到满足约束条件的最优解。
Matlab符号工具箱在实时电价中的应用
０ 引言
近年来 ， 随着 电力 市场 的深 入研究 ， 电力 系统 实时 电价 也越来越受 到电力工 作者 和学 者 的关 注 。电力 系统最 优潮 流 （ Ｐ ） 一 个 极 为 复 杂 的非 线 性 规 划 问 题 ， 满 足 特 定 的 ＯＦ 是 在
， ） 发 电机 有 功 生 产 成 本 （ 用 二 次 曲线 ） 无 功 机 会 成 本 （ 是 采 、
线计算 ， 编写程 序 简单 ， 用来校 验 一种 算 法和其他 程 序 的有 效性 。 但 可 关键 词 ： 实时 电价 ； 最优 潮 流 ； 原对偶 内点算 法 ； 号计 算 ； ｔｂ 符 Ｍａａ ｌ
中 图分类 号 ：Ｍ７ ４ Ｔ ４ 文 献标 识码 ： Ｂ 文章 编 号 ：６ ３— ０ ７ ２ ０ ） ２— ０ ｌ ３ １７ ０ ９ （ ０ ８ ０ ０ ３ 一０
参见文献 ； ） （ 为等式约束 向量 ， 这里为潮流平衡 方程 （ ｍ 维 ） ｇ ） ； （ 为不等式约束 向量 ， 这里 主要有 容量约束 、 电压约 束、 线路潮流约束 等（ 维 ） ｒ 。 １２ 原 对 偶 内点 算 法 ｎ 。 ． ］ 弓入松 弛变量将不等式约束化 为等式约束 ， 弓入 障碍 Ｉ 并 Ｉ 函数项 ， 根据实时 电价模型 可以定义 Ｌ ｍ ｎ函数 ， ａ ａ 如下 ： ， ， ｚ Ｗ ＝ （ （ Ａ，， ， ） ， ）一Ａ ｈ ｘ ， （ ）＋Ｚ ［ （ Ｔ ｇ ）一ｌ ］＋ —ｇ
和 Ｑ 的偏 导数 ， 即等 于最优 解处 各节 点 的拉格 朗 日乘子
的值 。
即：
＝
簧＝ ； ＝ Ａ ＝
＝
（ ５ ）
（） ６
（） ７ （） ８ （） ９ （ｏ １） （ １ １）