Matlab牛顿-拉夫森法计算潮流的方法

许多朋友报告该程序报告了错误:

???在调用过程中未分配输出参数“总线”（可能还有其他参数）

“ E: \ MAT \最佳\ OpDF_.m> OpDF _”.

由于输入数据的格式，输入数据包含总线和线路，而不仅仅是两个矩阵，例如:

bus = [

11.00 0.00-0.30 -0.181;

21.00 0.00-0.55 -0.131;

31.100.000.500.00 2;

41.050.000.000.00 3];

line = [

12 0.10 0.400.00.015280;

42 0.08 0.400.00.014130;

14 0.12 0.500.00.01920;

31 0.000.3 0.00.0-1.1];

所有蓝色文本是数据文件的内容

---------------------分钟---------------------剪切--- ------------------ LINE ---------------------

【注意】

首先，该程序仅适用于求解以节点电压的极坐标表示的潮流方程，而无需考虑节点优化数

第二，该程序通过了Matlab 6.5的测试，应该适用于其所有后续版本

三，程序主要使用文件的输入输出，对输入文件的格式有严格的要求，如下:

1. 输入文件可以通过创建新的m文件直接在Matlab中编写，也可以将内容以文本形式编写，但是后缀必须更改为“ .m”文件. 文件名的第一个字符必须是字母，然后可以是字母，数字和下划线的任意组合，但不能与现有文件和功能冲突，并且不能包含中文.

2. 输入文件内容格式:

内容包括``总线''（节点数据）和``线路''（线路数据）两个矩阵

“公共汽车”的格式为:

bus = [节点编号节点电压节点相角（弧度系统）有功注入无功注入节点类型];

“行”的格式为:

line = [节点i节点j的线电阻电抗电导变压器比（公共线为零）];

第四，程序分为主程序和每个子程序，所有都需要放在Matlab的当前目录中进行调用，计算开始后，将在当前目录中生成文件“ Result.m”保存操作结果.

---------------------分钟---------------------剪切--- ------------------ LINE ---------------------

以下是程序代码（红色部分，在“％”之后的注释）:

主程序“ PowerFlow_NR.m”

函数[bus_res，S_res] = PowerFlow_NR_2％牛顿-拉夫森法求解潮流方程的主程序

%%%%%%%%%%%%%%%%%%%%%% bylongdinhohe %%%%%%%%%%%%%%%%%%%%%%%%% %% %%%%%%%%%%%%%%%%%%%

[总线，行] = OpDF _;％打开数据文件的子例程，将总线（节点数据）和行（行数据）返回到主程序

[nb，mb] =大小（总线）；

[nl，ml] =大小（行）；％计算总线和行矩阵的行数和列数

[总线，线路，nPQ，nPV，节点编号] = Num_（总线，线路）；％子例程以对节点重新排序

Y = y_（总线，线路）；％子例程，用于计算节点导纳矩阵

myf = fopen（'Result.m'，'w'）;

fprintf（myf，'---------------由longdinhohe ------------------------- ----- \ n \ n'）;

fclose（myf）;％在当前目录中生成一个“ Result.m”文件并写入节点导纳矩阵

长格式

EPS = 1.0e-10;设置误差的准确度

对于t = 1: 100％开始迭代计算，将最大迭代次数设置为100，以便及时跳出而不会收敛

[dP，dQ] = dPQ_（Y，总线，nPQ，nPV）；％子程序，用于计算功率偏差dP和dQ

J = Jac_（bus，Y，nPQ）;用于计算雅可比矩阵的％子程序

UD =零（nPQ，nPQ）;

对于i = 1: nPQ

UD（i，i）=总线（i，2）;％生成对角电压矩阵

结束

dAngU = J \ [dP; dQ];

dAng = dAngU（1: nb-1,1）;％计算相位角校正量

dU = UD *（dAngU（nb: nb + nPQ-1,1））;％计算电压校正量

总线（1: nPQ，2）=总线（1: nPQ，2）-dU;％校正电压

