前言
本文章是根据法国国立高等电力技术、电子学、计算机、水力学与电信学校 (E.N.S.E.E.I.H.T.) 第七学期课程*“Recherche Operation”* 总结而来的【部分课程笔记】。碍于本人学识有限,部分叙述难免存在纰漏,请读者注意甄别。
总而言之,运筹学就是从真实系统中建立模型,用数学的形式表示出来。
附录一:KKT 和 拉格朗日乘子
关系:KKT ,拉格朗日乘子 和 线性规划
本科阶段,我们在高等数学中认识了拉格朗日乘子。其研究范围可以认为是
拉格朗日乘子∼在约束条件为强约束(即约束条件都为等式)的情况下,非线性规划的问题中求解最优值。
而 KKT 研究范围可以认为是
KKT∼在约束条件为弱约束(即约束条件有不等式)的情况下,非线性规划的问题中求解最优值。
线性规划问题的研究范围可以认为是
LP(线性规划)∼在约束条件下,线性规划的问题中求解最优值。
所以我们可以认为 KKT>拉格朗日乘子,KKT>LP>整数线性规划
拉格朗日乘子
引例:我们有如下问题
W=f(x,y,z)=3x2+2y2−4z2s.t.{3x+4y−z=06x2+y−z2=0
由于对于拉格朗日乘子法的证明过程并不感兴趣,所以在此不加证明地直接给出其计算方法:
-
已知我们想要求解目标函数 W=f(x,y,z)=3x2+2y2−4z2 的最大值maxW 或最小值minW,我们有两个约束条件g1,g2: s.t.{g1:3x+4y−z=0g2:6x2+y−z2=0
-
我们引入拉格朗日乘子 λi 辅助求解该问题。拉格朗日乘子 λi 引入的数量由约束条件的数量决定。所以本引例中需要引入两个拉格朗日乘子 λ1,λ2 。
-
构造一个包含拉格朗日乘子 λ1,λ2 的新函数 F(x,y,z,λi)=f(x,y,z)+λ1g1+λ2g2+...+λngn 。所以原函数可转化为以下的拉格朗日函数:
F(x,y,z,λ1,λ2)=f(x,y,z)+λ1g1+λ2g2=3x2+2y2−4z2+λ1(3x+4y−z)+λ2(6x2+y−z2)
-
则我们就可以通过导数的方法求解此问题 F(x,y,z,λ1,λ2) 。对 F 求偏导,有几个自变量,就求几次偏导,令其都等于 0。方程 ①②③④⑤ 一定能得到不止一组的解 (xi,yi,zi,λ1i,λ2i) ,将解得的解 (xi,yi,zi,λ1i,λ2i) 代入原来的目标函数 W=f(x,y,z) 中,算出的 Wmax 和 Wmin。
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧∂x∂F=0①∂y∂F=0②∂z∂F=0③∂λ1∂F=0④∂λ2∂F=0⑤⟹解⎝⎜⎜⎜⎜⎛xiyiziλ1iλ2i⎠⎟⎟⎟⎟⎞⟹代入原来的目标函数f(xi,yi,zi)⟹找Min和Max(Wmax,Wmin)
KKT
KKT 本质上是拉格朗日乘数法的推广,即推广到弱约束(即约束条件有不等式)条件。
【拉格朗日乘数法】:
Zmax或Zmin=f(x,y,z,...)g1=0,g2=0,g3=0①观察约束条件,有几个约束,方程就引入几个拉格朗日乘子λ,对于λ无限制②构建F(x,y,z,λ1,λ2,λ3)=f(x,y,z)+λ1g1+λ2g2+λ3g3③⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧∂x∂F=0⟹本质是F对目标函数的自变量求偏导∂y∂F=0⟹本质是F对目标函数的自变量求偏导∂z∂F=0⟹本质是F对目标函数的自变量求偏导∂λ1∂F=0⟹g1=0∂λ2∂F=0⟹g2=0∂λ3∂F=0⟹g3=0⟹约束gi=0④解方程组⑤将所有解代入目标函数f(x,y,z),判断最大值和最小值⑥结束
【KKT】:
Zmax或Zmin=f(x,y,z,...)g1≤0,g2≥0,gn=0⓪化为标准型∗①观察约束条件,有几个约束,方程就引入几个拉格朗日乘子λ∗,λ∗≥0②构建L(x,y,z,λ1∗,λ2∗,λ3∗,λ4∗)=f(x,y,z)+λ1∗g1′+λ2∗g2′+λ3∗g3′+λ4∗g4′=f(x,y,z)+λ1∗g1−λ2∗g2+λ3∗g3−λ4∗g3③⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧∂x∂L=0⟹本质是F对目标函数的自变量求偏导∂y∂L=0⟹本质是F对目标函数的自变量求偏导∂z∂L=0⟹本质是F对目标函数的自变量求偏导λ1∗g1=0λ2∗(−g2)=0λ3∗(g3)=0λ4∗(−g3)=0⟹广义乘子λi∗×约束gi=0④解方程组⑤将所有解代入目标函数f(x,y,z),判断最大值和最小值⑥结束
-
化为【标准型】:
对于上述例子,我们可以看到【约束条件】很“混乱”,即符号不同,有≤,≥,= ,我们需要将它们化为统一标准形式(“代数式≥0”)。
- 对于 = :g3=0⇔g3≤0且g3≥0
- 对于 ≤ :g1≤0⇔−g1≥0
所以我们可以按上述规则将约束条件改写为
s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧−g1≥0g2≥0g3≥0−g3≥0⟹看作整体s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧g1′=−g1≥0g2′=g2≥0g3′=g3≥0g4′=−g3≥0
不过更建议化成 “代数式≤0” 的形式,因为可以更好的和拉格朗日乘数法结合观察:
s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧−g1≥0g2≥0g3≥0−g3≥0⟹看作整体s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧g1′=g1≤0g2′=−g2≤0g3′=g3≤0g4′=−g3≤0
-
观察约束条件,有几个约束,方程就引入几个广义拉格朗日乘子 λ∗:λ1∗,λ2∗,λ3∗,λ4∗ 。注意,在拉格朗日乘数法中,引入拉格朗日乘子 λ 后并无任何限制,但是在 KKT 中,我们要求所有的广义拉格朗日乘子
λ∗≥0
-
构造一个包含广义拉格朗日乘子 λ1∗,λ2∗,λ3∗,λ4∗ 的新函数 L(x,y,z,λi∗)=f(x,y,z)+λ1∗g1+λ2∗g2+...+λn∗gn 。所以原函数可转化为以下的广义拉格朗日函数:
L(x,y,z,λ1∗,λ2∗,λ3∗,λ4∗)=f(x,y,z)+λ1∗g1′+λ2∗g2′+λ3∗g3′+λ4∗g4′=f(x,y,z)+λ1∗g1−λ2∗g2+λ3∗g3−λ4∗g3
-
则我们就可以通过导数的方法求解此问题 L(x,y,z,λ1∗,λ2∗,λ3∗,λ4∗) 。写出如下形式的方程。方程定能得到不止一组的解 (xi,yi,zi,λ1∗i,λ2∗i,λ3∗i,λ4∗i ,将解得的解 (xi,yi,zi,λ1∗i,λ2∗i,λ3∗i,λ4∗i 代入原来的目标函数 Z=f(x,y,z) 中,算出的 Zmax 和 Zmin。
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧∂x∂L=0∂y∂L=0∂z∂L=0λ1∗g1=0λ2∗(−g2)=0λ3∗(g3)=0λ4∗(−g3)=0⟹解⎝⎜⎜⎜⎜⎜⎜⎜⎜⎛xiyiziλ1∗iλ2∗iλ3∗iλ4∗i⎠⎟⎟⎟⎟⎟⎟⎟⎟⎞⟹代入原来的目标函数f(xi,yi,zi)⟹找Min和Max(Zmax,Zmin)
例题
f(x)=(x−3)2,0≤x≤5,求minf(x)
-
规范化(化标准型“代数式≤0”):
0≤x≤5⇒x≥0x≤5⇒−x≤0x−5≤0⇒g1=−x≤0g2=x−5≤0
-
观察约束,有 2 个约束,所以引入 2 个广义拉格朗日乘子 λ1∗,λ2∗,
λ1∗≥0λ2∗≥0
-
构建广义拉格朗日函数 L(x,λ1∗,λ2∗) :
L(x,λ1∗,λ2∗)=f(x)+λ1∗g1+λ2∗g2=(x−3)2+λ1∗(−x)+λ2∗(x−5)
-
联立方程组:
⎩⎪⎪⎪⎨⎪⎪⎪⎧∂x∂F=2(x−3)−λ1∗+λ2∗=0①λ1∗g1=λ1∗(−x)=0②λ2∗g2=λ2∗(x−5)=0③
-
解方程组。(【注意⚠️】:一般需要详细讨论解的情况。)
【1】我们先从 ② 式开始 :②⇒λ1∗x=0⟹两种情况λ1∗=0或x=0
-
当 λ1∗=0 时,我们可以得到 x=0 。我们再结合 ③λ2∗(x−5)=0,我们得到 λ2∗=0 。我们再结合 $ ① ; 2(x-3) - \lambda_1^* + \lambda_2^* = 0$ ,我们可以解得 λ1∗=−6<0 。因为我们有约束条件 λ1∗≥0,所以我们要舍弃这组解。
-
当 λ1∗=0 时,在 ② 式中 x 可以是任意值,所以我们要结合 ① ③ 式:
{2(x−3)+λ2∗=0①λ2∗(x−5)=0③⟹解得{x=3λ2∗=0或{x=5λ2∗=−4<0
我们可以看出,其中存在一组解
⎩⎪⎨⎪⎧x=3λ1∗=0λ2∗=0
【2】我们再回到 ③ 式重新讨论:$ ③ ; \lambda_2^* (x-5) = 0 \quad \stackrel{两种情况}{\Longrightarrow} \quad \lambda_2^*=0 ; 或; x=5$
-
x=5,λ2∗=0 的情况已经在 ② 式的第 2 种情况中讨论过了,所以跳过。
-
当 λ2∗=0 时,在 ③ 式中 x 可以是任意值,所以我们要结合 ① ② 式:
{2(x−3)−λ1∗=0①λ1∗(−x)=0②⟹解得{x=0λ1∗=−6<0或{x=3λ1∗=0
我们可以看出,其中存在一组解
⎩⎪⎨⎪⎧x=3λ1∗=0λ2∗=0
-
所以 minf(x)=f(3)=(3−3)2=0 。
【一点想法】
当我们使用KKT得到一个解时,若该解中某个广义拉格朗日乘子 λi∗>0,则说明 λi 所对应的约束 gi=0,即这个解所对应的点恰好落在了曲线 gi=0 上。
附录二:线性代数回顾
简言之,线性代数的目的就是求解线性方程组。经典的理论是与几何联系起来。
方法论
graph LR
提取系统信息 --经过运算--> 判断解的情况
提取系统信息 --经过运算--> 解方程
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧a1,1x1+a1,2x2+...a1,nxn=b1a2,1x1+a2,2x2+...a2,nxn=b2............ai,1x1+ai,2x2+...ai,nxn=bi............an,1x1+an,2x2+...an,nxn=bn⟹提取系统信息⎝⎜⎜⎜⎛a1,1a2,1⋮an,1a1,2a2,2⋮an,2⋯⋯⋮⋯a1,na2,n⋮an,n⎠⎟⎟⎟⎞⎝⎜⎜⎜⎛x1x2⋮xn⎠⎟⎟⎟⎞=⎝⎜⎜⎜⎛b1b2⋮bn⎠⎟⎟⎟⎞
【秩】Rank:矩阵的本质属性。如果把矩阵看成一个个行向量或者列向量,如果其中有一个向量能被另一个向量线性表示的话,那么该矩阵表示的方程组有无穷多解;如果没有任何一个向量能被另一个向量线性表示的话,那么该矩阵表示的方程组有唯一解。【秩】的本质就是矩阵中独立向量的个数。
【行列式】:行列式的本质是代表一种算法的运算符,和 +,−,×,÷ 类似。用 A2×2 举例:
A2×2=∣∣∣∣xzyv∣∣∣∣=xv−yz
【行列式的几何意义】
- 二阶:是这 2 个向量围成的平行四边形的面积
- 三阶:是这 3 个向量围成的平行六面体的体积
- n阶:n维超立方体的体积
所以,行列式一旦 =0 ,说明行列式对应的矩阵里有线性相关的向量。当行列式 =0 时,说明行列式对应的矩阵是满秩的。
第一部分:线性规划
1.1 线性规划问题及数学模型
1.1.1 问题的提出
我们可以考虑以下一个简单的例子:
例1,某工厂在计划期内要安排生产甲、乙两种产品,已知生产单位产品所需的设备合时及 A、B 两种原材料的消耗如表所示。该厂每生产一件产品甲可获利2元,一件产品乙可获利3元,问该工厂如何安排生产计划使利润最大?
|
产品甲 |
产品乙 |
资源限量 |
设备 |
1 台时/件 |
2 台时/件 |
8 台时 |
原材料 A |
4 kg/件 |
0 |
16 kg |
原材料 B |
0 |
4 kg/间 |
12 kg |
利润 |
2 元 |
3 元 |
|
我们可以根据以上描述建立一个数学模型:
Profitmax=2x+3ys.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧x1+2x2≤84x1≤164x2≤12x1,x2≥0
其中,
- x1 和 x2 为【决策变量】;
- Profitmax=2x+3y 是【决策函数】,是决策者希望达到的计划目标;
- s.t. 为“subject to”,为【约束条件】,为决策变量取值时受到的限制。
这类优化问题的共同特征:
- 每一个问题都用一组决策变量表示某一方案,这组决策变量的值就代表一个具体方案。一般取值非负且连续。
- 要有建模的相关数据,如资源拥有量、消耗资源定额、创造的新价值量等并构造不矛盾的约束条件,由线性等式或不等式来表示。
- 都有一个要达到的目标,可用决策变量及其有关的价值系数构成的线性函数来表示,一般要求最大或最小
线性规划问题模型的【一般表达式形式】为:
max(min)Z=c1x1+c2x2+c3x3+...+cnxns.t.⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧a1,1x1+a1,2x2+...a1,nxn≤(≥,=)b1a2,1x1+a2,2x2+...a2,nxn≤(≥,=)b2............ai,1x1+ai,2x2+...ai,nxn≤(≥,=)bi............an,1x1+an,2x2+...an,nxn≤(≥,=)bnx1,x2,...xn≥0
其中,
- cn 为价值系数
- aij 为技术系数
- bi 为资源限制
- x1,x2,...xn≥0 为决策变量约束
1.1.2 图解法
图解法求解步骤(只适用于两个决策变量的问题 - 二维):
- 由全部约束条件作图求出【可行域】;
- 作目标函数【等值线】,确定使目标函数最优的移动方向;
- 平移目标函数的等值线,找出【最优点】,算出最优值。
1.1.3 化标准型
参考以上线性规划问题模型的一般表达式形式,我们发现线性规划的共同特证:
-
决策变量:每个问题都用一组决策变量表示某个方案
决策变量:决策变量的取值一般都是非负且连续的
-
约束条件:与决策变量不矛盾的条件,用线性等式或不等式表示
-
目标函数:决策变量与价值系数组成,一般要求实现最大或最小
所以我们建模的一般思路为:确定决策变量,写出目标函数,找出约束条件。
如果我们想用单纯形法求解线性规划问题,第一步是将一般表达式为以下的标准型:
maxZ=c1x1+c2x2+c3x3+...+cnxns.t.⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧a1,1x1+a1,2x2+...a1,nxn=b1a2,1x1+a2,2x2+...a2,nxn=b2............ai,1x1+ai,2x2+...ai,nxn=bi............am,1x1+am,2x2+...am,nxn=bmx1,x2,...xn≥0b1,b2,...bn≥0
简言之:
- 目标函数 Z 最大;
- 约束条件 ai,1x1+ai,2x2+...ai,nxn=bi 为等式;
- 决策条件 x1,x2,...xn 非负
- 资源限量 b1,b2,...bn 非负
上方的单纯形法的标准型可以简写成以下形式
maxZ=j=1∑ncjxjs.t.⎩⎪⎪⎨⎪⎪⎧j=1∑nai,jxj=bji=1,2,3,...,mxj≥0j=1,2,3,...,n
或写成向量的形式
maxZ=CXs.t.⎩⎪⎪⎨⎪⎪⎧j=1∑nPjxj=bxj≥0j=1,2,3,...,n价值向量C=(c1,c2,...,cn);决策变量的向量X=(x1,x2,...,xn)T资源向量b=(b1,b2,...,bm)T;约束条件的系数列向量Pj=(a1,j,a2,j,...,am,j)
或写成矩阵的形式
maxZ=CXs.t.{AX=bX≥O其中,A=⎝⎜⎜⎜⎛a1,1a2,1⋮am,1a1,2a2,2⋮am,2⋯⋯⋮⋯a1,na2,n⋮am,n⎠⎟⎟⎟⎞,O=⎝⎜⎜⎜⎛00⋮0⎠⎟⎟⎟⎞,X=⎝⎜⎜⎜⎛x1x2⋮xn⎠⎟⎟⎟⎞
一般形式化成标准型
例2:我们需要将以下一般形式化成标准型
minZ=x1+2x2+3x3①s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧−2x1+x2+x3≤9②−3x1+x2+2x3≥4③4x1−2x2−3x3=−6④x1≤0,x2≥0,x3无约束
-
【目标函数最大】
如果目标函数的求解是最小值问题,我们需要将其转换成最大值问题,方法为取负数:
minZ=CX⟹Z′=−ZmaxZ′=−CX
所以在例题中:
minZ=x1+2x2+3x3⟹Z′=−ZmaxZ′=−x1−2x2−3x3
-
【资源限量非负】
bi<0⟹两端同乘(−1)右端项非负bi=0的情况日后再议
我们可以发现,例中的 ④ 式的资源限量为 -6,所以
④4x1−2x2−3x3=−6⟹两端同乘(−1)−4x1+2x2+3x3=6
-
【约束条件化为等式】
将约束条件中的≤⟹【加上】松弛变量变成=将约束条件中的≥⟹【减掉】剩余变量变成=
例如例中的 ② 式,我们需要将其化为等式,则
②−2x1+x2+x3≤9⟹【加上】松弛变量−2x1+x2+x3+x4=9
注意⚠️:此时的 x4 实际上是无意义的,只是我们加入的使等式成立的松弛变量。
同理,对于 ③ 式,我们有
③−3x1+x2+2x3≥4⟹【减掉】剩余变量−3x1+x2+2x3−x5=4
松弛变量与剩余变量在实际问题中分别表示未被充分利用的资源和超出的资源,均未转化为价值和利润,所以引进模型后它们**在目标函数中的系数均为 0 **。
-
【决策变量非负】
x≤0⟹x′=−xx′≥0
例如例题中的 x1≤0,我们需要将其化为大于 0 的形式,则
x1≤0⟹x1′=−x1x1′≥0
-
【处理无约束变量】
x无约束⟹x=x′−x′′x′≥0,x′′≥0
意为,一个无约束的数 x 可以等于两个非负数 x′,x′′ 的差,其结果依然无约束。但是 x′,x′′ 是存在约束的。
例如例题中的 x3无约束,我们将其变化为
x3无约束⟹x3=x3′−x3′′x3′≥0,x3′′≥0
-
最后请注意替换后的变量代入原式时的替换问题
所以,通过以上步骤,例题中的原式可化为如下形式:
maxZ=x1′−2x2−3x3′+3x3′′+0x4+0x5①s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧2x1′+x2+x3′−x3′′+x4=9②3x1′+x2+2x3′−2x3′′−x5=4③−4x1′+2x2+3x3′−3x3′′=6④x1′,x2,x3′,x3′′,x4,x5≥0
1.1.4 解的概念
可行解 和 最优解
在以下线性规划问题的一般形式中,
maxZ=j=1∑ncjxj①s.t.⎩⎪⎪⎨⎪⎪⎧j=1∑nai,jxj=bji=1,2,3,...,m②xj≥0j=1,2,3,...,n③
- 【可行解】:满足约束条件 ② 和 ③ 的解。全部可行解的集合称为可行域
- 【最优解】:是目标函数 ① 最大的可行解
基解 和 基可行解
maxZ=CXs.t.{AX=bX≥O其中,C=⎝⎜⎜⎜⎛c1c2⋮cn⎠⎟⎟⎟⎞,X=⎝⎜⎜⎜⎛x1x2⋮xn⎠⎟⎟⎟⎞,Am<n=⎝⎜⎜⎜⎛a1,1a2,1⋮am,1a1,2a2,2⋮am,2⋯⋯⋮⋯a1,na2,n⋮am,n⎠⎟⎟⎟⎞,b=⎝⎜⎜⎜⎛b1b2⋮bn⎠⎟⎟⎟⎞,O=⎝⎜⎜⎜⎛00⋮0⎠⎟⎟⎟⎞
- 【基】:设 A 是约束方程组 ② 的 m×n 阶系数矩阵(设 n>m,变量的个数大于方程的个数,即方程组一定有无穷多个解),其秩(rank) 为 m,B 是 A 中的一个 m×m 阶的满秩子矩阵(行列式 ∣B∣=0 的非奇异子矩阵),称 B 是线性规划问题的一个基。
R(A)=m⇒⎝⎜⎜⎜⎛a1,1a2,1⋮am,1a1,2a2,2⋮am,2⋯⋯⋮⋯a1,ma2,m⋮am,m⎠⎟⎟⎟⎞⇒(P1P2⋯Pm),其中Pi=⎝⎜⎜⎜⎛a1,ia2,i⋮am,i⎠⎟⎟⎟⎞
- 【基向量组 B】(极大线性无关组):向量组 A 中有 m 个线性无关的向量,且其余向量均可由这 m 个向量来线性表示
⎝⎜⎜⎜⎛a1,1a2,1⋮am,1a1,2a2,2⋮am,2⋯⋯⋮⋯a1,na2,n⋮am,n⎠⎟⎟⎟⎞×X⇒⎝⎜⎜⎜⎛a1,1a2,1⋮am,1⎠⎟⎟⎟⎞x1+⎝⎜⎜⎜⎛a1,2a2,2⋮am,2⎠⎟⎟⎟⎞x2+...+⎝⎜⎜⎜⎛a1,na2,n⋮am,n⎠⎟⎟⎟⎞xn=⎝⎜⎜⎜⎛b1b2⋮bm⎠⎟⎟⎟⎞
设方程组前 m 个变量的系数列向量就是它的基向量(极大线性无关组),则
⎝⎜⎜⎜⎛a1,1a2,1⋮am,1⎠⎟⎟⎟⎞x1+⎝⎜⎜⎜⎛a1,2a2,2⋮am,2⎠⎟⎟⎟⎞x2+...+⎝⎜⎜⎜⎛a1,na2,n⋮am,n⎠⎟⎟⎟⎞xn=⎝⎜⎜⎜⎛b1b2⋮bm⎠⎟⎟⎟⎞⇒
⎝⎜⎜⎜⎛a1,1a2,1⋮am,1⎠⎟⎟⎟⎞x1+⎝⎜⎜⎜⎛a1,2a2,2⋮am,2⎠⎟⎟⎟⎞x2+...+⎝⎜⎜⎜⎛a1,ma2,m⋮am,m⎠⎟⎟⎟⎞xm=⎝⎜⎜⎜⎛b1b2⋮bm⎠⎟⎟⎟⎞−⎝⎜⎜⎜⎛a1,m+1a2,m+1⋮am,m+1⎠⎟⎟⎟⎞xm+1−...−⎝⎜⎜⎜⎛a1,na2,n⋮am,n⎠⎟⎟⎟⎞xn
其中,P1,P2,...,Pm 被称为【基向量】,它们对应的变量 x1,x2,...,xm 被称作【基变量】;Pm+1,...,Pn 被称为【非基向量】,它们对应的变量 xm+1,...,xn 被称作【非基变量】。
令所有的非基变量 xm+1=xm+2=...=xn=0,又因为 ∣B∣=0,根据克莱姆法则,可以求出一个唯一解:
⎝⎜⎜⎜⎛a1,1a2,1⋮am,1⎠⎟⎟⎟⎞x1+⎝⎜⎜⎜⎛a1,2a2,2⋮am,2⎠⎟⎟⎟⎞x2+...+⎝⎜⎜⎜⎛a1,ma2,m⋮am,m⎠⎟⎟⎟⎞xm=⎝⎜⎜⎜⎛b1b2⋮bm⎠⎟⎟⎟⎞⟹解得XB⎝⎜⎜⎜⎛x1x2⋮xm⎠⎟⎟⎟⎞
【基解】:根据基 B 求得的解 X 被称作基解。
X=(XB+XN)=(x1,x2,...,xm,0,0,...,0)T
其中,
- 基解中非零分量的数目不大于 m(方程个数);
- 有一个基,就能求得一组基解。
【基可行解】:基解中所有分量都满足非负条件的解,即
X=(XB+XN)=(x1,x2,...,xm,0,0,...,0)T,xi≥0,i∈1,2,...,n
【可行即】:对应于基可行解的基,称为可行基。
四个解的关系
- 当最优解唯一时,最优解也是基可行解;
- 当最优解不唯一时,最优解不一定是基可行解。
1.2 单纯形法
单纯形法的思想
我们继续探究例2,用基的方法求解最优解。在 1.1.3 节的末尾,我们得到了化为标准型的线性规划问题模型:
maxZ=x1′−2x2−3x3′+3x3′′+0x4+0x4①s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧2x1′+x2+x3′−x3′′+x4=9②3x1′+x2+2x3′−2x3′′−x5=4③−4x1′+2x2+3x3′−3x3′′=6④x1′,x2,x3′,x3′′,x4,x5≥0
-
首先,我们可以将约束条件中的系数矩阵 A 写出来:
A=⎝⎛23−4112123−1−2−3100010⎠⎞⇒(P1P2P3P4P5P6)
-
我们可以看出,这是 A 是一个 3×6 ,秩为 3 的矩阵,则我们需要找到一个行列式不为 0 的 3×3 的基矩阵 B
∣B1∣=∣∣P1P2P3∣∣=∣∣∣∣∣∣23−4112123∣∣∣∣∣∣=−9
-
因为 P1,P2,P3 线性无关,令 (P1P2P3) 为基,则
⎝⎛23−4112123−1−2−3100010⎠⎞×⎝⎜⎜⎜⎜⎜⎜⎛x1′x2x3′000⎠⎟⎟⎟⎟⎟⎟⎞=⎝⎛946⎠⎞⇒X(1)=⎝⎜⎜⎜⎜⎜⎜⎛x1′x2x3′000⎠⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎛0.777813.2222−5.77778000⎠⎟⎟⎟⎟⎟⎟⎞
-
求解出的基解有元素小于 0,所以这个基解 X(1)不是可行解。
-
重复第 2 步,找到一个行列式不为 0 的 3×3 的基矩阵 Bi,解得它的基解 X(i) 。若 X(i) 的所有元素都满足非负约束(≥0),则 X(i) 就是可行解。
- 【出基】:将列向量 Pi 从基里移除,即
Bi=(PaPiPb)=⎝⎛a1,aa2,aa3,aa1,ia2,ia3,ia1,ba2,ba3,b⎠⎞⟹Pi出基(PaPb)=⎝⎛a1,aa2,aa3,aa1,ba2,ba3,b⎠⎞
【进基】:将列向量 Pj 移进基中,即Bj=(PaPb)=⎝⎛a1,aa2,aa3,aa1,ba2,ba3,b⎠⎞⟹Pj进基(PaPjPb)=⎝⎛a1,aa2,aa3,aa1,ja2,ja3,ja1,ba2,ba3,b⎠⎞
尽量选择目标函数里对结果影响大的向量【进基】,因为我们不希望它是非基向量而使其基变量为 0;尽量选择“容量”最小的列向量【出基】
-
解出所有的基解 X(n),其中基解的数目 n≤C64=15 个。找出所有基解中的可行解
-
再将所有可行解代入目标函数 ① 中,结果最大的就是【最优解】
单纯形法的计算步骤
- 为了书写规范和便于计算,对单纯形法的计算设计了【单纯形表】
- 每一次迭代对应一张单纯形表
- 含初始基可行解的单纯形表称为【初始单纯形表】
- 含最优解的单纯形表称为【最终单纯形表】
单纯形表
单纯形表解题方法:
- 将题目给出的线性规划问题化为标准型
- 根据标准型填写初始单纯形表
- 列数: 未知数xi的个数+cB+XB+b+θi
- 行数:约束方程的个数+cj+σj
- cj:目标函数中各个变量的系数,对应地填入表内
- 每个变量所在列的数字根据大括号中约束方程组中该变量的系数
- 【基】找出变量下放已填写数字构成的矩阵中的单位矩阵,依次将该单位矩阵对应的变量写在 XB 的下面
- 将上一步骤中找到的变量上方的数字 cj 依次对应地填入 cB 中
- 将大括号中约束方程组等号右边的 bi 依次填到 b 下面
- 找出可行解
- 令 XB 所在列的变量与**b 所在列的数字**对应相等,再令其他变量等于 0
- 求出检验数 σj=cj−(cB1xj1+cB2xj2+...)
- XB 下面的变量对应的 σ 的值为 0
- 【初始单纯形表】中 cj 的数字与 σj 行的数字完全相同
- 观察一下 σj 这一行的数字看一下是否都 ≤0
- 若这些数字都 ≤0,则该可行解就是最优解
- 若这些数字有>0的,则该可行解不是最优解,继续进行步骤 6
- 找到 σj 行最大的数字 max(σj) 那一列对应的变量 xa【进基变量】,求出 θi=bi/xai:
- 若 θ≥0,则把 θ 的值填到表中
- 若 θ<0,则不用把 θ 的值填到表中
- 若 θ无意义,则不用把 θ 的值填到表中
- 找到表中 θi 最小值对应 XB 列的变量 xb(出基变量),找到【 xa 的列】和 【xb 的行】交叉的数字 m
- 用 【xa 上面的数字 ca】 替代 【xb 前面的数字 cb】,用 xa 替代 xb,清空 σj 行与 θi 列
- 对 x1,x2…xn 与 b 列组成的矩阵进行运算,【将 m 变成 1】(矩阵的初等行变换),同列其他元素变成 0,形成一个新的矩阵,将该矩阵中的数字填入表格中对应的位置形成新的单纯形表并进行步骤 3
- 如何检验所求的最优解是不是唯一最优解:
- 找出 XB 列变量以外的变量
- 找出 XB 列变量以外的变量对应的 σi
- 若 XB 列变量以外的变量对应的 σi 都 <0,则 X∗ 为唯一最优解
- 若 XB 列变量以外的变量对应的 σi 中有 =0 的,则 X∗ 是无穷多最优解中的一个
【例】给出下列问题:
maxZ=2x1+3x2+0x3+0x4+0x5①s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧x1+2x2+x3+0x4+0x5=8②4x1+0x2+0x3+x4+0x5=16③0x1+4x2+0x3+0x4+x5=12④x1,x2,x3,x4,x5≥0
-
该线性规划问题已化成标准形式;
-
根据标准型填写如下的初始单纯形表。
\ |
cj: |
2 |
3 |
0 |
0 |
0 |
\ |
\ |
cB |
XB |
x1 |
x2 |
x3 |
x4 |
x5 |
b |
θi |
0 |
x3 |
1 |
2 |
1 |
0 |
0 |
8 |
4 |
0 |
x4 |
4 |
0 |
0 |
1 |
0 |
16 |
❌ |
0 |
x5 |
0 |
4 |
0 |
0 |
1 |
12 |
3 |
\ |
σi: |
2 |
3 |
0 |
0 |
0 |
\ |
\ |
-
找出可行解:令 XB 所在列的变量与**b 所在列的数字**对应相等,再令其他变量等于 0。故
X(0)=(x3=8,x4=16,x5=12,x1=0,x2=0)⟹整理X(0)=(0,0,8,16,12)T
-
求出检验数 σj=cj−(cB1xj1+cB2xj2+...) :
- σ1=c1−(cB1x11+cB2x12+cB3x13)=2−(0×1+0×4+0×0)=2
- σ2=c2−(cB1x21+cB2x22+cB3x23)=3−(0×2+0×0+0×4)=3
- σ3=c3−(cB1x31+cB2x32+cB3x33)=0−(0×1+0×0+0×0)=0
- σ4=c4−(cB1x41+cB2x42+cB3x43)=0−(0×0+0×1+0×0)=0
- σ5=c5−(cB1x51+cB2x52+cB3x53)=0−(0×0+0×0+0×1)=0
-
观察一下 σj 这一行的数字看一下是否都 ≤0 ?否,继续进行下一步。
-
找到 σj 行最大的数字 max(σj) 那一列对应的变量 xa【进基变量】,求出 θi=bi/xai:
max(σj)=3=θ2⟹对应的进基变量xa=x2
- θ1=b1/x21=8/2=4>0,填入表中
- θ2=b2/x22=16/0:无意义
- θ3=b3/x31=12/4=3>0,填入表中
-
找到表中 θi 最小值对应 XB 列的变量 xb【出基变量】,找到 xa 的列和 xb 的行交叉的数字 m:
min(θi)=3=θ3⟹对应的出基变量xb=x5。xa=x2的列 与 xb=x5 的行交叉的数字 m=4
-
用 xa 上面的数字替代 xb 前面的数字,用 xa 替代 xb,清空 σj 行与 θi 列:
\ |
cj: |
2 |
3 |
0 |
0 |
0 |
\ |
\ |
cB |
XB |
x1 |
x2 |
x3 |
x4 |
x5 |
b |
θi |
0 |
x3 |
1 |
2 |
1 |
0 |
0 |
8 |
|
0 |
x4 |
4 |
0 |
0 |
1 |
0 |
16 |
|
3 |
x2 |
0 |
4 |
0 |
0 |
1 |
12 |
|
\ |
σi: |
|
|
|
|
|
\ |
\ |
-
对 x1,x2…xn 与 b 列组成的矩阵进行运算,将 m 变成 1,同列其他元素变成 0,形成一个新的矩阵,将该矩阵中的数字填入表格中对应的位置形成新的单纯形表并进行步骤 3:
M=(x1,x2…xn∣b)=(A∣b)=⎝⎛14020410001000181612⎠⎞m=1的初等行变换1/4×第3行⎝⎛140201100010001/48163⎠⎞m所在列其他元素为0第1行+(−2)×第3行⎝⎛140001100010−1/201/42163⎠⎞
\ |
cj: |
2 |
3 |
0 |
0 |
0 |
\ |
\ |
cB |
XB |
x1 |
x2 |
x3 |
x4 |
x5 |
b |
θi |
0 |
x3 |
1 |
0 |
1 |
0 |
-1/2 |
2 |
|
0 |
x4 |
4 |
0 |
0 |
1 |
0 |
16 |
|
3 |
x2 |
0 |
1 |
0 |
0 |
1/4 |
3 |
|
\ |
σi: |
|
|
|
|
|
\ |
\ |
重复步骤 3 …
-
经过多次迭代,列出每次迭代过后的单纯形表:
-
第一次迭代:
\ |
cj: |
2 |
3 |
0 |
0 |
0 |
\ |
\ |
cB |
XB |
x1 |
x2 |
x3 |
x4 |
x5 |
b |
θi |
0 |
x3 |
1 |
2 |
1 |
0 |
0 |
8 |
4 |
0 |
x4 |
4 |
0 |
0 |
1 |
0 |
16 |
❌ |
0 |
x5 |
0 |
4 |
0 |
0 |
1 |
12 |
3 |
\ |
σi: |
2 |
3 |
0 |
0 |
0 |
\ |
\ |
X(0)=(0,0,8,16,12)T
-
第二次迭代:
\ |
cj: |
2 |
3 |
0 |
0 |
0 |
\ |
\ |
cB |
XB |
x1 |
x2 |
x3 |
x4 |
x5 |
b |
θi |
0 |
x3 |
1 |
0 |
1 |
0 |
-1/2 |
2 |
2 |
0 |
x4 |
4 |
0 |
0 |
1 |
0 |
16 |
4 |
3 |
x2 |
0 |
1 |
0 |
0 |
1/4 |
3 |
❌ |
\ |
σi: |
2 |
0 |
0 |
0 |
-3/4 |
\ |
\ |
X(1)=(0,3,2,16,0)T
-
第三次迭代:
\ |
cj: |
2 |
3 |
0 |
0 |
0 |
\ |
\ |
cB |
XB |
x1 |
x2 |
x3 |
x4 |
x5 |
b |
θi |
2 |
x1 |
1 |
0 |
1 |
0 |
-1/2 |
2 |
❌ |
0 |
x4 |
0 |
0 |
-4 |
1 |
2 |
8 |
4 |
3 |
x2 |
0 |
1 |
0 |
0 |
1/4 |
3 |
12 |
\ |
σi: |
0 |
0 |
-2 |
0 |
1/4 |
\ |
\ |
X(2)=(2,3,0,8,0)T
-
第四次迭代:
\ |
cj: |
2 |
3 |
0 |
0 |
0 |
\ |
\ |
cB |
XB |
x1 |
x2 |
x3 |
x4 |
x5 |
b |
θi |
2 |
x1 |
1 |
0 |
0 |
0 |
1/4 |
4 |
|
0 |
x5 |
0 |
0 |
-2 |
1/2 |
1 |
4 |
|
3 |
x2 |
0 |
1 |
1/2 |
-1/8 |
0 |
2 |
|
\ |
σi: |
0 |
0 |
-3/2 |
-1/8 |
0 |
\ |
\ |
X(3)=(4,2,0,0,4)T
我们观察观察一下 σj 这一行的数字看一下是否都 ≤0?是,则该可行解就是最优解,即X∗=X(3)=(4,2,0,0,4)T
所以,maxZ=2x1+3x2+0x3+0x4+0x5=14
第二部分:整数线性规划
整数线性规划模型是在一般线性规划模型的基础上加上决策变量取整数值的约束,由此,我们可以得到整数线性规划问题的一般形式如下:
maxZ=j=1∑ncjxj①s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧j=1∑nai,jxj=bji=1,2,3,...,m②xj≥0j=1,2,3,...,n③x1,x2,x3,...,xn∈Z④
- 先不考虑 x∈Z,求解问题;
- 再不考虑 x∈Z,结合 1 中的解,找到最优解。
实际上,整数线性规划根据决策变量取整数的情况可以进行如下区分:
-
混合整数规划 & 全整数规划
-
0-1型整数规划 & 非0-1型整数规
含 0-1 变量的建模
对不含0-1变量的整数规划问题,我们只需要按照一般线性规划问题的形式进行建模,然后综合问题的性质给相应变量增加取整的约束即可;但对于包含0-1变量的问题,由于0-1变量作为二进制变量(或逻辑量),其应用性很强,在实际问题的建模中常常需要灵活运用。
-
表示系统是否处于某个特定状态,或决策时是否取某个特定方案,如
x={1决策取方案P时0决策不取方案P时(Pˉ)
-
当问题含多项要素,每项要素面临两种选择时,可用0-1变量来描述。
设问题有有限项要素 E1,E2,...,En,其中每项E,有两种选择 Aj 和 Aˉj(j=1,2,...,n),可令:
xj={1若Ej选择Aj0若Ej选择Aˉj
含有相互排斥的约束条件的问题
我们可以设两个 0-1 变量 y1,y2:
y1={1若Ej选择Aj0若Ej选择Aˉj
分支定界法
可以理解为“二分比较”,对可行域进行二分,如果在当前可行域内有最优解
maxZ=40x1+90x2①s.t.⎩⎪⎨⎪⎧9x1+7x2≤56②7x1+20x2≤70③x1,x2≥0x1,x2∈Z④
-
我们先设不考虑 x1,x2∈Z 的情况为B问题:
maxZ=40x1+90x2①s.t.⎩⎪⎨⎪⎧9x1+7x2≤56②7x1+20x2≤70③x1,x2≥0④
-
运用单纯形法,我们可以得出B问题的最优解 XB∗=(4.81,1.82)T,maxZB=356
-
我们可以根据整数对可行域进行二分,然后再讨论:
-
不考虑 x1,x2∈Z 的情况下 x1∗=4.81,所以我们可以在左右两边对区域进行划分,即 x1≤4 和 x1≥5 ,将新得到的约束分别加入到 B 问题中
(1) 问题 B1
maxZ=40x1+90x2①s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧9x1+7x2≤56②7x1+20x2≤70③x1,x2≥0④x1≤4⑤
运用单纯形法,我们可以得出 B1 问题的最优解 XB1∗=(4,2.1)T,maxZB1=349
我们可以在 x2 左右两边对区域进行划分,即 x2≤2 和 x2≥3 ,将新得到的约束分别加入到 B1 问题中:
(i) 问题 B11
maxZ=40x1+90x2①s.t.⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧9x1+7x2≤56②7x1+20x2≤70③x1,x2≥0④x1≤4⑤x2≤2⑥
运用单纯形法,我们可以得出 B1 问题的最优解 XB11∗=(4,2)T,maxZB11=340
(ii) 问题 B12
maxZ=40x1+90x2①s.t.⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧9x1+7x2≤56②7x1+20x2≤70③x1,x2≥0④x1≤4⑤x2≥3⑥
运用单纯形法,我们可以得出 B1 问题的最优解 XB12∗=(1.42,3)T,maxZB12=327
(2) 问题 B2
maxZ=40x1+90x2①s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧9x1+7x2≤56②7x1+20x2≤70③x1,x2≥0④x1≥5⑤
运用单纯形法,我们可以得出 B2 问题的最优解 XB2∗=(5,1.57)T,maxZB2=341
我们可以在 x2 左右两边对区域进行划分,即 x2≤1 和 x2≥2 ,将新得到的约束分别加入到 B2 问题中:
(i) 问题 B21
maxZ=40x1+90x2①s.t.⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧9x1+7x2≤56②7x1+20x2≤70③x1,x2≥0④x1≥5⑤x2≤1⑥
运用单纯形法,我们可以得出 B1 问题的最优解 XB21∗=(5.44,1)T,maxZB21=308
(ii) 问题 B22
maxZ=40x1+90x2①s.t.⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧9x1+7x2≤56②7x1+20x2≤70③x1,x2≥0④x1≥5⑤x2≥2⑥
运用单纯形法,我们可以得出 B1 问题的最优解 XB22∗无可行解
第三部分:动态规划
我们有如下问题,给定一个线路网络,两点之间连线上的数字表示两点间的距离(或费用),试求一段由 A 到 E 的管线铺路,使总距离为最短(或总费用最小)。
其实从A到E并不是一次性的决策,A→B,B→C,C→D,D→E。
由此我们可见,将问题分成若干相互联系的阶段,在每个阶段都需要做出决策,使全过程达到整体最优,这类问题我们称为多阶段决策问题(资源分配、生产存储)。最大特点就是可以把问题拆分成很多阶段。
基本概念
-
阶段:将问题划分成若干个相互联系的阶段。k 表示为阶段变量。
上图有四个阶段,分别是:A→B,B→C,C→D,D→E
-
状态:每个阶段开始时所处的自然状态。sk 表示为状态变量。
上图有10个状态,分别是:A,B1,B2,B3C1,C2D1,D2,D3E
-
决策:当过程处于某一阶段的某一状态时作出的决定,从而确定下一阶段的状态。xk 表示为决策变量。
- 允许决策集合:决策变量的取值范围。Dk(sk) 表示第k阶段从状态sk 出发的允许决策集合。xk∈Dk(sk)
-
策略:按顺序排列的决策集合。
-
状态转移方程:确定过程由一个状态到另一个状态的演变过程。Sk+1=Tk(Sk,xk)。
-
指标函数:衡量所实现过程优劣的数值指标
- 阶段指标函数 vk(Sk,xk) :从k阶段状态 sk 出发,选择决策 xk 所产生的第k阶段指标
- 过程指标函数 Vk(Sk,xk,xk+1,...,xn):从k阶段状态 Sk 出发,选择决策 xk,xk+1,...,xn 所产生的过程指标
- Vk(sk,xk,xk+1,...,xn)=V[vk(sk,xk),Vk+1(Sk+1,xx+1,...,xn)]
-
最优值函数 fx(Sk)=min/maxxk∈Dk(Sk)Vk(sk,xk,xk+1,...,xn)
|
连和形式 |
连乘形式 |
指标函数 |
Vk(Sk,xk,xk+1,...,xn)=∑j=kn−1vj(sj,xj)+Vn |
Vk(Sk,xk,xk+1,...,xn)=∏j=kn−1vj(sj,xj)+Vn |
最优值函数 |
fx(Sk)=min/maxxk∈Dk(Sk)Vk(sk,xk,xk+1,...,xn) |
|
动态规划模型 |
|
|
基本原理
- 最优性原理:一个最优策略的子策略也是最优的。
例如:已知最优策略为 A→B1→C1→D1→E,则 C→E 的最优策略为 C1→D1→E 。
- 无后效性原理:如果某阶段的状态给定之后,在此阶段以后过程的发展不受这个阶段以前各个状态的影响
我们以上图为例:
阶段 k:k=1,2,3,4
状态变量 Sk:S1={A}S1={A}S2={B1,B2,B3}S3={C1,C2}S4={D1,D2,D3}S5={E}
决策变量 xk :Sk 状态时选择到达下一阶段的点
允许决策集合:xk∈Dk(Sk)
状态指标函数 vk(Sk,xk):状态点 Sk 到决策点 xk 间的最短距离
最优指标函数 fk(Sk):k阶段状态为 Sk 时到终点 E 间的最短距离
递推方程:{fx(Sk)=minxk∈Dk(Sk){vk(sk,xk)+fk+1(Sk+1)},k=1,2,3,4f5(S5)=0