牛顿拉夫逊法潮流计算(直角坐标)

%节点标号无特殊要求

clc

if nargin < 1

% linedata = [

% 1 2 0.1 0.4 0.01528 1;

% 1 3 0 0.3 0 1.1;

% 1 4 0.12 0.5 0.01920 1;

% 2 4 0.08 0.40 0.01413 1];

% busdata = [

% 1 0 -0.30 -0.18 1 0;

% 2 0 -0.55 -0.13 1 0;

% 3 2 0.5 0 1.10 0;

% 4 1 0 0 1.05 0;];

linedata = [

1 4 0.1 0.4 0.01528 1;

1 3 0 0.3 0 1.1;

1 2 0.12 0.5 0.01920 1;

4 2 0.08 0.40 0.01413 1];

busdata = [

1 0 -0.30 -0.18 1 0;

4 0 -0.55 -0.13 1 0;

3 2 0.5 0 1.10 0;

2 1 0 0 1.05 0;];

end

if nargin < 3

error = 10^-5;

end

clc

%优化节点排序

pqbus = find(busdata(:2)==0);

pvbus = find(busdata(:2)==2);

baxxxxsebus = find(busdata(:2)==1);

bussort = [pqbus' pvbus' baxxxxsebus'];

Num = busdata(bussort1);

%对节点进行重新排列，使其排列顺序为 PQ节点->PV节点->平衡节点

linedata = NumResort(linedataNum);

P = busdata(:3);

Q = busdata(:4);

U = busdata(:5);

Us = U;

deta = busdata(:6);

e = U.*cos(deta);

f = U.*sin(deta);

busnum = length(busdata(bussort1));

pqnum = length(pqbus);

%生成节点导纳矩阵

Y = BuildY(linedata);

for k = 1:100

[dPdQdU2] = CalcCollection(PQUsefYbusnumpqnum);

if max(abs([dP dQ]))break;

end

J = CalcJ(efYbusnumpqnum)

dPQU = zeros(busnum-11);

for ii = 1:pqnum

dPQU(2*ii-1) = dP(ii);

dPQU(2*ii) = dQ(ii);

end

for ii = pqnum+1:busnum-1

dPQU(2*ii-1) = dP(ii);

dPQU(2*ii) = dU2(ii);

end

def = Jordan(JdPQU);

for ii = 1:busnum-1

f(ii) = f(ii)+def(2*ii-1);

e(ii) = e(ii)+def(2*ii);

end

end

if k >= 100

disp('不收敛!');

return;

end

[SbusSijSvc] = CalcS(linedataefYbusnum);

disp('编号:');

disp(Num');

disp('节点电压:');

disp(e'+1j*f');

disp('平衡节点功率:');

disp(Svc);

disp('线路功率');

for n = 1:length(Sij(:1))

fprintf('%o %o 'Num(Sij(n1))Num(Sij(n2)))

disp(Sij(n3))

end

disp('节点功率');

for n = 1:busnum

fprintf('%o 'Num(n))

disp(Sbus(n))

end

end

function [SbusSijSvc] = CalcS(linedataefYbusnum)

Sbus = zeros(1busnum);

Qi = zeros(1busnum);

Ud = e+1j*f;

%公式11-57

for ii = 1:busnum

for jj = 1:busnum

Sbus(ii) = Sbus(ii)+Ud(ii)*conj(Y(iijj))*conj(Ud(jj));

end

end

Svc = Sbus(busnum);

Sij = zeros(length(linedata(:1))3);

for n=1:length(linedata(:1))

ii=linedata(n1);

jj=linedata(n2);

Sij(2*n-11) = ii;

Sij(2*n-12) = jj;

Sij(2*n-13) = Ud(ii)*(conj(Ud(ii))...

*conj(1j*linedata(n5)+linedata(n6)*(linedata(n6)-1)/(linedata(n3)+1j*linedata(n4)))...

+(conj(Ud(ii))-conj(Ud(jj)))*conj(-Y(iijj)));

ii=linedata(n2);

jj=linedata(n1);

Sij(2*n1) = ii;

Sij(2*n2) = jj;

Sij(2*n3) = Ud(ii)*(conj(Ud(ii))...

*conj(1j*linedata(n5)+(1-linedata(n6))/(linedata(n3)+1j*linedata(n4)))...

+(conj(Ud(ii))-conj(Ud(jj)))*conj(-Y(iijj)));

end

end

function linedata1 = NumResort(linedataNum)

%对节点进行重排

linedata1 = linedata;

for n = 1:length(Num)

n1 = find(linedata(:1)==Num(n));

n2 = find(linedata(:2)==Num(n));

linedata1(n11) = n*ones(length(n1)1);

linedata1(n22) = n*ones(length(n2)1);

end

end

function J = CalcJ(efYbusnumpqnum)

%计算雅可比矩阵

G = real(Y);

B = imag(Y);

J = zeros(2*(busnum-1));

H = zeros(pqnumpqnum);

N = H;

M = H;

L = H;

for ii = 1:pqnum

for jj = 1:pqnum

if ii~=jj

H(iijj) = -B(iijj)*e(ii)+G(iijj)*f(ii);

N(iijj) = G(iijj)*e(ii)+B(iijj)*f(ii);

M(iijj) = -N(iijj);

L(iijj) = H(iijj);

J(2*ii-12*jj-1) = H(iijj);

J(2*ii-12*jj) = N(iijj);

J(2*ii2*jj-1) = M(iijj);

J(2*ii2*jj) = L(iijj);

end

end

Iii = G(iiii)*e(ii)-B(iiii)*f(ii)...

+1j*(G(iiii)*f(ii)+B(iiii)*e(ii));

for jj = 1:busnum

if ii~=jj

Iii = Iii+G(iijj)*e(jj)-B(iijj)*f(jj)...

+1j*(G(iijj)*f(jj)+B(iijj)*e(jj));

end

end

H(iiii) = -B(iiii)*e(ii)+G(iiii)*f(ii)+imag(Iii);

N(iiii) = G(iiii)*e(ii)+B(iiii)*f(ii)+real(Iii);

M(iiii) = -G(iiii)*e(ii)-B(iiii)*f(ii)+real(Iii);

L(iiii) = -B(iiii)*e(ii)+G(iiii)*f(ii)-imag(Iii);

J(2*ii-12*ii-1) = H(iiii);

J(2*ii-12*ii) = N(iiii);

J(2*ii2*ii-1) = M(iiii);

J(2*ii2*ii) = L(iiii);

end

H = zeros(pqnumbusnum-pqnum-2);

N = H;

M = H;

L = H;

for ii = 1:pqnum

for jj = pqnum+1:busnum-1

if ii~=jj

H(iijj) = -B(iijj)*e(ii)+G(iijj)*f(ii);

N(iijj) = G(iijj)*e(ii)+B(iijj)*f(ii);

M(iijj) = -N(iijj);

L(iijj) = H(iijj);

J(2*ii-12*jj-1) = H(iijj);

J(2*ii-12*jj) = N(iijj);

J(2*ii2*jj-1) = M(iijj);

J(2*ii2*jj) = L(iijj);

end

end

end

H = zeros(busnum-pqnum-2pqnum);

N = H;

R = H;

S = H;

for ii = pqnum+1:busnum-1

for jj = 1:pqnum

if ii~=jj

H(iijj) = -B(iijj)*e(ii)+G(iijj)*f(ii);

N(iijj) = G(iijj)*e(ii)+B(iijj)*f(ii);

R(iijj) = 0;

S(iijj) = 0;

J(2*ii-12*jj-1) = H(iijj);

J(2*ii-12*jj) = N(iijj);

J(2*ii2*jj-1) = R(iijj);

J(2*ii2*jj) = S(iijj);

end

end

end

H = zeros(busnum-pqnum-2);

N = H;

R = H;

S = H;

for ii = pqnum+1:busnum-1

for jj = pqnum+1:busnum-1

if ii~=jj

H(iijj) = -B(iijj)*e(ii)+G(iijj)*f(ii);

N(iijj) = G(iijj)*e(ii)+B(iijj)*f(ii);

R(iijj) = 0;

S(iijj) = 0;

J(2*ii-12*jj-1) = H(iijj);

J(2*ii-12*jj) = N(iijj);

J(2*ii2*jj-1) = R(iijj);

J(2*ii2*jj) = S(iijj);

end

end

Iii = G(iiii)*e(ii)-B(iiii)*f(ii)...

+1j*(G(iiii)*f(ii)+B(iiii)*e(ii));

for jj = 1:busnum

if ii~=jj

Iii = Iii+G(iijj)*e(jj)-B(iijj)*f(jj)...

+1j*(G(iijj)*f(jj)+B(iijj)*e(jj));

end

end

H(iiii) = -B(iiii)*e(ii)+G(iiii)*f(ii)+imag(Iii);

N(iiii) = G(iiii)*e(ii)+B(iiii)*f(ii)+real(Iii);

R(iiii) = 2*f(ii);

S(iiii) = 2*e(ii);

J(2*ii-12*ii-1) = H(iiii);

J(2*ii-12*ii) = N(iiii);

J(2*ii2*ii-1) = R(iiii);

J(2*ii2*ii) = S(iiii);

end

end

function [dPdQdU2] = CalcCollection(PQUsefYbusnumpqnum)

%计算修正值

G = real(Y);

B = imag(Y);

Pt = P;

Qt = Q;

dP = zeros(1busnum);

dQ = zeros(1busnum);

dU2 = zeros(1busnum);

for ii = 1:pqnum

Pi = 0;

Qi = 0;

for jj = 1:busnum

Pi = Pi+e(ii)*(G(iijj)*e(jj)-B(iijj)*f(jj))...

+f(ii)*(G(iijj)*f(jj)+B(iijj)*e(jj));

Qi = Qi+f(ii)*(G(iijj)*e(jj)-B(iijj)*f(jj))...

-e(ii)*(G(iijj)*f(jj)+B(iijj)*e(jj));

end

dP(ii) = P(ii)-Pi;

dQ(ii) = Q(ii)-Qi;

Pt(ii) = Pi;

Qt(ii) = Qi;

end

for ii = pqnum+1:busnum-1

Pi = 0;

for jj = 1:busnum

Pi = Pi+e(ii)*(G(iijj)*e(jj)-B(iijj)*f(jj))...

+f(ii)*(G(iijj)*f(jj)+B(iijj)*e(jj));

end

dP(ii) = P(ii)-Pi;

Pt(ii) = Pi;

dU2(ii) = Us(ii)^2-(e(ii)^2+f(ii)^2);

end

end

function Y = BuildY(linedata)

if nargin < 1

linedata = [

1 4 0.1 0.4 0.01528 1; 1 3 0 0.3 0 1.1; 1 2 0.12 0.5 0.01920 1; 4 2 0.08 0.40 0.01413 1];end

nf = linedata(:1);

nt = linedata(:2);r = linedata(:3);

x = linedata(:4);

b = linedata(:5);

k = linedata(:6);

branchnum = length(nf);

busnum = max([nf'nt']);

y = ones(branchnum1)./(r+1j*x);

Y = zeros(busnum);

for n = 1:branchnum

Y(nf(n)nt(n)) = Y(nf(n)nt(n))-k(n)*y(n);

Y(nt(n)nf(n)) = Y(nf(n)nt(n));

Y(nf(n)nf(n)) = Y(nf(n)nf(n))+k(n)^2*y(n)+1j*b(n);

Y(nt(n)nt(n)) = Y(nt(n)nt(n))+y(n)+1j*b(n);

end

end

function s = Jordan(Ab)

%约当消元

if nargin < 1

A=[1 2 3;2 7 5;1 4 9];

b=[1 6 -3]';

end

if rank(A)~=rank([Ab])

error('A矩阵的秩和增广矩阵的秩不相同方程不存在唯一解');

end

[~n]=size(A);

A(:n+1) = b;

for k = 1:n

A(kk:end) = A(kk:end)/A(kk);

% 规格化

for q = 1:k-1

% 消上三角

if A(qk)~=0

A(qk:end) = A(qk:end)-A(kk:end).*A(qk);

end

end

for p = k+1:n

% 消下三角

if A(pk)~=0

A(pk:end) = A(pk:end)-A(kk:end).*A(pk);

end

end

end

s = A(:end);

end