Lecture5-6非线性规划及Matlab实现课件

第三章非线性规划及Matlab实现华南理工大学模具研究室mesjzhang@非线性规划基本概念及分类当目标函数或约束条件中有一个或多个为非线性函数，则称这样的规划问题为非线性规划（NonlinearProgramming）。其数学模型为：非线性规划无约束非线性规划约束非线性规划梯度(gradient)设n元函数f(x)在点x处可微，则称以下向量为f(x)在点x处的梯度梯度的几何意义：如果函数f(x)在点x处的梯度是非零向量，那么就是f(x)的等值面在点x处的法向量，垂直于等值面在点x处的切平面，且指向f(x)的函数值增大的方向。Hessian矩阵(HessianMatrix)设n元函数f(x)在点x处二次可微，将f(x)在x点的二阶偏导数按如下组成，则称是函数f(x)在x点的Hessian矩阵。Hessian矩阵和梯度的内在关系方向导数(directionderivative)设f(x)在点x处可微，d是给定的非零向量，如果极限：存在，则称此极限为函数f(x)在点x处沿着方向d的方向导数，记：下降方向(descentdirection)设d是非零向量，若存在正数ε>0，当tЄ(0,ε)时，必有：则称方向d是函数f(x)在点x处的一个下降方向。如果函数f(x)在点x沿着方向d的方向导数满足条件：那么方向d必是函数f(x)在点x处的一个下降方向。负梯度方向称为最速下降方向。参看：Lecture5Non-linearProgramming.pdf3.2.1节判断是否下降方向的一个简单方法就是看该方向与负梯度方向之间的夹角。正定矩阵考虑二次型

，z为n维向量正定的二次型：若对于任意

，有

；半正定的二次型：若对于任意

，有

；负定的二次型：若对于任意

，有

；半负定的二次型：若对于任意

，有

；不定二次型：

，有

，又

，有.二次型

为正定的充要条件是它的矩阵

的左上角各阶主子式都大于零.矩阵M的所有的特征值λi都是正的。函数Taylor展开相关数学基础知识参看：Lecture4MathematicalPreliminaries.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：利用极值条件解优化问题算法及相关概念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+1yesno3.终止条件b.d.a.c.xkxk+1xkxk+1x*一般而言，线性收敛速度是比较慢的，超线性收敛速度相对较快，而二阶收敛速度则相当快。如果一个算法具有超线性收敛速度，那么从计算的角度就可以认为是一个比较好的算法则称的收敛阶为。a.设算法A所得的点列为，如果b.4.收敛速度单谷函数（单峰函数）：设f(t)是定义在区间[a,b]上的一元函数，t*是f(t)在[a,b]上的全局极小点。如果在[a,b]上任取两个点t1<t2，当t2≤t*时，必有f(t1)>f(t2)；而当t1≥t*时，必有f(t1)<f(t2)，则称f(t)在[a,b]上是一个单谷函数。利用上述定理，可以缩短搜索区间设[a,b]是单谷函数f(t)的一个已知搜索区间，在[a,b]中任取两个点t1<t2：

如果f(t1)>f(t2)，则搜索区间可缩短为[t1,b]；

如果f(t1)≤f(t2)，则搜索区间可缩短为[a,t2]；无约束优化：线搜索法给定初始估计x(0)，设x(k)处有g(k)

≠0，则第k次迭代：(a)中的搜索方向有许多种不同的确定方法！根据某种模型函数确定x(k)处的搜索方向d(k)线搜索：确定的极小点置(b)中的确定步长：精确线搜索、非精确线搜索一维收索的方法很多，大致可以分为两类：一为不使用导数的方法，如进退法、Fibonacci法、0.618法；二为使用导数的方法，如对分法、牛顿法、抛物线法、三次插值法等。确定的极小点确定步长方法：分析法（最优步长）最优步长参看：Lecture5Non-linearProgramming.pdf节进退法求初始不确定区间

