查看: 724|回复: 10|关注: 0

[已解决] matlab索引超出数组边界问题

[复制链接]

新手

7 麦片

财富积分


050


1

主题

5

帖子

0

最佳答案
本帖最后由 居然然然 于 2019-12-25 16:56 编辑
  1. %本程序的功能是用牛顿——拉夫逊法进行潮流计算
  2. clear;
  3. n=6;
  4. nl=6;
  5. isb=1;
  6. pr=0.0001;
  7. B1=[1 2 (0.065+0.208i)*27/484/2 0.00000364*27*484*2 1 0;
  8.     2 3 (0.065+0.208i)*17/484/2 0.00000364*17*484*2 1 0;
  9. 2 4 (0.065+0.208i)*6/484/2 0.00000364*6*484*2 1 0;
  10. 3 6 (0.065+0.208i)*18/484/2 0.00000364*18*484*2 1 0;
  11. 4 5 (0.065+0.208i)*5/484/2 0.00000364*5*484*2 1 0;
  12. 5 6 (0.065+0.208i)*24/484/2 0.00000364*24*484*2 1 0];
  13. B2=[0 0 1.0455 1.05 0 1;
  14. 0 (250+120i)/100 1 1.05 0 2;
  15. 0 (330+190i)/100 1 1.05 0 2;
  16. 0 (250+80i)/100 1 1.05 0 2;
  17. 0 (400+200i)/100 1 1.05 0 2;
  18. (1+0.6197i)*251.1/100 0 1.05 1.05 0 3];
  19. %?B1 矩阵: 1、 支路首端号; 2、 末端号; 3、 支路阻抗; 4、 支路对地电纳
  20. %?????????5、 支路的变比; 6、 支路首端处于 K 侧为 1, 1 侧为 0
  21. %?B2 矩阵: 1、 该节点发电机功率; 2、 该节点负荷功率; 3、 节点电压初始值
  22. %?????????4、 PV 节点电压 V 的给定值; 5、 节点所接的无功补偿设备的容量
  23. %?????????6、 节点分类标号: 1 为平衡节点(应为 1 号节点); 2 为 PQ 节点;
  24. %????????????3 为 PV 节点;
  25. isn_1=input('是否在进行n‐1校验,是为1,否为0:');
  26. if isn_1==0
  27. isge=input('请输入发电机出力状态,满载为 1,轻载为 0:');
  28. if isge==0,B2(6,1)=(1+0.6197i)*251.1/100;
  29. else
  30.     if isge==1,B2(6,1)=(1+0.6197i)*753.3/100;
  31.     end
  32. end
  33. else
  34. if isn_1==1,B2(6,1)=(1+0.6197i)*753.3/100;
  35. nn_1=input('请输入n‐1校验断线序号,从 1‐6:');
  36. B1(nn_1,3)=B1(nn_1,3)*2;
  37. B1(nn_1,4)=B1(nn_1,4)/2;
  38. end
  39. end
  40. Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);
  41. %?%?%‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐
  42. for i=1:nl   %支路数
  43. if B1(i,6)==0   %左节点处于 1 侧
  44. p=B1(i,1);q=B1(i,2);
  45.     else  %左节点处于 K 侧
  46.     p=B1(i,2);q=B1(i,1);
  47. end
  48. Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5));  %非对角元
  49. Y(q,p)=Y(p,q);  %非对角元
  50. Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2;   %对角元 K 侧
  51. Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2;  %对角元 1 侧
  52. end
  53. %求导纳矩阵
  54. disp('导纳矩阵?Y=');
  55. disp(Y)
  56. %‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐
  57. G=real(Y);B=imag(Y);  %分解出导纳阵的实部和虚部
  58. for i=1:n  %给定各节点初始电压的实部和虚部
  59. e(i)=real(B2(i,3));
  60. f(i)=imag(B2(i,3));
  61. V(i)=B2(i,4);   %PV 节点电压给定模值
  62. end
  63. for i=1:n   %给定各节点注入功率
  64. S(i)=B2(i,1)-B2(i,2);  %i 节点注入功率 SG‐SL??
  65. B(i,i)=B(i,i)+B2(i,5);  %i 节点无功补偿量
  66. end
  67. %===================================================================
  68. P=real(S);Q=imag(S);  %分解出各节点注入的有功和无功功率
  69. ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0;  %迭代次数 ICT1、 a; 不满足收敛要求的节点数 IT2
  70. while IT2~=0  %?N0=2*n?雅可比矩阵的阶数; N=N0+1 扩展列
  71. IT2=0;a=a+1;
  72. for i=1:n
  73.   if i~=isb   %非平衡节点
  74.    C(i)=0;D(i)=0;
  75.    for j1=1:n
  76.     C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ (Gij*ej‐Bij*fj)
  77. D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ (Gij*fj+Bij*ej)
  78.    end
  79. P1=C(i)*e(i)+f(i)*D(i);%节点功率 P 计算 eiΣ (Gij*ej‐Bij*fj)+fiΣ (Gij*fj+Bij*ej)
  80. Q1=C(i)*f(i)-e(i)*D(i);%节点功率 Q 计算 fiΣ (Gij*ej‐Bij*fj)‐eiΣ (Gij*fj+Bij*ej)
  81. %求 i 节点有功和无功功率 P',Q'的计算值
  82. V2=e(i)^2+f(i)^2;   %电压模平方
  83. %=========?以下针对非 PV 节点来求取功率差及 Jacobi 矩阵元素?=========
  84.   if B2(i,6)~=3   %非 PV 节点
  85. DP=P(i)-P1;   %节点有功功率差
  86. DQ=Q(i)-Q1;   %节点无功功率差
  87. %===============?以上为除平衡节点外其它节点的功率计算?=================
  88. %=================?求取 Jacobi 矩阵?===================
  89.    for j1=1:n
  90.     if j1~=isb&j1~=i  %非平衡节点&非对角元
  91. X1=-G(i,j1)*e(i)-B(i,j1)*f(i);  %?dP/de=‐dQ/df
  92. X2=B(i,j1)*e(i)-G(i,j1)*f(i); %?dP/df=dQ/de????
  93. X3=X2; %?X2=dp/df??X3=dQ/de
  94. X4=-X1; %?X1=dP/de??X4=dQ/df
  95. p=2*i-1;q=2*j1-1;
  96. J(p,q)=X3;J(p,N)=DQ;m=p+1; %?X3=dQ/de??J(p,N)=DQ 节点无功功率差
  97. J(m,q)=X1;J(m,N)=DP;q=q+1; %?X1=dP/de??J(m,N)=DP 节点有功功率差
  98. J(p,q)=X4;J(m,q)=X2; %?X4=dQ/df??X2=dp/df
  99.     else
  100.     if j1==i&j1~=isb  %非平衡节点&对角元
  101. X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%?dP/de
  102. X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%?dP/df
  103. X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);%?dQ/de
  104. X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);%?dQ/df
  105. p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q
  106. m=p+1;
  107. J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P
  108. J(m,q)=X2;
  109.     end
  110.     end
  111.    end
  112. %===============?下面是针对 PV 节点来求取 Jacobi 矩阵的元素?===========
  113. DP=P(i)-P1; %?PV 节点有功误差
  114. DV=V(i)^2-V2; %?PV 节点电压误差
  115.    for j1=1:n
  116.     if j1~=isb&j1~=i  %非平衡节点&非对角元
  117. X1=-G(i,j1)*e(i)-B(i,j1)*f(i); %?dP/de
  118. X2=B(i,j1)*e(i)-G(i,j1)*f(i); %?dP/df
  119. X5=0;X6=0;
  120. p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; %?PV 节点电压误差
  121. m=p+1;
  122. J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; %?PV 节点有功误差
  123. J(m,q)=X2;
  124.     else
  125.     if j1==i&j1~=isb  %非平衡节点&对角元
  126. X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%?dP/de
  127. X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%?dP/df
  128. X5=-2*e(i);
  129. X6=-2*f(i);
  130. p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; %?PV 节点电压误差
  131. m=p+1;
  132. J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; %?PV 节点有功误差
  133. J(m,q)=X2;
  134.     end
  135.     end
  136.    end
  137.   end
  138.   end
  139. end
  140. %=========?以 上 为 求 雅 可 比 矩 阵 的 各 个 元 素 及 扩 展 列 的 功 率 差 或 电 压 差===================
  141. for k=3:N0 %?N0=2*n?(从第三行开始, 第一、 二行是平衡节点)
  142.     k1=k+1;N1=N; %?N=N0+1?即?N=2*n+1 扩展列△P、 △Q?或△U
  143.     for k2=k1:N1  %?从 k+1 列的 Jacobi 元素到扩展列的△P、 △Q?或△U??
  144.         J(k,k2)=J(k,k2)./J(k,k);%?用 K 行 K 列对角元素去除 K 行 K 列后的非对角元素进行规格化
  145.     end
  146.     J(k,k)=1; %?对角元规格化 K 行 K 列对角元素赋 1
  147. %====================? 回 代 运 算?================================
  148.     if k~=3 %?不是第三行??k?>?3
  149.         k4=k-1;
  150.         for k3=3:k4  %?用 k3 行从第三行开始到当前行的前一行 k4 行消去
  151.              for k2=k1:N1  %?k3 行后各行上三角元素
  152.               J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算(当前行 k 列元素消为 0)
  153.              end %用当前行 K2 列元素减去当前行 k 列元素乘以第 k 行 K2 列元素
  154. J(k3,k)=0; %当前行第 k 列元素已消为 0
  155.         end
  156.         if k==N0 %若已到最后一行
  157.         break;
  158.         end
  159. %==================?前代运算??==================================
  160.         for k3=k1:N0 %?从 k+1 行到 2*n 最后一行
  161.             for k2=k1:N1 %??从 k+1 列到扩展列消去 k+1 行后各行下三角元素
  162.                 J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算
  163.             end %用当前行 K2 列元素减去当前行 k 列元素乘以第 k 行 K2 列元素
  164.             J(k3,k)=0; %当前行第 k 列元素已消为 0
  165.         end
  166.      else %是第三行 k=3
  167. %======================?第三行 k=3 的前代运算?========================??????
  168.        for k3=k1:N0 %从第四行到 2n 行(最后一行)
  169.            for k2=k1:N1 %从第四列到 2n+1 列(即扩展列)
  170.                J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算(当前行 3 列元素消为 0)
  171.            end %用当前行 K2 列元素减去当前行 3 列元素乘以第三行 K2 列元素
  172.            J(k3,k)=0; %当前行第 3 列元素已消为 0
  173.        end
  174.     end
  175. end
  176. %====上面是用线性变换方式高斯消去法将 Jacobi 矩阵化成单位矩阵=====
  177. for k=3:2:N0-1
  178.           L=(k+1)./2;
  179.           e(L)=e(L)-J(k,N); %修改节点电压实部
  180.           k1=k+1;
  181.           f(L)=f(L)-J(k1,N); %修改节点电压虚部
  182. end
  183.   %‐‐‐‐‐‐修改节点电压‐‐‐‐‐‐‐‐‐‐‐
  184. for k=3:N0
  185.           DET=abs(J(k,N));
  186.           if DET>=pr  %电压偏差量是否满足要求
  187.               IT2=IT2+1; %不满足要求的节点数加 1
  188.           end
  189. end
  190. ICT2(a)=IT2; %不满足要求的节点数
  191. ICT1=ICT1+1; %迭代次数
  192. end
  193. %用高斯消去法解"w=‐J*V"
  194. disp('迭代次数: ');
  195. disp(ICT1);
  196. disp('没有达到精度要求的个数: ');
  197. disp(ICT2);
  198. for k=1:n
  199.      V(k)=sqrt(e(k)^2+f(k)^2); %计算各节点电压的模值
  200. sida(k)=atan(f(k)./e(k))*180./pi; %计算各节点电压的角度
  201.      E(k)=e(k)+f(k)*j; %将各节点电压用复数表示
  202. end
  203. %===============?计算各输出量?===========================
  204. disp('各节点的实际电压标幺值 E 为(节点号从小到大排列): ');
  205. disp(E); %显示各节点的实际电压标幺值 E 用复数表示
  206. disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  207. disp('各节点的电压大小 V 为(节点号从小到大排列): ');
  208. disp(V); %显示各节点的电压大小 V 的模值
  209. disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  210. disp('各节点的电压相角 sida 为(节点号从小到大排列): ');
  211. disp(sida); %显示各节点的电压相角
  212. for p=1:n
  213. C(p)=0;
  214.   for q=1:n
  215.       C(p)=C(p)+conj(Y(p,q))*conj(E(q)); %计算各节点的注入电流的共轭值
  216.   end
  217.   S(p)=E(p)*C(p); %计算各节点的功率?S?=?电压?X?注入电流的共轭值
  218. end
  219. disp('各节点的功率 S 为(节点号从小到大排列): ');
  220. disp(S); %显示各节点的注入功率
  221. disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  222. disp('各条支路的首端功率 Si 为(顺序同您输入 B1 时一致): ');
  223. for i=1:nl
  224.     p=B1(i,1);q=B1(i,2);
  225. if B1(i,6)==0
  226.    Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...
  227.    -conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
  228.    Siz(i)=Si(p,q);
  229. else
  230.    Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...
  231.    -conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
  232.    Siz(i)=Si(p,q);
  233. end
  234. disp(Si(p,q));
  235. SSi(p,q)=Si(p,q);
  236.    ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
  237. disp(ZF);
  238. disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  239. end
  240. disp('各条支路的末端功率 Sj 为(顺序同您输入 B1 时一致): ');
  241. for i=1:nl
  242.     p=B1(i,1);q=B1(i,2);
  243. if B1(i,6)==0
  244.     Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...
  245.     -conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
  246.     Sjy(i)=Sj(q,p);
  247. else
  248.     Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...
  249.     -conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
  250.     Sjy(i)=Sj(q,p);
  251. end
  252. disp(Sj(q,p));
  253. SSj(q,p)=Sj(q,p);
  254.    ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
  255. disp(ZF);
  256. disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  257. end
  258. disp('各条支路的功率损耗 DS 为(顺序同您输入 B1 时一致): ');
  259. for i=1:nl
  260. p=B1(i,1);q=B1(i,2);
  261. DS(i)=Si(p,q)+Sj(q,p);
  262. disp(DS(i));
  263. DDS(i)=DS(i);
  264. ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
  265. disp(ZF);
  266. disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  267. end
  268. xsl=[];
  269. for i=1:nl
  270. disp('线路线损率');%输出线损率
  271. DS1(i)=real(DS(i));
  272. Si1(i)=real(Siz(i));
  273. xsl=[xsl,DS1(i)/Si1(i)];
  274. ZF=['线损率(',num2str(B1(i,1)),'‐',num2str(B1(i,2)),')=',num2str(abs(xsl(i)))];
  275. disp(ZF);
  276. disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  277. end
