线性规划——单纯形法

最新推荐文章于 2024-09-29 08:00:00 发布

CaptainHarryChen

最新推荐文章于 2024-09-29 08:00:00 发布

阅读量2.2k

点赞数 3

分类专栏： OI知识总结计算机数学文章标签：数学

本文链接：https://blog.csdn.net/can919/article/details/88567924

版权

OI知识总结同时被 2 个专栏收录

37 篇文章 3 订阅

订阅专栏

计算机数学

29 篇文章 0 订阅

订阅专栏

问题

标准型问题：

最大化 $\sum\limits_{j=1}^{n}c_jx_j$
满足约束 $\sum\limits_{j=1}^{n}a_{ij}\cdot x_j\le b_i\ (i=1,2,...,m)$
且 $x_j\geq 0\ (j=1,2,3,...m)$

可以转换为松弛型问题：

最大化 $\sum\limits_{j=1}^{n}c_jx_j$
满足约束 $x_{i+n}=b_i-\sum\limits_{j=1}^{n}a_{ij}\cdot x_j\ (i=1,2,...,m)$
且 $x_j\geq 0\ (j=1,2,3,...n+m)$

把约束中 $x_{i+n}$ 称为非基变量， $x_j$ 称为基变量

算法

基本思路

初始找出一个满足所有约束的初始解，通过这个解不断的更新，找到更优的解，直到找不到位置（类似爬山算法）。

可以发现这样一定能找到最优解，因为如果把求解的x当做n维空间的一个点，约束条件一定在n维空间中框出一个n维凸“多面体”，一个解为凸“多面体”的一个顶点。凸面体从哪一个发现看都是单峰的，用爬山的算法是肯定能找到最优解的。

例子：
最大化： $S=3x_1+x_2+2x_3$
约束：
$\begin{aligned} x_4 &= 30 − x_1 − x_2 − 3x_3\\ x_5 &= 24 − 2x_1 − 2x_2 − 5x_3\\ x_6 &= 36 − 4x_1 − x_2 − 2x_3\\ x_j&\geq 0\ (j=1,2,3) \end{aligned}$
显然可以看出初始解(0,0,0,30,24,36)

把第三个式子化为：
$x_1=9-\frac 1 4 x_6-\frac 1 4 x_2-\frac 1 2 x_3$
再将 $x_1$ 代入其它每个式子中去，可以发现 $S=27-\frac 3 4 x_6+\frac 3 4 x_2+\frac 3 2x_3$
如此转换，如果我们把S式中所有的x系数转为负数，就可以把常数项当做最优解。

这种把一个非基变量和基变量交换的操作叫做转轴（pivot）。

转轴

转轴操作可以把交换一个非基变量和一个基变量
设交换非基变量 $x_N$ 和基变量 $x_B$
$x_B=b_i-\sum\limits_{j=1}^na_{ij}\cdot x_j$
则转为 $x_N=(b_i-x_B-\sum\limits_{j=1,j\ne N}^na_{ij}\cdot x_j)/a_{iN}$

上面的例子转轴 $x_1,x_6$ 之后为
最大化： $S=27-\frac 3 4 x_6+\frac 3 4 x_2+\frac 3 2x_3$
约束：
$\begin{aligned} x_4 &= 21 − \frac 3 4 x_6 + \frac 1 4 x_2 + \frac 1 2 x_3\\ x_5 &= 6 +\frac 1 2 x_6− \frac 3 2 x_2 − 4x_3\\ x_1&=9-\frac 1 4 x_6-\frac 1 4 x_2-\frac 1 2 x_3\\ x_j&\geq 0\ (j=1,2,3) \end{aligned}$

可以发现S至少为27了

simplex操作

即不停的转轴，使得S式中所有的参数都为负。

每次寻找S中参数最大的 $x_j$ （贪心的想，用最大的参数可以使答案扩得更大），然后从约束条件中，找到对 $x_j$ 限制最大的约束 $i$ ，即 $\frac {b_i} {a_{ij}}$ 最小，转轴 $x_j$ 与 $x_{i+n}$ （选择限制最大的约束，使得转轴后，其余约束条件中的 $b_i$ 非负）

寻找初始解

当满足所有 $b_i\geq 0$ ，显然可以把 $x_j=0,x_{i+n}=b_i$ 当做一个可行的初始解。

如果存在 $b_i<0$ ，有两种方法

法一

论文方法：

