线性规划——单纯形法学习笔记

发布于 2020-12-12  34 次阅读



首先把板子放在这儿......不过这个板子好像还是有点不太对,后面有时间再改改吧(

简化版:

namespace Simplex //单纯形法,求得目标函数的max值
{
    const int N = 150, M = 150; //N是变量个数,M是条件限制个数
    const double INF = 0x3f3f3f3f3f3f3f3f;
    const double eps = 1e-7; //精度误差
    int n, m;
    double a[M][N], b[M], c[N], v;
    //max cX
    //s.t aX<=b
    void pivot(int l, int e)
    {
        b[l] /= a[l][e];
        for (int j = 1; j <= n; j++)
            if (j != e)
                a[l][j] /= a[l][e];
        a[l][e] = 1 / a[l][e];
        for (int i = 1; i <= m; i++)
            if (i != l && fabs(a[i][e]) > 0)
            {
                b[i] -= a[i][e] * b[l];
                for (int j = 1; j <= n; j++)
                    if (j != e)
                        a[i][j] -= a[i][e] * a[l][j];
                a[i][e] = -a[i][e] * a[l][e];
            }
        v += c[e] * b[l];
        for (int j = 1; j <= n; j++)
            if (j != e)
                c[j] -= c[e] * a[l][j];
        c[e] = -c[e] * a[l][e];
    }
    double simplex()
    {
        while (true)
        {
            int e = 0, l = 0;
            for (e = 1; e <= n; e++)
                if (c[e] > eps)
                    break;
            if (e == n + 1)
                return v;
            double mn = INF;
            for (int i = 1; i <= m; i++)
                if (a[i][e] > eps && mn > b[i] / a[i][e])
                    mn = b[i] / a[i][e], l = i;
            if (mn == INF)
                return INF; // unbounded
            pivot(l, e);
        }
    }
} // namespace Simplex

几乎完全版

namespace Simplex
{
    const int N = 150, M = 150; //N是变量个数,M是条件限制个数
    const double INF = 3e9;
    const double eps = 1e-7; //精度误差
    int n, m;
    double a[M][N], b[M], c[N], v, A[M][N];
    int idx[N], idy[N];
    void pivot(const int &e, const int &l)
    {
        double tmp = -A[l][e];
        for (int i = 0; i <= n; i++)
            if (i != e)
                A[l][i] /= tmp;
        A[l][e] = -1 / tmp;
        for (int i = 0; i <= m; i++)
        {
            if (i == l)
                continue;
            for (int j = 0; j <= n; j++)
            {
                if (j == e)
                    continue;
                A[i][j] += A[i][e] * A[l][j];
            }
            A[i][e] *= A[l][e];
        }
        swap(idx[e], idy[l]);
    }
    double ans[N];
    double simplex()
    {
        for (int i = 1; i <= n; ++i)
            A[0][i] = c[i];
        for (int i = 1; i <= m; ++i)
        {
            for (int j = 1; j <= n; ++j)
                A[i][j] = a[i][j];
            A[i][0] = b[i];
        }
        for (int i = 1; i <= n; i++)
            idx[i] = i;
        for (int i = 1; i <= m; i++)
            idy[i] = n + i;
        int x, y;
        while (1)
        {
            x = y = 0;
            for (int i = 1; i <= m; i++)
                if (A[i][0] < -eps)
                {
                    y = i;
                    if (rand()
                        break;
                }
            if (!y)
                break;
            for (int i = 1; i <= n; i++)
                if (A[y][i] > eps)
                {
                    x = i;
                    if (rand()
                        break;
                }
            if (!x)
                return -INF;
            pivot(x, y);
        }
        while (1)
        {
            x = y = 0;
            for (int i = 1; i <= n; i++)
                if (A[0][i] > eps)
                {
                    x = i;
                    break;
                }
            if (!x)
                break;
            double temp = INF;
            for (int i = 1; i <= m; i++)
                if (A[i][x] < -eps)
                {
                    double nt = A[i][0] / -A[i][x];
                    if (nt < temp)
                        temp = nt, y = i;
                }
            if (!y)
                return INF;
            pivot(x, y);
        }
        // printf(
        return A[0][0];
        // if (!t)
        //     return;
        // for (int i = 1; i <= m; i++)
        //     if (idy[i] <= n)
        //         ans[idy[i]] = A[i][0];
        // for (int i = 1; i <= n; i++)
        //     printf(
    }
} // namespace Simplex

怎么用呢?这个板子实际上求解的是:
$$
\max\space\space c^TX \\
s.t\space\space AX\le b
$$
如果要求 $\min$ 的话,转一下对偶即可:
$$
\max\space\space b^TY \\
s.t\space\space A^TY\ge c
$$
总之要注意所有的条件一个都不能漏!
然后接下来开始线性规划的学习吧~


线性规划:将不等式约束通过引入松弛变量和剩余变量转化为等式约束,后者称为标准型

最优解一定在多面体端点上,检查所有端点的复杂度为指数级别。


松弛变量:$\sum_{j=1}^na_{pj}x_j+x_{n+p}=b_p(x_{n+p\ge 0})$

剩余变量:$\sum_{j=1}^na_{pj}x_j-x_{n+p}=b_p(x_{n+p\ge 0})$

自由变量:$x_i=x_i^{+}-x_i^{-}(x_i^{+},x_i^{-}\ge 0, x_i\in R)$

将目标函数求 max 转为求 min

$$
AX=b(b\ge 0) \
A=[B,N],B可逆(A的基向量组)\
基本解:x=(B^{-1}b,0)^T\ge0\space\space(令N=0)
$$
基可行解:既是基本解,又是可行解(可行解:满足方程组变量条件,$X\ge 0$)

不同的基最多有 $\binom{n}{m}$ 个,而一个基最多对应一个基可行解。


顶点(极点)

设 $X^0\in D\space R^n,D是凸集$,$X^0$ 不能表示为 $D$ 中其他任意两不同点的凸组合,则 $X^0$ 为 $D$ 的顶点(极点)。

在 $D^2$ 中,多边形的顶点就是那些端点。半圆的话,圆弧上的点都是顶点。

或者说这个点不在 $D$ 内的任意一条线段上。

设 $X$ 是线性规划的可行解,则 $X$ 为基可行解充要条件为:$X$ 的非零向量在 $A$ 中对应的列向量组线性无关。

一组基最多对应一个顶点。

线性规划、凸优化的两个最优解的凸组合也是最优解(凹优化也是)


单纯形法

从多面体某一个顶点出发,沿着目标函数值下降的方向寻找下一个顶点(沿着端点挨着跳,一边跳一边检查)。

一般来说 $AX=B$ 好解,$AX\le B$ 很困难。

准备工作

最优解判别准则:何时迭代终止,找到最优解

$$
\min\space f(X)=C^TX \
s.t.\space AX=b \
X\ge 0
$$

设 $A=[I\space\space N]$ 转化为
$$
\min\space f(X)=C_I^TX_I+C_N^TX_N \
s.t.\space [I\space\space N][X_I\space\space X_N]^T=b \
X_I,X_N\ge 0 \
I\space 是单位矩阵
$$
$\bar{x}=\binom{b}{0}\ge 0\space\space(基本可行解,即顶点),\space 则\space f(\bar{x})=C_{I}^Tb\space 为最小值的充要条件为\space C_I^TN-C_N^T\le 0$

证明:$f(x)=C_I^TX_I+C_N^TX_N=C_I^T(b-NX_N)+C_N^TX_N=f(\bar{x})-[C_I^TN-C_N^T]X_N\ge f(\bar{x}) \space (C_I^TN-C_N^T\le 0)$

通常叫 $\space C_I^TN-C_N^T$ 为判别式。

若判别式中某个判别数 $\sigma_j>0$ 但相应的列向量 $P_j\le 0$,则该线性规划无最优解(即没有下界)。

换基运算:从一个顶点转到另一个顶点

利用不交换行的行初等变换,选定其中一个位置(该位置不等于 $0$),该位置变为 $1$,这一列其他位置变为 $0$.

选择 $\sigma_j$ 最大的一列来换,这样函数值下降最多。$b_k'=\frac{b_k}{a_{kl}},\space b_i'=b_i-b_k'a_{il}$,故应选择这一列的正数 $a_{kl}$ 且使得 $\frac{b_k}{a_{kl}}$ 最小。$\frac{b_k}{a_{kl}}$ 称为转轴元素。这一列称为进基列。

单纯形表

用分块矩阵的运算来解释,进行行初等变换即可。
$$
\begin{bmatrix}
A&b \\
-C^T&0
\end{bmatrix}
\to
\begin{bmatrix}
I&N &b \\
-C_I^T&-C_N^T &0
\end{bmatrix}
\to
\begin{bmatrix}
I&N &b \\
0&C_I^TN-C_N^T &C_I^Tb
\end{bmatrix}
$$
最后一个称为初始单纯形表。如下:
$$
\begin{bmatrix}
P_1& P_2&...&P_m&P_{m+1}&...&P_l&...&P_n&b \\
1 & & & &a_{1,m+1}&...&a_{1l}&...&a_{1n}&b_1 \\
...\\
0&0&...&0&\sigma_{m+1}&...&0&...&\sigma_n&f
\end{bmatrix}
$$


两阶段方法

如果 $A=[B\space\space N]$,即不显含单位矩阵 $I$,则构造人工变量 $Y=(y_1,y_2,...,y_m)^T$,有
$$
\min\space\space \sum y_i \\
s.t\space\space Y+AX=b
$$
注意到含 $y_i$ 的系数都是 $1$,此时有一个显含的单位矩阵 $I$,进行 $m$ 次行初等变换之后就变成一个初始单纯形表。

上述问题为辅助线性规划,求得最优解。

若 $\sum Y>0$,则原问题无可行解。

证明:若有可行解,则 $\binom{0}{\bar{x}}$ 一定是辅助线性规划的一个特解(因为这是原问题的一个特解)。此时一定有 $\sum Y=0$。 证毕。

若 $\sum Y=0$,则 $y_1=y_2=...=0$。这是我们希望的。此时约束条件就是原问题的约束条件。

转轴之后单位矩阵的列会转到所有 $x_i$ 列上,此时就找到一个极大无关组,即 $B$。

这种做法比直接行初等变换找极大无关组速度更快。


对偶问题

$$
\max\space\space c^TX \\
s.t\space\space AX\le b
$$
这类问题可以用对偶原则转化为:
$$
\min\space\space b^TY \\
s.t.\space\space A^TY\ge c
$$
证明:写出拉格朗日函数:$L(x,\lambda,y)=c^TX+\lambda^T(AX-b)-\mu^TX $

$q(\lambda,\mu)=\min\limits_{X\in \mathcal{R^{n}}} L(x,\lambda,\mu)=\min\limits_{X\in \mathcal{R^n}}{(c^T-\lambda^TA-\mu^T)X}+\lambda^Tb=\lambda^Tb+ \min\limits_{X\in \mathcal{R^n}}(c-A^T\lambda-\mu)^TX=\lambda^Tb\space\space(c-A^T\lambda-\mu=0) $

故对偶问题可以表示为:
$$
\max\space\space \lambda^Tb \\
s.t\space\space c-A^T\lambda \ge0
$$


至于线性规划如何转化成网络流问题,这个坑先留在这里吧(


下面是一些题目

方格取数问题

这个题用费用流解蛮经典的了,但是实际上用线性规划来表示的话也很直观:
$$
\max\space\space \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m} c_{ij}x_{ij} \\
s.t\space\space x_{ij}+x_{i(j+1)}\le 1 \\
x_{ij}+x_{(i+1)j}\le 1
$$
一共 $n\cdot m$ 个变量,$(n-1)(m-1)$ 个约束,实际上直接求解即可......可能是因为单纯形法的复杂度一直都很玄学吧(

#include <bits/stdc++.h>
#define ll long long
#define ls id << 1
#define rs id << 1 | 1
#define mem(array, value, size, type) memset(array, value, ((size) + 5) * sizeof(type))
#define memarray(array, value) memset(array, value, sizeof(array))
#define fillarray(array, value, begin, end) fill((array) + (begin), (array) + (end) + 1, value)
#define fillvector(v, value) fill((v).begin(), (v).end(), value)
#define pb(x) push_back(x)
#define st(x) (1LL << (x))
#define pii pair<int, int>
#define mp(a, b) make_pair((a), (b))
#define Flush fflush(stdout)
#define vecfirst (*vec.begin())
#define veclast (*vec.rbegin())
#define vecall(v) (v).begin(), (v).end()
#define vecupsort(v) (sort((v).begin(), (v).end()))
#define vecdownsort(v, type) (sort(vecall(v), greater<type>()))
#define veccmpsort(v, cmp) (sort(vecall(v), cmp))
using namespace std;
const int N = 150;
const int inf = 0x3f3f3f3f;
const ll llinf = 0x3f3f3f3f3f3f3f3f;
const int mod = 998244353;
const int MOD = 1e9 + 7;
const double PI = acos(-1.0);
clock_t TIME__START, TIME__END;
void program_end()
{
#ifdef ONLINE
    printf("\n\nTime used:
    system("pause");
#endif
}
namespace Simplex //单纯形法,求得目标函数的max值
{
    const int N = 10001, M = 9900; //N是变量个数,M是条件限制个数
    const double INF = 2e9;
    const double eps = 1e-9; //精度误差
    int n, m;
    double a[M][N], b[M], c[N], v;
    //max cX
    //s.t aX<=b
    void pivot(int l, int e)
    {
        b[l] /= a[l][e];
        for (int j = 1; j <= n; j++)
            if (j != e)
                a[l][j] /= a[l][e];
        a[l][e] = 1 / a[l][e];
        for (int i = 1; i <= m; i++)
            if (i != l && fabs(a[i][e]) > 0)
            {
                b[i] -= a[i][e] * b[l];
                for (int j = 1; j <= n; j++)
                    if (j != e)
                        a[i][j] -= a[i][e] * a[l][j];
                a[i][e] = -a[i][e] * a[l][e];
            }
        v += c[e] * b[l];
        for (int j = 1; j <= n; j++)
            if (j != e)
                c[j] -= c[e] * a[l][j];
        c[e] = -c[e] * a[l][e];
    }
    double simplex()
    {
        while (true)
        {
            int e = 0, l = 0;
            for (e = 1; e <= n; e++)
                if (c[e] > eps)
                    break;
            if (e == n + 1)
                return v;
            double mn = INF;
            for (int i = 1; i <= m; i++)
                if (a[i][e] > eps && mn * a[i][e] > b[i])
                    mn = b[i] / a[i][e], l = i;
            if (mn == INF)
                return INF; // unbounded
            pivot(l, e);
        }
    }
} // namespace Simplex
int n, m, a[N][N];

inline void solve()
{
    scanf(
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= m; ++j)
            scanf(
    Simplex ::n = n * m;
    for (int i = 1; i <= n; ++i)
        for (int j = 1; j <= m; ++j)
            Simplex::c[(i - 1) * m + j] = a[i][j];
    int k = 0;
    for (int i = 1; i <= n; ++i)
    {
        for (int j = 1; j < m; ++j)
        {
            k++;
            Simplex::b[k] = 1;
            Simplex::a[k][(i - 1) * m + j] = 1;
            Simplex::a[k][(i - 1) * m + j + 1] = 1;
        }
    }
    for (int j = 1; j <= m; ++j)
    {
        for (int i = 1; i < n; ++i)
        {
            k++;
            Simplex::b[k] = 1;
            Simplex::a[k][(i - 1) * m + j] = 1;
            Simplex::a[k][i * m + j] = 1;
        }
    }
    Simplex::m = k;
    int ans = Simplex::simplex() + 0.5;
    printf(
}

int main()
{
    TIME__START = clock();
    int Test = 1;
    // scanf(
    while (Test--)
    {
        solve();
        // if (Test)
        //     putchar('\n');
    }
    TIME__END = clock();
    program_end();
    return 0;
}

Codeforces 1061E - Politics

这题时限 4s,我愣是 3650ms 卡过去了呃......

反正线性规划式子也很好推,注意要限制每一个点最多被选一次,还有判断无解的情况。

然后直接按照题意来设限制即可,复杂度玄学(

#include <bits/stdc++.h>
#define ll long long
#define ls id << 1
#define rs id << 1 | 1
#define mem(array, value, size, type) memset(array, value, ((size) + 5) * sizeof(type))
#define memarray(array, value) memset(array, value, sizeof(array))
#define fillarray(array, value, begin, end) fill((array) + (begin), (array) + (end) + 1, value)
#define fillvector(v, value) fill((v).begin(), (v).end(), value)
#define pb(x) push_back(x)
#define st(x) (1LL << (x))
#define pii pair<int, int>
#define mp(a, b) make_pair((a), (b))
#define Flush fflush(stdout)
#define vecfirst (*vec.begin())
#define veclast (*vec.rbegin())
#define vecall(v) (v).begin(), (v).end()
#define vecupsort(v) (sort((v).begin(), (v).end()))
#define vecdownsort(v, type) (sort(vecall(v), greater<type>()))
#define veccmpsort(v, cmp) (sort(vecall(v), cmp))
using namespace std;
const int N = 505;
const int inf = 0x3f3f3f3f;
const ll llinf = 0x3f3f3f3f3f3f3f3f;
const int mod = 998244353;
const int MOD = 1e9 + 7;
const double PI = acos(-1.0);
clock_t TIME__START, TIME__END;
void program_end()
{
#ifdef ONLINE
    printf("\n\nTime used:
    system("pause");
#endif
}
namespace Simplex
{
    const int N = 505, M = 2550; //N是变量个数,M是条件限制个数
    const double INF = 3e9;
    const double eps = 1e-6; //精度误差
    int n, m;
    double a[M][N], b[M], c[N], v, A[M][N];
    int idx[N], idy[M];
    void pivot(const int &e, const int &l)
    {
        double tmp = -A[l][e];
        for (int i = 0; i <= n; i++)
            if (i != e)
                A[l][i] /= tmp;
        A[l][e] = -1 / tmp;
        for (int i = 0; i <= m; i++)
        {
            if (i == l)
                continue;
            for (int j = 0; j <= n; j++)
            {
                if (j == e)
                    continue;
                A[i][j] += A[i][e] * A[l][j];
            }
            A[i][e] *= A[l][e];
        }
        swap(idx[e], idy[l]);
    }
    double ans[N];
    double simplex()
    {
        srand(1);
        for (int i = 1; i <= n; ++i)
            A[0][i] = c[i];
        for (int i = 1; i <= m; ++i)
        {
            for (int j = 1; j <= n; ++j)
                A[i][j] = -a[i][j];
            A[i][0] = b[i];
        }
        for (int i = 1; i <= n; i++)
            idx[i] = i;
        for (int i = 1; i <= m; i++)
            idy[i] = n + i;
        int x, y;
        while (1)
        {
            x = y = 0;
            for (int i = 1; i <= m; i++)
                if (A[i][0] < -eps)
                {
                    y = i;
                    if (rand()
                        break;
                }
            if (!y)
                break;
            for (int i = 1; i <= n; i++)
                if (A[y][i] > eps)
                {
                    x = i;
                    if (rand()
                        break;
                }
            if (!x)
                return -INF;
            pivot(x, y);
        }
        while (1)
        {
            x = y = 0;
            for (int i = 1; i <= n; i++)
                if (A[0][i] > eps)
                {
                    x = i;
                    break;
                }
            if (!x)
                break;
            double temp = INF;
            for (int i = 1; i <= m; i++)
                if (A[i][x] < -eps)
                {
                    double nt = A[i][0] / -A[i][x];
                    if (nt < temp)
                        temp = nt, y = i;
                }
            if (!y)
                return INF;
            pivot(x, y);
        }
        return A[0][0];
    }
} // namespace Simplex
int n, rx, ry;
vector<int> e1[N], e2[N];
int a[N];
int depx[N], depy[N];
void dfsdepx(int u, int fa)
{
    depx[u] = depx[fa] + 1;
    for (auto v : e1[u])
    {
        if (v == fa)
            continue;
        dfsdepx(v, u);
    }
}
void dfsdepy(int u, int fa)
{
    depy[u] = depy[fa] + 1;
    for (auto v : e2[u])
    {
        if (v == fa)
            continue;
        dfsdepy(v, u);
    }
}
void dfsx(int u, int fa, int k, int f)
{
    Simplex::a[k][u] = f;
    for (auto v : e1[u])
    {
        if (v == fa || depx[v] < depx[u])
            continue;
        dfsx(v, u, k, f);
    }
}
void dfsy(int u, int fa, int k, int f)
{
    Simplex::a[k][u] = f;
    for (auto v : e2[u])
    {
        if (v == fa || depy[v] < depy[u])
            continue;
        dfsy(v, u, k, f);
    }
}

inline void solve()
{
    scanf(
    Simplex::n = n;
    for (int i = 1; i <= n; ++i)
    {
        scanf(
        Simplex::c[i] = a[i];
    }
    for (int i = 1; i < n; ++i)
    {
        int x, y;
        scanf(
        e1[x].push_back(y);
        e1[y].push_back(x);
    }
    for (int i = 1; i < n; ++i)
    {
        int x, y;
        scanf(
        e2[x].push_back(y);
        e2[y].push_back(x);
    }
    dfsdepx(rx, 0);
    dfsdepy(ry, 0);
    int K = 0;
    for (int i = 1; i <= n; ++i)
    {
        K++;
        Simplex::b[K] = 1;
        Simplex::a[K][i] = 1;
    }
    int qx, qy;
    scanf(
    while (qx--)
    {
        int k, x;
        scanf(
        K++;
        Simplex::b[K] = x;
        dfsx(k, 0, K, 1);
        K++;
        Simplex::b[K] = -x;
        dfsx(k, 0, K, -1);
    }
    scanf(
    while (qy--)
    {
        int k, x;
        scanf(
        K++;
        Simplex::b[K] = x;
        dfsy(k, 0, K, 1);
        K++;
        Simplex::b[K] = -x;
        dfsy(k, 0, K, -1);
    }
    Simplex::m = K;
    double ans = Simplex::simplex();
    if (ans == -Simplex::INF)
        return puts("-1"), void();
    // printf(
    int ans2 = ans + 0.5;
    cout << ans2 << endl;
}

int main()
{
    TIME__START = clock();
    int Test = 1;
    // scanf(
    while (Test--)
    {
        solve();
        // if (Test)
        //     putchar('\n');
    }
    TIME__END = clock();
    program_end();
    return 0;
}