fm_pmuloc.m

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

    for i = 1:length(pmu_test2)
      numpmu = min(numpmu, length(pmu_test2{i,1}));
    end
    for i = 1:length(pmu_test2),
      if numpmu == length(pmu_test2{i,1})
        index_pmu = [index_pmu; i];
      end
    end
    pmu_test = pmu_test2(index_pmu);
    I_idx = I_idx_test(index_pmu);
    pseudoi = pseudi_test(index_pmu);
    if Fig.pmu
      set(hdl_pmu,'String',int2str(pmunum(index_pmu(1))));
      set(hdl_nob,'String',int2str(0));
      drawnow
    end
    pmunum = numpmu;
    Line.con = Line_old;
    Bus.con(:,1) = bus_old;
    Bus.int = busintold;
    Line.Y = Y_old;

   case 7

    metodo = 'Direct (N-1) Spanning Tree';
    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);
      coll_idx = find(coll == j & pesi < 100);
      for i = 1:length(coll_idx)
        pmu_idx = find(A(coll_idx(i),:));
        a = find(pesi(pmu_idx) >= 1);
        b = find(pesi(pmu_idx) > 99);
        if length(a) >= 1 & length(b) < round(2*j/3)
          %isempty(find(pesi(pmu_idx(a)) > 99))
          pmu_test = [pmu_test, Bus.con(coll_idx(i),1)];
          pesi(pmu_idx) = pesi(pmu_idx) + 1;
          pesi(coll_idx(i)) = 100;
        end
      end
      j = j + 1;
      if j > max_coll, break, end
    end

    % assenazione di un PMU ai nodi non raggiungibili due volte
    while 1
      a = find(pesi < 99);
      b = coll(a);
      [b, c] = sort(b);
      a = a(c);
      for i = length(a):-1:1
        coll_idx = find(A(a(i),:));
        stoppa = 1;
        if length(find(pesi(coll_idx) > 99)) < 2 & length(coll_idx) > 1
          pesi(a(i)) = pesi(a(i))+100;
          pmu_test = [pmu_test, Bus.con(a(i),1)];
          stoppa = 0;
          break
        end
      end
      if stoppa, break, end
    end

    % determinazione della matrice di incidenza delle correnti misurate
    ntest = length(pmu_test);
    if Fig.pmu
      set(hdl_pmu,'String',int2str(ntest));
      set(hdl_nob,'String',int2str(length(find(pesi == 0))));
      drawnow
    end
    linee = Line.con(:,[1 2]);
    I_idx = [];
    nodi = [];
    pseudoi = 0;

    for ijk = 1:ntest
      i_idx = fm_iidx(pmu_test(ijk),linee);
      I_idx = [I_idx; i_idx];
      nodi_oss = [i_idx(:,4);pmu_test(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

    index_pmu = 1;

  end

  if iscell(pmu_test)
    I_idx{1}(:,1) = [1:length(I_idx{1}(:,1))]';
    [Vest, angest, rangoH] = fm_lssest(pmu_test{1}, I_idx{1},1);
  else
    I_idx(:,1) = [1:length(I_idx(:,1))]';
    [Vest, angest, rangoH] = fm_lssest(pmu_test, I_idx,1);
  end

  pmucell = cell(Bus.n,1);
  vcell = cell(Bus.n,1);
  angcell = cell(Bus.n,1);

  for i = 1:Bus.n
    angolo = round(angest(i)/Settings.lftol)*Settings.lftol;
    if angolo < 0
      segno = '-';
    else
      segno = ' ';
    end
    angcell{i,1} = [segno, num2str(abs(angolo))];
  end

  if ~iscell(pmu_test)
    pmunum = length(pmu_test);
  end
  pmucif = ceil(log10(pmunum))+2;
  b = '                             ';
  b = b(1:pmucif);

  if type < 4, index_pmu = 1; end
  Matrix = cell(3,1);
  if length(index_pmu) > 1
    Matrix{3,1} = zeros(Bus.n+2,length(index_pmu));
  else
    Matrix{3,1} = zeros(Bus.n,1);
  end

  for i = 1:Bus.n
    vcell{i,1} = num2str(Vest(i));
    if iscell(pmu_test)
      pmuyes = find(pmu_test{1} == Bus.con(i,1));
    else
      pmuyes = find(pmu_test == Bus.con(i,1));
    end
    pmuname = Varname.bus{Bus.int(Bus.con(i,1))};
    if ~isempty(pmuyes)
      a = [int2str(pmuyes(1)),'      '];
      a = a(1:pmucif);
      vcell{i,1} = ['PMU ',a,' # ', vcell{i,1}];
      angcell{i,1} = ['PMU ',a,' # ', angcell{i,1}];
      pmucell{i,1} = ['PMU ',a,' # ',pmuname];
      Matrix{3,1}(i,1) = 1;
    else
      vcell{i,1} = [b,'       ', vcell{i,1}];
      angcell{i,1} = [b,'       ', angcell{i,1}];
      pmucell{i,1} = [b,'       ',pmuname];
    end
  end

  if type >= 4 & length(index_pmu) > 1
    for k = 2:length(index_pmu)
      pmu_test = pmu_test2{index_pmu(k),1};
      for i = 1:Bus.n
        pmuyes = find(pmu_test == Bus.con(i,1));
        pmuname = Varname.bus{Bus.int(Bus.con(i,1))};
        if ~isempty(pmuyes)
          a = [int2str(pmuyes(1)),'      '];
          a = a(1:pmucif);
          pmucell{i,k} = ['PMU ',a,' # ',pmuname];
          Matrix{3,1}(i,k) = 1;
        else
          pmucell{i,k} = [b,'       ',pmuname];
        end
      end
    end
  end

  pmunum = pmunum(1);
  if iscell(pseudoi)
    pseudi = pseudoi{1};
  else
    pseudi = pseudoi;
  end
  if iscell(I_idx)
    measui = length(I_idx{1}(:,1)) - pseudoi{1};
  else
    measui = length(I_idx(:,1)) - pseudoi;
  end

  if Fig.pmu
    set(hdl_pmv,'String',pmucell,'Value',1);
    set(hdl_ang,'String',angcell,'Value',1);
    set(hdl_V,'String',vcell,'Value',1);
    set(hdl_pmc,'String',pseudi);
    set(hdl_mv,'String',num2str(pmunum));
    set(hdl_mc,'String',measui);
  end

  PMU.location = pmucell;
  PMU.angle = angcell;
  PMU.voltage = vcell;
  PMU.measv = pmunum;
  PMU.measc = measui;
  PMU.pseudo = pseudi;
  PMU.number = pmunum;
  PMU.noobs = 0;

  tempo = toc;
  secondi = rem(tempo,60);
  tempo = (tempo - secondi)/60;
  minuti = rem(tempo,60);
  ore = (tempo - minuti)/60;
  durata = [int2str(ore),'h   ',int2str(minuti),'m   ',num2str(secondi),'s'];


  % PMU report
  % ----------------------------------------------------------

  % Headings
  Header{1,1}{1,1} = 'PMU PLACEMENT REPORT';
  Header{1,1}{2,1} = ' ';
  Header{1,1}{3,1} = ['P S A T  ',Settings.version];
  Header{1,1}{4,1} = ' ';
  Header{1,1}{5,1} = 'Author:  Federico Milano, (c) 2002-2006';
  Header{1,1}{6,1} = 'e-mail:  [e-mail]';
  Header{1,1}{7,1} = 'website: http://thunderbox.uwaterloo.ca/~fmilano';
  Header{1,1}{8,1} = ' ';
  Header{1,1}{9,1} = ['File:  ', Path.data,strrep(File.data,'(mdl)','.mdl')];
  Header{1,1}{10,1} = ['Date:  ',datestr(now,0)];
  Header{1,1}{11,1} = ' ';
  Header{1,1}{12,1} = ['Placement Method:  ', metodo];
  Header{1,1}{13,1} = ['Elapsed Time:      ', durata];

  Matrix{1,1} = [];
  Cols{1,1} = '';
  Rows{1,1} = '';

  % Network and PMU statistics
  Header{2,1} = 'STATISTICS';
  Cols{2,1} = '';
  Rows{2,1}{1,1} = 'Buses';
  Matrix{2,1}(1,1) = Bus.n;
  Rows{2,1}{2,1} = 'Lines';
  Matrix{2,1}(2,1) = Line.n;
  Rows{2,1}{3,1} = 'PMUs';
  Matrix{2,1}(3,1) = PMU.number;
  Rows{2,1}{4,1} = 'PMU Sets';
  Matrix{2,1}(4,1) = length(pseudoi);
  if type < 4 | length(index_pmu) == 1
    Rows{2,1}{5,1} = 'Meas. Currents';
    Matrix{2,1}(5,1) = PMU.measc;
    Rows{2,1}{6,1} = 'Pseudo-Meas. Currents';
    Matrix{2,1}(6,1) = PMU.pseudo;
  end

  % PMU Placement
  nconf = length(index_pmu);
  Header{3,1} = 'PMU PLACEMENT';
  Cols{3,1}{1,1} = 'Bus Name';
  for kk = 1:nconf
    Cols{3,1}{1,kk+1} = ['Set ',num2str(kk)];
  end
  Rows{3,1} = Varname.bus;
  if nconf > 1
    Rows{3,1}{Bus.n+1,1} = 'MC';
    Rows{3,1}{Bus.n+2,1} = 'PMC';
    for k = 1:nconf
      Matrix{3,1}(Bus.n+1,k) = pseudoi{k,1};
      Matrix{3,1}(Bus.n+2,k) = length(I_idx{k}(:,1) - pseudoi{k,1});
    end
  end

  if nconf > 7
    uno = fix(nconf/7);
    due = rem(nconf,7);
    for i = 2:uno;
      Header{2+i,1} = 'PMU PLACEMENT';
      Rows{2+i,1} = Varname.bus;
      Rows{2+i,1}{Bus.n+1,1} = 'MC';
      Rows{2+i,1}{Bus.n+2,1} = 'PMC';
      Cols{2+i,1}{1,1} = 'Bus Name';
      idx1 = (i-1)*7+1;
      idx2 = i*7;
      for kk = 1:7
        Cols{2+i,1}{1,kk+1} = ['Set ',num2str(kk+idx1-1)];
      end
      Matrix{2+i,1} = Matrix{3,1}(:,idx1:idx2);
    end
    if due
      Header{3+uno,1} = 'PMU PLACEMENT';
      Rows{3+uno,1} = Varname.bus;
      Rows{3+uno,1}{Bus.n+1,1} = 'MC';
      Rows{3+uno,1}{Bus.n+2,1} = 'PMC';
      Cols{3+uno,1}{1,1} = 'Bus Name';
      idx1 = uno*7+1;
      idx2 = uno*7+due;
      for kk = 1:due
        Cols{3+uno,1}{1,kk+1} = ['Set ',num2str(kk+uno*7)];
      end
      Matrix{3+uno,1} = Matrix{3,1}(:,idx1:idx2);
    end
    Cols{3,1} = Cols{3,1}(1:8);
    Matrix{3,1} = Matrix{3,1}(:,1:7);
  end


  PMU.report.Matrix = Matrix;
  PMU.report.Header = Header;
  PMU.report.Cols = Cols;
  PMU.report.Rows = Rows;

  % End PMU report
  % -----------------------------------------------------------

  fm_disp(['PMU Placement Routine successfully completed in ',durata],1)
  if Settings.beep, beep, end