function Ybus = makeYbus(baseMVA, bus, branch)%MAKEYBUS   Builds the bus admittance matrix and branch admittance matrices.%   [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.%   MATPOWER%   $Id: makeYbus.m,v 1.8 2007/06/26 16:16:17 ray Exp $%   by Ray Zimmerman, PSERC Cornell%   Copyright (c) 1996-2005 by Power System Engineering Research Center (PSERC)%   See http://www.pserc.cornell.edu/matpower/ for more info.%% constantsj = sqrt(-1);nb = size(bus, 1);          %% number of busesnl = 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('buses must appear in order by bus number')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 branchesYs = stat ./ (branch(:, BR_R) + j * branch(:, BR_X));   %% series admittanceBc = stat .* branch(:, BR_B);                           %% line charging susceptancetap = ones(nl, 1);                              %% default tap ratio = 1i = find(branch(:, TAP));                       %% indices of non-zero tap ratiostap(i) = branch(i, TAP);                        %% assign non-zero tap ratiostap = tap .* exp(j*pi/180 * branch(:, SHIFT));  %% add phase shiftersYtt = Ys + j*Bc/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) + j * bus(:, BS)) / baseMVA;  %% vector of shunt admittances%% build Ybusf = branch(:, F_BUS);                           %% list of "from" busest = branch(:, T_BUS);                           %% list of "to" busesCf = sparse(f, 1:nl, ones(nl, 1), nb, nl);      %% connection matrix for line & from busesCt = sparse(t, 1:nl, ones(nl, 1), nb, nl);      %% connection matrix for line & to busesYbus = spdiags(Ysh, 0, nb, nb) + ...            %% shunt admittance    Cf * spdiags(Yff, 0, nl, nl) * Cf' + ...    %% Yff term of branch admittance    Cf * spdiags(Yft, 0, nl, nl) * Ct' + ...    %% Yft term of branch admittance    Ct * spdiags(Ytf, 0, nl, nl) * Cf' + ...    %% Ytf term of branch admittance    Ct * spdiags(Ytt, 0, nl, nl) * Ct';         %% Ytt term of branch admittance%% Build Yf and Yt such that Yf * V is the vector of complex branch currents injected%% at each branch's "from" bus, and Yt is the same for the "to" bus endreturn;