非线性规划的算法研究

论文 关键词：非线性规划 最优决策 初值依赖 matlab

论文摘要：本课题主要研究非线性规划的算法。非线性规划在军事， 经济 ，管理，生产过程自动化，工程设计和产品优化设计等方面都有着重要的应用。但非线性规划的研究目前还不成熟，有许多问题需要进一步完善。非线性规划不像线性规划有统一的算法，对于不同的问题需要用不同的算法处理，现阶段各种算法都有一定的局限性，只有对各种算法加以改正，才能有效地解决人们在日常的生产、生活中遇到的优化问题，做出最优决策。 本文主要是对现有的各种算法加以测试，指出各种算法的优缺点，寻找一种不受初值依赖，收敛更快的最优算法。首先介绍了非线性规划研究的背景和国内外研究状况，然后论述了方案的选取过程，重点描实验过程，主要是对各种非线性最优 计算 方法用matlab 软件编程，给出一个在工程中具有代表性的最优函数实例，经过大量的测试，并给出了结果分析。最后给出了整个实验的 总结 和由此对未来的展望。

1 选题背景

1.1 课题背景

1.2 课题研究的目的和意义 一般来说，解非线性规划问题要比求解线性规划问题困难得多，而且也不像线性规划那样有统一的数学模型及如单纯形法这一通用解法。非线性规划的各种算法大都有自己特定的适用范围。都有一定的局限性，到目前为止还没有适合于各种非线性规划问题的一般算法。这正是需要人们进一步研究的课题。非线性规划在工程、管理、经济、科研、军事等方面都有广泛的应用，为最优设计提供了有力的工具。例如：如何在现有人力、物力、财力条件下合理安排产品生产，以取得最高的利润；如何设计某种产品，在满足规格、性能要求的前提下，达到最低的成本；如何确定一个自动控制系统的某些参数，使系统的工作状态最佳；如何分配一个动力系统中各电站的负荷，在保证一定指标要求的前提下，使总耗费最小；如何安排库存储量，既能保证供应，又使储存费用最低；如何组织货源，既能满足顾客需要，又使资金周转最快等。对于静态的最优化问题，当目标函数或约束条件出现未知量的非线性函数，且不便于线性化，或勉强线性化后会招致较大误差时，就可应用非线性规划的方法去处理。

随着社会的发展，实际问题越来越复杂，例如全局最优化问题。经典算法一般都用得局部信息，如单个初始点及所在点的导数等，这使得经典算法无法避免局部极小问题。全局最优化是np-hard 问题, 所以原有的经典算法不再使用，必须对其进行改进，或将其与启发式算法结合。启发式算法是受大 自然 的启发，人们从大自然的运行 规律 中找到了许多解决实际问题的方法。启发式算法的计算量都比较大，所以启发式算法伴随着计算机技术的发展，取得了巨大的成就。

40年代：由于实际需要，人们已经提出了一些解决实际问题快速有效的启发式算法。 50年代：启发式算法的研究逐步繁荣起来。随后，人们将启发式算法的思想和人工智能领域中的各种有关问题的求解的收缩方法相结合，提出了许多启发式的搜索算法。其中贪婪算法和局部搜索等到人们的关注。

60年代: 随着人们对数学模型和优化算法的研究越来越重视，发现以前提出的启发式算法速度很快，但是解得质量不能保证。虽然对优化算法的研究取得了很大的进展，但是较大规模的问题仍然无能为力（计算量还是太大）。

70年代：计算复杂性理论的提出。np 完全理论告诉我们，许多实际问题不可能在合理的时间范围内找到全局最优解。发现贪婪算法和局部搜索算法速度快，但解不好的原因主要是他们只是在局部的区域内找解，得到的解不能保证全局最优性。由此必须引入新的搜索机制和策略，才能有效地解决这些困难问题。

优胜劣汰是大自然的普遍规律，它主要通过选择和变异来实现。选择是优化的基本思想，变异（多样化）是随机搜索或非确定搜索的基本思想。“优胜劣汰”是算法搜索的核心，根据

“优胜劣汰”策略的不同，可以获得不同的超启发式算法。超启发式算法的主要思想来自于人类经过长期对物理、生物、社会的自然现象仔细的观察和实践，以及对这些自然现象的深刻理解，逐步向大自然学习，模仿其中的自然现象的运行机制而得到的。

遗传算法：是根据生物演化，模拟演化过程中基因染色体的选择、交叉和变异得到的算法。在进化过程中，较好的个体有较大的生存几率。

