运筹学 matlab实现单纯形法

weixin_52624015

已于 2022-03-22 22:35:29 修改

阅读量4.7k

点赞数 3

分类专栏：运筹学文章标签： matlab

于 2022-03-22 15:50:26 首次发布

本文链接：https://blog.csdn.net/weixin_52624015/article/details/123662712

版权

运筹学专栏收录该内容

3 篇文章 5 订阅

订阅专栏

文章目录

一,实验内容

$[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-8EH3ZUjI-1647959486435)(C:\Users\86150\AppData\Roaming\Typora\typora-user-images\image-20220321213956801.png)]$

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

$max 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$

$\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.$

三，实验方法与步骤：

1，化为标准型：

$max 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$

$\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.$

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:mA
    if i~=outb   %~=为不等号
       A(i,:)=A(i,:)-A(outb,:)*A(i,inb);
     end
end

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

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

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

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

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

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

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

c,实现matlab程序 ssimplex.m：

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


val =

        6160


kk =

     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 12
	2 3 0 0 0  0];
N=[3 4 5];
[sol,val,kk]=ssimplex(A,N)