cdf2matp.m

关于电力系统的计算程序!最新的改进法计算潮流程序！用matlab实现的！ 适用于电力行业的相关从业人员！

function cdf2matp
% CDF2MATP reads in load flow data written in IEEE Common Format, outputs 
% load flow data and OPF data suitable for matpower. 
%
% Since a typical set of system data in IEEE Common Format does not 
% contain OPF data format, it sets up generator cost data.
% Cost curves are modeled by quadratic functions with default coefficients. 
% You can change the default parameters if you want. Cost curve can 
% also be modelled by piece-wise linear function. The formats of 
% cost curve data are given in CASE.M.
%
% CDF2MATP also sets up the following data which are necessary to 
% run optimal power flow:
% 
% generation upper bound = Pg + baseMVA
% generation lower bound = 0
% voltage upper bound = 1.06
% lower bounds = 0.94
%
% CDF2MATP may modify some of the data which are "infeasible" for 
% running optimal power flow. If so, warning information will be 
% printed out on screen.
%
% Note: since our code can not handle transformers with variable tap, 
% you may not expect to get exactly the same power flow solution 
% using converted data. This is the case when we converted
% ieee300.cdf

%   MATPOWER Version 2.0
%   by Deqiang (David) Gan, PSERC Cornell    12/22/97
%   Copyright (c) 1996, 1997 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/ for more info.

% ----- get the original load flow data -----

cdffilename =  input('please enter IEEE Common Data Format file name --- ', 's');
matpfilename = input('please enter MATPOWER file name                --- ', 's');

fid = fopen(cdffilename, 'r');


% get baseMVA

title = fgetl(fid);
baseMVA = str2num(title(32:37));



% find string 'BUS DATA FOLLOWS'
while 1	
	line = fgetl(fid);
        if line(1:16) == 'BUS DATA FOLLOWS', break, end;
end;

% ----- get bus data, feed them into matrix bus, gen, gencost, and area

ibus = 0;
igen = 0;
iarea = 0;
area=[];

while 1
	line = fgetl(fid);
	if line(1:4) == '-999', break, end;

	% feed bus data
	ibus = ibus + 1;
	bus(ibus, 1) = str2num(line(1:4)); 							% bus number
	bus(ibus, 2) = str2num(line(25:26)); if bus(ibus, 2) == 0, bus(ibus, 2) = 1; end;	% bus type
	bus(ibus, 3) = str2num(line(41:49));							% Pd
	bus(ibus, 4) = str2num(line(50:59));							% Qd
	bus(ibus, 5) = str2num(line(107:114));							% Gs
	bus(ibus, 6) = str2num(line(115:122));							% Bs
	if abs(bus(ibus, 6)) > 0.0000001, bus(ibus, 6) = baseMVA*bus(ibus, 6); end;
	bus(ibus, 7) = str2num(line(19:20));							% area
	bus(ibus, 8) = str2num(line(28:33));							% Vm
	bus(ibus, 9) = str2num(line(34:40));							% Va
	bus(ibus, 10) = str2num(line(77:83));							% baseKV
	bus(ibus, 11) = str2num(line(21:23));							% zone
	bus(ibus, 12) = 1.06;									% default voltage upper limit
	bus(ibus, 13) = 0.94;									% default voltage lower limit

	% feed gen and gencost
	Pg = str2num(line(60:67)); 
	Qg = str2num(line(68:75));
	Qmax = str2num(line(91:98)); 
	Qmin = str2num(line(99:106));
	if Pg ~=0 | Qg ~=0 | Qmax ~=0 | Qmin ~=0 | bus(ibus, 2) == 3
		igen = igen + 1;
		if bus(ibus, 2) == 3, refgen = igen; end;
		gen(igen, 1) = bus(ibus, 1);						% bus number
		gen(igen, 2) = Pg;							% Pg
		if gen(igen, 2) < 0
			bus(ibus, 3) = bus(ibus, 3) - gen(igen, 2);			% negtive Pg is transformed as load 
			fprintf('\n ***** negtive Pg at bus %g treated as Pd', bus(ibus, 1));
