fm_atc.m

电力系统的psat

function fm_atc
%FM_ATC determine the Available Transfer Capability (ATC)
%       by means of an iterative OPF/CPF method or a power
%       flow sensitivity analysis
%
%FM_ATC
%
%see OPF and CPF structures for settings
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Version:   1.0.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('ATC analysis')
  return
end

fm_disp

if Demand.n == 0
  fm_disp('ATC computations requires the definition of "Demand" data',2)
  return
end
if Supply.n == 0
  fm_disp('ATC computations requires the definition of "Supply" data',2)
  return
end
switch OPF.type
 case 4
  fm_disp('Determination of ATC by means of an iterative OPF-CPF method')
 case 5
  fm_disp('Determination of ATC by means of a power flow sensitivity analysis')
end
tempo1 = clock;
OPF.show = 0;
Settings.show = 0;
CPFold = CPF;
CPF.show = 0;

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

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

% ===================================================================
%  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:')
  rows = int2str(antennas);
  busf = int2str(Line.from(antennas));
  bust = int2str(Line.to(antennas));
  fm_disp(strcat('Line #',rows,' from bus #',busf,' to bus #',bust))
  fm_disp('These lines are not used for (N-1) contingency evaluations.')
else
  fm_disp('All lines are used for (N-1) contingency evaluations.')
  antennas = 0;
end

if Fig.main
  fm_disp
  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 = ' * ';

CPFold = CPF;

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

Vmin = zeros(Bus.n,1);                 Vmax = zeros(Bus.n,1);
Vmin(PQ.bus) = PQ.con(:,7);            Vmax(PQ.bus) = PQ.con(:,6);
Vmin(PV.bus) = PV.con(:,9);            Vmax(PV.bus) = PV.con(:,8);
Vmin(SW.bus) = SW.con(:,9);            Vmax(SW.bus) = SW.con(:,8);
Vmin(find(Vmin == 0)) = 0.8;           Vmax(find(Vmax == 0)) = 1.2;

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

% =======================================================================
% First OPF solution without contingencies
% =======================================================================

fm_disp('First OPF solution without faults.')
fm_opfsdr
if ~OPF.conv & Fig.main
  set(Fig.main,'UserData',0)
end
if strcmp(History.text{end},'Optimization routine interrupted.')
  return
end
fm_disp([sp,'ATC = ', num2str(OPF.atc)])

line_orig = Line.con;
Y_orig = Line.Y;
idx_old = 0;
atc_old = 0;
%atc = 0;
V = DAE.V;
ang = DAE.a;
MVA = Settings.mva;

sw = SW;
pv = PV;
pq = PQ;
branch = Line;
pq_idx = [];
for i = 1:Demand.n
  pq_idx = [pq_idx, find(PQ.bus == Demand.bus(i))];
end

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

