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

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

  for i = 1:Ddsg.n

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

  end

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