📄 fm_pod.m
字号:
function fm_pod(flag)
% FM_POD defines a Supplementary Stabilizing Control Loop for SVCs
%
% FM_POD(FLAG)
% FLAG = 0 initialization
% FLAG = 3 differential equations
% FLAG = 4 state matrix
%
global Pod DAE Bus Settings Line
% Define input signal indices
VSI = zeros(Pod.n,1);
SIv = find(Pod.type == 1); % V
SIPs = find(Pod.type == 2); % Pij
SIPr = find(Pod.type == 3); % Pji
SIIs = find(Pod.type == 4); % Iij
SIIr = find(Pod.type == 5); % Iji
SIQs = find(Pod.type == 6); % Qij
SIQr = find(Pod.type == 7); % Qji
if SIPs, [Ps, JPs] = fm_pflows(1,Pod.idx(SIPs),flag); end
if SIPr, [Pr, JPr] = fm_pflows(2,Pod.idx(SIPr),flag); end
if SIIs, [Is, JIs] = fm_pflows(3,Pod.idx(SIIs),flag); end
if SIIr, [Ir, JIr] = fm_pflows(4,Pod.idx(SIIr),flag); end
if SIQs, [Qs, JQs] = fm_pflows(5,Pod.idx(SIQs),flag); end
if SIQr, [Qr, JQr] = fm_pflows(6,Pod.idx(SIQr),flag); end
if SIv , VSI(SIv) = DAE.V(Pod.idx(SIv)); end
if SIPs, VSI(SIPs) = Ps; end
if SIPr, VSI(SIPr) = Pr; end
if SIIs, VSI(SIIs) = Is; end
if SIIr, VSI(SIIr) = Ir; end
if SIQs, VSI(SIQs) = Qs; end
if SIQr, VSI(SIQr) = Qr; end
vsmax = Pod.con(:,5);
vsmin = Pod.con(:,6);
Kw = Pod.con(:,7);
Tw = Pod.con(:,8);
T1 = Pod.con(:,9);
T2 = Pod.con(:,10);
T3 = Pod.con(:,11);
T4 = Pod.con(:,12);
Tr = Pod.con(:,13);
switch flag
case 0
idx = find(Tw == 0);
if idx,
Pod.con(idx,7) = 1;
podwarn(idx,' Tw cannot be zero. Default value Tw = 1 will be used.')
end
idx = find(T2 == 0);
if idx,
Pod.con(idx,10) = 0.01;
podwarn(idx,' T2 cannot be zero. Default value T2 = 0.01 will be used.')
end
idx = find(T4 == 0);
if idx,
Pod.con(idx,12) = 0.01;
podwarn(idx,' T4 cannot be zero. Default value T4 = 0.01 will be used.')
end
DAE.x(Pod.v1) = -Kw.*VSI;
DAE.x(Pod.v2) = 0;
DAE.x(Pod.v3) = 0;
DAE.x(Pod.Vs) = 0;
Pod.u = ones(Pod.n,1);
fm_disp('Initialization of Power Oscillation Damper completed.')
case 3
DAE.f(Pod.v1) = -(Kw.*VSI+DAE.x(Pod.v1))./Tw;
A = T1./T2;
B = 1-A;
C = T3./T4;
D = 1-C;
% output of wash-out block
S1 = Kw.*VSI+DAE.x(Pod.v1);
DAE.f(Pod.v2) = (B.*S1 - DAE.x(Pod.v2))./T2;
DAE.f(Pod.v3) = (D.*(DAE.x(Pod.v2)+A.*S1)-DAE.x(Pod.v3))./T4;
VS = DAE.x(Pod.v3)+C.*(DAE.x(Pod.v2)+A.*S1);
DAE.f(Pod.Vs) = (VS-DAE.x(Pod.Vs))./Tr;
% non-windup limits
idx = find(VS >= vsmax & DAE.f(Pod.Vs) > 0);
if idx, DAE.f(Pod.Vs(idx)) = 0; end
idx = find(VS <= vsmin & DAE.f(Pod.Vs) < 0);
if idx, DAE.f(Pod.Vs(idx)) = 0; end
DAE.x(Pod.Vs) = min(DAE.x(Pod.Vs),vsmax);
DAE.x(Pod.Vs) = max(DAE.x(Pod.Vs),vsmin);
Pod.u = DAE.x(Pod.Vs) <= vsmax & DAE.x(Pod.Vs) >= vsmin;
case 4
A = T1./T2;
B = 1-A;
C = T3./T4;
D = 1-C;
E = C.*A;
F = D./T4;
G = B./T2;
DAE.Fx = DAE.Fx - sparse(Pod.v1,Pod.v1,1./Tw,DAE.n,DAE.n); % df1/dv1
DAE.Fx = DAE.Fx + sparse(Pod.v2,Pod.v1,G,DAE.n,DAE.n); % df2/dv1
DAE.Fx = DAE.Fx - sparse(Pod.v2,Pod.v2,1./T2,DAE.n,DAE.n); % df2/dv2
DAE.Fx = DAE.Fx + sparse(Pod.v3,Pod.v1,F.*A,DAE.n,DAE.n); % df3/dv1
DAE.Fx = DAE.Fx + sparse(Pod.v3,Pod.v2,F,DAE.n,DAE.n); % df3/dv2
DAE.Fx = DAE.Fx - sparse(Pod.v3,Pod.v3,1./T4,DAE.n,DAE.n); % df3/dv3
DAE.Fx = DAE.Fx + sparse(Pod.Vs,Pod.v1,Pod.u.*E./Tr,DAE.n,DAE.n); % df4/dv1
DAE.Fx = DAE.Fx + sparse(Pod.Vs,Pod.v2,Pod.u.*C./Tr,DAE.n,DAE.n); % df4/dv2
DAE.Fx = DAE.Fx + sparse(Pod.Vs,Pod.v3,Pod.u./Tr,DAE.n,DAE.n); % df4/dv3
DAE.Fx = DAE.Fx - sparse(Pod.Vs,Pod.Vs,1./Tr,DAE.n,DAE.n); % df4/dVs
% Algebraic Jacobians
K1 = -Kw./Tw;
K2 = Kw.*G;
K3 = Kw.*D.*A./T4;
K4 = Pod.u.*Kw.*E./Tr;
if SIv
% Jacobians for Bus voltage input signals
DAE.Fy = DAE.Fy + sparse(Pod.v1,Pod.con(SIv,1)+Bus.n,K1,DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v2,Pod.con(SIv,1)+Bus.n,K2,DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v3,Pod.con(SIv,1)+Bus.n,K3,DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.Vs,Pod.con(SIv,1)+Bus.n,K4,DAE.n,2*Bus.n);
end
S = find(Pod.type > 1);
if ~isempty(S)
J = zeros(Pod.n,4);
if SIPs, J(SIPs,:) = JPs; end
if SIPr, J(SIPr,:) = JPr; end
if SIQs, J(SIQs,:) = JQs; end
if SIQr, J(SIQr,:) = JQr; end
if SIIs, J(SIIs,:) = JIs; end
if SIIr, J(SIIr,:) = JIr; end
L = Pod.idx(S);
DAE.Fy = DAE.Fy + sparse(Pod.v1(S),Line.from(L),K1(S).*J(S,1),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v1(S),Line.from(L)+Bus.n,K1(S).*J(S,2),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v1(S),Line.to(L),K1(S).*J(S,3),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v1(S),Line.to(L)+Bus.n,K1(S).*J(S,4),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v2(S),Line.from(L),K2(S).*J(S,1),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v2(S),Line.from(L)+Bus.n,K2(S).*J(S,2),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v2(S),Line.to(L),K2(S).*J(S,3),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v2(S),Line.to(L)+Bus.n,K2(S).*J(S,4),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v3(S),Line.from(L),K3(S).*J(S,1),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v3(S),Line.from(L)+Bus.n,K3(S).*J(S,2),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v3(S),Line.to(L),K3(S).*J(S,3),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.v3(S),Line.to(L)+Bus.n,K3(S).*J(S,4),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.Vs(S),Line.from(L),K4(S).*J(S,1),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.Vs(S),Line.from(L)+Bus.n,K4(S).*J(S,2),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.Vs(S),Line.to(L),K4(S).*J(S,3),DAE.n,2*Bus.n);
DAE.Fy = DAE.Fy + sparse(Pod.Vs(S),Line.to(L)+Bus.n,K4(S).*J(S,4),DAE.n,2*Bus.n);
end
case 5 % non-windup limiters
fm_windup(Pod.Vs,Pod.con(:,5),Pod.con(:,6));
end
% ----------------------------------------------------------------------------
% function for creating warning messages
% ----------------------------------------------------------------------------
function podwarn(idx, msg)
global Pod
fm_disp(strcat('Warning: POD #',int2str(idx), ...
' at bus #',int2str(Pod.svc(idx)),msg))
% ----------------------------------------------------------------------------
function [F, J] = fm_pflows(type,pfl,flag)
% ----------------------------------------------------------------------------
% FM_LFLOWS compute current and power in transmission lines
% and the associated Jacobian
%
% NOTE: Flows are Jacobiasn computed by this function are only valid for
% Tansmsmission lines whith tap ratio = 1.
%%
% Reference: PSAT function fm_lines
%
% pfl: lines for wich flows and jacobians are requeried
% nsl=length(pfl)
%
% Fij -> column vector (nsl,1) of flows from bus "i" to bus "j"
% Jij -> Jacobian matrix (nsl,Bus.n) of flows from bus "i" to bus "j"
%
% Fji -> column vector (nsl,1) of flows from bus "j" to bus "j"
% Jji -> Jacobian matrix (nsl,Bus.n) of flows from bus "j" to bus "i"
%
%
global Line Bus DAE jay
% ==========================================================================
% line flow solution
% ==========================================================================
V1 = DAE.V(Line.from(pfl));
V2 = DAE.V(Line.to(pfl));
theta1 = DAE.a(Line.from(pfl));
theta2 = DAE.a(Line.to(pfl));
rl = Line.con(pfl,8);
xl = Line.con(pfl,9);
bl = Line.con(pfl,10);
z = rl + jay*xl;
y = 1./z;
g12 = real(y);
b12 = imag(y);
bl0 = 0.5*bl;
J = zeros(length(pfl),4);
switch type
case 1 % P12
cc = cos(theta1-theta2);
ss = sin(theta1-theta2);
F = V1.^2.*g12 - V1.*V2.*(g12.*cc+b12.*ss);
if flag == 4
J(:,1) = -V1.*V2*(-g12.*ss+b12.*cc); % dP12/d(theta1)
J(:,2) = 2.*V1.*g12-V2.*(g12.*cc+b12.*ss); % dP12/dV1
J(:,3) = -V1.*V2.*(g12.*ss-b12.*cc); % dP12/d(theta2)
J(:,4) = -V1.*(g12.*cc+b12.*ss); % dP12/dV2
end
case 2 % P21
cc = cos(theta1-theta2);
ss = sin(theta1-theta2);
F = V2.^2.*g12 - V1.*V2.*(g12.*cc-b12.*ss);
if flag == 4
J(:,1) = -V1.*V2.*(-g12.*ss-b12.*cc); % dP21/d(theta1)
J(:,2) = -V2.*(g12.*cc-b12.*ss); % dP21/dV1
J(:,3) = 2*V2.*g12-V1.*(g12.*cc-b12.*ss); % dP21/d(theta2)
J(:,4) = -V1.*V2.*(g12.*ss+b12.*cc); % dP21/dV2
end
case 3 % I12
V11 = V1.*exp(jay*theta1);
V22 = V2.*exp(jay*theta2);
I12 = (V11-V22).*y + V11.*jay*bl0; % Is
Ir12=real(I12); Im12=imag(I12);
F = abs(I12);
if flag == 4
cs1 = cos(theta1);
cs2 = cos(theta2);
sn1 = sin(theta1);
sn2 = sin(theta2);
JIr12_1=g12.*cs1-b12.*sn1-bl0.*sn1; % dIr12/dV1
JIr12_2=-V1.*g12.*sn1-V1.*b12.*cs1-V1.*bl0.*cs1; % dIr12/theta1
JIr12_3=-g12.*cs2+b12.*sn2; % dIr12/dV2
JIr12_4=V2.*g12.*sn2+V2.*b12.*cs2; % dIr12/theta2
JIm12_1=b12.*cs1+g12.*sn1+bl0.*cs1; % dIm12/dV1
JIm12_2=-V1.*b12.*sn1+V1.*g12.*cs1-V1.*bl0.*sn1; % dIm12/theta1
JIm12_3=-b12.*cs2-g12.*sn2; % dIm12/dV2
JIm12_4=V2.*b12.*sn2-V2.*g12.*cs2; % dIm12/theta2
J(:,1) = 0.5*(2*Ir12.*JIr12_2+2*Im12.*JIm12_2)./abs(I12); % dI12/theta1
J(:,2) = 0.5*(2*Ir12.*JIr12_1+2*Im12.*JIm12_1)./abs(I12); % dI12/dV1
J(:,3) = 0.5*(2*Ir12.*JIr12_4+2*Im12.*JIm12_4)./abs(I12); % dI12/theta2
J(:,4) = 0.5*(2*Ir12.*JIr12_3+2*Im12.*JIm12_3)./abs(I12); % dI12/dV2
end
case 4 % I21
V11 = V1.*exp(jay*theta1);
V22 = V2.*exp(jay*theta2);
I21 = abs((V22-V11).*y + V22.*jay*bl0); % Ir
Ir12=real(I12); Im12=imag(I12);
F = abs(I21);
if flag == 4
cs1 = cos(theta1);
cs2 = cos(theta2);
sn1 = sin(theta1);
sn2 = sin(theta2);
JIr21_1=-g12.*cs1+b12.*sn1; % dIr21/dV1
JIr21_2=V1.*b12.*cs1+V1.*g12.*sn1; % dIr21/theta1
JIr21_3=g12.*cs2-b12.*sn2-bl0.*sn2; % dIr21/dV2
JIr21_4=-V2.*g12.*sn2-V2.*b12.*cs2-V2.*bl0.*cs2; % dIr21/theta2
JIm21_1=-b12.*cs1-g12.*sn1; % dIm21/dV1
JIm21_2=V1.*b12.*sn1-V1.*g12.*cs1; % dIm21/theta1
JIm21_3=g12.*sn2+b12.*cs2+bl0.*cs2; % dIm21/dV2
JIm21_4=V2.*g12.*cs2-V2.*b12.*sn2-V2.*bl0.*sn2; % dIm21/theta2
J(:,1) = 0.5*(2*Ir21.*JIr21_2+2*Im21.*JIm21_2)./abs(I21); % dI21/theta1
J(:,2) = 0.5*(2*Ir21.*JIr21_1+2*Im21.*JIm21_1)./abs(I21); % dI21/dV1
J(:,3) = 0.5*(2*Ir21.*JIr21_4+2*Im21.*JIm21_4)./abs(I21); % dI21/theta2
J(:,4) = 0.5*(2*Ir21.*JIr21_3+2*Im21.*JIm21_3)./abs(I21); % dI21/dV2
end
case 5 % Q12
cc = cos(theta1-theta2);
ss = sin(theta1-theta2);
F = -V1.^2.*(b12+bl0)-V1.*V2.*(g12.*ss-b12.*cc);
if flag == 4;
J(:,1) = -V1.*V2.*(g12.*cc+b12.*ss); % dQ12/d(theta1)
J(:,2) = -2.*V1.*(b12+bl0)-V2.*(g12.*ss-b12.*cc); % dQ12/dV1
J(:,3) = -V1.*V2.*(-g12.*cc-b12.*ss); % dQ12/d(theta2)
J(:,4) = -V1.*(g12.*ss-b12.*cc); % dQ12/dV2
end
case 6 % Q21
cc = cos(theta1-theta2);
ss = sin(theta1-theta2);
F = -V2.^2.*(b12+bl0)+V1.*V2.*(g12.*ss+b12.*cc);
if flag == 4;
J(:,1) = V1.*V2.*(g12.*cc-b12.*ss); % dQ12/d(theta1)
J(:,2) = V2.*(g12.*ss+b12.*cc); % dQ12/dV1
J(:,3) = V1.*V2.*(-g12.*cc+b12.*ss); % dQ12/d(theta2)
J(:,4) = -2*V2.*(b12+bl0)+V1.*(g12*ss+b12.*cc); % dQ12/dV2
end
end
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -