fm_ncomp.m

电力系统的psat

function check = fm_ncomp
% FM_NCOMP search components used in the current data
%          file and initializes fields used for power
%          flow computations.
%
% CHECK = FM_NCOMP
%
%see also FM_SLF
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    22-Feb-2003
%Version:   1.0.1
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2005 Federico Milano
%
% This toolbox is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2.0 of the License, or
% (at your option) any later version.
%
% This toolbox is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANDABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this toolbox; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
% USA.

global DAE Bus Settings clpsat

% component_con
global SW PV PQ Mn Pl Mot Svc Shunt Line Twt
global Syn Tap SAE1 SAE2 SAE3 Breaker Fl Varname
global Pss Ltc Oxl Exc Tg Fault Mass Thload

check = 0;

% check buses
if isempty(Bus.con)
  fm_disp(['The data file does not seem to be in a valid ', ...
           'format: no bus found.'])
  return
end
Bus.n = length(Bus.con(:,1));
busmax = max(Bus.con(:,1));
% set up internal bus numbers for second indexing of buses
Bus.int = zeros(busmax,1);
for i = 1:Bus.n
  Bus.int(round(Bus.con(i,1))) = i;
end

% check three-winding transformers
if ~isempty(Twt.con)

  % check data
  [n_twt_rows, n_twt_cols] = size(Twt.con);
  if n_twt_cols < 14
    fm_disp(['Three-winding transformer data does not ', ...
             'seems in a valid format',2])
    return
  elseif n_twt_cols < 15
    Twt.con = [Twt.con, ones(n_twt_rows,1), zeros(n_twt_rows,9)];
  elseif n_twt_cols < 24
    Twt.con = [Twt.con, zeros(n_twt_rows,24-n_twt_cols)];
  end
  
  % modify Bus.con
  twt_bus = busmax+[1:n_twt_rows]';
  twt_Vn = Twt.con(:,7);
  twt_ar = Bus.con(Bus.int(Twt.con(:,2)),[4 5]);
  Bus.con = [Bus.con; [twt_bus, twt_Vn, ...
                      ones(n_twt_rows,1), ...
                      zeros(n_twt_rows,1), twt_ar]];
  Bus.int = [Bus.int; Bus.n+[1:n_twt_rows]'];
  Bus.n = Bus.n + n_twt_rows;
  busmax = busmax + n_twt_rows;
  
  % check Varname
  if ~isempty(Varname.bus)
    for i = 1:n_twt_rows
      Varname.bus{end+1,1} = ...
          [Varname.bus{Bus.int(Twt.con(i,1))},'(twt)'];
    end
  end
  
  % compute branch impedances
  r12 = Twt.con(:,9);
  r13 = Twt.con(:,10);
  r23 = Twt.con(:,11);
  r1 = 0.5*(r12+r13-r23);
  r2 = 0.5*(r12+r23-r13);
  r3 = 0.5*(r23+r13-r12);
  x12 = Twt.con(:,12);
  x13 = Twt.con(:,13);
  x23 = Twt.con(:,14);
  x1 = 0.5*(x12+x13-x23);
  x2 = 0.5*(x12+x23-x13);
  x3 = 0.5*(x23+x13-x12);
  
  % check for zero impedances in the equivalent star
  idx = find(abs(x1) < 1e-5);
  if ~isempty(idx), x1(idx) = 0.0001; end
  idx = find(abs(x2) < 1e-5);
  if ~isempty(idx), x2(idx) = 0.0001; end
  idx = find(abs(x3) < 1e-5);
  if ~isempty(idx), x3(idx) = 0.0001; end

  % check Line.con column number
  if ~isempty(Line.con)
    ncol = length(Line.con(1,:));
    if ncol < 15,
      Line.con = [Line.con, zeros(Line.n,15-ncol)];
    end
    if ncol > 15,
      Line.con = Line.con(:,[1:15]);
    end
  end
  
  % modify Line.con
  line1 = [Twt.con(:,1),twt_bus, Twt.con(:,[4,6,5]),zeros(n_twt_rows,1), ...
           Twt.con(:,6)./Twt.con(:,7), r1,x1,zeros(n_twt_rows,1), ...
           Twt.con(:,15), zeros(n_twt_rows,1),Twt.con(:,[17,19,22])];
  line2 = [twt_bus,Twt.con(:,[2,4,7,5]),zeros(n_twt_rows,2), r2,x2, ...
           zeros(n_twt_rows,1), ones(n_twt_rows,1), ...
           zeros(n_twt_rows,1),Twt.con(:,[18,20,23])];
  line3 = [twt_bus,Twt.con(:,[3,4,7,5]),zeros(n_twt_rows,1), ...
           Twt.con(:,7)./Twt.con(:,8),r3,x3, zeros(n_twt_rows,1), ...
           ones(n_twt_rows,1), zeros(n_twt_rows,1),Twt.con(:,[19,21,24])];
  Line.con = [Line.con; line1; line2; line3];
  
  % clear Twt.con
  Twt.con = [];
end

% check lines
if ~isempty(Line.con)
  Line.n = length(Line.con(:,1));
end
% process line data
if Line.n > 0
  from_bus = Line.con(:,1);
  Line.from = Bus.int(round(from_bus));
  to_bus = Line.con(:,2);
  Line.to = Bus.int(round(to_bus));
  ncol = length(Line.con(1,:));
  if ncol < 15,
    Line.con = [Line.con, zeros(Line.n,15-ncol)];
  end
  Line.con(find(abs(Line.con(:,11))==0),11) = 1;
  % check for parameter units (if the length is > 0, unit is km)
  idx = find(Line.con(:,6));
  if ~isempty(idx)
    XB = Line.con(idx,4).*Line.con(idx,4)./Line.con(idx,3);
    Line.con(idx,8)  = Line.con(idx,6).*Line.con(idx,8)./XB;
    Line.con(idx,9)  = Line.con(idx,6).*Line.con(idx,9)./XB*Settings.rad;
    Line.con(idx,10) = Line.con(idx,6).*Line.con(idx,10).*XB/Settings.rad;
    bidx = find(Line.con(idx,10) > 10);
    if ~isempty(bidx)
      fm_disp('Warning: Some line susceptances are too high...')
      fm_disp(['         microFarad are assumed for those susceptance' ...
	       ' values'])
      Line.con(idx(bidx),10) = Line.con(idx(bidx),10)/1e6;
    end
  end
end

% set up shunt admittances
Shunt.g = zeros(Bus.n,1);     % bus conductance
Shunt.b = zeros(Bus.n,1);     % bus susceptance
if ~isempty(Shunt.con)
  Shunt.bus = Bus.int(Shunt.con(:,1));
  Shunt.con(:,5) = Settings.mva*Shunt.con(:,5)./Shunt.con(:,2);
  Shunt.con(:,6) = Settings.mva*Shunt.con(:,6)./Shunt.con(:,2);
  if Settings.octave
    Shunt.g = zeros(Bus.n,1);
    Shunt.b = zeros(Bus.n,1);
    Shunt.g(Shunt.bus) = Shunt.con(:,5);
    Shunt.b(Shunt.bus) = Shunt.con(:,6);
  else
    Shunt.g = sparse(Shunt.bus,1,Shunt.con(:,5),Bus.n,1);
    Shunt.b = sparse(Shunt.bus,1,Shunt.con(:,6),Bus.n,1);
  end
end

if ~isempty(Breaker.con)
  Breaker.line = Breaker.con(:,1);
  Breaker.bus = Bus.int(Breaker.con(:,2));
  Breaker.n = length(Breaker.con(:,1));
  % store line data
  Breaker.con(:,[9 10 11 12 13]) = ...
      Line.con(Breaker.line,[8 9 10 11 12]);
  idx = find(~Breaker.con(:,6));
  % change status of the line in case breakers are open
  if idx & ~isempty(Line.con)
    Line.con(Breaker.line(idx),[8 9 10 11 12]) = ...
        [zeros(length(idx),1), 1e6*ones(length(idx),1), ...
         zeros(length(idx),1), ones(length(idx),1), ...
         zeros(length(idx),1)]; % + Breaker.con(idx,[9 10 11 12 13]);
    % swap intervention time so that first
    % intervention will close the breaker
    Breaker.con(idx,[7 8]) = Breaker.con(idx,[8 7]);
  end
  % intervention times t = 0 are not allowed
  idx = find(Breaker.con(:,7) == 0);
  if ~isempty(idx)
    Breaker.con(idx,7) = 1e-6;
    Breaker.con(idx,8) = Breaker.con(idx,8)+1e-6;
  end
end

if Settings.octave
  DAE.J11 = zeros(Bus.n,Bus.n);
  DAE.J21 = zeros(Bus.n,Bus.n);
  DAE.J12 = zeros(Bus.n,Bus.n);
  DAE.J22 = zeros(Bus.n,Bus.n);
  DAE.Jlf = zeros(2*Bus.n,2*Bus.n);
  DAE.Jlfv = zeros(2*Bus.n,2*Bus.n);
  DAE.Jlfd = zeros(2*Bus.n,2*Bus.n);
else
  DAE.J11 = sparse(Bus.n,Bus.n);
  DAE.J21 = sparse(Bus.n,Bus.n);
  DAE.J12 = sparse(Bus.n,Bus.n);
  DAE.J22 = sparse(Bus.n,Bus.n);
  DAE.Jlf = sparse(2*Bus.n,2*Bus.n);
  DAE.Jlfv = sparse(2*Bus.n,2*Bus.n);
  DAE.Jlfd = sparse(2*Bus.n,2*Bus.n);
end
DAE.gp = zeros(Bus.n,1);
DAE.glfp = zeros(Bus.n,1);
DAE.gq = zeros(Bus.n,1);
DAE.glfq = zeros(Bus.n,1);
if length(Bus.con(1,:)) >= 4,
  Vlow  = find(Bus.con(:,3) < 0.5);
  Vhigh = find(Bus.con(:,3) > 1.5);
  if ~isempty(Vlow),
    fm_disp(['Warning: some initial guess voltage amplitudes are ' ...
             'too low.']),
  end
  if ~isempty(Vhigh),
    fm_disp(['Warning: some initial guess voltage amplitudes are ' ...
             'too high.']),
  end
  DAE.V = Bus.con(:,3);
  aref = min(abs(Bus.con(:,4)));
  alow  = find(Bus.con(:,4)-aref < -0.6);
  ahigh = find(Bus.con(:,4)-aref >  0.6);
  if ~isempty(alow),
    fm_disp(['Warning: some initial guess voltage phases are too ' ...
             'low.']),
  end
  if ~isempty(ahigh),
    fm_disp(['Warning: some initial guess voltage phases are too ' ...
             'high.']),
  end
  DAE.a = Bus.con(:,4);
else
  DAE.V = ones(Bus.n,1);
  DAE.a = zeros(Bus.n,1);
end
Bus.Pl = zeros(Bus.n,1);
Bus.Ql = zeros(Bus.n,1);
Bus.Pg = zeros(Bus.n,1);
Bus.Qg = zeros(Bus.n,1);

% count buses of different types
if (~isempty(SW.con))
  if clpsat.init, SW.store = SW.con; end
  SW.n = length(SW.con(:,1));
  SW.bus = Bus.int(round(SW.con(:,1)));
  for i = 1:SW.n
    if length(find(SW.bus == SW.bus(i))) > 1
      fm_disp(['Error: More than one slack generator ', ...
	       'connected to the same bus.'],2)
      return
    end
  end
  if length(SW.con(1,:)) < 8
    SW.con = [SW.con(:,1),100*ones(SW.n,1), ...
	      ones(SW.n,1),SW.con(:,[2:end])];
  end
  if length(SW.con(1,:)) == 5
    SW.con = [SW.con, 999*ones(SW.n,1), -999*ones(SW.n,1), ...
	      1.1*ones(SW.n,1), 0.9*ones(SW.n,1)];
  end
  if length(SW.con(1,:)) == 10
    SW.con = [SW.con, ones(SW.n,1)];
  end
  DAE.V(SW.bus) = SW.con(:,4);
  if ~sum(DAE.a) & SW.n == 1
    DAE.a(:) = SW.con(1,5);
  else
    DAE.a(SW.bus) = SW.con(:,5);
  end
else
  fm_disp('Error: No slack bus found.',2)
  return
end

if (~isempty(PV.con))
  if clpsat.init, PV.store = PV.con; end
  PV.n = length(PV.con(:,1));
  PV.bus = Bus.int(round(PV.con(:,1)));
  for i = 1:PV.n
    if length(find(PV.bus == PV.bus(i))) > 1
      fm_disp(['Error: More than one PV generator ', ...
	       'connected to the same bus.'],2)
      return
    end
  end
  DAE.V(PV.bus) = PV.con(:,5);
  if length(PV.con(1,:)) == 5
    PV.con = [PV.con, 999*ones(PV.n,1), ...
              -999*ones(PV.n,1), 1.1*ones(PV.n,1), ...
              0.9*ones(PV.n,1), ones(PV.n,1)];
  end
  if length(PV.con(1,:)) == 9; PV.con = [PV.con, ones(PV.n,1)]; end
end

% check of distributed slack bus model consistency
if ~isempty(PV.con)
  idx = find(~[PV.con(:,10); SW.con(:,11)]);
else
  idx = find(~SW.con(:,11));
end
if length(idx) == PV.n+SW.n, SW.con(:,11) = 1; end

if (~isempty(PQ.con))
  if clpsat.init, PQ.store = PQ.con; end
  PQ.n = length(PQ.con(:,1));
  PQ.bus = Bus.int(round(PQ.con(:,1)));
  for i = 1:PQ.n
    if length(find(PQ.bus == PQ.bus(i))) > 1
      fm_disp('Error: More than one PQ load connected to the same bus.',2)
      return
    end
  end
  if length(PQ.con(1,:)) == 5
    PQ.con = [PQ.con, 1.2*ones(PQ.n,1), ...
              0.8*ones(PQ.n,1), zeros(PQ.n,1)];
  elseif length(PQ.con(1,:)) == 7
    PQ.con = [PQ.con, zeros(PQ.n,1)];
  end
  for i = 1:PQ.n
    if PQ.con(i,6) == 0, PQ.con(i,6) = 1.2; end
    if PQ.con(i,7) == 0, PQ.con(i,7) = 0.8; end
  end
  PQ.P0 = PQ.con(:,4);
  PQ.Q0 = PQ.con(:,5);
end

if ~isempty(Mn.con)
  if clpsat.init, Mn.store = Mn.con; end
  Mn.n = length(Mn.con(:,1));
  Mn.bus = Bus.int(round(Mn.con(:,1)));
  Mn.init = find(~Mn.con(:,8));
end

if ~isempty(Pl.con)
  Pl.n = length(Pl.con(:,1));
  Pl.bus = Bus.int(round(Pl.con(:,1)));
  Pl.init = find(~Pl.con(:,11));
end

if ~isempty(Fl.con)
  Fl.n = length(Fl.con(:,1));
  Fl.bus = Bus.int(round(Fl.con(:,1)));
end

if ~isempty(Mot.con)
  Mot.n = length(Mot.con(:,1));
  Mot.bus = Bus.int(round(Mot.con(:,1)));
  Mot.dat = ones(Mot.n,11);
  Mot.slip = ones(Mot.n,1);
  Mot.e1r = ones(Mot.n,1);
  Mot.e1m = ones(Mot.n,1);
  Mot.e2r = ones(Mot.n,1);
  Mot.e2m = ones(Mot.n,1);
  DAE.x = zeros(Mot.n,1);
  fm_mot(0)
end