模拟退火：是模拟统计物理中固体物质的结晶过程。在退火的过程中，如果搜索到好的解接受；否则，以一定的概率接受不好的解（即实现多样化或变异的思想），达到跳出局部最优解得目的。

神经网络：模拟大脑神经处理的过程，通过各个神经元的竞争和协作，实现选择和变异的过程。

禁忌搜索：模拟人的经验，通过禁忌表记忆最近搜索过程中的 历史 信息，禁忌某些解，以避免走回头路，达到跳出局部最优解的目的。

蚂蚁算法：模拟蚂蚁的行为，拟人拟物，向蚂蚁的协作方式学习。

虽然人们研究对启发式算法的研究将近５０年，但它还有很多不足：

1）启发式算法目前缺乏统一、完整的理论体系。

2）由于np 理论，各种启发式算法都不可避免的遭遇到局部最优的问题，如何判断已经找到的最优值是全局最优值。

3）各种启发式算法都有各自优点，如何完美结合。

4）启发式算法中的参数对算法的效果起着至关重要的作用，如何有效设置参数。

5）启发算法缺乏有效的迭代停止条件。

6）启发式算法收敛速度的研究等。

2 方案论证

2.1 设计原理

对实际规划问题作定量分析，必须建立数学模型。建立数学模型首先要选定适当的目标变量和决策变量，并建立起目标变量与决策变量之间的函数关系, 称之为目标函数。然后将各种限制条件加以抽象, 得出决策变量应满足的一些等式或不等式，称之为约束条件。一般非线性规划的数学模型可表示为：

有约束问题与无约束问题是非线性规划的两大类问题，它们在处理方法上有明显的不同。实际问题中，大多数都是有约束条件的问题。求解带有约束条件的问题比起无约束问题要困难得多，也复杂得多。在每次迭代时，不仅要使目标函数值有所下降，而且要使迭代点都落在可行域内（个别算法除外）。求解带有约束的极值问题常用方法是：将约束问题化为一个或一系列的无约束极值问题；将非线性规划化为近似的线性规划；将复杂问题变为较简单问题等等。

关于求解非线性规划的迭代法有二分法、简单迭代法、牛顿迭代法与拟牛顿迭代法、同论延拓法、单纯形法等。以上方法都有一定的局限性。目前需寻找一种初值不受限制，以较快的是速度收敛到最优解的迭代法。

2.2方案选择

求解非线性规划需要编写各种算法程序，可供选择的高级语言有很多种，如fortran 和c 在内的多种高级语言均可以实现这些非线性最优算法，这里选择用matlab 语言实现，为什要选择matlab ？在通过一定的比较和分析之后，不难发现，matlab 语言应该是一个最佳方案。 matlab 最突出的特点就是简洁。matlab 用更直观的，符合人们思维习惯的代码，代替了c 和 fortran语言的冗长代码。matlab 给用户带来的是最直观，最简洁的程序开发环境。如

果用fortran 或c 语言去为实现平台功能而编写程序，尤其当涉及矩阵运算和画图时，编程会很麻烦，如下文中的一个二元函数的三维立体图就是通过matlab 的一个简单的命令绘出来的。

以fortran 和c 这样的高级语言为例，如果用户想求解一个非线性代数方程，就得编写一个程序块读入数据，然后再使用一种求解非线性方程的算法（例如牛顿法）编写一个程序块来求解方程，最后再输出计算结果。在求解过程中，最麻烦的要算第二部分。解非线性方程的麻烦在于要对矩阵的元素作循环，选择稳定的算法以及代码的调试都不容易。即使有部分源代码，用户也会感到麻烦，且不能保证运算的稳定性。解非线性方程的程序用fortran 和c 这样的高级语言编写，估计要编写代码上百行，调试这种几百行的计算程序可以说很困难。而matlab 语言强大的库函数和优化工具箱使得这些问题迎刃而解，其优越性不言自明。 以下简单介绍一下matlab 相较于其他语言（主要指fortran 和c 语言）的主要优越性：

1）语言简洁紧凑，使用方便灵活，库函数极其丰富。matlab 程序书写形式自由，利用丰富的库函数避开繁杂的子程序编程任务，压缩了一切不必要的编程工作。由于库函数都由本领域的专家编写，用户不必担心函数的可靠性。如果用fortran 或c 语言去编写程序，尤其当涉及矩阵运算和画图时，编程会很麻烦，而matlab 的程序是极其简短的。更为难能可贵的是，matlab 甚至具有一定的智能水平，比如解方程时，matlab 会根据矩阵的特性选择方程的求解方法，所以用户根本不用怀疑matlab 的准确性。

