fm_out.m

一个较好的MATLAB潮流程序

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-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 Settings Line Ltc Bus Syn Exc Oxl RLC
global DAE jay Shunt Varout

switch flag
 case 0

  chunk = Settings.chunk;
  Varout.t = zeros(chunk,1);
  if DAE.n
    Varout.x = zeros(chunk, DAE.n);
  end
  Varout.V = zeros(chunk,Bus.n);
  Varout.ang = zeros(chunk,Bus.n);
  Varout.p = zeros(chunk,Bus.n);
  Varout.q = zeros(chunk,Bus.n);
  if Syn.n
    Varout.Pm = zeros(chunk,Syn.n);
    Varout.Vf = zeros(chunk,Syn.n);
  end
  if Exc.n
    Varout.Vref = zeros(chunk,Exc.n);
  end
  if Oxl.n
    Varout.If = zeros(chunk,Oxl.n);
  end
  Varout.psflow = zeros(chunk,Line.n+Ltc.n);
  Varout.qsflow = zeros(chunk,Line.n+Ltc.n);
  Varout.prflow = zeros(chunk,Line.n+Ltc.n);
  Varout.qrflow = zeros(chunk,Line.n+Ltc.n);

 case 1

  chunk = Settings.chunk;
  Varout.t = [Varout.t; zeros(chunk,1)];
  if DAE.n
    Varout.x = [Varout.x; zeros(chunk, DAE.n)];
  end
  Varout.V = [Varout.V; zeros(chunk,Bus.n)];
  Varout.ang = [Varout.ang; zeros(chunk,Bus.n)];
  Varout.p = [Varout.p; zeros(chunk,Bus.n)];
  Varout.q = [Varout.q; zeros(chunk,Bus.n)];
  if Syn.n
    Varout.Pm = [Varout.Pm; zeros(chunk,Syn.n)];
    Varout.Vf = [Varout.Vf; zeros(chunk,Syn.n)];
  end
  if Exc.n
    Varout.Vref = [Varout.Vref; zeros(chunk,Exc.n)];
  end
  if Oxl.n
    Varout.If = [Varout.If; zeros(chunk,Oxl.n)];
  end
  Varout.psflow = [Varout.psflow; zeros(chunk,Line.n+Ltc.n)];
  Varout.qsflow = [Varout.qsflow; zeros(chunk,Line.n+Ltc.n)];
  Varout.prflow = [Varout.prflow; zeros(chunk,Line.n+Ltc.n)];
  Varout.qrflow = [Varout.qrflow; zeros(chunk,Line.n+Ltc.n)];

 case 2

 Varout.t(k) = t;
 if DAE.n
   Varout.x(k,:) = DAE.x.';
 end
 Varout.V(k,:) = DAE.V.';
 Varout.ang(k,:) = DAE.a.';

 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
 if RLC.n > 0
   P(RLC.bus1) = P(RLC.bus1) + RLC.dat(:,1);
   P(RLC.bus2) = P(RLC.bus2) + RLC.dat(:,3);
   Q(RLC.bus1) = Q(RLC.bus1) + RLC.dat(:,2);
   Q(RLC.bus2) = Q(RLC.bus2) + RLC.dat(:,4);
 end

 Varout.p(k,:) = P.';
 Varout.q(k,:) = Q.';

 if Syn.n, Varout.Pm(k,:) = Syn.pm.'; Varout.Vf(k,:) = Syn.vf.'; end
 if Exc.n, Varout.Vref(k,:) = Exc.vrif.'; end
 if Oxl.n, Varout.If(k,:) = Oxl.If.'; end

 % flussi di potenza nelle linee e nei trasformatori
 if Line.n > 0
   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));
   Varout.psflow(k,1:Line.n) = real(MW_s).';
   Varout.qsflow(k,1:Line.n) = imag(MW_s).';

   MW_r = VV(Line.to).* ...
          conj((VV(Line.to)-VV(Line.from)./tps).*y + ...
               VV(Line.to).*(jay*chrg2));
   Varout.prflow(k,1:Line.n) = real(MW_r).';
   Varout.qrflow(k,1:Line.n) = imag(MW_r).';
 end
 if Ltc.n > 0
   % active power sent out by from_bus to to_bus
   Varout.psflow(k,Line.n+1:Line.n+Ltc.n) = [-real(Ltc.dat(:,5))].';
   % reactive power sent out by from_bus to to_bus
   Varout.qsflow(k,Line.n+1:Line.n+Ltc.n) = [imag(Ltc.dat(:,5))].';
   % active power received by to_bus from from_bus
   Varout.prflow(k,Line.n+1:Line.n+Ltc.n) = [-real(Ltc.dat(:,6))].';
   % reactive power received by to_bus from from_bus
   Varout.qrflow(k,Line.n+1:Line.n+Ltc.n) = [imag(Ltc.dat(:,6))].';
 end

case 3

 Varout.t = Varout.t(1:k);
 if DAE.n
   Varout.x = Varout.x(1:k,:);
 end
 Varout.V = Varout.V(1:k,:);
 if Syn.n
   Varout.Pm = Varout.Pm(1:k,:);
   Varout.Vf = Varout.Vf(1:k,:);
 end
 if Exc.n
   Varout.Vref = Varout.Vref(1:k,:);
 end
 if Oxl.n
   Varout.If = Varout.If(1:k,:);
 end
 Varout.ang = Varout.ang(1:k,:);
 Varout.p = Varout.p(1:k,:);
 Varout.q = Varout.q(1:k,:);
 Varout.psflow = Varout.psflow(1:k,:);
 Varout.qsflow = Varout.qsflow(1:k,:);
 Varout.prflow = Varout.prflow(1:k,:);
 Varout.qrflow = Varout.qrflow(1:k,:);

end