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:    fmilano@thunderbox.uwaterloo.ca%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_varif ~autorun('LIB Direct Method')  returnendtype   = LIB.type;slack  = LIB.slack;bus_no = LIB.selbus;% check for loaded componentsif Mn.n, mn = sum(~Mn.con(:,8));  else, mn = 0; endif Pl.n, pl = sum(~Pl.con(:,11)); else, pl = 0; endncload = mn + pl + Lines.n + Tcsc.n;if ncload | DAE.n  fm_disp('only PV, PQ and SW buses are allowed for LIB computations.')  returnendfm_disp('  ',1)fm_disp('Newton-Rapshon Method for LIB Computation - Distribuited Slack Bus',1)fm_disp(['Data file "',Path.data,File.data,'"'],1)fm_displength(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;endif bus_no < 1  fm_disp('Bus_no should be an integer > 0',2)  failed = 1;endbus_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;endif failed  DAE.n = dynordold;  PQ = PQold;  PV = PVold;  SW = SWold;  returnend% for the slack bus, it is used the power injection obtained by the power flowfor 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  endend% if no Demand.con is imposed, the load power direction% is assumed to be equal to the PQ oneif 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  endelse  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 oneif 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  endelse  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;  endend% size Jacobian matricesif DAE.n ==0  DAE.f = 0;  DAE.x = 1;  DAE.Fx = 1;endif isempty(DAE.f), DAE.f = 0; endif DAE.n == 0, DAE.n = 1; endif 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);endCjac(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 limittic;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)endfm_dispfm_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;endfor i = 1:Demand.n  k = Demand.bus(i);  Dpf1c(k,Supply.n+i) = lambda;endDxf2c = 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, endend% 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 = [];endif noSup  Supply.con = [];  Supply.n = 0;  Supply.bus = [];endSNB.init = 0;LIB.init = 1;CPF.init = 0;OPF.init = 0;