PLOT_MPC Plot an electrically meaningful drawing of a MATPOWER case. HANDLE = PLOT_MPC(MPC) HANDLE = PLOT_MPC(MPC, <option1_name>, <option1_value>, ...) Examples: plot_mpc('case118', 'MaxBusLabels', 40) Plots an electrically meaningful drawing of a MATPOWER case. Higher voltage branches will be drawn using thicker lines. Systems of more than ~300 buses may be very slow to draw. Inputs: MPC: a MATPOWER case struct or string with the name of the case file <options>: name, value pairs of options Three optional name/value attribute pairs are available: 'MaxBusLabels' / integer : This integer sets how many bus numbers to display. Labels are given to the highest voltage level buses first, and progressing down as space permits. 'DoThevLayout' / boolean : This sets which electrical distance metric to use as a basis for the layout. The default, 0, uses a distance metric based on power transfer distances, which tends to give attractive, legible layouts. Setting this to 1 will give a diagram which explicitly shows the effective Thevenin impedance between all bus pairs 'FigSize' / scalar : The basic image size is 8.89cm by 6.27 cm, which fits the column width of an IEEE journal. This can be scaled up by setting a scaler here, for larger or higher resolution images. The default is 3. This function implements visualization techniques described in this paper: P. Cuffe; A. Keane, "Visualizing the Electrical Structure of Power Systems," in IEEE Systems Journal, vol.PP, no.99, pp.1-12. DOI: 10.1109/JSYST.2015.2427994 http://dx.doi.org/10.1109/JSYST.2015.2427994 Figures produced using this script can be published freely, however a citation of the above paper would be appreciated. Contact: Paul Cuffe, University College Dublin, March 2016 paul.cuffe@ucd.ie
0001 function [handle] = plot_mpc(mpc, varargin) 0002 %PLOT_MPC Plot an electrically meaningful drawing of a MATPOWER case. 0003 % HANDLE = PLOT_MPC(MPC) 0004 % HANDLE = PLOT_MPC(MPC, <option1_name>, <option1_value>, ...) 0005 % 0006 % Examples: 0007 % plot_mpc('case118', 'MaxBusLabels', 40) 0008 % 0009 % Plots an electrically meaningful drawing of a MATPOWER case. 0010 % Higher voltage branches will be drawn using thicker lines. Systems of 0011 % more than ~300 buses may be very slow to draw. 0012 % 0013 % Inputs: 0014 % MPC: a MATPOWER case struct or string with the name of the case file 0015 % <options>: name, value pairs of options 0016 % Three optional name/value attribute pairs are available: 0017 % 0018 % 'MaxBusLabels' / integer : This integer sets how many bus numbers 0019 % to display. Labels are given to the highest voltage level 0020 % buses first, and progressing down as space permits. 0021 % 0022 % 'DoThevLayout' / boolean : This sets which electrical distance 0023 % metric to use as a basis for the layout. The default, 0, uses 0024 % a distance metric based on power transfer distances, which 0025 % tends to give attractive, legible layouts. Setting this to 1 0026 % will give a diagram which explicitly shows the effective 0027 % Thevenin impedance between all bus pairs 0028 % 0029 % 'FigSize' / scalar : The basic image size is 8.89cm by 6.27 cm, 0030 % which fits the column width of an IEEE journal. This can be 0031 % scaled up by setting a scaler here, for larger or higher 0032 % resolution images. The default is 3. 0033 % 0034 % This function implements visualization techniques described in this paper: 0035 % 0036 % P. Cuffe; A. Keane, "Visualizing the Electrical Structure of Power 0037 % Systems," in IEEE Systems Journal, vol.PP, no.99, pp.1-12. 0038 % DOI: 10.1109/JSYST.2015.2427994 0039 % http://dx.doi.org/10.1109/JSYST.2015.2427994 0040 % 0041 % Figures produced using this script can be published freely, however a 0042 % citation of the above paper would be appreciated. 0043 % 0044 % Contact: 0045 % Paul Cuffe, University College Dublin, March 2016 0046 % paul.cuffe@ucd.ie 0047 0048 % MATPOWER 0049 % Copyright (c) 2016, Power Systems Engineering Research Center (PSERC) 0050 % by Paul Cuffe 0051 % 0052 % This file is part of MATPOWER. 0053 % Covered by the 3-clause BSD License (see LICENSE file for details). 0054 % See http://www.pserc.cornell.edu/matpower/ for more info. 0055 0056 0057 MaxBusLabels=[]; %Maximum number of bus labels to plot - need to limit this to stop clutter 0058 DoThevLayout=[]; %Whether to use Thevenin impedance as the inter-node distance measure, or the power transfer distance measure 0059 FigSize=[]; %The multipier on figure size. Setting = 1 gives a figure with width 8.89 cm, which is equal to the column width for an IEEE publication. Height is set to width / sqrt(2) to maintain pleasing proportions 0060 0061 0062 i=1; 0063 while i <= length(varargin); 0064 switch (varargin{i}); 0065 case 'MaxBusLabels'; 0066 MaxBusLabels=varargin{i+1}; i=i+2; 0067 case 'DoThevLayout'; 0068 DoThevLayout = 1; i=i+2; 0069 case 'FigSize'; 0070 FigSize=varargin{i+1}; i=i+2; 0071 otherwise; 0072 error('Unknown option: %s\n',varargin{i}); i=i+1; 0073 end 0074 end 0075 0076 % If the variables remains empty after the above assignment give them 0077 % defaults 0078 if isempty(MaxBusLabels), MaxBusLabels = 50; end %50 seems a sensible default so things don't get too cluttered 0079 if isempty(DoThevLayout), DoThevLayout =0; end %0 means do power transfer layout 0080 if isempty(FigSize), FigSize=3; end 0081 0082 define_constants; 0083 0084 NumberedSystem = ext2int(runpf(mpc)); %Just in case bus numbers are all out of whack 0085 0086 NumBuses = size(NumberedSystem.bus,1); 0087 0088 if(NumBuses > 300) 0089 warning('Large system: layout calculation may be quite slow') 0090 end 0091 0092 if (DoThevLayout == 1) %Use internode effective impedances as the basis for the projection 0093 0094 [Ybuspu, Yf, Yt] = makeYbus(NumberedSystem); 0095 0096 Ybuspu = full(Ybuspu); 0097 Zbuspu = pinv(Ybuspu); 0098 0099 ThevMatrix = repmat(diag(Zbuspu),1, NumBuses) + repmat(diag(Zbuspu),1, NumBuses).' - 2*Zbuspu; %Klein resistance distance formula 0100 0101 ThevMatrix = ThevMatrix - tril(ThevMatrix); %Force symmetry in case a little numerical noise is preventing us being a strict distance matrix 0102 ThevMatrix = ThevMatrix + ThevMatrix.'; 0103 0104 [Layout, stressThev] = mdscale(abs(ThevMatrix),2,'criterion','sammon','start','random'); %apply multidimensional scaling to project into two dimensions 0105 0106 else %Use internode "power transfer distances" as the basis for the projection 0107 0108 TransMatrixMW = zeros(NumBuses); %init with zeros 0109 0110 for ThisBus = 1:NumBuses % loop through each bus 0111 0112 CurPTDFMatrix = makePTDF(NumberedSystem.baseMVA, NumberedSystem.bus, NumberedSystem.branch, ThisBus); %calculate the bus-to-branch PTDFs for an injection at every bus absorbed at ThisBus 0113 0114 TotalBusFlows = sum(abs(CurPTDFMatrix),1); %take ABS value as we're interested in total MW flow shift, regardless of direction 0115 0116 TransMatrixMW(ThisBus,:) = TotalBusFlows; 0117 0118 end 0119 0120 TransMatrixMW = TransMatrixMW - tril(TransMatrixMW); %Force symmetry in case a little numerical noise is preventing us being a strict distance matrix 0121 TransMatrixMW = TransMatrixMW + TransMatrixMW'; 0122 TransMatrixMW(isnan(TransMatrixMW)) = mean(TransMatrixMW(~isnan(TransMatrixMW))); %Just a hack in case some PTDF calculations haven't been succesful 0123 0124 [Layout, stressMW] = mdscale(TransMatrixMW,2,'criterion','sammon','start','random'); %apply multidimensional scaling to project into two dimensions 0125 0126 end 0127 0128 BranchVoltages = ( NumberedSystem.bus(NumberedSystem.branch(:, F_BUS), BASE_KV) + NumberedSystem.bus(NumberedSystem.branch(:, T_BUS), BASE_KV) ) / 2; %Divide by two so transformers take on intermediate values 0129 0130 DistinctVoltageLevels = unique([BranchVoltages; NumberedSystem.bus(:, BASE_KV)]); 0131 0132 handle = figure; 0133 0134 cmap = colormap(winter(size(DistinctVoltageLevels,1) + 5)); 0135 0136 AllHandles = []; 0137 0138 for ThisVoltageLevel = 1:size(DistinctVoltageLevels,1) %Run through and draw all the branches first 0139 0140 BranchesAtThisLevel = find(BranchVoltages == DistinctVoltageLevels(ThisVoltageLevel)); 0141 0142 if (~isempty(BranchesAtThisLevel)) 0143 0144 FromBuses = NumberedSystem.branch(BranchesAtThisLevel, F_BUS); 0145 ToBuses = NumberedSystem.branch(BranchesAtThisLevel, T_BUS); 0146 0147 LineSegments = [Layout(FromBuses,:), Layout(ToBuses,:)]; 0148 0149 CurrentThickness = ThisVoltageLevel/size(DistinctVoltageLevels,1) * 3; 0150 0151 ThisHandle = line([LineSegments(:,1),LineSegments(:,3)]' , [LineSegments(:,2), LineSegments(:,4)]', 'color', cmap(ThisVoltageLevel + 1, :), 'LineWidth', CurrentThickness*1.5, 'DisplayName', [int2str(DistinctVoltageLevels(ThisVoltageLevel)), ' kV'] ); 0152 0153 AllHandles = [AllHandles; ThisHandle(1)]; %Just take the first sample line from each 0154 0155 hold on; 0156 end 0157 0158 end 0159 0160 0161 for ThisVoltageLevel = 1:size(DistinctVoltageLevels,1) %Then draw in buses circles. Doing this in two steps stops unwanted overlaps 0162 0163 BusesAtThisLevel = find(NumberedSystem.bus(:, BASE_KV) == DistinctVoltageLevels(ThisVoltageLevel)); 0164 0165 CurrentThickness = ThisVoltageLevel/size(DistinctVoltageLevels,1) * 3; 0166 0167 hold on; 0168 0169 scatter(Layout(BusesAtThisLevel,1), Layout(BusesAtThisLevel,2), CurrentThickness * 12, cmap((ThisVoltageLevel + 1),:),'fill','MarkerEdgeColor','w'); 0170 0171 end 0172 0173 RemainingLabelSlots = MaxBusLabels; 0174 0175 for ThisVoltageLevel = size(DistinctVoltageLevels,1):-1:1 %Now start at the highest voltage level and work our way down, labelling bus nodes until we are too cluttered. 0176 0177 BusesAtThisLevel = find(NumberedSystem.bus(:, BASE_KV) == DistinctVoltageLevels(ThisVoltageLevel)); 0178 0179 if(size(BusesAtThisLevel,1) < RemainingLabelSlots) %For consistency, we should label all buses at a voltage level together, or not at all 0180 0181 text(Layout(BusesAtThisLevel,1), Layout(BusesAtThisLevel,2), int2str(NumberedSystem.order.bus.i2e(BusesAtThisLevel)), 'FontName', 'Times Roman', 'FontSize', 18, 'HorizontalAlignment','center','VerticalAlignment','top') 0182 0183 RemainingLabelSlots = RemainingLabelSlots - size(BusesAtThisLevel,1); 0184 0185 else 0186 break %Needs to be only the highest voltage buses we label 0187 end 0188 0189 end 0190 0191 FigWidth = 8.89 * FigSize; %Size of a column in an IEEE paper 0192 FigHeight = 6.27 * FigSize; %Height/sqrt(2) which is a classicaly nice proportion 0193 0194 set(gcf, 'Units','centimeters', 'Position',[0 0 FigWidth FigHeight]) 0195 set(gca, 'Units','centimeters') 0196 set(gca, 'Position',[0 0 FigWidth FigHeight]) 0197 0198 set(gcf, 'PaperUnits','centimeters', 'PaperPosition',[(FigWidth * 0.02) (FigHeight * 0.02) (FigWidth * 0.98) (FigHeight * 0.98)]) 0199 0200 axis([min(Layout(:,1))*1.05 max(Layout(:,1))*1.05 min(Layout(:,2))*1.1 max(Layout(:,2))*1.02]) %Just give a little breathing room around the edges 0201 0202 axis equal; %it's an electrical map so need consistent scale 0203 0204 set(gca,'YTick',get(gca,'XTick')) %consistent tick interval 0205 set(gca, 'XTickLabelMode', 'manual', 'XTickLabel', []); %no numbers just ticks 0206 set(gca, 'YTickLabelMode', 'manual', 'YTickLabel', []); %no numbers just ticks 0207 0208 0209 %Now put in a legend using the line handles we recorded earlier 0210 0211 NumHandles = size(AllHandles,1); %Just in case we have too many voltage levels 0212 0213 if(NumHandles > 5) %Limit it to five for clarity's sake 0214 0215 TrimHandles = [AllHandles(1); AllHandles(floor(NumHandles*0.25));AllHandles(floor(NumHandles*0.5));AllHandles(floor(NumHandles*0.75)) ;AllHandles(end)]; %Top and bottom and a few intermediate 0216 0217 else 0218 TrimHandles = AllHandles; 0219 end 0220 0221 VoltageLegend = legend(flipud(TrimHandles)); 0222 legend('boxoff'); 0223 0224 set(VoltageLegend, 'FontName', 'Times Roman', 'FontSize', 22); %To be consistent with node numbers 0225 set(VoltageLegend, 'Location', 'best'); %Locate it to not obscure data