psat2epri.m

这是一个很适合研究和学习用的电力系统仿真软件

function check = psat2epri(filename, pathname)
% PSAT2EPRI converts PSAT data file into EPRI Data Format
%
% CHECK = PSAT2EPRI(FILENAME,PATHNAME)
%       FILENAME name of the file to be converted
%       PATHNAME path of the file to be converted
%
%       CHECK = 1 conversion completed
%       CHECK = 0 problem encountered (no data file created)
%
%Author:    Federico Milano
%Date:      06-Oct-2003
%Version:   1.0.0
%
%E-mail:    [e-mail]
%Web-site:  http://thunderbox.uwaterloo.ca/~fmilano
%
% Copyright (C) 2002-2006 Federico Milano

global Settings

if strcmp(pathname(end),filesep)
  pathname = pathname(1:end-1);
end
if ~strcmp(pathname,pwd)
  cd(pathname)
end

fm_disp
fm_disp(['Opening PSAT file "',filename,'"...'])

% General Settings
% -----------------------------------------------------------

check = 1;
b128 = [blanks(128),'\n'];
b12 = blanks(12);

% Defining local data structures
% -----------------------------------------------------------

DAE = struct('a',[],'V',[]);
Bus = struct('con',[],'n',0,'int',[]);
Twt = struct('con',[],'n',0);
Line = struct('con',[],'n',0);
Ltc = struct('con',[],'n',0);
Phs = struct('con',[],'n',0);
Shunt = struct('con',[],'bus',[]);
SW = struct('con',[],'n',0,'bus',[]);
PV = struct('con',[],'n',0,'bus',[]);
PQ = struct('con',[],'n',0,'bus',[]);
Varname = struct('bus','');

% Reading Data from PSAT Data File
% -----------------------------------------------------------

a = exist(filename);
if a == 2,
  eval(filename(1:end-2))
else,
  fm_disp(['File "',pathname,filesep,filename, ...
           '" not found or not an m-file'],2)
  check = 0;
  return
end

% Completing data settings
% -----------------------------------------------------------

% Bus
Bus.n = length(Bus.con(:,1));
busmax = max(Bus.con(:,1));
Bus.int = zeros(busmax,1);
for i = 1:Bus.n, Bus.int(round(Bus.con(i,1))) = i; end
DAE.V = ones(Bus.n,1);
DAE.a = zeros(Bus.n,1);

% Bus names
if isempty(Varname.bus)
  Varname.bus = cellstr(strcat('Bus',num2str(Bus.con(:,1))));
end

% Shunt
if ~isempty(Shunt.con)
  Shunt.bus = Bus.int(Shunt.con(:,1));
  Shunt.con(:,5) = 100*Shunt.con(:,5)./Shunt.con(:,2);
  Shunt.con(:,6) = 100*Shunt.con(:,6)./Shunt.con(:,2);
end