找三点使两端点的函数值大于中间点的函数值。思路：任取λ0，步长δ>0，取λ1=λ0+δ

，

1°若ф(λ0)<ф(λ1),令δ=2δ（步长加倍），λ2=λ0-δ

，

若ф(λ2)<ф(λ0),则令λ1=λ0

，λ0=λ2,重复1°

若ф(λ2)>ф(λ0),则停，a=λ2，b=λ1(图1)

2°若ф(λ0)>ф(λ1),令δ=2δ,λ2=λ1+δ

，若ф(λ2)<ф(λ1),则令λ0=λ1，λ1=λ2,重复2°

若ф(λ2)>ф(λ1),则停，a=λ0

，b=λ2(图2)λ2λ0λ1

向左找λ2

图1

λ0λ1λ2

向右找λ2

图2区间法参看：Lecture5Non-linearProgramming.pdf节0.618法和Fibonacci法0.618法和Fibonacci法都是分割方法，其基本思想是通过取试探点和函数值进行比较，使包含极小点的搜索区间不断缩短，当区间长度缩短到一定程度时，区间上各点的函数值均接近极小值，从而各点可以看作是极小点的近似。这类方法仅需计算函数值，不涉及导数，又称直接法，用途很广，尤其适用于非光滑及导数表达式复杂或写不出的情况。这类方法要求所考虑区间上的目标函数是单峰函数。如果不满足这个条件，可以考虑把区间分成若干个小区间，使得每个小区间上函数都是单峰，分别求极小值0.618法0.618法的基本思想是在搜索区间[a,b]上选取两个对称点λ1和λ2，且λ1<λ2，通过比较这两个点处的函数值φ(λ1)和φ(λ2)的大小来决定删除左半区间[a,λ1)，还是右半区间(λ2,b]。删除后的新区间长度是原来区间长度的0.618倍。新区间包含原区间中两个对称点中的一个，只要再选一个对称点，并利用这两个新对称点的函数值继续比较。重复这个过程，最后确定出极小点。黄金分割法（0.618法）（算法）初始[α,β],ε>0

λ=α+（1-t）(β-α)μ=α+t(β-α)β-α<ε?

STOP;λ*

=(α+β)/2yesФ（λ）-Ф（μ）>0?Noα=λ,λ=μ

μ=α+t(β-α)yesβ=μ,μ=λλ=α+（1-t）(β-α)Noαλμβμβ

αλμβαλ黄金分割法Matlab实现程序，参看golden_sect.pdfFibonacci法0.618法和Fibonnacci法都用于单谷函数，通过不断缩小搜索区间来获得极小点的近似值。0.618法的区间长度缩短比率是常数，而Fibonacci法的区间长度缩短比率有Fibonnacci数确定。当所需插入点的个数n充分大时，Fibonacci法的搜索区间缩短比率趋近于0.618法的缩短比率。最速下降法令

最速下降法参看：Lecture5Non-linearProgramming.pdf节最优解最优方向克服负梯度方向缺陷的途径用适当的正定矩阵（尺度矩阵）乘负梯度方向，其作用是对后者进行适当的旋转，以获得更好的方向牛顿法在局部用一个二次函数g(x)近似地代替目标函数，然后用g(x)的极小值作为f(x)的近似极小点牛顿方向牛顿法特点牛顿法的改进:修正牛顿法/阻尼牛顿法牛顿法参看：Lecture5Non-linearProgramming.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..拟牛顿法参看：“AnIntroductiontoOptimization”第11章，以及Lecture5Non-linearProgramming.pdf节共轭方向法最速下降法计算步骤简单，但收敛速度太慢；牛顿法和阻尼牛顿法收敛速度快，但要计算二阶偏导数矩阵及其逆矩阵，计算量太大。共轭方向法是兼有这两种方法优点的一类方法，比最速下降法的收敛速度快得多，但同时又避免了牛顿法所要求的Hessian矩阵的计算和求逆。共轭方向法，主要是其中的共轭梯度法（conjugategradientmethod）。1952年Hesteness和Stiefel为求解线性方程组提出了共轭梯度法，1964年Fletcher和Reeves提出了F-R共轭梯度法（简称FR法），用来求解无约束优化问题。共轭梯度法FR法的基本思想是把共轭性与最速下降法相结合，利用已知点处的梯度构造一组共轭方向，并沿着这组方向进行搜索，求出目标函数的极小值。其中,是对称正定矩阵,是常数.我们先讨论对于二次凸函数的共轭梯度法,然后再把这种方法推广到极小化一般函数的情形.

