fcall.m

电力系统分析计算程序

function fcall(p)

global DAE 

if ~p.n, return, end

mltc = DAE.x(p.m);

h = p.con(:,7);
k = p.con(:,8);
mmax = p.con(:,9);
mmin = p.con(:,10);
mstep = p.con(:,11);
ref = p.con(:,12);

error = zeros(p.n,1);
ty1 = find(p.con(:,16) == 1 | p.con(:,16) == 3);
ty2 = find(p.con(:,16) == 2);
if ty1
  error(ty1) = DAE.y(p.vr(ty1))-ref(ty1);
end
if ty2
  Vf = p.u.*DAE.y(p.v1).*exp(i*DAE.y(p.bus1));
  Vt = p.u.*DAE.y(p.v2).*exp(i*DAE.y(p.bus2));  
  y = admittance(p);
  q2 = imag(Vt.*conj((Vt.*mltc-Vf).*y.*mltc));
  error(ty2) = -ref(ty2)-q2(ty2);
end
errel = error./ref;

DAE.f(p.m) = p.u.*(p.con(:,11) <= 0).*(k.*error-h.*mltc);
relay = (errel > mstep) - (errel < -mstep);

% delay tap change for NR stability
p.iters = p.iters + abs(relay);

% delay tap changer during TDS
if DAE.t > 0
  % fixed delay of 0.5 s
  idx1 = find((DAE.t-p.delay) > 0.5);
  idx2 = find((DAE.t-p.delay) <= 0.5);
  if ~isempty(idx1) % let state variables move
    p.delay(idx2) = DAE.t;
  end
  if ~isempty(idx2) % freeze state variables  
    relay(idx2) = 0;
  end
end

idx = find(p.iters == 2);
if ~isempty(idx)
  DAE.x(p.m(idx)) = DAE.x(p.m(idx))+relay(idx).*p.con(idx,11).*p.u(idx);
  p.iters(idx) = 0;
end

% non-windup limits
idx = find(mltc >= mmax & DAE.f(p.m) > 0);
if idx, DAE.f(p.m(idx)) = 0; end
idx = find(mltc <= mmin & DAE.f(p.m) < 0);
if idx, DAE.f(p.m(idx)) = 0; end
DAE.x(p.m) = min(DAE.x(p.m),mmax);
DAE.x(p.m) = max(DAE.x(p.m),mmin);