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:    Federico Milano
%Date:      11-Nov-2002
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2005 Federico Milano
%
% This toolbox is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2.0 of the License, or
% (at your option) any later version.
%
% This toolbox is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANDABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this toolbox; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307,
% USA.

global Sssc Bus DAE jay

Vdc  = DAE.x(Sssc.Vdc);
beta = DAE.x(Sssc.beta);
mse  = DAE.x(Sssc.m);
phi  = DAE.x(Sssc.Pmac);
Vmdc = DAE.x(Sssc.Vmdc);

bus1 = Sssc.bus1;
bus2 = Sssc.bus2;

V1 = DAE.V(bus1);
theta1 = DAE.a(bus1);
V2 = DAE.V(bus2);
theta2 = DAE.a(bus2);

Rdc = Sssc.con(:,7);
Cdc = Sssc.con(:,8);
Rac = Sssc.con(:,9);
Xac = Sssc.con(:,10);
K  = Sssc.con(:,11);
Kd = Sssc.con(:,12);
T1 = Sssc.con(:,13);
T2 = Sssc.con(:,14);
Kp = Sssc.con(:,15);
Ki = Sssc.con(:,16);
Kmac = Sssc.con(:,17);
Tmac = Sssc.con(:,18);
Kmdc = Sssc.con(:,19);
Tmdc = Sssc.con(:,20);
Vrefdc = Sssc.con(:,21);
phi_ref = Sssc.con(:,22);
Rt = Sssc.con(:,23);
Xt = Sssc.con(:,24);
b_max = Sssc.con(:,25);
b_min = Sssc.con(:,26);
m_max = Sssc.con(:,27);
m_min = Sssc.con(:,28);
c_type = Sssc.con(:,3);

a0r = Sssc.dat(:,1);
a1r = Sssc.dat(:,2);
a2r = Sssc.dat(:,3);
a0i = Sssc.dat(:,4);
a1i = Sssc.dat(:,5);
a2i = Sssc.dat(:,6);

