基于MATLAB的可行方向法求极值问题.doc

机械优化设计报告（5）基于matlab可行方向法求极值的实现姓名：xxx学号：xxx（北京理工大学机械与车辆学院车辆工程，北京 100081）摘要：在工程实际的优化设计中，随着设计变量数和约束条件数的增加，随机方向搜索法和复合形法等直接优化解法的求解效率会偏低。可行方向法，顾名思义，一种始终在可行域内寻找下降方向的搜索法，以其收敛速度快、效果好的优点已成为求解约束非线性问题的一种有代表性的直接解法，同时也是求解大型约束优化问题的主要方法之一。本文将简单介绍可行方向法的数学思想，采用线性规划法和约束最优步长法编写matlab程序，最后通过算例完成对优化问题的求解。关键字：可行方向法；matlab；优化方法。1. 可行方向法的基本数学思想1.1可行方向法的搜索策略可行方向法迭代的第一步都是从可行域的某一初始点出发，沿负梯度方向移至某一个或j个起作用约束面的交集上。以后的搜索路线和迭代计算可根据约束函数和目标函数的不同性状，分别采用以下三种不同策略继续搜索。1） 由点出发，沿可行方向作一维最优化搜索，若所得新点在可行域内，则再沿方向作一维最优化搜索；若所得的新点不在可行域内，则将它移至约束面上再反复重复上述步骤，若，则停止迭代，如图1.1所示。2） 由点出发，沿可行方向作一维最优化搜索，若所得新点在可行域外，则沿可行方向以最大步长到达另一个约束面上一点，将该点作为迭代点进行反复搜索，直至满足给出的k-t条件，如图1.2所示。3） 沿着约束面进行搜索。对于只具有线性约束的非线性规划问题，如图1.3所示，从点出发，沿约束面移动，在有限的几步内即可搜索到约束最优点；对于非线性约束函数，如图1.4所示情况，就是沿着约束面的切线移动，但这样将会进入非可行域，使问题变得复杂。此时，应将进入非可行域的新点x设法调整回约束边界。调整方法为：先规定约束面容差，建立新的约束边界，然后将已经离开约束面的x点沿起作用约束的负梯度方向返回约束面上，计算公式为式中，为调整步长，可以用试探法确定，也可以用下式估计，使其回到约束边界上。 (1-1)图1.1图1.2图1.3图1.4不管采用哪种搜索方法都要进行两条决策：一是产生一个适用的可行方向；二是沿方向确定一个不会越出可行域的适当的步长因子。1.2产生可行方向和步长的数学原理可行方向要保证沿该方向作微小移动后所得的新点是可行点，且目标函数值有所下降。显然可行方向应满足可行和下降两个条件。这里从数学原理的角度分析可行性和下降性。1.2.1可行性设为可行域d中的一个点，即。对于某一方向来说，若存在数，使其对于任意的a，式1-2均成立，则称方向对点满足可行性条件。 （1-2）若点在可行域内，即为内点，则任何方向都满足可行性条件。当点在某一约束边界上时，如图1.5所示，该约束面就是起作用约束。当点位于一个约束面上时，如图1.6所示，作出初起作用约束的梯度和切线，可见只要与起作用约束的梯度成直角或钝角，则可指向可行域内，满足可行性条件。当同时处于j个约束面时，就要求与这j个约束面的梯度均垂直或交成钝角，如图1.6所示，写成数学表达式如下：式中，j为起作用约束的个数。上式称为可行性条件。图1.5图1.61.2.2下降性沿一方向搜索时，要求下降愈快愈好。对某一点来说，负梯度方向为最速下降方向。如果负梯度方向是可行方向，那么负梯度方向是最有利的方向；否则为了保证目标函数值有所下降，至少要使搜索方向与目标函数的负梯度方向成锐角或与梯度方向成钝角，写成数学表达式为该式称为下降性条件。如图1.7所示，可行下降方向显然位于点约束面的切线与目标函数等值线的切线所围成的扇形区域内，推广到一般的情况就是可行方向在目标函数超等值面的超切面和j个起作用约束的超切面所围成的超锥体内。该区域称为可行下降方向区。图1.7综上所述，当点位于j个起作用的约束面上时，满足 （1-3）式中，即为可行方向。1.3可行方向的产生本文采用线性规划法来产生可行方向，用约束最优步长确定最优步长，为后面编制matlab程序做理论基础。本节先阐述线性规划法，下节再对约束最优步长法进行描述。线性规划法对包含线性和非线性的不等式约束的最优化问题都适用，但不允许有等式约束。其基本原理是将具有一阶连续偏导数的目标函数和约束条件在点用taylor展开式展成线性近似函数（一次项），并用这些线性近似函数代替目标函数和它的约束条件，使得问题线性化。这样成为约束条件变为用代替上式中的x，得到s.t. 式中，、为常数，只有搜索方向和是未知量。当迭代点位于j个起作用约束边界的交点上时，问题就变成求解线性规划问题：s.t. 式中，只起到方向作用，故规定其向量的模是有界的，也就是规定每个分量的绝对值不大于1，即由于可行方向还得满足式1-3的可行性和下降性条件，故确定可行方向的线性规划数学模型如下：s.t. (1-4)求解后，可得到迭代方向，作为下一步的搜索方向。1.4最优步长的确定由线性规划法求的可行方向后，由出发，找到该方向上的可行新点。新点必须满足两个条件，一是新点必须还是可行点，二是目标函数值具有最大的下降量。约束最优步长即可满足这两点，可用以下一维搜索法求解：s.t. 式中，。2. 可行方向法的迭代过程和算法流程图2.1迭代过程1)输入收敛精度；2)形成可行初始点，令k=1;3)确定起作用约束下标集合，；4)判断e是否为空集，若是空集则说明在可行域内，转步骤5）；否则说明在约束边界上，转步骤6)；5)判断是否满足终止条件，若满足，转步骤10)；若不满足，则取搜索方向为，转步骤8)；6)采用线性规划法，用式1-4求解后，得到可行方向；7)判断是否成立，若成立，则转步骤10)输出，停机结束；若不成立，则继续下一步；8)采用约束最优步长法进行一维搜索，得到新点;9)令k=k+1，转步骤3)；10)输出，停机结束。2.2算法流程图算法流程图如图2-1。图 2-1 可行方向法的算法流程图3. 可行方向法的matlab程序function kxfxf(x0)eps=1.0e-6;%刚开始给的x0为行向量x0=transpose(x0);funcsz=length(x0);g1 = g;for k=1:1:50 %f(x)梯度： sf=diff_val(x0); sf=eval(sf); %f(x)梯度范数： norm_s=norm(diff_val(x0); gl=length(g); g1 = g; g_copy=g; for i=1:1:gl g=subs(g1(i,:),x,x0); g1(i,:)=g; end g_zero=eval(g1); for i=gl:-1:1 if abs(g_zero(i,:)0.1 g_zero(i,:)=; g_copy(i,:)=; end end i=size(g_zero,1); add=-ones(i,1); %根据e是否为空集(i=0)确定不同情况下的搜索方向 if i=0 if norm_s=3 最优解： 自变量值: x0 目标函数值： f=fval(x0) 迭代次数： k= return else d0=-diff_val(x0); %线性搜索 d0=transpose(d0); vmax = gmin(x0,d0); h=fmin(x0,d0,vmax); x0=x0+h*d0; end else %线性规划 grad=jacobian(g_copy,x); g_zero=subs(grad,x,x0); ac=sf;g_zero; lb=-1*ones(sz,1); ub=ones(sz,1); p=zeros(i+1,1); dz = size(1,sz); c=sf*dz; dz,mn,m1,m2,m3=linprog(c,ac,p,lb,ub); z0 = sf*dz; z0=abs(z0); if z0 x0=3 2 kxfxf(x0)计算结果为：最优解：自变量值:x0 = 1.0e-03 * 0.066121027424648 0.212777459551459ans =目标函数值：f = -11.999735470616054ans =迭代次数：k = 3与书中提供的精确解基本一致，最优解的相对误差为， 说明此方法有效可行。4.2实例结果分析此处主要从以下三个方面来分析此算法及计算结果。1) 可靠性。通过对其他多元多次函数（教科书例5.1）进行求极值计算，结果都比较精确，说明此算法可靠性、通用性好。通过对教课书例5.1进行求解，可得以下结果：x0 = 6.000146685451760 6.999853314348210目标函数值：f = 2.000586785240415迭代次数：k = 3与课本上精确解相比，几乎完全吻合。2) 有效性。对于本例，目标函数不是很复杂，其求解速率在1s内可完成，迭代次数为k=3次。3) 简便性。此程序逻辑清晰，简明易懂，对计算机的要求也不是很高，其简便性不言而明。参考文献1 李志锋。机械优化设计。高等教育出版社。2 孙靖民。机械优化设计。机械工业出版社。1999.5。3 蒲俊，吉家锋，伊良忠。matlab 6.0数学手册。浦东电子出版社。2002.1。附录一 回调函数内容1. func 记录函数及其自变量信息根据所求解的目标函数和相应的约束条件不同而不同，由用户进行手动输入。示例见4.1。2. fval(x0) 计算函数在x0处得函数值function f_val=fval(x0)x0=transpose(x0);func;f_val=subs(f,x,x0);3. gmin(x0,d0) 求函数在满足约束条件下的最小值function h=gmin(x0,d0)syms h;funcn = length(g);alpha = zeros(1,n);a=x0+h*d0;g_val=subs(g,x,a);k = 1;for i=1:1:n alpha_temp = solve(g_val(i); if length(alpha_temp)1 if imag(alpha_temp) = 0 for j=1:1:length(alpha_temp) if alpha_temp(j,1)0 | alpha_temp(j,1) = 0 alpha(1,k) = alpha_temp(j,1); k = k+1; end end else end else if alpha_temp 0 | alpha_temp = 0 alpha(1,k) = alpha_temp; k = k+1; else alpha(:,k) = ; end endendh_temp = alpha(1,1);len = length(alpha);for i=1:1:len if h_temp alpha(1,i) h_temp = alpha(1,i); endendh=h_temp;4. fmin(x0,d0,vmax)求函数在0,vmax上的最小值function h=