fm_sssc.m

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

function fm_sssc(flag)
% FM_SSSC define Static Synchronous Series Compensators - SSSC
%
% FM_SSSC(FLAG)
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%       FLAG = 5 non-windup limits
%
%Author:    Hugo M. Ayres (revised by Federico Milano)
%Date:      18-Feb-2006
%Version:   1.0.0
%
%E-mail:    [e-mail]
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2006 Hugo M. Ayres & Federico Milano

global Sssc Syn Pod Bus DAE PV Line jay

bus1 = Sssc.bus1;
bus2 = Sssc.bus2;
V1 = DAE.V(bus1);
V2 = DAE.V(bus2);
a1 = DAE.a(bus1);
a2 = DAE.a(bus2);
ss = sin(a1-a2);
cc = cos(a1-a2);
den = max(sqrt(V1.^2+V2.^2-2.*V1.*V2.*cc),1e-6*ones(Sssc.n,1));
kp = Sssc.Cp./(1-Sssc.Cp);
vcs = DAE.x(Sssc.vcs);
vpi = DAE.x(Sssc.vpi);
Tr = Sssc.con(:,7);
vcs_max = Sssc.con(:,8);
vcs_min = Sssc.con(:,9);

switch flag

 case 0
  
  ty3 = Sssc.con(:,2) == 3;
  typwr = Sssc.con(:,10) == 1;
  tyang = Sssc.con(:,10) == 2;

  DAE.x(Sssc.vcs) = kp.*den;
  DAE.x(Sssc.vpi) = DAE.x(Sssc.vcs(find(ty3)));

  % reset transmission line reactance and admittance matrix
  Line.con(Sssc.line,9) = Line.con(Sssc.line,9) + Sssc.xcs;
  fm_y;

  [Ps0,Qs0,Pr0,Qr0,Is0,Ir0] = fex_lineflows;
  Sssc.Pref = Ps0(Sssc.line) + ty3.*tyang.*Pr0(Sssc.line);

  idx = find(DAE.x(Sssc.vcs) > vcs_max);
  if idx, ssscwarn(idx,' Vs is over its max limit.'), end
  idx = find(DAE.x(Sssc.vcs) < vcs_min);
  if idx, ssscwarn(idx,' Vs is under its min limit.'), end
  DAE.x(Sssc.vcs) = max(DAE.x(Sssc.vcs),vcs_min);
  DAE.x(Sssc.vcs) = min(DAE.x(Sssc.vcs),vcs_max);
  
  % reference voltage signal
  Sssc.V0 = DAE.x(Sssc.vcs);
  
  fm_disp('Initialization of SSSC completed.')
  
 case 1 % algebraic equations
   
  c1 = vcs.*Sssc.y./den;  
  P1 = c1.*V1.*V2.*ss;  
  Q1 = c1.*(V1.^2-V1.*V2.*cc);
  Q2 = c1.*(V2.^2-V1.*V2.*cc);
  
  DAE.gp = DAE.gp + sparse(bus1,1,P1,Bus.n,1) - sparse(bus2,1,P1,Bus.n,1);
  DAE.gq = DAE.gq + sparse(bus1,1,Q1,Bus.n,1) + sparse(bus2,1,Q2,Bus.n,1);  
  
 case 2 % Jacobians of active & reactive powers
     
  c1 = vcs.*Sssc.y./den;
  c2 = vcs.*Sssc.y./(den.^3);  
  M2 = c2.*V1.*V2.*ss;  
  M1 = -M2;
  M3 = c2.*(V2.*cc-V1);
  M4 = c2.*(V1.*cc-V2);
  
  P1a1 = V1.*V2.*(M1.*ss+c1.*cc);
  P1a2 = V1.*V2.*(M2.*ss-c1.*cc);
  P1v1 = V2.*ss.*(M3.*V1+c1);
  P1v2 = V1.*ss.*(M4.*V2+c1);

  Q1a1 = -V1.*V2.*(M1.*cc-c1.*ss)+M1.*V1.^2;
  Q1a2 = -V1.*V2.*(M2.*cc+c1.*ss)+M2.*V1.^2;
  Q2a1 = -V1.*V2.*(M1.*cc-c1.*ss)+M1.*V2.^2;  
  Q2a2 = -V1.*V2.*(M2.*cc+c1.*ss)+M2.*V2.^2;  

  Q1v1 = M3.*V1.*(V1-V2.*cc)+c1.*(2.*V1-V2.*cc);
  Q1v2 = M4.*V1.*(V1-V2.*cc)-c1.*V1.*cc;
  Q2v1 = M3.*V2.*(V2-V1.*cc)-c1.*V2.*cc;
  Q2v2 = M4.*V2.*(V2-V1.*cc)+c1.*(2.*V2-V1.*cc);

  DAE.J11 = DAE.J11 ...
            + sparse(bus1,bus1,P1a1,Bus.n,Bus.n) ...
            + sparse(bus1,bus2,P1a2,Bus.n,Bus.n) ...
            - sparse(bus2,bus1,P1a1,Bus.n,Bus.n) ...
            - sparse(bus2,bus2,P1a2,Bus.n,Bus.n);
  
  DAE.J12 = DAE.J12 ...
            + sparse(bus1,bus1,P1v1,Bus.n,Bus.n) ...
            + sparse(bus1,bus2,P1v2,Bus.n,Bus.n) ...
            - sparse(bus2,bus1,P1v1,Bus.n,Bus.n) ...
            - sparse(bus2,bus2,P1v2,Bus.n,Bus.n);  
  
  DAE.J21 = DAE.J21 ...
            + sparse(bus1,bus1,Q1a1,Bus.n,Bus.n) ...
            + sparse(bus1,bus2,Q1a2,Bus.n,Bus.n) ...
            + sparse(bus2,bus1,Q2a1,Bus.n,Bus.n) ...
            + sparse(bus2,bus2,Q2a2,Bus.n,Bus.n);   
  
  DAE.J22 = DAE.J22 ...
            + sparse(bus1,bus1,Q1v1,Bus.n,Bus.n) ...
            + sparse(bus1,bus2,Q1v2,Bus.n,Bus.n) ...
            + sparse(bus2,bus1,Q2v1,Bus.n,Bus.n) ...
            + sparse(bus2,bus2,Q2v2,Bus.n,Bus.n);
 
 case 3 % differential equations

  ty2 = find(Sssc.con(:,2) == 2);
  ty3 = find(Sssc.con(:,2) == 3);
  
  V0 = Sssc.V0;
  if ty2, V0(ty2) = kp(ty2).*den(ty2); end

  if ty3    
    [Ps,Qs,Pr,Qr,Is,Ir] = fex_lineflows;
    Kpr = Sssc.con(:,11);
    Kin = Sssc.con(:,12);
    tp = ty3(find(Sssc.con(ty3,10) == 1));
    ta = ty3(find(Sssc.con(ty3,10) == 2));    
    if tp
      V0(tp) = Kpr(tp).*(Sssc.Pref(tp)-Ps(Sssc.line(tp))) + DAE.x(Sssc.vpi(tp));
      DAE.f(Sssc.vpi(tp)) = Kin(tp).*(Sssc.Pref(tp)-Ps(Sssc.line(tp)));
    end
    if ta
      V0(ta) = Kpr(ta).*(Sssc.Pref(ta)-Ps(Sssc.line(ta))-Pr(Sssc.line(ta))) ...
               + DAE.x(Sssc.vpi(ta));
      DAE.f(Sssc.vpi(ta)) = Kin(ta).*(Sssc.Pref(ta)-Ps(Sssc.line(ta))-Pr(Sssc.line(ta)));
    end
  end

  % Updating signal V0
  V0(Sssc.pod) = V0(Sssc.pod) + DAE.x(Pod.Vs(Pod.sssc));

  u = ~(vcs >= vcs_max & DAE.f(Sssc.vcs) > 0) & ...
      ~(vcs <= vcs_min & DAE.f(Sssc.vcs) < 0);
  DAE.f(Sssc.vcs) = u.*(V0-vcs)./Tr;
  DAE.x(Sssc.vcs) = min(vcs_max,DAE.x(Sssc.vcs));
  DAE.x(Sssc.vcs) = max(vcs_min,DAE.x(Sssc.vcs));

 case 4 % Jacobians of state variables
 
  Ms = Sssc.y./den;
  
  P1vs = Ms.*V1.*V2.*ss;  
  P2vs = -P1vs;
  Q1vs = Ms.*V1.*(V1-V2.*cc);
  Q2vs = Ms.*V2.*(V2-V1.*cc);
  
  DAE.Fx = DAE.Fx - sparse(Sssc.vcs,Sssc.vcs,1./Tr,DAE.n,DAE.n);
  
  u = (vcs <= vcs_max & vcs >= vcs_min);
  ty2 = Sssc.con(:,2) == 2;
  ty3 = Sssc.con(:,2) == 3;

  if Sssc.pod
    DAE.Fx(Sssc.vcs(Sssc.pod),Pod.Vs(Pod.sssc)) = Pod.u(Pod.sssc).*u(Sssc.pod)./Tr(Sssc.pod);
  end
    
  DAE.Gx = DAE.Gx + sparse(bus1,Sssc.vcs,u.*P1vs,2*Bus.n,DAE.n);  
  DAE.Gx = DAE.Gx + sparse(bus2,Sssc.vcs,u.*P2vs,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(bus1+Bus.n,Sssc.vcs,u.*Q1vs,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(bus2+Bus.n,Sssc.vcs,u.*Q2vs,2*Bus.n,DAE.n);    
  
  a = find(u.*ty2);
  if a
    k = kp(a)./den(a)./Tr(a);
    F1 = k.*(V1(a).*V2(a).*ss(a));   
    F3 = k.*(V1(a)-V2(a).*cc(a));
    F4 = k.*(V2(a)-V1(a).*cc(a));     
    DAE.Fy = DAE.Fy + sparse(Sssc.vcs(a),bus1(a),F1,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy - sparse(Sssc.vcs(a),bus2(a),F1,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(Sssc.vcs(a),bus1(a)+Bus.n,F3,DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy + sparse(Sssc.vcs(a),bus2(a)+Bus.n,F4,DAE.n,2*Bus.n);
  end

  a = find(u.*ty3);
  if a

    Kpr = Sssc.con(a,11);
    Kin = Sssc.con(a,12);
    c1 = Sssc.y(a).*(1+vcs(a)./den(a));
    c2 = Sssc.y(a).*vcs(a)./(den(a).^3);    
    Ms = Sssc.y(a)./den(a);
    M2 = c2.*V1(a).*V2(a).*ss(a);
    M1 = -M2;
    M3 = c2.*(V2(a).*cc(a)-V1(a));
    M4 = c2.*(V1(a).*cc(a)-V2(a));              
    Jps1 = V1(a).*V2(a).*(M1.*ss(a)+c1.*cc(a));
    Jps2 = V1(a).*V2(a).*(M2.*ss(a)-c1.*cc(a));
    Jps3 = V2(a).*ss(a).*(M3.*V1(a)+c1);
    Jps4 = V1(a).*ss(a).*(M4.*V2(a)+c1);
    Jps5 = Ms.*V1(a).*V2(a).*ss(a);

    DAE.Fx = DAE.Fx - sparse(Sssc.vcs(a),Sssc.vcs(a),Jps5.*Kpr(a)./Tr(a),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx + sparse(Sssc.vcs(a),Sssc.vpi(a),1./Tr(a),DAE.n,DAE.n);
    DAE.Fx = DAE.Fx - sparse(Sssc.vpi(a),Sssc.vcs(a),Jps5.*Kin(a),DAE.n,DAE.n);

    DAE.Fy = DAE.Fy - sparse(Sssc.vcs(a),bus1(a),Jps1.*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy - sparse(Sssc.vcs(a),bus2(a),Jps2.*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy - sparse(Sssc.vcs(a),bus1(a)+Bus.n,Jps3.*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy - sparse(Sssc.vcs(a),bus2(a)+Bus.n,Jps4.*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
    
    DAE.Fy = DAE.Fy - sparse(Sssc.vpi(a),bus1(a),Jps1.*Kin(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy - sparse(Sssc.vpi(a),bus2(a),Jps2.*Kin(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy - sparse(Sssc.vpi(a),bus1(a)+Bus.n,Jps3.*Kin(a),DAE.n,2*Bus.n);
    DAE.Fy = DAE.Fy - sparse(Sssc.vpi(a),bus2(a)+Bus.n,Jps4.*Kin(a),DAE.n,2*Bus.n);
    
    c = find(Sssc.con(a,10) == 2);      
    if c
      a = a(c);
    
      DAE.Fy = DAE.Fy + sparse(Sssc.vcs(a),bus1(a),Jps1(c).*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.vcs(a),bus2(a),Jps2(c).*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.vcs(a),bus1(a)+Bus.n,Jps3(c).*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.vcs(a),bus2(a)+Bus.n,Jps4(c).*Kpr(a)./Tr(a),DAE.n,2*Bus.n);
      
      DAE.Fy = DAE.Fy + sparse(Sssc.vpi(a),bus1(a), Jps1(c).*Kin(a),DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.vpi(a),bus2(a), Jps2(c).*Kin(a),DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.vpi(a),bus1(a)+Bus.n,Jps3(c).*Kin(a),DAE.n,2*Bus.n);
      DAE.Fy = DAE.Fy + sparse(Sssc.vpi(a),bus2(a)+Bus.n,Jps4(c).*Kin(a),DAE.n,2*Bus.n);

    end
    
  end    
  
 case 5 % non-windup limiters

  a = find((vcs >= vcs_max | vcs <= vcs_min) & DAE.f(Sssc.vcs) == 0);
  if ~isempty(a)
    DAE.tn(Sssc.vcs(a)) = 0;
    DAE.Ac(Sssc.vcs(a),:) = 0;
    DAE.Ac(:,Sssc.vcs(a)) = 0;
    DAE.Ac(Sssc.vcs(a),Sssc.vcs(a)) = speye(length(a));
  end

end  

% function for creating warning messages
function ssscwarn(idx, msg)
global Sssc Varname
fm_disp(strcat('Warning: SSSC #',int2str(idx), ...
	       ' at bus <',Varname.bus(Sssc.bus1(idx)), ...
	       '>: ',msg))