fm_limit.m

电力系统的psat

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-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('LIB Direct Method')
  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 + Tcsc.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;

PQold = PQ;
PVold = PV;
SWold = SW;
noDem = 0;
noSup = 0;
Vmax = PQ.con(:,6);
Vmin = PQ.con(:,7);
PQ.con(:,6) = 10000*ones(PQ.n,1);
PQ.con(:,7) = -10000*ones(PQ.n,1);
n_gen = PV.n+SW.n;
bus_gen = sort([SW.bus; PV.bus]);
Qmin = [SW.con(:,7); PV.con(:,7)];
Qmax = [SW.con(:,6); PV.con(:,6)];
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 = find(PQ.bus(:,1) == bus_no);
  if isempty(a)
    fm_disp('No PQ load found for the specified bus number',2)
    failed = 1;
  end
  eta = min(Vmax(a));
 case 2
  a = find(PQ.bus(:,1) == bus_no);
  if isempty(a)
    fm_disp('No PQ load found for the specified bus number',2)
    failed = 1;
  end
  eta = max(Vmin(a));
 case 3
  a = find(PV.bus == bus_no);
  b = find(SW.bus == bus_no);
  if isempty(a) & isempty(b),
    fm_disp('No generator found for the specified bus number',2)
    failed = 1;
  end
  if a
    PQ.con = [PQ.con; PV.con(a,[1 2 3]),-PV.con(a,[4 6]), ...
              10000, -10000, 0];
    PQ.bus = [PQ.bus; PV.bus(a)];
    PQ.n = PQ.n+1;
    eta = PV.con(a,5);
    PV.con(a,:) = [];
    PV.bus(a) = [];
    PV.n = PV.n-1;
  end
  if b
    PQ.con = [PQ.con; SW.con(b,[1 2 3]),-Bus.Pg(SW.bus(b)), ...
              -SW.con(b,6), 10000, -10000, 0];
    PQ.bus = [PQ.bus; SW.bus(b)];
    PQ.n = PQ.n+1;
    eta = SW.con(b,4);
    SW.con(b,:) = [];
    SW.bus(b) = [];
    SW.n = SW.n-1;
  end
 case 4
  a = find(PV.bus == bus_no);
  b = find(SW.bus == bus_no);
  if isempty(a) & isempty(b),
    fm_disp('No generator found for the specified bus number',2)
    failed = 1;
  end
  if a
    PQ.con = [PQ.con; PV.con(a,[1 2 3]),-PV.con(a,[4 7]), ...
              10000, -10000, 0];
    PQ.bus = [PQ.bus; PV.bus(a)];
    PQ.n = PQ.n+1;
    eta = PV.con(a,5);
    PV.con(a,:) = [];
    PV.bus(a) = [];
    PV.n = PV.n-1;
  end
  if b
    PQ.con = [PQ.con; SW.con(b,[1 2 3]),-Bus.Pg(SW.bus(b)), ...
              -SW.con(b,7), 10000, -10000, 0];
    PQ.bus = [PQ.bus; SW.bus(b)];
    PQ.n = PQ.n+1;
    eta = SW.con(b,4);
    SW.con(b,:) = [];
    SW.bus(b) = [];
    SW.n = SW.n-1;
  end
 otherwise
  fm_disp('ERROR: option "',num2str(type),'" is not defined.',2)
  failed = 1;
end

if failed
  DAE.n = dynordold;
  PQ = PQold;
  PV = PVold;
  SW = SWold;
  return
end

% for the slack bus, it is used the power injection obtained by the power flow
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 equal to the PQ one
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

% 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

if Settings.octave
  Flambda = zeros(DAE.n,1);
  Fk = zeros(DAE.n,1);
  Fc = zeros(DAE.n,1);
  Kjac = zeros(1,2*Bus.n+DAE.n+2);
  Cjac = zeros(1,2*Bus.n+DAE.n+2);
  Gl = zeros(2*Bus.n,1);
  Gk = zeros(2*Bus.n,1);
else
  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);
  Gl = sparse(2*Bus.n,1);
  Gk = sparse(2*Bus.n,1);
end

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

  % call components functions and models
  fm_lf(1); fm_pq(1); fm_pv(1); fm_lf(2); fm_pv(2);
  for i = 1:SW.n,
    k = Settings.refbus(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];
  DAE.Jlfv(:,Settings.refbus) = 0;

  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);
      Gl(k,1) = Demand.con(i,3);
    end
    if isempty(find(bus_gen == k)),
      DAE.gq(k) = lambda*Demand.con(i,4) + DAE.gq(k);
      Gl(k+Bus.n,1) = Demand.con(i,4);
    end
    %DAE.gp(k) = lambda*Demand.con(i,3) + DAE.gp(k);
    %DAE.gq(k) = lambda*Demand.con(i,4) + DAE.gq(k);
    %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);
    DAE.gp(k) = DAE.gp(k) - lambda*Supply.con(i,3);
    Gl(k,1) = -Supply.con(i,3);
  end
  if slack
    for i = 1:Supply.n
      k = Supply.bus(i);
      DAE.gp(k) = DAE.gp(k) - kg*Supply.con(i,3);
      Gk(k,1) = -Supply.con(i,3);
    end
  else
    DAE.gp(Settings.refbus) = DAE.gp(Settings.refbus) - ...
        kg*Supply.con(find(Supply.bus==Settings.refbus),3);
    Gk(Settings.refbus,1) = -Supply.con(find(Supply.bus==Settings.refbus),3);
  end

  inc = -[DAE.Fx, DAE.Fy, Flambda, Fk; DAE.Gx, DAE.Jlfv, Gl, Gk; Cjac; Kjac]\ ...
        [DAE.f; DAE.gp; DAE.gq; DAE.V(bus_no)-eta; 0];
  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,Gk;k_jac];
Dlf1c = [Gl;0];
Dpf1c = sparse(2*Bus.n+1,Demand.n+Supply.n);
for i = 1:Supply.n
  k = Supply.bus(i);
  Dpf1c(k,i) = -lambda;
  Dpf1c(end,i) = kg;
end
for i = 1:Demand.n
  k = Demand.bus(i);
  Dpf1c(k,Supply.n+i) = lambda;
end

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

% 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 = 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 = 0;
LIB.init = 1;
CPF.init = 0;
OPF.init = 0;