CALC_V_PQ_SUM Solves the power flow using the admittance summation method. [V, Qpv, Sf, St, Sslack, iter, success] = calc_v_y_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 admittance summation taken from: Dragoslav Rajičić, Rubin Taleski, Two novel methods for radial and weakly meshed network analysis, Electric Power Systems Research, Volume 48, Issue 2, 15 December 1998, Pages 79-87 https://doi.org/10.1016/S0378-7796(98)00067-4 See also RADIAL_PF.
0001 function [V, Qpv, Sf, St, Sslack, iter, success] = calc_v_y_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 admittance summation method. 0003 % 0004 % [V, Qpv, Sf, St, Sslack, iter, success] = calc_v_y_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 admittance summation taken from: 0017 % Dragoslav Rajičić, Rubin Taleski, Two novel methods for radial and 0018 % weakly meshed network analysis, Electric Power Systems Research, 0019 % Volume 48, Issue 2, 15 December 1998, Pages 79-87 0020 % https://doi.org/10.1016/S0378-7796(98)00067-4 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 Ye = conj(Sdz) + Yd; 0063 D = zeros(nl,1); 0064 for k = nl:-1:2 0065 D(k) = 1 / (1 + Zb(k)*Ye(k)); 0066 i = f(k); 0067 Ye(i) = Ye(i) + D(k)*Ye(k); 0068 end 0069 while success == 0 && iter < iter_max 0070 iter = iter + 1; 0071 % calculate load demand using actual voltages (Sdz and Yd are already 0072 % included in Ye) 0073 Vm = abs(V); 0074 S = Sdp + Sdi.*Vm; 0075 % backward sweep 0076 Je = conj(S./V); 0077 for k = nl:-1:2 0078 i = f(k); 0079 Je(i) = Je(i) + D(k)*Je(k); 0080 end 0081 % forward sweep 0082 for k = 2:nl 0083 i = f(k); 0084 V(k) = D(k) * (V(i)-Zb(k)*Je(k)); 0085 end 0086 % check for convergence 0087 DU = abs(V - Vold); 0088 DU(isnan(DU)) = inf; 0089 if mpopt.verbose > 1 0090 fprintf('\n%3d %10.3e', iter, max(DU)); 0091 end 0092 if max(DU) > tol 0093 Vold = V; 0094 % update PV generators reactive powers 0095 if ~isempty(pv) 0096 DE = (Vg./abs(V(pv))-1).*real(V(pv)); % Rajicic (VCPF) 0097 DD = Bpv * DE; 0098 DC = DD .* imag(V(pv))./real(V(pv)); 0099 if vcorr 0100 V_corr = make_vcorr(DC+1j*DD,pv,nb,nl,f,Zb); 0101 V = V + V_corr; 0102 end 0103 DQ = DD .* abs(V(pv)).^2 ./ real(V(pv)); 0104 Qpv = Qpv + DQ; 0105 Sdp(pv) = Sdp(pv) - 1j*DQ; 0106 end 0107 else 0108 success = 1; 0109 end 0110 end 0111 if mpopt.verbose 0112 if success 0113 fprintf('\nAdmittance summation converged in %d iterations.\n', iter); 0114 else 0115 fprintf('\nAdmittance summation did not converge in %d iterations.\n', iter); 0116 end 0117 end 0118 %% calculate branch flows 0119 % take out the first artificial branch 0120 Sslack = V(1)*conj(Je(1)) + conj(Ye(1))*abs(V(1))^2; 0121 f = f(2:end); 0122 Zb = Zb(2:end); 0123 % calculate branch flows 0124 I = (V(f) - V(2:end))./Zb; 0125 Sf = V(f).*conj(I) + conj(Ybf) .* abs(V(f)).^2; 0126 St = V(2:end).*conj(I) - conj(Ybt) .* abs(V(2:end)).^2;