非线性规划的算法研究
关键词：非线性规划 最优决策 初值依赖 matlab
论文摘要：本课题主要研究非线性规划的算法。非线性规划在军事，，管理，生产过程自动化，工程设计和产品优化设计等方面都有着重要的应用。但非线性规划的研究目前还不成熟，有许多问题需要进一步完善。非线性规划不像线性规划有统一的算法，对于不同的问题需要用不同的算法处理，现阶段各种算法都有一定的局限性，只有对各种算法加以改正，才能有效地解决人们在日常的生产、生活中遇到的优化问题，做出最优决策。
本文主要是对现有的各种算法加以测试，指出各种算法的优缺点，寻找一种不受初值依赖，收敛更快的最优算法。首先介绍了非线性规划研究的背景和国内外研究状况，然后论述了方案的选取过程，重点描实验过程，主要是对各种非线性最优方法用matlab软件编程，给出一个在工程中具有代表性的最优函数实例，经过大量的测试，并给出了结果分析。最后给出了整个实验的和由此对未来的展望。
1 选题背景
1.1 课题背景
非线性规划的一个重要理论是1951年Kuhn-Tucker最优条件(简称KT条件)的建立.此后的50年代主要是对梯度法和牛顿法的研究.以Davidon(1959), Fletcher和Powell(1963)提出的DFP方法为起点, 60年代是研究拟牛顿方法活跃时期, 同时对共轭梯度法也有较好的研究. 在1970年由Broyden,Fletcher,Goldfarb 和Shanno从不同的角度共同提出的BFGS方法是目前为止最有效的拟牛顿方法. 由于Broyden, Dennis 和More的工作使得拟牛顿方法的理论变得很完善. 70年代是非线性规划飞速时期, 约束变尺度(SQP)方法(Han和Powell为代表)和Lagrange乘子法(代表人物是Powell 和Hestenes)是这一时期主要研究成果.计算机的飞速发展使非线性规划的研究如虎添翼.80年代开始研究信赖域法、稀疏拟牛顿法、大规模问题的方法和并行计算, 90年代研究解非线性规划问题的内点法和有限储存法. 可以毫不夸张的说, 这半个世纪是最优化发展的黄金时期. 以后随着计算机的普遍使用，非线性规划的理论和方法有了很大的发展，其应用的领域也越来越广泛，特别是在军事，经济，管理，生产过程自动化，工程设计和产品优化设计等方面都有着重要的应用。
1.2 课题研究的目的和意义
一般来说，解非线性规划问题要比求解线性规划问题困难得多，而且也不像线性规划那样有统一的数学模型及如单纯形法这一通用解法。非线性规划的各种算法大都有自己特定的适用范围。都有一定的局限性，到目前为止还没有适合于各种非线性规划问题的一般算法。这正是需要人们进一步研究的课题。非线性规划在工程、管理、经济、科研、军事等方面都有广泛的应用，为最优设计提供了有力的工具。例如：如何在现有人力、物力、财力条件下合理安排产品生产，以取得最高的利润；如何设计某种产品，在满足规格、性能要求的前提下，达到最低的成本；如何确定一个自动控制系统的某些参数，使系统的工作状态最佳；如何分配一个动力系统中各电站的负荷，在保证一定指标要求的前提下，使总耗费最小；如何安排库存储量，既能保证供应，又使储存费用最低；如何组织货源，既能满足顾客需要，又使资金周转最快等。对于静态的最优化问题，当目标函数或约束条件出现未知量的非线性函数，且不便于线性化，或勉强线性化后会招致较大误差时，就可应用非线性规划的方法去处理。
1.3 国内外概况
解决最优化问题的算法分为经典优化算法和启发式优化算法。经典算法的确立可以从（G.B.Dantzig 1947）提出解决线形规划的单纯形（Simplex method）开始，但单纯形不是多项式算法。随后，Kamaka提出了椭球算法(多项式算法)，内点法。对于非线性问题，起初人们试图用线性优化理论去逼近求解非线性问题，但效果并不理想。后来的非线性理论大多都建立在二次（凸）函数的基础上，也就用二次函数去逼近其他非线性函数。在此基础上提出许多优化经典的优化算法。无约束的优化算法包括：最速下降法（steepest）、共轭梯度法、牛顿法（Newton Algorithm）、拟牛顿法（pseudo Newton Algorithms）、信赖域法。约束优化算法包括：拉格朗日乘子法（Augmented Lagrangian Algorithms），序列二次规划（ＳＱＰ）等。
随着社会的发展，实际问题越来越复杂，例如全局最优化问题。经典算法一般都用得局部信息，如单个初始点及所在点的导数等，这使得经典算法无法避免局部极小问题。全局最优化是NP-Hard问题,所以原有的经典算法不再使用，必须对其进行改进，或将其与启发式算法结合。启发式算法是受大的启发，人们从大自然的运行中找到了许多解决实际问题的方法。启发式算法的计算量都比较大，所以启发式算法伴随着计算机技术的发展，取得了巨大的成就。
40年代：由于实际需要，人们已经提出了一些解决实际问题快速有效的启发式算法。
50年代：启发式算法的研究逐步繁荣起来。随后，人们将启发式算法的思想和人工智能领域中的各种有关问题的求解的收缩方法相结合，提出了许多启发式的搜索算法。其中贪婪算法和局部搜索等到人们的关注。
60年代: 随着人们对数学模型和优化算法的研究越来越重视，发现以前提出的启发式算法速度很快，但是解得质量不能保证。虽然对优化算法的研究取得了很大的进展，但是较大规模的问题仍然无能为力（计算量还是太大）。
70年代：计算复杂性理论的提出。NP完全理论告诉我们，许多实际问题不可能在合理的时间范围内找到全局最优解。发现贪婪算法和局部搜索算法速度快，但解不好的原因主要是他们只是在局部的区域内找解，得到的解不能保证全局最优性。由此必须引入新的搜索机制和策略，才能有效地解决这些困难问题。
80年代以后： 模拟退火算法（Simulated Annealing Algorithm），人工神经（Artificial Neural Network），禁忌搜索（Tabu Search）相继出现。最近，演化算法（Evolutionary Algorithm）, 蚁群算法（Ant Algorithms）， 拟人拟物算法，量子算法等油相继兴起，掀起了研究启发式算法的高潮。由于这些算法简单和有效，而且具有某种智能，因而成为计算和人类之间的桥梁。
优胜劣汰是大自然的普遍规律，它主要通过选择和变异来实现。选择是优化的基本思想，变异（多样化）是随机搜索或非确定搜索的基本思想。“优胜劣汰”是算法搜索的核心，根据“优胜劣汰”策略的不同，可以获得不同的超启发式算法。超启发式算法的主要思想来自于人类经过长期对物理、生物、社会的自然现象仔细的观察和实践，以及对这些自然现象的深刻理解，逐步向大自然学习，模仿其中的自然现象的运行机制而得到的。
遗传算法：是根据生物演化，模拟演化过程中基因染色体的选择、交叉和变异得到的算法。在进化过程中，较好的个体有较大的生存几率。
模拟退火：是模拟统计物理中固体物质的结晶过程。在退火的过程中，如果搜索到好的解接受；否则，以一定的概率接受不好的解（即实现多样化或变异的思想），达到跳出局部最优解得目的。
神经网络：模拟大脑神经处理的过程，通过各个神经元的竞争和协作，实现选择和变异的过程。
禁忌搜索：模拟人的经验，通过禁忌表记忆最近搜索过程中的信息，禁忌某些解，以避免走回头路，达到跳出局部最优解的目的。
蚂蚁算法：模拟蚂蚁的行为，拟人拟物，向蚂蚁的协作方式学习。
虽然人们研究对启发式算法的研究将近５０年，但它还有很多不足：
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的运算符将使程序变得极为简短。
3）MATLAB既具有结构化的控制语句（如for循环，while循环，break语句和if语句），又有面向对象编程的特性。
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。建立的数学模型就是指编写的各种算法，选用的实验方案就是对于需要输入初值的算法通过多组初值的输入，对于需要修改参数的算法，更换参数。运行程序，得要实验结果，并进行对比分析。
以下是对实验结果分析的依据：
1）初值依赖的分析：大多数算法都需要输入初值（如梯度法、newton迭代法、拟newton迭代法、方向加速法、单纯形法），。如果不同的初值能得到不同的迭代结果，就说明改算法对初值依赖，如果不同的初值都能收敛到同一个全局最优解，就说明改算法不存在初值赖；另一种情况是不需要手动输入初值（如粒子群法），初值是由一些随机函数产生的，这些初值点在定义域内产生是等概率的，多次实验都能收敛到同一结果，则说明该算法不存在初值依赖。
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）
总是考虑从点 出发沿哪一个方向 ，使目标函数 下降得最快。微积分的知识告诉大家，点 的负梯度方向
根据上述步骤，编写newton.m、nwfun.m和detaf.m文件。（具体源程序请参见附录）
3.3.3.3拟牛顿法
变尺度法（Variable Metric Algorithm）是近20多年来发展起来的，它不仅是求解无约束极值问题非常有效的算法，而且也已被推广用来求解约束极值问题。由于它既避免了计算二阶导数矩阵及其求逆过程，又比梯度法的收敛速度快，特别是对高维问题具有显著的优越性，因而使变尺度法获得了很高的声誉。下面简要地介绍一种变尺度法—DFP法的基本原理及其计算过程。这一方法首先由Davidon在1959年提出，后经Fletcher和Powell加以改进。
这就是常说的拟Newton条件。
DFP变尺度法的计算步骤总结如下：
根据上面的步骤编写matlab程序dfp.m和dfun.m（具体源程序请参见附录）
3.3.3.4 Powell法
方向加速法是有powell（1964年）提出来的，因此又称powell方法。迄今为止，它是直接方法中最有效的方法。
对于正定二次函数
在不用导数的前提下，powell方法逐次迭代够造A共轭方向系，因此powell方法本质上属于共轭方向法，具有二次收敛性。用于非二次函数时一般也具有较快的收敛速度。
这个方法主要由所谓基本搜索、加速搜索和调整搜索方向三部分组成。
具体步骤如下：
根据以上步骤，编写matlab程序powell.m（具体源程序请参见附录）
3.3.3.5单纯型法
单纯形指的是ｎ维空间中的具有ｎ＋１个顶点的多面体。通过计算单纯形各顶点的目标函数值并加以比较，从中确定有利的搜索方向和步长，找到一个较好的点取代单纯形中较差的点，组成新的单纯形来代替原来的单纯形，使单纯形不断向目标函数的极小点靠近，直到搜索到极小点为止。这就是单纯形法的基本思想。
单纯形法属于最优化方法中的无约束类中的直接搜索法，算法简单、收敛速度较快，适用范围广泛，无需通过计算目标函数的导数或梯度来寻求目标下降的最优方向。但单纯形法的求解性能依赖于初始点的选择，同时采用的是确定性搜索机制，算法陷入局部极值点后难以跳出。
相对于newton法、powell法等其他常规经典优化算法而言，单纯形法具有群体搜索的特性，算法在每一步的迭代过程中，针对的是单纯形这一群体，而不是单独的某个点，因此具有很好的并行求解特性。
单纯形算法计算步骤如下：
根据以上步骤编写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
algorithm: 'medium-scale: SQP, Quasi-Newton, line-search'
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函数求最优值时，对初值的依赖很严重，只有在靠近全局极值点附近的初值才能收敛到最优极值点。可见这种求解最优值的方法很不稳定。
2）梯度法求解非线性函数最优值的实验结果为：在matlab的命令窗口输入zuisu并回车，
初值x=[-2;0]时，最后得到：
x =
0.2288
-1.6261
k =
11
f0 =
-6.5459
g =
1.0e-006 *
-0.2747 0.7430
当令重设初值x=[-2;0]时，重复上述操作，得到的结果为：
x =
-1.3597
0.2145
k =
14
f0 =
-2.7683
g =
1.0e-006 *
-0.6480 0.2380
当初值为x=[0,-2]时，经过11次迭代，求得的极值点为全局极值点（0.2288， -1.6261），极值为f0 = -6.5459。令重设初值x=[-2;0]时，经过的14次迭代，求得的极值点为另外一个局部极值点（-1.3597，0.2145），极值为f0 = -2.7683，，说明梯度法求非线性最优值受初值的影响很大，速度也较慢，只有在全局最优解附近的初值才有可能收敛到最终要求的极值点，当一个函数有多个极值点的时候，就有可能陷入局部极值点而无法跳出。需要寻找一种更好的方法。
3）牛顿法求解非线性函数最优值的实验结果为：
在matlab的命令窗口输入newton并回车
初值x=[-2;0]时，最后得到：
x =
0.2288
-1.6261
k =
11
f0 =
-6.5459
g1 =
1.0e-006 *
-0.0391 0.2705
g2 =
19.2367 -3.8092
-3.8092 28.9026
当令初值x=[-2;0]时，重复上述操作，运行1001次后停止。
x =
-2.0000
-0.0000
k =
1001
f0 =
-1.2101
g1 =
-2.8003 -0.3638
g2 =
-2.7786 -1.2128
-1.2128 3.1478
当令初值x=[-1.5;0]时，重复上述操作，得到的结果为：
x =
-1.3597
0.2145
k =
11
f0 =
-2.7683
g1 =
1.0e-007 *
0.7809 -0.8835
g2 =
13.9596 -2.2092
-2.2092 8.7492
当令newton法的迭代初值为x=[2;-2]时，运行1001次后停止:
x =
2
-2
k =
1001
f0 =
-0.0616
g1 =
0.3266 -0.0186
g2 =
-1.4224 0.3136
0.3136 0.5047
当令初值为x=[-2;0]和x=[2;-2]时，程序运行了1001次后停止。经对newton法分析，其迭代步长t最终因为不停减半变得太小导致计算机将其近视为零处理，导致每一次迭代完之后仍然停留在原初值点，迭代失败，该算法需要进一步改进。
当初值为（-1.5，0）时，求得的极值点为另外一个局部极值点（-1.3597，0.2145），极值为f0 = -2.7683，迭代次数为11次。当初值为（0，-2）时，求得的极值点为全局极值点（0.2288， -1.6261），极值为f0 = -6.5459，迭代次数为11次。说明newton法求非线性最优值受初值的影响非常大，只有在全局最优解很近的初值才有可能收敛到最终要求的极值点，当初值设置不好有可能无法收敛。当一个函数有多个极值点的时候，就有可能陷入局部极值点而无法跳出。newton法的优点是收敛速度快；缺点是有时不好用而需采取改进措施，此外，当维数较高时，计算 的工作量很大。需要寻找一种比newton更好的计算方法。
4）拟牛顿法求解非线性函数最优值的实验结果为：
在matlab的命令窗口输入dfp并回车
初值x=[-2;0]时，最后得到：
k =
4
opt =
0.2288
-1.6261
f =
-6.5459
当令初值x=[-2;0]时，重复上述操作，得到的结果为：
k =
5
opt =
-1.3597
0.2145
f =
-2.7683
当令初值x=[2;-2]时，重复上述操作，得到的结果为：
k =
7
opt =
0.2288
-1.6261
f =
-6.5459
当初值为x=[0,-2]和x=[2;-2]时，均能收敛到全局最优值，由此可以得出结论：拟牛顿法中的dfp法相对于没有改进前的牛顿法对初值的依赖程度相对减弱，初值的选择范围变大；对上面的选择的初值分别经过4次、7次收敛到全局最优解，相对于牛顿法收敛速度也有加快。但是它的初值如果选择不好仍然会收敛到局部最优解最终无法跳出，得不到全局最优解。但此方法相对于牛顿法已有明显改善。
5）Powell法求解非线性函数最优值的实验结果为：
在matlab的命令窗口输入powell并回车
初值x=[-2;0]时，最后得到：
count =
4
lastx =
0.2288 -1.6261
lastf =
-6.5459
当令初值x=[-2;0]时，重复上述操作，得到的结果为：
count =
4
lastx =
-0.8060 0
lastf =
-2.7560
当令初值x=[2;-2]时，重复上述操作，得到的结果为：
count =
4
lastx =
-1.6420 -2.0000
lastf =
-6.5227
可见powell法所迭代出的最优值也跟初值有关，甚至可能求得的点不是局部最优点，只有再初值相当接近全局最优点的时候才有可能收敛到全局最优点，可见powell法也很不稳定，需要继续寻求一种好的最优算法。
6）单纯形法求解非线性函数最优值的实验结果为：
在matlab命令窗口输入dxc并回车
初值x=[-2;0]时，得到的最后结果为：
xo =
0.228846942455879 -1.626050258682904
xr =
0.228846938467961 -1.626050254103438
N =
1000
ans =
-6.545892011017174
经过一千次迭代后，收敛到的全局最优点（0.228846942455879，-1.626050258682904）
最优值为-6.545892011017174，若将循环条件修改为N<30，则得到的结果为：
xc =
0.234177084377734 -1.625755309876695
N =
30
ans =
-6.545115886509709
如果将初值迭代点x1改为（-1.3，0.2），此时三个初值点构成的三角区域将包括另外一个局部极值点（-1.3597，0.2145），，运行程序，得到的结果为：
xo =
0., , 228846962960411 -1.626050259455931
xr =
0.228846962960411 -1.626050259455933
N =
1000
ans =
-6.545892011017175
依然能很好的收敛到全局最优解，当将迭代初值设置为x1（-2，4）、x2（-3，-4）、x1（3，-4），运行程序，得到的结果为：
xo =
10.250000000000000 -24.000000000000000
xr =
17.500000000000000 -44.000000000000000
N =
1000
ans =
0
程序不能收敛，迭代失败。
由此可以得出结论：单纯形算法的收敛速度相对较慢，迭代过程中跳出局部最优解的能力也有所加强，但是受初值的影响依然很大，如果初值选择不当，将导致初值点发散，无法收敛，最终迭代失败，所以单纯形算法还算不上理想的寻优算法。
7）模拟退火算法求解非线性函数最优值的实验结果为：
在matlab命令窗口输入以下命令：
>> loss= @(x) 3*(1-x(1))^2*exp(-x(1)^2-(x(2)+1)^2)-10*(x(1)/5-x(1)^3-x(2)^5)*exp(-x(1)^2-x(2)^2);
>> [x f] = mnth(loss,[0 -2])
得到的最后结果为：
Initial temperature: 1
Final temperature: 3.6185e-007
Consecutive rejections: 1254
Number of function calls: 7470
Total final loss: -6.54589
x =
0.2289 -1.6261
f =
-6.5459
当输入>> [x f] = mnth(loss,[-2 0])
得到的最后结果为：
Initial temperature: 1
Final temperature: 2.05688e-007
Consecutive rejections: 1058
Number of function calls: 6634
Total final loss: -2.76829
x =
-1.3599 0.2144
f =
-2.7683
当输入>> [x f] = mnth(loss,[2 -2])
得到的最后结果为：
Initial temperature: 1
Final temperature: 1.6455e-007
Consecutive rejections: 1117
Number of function calls: 8231
Total final loss: -6.54589
x =
0.2289 -1.6260
f =
-6.5459
当初始点为（0，-2）和点（2，-2）时，都能收敛到全局最优解，当初始点为[-2,0]时，收敛到局部最优点（-1.3599，0.2144），说明需要调节一些参数才能使得迭代跳出局部最优值，这需要大量的数值模拟计算来选择好参数。
4.2 结果
对各种非线性最优算法的有缺点列表进行对比：
算法 | 优点 | 缺点 |
梯度法 | 计算量小，存储变量较少,初始点要求不高 | 初值依赖，收敛慢，最速下降法适用于寻优过程的前期迭代或作为间插步骤，越接近极值点时，收敛熟读越慢，后期宜选用收敛快的算法 |
牛顿法 | 收敛速度很快 | 当维数较高时，计算 |
拟牛顿法 | 收敛速度快，避免牛顿矩阵求逆运算，算法更稳定 | 初值依赖程度相对牛顿法减弱，但仍然存在 |
方向加速法 | 收敛速度较快 | 初值依赖，算法不稳定，有可能迭代到一个一般的点（非最优点）就停止迭代 |
单纯形法 | 算法简单、收敛速度较快，适用范围广泛，无需通过计算目标函数的导数或梯度来寻求目标下降的最优方向，具有群体搜索的特性 | 依赖于初始点的选择，同时采用的是确定性搜索机制，算法陷入局部极值点后难以跳出。当初值选择不好时，经过几次循环后初值点坑发散开 |
模拟退火算法 | 初值依赖小，不用求目标函数的偏导数及解大型矩阵方程组，即能找到一个全局最优解，而且易于加入约束条件，编写程序简单 | 算法通常需要较高的初温、较慢的降温速率、较低的终止温度以及各温度下足够多的抽样次数，因而模拟退火算法往往过程较长，实现效率较低，参数调节麻烦 |
的工作量很大，初值依赖，当初值选择不好时，有可能计算出现异常，导致迭代无法进行，该法需要修正
5 总结与展望
5.1 设计总结
总体来说，这次毕业设计的任务是基本完成的，这次毕业设计的作品——氧化铝物料平衡计算软件开发。我所做的工作主要是最优算法这一部分，从开始的用matlab工具箱的优化函数到梯度法、牛顿法、拟牛顿法、powell法、单纯形法、模拟退火算法，这些方法都通过了
这个实例的一一验证，分别得出了各种方法的优缺点，最终找到一种行之有效的方法。
但是这个验证过程还存在一些不完善的地方，测试的函数只有一个实例，这个实例只有两个局部最小值，可能在实际的工程项目中所遇到的函数更为复杂，局部最小值也更多，粒子群法还能否稳定快速地收敛到全局最优值需要进一步加以验证，如果不能就需要进一步对该算法加以改正，或者继续验证其它的最优算法。
无论如何，这次毕业设计都给了我很大的锻炼，从进行一项设计本身开始对论题背景、目的和意义的分析，到对论题方案的建立、选取和论证，再到方案的施行和最终的总结和展望，这一切都是对我这大学四年中所学到的知识和能力的巨大考验和提高。在这个过程中，从最开始的线性规划、非线性规划的各种迭代算法的学习，到掌握matlab编程。虽然只是整个设计中的一部分工作，但是，很大地提高了我分析问题、解决问题的能力。
值得一提的是，这次毕业设计中的锻炼使我对MATLAB语言有了较深刻的了解和认识：MATLAB有很多显而易见的优越性，但是MATLAB也有它的缺点，即由于MATLAB的程序不用编译等预处理，也不生成可执行文件，程序为解释执行，所以和其他高级程序相比起来，程序的执行速度较慢。不过这一点小小的不足并不能影响MATLAB将来在更多的领域中有更大的作为。
5.2 设计展望
通过这次毕业设计我深刻的感受到非线性最优问题的研究仍然任重道远，非线性规划还不成熟，不像线性规划那样有一般的求解方法，要对具体问题具体分析，选择合适的最优算法求解。需要更多的科研人员投身到非线性规划的研究，不断的对各种算法加以改进，不断的设计更加优秀的算法，用来解决各个领域的最优规划问题。
此外，MATLAB作为一个非常优秀的软件，其巨大的优越性也越来越突现出来，如果简单的做一个展望的话，MATLAB今后必将在社会生产的各个领域有更多更广运用，并将产生更深远的影响。
所以在我看来，MATLAB作为一个今后必备的数学工具应用软件，完全可以在各大专院校的理工科高年级或有高级编程语言的基础之后开必修课，以更好的完善我国理工科高等人才的综合能力，提高他们思考、分析和解决问题的方式和能力，使得他们在社会中更好更高效率的学习和工作。
[1]黄象鼎.非线性数值分析.武汉大学出版社.2000
[2]张光澄.非线性最优化计算方法.高等出版社.2005
[3]魏诺.非线性基础与应用.科学出版社.2004
[4]王凌.智能优化算法及其应用.北京:清华大学出版社.2000
[5]曹建福,韩嵩昭.非线性系统理论及应用.西安大学出版社.2002
[6]韩茂安,顾圣土.非线性系统的理论和方法.科学出版社.2004
[7]张圣勤.MATLAB 7.0实用敎程.机械出版社.2006
[8]王家文.MATLAB 7.0编程基础.机械工业出版社.2005
[9]田宝国,谷可,姜璐. 从线性到非线性——科学的历程. 系统辩证学学报 .第9卷第3期.2001年7月
[10]李志斌.求解非线性规划问题的一种迭代法 [J] 大连铁道学院学报,第19卷.第4期.1998年12月
[11]黄平，沈沉，陈颖.非线性代数方程组异步迭代算法研究及其在潮流计算中的应用.
[12]沈飞.非线性规划单纯形算法的改进算法[J] 测绘信息与工程.Journal of Geomatics.2004年8月
[13]ZHIYOU WU，LIANSHENG ZHANG.一些类型的数学规划问题的全局最优解. 运筹学学报.第7卷第2期.2003年5月
[14]胡运红.非线性最优化问题及其算法研究[J] 运城学院学报.第2l卷第3期 .2OO3年6月
[15]简金宝.非线性规划问题的一个全局收敛的次可行方法[J] 曲阜师范大学学报.第18卷第4期.1992年10月
[16]陈加民，王希云.求解非线性最优化问题的Newton算法研究[J] 井冈山学院学报(科学).第28卷第8期.2007年8月
[17]夏念凌.非线性规划问题的实用优化算法[J] 安徽省水利水电设计院.第10卷第4-5期.1989年4月
[18]吴至友.目标混合型且带约束的多目标最优决策及其应用[J].重庆师范学院学报.第16巷第3期.1999年9月
[19]童小娇 何伟. 解非线性约束方程的拉格朗日全局投影方法.数学物报.2008年8月
[20]闫利军,李宗斌,卫军胡.模拟退火算法的一种参数设定方法研究.系统仿真学报[J].第20卷第1期.2008年
[21]王福昌，张艳芳.一种改进模拟退火算法在非线性方程组求解中的应用.航空计算技术.第37卷第6期.2007年
[22]陈利民，苏宏业，牟盛静，褚健.基于有界变量单纯形法的改进区间牛顿法.浙江大学学报.第37卷第3期.2003年
