First, you will need to become familiar with basic Matlab programming, if you are not already.
The only difference is that in the last line of the script you will need to provide and index vector of the load buses (e.g. load_bus_idx) to index the rows of the bus matrix. sum_of_V_deviations = sum(abs(results.bus(load_bus_idx, VM) - 1.0)) — Ray > On Aug 31, 2018, at 9:03 AM, Sabhan Kanata <[email protected]> wrote: > > Dear Prof Ray Zimmerman > I have run this program. This program has been successful and result below: > > define_constants; > mpc = loadcase('case14'); > results = runpf(mpc); > total_losses = sum(get_losses(results)) > sum_of_V_deviations = sum(abs(results.bus(:, VM) - 1.0)) > > result... > total_losses = > 13.3933 +54.5383i > sum_of_V_deviations = > 0.6786 > > But after I calculate the amount of voltage deviation, the output of 0.6786 > is the voltage deviation for all buses. If I just want to calculate for the > load bus only. How does channge the script ? > > Thanks. > > > > > > Pada Jumat, 31 Agustus 2018 04.51.06 WIB, Ray Zimmerman <[email protected]> > menulis: > > > Sorry for the late response. Just catching up on e-mails after being away for > a while. > > To minimize power losses simply run an OPF with equal linear generator costs > for all generators, as described in FAQ #1 > <http://www.pserc.cornell.edu/matpower/#lossminimization>. To minimize > voltage deviations from some voltage set points, the best way is probably to > use an OPF with the generator costs set to zero and some extensions to the > standard OPF to put a cost on the voltage deviations. See Section 6.3 and > chapter 7 for details on extending the OPF. There is also a relevant example > starting at line 177 in t_opf_mips.m > <http://www.pserc.cornell.edu/matpower/docs/ref/matpower6.0/t/t_opf_mips.html>. > > Now, it is also possible to simply evaluate a power flow using a call to an > unmodified runpf(), followed by some additional code to compute the losses > (e.g. using get_losses()) or voltage deviations from the result. For example … > > define_constants; > mpc = loadcase('mycase'); > results = runpf(mpc); > total_losses = sum(get_losses(results)); > sum_of_V_deviations = sum(abs(results.bus(:, VM) - 1.0)); > > Hope this helps, > > Ray > > >> On Jul 28, 2018, at 12:03 AM, Sabhan Kanata <[email protected] >> <mailto:[email protected]>> wrote: >> >> Dear Mr. Ray Zimmerman, >> >> I am Sabhan Kanata, a Doctoral Student of Bandung Institute of Technology, >> Indonesia. I am taking a reseach on Optimazation of Power System to >> investigate the way to minimalize the power losses and voltage deviation. >> >> I have a big problem when dealing with the calculation of voltage deviation >> using MAT-POWER. Through this email, there is an attachment of runpf.m >> coding to calculate the losses. I need your assistance in oder to >> calculate the voltage deviation by employing the following formula.: >> >> min VD = min (Vi - 1.0) to minimize the deviations in voltage magnitud at >> load bus. >> >> I want to change objective function of losses to deviations in voltage. >> ====================================== >> function [PLOSS] = runpf(casename,AA) >> %RUNPF Runs a power flow. >> % MATPOWER >> % $Id: runpf.m,v 1.14 2006/09/29 19:23:07 ray Exp $ >> % by Ray Zimmerman, PSERC Cornell >> % Enforcing of generator Q limits inspired by contributions >> % from Mu Lin, Lincoln University, New Zealand (1/14/05). >> % Copyright (c) 1996-2005 by Power System Engineering Research Center >> (PSERC) >> % See http://www.pserc.cornell.edu/matpower/ >> <http://www.pserc.cornell.edu/matpower/> for more info. >> >> %%----- initialize ----- >> %% define named indices into bus, gen, branch matrices >> [PQ, PV, REF, NONE, BUS_I, BUS_TYPE, PD, QD, GS, BS, BUS_AREA, VM, ... >> VA, BASE_KV, ZONE, VMAX, VMIN, LAM_P, LAM_Q, MU_VMAX, MU_VMIN] = idx_bus; >> [F_BUS, T_BUS, BR_R, BR_X, BR_B, RATE_A, RATE_B, RATE_C, ... >> TAP, SHIFT, BR_STATUS, PF, QF, PT, QT, MU_SF, MU_ST, ... >> ANGMIN, ANGMAX, MU_ANGMIN, MU_ANGMAX] = idx_brch; >> [GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, GEN_STATUS, PMAX, PMIN, ... >> MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN, PC1, PC2, QC1MIN, QC1MAX, ... >> QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen; >> >> %% default arguments >> solvedcase = ''; %% don't save solved case >> fname = ''; %% don't print results to a file >> mpopt = mpoption; %% use default options >> >> >> %% options >> verbose = mpopt(31); >> qlim = mpopt(6); %% enforce Q limits on gens? >> dc = mpopt(10); %% use DC formulation? >> >> %% read data & convert to internal bus numbering >> [baseMVA, bus, gen, branch] = loadcase(casename); >> >> %Untuk sistem tenaga ieee 14 bus >> gen(1,6)=AA(1); >> gen(2,6)=AA(2); >> gen(3,6)=AA(3); >> gen(4,6)=AA(4); >> gen(5,6)=AA(5); >> bus(9,6)=2*AA(6); >> bus(14,6)=2*AA(7); >> branch(8,9)=1+0.025*AA(8); >> branch(9,9)=1+0.025*AA(9); >> branch(10,9)=1+0.025*AA(10); >> >> [i2e, bus, gen, branch] = ext2int(bus, gen, branch); >> >> %% get bus index lists of each type of bus >> [ref, pv, pq] = bustypes(bus, gen); >> >> %% generator info >> on = find(gen(:, GEN_STATUS) > 0); %% which generators are on? >> gbus = gen(on, GEN_BUS); %% what buses are they at? >> >> %%----- run the power flow ----- >> t0 = clock; >> if dc %% DC formulation >> %% initial state >> Va0 = bus(:, VA) * (pi/180); >> >> %% build B matrices and phase shift injections >> [B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA, bus, branch); >> >> %% compute complex bus power injections (generation - load) >> %% adjusted for phase shifters and real shunts >> Pbus = real(makeSbus(baseMVA, bus, gen)) - Pbusinj - bus(:, GS) / >> baseMVA; >> >> %% "run" the power flow >> Va = dcpf(B, Pbus, Va0, ref, pv, pq); >> >> %% update data matrices with solution >> branch(:, [QF, QT]) = zeros(size(branch, 1), 2); >> branch(:, PF) = (Bf * Va + Pfinj) * baseMVA; >> branch(:, PT) = -branch(:, PF); >> bus(:, VM) = ones(size(bus, 1), 1); >> bus(:, VA) = Va * (180/pi); >> %% update Pg for swing generator (note: other gens at ref bus are >> accounted for in Pbus) >> %% Pg = Pinj + Pload + Gs >> %% newPg = oldPg + newPinj - oldPinj >> refgen = find(gbus == ref); %% which is(are) the reference >> gen(s)? >> gen(on(refgen(1)), PG) = gen(on(refgen(1)), PG) + (B(ref, :) * Va - >> Pbus(ref)) * baseMVA; >> >> success = 1; >> else %% AC formulation >> %% initial state >> % V0 = ones(size(bus, 1), 1); %% flat start >> V0 = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA)); >> V0(gbus) = gen(on, VG) ./ abs(V0(gbus)).* V0(gbus); >> >> if qlim >> ref0 = ref; %% save index and angle of >> Varef0 = bus(ref0, VA); %% original reference bus >> limited = []; %% list of indices of gens @ Q >> lims >> fixedQg = zeros(size(gen, 1), 1); %% Qg of gens at Q limits >> end >> repeat = 1; >> while (repeat) >> %% build admittance matrices >> [Ybus, Yf, Yt] = makeYbus(baseMVA, bus, branch); >> >> %% compute complex bus power injections (generation - load) >> Sbus = makeSbus(baseMVA, bus, gen); >> >> %% run the power flow >> alg = mpopt(1); >> if alg == 1 >> [V, success, iterations] = newtonpf(Ybus, Sbus, V0, ref, pv, pq, >> mpopt); >> elseif alg == 2 | alg == 3 >> [Bp, Bpp] = makeB(baseMVA, bus, branch, alg); >> [V, success, iterations] = fdpf(Ybus, Sbus, V0, Bp, Bpp, ref, >> pv, pq, mpopt); >> elseif alg == 4 >> [V, success, iterations] = gausspf(Ybus, Sbus, V0, ref, pv, pq, >> mpopt); >> else >> error('Only Newton''s method, fast-decoupled, and Gauss-Seidel >> power flow algorithms currently implemented.'); >> end >> >> %% update data matrices with solution >> [bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, >> V, ref, pv, pq); >> >> if qlim %% enforce generator Q limits >> %% find gens with violated Q constraints >> mx = find( gen(:, GEN_STATUS) > 0 & gen(:, QG) > gen(:, QMAX) ); >> mn = find( gen(:, GEN_STATUS) > 0 & gen(:, QG) < gen(:, QMIN) ); >> >> if ~isempty(mx) | ~isempty(mn) %% we have some Q limit >> violations >> if isempty(pv) >> if verbose >> if ~isempty(mx) >> fprintf('Gen %d (only one left) exceeds upper Q >> limit : INFEASIBLE PROBLEM\n', mx); >> else >> fprintf('Gen %d (only one left) exceeds lower Q >> limit : INFEASIBLE PROBLEM\n', mn); >> end >> end >> success = 0; >> break; >> end >> if verbose & ~isempty(mx) >> fprintf('Gen %d at upper Q limit, converting to PQ >> bus\n', mx); >> end >> if verbose & ~isempty(mn) >> fprintf('Gen %d at lower Q limit, converting to PQ >> bus\n', mn); >> end >> >> %% save corresponding limit values >> fixedQg(mx) = gen(mx, QMAX); >> fixedQg(mn) = gen(mn, QMIN); >> mx = [mx;mn]; >> >> %% convert to PQ bus >> gen(mx, QG) = fixedQg(mx); %% set Qg to binding limit >> gen(mx, GEN_STATUS) = 0; %% temporarily turn off gen, >> for i = 1:length(mx) %% (one at a time, since >> bi = gen(mx(i), GEN_BUS); %% they may be at same bus) >> bus(bi, [PD,QD]) = ... %% adjust load accordingly, >> bus(bi, [PD,QD]) - gen(mx(i), [PG,QG]); >> end >> bus(gen(mx, GEN_BUS), BUS_TYPE) = PQ; %% & set bus type to >> PQ >> >> %% update bus index lists of each type of bus >> ref_temp = ref; >> [ref, pv, pq] = bustypes(bus, gen); >> if verbose & ref ~= ref_temp >> fprintf('Bus %d is new slack bus\n', ref); >> end >> limited = [limited; mx]; >> else >> repeat = 0; %% no more generator Q limits violated >> end >> else >> repeat = 0; %% don't enforce generator Q limits, once is >> enough >> end >> end >> if qlim & ~isempty(limited) >> %% restore injections from limited gens (those at Q limits) >> gen(limited, QG) = fixedQg(limited); %% restore Qg value, >> for i = 1:length(limited) %% (one at a time, since >> bi = gen(limited(i), GEN_BUS); %% they may be at same bus) >> bus(bi, [PD,QD]) = ... %% re-adjust load, >> bus(bi, [PD,QD]) + gen(limited(i), [PG,QG]); >> end >> gen(limited, GEN_STATUS) = 1; %% and turn gen back on >> if ref ~= ref0 >> %% adjust voltage angles to make original ref bus correct >> bus(:, VA) = bus(:, VA) - bus(ref0, VA) + Varef0; >> end >> end >> end >> et = etime(clock, t0); >> >> %%----- output results ----- >> %% convert back to original bus numbering & print results >> [bus, gen, branch] = int2ext(i2e, bus, gen, branch); >> if fname >> [fd, msg] = fopen(fname, 'at'); >> if fd == -1 >> error(msg); >> else >> printpf(baseMVA, bus, gen, branch, [], success, et, fd, mpopt); >> fclose(fd); >> end >> end >> [PLOSS]=printpf(baseMVA, bus, gen, branch, [], success, et, 1, mpopt); >> >> % % save solved case >> % if solvedcase >> % savecase(solvedcase, baseMVA, bus, gen, branch); >> % end >> >> %% this is just to prevent it from printing baseMVA >> %% when called with no output arguments >> if nargout, MVAbase = baseMVA; end >> >> return; >> >> ============================================================== >> >> >> Hopefully I can get your answer to cope with the problem i am now facing. >> I would like to thank you in advance for the help at your convenience. >> >> >> Yours Faithfully, >> >> >> Sabhan Kanata >> Thanks for your helping >> >> >> >
