前言:这些是一些大神的博客
https://www.cnblogs.com/ECJTUACM-873284962/p/7097864.html
https://ldxcaicai.github.io/LP/
线性规划
标准型:
max(∑i=1ncixi)max(\sum_{i=1}^n c_ix_i)max(i=1∑ncixi)
limit:∑j=1nai,j∗xj≤bi, x≥0limit:\sum_{j=1}^na_{i,j}*x_j\le b_i,\ x\ge 0limit:j=1∑nai,j∗xj≤bi, x≥0
松弛型:
max(∑i=1ncixi)max(\sum_{i=1}^n c_ix_i)max(i=1∑ncixi)
limit:∑j=1nai,j∗xj=bi, x≥0limit:\sum_{j=1}^na_{i,j}*x_j= b_i,\ x\ge 0limit:j=1∑nai,j∗xj=bi, x≥0
标准型与松弛型的转换:
令 xi+m=bi−∑j=1nai,j∗xjx_{i+m}=b_i-\sum_{j=1}^na_{i,j}*x_jxi+m=bi−∑j=1nai,j∗xj,那么有 ∑j=1nai,j∗xj+xi+m=bi\sum_{j=1}^na_{i,j}*x_j+x_{i+m}= b_i∑j=1nai,j∗xj+xi+m=bi
单纯形求解线性规划:
一直重复以下操作:
- 在目标函数中找到一个 ce≥0c_e\ge 0ce≥0 的 eee
- 在约束条件中找到一个al,e>0a_{l,e}>0al,e>0且使blal,e\frac{b_l}{a_{l,e}}al,ebl最大的 lll
- 将第 lll 个等式进行移项,xe=b[l]−∑j!=ea[l][j]−xl+ma[l][e]x_e=\frac{b[l]-\sum_{j!=e}a[l][j]-x_{l+m}}{a[l][e]}xe=a[l][e]b[l]−∑j!=ea[l][j]−xl+m
- 带到其他行,消去 xex_exe,新代入 xl+mx_{l+m}xl+m
- 当无法找到 eee 时,说明已经最优
- 当无法找到合法的 lll 时,说明该线性规划无界
令:当存在 b[i]<0b[i]<0b[i]<0 的时候可能会出锅,可以选择对 b[i]<0b[i]<0b[i]<0 的列随机扰动
具体实现参考代码
满足线性规划问题约束条件的所有点组成的集合就是线性规划的可行域。若可行域有界(以下主要考虑有界可行域),线性规划问题的目标函数最优解必然在可行域的顶点上达到最优。
单纯形法就是通过设置不同的基向量,经过矩阵的线性变换,求得基可行解(可行域顶点),并判断该解是否最优,否则继续设置另一组基向量,重复执行以上步骤,直到找到最优解。所以,单纯形法的求解过程是一个循环迭代的过程
对偶原理:
最大化和最小化可以通过正负来转换,而大于小于号则不行
考虑如下的一个线性规划:
min(∑i=1ncixi)min(\sum_{i=1}^n c_ix_i)min(i=1∑ncixi)
limit:∑j=1nai,j∗xj≥bi, x≥0limit:\sum_{j=1}^na_{i,j}*x_j\ge b_i,\ x\ge 0limit:j=1∑nai,j∗xj≥bi, x≥0
利用待定系数法的思想可以知道,该线性规划的任意可行解 {x1,x2,...,xn}\{x_1,x_2,...,x_n\}{x1,x2,...,xn} 一定能表示成如下形式:
∑i=1ncixi≥∑i=1myi(∑j=1nai,j∗xj)≥∑i=1myi∗bi\sum_{i=1}^nc_ix_i\ge \sum_{i=1}^my_i(\sum_{j=1}^na_{i,j}*x_j)\ge\sum_{i=1}^my_i*b_ii=1∑ncixi≥i=1∑myi(j=1∑nai,j∗xj)≥i=1∑myi∗bi
∑i=1myi∗aj,i≤cj,yi≥0\sum_{i=1}^my_i*a_{j,i}\le c_j,y_i\ge 0i=1∑myi∗aj,i≤cj,yi≥0
于是可以转换为一个新的线性规划
找到 ∑i=1myi∗bi\sum_{i=1}^my_i*b_i∑i=1myi∗bi 的上界就相当于找到 ∑i=1nci∗xi\sum_{i=1}^nc_i*x_i∑i=1nci∗xi
但是我们还要证明当 ∑i=1myi∗bi\sum_{i=1}^my_i*b_i∑i=1myi∗bi 取到最优解即最大值时,∑i=1nci∗xi\sum_{i=1}^nc_i*x_i∑i=1nci∗xi 也取到最优解即最小值
证明:
可以参考
https://blog.youkuaiyun.com/itnerd/article/details/83352441
将线性规划表示成矩阵的形式就是要证这么一个东西:
min{cTx∣Ax≥b,x≥0}⇔max{bTy∣ATy≤c,y≥0}min\{c^Tx|Ax\ge b,x\ge 0\}\Leftrightarrow \max\{b^Ty|A^Ty\le c,y\ge 0\}min{cTx∣Ax≥b,x≥0}⇔max{bTy∣ATy≤c,y≥0}
对于第一个线性规划,构造一个函数 f(y)=cTx+yT(b−Ax),y≥0f(y)=c^Tx+y^T(b-Ax),y\ge 0f(y)=cTx+yT(b−Ax),y≥0
令 g(y)=min{f(y)∣Ax≥b,x≥0}g(y)=min\{f(y)|Ax\ge b,x\ge0\}g(y)=min{f(y)∣Ax≥b,x≥0}
∀g(y)≤cTx\forall g(y)\le c^Tx∀g(y)≤cTx, max(g(y))=cTxmax(g(y))=c^T xmax(g(y))=cTx
于是有
min(cTx)=max(g(y))min(c^Tx)=max(g(y))min(cTx)=max(g(y))
=max(yTb+min((cT−yTA)x))=max(y^Tb+min((c^T-y^TA)x))=max(yTb+min((cT−yTA)x))
=max(bTy+min((cT−yTA)x))=max(b^Ty+min((c^T-y^TA)x))=max(bTy+min((cT−yTA)x))
当 cT−yTA≥0c^T-y^TA\ge 0cT−yTA≥0,那么 min((cT−yTA)x)=0min((c^T-y^TA)x)=0min((cT−yTA)x)=0
反之 min((cT−yTA)x)=−∞min((c^T-y^TA)x)=-\inftymin((cT−yTA)x)=−∞
所以
=max(bTy+min((cT−yTA)x))=max(bTy)=max(b^Ty+min((c^T-y^TA)x))=max(b^Ty)=max(bTy+min((cT−yTA)x))=max(bTy)
令:翻到一个感性的证明,有点像经济学上买家卖家的博弈
线性规划互相松弛定理
这里给出 NOI 2008 志愿者招募的题解及代码
令 ai,ja_{i,j}ai,j 表示第 iii 类志愿者在第 jjj 天是否工作,xix_ixi 表示第 iii 类志愿者的个数
那么有线性规划
min(∑i=1mxi∗ci)min(\sum_{i=1}^mx_i*c_i)min(i=1∑mxi∗ci)
limit:∑j=1naj,i∗xj≥bi, xi≥0limit:\sum_{j=1}^na_{j,i}*x_j\ge b_i,\ x_i\ge 0limit:j=1∑naj,i∗xj≥bi, xi≥0
对偶过后发现
max(∑i=1nyi∗bi)max(\sum_{i=1}^ny_i*b_i)max(i=1∑nyi∗bi)
limit:∑j=1mai,j∗yj≤cilimit:\sum_{j=1}^ma_{i,j}*y_j\le c_ilimit:j=1∑mai,j∗yj≤ci
已经是标准形式了
#include<bits/stdc++.h>
#define cs const
using namespace std;
int read(){
int cnt = 0, f = 1; char ch = 0;
while(!isdigit(ch)){ ch = getchar(); if(ch == '-') f = -1; }
while(isdigit(ch)) cnt = cnt*10 + (ch-'0'), ch = getchar();
return cnt * f;
}
cs int N = 1e3 + 5, M = 1e4 + 5;
cs double eps = 1e-8, INF = 1e15;
int n, m; double a[M][N];
void pivot(int l, int e){
double t = a[l][e]; a[l][e] = 1;
for(int i = 0; i <= n; i++) a[l][i] /= t;
for(int i = 0; i <= m; i++) if(l!=i && fabs(a[i][e]) > eps){
t = a[i][e]; a[i][e] = 0;
for(int j = 0; j <= n; j++) a[i][j] -= a[l][j] * t;
}
}
bool simplex(){
while(true){
int l = 0, e = 0; double mi = INF;
for(int i = 1; i <= n; i++) if(a[0][i] > eps) { e = i; break; } if(!e) break;
for(int i = 1; i <= m; i++) if(a[i][e] > eps && a[i][0] / a[i][e] < mi) mi = a[i][0] / a[i][e], l = i;
if(!l) return false; pivot(l, e);
} return true;
}
int main(){
n = read(), m = read();
for(int i = 1; i <= n; i++) a[0][i] = read();
for(int i = 1; i <= m; i++){
int l = read(), r = read(), c = read();
for(int j = l; j <= r; j++) a[i][j] = 1;
a[i][0] = c;
} if(simplex()) cout << (long long)(-a[0][0] + 0.5);
return 0;
}
附:随机扰动的代码
bool init(){
while(true){
int l = 0, e = 0;
for(int i = 1; i <= m; i++)
if(a[i][0] < -eps && (!l || (rand()&1))) l = i;
if(!l) break;
for(int j = 1; j <= n; j++)
if(a[l][j] < -eps && (!e || (rand()&1))) e = j;
if(!e) return false; pivot(l, e);
} return true;
}