fm_n1cont.m

电力系统的psat

function Pout = fm_n1cont
% FM_N1CONT compute power limits in transmission lines
%           with an N-1 contingency criterion
%
%PMAX = FM_N1CONT (output stored in CPF.pmax)
%
%see also FM_SNB
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    05-July-2004
%Version:   1.1.0
%
%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.

fm_var

if ~autorun('(N-1) Contingency Analysis'), return, end

fm_disp
fm_disp('N-1 Contingency Computation')

tempo1 = clock;
Settings.show = 0;
fm_set('lf')

CPFold = CPF;
CPF.show = 0;
length(Snapshot.V);

Demand_old = Demand;
Supply_old = Supply;
PV_old = PV;
SW_old = SW;
PQ_old = PQ;

% Continuation Power Flow settings
% ------------------------------------------------------------------

CPF.method = 1;
CPF.flow = 1;
CPF.type = 3;
CPF.sbus = 0;
CPF.vlim = 1;
CPF.ilim = 1;
CPF.qlim = 1;
CPF.linit = 0;

% ==================================================================
% Determination of "antennas", i.e. lines that connects a
% PQ or a PV bus to the rest of the network
% ==================================================================

antennas = [];
fromto = [];
for i = 1:Bus.n
  Iidx = fm_iidx(i,Line.con);
  if length(Iidx(:,1)) == 1,
    antennas = [antennas; Iidx(2)];
    if Iidx(5) == 1
      fromto = [fromto; Iidx(4)];
    else
      fromto = [fromto; Iidx(3)];
    end
  end
end

fm_disp
if ~isempty(antennas)
  fm_disp('Detected the following antennas:')
  for i = 1:length(antennas)
    k = antennas(i);
    fm_disp(['Line #',fvar(k,4),' from ', ...
             fvar(Varname.bus{Line.from(k)},12), ...
	     ' to ',fvar(Varname.bus{Line.to(k)},12)])
  end
  fm_disp(['When these lines are out, connected generators ', ...
           'and/or loads will be eliminated.'])
else
  fm_disp('All lines are used for (N-1) contingency evaluations.')
  antennas = 0;
end

fm_disp
if Fig.main
  if Hdl.status ~= 0, delete(Hdl.status), Hdl.status = 0; end
  hdl = findobj(Fig.main,'Tag','PushClose');
  set(hdl, ...
      'String','STOP', ...
      'ForegroundColor','w', ...
      'FontWeight','bold', ...
      'BackgroundColor',Theme.color08);
  set(Fig.main,'UserData',1);
end
sp = ' * ';

% ====================================================================
% Saving complete impedance matrix and voltages
% ====================================================================

line_orig = Line.con;
Y_orig = Line.Y;
idx_old = 0;
V = DAE.V;
ang = DAE.a;

% ====================================================================
% Continuation Power Flow loop for the (N-1) contingency evaluations
% ====================================================================

lcrit = [];

if Supply.n
  if Supply.con(:,3) == 0
    Supply.con(:,3) = Supply.con(:,6);
  end
end
if Demand.n
  if Demand.con(:,3) == 0,
    Demand.con(:,3) = Demand.con(:,7);
    Demand.con(:,4) = Demand.con(:,7).*PQ.con(:,5)./PQ.con(:,4);
  end
end

fm_disp('Continuation Power Flow Computations')
Pij = [];
Pji = [];

tps = Line.con(:,11).*exp(jay*Line.con(:,12)*pi/180);
r = Line.con(:,8);
rx = Line.con(:,9);
chrg = Line.con(:,10)/2;
z = r + jay*rx;
y = 1./z;
g = real(y);
b = imag(y);
VV = DAE.V.*exp(jay*DAE.a);

Fijbc = VV(Line.from).*conj(((VV(Line.from) - tps.*VV(Line.to)).*y + ...
    VV(Line.from).*(jay*chrg))./(tps.*conj(tps)));
Fjibc = VV(Line.to).*conj((VV(Line.to) - VV(Line.from)./tps).*y + ...
    VV(Line.to).*(jay*chrg));

Fijbc = abs(real(Fijbc));
Fjibc = abs(real(Fjibc));
Fmax = max(Fijbc,Fjibc);

