fxcall.m

电力系统分析计算程序

function Fxcall(a)

global DAE Settings

if ~a.n, return, end

ord = a.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  = a.bus(is2);
bs3  = a.bus(is3);
bs4  = a.bus(is4);
bs51 = a.bus(is51);
bs52 = a.bus(is52);
bs53 = a.bus(is53);
bs6  = a.bus(is6);
bs8  = a.bus(is8);

vs2  = a.vbus(is2);
vs3  = a.vbus(is3);
vs4  = a.vbus(is4);
vs51 = a.vbus(is51);
vs52 = a.vbus(is52);
vs53 = a.vbus(is53);
vs6  = a.vbus(is6);
vs8  = a.vbus(is8);

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

ag = DAE.y(a.bus);
vg = a.u.*DAE.y(a.vbus);
ss = sin(delta-ag);
cc = cos(delta-ag);

iM = a.u./a.con(:,18);
D = a.con(:,19);

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

M1 =  vg.*(a.c1.*cc-a.c3.*ss);
M2 = -vg.*(a.c2.*cc+a.c1.*ss);
M3 = -(a.c1.*ss+a.c3.*cc).*a.u;
M4 =  (a.c2.*ss-a.c1.*cc).*a.u;

rad = 2*pi*Settings.freq;
Wn = rad*a.u;

