LoadRedistribution

 Subroutine LoadRedistribution moves loads in reduced model in order to
 make the dcpf solution on reduced model identical to the full model with
 external generator placed in subroutine MoveExGen.

   [mpcreduced,BCIRCr]=LoadRedistribution(mpcfull,mpcreduced,BCIRCr,Pf_flag)

 INPUT DATA:
   mpcfull: struct, full model data in MATPOWER case format
   mpcreduced: struct, reduced model data with all external buses
       eliminated in MATPOWER case format
   BCIRCr: 1*n array, includes circuit number of branches in reduced model
   Pf_flag: scalar, indicate if dc power flow need to be solved before
   load redistribution.

 OUTPUT DATA:
   mpcreduced: struct, reduced model data with load redistributed
   BCIRCr: 1*n array, includes reordered branch circuit number in reduced
       model

 NOTE:
   The subroutine will first run a dc power flow on the full model. If the
   dc power flow can not be solved on the full model, the subroutine will
   be terminated and an error will be returned.

0001 function [mpcreduced,BCIRCr]=LoadRedistribution(mpcfull,mpcreduced,BCIRCr,Pf_flag)
0002 % Subroutine LoadRedistribution moves loads in reduced model in order to
0003 % make the dcpf solution on reduced model identical to the full model with
0004 % external generator placed in subroutine MoveExGen.
0005 %
0006 %   [mpcreduced,BCIRCr]=LoadRedistribution(mpcfull,mpcreduced,BCIRCr,Pf_flag)
0007 %
0008 % INPUT DATA:
0009 %   mpcfull: struct, full model data in MATPOWER case format
0010 %   mpcreduced: struct, reduced model data with all external buses
0011 %       eliminated in MATPOWER case format
0012 %   BCIRCr: 1*n array, includes circuit number of branches in reduced model
0013 %   Pf_flag: scalar, indicate if dc power flow need to be solved before
0014 %   load redistribution.
0015 %
0016 % OUTPUT DATA:
0017 %   mpcreduced: struct, reduced model data with load redistributed
0018 %   BCIRCr: 1*n array, includes reordered branch circuit number in reduced
0019 %       model
0020 %
0021 % NOTE:
0022 %   The subroutine will first run a dc power flow on the full model. If the
0023 %   dc power flow can not be solved on the full model, the subroutine will
0024 %   be terminated and an error will be returned.
0025 
0026 %   MATPOWER
0027 %   Copyright (c) 2014-2016, Power Systems Engineering Research Center (PSERC)
0028 %   by Yujia Zhu, PSERC ASU
0029 %
0030 %   This file is part of MATPOWER.
0031 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0032 %   See http://www.pserc.cornell.edu/matpower/ for more info.
0033 
0034 if Pf_flag==1
0035 % OPT=mpoption('out.all',0);
0036 [resultfull,successfull]=rundcpf(mpcfull);
0037 if ~successfull
0038     error('unable to solve dc powerflow with original full model, load cannot be redistributed')
0039 end
0040 else resultfull=mpcfull;
0041     successfull=1;
0042 end
0043     
0044 %% Read the full model bus data
0045 %  [BusID, V_mag, V_angle
0046 OrigBusRec=resultfull.bus(:,[1,8,9]);
0047 OrigBusRec=sortrows(OrigBusRec,1); %reorder bus records
0048 %% Read Bus Data
0049 BusRec = mpcreduced.bus;
0050 BusRec = sortrows(BusRec,1);
0051 BusNo = BusRec(:,1);
0052 
0053 %% Use original bus voltage
0054 [ignore, ind]=ismember(BusNo,OrigBusRec(:,1));
0055 BusRec(:,8)=OrigBusRec(ind,2); % Vm
0056 BusRec(:,9)=OrigBusRec(ind,3); % Vang
0057 %% Read Branch DATA
0058 branchdata = mpcreduced.branch;
0059 branchdata = branchdata(branchdata(:,11)==1,:);
0060 [branchdata,braindex] = sortrows(branchdata,[1,2]);
0061 BranchRec=branchdata;
0062 %% Renumber the branch terminal buses
0063 Sbase=100; %MVA
0064 NewBusNo = (1:length(BusNo))';
0065 BusRec(:,1) = NewBusNo;
0066 BranchRec(:,1) = interp1(BusNo,NewBusNo,BranchRec(:,1));
0067 BranchRec(:,2) = interp1(BusNo,NewBusNo,BranchRec(:,2));
0068 %% read phase shifter information
0069 ind=find( abs(BranchRec(:,10)) );
0070 flag=0;
0071 if isempty(ind)==0
0072     flag=1;
0073     phase_shifter=BranchRec(ind,:);
0074 end
0075 %% Form complex voltage vector
0076 Bus_V_Mag_PU=BusRec(:,8);
0077 Bus_V_Pha=BusRec(:,9)/180*pi;
0078 
0079 %% Form Y Matrix
0080 BB=zeros(length(BusNo),length(BusNo));
0081 bb=BranchRec(:,4);
0082 BranchRec(BranchRec(:,9)==0,9)=1;
0083 bb=bb.*(BranchRec(:,9)); % x/tap
0084 bb=1./bb;
0085 for i=1:length( BranchRec(:,4) )
0086 %     i
0087     m=BranchRec(i,1);
0088     n=BranchRec(i,2);
0089     
0090     BB(m,m)=BB(m,m)+bb(i);
0091     BB(n,n)=BB(n,n)+bb(i);
0092     
0093     BB(m,n)=BB(m,n)-bb(i);
0094     BB(n,m)=BB(n,m)-bb(i);  
0095 end
0096 
0097 P_injected2=BB*Bus_V_Pha*Sbase;
0098 
0099 if flag==1
0100     %% phase_shifter
0101     B_fix=zeros(length(BB(:,1)),1);
0102     for i=1:length( phase_shifter(:,1) )
0103         B_fix( phase_shifter(i,1) )= B_fix( phase_shifter(i,1) ) - phase_shifter(i,10)*pi/180/phase_shifter(i,4);
0104         B_fix( phase_shifter(i,2) )= B_fix( phase_shifter(i,2) ) + phase_shifter(i,10)*pi/180/phase_shifter(i,4);
0105     end
0106     B_fix=B_fix*Sbase;
0107 end
0108 
0109 gen = mpcreduced.gen;
0110 gen(:,2)=resultfull.gen(:,2); % use the full model solution
0111 gen(:,1) = interp1(BusNo,NewBusNo,gen(:,1));
0112 Generation = zeros(size(mpcreduced.bus,1),2);
0113 Generation(:,1) = NewBusNo;
0114 for i = 1:size(gen,1)
0115     Generation(gen(i,1),2) = Generation(gen(i,1),2)+gen(i,2);
0116 end
0117 gen(:,1) = interp1(NewBusNo,BusNo,gen(:,1));
0118 %% fix the phase shifter
0119 if flag==1
0120     P_injected2=P_injected2+B_fix;
0121 end
0122 P_L_should=Generation(:,2)-P_injected2;
0123 %% dealing with HVDC lines
0124 if isfield(mpcreduced,'dcline')
0125     dcline = mpcfull.dcline;
0126     HVDC_Line=[dcline(:,1),dcline(:,2),dcline(:,4),dcline(:,5)];  
0127     HVDC_Line=sortrows(HVDC_Line,[1 2]);    
0128     HVDC_Line(:,1)=interp1(BusNo,NewBusNo,HVDC_Line(:,1));
0129     HVDC_Line(:,2)=interp1(BusNo,NewBusNo,HVDC_Line(:,2));
0130     % for HVDC lines if one bus of a line is isolated then the buses on the other end
0131     % of the line will be ignored in the inverse power flow program
0132     for i=1:length(HVDC_Line(:,1))
0133         if (BusRec(HVDC_Line(i,1),2)~=4)&&(BusRec(HVDC_Line(i,2),2)~=4)
0134             P_L_should(HVDC_Line(i,1))=P_L_should(HVDC_Line(i,1))-HVDC_Line(i,3);
0135             P_L_should(HVDC_Line(i,2))=P_L_should(HVDC_Line(i,2))+HVDC_Line(i,4); % YZ compensate HVDC line by adding/reducing the loads from the HVDC flows
0136         end
0137     end
0138 end
0139 %% Plug in the results
0140 mpcreduced.bus(:,3)=P_L_should;
0141 mpcreduced.gen = gen;
0142 end