算法实验之线性规划解决配料问题

公式和图片不显示...

1 问题描述

    食谱/配料问题：设有n种配料,每种配料含有m种营养成分,用 表示一个单位的第 j (j<=n)种配料中含有的第 i (i<=m)种营养成分数量,用 表示对第 i 种营养成分的最低需要量,表示第 j 种食品的单价,表示所用的第 j 种食品的数量，配料问题的求解即解决应如何搭配食物（即求解）, 一方面满足种营养成分的需要，同时使食物的总成本最低的问题。

 

2 问题的相关背景描述

这类问题其实是在计划任务确定的情况下，如何统筹安排，用作少的资源来实现任务要求，可以抽象为一个求极小值的问题。对于社会实际问题而言，这类问题研究和应用的内容是实现系统的投入产出的问题，用最少的劳力和物力等投资消耗，获利更多更好的社会需求产品。

解决配料问题常用的方法是图解法和单纯形法。图解法虽然清晰直观、方便简单，但是只适用于二维的线性规划问题，在我们的高中时期经常用到。单纯形法的优点是可以适用于所有的线性规划问题，但是单纯形法中涉及大量不同的算法，为了针对不同的线性规划问题，计算量大，复杂繁琐。除此之外还有内点法、椭球法等求解线性规划问题的算法。

 

参考文献：第1卷第1期 西昌学院学报·自然科学版 Vo1．1，No．12013年06月 Journal of Xichang College·Natural Science Edition Jun．，2013

3.问题分析

解决线性规划问题，首先要对线性规划问题建立数学模型，确定目标函数和相应的约束条件。

从实际问题中建立数学模型一般有以下步骤：

（1）、根据所求目标的影响因素找到决策变量；

（2）、由决策变量和所求目标的函数关系确定目标函数；

（3）、由决策变量所受的限制条件确定决策变量所需要满足的约束条件。

3.1单纯形法（主要）

算法步骤：

（1）、将非标准型线性规划化成标准型，加入松弛变量；

（2）、求出线性规划的初始基可行解，列出初始单纯形表；

（3）、进行最优解检验。如果表中所有检验数小于等于零，则表中的基可行解就是问题的最优解，计算停止，否则进行下一步。

（4）、从一个基可行解转换到另一个目标值更大的基可行解，列出新的单纯形表。确定换入基的变量，当有一个以上检验数大于0时，一半选择最大的一个检验数，其对应的作为换入变量；确定换出变量，根据计算并选择，选最小的对应基变量作为换出变量；用换入变量替换基变量中的换出变量，得到一个新的基。对应新的基可以找出一个新的基可行解，并相应地可以画出一个新的单纯形表。

（5）、重复（3）（4）步骤直到计算结束为止。

算法评估：平均计算量统计O(n)

3.2椭球法

算法步骤：

（1）、变换问题提法：原问题：,    对偶问题：, ，于是知，若有最优解则构造的下述复合不等式必成立：,;

（2）、变换上述不等式并试图求解，然后构造一个大的球体，使其必包含不等式可行解（若存在的话），对球心判断是否为可行解，若是，结束，否则，切割球体（切去肯定不包含可行解的部分）直至找到可行解，或证明无可行解。

算法评估：并不实用，计算量

3.3内点法

算法步骤：

算法评估：计算量 平均计算量统计O(lgn)

3.4蛮力法

算法步骤：根据条件将解空间缩小后，在解空间中搜索最优解

算法评估：在特殊情况下可能会表现突出，但是更普遍的情况下并不适用。

 

 

由于单纯形法更具有普遍性，在解决配料问题等各类线性规划问题方面也表现突出，所以单纯形法是我们组重点研究的对象。

 

 

4.计算模型

1、  化为标准形式

…

 

 

2、  化为松弛形式

此时依然成立, 我们将这些变量称为非基变量，它们构成的集合记为N。将这些变量称为基变量，它们构成的集合记为B。简单地理解，非基变量能够由基变量唯一确定。

s.t.

3、  初始化过程

如果b向量并不全为非负，则我们需要通过初始化过程来找到一个可行解，然后才可以使用最优化过程进行优化。当然，此时原线性规划不一定存在可行解。具体的方法是，加入一个新的非基变量，并在原线性规划的基础上构造一个新的辅助的线性规划。

 注意这里N集合并不包含。然后，选择一个基变量使得最小，执行转轴操作pivot(d, 0)。不难证明该操作过后所有的b值全部非负。然后，使用前文中所述的最优化过程求解该辅助线性规划。由于非负，因此该线性规划的答案非正。如果答案为负数，则说明原线性规划不可能让所有的基变量都非负，因此原线性规划无可行解。否则，只要令所有变量为0，并去掉变量，就可以得到可行解。在从辅助线性规划转化到原来的线性规划的过程中，如果已经是非基变量，则可以将其从约束条件和目标函数中直接去掉。否则，需要任取一个非基变量执行pivot(0，e)，将变为非基变量。由于此时是基变量且=0，故一定成立，因此这个转轴操作不会破坏b向量的非负性。

4、  最优化过程

如果b向量所有元素非负，则显然我们只需要令所有的变量等于0，就可以得到一个可行解。在这种情况下，通过下述最优化过程，我们可以得到该线性规划的最优解，或者指出该线性规划的最优解为无穷大（不存在）。

1、任取一个非基变量，使得。

2、选取一个基变量，使得，且最小化

3、执行转轴操作pivot(d, e)，并转到第一步继续算法。

根据的最小性不难证明pivot(d, e)不会破坏b的非负性。因此将所有变量取0值仍然是可行解。同时，根据，我们发现v一定是不降的。这就达到了更新解的目的。不难发现，算法终止有两种情况：

1、对于所有的非基变量，c均非正。

2、对于某一个e，所有的均非正。

可以证明，对于第一种情况，我们已经得到了该线性规划的最优解。当前的v即为答案。严格证明比较复杂，但是直观上是很容易理解的。因为所有的非基变量都是非负的，而所有的c都是非正的，因此只要某个非基变量不为0，就会使得目标函数更小。对于第二种情况来说，很容易证明此时线性规划的最优解是无穷大。只要让其他所有变量均为0，变量为正无穷。由于所有的都非正，因此非基变量的非负性得到保证。同时由于，目标函数值为正无穷。

5、  选取入基变量离基变量

步骤1：选入基变量。如果所有≤0，则当前基本可行解为最优解，计算结束。否则，取>0，相应的非基本变量为入基变量。

步骤2：选离基变量。 对于入基变量，如果所有≤0(i=1,2, …, m)，则最优解无界，计算结束。否则，计算

选取基本变量x_k为离变量。

6、  转轴操作

如果,我们有

 

 

实例1、一饲养场饲养供实验用的动物,已知动物生长对饲料中的三种营养成分—蛋白质、

矿物质和维生素特别敏感。每个动物每天至少需要蛋白质70克,矿物质3克,维生素10毫克,该场能搞到五种饲料,每种饲料的成本如表4.2.1

饲料

A1

A2

A3

A4

A5

成本（元）

2

7

4

3

5

每一千克饲料中所含的营养成分如表4.2.2,我们希望确定既能满足动物需要,又使成本最低的饲料配方。

饲料

蛋白质

矿物质

维生素

A1

0.30

0.10

0.05

A2

2.00

0.05

0.10

A3

1.00

0.02

0.02

A4

0.60

0.20

0.20

A5

1.80

0.05

0.08

本题要求我们给出一个饲料配方,即确定在混合饲料中每种饲料的重量,因此变量应该是:

混合饲料中所含的第j种饲料的数量.

我们的目标是使得总成本最低。成本=0.1(2 +7 +4 +3 +5)。动物对蛋白质、矿物质及维生素的需求构成本问题的约束条件:

 

0.30 +2 + +0.6 +1.8 ≥70

0.1 +0.05 +0.02 +0.2 +0.05 ≥3

0.05 +0.12 +0.02 +0.2+0.08 ≥10

≥0, j=1,2,...,5

优化模型：

min 0.1(2+7+4+3+5)

0.1 +0.05 +0.02 +0.2 +0.05 ≥3

0.05 +0.12 +0.02 +0.2 +0.08 ≥10

≥0, j=1,2,...,5

 

实例2、某养猪场所用的混合饲料由 A 、 B 、 C三种配料组成，表5.17给出1单位各种配料所含的营养成分，单位成本以及1份混合饲料必须含有的各种营养成分，问如何配制饲料时成本最小？

配料

营养成分

单位成本

D

E

F

A

1

1/2

2

6

B

1

1/2

1

3

C

1

1/4

1

2

1份饲料应含量

20

6

10

 

解   设为混合饲料中第j中配料的含量，j=A,B,C,

目标函数,

满足:  

  

  

 

 

 

5.算法设计与描述

// 通用算法伪代码描述

// 目标函数的系数c[j], j from 1 to n

// 约束条件的系数a[i][j], i from 1 to m

// 约束条件的常数b[i]

 

 

InitializeSimple()

{

    输入变量数目 n 和约束数目 m;

    输入目标函数的系数向量 (c[j]);

    输入约束条件的系数矩阵 (a[i][j]);

    输入约束条件的常数项 (b[i]);

}

 

 

Simplex()

{

    InitializeSimple(); //初始化单纯形表

   

    while ( 目标函数中，存在系数 c[j] > 0)

       

        // 从非基本变量中选择入基变量

        for ( e in all j ) {

            if ( x[e]为非基本变量&& c[e] > 0)

                选择对应的 x[e] ;

                记录该入基变量下标e;

        }

           

        // 从非基本变量中选择离基变量

        for ( i <- 1 to m ) {

            if ( x[i] 为基本变量&& 该行中离基变量 x[e] 系数a[i][e] > 0 )

                找到使 min <- b[i] /a[i][e] 的值最小的 i, 即找到离基变量 x[i];

                记录该离基变量下标l;

        }

               

        if ( min 在精度允许范围内为无穷 )

            return "该线性规划无界";

        else

            //转轴变换

            GetPivot(l, e);

 

        // 当前解的下标 ans

        ans <- 1

        for ( i <- 1 to n )

            if ( x[i] 是基本变量 )

                更新当前解 x[ans] <-b[i];

}

 

 

GetPivot()

{

    // 为基本变量x[e]计算约束方程常数项

    b1[l] = b[l] / a[l][e];

    // 对基本变量 x[e] 所在方程中，其他非基本变量且非入基变量的系数 a[e][i] 更新

    for ( i <- 1 to n+m ) {

        if ( x[i] 为非基本变量&& x[i] 不是入基变量 )

        a[e][i] <- a[l][i] / a[l][e];

    }       

    //更新入基变量系数

    a[e][l] 变成原来的倒数;

 

 

 

    // 计算剩余的约束系数

    for ( i <- 1 to n+m ) {

        if ( x[i] 是基本变量&& x[i] 不是离基变量 ) {

            //更新常数项

            b1[i] <- b[i] - a[i][e] * b1[e];

            for ( j <- 1 to n+m ) {

                if ( x[j] 不是基本变量&& x[j] 不是入基变量 )

                    //更新系数

                    a1[i][j] <- a[i][j] -a[i][e] * a1[e][j];

            }

            // 更新离基变量的系数

            a1[i][l] <- (-1) * a[i][e] * a1[e][l];

        }

    }

 

 

 

    // 计算目标函数

    v1 <- v + c[e] * b1[e];

    for ( i <- 1 to n+m ) {

        if ( x[i] 不是基本变量&& x[i] 不是入基变量 )

            // 更新目标函数中非离基变量的系数

            c1[i] <- c[i] - c[e] * a1[e][i];

    }

    // 更新目标函数中离基变量的系数

    c1[l] <- (-1) * c[e] * a1[e][l];

 

 

   

    // 更新基本变量和非基本变量

    for ( i <- 1 to n+m ) {

        // 更新目标函数的

        c[i] <- c1[i];

        b[i] <- b1[i];

        // 更新约束方程系数

        for (j <- 1 to n+m )

            a[i][j] <- a1[i][j];

    }

}

 

6.算法分析

6.1实例分析：

在实例2中，转换为松弛型后，变量数目 n＝6，约束数目 m＝3，目标函数的系数向量(c[j])＝｛6 3 2｝；

约束条件的系数矩阵 a[i][j]＝｛-20 1 1 1，－6 1/2 1/2 1/4，－10 2 1 1｝；

约束条件的常数项 b[i]＝｛－20 －6 －10｝；

每次迭代的时间复杂度上界，在次迭代下，总的时间复杂度

＝

 

，总的空间复杂度为＝27

 

 

6.2抽象分析：

时间复杂度：

单纯形法的时间复杂度由两个因素决定：一是每次迭代中的步数，二是迭代次数。

 

算法分三个部分，初始化的InitializeSimplex()，转轴操作 GetPivot() 和函数Simplex()，根据伪代码，将每行语句（除了for循环的循环条件部分和函数调用）看作单位时间。

 

·InitializeSimplex:

由于只是输入的初始化，根据每行的输入，

 

·GetPivot()

对伪代码进行分析，共出现6次for循环，且有四个出现在两个嵌套循环，且分为四部分：

1、为基本变量计算约束方程常数项和约束方程系数

 

       2、计算剩余的约束系数

 

       3、计算目标函数值和系数

 

4、更新变量系数

 

可得GetPivot()运行一次的总时间

 

其中关键步骤为1和2，不难看出每次迭代的时间复杂度上界。

 

·Simplex()

设每次迭代过程里，目标函数中仍有c个系数严格大于0。

 

第次迭代（精度范围内有界）的复杂度如下：

 

故每次迭代时间复杂度如下：

 

在次迭代下，总的时间复杂度

 

由此推导与前人的证明里，单纯形法的最终时间复杂度上界为，根据经验得为关于n的指数级别。

 

空间复杂度：

在本算法的实现中，利用的辅助变量有一个长度为n的一维数组，一个长为m的一维数组，以及一个m*n的二维数组，故

得空间复杂度上界为。

7.实验结果分析

实例1:目标函数minZ=0.1(2 +7 +4 +3 +5)

约束条件:

0.30 +2 + +0.6 +1.8 ≥70

0.1 +0.05 +0.02 +0.2 +0.05 ≥3

0.05 +0.12 +0.02 +0.2+0.08 ≥10

≥0, j=1,2,...,5

解得    

实例2:目标函数,

约束条件:

  

  

  

 

解得 

 

8.软件使用说明

请输入线性规划问题的变量个数N：…

请输入线性规划问题的约束条件个数M：…

请输入目标函数各个变量的系数所构成的系数阵C：…

输入矩阵：…

参考文献

[1] 维基百科 https://zh.m.wikipedia.org/zh-cn/单纯形法

[2] 百度百科 内点法http://baike.baidu.com/

[3] 百度百科 椭球法http://baike.baidu.com/