1、 -32- 第三章 非线性规划 1 非线性规划 1.1 非线性规划的实例与定义 如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不象线性规划有单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。 下面通过实例归纳出非线性规划数学模型的一般形式， 介绍有关非线性规划的基本概念。 例 1 （投资决策问题）某企业有 n 个项目可供选择投资，并且至少要对其中一个项目投资。已知该企业拥有总资金 A元，投资于第 ),1( nii L= 个项目需花资金ia 元，并预计可收益ib
2、元。试选择最佳投资方案。 解 设投资决策变量为 =个项目决定不投资第，个项目决定投资第iixi0,1， ni ,1 L= ， 则投资总额为=niiixa1，投资总收益为=niiixb1。因为该公司至少要对一个项目投资，并且总的投资金额不能超过总资金 A，故有限制条件 =niiiAxa10 另外，由于 ),1( nixiL= 只取值 0 或1，所以还有 .,1,0)1( nixxiiL= 最佳投资方案应是投资额最小而总收益最大的方案， 所以这个最佳投资决策问题归结为总资金以及决策变量（取 0 或 1）的限制条件下，极大化总收益和总投资之比。因此，其数学模型为： =niiiniiixaxbQ11m
3、ax s.t. =niiiAxa10 .,1,0)1( nixxiiL= 上面例题是在一组等式或不等式的约束下，求一个函数的最大值（或最小值）问题，其中至少有一个非线性函数，这类问题称之为非线性规划问题。可概括为一般形式 )(min xf qjxhj,1,0)(s.t. L= (NP) pixgi,1,0)( L= 各位亲，插个小小的广告，如果不感兴趣，请跳过这一页 代写数学建模论文 数学建模意义非凡！众所周知，在数模竞赛中获奖意义重大。它是对个人能力的肯定，是包括计算机、数学、英语等知识和团队合作等多项能力的体现，是保研或出国留学的通行证！ 但是获奖是不容易的！本人系清华大学数学系保送研究生
4、，本科期间获得过国家赛一等奖和国际赛一等奖。本团队提供代写数模论文服务，我团队已有6年代写经验，已经为全国范围内的参赛考生获取过8次国家一等奖，14次国家二等奖，获得美赛一等奖6次，二等奖18次。本人QQ号1105758397，如有购需要，可以和我联系，价格QQ商议。 同时本人还代做毕业论文，算法设计和实现，matlab实现算法等，有需要可加我QQ。 -33- 其中Tnxxx 1L= 称为模型 （NP） 的决策变量， f 称为目标函数，ig ),1( pi L=和 ),1( qjhjL= 称为约束函数。另外， 0)( =xgi ),1( pi L= 称为等式约束，0)( xhj ),1( qj
5、 L= 称为不等式的约束。 对于一个实际问题，在把它归结成非线性规划问题时，一般要注意如下几点： （i）确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基础上，确认什么是问题的可供选择的方案，并用一组变量来表示它们。 （ii）提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。 （iii）给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或“坏”的价值标准，并用某种数量形式来描述它。 （iv）寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果，因此还需要寻找出问题
6、的所有限制条件，这些条件通常用变量之间的一些不等式或等式来表示。 1.2 线性规划与非线性规划的区别 如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到） ；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。 1.3 非线性规划的 Matlab 解法 Matlab 中非线性规划的数学模型写成以下形式 )(min xf =0)(0)(xCeqxCBeqxAeqBAx, 其中 )(xf 是标量函数， BeqAeqBA , 是相应维数的矩阵和向量， )(),( xCeqxC 是非线性向量函数。 Matlab 中的命令是 X=FMINCON(FU
7、N,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS) 它的返回值是向量 x ，其中 FUN 是用 M 文件定义的函数 )(xf ；X0 是 x 的初始值；A,B,Aeq,Beq 定义了线性约束 BeqXAeqBXA = *,* ，如果没有线性约束，则A=,B=,Aeq=,Beq=；LB 和 UB 是变量 x 的下界和上界，如果上界和下界没有约束，则LB=，UB=，如果 x 无下界，则 LB 的各分量都为-inf，如果 x 无上界，则 UB的各分量都为 inf； NONLCON 是用M 文件定义的非线性向量函数 )(),( xCeqxC ； OPTIONS定义了优化参数
8、，可以使用 Matlab缺省的参数设置。 例2 求下列非线性规划 8)(min 232221+= xxxxf s.t. 023221+ xxx 2033221+ xxx 02221=+ xx -34- 32232=+ xx 0,321xxx 解 （i）编写 M 文件fun1.m 定义目标函数 function f=fun1(x); f=sum(x.2)+8; （ ii）编写 M文件 fun2.m定义非线性约束条件 function g,h=fun2(x); g=-x(1)2+x(2)-x(3)2 x(1)+x(2)2+x(3)3-20; %非线性不等式约束 h=-x(1)-x(2)2+2 x(
9、2)+2*x(3)2-3; %非线性等式约束 （ iii）编写主程序文件 example2.m 如下： options=optimset(largescale,off); x,y=fmincon(fun1,rand(3,1),zeros(3,1), . fun2, options) 就可以求得当 9478.0,2033.1,5522.0321= xxx 时，最小值 6511.10=y 。 1.4 求解非线性规划的基本迭代格式 记（ NP）的可行域为 K 。 若 Kx *，并且 Kxxfxf ),()(* 则称*x 是（ NP）的整体最优解， )(*xf 是 (NP)的整体最优值。如果有 *,)
10、,()( xxKxxfxf 则称*x 是（ NP）的严格整体最优解， )(*xf 是 (NP)的严格整体最优值。 若 Kx *，并且存在*x 的邻域 )(*xN，使 KxNxxfxf I)(),()(* ， 则称*x 是（ NP）的局部最优解， )(*xf 是 (NP)的局部最优值。如果有 KxNxxfxf I)(),()(* ，使 ),0(),()( t ，使 Ktpx + ， 称向量 p 是点 x 处关于 K 的可行方向。 一个向量 p ，若既是函数 f 在点 x 处的下降方向，又是该点关于区域 K 的可行方向，称之为函数 f 在点 x 处关于 K 的可行下降方向。 现在，我们给出用基本迭
11、代格式（ 1）求解 (NP)的一般步骤如下： 0 选取初始点0 x ，令 0:=k 。 1 构造搜索方向，依照一定规划，构造 f 在点kx 处关于 K 的可行下降方向作为搜索方向kp 。 2 寻求搜索步长。以kx 为起点沿搜索方向kp 寻求适当的步长kt ，使目标函数值有某种意义的下降。 3 求出下一个迭代点。按迭代格式（ 1）求出 kkkkptxx +=+1。 若1+kx 已满足某种终止条件，停止迭代。 4 以1+kx 代替kx ，回到 1步。 1.5 凸函数、凸规划 设 )(xf 为定义在 n 维欧氏空间)(nE 中某个凸集 R 上的函数，若对任何实数)10( 以及 R 中的任意两点)1(
12、x 和)2(x ，恒有 )()1()()1()2()1()2()1(xfxfxxf + 则称 )(xf 为定义在 R 上的凸函数。 若对每一个 )10( 和 Rxx )2()1(恒有 )()1()()1()2()1()2()1(xfxfxxf + 则称 )(xf 为定义在 R 上的严格凸函数。 考虑非线性规划 =,2,1,0)(|)( minljxgxRxfjRxL 假定其中 )(xf 为凸函数， ),2,1)( ljxgjL= 为凸函数，这样的非线性规划称为凸规划。 可以证明，凸规划的可行域为凸集，其局部最优解即为全局最优解，而且其最优解的集合形成一个凸集。当凸规划的目标函数 )(xf 为严
13、格凸函数时，其最优解必定唯一（假定最优解存在） 。由此可见，凸规划是一类比较简单而又具有重要理论意义的非 -36- 线性规划。 2 无约束问题 2.1 一维搜索方法 当用迭代法求函数的极小点时，常常用到一维搜索，即沿某一已知方向求目标函数的极小点。一维搜索的方法很多，常用的有： （ 1）试探法（ “成功失败” ，斐波那契法，0.618 法等） ； （ 2）插值法（抛物线插值法，三次插值法等） ； （ 3）微积分中的求根法（切线法，二分法等） 。 考虑一维极小化问题 )(min tfbta （ 2） 若 )(tf 是 , ba 区间上的下单峰函数，我们介绍通过不断地缩短 , ba 的长度，来搜索
14、得（ 2）的近似最优解的两个方法。 为了缩短区间 , ba ，逐步搜索得（ 2）的最优解*t 的近似值，我们可以采用以下途径：在 , ba 中任取两个关于 , ba 是对称的点1t 和2t （不妨设12tt ，并把它们叫做搜索点） ，计算 )(1tf 和 )(2tf 并比较它们的大小。对于单峰函数，若 )()(12tftf ，则必有 ,1*tat ，因而 ,1ta 是缩短了的单峰区间；若 )()(21tftf ，则有,2*btt ，故 ,2bt 是缩短了的单峰区间；若 )()(12tftf = ，则 ,1ta 和 ,2bt 都是缩短了的单峰。因此通过两个搜索点处目标函数值大小的比较，总可以获得
15、缩短了的单峰区间。对于新的单峰区间重复上述做法，显然又可获得更短的单峰区间。如此进行，在单峰区间缩短到充分小时，我们可以取最后的搜索点作为（ 2）最优解的近似值。 应该按照怎样的规则来选取探索点，使给定的单峰区间的长度能尽快地缩短？ 2.1.1 Fibonacci 法 如用nF 表示计算 n 个函数值能缩短为单位长区间的最大原区间长度， 可推出nF 满足关系 110= FF ,3,2,12L=+=nFFFnnn 数列 nF 称为 Fibonacci 数列，nF 称为第 n 个 Fibonacci 数，相邻两个 Fibonacci 数之比nnFF1称为 Fibonacci 分数。 当用斐波那契法
16、以 n 个探索点来缩短某一区间时，区间长度的第一次缩短率为nnFF1，其后各次分别为212312,FFFFFFnnnnL。由此，若1t 和 )(122ttt ，这就要求最后区间的长度不超过 ，即 nFab （ 4） 据此，我们应按照预先给定的精度 ，确定使（ 4）成立的最小整数 n 作为搜索次数，直到进行到第 n 个探索点时停止。 用上述不断缩短函数 )(tf 的单峰区间 , ba 的办法，来求得问题（ 2）的近似解，是 Kiefer(1953 年 )提出的，叫做 Finbonacci 法，具体步骤如下： o1 选取初始数据，确定单峰区间 ,00ba ，给出搜索精度 0 ，由（ 4）确定搜索次
17、数 n 。 o2 00,1 bbaak = ，计算最初两个搜索点，按（ 3）计算1t 和2t 。 o3 while 1 nk )(),(2211tfftff = if 21ff ，满足比例关系 -38- =11 称之为黄金分割数，其值为 L6180339887.0215= 。 黄金分割数 和 Fibonacci 分数之间有着重要的关系 nnnFF1lim= 。 现用不变的区间缩短率 0.618，代替斐波那契法每次不同的缩短率，就得到了黄金分割法（ 0.618 法） 。这个方法可以看成是斐波那契法的近似，实现起来比较容易，效果也相当好，因而易于为人们所接受。 用 0.618 法求解，从第 2 个
18、探索点开始每增加一个探索点作一轮迭代以后，原单峰区间要缩短 0.618 倍。计算 n 个探索点的函数值可以把原区间 ,00ba 连续缩短 1n次，因为每次的缩短率均为 ，故最后的区间长度为 100)(nab 这就是说，当已知缩短的相对精度为 时，可用下式计算探索点个数 n ： 1n 当然，也可以不预先计算探索点的数目 n ，而在计算过程中逐次加以判断，看是否已满足了提出的精度要求。 0.618 法是一种等速对称进行试探的方法，每次的探索点均取在区间长度的 0.618倍和 0.382 倍处。 2.2 二次插值法 对极小化问题（ 2） ，当 )(tf 在 , ba 上连续时，可以考虑用多项式插值来
19、进行一维搜索。它的基本思想是：在搜索区间中，不断用低次（通常不超过三次）多项式来近似目标函数，并逐步用插值多项式的极小点来逼近（ 2）的最优解。 2.3 无约束极值问题的解法 无约束极值问题可表述为 )(),( minnExxf （ 5） 求解问题（ 5）的迭代法大体上分为两点：一是用到函数的一阶导数或二阶导数，称为解析法。另一是仅用到函数值，称为直接法。 2.3.1 解析法 2.3.1.1 梯度法（最速下降法） 对基本迭代格式 kkkkptxx +=+1 （ 6） 我们总是考虑从点kx 出发沿哪一个方向kp ，使目标函数 f 下降得最快。微积分的知识告诉我们，点kx 的负梯度方向 )(kkx
20、fp = ， 是从点kx 出发使 f 下降最快的方向。为此，称负梯度方向 )(kxf 为 f 在点kx 处的最速下降方向。 -39- 按基本迭代格式（ 6） ，每一轮从点kx 出发沿最速下降方向 )(kxf 作一维搜索，来建立求解无约束极值问题的方法，称之为最速下降法。 这个方法的特点是，每轮的搜索方向都是目标函数在当前点下降最快的方向。同时，用 0)( =kxf 或 )(kxf 作为停止条件。其具体步骤如下： 1选取初始数据。选取初始点0 x ，给定终止误差，令 0:=k 。 2求梯度向量。计算 )(kxf ， 若 )(kxf ，停止迭代，输出kx 。否则，进行 3。 3 构造负梯度方向。取
21、 )(kkxfp = . 4 进行一维搜索。求kt ，使得 )(min)(0kktkkktpxfptxf +=+ 令 ,1:,1+=+=+kkptxxkkkk转 2。 例 4 用最速下降法求解无约束非线性规划问题 222125)(min xxxf += 其中Txxx ),(21= ，要求选取初始点Tx )2,2(0= 。 解 （ i）Txxxf )50,2()(21= 编写 M 文件 detaf.m，定义函数 )(xf 及其梯度列向量如下 function f,df=detaf(x); f=x(1)2+25*x(2)2; df=2*x(1) 50*x(2); （ii）编写主程序文件zuisu.
22、m如下： clc x=2;2; f0,g=detaf(x); while norm(g)0.000001 p=-g/norm(g); t=1.0;f=detaf(x+t*p); while ff0 t=t/2; f=detaf(x+t*p); end x=x+t*p; f0,g=detaf(x); end x,f0 2.3.1.2 Newton 法 考虑目标函数 f 在点kx 处的二次逼近式 )()(21 )()()()()(2 kkTkkTkkxxxfxxxxxfxfxQxf += 假定 Hesse 阵 -40- =2212122122)()()()()(nknknkkkxxfxxxfxxx
23、fxxfxfLMLML 正定。 由于 )(2 kxf 正定，函数 Q 的驻点1+kx 是 )(xQ 的极小点。为求此极小点，令 0)()()(121=+=+ kkkkkxxxfxfxQ ， 即可解得 )()(121 kkkkxfxfxx =+. 对照基本迭代格式（ 1） ，可知从点kx 出发沿搜索方向。 )()(12 kkkxfxfp = 并取步长 1=kt 即可得 )(xQ 的最小点1+kx 。通常，把方向kp 叫做从点kx 出发的Newton 方向。从一初始点开始，每一轮从当前迭代点出发，沿 Newton 方向并取步长为 1 的求解方法，称之为 Newton 法。其具体步骤如下： 1选取初
24、始数据。选取初始点0 x ，给定终止误差 0 ，令 0:=k 。 2求梯度向量。计算 )(kxf ，若 )(kxf ，停止迭代，输出kx 。否则，进行 3。 3构造 Newton 方向。计算12)(kxf ，取 )()(12 kkkxfxfp =. 4 求下一迭代点。令 1:,1+=+=+kkpxxkkk，转 2。 例 5 用 Newton 法求解， 2221424125)(min xxxxxf += 选取Tx )2,2(0= 。 解： （ i）Txxxxxxxf 210024)(2213222131+= +=2122212122212230044212xxxxxxxxf （ii）编写 M 文
25、件 nwfun.m 如下： function f,df,d2f=nwfun(x); f=x(1)4+25*x(2)4+x(1)2*x(2)2; df=4*x(1)3+2*x(1)*x(2)2;100*x(2)3+2*x(1)2*x(2); d2f=2*x(1)2+2*x(2)2,4*x(1)*x(2) 4*x(1)*x(2),300*x(2)2+2*x(1)2; （III）编写主程序文件 example5.m 如下： clc x=2;2; f0,g1,g2=nwfun(x); while norm(g1)0.00001 p=-inv(g2)*g1; x=x+p; -41- f0,g1,g2=n
26、wfun(x); end x, f0 如果目标函数是非二次函数，一般地说，用 Newton 法通过有限轮迭代并不能保证可求得其最优解。 为了提高计算精度，我们在迭代时可以采用变步长计算上述问题，编写主程序文件example5_2如下： clc,clear x=2;2; f0,g1,g2=nwfun(x); while norm(g1)0.00001 p=-inv(g2)*g1;p=p/norm(p); t=1.0;f=nwfun(x+t*p); while ff0 t=t/2;f=nwfun(x+t*p); end x=x+t*p; f0,g1,g2=nwfun(x); end x,f0 Ne
27、wton 法的优点是收敛速度快；缺点是有时不好用而需采取改进措施，此外，当维数较高时，计算12)(kxf 的工作量很大。 2.3.1.3 变尺度法 变尺度法（ Variable Metric Algorithm）是近 20 多年来发展起来的，它不仅是求解无约束极值问题非常有效的算法，而且也已被推广用来求解约束极值问题。由于它既避免了计算二阶导数矩阵及其求逆过程，又比梯度法的收敛速度快，特别是对高维问题具有显著的优越性，因而使变尺度法获得了很高的声誉。下面我们就来简要地介绍一种变尺度法 DFP 法的基本原理及其计算过程。这一方法首先由 Davidon 在 1959 年提出，后经 Fletcher
28、 和 Powell 加以改进。 我们已经知道，牛顿法的搜索方向是 )()(12 kkxfxf ，为了不计算二阶导数矩阵 )(2 kxf 及其逆阵，我们设法构造另一个矩阵，用它来逼近二阶导数矩阵的逆阵12)(kxf ，这一类方法也称拟牛顿法（ Quasi-Newton Method） 。 下面研究如何构造这样的近似矩阵，并将它记为)(kH 。我们要求：每一步都能以现有的信息来确定下一个搜索方向；每做一次选代，目标函数值均有所下降；这些近似矩阵最后应收敛于解点处的 Hesse 阵的逆阵。 当 )(xf 是二次函数时， 其 Hesse 阵为常数阵 A， 任两点kx 和1+kx 处的梯度之差为 )()
29、()(11 kkkkxxAxfxf =+ 或 )()(111 kkkkxfxfAxx =+ 对于非二次函数，仿照二次函数的情形，要求其 Hesse 阵的逆阵的第 1+k 次近似矩阵)1( +kH 满足关系式 )()(1)1(1 kkkkkxfxfHxx =+ （ 7） -42- 这就是常说的拟 Newton 条件。 若令 =+kkkkkkxxxxfxfG11)()()( （ 8） 则式（ 7）变为 )()1( kkkGHx =+， （ 9） 现假定)(kH 已知，用下式求)1( +kH （设)(kH 和)1( +kH 均为对称正定阵） ； )()()1( kkkHHH +=+ （ 10） 其中
30、)(kH 称为第 k 次校正矩阵。显然，)1( +kH 应满足拟 Newton 条件（ 9） ，即要求 )()()()(kkkkGHHx += 或 )()()()( kkkkkGHxGH = （ 11） 由此可以设想， )(kH 的一种比较简单的形式是 TkkkTkkkWGHQxH )()()()()()()(= (12) 其中)(kQ 和)(kW 为两个待定列向量。 将式（ 12）中的)(kH 代入（ 11） ，得 )()()()()()()()()()(kkkkTkkkkTkkGHxGWGHGQx = 这说明，应使 1)()()()()()(=kTkkTkGWGQ （ 13） 考虑到)(k
31、H 应为对称阵，最简单的办法就是取 =)()()()(kkkkkkkGHWxQ （ 14） 由式（ 13）得 1)()()()()()(=kkTkkkTkkGHGGx （ 15） 若)()(kTkGx 和)()()()(kkTkGHG 不等于零，则有 =)()()()()()(1)(1)(1kkTkkkTkkTkkGHGxGGx （ 16） 于是，得校正矩阵 )()()()()()()()()()()()()(kkTkkTkkkkTkTkkkGHGHGGHxGxxH= (17) 从而得到 )()()()()()()()()()1()()()()(kkTkkTkkkkTkTkkkkGHGHGGH
32、xGxxHH+=+ （ 18） 上述矩阵称为尺度矩阵。通常，我们取第一个尺度矩阵)0(H 为单位阵，以后的尺度矩 -43- 阵按式（ 18）逐步形成。可以证明： （ i）当kx 不是极小点且)(kH 正定时，式（ 17）右端两项的分母不为零，从而可按式（ 18）产生下一个尺度矩阵)1( +kH ； （ ii）若)(kH 为对称正定阵，则由式（ 18）产生的)1( +kH 也是对称正定阵； （ iii）由此推出 DFP 法的搜索方向为下降方向。 现将 DFP 变尺度法的计算步骤总结如下。 1给定初始点0 x 及梯度允许误差 0 。 2若 )(0 xf ，则0 x 即为近似极小点，停止迭代，否则，
33、转向下一步。 3令 IH =)0(（单位矩阵） ， )(0)0(0 xfHp = 在0p 方向进行一维搜索，确定最佳步长0 ： )()(min00000pxfpxf +=+ 如此可得下一个近似点 0001pxx += 4一般地，设已得到近似点kx ，算出 )(kxf ，若 )(0 xf 则kx 即为所求的近似解，停止迭代；否则，计算)(kH ： )1()1()1()1()1()1()1(1)1(11)1()()()( )()(+=kkTkkTkkkkTkTkkkkGHGHGGHxGxxHH 并令 )()( kkkxfHp = ，在kp 方向上进行一维搜索，得k ，从而可得下一个近似点 kkkk
34、pxx +=+1 5若1+kx 满足精度要求，则1+kx 即为所求的近似解，否则，转回 4，直到求出某点满足精度要求为止。 2.3.2 直接法 在无约束非线性规划方法中，遇到问题的目标函数不可导或导函数的解析式难以表示时，人们一般需要使用直接搜索方法。同时，由于这些方法一般都比较直观和易于理解，因而在实际应用中常为人们所采用。下面我们介绍 Powell 方法。 这个方法主要由所谓基本搜索、加速搜索和调整搜索方向三部分组成，具体步骤如下： 1 选取初始数据。选取初始点0 x ， n 个线性无关初始方向，组成初搜索方向组,110 nppp L 。给定终止误差 0 ，令 0:=k 。 2进行基本搜索
35、。令kxy =:0，依次沿 ,110 nppp L 中的方向进行一维搜索。对应地得到辅助迭代点nyyy ,21L ，即 111 +=jjjjptyy njtpyfptyfjjtjjj,1),(min)(110111L=+=+ -44- 3构造加速方向。令0yypnn= ，若 np ，停止迭代，输出nkyx =+1。否则进行 4。 4确定调整方向。按下式 1|)()(max)()(11njyfyfyfyfjjmm= 找出 m 。若 )()(2)2()(2)(100 mmnnyfyfyyfyfyf M ，构造函数 =+=tiisiiriixkMxhMxgMxfMxP111|)(|)0),(min(
36、)0),(max()(),( （或 |)(|0)(minsum0)(maxsum)(),( xKMxHMxGMxfMxP += 这里 )()()(1xgxgxGrL= ， )()()(1xhxhxHsL= ， )()()(1tkxkxKtL= ， Matlab 中可以直接利用 max 、 min 和 sum 函数。 ）则以增广目标函数 ),( MxP 为目标函数的无约束极值问题 ),(min MxP 的最优解 x 也是原问题的最优解。 例 9 求下列非线性规划 =+=.0, 02 0 8)(min 212212212221xxxxxxxxxf 解 (i)编写 M 文件 test.m funct
37、ion g=test(x); M=50000; f=x(1)2+x(2)2+8; g=f-M*min(x(1),0)-M*min(x(2),0)-M*min(x(1)2-x(2),0)+. -47- M*abs(-x(1)-x(2)2+2); 或者是利用 Matlab的求矩阵的极小值和极大值函数编写 test.m如下： function g=test(x); M=50000; f=x(1)2+x(2)2+8; g=f-M*sum(min(x;zeros(1,2)-M*min(x(1)2-x(2),0)+. M*abs(-x(1)-x(2)2+2); 我们也可以修改罚函数的定义，编写 test.
38、m如下： function g=test(x); M=50000; f=x(1)2+x(2)2+8; g=f-M*min(min(x),0)-M*min(x(1)2-x(2),0)+M*(-x(1)-x(2)2+2)2; (ii)在 Matlab 命令窗口输入 x,y=fminunc(test,rand(2,1) 即可求得问题的解。 3.3 Matlab 求约束极值问题 在 Matlab 优化工具箱中，用于求解约束最优化问题的函数有： fminbnd、 fmincon、quadprog、 fseminf、 fminimax，上面我们已经介绍了函数 fmincon 和 quadprog。 3.3
39、.1 fminbnd 函数 求单变量非线性函数在区间上的极小值 , , )(min21xxxxfx Matlab 的命令为 X,FVAL = FMINBND(FUN,x1,x2,OPTIONS)， 它的返回值是极小点 x 和函数的极小值。这里 fun 是用 M 文件定义的函数或 Matlab 中的单变量数学函数。 例 10 求函数 5,0 ,1)3()(2= xxxf 的最小值。 解 编写 M 文件 fun5.m function f=fun5(x); f=(x-3)2-1; 在 Matlab 的命令窗口输入 x,y=fminbnd(fun5,0,5) 即可求得极小点和极小值。 3.3.2 f
40、seminf 函数 求 0),(,0)(,0)(|)(min = wxPHIxCeqxCxFx =BeqxAeqBxA* s.t. 其中 ),(),(),( wxPHIxCeqxC 都是向量函数； w是附加的向量变量， w的每个分量都限定在某个区间内。 上述问题的 Matlab 命令格式为 X=FSEMINF(FUN,X0,NTHETA,SEMINFCON,A,B,Aeq,Beq) 其中 FUN 用于定义目标函数 )(xF ； X0 为 x 的初始值； NTHETA 是半无穷约束),( wxPHI 的个数；函数 SEMINFCON 用于定义非线性不等式约束 )(xC ，非线性等 -48- 式约
41、束 )(xCeq 和半无穷约束 ),( wxPHI 的每一个分量函数， 函数 SEMINFCON 有两个输入参量 X 和 S， S 是推荐的取样步长，也许不被使用。 例 11 求函数232221)5.0()5.0()5.0()( += xxxxf 取最小值时的 x 值 ,约束为： 1)sin()50(10001)cos()sin(),(33121211111= xxwwxwxwwxK 1)sin()50(10001)cos()sin(),(33222122222= xxwwxwxwwxK 10011 w ， 10012 w 解 （ 1）编写 M 文件 fun6.m 定义目标函数如下： func
42、tion f=fun6(x,s); f=sum(x-0.5).2); （ 2）编写 M 文件 fun7.m 定义约束条件如下： function c,ceq,k1,k2,s=fun7(x,s); c=;ceq=; if isnan(s(1,1) s=0.2,0;0.2 0; end %取样值 w1=1:s(1,1):100; w2=1:s(2,1):100; %半无穷约束 k1=sin(w1*x(1).*cos(w1*x(2)-1/1000*(w1-50).2-sin(w1*x(3)-x(3)-1; k2=sin(w2*x(2).*cos(w2*x(1)-1/1000*(w2-50).2-si
43、n(w2*x(3)-x(3)-1; %画出半无穷约束的图形 plot(w1,k1,-,w2,k2,+); （ 3）调用函数 fseminf 在 Matlab 的命令窗口输入 x,y=fseminf(fun6,rand(3,1),2,fun7) 即可。 3.3.3 fminimax 函数 求解)(maxmin xFiFx =UBxLBxCeqxCBeqxAeqbxA0)(0)(* .s.t 其中 )(,),()(1xFxFxFmL= 。 上述问题的 Matlab 命令为 X=FMINIMAX(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON) -49- 例 12 求函数族 )()
44、,(),(),(),(54321xfxfxfxfxf 取极大极小值时的 x 值。 其中： +=+=+=8)()(183)(3)(30440482)(215214213222122122211xxxfxxxfxxxfxxxfxxxxxf 解 （ 1）编写 M 文件 fun8.m 定义向量函数如下： function f=fun8(x); f=2*x(1)2+x(2)2-48*x(1)-40*x(2)+304 -x(1)2-3*x(2)2 x(1)+3*x(2)-18 -x(1)-x(2) x(1)+x(2)-8; （ 2）调用函数 fminimax x,y=fminimax(fun8,rand(
45、2,1) 3.3.4 利用梯度求解约束优化问题 例 13 已知函数 )12424()(22122211+= xxxxxexfx，且满足非线性约束： 105.1212121xxxxxx 求 )(min xfx 分析：当使用梯度求解上述问题时，效率更高并且结果更准确。 题目中目标函数的梯度为： +)244()168424(212121222111xxexxxxxxexx 解 （ 1）编写 M 文件 fun9.m 定义目标函数及梯度函数： function f,df=fun9(x); f=exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1); df=exp(
46、x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+8*x(1)+6*x(2)+1);exp(x(1)*(4*x(2)+4*x(1)+2); （ 2）编写 M 文件 fun10.m 定义约束条件及约束条件的梯度函数： function c,ceq,dc,dceq=fun10(x); c=x(1)*x(2)-x(1)-x(2)+1.5;-x(1)*x(2)-10; dc=x(2)-1,-x(2);x(1)-1,-x(1); ceq=;dceq=; （ 3）调用函数 fmincon，编写主函数文件 example13.m 如下： %采用标准算法 options=optimset(l
47、argescale,off); %采用梯度 options=optimset(options,GradObj,on,GradConstr,on); x,y=fmincon(fun9,rand(2,1),fun10,options) -50- 3.4 Matlab 优化工具箱的用户图形界面解法 Matlab 优化工具箱中的 optimtool 命令提供了优化问题的用户图形界面解法。optimtool 可应用到所有优化问题的求解，计算结果可以输出到 Matlab 工作空间中。 图 1 优化问题用户图形界面解法示意图 例 14 用 optimtool 重新求解例 1。 利用例 1 已经定义好的函数
48、fun1 和 fun2。在 Matlab 命令窗口运行 optimtool，就打开图形界面，如图 1 所示，填入有关的参数，未填入的参数取值为空或者为默认值，然后用鼠标点一下“ start”按钮，就得到求解结果，再使用“ file”菜单下的“ Export to Workspace”选项，把计算结果输出到 Matlab 工作空间中去。 4 飞行管理问题 在约 10， 000m 高空的某边长 160km 的正方形区域内，经常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区
49、域内的飞机发生碰撞。如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞。现假定条件如下： 1）不碰撞的标准为任意两架飞机的距离大于 8km; 2）飞机飞行方向角调整的幅度不应超过 30 度； 3）所有飞机飞行速度均为每小时 800km; 4）进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在 60km 以上； 5）最多需考虑 6 架飞机； 6）不必考虑飞机离开此区域后的状况。 请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过 0.01 度） ，要求飞机飞行方向角调整的幅度尽量小。 设该区域 4 个顶点的座标为 (0
50、,0)， (160,0)， (160,160)， (0,160)。记录数据见表 1。 -51- 表 1 飞行记录数据 飞机编号 横座标 x 纵座标 y 方向角（度） 1 150 140 243 2 85 85 236 3 150 155 220.5 4 145 50 159 5 130 150 230 新进入 0 0 52 注：方向角指飞行方向与 x 轴正向的夹角。 为方便以后的讨论，我们引进如下记号： D 为飞行管理区域的边长； 为飞行管理区域，取直角坐标系使其为 ,0,0 DD ； a 为飞机飞行速度， 800=a km/h； ),(00iiyx 为第 i 架飞机的初始位置； )(),(
