fm_cswt.m

一个较好的MATLAB潮流程序

function fm_cswt(flag)
% FM_CSWT define Constant Speed Wind Turbine
%
% FM_CSWT(FLAG)
%       FLAG = 0 initializations
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%
%see also FM_WIND
%
%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 Cswt Wind DAE Settings Varname PV SW

Wn = Settings.rad;

omega_wr = DAE.x(Cswt.omega_wr);
omega_m = DAE.x(Cswt.omega_m);
gamma = DAE.x(Cswt.gamma);
e1r = DAE.x(Cswt.e1r);
e1m = DAE.x(Cswt.e1m);
vw = DAE.x(Wind.vw(Cswt.wind));
rho = Wind.con(Cswt.wind,3);

V = DAE.V(Cswt.bus);
t = DAE.a(Cswt.bus);
st = sin(t);
ct = cos(t);
Vr = -V.*st;
Vm =  V.*ct;

r1 = Cswt.con(:,6);
rr = Cswt.con(:,8);

x0 = Cswt.dat(:,1);
x1 = Cswt.dat(:,2);
T10 = Cswt.dat(:,3);
i2Hwr = Cswt.dat(:,4);
i2Hm = Cswt.dat(:,5);
Ks = Cswt.con(:,13);
R = Cswt.dat(:,6);
A = Cswt.dat(:,7);
B = Cswt.dat(:,8);

a = r1.^2+x1.^2;
a13 = r1./a;
a23 = x1./a;
a33 = x0-x1;

Im = -a23.*(e1r-Vr)+a13.*(e1m-Vm);
Ir =  a13.*(e1r-Vr)+a23.*(e1m-Vm);

