fm_pmuloc.m

这是一个很适合研究和学习用的电力系统仿真软件

        if ~isempty(i_idx)
          I_idx = [I_idx; i_idx];
          % nodi osservabili dal primo PMU
          nodi_oss = [i_idx(:,4);pmu_test(ii)];
          nodi = [nodi; nodi_oss];
          linee(i_idx(:,2),[1 2]) = zeros(length(nodi_oss)-1,2);
        end

        % determinazione delle pseudo-correnti nelle linee ai cui
        % estremi sono note le tensioni
        pi_idx = [];
        for i = 1:length(nodi)
          I_idx_from = find(linee(:,1) == nodi(i));
          I_idx_to = [];
          for j = 1:length(nodi)
            ifrom = find(linee(I_idx_from,2) == nodi(j));
            I_idx_to = [I_idx_to; I_idx_from(ifrom)];
          end
          if ~isempty(I_idx_to);
            n_current = length(I_idx_to);
            api = [[1:n_current]', I_idx_to];
            bpi = linee(I_idx_to,[1 2]);
            cpi = ones(length(I_idx_to),1);
            pi_idx = [pi_idx; [api, bpi, cpi]];
          end
        end

        if ~isempty(pi_idx)
          % pmu_con(length(pmu_con));
          linee(pi_idx(:,2),[1 2]) = zeros(length(pi_idx(:,1)),2);
          I_idx = [I_idx; pi_idx];
          pseudoi = pseudoi + length(pi_idx(:,1));
        end

        % determinazione delle pseudo-correnti  determinate con la
        % legge di Kirchhoff per le correnti
        % ed eliminazione dei nodi di cui si può determinare la
        % tensione con la legge di Ohm
        count = 1;
        while count < length(nodi)
          %for ii = 1:length(nodi)
          if zeroinj(Bus.int(nodi(count))) == 0
            I_idx_from = find(linee(:,1) == nodi(count));
            I_idx_to = find(linee(:,2) == nodi(count));
            ncfrom = length(I_idx_from);
            ncto = length(I_idx_to);
            nc = ncfrom + ncto;
            if nc == 1
              if ncfrom == 1
                ki_idx = [1, I_idx_from, linee(I_idx_from,[1 2]), 1];
              else
                ki_idx = [1, I_idx_to, linee(I_idx_to,[2 1]), -1];
              end
              linee(ki_idx(2),[1 2]) = zeros(length(ki_idx(1)),2);
              I_idx = [I_idx; ki_idx];
              pseudoi = pseudoi + length(ki_idx(1));
              nodi_oss = [nodi_oss; ki_idx(4)];
              nodi = [ki_idx(4); nodi];
              count = 1;
            else
              count = count + 1;
            end
          else
            count = count + 1;
          end
        end

      end

      I_idx(:,1) = [1:length(I_idx(:,1))]';
      [Vest, angest, rangoH] = fm_lssest(pmu_test, I_idx,0);
      [nu, nl, pmu_test, I_idx, pseudoi, ok] = ...
          fm_annealing(rangoH, ntest, pmu_test, nu, nl, hdl_nob, zeroinj);
      if ok
        pmu_test_old = pmu_test;
        I_idx_old = I_idx;
        pseudoi_old = pseudoi;
        if Fig.pmu
          set(hdl_pmu,'String',int2str(nu));
          set(hdl_nob,'String',int2str(0));
          drawnow
        end
      end
    end

    if ~ok
      pmu_test = pmu_test_old;
      I_idx = I_idx_old;
      pseudoi = pseudoi_old;
    end
    if Fig.pmu, set(hdl_nob,'String',int2str(0)); end

   case 4

    metodo = 'Minimum Spanning Tree';
    fm_disp([metodo, ' Method '],1)
    [pmuloc, pmunum] = fm_mintree(zeroinj,hdl_pmu,hdl_nob);
    [pmu_test, pmu_test2, I_idx, pseudoi, index_pmu, pmunum] = ...
        fm_spantree(zeroinj, pmuloc, pmunum, hdl_pmu, hdl_nob);

   case 5

    metodo = 'Direct Spanning Tree';
    fm_disp([metodo, ' Method '],1)
    pesi = zeros(1,Bus.n);
    coll = zeros(1,Bus.n);

    % costruzione della matrice di adiacenza di rete
    A = zeros(Bus.n,Bus.n);
    for i = 1:Bus.n
      nonzero = find(Line.Y(i,:));
      A(i,nonzero) = ones(1,length(nonzero));
      A(i,i) = 0;
      coll(i) = sum(A(i,:));
    end

    % posizionamento dei PMU nei nodi collegati alle
    % antenne (interconnessione singola)
    pmu_test = [];
    coll1_idx = find(coll == 1);
    for i = 1:length(coll1_idx)
      pmu_idx = find(A(coll1_idx(i),:));
      pmu_test = [pmu_test, Bus.con(pmu_idx,1)];
      if Fig.pmu
        set(hdl_pmu,'String',int2str(length(pmu_test)));
        set(hdl_nob,'String',int2str(length(find(pesi == 0))));
        drawnow
      end
      pesi(pmu_idx) = 100;
      pmu_idx = find(A(pmu_idx,:));
      pesi(pmu_idx) = pesi(pmu_idx) + ones(1,length(pmu_idx));
    end

    % posizionamento dei PMU, nodi ad interconnessione maggiore di 1
    max_coll = max(coll);
    j = 2;
    while 1
      n_pmu = length(pmu_test);
      basta = 1;
      coll_idx = find(coll == j & pesi == 0);
      for i = 1:length(coll_idx)
        pmu_idx = find(A(coll_idx(i),:));
        a = find(pesi(pmu_idx) >= 1);
        b = find(pesi(pmu_idx) == 0);
        if length(b) == 1 & isempty(find(pesi(pmu_idx(a)) > 99))
          basta = 0;
          break
        end
      end
      if basta,
        j = j + 1;
        if j > max_coll, break, end
      else
        pesi(pmu_idx(a)) = pesi(pmu_idx(a)) + ones(1,length(pmu_idx(a)));
        pesi(pmu_idx(b)) = 100;
        pmu_test = [pmu_test, Bus.con(pmu_idx(b),1)];
        if Fig.pmu
          set(hdl_pmu,'String',int2str(length(pmu_test)));
          set(hdl_nob,'String',int2str(length(find(pesi == 0))));
          drawnow
        end
        pmu_idx2 = find(A(pmu_idx(b),:));
        pesi(pmu_idx2) = pesi(pmu_idx2) + ones(1,length(pmu_idx2));
      end
      if n_pmu < length(pmu_test), j = 2; end
    end

    % assegnazione di un PMU ai nodi non ancora osservabili
    while 1
      a = find(pesi == 0);
      if isempty(a), break, end
      pmu_test = [pmu_test, Bus.con(a(1),1)];
      pesi(a(1)) = 100;
      pmu_idx2 = find(A(a(1),:));
      pesi(pmu_idx2) = pesi(pmu_idx2) + ones(1,length(pmu_idx2));
      if Fig.pmu
        set(hdl_pmu,'String',int2str(length(pmu_test)));
        set(hdl_nob,'String',int2str(length(find(pesi == 0))));
        drawnow
      end
    end

    pmuloc = zeros(1,Bus.n);
    pmuloc(Bus.int(pmu_test)) = 1; %ones(1,pmunum);
    pmunum = length(find(pmuloc));
    [pmu_test, pmu_test2, I_idx, pseudoi, index_pmu, pmunum] = ...
        fm_spantree(zeroinj, pmuloc, pmunum, hdl_pmu, hdl_nob);

   case 6

    metodo = 'Minimum (N-1) Spanning Tree';
    fm_disp([metodo, ' Method '],1)

    Line_old = Line.con;
    bus_old = Bus.con(:,1);
    busintold = Bus.int;
    Y_old = Line.Y;
    %p = symrcm(abs(Line.Y));
    p = 1:Bus.n;
    Bus.con(:,1) = Bus.con(p,1);
    Line.con(:,[1 2]) = [Bus.con(p(Bus.int(Line.con(:,1))),1), ...
                        Bus.con(p(Bus.int(Line.con(:,2))),1)];
    Line.Y = Line.Y(p,p);
    busmax = max(Bus.con(:,1));
    Bus.int = zeros(busmax,1);
    for i = 1:Bus.n
      Bus.int(round(Bus.con(i,1))) = i;
    end

    bus20 = zeros(Bus.n,1);
    bus0 = bus20;

    for i = 1:Bus.n
      bus2 = bus20;
      bus2(i) = 1;
      pmu_test = Bus.con(i,1);
      bus_pmu = Bus.con(i,1);
      while ~isempty(find(~bus2))
        buspmu_new = [];
        for j = 1:length(bus_pmu)
          [buspmu, bus2]= fm_pmurec(bus_pmu(j),bus2,bus0);
          buspmu_new = [buspmu_new; buspmu];
          pmu_test = [pmu_test; buspmu];
        end
        bus_pmu = buspmu_new;
      end
      pmu_test = sort(pmu_test);
      num_bus = length(pmu_test);
      nodi_el = [];
      antenne = [];
      for jjj = 1:length(pmu_test)
        nodi_el = [nodi_el; jjj+find(pmu_test([jjj+1:length(pmu_test)]) ...
                                     == pmu_test(jjj))];
      end

      pmu_test(nodi_el) = [];
      for jjj = 1:length(pmu_test)
        a = length(fm_busidx(pmu_test(jjj)));
        if a == 1
          antenne = [antenne, jjj];
        end
      end
      pmu_test(antenne) = [];
      for jjj = 1:Bus.n
        if isempty(find(pmu_test == Bus.con(jjj,1)))
          a = fm_busidx(Bus.con(jjj,1));
          if length(a) > 1
            a_idx = [];
            for kkk = 1:length(a), a_idx = [a_idx; find(pmu_test == a(kkk))]; end
            if length(a_idx) < 2, pmu_test = [pmu_test; a]; end
          end
        end
      end
      nodi_el = [];
      for jjj = 1:length(pmu_test)
        nodi_el = [nodi_el; jjj+find(pmu_test([jjj+1:length(pmu_test)]) ...
                                     == pmu_test(jjj))];
      end

      pmu_test(nodi_el) = [];
      pmutest{i,1} = pmu_test;
      if Fig.pmu
        set(hdl_pmu,'String',int2str(length(pmu_test)));
        set(hdl_nob,'String',int2str(0));
        drawnow
      end
    end

    for i = 1:length(pmutest)

      pmu_test_new = pmutest{i};
      ntest = length(pmu_test_new);
      pmunum(i) = ntest;
      if Fig.pmu
        set(hdl_pmu,'String',int2str(ntest));
        drawnow
      end
      linee = Line.con(:,[1 2]);

      I_idx = [];
      nodi = [];
      pseudoi = 0;

      % nodi osservabili dai PMU posizionati con l'algoritmo
      % di spanning tree
      for ijk = 1:ntest
        i_idx = fm_iidx(pmu_test_new(ijk),linee);
        I_idx = [I_idx; i_idx];
        nodi_oss = [i_idx(:,4);pmu_test_new(ijk)];
        nodi = [nodi; nodi_oss];
      end

      nodi = sort(nodi);
      num_nodi = length(nodi);
      nodi_el = [];
      for jjj = 1:num_nodi
        nodi_el = [nodi_el; jjj+find(nodi([jjj+1:num_nodi]) == ...
                                     nodi(jjj))];
      end
      nodi(nodi_el) = [];

      % determinazione delle pseudo-correnti nelle linee ai cui
      % estremi sono note le tensioni
      pi_idx = [];
      for ii = 1:length(nodi)
        I_idx_from = find(linee(:,1) == nodi(ii));
        I_idx_to = [];
        for jj = 1:length(nodi)
          ifrom = find(linee(I_idx_from,2) == nodi(jj));
          I_idx_to = [I_idx_to; I_idx_from(ifrom)];
        end
        if ~isempty(I_idx_to);
          n_current = length(I_idx_to);
          api = [[1:n_current]', I_idx_to];
          bpi = linee(I_idx_to,[1 2]);
          cpi = ones(length(I_idx_to),1);
          pi_idx = [pi_idx; [api, bpi, cpi]];
        end
      end

      if ~isempty(pi_idx)
        linee(pi_idx(:,2),[1 2]) = zeros(length(pi_idx(:,1)),2);
        I_idx = [I_idx; pi_idx];
        pseudoi = pseudoi + length(pi_idx(:,1));
      end

      % determinazione delle pseudo-correnti  determinate con
      % la legge di Kirchhoff per le correnti
      % ed eliminazione dei nodi di cui si può determinare la
      % tensione con la legge di Ohm
      count = 1;
      while count < length(nodi)
        if zeroinj(Bus.int(nodi(count))) == 0
          I_idx_from = find(linee(:,1) == nodi(count));
          I_idx_to = find(linee(:,2) == nodi(count));
          ncfrom = length(I_idx_from);
          ncto = length(I_idx_to);
          nc = ncfrom + ncto;
          if nc == 1
            if ncfrom == 1
              ki_idx = [1, I_idx_from, linee(I_idx_from,[1 2]), 1];
            else
              ki_idx = [1, I_idx_to, linee(I_idx_to,[2 1]), -1];
            end
            linee(ki_idx(2),[1 2]) = zeros(length(ki_idx(1)),2);
            I_idx = [I_idx; ki_idx];
            pseudoi = pseudoi + length(ki_idx(1));
            nodi_oss = [nodi_oss; ki_idx(4)];
            nodi = [ki_idx(4); nodi];
            count = 1;
          else
            count = count + 1;
          end
        else
          count = count + 1;
        end
      end

      pmu_test2{i,1} = pmu_test_new;
      I_idx_test{i,1} = I_idx;
      pseudi_test{i,1} = pseudoi;
      if Fig.pmu
        set(hdl_nob,'String',int2str(0));
        drawnow
      end
    end
    numpmu = Bus.n;
    numpmuold = Bus.n;
    index_pmu = [];