母线（1: nb-1,3）=母线（1: nb-1,3）-dAng;校正相角百分比

如果（max（abs（dU））

休息

end％判断精度误差是否满足，如果满足，将跳出，否则返回继续迭代

结束

bus = PQ_（总线，Y，nPQ，nPV）；％子例程，用于计算每个节点的有功和无功注入

[总线，线路] = ReNum_（总线，线路，nodenum）；％子例程，用于还原节点号

YtYm = YtYm_（线）；％子程序，用于计算线的等效Yt和Ym来计算线功率潮流

bus_res = bus_res_（bus）;用于计算节点数据结果的％子例程

S_res = S_res_（总线，线路，YtYm）；％子例程来计算线路潮流

myf = fopen（'Result.m'，'a'）;

fprintf（myf，'---------------牛顿-拉夫森法潮流计算结果---------- \ n \ nNode计算结果: \ n节点节点电压节点相角（角度）节点注入功率\ n'）;

对于i = 1: nb

fprintf（myf，'％2.0f'，bus_res（i，1））;

fprintf（myf，'. 6f'，bus_res（i，2））;

fprintf（myf，'.6f'，bus_res（i，3））;

fprintf（myf，'. 6f + j.6f \ n'，实数（bus_res（i，4）），imag（bus_res（i，4）））;

结束

fprintf（myf，'\ n线路计算结果: \ nNode I节点J线路功率S（I，J）线路功率S（J，I）线损dS（I，J）\ n'）; <

对于i = 1: nl

fprintf（myf，'％2.0f'，S_res（i，1））;

fprintf（myf，'％2.0f'，S_res（i，2））;

fprintf（myf，'. 6f + j.6f'，实数（S_res（i，3）），imag（S_res（i，3）））;

fprintf（myf，'. 6f + j.6f'，实数（S_res（i，4）），imag（S_res（i，4）））;

fprintf（myf，'. 6f + j.6f \ n'，实数（S_res（i，5）），imag（S_res（i，5）））;

结束

fclose（myf）;％迭代结束后继续将节点计算结果和行计算结果写入“ Result.m”中

子例程1“ OpDF_.m”用于打开数据文件

功能[总线，线路] = OpDF _

[dfile，pathname] = uigetfile（'*. m'，'Select DataFile'）;％数据文件类型为m文件，窗口标题为“ Select Data File”

如果路径名== 0

错误（“您必须选择一个有效的数据文件”）％如果未选择任何有效的文件，则会出现错误消息

其他

lfile =长度（dfile）;％读取文件字符串长度

eval_r（dfile（1: 1file-2））;％删除后缀并打开文件！注意！新浪博客中的“ eval”功能会自动添加“ _r”后缀，此处的功能名称为“ eval”而不是“ eval_r”，复制后请删除“ _r”后缀

子例程2“ Num.m”用于对节点重新排序并修改相应的行数据

功能[总线，线路，nPQ，nPV，节点编号] = Num_（总线，线路）

[nb，mb] =大小（总线）；

[nl，ml] =尺寸（线）；

nSW = 0；％nSW是平衡节点数

nPV = 0；％nPV是PV节点数

nPQ = 0；％nPQ是PQ节点数

对于i = 1: nb，％nb是节点总数

type =总线（i，6）;

如果类型== 3，

nSW = nSW +1;

SW（nSW，:) =总线（i，: ）;％计算并存储余额节点

elseif类型== 2，

nPV = nPV +1;

PV（nPV，:) =总线（i，: ）;％计算并存储PV节点

其他

nPQ = nPQ +1;

PQ（nPQ，:) =总线（i，: ）;％计算并存储PQ节点

结束

结束

bus = [PQ; PV; SW];％按照PQ，PV和平衡节点的顺序对总线矩阵进行重新排序

