最优化实验报告

最优化方法

课程设计报告班级： 学号：===^成绩：2017年5月21日TOC\o"1-5"\h\z一、摘要 1二、单纯形算法 2单纯形算法的基本思路 2\o"CurrentDocument"算法流程图 3用matlab编写源程序 4黄金分割法 7黄金分割法的基本思路 7算法流程图 8用matlab编写源程序 9黄金分割法应用举例 11最速下降法 11最速下降法的基本思路 11算法流程图 13用matlab编写源程序 13最速下降法应用举例 13惩罚函数法 17惩罚函数法的基本思路 17算法流程图 18用matlab编写源程序 18惩罚函数法应用举例 19自我总结 20六、参考文献 20#用matlab编写源程序Matlab程序源代码： simplexTab.m子函数 functionsimplexTab(mat,numFreeVar)maxRow=length(mat(:,1));maxCol=length(mat(1,:));objEntryExcludingMaxPayOff=mat(maxRow,1:maxCol-2);[objEntbestColToPivot]=min(objEntryExcludingMaxPayOff);while(objEnt<0)lastColExcludingObjEnty=mat(1:(maxRow-1),maxCol);ithColExcludingObjEnty=mat(1:(maxRow-1),bestColToPivot);a=lastColExcludingObjEnty./ithColExcludingObjEnty;[valbestRowToPivot]=min(a);sprintf('thebestPivotis%drowand%dcol',bestRowToPivot,bestColToPivot)disp('单纯形表化为：’)；[mat,[a;0]]disp(按任意键继续’)；pause；if(val<0)[sindices]=sort(a)；if(max(a)>0)count=1；while(s(count)<0)count=count+1；endbestRowToPivot=indices(count)；endendif(length(a)==0)length(a)returnendmat=pivot(mat,bestRowToPivot,bestColToPivot)；objEntryExcludingMaxPayOff=mat(maxRow,1:maxCol-2)；[objEntbestColToPivot]=min(objEntryExcludingMaxPayOff)；endsprintf('thebestPicotis%d行and%d歹U',bestRowToPivot,bestColToPivot)disp('单纯形表化为：’)；[mat,[a；0]]disp('运行结束!’)；interChange.m子函数functionnewMat=interChange(mat,row1,row2)temp=mat(row1,:);mat(row1,:)=mat(row2,:);mat(row2,:)=temp;newMat=mat; multiFromRowToRow.m子函数 functionnewMat=multiFromRowToRow(mat,fromRow,toRow,multiplier)rG=mat(fromRow,:)*multiplier;mat(toRow,:)=rG+mat(toRow,:);newMat=mat; multMat.m子函数 fuctionnewMat=multMat(mat,row,mult)%multiplyarowofamatrixbyanonzeroconstantmat(row,:)=mat(row,:)*mult;newMat=mat; pivot.m子函数 functionnewMat=pivot(mat,row,col)%normalizethisrowmat(row,:)=mat(row,:)./mat(row,col);%maketheleadinganumbera1forr=1:length(mat(:,1))if(r~=row)mat=multiFromRowToRow(mat,row,r,-mat(r,col));endendnewMat=mat;单纯形算法应用举例题目：使用单纯形法解下面线性规划问题：目标函数为：maxf=—x—x—x；123约束条件是： /7x+3x+9x<1TOC\o"1-5"\h\z1 2 3st. 8x+5x+4x<11 2 36x+9x+5x<11 2 3Ix,x,x>01 2 3

解：化为标准形式：TOC\o"1-5"\h\zmaxg=x+x+x+0x+0x+0x

123 4 5 67x+3x+9x+x+0x+0x =11 2 34 5 6st. 8x+5x+4x+0x+x+0x=11 2 3 4 5 66x+9x+5x+0x+0x+x=1

