无约束非线性规划求解方法及其实现

无约束非线性规划求解方法及其实现

作者：杨玲指导老师：陈素根

摘要:

非线性规划是具有非线性约束条件或目标函数的数学规划，是运筹学的一个重要分支。非线性规划属于最优化方法的一种，是线性规划的延伸。非线性规划研究一个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维向量空间：n

错误！未找到引用源。,则称

x R

这样的非线性规划问题为无约束问题：错误！未找到引用源。或错误！未找到引用源。 。

一般我们研究的无约束非线性规划问题大都可以归结为求无约束最优化问题。

2.2 介绍无约束非线性规划的几种主要的求解方法及其实现 求解无约束非线性规划问题就是求解无约束非线性规划最优化

的问题，可以表述为()min ,n

f x x R ∈。它的求解方法有许多种，

大体上可以概括为两大类，一是直接法，二是解析法。解析法又被称为代数法，值得是通过计算

()f x 的一阶，二阶偏导数及其

函数的解析性质来实现极值的求解方法。相应的，不必计算

()f x 的一阶、二阶偏导数及其函数的解析性质，仅用到函数值

来实现近似值的求解方法叫直接法。

1. 先介绍直接法中的一维搜索方法，包括Fibonacci 法和0.618法。

一维搜索方法就是在用迭代法沿某一已知方向求目标函数极小点的方法，常用的由斐波那契法和黄金分割法。

考虑一维极小值问题()min a t b f t ≤≤，若()f t 是[],a b 区间上的下单峰函数，我们将通过不断的缩短[],a b 的长度，来探索()min a t b f t ≤≤的近似最优解。在[],a b 中任意取两个关于[],a b 是对称的点1t 和1t (不妨设,2

1t t <并称它们为搜索点)，计算()1f t 与()2f t 并

比较它们的大小。对于单峰函数，若

()()12f t f t <,则必有

[]

1*,t a t ∈，

因而[]1,a t 是缩短了的单峰区间，若()()21f t f t <,

则有

[]

2*,t t b ∈，故[]2

,t b 是缩短了的单峰区间，若

()()21f t f t =，则[]1,a t 和[]2,t b 都是缩短了的单峰。因而通

过两个搜索点处目标函数值大小的比较，总可以获得缩短了的单峰区间。对于新的单峰区间重复上述做法，又可以获得更短的单峰区间。如此下去，在单峰区间缩短到充分小时，可以取最后的

搜索点作为()min a t b

f t ≤≤最优解的近似值，下面介绍斐波那契法来选取搜索点，使给定的单峰区间的长度能尽快缩短。 Fibonacci 法： 若数列

{}

n F 满足关系：011F F ==，21n n n F F F --=+，2,3,

n =，则称n F 为Fibonacci 数列，n F 称为第n 个

Fibonacci 数，称相邻两个Fibonacci 数之比1

n n

F F -为Fibonacci 分

数。当用斐波那契法以n 个探索点来缩短某一区间时，区间长度

的第一次缩短率为1

n n

F F -，其后各次分别为

231

122

,,,n n n n F F F F F F ----，由此，若1t 和2t ，

()21t t <单峰区间[],a b 中的第1个和第2

个探索点的话，则应有比例关系1

1n n F t a b a F --=-，

22a n n t F b a F --=-，从

而()11n n F t a b a F -=+-，

()2

2n n

F t a b a F -=+-，它们关于[]

,a b 是对称的点。

如果要求经过一系列探索点搜索之后，使最后的探索点和最优解之间的距离不超过精度0δ>，这也要求最后区间的长度不超过

δ

，即

n

b a

F δ-≤。由此，按照预先给定的精度δ，确定使n b a F δ-≤成立的最小整数n 作为搜索次数，直到进行第n 次探索点为止。

用上述不断缩短函数()f x 单峰区间的方法来求()min a t b

f t ≤≤的近似解是Kiefer(1953年)提出的，叫Fibonacci 法。具体步骤如下：

