建模仿真与优化设计.ppt

,建模仿真与优化设计 第一部分 中北大学 张保成,掌握优化问题的数学描述方法； 2. 熟练掌握常用优化算法。,一 优化模型的一般意义 二 非线性规划建模无约束最优化 非线性规划建模有约束非线性规划 四 连续结构体建模与优化设计,学 习 内 容,目 的,（一）优化模型的数学描述,一 优化模型的一般意义,“受约束于”之意,（二）优化模型的分类,1.根据是否存在约束条件 有约束问题和无约束问题。 2.根据设计变量的性质 静态问题和动态问题。,3.根据目标函数和约束条件表达式的性质 线性规划，非线性规划，二次规划，多目标规划等。,（1）非线性规划 目标函数和约束条件中，至少有一个非线性函数。,（2）线性规划（LP） 目标函数和所有的约束条件都是设计变量 的线性函数。,（3）二次规划问题 目标函数为二次函数，约束条件为线性约束,5. 根据变量具有确定值还是随机值 确定规划和随机规划。,4. 根据设计变量的允许值,整数规划（0-1规划）和实数规划。,（三）建立优化模型的一般步骤,1.确定设计变量和目标变量； 2.确定目标函数的表达式； 3.寻找约束条件。,二、优化求解的数学基础,函数的梯度,泰勒展开,二阶导数矩阵,矢量的概念、运算和点积,矩阵的运算和逆矩阵,（一）方向导数,和 分别是 函数 f( x1 , x2 )在x0点处沿坐标轴x1和x2方向的变化率,故函数 f( x1 , x2 )在x0(x10 , x20)点处沿某一方向S的变化率为： 称为该函数沿此方向的方向导数 偏导数可以看作是函数沿坐标轴 方向的方向导数，并有 （二）梯度 二元函数在点x0的梯度是由函数在该点的各一阶偏导数组成的向量。即 :,二、优化求解的数学基础,设S方向单位向量 则有 函数的梯度具有以下性质： （1）函数在一点的梯度是一个向量。梯度的方向是该点函数值上升最快的方向，与梯度相反的方向是该点函数值下降的最快的方向，梯度的大小就是它的模长。 （2）一点的梯度方向是与过该点的等值线或等值面的切线或切平面相垂直的方向，或者说是该点等值线或等值面的法线方向。 （3）梯度是函数在一点邻域内局部形态的描述。在一点上升得快的方向，离开该领域后就不一定上升得快，甚至可能下降。,二、优化求解的数学基础,（三）泰勒展开 为了便于数学问题的分析和求解，往往需要将一个复杂的非线性函数简化成线性函数或二次函数。简化的方法可以采用泰勒展开式。 由高等数学可知，一元函数 f(x)若在点x0的邻域内n阶可导，则函数可在该点邻域内作泰勒展开： 二元函数 f(x)在点x0(x10 , x20)也可以作泰勒展开，展开式一般取前三项，即：,二、优化求解的数学基础,将上式写为矩阵形式 其中，G(x0)称为函数 f( x1 , x2 )在点x0处的二阶导数矩阵或海赛矩阵。,二、优化求解的数学基础,在优化计算中，当某点附近的函数值采用泰勒展开式作近似表达时，研究该点邻域的极值问题需要分析二次型函数是否正定。当对任何非零向量x使 则二次型函数正定，G为正定矩阵,二、优化求解的数学基础,（四）二次函数 当将函数的泰勒展开式取到二次项时得到二次函数形式。优化计算经常把目标函数表示为二次函数以使问题分析得以简化。在线性代数中将二次齐次函数称作二次型，其矩阵形式 在优化计算中，当某点附近的函数值采用泰勒展开式作近似表达时，研究该点邻域的极值问题需要分析二次型函数是否正定。当对任何非零向量x使 则二次型函数正定，G为正定矩阵,二、优化求解的数学基础,对于一般二次函数 矩阵有正定和负定之分。对于所有非零向量： （1）若有 ，则称矩阵是正定的； （2）若有 ，则称矩阵是半正定的； （3）若有 ，则称矩阵是负定的； （4）若有 ，则称矩阵是半负定的； （5）若有 ，则称矩阵是不定的,二、优化求解的数学基础,可以证明，正定二次函数具有以下性质： （1）正定二次函数的等值线或等值面是一簇同心圆或同心椭球。椭圆簇或椭球簇的中心就是该二次函数的极小点。 （2）非正定二次函数在极小点附近的等值线或等值面近似于椭圆或椭球。,二、优化求解的数学基础,（五）无约束优化问题的极值条件 二次函数 f(x1 , x2)在x0取得极值的必要条件为 充分条件为：该点处海赛矩阵正定，即,二、优化求解的数学基础,正定：各阶主子式均大于零。,现代设计方法概论课程教案,（六） 下降迭代算法及其收敛性,无约束最优化问题求优过程的求解方法大致分为两类。 （1）解析法,二、优化求解的数学基础,现代设计方法概论课程教案,（2）数值迭代计算,（六） 下降迭代算法及其收敛性,二、优化求解的数学基础,三、优化设计的迭代计算,（一）优化问题的求解方法,1、优化问题的本质,优化问题的本质求极值的数学问题。,2、优化问题的求解方法,理论上,解析法,数值计算法,能求解吗？,由于实际优化数学模型的目标函数及约束函数往往是非线性的,解析法求解非常困难，甚至无法实现,（二）数值计算法的迭代方法,1、数值计算法的数学基础,计算方法,2、数值计算法的迭代方法,数值计算法可以较好地解决这类问题,三、优化设计的迭代计算,从目标函数出发,构造一种使目标函数值逐次下降的数值计算方法,利用计算机进行反复迭代运算,一步步搜索、调优逐步逼近函数极值点或最优点,直到满足一定的精度时终止迭代计算,最后所逼近的设计点即最优点，所得到的解即一定精度下的近似解,三、优化设计的迭代计算,（三）迭代计算的迭代过程,由选定的初始点x(0)出发,使满足,由于各设计点的函数值依次下降，可见迭代点不断向理论最优点逼近，最后可得到一定精度下的近似最优点，记作x*。,三、优化设计的迭代计算,迭代过程图示：,X(0),X(1),X(k),s(0),a(0),X(2),x*,三、优化设计的迭代计算,（四）迭代计算的终止准则,由于数值迭代是逐步逼近最优点而获得近似解的，所以要考虑优化问题解的收敛性及迭代过程的终止条件。,1、迭代的收敛性,指某种迭代程序产生的序列xk (k=0,1,2,)收敛于,2、通常采用的终止准则,（1）点距准则,|xk+1-xk| 1,相邻两个设计点的距离已达到充分小,或,两迭代点的坐标分量之差,三、优化设计的迭代计算,（2）函数下降量准则,相邻迭代点的函数值下降量达到充分小,或,相邻迭代点的函数值的相对值达到充分小,（3）梯度准则,无约束最优化,非线性规划建模,教学目的,教学内容,2、掌握用数学软件包求解无约束最优化问题。,1、了解无约束最优化基本算法。,1.无约束优化基本思想及基本算法。,3. 用MATLAB求解无约束优化问题。,2. MATLAB优化工具箱简介,无约束最优化问题,求解无约束最优化问题的的基本思想,*无约束最优化问题的基本算法,标准形式：,求解无约束最优化问题的基本思想,求解的基本思想 ( 以二元函数为例 ),5,3,1,连续可微,多局部极小,唯一极小 (全局极小),搜索过程,最优点 (1 1) 初始点 (-1 1),-1,1,4.00,-0.79,0.58,3.39,-0.53,0.23,2.60,-0.18,0.00,1.50,0.09,-0.03,0.98,0.37,0.11,0.47,0.59,0.33,0.20,0.80,0.63,0.05,0.95,0.90,0.003,0.99,0.99,1E-4,0.999,0.998,1E-5,0.9997,0.9998,1E-8,坐标轮换法,一、坐标轮换法的迭代过程 如图，以二次函数为例。,x2,x0 x01,x02 x21,x11,x12,x1,x21,坐标轮换法,任取一初始点x0作为第一轮的始点x01,先沿第一坐标轴的方向e1作一维搜索，用一维优化方法确定最优步长11，得第一轮的第一个迭代点x11=x01+ 11 e1，然后以x11为新起点，沿第二坐标轴的方向e2作一维搜索，确定步长21 ，得第一轮的第二个迭代点x21=x11+ 11 e2 第二轮迭代，需要 x11x01 x12 x02+ 12 e1 x22=x12+ 22 e2 依次类推，不断迭代，目标函数值不断下降，最后逼近该目标函数的最优点。,坐标轮换法,二、终止准则 可以采用点距准则或者其它准则。 注意： 若采用点距准则或函数值准则，其中采用的点应该是一轮迭代的始点和终点，而不是某搜索方向的前后迭代点。,坐标轮换法,三、坐标轮换法的流程图,入口,给定：x0，,K=1,i=1,Xik=x0,沿ei方向一维搜索求i xik=xi-1k+ ikei x=xk f=f(x),i=n?,|xnk-x0k|?,x*=x f*=f(x*),出口,i=i+1,x0=x0k,k=k+1,N,Y,N,Y,坐标轮换法,四、例题 五、小结 坐标轮换法程序简单，易于掌握。但是计算效率比较低，尤其是当优化问题的维数较高时更为严重。一般把此种方法应用于维数小于10的低维优化问题。 对于目标函数存在 “脊线”的情况，在脊线 的尖点处没有一个坐标方 向可以使函数值下降，只 有在锐角所包含的范围搜 索才可以达到函数值下降 的目的，故坐标轮换法对 此类函数会失效。,x2,x1,脊线,鲍威尔方法方向加速法,鲍威尔方法是直接利用函数值来构造共轭方向的一种共轭方向法。这种方法是在研究具有正定矩阵G的二次函数 的极小化问题时形成的。其基本思想是在不用导数的前提下，在迭代中逐次构造G的共轭方向。 一、共轭方向的概念 设G为一正定对称矩阵，若有一组非零向量S1，S2， ，Sn满足 SiTGSj=0 (i j )，则称这组向量关于矩阵G共轭。 共轭方向对于构造一种有效的算法是很重要的。以正定二元二次函数为例，我们进行探讨。,鲍威尔方法,正定的二元二次函数的等值线为一组椭圆，任选初始点x0沿某个下降方向S0作一维搜索，得x1 x1=x0+ 0S0此时，点x1的梯度必然与方向S0垂直，即有 f(x1)TS0=0S0与某一等值线相切于x1点。下一次的迭代，若选择负梯度方向为搜索方向，将产生锯齿现象。为避免锯齿的产生，我们取迭代方向S1直指极小点x*，如图所示。,x0,x*,x1,1S1,- f(x1),S1,0S0,鲍威尔方法,若选定这样的方向，则对于二元二次函数只需进行S0 和S1两次搜索就可以求得极小点x*，即 x*=x1+ 1S1由于f(x1)=Gx1+b，当x1 x*时 1S1 0，由于x*是f （x）的极小点，应满足极值必要条件，故f(x*)=Gx*+b=0，即 f(x*)=G（x1+ 1S1）+b= f(x1)+ 1GS1 =0 等式两端同乘以(S0)T,得（S0）TGS1=0 ，故两个向量S0、S1是G的共轭向量。 因此，若要使第二个迭代点成为正定二元二次函数的极小点，只要使两次一维搜索的方向S0、S1关于函数的二阶导数矩阵G共轭就可以了。,鲍威尔方法,二、共轭 方向的生成 设xk、xk+1为从不同点出发，沿同一方向Sj 进行一维搜索得到的两个极小点，如图所示。根据梯度和等值面的性质， Sj 和xk、xk+1两点处的梯度gk、gk+1之间存在如下关系 (Sj)Tgk=0 (Sj)Tgk+1=0又因为xk、xk+1两点处的梯度 可表示为 gk=Gxk+b gk+1=Gxk+1+b 两式相减，得 gk+1-gk=G(xk+1-xK),鲍威尔方法,因此有 (Sj)T (gk+1-gk)= (Sj)T G(xk+1-xK)=0 若取方向Sk= xk+1-xK，则Sk和Sj对G共轭。这说明只要沿方向Sj 分别对函数作两次一维搜索，得到两个极小点，则这两点的连线方向就是与Sj 共轭的方向。,鲍威尔方法,三、鲍威尔基本算法 如图所示，以三维二次目标函数的无约束优化问题为例。,x1,x3,x2,e1,e2,e3,e2,e3,e3,x01,x11,x21,x31,x1,x12,x22,x32,x13,x23,x33,x2,x3,S3,-S1,S2,S2,S1,鲍威尔方法,鲍威尔基本算法的步骤： 第一环基本方向组取单位坐标矢量系e1、 e2、 e3 、 en，沿这些方向依次作一维搜索，然后将始末两点相连作为新生方向，再沿新生方向作一维搜索，完成第一环的迭代。以后每环的基本方向组是将上环的第一个方向淘汰，上环的新生方向补入本环的最后而构成。n维目标函数完成n环的迭代过程称为一轮。从这一轮的终点出发沿新生方向搜索所得到的极小点，作为下一轮迭代的始点。这样就形成了算法的循环。,鲍威尔方法,鲍威尔基本算法的缺陷： 可能在某一环迭代中出现基本方向组为线性相关的矢量系的情况。如第k环中，产生新的方向： Sk=xnk-x0k=1kS1k+ 2kS2k + + nkSnk 式中， S1k、S2k 、 、 Snk为第k环基本方向组矢量， 1k 、 2k 、 、 nk为个方向的最优步长。 若在第k环的优化搜索过程中出现1k =0，则方向Sk表示为S2k 、 S3k 、 、 Snk的线性组合，以后的各次搜索将在降维的空间进行，无法得到n维空间的函数极小值，计算将失败。,鲍威尔方法,鲍威尔基本算法的退化,如图所示为一个三维优化问题的示例，设第一环中1 =0 ，则新生方向与e2 、e3共面，随后的各环方向组中，各矢量必在该平面内，使搜索局限于二维空间，不能得到最优解。,鲍威尔方法,四、鲍威尔修正算法 在某环已经取得的n+1各方向中，选取n个线性无关的并且共轭程度尽可能高的方向作为下一环的基本方向组。 鲍威尔修正算法的搜索方向的构造： 在第k环的搜索中， x0k 为初始点，搜索方向为S1k、S2k 、 、 Snk，产生的新方向为Sk ，此方向的极小点为xk。点 xn+1k=2xnk-x0k ， 为x0k对xnk的映射点。 计算x0k 、 x1k 、 、 xnk 、 xk、 x n+1k 各点的函数值，记作： F1=F(x0k) F2=F(xnk) F3=F(xn+1k) = F(xmk) -F(xm-1k) 是第k环方向组中，依次沿各方向搜索函数值下降最大值，即Smk方向函数下降最大。,鲍威尔方法,为了构造第k+1环基本方向组，采用如下判别式： 按照一下两种情况处理： 1、上式中至少一个不等式成立，则第k+1环的基本方向仍用老方向组S1k、S2k 、 、 Snk。 k+1的环初始点取 x0k+1=xnk F2F3 x0k+1=xn+1k F2F3 2、两式均不成立，则淘汰函数值下降最大的方向，并用第k环的新生方向补入k+1环基本方向组的最后，即k+1环的方向组为S1k、S2k 、 、 Sm-1k、Sm+1k 、 Snk 、 Sn+1k 。k+1环的初始点取 x0k+1=xk xk是第k环沿Sn+1k方向搜索的极小点。,鲍威尔方法,x0k,x1k,x2k,x3k,xm-1k,xmk,Snk,xnk,Smk,S3k,S2k,Sk,xk,xn+1k,F1,F2,F3,鲍威尔方法,鲍威尔算法的终止条件： | xk-x0k | 五、鲍威尔算法的迭代步骤及流程图 六、例题,梯度法,优化设计是追求目标函数值最小，因此，自然可以设想从某点出发，其搜索方向取该点的负梯度方向，使函数值在该点附近下降最快。这种方法也称为最速下降法。 一、基本原理梯度法的迭代公式为： x(k+1)=x(k)-(k)g(k)其中g(k)是函数F(x)在迭代点x(k)处的梯度f(xk) ， (k)一般采用一维搜索的最优步长，即 f(x(k+1)=f(x(k)-(k)g(k) =min f(x(k)-(k)g(k)=min()据一元函数极值条件和多元复合函数求导公式，得 ()= -( f(x(k)-(k)g(k)T g(k) =0即 （ f(x(k+1)T g(k) =0或 (g(k+1)Tg(k)=0,梯度法,此式表明，相邻的两个迭代点的梯度是彼此正交的。也即在梯度的迭代过程中，相邻的搜索方向相互垂直。梯度法向极小点的逼近路径是锯齿形路线，越接近极小点，锯齿越细，前进速度越慢。 这是因为，梯度是函数的局部性质，从局部上看，在该点附近函数的下降最快，但从总体上看则走了许多弯路，因此函数值的下降并不快。,梯度法,二、迭代终止条件 采用梯度准则： | g(k) | 三、迭代步骤（1）任选初始迭代点x(0)，选收敛精度 。（2）确定x(k)点的梯度（开始k=0)（3）判断是否满足终止条件| g(k) | ？若满足输出最优解，结束计算。否则转下步。（4）从x(k)点出发，沿-g(k)方向作一维搜索求最优步长(k)。得下一迭代点 x(k+1)=x(k)-(k)g(k) ，令k=k+1 返回步骤（2）。,梯度法,四、梯度法流程图,入口,给定： x(0)，,k=0,|g(k)| ？,x*=x(k) f*=f(x(k),出口,x(k)= x(0,计算：g(k),k=k+1,沿g(k)方向一维搜索， 求最优步长(k)。,N,Y,共轭梯度法,共轭梯度法是共轭方向法的一种，因为该方法中每一个共轭向量都是依赖于迭代点处的负梯度而构造出来的，所以称作共轭梯度法。 一、共轭梯度法的搜索方向 共轭梯度法的搜索方向采用梯度法基础上的共轭方向，如图所示，目标函数F(x)在迭代点xk+1处的负梯度为- f(xk+1)，该方向与前一搜索方向Sk互为正交，在此基础上构造一种具有较高收敛速度的算法，该算法的搜索方向要满足以下两个条件： （1）以xk+1点出发的搜索方向 Sk+1是- f(xk+1)与Sk的线性组合。即,xk,x*,xk+1,- f(xk+1),Sk+1,Sk,共轭梯度法,Sk+1 = - f(xk+1) + kSk （2）以Sk与Sk+1为基底的子空间中，矢量Sk与Sk+1相共轭，即满足 Sk+1T G Sk = 0 二、 k的确定 确定方法不作要求。记住 三、 共轭梯度法的算法 （1）选初始点x0和收敛精度。 （2）令k=0，计算S0 = - f(x0) 。 （3）沿Sk方向进行一维搜索求(k)，得 x(k+1)=x(k)+(k)S(k) （4）计算 f(xk+1) ，若| f(xk+1)| ，则终止迭代，取x*=xk+1；否则进行下一步。,共轭梯度法,（5）检查搜索次数，若k=n，则令x0=xk+1，转（2）， 否则，进行下一步。 （6）构造新的共轭方向 Sk+1 = - f(xk+1) + kSk 令k=k+1，转（3）,共轭梯度法,四、共轭梯度法流程图,入口,k=0，计算：- f(x0),| f(xk+1) | ？,出口,求(k) ，x(k+1)= x(k) +(k)S(k),计算： f(xk+1),x*=xk+1 f(x*)=f(xk+1),Y,N,给定： x(0)，,k n ？,x0=xk+1,N,Y,Sk+1 = - f(xk+1) + kSk,K=K+1,共轭梯度法,五、共轭梯度法的特点 共轭梯度法属于解析法，其算法需求一阶导数，所用公式及算法简单，所需存储量少。该方法以正定二次函数的共轭方向理论为基础，对二次型函数可以经过有限步达到极小点，所以具有二次收敛性。但是对于非二次型函数，以及在实际计算中由于计算机舍入误差的影响，虽然经过n次迭代，仍不能达到极小点，则通常以重置负梯度方向开始，搜索直至达到预定精度，其收敛速度也是较快的。 六、例题,牛顿法,牛顿法是求无约束最优解的一种古典解析算法。牛顿法可以分为原始牛顿法和阻尼牛顿法两种。实际中应用较多的是阻尼牛顿法。,原始牛顿法,一、原始牛顿法的基本思想 在第k次迭代的迭代点x(k)邻域内，用一个二次函数去近似代替原目标函数f(x)，然后求出该二次函数的极小点作为对原目标函数求优的下一个迭代点，依次类推，通过多次重复迭代，是迭代点逐步逼近原目标函数的极小点。如图所示。,F(x),1(x),0(x),x0,x1,x2,x*,原始牛顿法,二、原始牛顿法的迭代公式 设目标函数f(x)具有连续的一、二阶导数，在x(k)点邻域内取f(x)的二次泰勒多项式作近似式，即取逼近函数(x)为设xk+1为(x)极小点，根据极值的必要条件，应有(xk+1)=0，即 f(xk)+ G(xk)x=0 又 x= xk+1 - xk 得 xk+1 = xk- G(xk)-1 f(xk) 即牛顿法迭代公式，方向- G(xk)-1 f(xk)称为牛顿方向,原始牛顿法,三、原始牛顿法的特点 若用原始牛顿法求某二次目标函数的最优解，则构造的逼近函数与原目标函数是完全相同的二次式，其等值线完全重合，故从任一点出发，一定可以一次达到目标函数的极小点。 因此，牛顿法是具有二次收敛性的算法。其优点是：对于二次正定函数，迭代一次即可以得到最优解，对于非二次函数，若函数二次性较强或迭代点已经进入最优点的较小邻域，则收敛速度也很快。 原始牛顿法的缺点是：由于迭代点的位置是按照极值条件确定的，并未沿函数值下降方向搜索，因此，对于非二次函数，有时会使函数值上升，即 f(xk+1) f(xk)，而使计算失败。,阻尼牛顿法,一、对原始牛顿法的改进 为解决原始牛顿法的不足，加入搜索步长(k)因此，迭代公式变为： xk+1 = xk - (k)G(xk)-1 f(xk) 这就是阻尼牛顿法的迭代公式，最优步长(k)也称为阻尼因子，是沿牛顿方向一维搜索得到的最优步长。 二、阻尼牛顿法的迭代步骤 （1）给定初始点和收敛精度 （2）计算 f(xk) 、 G(xk)、 G(xk)-1 （3）求xk+1 = xk - (k)G(xk)-1 f(xk) （4）检查收敛精度，若| xk+1- xk |则x*=xk+1，停止，否则 k=k+1，返回（2）继续,阻尼牛顿法,三、阻尼牛顿法的特点 优点： 由于阻尼牛顿法每次迭代都在牛顿方向进行一维搜索，避免了迭代后函数值上升的现象，从而保持了牛顿法的二次收敛性，而对初始点的选择没有苛刻的要求。 缺点： 1、对目标函数要求苛刻，要求函数具有连续的一、二阶导数；为保证函数的稳定下降，海赛矩阵必须正定；为求逆阵要求海赛矩阵非奇异。 2、计算复杂且计算量大，存储量大,变尺度法,Matlab优化工具箱简介,1.MATLAB求解优化问题的主要函数,2. 优化函数的输入变量,使用优化函数或优化工具箱中其它优化函数时, 输入变量见下表:,3. 优化函数的输出变量下表:,4控制参数options的设置,(3) MaxIter: 允许进行迭代的最大次数,取值为正整数.,Options中常用的几个参数的名称、含义、取值如下:,(1)Display: 显示水平.取值为off时,不显示输出; 取值为iter时,显示每次迭代的信息;取值为final时,显示最终结果.默认值为final.,(2)MaxFunEvals: 允许进行函数评价的最大次数,取值为正整数.,例：opts=optimset(Display,iter,TolFun,1e-8) 该语句创建一个称为opts的优化选项结构,其中显示参数设为iter, TolFun参数设为1e-8.,控制参数options可以通过函数optimset创建或修改。命令的格式如下：,(1) options=optimset(optimfun) 创建一个含有所有参数名,并与优化函数optimfun相关的默认值的选项结构options.,（2）options=optimset(param1,value1,param2,value2,.) 创建一个名称为options的优化选项参数,其中指定的参数具有指定值,所有未指定的参数取默认值.,(3)options=optimset(oldops,param1,value1,param2, value2,.) 创建名称为oldops的参数的拷贝,用指定的参数值修改oldops中相应的参数.,用Matlab解无约束优化问题,其中（3）、（4）、（5）的等式右边可选用（1）或（2）的等式右边。 函数fminbnd的算法基于黄金分割法和二次插值法，它要求目标函数必须是连续函数，并可能只给出局部最优解。,常用格式如下： （1）x= fminbnd (fun,x1,x2) （2）x= fminbnd (fun,x1,x2 ，options) （3）x，fval= fminbnd（.） （4）x，fval，exitflag= fminbnd（.） （5）x，fval，exitflag，output= fminbnd（.）,主程序为wliti1.m: f=2*exp(-x).*sin(x); fplot(f,0,8); %作图语句 xmin,ymin=fminbnd (f, 0,8) f1=-2*exp(-x).*sin(x); xmax,ymax=fminbnd (f1, 0,8),例2 对边长为3米的正方形铁板，在四个角剪去相等的正方形以制成方形无盖水槽，问如何剪法使水槽的容积最大？,解,先编写M文件fun0.m如下: function f=fun0(x) f=-(3-2*x).2*x;,主程序为wliti2.m: x,fval=fminbnd(fun0,0,1.5); xmax=x fmax=-fval,运算结果为: xmax = 0.5000,fmax =2.0000.即剪掉的正方形的边长为0.5米时水槽的容积最大,最大容积为2立方米.,命令格式为: （1）x= fminunc（fun,X0 ）；或x=fminsearch（fun,X0 ） （2）x= fminunc（fun,X0 ，options）； 或x=fminsearch（fun,X0 ，options） （3）x，fval= fminunc（.）； 或x，fval= fminsearch（.） （4）x，fval，exitflag= fminunc（.）； 或x，fval，exitflag= fminsearch （5）x，fval，exitflag，output= fminunc（.）； 或x，fval，exitflag，output= fminsearch（.）,2、多元函数无约束优化问题,标准型为：min F(X),3 fminunc为中型优化算法的步长一维搜索提供了两种算法， 由options中参数LineSearchType控制： LineSearchType=quadcubic(缺省值)，混合的二次和三 次多项式插值； LineSearchType=cubicpoly，三次多项式插,使用fminunc和 fminsearch可能会得到局部最优解.,说明:,fminsearch是用单纯形法寻优. fminunc的算法见以下几点说明：,1 fminunc为无约束优化提供了大型优化和中型优化算法。由options中的参数LargeScale控制： LargeScale=on(默认值),使用大型算法 LargeScale=off(默认值),使用中型算法,2 fminunc为中型优化算法的搜索方向提供了4种算法，由 options中的参数HessUpdate控制： HessUpdate=bfgs（默认值），拟牛顿法的BFGS公式； HessUpdate=dfp，拟牛顿法的DFP公式； HessUpdate=steepdesc，最速下降法,例3 min f(x)=(4x12+2x22+4x1x2+2x2+1)*exp(x1),1、编写M-文件 fun1.m: function f = fun1 (x) f = exp(x(1)*(4*x(1)2+2*x(2)2+4*x(1)*x(2)+2*x(2)+1); 2、输入M文件wliti3.m如下: x0 = -1, 1; x=fminunc(fun1,x0); y=fun1(x),3、运行结果: x= 0.5000 -1.0000 y = 1.3029e-10,3.用fminsearch函数求解,输入命令: f=100*(x(2)-x(1)2)2+(1-x(1)2; x,fval,exitflag,output=fminsearch(f, -1.2 2),运行结果: x =1.0000 1.0000 fval =1.9151e-010 exitflag = 1 output = iterations: 108 funcCount: 202 algorithm: Nelder-Mead simplex direct search,4. 用fminunc 函数,(1)建立M-文件fun2.m function f=fun2(x) f=100*(x(2)-x(1)2)2+(1-x(1)2,(2)主程序wliti44.m,Rosenbrock函数不同算法的计算结果,可以看出，最速下降法的结果最差.因为最速下降法特别不适合于从一狭长通道到达最优解的情况.,例5 产销量的最佳安排 某厂生产一种产品有甲、乙两个牌号，讨论在产销平衡的情况下如何确定各自的产量，使总利润最大. 所谓产销平衡指工厂的产量等于市场上的销量.,基本假设,1价格与销量成线性关系,2成本与产量成负指数关系,模型建立,若根据大量的统计数据,求出系数b1=100,a11=1,a12=0.1,b2=280, a21=0.2,a22=2,r1=30,1=0.015,c1=20, r2=100,2=0.02,c2=30,则 问题转化为无约束优化问题：求甲,乙两个牌号的产量x1，x2，使 总利润z最大.,为简化模型,先忽略成本,并令a12=0,a21=0,问题转化为求: z1 = ( b1 - a11x1 ) x1 + ( b2 - a22x2 ) x2 的极值. 显然其解为x1 = b1/2a11 = 50, x2 = b2/2a22 = 70, 我们把它作为原问题的初始值.,总利润为： z(x1,x2)=(p1-q1)x1+(p2-q2)x2,模型求解,1.建立M-文件fun.m: function f = fun(x) y1=(100-x(1)- 0.1*x(2)-(30*exp(-0.015*x(1)+20)*x(1); y2=(280-0.2*x(1)- 2*x(2)-(100*exp(-0.02*x(2)+30)*x(2); f=-y1-y2;,2.输入命令: x0=50,70; x=fminunc(fun,x0), z=fun(x),3.计算结果: x=23.9025, 62.4977, z=6.4135e+003 即甲的产量为23.9025,乙的产量为62.4977,最大利润为6413.5.,非线性规划建模,有约束非线性规划,教学目的,教学内容,2、掌握用数学软件求解优化问题。,1、直观了解非线性规划的基本内容。,1.非线性规划的基本理论。,2. 用数学软件求解非线性规划。,3. 钢管订购及运输优化模型,*非线性规划的基本解法,非线性规划,非线性规划的基本概念,定义 如果目标函数或约束条件中至少有一个是非线性函数时的最优化问题就叫做非线性规划问题,非线性规划的基本概念,一般形式: （1） 其中 ， 是定义在 En 上的实值函数，简记:,其它情况: 求目标函数的最大值或约束条件为小于等于零的情况，都可通过取其相反数化为上述一般形式,定义1 把满足问题（1）中条件的解 称为可行解（或可行点），所有可行点的集合称为可行集（或可行域）记为D即 问题(1)可简记为 ,定义2 对于问题(1)，设 ， 若存在 ,使得对一切 ，且 ，都有 ，则称X*是f(X)在D上的局部极小值点（局部最优解）特别地当 时,若 ，则称X*是f(X)在D上的严格局部极小值点（严格局部最优解）,定义3 对于问题(1),设 ,对任意的 ,都有 则称X*是f(X)在D上的全局极小值点（全局最优解）特别地当 时，若 ，则称X*是f(X)在D上的严格全局极小值点（严格全局最优解）,非线性规划的基本解法,SUTM外点法,SUTM内点法（障碍罚函数法）,1. 罚函数法,2. 近似规划法,罚函数法,罚函数法基本思想是通过构造罚函数把约束问题转化为一系列无约束最优化问题，进而用无约束最优化方法去求解这类方法称为序列无约束最小化方法.简称为SUMT法 其一为SUMT外点法,其二为SUMT内点法,其中T(X,M)称为罚函数，M称为罚因子，带M的项称为罚项，这里的罚函数只对不满足约束条件的点实行惩罚:当 时，满足各 ，故罚项=0，不受惩罚当 时，必有 的约束条件，故罚项0，要受惩罚,SUTM外点法,罚函数法的缺点是：每个近似最优解Xk往往不是容许解，而只能近似满足约束，在实际问题中这种结果可能不能使用；在解一系列无约束问题中，计算量太大，特别是随着Mk的增大，可能导致错误,1、任意给定初始点X0，取M11，给定允许误差 ，令k=1； 2、求无约束极值问题 的最优解，设为Xk=X(Mk)，即 ； 3、若存在 ，使 ，则取MkM( )令k=k+1返回（2），否则，停止迭代得最优解 . 计算时也可将收敛性判别准则 改为 .,SUTM外点法(罚函数法)的迭代步骤,SUTM内点法（障碍函数法）,内点法的迭代步骤,近似规划法的基本思想：将问题(3)中的目标函数 和约束条件 近似为线性函数，并对变量的取值范围加以限制，从而得到一个近似线性规划问题,再用单纯形法求解之，把其符合原始条件的最优解作为(3)的解的近似,近似规划法,每得到一个近似解后，都从这点出发，重复以上步骤,这样，通过求解一系列线性规划问题，产生一个由线性规划最优解组成的序列，经验表明，这样的序列往往收敛于非线性规划问题的解。,近似规划法的算法步骤如下,用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、二次规划,例1 min f(x1,x2)=-2x1-6x2+x12-2x1x2+2x22 s.t. x1+x22 -x1+2x22 x10, x20,1、写成标准形式：,2、 输入命令： H=1 -1; -1 2; 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 =0.6667 1.3333 z = -8.2222,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法更新拉格朗日Hessian矩阵。 3 fmincon函数可能会给出局部最优解，这与初值X0的选取有关。,1、写成标准形式： s.t.,2x1+3x2 6 s.t x1+4x2 5 x1,x2 0,例2,2、先建立M-文件 fun3.m: function f=fun3(x); f=-x(1)-2*x(2)+(1/2)*x(1)2+(1/2)*x(2)2,3、再建立主程序youh2.m： 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 = 0.7647 1.0588 fval = -2.0294,1先建立M文件 fun4.m,定义目标函数: function f=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 s.t. 1.5+x1x2 - x1 - x2 0 -x1x2 10 0,例3,2再建立M文件mycon.m定义非线性约束： function g,ceq=mycon(x) g=x(1)+x(2);1.5+x(1)*x(2)-x(1)-x(2);-x(1)*x(2)-10;,3主程序youh3.m为: x0=-1;1; A=;b=; Aeq=1 1;beq=0; vlb=;vub=; x,fval=fmincon(fun4,x0,A,b,Aeq,beq,vlb,vub,mycon),3. 运算结果为： x = -1.2250 1.2250 fval = 1.8951,例4,1先建立M-文件fun.m定义目标函数: function f=fun(x); f=-2*x(1)-x(2);,2再建立M文件mycon2.m定义非线性约束： function g,ceq=mycon2(x) g=x(1)2+x(2)2-25;x(1)2-x(2)2-7;,3. 主程序fxx.m为: x0=3;2.5; VLB=0 0;VUB=5 10; x,fval,exitflag,output =fmincon(fun,x0,VLB,VUB,mycon2),4. 运算结果为: x = 4.0000 3.0000 fval =-11.0000 exitflag = 1 output = iterations: 4 funcCount: 17 stepsize: 1 algorithm: 1x44 char firstorderopt: cgiterations: ,应用实例： 供应与选址,某公司有6个建筑工地要开工，每个工地的位置（用平面坐标系a，b表示，距离单位：千米 ）及水泥日用量d(吨)由下表给出。目前有两个临时料场位于A(5,1)，B(2,7)，日储量各有20吨。假设从料场到工地之间均有直线道路相连。 （1）试制定每天的供应计划，即从A，B两料场分别向各工地运送多少吨水泥，使总的吨千米数最小。 （2）为了进一步减少吨千米数，打算舍弃两个临时料场，改建两个新的，日储量各为20吨，问应建在何处，节省的吨千米数有多大？,（一）、建立模型,记工地的位置为(ai，bi)，水泥日用量为di，i=1,6;料场位置为(xj，yj)，日储量为ej，j=1,2；从料场j向工地i的运送量为Xij。,当用临时料场时决策变量为：Xij， 当不用临时料场时决策变量为：Xij，xj，yj。,（二）使用临时料场的情形,使用两个临时料场A(5,1)，B(2,7).求从料场j向工地i的运送量为Xij，在各工地用量必须满足和各料场运送量不超过日储量的条件下，使总的吨千米数最小，这是线性规划问题. 线性规划模型为：,设X11=X1, X21= X 2, X31= X 3, X41= X 4, X51= X 5, X61= X 6 X12= X 7, X22= X 8, X32= X 9, X42= X 10, X52= X 11, X62= X 12 编写程序gying1.m,计算结果为：,x = 3.0000 5.0000 0.0000 7.0000 0.0000 1.0000 0.0000 0.0000 4.0000 0.0000 6.0000 10.0000 fval = 136.2275,（三）改建两个新料场的情形,改建两个新料场,要同时确定料场的位置(xj,yj)和运送量Xij，在同样条件下使总吨千米数最小.这是非线性规划问题。非线性规划模型为：,设 X11=X1, X21= X 2, X31= X 3, X41= X 4, X51= X 5, X61= X 6 X12= X 7, X22= X 8, X32= X 9, X42= X 10, X52= X 11, X62= X 12 x1=X13, y1=X14, x2=X15, y2=X16,（1）先编写M文件liaoch.m定义目标函数。,(2) 取初值为线性规划的计算结果及临时料场的坐标: x0=3 5 0 7 0 1 0 0 4 0 6 10 5 1 2 7; 编写主程序gying2.m.,(3) 计算结果为：,x= 3.0000 5.0000 0.0707 7.0000 0 0.9293 0 0 3.9293 0 6.0000 10.0707 6.3875 4.3943 5.7511 7.1867 fval = 105.4626 exitflag = 1,(4) 若修改主程序gying2.m, 取初值为上面的计算结果: x0= 3.0000 5.0000 0.0707 7.0000 0 0.9293 0 0 3.9293 0 6.0000 10.0707 6.3875 4.3943 5.7511 7.1867,得结果为: x=3.0000 5.0000 0.3094 7.0000 0.0108 0.6798 0 0 3.6906 0 5.9892 10.3202 5.5369 4.9194 5.8291 7.2852 fval =103.4760 exitflag = 1,总的吨千米数比上面结果略优.,(5) 若再取刚得出的结果为初值, 却计算不出最优解.,(6) 若取初值为: x0=3 5 4 7 1 0 0 0 0 0 5 11 5.6348 4.8687 7.2479 7.7499, 则计算结果为: x=3.0000 5.0000 4.0000 7.0000 1.0000 0 0 0 0 0 5.0000 11.0000 5.6959 4.9285 7.2500 7.7500 fval =89.8835 exitflag = 1 总的吨千米数89.8835比上面结果更好.,通过此例可看出fmincon函数在选取初值上的重要性.,连续结构体建模与优化设计,连续体形状优化例一,连续体形状优化例二,注：并非结 构变形动画， 而是结构受 力状况下自 动寻优过程 注意：网格 变化过程 （网格拓扑 不改变，类 似蜘蛛网受 力变化）,离散模型敏度分析方法,对E求导,？,对E求导,离散法主要为有限 元程序开发者应用,代入已知得所求,连续模型敏度分析方法,敏度表达为位移、 应力等和区域形状 变化量的边界积分 或区域积分的形式,图 连续体变形过程的映射示意图,网格均匀化(Mesh Smoothing),启发式算法,启发式算法(heuristic algorithm)是相对于最优化算法提出的,启发式算法定义：一种技术。这种技术使得在可接受的计算 的计算费用内去寻找最好的解，但不一定能保证所得解的可 行性和最优性，甚至在多数情况下，无法阐述所得解同最优 解的近似程度,结构优化之满应力法,满应力法理论上受启发于生物进化方式,满应力法是一种经过实践检验的较好的优化方法,满应力法数学模型可以看成是最优问题的近似模型,连续体拓扑优化例一,注意：图片中的孤立块和锯齿边界形状,连续体拓扑优化例二,注意：图片中的棋盘图案,连续体拓扑优化例三,注意：拓扑优化是 一种概念设计,映射反演方法例一,微分方程,代数方程,求解代数方程,方程解,拉氏变换,反拉氏变换,RMI方法使用条件：,映射须是两类数学对象之间 的一一对应关系； 映射须是可定映的； 逆映射（反演）具有能行性,授课结束 谢 谢 大 家！,