fm_syn.m

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

function fm_syn(flag)
% FM_SYN defines Synchronous Machines
%
% FM_SYN(FLAG)
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%       FLAG = 5 non-windup limits
%
%see also FM_SYNIT
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderboxx.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global Bus Syn DAE Exc Settings Tg

ord = Syn.con(:,5);

is2  = find(ord == 2);
is3  = find(ord == 3);
is4  = find(ord == 4);
is51 = find(ord == 5.1);
is52 = find(ord == 5.2);
is53 = find(ord == 5.3);
is6  = find(ord == 6);
is8  = find(ord == 8);

bs2  = Syn.bus(is2);
bs3  = Syn.bus(is3);
bs4  = Syn.bus(is4);
bs51 = Syn.bus(is51);
bs52 = Syn.bus(is52);
bs53 = Syn.bus(is53);
bs6  = Syn.bus(is6);
bs8  = Syn.bus(is8);

delta = DAE.x(Syn.delta);
omega = DAE.x(Syn.omega);
e1q = zeros(Syn.n,1);
e1d = zeros(Syn.n,1);
e2q = zeros(Syn.n,1);
e2d = zeros(Syn.n,1);
psiq = zeros(Syn.n,1);
psid = zeros(Syn.n,1);

if (~isempty(is3)),
  e1q(is3) = DAE.x(Syn.e1q(is3));
end
if (~isempty(is4)),
  e1d(is4) = DAE.x(Syn.e1d(is4));
  e1q(is4) = DAE.x(Syn.e1q(is4));
end
if (~isempty(is51)),
  e1d(is51) = DAE.x(Syn.e1d(is51));
  e1q(is51) = DAE.x(Syn.e1q(is51));
  e2d(is51) = DAE.x(Syn.e2d(is51));
end
if (~isempty(is52)),
  e1q(is52) = DAE.x(Syn.e1q(is52));
  e2q(is52) = DAE.x(Syn.e2q(is52));
  e2d(is52) = DAE.x(Syn.e2d(is52));
end
if (~isempty(is53)),
  e1q(is53) = DAE.x(Syn.e1q(is53));
  psid(is53) = DAE.x(Syn.psid(is53));
  psiq(is53) = DAE.x(Syn.psiq(is53));
end
if (~isempty(is6)),
  e1d(is6) = DAE.x(Syn.e1d(is6));
  e1q(is6) = DAE.x(Syn.e1q(is6));
  e2d(is6) = DAE.x(Syn.e2d(is6));
  e2q(is6) = DAE.x(Syn.e2q(is6));
end
if (~isempty(is8)),
  e1d(is8) = DAE.x(Syn.e1d(is8));
  e1q(is8) = DAE.x(Syn.e1q(is8));
  e2d(is8) = DAE.x(Syn.e2d(is8));
  e2q(is8) = DAE.x(Syn.e2q(is8));
  psid(is8) = DAE.x(Syn.psid(is8));
  psiq(is8) = DAE.x(Syn.psiq(is8));
end

ag = DAE.a(Syn.bus);
ss = sin(delta-ag);
cc = cos(delta-ag);

iM = 1./Syn.con(:,18);
D = Syn.con(:,19);

ra   = Syn.con(:,7);    xl   = -Syn.con(:,6);
xd   = Syn.con(:,8);    xq   = Syn.con(:,13);
xd1  = Syn.con(:,9);    xq1  = Syn.con(:,14);
xd2  = Syn.con(:,10);   xq2  = Syn.con(:,15);
Td10 = Syn.con(:,11);   Tq10 = Syn.con(:,16);
Td20 = Syn.con(:,12);   Tq20 = Syn.con(:,17);
Kw   = Syn.con(:,20);   Kp   = Syn.con(:,21);

if ~isempty(is3)
  a34 = 1./Td10(is3);
  a35 = a34.*(xd(is3)-xd1(is3));
end
if ~isempty(is4)
  a44 = 1./Td10(is4);
  a45 = a44.*(xd(is4)-xd1(is4));
  b43 = 1./Tq10(is4);
  b44 = b43.*(xq(is4)-xq1(is4));