2）运算符丰富。由于matlab 是用c 语言编写的，matlab 提供了和c 语言几乎一样多的运算符，灵活使用matlab 的运算符将使程序变得极为简短。

4）程序限制不严格，程序设计自由度大。例如，在matlab 里，用户无需对矩阵预定义就可使用。

5）程序的可移植性很好，基本上不做修改就可以在各种型号的计算机和操作系统上运行。

6）matlab 的图形功能强大。在fortran 和c 语言里，绘图都很不容易，但在matlab 里，数据的可视化非常简单。matlab 还具有较强的编辑图形界面的能力。

7）功能强大的工具箱是matlab 的另一特色。matlab 包含两个部分：核心部分和各种可选的工具箱。核心部分中有数百个核心内部函数。其工具箱又分为两类：功能性工具箱和学科性工具箱。功能性工具箱主要用来扩充其符号计算功能，图示建模仿真功能，文字处理功能以及与硬件实时交互功能。功能性工具箱用于多种学科。而学科性工具箱都是由该领域内学术水平很高的专家编写的，所以用户无需编写自己学科范围内的基础程序，而直接进行高，精，尖的研究。

8）源程序的开放性。开放性也许是matlab 最受人们欢迎的特点。除内部函数以外，所有matlab 的核心文件和工具箱文件都是可读可改的源文件，用户可通过对源文件的修改以及加入自己的文件构成新的工具箱。

综上所述，matlab 这些突出的优越性使得其成为目前实现求解非线性规划的不二选择。 3 过程论述

3.1 关键技术

整个过程中需要掌握的相关知识有线性规划、非线性规划，最优计算方法、matlab 优化工具箱的使用、matlab 编程。非线性规划求最优值的算法包括梯度法、牛顿法、拟牛顿法、方向加速法、单纯型法、粒子群法„需要对这些算法进行深入的学习与研究。并根据算法步骤编写程序。这里具体要做的工作就是利用matlab 的绘图函数绘出实例函数的立体图，利用matlab 的优化工具箱求解函数的最优值，利用matlab 语言完成非线性最优算法编程，通过实例测试各种算法的优缺点，最终寻找一种具有较强的全局搜索能力和较高的收敛精度的计算方法，能够不受初值条件限制，跳出局部最优值，以较快的速度找出全局最优解。用于解决非线性规划问题。

3.2 实验方案 整个的实验过程是一个非线性函数寻优过程，用matlab 软件编写各种最优算法程序。用最优化方法解决最优化问题的技术称为最优化技术，它包含两个方面的内容：

1） 建立数学模型 即用数学语言来描述最优化问题。模型中的数学关系式反映了最优化问题所要达到的目标和各种约束条件。

2） 数学求解 数学模型建好以后，选择合理的最优化方法进行求解。

这里选择的编程语言是matlab 。建立的数学模型就是指编写的各种算法，选用的实验方案就是对于需要输入初值的算法通过多组初值的输入，对于需要修改参数的算法，更换参数。运行程序，得要实验结果，并进行对比分析。

以下是对实验结果分析的依据：

2）迭代速度的分析：在一个算法能成功收敛到最优值的情况下，并且是同样的迭代误差范围，经过迭代次数越少的算法就是迭代速度越快的算法；或者是同样的初值，经过了同样的步数，迭代出的结果更接近最优值的算法是速度更快的算法。

判断一个算法的优劣就是看这个算法是否初值依赖，是否能够快速的收敛，一个不存在初值依赖、能快速稳定的收敛到全局最优解的算法才是好算法。

本实验通过逐一验证的方法寻找一种最优的计算方法，通过一个有代表性的测试函数对这些最优化方法进行测试，对各组结果对比分析，验证这些方法是否符合要求，直到寻找到一种符合要求的迭代方法为止。

3.3 实验过程

3.3.1 一个测试实例的初步分析

首先给出一个测试实例：

;

选择的实例是工程上所遇到的求最优问题最具代表性的例子，这是对一个二元非线性方程求其最小值，对于实际问题的应用中，这两个变量可能存在一定的线性或非线性约束。 利用matlab 的绘图功能，

在matlab 的命令窗口中输入：

f=@(x,y) 3*(1-x)^2*exp(-x^2-(y+1)^2)-10*(x/5-x^3-y^5)*exp(-x^2-y^2);

ezsurf(f,[-6 6])

可以绘出如下三维图：

