龚纯《精通MATLAB最优化计算》随书源码(M文件)

function [intx,intf] = DividePlane(A,c,b,baseVector) sz = size(A); nVia = sz(2); n = sz(1); xx = 1:nVia; if length(baseVector) ~= n disp('基变量的个数要与约束矩阵的行数相等！'); mx = NaN; mf = NaN; return; end M = 0; sigma = -[transpose(c) zeros(1,(nVia-length(c)))]; xb = b; while 1 [maxs,ind] = max(sigma); if maxs <= 0 vr = find(c~=0 ,1,'last'); for l=1:vr ele = find(baseVector == l,1); if(isempty(ele)) mx(l) = 0; else mx(l)=xb(ele); end end if max(abs(round(mx) - mx))<1.0e-7 intx = mx; intf = mx*c; return; else sz = size(A); sr = sz(1); sc = sz(2); [max_x, index_x] = max(abs(round(mx) - mx)); [isB, num] = find(index_x == baseVector); fi = xb(num) - floor(xb(num)); for i=1:(index_x-1) Atmp(1,i) = A(num,i) - floor(A(num,i)); end for i=(index_x+1):sc Atmp(1,i) = A(num,i) - floor(A(num,i)); end Atmp(1,index_x) = 0; A = [A zeros(sr,1);-Atmp(1,:) 1]; xb = [xb;-fi]; baseVector = [baseVector sc+1]; sigma = [sigma 0]; %对偶单纯性 while 1 if xb >= 0 if max(abs(round(xb) - xb))<1.0e-7 vr = find(c~=0 ,1,'last'); for l=1:vr ele = find(baseVector == l,1); if(isempty(ele)) mx_1(l) = 0; else mx_1(l)=xb(ele); end end intx = mx_1; intf = mx_1*c; return; else sz = size(A); sr = sz(1); sc = sz(2); [max_x, index_x] = max(abs(round(mx_1) - mx_1)); [isB, num] = find(index_x == baseVector); fi = xb(num) - floor(xb(num)); for i=1:(index_x-1) Atmp(1,i) = A(num,i) - floor(A(num,i)); end for i=(index_x+1):sc Atmp(1,i) = A(num,i) - floor(A(num,i)); end Atmp(1,index_x) = 0; A = [A zeros(sr,1);-Atmp(1,:) 1]; xb = [xb;-fi]; baseVector = [baseVector sc+1]; sigma = [sigma 0]; continue; end else minb_1 = inf; chagB_1 = inf; sA = size(A); [br,idb] = min(xb); for j=1:sA(2) if A(idb,j)<0 bm = sigma(j)/A(idb,j); if bm<minb_1 minb_1 = bm; chagB_1 = j; end end end sigma = sigma -A(idb,:)*minb_1; xb(idb) = xb(idb)/A(idb,chagB_1); A(idb,:) = A(idb,:)/A(idb,chagB_1); for i =1:sA(1) if i ~= idb xb(i) = xb(i)-A(i,chagB_1)*xb(idb); A(i,:) = A(i,:) - A(i,chagB_1)*A(idb,:); end end baseVector(idb) = chagB_1; end end end else minb = inf; chagB = inf; for j=1:n if A(j,ind)>0 bz = xb(j)/A(j,ind); if bz<minb minb = bz; chagB = j; end end end sigma = sigma -A(chagB,:)*maxs/A(chagB,ind); xb(chagB) = xb(chagB)/A(chagB,ind); A(chagB,:) = A(chagB,:)/A(chagB,ind); for i =1:n if i ~= chagB xb(i) = xb(i)-A(i,ind)*xb(chagB); A(i,:) = A(i,:) - A(i,ind)*A(chagB,:); end end baseVector(chagB) = ind; end M = M + 1; if (M == 1000000) disp('找不到最优解！'); mx = NaN; minf = NaN; return; end end