fm_snb.m

电力系统的psat

function fm_snb
% FM_SNB compute Saddle-node Bifurcation
%        by means of a Direct Method
%
% FM_SNB
%
%System Equations:
%       g(y,kg,lambda) = 0   ->  operating point
%       dg/dy'*ro = 0        ->  transversality condition
%       |ro| - 1 = 0         ->  non-trivial condition
%
%       where:  g = load flow equations
%               y = [theta; V]
%               ro = right eigenvector
%               lambda = loading parameter
%
% Options:
%          SNB.slack:  1  ->  distribuited slack bus
%                      0  ->  single slack bus
%
% Output:  SNB.lambda: saddle node value of lambda
%          SNB.dldp:   sensitivity factors
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    11-Feb-2003
%Version:   1.0.2
%
%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('SNB Direct Method')
  return
end

% check for continuation power flow solution
% lambda the maximum lambda obtained with the CPF analysis
% or 1e-3 if no CPF analysis has been performed.
if ~CPF.init & ~clpsat.init
  Settings.ok = 0;
  uiwait(fm_choice(['It is strongly recommended running a CPF to get ' ...
      'a solution close to the SNB point. Run CPF?'],1))
  if Settings.ok
    cpf_old = CPF;
    CPF.vlim = 0;
    CPF.ilim = 0;
    CPF.qlim = 0;
    CPF.type = 2;
    CPF.linit = 0;
    CPF.nump = 50;
    fm_cpf('snb');
    lambda = CPF.lambda; 
    CPF = cpf_old;
  else
    lambda = 1e-3;
  end
else
  lambda = CPF.lambda;   
end

distr = SNB.slack;

% check for loaded components
if Mn.n, mn = sum(~Mn.con(:,8));  else, mn = 0; end
if Pl.n, pl = sum(~Pl.con(:,11)); else, pl = 0; end
ncload = mn + pl + Lines.n + Tcsc.n;
if ncload | DAE.n
  fm_disp('only PV, PQ and SW buses are allowed for SNB computations.')
  return
end

fm_disp
if distr
  fm_disp('Direct Method for SNB Computation - Distribuited Slack Bus')
else
  fm_disp('Direct Method for SNB Computation - Single Slack Bus')
end
fm_disp('                                  - Right Eigenvector')
fm_disp
fm_disp(['Data file "',Path.data,File.data,'"'])
fm_disp

length(Snapshot.V);
%DAE.V = Snapshot(1).V;
%DAE.a = Snapshot(1).ang;
%DAE.x = Snapshot(1).x;
%DAE.Jlfv = Snapshot(1).Jlfv;

Vmax_old = PQ.con(:,6);
Vmin_old = PQ.con(:,7);
PQold = PQ;
PVold = PV;
SWold = SW;
noDem = 0;
noSup = 0;
PQ.con(:,6) = 10000*ones(PQ.n,1);
PQ.con(:,7) = -10000*ones(PQ.n,1);
bus_gen = [SW.bus; PV.bus];
n_gen = PV.n+SW.n;

% for the slack bus, the power injection obtained by the power flow is used
for i = 1:SW.n,
  k = SW.bus(i);
  if ~isempty(Snapshot(1).Pg)
    pg = Snapshot(1).Pg;
    if pg(k) ~= 0
      SW.con(i,10) = pg(k);
    end
  end
end

% if no Demand.con is imposed, the load power direction
% is assumed to be the one defined by the PQ loads
if Demand.n
  no_dpq = find(Demand.con(:,3) == 0 & Demand.con(:,4) == 0);
  if ~isempty(no_dpq),
    fm_disp
    for i = 1:length(no_dpq)
      fm_disp(['No power direction found in "Demand.con" for Bus ', ...
               Varname.bus{Demand.bus(no_dpq(i))}])
    end
    fm_disp('Continuation load flow routine may have convergence problems.',2)
    fm_disp
  end
else
  noDem = 1;
  Ppos = find(PQ.con(:,4) > 0);
  Demand.con = zeros(length(Ppos),14);
  Demand.n = length(Ppos);
  Demand.bus = PQ.bus(Ppos);
  Demand.con(:,[1 2 3 4]) = PQ.con(Ppos,[1 2 4 5]);
  Demand.con(:,14) = 1;
  PQ.con(Ppos,[4 5]) = 0;
