CALC_V_PQ_SUM Solves the power flow using the power summation method. [V, Qpv, Sf, St, Sslack, iter, success] = calc_v_pq_sum(Vslack,nb,nl,f,Zb,Ybf,Ybt,,Yd,Sd,pv,Pg,Vg,tol,iter_max) Solves for bus voltages, generator reactive power, branch active and reactive power flows and slack bus active and reactive power. The input data consist of slack bus voltage, vector "from bus" indices, branch impedance and shunt admittance, vector of bus shunt admittances and load demand, as well as vectors with indicies of PV buses with their specified voltages and active powers. It is assumed that the branches are ordered using the principle of oriented ordering: indicies of sending nodes are smaller then the indicies of the receiving nodes. The branch index is equal to the index of their receiving node. Branch addmittances are added in Yd and treated as constant admittance bus loads. The applied method is Voltage correction power flow (VCPF) taken from: D. Rajicic, R. Ackovski and R. Taleski, "Voltage correction power flow," IEEE Transactions on Power Delivery, vol. 9, no. 2, pp. 1056-1062, Apr 1994. https://doi.org/10.1109/61.296308 See also RADIAL_PF.
0001 function [V, Qpv, Sf, St, Sslack, iter, success] = calc_v_pq_sum(Vslack,nb,nl,f,Zb,Ybf,Ybt,Yd,Sd,pv,Pg,Vg,mpopt) 0002 %CALC_V_PQ_SUM Solves the power flow using the power summation method. 0003 % 0004 % [V, Qpv, Sf, St, Sslack, iter, success] = calc_v_pq_sum(Vslack,nb,nl,f,Zb,Ybf,Ybt,,Yd,Sd,pv,Pg,Vg,tol,iter_max) 0005 % 0006 % Solves for bus voltages, generator reactive power, branch active and 0007 % reactive power flows and slack bus active and reactive power. The input 0008 % data consist of slack bus voltage, vector "from bus" indices, branch 0009 % impedance and shunt admittance, vector of bus shunt admittances and 0010 % load demand, as well as vectors with indicies of PV buses with their 0011 % specified voltages and active powers. It is assumed that the branches 0012 % are ordered using the principle of oriented ordering: indicies of 0013 % sending nodes are smaller then the indicies of the receiving nodes. The 0014 % branch index is equal to the index of their receiving node. Branch 0015 % addmittances are added in Yd and treated as constant admittance bus 0016 % loads. The applied method is Voltage correction power flow (VCPF) taken 0017 % from: 0018 % D. Rajicic, R. Ackovski and R. Taleski, "Voltage correction power flow," 0019 % IEEE Transactions on Power Delivery, vol. 9, no. 2, pp. 1056-1062, Apr 1994. 0020 % https://doi.org/10.1109/61.296308 0021 % 0022 % See also RADIAL_PF. 0023 0024 %% initialize 0025 tol = mpopt.pf.tol; 0026 iter_max = mpopt.pf.radial.max_it; 0027 vcorr = mpopt.pf.radial.vcorr == 1; 0028 Sd(pv) = Sd(pv) - Pg; 0029 V = Vslack * ones(nb,1); 0030 Vold = V; 0031 iter = 0; 0032 success = 0; 0033 % ZIP load model 0034 pw = mpopt.exp.sys_wide_zip_loads.pw; 0035 qw = mpopt.exp.sys_wide_zip_loads.qw; 0036 if isempty(pw) 0037 pw = [1 0 0]; 0038 end 0039 if isempty(qw) 0040 qw = pw; 0041 end 0042 Sdz = real(Sd) * pw(3) + 1j * imag(Sd) * qw(3); % constant impedance 0043 Sdi = real(Sd) * pw(2) + 1j * imag(Sd) * qw(2); % constant current 0044 Sdp = real(Sd) * pw(1) + 1j * imag(Sd) * qw(1); % constant power 0045 % Add artificial branch at the top of the branch list, so that the branch 0046 % index for other branches is equal to the index of their receiving node. 0047 f = [0; f]; 0048 Zb = [0; Zb]; 0049 nl = nl + 1; 0050 %% make Zpv matrix, for calculation of the PV generators reactive powers 0051 if ~isempty(pv) 0052 Zpv = make_zpv(pv,nb,nl,f,Zb,Yd); 0053 Bpv = (imag(Zpv))^-1; 0054 end 0055 npv = length(pv); 0056 Qpv = zeros(npv,1); 0057 %% do brackward-forward iterations 0058 if mpopt.verbose > 1 0059 fprintf('\n it max V mismatch (p.u.)'); 0060 fprintf('\n---- ----------------------'); 0061 end 0062 while success == 0 && iter < iter_max 0063 iter = iter + 1; 0064 % calculate load demand using actual voltages 0065 Vm = abs(V); 0066 S = Sdp + Sdi.*Vm + Sdz.*Vm.^2 + conj(Yd).*Vm.^2; 0067 % backward sweep 0068 St = S; 0069 Sf = St; 0070 for k = nl:-1:2 0071 i = f(k); 0072 Sf(k) = St(k) + Zb(k) * abs(St(k)/V(k))^2; 0073 St(i) = St(i) + Sf(k); 0074 end 0075 % forward sweep 0076 for k = 2:nl 0077 i = f(k); 0078 V(k) = V(i) - Zb(k) * conj(Sf(k)/V(i)); 0079 end 0080 % check for convergence 0081 DU = abs(V - Vold); 0082 DU(isnan(DU)) = inf; 0083 if mpopt.verbose > 1 0084 fprintf('\n%3d %10.3e', iter, max(DU)); 0085 end 0086 if max(DU) > tol 0087 Vold = V; 0088 % update PV generators reactive powers 0089 if ~isempty(pv) 0090 DE = (Vg./abs(V(pv))-1).*real(V(pv)); % Rajicic (VCPF) 0091 DD = Bpv * DE; 0092 DC = DD .* imag(V(pv))./real(V(pv)); 0093 if vcorr 0094 V_corr = make_vcorr(DC+1j*DD,pv,nb,nl,f,Zb); 0095 V = V + V_corr; 0096 end 0097 DQ = DD .* abs(V(pv)).^2 ./ real(V(pv)); 0098 Qpv = Qpv + DQ; 0099 Sdp(pv) = Sdp(pv) - 1j*DQ; 0100 end 0101 else 0102 success = 1; 0103 end 0104 end 0105 if mpopt.verbose 0106 if success 0107 fprintf('\nPower summation converged in %d iterations.\n', iter); 0108 else 0109 fprintf('\nPower summation did not converge in %d iterations.\n', iter); 0110 end 0111 end 0112 %% calculate branch flows 0113 % take out the first artificial branch 0114 Sslack = St(1); 0115 Sf = Sf(2:end); 0116 St = St(2:end); 0117 f = f(2:end); 0118 % correct branch flows to account for branch shunt admittances 0119 Sf = Sf + conj(Ybf) .* abs(V(f)).^2; 0120 St = St - conj(Ybt) .* abs(V(2:end)).^2;