Dear Prof Ray Zimmerman
I have run thisprogram. 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. 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.
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]> 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 rayExp $
% by Ray Zimmerman, PSERC Cornell
% Enforcing of generator Q limits inspired bycontributions
% from Mu Lin, Lincoln University, New Zealand(1/14/05).
% Copyright (c) 1996-2005 by Power SystemEngineering Research Center (PSERC)
% See http://www.pserc.cornell.edu/matpower/for more info.
%%----- initialize -----
%% define named indicesinto 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 afile
mpopt = mpoption; %% use default options
%% options
verbose = mpopt(31);
qlim = mpopt(6); %% enforceQ limits on gens?
dc = mpopt(10); %% use DCformulation?
%% read data & convertto internal bus numbering
[baseMVA, bus, gen, branch] =loadcase(casename);
%Untuk sistem tenaga ieee14 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 ofeach 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 %% DCformulation
%% initial state
Va0 = bus(:, VA) * (pi/180);
%% build B matrices and phase shiftinjections
[B, Bf, Pbusinj, Pfinj] = makeBdc(baseMVA,bus, branch);
%% compute complex bus power injections(generation - load)
%% adjusted for phase shifters and realshunts
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) thereference
gen(s)?
gen(on(refgen(1)), PG) = gen(on(refgen(1)),PG) + (B(ref, :) * Va -
Pbus(ref)) * baseMVA;
success = 1;
else %% ACformulation
%%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 ofindices 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('OnlyNewton''s method, fast-decoupled, and Gauss-Seidel power
flow algorithmscurrently 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 Qlimits
%% find gens with violated Qconstraints
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 (onlyone left) exceeds upper Q
limit : INFEASIBLE PROBLEM\n', mx);
else
fprintf('Gen %d (onlyone left) exceeds lower Q
limit : INFEASIBLE PROBLEM\n', mn);
end
end
success = 0;
break;
end
if verbose & ~isempty(mx)
fprintf('Gen %d atupper Q limit, converting to PQ bus\n',
mx);
end
if verbose & ~isempty(mn)
fprintf('Gen %d atlower Q limit, converting to PQ bus\n',
mn);
end
%% save corresponding limitvalues
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 ofeach type of bus
ref_temp = ref;
[ref, pv, pq] = bustypes(bus,gen);
if verbose & ref ~= ref_temp
fprintf('Bus %d is newslack bus\n', ref);
end
limited = [limited; mx];
else
repeat = 0; %% no moregenerator Q limits violated
end
else
repeat = 0; %% don't enforce generator Qlimits, once is enough
end
end
if qlim & ~isempty(limited)
%% restore injections from limited gens(those at Q limits)
gen(limited, QG) =fixedQg(limited); %% restoreQg 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 makeoriginal ref bus correct
bus(:, VA) = bus(:, VA) - bus(ref0,VA) + Varef0;
end
end
end
et = etime(clock, t0);
%%----- output results -----
%% convert back to originalbus 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 preventit from printing baseMVA
%% when called with nooutput 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
