www.gusucode.com > matpower工具箱源码程序 > matpower工具箱源码程序/MP2_0/LPsetup.m
function [x, duals, idx_workc, idx_bindc] = LPsetup(a, f, b, nequs, vlb, vub, idx_workc, mpopt) % LPSOLVER solves a LP problem using a callable LP routine % The LP problem is defined as follows: % % min f' * x % S.T. a * x =< b % vlb =< x =< vub % % All of the equality constraints must appear before inequality constraints. % nequs specifies how many of the constraints are equality constraints. % % The algorithm (set in mpopt) can be set to the following options: % % 220 - solve LP using ICS (equality constraints are eliminated) % 240 - solve LP using Iterative Constraint Search (ICS) % (equality constraints are preserved, typically superior to 220 and 250) % 250 - solve LP with full set of constraints % MATPOWER Version 2.0 % by Deqiang (David) Gan, PSERC Cornell 12/12/97 % Copyright (c) 1996, 1997 by Power System Engineering Research Center (PSERC) % See http://www.pserc.cornell.edu/ for more info. %% options alg = mpopt(11); % ----- solve LP directly ----- if opf_slvr(alg) == 3 %% sparse LP with full constraints a = full(a); [x, duals] = lp(f, a, b, vlb, vub, [], nequs, -1); duals = duals(1:length(b)); % MATLAB built-in LP solver has more elements in duals than we want idx_workc = []; idx_bindc = []; return; end % ----- solve LP using constraint relaxation (equality constraints are preserved) ------ if opf_slvr(alg) == 2 %% sparse LP with relaxed constraints if isempty(idx_workc) == 1 idx_workc = find(b < 1.0e-5); end [x, duals, idx_workc, idx_bindc] = LPrelax(a, f, b, nequs, vlb, vub, idx_workc, mpopt); return; end % ----- solve LP using constraint relaxation (equality constraints are eliminated) ------ % so opf_slvr(alg) == 1 %% dense LP % set up the indicies of variables and constraints idx_x1 = 1:nequs-1; idx_x2 = nequs : length(f); idx_c1 = 1:nequs-1; idx_c2 = nequs : length(b); % eliminate equality constraints b1 = b(idx_c1); b2 = b(idx_c2); a11 = a(idx_c1, idx_x1); a12 = a(idx_c1, idx_x2); a21 = a(idx_c2, idx_x1); a22 = a(idx_c2, idx_x2); a11b1 = a11 \ b1; a11a12 = a11 \ a12; % set up the reduced LP fred = -((f(idx_x1))' * a11a12)' + f(idx_x2); ared = [-a21 * a11a12 + a22 -a11a12 a11a12]; bred = [ b2 - a21 * a11b1 vub(idx_x1) - a11b1 a11b1 - vlb(idx_x1)]; vubred = vub(idx_x2); vlbred = vlb(idx_x2); nequsred = nequs - length(idx_x1); % solve the reduced LP problem using constraint relaxation if isempty(idx_workc) == 1 idx_workc = find(b2< 1.0e-5); end [x2, dualsred, idx_workc, idx_bindc] = LPrelax(ared, fred, bred, nequsred, vlbred, vubred, idx_workc, mpopt); % parse the solution of the reduced LP to get the solution of the original LP x(idx_x1) = a11b1 - a11a12 * x2; x(idx_x2) = x2; x = x'; dualsc2 = dualsred(1:length(idx_c2)); temp = find(dualsc2); dualsc1 = a11' \ ( -f(idx_x1) - (a21(temp, :))' * dualsc2(temp) ); duals(idx_c1) = dualsc1; duals(idx_c2) = dualsc2; duals = duals'; return;