实验七(无约束优化)

1、数学实验 实验七 无约束优化化21 张腾2012011864 2014-4-18大学数学实验 实验报告无约束优化一、 实验目的1、掌握MATLAB优化工具箱的基本用法，对不同算法进行初步分析、比较。2、练习用无约束优化方法建立和求解实际问题的模型（包括最小二乘拟合）。二、 实验内容项目一：某分子由25个原子组成，并且已经通过实验测量得到了其中某些原子队之间的距离（假设在平面结构上讨论），如表所示。请你确定每个原子的位置关系。原子对距离原子对距离原子对距离原子对距离(4,1)0.9607(5,4)0.4758(18,8)0.8363(15,13)0.5725(12,1)0.4399(12,4)1

2、.3402(13,9)0.3208(19,13)0.7660(13,1)0.8143(24,4)0.7006(15,9)0.1574(15,14)0.4394(17,1)1.3765(8,6)0.4945(22,9)1.2736(16,14)1.0952(21,1)1.2722(13,6)1.0559(11,10)0.5781(20,16)1.0422(5,2)0.5294(19,6)0.6810(13,10)0.9254(23,16)1.8255(16,2)0.6144(25,6)0.3587(19,10)0.6401(18,17)1.4325(17,2)0.3766(8,7)0.3351(

3、20,10)0.2467(19,17)1.0851(25,2)0.6893(14,7)0.2878(22,10)0.4727(20,19)0.4995(5,3)0.9488(16,7)1.1346(18,11)1.3840(23,19)1.2277(20,3)0.8000(20,7)0.3870(25,11)0.4366(24,19)1.1271(21,3)1.1090(21,7)0.7511(15,12)1.0307(23,21)0.7060(24,3)1.1432(14,8)0.4439(17,12)1.3904(23,22)0.8025问题分析：每个原子的位置都是未知的，在坐标系中只有相

4、对的位置参数，不妨固定原子1的坐标为（0，0）。并且分子可以在平面内任意旋转，很难确定每个原子的绝对位置，表中给出了52组数据，可以得到52个方程，而未知数个数为48个（每个原子的x,y坐标）。属于超定方程组，没有确定的解，只能求得最优解。故我采用最小方差标准来求最小二乘解得到最优解。模型建立：设每个原子的坐标为（xi，yi），则任意两个原子之间的距离可表示为：dij=xi-xj2-yi-yj2(i,j=1,225)则本题所要求的就是使方差：ijxi-xj2-yi-yj2-dij2最小的一组坐标值。所的到的坐标值即为最优的解。解决方案：根据题目中所给的数据，利用matlab编写函数M文件如下：

5、function s=yuanzif(x,d)x(1,1)=0;x(2,1)=0;d=0.9607 0.4399 0.8143 1.3765 1.2722 0.5294 0.6144 0.3766 0.6893 0.9488 0.8000 1.1090 1.1432 0.4758 1.3402 0.7006 0.4945 1.0559 0.6810 0.3587 0.3351 0.2878 1.1346 0.3870 0.7511 0.4439 0.8363 0.3208 0.1574 1.2736 0.5781 0.9254 0.6401 0.2467 0.4727 1.3840 0.436