从该图中可以看出：该函数有两个局部极小值，一个在点（-1.5，0）附近，一个在点（0，-2）附近，比较两点函数值的大小，得出在点（0，-2）附近的极小值是全局最优点。因为局部最小值不唯一，这就能检测出一种算法是否能跳出局部最优点，通过迭代次数能看出算法是否能快速有效地收敛到全局最优点，以下是用不同的方法求解的过程。

3.3.2 利用matlab 优化工具箱求解

function [c,ceq] = mycon(x)

c = ... % 计算x 处的非线性不等式。

ceq = ... % 计算x 处的非线性等式。

若还计算了约束的梯度，即

options = optimset('gradconstr','on')

则nonlcon 函数必须在第三个和第四个输出变量中返回c(x)的梯度gc 和ceq(x)的梯度gceq 。当被调用的nonlcon 函数只需要两个输出变量（此时优化算法只需要c 和ceq 的值，而不需要gc 和gceq ）时，可以通过查看nargout 的值来避免计算gc 和gceq 的值。 function [c,ceq,gc,gceq] = mycon(x)

c = ... % 解x 处的非线性不等式。

ceq = ... % 解x 处的非线性等式。

if nargout & 2 % 被调用的nonlcon 函数，要求有4个输出变量。

gc = ... % 不等式的梯度。

gceq = ... % 等式的梯度。

end

若nonlcon 函数返回m 元素的向量c 和长度为n 的x ，则c(x)的梯度gc 是一个n*m的矩阵，其中gc(i,j)是c(j)对x(i)的偏导数。同样，若ceq 是一个p 元素的向量，则ceq(x)的梯度gceq 是一个n*p的矩阵，其中gceq(i,j)是ceq(j)对x(i)的偏导数。其它参数意义同前。

编写目标函数objfun.m 和非线性约束函数文件confun.m （具体源程序请参见附录）

3.3.3 各种优化算法的matlab 编程

3.3.3.1梯度法

最速下降法（steepest descent method ）由法国数学家cauchy 于1847年首先提出。在每次迭代中，沿最速下降方向（负梯度方向）进行搜索，每步沿负梯度方向取最优步长，因此这种方法称为最优梯度法。最速下降法是一种最基本的算法，它在最优化方法中占有重要地位. 最速下降法的优点是工作量小，存储变量较少, 初始点要求不高；缺点是收敛慢，最速下降法适用于寻优过程的前期迭代或作为间插步骤，当接近极值点时，宜选用别种收敛快的算法.

用最速下降法（梯度法）求解该非线性最优规划：

其基本迭代格式

（6）

总是考虑从点 出发沿哪一个方向 ，使目标函数 下降得最快。微积分的知识告诉大家，点 的负梯度方向

dfp 变尺度法的计算步骤总结如下：

根据上面的步骤编写matlab 程序dfp.m 和dfun.m （具体源程序请参见附录）

对于正定二次函数

这个方法主要由所谓基本搜索、加速搜索和调整搜索方向三部分组成。

具体步骤如下：

3.3.3.5单纯型法 单纯形指的是ｎ维空间中的具有ｎ＋１个顶点的多面体。通过计算单纯形各顶点的目标函数值并加以比较，从中确定有利的搜索方向和步长，找到一个较好的点取代单纯形中较差的点，组成新的单纯形来代替原来的单纯形，使单纯形不断向目标函数的极小点靠近，直到搜索到极小点为止。这就是单纯形法的基本思想。

单纯形法属于最优化方法中的无约束类中的直接搜索法，算法简单、收敛速度较快，适用范围广泛，无需通过计算目标函数的导数或梯度来寻求目标下降的最优方向。但单纯形法的求解性能依赖于初始点的选择，同时采用的是确定性搜索机制，算法陷入局部极值点后难

以跳出。

单纯形算法计算步骤如下：

根据以上步骤编写matlab 程序dcx.m （具体源程序请参见附录）

3.3.3.6模拟退火算法 模拟退火算法是一种通用的随机搜索算法，是局部搜索算法的扩展。它的思想是再1953年由metropolis 提出来的，到1983年由kirkpatrick 等人成功地应用在组合优化问题中。

模拟退火算法来源于固体退火原理，将固体加温至充分高，再让其徐徐冷却，加温时，固体内部粒子随温升变为无序状，内能增大，而徐徐冷却时粒子渐趋有序，在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。根据metropolis 准则，粒子在温度t 时趋于平衡的概率为e-δe/(kt)，其中e 为温度t 时的内能，δe 为其改变量，k 为boltzmann 常数。用固体退火模拟组合优化问题，将内能e 模拟为目标函数值f ，温度t 演化成控制参数t ，即得到解组合优化问题的模拟退火算法：由初始解i 和控制参数初值t 开始，对当前解重复“产生新解→计算目标函数差→接受或舍弃”的迭代，并逐步衰减t 值，算法终止时的当前解即为所得近似最优解，这是基于蒙特卡罗迭代求解法的一种启发式随机搜索过程。退火过程由冷却进度表(cooling schedule)控制，包括控制参数的初值t 及其衰减因子δt 、每个t 值时的迭代次数l 和停止条件s 。

