最速下降法无约束最优化.doc
- 文档编号:2505063
- 上传时间:2023-05-03
- 格式:DOC
- 页数:17
- 大小:440KB
最速下降法无约束最优化.doc
《最速下降法无约束最优化.doc》由会员分享,可在线阅读,更多相关《最速下降法无约束最优化.doc(17页珍藏版)》请在冰点文库上搜索。
《MATLAB程序设计实践》课程考核
实践一、编程实现以下科学计算法,并举一例应用之。
(参考书籍《精通MATLAB科学计算》,王正林等著,电子工业出版社,2009年)
“最速下降法无约束最优化”
最速下降法:
解:
算法说明:
最速下降法是一种沿着N维目标函数的负梯度方向搜索最小值的方法。
原理:
由高等数学知识知道任一点的负梯度方向是函数值在该点下降最快的方向,那么利用负梯度作为极值搜索方向,达到搜寻区间最速下降的目的。
而极值点导数性质,知道该点的梯度=0,故而其终止条件也就是梯度逼近于0,也就是当搜寻区间非常逼近极值点时,即:
当▽f(a)→0推出f(a)→,f(a)即为所求。
该方法是一种局部极值搜寻方法。
函数的负梯度表示如下:
-g(x)=-▽f(x)=-…
搜索步长可调整,通常记为αk(第k次迭代中的步长)。
该算法利用一维的线性搜索方法,如二次逼近法,沿着负梯度方向不断搜索函数的较小值,从而找到最优解。
方法特点
(1)初始值可任选,每次迭代计算量小,存储量少,程序简短。
即使从一个不好的初始点出发,开始的几步迭代,目标函数值下降很快,然后慢慢逼近局部极小点。
(2)任意相邻两点的搜索方向是正交的,它的迭代路径胃绕道逼近极小点。
当迭代点接近极小点时,步长变得很小,越走越慢。
(3)全局收敛,线性收敛,易产生扭摆现象而造成早停。
算法步骤:
最速下降法的基本求解流程如下:
第一步
迭代次数初始化为k=0,求出初始点的函数值f=f()。
第二步
迭代次数加1,即k=k+1,用一维线性搜索方法确定沿负梯度方向-的步长,其中=ArgMinaf()。
第三步
沿着负梯度方向寻找下一个接近最小值的点,其中步长为,得到下一点的坐标为:
。
第四步
如果≈,且f()≈f(),那么就认为为所求的最小值点,并结束循环;否则,就跳到步骤二。
流程图:
-g(x)=-▽f(x)
给定,e
k=0
开始
。
:
minf()
=
结束
=ArgMinaf()
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处的函数值。
%%%%%%判断输入的变量数,设定一些变量为默认值
ifnargin<7
MaxIter=100;%最大的迭代次数默认为100
end
ifnargin<6
dist0=10;%初始步长默认为10
end
ifnargin<5
TolFun=1e-8;%函数值误差为1e-8
end
ifnargin<4
TolX=1e-6;%自变量距离误差
end
x=x0;
fx0=feval(f,x0);
fx=fx0;
dist=dist0;
kmax1=25;%线性搜索法确定步长的最大搜索次数
warning=0;
%%%%%迭代计算求最优解
fork=1:
MaxIter
g=feval(grad,x);
g=g/norm(g);%求点x处的梯度
%%线性搜索方法确定步长
dist=dist*2;
fx1=feval(f,x-dist*2*g);
fork1=1:
kmax1
fx2=fx1;
fx1=feval(f,x-dist*g);
iffx0>fx1+TolFun&&fx1 den=4*fx1-2*fx0-2*fx2;num=den-fx0+fx2;%二次逼近法 dist=dist*num/den; x=x-dist*g;fx=feval(f,x);%确定下一点 break; else dist=dist/2; end end ifk1>=kmax1 warning=warning+1;%无法确定最优步长 else warning=0; end ifwarning>=2||(norm(x-x0) break; end x0=x; fx0=fx; end xo=x;fo=fx; ifk==MaxIter fprintf('Justbestin%diteration',MaxIter); end Q1.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)/L Na=q*L/2-Nb 第二步: 建立弯矩方程 以截面c,d为分界面,将梁划分为ac,cd,db三段 分别建立ac,cd,db三段对应的弯矩方程: M1=Na*x-q*x.^2/2;0≦x≦L/2 M2=Nb*(L-x)-M0;L/2≦x≦3L/4 M3=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=0 Y(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积分求转角A Y0=cumtrapz(A0)*dx;%由转角积分求挠度Y C=[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= Columns1through11 00.01750.03490.05240.06980.08720.10450.12190.13920.15640.1737 Columns12through22 0.19080.20790.22490.24190.25880.27560.29230.30900.32550.34200.3583 Columns23through33 0.37460.39070.40670.42260.43840.45400.46950.48480.50000.51500.5299 Columns34through44 0.54460.55920.57360.58780.60180.61570.62930.64280.65610.66910.6820 Columns45through55 0.69470.70710.71930.73130.74310.75470.76600.77710.78800.79860.8090 Columns56through66 0.81910.82900.83870.84800.85710.86600.87460.88290.89100.89870.9062 Columns67through77 0.91350.92040.92710.93350.93960.94540.95100.95630.96120.96590.9703 Columns78through88 0.97440.97820.98170.98490.98780.99040.99270.99460.99630.99770.9987 Columns89through91 0.99950.99991.0000 yj= Columns1through11 0.00000.01740.03490.05230.06970.08710.10450.12180.13910.15640.1736 Columns12through22 0.19080.20790.22490.24190.25880.27560.29240.30900.32560.34200.3584 Columns23through33 0.37460.39070.40670.42260.43840.45400.46950.48480.50000.51500.5299 Columns34through44 0.54460.55920.57360.58780.60180.61570.62930.64280.65610.66910.6820 Columns45through55 0.69460.70710.71930.73130.74310.75470.76600.77710.78800.79860.8090 Columns56through66 0.81910.82900.83860.84800.85710.86600.87460.88290.89100.89880.9063 Columns67through77 0.91350.92050.92720.93360.93970.94550.95100.95630.96120.96590.9703 Columns78through88 0.97430.97810.98160.98480.98770.99020.99250.99450.99620.99750.9986 Columns89through91 0.99940.99981.0000 通过比较,两种方法得到的结果近似相等。 B.正切值算法: x=0: pi/12: 5*pi/12; y=[00.26790.57741.00001.73203.7320]; xi=0: pi/180: 5*pi/12;%三次样条差值 yi=interp1(x,y,xi,'spline') %五次多项式拟合 A=polyfit(x,y,5); yj=polyval(A,xi) 运行结果: yi= Columns1through11 00.01840.03650.05450.07240.09020.10790.12550.14310.16070.1784 Columns12through22 0.19610.21380.23170.24970.26790.28630.30480.32360.34270.36200.3817 Columns23through33 0.40170.42210.44290.46410.48580.50790.53050.55370.57740.60170.6266 Columns34through44 0.65200.67800.70460.73170.75930.78760.81630.84560.87540.90580.9367 Columns45through55 0.96811.00001.03251.06581.10031.13641.17431.21451.25721.30281.3516 Columns56through66 1.40411.46041.52111.58631.65651.73201.81311.90021.99362.09372.2008 Columns67through76 2.31522.43742.56752.70602.85323.00953.17523.35063.53613.7320 yj= Columns1through11 -0.00000.02350.04540.06580.08500.10320.12060.13750.15400.17010.1862 Columns12through22 0.20220.21830.23450.25110.26790.28510.30280.32080.33940.35850.3781 Columns23through33 0.39820.41880.44000.46160.48380.50650.52970.55330.57740.60200.6270 Columns34through44 0.65240.67830.70470.73150.75880.78670.81500.84400.87360.90390.9351 Columns45through55 0.96701.00001.03411.06931.10601.14421.18411.22591.26991.31621.3652 Columns56through66 1.41711.47231.53101.59351.66041.73201.80871.89101.97932.07422.1762 Columns67through76 2.28602.40402.53102.66772.81472.97273.14273.32533.52143.7320 通过比较知,角度较小时五次多项式算得的值较大,角度增大则两种方法得到的结果近似相等。 (2) %三次多项式插值 x=[149162536496481100]; y=[12345678910]; xi=1: 100; f=interp1(x,y,xi,'cubic') 运行结果: f= Columns1through11 1.00001.37291.71252.00002.24052.45512
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 下降 无约束 优化
![提示](https://static.bingdoc.com/images/bang_tan.gif)