fm_sofc.m

这是一个很适合研究和学习用的电力系统仿真软件

function  fm_sofc(flag)
% FM_SOFC define Solyd Oxyde Fuel Cells
%
% FM_SOFC(FLAG)
%       FLAG = 0 -> initialization
%       FLAG = 1 -> algebraic equations
%       FLAG = 2 -> algebraic Jacobians
%       FLAG = 3 -> differential equations
%       FLAG = 4 -> state Jacobians
%       FLAG = 5 -> non-windup limiters
%
%Author:    Federico Milano
%Date:      25-Nov-2003
%Version:   2.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global Sofc DAE Bus PV SW

Ik = DAE.x(Sofc.Ik);
Vk = DAE.x(Sofc.Vk);
Vb = DAE.V(Sofc.bus);
pH2 = DAE.x(Sofc.pH2);
pH20 = DAE.x(Sofc.pH20);
pO2 = DAE.x(Sofc.pO2);
qH2 = DAE.x(Sofc.qH2);
m = DAE.x(Sofc.m);
Sn = Sofc.con(:,2);
Vn = Sofc.con(:,3);
Te = Sofc.con(:,4);
TH2 = Sofc.con(:,5);
KH2 = Sofc.con(:,6);
Kr = Sofc.con(:,7);
TH20 = Sofc.con(:,8);
KH20 = Sofc.con(:,9);
TO2 = Sofc.con(:,10);
KO2 = Sofc.con(:,11);
rHO = Sofc.con(:,12);
Tf = Sofc.con(:,13);
Uopt = Sofc.con(:,14);
r = Sofc.con(:,17);
N0 = Sofc.con(:,18);
E0 = Sofc.con(:,19);
RTon2F = 0.5*8.314*Sofc.con(:,20)/96487;
Pref = Sofc.con(:,21);
Vref = Sofc.con(:,22);
xt = Sofc.con(:,26);
Km = Sofc.con(:,27);
Tm = Sofc.con(:,28);
m_max = Sofc.con(:,29);
m_min = Sofc.con(:,30);

