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-2015 by Power System Engineering Research Center (PSERC) 0037 % by Yujia Zhu, PSERC ASU 0038 % 0039 % $Id: PartialNumLU.m 2655 2015-03-18 16:40:32Z ray $ 0040 % 0041 % This file is part of MATPOWER. 0042 % Covered by the 3-clause BSD License (see LICENSE file for details). 0043 % See http://www.pserc.cornell.edu/matpower/ for more info. 0044 0045 % do numerical factorization on input data 0046 numRow = dim; % number of rows of given data matrix 0047 ICPL = ERPU+1; % the initial column pointer equal to the last end of row pointer+1 0048 ICPL(end) = []; 0049 ICPL = [0,ICPL]; 0050 RX = 0; % Initiate the RX value; 0051 Link = zeros(1,dim); 0052 ExAcum = Link; % Initiate the ExAcum; 0053 Diag = zeros(1,numRow); 0054 DataEQ=zeros(size(CIndxEQ)); 0055 DataShunt=zeros(1,length(BoundBus)); 0056 % Initialization based on ERPU, CindU; 0057 % Sort the CIndxU to make it Ordered CindUU->CindUO 0058 for i = 2:numRow 0059 RowColInd = ERPU(i-1)+1:ERPU(i); % for every row the pointer of the column index 0060 [CIndxU(RowColInd)] = sort(CIndxU(RowColInd)); % sort the CIndx of every row in ascending order 0061 end 0062 %% begin Numerical Factorization 0063 RowIndex = 1; % Start at row 1 0064 while RowIndex<=numRow 0065 %% step 1 a,b 0066 % zero ExAcum using Link and CIndxUO; 0067 % This give the active element of current row 0068 % get the array from the self referential link 0069 if RowIndex>stop 0070 ExAcumEQ=zeros(size(ExAcum)); 0071 end 0072 if Link(RowIndex)~=0 0073 [LinkPos,LinkArray,LinkCounter] = SelfLink(Link,RowIndex); 0074 else LinkCounter=0; 0075 end 0076 if LinkCounter~=0 0077 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 0078 ExAcum([LinkArray,CIndxU(ERPU(RowIndex)+1:ERPU(RowIndex+1))])=0;% zero non-zero position from both native and fill 0079 else 0080 ExAcum(LinkArray)=0; % for last row, fill only 0081 end 0082 else 0083 ExAcum(CIndxU(ERPU(RowIndex)+1:ERPU(RowIndex+1)))=0; 0084 end 0085 %% step 1c 0086 % load corresponding values to ExAcum 0087 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)) 0088 %% step 2a 0089 RX = 0; % initiate RX 0090 %% step 2b,c 0091 if LinkCounter~=0 0092 [LinkArray,LinkPosInd]=sort(LinkArray); 0093 LinkPos = LinkPos(LinkPosInd); 0094 Link(LinkPos)=0; 0095 0096 for i = 1:LinkCounter 0097 RX = LinkArray(i); % assign RX value to current fill generating row 0098 %% step 2d 0099 if RowIndex>stop 0100 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)); 0101 end 0102 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)); 0103 0104 %% step 2ef 0105 LCO(ICPL(RX+1))=ExAcum(RX)*Diag(RX); 0106 ICPL(RX+1) = ICPL(RX+1)+1; 0107 % Link(LinkPos(i))=0; 0108 %% step 2g 0109 if ICPL(RX+1)<=ERPU(RX+1) 0110 SelfRef=CIndxU(ICPL(RX+1)); 0111 while Link(SelfRef)~=0 0112 SelfRef = Link(SelfRef); % exhaust the link list and find a 0 position to store RX 0113 end 0114 Link(SelfRef)=RX; 0115 end 0116 end 0117 end 0118 if RowIndex>stop 0119 DataEQ(ERPEQ(RowIndex)+1:ERPEQ(RowIndex+1))=1./ExAcumEQ(CIndxEQ(ERPEQ(RowIndex)+1:ERPEQ(RowIndex+1))); 0120 DataShunt(RowIndex-stop)=ExAcumEQ(RowIndex);%+sum(ExAcumEQ(CIndxEQ(ERPEQ(RowIndex)+1:ERPEQ(RowIndex+1)))); 0121 end 0122 %% step 4 0123 if RowIndex<=stop 0124 Diag(RowIndex)=1/ExAcum(RowIndex); % Invert the diagonal value 0125 %% step 5 0126 if RowIndex<numRow 0127 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 0128 %% step 6 0129 SelfRef=CIndxU(ICPL(RowIndex+1)); 0130 while Link(SelfRef)~=0 0131 SelfRef=Link(SelfRef); 0132 end 0133 Link(SelfRef)=RowIndex; 0134 end 0135 0136 elseif sum(Link)==0 0137 break; 0138 end 0139 %% Prepare for next loop 0140 RowIndex = RowIndex+1; % Increment the RowIndex 0141 end 0142 end