clc
clear
load('Feeder.mat')
NumN = Feeder.NumN;NumL = Feeder.NumL;
NODES = zeros(NumN,6);
LINES = zeros(NumL,3);
for k = 1:NumL
N1 = Feeder.Topology(k,1);
N2 = Feeder.Topology(k,2);
Rkm = real(sum(diag(Feeder.ZLIN(:,:,1)))/3)*Feeder.Topology(k,3);
LINES(k,:) = [N1 N2 Rkm];
end
rand('twister',12345)
Loadpropose1 = 0.8 + (0.4)*rand(NumN,3);
Loadpropose2 = randi([0,1],NumN,3);
LOAD = Loadpropose1.*Loadpropose2;
LOAD = 2*3490*LOAD/(sum(sum(LOAD)));
Imax = 1.1*LOAD/4.16;
Nodes_ID = unique(LINES(:,[1 2]));
NODES= [Nodes_ID LOAD Imax];
Solid = 0;
Vary = 1;
Vb = 4.16;
Pb = 1000;
Rb = (Vb^2)/(Pb/1000);
Ibase = Pb/Vb;
Vref = 1;
Damping = [20 30];
Nodog = [80 120];
tic
NE = 200;
LINES(:,3) = LINES(:,3)/Rb;
NODES(:,2:4) = NODES(:,2:4)/Pb;
NODES(:,5:end) = NODES(:,5:end)/Ibase;
for kk = 1:NE
if Vary == 1
NODESin(:,2:4) = NODES(:,2:4)*(1.25*rand+0.75);
else
NODESin(:,2:4) = NODES(:,2:4);
end
N = size(NODES,1);
L = size(LINES,1);
A = zeros(N,L);
for k = 1:L
n = LINES(k,1); m = LINES(k,2);
A(n,k) = 1; A(m,k) = -1;
end
A3p = zeros(3*N,3*L);
for m = 1:L
for k = 1:N
if A(k,m) == 1
A3p(3*k-2:3*k,3*m-2:3*m) = eye(3);
elseif A(k,m) == -1
A3p(3*k-2:3*k,3*m-2:3*m) = -eye(3);
end
end
end
G = diag(1./LINES(:,3));
G3p = zeros(3*L,3*L);
for k = 1:L
G3p(3*k-2:3*k,3*k-2:3*k) = G(k,k)*eye(3);
end
G3bus = A3p*G3p*A3p.';
G3bds = G3bus(4:end,1:3);
G3bdd = G3bus(4:end,4:end);
Z3bdd = inv(G3bdd);
Vdpn = ones(3*N-3,1);
Vdpn0 = ones(3*N-3,1);
PdpnY = ones(2*N-2,1);
PdpnD = ones(2*N-2,1);
Vspn = [1;0;-1];
for k = 1:N-1
Vdpn0(3*k-2:3*k,1) = Vspn;
if Vary == 0
Vdpn0(3*k-2:3*k,1) = (1.1*rand+0.2)*Vspn;
end
PdpnY(2*k-1:2*k,1) = [NODESin(k+1,2);NODESin(k+1,3)];
PdpnD(2*k-1:2*k,1) = [NODESin(k+1,4);NODESin(k+1,4)];
end
z = zeros(3*(N-1),1); k2 = 1;
for k = 2:N
for k1 = 1:3
if NODES(k,k1+1) == 0
else
z(k2) = NODES(k,k1+4)^2/abs(NODESin(k,k1+1));
end
k2 = k2 +1;
end
end
L = norm(z)*norm(Z3bdd);
Lkk(kk) = L;
itermax = 1000;
e = 1e-10;
IdpnY = zeros(size(Vdpn0));
IdpnD = zeros(size(Vdpn0));
for t = 1:itermax
ll = 0;
for k = 2:N
if any(k == Nodog)
ll = 1+ ll;
IdpnY((3*k-5),1) = (PdpnY(2*k-3)-Damping(ll)*(Vref-Vdpn0(3*k-5,1)))/(Vdpn0(3*k-5,1)-Vdpn0(3*k-4,1));
IdpnY((3*k-4),1) = PdpnY(2*k-3)/(Vdpn0(3*k-4,1)- Vdpn0(3*k-5,1)) + PdpnY(2*k-2)/(Vdpn0(3*k-4,1)- Vdpn0(3*k-3,1));
IdpnY((3*k-3),1) = (PdpnY(2*k-2)+Damping(ll)*(-Vref-Vdpn0(3*k-3,1)))/(Vdpn0(3*k-3,1)- Vdpn0(3*k-4,1));
else
IdpnY((3*k-5),1) = PdpnY(2*k-3)/(Vdpn0(3*k-5,1)-Vdpn0(3*k-4,1));
IdpnY((3*k-4),1) = PdpnY(2*k-3)/(Vdpn0(3*k-4,1)- Vdpn0(3*k-5,1)) + PdpnY(2*k-2)/(Vdpn0(3*k-4,1)- Vdpn0(3*k-3,1));
IdpnY((3*k-3),1) = PdpnY(2*k-2)/(Vdpn0(3*k-3,1)- Vdpn0(3*k-4,1));
end
IdpnD((3*k-5),1) = PdpnD(2*k-3)/(Vdpn0(3*k-5,1) - Vdpn0(3*k-3,1));
IdpnD((3*k-4),1) = 0;
IdpnD((3*k-3),1) = PdpnD(2*k-2)/(Vdpn0(3*k-3,1) - Vdpn0(3*k-5,1));
end
Idpn = IdpnY + IdpnD;
Vdpn = -Z3bdd*(G3bds*Vspn + Idpn);
if Solid == 0
for m = 1:N-1
Vdpn(3*m-1,1) = 0;
end
end
error(t,kk) = max(abs(abs(Vdpn)-abs(Vdpn0)));
if error(t,kk) < e
Vpn = [Vspn;Vdpn];
break;
else
Vdpn0 = Vdpn;
end
end
end
t1 = toc;
Vall = zeros(NumN,4);
for k = 2:NumN+1
Vall(k-1,:) = [k-1, Vpn(3*k-5), Vpn(3*k-4), Vpn(3*k-3)] ;
end