fm_ddsg.m

基于PSAT 软件的多目标最优潮流计算用于中小型电力系统的分析和管理

function fm_ddsg(flag)
% FM_DDSG define a direct drive synchronous generator wind turbine
%
% FM_DDSG(FLAG)
%       FLAG = 0 initializations
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%       FLAG = 5 non-windup limiters
%
%see also FM_WIND, FM_CSWT and FM_DFIG
%
%Author:    Federico Milano
%Date:      12-May-2004
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global Bus Ddsg Wind DAE Settings Varname PV SW

omega_m = DAE.x(Ddsg.omega_m);
theta_p = DAE.x(Ddsg.theta_p);
ids = DAE.x(Ddsg.ids);
iqs = DAE.x(Ddsg.iqs);
idc = DAE.x(Ddsg.idc);
vw = DAE.x(Wind.vw(Ddsg.wind));

rho = Wind.con(Ddsg.wind,3);

V = DAE.V(Ddsg.bus);
th = DAE.a(Ddsg.bus);
st = sin(th);
ct = cos(th);

rs = Ddsg.con(:,6);
xd = Ddsg.con(:,7);
xq = Ddsg.con(:,8);
psip = Ddsg.con(:,9);

Hm = Ddsg.con(:,10);
Kp = Ddsg.con(:,11);
Tp = Ddsg.con(:,12);
Kv = Ddsg.con(:,13);
Tv = Ddsg.con(:,14);
Tep = Ddsg.con(:,15);
Teq = Ddsg.con(:,16);
R = Ddsg.dat(:,4);
A = Ddsg.dat(:,5);
Vref = Ddsg.dat(:,1);
Pref = Ddsg.dat(:,6);
Qref = Ddsg.dat(:,7);

ids_max = -Ddsg.con(:,24);
ids_min = -Ddsg.con(:,23);
idc_max = -Ddsg.con(:,24);
idc_min = -Ddsg.con(:,23);
iqs_max = -Ddsg.con(:,22);
iqs_min = -Ddsg.con(:,21);

