numjacs.m

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

fm_call('i')

ffn  = DAE.f;
xa   = DAE.x;
anga = DAE.a;
Va   = DAE.V;   
tol  = NLA.tol;

A = zeros(DAE.n,Bus.n);
B = zeros(DAE.n,Bus.n);
Gx = zeros(2*Bus.n,DAE.n);
Fx = zeros(DAE.n,DAE.n);

fm_call('i')

for j=1:Bus.n;
  deltaa = max(tol,abs(tol*DAE.a(j)));
  DAE.a(j)=DAE.a(j)+deltaa;
  
  fm_call('i')

  A(:,j)=(DAE.f-ffn)./deltaa;

  DAE.x  = xa;
  DAE.a  = anga;
  DAE.V  = Va;   

end

fm_call('i')

for j=1:Bus.n;
  deltaV = max(tol,abs(tol*DAE.V(j)));
  DAE.V(j)=DAE.V(j)+deltaV;
  
  fm_call('i')

  B(:,j)=(DAE.f-ffn)./deltaV;

  DAE.x  = xa;
  DAE.a  = anga;
  DAE.V  = Va;   

end

Fy = [A,B];

fm_call('i')

ggn = DAE.g;

A = zeros(2*Bus.n,Bus.n);
B = zeros(2*Bus.n,Bus.n);

for j=1:Bus.n;
  deltaa = max(tol,abs(tol*DAE.a(j)));
  DAE.a(j)=DAE.a(j)+deltaa;
  
  fm_call('i')

  A(:,j)=(DAE.g-ggn)./deltaa;

  DAE.x  = xa;
  DAE.a  = anga;
  DAE.V  = Va;   

end

fm_call('i')

for j=1:Bus.n;
  deltaV = max(tol,abs(tol*DAE.V(j)));
  DAE.V(j)=DAE.V(j)+deltaV;
  
  fm_call('i')

  B(:,j)=(DAE.g-ggn)./deltaV;

  DAE.x  = xa;
  DAE.a  = anga;
  DAE.V  = Va;   

end

Jlfv = [A,B];

fm_call('i')

for j=1:DAE.n;
  deltaX = max(tol,abs(tol*DAE.x(j)));
  DAE.x(j)=DAE.x(j)+deltaX;
  
  fm_call('i')

  Gx(:,j)=(DAE.g-ggn)./deltaX;

  DAE.x  = xa;
  DAE.a  = anga;
  DAE.V  = Va;   
  fm_call('i')

end

fm_call('i')

for j=1:DAE.n;
  deltaX = max(tol,abs(tol*DAE.x(j)));
  DAE.x(j)=DAE.x(j)+deltaX;
  if Oxl.n
    idx = find(Oxl.v == j);
    if idx, DAE.x(j)=DAE.x(j)-deltaX; end
  end
  
  fm_call('i')

  Fx(:,j)=(DAE.f-ffn)./deltaX;

  DAE.x  = xa;
  DAE.a  = anga;
  DAE.V  = Va;   
  fm_call('i')

end

Gbus = Bus.n+[SW.bus;PV.bus];
Jlfv(Gbus,:) = 0;
Jlfv(:,Gbus) = 0;
Jlfv(Gbus,Gbus) = eye(PV.n+SW.n);
Jlfv(:,SW.bus) = 0;
Jlfv(SW.bus,:) = 0;
Jlfv(SW.bus,SW.bus) = eye(SW.n);
Fy(:,SW.bus) = 0;
Gx(SW.bus,:) = 0;
Fy(:,Gbus) = 0;
Gx(Gbus,:) = 0;