复制代码


运行后出错显示为
位置 1 处的索引超出数组边界(不能超出 10)。
出错 chaoliujisuan1 (line 170)
               J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算(当前行 3 列元素消为 0)

用的是R2018b、win10系统

论坛优秀回答者

18

主题

2237

帖子

453

最佳答案
  • 关注者: 89
发表于 2019-12-25 10:18:06 | 显示全部楼层
你的代码没办法独立运行,所以没办法复现你的问题,加上你的代码直接放在正文中,没有放在代码框中,给人的感觉就是,你就是随意将问题贴过来了,没有考虑别人帮你解决问题时的感受,在我看来,基本没人会帮你解决。

给你一个建议,在错误的前一句打断点,然后运行看看哪里出问题了。matlab的变量都能直接在工作区中看到,调试就行了。
多看帮助文档
说明你的matlab版本

新手

7 麦片

财富积分


050


1

主题

5

帖子

0

最佳答案
 楼主| 发表于 2019-12-25 16:50:02 | 显示全部楼层
深蓝孩童 发表于 2019-12-25 10:18
你的代码没办法独立运行,所以没办法复现你的问题,加上你的代码直接放在正文中,没有放在代码框中,给人的 ...

抱歉第一次发贴不太懂,谢谢你呀

新手

7 麦片