1 2 3 4 56x,x,x,x,x,x>0123456利用Matlab程序计算：命令窗口输入：mat=[7391001;8540101;6950011;-1-1-10000];numFreeVar=3;%自由变量个数simplexTab(mat,numFreeVar)结果输出：初始结果：ans=thebestPivotis2rowand1col单纯形表化为：-1.0000 -1.0000 -1.0000 0第一次转换结果ans=thebestPivotis1rowand3col单纯形表化为：00000 -0.3750 -0.5000 0第二次转换结果：ans=thebestPivotis1rowand2col单纯形表化为：0.125000.125000 -0.5000 0 0.0909第三次转换结果ans=thebestPicotis3rowand1col0.045500.13640单纯形表化为0 0 0 0.0593 0.0079 0.0870 0.1542 0运行结束！〔并表示程序运行完毕〕最后单纯形表为xx1x2x3x4x5 (bthetaA000000判别数0000x=(0.0316,0.0870,0.0356),g=0.1542,f=—x—x—x=—0.15421 2 3所以该线性规划的最优解是g=0.1542。二、黄金分割法黄金分割法的基本思路一维搜索是解函数极小值的方法之一,其解法思想为沿某一已知方向求目标函数的极小值点。一维搜索的解法很多,这里主要采用黄金分割法〔0.618法〕。该方法用不变的区间缩短率0.618代替斐波那契法每次不同的缩短率,从而可以看成是斐波那契法的近似,实现起来比较容易,也易于人们所接受。黄金分割法是用于一元函数f(x)在给定初始区间［a,b］内搜索极小点x的min一种方法。它是优化计算中的经典算法,以算法简单、收敛速度均匀、效果较好而著称,是许多优化算法的基础,但它只适用于一维区间上的凸函数,即只在单峰区间内才能进行一维寻优,其收敛效率较低。其基本原理是：依照“去劣存优”原则、对称原则、以及等比收缩原则来逐步缩小搜索区间。具体步骤是：在区间［a,b］内取点：a,a把［a,b］分为三段。12①如果f(a)>f(a)，令a=a,a1=a,a=a+0.618*(b-a)；1 2 1 22

②如果f(a)<f(a),令b=a,a=a,a=b-0.618^(b-a)；1 2 2 2 1 1如果(b-a)/b|和(y1-y2)/y2|都大于收敛精度£重新开始循环。因为［a,b］为单峰区间，这样每次可将搜索区间缩小倍，处理后的区间都将包含极小点的区间缩小，然后在保留下来的区间上作同样的处理，如此迭代下去，将使搜索区［a,b］逐步缩小，直到满足预先给定的精度时，即获得一维优化问题的近似最优解。插入点原理图如下：1—九二九2"九2十九一1=0"1二0.618所谓的黄金分割法，是指将一线段分成两段的方法，使整段的长与较长段的长的长度比值等于较长段与较短段的比值，即：1:九二九:(1—X)算法流程图黄金分割法的算法步骤：给定a,b〔a>b〕，控制误差£>0。一般取［a,b］为［0,1］，当然也可以根据自己需要调整。如果目标函数对变量非常灵敏，可以将补偿的可行范围进行范围缩小，以提高搜索精度；反之，可以将步长的可行范围放大，以提高算法的计算速度。

用matlab编写源程序Matlab程序源代码： golden.m子函数 functionPa=golden(x,p,LowBound,UpBound)%goldenisfunctionforgoldenlinearsearch%x(k+1)=x(k)+Pa*p;%input:% x:搜索初始点% p:搜索方向% LowBound,UpBound:搜索区间%output:% Pa:搜索步长%disp('itisthegoldenlinereasch');% 搜索步长范围［0,1］可以根据需要修改gold_a=LowBound;gold_b=UpBound;z2=gold_a+0.618*(gold_b-gold_a);%%step1f2=objfun(x+z2*p);z1=gold_a+0.382*(gold_b-gold_a);%%step2f1=objfun(x+z1*p);in_itera=0;whileabs(gold_b-gold_a)>0.0000000001%absa=0,b=1goldenlinesearchin_itera=in_itera+1;%f1,f2iff1<f2gold_b=z2;z2=z1;f2=f1;z1=gold_a+0.382*(gold_b-gold_a);f1=objfun(x+z1*p);continueelseiff1==f2gold_a=z1;gold_b=z2;z2=gold_a+0.618*(gold_b-gold_a);%%step2f2=objfun(x+z2*p);z1=gold_a+0.382*(gold_b-gold_a);%%step2f1=objfun(x+z1*p);continueelsegold_a=z1;z1=z2;f1=f2;z2=gold_a+0.618*(gold_b-gold_a);f2=objfun(x+z2*p);endendPa=(gold_a+gold_b)/2;%Pax(k+1)=x(k)+Pa*p;dispC最优解是：')%Pa,in_itera%obj(x+Pa*p)objfun.m子函数functionf=objfun(x)f=x(1)A3+x(2)A2-10*x(1)*x(2)+1;

黄金分割法应用举例题目：利用黄金分割法求解下面函数的最优解：f=x3+x2—10xx+11 2 12解：利用Matlab程序计算：命令窗口输入：X=[0,0];P=[1,1];Pa=golden（X,P,0,1）结果输出：最优解是：三、最速下降法3.1最速下降法的基本思路最速下降法的搜索法向是目标函数的负梯度方向，最速下降法从目标函数的负梯度方向一直前进，直到到达目标函数的最低点。已知目标函数在X（k）点的梯度为：vf（x（k））1±q,..q-'dx dx dx_ 1 2 n-当求目标函数的最小点时，由于函数沿负梯度方向下降最快，故在X（k）点的探索方向应取该点的负梯度方向，即V探索方向应取该点的负梯度方向，即Vf(x(k))显然，s（k）为单位向量。这样第k+1次迭代计算所得的新点为X(k+1)=X(k)+a(k)s(k)=X(k)X(k+1)=X(k)+a(k)s(k)=X(k)—负梯度仅给出了最优化方向，而没有给出步长的大小，所以可能有各种各样的最速下降的过程，它们依赖于a(k速下降的过程，它们依赖于a(k)VfIX(k)的大小。步长a(k)有两种取法：一种方法是任意给定一个初始步长，使满足条件：f(X(k)+a(k)S(k))<f(X(k))另外一种方法是沿负梯度方向做一维探索，以求解一维最优化问题的最优步长a，即对目标函数极小，以得到最优步长：minf(X(k)+aS(k))=f(X(k)+a(k)S(k))a〉0以此最优步长作为由X(k)点出发沿该点的负梯度方向探索的步长a(k)。这种方法的迭代计算的收敛性，可用以下三式中的任一式或二式作为准则来进行判断：卜f(X(k)"1||f(X(k))—f(X(k一1/：|f(X(k)匕|X(k)-X(k-1)||<£用最速下降法求无约束多维极值问题minf(x),xeRn的算法步骤如下：〔1〕、取初始点x(0)，精度e>0，令k=0〔2〕、计算搜索方向v(k)=-Vf(x(k))，其中Vf(x(k))表示函数f(x)在点x(k)处的梯度；〔3〕、假设|V(k)||<£，则停止计算；否则，从x(k)出发，沿v(k)进行一维搜索，即求入，使得f(x(k)+入v(k))=minf(x(k)+入v(k))。此处的一维搜索可以用k k X>0MATLAB的fminunc函数；令x(k+i)=x(k)+九v(k),k=k+1，转步骤〔2〕。k

算法流程图用matlab编写源程序 BanaFunWithGrad.m子函数 function[f,g]=BanaFunWithGrad(x)%(不含导数解析式)f=100*(x(2)-x(1)A2)A2+(1-x(1))A2;g=[100*(4*x(1)A3-4*x(1)*x(2))+2*x(1)-2;100*(2*x(2)-2*x(1)A2)]; BanaFun.m子函数 functionf=BanaFun(x)%(含导数解析式)f=100*(x(2)-x(1)A2)A2+(1-x(1))A2;最速下降法应用举例题目：求解函数f=100(X—X2)2+(1—X)22 1 1

解：利用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)结果输出：迭代次数函数计算次数函数值止工步长一阶导数最优性Gradient'sIterationFunc-countf(x)Step-sizeinfinity-norm01267.621.23e+00312214.4160.000813405519295.836390.0009849313.43155.782920.0005673051.584185.731270.038797912.95245.680230.0005837361.596275.630810.036759912.57335.58190.0006002911.68365.534440.034948212.19425.487430.0006170091.6110455.441730.033324511.711515.396410.0006339251.6212545.352280.031858711.313605.308490.0006510741.6314635.265790.03052711115695.223380.0006684871.6416725.181970.029310610.717785.140820.0006861941.6418815.100590.028193710.40.0007042251.6520905.021440.027163510.121964.982490.0007226091.6622994.944330.02620959.89231054.906350.0007413761.67241084.869120.02532289.65251144.832040.0007605541.68261174.795660.02449589.42271234.75940.0007801741.69281264.723810.0237229.2291324.688330.0008002671.7301354.653470.0229968.99311414.618710.0008208661.7

