fm_dfig.m

一个较好的MATLAB潮流程序

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-2005 Federico Milano
%
% This toolbox is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2.0 of the License, or
% (at your option) any later version.
%
% This toolbox is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANDABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this toolbox; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
% USA.

global Bus Dfig Wind DAE Settings Varname PV SW

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;

  % 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;

    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) =  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;
      %disp(x')
      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);

    % 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
        fm_disp(['Initialization of wind speed #', ...
                 num2str(Dfig.wind(i)), ...
                 ' failed (convergence problem).'])
       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);

  end

  % find and delete PV generators
  for j = 1:Dfig.n
    idx_pv = [];
    idx_sw = [];
    if PV.n
      idx_pv = find(PV.bus == Dfig.bus(j));
    end
    if SW.n
      idx_sw = find(SW.bus == Dfig.bus(j));
    end
    if isempty(idx_pv) & isempty(idx_sw)
      fm_disp(['Error: No generator found at bus #', ...
	       int2str(Dfig.bus(j))])
      check = 0;
    elseif ~isempty(idx_pv)
      PV.con(idx_pv,:) = [];
      PV.bus(idx_pv) = [];
      PV.n = PV.n - 1;
    else
      fm_disp(['Warning: Doubly fed induction generator cannot ', ...
	       'substitute slack buses.'])
      fm_disp(['Doubly fed induction generator will be initialized ', ...
               'but eigenvalue analyses and time domain'])
      fm_disp(['simulations will produce errors ', ...
	       '(singularity in the Jacobian matrix).'])
    end
  end