fm_dfig.m

电力系统的psat

  if ~check
    fm_disp('Doubly fed induction generator cannot be properly initialized.')
  else
    fm_disp(['Doubly fed induction generators initialized.'])
  end

 case 1

  p = vds.*ids+vqs.*iqs+vdr.*idr+vqr.*iqr;
  q = vqs.*ids-vds.*iqs;

  DAE.gp = DAE.gp - sparse(Dfig.bus,1,p,Bus.n,1);
  DAE.gq = DAE.gq - sparse(Dfig.bus,1,q,Bus.n,1);

 case 2

  pv = -ids.*st + iqs.*ct;
  pt = -ids.*vqs + iqs.*vds;
  qv =  ids.*ct + iqs.*st;
  qt =  ids.*vds + iqs.*vqs;

  idsv = a13.*st - a23.*ct;
  idst = a13.*vqs - a23.*vds;
  iqsv = -a23.*st - a13.*ct;
  iqst = -a23.*vqs - a13.*vds;

  k = (1-omega_m).*xm;

  vdrv = k.*iqsv;
  vdrt = k.*iqst;
  vqrv = -k.*idsv;
  vqrt = -k.*idst;

  j11 = vds.*idst + vqs.*iqst + vdrt.*idr + vqrt.*iqr + pt;
  j12 = vds.*idsv + vqs.*iqsv + vdrv.*idr + vqrv.*iqr + pv;
  j21 = vqs.*idst - vds.*iqst + qt;
  j22 = vqs.*idsv - vds.*iqsv + qv;

  DAE.J11 = DAE.J11 - sparse(Dfig.bus,Dfig.bus,j11,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 - sparse(Dfig.bus,Dfig.bus,j12,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 - sparse(Dfig.bus,Dfig.bus,j21,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 - sparse(Dfig.bus,Dfig.bus,j22,Bus.n,Bus.n);

 case 3

  % wind speed in m/s
  Vw = vw.*Wind.con(Dfig.wind,2);
  Pw = windpower(rho,Vw,A,R,omega_m,theta_p,1)/Settings.mva/1e6;
  % mechanical torque
  Tm = Pw./omega_m;

  % motion equation
  % ---------------
  DAE.f(Dfig.omega_m) = (Tm-xm.*(iqr.*ids-idr.*iqs)).*i2Hm;

  % speed control equations
  % -----------------------
  Pm = Dfig.con(:,3).*max(min(2*omega_m-1,1),0)/Settings.mva;
  DAE.f(Dfig.iqr) = (-(xs+xm).*Pm./V./xm./omega_m-iqr-Dfig.dat(:,7))./Te;
  % non-windup limiter
  iqr_max = -Dfig.con(:,21);
  iqr_min = -Dfig.con(:,20);
  idx = find(iqr <= iqr_min & DAE.f(Dfig.iqr) < 0);
  if idx, DAE.f(Dfig.iqr(idx)) = 0; end
  DAE.x(Dfig.iqr) = max(DAE.x(Dfig.iqr),iqr_min);
  idx = find(iqr >= iqr_max & DAE.f(Dfig.iqr) > 0);
  if idx, DAE.f(Dfig.iqr(idx)) = 0; end
  DAE.x(Dfig.iqr) = min(DAE.x(Dfig.iqr),iqr_max);

  % voltage control equation
  % ------------------------
  DAE.f(Dfig.idr) = Kv.*(V-Dfig.dat(:,6))-V./xm-idr;
  % non-windup limiter
  idr_max = -Dfig.con(:,23);
  idr_min = -Dfig.con(:,22);
  idx = find(idr <= idr_min & DAE.f(Dfig.idr) < 0);
  if idx, DAE.f(Dfig.idr(idx)) = 0; end
  DAE.x(Dfig.idr) = max(DAE.x(Dfig.idr),idr_min);
  idx = find(idr >= idr_max & DAE.f(Dfig.idr) > 0);
  if idx, DAE.f(Dfig.idr(idx)) = 0; end
  DAE.x(Dfig.idr) = min(DAE.x(Dfig.idr),idr_max);

  % pitch control equation
  % ----------------------
  % vary the pitch angle only by steps of 1% of the fn
  phi = round(1000*(omega_m-1))/1000;
  DAE.f(Dfig.theta_p) = (Kp.*phi-theta_p)./Tp;
  % non-windup limiter
  idx = find(theta_p <= 0 & DAE.f(Dfig.theta_p) < 0);
  if idx, DAE.f(Dfig.theta_p(idx)) = 0; end
  DAE.x(Dfig.theta_p) = max(DAE.x(Dfig.theta_p),0);

 case 4

  Vw = vw.*Wind.con(Dfig.wind,2);
  dPwdx = windpower(rho,Vw,A,R,omega_m,theta_p,2)./Settings.mva/1e6;
  Pw = windpower(rho,Vw,A,R,omega_m,theta_p,1)/Settings.mva/1e6;
  Tm = Pw./omega_m;
  Tsp = Dfig.con(:,3).*min(2*omega_m-1,1)./omega_m/Settings.mva;
  Tsp = max(Tsp,0);

  slip = 1-omega_m;
  iqr_min = -Dfig.con(:,3)/Settings.mva;

  % d f / d y
  % -----------

  idsv =  a13.*st  - a23.*ct;
  idst =  a13.*vqs - a23.*vds;
  iqsv = -a23.*st  - a13.*ct;
  iqst = -a23.*vqs - a13.*vds;

  ot = xm.*(idr.*iqst-iqr.*idst).*i2Hm;
  ov = xm.*(idr.*iqsv-iqr.*idsv).*i2Hm;
  iqrv = (xs+xm).*Tsp./Te./V./V./xm;

  DAE.Fy = DAE.Fy + sparse(Dfig.omega_m,Dfig.bus,ot,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Dfig.omega_m,Dfig.bus+Bus.n,ov,DAE.n,2*Bus.n);

  % d g / d x
  % -----------

  idsidr = -a23.*xm;
  idsiqr =  a13.*xm;
  iqsidr = -a13.*xm;
  iqsiqr = -a23.*xm;

  vdridr = -rr + slip.*xm.*iqsidr;
  vdriqr = slip.*(a33 + xm.*iqsiqr);
  vqriqr = -rr - slip.*xm.*idsiqr;
  vqridr = -slip.*(a33 + xm.*idsidr);

  vdrom = -(a33.*iqr+xm.*iqs);
  vqrom = a33.*idr+xm.*ids;

  pidr = vds.*idsidr + vqs.*iqsidr + idr.*vdridr + vdr + iqr.*vqridr;
  piqr = vds.*idsiqr + vqs.*iqsiqr + idr.*vdriqr + vqr + iqr.*vqriqr;
  pom = idr.*vdrom + iqr.*vqrom;
  qidr = vqs.*idsidr - vds.*iqsidr;
  qiqr = vqs.*idsiqr - vds.*iqsiqr;

  DAE.Gx = DAE.Gx - sparse(Dfig.bus,Dfig.omega_m,pom,2*Bus.n,DAE.n);

  % d f / d x
  % -----------

  oidr = xm.*(idr.*iqsidr+iqs-iqr.*idsidr).*i2Hm;
  oiqr = xm.*(idr.*iqsiqr-ids-iqr.*idsiqr).*i2Hm;

  % mechanical equation
  % -------------------
  DAE.Fx = DAE.Fx + sparse(Dfig.omega_m,Dfig.omega_m, ...
                           (dPwdx(:,1)-Tm).*i2Hm./omega_m,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Dfig.omega_m,Wind.vw(Dfig.wind), ...
                           dPwdx(:,2).*i2Hm./omega_m,DAE.n,DAE.n);

  % pitch angle control equation
  % ----------------------------
  DAE.Fx = DAE.Fx + sparse(Dfig.theta_p,Dfig.theta_p,-1./Tp,DAE.n,DAE.n);
  idx = find(theta_p > 0 & DAE.f(Dfig.theta_p) ~= 0);
  if ~isempty(idx)
    DAE.Fx = DAE.Fx + sparse(Dfig.theta_p(idx),Dfig.omega_m(idx), ...
                             -Kp(idx)./Tp(idx),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Dfig.omega_m(idx),Dfig.theta_p(idx), ...
                             dPwdx(idx,3).*i2Hm(idx)./omega_m(idx),DAE.n,DAE.n);
  end

  % speed control equation
  % ----------------------
  kiqr = -(xs+xm)./V./xm./Te;
  tspo = 2*Dfig.con(:,3)/Settings.mva;
  idx = find(Tsp == 0 & (2*omega_m-1) >= 1);
  if ~isempty(idx), tspo(idx) = 0; end
  iqr_max = -Dfig.con(:,21);
  iqr_min = -Dfig.con(:,20);

  idx = find(iqr >= iqr_min & iqr <= iqr_max & DAE.f(Dfig.iqr) ~= 0);
  if ~isempty(idx)
    ki = Dfig.iqr(idx);
    ko = Dfig.omega_m(idx);
    kb = Dfig.bus(idx);

    DAE.Fx = DAE.Fx + sparse(ki,ko,kiqr(idx).* ...
                             (tspo(idx)-Tsp(idx))./omega_m(idx),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(ko,ki,oiqr(idx),DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(ki,kb+Bus.n,iqrv(idx),DAE.n,2*Bus.n);
    DAE.Gx = DAE.Gx - sparse(kb+Bus.n,ki,qiqr(idx),2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx - sparse(kb,ki,piqr(idx),2*Bus.n,DAE.n);
  end

  DAE.Fx = DAE.Fx + sparse(Dfig.iqr,Dfig.iqr,-1./Te,DAE.n,DAE.n);

  % voltage control equation
  % ------------------------
  idr_max = -Dfig.con(:,23);
  idr_min = -Dfig.con(:,22);
  idx = find(idr >= idr_min & idr <= idr_max & DAE.f(Dfig.idr) ~= 0);
  if ~isempty(idx)
    ki = Dfig.idr(idx);
    ko = Dfig.omega_m(idx);
    kb = Dfig.bus(idx);

    DAE.Fx = DAE.Fx - sparse(ko,ki,oidr(idx),DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(ki,kb+Bus.n,Kv(idx)-1./xm(idx),DAE.n,2*Bus.n);
    DAE.Gx = DAE.Gx - sparse(kb,ki,pidr(idx),2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx - sparse(kb+Bus.n,ki,qidr(idx),2*Bus.n,DAE.n);
  end
  DAE.Fx = DAE.Fx + sparse(Dfig.idr,Dfig.idr,-1,DAE.n,DAE.n);

 case 5 % non-windup limiters

  idx = find(theta_p <= 0 & DAE.f(Dfig.theta_p) <= 0);
  if ~isempty(idx)
    k = Dfig.theta_p(idx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    if Settings.octave
      DAE.Ac(k,k) = eye(length(idx));
    else
      DAE.Ac(k,k) = speye(length(idx));
    end
  end

  iqr_max = -Dfig.con(:,21);
  iqr_min = -Dfig.con(:,20);
  idx = find((iqr <= iqr_min | iqr >= iqr_max) & DAE.f(Dfig.iqr) == 0);
  if ~isempty(idx)
    k = Dfig.iqr(idx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    if Settings.octave
      DAE.Ac(k,k) = eye(length(idx));
    else
      DAE.Ac(k,k) = speye(length(idx));
    end
  end

  idr_max = -Dfig.con(:,23);
  idr_min = -Dfig.con(:,22);
  idx = find((idr <= idr_min | idr >= idr_max) & DAE.f(Dfig.idr) == 0);
  if ~isempty(idx)
    k = Dfig.idr(idx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    if Settings.octave
      DAE.Ac(k,k) = eye(length(idx));
    else
      DAE.Ac(k,k) = speye(length(idx));
    end
  end

end

%--------------------------------------------------------
function output = windpower(rho,vw,Ar,R,omega,theta,type)
%--------------------------------------------------------

lambda = omega.*R./vw;
lambdai = 1./(1./(lambda+0.08*theta)-0.035./(theta.^3+1));

switch type
 case 1 % Pw

  cp = 0.22*(116./lambdai-0.4*theta-5).*exp(-12.5./lambdai);
  output = 0.5*rho.*cp.*Ar.*vw.^3;

 case 2 % d Pw / d x

  output = zeros(length(omega),3);

  a1 = exp(-12.5./lambdai);
  a2 = (lambda+0.08*theta).^2;
  a3 = 116./lambdai-0.4*theta-5;
  a4 = -9.28./(lambda+0.08*theta).^2 + ...
       12.180*theta.*theta./(theta.^3+1).^2-0.4;
  a5 = 1.000./(lambda+0.08*theta).^2 - ...
       1.3125*theta.*theta./(theta.^3+1).^2;

  % d Pw / d omega_m
  output(:,1) = R.*a1.*rho.*vw.*vw.*Ar.*(-12.760+1.3750*a3)./a2;

  % d Pw / d vw
  output(:,2) = (omega.*R.*(12.760-1.3750*a3)./a2 ...
      + 0.330*a3.*vw).*vw.*Ar.*rho.*a1;

  % d Pw / d theta_p
  output(:,3) = 0.110*rho.*(a4 + a3.*a5).*a1.*Ar.*vw.^3;

end