fm_ltc.m

这是一个很适合研究和学习用的电力系统仿真软件

function fm_ltc(flag)
% FM_LTC defines an Load Tap Changer with a continuous or discrete model
%        and three different control models:
%        (a) secondary winding voltage
%        (b) reactive power
%        (c) remote voltage
%
% FM_LTC(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
%Update:    07-Oct-2003
%Version:   1.0.1
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global Bus Ltc DAE jay Settings

persistent iters
persistent delay

%Ltc.dat:
%      1.  longitudinal admittance
%      2.  y11
%      3.  y12 == y21
%      4.  y22
%      5.  complex power at bus 1
%      6.  complex power at bus 2
%      7.  regulated bus

mltc = DAE.x(Ltc.m);

bus1 = Ltc.bus1;
bus2 = Ltc.bus2;
busr = Ltc.dat(:,7);
v1 = DAE.V(bus1);
v2 = DAE.V(bus2);
ang1 = DAE.a(bus1);
ang2 = DAE.a(bus2);
hltc = Ltc.con(:,7);
kltc = Ltc.con(:,8);
mmax = Ltc.con(:,9);
mmin = Ltc.con(:,10);
if ~isempty(mltc), mltc = max(mmin,mltc); end
mstep = Ltc.con(:,11);
vrif = Ltc.con(:,12);
im1 = 1./mltc;
im2 = im1./mltc;
im3 = im2./mltc;

switch flag
 case 0

  iters = 0;
  delay = 0;
  xl = Ltc.con(:,13);
  rl = Ltc.con(:,14);
  zl = rl + jay*xl;
  yl= 1./zl;
  Ltc.dat(:,1) = yl;
  Ltc.con(find(Ltc.con(:,16) == 3),16) = 1;
  Ltc.dat(:,7) = Ltc.bus2;
  idx = find(Ltc.con(:,15));
  if ~isempty(idx)
    Ltc.dat(idx,7) = Bus.int(round(Ltc.con(idx,15)));
  end
  idx = find(Ltc.con(:,12) <= 0);
  if ~isempty(idx), Ltc.con(:,12) = 1; end
  %DAE.x(Ltc.m) = ones(Ltc.n,1);

 case 1

  for i = 1:Ltc.n

    % Voltage in rectangular coordinates
    V_rect = [v1(i);v2(i)].*exp(jay*[ang1(i);ang2(i)]);

    % LTC admittance matrix
    a = Ltc.dat(i,1)*im1(i);
    b = Ltc.dat(i,1)*(im2(i)-im1(i));
    c = Ltc.dat(i,1)*(1-im1(i));
    Ltc.dat(i,2) = (a+b);
    Ltc.dat(i,3) = -a;
    Ltc.dat(i,4) = (a+c);

    idx = [bus1(i),bus2(i)];

    % complex power vector
    S = V_rect.*conj([Ltc.dat(i,2), Ltc.dat(i,3); ...
		      Ltc.dat(i,3), Ltc.dat(i,4)]*V_rect);
    DAE.gp(idx) = DAE.gp(idx) + real(S);
    DAE.gq(idx) = DAE.gq(idx) + imag(S);
    Ltc.dat(i,5) = S(1);
    Ltc.dat(i,6) = S(2);

  end

 case 2

  n = [1 2];
  for i = 1:Ltc.n

    Y = [Ltc.dat(i,2), Ltc.dat(i,3); Ltc.dat(i,3), Ltc.dat(i,4)];
    E = exp(jay*[ang1(i);ang2(i)]);
    V = [v1(i);v2(i)].*E;
    I = Y*V;

    diagVc = [V(1), 0; 0, V(2)];
    diagVn = [E(1), 0; 0, E(2)];
    diagIc = [I(1), 0; 0, I(2)];

    dS = diagVc*conj(Y*diagVn)+conj(diagIc)*diagVn;
    J12 = real(dS);
    J22 = imag(dS)+1e-6*[1 0; 0 1];

    dS = j*diagVc*conj(diagIc-Y*diagVc);
    J11 = real(dS);
    J21 = imag(dS)+1e-6*[1 0; 0 1];

    idx = [bus1(i),bus2(i)];

    DAE.J11(idx,idx) = DAE.J11(idx,idx) + J11;
    DAE.J12(idx,idx) = DAE.J12(idx,idx) + J12;
    DAE.J21(idx,idx) = DAE.J21(idx,idx) + J21;
    DAE.J22(idx,idx) = DAE.J22(idx,idx) + J22;

  end

 case 3

  error = zeros(Ltc.n,1);
  ty1 = find(Ltc.con(:,16) == 1);
  ty2 = find(Ltc.con(:,16) == 2);
  if ty1
    error(ty1) = DAE.V(busr(ty1))-vrif(ty1);
  end
  if ty2
    error(ty2) = -vrif(ty2)-imag(Ltc.dat(ty2,6));
  end
  errel = error./vrif;

  DAE.f(Ltc.m) = (Ltc.con(:,11) <= 0).*(kltc.*error-hltc.*mltc);
  relay = (errel > mstep) - (errel < -mstep);

  % delay tap change for NR stability
  if ~isempty(find(relay))
    iters = iters + 1;
  end

  % delay tap change during TDS
  if DAE.t > 0
    if DAE.t-delay > 0.05 % fixed delay of 0.05 s
      delay = DAE.t;
    else
      relay = 0*relay;
    end
  end

  if iters == 2
    DAE.x(Ltc.m) = DAE.x(Ltc.m)+relay.*Ltc.con(:,11);
    iters = 0;
  end

  idx = find(DAE.x(Ltc.m) >= mmax & DAE.f(Ltc.m) > 0);
  if idx
    DAE.f(Ltc.m(idx)) = 0;
  end
  idx = find(DAE.x(Ltc.m) <= mmin & DAE.f(Ltc.m) < 0);
  if idx
    DAE.f(Ltc.m(idx)) = 0;
  end
  DAE.x(Ltc.m) = min(DAE.x(Ltc.m),mmax);
  DAE.x(Ltc.m) = max(DAE.x(Ltc.m),mmin);

 case 4

  gl = real(Ltc.dat(:,1));
  bl = imag(Ltc.dat(:,1));
  c12 = cos(ang1-ang2);
  s12 = sin(ang1-ang2);
  c21 = c12;
  s21 = -s12;

  for i=1:Ltc.n
    k=Ltc.m(i);
    if ((mltc(i) >= mmax(i) | mltc(i) <= mmin(i)) & ...
        DAE.f(Ltc.m(i)) == 0) | Ltc.con(i,11)
      DAE.Fx(k,k) = -1;
      DAE.Gx(:,k) = zeros(2*Bus.n,1);
      DAE.Fy(k,:) = zeros(1,2*Bus.n);
    else
      DAE.Gx(bus1(i),k) = ...
          v1(i)*(-2*im3(i)*gl(i)*v1(i) + ...
                 v2(i)*im2(i)*(gl(i)*c12(i) + bl(i)*s12(i)));
      DAE.Gx(bus1(i)+Bus.n,k) = ...
          v1(i)*(2*im3(i)*bl(i)*v1(i) - ...
                 v2(i)*im2(i)*(-gl(i)*s12(i) + bl(i)*c12(i)));
      DAE.Gx(bus2(i),k) = v1(i)*v2(i)*im2(i)*(gl(i)*c21(i) + bl(i)*s21(i));
      DAE.Gx(bus2(i)+Bus.n,k) = v1(i)*v2(i)*im2(i)*(gl(i)*s21(i) - bl(i)*c21(i));
      switch Ltc.con(i,16)
       case 1
	DAE.Fx(k,k) = -hltc(i);
	if DAE.Fx(k,k) == 0
	  DAE.Fx(k,k) = -0.0001;
	end
	DAE.Fy(k,busr(i)+Bus.n)= kltc(i);
       case 2
	Y = [Ltc.dat(i,2), Ltc.dat(i,3); Ltc.dat(i,3), Ltc.dat(i,4)];
	G = real(Y);
	B = imag(Y);
	gp = -real(Ltc.dat(i,6));
	gq =  imag(Ltc.dat(i,6));
	DAE.Fx(k,k) = -hltc(i) - kltc(i)*(DAE.Gx(bus2(i)+Bus.n,k));
	DAE.Fy(k,bus1(i)) = kltc(i)*v1(i)*v2(i)*(G(2,1)*c12(i)-B(2,1)*s12(i));
	DAE.Fy(k,bus1(i)+Bus.n) = kltc(i)*v2(i)*(G(2,1)*s12(i)+B(2,1)*c12(i));
	DAE.Fy(k,bus2(i)) = kltc(i)*(gp+G(2,2)*v2(i)*v2(i));

        V2 = v2(i);
        if V2 <= 0
          V2 = 1;
        end

        DAE.Fy(k,bus2(i)+Bus.n) = -kltc(i)*(gq/V2-B(2,2)*v2(i));
      end
    end
  end

 case 5

  % hard limits (non-windup)
  idx = find(((mltc >= mmax | mltc <= mmin) & DAE.f(Ltc.m)== 0) | Ltc.con(:,11));
  if ~isempty(idx)
    k = Ltc.m(idx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    DAE.Ac(k,k) = speye(length(idx));
  end

end