switch flag
 case 0

  check = 1;

  %check time constants
  idx = find(Hm == 0);
  if idx
    ddsgwarn(idx, 'Inertia Hm cannot be zero. Hm = 2.5 s will be used.'),
  end
  Ddsg.con(idx,5) = 2.5;
  idx = find(Tp == 0);
  if idx
    ddsgwarn(idx, 'Time constant Tp cannot be zero. Tp = 0.001 s will be used.'),
  end
  Ddsg.con(idx,6) = 0.001;
  idx = find(Tv == 0);
  if idx
    ddsgwarn(idx, 'Time constant Tv cannot be zero. Tv = 0.001 s will be used.'),
  end
  Ddsg.con(idx,9) = 0.001;
  idx = find(Teq == 0);
  if idx
    ddsgwarn(idx, 'Time constant Teq cannot be zero. Teq = 0.001 s will be used.'),
  end
  Ddsg.con(idx,15) = 0.001;
  idx = find(Tep == 0);
  if idx
    ddsgwarn(idx, 'Time constant Tep cannot be zero. Tep = 0.001 s will be used.'),
  end
  Ddsg.con(idx,17) = 0.001;

  % Constants
  % etaGB*4*R*pi*f/p
  Ddsg.dat(:,4) = 2*Settings.rad*Ddsg.con(:,17)./ ...
      Ddsg.con(:,18).*Ddsg.con(:,20);
  % A
  Ddsg.dat(:,5) = pi*Ddsg.con(:,17).*Ddsg.con(:,17);

  % Initialization of state variables
  Pc = Bus.Pg(Ddsg.bus);
  Qc = Bus.Qg(Ddsg.bus);
  Vc = DAE.V(Ddsg.bus);
  ac = DAE.a(Ddsg.bus);

  vdc = -Vc.*sin(ac);
  vqc =  Vc.*cos(ac);

  for i = 1:Ddsg.n

    % taking into account PQ generators
    if abs(Pc(i)) < 1e-5 & Bus.Pl(Ddsg.bus(i)) < 0
      Pc(i) = -Bus.Pl(Ddsg.bus(i));
      Qc(i) = -Bus.Ql(Ddsg.bus(i));
    end

    % idc
    DAE.x(Ddsg.idc(i)) = cos(ac(i))*(Qc(i)-Pc(i)*tan(ac(i)))/Vc(i);
    % vref
    Ddsg.dat(i,1) =  DAE.x(Ddsg.idc(i))/Kv(i)+Vc(i);

    jac = zeros(4,4);
    jac(3,1) = -1;
    jac(4,2) = 1;
    jac(3,3) = -rs(i);
    jac(3,4) = xq(i);
    jac(4,3) = xd(i);
    jac(4,4) = rs(i);

    x = zeros(4,1);
    x(2) = 1;
    x(3) = -1/xd(i)-psip(i);

    iter = 0;
    inc = 1;

    while max(abs(inc)) > 1e-8

      if iter > 20
        fm_disp(['Initialization of direct drive syn. gen. #', ...
                 num2str(i),' failed.'])
        check = 0;
        break
      end

      eqn(1,1) = x(1)*x(3)+x(2)*x(4)-Pc(i);
      eqn(2,1) = x(2)*x(3)-x(1)*x(4);
      eqn(3,1) = -x(1)-rs(i)*x(3)+xq(i)*x(4);
      eqn(4,1) = x(2)+rs(i)*x(4)+xd(i)*x(3)-psip(i);

      jac(1,1) = x(3);
      jac(1,2) = x(4);
      jac(1,3) = x(1);
      jac(1,4) = x(2);
      jac(2,1) = -x(4);
      jac(2,2) = x(3);
      jac(2,3) = x(2);
      jac(2,4) = -x(1);

      inc = -jac\eqn;
      x = x + inc;
      %disp(x')
      iter = iter + 1;

    end

    vds = x(1);
    vqs = x(2);
    ids = x(3);
    iqs = x(4);

    omega_m = 1;

    % theta_p
    theta_p = Kp(i)*round(1000*(omega_m-1))/1000;
    theta_p = max(theta_p,0);

    % wind turbine state variables and constants
    DAE.x(Ddsg.ids(i)) = ids; % ids
    Ddsg.dat(i,2) =  ids;     % ids_ref
    DAE.x(Ddsg.iqs(i)) = iqs; % iqs
    Ddsg.dat(i,3) =  iqs;     % iqs_ref
    DAE.x(Ddsg.omega_m(i)) = omega_m;
    DAE.x(Ddsg.theta_p(i)) = theta_p;

    % electrical torque
    Tel = ((xq(i)-xd(i))*ids+psip(i))*iqs;

    % Ps_ref
    Ddsg.dat(i,6) = iqs*(psip(i)-xd(i)*ids);
    % Qs_ref
    Ddsg.dat(i,7) = xd(i)*((psip(i)/xd(i))^2-(ids-psip(i)/xd(i))^2);

    % wind power [MW]
    Pw = Tel*omega_m*Settings.mva*1e6;
    if Pc(i) < 0
      fm_disp([' * * Turbine power is negative at bus <',Varname.bus{Ddsg.bus(i)},'>.'])
      fm_disp(['     Wind speed <',num2str(Ddsg.wind(i)),'> cannot be initilized.'])
      DAE.x(Wind.vw(Ddsg.wind(i))) = 1;
      continue
    end
    % wind speed
    iter = 0;
    incvw = 1;
    eqnvw = 1;
    R = Ddsg.dat(i,4);
    A = Ddsg.dat(i,5);
    % initial guess for wind speed
    vw = 0.9*Wind.con(Ddsg.wind(i),2);
    while abs(eqnvw) > 1e-7
      if iter > 50
        wspeed = num2str(Ddsg.wind(i));
        fm_disp([' * * Initialization of wind speed <', ...
                 wspeed,'> failed (convergence problem).'])
        fm_disp(['     Tip: Try increasing the nominal wind speed <',wspeed,'>.'])
        check = 0;
        break
      end
      eqnvw = windpower(rho(i),vw,A,R,omega_m,theta_p,1)-Pw;
      jacvw = windpower(rho(i),vw,A,R,omega_m,theta_p,2);
      incvw = -eqnvw/jacvw(2);
      vw = vw + incvw;
      iter = iter + 1;
    end
    % average initial wind speed [p.u.]
    DAE.x(Wind.vw(Ddsg.wind(i))) = vw/Wind.con(Ddsg.wind(i),2);
    % find & delete static generators
    if ~fm_rmgen(Ddsg.bus(i)), check = 0; end
  end

  if ~check
    fm_disp(['Direct drive synchronous generator cannot be properly ' ...
             'initialized.'])
  else
    fm_disp(['Direct drive synchronous generators initialized.'])
  end

 case 1

  vds = -rs.*ids+omega_m.*xq.*iqs;
  vqs = -rs.*iqs-omega_m.*(xd.*ids-psip);
  ps = vds.*ids+vqs.*iqs;

  DAE.gp = DAE.gp - sparse(Ddsg.bus,1,ps,Bus.n,1);
  DAE.gq = DAE.gq - sparse(Ddsg.bus,1,V.*idc./cos(th)+tan(th).*ps,Bus.n,1);

 case 2

  c1 = cos(th);
  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;
    DAE.Ac(k,k) = speye(length(idx));
  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;
    DAE.Ac(k,k) = speye(length(idx));
  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;
    DAE.Ac(k,k) = speye(length(idx));
  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;
    DAE.Ac(k,k) = speye(length(idx));
  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))