321444.584530.0223138.79331504.550430.0008420021.71341534.51690.02166898.6351594.483430.0008637111.72361624.450490.02106018.42371684.41760.0008860291.73381714.385220.02048348.240.0009089951.74401804.321030.0199368.07411864.28920.0009326491.75421894.257840.01941547.91431954.22650.0009570331.76441984.195590.01891947.75452044.16470.0009821941.77462074.134230.01844627.6472134.103770.001008181.78482164.073710.01799397.45492224.043640.001035051.79502254.013960.01756097.3512313.984270.001062841.8522343.954950.01714587.16532403.925610.001088861.8542433.894410.01810857.29552493.863210.001117641.82Maximumnumberoffunctionevaluationsexceeded;〔到达最大的函数值计算次数〕〔建议增加最大的函数值计算次数〕x=〔最后迭代点〕fval=〔最优点对应的函数值〕exitflag=0output=(函数基本信息)iterations:56〔迭代次数〕funcCount:250〔目标函数最大计算次数〕stepsize: 〔步长〕algorithm:'medium-scale:Quasi-Newtonlinesearch'message:[1x78char]从结果可见：在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)结果输出：

迭代次数函数计算次数函数值止工步长一阶导数最优性Gradient'sIterationFunc-countf(x)Step-sizeinfinity-norm01267.621.23e+00312214.4160.000813405519295.836390.0009849313.43155.782920.0005673051.584185.731270.038797912.95245.680230.0005837361.596275.630810.036759912.57335.58190.0006002911.68365.534440.034948212.19425.487430.0006170091.6110455.441730.033324511.711515.396410.0006339251.6212545.352280.031858711.313605.308490.0006510741.6314635.265790.03052711115695.223380.0006684871.6416725.181970.029310610.717785.140820.0006861941.6418815.100590.028193710.419875.060580.0007042251.6540014590.0005415110.001117680.016240114620.0005380840.01415250.0574Maximumnumberofiterationsexceeded;increaseoptions.MaxIter.exitflag= 0output=iterations:401funcCount:1462algorithm:'medium-scale:Quasi-Newtonlinesearch'message:[1x67char]从以上两个算法的结果可以看出，最速下降法对bana函数不是十分有效。当目标函数不可微或者导数求解复杂时，可以无需提供目标函数的导函数的解析形式，MATLAB可以使用差分的方法求导〔下降方向〕。命令窗口输入对应的代码为：OPTIONS=optimset('LargeScale','off','HessUpdate','steepdesc','MaxFunEvals',250);x=[-1.9,2];[x,fval,exitflag,output]=fminunc(@BanaFun,x,OPTIONS)结果输出：Maximumnumberoffunctionevaluationsexceeded;x=fval=exitflag= 0output=iterations:19funcCount:252stepsize:algorithm:'medium-scale:Quasi-Newtonlinesearch'message:[1x78char]从计算结果可以明显看出，这对于最速下降发，有导数解析式会比没有导数解析式时，计算质量会有所提高。四、惩罚函数法4.1惩罚函数法的基本思路惩罚函数法是应用广泛,非常有效的间接解法，又称为序列无约束极小化方法(SUMT法)。该方法通过将原约束优化问题中的等式和不等式约束函数加权处理后与原目标函数结合,得到新的目标函数(惩罚函数)。原问题转化为新的无约束优化问题,求解该新的无约束优化问题,间接得到原约束优化问题的最优解。程序步骤：①选择适当的初始罚因子M(0)、初始点X(0)、收敛精度E和罚因子系数c。在本程序中分别取令迭代步数k=0。②采用牛顿法求无约束问题min巾(X,M(k))的极值点X*(M(k))。③检验迭代终止准则，假设满足f[X*(M(k))]-f[X*(M(k-D)]»||X*(M(k))-X(M(k-i))IIV£及|- - —1<£f[X*(M(k-1))]则停止迭代计算，输出最优,“*=X*(M(k))；否则，转入步骤④。④取M(k+1)=cM(k)，X(o)=X*(M(k))，卜=卜+1，转入步骤②继续迭代。

算法流程图用matlab编写源程序Matlab程序源代码： 主程序 symsx1x2M; %M为罚因子。m⑴=1;a⑴=20;b⑴=20;c=8;%c为递增系数。赋初值。f=(x1-2厂2+(x2-1厂2+M*((x「2-x2『2+(x1+x2-2厂2);%外点罚函数f0(1)=500;%求偏导、Hessian元素fx1=diff(f,'x1');fx2=diff(f,'x2');fx1x1=diff(fx1,'x1');fx1x2=diff(fx1,'x2');fx2x1=diff(fx2,'x1');fx2x2=diff(fx2,'x2');%外点法M迭代循环fork=1:100x1=a(k);x2=b(k);M=m(k);%牛顿法求最优值forn=1:100f1=subs(fx1);%求解梯度值和Hessian矩阵f2=subs(fx2);f11=subs(