fm_phs.m

电力系统的psat

function  fm_phs(flag)
%FM_PHS defines Phase Shifting Transformer
%
%Data Format Phs.con:
%          col #1: Bus 1 number
%          col #2: Bus 2 number
%          col #3: Power rate [MVA]
%          col #4: Bus 1 Voltage Rate [kV]
%          col #5: Bus 2 Voltage Rate [kV]
%          col #6: Frequency rate [Hz]
%          col #7: Tm Measurement Time Constant [Second]
%          col #8: Kp None [None]
%          col #9: Ki None [None]
%          col #10: Pref None [None]
%          col #11: rt None [None]
%          col #12: xt None [None]
%          col #13: alpha_max
%          col #14: alpha_max
%
%Author:    Federico Milano
%Date:      16-May-2003
%Update:    01-Aug-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 Phs DAE Bus Settings

a = DAE.x(Phs.alpha);
Pm = DAE.x(Phs.Pm);
V1 = DAE.V(Phs.bus1);
V2 = DAE.V(Phs.bus2);
t1 = DAE.a(Phs.bus1);
t2 = DAE.a(Phs.bus2);
Tm = Phs.con(:,7);
Kp = -Phs.con(:,8);
Ki = -Phs.con(:,9);
Pref = Phs.con(:,10);
rt = Phs.con(:,11);
xt = Phs.con(:,12);
zt = rt+i*xt;
yt = 1./zt;
gt = real(yt);
bt = imag(yt);
a_max = Phs.con(:,13);
a_min = Phs.con(:,14);

s12a = sin(t1-t2-a);
c12a = cos(t1-t2-a);
V12 = V1.*V2;

switch flag
 case 1 % algebraic equations

  a1 = bt.*s12a+gt.*c12a;
  a2 =-bt.*s12a+gt.*c12a;
  a3 = bt.*c12a+gt.*s12a;
  a4 =-bt.*c12a+gt.*s12a;
  V11 = V1.*V1;
  V22 = V2.*V2;
  DAE.gp = DAE.gp + sparse(Phs.bus1,1, V11.*gt-V12.*a1,Bus.n,1);
  DAE.gq = DAE.gq + sparse(Phs.bus1,1,-V11.*bt-V12.*a4,Bus.n,1);
  DAE.gp = DAE.gp + sparse(Phs.bus2,1, V22.*gt-V12.*a2,Bus.n,1);
  DAE.gq = DAE.gq + sparse(Phs.bus2,1,-V22.*bt+V12.*a3,Bus.n,1);

