fm_limit.m

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

function fm_limit
% FM_LIMIT compute Limit-Induced Bifurcation (LIB)
%          by means of a Newton-Raphson routine.
%
% FM_LIMIT
%
%     LIB.type:  1 = 'Vmax' for maximum voltage limit
%                2 = 'Vmin' for minimum voltage limit
%                3 = 'Qmax' for maximum reactive power limit
%                4 = 'Qmin' for minimum reactive power limit
%     LIB.slack:  0   ->   single slack bus
%                 1   ->   distribuited slack bus
%     LIB.bus:    bus number at which the limit will be applied
%     LIB.lambda: the critical value of lambda at the LIB eq. point
%
%                                              d P   |
%     LIB.dpdl:  sensitivity coefficient    -------- |
%                                           d lambda |0
%
%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-2006 Federico Milano

fm_var

if ~autorun('LIB Direct Method',0), return, end

type   = LIB.type;
slack  = LIB.slack;
bus_no = LIB.selbus;

% 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 LIB computations.')
  return
end

fm_disp('  ',1)
fm_disp('Newton-Rapshon Method for LIB Computation - Distribuited Slack Bus',1)
fm_disp(['Data file "',Path.data,File.data,'"'],1)
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;

dynordold = DAE.n;

noDem = 0;
noSup = 0;
PQ = novlim(PQ,'all');
n_gen = PV.n+SW.n;
bus_gen = sort([SW.bus; PV.bus]);
[Qmax,Qmin] = fm_qlim(99,-99,'gen');
[Vmax,Vmin] = fm_vlim(1.2,0.8);
failed = 0;

if bus_no > Bus.n
  fm_disp('Bus_no exceeds bus number',2)
  failed = 1;
end
if bus_no < 1
  fm_disp('Bus_no should be an integer > 0',2)
  failed = 1;
end
bus_no = Bus.int(round(bus_no));

switch type
 case 1
  a = findbus(PQ,bus_no);
  if isempty(a)
    fm_disp('No PQ load found for the specified bus number',2)
    failed = 1;
  end
  eta = Vmax(a);
 case 2
  a = findbus(PQ,bus_no);
  if isempty(a)
    fm_disp('No PQ load found for the specified bus number',2)
    failed = 1;
  end
  eta = Vmin(a);
 case 3
  a = findbus(PV,bus_no);
  b = findbus(SW,bus_no);
  if isempty(a) & isempty(b),
    fm_disp('No generator found for the specified bus number',2)
    failed = 1;
  end
  if a
    PQ = add(PQ,[PV.con(a,[1 2 3]),-PV.con(a,[4 6]),10,-10,0]);
    eta = getvg(PV,a);
    PV = remove(PV,a);
  end
  if b
    PQ = add(PQ,[SW.con(b,[1 2 3]),-Bus.Pg(SW.bus(b)),-SW.con(b,6),10,-10,0]);
    eta = getvg(SW,b);
    SW = remove(SW,b);
    SW = add(SW,move2sw(PV));
  end
 case 4
  a = findbus(PV,bus_no);
  b = findbus(SW,bus_no);
  if isempty(a) & isempty(b),
    fm_disp('No generator found for the specified bus number',2)
    failed = 1;
  end
  if a
    PQ = add(PQ,[PV.con(a,[1 2 3]),-PV.con(a,[4 7]),10,-10,0]);
    eta = getvg(PV,a);
    PV = remove(PV,a);
  end
  if b
    PQ = add(PQ,[SW.con(b,[1 2 3]),-Bus.Pg(SW.bus(b)),-SW.con(b,7),10,-10,0]);
    eta = getvg(SW,b);
    SW = remove(SW,b);
    SW = add(SW,move2sw(PV));
  end
 otherwise
  fm_disp('ERROR: option "',num2str(type),'" is not defined.',2)
  failed = 1;
end

if failed
  DAE.n = dynordold;
  PQ = restore(PQ);
  PV = restore(PV);
  SW = restore(SW);
  return
end

% if no Demand.con is imposed, the load power direction
% is assumed to be equal to the PQ one
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;
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 = 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;
  if slack, SW = move2sup(SW); end 
