Lecture5-6非线性规划及Matlab实现_修改.ppt

华南理工大学 模具研究室 非线性规划基本概念及分类 当目标函数或约束条件中有一个或多个为非线性函数，则称 这样的规划问题为非线性规划（Nonlinear Programming） 。其数学模型为： 非线性规划 无约束非线性规划约束非线性规划 梯度(gradient) 设n元函数f(x)在点x处可微，则称以下向量为f(x)在点x处的梯度 梯度的几何意义：如果函数f(x)在点x处的梯度 是非零向量 ，那么 就是f(x)的等值面在点x处的法向量，垂直于等值面 在点x处的切平面，且指向f(x)的函数值增大的方向。 Hessian矩阵(Hessian Matrix) 设n元函数f(x)在点x处二次可微，将f(x)在x点的二阶偏导数按如 下组成 ，则称 是函数f(x)在x点的Hessian矩阵。 Hessian矩阵和梯度的内在关系 方向导数(direction derivative) 设f(x)在点x处可微，d是给定的非零向量，如果极限： 存在，则称此极限为函数f(x)在点x处沿着方向d的方向导数，记 ： 下降方向(descent direction) 设d是非零向量，若存在正数 0，当t(0, )时，必有： 则称方向d是函数f(x)在点x处的一个下降方向。 如果函数f(x)在点x沿着方向d的方向导数满足条件： 那么方向d必是函数f(x)在点x处的一个下降方向。 负梯度方向称为最速下降方向。 参看： Lecture 5 Non-linear Programming.pdf 3.2.1节 判断是否下降方向的一个简单方法就是看该方向与负梯度 方向之间的夹角。 正定矩阵 考虑二次型 ，z 为n维向量 正定的二次型：若对于任意 ，有 ； 半正定的二次型：若对于任意 ，有 ； 负定的二次型：若对于任意 ，有 ； 半负定的二次型：若对于任意 ，有 ； 不定二次型： ，有 ，又 ，有 . 二次型 为正定的充要条件是它的矩阵 的左上角各阶主 子式都大于零. 矩阵M的所有的特征值i都是正的。 函数Taylor展开 相关数学基础知识参看： Lecture 4 Mathematical Preliminaries.pdf 无约束非线性规划最优性条件 局部极小点、全局极小点、非光滑的极小点 局部极小的条件极小点的类型 无约束问题的最优性条件 无约束极小点的一阶必要条件： 设f(x)在点x处可微，若x是f(x)的无约束局部极小点，则 必有梯度 无约束极小点的二阶充分条件： 设f(x)在点x处二次可微，若x是f(x)的无约束局部极小点，则必有 梯度 ，且f(x)在x处的Hessian矩阵 半正定。 严格无约束局部极小点充分条件： 设f(x)在x处二次可微，若梯度 ，且f(x)在x处的 Hessian矩阵 正定，则可断言x是严格局部极小点 例：利用极值条件解优化问题 算法及相关概念 1、迭代算法 集合 D上的迭代算法A： （1）初始点； （2）按照某种规则A产生下一个迭代点。 （i）如果点列收敛于最优解，则称算法A收敛。 （ii）如果 ，则称算法A为 下降迭代算法。 . . 2.下降迭代算法步骤 （1）给出初始点，令； （2）按照某种规则确定下降搜索方向； （3）按照某种规则确定搜索步长，使得 ； （4）令， ； （5）判断是否满足停止条件。是则停止，否则转第2步。 搜索步长确定方法： 称。为最优步长，且有 模型算法 线性搜索求 ， 新点 使x(k+1)S 初始x(1) S, k =1 对x(k)点选择下降 可行方向d(k) 是否满足停机条件？ 停 k=k+1 yes no 3.终止条件 b. d. a. c. xk xk+1 xk xk+1 x* 一般而言，线性收敛速度是比较慢的，超线性收敛 速度相对较快，而二阶收敛速度则相当快。如果一 个算法具有超线性收敛速度，那么从计算的角度就 可以认为是一个比较好的算法 则称 的收敛阶为 。 a. 设算法A所得的点列为 ，如果 b. 4.收敛速度 单谷函数（单峰函数） ： 设f(t)是定义在区间a,b上的一元函数，t*是f(t)在a,b 上的全局极小点。如果在a,b上任取两个点t1f(t2)；而当t1t*时，必有f(t1)f(t2)，则搜索区间可缩短为t1,b； 如果f(t1)f(t2)，则搜索区间可缩短为a,t2； 无约约束优优化：线线搜索法 给定初始估计x(0)，设x(k)处有g(k) 0，则第 k 次迭代： (a)中的搜索方向有许多种不同的确定方法！ (a)根据某种模型函数确定x(k)处的搜索方向d(k) (b)线搜索：确定 的极小点 (c) 置 (b)中的确定步长：精确线搜索、非精确线搜索 一维收索的方法很多，大致可以分为两类：一为不使 用导数的方法，如进退法、Fibonacci法、0.618法；二 为使用导数的方法，如对分法、牛顿法、抛物线法、 三次插值法等。 确定 的极小点 确定步长方法：分析法（最优步长） 最优步长参看： Lecture 5 Non-linear Programming.pdf 节 进退法求初始不确定区间 找三点使两端点的函数值大于中间点的函数值。 思路：任取0，步长 0，取1=0 + ， 1若(0 )(0),则停，a= 2，b= 1 (图1) 2若(0 )(1), 令=2 , 2=1 + ， 若(2 )(1),则停，a= 0 ，b= 2 (图2) 2 0 1 向左找 2 图1 0 1 2 向右找2 图2 区间法参看： Lecture 5 Non-linear Programming.pdf 节 0.618法和Fibonacci法 0.618法和Fibonacci法都是分割方法，其基本思想是通 过取试探点和函数值进行比较，使包含极小点的搜索区 间不断缩短，当区间长度缩短到一定程度时，区间上各 点的函数值均接近极小值，从而各点可以看作是极小点 的近似。 这类方法仅需计算函数值，不涉及导数，又称直接法， 用途很广，尤其适用于非光滑及导数表达式复杂或写不 出的情况。 这类方法要求所考虑区间上的目标函数是单峰函数。如 果不满足这个条件，可以考虑把区间分成若干个小区间 ，使得每个小区间上函数都是单峰，分别求极小值 0.618法 0.618法的基本思想是在搜索区间a,b上选取两个对称 点1和2 ，且10 = + （1-t）( - ) = +t( - ) - 0? No = , = = +t( - ) yes = , = = + （1-t）( - ) No 黄金分割法Matlab实现程序，参看golden_sect.pdf Fibonacci法 0.618法和Fibonnacci法都用于单谷函数，通过不断缩小搜索 区间来获得极小点的近似值。0.618法的区间长度缩短比率是 常数，而Fibonacci法的区间长度缩短比率有Fibonnacci数确 定。当所需插入点的个数n充分大时，Fibonacci法的搜索区 间缩短比率趋近于0.618法的缩短比率。 最速下降法 令 最速下降法参看： Lecture 5 Non-linear Programming.pdf 节 最优解 最优方向 克服负梯度方向缺陷的途径 用适当的正定矩阵（尺度矩阵）乘负梯度方向，其 作用是对后者进行适当的旋转，以获得更好的方向 牛顿法 在局部用一个二次函数 g(x)近似地代替目标函数 ，然后用g(x)的极小值作 为f(x)的近似极小点 牛顿方向 牛顿法特点 牛顿法的改进: 修正牛顿法/ 阻尼牛顿法 牛顿法参看： Lecture 5 Non-linear Programming.pdf 节 变尺度法/拟牛顿法 前面介绍了牛顿法,它的突出优点是收敛很快.但是,运用 牛顿法需要计算二阶导数,而且目标函数的Hessian矩阵可 能非正定.为了克服牛顿法的缺点,人们提出了拟牛顿法. 它的基本思想是用不包含二阶导数的矩阵近似牛顿法中 的Hessian矩阵的逆矩阵.由于构造近似矩阵的方法不同,因 而出现不同的拟牛顿法.经理论证明和实践检验,拟牛顿法 已经成为一类公认的比较有效的算法. 下面分析怎样构造近似矩阵并用它取代牛顿法中的Hessian 矩阵的逆. 前面已经给出牛顿法的迭代公式,即 其中 是在点 处的牛顿方向: 是从 出发沿牛顿方向搜索的最优步长. 为构造 的近似矩阵 ,先分析 与一 阶导数的关系. 设在第 次迭代后,得到点 ,我们将目标函数 在点 展成Taylor级数,并取二阶近似,得到 由此可知,在 附近有 令 ,则 记作 则有 又设Hessian矩阵 可逆,则 这样,计算出 后,可以根据上式估计在 处的Hessian矩 阵的逆.因此,为了用不包含二阶导数的矩阵 取代牛顿法中的 Hessian矩阵 的逆矩阵,有理由令 满足 该式有时称为拟牛顿条件. 下面就来研究怎样确定满足这个条件的矩阵 . DFP算法 著名的 DFP 方法是 Davidon 首先提出, 后来又被 Fletcher 和 Powell 改进的算法,又称为变尺度法.在这种方法中,定义校正矩阵 为 容易验证,这样定义校正矩阵 ,得到的矩阵 DFP方法计算步骤如下: 1.给定初始点 ,允许误差 . 2.置 (单位矩阵),计算出在 处的梯度 3.令 4.从 出发,沿方向 搜索,求步长 ,使它满足 令 5.检验是否满足收敛准则,若 则停止迭代,得到点 ;否则,进行步6 6.若 ,则令 ,返回步2.;否则,进行步7. （即迭代一轮n次仍没求得最优解，以新的 为起点重新开始一轮新的迭代） 7.令 利用公式计算 置 ,返回步3 拟牛顿法参看： “An Introduction to Optimization” 第11章，以及Lecture 5 Non-linear Programming.pdf 节 共轭方向法 最速下降法计算步骤简单，但收敛速度太慢； 牛顿法和阻尼牛顿法收敛速度快，但要计算二阶偏导数 矩阵及其逆矩阵，计算量太大。 共轭方向法是兼有这两种方法优点的一类方法， 比最速下降法的收敛速度快得多，但同时又避免了 牛顿法所要求的Hessian矩阵的计算和求逆。 共轭方向法，主要是其中的共轭梯度法（conjugate gradient method）。1952年Hesteness和Stiefel为求解线性方程组提出 了共轭梯度法，1964年Fletcher和Reeves提出了F-R共轭梯度 法（简称FR法），用来求解无约束优化问题。 共轭梯度法 FR法的基本思想是把共轭性与最速下降法相结合，利用已 知点处的梯度构造一组共轭方向，并沿着这组方向进行搜 索，求出目标函数的极小值。 其中 , 是对称正定矩阵, 是常数. 我们先讨论对于二次凸函数的共轭梯度法,然后再把这种方 法推广到极小化一般函数的情形. 考虑问题 求解方法 ： Fletcher-Reeves 公式 拟牛顿法参看： “An Introduction to Optimization” 第10章，以及Lecture 5 Non-linear Programming.pdf 节 约束极值的最优化方法 约束问题的最优性条件 研究一点处的可行方向时，只需考虑在这一点的积极约束，可以 暂时不管那些不积极约束 可行方向 Lagrange必要条件 Lagrange定理实 质是将等式约束 问题转化为无约 束问题来处理 等式约束问题的最优性条件 Lagrange充分条件 不等式约束问题的最优性条件 几何最优性条件 几何最优性 条件是不等 式约束条件 下极小点的 必要条件 可行方向代数条件 下降方向代数条件 可行下降方向 KuhnTucker条件 Kuhn-Tucker定理简称为K- T定理，称所导出的有关极 小点的必要条件为K-T条件 ，满足K-T条件的点称为K-T 点 KT条件几何解释参看： Lecture 6 Mathematical Preliminaries2.pdf Zoutendijk可行方向法 可行方向法基本原理 仅带线性约束的Zoutendijk可行方向法 带非线性约束的Zoutendijk可行方向法 求出线性规划问题的最优解 ，若 z=0，说明在 点处不存 在可行下降方向.在 为正则点 时， x为KT点.若z0 ，则 为可行下降方向. 制约函数法 这类方法是将约束条件随同一个参数并入到目标函 数而形成一个新的无约束问题.当参数变化时就得到 了一系列无约束问题，用求解这一系列无约束问题 而得到的解来逼近原有约束问题的解，也称这类方 法为序列无约束极小化方法. 罚函数法(外点法 ) 障碍函数法（内点法） 罚函数法(外点法 ) 障碍函数法（内点法） 用MATLAB软件求解,其输入格式如下: 1x=quadprog(H,C,A,b); 2x=quadprog(H,C,A,b,Aeq,beq); 3x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB); 4x=quadprog(H,C,A,b, Aeq,beq ,VLB,VUB,X0); 5x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB,X0,options); 6x,fval=quaprog(); 7x,fval,exitflag=quaprog(); 8x,fval,exitflag,output=quaprog(); 1二次规划 例1 min f(x1,x2)=-2x1-6x2+x12-2x1x2+2x22 s.t. x1+x22 -x1+2x22 x10, x20 MATLAB（youh1） 1写成标准形式： 2输入命令： H=2 -2; -2 4; c=-2 ;-6;A=1 1; -1 2;b=2;2; Aeq=;beq=; VLB=0;0;VUB=; x,z=quadprog(H,c,A,b,Aeq,beq,VLB,VUB) 3运算结果为： x =4/5 6/5 z = -36/5 s.t. 1 首先建立M文件fun.m,用来定义目标函数F（X）: function f=fun(X); f=F(X); 2一般非线性规划 其中X为n维变元向量，G(X)与Ceq(X)均为非线性函数组成的向量，其 他变量的含义与线性规划、二次规划中相同用MATLAB求解上述问题， 基本步骤分三步： 3 建立主程序.求解非线性规划的函数是fmincon,命令的基本格式 如下： (1) x=fmincon(fun,X0,A,b) (2) x=fmincon(fun,X0,A,b,Aeq,beq) (3) x=fmincon(fun,X0,A,b, Aeq,beq,VLB,VUB) (4) x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon) (5)x=fmincon(fun,X0,A,b,Aeq,beq,VLB,VUB,nonlcon,options) (6) x,fval= fmincon() (7) x,fval,exitflag= fmincon() (8)x,fval,exitflag,output= fmincon() 输出极值点 M文件迭代的初值参数说明变量上下限 注意： 1 fmincon函数提供了大型优化算法和中型优化算法默认 时: 若在fun函数中提供了梯度（options参数的GradObj设置 为on），并且只有上下界存在或只有等式约束，fmincon函 数将选择大型算法当既有等式约束又有梯度约束时，使用中型 算法 2 fmincon函数的中型算法使用的是序列二次规划法在每 一步迭代中求解二次规划子问题，并用BFGS法更新拉格朗日 Hesse矩阵 3 fmincon函数可能会给出局部最优解，这与初值X0的选取 有关 1写成标准形式： s.t. 2x1+3x2 6 s.t. x1+4x2 5 x1,x2 0 例2 2先建立M-文件 fun3m: function f=fun3(x); f=-x(1)-2*x(2)+(1/2)*x(1)2+(1/2)*x(2)2 MATLAB(youh2) 3再建立主程序youh2m： x0=1;1; A=2 3 ;1 4; b=6;5; Aeq=;beq=; VLB=0;0; VUB=; x,fval=fmincon(fun3,x0,A,b,Aeq,beq,VLB,VUB ) 4运算结果为： x = 07647 10588 fval = -20294 1先建立M文件fun4m定义目标函数: function f=f