fcall.m

电力系统分析计算程序

function p = fcall(p)

global DAE

if ~p.n, return, end

V = DAE.y(p.vbus);
Vref = DAE.y(p.vref);

if ~isempty(p.ty1)
  bcv = DAE.x(p.bcv);
  Tr = p.con(p.ty1,6);
  Kr = p.con(p.ty1,7);
  bcv_max = p.u(p.ty1).*p.con(p.ty1,9);
  bcv_min = p.u(p.ty1).*p.con(p.ty1,10);
  DAE.f(p.bcv) = p.u(p.ty1).*(Kr.*(Vref(p.ty1)-V(p.ty1))-bcv)./Tr;
  a = find(bcv >= bcv_max & DAE.f(p.bcv) > 0);
  if ~isempty(a), DAE.f(p.bcv(a)) = 0; end
  a = find(bcv <= bcv_min & DAE.f(p.bcv) < 0);
  if ~isempty(a), DAE.f(p.bcv(a)) = 0; end
  DAE.x(p.bcv) = min(bcv_max,DAE.x(p.bcv));
  DAE.x(p.bcv) = max(bcv_min,DAE.x(p.bcv));
end

if ~isempty(p.ty2)
  alpha = DAE.x(p.alpha);
  vm = DAE.x(p.vm);
  a_max = p.u(p.ty2).*p.con(p.ty2,9);
  a_min = p.u(p.ty2).*p.con(p.ty2,10);
  T2 = p.con(p.ty2,6);
  K = p.con(p.ty2,7);
  Kd = p.con(p.ty2,11);
  T1 = p.con(p.ty2,12);
  Km = p.con(p.ty2,13);
  Tm = p.con(p.ty2,14);
  xl = p.con(p.ty2,15);
  xc = p.con(p.ty2,16);
  DAE.f(p.vm) = p.u(p.ty2).*(Km.*V(p.ty2)-vm)./Tm;
  DAE.f(p.alpha) = p.u(p.ty2).*(-Kd./T2.*alpha + ...
                                K.*T1./T2./Tm.*(vm-Km.*V(p.ty2)) ...
                                + K./T2.*(Vref(p.ty2)-vm));
  a = find(alpha >= a_max & DAE.f(p.alpha) > 0);
  if ~isempty(a), DAE.f(p.alpha(a)) = 0; end
  a = find(alpha <= a_min & DAE.f(p.alpha) < 0);
  if ~isempty(a), DAE.f(p.alpha(a)) = 0; end
  DAE.x(p.alpha) = min(a_max,DAE.x(p.alpha));
  DAE.x(p.alpha) = max(a_min,DAE.x(p.alpha));
end

p.Be = bsvc(p);