fm_synit.m

电力系统的psat

function fm_synit
% FM_SYNIT initialize Synchornous Machines state variables,
%         field voltages and mechanical torques
%
% FM_SYNIT
%
%see also FM_SYN
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%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 DAE Syn Bus Varname Settings SW PV jay

if Syn.n == 0; return; end

% machine dynamic orders
ord = Syn.con(:,5);

% indexes of machine orders
is2 = find(ord == 2);
is3 = find(ord == 3);
is4 = find(ord == 4);
is51 = find(ord == 5.1);
is52 = find(ord == 5.2);
is53 = find(ord == 5.3);
is6 = find(ord == 6);
is8 = find(ord == 8);

% indexes of machine buses
bs2 = Syn.bus(is2);
bs3 = Syn.bus(is3);
bs4 = Syn.bus(is4);
bs51 = Syn.bus(is51);
bs52 = Syn.bus(is52);
bs53 = Syn.bus(is53);
bs6 = Syn.bus(is6);
bs8 = Syn.bus(is8);

% constants
Syn.c1 = zeros(Syn.n,1);
Syn.c2 = zeros(Syn.n,1);
Syn.c3 = zeros(Syn.n,1);
Syn.con([is2;is53;is8],20) = 0; % Kw = 0 for these machines
Syn.con([is2;is53;is8],21) = 0; % Kp = 0 for these machines

% check parameters
% -------------------------------------------------------------------
% check inertia M and x'd
idx = find(Syn.con(:,18) <= 0);
if ~isempty(idx)
  Syn.con(idx,18) = 10;
  if Settings.conv
    Syn.con(idx,18) = Syn.con(idx,18).*Syn.con(idx,2)/Settings.mva;
  end
  synwarn(idx,' Inertia cannot be <= 0. M = 10 [kWs/kVa] will be used.')
end
idx = find(Syn.con(:,9) <= 0);
if ~isempty(idx)
  Syn.con(idx,9) = 0.302;
  if Settings.conv
    Syn.con(idx,9)= (Syn.con(idx,9)./Syn.con(idx,2))*Settings.mva;
  end
  synwarn(idx,' x''d cannot be <= 0. x''d = 0.302 [p.u.] will be used.')
end
% check xd, T'd0 and T"d0
is = [is3; is4; is51; is52; is53; is6; is8];
idx = find(Syn.con(is,8) <= 0);
if ~isempty(idx)
  idx = is(idx);
  Syn.con(idx,8) = 1.90;
  if Settings.conv
    Syn.con(idx,8)= (Syn.con(idx,8)./Syn.con(idx,2))*Settings.mva;
  end
  synwarn(idx,' xd cannot be <= 0. xd = 1.90 [p.u.] will be used.')
end
idx = find(Syn.con(is,11) <= 0);
if ~isempty(idx)
  idx = is(idx);
  Syn.con(idx,11) = 8.00;
  synwarn(idx,' T''d0 cannot be <= 0. T''d0 = 8.00 [s] will be used.')
end
% check T"d0, x"d and x"q
is = [is52; is6; is8];
idx = find(Syn.con(is,12) <= 0);
if ~isempty(idx)
  idx = is(idx);
  Syn.con(idx,12) = 0.04;
  synwarn(idx,' T"d0 cannot be <= 0. T"d0 = 0.04 [s] will be used.')
end
idx = find(Syn.con(is,10) <= 0);
if ~isempty(idx)
  idx = is(idx);
  Syn.con(idx,10) = 0.204;
  if Settings.conv
    Syn.con(idx,10)= (Syn.con(idx,10)./Syn.con(idx,2))*Settings.mva;
  end
  synwarn(idx,' x"d cannot be <= 0. x"d = 0.204 [p.u.] will be used.')
end
idx = find(Syn.con(is,15) <= 0);
if ~isempty(idx)
  idx = is(idx);
  Syn.con(idx,15) = 0.30;
  if Settings.conv
    Syn.con(idx,15)= (Syn.con(idx,15)./Syn.con(idx,2))*Settings.mva;
  end
  synwarn(idx,' x"q cannot be <= 0. x"q = 0.30 [p.u.] will be used.')
end
% check T'q0
is = [is4; is51; is6; is8];
idx = find(Syn.con(is,16) <= 0);
if ~isempty(idx)
  idx = is(idx);
  Syn.con(idx,16) = 0.80;
  synwarn(idx,' T''q0 cannot be <= 0. T''q0 = 0.80 [s] will be used.')
