数学建模竞赛A题论文-建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位.doc

建立模型计算太阳影子长度和梯度下降法拟合太阳影子数据定位1 摘要本文建立了一个较为理想的几何立体模型来计算某个日期，某个时间，某个纬度的太阳高度角，从而计算太阳照射直杆的影长。并且通过对一系列按时间变化记录的数据以及部分已知参数，利用梯度下降法进行曲线拟合样本数据点，以此求出相应的未知的参数，算出数据记录时所在的可能的地点经纬度和日期。关键词：几何立体模型 拟合曲线 梯度下降法 2 正文1. 建立太阳光照射地球，地球上垂直地面杆的影长模型 假设太阳光线完全平行，地球是完全标准的球体，地球公转轨道是标准圆形，且随日期变化匀速运转。1.1 建立太阳光线和地球上某一纬度的上某点X切面所形成的太阳高度角模型O为地球球心。设太阳光线和地球的赤道面形成的角度为EOC = ，纬度即BOA=，根据纬度的定义显然OB是纬度上X点的切面（即X点地面）的法向量方向。X点的经度和太阳垂直照射点所在的经度的夹角即AOC为（可以由该点的真实时间计算得到，地球每小时自转15，=（该点真实时间-12） * 15 ）。欲求太阳在X点的高度角，相当于求X点地面法向量方向OB和太阳光线的夹角BOE的余角。假设太阳高度角为，则sin = cos BOE。如图作ACOC，作EC和BA垂直于赤道面，BDCE于D。设OA长为z，则OB = z / cos，OE=z cos / cos，BD=AC=z sin，CD=AB=z tan，CE=z cos tan。DE = CE - CD =z cos tan - z tan。BE = BD + DE=（z sin）+（z cos tan - z tan）。根据余弦定理可得，经化简可得：cosBOE = sinsin + coscoscos。即sin = sinsin + coscoscos用反三角函数即可知太阳高度角。1.2 建立太阳光线和地球赤道面所成角度随日期变化的模型如图，如果不考虑地球自转，则随着公转，太阳直射点将在地球上画出一个圆，这个圆所在的平面即为回归面。设A点为夏至日（6月22日）的直射点，此时AOB=N=2326，随着日期推移，直射点移到点C，则此时COD=正是我们所要求的太阳光线和地球赤道面所成的角，而AOC则是日期相对6月22日（公转导致）移动所形成的夹角，所以AOC=相对6月22日天数差距/365 * (闰年为366)，其中比6月22日小则天数差距为负值。易证，赤道面和回归面的二面角大小为N=2326，平面AOB赤道面，平面COD赤道面。所以该问题可以转化为如下图所示的情形。即一个二面角AOB大小为N，垂直于二面角AOB其中一个面BOD的面COD，求COD被这个二面角所夹出来的COD（即）大小。图中AOB，AOC，COD为直角三角形。CD = AB = OA sinN，OC=OA / cosAOC，sin=CD / OC = sinN * cosAOC用反三角函数即可求得。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（sinsin + coscoscos）可由日期2015年4月18日计算得到。每个数据样本都可以计算出对应的。设数据样本有n个，第i个样本计算得到的记为i，第i个样本的影子长度记为Lexam i。首先定义一个代价函数：拟合即找到使代价函数达到尽可能小的值的参数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毫米。并且考虑了经度对时间的影响，试用了不同经度来拟合曲线。得到了拟合结果类似如下：由于差距很小，所以两条线几乎合并在一起，放大后如下图：最终将拟合后得到的一些可能的经纬度制成下图：如图，在东经100125，北纬23.55都是可能的地点。杆长为2.6073米。3. 根据某固定直杆在水平地面上的太阳影子顶点坐标数据，建立数学模型确定直杆所处的地点和日期。将你们的模型分别应用于附件2和附件3的影子顶点坐标数据，给出若干个可能的地点与日期。由于数据的进一步缺失，没有交代日期，所以依旧采用梯度下降的方法拟合数据。但是多了一个变量需要拟合就是太阳光和赤道面的夹角。对比上一题多了一个需拟合的变量，同样地求出cost关于的偏导数cost / （具体看MATLAB代码文件gradientdescent.m）。迭代过程中增加每次对的修正：:= - stepa * cost / 。由于多出了变量，导致拟合的可能性变多，所以对结果进行筛选，需要满足cost0.008且拟合最终得到的参数|2326；| 23+26/60 | alpha1d 66 + 34/60 | beta1d 23+26/60 | alpha2d 66 + 34/60 | beta2d M(month) day = day-M(month); month = month + 1;end month = month - 1; end Descent2.mfunction b l = Descent2( l,a,b,e,shadowlexample,stepb,stepl,itertimes )%UNTITLED5 Summary of this function goes here% Detailed explanation goes herefor i = 1:itertimes cost = costFunc(l,a,b,e,shadowlexample); % to show cost in the descent process when debugging . %and deciding the parameters of gradient descent. dl da db = gradientdescent(l,a,b,e,shadowlexample); b = b - db*stepb; l = l - dl*stepl;end end Descent3.mfunction a b l goodflag = Descent3( l,a,b,e,shadowlexample,stepa,stepb,stepl,itertimes )%UNTITLED5 Summary of this function goes here% Detailed explanation goes herefor i = 1:itertimes cost = costFunc(l,a,b,e,shadowlexample); % to show cost in the descent process when debugging . %and deciding the parameters of gradient descent. dl da db = gradientdescent(l,a,b,e,shadowlexample); b = b - db*stepb; l = l - dl*stepl; a = a - da*stepa; end if(cost 0.0082) goodflag = 1; else goodflag = 0; end end gradientdescent.mfunction gradl grada gradb = gradientdescent( l,a,b,e,lexample ) %length,alpha,beta,theta.theta is calculated by timeexample.%UNTITLED2 Summary of this function goes here% Detailed explanation goes hereexamplenumber = size(lexample,1);shadowlength = l*cot(asin(sinr(a,b,e); dl = (1 - (sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2).(1/2). ./(sin(a)*sin(b) + cos(a)*cos(b)*cos(e); % dshadowlength/dl gradl = (1/examplenumber) * sum( 2*(shadowlength - lexample) .* dl , 1); % dcost/dl da = - (l*(cos(a)*sin(b) - cos(b)*cos(e)*sin(a). ./(1 - (sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2).(1/2). - (l*(cos(a)*sin(b) - cos(b)*cos(e)*sin(a).*(1 - (sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2).(1/2). ./(sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2; % dshadowlength/da grada = (1/examplenumber) * sum( 2*(shadowlength - lexample) .* da , 1); % dcost/da db = - (l.*(cos(b)*sin(a) - cos(a)*cos(e)*sin(b). ./(1 - (sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2).(1/2) . - (l*(cos(b)*sin(a) - cos(a)*cos(e)*sin(b).*(1 - (sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2).(1/2). ./(sin(a)*sin(b) + cos(a)*cos(b)*cos(e).2; % dshadowlength/db gradb = (1/examplenumber) * sum( 2*(shadowlength - lexample) .* db , 1); % dcost/db end sinr.mfunction sinr = sinr( a,b,e )%alpha,beta,theta%alpha sun vertical angle on equatorial planebeta latitude angletheta time angle% %UNTITLED Summary of this function goes here% Detailed explanation goes here%cosr = (cos(e)2/cos(a)2 + 1/cos(b)2 - (tan(b) - cos(e)*tan(a)2 - sin(e)2 ). /(2 * cos(e)/cos(a) * 1/cos(b)%cosr = (cos(e)2*cos(b)2 + cos(a)2 -(sin(b)*cos(a)-cos(e)*sin(a)*cos(b)2 - sin(