lpeqslvr.m

关于电力系统的计算程序!最新的改进法计算潮流程序！用matlab实现的！ 适用于电力行业的相关从业人员！

function [x, success] = LPeqslvr(x, baseMVA, bus, gen, gencost, branch, Ybus, Yf, Yt, V, ref, pv, pq, mpopt);

%   MATPOWER Version 2.0
%   by Deqiang (David) Gan, PSERC Cornell    12/16/97
%   Copyright (c) 1996, 1997 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/ for more info.


%% define named indices into data 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;
[GEN_BUS, PG, QG, QMAX, QMIN, VG, MBASE, ...
        GEN_STATUS, PMAX, PMIN, MU_PMAX, MU_PMIN, MU_QMAX, MU_QMIN] = idx_gen;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, N, COST] = idx_cost;

%% options
alg		= mpopt(11);			%% algorithm
verbose	= mpopt(31);			%% verbose

%% set up constants
nb = size(bus, 1);
nl = size(branch, 1);
npv     = length(pv);
npq     = length(pq);
ng = npv + 1;					%% number of generators that are turned on
on = find(gen(:, GEN_STATUS));	%% which generators are on?
pvpq = [pv; pq];
gbus = gen(on, GEN_BUS);		%% what buses are they at?  

%% check for costs for Qg
[pcost, qcost] = pqcost(gencost, size(gen, 1), on);

%parse x, update bus, gen
bus(pvpq, VA) = x(1:nb-1) * 180/pi;
bus(:, VM)= x(nb:nb + nb-1);
gen(on, PG) = x(2*nb : 2*nb-1+ng) * baseMVA;
gen(on, QG) = x(2*nb+ng : 2*nb-1+2*ng) * baseMVA;


% run load flow
V       = bus(:, VM) .* exp(sqrt(-1) * pi/180 * bus(:, VA));

success = 1;
Sbus = makeSbus(baseMVA, bus, gen);											%% build Sg-Sl

%% turn down verbosity one level for call to power flow
if verbose
	mpopt = mpoption(mpopt, 'VERBOSE', verbose-1);
end
[V, converged, iterations] = newtonpf(Ybus, Sbus, V, ref, pv, pq, mpopt);	%% do NR iteration
if converged == 0, success = 0; end
[bus, gen, branch] = pfsoln(baseMVA, bus, gen, branch, Ybus, Yf, Yt, V, ref, pv, pq);	%% post-processing
et = 1;
% printpf(baseMVA, bus, gen, branch, [], success, et, 1, mpopt);

Pg = gen(on, PG)/baseMVA;
Qg = gen(on, QG) / baseMVA;

%% set up x
if opf_form(alg) == 1
	Cp = [];
	Cq = [];
else
	Cp = totcost(pcost, Pg * baseMVA);
	Cq = totcost(qcost, Qg * baseMVA);		%% empty if qcost is empty
end
x = [angle(V(pvpq)); abs(V); Pg; Qg; Cp; Cq];

return;