% Three-Winding Transformers
if ~isempty(Twt.con)

  % check data
  [n_twt_rows, n_twt_cols] = size(Twt.con);
  if n_twt_cols < 14
    fm_disp(['Three-winding transformer data does not ', ...
             'seems in a valid format',2])
    return
  elseif n_twt_cols < 15
    Twt.con = [Twt.con, ones(n_twt_rows,1), zeros(n_twt_rows,9)];
  elseif n_twt_cols < 24
    Twt.con = [Twt.con, zeros(n_twt_rows,24-n_twt_cols)];
  end

  % modify Bus.con
  twt_bus = busmax+[1:n_twt_rows]';
  twt_Vn = Twt.con(:,7);
  twt_ar = Bus.con(Bus.int(Twt.con(:,2)),[4 5]);
  Bus.con = [Bus.con; [twt_bus, twt_Vn, ...
                      ones(n_twt_rows,1), ...
                      zeros(n_twt_rows,1), twt_ar]];
  Bus.int = [Bus.int; Bus.n+[1:n_twt_rows]'];
  Bus.n = Bus.n + n_twt_rows;
  busmax = busmax + n_twt_rows;

  % check Varname
  if ~isempty(Varname.bus)
    for i = 1:n_twt_rows
      Varname.bus{end+1,1} = ...
          [Varname.bus{Bus.int(Twt.con(i,1))},'(twt)'];
    end
  end

  % compute branch impedances
  r12 = Twt.con(:,9);
  r13 = Twt.con(:,10);
  r23 = Twt.con(:,11);
  r1 = 0.5*(r12+r13-r23);
  r2 = 0.5*(r12+r23-r13);
  r3 = 0.5*(r23+r13-r12);
  x12 = Twt.con(:,12);
  x13 = Twt.con(:,13);
  x23 = Twt.con(:,14);
  x1 = 0.5*(x12+x13-x23);
  x2 = 0.5*(x12+x23-x13);
  x3 = 0.5*(x23+x13-x12);

  % check for zero impedances in the equivalent star
  idx = find(abs(x1) < 1e-5);
  if ~isempty(idx), x1(idx) = 0.0001; end
  idx = find(abs(x2) < 1e-5);
  if ~isempty(idx), x2(idx) = 0.0001; end
  idx = find(abs(x3) < 1e-5);
  if ~isempty(idx), x3(idx) = 0.0001; end

  % check Line.con column number
  if ~isempty(Line.con)
    ncol = length(Line.con(1,:));
    if ncol < 15,
      Line.con = [Line.con, zeros(Line.n,15-ncol)];
    end
    if ncol > 15,
      Line.con = Line.con(:,[1:15]);
    end
  end

  % modify Line.con
  line1 = [Twt.con(:,1),twt_bus, Twt.con(:,[4,6,5]),zeros(n_twt_rows,1), ...
           Twt.con(:,6)./Twt.con(:,7), r1,x1,zeros(n_twt_rows,1), ...
           Twt.con(:,15), zeros(n_twt_rows,1),Twt.con(:,[17,19,22])];
  line2 = [twt_bus,Twt.con(:,[2,4,7,5]),zeros(n_twt_rows,2), r2,x2, ...
           zeros(n_twt_rows,1), ones(n_twt_rows,1), ...
           zeros(n_twt_rows,1),Twt.con(:,[18,20,23])];
  line3 = [twt_bus,Twt.con(:,[3,4,7,5]),zeros(n_twt_rows,1), ...
           Twt.con(:,7)./Twt.con(:,8),r3,x3, zeros(n_twt_rows,1), ...
           ones(n_twt_rows,1), zeros(n_twt_rows,1),Twt.con(:,[19,21,24])];
  Line.con = [Line.con; line1; line2; line3];

  % clear Twt.con
  Twt.con = [];
end

% take into account initial guess
if length(Bus.con(1,:) > 2) & max(Bus.con(:,3) < 2)
  DAE.V = Bus.con(:,3);
end
if length(Bus.con(1,:) > 3) & max(Bus.con(:,4) < pi)
  DAE.a = Bus.con(:,4);
end

%  Line
if isempty(Line.con)
  fm_disp(['The selected data file has no Lines, which ', ...
           'are required in the WSCC format.'],2)
  check = 0;
  return
else
  Line.n = length(Line.con(:,1));
  Line.con(:,1) = Bus.int(round(Line.con(:,1)));
  Line.con(:,2) = Bus.int(round(Line.con(:,2)));
  % parameters in p.u. with respect to system base
  Line.con(:,8)  = Line.con(:,8)./Line.con(:,3)*Settings.mva;
  Line.con(:,9)  = Line.con(:,9)./Line.con(:,3)*Settings.mva;
  Line.con(:,10) = Line.con(:,10).*Line.con(:,3)/Settings.mva;
  % check for paramter units (if the length is > 0, unit is km)
  idx = find(Line.con(:,6));
  if ~isempty(idx)
    XB = Line.con(idx,4).*Line.con(idx,4)./Line.con(idx,3);
    Line.con(idx,8)  = Line.con(idx,6).*Line.con(idx,8)./XB;
    Line.con(idx,9)  = Line.con(idx,6).*Line.con(idx,9)./XB;
    Line.con(idx,10) = Line.con(idx,6).*Line.con(idx,10).*XB;
    bidx = find(Line.con(idx,10) > 10);
    if ~isempty(bidx)
      fm_disp('Warning: Some line susceptances are too high...')
      fm_disp(['         microFarad are assumed for those ', ...
               'susceptance values'])
      Line.con(idx(bidx),10) = Line.con(idx(bidx),10)/1e6;
    end
  end
end

% Swing generators
if isempty(SW.con)
  fm_disp(['The selected data file has no slack bus, which ', ...
           'is required in the WSCC fromat.'],2)
  check = 0;
  return
else
  SW.n = length(SW.con(:,1));
  SW.bus = Bus.int(round(SW.con(:,1)));
  if length(SW.con(1,:)) < 8
    SW.con = [SW.con(:,1),100*ones(SW.n,1),ones(SW.n,1),SW.con(:,[2:end])];
  end
  if length(SW.con(1,:)) == 5
    SW.con = [SW.con, 999*ones(SW.n,1), -999*ones(SW.n,1), ...
	      1.1*ones(SW.n,1), 0.9*ones(SW.n,1)];
  end
  DAE.V(SW.bus) = SW.con(:,4);
  DAE.a(SW.bus) = SW.con(:,5);
  SW.con(:,6)  = SW.con(:,6).*SW.con(:,2)/100;
  SW.con(:,7)  = SW.con(:,7).*SW.con(:,2)/100;
  SW.con(:,10) = SW.con(:,10).*SW.con(:,2)/100;
end

% PV generators
if ~isempty(PV.con)
  PV.n = length(PV.con(:,1));
  PV.bus = Bus.int(round(PV.con(:,1)));
  DAE.V(PV.bus) = PV.con(:,5);
  if length(PV.con(1,:)) == 5
    PV.con = [PV.con, zeros(PV.n,2), ...
              1.1*ones(PV.n,1), 0.9*ones(PV.n,1), ones(PV.n,1)];
  end
  PV.con(:,4) = PV.con(:,4).*PV.con(:,2)/100;
  PV.con(:,6) = PV.con(:,6).*PV.con(:,2)/100;
  PV.con(:,7) = PV.con(:,7).*PV.con(:,2)/100;
end

% PQ loads
if ~isempty(PQ.con)
  PQ.n = length(PQ.con(:,1));
  PQ.bus = Bus.int(round(PQ.con(:,1)));
  if length(PQ.con(1,:)) == 5
    PQ.con = [PQ.con, zeros(PQ.n,3)];
  elseif length(PQ.con(1,:)) == 7
    PQ.con = [PQ.con, zeros(PQ.n,1)];
  end
  for i = 1:PQ.n
    if PQ.con(i,6) == 0, PQ.con(i,6) = 1.2; end
    if PQ.con(i,7) == 0, PQ.con(i,7) = 0.8; end
  end
  PQ.con(:,4) = PQ.con(:,4).*PQ.con(:,2)/100;
  PQ.con(:,5) = PQ.con(:,5).*PQ.con(:,2)/100;
end

% Regulating transformers
if ~isempty(Ltc.con)
  Ltc.n = length(Ltc.con(:,1));
  Ltc.bus1 = Bus.int(round(Ltc.con(:,1)));
  Ltc.bus2 = Bus.int(round(Ltc.con(:,2)));
  Ltc.busc = Ltc.bus2;
  idx = find(Ltc.con(:,15));
  if ~isempty(idx)
    Ltc.busc(idx) = Bus.int(round(Ltc.con(idx,15)));
  end
end
if ~isempty(Phs.con)
  Phs.n = length(Phs.con(:,1));
  Phs.bus1 = Bus.int(round(Phs.con(:,1)));
  Phs.bus2 = Bus.int(round(Phs.con(:,2)));
end

% Opening File
% -----------------------------------------------------------

newfile = [filename(1:end-2),'.wsc'];
fm_disp(['Writing WSCC file "',newfile,'"...'])
fid = fopen([pathname,filesep, newfile], 'wt');
comment = ['C\nC',repmat('*',1,79),'\nC\n'];
count = fprintf(fid,comment);

% Header and Title
% -----------------------------------------------------------

count = fprintf(fid,'HDG\n');
count = fprintf(fid,['PSAT ARCHIVE\n']);
count = fprintf(fid,[num2str(Bus.n),'-Bus ', ...