fm_simrep.m

电力系统分析计算程序

        zpeak = DAE.x(Syn.omega(idx_gen));
      end
      Varout.hdl(idx_gen) = plot3([xmean xmean],[ymean ymean], ...
                              [zlevel,zpeak],'k:');
    end

    switch borien
     case 'right'
      len = yscale*(ymax-ymin);
      if bin, in_off = len/bin; end
      if bou, ou_off = len/bou; end
      for j = 1:bin
        yi = bpos(4)-ycorr-yscale*ymin-in_off/2-in_off*(j-1);
        xi = bpos(1)+xcorr+xscale*xmin;
        [yf,xf] = closerline(yi,xi-5,yl,xl);
        plot3([xi, xf],[yf, yf],[zlevel, zlevel],'k')
      end
      for j = 1:bou
        yi = bpos(4)-ycorr-yscale*ymin-ou_off/2-ou_off*(j-1);
        xi = bpos(1)+xcorr+xscale*xmax;
        [yf,xf] = closerline(yi,xi+5,yl,xl);
        plot3([xi, xf],[yf, yf],[zlevel, zlevel],'k')
      end
     case 'left'
      len = yscale*(ymax-ymin);
      if bin, in_off = len/bin; end
      if bou, ou_off = len/bou; end
      for j = 1:bin
        yi = bpos(4)-ycorr-yscale*ymin-in_off/2-in_off*(j-1);
        xi = bpos(1)+xcorr+xscale*xmax;
        [yf,xf] = closerline(yi,xi+5,yl,xl);
        plot3([xi, xf],[yf, yf],[zlevel, zlevel],'k')
      end
      for j = 1:bou
        yi = bpos(4)-ycorr-yscale*ymin-ou_off/2-ou_off*(j-1);
        xi = bpos(1)+xcorr+xscale*xmin;
        [yf,xf] = closerline(yi,xi-5,yl,xl);
        plot3([xi, xf],[yf, yf],[zlevel, zlevel],'k')
      end
     case 'up'
      len = xscale*(xmax-xmin);
      if bin, in_off = len/bin; end
      if bou, ou_off = len/bou; end
      for j = 1:bin
        yi = bpos(4)-ycorr-yscale*ymin;
        xi = bpos(1)+xcorr+xscale*xmin+in_off/2+in_off*(j-1);
        [xf,yf] = closerline(xi,yi+5,xl,yl);
        plot3([xf, xf],[yi, yf],[zlevel, zlevel],'k')
      end
      for j = 1:bou
        yi = bpos(4)-ycorr-yscale*ymax;
        xi = bpos(1)+xcorr+xscale*xmin+ou_off/2+ou_off*(j-1);
        [xf,yf] = closerline(xi,yi-5,xl,yl);
        plot3([xf, xf],[yi, yf],[zlevel, zlevel],'k')
      end
     case 'down'
      len = xscale*(xmax-xmin);
      if bin, in_off = len/bin; end
      if bou, ou_off = len/bou; end
      for j = 1:bin
        yi = bpos(4)-ycorr-yscale*ymax;
        xi = bpos(1)+xcorr+xscale*xmin+in_off/2+in_off*(j-1);
        [xf,yf] = closerline(xi,yi-5,xl,yl);
        plot3([xf, xf],[yi, yf],[zlevel, zlevel],'k')
      end
      for j = 1:bou
        yi = bpos(4)-ycorr-yscale*ymin;
        xi = bpos(1)+xcorr+xscale*xmin+ou_off/2+ou_off*(j-1);
        [xf,yf] = closerline(xi,yi+5,xl,yl);
        plot3([xf, xf],[yi, yf],[zlevel, zlevel],'k')
      end
    end

    for j = 1:length(x)
      z = zlevel*ones(length(x{j}),1);
      xx = bpos(1)+xcorr+xscale*(x{j});
      yy = bpos(4)-yscale*(y{j})-ycorr;
      if ~type
        if ~isempty(xx)
          x_max = max(x_max,max(xx));
          x_min = min(x_min,min(xx));
        end
        if ~isempty(yy)
          y_max = max(y_max,max(yy));
          y_min = min(y_min,min(yy));
        end
      end
      plot3(xx,yy,z,s{j})
    end
    x1 = [xlim]';
    x2 = [ylim]';
    x1mean = 0.5*(x1(1)+x1(2));
    x2mean = 0.5*(x2(1)+x2(2));
    xba = [xb;x1(1);x1(1);x1(2);x1(2);x1mean;x1mean;x1(1);x1(2)];
    yba = [yb;x2(1);x2(2);x2(1);x2(2);x2(1);x2(2);x2mean;x2mean];
    Varout.xb = xba;
    Varout.yb = yba;
  end

  if ~type % draw only one-line diagram
    xframe = 0.05*(x_max-x_min);
    yframe = 0.05*(y_max-y_min);
    set(hdla,'XLim',[x_min-xframe x_max+xframe],'YLim',[y_min-yframe y_max+yframe])
    hold off
    return
  end

  x1 = [xlim]';
  x2 = [ylim]';

  switch type
   case 1
    zb = formz(DAE.y(Bus.v),Bus.n);
    if abs(mean(DAE.y(Bus.v))-1) < 1e-3
      zba = [zb; 0.9999*ones(8,1)];
    else
      zba = [zb; ones(8,1)];
    end
   case 2
    zb = formz(180*DAE.y(Bus.a)/pi,Bus.n);
    zba = [zb; mean(180*DAE.y(Bus.a)/pi)*ones(8,1)];
   case 3
    [sij,sji] = flows(Line,3);
    zb = formz(sij,Line.n);
    zba = [zb; mean(sij)*ones(8,1)];
   case 4
    zb = formz(180*DAE.x(Syn.delta)/pi,Syn.n);
    zba = [zb; mean(180*DAE.x(Syn.delta)/pi)*ones(8,1)];
   case 5
    zb = formz(DAE.x(Syn.omega),Syn.n);
    zba = [zb; 0.999*ones(8,1)];
   case 6
    zb = formz(LMP,Bus.n);
    zba = [zb; mean(LMP)*ones(8,1)];
   case 7
    zb = formz(NCP,Bus.n);
    zba = [zb; mean(NCP)*ones(8,1)];
  end
  [XX,YY] = meshgrid(x1(1):5:x1(2),x2(1):5:x2(2));
  ZZ = griddata(Varout.xb,Varout.yb,zba,XX,YY,'cubic',{''});
  if method
    zlevel = Varout.zlevel;
    switch type
     case 1
      for i = 1:Bus.n
        set(Varout.hdl(i),'ZData',[zlevel,DAE.y(Bus.v(i))]);
      end
     case 2
      for i = 1:Bus.n
        set(Varout.hdl(i),'ZData',[zlevel,180*DAE.y(Bus.a(i))/pi]);
      end
     case 3
      [sij,sji] = flows(Line,3);
      for i = 1:Line.n
        set(Varout.hdl(i),'ZData',[zlevel,sij(i)]);
      end
     case 4
      for i = 1:Syn.n
        set(Varout.hdl(i),'ZData',[zlevel,180*DAE.x(Syn.delta(i))/pi]);
      end
     case 5
      for i = 1:Syn.n
        set(Varout.hdl(i),'ZData',[zlevel,zlevel,DAE.x(Syn.omega(i))]);
      end
     case 6
      for i = 1:Bus.n
        set(Varout.hdl(i),'ZData',[zlevel,zlevel,LMP(i)]);
      end
     case 7
      for i = 1:Bus.n
        set(Varout.hdl(i),'ZData',[zlevel,zlevel,NCP(i)]);
      end
    end
    delete(Varout.surf)
    if strcmp(Settings.xlabel,'Loading Parameter \lambda (p.u.)')
      xlabel([Settings.xlabel,' = ',sprintf('%8.4f',DAE.lambda)])
    else
      xlabel([Settings.xlabel,' = ',sprintf('%8.4f',DAE.t)])
    end
    Varout.surf = surf(XX,YY,ZZ);
    alpha(Varout.alpha)
    shading interp
    axis manual
    Varout.movie(end+1) = getframe(findobj(Fig.threed,'Tag','Axes1'));
    %Varout.movie(end+1) = getframe(Fig.threed);
  else
    if type == 1 & Varout.caxis
      caxis([0.9 1.1])
    else
      caxis('auto')
    end
    Varout.surf = surf(XX,YY,ZZ);
    axis auto
    shading interp
    alpha(Varout.alpha)
    xlabel('')
    set(gca,'YDir','reverse')
    set(gca,'XTickLabel',[])
    set(gca,'YTickLabel',[])
    set(gca,'XTick',[])
    set(gca,'YTick',[])
    switch type
     case 1
      zlabel('Voltage Magnitudes [p.u.]')
     case 2
      zlabel('Voltage Angles [deg]')
     case 3
      zlabel('Line Flows [p.u.]')
     case 4
      zlabel('Gen. Rotor Angles [deg]')
     case 5
      zlabel('Gen. Rotor Speeds [p.u.]')
     case 6
      zlabel('Locational Marginal Prices [$/MWh]')
     case 7
      zlabel('Nodal Congestion Prices [$/MWh]')
    end
    colormap(maps{maptype})
    colorbar('EastOutside')
    xlim(x1)
    box on
  end

  hold off

 case 'VoltageReport'

  busidx = find(strcmp(masks,'Bus'));
  for i = 1:Bus.n
    valore = ['|V| = ', ...
              fvar(DAE.y(i+Bus.n),7), ...
              ' p.u.\n<V  = ', ...
              fvar(DAE.y(i),7), ...
              ' rad '];
    set_param(blocks(busidx(i)),'AttributesFormatString',valore);
  end

 case 'WipeVoltageReport'

  busidx = find(strcmp(masks,'Bus'));
  for i = 1:Bus.n,
    set_param(blocks(busidx(i)),'AttributesFormatString','');
  end

 case 'PowerFlowReport'

  simrep(Line, blocks, masks, lines)
  simrep(Hvdc, blocks, masks, lines)
  simrep(Ltc, blocks, masks, lines)
  simrep(Phs, blocks, masks, lines)
  simrep(Lines, blocks, masks, lines)
  simrep(Tg, blocks, masks, lines)
  simrep(Exc, blocks, masks, lines)
  simrep(Oxl, blocks, masks, lines)
  simrep(Pss, blocks, masks, lines)
  simrep(Syn, blocks, masks, lines)
  simrep(PQ, blocks, masks, lines)
  simrep(PV, blocks, masks, lines)
  simrep(SW, blocks, masks, lines)

  %for i = 1:Comp.n
  %  simrep(eval(Comp.names{i},blocks,masks,lines)
  %end

 case 'WipePowerFlowReport'

  for i = 1:length(lines)
    set_param(lines(i),'Name','')
  end

 case 'HideNames'

  nobusidx = find(~strcmp(masks,'Bus'));
  for i = 1:length(nobusidx)
    set_param(blocks(nobusidx(i)),'ShowName','off')
  end
  nobusidx = find(strcmp(masks,''));
  for i = 1:length(nobusidx)
    set_param(blocks(nobusidx(i)),'ShowName','on')
  end

 case 'ShowNames'

  for i = 1:nblock
    set_param(blocks(i),'ShowName','on')
  end

 case 'HideBusNames'

  busidx = find(strcmp(masks,'Bus'));
  for i = 1:Bus.n
    set_param(blocks(busidx(i)),'ShowName','off')
  end

 case 'ShowBusNames'

  busidx = find(strcmp(masks,'Bus'));
  for i = 1:Bus.n
    set_param(blocks(busidx(i)),'ShowName','on')
  end

 case 'FontType'

  for i = 1:nblock
    set_param(blocks(i),'FontName',fontname)
  end
  for i = 1:length(lines)
    set_param(lines(i),'FontName',fontname)
  end

 case 'FontSize'

  for i = 1:nblock,
    set_param(blocks(i),'FontSize',fontsize);
  end
  for i = 1:length(lines),
    set_param(lines(i),'FontSize',fontsize);
  end

 case 'mdl2eps'

  % fontsize == 0:  export to color eps
  % fontsize == 1:  export to grey-scale eps
  grey_eps = fontsize;
  pat = '^/c\d+\s\{.* sr\} bdef';

  cd(Path.data)
  fileeps = [filedata,'.eps'];
  a = dir(fileeps);
  Settings.ok = 1;
  if ~isempty(a)
    uiwait(fm_choice(['Overwrite "',fileeps,'" ?']))
  end
  if ~Settings.ok, return, end
  orient portrait
  print('-s','-depsc',fileeps)

  if ~Settings.noarrows
    cd(Path.local)
    fm_disp(['PSAT model saved in ',Path.data,fileeps])
    return
  end

  file = textread(fileeps,'%s','delimiter','\n');

  idx = [];
  d2 = zeros(1,4);

  for i = 1:length(file)
    if grey_eps & ~isempty(regexp(file{i},pat))
      matchexp = regexp(file{i},'^/c\d+\s\{','match');
      colors = strrep(file{i},matchexp{1},'[');
      colors = regexprep(colors,' sr\} bdef',']');
      rgb = sum(str2num(colors));
      colors = num2str([rgb rgb rgb]/3);
      file{i} = [matchexp{1},colors,' sr} bdef'];
    end
    if strcmp(file{i},'PP')
      if strcmp(file{i-1}(end-1:end),'MP')
        d1 = str2num(strrep(file{i-1},'MP',''));
        if ~d1(3)
          d2(2) = d1(6)-d1(2);
          d2(4) = d1(6)-d1(2);
          d2(1) = d1(5)-d1(2);
          d2(3) = d1(5)+d1(1);
        else
          d2(2) = d1(6)+d1(2);
          d2(4) = d1(6)-d1(1);
          d2(1) = d1(5)-d1(1);
          d2(3) = d1(5)-d1(1);
        end
        file{i-1} = sprintf('%4d %4d mt %4d %4d L',d2);
        idx = [idx, i];
      end
    end
    if ~isempty(findstr(file{i}(max(1,end-1):end),'PO'))
      d1 = str2num(strrep(file{i},'PO',''));
      if ~isempty(findstr(file{i+1},'L'))
        nextline = strrep(file{i+1},'L','');
        d2 = str2num(strrep(nextline,'mt',''));
        if d1(4) == d2(2)
          if d2(1) > d1(3);
            d2(1) = d1(3)-d1(1);
          else
            d2(1) = d2(1)+d1(1)+abs(d1(3)-d2(1));
          end
        else
          if d2(2) > d1(4)
            d2(2) = d1(4)-d1(2);
          else
            d2(2) = d2(2)+d1(2)+abs(d1(4)-d2(2));
          end
        end
        file{i+1} = sprintf('%4d %4d mt %4d %4d L',d2);
      elseif ~isempty(findstr(file{i+2},'MP stroke'))
        d2 = str2num(strrep(file{i+2},'MP stroke',''));
        if d2(1)
          d2(3) = d2(3)-d2(1);
          if d2(1) < 0
            d2(3) = d2(3) - 4;
          end
        end
        if d2(2)
          d2(4) = d2(4)-d2(2);
          if d2(2) < 0
            d2(4) = d2(4) - 4;
          end
        end
        file{i+1} = sprintf('%4d %4d %4d %4d %4d MP stroke',d2);
      end
      idx = [idx, i];
    end
  end

  file(idx) = [];

  fid = fopen(fileeps,'wt+');
  for i = 1:length(file)
    fprintf(fid,'%s\n',file{i});
  end
  fclose(fid);
  cd(Path.local)
  fm_disp(['PSAT model saved in ',Path.data,fileeps])

 otherwise

  fm_disp('Unknown command for Simulink Settings GUI...',2)

end

cd(Path.local)

% -------------------------------------------------------------------------
function [xa,ya] = closerline(xa,ya,xl,yl)

[err,idx] = min(abs(sqrt((xl-xa).^2+(yl-ya).^2)));
if err < 15
  xa = xl(idx);
  ya = yl(idx);
end

% -------------------------------------------------------------------------
function zb = formz(vec,n)

zb = zeros(3*n,1);
zb(1:3:3*n) = vec;
zb(2:3:3*n) = vec;
zb(3:3:3*n) = vec;