1选取初始数据，确定单峰区间[]00,a b ,给出搜索精度0δ>，由

n

b a

F δ-≤确定搜索次数n ； 21k =，0a a =，0b b =，计算最初两个搜索点，按

()11n n F t a b a F -=+-，

()2

2n n F t a b a F -=+-计算1t ，2t ；

3while 1k n <-

()11f f t =,

()22f f t =

if 12f f <

2a t =；21t t =；()11n k

n k

F t a b a F ---=+-；

else

1b t =；12t t =；()12n k

n k

F t b b a F ---=+

-；

end

1k k =+

end

4当进行至1k n =-时，()121

2t t a b ==+，此时无法比较

()1f t 与()2f t 的大小来确定最终区间，为此，取

()()211212t a b t a b a ε?

=+???

???=++- ?????

，其中ε为任意小点的数，在1t 和2t 这两点中，以函数值较小者为近似极小值，相应的函数值为近似极小值，并得最终区间

[]1,a t 或[]2,t b 。

由上述分析可知，斐波那契法使用对称搜索方法，逐步缩短所考察的区间，它能以尽量少的函数求值次数达到预定的某一缩短率。

2. 下面介绍解析法中的最速下降法、牛顿法、共轭梯度法和变尺度法。

（1)最速下降法：

对基本迭代格式

1k k k

k x x t p

+=+，我们一般考虑从点

k

x

出发沿哪一个方向

k

p

，使目标函数下

f 降得最快。

由微积分的相关知识我们可以知道，点

k

x

的负梯度方向

()k k

p f x =-?是从点

k

x

出发使

f

下降最快的方向。为此，

负梯度方向

()k

f x -?为

f

在点

k

x

处的最速下降

方向，按基本迭代式1

k k

k

k x x t p

+=+每一轮从点

k

x

出发沿

最速下降方向

()k

f x -?作一维搜索，来建立求解无约束

极值问题的方法称之为最速下降法。该方法的特点是每轮的搜索方向都是目标函数在当前点下降最快的方向。同时，

()0k f x ?=或

()k

f x ε?≤作为停止条件。其具体的步骤为：(a).选取初始数据。选取初始点0

x

，给定终止误差，

令k:=0. (b).求梯度向量。计算()k

f x ?，若()k

f x ε?≤，停

止迭代，输出

k

x

。否则进行(c)。 (c). 构造负梯度方向。取()k

k

p f x =-?。

(d). 进行一维搜索。求

k

t ，使得

()m i n ()k k k k

k t f x t p f x t p ≥+=+

，令

1k k k

k x x t p +=+，

=1k k +：进行(b).

例题：用最速下降法求解无约束非线性规划问题，

()2

21

2

min 25f x x x

=+，其中

()

12,T

x x x =，要求选取初始点

()

02,2T

x =，终止误差

6

10ε-=。 解：1）

()()122,50T

f x x x ?=,编写M 文件deta f.m 如下：

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

(2). 牛顿法：

牛顿法是求无约束最优解的一种古典解析算法，其基本思想是在寻找收敛速度最快的无约束最优化方法中，在每次迭代时，用适当的二次函数去近似目标函数

f

，并用迭代点指向近似二次

函数极小点的方向来构造搜索方向，然后精确地求近似二次函数的极小点，以这个极小点作为f

的极小点的近似值。 牛顿法的迭代步骤： 1）给定初始点和收敛精度；

2）计算

()k f x ?，

()

k H x ，

()1

k H x -????

3）求

()()

1

1k k k k k x x K H x f x -+=-?????；

4）检查收敛精度，若

1k k x x ε

+-<，则

1

*k x x +=，停止；否则

1

k k =+，返回2）继续。

牛顿法的优点是每次迭代都在牛顿方向进行一维搜索，避免了迭代后函数值变大的现象，从而保持了牛顿法的二次收敛性，而对初始点的选择没有严格要求。但是牛顿法也有缺点，它对目标函数要求严格，函数必须具有连续的一、二阶导数；海赛矩阵必须正定且非奇异。还有牛顿法的计算复杂，存储量大。

