pfinpt.m

电力系统Matlab潮流程序。 根据给定输入数据

clear all, clc, close all

%      PF data in  MATLAB format
%			IEEE 14 Bus  Test System
Sb = 100; 
econv = 0.001; 
mxitr = 10;
%			Bus Data
%	Bus   Pg (MW)   Qg(Mvar)   Pl (MW)  Ql (Mvar)    Kv (nom.)  kV(reg) 
%   1     2         3          4        5            6          7 
bus = [
	 1 	0	0	0	0	230	240
	12	40	0	20	10	230	230
	13	60	0	40	20	230	230
	14	35	0	30	10	230	230
	15	20	0	5	3	130	130
	2	0	0	8	2	230	0
	3	0	0	0	0	230	0
	4	0	0	9	3	130	0
	5	0	0	15	10	130	0
	6	0	0	25	15	130	0
	7	0	0	0	0	130	0
	8	0	0	10	6	130	0
	9	0	0	0	0	130	0
	10	0	0	15	6	130	0
	11	0	0	0	0	230	0   ];

%			Branch Data
%	Fr	To        r (p.u)    x (p.u)        b (p.u)    t?(tap ratio)
%   1   2         3          4              5          6
brnch = [
	1	2	.005	.010	0.02	0
	1	12	.002	.006	0.03	0
	2	3	.002	.004	0	    0
	2	12	.016	.120	0.03	0
	2	15	.000	.085	0	    1.00
	3	4	.000	.080	0	    1.03
	3	12	.016	.080    0       0
	3	13	.004	.010    0       0
	3	14	.001	.004    0       0
	4	6	.005	.040    0       0
	4	7	.003	.050    0       0
	5	8	.015	.120    0       0
	5	6	.045	.080    0       0
	5	15	.010	.060    0       0
	7	8	.070	.120    0       0
	8	9	.015	.060    0       0
	9	10	.000	.030    0       0
	11	10	.000	.060	0	    1.03
	12	13	.016	.087	0.02    0
	14	11	.036	.130	0       0	];


% Begin program %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define main variable
% BusNum,1 
% Pg (MW),2
% Qg(Mvar),3
% Pl(MW),4
% Ql(Mvar),5
% Kv(nom.),6
% kV(reg),7
% Type,8
% v,9
% delta,10
% p,11
% q,12
Bus = []; % processed bus data
Branch = []; % processed branch data
YBus = []; % YBus
Bprime = [];
Bprimeprime = [];
Tolerance = econv;
Branch = brnch;
NumBus = size(bus,1);
NumPQ = 0;
NumPV = 0;
MaxNumIter = mxitr;

'preprocess data ...'
for i=1:size(bus,1)
    for j=1:size(bus,1)
        if (bus(j,1)==i)
            Bus(i,:)=bus(j,:);
        end
    end
end

% define bus type and get per-unit value
for i=1:size(Bus,1)
    % get per-unit
    Bus(i,2) = Bus(i,2)/Sb;
    Bus(i,3) = Bus(i,3)/Sb;
    Bus(i,4) = Bus(i,4)/Sb;
    Bus(i,5) = Bus(i,5)/Sb;
    
    % bus type, 8th colunm
    if (i==1)
        Bus(i,8)=3; %swing bus
        Bus(i,9)=Bus(i,7)/Bus(i,6); % v given
        Bus(i,10)=0; % delta given
        Bus(i,11)=0;
        Bus(i,12)=0;
    else
        if (Bus(i,7) ~=0)
            Bus(i,8)=2; %PV bus
            Bus(i,9)=Bus(i,7)/Bus(i,6); %v
            Bus(i,10)=0; % delta 
            Bus(i,11)=Bus(i,2)-Bus(i,4); %p
            Bus(i,12)=0; %q
            NumPV = NumPV + 1;
        else
            Bus(i,8)=1; %PQ bus
            Bus(i,9)=1; %v
            Bus(i,10)=0; % delta 
            Bus(i,11)=Bus(i,2)-Bus(i,4); %p
            Bus(i,12)=Bus(i,3)-Bus(i,5); %q
            NumPQ = NumPQ + 1;
        end
    end 
end

% Build YBus %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%begin
'Building YBus ...'
% initialize
for i=1:size(Bus,1)
    for j=1:size(Bus,1)
        YBus(i,j)=0;
    end
end 

% find diagonal elemen, YY
for i=1:size(Bus,1)
    curbus=Bus(i,1);
    for j=1:size(Branch,1) 
        curfrombus=Branch(j,1);
        curtobus=Branch(j,2);
        % either from bus or to bus
        if ((Bus(i,1)==Branch(j,1)) | (Bus(i,1)==Branch(j,2)))
            y = 1/complex(Branch(j,3),Branch(j,4));
            if (Branch(j,6)~=0) % tap changing xfrm
                if (Bus(i,1)==Branch(j,1)) % from bus
                    y = y * (Branch(j,6)*Branch(j,6)); % t*t*y
                end
            end
            
            YBus(i,i) = YBus(i,i) + y + complex(0,Branch(j,5)/2); % shunt
        end
    end
end

% form non-diagonal element
for j=1:size(Branch,1)
    y = -1/complex(Branch(j,3),Branch(j,4));
    if (Branch(j,6)~=0) 
        y = y*Branch(j,6);
    end
    YBus(Branch(j,1),Branch(j,2)) = y;
    YBus(Branch(j,2),Branch(j,1)) = y;
end
% Build YBus %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end


'Form B prime and B double prime ...'
Bprime=imag(YBus);
Bprimeprime=imag(YBus);
% B prime includes n-1 dimension, no swing bus
for i=1:size(Bus,1)
    if (Bus(i,8)==3) % swing bus
        Bprime(Bus(i,1),:)=[];
        Bprime(:,Bus(i,1))=[];
        break;
    end
