那位大神帮我用matlab写个运筹学单纯性算法的程序,最好带注释,急求!!!

2024-12-31 15:51:37
推荐回答(1个)
回答1:

下面是我多年前学习最优化课程时编写的单纯形法程序,100%原创。

单纯形算法的基本思想是,从多面体的某个顶点出发,移动到使得目标函数有所改进的相邻顶点;然后,从相邻顶点出发,移动到另一个更好的顶点,直至到达最优的顶点。从代数的观点看,也就是从一个基本可行解出发,求出使目标函数得到改进的相邻的基本可行解,循环往复,直至求得最优的基本可行解,或者判断最优解不可能存在。

从一个基本可行解得到另一个与之相邻的基本可行解的过程称为转轴过程(pivoting process,简称转轴),就是以某个非基变量代替另一个基变量。如何确定换入和换出变量(即所谓的转轴法则)是不同单纯形方法的主要差别。

我所编写程序的主要特色在于,不仅能够给出最优解,而且还能够给出迭代过程中各步的单纯形表(simplex tableau),这对于理解单纯形法的求解过程非常有帮助。

程序自带调用实例,使用了符号数学工具箱,所以显得有点繁琐,好好读一读对你会有帮助。

function [G, BasisSet] = simplex_table(A, b, c, idxB)

if ~nargin
A=[
3 3 1 0 0;
4 -4 0 1 0;
2 -1 0 0 1
];
b=[30 16 12]';
c=[-3 -1 0 0 0];
idxB = 3 : 5;
simplex_table(A, b, c, idxB)
return
end

if nargin == 2
G = A;
BasisSet = b;
m = size(G, 1) - 1;
n = size(G, 2) - 1;
else
m = size(A,1);
n = size(A,2);
BasisSet = idxB;
idxN = true(1, n);
idxN(idxB) = ~idxN(idxB);
idxB = ~idxN;
B = sym( A(:, BasisSet) );
cB = c(BasisSet);
w = cB * B^-1;
G = sym( zeros( size(A)+1 ) );
G(1:end-1, 1:end-1) = B^-1 * A;
G(1:end-1, end) = B^-1 * b;
G(end, 1:end-1) = w * A -c;
G(end, end) = w * b;
end
output_stars(n*16);
count = 0;
fprintf('\n initial table:\n');
output_table(G, BasisSet);

while count <= 100,
last_row = G(end, 1:end-1);
try
[d, in] = max( double( last_row ) );
catch
syms M
lim = limit(last_row, M, inf);
[d, in] = max( double(lim) );
if d ==inf
lim = limit(last_row / M, M, inf);
[d, in] = max( double(lim) );
end
end
if d <= 0,
disp(' ***** success');
x = sym(zeros(1, n));
x(BasisSet) = G(1:end-1, end)';
fprintf(' Optimal solution is: x = %s, ', symvec2str(x));
fprintf(' fmin = %s', char( G(end,end) ) );
break,
end

b_ = G(1:end-1, end);
yi = G(1:end-1, in);
if max( double(yi) ) <= 0, disp(' ***** No solution'); break, end
r = double(b_ ./ yi);
r( double(yi) <= 0 ) = NaN;
if 1
[yir, out] = min( r );
else
[yir, out] = min( r(end:-1:1) );
out = length(r) + 1 - out;
end

yir = yi(out);
count = count + 1;
fprintf(' The #%i iteration:\n', count);
fprintf(' Pivoting: x%i out, x%i in\n', BasisSet(out), in);
BasisSet(out) = in;

G(out, :) = G(out, :) / yir;
for i = 1 : m + 1
if i ~= out
G(i, :) = G(i, :) - G(out, :) * G(i, in);
end
end
G = simple(G);
output_table(G, BasisSet);
end
output_stars(n*16);

function output_table(G, BasisSet)
n = size(G, 2);
m = size(G, 1);
fprintf('\n ');
for j = 1 : n - 1
fprintf('\t%9s%i', 'x', j );
end
fprintf('\n');
for i = 1 : m
if i < m,
fprintf(' x%i', BasisSet(i));
else
fprintf(' ');
end
for j = 1 : n
fprintf('\t%10s', char( G(i,j) ) );
end
fprintf('\n');
end
fprintf('\n');

function output_stars(n)
fprintf('\n');
for i=1:n, fprintf('-'); end
fprintf('\n');

function s = symvec2str(x)
n=length(x);
s = '[ ';
for i=1:n-1
s = [s char(x(i)) ', ' ];
end
s = [s char(x(n)) ' ]'];