switch flag
 case 0

  check = 1;

  % Constants
  % x0
  Cswt.dat(:,1) =  Cswt.con(:,7) + Cswt.con(:,10);
  % x'
  Cswt.dat(:,2) =  Cswt.con(:,7) + Cswt.con(:,10).*Cswt.con(:,9)./ ...
      (Cswt.con(:,10)+Cswt.con(:,9));
  % T'0
  Cswt.dat(:,3) = (Cswt.con(:,10)+Cswt.con(:,9))./Cswt.con(:,8)/Wn;
  % 1/(2*Hwr)
  Cswt.dat(:,4) = 1./(2*Cswt.con(:,11));
  % 1/(2*Hm)
  Cswt.dat(:,5) = 1./(2*Cswt.con(:,12));
  % 4*R*pi*f/p
  Cswt.dat(:,6) = 2*Wn*Cswt.con(:,14)./Cswt.con(:,15).*Cswt.con(:,17);
  % A
  Cswt.dat(:,7) = pi*Cswt.con(:,14).*Cswt.con(:,14);

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

  vr = -Vc.*sin(ac);
  vm =  Vc.*cos(ac);

  for i = 1:Cswt.n

    % parameters
    Vr = vr(i);
    Vm = vm(i);
    Rs = r1(i);
    X0 = Cswt.dat(i,1);
    X1 = Cswt.dat(i,2);
    T10 = Cswt.dat(i,3);
    Pg = Pc(i);

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

    jac(1,1) = Vr;
    jac(1,2) = Vm;
    jac(2,1) = -Rs;
    jac(2,2) = X1;
    jac(2,3) = 1;
    jac(3,1) = -X1;
    jac(3,2) = -Rs;
    jac(3,4) = 1;
    jac(4,2) = (X0-X1)/T10;
    jac(4,3) = -1/T10;
    jac(5,1) = -(X0-X1)/T10;
    jac(5,4) = -1/T10;

    % variables: x(1) = ir
    %            x(2) = im
    %            x(3) = e'r
    %            x(4) = e'm
    %            x(5) = sigma
    % first guess
    x = jac([1:5],[1:4])\[Pg;Vr;Vm;0;0];
    x(5) = 0;

    iter = 0;

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

      if iter > 20
        fm_disp(['Initialization of constant speed wind turbine #',...
                 num2str(i),' failed (convergence problem).'])
        break
      end

      eqn(1) = Vr*x(1) + Vm*x(2) - Pg;
      eqn(2) = x(3) - Vr - Rs*x(1) + X1*x(2);
      eqn(3) = x(4) - Vm - Rs*x(2) - X1*x(1);
      eqn(4) =  Wn*x(5)*x(4) - (x(3)-(X0-X1)*x(2))/T10;
      eqn(5) = -Wn*x(5)*x(3) - (x(4)+(X0-X1)*x(1))/T10;

      jac(4,4) =  Wn*x(5);
      jac(4,5) =  Wn*x(4);
      jac(5,3) = -Wn*x(5);
      jac(5,5) = -Wn*x(3);

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

    end

    Qg = Vm*x(1)-Vr*x(2);
    Te = x(3)*x(1)+x(4)*x(2);

    % shunt capacitor
    Cswt.dat(i,8) = (Qc(i)-Qg)/V(i)/V(i);

    % wind turbine state variables
    DAE.x(Cswt.e1r(i)) = x(3);
    DAE.x(Cswt.e1m(i)) = x(4);
    DAE.x(Cswt.omega_m(i)) = 1-x(5);
    DAE.x(Cswt.gamma(i)) = Te/Cswt.con(i,13);
    DAE.x(Cswt.omega_wr(i)) = DAE.x(Cswt.omega_m(i));

    % wind power
    Pw = Te*DAE.x(Cswt.omega_m(i));
    % wind speed
    iter = 0;
    incvw = 1;
    eqnvw = 1;
    R = Cswt.dat(i,6);
    A = Cswt.dat(i,7);

    vw = 0.9*Wind.con(Cswt.wind(i),2);

    while abs(eqnvw) > 1e-8
      if iter > 20
        fm_disp(['Initialization of wind speed #',...
                 num2str(Cswt.wind(i)),' failed (convergence problem).'])
        break
      end
      eqnvw = windpower(rho(i),vw, ...
                        A,R,DAE.x(Cswt.omega_m(i)),1)-Pw*Settings.mva*1e6;
      jacvw = windpower(rho(i),vw, ...
                        A,R,DAE.x(Cswt.omega_m(i)),2);
      incvw = -eqnvw/jacvw(2);
      %disp([eqnvw+Pw*Cswt.con(i,3),incvw])
      vw = vw + incvw;
      %disp([vw,iter]);
      iter = iter + 1;
    end
    DAE.x(Wind.vw(Cswt.wind(i))) = vw/Wind.con(Cswt.wind(i),2);

  end

  % find and delete PV generators
  check = 1;
  for j = 1:Cswt.n
    idx_pv = [];
    idx_sw = [];
    if PV.n
      idx_pv = find(PV.bus == Cswt.bus(j));
    end
    if SW.n
      idx_sw = find(SW.bus == Cswt.bus(j));
    end
    if isempty(idx_pv) & isempty(idx_sw)
      fm_disp(['Error: No generator found at bus #', ...
	       int2str(Cswt.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: Constant speed wind turbine cannot ', ...
	       'substitute slack buses.'])
      fm_disp(['Constant speed wind turbine will be initialized ', ...
               'but eigenvalue analyses and time domain'])
      fm_disp(['simulations will produce errors ', ...
	       '(singularity in the Jacobian matrix).'])
    end
  end

  % random initial phase angle for shadow tower effect
  %Cswt.dat(:,9) = 2*pi*rand(Cswt.n,1);

  if ~check
    fm_disp('Constant speed wind turbine cannot be properly initialized.')
  else
    fm_disp(['Constant speed wind turbines initialized.'])
  end

 case 1

  DAE.gp = DAE.gp - sparse(Cswt.bus,1,Vr.*Ir+Vm.*Im,Bus.n,1);
  DAE.gq = DAE.gq - sparse(Cswt.bus,1,Vm.*Ir-Vr.*Im+B.*V.*V,Bus.n,1);

 case 2

  Pv = -Ir.*st + Im.*ct;
  Pt = -Ir.*Vm + Im.*Vr;
  Qv =  Ir.*ct + Im.*st;
  Qt =  Ir.*Vr + Im.*Vm;

  IrV =  a13.*st-a23.*ct;
  ImV = -a23.*st-a13.*ct;
  Irt =  a13.*Vm-a23.*Vr;
  Imt = -a23.*Vm-a13.*Vr;

  DAE.J11 = DAE.J11 - sparse(Cswt.bus,Cswt.bus,Pt+Vr.*Irt+Vm.*Imt,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 - sparse(Cswt.bus,Cswt.bus,Pv+Vr.*IrV+Vm.*ImV,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 - sparse(Cswt.bus,Cswt.bus,Qt+Vm.*Irt-Vr.*Imt,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 - sparse(Cswt.bus,Cswt.bus,Qv+Vm.*IrV-Vr.*ImV+2*B.*V,Bus.n,Bus.n);

 case 3

  slip = 1-omega_m;
  % periodic torque pulsation due to "tower shadow" phenomenon
  omega_r = 2*Wn*omega_wr.*Cswt.con(:,17)./Cswt.con(:,15)./Cswt.con(:,16);
  shadow = 0.025*sin(max(DAE.t,0)*omega_r);

  Twr = windpower(rho,vw.*Wind.con(Cswt.wind,2), ...
                  A,R,omega_wr,1)./omega_wr./Settings.mva/1e6;
  DAE.f(Cswt.omega_wr) = (Twr.*(1+shadow)-Ks.*gamma).*i2Hwr;
  DAE.f(Cswt.omega_m) = (Ks.*gamma-e1r.*Ir-e1m.*Im).*i2Hm;
  DAE.f(Cswt.gamma) = Wn*(omega_wr-omega_m);
  DAE.f(Cswt.e1r)  =  Wn*slip.*e1m-(e1r-a33.*Im)./T10;
  DAE.f(Cswt.e1m)  = -Wn*slip.*e1r-(e1m+a33.*Ir)./T10;

 case 4

  dPwdx = windpower(rho,vw.*Wind.con(Cswt.wind,2), ...
                    A,R,omega_wr,2)./Settings.mva/1e6;
  Twr = windpower(rho,vw.*Wind.con(Cswt.wind,2), ...
                  A,R,omega_wr,1)./omega_wr./Settings.mva/1e6;
  slip = 1-omega_m;

  IrV =  a13.*st-a23.*ct;
  ImV = -a23.*st-a13.*ct;
  Irt =  a13.*Vm-a23.*Vr;
  Imt = -a23.*Vm-a13.*Vr;

  DAE.Fx = DAE.Fx + sparse(Cswt.omega_wr,Cswt.gamma,-Ks.*i2Hwr,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Cswt.omega_m,Cswt.gamma,Ks.*i2Hm,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Cswt.gamma,Cswt.omega_wr,Wn,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Cswt.gamma,Cswt.omega_m,-Wn,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx + sparse(Cswt.omega_wr,Cswt.omega_wr,(dPwdx(:,1)-Twr).*i2Hwr./omega_wr,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Cswt.omega_wr,Wind.vw(Cswt.wind),dPwdx(:,2).*i2Hwr./omega_wr,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx + sparse(Cswt.omega_m,Cswt.e1r,(-Ir-e1r.*a13+e1m.*a23).*i2Hm,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Cswt.omega_m,Cswt.e1m,(-Im-e1r.*a23-e1m.*a13).*i2Hm,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx + sparse(Cswt.e1r,Cswt.omega_m,-Wn.*e1m,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Cswt.e1r,Cswt.e1r,-(1+a33.*a23)./T10,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Cswt.e1r,Cswt.e1m,Wn.*slip+a33.*a13./T10,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx + sparse(Cswt.e1m,Cswt.omega_m,Wn.*e1r,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Cswt.e1m,Cswt.e1r,-Wn.*slip-a33.*a13./T10,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Cswt.e1m,Cswt.e1m,-(1+a33.*a23)./T10,DAE.n,DAE.n);

  DAE.Gx = DAE.Gx - sparse(Cswt.bus,Cswt.e1r, a13.*Vr-a23.*Vm,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx - sparse(Cswt.bus,Cswt.e1m, a23.*Vr+a13.*Vm,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx - sparse(Cswt.bus+Bus.n,Cswt.e1r, a23.*Vr+a13.*Vm,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx - sparse(Cswt.bus+Bus.n,Cswt.e1m,-a13.*Vr+a23.*Vm,2*Bus.n,DAE.n);

  DAE.Fy = DAE.Fy + sparse(Cswt.omega_m,Cswt.bus,(-e1r.*Irt-e1m.*Imt).*i2Hm,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Cswt.e1r,Cswt.bus,a33.*Imt./T10,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Cswt.e1m,Cswt.bus,-a33.*Irt./T10,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Cswt.omega_m,Cswt.bus+Bus.n,(-e1r.*IrV-e1m.*ImV).*i2Hm,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Cswt.e1r,Cswt.bus+Bus.n,a33.*ImV./T10,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Cswt.e1m,Cswt.bus+Bus.n,-a33.*IrV./T10,DAE.n,2*Bus.n);

end

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

lambda = omega.*R./vw;
lambdai = 1./(1./lambda+0.002);

switch type
 case 1 % Pw
  cp = 0.44*(125./lambdai-6.94).*exp(-16.5./lambdai);
  output = 0.5*rho.*cp.*Ar.*vw.^3;
 case 2
  output = zeros(length(omega),2);
  a1 = exp(-16.5./lambda-0.0330);
  a2 = 125./lambda-6.690;
  % d Pw / d omega
  output(:,1) = rho.*(-27.5+3.63*a2).*a1.*Ar.*vw.^4./omega./omega./R;
  % d Pw / d vw
  output(:,2) = rho.*((27.5-3.63*a2)./lambda+0.66*a2).*a1.*Ar.*vw.^2;

  % > diff(Pw,vw);
  %                      16.5 vw                3
  %27.50000000 rho exp(- ------- - 0.0330) Ar vw
  %                      omega R
  %----------------------------------------------
  %                   omega R
  %
  %                       /125 vw         \       16.5 vw                3
  %       3.630000000 rho |------- - 6.690| exp(- ------- - 0.0330) Ar vw
  %                       \omega R        /       omega R
  %     - ----------------------------------------------------------------
  %                                   omega R
  %
  %                        /125 vw         \       16.5 vw                2
  %     + 0.6600000000 rho |------- - 6.690| exp(- ------- - 0.0330) Ar vw
  %                        \omega R        /       omega R
  %
  %> diff(Pw,omega);
  %                    4       16.5 vw
  %  27.50000000 rho vw  exp(- ------- - 0.0330) Ar
  %                            omega R
  %- ----------------------------------------------
  %                          2
  %                     omega  R
  %
  %                       /125 vw         \   4       16.5 vw
  %       3.630000000 rho |------- - 6.690| vw  exp(- ------- - 0.0330) Ar
  %                       \omega R        /           omega R
  %     + ----------------------------------------------------------------
  %                                        2
  %                                   omega  R

end