本科毕业设计建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位全国大学生数学建模竞赛a题论文文档格式.docx
《本科毕业设计建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位全国大学生数学建模竞赛a题论文文档格式.docx》由会员分享，可在线阅读，更多相关《本科毕业设计建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位全国大学生数学建模竞赛a题论文文档格式.docx（31页珍藏版）》请在冰点文库上搜索。
1.3建立某太阳高度角γ下，杆影长的模型
由于杆影长大部分时候远小于地球的半径，所以不考虑地球地面的弯曲，当做地面为平面。
则影长shadowlength=L*cotγ。
其中L为杆长。
1.4北京时间和某一经度longitude真实时间的关系。
由于北京时间是以东经120°
为基准计算的，所以该经度longitude上的真实时间应该为realtime=（longitude-120）/15+beijingtime（小时），其中单位为小时和度，除以15因为地球每小时自转15°
2.画出2015年10月22日北京时间9:
00-15:
00之间天安门广场（北纬39度54分26秒,东经116度23分29秒）3米高的直杆的太阳影子长度的变化曲线。
根据上述模型，在MATLAB上编码得到该变化曲线。
编码参照附录
3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点。
将你们的模型应用于附件1的影子顶点坐标数据，给出若干个可能的地点。
附件1得到不同时刻的影长数据，我们采用梯度下降的方法，去拟合出一条曲线尽可能符合附件1的数据。
根据模型，影长函数
shadowlength（L，α，β，Θ）=L*cotγ
=L*cot（arcsin（sinαsinβ+cosαcosβcosΘ））
α可由日期2015年4月18日计算得到。
每个数据样本都可以计算出对应的Θ。
设数据样本有n个，第i个样本计算得到的Θ记为Θi，第i个样本的影子长度记为Lexami。
首先定义一个代价函数：
拟合即找到使代价函数达到尽可能小的值的参数L和β。
让cost对L和β分别求偏导数记为∂cost/∂L和∂cost/∂β。
（具体求导由MATLAB完成，在代码文件gradientdescent.m中实现）
初始化L，β为随机值，进行多次迭代每次迭代都修正参数L和β的值，以减小cost的值。
每次迭代修正方法：
L:
=L-stepl*∂cost/∂L；
β:
=β-stepb*∂cost/∂β；
则cost将随着L和β的修正逐渐变小，即拟合程度变好。
关于梯度下降法原理：
详见[1]。
其中stepl和stepb是控制修正速度的值，在机器学习中也叫学习速率，太大容易造成略过最优点，太小造成修正缓慢，效率低下。
此次计算中选用速率都为0.03，效果较好。
迭代次数为40000次。
并且为了防止出现局部最优解，试用了多次不同初始值。
最终得到的拟合结果，每个样本点影长差距远小于1毫米。
并且考虑了经度对时间的影响，试用了不同经度来拟合曲线。
得到了拟合结果类似如下：
由于差距很小，所以两条线几乎合并在一起，放大后如下图：
最终将拟合后得到的一些可能的经纬度制成下图：
如图，在东经100°
~125°
，北纬23.55°
都是可能的地点。
杆长为2.6073米。
3.根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。
将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。
由于数据的进一步缺失，没有交代日期，所以依旧采用梯度下降的方法拟合数据。
但是多了一个变量需要拟合就是太阳光和赤道面的夹角α。
对比上一题多了一个需拟合的变量，同样地求出cost关于α的偏导数∂cost/∂α（具体看MATLAB代码文件gradientdescent.m）。
迭代过程中增加每次对α的修正：
α:
=α-stepa*∂cost/∂α。
由于多出了变量，导致拟合的可能性变多，所以对结果进行筛选，需要满足cost<
0.008²
且拟合最终得到的参数|α|<
23°
26′；
|β|<
66°
34′，因为太阳光和赤道面所成角的范围所致。
考虑到附件2中影子在北京时间12点41分后直到13点41分仍慢慢减短，说明其只可能处于经度小于95度的范围内（这些地点的真实时间还没到12点，正在向12点靠近，所以影子在缩短）。
所以在东经50~95°
范围每隔5°
进行一次拟合。
拟合结果：
附件2可能的地点与日期以及杆长：
（7月1日，50°
E，42.71°
N，杆长0.846917m），
（7月13日，55°
E，40.19°
N，杆长0.90609m），
（7月24日，60°
E，37.15°
N，杆长0.980429m）
还有一些，具体看MATLAB代码的运行结果。
经度越靠近90度拟合效果越差。
同样地，附件3影子在13点09后影长数据不断增大，所以考虑拟合经度范围在东经110~120°
范围每隔1°
进行一次拟合（原本在105~150°
每隔5°
进行一次拟合，只有在110~120°
之间有两组拟合成功数据，所以改为110~120°
）。
拟合结果:
附件3可能的地点与日期以及杆长：
（9月18日，111°
E，37.95°
N，杆长3.059654m）
（9月7日，112°
E，42.47°
N，杆长3.027828m）
（8月27日，113°
E，46.27°
N，杆长3.005987m）
拟合较好的是在111°
E到113°
E之间。
4.附件4为一根直杆在太阳下的影子变化的视频，并且已通过某种方式估计出直杆的高度为2米。
请建立确定视频拍摄地点的数学模型，并应用你们的模型给出若干个可能的拍摄地点。
如果拍摄日期未知，你能否根据视频确定出拍摄地点与日期？
视频中给出了时间即相当于给出Θ，给出了杆长L，给出了日期相当于给出了α，所以唯一需要的就是纬度β。
用尺子粗略地测量了视频中影子的长度，得到以下数据：
杆长11厘米。
影长的变化情况：
8:
55、影长13厘米；
57、12.9厘米；
59、12.8厘米9:
01、12.6厘米；
9:
03、12.4厘米；
05、12.2厘米；
07、12厘米；
09、11.8厘米；
11、11.7厘米；
13、11.6厘米；
15、11.5厘米；
17、11.3厘米；
19、11.2厘米；
21、11厘米；
23、10.8厘米；
25、10.7厘米；
27、10.55厘米；
9：
29、10.4厘米；
31、10.3厘米；
33、10.2厘米
将数据输入程序，在经度为90°
对曲线拟合，得到β的情况，即可能的地点如下：
Possibleplace:
longitude:
90,latitude:
52.383329
91,latitude:
50.940848
92,latitude:
49.560709
93,latitude:
48.241217
94,latitude:
46.980504
95,latitude:
45.776579
96,latitude:
44.627373
97,latitude:
43.530781
98,latitude:
42.484696
99,latitude:
41.487046
100,latitude:
40.535827
101,latitude:
39.629140
102,latitude:
38.765238
103,latitude:
37.942608
104,latitude:
37.160098
105,latitude:
36.417195
106,latitude:
35.714660
107,latitude:
35.056272
108,latitude:
34.454905
109,latitude:
33.964138
110,latitude:
34.130393
111,latitude:
40.322309
112,latitude:
44.391538
113,latitude:
47.150567
114,latitude:
49.289467
115,latitude:
51.049417
116,latitude:
52.547724
117,latitude:
53.851881
118,latitude:
55.004924
119,latitude:
56.036277
120,latitude:
56.967101
121,latitude:
57.813207
122,latitude:
58.586770
123,latitude:
59.297392
124,latitude:
59.952799
125,latitude:
60.559310
其中数据都为东经和北纬。
如果没有给出日期，则相当于对两个参数α和β进行梯度下降逼近。
依然可以计算出可能的地点（纬度β和经度决定）和日期（太阳和赤道面所成角α决定）。
参考文献：
[1]woxincd，梯度下降法，
网址：
访问时间：
2015年9月13日
附录：
使用软件：
MATLAB2012a
以下A1、A2、A3、A4为各个小题的解决程序脚本，其余为脚本调用的各个函数。
A1.m
clearall;
clc;
%betaislatitude.
beta=（（39*60+54）*60+26）/（180*60*60）*pi;
%alphaissunverticleangleonequatorialplane.
alpha=alphaondate（10,22,0）;
%bjtime=[8:
0.05:
16];
bjtime=[9:
15];
%thetaisthetimeangle.
theta=thetafuncontime（BeijingChange2realTime（（116+（23+29/60）/60）,bjtime））;
%gammaissunverticleangleonearth.
gamma=asin（sinr（alpha,beta,theta））;
shadowlength=cot（gamma）*3;
plot（bjtime,shadowlength）;
xlabel（'
time'
）;
ylabel（'
lengthofstick'
'
sshadow'
A2.m
clc;
alpha=alphaondate（4,18,0）;
%inputexample
x=[1.0365
1.0699
1.1038
1.1383
1.1732
1.2087
1.2448
1.2815
1.3189
1.3568
1.3955
1.4349
1.4751
1.516
1.5577
1.6003
1.6438
1.6882
1.7337
1.7801
1.8277
];
y=[0.4973
0.5029
0.5085
0.5142
0.5198
0.5255
0.5311
0.5368
0.5426
0.5483
0.5541
0.5598
0.5657
0.5715
0.5774
0.5833
0.5892
0.5952
0.6013
0.6074
0.6135
bjtime=[14+42/60:
3/60:
15+42/60]'
;
%size（bjtime）
shadowlexample=（x.^2+y.^2）.^（1/2）;
%size（lexample）
beta=60/180*pi;
%40
l=2;
%31.52
longtitude=100;
%120116
i=1;
stepb=0.03;
stepl=0.03;
%（0.03,0.1）（0.03,0.07）alternative
forlontitude=100:
5:
125
%initializationdone.
realtime=BeijingChange2realTime（longtitude,bjtime）;
theta=thetafuncontime（realtime）;
[betal]=Descent2（l,alpha,beta,theta,shadowlexample,stepb,stepl,40000）;
latitude（i）=（beta/pi）*180;
%recordthedifferentlatitudecombinedwithdifferentlontitude
i=i+1;
l%lshowthelengthofstick.
shadowlength=l*cot（asin（sinr（alpha,beta,theta）））;
figure,plot（bjtime,shadowlength）;
%showthemodelpredictedlengthofshadow.
xlabel（'
beijingtime'
ylabel（'
shadowlength'
holdon;
plot（bjtime,shadowlexample,'
r'
end
%shadowlexample%showthelengthmeasurementofshadow.
figure,plot（[100:
125],latitude）;
%showthepossibleplace.
longtitude'
latitude'
A3.m
x1=[-1.2352
-1.2081
-1.1813
-1.1546
-1.1281
-1.1018
-1.0756
-1.0496
-1.0237
-0.998
-0.9724
-0.947
-0.9217
-0.8965
-0.8714
-0.8464
-0.8215
-0.7967
-0.7719
-0.7473
-0.7227
y1=[0.173
0.189
0.2048
0.2203
0.2356
0.2505
0.2653
0.2798
0.294
0.308
0.3218
0.3354
0.3488
0.3619
0.3748
0.3876
0.4001
0.4124
0.4246
0.4366
0.4484
shadowlexample1=（x1.^2+y1.^2）.^（1/2）;
bjtime1=[12+41/60:
13+41/60]'
x2=[1.1637
1.2212
1.2791
1.3373
1.396
1.4552
1.5148
1.575
1.6357
1.697
1.7589
1.8215
1.8848
1.9488
2.0136
2.0792
2.1457
2.2131
2.2815
2.3508
2.4213
y2=[3.336
3.3299
3.3242
3.3188
3.3137
3.3091
3.3048
3.3007
3.2971
3.2937
3.2907
3.2881
3.2859
3.284
3.2824
3.2813
3.2805
3.2801
3.2804
3.2812
shadowlexample2=（x2.^2+y2.^2）.^（1/2）;
bjtime2=[13+09/60:
14+09/60]'
alpha1=0/180*pi;
alpha2=0/180*pi;
beta1=60/180*pi;
beta2=45/180*pi;
l1=1;
l2=3;
longtitude=120;
%longtitude2shouldbelargerthan（120-41/60*15）degree.
%Becausethedataisturnlarger,therealtimeontheplace..
%shouldbeafter12o'
clock.
%whilelongitude1shouldbesmallerthan（120-15-41/60*15）
%forlongtitude=:
:
%initializationdone.
forlongtitude1=50:
5:
95
realtime1=BeijingChange2realTime（longtitude1,bjtime1）;
theta1=thetafuncontime（realtime1）;
stepa=0.003;
stepb=0.003;
stepl=0.003;
%%{
%[alpha1alpha2betal]=Descent3（l,alpha1,alpha2,beta,theta1,theta2,shadowlexample1,shadowlexample2,stepa,stepb,stepl,100000）;
[alpha1beta1l1goodflag]=Descent3（l1,alpha1,beta1,theta1,shadowlexample1,
升级会员 