CALC_V_I_SUM Solves the power flow using the current summation method. [V, Qpv, Sf, St, Sslack, iter, success] = calc_v_i_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 current summation taken from: D. Shirmohammadi, H. W. Hong, A. Semlyen and G. X. Luo, "A compensation-based power flow method for weakly meshed distribution and transmission networks," IEEE Transactions on Power Systems, vol. 3, no. 2, pp. 753-762, May 1988. https://doi.org/10.1109/59.192932 and G. X. Luo and A. Semlyen, "Efficient load flow for large weakly meshed networks," IEEE Transactions on Power Systems, vol. 5, no. 4, pp. 1309-1316, Nov 1990. https://doi.org/10.1109/59.99382 See also RADIAL_PF.
0001 function [V, Qpv, Sf, St, Sslack, iter, success] = calc_v_i_sum(Vslack,nb,nl,f,Zb,Ybf,Ybt,Yd,Sd,pv,Pg,Vg,mpopt) 0002 %CALC_V_I_SUM Solves the power flow using the current summation method. 0003 % 0004 % [V, Qpv, Sf, St, Sslack, iter, success] = calc_v_i_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 current summation taken from: 0017 % D. Shirmohammadi, H. W. Hong, A. Semlyen and G. X. Luo, 0018 % "A compensation-based power flow method for weakly meshed distribution 0019 % and transmission networks," IEEE Transactions on Power Systems, 0020 % vol. 3, no. 2, pp. 753-762, May 1988. https://doi.org/10.1109/59.192932 0021 % and 0022 % G. X. Luo and A. Semlyen, "Efficient load flow for large weakly meshed 0023 % networks," IEEE Transactions on Power Systems, vol. 5, no. 4, 0024 % pp. 1309-1316, Nov 1990. https://doi.org/10.1109/59.99382 0025 % 0026 % See also RADIAL_PF. 0027 0028 %% initialize 0029 tol = mpopt.pf.tol; 0030 iter_max = mpopt.pf.radial.max_it; 0031 vcorr = mpopt.pf.radial.vcorr == 1; 0032 Sd(pv) = Sd(pv) - Pg; 0033 V = Vslack * ones(nb,1); 0034 Vold = V; 0035 iter = 0; 0036 success = 0; 0037 % ZIP load model 0038 pw = mpopt.exp.sys_wide_zip_loads.pw; 0039 qw = mpopt.exp.sys_wide_zip_loads.qw; 0040 if isempty(pw) 0041 pw = [1 0 0]; 0042 end 0043 if isempty(qw) 0044 qw = pw; 0045 end 0046 Sdz = real(Sd) * pw(3) + 1j * imag(Sd) * qw(3); % constant impedance 0047 Sdi = real(Sd) * pw(2) + 1j * imag(Sd) * qw(2); % constant current 0048 Sdp = real(Sd) * pw(1) + 1j * imag(Sd) * qw(1); % constant power 0049 % Add artificial branch at the top of the branch list, so that the branch 0050 % index for other branches is equal to the index of their receiving node. 0051 f = [0; f]; 0052 Zb = [0; Zb]; 0053 nl = nl + 1; 0054 %% make Zpv matrix, for calculation of the PV generators reactive powers 0055 if ~isempty(pv) 0056 Zpv = make_zpv(pv,nb,nl,f,Zb,Yd); 0057 Bpv = (imag(Zpv))^-1; 0058 end 0059 npv = length(pv); 0060 Qpv = zeros(npv,1); 0061 %% do backward-forward iterations 0062 if mpopt.verbose > 1 0063 fprintf('\n it max V mismatch (p.u.)'); 0064 fprintf('\n---- ----------------------'); 0065 end 0066 while success == 0 && iter < iter_max 0067 iter = iter + 1; 0068 % calculate load demand using actual voltages 0069 Vm = abs(V); 0070 S = Sdp + Sdi.*Vm + Sdz.*Vm.^2 + conj(Yd).*Vm.^2; 0071 % backward sweep 0072 I = conj(S./V); 0073 for k = nl:-1:2 0074 i = f(k); 0075 I(i) = I(i) + I(k); 0076 end 0077 % forward sweep 0078 for k = 2:nl 0079 i = f(k); 0080 V(k) = V(i) - Zb(k) * I(k); 0081 end 0082 % check for convergence 0083 DU = abs(V - Vold); 0084 DU(isnan(DU)) = inf; 0085 if mpopt.verbose > 1 0086 fprintf('\n%3d %10.3e', iter, max(DU)); 0087 end 0088 if max(DU) > tol 0089 Vold = V; 0090 % update PV generators reactive powers 0091 if ~isempty(pv) 0092 DE = (Vg./abs(V(pv))-1).*real(V(pv)); % Rajicic (VCPF) 0093 DD = Bpv * DE; 0094 if vcorr 0095 DC = DD .* imag(V(pv))./real(V(pv)); 0096 V_corr = make_vcorr(DC+1j*DD,pv,nb,nl,f,Zb); 0097 V = V + V_corr; 0098 end 0099 DQ = DD .* abs(V(pv)).^2 ./ real(V(pv)); 0100 Qpv = Qpv + DQ; 0101 Sdp(pv) = Sdp(pv) - 1j*DQ; 0102 end 0103 else 0104 success = 1; 0105 end 0106 end 0107 if mpopt.verbose 0108 if success 0109 fprintf('\nCurrent summation converged in %d iterations.\n', iter); 0110 else 0111 fprintf('\nCurrent summation did not converge in %d iterations.\n', iter); 0112 end 0113 end 0114 %% calculate branch flows 0115 % take out the first artificial branch 0116 Sslack = V(1)*conj(I(1)); 0117 I = I(2:end); 0118 f = f(2:end); 0119 % calculate branch flows 0120 Sf = V(f).*conj(I) + conj(Ybf) .* abs(V(f)).^2; 0121 St = V(2:end).*conj(I) - conj(Ybt) .* abs(V(2:end)).^2;