6、6 1.0307 1.3904 0.5725 0.7660 0.4394 1.0952 1.0422 1.8255 1.4325 1.0851 0.4995 1.2277 1.1271 0.7060 0.8052;s(1)=(d(1)-(x(1,1)-x(1,4)2+(x(2,1)-x(2,4)2)0.5);s(2)=(d(2)-(x(1,1)-x(1,12)2+(x(2,1)-x(2,12)2)0.5);s(3)=(d(3)-(x(1,1)-x(1,13)2+(x(2,1)-x(2,13)2)0.5);s(4)=(d(4)-(x(1,1)-x(1,17)2+(x(2,1)-x(2,17)2)

7、0.5);s(5)=(d(5)-(x(1,1)-x(1,21)2+(x(2,1)-x(2,21)2)0.5);s(6)=(d(6)-(x(1,2)-x(1,5)2+(x(2,2)-x(2,5)2)0.5);s(7)=(d(7)-(x(1,2)-x(1,16)2+(x(2,2)-x(2,16)2)0.5);s(8)=(d(8)-(x(1,2)-x(1,17)2+(x(2,2)-x(2,17)2)0.5);s(9)=(d(9)-(x(1,2)-x(1,25)2+(x(2,2)-x(2,25)2)0.5);s(10)=(d(10)-(x(1,3)-x(1,5)2+(x(2,3)-x(2,5)2)0.

8、5);s(11)=(d(11)-(x(1,3)-x(1,20)2+(x(2,3)-x(2,20)2)0.5);s(12)=(d(12)-(x(1,3)-x(1,21)2+(x(2,3)-x(2,21)2)0.5);s(13)=(d(13)-(x(1,3)-x(1,24)2+(x(2,3)-x(2,24)2)0.5);s(14)=(d(14)-(x(1,4)-x(1,5)2+(x(2,4)-x(2,5)2)0.5);s(15)=(d(15)-(x(1,4)-x(1,12)2+(x(2,4)-x(2,12)2)0.5);s(16)=(d(16)-(x(1,4)-x(1,24)2+(x(2,4)-x

9、(2,24)2)0.5);s(17)=(d(17)-(x(1,6)-x(1,8)2+(x(2,6)-x(2,8)2)0.5);s(18)=(d(18)-(x(1,6)-x(1,13)2+(x(2,6)-x(2,13)2)0.5);s(19)=(d(19)-(x(1,6)-x(1,19)2+(x(2,6)-x(2,19)2)0.5);s(20)=(d(20)-(x(1,6)-x(1,25)2+(x(2,6)-x(2,25)2)0.5);s(21)=(d(21)-(x(1,7)-x(1,8)2+(x(2,7)-x(2,8)2)0.5);s(22)=(d(22)-(x(1,7)-x(1,14)2+(

10、x(2,7)-x(2,14)2)0.5);s(23)=(d(23)-(x(1,7)-x(1,16)2+(x(2,7)-x(2,16)2)0.5);s(24)=(d(24)-(x(1,7)-x(1,20)2+(x(2,7)-x(2,20)2)0.5);s(25)=(d(25)-(x(1,7)-x(1,21)2+(x(2,7)-x(2,21)2)0.5);s(26)=(d(26)-(x(1,8)-x(1,14)2+(x(2,8)-x(2,14)2)0.5);s(27)=(d(27)-(x(1,8)-x(1,18)2+(x(2,8)-x(2,18)2)0.5);s(28)=(d(28)-(x(1,9

11、)-x(1,13)2+(x(2,9)-x(2,13)2)0.5);s(29)=(d(29)-(x(1,9)-x(1,15)2+(x(2,9)-x(2,15)2)0.5);s(30)=(d(30)-(x(1,9)-x(1,22)2+(x(2,9)-x(2,22)2)0.5);s(31)=(d(31)-(x(1,10)-x(1,11)2+(x(2,10)-x(2,11)2)0.5);s(32)=(d(32)-(x(1,10)-x(1,13)2+(x(2,10)-x(2,13)2)0.5);s(33)=(d(33)-(x(1,10)-x(1,19)2+(x(2,10)-x(2,19)2)0.5);s

12、(34)=(d(34)-(x(1,10)-x(1,20)2+(x(2,10)-x(2,20)2)0.5);s(35)=(d(35)-(x(1,10)-x(1,22)2+(x(2,10)-x(2,22)2)0.5);s(36)=(d(36)-(x(1,11)-x(1,18)2+(x(2,11)-x(2,18)2)0.5);s(37)=(d(37)-(x(1,11)-x(1,25)2+(x(2,11)-x(2,25)2)0.5);s(38)=(d(38)-(x(1,12)-x(1,15)2+(x(2,12)-x(2,15)2)0.5);s(39)=(d(39)-(x(1,12)-x(1,17)2+

13、(x(2,12)-x(2,17)2)0.5);s(40)=(d(40)-(x(1,13)-x(1,15)2+(x(2,13)-x(2,15)2)0.5);s(41)=(d(41)-(x(1,13)-x(1,19)2+(x(2,13)-x(2,19)2)0.5);s(42)=(d(42)-(x(1,14)-x(1,15)2+(x(2,14)-x(2,15)2)0.5);s(43)=(d(43)-(x(1,14)-x(1,16)2+(x(2,14)-x(2,16)2)0.5);s(44)=(d(44)-(x(1,16)-x(1,20)2+(x(2,16)-x(2,20)2)0.5);s(45)=(

14、d(45)-(x(1,16)-x(1,23)2+(x(2,16)-x(2,23)2)0.5);s(46)=(d(46)-(x(1,17)-x(1,18)2+(x(2,17)-x(2,18)2)0.5);s(47)=(d(47)-(x(1,17)-x(1,19)2+(x(2,17)-x(2,19)2)0.5);s(48)=(d(48)-(x(1,19)-x(1,20)2+(x(2,19)-x(2,20)2)0.5);s(49)=(d(49)-(x(1,19)-x(1,23)2+(x(2,19)-x(2,23)2)0.5);s(50)=(d(50)-(x(1,19)-x(1,24)2+(x(2,1

15、9)-x(2,24)2)0.5);s(51)=(d(51)-(x(1,21)-x(1,23)2+(x(2,21)-x(2,23)2)0.5);s(52)=(d(52)-(x(1,22)-x(1,23)2+(x(2,22)-x(2,23)2)0.5);end利用lsqnonlin函数编写主程序求最优解并绘图进行直观的观察：clear all;d=0.9607 0.4399 0.8143 1.3765 1.2722 0.5294 0.6144 0.3766 0.6893 0.9488 0.8000 1.1090 1.1432 0.4758 1.3402 0.7006 0.4945 1.0559 0

16、.6810 0.3587 0.3351 0.2878 1.1346 0.3870 0.7511 0.4439 0.8363 0.3208 0.1574 1.2736 0.5781 0.9254 0.6401 0.2467 0.4727 1.3840 0.4366 1.0307 1.3904 0.5725 0.7660 0.4394 1.0952 1.0422 1.8255 1.4325 1.0851 0.4995 1.2277 1.1271 0.7060 0.8052;x0=0*ones(2,25);x,norms=lsqnonlin(yuanzif,x0,d);for i=1:1:25 pl

17、ot(x(1,i),x(2,i),'bo'); hold on;end运行结果如下：原子xy10.00000.000020.82910.307130.40980.048541.00150.251351.09520.730560.03421.015370.22460.463380.42720.733390.6652-0.0238100.12790.948311-0.03310.357312-0.26510.3835130.80720.2967140.37440.2402150.5733-0.1819161.35930.6456171.19060.4671180.39931.629

18、4190.09110.3446200.32770.8217210.39531.1878220.51031.235223-0.23431.5277241.05360.9621250.24430.7030原子分布图像如下（蓝色）：改变初值x0=1*ons(2,25)得到另一组数据（不再一一列出），图像如下（红色）：将两组数据在图像上对比发现二者差距很大，说明改变迭代初值将改变结果。这也是合理的，首先原子的位置信息需要至少两两原子间的距离数据才能准确确定，也就是需要25*24/2组相对距离数据，而本题中只有52组，显然解是有无穷多个的，再加上实际原子无论怎样旋转都不会改变相对位置，确定这些原子的绝对

19、位置至少还需要三个原子的绝对位置信息才能确定出所有原子的绝对位置信息。项目二：经济学中著名的Cobb-Douglas生产函数的一般形式为QK，L=aKL，0<，<1其中Q，K，L分别表示产值、资金、劳动力，式中，要由经济统计数据确定。现有中国统计年鉴(2003)给出的统计数据如表2所示，请分别用线性和非线性最小二乘拟合求出式中的，并解释，的含义。年份总产值/万亿元资金/万亿元劳动力/亿人19840.71710.09104.817919850.89640.25434.987319861.02020.31215.128219871.19620.37925.278319881.49280

20、.47545.433419891.69090.44105.532919901.85480.45176.474919912.16180.55956.549119922.66380.80806.615219933.46341.30726.680819944.67591.70426.745519955.84782.00196.806519966.78852.29146.895019977.44632.49416.982019987.83452.84067.063719998.20682.98547.139420009.94683.29187.208520019.73153.73147.3025200

21、210.47914.35007.3740问题分析：本题目给出了19组数据，根据给出的关系式可以列出19个方程，其中有3个未知量，属于超定方程组，可以求出最优解。题目要求用线性最小二乘和非线性最小二乘寻找最优解，故用两种不同的方法求解，可以得到两组不同的答案。模型建立：对于线性最小二乘最优解，由于原方程是非线性的，故首先对原等式进行数学处理使其关于未知参数线性化：QK，L=aKL同时取自然对数得：lnQK，L=lna+lnK+lnL由此继而得到方程组：lnQK(i)，L(i)=lnai+lnKi+lnLi; i=1,2,3,19对于非线性最小二乘求解，则可以直接列出方程：QK(i)，L(i)=a

22、(i)K(i)L(i); i=1,2,3,19解决方案：对于线性化最小二乘求最优解，直接利用matlab中的左除命令编写程序运算：q=0.7171,0.8964,1.0202,1.1962,1.4928,1.6909,1.8548,2.1618,2.6638,3.4634,4.6759,5.8478,6.7885,7.4463,7.8345,8.2068,9.9468,9.7315,10.4791'k=0.0910,0.2543,0.3121,0.3792,0.4754,0.4410,0.4517,0.5595,0.8080,1.3072,1.7042,2.0019,2.2914,2.

23、4941,2.8406,2.9854,3.2918,3.7314,4.3500'l=4.8179,4.9873,5.1282,5.2783,5.4334,5.5329,6.4749,6.5491,6.6152,6.6808,6.7455,6.8065,6.8950,6.9820,7.0637,7.1394,7.2085,7.3025,7.3740'lnq=log(q);lnk=log(k);lnl=log(l);A=ones(19,1) lnk lnl;b=lnq;x=Ab;ans=exp(x(1) x(2) x(3)运行得到结果：a=0.3148，=0.6600，=1.2677 但是其中=1.2677不满足题中对的要求，因此这组解没有实际意义。对于非线性最小二乘求最优解，首先编写函数M文件如下：function f=shengchan(x,k,l,q)f=x(1).*k.(x(2).*l.(x(3)-q;end编写主脚本文件如下：clear all;q=0.7171,0.8964,1.0202,1.1962,1.4928,1.6909,1.8548,2.1618,2.6638,3.4634,4.6759,5.8478,6.7885,7.4463,