清华大学谢金星数学实验作业4.docx

清华大学谢金星数学实验作业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（来自‘维基百科’），而从本题给出的数据看，中国的劳动力产值弹性系数超过一般水平，猜测原因是中国的产值严重依赖于劳动力。

中国在初期的经济发展主要依靠劳动力密集型企业，本题的结论印证了这个结果。

【实验总结】

最小二乘优化是处理优化问题的重要方法，在工程，生产等领域的决策中十分方便，本次实验让我对优化问题有了初步认识，相信通过后续的学习会有更多收获。