fm_snb.m

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

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

fm_var

if ~autorun('SNB Direct Method',0), 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.show = 0;
    CPF.nump = 200;
    CPF.step = 0.25;
    fm_cpf('snb');
    lambda = CPF.lambda;
    kg = CPF.kg;
    CPF = cpf_old;
  else
    lambda = 1;
    kg = 1e-3;
  end
else
  lambda = CPF.lambda;
  kg = CPF.kg;
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;
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

PQ = novlim(PQ,'all');
noDem = 0;
noSup = 0;
bus_gen = [SW.bus; PV.bus];
n_gen = PV.n+SW.n;
Settings.pv2pq = 0;

% if no Demand.con is used, the load power direction
% is assumed to be the one defined by the PQ loads
if Demand.n
  no_dpq = findzero(Demand);
  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 power flow routine may have convergence problems.',2)
    fm_disp
  end
else
  noDem = 1;
end

% if no Supply.con is used, the generator power direction
% is assumed to be equal to the PV one
if Supply.n
  no_sp = findzero(Supply);
  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;
end

% vector and matrix initialization
DAE.Gl = sparse(2*Bus.n,1);
DAE.Gk = sparse(2*Bus.n,1);
Nrow = sparse(1,2*Bus.n+1);
Ncol = sparse(2*Bus.n,1);
Nmat = sparse(2*Bus.n+1,2*Bus.n+1);
Kjac = sparse(1,2*Bus.n+1);
Kjac(1,SW.bus) = 1;

if ~distr, SW = setgamma(SW); end

if noDem
  Glcall(PQ);
else
  Glcall(Demand);
end

Glcall(Supply);
if distr, Gkcall(Supply); end

if noSup
  Glcall(PV);
  if distr, Gkcall(PV); end
  Glcall(SW);
end
Gkcall(SW);

if SW.n, DAE.Gl(SW.bus+Bus.n) = 0; end
if PV.n, DAE.Gl(PV.bus+Bus.n) = 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;
Z = 0;

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'})

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

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

  %DAE.Gl = sparse(2*Bus.n,1);
  %DAE.Gk = sparse(2*Bus.n,1);

  % call components functions and models
  fm_lf(1);
  fm_lf(2);
  if noDem
    glambda(PQ,lambda);
    %Glcall(PQ);
  else
    gcall(PQ);
    glambda(Demand,lambda);
    %Glcall(Demand);
  end

  glambda(Supply,lambda,distr*kg);
  %Glcall(Supply);
  %if distr, Gkcall(Supply); end

  if noSup
    glambda(PV,lambda,distr*kg)
    glambda(SW,lambda,kg)
    greactive(PV)
    %if distr, Gkcall(PV), end
    %Glcall(PV)
    Gyreactive(PV)
    %Glcall(SW)
  else
    gcall(PV);
    Gycall(PV);
    glambda(SW,1,kg)
  end
  %Gkcall(SW)
  greactive(SW)
  Gyreactive(SW)
  
  DAE.Jlfv = [DAE.J11, DAE.J12; DAE.J21, DAE.J22];

  % non-trivial condition ||w|| = 1
  norm_eq = norm([M', W', Z]);
  Vjac = [M', W', Z]/norm_eq;
  gy = [DAE.Jlfv, DAE.Gk; Kjac]'*[M; W; Z];

  % complete Jacobian matrix
  AC = [[fm_hessian([M;W]),Ncol;Nrow],[DAE.Jlfv,DAE.Gk; Kjac]',[Ncol;0]; ...
        [DAE.Jlfv,DAE.Gk;Kjac],Nmat,[DAE.Gl;0];Nrow,Vjac,0];
  
  inc = -robust*(AC\[gy; DAE.gp; DAE.gq; DAE.a(SW.bus);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;
    Z  = z_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;
    z_old = Z;

    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;
    kg = kg + inc(idx+1);                 idx = idx + 1;
    M = M + inc(idx+1:idx+Bus.n);         idx = idx + Bus.n;
    W = W + inc(idx+1:idx+Bus.n);         idx = idx + Bus.n;
    Z = Z + inc(idx+1);                   idx = idx + 1;
    lambda = lambda + inc(idx+1);

    err_max_old = err_max;
    iteration = iteration + 1;
    if Settings.show
      fm_disp(['iteration = ', int2str(iteration), ...
               ';   lambda = ',num2str(lambda), ...
               '   kg = ',num2str(kg), ...
               '   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; Z];
d1 = 2*Bus.n+1;
d2 = Supply.n;
d3 = Demand.n;
Sbus = double(Supply.bus);
Dbus = Demand.bus;
if noDem
  d3 = PQ.n;
  Dbus = double(PQ.bus);
end
if noSup
  d2 = SW.n + PV.n;
  Sbus = double([PV.bus;SW.bus]);
end
Dpfc = sparse(d1,d2+d3);
Dpfc = Dpfc + sparse(Sbus,[1:d2],-(lambda+kg),d1,d2+d3);
Dpfc = Dpfc + sparse(Dbus,d2+[1:d3],lambda,d1,d2+d3);

dl_dp = full(-W'*Dpfc/(W'*[DAE.Gl;0]))';

% Update Pg, Qg, Pl and Ql
% ===========================================================================

DAE.gp = zeros(Bus.n,1);
DAE.gq = zeros(Bus.n,1);
fm_call('pq');
glambda(Demand,lambda)
Bus.Pl = DAE.gp;
Bus.Ql = DAE.gq;
if Line.n, fm_lf(1), end
DAE.gp = zeros(Bus.n,1);
DAE.gq = zeros(Bus.n,1);
fm_call('1');
Bus.Pg = DAE.gp + DAE.glfp;
Bus.Qg = DAE.gq + DAE.glfq;

% Display results
% ===========================================================================
serv_test{1,1} = ['lambda = ',fvar(lambda,8)];
serv_test{2,1} = ['kg = ',fvar(kg,8)];
serv_test{3,1} = '  ';

if d2
  serv_test = [serv_test;strcat('d lambda / d Ps_',Varname.bus(Sbus), ...
                                '=',num2str(dl_dp(1:d2)))];
end
serv_test{end+1,1} = ' ';
if d3
  serv_test = [serv_test;strcat('d lambda / d Pd_',Varname.bus(Dbus), ...
                                '=',num2str(dl_dp(d2+[1:d3])))];
end
fm_disp
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(['Direct Method for SNB did not converge.'],2)
else
  fm_disp(['Direct Method for SNB 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 = resetvlim(PQ);
PQ = pqreset(PQ,'all');
SW = restore(SW);

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