fm_ddsg.m

电力系统的psat

  s1 = sin(th);
  t1 = s1./c1;
  vds = -rs.*ids+omega_m.*xq.*iqs;

  DAE.J22 = DAE.J22 - sparse(Ddsg.bus,Ddsg.bus,idc./c1,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 - sparse(Ddsg.bus,Ddsg.bus,t1.*V.*idc./c1+(1+t1.^2).*vds.*ids,Bus.n,Bus.n);

 case 3

  % wind speed in m/s
  % -----------------
  Vw = vw.*Wind.con(Ddsg.wind,2);

  % mechanical torque
  % -----------------
  Pw = windpower(rho,Vw,A,R,omega_m,theta_p,1)/Settings.mva/1e6;
  Tm = Pw./omega_m;

  % motion equation
  % ---------------
  DAE.f(Ddsg.omega_m) = (0.5.*(Tm-(-xd.*ids+psip).*iqs-xq.*iqs.*ids))./Hm;

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

  % voltage control equation
  % ------------------------
  DAE.f(Ddsg.idc) = (-idc+Kv.*(Vref-V))./Tv;
  idx = find(idc >= idc_max & DAE.f(Ddsg.idc) > 0);
  if idx, DAE.f(Ddsg.idc(idx)) = 0; end
  DAE.x(Ddsg.idc) = min(idc,idc_max);
  idx = find(idc <= idc_min & DAE.f(Ddsg.idc) < 0);
  if idx, DAE.f(Ddsg.idc(idx)) = 0; end
  DAE.x(Ddsg.idc) = max(idc,idc_min);

  % generator reactive power control equations
  % ------------------------------------------
  DAE.f(Ddsg.ids) = (psip./xd-sqrt(psip.*psip./xd./xd-Qref./omega_m./xd)-ids)./Teq;
  idx = find(ids >= ids_max & DAE.f(Ddsg.ids) > 0);
  if idx, DAE.f(Ddsg.ids(idx)) = 0; end
  DAE.x(Ddsg.ids) = min(ids,ids_max);
  idx = find(ids <= ids_min & DAE.f(Ddsg.ids) < 0);
  if idx, DAE.f(Ddsg.ids(idx)) = 0; end
  DAE.x(Ddsg.ids) = max(ids,ids_min);

  % speed control equations
  % -----------------------
  Pm = Ddsg.con(:,3).*max(min(2*omega_m-1,1),0)/Settings.mva;
  DAE.f(Ddsg.iqs) = (Pm./(psip-xd.*ids)./omega_m-iqs)./Tep;
  idx = find(iqs >= iqs_max & DAE.f(Ddsg.iqs) > 0);
  if idx, DAE.f(Ddsg.iqs(idx)) = 0; end
  DAE.x(Ddsg.iqs) = min(iqs,iqs_max);
  idx = find(iqs <= iqs_min & DAE.f(Ddsg.iqs) < 0);
  if idx, DAE.f(Ddsg.iqs(idx)) = 0; end
  DAE.x(Ddsg.iqs) = max(iqs,iqs_min);

 case 4

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

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

  DAE.Fy = DAE.Fy + sparse(Ddsg.idc,Ddsg.bus+Bus.n,-Kv./Tv,DAE.n,2*Bus.n);

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

  w = omega_m;
  t1 = tan(th);
  c1 = cos(th);
  k1 = psip-xd.*ids;
  k2 = xq.*iqs;
  k3 = xd.*iqs;
  k4 = xq.*ids;
  k5 = rs.*ids;
  k6 = rs.*iqs;

  DAE.Gx = DAE.Gx - sparse(Ddsg.bus,Ddsg.omega_m,k2.*ids+k1.*iqs,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx - sparse(Ddsg.bus+Bus.n,Ddsg.omega_m,t1.*k2.*ids+k1.*iqs,2*Bus.n,DAE.n);

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

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


  % voltage control
  % ---------------
  DAE.Fx = DAE.Fx + sparse(Ddsg.idc,Ddsg.idc,-1./Tv,DAE.n,DAE.n);
  idx = find(ids >= ids_min & ids <= ids_max & DAE.f(Ddsg.ids) ~= 0);
  if ~isempty(idx)
    DAE.Gx = DAE.Gx - sparse(Ddsg.bus(idx)+Bus.n, ...
                             Ddsg.idc(idx),V(idx)./c1(idx),2*Bus.n,DAE.n);
  end

  % reactive power control
  % ----------------------
  DAE.Fx = DAE.Fx + sparse(Ddsg.ids,Ddsg.ids,-1./Teq,DAE.n,DAE.n);

  idx = find(ids >= ids_min & ids <= ids_max & DAE.f(Ddsg.ids) ~= 0);
  if ~isempty(idx)

    i1 = Ddsg.ids(idx);
    i2 = Ddsg.omega_m(idx);
    i3 = Ddsg.bus(idx);

    o1 = (-0.5./((psip(idx)./xd(idx)).^2-Qref(idx)./w(idx)./xd(idx)).^0.5 ...
          .*Qref(idx)./w(idx).^2./xd(idx))./Teq(idx);
    o2 = (xd(idx)-xq(idx)).*iqs(idx).*i2Hm(idx);
    o3 = -2.*k5(idx)+w(idx).*k2(idx)-w(idx).*k3(idx);
    o4 = -t1(idx).*k5(idx)+t1(idx).*(-k5(idx)+w(idx).*k2(idx))-w(idx).*k3(idx);

    DAE.Fx = DAE.Fx + sparse(i1,i2,o1,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(i2,i1,o2,DAE.n,DAE.n);
    DAE.Gx = DAE.Gx - sparse(i3,i1,o3,2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx - sparse(i3+Bus.n,i1,o4,2*Bus.n,DAE.n);

  end

  % speed control
  % -------------
  DAE.Fx = DAE.Fx + sparse(Ddsg.iqs,Ddsg.iqs,-1./Tep,DAE.n,DAE.n);

  tspo = 2*Ddsg.con(:,3)/Settings.mva;
  idx = find(Tsp == 0 & (2*omega_m-1) >= 1);
  if ~isempty(idx), tspo(idx) = 0; end

  idx = find(iqs >= iqs_min & iqs <= iqs_max & DAE.f(Ddsg.iqs) ~= 0);
  if ~isempty(idx)

    i1 = Ddsg.iqs(idx);
    i2 = Ddsg.omega_m(idx);
    i3 = Ddsg.bus(idx);

    o1 = (tspo(idx)-Tsp(idx))./w(idx)./ ...
         (psip(idx)-xd(idx).*ids(idx))./Teq(idx);
    o4 = Pw(idx).*xd(idx)./w(idx)./(psip(idx)-xd(idx).*ids(idx)).^2./Tep(idx);
    o5 = ((xd(idx)-xq(idx)).*ids(idx)-psip(idx)).*i2Hm(idx);
    o6 = w(idx).*k4(idx)-2.*k6(idx)+w(idx).*k1(idx);
    o7 = t1(idx).*w(idx).*k4(idx)-2.*k6(idx)+w(idx).*k1(idx);

    DAE.Fx = DAE.Fx + sparse(i1,i2,o1,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(i1,Ddsg.theta_p(idx),o2,DAE.n,DAE.n);
    %DAE.Fx = DAE.Fx + sparse(i1,Wind.vw(Ddsg.wind(idx)),o3,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(i1,Ddsg.ids(idx),o4,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(i2,i1,o5,DAE.n,DAE.n);
    DAE.Gx = DAE.Gx - sparse(i3,i1,o6,2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx - sparse(i3+Bus.n,i1,o7,2*Bus.n,DAE.n);

  end

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

 case 5 % non-windup limiters

  idx = find(theta_p <= 0 & DAE.f(Ddsg.theta_p) <= 0);
  if ~isempty(idx)
    k = Ddsg.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

  idx = find((idc >= idc_max | idc <= idc_min) & DAE.f(Ddsg.idc) == 0);
  if ~isempty(idx)
    k = Ddsg.idc(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

  idx = find((ids >= ids_max | ids <= ids_min) & DAE.f(Ddsg.ids) == 0);
  if ~isempty(idx)
    k = Ddsg.ids(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

  idx = find((iqs >= iqs_max | iqs <= iqs_min) & DAE.f(Ddsg.iqs) == 0);
  if ~isempty(idx)
    k = Ddsg.iqs(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

% -------------------------------------------------------------------
% function for creating warning messages
function ddsgwarn(idx, msg)
fm_disp(strcat('Warning: DDSG #',int2str(idx),msg))