Description of plot

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 %     https://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 Extras.
0053 %   Covered by the 3-clause BSD License (see LICENSE file for details).
0054 %   See https://github.com/MATPOWER/matpower-extras 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
plot_mpc

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