end
% check T"q0
is = [is51; is52; is6; is8];
idx = find(Syn.con(is,17) <= 0);
if ~isempty(idx)
  idx = is(idx);
  Syn.con(idx,17) = 0.02;
  synwarn(idx,' T"q0 cannot be <= 0. T"q0 = 0.02 [s] will be used.')
end
% check xq
is = [is3; is4; is51; is52; is53; is6; is8];
idx = find(Syn.con(is,13) <= 0);
if ~isempty(idx)
  idx = is(idx);
  Syn.con(idx,13) = 1.70;
  if Settings.conv
    Syn.con(idx,13)= (Syn.con(idx,13)./Syn.con(idx,2))*Settings.mva;
  end
  synwarn(idx,' xq cannot be <= 0. xq = 1.70 [p.u.] will be used.')
end
% check x'q
is = [is4; is51; is6; is8];
idx = find(Syn.con(is,14) <= 0);
if ~isempty(idx)
  idx = is(idx);
  Syn.con(idx,14) = 0.50;
  if Settings.conv
    Syn.con(idx,14)= (Syn.con(idx,14)./Syn.con(idx,2))*Settings.mva;
  end
  synwarn(idx,' x''q cannot be <= 0. x''q = 0.50 [p.u.] will be used.')
end
% check Taa and saturation factors
% if saturation factors are defined, Taa = 0 is used.
idx = find(Syn.con(:,25) | Syn.con(:,26));
if ~isempty(idx)
  Syn.con(idx,24) = 0;
end

% parameters
xl = -Syn.con(:,6);
ra = Syn.con(:,7);
xd = Syn.con(:,8);
xd1 = Syn.con(:,9);
xd2 = Syn.con(:,10);
xq = Syn.con(:,13);
xq1 = Syn.con(:,14);
xq2 = Syn.con(:,15);
Td10 = Syn.con(:,11);
Tq10 = Syn.con(:,16);
Td20 = Syn.con(:,12);
Tq20 = Syn.con(:,17);
K0 = Syn.con(:,20);
Kp = Syn.con(:,21);
Taa = Syn.con(:,24);

% adjusting parameters for initialization
if ~isempty(is2)
  xq(is2) = xd1(is2);
end

% rotor speeds
DAE.x(Syn.omega) = ones(Syn.n,1);

