gusucode.com > matpower工具箱源码程序 > matpower工具箱源码程序/MP2_0/LPeqslvr.m
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;