% function [loss22,Vdev,stc,V3,delta,fitness]=chaoliu(x)
% function [loss22,V3,delta]=chaoliu(x)
% clear all
% clc
%
% load node30x
% bus=node30x(1:30,:);
% branch=node30x(31:end,:);
% load IEEE30
% bus=IEEE30(1:30,:);
% branch=IEEE30(31:end,:);
load IEEE14
bus=IEEE14(1:14,:);
branch=IEEE14(15:end,:);
% load node14x
% bus=node14x(1:14,:);
% branch=node14x(15:end,:);
%对节点进行排序,按照先PQ,再PV,再平衡节点的顺序
[ nn,pp]=size(bus);
stc=0;
% % x=zeros(1,14);x(2)=200;
% % for i=1:nn
% % bus(i,4)=bus(i,4)+x(i);
% % stc=stc+x(i);
% % end
% for i=1:nn
% if x(i)~=0
% bus(i,4)=bus(i,4)+x(nn+i);
% stc=stc+abs(x(nn+i));
% end
% end
%nq为PQ节点数
nq=0;
for i=1:nn
if bus(i,10)==1
nq=nq+1;
end
end
%二进制无功补偿
% %对PQ节点进行无功补偿
% t=1;
% stc=0;
% for i=1:nn
% if bus(i,10)==1
% bus(i,4)=bus(i,4)+x(nq+t)*x(t);
% stc=stc+abs(x(nq+t))*x(t);
% t=t+1;
% end
% end
%十进制无功补偿
%对PQ节点进行无功补偿,对应nstcpso程序
% x=[5 0.5];
%
% [ro,cl]=size(x);
% for i=1:nn
% for j=1:cl/2
% if x(j)==i
% bus(i,4)=bus(i,4)+x(cl/2+j);
% end
% end
% end
% n_num=[3 4 6 7 9 10 12 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30];%IEEE30节点系统PQ节点编号
%%1-PQ节点,2-PV节点,3-平衡节点
bus1=[];
k=1;
for j=1:nn
if bus(j,10)==1 %先找出PQ节点排在前列
bus1(k,:)=bus(j,:);
k=k+1;
end
end
for j=1:nn
if bus(j,10)==2 %找出PV节点排在中间
bus1(k,:)=bus(j,:);
k=k+1;
end
end
for j=1:nn
if bus(j,10)==3 %平衡节点排在末尾
bus1(k,:)=bus(j,:);
k=k+1;
end
end
bus1;
[b,pp]=size(branch);
bran=[1:nn]';
bran=[bran,bus1];
branch1=branch;
for i=1:nn
vv1=find(branch1(:,2)==bran(i,2));
branch(vv1,2)=bran(i,1);
vv2=find(branch1(:,3)==bran(i,2));
branch(vv2,3)=bran(i,1);
end
%branch=input('branch=');
%将STATCOM的注入功率加入QG中,以便计算潮流
k=1;
%根据支路信息形成支路的关联矢量
M=zeros(nn,b);
for i=1:b
%判断为普通支路时,对From bus、To bus行分别赋值1、-1
if(branch(i,7)==1)
M(branch(i,2),i)=1;
M(branch(i,3),i)=-1;
%判断为变压器支路时,对From bus、To bus行分别赋值1、-1/t
else if(branch(i,7)==2)
M(branch(i,2),i)=1;
M(branch(i,3),i)=-1/branch(i,8);
k=k+1;
end
end
end
%由支路导纳及支路关联矢量生成节点导纳矩阵
j=sqrt(-1);
Y=zeros(nn,nn);
for i=1:b
Y=Y+M(:,i)*(1/(branch(i,4)+branch(i,5)*j))*M(:,i)';
end
%输出导纳矩阵
Y;
G=real(Y);B=imag(Y);
Pi1=find(bus1(:,10)~=3);%除平衡节点外的节点编号
PQi1=find(bus1(:,10)==1);%PQ节点编号
n=length(Pi1);
m=length(PQi1);
PS=(bus1(1:n,2)-bus1(1:n,4));
QS=(bus1(1:m,3)-bus1(1:m,5));
V=(bus1(:,6))';
delta=(bus1(:,7))';
for i=1:n
for j=1:n+1
Pn(j)=V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
Qn(j)=V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
end
P(i)=sum(Pn(:));
if i<=m
Q(i)=sum(Qn(:));
end
end
QS;
Q;
dQS=QS-Q';dPS=PS-P';
dPQ=[dPS;dQS];
dPQmax=max(abs(dPQ));
for k1=1:10
if dPQmax<=0.0001
break
end
n1=length(V);
H=zeros(n1);
for i=1:n+1
for j=1:n+1
if i~=j
H(i,j)=-V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
else
t=setdiff(1:(n+1),i);
H(i,i)=V(i)*sum(V(t).*(G(i,t).*sin(delta(repmat(i,1,n))-delta(t))-B(i,t).*cos(delta(repmat(i,1,n))-delta(t))));
end
end
end
H=H(1:n,1:n);
N=zeros(n1);
for i=1:n+1
for j=1:m
if i~=j
N(i,j)=-V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
else
t=setdiff(1:n+1,i);
N(i,i)=-V(i)*sum(V(t).*(G(i,t).*cos(delta(repmat(i,1,n))-delta(t))+B(i,t).*sin(delta(repmat(i,1,n))-delta(t))))-2*V(i)^2*G(i,i);
end
end
end
N=N(1:n,1:m);
M=zeros(n1);
for i=1:m
for j=1:n+1
if i~=j
M(i,j)=V(i)*V(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
else
t=setdiff(1:n+1,i);
M(i,i)=-V(i)*sum(V(t).*(G(i,t).*cos(delta(repmat(i,1,n))-delta(t))+B(i,t).*sin(delta(repmat(i,1,n))-delta(t))));
end
end
end
M=M(1:m,1:n);
L=zeros(m);
for i=1:m
for j=1:m
if i~=j
L(i,j)=-V(i)*V(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
else
t=setdiff(1:n+1,i);
L(i,i)=-V(i)*sum(V(t).*(G(i,t).*sin(delta(repmat(i,1,n))-delta(t))-B(i,t).*cos(delta(repmat(i,1,n))-delta(t))))+2*V(i)^2*B(i,i);
end
end
end
Jacobi=[H N;M L];
DeltaV=-Jacobi\dPQ;
delta(Pi1) =delta(Pi1)'+DeltaV(1:length(Pi1));
V(PQi1)=V(PQi1)+V(PQi1).*DeltaV(n+1:end)';
V1=V;
V2=V.*exp(delta*1j);
for i=1:n
for j=1:n+1
Pn(j)=V1(i)*V1(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
Qn(j)=V1(i)*V1(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
end
P(i)=sum(Pn(:));
if i<=m
Q(i)=sum(Qn(:));
end
end
dQS=QS-Q';dPS=PS-P';
dPQ=[dPS;dQS];
dPQmax=max(abs(dPQ));
end
k1;%迭代次数
S=zeros(length(V2),1);
for k=1:length(V2)
S(k)=V2(k)'*sum(conj(Y(k,:)).*conj(V2));
end
for k=1:b
m=branch(k,2);n=branch(k,3);%frombus,tobus
z=branch(k,4)+branch(k,5)*sqrt(-1);%R+jX
if branch(k,7)==2%变压器支路
loss1(k)=V2(m)*V2(m)*conj((branch(k,8)-1)/(branch(k,8)*z))+conj(V2(m))*(conj(V2(m))-conj(V2(n)))*conj(1/(z*branch(k,8)));
loss2(k)=V2(n)*V2(n)*conj((1-branch(k,8))/(branch(k,8)^2*z))+conj(V2(n))*(conj(V2(n))-conj(V2(m)))*conj(1/(z*branch(k,8)));
else
loss1(k)=V2(m)*V2(m)*conj(branch(k,6)*sqrt(-1)/2)+conj(V2(m))*(conj(V2(m))-conj(V2(n)))*conj(1/(z));
loss2(k)=V2(n)*V2(n)*conj(branch(k,6)*sqrt(-1)/2) +conj(V2(n))*(conj(V2(n))-conj(V2(m)))*conj(1/(z));
end
end
loss=(loss1+loss2)';%功率损耗
for i=1:nn %将节点编号还原,相应的电压、相角、功率编号也还原
for j=i:nn
if bus1(i,1)>bus1(j,1)
pp=V2(i);qq=delta(i);rr=S(i);ss=bus1(i,:);
V2(i)=V2(j);delta(i)=delta(j);S(i)=S(j);bus1(i,:)=bus1(j,:);
V2(j)=pp;delta(j)=qq;S(j)=rr;bus1(j,:)=ss;
end
end
end
Vdev=0;
for i=1:nn
Vdev=Vdev+(abs(V2(i))-1)^2;
end
Vdev=sqrt(Vdev);%总的电压偏差
% loss=abs(sum(loss));%总的功率损耗
stc;%接入STATCOM的总容量,标幺值
loss22=sum(abs(real(loss)));
fitness=1/3*(loss22+Vdev+stc);
% disp('有功损耗大小')
% disp(loss22)
% loss3=sum(abs(loss));
% disp('总的电压偏差')
% disp(Vdev)
% aa=[abs(V2)',delta'];
% loss4=[loss1',loss2',loss,abs(real(loss)),abs(imag(loss))];
% disp('电压大小')
% disp(abs(V2)')
% disp('电压相角(弧度)')
% disp(angle(V2)')
% AA=abs(V2');BB=angle(V2');
V3=abs(V2);
% flag=0;
% for i=1:length(n_num)
% if V3(n_num(i))<0.95|V3(n_num(i))>1.05
% flag=1;
% end
% end
% flag