fm_ssr.m

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

function fm_ssr(flag)
% FM_SSR Model of the synchronous machine with shaft dynamics and compensated line
%
% FM_SSR(FLAG)
%       FLAG = 0 initializations
%       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-2006 Federico Milano

global SSR DAE Bus Settings PV jay

id = DAE.x(SSR.Id);
iq = DAE.x(SSR.Iq);
If = DAE.x(SSR.If);
ed = DAE.x(SSR.Edc);
eq = DAE.x(SSR.Eqc);
w1 = DAE.x(SSR.omega_HP);
w2 = DAE.x(SSR.omega_IP);
w3 = DAE.x(SSR.omega_LP);
w4 = DAE.x(SSR.omega);
w5 = DAE.x(SSR.omega_EX);
t1 = DAE.x(SSR.delta_HP);
t2 = DAE.x(SSR.delta_IP);
t3 = DAE.x(SSR.delta_LP);
t4 = DAE.x(SSR.delta);
t5 = DAE.x(SSR.delta_EX);

xd  = SSR.con(:,5);
xq  = SSR.con(:,6);
xad = SSR.con(:,8);
xf  = SSR.con(:,13);

xl  = SSR.con(:,10);
xc  = SSR.con(:,11);
xeq = xad.*xad-xf.*(xd+xl);

rf  = SSR.con(:,12);
ra  = SSR.con(:,7);
r   = SSR.con(:,9);

D1 = SSR.con(:,19);
D2 = SSR.con(:,20);
D3 = SSR.con(:,21);
D4 = SSR.con(:,22);
D5 = SSR.con(:,23);

M1 = SSR.con(:,14);
M2 = SSR.con(:,15);
M3 = SSR.con(:,16);
M4 = SSR.con(:,17);
M5 = SSR.con(:,18);

k12 = SSR.con(:,24);
k23 = SSR.con(:,25);
k34 = SSR.con(:,26);
k45 = SSR.con(:,27);

Tm  = SSR.Tm;
Efd = SSR.Efd;
Wb = Settings.rad;

V = DAE.V(SSR.bus);
a = DAE.a(SSR.bus);
cdt = cos(t4-a);
sdt = sin(t4-a);