end

% B prime prime includes n-m-1 dimension, no swing bus and PV bus
for i=1:size(Bus,1)
    if (Bus(i,8)~=1) % swing bus or PV bus
        Bprimeprime(Bus(i,1),:)=0;
        Bprimeprime(:,Bus(i,1))=0;
    end
end
Temp=diag(Bprimeprime);
while (size(find(Temp),1) ~= size(Temp,1))
    for i=1:size(Bprimeprime,1)
        if (Bprimeprime(i,i)==0) % swing bus or PV bus
            Bprimeprime(i,:)=[];
            Bprimeprime(:,i)=[];
            Temp=diag(Bprimeprime);
            break;
        end
    end
end
clear Temp;

'Begin Decoupled Power Flow ...'
NumIter=1;
Kp = true;
Kq = true;

DeltaPoverV = [];


DeltaQoverV = [];


DeltaAngle = [];
DeltaV = [];

Kv=false;
Kangle=false;
IterNum = 1;
% either p mismatch or q mismatch meets the minimum requirement
while( (Kp == true) | (Kq == true) )

    % calculate P mismatch for PQ bus and PV bus, no swing bus
    MaxDeltaP = 0;
    for i=1:(NumPQ+NumPV)
        DeltaPoverV(i)=0;
    end
    iDeltaPoverV = 1;
    for i=1:size(Bus,1)
        if (Bus(i,8)==3) % swing bus
            continue;
        end
        
        for j=1:size(Bus,1)
            Gij = real(YBus(i,j));
            Bij = imag(YBus(i,j));
            DELTAij = Bus(i,10) - Bus(j,10);
            DeltaPoverV(iDeltaPoverV) = DeltaPoverV(iDeltaPoverV) + Bus(j,9) * (Gij*cos(DELTAij) + Bij*sin(DELTAij));
        end
        % get the max delta p
        DeltaP = DeltaPoverV(iDeltaPoverV)*Bus(i,9);
        if (abs(DeltaP)>MaxDeltaP)
            MaxDeltaP = abs(DeltaP);
        end
        % (P-P(x))/V
        DeltaPoverV(iDeltaPoverV) = (Bus(i,11) - DeltaPoverV(iDeltaPoverV)*Bus(i,9)) / Bus(i,9);
        
        iDeltaPoverV = iDeltaPoverV + 1;
    end
    if (MaxDeltaP<Tolerance)
        Kp == false;
    else
        Kp == true;
    end
    if(Kp == true)
        % calculate delta angle
        DeltaAngle = (-1)* eye(NumPQ+NumPV) * inv(Bprime) * (DeltaPoverV');
        if (max(abs(DeltaAngle))<Tolerance)
            Kangle = true;
        end
        % update Angle
        iDeltaAngle=1;
        for i=1:NumBus
            if (Bus(i,8)==3) % swing bus
                continue;
            end
            Bus(i,10) = Bus(i,10) + DeltaAngle(iDeltaAngle);
            iDeltaAngle = iDeltaAngle + 1;
        end
    end
    
    
    % calculate P mismatch for PQ bus, no swing bus and PV bus
    MaxDeltaQ = 0;
    for i=1:NumPQ
        DeltaQoverV(i)=0;
    end
    iDeltaQoverV = 1;
    for i=1:size(Bus,1)
        if (Bus(i,8)==3 | Bus(i,8)==2) % swing bus or PV bus
            continue;
        end
                
        for j=1:size(Bus,1)
            Gij = real(YBus(i,j));
            Bij = imag(YBus(i,j));
            DELTAij = Bus(i,10) - Bus(j,10);
            DeltaQoverV(iDeltaQoverV) = DeltaQoverV(iDeltaQoverV) + Bus(j,9) * (Gij*sin(DELTAij) - Bij*cos(DELTAij));
        end
        
         % get the max delta q
        DeltaQ = DeltaQoverV(iDeltaQoverV) * Bus(i,9);
        if (abs(DeltaQ)>MaxDeltaQ)
            MaxDeltaQ = abs(DeltaQ);
        end
        
        % (Q-Q(x))/V
        DeltaQoverV(iDeltaQoverV) = (Bus(i,12) - DeltaQoverV(iDeltaQoverV) * Bus(i,9)) / Bus(i,9);
        
        iDeltaQoverV = iDeltaQoverV + 1;
    end
    % evaluate the mismatch
    if (DeltaQ<Tolerance)
        Kq == false;
    else
        Kq == true;
    end
    if(Kq == true)
        % calculate delta V
        DeltaV = (-1)* eye(NumPQ) * inv(Bprimeprime) * (DeltaQoverV');
        if (max(abs(DeltaV))<Tolerance)
            Kv = true;
        end
        % update Angle
        iDeltaV=1;
        for i=1:NumBus
            if (Bus(i,8)==3 | Bus(i,8)==2) % swing bus
                continue;
            end
            Bus(i,9) = Bus(i,9) + DeltaV(iDeltaV);
            iDeltaV = iDeltaV + 1;
        end
    end
    
    if (Kv==true & Kangle==true)
        break;
    end

   if (IterNum>MaxNumIter)
       break;
   end
   
   IterNum = IterNum + 1;
end

fprintf('After %d iteration, power flow mismatch meets the tolerance.\n',IterNum)
fprintf('BusID, Vmag, Vangle\n',IterNum)
for i=1:NumBus
    fprintf('%d, %f, %f\n',Bus(i,1),Bus(i,9),Bus(i,10))
end