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
buses0001 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. 0040 % Covered by the 3-clause BSD License (see LICENSE file for details). 0041 % See http://www.pserc.cornell.edu/matpower/ 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