switch flag

 case 0 % Initialization

  a0 = 1./(Rt+Rac + jay*(Xt+Xac));
  a1 = (Rac+jay*Xac).*a0;
  a2 = (Rt+jay*Xt).*a0;

  Sssc.dat(:,1) = real(a0);
  Sssc.dat(:,2) = real(a1);
  Sssc.dat(:,3) = real(a2);
  Sssc.dat(:,4) = imag(a0);
  Sssc.dat(:,5) = imag(a1);
  Sssc.dat(:,6) = imag(a2);

  idx = find(c_type == 2);
  if ~isempty(idx)
    Sssc.con(idx,22) = Sssc.con(idx,22).^2;
  end
  idx = find(c_type == 1);
  if ~isempty(idx)
    Sssc.con(idx,11) = sign(Sssc.con(idx,22)).*Sssc.con(idx,11);
  end
  Sssc.con(:,15) = -abs(Sssc.con(:,15));
  Sssc.con(:,16) = -abs(Sssc.con(:,16));

 case 1 % Algebraic equations

  t1 = cos(theta1);
  t2 = V1.*t1;
  t3 = cos(theta2);
  t4 = V2.*t3;
  t5 = sqrt(6);
  t6 = t5.*mse;
  t7 = cos(beta);
  t11 = t2-t4-t6.*Vdc.*t7./4;
  t13 = sin(theta1);
  t14 = V1.*t13;
  t15 = sin(theta2);
  t16 = V2.*t15;
  t17 = sin(beta);
  t21 = t14-t16-t6.*Vdc.*t17./4;
  t23 = a0r.*t11-a0i.*t21;
  t27 = a0i.*t11+a0r.*t21;

  DAE.gp = DAE.gp + sparse(bus1,1,t2.*t23+t14.*t27,Bus.n,1);
  DAE.gq = DAE.gq + sparse(bus1,1,t14.*t23-t2.*t27,Bus.n,1);
  DAE.gp = DAE.gp + sparse(bus2,1,-t4.*t23-t16.*t27,Bus.n,1);
  DAE.gq = DAE.gq + sparse(bus2,1,-t16.*t23+t4.*t27,Bus.n,1);

 case 2  % Algebraic Jacobians

  t1 = sin(theta1);
  t2 = V1.*t1;
  t3 = cos(theta1);
  t4 = V1.*t3;
  t5 = cos(theta2);
  t6 = V2.*t5;
  t7 = sqrt(6);
  t8 = t7.*mse;
  t9 = cos(beta);
  t13 = t4-t6-t8.*Vdc.*t9./4;
  t15 = sin(theta2);
  t16 = V2.*t15;
  t17 = sin(beta);
  t21 = t2-t16-t8.*Vdc.*t17./4;
  t23 = a0r.*t13-a0i.*t21;
  t25 = a0r.*V1;
  t27 = a0i.*V1;
  t29 = -t25.*t1-t27.*t3;
  t33 = a0i.*t13+a0r.*t21;
  t37 = -t27.*t1+t25.*t3;
  t43 = a0r.*t3-a0i.*t1;
  t48 = a0i.*t3+a0r.*t1;
  t51 = a0r.*V2;
  t53 = a0i.*V2;
  t55 = t51.*t15+t53.*t5;
  t59 = t53.*t15-t51.*t5;
  t64 = -a0r.*t5+a0i.*t15;
  t68 = -a0i.*t5-a0r.*t15;

  DAE.J11 = DAE.J11 + sparse(bus1,bus1,-t2.*t23+t4.*t29+t4.*t33+t2.*t37,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(bus1,bus1,t3.*t23+t4.*t43+t1.*t33+t2.*t48,Bus.n,Bus.n);
  DAE.J11 = DAE.J11 + sparse(bus1,bus2,t4.*t55+t2.*t59,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(bus1,bus2,t4.*t64+t2.*t68,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(bus1,bus1,t4.*t23+t2.*t29+t2.*t33-t4.*t37,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(bus1,bus1,t1.*t23+t2.*t43-t3.*t33-t4.*t48,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(bus1,bus2,t2.*t55-t4.*t59,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(bus1,bus2,t2.*t64-t4.*t68,Bus.n,Bus.n);
  DAE.J11 = DAE.J11 + sparse(bus2,bus1,-t6.*t29-t16.*t37,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(bus2,bus1,-t6.*t43-t16.*t48,Bus.n,Bus.n);
  DAE.J11 = DAE.J11 + sparse(bus2,bus2,t16.*t23-t6.*t55-t6.*t33-t16.*t59,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 + sparse(bus2,bus2,-t5.*t23-t6.*t64-t15.*t33-t16.*t68,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(bus2,bus1,-t16.*t29+t6.*t37,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(bus2,bus1,-t16.*t43+t6.*t48,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(bus2,bus2,-t6.*t23-t16.*t55-t16.*t33+t6.*t59,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(bus2,bus2,-t15.*t23-t16.*t64+t5.*t33+t6.*t68,Bus.n,Bus.n);

 case 3 % Differential Equations

  t1 = cos(theta1);
  t2 = V1.*t1;
  t3 = cos(theta2);
  t4 = V2.*t3;
  t5 = t2-t4;
  t7 = sin(theta1);
  t8 = V1.*t7;
  t9 = sin(theta2);
  t10 = V2.*t9;
  t11 = t8-t10;
  t13 = sqrt(6);
  t14 = a2r.*t13;
  t15 = mse.*Vdc;
  t16 = cos(beta);
  t17 = t15.*t16;
  t20 = a2i.*t13;
  t21 = sin(beta);
  t22 = t15.*t21;
  t26 = t13.*mse;
  t30 = t2-t4-t26.*Vdc.*t16./4;
  t35 = t8-t10-t26.*Vdc.*t21./4;
  t37 = a0r.*t30-a0i.*t35;
  t48 = a0i.*t30+a0r.*t35;
  t51 = 1./Cdc;
  t53 = 1./Vdc;
  t58 = t37.^2;
  t59 = t48.^2;
  t60 = t58+t59;
  t70 = 1./Tmac;
  t71 = (-phi+Kmac.*(t2.*t37+t8.*t48)).*t70;
  t72 = Kd.*mse;
  t74 = K.*(phi_ref-phi);
  t75 = K.*T1;
  t78 = 1./T2;
  t82 = (-phi+Kmac.*t60).*t70;
  t88 = 1./Tmdc;

  DAE.f(Sssc.Vdc) = ((a1r.*t5-a1i.*t11+t14.*t17./4-t20.*t22./4).* ...
                     t37+(a1i.*t5+a1r.*t11+t14.*t22./4+t20.*t17./4).* ...
                     t48).*t51.*t53-Vdc./Rdc.*t51-Rac.*t60.*t51.*t53;

  c1 = find(c_type == 1);
  if ~isempty(c1)
    DAE.f(Sssc.Pmac(c1)) = t71(c1);
    DAE.f(Sssc.m(c1)) = (-t72(c1)+t74(c1)-t75(c1).*t71(c1)).*t78(c1);
  end

  c2 = find(c_type == 2);
  if ~isempty(c2)
    DAE.f(Sssc.Pmac(c2)) = t82(c2);
    DAE.f(Sssc.m(c2)) = (-t72(c2)+t74(c2)-t75(c2).*t82(c2)).*t78(c2);
  end

  DAE.f(Sssc.Vmdc) = (-Vmdc+Kmdc.*Vdc).*t88;
  DAE.f(Sssc.beta) = (Kp.*t88-Ki).*Vmdc+Ki.*Vrefdc-Kp.*Kmdc.*t88.* ...
      Vdc;

  % m limits
  idx = find(mse <= m_min & DAE.f(Sssc.m) < 0);
  if idx, DAE.f(Sssc.beta(idx)) = 0; end
  DAE.x(Sssc.m) = max(DAE.x(Sssc.m),m_min);
  idx = find(mse >= m_max & DAE.f(Sssc.m) > 0);
  if idx, DAE.f(Sssc.m(idx)) = 0; end
  DAE.x(Sssc.m) = min(DAE.x(Sssc.m),m_max);

  % beta limits
  idx = find(beta <= b_min & DAE.f(Sssc.beta) < 0);
  if idx, DAE.f(Sssc.beta(idx)) = 0; end
  DAE.x(Sssc.beta) = max(DAE.x(Sssc.beta),b_min);
  idx = find(beta >= b_max & DAE.f(Sssc.beta) > 0);
  if idx, DAE.f(Sssc.beta(idx)) = 0; end
  DAE.x(Sssc.beta) = min(DAE.x(Sssc.beta),b_max);

 case 4 % Jacobians of state and service variables

  t1 = sqrt(6);
  t2 = a2r.*t1;
  t3 = cos(beta);
  t4 = mse.*t3;
  t6 = a2i.*t1;
  t7 = sin(beta);
  t8 = t7.*mse;
  t11 = cos(theta1);
  t12 = V1.*t11;
  t13 = cos(theta2);
  t14 = V2.*t13;
  t15 = t1.*mse;
  t16 = Vdc.*t3;
  t19 = t12-t14-t15.*t16./4;
  t21 = sin(theta1);
  t22 = V1.*t21;
  t23 = sin(theta2);
  t24 = V2.*t23;
  t25 = Vdc.*t7;
  t28 = t22-t24-t15.*t25./4;
  t30 = a0r.*t19-a0i.*t28;
  t32 = t12-t14;
  t34 = t22-t24;