nbus=5; nlines=5; tolerance=0.00001; tap=[1;1;1;1.05;1.05]; x=[1;1;2;2;3]; y=[2;3;3;4;5]; zreal=[0.04;0.1;0.08;0;0]; zimag=[0.25;0.35;0.3;0.015;0.03]; y1=[0.25;0;0.25;0;0]; Ybus=zeros(nbus); for i=1:length(x) tap_effect11=(1/(zreal(i)+zimag(i)*j))/(tap(i)^2); tap_effect12=(-1/(zreal(i)+zimag(i)*j))/(tap(i)); tap_effect21=(-1/(zreal(i)+zimag(i)*j))/(tap(i)); tap_effect22=(1/(zreal(i)+zimag(i)*j)); Ybus(x(i),x(i))=Ybus(x(i),x(i))+tap_effect11+(y1(i))*j; Ybus(y(i),y(i))=Ybus(y(i),y(i))+tap_effect22+(y1(i))*j; Ybus(x(i),y(i))=Ybus(x(i),y(i))+tap_effect12; Ybus(y(i),x(i))=Ybus(y(i),x(i))+tap_effect21; end Ybus; G=real(Ybus); B=imag(Ybus); v_bus=[1.05;1.05;1.05;1.05;1.05]; ang_bus=[0;0;0;0;0]; p_real=[-1.6;-2;-3.7;5]; q_real=[-0.8;-1;-1.3]; count=0; flag=1; while flag>0 p_calc=[]; q_calc=[]; for i=1:nbus p_calc(i)=0; q_calc(i)=0; for j=1:nbus p_calc(i)=p_calc(i)+v_bus(i)*v_bus(j)*(G(i,j)*cos(ang_bus(i)-ang_bus(j))+B(i,j)*sin(ang_bus(i)- ang_bus(j))); q_calc(i)=q_calc(i)+v_bus(i)*v_bus(j)*(G(i,j)*sin(ang_bus(i)-ang_bus(j))-B(i,j)*cos(ang_bus(i)-a ng_bus(j))); end 1 

end p_calcu=p_calc; p_calcu(5)=[]; mmp=p_real'-p_calcu; q_calcu=q_calc; q_calcu(4:5)=[]; mmq=q_real'-q_calcu; mm=[mmp mmq]'; nstate=2*nbus; Hp=zeros(nbus,nstate); Hq=zeros(nbus,nstate); bus_number=[]; for i=1:nbus bus_number=[bus_number i]; end for i=1:length(mmp) for j=1:nstate if i==j Hp(i,j)=-q_calc(i)-v_bus(i)*v_bus(i)*B(i,i); end if j~=i & j<=nbus Hp(i,j)=v_bus(i)*v_bus(j)*(G(i,j)*sin(ang_bus(i)-ang_bus(j))-B(i,j)*cos(ang_bus(i)-ang_bus(j))); end if j==i+nbus Hp(i,j)=p_calc(i)/v_bus(i)+ v_bus(i)*G(i,i); end if j~=i+nbus & j>nbus Hp(i,j)=v_bus(i)*(G(i,j-nbus)*cos(ang_bus(i)-ang_bus(j-nbus))+B(i,j-nbus)*sin(ang_bus(i)-ang_b us(j-nbus))); end end end for i=1:length(mmq) for j=1:nstate if i==j Hq(i,j)=p_calc(i)-v_bus(i)*v_bus(i)*G(i,i); end if j~=i & j<=nbus 2 

