%%% 單位:電流單位為A,電壓單位為kV，功率單位為kW，阻抗單位為歐姆
function [U,iterations,error,costime]=formback
t1=clock;
bus16;%% 導出networkdata文件中的data矩陣,和網絡的輸入端的額定電壓為Un，網絡的參數都存儲在data中
%%global A;
global Ibr;
A=data;  %%% 輸入參數矩陣A說明：矩陣A的行數為配電網的支路數，共6列，第1列
[Rower,Column]=size(A);
Nodes=Rower+1;
Zbr=A(:,3)+j*A(:,4);
Srnode=A(:,5)+j*A(:,6);
fork=0;%% 定義fork為存儲分叉節點的行向量，其初始值為一個數0，即行向量fork的第一列為0

%% 搜索分叉節點，并將所有的分叉節點存在向量fork里面，除第一列外
%% 搜索分叉節點，并將所有的分叉節點存在向量fork里面，除第一列外
for ii=1:Rower-1
    k=0;
    [r,fsize]=size(fork); %% 取行向量fork的長度為fsize
    for jj=ii+1:Rower
        if A(ii,1)==A(jj,1)
            %% 判斷分叉節點A(i,1)是否已經被保存過？
            for m=1:fsize
                if A(jj,1)==fork(m)
                    k=1;  %% k=1，表示分叉點已經被保存過，k=0,表示沒有保存過
                end
            end
            %% 判斷分叉節點A(i,1)是否已經被保存過？
            if k==0
                fork=[fork,A(ii,1)]; %% 保存分叉節點A(i,1)
            end
            break  %% 當搜索到分叉點后立即退出搜索,返回重新搜索
        end
    end
end
%% 搜索分叉節點，并將所有的分叉節點存在向量fork里面，除第一列外
%% 搜索分叉節點，并將所有的分叉節點存在向量fork里面，除第一列外

%% 搜索末節點
%% 搜索末節點
End=[0,0]; %% 定義End為存儲末端節點的兩列矩陣，第一列為末端節點所在的行,第二列為末端節點編號;其初始值為[0,0]，即End的第一行為0
for ii=1:Rower
    k=0;
    for jj=1:Rower
        if A(ii,2)==A(jj,1)
            k=1; %% 搜索有沒有與受端節點相同的送端節點,如果有使k=1,說明該受端節點不是末節點
            break
        end
    end
    if k==0
        End=[End;[ii,A(ii,2)]]; %% k=0,表示受端節點是末節點，將A(ii,2)存儲在End行向量中
    end
end
%% 搜索末節點
%% 搜索末節點

[r,fsize]=size(fork);[esize,c]=size(End);
forknods=fork(2:fsize); %% 取出所有的分叉節點存放在行向量forknods中
endnods=End(2:esize,:); %% 取出所有的末端節點存放在矩陣endnods中,第一列表示末端節點所在支路在矩陣A中的行數,第二列表示末端節點的編號
%%global Unods;
Unods(Rower,3)=0; 
Unods(:,1:2)=A(:,1:2); %% Unods第一列表示支路的送端節點編號,第二列表示受端節點編號,第三列表示受端節點電壓
Unods(:,3)=Un/Ubse; %% 定義所有受端節點的電壓的初始值均為額定電壓值Un,Ubse為基準值

accuracy=1.0e-8;
error=Un;
iterations=0;

while error>accuracy
      iterations=iterations+1;
      Uabs1=abs(Unods(:,3));
