建模仿真与优化设计.pptVIP

建模仿真与优化设计 第一部分中北大学张保成 掌握优化问题的数学描述方法 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 连续可微 多局部极小 唯一极小 全局极小 搜索过程 最优点 11 初始点 11 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 11e1 然后以x11为新起点 沿第二坐标轴的方向e2作一维搜索 确定步长 21 得第一轮的第二个迭代点x21 x11 11e2第二轮迭代 需要x11 x01x12 x02 12e1x22 x12 22e2依次类推 不断迭代 目标函数值不断下降 最后逼近该目标函数的最优点 坐标轮换法 二 终止准则可以采用点距准则或者其它准则 注意 若采用点距准则或函数值准则 其中采用的点应该是一轮迭代的始点和终点 而不是某搜索方向的前后迭代点 坐标轮换法 三 坐标轮换法的流程图 入口 给定 x0 K 1 i 1 Xik x0 沿ei方向一维搜索求 ixik xi 1k ikeix xkf f x i n xnk x0k x xf 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作一维搜索 得x1x1 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 bgk 1 Gxk 1 b两式相减 得gk 1 gk G xk 1 xK 鲍威尔方法 因此有 Sj T gk 1 gk Sj TG 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 xn 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 xnkF2 F3x0k 1 xn 1kF2 F32 两式均不成立 则淘汰函数值下降最大的方向 并用第k环的新生方向补入k 1环基本方向组的最后 即k 1环的方向组为S1k S2k Sm 1k Sm 1k Snk Sn 1k k 1环的初始点取x0k 1 xkxk是第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 minf x k k g k min 据一元函数极值条件和多元复合函数求导公式 得 f x k k g k Tg k 0即 f x k 1 Tg 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 1 TGSk 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 1f 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如下 functionf fun0 x f 3 2 x 2 x 主程序为wliti2 m x fval fminbnd fun0 0 1 5 xmax xfmax 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 多元函数无约束优化问题 标准型为 minF 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 最速下降法 例3minf x 4x12 2x22 4x1x2 2x2 1 exp x1 1 编写M 文件fun1 m functionf 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 0000y 1 3029e 10 3 用fminsearch函数求解 输入命令 f 100 x 2 x 1 2 2 1 x 1 2 x fval exitflag output fminsearch f 1 22 运行结果 x 1 00001 0000fval 1 9151e 010exitflag 1output iterations 108funcCount 202algorithm Nelder Meadsimplexdirectsearch 4 用fminunc函数 1 建立M 文件fun2 mfunctionf 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 functionf 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 取M1 1 给定允许误差 令k 1 2 求无约束极值问题的最优解 设为Xk X Mk 即 3 若存在 使 则取Mk M 令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 二次规划 例1minf x1 x2 2x1 6x2 x12 2x1x2 2x22s t x1 x2 2 x1 2x2 2x1 0 x2 0 1 写成标准形式 2 输入命令 H 1 1 12 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 0 66671 3333z 8 2222 s 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法更新拉格朗日Hessian矩阵 3 fmincon函数可能会给出局部最优解 这与初值X0的选取有关 1 写成标准形式 s t 2x1 3x26s tx1 4x25x1 x20 例2 2 先建立M 文件fun3 m functionf 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 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 0294 1 先建立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 0s t 1 5 x1x2 x1 x20 x1x2 100 例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 11 beq 0 vlb vub x fval fmincon fun4 x0 A b Aeq beq vlb vub mycon 3 运算结果为 x 1 22501 2250fval 1 8951 例4 1 先建立M 文件fun m定义目标函数 functionf 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 00 VUB 510 x fval exitflag output fmincon fun x0 VLB VUB mycon2 4 运算结果为 x 4 00003 0000fval 11 0000exitflag 1output iterations 4funcCount 17stepsize 1algorithm 1x44char 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