% 			fprintf('\n ***** attention: negtive Pg at bus %g has been transformed into Pd \', bus(ibus, 1));
% 			fprintf('\n       the transformation should not affect the solution\n');
			gen(igen, 2) = 0;
		end;
		gen(igen, 3) = Qg;							% Qg
		gen(igen, 4) = Qmax;							% Qmax
		gen(igen, 5) = Qmin;							% Qmin
		if Qmax - Qmin < 0.01, 
			gen(igen, 4) = Qmin + 0.1 * baseMVA;				% Qmax is modified
			fprintf('\n ***** Qmax = Qmin at generator at bus %4i (Qmax set to Qmin + %g)', bus(ibus, 1), baseMVA/10); 
% 			fprintf('\n ***** attention: Qmax - Qmin = 0 at generator connected to  bus %4i \n', bus(ibus, 1)); 
% 			fprintf('\n       Qmax has been set to Qmin + 0.1 * %g by this converter', baseMVA);
% 			fprintf('\n       you can open file %s, change Qmax in matrix "gen" \n ', matpfilename);
		end;
		gen(igen, 6) = str2num(line(85:90));					% specified voltage
		gen(igen, 7) = baseMVA;							% baseMVA
		gen(igen, 8) = 1;							% default status is 'on'
		gen(igen, 9) = gen(igen, 2) + baseMVA;					% Pmax
		gen(igen, 10) = 0;							% Pmin = 0 by default




		gencost(igen, 1) = 2;						% by default, sets the model as polynomial
		gencost(igen, 2) = 0; 						% start up cost is zero by default
		gencost(igen, 3) = 0;						% shut down cost is zero by default
		gencost(igen, 4) = 3;						% number of coefficients in polynomial cost
%		gencost(igen, 5) = 0.01;					% default c2
%		gencost(igen, 6) = 0.3;						% default c1
%		gencost(igen, 7) = 0.2;						% default c0


	end;

	% feed area data                   ***** note: this part of code is under construction
	if isempty(area) == 1
		fprintf('\n ***** area data conversion not yet implemented (creating dummy area data)');
% 		fprintf('\n ***** attention: the area data can not be properly processed by this version of converter\');
		iarea = iarea +1;
		area(iarea, 1) = bus(ibus, 7);
		area(iarea, 2) = bus(ibus, 1);
	end;

end;

totload = sum(bus(:, 3)); 
totgen = sum(gen(:, 2));
if totgen < 1.04 * totload
	gen(refgen, 9) = gen(refgen, 2) + 1.1 * totload - totgen;			% Pg at slack bus is modified
	fprintf('\n ***** Insufficient generation, setting Pmax at slack bus (bus %d) to %g', gen(refgen, [1, 9]));
% 	fprintf('\n ***** attention: Pmax at slack bus %i is overly small to get feasible solutions\n', gen(refgen, 1));
% 	fprintf('\n       Pg has been set to %g by this converter', gen(refgen, 9));
% 	fprintf('\n       you can open file %s, change Pg in matrix "gen" \n', matpfilename);
end;



% ----- set up the cost coefficients of generators

gencost(:, 5) = zeros(size(gen, 1), 1);
gencost(:, 6) = 100*ones(size(gen, 1), 1) ./ (gen(:, 2) + 10*ones(size(gen, 1), 1));
gencost(:, 7) = 100*ones(size(gen, 1), 1) ./ (gen(:, 2) + 10*ones(size(gen, 1), 1));





% find string 'BRANCH DATA FOLLOWS'

while 1	
	line = fgetl(fid);
        if line(1:19) == 'BRANCH DATA FOLLOWS', break, end;
end;

% ----- get branch data, feed them into matrix branch
k = 0;
while 1
	line = fgetl(fid);
	if line(1:4) == '-999', break, end;

	k = k + 1;
	branch(k, 1) = str2num(line(1:4)); 							% fbus (also the tap bus)
	branch(k, 2) = str2num(line(6:9)); 							% tbus
	branch(k, 3) = str2num(line(20:29)); 							% R
	branch(k, 4) = str2num(line(30:40)); 							% X
	branch(k, 5) = str2num(line(41:50)); 							% B
	branch(k, 6) = str2num(line(51:55));							% RATE A
	if branch(k, 6) < 0.000001, 
		branch(k, 6) = 99 * baseMVA; 							% RATE A is modified
		fprintf('\n ***** MVA limit of branch %d - %d not given, set to %g', branch(k, [1, 2, 6]));
% 		fprintf('\n ***** attention: thermal limit of branch %g -- %g is not give\', branch(k, 1), branch(k, 2));
% 		fprintf('\n       it has been set to 99 * %g by this converter\', baseMVA);
% 		fprintf('\n       you can open %s, change the branch data if you want\', matpfilename);
	end;
	branch(k, 7) = str2num(line(57:61));							% RATE B
	branch(k, 8) = str2num(line(63:67));							% RATE C
	branch(k, 9) = str2num(line(77:82));							% transformer turns ratio
	branch(k, 10) = 0;									% phase shifter can not be modelled
	branch(k, 11) = 1;									% by default, branch is on
end;
fprintf('\n');
fclose(fid);

% ----- write data in MATPOWER format -----
fid = fopen(matpfilename, 'w');

% print out function title, baseMVA, etc.
temp = matpfilename(1:length(matpfilename)-2);
fprintf(fid, 'function [baseMVA, bus, gen, branch, area, gencost] = %s', temp);
fprintf(fid, '\n\nbaseMVA = %5.1f; \n', baseMVA);

% print out bus matrix
fprintf(fid, '\n%%bus type  Pd      Qd       Gs      Bs    area  Vm      Va      baseKV zone   Vmax   Vmin');
fprintf(fid, '\nbus = [');
for i = 1:size(bus, 1)
		%       bus type Pd    Qd    Gs    Bs    area Vm    Va    baseKV zone Vmax  Vmin
	fprintf(fid, '\n%4i %2i  %8.2f %8.2f %8.3f %8.3f %3i  %8.5f %8.3f %8.2f  %i   %8.5f %8.5f', bus(i,1:13));
end
fprintf(fid, '\n];\n');

% print out gen matrix
fprintf(fid, '\n%%bus  Pg       Qg      Qmax    Qmin   Vsp      base  status   Pmax    Pmin');
fprintf(fid, '\ngen = [');
for i = 1:size(gen, 1)
		%	bus Pg    Qg    Qmax  Qmin  Vsp   base  status Pmax  Pmin
	fprintf(fid, '\n%4i %8.2f %8.2f %8.2f %8.2f %8.5f %8.2f %4i    %8.2f %8.2f', gen(i,1:10));
end
fprintf(fid, '\n];\n');

% print out branch matrix
fprintf(fid, '\n%%fbus tbus      r     x     b       ratea   rateb   ratec ratio   angle  statu');
fprintf(fid, '\nbranch = [');
for i = 1:size(branch, 1)
		%	f   t    r     x     b     ratea rateb  ratec ratio angle  status
	fprintf(fid, '\n%4i %4i  %8.5f %8.5f %8.5f %8.2f %8.2f  %8.2f %8.5f %8.3f  %4i', branch(i, 1:11));
end
fprintf(fid, '\n];\n');

% print out area matrix
if isempty(area) ~= 1
	fprintf(fid, '\narea = [');
	for i = 1:size(area, 1)
		fprintf(fid, '\n%4i   %4i', area(i, 1:2));
	end;
	fprintf(fid, '\n];	\n');
end

% print out gencost matrix
if isempty(gencost) ~= 1
	fprintf(fid, '\ngencost = [');
	for i = 1:size(gen, 1)
			%       i   startup  shut down  n   c2    c1    c0
		fprintf(fid, '\n%4i %8.2f    %8.2f      %4i %g    %g    %g', gencost(i,1:7));
	end;
	fprintf(fid, '\n];\n');
end;

% print out ancillary stuff
fprintf(fid, '\nreturn;	\n');

fclose(fid);

% things need to be done

% area is not properly set up, this should not affect most of applications
% PTI format can not be handled