财富积分


050


1

主题

5

帖子

0

最佳答案
 楼主| 发表于 2019-12-25 16:57:07 | 显示全部楼层
深蓝孩童 发表于 2019-12-25 10:18
你的代码没办法独立运行,所以没办法复现你的问题,加上你的代码直接放在正文中,没有放在代码框中,给人的 ...

我重新编辑了一下,你能帮忙看一下嘛,拜托啦

论坛优秀回答者

18

主题

2237

帖子

453

最佳答案
  • 关注者: 89
发表于 2019-12-25 17:06:23 | 显示全部楼层 |此回复为最佳答案
居然然然 发表于 2019-12-25 16:57
我重新编辑了一下,你能帮忙看一下嘛,拜托啦

你的J是 10x13的,但是你取了第11行,所以必然错误。具体为啥要取到11行,这是你的算法逻辑决定的,你自己决定要怎么更改。
多看帮助文档
说明你的matlab版本

论坛优秀回答者

18

主题

2237

帖子

453

最佳答案
  • 关注者: 89
发表于 2019-12-25 17:09:04 | 显示全部楼层
居然然然 发表于 2019-12-25 16:57
我重新编辑了一下,你能帮忙看一下嘛,拜托啦

还有,你有没有发现,复制你的代码之后,再拷贝到matlab编辑器中,是有很多很多很多空格错误的,你自己尝试一下,你是从之前的帖子中直接复制到代码框中的么?

