fm_wscc.m

电力系统的psat

function tutto_ok = fm_wscc(dfile,pathname,flag)
% FM_WSCC convert BPA-WSCC data files in PSAT data files
%
% OK = FM_WSCC(FILE,PATH,FLAG)
%      FILE name of the file to be converted
%      PATH path of the file to be converted
%      FLAG append orignal data file [1,0]
%      OK check  1 -> conversion successfully completed
%                0 -> conversion failed
%
%Author:    Federico Milano
%Date:      11-Nov-2002
%Update:    08-Oct-2003
%Update:    23-July-2004
%Version:   1.1.1
%
%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 Settings Path Varname

fm_disp(' ',1);
fm_disp('Conversion From BPA-WSCC Data Format...',1);
fm_disp(['Source Data File "',pathname,filesep,dfile,'"'],1)

tipi = length(Varname.fcomp);
type_con = cell(tipi,1);
for i = 1:tipi
  type_con{i,1} = [Varname.fcomp{i},'_con'];
  eval([type_con{i,1},' = [];'])
end

tutto_ok = 1;
fidtxt = fopen([pathname, filesep, dfile]);
b80 = blanks(80);


% header
while 1
  sline = fgetl(fidtxt);
  if sline == -1
    fm_disp(['The Data File is not in a correct WSCC Format. ', ...
	     'Conversion Process Aborted'],2);
    tutto_ok = 0;
    break
    return
  end
  % EPRI format
  if ~isempty(strfind(sline,'HDG'))
    intestazione = fgetl(fidtxt);
    intestazione = [intestazione,'\n%%',fgetl(fidtxt)];
    intestazione = [intestazione,'\n%%',fgetl(fidtxt)];
    break
  end
  % BPA/China format
  if ~isempty(strfind(sline,'CASEID='))
    idxeq = findstr(sline,'=');
    idxcom = findstr(sline,',');
    intestazione = ['CASE: ',sline(idxeq(1):idxcom(2)), ...
                    '\n%% PROJECT: ',sline(idxeq(2):end)];
    break
  end
end

