fm_cpf.m

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

function lambda_crit = fm_cpf(fun)
% FM_CPF continuation power flow routines for computing nose curves
%        and determining critical points (saddle-node bifurcations)
%
% [LAMBDAC] = FM_CPF
%
%       LAMBDAC: loading paramter at critical or limit point
%
%see also CPF structure and FM_CPFFIG
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    16-Sep-2003
%Version:   1.0.1
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

fm_var

if ~autorun('Continuation Power Flow',0), return, end

lambda_crit = [];

%  Settings
if CPF.ilim
  ilim = CPF.flow;
else
  ilim = 0;
end
type = double(~CPF.sbus);
if type & SW.n == 1 & Supply.n == 1
  if Supply.bus ~= SW.bus
    fm_disp(' * * ')
    fm_disp('Only one Supply found. Single slack bus model will be used.')
    type = 0;
  end
end
if ~PQ.n & ~Demand.n
  fm_disp('* * * The systems does not present PQ load or Demand data.')  
  fm_disp('* * * Continuation Power Flow analysis interrupted.')  
  return
end
stop = CPF.type - 1;
vlim = CPF.vlim;
qlim = CPF.qlim;
perp = ~(CPF.method - 1);

if strcmp(fun,'atc') | strcmp(fun,'gams')
  one = 1;
else
  one = 0;
end

%if PV.n+SW.n == 1, type = 0; end

if CPF.show
  fm_disp
  if type
    fm_disp('Continuation Power Flow - Distribuited Slack Bus')
  else
    fm_disp('Continuation Power Flow - Single Slack Bus')
  end
  fm_disp(['Data file "',Path.data,File.data,'"'])
  fm_disp
  if perp
    fm_disp('Corrector Method:         Perpendicular Intersection')
  else
    fm_disp('Corrector Method:         Local Parametrization')
  end
  if vlim
    fm_disp('Check Voltage Limits:     Yes')
  else
    fm_disp('Check Voltage Limits:     No')
  end
  if qlim
    fm_disp('Check Generator Q Limits: Yes')
  else
    fm_disp('Check Generator Q Limits: No')
  end
  if ilim
    fm_disp('Check Flow Limits:        Yes')
  else
    fm_disp('Check Flow Limits:        No')
  end
  fm_disp
end

% Initializations of vectors and matrices
% ----------------------------------------------------------------------

% disable conversion to impedance for PQ loads
PQ = novlim(PQ,'all');

nodyn = 0;

if Syn.n
  Pm_syn = Syn.pm;
  pm_idx = Syn.omega;
  syn_idx = [1:Syn.n]';
  if Tg.n
    pm_idx(Tg.syn) = [];
    syn_idx(Tg.syn) = [];
  end
end

% initialization of the state vector and Jacobian matices
if DAE.n > 0
  DAE.f = ones(DAE.n,1);  % state variable time derivatives
  DAE.x = ones(DAE.n,1);  % state variables
  fm_xfirst;
  DAE.Fx = sparse(DAE.n,DAE.n);   % Dx(f)
  DAE.Fy = sparse(DAE.n,2*Bus.n); % Dy(f)
  DAE.Gx = sparse(2*Bus.n,DAE.n); % Dx(g)
  Flambda = sparse(DAE.n,1);
  Fk = sparse(DAE.n,1);
else  % no dynamic components
  nodyn = 1;
  DAE.n = 1;
  DAE.f = 0;
  DAE.x = 0;
  DAE.Fx = 1;
  DAE.Fy = sparse(1,2*Bus.n);
  DAE.Gx = sparse(2*Bus.n,1);
  Flambda = 0;
  Fk = 0;
end

PV2PQ = Settings.pv2pq;
noDem = 0;
noSup = 0;
sp = ' * ';

Settings.pv2pq = 0;

