fex_abcd.m

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

function fex_abcd
% This function computes matrices (A, B, C, D) for linear analysis by
% means of numerical differentiation.  It has to be run after state
% initialization performed by power flow PSAT functions in Command
% Line mode.
%
% Structure NLA contains Jacobian matrices and settings for the
% numeric linear analysis library.  NLA.a_sys is the state matrix, and
% can be compared with the analytical state matrix calculated by PSAT
% functions.
%
% Author:  Alberto Del Rosso
% Date:    June 2004
% Update:
% Version: preliminary version
% E-mail:  [e-mail]
% 
% I M P O R T A N T   N O T E 
% =========================== 
% This external function is not part of the original PSAT software
% package although it uses many of PSAT functions. It has not been
% approved yet by Dr. Federico Milano to be distributed along with
% PSAT files.

fm_var

fm_disp
% Check for dynamic components
if ~DAE.n 
  fm_disp('No dynamic component is loaded.')
  fm_disp('Computation of linear dynamic matrices aborted.')
  return
end
fm_disp('Computation of linear analysis matrices A,B,C,D:')
fm_disp
fm_disp('   [Delta dx/dt] = [A] [Delta x] + [B] [Delta u]')
fm_disp('       [Delta y] = [C] [Delta x] + [D] [Delta u]')
fm_disp
fm_disp('Input and Output matrices are stored in structure NLA')

% Excitation  systems
if ~Exc.n, fm_disp('No excitation systems found'), end
% Turbine governors
if ~Tg.n, fm_disp('No turbine governor found'), end
% Static var systems
if ~Svc.n, fm_disp('No Static Var Compensators found'), end

% Define matrices NLA structure
% State matrix
nstates = DAE.n;
NLA.a_sys = zeros(nstates);

% -------------------------------------------------------------------
% Output matrices  
% -------------------------------------------------------------------

% output matrix corresponding to bus voltage magnitudes
NLA.c_V  = zeros(Bus.n, DAE.n);  

% output matrix corresponding to generator active powers
NLA.c_Pg = zeros(Syn.n, DAE.n);  

% output matrix for the line power flow at the from buses Ps
NLA.c_ps = zeros(Line.n, DAE.n);  

% output matrix for the line reactive power flow at the from buses
% Qs
NLA.c_qs = zeros(Line.n, DAE.n);  

% output matrix for the line power flow at the to buses Pr
NLA.c_pr = zeros(Line.n, DAE.n);  

% output matrix for the line reactive power flow at the to buses Qr
NLA.c_qr = zeros(Line.n, DAE.n);  

% output matrix for the line current magnitude at the from buses Is
NLA.c_Is = zeros(Line.n, DAE.n);  

% output matrix for the line current magnitude at the to buses Ir
NLA.c_Ir = zeros(Line.n, DAE.n);  

% -------------------------------------------------------------------
% Input matrices
% -------------------------------------------------------------------

% input matrix corresponding to the exciter reference Vrif0 
NLA.b_Vr = zeros(DAE.n, Exc.n);  

% input matrix corresponding to governor reference
NLA.b_Tr = zeros(DAE.n, Tg.n);

% input matrix corresponding to the svc reference input
NLA.b_svc= zeros(DAE.n, Svc.n);  

% Set tolerance for variable perturbation
% ----------------------------------------------------------------------
tol = NLA.tol;

% Backup of actual variables
% ----------------------------------------------------------------------
if isempty(DAE.x), DAE.x = 0; end
if isempty(DAE.f), DAE.f = 0; end

% refresh equations
fm_call('i')
fm_nrlf(Settings.lfmit,Settings.lftol,0);

xa    = DAE.x;
anga  = DAE.a;
Va    = DAE.V;
ffn   = DAE.f;
ga    = DAE.g;
Pg0   = Syn.Pg;
Vrif0 = Exc.vrif0;

% Initial line power flows
[Ps0,Qs0,Pr0,Qr0,Is0,Ir0] = fex_lineflows;

% Comuputation of state matrix "NLA.a_sys" and output matrices

% Sequentially perturb differential and algebraic equations 

% A matrix
% ----------------------------------------------------------------------
for j=1:nstates;
  deltaX = max(tol,abs(tol*DAE.x(j)));
  DAE.x(j)=DAE.x(j)+deltaX;
  % --------------------------------------------------------------------
  % Solve algebraic equations  g(x0,y) for algebraic variables y,
  % where x0 is the perturbed state vector which is
  % kept constant during the NR loop.
  % --------------------------------------------------------------------
  fm_nrlf(Settings.lfmit,Settings.lftol,0);
  % Differential equations & state matrix calculation
  ff = DAE.f;
  NLA.a_sys(:,j)=(ff-ffn)./deltaX;
  % Power flows and current magnitude in  Lines 
  [Ps,Qs,Pr,Qr,Is,Ir] = fex_lineflows;
  % Output matrices matrices calculation 
  NLA.c_V(:,j) =(DAE.V-Va)/deltaX;
  NLA.c_ps(:,j)=(Ps-Ps0)/deltaX; 
  NLA.c_qs(:,j)=(Qs-Qs0)/deltaX; 
  NLA.c_pr(:,j)=(Pr-Pr0)/deltaX; 
  NLA.c_qr(:,j)=(Qr-Qr0)/deltaX; 
  NLA.c_Is(:,j)=(Is-Is0)/deltaX; 
  NLA.c_Ir(:,j)=(Ir-Ir0)/deltaX; 
  if Syn.n, NLA.c_Pg(:,j)=(Syn.Pg-Pg0)/deltaX; end
  DAE.x  = xa;
  DAE.a  = anga;
  DAE.V  = Va;   
  Syn.Pg = Pg0;
end

% -----------------------------------------------------------------------
% Computation of input matrices B
% -----------------------------------------------------------------------

% B matrix for Exciters reference voltages 
% -----------------------------------------------------------------------
fm_nrlf(Settings.lfmit,Settings.lftol,0);
if Exc.n
  incY=0; 
  nexcs=Exc.n;
  for j=1:nexcs;
    deltaU = max(tol,abs(tol*Exc.vrif0(j)));
    Exc.vrif0(j) = Exc.vrif0(j)+deltaU;

    % --------------------------------------------------------------------
    % Solve algebraic equations -- solve equations g(x,y,U0) for algebraic
    % variables y, where U0 is the perturbed input vector which is
    % kept constant during the NR loop.  
    % --------------------------------------------------------------------
    
    % Initialize NR loop
    fm_nrlf(Settings.lfmit,Settings.lftol,0);
    % Differential equations & state matrix formation    
    ff=DAE.f;
    NLA.b_Vr(:,j)=(ff-ffn)./deltaU;
    DAE.x = xa;
    DAE.a = anga;
    DAE.V = Va;   
    Exc.vrif0 = Vrif0;
  end
end

% B matrix for Governor reference  
numdiff(Tg,ffn,xa,Va,anga);

% B matrix SVC voltage reference  
numdiff(Svc,ffn,xa,Va,anga);

fm_disp('Numerical Linear Analysis completed.')