fxcall.m

来自「基于PSAT 软件的多目标最优潮流计算用于中小型电力系统的分析和管理」· M 代码 · 共 69 行

function Fxcall(a)

global DAE Bus Pod Cluster

if ~a.n, return, end

V = DAE.V(a.bus);

if ~isempty(a.ty1)
  bcv = DAE.x(a.bcv);
  Tr = a.con(a.ty1,6);
  Kr = a.con(a.ty1,7);
  bcv_max = a.con(a.ty1,9);
  bcv_min = a.con(a.ty1,10);
  DAE.Fx = DAE.Fx + sparse(a.bcv,a.bcv,-1./Tr,DAE.n,DAE.n);
  u = (bcv <= bcv_max & bcv >= bcv_min);
  DAE.Gx = DAE.Gx - sparse(Bus.n+a.bus(a.ty1),a.bcv,u.*V(a.ty1).^2,2*Bus.n,DAE.n);
  DAE.Fy = DAE.Fy - sparse(a.bcv,Bus.n+a.bus(a.ty1),u.*Kr./Tr,DAE.n,2*Bus.n);
  z = a.ty1(find(u));
  if z & a.pod
    h = a.pod(z);
    k = Pod.svc(z);
    DAE.Fx(a.bcv(h),Pod.Vs(k)) = Pod.u(k).*Kr(h)./Tr(h);
  end
  if z & a.cluster
    h = a.cluster(z);
    k = Cluster.svc(z);
    DAE.Fx(a.bcv(h),Cluster.Vs(k)) = Cluster.u(k).*Kr(h)./Tr(h);
    DAE.Fx(Cluster.Vs(k),a.bcv(h)) = -Cluster.u(k).*Cluster.dVsdQ(k).*V(z).^2;
  end
end

if ~isempty(a.ty2)
  alpha = DAE.x(a.alpha);
  vm = DAE.x(a.vm);
  a_max = a.con(a.ty2,9);
  a_min = a.con(a.ty2,10);
  T2 = a.con(a.ty2,6);
  K = a.con(a.ty2,7);
  Kd = a.con(a.ty2,11);
  T1 = a.con(a.ty2,12);
  Km = a.con(a.ty2,13);
  Tm = a.con(a.ty2,14);
  xl = a.con(a.ty2,15);
  xc = a.con(a.ty2,16);
  DAE.Fx = DAE.Fx + sparse(a.vm,a.vm,-1./Tm,DAE.n,DAE.n);
  DAE.Fy = DAE.Fy + sparse(a.vm,Bus.n+a.bus(a.ty2),Km./Tm,DAE.n,2*Bus.n);
  DAE.Fx = DAE.Fx + sparse(a.alpha,a.alpha,-Kd./T2,DAE.n,DAE.n);
  u = (alpha < a_max & alpha > a_min);
  k1 = -u.*K.*T1./T2./Tm.*Km;
  k2 = u.*K.*T1./T2./Tm - K./T2;
  k3 = -2*u.*V(a.ty2).*V(a.ty2).*(1-cos(2*alpha))./xl/pi;
  DAE.Fy = DAE.Fy + sparse(a.alpha,a.bus(a.ty2)+Bus.n,k1,DAE.n,2*Bus.n);
  DAE.Fx = DAE.Fx + sparse(a.alpha,a.vm,k2,DAE.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(Bus.n+a.bus(a.ty2),a.alpha,k3,2*Bus.n,DAE.n);
  z = a.ty2(find(u));
  if z & a.pod
    h = a.pod(z);
    k = Pod.svc(z);
    DAE.Fx(a.alpha(h),Pod.Vs(k)) = Pod.u(k).*K(h)./T2(h);
  end
  if z & a.cluster
    h = a.cluster(z);
    k = Cluster.svc(z);
    DAE.Fx(a.alpha(h),Cluster.Vs(k)) = Cluster.u(k).*K(h)./T2(h);
    DAE.Fx(Cluster.Vs(k),a.alpha(h)) = -Cluster.u(k).*k3(z).*Cluster.dVsdQ(k);
  end
end