matlab数据边界不对,matlab索引超出数组边界问题

%本程序的功能是用牛顿——拉夫逊法进行潮流计算

clear;

n=6;

nl=6;

isb=1;

pr=0.0001;

B1=[1 2 (0.065+0.208i)*27/484/2 0.00000364*27*484*2 1 0;

2 3 (0.065+0.208i)*17/484/2 0.00000364*17*484*2 1 0;

2 4 (0.065+0.208i)*6/484/2 0.00000364*6*484*2 1 0;

3 6 (0.065+0.208i)*18/484/2 0.00000364*18*484*2 1 0;

4 5 (0.065+0.208i)*5/484/2 0.00000364*5*484*2 1 0;

5 6 (0.065+0.208i)*24/484/2 0.00000364*24*484*2 1 0];

B2=[0 0 1.0455 1.05 0 1;

0 (250+120i)/100 1 1.05 0 2;

0 (330+190i)/100 1 1.05 0 2;

0 (250+80i)/100 1 1.05 0 2;

0 (400+200i)/100 1 1.05 0 2;

(1+0.6197i)*251.1/100 0 1.05 1.05 0 3];

%?B1 矩阵： 1、 支路首端号； 2、 末端号； 3、 支路阻抗； 4、 支路对地电纳

%?????????5、 支路的变比； 6、 支路首端处于 K 侧为 1， 1 侧为 0

%?B2 矩阵： 1、 该节点发电机功率； 2、 该节点负荷功率； 3、 节点电压初始值

%?????????4、 PV 节点电压 V 的给定值； 5、 节点所接的无功补偿设备的容量

%?????????6、 节点分类标号： 1 为平衡节点(应为 1 号节点)； 2 为 PQ 节点；

%????????????3 为 PV 节点；

isn_1=input('是否在进行n‐1校验，是为1，否为0:');

if isn_1==0

isge=input('请输入发电机出力状态，满载为 1，轻载为 0:');

if isge==0,B2(6,1)=(1+0.6197i)*251.1/100;

else

if isge==1,B2(6,1)=(1+0.6197i)*753.3/100;

end

end

else

if isn_1==1,B2(6,1)=(1+0.6197i)*753.3/100;

nn_1=input('请输入n‐1校验断线序号，从 1‐6:');

B1(nn_1,3)=B1(nn_1,3)*2;

B1(nn_1,4)=B1(nn_1,4)/2;

end

end

Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);

%?%?%‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐

for i=1:nl   %支路数

if B1(i,6)==0   %左节点处于 1 侧

p=B1(i,1);q=B1(i,2);

else  %左节点处于 K 侧

p=B1(i,2);q=B1(i,1);

end

Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));  %非对角元

Y(q,p)=Y(p,q);  %非对角元

Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;   %对角元 K 侧

Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;  %对角元 1 侧

end

%求导纳矩阵

disp('导纳矩阵?Y=');

disp(Y)

%‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐

G=real(Y);B=imag(Y);  %分解出导纳阵的实部和虚部

for i=1:n  %给定各节点初始电压的实部和虚部

e(i)=real(B2(i,3));

f(i)=imag(B2(i,3));

V(i)=B2(i