运筹学 matlab实现单纯形法

文章目录

一,实验内容
二，建立的线性规划模型：
三，实验方法与步骤：
- 1，化为标准型：
- 2，在matlab中实现求解标准型的程序。
- - a,定义输入，用单纯形表计算：
  - b,计算步骤：
  - - 1，根据数学模型确定初始可行基和初始基可行解，建立初始单纯形表。
    - 2，计算各非基变量检验数，并判断是否是最优解：
    - 3，判断此问题是否无界：
    - 4，确定换入换出变量并记录位置以及更新N矩阵：（N矩阵用来存取的基变量的位置）
    - 5，以换入换出变量找到的alka_{lk}alk位置进行迭代即利用高斯消去法或称为旋转运算。
  - c,实现matlab程序 ssimplex.m：
四，实验结果：
五，实验结果分析：
六，附录：部分其余样例

一,实验内容

二，建立的线性规划模型：

maxz=0.9x1+1.4x2+1.9x3+0.45x4+0.95x5+1.45x6−0.05x7+0.45x8+0.95x9max z=0.9x_1+1.4x_2+1.9x_3+0.45x_4+0.95x_5+1.45x_6-0.05x_7+0.45x_8+0.95x_9 maxz=0.9x1+1.4x2+1.9x3+0.45x4+0.95x5+1.45x6−0.05x7+0.45x8+0.95x9

