fm_tcsc.m

电力系统的psat

function  fm_tcsc(flag)
% FM_TCSC defines Synchronous Machines
%
% FM_TCSC(FLAG)
%       FLAG = 0 initializations
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%       FLAG = 5 non-windup limits
%
%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 Tcsc Bus DAE Settings

type = Tcsc.con(:,6);
ty1 = find(type == 1);
ty2 = find(type == 2);
V1 = DAE.V(Tcsc.bus1);
V2 = DAE.V(Tcsc.bus2);
t1 = DAE.a(Tcsc.bus1);
t2 = DAE.a(Tcsc.bus2);
ss = sin(t1-t2);
cc = cos(t1-t2);
Pref = Tcsc.con(:,8);

if Settings.init
  Kr = Tcsc.con(:,7);
  Tw = Tcsc.con(:,13);
  T1 = Tcsc.con(:,14);
  T2 = Tcsc.con(:,15);
  T3 = Tcsc.con(:,16);
  x_max = Tcsc.con(:,9);
  x_min = Tcsc.con(:,10);
end

switch flag
 case 0
  Tcsc.B = Pref./V1./V2./ss;
  Tcsc.SI = zeros(Tcsc.n,1);
  if ty1
    Tcsc.SI(ty1) = -1./Tcsc.B(ty1);
    jdx = find(Tcsc.SI(ty1) > x_max(ty1));
    if jdx
      tcscwarn(ty1(jdx),' Xc is over its maximum limit.')
    end
    jdx = find(Tcsc.SI(ty1) < x_min(ty1));
    if jdx
      tcscwarn(ty1(jdx),' Xc is under its minimum limit.')
    end
  end
  if ty2
    xC = Tcsc.con(ty2,12);
    xL = Tcsc.con(ty2,11);
    % initialization of alpha
    ar = -pi:0.001:pi;
    for i = 1:length(ty2)
      % find a "good" initial guess in a brutal way...
      B = tcsca(ar,xC(i),xL(i));
      [val,idx] = min(abs(B-Tcsc.B(ty2(i))));
      a = ar(idx);
      err = 1;
      iter = 0;
      while abs(err) > Settings.lftol
	if iter > 20, break, end
	err = (tcsca(a,xC(i),xL(i)) - Tcsc.B(ty2(i)))/tcscda(a,xC(i),xL(i));
	a = a + err;
	iter = iter + 1;
      end
      if iter > 20
	tcscwarn(ty2(i),[' convergence not reached when initializing' ...
			 ' the firing angle alpha.'])
      end
      Tcsc.SI(ty2(i)) = a;
    end
    jdx = find(Tcsc.SI(ty2) > x_max(ty2));
    if jdx
      tcscwarn(ty2(jdx),' alpha is over its maximum limit.')
    end
    jdx = find(Tcsc.SI(ty2) < x_min(ty2));
    if jdx
      tcscwarn(ty2(jdx),' alpha is under its minimum limit.')
    end
    %Tcsc.con(ty2,7) = -Tcsc.con(ty2,7);
  end
  DAE.x(Tcsc.x1) = 0;
  DAE.x(Tcsc.x2) = 0;
  DAE.x(Tcsc.x3) = 0;

  %SI.*(1-Tcsc.con(:,15)./Tcsc.con(:,16));
  %Tcsc.con(:,8) = Pref - DAE.x(Tcsc.x1)./Tcsc.con(:,7)./Tcsc.con(:,13);

 case 1 % algebraic equations

  if Settings.init % complete model
    Tcsc.Pe = V1.*V2.*ss.*Tcsc.B;
    DAE.gp = DAE.gp + sparse(Tcsc.bus1,1, Tcsc.Pe,Bus.n,1);
    DAE.gq = DAE.gq + sparse(Tcsc.bus1,1, V1.*(V1-V2.*cc).*Tcsc.B,Bus.n,1);
    DAE.gp = DAE.gp + sparse(Tcsc.bus2,1,-Tcsc.Pe,Bus.n,1);
    DAE.gq = DAE.gq + sparse(Tcsc.bus2,1, V2.*(V2-V1.*cc).*Tcsc.B,Bus.n,1);
  else % static model used for power flow solution
    DAE.gp = DAE.gp + sparse(Tcsc.bus1,1, Pref,Bus.n,1);
    DAE.gq = DAE.gq + sparse(Tcsc.bus1,1, Pref.*(V1./V2./ss-cc./ss),Bus.n,1);
    DAE.gp = DAE.gp + sparse(Tcsc.bus2,1,-Pref,Bus.n,1);
    DAE.gq = DAE.gq + sparse(Tcsc.bus2,1, Pref.*(V2./V1./ss-cc./ss),Bus.n,1);
  end

 case 2 % algebraic Jacobians

  if Settings.init  % complete model
    a1 = ss.*Tcsc.B;
    a2 = cc.*Tcsc.B;
    a3 = V1.*V2;
    a4 = a3.*a1;
    a5 = a3.*a2;
    % dP1/dy
    DAE.J12 = DAE.J12 + sparse(Tcsc.bus1,Tcsc.bus1, V2.*a1,Bus.n,Bus.n);
    DAE.J12 = DAE.J12 + sparse(Tcsc.bus1,Tcsc.bus2, V1.*a1,Bus.n,Bus.n);
    DAE.J11 = DAE.J11 + sparse(Tcsc.bus1,Tcsc.bus1, a5,Bus.n,Bus.n);
    DAE.J11 = DAE.J11 + sparse(Tcsc.bus1,Tcsc.bus2,-a5,Bus.n,Bus.n);
    % dP2/dy
    DAE.J22 = DAE.J22 + sparse(Tcsc.bus1,Tcsc.bus1,(2*V1-V2.*cc).*Tcsc.B,Bus.n,Bus.n);
    DAE.J22 = DAE.J22 + sparse(Tcsc.bus1,Tcsc.bus2,-V1.*a2,Bus.n,Bus.n);
    DAE.J21 = DAE.J21 + sparse(Tcsc.bus1,Tcsc.bus1, a4,Bus.n,Bus.n);
    DAE.J21 = DAE.J21 + sparse(Tcsc.bus1,Tcsc.bus2,-a4,Bus.n,Bus.n);
    % dQ1/dy
    DAE.J12 = DAE.J12 + sparse(Tcsc.bus2,Tcsc.bus1,-V2.*a1,Bus.n,Bus.n);
    DAE.J12 = DAE.J12 + sparse(Tcsc.bus2,Tcsc.bus2,-V1.*a1,Bus.n,Bus.n);
    DAE.J11 = DAE.J11 + sparse(Tcsc.bus2,Tcsc.bus1,-a5,Bus.n,Bus.n);
    DAE.J11 = DAE.J11 + sparse(Tcsc.bus2,Tcsc.bus2, a5,Bus.n,Bus.n);
    % dQ2/dy
    DAE.J22 = DAE.J22 + sparse(Tcsc.bus2,Tcsc.bus2,(2*V2-V1.*cc).*Tcsc.B,Bus.n,Bus.n);
    DAE.J22 = DAE.J22 + sparse(Tcsc.bus2,Tcsc.bus1,-V2.*a2,Bus.n,Bus.n);
    DAE.J21 = DAE.J21 + sparse(Tcsc.bus2,Tcsc.bus1, a4,Bus.n,Bus.n);
    DAE.J21 = DAE.J21 + sparse(Tcsc.bus2,Tcsc.bus2,-a4,Bus.n,Bus.n);
  else %static model used for power flow solution
    ct = cot(t1-t2);
    % dQ1/dy
    DAE.J22 = DAE.J22 + sparse(Tcsc.bus1,Tcsc.bus1, Pref./V2./ss,Bus.n,Bus.n);
    DAE.J22 = DAE.J22 + sparse(Tcsc.bus1,Tcsc.bus2,-V1.*Pref./V2./V2./ss,Bus.n,Bus.n);
    DAE.J21 = DAE.J21 + sparse(Tcsc.bus1,Tcsc.bus1, Pref.*(1+ct.*ct-V1.*ct./ss./V2),Bus.n,Bus.n);
    DAE.J21 = DAE.J21 + sparse(Tcsc.bus1,Tcsc.bus2,-Pref.*(1+ct.*ct-V1.*ct./ss./V2),Bus.n,Bus.n);
    % dQ2/dy
    DAE.J22 = DAE.J22 + sparse(Tcsc.bus2,Tcsc.bus2, Pref./V1./ss,Bus.n,Bus.n);
    DAE.J22 = DAE.J22 + sparse(Tcsc.bus2,Tcsc.bus1,-V2.*Pref./V1./V1./ss,Bus.n,Bus.n);
    DAE.J21 = DAE.J21 + sparse(Tcsc.bus2,Tcsc.bus1, Pref.*(1+ct.*ct-V2.*ct./ss./V1),Bus.n,Bus.n);
    DAE.J21 = DAE.J21 + sparse(Tcsc.bus2,Tcsc.bus2,-Pref.*(1+ct.*ct-V2.*ct./ss./V1),Bus.n,Bus.n);
  end

 case 3 % differential equations

  if Settings.init
    DAE.f(Tcsc.x1) = (-Kr.*(Pref-Tcsc.Pe)-DAE.x(Tcsc.x1))./Tw;
    DAE.f(Tcsc.x2) = (DAE.x(Tcsc.x1)-DAE.x(Tcsc.x2)+Kr.*(Pref-Tcsc.Pe))./T1;
    DAE.f(Tcsc.x3) = ((1-T2./T3).*DAE.x(Tcsc.x2)-DAE.x(Tcsc.x3))./T3;
    SI = max(DAE.x(Tcsc.x3) + DAE.x(Tcsc.x2).*T2./T3+Tcsc.SI,x_min);
    SI = min(SI,x_max);
    % compute susceptances Tcsc.B
    if ~isempty(ty1), Tcsc.B(ty1) = -1./SI(ty1); end
    if ~isempty(ty2), Tcsc.B(ty2) = -tcsca(SI(ty2),Tcsc.con(ty2,12),Tcsc.con(ty2,11)); end
  end

 case 4 % state Jacobians

  if Settings.init
    DB = zeros(Tcsc.n,1);
    SI = DAE.x(Tcsc.x3) + DAE.x(Tcsc.x2).*T2./T3+Tcsc.SI;
    if ~isempty(ty1)
      DB(ty1) = -1./SI(ty1)./SI(ty1);
    end
    if ~isempty(ty2)
      DB(ty2) = tcscda(SI(ty2),Tcsc.con(ty2,12),Tcsc.con(ty2,11));
    end

    DAE.Fx = DAE.Fx + sparse(Tcsc.x1,Tcsc.x1,-1./Tw,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tcsc.x2,Tcsc.x1, 1./T1,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tcsc.x2,Tcsc.x2,-1./T1,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tcsc.x3,Tcsc.x2,(1-T2./T3)./T3,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Tcsc.x3,Tcsc.x3,-1./T3,DAE.n,DAE.n);

    a0 = V1.*V2;
    a1 = a0.*ss;
    a2 = V1-V2.*cc;
    a3 = Kr.*ss.*Tcsc.B;
    a4 = a0.*Kr.*cc.*Tcsc.B;

    DAE.Fy = DAE.Fy + sparse(Tcsc.x1,Bus.n+Tcsc.bus1,V2.*a3./Tw,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(Tcsc.x1,Bus.n+Tcsc.bus2,V1.*a3./Tw,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(Tcsc.x1,Tcsc.bus1, a4./Tw,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(Tcsc.x1,Tcsc.bus2,-a4./Tw,DAE.n,2*Bus.n);

    DAE.Fy = DAE.Fy - sparse(Tcsc.x2,Bus.n+Tcsc.bus1,V2.*a3./T1,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy - sparse(Tcsc.x2,Bus.n+Tcsc.bus2,V1.*a3./T1,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy - sparse(Tcsc.x2,Tcsc.bus1, a4./T1,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy - sparse(Tcsc.x2,Tcsc.bus2,-a4./T1,DAE.n,2*Bus.n);

    % windup limiter
    a = find(SI < x_max & SI > x_min);

    if ~isempty(a)

      DAE.Fx = DAE.Fx - sparse(Tcsc.x1(a),Tcsc.x3(a),Kr(a).*a1(a).*DB(a)./Tw,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx - sparse(Tcsc.x1(a),Tcsc.x2(a),Kr(a).*a1(a).*DB(a).*T2(a)./T3(a)./Tw,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Tcsc.x2(a),Tcsc.x3(a),Kr(a).*a1(a).*DB(a)./T1,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Tcsc.x2(a),Tcsc.x2(a),Kr(a).*a1(a).*DB(a).*T2(a)./T3(a)./T1,DAE.n,DAE.n);

      DAE.Gx = DAE.Gx + sparse(Tcsc.bus1(a),Tcsc.x3(a),a1(a).*DB(a),2*Bus.n,DAE.n);
      DAE.Gx = DAE.Gx + sparse(Bus.n+Tcsc.bus1(a),Tcsc.x3(a),V1(a).*a2(a).*DB(a),2*Bus.n,DAE.n);
      DAE.Gx = DAE.Gx + sparse(Tcsc.bus2(a),Tcsc.x3(a),-a1(a).*DB(a),2*Bus.n,DAE.n);
      DAE.Gx = DAE.Gx + sparse(Bus.n+Tcsc.bus2(a),Tcsc.x3(a),V2(a).*a2(a).*DB(a),2*Bus.n,DAE.n);

      DAE.Gx = DAE.Gx + sparse(Tcsc.bus1(a),Tcsc.x2(a),a1(a).*DB(a).*T2(a)./T3(a),2*Bus.n,DAE.n);
      DAE.Gx = DAE.Gx + sparse(Bus.n+Tcsc.bus1(a),Tcsc.x2(a),V1(a).*a2(a).*DB(a).*T2(a)./T3(a),2*Bus.n,DAE.n);
      DAE.Gx = DAE.Gx + sparse(Tcsc.bus2(a),Tcsc.x2(a),-a1(a).*DB(a).*T2(a)./T3(a),2*Bus.n,DAE.n);
      DAE.Gx = DAE.Gx + sparse(Bus.n+Tcsc.bus2(a),Tcsc.x2(a),V2(a).*a2(a).*DB(a).*T2(a)./T3(a),2*Bus.n,DAE.n);

    end
  end

end

% --------------------------------------------------------------------------------------
function Be = tcsca(af,xC,xL) % Be(alpha)

kx1 = sqrt(xC./xL);
kx2 = kx1.*kx1;
kx3 = kx2.*kx1;
kx4 = kx3.*kx1;
ckf = cos(kx1.*(pi-af));
skf = sin(kx1.*(pi-af));
s2a = sin(2*af);
caf = cos(af);
saf = sin(af);
Be = (pi*(kx4-2*kx2+1).*ckf)./(xC.*((pi*kx4-pi-2*kx4.*af+2*af.*kx2- ...
      kx4.*s2a+kx2.*s2a-4*kx2.*caf.*saf).*ckf -4*kx3.*caf.*caf.*skf));

% --------------------------------------------------------------------------------------
function DB = tcscda(af,xC,xL) % -dBe/d(alpha)

kx1 = sqrt(xC./xL);
kx2 = kx1.*kx1;
kx3 = kx2.*kx1;
kx4 = kx3.*kx1;
ckf = cos(kx1.*(-pi+af));
skf = sin(kx1.*(-pi+af));
ck2 = ckf.*ckf;
c2a = cos(2*af);
s2a = sin(2*af);
caf = cos(af);
saf = sin(af);
ca2 = caf.*caf;
DB = -2*pi*(-kx4+2*kx2-1).*(2*skf.*skf.*kx1.*kx3.*ca2-ck2.*kx4+ck2.*kx2- ...
      ck2.*kx4.*c2a+ck2.*kx2.*c2a+2*ck2.*kx2.*saf.*saf-2*ck2.*kx2.*ca2- ...
      4*ckf.*kx3.*caf.*skf.*saf+2*kx3.*ca2.*ck2.*kx1)./xC./(ckf.*pi.*kx4- ...
      ckf.*pi-2*ckf.*kx4.*af+2*ckf.*af.*kx2-ckf.*kx4.*s2a+ckf.*kx2.*s2a- ...
      4*ckf.*kx2.*caf.*saf+4*kx3.*ca2.*skf).^2;

% --------------------------------------------------------------------------------------
% function for creating warning messages
function tcscwarn(idx, msg)
global Tcsc
fm_disp(strcat('Warning: TCSC #',int2str(idx),' between buses #', ...
	       int2str(Tcsc.bus1(idx)),' and #',int2str(Tcsc.bus2(idx)),msg))