最速下降法无约束最优化

《MATLAB程序设计实践》课程考核实践一、编程实现以下科学计算法，并举一例应用之。（参考书籍《精通MATLAB科学计算》，王正林等著，电子工业出版社，2009年）“最速下降法无约束最优化”最速下降法：解：算法说明：最速下降法是一种沿着N维目标函数的负梯度方向搜索最小值的方法。原理：由高等数学知识知道任一点的负梯度方向是函数值在该点下降最快的方向，那么利用负梯度作为极值搜索方向，达到搜寻区间最速下降的目的。而极值点导数性质，知道该点的梯度=0，故而其终止条件也就是梯度逼近于0，也就是当搜寻区间非常逼近极值点时，即：当▽f(a）→0推出f(a）→极值）（xf，f(a）即为所求。该方法是一种局部极值搜寻方法。函数的负梯度表示如下：-g(x）=-▽f(x)=-1)(xxf2)(xxf…TNxxf)(搜索步长可调整，通常记为αk（第k次迭代中的步长）。该算法利用一维的线性搜索方法，如二次逼近法，沿着负梯度方向不断搜索函数的较小值，从而找到最优解。方法特点（1）初始值可任选，每次迭代计算量小，存储量少，程序简短。即使从一个不好的初始点出发，开始的几步迭代，目标函数值下降很快，然后慢慢逼近局部极小点。（2）任意相邻两点的搜索方向是正交的，它的迭代路径胃绕道逼近极小点。当迭代点接近极小点时，步长变得很小，越走越慢。（3）全局收敛，线性收敛，易产生扭摆现象而造成早停。算法步骤：最速下降法的基本求解流程如下：第一步迭代次数初始化为k=0，求出初始点0x的函数值f0=f（0x）。第二步迭代次数加1，即k=k+1，用一维线性搜索方法确定沿负梯度方向-1kg的步长1k，其中1k=ArgMinaf（111k/kkggx）。第三步沿着负梯度方向寻找下一个接近最小值的点，其中步长为1k，得到下一点的坐标为：1111/kkkkkggxx。第四步如果kx≈1kx，且f（kx）≈f（1kx），那么就认为kx为所求的最小值点，并结束循环；否则，就跳到步骤二。流程图：-g(x）=-▽f(x)给定0x，ek=0开始1111/kkkkkggxx。1kx：minf（111k/kkggx）?1exxkk*x=kx结束1k=ArgMinaf（111k/kkggx）k=k+1是否题目：最速下降法求解无约束最优化问题实例。采用最速下降法求如下函数的最小值问题：f（x，y）=x(x-5-y)+y(y-4)即用最速下降法求解函数的最小值问题。解：需先求出该函数的梯度函数。可知其梯度函数为：g(x)=(2x-5-y,-x+2y-4)。源程序代码如下：Opt_Steepest.m文件%用最速下降法求最优化解；function[xo,fo]=Opt_Steepest(f,grad,x0,TolX,TolFun,dist0,MaxIter)%f:函数名；%grad:梯度函数；%x0：搜索初始值；%TolX:最优值点间的误差阈值；%TolFun:函数的误差阈值；%dist0:初始步长；%MaxIter:最大的迭代次数；%xo：最优化点值；%fo:函数在点xo处的函数值。%%%%%%判断输入的变量数，设定一些变量为默认值ifnargin7MaxIter=100;%最大的迭代次数默认为100endifnargin6dist0=10;%初始步长默认为10endifnargin5TolFun=1e-8;%函数值误差为1e-8endifnargin4TolX=1e-6;%自变量距离误差endx=x0;fx0=feval(f,x0);fx=fx0;dist=dist0;kmax1=25;%线性搜索法确定步长的最大搜索次数warning=0;%%%%%迭代计算求最优解fork=1:MaxIterg=feval(grad,x);g=g/norm(g);%求点x处的梯度%%线性搜索方法确定步长dist=dist*2;fx1=feval(f,x-dist*2*g);fork1=1:kmax1fx2=fx1;fx1=feval(f,x-dist*g);iffx0fx1+TolFun&&fx1fx2-TolFun%fx0fx1fx2,den=4*fx1-2*fx0-2*fx2;num=den-fx0+fx2;%二次逼近法dist=dist*num/den;x=x-dist*g;fx=feval(f,x);%确定下一点break;elsedist=dist/2;endendifk1=kmax1warning=warning+1;%无法确定最优步长elsewarning=0;endifwarning=2||(norm(x-x0)TolX&&abs(fx-fx0)TolFun)break;endx0=x;fx0=fx;endxo=x;fo=fx;ifk==MaxIterfprintf('Justbestin%diteration',MaxIter);endQ1.m文件f1004=inline('[x(1)*(x(1)-5-x(2))+x(2)*(x(2)-4)]','x');%目标函数grad=inline('[2*x(1)-5-x(2),-x(1)+2*x(2)-4]','x');%目标函数的梯度函数x0=[14];TolX=1e-4;TolFun=1e-9;MaxIter=100;dist0=1;[xo,fo]=Opt_Steepest(f1004,grad,x0,TolX,TolFun,dist0,MaxIter)运行结果如下：由计算结果可知，当x=4.6667,y=4.3333时，函数f(x,y)=x(x-5-y)+y(y-4)取得最小值-20.3333。二．编程解决以下科学计算和工程实际问题。简支梁受左半均匀分布载荷q及右边L/4处集中力偶M0作用（如下图1-1），求其弯矩、转角和挠度。设L=2m，q=1000N/m，M0=900N*m，E=200*109N/m2,I=2*10-6m4.图1-1①解题思路：首先对简支梁进行受力分析，受力分析图（如下图1-2）所示：图1-2从材料力学的知识可知道，由弯矩求转角要经过一次不定积分，而由转角求挠度又要经过一次不定积分，通常这是很麻烦而且容易出错的，而在MATLAB中，可用cumsum函数或cumtrapz函数作近似的不定积分，只要x取得足够密，其结果将相当准确，而程序非常简单。第一步：计算支反力设支座a和b处的支反力分别为Na和Nb，则据∑Ma=0，∑Fy=0得到平衡方程为：Nb=(q*L^2/8+M0)/LNa=q*L/2-Nb第二步：建立弯矩方程以截面c，d为分界面，将梁划分为ac，cd，db三段分别建立ac，cd，db三段对应的弯矩方程：M1=Na*x-q*x.^2/2;0≦x≦L/2M2=Nb*(L-x)-M0;L/2≦x≦3L/4M3=Nb*(L-x);3L/4≦x≦L第三步：建立挠曲轴近似微分方程并积分建立挠曲轴近似微分方程d2Y/dx2=M(x)/EI对M/EI积分，得转角A，再做一次积分，得挠度Y，每次积分都有一个待定的积分常数。A=∫(M)*dx/(E*I)+Ca=A0（x）+Ca，此处设A0（x）=cumtrapz(M)*dx/(E*I)Y=∫(A)*dx+Cy=∫A0(x)*dx+Ca*x+Cy,此处设Y0（x）=cumtrapz(A0)*dx第四步：确定相应的积分常数Ca，Cy由边界条件Y(0)=0，Y(L)=0确定Y(0)=Y0(0)+Cy=0Y(L)=Y0(L)+Ca*L+Cy=0即[01；L1][Ca；Cy]=[-Y0（0）；-Y0（L）][Ca；Cy]=[0,1;L,1]\[-Y0(0);-Y0（L)];第五步：根据计算结果绘制弯矩、转角以及挠度图形②源程序：L=2;q=1000;M0=900;E=200e9;I=2e-6;%输入已知参数Nb=(q*L^2/8+M0)/L;Na=q*L/2-Nb;%求支反力x=linspace(0,L,101);dx=L/100;%linspace是线性空间向量M1=Na*x(1:51)-q*x(1:51).^2/2;%分三段用数组列出M数组M2=Nb*(L-x(52:76))-M0;M3=Nb*(L-x(77:101));M=[M1,M2,M3];%写出完整的M数组A0=cumtrapz(M)*dx/(E*I);%由M积分求转角AY0=cumtrapz(A0)*dx;%由转角积分求挠度YC=[0,1;L,1]\[-Y0(1);-Y0(101)];%由边界条件求积分常数Ca=C(1),Cy=C(2),A=A0+Ca;Y=Y0+Ca*x+Cy;%A为转角，Y为挠度，求出转角和挠度的完整数值subplot(3,1,1),plot(x,M),grid%绘制弯矩图形subplot(3,1,2),plot(x,A),grid%绘制转角图形subplot(3,1,3),plot(x,Y),grid%绘制挠度图形③运行结果：④流程图：开始输入已知参数L，q，M0，E，I求支反力Nb，Na使x=linspace(0,L,101)，dx=L/100，划分x为空间线性向量分三段用数组列出M数组，写出完整的M数组M=[M1,M2,M3]由M积分求转角A，由转角积分求挠度Y（用cumtrapz积分），由边界条件求积分常数求出转角和挠度的完整数值，A=A0+Ca;Y=Y0+Ca*x+Cy;分别绘制弯矩，转角，挠度的图形用subplot分别绘制弯矩，转角，挠度的图形并输出输出积分常数Ca,Cy结束第三题流程图：开始输入已知的数据表作为样本；设置插值节点针对不同的方法选用相应的函数及格式将已知数据和插值节点代入求得插值节点处的函数值（1）A．正弦值算法：x=0:pi/12:pi/2;y=[00.25880.50000.70710.86600.96591.0000];xi=0:pi/180:pi/2;%三次样条差值yi=interp1(x,y,xi,'spline')%五次多项式拟合A=polyfit(x,y,5);yj=polyval(A,xi)运行结果：yi=Columns1through1100.01750.03490.05240.06980.08720.10450.12190.13920.15640.1737Columns12through220.19080.20790.22490.24190.25880.27560.29230.30900.32550.34200.3583Columns23through330.37460.39070.40670.42260.43840.45400.46950.48480.50000.51500.5299Columns34through440.54460.55920.57360.58780.60180.61570.62930.64280.65610.66910.6820Columns45through550.69470.70710.71930.73130.74310.75470.76600.77710.78800.79860.8090Columns56through660.81910.82900.83870.84800.85710.86600.87460.88290.89100.89870.9062Columns67through770.91350.92040.92710.93350.93960.94540.95100.95630.96120.96590.9703Columns78through880.97440.97820.98170.98490.98780.99040.99270.99460.99630.99770.9987Columns89through910.99950.99991.0000yj=Columns1through110.00000.01740.03490.05230.06970.08710.10450.12180.1391