while 1

  if Fig.main
    if ~get(Fig.main,'UserData'), break, end
  end

  lcrit = [];

  Supply.con(:,3) = Supply.con(:,6);
  Demand.con(:,3) = Demand.con(:,7);
  Demand.con(:,4) = Demand.con(:,7).*PQ.con(pq_idx,5)./PQ.con(pq_idx,4);

  %fm_disp('Continuation Power Flow Computations')

  switch OPF.type

   case 4  % (N-1) contingency analysis
    for i = 1:Line.n
      if Fig.main
        if ~get(Fig.main,'UserData'), break, end
      end

      fm_disp('Continuation Power Flow Computations')

      j = find(antennas == i);

      if isempty(j)
	fm_disp(['Line #',fvar(i,4)])
	fm_disp([sp,'from bus #',fvar(Line.from(i),5), ...
		 ' to bus #',fvar(Line.to(i),5)])
	DAE.V = OPF.Vc;
	DAE.a = OPF.ac;
        length(Snapshot.V);
	DAE.x = Snapshot(1).x;
	Line.con = line_orig;
	Line.con(i,[8 9 10]) = [0, 1e6, 0];
	fm_y;
	idx = find(abs(DAE.V-Vmax) < 1e-5);
	if ~isempty(idx), DAE.V(idx) = Vmax(idx)-1e-2; end
	CPF.init = 0;
	OPF.init = 0;
	a = fm_cpf('atc');
        %a = a+1;
	%disp(a)
	if ~isempty(a)
	  lcrit = [lcrit; [a*sum(Demand.con(:,7))*MVA+ ...
			   sum(PQ.con(:,4))*MVA, i]];
	else
	  lcrit = [lcrit; [NaN, 0]];
	end
	fm_disp([sp,'ATC = ', num2str(lcrit(end,1))])
      end
    end
    Line.con = line_orig;
    Line.Y = Y_orig;

   case 5  % sensitivity analysis: ranking (dPij/dlambda) & (Pij)

    fm_disp
    fm_disp('Sensitivity analysis.')

    lambda = OPF.guess(end-4*Bus.n);
    kg = OPF.guess(end-4*Bus.n-1-pv.n-sw.n);
    lambda = lambda - 1e-6;
    VV = OPF.Vc.*exp(jay*OPF.ac);
    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);
    Fij = abs(real(VV(Line.from).* ...
	  conj(((VV(Line.from) - tps.*VV(Line.to)).*y + ...
	  VV(Line.from).*(jay*chrg))./(tps.*conj(tps)))));
    Fji = abs(real(VV(Line.to).* ...
	  conj((VV(Line.to) - VV(Line.from)./tps).*y + ...
	  VV(Line.to).*(jay*chrg))));
    Pflow1 = max(Fij,Fji);

    PQold = PQ;
    PVold = PV;
    PQ.con(:,4) = (1+lambda)*PQ.con(:,4);
    PQ.con(:,5) = (1+lambda)*PQ.con(:,5);
    for i = 1:Demand.n
      a = find(PQ.bus == Demand.bus(i));
      PQ.con(a,4) = PQ.con(a,4) + (1+lambda)*Demand.con(i,7);
      PQ.con(a,5) = PQ.con(a,5) + ...
	  (1+lambda)*Demand.con(i,4)*Demand.con(i,7)/Demand.con(i,3);
    end
    PV.con(:,4) = (1+lambda+kg)*PV.con(:,4);
    for i = 1:Supply.n
      a = find(PV.bus == Supply.bus(i));
      PV.con(a,4) = PV.con(a,4) + (1+lambda+kg)*Supply.con(i,6);
    end
    PV.con(:,5) = OPF.Vc(PV.bus);
    SW.con(:,4) = OPF.Vc(SW.bus);
    DAE.V = OPF.Vc;
    DAE.a = OPF.ac;
    fm_spf
    fm_disp

    VV = DAE.V.*exp(jay*DAE.a);
    Fij = abs(real(VV(Line.from).* ...
	  conj(((VV(Line.from) - tps.*VV(Line.to)).*y + ...
	  VV(Line.from).*(jay*chrg))./(tps.*conj(tps)))));
    Fji = abs(real(VV(Line.to).* ...
	  conj((VV(Line.to) - VV(Line.from)./tps).*y + ...
	  VV(Line.to).*(jay*chrg))));
    Pflow2 = max(Fij,Fji);

    % sensitivity coefficients
    dPdl = (Pflow1-Pflow2)/1e-6;
    fm_disp(['Line                 |Pij(lambda)|      ', ...
	     '|Pij(lambda-dl)|   |Pij|*|d Pij/d lambda|'])
    for i = 1:Line.n
      fm_disp([fvar(i,4),'   ',fvar(Line.from(i),4),'  -> ', ...
	       fvar(Line.to(i),4),' ',fvar(Pflow1(i),19), ...
	       fvar(Pflow2(i),19),fvar(Pflow1(i)*abs(dPdl(i)),19)])
    end
    PQ = PQold;
    PV = PVold;
  end

  fm_disp
  switch OPF.type
   case 4
    [atc,idx] = min(lcrit(:,1));
    OPF.line = lcrit(idx,2);
    fm_disp(['Critical Line #',num2str(idx), ...
	     ' from bus #',fvar(Line.from(idx),5), ...
	     ' to bus #',fvar(Line.to(idx),5)])
    fm_disp([sp,'ATC = ', num2str(lcrit(idx,1))])
   case 5
    if antennas, dPdl(antennas) = 0; end
    pfactor = Pflow1.*abs(dPdl);
    [pfactor, pfactor_idx] = sort(pfactor);
    %pfactor = pfactor(end:-1:1);
    %pfactor_idx = pfactor_idx(end:-1:1);
    %disp(num2str([pfactor_idx,pfactor]))
    CPF.Pij = Pflow1;
    CPF.dPdl = dPdl;
    %return
    fm_disp('OPF solution without contingencies:')
    fm_disp([sp,'max(Pij*|dPij/dlambda)| = ', num2str(pfactor(end)), ...
	     '  (Line # ',fvar(pfactor_idx(end),4),')'])
    fm_disp([sp,'lambda = ',num2str(lambda)])
    fm_disp([sp,'ATC = ',num2str(OPF.atc)])
  end
  fm_disp

  % ============================================================================
  %  Optimal Power computation with inclusion of the fault on the critical line
  % ============================================================================

  %SW.con(:,8) = sw.con(:,8);
  %SW.con(:,9) = sw.con(:,9);
  %PV.con(:,8) = pv.con(:,8);
  %PV.con(:,9) = pv.con(:,9);
  %PQ.con(:,6) = pq.con(:,6);
  %PQ.con(:,7) = pq.con(:,7);

  fm_disp('OPF Computation with contingencies')

  switch OPF.type
   case 4
    CPF.init = 0;
    DAE.V = Snapshot(1).V;
    DAE.a = Snapshot(1).ang;
    DAE.x = Snapshot(1).x;
    fm_spf
    fm_opfsdr
    if Fig.main & ~OPF.conv, set(Fig.main,'UserData',0), end
    %fm_disp([sp,'ATC = ', num2str(OPF.atc)])
   case 5
    nopf = min(5,length(pfactor_idx));
    atc = zeros(nopf,1);
    rep = cell(nopf,1);
    for i = 1:nopf
      OPF.init = 0;
      CPF.init = 0;
      if Fig.main
        if ~get(Fig.main,'UserData'), break, end
      end
      OPF.line = pfactor_idx(i);
      DAE.V = Snapshot(1).V;
      DAE.a = Snapshot(1).ang;
      DAE.x = Snapshot(1).x;
      fm_spf
      fm_disp(['OPF computation with contingency on Line # ', ...
	       num2str(pfactor_idx(end-(i-1)))])
      fm_opfsdr
      if Fig.main & ~OPF.conv, set(Fig.main,'UserData',0), end
      atc(i) = OPF.atc;
      rep{i} = OPF.report;
      fm_disp([sp,'ATC = ', num2str(OPF.atc)])
    end
    [atc, idx_atc] = min(atc);
    OPF.report = rep{idx_atc};
    fm_disp
    idxc = pfactor_idx(end-(idx_atc-1));
    fm_disp(['Critical Line #',num2str(idxc),' from bus #', ...
	     fvar(Line.from(idxc),5), ...
	     ' to bus #',fvar(Line.to(idxc),5)])
    fm_disp([sp,'ATC = ',num2str(atc)])
  end

  % ====================================================================
  %  Convergency Criteria
  % ====================================================================

  % stop if the method uses sensitivity analysis
  if OPF.type == 5, break, end
  % stop if ATC level has not changed
  if abs(OPF.atc - atc) < 1e-2, break, end
  % stop if the worst line is always the same
  if idx == idx_old(end), break, end
  idx_old = [idx_old, idx];
  if OPF.type == 4,
    if ~isempty(OPF.report)
      rep_old(:,length(idx_old)-1) = OPF.report;
    end
    atc_old = [atc_old, OPF.atc];
  end
  % stop if the worst lines are always the same two
  % and chose the solution which has the lowest ATC
  if length(idx_old) > 3
    if idx_old(end) == idx_old(end-2),
      fm_disp
      if atc_old(end) > atc_old(end-1),
	OPF.report = rep_old(:,end-1);
	idxc = idx_old(end-1);
      else
	idxc = idx_old(end);
      end
      fm_disp(['Critical Line #',num2str(idxc), ...
	       ' from bus #',fvar(Line.from(idxc),5), ...
	       ' to bus #',fvar(Line.to(idxc),5)])
      break
    end
  end
end

Line = branch;
SW = sw;
PV = pv;
PQ = pq;
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;

if Fig.main, set(hdl,'String','Close'); end
Settings.show = 1;
OPF.show = 1;
OPF.line = 0;
CPF = CPFold;
CPF.show = 1;
CPF.init = 2;
OPF.init = 0;
LIB.init = 0;
SNB.init = 0;

if Fig.main
  if get(Fig.main,'UserData'),
    History.text = [History.text; ...
                    {' '; 'Final OPF Solution:'; ' '}; ...
                    OPF.report];
    fm_disp(['ATC = ',num2str(atc)])
    fm_disp(['ATC computation completed in ', ...
             num2str(etime(clock,tempo1)),' s'])
  else
    fm_disp('ATC computation interrupted.')
  end
end