lpsetup.m

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

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;