function [ref, pv, pq] = bustypes(bus, gen)%BUSTYPES   Builds lists of each type of bus (ref, pv, pq).%   [ref, pv, pq] = bustypes(bus, gen)%   Generators with "out-of-service" status are treated as PQ buses with%   zero generation (regardless of Pg/Qg values in gen).%   MATPOWER%   $Id: bustypes.m,v 1.4 2004/08/23 20:55:59 ray Exp $%   by Ray Zimmerman, PSERC Cornell%   Copyright (c) 1996-2004 by Power System Engineering Research Center (PSERC)%   See http://www.pserc.cornell.edu/matpower/ for more info.%% constants[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;[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...    GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;%% get generator status% bus_gen_status = zeros(size(bus, 1), 1);% bus_gen_status(gen(:, GEN_BUS)) = gen(:, GEN_STATUS) > 0;nb = size(bus, 1);ng = size(gen, 1);Cg = sparse(gen(:, GEN_BUS), [1:ng]', gen(:, GEN_STATUS) > 0, nb, ng);  %% gen connection matrix                                        %% element i, j is 1 if, generator j at bus i is ONbus_gen_status = Cg * ones(ng, 1);      %% number of generators at each bus that are ON%% form index lists for slack, PV, and PQ busesref = find(bus(:, BUS_TYPE) == REF & bus_gen_status);   %% reference bus indexpv  = find(bus(:, BUS_TYPE) == PV  & bus_gen_status);   %% PV bus indicespq  = find(bus(:, BUS_TYPE) == PQ | ~bus_gen_status);   %% PQ bus indices%% pick a new reference bus if for some reason there is none (may have been shut down)if isempty(ref)    ref = pv(1);                %% use the first PV bus    pv = pv(2:length(pv));      %% take it off PV listendreturn;