1、第 0 页 共 23 页最优化方法实验报告Numerical Linear Algebra And Its Applications学生所在学院：理 学 院学生所在班级：计算数学 10-1学 生 姓 名：甘 纯指 导 教 师：单 锐教 务 处2013 年 5 月第 1 页 共 23 页实验三实验名称： 无约束最优化方法的 MATLAB 实现 实验时间: 2013 年 05 月 10 日 星期三 实验成绩： 一、实验目的：通过本次实验的学习，进一步熟悉掌握使用 MATLAB 软件，并能利用该软件进行无约束最优化方法的计算。二、实验背景：（一）最速下降法1、算法原理最速下降法的搜索方向是目标函数的
2、负梯度方向，最速下降法从目标函数的负梯度方向一直前进，直到到达目标函数的最低点。2、算法步骤用最速下降法求无约束问题 的算法步骤如下：nRxf，)(mina）给定初始点 ，精度 ，并令 k=0；)0(x0b）计算搜索方向 ，其中 表示函数 在点)()(kkxfv)(kxf )(xf处的梯度；)(kxc）若 ，则停止计算；否则，从 出发，沿 进行一维搜)(kv )(kx)(kv索，即求 ，使得 ；kmin)( )()(0()kkfvxf 第 2 页 共 23 页d）令 ，转 b） 。1,)()()1( kvxkk（二）牛顿法1、算法原理牛顿法是基于多元函数的泰勒展开而来的，它将作为搜索方向，因此
3、它的迭代公式可直接写出)()(-12kkxff来： )()(12)()( kkk xffx2、算法步骤用牛顿法求无约束问题 的算法步骤如下：nRxf),(mia）给定初始点 ,精度 ，并令 k=0；)0(xb）若 ，停止，极小点为 ，否则转 c） ；()kf )(kxc）计算 ；, )(1)(2)(1)(2 kkffpx令d）令 ，转 b） 。)()()1( kk（三）共轭梯度法1、算法原理共轭梯度法是利用目标函数梯度逐步产生共轭方向作为线搜索方向的方法，每次搜索方向都是在目标函数梯度的共轭方向，搜索步长通过一维极值算法确定。2、算法步骤第 3 页 共 23 页a）给定初始点 ,精度 ；)0(
4、xb）若 ，停止，极小点为 ，否则转 c） ；()f )0(xc） ；)0()0( kxp） ， 且 置取d）用一维搜索方法求 ，使得kt )(min)( ()0() ktktpxfpxf 令 ，转 e） ；)()()1( kkkptxe）若 ，停止，极小值为 ，否则转 f） ；)1(f )1(kxf）若 转 c） ，否则转 g） ；,)()0(nxk令g）令 ，2)(1)1()1( kkkk xffp三、实验内容：1最速下降法的 MATLAB 实现2牛顿法的 MATLAB 实现3共轭梯度法的 MATLAB 实现 四、实验过程：1最速下降法的函数：function x,minf = minFD
5、(f,x0,var,eps)%最速下降法主函数if nargin = 3eps = 1.0e-6;end第 4 页 共 23 页syms l;tol = 1;gradf = - jacobian(f,var);while tolepsv = Funval(gradf,var,x0);tol = norm(v);y = x0 + l*v;yf = Funval(f,var,y);a,b = minJT(yf,0,0.1);%进退法求区间xm = minHJ(yf,a,b);%黄金分割法x1 = x0 + xm*v;x0 = x1;endx = x1;minf = Funval(f,var,x);
6、%进退法函数function minx,maxx = minJT(f,x0,h0,eps)if nargin = 3eps = 1.0e-6;endx1 = x0;k = 0;h = h0;while 1x4 = x1 + h;k = k+1;f4 = subs(f, findsym(f),x4);f1 = subs(f, findsym(f),x1);if f4 eps l = u;u = a + 0.618*(b - a);elseb = u;u = l;第 6 页 共 23 页l = a + 0.382*(b-a);endk = k+1;tol = abs(b - a);endif k
7、= 100000disp(找不到最小值！);x = NaN;minf = NaN;return;endx = (a+b)/2;minf = subs(f, findsym(f),x);2牛顿法的函数：function x,minf = minNT(f,x0,var,eps)if nargin = 3eps = 1.0e-6;endtol = 1;x0 = transpose(x0);gradf = jacobian(f,var);jacf = jacobian(gradf,var);while tolepsv = Funval(gradf,var,x0);tol = norm(v);pv =
8、Funval(jacf,var,x0);p = -inv(pv)*transpose(v);p = double(p);x1 = x0 + p;x0 = x1;endx = x1;minf = Funval(f,var,x);第 7 页 共 23 页3. 共轭梯度法的函数：function x,minf = minGETD(f,x0,var,eps)if nargin = 3eps = 1.0e-6;endx0 = transpose(x0);n = length(var);syms l;gradf = jacobian(f,var);v0 = Funval(gradf,var,x0);p =
9、 -transpose(v0);k = 0;while 1v = Funval(gradf,var,x0);tol = norm(v);if toleps l = u;u = a + 0.618*(b - a);elseb = u;u = l;l = a + 0.382*(b-a);endk = k+1;tol = abs(b - a);endif k = 100000disp(找不到最小值！);x = NaN;minf = NaN;return;end第 10 页 共 23 页x = (a+b)/2;minf = subs(f, findsym(f),x);五、实验结果(总结/方案)1、最速
10、下降法：用最速下降法求函数 的极小值，初始点取1)2()4(, stsf。)3,(0x在 command window 中输入：syms t s;f=(t-4)2+(s+2)2+1;t,mf=minFD(f,1 -3,t s)输出结果：x= 4.0000 -2.0000mf= 12、牛顿法：用牛顿法求函数 的极小值，其中初始点取为1)2()4(, stsf。)0,(x在 command window 中输入：syms t s;f=(t-4)2+(s+2)2+1;t,mf=minNT(f,0 0,t s)输出结果：第 11 页 共 23 页x= 4-2mf= 13、共轭梯度法：用共轭梯度法求函数
11、 的极小值，其中初始值取2)3(,stsf。)6,2(0x在 command window 中输入：syms t s;f=(t-3)2+s2;t,mf=minGETD(f,-2,6,t s)输出结果：x= 3.00000.0000mf= 2.00116e-037第 12 页 共 23 页实验四实验名称： 约束最优化方法的 MATLAB 实现 实验时间: 2013 年 05 月 10 日 星期三 实验成绩： 一、 实验目的：通过本次实验使学生较为熟练使用 MATLAB 软件，并能利用该软件进行约束最优化方法的计算。2、实验内容： 1罚函数法的 MATLAB 实现2可行方向法的 MATLAB 实现
12、三、实验背景：（一）罚函数1、算法原理外罚点函数是通过一系列罚因子 ，求罚函数的极小点来逼ic近原约束问题的最优点。之所以称为外点罚函数法，是因为他是从可行域外部向约束边界逐步靠拢的。2、算法步骤用外点罚函数法求解线性约束优化问题 的算法过程如下：bAxfmin1) 给定初始点 ，罚参数列 及精度 ，置 k=1；0xic0第 13 页 共 23 页2）构造罚函数 ；2bAxcfxF3）用某种无约束非线性规划，以 为初始点求解 ；1k xFmin4）设最优解为 ，若 满足某种终止条件，则停止迭代输出 ，kxk k否则令 ，转 2） 。1罚函数列 的选法：通常先选定一个初始常数 和一个比例系数ic
13、 1c，则其余的可表示为 。终止条件可采用 ，其中21iicxS。2bAxcS（二）可行方向法1、算法原理可行方向法是求解如下约束最优化问题 的算法，其中bAxfmin约束条件为线性约束，A 为约束系数矩阵，b 为约束向量。其基本思路是从可行点出发，沿着目标函数值减小的方向搜索求出新的可行点，如此迭代下去。可行点是满足约束条件的点。2、算法步骤 可行方向法的算法过程如下：1）给定初始可行点 ，使其满足约束条件 ，置 ；1xbAx1k2）在 处，将 A，b 分解成 和 ，使 ， ；kx21A21b11k22bx3）如果 是空的，则令 p=1，否则令 ;1 11AIPT第 14 页 共 23 页4
14、)计算 ，若 则转 6） ，否则转 5） ；kkxfPd0kd5）若 是空的，则停止计算，输出 ，否则计算1Akx；如果 ，则停止计算，输出 ，若 包kTxfw1wkxw含负的分量，则选择一个负分量，去掉 对应的行，转 3） ；1A6）求一维约束问题 ，其中 的计算方max0.,mintsdxfkk max法如下： ， 求出最优解，2 bAdk其 他,inaxiib设为 ，令 ，置 ，转 2） 。kkkkdx1 1k四、实验过程：1、罚函数法的函数function x,minf = minPF(f,x0,A,b,c1,p,var,eps)format long;if nargin = 7eps
15、 = 1.0e-4;endk = 0;FE = 0;for i=1:length(b)FE = FE + (var*transpose(A(1,:) - b(i)2;endx1 = transpose(x0);x2 = inf;while 1M = c1*p;FF = M*FE;SumF = f + FF;x2,minf = minNT(SumF,transpose(x1),var);%牛顿法函数第 15 页 共 23 页if norm(x2 - x1)epsv = Funval(gradf,var,x0);tol = norm(v);pv = Funval(jacf,var,x0);p =
16、-inv(pv)*transpose(v);p = double(p);x1 = x0 + p;x0 = x1;endx = x1;minf = Funval(f,var,x);2、可行方向法的函数function x,minf = minRosen(f,A,b,x0,var,eps)第 16 页 共 23 页if nargin = 5eps = 1.0e-6;endsyms l;x0 = transpose(x0);n = length(var);sz = size(A);m = sz(1);gf = jacobian(f,var);bConti = 1;while bContik = 0;
17、s = 0;A1 = A;A2 = A;b1 = b;b2 = b;for i=1:mdfun = A(i,:)*x0 - b(i);if abs(dfun) 0A1 = A1(1:k,:);b1 = b1(1:k,:);endif s 0A2 = A2(1:s,:);b2 = b2(1:s,:);第 17 页 共 23 页endwhile 1P = eye(n,n);if k 0tM = transpose(A1);P = P - tM*inv(A1*tM)*A1;endgv = Funval(gf, var, x0);gv = transpose(gv);d = -P*gv;if d =
18、0if k = 0x = x0;bConti = 0;break;elsew = inv(A1*tM)*A1*gv;if w=0x = x0;bConti = 0;break;elseu,index = min(w);sA1 = size(A1);if sA1(1) = 1k = 0;elsek = sA1(2);A1 = A1(1:(index-1),:); A1(index+1):sA1(2),:);endendendelsebreak;endendyl = x0 + l*d;第 18 页 共 23 页tmpf = Funval(f,var,yl);bb = b2 - A2*x0;dd =
19、 A2*d;if dd = 0tmpI,lm = minJT(tmpf,0,0.1);%进退法函数elselm = inf;for i=1:length(dd)if dd(i) eps l = u;u = a + 0.618*(b - a);elseb = u;u = l;l = a + 0.382*(b-a);endk = k+1;tol = abs(b - a);endif k = 100000disp(找不到最小值！);x = NaN;minf = NaN;return;endx = (a+b)/2;minf = subs(f, findsym(f),x);五、实验结果(总结/方案)1、
20、罚函数法：用外点罚函数求下面的优化问题： 其中取1.25.0.),(minststf，初始点取 。2,05.1c )0,(x在 command window 中输入：f=0.5*t2+s2/4;A=1,1;b=1;c1=0.05;p=2;第 21 页 共 23 页x.minf=minPF(f,0 0,A,b,c1,p,t,s) 输出结果：x= 0.3333 0.6666mf= 0.16662、可行方向法：用 Rosen 梯度投影法求下面函数的极小值：，其中初始可行点取8534. 2832),(minsts stttf )0,(x在 command window 中输入：syms t s;f=2*t2+s2-2*t*s+3*t-8*s+2;x,mf=minRosen(f,-1 1;-3 -5,-4;-8,0,0,t,s)输出结果：x= -0.37641.8258mf= -8.7444第 22 页 共 23 页