% (N-1) contingency analysis
for i = 1:Line.n
  if Fig.main
    if ~get(Fig.main,'UserData'), break, end
  end
  SW = SW_old;
  PV = PV_old;
  PQ = PQ_old;
  Demand = Demand_old;
  Supply = Supply_old;

  j = find(antennas == i);

  if ~isempty(j)

    SWbus = find(SW.bus == fromto(j));
    PVbus = find(PV.bus == fromto(j));
    PQbus = find(PQ.bus == fromto(j));
    Dbus  = find(Demand.bus == fromto(j));
    Sbus  = find(Supply.bus == fromto(j));

    if ~isempty(SWbus)
      SW.con(SWbus,:) = [];
      [pvmax, pvidx] = max(PV.con(:,4));
      SW.con(SWbus,:) = [PV.con(pvidx,[1 2 3 5]), ang(PV.bus(pvidx)), ...
			 PV.con(pvidx,[6 7 8 9 4 10])];
      PV.n = PV.n - 1;
      SW.bus = PV.bus(pvidx);
      PV.bus(pvidx) = [];
      PV.con(pvidx,:) = [];
    end
    if ~isempty(PVbus)
      PV.con(PVbus,:) = [];
      PV.bus(PVbus) = [];
      PV.n = PV.n - 1;
    end
    if ~isempty(PQbus)
      PQ.con(PQbus,:) = [];
      PQ.bus(PQbus) = [];
      PQ.n = PQ.n - 1;
    end
    if ~isempty(Dbus)
      Demand.con(Dbus,:) = [];
      Demand.bus(Dbus,:) = [];
      Demand.n = Demand.n - 1;
    end
    if ~isempty(Sbus)
      Supply.con(Sbus,:) = [];
      Supply.bus(Sbus,:) = [];
      Supply.n = Supply.n - 1;
    end
  end

  fm_disp(['Line #',fvar(i,4)])
  fm_disp([sp,'from ',fvar(Varname.bus{Line.from(i)},12),' to ', ...
	   fvar(Varname.bus{Line.to(i)},12)])
  Line.con = line_orig;
  Line.Y = Y_orig;
  Line.con(i,[8 9 10]) = [0, 1e6, 0];
  fm_y;
  a = [];
  guess = 0;
  while isempty(a)
    if guess > 20, break, end
    DAE.V = V;
    DAE.a = ang;
    DAE.x = Snapshot(1).x;
    CPF.init = 0;
    a = fm_cpf('n1cont');
    CPF.linit = CPF.linit-0.1;
    guess = guess + 1;
  end
  CPF.linit = 0;
  if ~isempty(a) & ~isempty(DAE.V)
    lcrit = [lcrit; [a, i]];
    % power flows in transmission lines
    VV = DAE.V.*exp(jay*DAE.a);
    Fij = VV(Line.from).*conj(((VV(Line.from) - tps.*VV(Line.to)).*y + ...
			       VV(Line.from).*(jay*chrg))./(tps.*conj(tps)));
    Fji = VV(Line.to).*conj((VV(Line.to) - VV(Line.from)./tps).*y + ...
			    VV(Line.to).*(jay*chrg));
    % take out the power flow in the line with contingency
    Fij(i) = NaN;
    Fji(i) = NaN;

    Fij = max(abs(real(Fij)),Fmax);
    b = find(Fij < 1.2*Fijbc);
    if ~isempty(b), Fij(b) = NaN; end

    Fji = max(abs(real(Fji)),Fmax);
    b = find(Fji < 1.2*Fjibc);
    if ~isempty(b), Fji(b) = NaN; end

    Pij = [Pij, Fij];
    Pji = [Pji, Fji];
  else
    lcrit = [lcrit; [NaN, 0]];
    Pij = [Pij, NaN*zeros(Line.n,1)];
    Pji = [Pji, NaN*zeros(Line.n,1)];
  end
  fm_disp([sp,'ATC = ', num2str(lcrit(end,1))])

end

Pmax = min(Pij',Pji');
if nargout == 1, Pout = Pmax; end
[Pmax, Pidx] = min(Pmax);

Header{1,1}{1,1} = 'N-1 CONTINGENCY ANALYSIS';
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-2005';
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)];

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

Header{2,1} = 'POWER FLOW LIMITS';
Cols{2,1} = {' Line',' Line',' Pij',' Pij max'; ...
             ' ',' Outage',' [p.u.]',' [p.u.]'};
for i = 1:Line.n
  Rows{2,1}{i,1} = [num2str(Bus.int(Line.from(i))),'-', ...
                    num2str(Bus.int(Line.to(i)))];
  Rows{2,1}{i,2} = [num2str(Bus.int(Line.from(Pidx(i)))),'-', ...
                    num2str(Bus.int(Line.to(Pidx(i))))];
end
Matrix{2,1} = [Fmax,Pmax'];

% writing data...
filename = [fm_filenum(Settings.export),['.',Settings.export]];
switch Settings.export
 case 'txt'
  fm_writetxt(Matrix,Header,Cols,Rows,filename)
 case 'xls'
  fm_writexls(Matrix,Header,Cols,Rows,filename)
 case 'tex'
  fm_writetex(Matrix,Header,Cols,Rows,filename)
end

if Fig.main
  set(hdl,'BackgroundColor',Theme.color02, ...
          'String','Close', ...
          'ForegroundColor','k', ...
          'FontWeight','normal');
end
Settings.show = 1;

CPF = CPFold;
Line.con = line_orig;
Line.Y = Y_orig;
Demand = Demand_old;
Supply = Supply_old;
PV = PV_old;
SW = SW_old;
PQ = PQ_old;
CPF.pmax = Pmax';
CPF.init = 3;

if Fig.main
  if ~get(Fig.main,'UserData'),
    fm_disp('N-1 contingency computation interrupted.')
  else
    fm_disp(['N-1 contingency computation completed in ', ...
             num2str(etime(clock,tempo1)),' s'])
  end
else
  fm_disp(['N-1 contingency computation completed in ', ...
           num2str(etime(clock,tempo1)),' s'])
end