nodenum = [[1: nb]'bus（: ，1）];％生成新旧节点比较表

fori = 1: nl

forj = 1: 2

叉子= 1: nb

如果第（i，j）行==节点编号（k，2）

line（i，j）= nodenum（k，1）;

休息

结束

结束

结束

end％排序后按节点顺序对行矩阵重新编号

子程序3“ y_.m”充当节点导纳矩阵

函数Y = y_（总线，线路）

[nb，mb] =大小（总线）；

[nl，ml] =尺寸（线）；

Y =零（nb，nb）;％为导纳矩阵分配初始值0

叉子= 1: nl

I =行（k，1）;

J =行（k，2）;

Zt =第（k，3）行+ j *第（k，4）行;％读取行参数

ifZt〜= 0

Yt = 1 / Zt；未计算接地分支的％Yt

结束

Ym =行（k，5）+ j *行（k，6）;

K =行（k，7）;

如果（K == 0）＆（J〜= 0）％普通行: K = 0

Y（I，I）= Y（I，I）+ Yt + Ym;

Y（J，J）= Y（J，J）+ Yt + Ym;

Y（I，J）= Y（I，J）-Yt;

Y（J，I）= Y（I，J）;

结束

如果（K == 0）＆（J == 0）％分支到地: K = 0，J = 0，R = X = 0

Y（I，I）= Y（I，I）+ Ym;

结束

ifK> 0％变压器线路: Zt和Ym转换为i侧值，K在j侧

Y（I，I）= Y（I，I）+ Yt + Ym;

Y（J，J）= Y（J，J）+ Yt / K / K;

Y（I，J）= Y（I，J）-Yt / K;

Y（J，I）= Y（I，J）;

结束

ifK <0％变压器行: Zt和Ym是转换为K侧的值，K在i侧

Y（I，I）= Y（I，I）+ Yt + Ym;

Y（J，J）= Y（J，J）+ K * K * Yt;

Y（I，J）= Y（I，J）+ K * Yt;

Y（J，I）= Y（I，J）;

结束

结束

子程序4“ dPQ_.m”用于计算功率偏差

函数[dP，dQ] = dPQ_（Y，总线，nPQ，nPV）％nPQ和nPV是相应节点的数量

n = nPQ + nPV +1;节点总数的百分比

dP =总线（1: n-1,4）;

dQ =总线（1: nPQ，5）;％为dP和dQ PV节点分配初始值不需要计算dQ Balance节点不参与计算

对于i = 1: n-1

对于j = 1: n

dP（i，1）= dP（i，1）-总线（i，2）*总线（j，2）*（实数（Y（i，j）））* cos（总线（i，3） -bus（j，3））+ imag（Y（i，j））* sin（公共汽车（i，3）-bus（j，3）））;

如果我

dQ（i，1）= dQ（i，1）-总线（i，2）*总线（j，2）*（实数（Y（i，j）））* sin（总线（i，3） -bus（j，3））-imag（Y（i，j））* cos（公共汽车（i，3）-bus（j，3）））;

结束

结束

end％使用循环计算来找到dP和dQ

子程序5“ Jac_.m”用于计算雅可比矩阵

函数J = Jac_（总线，Y，nPQ）

[nb，mb] =大小（总线）；

H =零（nb-1，nb-1）;

N =零（nb-1，nPQ）;

K =零（nPQ，nb-1）;

L =零（nPQ，nPQ）；％除雅可比矩阵，即: J = [H N; K L]，然后初始化

Qi =零（nb-1,1）;

Pi =零（nb-1,1）;

对于i = 1: nb-1

对于j = 1: nb

Pi（i，1）= Pi（i，1）+总线（i，2）*总线（j，2）*（实数（Y（i，j））* cos（总线（i，3） -bus（j，3））+ imag（Y（i，j））* sin（公共汽车（i，3）-bus（j，3）））;

Qi（i，1）= Qi（i，1）+总线（i，2）*总线（j，2）*（真实（Y（i，j）））* sin（总线（i，3） -bus（j，3））-imag（Y（i，j））* cos（公共汽车（i，3）-bus（j，3）））;

结束

end％初始化并计算Qi和Pi

对于i = 1: nb-1

对于j = 1: nb-1

如果i〜= j

H（i，j）=-总线（i，2）*总线（j，2）*（实数（Y（i，j））* sin（总线（i，3）-总线（j，3 ））-imag（Y（i，j））* cos（公共汽车（i，3）-公共汽车（j，3）））;

其他

H（i，j）= Qi（i，1）+ imag（Y（i，j））*（（总线（i，2））^ 2）;

end％分别计算H矩阵的对角线元素和非对角线元素

如果j

如果i〜= j

N（i，j）=-总线（i，2）*总线（j，2）*（实数（Y（i，j））* cos（总线（i，3）-总线（j，3 ））+ imag（Y（i，j））* sin（公共汽车（i，3）-公共汽车（j，3）））;

其他

N（i，j）= -Pi（i，1）-real（Y（i，j））*（（总线（i，2））^ 2）;

结束

end％计算N矩阵的对角线元素和非对角线元素

如果我

如果i〜= j

K（i，j）=总线（i，2）*总线（j，2）*（实数（Y（i，j））* cos（总线（i，3）-总线（j，3） ）+ imag（Y（i，j））* sin（公共汽车（i，3）-公共汽车（j，3）））;

其他

K（i，j）= -Pi（i，1）+实数（Y（i，j））*（（总线（i，2））^ 2）;

end％计算K矩阵的对角线元素和非对角线元素

如果j

如果i〜= j

L（i，j）=-总线（i，2）*总线（j，2）*（实数（Y（i，j））* sin（总线（i，3）-总线（j，3 ））-imag（Y（i，j））* cos（公共汽车（i，3）-公共汽车（j，3）））;

其他

L（i，j）= -Qi（i，1）+ imag（Y（i，j））*（（总线（i，2））^ 2）;

结束

end％分别计算L矩阵的对角线元素和非对角线元素

结束

结束

结束

J = [H N; KL];％生成雅可比矩阵

子程序6“ PQ_.m”用于计算每个节点的功率注入

功能总线= PQ_（总线，Y，nPQ牛顿-拉夫逊法潮流计算，nPV）

n = nPQ + nPV + 1；％n是节点总数

对于i = nPQ + 1: n-1

对于j = 1: n

总线（i，5）=总线（i，5）+总线（i，2）*总线（j，2）*（实数（Y（i，j））* sin（总线（i，3） -bus（j，3））-imag（Y（i，j））* cos（公共汽车（i，3）-bus（j，3）））;

结束

end％使用公式来计算光伏节点的无功功率注入

对于j = 1: n

总线（n，4）=总线（n，4）+总线（n，2）*总线（j，2）*（实数（Y（n，j））* cos（总线（n，3） -bus（j，3））+ imag（Y（n，j））* sin（bus（n，3）-bus（j，3）））;

总线（n，5）=总线（n，5）+总线（n牛顿-拉夫逊法潮流计算，2）*总线（j，2）*（实数（Y（n，j））* sin（总线（n，3） -bus（j，3））-imag（Y（n，j））* cos（总线（n，3）-bus（j，3）））;

平衡节点的反应注入和主动注入的结束％计算

子例程7“ ReNum_.m”用于恢复节点和行数据的数量

功能[总线，线路] = ReNum_（总线，线路，节点编号）

[nb，mb] =大小（总线）；

[nl，ml] =尺寸（线）；

bus_temp =零（nb，mb）;％bus_temp矩阵用于临时存储总线矩阵数据

k = 1;

对于i = 1: nb

对于j = 1: nb

如果总线（j，1）== k

bus_temp（k，:) =总线（j，: ）;