% data matrices
Bus_n = 0;
Line_n = 0;
ntrasf = 0;
pv_idx = [];
while 1
  sline = fgetl(fidtxt);
  if sline == -1, break, end
  if ~isempty(findstr(sline,'BAS')) | ...
        ~isempty(findstr(sline,'/NETWORK_DATA\'))
    break
  end
end
while 1
  sline = fgetl(fidtxt);
  if isempty(sline), sline = ' '; end
  if sline == -1, break, end
  sline = [sline, b80];

  if strcmp(sline(1),'C')
    % comment line in EPRI format
  elseif strcmp(sline(1),'.')
    % comment line in BPA/China format
  elseif strcmp(sline(1:4),'-999')
    % end of bus data in BPA/China format
  elseif strcmp(sline(1),'B') & ~strcmp(sline([1 2]),'BD') ...
	& ~strcmp(sline([1 2]),'BZ')
    % bus data
    Bus_n = Bus_n + 1;
    busname{Bus_n,1} = sline(7:14);
    kV(Bus_n,1) = str2num(sline(15:18));
    Z(Bus_n,1) =  getv(sline,19,20);
    PL(Bus_n,1) = getv(sline,21,25);
    QL(Bus_n,1) = getv(sline,26,30);
    MWshunt(Bus_n,1) =  getv(sline,31,34);
    Mvashunt(Bus_n,1) = getv(sline,35,38);
    PM(Bus_n,1) =  getv(sline,39,42);
    P(Bus_n,1) =   getv(sline,43,47);
    QM(Bus_n,1) =  getv(sline,48,52,100);
    Qm(Bus_n,1) =  getv(sline,53,57,100);
    Vpu(Bus_n,1) = getv(sline,58,61,1000,10);
    Vm(Bus_n,1) =  getv(sline,62,65);
    RemoteBusName(Bus_n,1) = getv(sline,66,73);
    RemoteBuskV(Bus_n,1) = getv(sline,74,77);
    Qperc(Bus_n,1) = getv(sline,78,80);

    if strcmp(sline(2),'Q')
      % PV generator with Q limits
      pv_idx = [pv_idx; Bus_n];
    elseif strcmp(sline(2),'E')
      % PV generator with no Q limits
      pv_idx = [pv_idx; Bus_n];
    elseif strcmp(sline(2),'V')
      % PQ generator bus with DAE.V limits
    elseif strcmp(sline(2),'G')
      % PQ generator bus with Q limits controlling
      % voltage on a remote PV load bus
    elseif strcmp(sline(2),'C')
      % PV load bus with remote voltage control
    elseif strcmp(sline(2),'T')
      % PQ load bus with voltage controled by LTC transformer
    elseif strcmp(sline(2),'S')
      % Swing bus
      slack_name = sline(7:14);
      slack_kV = str2num(sline(15:18));
      slack_angle = 0.00;
    end

  elseif strcmp(sline(1),'L') & ~strcmp(sline([1 2]),'LD')

    % line data
    Line_n = Line_n + 1;
    Lbus1name{Line_n,1} = sline(7:14);
    LkV1(Line_n,1) = str2num(sline(15:18));
    Lbus2name{Line_n,1} = sline(20:27);
    LkV2(Line_n,1) = str2num(sline(28:31));
    LIn(Line_n,1) = getv(sline,34,37,1e5,1,1);
    LR(Line_n,1)   = getv(sline,39,44,1e5,1,1);
    LX(Line_n,1)   = getv(sline,45,50,1e5,1,1);
    LG2(Line_n,1)  = getv(sline,51,56,1e5,3,1);
    LB2(Line_n,1)  = getv(sline,57,62,1e5,10,1);
    LMil(Line_n,1) = getv(sline,63,66);
    if ~LMil(Line_n,1), LMil(Line_n,1) = 1; end

  elseif strcmp(sline(1),'T')

    % transformer data
    ntrasf = ntrasf + 1;
    Tbus1name{ntrasf,1} = sline(7:14);
    TkV1(ntrasf,1) = str2num(sline(15:18));
    Tbus2name{ntrasf,1} = sline(20:27);
    TkV2(ntrasf,1) = str2num(sline(28:31));

    TSn(ntrasf,1) = getv(sline,34,38);
    TR(ntrasf,1) = getv(sline,39,44,1e5,1,1);
    TX(ntrasf,1) = getv(sline,45,50,1e5,1,1);
    TG(ntrasf,1) = getv(sline,51,56);
    TB(ntrasf,1) = getv(sline,57,62);

    % Tap ratio 1
    a = sline(63:67);
    b = strrep(a(1:3),'.','');
    c = strrep(a(4:5),'.','');
    a = [b,'.',c];
    a = str2num(a);
    if isempty(a), a = 1; end
    Tap1(ntrasf,1) = a;

    % Tap ratio 2
    a = sline(68:72);
    b = strrep(a(1:3),'.','');
    c = strrep(a(4:5),'.','');
    a = [b,'.',c];
    a = str2num(a);
    if isempty(a), a = 1; end
    Tap2(ntrasf,1) = a;

  elseif strcmp(sline(1:2),'RT')
    % ULTC data
    fm_disp('Skipping regulating transformer data ...')

  elseif strcmp(sline(1:3),'SOL')
    % slack bus data (EPRI format only)
    slack_name = sline(31:38);
    slack_kV = str2num(sline(39:45));
    slack_angle = str2num(sline(46:55));
  elseif strcmp(sline(1:5),'(END)')
    % end of data (BPA/China only)
  end
end

if flag == 1
  status = fseek(fidtxt,0,'bof');
  orig_file = cell(1,1);
  linea = 1;
  while 1
    riga = fgetl(fidtxt);
    if ischar(riga)
      orig_file{linea,1} = ['%% ', riga, '\n'];
      linea = linea + 1;
    else
      break
    end
  end
end

fclose(fidtxt);

% build matrix "Bus.con"
bus = [[1:Bus_n]',kV];

% build matrix "SW.con"
for i = 1:Bus_n
  if strcmp(busname{i},slack_name)
    sw_idx = i;
    break
  end
end
SW_con = [sw_idx, 100*ones(length(sw_idx),1), kV(sw_idx), ...
	  Vpu(sw_idx), slack_angle, QM(sw_idx), Qm(sw_idx), ...
	  1.2*Vpu(sw_idx),Vm(sw_idx), P(sw_idx)/100, ...
	  ones(length(sw_idx),1)];

% build matrix "PQ.con"
pq_idx = find(PL | QL);

PQ_con = [pq_idx, 100*ones(length(pq_idx),1), ...
	  kV(pq_idx), PL(pq_idx)/100, ...
	  QL(pq_idx)/100, Vpu(pq_idx), Vm(pq_idx), ...
	  zeros(length(pq_idx),1)];

% build matrix "PV.con"
PV_con = [pv_idx, 100*ones(length(pv_idx),1), ...
	  kV(pv_idx), P(pv_idx)/100,Vpu(pv_idx), ...
	  QM(pv_idx), Qm(pv_idx),1.2*Vpu(pv_idx), ...
	  Vm(pv_idx), ones(length(pv_idx),1)];

% build matrix "Shunt.con"
shunt_idx = find(MWshunt | Mvashunt);

Shunt_con = [shunt_idx, 100*ones(length(shunt_idx),1), ...
	     kV(shunt_idx),60*ones(length(shunt_idx),1), ...
	     MWshunt(shunt_idx)/100, Mvashunt(shunt_idx)/100];

% build matrix "Line.con"
from_idx = [];
for i = 1:Line_n
  for j = 1:Bus_n
    if strcmp(Lbus1name{i},busname{j})
      from_idx = [from_idx; j];
      break
    end
  end
end
to_idx = [];
for i = 1:Line_n
  for j = 1:Bus_n
    if strcmp(Lbus2name{i},busname{j})
      to_idx = [to_idx; j];
      break
    end
  end
end
Line_con = [from_idx, to_idx, 100*ones(Line_n,1), LkV1, ...
            60*ones(Line_n,1), zeros(Line_n,1), zeros(Line_n,1), ...
            LR.*LMil, LX.*LMil, 2*LB2.*LMil, zeros(Line_n,2), ...
            LIn./(100000./(sqrt(3)*LkV1))];
from_idx = [];
for i = 1:ntrasf
  for j = 1:Bus_n
    if strcmp(Tbus1name{i},busname{j})
      from_idx = [from_idx; j];
      break
    end
  end
end
to_idx = [];
for i = 1:ntrasf
  for j = 1:Bus_n
    if strcmp(Tbus2name{i},busname{j})
      to_idx = [to_idx; j];
      break
    end
  end
end
if ntrasf > 0
  tap_ratio = 1./(TkV1./TkV2./Tap1.*Tap2);
  Line_con = [Line_con; ...
	      [from_idx, to_idx, 100*ones(ntrasf,1),TkV1, ...
               60*ones(ntrasf,1), zeros(ntrasf,1),TkV1./TkV2, ...
               TR, TX, TB, tap_ratio, zeros(ntrasf,1),TSn]];
end

lfile = length(dfile)-4;
nfile = dfile(1:lfile);
newfile = ['d_', nfile, '.m'];
newfile = strrep(newfile,'d_d_','d_');
newfile = strrep(newfile,'d_d','d_');
newpath = pathname;

finenf = newfile(length(newfile)-1:length(newfile));
if ~strcmp(finenf, '.m')
  newfile = [newfile, '.m'];
end
if ~strcmp(newfile(1), 'd')
  newfile(1) = 'd';
end
fid = fopen([newpath,filesep, newfile], 'wt');

% write PSAT data file
count = fprintf(fid, ['%%  ', intestazione, '\n\n']);
count = fprintf(fid,'Bus.con = [ ...\n      ');
for i = 1:Bus_n-1
  count = fprintf(fid, '%4d %8.4g; ',bus(i,:));
  if rem(i,5) == 0
    count = fprintf(fid,'\n      ');
  end
end
count = fprintf(fid, '%4d %8.4g]; \n\n\n', bus(end,:));

for j = 2:length(type_con)
  typecon = type_con{j,1};
  format = Varname.format{j};
  if eval(['~isempty(',typecon,')'])
    count = fprintf(fid, [strrep(typecon,'_','.'), ' = [ ...\n']);
    eval(['quali = ',typecon,';']);
    if ~isempty(quali)
      count = fprintf(fid, [format,';\n'],quali(1:end-1,:)');
      count = fprintf(fid, [format,'];\n\n\n'],quali(end,:));
    end
  end
end

% write "Varname.bus"
count = fprintf(fid, 'Varname.bus = {...\n      ');
for i = 1:Bus_n-1
  count = fprintf(fid, ['''', busname{i,1},'''; ']);
  if rem(i,5) == 0
    count = fprintf(fid,'\n      ');
  end
end
if iscell(busname)
  count = fprintf(fid, ['''', busname{length(busname),1},'''};\n\n']);
else
  count = fprintf(fid, ['''', busname,'''};\n\n']);
end

if flag == 1
  count = fprintf(fid, '%% ORIGINAL DATA FILE: \n\n');
  for i = 1:length(orig_file)
    count = fprintf(fid, '%s\n', orig_file{i,1});
  end
end
count = fclose(fid);

fm_disp(['Conversion into data file "', ...
	 pathname,filesep,newfile,'" completed.'],1)
if Settings.beep, disp(char(7)), end


%---------------------------------------------------------------------
function value = getv(sline,int1,int2,base,check,cabs)

value = str2num(sline(int1:int2));

if isempty(value), value = 0; end

switch nargin
 case 6, if abs(value) > check, value = value/base; end
 case 5, if value > check, value = value/base; end
 case 4, value = value/base;
end