1、 2012 年 最优化理论与方法 上机报告 学院： 电子工程 专业： 智能科学与技术 文章详细介绍了 使用 matlab 实现无约束和约束最优化问题 的方法和程序 ， 及使用方法 ， 包含黄金分割法、最速下降法、牛顿法、外点法、内点法 ， 例子丰富 . 最优化理论与方法上机报告 1 / 25 最优化理论与方法上机报告 虽然 上机题仅有两个 ， 但我又增 加了 一 个 ， 以验证程序的正确性 ！ 报告 按照优化方法 ： 黄金分割法、最速下降法、牛顿法、外点法、内点法 的顺序 叙述 ， 分为两篇 ， 第一篇为 无约束优化 ， 第二篇为约束优化 . 第一篇 无约束优化 1. 黄金分割法 黄金分割法又
2、称 0.618 法 ， 其基本思想是试探点函数值的比较 ， 使包含极小点的搜索区间不断缩小 .该方法仅需计算函数值 ， 适用范围广 ， 使用方便 .由于该方法思想 简单 ， 在此不给出具体算法流程 ， 仅介绍思想和 matlab 实现的程序 . 1.1 思想 在搜索区间 , ab 内插入点 12,xx， 其中 ， 120 .3 8 2 ( )0 .6 1 8 ( )x a b ax a b a ， 把 , ab 区间分为了三段 ， 比较 12( ) ( )f x f x与 的大小 ， 若 12( ) ,f x f x 则 * 1 , ,x x b 去掉 1, ;ax 若 12( ) ,f x
3、f x 则 * 2 , ,x a x 去 掉 2 , xb， 然后在余下的区间上根据对称原则再计算一个对称点的函数值 ， 再重复上述操作 ， 直到满足精度 . 1.2 实现 及结果 采用 matlab 2011a 实现 ， 命名为 golden， 具体程序见附录 golden.m 文件 . 该程序的使用方法 ， 请参见 附录中 golden.m文件中 绿色 部分 ， 下面给出 程序运行结果 ， 其中目标函数为： 211( ) 6f x x x ， 精度采用默认： 610 . 最优化理论与方法上机报告 2 / 25 图 1-1-1 黄金分割法运行结果 由 211( ) 6f x x x ， 求导
4、得 1( ) 2 1,f x x 易知 1 0.5x 为其极小点 ， 可见程序正确 . 2. 最速下降法 最速下降法 是求解无约束优化问题最简单和最古老的方法之一 ， 时至今日已不再具有实用性 ， 但它是研究其它无约束优化算法的基础 . 最速下降法适合解决下述优化问题 ： min ( )nxRfx 其中 ()fx是 nR 上的连续可谓函数 . 2.1 思想 无约束优化问题下降类算法的一般框架是 ， 用不同的方式确定搜索方向和搜索步长 ， 就得到不同的算法 ， 最速下降法采用负梯度方向 ()kkd f x 最优化理论与方法上机报告 3 / 25 作为搜索方向 ， 由于该方向是 函数值 下降最快的
5、方向 ， 故称为最速下降法 . 2.2 算法 最速下降法的算法如下 : Step1 取初始点 0 ,nxR 容许误差 01 ， 令 k:=0. Step2 计算 ()kkd f x . Step3 若 kd ， 则停止迭代 ， 取 * kxx ， 否则 ， 转下一步 . Step4 求 0k ， 使0( ) m i n ( )k k k k kf x d f x d . Step5 令 1 , : 1k k kx x d k k 转 Step2. 2.3 实现及结果 采用 matlab 2011a 实现 ， 具体程序见附录 fast.m 文件 . 下面以上机题第一题为例介绍 具体使用方法 :
6、上机题 1的 优化目标函数为 2 2 4 41 2 3 4 2 3 1 4( ) ( 1 0 ) 5 ( ) ( 2 ) 1 0 ( ) ,f x x x x x x x x x 首先分别建立两个文件即目标函数文件 func.m和 对应梯度文件 gfunc.m内容分别如下 : func.m文件 function f = func( x ) % -最速下降法、牛顿法用 - f=(x(1)+10*x(2)2+5*(x(3)-x(4)2+(x(2)-2*x(3)4+10*(x(1)-x(4)4; gfunc.m文件 function gf = gfunc( x ) %目标函数的梯度 gf=2*(x(
7、1)+10*x(2)+40*(x(1)-x(4)3 20*(x(1)+10*x(2)+4*(x(2)-2*x(3)3 10*(x(3)-x(4)-8*(x(2)-2*x(3)3 -10*(x(3)-x(4)-40*(x(1)-x(4)3; 最优化理论与方法上机报告 4 / 25 然后 在命令窗口输入 :X0=1 2 3 4， 再 输入 :fast(func， gfunc， X0);即可 (加粗部分 ， X0为初始点 ， 任意 )， 运行结果如下 : (a)精度较 高 时的情况 (b)精度较低时的情况 图 1-2-1 最速下降法 运行结果 由于最速下降法收敛速度 很缓慢 ， 要达到较高精度 ，
8、迭代的次数要很大 ， 正如上图 (a)所示 ， 当设置的精度较高时 ( 610 )， 迭代了 50000 次也未达到 ， 实际上需要 366304 次才可以 ，给出的最优解 为 * 0 . 0 0 4 9 , 0 . 0 0 0 5 , 0 . 0 0 2 4 , 0 . 0 0 2 4 Tx .但精度不是很高时 ， 得到的结果也是令人满意的 ， 如图 (b)所示 . 2.4 体会 Matlab 的运算速度不太令人满意，当然有时程序的运算速度是和算法及编程实现者的水平密不可分的！ 在最速下降法的实现过程中 ， 刚开始使用 subs 函数计算函数值 和一般的函数调用方式 会 大大降低程序的速度
9、， 几经修改后 ， 得到了上面的程序 ， 原来的程序迭代 5000 次需要 30 分钟不 只 ， 改进后的程序仅需不到 10 秒的时间 ， 即使是要达到 610 的精度，也只 需要 4 分钟， 可见效果是 相当 可观的 ! 这其中 ， 应用了函数句柄和匿名函数 ， 其实它们就是 为了使 feval 及借助于它的泛函指令工作更可靠；使“函数调用”像“变量调用”一样方便灵活；提高函数调用速度，特别在反复调用情况下更显效率；提高软件重用性，扩大子函数和私用函数的可调用范围； 迅速获得同名重载函数的位置、类型信息 . 最优化理论与方法上机报告 5 / 25 函数句柄和匿名函数 的 操作符 均 为 ,匿
10、名函数是函数句柄的一种特殊用法，这里所得到的函数句柄变量不指向特定的函数（即不指向函数 M文件中的函数名 ），而是指向一个函数表达式（具体表达式） .在命令窗口中输入 func就是用了函数句柄将其传给 objf变量 ,而 step=golden(x)objf(X0+x*d)， 0， 1，epsilon);则是 通过 匿名函数 将函数传给 golden程序 ,此两处的改变就大大提高了程序的运行速度 ! 3. 牛顿法 牛顿法也是一种经典的无约束优化方法 ， 并且其收敛速度快 ， 具有自适应的特点 ， 至今仍受到人们的青睐 . 3.1 思想 牛顿法的基本思想是用迭代点 kx 处的一阶导数（梯度）和二
11、阶导数（ Hesse阵）对目标函数进行二次函数近似，然后把二次模型的极小点作 为新的迭代点，并不断重复这一过程，直至达到精度 . 3.2 算法 牛顿法 的算法如下 : Step1 取初始点 0 ,nxR 容许误差 01 ，令 k:=0. Step2 计算 ()kkg f x . Step3 若 kg ，则停止迭代，取 * kxx ， 否则 ，转下一步 . Step4 计算 2 ( ),kkH f x 并由 ,kkHd g 求 得 k k kd H g . Step5 令 1 , : 1k k kx x d k k 转 Step2. 3.3 实现及结果 采用 matlab 2011a 实现， 具
12、体程序见附录 Newton.m 文件 . 下面以上机题第一题为例介绍具体使用方法 : 上机题 1的优化目标函数为 2 2 4 41 2 3 4 2 3 1 4( ) ( 1 0 ) 5 ( ) ( 2 ) 1 0 ( ) ,f x x x x x x x x x 仍然使用 最速下降法中建立的 func.m文件 ,在命令窗口中输 入 (加粗 部分 )：最优化理论与方法上机报告 6 / 25 X0=1 2 3 4;Newton(func,X0);即可按默认的 精度 610 计算 ，若想改变精度，如 310 则输入 ： Newton(func,X0,0.001);对应的结果 分别 如下 图的 (a)
13、、 (b)所示 : (a)采用默认 精度 (b)精度 为 0.001 图 1-3-1 牛顿法 运行结果 牛顿法收敛速度快，正 如上图（ a）所示，达到 610 精度只需 18 步 ，而最速下降法却需要 366304 步，两者的收敛速度由此可见一斑 ，可谓有天壤之别 ！ 3.4 体会 由于牛顿法收敛速度快， 且人工求目标函数的梯度、 Hesse 矩阵，有时较为繁杂， 故算法实现时 ，实现了自动求梯度和 Hesse 矩阵 ，这样便在 程序的易用性和时间复杂度两者间 得到了很好的折中 . 第 二 篇 约束优化 求解约束优化问题的经典算法是罚函数法， 其基本思想是：根据约束条件的特点将其转化为某种惩罚
14、函数加到目标函数中去，从而将约束优化问题 转化为一系列的无约束优化问题来求解。 1. 外点法 外点法适合求解 如下约束优化问题 ： m i n ( ) , ,. . ( ) 0 , 1 , , ( ) 0 , 1 , , .niif x x Rs t h x i E pg x i I m ，最优化理论与方法上机报告 7 / 25 但本次实验仅需求解 如下约束优化问题 ： m i n ( ) , ,. . ( ) 0 , 1 , , .nif x x Rs t g x i I m 而且，等式约束也可以转化为不等式约束！ 记可行域为 | ( ) 0 ( ) .n iD x R g x i I 1.
15、1 思想 方法 构造罚函数 21( ) m in 0 , ( ) ,miiP x g x 和增广目标函数 ( , ) ( ) ( ) ,T x M f x M P x 其中 0M 为惩罚因子 ，不难发现当 xD 时，即 x 为可行点时 ，( , ) ( ),T x M f x 此时 目标函数没有受到惩罚 ； 当 xD 时， ( , ) ( ),T x M f x 目标函数受到惩罚 ， M 越大 ，受到的惩罚越重 ， 当 M 足够大时， 要使 ()fx取得极小 ，罚函数 ()Px 应充分小， 从而 ( , )TxM 的极小点充分逼近可行域 D， 这样求解一般约束优化问题就化为求解一系列无约束的优
16、化问题 min ( , )T x M . 1.2 算法 外点法 的算法如下 : Step1 取初始点 0 ,nxR 容许误差 01 ， 惩罚因子 1 0 , 1,Mr令 k:=1. Step2 以 1kx 为初始点求解子问题 m i n ( , ) ( ) ( )n kkxR T x M f x M P x ，设其极小点为 kx . Step3 若 (1 ),j j m 使 ( ) ,jkgx 则 令 1 , 1,kkM rM k k 转 Step2,否则 停止迭代，取 * kxx . 1.3 实现及结果 最优化理论与方法上机报告 8 / 25 采用 matlab 2011a 实现， 无约束优
17、化采用牛顿法， 即直接调用第一篇中的Newton 函数 ， 具体程序见附录 EPM.m 文件 . 程序中 的 Min=Newton(augf,Xk,epsilon);是无约束优化函数的调用 ，参数augf为增广目标函数 文件 名 ，故还需 建立 该文件 ， 详 见附录 augf.m文件 . augf.m中还涉及到（调用）了两个函数： func.m和 constrains.m，前者 为目标函数文件， 与 前面的 func.m一样 ， 只是目标函数不同而已， 后者 为约束函数 ，详 见附录 constrains.m文件 . 现 举 两个例子 ，来说明程序运行方法及运行结果 . 例 1 求解以下约束
18、优化问题 31211221m i n ( ) ( 1 ) ,3. . ( ) 1 0 ,( ) 0 .f x x xs t g x xg x x 建立 func.m文件和 constrains.m文件如下： func.m文件 function f = func( x ) f=(1/3)*(x(1)+1)3+x(2); end constrains.m文件 function g= constrains( x ) g(1)=x(1)-1;%约束 1 g(2)=x(2);%约束 2 end 然后 修改 augf.m文件中 不等式约束的个数变量 inequN=2、去掉外点法罚函数部分的注释并注释掉内点
19、法对应部分（ augf.m文件中有详细说明 ） ;保存后在 命令窗口输入 （粗体部分） ： EPM(-1 -1);或 EPM(-1 -1,10,2);运行结果分别如下图（ a）、（ b）所示： 最优化理论与方法上机报告 9 / 25 (a)采用默认 值的结果 (b)M=10， r=2时的结果 图 2-1-1 例 1 外点法 运行结果 本例 的实际最优解 就是 1,0，可见程序运行正确 ，由上图可知增大 M 和 r的初值可加快收敛速度 . 例 2 求解以下约束优化问题 321212m i n ( ) ,. . ( ) 1 0 .f x x xs t g x x x 建立 func.m文件和 co
20、nstrains.m文件如下： func.m文件 function f = func( x ) f=x(1)3+x(2)2; end constrains.m文件 function g= constrains( x ) g=x(1)+x(2)-1; end 然后修改 augf.m 文件中 不等式约束的个数变量 inequN=1、去掉外点法罚函数部分的注释并注释 掉内点法对应部分（ augf.m 文件中有详细说明 ） ;保存后在命令窗口输入（粗体部分）： EPM(2 -2);或 EPM(2 -2,10,2);运行结果分别如下图（ a）、（ b）所示： 最优化理论与方法上机报告 10 / 25 (
21、a)采用默认 值的结果 (b)M=10， r=2时的结果 图 2-1-2 例 2 外点法 运行结果 本例 没有 实际最优解， 但有局部最优解， 可见程序运行正确，由上图可知增大 M 和 r 的初值可加快收敛速度 . 2. 内点 法 内 点法 一般只适用于 求解如下约束优化 问题： m i n ( ) , ,. . ( ) 0 , 1 , , .nif x x Rs t g x i I m 记 可行域 | ( ) 0 ( ) .n iD x R g x i I 2.1 思想 内点法也属于罚方法的范畴， 其基本思想是保持每一个迭代点 kx 是可行域 D内的点 ，可行域的边界被筑起一道很高的“围墙”
22、作为障碍， 当迭代点靠近边界时，增广目标函数 的值骤然增大，以示 “ 惩罚 ” ，并阻止迭代点穿过边界 .因此内点法也成为内罚函数法或障碍函数法 . 构造高墙 ： 引入 ( , ) ( ) ( ) ,I x M f x M B x其中 ()Bx 满足下列条件 ： 1) ()Bx 在 D 内连续； 2) 在 D 内 当 x 趋于边界点时 ， ()Bx ，当 x 远离边界点时 ( ) 0Bx . 最优化理论与方法上机报告 11 / 25 可取11( ) ,()mi iBx gx或 者21 1( ) ,()mi iBx gx或 者1( ) ln ( ),miiB x g x这样原问题就 转化为求解无
23、约束优化子问题 m i n ( , ) ( ) ( ) .kkI x M f x M B x 2.2 算法 内点法 的算法如下 : Step1 取初始点 0 ,xD 容许误差 01 ， 惩罚因子 1 0 , 0 1,Mc 令k:=1. Step2 以 1kx 为初始点 ， 求解子问题 m i n ( , ) ( ) ( )n kkxR T x M f x M P x ，设其极小点为 kx . Step3 若 ( ) ,kkM B x 则停止迭代，取 * kxx ，否则， 令 1 , 1,kkM cM k k 转 Step2. 2.3 结果 采用 matlab 2011a 实现，无约束优化采用牛
24、顿法，即直接调用第一篇 中的Newton 函数， 具体程序见附录 IPM.m 文件 . 程序中用到的 augf.m文件 、 func.m文件、 constrains.m文件 与外点法用到的为同一文件 ，使用时只需注释掉 另外一部分并取消内点法部分的注释即可 . 仍以外点法中的 两个例子 为例 ，来说明程序运行方法及运行结果 . 只需 修改 augf.m文件中 不等式约束的个数变量 inequN=2、去掉 内 点法罚函数部分的注释并注释掉 外 点法对应部分（ augf.m文件中有详细说明 ） ;保存 对于例 1， 在命令窗口输入（粗体部分）： IPM(2 2);或 IPM(2 2,20,0.5)
25、;运行结果分别如下图（ a）、（ b）所示： 最优化理论与方法上机报告 12 / 25 (a)采用默认 值的结果 (b)M=20， c=0.5时的结果 图 2-2-1 例 1外点法 运行结果 本例的实际最优解就是 1,0，可见程序运行正确，由上图可知增大 M 和 c的初值可加快收敛速度 . 对于例 2，在命令窗口输入（粗体部分）： IPM(1 2);或 IPM(1 2,20,0.3);运行结果分别如下图（ a）、（ b）所示： (a)采用默认 值的结果 (b)M=20， c=0.3时的结果 图 2-2-2 例 2 外点法 运行结果 本例没有实际最优解，但有局部最优解，可见程序运行正确，由上图可
26、知增大 M和 c 的初值可加快收敛速度 . 最优化理论与方法上机报告 13 / 25 附录 源程序清单 1. func.m 文件： function f = func( x ) %-目标函数 - % -刘志 - %-2012/11/12- % 输出：目标函数表达式 % 提示：注释： Ctrl + R 取消注释： Ctrl + T % 提示：注释（取消注释）一行，请将光标置于该行任意位置，使用上述快捷键 % 提示：注释（取消注释）多行，请选中多行 (不需全选 )，使用上述快捷键 % 提示：注意保存！ % -以下为 无约束优化 部分的目标函数 - % -黄金分割法用 - % f=x(1)2-x(1
27、)+6; % -最速下降法、牛顿法用 - % f=(x(1)+10*x(2)2+5*(x(3)-x(4)2+(x(2)-2*x(3)4+10*(x(1)-x(4)4; % -以下为 约束优化 部分的目标函数 - % -对应 约束函数 见 constrains.m 文件 - f=(1/3)*(x(1)+1)3+x(2); % f=2*x(1)+3*x(2); % f=x(1)3+x(2)2; end 2. gfunc.m文件： function gf = gfunc( x ) %目标函数的梯度 % -刘志 - %-2012/11/12- %输出：目标函数的梯度表达式 gf=2*(x(1)+10*
28、x(2)+40*(x(1)-x(4)3 20*(x(1)+10*x(2)+4*(x(2)-2*x(3)3 10*(x(3)-x(4)-8*(x(2)-2*x(3)3 -10*(x(3)-x(4)-40*(x(1)-x(4)3; end 最优化理论与方法上机报告 14 / 25 3. constrains.m 文件： function g= constrains( x ) %-约束函数 - % 不等式约束，也很容易加入等式约束 % 提示：注释： Ctrl + R 取消注释： Ctrl + T % 提示：注释（取消注释）一行，请将光标置于该行任意位置，使用上述快捷键 % 提示：注释（取消注释）多行
29、，请选中多行 (不需全选 )，使用上述快捷键 % 提示：注意保存！ %-对应目标函数 f=(1/3)*(x(1)+1)3+x(2) 的不等式约束 - g(1)=x(1)-1;%约束 1 g(2)=x(2);%约束 2 %-对应目标函数 f=2*x(1)+3*x(2); 的不等式约束 - % g=1-2*x(1)2-x(2)2; %-对应目标函数 f=x(1)3+x(2)2; 的不等式约束 - % g=x(1)+x(2)-1; end 4. augf.m 文件： function af = augf( x ) %-增广目标函数 - % -刘志 - %-2012/12/12- % 此函数用来得到增
30、广目标函数 % 使用 外点法 时，请 注释掉 内点法 部分 % 使用 内点法 时，请注释掉 外点法 部分 % 提示：注释： Ctrl + R 取消注释： Ctrl + T % 提示：注释（取消注释）一行，请将光标置于该行任意位置，使用上述快捷键 % 提示：注释（取消注释）多行，请选中多行 (不需全选 )，使用上述快捷键 % 提示：注意保存！注意修改 inequN 的值！ global M inequN penf; %penalty function %注意修改 inequN 的值！ % M：惩罚 (/障碍 )因子； inequN：不等式 约束个数； penf：罚 (/障碍 )函数 inequN
31、=2; %注意修改 inequN 的值！ penf=0; 最优化理论与方法上机报告 15 / 25 g=constrains(x);%不等式约束 %-以下为外点罚函数 - % global Xk; % G=constrains(Xk);%计算各约束条件在 Xk 处的值 ,以便得到罚函数 % for i=1:inequN % if(G(i)epsilon if objf1objf2 %在 x1,b a=x1; x1=x2; objf1=objf2; x2=a+0.618*(b-a); objf2=objf(x2); else %在 a,x2 b=x2; x2=x1; objf2=objf1; 最
32、优化理论与方法上机报告 17 / 25 x1=a+0.382*(b-a); objf1=objf(x1); end k=k+1; I(k,:)=a,x1,x2,b; end dx=abs(b-a); min=(a+b)/2; fmin=objf(min); if nargout=0 disp( 极小点 :); disp(min); disp( 极小点对应函数值 :); disp(fmin); disp( 迭代次数 :); disp(k); disp( 迭代 区间 :); disp( a x1 x2 b); disp(I); disp( 迭代终止时的区间长度 :); disp(dx); end
33、end 6. fast.m文件： %-最速下降法 (步长采用 0.618法 )- % -刘志 - %-2012/11/12- %- % % function: Min,Fmin,k =golden(objf,gobjf,X0,epsilon) % %输入： objf:目标函数 ; gobjf:目标函数的梯度函数 , % X0：搜索起始点 ; epsilon:精度 (默认为 :10-6). %输出： Min:极小点 , % Fmin:极小点对应函数值矩阵 , % k:迭代次数 . %-% % %使用方法： 最优化理论与方法上机报告 18 / 25 % 1.先建立名为 func的 .m文件，在里面
34、写入表达式，举例如下： % function f = func( x ) % f=100*(x(1)2-x(2)2+(x(1)-1)2; % end % 再建立名为 gfunc的 .m文件，在里面写入 func的梯度表达式 ，如下： % function f = gfunc( x ) % gf=400*x(1)*(x(1)2-x(2)+2*(x(1)-1) % -200*(x(1)2-x(2); % end % 2.在命令窗口中 % 首先输入： X0=1,2 接着 % (1)输入： % fast(func,gfunc,X0， 0.0001); % 程序自动输出所有输出参数，如下： % 极小点
35、: % 1.0001 % 1.0002 % % 极小点对应函数值 : % 1.1837e-008 % % 迭代次数 : % 1452 % (2)输入： % fast(func,gfunc,X0); % 程序按照默认的 epsilon计算，并自动输出所有输出参数 ; % (3) % 如果输入的命令中 ,输出参数不为 0，则按照从左至右的顺序输 % 出 n个参数，但命令结尾不能加 ;号 . % 如输入： min=fast(func,gfunc,X0) 则输出： % Not enough input arguments! % Use default: epsilon=10-6 ! %min = %
36、1.0000 % 1.0000 %function Min,fmin,k = fast( objf,gobjf,X0,epsilon ) %- function Min,fmin,k = fast( objf,gobjf,X0,epsilon ) disp( -最速下降法 (步长采用 0.618法 )-); if nargin=3 epsilon = 0.000001; disp( Not enough input arguments!); 最优化理论与方法上机报告 19 / 25 disp( Use default: epsilon=10-6 !); elseif narginepsilon
37、) step=golden(x)objf(X0+x*d),0,1,epsilon); %d=d/norm(d); X0=X0+step*d; k=k+1; d=-gobjf(X0); end Min=X0; fmin=objf(X0); if norm(d)epsilon disp( 未达到精度 !); end if nargout=0 disp( 极小点 :); disp(Min); disp( 极小点对应函数值 :); disp(fmin); disp( 迭代次数 :); disp(k); end end 7. Newton.m 文件： %-牛顿法 -% % -刘志 - %-2012/11
38、/16-% %-% %输入： objf:目标函数， X0：搜索起始点， % epsilon:精度 (默认为 :10-6); %输出： Min:极小点， fmin极小点对应函数值， % k:迭代次数 %-% function Min,fmin,k = Newton( objf,X0,epsilon ) 最优化理论与方法上机报告 20 / 25 if nargin=2 epsilon = 0.000001; disp( Not enough input arguments! epsilon=10-6); elseif narginepsilon) H=subs(Hess,v,X0); delta=
39、-inv(H)*G; X0=X0+delta; G=subs(gobjf,v,X0); k=k+1; end Min=X0;fmin=objf(X0); if nargout=0 disp( 极小点 :); disp(Min); disp( 极小点对应函数值 :); disp(fmin); disp( 迭代次数 :); disp(k); end end 8. EPM.m 文件： % -外点法 -% % -刘志 - % -2012/12/12-% % -% % function Min,fmin,k = EPM( X0,M1,r,epsilon ) % 输入： X0：搜索起始点， M1：初始惩罚
40、因子（默认为： 1）， 最优化理论与方法上机报告 21 / 25 % r：惩罚因子变化的倍数（默认为： 2）， % epsilon:精度（默认为： 10-6） ; % 输出： Min：极小点， fmin：极小点对应函数值， % k：迭代次数 function Min,fmin,k = EPM( X0,M1,r,epsilon ) global M inequN Xk; % M：惩罚因子； inequN：不等式约束个数 disp( -外点法 -); if nargin=3 epsilon = 0.000001; disp( Not enough input arguments! epsilon=
41、10-6); elseif nargin=1 | cepsilon) Min=Newton(augf,X0,epsilon ); X0=Min; M=M*c; k=k+1; stop=M*penf;%内点法中止条件 :障碍函数准则 end fmin=func(Min); if nargout=0 disp( 极小点 :); disp(Min); disp( 极小点对应函数值 :); disp(fmin); disp( 迭代次数 :); disp(k); end end 所给题为： 321212min. . 1 0xxs t x x 如果使 211xx即 2x 始终比 1x 大 1， 则 121
42、xx， 满足约束条件 ， 那么 1 1 111113 2 3 21 2 1 13211321 1 121 1 1111l i m l i m ( ) l i m ( ( 1 ) )l i m ( ( 1 ) )l i m ( 2 1 )l i m ( ( 2 ) 1 )l i m ( ( ( 1 ) 2 ) 1 )x x xxxxxf x x x xxxx x xx x xxxx 故此函数无最小值。 最优化理论与方法上机报告 24 / 25 举例来说，取 121 8 0 0, 1 8 0 1,xx 则 211 0,xx 满足约束条件，那么 3 2 3 212 ( 1 8 0 0 ) (1 8 0 1 )f x x 很显然上式的值很小，计算器求得： -5828756399， 其实也可以从增长速度上理解，那么 当 1 0x 时， 3 2 3 2 3 21 2 1 1 1 1( 1 ) ( 1 )f x x x x x x 又1 1 1221131 1 11( 1 )( 1 ) 1l i m l i m l i m 0x x xxxx x x ，即分母比分子增长的快 所以 3 2 3 2 3 21 2 1 1 1 1331 1 1( 1 ) ( 1 )( ) ( )f x x x x x xx x x
