fxcall.m

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

function Fxcall(a)

global DAE Bus Pod

if ~a.n, return, end

bus1 = a.bus1;
bus2 = a.bus2;
V1 = DAE.V(bus1);
V2 = DAE.V(bus2);
a1 = DAE.a(bus1);
a2 = DAE.a(bus2);

vp = DAE.x(a.vp);
vq = DAE.x(a.vq);
iq = DAE.x(a.iq);

Kr = a.con(:,7);
Tr = a.con(:,8);

ss = sin(a1-a2+a.gamma);
cc = cos(a1-a2+a.gamma);

c1 = a.y.*V1;
vs = sqrt(vp.*vp+vq.*vq);
c2 = a.y.*vs;
U = max(V1.^2+V2.^2-2.*V1.*V2.*cos(a1-a2),1e-6);

P1r =  c1.*V2.*ss;
P2r = -c1.*V2.*ss;
Q1r =  c1.*V1.*cos(a.gamma);
Q2r = -c1.*V2.*cc;

P1a =  c2.*V2.*cc;
P2a = -c2.*V2.*cc;
Q1a = -c2.*V1.*sin(a.gamma);  
Q2a =  c2.*V2.*ss;  

P1vp = P1r.*vp./(V1.*vs) - P1a.*vq./vs.^2;
P2vp = P2r.*vp./(V1.*vs) - P2a.*vq./vs.^2;
Q1vp = Q1r.*vp./(V1.*vs) - Q1a.*vq./vs.^2;
Q2vp = Q2r.*vp./(V1.*vs) - Q2a.*vq./vs.^2;  

P1vq = P1r.*vq./(V1.*vs) + P1a.*vp./vs.^2;
P2vq = P2r.*vq./(V1.*vs) + P2a.*vp./vs.^2;
Q1vq = Q1r.*vq./(V1.*vs) + Q1a.*vp./vs.^2;
Q2vq = Q2r.*vq./(V1.*vs) + Q2a.*vp./vs.^2;

ty1 = a.con(:,2) == 1;
ty2 = a.con(:,2) == 2;

tyvp = a.con(:,15);
tyvq = a.con(:,16);
tyiq = a.con(:,17);

DAE.Fx = DAE.Fx + sparse(a.vp,a.vp,-1./Tr,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(a.vq,a.vq,-1./Tr,DAE.n,DAE.n);
DAE.Fx = DAE.Fx + sparse(a.iq,a.iq,-1./Tr,DAE.n,DAE.n);

ap = (vp <= a.con(:,9)  & vp >= a.con(:,10));
aq = (vq <= a.con(:,11) & vq >= a.con(:,12));
ac = (iq <= a.con(:,13) & iq >= a.con(:,14));

if a.pod
  z = Pod.u(Pod.upfc)./Tr(a.pod);
  DAE.Fx(a.vp(a.pod),Pod.Vs(Pod.upfc)) = z.*tyvp(a.pod).*ap(a.pod);
  DAE.Fx(a.vq(a.pod),Pod.Vs(Pod.upfc)) = z.*tyvq(a.pod).*aq(a.pod);
  DAE.Fx(a.iq(a.pod),Pod.Vs(Pod.upfc)) = z.*Kr(a.pod).*tyiq(a.pod).*ac(a.pod);
end

b = find(ap);
if b
  DAE.Gx = DAE.Gx + sparse(bus1(b),a.vp(b),P1vp(b),2*Bus.n,DAE.n);  
  DAE.Gx = DAE.Gx + sparse(bus2(b),a.vp(b),P2vp(b),2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(bus1(b)+Bus.n,a.vp(b),Q1vp(b),2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(bus2(b)+Bus.n,a.vp(b),Q2vp(b),2*Bus.n,DAE.n);
end

b = find(aq);
if b
  c = find(aq.*ty2);
  if c
    Factor(c) = a.Cp(c)./(1-a.Cp(c));
    F1(c) = Factor(c).*V1(c).*V2(c).*sin(a1(c)-a2(c))./U(c).^.5;
    F2(c)=-F1(c);
    F3(c) = Factor(c).*(V1(c)-V2(c).*cos(a1(c)-a2(c)))./U(c).^.5;
    F4(c) = Factor(c).*(V2(c)-V1(c).*cos(a1(c)-a2(c)))./U(c).^.5;     
    DAE.Fy = DAE.Fy + sparse(a.vq(c),bus1(c), F1(c)./Tr(c),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(a.vq(c),bus2(c), F2(c)./Tr(c),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(a.vq(c),bus1(c)+Bus.n, F3(c)./Tr(c),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(a.vq(c),bus2(c)+Bus.n, F4(c)./Tr(c),DAE.n,2*Bus.n);
  end
  DAE.Gx = DAE.Gx + sparse(bus1(b),a.vq(b),P1vq(b),2*Bus.n,DAE.n);  
  DAE.Gx = DAE.Gx + sparse(bus2(b),a.vq(b),P2vq(b),2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(bus1(b)+Bus.n,a.vq(b),Q1vq(b),2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(bus2(b)+Bus.n,a.vq(b),Q2vq(b),2*Bus.n,DAE.n); 
end

b = find(ac);
if b
  DAE.Gx = DAE.Gx - sparse(Bus.n+bus1(b),a.iq(b),V1(b),2*Bus.n,DAE.n);
  DAE.Fy = DAE.Fy - sparse(a.iq(b),Bus.n+bus1(b),Kr(b)./Tr(b),DAE.n,2*Bus.n);
end