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