fm_mixed.m

一个较好的MATLAB潮流程序

function  fm_mixed(flag)
% FM_MIXED defines Mixed Loads
%
% FM_MIXED(FLAG)
%       FLAG = 0 initialization
%       FLAG = 1 algebraic equations
%       FLAG = 2 algebraic Jacobians
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%
%Author:    Federico Milano
%Date:      23-May-2005
%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 Mixed DAE Bus Settings PQ

x = DAE.x(Mixed.x);
y = DAE.x(Mixed.y);
V1 = DAE.V(Mixed.bus);
t1 = DAE.a(Mixed.bus);
Kpf = Mixed.con(:,5);
Kpv = Mixed.con(:,6);
alpha = Mixed.con(:,7);
Tpv = Mixed.con(:,8);
Kqf = Mixed.con(:,9);
Kqv = Mixed.con(:,10);
beta = Mixed.con(:,11);
Tqv = Mixed.con(:,12);
Tfv = Mixed.con(:,13);
Tft = Mixed.con(:,14);
k = 0.5/pi/Settings.freq;
V0 = Mixed.dat(:,1);
t0 = Mixed.dat(:,2);

switch flag
 case 0 % initialization

  for k = 1:length(Mixed.n)
    idx = find(PQ.bus == Mixed.bus(k));
    if isempty(idx)
      fm_disp(['No PQ load found for initializing mixed load ', ...
               'at bus #',int2str(Mixed.bus(k))])
      Mixed.con(k,6) = 0;
      Mixed.con(k,10) = 0;
    else
      PQ.con(idx,4) = PQ.con(idx,4) - ...
          PQ.P0(idx)*sum(Mixed.con(k,[6:8]))/100;
      PQ.con(idx,5) = PQ.con(idx,5) - ...
          PQ.Q0(idx)*sum(Mixed.con(k,[9:11]))/100;
      Mixed.con(k,6) = Mixed.con(k,6)*PQ.P0(idx)/100;
      Mixed.con(k,10) = Mixed.con(k,10)*PQ.Q0(idx)/100;
      if PQ.con(idx,4) == 0 & PQ.con(idx,5) == 0,
        PQ.con(idx,:) = [];
        PQ.bus(idx) = [];
        PQ.n = PQ.n - 1;
        PQ.P0(idx) = [];
        PQ.Q0(idx) = [];
      end
    end
  end

  %check time constants
  idx = find(Tfv == 0);
  if idx
    mixedwarn(idx, 'Time constant Tfv cannot be zero. Tfv = 0.001 s will be used.'),
  end
  Mixed.con(idx,13) = 0.001;
  idx = find(Tft == 0);
  if idx
    mixedwarn(idx, 'Time constant Tft cannot be zero. Tft = 0.001 s will be used.'),
  end
  Mixed.con(idx,14) = 0.001;

  %variable initialization
  DAE.x(Mixed.x) = -V1./Tfv;
  x = DAE.x(Mixed.x);
  DAE.x(Mixed.y) = 0;
  y = DAE.x(Mixed.y);
  Mixed.dat(:,1) = V1;
  Mixed.dat(:,2) = t1;

  %check limits
  fm_disp('Initialization of mixed loads completed.')

 case 1 % algebraic equations

  DAE.gp = DAE.gp + sparse(Mixed.bus,1,Kpf.*(y+k.*(t1-t0)./Tft)+Kpv.*((V1./V0).^alpha+Tpv.*(x+V1./Tfv)),Bus.n,1);
  DAE.gq = DAE.gq + sparse(Mixed.bus,1,Kqf.*(y+k.*(t1-t0)./Tft)+Kqv.*((V1./V0).^beta+Tqv.*(x+V1./Tfv)),Bus.n,1);

 case 2 % algebraic Jacobians

  DAE.J12 = DAE.J12 + sparse(Mixed.bus,Mixed.bus,Kpv.*((V1./V0).^alpha.*alpha./V1+Tpv./Tfv),Bus.n,Bus.n);
  DAE.J11 = DAE.J11 + sparse(Mixed.bus,Mixed.bus,Kpf.*k./Tft,Bus.n,Bus.n);
  DAE.J22 = DAE.J22 + sparse(Mixed.bus,Mixed.bus,Kqv.*((V1./V0).^beta.*beta./V1+Tqv./Tfv),Bus.n,Bus.n);
  DAE.J21 = DAE.J21 + sparse(Mixed.bus,Mixed.bus,Kqf.*k./Tft,Bus.n,Bus.n);

 case 3 % differential equations

  DAE.f(Mixed.x) = (-V1./Tfv-x)./Tfv;
  DAE.f(Mixed.y) = (-k.*(t1-t0)./Tft-y)./Tft;

 case 4 % state variable Jacobians

  % DAE.Fx
  DAE.Fx = DAE.Fx + sparse(Mixed.x,Mixed.x,-1./Tfv,DAE.n,DAE.n);
  DAE.Fx = DAE.Fx + sparse(Mixed.y,Mixed.y,-1./Tft,DAE.n,DAE.n);

  % DAE.Fy
  DAE.Fy = DAE.Fy + sparse(Mixed.x,Mixed.bus+Bus.n,-1./Tfv./Tfv,DAE.n,2*Bus.n);
  DAE.Fy = DAE.Fy + sparse(Mixed.y,Mixed.bus,-k./Tft./Tft,DAE.n,2*Bus.n);

  % DAE.Gx
  DAE.Gx = DAE.Gx + sparse(Mixed.bus,Mixed.x,Kpv.*Tpv,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(Mixed.bus,Mixed.y,Kpf,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(Mixed.bus+Bus.n,Mixed.x,Kqv.*Tqv,2*Bus.n,DAE.n);
  DAE.Gx = DAE.Gx + sparse(Mixed.bus+Bus.n,Mixed.y,Kqf,2*Bus.n,DAE.n);

end


% -------------------------------------------------------------------
% function for creating warning messages
function mixedwarn(idx, msg)
fm_disp(strcat('Warning: Mixed Load #',int2str(idx),msg))