线性规划（单纯型法）

最新推荐文章于 2024-08-07 09:24:57 发布

JK Chen

最新推荐文章于 2024-08-07 09:24:57 发布

阅读量623

点赞数

分类专栏：数论/数学知识点 ACM中的数学问题合集 All

本文链接：https://blog.csdn.net/jk_chen_acmer/article/details/101603156

版权

All 同时被 3 个专栏收录

629 篇文章 6 订阅

订阅专栏

数论/数学知识点

48 篇文章 6 订阅

订阅专栏

ACM中的数学问题合集

44 篇文章 33 订阅

订阅专栏

问题模型

有 $n$ 个变量 $x_1,x_2,...,x_n$ ， $m$ 个约束条件 $\sum k_ix_i= b_i$ 。

题目要求所有的解非负，即 $x_i\geq 0$ 。

求解目标方程 $\sum k_ix_i$ 的最大值。

模型转化

遇到 $x_1+x_2\leq 7$ ，可以增加一个非负松弛变量 $x_3$ ，变为： $x_1+x_2+x_3=7,x_3\geq 0$

遇到 $x_1+x_2\geq 7$ ，可以增减去一个非负松弛变量 $x_3$ ，变为： $x_1+x_2-x_3=7,x_3\geq 0$

求 $x_1+x_2$ 的最小值，可以取反，变为求 $x_1-x_2$ 的最大值。

无约束的变量 $-\infty\leq x_1\leq\infty$ ，转化为两个变量的差： $x_1=x_2-x_3,x_2\geq 0,x_3\geq 0$

求解

先找出 $m * n$ 矩阵的基（行列式非0的 $m * m$ 子矩阵），令其它变量为0，若解出的基解满足约束条件（ $\geq 0$ ），则为基可行解。有定理：最优解一定是某个基可行解。

代码

original link - http://uoj.ac/problem/179

#include <bits/stdc++.h>
#define fo(i,a,b) for(int i=a;i<=b;++i)
#define fod(i,a,b) for(int i=a;i>=b;--i)
#define N 50
using namespace std;
const double eps=1e-9;
double a[N][N],ans[N];
int n,m,t,id[N],d[N];
int rd(int k) {
    return rand()%k;
}
void pivot(int l,int e) {
    swap(id[l+n],id[e]);
    double t=a[l][e];
    a[l][e]=1;
    fo(i,0,n) a[l][i]/=t;
    fo(i,0,m) {
        if(i!=l&&abs(a[i][e])>eps) {
            double t=a[i][e];
            a[i][e]=0;
            fo(j,0,n) a[i][j]-=t*a[l][j];
        }
    }
}
bool prp() {
    while(1) {
        int l=0,e=0;
        d[0]=0;
        fo(i,1,m) if(a[i][0]<-eps)
            d[++d[0]]=i;
        if(!d[0])
            return 1;
        l=d[rd(d[0])+1],d[0]=0;
        fo(i,1,n) if(a[l][i]<-eps)
            d[++d[0]]=i;
        if(!d[0])
            return 0;
        e=d[rd(d[0])+1];
        pivot(l,e);
    }
}
bool solve() {
    while(1) {
        int l=0,e=0;
        double mi=1e18;
        fo(i,1,n)
        if(a[0][i]>eps&&(!e||a[0][i]>a[0][e]))
            e=i;
        if(!e)
            return 1;
        fo(i,1,m) {
            if(a[i][e]>eps&&a[i][0]/a[i][e]<mi)
                mi=a[i][0]/a[i][e],l=i;
        }
        if(!l)
            return 0;
        pivot(l,e);
    }
}
int main() {
    srand(time(0));
    cin>>n>>m>>t; // n变量 m约束
    fo(i,1,n) {
        int x;
        scanf("%d",&x);
        a[0][i]=x; // 目标方程
    }
    fo(i,1,m) { // 约束 aij xj <= bi
        int x;
        fo(j,1,n) scanf("%d",&x),a[i][j]=x;
        scanf("%d",&x),a[i][0]=x;
    }
    fo(i,1,n+m) id[i]=i;
    if(!prp()) { // 无可行域
        printf("Infeasible\n");
        return 0;
    }
    if(!solve()) { // 无界，取得inf
        printf("Unbounded\n");
        return 0;
    }
    printf("%.10lf\n",-a[0][0]);
    if(t==1) {
        fo(i,1,m) ans[id[i+n]]=a[i][0];
        fo(i,1,n) printf("%.10lf ",ans[i]);
    }
}