lpsetup.m

电力系统计算软件

function [x, duals, idx_workc, idx_bindc] = LPsetup(a, f, b, nequs, vlb, vub, idx_workc, mpopt)
% LPSETUP 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
%   $Id: LPsetup.m,v 1.8 2004/11/02 15:30:59 ray Exp $
%   by Deqiang (David) Gan, PSERC Cornell & Zhejiang University
%   Copyright (c) 1996-2004 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ for more info.

%% options
alg = mpopt(11);

% ----- solve LP directly -----

if opf_slvr(alg) == 3           %% sparse LP with full constraints
    [x, duals] = mp_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;