fm_ncomp.m

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

function 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-2006 Federico Milano

fm_var

Settings.ok = 1;

% 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));
% set up internal bus numbers for second indexing of buses
Bus.int(round(Bus.con(:,1)),1) = [1:Bus.n];
busmax = length(Bus.int);

% check bus voltage rates
if length(Bus.con(1,:)) < 2
  fm_disp('No voltage rates found in Bus data.',2)
  return
end
idx = find(Bus.con(:,2) == 0);
if ~isempty(idx)
  fm_disp('Some Bus voltage rate is zero! 1 kV will be used.')
  Bus.con(idx,2) = 1;
end

% check or create bus names
if length(Varname.bus) ~= Bus.n
  fm_disp('Bus names in Varname.bus does not match bus number.',2)
  Varname.bus = '';
end
if isempty(Varname.bus)
  Varname.bus = cell(Bus.n,1);
  ibusmax = fix(log10(max(Bus.con(:,1))));
  for j = 1:Bus.n
    b = int2str(Bus.con(j,1));
      zerib = [];
      for i = 1:ibusmax-length(b)
        zerib = [zerib, ' '];
      end
      Varname.bus{j} = ['Bus ', zerib, b];
  end
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(:,6);
  twt_ar = Bus.con(Bus.int(Twt.con(:,2)),[5 6]);
  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), ...
           ones(n_twt_rows,1),r1,x1,zeros(n_twt_rows,1), ...
           Twt.con(:,15), zeros(n_twt_rows,1),Twt.con(:,[16,19,22])];
  line2 = [twt_bus,Twt.con(:,[2,4,6,5]),zeros(n_twt_rows,1), ...
           Twt.con(:,6)./Twt.con(:,7),r2,x2,zeros(n_twt_rows,1), ...
           ones(n_twt_rows,1),zeros(n_twt_rows,1),Twt.con(:,[17,20,23])];
  line3 = [twt_bus,Twt.con(:,[3,4,6,5]),zeros(n_twt_rows,1), ...
           Twt.con(:,6)./Twt.con(:,8),r3,x3, zeros(n_twt_rows,1), ...
           ones(n_twt_rows,1), zeros(n_twt_rows,1),Twt.con(:,[18,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.g = sparse(Shunt.bus,1,Shunt.con(:,5),Bus.n,1);
  Shunt.b = sparse(Shunt.bus,1,Shunt.con(:,6),Bus.n,1);
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 ~isempty(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

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);
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 magnitudes are too low.'])
  end
  if ~isempty(Vhigh),
    fm_disp(['Warning: some initial guess voltage magnitudes are too high.'])
  end
  DAE.V = Bus.con(:,3);
  % check voltage phases
  aref = min(abs(Bus.con(:,4)));
  alow  = find(Bus.con(:,4)-aref < -1.5708);
  ahigh = find(Bus.con(:,4)-aref >  1.5708);
  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);

% setup components

SW = setup(SW);

% check areas
if ~isempty(Area.con)
  Area.n = length(Area.con(:,1));
  Area.slack = Bus.int(round(Area.con(:,2)));
else
  Area.con = [1, SW.con(1,1), 100, 0, 0];
  Area.n = 1;
  Area.slack = SW.bus(1);
end

PV = setup(PV);

% check of distributed slack bus model consistency
idxpv = pvgamma(PV,'sum');
idxsw = swgamma(SW,'sum');
if ~idxpv & ~idxsw, SW = setgamma(SW); end

PQ = setup(PQ);
PQgen = setup(PQgen);
PQ = addgen(PQ,PQgen);
PQgen = init(PQgen);

if ~isempty(Mn.con)
  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

Fl = setup(Fl);

if ~isempty(Mot.con)
  Mot.n = length(Mot.con(:,1));
  Mot.bus = Bus.int(round(Mot.con(:,1)));
  Mot.dat = ones(Mot.n,12);
  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

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

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