fm_simrep.m

电力系统分析计算程序

function fm_simrep(varargin)
% FM_SIMREP generate data report in Simulink models
%
% FM_SIMREP(FLAG,FONTSIZE,FONTNAME)
%       FLAG: 1 - set up voltage report on current loaded network
%             2 - wipe voltage report on current loaded network
%             3 - set up power flows on current loaded network
%             4 - wipe power flows on current loaded network
%             5 - hide component names except for buses
%             6 - show component names
%             7 - hide bus names
%             8 - show bus names
%             9 - set font name
%            10 - set font size
%            11 - save model diagram to eps file
%       FONTSIZE: font size (integer)
%       FONTNAME: font name (string)
%
%see also FM_LIB, FM_SIMSET
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://www.uclm.es/area/gsee/Web/Federico
%
% Copyright (C) 2002-2008 Federico Milano

fm_var

type = 1;

switch nargin
 case 1
  flag = varargin{1};
 case 3
  flag = varargin{1};
  fontsize = varargin{2};
  fontname = varargin{3};
 case 4
  flag = varargin{1};
  maptype = varargin{2};
  type = varargin{3};
  method = varargin{4};
end

if isempty(File.data) & type,
  fm_disp('No loaded system is present at the moment.',2),
  return
end
if isempty(strfind(File.data,'(mdl)')) & type
  fm_disp('The actual data file is not generated from a Simulink model.',2),
  return
end
if ~Settings.init & ~strcmp(flag,{'mdl2eps','ViewModel','DrawModel'}) & type
  fm_disp('Perform Power Flow before using this utility.',2),
  return
end

lasterr('');

% load Simulink model
cd(Path.data);
filedata = File.data(1:end-5);
open_sys = find_system('type','block_diagram');
if ~sum(strcmp(open_sys,filedata))
  try
    if nargin > 3
      load_system(filedata)
    else
      open_system(filedata);
    end
  catch
    fm_disp(lasterr,2)
    return
  end
end
cur_sys = get_param(filedata,'Handle');
if ~strcmp(flag,'DrawModel'), set_param(cur_sys,'Open','on'), end
blocks = find_system(cur_sys,'Type','block');
lines = find_system(cur_sys, ...
                    'FindAll','on', ...
                    'type','line');
masks = get_param(blocks,'Masktype');
nblock = length(blocks);