{−25x1+35x2+35x3≤0−15x1−15x2+45x3≤0−1720x4+320x5+320x6≤0−35x4−35x5+25x6≤0−12x7−12x8+12x9≤0x1+x4+x7≤2000x2+x5+x8≤2500x3+x6+x9≤1200x1,x2,x3,x4,x5,x6,x7,x8,x9≥0\left\{\begin{matrix} -\frac{2}{5}x_1+\frac{3}{5}x_2+\frac{3}{5}x_3&\le 0\\ -\frac{1}{5}x_1-\frac{1}{5}x_2+\frac{4}{5}x_3&\le 0\\ -\frac{17}{20}x_4+\frac{3}{20}x_5+\frac{3}{20}x_6&\le 0\\ -\frac{3}{5}x_4-\frac{3}{5}x_5+\frac{2}{5}x_6&\le 0\\ -\frac{1}{2}x_7-\frac{1}{2}x_8+\frac{1}{2}x_9&\le 0\\ x_1+x_4+x_7 &\le 2000\\ x_2+x_5+x_8&\le 2500\\ x_3+x_6+x_9&\le 1200\\ x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9 &\ge 0 \end{matrix}\right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧−52x1+53x2+53x3−51x1−51x2+54x3−2017x4+203x5+203x6−53x4−53x5+52x6−21x7−21x8+21x9x1+x4+x7x2+x5+x8x3+x6+x9x1,x2,x3,x4,x5,x6,x7,x8,x9≤0≤0≤0≤0≤0≤2000≤2500≤1200≥0

三，实验方法与步骤：

1，化为标准型：

{−25x1+35x2+35x3+x10=0−15x1−15x2+45x3+x11=0−1720x4+320x5+320x6+x12=0−35x4−35x5+25x6+x13=0−12x7−12x8+12x9+x14=0x1+x4+x7+x15=2000x2+x5+x8+x16=2500x3+x6+x9+x17=1200x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17≥0\left\{\begin{matrix} -\frac{2}{5}x_1+\frac{3}{5}x_2+\frac{3}{5}x_3+x_{10}&= 0\\ -\frac{1}{5}x_1-\frac{1}{5}x_2+\frac{4}{5}x_3+x_{11}&= 0\\ -\frac{17}{20}x_4+\frac{3}{20}x_5+\frac{3}{20}x_6+x{12}&= 0\\ -\frac{3}{5}x_4-\frac{3}{5}x_5+\frac{2}{5}x_6+x_{13}&= 0\\ -\frac{1}{2}x_7-\frac{1}{2}x_8+\frac{1}{2}x_9+x_{14}&= 0\\ x_1+x_4+x_7+x_{15} &= 2000\\ x_2+x_5+x_8+x_{16}&=2500\\ x_3+x_6+x_9+x_{17}&= 1200\\ x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9,x_{10},x_{11},x_{12},x_{13},x_{14},x_{15},x_{16},x_{17} &\ge 0 \end{matrix}\right. ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧−52x1+53x2+53x3+x10−51x1−51x2+54x3+x11−2017x4+203x5+203x6+x12−53x4−53x5+52x6+x13−21x7−21x8+21x9+x14x1+x4+x7+x15x2+x5+x8+x16x3+x6+x9+x17x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17=0=0=0=0=0=2000=2500=1200≥0

2，在matlab中实现求解标准型的程序。

a,定义输入，用单纯形表计算：

% 单纯形法matlab程序-ssimplex
% 求解标准型线性规划:max c*x; s.t. A*x=b; x>=0
% 本函数中的A是单纯初始表，包括:最后一行是初始的检验数，最后一列是资源向量b
% N是初始的基变量的下标
% 输出变量sol是最优解, 其中松弛变量（或剩余变量）可能不为0
% 输出变量val是最优目标值，kk是迭代次数

b,计算步骤：

1，根据数学模型确定初始可行基和初始基可行解，建立初始单纯形表。

即初始输入的A矩阵和N矩阵。

A=[ -0.4  0.6  0.6  0    0    0    0    0    0   1 0 0 0 0 0 0 0    0;-0.2 -0.2  0.8  0    0    0    0    0    0   0 1 0 0 0 0 0 0  0;0     0    0   -0.85 0.15 0.15 0    0    0   0 0 1 0 0 0 0 0  0; 0     0    0   -0.6 -0.6  0.4  0    0    0   0 0 0 1 0 0 0 0 0;0     0    0    0    0    0   -0.5 -0.5  0.5 0 0 0 0 1 0 0 0  0;1     0    0    1    0    0    1    0    0   0 0 0 0 0 1 0 0  2000;0     1    0    0    1    0    0    1    0   0 0 0 0 0 0 1 0   2500;0     0    1    0    0    1    0    0    1   0 0 0 0 0 0 0 1   1200;0.9   1.4  1.9  0.45 0.95 1.45 -0.05 0.45 0.95 0 0 0 0 0 0 0 0    0];
N=[10 11 12 13 14 15 16 17];

2，计算各非基变量检验数，并判断是否是最优解：

初始时检验数即为开始的max z里面的系数，后面在转轴运算的时候通过如下代码计算。

%i=mA时此时就是更新检验数
for i=1:mAif i~=outb   %~=为不等号A(i,:)=A(i,:)-A(outb,:)*A(i,inb);end
end

3，判断此问题是否无界：

        for i=1:nA-1if A(mA,i)>0&A(1:mA-1,i)<=0 % 第i个为换入变量但是无相应的换出变量。问题有无界解disp('have infinite solution!');flag=0;break;endend

4，确定换入换出变量并记录位置以及更新N矩阵：（N矩阵用来存取的基变量的位置）

代码较长且简单直接看ssimplex.m文件即可。

5，以换入换出变量找到的alka_{lk}alk位置进行迭代即利用高斯消去法或称为旋转运算。

            % 以下进行转轴运算A(outb,:)=A(outb,:)/A(outb,inb);for i=1:mAif i~=outb   %~=为不等号A(i,:)=A(i,:)-A(outb,:)*A(i,inb);endend

6，重复1~5的直到得到最优解。

c,实现matlab程序 ssimplex.m：

function [sol,val,kk]=ssimplex(A,N)
[mA,nA]=size(A);  %返回矩阵 行数和列数
kk=0; % 迭代次数
flag=1;
while flagkk=kk+1;if A(mA,:)<=0 % 已找到最优解  :表示所有的行或者列（取决于位置在哪）flag=0;sol=zeros(1,nA-1);  %生成一个1行nA-1列的矩阵for i=1:mA-1sol(N(i))=A(i,nA);  %这是将单纯形表的最右边一列的b赋给基变量endval=-A(mA,nA);    %保存结果z的值elsefor i=1:nA-1if A(mA,i)>0 & A(1:mA-1,i)<=0 % 第i个为换入变量但是无相应的换出变量。问题有无界解disp('have infinite solution!');flag=0;break;endendif flag % 还不是最优表，进行转轴运算temp=0;for i=1:nA-1if A(mA,i)>temptemp=A(mA,i);inb=i; % 进基变量的下标endend  %寻找检验数最大的为换入变量，记录值，并记录下标sita=zeros(1,mA-1)*nan;  %寻找换出变量for i=1:mA-1sita(i)=-1;endfor i=1:mA-1if A(i,inb)>0sita(i)=A(i,nA)/A(i,inb);endendtemp=inf;  %寻找最小的sita来找出换出变量for i=1:mA-1if sita(i)>=0&sita(i)<temp                    temp=sita(i);outb=i; % 出基变量下标endend% 以下更新NN(outb)=inb;% 以下进行转轴运算A(outb,:)=A(outb,:)/A(outb,inb);for i=1:mAif i~=outb   %~=为不等号A(i,:)=A(i,:)-A(outb,:)*A(i,inb);endendendend
end

输入可以为：

A=[ -0.4  0.6  0.6  0    0    0    0    0    0   1 0 0 0 0 0 0 0    0;-0.2 -0.2  0.8  0    0    0    0    0    0   0 1 0 0 0 0 0 0  0;0     0    0   -0.85 0.15 0.15 0    0    0   0 0 1 0 0 0 0 0  0; 0     0    0   -0.6 -0.6  0.4  0    0    0   0 0 0 1 0 0 0 0 0;0     0    0    0    0    0   -0.5 -0.5  0.5 0 0 0 0 1 0 0 0  0;1     0    0    1    0    0    1    0    0   0 0 0 0 0 1 0 0  2000;0     1    0    0    1    0    0    1    0   0 0 0 0 0 0 1 0   2500;0     0    1    0    0    1    0    0    1   0 0 0 0 0 0 0 1   1200;0.9   1.4  1.9  0.45 0.95 1.45 -0.05 0.45 0.95 0 0 0 0 0 0 0 0    0];
N=[10 11 12 13 14 15 16 17];
[sol,val,kk]=ssimplex(A,N)

或者将A,N矩阵加入在函数中求得结果。

四，实验结果：

结果为：

sol =1.0e+03 *1.5267    1.0178         0    0.4733    1.4822    1.2000         0         0         0         0    0.5089         0    0.6933         0         0         0         0val =6160kk =9

以下是matlab截图结果：

五，实验结果分析：

经过计算，实验结果正确，并使用多个样例测试程序皆有正确结果。

六，附录：部分其余样例

%样例1    4360% 结果为  x = 320.0000  360.0000   0   0  20.0000 390.000 ]
A=[ 2 1  1 0 0 0 1000;3 4  0 1 0 0 2400;1 1  0 0 1 0 700;1 -1 0 0 0 1 350;8 5  0 0 0 0 0];
N=[3 4 5 6];
[sol,val,kk]=ssimplex(A,N)   /不用加分号

书本例题：
%样例3
A=[ 1 2 1 0 0 8;4 0 0 1 0 16;0 4 0 0 1 122 3 0 0 0  0];
N=[3 4 5];
[sol,val,kk]=ssimplex(A,N)