fxcall.m

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

function Fxcall(p)

global DAE Bus Pod

if ~p.n, return, end

bus1 = p.bus1;
bus2 = p.bus2;
V1 = DAE.V(bus1);
V2 = DAE.V(bus2);
a1 = DAE.a(bus1);
a2 = DAE.a(bus2);
ss = sin(a1-a2);
cc = cos(a1-a2);
den = ssscden(p);
kp = p.Cp./(1-p.Cp);
vcs = DAE.x(p.vcs);
vpi = DAE.x(p.vpi);
Tr = p.con(:,7);
vcs_max = p.con(:,8);
vcs_min = p.con(:,9);

Ms = p.y./den;

P1vs = Ms.*V1.*V2.*ss;  
P2vs = -P1vs;
Q1vs = Ms.*V1.*(V1-V2.*cc);
Q2vs = Ms.*V2.*(V2-V1.*cc);

DAE.Fx = DAE.Fx - sparse(p.vcs,p.vcs,1./Tr,DAE.n,DAE.n);

u = (vcs <= vcs_max & vcs >= vcs_min);
ty2 = p.con(:,2) == 2;
ty3 = p.con(:,2) == 3;

if p.pod
  DAE.Fx(p.vcs(p.pod),Pod.Vs(Pod.sssc)) = Pod.u(Pod.sssc).*u(p.pod)./Tr(p.pod);
end

DAE.Gx = DAE.Gx + sparse(bus1,p.vcs,u.*P1vs,2*Bus.n,DAE.n);  
DAE.Gx = DAE.Gx + sparse(bus2,p.vcs,u.*P2vs,2*Bus.n,DAE.n);
DAE.Gx = DAE.Gx + sparse(bus1+Bus.n,p.vcs,u.*Q1vs,2*Bus.n,DAE.n);
DAE.Gx = DAE.Gx + sparse(bus2+Bus.n,p.vcs,u.*Q2vs,2*Bus.n,DAE.n);    

a = find(u.*ty2);
if a
  k = kp(a)./den(a)./Tr(a);
  F1 = k.*(V1(a).*V2(a).*ss(a));   
  F3 = k.*(V1(a)-V2(a).*cc(a));
  F4 = k.*(V2(a)-V1(a).*cc(a));     
  DAE.Fy = DAE.Fy + sparse(p.vcs(a),bus1(a),F1,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy - sparse(p.vcs(a),bus2(a),F1,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(p.vcs(a),bus1(a)+Bus.n,F3,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(p.vcs(a),bus2(a)+Bus.n,F4,DAE.n,2*Bus.n);
end

a = find(u.*ty3);
if a
  
  Kpr = p.con(a,11);
  Kin = p.con(a,12);
  c1 = p.y(a).*(1+vcs(a)./den(a));
  c2 = p.y(a).*vcs(a)./(den(a).^3);    
  Ms = p.y(a)./den(a);
  M2 = c2.*V1(a).*V2(a).*ss(a);
  M1 = -M2;
  M3 = c2.*(V2(a).*cc(a)-V1(a));
  M4 = c2.*(V1(a).*cc(a)-V2(a));              
  Jps1 = V1(a).*V2(a).*(M1.*ss(a)+c1.*cc(a));
  Jps2 = V1(a).*V2(a).*(M2.*ss(a)-c1.*cc(a));
  Jps3 = V2(a).*ss(a).*(M3.*V1(a)+c1);
  Jps4 = V1(a).*ss(a).*(M4.*V2(a)+c1);
  Jps5 = Ms.*V1(a).*V2(a).*ss(a);
  
  DAE.Fx = DAE.Fx - sparse(p.vcs(a),p.vcs(a),Jps5.*Kpr(a)./Tr(a),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(p.vcs(a),p.vpi(a),1./Tr(a),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(p.vpi(a),p.vcs(a),Jps5.*Kin(a),DAE.n,DAE.n);
  
  DAE.Fy = DAE.Fy - sparse(p.vcs(a),bus1(a),Jps1.*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy - sparse(p.vcs(a),bus2(a),Jps2.*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy - sparse(p.vcs(a),bus1(a)+Bus.n,Jps3.*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy - sparse(p.vcs(a),bus2(a)+Bus.n,Jps4.*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
  
  DAE.Fy = DAE.Fy - sparse(p.vpi(a),bus1(a),Jps1.*Kin(a),DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy - sparse(p.vpi(a),bus2(a),Jps2.*Kin(a),DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy - sparse(p.vpi(a),bus1(a)+Bus.n,Jps3.*Kin(a),DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy - sparse(p.vpi(a),bus2(a)+Bus.n,Jps4.*Kin(a),DAE.n,2*Bus.n);
  
  c = find(p.con(a,10) == 2);      
  if c
    a = a(c);
    
    DAE.Fy = DAE.Fy + sparse(p.vcs(a),bus1(a),Jps1(c).*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(p.vcs(a),bus2(a),Jps2(c).*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(p.vcs(a),bus1(a)+Bus.n,Jps3(c).*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(p.vcs(a),bus2(a)+Bus.n,Jps4(c).*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
    
    DAE.Fy = DAE.Fy + sparse(p.vpi(a),bus1(a),Jps1(c).*Kin(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(p.vpi(a),bus2(a),Jps2(c).*Kin(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(p.vpi(a),bus1(a)+Bus.n,Jps3(c).*Kin(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(p.vpi(a),bus2(a)+Bus.n,Jps4(c).*Kin(a),DAE.n,2*Bus.n);

  end
  
end