MATPOWER / matpower

MATPOWER – steady state power flow simulation and optimization for MATLAB and Octave
https://matpower.org
Other
423 stars 151 forks source link

an improvement suggestion for makeYbus.m #52

Closed AsprinChina closed 5 years ago

AsprinChina commented 5 years ago

In makeYbus.m, we use the multiplication of the connection matrices and branch admittance matrices to get the node admittance matrices, like this:

Ybus = Cf' Yf + Ct' Yt

 In fact, since we have get the lists of "from" buses and "to" buses, we could from the node admittance matrices directly without multiplication, like this

Ybus = sparse([f,f,t,t], [f,t,f,t], [Yff,Yft,Ytf,Ytt], nb, nb)

 If we use the latter to replace the former, the time to form the node admittance matrix could reduce 70%. And the performance of this function could have of lift of 18%, which has been tested  on the forming of the case39.

mpc = case39(); for i=1:100000 makeYbus(mpc); end for i=1:100000 makeYbus_modified(mpc); end

The time consumed are shown in this figure. time

The makeYbus_modified.m is `function [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch) %MAKEYBUS Builds the bus admittance matrix and branch admittance matrices. % [YBUS, YF, YT] = MAKEYBUS(MPC) % [YBUS, YF, YT] = MAKEYBUS(BASEMVA, BUS, BRANCH) %
% Returns the full bus admittance matrix (i.e. for all buses) and the % matrices YF and YT which, when multiplied by a complex voltage vector, % yield the vector currents injected into each line from the "from" and % "to" buses respectively of each line. Does appropriate conversions to p.u. % Inputs can be a MATPOWER case struct or individual BASEMVA, BUS and % BRANCH values. Bus numbers must be consecutive beginning at 1 % (i.e. internal ordering). % % See also MAKEJAC, MAKESBUS, EXT2INT.

% MATPOWER % Copyright (c) 1996-2016, Power Systems Engineering Research Center (PSERC) % by Ray Zimmerman, PSERC Cornell % % This file is part of MATPOWER. % Covered by the 3-clause BSD License (see LICENSE file for details). % See http://www.pserc.cornell.edu/matpower/ for more info.

%% extract from MPC if necessary if nargin < 3 mpc = baseMVA; baseMVA = mpc.baseMVA; bus = mpc.bus; branch = mpc.branch; end

%% constants nb = size(bus, 1); %% number of buses nl = size(branch, 1); %% number of lines

%% define named indices into bus, branch matrices [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;

%% check that bus numbers are equal to indices to bus (one set of bus numbers) if any(bus(:, BUS_I) ~= (1:nb)') error('makeYbus: buses must be numbered consecutively in bus matrix; use ext2int() to convert to internal ordering') end

%% for each branch, compute the elements of the branch admittance matrix where %% %% | If | | Yff Yft | | Vf | %% | | = | | | | %% | It | | Ytf Ytt | | Vt | %% stat = branch(:, BR_STATUS); %% ones at in-service branches Ys = stat ./ (branch(:, BR_R) + 1j branch(:, BR_X)); %% series admittance Bc = stat . branch(:, BR_B); %% line charging susceptance tap = ones(nl, 1); %% default tap ratio = 1 i = find(branch(:, TAP)); %% indices of non-zero tap ratios tap(i) = branch(i, TAP); %% assign non-zero tap ratios tap = tap . exp(1jpi/180 branch(:, SHIFT)); %% add phase shifters Ytt = Ys + 1jBc/2;
Yff = Ytt ./ (tap .
conj(tap)); %
Yft = - Ys ./ conj(tap); Ytf = - Ys ./ tap;

%% compute shunt admittance %% if Psh is the real power consumed by the shunt at V = 1.0 p.u. %% and Qsh is the reactive power injected by the shunt at V = 1.0 p.u. %% then Psh - j Qsh = V conj(Ysh V) = conj(Ysh) = Gs - j Bs, %% i.e. Ysh = Psh + j Qsh, so ... Ysh = (bus(:, GS) + 1j * bus(:, BS)) / baseMVA; %% vector of shunt admittances

%% build connection matrices f = branch(:, F_BUS); %% list of "from" buses t = branch(:, T_BUS); %% list of "to" buses

%% build Ybus i = [1:nl; 1:nl]';
Yf = sparse(i, [f; t], [Yff; Yft], nl, nb); Yt = sparse(i, [f; t], [Ytf; Ytt], nl, nb); Ybus = sparse([f,f,t,t], [f,t,f,t], [Yff,Yft,Ytf,Ytt], nb, nb) + ... %% branch admittances sparse(1:nb, 1:nb, Ysh, nb, nb); %% shunt admittance ` Hope you could adopt this suggestion, and thanks for providing such a powerful instrument!

rdzman commented 5 years ago

Thanks for this suggestion.

Unfortunately, it seems the improved performance is not so consistent. In particular, for me it looks like the improvement only applies to very small systems and the change actually slows things down for larger systems, but only on MATLAB, not on Octave.

Running time_makeYbus.m, which also requires makeYbus_new.m, results in output like the following on my machine.