end

% if no Supply.con is imposed, the generator power direction
% is assumed to be equal to the PV one
if Supply.n
  no_sp = find(Supply.con(:,3) == 0);
  if ~isempty(no_sp),
    fm_disp
    if length(no_sp) == Supply.n
      fm_disp(['No power directions found in "Supply.con" for all buses.'])
      fm_disp('Remove "Supply" components or set power directions.')
      fm_disp('Continuation power flow interrupted',2)
      if CPF.show, set(Fig.main,'Pointer','arrow'); end
      return
    else
      for i = 1:length(no_sp)
        fm_disp(['No power direction found in "Supply.con" for Bus ', ...
                 Varname.bus{Supply.bus(no_sp(i))}])
      end
      fm_disp('Continuation power flow routine may have convergence problems.',2)
      fm_disp
    end
  end
else
  noSup = 1;
  if PV.n
    Ppos = find(PV.con(:,4) > 0);
    Supply.con = zeros(length(Ppos),14);
    Supply.n = length(Ppos);
    Supply.bus = PV.bus(Ppos);
    Supply.con(:,[1 2 3]) = PV.con(Ppos,[1 2 4]);
    Supply.con(:,14) = 1;
    PV.con(Ppos,4) = 0;
  end
  if SW.n
    Ppos = find(SW.con(:,10) > 0);
    Supply.n = Supply.n+length(Ppos);
    Supply.bus = [Supply.bus; SW.bus(Ppos)];
    Supply.con(end+[1:length(Ppos)],[1 2 3]) = SW.con(Ppos,[1 2 10]);
    Supply.con((end-length(Ppos)+1):end,14) = 1;
    SW.con(Ppos,10) = 0;
  end
end

% vector and matrix initialization
if Settings.octave
  Gl = zeros(2*Bus.n,1);
  Gk = zeros(2*Bus.n,1);
  Nrow = zeros(1,2*Bus.n);
  Ncol = zeros(2*Bus.n,1);
  Nmat = zeros(2*Bus.n,2*Bus.n);
else
  Gl = sparse(2*Bus.n,1);
  Gk = sparse(2*Bus.n,1);
  Nrow = sparse(1,2*Bus.n);
  Ncol = sparse(2*Bus.n,1);
  Nmat = sparse(2*Bus.n,2*Bus.n);
end

for i = 1:Demand.n
  k = Demand.bus(i);
  Gl(k,1) = Demand.con(i,3);
  Gl(k+Bus.n,1) = Demand.con(i,4);
end
for i = 1:Supply.n
  k = Supply.bus(i);
  Gl(k,1) = -Supply.con(i,3);
  Gk(k,1) = -Supply.con(i,3);
end
if distr == 0, Gl(SW.bus) = 0; end

[Vect,D] = eig(full(DAE.Jlfv));
diagD = diag(D);
idx = find(diagD == 1);
if ~isempty(idx)
  diagD(idx) = 1e6;