% check active and reactive power ratios
synbus = sort(Syn.bus);
n_old = -1;
for i = 1:Syn.n
  n_new = synbus(i);
  if n_new ~= n_old
    idx = find(Syn.bus == n_new);
    if length(idx) == 1
      if Syn.con(idx,22) ~= 1
        fm_disp(['Warning: Active power ratio of ', ...
                 'generator #', ...
                 num2str(idx),' must be 1'])
        Syn.con(idx,22) = 1;
      end
      if Syn.con(idx,23) ~= 1
        fm_disp(['Warning: Reactive power ratio of ', ...
                 'generator #', ...
                 num2str(idx),' must be 1'])
        Syn.con(idx,23) = 1;
      end
    elseif length(idx) > 1
      numsyn = length(idx);
      ratiop = sum(Syn.con(idx,22));
      ratioq = sum(Syn.con(idx,23));
      if abs(ratiop-1) > 1e-5
        fm_disp(['Warning: The sum of active power ', ...
                 'ratios of generators #', ...
                 num2str(idx'),' must be 1'])
        Syn.con(idx,22) = 1/numsyn;
      else
        Syn.con(idx(1),22) = Syn.con(idx(1),22)-(ratiop-1);
      end
      if abs(ratioq-1) > 1e-5
        fm_disp(['Warning: The sum of reactive power ', ...
                 'ratios of generators #', ...
                 num2str(idx'),' must be 1'])
        Syn.con(idx,23) = 1/numsyn;
      else
        Syn.con(idx(1),23) = Syn.con(idx(1),23)-(ratioq-1);
      end
    end
    n_old = n_new;
  end
end

% powers and rotor angles
Syn.Pg = -Bus.Pg(Syn.bus).*Syn.con(:,22);
Syn.Qg = -Bus.Qg(Syn.bus).*Syn.con(:,23);
Vg = DAE.V(Syn.bus);
ag = DAE.a(Syn.bus);
V =  Vg.*exp(jay*ag);
S = -Syn.Pg + jay*Syn.Qg;
I = S./conj(V);
delta = angle(V + (ra + jay*(xq+xl)).*I);
DAE.x(Syn.delta) = delta;

% d and q-axis voltages and currents
Vdq = V.*exp(-jay*(delta-pi/2));
Idq = I.*exp(-jay*(delta-pi/2));
vd = real(Vdq);
vq = imag(Vdq);
Syn.Id = real(Idq);
Syn.Iq = imag(Idq);

% mechanical torques/powers
Syn.pm = (vq+ra.*Syn.Iq).*Syn.Iq+(vd+ra.*Syn.Id).*Syn.Id;

% remaining state variables and field voltages
if ~isempty(is2)
  K = 1./(ra(is2).^2+xd1(is2).*xd1(is2)+xd1(is2).*xl(is2)+ ...
          xl(is2).*xd1(is2)+xl(is2).^2);
  Syn.c1(is2) = ra(is2).*K;
  Syn.c2(is2) = (xd1(is2)+xl(is2)).*K;
  Syn.c3(is2) = (xd1(is2)+xl(is2)).*K;
  Syn.vf(is2) = vq(is2)+ra(is2).*Syn.Iq(is2)+ ...
      (xd1(is2)+xl(is2)).*Syn.Id(is2);
end
if ~isempty(is3)
  K = 1./(ra(is3).^2+xq(is3).*xd1(is3)+ ...
          xq(is3).*xl(is3)+xl(is3).*xd1(is3)+xl(is3).^2);
  Syn.c1(is3) = ra(is3).*K;
  Syn.c2(is3) = (xd1(is3)+xl(is3)).*K;
  Syn.c3(is3) = (xq(is3)+xl(is3)).*K;
  DAE.x(Syn.e1q(is3)) = vq(is3)+ra(is3).*Syn.Iq(is3)+ ...
      (xd1(is3)+xl(is3)).*Syn.Id(is3);
  Syn.vf(is3) = Sat(1,DAE.x(Syn.e1q(is3)),is3)+ ...
      (xd(is3)-xd1(is3)).*Syn.Id(is3);
end
if ~isempty(is4)
  K = 1./(ra(is4).^2+xq1(is4).*xd1(is4)+ ...
          xq1(is4).*xl(is4)+xl(is4).*xd1(is4)+xl(is4).^2);
  Syn.c1(is4) = ra(is4).*K;
  Syn.c2(is4) = (xd1(is4)+xl(is4)).*K;
  Syn.c3(is4) = (xq1(is4)+xl(is4)).*K;
  DAE.x(Syn.e1q(is4)) = vq(is4)+ra(is4).*Syn.Iq(is4)+ ...
      (xd1(is4)+xl(is4)).*Syn.Id(is4);
  DAE.x(Syn.e1d(is4)) = vd(is4)+ra(is4).*Syn.Id(is4)- ...
      (xq1(is4)+xl(is4)).*Syn.Iq(is4);
  Syn.vf(is4) = Sat(1,DAE.x(Syn.e1q(is4)),is4)+ ...
      (xd(is4)-xd1(is4)).*Syn.Id(is4);
end
if ~isempty(is51)
  K = 1./(ra(is51).^2+xd1(is51).*xd1(is51)+ ...
          xd1(is51).*xl(is51)+xl(is51).*xd1(is51)+xl(is51).^2);
  Syn.c1(is51) = ra(is51).*K;
  Syn.c2(is51) = (xd1(is51)+xl(is51)).*K;
  Syn.c3(is51) = (xd1(is51)+xl(is51)).*K;
  DAE.x(Syn.e1q(is51)) = vq(is51)+ra(is51).*Syn.Iq(is51)+ ...
      (xd1(is51)+xl(is51)).*Syn.Id(is51);
  DAE.x(Syn.e2d(is51)) = vd(is51)+ra(is51).*Syn.Id(is51)- ...
      (xd1(is51)+xl(is51)).*Syn.Iq(is51);
  DAE.x(Syn.e1d(is51)) = (xq(is51)-xq1(is51)- ...
                          Tq20(is51).*xd1(is51).* ...
                          (xq(is51)-xq1(is51))./Tq10(is51)./ ...
                          xq1(is51)).*Syn.Iq(is51);
  Syn.vf(is51) = Sat(1,DAE.x(Syn.e1q(is51)),is51)+ ...
      (xd(is51)-xd1(is51)).*Syn.Id(is51);
end
if ~isempty(is52)
  K = 1./(ra(is52).^2+(xq2(is52)+xl(is52)).*(xd2(is52)+xl(is52)));
  Syn.c1(is52) = ra(is52).*K;
  Syn.c2(is52) = (xd2(is52)+xl(is52)).*K;
  Syn.c3(is52) = (xq2(is52)+xl(is52)).*K;
  DAE.x(Syn.e2q(is52)) = vq(is52)+ra(is52).*Syn.Iq(is52)+ ...
      (xd2(is52)+xl(is52)).*Syn.Id(is52);
  DAE.x(Syn.e2d(is52)) = vd(is52)+ra(is52).*Syn.Id(is52)- ...
      (xq2(is52)+xl(is52)).*Syn.Iq(is52);
  k1 = xd(is52)-xd1(is52)-Td20(is52).*xd2(is52).* ...
       (xd(is52)-xd1(is52))./Td10(is52)./xd1(is52);
  k2 = xd1(is52)-xd2(is52)+Td20(is52).*xd2(is52).* ...
       (xd(is52)-xd1(is52))./Td10(is52)./xd1(is52);
  %Syn.vf(is52) = (k1+k2).*Syn.Id(is52)+DAE.x(Syn.e2q(is52));
  DAE.x(Syn.e1q(is52)) = -k1.*Taa(is52)./Td10(is52).* ...
      Syn.Id(is52)+(1-Taa(is52)./Td10(is52)).* ...
      (DAE.x(Syn.e2q(is52)) + k2.*Syn.Id(is52));
  Syn.vf(is52) = (k1.*Syn.Id(is52)+Sat(1,DAE.x(Syn.e1q(is52)),is52))./ ...
      (1-Taa(is52)./Td10(is52));
end
if ~isempty(is53)
  DAE.x(Syn.psid(is53)) = ra(is53).*Syn.Iq(is53) + vq(is53);
  DAE.x(Syn.psiq(is53)) = -(ra(is53).*Syn.Id(is53) + vd(is53));
  DAE.x(Syn.e1q(is53)) = DAE.x(Syn.psid(is53)) + ...
      (xd(is53)+xl(is53)).*Syn.Id(is53);
  Syn.vf(is53) = DAE.x(Syn.e1q(is53));
end
if ~isempty(is6)
  K = 1./(ra(is6).^2+xq2(is6).*xd2(is6)+xq2(is6).* ...
          xl(is6)+xl(is6).*xd2(is6)+xl(is6).^2);
  Syn.c1(is6) = ra(is6).*K;
  Syn.c2(is6) = (xd2(is6)+xl(is6)).*K;
  Syn.c3(is6) = (xq2(is6)+xl(is6)).*K;
  DAE.x(Syn.e2q(is6)) = vq(is6)+ra(is6).*Syn.Iq(is6)+ ...
      (xd2(is6)+xl(is6)).*Syn.Id(is6);
  DAE.x(Syn.e2d(is6)) = vd(is6)+ra(is6).*Syn.Id(is6)- ...
      (xq2(is6)+xl(is6)).*Syn.Iq(is6);
  DAE.x(Syn.e1d(is6)) = (xq(is6)-xq1(is6)-Tq20(is6).* ...
                         xq2(is6).*(xq(is6)-xq1(is6))./ ...
                         Tq10(is6)./xq1(is6)).*Syn.Iq(is6);
  k1 = xd(is6)-xd1(is6)-Td20(is6).*xd2(is6).* ...
       (xd(is6)-xd1(is6))./Td10(is6)./xd1(is6);
  k2 = xd1(is6)-xd2(is6)+Td20(is6).*xd2(is6).* ...
       (xd(is6)-xd1(is6))./Td10(is6)./xd1(is6);
  %Syn.vf(is6) = (k1+k2).*Syn.Id(is6)+DAE.x(Syn.e2q(is6));
  DAE.x(Syn.e1q(is6)) = DAE.x(Syn.e2q(is6))+k2.*Syn.Id(is6)- ...
      Taa(is6)./Td10(is6).*((k1+k2).*Syn.Id(is6)+DAE.x(Syn.e2q(is6)));
  Syn.vf(is6) = (k1.*Syn.Id(is6)+Sat(1,DAE.x(Syn.e1q(is6)),is6))./ ...
      (1-Taa(is6)./Td10(is6));
end
if ~isempty(is8)
  DAE.x(Syn.psid(is8)) = ra(is8).*Syn.Iq(is8) + vq(is8);
  DAE.x(Syn.psiq(is8)) = -(ra(is8).*Syn.Id(is8) + vd(is8));
  DAE.x(Syn.e2d(is8)) = -(DAE.x(Syn.psiq(is8)) + ...
                          (xq2(is8)+xl(is8)).*Syn.Iq(is8));
  DAE.x(Syn.e2q(is8)) = DAE.x(Syn.psid(is8)) + ...
      (xd2(is8)+xl(is8)).*Syn.Id(is8);
  DAE.x(Syn.e1d(is8)) = (xq(is8)-xq1(is8)-Tq20(is8).* ...
                         xq2(is8).*(xq(is8)-xq1(is8))./ ...
                         Tq10(is8)./xq1(is8)).*Syn.Iq(is8);
  k1 = xd(is8)-xd1(is8)-Td20(is8).*xd2(is8).* ...
       (xd(is8)-xd1(is8))./Td10(is8)./xd1(is8);
  k2 = xd1(is8)-xd2(is8)+Td20(is8).*xd2(is8).* ...
       (xd(is8)-xd1(is8))./Td10(is8)./xd1(is8);
  %Syn.vf(is8) = (k1+k2).*Syn.Id(is8)+DAE.x(Syn.e2q(is8));
  DAE.x(Syn.e1q(is8)) = -k1.*Taa(is8)./Td10(is8).* ...
      Syn.Id(is8)+(1-Taa(is8)./Td10(is8)).* ...
      (DAE.x(Syn.e2q(is8)) + k2.*Syn.Id(is8));
  Syn.vf(is8) = (k1.*Syn.Id(is8)+Sat(1,DAE.x(Syn.e1q(is8)),is8))./ ...
      (1-Taa(is8)./Td10(is8));
end

% adjusting field voltages in case of power feedback
Syn.vf = Syn.vf+Kp.*(Syn.pm+Syn.Pg);

% check if a PV or Slack generator are connected at the machine buses
for i = 1:Syn.n
  ok = find(SW.bus == Syn.bus(i));
  if isempty(ok)
    ok = find(PV.bus == Syn.bus(i));    
  end
  if isempty(ok)
    fm_disp(['No PV or Slack generator found at the machine bus: ', ...
            Varname.bus{Syn.bus(i)}],2)
  end
end

% search generators connected to SW and PV buses
for i = 1:SW.n
  SW_bus = find(Syn.bus(:,1) == SW.bus(i));
  if (~isempty(SW_bus))
    SW.bus(i) = 0;
  end
end
for i = 1:PV.n
  PV_bus = find(Syn.bus(:,1) == PV.bus(i));
  if (~isempty(PV_bus))
    PV.bus(i) = 0;
  end
end

if SW.n > 0
  SW.con(SW.bus == 0,:) = [];
  SW.bus = SW.bus(SW.bus>0);
  SW.n = length(SW.bus);
end
if PV.n > 0
  PV.con(PV.bus == 0,:) = [];
  PV.bus = PV.bus(PV.bus>0);
  PV.n = length(PV.bus);
end

fm_disp('Initialization of Synchronous Machines completed.')

% -----------------------------------------------------------
% function for creating warning messages
% -----------------------------------------------------------
function synwarn(idx, msg)
global Syn
fm_disp(strcat('Warning: Synchronous Machine #', ...
               int2str(idx),'(model ',num2str(Syn.con(idx,5)), ...
               ') at bus #',int2str(Syn.bus(idx)),msg))

% --------------------------------------------------------------
% Saturation function
% --------------------------------------------------------------
function output = Sat(flag,e1q,isx)
global Syn

b = [0.8*ones(length(isx),1), 1-Syn.con(isx,25), 1.2*(1-Syn.con(isx,26))];
c2 = b*[12.5; -25; 12.5];
c1 = b*[-27.5; 50; -22.5];
c0 = b*[15; -24; 10];

switch flag

 case 1 % saturation function
  output = e1q;
  idx = find(output > 0.8);
  if ~isempty(idx)
    output(idx)=c2(idx).*e1q(idx).^2+c1(idx).*e1q(idx)+c0(idx);
  end

 case 2 % Jacobian
  output = ones(length(isx),1);
  idx = find(e1q(isx) > 0.8);
  if ~isempty(idx)
    output(idx)=2*c2(idx).*e1q(idx) + c1(idx);
  end

end