可以通过引入辅助变量 $x_0$ ，经过转轴操作将 $b_i$ 转换为 $\geq 0$ 。
最大化 $x_0$
满足约束 $x_{i+n}=b_i-\sum\limits_{j=1}^{n}a_{ij}\cdot x_j+x_0\ (i=1,2,...,m)$
如果得到最优解 $x_0=0$ ，则这个 $x_0$ 对原来问题的答案没有影响，否则原问题无解。

选择最小的 $b_i$ ，设为第k个约束，把那个约束中 $x_0$ 进行转轴操作换出来，则
$x_0=-b_k+x_{k+n}+\sum\limits_{j=1}^{n}a_{kj}\cdot x_j$
其它约束变为
$x_{i+n}=b_i-b_k+x_{k+n}+\sum\limits_{j=1}^{n}(a_{kj}-a_{ij})\cdot x_j$
因为 $b_k$ 为最小的b，所以这个式子中 $b_i-b_k$ 一定为非负数，从而把所有的约束中常数项变为非负。

此时再从要求最大化的式子中，选一个参数为正的 $x_p$ ，将其与 $x_0$ 转轴。

重复这样的操作，直到所有 $b_i$ 全为正。

实现简单的法二

每次随机找到一个<0的 $b_i$ ，在从 $b_i$ 的约束条件中找一个系数大于0的 $x_j$ ，转轴 $x_{i+n},x_j$ ，就可以把当前 $b_i$ 变为正的（虽然有可能把其它条件的 $b_i$ 变为负数，所以要随机。。。）
多次执行随机操作，基本上可以把所有 $b_i$ 变为非负数。（UOJ拿97分）

法三

把法二的 $b_i$ ，每次选择最小的负数 $b_i$ 来进行操作也能过，并不知道为什么。。。

代码

UOJ179

#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
#include<algorithm>
using namespace std;
const int MAXN=55;
const double EPS=1e-8;

int dcmp(double a,double b)
{
    if(fabs(a-b)<=EPS)
        return 0;
    return a<b?-1:1;
}

int n,m,type;
double A[MAXN][MAXN];
int id[MAXN],tp[MAXN];

void pivot(int r,int c)
{
    swap(id[r+n],id[c]);
    double tmp=-A[r][c];
    A[r][c]=-1;
    for(int j=0;j<=n;j++)
        A[r][j]/=tmp;
    for(int i=0;i<=m;i++)
        if(i!=r&&dcmp(A[i][c],0)!=0)
        {
            tmp=A[i][c];
            A[i][c]=0;
            for(int j=0;j<=n;j++)
                A[i][j]+=tmp*A[r][j];
        }
}

void solve()
{
    for(int i=1;i<=n;i++)
        id[i]=i;
    for(;;)
    {
        int i=0,j=0;
        for(int k=1;k<=m;k++)
            if(dcmp(A[k][0],0)==-1&&(!i||rand()&1))
            	i=k;
        if(i==0)
            break;
        for(int k=1;k<=n;k++)
            if(dcmp(A[i][k],0)==1&&(!j||rand()&1))
           		j=k;
        if(j==0)
        {
            puts("Infeasible");
            return;
        }
        pivot(i,j);
    }
    for(;;)
    {
        int i=0,j=0;
        double mx=0,mn=1e9;
        for(int k=1;k<=n;k++)
            if(dcmp(A[0][k],mx)==1)
                mx=A[0][k],j=k;
        if(j==0)
            break;
        for(int k=1;k<=m;k++)
            if(dcmp(A[k][j],0)==-1&&dcmp(-A[k][0]/A[k][j],mn)==-1)
                mn=-A[k][0]/A[k][j],i=k;
        if(i==0)
        {
            puts("Unbounded");
            return;
        }
        pivot(i,j);
    }
    printf("%.9f\n",A[0][0]);
    for(int i=n+1;i<=n+m;i++)
        tp[id[i]]=i-n;
    if(type)
        for(int i=1;i<=n;i++)
            printf("%.9f%c",tp[i]?A[tp[i]][0]:0," \n"[i==n]);
}

int main()
{
    scanf("%d%d%d",&n,&m,&type);
    for(int i=1;i<=n;i++)
        scanf("%lf",&A[0][i]);
    for(int i=1;i<=m;i++)
    {
        for(int j=1;j<=n;j++)
        {
            scanf("%lf",&A[i][j]);
            A[i][j]*=-1;
        }
        scanf("%lf",&A[i][0]);
    }

    solve();

    return 0;
}