end
if ~isempty(is51)
  gq = xd1(is51)./xq1(is51).*Tq20(is51)./Tq10(is51).*(xq(is51)-xq1(is51));
  a514 = 1./Td10(is51);
  a515 = a514.*(xd(is51)-xd1(is51));
  b511 = 1./Tq20(is51);
  b512 = b511.*(xq1(is51)-xd1(is51)+gq);
  b513 = 1./Tq10(is51);
  b514 = b513.*(xq(is51)-xq1(is51)-gq);
end
if ~isempty(is52)
  Taa = Syn.con(:,24);
  gd = xd2(is52)./xd1(is52).*Td20(is52)./Td10(is52).*(xd(is52)-xd1(is52));
  a521 = 1./Td20(is52);
  a522 = a521.*(xd1(is52)-xd2(is52)+gd);
  a523 = Taa(is52)./Td10(is52)./Td20(is52);
  a524 = 1./Td10(is52);
  a525 = a524.*(xd(is52)-xd1(is52)-gd);
  a526 = a524.*(1-Taa(is52)./Td10(is52));
  b521 = 1./Tq20(is52);
  b522 = b521.*(xq(is52)-xq2(is52));
end
if ~isempty(is53)
  a531 = (xd(is53)-xd1(is53))./(xd(is53)+xl(is53));
  a532 = 1./(1-a531);
  a534 = 1./Td10(is53);
  c534 = 1./(xd(is53)+xl(is53));
  c535 = 1./(xq(is53)+xl(is53));
end
if ~isempty(is6)
  Taa = Syn.con(:,24);
  gd = xd2(is6)./xd1(is6).*Td20(is6)./Td10(is6).*(xd(is6)-xd1(is6));
  gq = xq2(is6)./xq1(is6).*Tq20(is6)./Tq10(is6).*(xq(is6)-xq1(is6));
  a1 = 1./Td20(is6);
  a2 = a1.*(xd1(is6)-xd2(is6)+gd);
  a3 = Taa(is6)./Td10(is6)./Td20(is6);
  a4 = 1./Td10(is6);
  a5 = a4.*(xd(is6)-xd1(is6)-gd);
  a6 = a4.*(1-Taa(is6)./Td10(is6));
  b1 = 1./Tq20(is6);
  b2 = b1.*(xq1(is6)-xq2(is6)+gq);
  b3 = 1./Tq10(is6);
  b4 = b3.*(xq(is6)-xq1(is6)-gq);
end
if ~isempty(is8)
  Taa = Syn.con(:,24);
  gd = xd2(is8)./xd1(is8).*Td20(is8)./Td10(is8).*(xd(is8)-xd1(is8));
  gq = xq2(is8)./xq1(is8).*Tq20(is8)./Tq10(is8).*(xq(is8)-xq1(is8));
  a18 = 1./Td20(is8);
  a28 = a18.*(xd1(is8)-xd2(is8)+gd);
  a38 = Taa(is8)./Td10(is8)./Td20(is8);
  a48 = 1./Td10(is8);
  a58 = a48.*(xd(is8)-xd1(is8)-gd);
  a68 = a48.*(1-Taa(is8)./Td10(is8));
  b18 = 1./Tq20(is8);
  b28 = b18.*(xq1(is8)-xq2(is8)+gq);
  b38 = 1./Tq10(is8);
  b48 = b38.*(xq(is8)-xq1(is8)-gq);
  xd2(is8) = xd2(is8)+xl(is8);
  xq2(is8) = xq2(is8)+xl(is8);
end

