fm_out.m

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

function  fm_out(flag,t,k)
% FM_OUT define output variable vector during time domain simulations
%
% FM_OUT(FLAG,T,K)
%   T = actual time
%   K = time step number
%   FLAG = 0 ->  Output Structure Initialization
%	   1 ->  Memory Allocation
%          2 ->  k-th Step Assignment
%          3 ->  Vector Redimensioning
%
%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 Settings Line Ltc Bus Syn Exc Oxl
global DAE jay Shunt Varout Varname Sssc Upfc Tcsc
persistent compute_pijs

switch flag

 case 0

  if DAE.n > 5000 | Bus.n > 5000
    Settings.chunk = 10;
  end
  if length(Varname.idx) > Settings.maxvar
    fm_disp('No more than 1500 variables will be stored for plotting.')
    Varname.idx = Varname.idx(1:Settings.maxvar);
  end        
  
  Varout.t = zeros(Settings.chunk,1);
  Varout.idx = Varname.idx;
  idx0 = length(Varout.idx);
  Varout.vars = zeros(Settings.chunk,idx0);
  idx1 = DAE.n+4*Bus.n+2*Syn.n+Exc.n+Oxl.n;
  if ~isempty(find(Varout.idx > idx1))
    compute_pijs = 1;
  else
    compute_pijs = 0;
  end

 case 1

  Varout.t = [Varout.t; zeros(Settings.chunk,1)];
  idx0 = length(Varout.idx);
  Varout.vars = [Varout.vars; zeros(Settings.chunk,idx0)];

 case 2

  Varout.t(k) = t;

  P = DAE.glfp;
  Q = DAE.glfq;
  if Ltc.n > 0
    P(Ltc.bus1) = P(Ltc.bus1) - real(Ltc.dat(:,5));
    P(Ltc.bus2) = P(Ltc.bus2) - real(Ltc.dat(:,6));
    Q(Ltc.bus1) = Q(Ltc.bus1) + imag(Ltc.dat(:,5));
    Q(Ltc.bus2) = Q(Ltc.bus2) + imag(Ltc.dat(:,6));
  end

  % power flows in transmission lines and transformers
  if compute_pijs & Line.n
    tps = Line.con(:,11).*exp(jay*Line.con(:,12)*pi/180);
    tps(find(tps == 0)) = 1;
    VV = DAE.V.*exp(jay*DAE.a);
    rx = Line.con(:,9);
    rr = Line.con(:,8);
    chrg = Line.con(:,10)/2;
    chrg1 = chrg - jay*Shunt.g(Line.from) + Shunt.b(Line.from);
    chrg2 = chrg - jay*Shunt.g(Line.to) + Shunt.b(Line.to);
    z = rr + jay*rx;
    y = ones(Line.n,1)./z;

    MW_s = VV(Line.from).* ...
           conj((VV(Line.from)-tps.*VV(Line.to)).*y + ...
                VV(Line.from).*(jay*chrg1))./(tps.*conj(tps));
    ps = real(MW_s).';
    qs = imag(MW_s).';

    MW_r = VV(Line.to).* ...
           conj((VV(Line.to)-VV(Line.from)./tps).*y + ...
                VV(Line.to).*(jay*chrg2));
    pr = real(MW_r).';
    qr = imag(MW_r).';
  end
  
  if compute_pijs
    [pr,ps,qr,qs] = flows(Tcsc,pr,ps,qr,qs);
    [pr,ps,qr,qs] = flows(Sssc,pr,ps,qr,qs);
    [pr,ps,qr,qs] = flows(Upfc,pr,ps,qr,qs);
  end
    
  if compute_pijs & Ltc.n > 0
    ps(1,Line.n+1:Line.n+Ltc.n) = [-real(Ltc.dat(:,5))].';
    qs(1,Line.n+1:Line.n+Ltc.n) = [imag(Ltc.dat(:,5))].';
    pr(1,Line.n+1:Line.n+Ltc.n) = [-real(Ltc.dat(:,6))].';
    qr(1,Line.n+1:Line.n+Ltc.n) = [imag(Ltc.dat(:,6))].';
  end

  idx0 = DAE.n+4*Bus.n+2*Syn.n+Exc.n+Oxl.n+4*Line.n+Ltc.n;
  if compute_pijs
    vars = [DAE.x.', DAE.V.', DAE.a.', ...
            P.', Q.', Syn.pm.', Syn.vf.', ...
            Exc.vrif.', Oxl.If.', ps, pr, qs, qr];
  else
    vars = [DAE.x.', DAE.V.', DAE.a.', ...
            P.', Q.', Syn.pm.', Syn.vf.', ...
            Exc.vrif.', Oxl.If.'];
  end
  if Settings.vs
    Varout.vars(k,[1:length(Varout.idx)]) = ...
        vars(Varout.idx(find(Varout.idx <= idx0)));
    fm_vs(k,idx0)
  else
    Varout.vars(k,:) = vars(Varout.idx);
  end

 case 3

  Varout.t = Varout.t(1:k);
  Varout.vars = Varout.vars(1:k,:);

end