fm_vs.m

电力系统的psat

function    fm_vs(flag,k)
% FM_VS save Jacobian Matrices during time domain integrations
%
% FM_VS(FLAG,K)
%      FLAG   0 -> initialization and memory allocation
%             1 -> enlarge memory allocation
%             2 -> computes determinants, eigenvalues
%                  and participation factors
%             3 -> adjust memory allocation at end of time domain
%                  integrations
%      K integration step number
%
%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.

global Settings Bus DAE Varout Exc

if flag == 0

   chunk = Settings.chunk;
   Varout.detJlf    = zeros(chunk,1);
   Varout.detJlfv   = zeros(chunk,1);
   Varout.detJlfd   = zeros(chunk,1);
   Varout.AutoState = zeros(chunk,DAE.n);
   Varout.AutoJlfr  = zeros(chunk,Bus.n);
   Varout.AutoJlfdr = zeros(chunk,Bus.n);
   Varout.pf_lf     = zeros(chunk,Bus.n,Bus.n);
   Varout.pf_lfd    = zeros(chunk,Bus.n,Bus.n);

elseif flag == 1

   % ampliamento dei vettori utilizzati per la valutazione della stabilità di tensione.

   chunk = Settings.chunk;
   Varout.detJlf    = [Varout.detJlf; zeros(chunk,1)];
   Varout.detJlfv   = [Varout.detJlfv; zeros(chunk,1)];
   Varout.detJlfd   = [Varout.detJlfd; zeros(chunk,1)];
   Varout.AutoState = [Varout.AutoState; zeros(chunk,DAE.n)];
   Varout.AutoJlfr  = [Varout.AutoJlfr; zeros(chunk,Bus.n)];
   Varout.AutoJlfdr = [Varout.AutoJlfdr; zeros(chunk,Bus.n)];
   Varout.pf_lf     = [Varout.pf_lf; zeros(chunk,Bus.n,Bus.n)];
   Varout.pf_lfd    = [Varout.pf_lfd; zeros(chunk,Bus.n,Bus.n)];

elseif flag == 2

   Fxr = DAE.Fx;
   Fyr = DAE.Fy;
   Gxr = DAE.Gx;

   pf = zeros(Bus.n,Bus.n);
   As = DAE.Fx - DAE.Fy*(DAE.Jlfv\DAE.Gx);
   Varout.AutoState(k,:) = sort([round(eig(full(As))*1e5)/1e5]');
   %[c, vjlf] = condest(DAE.Jlf);
   Jlfcond = DAE.Jlfv; % + sparse(1:2*Bus.n,1:2*Bus.n,vjlf);
   Varout.detJlf(k) = det(DAE.Jlf);
   Varout.detJlfv(k) = det(DAE.Jlfv);
   Jlfr = (Jlfcond(Bus.n+1:2*Bus.n,Bus.n+1:2*Bus.n) - Jlfcond(Bus.n+1:2*Bus.n, 1:Bus.n)*...
      inv(Jlfcond(1:Bus.n,1:Bus.n))*Jlfcond(1:Bus.n,Bus.n+1:2*Bus.n));
   [Vlfr,Dlfr] = eig(full(Jlfr));
   Dlfr_vec = diag(Dlfr)';
   [Varout.AutoJlfr(k,:),Ilfr] = sort(Dlfr_vec);
   pf = Vlfr.*(inv(Vlfr))';
   Varout.pf_lf(k,:,:) = pf(:,Ilfr);
   DAE.Jlfd = (DAE.Jlfv - Gxr*(Fxr\Fyr));
   %[c, vjlfd] = condest(DAE.Jlfd);
   %Jlfd = Jlfd; % + sparse(1:2*Bus.n,1:2*Bus.n,vjlfd);
   Varout.detJlfd(k) = det(DAE.Jlfd);
   Jlfdr = (DAE.Jlfd(Bus.n+1:2*Bus.n,Bus.n+1:2*Bus.n) - DAE.Jlfd(Bus.n+1:2*Bus.n, 1:Bus.n)*...
      inv(DAE.Jlfd(1:Bus.n,1:Bus.n))*DAE.Jlfd(1:Bus.n,Bus.n+1:2*Bus.n));
   [Vlfdr,Dlfdr] = eig(full(Jlfdr));
   Dlfdr_vec = diag(Dlfdr)';
   [Varout.AutoJlfdr(k,:),Ilfdr] =  sort(Dlfdr_vec);
   pf = Vlfdr.*(inv(Vlfdr))';
   Varout.pf_lfd(k,:,:) = pf(:,Ilfdr);

elseif  flag == 3

   Varout.detJlf = Varout.detJlf(1:k);
   Varout.detJlfv = Varout.detJlfv(1:k);
   Varout.detJlfd = Varout.detJlfd(1:k);
   Varout.AutoState = Varout.AutoState(1:k,:);
   Varout.AutoJlfr = Varout.AutoJlfr(1:k,:);
   Varout.AutoJlfdr = Varout.AutoJlfdr(1:k,:);
   Varout.pf_lf = Varout.pf_lf(1:k,:,:);
   Varout.pf_lfd = Varout.pf_lfd(1:k,:,:);;

end