end

% size Jacobian matrices
if DAE.n ==0
  DAE.f = 0;
  DAE.x = 1;
  DAE.Fx = 1;
end
if isempty(DAE.f), DAE.f = 0; end
if DAE.n == 0, DAE.n = 1; end

Flambda = sparse(DAE.n,1);
Fk = sparse(DAE.n,1);
Fc = sparse(DAE.n,1);
Kjac = sparse(1,2*Bus.n+DAE.n+2);
Cjac = sparse(1,2*Bus.n+DAE.n+2);
Cjac(DAE.n+Bus.n+bus_no) = 1;
Kjac(1,DAE.n+Settings.refbus) = 1;

iter_max = Settings.lfmit;
iterazione = 0;
tol = Settings.lftol;
err_max = tol+1;

%  Power Flow Routine with inclusion of limit
tic;

fm_status('lib','init',iter_max,{'b'},{'-'},{'y'},[-1 5])

l_vect = [];
lambda = 1;
kg = 0;

while err_max > tol
  if (iterazione >= 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,slack*kg);
  Glcall(Supply);
  if slack, Gkcall(Supply); end
  
  if noSup
    glambda(PV,lambda,slack*kg)
    glambda(SW,lambda,kg)
    greactive(PV)
    if slack, 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];

  inc = -[DAE.Fx, DAE.Fy, Flambda, Fk; DAE.Gx, DAE.Jlfv, DAE.Gl, DAE.Gk; Cjac; Kjac]\ ...
        [DAE.f; DAE.gp; DAE.gq; DAE.V(bus_no)-eta; DAE.a(Settings.refbus)];
  DAE.x = DAE.x + inc(1:DAE.n);
  DAE.a = DAE.a + inc(1+DAE.n: Bus.n+DAE.n);
  DAE.V = DAE.V + inc(DAE.n+Bus.n+1:end-2);
  lambda = lambda + inc(end-1);
  kg = kg + inc(end);

  err_max = max(abs(inc));
  iterazione = iterazione + 1;
  fm_status('lib','update',[iterazione, err_max],iterazione)
  fm_disp(['iteration = ',int2str(iterazione), ...
           '    lambda = ',num2str(lambda), ...
           '    kg = ',num2str(kg), ...
           '    err = ',num2str(err_max)],1)
end

fm_disp
fm_disp(['lambda critical = ',num2str(lambda)])

% sensitivity coefficients
% ===========================================================================
k_jac = Kjac([DAE.n+1:end-2, end]);
Dxf1c = [DAE.Jlfv,DAE.Gk;k_jac];
Dlf1c = [DAE.Gl;0];
d1 = 2*Bus.n+1;
d2 = Demand.n+Supply.n;
Dpf1c = sparse(d1,d2);
Dpf1c = Dpf1c + sparse(Supply.bus,[1:Supply.n],-(lambda+kg),d1,d2);
Dpf1c = Dpf1c + sparse(Demand.bus,Supply.n+[1:Demand.n],lambda,d1,d2);

Dxf2c = Dxf1c;
Dxf2c(bus_no,:) = zeros(1,2*Bus.n+1);
Dxf2c(:,bus_no) = zeros(2*Bus.n+1,1);
Dxf2c(bus_no,bus_no) = 1;
Dlf2c = Dlf1c;
Dpf2c = Dpf1c;

mu = Dlf2c - Dxf2c*(Dxf1c\Dlf1c);

dl_dp = full((mu')*(Dxf2c*(Dxf1c\Dpf1c) - Dpf2c)/(mu'*mu))';

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

% 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
% ===========================================================================
Settings.lftime = toc;
Settings.iter = iterazione;
if iterazione >= iter_max
  fm_disp(['Reached Maximum Number of Iterations for ', ...
           'LIB computation without Convergence'],2)
else
  fm_disp(['Limit Induced Bifurcation computed in ', ...
           num2str(Settings.lftime),' s'],1)
  if Settings.showlf == 1, fm_stat, end
end

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

DAE.n = dynordold;
PQ = restore(PQ);
PV = restore(PV);
SW = restore(SW);

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