fm_dfig.m

这是一个很适合研究和学习用的电力系统仿真软件

function fm_dfig(flag)
% FM_DFIG define Constant Speed Wind Turbine
%
% FM_DFIG(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 and FM_CSWT
%
%Author:    Federico Milano
%Date:      30-Nov-2003
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global Bus Dfig Wind DAE Settings Varname PV SW PQ

omega_m = DAE.x(Dfig.omega_m);
theta_p = DAE.x(Dfig.theta_p);
idr = DAE.x(Dfig.idr);
iqr = DAE.x(Dfig.iqr);
vw = DAE.x(Wind.vw(Dfig.wind));

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

V = DAE.V(Dfig.bus);
t = DAE.a(Dfig.bus);
st = sin(t);
ct = cos(t);

rs = Dfig.con(:,6);
xs = Dfig.con(:,7);
rr = Dfig.con(:,8);
xr = Dfig.con(:,9);
xm = Dfig.con(:,10);

i2Hm = Dfig.dat(:,3);
Kp = Dfig.con(:,12);
Tp = Dfig.con(:,13);
Kv = Dfig.con(:,14);
Te = Dfig.con(:,15);
R = Dfig.dat(:,4);
A = Dfig.dat(:,5);

a = rs.^2+Dfig.dat(:,1).^2;
a13 = rs./a;
a23 = Dfig.dat(:,1)./a;
a33 = Dfig.dat(:,2);

vds = -V.*st;
vqs =  V.*ct;

ids =  -a13.*(vds-xm.*iqr)-a23.*(vqs+xm.*idr);
iqs =   a23.*(vds-xm.*iqr)-a13.*(vqs+xm.*idr);

vdr = -rr.*idr+(1-omega_m).*(a33.*iqr+xm.*iqs);
vqr = -rr.*iqr-(1-omega_m).*(a33.*idr+xm.*ids);

switch flag
 case 0

  check = 1;

  % Constants
  % xs + xm
  Dfig.dat(:,1) =  Dfig.con(:,7) + Dfig.con(:,10);
  % xr + xm
  Dfig.dat(:,2) =  Dfig.con(:,9) + Dfig.con(:,10);
  % 1/(2*Hm)
  Dfig.dat(:,3) = 1./(2*Dfig.con(:,11));
  % etaGB*4*R*pi*f/p
  Dfig.dat(:,4) = 2*Settings.rad*Dfig.con(:,16)./ ...
      Dfig.con(:,17).*Dfig.con(:,19);
  % A
  Dfig.dat(:,5) = pi*Dfig.con(:,16).*Dfig.con(:,16);
  % Vref
  Dfig.dat(:,6) = V;
  % iqr max & min
  Dfig.dat(:,8) = -Dfig.con(:,21).*Dfig.con(:,10)./Dfig.dat(:,1);
  Dfig.dat(:,9) = -Dfig.con(:,20).*Dfig.con(:,10)./Dfig.dat(:,1);
  % idr max & min
  Dfig.dat(:,10) = -(1+Dfig.con(:,23).*Dfig.con(:,10))./Dfig.dat(:,1);
  Dfig.dat(:,11) = -(1+Dfig.con(:,22).*Dfig.con(:,10))./Dfig.dat(:,1);

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

  vds = -Vc.*sin(ac);
  vqs =  Vc.*cos(ac);

  for i = 1:Dfig.n

    % parameters
    Vds = vds(i);
    Vqs = vqs(i);
    Rs = rs(i);
    Rr = rr(i);
    Xm = xm(i);
    x1 = Dfig.dat(i,1);
    x2 = Dfig.dat(i,2);
    Pg = Pc(i);
    Qg = Qc(i);

    % Initial guess
    Pb = Settings.mva*Pg/Dfig.con(i,3);
    Pb = max(min(Pb,1),0);
    omega = (Pb+1)/2;
    slip = 1-omega;

    iqr = -x1*Dfig.con(i,3)*(2*omega-1)/Vc(i)/Xm/Settings.mva/omega;
    A = [-Rs  x1; Vqs  -Vds];
    B = [Vds-Xm*iqr; Qg];
    Is = A\B;
    ids = Is(1);
    iqs = Is(2);
    idr = -(Vqs+Rs*iqs+x1*ids)/Xm;
    vdr = -Rr*idr+slip*(x2*iqr+Xm*iqs);
    vqr = -Rr*iqr-slip*(x2*idr+Xm*ids);

    jac = zeros(7,7);
    eqn = zeros(7,1);
    inc = ones(7,1);

    x = zeros(7,1);
    x(1) = ids;
    x(2) = iqs;
    x(3) = idr;
    x(4) = iqr;
    x(5) = vdr;
    x(6) = vqr;
    x(7) = 0.51;

    jac(1,1) = -Rs;
    jac(1,2) =  x1;
    jac(1,4) =  Xm;
    jac(2,1) = -x1;
    jac(2,2) = -Rs;
    jac(2,3) = -Xm;
    jac(3,3) = -Rr;
    jac(3,5) = -1;
    jac(4,4) = -Rr;
    jac(4,6) = -1;
    jac(5,1) =  Vds;
    jac(5,2) =  Vqs;
    %jac(6,1) =  Vqs;
    %jac(6,2) = -Vds;
    jac(6,3) = -Xm*Vc(i)/x1;

    k = x1*Dfig.con(i,3)/Vc(i)/Xm/Settings.mva;

    iter = 0;

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

      if iter > 20
        fm_disp([' * * Initialization of doubly fed ind. gen. #', ...
                 num2str(i),' failed.'])
        check = 0;
        break
      end

      eqn(1) = -Rs*x(1)+x1*x(2)+Xm*x(4)-Vds;
      eqn(2) = -Rs*x(2)-x1*x(1)-Xm*x(3)-Vqs;
      eqn(3) = -Rr*x(3)+(1-x(7))*(x2*x(4)+Xm*x(2))-x(5);
      eqn(4) = -Rr*x(4)-(1-x(7))*(x2*x(3)+Xm*x(1))-x(6);
      eqn(5) =  Vds*x(1)+Vqs*x(2)+x(5)*x(3)+x(6)*x(4)-Pg;
      eqn(6) = -Xm*Vc(i)*x(3)/x1-Vc(i)*Vc(i)/x1 - Qg;
      %eqn(6) =  Vqs*x(1)-Vds*x(2)-Qg;
      eqn(7) = -k*(2*x(7)-1)-x(4)*x(7);

      jac(3,2) = (1-x(7))*Xm;
      jac(3,4) = (1-x(7))*x2;
      jac(3,7) = -(x2*x(4)+Xm*x(2));
      jac(4,1) = -(1-x(7))*Xm;
      jac(4,3) = -(1-x(7))*x2;
      jac(4,7) = x2*x(3)+Xm*x(1);
      jac(5,3) = x(5);
      jac(5,4) = x(6);
      jac(5,5) = x(3);
      jac(5,6) = x(4);
      jac(7,4) = -x(7);
      jac(7,7) = -x(4)-2*k;

      inc = -jac\eqn;
      x = x + inc;
      iter = iter + 1;

    end

    ids = x(1);
    iqs = x(2);
    idr = x(3);
    iqr = x(4);
    vdr = x(5);
    vqr = x(6);
    omega = x(7);

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

    % wind turbine state variables
    DAE.x(Dfig.idr(i)) = idr;
    DAE.x(Dfig.iqr(i)) = iqr;
    DAE.x(Dfig.omega_m(i)) = omega;
    DAE.x(Dfig.theta_p(i)) = theta;
    % Vref
    KV = Kv(i);
    if KV == 0 % no voltage control
      Dfig.dat(i,6) = 0;
    else
      Dfig.dat(i,6) = Vc(i)-(idr+Vc(i)/Xm)/Kv(i);
    end
    % iqr offset (~= 0 if Pg = 1 p.u.)
    Dfig.dat(i,7) = max(2*k*(omega-1)/omega,0);

    % electrical torque
    Tel = Xm*(iqr*ids-idr*iqs);
    if Pg < 0
      fm_disp([' * * Turbine power is negative at bus <',Varname.bus{Dfig.bus(i)},'>.'])
      fm_disp(['     Wind speed <',num2str(Dfig.wind(i)),'> cannot be initilized.'])
      DAE.x(Wind.vw(Dfig.wind(i))) = 1;
      continue
    end
    % wind power [MW]
    Pw = Tel*omega*Settings.mva*1e6;
    % wind speed
    iter = 0;
    incvw = 1;
    eqnvw = 1;
    R = Dfig.dat(i,4);
    A = Dfig.dat(i,5);
    % initial guess for wind speed
    vw = 0.9*Wind.con(Dfig.wind(i),2);
    while abs(eqnvw) > 1e-7
      if iter > 50
        wspeed = num2str(Dfig.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,theta,1)-Pw;
      jacvw = windpower(rho(i),vw,A,R,omega,theta,2);
      incvw = -eqnvw/jacvw(2);
      vw = vw + incvw;
      iter = iter + 1;
    end
    % average initial wind speed [p.u.]
    DAE.x(Wind.vw(Dfig.wind(i))) = vw/Wind.con(Dfig.wind(i),2);
    % find & delete static generators
    if ~fm_rmgen(Dfig.bus(i)), check = 0; end
  end

  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;
  q = -V.*(Dfig.con(:,10).*idr+V)./Dfig.dat(:,1);

  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