fm_wind.m

一个较好的MATLAB潮流程序

function fm_wind(flag)
% FM_WIND define wind speed measures
%
% FM_WIND(FLAG)
%       FLAG = 0 initializations
%       FLAG = 3 differential equations
%       FLAG = 4 state matrix
%
%see also FM_CSWT
%
%Author:    Federico Milano
%Date:      01-Dic-2003
%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 Wind Bus DAE Settings

switch flag
 case -1

  if ~Settings.init
    fm_choice('No data found. Solve Power Flow first.',2)
    return
  end
  if ~Wind.n
    fm_choice('No Wind data found!',2)
    return
  end

  % plot Wind speeds
  colors = {'b','g','r','c','m','y','k'};
  figure
  hold on
  for i = 1:Wind.n
    leg{i} = ['v_{w',num2str(i),'}'];
    plot(Wind.speed(i).time,Wind.speed(i).vw*Wind.con(i,2),colors{rem(i-1,7)+1})
  end
  hold off
  legend(leg)
  title('Wind Speeds')
  xlabel('time [s]')
  ylabel('v_w [m/s]')
  box('on')

 case 0

  % initialize average wind speed;
  Wind.vwa = DAE.x(Wind.vw);

  for i = 1:Wind.n
    t0 = Settings.t0;
    tf = Settings.tf;
    dt = Wind.con(i,5);
    Wind.speed(i).time = [t0:dt:tf]';
    type = Wind.con(i,1);

    switch type

     case 1 % resample data in case of measurement data

      if ~isempty(Wind.speed(i).vw)
        told = Wind.speed(i).vw(:,1);
        vold = Wind.speed(i).vw(:,2);
        Wind.speed(i).vw = interp1(told,vold,Wind.speed(i).time)/Wind.con(i,2);
        Wind.speed(i).vw(1) = Wind.vwa(i);
      else % if no data is found, Weibull distribution is used
        n = length(Wind.speed(i).time);
        c = Wind.con(i,6);
        if c <= 0, c = 5; end
        k = Wind.con(i,7);
        if k <= 0, k = 2; end
        Wind.speed(i).vw = (-log(rand(n,1))/c).^(1/k);
        % set average wind speed as initial wind speed
        mean_vw = mean(Wind.speed(i).vw);
        Wind.speed(i).vw(1) = Wind.vwa(i);
        Wind.speed(i).vw(2:end) = abs(1+Wind.speed(i).vw(2:end)-mean_vw)*Wind.vwa(i);
       end

     case 2 % Weibull random distribution

      n = length(Wind.speed(i).time);
      c = Wind.con(i,6);
      if c <= 0, c = 5; end
      k = Wind.con(i,7);
      if k <= 0, k = 2; end
      Wind.speed(i).vw = (-log(rand(n,1))/c).^(1/k);
      % set average wind speed as initial wind speed
      mean_vw = mean(Wind.speed(i).vw);
      Wind.speed(i).vw(1) = Wind.vwa(i);
      Wind.speed(i).vw(2:end) = abs(1+Wind.speed(i).vw(2:end)-mean_vw)*Wind.vwa(i);

     case 3 % Composite wind model

      n = length(Wind.speed(i).time);
      % average wind speed
      vwa = Wind.vwa(i);

      % wind ramp
      Tsr = Wind.con(i,8);
      Ter = Wind.con(i,9);
      Awr = Wind.con(i,10);
      if Tsr > Ter
        fm_disp(['Start ramp time Tsr cannot be greater than end ramp ' ...
                 'time Ter'],2)
        fm_disp('Ter = Tsr + 10s will be used.')
        Ter = Tsr + 10;
      end
      vwr = zeros(n,1);
      idxmax = find(Wind.speed(i).time > Ter);
      if ~isempty(idxmax)
        vwr(idxmax) = Awr;
      end
      idxramp = find(Wind.speed(i).time <= Ter & ...
                     Wind.speed(i).time >= Tsr);
      if ~isempty(idxramp)
        vwr(idxramp) = Awr*(Wind.speed(i).time(idxramp)-Tsr)/(Ter-Tsr);
      end

      % wind gust
      Tsg = Wind.con(i,11);
      Teg = Wind.con(i,12);
      Awg = Wind.con(i,13);
      if Tsg > Teg
        fm_disp(['Start gust time Tsg cannot be greater than end gust ' ...
                 'time Teg'],2)
        fm_disp('Teg = Tsg + 10s will be used.')
        Teg = Tsg + 10;
      end
      vwg = zeros(n,1);
      idxmax = find(Wind.speed(i).time > Teg);
      if ~isempty(idxmax)
        vwg(idxmax) = Awg;
      end
      idxgust = find(Wind.speed(i).time <= Teg & ...
                     Wind.speed(i).time >= Tsg);
      if ~isempty(idxgust)
        vwg(idxgust) = 0.5*Awg*(1-cos(2*pi*(Wind.speed(i).time(idxgust)-Tsr)/(Ter-Tsr)));
      end

      % wind turbolence
      h  = Wind.con(i,14);
      z0 = Wind.con(i,15);
      df = Wind.con(i,16);
      nh = round(Wind.con(i,17));
      if h < 30
        el = 20*h;
      else
        el = 600;
      end
      f = 0;
      vwt = zeros(n,1);
      for hh = 1:nh
        f = f + df;
        phi = 2*pi*rand;
        Swt = (el*vwa/(log(h/z0))^2)/(1+1.5*f*el/vwa)^(5/3);
        vwt = vwt + sqrt(Swt*df)*cos(2*pi*Wind.speed(i).time*f+phi+0.05*pi*rand(n,1));
      end

      % composite wind
      Wind.speed(i).vw = vwa + vwr + vwg + vwt;
      Wind.speed(i).vw(1) = vwa;

    end

  end

 case 3

  % define wind speed
  Vw = zeros(Wind.n,1);
  t = DAE.t;
  if t < 0
    t = Settings.t0;
  end
  for i = 1:Wind.n
    Vw(i) = interp1(Wind.speed(i).time,Wind.speed(i).vw,t);
  end
  % filter wind speed
  DAE.f(Wind.vw) = (Vw-DAE.x(Wind.vw))./Wind.con(:,4);

 case 4

  DAE.Fx = DAE.Fx + sparse(Wind.vw,Wind.vw,-1./Wind.con(:,4),DAE.n,DAE.n);

end