case 2 % algebraic Jacobians

  b1 =  bt.*s12a+gt.*c12a;
  b2 =  bt.*c12a-gt.*s12a;
  b3 = -bt.*s12a+gt.*c12a;
  b4 = -bt.*c12a-gt.*s12a;

  DAE.J12 = DAE.J12 + sparse(Phs.bus1,Phs.bus1,2.*V1.*gt-V2.*b1,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(Phs.bus1,Phs.bus2,-V1.*b1,Bus.n,Bus.n);
  DAE.J11 = DAE.J11 + sparse(Phs.bus1,Phs.bus1,-V12.*b2,Bus.n,Bus.n);
  DAE.J11 = DAE.J11 + sparse(Phs.bus1,Phs.bus2, V12.*b2,Bus.n,Bus.n);

  DAE.J22 = DAE.J22 + sparse(Phs.bus1,Phs.bus1,-2.*V1.*bt+V2.*b2,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(Phs.bus1,Phs.bus2, V1.*b2,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(Phs.bus1,Phs.bus1,-V12.*b1,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(Phs.bus1,Phs.bus2, V12.*b1,Bus.n,Bus.n);

  DAE.J12 = DAE.J12 + sparse(Phs.bus2,Phs.bus1,-V2.*b3,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(Phs.bus2,Phs.bus2,2.*V2.*gt-V1.*b3,Bus.n,Bus.n);
  DAE.J11 = DAE.J11 + sparse(Phs.bus2,Phs.bus1,-V12.*b4,Bus.n,Bus.n);
  DAE.J11 = DAE.J11 + sparse(Phs.bus2,Phs.bus2, V12.*b4,Bus.n,Bus.n);

  DAE.J22 = DAE.J22 + sparse(Phs.bus2,Phs.bus1,-V2.*b4,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(Phs.bus2,Phs.bus2,-2.*V2.*bt-V1.*b4,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(Phs.bus2,Phs.bus1, V12.*b3,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(Phs.bus2,Phs.bus2,-V12.*b3,Bus.n,Bus.n);

 case 3 % differential equations

  errP = V1.*V1.*gt-V12.*(bt.*s12a+gt.*c12a)-Pm;

  DAE.f(Phs.alpha) = -Kp.*errP./Tm+Ki.*(Pref-Pm);
  DAE.f(Phs.Pm) = errP./Tm;

  % non-windup limits
  idx = find(a >= a_max & (DAE.f(Phs.alpha) > 0 | ~Settings.init));
  if idx, DAE.f(Phs.alpha(idx)) = 0; end
  DAE.x(Phs.alpha) = min(DAE.x(Phs.alpha),a_max);
  idx = find(a <= a_min & (DAE.f(Phs.alpha) < 0 | ~Settings.init));
  if idx, DAE.f(Phs.alpha(idx)) = 0; end
  DAE.x(Phs.alpha) = max(DAE.x(Phs.alpha),a_min);

 case 4 % state variable Jacobians

  V11 = V1.*V1;
  V22 = V2.*V2;
  b1 =  bt.*s12a+gt.*c12a;
  b2 =  V1.*b1;
  b3 =  2.*V1.*gt-V2.*b1;
  b4 =  V12.*(bt.*c12a-gt.*s12a);
  b5 =  V12.*(bt.*c12a+gt.*s12a);
  b6 =  V12.*(gt.*c12a-bt.*s12a);

  idx = find(a <= a_max & a >= a_min & DAE.f(Phs.alpha) ~= 0);
  if ~isempty(idx)

    ai  = Phs.alpha(idx);
    b1i = Phs.bus1(idx);
    b2i = Phs.bus2(idx);
    pmi = Phs.Pm(idx);

    % DAE.Fx
    DAE.Fx = DAE.Fx + sparse(ai,pmi,Kp(idx)./Tm(idx)-Ki(idx),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(pmi,ai,b4(idx)./Tm(idx),DAE.n,DAE.n);

    % DAE.Fy
    DAE.Fy = DAE.Fy + sparse(ai,b1i+Bus.n,-Kp(idx).*b3(idx)./Tm(idx),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(ai,b2i+Bus.n,Kp(idx).*b2(idx)./Tm(idx),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(ai,b1i,Kp(idx).*b4(idx)./Tm(idx),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(ai,b2i,-Kp(idx).*b4(idx)./Tm(idx),DAE.n,2*Bus.n);

    % DAE.Gx
    DAE.Gx = DAE.Gx + sparse(b1i,ai,b4(idx),2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(b1i+Bus.n,ai,-V12(idx).*b1(idx),2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(b2i,ai,-b5(idx),2*Bus.n,DAE.n);
    DAE.Gx = DAE.Gx + sparse(b2i+Bus.n,ai,-b6(idx),2*Bus.n,DAE.n);

  end

  % DAE.Fx
  DAE.Fx = DAE.Fx + sparse(Phs.alpha,Phs.alpha,-Kp.*b4./Tm,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Phs.Pm,Phs.Pm,-1./Tm,DAE.n,DAE.n);

  % DAE.Fy
  DAE.Fy = DAE.Fy + sparse(Phs.Pm,Phs.bus1+Bus.n,b3./Tm,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Phs.Pm,Phs.bus2+Bus.n,-b2./Tm,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Phs.Pm,Phs.bus1,-b4./Tm,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Phs.Pm,Phs.bus2,b4./Tm,DAE.n,2*Bus.n);

case 5 % non-windup limiters

  idx = find((a >= a_max | a <= a_min) & DAE.f(Phs.alpha) == 0);
  if ~isempty(idx)
    k = Phs.alpha(idx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    if Settings.octave
      DAE.Ac(k,k) = -eye(length(idx));
    else
      DAE.Ac(k,k) = -speye(length(idx));
    end
  end

end

% -------------------------------------------------------------------
% function for creating warning messages
function phshifterwarn(idx, msg)
fm_disp(strcat('Warning: Phase shifter #',int2str(idx),msg))