switch flag
 case 0 % initialization

  VV =  V.*exp(jay*a);
  S = Bus.Pg(SSR.bus) - jay*Bus.Qg(SSR.bus);
  I = S./conj(VV);
  delta = angle(VV + (ra+r + jay*(xq+xl-xc)).*I);
  cdt = cos(delta-a);
  sdt = sin(delta-a);
  Idq = I.*exp(-jay*(delta-pi/2));
  Id = real(Idq);
  Iq = imag(Idq);
  If = (xl.*Id+(ra+r).*Iq+xd.*Id-xc.*Id+V.*cdt)./xad;
  Te = (xq-xd).*Id.*Iq+xad.*If.*Iq;
  delta_LP = Te./k34+delta;
  delta_IP = ((k23+k34).*delta_LP-k34.*delta)./k23;
  delta_HP = ((k12+k23).*delta_IP-k23.*delta_LP)./k12;
  B = (ra+r).*Id-xl.*Iq-xq.*Iq+xc.*Iq+V.*sdt;

  SSR.Tm = k12.*delta_HP-k12.*delta_IP;
  SSR.Efd = xf.*B./rf + xad.*If;

  DAE.x(SSR.Id) = Id;
  DAE.x(SSR.Iq) = Iq;
  DAE.x(SSR.If) = If;

  DAE.x(SSR.omega_HP) = 1;
  DAE.x(SSR.omega_IP) = 1;
  DAE.x(SSR.omega_LP) = 1;
  DAE.x(SSR.omega) = 1;
  DAE.x(SSR.omega_EX) = 1;

  DAE.x(SSR.Eqc) = -xc.*Id;
  DAE.x(SSR.Edc) =  xc.*Iq;

  DAE.x(SSR.delta) = delta;
  DAE.x(SSR.delta_EX) = delta;
  DAE.x(SSR.delta_HP) = delta_HP;
  DAE.x(SSR.delta_IP) = delta_IP;
  DAE.x(SSR.delta_LP) = delta_LP;

  % find & delete static generators
  check = 1;
  for j = 1:SSR.n
    if ~fm_rmgen(SSR.bus(j)), check = 0; end
  end

  if ~check
    fm_disp('SSR models cannot be properly initialized.')
  else
    fm_disp('Initialization of SSR models completed.')
  end

 case 1 % algebraic equations

  DAE.gp = DAE.gp - sparse(SSR.bus,1,V.*sdt.*id+V.*cdt.*iq,Bus.n,1);
  DAE.gq = DAE.gq - sparse(SSR.bus,1,V.*cdt.*id-V.*sdt.*iq,Bus.n,1);

 case 2 % algebraic Jacobians

  DAE.J11 = DAE.J11 + sparse(SSR.bus,SSR.bus,V.*cdt.*id-V.*sdt.*iq,Bus.n,Bus.n);
  DAE.J12 = DAE.J12 - sparse(SSR.bus,SSR.bus,sdt.*id+cdt.*iq,Bus.n,Bus.n);
  DAE.J21 = DAE.J21 - sparse(SSR.bus,SSR.bus,V.*sdt.*id+V.*cdt.*iq,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(SSR.bus,SSR.bus,sdt.*iq-cdt.*id,Bus.n,Bus.n);

 case 3 % differential equations

  Te = ((xq-xd).*id+xad.*If).*iq;
  B = Wb*((ra+r).*id-(xl+w4.*xq).*iq+ed+V.*sdt);
  C = Wb*(rf.*Efd./xad-rf.*If);

  DAE.f(SSR.Id) = (xf.*B-xad.*C)./xeq;
  DAE.f(SSR.Iq) = -Wb*((xl+w4.*xd).*id+(ra+r).*iq-w4.*xad.*If+eq+V.*cdt)./(xq+xl);
  DAE.f(SSR.If) = (xad.*B-(xd+xl).*C)./xeq;
  DAE.f(SSR.Edc) = Wb*(xc.*id+eq);
  DAE.f(SSR.Eqc) = Wb*(xc.*iq-ed);
  DAE.f(SSR.omega_HP) = (Tm-D1.*(w1-1)+k12.*(t2-t1))./M1;
  DAE.f(SSR.omega_IP) = (-D2.*(w2-1)+k12.*(t1-t2)+k23.*(t3-t2))./M2;
  DAE.f(SSR.omega_LP) = (-D3.*(w3-1)+k23.*(t2-t3)+k34.*(t4-t3))./M3;
  DAE.f(SSR.omega) = (-Te-D4.*(w4-1)+k34.*(t3-t4)+k45.*(t5-t4))./M4;
  DAE.f(SSR.omega_EX) = (-D5.*(w5-1)+k45.*(t4-t5))./M5;

  DAE.f(SSR.delta_HP) = Wb*(w1-1);
  DAE.f(SSR.delta_IP) = Wb*(w2-1);
  DAE.f(SSR.delta_LP) = Wb*(w3-1);
  DAE.f(SSR.delta)    = Wb*(w4-1);
  DAE.f(SSR.delta_EX) = Wb*(w5-1);

 case 4 % Jacobians of state variables

  DAE.Fx = DAE.Fx + sparse(SSR.delta_HP,SSR.omega_HP,Wb,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.delta_IP,SSR.omega_IP,Wb,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.delta_LP,SSR.omega_LP,Wb,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.delta,SSR.omega,Wb,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.delta_EX,SSR.omega_EX,Wb,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx - sparse(SSR.omega_HP,SSR.omega_HP,D1./M1,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.omega_HP,SSR.delta_HP,k12./M1,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.omega_HP,SSR.delta_IP,k12./M1,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx - sparse(SSR.omega_IP,SSR.omega_IP,D2./M2,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.omega_IP,SSR.delta_HP,k12./M2,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.omega_IP,SSR.delta_IP,(k12+k23)./M2,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.omega_IP,SSR.delta_LP,k23./M2,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx - sparse(SSR.omega_LP,SSR.omega_LP,D3./M3,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.omega_LP,SSR.delta_IP,k23/M3,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.omega_LP,SSR.delta_LP,(k23+k34)./M3,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.omega_LP,SSR.delta,k34./M3,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx - sparse(SSR.omega_EX,SSR.omega_EX,D5./M5,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.omega_EX,SSR.delta,k45./M5,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.omega_EX,SSR.delta_EX,k45./M5,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx + sparse(SSR.omega,SSR.delta_LP,k34./M4,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.omega,SSR.delta,(k34+k45)./M4,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.omega,SSR.delta_EX,k45./M4,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.omega,SSR.omega,D4./M4,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.omega,SSR.Id,(xq-xd).*iq./M4,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.omega,SSR.Iq,((xq-xd).*id+xad.*If)./M4,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.omega,SSR.If,xad.*iq./M4,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx + sparse(SSR.Edc,SSR.Id,Wb*xc,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.Edc,SSR.Eqc,Wb,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.Eqc,SSR.Iq,Wb*xc,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.Eqc,SSR.Edc,Wb,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx + sparse(SSR.Id,SSR.Id,Wb*xf.*(ra+r)./xeq,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.Id,SSR.Iq,Wb*xf.*(xl+w4.*xq)./xeq,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.Id,SSR.If,Wb*xad.*rf./xeq,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.Id,SSR.Edc,Wb*xf./xeq,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.Id,SSR.omega,Wb*xf.*iq.*xq./xeq,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.Id,SSR.delta,Wb*V.*xf.*cdt./xeq,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx + sparse(SSR.If,SSR.Id,Wb*xad.*(ra+r)./xeq,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.If,SSR.Iq,Wb*xad.*(xl+w4.*xq)./xeq,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.If,SSR.If,Wb*(xd+xl).*rf./xeq,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.If,SSR.Edc,Wb*xad./xeq,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.If,SSR.omega,Wb*xad.*iq.*xq./xeq,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.If,SSR.delta,Wb*xad.*V.*cdt./xeq,DAE.n,DAE.n);

  DAE.Fx = DAE.Fx - sparse(SSR.Iq,SSR.Id,Wb*(xl+w4.*xd)./(xq+xl),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.Iq,SSR.Iq,Wb*(ra+r)./(xq+xl),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.Iq,SSR.If,Wb*(w4.*xad)./(xq+xl),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx - sparse(SSR.Iq,SSR.Eqc,Wb./(xq+xl),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.Iq,SSR.omega,Wb*(If.*xad-xd.*id)./(xq+xl),DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(SSR.Iq,SSR.delta,Wb*V.*sdt./(xq+xl),DAE.n,DAE.n);

  DAE.Fy = DAE.Fy - sparse(SSR.Id,SSR.bus,Wb*V.*cdt.*xf./xeq,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(SSR.Id,SSR.bus+Bus.n,Wb*sdt.*xf./xeq,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy - sparse(SSR.If,SSR.bus,Wb*V.*cdt.*xad./xeq,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(SSR.If,SSR.bus+Bus.n,Wb*sdt.*xad./xeq,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy - sparse(SSR.Iq,SSR.bus,Wb*V.*sdt./(xq+xl),DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy - sparse(SSR.Iq,SSR.bus+Bus.n,Wb*cdt./(xq+xl),DAE.n,2*Bus.n);

  DAE.Gx = DAE.Gx - sparse(SSR.bus,SSR.Id,V.*sdt,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx - sparse(SSR.bus+Bus.n,SSR.Id,V.*cdt,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx - sparse(SSR.bus,SSR.Iq,V.*cdt,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(SSR.bus+Bus.n,SSR.Iq,V.*sdt,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(SSR.bus,SSR.delta,-V.*cdt.*id+V.*sdt.*iq,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(SSR.bus+Bus.n,SSR.delta,V.*sdt.*id+V.*cdt.*iq,2*Bus.n,DAE.n);

end