fm_pod.m

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

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