MATLAB实验
实验一
1.判断是否存在单位矩阵，并输出列
m=3;n=8; A = 10*rand(m,n);I = eye(m,m);
randIndex = randperm(size(A,2));
A(:,randIndex(1:m))=I;
%准备数据，想了好久这是在干嘛
count=0;
C = nchoosek(randIndex,m);
%R = rref(A)
%%%
for i=1:size(C,1)
if (A(:,C(i,:))==I)
%disp("find eye!");
count = count + 1;
Column = C(i,:);
end
end
%%%
if (count == 0)
disp("no eye!");
else
disp("find eye!");
Column
end
组合排列函数
nchoosek(n,m)
含义：从n个元素中取出m个元素的所有组合。
将矩阵化成行最简形的命令是rref或rrefmovie。
函数 rref或rrefmovie格式
R = rref(A) %用高斯—约当消元法和行主元法求A的行最简行矩阵R
2.高斯消元解方程
function x=gauss_elim(A,b)
%参量说明：A,系数矩阵；B，常数列向量；zg，增广矩阵
%将增广矩阵化为上三角，再回带求解x
%此方法较为常规，将zg(k,k)元素乘以-zg(i,k)/zg(k,k)加到第i行
%从1:n-1列，主对角元素的以下行，通过两层循环来遍历
zg=[A,b];
zj=rref(zg);
n=length(b);
ra=rank(A);
rz=rank(zg);
temp=rz-ra;
if temp>0
disp('无解');
return;
end
if ra==rz
if ra==n
x=zeros(n,1);
for p=1:n-1
for k=p+1:n
m=zg(k,p)/zg(p,p);
zg(k,p:n+1)=zg(k,p:n+1)-m*zg(p,p:n+1);
end
end
%用第p层（从1到n-1）消元
b=zg(1:n,n+1);
A=zg(1:n,1:n);
x(n)=b(n)/A(n,n);
for q=n-1:-1:1
x(q)=(b(q)-sum(A(q,q+1:n)*x(q+1:n)))/A(q,q);
end
end
endfind（）
查找非零元素的索引和值
k = find(X)
k = find(X,n)
k = find(X,n,direction)返回一个包含数组 X 中每个非零元素的线性索引的向量。
实验二
实现函数function [xs,Bs,x_num]=BFS(A,b)，求基本可行解和对应的基（矩阵）。
测试程序
%BFS_test.m
A=[2,1,1,0,0;
1,1,0,1,0;
0,1,0,0,1;]
b=[10,8,7]'
[xs,Bs,x_num]=BFS(A,b)function [xs,Bs,x_num]=BFS(A,b)
%参量说明：A,系数矩阵；b，常数列向量；xs,基本可行解，Bs对应基矩阵，x_num，个数
clc %清屏
%clear all % 清除变量
close all % 关掉图像窗口
format RAT % 显示分数形式
%m = 3; %矩阵A的行数
%n =5; %矩阵A的列数% 随机生成一个m*n行满秩矩阵A测试,这里按照ppt给的例子设置
[m,n]=size(A);
x_num = 0;
idxs = nchoosek(1:n,m); %所有可能的基阵的下标组合
count = size(idxs,1);
for k = 1:1:count
idx = idxs(k,:);
Bs_ = A(:,idx);
if det(Bs_) == 0
%fprintf(['\n1.B不能构成基阵.\n']); % B是奇异矩阵，不能构成基阵
else
%fprintf(['\n1. B能构成基阵.\n']); % B不是奇异矩阵，能构成基阵
% 令X_N=0，解方程B_X*B+N*X_N=b,计算基本解为X_B=B^(-1)b,XN=0;
xs_ = zeros(1,n);
X_B = inv(Bs_)*b;
X_N = 0;
xs_(idx)=X_B;
%fprintf('2. 该基阵对应的基本解：xs= [%s].\n',rats(xs));% 判断是否为基本可行解
[r,c]=size(xs_);
flag=0;%判断基本解是否为基本可行解
for i = 1:c
if (xs_(1,i)<0)
flag=flag+1;
end
end
if (flag == 0)
%fprintf(['\n3. xs是基本可行解.\n']);
x_num =x_num + 1;
xs(x_num,:)=xs_;
Bs(:,:,x_num)=Bs_;
else
%fprintf(['\n3. xs不是基本可行解.\n']);
end
end
end实验三
测试程序
%test
A=[ 2 -3 2 1 0;
1/3 1 5 0 1];
b=[15 20]';
c=[1 2 1 0 0];
format rat %元素使用分数表示
% [x_opt,fx_opt,iter] = Simplex_eye(A,b,c);函数实现
function [x_opt,fx_opt,iter] = Simplex_eye(A,b,c)
%参量说明：x_opt为最求解，fx_opt为最优函数值，iter为迭代次数
% max z=(c)'x
% Ax=b,x>=0
% r(A)=m,A大小(m*n),(n>=m),b>=0(m*1),c(1*n)这里有问题输入可以是，
%可能是方便输入,直接输入的c其实是是转置后（c'）的行向量
%A是化作标准型后的系数矩阵，
format rat %元素使用分数表示
[m,n]=size(A)%m约束条件个数, n决策变量数
I = eye(m,m);
randIndex = randperm(size(A,2));%A的列数整数的随机排列
C = nchoosek(randIndex,m);%列数的排列组合的集合
ind_B=[];
for i=1:size(C,1)%遍历每一种组合
if (A(:,C(i,:))==I)%找到单位矩阵
disp("find eye!");
ind_B = C(i,:);%找到基变量的列数序号
end
end
ind_N = setdiff(1:n, ind_B); %非基变量的索引，从小到大排序
iter=0;%记录迭代次数
%开始循环迭代
while(true)
%取基变量，非基变量取值为0
%选择初始可行基
x0=zeros(n,1);
x0(ind_B) = b;%初始基可行解
cB = c(ind_B);%计算检验数，选一个检验数最大的非基变量作为换入变量
cN=c(ind_N)
%注意基变量检验数为0
B=A(:,ind_B)%基本矩阵
Pjlist=A(:,ind_N)
wlist=inv(B)*Pjlist;%wlist=y_k
zlist=(cB)*inv(B)*Pjlist;
clist=cN-zlist;
clist;
[maxVal maxInd] = max(clist);%找到换入变量xk，k=maxInd
k=maxInd;
%判断是否已经是最优解
if all(clist(:)<=0)%clist中检验值是否都<=0
x_opt = x0;
fx_opt = c*x_opt;
return
end
%判断是否已经是无界
if all(wlist(:)<=0)%clist中元素是否都<=0
x_opt = [];
fx_opt = [];
disp('问题无界')
return
end
%选一个基变量作为换出变量
Pk=A(:,k)
yk=inv(B)*Pk
% yk(yk<=0)=0.0000001;
thetalist=b./yk;
%a_ik>0，b一定大于0，如果a_ik<=0就让b_i/a_ik变得无穷大
thetalist(thetalist<=0) = 10000;
[~,minId] = min(thetalist);%确定出基变量xr索引r,这里写错好多次改好久
r=ind_B(minId)
% 换基
ind_B(ind_B == r) = k; %新的基变量索引
ind_N = setdiff(1:n, ind_B); %非基变量索引
% 更新A和b,化最简型
A(:,ind_N) = A(:,ind_B) \ A(:,ind_N)%基矩阵的逆乘以非基矩阵
b = A(:,ind_B) \ b %基矩阵的逆乘以b
A(:,ind_B) = eye(m,m) %基矩阵更新为单位矩阵
iter=iter+1;
end
end
花了好久把单纯形法和单纯形表代码写出来，结果正确！
1.setdiff函数
set difference.
C=setdiff(A,B) for vector A and B, return the values in A that are not in B with no repetitions. C will be sorted.
对于向量A,向量B，C=setdiff(A,B)函数返回在向量A中却不在向量B中的元素，并且C中不包含重复元素，并且从小到大排序。
————————————————
版权声明：本文为CSDN博主「风景不在对岸wj」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/u011089523/article/details/79986511
参考资料
CSDN一搜全是上大智科的学长，学长太强了。我的代码水平还是不够。
以收集上大智科学长的运筹学博客为乐，这些代码可以一直传承下去。
“最近在上《运筹与优化》这门课，讲到了单纯形法这部分，老师让我们上机用Matlab实现单纯形法，数学公式实现过程看不懂的我就暴力模拟了单纯形法的整个过程，包含无界解等情况。”——努力变成大白的小白学长
学长写的代码好复杂，不如上面18级学长写的清楚。
这个17级的学长写了120行，其实只要写50行不到（20行注释），可以年年优化，减少代码行数。
纸质作业周四交
