《数学建模》算法全收录（算法大全）：第03章非线性规划.pdf_大学文库

其中x=[x1…xn称为模型(NP)的决策变量,∫称为目标函数,g;(i=1…,pP) 和h(=1,…,q)称为约束函数。另外,g;(x)=0(=l…,P)称为等式约束, h(x)≤0G=1,…q)称为不等式的约束。对于一个实际问题,在把它归结成非线性规划问题时,一般要注意如下几点 (i)确定供选方案:首先要收集同问题有关的资料和数据,在全面熟悉问题的基础上,确认什么是问题的可供选择的方案,并用一组变量来表示它们。 ⅱ)提出追求目标:经过资料分析,根据实际需要和可能,提出要追求极小化或极大化的目标。并且,运用各种科学和技术原理,把它表示成数学关系式。你>(ii)给出价值标准:在提出要追求的目标之后,要确立所考虑目标的“好”或的价值标准,并用某种数量形式来描述它。 (iv)寻求限制条件:由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果,因此还需要寻找出问题的所有限制条件,这些条件通常用变量之间的一些不等式或等式来表示 12线性规划与非线性规划的区别如果线性规划的最优解存在,其最优解只能在其可行域的边界上达到(特别是可行域的顶点上达到);而非线性规划的最优解(如果最优解存在)则可能在其可行域的任意一点达到。 1.3非线性规划的 Matlab解法 Matlab中非线性规划的数学模型写成以下形式 min f(x) Ax< B Aeg·x=Beq (x)≤0 ICeq(x)=0 其中f(x)是标量函数,A,B,Aeqg,Beq是相应维数的矩阵和向量,C(x),Ceq(x)是非线性向量函数。 Matlab中的命令是 X=FMINCON (FUN, XO, A, B, Aeg, Beg, 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 的各分量都为inf:; NONLCON是用M文件定义的非线性向量函数C(x),Ceq(x): OPTIONS 定义了优化参数,可以使用 Matlab缺省的参数设置例2求下列非线性规划 min f(x)=x+x2+x3+8 x2-x2+x2≥0 x1+x2+x3≤20 -x,-x,+2=0

-33- 其中 T n x [x x ] = 1 L 称为模型（NP）的决策变量，f 称为目标函数，gi (i = 1,L, p) 和 h ( j 1, ,q) j = L 称为约束函数。另外， gi(x) = 0 (i = 1,L, p) 称为等式约束， hj(x) ≤ 0 ( j = 1,L,q) 称为不等式的约束。对于一个实际问题，在把它归结成非线性规划问题时，一般要注意如下几点：（i）确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基础上，确认什么是问题的可供选择的方案，并用一组变量来表示它们。（ii）提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。（iii）给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或 “坏”的价值标准，并用某种数量形式来描述它。（iv）寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果，因此还需要寻找出问题的所有限制条件，这些条件通常用变量之间的一些不等式或等式来表示。 1.2 线性规划与非线性规划的区别如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到）；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。 1.3 非线性规划的 Matlab 解法 Matlab 中非线性规划的数学模型写成以下形式 min f (x) ⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ = ≤ ⋅ = ≤ ( ) 0 ( ) 0 Ceq x C x Aeq x Beq Ax B , 其中 f (x)是标量函数， A, B, Aeq, Beq是相应维数的矩阵和向量，C(x),Ceq(x) 是非线性向量函数。 Matlab 中的命令是 X=FMINCON(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 的各分量都为 inf；NONLCON 是用 M 文件定义的非线性向量函数C(x),Ceq(x) ；OPTIONS 定义了优化参数，可以使用 Matlab 缺省的参数设置。例2 求下列非线性规划 min ( ) 8 2 3 2 2 2 f x = x1 + x + x + s.t. 0 2 2 3 2 x1 − x + x ≥ 20 3 3 2 x1 + x2 + x ≤ 2 0 2 − x1 − x2 + =

-34- 2 3 2 x2 + x3 = x1, x2 , x3 ≥ 0 解（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(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) 就可以求得当 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 ) I K * * ≤ ∀ ∈ δ ，则称 * x 是（NP）的局部最优解， ( ) * f x 是(NP)的局部最优值。如果有 f (x ) f (x), x N (x ) I K * * < ∀ ∈ δ 则称 * x 是（NP）的严格局部最优解， ( ) * f x 是(NP)的严格局部最优值。由于线性规划的目标函数为线性函数，可行域为凸集，因而求出的最优解就是整个可行域上的全局最优解。非线性规划却不然，有时求出的某个解虽是一部分可行域上的极值点，但却并不一定是整个可行域上的全局最优解。对于非线性规划模型(NP)，可以采用迭代方法求它的最优解。迭代方法的基本思想是：从一个选定的初始点 n x ∈ R 0 出发，按照某一特定的迭代规则产生一个点列 { } k x ，使得当{ } k x 是有穷点列时，其最后一个点是(NP)的最优解；当{ } k x 是无穷点列时，它有极限点，并且其极限点是(NP)的最优解。设 k n x ∈ R 是某迭代方法的第k 轮迭代点， k n x ∈ R +1 是第k +1轮迭代点，记 k k k k x = x + t p +1 （1）这里 , , 1 1 ∈ ∈ = k n k tk R p R p ，并且 k p 的方向是从点 k x 向着点 k +1 x 的方向。式（1）就是求解非线性规划模型(NP)的基本迭代格式

-35- 通常，我们把基本迭代格式（1）中的 k p 称为第k 轮搜索方向， kt 为沿 k p 方向的步长，使用迭代方法求解(NP)的关键在于，如何构造每一轮的搜索方向和确定适当的步长。设 x ∈ R , p ≠ 0 n ，若存在δ > 0 ，使 f (x + tp) 0，使 x + tp ∈ K ，称向量 p 是点 x 处关于 K 的可行方向。一个向量 p ，若既是函数 f 在点 x 处的下降方向，又是该点关于区域 K 的可行方向，称之为函数 f 在点 x 处关于 K 的可行下降方向。现在，我们给出用基本迭代格式（1）求解(NP)的一般步骤如下： 0° 选取初始点 0 x ，令k := 0 。 1° 构造搜索方向，依照一定规划，构造 f 在点 k x 处关于 K 的可行下降方向作为搜索方向 k p 。 2° 寻求搜索步长。以 k x 为起点沿搜索方向 k p 寻求适当的步长 kt ，使目标函数值有某种意义的下降。 3° 求出下一个迭代点。按迭代格式（1）求出 k k k k x = x + t p +1 。若 k +1 x 已满足某种终止条件，停止迭代。 4° 以 k +1 x 代替 k x ，回到 1°步。 1.5 凸函数、凸规划设 f (x) 为定义在 n 维欧氏空间 (n) E 中某个凸集 R 上的函数，若对任何实数 α(0 <α <1) 以及 R 中的任意两点 (1) x 和 (2) x ，恒有 ( (1 ) ) ( ) (1 ) ( ) (1) (2) (1) (2) f αx + −α x ≤αf x + −α f x 则称 f (x)为定义在 R 上的凸函数。若对每一个α(0 <α <1) 和 x ≠ x ∈R (1) (2) 恒有 ( (1 ) ) ( ) (1 ) ( ) (1) (2) (1) (2) f αx + −α x <αf x + −α f x 则称 f (x)为定义在 R 上的严格凸函数。考虑非线性规划 ⎪⎩ ⎪ ⎨ ⎧ = ≤ = ∈ { | ( ) 0, 1,2, , } min ( ) R x g x j l f x j x R L 假定其中 f (x)为凸函数，g (x)( j 1,2, ,l) j = L 为凸函数，这样的非线性规划称为凸规划。可以证明，凸规划的可行域为凸集，其局部最优解即为全局最优解，而且其最优解的集合形成一个凸集。当凸规划的目标函数 f (x)为严格凸函数时，其最优解必定唯一（假定最优解存在）。由此可见，凸规划是一类比较简单而又具有重要理论意义的非

线性规划 §2无约束问题 2.1一维搜索方法当用迭代法求函数的极小点时,常常用到一维搜索,即沿某一已知方向求目标函数的极小点。一维搜索的方法很多,常用的有:(1)试探法(“成功一失败”,斐波那契法 0618法等);(2)插值法(抛物线插值法,三次插值法等);(3)微积分中的求根法(切线法,二分法等)。考虑一维极小化问题 min f(o) (2) ast≤b 若∫()是[an,b]区间上的下单峰函数,我们介绍通过不断地缩短[a,b]的长度,来搜索得(2)的近似最优解的两个方法为了缩短区间[a,b],逐步搜索得(2)的最优解r的近似值,我们可以采用以下途径:在[a,b]中任取两个关于[a,b]是对称的点1和l2(不妨设l2<l1,并把它们叫做搜索点),计算∫(1)和∫(l2)并比较它们的大小。对于单峰函数,若f(2)<f(1), 则必有r∈[a,l],因而[a,1]是缩短了的单峰区间;若f(t1)<f(l2),则有 t∈[l2,b,故[t2,b是缩短了的单峰区间;若∫(t2)=f(1),则[a,41]和[t2,b都是缩短了的单峰。因此通过两个搜索点处目标函数值大小的比较,总可以获得缩短了的单峰区间。对于新的单峰区间重复上述做法,显然又可获得更短的单峰区间。如此进行, 在单峰区间缩短到充分小时,我们可以取最后的搜索点作为(2)最优解的近似值应该按照怎样的规则来选取探索点,使给定的单峰区间的长度能尽快地缩短? 2.1.1 Fibonacci法如用F表示计算n个函数值能缩短为单位长区间的最大原区间长度,可推出F满足关系 Fo=F Fn=Fn-2+F1,n=2,3 数列{F}称为 Fibonacci数列,F称为第n个 Fibonacci数,相邻两个 Fibonacci数之比n称为 Fibonacci分数 F 当用斐波那契法以n个探索点来缩短某一区间时,区间长度的第一次缩短率为其后各次分别为 F 由此,若t1和t2(12<1)是单峰区间[a,b] 中第1个和第2个探索点的话,那么应有比例关系 F 从而 t1=a+2(b-a),l2=a+-n2(b-a) (3) 它们关于[a,b]确是对称的点 -36

-36- 线性规划。 §2 无约束问题 2.1 一维搜索方法当用迭代法求函数的极小点时，常常用到一维搜索，即沿某一已知方向求目标函数的极小点。一维搜索的方法很多，常用的有：（1）试探法（“成功—失败”，斐波那契法， 0.618 法等）；（2）插值法（抛物线插值法，三次插值法等）；（3）微积分中的求根法（切线法，二分法等）。考虑一维极小化问题 min f (t) a≤t≤b （2）若 f (t) 是[a,b] 区间上的下单峰函数，我们介绍通过不断地缩短[a,b] 的长度，来搜索得（2）的近似最优解的两个方法。为了缩短区间[a,b] ，逐步搜索得（2）的最优解 * t 的近似值，我们可以采用以下途径：在[a,b] 中任取两个关于[a,b] 是对称的点 1 t 和 2t （不妨设 2 1 t < t ，并把它们叫做搜索点），计算 ( )1 f t 和 ( ) 2 f t 并比较它们的大小。对于单峰函数，若 ( ) ( ) 2 1 f t < f t ，则必有 [ , ]1 * t ∈ a t ，因而 [ , ]1 a t 是缩短了的单峰区间；若 ( ) ( ) 1 2 f t < f t ，则有 [ , ] 2 * t ∈ t b ，故[ , ] 2t b 是缩短了的单峰区间；若 ( ) ( ) 2 1 f t = f t ，则[ , ]1 a t 和[ , ] 2t b 都是缩短了的单峰。因此通过两个搜索点处目标函数值大小的比较，总可以获得缩短了的单峰区间。对于新的单峰区间重复上述做法，显然又可获得更短的单峰区间。如此进行，在单峰区间缩短到充分小时，我们可以取最后的搜索点作为（2）最优解的近似值。应该按照怎样的规则来选取探索点，使给定的单峰区间的长度能尽快地缩短？ 2.1.1 Fibonacci 法如用 Fn 表示计算n 个函数值能缩短为单位长区间的最大原区间长度，可推出 Fn 满足关系 F0 = F1 = 1 , 2,3, , Fn = Fn−2 + Fn−1 n = L 数列{ } Fn 称为 Fibonacci 数列， Fn 称为第n 个 Fibonacci 数，相邻两个 Fibonacci 数之比 n n F F −1 称为 Fibonacci 分数。当用斐波那契法以 n 个探索点来缩短某一区间时，区间长度的第一次缩短率为 n n F F −1 ，其后各次分别为 2 1 2 3 1 2 , , , F F F F F F n n n n L − − − − 。由此，若 1 t 和 ( ) 2 2 1 t t < t 是单峰区间[a,b] 中第 1 个和第 2 个探索点的话，那么应有比例关系 n n F F b a t1 a −1 = − − ， n n F F b a t2 a −2 = − − 从而 ( ) 1 1 b a F F t a n n = + − − ， ( ) 2 2 b a F F t a n n = + − − （3）它们关于[a,b] 确是对称的点

-38- ω ω −ω = 1 1 称之为黄金分割数，其值为 0.6180339887L 2 5 1 = − ω = 。黄金分割数ω 和 Fibonacci 分数之间有着重要的关系 n n n F F 1 lim − →∞ ω = 。现用不变的区间缩短率 0.618，代替斐波那契法每次不同的缩短率，就得到了黄金分割法（0.618 法）。这个方法可以看成是斐波那契法的近似，实现起来比较容易，效果也相当好，因而易于为人们所接受。用 0.618 法求解，从第 2 个探索点开始每增加一个探索点作一轮迭代以后，原单峰区间要缩短 0.618 倍。计算n 个探索点的函数值可以把原区间[ , ] a0 b0 连续缩短n −1 次，因为每次的缩短率均为 μ ，故最后的区间长度为 1 0 0 ( ) − − n b a μ 这就是说，当已知缩短的相对精度为δ 时，可用下式计算探索点个数n ： μ ≤ δ n−1 当然，也可以不预先计算探索点的数目n ，而在计算过程中逐次加以判断，看是否已满足了提出的精度要求。 0.618 法是一种等速对称进行试探的方法，每次的探索点均取在区间长度的 0.618 倍和 0.382 倍处。 2.2 二次插值法对极小化问题（2），当 f (t) 在[a,b] 上连续时，可以考虑用多项式插值来进行一维搜索。它的基本思想是：在搜索区间中，不断用低次（通常不超过三次）多项式来近似目标函数，并逐步用插值多项式的极小点来逼近（2）的最优解。 2.3 无约束极值问题的解法无约束极值问题可表述为 ( ) min ( ), n f x x ∈ E （5）求解问题（5）的迭代法大体上分为两点：一是用到函数的一阶导数或二阶导数，称为解析法。另一是仅用到函数值，称为直接法。 2.3.1 解析法 2.3.1.1 梯度法（最速下降法）对基本迭代格式 k k k k x = x + t p +1 （6）我们总是考虑从点 k x 出发沿哪一个方向 k p ，使目标函数 f 下降得最快。微积分的知识告诉我们，点 k x 的负梯度方向 ( ) k k p = −∇f x ，是从点 k x 出发使 f 下降最快的方向。为此，称负梯度方向 ( ) k − ∇f x 为 f 在点 k x 处的最速下降方向

-41- [f0,g1,g2]=nwfun(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 f>f0 t=t/2;f=nwfun(x+t*p); end x=x+t*p; [f0,g1,g2]=nwfun(x); end x,f0 Newton 法的优点是收敛速度快；缺点是有时不好用而需采取改进措施，此外，当维数较高时，计算 2 1 [ ( )]− − ∇ k f x 的工作量很大。 2.3.1.3 变尺度法变尺度法（Variable Metric Algorithm）是近 20 多年来发展起来的，它不仅是求解无约束极值问题非常有效的算法，而且也已被推广用来求解约束极值问题。由于它既避免了计算二阶导数矩阵及其求逆过程，又比梯度法的收敛速度快，特别是对高维问题具有显著的优越性，因而使变尺度法获得了很高的声誉。下面我们就来简要地介绍一种变尺度法—DFP 法的基本原理及其计算过程。这一方法首先由 Davidon 在 1959 年提出，后经 Fletcher 和 Powell 加以改进。我们已经知道，牛顿法的搜索方向是 [ ( )] ( ) 2 k 1 k − ∇ f x ∇f x − ，为了不计算二阶导数矩阵[ ( )] 2 k ∇ f x 及其逆阵，我们设法构造另一个矩阵，用它来逼近二阶导数矩阵的逆阵 2 1 [ ( )]− ∇ k f x ，这一类方法也称拟牛顿法（Quasi-Newton Method）。下面研究如何构造这样的近似矩阵，并将它记为 (k ) H 。我们要求：每一步都能以现有的信息来确定下一个搜索方向；每做一次选代，目标函数值均有所下降；这些近似矩阵最后应收敛于解点处的 Hesse 阵的逆阵。当 f (x)是二次函数时，其 Hesse 阵为常数阵 A ，任两点 k x 和 k +1 x 处的梯度之差为 ( ) ( ) ( ) k 1 k k 1 k ∇f x − ∇f x = A x − x + + 或 [ ( ) ( )] k 1 k 1 k 1 k x − x = A ∇f x − ∇f x + − + 对于非二次函数，仿照二次函数的情形，要求其 Hesse 阵的逆阵的第k +1次近似矩阵 (k +1) H 满足关系式 [ ( ) ( )] k 1 k (k 1) k 1 k x − x = H ∇f x − ∇f x + + + （7）

《数学建模》算法全收录（算法大全）：第03章 非线性规划

《数学建模》算法全收录（算法大全）：第03章非线性规划