📄 fm_pss.m
字号:
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: fmilano@thunderbox.uwaterloo.ca
%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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -