ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

Note -「线性规划」学习笔记

2021-05-21 12:34:43  阅读:208  来源: 互联网

标签:frac 线性规划 rep boldsymbol 笔记 Note int cdots operatorname


\(\mathcal{Definition}\)

  线性规划(Linear Programming, LP)形式上是对如下问题的描述:

\[\operatorname{maximize}~~~~z=\sum_{i=1}^nc_ix_i\\\operatorname{s.t.}\begin{cases} \sum_{j=1}^na_{ij}x_j\ge b_i&i=1,2,\cdots,m\\ x_i\ge0&i=1,2,\cdots,n\end{cases} \]

其中,\(\operatorname{maximize} z\) 即最大化 \(z\) 的值,亦可作 \(\operatorname{minimum} z\)(最小化),我们称其为目标函数。\(\operatorname{s.t.}\)(subject to,服从)称为约束条件。像这样的问题形式是线性规划的标准型。对于以下问题:

\[\operatorname{maximize}~~~~z=\sum_{i=1}^nc_ix_i\\\operatorname{s.t.}\begin{cases}\sum_{j=1}^na_{ij}x_j=b_i&i=1,2,\dots,m\\x_i\ge0&i=1,2,\dots,n\end{cases} \]

则是所谓松弛型,因为我们可以通过引入松弛变量将标准型转化为松弛型:

\[\operatorname{maximize}~~~~z=\sum_{i=1}^nc_ix_i\\\operatorname{s.t.}\begin{cases} x_{i+n}=b_i-\sum_{j=1}^na_{ij}x_j&i=1,2,\cdots,m\\ x_i\ge0&i=1,2,\cdots,m+n\end{cases} \]

所以数学考试大方程题的第一问也是线性规划。此时,原有的变量 \(x_1,x_2,\cdots,x_n\) 称为非基变量,松弛用的变量 \(x_{n+1},x_{n+2},\cdots,x_{n+m}\)(\(m\) 为约束个数)称为基变量。更简便地,我们用矩阵来描述标准形式:

\[\operatorname{maximize}~~~~z=\boldsymbol{c}^T\boldsymbol x\\\operatorname{s.t.}\begin{cases} \boldsymbol{Ax}=\boldsymbol b\\ x_i\ge0\end{cases}\\ \]

其中,

\[\boldsymbol c=\begin{pmatrix}c_1\\c_2\\\vdots\\c_n\end{pmatrix}\\\boldsymbol A=\begin{pmatrix} a_{11}&a_{12}&\cdots&a_{1n}\\ a_{21}&a_{22}&\cdots&a_{2n}\\ \vdots&\vdots&\ddots&\vdots\\ a_{m1}&a_{m2}&\cdots&a_{mn}\end{pmatrix} \]

并称 \(\boldsymbol A\) 为约束矩阵,\(\boldsymbol b\) 为资源向量,\(\boldsymbol c\) 为价值向量,\(\boldsymbol x\) 为决策价值变量向量

\(\mathcal{Algorithm}\)

  说了那么多定义,接下来我们来设计一个算法求解线性规划问题。算法的基本思路为:

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

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

  以

\[\operatorname{maximize}~~~~z=3x_1+x_2+2x_3+0x_4+0x_5+0x_6\\ \operatorname{s.t.}\begin{cases} 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_1,x_2,\cdots,x_6\ge0 \end{cases} \]

为例,不难看出初始解可为 \(\boldsymbol x=\begin{pmatrix}0,0,0,30,24,36\end{pmatrix}\)。接着变换第三个式子:

\[x_6=36-(4x_1+x_2+2x_3)\\ \Rightarrow x_1=9-\frac{1}4x_2-\frac{1}2x_3-\frac{1}4x_6 \]

代入 \(z\),得到:

\[z=27+\frac{3}4x_2+\frac{3}2x_3-\frac{3}4x_6 \]

  如此,交换 \(z\) 表达式中的一个非基变量和一个基变量的操作称为转轴(pivot)。形式地,设交换非基变量 \(x_u\) 与基变量 \(x_{v+n}\):

\[x_{v+n}=b_v-\sum_{i=1}^na_{vi}x_i\\\Rightarrow~~~~x_u=\frac{b_v-x_{v+n}-\sum_{i=1}^n[i\not=u]a_{vi}x_i}{a_{vu}} \]

