两阶段单纯形法的Matlab实现

function x=Simplex2phase(A,b,c) % 两阶段单纯形法求线性规划最优解，针对于A中未必刚好有一个m阶单位矩阵 % min cx s.t Ax=b,x>=0,注意A是m*n矩阵，b是m维列向量! c是n维行向量 % 第一阶段求解最优解 [m,n]=size(A); % 列单纯形表 subscript=((n+1):(n+m))'; % 基变量的下标 cB=ones(1,m); % 基变量对应的cB zj_cj=cB*A; % 判别式 f=cB*b; % 当前的函数值 whole=[subscript A eye(m) b;0 zj_cj zeros(1,m) f]; % whole和下面的whole2是单纯形表 % Whole 是(m+1)*(m+n+2)的矩阵 % Whole(m+1,2:(m+n+1)) 是判别数向量 while any(whole(m+1,2:(m+n+1))>0.00001)% m+1行，2到m+n+1列中的判别数>0 % >0.00001 是 >0， 如果判别式中有>0的，表明把这个>0的分量作为进基变量，函数值还可以更小，没有找到最优值 zj_cj=whole(m+1,2:(m+n+1)); % 把所有的判别数的值赋给zj_cj [maxzj_cj , k]=max(zj_cj); % 求得所有判别式中最大的那个数的下标是k,也就是进基变量的下标 % 从左往右进行比较，如果有两个一样的最大值，则选择最左边的那个 yk=whole(1:m,k+1); % 但是由于whole矩阵还有第一列 % 是基变量下标subscript，所以如果在whole矩阵中引用这一列，标号应该是k+1 % yk位1到m行，k+1列的数 [yk1 ,J]=find(yk'>0); % 找到yk中>0的那些数 if isempty(yk1) error('问题不存在最优解'); % 如果没有大于0的，表明不存在最优解 end b_=whole(1:m,m+n+2); %对应的是b的一个列向量 [b_yk1 , J1]=min(b_(J)./yk(J)); %b_yk1表示最小的那个数，J1表示b_(J)./yk(J)中最小数所在列位置(./表示的是对应元素的除法运算) r=J(J1); % 找到进基列的那个变化后非零分量的行数r y=whole(:,k+1); % y是进基变量对应的整一列向量 E=invmatrix(y,r); % 下面要对整个单纯性表（不包括subscript)做初等行变换，使得y向量变成(0,0,...1,0,...0)T,1在r个分量位置 % 由于初等行变换相当于左乘一个可逆矩阵，这一个子函数invmarrix就是求出这一个可逆矩阵 whole1=E*whole(:,2:(m+n+2)); % 做初等行变换（左乘可逆矩阵） subscript(r)=k; % 更新subscript whole=[[subscript;0] whole1]; % 更新whole end % 第一阶段求解最优解结束，进入判断三种情况 %情况一：最后的最优值大于0，无可行解 %情况二：最后的最优解等于0，基变量全在x1,x2,...,xn中，则x=(x1,x2,...,xn)T就是初始基可行解 %情况三：最后的最优解等于0，基变量不全在x1,x2,...,xn中，yk仍然是基变量，则要把yk变为非基变量 if abs(whole(end,end))<0.00001 % 如果最优值=0，表明原问题存在可行解，否则不存在 [s1,s2]=find(subscript'>n); % s2保存subscript中大于n的下标，由于大于n的基变量为人工变量，要用原有的变量将其驱赶 if ~isempty(s2) % 如果s2是空，则不存在大于n的下标，直接把whole的最后一行和人工变量所在列去掉就行，见下面 p=size(s2,2); % p是s2中的列数，在这相当于求s2元素个数，表示基变量的选取选择人工变量的有s2个 del=zeros(1,p); % 如果有一个人工变量是基变量但是这一行对应的所有原来的变量的值均为0，即无法驱赶，则将这一行整个删掉,用del记录 for i=1:p ind=find(whole(s2(i),2:n+1)); % 寻找这一行中x1到xn列中非零元素 if ~isempty(ind) % 如果没找到，则删除这一行，见下面，找到了，就取第一个 k1=ind(1);% 在whole中位置是k1+1; E1=invmatrix(whole(:,k1+1),s2(i)); % 驱赶的过程和上面进基变量替换离基变量过程类似，初等行变换， whole1=E1*E*whole(:,2:(m+n+2)); % 相当于左乘一个矩阵 subscript(s2(i))=k; % 更新subscript whole=[[subscript;0] whole1]; % 更新whole else del(i)=1; % 记录要删除的这一行的位置 end end ind2=setdiff(1:m,s2(logical(del))); % setdiff(a,b) 相当于集合减法，a-b=a-(a交b) a=1:m而whole有m+1列，相当于同时把whole的最后一行去掉 % ind2得到删除后剩下的下标 whole=whole(ind2,[1:(n+1) m+n+2]); % 新的whole矩阵 else whole=whole(1:(end-1),[1:(n+1) m+n+2]); % 如果不需要驱赶，则直接把whole的最后一行和人工变量所在列去掉 end else error('原问题不存在可行解'); end % 已经通过两阶段法的第一阶段找到初始基本可行解 % 把whole矩阵的判别式列补上 [m1 , n1]=size(whole); zj_cj2=c(whole(:,1))*whole(:,2:(n1-1))-c; % 计算判别式，注意这里把基变量和非基变量放在一起算，如果是基变量，zj-cj2这样算的结果就等于0 b_2=c(whole(:,1))*whole(:,n1); % 计算函数值，这里选出基变量所对应的c的值c(whole(:,1))来求判别式 whole2=[whole;[0 zj_cj2 b_2]]; % whole2是第二阶段的单纯形表 subscript=whole(1:m1,1); % 第二阶段求最优解，和第一阶段类似 % whole2 的大小: m1+1 * n1 while any(whole2(m1+1,2:(n1-1))>0.00001) zj_cj=whole2(m1+1,2:(n1-1)); [maxzj_cj , k]=max(zj_cj); % 在whole里面进基变量的位置是k+1; yk=whole2(1:m1,k+1); [yk1 , J]=find(yk'>0); % 注意 find函数中 行向量和列向量的返回值不同，如果是列向量，则 % 第一个返回值是下标，行向量是第二个返回值是下标 if isempty(yk1) error('问题不存在最优解'); end b_=whole2(1:m1,n1); [b_yk1 , J1]=min(b_(J)./yk(J)); r=J(J1); y=whole2(:,k+1); E=invmatrix(y,r); whole3=E*whole2(:,2:(n1)); % 更新subscript subscript(r)=k; whole2=[[subscript;0] whole3] end %求最优解所对应的x [m3,n3]=size(whole2); subscript=whole2(1:m3-1,1); s=n3-2; x=zeros(1,s); for i=1:m3-1 x(subscript(i))=whole2(i,n3); end end