无约束非线性规划求解方法及其实现
无约束非线性规划求解方法及其实现
作者：杨玲 指导老师：陈素根 摘要:
非线性规划是具有非线性约束条件或目标函数的数学规划，是运筹学的一个重要分支。非线性规划属于最优化方法的一种，是线性规划的延伸。非线性规划研究一个n元实函数在一组灯饰或不等式的约束条件下的极值问题，且目标函数和约束条件至少有一个是未知量的非线性函数。目标函数和约束条件都是线性函数的情形则属于线性规划。非线性规划是20世纪50年代才形成的一门新兴学科。1951年H.W库恩和A.W塔克发表的关于最优性条件的论文是非线性规划正是诞生的一个重要标志。在50年代还得出了可分离规划和二次规划的n种解法，它们大都是以G.B.丹齐克提出的解线性规划的单纯形法为基础的。50年代末到60年代末出现了许多解线性规划问题的有效的算法，70年代又得到进一步的发展。非线性规划在工程，管理，经济，科研，军事等发面都有广泛的应用，为最优设计提供了有力的工具。20世纪80年代以来，随着计算机技术的快速发展，非线性规划在信赖域法、稀疏牛顿法、并行计算、内点法和有限存储法等领域取得了丰硕的成果，无约束非线性规划问题是非线性规划的一个重要内容，很多学者对非线性规划问题进行了深入且系统的研究，研究成果丰硕。
关键词 最优化 共轭梯度法 非线性 无约束
1 引言
1.1 无约束非线性规划问题是最基本的非线性规划问题，在1959~1963年幼三位数学家共同研究成功求解无约束问题的DFP变尺度法，该算法的研究成功是无约束优化算法的一个大飞跃，引起了一系列的理论工作，并陆续出现了许多新的算法。20世纪80年代以来，随着计算机技术的快速发展，非线性规划在信赖域法、稀疏牛顿法、并行计算、内点法和有限存储法等领域取得了丰硕的成果。无约束非线性规划问题是非线性规划的一个重要内容，很多学者对非线性规划问题进行了深入且系统的研究，研究成果丰硕。
1.2 本文主要研究无约束非线性规划问题，将文章分成四个部分，首先会具体介绍无约束非线性规划的相关概念，并在此基础上研究非线性规划的相关理论与基本算法问题，接着详细介绍无约束非线性规划的几种主要的求解方法，最后举例说明他在实际生活中的应用，并编程实现它。 2 正文
2.1主要介绍无约束非线性规划的相关概念
一个非线性规划问题的自变量x没有任何约束，或说可行域即是整个n维向量空间：x Rn错误！未找到引用源。,则称
这样的非线性规划问题为无约束问题：错误！未找到引用源。或错误！未找到引用源。 。
一般我们研究的无约束非线性规划问题大都可以归结为求无约束最优化问题。
2.2 介绍无约束非线性规划的几种主要的求解方法及其实现 求解无约束非线性规划问题就是求解无约束非线性规划最优化
n
minfx,x∈R()的问题，可以表述为。它的求解方法有许多种，
大体上可以概括为两大类，一是直接法，二是解析法。解析法又被称为代数法，值得是通过计算
f(x)的一阶，二阶偏导数及其
函数的解析性质来实现极值的求解方法。相应的，不必计算
f(x)的一阶、二阶偏导数及其函数的解析性质，仅用到函数值
来实现近似值的求解方法叫直接法。
1. 先介绍直接法中的一维搜索方法，包括Fibonacci法和0.618法。
一维搜索方法就是在用迭代法沿某一已知方向求目标函数极小点的方法，常用的由斐波那契法和黄金分割法。
f(t)，若f(t)是[a,b]区间上的下单考虑一维极小值问题mina≤t≤bf(t)峰函数，我们将通过不断的缩短[a,b]的长度，来探索mina≤t≤b的近似最优解。在[a,b]中任意取两个关于[a,b]是对称的点t1和t1(不妨设,t2
比较它们的大小。对于单峰函数，若
f(t1)
t*∈[a,t1]
则有
因而，
[a,t1]是缩短了的单峰区间，若f(t2)
[t,b]是缩短了的单峰区间，若
，故2
t*∈[t2,b]
f(t2)=f(t1)，则[a,t1]和[t2,b]都是缩短了的单峰。因而通
过两个搜索点处目标函数值大小的比较，总可以获得缩短了的单峰区间。对于新的单峰区间重复上述做法，又可以获得更短的单峰区间。如此下去，在单峰区间缩短到充分小时，可以取最后的
f(t)最优解的近似值，下面介绍斐波那契法搜索点作为mina≤t≤b
来选取搜索点，使给定的单峰区间的长度能尽快缩短。 Fibonacci法： 若数列
{Fn}
满足关系：F0=F1=1，Fn=Fn-2+Fn-1，，则称Fn为Fibonacci数列，Fn称为第n个
n=2,3,
Fn-1
Fibonacci数，称相邻两个Fibonacci数之比F为Fibonacci分
n
数。当用斐波那契法以n个探索点来缩短某一区间时，区间长度
Fn-2Fn-3Fn-1
,,的第一次缩短率为，其后各次分别为
Fn-1Fn-2Fn
由此，若t1和t2，
F1
,，F2
(t2
t2-aFn-2t1-aFn-1
=个探索点的话，则应有比例关系b-a=F，从b-aFnn，
而t1=a+的点。
Fn-1Fn-2
b-a()t2=a+(b-a)，它们关于a,b是对称FnFn
，
[]
如果要求经过一系列探索点搜索之后，使最后的探索点和最优解之间的距离不超过精度δ>0，这也要求最后区间的长度不超过
δ
b-ab-a
≤δ≤δ。由此，按照预先给定的精度δ，确定使，即
FnFn
成立的最小整数n作为搜索次数，直到进行第n次探索点为止。
f(t)的近似用上述不断缩短函数f(x)单峰区间的方法来求mina≤t≤b
解是Kiefer(1953年)提出的，叫Fibonacci法。具体步骤如下：
1选取初始数据，确定单峰区间[a0,b0],给出搜索精度δ>0，由
b-a
≤δ确定搜索次数n； Fn
2k=1，a=a0，b=b0，计算最初两个搜索点，按
Fn-1Fn-2
t1=a+b-a()t2=a+(b-a)计算t1t2；
，FnFn，
3while k
f1=f(t1)f2=f(t2)
,
if f1
Fn-1-k
t=a+(b-a) a=t2；t2=t1；1Fn-k
；
else
Fn-1-k
(b-a) b=t1；t1=t2；t2=b+
Fn-k
；
end
k=k+1
end
1
4当进行至k=n-1时，t1=t2=2(a+b)，此时无法比较
f(t1)与f(t2)的大小来确定最终区间，为此，取
1⎧
t2=(a+b)⎪2⎪⎨
1⎫ε为任意小点的数，在t1和t2，其中⎪t1=a+⎛+ε⎪(b-a) ⎪⎝2⎭⎩
这两点中，以函数值较小者为近似极小值，相应的函数值为近似极小值，并得最终区间
[a,t1]或[t2,b]。
由上述分析可知，斐波那契法使用对称搜索方法，逐步缩短所考察的区间，它能以尽量少的函数求值次数达到预定的某一缩短率。
2. 下面介绍解析法中的最速下降法、牛顿法、共轭梯度法和变尺度法。
（1)最速下降法：
对基本迭代格式
x=x+tkp
k+1kk
，我们一般考虑从点
x
k
出发沿哪一个方向
p
k
f，使目标函数下
x
k
降得最快。
由微积分的相关知识我们可以知道，点
的负梯度方向
p=-∇f(x)是从点
负梯度方向
kk
x
k+1
k
出发使
k
f
下降最快的方向。为此，
-∇f(x)为
k
f
k
在点
x
k
处的最速下降
方向，按基本迭代式x最速下降方向
=x+tkp
k
每一轮从点
x
k
出发沿
-∇f(x)作一维搜索，来建立求解无约束
极值问题的方法称之为最速下降法。该方法的特点是每轮的搜索方向都是目标函数在当前点下降最快的方向。同时，
k
∇f(xk)=0∇f(x)≤ε作为停止条件。其具体的步或
骤为：(a).选取初始数据。选取初始点令k:=0. (b).求梯度向量。计算∇f(x止迭代，输出
k
x
k
，给定终止误差，
k
)，若∇f(x)≤ε，停
x
k
k
。否则进行(c)。 (c). 构造负梯度方向。取(d). 进行一维搜索。求
k
p=-∇f(x)。
t≥0
tk
，使得
，
kkf(x+ktkp)=minf(x+
k
tp)xk+1=xk+ktp，令
k：=k+1进行(b).
例题：用最速下降法求解无约束非线性规划问题，
minf(x)=x+25x
2
1
22
，其中
x=(x1,x2)
T
，要求选取初始点
x0=(2,2)
T
-6
ε=10。 ，终止误差
∇fx=2x,50x()(12)解：1）
function.[f,df]=deta f(x); f=x(1)^2+25*x(2)^2; df(1)=2*x(1); df(2)=50*x(2)^2; 2)编写M文件zuisu.m clc x=[2,2]; [f0,g]=deta f(x);
while norm(g)>0.000001 p=-g/norm(g); t=1.0;f=deta f(x+t*p); end x=x+t*p [f0,g]=deta f(x) end
T
,编写M文件deta f.m如下：
(2). 牛顿法：
牛顿法是求无约束最优解的一种古典解析算法，其基本思想是在寻找收敛速度最快的无约束最优化方法中，在每次迭代时，用适当的二次函数去近似目标函数
f
，并用迭代点指向近似二次
函数极小点的方向来构造搜索方向，然后精确地求近似二次函数的极小点，以这个极小点作为的极小点的近似值。 牛顿法的迭代步骤： 1）给定初始点和收敛精度；
f
∇fx(k)2）计算
3）求
，
H(xk)
-1
H⎡⎣，
(xk)⎤⎦
-1
xk+1=xk-Kk⎡⎣H(xk)⎤⎦∇f(xk)
；
4）检查收敛精度，若
xk+1-xk
，则
x*=xk+1
，停止；否则
k=k+1
，返回2）继续。
牛顿法的优点是每次迭代都在牛顿方向进行一维搜索，避免了迭代后函数值变大的现象，从而保持了牛顿法的二次收敛性，而对初始点的选择没有严格要求。但是牛顿法也有缺点，它对目标函数要求严格，函数必须具有连续的一、二阶导数；海赛矩阵必须正定且非奇异。还有牛顿法的计算复杂，存储量大。
(3). 共轭梯度法：
共轭梯度法最早是由计算数学家Hestenes和几何学家Stiefel
n
Ax=b,x∈R在20世纪50年代初为求解线性方程组而各
自独立提出的。在
A
为对称正定阵时，方程组
1TTminxAxbxn
nAx=b,x∈R等价于最优化问题x∈R2
，
由此，Hestenes和Stief的方法也可以看做是求二次函数极小值得共轭梯度法。在1964年Fetcher和Reeves将这种方法推广到了非线性优化，得到了求一般函数极小值的共轭梯度法。对于无
约束最优化问题
minf(x)n
x∈R
，其中
f:R→R连续可微有下
n
界，共轭梯度法是解决这类问题中的最有效的数值方法之一。特别是在大规模问题上，共轭梯度法因为算法简便、所需存储量小、收敛速度快等特性而在许多工程科学领域采用。
x 对于无约束优化问题，给出一个初始值0，算法迭代产生点
xx1,x2,}{列。假设某一k是无约束优化问题的解，或者该
k+1次迭代中，当前迭代点为
点列收敛于最优解。在第
xk，
nx=x=add∈Rkkk，其中k
产生下一个迭代点k+1是搜索方向，
ak>0是步长因子，它满足某线搜索终止条件。显然，每步迭
dak代主要由两部分组成：一是搜索方向；另一是步长因子k。
求解无约束优化问题的共轭梯度法是从求解线性方程组的线性共轭梯度法推广而来的，其搜索方向是负梯度方向与上一次迭代的搜索方向的线性组合，它表示为
d0=-g0,dk=-gk+βk-1dk-1，关于参数βk
的不同取
法对应于不同的共轭梯度法，著名的有HS方法,FR方法,PRP方法,CD方法,LS方法，DY方法。其中FR方法、DY方法和CD方法具有很好的收敛性质，但数值表现结果却差强人意，而PRP方法、HS方法和LS方法对一般函数不具备收敛性，但当收敛时，往往数值表现却很好。因此近年来研究出了混合共轭梯度法，许多学者作出了尝试。焦宝聪、陈兰平和李娟对求解无约束优化问题提出一类三项混合共轭梯度算法，新算法将HS算法与DY方法相结合，并在不需给定下降条件的情况下，证明了算法在Wolfe线搜索原则下的收敛性，数值试验亦显示出这种混合共轭梯度算法较之HS和PRP的优势。焦宝聪、陈兰平和潘翠英Ⅲo结合FR算法和DY算法，给出了一类新的混合共轭梯度算法，并结合Goldstein线搜索，在较弱的条件下证明了算法的收敛性。
共轭梯度法是一种很有效求解无约束优化的方法，共轭梯度法是根据共轭方向去搜索，可以由较快的收敛速度找到最优解求得极小点。
(4).变尺度法：
变尺度法也称为拟牛顿法，它是基于牛顿法的思想而又作出了重大改进的一种方法，是由Davidon提出，经过Fletcher和Powell加以发展和完善的，称为DFP变尺度法。
拟牛顿法的一般步骤为：
a. 给定初始点x ，初始对称正定矩阵H0，
(0)
g0=gx()
()及
ε>0； 精度
(k)p=-Hkgk； b. 计算搜索方向
c. 作直线搜索
x
(k+1)
=Fx,p
(k+1)
(
(k)(k)
)，计算f
k+1
=fx
()
(k+1)
，
gk+1=gx
()
(k+1)
，
Sk=x
-x
(k)
y=g-gkk+1k； ，
d. 判断终止准则是否满足； e. 令
HK+1=Hk+Ek置k=k+1，转步骤(b).(不同的拟牛顿
Ek
)
法对应不同的
DFP变尺度法的算法原理：
DFP
SSHyyH
Hk+1=Hk+-
SyyHy算法中校正公式为
TkkTkk
TkkkkT
kkk，
为了保证
Hk的正定性，在下面的算法迭代一定的次数后，
重置初始点和迭代矩阵再进行迭代。 DFP变尺度法的算法步骤： 1） 给定初始点x
k
(0)
，初始矩阵H0=In及精度ε
>0；
2） 若
∇f(x)≤ε，停止，极小点为x(0)；否则转步骤3）；
(0)
3） 取p4） 用
=-H0∇f(x),令k=0；
维
搜
索
法
求
(0)
一
tk
，使得
tkfx+tkp
x
(k+1)
(k)
(
(k)(k)
)=minf(x
α≥0
(k)
(k)
+tp
(k)
)
，令
=x+tp
∇f(x
(k+1)
，转步骤5）；
5）
)≤ε，停止，极小值点为x
(0)
(k+1)
；否则转步
骤6）；
6） 若k+1=n，令x7） 令
=x
(n)
，转步骤3);否则转步骤7);
Hk+1=Hk+
(x(x
(k+1)
(k+1)
-x
T
(k)
)(x
(k+1)(k+1)
-x
(k)
)
T
-x
(k)
)(∇f(x
)-∇f(x)
(k)
)
-
Hk∇f(x
(
(k+1)
)-∇f(x)∇f(x
(k)
(k)
T
)(
(k+1)
)-∇f(x)
(k)
(k)
)
T
(∇f(x
(k)
(k+1)
)-∇f(x)Hk∇f(x
(k+1)
)(
(k+1)
)-∇f(x)
)
，取
p=-Hk+1∇f(x
)，置k=k+1，步骤4)。
DFP变尺度法的算法MATLAB程序： 调用格式：[x,min f]=min DFP(f,x0,var,eps) 其中， f:目标函数 x0:初始点 var自变量向量 eps精度
x：目标函数取最小值时的自变量值 min f :目标函数的最小值
DFP的MATLAB程序代码如下： fuction [x,min f]=minDFP(f,x0,var,eps) %目标函数：f; %初始点：x0; %自变量向量：var;
%目标函数取最小值时的自变量值:x; %目标函数的最小值:min f; format long; If nargin==3
eps=1.0e-6; end
x0=transpose(x0); n=length(var); syms l; H=eye(n,n);
grad f=jacobian(f,var); v0=Funval(gradf,var,x0); p=-H*transpose(v0); k=0; while 1
v=Funval(gradf,var,x0); tol =norm(v); if tol
yf=Funval(f,var,y); [a,b]=minJT(yf,0,0.1); xm=minHJ(yf,a,b); x1=x0+xm*p;
vk=Funval(gradf,var,x1); tol=norm(vk); if tol
if k+1==n x0=x1; continue; else
dx=x1-x0; dgf=vk-v;
dgf=transpose(dgf); dxT=transpose(dx); dgfT=transpose(dgf); mdx=dx*dxT; mdgf=dgf*dgfT; fz=H*(dgf*(dgfT*H));
H=H+mdx/(dxT*dgf)-inv(dgfT*(H*dgf))*fz; p=-H*transpose(vk); k=k+1; x0=x1;
end end
minf=Funval(f,var,x); format short;