function  fm_exc(flag)% FM_EXC define Automatic Voltage Regulators%% FM_EXC(FLAG)%       FLAG = 3 differential equations%       FLAG = 4 state matrix%       FLAG = 5 non-windup limits%%see also FM_EXCIN%%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.global Bus Exc Syn DAE Oxl Pss Clustertype = Exc.con(:,2);ty1 = find(type == 1);ty2 = find(type == 2);ty3 = find(type == 3);vg = DAE.V(Exc.bus);vrmax = Exc.con(:,3);vrmin = Exc.con(:,4);Td = Exc.con(:,10);Tr = Exc.con(:,11);A = Exc.con(:,12);B = Exc.con(:,13);if ty1  vm_1  = DAE.x(Exc.vm(ty1));  vr1_1 = DAE.x(Exc.vr1(ty1));  vr2_1 = DAE.x(Exc.vr2(ty1));  efd_1 = DAE.x(Exc.vf(ty1));  m0 = Exc.con(ty1,5);  T1 = Exc.con(ty1,6);  T2 = Exc.con(ty1,7);  T3 = Exc.con(ty1,8);  T4 = Exc.con(ty1,9);  K1 = m0.*T2./T1;  K2 = m0 - K1;  K3 = T3./T4;  K4 = 1 - K3;  vr = m0.*vr2_1 + K3.*(K1.*(Exc.vrif(ty1) - vm_1) + vr1_1);endif ty2  vm_2  = DAE.x(Exc.vm(ty2));  vr1_2 = DAE.x(Exc.vr1(ty2));  vr2_2 = DAE.x(Exc.vr2(ty2));  efd_2 = DAE.x(Exc.vf(ty2));  Ka = Exc.con(ty2,5);  Ta = Exc.con(ty2,6);  Kf = Exc.con(ty2,7);  Tf = Exc.con(ty2,8);  K5 = Kf./Tf;endif ty3  vm_3  = DAE.x(Exc.vm(ty3));  vr3_3 = DAE.x(Exc.vr3(ty3));  efd_3 = DAE.x(Exc.vf(ty3));  Kr  = Exc.con(ty3,5);  T2r = Exc.con(ty3,6);  T1r = Exc.con(ty3,7);  Kr1 = Kr.*T1r./T2r;  Kr2 = Kr - Kr1;  vf0 = Exc.con(ty3,8);  v0  = Exc.con(ty3,9);endswitch flag case 3  % Updating of reference voltages  Exc.vrif = Exc.vrif0;  if Oxl.n, Exc.vrif(Oxl.exc) = Exc.vrif(Oxl.exc) - DAE.x(Oxl.v); end  if Pss.n, Exc.vrif(Pss.exc) = Exc.vrif(Pss.exc) + DAE.x(Pss.vss); end  if Cluster.n,    type = find(Cluster.con(:,2) == 1);    if ~isempty(type)      Exc.vrif(Cluster.exc(type)) = Exc.vrif(Cluster.exc(type)) + ...          DAE.x(Cluster.Vs(type));    end  end  DAE.f(Exc.vm)  = (vg - DAE.x(Exc.vm))./Tr;  % differential equations  if ty1    DAE.f(Exc.vr1(ty1)) = (K2.*(Exc.vrif(ty1) - vm_1) - vr1_1)./T1;    DAE.f(Exc.vr2(ty1)) = (K4.*(vr1_1 + K1.*(Exc.vrif(ty1) - vm_1)) ...                           - m0.*vr2_1)./(T4.*m0);    % windup limiter    vr = min(vr,vrmax(ty1));    vr = max(vr,vrmin(ty1));    DAE.f(Exc.vf(ty1)) = ...        (-efd_1+vr-ceiling(efd_1,A(ty1),B(ty1),1))./Td(ty1);  end  if ty2    DAE.f(Exc.vr1(ty2)) = (Ka.*(Exc.vrif(ty2)-vm_2-vr2_2- ...                                K5.*(efd_2))-vr1_2)./Ta;    DAE.f(Exc.vr2(ty2)) = (-K5.*efd_2-vr2_2)./Tf;    % non-windup limiter    idx = find(vr1_2 > vrmax(ty2) & DAE.f(Exc.vr1(ty2)) > 0);    DAE.f(Exc.vr1(ty2(idx))) = 0;    idx = find(vr1_2 < vrmin(ty2) & DAE.f(Exc.vr1(ty2)) < 0);    DAE.f(Exc.vr1(ty2(idx))) = 0;    vr1_2 = min(vr1_2,vrmax(ty2));    vr1_2 = max(vr1_2,vrmin(ty2));    DAE.x(Exc.vr1(ty2)) = vr1_2;    DAE.f(Exc.vf(ty2)) = (-efd_2 + vr1_2 - ...                           ceiling(efd_2,A(ty2),B(ty2),1))./Td(ty2);  end  if ty3    DAE.f(Exc.vr3(ty3)) = (Kr2.*(Exc.vrif(ty3) - vm_3) - vr3_3)./T2r;    DAE.f(Exc.vf(ty3)) = 1000*((vr3_3+Kr1.*(Exc.vrif(ty3)-vm_3)+vf0).*(vg(ty3)./v0)-efd_3);    % non-windup limiter    idx = find(efd_3 >= vrmax(ty3) & DAE.f(Exc.vf(ty3)) > 0);    if idx, DAE.f(Exc.vf(ty3)) = 0; end    DAE.x(Exc.vf(ty3)) = min(DAE.x(Exc.vf(ty3)),vrmax(ty3));    idx = find(efd_3 <= vrmin(ty3) & DAE.f(Exc.vf(ty3)) < 0);    if idx, DAE.f(Exc.vf(ty3)) = 0; end    DAE.x(Exc.vf(ty3)) = max(DAE.x(Exc.vf(ty3)),vrmin(ty3));  end case 4  % common AVR Jacobians  DAE.Fx = DAE.Fx + sparse(Exc.vm,Exc.vm,-1./Tr,DAE.n,DAE.n);  DAE.Fy = DAE.Fy + sparse(Exc.vm,Bus.n+Exc.bus,1./Tr,DAE.n,2*Bus.n);  if ty1    % AVR jacobians    DAE.Fx = DAE.Fx + sparse(Exc.vr1(ty1),Exc.vm(ty1),-K2./T1,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vr1(ty1),Exc.vr1(ty1),-1./T1,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vr2(ty1),Exc.vm(ty1),-K4.*K1./T4./m0,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vr2(ty1),Exc.vr1(ty1),K4./T4./m0,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vr2(ty1),Exc.vr2(ty1),-1./T4,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vf(ty1),Exc.vf(ty1), ...                             -(1+ceiling(efd_1,A(ty1),B(ty1),2))./ ...                             Td(ty1),DAE.n,DAE.n);    idx = find(vr < vrmax(ty1) & vr > vrmin(ty1));    if idx,      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty1(idx)),Exc.vr1(ty1(idx)), ...                               K3(idx)./Td(ty1(idx)),DAE.n,DAE.n);      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty1(idx)),Exc.vm(ty1(idx)), ...                               -K3(idx).*K1(idx)./Td(ty1(idx)),DAE.n,DAE.n);      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty1(idx)),Exc.vr2(ty1(idx)), ...                               m0(idx)./Td(ty1(idx)),DAE.n,DAE.n);    end    for i = 1:length(ty1), excjac(Exc.syn(ty1(i)),Exc.vf(ty1(i))); end  end  if ty2    % AVR jacobians    DAE.Fx = DAE.Fx + sparse(Exc.vr1(ty2),Exc.vm(ty2),-Ka./Ta,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vr1(ty2),Exc.vr1(ty2),-1./Ta,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vr1(ty2),Exc.vr2(ty2),-Ka./Ta,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vr1(ty2),Exc.vf(ty2),-K5.*Ka./Ta,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vr2(ty2),Exc.vr2(ty2),-1./Tf,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vr2(ty2),Exc.vf(ty2),-K5./Tf,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vf(ty2),Exc.vf(ty2), ...                             -(1+ceiling(efd_2,A(ty2),B(ty2),2))./ ...                             Td(ty2),DAE.n,DAE.n);    % windup limiter    idx = find(vr1_2 < vrmax(ty2) & vr1_2 > vrmin(ty2));    if idx      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty2(idx)), ...                               Exc.vr1(ty2(idx)),1./Td(ty2(idx)), ...                               DAE.n,DAE.n);    end    for i = 1:length(ty2), excjac(Exc.syn(ty2(i)),Exc.vf(ty2(i))); end  end  if ty3    % AVR jacobians    DAE.Fx = DAE.Fx + sparse(Exc.vr3(ty3),Exc.vr3(ty3),-1./T2r,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vr3(ty3),Exc.vm(ty3),-Kr2./T2r,DAE.n,DAE.n);    DAE.Fx = DAE.Fx + sparse(Exc.vf(ty3),Exc.vf(ty3),-1000,DAE.n,DAE.n);    idx = find(efd_3 < vrmax(ty3) & efd_3 > vrmin(ty3));    if ~isempty(idx)      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty3(idx)),Exc.vr3(ty3(idx)), ...                               1000.*vg(ty3(idx))./v0(idx),DAE.n,DAE.n);      DAE.Fx = DAE.Fx + sparse(Exc.vf(ty3(idx)),Exc.vm(ty3(idx)), ...                               -1000*Kr1(idx).*vg(ty3(idx))./v0(idx),DAE.n,DAE.n);      DAE.Fy = DAE.Fy + ...               sparse(Exc.vf(ty3(idx)),Bus.n+Exc.bus(ty3(idx)),1000*(vr3_3(idx)+Kr1(idx).* ...               (Exc.vrif(ty3(idx))-vm_3(idx))+vf0(idx))./v0(idx),DAE.n,2*Bus.n);      for i = 1:length(idx), excjac(Exc.syn(ty3(idx(i))),Exc.vf(ty3(idx(i)))); end    end  end case 5  global Settings  if ty2    idx = find((vr1_2 >= vrmax(ty2) | vr1_2 <= vrmin(ty2)) & ...               DAE.f(Exc.vr1(ty2)) == 0);    if ~isempty(idx)      k = Exc.vr1(ty2(idx));      DAE.tn(k) = 0;      DAE.Ac(:,k) = 0; %zeros(DAE.n + 2*Bus.n, 1);      DAE.Ac(k,:) = 0; %zeros(1, DAE.n + 2*Bus.n);      if Settings.octave        DAE.Ac(k,k) = eye(length(idx));      else        DAE.Ac(k,k) = speye(length(idx));      end    end  end  if ty3    idx = find((efd_3 >= vrmax(ty3) | efd_3 <= vrmin(ty3)) & ...               DAE.f(Exc.vf(ty3)) == 0);    if ~isempty(idx)      k = Exc.vf(ty3(idx));      DAE.tn(k) = 0;      DAE.Ac(:,k) = 0; %zeros(DAE.n + 2*Bus.n, 1);      DAE.Ac(k,:) = 0; %zeros(1, DAE.n + 2*Bus.n);      if Settings.octave        DAE.Ac(k,k) = eye(length(idx));      else        DAE.Ac(k,k) = speye(length(idx));      end    end  endend%===================================================================function output = ceiling(vf,A,B,flag)Se = A.*(exp(B.*abs(vf))-1);switch flag case 1, output = Se.*vf; case 2, output = Se+A.*B.*abs(vf).*exp(B.*abs(vf));end%===================================================================function excjac(h,k)global Syn Exc DAE% updating Syncronous Machine JacobiansTd10 = Syn.con(h,11);ord = Syn.con(h,5);switch ord case 3,  DAE.Fx(Syn.e1q(h),k) = 1/Td10; case 4,  DAE.Fx(Syn.e1q(h),k) = 1/Td10; case 5.1,  DAE.Fx(Syn.e1q(h),k) = 1/Td10; case 5.3,  xd   = Syn.con(h,8);  x1d  = Syn.con(h,9);  xL   = Syn.con(h,6);  DAE.Fx(Syn.e1q(h),k) = (xd+xL)/(x1d+xL)/Td10; otherwise,  Td20 = Syn.con(h,12);  Taa  = Syn.con(h,24);  DAE.Fx(Syn.e1q(h),k) = (1-Taa/Td10)/Td10;  DAE.Fx(Syn.e2q(h),k) = Taa/Td20/Td10;end