fm_pss.m

基于PSAT 软件的多目标最优潮流计算用于中小型电力系统的分析和管理

function  fm_pss(flag)
% FM_PSS define Power System Stabilizers
%
% FM_PSS(FLAG)
%       FLAG = 0 initialization
%       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-2006 Federico Milano

global Exc Syn DAE Bus Pss Settings

type = Pss.con(:,2);
ty1 = find(type == 1);
ty2 = find(type == 2 | type == 4);
ty3 = find(type == 3 | type == 5);
tya = find(type > 1);
tyb = find(type > 3);

VSI = zeros(Pss.n,1);
SIw = find(Pss.con(:,3) == 1);
SIp = find(Pss.con(:,3) == 2);
SIv = find(Pss.con(:,3) == 3);
if SIw, VSI = DAE.x(Syn.omega(Pss.syn(SIw))); end
if SIp, VSI = -Syn.Pg(Pss.syn(SIp)); end
if SIv, VSI = DAE.V(Pss.bus(SIv)); end

Tw = Pss.con(:,7);
T1 = Pss.con(:,8);
T2 = Pss.con(:,9);
T3 = Pss.con(:,10);
T4 = Pss.con(:,11);
Ta = Pss.con(:,13);

Kw = Pss.con(:,6);
Ka = Pss.con(:,12);
Kp = Pss.con(ty1,14);
Kv = Pss.con(ty1,15);

vsmax = Pss.con(:,4);
vsmin = Pss.con(:,5);
vamax = Pss.con(:,16);
vathr = Pss.con(:,17);
v3max = Pss.con(:,18);
v3min = Pss.con(:,19);

ETHR = Pss.con(:,20);
WTHR = Pss.con(:,21);

S2 = Pss.con(:,22);