%% 利用前推法求解各支路的電流值
%% 利用前推法求解各支路的電流值
      Inlp(Rower,3)=0;  %% 定義節點負荷電流存儲矩陣
      Inlp(:,1:2)=A(:,1:2); %% Inlp的第一列存儲送端節點的編號,第二列存儲受端節點的編號,第三列存儲受端節點的負荷電流(初始值為0)
      Inlp(:,3)=conj(Srnode./Unods(:,3)); %% 求解受端節點的負荷電流

      [r,fsize]=size(forknods);
      [endsize,c]=size(endnods);

      Ibr(Rower,3)=0;  %% 定義支路負荷電流存儲矩陣
      Ibr(:,1:2)=A(:,1:2);  %% Ibr的第一列存儲送端節點的編號,第二列存儲受端節點的編號,第三列存儲支路的負荷電流(初始值為0)
      fork1=0;
      for ii=1:endsize
          Ibr(endnods(ii,1),3)=Inlp(endnods(ii,1),3); %% 將末端節點的負荷電流賦給所連分支的分支電流
          Sendnod=Ibr(endnods(ii,1),1); %% Sendnod為之路的送端節點編號，它隨之路的不同而變化的
          I=Ibr(endnods(ii,1),3);  %% 定義臨時支路電流，它也隨之路的不同而變化的
          k=0;
          kk=0;
          while kk<=Rower
                kk=kk+1;
                for jj=1:fsize
                    %%if Ibr(Sendnod,1)==forknods(jj)
                    if Sendnod==forknods(jj)
                       k=1;  %% 如果支路送端節點Sendnod是分叉節點的話，則賦給k=1,并退出本層循環體
                       break
                    end
                end
                if k==1
                    i1=1;
                    [r7,fsize1]=size(fork1);
                    for mm=1:fsize1
                        if Sendnod==fork1(mm)
                            i1=i1+1;
                        end
                    end
                    if i1==1
                        fork1=[fork1,Sendnod];
                    end
                    break   %% 如果k等于1,則表示支路的送端節點為分叉節點,退出本層循環體
                end
                %% 尋找受端節點是節點Sendnod的支路
                if Ibr(kk,2)==Sendnod
                   Ibr(kk,3)=I+Inlp(kk,3); %% 受端節點是節點Sendnod的支路的電流=送端節點是Sendnod的支路的電流+受端節點Sendnod處的負荷電流
                   Sendnod=Ibr(kk,1); %% 更新送端節點編號Sendnod為新支路的送端節點編號
                   I=Ibr(kk,3); %% 更新支路電流為新支路的支路電流
                   kk=0;
                end
                %% 尋找受端節點是節點Sendnod的支路
          end
      end
%% 由于在節點編號時有意將分叉點由小到大排列，由大的分叉點向小的分叉點聚攏，最后聚攏到一條支路上
      [r6,fsize1]=size(fork1);
      fork1=fork1(:,2:fsize1);;
      [r6,fsize1]=size(fork1);
      for n=1:fsize1             %% 交叉點由大到小循環
          Sendnod=fork1(n);
          I=0;
          for m=1:Rower
              if Ibr(m,1)==Sendnod  %% 尋找送端節點是交叉節點的支路
                 I=I+Ibr(m,3); %% 計算從交叉點流向其它分支的的總的電流，不包括節點Sendnod的負荷電流
              end
          end
          k=0;
          kk=0;
          while kk<Rower
                kk=kk+1;
                if Ibr(kk,2)==Sendnod  
                   Ibr(kk,3)=I+Inlp(kk,3); %% 計算受端節點為Sendnod的支路的電流 
                   Sendnod=Ibr(kk,1);
                   kk=0;
                   for nn=1:fsize1
                       if Sendnod==fork1(nn)
                          k=1; %% 如果受端節點為Sendnod的支路的送端節點編號是交叉點，則退出
                          break
                       end
                   end
                end
                if k==1
                   break
                end
          end
      end
%% 利用前推法求解各支路的電流值
%% 利用前推法求解各支路的電流值

%% 利用回代法求解各支路的電流值
%% 利用回代法求解各支路的電流值            
      Sendnode=A(1,1); %Ufn=Un; %% Sendnode為首端節點編號，Ufn為首端節點電壓，一般為額定值不變
      %Uupnode=Ufn;
      for ii=1:Rower               %%有問題
          %if Ibr(ii,1)==Sendnode
              %Unods(ii,3)=Un;
          %end
          for jj=1:Rower
              if Ibr(jj,1)==Sendnode
                  Unods(jj,3)=Un/Ubse-Zbr(jj)*Ibr(jj,3);
              elseif Ibr(jj,1)==Ibr(ii,2)
                  Unods(jj,3)=Unods(ii,3)-Zbr(jj)*Ibr(jj,3);
              end
          end
      end   
      %%% 利用回代法求解各支路的電流值
      %%% 利用回代法求解各支路的電流值 
      Uabs2=abs(Unods(:,3));
      error=max(Uabs1-Uabs2);
