企业官网建设流程全解析

网站建设分金手指排名十四,头皮痒 脱发严重怎么办,商丘网站开发公司,游戏开发团队第三章 非线性规划1 非线性规划 1.1 非线性规划的实例与定义 如果目标函数或约束条件中包含非线性函数#xff0c;就称这种规划问题为非线性规划问题。一般说来#xff0c;解非线性规划要比解线性规划问题困难得多。而且#xff0c; 也不象线性规划有单纯形法这一通…第三章 非线性规划§ 1 非线性规划1.1 非线性规划的实例与定义如果目标函数或约束条件中包含非线性函数就称这种规划问题为非线性规划问题。一般说来解非线性规划要比解线性规划问题困难得多。而且 也不象线性规划有单纯形法这一通用方法非线性规划目前还没有适于各种问题的一般算法各个方法都有自己特定的适用范围。下面通过实例归纳出非线性规划数学模型的一般形式介绍有关非线性规划的基本概念。例 1 投资决策问题某企业有n 个项目可供选择投资并且至少要对其中一个项目投资。已知该企业拥有总资金 A 元投资于第i(i 1, … , n) 个项目需花资金ai 元并预计可收益bi 元。试选择最佳投资方案。解 设投资决策变量为xi 目i 1, … , n 则投资总额为 ai xi 投资总收益为 bi xi 。因为该公司至少要对一个项目投资并且总的投资金额不能超过总资金 A 故有限制条件另外由于xi (i 1, … , n) 只取值 0 或 1所以还有xi (1− xi ) 0, i 1, … , n.最佳投资方案应是投资额最小而总收益最大的方案所以这个最佳投资决策问题归结为总资金以及决策变量取 0 或 1的限制条件下极大化总收益和总投资之比。因此其数学模型为s.t. ai xi ≤ Axi (1− xi ) 0, i 1, … , n.上面例题是在一组等式或不等式的约束下求一个函数的最大值或最小值问题其中至少有一个非线性函数这类问题称之为非线性规划问题。可概括为一般形式min f (x)s.t. hj (x) ≤ 0, j 1, … , q (NP)gi (x) 0, i 1, … , p-32-其中x [x1 … xn ]T 称为模型NP的决策变量f 称为目标函数gi (i 1, … , p)和 hj (j 1, … , q) 称为约束函数。另外 gi (x) 0 (i 1, … , p) 称为等式约束 hj (x) ≤ 0 (j 1, … , q) 称为不等式的约束。对于一个实际问题在把它归结成非线性规划问题时一般要注意如下几点i确定供选方案首先要收集同问题有关的资料和数据在全面熟悉问题的基础上确认什么是问题的可供选择的方案并用一组变量来表示它们。ii提出追求目标经过资料分析根据实际需要和可能提出要追求极小化或极大化的目标。并且运用各种科学和技术原理把它表示成数学关系式。iii给出价值标准在提出要追求的目标之后要确立所考虑目标的“好”或“坏”的价值标准并用某种数量形式来描述它。iv寻求限制条件由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果因此还需要寻找出问题的所有限制条件这些条件通常用变量之间的一些不等式或等式来表示。1.2 线性规划与非线性规划的区别如果线性规划的最优解存在其最优解只能在其可行域的边界上达到特别是可行域的顶点上达到而非线性规划的最优解如果最优解存在则可能在其可行域的任意一点达到。1.3 非线性规划的 Matlab 解法Matlab 中非线性规划的数学模型写成以下形式min f (x)其中f (x) 是标量函数A, B, Aeq, Beq 是相应维数的矩阵和向量C(x), Ceq(x) 是非线性向量函数。Matlab 中的命令是XFMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)它的返回值是向量x 其中 FUN 是用 M 文件定义的函数f(x) X0 是x 的初始值 A,B,Aeq,Beq 定义了线性约束 A * X ≤ B, Aeq * X Beq 如果没有线性约束 则A[],B[],Aeq[],Beq[]LB 和 UB 是变量x 的下界和上界如果上界和下界没有约束则 LB[]UB[]如果x 无下界则 LB 的各分量都为-inf如果x 无上界则 UB的各分量都为 infNONLCON 是用 M 文件定义的非线性向量函数 C(x), Ceq(x) OPTIONS定义了优化参数可以使用 Matlab 缺省的参数设置。例2 求下列非线性规划min f (x) x12 x22 x 8s.t. x − x2 x ≥ 0x1 x x ≤ 20− x1 − x22 2 0-33-x2 2x 3x1 , x2 , x3 ≥ 0解 i编写 M 文件 fun1.m 定义目标函数function ffun1(x);fsum(x.^2)8;ii编写M文件fun2.m定义非线性约束条件function [g,h]fun2(x);g[-x(1)2x(2)-x(3)2x(1)x(2)2x(3)3-20]; %非线性不等式约束h[-x(1)-x(2)^22x(2)2x(3)^2-3]; %非线性等式约束iii编写主程序文件 example2.m 如下optionsoptimset(‘largescale’,‘off’);[x,y]fmincon(‘fun1’,rand(3,1),[],[],[],[],zeros(3,1),[], … ‘fun2’, options)就可以求得当x1 0.5522, x2 1.2033, x3 0.9478 时最小值y 10.6511。1.4 求解非线性规划的基本迭代格式记NP的可行域为K 。若x * ∈ K 并且f (x * ) ≤ f (x), ∀x ∈ K则称x是NP的整体最优解f (x * ) 是(NP)的整体最优值。如果有f (x * ) f (x), ∀x ∈ K, x ≠ x *则称x* 是NP的严格整体最优解f (x * ) 是(NP)的严格整体最优值。若x * ∈ K 并且存在 x* 的邻域Nδ (x * ) 使f (x* ) ≤ f (x), ∀x ∈ Nδ (x* ) ∩ K 则称x* 是NP的局部最优解f (x * ) 是(NP)的局部最优值。如果有f (x* ) f (x), ∀x ∈ Nδ (x* ) ∩ K则称x* 是NP的严格局部最优解f (x * ) 是(NP)的严格局部最优值。由于线性规划的目标函数为线性函数可行域为凸集因而求出的最优解就是整个可行域上的全局最优解。非线性规划却不然 有时求出的某个解虽是一部分可行域上的极值点但却并不一定是整个可行域上的全局最优解。对于非线性规划模型(NP)可以采用迭代方法求它的最优解。迭代方法的基本思想是从一个选定的初始点 x0 ∈ R n 出发按照某一特定的迭代规则产生一个点列{xk } 使得当{xk } 是有穷点列时其最后一个点是(NP)的最优解当 {xk } 是无穷点列时它有极限点并且其极限点是(NP)的最优解。设x k ∈ Rn 是某迭代方法的第k 轮迭代点x k 1 ∈ Rn 是第k 1 轮迭代点记x k 1 x k tk pk 1这里tk ∈ R1 , pk ∈ R n , pk 1 并且pk 的方向是从点x k 向着点xk 1 的方向。式1就是求解非线性规划模型(NP)的基本迭代格式。-34-通常我们把基本迭代格式1中的pk 称为第k 轮搜索方向tk 为沿pk 方向的步长使用迭代方法求解(NP)的关键在于如何构造每一轮的搜索方向和确定适当的步长。设x ∈ R n , p ≠ 0 若存在δ 0 使f (x tp) f (x ), ∀t ∈(0, δ) 称向量p 是f 在点x 处的下降方向。设x ∈ R n , p ≠ 0 若存在 t 0 使x tp ∈ K 称向量p 是点x 处关于K 的可行方向。一个向量p 若既是函数f 在点x 处的下降方向又是该点关于区域K 的可行方向称之为函数f 在点x 处关于K 的可行下降方向。现在我们给出用基本迭代格式1求解(NP)的一般步骤如下0° 选取初始点x0 令 k : 0 。1° 构造搜索方向依照一定规划构造f 在点 x k 处关于K 的可行下降方向作为搜索方向pk 。2° 寻求搜索步长。以xk 为起点沿搜索方向pk 寻求适当的步长tk 使目标函数值有某种意义的下降。3° 求出下一个迭代点。按迭代格式1求出x k 1 x k tk pk 。若x k 1 已满足某种终止条件停止迭代。4° 以xk 1 代替x k 回到 1°步。1.5 凸函数、凸规划设 f(x) 为定义在 n 维欧氏空间 E(n ) 中某个凸集 R 上的函数若对任何实数α(0 α 1) 以及R 中的任意两点x( 1) 和x(2) 恒有f ( αx( 1) (1 − α)x( 2) ) ≤ αf (x( 1) ) (1 − α)f (x( 2) )则称f(x) 为定义在R 上的凸函数。若对每一个 α(0 α 1) 和 x( 1) ≠ x( 2) ∈ R 恒有f ( αx( 1) (1 − α)x( 2) ) αf (x( 1) ) (1 − α)f (x( 2) )则称f(x) 为定义在R 上的严格凸函数。考虑非线性规划假定其中f (x) 为凸函数gj (x)(j 1,2, … , l) 为凸函数这样的非线性规划称为凸规划。可以证明凸规划的可行域为凸集其局部最优解即为全局最优解而且其最优解的集合形成一个凸集。当凸规划的目标函数f(x) 为严格凸函数时其最优解必定唯一假定最优解存在。由此可见 凸规划是一类比较简单而又具有重要理论意义的非-35-线性规划。§2 无约束问题2.1 一维搜索方法当用迭代法求函数的极小点时常常用到一维搜索即沿某一已知方向求目标函数的极小点。一维搜索的方法很多常用的有1试探法“成功—失败”斐波那契法 0.618 法等2插值法抛物线插值法三次插值法等3微积分中的求根法切线法二分法等。考虑一维极小化问题min f (t) 2若f (t) 是 [a, b] 区间上的下单峰函数我们介绍通过不断地缩短[a, b] 的长度来搜索得2的近似最优解的两个方法。为了缩短区间 [a, b] 逐步搜索得2的最优解 t* 的近似值我们可以采用以下途径在 [a, b] 中任取两个关于 [a, b] 是对称的点 t1 和 t2 不妨设 t2 t1 并把它们叫做搜索点计算f (t1 ) 和f (t2 ) 并比较它们的大小。对于单峰函数 若f (t2 ) f (t1 ) 则必有 t* ∈ [a, t1 ] 因 而 [a, t1 ] 是缩短了 的单峰区间 若 f (t1 ) f (t2 ) 则有t* ∈ [t2, b] 故 [t2, b] 是缩短了的单峰区间若f (t2 ) f (t1 ) 则 [a, t1 ] 和 [t2, b] 都是缩短了的单峰。因此通过两个搜索点处目标函数值大小的比较 总可以获得缩短了的单峰区间。对于新的单峰区间重复上述做法显然又可获得更短的单峰区间。如此进行在单峰区间缩短到充分小时我们可以取最后的搜索点作为2最优解的近似值。应该按照怎样的规则来选取探索点使给定的单峰区间的长度能尽快地缩短2.1.1 Fibonacci 法如用Fn 表示计算n 个函数值能缩短为单位长区间的最大原区间长度可推出Fn 满足关系F0 F1 1Fn Fn −2 Fn −1 , n 2,3, … ,数列{Fn } 称为 Fibonacci 数列Fn 称为第n 个 Fibonacci 数相邻两个 Fibonacci 数之比 称为 Fibonacci 分数。当用斐波那契法以 n 个探索点来缩短某一区间时区间长度的第一次缩短率为 , 其后各次分别为 。由此 若t1 和t2 (t2 t1 ) 是单峰区间 [a, b]中第 1 个和第 2 个探索点的话那么应有比例关系从而它们关于 [a, b] 确是对称的点。-36-如果要求经过一系列探索点搜索之后使最后的探索点和最优解之间的距离不超过精度δ 0 这就要求最后区间的长度不超过δ ,即4据此我们应按照预先给定的精度 δ ,确定使4成立的最小整数n 作为搜索次数直到进行到第n 个探索点时停止。用上述不断缩短函数f(t) 的单峰区间 [a, b] 的办法来求得问题2的近似解是 Kiefer(1953 年)提出的叫做 Finbonacci 法具体步骤如下1o 选取初始数据确定单峰区间 [a0, b0 ] 给出搜索精度 δ 0 由4确定搜索次数n 。2o k 1, a a0, b b0 计算最初两个搜索点按3计算t1 和 t2 。3o while k n −1f 1 f (t1 ), f2 f (t2 )if f1 f2elseendk k 1end4o 当进行至k n −1时这就无法借比较函数值f (t1 ) 和f(t2 ) 的大小确定最终区间为此取其中ε 为任意小的数。在t1 和t2 这两点中以函数值较小者为近似极小点相应的函数值为近似极小值。并得最终区间 [a, t1 ] 或 [t2, b] 。由上述分析可知斐波那契法使用对称搜索的方法逐步缩短所考察的区间它能以尽量少的函数求值次数达到预定的某一缩短率。例 3 试用斐波那契法求函数f(t) t2 − t 2 的近似极小点要求缩短后的区间不大于区间[−1,3] 的 0.08 倍。程序留作习题。2.1.2 0.618 法若 ω 0 满足比例关系-37-称之为黄金分割数其值为 黄金分割数 ω 和 Fibonacci 分数之间有着重要的关系现用不变的区间缩短率 0.618代替斐波那契法每次不同的缩短率就得到了黄金分割法0.618 法。这个方法可以看成是斐波那契法的近似实现起来比较容易效果也相当好因而易于为人们所接受。用 0.618 法求解从第 2 个探索点开始每增加一个探索点作一轮迭代以后原单峰区间要缩短 0.618 倍。计算n 个探索点的函数值可以把原区间 [a0, b0 ] 连续缩短n −1次因为每次的缩短率均为 μ ,故最后的区间长度为(b0 − a0 ) μ n −1这就是说当已知缩短的相对精度为 δ 时可用下式计算探索点个数n μ n −1 ≤ δ当然也可以不预先计算探索点的数目n 而在计算过程中逐次加以判断看是否已满足了提出的精度要求。0.618 法是一种等速对称进行试探的方法每次的探索点均取在区间长度的 0.618倍和 0.382 倍处。2.2 二次插值法对极小化问题2当 f(t) 在 [a, b] 上连续时可以考虑用多项式插值来进行一维搜索。它的基本思想是在搜索区间中不断用低次通常不超过三次多项式来近似目标函数并逐步用插值多项式的极小点来逼近2的最优解。2.3 无约束极值问题的解法无约束极值问题可表述为min f (x), x ∈ E(n ) 5求解问题5的迭代法大体上分为两点一是用到函数的一阶导数或二阶导数称为解析法。另一是仅用到函数值称为直接法。2.3.1 解析法2.3.1.1 梯度法最速下降法对基本迭代格式x k 1 x k tk pk 6我们总是考虑从点xk 出发沿哪一个方向pk 使目标函数f 下降得最快。微积分的知识告诉我们点xk 的负梯度方向pk −▽f (xk ) 是从点xk 出发使f 下降最快的方向。为此称负梯度方向 − ▽f(xk ) 为f 在点xk 处的最速下降方向。-38-按基本迭代格式6每一轮从点xk 出发沿最速下降方向 − ▽f(xk ) 作一维搜索来建立求解无约束极值问题的方法称之为最速下降法。这个方法的特点是每轮的搜索方向都是目标函数在当前点下降最快的方向。同时用 ▽f (xk ) 0 或作为停止条件。其具体步骤如下1°选取初始数据。选取初始点x0 给定终止误差令 k : 0 。2°求梯度向量。计算 ▽f (xk ) 若 ▽f (xk ) ≤ ε ,停止迭代输出 x k 。否则进行 3°。3° 构造负梯度方向。取pk −▽f (xk ) .4° 进行一维搜索。求tk 使得令xk 1 xk tk pk , k : k 1, 转 2°。例 4 用最速下降法求解无约束非线性规划问题min f (x) x 25x其中x (x1, x2 )T 要求选取初始点 x0 (2,2)T 。解 i ▽f (x) (2x1,50x2 )T编写 M 文件 detaf.m定义函数f (x) 及其梯度列向量如下function [f,df]detaf(x);fx(1)225*x(2)2;df[2x(1)50x(2)];ii编写主程序文件zuisu.m如下clcx[2;2];[f0,g]detaf(x);while norm(g)0.000001p-g/norm(g);t1.0;fdetaf(xtp);while ff0tt/2;fdetaf(xtp);endxxt*p;[f0,g]detaf(x);endx,f02.3.1.2 Newton 法考虑目标函数f 在点xk 处的二次逼近式假定 Hesse 阵-39-正定。由于 ▽2 f (xk ) 正定函数Q 的驻点xk 1 是Q(x) 的极小点。为求此极小点令▽Q(xk 1) ▽f (xk ) ▽2 f (xk )(xk 1 − xk ) 0 即可解得x k 1 xk − [▽2 f (xk )]−1 ▽f (xk ) .对照基本迭代格式1可知从点 xk 出发沿搜索方向。pk −[▽2 f (xk )]−1 ▽f (xk )并取步长 tk 1 即可得 Q(x) 的最小点 xk 1 。通常 把方向 pk 叫做从点 xk 出发的Newton 方向。从一初始点开始每一轮从当前迭代点出发沿 Newton 方向并取步长为 1 的求解方法称之为 Newton 法。其具体步骤如下1°选取初始数据。选取初始点x0 给定终止误差 ε 0 令 k : 0 。2°求梯度向量。计算 ▽f (xk ) 若▽f (xk ) ≤ ε ,停止迭代输出 xk 。否则进行 3°。3°构造 Newton 方向。计算 [▽2 f (xk )]−1 取pk −[▽2 f (xk )]−1 ▽f (xk ) .4° 求下一迭代点。令x k 1 x k pk , k : k 1转 2°。例 5 用 Newton 法求解min f (x) x 25x xx选取x0 (2,2)T 。解i ▽f (x) [4x 2x1x 100x 2xx2 ]Tii编写 M 文件 nwfun.m 如下function [f,df,d2f]nwfun(x);fx(1)425*x(2)4x(1)2*x(2)2;df[4x(1)32*x(1)*x(2)2;100x(2)32*x(1)2x(2)];d2f[2x(1)22*x(2)2,4*x(1)x(2)4x(1)x(2),300x(2)22*x(1)2];III编写主程序文件 example5.m 如下clcx[2;2];[f0,g1,g2]nwfun(x);while norm(g1)0.00001p-inv(g2)*g1;xxp;-40-[f0,g1,g2]nwfun(x);endx, f0如果目标函数是非二次函数一般地说用 Newton 法通过有限轮迭代并不能保证可求得其最优解。为了提高计算精度我们在迭代时可以采用变步长计算上述问题编写主程序文件example5_2 如下clc,clearx[2;2];[f0,g1,g2]nwfun(x);while norm(g1)0.00001p-inv(g2)g1;pp/norm§;t1.0;fnwfun(xtp);while ff0tt/2;fnwfun(xtp);endxxtp;[f0,g1,g2]nwfun(x);endx,f0Newton 法的优点是收敛速度快缺点是有时不好用而需采取改进措施此外当维数较高时计算 − [▽2 f (xk )]−1 的工作量很大。2.3.1.3 变尺度法变尺度法Variable Metric Algorithm是近 20 多年来发展起来的它不仅是求解无约束极值问题非常有效的算法而且也已被推广用来求解约束极值问题。由于它既避免了计算二阶导数矩阵及其求逆过程又比梯度法的收敛速度快特别是对高维问题具有显著的优越性因而使变尺度法获得了很高的声誉。下面我们就来简要地介绍一种变尺度法—DFP 法的基本原理及其计算过程。这一方法首先由Davidon 在 1959 年提出后经 Fletcher 和 Powell 加以改进。我们已经知道牛顿法的搜索方向是 − [▽2 f (xk )]−1▽f (xk ) 为了不计算二阶导数矩阵 [▽2 f (xk )] 及其逆阵我们设法构造另一个矩阵用它来逼近二阶导数矩阵的逆阵 [▽2 f (xk )]−1 这一类方法也称拟牛顿法Quasi-Newton Method。下面研究如何构造这样的近似矩阵并将它记为 。我们要求每一步都能以现有的信息来确定下一个搜索方向每做一次选代目标函数值均有所下降这些近似矩阵最后应收敛于解点处的 Hesse 阵的逆阵。当f (x) 是二次函数时其 Hesse 阵为常数阵 A 任两点xk 和xk 1 处的梯度之差为▽f (xk 1) − ▽f (xk ) A(xk 1 − xk )或对于非二次函数仿照二次函数的情形要求其 Hesse 阵的逆阵的第k 1 次近似矩阵H (k 1) 满足关系式-41-这就是常说的拟 Newton 条件。若令则式7变为现假定H (k ) 已知用下式求H (k 1) 设 H (k ) 和H (k 1) 均为对称正定阵其中 称为第k 次校正矩阵。显然H (k 1) 应满足拟 Newton 条件9即要求或由此可以设想 ΔH (k ) 的一种比较简单的形式是ΔH (k ) Δxk (Q(k ) )T - H (k ) ΔG(k ) (W (k ) )T (12)其中Q(k ) 和W(k ) 为两个待定列向量。将式12中的 ΔH (k ) 代入11得这说明应使(Q(k ) )T ΔG(k ) (W (k ) )T ΔG(k ) 1 13考虑到 ΔH (k ) 应为对称阵最简单的办法就是取由式13得若(Δxk )T ΔG(k ) 和 不等于零则有于是得校正矩阵从而得到上述矩阵称为尺度矩阵。通常我们取第一个尺度矩阵H(0) 为单位阵以后的尺度矩-42-阵按式18逐步形成。可以证明i当 xk 不是极小点且H (k ) 正定时式17右端两项的分母不为零从而可按式18产生下一个尺度矩阵ii若H (k ) 为对称正定阵则由式18产生的 H (k 1) 也是对称正定阵iii由此推出 DFP 法的搜索方向为下降方向。现将 DFP 变尺度法的计算步骤总结如下。1°给定初始点x0 及梯度允许误差 ε 0 。2°若 则 x0 即为近似极小点停止迭代否则转向下一步。 3°令H ( 0) I 单位矩阵p0 −H ( 0) ▽f (x0 )在p0 方向进行一维搜索确定最佳步长 λ0 如此可得下一个近似点x1 x0 λ0 p04 °一般地设已得到近似点xk 算出 ▽f (xk ) 若▽f (x0 ) ≤ ε则x k 即为所求的近似解停止迭代否则计算H (k ) 并令pk −H(k ) ▽f (xk ) 在 pk 方向上进行一维搜索得 λk 从而可得下一个近似点x k 1 xk λk pk5°若xk 1 满足精度要求则xk 1 即为所求的近似解否则转回 4°,直到求出某点满足精度要求为止。2.3.2 直接法在无约束非线性规划方法中遇到问题的目标函数不可导或导函数的解析式难以表示时人们一般需要使用直接搜索方法。同时 由于这些方法一般都比较直观和易于理解因而在实际应用中常为人们所采用。下面我们介绍 Powell 方法。这个方法主要由所谓基本搜索、加速搜索和调整搜索方向三部分组成 具体步骤如下1° 选取初始数据。选取初始点x0 n 个线性无关初始方向组成初搜索方向组{p0 , p1 , … , pn −1 } 。给定终止误差 ε 0 令 k : 0 。2°进行基本搜索。令 y0 : xk 依次沿 {p0 , p1 , … , pn −1 } 中的方向进行一维搜索。对应地得到辅助迭代点y1 , y2 , … , y n 即yj yj −1 tj −1pj−1f j 1, … , n-43-3°构造加速方向。令p n y n − y 0 若 p n ≤ ε ,停止迭代输出 x k 1 y n 。否则进行 4°。4°确定调整方向。按下式f (y m−1) − f (y m ) max{f (yj −1) − f (yj ) | 1 ≤ j ≤ n}找出m 。若f (y 0 ) − 2f (y n) f (2y n − y 0 ) 2[f (y m−1) − f (y m )]成立进行 5°。否则进行 6°。5°调整搜索方向组。令x k 1 y n tn p n : f (y n tn p n) min f (y n tpn ) .t ≥0同时令{p0 , p1 , … , pn 1 }k 1 : {p0 , … , pm−1 , pm1 , … , pn −1 , pn} k : k 1转 2°。6°不调整搜索方向组。令x k 1 : y n , k : k 1转 2°。2.4 Matlab 求无约束极值问题在 Matlab 工具箱中用于求解无约束极值问题的函数有 fminunc 和 fminsearch用法介绍如下。求函数的极小值其中x 可以为标量或向量。Matlab 中 fminunc 的基本命令是[X,FVAL]FMINUNC(FUN,X0,OPTIONS,P1,P2, …)其中的返回值 X 是所求得的极小点FVAL 是函数的极小值其它返回值的含义参见相关的帮助。FUN 是一个 M 文件当 FUN 只有一个返回值时它的返回值是函数f(x)当 FUN 有两个返回值时它的第二个返回值是f(x) 的梯度向量当FUN 有三个返回值时它的第三个返回值是f (x) 的二阶导数阵Hessian 阵。X0 是向量x 的初始值 OPTIONS 是优化参数可以使用缺省参数。P1 P2 是可以传递给 FUN 的一些参数。例 6 求函数f (x) 100(x2 − x12 )2 (1 − x1 )2 的最小值。解编写 M 文件 fun2.m如下function [f,g]fun2(x);f100*(x(2)-x(1)2)2(1-x(1))^2;g[-400x(1)(x(2)-x(1)2)-2*(1-x(1));200*(x(2)-x(1)2)];编写主函数文件example6.m如下options optimset(‘GradObj’,‘on’);[x,y]fminunc(‘fun2’,rand(1,2),options)即可求得函数的极小值。在求极值时也可以利用二阶导数编写 M 文件 fun3.m如下function [f,df,d2f]fun3(x);-44-编写主函数文件example62.m如下options optimset(‘GradObj’,‘on’,‘Hessian’,‘on’);[x,y]fminunc(‘fun3’,rand(1,2),options)即可求得函数的极小值。求多元函数的极值也可以使用 Matlab 的 fminsearch 命令其使用格式为[X,FVAL,EXITFLAG,OUTPUT] FMINSEARCH(FUN,X0,OPTIONS,P1,P2,…)例 7 求函数f (x) sin(x) 3 取最小值时的x 值。解 编写f (x) 的 M 文件 fun4.m 如下function ffun4(x);fsin(x)3;编写主函数文件example7.m如下x02;[x,y]fminsearch(fun4,x0)即求得在初值 2 附近的极小点及极小值。§3 约束极值问题带有约束条件的极值问题称为约束极值问题也叫规划问题。求解约束极值问题要比求解无约束极值问题困难得多。为了简化其优化工作 可采用以下方法将约束问题化为无约束问题将非线性规划问题化为线性规划问题以及能将复杂问题变换为较简单问题的其它方法。库恩—塔克条件是非线性规划领域中最重要的理论成果之一是确定某点为最优点的必要条件但一般说它并不是充分条件对于凸规划它既是最优点存在的必要条件同时也是充分条件。3.1 二次规划若某非线性规划的目标函数为自变量x 的二次函数约束条件又全是线性的就称这种规划为二次规划。Matlab 中二次规划的数学模型可表述如下这里H是实对称矩阵 f , b 是列向量 A 是相应维数的矩阵。Matlab 中求解二次规划的命令是[X,FVAL] QUADPROG(H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS)返回值 X 是决策向量x 的值返回值 FVAL 是目标函数在x 处的值。具体细节可以参看在 Matlab 指令中运行 help quadprog 后的帮助。例 8 求解二次规划解 编写如下程序-45-h[4,-4;-4,8];f[-6;-3];a[1,1;4,1];b[3;9];[x,value]quadprog(h,f,a,b,[],[],zeros(2,1))求得3.2 罚函数法利用罚函数法可将非线性规划问题的求解转化为求解一系列无约束极值问题因而也称这种方法为序列无约束最小化技术简记为 SUMT (Sequential Unconstrained Minization Technique)。罚函数法求解非线性规划问题的思想是利用问题中的约束函数作出适当的罚函数由此构造出带参数的增广目标函数把问题转化为无约束非线性规划问题。主要有两种形式一种叫外罚函数法另一种叫内罚函数法下面介绍外罚函数法。考虑问题min f (x)取一个充分大的数 M 0 构造函数这里 G(x) [g1 (x) … gr (x)] H(x) [h1 (x) … hs (x)]K(x) [k1 (x) … kt (t)]Matlab 中可以直接利用 max 、min 和 sum 函数。则以增广目标函数P(x, M) 为目标函数的无约束极值问题min P(x, M)的最优解x 也是原问题的最优解。例 9 求下列非线性规划解 (i)编写 M 文件 test.mfunction gtest(x);M50000;fx(1)2x(2)28;gf-Mmin(x(1),0)-Mmin(x(2),0)-M*min(x(1)^2-x(2),0)…-46-Mabs(-x(1)-x(2)^22);或者是利用Matlab的求矩阵的极小值和极大值函数编写test.m如下function gtest(x);M50000;fx(1)2x(2)28;gf-Msum(min([x’;zeros(1,2)]))-Mmin(x(1)^2-x(2),0)…Mabs(-x(1)-x(2)^22);我们也可以修改罚函数的定义编写test.m如下 function gtest(x);M50000;fx(1)2x(2)28;gf-Mmin(min(x),0)-Mmin(x(1)2-x(2),0)M*(-x(1)-x(2)22)^2;(ii)在 Matlab 命令窗口输入[x,y]fminunc(‘test’,rand(2,1))即可求得问题的解。3.3 Matlab 求约束极值问题在 Matlab 优化工具箱中用于求解约束最优化问题的函数有fminbnd 、fmincon、 quadprog 、fseminf、fminimax上面我们已经介绍了函数 fmincon 和 quadprog。3.3.1 fminbnd 函数求单变量非线性函数在区间上的极小值Matlab 的命令为[X,FVAL] FMINBND(FUN,x1,x2,OPTIONS)它的返回值是极小点x 和函数的极小值。这里fun 是用 M 文件定义的函数或 Matlab 中的单变量数学函数。例 10 求函数 f (x) (x − 3)2 −1, x ∈ [0,5] 的最小值。解 编写 M 文件 fun5.m function ffun5(x);f(x-3)^2-1;在 Matlab 的命令窗口输入[x,y]fminbnd(‘fun5’,0,5)即可求得极小点和极小值。3.3.2 fseminf 函数求min {F(x) | C(x) ≤ 0, Ceq(x) 0, PHI(x, w) ≤ 0}x其中 C(x),Ceq(x),PHI(x, w) 都是向量函数 w 是附加的向量变量 w 的每个分量都限定在某个区间内。上述问题的 Matlab 命令格式为XFSEMINF(FUN,X0,NTHETA,SEMINFCON,A,B,Aeq,Beq)其中 FUN 用于定义目标函数 F(x) X0 为 x 的初始值 NTHETA 是半无穷约束PHI(x, w) 的个数函数 SEMINFCON 用于定义非线性不等式约束 C(x) 非线性等-47-式约束 Ceq(x) 和半无穷约束PHI(x, w) 的每一个分量函数函数 SEMINFCON 有两个输入参量 X 和 S S 是推荐的取样步长也许不被使用。例 11 求函数 f (x) (x1 - 0.5)2 (x2 - 0.5)2 (x3 - 0.5)2 取最小值时的 x 值,约束为1 ≤ w1 ≤ 100 1 ≤ w2 ≤ 100解 1编写 M 文件 fun6.m 定义目标函数如下function ffun6(x,s);fsum((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%取样值w11:s(1,1):100;w21:s(2,1):100;%半无穷约束k1sin(w1x(1)).cos(w1x(2))-1/1000(w1-50).^2-sin(w1x(3))-x(3)-1;k2sin(w2x(2)).cos(w2x(1))-1/1000*(w2-50).^2-sin(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 函数求解 x F其中F(x) {F1 (x),… , Fm (x)}。上述问题的 Matlab 命令为XFMINIMAX(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON)-48-例 12 求函数族{f1 (x), f2 (x), f3 (x), f4 (x), f5 (x)} 取极大极小值时的x 值。其中解 1编写 M 文件 fun8.m 定义向量函数如下function ffun8(x);2调用函数 fminimax[x,y]fminimax(fun8,rand(2,1))3.3.4 利用梯度求解约束优化问题例 13 已知函数f (x) ex1 (4x 2x 4x1x2 2x2 1) 且满足非线性约束求min f (x)x分析当使用梯度求解上述问题时效率更高并且结果更准确。题目中目标函数的梯度为解 1编写 M 文件 fun9.m 定义目标函数及梯度函数function [f,df]fun9(x);fexp(x(1))(4x(1)22*x(2)24x(1)x(2)2x(2)1);df[exp(x(1))(4x(1)22*x(2)24x(1)x(2)8x(1)6x(2)1);exp(x(1))(4x(2) 4x(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 如下%采用标准算法optionsoptimset(‘largescale’,‘off’);%采用梯度optionsoptimset(options,‘GradObj’,‘on’,‘GradConstr’,‘on’);[x,y]fmincon(fun9,rand(2,1),[],[],[],[],[],[],fun10,options)-49-3.4 Matlab 优化工具箱的用户图形界面解法Matlab 优化工具箱中的 optimtool 命令提供了优化问题的用户图形界面解法。 optimtool 可应用到所有优化问题的求解计算结果可以输出到 Matlab 工作空间中。图 1 优化问题用户图形界面解法示意图例 14 用 optimtool 重新求解例 1。利用例 1 已经定义好的函数 fun1 和 fun2。在 Matlab 命令窗口运行 optimtool就打开图形界面如图 1 所示填入有关的参数未填入的参数取值为空或者为默认值然后用鼠标点一下“start ”按钮就得到求解结果再使用“file”菜单下的“Export to Workspace…”选项把计算结果输出到 Matlab 工作空间中去。§4 飞行管理问题在约 10 000m 高空的某边长 160km 的正方形区域内经常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据以便进行飞行管理。当一架欲进入该区域的飞机到达区域边缘时记录其数据后要立即计算并判断是否会与区域内的飞机发生碰撞。如果会碰撞 则应计算如何调整各架包括新进入的飞机飞行的方向角以避免碰撞。现假定条件如下1不碰撞的标准为任意两架飞机的距离大于 8km;2飞机飞行方向角调整的幅度不应超过 30 度3所有飞机飞行速度均为每小时 800km;4进入该区域的飞机在到达区域边缘时与区域内飞机的距离应在 60km 以上5最多需考虑 6 架飞机6不必考虑飞机离开此区域后的状况。请你对这个避免碰撞的飞行管理问题建立数学模型列出计算步骤对以下数据进行计算方向角误差不超过 0.01 度要求飞机飞行方向角调整的幅度尽量小。设该区域 4 个顶点的座标为(0,0) (160,0) (160,160) (0,160)。记录数据见表 1。-50-表 1 飞行记录数据飞机编号 横座标 x 纵座标 y 方向角度1 150 140 2432 85 85 2363 150 155 220.54 145 50 1595 130 150 230新进入 0 0 52注方向角指飞行方向与x 轴正向的夹角。为方便以后的讨论我们引进如下记号D 为飞行管理区域的边长Ω 为飞行管理区域取直角坐标系使其为[0, D]×[0, D]a 为飞机飞行速度a 800 km/h(x, y) 为第i 架飞机的初始位置(xi (t), yi (t)) 为第i 架飞机在t 时刻的位置θi0 为第i 架飞机的原飞行方向角即飞行方向与x 轴夹角0 ≤ θi0 2π ;Δ θi 为第i 架飞机的方向角调整θi θi0 Δθi 为第i 架飞机调整后的飞行方向角。4.1 模型一根据相对运动的观点在考察两架飞机i 和 j 的飞行时可以将飞机i 视为不动而飞机j 以相对速度vij vj - vi (a cos θj - a cos θi, a sinθj - a sinθi ) 19相对于飞机i 运动对19式进行适当的化约可得不妨设θj ≥ θi 此时相对飞行方向角为 见图 2。图 2 相对飞行方向角由于两架飞机的初始距离为-51-则只要当相对飞行方向角 βij 满足α βij 2π − α 23时两架飞机不可能碰撞见图 2。记 β为调整前第 j 架飞机相对于第i 架飞机的相对速度矢量与这两架飞机连线从j 指向i 的矢量的夹角以连线矢量为基准逆时针方向为正顺时针方向为负。则由式23知两架飞机不碰撞的条件为其中 相对速度vmn 的幅角从n 指向m 的连线矢量的幅角注意 βn 表达式中的i 表示虚数单位这里我们利用复数的幅角可以很方便地计算角度 βn m, n 1,2, … ,6 。本问题中的优化目标函数可以有不同的形式如使所有飞机的最大调整量最小所有飞机的调整量绝对值之和最小等。这里以所有飞机的调整量绝对值之和最小为目标函数可以得到如下的数学规划模型Δ θi ≤ 30o i 1,2, … ,6利用如下的程序clc,clearx0[150 85 150 145 130 0];y0[140 85 155 50 150 0];q[243 236 220.5 159 230 52];xy0[x0; y0];d0dist(xy0); %求矩阵各个列向量之间的距离d0(find(d00))inf;a0asind(8./d0) %以度为单位的反函数xy1x0iy0xy2exp(iq*pi/180)for m1:6for n1:6if n~mb0(m,n)angle((xy2(n)-xy2(m))/(xy1(m)-xy1(n)));end-52-endendb0b0*180/pi;dlmwrite(‘txt1.txt’,a0,‘delimiter’, ‘\t’,‘newline’,‘PC’);fidfopen(‘txt1.txt’,‘a’);fwrite(fid,‘~’,‘char’); %往纯文本文件中写 LINGO 数据的分割符dlmwrite(‘txt1.txt’,b0,‘delimiter’, ‘\t’,‘newline’,‘PC’,‘-append’,‘roffset’, 1)求得 α 的值如表 2 所示。表 2 α 的值1 2 3 4 5 61 0 5.39119 32.23095 5.091816 20.96336 2.2345072 5.39119 0 4.804024 6.61346 5.807866 3.8159253 32.23095 4.804024 0 4.364672 22.83365 2.1255394 5.091816 6.61346 4.364672 0 4.537692 2.9898195 20.96336 5.807866 22.83365 4.537692 0 2.3098416 2.234507 3.815925 2.125539 2.989819 2.309841 0求得 β 的值如表 3 所示。表 3 β 的值1 2 3 4 5 61 0 109.26 -128.25 24.18 173.07 14.4752 109.26 0 -88.871 -42.244 -92.305 93 -128.25 -88.871 0 12.476 -58.786 0.310814 24.18 -42.244 12.476 0 5.9692 -3.52565 173.07 -92.305 -58.786 5.9692 0 1.91446 14.475 9 0.31081 -3.5256 1.9144 0上述飞行管理的数学规划模型可如下输入 LINGO 求解model:sets:plane/1…6/:delta;link(plane,plane):alpha,beta;endsets data:alphafile(‘txt1.txt’); !需要在alpha的数据后面加上分隔符~;betafile(‘txt1.txt’);enddataminsum(plane:abs(delta));for(plane:bnd(-30,delta,30));for(link(i,j) |i#ne#j:abs(beta(i,j)0.5delta(i)0.5delta(j))alpha(i,j));end求得的最优解为 Δθ3 2.83858o Δ θ5 -21.0138o Δ θ6 0.7908o 其它调整角度为 0。4.2 模型二-53-两架飞机i, j 不发生碰撞的条件为(xi (t) − xj (t))2 (yi (t) − yj (t))2 64 251 ≤ i ≤ 5 i 1 ≤ j ≤ 6 0 ≤ t ≤ min{Ti, Tj }其中 Ti, Tj 分别表示第i, j 架飞机飞出正方形区域边界的时刻。这里xi (t) x at cos θi yi (t) y at sinθi i 1,2, … , n θi θi0 Δθi | Δ θi i 1,2, … , n 下面我们把约束条件25加强为对所有的时间t 都成立记li ,j (xi (t) − xj (t))2 (yi (t) − yj (t))2 − 64 (i, j)t 2 (i, j)t (i, j)其中 i, j) 4a2 sin(i, j) 2a[(xi (0) − xj (0))(cosθi − cos θj ) (yi (0) − yj (0))(sinθi − sinθj )](i, j) (xi (0) − xj (0))2 (yi (0) − yj (0))2 − 64则两架i, j 飞机不碰撞的条件是Δ (i, j) (i, j)2 − 4(i, j)(i, j) 0 26这样我们建立如下的非线性规划模型s.t. Δ (i, j) 0 1 ≤ i ≤ 5 i 1 ≤ j ≤ 6i 1,2, … ,6习 题 三用最速下降法梯度法求函数f (x) 4x1 6x2 − 2x − 2x1x2 − 2x的极大点。给定初始点x 0 (1,1)T 。试用牛顿法求解取初始点x (0) (4,0)T 并将采用变步长和采用固定步长 λ 1.0 时的情形做比较。3. 某工厂向用户提供发动机按合同规定其交货数量和日期是第一季度末交40 台第二季末交 60 台第三季末交 80 台。工厂的最大生产能力为每季 100 台每季的生产费用是f (x) 50x 0.2x2 元此处x 为该季生产发动机的台数。若工厂生产的多多余的发动机可移到下季向用户交货这样工厂就需支付存贮费每台发动机每季的存贮费为 4 元。问该厂每季应生产多少台发动机 才能既满足交货合同又使工厂所花费的费用最少假定第一季度开始时发动机无存货。4. 用 Matlab 的非线性规划命令 fmincon 求解飞行管理问题的模型二。5. 用罚函数法求解飞行管理问题的模型二。-54-求下列问题的解max f (x) 2x1 3x 3x2 x x3s.t. x1 2x x2 2x x3 ≤ 10x1 x x2 x − x3 ≤ 50 2x1 x 2x2 x3 ≤ 40 x x3 2x1 2x2 ≥ 1x1 ≥ 0, x2, x3不约束7图 3 影院剖面示意图图 3 为影院的剖面示意图座位的满意程度主要取决于视角 α 和仰角 β 。视角 α 是观众眼睛到屏幕上、下边缘视线的夹角 α 越大越好仰角 β 是观众眼睛到屏幕上边缘视线与水平线的夹角 β 太大使人的头部过分上仰引起不舒服感一般要求 β 不超过30o 。记影院屏幕高h 上边缘距地面高H 地板线倾角θ , 第一排和最后一排座位与屏幕水平距离分别为d 和D 观众平均座高为c指眼睛到地面的距离。已知参数h 1.8 H 5 d 4.5 D 19 c 1.1 单位m。1地板线倾角 θ 10o 问最佳座位在什么地方2求地板线倾角一般不超过20o 使所有观众的平均满意程度最大。3地板线设计成什么形状可以进一步提高观众的满意程度。