//解线性规划的用乘数的单纯形法
//问题:x0 + Sum{j = 1..n;
a[0, j] * x[j]} = a[0, 0] (1)
// Sum{a[i, j] * x[j] = a[i, 0] (i = 1..m) (2)
// x[j] >= 0 (j = 1..n) (3)
//即求出满足约束条件(2)和(3)、且使目标函数x0为极小或极大的
//x = (x[j1], x[j2], ..., x[jm], 0, ..., 0)。
//其中，a[i, j] (i = 0..m;
j = 0..n)开始时给定，在求解过程中
//被不断地改变。
//
//参数：
//d 指示值。当d = -1时解极大问题，d = 1时解极小问题
//m 线性规划问题的阶数（即约束条件数）
//n 线性规划问题的维数（为原给的线性规划问题的变量数加m
// 详见例子
//id 数组id[0..m]，存放基本变量的下标，开始时id[0] = 0,
// id = n - m + i (i = 1..m)
//a 数组a[0..m, 0..n]，开始时存放线性规划问题的全部给定
// 数据，在求解过程中其值不断改变（具体存放方式见例子）
//k 指示量，计算结束时如果k = 0则表示线性线性规划问题有
// 有限最优解。
//常数说明：
//PSUMChangeBaseZero 规定老基变新基运算中的无穷小量
//PSUMAccumulateZero 为减小由于舍入误差的积累造成错误允许值，
// 用于引进新基本变量
//PSUMZero 在决定应从原来的基本变量中去掉的基本变量的
// 下标时，为减小舍入误差的影响而使用的无穷小。
//例子：
// x0 -1.1 x1 - 2.2 x2 + 3.3 x3 - 4.4 x4 = 0
// x1 + x2 + 2 x3 = 4
// x1 + 2 x2 + 2.5 x3 + 3 x4 = 5
//d = -1, m = 2, 变量数 = 4，于是：
//m = 2, n = 4 + m = 6,
//a = 0, -1.1, -2.2, 3.3, -4.4, (0), (0) //后面补m个0
// 4, 1, 1, 2, 0, (1), (0) //后面形成单位阵
// 5, 1, 2, 2.5, 3, (0), (1)
type
TMatrix = array of array of Extended;
const
PSUMInfinity: Extended = 1E30;
PSUMChangeBaseZero: Extended = 1E-15;
PSUMAccumulateZero: Extended = 1E-15;
PSUMZero: Extended = 1E-15;
procedure PSUM(d: Integer;
m, n: Integer;
var id: array of Integer;
a: TMatrix;
var k: Integer);
implementation
procedure PSUM(d: Integer;
m, n: Integer;
var id: array of Integer;
a: TMatrix;
var k: Integer);
var
i, j, L: Integer;
f, exm: Extended;
function SI(c, cr: Extended): Extended;
begin
if Abs(c) > cr then
Result := c else
Result := 0;
end;
begin
repeat
exm := PSUMInfinity * d;
k := -1;
for j := 1 to ndo
if (d * a[0, j] > PSUMAccumulateZero)
and (-d * (a[0, j] - exm) > 0) then
begin
exm := a[0, j];
k := j;
end;
if k = -1 then
begin
k := 0;
Break;
end;
if k > n - m then
Break;
exm := PSUMInfinity;
L := -1;
for i := 1 to mdo
begin
f := a[i, 0] / a[i, k];
if (a[i, k] > PSUMChangeBaseZero)
and (exm > f) then
begin
exm := f;
L := i;
end;
end;
if L = -1 then
Break;
id[L] := k;
for j := 0 to ndo
if j <> k then
a[L, j] := a[L, j] / a[L, k];
a[L, k] := 1;
for i := 0 to mdo
if i <> L then
begin
for j := 0 to ndo
if j <> k then
a[i, j] := SI(a[i, j] - a[L, j] * a[i, k], PSUMChangeBaseZero);
a[i, k] := SI(a[i, k] * (1 - a[L, k]), PSUMChangeBaseZero);
end;
until False;
end;
例子：
var
a: TMatrix;
id: array [0..2] of Integer;
k: Integer;
begin
SetLength(a, 3, 7);
a[0, 0] := 0;
a[0, 1] := -1.1;
a[0, 2] := -2.2;
a[0, 3] := 3.3;
a[0, 4] := -4.4;
a[0, 5] := 0;
a[0, 6] := 0;
a[1, 0] := 4;
a[1, 1] := 1;
a[1, 2] := 1;
a[1, 3] := 2;
a[1, 4] := 0;
a[1, 5] := 1;
a[1, 6] := 0;
a[2, 0] := 5;
a[2, 1] := 1;
a[2, 2] := 2;
a[2, 3] := 2.5;
a[2, 4] := 3;
a[2, 5] := 0;
a[2, 6] := 1;
id[0] := 0;
id[1] := 4;
id[2] := 5;
PSUM(-1, 2, 6, id, a, k);
if k = 0 then
ShowMessage('failed')
else
ShowMessage(Format('Base variables: x%d = %f, x%d = %f, max x0 = %f',
[id[1], a[1, 0], id[2], a[2, 0], a[0, 0]]));
end;
注：如果要改变上界或下界，那么要添加变量数来调整矩阵A
如：下界：x2 > 3
则要添加一个变量和一个约束：
x2 - x[k] = 3
上界：x2 < 10
则要添加一个变量和一个约束：
x2 + x[k] = 10
其中k为原变量数 + 1