end
[val,idx] = min(abs(diagD));
Vect = inv(Vect.');
M = real(Vect(1:Bus.n,idx))+1e-4;
W = real(Vect(Bus.n+1:2*Bus.n,idx))+1e-4;

iter_max = Settings.lfmit;
iteration = 0;
tol = Settings.lftol;
g1=ones(2*Bus.n,1);
if DAE.n == 0, dynordold = 0; DAE.n = 1; end

%  Direct Method
tic

fm_status('snb','init',iter_max,{'r'},{'-'},{'g'})

G = real(Line.Y);
B = imag(Line.Y);

%lambda = CPF.lambda; 
kg = CPF.kg;
if distr, j = 1; else, j = 0; end

a_old = DAE.a;
v_old = DAE.V;
l_old = lambda;
k_old = kg;
m_old = M;
w_old = W;

err_max = tol+1;
err_max_old = err_max + 1;
robust = 1;

while err_max > tol
  if (iteration > iter_max), break, end
  if Fig.main
    if ~get(Fig.main,'UserData'), break, end
  end

  if isempty(Line.Y)
    DAE.gp = zeros(Bus.n,1);
    DAE.gq = zeros(Bus.n,1);
    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);
    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);
    end
  end

  % Jlfv computation
  if distr
    fm_lf(1); fm_pq(1); fm_pv(1); fm_lf(2); fm_pv(2);
    for i = 1:SW.n,
      k = SW.bus(i);
      DAE.gp(k) = DAE.gp(k) - SW.con(i,10);
      DAE.gq(k) = 0;
      DAE.J21(k,:)=zeros(1,Bus.n);
      DAE.J22(k,:)=zeros(1,Bus.n);
      DAE.J12(:,k)=zeros(Bus.n,1);
      DAE.J22(:,k)=zeros(Bus.n,1);
      DAE.J22(k,k)=1;
    end
    DAE.Jlfv = [DAE.J11, DAE.J12; DAE.J21, DAE.J22];
  else
    fm_call('l')
  end

  % dF/dlambda computation
  for i = 1:Demand.n
    k = Demand.bus(i);
    if isempty(find(SW.bus == k))
      DAE.gp(k) = lambda*Demand.con(i,3) + DAE.gp(k);
    end
    if isempty(find(bus_gen == k))
      DAE.gq(k) = lambda*Demand.con(i,4) + DAE.gq(k);
    end
  end
  if distr
    for i = 1:Supply.n
      k = Supply.bus(i);
      DAE.gp(k) = DAE.gp(k) - (lambda+kg)*Supply.con(i,3);
    end
  else
    for i = 1:Supply.n
      k = Supply.bus(i);
      if isempty(find(SW.bus == k))
        DAE.gp(k) = DAE.gp(k) - lambda*Supply.con(i,3);
      end
    end
  end

  % non-trivial condition ||w|| = 1
  norm_eq = norm([M', W']);
  Vjac = [M', W']/norm_eq;
  if distr
    gy = [DAE.Jlfv, Gk]'*[M; W];
  else
    gy = DAE.Jlfv'*[M; W];
  end
  DAE.gp(SW.bus) = 0;

  % complete Jacobian matrix
  if distr
    AC = [[fm_hessian([M;W]),Ncol;Nrow,0],[DAE.Jlfv,Gk]', ...
          [Ncol;0];DAE.Jlfv,Gk,Nmat,Gl;Nrow,0,Vjac,0];
  else
    AC = [fm_hessian([M;W]),DAE.Jlfv',Ncol; ...
          DAE.Jlfv,Nmat,Gl;Nrow,Vjac,0];
  end

  AC(SW.bus,:) = 0;
  AC(:,SW.bus) = 0;
  AC(SW.bus,SW.bus) = 1;
  gy(SW.bus) = 0;

  for i = 1:n_gen
    k = bus_gen(i);
    AC(k+Bus.n,:) = 0;
    AC(:,k+Bus.n) = 0;
    AC(k+Bus.n,k+Bus.n) = 1;

    AC(k+3*Bus.n+j,:) = 0;
    AC(:,k+3*Bus.n+j) = 0;
    AC(k+3*Bus.n+j,k+3*Bus.n+j) = 1;
  end

  AC(SW.bus+2*Bus.n+j,:) = 0;
  AC(:,SW.bus+2*Bus.n+j) = 0;
  AC(SW.bus+2*Bus.n+j,SW.bus+2*Bus.n+j) = 1;

  gy(Bus.n+bus_gen) = 0;
  DAE.gp(SW.bus) = 0;
  inc = -robust*(AC\[gy; DAE.gp; DAE.gq; norm_eq-1]);
  err_max = max(abs(inc));

  if err_max > 2*err_max_old & iteration

    DAE.a  = a_old;
    DAE.V  = v_old;
    lambda = l_old;
    kg = k_old;
    M  = m_old;
    W  = w_old;

    robust = 0.5*robust;
    if robust < tol, iteration = iter_max+1; break, end

  else

    a_old = DAE.a;
    v_old = DAE.V;
    l_old = lambda;
    k_old = kg;
    m_old = M;
    w_old = W;

    DAE.a = DAE.a + inc(1:Bus.n);         idx = Bus.n;
    DAE.V = DAE.V + inc(idx+1:idx+Bus.n); idx = idx + Bus.n;
    if distr, kg = kg + inc(idx+1);       idx = idx + 1; end
    M = M + inc(idx+1:idx+Bus.n);         idx = idx + Bus.n;
    W = W + inc(idx+1:idx+Bus.n);         idx = idx + Bus.n;
    lambda = lambda + inc(idx+1);

    err_max_old = err_max;
    iteration = iteration + 1;
    if Settings.show
      fm_disp(['iteration = ', int2str(iteration), ...
               ';   lambda = ',num2str(lambda), ...
               '   err = ',num2str(err_max)])
    end
    fm_status('snb','update',[iteration, err_max],iteration)
  end