k = k +1;

结束

结束

end％使用总线矩阵的第一列号对总线矩阵进行重新排序并将其存储在bus_temp矩阵中

bus = bus_temp;％重新分配回总线，完成总线矩阵的重新编号

fori = 1: nl

forj = 1: 2

叉子= 1: nb

如果第（i，j）行==节点编号（k，1）

line（i，j）= nodenum（k，2）;

休息

结束

结束

结束

结束％恢复行号

子程序8“ YtYm_.m”用于计算线路的等效Yt和Ym，以计算线路功率潮流

函数YtYm = YtYm_（行）

[nl，ml] =尺寸（线）；

YtYm =零（nl，5）;％为YtYm矩阵分配初始值0

YtYm（: ，1: 2）= line（: ，1: 2）;％矩阵的前两列是该行的两个节点号，后三列是与行等效的Yt， i侧的等效Ym，j侧的等效Ym

j = sqrt（-1）;

叉子= 1: nl

I =行（k，1）;

J =行（k，2）;

Zt =行（k，3）+ j *行（k，4）;

ifZt〜= 0

Yt = 1 / Zt;

结束

Ym =行（k，5）+ j *行（k，6）;

K =行（k，7）;

如果（K == 0）＆（J〜= 0）％普通行: K = 0

YtYm（k，3）= Yt;