上例中,转轴 \(x_1-x_6\),再将新的等式依次代入其他等式和目标函数,得到了 \(z=27+\frac{3}4x_2+\frac{3}2x_3-\frac{3}4x_6\)。如果我们通过若干次转轴使得 \(z\) 表达式中所有变量系数非正且 \(b_i\) 非负,就能得到此时的常数项为 \(z\) 的最值。总结一下,我们的算法分为两步:

  1. 找初始解。一种比较好写的随机算法:随机取 \(b_i<0\),再随机取 \(a_{ij}<0\),转轴 \(x_j-x_{i+n}\)。不存在 \(b_i\) 时结束。虽然这样仅保证 \(b_i\ge 0\)。更为精准的算法见学长的博客
  2. 转轴使得 \(c_j\le0\)。取 \(j=\arg\max\{c_j|c_j>0\}\),找出 \(i=\arg\min\{\frac{b_i}{a_{ij}}|a_{ij}>0\}\),最后转轴 \(x_j-x_{i+n}\)。如此仍能保持 \(b_i\ge0\)。

  无解当且仅当找不到初始解,这显而易见;无界当且仅当第二步找不到 \(i\),因为此时任意增大 \(c_j>0\) 的非基变量 \(x_j\),只会使 \(z\) 和 \(\boldsymbol b\) 中某些元素增大而不可能减少,一定能够保持合法性。

\(\mathcal{Example}\)

  UOJ #179 \(97\) 分代码:

/* Clearink */

#include <ctime>
#include <cstdio>
#include <random>

#define rep( i, l, r ) for ( int i = l, rpbound##i = r; i <= rpbound##i; ++i )
#define per( i, r, l ) for ( int i = r, rpbound##i = l; i >= rpbound##i; --i )

const int MAXN = 50;
const double EPS = 1e-9;
int n, m, type, id[MAXN + 5], ref[MAXN + 5];
double a[MAXN + 5][MAXN + 5];

inline void iswp ( int& a, int& b ) { a ^= b ^= a ^= b; }
inline double dabs ( const double a ) { return a < 0 ? -a : a; }
inline int dcmp ( const double a, const double b ) {
	return dabs ( a - b ) < EPS ? 0 : ( a < b ? -1 : 1 );
}

inline bool halfp () {
	static std::mt19937 rnd ( time ( 0 ) ^ 20120712 );
	return rnd () & 1;
}

inline void pivot ( const int r, const int c ) {
	iswp ( id[r + n], id[c] );
	double tmp = -a[r][c]; a[r][c] = -1;
	rep ( i, 0, n ) a[r][i] /= tmp;
	rep ( i, 0, m ) if ( i != r && dcmp ( a[i][c], 0 ) ) {
		tmp = a[i][c], a[i][c] = 0;
		rep ( j, 0, n ) a[i][j] += tmp * a[r][j];
	}
}

inline int simplex () {
	rep ( i, 1, n ) id[i] = i;
	while ( true ) {
		int i = 0, j = 0;
		rep ( k, 1, m ) {
			if ( !~dcmp ( a[k][0], 0 ) && ( !i || halfp () ) ) {
				i = k;
			}
		}
		if ( !i ) break;
		rep ( k, 1, n ) {
			if ( dcmp ( a[i][k], 0 ) == 1 && ( !j || halfp () ) ) {
				j = k;
			}
		}
		if ( !j ) return 2;
		pivot ( i, j );
	}
	while ( true ) {
		int i = 0, j = 0;
		double mx = 0, mn = 1e9;
		rep ( k, 1, n ) if ( dcmp ( a[0][k], mx ) == 1 ) mx = a[0][j = k];
		if ( !j ) break;
		rep ( k, 1, m ) {
			if ( !~dcmp ( a[k][j], 0 ) && dcmp ( -a[k][0] / a[k][j], mn ) == -1 ) {
				mn = -a[i = k][0] / a[k][j];
			}
		}
		if ( !i ) return 1;
		pivot ( i, j );
	}
	return 0;
}

int main () {
	scanf ( "%d %d %d", &n, &m, &type );
	rep ( i, 1, n ) scanf ( "%lf", &a[0][i] );
	rep ( i, 1, m ) {
		rep ( j, 1, n ) scanf ( "%lf", &a[i][j] ), a[i][j] *= -1;
		scanf ( "%lf", &a[i][0] );
	}
	int res = simplex ();
	if ( res == 2 ) puts ( "Infeasible" );
	else if ( res ) puts ( "Unbounded" );
	else {
		printf ( "%.12f\n", a[0][0] );
		if ( type ) {
			rep ( i, 1, m ) ref[id[i + n]] = i;
			rep ( i, 1, n ) {
				printf ( "%.12f%c", ref[i] ? a[ref[i]][0] : 0, i ^ n ? ' ' : '\n' );
			}
		}
	}
	return 0;
}

标签:frac,线性规划,rep,boldsymbol,笔记,Note,int,cdots,operatorname
来源: https://www.cnblogs.com/rainybunny/p/14230790.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有