end
%% 利用最終計算所得到的節點電壓值重新計算各分支電流
%% 利用最終計算所得到的節點電壓值重新計算各分支電流
Inlp(Rower,3)=0;  %% 定義節點負荷電流存儲矩陣
      Inlp(:,1:2)=A(:,1:2); %% Inlp的第一列存儲送端節點的編號,第二列存儲受端節點的編號,第三列存儲受端節點的負荷電流(初始值為0)
      Inlp(:,3)=conj(Srnode./Unods(:,3)); %% 求解受端節點的負荷電流

      [r,fsize]=size(forknods);
      [endsize,c]=size(endnods);
      Ibr(Rower,3)=0;  %% 定義支路負荷電流存儲矩陣
      Ibr(:,1:2)=A(:,1:2);  %% Ibr的第一列存儲送端節點的編號,第二列存儲受端節點的編號,第三列存儲支路的負荷電流(初始值為0)

      for ii=1:endsize
          Ibr(endnods(ii,1),3)=Inlp(endnods(ii,1),3); %% 將末端節點的負荷電流賦給所連分支的分支電流
          Sendnod=Ibr(endnods(ii,1),1); %% Sendnod為之路的送端節點編號，它隨之路的不同而變化的
          I=Ibr(endnods(ii,1),3);  %% 定義臨時支路電流，它也隨之路的不同而變化的
          k=0;
          kk=0;
          while kk<=Rower
                kk=kk+1;
                for jj=1:fsize
                   %% if Ibr(Sendnod,1)==forknods(jj)
                    if Sendnod==forknods(jj)
                       k=1;  %% 如果支路送端節點Sendnod是分叉節點的話，則賦給k=1,并退出本層循環體
                       break
                    end
                end
                if k==1
                   break   %% 如果k等于1,則表示支路的送端節點為分叉節點,退出本層循環體
                end
                %% 尋找受端節點是節點Sendnod的支路
                if Ibr(kk,2)==Sendnod
                   Ibr(kk,3)=I+Inlp(kk,3); %% 受端節點是節點Sendnod的支路的電流=送端節點是Sendnod的支路的電流+受端節點Sendnod處的負荷電流
                   Sendnod=Ibr(kk,1); %% 更新送端節點編號Sendnod為新支路的送端節點編號
                   I=Ibr(kk,3); %% 更新支路電流為新支路的支路電流
                   kk=0;
                end
                %% 尋找受端節點是節點Sendnod的支路
          end
      end

%% 由于在節點編號時有意將分叉點由小到大排列，由大的分叉點向小的分叉點聚攏，最后聚攏到一條支路上

      for n=1:fsize1             %% 交叉點由大到小循環
          Sendnod=fork1(n);
          I=0;
          for m=1:Rower
              if Ibr(m,1)==Sendnod  %% 尋找送端節點是交叉節點的支路
                 I=I+Ibr(m,3); %% 計算從交叉點流向其它分支的的總的電流，不包括節點Sendnod的負荷電流
              end
          end
          k=0;
          kk=0;
          while kk<Rower
                kk=kk+1;
                if Ibr(kk,2)==Sendnod  
                   Ibr(kk,3)=I+Inlp(kk,3); %% 計算受端節點為Sendnod的支路的電流 
                   Sendnod=Ibr(kk,1);
                   kk=0;
                   for nn=1:fsize1
                       if Sendnod==fork1(nn)
                          k=1; %% 如果受端節點為Sendnod的支路的送端節點編號是交叉點，則退出
                          break
                       end
                   end
                end
                if k==1
                   break
                end
          end
      end
%% 利用最終計算所得到的節點電壓值重新計算各分支電流
%% 利用最終計算所得到的節點電壓值重新計算各分支電流

U(:,1:2)=Unods(:,1:2);
U(:,3)=abs(Unods(:,3));
U(:,4)=angle(Unods(:,3));
fork1

%% 計算各支路的功率損耗
%% 計算各支路的功率損耗
IIbr=(Ibr(:,3).*conj(Ibr(:,3)));
Sloss=IIbr.*Zbr*Sbse*10^6; %% 單位為W
%% 計算各支路的功率損耗
%% 計算各支路的功率損耗
t2=clock;
t=t2-t1;
costime=t;
return