3机9节点潮流计算编程

matlab潮流计算编程。可得出正确的可运行的雅可比矩阵和节点导纳矩阵。

clear;clc;
a=[%首节点t1	末节点t2	电阻r     电抗x   容纳之半ban  变比k
1		    	2		0.010		0.085		0.088	1;
1				3		0.017		0.092		0.079	1;
2				4		0.032		0.161		0.153	1;
3				6		0.039		0.170		0.179	1;
4				5		0.0085		0.072		0.0745  1;
5		    	6		0.0119		0.1008		0.1045	1;
9				1		0			0.0576		0		1;
7				4		0			0.0625		0		1;
8				6		0			0.0586		0		1;
];
b=[%节点	电压幅值 电压相角	 节点有功	节点无功	节点类型 /“1”为PQ节点，“2”为PV节点，“3”为平衡节点
1		1.0	 		0		0			0		1;
2		1.0	 		0		-1.2500		-0.5	1;
3		1.0	 		0		-0.9		-0.3	1;
4		1.0	 		0		0			0		1;
5		1.0	 		0		-1			-0.35	1;
6		1.0      	0		0			0		1;
7		1.025	 	0		1.6300		0		2;
8		1.025	 	0		0.8500		0		2;
9		1.040	 	0		0			0		3;
];

% a=load（‘zhilu.txt‘）;
% b=load（‘jiedian.txt‘）;
N=size（a1）;%支路数
Nbus=size（b1）;%节点数
pq=0;
for k=1:Nbus %PQ节点个数
    if b（k6）==1 
        pq=pq+1;
    end
end
Y=zeros（Nbus）;%节点导纳矩阵（生成一个Nbus*Nbus的零矩阵）
for k=1:N
    t1=a（k1）;t2=a（k2）;r=a（k3）;x=a（k4）;ban=a（k5）;K=a（k6）;
    Y（t1t1）=Y（t1t1）+1/（r+j*x）+j*ban; %电力工程P175
    Y（t1t2）=Y（t1t2）-1/（K*（r+j*x））;
    Y（t2t1）=Y（t2t1）-1/（K*（r+j*x））;
    Y（t2t2）=Y（t2t2）+1/（K*K*（r+j*x））+j*ban;
end
G=real（Y）;B=imag（Y）;      %提取导纳矩阵实部虚部
 
precision=1;           %误差精度
t=0 ;%存储迭代次数
while precision>0.00001   %只要精度一直大于0.00001就循环
    P=zeros（Nbus1）;%存储所有节点的有功
    Q=zeros（Nbus1）;%存储所有节点的无功
    for m=1:Nbus %求pv、pq和平衡节点的有功
        for n=1:Nbus
            P（m1）=P（m1）+b（m2）*b（n2）*（G（mn）*cos（b（m3）-b（n3））+B（mn）*sin（b（m3）-b（n3）））;
        %Pm=Um*E（1:n）Un[Gmn*cos（&m-&n）+Bmn*sin（&m-&n）]
        end
    end
    for m=1:Nbus %求pq、pv和平衡节点的无功
        for n=1:Nbus
            Q（m1）=Q（m1）+b（m2）*b（n2）*（G（mn）*sin（b（m3）-b（n3））-B（mn）*cos（b（m3）-b（n3）））;
       end
    end
    deltp=b（1:Nbus-14）-P（1:Nbus-11）;%pq和pv节点有功差（除去第9行平衡节点）
    deltq=b（1:pq5）-Q（1:pq1）;%pq节点无功差（只有pq有Q）
    deltPQ=[deltp;deltq]; %功率不平衡量
    H=zeros（Nbus-1）;%极坐标形式不用计算平衡节点
    for m=1:Nbus-1 %求H矩阵（雅可比HNJL）
        for n=1:Nbus-1
            if m~=n  %不等于，非对角元素，对&求偏导
                H（mn）=-b（m2）*b（n2）*（G（mn）*sin（b（m3）-b（n3））-B（mn）*cos（b（m3）-b（n3）））;
            else
                H（mm）=b（m2）*b（m2）*B（mm）+Q（m1）;%陈珩书P84
            end
        end
    end
    N=zeros（Nbus-1pq）;%求N矩阵
    for m=1:Nbus-1 
        for n=1:pq
            if m~=n
                N（mn）=-b（m2）*b（n2）*（G（mn）*cos（b（m3）-b（n3））+B（mn）*sin（b（m3）-b（n3）））;
            else
                N（mm）=-b（m2）*b（m2）*G（mm）-P（m1）;
            end
        end
    end
    J=zeros（pqNbus-1）;
    for m=1:pq %求J矩阵
        for n=1:Nbus-1
            if m~=n
                J（mn）=b（m2）*b（n2）*（G（mn）*cos（b（m3）-b（n3））+B（mn）*sin（b（m3）-b（n3）））;
            else
                J（mm）=b（m2）*b（m2）*G（mm）-P（m1）;
            end
        end
    end
    L=zeros（pqpq）;
    for m=1:pq %求L矩阵
        for n=1:pq
            if m~=n
                L（mn）=-b（m2）*b（n2）*（G（mn）*sin（b（m3）-b（n3））-B（mn）*cos（b（m3）-b（n3）））;
            else
 

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       3731  2018-08-06 16:34  3机9节点\motor3bus9.m

     文件     418408  2018-08-06 16:48  3机9节点\牛拉法潮流计算程序（附3机9节点结果对比）.docx

     目录          0  2018-08-06 16:48  3机9节点

----------- ---------  ---------- -----  ----

               422139                    3