(3). 共轭梯度法：

共轭梯度法最早是由计算数学家Hestenes 和几何学家Stiefel

在20世纪50年代初为求解线性方程组,n

Ax b x R =∈而各

自独立提出的。在

A

为对称正定阵时，方程组

,n

Ax b x R =∈等价于最优化问题1min 2

n T T x R x Ax b x ∈，

由此，Hestenes 和Stief 的方法也可以看做是求二次函数极小值得共轭梯度法。在1964年Fetcher 和Reeves 将这种方法推广到了非线性优化，得到了求一般函数极小值的共轭梯度法。对于无

约束最优化问题

()min n

x R

f x ∈，其中

:n

f R R →连续可微有下

界，共轭梯度法是解决这类问题中的最有效的数值方法之一。特别是在大规模问题上，共轭梯度法因为算法简便、所需存储量小、收敛速度快等特性而在许多工程科学领域采用。

对于无约束优化问题，给出一个初始值0x ，算法迭代产生点

列{}12,,x x 。假设某一k x 是无约束优化问题的解，或者该

点列收敛于最优解。在第

1k +次迭代中，当前迭代点为

k x ，

产生下一个迭代点1k k k k x x a d +==，其中n k

d R ∈是搜索方向，

0k a >是步长因子，它满足某线搜索终止条件。显然，每步迭

代主要由两部分组成：一是搜索方向k d ；另一是步长因子k a 。

求解无约束优化问题的共轭梯度法是从求解线性方程组的线性共轭梯度法推广而来的，其搜索方向是负梯度方向与上一次迭代的搜索方向的线性组合，它表示为

0011,k k k k d g d g d β--=-=-+，关于参数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. 给定初始点()

0x ，初始对称正定矩阵0H ，

()

()0

0g g x =及

精度

0ε>； b. 计算搜索方向

()k k k p H g =-； c. 作直线搜索

()

()()

(

)1,k k k x

F x p

+=，计算()

()

11

k k f

f x

++=，

()

()

11k k g g x

++=，

()

()

1k k k S x

x

+=-，

1k k k y g g +=-； d. 判断终止准则是否满足； e. 令

1K k k H H E +=+置1k k =+，转步骤(b).(不同的拟牛顿

法对应不同的

k

E )

DFP 变尺度法的算法原理：

DFP

算法中校正公式为

1T T k k k k k k k k T T

k k

k k k S S H y y H

H H S y y H y +=+-

，

为了保证

k H 的正定性，在下面的算法迭代一定的次数后，

重置初始点和迭代矩阵再进行迭代。 DFP 变尺度法的算法步骤： 1） 给定初始点()

0x

，初始矩阵0n H I =及精度0ε

>；

2） 若

()k

f x ε?≤，停止，极小点为()0x ；否则转步骤3）；

3） 取()

()

000()p H f x =-?,令0k =；

4） 用

一

维

搜

索

法

求

k

t ，使得

()()

(

)()

()

()

min k k k k k k t f x t p

f x

tp

α≥+=+，令

()

()

()

1k k k x

x tp

+=+，转步骤5）；

5）

()

1()k f x

ε+?≤，停止，极小值点为()

1k x

+；否则转步

骤6）；

6） 若1k n +=，令()

()

0n x x

=，转步骤3);否则转步骤7);

7） 令

()

()

()()()

()

()

()

()()

()

()

11111()()

T

k k k k k k T

k k k k x x

x

x

H H x

x

f x

f x +++++--=+

-?-?

()

()

(

)()

()

(

)

()

()

()()

()

(

)

1111()()()()

()()()()

T

k k k k k T

k k k k k H f x

f x f x

f x f x

f x H f x

f x ++++?-??-?-

?-??-?，取

()

()

11()k k k p H f x

++=-?，置1k k =+，步骤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<=eps

x=x0;

break;

end

y=x0+l*p;

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<=eps

x=x1;

break;

end

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;