IEEE33节点前推回代程序

clc
clear all;
format short;
Z=[0  1  0.0922+j*0.047   100+j*60
   1  2  0.4930+j*0.2511  90+j*40
   2  3  0.3660+j*0.1864  120+j*80
   3  4  0.3811+j*0.1941  60+j*30 
   4  5  0.8190+j*0.7070  60+j*20
   5  6  0.1872+j*0.6188  200+j*100
   6  7  0.7114+j*0.2351  200+j*100
   7  8  1.0300+j*0.7400  60+j*20
   8  9  1.0440+j*0.7400  60+j*20
   9  10 0.1966+j*0.0650  45+j*30
   10 11 0.3744+j*0.1238  60+j*35
   11 12 1.4680+j*1.1550  60+j*35
   12 13 0.5416+j*0.7129  120+j*80
   13 14 0.5910+j*0.5260  60+j*10
   14 15 0.7463+j*0.5450  60+j*20
   15 16 1.2890+j*1.7210  60+j*20
   16 17 0.3720+j*0.5740  90+j*40
   1  18 0.1640+j*0.1565  90+j*40
   18 19 1.5042+j*1.3554  90+j*40
   19 20 0.4095+j*0.4784  90+j*40
   20 21 0.7089+j*0.9373  90+j*40
   2  22 0.4512+j*0.3083  90+j*50
   22 23 0.8980+j*0.7091  420+j*200
   23 24 0.8960+j*0.7011  420+j*200
   5  25 0.2030+j*0.1034  60+j*25
   25 26 0.2842+j*0.1447  60+j*25
   26 27 1.0590+j*0.9337  60+j*20
   27 28 0.8042+j*0.7006  120+j*70
   28 29 0.5075+j*0.2585  200+j*600
   29 30 0.9744+j*0.9630  150+j*70
   30 31 0.3105+j*0.3619  210+j*100
   31 32 0.3410+j*0.5362  60+j*40];%根据网络拓扑结构，第一列为支路首节点，第二列为支路末节点，第三列为支路阻抗，第四列为节点负荷功率
%将参数化为标幺值进行计算（Sb=10MVA;Vb=12.66KV）
Sb=10;  
Vb=12.66;
Z（:4）=Z（:4）/（Sb*1000）;
Z（:3）=Z（:3）*Sb/Vb^2;

                       %Start
[mn]=size（Z）;
V=ones（m+11）;
dU=ones（m1）;
S=zeros（m+11）;
Sr=Z（:4）;
dSr=zeros（m1）;