From c4db78b9251bc5f97faaa0ff39c586684d99e47a Mon Sep 17 00:00:00 2001 From: huhaoo Date: Sun, 31 Mar 2019 15:40:58 +0800 Subject: [PATCH] Update linear-programming.md --- docs/math/linear-programming.md | 174 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 174 insertions(+) diff --git a/docs/math/linear-programming.md b/docs/math/linear-programming.md index 66873dee..7dda7608 100644 --- a/docs/math/linear-programming.md +++ b/docs/math/linear-programming.md @@ -106,3 +106,177 @@ 显而易见的，打了蓝色斜线的区域为三块区域的交集，这就是这个线性规划的所有可行解。因为题目中说明，需要最小化 x 和 y，观察图像可知，点（2.3）为可行解中 x 和 y 最小的一个。因此， $x_{min}=2,y_{min}=3$ 把求解线性规划的图解法总结起来，就是先在坐标系中作出所有的约束条件，然后作出需要求极值的线性函数的定义域。定义域与约束条件的交集就是这个线性规划的解集，而所需求的极值由观察可以轻易得出。 + +### 单纯形法 + +> $$ +> \forall k,\sum_{i}a_{k,i}x_i\le b_{k}\\ +> \forall i,x_i\ge 0 +> $$ +> +> 求 +> $$ +> \max(\sum_{i}c_{i}x_i) +> $$ + +我们考虑一个下标从$0$开始的矩阵$A_{0\dots m,0\dots n+m}$ +$$ +A= +\begin{bmatrix} +0&c_1&\ldots&c_n&0&0&\ldots&0\\ +b_1&a_{1,1}&\ldots&a_{1,n}&1&0&\dots&0\\ +b_2&a_{2,1}&\ldots&a_{2,n}&0&1&\dots&0\\ +\vdots&\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\ +b_m&a_{m,1}&\ldots&a_{m,n}&0&0&\dots&1 +\end{bmatrix} +$$ +其中$A_{0,0}$为这个状态和初始时答案之差，新增的$m$行相当于新建$m$个变量，编号为$x_{n+1\dots n+m}$，要求它们不小于$0$ + +于是 +$$ +\forall i\in[1,m]\\ +\sum_{1\le j\le n+m}a_{i,j}x_j=b_i +$$ +显然，要先求出一套可行的方案 + +如果 +$$ +x_i=\cases{0&$i\in[1,n]$\\b_{i-n}&$i\in[n+1,n+m]$} +$$ +那么显然满足，但是由于$b_i$不一定不小于$0$，所以这样可能导致$x_{i}$小于$0$ + +然后，介绍一个魔法，通过矩阵变化改变$b$的值来求合法解/最优解 + +#### 转轴 + +~~之所以叫这个名称，是因为可行域是一个超立方体，这个魔法是经过一条轴，但是我不讲~~ + +这个魔法是指： + +选择两个数$u,v$，考虑将$x_v$变为非$0$数而$x_{u+n}$变为$0$ + +也就是让$x_v$来代替$x_{u+n}$作为等于$b_u$的位置 + +如果要实现，就要运用矩阵变换 + +我们可以先让第$u$行全部除一个$A_{u,v}$ + +然后把其它所有行$i$的位置$v$变为$0$（即使$A_{i,v}=0$），可以直接把第$i$行减去$A_{i,v}$倍第$u$行 + +然后就可以了！~~（简单吧）~~ + +因为（不算$u+n,v$列的话）后面$m$列只有$A_{i,i+n}=1$，前面$n$列都是$0$，所以矩阵变化不会影响其它的$x$ + +另外，这个会影响 + +当然，最好把第$u+n,v​$行换一下，这样会方便很多 + +（还有，因为矩阵后$m$列值的特殊性，一般不存这$m$列） + +```cpp +int id[N<<1]; +double a[N][N]; +void pivot(int u,int v) +{ + double k=a[u][v]; + swap(id[u+n],id[v]); + a[u][v]=1.;//已经知道第v列要被消掉了，就可以直接变为u+n列了（尽管要用一遍，就临时存一下） + for(int i=0;i<=n;i++) + a[u][i]/=k; + for(int i=0;i<=m;i++) + if(i!=u&&fabs(a[i][v])>eps) + { + k=a[i][v]; + a[i][v]=0;//同上 + for(int j=0;j<=n;j++) + a[i][j]-=k*a[u][i]; + } +} +``` + +知道了基本的变（魔）化（法），就来看怎么求一个合法解 + +#### 初始化 + +引用一下Candy?的博客 + +> 算法导论上有一个辅助线性规划的做法 +> +> 但我发现好多人都用了**随机初始化**的黑科技 +> +> 在所有$a_{i,0}<0$的约束中随机选一个作为$u$，再随机选一个$a_{u,v}<0$作为$v$，然后`pivot(u,v)`后$a_{i,0}$就变正了... + +但这可能出现循环，然后不断损精度，然后WA/TLE，总之很玄学就对了 + +但是有『优先前面/后面』的选择倾向在加以小随机效果似乎不错 + +```cpp +int init() +{ + while(1) + { + int u=0,v=0; + for(int i=1;i<=m;i++) + if(a[i][0]<-eps&&(!u||rand()%3==0)) + u=i; + if(!u)//说明这个解是可行解 + return 1; + for(int i=1;i<=n;i++) + if(a[u][i]<-eps&&(!v||rand()%3==0)) + v=i; + if(!v)//说明无解 + { + printf("Infeasible\n"); + return 0; + } + pivot(u,v); + } +} +``` + +#### 求解 + +每次找一个$v$使得$A_{0,v}>0$再找一个$u$使得$A_{u,v}>0$，然后`pivot(u,v)` + +##### 退化 + +在旋转的过程中，可能会存在保持答案不变的情况，这种现象称为退化。 + +##### Bland规则 + +在选择替入变量和替出变量的时候，我们总是选择满足条件的下标最小值。 + +- 替入变量$v$：目标条件中，系数为正数（如果求最小值就是负数）的第一个作为替入变量 +- 替出变量$u$：对所有的约束条件中，选择对$u$约束最紧的第一个（也就是$\dfrac{A_{u,0}}{A_{u,v}}$最小的那个） + +```cpp +int simplex() +{ + while(1) + { + int u=0,v=0; + double mi=inf; + for(int i=1;i<=n;i++) + if(a[0][i]>eps) + { + v=i; + return; + } + if(!v) + return 1; + for(int i=1;i<=m;i++) + if(a[i][v]>eps&&a[i][0]/a[i][v]