一、最优化问题与预备知识
1. 预备知识
范数
向量范数满足如下性质:
- 非负性:当α=0,∣∣α∣∣>0,当且仅当α=0时∣∣α∣∣=0.
- 齐次性:∣∣kα∣∣=∣k∣ ∣∣α∣∣,k为任意数.
- 三角不等式:任取α,β∈Rn,都有∣∣α+β∣∣≤∣∣α∣∣+∣∣β∣∣.
常见向量范数:
- 1-范数:∣∣α∣∣1=∑i=1n∣ai∣.
- 2-范数:∣∣α∣∣2=(∑i=1n∣ai∣2)21=αHα,又称为欧氏范数.
- ∞-范数:∣∣α∣∣∞=1≤i≤nmax∣ai∣.
- p-范数:∣∣α∣∣p=(∑i=1n∣ai∣p)p1,p∈[1,∞).
- 椭球范数:∣∣α∣∣A=αTAα,∀x∈Rn,A为正定矩阵.
这些向量范数之间的关系:∣∣x∣∣∞≤∣∣x∣∣2≤∣∣x∣∣1≤n∣∣x∣∣∞.
向量范数的等价性:设∣∣⋅∣∣a、∣∣⋅∣∣b是定义在Rn上的两种向量范数,那么总存在两个正数d1、d2使得对任意x∈Rn都有d1∣∣x∣∣b≤∣∣x∣∣a≤d2∣∣x∣∣b.
矩阵范数满足如下性质:
- 非负性:当A=0,∣∣A∣∣>0,当且仅当A=0时∣∣A∣∣=0.
- 齐次性:∣∣kA∣∣=∣k∣ ∣∣A∣∣,其中k为任意复数.
- 三角不等式:任取A,B∈Rm×n,都有∣∣A+B∣∣≤∣∣A∣∣+∣∣B∣∣.
常见矩阵范数:
- Frobenious范数,∣∣A∣∣F=(∑i=1m∑j=1n∣aij∣2)21=tr(AHA)=∑i=1nλi(AHA),其中λj(AHA)表示λ是AHA的特征值.
- 列和范数:∣∣A∣∣1=jmax(∑i=1m∣aij∣),j=1,2,⋯,n.(即每一列各分量模长和的最大值)
- 谱范数:∣∣A∣∣2=jmax(λj(AHA)).
- 行和范数:∣∣A∣∣∞=imax(∑j=1n∣aij∣),i=1,2,⋯,n.
- 诱导范数:∣∣A∣∣i=X=0max∣∣X∣∣α∣∣AX∣∣α.
范数的相容性:
- 矩阵范数的相容性:∣∣AB∣∣≤∣∣A∣∣ ∣∣B∣∣,其中∣∣⋅∣∣是矩阵范数。
- 矩阵范数与向量范数之间的相容性:∣∣Ax∣∣≤∣∣A∣∣∗∣∣x∣∣,其中∣∣⋅∣∣∗是矩阵范数,∣∣⋅∣∣是向量范数,A是矩阵,x是向量。
常用范数不等式:
Cauchy-Schwarz不等式:∀x,y∈Rn,∣xTy∣≤∣∣x∣∣2∣∣y∣∣2,当且仅当x和y线性相关时等号成立。
广义Cauchy-Schwarz不等式:设A∈Rn×n正定,则∀x,y∈Rn,∣xTy∣≤∣∣x∣∣A∣∣y∣∣A−1.
Young不等式:设p,q>1且p1+q1=1,如果a,b∈R,则∣ab∣≤p∣a∣p+q∣b∣q,当且仅当∣a∣p=∣b∣q时等号成立。
Holder不等式:设x,y∈Rn,则∣xTy∣≤∣∣x∣∣p∣∣y∣∣q=(∑i=1n∣xi∣p)p1(∑i=1n∣y∣q)q1,其中p>1、q>1,并且p1+q1=1.
Minkowski不等式:设x,y∈Rn,p∈[1,∞),则 ∣∣x+y∣∣p≤∣∣x∣∣p+∣∣y∣∣p=(∑i=1n∣xi∣p)p1+(∑i=1n∣yi∣p)p1.
函数可微性
连续:若给定函数f:F⊆Rn,函数f在F上连续⇔f在每一点x∈F连续;记作f∈C.
连续可微:若F为开集且在每一点x∈F处,f在F上连续可微⇔每一个偏导数∂xi∂f(x)(i=1,2,⋯,n)存在且连续;记作f∈C1.
二次连续可微:若在每一点x∈F处,f在F上二次连续可微⇔每一个二阶偏导数∂xi∂xj∂2f(x)(i,j=1,2,⋯,n)存在且连续;记作f∈C2.
梯度、海塞阵(Hesse)、Jacobi阵
梯度:设f:F⊆Rn→R是一阶连续可微的,则f在x处的一阶偏导数(即f在x处的梯度)为∇f(x)=(∂x1∂f(x),∂x2∂f(x),⋯,∂xn∂f(x))T.
海塞阵(Hesse):设f是二阶连续可微的,则f在x处的二阶导数(即f在x处的Hesse阵)为H=∇2f(x)=∂x12∂2f∂x2∂x1∂2f⋮∂xn∂x1∂2f∂x1∂x2∂2f∂x22∂2f⋮∂xn∂x2∂2f⋯⋯⋱⋯∂x1∂xn∂2f∂x2∂xn∂2f⋮∂xn2∂2f=(∂xi∂xj∂2f)n×n.(由于f二阶连续可微,所以混合偏导与顺序无关,该矩阵是对称的)
比如:设A对称,二次函数f(x)=21xTAx+bTx+c,其中A∈Rn×n,b∈Rn,c∈R,则∇f(x)=Ax+b,∇2f(x)=A.
多变量向量值函数的Jacobi阵:设多变量向量值函数F:F⊆Rn→Rm在x∈F处连续可微,则F在x∈F处的一阶导数为J=F′(x)=∂x1∂F1(x)∂x1∂F2(x)⋮∂x1∂Fm(x)∂x2∂F1(x)∂x2∂F2(x)⋮∂x2∂Fm(x)∂x3∂F1(x)∂x3∂F2(x)⋮∂x3∂Fm(x)⋯⋯⋱⋯∂xn∂F1(x)∂xn∂F2(x)⋮∂xn∂Fm(x),称为F在x处的Jacobi阵。
比如:设多变量向量值函数F(x)=Ax,其中A∈Rm×n,则Jacobi阵为F′(x)=A.
例1:求函数f(x)=x12+2x22−2x1x2−4x1+3的梯度与Hesse阵
解:方法一,直接求导写在对应位置上:∇f(x)=(2x1−2x2−44x2−2x1),∂2f(x)=(2−2−24).
方法二:f(x)=21xT(2−2−24)x+(−40)x+3,∇f(x)=Ax+b=(2−2−24)x+(−40)x=(2x1−2x2−44x2−2x1),∂2f(x)=A=(2−2−24).
例2:求向量值函数F(x)=(ex1x2+3sinx2,1+x22cosx1)T的Jacobi阵
解:F′(x)=(∂x1∂F1(x)∂x1∂F2(x)∂x2∂F1(x)∂x2∂F2(x))=(x2ex1x2−x22sinx1x1ex1x2+3cosx22x2cosx1).
多变量实值函数的中值定理、泰勒公式
定理:设f:F⊆Rn→R,且x,x∗∈F,如果f是一阶连续可微的,则
(1)存在α∈(0,1)使得f(x)=f(x∗)+∇f(ξ)T(x−x∗),其中ξ=x∗+α(x−x∗).
(2)f在x∗处有一阶Taylor公式:f(x)=f(x∗)+∇f(x∗)T(x−x∗)+o(∣∣x−x∗∣∣).
如果f在F上是二阶连续可微的,则
(3)存在α∈(0,1)使得f(x)=f(x∗)+∇f(x∗)T(x−x∗)+21(x−x∗)T∇f2(ξ)(x−x∗),其中ξ=x∗+α(x−x∗).
(4)f在x∗处有二阶Taylor公式:f(x)=f(x∗)+∇f(x∗)T(x−x∗)+21(x−x∗)T∇f2(x∗)(x−x∗)+o(∣∣x−x∗∣∣2).
2. 凸集
凸集:给定非空集合F⊆Rn,如果∀x,y∈F,α∈[0,1],都有αx+(1−α)y∈F,那么称F为Rn中的一个凸集,如果凸集为开集,则称为开凸集;若凸集为闭集,则称为闭凸集。并且我们规定,空集为凸集。
如果α给定,那么αx+(1−α)y就是一条线,因此有凸集的另外一种定义:F为凸集⇔F中任意两点的连线段仍然属于F。
易证,若ω∈Rn\{0}、β∈R,下面的集合都是凸集
- 单点集、Rn.
- 超平面:H={x∈Rn∣ωTx=β}.
- 闭半空间:{x∈Rn∣ωTx≥β}和{x∈Rn∣ωTx≤β}.
- 开半空间:{x∈Rn∣ωTx>β}和{x∈Rn∣ωTx<β}.
- 超球:{x∈Rn∣∥x∥≤β},其中β≥0.
命题:假设Fi⊆Rn为凸集、βi∈R(i∈{1,2,⋯,p}),则下列集合均为凸集
- 交集F:=F1∩F2∩⋯∩Fp.(有限个凸集的交仍为凸集)
- 集合β1F1:={β1x∣x∈F1}.(凸集的数乘仍为凸集)
- 和集F1+F2:={x+y∣x∈F1,y∈F2}.(凸集的和仍为凸集)
- 集合∑i=1pβiFi.(凸集数乘后的和仍为凸集)
定理:非空集合F⊆Rn为凸集⇔∀xi∈F及任意满足∑i=1pαi=1的非负实数αi(i∈{1,2,⋯,p}且p≥2为正整数)都有∑i=1pαixi∈F.
证明:根据凸集定义可得定理的充分性,下面用归纳法证明必要性
p=2时由定义知结论成立;假设当p=k时结论成立,下面证明p=k+1时结论也成立:
对∀xi∈F,∑i=1k+1αi=1,则∑i=1k+1αixi=∑i=1kαixi+αk+1xk+1=(1−αk+1)(∑i=1k1−αk+1αixi)+αk+1xk+1,由∑i=1k1−αk+1αi=1且1−αk+1αi≥0知,归纳假设∑i=1k1−αk+1αixi∈F,于是上式可看作(1−αk+1)x+αk+1y,根据p=2时结论成立知,(1−αk+1)x+αk+1y∈F,即∑i=1k+1αixi∈F,由归纳原理知结论成立。
集合的(严格)分离:设F1,F2⊆Rn为两个非空凸集,如果存在非零向量ω∈Rn和实数t,使得
(1)对任意x∈F1、y∈F2都有ωTx≥t且ωTy≤t,则称超平面π={x∈Rn∣ωTx=t}分离集合F1、F2.
(2)对任意x∈F1、y∈F2都有ωTx>t且ωTy<t,则称超平面π={x∈Rn∣ωTx=t}严格分离集合F1、F2.
点与凸集分离定理:设F为Rn中的非空闭凸集,x0∈Rn且x0∈/F,则存在Rn中的超平面严格分离集合F与{x0}.
3. 凸函数
Farkas引理:设A∈Rm×n、b∈Rn,则不等式组(1):{Ax≤0bTx>0、(2):{ATy=by≥0有且仅有一组有解(即(1)和(2)不能同时有解)。
证明:假设(2)有解,即∃y∈Rn使得ATy=b,y≥0,若(1)也有解则∃x∈Rn使得Ax≤0,此时bTx=(ATy)Tx=yTAx≤0
假设(2)无解,记Ω={z:ATy∣y≥0},则Ω⊆Rn为非空闭凸集,且b∈/Ω,根据点与凸集分离定理知,存在ω∈Rn、t∈R使得对所有z∈Ω有ωTb>t、ωTz<t,并且由0∈Ω知t>0,t>ωTz=ωT(ATy)=yTAω,y≥0,由y的任意性知,Aω≤0,bTω>0,从而ω为(1)的解。
定理:设p、q是两个非负整数,u0,u1,⋯,up,v1,v2,⋯,vq∈Rn,则等式与不等式组(1):⎩⎨⎧dTu0<0dTui=0dTvi≥0,i∈{1,2,⋯,p},i∈{1,2,⋯,q}无解⇔存在实数αi(i∈{1,2,⋯,p})和非负实数βi(i∈{1,2,⋯,q})使得(2):u0=∑i=1pαiui+∑i=1qβivi.
凸函数
凸函数的定义:设F⊆Rn为非空凸集,给定函数f:F→R:
(1)若对任意的x,y∈F及任意的α∈[0,1],有f(αx+(1−α)y)≤αf(x)+(1−α)f(y),则称函数f为凸集F上的凸函数。
x与y凸组合的函数小于等于x与y的函数的凸组合叫凸函数。
(2)若对任意的x,y∈F且x=y,及任意的α∈(0,1),有f(αx+(1−α)y)<αf(x)+(1−α)f(y),则称函数f为凸集F上的严格凸函数。
(3)若存在常数c>0,使得对任意的x,y∈F及任意的α∈(0,1),有f(αx+(1−α)y)≤αf(x)+(1−α)f(y)−cα(1−α)∥x−y∥2,则称函数f为凸集F上的强凸函数(或一致凸函数)。
对应凹函数的定义:
(1)f(αx+(1−α)y)≥αf(x)+(1−α)f(y)-凹函数。
(2)f(αx+(1−α)y)>αf(x)+(1−α)f(y)-严格凹函数。
(3)f(αx+(1−α)y)≥αf(x)+(1−α)f(y)+cα(1−α)∥x−y∥2-强凹函数。
例如,给定向量c∈Rn,则线性函数f(x)=cTx既是凸函数又是凹函数,因为此时f(αx+(1−α)y)=αf(x)+(1−α)f(y).
凸函数的性质(设F⊆Rn是凸集):
- 函数f是F上的(严格)凸函数⇔−f是F上的(严格)凹函数。
- 设函数f1、f2是F上的凸函数,实数α1,α2≥0,则函数α1f1+α2f2也是F上的凸函数。
- 设函数f是F上的凸函数,t为实数,则水平集L1(f)={x∈F:f(x)≤t}是凸集。
一阶判别定理(凸函数的判别准则1):设函数f在凸集F⊆Rn上可微,则:
(1)f在F上为凸函数⇔对任意的x,y∈F,有f(y)−f(x)≥∇f(x)T(y−x).
(2)f在F上为严格凸函数⇔对任意不同的x,y∈F,有f(y)−f(x)>∇f(x)T(y−x).
第一条的证明:
必要性⇒:由于f是F上的凸函数,则∀x,y∈F、∀α∈[0,1],有f(αy+(1−α)x)≤αf(y)+(1−α)f(x),整理得f(x+α(y−x))≤f(x)+α[f(y)−f(x)],由一阶泰勒展开得f(x+α(y−x))=f(x)+α∇f(x)T(y−x)+o(∥α(y−x)∥),所以f(y)−f(x)≥∇f(x)T(y−x)+αo(∥α(y−x)∥),令α→0得f(y)−f(x)≥∇f(x)T(y−x).
充分性⇐:∀x,y∈F、∀α∈[0,1],有z=αx+(1−α)y∈F,于是{f(x)−f(z)≥∇f(z)T(x−z)f(y)−f(z)≥∇f(z)T(y−z),两式分别乘以α、(1−α)再相加,得αf(x)+(1−α)f(y)−f(z)≥0,即f(αx+(1−α)y)≤αf(x)+(1−α)f(y),因此f是F上的凸函数。
第二条的证明:
必要性⇒:∀x,y∈F且x=y,令z=21(x+y)此时z∈F,根据第一条的必要性证明可以得到f(z)−f(x)≥∇f(x)T(z−x),并且根据凸函数定义有f(z)=f(21x+21y)<21f(x)+21f(y),所以21f(x)+21f(y)>f(x)+∇f(x)T(z−x),整理得f(y)−f(x)>∇f(x)T(y−x).
充分性⇐:证明和第一条类似。
二阶判别定理(凸函数的判别准则2):设在开凸集F⊆Rn内函数f二阶可微,则:
(1)f在F内为凸函数⇔对任意的x∈F,∇2f(x)半正定.
(2)若∀x∈F,∇2f(x)正定,则f在F内为严格凸函数。
第一条的证明:
必要性⇒:∀x∈F,0=y∈Rn,由F是开集知,存在ϵ>0使得∀α∈(−ϵ,ϵ),x+αy∈F,由一阶判别定理有f(x+αy)≥f(x)+α∇f(x)Ty,由Taylor展开式得f(x+αy)=f(x)+α∇f(x)Ty+21α2yT∇2f(x)y+o(∥αy∥2),所以yT∇2f(x)y+2α2o(∥αy∥2)≥0,令α→0得yT∇2f(x)y≥0,由y的任意性得∇2f(x)半正定。
充分性⇐:∀x,y∈F及∇2f(x)对称半正定,对f(y)做代拉格朗日余项的Taylor展开,并由于XT∇2f(ξ)X≥0得f(y)=f(x)+∇f(x)T(y−x)+21(y−x)T∇2f(ξ)(y−x)≥f(x)+∇f(x)T(y−x),其中ξ=x+t(y−x),t∈(0,1),因此f是F上的凸函数。
强凸函数的判别定理:设f:Rn→R二次连续可微,则:f是强凸函数⇔∇2f(x)一致正定。(注:一致正定即存在c>0对所有x,d∈Rn都有dT∇2f(x)d≥c∥d∥,也就是∇2f(x)的最小特征值是以c为下界的)
4. 凸规划
凸规划问题:设F⊆Rn为凸集,f:F→R为凸函数,则称x∈Fminf(x)为凸规划问题。
凸规划的性质:
- 凸规划问题的任一局部最优解x为其全局最优解.
- 凸规划问题的最优解集S为凸集.
- 若函数f为非空凸集F上的严格凸函数,且凸规划问题存在全局最优解,则其全局最优解唯一.
第一点的证明(反证法):假设x∗是凸规划问题的局部最优解但非全局最优解,则至少存在一个y∗∈F使得f(y∗)<f(x∗),由函数f为凸函数、F为凸集得,对任意α∈[0,1]有{αy∗+(1−α)x∗∈Ff(αy∗+(1−α)x∗)≤αf(y∗)+(1−α)f(x∗)<f(x∗),当α→0+时,αy∗+(1−α)x∗→x∗,因此在充分接近x∗时,f(αy∗+(1−α)x∗)<f(x∗),这与x∗是局部最优解矛盾。
第二点的证明:由空集与单元素集为凸集,不妨设凸规划问题的最优解集S至少含两个元素,假设x∗,y∗∈S,则∀α∈[0,1]有x∗,y∗∈F且f(αx∗+(1−α)y∗)≤αf(x∗)+(1−α)f(y∗)=f(x∗)=x∈Fminf(x)(因为x∗、y∗都是最优解,所以f(x∗)=f(y∗)),由于x∈Fminf(x)是最优值,根据定义它是最小值,因此小于等于只能取等于号,即αx∗+(1−α)y8∈S.
第三点的证明(反证法):假设全局最优解不唯一,则∃x∗,y∗∈S且x∗=y∗,显然x∗,y8∈F且f(x∗)=f(y∗),∀α∈[0,1],αx+(1−α)y∗∈S且f(αx∗+(1−α)y∗)<αf(x∗)+(1−α)f(y∗)=f(x∗),这与x∗是全局最优解矛盾,所以凸规划问题有唯一最优解。
5. 线搜索迭代算法
线搜索算法的一般框架
基本思路:从某个点x0∈Rn出发,按照某种规则产生一个迭代点序列{xk},直到算法终止。
算法分类:
(1)按照迭代点是否可行,分为 ①可行算法-所有迭代点都是可行点
②不可行算法-迭代点中至少存在一个不可行点
(2)按照目标函数值是否下降分为 ①单调下降算法f(xk+1)≤f(xk),∀k
②非单调下降算法
一般框架:
步1:选择一个初始点x0∈Rn,令k=0.
步2:判断当前的迭代点xk是否满足终止条件.
步3:从当前点出发,选择沿什么方向进行迭代(记迭代方向为dk)以及沿该方向走多远(记迭代步长为λk),确定下一个迭代点xk+1=xk+λkdk.
步4:令k=k+1,转步2.
注1:关于终止规则
最理想的终止规则:∣f(xk)−f(x∗)∣≤ϵ;∥xk−x∗∥≤ϵ,x∗∈S.
常用的终止规则:
(1)∥xk+1−xk∥<ϵ (2)∣f(xk)−f(x∗)∣<ϵ (3)∥xk∥∥xk+1−xk∥<ϵ (4)∣f(xk)∣∣f(xk+1)−f(xk)∣<ϵ (5)∥∇f(xk)∥≤ϵ
注2:迭代方向
迭代方向定义:在点xk处,对于向量dk∈Rn\{0},若存在λ>0使得对任意λ∈(0,λ)都有f(xk+λdk)<f(xk),则称dk为函数f在xk处的一个下降方向。
命题:假设函数f一阶连续可微,那么dk为f在xk处的下降方向当且仅当∇f(xk)Tdk<0.
证明:由Taylor展开式得f(xk+λdk)=f(xk)+λ∇f(xk)Tdk+o(∥dk∥),所以f(xk+λdk)−f(xk)=λ[∇f(xk)Tdk+o(∥dk∥)],从而结论成立。
注3:可行方向
可行方向定义:给定非空集合F⊆Rn和点xk∈F,对于向量dk∈Rn\{0},若存在λ>0使得对任意λ∈(0,λ)都有xk+λdk∈F,则称dk为xk处关于F的一个可行方向。
可行下降方向定义:如果dk为点xk处的可行下降方向,那么存在λk使得xk+λkdk∈F,且f(xk+λkdk)<f(xk).
注4:迭代步长的选取
1、精确一维线搜索:选取λk=argλ>0minφ(λ)=λ>0minf(xk+λdk),f(xk+1)=f(xk+λkdk).
2、非精确一维线搜索
(1)Goldstein型线搜索(1965):选取λk>0使得{f(xk+λkdk)−f(xk)≤σλk∇f(xk)Tdkf(xk+λkdk)−f(xk)≥(1−σ)λk∇f(xk)Tdk,其中σ∈(0,21).
(2)Armijo型线搜索(1966):选取λk=ργmk使得mk为满足下式的最小非负整数:f(xk+ργmkdk)−f(xk)≤σργmk∇f(xk)Tdk,其中ρ>0,γ∈(0,1).
(3)Wolfe型线搜索:选取λk>0使得{f(xk+λkdk)−f(xk)≤σλk∇f(xk)Tdk∇f(xk+λkdk)Tdk≥δ∇f(xk)Tdk,其中σ∈(0,21),δ∈(σ,1).
3、非单调一维线搜索
(1)Grippo-Lampariello-Lucidi非单调线搜索(1986):寻找步长λk=ργhk使得hk是满足下式的最小非负整数:f(xk+λkdk)≤0≤j≤mkmaxf(xk−j)+δλk∇f(xk)Tdk,其中mk是一个正整数且mk≤M(M为某一正整数),λk0是一个给定的小实数,σ,ρ∈(0,1).
注:当f(xk)=0≤j≤mkmaxf(xk−j)时,即为单调的Armijo型线搜索。
(2)Zhang-Hager非单调线搜索(2004):寻找步长λk>0使得{f(xk+λkdk)≤Ck+δλk∇f(xk)Tdk∇f(xk+λkdk)Tdk≥σ∇f(xk)Tdk,其中0<δ<σ<1,Ck按如下方式选取:Ck+1=(ηkQkCk+f(xk+1))\Qk+1,这里Qk+1=ηkQk+1,C0=f(x0),Q0=1,且ηk∈[ηmin,ηmax],0≤ηmin≤ηmax≤1.
注:Ck是f(x0),f(x1),⋯,f(xk)的凸组合,含有k步迭代之前所有迭代点函数值的信息;ηk的选取控制着”非单调性”的度,如果ηk=0,k∈{0,1,2,⋯},则即为单调的Wolfe搜索。
(3)Hu-Huang-Lu非单调线搜索(2010):寻找步长λk>0使得{f(xk+λkdk)≤Ck+δλk∇f(xk)Tdk∇f(xk+λkdk)Tdk≥σ∇f(xk)Tdk,其中0<δ<σ<1,Ck按如下方式选取:Ck+1=(ηk∑l=1mk−1ηk−lf(xk−l)+f(xk))\Qk,且Qk=ηk∑l=1mk−1ηk−l+1,其中ηk∈[ηmin,ηmax],0≤ηmin≤ηmax≤1,f(xk)−Ck≤2δlk∣∇f(xk)Tdk∣,这里,若∇f是Lipschitz连续的,则令lk=L∥dk∥2(1−σ)∣∇f(xk)Tdk∣,其中L为Lipschitz常数;否则,令lk=0.
注一个特例:若∑l=1mk−1ηk−lf(xk−l)≥∑l=1mk−1ηk−lf(xk),则ηk∈(0,1];否则,选择ηk=0。因此,每一步要么使用单调线搜索,要么使用非单调线搜索,是一种混合线搜索方法。
算法收敛性
1、适定性(well-defined):如果算法的每一步是适定的,则这个算法是适定的。
2、算法是收敛的:如果算法是有限终止的或所产生迭代点列{xk}的每个聚点是优化问题的最优解,则称该算法是收敛的。
3、全局收敛(globally convergent):局部收敛(locally
convergent):如果对于任意的初始点x0算法是收敛的,则称该算法是全局收敛的;如果只有初始点x0充分靠近最优解时,算法是收敛的,则称该算法是局部收敛的。
4、收敛率(rate of convergence):假设算法生产无穷点列{xk},且算法收敛到最优解x∗,不妨设{xk}收敛到x∗,即k→∞lim∥xk−x∗∥=0,则:
(1)全局Q-线性收敛:如果算法初始迭代点的选取无关于最优解的信息,且存在β∈(0,1)使得∥xk+1−x∗∥≤β∥xk−x∗∥,∀k∈{0,1,2,⋯},则称该算法是全局Q-线性收敛的。
(2)局部Q-收敛率:如果存在β>0、ζ≥0使得k→∞lim∥xk−x∗∥β∥xk+1−x∗∥=ζ,则称该算法是局部Q-β阶收敛的。
注:特别地,
若β=1,ζ∈(0,1),则称该算法是局部Q-线性收敛的.
若β∈(1,2)、ζ>0或β=1、ζ=0,则称该算法是局部Q-超线性收敛的.
若β=2,ζ>0,则称算法是Q-二阶收敛的.
二、线性规划问题
1. 线性规划问题的数学模型
例1:某制药厂生产甲、乙两种药品,生产这两种药品都有消耗某种原料,生产每吨药品所需的原料量及所占设备时间见下表
项目\消耗量 |
每吨产品的消耗 |
|
每周资源总量 |
|
甲 |
乙 |
|
原料/kg |
30 |
20 |
160 |
设备/台班 |
5 |
1 |
15 |
该厂每周所能得到的原料量为160kg,每周设备最多能开15个台班,且根据市场需求,甲种产品每周产量不应超过4t。已知该厂生产每吨甲、乙两种产品的利润分别为5万元和2万元,问该厂应如何安排两种产品的产量才能使每周获得的利润最大?
解:设该厂每周安排生产甲、乙两种药品的产量分别为x1、x2,则x1,x2需要满足不等式条件s.t.⎩⎨⎧30x1+20x2≤1605x1+x2≤15x1≤4x1≥0x2≥0,并且该问题变为maxz=5x1+2x2.
例2:某铁器加工厂需要制作100套钢架,每套要用长分别为2.9m、2.1m和1.5m的圆钢各一根,已知原料长为7.4m。问应如何下料,可使用材料最省?
解:通过分析,可知有如下几种切割方式
长度(数量)\方案编号 |
Ⅰ |
Ⅱ |
Ⅲ |
Ⅳ |
Ⅴ |
2.9m |
1 |
2 |
0 |
1 |
0 |
2.1m |
0 |
0 |
2 |
2 |
1 |
1.5m |
3 |
1 |
2 |
0 |
3 |
料头/m |
0 |
0.1 |
0.2 |
0.3 |
0.8 |
假设按第i种方案下料的原料根数为xi(i=1,2,3,4,5),则求minz=0x1+0.1x2+0.2x3+0.3x4+0.8x5,s.t.⎩⎨⎧x1+2x2+x4=1002x3+2x4+x5=1003x1+x2+2x3+3x5=100x1,x2,x3,x4,x5≥0且为整数.
线性规划问题的数学模型:像上面那样,每个问题都有一组决策变量x1,x2,⋯,每个问题都有这些决策变量需要满足的一组约束条件(一般为线性等式或不等式),每个问题都有一个关于决策变量的线性目标函数并且该函数要在满足约束条件的条件下实现最大化或最小化,那么我们将约束条件和目标函数都是决策变量的线性函数的规划问题称为线性规划。
2. 图解法(两个决策变量)
对于只有两个决策变量的线性规划问题,一般采用图解法
例:求解maxz=5x1+2x2,s.t.⎩⎨⎧30x1+20x2≤1605x1+x2≤15x1≤4x1≥0,x2≥0.
解:以x1、x2为坐标轴建立直角坐标系,分别作出约束条件所对应的半平面(图中的蓝色区域就是可行域)
两个决策变量的线性规划问题-例
求出z=5x1+2x2的梯度t=(∂x1∂z,∂x2∂z)T=(5,2)T,将任一等值线l沿t的方向平移,直到将离开还未离开可行域时的一根等值线为l∗,即为过图中A(2,5)点的等值线,故(2,5)为唯一的最优解。
关于可行域与解的一些结论
- 若可行域非空且有界,则线性规划问题必有最优解;若可行域无界,则线性规划问题可能有最优解,也可能无最优解。
- 若线性规划问题有最优解,其最优解必在某个顶点处达到,最优解的个数或是唯一的,或有无穷多个。
3. 线性规划问题的标准化
标准的线性规划问题要求满足:
- 目标函数一律改为最大化.
- 约束条件一律改为等式(非负约束条件除外).
- 等式右端项大于等于零.
线性规划问题的矩阵形式:maxz=cx,s.t.{Ax=bx≥0,其中c=(c1,⋯,cn),x=x1⋮xn,b=b1⋮bm,A=a11a21⋮am1a12a22⋮am2⋯⋯⋯a1na2n⋮amn.
线性规划问题的向量形式:maxz=cx,s.t.{∑j=1nxjpj=bx≥0,其中c,x,b同矩阵形式,而pj=a1j⋮amj,j=1,2,⋯,n,即A=(p1,⋯,pn).
如何将线性规划问题改写为标准形式:
(1)若目标函数为最小化minz=cx,则作函数z′=−cx,对z′实现最大化,即maxz′=−cx.
(2)若约束条件是小于等于型,则在该约束不等式左边加上一个松弛变量,将不等式改为等式.
比如:x1−2x2+3x3≤8⇒x1−2x2+3x3+x4=8.
(3)若约束条件是大于等于型,则在该约束不等式左边减去一个剩余变量,将不等式改为等式。
比如:2x1−3x2−4x3≥5⇒2x1−3x2−4x3−x4=5.
(4)若某个约束方程右端项bi<0,则在约束方程两端乘以−1,不等式改变方向,然后再将不等式转化为等式.
(5)若决策变量xk无非负要求,则可令两个新变量xk′≥0和xk′′≥0,作xk=xk′−xk′′,在原有的数学模型中,xk均用xk′−xk′′来替代,而在非负约束中,增加xk′≥0、xk′′≥0.
即一般的线性规划问题都要求决策变量x≥0,如果实际问题中x可以小于零,就要转化为两个非负的新的变量相减。
例:将下列线性规划问题化为标准形式:minz=x1−2x2+3x3,s.t.⎩⎨⎧x1+x2+x3≤7x1−x2+x3≥2−3x1+x2+2x3=−5x1≤0,x2≥0,x3无约束
解:令z′=−z,再令x3=x4−x5,并在第一个约束条件中加入松弛变量x6,在第二个约束条件左边减去剩余变量x7,第三个约束条件两边同乘−1,得maxz′=−x1+2x2−3(x4−x5)+0⋅x6+0⋅x7,s.t.⎩⎨⎧x1+x2+x4−x5+x6=7x1−x2+x4−x5−x7=23x1−x2−2x4+2x5=5x1,x2,x3,x4,x5,x6,x7≥0.
4.线性规划问题的解
可行解、最优解:若数学模型为maxz=cx,s.t.{Ax=bx≥0,其中A=(aij)m×n,r(A)=m<n,x∈Rn,b∈Rm,b≥0,c∈Rn,我们把凡是满足Ax=b及x≥0的解x=x1⋮xn称为线性规划问题的可行解,同时满足maxz=cx的可行解为最优解。
基矩阵、基向量:设线性规划问题约束方程组的系数矩阵Am×n的秩为m,则A中某m列组成的任一m阶可逆阵B称为该规划问题的基矩阵,简称基,若把基记为B=(p1,p2,⋯,pm),则称pk为基B中的一个基向量,而A中其余n−m个列向量称为非基向量。
基变量:当确定了Ax=b的一个基B后,与基向量pk相对应的决策变量xk称为关于基B的一个基变量,而与非基向量对应的决策变量称为非基变量。
基本解、基本可行解:若确定了基B及其对应的基变量xk1,⋯,xkm,我们称非基变量取值均为零且满足约束条件的一个解x,为关于基B的一个基本解;如果该基本解同时满足非负条件x≥0,则称基本解为基本可行解。
基本解的确定:设B为一个基矩阵,xB为对应的基变量,N为非基矩阵,xN为对应的非基变量,那么Ax=b可以改写为BxB+NxN=b,因此xB=B−1b−B−1NxN,如果令非基变量xN=0,则基变量xB=B−1b,故[xBxN]=[B−1b0]就是关于基B的一个基本解。
例:求{x1+2x2≤8x2≤2的基本解和基本可行解
解:化为标准形式s.t.{x1+2x2+x3=8x2+x4=2,于是A=[10211001]=(p1,p2,p3,p4),设B1=(p1,p2)=[1021]是一个基,且x1,x2是相应的基变量,而x3,x4是非基变量,令x3=x4=0,算出基本解为x1=4200,又因为x1≥0,那么x1也是基本可行解。
另外,如果选取B2=(p1,p3),得到x2=8002,此时x3也是基本可行解。其他选法类似。
5. 线性规划问题基本理论
引入:通过两个变量的线性规划问题的图解法,可以得到如下结论:
- 线性规划问题的可行域是一个有界或无界的凸多边形,其顶点的个数是有限个。
- 若线性规划问题有最优解,则其最优解必可在某顶点处达到。
下面考虑对于多变量线性规划问题,讨论以上结论是否仍然成立。
凸集:设集合S∈Rn,若对∀x1,x2∈S及每一个数λ∈[0,1],都有λx1+(1−λ)x2∈S成立,则称S为凸集。
从直观上看,没有凹入部分,或没有空洞的是凸集。若集合S中任意两点连线上的每一点仍在S中,则S为凸集。
凸组合:设x1,x2,⋯,xk是n维欧式空间Rn中的k个点,若存在实数λ1,λ2,⋯,λk,其中λi∈[0,1]且∑i=1kλi=1,有x=λ1x1+λ2x2+⋯+λkxk,则称x为x1,x2,⋯,xk的一个凸组合。
严格凸组合:若存在实数λ1,λ2,⋯,λk,其中λi∈(0,1)且∑i=1kλi=1,有x=λ1x1+λ2x2+⋯+λkxk,则称x为x1,x2,⋯,xk的一个严格凸组合。
极点:设S是凸集,x∈S,若x不能表示为S中任意两个不同点x1,x2的一个严格凸组合,则称x为S的一个极点。
若x是S中的一个极点,且有x=λx1+(1−λ)x2,λ∈(0,1),x1,x2∈S,则必有x=x1=x2.
五边形的每个顶点都是极点,而圆域边界上任意一点都是极点。
定理:若线性规划问题存在可行域,则可行域S={x∣Ax=b,x≥0,A∈Mm×n,b∈Rm,x∈Rn}为一凸集。
定理:线性规划问题的可行解x=(x1,x2,⋯,xn)T为基本可行解当且仅当x中正分量所对应的系数列向量线性无关。
定理:线性规划问题的每个基本可行解x对应于可行域S的一个极点。
证明:设基本可行解x的前m个分量为基变量,即∑j=1mpjxj=b,如果x不是可行域S中的极点,则x一定可以表示为可行域S中两个不同的点x(1)和x(2)的一个严格凸组合,即x=λx(1)+(1−λ)x(2),0<λ<1,若记x(1)=(x1(1),x2(1),⋯,xn(2))T,x(2)=(x1(2),x2(2),⋯,xn(2))T,则有x1⋮xm0⋮0=λx1(1)⋮xm(1)xm+1(1)⋮xn(1)+(1−λ)x1(2)⋮xm(2)xm+1(2)⋮xn(2),故λxk(1)+(1−λ)xk(2)=0(k=m+1,⋯,n),根据基本可行解的非负条件知,xk(1)=0、xk(2)=0,因为x(1),x(2)∈S是可行解,有{x1(1)p1+x2(1)p2+⋯+xm(1)pm=bx1(2)p1+x2(2)p2+⋯+xm(2)pm=b,两式相减得(x1(1)−x1(2))p1+(x2(1)−x2(2))p2+⋯+(xm(1)−xm(2))pm=0,由于x(1)=x(2),所以至少有一个xj(1)−xj(2)=0,其中1≤j≤m,也就是说p1,p2,⋯,pm线性相关,这与x是基本可行解矛盾,故假设不成立。
定理:有界凸集S上任一点x都可以表示为S的极点的凸组合。
定理:对线性规划问题的标准形式,若可行域有界,且存在最优解,则目标函数必可在其可行域S的某个顶点达到最优值。
定理:对于线性规划问题的标准形式,若存在可行解,则必存在基本可行解。其中约束矩阵为秩是m的m×n型矩阵A。
通过以上分析,可得到以下结论:
- 线性规划问题的可行域是一个凸集,可行域可能有界也可能无界,但其顶点数是有限个。
- 线性规划问题每个基本可行解对应于可行域的一个极点(顶点)。
- 若线性规划问题有最优解,则必可在其可行域的某个(或多个)极点上达到最优值。
于是,我们就可以用代数方法寻找最优解:求最优解先找基本可行解,而基本可行解的个数不会超过基本解的个数Cnm。
6.单纯形法
单纯形法概述
引入:前面我们已经知道了要求线性规划问题的最优解可以先找基本可行解再逐个比较,但是往往基本可行解的个数Cnm会随着m、n增大而迅速增大,所以这种方法往往行不通。于是我们设想如果可以从某一基本可行解(初始基本可行解)出发,每次寻求比上次更“好”的基本可行解,这样就大大减少了计算量。
要使用这种逐步改善的求解方法,要解决如下几个问题:(1)如何得到一个初始基本解
(2)如何判别是否达到了最优解
(3)若当前不是最优解,如何寻找一个更好的基本可行解
下面利用例题来说明单纯性法的思路:
例:maxz=5x1+2x2,s.t.⎩⎨⎧30x1+20x2≤1605x1+x2≤15x1≤4x1,x2≥0.
解: 将问题标准化得,maxz=5x1+2x2+0x3+0x4+0x5,s.t.⎩⎨⎧30x1+20x2+x3=1605x1+x2+x4=15x1+x5=4x1,x2,x3,x4,x5≥0,写出约束矩阵为A=30512010100010001。
步骤一:选取一个可逆的子阵作为初始基矩阵,比如B(0)=100010001,则相应的基变量为xB=x3x4x5,于是初始基本可行解为[xBxN]=[B−1b0],带入计算得初始基本可行解为(0,0,160,15,4)T,相应的目标函数值为maxz=0。
步骤二:判断当前解是否为最优解。将基变量用非基变量表示⎩⎨⎧x3=160−30x1−20x2x4=15−5x1−x2x5=4−x1,用非基变量表示目标函数z=0+5x1+2x2,从数学角度看,目标函数中非基变量x1,x2的系数为正数,故若让非基变量x1(或x2)的取值从零增加,相应的目标函数值z也将增加,因此就有可能找到一个新的基本可行解,使目标函数值比x(0)的更“好”。因此,这个不是最优解。
步骤三:解的改进。在z=5x1+2x2中,x1前的系数必比x2前的系数大,即x1每增加一个单位对z的贡献比x2大,故让x1的取值从零变成正值。这样,x1从非基变量转为基变量,我们称之为进基变量。但对于本例,任意一个基本可行解中只能有3个基变量,因此必须从原有的基变量x3,x4,x5中选择一个离开基变量,我们称之为离基变量(或出基变量),为了选取合适的离基变量,观察⎩⎨⎧x3=160−30x1−20x2x4=15−5x1−x2x5=4−x1,此时我们已经选取了x1为进基变量,而x2仍为非基变量取值仍为0,我们选取在x1从零增加直到使x3,x4,x5的取值减少到零时停止,第一个变为零的基变量为离基变量。显然x1=min{30160,515,14}=515=3,因此x4为离基变量(即令x4=0),将x1=3代入前面等式,可得到新的基本可行解x(1)。
经过以上三步,我们得到了一个新的基本可行解x(1)=(3,0,70,0,1)T,然后再重复上述步骤。(1)算出新的基本可行解后再重新写出基矩阵、基变量和目标函数:相应的基矩阵为B(1)=(p3,p1,p5),基变量为x3,x1,x5,非基变量为x2,x4,对应的目标函数值是z(x(1))=15>z(x(0))=0。(2)判断当前解是否为最优解:将基变量用非基变量表示⎩⎨⎧x3=70−14x2+6x4x1=3−51x2−51x4x5=1+51x2+51x4,用非基变量表示目标函数z=15+x2−x4,此时目标函数中x2的系数仍为正数,于是选x2为进基变量,迭代到一个新的基本可行解x(2),就有可能使目标函数值再增加,因此x(1)仍不时最优解。(3)解的改进:在基变量用非基变量表示的式子中,x4仍为非基变量取值为0,在x2从零增加的过程中,x3首先到达零,故取x2=min{1470,513}=1470=5,则x3=0,x1=2>0,x5=2>0。
经过两次迭代,可得到新的基本可行解x(2)=(2,5,0,0,2)T,再次判断此解是否为最优解:将基变量用非基变量表示⎩⎨⎧x2=5−141x3+73x4x1=2+701x3−702x4x5=2−701x3+72x4,用非基变量表示目标函数z=20−141x3−74x4,这时x3,x4中任意一个从零开始增加,都会使z减少,因此当前解x(2)=(2,5,0,0,2)T就是最优解了,此时最优值为z∗=z(x(2))=20。
通过上面例子,总结出单纯性法的求解步骤:
(1)在线性规划问题的标准形式中,设法在约束矩阵Am×n中构造出一个m阶单位矩阵作为初始可行基,获得一个初始基本可行解;
(2)判断当前基本可行解是否为最优解。判断方法为:求出用非基变量表示基变量及目标函数的表示式,称为线性规划问题的典式(或称为规范式),在目标函数的典式中,若至少有一个非基变量前的系数为正数,则当前解不是最优解;而如果所有非基变量的系数均为非正数,则当前解就是最优解。因此我们把目标函数的典式中非基变量前的系数称为检验数,对于最大化问题,当所有的检验数≤0时,当前解为最优解;对于最小化问题,当所有的检验数≥0时,当前解为最优解;
个人理解:随着学习的深入,最优解的意义不同。比如学到对偶单纯形法之后,最大化问题的所有检验数≤0只能叫做正则解,而不能再叫做最优解,因为后面的最优解还要求b≥0,这样才能利用互补松弛性定理写出对偶问题的解,但是如果不考虑对偶问题的话,这里管检验数≤0的解叫最优解也没错,只是学的越后面,对最优解的要求变得越多。
(3)若当前解不是最优解,则要进行基变换迭代到下一个基本可行解。变换方式为:首先从当前解的非基变量中选一个作进基变量(选择的原则一般为目标函数的典式中,最大的正检验数所属的非基变量作为进基变量);再从当前解的基变量中选一个作为离基变量(选择的方法是在用非基变量表示基变量的典式中,除进基变量外,让其余非基变量取值为零,再按最小比值准则确定离基变量)。选好后再继续重复上面步骤,直至找到最优解。
单纯形法的矩阵描述
考虑线性规划问题的标准模型maxz=cx,s.t.{Ax=bx≥0,其中A=(aij)m×n,r(A)=m<n,x∈Rn,b∈Rm,cT∈Rn,在下面的讨论中,假设:(1)此线性规划问题的可行域非空
(2)所有的基本可行解不退化。
步骤一:将A分解为A=(B,N),其中B为基矩阵,N为非基矩阵,相应地,xB为基变量,xN为非基变量,那么解向量可写为x=[xBxN],根据前面的结论知,当前的基本可行解[xBxN]=[B−1b0],相应地,将价值系数向量也分为两块c=(cB,cN);
步骤二:判断当前解是否为最优解,即求检验数。①先用非基变量表示基变量,由Ax=b得BxB+NxN=b,又因为B可逆,则有xB=B−1b−B−1NxN。②再用非基变量表示目标函数,z=cB(B−1b−B−1NxN)+cNxN=cBB−1b+∑k∈N(ck−cBB−1pk)xk,故对于每个非基变量xk,其对应的检验数为σk=ck−cBB−1pk,其中pk是xk在约束矩阵中对应的列向量,而基变量的检验数ck−cBB−1pk都为零;
步骤三:若当前解不是最优解,则要进行基变换迭代到下一个基本可行解。①若某个非基变量的检验数σk>0,则当前解不是最优解,而且σk越大,目标函数值增加越多,选择进基变量xi使σi=maxk∈N{σk}。②选择离基变量,因为xB=B−1b−B−1NxN=B−1b−∑k∈NB−1pkxk,除xi以外,其他的非基变量取值仍为零,所以xB=B−1b−B−1pixi=b−Yixi,其中b=B−1b,Yi=B−1pi都为m维向量,于是xB=xB1xB2⋮xBm=b1b2⋮bm−y1iy2i⋮ymixi,由最小比值准则知,xi=min{ylibl∣yli>0}=yribr,故当前的第r个基标量xBr为离基变量。于是新的基矩阵为B=(pB1,pB2,⋯,pBr−1,pi,pBr+1,⋯,pBm)。
定理:根据上面分析,可以看出,对于最大化线性规划问题,若存在一个检验数σk>0,且yli≤0,l=1,2,⋯,m,即B−1pi≤0,则线性规划问题有无界解(无最优解)。
定理:对于最大化线性规划问题,若所有非基变量的检验数σk≤0,则当前解为最优解;对于最小化线性规划问题,若所有非基变量的检验数σk>0,则当前解为最优解。
例:用单纯形法求解max4x1+x2,s.t.⎩⎨⎧−x1+2x2≤42x1+3x2≤12x1−x2≤3x1,x2≥0.
解:将问题化为标准型得max4x1+x2+0x3+0x4+0x5,s.t.⎩⎨⎧−x1+2x2+x3=42x1+3x2+x4=12x1−x2+x5=3x1,x2,x3,x4,x5≥0,写出系数矩阵A=(p1,p2,p3,p4,p5)=−12123−1100010001。
第一次迭代:步骤一:选择初始可行基B=(p3,p4,p5)=100010001,则xB=(x3,x4,x5)T=B−1b=4123,故基本可行解为x(0)=(0,0,4,12,3)T,相应的目标函数值z(0)=cBxB=0。步骤二:计算非基变量的检验数σk=ck−cBB−1pk,对于非基变量x1,其检验数σ1=c1−cBB−1p1=4−[000]100010001−121=4,同理,x2的检验数为σ2=c2−cBB−1p2=1−[000]10001000123−1=1。步骤三:由于σ1,σ2都大于0,选取σi=max{σ1,σ2}=σ1,其对应的非基变量x1为进基变量;为选取离基变量,先求出Y1=B−1p1=100010001−121=−121,利用最小比值准则,找出min{ylibl∣yli>0}=min{y21b2,y31b3}=min{212,13}=13=y31b3,其中b=B−1b,,故xB中的第三个分量x5为离基变量,相应的新的基矩阵为B=(p3.p4,p1)。
第二次迭代:步骤一:新可行基的逆B−1=1000101−21,则基变量xB=(x3,x4,x1)T=B−1b=1000101−214123763,非基变量为xN=(x2,x5)T,于是初始基本可行解为x(1)=(3,0,7,6,0)T,相应的目标函数值z(1)=z(0)+(c1−z1)x1=0+4⋅3=12。步骤二:计算非基变量的检验数,公式为σk=ck−cBB−1pk,算出σ2=c2−cBB−1p2=1−[004]1000101−2123−1=5,同理σ5=0−1000101−21001=−4。步骤三:σi=max{σ2,σ5}=σ2,故对应的x2为进基变量;此时Y2=B−1p2=1000101−2123−1=15−1,找出min{ylibl∣yli>0}=min{17,56}=y22b2,其中b=B−1b,故xB中的第二个分量x4为离基变量,相应的新的基矩阵为B=(p3,p2,p1)=10023−1−121。
第三次迭代:新可行基的逆为B−1=100−51515157−5253,则xB=(x3,x2,x1)T=B−1b=100−51515157−52534123=52956521,xN=(x4,x5)T,故初始基本可行解为x(2)=(521,56,529,0,0)T,目标函数值为z(2)=z(1)+(c2−z2)x2=18。计算非基变量的检验数σ4=c4−cBB−1p4=0−[014]100−51515157−5253010=−1,同理σ5=−2,由前面定理知,所有非基变量的检验数σk≤0,则当前解为最优解。
综上,最优解x∗=x(2)=(521,56,529,0,0)T,最优值z∗=z(2)=18。
单纯形表
考虑线性规划问题的标准模型maxz=cx,s.t.{Ax=bx≥0,其中A=(aij)m×n,r(A)=m<n,x∈Rn,b∈Rm,cT∈Rn,假设B为初始的可行基,xB为基变量,xN为非基变量,则A=(B,N),x=[xBxN],c=(cB,cN),于是上述线性规划问题可以写成如下等价形式,maxz,s.t.⎩⎨⎧BxB+NxN=b−z+cBxB+cNxN=0xB≥0,xN≥0,进一步化为maxz,s.t.⎩⎨⎧xB+B−1NxN=B−1b−z+0⋅xB+(cN−cBB−1N)xN=−cBB−1bxB≥0,xN≥0(其中第二条用到了xB=B−1b−B−1NxN),该约束方程可以看作是以xB,xN为变量的方程组,把其系数列成表格,得到单纯形表:
z |
xB |
xN |
右端项 |
0 |
Em(单位阵) |
B−1N |
B−1b |
-1 |
0 |
cN−cBB−1N |
−cBB−1b |
若假设x1,x2,⋯,xm为基变量,可设计如下单纯形表:
|
cj |
|
c1 |
c2 |
⋯ |
cm |
cm+1 |
⋯ |
cn |
|
cB |
xB |
b |
x1 |
x2 |
⋯ |
xm |
xm+1 |
⋯ |
xn |
θ |
c1 |
x1 |
b1 |
1 |
0 |
⋯ |
0 |
a1,m+1 |
⋯ |
a1,n |
|
c2 |
x2 |
b2 |
0 |
1 |
⋯ |
0 |
a2,m+1 |
⋯ |
a2,n |
|
⋮ |
⋮ |
⋮ |
|
|
|
|
|
|
|
|
cm |
xm |
bm |
0 |
0 |
⋯ |
1 |
am,m+1 |
⋯ |
am,n |
|
|
−z |
−z0 |
0 |
0 |
⋯ |
0 |
σm+1 |
⋯ |
σn |
|
原表格
若按行看,c1,⋯,cm,cm+1,⋯,cn为价值系数,x1,⋯,xm,xm+1,⋯,xn为决策变量,中间的m×n矩阵包含着一个m阶方阵,该方阵对应着m个基变量,而0,⋯,0,σm+1,⋯,σn是检验行(其中基变量的检验数都是零);若按列看,cB,c1,⋯,cm是基变量的价值系数,xB,x1,⋯,xm是基变量,b,b1,⋯,bm是基变量取值(因为b=B−1b=xB),而最右边的θ这一列是最小比值列。
根据单纯形表,可以:
(1)检验当前基本可行解是否为最优解:若检验行所有的σj≤0,则已获得最优解,停止计算,否则要进行下一次步计算。
(2)检验是否为无界解:在σj>0(j∈JN)中,若有一个σm+t>0,而在单纯形表中σm+t所在列的Ym+t=B−1pm+t所有元素都≤0,则该问题无最优解,停止计算,否则进行下一步计算。
(3)选择进基变量:由进基变量选择准则maxj∈JNσj=σm+t(σj>0)(即在所有检验数为正的当中选取最大的检验数),选择进基变量xm+t,相应的列pm+t为进基向量,称表中pm+t所在列为主列。
(4)选择离基变量:由离基变量的最小比值准则θ=min{ai,m+tbi∣ai,m+t>0}=ar,m+tbr,则称第r行为主行,与主行对应的基变量xBr为离基变量。
(5)基变换:将可行基由(p1,p2,⋯,pr,⋯,pm)→(p1,⋯,pl−1,pm+t,⋯,pm),且将主列pm+t化为单位列向量er,即pm+t=a1,m+1a2,m+t⋮am,m+t⇒pr=0⋮1r⋮0。
例:用单纯形表解下列问题,max5x1+2x2+0x3+0x4+0x5,s.t.⎩⎨⎧30x1+20x2+x3=1605x1+x2+x4=15x1+x5=4x1,x2,x3,x4,x5≥0.
解:系数矩阵为A=(p1,p2,p3,p4,p5)=30512010100010001。
第一次迭代:步骤一,选取初始可行基B=100010001,B−1b=160154,其中xB=(x3,x4,x5)T,xN=(x1,x2)T,故初始基本可行解为x(0)=(0,0,160,15,4)T,z(0)=0。步骤二,计算检验数σj=cj−cBB−1pj,σ1=c1−cBB−1p1=5−[000]1000100013051=5,σ2=c2−cBB−1p2=2−[000]1000100012010=2,建立初始单纯形表
线性规划问题-单纯形法-单纯形表-例题-01
步骤三,选取检验数最大的5所在列的x1作为进基变量,此时x1所在列称为主列,再用b列对应元素除以主列对应元素(根据最小比值准则知,主列为负的不用除),即可得到θ列(最小比值列),最小比值3所在的行对应x4,于是选取x4为离基变量,这时检验数5所在列与最小比值3所在行交叉处的元素5(用中括号括起来),称这个5为主元素。
更新单纯形表:我们要将主元素5变为1,再把主元素所在列其他的元素变为0,这样主列才能变成单位向量。于是把主元素所在行所有元素都乘以51(中间的大矩阵及b列),把主元素变为1,同时把主元素适当倍数加到除主元素所在行以外的其他行(中间的大矩阵及b列),就能把主元素所在列其他的元素变为0。按顺序写出新的基变量xB及其价值函数系数cB,再次计算−z、检验数和最小比值,于是得到了新的单纯形表如下
线性规划问题-单纯形法-单纯形表-例题-02
此时检验数仍然有大于零的数,所以当前解还不是最优解,再次重复上面步骤,得到新的单纯形表
线性规划问题-单纯形法-单纯形表-例题-03
此时所有检验数都≤0,即当前解x∗=x(2)=(2,5,0,0,2)T是最优解,z∗=20。
这种方法的好处是,每一次都令B是单位阵,在计算检验数的过程中就不用反复计算B−1,并且也不需要计算Yi,直接用b除以主列对应元素即可得到θ。
线性规划问题解的数目
定理:在最大化问题中,若当前基本可行解的非基变量的检验数满足最优准则σj≤0(j∈JN),且其中有一个非基变量的检验数σj=0,则该线性规划问题有无穷多个最优解。
证明:设x(1)为当前的基本可行解,σj≤0(j∈JN),满足最优准则,故此时x∗=x(1),z∗=z(x(1))。不妨设x1为非基变量,且检验数σ1=0,选x1为进基变量,则x1对应的数列Y1=B−1p1有以下两种情况:
(1)若Y1中至少有一个分量yi1>0,则按最小比值准则选择离基变量,可由当前解迭代到下一个基本可行解x(2),而对应的目标函数值为z(x(2))=z(x(1))+(c1−z1)x1=z(x(1))(因为其中的c1−z1就是x1的检验数σ1),因此x(2)也是最优解,由凸集的性质可知,x(1)和x(2)连线上的任一点x^都是方程的可行解,且z(x^)=z(x(1))=z(x(2)),因此有无穷多个最优解。
(2)若Y1中所有分量yi1≤0,则不能运用最小比值准则,但仍然可以构造其他不同的最优解x(2)。把xB写成xB=xB1xB2⋮xBm=b1b2⋮bm−y11y21⋮ym1x1,令x1(2)=θ>0,xBi(2)=bi−yi1θ,i=1,2,⋯,m,其余xj(2)=0,不难发现,x(2)必为一个可行解(但不是基本可行解),由于θ可取任意正值,则相应的目标函数值z(x(2))=z(x(1))+(c1−z1)θ=z(x(1)),故x(2)也是最优解,也就是说该问题有无穷多个最优解。
定理:对线性规划问题的某一个基本可行解x(1)=[x1(1),x2(1),⋯,xm(1),0,⋯,0]T=[b1(1),b2(1),⋯,bm(1),0,⋯,0]T,若第t个非基变量xm+t的检验数σm+t>0,而xm+t对应的系数列向量pm+t≤0(即单纯形表中xm+t所在列的各元素都小于等于零),则该线性规划问题没有最优解。
例:求下列问题的最优解,maxx1+6x2+4x3,s.t.⎩⎨⎧−x1+2x2+2x3≤134x1−4x2+x3≤20x1+2x2+x3≤17x1≥1,x2≥2,x3≥3.
解:把问题标准为maxx1′+6x2′+4x3′+0x4+0x5+0x6,s.t.⎩⎨⎧−x1′+2x2′+2x3′+x4=44x1′−4x2′+x3′+x5=21x1′+2x2′+x3′+x6=9x1′,x2′,x3′,x4,x5,x6≥0,列出初始单纯形表
线性规划问题-单纯形法-线性规划问题解的数目-01
这时当前解还不是最优解,经过计算,得到最优单纯形表为:
线性规划问题-单纯形法-线性规划问题解的数目-02
由于有一个非基变量的检验数为零(蓝色那一列),根据定理,该问题一定有无穷多个最优解。验证:选x3′为进基变量,x5为离基变量,得到第二个最优单纯形表:
线性规划问题-单纯形法-线性规划问题解的数目-03
根据凸集的性质,这两个最优解连线上的任意一点都是最优解,即该问题有无穷多个最优解。
7. 大M法与两阶段法
人工变量及其处理方法:前面的单纯形法告诉我们,选取单位阵作为初始可行基可以很大程度地简化计算。为了便于计算及寻找计算方法上的规律性,希望在化线性规划为标准形式时,约束矩阵Am×n中含有一个m阶单位阵作为初始可行基,这有以下几种情况:
- 若在化标准形式前,m个约束条件都是“≤型”,则在化标准形式时,每一个约束条件左边都加上一个松弛变量xn+i,该松弛变量所对应的系数列向量就是单位向量ei,且这m个松弛变量所对应的系数列向量恰好组成了一个m阶单位阵。
- 若在化标准形式前,约束条件中有≥的不等式,则在约束条件左端减去剩余变量化为标准形式后,再加上一个非负的新变量,称为人工变量,显然该人工变量的系数列向量也为单位向量。(比如:2x1−x2+3x3≥30⇒2x1−x2+3x3−x5+x6=30)
- 若在化标准形式前,约束条件中有等式方程,则直接在该条件左端添加人工变量。
比如,考虑线性规划maxz=∑j=1ncjxj,s.t.{∑j=1nxjpj=bxj≥0,j=1,2,⋯,n,在每个约束方程左边加上一个人工变量xn+i(i=1,2,⋯,m)有⎩⎨⎧a11x1+a12x2+⋯+a1nxn+xn+1=b1a21x1+a22x2+⋯+a2nxn+xn+2=b2⋯am1x1+am2x2+⋯+amnxn+xn+m=bm,这样就让其中含有了一个m阶单位矩阵,以xn+1,xn+2,⋯,xn+m为基变量,得初始基本可行解x(0)=(0,⋯,0,b1,⋯,bm)T。
加入人工变量后的数学模型与未加人工变量的数学模型一般不是等价的。因此,人工变量与松弛变量或剩余变量是不同的,松弛变量或剩余变量只是将不等式改写为等式,而改写前后,两个约束是等价的。虽然加入人工变量后的新问题与原问题不等价,但它们有如下关系:
- 新问题的最优解中,若人工变量都处在非基变量的位置,则原问题有最优解,且将新问题最优解中去掉人工变量部分即为原问题的最优解。
- 新问题的最优解中,包含非零人工变量,则原问题无可行解。
- 新问题最优解的基变量中包含人工变量,但该人工变量取值为0,这时可将某个非基变量引入基变量中来替换该人工变量,从而得到原问题的最优解。
当以新问题的x(0)作为初始基本可行解进行迭代时,要怎样才能将所有的人工变量从基变量中“赶”出去,通常有大M法与两阶段法两种。
大M法
对于最大化问题,当以下式作为约束方程组(每个方程加一个人工变量)⎩⎨⎧a11x1+a12x2+⋯+a1nxn+xn+1=b1a21x1+a22x2+⋯+a2nxn+xn+2=b2⋯am1x1+am2x2+⋯+amnxn+xn+m=bm,若将目标函数修改为maxz=∑j=1ncjxj−Mxn+1−Mxn+2−⋯−Mxn+m,其中M是一个很大的正数,因为是对目标函数实现最大化,因此人工变量必须从基变量中迅速换出去,否则目标函数不可能实现最大化。
注:如果是最小化问题要用大M法,则把最小化目标函数化为minz=∑j=1ncjxj+Mxn+1+Mxn+2+⋯+Mxn+m。
例:求解线性规划问题maxz=3x1−x2−x3,s.t.⎩⎨⎧x1−2x2+x3≤11−4x1+x2+2x3≥3−2x1+x3=1x1,x2,x3≥0.
解:将问题化为maxz=3x1−x2−x3−Mx6−Mx7,s.t.⎩⎨⎧x1−2x2+x3+x4=11−4x1+x2+2x3−x5+x6=3−2x1+x3+x7=1xj≥0,j=1,2,⋯,7,写出初始单纯形表
线性规划问题-大M法-例题-01
这里的M看作是一个很大的正数即可(因此选取x3所在列为主列),第一次迭代:
线性规划问题-大M法-例题-02
第二次迭代:
线性规划问题-大M法-例题-03
第三次迭代:
线性规划问题-大M法-例题-04
对于最大化问题,检验数全部≥0时取得最优解,于是当前解为新问题的最优解。此时基变量中不含人工变量,且人工变量取值为零,所以当前解去掉x6,x7的部分就是原问题的最优解,即原问题的最优解为x∗=(4,1,9,0,0)T,z∗=2。
两阶段法
当线性规划问题添加人工变量后,将问题拆成两个线性规划问题:
(1)第一阶段:求解第1个线性规划minw=∑i=1mxn+i,s.t.⎩⎨⎧a11x1+a12x2+⋯+a1nxn+xn+1=b1a21x1+a22x2+⋯+a2nxn+xn+2=b2⋯am1x1+am2x2+⋯+amnxn+xn+m=bmxj≥0,j=1,2,⋯,n;i=1,2,⋯,m,即第1个线性规划的目标函数是对所有人工变量之和求最小,分以下情况进行讨论:
- 若求得的最优解中,所有人工变量都处在非基变量的位置,即xn+i=⋯=xn+m=0及w∗=0,则从第一阶段的最优解中去掉人工变量,得原问题的一个基本可行解。以其为原问题的初始基本解,进入第二阶段求解原问题的最优解。
- 若求得的最优解中,至少有一个人工变量不为零值,则说明添加人工变量之前的原问题无可行解,无需再进行第二阶段。
(2)第二阶段:得出原问题的一个基本可行解后,通过单纯形法对原问题进行计算。
例:考虑线性规划问题maxz=3x1−x2−x3,s.t.⎩⎨⎧x1−2x2+x3≤11−4x1+x2+2x3≥3−2x1+x3=1x1,x2,x3≥0.
解:建立第一阶段的线性规划问题minz=x6+x7,s.t.⎩⎨⎧x1−2x2+x3+x4=11−4x1+x2+2x3−x5+x6=3−2x1+x3+x7=1xj≥0,j=1,2,⋯,7,令B(0)=(p4,p6,p7)作初始可行基,作初始单纯形表:
线性规划问题-两阶段法-例题-01
由于这里是最小化问题,当所有检验数都≥0时才取得最优解,选取第三列为主列,进行第一次迭代:
线性规划问题-两阶段法-例题-02
第二次迭代:
线性规划问题-两阶段法-例题-03
此时已达到最优值,因此第一阶段的最优解为x=(0,1,1,12,0,0,0)T,将该表中人工变量列划去,即可得到第二阶段的初始单纯形表(由于第二阶段目标函数不同,涉及目标函数的部分需重新计算),且x(0)=(0,1,1,12,0)T为第二阶段的初始基本可行解。
建立第二阶段的线性规划问题maxz=3x1−x2−x3,s.t.⎩⎨⎧x1−2x2+x3+x4=11−4x1+x2+2x3−x5=3−2x1+x3=1xj≥0,j=1,2,⋯,5,写出初始单纯形表:
线性规划问题-两阶段法-例题-04
第一次迭代:
线性规划问题-两阶段法-例题-05
因此原问题的最优解x∗=(4,1,9,0,0)T,z∗=2。
8. 对偶问题
对偶问题概述
引入:一个线性规划问题往往伴随着与之配对的,两者有密切联系的另一个线性规划问题,我们将其中一个称为原问题,另一个称为对偶问题。对偶问题在经济学上有重要意义,下面通过一个问题来说明。
问题:某家具长木器车间生产木门和木窗两种产品,加工木门收入为56元/扇,加工木窗收入为30元/扇,生产一扇木门需要木工4小时、油漆工2小时,生产一扇木窗需要木工3小时、油漆工1小时。该车间每日可用木工总工时为120小时,油漆工总工时为50小时,问该车间应如何安排生产才能使每日收入最大?
解:令该车间每日安排生产木门x1扇、木窗x2扇,则数学模型为maxz=56x1+30x2,s.t.⎩⎨⎧4x1+3x2≤1202x1+x2≤50x1,x2≥0,用图解法或单纯形法,可求得最优解x∗=(x1,x2)T=(15,20)T,z∗=1440元。
现在从另一个角度来考虑车间的生产问题。假设有一个个体经营者,手中有一批木器家去生产订单,他想利用该木器生产车间的木工与油漆工来完成他的订单,他就要实现考虑付给该车间每个工时的价格。他可以构造一个数学模型来研究如何定价才能:(1)木器生产车间觉得有利可图从而愿意替他加工这批订单(2)他自己所付的工时费用总数最小。设w1为付给木工每个工时的价格,w2为付给油漆工每个工时的价格,则该个体经营者的目标函数为每日所付工时总费用最小,即minf=120w1+50w2,但该个体经营者所付的价格不能太低,至少不能低于该车间生产木门、木窗时所得到的收入,否则车间觉得无利可图就不会替他加工这批订单,因此w1,w2应满足s.t.⎩⎨⎧4w1+2w2≥563w1+w2≥30w1,w2≥0,解出w∗=(2,24)T,f∗=1440。
对比两个问题,maxs.t.原问题z=56x1+30x24x1+3x2≤1202x1+x2≤50x1,x2≥0x∗=(15,20)Tz∗=1440与mins.t.对偶问题f=120w1+50w24w1+2w2≥563w1+w2≥30w1,w2≥0w∗=(2,24)Tf∗=1440,他们有某些共同点:原问题的价值系数是对偶问题约束条件的右端项,对偶问题价值系数是原问题约束条件右端项,原问题约束不等式系数的行是对偶问题约束不等式系数的列,原问题约束不等式系数的列是对偶问题约束不等式系数的行,…
原问题与对偶问题:原问题与对偶问题一般表示为maxs.t.LPz=cxAx≤bx≥0n×1,mins.t.DPf=wbwA≥cw≥01×m,其中c1×n,xn×1,w1×m,Am×n,bm×1,两个问题之间的对应关系总结如下表(如果原问题是最小化问题,则从右侧对应左侧):
线性规划问题-对偶问题-原问题与对偶问题的关系-01
比如:写出下列问题的对偶问题:mins.t.原问题z=25x1+2x2+3x3−x1+x2−x3≤1x1+2x2−x3≥12x1−x2+x3=1x1≥0,x2≤0,x3无限制,maxs.t.对偶问题w1+w2+w3−w1+w2+2w3≤25w1+2w2−w3≥2−w1−w2+w3=3w1≤0,w2≥0,w3无限制.
对偶理论
对于maxs.t.LPz=cxAx≤bx≥0n×1与mins.t.DPf=wbwA≥cw≥01×m,其中c1×n,xn×1,w1×m,Am×n,bm×1,有如下定理:
弱对偶性定理:设x(0)及w(0)分别是LP和DP的任一可行解,则恒有cx(0)≤w(0)b。
证明:cx(0)≤w(0)Ax(0)≤w(0)b.
该定理说明了,最大化问题的任一可行解的目标函数值都是其对偶最小化问题目标函数的下界。最小化问题的任一可行解的目标函数值都是其对偶最大化问题目标函数值的上界。
定理:设x(0)及w(0)分别是LP及DP问题的可行解,则当cx(0)=w(0)b时,x(0)与w(0)必分别是各自问题的最优解。
证明:设x是LP问题的任意一个可行解,则由弱对偶性定理知,必有cx≤w(0)b=cx(0),说明x(0)是LP问题的最优解;同理可证w(0)是DP问题的最优解。
定理:若原问题LP与对偶问题DP同时有可行解,则它们必都有最优解。
定理:若原问题的目标函数无界,则其对偶问题必无可行解。
强对偶定理:设LP与DP中有一个最优解,则另一个问题也必存在最优解,且两个问题最优解的目标函数值必相等。
证明:将原问题标准化为maxs.t.z′=cx+caxaAx+Exa=bx≥0,xa≥0,其中xa为松弛变量,ca为松弛变量的系数(显然ca=0),E为单位矩阵,设其相应的最优基为B,则所有变量的检验数σj≤0。考虑变量x1,x2,⋯,xn的检验数(c1,c2,⋯,cn)−cBB−1(p1,p2,⋯,pn)≤0,可以写成c−cBB−1A≤0,若记w(0)=cBB−1,代入上式得w(0)A≥c,因此w(0)满足对偶问题的约束条件;再考虑xn+1,xn+2,⋯,xn+m的检验数(0,0,⋯,0)−cBB−1(pn+1,pn+2,⋯,pn+m)≤0,可以写成0−cBB−1E≤0,那么w(0)=cBB−1≥0,因此w(0)=cBB−1是对偶问题的可行解,此时原问题与对偶问题的目标函数值有w(0)b=(cBB−1)b=cB(B−1b),故w(0)=cBB−1为对偶问题的最优解(前面的定理,目标函数值相等即为各自最优解)。
对称形式的互补松弛性定理:设x(0)与w(0)分别是对称形式的原问题及其对偶问题的两个可行解,则x(0)与w(0)分别是各自问题的最优解的充分必要条件为:对所有的i,j,下列各式都成立:
(1)若xj(0)>0,必有w(0)pj=cj;
(2)若w(0)pj>cj,必有xj(0)=0;
(3)若wi(0)>0,必有Aix(0)=bi;
(4)若Aix(0)<bi,必有wi(0)=0。
其中Ai为原问题约束矩阵A的第i行,pj为A的第j列,w(0)=(w1(0),w2(0),⋯,wm(0))是对偶变量,x(0)=(x1(0),x2(0),⋯,xn(0))T是原问题的决策变量,cj是原问题的价值系数,bi是原问题约束方程右端项。
例:若已知maxx1+2x2,s.t.⎩⎨⎧3x1+x2≤2−x1+2x2≤3x1−3x2≤1x1,x2≥0的最优解为x∗=(71,711)T,求其对偶问题的最优解
解:对偶问题为min2w1+3w2+w3,s.t.⎩⎨⎧3w1−w2+w3≥1w1+2w2−3w3≥2w1,w2,w3≥0,根据互补松弛定理的第(1)条,因为所有的x1∗,x2∗>0,列出w∗的方程得{3w1∗−w2∗+w3∗=1w1∗+2w2∗−3w3∗=2,此时方程还无法求解,于是再根据互补松弛定理的第(4)条,再将x∗=(71,711)T代入原问题的约束条件中,得⎩⎨⎧3x1∗+x2∗=2−x1∗+2x2∗=3x1∗−3x2∗=−732<1,所以w3∗=0,最后解得对偶问题最优解为w∗=(74,75,0)T.
对偶单纯形法
引入:单纯形法从一个基本可行解迭代到下一个基本可行解时,总是保持解的可行性不变,变化的只是检验数向量σ,对于最大化问题,当σ≤0时最优解诞生,此时由σ=c−cBB−1A=c−w(0)A≤0及w(0)=cBB−1≥0知,w(0)=cBB−1是对偶问题的可行解。因此,我们可以这样来描述单纯形法的过程:从基本解向最优解迭代时,是始终保持原问题可行性的条件下,其对偶问题由不可行(c≰w(0)A)向可行(c≤w(0)A)转化,一旦其对偶问题也成为可行解时,根据前面定理,原问题的可行解即成为最优解。而对偶单纯形法是把上面的过程反过来,在迭代的过程中,始终保持对偶问题解的可行性,而使原问题的解由不可行逐渐向可行转化,一旦原问题的解满足了可行性,对偶问题也就达到了最优解。
正则解与正则基的定义:设x(0)是原问题的一个基本解,对应的基是B,若它所对应的检验数向量σ=c−cBB−1A≤0成立,则称x(0)是原问题的一个正则解,对应的基矩阵B称为正则基。
正则解对原问题而言只是一个基本解,并不要求是可行解,而其σ≤0成立,即满足了对偶问题的可行性,称之为正则。因此用正则性语言来描述对偶单纯形法的思路是:在保持正则解的正则性不变的条件下,在迭代过程中,使原问题的不可行性逐步消失,一旦迭代到可行解时,即达到了最优解。
对偶单纯形法计算步骤:
(1)给定一个初始正则解x(0),对应的正则基为B;
(2)计算b=B−1b,若b=B−1b≥0时,则停止计算,当前的正则解x=B−1b便是最优解,否则继续进入到下一步;
(3)确定离基变量:令br=min{bi∣bi<0},则xBr为离基变量(此时称第r行称为主行);
(4)检查单纯形表中第r行的系数:Ar′=(ar1′,ar2′,⋯,arm′),若所有的arj′≥0,则原问题无可行解,否则转下一步(理由如下:对于第r个约束方程xBr=br′−∑j∈JNarj′xj,因为所有的arj′≥0,因此不论xj(j∈JN)改为怎样的正值,都无法使xBr的值转化为正数,故本问题无可行解);
(5)确定进基变量:取ark′σk=min{arj′σj∣arj′<0,j∈JN},则xk是进基变量,pk是主列;
(6)迭代:以ark′为主元素进行消元变换,迭代到下一个正则解。
对偶单纯形法与单纯形法的不同之处在于,先根据b确定离基变量得到主行,再用检验数除以主行中的负元素得到比值确定进基变量。
例:用对偶单纯形法求解minz=2x1+3x2+4x3,s.t.⎩⎨⎧x1+2x2+x3≥32x1−x2+3x3≥4x1,x2,x3≥0.
解:将问题标准化为maxz=−2x1−3x2−4x3+0x4+0x5,s.t.⎩⎨⎧x1+2x2+x3−x4=32x1−x2+3x3−x5=4x1,x2,x3,x4,x5≥0,为了利用p4,p5作为初始基,将两个约束方程左右两端各乘−1,得maxz=−2x1−3x2−4x3+0x4+0x5,s.t.⎩⎨⎧−x1−2x2−x3+x4=−3−2x1+x2−3x3+x5=−4x1,x2,x3,x4,x5≥0,写出初始单纯形表:
线性规划问题-对偶问题-对偶单纯形法-例题-01
第一次迭代:
线性规划问题-对偶问题-对偶单纯形法-例题-02
第二次迭代:
线性规划问题-对偶问题-对偶单纯形法-例题-03
此时达到了最优解。
9. 单纯形法的灵敏度分析
引入:之前的计算中,总是把c(价值系数)、A(约束矩阵)、b(约束方程右端项)视作常数,而灵敏度分析是讨论他们如果变化,会对解产生什么影响。
灵敏度分析通常有两类问题:
- 当c、A、b中某一部分数据发生变化时,最优解与最优值怎么变?
- 研究c、A、b中数据在多大范围内波动时,使原有最优解仍为最优解,同时讨论此时最优值如何变动?
设考虑的问题为maxz=cx,s.t.{Ax=bx≥0,相应的最优单纯形表为:
线性规划问题-灵敏度分析-01
下面分情况讨论(注意计算时由于题目数据发生了改变,大部分情况下都要用题目的数据进行重算,而不能用单纯形表中的数据计算):
价值系数c发生改变
当c由c→c+Δc时,最优表会发生改变的只是最后一行(检验行)。设c=(c1,⋯,cn),则c+Δc=(c1+Δc1,⋯,cn+Δcn),则检验数应修改为(c+Δc)−(cB+ΔcB)B−1A,而目标函数值应修改为(cB+ΔcB)B−1b。此时:
- 若检验数仍保持≤0,则原最优解仍为最优解,但目标函数值已变;
- 若检验数不满足最优性条件,则当前解已不是最优解,要从修改后的单纯形表出发,重新进行迭代。
例1:已知标准形式的线性规划问题maxz=−x1+2x2+x3,s.t.⎩⎨⎧x1+x2+x3+x4=62x1−x2+x5=4x1,x2,⋯,x5≥0,其最优单纯形表如下
线性规划问题-灵敏度分析-例题1-01
问:(1)当c1由−1变为c1+Δc1=4时,求新问题的最优解(2)讨论c2在什么范围内变化时,原有的最优解仍是最优解
解:(1)因为只有非基变量的x1系数c1发生了改变,只需重新计算σ1′=(c1+Δc1)−cBB−1p1=σ1+5=(−3)+5=2>0,可见最优性准则已不满足,修改上表后重新迭代可得问题的最优解:
线性规划问题-灵敏度分析-例题1-02
(2)要使原最优解仍为最优解,须在新条件下仍满足σ≤0,记c2′=c2+Δc2,由于x2为基变量,cB发生了改变全部检验数都要重算,于是σ=c−cBB−1A=(c1,c2′,c3,c4,c5)−cB′B−1(p1,p2,p3,p4,p5)=(−1,c2′,1,0,0)−(c2′,c5)[1−101]−1[121−1101001]=(−1,c2′,1,0,0)−(c2′,0)[1101][121−1101001]=(−1,c2′,1,0,0)−(c2′,c2′,c2′,c2′,0)=(−1−c2′,0,1−c2′,−c2′,0)≤0,即当c2′≥1或写成Δc2≥−1时,原最优解仍为最优解。
右端项b发生改变
当右端列向量b→b+Δb,改变的只是表中第三列(右端列),即基变量的取值由xB=B−1b→xB′=B−1(b+Δb),而目标函数值由−z=−cBB−1→−z′=−cBB−1(b+Δb)。此时:
- 若xB′=B−1(b+Δb)≥0仍成立,则因为σj(j∈JN)没有改变,此时xB′为新问题的最优解,z′为新问题的最优值;
- 若xB′=B−1(b+Δb)≱0,但因为σj≤0仍成立,故[xB′xN′]=[B−1(b+Δb)0]是一个正则解,故可用对偶单纯形法再次进行迭代,直到求得新的最优解。
例2:已知线性规划问题maxz=−x1−x2+4x3,s.t.⎩⎨⎧x1+x2+2x3+x4=9x1+x2−x3+x5=2−x1+x2+x3+x6=4x1,x2,⋯,x6≥0,若右端列向量从b=(9,2,4)T→(3,2,3)T,求新问题的最优解。
解:原问题的最优单纯形表如下:
当b由(9,2,4)T改变为(3,2,3)T时,xB′=B−1(b+Δb)=11−10102−11−1323=31031010−32131323=−152≱0,−z′=−cBB−1(b+Δb)=−[−104]−152=−9,当前解只是一个正则解,还要用对偶单纯形法进行迭代:
线性规划问题-灵敏度分析-例题2-02
第一次迭代:
线性规划问题-灵敏度分析-例题2-03
故新问题的最优解是x∗=(0,0,23,0,27,23)T,最优值为z∗=6。
约束系数列向量pk发生改变
设系数列向量由pk→pk+Δpk=(a1k+Δa1k,a2k+Δa2k,⋯,amk+Δamk)。此时:
- 若pk不属于最优基B,即xk不是基变量,这时非基变量xk的检验数为σk′=ck−cBB−1(pk+Δpk),那么
- σk′≤0,则原最优解仍为最优解,最优值也不变;
- σk′>0,则最优性准则已不满足,修改最优单纯形表中的第k列,即B−1pk→B−1(pk+Δpk),以xk作进基变量进行迭代,求出新问题的最优解;
- 若pk属于最优基B,即xk为基变量,此时表中所有的系数都要发生改变,因此还是用单纯形法重新计算比较方便。
例3:已知线性规划问题maxz=2x1+3x2+x3,s.t.⎩⎨⎧31x1+31x2+31x3+x4=131x1+34x2+37x3+x5=3x1,x2,⋯,x5≥0,若p3由原来的(31,37)T→(101,31)T,最优解将如何变化?
解:原问题的最优单纯形表如下
线性规划问题-灵敏度分析-例题3-01
若p3发生改变,此时σ3=c3−cBB−1p3′=1−[23][31313134]−1[10131]=1−[23][4−1−11][10131]=61>0,最优准则已不满足,修改最优单纯形表中的第三列为[p3′σ3]=[B−1p3σ3]=15130761,选出主元素再进行迭代
线性规划问题-灵敏度分析-例题3-02
第一次迭代:
线性规划问题-灵敏度分析-例题3-03
故新问题的最优解是x∗=(73,0,760,0,0)T。
增加一个新变量xn+1
考虑线性规划问题maxz=cx,s.t.{Ax=bx≥0,设该问题已经得到了最优解x∗=(x1∗,x2∗,⋯,xn∗)T,最优基为B,现追加一个新变量xn+1,其价值系数为cn+1,系数列向量为pn+1=(a1,n+1,a2,n+1,⋯,am,n+1)T,得到新问题max=f=c1x1+c2x2+⋯+cnxn+cn+1xn+1,s.t.{ai1x1+⋯+ainxn+ai,n+1xn+1=bi,i=1,2,⋯,mxj≥0,j=1,2,⋯,n,n+1,这时原问题的最优基B是新问题的可行基,原由变量的检验数都保持不变,而σn+1=cn+1−cBB−1pn+1。此时:
- 若σn=1≤0,则新问题的最优性准则仍满足,故x∗=(x1∗,x2∗,⋯,xn∗,0)T是新问题的最优解,此时xn+1=0,说明所追加的新变量xn+1对最优解没有影响,或说新增加的内容对总的结果是不利的(因为若xn+1>0则达不到最优解);
- 若σn+1>0,说明新增加的内容对总的结果有利,故x∗=(x1∗,x2∗,⋯,xn∗,0)T一定不是新问题的最优解,这时在原问题的最优单纯形表上增加一列[pn+1′σn+1]=[B−1pn+1cn+1−cBB−1pn+1],并以xn+1作为进基变量继续迭代。
例4:已知线性规划问题maxz=−x1−x2+4x3,s.t.⎩⎨⎧x1+x2+2x3+x4=9x1+x2−x3+x5=2−x1+x2+x3+x6=4x1,x2,⋯,x6≥0,现增加新变量x7,且c7=3,p7=(3,1,−3)T,求新问题的最优解。
解:原问题的最优单纯形表如下
线性规划问题-灵敏度分析-例题4-01
计算决策变量x7的检验数σ7=c7−cBB−1p7=3−[−104]11−10102−11−131−3=3−[−104]31031010−3213131−3=6≥0,说明新增加的x7对总的结果有利,目标函数肯定有所增加。在原最优单纯形表加入列[p7′σ7]=[B−1p7σ7]=(3,−2,0,6)T,然后以p7作为进基变量进行迭代
线性规划问题-灵敏度分析-例题4-02
第一次迭代:
线性规划问题-灵敏度分析-例题4-03
故新问题的最优解是x∗=(0,0,313,0,956,0,91)T。
10. 运输问题
运输问题的数学模型
运输问题是一种典型的线性规划问题,所以它可以用单纯形法求解,但是其数学模型有特殊的结构,存在一种比单纯形法更简便的计算方法–表上作业法(表上作业法本质上仍是单纯形法)。下面通过例子来了解运输问题:
问题:设有m个生产点A1,A2,⋯,Am可供应某种物质,其生产量分别为a1,a,1,⋯,am,另有n销售点B1,B2,⋯,Bn,其销售量分别为b1,b2,⋯,bn。如果从Ai到Bj运输单位物质的运价为cij,那么在产销平衡的条件下(∑i=1mai=∑j=1nbj),如何设计调运方案才使得运费最小?
解:将问题通过表格来描述,其中xij表示从Ai到Bj的发运量:
|
B1 |
B2 |
B3 |
⋯ |
Bn |
A1 |
x11[c11] |
x12[c12] |
x13[c13] |
⋯ |
x1n[c1n] |
A2 |
x21[c21] |
x22[c22] |
x23[c23] |
⋯ |
x2n[c2n] |
⋯ |
⋯ |
⋯ |
⋯ |
⋯ |
⋯ |
Am |
xm1[cm1] |
xm2[cm2] |
xm3[cm3] |
⋯ |
xmn[cmn] |
于是可以得到数学模型:min∑i=1m∑j=1ncijxij,s.t.⎩⎨⎧∑j=1nxij=ai∑i=1mxij=bjxij≥0i=1,2,⋯,mj=1,2,⋯,ni=1,2,⋯,m和j=1,2,⋯,n,把前两个约束展开,得到⎩⎨⎧x11+x12+⋯+x1n=a1x21+x22+⋯+x2n=a2⋯xm1+xm2+⋯+xmn=amx11+x21+⋯+xm1=b1x12+x22+⋯+xm2=b2⋯x1n+x2n+⋯+xmn=bn,其系数矩阵为A=1111⋯⋱111111⋯⋱11⋱⋯1111⋯⋱11(m+n)×(mn),观察可知:
A是一个结构特殊的稀疏矩阵,每列的(m+n)个元素中只有两个元素为1,其余的为0;
变量xij所对应的系数列向量pij=ei+em+j=(0,⋯,1,⋯,0,⋯,1,⋯,0)T,即只有第i个和第m+j位置元素为1,其余为0;
矩阵A的秩为(m+n−1);
证明:已知A是一个(m+n)×mn型矩阵,并假设mn>m+n,当产销平衡时前m个方程之和与后n个方程之和相等,即产销平衡的约束方程组中有一个方程是多余的,亦即矩阵A的前m行之和与后n行之和相等,故矩阵A的任一(m+n)阶子式比为零,即R(A)≤m+n−1;将A的第一行删去得A=[p11′,p12′,⋯,p1n′,p21′,p22′,⋯,p2n′,⋯,pm1′,pm2′,⋯,pmn′],取A的前n列,第n+1列、第2n+1列,⋯和第(m−1)n+1列构成一个(m+n−1)阶子矩阵A=0101⋱⋱01110110⋱⋯⋯110=[p11′,⋯,p1n′,p21′,p31′,⋯,pm1′],把它看作是2×2的分块矩阵,容易算出∣A∣=±1=0,也就是说R(A)=m+n−1,即向量组p11′,⋯,p1n′,p21′,p31′,⋯,pm1′线性无关,即向量组p11,⋯,p1n,p21,p31,⋯,pm1也线性无关,故R(A)=m+n−1。
定理(运输问题解的存在性):产销平衡问题存在最优基本可行解。
证明:记∑i=1mai=∑j=1nbj=k,若构造xij=kaibj(i=1,2,⋯,m,j=1,2,⋯,n),将xij代入约束方程得{∑j=1nxij=∑j=1nkaibj=kai∑j=1nbj=ai,i=1,2⋯,m∑i=1mxij=∑i=1mkaibj=kbj∑i=1mai=bj,j=1,2,⋯,n,因为ai≥0,bj≥0,所以xij≥0是产销平衡问题的可行解,故原问题有基本可行解。
再由xij的定义知,xij≤min{ai,bj},即产销问题的任意可行解是有限的,故其可行域是有界的,因此原问题必有最优解。
(补充)产销不平衡问题的数学模型:当产量大于销量时,有∑i=1mai>∑j=1nbj,相应的数学模型为min∑i=1m∑j=1ncijxij,s.t.⎩⎨⎧∑j=1nxij≤ai∑i=1mxij=bjxij≥0i=1,2,⋯,mj=1,2,⋯,ni=1,2,⋯,m和j=1,2,⋯,n,
表上作业法:要求解产销平衡时的运输问题,可以用西北角法或最小元素法确定初始可行解,然后用位势法求检验数,再用闭合回路法调整基本可行解,由于以上都在表(价格表或调运表)中进行,所以它们都是表上作业法。
首先了解几个相关的定义:将变量xij在调运表中对应的空格记作(i,j),称为格点(i,j),将xij的系数列向量pij称为格点(i,j)所对应的系数列向量。若xij为基变量,则(i,j)称为基格,否则称为非基格。
西北角法
下面通过例题了解西北角法如何确定初始可行解:
例:某公司生产糖果,它有3个加工厂A1,A2,A3,每月产量分别为7t,4t,9t,该公司把产品分别运往4个销售店B1,B2,B3,B4,每月销量分别为3t,6t,5t,6t。已知从第i个加工厂到第j个销售店的每吨糖果的运价cij见下表,试设计在满足销售店需求量的前提下,个加工厂到每个销售店的调运方案,使该公司所花总的运费最小。
价格表:
|
B1 |
B2 |
B3 |
B4 |
A1 |
3 |
11 |
3 |
10 |
A2 |
1 |
9 |
2 |
8 |
A3 |
7 |
4 |
10 |
5 |
解:画出调运表如下:
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
|
|
|
|
7 |
A2 |
|
|
|
|
4 |
A3 |
|
|
|
|
9 |
收量 |
3 |
6 |
5 |
6 |
|
要填写该表,可以用西北角法,观察最西北角处的格子(A1行与B1列交界处),此时发量(即产量)为7,而收量(即需求量)为3,因此可以满足,B1不再需要来自A2、A3的货物,画斜线,于是填写如下
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
3 |
|
|
|
7-3 |
A2 |
/ |
|
|
|
4 |
A3 |
/ |
|
|
|
9 |
收量 |
3-3 |
6 |
5 |
6 |
|
那么剩下的格子的西北角是A1与B2交界处位置,此时收量为6,而发量只剩下4,于是剩下的B3、B4不可能再接收来自A1的货物,划斜线,而剩下的收量2需要来自A2,于是填写如下
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
3 |
4 |
/ |
/ |
7-3-4 |
A2 |
/ |
2 |
|
|
4-2 |
A3 |
/ |
/ |
|
|
9 |
收量 |
3-3 |
6-4-2 |
5 |
6 |
|
以此类推,最后填写如下
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
3 |
4 |
/ |
/ |
7-3-4 |
A2 |
/ |
2 |
2 |
/ |
4-2-2 |
A3 |
/ |
/ |
3 |
6 |
9-3-6 |
收量 |
3-3 |
6-4-2 |
5 |
6 |
|
故原问题的一个初始可行解x(0)=(3,4,0,0,0,2,2,0,0,0,3,6)T,相应的总费用为1330元。
最小元素法
下面通过例题了解最小元素法如何确定初始可行解:
例题同西北角法,但是这次使用最小元素法确定初始可行解,价格表如下:
|
B1 |
B2 |
B3 |
B4 |
A1 |
3 |
11 |
3 |
10 |
A2 |
1 |
9 |
2 |
8 |
A3 |
7 |
4 |
10 |
5 |
解:画出调运表如下
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
|
|
|
|
7 |
A2 |
|
|
|
|
4 |
A3 |
|
|
|
|
9 |
收量 |
3 |
6 |
5 |
6 |
|
要填写该表,根据最小元素法,从价格表中找出运费最低的是1,即从A2运往B1,于是在调运表对应位置尽可能发更多的货,对应位置发货量为3可以满足,于是此时B1达到饱和,不需要再来自A1、A3的货物,画斜线
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
/ |
|
|
|
7 |
A2 |
3 |
|
|
|
4-3 |
A3 |
/ |
|
|
|
9 |
收量 |
3-3 |
6 |
5 |
6 |
|
再找下一个最小的元素,是从A2运往B3的2,A2给B3发完1的货物后,已经不可能再给B2、B4发货了,于是填写如下
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
/ |
|
|
|
7 |
A2 |
3 |
/ |
1 |
/ |
4-3-1 |
A3 |
/ |
|
|
|
9 |
收量 |
3-3 |
6 |
5-1 |
6 |
|
再找下一个最小的元素,是从A1运往B3的3,以此类推,最后填写如下
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
/ |
/ |
4 |
3 |
7-4-3 |
A2 |
3 |
/ |
1 |
/ |
4-3-1 |
A3 |
/ |
6 |
/ |
3 |
9-6-3 |
收量 |
3-3 |
6-6 |
5-1-4 |
6 |
|
故原问题的一个初始可行解x(0)=(0,0,4,3,3,0,1,0,0,6,0,3)T。
位势法
位势法的原理是利用基变量检验数为零,得到公式cij=(ui+vj)((i,j)∈JB,ui,vj是其对偶问题的决策变量),依此得到方程组解出ui,vj,再利用非基变量检验数公式σij=cij−(ui+vj)算出非基变量的检验数。
推导:在原来检验数公式的基础上修改,得到双下标形式的检验数公式σij=cij−zij=cij−cBB−1pij,若记cBB−1=w,则w就是对偶问题的决策变量,设w=(u1,⋯,um,v1,⋯,vn),则运输问题的对偶问题为maxf=∑i=1maiui+∑j=1nbjvj,s.t.{ui+vj≤cij,i=1,2⋯,m;j=1,2,⋯,nui,vj无正负限制,决策变量的检验数为σij=cij−wpij=cij−(u1,⋯,um,v1,⋯,vn)(ei+em+j)=cij−(ui+vj),(i=1,2,⋯,m;j=1,2,⋯,n),因为基变量的检验数σij=0,故当(i,j)∈JB时cij−(ui+vj)=0,这是一个含有(m+n−1)个方程、(m+n)个未知量的线性方程组,为求解方程组,我们人为规定令vn=0,这样就可以求出所有的ui,vj,以及所有非基变量的检验数。
下面通过例题了解位势法如何求检验数(注意位势法是在价格表上进行的):
同样是前面的问题,用西北角法求出的可行解x(0)=(3,4,0,0,0,2,2,0,0,0,3,6)T为例,下面是其对应的调运表:
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
3 |
4 |
/ |
/ |
7-3-4 |
A2 |
/ |
2 |
2 |
/ |
4-2-2 |
A3 |
/ |
/ |
3 |
6 |
9-3-6 |
收量 |
3-3 |
6-4-2 |
5 |
6 |
|
在例题的价格表上加上一行一列,其中ui、vj是其对偶问题的决策变量,在对应的调运表中非零元素的位置作加粗标识(其对应的xij是基变量),非基变量用中括号括住,如下表
|
B1 |
B2 |
B3 |
B4 |
ui |
A1 |
3 |
11 |
[3] |
[10] |
|
A2 |
[1] |
9 |
2 |
[8] |
|
A3 |
[7] |
[4] |
10 |
5 |
|
vj |
|
|
|
|
|
由于检验数公式为cij−(ui+vj),所以我们应将ui和vj确定下来,首先根据前面的分析,vn人为规定为0,于是填写如下
|
B1 |
B2 |
B3 |
B4 |
ui |
A1 |
3 |
11 |
[3] |
[10] |
|
A2 |
[1] |
9 |
2 |
[8] |
|
A3 |
[7] |
[4] |
10 |
5 |
|
vj |
|
|
|
0 |
|
以某个基格为标准,结合公式cij=(ui+vj)(仅当(i,j)∈JB时成立)计算运输问题的对偶问题的基变量ui和vj,比如c34=u3+v4⇒5=u3+0⇒u3=5,(注意只看基变量对应的格子,即加粗元素的位置),最后填写如下
|
B1 |
B2 |
B3 |
B4 |
ui |
A1 |
3 |
11 |
[3] |
[10] |
-1 |
A2 |
[1] |
9 |
2 |
[8] |
-3 |
A3 |
[7] |
[4] |
10 |
5 |
5 |
vj |
4 |
12 |
5 |
0 |
|
对于非基变量的检验数,通过σij=cij−(ui+vj)计算,比如σ31=c31−(u3+v1)=7−(5+4)=−2,把检验数写在价格表上,如下表(原非基变量的cij用中括号括住,而基变量的检验数一定为0就不填写了)
|
B1 |
B2 |
B3 |
B4 |
ui |
A1 |
3 |
11 |
[3]-1 |
[10]11 |
-1 |
A2 |
[1]0 |
9 |
2 |
[8]11 |
-3 |
A3 |
[7]-2 |
[4]-13 |
10 |
5 |
5 |
vj |
4 |
12 |
5 |
0 |
|
由于问题是最小化问题,当所有检验数都≥0时得到最优解,显然当前解还不是最优解,还要再用闭合回路法对基本可行解进行调整。
闭合回路法
定义:若一组格点经过适当的排序后,能写成以下形式:(i1,j1)、(i1,j2)、(i2,j2)、(i2,j3)、(i3,j3)、⋯、(is,js)、(is,j1),则称这组格点构成了闭合回路。
比如以下格点:
线性规划问题-运输问题-闭合回路法-01
由梨花组成的(1,1),(1,2),(3,2),(3,1)就构成了闭合回路;再比如心形组成的(3,3),(1,3),(1,4),(2,4),(2,5),(3,5)也构成了闭合回路,…
下面通过例题了解闭合回路法如何调整基本可行解:
前面通过位势法得到了下面价格表:
|
B1 |
B2 |
B3 |
B4 |
ui |
A1 |
3 |
11 |
[3]-1 |
[10]11 |
-1 |
A2 |
[1]0 |
9 |
2 |
[8]11 |
-3 |
A3 |
[7]-2 |
[4]-13 |
10 |
5 |
5 |
vj |
4 |
12 |
5 |
0 |
|
同单纯形法选择进基变量的准则,对于最小化问题,找出最小检验数对应变量作为进基变量,σkl=mini,j∈JN{σij∣σij<0}=σ32=−13,即x32为进基变量,也就是说(3,2)为进基格,为确定离基变量,以进基格为起点作一个闭合回路,要求除起始格为进基格外,该闭合回路的其余顶点均为基格。
我们在调运表中进行操作(进基格用∙标注):
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
3 |
4 |
|
|
7 |
A2 |
|
2 |
2 |
|
4 |
A3 |
|
∙ |
3 |
6 |
9 |
收量 |
3 |
6 |
5 |
6 |
|
找出任意一个闭合回路,比如(3,2),(3,3),(2,3),(2,2),考虑(闭合回路用∙标注):
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
3 |
4 |
|
|
7 |
A2 |
|
2−θ ∙ |
2+θ ∙ |
|
4 |
A3 |
|
0+θ ∙ |
3−θ ∙ |
6 |
9 |
收量 |
3 |
6 |
5 |
6 |
|
其中θ=min{x22=2,x33=3∣(2,2)与(3,3)为闭合回路中排列的序号为偶数的格点}=2。因为根据产销平衡原则,进基变量要增大(从0变为正数),所以销售点B2从A3进更多的货(θ吨)就要从A2少进同等数量的货(θ吨),同时闭合回路上的B3也从A3进了货,由于A2、A3的发量是固定的,所以销售点B3要从A3少进θ吨货,同时从A2进θ吨货以保持收量;而θ的取值准则是尽可能大,由于最多只能减到零,所以找(2,2)和(3,3)(也就是θ符号为负)中的最小值。
此时显然易见,选择θ=2后(2,2)变为了零,于是x22就是离基变量。因此改善后的基本可行解为
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
3 |
4 |
|
|
7 |
A2 |
|
|
4 |
|
4 |
A3 |
|
2 |
1 |
6 |
9 |
收量 |
3 |
6 |
5 |
6 |
|
此时再通过位势法检验是否为最优解,若不是,再重复前面的步骤直至得到最优解。
后续过程
第二次迭代
用位势法得到的价格表:
|
B1 |
B2 |
B3 |
B4 |
ui |
A1 |
3 |
11 |
[3]-14 |
[10]-2 |
12 |
A2 |
[1]13 |
[9]13 |
2 |
[8]11 |
-3 |
A3 |
[7]11 |
4 |
10 |
5 |
5 |
vj |
-9 |
-1 |
5 |
0 |
|
检验数还有小于零的部分,继续迭代,对应的调运表:
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
3 |
4−θ ∙ |
0+θ ∙ |
|
7 |
A2 |
|
|
4 |
|
4 |
A3 |
|
2+θ ∙ |
1−θ ∙ |
6 |
9 |
收量 |
3 |
6 |
5 |
6 |
|
用闭合回路法得到的调运表:
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
3 |
3 |
1 |
|
7 |
A2 |
|
|
4 |
|
4 |
A3 |
|
3 |
|
6 |
9 |
收量 |
3 |
6 |
5 |
6 |
|
第三次迭代
再用位势法检验是否为最优解,得到的价格表为:
|
B1 |
B2 |
B3 |
B4 |
ui |
A1 |
3 |
11 |
3 |
[10]-2 |
12 |
A2 |
[1]-1 |
[9]-1 |
2 |
[8]-3 |
11 |
A3 |
[7]11 |
4 |
[10]14 |
5 |
5 |
vj |
-9 |
-1 |
-9 |
0 |
|
检验数还有小于零的部分,继续迭代,对应的调运表:
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
3 |
3−θ ∙ |
1+θ ∙ |
|
7 |
A2 |
|
|
4−θ ∙ |
0+θ ∙ |
4 |
A3 |
|
3+θ ∙ |
|
6−θ ∙ |
9 |
收量 |
3 |
6 |
5 |
6 |
|
…
经过有限次运算得到最优解,且最优解有无穷多个
|
B1 |
B2 |
B3 |
B4 |
发量 |
A1 |
2 |
|
5 |
|
7 |
A2 |
1 |
|
|
3 |
4 |
A3 |
|
6 |
|
3 |
9 |
收量 |
3 |
6 |
5 |
6 |
|
相应的最优值为z=85.
三、整数规划问题
1. 整数规划问题的数学模型
在许多线性规划问题中,决策变量具有不可分割的性质,如人数、设备的台数、车辆船只数、项目的个数等
- 若一规划所有的决策变量都取整数,则称其为纯整数规划问题;
- 若有部分决策变量要求取整数,则称其为混合整数规划问题;
- 若决策变量只能取0或者1,则称其为0-1型整数线性规划问题。
例1:设有一背包的最大容量V0,现有n种物品可供选择装入背包,每种物品数量不限,设第种物品的体积为Vi(i=1,2,⋯,n),相应的价值为ci,问每种物品各取多少件装入,可使背包中价值最大?
解:设第i种物品有xi件装入背包,则maxz=c1x1+c2x2+⋯+cnxn,s.t.{V1x1+V2x2+⋯+Vnxn≤V0x1,x2,⋯,xn≥0且均为整数,这是一个纯整数规划问题。
通常情况下,整数规划无法通过单纯形法求解然后对结果舍零取整得到解,因为取整后的解不一定仍然是可行解,即使是可行解也不一定是最优解。我们把整数规划问题中的整数要求去掉后的问题称为原问题的伴随规划。
例2:某厂计划用集装箱托运甲乙两种货物,每箱的体积、重量、可获利润及托运所受限制见下表,如何设计装运方案可使所获利润最大?
货物/箱 |
体积/m2 |
重量/百斤 |
利润/百元 |
甲 |
5 |
2 |
20 |
乙 |
4 |
5 |
10 |
托运限制/集装箱 |
24 |
13 |
|
解:数学模型为maxz=20x1+10x2,s.t.⎩⎨⎧5x1+4x2≤242x1+5x2≤13x1,x2≥0且为整数,这是一个纯整数规划问题,其伴随规划为maxz=20x1+10x2,s.t.⎩⎨⎧5x1+4x2≤242x1+5x2≤13x1,x2≥0,伴随规划的最优解为x∗=(4.8,0)T,易知舍零取整后的解x1=(5,0)T不是原问题的整数规划的可行解;而x2=(4,0)T是可行解,但不是整数规划的最优解。
2. 分枝定界法
分枝定界法既可以用来求解整数规划,也可以用来求解混合整数规划。分枝定界法的主要思路是:首先求解整数规划的伴随规划,如果求得的最优解不符合整数条件,则增加新约束条件以便缩小可行域,将原整数规划问题分枝为两个子规划,再解子规划的伴随规划,通过求解一系列子规划的伴随规划问题,从而不断地定界,最后得到原整数规划问题的整数最优解。
原问题与伴随规划问题的解的关系如下:
- 如果伴随规划问题无可行解二,则原整数规划问题无可行解;
- 若伴随规划的最优解满足整数规划的整数条件,则该最优解也为原整数规划的最优解;
- 否则使用分支定界法进行计算。
分枝定界法的步骤如下(一般将原问题(整数规划问题)称为问题A0,将它的伴随问题称为B0):
步骤一:计算原问题A0目标函数的初始上界z。因为问题A0的可行域⊂问题B0的可行域,故一般用伴随问题的z∗作为初始上界;
步骤二:计算原问题A0目标函数的初始下界z。若能从问题A0的约束条件观察到一个整数可行解,则可将其目标函数值作为问题A0目标函数的初始下界,否则可以令初始下界为z=−∞。
步骤三:增加约束条件将原问题分枝。一般是分别采用向上、向下取整将伴随问题的小数部分去掉,从而分为两枝。
步骤四:分别求解每一分枝。对某分枝的伴随问题求解,可能出现以下情形:
- 无可行解:说明该分枝情况已查明,不需对其继续分枝,称该分枝为明枝;
- 得到整数最优解:这种情况也不需再进行分枝,也称为明枝;
- 得到不满足整数条件的最优解:
- 该最优解的目标函数值z小于当前下界z,该分枝内不可能含有原问题的整数最优解,故不需对其继续分枝,称为枯枝;
- 该最优解的目标函数值z大于当前下界z,仍需对该枝继续分枝,以查明分枝内是否有目标函数值比当前的z更好的整数最优解。(如果有多对分枝都需要继续分枝时,首先对目标函数值较优的分枝求解,待目标函数较优的那枝全部分解到不能再分,全部查清时为止,再回过头来考虑目标函数稍差的那枝。)
步骤五:修改上界z与下界z。修改的时机和原则如下:
修改下界z的时机:每求出一个整数可行解,须修改下界z。
修改下界z的原则:在至今所求的全部整数可行解中,选目标函数值最大的那个作为最新下界z。
修改上界z的时机:每求解完一对分枝,须考虑修改上界z。
修改上界z的原则:挑选目前为止所有未被分枝问题的目标函数值中最大的一个作为新的上界。
步骤六:结束准则:当所有分枝均已查明,且此时z=z,则已得到了原问题的整数最优解,即目标函数值为下界z的那个整数解。
下面通过例题了解分枝定界法
例:求解下列整数规划问题maxz=10x1+20x2,s.t.⎩⎨⎧5x1+8x2≤60x1≤8x2≤4x1,x2≥0且为整数.
解:为方便表示,这里将原问题称为问题A0,将它的伴随问题称为B0。通过之前的方法可以求出B0的最优解为x0∗=(5.6,4)T,最优值为f0∗=136,下面开始使用分枝定界法:
步骤一:计算原问题A0目标函数的初始上界z。因为问题B0的最优解不满足整数条件,故不为A0的最优解,又因为A0的可行域⊂问题B0的可行域,故有z∗≤f0∗=136,于是z∗的初始上界为z=136。
步骤二:计算原问题A0目标函数的初始下界z。对于本例,一个显然的可行解为零,故z=0。
步骤三:分枝。由于B0的最优解为x0∗=(5.6,4)T,可以将问题分枝为问题A1maxz=10x1+20x2s.t.⎩⎨⎧5x1+8x2≤60x1≤8x2≤4x1≤5x1,x2≥0且为整数和问题A2maxz=10x1+20x2s.t.⎩⎨⎧5x1+8x2≤60x1≤8x2≤4x1≥6x1,x2≥0且为整数。
步骤四、五、六:分别求解这一对分枝的伴随问题,并更新上、下界,一直进行分枝,直至上界与下界相等时得到整数最优解。
整数规划问题-分枝定界法-例题-01
3. 割平面法
割平面法主要用于求解纯整数规划。割平面法的主要思路是:不断切割原问题伴随规划的可行域,使它在不断缩小的过程中,将原问题的整数最优解逐渐暴露,且趋于可行域极点。
下面通过例题了解割平面法
例:求解下列整数规划问题原问题A0maxz=x1+x2s.t.⎩⎨⎧−x1+x2≤13x1+x2≤4x1,x2≥0且为整数、伴随规划B0的标准形maxz=x1+x2s.t.⎩⎨⎧−x1+x2+x3=13x1+x2+x4=4x1,x2,x3,x4≥0.
解:求出伴随问题的最优单纯形表如下
整数规划问题-割平面法-例题-01
割平面法步骤如下:
(1)割平面方程可以由上述最优单纯形表上任一含有不满足整数条件的基变量的约束方程演变得到,比如上表中x1对应的行(第4行),有x1−41x3+41x4=43;
(2)将所选的约束方程中非基变量的系数及常数项进行拆分处理,具体规则是:将上述系数和常数项均拆成一个整数加一个非负的真分数之和,即x1+(−1+43)x3+(0+41)x4=(0+43);
(3)将上述方程重新组合,组合的原则是:将非基变量系数及常数项中的非负真分数移到等号左端,将其他部分移到等式右端,即43x3+41x4−43=0−x1+x3+0⋅x4;
(4)求割平面方程,等式右端全部项为整数,并且松弛变量被认为是非负整数(当原问题A0的约束方程组中的系数或常数项中有非整数时,可以先将它们化为整数再标准化,所以这里只考虑整数的情况)。这里43x3+41x4=43+(x3−x1),故x3−x1≥0,否则43+(x3−x1)<0与x3,x4非负的条件矛盾(反证法:由于x1,x3,x4都是整数,如果x3−x1严格小于零,那么最大只能取−1,此时43−1<0,而如果x3,x4非负,43x3+41x4不可能等于一个小于零的数),因此割平面条件为:43x3+41x4≥43;
(5)把割平面条件化为−43x3−41x4≤−43,加入松弛变量得−43x3−41x4+x5=−43,将此方程加入到伴随规划B0的约束方程中,得到新问题maxz=x1+x2,s.t.⎩⎨⎧−x1+x2+x3=13x1+x2+x4=4−43x3−41x4+x5=−43x1,x2,x3,x4,x5≥0;
新问题的单纯形表可以由原来的单纯形表改写得到,直接添加一行一列如下
整数规划问题-割平面法-例题-02
此时还不是最优单纯形表,通过对偶单纯形法得到最优单纯形表如下
整数规划问题-割平面法-例题-03
可见此时已得到整数最优解,所以最优解为x∗=(1,1)T,且相应的目标函数值为z∗=2。
4. 指派问题与匈牙利算法
最小化的指派问题
指派问题:在实践中经常会遇到这样一种问题:有n项不同的工作或任务,需要n个人去完成,要求每人只完成一项工作,由于每人的知识、能力、经验等不同,故各人完成不同任务所需要的时间不同,问应指派何人完成何项工作,使完成n项工作总耗时最少,这种问题叫做指派问题,指派问题是一种整数规划问题,同时也是一类特殊的运输问题。
指派问题的数学模型为minz=∑i=1n∑j=1ncijxij,s.t.⎩⎨⎧∑j=1nxij=1,∑i=1nxij=1,xij=0或1,i=1,2,⋯,nj=1,2,⋯,ni,j=1,2,⋯,n,其中xij表示令第i个人去完成第j份工作,cij表示第i个人做第j份工作的时间(或资源浪费,匈牙利算法要求cij非负),约束条件表示一个人只能做一份工作,并且每份工作只由一个人来完成。将效率系数cij排成一个n×n矩阵,称为效率矩阵或价值系数矩阵,即C=c11c21⋮cn1c12c22⋮cn2⋯⋯⋯c1nc2n⋮cnn。
定理:设指派问题的效率矩阵为C=(cij)n×n,若将该矩阵的某一行或列的各元素都减去同一个常数t,得到新的效率矩阵C′=(cij′)n×n,则以C′为效率矩阵的新指派问题与原指派问题的最优解相同,但其最优值比原最优值减少t。
证明:根据指派问题定义,假设第k行减去了t,z′=∑i=1n∑j=1ncij′xij=∑i=1, i=kn∑j=1ncij′xij+∑j=1nckj′xkj=∑i=1, i=kn∑j=1ncijxij+∑j=1n(ckj−t)xkj=∑i=1, i=kn∑j=1ncijxij+∑j=1nckjxkj−t∑j=1nxkj=z−t.
定理:若将指派问题的效率矩阵每一行(或每一列)分别减去各自的最小元素,则得到的新指派问题与原指派问题有相同的最优解。
将效率矩阵的每一行减去各行最小元素,所得矩阵的每一列再减去当前列中最小元素,最后的新效率矩阵C′中必然会出现一些零元素,显然,元素cij′=0表示第i个人去干第j项工作效率最好,或表示第j项工作由第i个人来干效率最高。
定义:在效率矩阵C中,有一组处在不同行不同列的零元素,称为独立零元素组,此时其中每个元素称为独立零元素。
例1:由效率矩阵的独立零元素组确定最优指派问题:C=5204035820600070,其中加粗的零就组成了一个独立零元素组,此时将独立零元素组所在的位置赋1,其他为零,就得到了一个解x(1)=0010100000010100,即让第1个人去完成第2份工作、第2个人去完成第4份工作,…;独立零元素组并不唯一,比如x(2)=0010100001000001。
可见,当独立零元素组的个数m和矩阵的阶数n相等时,可以直接利用定理2及独立零元素组写出解。但是在有的问题中,效率矩阵C的独立零元素的个数不到n个,这样就无法求到最优指派方案,需要作进一步分析。
定理:效率矩阵C中独立零元素的最多个数等于能覆盖所有零元素的最少直线数。
比如:C1= 5204035820600070 是4个独立零元素,C2= 5204003586205030070620245 也是4个独立零元素。
圈零法:要找出矩阵的独立零元素,可以使用圈零法,具体步骤如下:
(1)进行行检验:对C进行逐行检查,每行只有一个未标记的零元素时,用O记号将该元素圈起,然后将被圈起的零元素所在列的其他未标记的零元素用记号×划去;
(2)进行列检验:与进行行检验类似;
(3)重复检验,直到每一行都没有未被标记的零元素或至少有两个未被标记的零元素为止。
根据圈零法标记的个数的不同,最后可能会有如下情况:
- 每一行均有一个圈零(即独立零元素)出现,圈零的个数m恰好等于矩阵阶数n,即前面的例1的情况,此时可以直接给出指派方案;
- 存在未标记过的零元素,但它们所在的行或列中,未标记过的零元素至少有两个,此时任选一个作为独立零元素(划圈),其他零元素划×标记(此时方案可能有多个);
- 不存在未被标记过的零元素,但圈零的个数m<n。
匈牙利算法的步骤:
(1)根据定理2,将各行或各列都减去各自的最小元素,得到新效率矩阵C1,此时新效率矩阵对应的最优解仍与原问题相同;
(2)用圈零法找出新矩阵C1中的独立零元素;
(3)进行试指派,根据前面圈零法讨论的不同情况:①情况1可以直接令圈零位置的决策变量取值1,其他决策变量取值0,得到一个最优指派方案;②情况2经过任选一个作为独立零元素的处理,再次进行圈零法的检验,会转变为情况1或情况3;③情况3,需进行第四、五步。
(4)增加独立零元素的个数。首先要作最少直线覆盖当前所有零元素:①对C1中所有不含圈零元素的行打;②对打的行中所有零元素所在列打;③对所有打列中圈零元素所在行打;④重复上述第②、③步,直到不能进一步打为止;然后对所有没有打的行划横线、对已打的列划竖线,便能用最少直线覆盖所有零元素。然后再在未被直线覆盖的元素中找出最小元素,将打行的各元素减去这个最小元素,将打列的各元素加上这个最小元素,这样就能增加独立零元素的个数。
(5)对已增加独立零元素的矩阵,重复(2)、(3)、(4)步,直至出现圈零数m等于矩阵阶数n为止。
下面通过例题了解匈牙利算法的完整求解步骤
例2:现有一个5×5的指派问题,其效率矩阵为C=1287154791714109612677614610969109,求该指派问题。
解:
步骤一:变换效率矩阵,根据定理2,将各行都减去各自的最小元素,得C1=52090031086205030070620245;
步骤二:用圈零法求出新矩阵C1中独立零元素,最后得到C1=52090×03108620×5030×0×70×620245,检查的过程为(这里用代替圈):先看第1、2行,都有多于一个的零,不作任何标记;到第3行时,只有一个零,对这个零作标记,此时同列的(5,1)位置的零划掉;再到后面4、5行都不符合,不作标记;行检查完一次,开始对列进行检查,第1列的所有零都已经有标记了,不作检查;第2列只有一个零,对这个零作标记,此时同行的(1,4)位置的零划掉;第3、4列都有多于一个的零,不作任何标记;第5列只有一个零,对这个零作标记,此时同行的(2,3)和(2,4)位置的零划掉;第二次行检查可以看到,只剩下(4,3)和(4,4)位置的零还没有被标记,任选一个画,另一个划×即可。
步骤三:进行试指派。此时显然独立零元素的个数4小于矩阵阶数5,属于情况3,进入步骤四。
步骤四:先作最少直线覆盖当前所有零元素,根据规则对行、列打,然后对所有没有打的行划横线、对已打的列划竖线,最后得到C1=52090×03108620×5030×0×70×620245,然后找出没有划线元素中最小元素为2,打行减去2、打列加上2,得到C2=74011003884203010050420043。
步骤五:找出C2的独立零元素,用圈零法得C2=740×1100388420×3010×050×420×043,此时可以直接写出最优解为x12=x24=x35=x43=x51=1,其余xij=0。
(补充)最大化的指派问题
对于最大化的指派问题,maxz=∑i=1n∑j=1ncijxij,s.t.⎩⎨⎧∑j=1nxij=1,∑i=1nxij=1,xij=0或1,i=1,2,⋯,nj=1,2,⋯,ni,j=1,2,⋯,n,不能使用将目标函数添加负号转化为最小化问题的方法,因为匈牙利算法要求每个元素都非负。
我们可以在效率矩阵C=(cij)n×n中,找出一个最大元素M,并令bij=M−cij(i,j=1,2,⋯,n),这样就能转化为求解最小化问题minz=∑i=1n∑j=1nbijxij,s.t.⎩⎨⎧∑j=1nxij=1,∑i=1nxij=1,xij=0或1,i=1,2,⋯,nj=1,2,⋯,ni,j=1,2,⋯,n。
这两个问题等价的证明:z′=∑i=1n∑j=1nbijxij=∑i=1n∑j=1n(M−cij)xij=M∑i=1n∑j=1nxij−∑i=1n∑j=1ncijxij=nM−z,即z′与z负相关,z′的最小值就是z的最大值。
四、非线性规划
1. 非线性规划的数学模型
非线性规划的概念:线性规划的一个明显的特点是其目标函数与约束条件都是决策变量的一次函数,如果在一个规划问题的目标函数和约束条件中,至少有一个方程是决策变量的非线性函数,则称这类规划问题为非线性规划。
例1(选址问题):设有n个市场,第k个市场的位置为(pk,qk),它对某种货物的需求量为bk(k=1,2,⋯,n),现计划建立m个仓库,第l个仓库的存储容量为al(l=1,2,⋯,m)。试确定仓库的位置使各仓库对各市场的运输量与路程乘积之和为最小。
解:设第l个仓库的位置为(xl,yl),它到第k个市场的货物供应量为zlk,则数学模型为min∑l=1m∑k=1nzlk(xl−pk)2+(yl−qk)2,s.t.⎩⎨⎧∑k=1nzlk≤al,∑l=1mzlk=bk,zlk≥0,l=1,2,⋯,mk=1,2,⋯,nl=1,2,⋯,m;k=1,2,⋯,n.
一般非线性规划的数学模型可表示为minf(x),s.t.{gi(x)≥0,hj(x)=0,i=1,2,⋯,mj=1,2,⋯,l,其中x=(x1,x2,⋯,xn)T∈Rn,而f,gi,hj都是Rn→R1的映射。
定义:我们称下列集合H={gi(x)≥0 (i=1,2,⋯,m),hj(x)=0 (j=1,2,⋯,l)}为非线性规划问题的可行域。当非线性规划问题的自变量x没有任何约束,或说可行域是整个n维向量空间H=Rn,则称该非线性规划问题为无约束问题,记作minf(x)或minx∈Rnf(x),有约束问题与无约束问题在处理方法上有明显的不同。
在微积分中曾讨论过极值问题,它们都属于无约束极值问题,即使有约束,也是等式约束,这样可利用拉格朗日乘子法。将等式约束极值问题化为无约束极值问题来求解,这种极值问题统称为经典极值问题。
定义:若x∗∈H,且满足minx∈Hf(x)=f(x∗),即对∀x∈H都有f(x∗)≤f(x),则称x∗为非线性规划的全局最优解。
定义:若x∗∈H,且存在一个δ>0,对于∀x∈N(x∗,δ)∩H(这里的N(x∗,δ)表示x∗的一个邻域),使minx∈N(x∗,δ)∩Hf(x)=f(x∗)成立,即∀x∈N(x∗,δ)∩H,都有f(x∗)≤f(x),则称x∗为非线性规划的局部最优解。
若线性规划问题有最优解,则其最优解一定在可行域的极点上,但非线性规划的最优解却可能是可行域上的任何一点。
线性规划问题的最优解一定是全局最优解,而非线性规划问题有局部最优解与全局最优解之分,一般的非线性规划算法往往求出的是局部最优解。
2. 无约束问题-最优性条件
定义:设u=f(x),x∈S⊆Rn,若在x0=(x1(0),x2(0),⋯,xn(0))T处对于自变量x=(x1,x2,⋯,xn)T的各分量的偏导数∂xi∂f(x0)都存在,则称函数u=f(x)在点x0处一阶可导,并称向量∇f(x0)=(∂x1∂f(x0),∂x2∂f(x0),⋯,∂xn∂f(x0))T是u=f(x)在点x0处的梯度或一阶导数。
定理(微分与梯度的关系):设u=f(x),x∈S⊆Rn,若f在点x0处可微,则f在点x0处的梯度存在,并且有d f(x0)=∇f(x0)TΔx,其中Δx=(x1−x1(0),x2−x2(0).⋯,xn−xn(0))T为x0处的增量。
定义:设u=f(x),x0∈S⊆Rn,若f在x0∈S处对于自变量x∈S的各分量的二阶偏导数∂xi∂xj∂2f(x0)(i,j=1,2,⋯,n)都存在,则称函数在f(x)在点x0处二阶可导,且称矩阵∇2f(x)=∂x12∂2f(x0)∂x2∂x1∂2f(x0)⋮∂xn∂x1∂2f(x0)∂x1∂x2∂2f(x0)∂x22∂2f(x0)⋮∂xn∂x2∂2f(x0)⋯⋯⋱⋯∂x1∂xn∂2f(x0)∂x2∂xn∂2f(x0)⋮∂xn2∂2f(x0)为f(x)在点x0处的二阶导数或黑塞矩阵,记作H(x0)。
定义:设f(x)是一个多元函数,x∈S⊆Rn,若∃x∗∈S,且∃δ>0,则:
- 如果对∀x∈N(x∗,δ)∩S,都有f(x∗)≤f(x),则称x∗为f(x)的局部极小点,称f(x∗)是局部极小值;
- 如果对∀x∈N(x∗,δ)∩S,x=x∗,都有f(x∗)<f(x),则称x∗为f(x)的严格局部极小点,称f(x∗)是严格局部极小值。
定义:设f(x)是一个多元函数,x∈S⊆Rn,若∃x∗∈S,则:
- 对∀x∈S,都有f(x∗)≤f(x),则称x∗为f(x)的整体极小点,称f(x∗)是整体极小值;
- 对∀x∈S,x=x∗,都有f(x∗)<f(x),则称x∗为f(x)的严格整体极小点,称f(x∗)是严格整体极小值。
多元函数局部极小点(极大点)的判定条件:
定理(一阶必要条件):设f(x)(x∈S⊆Rn)在点x∗∈S处可微,若x∗是f(x)的局部极值点,则∇f(x∗)=0。
注意∇f(x∗)=0只是函数取得极值点的必要条件;
我们将满足∇f(x)=0的点称为函数的平稳点;函数的平稳点或是其极小点,或是其极大点,或两者都不是。
定理(二阶必要条件):设f(x)(x∈S⊆Rn)在点x∗∈S处二次可微,若x∗是f(x)的局部极小点,则∇f(x∗)=0,且∇2f(x∗)半正定。
定理(二阶充分条件):设f(x)(x∈S⊆Rn)在点x∗∈S处二次可微,若∇f(x∗)=0,且∇2f(x∗)正定,则x∗是函数f(x)的严格局部极小点。
定理(二阶充分条件):设f(x)(x∈S⊆Rn)在点x∗∈Rn的一个邻域N(x∗,δ)内二次可微,若在点x∗处满足∇f(x∗)=0,且∀x∈N(x∗,δ)都有∇2f(x)半正定,则x∗是函数f(x)的局部极小点。
例1:求下列函数的极小点f(x)=(x12−1)2+x12+2x22−2x1.
解:先找出平稳点{fx1(x)=4x13−2x1−2=0fx2(x)=4x2=0,得出平稳点为x∗=(1,0)T,然后再判断黑塞阵的正定性∇2f(x)=[12x12−2004],在平稳点的黑塞阵为∇2f(x∗)=[10004],观察可知它是一个正定矩阵(因为特征值都大于零),根据前面的定理(二阶充分条件1)知,x∗=(1,0)T是严格局部极小点。
例2:求下列函数的极小点f(x)=(x1−2)4+(x1−2x2)2.
解:先找出平稳点{fx1(x)=4(x1−2)3+2(x1−2x2)=0fx2(x)=−4(x1−2x2)=0,得出平稳点为x∗=(2,1)T,黑塞阵为∇2f(x)=[12(x1−2)2+2−4−48],在平稳点的黑塞阵为∇2f(x)=[2−4−48],该矩阵是半正定的(一阶顺序主子式大于零,而二阶顺序主子式等于零),考虑x∗的邻域x=x∗+Δx=(2+Δx11+Δx2),则H(x)=[12Δx12+2−4−48],易证上述矩阵是半正定的,故x∗是局部极小点。
3. 解非线性规划的基本思路
引入:前面讨论了无约束问题的最优性条件,但是对大多数实际问题,要用最优性条件来求解非线性规划是很困难的:有的问题导数不存在,有的问题导数即使存在,计算也很麻烦;多数问题由条件∇f(x)=0得到的是一个非线性方程组,求解非常困难,甚至根本无法得到解析解。
求解非线性规划问题一般采用数值计算的选代方法。迭代就是从某个已知点x(k)出发,按照某种算法求出后继点x(k+1),用k+1代替k,重复以上过程,得到一个点列{x(k)}。非线性规划选代方法的基本思想是:使得要求的送代算法(1)当{x(k)}是穷点列时,其最后一点是该问题的最优解;(2)当{x(k)}是无穷点列时,它的极限点是该问题的最优解。
定义:设x(k),x(k+1)∈Rn是某种迭代算法的第k轮和第k+1轮迭代点,记Δxk为两者之差,即Δxk=x(k+1)−x(k),故x(k+1)=x(k)+Δxk,令p(k)是向量Δxk方向上的单位向量,则有Δxk=λkp(k),进而有x(k+1)=x(k)+λkp(k),其中λk是大于零的实数。这里的x(k+1)=x(k)+λkp(k)就是基本迭代格式,并且从该选代格式可知,求解非线性规划问题的关键在于:如何构造每一轮的搜索方向p(k)和步长因子λk。
根据基本选代格式,求解非线性规划问题的选代算法的关键在于构造搜索方向,构造不同的搜索方向,就是不同的算法。但这些不同的搜索方向有一个共同的要求:就是当算法从x(k)产生了搜索方向p(k)和步长λk,从而得到了x(k+1)=x(k)+λkp(k),要求f(x(k+1))=f(x(k)+λkp(k))<f(x(k)),也就是说,搜索方向p(k)应指向目标函数值减小的方向。
定义:设f:Rn→R1,点x∈Rn,向量p∈Rn(p=0),若存在一个数δ>0,使∀λ∈(0,δ)都有f(x+λp)<f(x),则称向量p是f(x)在点x处的下降方向。
一般来说,f(x)在点x处的下降方向不止一个,可能有无穷多个,那么如何在点x处选择使函数下降最快的方向是一个重要的问题。并且对于有约束的非线性规划问题minx∈Hf(x),算法不仅要保证搜索方向为下降方向,而且要保证x(k)+λkp(k)仍在可行域H内,因此对于有约束的非线性规划问题,要寻找可行下降方向。
计算的终止条件:常用的计算终止条件有:
- 当自变量的改变量充分小时终止计算∣∣x(k+1)−x(k)∣∣<ϵ1,或相继两次迭代的相对误差充分小时终止计算∣∣x(k)∣∣∣∣x(k+1)−x(k)∣∣<ϵ2;
- 当函数值的下降量充分小时终止计算f(x(k))−f(x(k+1))<ϵ3,或相继两次迭代的相对误差充分小时终止计算∣f(x(k))∣f(x(k))−f(x(k+1))<ϵ4;
- 对于无约束优化问题,当函数梯度的模充分小时终止计算∣∣∇f(x(k+1))∣∣<ϵ5;
非线性规划迭代算法的一般步骤:
(1)选取初始点x(0),令k:=0;
(2)构造搜索方向(下降方向或可行下降方向,构造的方向不同,就是不同的算法);
(3)确定步长因子λk,此时f(x(k)+λkp(k))只是λ的一元函数,求步长λk就相当于求一元函数f(x(k)+λkp(k))的极小点,即f(x(k)+λkp(k))=minλ>0f(x(k)+λp(k)),称上述λk为最优步长,并称求最优步长的过程为一维搜索;
(4)求出下一个迭代点x(k+1)=x(k)+λkp(k);若x(k+1)已满足规定的终止条件,停止迭代,否则令k:=k+1,转第二步。
4. 无约束问题-步长(一维搜索)
前面提到过,求解f(x(k)+λkp(k))=minλ>0f(x(k)+λp(k))的过程称作一维搜索,因为λk是φ(λ)=f(x(k)+λp(k))的极小点,故在理论上λk应满足dλdφ(λ)∣λ=λk=0,但是对许多问题来说,求导计算并不容易,且不便于利用计算机来解决,因此一般不是用求导的方法来求λk的精确值,而是求它的近似值。求解λk近似值的方法主要分为区间收缩法、函数逼近法。
后面介绍的黄金分割法属于区间收缩法,而牛顿法和抛物线法属于函数逼近法。
黄金分割法
黄金分割法也称0.618法,属于区间收缩法,其主要步骤为:①找出包含极小点的初始搜索区间;②按黄金分割点通过对函数值的比较不断缩小搜索区间(要保证极小点始终在搜索区间内);③当区间长度小到精度范围之内时,可以粗略地认为区间端点的平均值即为极小点的近似值。
定义:设函数φ(λ):R1→R1,闭区间[a0,b0]⊂R1,若存在一点λ∗∈[a0,b0],使φ(λ)在[a0,λ∗]上严格递减,在[λ∗,b0]上严格递增,则称φ(λ)为[a0,b0]的单谷函数,[a0,b0]为φ(λ)的单谷区间。
单谷函数及其性质::若φ(λ)是单谷区间[a0,b0]上的单谷函数,极小点为λ∗,在[a0,b0]中任取两点a1,b1,且a1<b1,则
- 当φ(a1)<φ(b1)时,λ∗∈[a0,b1];
- 当φ(a1)>φ(b1)时,λ∗∈[a1,b0];
非线性规划-一维搜索-黄金分割法-01
黄金分割法的步骤:
(1)在单谷区间[a0,b0]中找两个试点x1,x1′后,比较φ(x1)与φ(x1′)的大小,就可把搜索区域缩小至[a1,b1];
(2)在区间[a1,b1]中继续找试点x2,x2′,比较φ(x2)与φ(x2′)的大小,又可以把搜索区域从[a1,b1]缩小至[a2,b2];
(3)重复上述过程,直到[ak,bk]的区间长度足够小,并达到事先规定的精度时,取λ∗≈21[ak+bk]。
x1和x1′一般取黄金分割点及对称点。
非线性规划-一维搜索-黄金分割法-02
黄金分割点及对称点:若记∣a0b0∣=l,∣a0x1∣=c,则∣x1b0∣=l−c,若希望∣a0b0∣∣a0x1∣=a0x1x1b0⇔lc=ll−c⇔lc=25−1,因此x1=a0+0.618(b0−a0),x1′=a0+0.382(b0−a0)。
取黄金分割点及对称点的优势:
- 若φ(x1′)<φ(x1),搜索区间被保留下来的是[a0,x1],则x1′是[a0,x1]区间的黄金分割点;
- 若φ(x1)<φ(x1′),搜索区间被保留下来的是[x1′,b0],则其中的x1是[x1′,b0]区间的对称点。
黄金分割法的简化步骤:
(1)选取初始数据:确定初始搜索区间[a0,b0],给出最后的区间精度δ>0;
(2)计算初始的两个试点x1,x1′:在区间[ak,bk]上黄金分割点用xk+1表示,其对称点用xk+1′表示,计算{x1=a0+0.618(b0−a0)x1′=a0+0.382(b0−a0),并计算相应试点的函数值φ(x1)和φ(x1′),此时规定k:=0;
(3)比较目标函数值:
- 当φ(xk+1′)≤φ(xk+1)时,取{ak+1=akbk+1=xk+1以缩小搜索区间,此时:
- 若Δ=b0−a0bk+1−ak+1<δ,算法终止,则λ∗=21[ak+1+bk+1];
- 若Δ=b0−a0bk+1−ak+1>δ,则要计算新的试点{xk+2=xk+1′xk+2′=ak+1+0.382(bk+1−ak+1),计算φ(xk+2)和φ(xk+2′),算法在k:=k+1转至第三步。
- 当φ(xk+1′)>φ(xk+1)时,取{ak+1=xk+1′bk+1=bk以缩小搜索区间,此时:
- 若Δ=b0−a0bk+1−ak+1<δ,算法终止,则λ∗=21[ak+1+bk+1];
- 若Δ=b0−a0bk+1−ak+1>δ,则要计算新的试点{xk+2′=xk+1xk+2=ak+1+0.618(bk+1−ak+1),计算φ(xk+2′)和φ(xk+2),算法在k:=k+1转至第三步。
例:用黄金分割法求minx≥0φ(x)=x3−2x+1的近似最优解,设初始的搜索区间为[0,3],精度b0−a0bk−ak≤δ=0.238.
解:初始数据题目已经给出;接下来计算初始的两个试点x1,x1′,计算{x1=a0+0.618(b0−a0)=0+0.618(3−0)=1.854x1′=a0+0.382(b0−a0)=0+0.382(3−0)=1.146,计算得φ(x1)=3.6648和φ(x1′)=0.2131,此时规定k:=0;显然φ(x1′)≤φ(x1),取{a1=a0=0b1=x1=1.854,此时Δ=b0−a0b1−a1=3−01.854−0=0.618>0.238=δ,精度不满足要求,故需继续寻找试点{x2=x1′=1.146x2′=a1+0.382(b1−a1)=0+0.382(1.854−0)=0.708,计算得φ(x2)=0.2131和φ(x2′)=−0.0611,算法在k:=k+1转至第三步。
显然φ(x2′)=−0.0611≤0.2131=φ(x2),取{a2=a1=0b2=x2=1.146,此时Δ=b0−a0b2−a2=3−01.146−0=0.382>0.238=δ,精度不满足要求,故需继续寻找试点{x3=x2′=0.708x3′=a2+0.382(b2−a2)=0+0.382(1.146−0)=0.438,计算得φ(x3)=0.0611和φ(x3′)=0.2080,算法在k:=k+1转至第三步。
显然φ(x3′)=0.2080>−0.0611=φ(x3),取{a3=x3′=0.438b3=b2=1.146,此时Δ=b0−a0b3−a3=3−01.146−0.438=0.236>0.238=δ,精度满足要求,因此近似最优解为x∗=2b3+a3=21.146+0.438=0.792,相应的近似最优值为φ(x∗)=φ(0.792)≈−0.087。
牛顿法
牛顿法是一种函数通近法,它的基本思想是:在极小点附近用二阶泰勒多项式近似代替目标函数f(x),从而求出f(x)极小点的估计值。
如minf(x),x∈R1,现已有f(x)极小点的第k级估计值x(k),在x(k)点将f(x)作二阶泰勒展开:f(x)=f(x(k))+f′(x(k))(x−x(k))+21f′′(x(k))(x−x(k))2+o(∣x−x(k)∣2),记φ(x)=f(x(k))+f′(x(k))(x−x(k))+21f′′(x(k))(x−x(k))2,则在x(k)附近可用φ(x)来近似f(x),即φ(x)≈f(x),因此可以用φ(x)的极小点来近似f(x)的极小点(假定f(x)的极小点在x(k)附近)。
首先要找出φ(x)的驻点x(k+1),φ′(x)=f′(x(k))+f′′(x(k))(x−x(k)),令φ′(x)=0,得驻点为x(k+1)=x(k)−f′′(x(k))f′(x(k)),以φ(x)的驻点x(k+1)作为f(x)在x(k)附近的极小点的第k+1级估计值。
牛顿法的步骤:
(1)给定初始点x(1),给定精度ϵ>0,令k:=1;
(2)计算f′(x(k))及f′′(x(k)):若∣f′(x(k))∣<ϵ,则停止迭代,输出近似极小点x(k),否则进行第三步;
(3)用迭代公式计算x(k+1):x(k+1)=x(k)−f′′(x(k))f′(x(k)),令k:=k+1,返回第二步。
例:用牛顿法求函数minf(x)=∫0xarctantdt的极小点,取x(1)=1,ϵ=0.01.
解:f′(x)=arctanx,f′′(x)=(1+x2)1
k |
x(k) |
f′(x(k)) |
f′′(x(k))1 |
ϵ |
x(k+1) |
1 |
1 |
0.7854 |
2 |
0.7854 |
-0.5708 |
2 |
-0.5708 |
-0.5187 |
1.3258 |
0.5187 |
0.1169 |
3 |
0.1169 |
0.1164 |
1.0137 |
0.1164 |
-0.00106 |
4 |
-0.00106 |
0.0010 |
|
|
|
加步探索法
黄金分割法及其他一维搜索方法都要事先给定一个包含极小点的初始搜索区间,加步探索法能解决这个问题。
加步探索法本身是一种区间试探法,其主要思路就是从一点出发,按一定的步长,试图确定出函数值呈现“高-低-高”的三点。首先从一个方向去找,若不成功,就退回来,再沿相反方向寻找,若方向正确,则加大步长进行探索,最终找到x1,x2,x3三点,并满足x1<x2<x3,f(x1)>f(x2),f(x2)<f(x3)。
加步探索法的步骤:
(1)给定初始点x1,初始步长h0>0;
(2)用加倍步长的外推法找出初始区间,由初始点x1向某个方向,比如x向增大方向走一步,步长为h0,得x2=x1+h0,计算f(x1),f(x2)并比较函数值的大小:
- 若f(x2)<f(x1),说明方向选对,则步长加倍,继续往前走,有x3=x2+2h0,若仍有f(x3)<f(x2),则将步长再加倍,有x4=x3+4h0,直到xk点函数值刚刚变为增加为止,则得三点xk−2<xk−1<xk,且相应的函数值呈“高-低-高”的形式,即f(xk−2)>f(xk−1),f(xk−1)<f(xk),故极小点必在区间[xk−2,xk]上。
- 若f(x2)>f(x1),说明方向选错,因此仍由x1点出发向相反方向走一步,有x3=x2−h0,若此时有f(x3)<f(x2),说明方向对,仍需在此方向上加步迈进,取x4=x3−2h0,若仍有f(x4)<f(x3),再加倍步长继续同方向迈出,直到函数值刚刚变为增加为止,这样就得到了三点xk<xk−1<xk−2,且f(xk)>f(xk−1),f(xk−1)<f(xk−2),故极小点必在[xk,xk−2]上。
(3)再进一步缩小搜索区间,在上述三个点xk−2,xk−1,xk之间,步长是逐次加倍的,故有xk−xk−1=2(xk−1−xk−2),若在xk−1与xk之间再插入一点xk+1,令xk+1=21(xk−1+xk),得到等间距的四个点xk−2,xk−1,xk+1,xk,比较f(xk−2),f(xk−1),f(xk+1),f(xk),令函数值最小的点为x2,它的左、右邻点分为称为x1与x3,则得到更小的搜索区间[x1,x3],这三点有x1<x2<x3,且f(x1)>f(x2),f(x2)<f(x3)。
例:用加步探索法确定minx≥0f(x)=x3−2x+1的搜索区间,要求选取x1=0,h0=1,步长的倍数α=2.
解:取x1=0,令x2=x1+h0=0+1=1,由于f(x1)=1>f(x2)=0,说明方向正确,取x3=x2+2h0=1+2=3,则f(x3)=22>f(x2)=0,即对x1,x2,x3三点有x1<x2<x3且f(x1)>f(x2),f(x2)<f(x3),在x2与x3间插入点x4=21(x2+x3)=2,且f(x4)=5,比较f(x1),⋯,f(x4)知,f(x2)=0为最小值,故[x1,x4]=[0,2]为初始搜索区间。
抛物线法
抛物线法的步骤:
(1)在极小点附近,用二次三项式φ(x)逼近目标函数f(x),若f(x)有x1<x2<x3这三个点。且满足f(x1)>f(x2),f(x2)<f(x3),令二次三项式φ(x)=ax2+bx+c,设⎩⎨⎧φ(x1)=ax12+bx1+c=f(x1)φ(x2)=ax22+bx2+c=f(x2)φ(x3)=ax32+bx3+c=f(x3)⇒a,b,c⇒φ(x);
(2)解出φ(x)的表达式后,通过φ′(x)=2ax+b=0求得φ(x)的驻点x=−2ab=21(x2−x3)f(x1)+(x3−x1)f(x2)+(x1−x2)f(x3)(x22−x32)f(x1)+(x32−x12)f(x2)+(x12−x22)f(x3),用抛物线φ(x)的极小点x来近似f(x)的极小点x(k),即x(k)=x;
(3)然后从x1,x2,x3,x(k)中选出三个点,选择的原则是以目标函数值最小的点作为新的x2点,其左右两个邻点分别作为新的x1及x3点,将新得到的x1,x2,x3点及新的函数值f(x1),f(x2),f(x3)代入公式,可求出极小点新的估计值x(k+1),继续重复该过程,可以得到一个点列{x(k)},这个点列可收敛于原问题的极小点。
注意:三个初始点x1,x2,x3的函数值须满足“高-低-高”的形式,这样才能保证二次三项式φ(x)的二次项系数a>0,且f(x)及φ(x)的极小点都在区间[x1,x3]之内,而寻找初始的这三点,可以用加步探索法。
5.无约束问题-搜索方向
前面介绍了一维搜索的几种方法,本节介绍构造无约束问题搜索方向的方法:
- 一类是在计算过程中只用到目标函数值,不用计算导数,通常称为直接搜索法;
- 一类是要用到目标函数的导数计算,称为解析法。
后面要介绍的变量轮换法属于直接法,而最速下降法、牛顿法、共轭方向法等属于解析法。
变量轮换法
变量轮换法是把多变量函数的优化问题转化为一系列单变量函数的优化问题来求解。它的基本思路是:认为有利的搜索方向是各坐标轴的方向,因此它轮流按各坐标轴的方向搜索最优点。
从某一个给定点出发,按第i个坐标轴xi的方向搜索时,在n个变量中,只有单个变量xi在变化,其余n−1个变量都取给定点的值保持不变。依次从x1到xn做n次单变量的一维搜索,完成了变量轮换法的依次选代.
变量轮换法的步骤:
(1)给定初始点x(1)∈Rn;
(2)从x(1)出发,先沿第1个坐标轴方向e1进行一维搜索,记求得的最优步长为λ1,则可得到新点x(2),即{x(2)=x(1)+λ1e1f(x(2))=f(x(1)+λ1e1)=minλf(x(1)+λe1);再以x(2)为起点,沿第2个坐标轴方向e2进行一维搜索,求得的最优步长为λ2,则可得到新点x(3),{x(3)=x(2)+λ2e2f(x(3))=f(x(2)+λ2e2)=minλf(x(2)+λe2);依此类推直至n个坐标轴方向全部搜索一遍,最后可得到点x(n+1),{x(n+1)=x(n)+λnenf(x(n+1))=f(x(n)+λnen)=minλf(x(n)+λen),此时完成了变量轮换法的一次迭代;
(3)令x(1):=x(n+1),返回第二步,即以xn+1点作为起点,再沿着各坐标轴方向依次进行一维搜索,一直到所有最新点x(n+1)满足给定的精度为止。
例:用变量轮换法求解minf(x)=3x12+2x22+x32,已知初始点x(1)=(1,2,3)T,当∣∣x(n+1)−x(1)∣∣<0.01时停止迭代.
解:第一次迭代:(1)先从初始点x(1)出发,沿x1轴方向e1进行一维搜索:x(1)+λe1=123+λ100=1+λ23,f(x(1)+λe1)=3(1+λ)2+2×22+32=3(1+λ)2+17=φ1(λ),由φ′(λ)=6λ+6=0求得步长λ1=−1,故x(2)=023;(2)然后从x(2)出发,沿x2轴方向e2进行一维搜索:x(2)+λe2=023+λ010=02+λ3,f(x(2)+λe2)=2(2+λ)2+9=φ2(λ),求得步长λ2=−2,故x(3)=003;(3)再从x(3)出发,沿x3轴方向e3进行一维搜索,x(3)+λe3=003+λ001=003+λ,f(x(3)+λe3)=(3+λ)2,求得步长λ3=−3,故x(4)=000;因为∣∣x(1)−x(4)∣∣>0.01,故令x(1):=x(4),再进行新一轮迭代。
第二次迭代:(1)x(1)+λe1=000+λ100=λ00,由f(x(1)+λe1)=3λ2,求得步长λ1=0,故x(2)=000;(2)x(2)+λe2=000+λ010=0λ0,由f(x(2)+λe2)=2λ2,求得步长λ2=0,故x(3)=000;(3)同理,x(3)=000;因为∣∣x(1)−x(4)∣∣=0<0.001,故x(4)=(0,0,0)T为极小点。
变量轮换法的缺点是收敛速度较慢,搜索效率较低。只有对那些具有特殊结构的函数使用起来尚好,但变量轮换法的基本思路非常简单:沿各坐标轴的方向进行搜索。
此方法把一维搜索放到了确定搜索方向的前面,并没有按照前面非线性规划基本思路里的顺序来进行。
最速下降法
最速下降法又称梯度法,它是许多非线性规划算法的一个基础。
最速下降法的原理:考虑问题minf(x),x∈Rn,f(x)∈R1,假设式中f(x)有一阶连续偏导数,有极小点x∗,若现已求得x∗的第k此近似值x(k),欲求得度k+1次近似值x(k+1),需选定搜索方向p(k)。设x=x(k)+λp(k),对f(x)在x(k)点作一阶泰勒展开f(x)=f(x(k)+λp(k))=f(x(k))+λ∇f(x(k))Tp(k)+o(∣∣λp(k)∣∣),其中o(∣∣λp(k)∣∣)=o(λ)是比λ高阶的无穷小量,而∣∣p(k)∣∣=1,故f(x(k+1))−f(x(k))≈λ∇f(x(k))Tp(k),由于要找极小点,所以f(x(k+1))比f(x(k))小,所以左边一定是负数,而右边的步长λ>0,所以∇f(x(k))与p(k)夹角的cos值一定小于零,因为∇f(x(k))Tp(k)=∣∣∇f(x(k))∣∣ ∣∣p(k)∣∣cosθ,其中θ为向量∇f(x(k))与向量p(k)之间的夹角,显然只有当θ=180°(即p(k)的方向与∇f(x(k))相反)时目标函数值下降最快,因为∇f(x(k))Tp(k)=∣∣∇f(x(k))∣∣ ∣∣p(k)∣∣cosθ=∣∣∇f(x(k))∣∣ ⋅1⋅cos180=−∣∣∇f(x(k))∣∣=−∇f(x(k))T∇f(x(k)),故p(k)=−∇f(x(k))。
最速下降法的步骤:
(1)给定初始数据:起始点x(0),给定终止误差ϵ>0,令k:=0;
(2)求梯度向量模的值∣∣∇f(x(k))T∣∣:①若∣∣∇f(x(k))∣∣<ϵ,停止计算,此时x(k)就是极小点的近似值;②若∣∣∇f(x(k))T∣∣>ϵ,转下一步;
(3)构造负梯度方向:p(k)=−∇f(x(k));
(4)进行一维搜索:无论用哪种方法求得λk后,令x(k+1)=x(k)+λkp(k)=x(k)−λk∇f(x(k)),令k:=k+1,转第二步。
对于一些简单的函数,也可以通过直接求导获得最优步长因子;当搜索方向为最速下降方向(即−∇f(x(k)))时,以x(k)为起点的最优步长公式为λk=∇f(x(k))TH(x(k))∇f(x(k))∇f(x(k))T∇f(x(k))(其中H(x(k))是黑塞阵)。
最优步长的推导:设f(x)有二阶连续偏导数,将它在x(k)点作二阶泰勒展开:f(x(k)−λ∇f(x(k)))=f(x(k))+∇f(x(k))T(−λ∇f(x(k)))+21(−λ∇f(x(k)))TH(x(k))(−λ∇f(x(k)))+o(∣∣λ∇f(x(k))∣∣2),将上式中的主要部分记为H(λ),即H(λ)=f(x(k))−λ∇f(x(k))T∇f(x(k))+21λ2∇f(x(k))TH(x(k))∇f(x(k)),函数H(λ)的唯一极值点即为最优步长。
例:用最速下降法求下述函数的极小点f(x)=(x1−1)2+(x2−1)2,初始点x(0)=(0,0)T,ϵ=0.01.
解:梯度向量模的值∇f(x)=[2x1−22x2−2],在x(0)=[00]有∇f(x(0))=[−2−2],且∣∣∇f(x(0))∣∣=8>ϵ,精度不满足要求,要进行搜索,先确定搜索方向p(0)=−∇f(x(0))=−[−2−2]=[22],再确定步长,这里使用前面的公式计算,H(x(0))=∇(∇f(x))=[2002],故λ0=∇f(x(0))TH(x(0))∇f(x(0))∇f(x(0))T∇f(x(0))=[−2−2][2002][−2−2][−2−2][−2−2]=21,因此下一个近四点为x(1)=x(0)+λ0p(0)=[00]+21[22]=[11],此时∣∣∇f(x(1))∣∣=∣∣[00]∣∣=0<ϵ=0.01,故f(x)的极小点为x∗=x(1)=[11].
最速下降法对初始点的选择要求不高,每一轮选代工作量较少,它可以比较快地从初始点达到极小点附近。但在接近极小点时,最速下降法却会出现锯齿现象,选代产生的点列{x(k)}所循路径是之字形的,每次选代移动的步长很小,收敛速度很慢。因此常常将梯度法与其他方法结合起来使用(比如与牛顿法结合),前期用最速下降法,而当接近极小点时改用牛顿法。
牛顿法
牛顿法是一种函数通近法,它的基本思想是:在极小点附近用x(k)点的二阶泰勒多项式来近似目标函数f(x),并用选代点x(k)处指向近似二次函数的极小点方向作为搜索方向p(k)。
牛顿法的原理:设规划问题minf(x),x∈Rn,其中f(x)在点x(k)处具有连续偏导数,黑塞阵∇2f(x(k))正定,f(x)=f(x(k))+∇f(x(k))T(x−x(k))+21(x−x(k))T∇2f(x(k))(x−x(k))+o(∣∣x−x(k)∣∣2),记上式中的主要部分为Q(x)=f(x(k))+∇f(x(k))T(x−x(k))+21(x−x(k))T∇2f(x(k))(x−x(k)),在x(k)附近可用Q(x)来近似f(x),故可用Q(x)的极小点来近似f(x)的极小点,先求出Q(x)的驻点:∇Q(x)=∇f(x(k))+∇2f(x(k))(x−x(k)),由∇Q(x)=0得Q(x)的平稳点x(k+1)=x(k)−[∇2f(x(k))]−1∇f(x(k)),故p(k)=x(k+1)−x(k)=−[∇2f(x(k))]−1∇f(x(k))(由x(k)指向x(k+1)的向量就是x(k+1)−x(k)),λk=1(对比上式与基本迭代格式中的x(k+1)=x(k)+λkp(k)得出步长为1)。
牛顿法的步骤:
(1)选取初始数据:起始点x(0),终止条件ϵ>0,令k:=0;
(2)求梯度向量∇f(x(k)),并计算∣∣∇f(x(k))∣∣:若∣∣∇f(x(k))∣∣<ϵ,停止计算,此时x(k)为极小点;否则转下一步;
(3)构造牛顿方向:p(k)=−[∇2f(x(k))]−1∇f(x(k));
(4)算法迭代:计算x(k+1)=x(k)+p(k),以x(k+1)作为下一轮迭代点,令k:=k+1,转第二步。
例:用牛顿法求f(x)=x12+25x22的极小点,其中初始点为x(0)=(2,2)T,ϵ=10−6.
解:求梯度和黑塞阵∇f(x)=[2x150x2],∇2f(x)=[20050],于是∇f(x(0))=[4100],且∇2f(x(0))=[20050]是正定矩阵,故[∇2f(x(0))]−1=[2100501],牛顿方向为p(0)=−[∇2f(x(0))]−1∇f(x(0))=[−2−2],于是x(1)=x(0)+p(0)=[22]+[−2−2]=[00],此时相应的梯度为∇f(x(1))=[00],且∣∣∇f(x(1))∣∣<ϵ,故x(1)=[00]是所求的极小点。
牛顿法在极小点附近收敛性很好、速度快,而最速下降法在极小点附近收敛速度很差。牛顿法要求初始点离最优解不远,若初始点选得离最优解太远时,牛顿法并不能保证其收敛,甚至也不是下降方向,因此经常是将牛顿法与最速下降法结合起来使用,前期用最速下降法,当选代到一定程度后改用牛顿法,可得到较好的效果。
修正牛顿法
为了克服牛顿法的缺点,人们保留了从牛顿法中选取牛顿方向作为搜索方向,摒弃其步长恒取1的做法,而用一维搜索确定最优步长来构造算法,这种方法通常称做修正牛顿法。
修正牛顿法的步骤:
(1)选取初始数据:初始点x(0),终止条件ϵ>0,令k:=0;
(2)求梯度向量∇f(x(k)),若∣∣∇f(x(k))∣∣<ϵ,停止计算,此时x(k)为极小点;否则转下一步;
(3)构造牛顿方向:p(k)=−[∇2f(x(k))]−1∇f(x(k));
(4)进行一维搜索:求λk,使f(x(k)+λkp(k))=minλ≥0f(x(k)+λp(k)),令x(x+1)=x(k)+λkp(k),令k:=k+1,转第二步。
例:用修正牛顿法求f(x)=x1−x2+2x12+2x1x2+x22的极小点,初始点x(0)=(0,0)T,ϵ=10−6.
解:求梯度和黑塞阵∇f(x)=[1+4x1+2x2−1+2x1+2x2],∇2f(x)=[4222],在初始点x(0)处有∇f(x(0))=[1−1],∇2f(x(0))=[4222],牛顿方向为p(0)=−[∇2f(x(0))]−1∇f(x(0))=−[21−21−211][1−1]=[−123],然后沿牛顿方向作一维搜索:x(1)=x(0)+λp(0)=[−λ23λ],f(x(0)+λp(0))=45λ2−25λ,令f′(x(0)+λp(0))=0,得最优步长λ0=1,然后计算在新迭代点处的梯度∇f(x(1)):易知∇f(x(1))=[00],∣∣∇f(x(1))∣∣=0<ϵ,即x(1)为所求极小点。
修正牛顿法虽然比牛顿法有了改进,但也有不足之处:一是要计算黑塞阵及其逆矩阵,工作量较大;二是要求选代点x(k)处的黑塞阵正定,可是有些函数未必能满足,因而牛顿方向未必是下降方向,也有一些函数的黑塞阵不可逆,因此不能确定出后继点,这些都是修正牛顿法与牛顿法的局限性。
6. 约束极值问题-最优性条件
求解带有约束的极值问题常用的方法有:
- 将约束问题化为一个或一系列的无约束极值问题;
- 将非线性规划化为近似的线性规划;
- 将复杂的问题变为较简单问题等等。
考虑只含不等式约束条件下求极小值问题的数学模型:minf(x),s.t. gi(x)≥0,i=1,2,⋯,m,其中可行域H={x∣x∈Rn, 且gi(x)≥0,i=1,2,⋯,m}:
定义:设f:Rn→R1,点x∈Rn,向量p∈Rn(p=0),若存在一个数δ>0,使∀λ∈(0,δ),都有f(x+λp)<f(x),则称向量p是f(x)在点x处的下降方向。
定理:设f:Rn→R1在点x处可微,若存在向量p=x−x∈Rn使得∇f(x)Tp<0,则存在数δ>0,使得对每个λ∈(0,δ)都有f(x+λx)<f(x)。若f(x)在点x处可微,当向量p满足∇f(x)Tp<0,则p就是f(x)在点x处的下降方向。
定义:设x∈H,若有gi(x)=0,则称不等式约束gi(x)≥0为点x处起作用的约束,且将下标集I(x)={i∣gi(x)=0,1≤i≤m}称为点x的起作用下标集。若有gi(x)>0,则称不等式约束gi(x)≥0为点x的不起作用约束。(显然等式约束hj(x)=0都是起作用约束)
设p是点x处的一个可行方向,则存在实数δ>0,使得对任意的λ∈[0,δ],都有x+λp∈H,即gi(x+λp)≥0,i=1,2,⋯,m,将gi(x+λp)在x处作泰勒展开,得gi(x+λp)=gi(x)+∇gi(x)Tλp+o(∣∣λp∣∣):
- 对于起作用的约束,因为gi(x)=0,故当λ>0足够小时,∇gi(x)Tp>0⇒gi(x+λp)≥0,i∈I(x);
- 对于不起作用的约束,因为gi(x)>0,故当λ>0足够小时,由于函数gi(x)(i∈/I(x))的连续性,仍有gi(x+λp)≥0,i∈/I(x)。
定义:对于非线性规划问题,如果可行点x处,各起作用约束的梯度向量线性无关,则称x是约束条件的一个正则点。
定义:将既是点x的下降方向又是点x处的可行方向的向量p称为点x处的可行下降方向,即有{∇f(x)Tp<0∇gi(x)Tp>0,i∈I(x)。
容易验证,点x的可行下降方向p与该点处的目标函数负梯度向量之间夹角成锐角,与该点处所有起作用约束的梯度向量之间夹角也都成锐角。
定理:设x∗∈H,f(x)在x∗处可微,若x∗是局部最优解,则x∗点处必不存在可行下降方向。
库恩-塔克条件(Kuhn-Tucker条件)是确定某点为局部最优解的一阶必要条件,只要是最优点就必满足这个条件。但一般来说它不是充分条件,即满足这个条件的点并非最优点,但是对于凸规划,库恩-塔克条件是充要条件。
只含有不等式约束的库恩-塔克条件:设x∗是极小点,那么x∗可能在可行域H的内部,也可能在可行域的边界上:
- 若x∗在H的内部,实际上是个无约束问题,则x∗必满足条件:∇f(x∗)=0;
- 若x∗位于可行域的边界上,分为几种情况来讨论:
- 设x∗位于一个约束条件形成的边界上:不妨设g1(x)≥0是x∗点处起作用约束,即g1(x∗)=0,若x∗是局部最优解,则−∇f(x∗)与∇g1(x∗)必共线且方向相反(反证法,如果不共线如下左图,此时x∗不会是极值点,只有共线才会如右图),也就是说,必存在实数γ1≥0,使得{∇f(x∗)−γ1∇g1(x∗)=0γ1≥0;

- 设x∗同时位于两个约束条件形成的边界上:不妨设g1(x∗)=0、g2(x∗)=0,此时∇f(x∗)必位于∇g1(x∗)与∇g2(x∗)所形成的夹角内(如下图),若x∗是局部最优解,且∇g1(x∗)与∇g2(x∗)线性无关,则∇f(x∗)必可由∇g1(x∗)与∇g2(x∗)正线性表出,即存在γ1≥0、γ2≥0,使{∇f(x∗)−γ1∇g1(x∗)−γ2∇g2(x∗)=0γ1≥0,γ2≥0;

- 依此类推,有{∇f(x∗)−∑i∈I(x∗)γi∇gi(x∗)=0γi≥0,此处要求点x∗处的作用约束梯度向量组线性无关,即x∗同时也是一个正则点,上式也可以写成⎩⎨⎧∇f(x∗)−∑i∈I(x∗)mγi∇gi(x∗)=0γigi(x∗)=0,i=1,2,⋯,mγi≥0,i=1,2,⋯,m。
含等式和不等式约束的库恩-塔克条件:对于问题minf(x),s.t.{gi(x)≥0,i=1,2,⋯,mhj(x)=0,j=1,2,⋯,l,其中hj(x)=0⇔{hj(x)≥0−hj(x)≥0,从而转化为只含不等式约束的问题,其库恩-塔克条件为⎩⎨⎧∇f(x∗)−∑i=1mγi∇gi(x∗)−∑j=1lλj∇hj(x∗)=0γigi(x∗)=0,i=1,2,⋯,mγi≥0,i=1,2,⋯,m。
例:求下列非线性规划问题的K−T点(库恩-塔克点):minf(x)=2x12+2x1x2+x22−10x1−10x2,s.t.{x12+x22≤53x1+x2≤6.
解:将约束改写为gi(x)≥0的形式s.t.{g1(x)=−x12−x22+5≥0g2(x)=−3x1−x2+6≥0,设所求的K−T点为x∗=(x1,x2)T,则∇f(x∗)=[4x1+2x2−102x1+2x2−10],∇g1(x∗)=[−2x1−2x2],∇g2(x∗)=[−31],
写出库恩-塔克条件⎩⎨⎧4x1+2x2−10+2γ1x1+3γ2=02x1+2x2−10+2γ1x2+γ2=0γ1(5−x12−x22)=0γ2(6−3x1−x2)=0γ1≥0γ2≥0:
①如果两个约束均是x∗处不起作用约束,则由{g1(x)=5−x1∗2−x2∗2>0g2(x)=6−3x1∗−x2∗>0⇒{γ1=0γ2=0,代入库恩-塔克条件解出{x1=0x2=5,但该点不满足g1(x∗)≥0,故该点不是可行点;
②若g1(x)≥0是起作用约束,g2(x)≥0是不起作用约束,则由γ2(6−3x1−x2)=0知γ2=0,代入库恩-塔克条件解出x1=1,x2=2,γ1=1或γ1=0,可知x∗=(1,2)T是可行点,且满足库恩-塔克条件,又是一个正则点,故它是一个K−T点;又因为γ1=0也成立,此时和第一张情形一样,该点也不是可行点;
③若g1(x)≥0是不起作用约束,g2(x)≥0是起作用约束,则由γ1(5−3x12−x22)=0知γ1=0,代入库恩-塔克条件解出γ2=0或γ2=−52,当γ2=−52时不满足γ2≥0的条件,当γ2=0时与第一种情形一样;
④若两个约束都是起作用的,此时γ1>0、γ2>0,于是库恩-塔克条件化简为⎩⎨⎧4x1+2x2−10+2γ1x1+3γ2=02x1+2x2−10+2γ1x2+γ2=05−x12−x22=06−3x1−x2=0,但是得出的解不满足γ1≥0,γ2≥0,故舍去;
综上,本题目的库恩-塔克点为x∗=(1,2)T。
本题的目标函数f(x)是凸函数,而g1(x)≥0和g2(x)≥0都是凹函数,故本题是凸规划,而对于凸规划而言,库恩-塔克条件也是充分的,因此x∗=(1,2)T也是本题的全局极小点。
7. 约束极值问题-可行方向法
可行方向法是用一个线性规划来确定搜索方向的方法,我们将约束条件分为线性和非线性的两种情形来讨论。
约束条件为线性
考虑minf(x),s.t. {Ax≥bEx=e,其中f(x)是可微函数,A为m×n矩阵,E为l×n矩阵,x∈Rn,b∈Rm,e∈Rl,即问题有m个线性不等式约束,有l个线性等式约束。
定理:设x是问题的一个可行解,且A1x=b1,A2x>b2,其中A=[A1A2],b=[b1b2](即将x的起作用约束部分的系数矩阵(不一定是方阵)作为A1、右端项作为b1,不起作用约束的系数矩阵作为A2、右端项作为b2),则下列结论成立:
- 非零向量p∈Rn是点x处的可行方向,当且仅当A1p≥0且Ep=0;
- 若向量p在满足可行方向条件的前提下,同时又满足∇f(x)Tp<0,则p是一个可行下降方向。
可行方向法计算步骤(约束条件为线性函数时):
(1)给定初始可行点x(0)∈H,允许误差ϵ1>0、ϵ2>0,令k:=0;
(2)在点x(k)处把A与b分解成A=[A1A2]、b=[b1b2],使得A1x(k)=b1、A2x(k)>b2;
(3)判断x(k)是否为问题可行域的内点;
- 若x(k)是可行域的一个内点(此时问题没有等式约束,即E=0且A1=0),而且∣∣∇f(x(k))∣∣<ϵ1,停止迭代,得到近似极小点x(k);
- 若x(k)是可行域的一个内点,且∣∣∇f(x(k))∣∣≥ϵ1,则取搜索方向p(k)=−∇f(x(k)),然后转第五步,即用目标函数的负梯度方向作搜索方向来求步长,此时类似于无约束问题;
- 若x(k)不是可行域的一个内点(即x(k)在可行域的边界上),则要寻找可行下降方向,转第四步;
(4)求可行下降方向就是要求解问题minz=∇f(x)Tp,s.t. ⎩⎨⎧A1p≥0Ep=0−1≤d1,d2,⋯,dn≤1,设得到的最优解为(p(k),z(k)),若∣z(k)∣=∣∇f(x(k))Tp(k)∣<ϵ2,则停止迭代,否则转第五步;
(5)先计算λ的上限λ={min{(A2p(k))i(b2−A2x(k))i∣(A2p(k))i<0}∞,当A2p(k)≱0,当A2p(k)≥0,再作一维搜索minf(x(k)+λp(k)),s.t. 0≤λ≤λ,求得最优解λk,令x(k+1)=x(k)+λkp(k);
(6)令k:=k+1,返回第二步。
搜索方向和步长计算步骤的推导:
根据定理,要寻找问题可行点x的一个可行下降方向p=(d1,d2,⋯,dn)T,相当于求解如下的一个线性规划问题,问题一:minz=∇f(x)Tp,s.t. ⎩⎨⎧A1p≥0Ep=0−1≤d1,d2,⋯,dn≤1(其中把p的各分量限制在[−1,1]之间是为了讨论方便,我们只关心它的方向而非大小),显然p=0是其中一个可行解,故目标函数的最优值z∗≤0,当z∗<0时,则p为可行下降方向;当z∗=0时,则x为库恩-塔克点。
确定了搜索方向,然后就是步长的确定了。假定x∈H是第k次迭代点的出发点x(k),其可行下降方向为p(k),则后继点x(k+1)为x(k+1)=x(k)+λp(k),为使x(k+1)∈H,且f(xk+1)的值尽可能小,求解一维搜索问题,问题二:min0≤λ≤λf(x(k)+λp(k)),其中λ=max{λ∣x(k)+λp(k)∈H}。
考虑到线性约束Ax≥b,Ex=e,那么求解问题二时首先求解问题三:minf(x(k)+λp(k)),s.t. {A(x(k)+λp(k))≥bE(x(k)+λp(k))=e,因为p(k)是可行方向,则Ep(k)=0,而x(k)是可行点,则Ex(k)=e,说明问题三中的第2个约束必定满足(因为E(x(k)+λp(k))=Ex(k)+λEp(k)=e+0=e),可不再考虑它。
在问题三的基础上,如果在点x(k)处将不等式约束分为起作用约束和不起作用约束,得问题四:minf(x(k)+λp(k)),s.t. {A1x(k)+λA1p(k)≥b1A2x(k)+λA2p(k)≥b2,其中A=(A1,A2)T,b=(b1,b2)T,因为p(k)是可行方向,根据前面定理A1p(k)≥0,再设λ≥0及A1x(k)=b1,故问题四中的第1个条件自然满足,可不再考虑它。
根据问题三、四的讨论知,问题二可以被简化为问题五:minf(x(k)+λp(k)),s.t. A2x(k)+λA2p(k)≥b2,要使步长尽可能大,就是要求λ的上界,λ的上界由下面的公式给出:λ={min{(A2p(k))i(b2−A2x(k))i∣(A2p(k))i<0}∞,当A2p(k)≱0,当A2p(k)≥0,此时问题二再次被化简为问题六:minf(x(k)+λp(k)),s.t. 0≤λ≤λ。
综上,对于约束是线性函数的非线性规划问题,若已知一个迭代点x(k)后,可由求解问题一得到可行下降方向p(k),再由问题六求解该方向上的步长λk。
例:用可行方向法求解minf(x)=(x1−6)2+(x2−2)2,s.t. ⎩⎨⎧x1−2x2≥−4−3x1−2x2≥−12x1,x2≥0,取x(0)=(2,3)T,ϵ1=0.001,ϵ2=0.001.
解:
第一次迭代:
将x(0)=(2,3)T代入约束条件⎩⎨⎧g1(x(0))=2−2×3=−4g2(x(0))=−3×2−2×3=−12g3(x(0))=2>0g4(x(0))=3>0,因此第1、2个约束为点x(0)的起作用约束,故A1=[1−3−2−2]、b1=[−4−12]。
先确定搜索方向,令p=(d1,d2)T,求解线性规划minz=∇f(x(0))Tp,s.t. ⎩⎨⎧A1p≥0Ep=0−1≤d1,d2≤1,因为E=0,∇f(x(0))=(−8,2)T,故线性规划进一步简化为minz=−8d1+2d2,s.t. ⎩⎨⎧d1−2d2≥0−3d1−2d2≥0−1≤d1,d2≤1,用单纯形法解得d1=32、d2=−1,因此x(1)=x(0)+λp(0)=[2+32λ3−λ];
再作一维搜索(A2x(k)+λA2p(k)≥b2)来确定λ的上界λ={min{(A2p(k))i(b2−A2x(k))i∣(A2p(k))i<0}∞,当A2p(k)≱0,当A2p(k)≥0,求出b=b2−A2x(0)=[00]−[1001][23]=[−2−3],p=A2p(0)=[1001]⋅[32−1]=[32−1],故λ=min{−1−3}=3;
为确定步长λ,还要求解下述规划minz=f(x(1))=913(λ−1533)2+13100,s.t. 0≤λ≤λ,故最优解和最优值为λ1=1333、f∗(x(1))=13100,故x(1)=x(0)+λ1p(0)=[1348136],且∇f(x(1))=(−1360,−1340)T,但是∣∣∇f(x(1))∣∣=5.547>0.001=ϵ1,还要进行下一次迭代;
第二次迭代:
对x(1)寻找可行下降方向:⎩⎨⎧g1(x(1))=1348−2×136=1336>−4g2(x(1))=−3×1348−2×136=−12g3(x(1))=1348>0g3(x(1))=136>0,因此第2个约束是点x(1)的起作用约束,故A1=(−3,−2)、b1=−12;
确定搜索方向:令p(1)=(d1,d2)T,求解线性规划minz=∇f(x(1))Tp(1)=−1360d1−1340d2,s.t.⎩⎨⎧−3d1−2d2≥0−1≤d1≤1−1≤d2≤1,解得d1=32、d2=−1,相应的目标函数值z∗=0,因为∣z∗∣=∣∇f(x(1))Tp(1)∣=0<0.001=ϵ2,此时已经没有可行下降方向了,即已经达到极值点(由于z∗=0,该点也是库恩塔克点),故x∗=(1348,136)即为所求。
约束条件非线性
考虑非线性规划minf(x),s.t. gi(x)≥0,i=1,2,⋯,m,其中x∈Rn,f(x),gi(x)(i=1,2,⋯,m)均为可微函数。
定理:设x是问题的一个可行解,I={i∣gi(x)=0}是点x处起作用约束下标集,若f(x),gi(x)(i∈I)均为可微函数,gi(x)(i∈/I)在点x处连续,如果{∇f(x)Tp<0∇gi(x)Tp>0,i∈I(x),则p是可行下降方向。
上述定理中的条件等价于用方程组求向量p和η:⎩⎨⎧∇f(x)Tp≤η−∇gi(x)Tp≤η ,i∈I(x)η<0,满足这一条件的可行下降方向p和数η一般有很多个,我们希望求出能使目标函数值下降最多的方向p。
则求可行下降方向p=(d1,d2,⋯,dn)T的问题可转化为求下列线性规划问题:minη,s.t.⎩⎨⎧∇f(x)Tp≤η−∇gi(x)Tp≤η ,i∈I(x)−1≤d1,d2,⋯,dn≤1。
得到x(k)与可行下降方向p(k)后,沿方向p(k)进行一维搜索确定λk:minf(x(k)+λp(k)),s.t. 0≤λ≤λ,其中λ=sup{λ∣gi(x(k)+λp(k))≥0,i=1,2,⋯,m},
可行方向法计算步骤(约束条件为非线性函数时):
(1)给定初始可行点x(0)∈H,允许误差ϵ1>0、ϵ2>0,令k:=0;
(2)确定x(k)点处的起作用下标集I(x(k))={i∣gi(x(k))=0,1≤i≤m};
(3)判断x(k)是否是问题可行域的内点:
- 若I(x(k))=∅(即不在边界上,是内点),且∣∣∇f(x(k))∣∣<ϵ1,停止迭代,得极小点x(k);
- 若I(x(k))=∅,但∣∣∇f(x(k))∣∣≥ϵ1,则取p(k)=−∇f(x(k))为搜索方向,然后转第五步;
- 若I(x(k))=∅(在边界上,不是内点),转第四步;
(4)确定搜索方向,即求解线性规划问题minη,s.t.⎩⎨⎧∇f(x)Tp≤η−∇gi(x)Tp≤η ,i∈I(x)−1≤d1,d2,⋯,dn≤1,设得到的最优解为(p(k),η(k)),若∣η(k)∣<ϵ2,停止迭代,得到极小点x(k),否则以p(k)为搜索方向转第五步;
(5)先计算λ的上限λ=sup{λ∣gi(x(k)+λp(k))≥0,i=1,2,⋯,m},然后作一维搜索minf(x(k)+λp(k)),s.t. 0≤λ≤λ,求得最优解λk,令x(k+1)=x(k)+λkp(k);
(6)令k:=k+1,返回第二步。