fm_n1cont.m

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

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

fm_var

if ~autorun('(N-1) Contingency Analysis',0)
  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);

% 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')
  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 = [];

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 = 1.2*max(Fijbc,Fjibc);
lcrit = zeros(Line.n,2);
Pij = zeros(Line.n,Line.n);
Pji = zeros(Line.n,Line.n);

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

  j = find(antennas == i);

  if ~isempty(j)

    idx_sw = findbus(SW,fromto(j));
    idx_pv = findbus(PV,fromto(j));
    idx_pq = findbus(PQ,fromto(j));
    idx_de = findbus(Demand,fromto(j));
    idx_su = findbus(Supply,fromto(j));

    SW = remove(SW,idx_sw);
    if ~isempty(idx_sw), SW = add(SW,move2sw(PV)); end
    PV = remove(PV,idx_pv);
    PQ = remove(PQ,idx_pq,'force');
    Demand = remove(Demand,idx_de);
    Supply = remove(Supply,idx_su);

  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) & abs(a-CPF.linit) > 1e-4

    lcrit(i,:) = [a,1];

    % 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;

    Pij(:,i) = abs(real(Fij));
    Pji(:,i) = abs(real(Fji));
  else
    lcrit(i,:) = [NaN,0];
    Pij(:,i) = NaN*zeros(Line.n,1);
    Pji(:,i) = 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-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)];

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

Header{2,1} = 'POWER FLOW LIMITS';
Cols{2,1} = {' Line',' Outage of', ' Worst case',' Pij base',' Pij max'; ...
             ' ',' this line',' line outage',' [p.u.]',' [p.u.]'};
for i = 1:Line.n
  Rows{2,1}{i,1} = [num2str(Bus.con(Line.from(i),1)),'-', ...
                    num2str(Bus.con(Line.to(i),1))];
  Rows{2,1}{i,3} = [num2str(Bus.con(Line.from(Pidx(i)),1)),'-', ...
                    num2str(Bus.con(Line.to(Pidx(i)),1))];
  if lcrit(i,2)
    Rows{2,1}{i,2} = ' Feasible';
  else
    Rows{2,1}{i,2} = ' Unfeasible';    
  end
end
Matrix{2,1} = [Fmax,Pmax'];

% writing data...
fm_write(Matrix,Header,Cols,Rows)

if Fig.main, set(hdl,'String','Close'), end
Settings.show = 1;

CPF = CPFold;
Line.con = line_orig;
Line.Y = Y_orig;

SW = restore(SW);
PV = restore(PV);
PQ = restore(PQ);
Demand = restore(Demand);
Supply = restore(Supply);

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