考虑问题求解方法：Fletcher-Reeves公式拟牛顿法参看：“AnIntroductiontoOptimization”第10章，以及Lecture5Non-linearProgramming.pdf节约束极值的最优化方法约束极值的条件的可叠加特性约束问题的最优性条件研究一点处的可行方向时，只需考虑在这一点的积极约束，可以暂时不管那些不积极约束可行方向Lagrange必要条件Lagrange定理实质是将等式约束问题转化为无约束问题来处理等式约束问题的最优性条件Lagrange充分条件不等式约束问题的最优性条件几何最优性条件几何最优性条件是不等式约束条件下极小点的必要条件可行方向代数条件下降方向代数条件可行下降方向Kuhn－Tucker条件Kuhn-Tucker定理简称为K-T定理，称所导出的有关极小点的必要条件为K-T条件，满足K-T条件的点称为K-T点KT条件几何解释Zoutendijk可行方向法可行方向法基本原理仅带线性约束的Zoutendijk可行方向法制约函数法这类方法是将约束条件随同一个参数并入到目标函数而形成一个新的无约束问题.当参数变化时就得到了一系列无约束问题，用求解这一系列无约束问题而得到的解来逼近原有约束问题的解，也称这类方法为序列无约束极小化方法.罚函数法(外点法)障碍函数法（内点法）罚函数法(外点法)障碍函数法（内点法）

用MATLAB软件求解,其输入格式如下:

1．x=quadprog(H,C,A,b);2．x=quadprog(H,C,A,b,Aeq,beq);3．x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB);4．x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB,X0);5．x=quadprog(H,C,A,b,Aeq,beq,VLB,VUB,X0,options);6．[x,fval]=quaprog(…);7．[x,fval,exitflag]=quaprog(…);8．[x,fval,exitflag,output]=quaprog(…);1．二次规划如何转化为二次型标准型参考“Lecture4MathematicalPreliminaries”0.1.4QuadraticForms一节例1min

f(x1,x2)=-2x1-6x2+x12-2x1x2+2x22s.t.

x1+x2≤2-x1+2x2≤2

x1≥0,x2≥0MATLAB（youh1）1．写成标准形式：2．输入命令：

H=[2-2;-24];c=[-2;-6];A=[11;-12];b=[2;2];Aeq=[];beq=[];VLB=[0;0];VUB=[];[x,z]=quadprog(H,c,A,b,Aeq,beq,VLB,VUB)3．运算结果为：

x=4/56/5z=-36/5s.t.

1．首先建立M文件fun.m,用来定义目标函数F（X）:functionf=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+3x26

s.t.

x1+4x25

x1,x20例22．先建立M-文件fun3．m:

functionf=fun3(x);f=-x(1)-2*x(2)+(1/2)*x(1)^2+(1/2)*x(2)^2MATLAB(youh2)3．再建立主程序youh2．m：

x0=[1;1];A=[23;14];b=[6;5];Aeq=[];beq=[];VLB=[0;0];VUB=[];[x,fval]=fmincon('fun3',x0,A,b,Aeq,beq,VLB,VUB)4．运算结果为：

x=0．76471．0588fval=-2．02941．先建立M文件fun4．m定义目标函数:

functionf=fun4(x);f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);

x1+x2=0