switch flag

 case 'ViewModel'

  % view the current Simulink model.
  hilite_system(cur_sys,'none')

 case 'DrawModel'

  maps = {'jet';'hot';'gray';'bone';'copper';'pink';
    'hsv';'cool';'autumn';'spring';'winter';'summer'};

  if ~method
    switch type
     case 0 % only one-line diagram
      zlevel = 0;
     case 1 % voltage magnitudes
      zlevel = 1.025*max(DAE.y(Bus.v));
     case 2 % voltage phases
      zmax = max(180*DAE.y(Bus.a)/pi);
      if zmax >= 0
        zlevel = 1.1*zmax;
      else
        zlevel = 0.9*zmax;
      end
     case 3 % line flows
      [sij,sji] = flows(Line,3);
      zlevel = 1.1*max(sij);
     case 4 % generator rotor angles
      if ~Syn.n
        if Settings.static
          fm_disp('Currently, synchronous machines are not loaded.')
          fm_disp('Uncheck the option "Discard dynamic data" and try again.')
        else
          fm_disp('There are no synchronous machines in the current system.',2)
        end
        return
      end
      zmax = max(180*DAE.x(Syn.delta)/pi);
      if zmax >= 0
        zlevel = 1.1*zmax;
      else
        zlevel = 0.9*zmax;
      end
     case 5 % generator rotor speeds
      if ~Syn.n
        if Settings.static
          fm_disp('Currently, synchronous machines are not loaded.')
          fm_disp('Uncheck the option "Discard dynamic data" and try again.')
        else
          fm_disp('There are no synchronous machines in the current system.',2)
        end
        return
      end
      zlevel = 1.005*max(DAE.x(Syn.omega));
     case 6 % locational marginal prices
      if ~OPF.init
        fm_disp(['Run OPF before displaying Locational Marginal ' ...
                 'Prices'],2)
        return
      end
      LMP = OPF.LMP;
      zlevel = 1.025*max(LMP);
     case 7 % nodal congestion prices
      if ~OPF.init
        fm_disp(['Run OPF before displaying Nodal Congestion Prices'], ...
                2)
        return
      end
      NCP = OPF.NCP(1:Bus.n);
      zlevel = 1.025*max(NCP);
     otherwise
      zlevel = 0;
    end
    Varout.zlevel = zlevel;
  end

  if Fig.threed & type
    figure(Fig.threed)
  elseif Fig.dir & ~type
    figure(Fig.dir)
    hdla = findobj(Fig.dir,'Tag','Axes1');
    cla(hdla)
  else
    figure
  end
  hold on

  if ~method
    if type, colorbar, end
    if Fig.threed, cla(get(Fig.threed,'UserData')), end
    lines = get_param(cur_sys,'Lines');
  end

  pos = get_param(cur_sys,'Location');
  xl = [];
  yl = [];

  for i = 1:length(lines)*(~method)
    z = zlevel*ones(length(lines(i).Points(:,1)),1);
    plot3(lines(i).Points(:,1),lines(i).Points(:,2),z,'k')
    xl = [xl; lines(i).Points([1 end],1)];
    yl = [yl; lines(i).Points([1 end],2)];
  end

  xb = [];
  yb = [];
  zb = [];
  x_max = 0;
  x_min = 2000;
  y_max = 0;
  y_min = 2000;
  idx = 0;
  idx_line = 0;
  idx_gen = 0;
  Compnames = cell(0,0);
  Compidx = [];
  ncomp = 0;
  if ~method
    switch type
     case 3
      Varout.hdl = zeros(Line.n,1);
     case {4,5}
      Varout.hdl = zeros(Syn.n,1);
     otherwise
      Varout.hdl = zeros(Bus.n,1);
    end
  end

  for i = 1:length(blocks)*(~method)

    bmask = get_param(blocks(i),'MaskType');

    jdx = strmatch(bmask,Compnames,'exact');
    if isempty(jdx)
      ncomp = ncomp + 1;
      Compnames{ncomp,1} = bmask;
      Compidx(ncomp,1) = 0;
      jdx = ncomp;
    end
    Compidx(jdx) = Compidx(jdx)+1;

    bpos = get_param(blocks(i),'Position');
    bidx = Compidx(jdx); %str2double(get_param(blocks(i),'UserData'));
    borien = get_param(blocks(i),'Orientation');
    bports = get_param(blocks(i),'Ports');
    bvalues = get_param(blocks(i),'MaskValues');
    bin = sum(bports([1 6]));
    bou = sum(bports([2 7]));

    switch bmask
     case 'Varname'
      % nothing to do!
      x = cell(1,1);
      y = cell(1,1);
      s = cell(1,1);
      x{1} = 0;
      y{1} = 0;
      s{1} = 'k';
     case 'Ltc'
      bname = get_param(blocks(i),'NamePlacement');
      [x,y,s] = mask(eval(bmask),bidx,{borien;bname},bvalues);
     case 'Line'
      bdescr = get_param(blocks(i),'MaskDescription');
      [x,y,s] = mask(eval(bmask),bidx,borien,bdescr);
     case 'PQ'
      [x,y,s] = mask(eval(bmask),bidx,borien,'PQ');
     case 'PQgen'
      [x,y,s] = mask(eval(bmask),bidx,borien,'PQgen');
     case 'Link'

      rot = strcmp(get_param(blocks(i),'NamePlacement'),'normal');

      x = cell(6,1);
      y = cell(6,1);
      s = cell(6,1);

      x{1} = [0.45 0];
      y{1} = [0.5 0.5];
      s{1} = 'k';

      x{2} = [1 0.55];
      y{2} = [0.5 0.5];
      s{2} = 'k';

      x{3} = [0.45 0.55 0.55 0.45 0.45];
      y{3} = [0.45 0.45 0.55 0.55 0.45];
      s{3} = 'g';

      x{4} = [0.5 0.5];
      y{4} = rot*0.45+[0.1 0.45];
      s{4} = 'g';

      x{5} = [0.45 0.55 0.55 0.45 0.45];
      y{5} = rot*0.9+[0 0 0.1 0.1 0];
      s{5} = 'g';

      x{6} = 1-rot;
      y{6} = 1-rot;
      s{6} = 'w';

      [x,y] = fm_maskrotate(x,y,borien);

     case 'Link2'

      x = cell(10,1);
      y = cell(10,1);
      s = cell(10,1);

      x{1} = [0 0.45];
      y{1} = [0.5 0.5];
      s{1} = 'k';

      x{2} = [0.5 0.5];
      y{2} = [0 0.4];
      s{2} = 'k';

      x{3} = [0.5 0.5];
      y{3} = [0 0.4];
      s{3} = 'k';

      x{4} = [0.5 1];
      y{4} = [0 0];
      s{4} = 'k';

      x{5} = [0.5 0.5];
      y{5} = [0.6 1];
      s{5} = 'g';

      x{6} = [0.5 0.9];
      y{6} = [1 1];
      s{6} = 'g';

      x{7} = [0.9 0.985 0.985 0.9 0.9];
      y{7} = [0.1 0.1 -0.1 -0.1 0.1]+1;
      s{7} = 'g';

      x{8} = [0.45 0.55 0.55 0.45 0.45];
      y{8} = [0.6 0.6 0.4 0.4 0.6];
      s{8} = 'g';

      x{9} = 0;
      y{9} = -0.5;
      s{9} = 'w';

      x{10} = 0;
      y{10} = 1.5;
      s{10} = 'w';

      [x,y] = fm_maskrotate(x,y,borien);

     otherwise
      [x,y,s] = mask(eval(bmask),bidx,borien,bvalues);
    end

    xt = [];
    yt = [];
    for j = 1:length(x)
      xt = [xt, x{j}];
      yt = [yt, y{j}];
    end

    xmin = min(xt);
    xmax = max(xt);
    if xmax == xmin, xmax = xmin+1; end
    ymin = min(yt);
    ymax = max(yt);
    if ymax == ymin, ymax = ymin+1; end
    dx = bpos(3)-bpos(1);
    dy = bpos(4)-bpos(2);

    xscale = dx/(xmax-xmin);
    yscale = dy/(ymax-ymin);

    xcorr = -xscale*(xmax+xmin)/2 + 0.5*dx;
    ycorr = -yscale*(ymax+ymin)/2 + 0.5*dy;

    xmean = bpos(1)+xcorr+xscale*(xmax+xmin)/2;
    ymean = bpos(4)-ycorr-yscale*(ymax+ymin)/2;

    if strcmp(bmask,'Bus')
      idx = idx + 1;

      if type == 1 | type == 2 | type == 6 | type == 7
        xb = [xb; bpos(1)+xcorr+xscale*xmax; xmean];
        yb = [yb; bpos(4)-yscale*ymax-ycorr; ymean];
        xb = [xb; bpos(1)+xcorr+xscale*xmin];
        yb = [yb; bpos(4)-yscale*ymin-ycorr];
      end

      %bcolor = get_param(blocks(i),'BackgroundColor');
      bname = get_param(blocks(i),'Name');

      switch borien
       case {'right','left'}
        switch get_param(blocks(i),'NamePlacement')
         case 'normal'
          xtext = xmean;
          ytext = bpos(4)-ycorr-yscale*ymin+7;
          bha = 'center';
          bva = 'top';
         case 'alternate'
          xtext = xmean;
          ytext = bpos(4)-ycorr-yscale*ymax-7;
          bha = 'center';
          bva = 'bottom';
        end
       case {'up','down'}
        switch get_param(blocks(i),'NamePlacement')
         case 'normal'
          xtext = bpos(1)+xcorr+xscale*xmax+7;
          ytext = ymean;
          bha = 'left';
          bva = 'middle';
         case 'alternate'
          xtext = bpos(1)+xcorr+xscale*xmin-7;
          ytext = ymean;
          bha = 'right';
          bva = 'middle';
        end
      end

      % write bus names
      if type
        h = text(xtext,ytext,zlevel,bname);
        set(h,'HorizontalAlignment',bha,'VerticalAlignment',bva, ...
              'FontSize',8)
      end
      if type == 1 | type == 2 | type == 6 | type == 7
        switch type
         case 1, zpeak = DAE.y(Bus.v(idx));
         case 2, zpeak = 180*DAE.y(Bus.a(idx))/pi;
         case 6, zpeak = LMP(idx);
         case 7, zpeak = NCP(idx);
        end
        Varout.hdl(idx) = plot3([xmean xmean],[ymean ymean], ...
                              [zlevel,zpeak],'k:');
      end
    end
    if strcmp(bmask,'Line') & type == 3
      idx_line = idx_line + 1;

      xb = [xb; bpos(1)+xcorr+xscale*xmax; xmean];
      yb = [yb; bpos(4)-yscale*ymax-ycorr; ymean];
      xb = [xb; bpos(1)+xcorr+xscale*xmin];
      yb = [yb; bpos(4)-yscale*ymin-ycorr];

      Varout.hdl(idx_line) = plot3([xmean xmean],[ymean ymean],[zlevel,sij(idx_line)],'k:');

    elseif strcmp(bmask,'Syn') & (type == 4 | type == 5)

      idx_gen = idx_gen + 1;

      xb = [xb; bpos(1)+xcorr+xscale*xmax; xmean];
      yb = [yb; bpos(4)-yscale*ymax-ycorr; ymean];
      xb = [xb; bpos(1)+xcorr+xscale*xmin];
      yb = [yb; bpos(4)-yscale*ymin-ycorr];

      switch type
       case 4
        zpeak = 180*DAE.x(Syn.delta(idx_gen))/pi;
       case 5