% common Jacobians
DAE.Fx = DAE.Fx - sparse(a.delta,a.delta,(~a.u),DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(a.delta,a.omega,Wn,DAE.n,DAE.n);
DAE.Fx = DAE.Fx - sparse(a.omega,a.omega,iM.*D+(~a.u),DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(a.omega,a.delta,2*ra.*(a.Id.*M1+a.Iq.*M2).*iM,DAE.n,DAE.n);
DAE.Fy = DAE.Fy - sparse(a.omega,a.bus,2*ra.*(a.Id.*M1+a.Iq.*M2).*iM,DAE.n,DAE.m);
DAE.Fy = DAE.Fy - sparse(a.omega,a.vbus,2*ra.*(a.Id.*M3+a.Iq.*M4).*iM,DAE.n,DAE.m);
DAE.Fy = DAE.Fy - sparse(a.omega,a.p,iM,DAE.n,DAE.m);
DAE.Fy = DAE.Fy + sparse(a.omega,a.pm,iM,DAE.n,DAE.m);
DAE.Gx = DAE.Gx + sparse(a.p,a.delta,a.J11,DAE.m,DAE.n);
DAE.Gx = DAE.Gx + sparse(a.q,a.delta,a.J21,DAE.m,DAE.n);

Gp1 = vg.*(a.c3.*ss + a.c1.*cc);
Gp2 = vg.*(a.c1.*ss - a.c2.*cc);
Gq1 = vg.*(a.c3.*cc - a.c1.*ss);
Gq2 = vg.*(a.c1.*cc + a.c2.*ss);

N1 = -2*ra.*(a.Id.*a.c3+a.Iq.*a.c1).*iM;
N2 = -2*ra.*(a.Id.*a.c1-a.Iq.*a.c2).*iM;

% Model II
if ~isempty(is2)
  o2 = a.omega(is2);
  q1 = vg(is2).*(ss(is2).*a.c3(is2)+cc(is2).*a.c1(is2));
  q2 = vg(is2).*(cc(is2).*a.c3(is2)-ss(is2).*a.c1(is2));
  k2 = -2*ra(is2).*iM(is2).*(a.c3(is2).*a.Id(is2)+a.c1(is2).*a.Iq(is2));
  DAE.Fy = DAE.Fy + sparse(o2,a.vf(is2),k2,DAE.n,DAE.m);
  DAE.Fy = DAE.Fy - sparse(o2,a.vf(is2),Kp(is2).*k2,DAE.n,DAE.m);
  DAE.Fx = DAE.Fx + sparse(o2,o2,Kw(is2).*k2,DAE.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.p(is2),o2,Kw(is2).*q1,DAE.m,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.q(is2),o2,Kw(is2).*q2,DAE.m,DAE.n);
end

% Model III
if (~isempty(is3))
  o3 = a.omega(is3);
  e1q3 = a.e1q(is3);
  e1q(is3) = DAE.x(e1q3);
  a34 = a.u(is3)./Td10(is3);
  a35 = a34.*(xd(is3)-xd1(is3));
  DAE.Fx = DAE.Fx + sparse(o3,e1q3,N1(is3),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e1q3,a.delta(is3),a35.*M1(is3),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e1q3,o3,a34.*Kw(is3),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1q3,e1q3,a34.*synsat(a,2,e1q(is3),is3)+a35.*a.c3(is3)+(~a.u(is3)),DAE.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.p(is3),e1q3,Gp1(is3),DAE.m,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.q(is3),e1q3,Gq1(is3),DAE.m,DAE.n);
  DAE.Fy = DAE.Fy + sparse(e1q3,a.vf(is3),a34,DAE.n,DAE.m);
  DAE.Fy = DAE.Fy - sparse(e1q3,bs3,a35.*M1(is3),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy - sparse(e1q3,vs3,a35.*M3(is3),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy - sparse(e1q3,a.p(is3),a34.*Kp(is3),DAE.n,DAE.m);
end

% Model IV
if ~isempty(is4)
  o4 = a.omega(is4);
  e1q4 = a.e1q(is4);
  e1d4 = a.e1d(is4);
  e1d(is4) = DAE.x(e1d4);
  e1q(is4) = DAE.x(e1q4);
  a44 = a.u(is4)./Td10(is4);
  a45 = a44.*(xd(is4)-xd1(is4));
  b43 = a.u(is4)./Tq10(is4);
  b44 = b43.*(xq(is4)-xq1(is4));
  DAE.Fy = DAE.Fy + sparse(e1q4,a.vf(is4),a44,DAE.n,DAE.m);
  DAE.Fx = DAE.Fx + sparse(o4,e1q4,N1(is4),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(o4,e1d4,N2(is4),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e1q4,a.delta(is4),a45.*M1(is4),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e1q4,o4,a44.*Kw(is4),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1q4,e1q4,a44.*synsat(a,2,e1q(is4),is4)+a45.*a.c3(is4)+(~a.u(is4)), DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1q4,e1d4,a45.*a.c1(is4), DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1d4,a.delta(is4),b44.*M2(is4), DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e1d4,e1q4,b44.*a.c1(is4), DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1d4,e1d4,b43+b44.*a.c2(is4)+(~a.u(is4)), DAE.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.p(is4),e1q4,Gp1(is4),DAE.m,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.p(is4),e1d4,Gp2(is4),DAE.m,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.q(is4),e1q4,Gq1(is4),DAE.m,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.q(is4),e1d4,Gq2(is4),DAE.m,DAE.n);
  DAE.Fy = DAE.Fy - sparse(e1q4,bs4,a45.*M1(is4),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy - sparse(e1q4,vs4,a45.*M3(is4),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy + sparse(e1d4,bs4,b44.*M2(is4),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy + sparse(e1d4,vs4,b44.*M4(is4),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy - sparse(e1q4,a.p(is4),a44.*Kp(is4),DAE.n,DAE.m);
end

% Model V Type 1
if(~isempty(is51))
  o51 = a.omega(is51);
  e1q51 = a.e1q(is51);
  e1d51 = a.e1d(is51);
  e2d51 = a.e2d(is51);
  e1d(is51) = DAE.x(e1d51);
  e1q(is51) = DAE.x(e1q51);
  e2d(is51) = DAE.x(e2d51);
  gq = xd1(is51)./xq1(is51).*Tq20(is51)./Tq10(is51).*(xq(is51)-xq1(is51));
  a514 = a.u(is51)./Td10(is51);
  a515 = a514.*(xd(is51)-xd1(is51));
  b511 = a.u(is51)./Tq20(is51);
  b512 = b511.*(xq1(is51)-xd1(is51)+gq);
  b513 = a.u(is51)./Tq10(is51);
  b514 = b513.*(xq(is51)-xq1(is51)-gq);
  DAE.Fy = DAE.Fy + sparse(e1q51,a.vf(is51),a514,DAE.n,DAE.m);
  DAE.Fx = DAE.Fx + sparse(o51,e1q51,N1(is51),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(o51,e2d51,N2(is51),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e1q51,a.delta(is51),a515.*M1(is51),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e1q51,o51,a514.*Kw(is51),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1q51,e1q51,a514.*synsat(a,2,e1q(is51),is51)+a515.*a.c3(is51)+(~a.u(is51)),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1q51,e2d51,a515.*a.c1(is51),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1d51,a.delta(is51),b514.*M2(is51),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e1d51,e1q51,b514.*a.c1(is51),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1d51,e1d51,b513+(~a.u(is51)),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1d51,e2d51,b514.*a.c2(is51),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e2d51,a.delta(is51),b512.*M2(is51),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e2d51,e1q51,b512.*a.c1(is51),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e2d51,e1d51,b511,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e2d51,e2d51,b511+b512.*a.c2(is51)+(~a.u(is51)),DAE.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.p(is51),e1q51,Gp1(is51),DAE.m,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.p(is51),e2d51,Gp2(is51),DAE.m,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.q(is51),e1q51,Gq1(is51),DAE.m,DAE.n);
  DAE.Gx = DAE.Gx + sparse(a.q(is51),e2d51,Gq2(is51),DAE.m,DAE.n);
  DAE.Fy = DAE.Fy - sparse(e1q51,bs51,a515.*M1(is51),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy - sparse(e1q51,vs51,a515.*M3(is51),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy + sparse(e1d51,bs51,b514.*M2(is51),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy + sparse(e1d51,vs51,b514.*M4(is51),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy + sparse(e2d51,bs51,b512.*M2(is51),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy + sparse(e2d51,vs51,b512.*M4(is51),DAE.n,DAE.m);
  DAE.Fy = DAE.Fy - sparse(e1q51,a.p(is51),a514.*Kp(is51),DAE.n,DAE.m);
end

% Model V Type 2
if(~isempty(is52))
  o52 = a.omega(is52);
  e1q52 = a.e1q(is52);
  e2q52 = a.e2q(is52);
  e2d52 = a.e2d(is52);
  e1q(is52) = DAE.x(e1q52);
  e2q(is52) = DAE.x(e2q52);
  e2d(is52) = DAE.x(e2d52);
  Taa = a.con(:,24);
  gd = xd2(is52)./xd1(is52).*Td20(is52)./Td10(is52).*(xd(is52)-xd1(is52));
  a521 = a.u(is52)./Td20(is52);
  a522 = a521.*(xd1(is52)-xd2(is52)+gd);
  a523 = a.u(is52).*Taa(is52)./Td10(is52)./Td20(is52);
  a524 = a.u(is52)./Td10(is52);
  a525 = a524.*(xd(is52)-xd1(is52)-gd);
  a526 = a524.*(1-Taa(is52)./Td10(is52));
  b521 = a.u(is52)./Tq20(is52);
  b522 = b521.*(xq(is52)-xq2(is52));
  DAE.Fx = DAE.Fx + sparse(o52,e2q52,N1(is52),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(o52,e2d52,N2(is52),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e1q52,a.delta(is52),a525.*M1(is52),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e1q52,o52,a526.*Kw(is52),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1q52,e1q52,a524.*synsat(a,2,e1q(is52),is52)+(~a.u(is52)),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1q52,e2q52,a525.*a.c3(is52),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e1q52,e2d52,a525.*a.c1(is52),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e2q52,a.delta(is52),a522.*M1(is52),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e2q52,o52,a523.*Kw(is52),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(e2q52,e1q52,a521,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e2q52,e2q52,a521+a522.*a.c3(is52)+(~a.u(is52)),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(e2q52,e2d52,a522.*a.c1(is52),DAE.n,DAE.n);