if ilim
  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);
  Imax = Line.con(:,12+ilim);
  switch ilim
   case 1, fw = 'I ';
   case 2, fw = 'P ';
   case 3, fw = 'S ';
  end
  Imax(find(Imax <= 0)) = 99999;  % no limit for undefined limit
end

% if no Demand.con is found, the load direction
% is assumed to be the one of the PQ load
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 load flow routine may have convergence problems.',2)
    fm_disp
  end
else
  noDem = 1;
  idx = findpos(PQ);
  Demand = add(Demand,dmdata(PQ,idx));
  PQ = pqzero(PQ,idx);
end

% if no Supply.con is found, the load direction
% is assumed to be the one of the PV or the SW generator
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.'])
      if noDem
        fm_disp('Supply data will be ignored.')
        noSup = 1;
        Supply = remove(Supply,[1:Supply.n]);
        if qlim & type, SW = move2sup(SW); end
        %SW = move2sup(SW);
        %PV = move2sup(PV);        
      else
        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
      end
    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 qlim & type, SW = move2sup(SW); end
  %SW = move2sup(SW); 
  %PV = move2sup(PV);
end

% Newton-Raphson routine settings
iter_max = Settings.lfmit;
iterazione = 0;
tol =  CPF.tolc;
tolf = CPF.tolf;
tolv = CPF.tolv;
err_max = tol+1;
proceed = 1;
sigma_corr = 1;
if DAE.n == 0, DAE.n = 1; end

Kjac = sparse(1,2*Bus.n+DAE.n+2);
Ljac = sparse(1,2*Bus.n+DAE.n+2);
kjac = sparse(1,2*Bus.n+DAE.n+1);
Kjac(1,DAE.n+Settings.refbus) = 1;
kjac(1,DAE.n+Settings.refbus) = 1;

% chose a PQ to display the CPF progress in the main window
PQidx = pqdisplay(PQ);

if ~PQidx % absence of PQ buses
  PQidx = pqdisplay(Demand);
  if ~PQidx
    if ~perp
      perp = 1;
      fm_disp('* * * Switch to perpendicular intersection * * * ')
    end
    PQidx = 1;
  end
end

% ---------------------------------------------------------------------
% Continuation Power Flow Routine
% ---------------------------------------------------------------------

tic

fm_status('cpf','init',12,{'m'},{'-'},{Theme.color11},[0 1.2])

l_vect = [];

lambda = CPF.linit;
lambda_crit = CPF.linit;
kg = 0;
lambda_old = -1;
Jsign = 0;
Jsign_old = 0;
Sflag = 1;
ksign = 1;
Qmax_idx = [];
Qmin_idx = [];
Vmax_idx = [];
Vmin_idx = [];
Iij_idx = [];
Iji_idx = [];
Iij = ones(Line.n,1);
Iji = ones(Line.n,1);
Qswmax_idx = [];
Qswmin_idx = [];

pqlim(PQ,vlim,0,0,0);
fm_out(0,0,0)
if nodyn
  DAE.Fx = 1; 
  Varout.idx = Varout.idx+1; 
end
fm_out(1,0,0)

% initial Jacobian Gk 
DAE.Gk = sparse(2*Bus.n,1);
if type, Gkcall(Supply), end
if type, Gkcall(PV), end
Gkcall(SW)

if Syn.n
  Fk(pm_idx) = Pm_syn(syn_idx)./Syn.con(syn_idx,18);
end

DAE.kg = 0;
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);
fm_lf(2);
if noSup
  Gyreactive(PV)
else
  Gycall(PV);
end
Gyreactive(SW)
DAE.Jlfv = [DAE.J11, DAE.J12; DAE.J21, DAE.J22];

% First predictor step
% ---------------------------------------------------------------

Jsign = 1; 
count_qmin = 0;
kqmin = 1;
d_lambda = 0;
lambda_old = CPF.linit;
inc = zeros(DAE.n+2*Bus.n+2,1);
Ljac(end) = 1;