Hq(i,j)=-v_bus(i)*v_bus(j)*(G(i,j)*cos(ang_bus(i)-ang_bus(j))+B(i,j)*sin(ang_bus(i)-ang_bus(j))) ; end if j==i+nbus Hq(i,j)=q_calc(i)/v_bus(i)-v_bus(i)*B(i,i); end if j~=i+nbus & j>nbus Hq(i,j)=v_bus(i)*(G(i,j-nbus)*sin(ang_bus(i)-ang_bus(j-nbus))-B(i,j-nbus)*cos(ang_bus(i)-ang_b us(j-nbus))); end end end Hp(5,:)=[]; Hq(4:5,:)=[]; H=[Hp; Hq]; H(:,9: 10)=[]; H(:,5)=[]; U=chol(H); smm=U\(U'\mm); if max(abs(smm))>tolerance flag=1; smma=smm(1:nbus-1,:); smmv=smm(nbus:2*nbus-3,:); smma=[smma;0]; smmv=[smmv;0;0]; count=count+1; v_bus=v_bus+smmv; ang_bus=ang_bus+smma; else flag=0; end end%对应 while% fprintf('\nConvergence has occurred in %d interations.\n', count+1); fprintf('\nEstimated bus voltage.'); fprintf('\nBus_number for i=1:nbus voltage\n'); fprintf('%f %f\n', bus_number(i), v_bus(i)); end fprintf('\nEstimated bus angle.'); 3 

fprintf('\nBus_number for i=1:nbus angle\n'); fprintf('%f %f\n', bus_number(i), ang_bus(i)); end p_calc=[]; q_calc=[]; for i=1:nbus p_calc(i)=0; q_calc(i)=0; for j=1:nbus p_calc(i)=p_calc(i)+v_bus(i)*v_bus(j)*(G(i,j)*cos(ang_bus(i)-ang_bus(j))+B(i,j)*sin(ang_bus(i)- ang_bus(j))); q_calc(i)=q_calc(i)+v_bus(i)*v_bus(j)*(G(i,j)*sin(ang_bus(i)-ang_bus(j))-B(i,j)*cos(ang_bus(i)-a ng_bus(j))); end end fprintf('\nEstimated real power bus injection.'); fprintf('\nBus_number P\n'); for i=1:nbus fprintf('%f %f\n', bus_number(i), p_calc(i)); end fprintf('\nEstimated reactive power bus injection.'); fprintf('\nBus_number Q\n'); for i=1:nbus fprintf('%f %f\n', bus_number(i), q_calc(i)); end start_bus=x; end_bus=y; pij_calc=[]; qij_calc=[]; for i=1:length(x) q=start_bus(i); r=end_bus(i); pij_calc(i)=-((v_bus(q))^2)*G(q,r)/tap(i)+v_bus(q)*v_bus(r)*(G(q,r)*cos(ang_bus(q)-ang_bus(r)) +B(q,r)*sin(ang_bus(q)-ang_bus(r))); qij_calc(i)=((v_bus(q))^2)*(B(q,r)/tap(i)-y1(i))+v_bus(q)*v_bus(r)*(G(q,r)*sin(ang_bus(q)-ang_b us(r))-B(q,r)*cos(ang_bus(q)-ang_bus(r))); end fprintf('\nReal line flows.'); fprintf('\nFrom_bus To_bus for i=1:length(start_bus) Pij\n'); 4 

fprintf('%f %f %f\n', start_bus(i), end_bus(i), pij_calc(i)); end fprintf('\nReactive line flows.'); fprintf('\nFrom_bus To_bus for i=1:length(start_bus) Qij\n'); fprintf('%f %f %f\n', start_bus(i), end_bus(i), qij_calc(i)); end pji_calc=[]; qji_calc=[]; for i=1:length(x) q=start_bus(i); r=end_bus(i); pji_calc(i)=-((v_bus(r))^2)*tap(i)*G(r,q)+v_bus(r)*v_bus(q)*(G(r,q)*cos(ang_bus(r)-ang_bus(q)) +B(r,q)*sin(ang_bus(r)-ang_bus(q))); qji_calc(i)=((v_bus(r))^2)*(tap(i)*B(r,q)-y1(i))+v_bus(r)*v_bus(q)*(G(r,q)*sin(ang_bus(r)-ang_b us(q))-B(r,q)*cos(ang_bus(r)-ang_bus(q))); end fprintf('\nReal line flows.'); fprintf('\nFrom_bus To_bus for i=1:length(start_bus) Pji\n'); fprintf('%f %f %f\n', end_bus(i), start_bus(i), pji_calc(i)); end fprintf('\nReactive line flows.'); fprintf('\nFrom_bus To_bus for i=1:length(start_bus) Qji\n'); fprintf('%f %f %f\n', end_bus(i), start_bus(i), qji_calc(i)); end 5