پخش بار نیوتن رافسون در متلب — از صفر تا صد

1%% POWER FLOW ANALYSIS 
2
3% THIS CODE IS ANALYSES THE COMPLEX POWER FLOW FOR A NETWORK
4
5%Files needed in directory are:
6% BusInputData.csv
7% createYbus.m -forms Ybus from busInputData
8% data2bus.m - stores bus data in new structure for calcs
9% findJacob.m - finds jacobian
10% sparseAdd.m - used in Ybus
11% calcPQ.m - calculations PQ injections using calculated node V and Th
12
13% CODE BY CHRIS MUTZEL
14% CHRIS.MUTZEL@GMAIL.COM
15
16% EE197 - POWER SYSTEMS ANALYSIS
17% PROFESSOR A. STANKOVIC
18% DEPEARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING
19% TUFTS UNIVERSITY, MEDFORD, MA
20
21% DECEMEBER 21TH, 2011
22
23
24clc
25clear
26format short  % keep output to 4 decimal places
27disp('NR Power Flow')
28disp('Code by Chris Mutzel')
29disp('chris.mutze@gmail.com')
30Ttotal = tic;
31
32%load Ybus
33[Ybus, data, N, lines] = data2bus(); 
34save Ybus;
35
36%Create new storage for bustype, P, Q, V, theta
37busData = zeros(N,5);           
38for n=1:N    %populate dataBus
39    busData(data(n,1),1) = data(n,3);                   % bus type
40    busData(data(n,1),2) = (data(n,8) - data(n,6))/100; % Pspec [MW]
41    busData(data(n,1),3) = (data(n,9) - data(n,7))/100; % Qspec [MVAR]
42    busData(data(n,1),4) = data(n,4); % V 
43end                                 
44disp(' '); disp('Specified Injections'); disp('    #         Pspec     Qspec');
45disp([(1:N)' busData(:,2) busData(:,3)])
46busData(1,4) = 1.0; %Slack voltage is one, angle zero
47
48numPV=0;
49for i = 1:N  %find number of each bus type
50   if (busData(i,1) == 2)
51      numPV = numPV + 1; %m
52   end
53end
54m = numPV; % N-m-1 are PQ buses
55
56%Intialize N-R Variables
57ep  = .00000001;                %stopping criteria for bus power mismatches
58STOP = 0;                   %set to one when all deltaPQ < ep
59 
60X = [zeros(N-1,1); 1*ones(N-m-1,1)];    %flat start
61deltaPQ = zeros(2*N-m-2,1);             %Mismatch vector
62j = 0;                                  %iteration count
63
64Traph = tic;                       %Start Timer
65while STOP == 0;
66    if j == 0;
67        disp('######################################################')
68        disp('Begin Newton Raphson')
69        disp('Flat Start->')
70        X
71    end
72    
73  
74    itNum = int2str(j + 1); 
75    disp('######################################################')
76    disp(' '); disp(strcat('Begin Iteration: ',itNum)); disp(' ')
77    
78    % Update V in dataBus for all PQ , use voltages from here for PQcalc
79    c = 1;
80    for n=1:N
81        if busData(n,1) == 0; % if the bus is a PQ bus
82           busData(n,4) = X(N+c-1); % store its V along with knowns for calc
83        c = c+1;
84        end
85    end
86    
87    %Update thetaBus for all PQ, use for PQcalc
88     c = 1;
89     for n=1:N
90         if n == 1;
91            busData(n,5) = 0;
92         end
93         if (busData(n,1) == 2 || busData(n,1) == 0);
94            busData(n,5) = X(c); 
95         c = c+1;
96         end
97         
98     end
99         
100    % ###################################################### 
101    %Calculate PQcalc based on V and theta
102    [Pcalc, Qcalc] = calcPQ(Ybus, busData, N);
103    %disp('    #      Pspec    Qspec')
104    %disp([(1:N)' busData(:,2) busData(:,3)])
105    %disp('    #      Pcalc     Qcalc')
106    %disp([(1:N)' Pcalc Qcalc])
107   
108    % ######################################################
109    %Calculate bus power mismatches
110    
111    %P
112    deltaPQ(1:N-1,1) = busData(2:N, 2) - Pcalc(2:N,1);
113    
114    %Q
115    bc = 0; %counter for index in delta PQ
116    for n=1:N
117        if (busData(n,1) == 0) % only for PQ buses
118            bc = bc + 1;
119            deltaPQ(N+bc-1,1) = busData(n,3) - Qcalc(n,1); 
120        end
121    end
122    %clear bc
123    %disp('Mismatch vector =')
124    %disp(deltaPQ)
125    
126    % ######################################################
127    % Check for convergence
128    error = norm(deltaPQ,2);
129    disp('Error, given by 2-norm of Mismatch vector')
130    disp(error)
131          
132    if (error < ep);
133        STOP = 1;
134        disp('Power Flow Calculation stopped')
135        disp('Error Threshold Reached')
136        disp(' ')
137    end
138    
139    % ######################################################
140    % New iteration if error over threshold
141    if STOP==0;
142        
143        %Find new Jacobian
144        %if j <= 2;
145            Jacob = findJacob(Ybus, Pcalc, Qcalc, busData, N, m);
146        %end
147            
148    % Solve the mismatch equations
149    % Invert the Jacobian
150        invJ = inv(Jacob);
151    % Solve for deltaX 
152        dX = invJ*deltaPQ;
153        
154    %Update Theta and V    
155        X = X + dX;
156    end
157    j = j+1;
158    %if (STOP == 0)
159    %    disp('    #         V       Theta')
160    %    disp([(1:N)' busData(:,4) busData(:,5)*360/(2*pi)])
161    %    disp(strcat('Finished Iteration:', itNum))
162    %end
163    disp(' ')
164    if j >= 10
165        STOP = 1; disp(' ');
166        disp('Power Flow Calculation stopped')
167        disp('Maximum Iterations (10) Reached'); disp(' ')
168    end
169      
170end   %NR iteration
171Traph = toc(Traph); %Stop timer
172
173% ######################################################
174% Newton Raphson Finished
175% Calculate line flows
176Zbus = inv(Ybus);
177lineFlows = [lines zeros(length(lines), 1)];
178for n = 1:length(lineFlows)
179    lineFlows(n,3) =  (exp(1i*angle(Zbus(lineFlows(n,1),lineFlows(n,2))))/abs(Zbus(lineFlows(n,1),lineFlows(n,2))))*...
180        (abs(busData(lineFlows(n,1),4))^2 - abs(busData(lineFlows(n,1),4))*abs(busData(lineFlows(n,2),4))*...
181         exp(1i*(busData(lineFlows(n,1),5) - busData(lineFlows(n,2),5))));      
182end
183
184
185
186Ttotal = toc(Ttotal);
187disp(['Final Solution after ', int2str(j), ' iterations is:'])
188disp('    #         V       Theta')
189disp([(1:N)' busData(:,4) busData(:,5)*360/(2*pi)])
190disp('Line Flows are:')
191disp('  Bus n      to      Bus m              Snm (per unit)')
192disp(lineFlows) 
193disp('Computation Time for NR-iteration was:')
194disp([num2str(Traph), ' seconds'])
195disp('Computation Time for program was:')
196disp([num2str(Ttotal), ' seconds'])

1%% Create Ybus using bus impedance data as input
2
3%% get the data
4clc
5clear
6
7load BusInputData.csv;
8%load testBusData.csv;
9
10data = BusInputData;
11%data = testBusData;
12[a,b] = size(data);
13%data has format: first n rows list the bus numbers in the problem
14% row n+1 holds -999 in first column, marking end input bus #s
15% after -999 is a NaN row, next m entries are impedance data 
16
17% n determines the size of Ybus
18
19%% Read in bus #ing, map to internal numbering
20
21%Find the first row containing -999
22BusIndex = 1;
23stop = 0;
24
25%Find the end of the bus numbering
26while stop == 0
27    if data(BusIndex,1) == -999;
28        stop = 1; 
29    end
30    if (stop == 0);
31        BusIndex = BusIndex + 1;
32    end
33    if BusIndex > a;    
34        stop = 1;
35    end
36end
37
38%Read in the bus numbers and create an internal mapping scheme 
39minBus = min(data(1:BusIndex-1,1));
40maxBus = max(data(1:BusIndex-1,1));
41mapS = minBus - 1; % External number = internal number + mapS
42
43
44
45%%
46clear Ybus;
47%Initialize a sparse matrix element
48Ybus = [maxBus - mapS, maxBus - mapS, 0];
49
50% Find stop point in the data
51stop = 0;
52ConnEnd = BusIndex+1;
53while stop == 0
54    if data(ConnEnd, 1) == -999;
55        stop = 1;
56    end
57    if (stop == 0);
58        ConnEnd = ConnEnd + 1;
59    end
60end
61
62%for each connection, update Ybus
63for n = BusIndex+1:ConnEnd-1
64    NodeA = data(n,1) - mapS; NodeB = data(n,2) - mapS;
65    % Self terms = diagonal elements
66    yab = 1/(data(n,6) + 1i*data(n,7));
67    Ybus = sparseAdd(Ybus, NodeA, NodeA, yab +.5i*data(n,8));
68    Ybus = sparseAdd(Ybus, NodeB, NodeB, yab +.5i*data(n,8));
69    % Off diagonal terms, only fill the upper diagonal
70    if NodeA > NodeB
71        Ybus = sparseAdd(Ybus, NodeB, NodeA, -yab);
72    else
73        Ybus = sparseAdd(Ybus, NodeA, NodeB, -yab);
74    end
75end
76    
77
78%% Construct the real matrix Ybus
79
80YbusRe = zeros(Ybus(end,1:2));
81
82for m=1:length(Ybus)-1;
83    YbusRe(Ybus(m,1),Ybus(m,2)) = Ybus(m,3);
84    if Ybus(m,1) ~= Ybus(m,2)
85        YbusRe(Ybus(m,2),Ybus(m,1)) = Ybus(m,3);
86    end
87end
88
89YbusRe

1function [Pcalc,Qcalc] = calcPQ(Ybus, busData, N)
2% Summary: Calculates P and Q injections at each node
3% 
4% Inputs:  
5%   scalar N = numNonSlack is # of non slack busses in the power flow
6%   X is column vector with 1st N - 1  entries for non Slack buses [theta]
7%                           then N - m -1  entries for PQ buses [voltage]
8%
9% Outputs:  
10%   PQ injections at all buses
11
12%Create storage for output
13Pcalc = zeros(N,1); %top row will refer to slack bus
14Qcalc = zeros(N,1);
15
16%P
17for i=1:N
18    for k=1:N
19        Pcalc(i,1) = Pcalc(i,1) + busData(i,4)*busData(k,4)*(...
20            real(Ybus(i,k))*cos(busData(i,5)-busData(k,5)) +...
21            imag(Ybus(i,k))*sin(busData(i,5)-busData(k,5)));
22    end
23end
24
25
26for i=1:N
27    for k=1:N
28        Qcalc(i,1) = Qcalc(i,1) + busData(i,4)*busData(k,4)*(...
29            real(Ybus(i,k))*sin(busData(i,5)-busData(k,5)) -...
30            imag(Ybus(i,k))*cos(busData(i,5)-busData(k,5)));
31    end
32end
33
34end

1function [YbusRe,data,numBus, lines] = data2bus()
2
3load BusInputData.csv;
4%load Test3bus.csv;
5
6data = BusInputData;
7%data = Test3bus;
8[a,b] = size(data);
9%data has format: first n rows list the bus numbers in the problem
10% row n+1 holds -999 in first column, marking end input bus #s
11% after -999 is a NaN row, next m entries are impedance data 
12% n determines the size of Ybus
13
14%Find the first row containing -999
15BusIndex = 1;
16stop = 0;
17
18%Find the end of the bus numbering
19while stop == 0
20    if data(BusIndex,1) == -999;
21        stop = 1; 
22    end
23    if (stop == 0);
24        BusIndex = BusIndex + 1;
25    end
26    if BusIndex > a;   % end of file 
27        stop = 1;
28    end
29end
30
31%Read in the bus numbers and create an internal mapping scheme 
32minBus = min(data(1:BusIndex-1,1));
33maxBus = max(data(1:BusIndex-1,1));
34mapS = minBus - 1; % External number = internal number + mapS
35numBus = maxBus - mapS;
36
37clear Ybus;
38%Initialize a sparse matrix element
39Ybus = [maxBus - mapS, maxBus - mapS, 0];
40
41
42% Find stop point in the data
43stop = 0;
44ConnEnd = BusIndex+1;
45while stop == 0
46    if data(ConnEnd, 1) == -999;
47        stop = 1;
48    end
49    if (stop == 0);
50        ConnEnd = ConnEnd + 1;
51    end
52end
53lines = zeros(ConnEnd - BusIndex -1 ,2);
54%for each connection, update Ybus, also store connections
55%lines = zeros(length(BusIndex+1:ConnEnd-1),2);
56for n = BusIndex+1:ConnEnd-1
57    NodeA = data(n,1) - mapS; NodeB = data(n,2) - mapS;
58    lines(n - BusIndex,1) = NodeA;
59    lines(n - BusIndex,2) = NodeB;
60    % Self terms = diagonal elements
61    yab = 1/(data(n,6) + 1i*data(n,7));
62    Ybus = sparseAdd(Ybus, NodeA, NodeA, yab +.5i*data(n,8));
63    Ybus = sparseAdd(Ybus, NodeB, NodeB, yab +.5i*data(n,8));
64    % Off diagonal terms, only fill the upper diagonal
65    if NodeA > NodeB
66        Ybus = sparseAdd(Ybus, NodeB, NodeA, -yab);
67    else
68        Ybus = sparseAdd(Ybus, NodeA, NodeB, -yab);
69    end
70end 
71
72%Construct the real matrix Ybus
73YbusRe = zeros(Ybus(end,1:2));
74
75for m=1:length(Ybus)-1;
76    YbusRe(Ybus(m,1),Ybus(m,2)) = Ybus(m,3);
77    if Ybus(m,1) ~= Ybus(m,2)
78        YbusRe(Ybus(m,2),Ybus(m,1)) = Ybus(m,3);
79    end
80end
81
82
83end

1function newSparseM = sparseAdd(M, row, col, v)
2%Function for manipulating sparse matrices
3% Input:
4%       M is a sparse matrix object 
5%       row and col are location in a real matrix M to which v should be added
6% Output: 
7%       new sparse matrix element O
8
9% M is ordered by row then by column, third column holds value for row and
10% col.  Last row has holds dimensions of the matrix object
11
12    
13% See if the element exists
14loc = 0;
15found = 0;       % Loc is filled with new position for an element if found  = 0
16c1_index = 1;   % Starting row to check M for existing element
17endloop = 0;    % 
18[a,b] = size(M); % store dimensions of M
19
20while (endloop ~= 1);    
21    if M(c1_index, 1) == row;
22        if M(c1_index, 2) == col;   % the element exists, stores its location
23            loc = c1_index;
24            found = 1;
25            endloop = 1;
26        end
27        if M(c1_index, 2) > col;    % element for same row but higher column   
28            loc = c1_index;
29            found = 0;            
30            endloop = 1;
31        end
32    end
33    
34    if M(c1_index, 1) > row;
35            loc = c1_index;
36            found = 0;
37            endloop = 1; 
38    end
39    c1_index = c1_index + 1;
40    if (c1_index > a)  %Check that we are not out of bounds of the array
41       found = 0;
42       endloop = 1;
43    end
44end
45
46% if it does not exist, create it
47% array has same form up to newly added element, then new element, then
48% remaining contents
49if found == 0;
50    newSparseM = [M(1:loc-1,:); [row, col, v]; M(loc:end,:)];
51end
52
53%If does exist, v is added to element value
54if found == 1;
55    M(loc,3) = M(loc,3) + v;
56    newSparseM = M;
57end
58
59%insert code here to delete row if value is now zero, not important for
60%Ybus since all values remain non-zero once created.
61end

1function [newJacob] = findJacob2(Ybus,Pcalc, Qcalc, busData, N, m)
2% Summary: Finds the Jacobian for the power flow problem given
3%
4% Inputs:
5%   scalar N = numNonSlack is # of non slack busses in the power flow
6%   X is column vector with 1st N entries theta for non Slack buses
7%                           then N entries V for non Slack buses
8%
9% Outputs:
10%   Jacob 
11
12%Create storage for output
13%square matrix with n = N-1 + N - m - 1
14Jacob = zeros(2*N,2*N);
15
16
17sizeJ = 2*N;
18
19for Jr=1:2
20    for Jc=1:2
21        
22        %J11
23        if (Jr == 1 && Jc == 1)
24            for SJr=1:N
25                
26                for SJc=1:N
27                   
28                    if (SJr == SJc) % diagonal elements 
29                        Jacob(SJr,SJr) = -Qcalc(SJc,1) - imag(Ybus(SJc,SJc))*abs(busData(SJc,4))^2;
30                                
31                    else            % off-diagonal elements
32                        Jacob(SJr, SJc) = abs(busData(SJr,4))*abs(busData(SJc,4))*(...
33                            real(Ybus(SJr,SJc))*sin(busData(SJr,5)-busData(SJc,5))-...
34                            imag(Ybus(SJr,SJc))*cos(busData(SJr,5)-busData(SJc,5)));
35                        
36                    end
37                end
38            end
39            
40        end
41        
42        %J12
43        if (Jr == 1 && Jc == 2)
44            for SJr=1:N
45                for SJc=1:N
46                    if SJr == SJc % diagonal elements
47                        Jacob(SJr, N + SJc) = Pcalc(SJc)/abs(busData(SJc,4)) +...
48                            real(Ybus(SJc,SJc))*abs(busData(SJc,4));
49                    else          % off-diagaonal elements
50                        Jacob(SJr, N + SJc) = abs(busData(SJr,4))*(...
51                            real(Ybus(SJr, SJc))*cos(busData(SJr,5)-busData(SJc, 5))+...
52                            imag(Ybus(SJr, SJc))*sin(busData(SJr,5)-busData(SJc, 5)));
53                    end
54                end
55            end
56        end
57        
58        
59        %J21
60        if (Jr == 2 && Jc == 1)
61            for SJr=1:N
62                for SJc=1:N
63                    if SJr == SJc % diagonal elements
64                        Jacob(N + SJr, SJc) = Pcalc(SJc) - real(Ybus(SJc, SJc))*abs(busData(SJc, 4))^2;
65                    else          % off-diagonal elements
66                        Jacob(N + SJr, SJc) = -abs(busData(SJr,4))*abs(busData(SJc,4))*(...
67                            real(Ybus(SJr, SJc))*cos(busData(SJr,5) - busData(SJc,5)) +... 
68                            imag(Ybus(SJr, SJc))*sin(busData(SJr,5) - busData(SJc,5)));
69                    end
70                end
71            end
72        end
73        
74        %J22
75        if (Jr == 2 && Jc == 2)
76            for SJr=1:N
77                for SJc=1:N
78                    if SJr == SJc % diagonal elements
79                        Jacob(N + SJr, N +  SJc) = Qcalc(SJc,1)/abs(busData(SJc,4)) - ...
80                            imag(Ybus(SJc, SJc))*abs(busData(SJc,4));
81                    else          %off-diagonal elements
82                        Jacob(N + SJr, N +  SJc) = abs(busData(SJr, 4))*(...
83                            real(Ybus(SJr, SJc))*sin(busData(SJr,5)-busData(SJc,5))-...
84                            imag(Ybus(SJr, SJc))*cos(busData(SJr,5)-busData(SJc,5)));
85                            
86                    end
87                end
88            end
89        end
90        
91    end
92end
93
94% ######################################################
95%Just found the Jacob for all buses, all only for PQ + PV, PV
96newJacob = zeros(2*N-m-2,2*N-m-2);
97
98%J11
99for c = 2:N
100    newJacob(1:N-1,c-1) = Jacob(2:N, c);
101end
102
103%J12
104a=0;
105for b = 1:N
106    if (busData(b,1) == 0)
107        a = a + 1;
108        newJacob(1:N-1, N-1+a) = Jacob(2:N,b+N);
109    end
110end
111
112%J21
113a=0;
114for b=1:N
115    if (busData(b,1) == 0)
116        a = a + 1;
117        newJacob(N-1+a, 1:N-1) = Jacob(b+N, 2:N);
118    end
119end
120
121%J22
122%nest for loops do same thing as above
123
124a=0;
125for b=1:N
126    if (busData(b,1) == 0)
127        a=a+1;
128        c=0;
129        for d=1:N
130            if (busData(d,1) == 0)
131                c=c+1;
132                newJacob(N-1+a, N-1+c) = Jacob(b+N,d+N);
133            end
134        end
135    end
136end
137end