switch flag
 case 1 % active & reactive powers

  Syn.Id = -DAE.V(Syn.bus).*(Syn.c1.*ss+Syn.c3.*cc);
  Syn.Iq =  DAE.V(Syn.bus).*(Syn.c2.*ss-Syn.c1.*cc);
  if ~isempty(is2)
    Syn.Id(is2) = Syn.Id(is2) + Syn.c3(is2).*Syn.vf(is2);
    Syn.Iq(is2) = Syn.Iq(is2) + Syn.c1(is2).*Syn.vf(is2);
  end
  if ~isempty(is3)
    Syn.Id(is3) = Syn.Id(is3) + Syn.c3(is3).*e1q(is3);
    Syn.Iq(is3) = Syn.Iq(is3) + Syn.c1(is3).*e1q(is3);
  end
  if ~isempty(is4)
    Syn.Id(is4) = Syn.Id(is4) + Syn.c1(is4).*e1d(is4) + Syn.c3(is4).*e1q(is4);
    Syn.Iq(is4) = Syn.Iq(is4) - Syn.c2(is4).*e1d(is4) + Syn.c1(is4).*e1q(is4);
  end
  if ~isempty(is51)
    Syn.Id(is51) = Syn.Id(is51) + Syn.c1(is51).*e2d(is51) + Syn.c3(is51).*e1q(is51);
    Syn.Iq(is51) = Syn.Iq(is51) - Syn.c2(is51).*e2d(is51) + Syn.c1(is51).*e1q(is51);
  end
  if ~isempty(is52)
    Syn.Id(is52) = Syn.Id(is52) + Syn.c1(is52).*e2d(is52) + Syn.c3(is52).*e2q(is52);
    Syn.Iq(is52) = Syn.Iq(is52) - Syn.c2(is52).*e2d(is52) + Syn.c1(is52).*e2q(is52);
  end
  if ~isempty(is53)
    Syn.Id(is53) = (e1q(is53)-psid(is53))./(xd(is53)+xl(is53));
    Syn.Iq(is53) = -psiq(is53)./(xq(is53)+xl(is53));
  end
  if ~isempty(is6)
    Syn.Id(is6) = Syn.Id(is6) + Syn.c1(is6).*e2d(is6) + Syn.c3(is6).*e2q(is6);
    Syn.Iq(is6) = Syn.Iq(is6) - Syn.c2(is6).*e2d(is6) + Syn.c1(is6).*e2q(is6);
  end
  if ~isempty(is8)
    Syn.Id(is8) = (e2q(is8)-psid(is8))./xd2(is8);
    Syn.Iq(is8) = (-e2d(is8)-psiq(is8))./xq2(is8);
  end
  Syn.Pg = -DAE.V(Syn.bus).*(Syn.Id.*ss+Syn.Iq.*cc);
  Syn.Qg = -DAE.V(Syn.bus).*(Syn.Id.*cc-Syn.Iq.*ss);

  DAE.gp = DAE.gp + sparse(Syn.bus,1,Syn.Pg,Bus.n,1);
  DAE.gq = DAE.gq + sparse(Syn.bus,1,Syn.Qg,Bus.n,1);

 case 2 % Jacobians of active & reactive powers
  
  M1 = DAE.V(Syn.bus).*(Syn.c1.*cc-Syn.c3.*ss);
  M2 = -DAE.V(Syn.bus).*(Syn.c2.*cc+Syn.c1.*ss);
  M3 = -(Syn.c1.*ss+Syn.c3.*cc);
  M4 = Syn.c2.*ss-Syn.c1.*cc;

  Syn.J11 = DAE.V(Syn.bus).*((Syn.Id-M2).*cc-(M1+Syn.Iq).*ss);
  Syn.J12 = -Syn.Id.*ss-Syn.Iq.*cc-DAE.V(Syn.bus).*(M3.*ss+M4.*cc);
  Syn.J21 = DAE.V(Syn.bus).*((M2-Syn.Id).*ss-(M1+Syn.Iq).*cc);
  Syn.J22 = -Syn.Id.*cc+Syn.Iq.*ss-DAE.V(Syn.bus).*(M3.*cc-M4.*ss);

  DAE.J11 = DAE.J11 + sparse(Syn.bus,Syn.bus,Syn.J11,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(Syn.bus,Syn.bus,Syn.J12,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(Syn.bus,Syn.bus,Syn.J21,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(Syn.bus,Syn.bus,Syn.J22,Bus.n,Bus.n);

 case 3 % Differential equations

  % updating Vf and Pm 
  if Exc.n, Syn.vf(Exc.syn) = DAE.x(Exc.vf); end
  Vf = Syn.vf + Kw.*(omega-1) + Kp.*(Syn.Pg-Syn.Pg0);
  setpm(Tg);

  DAE.f(Syn.delta) = Settings.rad*(omega-1);
  DAE.f(Syn.omega) = (Syn.pm+Syn.Pg-ra.*(Syn.Id.^2+Syn.Iq.^2)-D.*(omega-1)).*iM;
  
  % Model III
  if (~isempty(is3))
    DAE.f(Syn.e1q(is3)) = -a34.*Sat(1,e1q,is3) - a35.*Syn.Id(is3) + a34.*Vf(is3);
  end
  % Model IV
  if (~isempty(is4))
    DAE.f(Syn.e1q(is4)) = -a44.*Sat(1,e1q,is4) - a45.*Syn.Id(is4) + a44.*Vf(is4);
    DAE.f(Syn.e1d(is4)) = -b43.*e1d(is4) + b44.*Syn.Iq(is4);
  end
  % Model V Type 1
  if (~isempty(is51))
    DAE.f(Syn.e1q(is51)) = -a514.*Sat(1,e1q,is51) - a515.*Syn.Id(is51) + a514.*Vf(is51);
    DAE.f(Syn.e1d(is51)) = -b513.*e1d(is51) + b514.*Syn.Iq(is51);
    DAE.f(Syn.e2d(is51)) = -b511.*e2d(is51) + b511.*e1d(is51) + b512.*Syn.Iq(is51);
  end
  % Model V Type 2
  if (~isempty(is52))
    DAE.f(Syn.e1q(is52)) = -a524.*Sat(1,e1q,is52) - a525.*Syn.Id(is52) + a526.*Vf(is52);
    DAE.f(Syn.e2q(is52)) = -a521.*e2q(is52) + a521.*e1q(is52) - ...
        a522.*Syn.Id(is52) + a523.*Vf(is52);
    DAE.f(Syn.e2d(is52)) = -b521.*e2d(is52) + b522.*Syn.Iq(is52);
  end
  % Model V Type 3
  if (~isempty(is53))
    DAE.f(Syn.psiq(is53)) = ...
        Settings.rad.*(DAE.V(bs53).*cc(is53) + ...
                       ra(is53).*Syn.Iq(is53) - omega(is53).*psid(is53));
    DAE.f(Syn.psid(is53)) = ...
        Settings.rad.*(DAE.V(bs53).*ss(is53) + ...
                       ra(is53).*Syn.Id(is53) + omega(is53).*psiq(is53));
    DAE.f(Syn.e1q(is53)) = (a534.*(Syn.vf(is53)-e1q(is53))- ...
                            a531.*DAE.f(Syn.psid(is53))).*a532;
  end
  % Model VI
  if (~isempty(is6))
    DAE.f(Syn.e1q(is6)) = -a4.*Sat(1,e1q,is6) - a5.*Syn.Id(is6) + a6.*Vf(is6);
    DAE.f(Syn.e1d(is6)) = -b3.*e1d(is6) + b4.*Syn.Iq(is6);
    DAE.f(Syn.e2q(is6)) = -a1.*e2q(is6) + a1.*e1q(is6) - a2.*Syn.Id(is6) + a3.*Vf(is6);
    DAE.f(Syn.e2d(is6)) = -b1.*e2d(is6) + b1.*e1d(is6) + b2.*Syn.Iq(is6);
  end
  % Model VIII
  if (~isempty(is8))
    DAE.f(Syn.e1q(is8)) = -a48.*Sat(1,e1q,is8) - a58.*Syn.Id(is8) + a68.*Vf(is8);
    DAE.f(Syn.e1d(is8)) = -b38.*e1d(is8) + b48.*Syn.Iq(is8);
    DAE.f(Syn.e2q(is8)) = -a18.*e2q(is8) + a18.*e1q(is8) - a28.*Syn.Id(is8) + a38.*Vf(is8);
    DAE.f(Syn.e2d(is8)) = -b18.*e2d(is8) + b18.*e1d(is8) + b28.*Syn.Iq(is8);
    DAE.f(Syn.psiq(is8)) = Settings.rad*(DAE.V(bs8).*cc(is8) + ...
                                         ra(is8).*Syn.Iq(is8) - omega(is8).*psid(is8));
    DAE.f(Syn.psid(is8)) = Settings.rad*(DAE.V(bs8).*ss(is8) + ...
                                         ra(is8).*Syn.Id(is8) + omega(is8).*psiq(is8));
  end

 case 4 % Jacobians of differential equations & state variables

  M1 = DAE.V(Syn.bus).*(Syn.c1.*cc-Syn.c3.*ss);
  M2 = -DAE.V(Syn.bus).*(Syn.c2.*cc+Syn.c1.*ss);
  M3 = -(Syn.c1.*ss+Syn.c3.*cc);
  M4 = Syn.c2.*ss-Syn.c1.*cc;

  Wn = Settings.rad;

  % common Jacobians