end

% ===========================================================================
%                                  d lambda |
% sensistivity coefficients        -------- |
%                                     d p   |c
% ===========================================================================

W = [M; W];
if Settings.octave
  Dpfc = zeros(2*Bus.n,Demand.n+Supply.n);
else
  Dpfc = sparse(2*Bus.n,Demand.n+Supply.n);
end
for i = 1:Supply.n
  k = Supply.bus(i);
  Dpfc(k,i) = -(lambda+j*kg);
end
for i = 1:Demand.n
  k = Demand.bus(i);
  Dpfc(k,Supply.n+i) = lambda;
end

dl_dp = full(-W'*Dpfc/(W'*Gl))';

% Update Demand.con & Supply.con
%____________________________________________________________________________

Demand.con(:,7) = lambda*Demand.con(:,3);
Supply.con(:,6) = lambda*Supply.con(:,3);

% Display results
% ===========================================================================
fm_disp
fm_disp(['Lambda Critical = ',num2str(lambda)],1)
fm_disp
fm_disp('Sensitivity Coefficients',1)
fm_disp
for i = 1:Supply.n
  fm_disp([fvar(['dLambda/dPs',int2str(Supply.con(i,1)),' = '],20), ...
           fvar(dl_dp(i),8)],1);
end
fm_disp
for i = 1:Demand.n
  fm_disp([fvar(['dLambda/dPd',int2str(Demand.con(i,1)),' = '],20), ...
           fvar(dl_dp(i+Supply.n),8)],1);
end
fm_disp

serv_test{1,1} = [fvar('Lambda Critical = ',20),fvar(lambda,8)];
serv_test{2,1} = '  ';
if distr,
  serv_test{3,1} = [fvar('Kg = ',20),fvar(lambda,8)];
  serv_test{4,1} = '  ';
end

for i = 1:Supply.n
  serv_test{i+2,1} = [fvar(['dLambda/dPs',int2str(Supply.con(i,1)),' = '],20),fvar(dl_dp(i),8)];
end
serv_test{end+1,1} = ' ';
for i = 1:Demand.n
  serv_test{i+3+Supply.n,1} = [fvar(['dLambda/dPd',int2str(Demand.con(i,1)),' = '],20),fvar(dl_dp(i+Supply.n),8)];
end
DAE.n = dynordold;
Settings.iter = iteration;

Settings.lftime = toc;
if robust < tol
  fm_disp('Increment step below threshold. Bad initial point. SNB routine stopped.',2)
elseif iteration >= iter_max
  fm_disp(['Reached maximum number of iterations for Direct Method for SNB computation without convergence'],2)
else
  fm_disp(['Direct Method for SNB computation completed in ',num2str(toc),' s'])
  if Settings.showlf == 1, fm_stat(serv_test), end
end

SNB.lambda = lambda;
SNB.dldp = dl_dp;
SNB.bus = strvcat(strcat('s_',num2str(Supply.bus)), strcat('d_',num2str(Demand.bus)));

% restore original data
% ===========================================================================
fm_status('snb','close')

PQ = PQold;
PV = PVold;
SW = SWold;
if noDem
  Demand.con = [];
  Demand.n = 0;
  Demand.bus = [];
end
if noSup
  Supply.con = [];
  Supply.n = 0;
  Supply.bus = [];
end

SNB.init = 1;
LIB.init = 0;
CPF.init = 0;
OPF.init = 0;