清华大学谢金星数学实验作业4.docx
《清华大学谢金星数学实验作业4.docx》由会员分享，可在线阅读，更多相关《清华大学谢金星数学实验作业4.docx（15页珍藏版）》请在冰豆网上搜索。
清华大学谢金星数学实验作业4
实验7无约束优化
【实验目的】
1）掌握Matlab优化工具箱的基本用法，对不同算法进行初步分析、比较。
2）练习用无约束优化方法建立和求解实际问题的模型（包括线性和非线性最小二乘拟合）。
【实验内容】
7.5原子问题：
●模型建立：
设第i个点的坐标
，并假设第一个点的坐标
，于是有24个未知点，共48个未知数，而题给条件共可以列出52个方程，所以一般没有精确解，通过最小二乘法寻找最接近的解答，即需要以下代数式达到最小：
其中，
表示两个原子的距离，所以本题的本质既是以下优化：
●matlab实现
用lsqnonlin命令实现：
首先建立函数m文件，代码如下所示：
functionf=yuanzifun（X,A,d）
X（1,1）=0;
X（2,1）=0;
fori=1:
52f（i）=sqrt（（X（1,A（1,i））-X（1,A（2,i）））^2+（X（2,A（1,i））-X（2,A（2,i）））^2）-d（i）;
end
说明：
1.该步定义了每次的残差，其中X是一个25×2的矩阵，第一行表示25个原子的x坐标，第二行表示25个原子的y坐标；
2.为了简化编程过程，建立2×52的A矩阵，A的第一行表示题中所给的52组原子对的第一个原子的编号，A的第二行表示题中所给的52组原子对的第二个原子的编号，这样可以简化在命令文件中的代码数量；
3.d是52个元素的距离向量组。
命令文件如下所示：
clc,clearall
A=[4,12,13,17,21,5,16,17,25,5,20,21,24,...
5,12,24,8,13,19,25,8,14,16,20,21,14,...
18,13,15,22,11,13,19,20,22,18,25,15,17,...
15,19,15,16,20,23,18,19,20,23,24,23,23;...
1,1,1,1,1,2,2,2,2,3,3,3,3,...
4,4,4,6,6,6,6,7,7,7,7,7,8,...
8,9,9,9,10,10,10,10,10,11,11,12,12,...
13,13,14,14,16,16,17,17,19,19,19,21,21];
d=[0.96070.43990.81431.37651.27220.52940.61440.37660.68930.94880.80001.10901.1432...
0.47581.34020.70060.49451.05590.68100.35870.33510.28781.13460.38700.75110.4439...
0.83630.32080.15741.27360.57810.92540.64010.24670.47271.38400.43661.03071.3904...
0.57250.76600.43941.09521.04221.82551.43251.08510.49951.22771.12710.70600.8052];
[X,norm,res,exit,out]=lsqnonlin（@yuanzifun,zeros（2,25）,[],[],[],A,d）;
zuobiao=X’
cancha2fanshu=norm
plot（X（1,:
）,X（2,:
）,'*'）;
说明：
1.输入相应的参量A矩阵和d矩阵；
2.用lsqnonlin命令进行非线性最小二乘优化求解。
运行结果如下所示：
zuobiao=
00
0.53700.8806
-0.3115-0.0244
-0.48480.8215
-0.00190.8718
-0.49830.5170
-0.26590.4437
-0.00020.6228
0.21740.7903
-0.43620.5057
-0.34511.0775
0.3505-0.2350
0.51960.6374
-0.44110.6609
0.00600.7545
0.30341.4440
0.92571.0172
0.1734-0.2007
-0.13081.0749
-0.49480.7518
-0.80400.9745
-0.90940.2980
-1.07650.2699
-1.18680.7078
-0.15070.6820
cancha2fanshu=
0.0245
相应的函数坐标散点图如下所示：
如果换用其他的初值，比如将zeros（2,25）改成[0,ones（1,24）;0,ones（1,24）]，运行结果如下：
zuobiao=
00
0.90720.3760
0.02651.2465
0.54640.8434
0.38660.3823
0.61151.3615
0.27760.7549
0.21151.0935
0.45920.5827
-0.22450.3405
0.32700.5493
0.2148-0.4132
0.71330.3426
0.58060.8428
0.24670.6092
0.8095-0.2404
1.18840.6251
0.51061.8956
0.12910.8813
0.01490.4320
1.02470.7844
-0.70540.1692
1.16941.5337
1.18131.1846
0.53420.9701
cancha2fanshu=
0.0338
对比两次结果，有一定的差异，说明不同的初值会产生不同的优化结果，可能的解释是求和表达式本身的极值可能就不是唯一的，那么从不同的点出发可以得到不同的结果和不同的精度（本例中当初值为zeros（2,25）时候的精度更高，即更可信），具体选择哪组结果时候，可以多取几次初值，然后选择残差2范数最小的结果作为最可靠结果。
7.7生产函数问题
●模型建立1
因为
，（i=1,2…,19）
如果要进行非线性最小二乘拟合，只需要：
●matlab实现1：
可以直接调用lsqnonlin命令求解非线性最小二乘拟合；
函数M文件如下所示：
functionf=chanzhifun（x,Q,K,L）;
fori=1:
19
f（i）=x
（1）*K（i）^（x
（2））*L（i）^（x（3））-Q（i）;
end
说明：
x向量为待确定的数组，x
（1）=a，x
（2）=α，x（3）=β
命令文件如下所示:
clc,clearall
Q=[0.71710.89641.02021.19621.49281.69091.85482.16182.66383.46344.6759...
5.84786.78857.44637.83458.20689.94689.731510.4791];
K=[0.09100.25430.31210.37920.47540.44100.45170.55950.80801.30721.7042...
2.00192.29142.49412.84062.98543.29183.73144.3500];
L=[4.81794.98735.12825.27835.43345.53296.47496.54916.61526.68086.7455...
6.80656.89506.98207.06377.13947.20857.30257.3740];
x0=[1,0.5,0.5];
[x,norm,res,exit,out]=lsqnonlin（@chanzhifun,x0,[],[],[],Q,K,L）
plot（（1984:
2002）,K,'blue'）;
holdon;
plot（（1984:
2002）,L,'black'）;
plot（（1984:
2002）,Q,'red'）;
holdoff
gtext（'K'）;gtext（'L'）;gtext（'Q'）;
运行结果为
x=
0.83370.77350.7317
norm=
2.7488
res=
Columns1through8
-0.30460.04030.09990.13380.1253-0.1437-0.0865-0.0576
Columns9through16
0.15290.65270.4133-0.0456-0.2861-0.4392-0.0194-0.0217
Columns17through19
-1.05670.15710.7350
exit=
1
out=
firstorderopt:
2.1488e-007
iterations:
18
funcCount:
76
cgiterations:
0
algorithm:
'large-scale:
trust-regionreflectiveNewton'
message:
[1x427char]
另外从原始数据可以得到K,L,Q与时间的关系如下所示：
●模型建立2
如果要进行线性最小二乘拟合，对
两端取对数，得到
所以基函数是
，
，1,可以列出矩阵表达式：
简写为
●matlab实现
clc,clearall
Q=[0.71710.89641.02021.19621.49281.69091.85482.16182.66383.46344.6759...
5.84786.78857.44637.83458.20689.94689.731510.4791]';
K=[0.09100.25430.31210.37920.47540.44100.45170.55950.80801.30721.7042...
2.00192.29142.49412.84062.98543.29183.73144.3500]';
L=[4.81794.98735.12825.27835.43345.53296.47496.54916.61526.68086.7455...
6.80656.89506.98207.06377.13947.20857.30257.3740]';
f=[log（L）,log（K）,ones（19,1）];%构造φ函数，a的三个元素分别为β，α，lna
a=f\log（Q）;
aa=inv（f'*f）*f'*log（Q）;
a（3）=exp（a（3））;
aa（3）=exp（aa（3））;
a'%三个元素分别对应β，α，a
aa'%三个元素分别对应β，α，a
运行结果为
ans=
1.26770.66000.3148
ans=
1.26770.66000.3148
比较两种结果几乎没有差异，进一步印证了，左除命令如果得不到一个精确解，则会采用最小二乘法求近似解。
但是结果β大于1，这与题给的约束并不符合，需要限定范围，如下所示：
clc,clearall
Q=[0.71710.89641.02021.19621.49281.69091.85482.16182.66383.46344.6759...
5.84786.78857.44637.83458.20689.94689.731510.4791]';
K=[0.09100.25430.31210.37920.47540.44100.45170.55950.80801.30721.7042...
2.00192.29142.49412.84062.98543.29183.73144.3500]';
L=[4.81794.98735.12825.27835.43345.53296.47496.54916.61526.68086.7455...
6.80656.89506.98207.06377.13947.20857.30257.3740]';
f=[log（L）,log（K）,ones（19,1）];%构造相应的φ函数
[x,norm,res,exit,out]=lsqlin（f,log（Q）,[],[],[],[],[0,0,-5]',[1,1,0]'）
%限定ln（a）的范围在-5到1之间是根据之前得到的结果给出的大致区间
a=exp（x（3））
alpha=x
（2）
beta=x
（1）
得到如下结果：
x=
0.6433
0.7349
-0.0067
norm=
0.3575
res=
0.4242
-0.1301
-0.1692
-0.1717
-0.1350
0.0332
0.0069
-0.0045
-0.0723
-0.1697
-0.0707
0.0289
0.0705
0.0926
0.0403
0.0433
0.1576
0.0353
-0.0097
exit=
3
out=
iterations:
6
algorithm:
'large-scale:
trust-regionreflectiveNewton'
firstorderopt:
0.0115
cgiterations:
5
message:
[1x123char]
a=
0.9934
alpha=
0.7349
beta=
0.6433
这与非线性拟合的结果有一定差异，可以通过绘图比较两个拟合结果：
%拟合对比
clc,clearall
Q=[0.71710.89641.02021.19621.49281.69091.85482.16182.66383.46344.6759...
5.84786.78857.44637.83458.20689.94689.731510.4791]';
K=[0.09100.25430.31210.37920.47540.44100.45170.55950.80801.30721.7042...
2.00192.29142.49412.84062.98543.29183.73144.3500]';
L=[4.81794.98735.12825.27835.43345.53296.47496.54916.61526.68086.7455...
6.80656.89506.98207.06377.13947.20857.30257.3740]';
fori=1:
19
Q1（i）=0.8337*K（i）^0.7735*L（i）^0.7317;
Q2（i）=0.9934*K（i）^0.7349*L（i）^0.6433;
end
plot（（1984:
2002）,Q1,'blue'）;
holdon
plot（（1984:
2002）,Q2,'yellow'）;
plot（（1984:
2002）,Q,'red'）;
holdoff
得到的结果如下：
说明：
1.在所给的数据范围内二者的拟合效果是差不多的，如果要进一步比较，需要更多的数据。
2.残差2范数比较，虽然非线性拟合的norm=2.7488，线性的norm=0.3575，但是并不能说明前者的拟合效果比后者差，因为二者的量级是不一样的，前者是未对数化结果，后者是直接对数化以后得到的结果，所以没有太大可比性，此时从图形上观察更直接。
3.通过查阅相关的资料，并结合该问的相关结果，可以知道ａ表示一个国家的总体技术发展水平大小，技术水平越高的国家产值越大。
α在经济学上代表劳动力产出的弹性系数，即劳动力投入的变化引起产值的变化的速率。
同理β表示资本产出的弹性系数，即资产投入的变化引起产值变化的速率，
这是因为根据
得到
这正是弹性的定义，同理也有：
4.另外国际上一般取α＝0.8～0.6,β=0.2～0.4（来自‘维基百科’），而从本题给出的数据看，中国的劳动力产值弹性系数超过一般水平，猜测原因是中国的产值严重依赖于劳动力。
中国在初期的经济发展主要依靠劳动力密集型企业，本题的结论印证了这个结果。
【实验总结】
最小二乘优化是处理优化问题的重要方法，在工程，生产等领域的决策中十分方便，本次实验让我对优化问题有了初步认识，相信通过后续的学习会有更多收获。
升级会员 