0001 function [MVAbase, bus, gen, branch, success, et] = runse(casedata, mpopt, fname, solvedcase)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041 [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ...
0042 VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus;
0043 [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ...
0044 TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ...
0045 ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch;
0046 [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ...
0047 MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ...
0048 QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
0049
0050
0051 if nargin < 4
0052 solvedcase = '';
0053 if nargin < 3
0054 fname = '';
0055 if nargin < 2
0056 mpopt = mpoption;
0057 if nargin < 1
0058 casedata = 'case9';
0059 end
0060 end
0061 end
0062 end
0063
0064
0065 dc = strcmp(upper(mpopt.model), 'DC');
0066
0067
0068 [baseMVA, bus, gen, branch] = loadcase(casedata);
0069 [i2e, bus, gen, branch] = ext2int(bus, gen, branch);
0070
0071
0072 [ref, pv, pq] = bustypes(bus, gen);
0073
0074
0075 on = find(gen(:, GEN_STATUS) > 0);
0076 gbus = gen(on, GEN_BUS);
0077
0078
0079 t0 = clock;
0080 if dc
0081
0082 Va0 = bus(:, VA) * (pi/180);
0083
0084
0085 [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch);
0086
0087
0088
0089 Pbus = real(makeSbus(baseMVA, bus, gen)) - Pbusinj - bus(:, GS) / baseMVA;
0090
0091
0092 Va = dcpf(B, Pbus, Va0, ref, pv, pq);
0093
0094
0095 branch(:, [QF, QT]) = zeros(size(branch, 1), 2);
0096 branch(:, PF) = (Bf * Va + Pfinj) * baseMVA;
0097 branch(:, PT) = -branch(:, PF);
0098 bus(:, VM) = ones(size(bus, 1), 1);
0099 bus(:, VA) = Va * (180/pi);
0100
0101
0102
0103 refgen = find(gbus == ref);
0104 gen(on(refgen(1)), PG) = gen(on(refgen(1)), PG) + (B(ref, :) * Va - Pbus(ref)) * baseMVA;
0105
0106 success = 1;
0107 else
0108
0109
0110 V0 = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA));
0111 V0(gbus) = gen(on, VG) ./ abs(V0(gbus)).* V0(gbus);
0112
0113
0114 [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch);
0115
0116
0117 Sbus = makeSbus(baseMVA, bus, gen);
0118
0119
0120 alg = upper(mpopt.pf.alg);
0121 switch alg
0122 case 'NR'
0123 [V, success, iterations] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, mpopt);
0124 case {'FDXB', 'FDBX'}
0125 [Bp, Bpp] = makeB(baseMVA, bus, branch, alg);
0126 [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, pv, pq, mpopt);
0127 case 'GS'
0128 [V, success, iterations] = gausspf(Ybus, Sbus, V0, ref, pv, pq, mpopt);
0129 otherwise
0130 error('Only Newton''s method, fast-decoupled, and Gauss-Seidel power flow algorithms currently implemented.');
0131 end
0132
0133
0134 [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq);
0135 end
0136 et = etime(clock, t0);
0137
0138
0139
0140 Pflf=branch(:,PF);
0141 Qflf=branch(:,QF);
0142 Ptlf=branch(:,PT);
0143 Qtlf=branch(:,QT);
0144 Sbuslf = V .* conj(Ybus * V);
0145 Vlf=V;
0146
0147
0148 [V, converged, i] = state_est(branch, Ybus, Yf, Yt, Sbuslf, Vlf, ref, pv, pq, mpopt);
0149
0150
0151
0152 Sbus = V .* conj(Ybus * V);
0153 bus(pq, PD) = -real(Sbus(pq)) * baseMVA;
0154 bus(pq, QD) = -imag(Sbus(pq)) * baseMVA;
0155
0156 on = find(gen(:, GEN_STATUS) > 0);
0157 gbus = gen(on, GEN_BUS);
0158 gen(on, PG) = real(Sbus(gbus)) * baseMVA + bus(gbus, PD);
0159
0160 [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq);
0161
0162
0163 Pfe=branch(:,PF);
0164 Qfe=branch(:,QF);
0165 Pte=branch(:,PT);
0166 Qte=branch(:,QT);
0167 nbr = length(Pfe);
0168 subplot(3,2,1), plot(180/pi*(angle(Vlf)-angle(V)),'.'), title('Voltage Angle (deg)');
0169 subplot(3,2,2), plot(abs(Vlf)-abs(V),'.'), title('Voltage Magnitude (p.u.)');
0170 subplot(3,2,3), plot((1:nbr),(Pfe-Pflf),'r.',(1:nbr),(Pte-Ptlf),'b.'), title('Real Flow (MW)');
0171 subplot(3,2,4), plot((1:nbr),(Qfe-Qflf),'r.',(1:nbr),(Qte-Qtlf),'b.'), title('Reactive Flow (MVAr)');
0172 subplot(3,2,5), plot(baseMVA*real(Sbuslf-Sbus), '.'), title('Real Injection (MW)');
0173 subplot(3,2,6), plot(baseMVA*imag(Sbuslf-Sbus), '.'), title('Reactive Injection (MVAr)');
0174
0175
0176
0177
0178 [bus, gen, branch] = int2ext(i2e, bus, gen, branch);
0179 if fname
0180 [fd, msg] = fopen(fname, 'at');
0181 if fd == -1
0182 error(msg);
0183 else
0184 if mpopt.out.all == 0
0185 printpf(baseMVA, bus, gen, branch, [], success, et, fd, ...
0186 mpoption(mpopt, 'out.all', -1));
0187 else
0188 printpf(baseMVA, bus, gen, branch, [], success, et, fd, mpopt);
0189 end
0190 fclose(fd);
0191 end
0192 end
0193 printpf(baseMVA, bus, gen, branch, [], success, et, 1, mpopt);
0194
0195
0196 if solvedcase
0197 savecase(solvedcase, baseMVA, bus, gen, branch);
0198 end
0199
0200
0201
0202 if nargout, MVAbase = baseMVA; end