% critical point
V_crit = DAE.V;
a_crit = DAE.a;
x_crit = DAE.x;
k_crit = kg;
l_crit = lambda;

while err_max > tol

  cpfmsg = [];

  % corrector step
  % ---------------------------------------------------------------
  
  V_old = DAE.V;
  ang_old = DAE.a;
  x_old = DAE.x;
  kg_old = kg;
  lambda_old = lambda;
  corrector = 1;
  if PV.n, inc(DAE.n+Bus.n+PV.bus) = 0; end
  if SW.n, inc(DAE.n+Bus.n+SW.bus) = 0; end

  while corrector

    modinc = norm(inc);
    if Fig.main
      if ~get(Fig.main,'UserData'), break, end
    end

    iter_corr = 0;
    err_corr = 1;

    while err_corr > tol

      if (iter_corr > 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 component functions and models
      fm_lf(1);
      fm_lf(2);
      gcall(PQ)
      DAE.Jlf = [DAE.J11, DAE.J12; DAE.J21, DAE.J22];

      glambda(Demand,lambda)
      glambda(Supply,lambda,type*kg)
      
      Glcall(Demand)
      Glcall(Supply)
      if type, Gkcall(Supply), end
      Gkcall(Supply)
      
      if Syn.n
        DAE.f(pm_idx) = DAE.f(pm_idx)+ ...
            (lambda+kg-1)*Pm_syn(syn_idx)./Syn.con(syn_idx,18);
        Flambda(pm_idx) = Pm_syn(syn_idx)./Syn.con(syn_idx,18);
        Fk(pm_idx) = Pm_syn(syn_idx)./Syn.con(syn_idx,18);
        DAE.Jlfv(Settings.refbus,:) = 0;
        DAE.Jlfv(:,Settings.refbus) = 0;
        DAE.Jlfv(Settings.refbus,Settings.refbus) = 1;
        DAE.gp(Settings.refbus) = 0;
      end
      
      if noSup
        glambda(PV,lambda,type*kg)
        glambda(SW,lambda,kg)
        greactive(PV)
        if type, 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];
      if SW.n
        DAE.Fy(:,Bus.n+SW.bus) = 0;
        DAE.Gx(Bus.n+SW.bus,:) = 0;
        DAE.Gl(SW.bus+Bus.n) = 0;
      end
      if PV.n
        DAE.Gl(PV.bus+Bus.n) = 0;
        DAE.Fy(:,Bus.n+PV.bus) = 0;
        DAE.Gx(Bus.n+PV.bus,:) = 0;
      end

      if perp*iterazione
        Cinc = sigma_corr*inc;        
        cont_eq = Cinc'*([DAE.x;DAE.a;DAE.V;kg;lambda]- ...
			 [x_old; ang_old; V_old;kg_old; ...
			  lambda_old]-Cinc);
        inc_corr = -[DAE.Fx, DAE.Fy, Fk, Flambda; DAE.Gx, DAE.Jlfv, ...
                     DAE.Gk, DAE.Gl; Cinc'; Kjac]\[DAE.f; DAE.gp; ...
                            DAE.gq; cont_eq; DAE.a(Settings.refbus)];
      else
        if iterazione
	  cont_eq = DAE.V(PQidx)-sigma_corr*inc(PQidx+DAE.n+Bus.n) - ...
		    V_old(PQidx);
        else
	  cont_eq = lambda - sigma_corr*d_lambda - lambda_old;
	end
        inc_corr = -[DAE.Fx, DAE.Fy, Fk, Flambda; DAE.Gx, DAE.Jlfv, ...
                     DAE.Gk, DAE.Gl; Ljac; Kjac]\[DAE.f; DAE.gp; ...
                            DAE.gq; cont_eq; DAE.a(Settings.refbus)];
      end
      
      DAE.x = DAE.x + inc_corr(1:DAE.n);
      DAE.a = DAE.a + inc_corr(1+DAE.n: Bus.n+DAE.n);