YtYm（k，4）= Ym;

YtYm（k，5）= Ym;

结束

如果（K == 0）＆（J == 0）％分支到地: K = 0，J = 0，R = X = 0

YtYm（k，4）= Ym;

结束

ifK> 0％变压器线路: Zt和Ym转换为i侧值，K在j侧

YtYm（k，3）= Yt / K;

YtYm（k，4）= Ym + Yt *（K-1）/ K;

YtYm（k，5）= Yt *（1-K）/ K / K;

结束

ifK <0％变压器行: Zt和Ym是转换为K侧的值，K在i侧

YtYm（k，3）= -Yt * K;

YtYm（k，4）= Ym + Yt *（1 + K）;

YtYm（k，5）= Yt *（K ^ 2 + K）;

结束

结束

子程序9“ bus_res_.m”计算并返回节点数据结果

功能bus_res = bus_res_（总线）;

[nb，mb] =大小（总线）；

bus_res =零（nb，4）;％bus_res矩阵存储节点计算结果

bus_res（: ，1: 2）=总线（: ，1: 2）;

bus_res（: ，3）=总线（: ，3）* 180 / pi;％相角采用角度系统

bus_res（: ，4）=总线（: ，4）+（sqrt（-1））*总线（: ，5）;％注入功率

子程序10“ S_res_.m”计算并返回线路功率流

函数S_res = S_res_（总线，线路，YtYm）

[nl，ml] =尺寸（线）；

S_res =零（nl，5）;％S_res矩阵存储潮流计算结果

S_res（: ，1: 2）= line（: ，1: 2）;％前两列是节点编号

对于k = 1: nl

I = S_res（k，1）;

J = S_res（k，2）;

如果（J〜= 0）＆（I〜= 0）

S_res（k，3）=总线（I，2）^ 2 *（conj（YtYm（k，3））+ conj（YtYm（k，4）））-总线（I，2）*总线（ J，2）*（cos（总线（I，3））+ j * sin（总线（I，3）））*（conj（cos（总线（J，3））+ j * sin（总线（J， 3））））* conj（YtYm（k，3））;

S_res（k，4）=总线（J，2）^ 2 *（conj（YtYm（k，3））+ conj（YtYm（k，5）））-总线（I，2）*总线（ J，2）*（conj（cos（公共汽车（I，3））+ j * sin（公共汽车（I，3））））*（cos（公共汽车（J，3））+ j * sin（公共汽车（J ，3）））* conj（YtYm（k，3））;

S_res（k，5）= S_res（k，3）+ S_res（k，4）;％使用公式计算不接地分支的功率流

elseifJ == 0

S_res（k，3）=总线（I，2）^ 2 * conj（YtYm（k，4））;

S_res（k，5）= S_res（k，3）+ S_res（k，4）;

其他

S_res（k，4）=总线（J，2）^ 2 * conj（YtYm（k，5））;

S_res（k，5）= S_res（k，3）+ S_res（k，4）;％使用公式计算接地分支的功率流

结束

结束

本文来自电脑杂谈，转载请注明本文网址：
http://www.pc-fly.com/a/tongxinshuyu/article-164677-1.html