fm_busfreq.m

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

function  fm_busfreq(flag)
%FM_BUSFREQ defines bus frequency measurement
%
%
%Author:    Federico Milano
%Date:      19-Dec-2003
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global Busfreq DAE Bus Settings

x = DAE.x(Busfreq.x);
w = DAE.x(Busfreq.w);
theta = DAE.a(Busfreq.bus);
Tf = Busfreq.con(:,2);
Tw = Busfreq.con(:,3);
theta0 = Busfreq.dat(:,1);
k = Busfreq.dat(:,2);

switch flag
 case 0 % initialization

  %check time constants
  idx = find(Tf == 0);
  if idx
    busfreqwarn(idx, ['Time constant Tf cannot be zero. Tf = 0.001 ' ...
                      's will be used.'])
  end
  Busfreq.con(idx,2) = 0.001;
  idx = find(Tw == 0);
  if idx
    busfreqwarn(idx, ['Time constant Tw cannot be zero. Tw = 0.001 ' ...
                      's will be used.'])
  end
  Busfreq.con(idx,3) = 0.001;

  %variable initialization
  DAE.x(Busfreq.x) = 0;
  x = DAE.x(Busfreq.x);
  DAE.x(Busfreq.w) = 1;
  w = DAE.x(Busfreq.w);
  Busfreq.dat(:,1) = theta;
  Busfreq.dat(:,2) = 1./Tf/Settings.rad;

  fm_disp('Initialization of frequency measurement components completed.')

 case 1 % algebraic equations


 case 2 % algebraic Jacobians


 case 3 % differential equations

  DAE.f(Busfreq.x) = (k.*(theta-theta0)-x)./Tf;
  DAE.f(Busfreq.w) = (-x+k.*(theta-theta0)+1-w)./Tw;
  % non-windoup limits

 case 4 % state variable Jacobians

  % DAE.Fx
  DAE.Fx = DAE.Fx + sparse(Busfreq.x,Busfreq.x,-1./Tf,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Busfreq.w,Busfreq.x,-1./Tw,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Busfreq.w,Busfreq.w,-1./Tw,DAE.n,DAE.n);

  % DAE.Fy
  DAE.Fy = DAE.Fy + sparse(Busfreq.x,Busfreq.bus,k./Tf,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Busfreq.w,Busfreq.bus,k./Tw,DAE.n,2*Bus.n);

end

% -------------------------------------------------------------------
% function for creating warning messages
function busfreqwarn(idx, msg)
fm_disp(strcat('Warning: BUSFREQ #',int2str(idx),msg))