fm_excin.m

电力系统的psat

function fm_excin
% FM_EXCIN initialize AVRs state variables and reference voltages
%
%see also FM_EXC
%
%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 DAE Syn Exc Varname Settings

if ~Exc.n, return, end

% AVR types
type = Exc.con(:,2);
ty1 = find(type == 1);
ty2 = find(type == 2);
ty3 = find(type == 3);

% bus voltages, field voltages & general parameters
vg = DAE.V(Exc.bus);
vf = Syn.vf(Exc.syn);
vrmax = Exc.con(:,3);
vrmin = Exc.con(:,4);
Tr = Exc.con(:,11);

% synchronous machine dynamic orders
ord = Syn.con(Exc.syn,5);

% general warnings
idx = find(ord == 2);
if idx
  excwarn(idx,' is connected to a generator of order 2.')
end
idx = find(Tr <= 0);
if idx
  Tr(idx) = 0.001;
  Exc.con(idx,11) = 0.001;
  excwarn(idx,[' Tr is negative or zero. ', ...
	       'Default value Tr = 0.001 will be used.'])
end

if ty1
  A = Exc.con(ty1,12);
  B = Exc.con(ty1,13);
  Ce = -vf(ty1).*(1+A.*(exp(B.*abs(vf(ty1)))-1));
  m0 = Exc.con(ty1,5);
  Td = Exc.con(ty1,10);
  T1 = Exc.con(ty1,6);
  T2 = Exc.con(ty1,7);
  T3 = Exc.con(ty1,8);
  T4 = Exc.con(ty1,9);
  idx = find(m0 <= 0);
  if idx
    m0(idx) = 400;
    Exc.con(ty1(idx),5) = 400;
    excwarn(idx, [' m0 cannot be zero. ', ...
	    'Default value m0 = 400 will be used.'])
  end
  idx = find(Td <= 0);
  if idx
    Td(idx) = 1;
    Exc.con(ty2(idx),10) = 1;
    excwarn(idx, [' Te cannot be zero. ', ...
	    'Default value Te = 1 will be used.'])
  end
  idx = find(T1 <= 0);
  if idx
    T1(idx) = 0.1;
    Exc.con(ty1(idx),6) = 0.1;
    excwarn(idx, [' T1 cannot be zero. ', ...
	    'Default value T1 = 0.1 will be used.'])
  end
  idx = find(T4 <= 0);
  if idx
    T4(idx) = 0.01;
    Exc.con(ty1(idx),9) = 0.01;
    excwarn(idx, [' T4 cannot be zero. ', ...
		  'Default value T4 = 0.01 will be used.'])
  end
  K1 = m0.*T2./T1;
  K2 = m0 - K1;
  K3 = T3./T4;
  K4 = 1 - K3;
  DAE.x(Exc.vm(ty1))  = vg(ty1);
  Exc.vrif(ty1) = vg(ty1)-Ce./m0;
  DAE.x(Exc.vr1(ty1)) = K2.*(Exc.vrif(ty1)-vg(ty1));
  DAE.x(Exc.vr2(ty1)) = K4.*(Exc.vrif(ty1)-vg(ty1));
  DAE.x(Exc.vf(ty1)) = vf(ty1);
  Exc.vrif0(ty1) = Exc.vrif(ty1);
  vr = m0.*DAE.x(Exc.vr2(ty1))+ ...
       K3.*(K1.*(Exc.vrif(ty1)-vg(ty1))+DAE.x(Exc.vr1(ty1)));
  idx = find(vr > vrmax(ty1));
  if idx
    excwarn(idx,' Vr is over its max limit.')
  end
  idx = find(vr < vrmin(ty1));
  if idx
    excwarn(idx,' Vr is under its min limit.')
  end
end
if ty2
  A = Exc.con(ty2,12);
  B = Exc.con(ty2,13);
  Ce = vf(ty2).*(1+A.*(exp(B.*abs(vf(ty2)))-1));
  Ka = Exc.con(ty2,5);
  Ta = Exc.con(ty2,6);
  Kf = Exc.con(ty2,7);
  Tf = Exc.con(ty2,8);
  Td = Exc.con(ty2,10);
  idx = find(Tf <= 0);
  if idx
    Tf(idx) = 0.1;
    Exc.con(ty2(idx),8) = 0.1;
    excwarn(idx,[' Tf cannot be zero. ', ...
	    'Default value Tf = 0.1 will be used.'])
  end
  idx = find(Td <= 0);
  if idx
    Td(idx) = 1;
    Exc.con(ty2(idx),10) = 1;
    excwarn(idx,[' Te cannot be zero. ', ...
	    'Default value Te = 1 will be used.'])
  end
  idx = find(Ta <= 0);
  if idx
    Exc.con(ty2(idx),6) = 0.1;
    excwarn(idx,[' Ta cannot be zero. ', ...
	    'Default value Ta = 0.1 will be used.'])
  end
  DAE.x(Exc.vm(ty2))  = vg(ty2);
  DAE.x(Exc.vr1(ty2)) = Ce;
  DAE.x(Exc.vr2(ty2)) = -Kf.*vf(ty2)./Tf;
  DAE.x(Exc.vf(ty2)) = vf(ty2);
  Exc.vrif(ty2) = Ce./Ka+vg(ty2);
  Exc.vrif0(ty2) = Exc.vrif(ty2);
  idx = find(Ce > vrmax(ty2));
  if idx
    excwarn(idx,' Vr1 is over its max limit.')
  end
  idx = find(Ce < vrmin(ty2));
  if idx
    excwarn(idx,' Vr1 is under its min limit.')
  end
end
if ty3
  T2 = Exc.con(ty3,6);
  Exc.con(ty3,8) = Syn.vf(Exc.syn(ty3)); % offset vf0
  Exc.con(ty3,9) = vg(ty3);              % initial voltage v0
  Exc.vrif(ty3)  = vg(ty3);
  Exc.vrif0(ty3) = vg(ty3);
  DAE.x(Exc.vm(ty3))  = vg(ty3);
  DAE.x(Exc.vr3(ty3)) = 0;
  idx = find(T2 <= 0);
  if idx
    Exc.con(ty3(idx),6) = 0.1;
    excwarn(idx,[' T2 cannot be zero. ', ...
	    'Default value T2 = 0.1 will be used.'])
  end
  idx = find(Syn.vf(Exc.syn(ty3)) > vrmax(ty3));
  if idx
    excwarn(idx,' Vf is over its max limit.')
  end
  idx = find(Syn.vf(Exc.syn(ty3)) < vrmin(ty3));
  if idx
    excwarn(idx,' Vf is under its min limit.')
  end
  DAE.x(Exc.vf(ty3)) = Syn.vf(Exc.syn(ty3));
end

fm_disp('Initialization of Automatic Voltage Regulators completed.')

% function for creating warning messages
function excwarn(idx, msg)
global Exc
fm_disp(strcat('Warning: AVR #',int2str(idx), ...
	       ' at bus #',int2str(Exc.bus(idx)),msg))