lpeqslvr.m

电力系统计算软件

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

%   MATPOWER
%   $Id: LPeqslvr.m,v 1.9 2006/03/08 15:16:19 ray Exp $
%   by Deqiang (David) Gan, PSERC Cornell & Zhejiang University
%   Copyright (c) 1996-2006 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ 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, PC1, PC2, QC1MIN, QC1MAX, ...
    QC2MIN, QC2MAX, RAMP_AGC, RAMP_10, RAMP_30, RAMP_Q, APF] = idx_gen;
[PW_LINEAR, POLYNOMIAL, MODEL, STARTUP, SHUTDOWN, NCOST, COST] = idx_cost;

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

%% generator info
on = find(gen(:, GEN_STATUS) > 0);      %% which generators are on?

%% set up constants
nb = size(bus, 1);
nl = size(branch, 1);
npv     = length(pv);
npq     = length(pq);
ng      = length(on);                   %% number of generators that are turned on
pvpq    = [pv; pq];

%% 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;