Note: You'll need to remove the .txt extension from each after downloading the two files linked above.

On MATLAB:

>> time_makeYbus
               case9 : 10000 reps : speedup = 24.6%
              case39 :  2329 reps : speedup = 7.0%
              case57 :  1603 reps : speedup = 6.3%
             case118 :   789 reps : speedup = 1.2%
             case300 :   327 reps : speedup = -1.5%
          case2383wp :    66 reps : speedup = -7.9%
         case6515rte :    42 reps : speedup = -17.7%
     case13659pegase :    35 reps : speedup = -10.6%
     case_ACTIVSg70k :    30 reps : speedup = -13.6%

On Octave:

octave:1> time_makeYbus
               case9 : 10000 reps : speedup = 4.6%
              case39 :  2329 reps : speedup = 5.9%
              case57 :  1603 reps : speedup = 7.1%
             case118 :   789 reps : speedup = 11.2%
             case300 :   327 reps : speedup = 16.2%
          case2383wp :    66 reps : speedup = 22.5%
         case6515rte :    42 reps : speedup = 9.6%
     case13659pegase :    35 reps : speedup = 11.0%
     case_ACTIVSg70k :    30 reps : speedup = 12.4%

Just out of curiosity, could you post the output of a typical run on your machine?

AsprinChina commented 5 years ago

Thanks for your reply!

This is the output of _timemakeYbus.m on my machine.

              case39 :  2329 reps : speedup = 7.2%
              case57 :  1603 reps : speedup = 6.3%
             case118 :   789 reps : speedup = -1.4%
             case300 :   327 reps : speedup = -2.9%
          case2383wp :    66 reps : speedup = -11.2%
         case6515rte :    42 reps : speedup = -13.8%
     case13659pegase :    35 reps : speedup = -14.1%

Exactly, this method only works effectively for small system. For a large system, it has a bad performance. Because in MATLAB, the sparse function occupies more time than matrix multiplication when the matrix is big enough, which I didn't notice before. Sorry to disturb you.

B.T.W. Your codes sets me a good example how to report a suggestion. Thanks for that again.

rdzman commented 5 years ago

Try running time_makeYbus again using this new version of makeYbus_new.m. It selects which method to use based on the size and whether it's running under MATLAB or Octave. For me it now results in a consistent performance improvement.

Here are typical results on my machine for MATLAB ...

>> time_makeYbus
               case9 : 10000 reps : speedup = 11.3%
              case39 :  2329 reps : speedup = 7.0%
              case57 :  1603 reps : speedup = 6.4%
             case118 :   789 reps : speedup = 1.3%
             case300 :   327 reps : speedup = 3.0%
          case2383wp :    66 reps : speedup = 19.1%
         case6515rte :    42 reps : speedup = 29.7%
     case13659pegase :    35 reps : speedup = 17.7%
     case_ACTIVSg70k :    30 reps : speedup = 8.4%

... and Octave ...

octave:1> time_makeYbus
               case9 : 10000 reps : speedup = 4.7%
              case39 :  2329 reps : speedup = 6.2%
              case57 :  1603 reps : speedup = 7.8%
             case118 :   789 reps : speedup = 12.7%
             case300 :   327 reps : speedup = 7.2%
          case2383wp :    66 reps : speedup = 26.3%
         case6515rte :    42 reps : speedup = 9.8%
     case13659pegase :    35 reps : speedup = 18.2%
     case_ACTIVSg70k :    30 reps : speedup = 20.0%

If you also see consistent improvements on your system, I think I'll go ahead and incorporate the change. If so, what name shall I use in the change log to thank you for the suggestion.

AsprinChina commented 5 years ago

After the method selecting, this is the result running on my machine.

>> time_makeYbus
               case9 : 10000 reps : speedup = 16.2%
              case39 :  2329 reps : speedup = 5.7%
              case57 :  1603 reps : speedup = 4.5%
             case118 :   789 reps : speedup = -1.8%
             case300 :   327 reps : speedup = 10.0%
          case2383wp :    66 reps : speedup = 24.0%
         case6515rte :    42 reps : speedup = 37.2%
     case13659pegase :    35 reps : speedup = 42.7%

...and if I change the threshold in line 78 of _makeYbusnew.m from if nb < 300 || have_fcn('octave') to if nb < 280 || have_fcn('octave'), then the performance improvement is more consistent:

>> time_makeYbus
               case9 : 10000 reps : speedup = 19.6%
              case39 :  2329 reps : speedup = 12.1%
              case57 :  1603 reps : speedup = 7.4%
             case118 :   789 reps : speedup = 0.0%
             case300 :   327 reps : speedup = 10.5%
          case2383wp :    66 reps : speedup = 23.4%
         case6515rte :    42 reps : speedup = 37.9%
     case13659pegase :    35 reps : speedup = 43.1%

Since there is MATLAB only in my machine, I cannot provide the test result of Octave on my machine. But I guess the result is similar.

I am Binbin Chen, a student in Tsinghua University majored in Electrical Engineering, and it's my glad to provide a suggestion for MatPower.