setx0.m

电力系统分析计算程序

function a = setx0(a)

global Bus DAE Settings Wind

if ~a.n, return, end

check = 1;

rho = getrho(Wind,a.wind);
Kv = a.con(:,13);
rs = a.con(:,6);
xd = a.con(:,7);
xq = a.con(:,8);
psip = a.con(:,9);

%check time constants
idx = find(a.con(:,10) == 0);
if idx
  warn(a,idx, 'Inertia Hm cannot be zero. Hm = 2.5 s will be used.'),
  a.con(idx,10) = 2.5;
end

idx = find(a.con(:,12) == 0);
if idx
  warn(a,idx, 'Time constant Tp cannot be zero. Tp = 0.001 s will be used.'),
  a.con(idx,12) = 0.001;
end

idx = find(a.con(:,14) == 0);
if idx
  warn(a,idx, 'Time constant Tv cannot be zero. Tv = 0.001 s will be used.'),
  a.con(idx,14) = 0.001;
end

idx = find(a.con(:,16) == 0);
if idx
  warn(a,idx, 'Time constant Teq cannot be zero. Teq = 0.001 s will be used.'),
  a.con(idx,16) = 0.001;
end

idx = find(a.con(:,15) == 0);
if idx
  warn(a,idx, 'Time constant Tep cannot be zero. Tep = 0.001 s will be used.'),
  a.con(idx,15) = 0.001;
end

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

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

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

for i = 1:a.n
  
  % idc
  DAE.x(a.idc(i)) = cos(ac(i))*(Qc(i)-Pc(i)*tan(ac(i)))/Vc(i);
  % Vref
  a.dat(i,3) =  DAE.x(a.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;
  k = a.con(i,3)/Settings.mva;
  omega_m = k/(2*k-iqs*(psip-xd.*ids));
  
  % theta_p
  theta_p = a.con(i,11)*round(1000*(omega_m-1))/1000;
  theta_p = max(theta_p,0);
  
  % wind turbine state variables and constants
  DAE.y(a.ids(i)) = ids; % ids
  DAE.x(a.iqs(i)) = iqs; % iqs
  DAE.x(a.omega_m(i)) = omega_m;
  DAE.x(a.theta_p(i)) = theta_p;
  
  % electrical torque
  Tel = ((xq(i)-xd(i))*ids+psip(i))*iqs;
  
  % wind power [MW]
  Pw = Tel*omega_m*Settings.mva*1e6;
  if Pc(i) < 0
    fm_disp([' * * Turbine power is negative at bus <',Bus.names{a.bus(i)},'>.'])
    fm_disp(['     Wind speed <',num2str(a.wind(i)),'> cannot be initilized.'])
    DAE.x(getidx(Wind,a.wind(i))) = 1;
    continue
  end
  % wind speed
  iter = 0;
  incvw = 1;
  eqnvw = 1;
  R = a.dat(i,1);
  A = a.dat(i,2);
  % initial guess for wind speed
  vw = 0.9*getvw(Wind,a.wind(i));
  while abs(eqnvw) > 1e-7
    if iter > 50
      wspeed = num2str(a.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(a,rho(i),vw,A,R,omega_m,theta_p,1)-Pw;
    jacvw = windpower(a,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(getidx(Wind,a.wind(i))) = vw/getvw(Wind,a.wind(i));
  % find & delete static generators
  if ~fm_rmgen(a.u(i)*a.bus(i)), check = 0; end
end

DAE.x(a.idc) = a.u.*DAE.x(a.idc);
DAE.y(a.ids) = a.u.*DAE.y(a.ids);
DAE.x(a.iqs) = a.u.*DAE.x(a.iqs);
DAE.x(a.omega_m) = a.u.*DAE.x(a.omega_m);
DAE.x(a.theta_p) = a.u.*DAE.x(a.theta_p);

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