switch flag
 case 0 % initialization

  % Check time constants
  idx = find(Te == 0);
  if idx
    fcwarn(idx, ['Time constant Te cannot be zero. ', ...
		 'Te = 0.001 s will be used.'])
    Sofc.con(idx,5) = 0.001;
  end

  idx = find(TH2 == 0);
  if idx
    fcwarn(idx, ['Time constant TH2 cannot be zero. ', ...
	   'TH2 = 0.001 s will be used.'])
    Sofc.con(idx,7) = 0.001;
  end

  idx = find(TH20 == 0);
  if idx
    fcwarn(idx, ['Time constant TH20 cannot be zero. ', ...
	   'TH20 = 0.001 s will be used.'])
    Sofc.con(idx,10) = 0.001;
  end

  idx = find(TO2 == 0);
  if idx
    fcwarn(idx, ['Time constant TO2 cannot be zero. ', ...
	   'TO2 = 0.001 s will be used.'])
    Sofc.con(idx,12) = 0.001;
 end

  idx = find(Tf == 0);
  if idx
    fcwarn(idx, ['Time constant Tf cannot be zero. ', ...
	   'Tf = 0.001 s will be used.'])
    Sofc.con(idx,15) = 0.001;
  end

  idx = find(Tm == 0);
  if idx
    fcwarn(idx, ['Time constant Tm cannot be zero. ', ...
		 'Tm = 0.001 s will be used.'])
    Sofc.con(idx,28) = 0.001;
  end

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

  % Get generator powers
  Pfc = Bus.Pg(Sofc.bus);
  Qfc = Bus.Qg(Sofc.bus);

  % Define reference active power
  Sofc.con(:,21) = Pfc;

  % Fuel Cell stae variable initialization
  DAE.x(Sofc.Ik) = Sn.*Pfc./Vb./Vn;
  Ik = DAE.x(Sofc.Ik);
  DAE.x(Sofc.qH2) = 2*Kr.*Ik./Uopt;
  qH2 = DAE.x(Sofc.qH2);
  DAE.x(Sofc.pO2) = (qH2./rHO-Kr.*Ik)./KO2;
  pO2 = DAE.x(Sofc.pO2);
  DAE.x(Sofc.pH2) = (qH2 - 2*Kr.*Ik)./KH2;
  pH2 = DAE.x(Sofc.pH2);
  DAE.x(Sofc.pH20) = 2*Kr.*Ik./KH20;
  pH20 = DAE.x(Sofc.pH20);
  DAE.x(Sofc.Vk) = -r.*Ik./Vn + ...
      N0.*(E0+RTon2F.*log(pH2.*sqrt(pO2)./pH20))./Vn;
  Vk = DAE.x(Sofc.Vk);

  % Base power and voltage
  Sofc.con(:,23) = (Ik.*Vk.*Vn./Sn)./Pfc;
  Sofc.con(:,24) = Vb./Vk;

  % Adjust value of control type
  Sofc.con(:,25) = ~Sofc.con(:,25);
  idx = find(Sofc.con(:,25));
  if idx, Sofc.con(idx,25) = Vk(idx); end

  % Initialize tap ratio
  DAE.x(Sofc.m) = sqrt(((xt./DAE.V(Sofc.bus)./Vk./Sofc.con(:,24)).^2) ...
		       .*(Pfc.^2+(Qfc+(DAE.V(Sofc.bus).^2)./xt).^2));

  % Define reference AC voltage
  Sofc.con(:,22) = DAE.V(Sofc.bus)+DAE.x(Sofc.m)./Km;

  if ~check
    fm_disp('Fuel cells cannot be properly initialized.')
  else
    fm_disp('Initialization of Solid Oxyde Fuel Cells completed.')
  end

 case 1 % algebraic equations

  DAE.gp = DAE.gp - sparse(Sofc.bus,1,Ik.*Vk.*Vn./Sn.*Sofc.con(:,24),Bus.n,1);
  Vt = m.*Vk.*Sofc.con(:,24);
  Vs = DAE.V(Sofc.bus);
  Q = -Vs.*(Vs - Vt.*sqrt(1-(xt.*Ik./Vs.*Vn./Sn./m).^2))./xt;
  DAE.gq = DAE.gq - sparse(Sofc.bus,1,Q,Bus.n,1);

 case 2

  Vt = m.*Vk.*Sofc.con(:,24);
  Vs = DAE.V(Sofc.bus);
  sq = sqrt(1-(xt.*Ik./Vs.*Vn./Sn./m).^2);
  dQdv = -2*Vs./xt+Vt./xt.*sq+0.5*Vs.*Vt./xt.*(2*((xt.*Ik.*Vn./Sn./m).^2)./(Vs.^3))./sq;
  DAE.J22 = DAE.J22 - sparse(Sofc.bus,Sofc.bus,dQdv,Bus.n,Bus.n);

 case 3 % differential equations

  Pref = Sofc.con(:,21).*Sofc.con(:,23);
  Umax = 0.5*Sofc.con(:,15).*qH2./Kr;
  Umin = 0.5*Sofc.con(:,16).*qH2./Kr;
  Input = Sn.*Pref./Vk./Vn;
  idx = find(Sofc.con(:,25));
  if idx, Input(idx) = Sn(idx).*Pref(idx)./Sofc.con(idx,25)./Vn(idx); end

  idx = find(Input > Umax);
  if idx, DAE.f(Sofc.Ik(idx)) = (Umax(idx) - Ik(idx))./Te(idx); end
  idx = find(Input < Umin);
  if idx, DAE.f(Sofc.Ik(idx)) = (Umin(idx) - Ik(idx))./Te(idx); end
  idx = find(Input <= Umax & Input >= Umin);
  if idx, DAE.f(Sofc.Ik(idx)) = (Input(idx) - Ik(idx))./Te(idx); end

  DAE.f(Sofc.Vk) = (-Vk-r.*Ik./Vn+N0.*(E0+RTon2F.*log(pH2.*sqrt(pO2)./pH20))./Vn)./1e-3;
  DAE.f(Sofc.pH2) = ((qH2-2.*Kr.*Ik)./KH2-pH2)./TH2;
  DAE.f(Sofc.pH20) = (2.*Kr.*Ik./KH20-pH20)./TH20;
  DAE.f(Sofc.pO2) = ((qH2./rHO-Kr.*Ik)./KO2-pO2)./TO2;
  DAE.f(Sofc.qH2) = (2.*Kr.*Ik./Uopt-qH2)./Tf;

  DAE.f(Sofc.m) = (Km.*(Vref-DAE.V(Sofc.bus))-m)./Tm;
  idx = find(m >= m_max & DAE.f(Sofc.m) > 0);
  if idx, DAE.f(Sofc.m(idx)) = 0; end
  DAE.x(Sofc.m) = min(DAE.x(Sofc.m),m_max);
  idx = find(m <= m_min & DAE.f(Sofc.m) < 0);
  if idx, DAE.f(Sofc.m(idx)) = 0; end
  DAE.x(Sofc.m) = max(DAE.x(Sofc.m),m_min);

 case 4 % state variable Jacobians

  Pref = Sofc.con(:,21).*Sofc.con(:,23);

  % DAE.Fx
  Umax = 0.5*Sofc.con(:,15).*qH2./Kr;
  Umin = 0.5*Sofc.con(:,16).*qH2./Kr;
  Input = Sn.*Pref./Vk./Vn;
  idxc = find(Sofc.con(:,25));
  if idxc
    Input(idxc) = Sn(idxc).*Pref(idxc)./Sofc.con(idxc,25)./Vn(idxc);
  end

  idx = find(Input > Umax);
  if idx
    DAE.Fx = DAE.Fx + sparse(Sofc.Ik(idx),Sofc.qH2(idx), ...
			     0.5*Sofc.con(idx,16)./Kr(idx),DAE.n,DAE.n);
  end
  idx = find(Input < Umin);
  if idx
    DAE.Fx = DAE.Fx + sparse(Sofc.Ik(idx),Sofc.qH2(idx), ...
			     0.5*Sofc.con(idx,16)./Kr(idx),DAE.n,DAE.n);
  end
  idx = find(Input <= Umax & Input >= Umin);
  idxc = find(~Sofc.con(idx,25));
  if idxc, idx = idx(idxc); end
  if idx
    DAE.Fx = DAE.Fx + sparse(Sofc.Ik(idx),Sofc.Vk(idx), ...
			     -Sn(idx).*Pref(idx)./Vk(idx).^2 ...
			     ./Vn(idx)./Te(idx),DAE.n,DAE.n);
  end

  DAE.Fx = DAE.Fx + sparse(Sofc.Ik,Sofc.Ik,-1./Te,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.Vk,Sofc.Vk,-1e3,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.Vk,Sofc.Ik,-r./Vn./1e-3,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.Vk,Sofc.pH2,N0.*RTon2F./pH2./Vn./1e-3,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.Vk,Sofc.pH20,-N0.*RTon2F./pH20./Vn./1e-3,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.Vk,Sofc.pO2,0.5*N0.*RTon2F./pO2./Vn./1e-3,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.pH2,Sofc.Ik,(-2.*Kr./KH2)./TH2,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.pH2,Sofc.pH2,(-1)./TH2,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.pH2,Sofc.qH2,(1./KH2)./TH2,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.pH20,Sofc.Ik,(2.*Kr./KH20)./TH20,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.pH20,Sofc.pH20,(-1)./TH20,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.pO2,Sofc.Ik,(-Kr./KO2)./TO2,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.pO2,Sofc.pO2,(-1)./TO2,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.pO2,Sofc.qH2,(1./rHO./KO2)./TO2,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.qH2,Sofc.Ik,(2.*Kr./Uopt)./Tf,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Sofc.qH2,Sofc.qH2,(-1)./Tf,DAE.n,DAE.n);

  Vt = m.*Vk.*Sofc.con(:,24);
  Vs = DAE.V(Sofc.bus);
  sq = sqrt(1-(xt.*Ik./Vs.*Vn./Sn./m).^2);
  dQdi = 0.5*Vs.*Vt./xt.*(-2*((xt.*Vn./Sn./m./Vs).^2).*Ik)./sq;
  dQdv = m.*Sofc.con(:,24).*Vs./xt.*sq;

  DAE.Gx = DAE.Gx - sparse(Sofc.bus+Bus.n,Sofc.Vk,dQdv,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx - sparse(Sofc.bus+Bus.n,Sofc.Ik,dQdi,2*Bus.n,DAE.n);

  idx = find(m <= m_max & m >= m_min & DAE.f(Sofc.m) ~= 0);
  if ~isempty(idx)
    DAE.Fx = DAE.Fx + sparse(Sofc.m,Sofc.m,-1./Tm,DAE.n,DAE.n);
    DAE.Fy = DAE.Fy + sparse(Sofc.m,Sofc.bus+Bus.n,-Km./Tm,DAE.n,2*Bus.n);
    dQdm = Vk.*Sofc.con(:,24).*Vs./xt.*sq + ...
	   0.5*Vs.*Vt./xt.*(2*((xt.*Ik.*Vn./Sn./Vs).^2)./(m.^3))./sq;
    DAE.Gx = DAE.Gx - sparse(Sofc.bus+Bus.n,Sofc.m,dQdm,2*Bus.n,DAE.n);
  end

  % DAE.Gx
  DAE.Gx = DAE.Gx - sparse(Sofc.bus,Sofc.Ik,Vk.*Vn./Sn,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx - sparse(Sofc.bus,Sofc.Vk,Ik.*Vn./Sn.*Sofc.con(:,24),2*Bus.n,DAE.n);

 case 5 % non-windup limiters

  idx = find((m >= m_max | m <= m_min) & DAE.f(Sofc.m) == 0);
  if ~isempty(idx)
    global Settings
    k = Sofc.m(idx);
    DAE.tn(k) = 0;
    DAE.Ac(:,k) = 0;
    DAE.Ac(k,:) = 0;
    DAE.Ac(k,k) = -speye(length(idx));
  end

end

% -------------------------------------------------------------------
% function for creating warning messages

function fcwarn(idx, msg)
fm_disp(strcat('Warning: FC #',int2str(idx),msg))