fm_out.m

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

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);
  end
  
  if compute_pijs & Sssc.n
    V1 = DAE.V(Sssc.bus1);
    V2 = DAE.V(Sssc.bus2);
    cc = cos(DAE.a(Sssc.bus1)-DAE.a(Sssc.bus2));
    u = DAE.x(Sssc.vcs)./sqrt(V1.^2+V2.^2-2.*V1.*V2.*cc);
    ps(Sssc.line) = (1+u).*ps(Sssc.line);
    pr(Sssc.line) = (1+u).*pr(Sssc.line);
    qs(Sssc.line) = (1+u).*qs(Sssc.line);
    qr(Sssc.line) = (1+u).*qr(Sssc.line);
  end
  
  if compute_pijs & Upfc.n
    V1 = DAE.V(Upfc.bus1);
    V2 = DAE.V(Upfc.bus2);
    theta = DAE.a(Upfc.bus1)-DAE.a(Upfc.bus2)+Upfc.gamma;
    ss = sin(theta);
    cc = cos(theta);
    c1 = sqrt(DAE.x(Upfc.vp).^2+DAE.x(Upfc.vq).^2).*Upfc.y;
    ps(Upfc.line) = ps(Upfc.line) + c1.*V2.*ss;
    pr(Upfc.line) = pr(Upfc.line) - c1.*V2.*ss;
    qs(Upfc.line) = qs(Upfc.line) + c1.*V1.*cos(Upfc.gamma)-DAE.x(Upfc.iq).*V1;
    qr(Upfc.line) = qr(Upfc.line) - c1.*V2.*cc;
  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