最优化实验报告.doc

最优化方法课程设计报告 班级：_ 姓名： _ 学号： _ 成绩： 2017年 5月 21 日 目 录一、摘 要1二、单纯形算法21.1 单纯形算法的基本思路21.2 算法流程图31.3 用matlab编写源程序3二、 黄金分割法72.1 黄金分割法的基本思路72.2 算法流程图82.3 用matlab编写源程序92.4 黄金分割法应用举例10三、 最速下降法103.1 最速下降法的基本思路103.2 算法流程图123.3 用matlab编写源程序123.4 最速下降法应用举例13四、 惩罚函数法164.1 惩罚函数法的基本思路164.2 算法流程图174.3 用matlab编写源程序174.4 惩罚函数法应用举例19五、 自我总结19六、参考文献19一、摘 要运筹学是一门以人机系统的组织、管理为对象，应用数学和计算机等工具来研究各类有限资源的合理规划使用并提供优化决策方案的科学。通过对数据的调查、收集和统计分析，以及具体模型的建立。收集和统计上述拟定之模型所需要的各种基础数据，并最终将数据整理形成分析和解决问题的具体模型。最优化理论和方法日益受到重视，已经渗透到生产、管理、商业、军事、决策等各个领域，而最优化模型与方法广泛应用于工业、农业、交通运输、商业、国防、建筑、通信、政府机关等各个部门及各个领域。伴随着计算机技术的高速发展，最优化理论与方法的迅速进步为解决实际最优化问题的软件也在飞速发展。其中，MATLAB软件已经成为最优化领域应用最广的软件之一。有了MATLAB这个强大的计算平台，既可以利用MATLAB优化工具箱（OptimizationToolbox）中的函数，又可以通过算法变成实现相应的最优化计算。关键词：优化、线性规划、黄金分割法、最速下降法、惩罚函数法 二、单纯形算法1.1 单纯形算法的基本思路线性规划问题的可行域是n维向量空间Rn中的多面凸集，其最优值如果存在必在该凸集的某顶点处达到。顶点所对应的可行解称为基本可行解。单纯形法的基本思想是：先找出一个基本可行解，对它进行鉴别，看是否是最优解；若不是，则按照一定法则转换到另一改进的基本可行解,再鉴别；若仍不是，则再转换,按此重复进行。因基本可行解的个数有限，故经有限次转换必能得出问题的最优解。如果问题无最优解也可用此法判别。单纯形法是从某一基可行解出发，连续地寻找相邻的基可行解，直到达到最优的迭代过程，其实质是解线性方程组。概述：根据单纯形法的原理，在线性规划问题中，决策变量（控制变量）x1，x2，xn的值称为一个解,满足所有的约束条件的解称为可行解。使目标函数达到最大值（或最小值）的可行解称为最优解。这样，一个最优解能在整个由约束条件所确定的可行区域内使目标函数达到最大值（或最小值）。求解线性规划问题的目的就是要找出最优解。用单纯形法求解线性规划问题所需的迭代次数主要取决于约束条件的个数。现在一般的线性规划问题都是应用单纯形法标准软件在计算机上求解，对于具有106个决策变量和104个约束条件的线性规划问题已能在计算机上解得。求解时可能出现下列情况之一：存在着一个最优解；存在着无穷多个最优解；不存在最优解，这只在两种情况下发生，即没有可行解或各项约束条件不阻止目标函数的值无限增大（或向负的方向无限增大）。要缩小对最优解的搜索范围，就必须认识最优解的一般性质，最优解如果存在的话，则它必然处于可行区域的边界上。任何一项约束条件的边界方程是用“”号来替换该约束条件中的“”或“”号而得到的。每一个边界方程确定一个超平面。因此，可行区域的边界是由那些满足一个或同时满足几个边界方程（即处在作为边界的一个或几个超平面上）的可行解所组成，而且最优解必在其中。最优解不仅是在可行区域的边界上，而且也在这个区域的一个隅角上。一个可行解，如果不处在由另两个可行解连接起来的任何线段上，它就是一个角点可行解。如果连接两个角点可行解的线段处在可行区域的边界上，这两个角点可行解就称为相邻的角点可行解。角点可行解具有下列三个重要性质：如果存在着一个最优解，那么它必定是角点可行解。如果存在有多个最优解，那么至少有两个最优解必定是相邻的角点可行解。只存在有限个数的角点可行解。如果一个角点可行解按目标函数值来衡量时比其所有的相邻角点可行解更好一些，那它就比所有其他角点可行解都更好，也就是最优解。上述这些性质构成单纯形法的原理基础。最后一个性质的重要性在于它为一个角点可行解是否是最优解提供了一种简便的检验标准，因而毋需列举所有的可行解。单纯形法正是利用了这个性质，只要检查少数的角点可行解，并且一旦这个最优性检验获得通过就可立即停止运算。1.2 算法流程图（1）、确定初始基可行解从线性规划标准形的系数矩阵中能直接找出m个线性独立的单位向量；对约束条件全为“=”连接的LP，化为标准形，左端添加松弛变量后即形成一个单位子矩阵；约束条件中含有“=”或“=”连接的方程，在插入剩余变量后找不到单位矩阵，则必须采用“人造基”法，（2）、单纯形法的运算步骤可归结为：起始步骤在一个角点可行解上开始。迭代步骤移动至一个更好一些的相邻角点可行解（根据需要反复进行这一步骤）。停止法则在当前角点可行解比所有相邻角点可行解都更好些时停止。当前角点可行解就是一个最优解。单纯形法的优点及其成功之处在于它只需要较少的有限次数的迭代，即可找到最优解。（3）、单纯形法的算法流程如下：把线性规划问题的约束方程组表达成典范型方程组，找出基本可行解作为初始基本可行解。若基本可行解不存在，即约束条件有矛盾，则问题无解。若基本可行解存在，从初始基本可行解作为起点，根据最优性条件和可行性条件，引入非基变量取代某一基变量，找出目标函数值更优的另一基本可行解。按步骤3进行迭代,直到对应检验数满足最优性条件（这时目标函数值不能再改善），即得到问题的最优解。若迭代过程中发现问题的目标函数值无界，则终止迭代。化标准形式、松弛变量求初始基本可行解选择非基本变量Y输出，终止程序转轴对换选择基本变量无解，终止程序NNY1.3 用matlab编写源程序Matlab程序源代码：- simplexTab.m子函数-function simplexTab(mat,numFreeVar)maxRow=length(mat(:,1);maxCol=length(mat(1,:);objEntryExcludingMaxPayOff=mat(maxRow,1:maxCol-2);objEnt bestColToPivot=min(objEntryExcludingMaxPayOff);while(objEnt0) lastColExcludingObjEnty=mat(1:(maxRow-1),maxCol); ithColExcludingObjEnty=mat(1:(maxRow-1),bestColToPivot); a=lastColExcludingObjEnty./ithColExcludingObjEnty; val bestRowToPivot=min(a); sprintf(the best Pivot is %d row and %d col ,bestRowToPivot,bestColToPivot) disp(单纯形表化为：); mat,a;0 disp(按任意键继续); pause; if(val0) count=1; while(s(count)，令a=, a1=, =a+0.618*(b-a)； 如果b），控制误差。一般取a, b为0,1，当然也可以根据自己需要调整。如果目标函数对变量非常灵敏，可以将补偿的可行范围进行范围缩小，以提高搜索精度；反之，可以将步长的可行范围放大，以提高算法的计算速度。2.3 用matlab编写源程序Matlab程序源代码：-golden.m子函数-function Pa=golden(x,p,LowBound,UpBound)%golden is function for golden linear search%x (k+1) =x (k) +Pa*p;%input:% x:搜索初始点% p:搜索方向% LowBound,UpBound:搜索区间%output:% Pa:搜索步长%disp (it is the golden linereasch);% 搜索步长范围0,1可以根据需要修改 gold_a=LowBound; gold_b=UpBound; z2=gold_a+0.618*(gold_b-gold_a); %step1 f2=objfun(x+z2*p); z1=gold_a+0.382*(gold_b-gold_a); %step2 f1=objfun(x+z1*p); in_itera=0; while abs(gold_b-gold_a)0.0000000001 % abs a=0,b=1 golden linesearch in_itera=in_itera+1; %f1, f2 if f1f2 gold_b=z2; z2=z1; f2=f1; z1=gold_a+0.382*(gold_b-gold_a); f1=objfun(x+z1*p); continue elseif f1=f2 gold_a=z1; gold_b=z2; z2=gold_a+0.618*(gold_b-gold_a); %step2 f2=objfun(x+z2*p); z1=gold_a+0.382*(gold_b-gold_a); %step2 f1=objfun(x+z1*p); continue else gold_a=z1; z1=z2; f1=f2; z2=gold_a+0.618*(gold_b-gold_a); f2=objfun(x+z2*p); end endPa= (gold_a+gold_b)/2; %Pa x(k+1)=x(k)+Pa*p;disp(最优解是：)%Pa,in_itera%obj(x+Pa*p)-objfun.m子函数-function f=objfun(x)f=x(1)3+x(2)2-10*x(1)*x(2)+1;-2.4 黄金分割法应用举例题目：利用黄金分割法求解下面函数的最优解：解： 利用Matlab程序计算：命令窗口输入：X=0,0;P=1,1;Pa=golden(X,P,0,1)结果输出：最优解是： Pa = 1.0000三、 最速下降法3.1 最速下降法的基本思路最速下降法的搜索法向是目标函数的负梯度方向，最速下降法从目标函数的负梯度方向一直前进，直到到达目标函数的最低点。已知目标函数在点的梯度为：当求目标函数的最小点时，由于函数沿负梯度方向下降最快，故在点的探索方向应取该点的负梯度方向，即显然，为单位向量。这样第次迭代计算所得的新点为负梯度仅给出了最优化方向，而没有给出步长的大小，所以可能有各种各样的最速下降的过程，它们依赖于的大小。步长有两种取法：一种方法是任意给定一个初始步长，使满足条件：另外一种方法是沿负梯度方向做一维探索，以求解一维最优化问题的最优步长，即对目标函数极小，以得到最优步长：以此最优步长作为由点出发沿该点的负梯度方向探索的步长。这种方法的迭代计算的收敛性，可用以下三式中的任一式或二式作为准则来进行判断：用最速下降法求无约束多维极值问题的算法步骤如下：（1）、取初始点，精度，令（2）、计算搜索方向，其中表示函数在点处的梯度；（3）、若，则停止计算；否则，从出发，沿进行一维搜索，即求，使得。此处的一维搜索可以用MATLAB的函数；令，转步骤（2）。3.2 算法流程图3.3 用matlab编写源程序-BanaFunWithGrad.m子函数-function f,g=BanaFunWithGrad(x) %(不含导数解析式)f=100*(x(2)-x(1)2)2+(1-x(1)2;g=100*(4*x(1)3-4*x(1)*x(2)+2*x(1)-2;100*(2*x(2)-2*x(1)2);-BanaFun.m子函数-function f=BanaFun(x) %(含导数解析式)f=100*(x(2)-x(1)2)2+(1-x(1)2;-3.4 最速下降法应用举例题目：求解函数解： 利用Matlab程序计算：命令窗口输入：OPTIONS=optimset(LargeScale,off,HessUpdate,steepdesc,gradobj,on,MaxFunEvals,250,display,iter);x=-1.9,2;x,fval,exitflag,output=fminunc(BanaFunWithGrad,x,OPTIONS)结果输出： 迭代次数 函数计算次数 函数值 步长 一阶导数最优性 Gradients Iteration Func-count f(x) Step-size infinity-norm 0 1 267.62 1.23e+003 1 2 214.416 0.000813405 519 2 9 5.83639 0.00098493 13.4 3 15 5.78292 0.000567305 1.58 4 18 5.73127 0.0387979 12.9 5 24 5.68023 0.000583736 1.59 6 27 5.63081 0.0367599 12.5 7 33 5.5819 0.000600291 1.6 8 36 5.53444 0.0349482 12.1 9 42 5.48743 0.000617009 1.61 10 45 5.44173 0.0333245 11.7 11 51 5.39641 0.000633925 1.62 12 54 5.35228 0.0318587 11.3 13 60 5.30849 0.000651074 1.63 14 63 5.26579 0.0305271 11 15 69 5.22338 0.000668487 1.64 16 72 5.18197 0.0293106 10.7 17 78 5.14082 0.000686194 1.64 18 81 5.10059 0.0281937 10.4 19 87 5.06058 0.000704225 1.65 20 90 5.02144 0.0271635 10.1 21 96 4.98249 0.000722609 1.66 22 99 4.94433 0.0262095 9.89 23 105 4.90635 0.000741376 1.67 24 108 4.86912 0.0253228 9.65 25 114 4.83204 0.000760554 1.68 26 117 4.79566 0.0244958 9.42 27 123 4.7594 0.000780174 1.69 28 126 4.72381 0.023722 9.2 29 132 4.68833 0.000800267 1.7 30 135 4.65347 0.022996 8.99 31 141 4.61871 0.000820866 1.7 32 144 4.58453 0.022313 8.79 33 150 4.55043 0.000842002 1.71 34 153 4.5169 0.0216689 8.6 35 159 4.48343 0.000863711 1.72 36 162 4.45049 0.0210601 8.42 37 168 4.4176 0.000886029 1.73 38 171 4.38522 0.0204834 8.24 39 177 4.35289 0.000908995 1.74 40 180 4.32103 0.019936 8.07 41 186 4.2892 0.000932649 1.75 42 189 4.25784 0.0194154 7.91 43 195 4.2265 0.000957033 1.76 44 198 4.19559 0.0189194 7.75 45 204 4.1647 0.000982194 1.77 46 207 4.13423 0.0184462 7.6 47 213 4.10377 0.00100818 1.78 48 216 4.07371 0.0179939 7.45 49 222 4.04364 0.00103505 1.79 50 225 4.01396 0.0175609 7.3 51 231 3.98427 0.00106284 1.8 52 234 3.95495 0.0171458 7.16 53 240 3.92561 0.00108886 1.8 54 243 3.89441 0.0181085 7.29 55 249 3.86321 0.00111764 1.82 Maximum number of function evaluations exceeded;（达到最大的函数值计算次数）increase options.MaxFunEvals（建议增加最大的函数值计算次数）x = -0.9634 0.9372（最后迭代点）fval = 3.8632（最优点对应的函数值）exitflag = 0output = (函数基本信息) iterations: 56 （迭代次数） funcCount: 250（目标函数最大计算次数） stepsize: 0.0011（步长） firstorderopt: 1.8155 algorithm: medium-scale: Quasi-Newton line search message: 1x78 char从结果可见：在250次最大的目标函数计算次数内，算法没有到最优点1,1。将最大目标函数计算次数调高到2500次，则：命令窗口输入：OPTIONS=optimset(LargeScale,off,HessUpdate,steepdesc,gradobj,on,MaxFunEvals,2500,display,iter);x=-1.9,2;x,fval,exitflag,output=fminunc(BanaFunWithGrad,x,OPTIONS)结果输出： 迭代次数 函数计算次数 函数值 步长 一阶导数最优性 Gradients Iteration Func-count f(x) Step-size infinity-norm 0 1 267.62 1.23e+003 1 2 214.416 0.000813405 519 2 9 5.83639 0.00098493 13.4 3 15 5.78292 0.000567305 1.58 4 18 5.73127 0.0387979 12.9 5 24 5.68023 0.000583736 1.59 6 27 5.63081 0.0367599 12.5 7 33 5.5819 0.000600291 1.6 8 36 5.53444 0.0349482 12.1 9 42 5.48743 0.000617009 1.61 10 45 5.44173 0.0333245 11.7 11 51 5.39641 0.000633925 1.62 12 54 5.35228 0.0318587 11.3 13 60 5.30849 0.000651074 1.63 14 63 5.26579 0.0305271 11 15 69 5.22338 0.000668487 1.64 16 72 5.18197 0.0293106 10.7 17 78 5.14082 0.000686194 1.64 18 81 5.10059 0.0281937 10.4 19 87 5.06058 0.000704225 1.65 400 1459 0.000541511 0.00111768 0.0162 401 1462 0.000538084 0.0141525 0.0574 Maximum number of iterations exceeded; increase options.MaxIter.x = 0.9770 0.9542fval = 5.3808e-004exitflag = 0output = iterations: 401 funcCount: 1462 stepsize: 0.0142 firstorderopt: 0.0574 algorithm: medium-scale: Quasi-Newton line search message: 1x67 char从以上两个算法的结果可以看出，最速下降法对bana函数不是十分有效。当目标函数不可微或者导数求解复杂时，可以无需提供目标函数的导函数的解析形式，MATLAB可以使用差分的方法求导（下降方向）。命令窗口输入对应的代码为：OPTIONS=optimset(LargeScale,off,HessUpdate,steepdesc,MaxFunEvals,250);x=-1.9,2;x,fval,exitflag,output=fminunc(BanaFun,x,OPTIONS)结果输出：Ma