模拟退火算法新解的产生和接受可分为如下四个步骤： 第一步是由一个产生函数从当前解产生一个位于解空间的新解；为便于后续的计算和接受，减少算法耗时，通常选择由当前新解经过简单地变换即可产生新解的方法，如对构成新解的全部或部分元素进行置换、互换等，注意到产生新解的变换方法决定了当前新解的邻域结构，因而对冷却进度表的选取有一定的影响。 第二步是计算与新解所对应的目标函数差。因为目标函数差仅由变换部分产生，所以目标函数差的计算最好按增量计算。事实表明，对大多数应用而言，这是计算目标函数差的最快方法。 第三步是判断新解是否被接受, 判断的依据是一个接受准则，最常用的接受准则是metropo1is 准则: 若δt ′<0则接受s ′作为新的当前解s ，否则以概率exp(-δt ′/t)接受s ′作为新的当前解s 。 第四步是当新解被确定接受时，用新解代替当前解，这只需将当前解中对应于产生新解时的变换部分予以实现，同时修正目标函数值即可。此时，当前解实现了一次迭代。可在此基础上开始下一轮试验。而当新解被判定为舍弃时，则在原当前解的基础上继续下一轮试验。 模拟退火算法与初始值无关，算法求得的解与初始解状态s(是算法迭代的起点) 无关；模拟退火算法具有渐近收敛性，已在理论上被证明是一种以概率l 收敛于全局最优解的全局优化算法；模拟退火算法具有并行性。 模拟退火算法的步骤： (1) 初始化：初始温度t(充分大) ，初始解状态s(是算法迭代的起点) ， 每个t 值的迭代次数l (2) 对k=1，„„，l 做第(3)至第6步： (3) 产生新解s ′ (4) 计算增量δt ′=c(s′)-c(s)，其中c(s)为评价函数 (5) 若δt ′<0则接受s ′作为新的当前解，否则以概率exp(-δt ′/t)接受s ′作为新的当前解.

(6) 如果满足终止条件则输出当前解作为最优解，结束程序。 终止条件通常取为连续若干个新解都没有被接受时终止算法。 (7) t逐渐减少，且t-&0，然后转第2步。

根据有关步骤编写matlab 程序mnth.m （请 参考 附录中的源程序）。

4 结果分析

4.1各种优化算法的测试 1）利用matlab 优化工具箱解函数最优值的结果分析： 在命令窗口中输入如下的命令：

&& clear

&& x0=[0 -2];

&& options=optimset('largescale','off','display','iter');

&&

[x,fval,exitflag,output]=fmincon(@objfun,x0,[],[],[],[],[],[],@confun,options) 输出结果如下：

max directional first-order

iter f-count f(x) constraint step-size derivative optimality procedure 1 11 -6.43488 -9.86 0.0625 5.74 1.88

2 19 -6.52097 -9.729 0.0625 1.4 0.711

3 23 -6.54274 -9.776 1 0.02 0.407

4 27 -6.54589 -9.772 1 8.77e-005 0.0112

5 31 -6.54589 -9.771 1 4.61e-007 0.000683

optimization terminated successfully:

magnitude of directional derivative in search direction

less than 2*options.tolfun and maximum constraint violation

is less than options.tolcon

no active constraints

x =

0.2289 -1.6261

fval =

-6.5459

exitflag =

1

output =

iterations: 5

funccount: 31

stepsize: 1

firstorderopt: 6.8256e-004

cgiterations: []

&&

通过matlab 优化工具箱提供的优化函数求得最终的极值点为（0.2289，-1.6261），极值为fval =-6.5459。误差值为6.8256e-004。

在上面的输入中，令x0=[-2 0]，重复上面的操作，输出的极值点为（-1.7775，1.1800），极值为fval =-0.3050。从前面的立体图上可以看出这是例外的一个局部极小值，显然利用matlab 自带的最优工具箱的fmincon 函数求最优值时，对初值的依赖很严重，只有在靠近全局极值点附近的初值才能收敛到最优极值点。可见这种求解最优值的方法很不稳定。