建议使用下边的格式,ctrl a 再 ctrl i 就行了。
  1. %本程序的功能是用牛顿——拉夫逊法进行潮流计算
  2. clear;
  3. n=6;
  4. nl=6;
  5. isb=1;
  6. pr=0.0001;
  7. B1=[1 2 (0.065+0.208i)*27/484/2 0.00000364*27*484*2 1 0;
  8.     2 3 (0.065+0.208i)*17/484/2 0.00000364*17*484*2 1 0;
  9.     2 4 (0.065+0.208i)*6/484/2 0.00000364*6*484*2 1 0;
  10.     3 6 (0.065+0.208i)*18/484/2 0.00000364*18*484*2 1 0;
  11.     4 5 (0.065+0.208i)*5/484/2 0.00000364*5*484*2 1 0;
  12.     5 6 (0.065+0.208i)*24/484/2 0.00000364*24*484*2 1 0];
  13. B2=[0 0 1.0455 1.05 0 1;
  14.     0 (250+120i)/100 1 1.05 0 2;
  15.     0 (330+190i)/100 1 1.05 0 2;
  16.     0 (250+80i)/100 1 1.05 0 2;
  17.     0 (400+200i)/100 1 1.05 0 2;
  18.     (1+0.6197i)*251.1/100 0 1.05 1.05 0 3];
  19. %?B1 矩阵: 1、 支路首端号; 2、 末端号; 3、 支路阻抗; 4、 支路对地电纳
  20. %?????????5、 支路的变比; 6、 支路首端处于 K 侧为 1, 1 侧为 0
  21. %?B2 矩阵: 1、 该节点发电机功率; 2、 该节点负荷功率; 3、 节点电压初始值
  22. %?????????4、 PV 节点电压 V 的给定值; 5、 节点所接的无功补偿设备的容量
  23. %?????????6、 节点分类标号: 1 为平衡节点(应为 1 号节点); 2 为 PQ 节点;
  24. %????????????3 为 PV 节点;
  25. isn_1=input('是否在进行n‐1校验,是为1,否为0:');
  26. if isn_1==0
  27.     isge=input('请输入发电机出力状态,满载为 1,轻载为 0:');
  28.     if isge==0,B2(6,1)=(1+0.6197i)*251.1/100;
  29.     else
  30.         if isge==1,B2(6,1)=(1+0.6197i)*753.3/100;
  31.         end
  32.     end
  33. else
  34.     if isn_1==1,B2(6,1)=(1+0.6197i)*753.3/100;
  35.         nn_1=input('请输入n‐1校验断线序号,从 1‐6:');
  36.         B1(nn_1,3)=B1(nn_1,3)*2;
  37.         B1(nn_1,4)=B1(nn_1,4)/2;
  38.     end
  39. end
  40. Y=zeros(n);e=zeros(1,n);f=zeros(1,n);V=zeros(1,n);sida=zeros(1,n);S1=zeros(nl);
  41. %?%?%‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐
  42. for i=1:nl %支路数
  43.     if B1(i,6)==0 %左节点处于 1 侧
  44.         p=B1(i,1);q=B1(i,2);
  45.     else %左节点处于 K 侧
  46.         p=B1(i,2);q=B1(i,1);
  47.     end
  48.     Y(p,q)=Y(p,q)-1./(B1(i,3)*B1(i,5)); %非对角元
  49.     Y(q,p)=Y(p,q); %非对角元
  50.     Y(q,q)=Y(q,q)+1./(B1(i,3)*B1(i,5)^2)+B1(i,4)./2; %对角元 K 侧
  51.     Y(p,p)=Y(p,p)+1./B1(i,3)+B1(i,4)./2; %对角元 1 侧
  52. end
  53. %求导纳矩阵
  54. disp('导纳矩阵?Y=');
  55. disp(Y)
  56. %‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐
  57. G=real(Y);B=imag(Y); %分解出导纳阵的实部和虚部
  58. for i=1:n %给定各节点初始电压的实部和虚部
  59.     e(i)=real(B2(i,3));
  60.     f(i)=imag(B2(i,3));
  61.     V(i)=B2(i,4); %PV 节点电压给定模值
  62. end
  63. for i=1:n %给定各节点注入功率
  64.     S(i)=B2(i,1)-B2(i,2); %i 节点注入功率 SG‐SL??
  65.     B(i,i)=B(i,i)+B2(i,5); %i 节点无功补偿量
  66. end
  67. %===================================================================
  68. P=real(S);Q=imag(S); %分解出各节点注入的有功和无功功率
  69. ICT1=0;IT2=1;N0=2*n;N=N0+1;a=0; %迭代次数 ICT1、 a; 不满足收敛要求的节点数 IT2
  70. while IT2~=0 %?N0=2*n?雅可比矩阵的阶数; N=N0+1 扩展列
  71.     IT2=0;a=a+1;
  72.     for i=1:n
  73.         if i~=isb %非平衡节点
  74.             C(i)=0;D(i)=0;
  75.             for j1=1:n
  76.                 C(i)=C(i)+G(i,j1)*e(j1)-B(i,j1)*f(j1);%Σ (Gij*ej‐Bij*fj)
  77.                 D(i)=D(i)+G(i,j1)*f(j1)+B(i,j1)*e(j1);%Σ (Gij*fj+Bij*ej)
  78.             end
  79.             P1=C(i)*e(i)+f(i)*D(i);%节点功率 P 计算 eiΣ (Gij*ej‐Bij*fj)+fiΣ (Gij*fj+Bij*ej)
  80.             Q1=C(i)*f(i)-e(i)*D(i);%节点功率 Q 计算 fiΣ (Gij*ej‐Bij*fj)‐eiΣ (Gij*fj+Bij*ej)
  81.             %求 i 节点有功和无功功率 P',Q'的计算值
  82.             V2=e(i)^2+f(i)^2; %电压模平方
  83.             %=========?以下针对非 PV 节点来求取功率差及 Jacobi 矩阵元素?=========
  84.             if B2(i,6)~=3 %非 PV 节点
  85.                 DP=P(i)-P1; %节点有功功率差
  86.                 DQ=Q(i)-Q1; %节点无功功率差
  87.                 %===============?以上为除平衡节点外其它节点的功率计算?=================
  88.                 %=================?求取 Jacobi 矩阵?===================
  89.                 for j1=1:n
  90.                     if j1~=isb&j1~=i %非平衡节点&非对角元
  91.                         X1=-G(i,j1)*e(i)-B(i,j1)*f(i); %?dP/de=‐dQ/df
  92.                         X2=B(i,j1)*e(i)-G(i,j1)*f(i); %?dP/df=dQ/de????
  93.                         X3=X2; %?X2=dp/df??X3=dQ/de
  94.                         X4=-X1; %?X1=dP/de??X4=dQ/df
  95.                         p=2*i-1;q=2*j1-1;
  96.                         J(p,q)=X3;J(p,N)=DQ;m=p+1; %?X3=dQ/de??J(p,N)=DQ 节点无功功率差
  97.                         J(m,q)=X1;J(m,N)=DP;q=q+1; %?X1=dP/de??J(m,N)=DP 节点有功功率差
  98.                         J(p,q)=X4;J(m,q)=X2; %?X4=dQ/df??X2=dp/df
  99.                     else
  100.                         if j1==i&j1~=isb %非平衡节点&对角元
  101.                             X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%?dP/de
  102.                             X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%?dP/df
  103.                             X3=D(i)+B(i,i)*e(i)-G(i,i)*f(i);%?dQ/de
  104.                             X4=-C(i)+G(i,i)*e(i)+B(i,i)*f(i);%?dQ/df
  105.                             p=2*i-1;q=2*j1-1;J(p,q)=X3;J(p,N)=DQ;%扩展列△Q
  106.                             m=p+1;
  107.                             J(m,q)=X1;q=q+1;J(p,q)=X4;J(m,N)=DP;%扩展列△P
  108.                             J(m,q)=X2;
  109.                         end
  110.                     end
  111.                 end
  112.                 %===============?下面是针对 PV 节点来求取 Jacobi 矩阵的元素?===========
  113.                 DP=P(i)-P1; %?PV 节点有功误差
  114.                 DV=V(i)^2-V2; %?PV 节点电压误差
  115.                 for j1=1:n
  116.                     if j1~=isb&j1~=i %非平衡节点&非对角元
  117.                         X1=-G(i,j1)*e(i)-B(i,j1)*f(i); %?dP/de
  118.                         X2=B(i,j1)*e(i)-G(i,j1)*f(i); %?dP/df
  119.                         X5=0;X6=0;
  120.                         p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; %?PV 节点电压误差
  121.                         m=p+1;
  122.                         J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; %?PV 节点有功误差
  123.                         J(m,q)=X2;
  124.                     else
  125.                         if j1==i&j1~=isb %非平衡节点&对角元
  126.                             X1=-C(i)-G(i,i)*e(i)-B(i,i)*f(i);%?dP/de
  127.                             X2=-D(i)+B(i,i)*e(i)-G(i,i)*f(i);%?dP/df
  128.                             X5=-2*e(i);
  129.                             X6=-2*f(i);
  130.                             p=2*i-1;q=2*j1-1;J(p,q)=X5;J(p,N)=DV; %?PV 节点电压误差
  131.                             m=p+1;
  132.                             J(m,q)=X1;J(m,N)=DP;q=q+1;J(p,q)=X6; %?PV 节点有功误差
  133.                             J(m,q)=X2;
  134.                         end
  135.                     end
  136.                 end
  137.             end
  138.         end
  139.     end
  140.     %=========?以 上 为 求 雅 可 比 矩 阵 的 各 个 元 素 及 扩 展 列 的 功 率 差 或 电 压 差===================
  141.     for k=3:N0 %?N0=2*n?(从第三行开始, 第一、 二行是平衡节点)
  142.         k1=k+1;N1=N; %?N=N0+1?即?N=2*n+1 扩展列△P、 △Q?或△U
  143.         for k2=k1:N1 %?从 k+1 列的 Jacobi 元素到扩展列的△P、 △Q?或△U??
  144.             J(k,k2)=J(k,k2)./J(k,k);%?用 K 行 K 列对角元素去除 K 行 K 列后的非对角元素进行规格化
  145.         end
  146.         J(k,k)=1; %?对角元规格化 K 行 K 列对角元素赋 1
  147.         %====================? 回 代 运 算?================================
  148.         if k~=3 %?不是第三行??k?>?3
  149.             k4=k-1;
  150.             for k3=3:k4 %?用 k3 行从第三行开始到当前行的前一行 k4 行消去
  151.                 for k2=k1:N1 %?k3 行后各行上三角元素
  152.                     J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算(当前行 k 列元素消为 0)
  153.                 end %用当前行 K2 列元素减去当前行 k 列元素乘以第 k 行 K2 列元素
  154.                 J(k3,k)=0; %当前行第 k 列元素已消为 0
  155.             end
  156.             if k==N0 %若已到最后一行
  157.                 break;
  158.             end
  159.             %==================?前代运算??==================================
  160.             for k3=k1:N0 %?从 k+1 行到 2*n 最后一行
  161.                 for k2=k1:N1 %??从 k+1 列到扩展列消去 k+1 行后各行下三角元素
  162.                     J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算
  163.                 end %用当前行 K2 列元素减去当前行 k 列元素乘以第 k 行 K2 列元素
  164.                 J(k3,k)=0; %当前行第 k 列元素已消为 0
  165.             end
  166.         else %是第三行 k=3
  167.             %======================?第三行 k=3 的前代运算?========================??????
  168.             for k3=k1:N0 %从第四行到 2n 行(最后一行)
  169.                 for k2=k1:N1 %从第四列到 2n+1 列(即扩展列)
  170.                     J(k3,k2)=J(k3,k2)-J(k3,k)*J(k,k2);%消去运算(当前行 3 列元素消为 0)
  171.                 end %用当前行 K2 列元素减去当前行 3 列元素乘以第三行 K2 列元素
  172.                 J(k3,k)=0; %当前行第 3 列元素已消为 0
  173.             end
  174.         end
  175.     end
  176.     %====上面是用线性变换方式高斯消去法将 Jacobi 矩阵化成单位矩阵=====
  177.     for k=3:2:N0-1
  178.         L=(k+1)./2;
  179.         e(L)=e(L)-J(k,N); %修改节点电压实部
  180.         k1=k+1;
  181.         f(L)=f(L)-J(k1,N); %修改节点电压虚部
  182.     end
  183.     %‐‐‐‐‐‐修改节点电压‐‐‐‐‐‐‐‐‐‐‐
  184.     for k=3:N0
  185.         DET=abs(J(k,N));
  186.         if DET>=pr %电压偏差量是否满足要求
  187.             IT2=IT2+1; %不满足要求的节点数加 1
  188.         end
  189.     end
  190.     ICT2(a)=IT2; %不满足要求的节点数
  191.     ICT1=ICT1+1; %迭代次数
  192. end
  193. %用高斯消去法解"w=‐J*V"
  194. disp('迭代次数: ');
  195. disp(ICT1);
  196. disp('没有达到精度要求的个数: ');
  197. disp(ICT2);
  198. for k=1:n
  199.     V(k)=sqrt(e(k)^2+f(k)^2); %计算各节点电压的模值
  200.     sida(k)=atan(f(k)./e(k))*180./pi; %计算各节点电压的角度
  201.     E(k)=e(k)+f(k)*j; %将各节点电压用复数表示
  202. end
  203. %===============?计算各输出量?===========================
  204. disp('各节点的实际电压标幺值 E 为(节点号从小到大排列): ');
  205. disp(E); %显示各节点的实际电压标幺值 E 用复数表示
  206. disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  207. disp('各节点的电压大小 V 为(节点号从小到大排列): ');
  208. disp(V); %显示各节点的电压大小 V 的模值
  209. disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  210. disp('各节点的电压相角 sida 为(节点号从小到大排列): ');
  211. disp(sida); %显示各节点的电压相角
  212. for p=1:n
  213.     C(p)=0;
  214.     for q=1:n
  215.         C(p)=C(p)+conj(Y(p,q))*conj(E(q)); %计算各节点的注入电流的共轭值
  216.     end
  217.     S(p)=E(p)*C(p); %计算各节点的功率?S?=?电压?X?注入电流的共轭值
  218. end
  219. disp('各节点的功率 S 为(节点号从小到大排列): ');
  220. disp(S); %显示各节点的注入功率
  221. disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  222. disp('各条支路的首端功率 Si 为(顺序同您输入 B1 时一致): ');
  223. for i=1:nl
  224.     p=B1(i,1);q=B1(i,2);
  225.     if B1(i,6)==0
  226.         Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)*B1(i,5))...
  227.             -conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
  228.         Siz(i)=Si(p,q);
  229.     else
  230.         Si(p,q)=E(p)*(conj(E(p))*conj(B1(i,4)./2)+(conj(E(p)./B1(i,5))...
  231.             -conj(E(q)))*conj(1./(B1(i,3)*B1(i,5))));
  232.         Siz(i)=Si(p,q);
  233.     end
  234.     disp(Si(p,q));
  235.     SSi(p,q)=Si(p,q);
  236.     ZF=['S(',num2str(p),',',num2str(q),')=',num2str(SSi(p,q))];
  237.     disp(ZF);
  238.     disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  239. end
  240. disp('各条支路的末端功率 Sj 为(顺序同您输入 B1 时一致): ');
  241. for i=1:nl
  242.     p=B1(i,1);q=B1(i,2);
  243.     if B1(i,6)==0
  244.         Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)./B1(i,5))...
  245.             -conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
  246.         Sjy(i)=Sj(q,p);
  247.     else
  248.         Sj(q,p)=E(q)*(conj(E(q))*conj(B1(i,4)./2)+(conj(E(q)*B1(i,5))...
  249.             -conj(E(p)))*conj(1./(B1(i,3)*B1(i,5))));
  250.         Sjy(i)=Sj(q,p);
  251.     end
  252.     disp(Sj(q,p));
  253.     SSj(q,p)=Sj(q,p);
  254.     ZF=['S(',num2str(q),',',num2str(p),')=',num2str(SSj(q,p))];
  255.     disp(ZF);
  256.     disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  257. end
  258. disp('各条支路的功率损耗 DS 为(顺序同您输入 B1 时一致): ');
  259. for i=1:nl
  260.     p=B1(i,1);q=B1(i,2);
  261.     DS(i)=Si(p,q)+Sj(q,p);
  262.     disp(DS(i));
  263.     DDS(i)=DS(i);
  264.     ZF=['DS(',num2str(p),',',num2str(q),')=',num2str(DDS(i))];
  265.     disp(ZF);
  266.     disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  267. end
  268. xsl=[];
  269. for i=1:nl
  270.     disp('线路线损率');%输出线损率
  271.     DS1(i)=real(DS(i));
  272.     Si1(i)=real(Siz(i));
  273.     xsl=[xsl,DS1(i)/Si1(i)];
  274.     ZF=['线损率(',num2str(B1(i,1)),'‐',num2str(B1(i,2)),')=',num2str(abs(xsl(i)))];
  275.     disp(ZF);
  276.     disp('‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐');
  277. end
