Home > matpower7.1 > extras > reduction > PartialNumLU.m

PartialNumLU

PURPOSE ^

FUNCTION SUMMARY:

SYNOPSIS ^

function [DataEQ,DataShunt] = PartialNumLU (CIndx,CIndxU,Data,dim,ERP,ERPU,stop,ERPEQ,CIndxEQ,BoundBus)

DESCRIPTION ^

FUNCTION SUMMARY:
   Subroutine PartialNumLU does partial numerical LU factorization to
   given full model bus addmittance matrix and calculate the equivalent
   branch reactance and the equivalent shunts (generated in the
   factorization process) added to the boundary buses.
  
   [DataEQ,DataShunt] = PartialNumLU (CIndx,CIndxU,Data,dim,ERP,ERPU,stop,ERPEQ,CIndxEQ,BoundBus)

 INPUT DATA:
   CIndx - N*1 array containning the column index of the rows, N is the
           number of non-zeros elements in the input data
   CIndxU - N*1 array containing the column index of off diagonal element
            in each row in the U matrix. The length N depends on the
            number of native plus filled non-zero elements in the off
            diagonal position in U matrix. CIndxU is unordered.
   Data -N*1 array containing the data of matrix element in the original
         input file.
   dim - scalar, dimension of the input matrix
   ERPU -   N_dim*1 array containing end of row pointer of all rows except
            the last row which doesn have any off diagonal element, N_dim is
            the dimension of the input matrix A in the original Ax = b
            problem
   ERP - (N_node+1)*1 array containning the end of row pointer data
   stop - scalar, equal to the number of external buses
   ERPEQ, CIndxEQ, 1*n arrays, together build the pointers of the
       equivalent branches
   BoundBus - 1*n array, list of boundary buses

 OUTPUT DATA:
   DataEQ: 1*n array, includes reactance value of the equivalent branches
   DataShunt: 1*n array, includes equivalent bus shunts of the boundary
   buses

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [DataEQ,DataShunt] = PartialNumLU (CIndx,CIndxU,Data,dim,ERP,ERPU,stop,ERPEQ,CIndxEQ,BoundBus)
0002 %FUNCTION SUMMARY:
0003 %   Subroutine PartialNumLU does partial numerical LU factorization to
0004 %   given full model bus addmittance matrix and calculate the equivalent
0005 %   branch reactance and the equivalent shunts (generated in the
0006 %   factorization process) added to the boundary buses.
0007 %
0008 %   [DataEQ,DataShunt] = PartialNumLU (CIndx,CIndxU,Data,dim,ERP,ERPU,stop,ERPEQ,CIndxEQ,BoundBus)
0009 %
0010 % INPUT DATA:
0011 %   CIndx - N*1 array containning the column index of the rows, N is the
0012 %           number of non-zeros elements in the input data
0013 %   CIndxU - N*1 array containing the column index of off diagonal element
0014 %            in each row in the U matrix. The length N depends on the
0015 %            number of native plus filled non-zero elements in the off
0016 %            diagonal position in U matrix. CIndxU is unordered.
0017 %   Data -N*1 array containing the data of matrix element in the original
0018 %         input file.
0019 %   dim - scalar, dimension of the input matrix
0020 %   ERPU -   N_dim*1 array containing end of row pointer of all rows except
0021 %            the last row which doesn have any off diagonal element, N_dim is
0022 %            the dimension of the input matrix A in the original Ax = b
0023 %            problem
0024 %   ERP - (N_node+1)*1 array containning the end of row pointer data
0025 %   stop - scalar, equal to the number of external buses
0026 %   ERPEQ, CIndxEQ, 1*n arrays, together build the pointers of the
0027 %       equivalent branches
0028 %   BoundBus - 1*n array, list of boundary buses
0029 %
0030 % OUTPUT DATA:
0031 %   DataEQ: 1*n array, includes reactance value of the equivalent branches
0032 %   DataShunt: 1*n array, includes equivalent bus shunts of the boundary
0033 %   buses
0034 
0035 %   MATPOWER
0036 %   Copyright (c) 2014-2016, Power Systems Engineering Research Center (PSERC)
0037 %   by Yujia Zhu, PSERC ASU
0038 %
0039 %   This file is part of MATPOWER/mx-reduction.
0040 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0041 %   See https://github.com/MATPOWER/mx-reduction/ for more info.
0042 
0043 % do numerical factorization on input data
0044 numRow = dim; % number of rows of given data matrix
0045 ICPL = ERPU+1; % the initial column pointer equal to the last end of row pointer+1
0046 ICPL(end) = [];
0047 ICPL = [0,ICPL]; 
0048 RX = 0; % Initiate the RX value;
0049 Link = zeros(1,dim);
0050 ExAcum = Link; % Initiate the ExAcum;
0051 Diag = zeros(1,numRow); 
0052 DataEQ=zeros(size(CIndxEQ));
0053 DataShunt=zeros(1,length(BoundBus));
0054 % Initialization based on ERPU, CindU;
0055 % Sort the CIndxU to make it Ordered CindUU->CindUO
0056 for i = 2:numRow
0057     RowColInd = ERPU(i-1)+1:ERPU(i); % for every row the pointer of the column index
0058     [CIndxU(RowColInd)] = sort(CIndxU(RowColInd)); % sort the CIndx of every row in ascending order
0059 end
0060 %% begin Numerical Factorization
0061 RowIndex = 1; % Start at row 1
0062 while RowIndex<=numRow
0063     %% step 1 a,b
0064     % zero ExAcum using Link and CIndxUO;
0065     % This give the active element of current row
0066     % get the array from the self referential link
0067      if RowIndex>stop
0068             ExAcumEQ=zeros(size(ExAcum));
0069         end
0070     if Link(RowIndex)~=0
0071         [LinkPos,LinkArray,LinkCounter] = SelfLink(Link,RowIndex);
0072     else LinkCounter=0;
0073     end
0074     if LinkCounter~=0
0075         if RowIndex<numRow % if this is the last row there will be nothing on the right of the diagonal in U only fills in row(numRow) of L
0076             ExAcum([LinkArray,CIndxU(ERPU(RowIndex)+1:ERPU(RowIndex+1))])=0;% zero non-zero position from both native and fill
0077         else
0078             ExAcum(LinkArray)=0; % for last row, fill only
0079         end
0080     else
0081         ExAcum(CIndxU(ERPU(RowIndex)+1:ERPU(RowIndex+1)))=0;
0082     end
0083     %% step 1c
0084     % load corresponding values to ExAcum
0085     ExAcum(CIndx(ERP(RowIndex)+1:ERP(RowIndex+1)))=Data(ERP(RowIndex)+1:ERP(RowIndex+1)); % Index in the original array is CIndx(ERP(RowIndex)+1:ERP(RowIndex+1))
0086     %% step 2a
0087     RX = 0; % initiate RX
0088     %% step 2b,c
0089     if LinkCounter~=0
0090         [LinkArray,LinkPosInd]=sort(LinkArray);
0091         LinkPos = LinkPos(LinkPosInd);
0092         Link(LinkPos)=0;
0093        
0094         for i = 1:LinkCounter
0095             RX = LinkArray(i); % assign RX value to current fill generating row
0096             %% step 2d
0097            if RowIndex>stop
0098             ExAcumEQ(CIndxU(ERPU(RX)+1:ERPU(RX+1)))=ExAcumEQ(CIndxU(ERPU(RX)+1:ERPU(RX+1)))-ExAcum(RX)*URO(ERPU(RX)+1:ERPU(RX+1));
0099            end
0100             ExAcum(CIndxU(ERPU(RX)+1:ERPU(RX+1)))=ExAcum(CIndxU(ERPU(RX)+1:ERPU(RX+1)))-ExAcum(RX)*URO(ERPU(RX)+1:ERPU(RX+1));
0101             
0102             %% step 2ef
0103             LCO(ICPL(RX+1))=ExAcum(RX)*Diag(RX);
0104             ICPL(RX+1) = ICPL(RX+1)+1;
0105 %             Link(LinkPos(i))=0;
0106             %% step 2g
0107             if ICPL(RX+1)<=ERPU(RX+1)
0108                 SelfRef=CIndxU(ICPL(RX+1));
0109                 while Link(SelfRef)~=0
0110                     SelfRef = Link(SelfRef); % exhaust the link list and find a 0 position to store RX
0111                 end
0112                 Link(SelfRef)=RX;
0113             end
0114         end
0115     end
0116     if RowIndex>stop
0117     DataEQ(ERPEQ(RowIndex)+1:ERPEQ(RowIndex+1))=1./ExAcumEQ(CIndxEQ(ERPEQ(RowIndex)+1:ERPEQ(RowIndex+1)));
0118     DataShunt(RowIndex-stop)=ExAcumEQ(RowIndex);%+sum(ExAcumEQ(CIndxEQ(ERPEQ(RowIndex)+1:ERPEQ(RowIndex+1))));
0119     end
0120     %% step 4
0121 if RowIndex<=stop
0122     Diag(RowIndex)=1/ExAcum(RowIndex); % Invert the diagonal value
0123     %% step 5
0124     if RowIndex<numRow
0125         URO(ERPU(RowIndex)+1:ERPU(RowIndex+1))=ExAcum(CIndxU(ERPU(RowIndex)+1:ERPU(RowIndex+1)))*Diag(RowIndex); % Multiply active ExAcum by Diag(1) & store in URO
0126         %% step 6
0127         SelfRef=CIndxU(ICPL(RowIndex+1));
0128         while Link(SelfRef)~=0
0129             SelfRef=Link(SelfRef);
0130         end
0131         Link(SelfRef)=RowIndex;
0132     end
0133    
0134 elseif sum(Link)==0
0135     break;
0136 end
0137  %% Prepare for next loop
0138 RowIndex = RowIndex+1; % Increment the RowIndex
0139 end
0140 end

Generated on Fri 09-Oct-2020 11:21:31 by m2html © 2005