switch flag
 case 0

  idx = find(Tw == 0);
  if idx,
    Pss.con(idx,6) = 0.01;
    psswarn(idx,[' Tw cannot be zero. Default value Tw = 0.01 will' ...
		 ' be used.'])
  end
  if ty1
    DAE.x(Pss.v1(ty1)) = -Kw(ty1)+Kp.*Syn.Pg(Pss.syn(ty1))-Kv.*DAE.V(Pss.bus(ty1));
  end
  if tya
    idx = find(T2(tya) == 0);
    if idx,
      Pss.con(tya(idx),9) = 0.01;
      psswarn(idx,[' T2 cannot be zero. Default value T2 = 0.01 will' ...
		   ' be used.'])
    end
    idx = find(T4(tya) == 0);
    if idx,
      Pss.con(tya(idx),11) = 0.01;
      psswarn(idx,[' T4 cannot be zero. Default value T4 = 0.01 will' ...
		   ' be used.'])
    end
    DAE.x(Pss.v1(tya)) = -Kw(tya).*VSI(tya);
    DAE.x(Pss.v2(tya)) = 0;
    DAE.x(Pss.v3(tya)) = 0;
  end
  if tyb
    idx = find(Ta(tyb) == 0);
    if idx,
      Pss.con(tyb(idx),13) = 0.01;
      psswarn(idx,[' Ta cannot be zero. Default value Ta = 0.01 will' ...
		   ' be used.'])
    end
    DAE.x(Pss.va(tyb)) = 0;
  end
  Pss.con(:,22) = ~Pss.con(:,22);
  DAE.x(Pss.vss) = 0;

  fm_disp('Initialization of Power System Stabilizers completed.')

 case 3

  if ty1
    VS = -Kw(ty1).*DAE.x(Syn.omega(Pss.syn(ty1))) + ...
	Kp.*Syn.Pg(Pss.syn(ty1))-Kv.*DAE.V(Pss.bus(ty1))- ...
	DAE.x(Pss.v1(ty1));
    %VS = round(VS/Settings.dyntol)*Settings.dyntol;
    DAE.f(Pss.v1(ty1)) = VS./Tw(ty1);
    Pss.Vs(ty1) = -VS;
  end
  if tya
    y = Kw.*VSI+DAE.x(Pss.v1);
    DAE.f(Pss.v1(tya)) = -y(tya)./Tw(tya);
  end
  if ty2
    A = T1(ty2)./T2(ty2);
    B = 1-A;
    C = T3(ty2)./T4(ty2);
    D = 1-C;
    DAE.f(Pss.v2(ty2)) = (B.*y(ty2)-DAE.x(Pss.v2(ty2)))./T2(ty2);
    DAE.f(Pss.v3(ty2)) = (D.*(DAE.x(Pss.v2(ty2))+A.*y(ty2))- ...
			  DAE.x(Pss.v3(ty2)))./T4(ty2);
    VS = DAE.x(Pss.v3(ty2))+C.*(DAE.x(Pss.v2(ty2))+A.*y(ty2));
    Pss.Vs(ty2) = VS; % round(VS/Settings.dyntol)*Settings.dyntol;
  end
  if ty3
    A = T3(ty3)./T4(ty3);
    B = 1-A;
    C = T1(ty3)-T2(ty3).*A;
    a12 = -1./T4(ty3);
    a13 = a12.*C;
    a22 = T2(ty3)./T4(ty3);
    a23 = a22.*C - B;
    DAE.f(Pss.v2(ty3)) = -a12.*DAE.x(Pss.v3(ty3))-a13.*y(ty3);
    DAE.f(Pss.v3(ty3)) = -a22.*DAE.x(Pss.v3(ty3))- ...
	DAE.x(Pss.v2(ty3))-a23.*y(ty3);
    VS = DAE.x(Pss.v2(ty3))+A.*y(ty3);
    Pss.Vs(ty3) = VS; %round(VS/Settings.dyntol)*Settings.dyntol;
  end
  if tyb
    Pss.s1 = (Pss.s1 | (Syn.vf(Pss.syn) < ETHR)) & ...
             (DAE.x(Syn.omega(Pss.syn)) >= WTHR);
    DAE.f(Pss.va(tyb)) = Pss.s1(tyb).*(Ka(tyb).*VSI(tyb)- ...
                                       DAE.x(Pss.va(tyb)))./Ta(tyb);
    idx = find(DAE.x(Pss.va(tyb)) > vamax(tyb) & DAE.f(Pss.va(tyb)) > 0);
    DAE.f(Pss.va(tyb(idx))) = 0;
    DAE.x(Pss.va(tyb)) = min(DAE.x(Pss.va(tyb)),vamax(tyb));
    s2 = ((DAE.x(Syn.omega(Pss.syn))-1) > 0) | S2;
    va = min(s2(tyb).*DAE.x(Pss.va(tyb)),vathr(tyb));
    va = max(va,0);
    Pss.Vs(tyb) = min(Pss.Vs(tyb),v3max(tyb));
    Pss.Vs(tyb) = max(Pss.Vs(tyb),v3min(tyb));
    Pss.Vs(tyb) = Pss.Vs(tyb)+va;
  end

  vss = DAE.x(Pss.vss);
  DAE.f(Pss.vss) = 1000*(Pss.Vs-vss);
  % non-windup limiter
  idx = find(vss >= vsmax & DAE.f(Pss.vss) > 0);
  if idx, DAE.f(Pss.vss) = 0; end
  DAE.x(Pss.vss) = min(DAE.x(Pss.vss),vsmax);
  idx = find(vss <= vsmin & DAE.f(Pss.vss) < 0);
  if idx, DAE.f(Pss.vss) = 0; end
  DAE.x(Pss.vss) = max(DAE.x(Pss.vss),vsmin);

 case 4

  Kw = zeros(Pss.n,1);
  Kp = zeros(Pss.n,1);
  Kv = zeros(Pss.n,1);
  Kw(SIw) = Pss.con(SIw,6);
  Kp(SIp) = Pss.con(SIp,6);
  Kv(SIv) = Pss.con(SIv,6);
  Kw(ty1) = Pss.con(ty1,6);
  Kp(ty1) = Pss.con(ty1,14);
  Kv(ty1) = Pss.con(ty1,15);
  
  % common Jacobians elements
  DAE.Fx = DAE.Fx + sparse(Pss.vss,Pss.vss,-1000,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Pss.v1,Pss.v1,-1./Tw,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Pss.v1,Syn.omega(Pss.syn),-Kw./Tw,DAE.n,DAE.n);
  DAE.Fy = DAE.Fy + sparse(Pss.v1,Pss.bus+Bus.n,-Kv./Tw,DAE.n,2*Bus.n);
  for i = 1:Pss.n, pssjac(Pss.syn(i),Pss.v1(i),Kp(i)/Tw(i)); end

  vss = DAE.x(Pss.vss);

  if ty1
    idx = find(vss(ty1) < vsmax(ty1) & vss(ty1) > vsmin(ty1));
    if ~isempty(idx)
      k = ty1(idx);
      h = Pss.syn(k);
      m = Pss.exc(k);
      DAE.Fx = DAE.Fx + sparse(Pss.vss(k),Syn.omega(h),1000*Kw(k),DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Pss.vss(k),Pss.v1(k),1000,DAE.n,DAE.n);
      DAE.Fy = DAE.Fy + sparse(Pss.vss(k),Pss.bus(k)+Bus.n,1000*Kv(k),DAE.n,2*Bus.n);
      vssexc(Pss.vss(k),Exc.vm(m));
      for i = 1:length(k), pssjac(k(i),h(i),-Kp(k(i))); end
    end
  end

  if ty2
    A = T1(ty2)./T2(ty2);
    B = 1-A;
    C = T3(ty2)./T4(ty2);
    D = 1-C;
    E = C.*A;
    F = D./T4(ty2);
    G = B./T2(ty2);

    DAE.Fx = DAE.Fx + sparse(Pss.v2(ty2),Pss.v1(ty2), G,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Pss.v2(ty2),Pss.v2(ty2),-1./T2(ty2),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Pss.v2(ty2),Syn.omega(Pss.syn(ty2)), ...
                             G.*Kw(ty2),DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(Pss.v2(ty2),Pss.bus(ty2)+Bus.n, ...
                             G.*Kv(ty2),DAE.n,2*Bus.n);
    for i = 1:length(ty2)
      pssjac(Pss.syn(ty2(i)),Pss.v2(ty2(i)),-G(i).*Kp(ty2(i)))
    end

    DAE.Fx = DAE.Fx + sparse(Pss.v3(ty2),Pss.v1(ty2), F.*A,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Pss.v3(ty2),Pss.v2(ty2), F,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Pss.v3(ty2),Pss.v3(ty2),-1./T4(ty2),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Pss.v3(ty2),Syn.omega(Pss.syn(ty2)), ...
                             F.*A.*Kw(ty2),DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(Pss.v3(ty2),Pss.bus(ty2)+Bus.n, ...
                             F.*A.*Kv(ty2),DAE.n,2*Bus.n);
    for i = 1:length(ty2)
      pssjac(Pss.syn(ty2(i)),Pss.v3(ty2(i)),-F(i).*A(i).*Kp(ty2(i)))
    end

    idx = find(vss(ty2) < vsmax(ty2) & vss(ty2) > vsmin(ty2));
    if ~isempty(idx)
      k = ty2(idx);
      h = Pss.syn(k);
      m = Pss.exc(k);
      DAE.Fx = DAE.Fx + sparse(Pss.vss(k),Pss.v3(k),1000,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Pss.vss(k),Pss.v2(k),1000*C(k),DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Pss.vss(k),Pss.v1(k),1000*E(k),DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Pss.vss(k),Syn.omega(h),1000*Kw(k).*E(k),DAE.n,DAE.n);
      DAE.Fy = DAE.Fy + sparse(Pss.vss(k),Pss.bus(k)+Bus.n,1000*Kv(k).*E(k),DAE.n,2*Bus.n);
      vssexc(Pss.vss(k),Exc.vm(m));
      for i = 1:length(k), pssjac(k(i),h(i),-E(k(i))*Kp(k(i))); end
    end
  end

  if ty3
    A = T3(ty3)./T4(ty3);
    B = 1-A;
    C = T1(ty3)-T2(ty3).*A;
    a12 = -1./T4(ty3);
    a13 = a12.*C;
    a22 = T2(ty3)./T4(ty3);
    a23 = a22.*C - B;

    DAE.Fx = DAE.Fx + sparse(Pss.v2(ty3),Pss.v1(ty3),-a13,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Pss.v2(ty3),Pss.v3(ty3),-a12,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Pss.v2(ty3),Syn.omega(Pss.syn(ty3)), ...
                             -a13.*Kw(ty3),DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(Pss.v2(ty3),Pss.bus(ty3)+Bus.n, ...
                             -a13.*Kv(ty3),DAE.n,2*Bus.n);
    for i = 1:length(ty3)
      pssjac(Pss.syn(ty3(i)),Pss.v2(ty3(i)),a13(i).*Kp(ty3(i)))
    end

    DAE.Fx = DAE.Fx + sparse(Pss.v3(ty3),Pss.v1(ty3),-a23,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Pss.v3(ty3),Pss.v2(ty3),-1,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Pss.v3(ty3),Pss.v3(ty3),-a22,DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Pss.v3(ty3),Syn.omega(Pss.syn(ty3)), ...
                             -a23.*Kw(ty3),DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(Pss.v3(ty3),Pss.bus(ty3)+Bus.n, ...
                             -a23.*Kv(ty3),DAE.n,2*Bus.n);
    for i = 1:length(ty3)
      pssjac(Pss.syn(ty3(i)),Pss.v3(ty3(i)),a23(i).*Kp(ty3(i)))
    end

    idx = find(vss(ty3) < vsmax(ty3) & vss(ty3) > vsmin(ty3));
    if ~isempty(idx)
      k = ty3(idx);
      h = Pss.syn(k);
      m = Pss.exc(k);
      DAE.Fx = DAE.Fx + sparse(Pss.vss(k),Pss.v2(k),1000,DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Pss.vss(k),Pss.v1(k),1000*A(k),DAE.n,DAE.n);
      DAE.Fx = DAE.Fx + sparse(Pss.vss(k),Syn.omega(h),1000*Kw(k).*A(k),DAE.n,DAE.n);
      DAE.Fy = DAE.Fy + sparse(Pss.vss(k),Pss.bus(k)+Bus.n,1000*Kv(k).*A(k),DAE.n,2*Bus.n);
      vssexc(Pss.vss(k),Exc.vm(m));
      for i = 1:length(k), pssjac(k(i),h(i),-A(k(i))*Kp(k(i))); end
    end
  end

  if tyb
    Kw(SIw) = Pss.con(SIw,12);
    Kp(SIp) = Pss.con(SIp,12);
    Kv(SIv) = Pss.con(SIv,12);

    DAE.Fx = DAE.Fx + sparse(Pss.va(tyb),Pss.va(tyb), ...
                             -Pss.s1(tyb)./Ta(tyb),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Pss.va(tyb),Syn.omega(Pss.syn(tyb)), ...
                             Pss.s1(tyb).*Kw(tyb)./Ta(tyb),DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(Pss.va(tyb),Pss.bus(tyb)+Bus.n, ...
                             Pss.s1(tyb).*Kv(tyb)./Ta(tyb),DAE.n,2*Bus.n);
    for i = 1:length(tyb)
      pssjac(Pss.syn(tyb(i)),Pss.va(tyb(i)),-Pss.s1(tyb(i))* ...
	      Kp(tyb(i))/Ta(tyb(i)))
    end

    s2 = ((DAE.x(Syn.omega(Pss.syn))-1) > 0) | S2;
    va = min(s2(tyb).*DAE.x(Pss.va(tyb)),vathr(tyb));

    idx = find(va(i) < vamax(tyb(i)) & va(i) < vathr(tyb(i)) & ...
            va(i) > 0 & Pss.vss(tyb(i)) > vsmin(tyb(i)) & ...
            Pss.vss(tyb(i)) < vsmax(tyb(i)));
    if ~isempty(idx)
      DAE.Fx = DAE.Fx + sparse(Pss.vss(tyb(idx)),Pss.va(tyb(idx)),1000,DAE.n,DAE.n);
    end
  end

 case 5

  % non-windup limter for the additional signal Va
  if tyb
    idx = find(DAE.x(Pss.va(tyb)) > vamax(tyb) & DAE.f(Pss.va(tyb)) == 0);
    if ~isempty(idx)
      k = Pss.va(tyb(idx));
      DAE.tn(k) = 0;
      DAE.Ac(:,k) = 0;
      DAE.Ac(k,:) = 0;
      DAE.Ac(k,k) = -speye(length(idx));
    end
  end

  vss = DAE.x(Pss.vss);
  idx = find((vss >= vsmax | vss <= vsmin) & ...
             DAE.f(Pss.vss) == 0);
  if ~isempty(idx)
    k = Pss.vss(idx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0; %zeros(DAE.n + 2*Bus.n, 1);
    DAE.Ac(k,:) = 0; %zeros(1, DAE.n + 2*Bus.n);
    DAE.Ac(k,k) = speye(length(idx));
  end

end

% ----------------------------------------------------------------------------
% function for creating warning messages
function psswarn(idx, msg)
global Pss
fm_disp(strcat('Warning: PSS #',int2str(idx),' at bus #',int2str(Pss.bus(idx)),msg))

% ----------------------------------------------------------------------------
% function for computing Jacobians due to active power of Synchronous Machines
function pssjac(m,k,A)
global Bus DAE Syn

DAE.Fy(k,Syn.bus(m)) = DAE.Fy(k,Syn.bus(m)) + A*Syn.J11(m);
DAE.Fy(k,Syn.bus(m)+Bus.n) = DAE.Fy(k,Syn.bus(m)+Bus.n) + A*Syn.J12(m);
DAE.Fx(k,Syn.delta(m)) = DAE.Fx(k,Syn.delta(m)) + A*Syn.Gp(m,1);
switch Syn.con(m,5)
 case 3
  DAE.Fx(k,Syn.e1q(m))  = DAE.Fx(k,Syn.e1q(m)) + A*Syn.Gp(m,3);
 case 4
  DAE.Fx(k,Syn.e1q(m))  = DAE.Fx(k,Syn.e1q(m)) + A*Syn.Gp(m,3);
  DAE.Fx(k,Syn.e1d(m))  = DAE.Fx(k,Syn.e1d(m)) + A*Syn.Gp(m,4);
 case 5.1
  DAE.Fx(k,Syn.e1q(m))  = DAE.Fx(k,Syn.e1q(m)) + A*Syn.Gp(m,3);
  DAE.Fx(k,Syn.e2d(m))  = DAE.Fx(k,Syn.e2d(m)) + A*Syn.Gp(m,6);
 case 5.2
  DAE.Fx(k,Syn.e2q(m))  = DAE.Fx(k,Syn.e2q(m)) + A*Syn.Gp(m,5);
  DAE.Fx(k,Syn.e2d(m))  = DAE.Fx(k,Syn.e2d(m)) + A*Syn.Gp(m,6);
 case 5.3
  DAE.Fx(k,Syn.e1q(m))  = DAE.Fx(k,Syn.e1q(m)) + A*Syn.Gp(m,3);
  DAE.Fx(k,Syn.psiq(m)) = DAE.Fx(k,Syn.psiq(m)) + A*Syn.Gp(m,7);
  DAE.Fx(k,Syn.psid(m)) = DAE.Fx(k,Syn.psid(m)) + A*Syn.Gp(m,8);
 case 6
  DAE.Fx(k,Syn.e2q(m))  = DAE.Fx(k,Syn.e2q(m)) + A*Syn.Gp(m,5);
  DAE.Fx(k,Syn.e2d(m))  = DAE.Fx(k,Syn.e2d(m)) + A*Syn.Gp(m,6);
 case 8
  DAE.Fx(k,Syn.e2q(m))  = DAE.Fx(k,Syn.e2q(m)) + A*Syn.Gp(m,5);
  DAE.Fx(k,Syn.e2d(m))  = DAE.Fx(k,Syn.e2d(m)) + A*Syn.Gp(m,6);
  DAE.Fx(k,Syn.psiq(m)) = DAE.Fx(k,Syn.psiq(m)) + A*Syn.Gp(m,7);
  DAE.Fx(k,Syn.psid(m)) = DAE.Fx(k,Syn.psid(m)) + A*Syn.Gp(m,8);
end

% ----------------------------------------------------------------------------
% function for computing Jacobians due to active the Vss signal in
% AVR equations
function vssexc(h,k)
global DAE

for i = 1:length(h)
  jac = DAE.Fx(:,k(i));
  jac(k(i)) = 0;
  idx = find(jac);
  if ~isempty(idx)
    DAE.Fx(idx,h(i)) = -jac(idx);
  end
end