复制代码
多看帮助文档
说明你的matlab版本

论坛优秀回答者

18

主题

2237

帖子

453

最佳答案
  • 关注者: 89
发表于 2019-12-25 17:09:47 | 显示全部楼层
居然然然 发表于 2019-12-25 16:57
我重新编辑了一下,你能帮忙看一下嘛,拜托啦

运行代码之前,先运行一句话,dbstop if error,然后就会停在运行错误的地方了,然后你就可以看变量了。
多看帮助文档
说明你的matlab版本

新手

7 麦片

财富积分


050


1

主题

5

帖子

0

最佳答案
 楼主| 发表于 2019-12-25 17:25:04 | 显示全部楼层
深蓝孩童 发表于 2019-12-25 17:09
还有,你有没有发现,复制你的代码之后,再拷贝到matlab编辑器中,是有很多很多很多空格错误的,你自己尝 ...

我是直接从m文件编辑框里复制的,不知道为什么这样,谢谢你,我去试啦!

新手

7 麦片

财富积分


050


1

主题

5

帖子

0

最佳答案
 楼主| 发表于 2019-12-25 18:14:42 | 显示全部楼层
深蓝孩童 发表于 2019-12-25 17:06
你的J是 10x13的,但是你取了第11行,所以必然错误。具体为啥要取到11行,这是你的算法逻辑决定的,你自 ...

那个,我想问问什么时候取到了第11行啊,我手写了一遍没发现呀= =

论坛优秀回答者

18

主题

2237

帖子

453

最佳答案
  • 关注者: 89
发表于 2019-12-25 19:20:38 | 显示全部楼层
居然然然 发表于 2019-12-25 18:14
那个,我想问问什么时候取到了第11行啊,我手写了一遍没发现呀= =

那就很令人好奇了,代码框中的竟然还有中文空格。

手写? 那是代码的问题了。你的N0是12,循环中肯定k3可以达到11,那为什么代码这么写呢。。。为什么会问我什么时候取到了11行。。
多看帮助文档
说明你的matlab版本
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

快速回复 返回顶部 返回列表