accuzusa
下面是我多年前学习最优化课程时编写的单纯形法程序，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) returnendif 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;endoutput_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; endyir = 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);endoutput_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 );endfprintf('/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');endfprintf('/n');function output_stars(n)fprintf('/n');for i=1:n， fprintf('-'); endfprintf('/n');function s = symvec2str(x)n=length(x);s = '[ ';for i=1:n-1 s = [s char(x(i)) '， ' ];ends = [s char(x(n)) ' ]']; 
