基于MATLAB的古塔变形情况数学模型全国大学生数学建模大赛.docx
《基于MATLAB的古塔变形情况数学模型全国大学生数学建模大赛.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的古塔变形情况数学模型全国大学生数学建模大赛.docx（79页珍藏版）》请在冰豆网上搜索。
基于MATLAB的古塔变形情况数学模型全国大学生数学建模大赛
2013高教社杯全国大学生数学建模竞赛
承诺书
我们仔细阅读了《全国大学生数学建模竞赛章程》和《全国大学生数学建模竞赛参赛规则》（以下简称为“竞赛章程和参赛规则”，可从全国大学生数学建模竞赛网站下载）。
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛章程和参赛规则的，如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛章程和参赛规则，以保证竞赛的公正、公平性。
如有违反竞赛章程和参赛规则的行为，我们将受到严肃处理。
我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。
我们参赛选择的题号是（从A/B/C/D中选择一项填写）：
C
我们的参赛报名号为（如果赛区设置报名号的话）：
5448
所属学校（请填写完整的全名）：
广东科学技术职业学院
参赛队员（打印并签名）：
指导教师或指导教师组负责人（打印并签名）：
康海刚
（论文纸质版与电子版中的以上信息必须一致，只是电子版中无需签名。
以上内容请仔细核对，提交后将不再允许做任何修改。
如填写错误，论文可能被取消评奖资格。
）
日期：
2013年9月13日
赛区评阅编号（由赛区组委会评阅前进行编号）：
2013高教社杯全国大学生数学建模竞赛
编号专用页
赛区评阅编号（由赛区组委会评阅前进行编号）：
赛区评阅记录（可供赛区评阅时使用）：
评
阅
人
评
分
备
注
全国统一编号（由赛区组委会送交全国前编号）：
全国评阅编号（由全国组委会评阅前进行编号）：
基于MATLAB的古塔变形情况数学模型
摘要
本文以MATLAB软件为核心工具，几何知识为基础，分别建立古塔的几何中心模型、变形情况模型以及变形预测模型。
在模型建立之前，首先对观察数据进行分析观察，以此确定建模目标与建模基础，通过观察发现，古塔塔身为近似正八棱台结构，从而为后续几何中心模型的建立提供了几何基础。
对数据整体观察时发现，对于86年与96年的第13层存在数据为空的现象，依据前面对古塔形状的定性，利用几何知识对缺失数据进行修复，便于后面的统计分析。
在之后对于观测数据的观察时，通过查阅相关资料文献，并结合实际判断，得出观测的数据点位置为古塔各层的层底数据，从而为后续确定各层的几何中心奠定基础。
对于几何中心模型，首先对于几何中心进行定义，对于该古塔，其形状为近似正八棱台，因此将其各层几何中心定义为各层的层底与层顶的几何中心点的连线中点为该层的几何中心点。
基于该定义，首先需要求出各层层底的几何中心，在求解的过程中，使用了LINGO软件，MATLAB中的直接搜索工具箱等予以辅助，以此为基础求解各层的结合中心。
对于变形情况模型，题目给出了三大方向，即分别从古塔的倾斜、弯曲和扭曲的角度进行考虑，因此在模型的建立过程中，分别对古塔的倾斜度、弯曲度、扭曲度进行合理的定义与相应的求解。
对于倾斜度，将其定义为古塔的倾斜位移值与塔高之比；对于弯曲度，定义为几何中心曲线与倾斜线的积分差，在具体建模时，对弯曲度模型进行扩展，考虑了各层的弯曲度，并利用SPSS软件采用聚类分析对弯曲度进行分级；对于扭曲度，定义其为各层几何中心散点的投影点与采用最小二乘法拟合的直线的距离平方和。
基于对三个指标的定义，分别对其建立相应的数学模型，并在求解模型结果后对结果进行分析与检验。
对于古塔变形预测模型，考虑到数据的小样本与不等时距的特点，建立不等间距灰色预测模型，弥补了小样本与不等时距的缺点。
模型建立时，首先利用前面问题的求解结果：
倾斜度等，并对数据进行归一化处理，以此建立模型后，再利用MATLAB软件编写相应代码对模型的结果进行求解，并对求解结果进行检验从而判断模型的合理性。
之后又采用BP神经网络算法同样求解模型结果，通过两组结果的相近性比较，从而验证灰色预测模型的合理性。
之后再利用模型对未来几年的变形情况进行预测，从而判断古塔的变形趋势。
在本文的最后，对于模型的改进和推广、优缺点分别进行了讨论，以此对模型进行评价与拓展。
关键词：
古塔变形，MATLAB，直接搜索工具箱，最小二乘法，不等间距灰色预测
1问题重述
C题古塔的变形
由于长时间承受自重、气温、风力等各种作用，偶然还要受地震、飓风的影响，古塔会产生各种变形，诸如倾斜、弯曲、扭曲等。
为保护古塔，文物部门需适时对古塔进行观测，了解各种变形量，以制定必要的保护措施。
某古塔已有上千年历史，是我国重点保护文物。
管理部门委托测绘公司先后于1986年7月、1996年8月、2009年3月和2011年3月对该塔进行了4次观测。
请你们根据附件1提供的4次观测数据，讨论以下问题：
1.给出确定古塔各层中心位置的通用方法，并列表给出各次测量的古塔各层中心坐标。
2.分析该塔倾斜、弯曲、扭曲等变形情况。
3.分析该塔的变形趋势。
2模型假设
1.假设2011年之后，不发生严重影响古塔变形的天灾人祸，从而保证预测模型的准确性。
2.附件提供的观测数据合理，即测量方法正确。
3.假设附件提供的不同年份观测数据观测点相同，即为同一空间坐标系。
3符号说明
符号
符号说明
e
观测点
min
各个散点到中心点的距离之和
d
各个散点到中心点的距离
U
层底的散点
δ
倾斜位移值
α
古塔倾斜度
H
塔高
λ
弯曲度
X
x轴坐标
Y
y轴坐标
Z
z轴坐标
4问题分析
4.1背景分析[2]
大型、特大型建筑物在其施工与运营过程中，由于受各种因素的影响，必将导致建筑物产生三维变形，如果变形超过允许的限度，就会影响建筑物的正常使用，严重的还会导致建筑物的破坏，给国家和人民的生命财产造成巨大损失。
本论文涉及的古塔，历经千年，是国家不可多得的宝贵财富，但由于多年来遭受气温、风力等作用使得古塔产生了变形，如若放任不管，将会有坍塌的危险，因此急需对古塔进行修缮工作，整修过程需要掌握现状的各种资料，例如各层的中心点、倾斜程度和扭曲程度等，本论文基于MATLAB等软件利用几何知识建立相应的数学模型，解决相关问题，从而保证重修工程的顺利实施。
4.2建模求解分析
4.2.1分析问题与建模前的准备——数据整理与分析
在对问题进行全面分析之前，需要对附件提供的数据进行整理与分析，从而保证分析问题时的数据基础可靠性。
4.2.1.1数据观察
附件提供了1~13层每一层其中8个点的观察数据，选取1986年第1层的8个点导入MALTAB软件中，并将各点连接起来，以俯视角度观察，得到图形如下：
图1-1：
测量数据点俯视图
观察图形，可以明显的看出，所观察的点构成的图形近似为一个正八边形，将其他各层的数据依次代入观察，发现图形结果相近，初步认定题目所涉及的古塔的形状应为近似八棱台的八角形塔。
在对整体数据进行观察时，发现1986年与1996年关于塔尖的观测点有4个，而到了2009年与2011年时仅剩下一个，考虑到产生这种结果的原因可能是由于在96年和09年之间由于塔尖的变形或缺失导致观测点无法再进行观测，从而在09年开始仅选择一点作为塔尖的观测数据，由于塔尖形状未知，无法对其进行修复等操作，故在考虑86年至96年的塔尖时仅取其中心点，并认定09年与11年的观测点即为其塔尖的中心点。
4.2.1.2数据修复
通过观察附件中的数据，发现1986年与1996年的观察数据中，关于第13层第5个点观测数据为空，而2009年与2001年则不为空，原因判断[1]：
在当时，测量塔的坐标数据时，需要在测量点安置反射片或棱镜，而第13层第5个点在当时可能不便于安置，导致无法测量，而后来（2009年之后）由于可采用无棱镜全站仪进行观测，才拥有测量数据结果。
虽然附件未能给出该点的精确测量结果，但通过上文描述即对数据的观察发现塔横截面近似为八边形，因此可利用几何知识拟合缺失点的X、Y坐标，这里是需要特别注意的，我们知道，通过将一层的8个点数据相连接所形成的图在三维空间里并不是一个平面图，而是一个凹凸不平的图，但其俯视图确为近似正八边形，因此我们通过八边形结合几何知识修复的仅是缺失点的X,Y坐标即平面坐标，对于高度即Z坐标，则通过第13层的其它各点与第12层的高度差关系，得出的缺失点的Z坐标，分别计算86年与96年的相应缺失点坐标，从而实现修复数据的目的。
具体步骤如下：
以86年为例，利用MATLAB软件作出第13层剩余7个点的坐标图的俯视图，并标号，结果如下：
图1-2：
第13层测量数据点俯视图
将补齐第5个点之后的图形理想化为正八边形，则缺失点5即为点1关于点3与点7之间连线所形成的直线的对称点，如下图：
图1-3：
数据修复操作图
点1的坐标为e1（566.308,525.092），点3的坐标为e3（564.481,521.521）,点7的坐标为e7（569.897,523.188），设点3与点7的连接直线为l,分别代入e3和e7，得直线l的方程可表示为
，设点5的坐标为（m,n），根据对称关系，（m,n）满足
（1）其中（a,b）为点1坐标，A表示直线l方程中的x的系数0.3078，B表示y的系数-1，C等于347.7784，代入相应数据，解得点5（m,n）=（567.9972,519.6039）
在修复了平面坐标后，还需要对Z坐标进行修复，依据上文提到的，分别计算第13层剩余7个点与第12层的对应点的Z坐标之差，所得的高度差平均值可作为第12层第5个点与第13层第5个点的高度差，由于第12层数据已知，从而得出缺失点的Z坐标。
具体计算如下：
第13层剩余7个点与第12层对应点的高度差如下表：
表1-1：
第13层剩余点与第12层高度差表
点
高度差
1
4.123
2
4.122
3
4.121
4
4.123
6
4.123
7
4.123
8
4.123
平均
4.123
通过上表可以看出，第13层各个点与第12层的对应点的高度差十分接近，符合实际情况，计算平均高度差与第12层第5个点之和即为第13层缺失点坐标，结果为48.713+4.123=52.836，从而得到修复点的三维坐标为（567.9972,519.6039，52.836），同理求得96年第13层的坐标为（568.0059,519.5896,52,.835）。
为了保证严谨性，还需对修复后的数据进行检验，判断其修复结果是否可信，将修复后的86年第13层的点与其他各点重新构图，分别计算两两各点之间的距离，如下表：
表1-2：
86年第13层各点距离图
边（点，点）
距离（m）
1,2
2.1710
2,3
2.1082
3,4
2.1663
4,5
2.1076
5,6
2.2387
6,7
2.1489
7,8
2.0976
8,1
2.2904
通过上表可以看出，本就完整的点7、点8、点1，这三个点之间的两两距离差为2.2904-2.0976=0.1928，反观修复的点5与点4、点3的距离为2.1076与2.2387，这两个距离与其他距离差距并不大，属于合理的范围，说明修复结果符合实际，可以采用。
4.2.1.3数据的定性
由于题目并没有明确指出各层的测量的数据具体是各层的哪个部分，因此在利用数据进行问题分析与计算之前，还需要对数据进行定性，即确定测的究竟是塔的哪些部分。
通过查阅参考文献，发现在对古塔的实际测量中，需要测量的是塔的各层的层底或层顶，因此需要引入层底和层顶的概念，具体如下图：
图2：
层底、层顶说明图
由此可以得出，附件给出的数据均为各层的层顶或层底的点的坐标，在给出的数据中，可以发现对于同一层，Z坐标相差并不大，所以可以认为实际上的层底或层顶是“很薄”的一层结果，但并不能因为高度很小而将其忽略从而将其视为平面，在本论文处理过程中，仍将其理解为测量数据是在一个三维空间中的测量，而不是一个平面。
针对第1层，假设其测量数据为第1层层顶的数据，那么在第1层所给数据中，Z坐标均为1.8m左右，说明第1层的层高仅为1.8m，这显然不符合实际，排除是层顶数据后，则证明附件给出的各层数据均为每层的层底数据，即1.8m为该塔的塔基，符合实际。
由于上一层的层底即可看作本层的层顶，因此测量数据仍可描述为层顶或层底的数据。
4.2.2分析问题一求各层的中心位置
题目要求确定古塔各层中心位置，并给出使用于所有其他层的通用方法，即要以一种计算模型的方式确定各层中心位置，首先通过数据整理得到的图，判断塔的横截面为近似正八边形，通过几何知识可知，正多边形任意两条边中垂线的交点即为正多边形的几何中心。
针对本题情况，由于本题涉及的是近似八棱台的三维立体的几何中心，因此分为二个步骤求解该立体的几何中心：
步骤一：
求解各层层底的中心点。
步骤二：
根据步骤一求得的两坐标，计算两坐标连线的中点即为该层的几何中心点。
对于步骤一，由于附件给出的对于一层塔的测量数据并非是在同一平面上的点，因此需要取任意三点构成一个平面，在这一平面中取任意两点的连线的中垂线与另外两点连线的中垂线的交点，重复遍历所有的任意三点，得到的交点集合在立体上的分布应是“一堆”散点，此时的每个散点都存在是该几何体（层）的几何中心点的可能，因此还需通过软件求解最优点，最终得出的最优中心点即为该层的几何中心点。
在求解最优点的问题上，我们采用了两种方法进行求解：
方法一：
利用LINGO软件，编写代码，求层底的三维空间中存在的一个点，使得这个点到其它中心点的距离最小，则该点即为所求层底的中心点。
方法二：
利用MATLAB直接搜索工具，采用模式搜索算法，计算中心点的散点集合中的最优点。
再根据两种方法求得的结果，进行比较。
对于步骤二：
结合上文数据定性可知，第二层的层底即为第一层的层顶，因此层底的中心点与层顶的中心的连线的中点极为该层的几何中心点，需要注意的是，其中对于第13层，由于没有第14层的存在，因此采用塔尖的中心点与第13层层底的中心点的连线中点为第13层的中心位置。
对于以上思路与操作步骤，作出流程图如下：
图3：
问题一求解流程图
4.2.3分析问题二塔的变形情况
题目提及的塔的变形情况有三种，分别是塔的倾斜情况、塔的弯曲情况、塔的扭曲情况。
对于塔的倾斜情况，基于问题一中的结果，得到塔底与塔顶的中心坐标的具体值，利用两坐标的平面X、Y坐标的距离计算倾斜位移值，倾斜位移值与塔顶高度之比即为该塔的倾斜度。
对于塔的弯曲情况的判断，是基于倾斜函数上进行的判断，可推论，若古塔仅发生倾斜并未发生弯曲，那么古塔各层的中心坐标的连线在三维空间中应是一条直线，而在实际情况中，由于古塔发生了弯曲，导致该连线是一条弯曲的曲线，基于此，对于古塔的弯曲程度可采用该曲线与倾斜直线之间的区别来说明古塔的弯曲度，对于弯曲度，定义其为：
以古塔倾斜方向为横轴，以Z坐标为纵轴形成的二维平面中，中心坐标所连曲线与倾斜线之间的差异，称为弯曲程度，差异越大说明弯曲越大，对于弯曲度的计算，择采用积分的思想，计算在二维图中几何中心曲线与倾斜线的积分之差，定义该差值为弯曲度。
对于塔的扭曲情况，不同于弯曲度，在几何中心坐标所在三维空间中所形成的散点图中，作出该空间图的俯视图即这些散点在平面上的投影图，所形成的图是分布在平面上的由中心点构成散点图，对这些散点利用最小二乘法的思想拟合出一条直线，求各个散点与这条拟合直线的离差平方和即为偏离程度，由此可判断塔的扭曲情况。
对于扭曲度的定义是：
在古塔各层中心坐标所在三维空间中，作出该三维空间的俯视图，得到以X,Y坐标为轴的平面图，判断各层中心坐标对由最小二乘法拟合的中心坐标点的直线的偏离程度，由此判断古塔的扭曲程度，可知，几何中心点与拟合越偏离，则说明古塔越扭曲，对于扭曲程度的计算，则计算各层中心点与拟合直线的距离的平方和，将其值定义为各层扭曲度，所有层的扭曲度之和则为整个塔的扭曲度。
对于古塔的倾斜、弯曲和扭曲程度的处理上，其基本思想都是将一个三维问题转化为一个二维问题进行求解，利用二维空间上的变化程度来描述三维空间的变化情况，即利用倾斜度，弯曲度，扭曲度三者结合来判断古塔的三维变形情况。
4.2.4分析问题三古塔的变形趋势
由题意可知，对于塔的变形趋势的分析即是要求建立一个预测未来塔的变形情况的预测模型，在建立预测模型时，需要特别考虑的是附件所提供数据的特点。
针对变形趋势的预测，理论上应该对于倾斜度，弯曲度以及扭曲度三者都进行预测处理，但通过对于问题一和问题二的求解结果的分析发现，对于扭曲度，存在特殊数据的可能：
86年与96年、09年与11年的扭曲度无太大变化，而96年与09年却存在十分大的区别，86年与96年相差10年，96年与09年相差13年，年份相差不大，但扭曲度的变化却出现很大差别，排除时间长度的影响，考虑扭曲的在96年与09年存在很大差别的原因可能是由于人工或天灾造成的，而这属于特殊情况，在有特殊情况影响的条件下作出对于扭曲度的变化趋势的预测显然是不合理的，因此对于古塔的变形趋势的考虑时，仅针对倾斜度和弯曲度进行处理。
由于附件仅提供了4个年份的数据，因此数据具有小样本数据的特点，同时4个年份1986、1996、2009、2011并非是等距的数据，所以数据具有不等距、小样本的特点。
针对样本数量小，建立预测模型时首先灰色模型进行预测，灰色模型具有应对小样本数据仍能较为精确的预测数据的特点，但缺点是灰色预测针对的是等距的数据样本，由于现实样本存在不等距的情况，所以需要对于样本数据化为等距数据进行处理，从而适应灰色模型的条件，得出利用灰色预测得到的结果。
考虑到模型的严谨性，还需对灰色模型求得的结果进行检验，因此可建立另一种预测模型——神经网络BP算法预测，比较两种模型求解的结果，检验灰色预测模型的合理性。
5模型建立与求解
5.1古塔几何中心模型的建立与求解:
依据上文的问题分析，可知问题一的目的是要求建立一个可以通用所有层的中心位置求解方法，即建立一个求解一层塔的中心位置的数学模型，通过上文分析可得，求解分为三大步骤：
一、构造各层层底几何中心的散点图；二、找出各层散点的最优点即为各层层底的几何中心坐标；三、连接各层层底的几何中心点，每两层层底之间连线的中点即是该层的几何中心点。
5.1.1构造各层层底几何中心散点图
结合上文问题分析的描述，可知附件提供的数据是各层层底的观测数据，但观测数据所连起来的图形的俯视图接近正八边形，因此对于各层层底的几何中心点的平面坐标可依据正八边形的几何关系计算，具体如下图：
图4-1：
二维平面八边形几何中心点关系图
图中有两个三角形，即通过选取六个点分别构造的两个三角形，在点1、2、3构造的三角形中，分别有两条虚线经过该三角形，两条虚线分别是边1,3与边2,3的中垂线，两条垂线相交于一点，该点即是三角形的外心，同时即是八边形的几何中心。
但在实际情况中，由于给出的点并不分布在同一水平面上，所以通过此法计算出来的点并非所有的两两中垂线都相交于同一点，所形成的是一个由
个点构成的散点图，构造该散点图的MATLAB代码详见附录，以86年第一层为例，散点图如下：
图4-2：
86年第一层层底散点图
5.1.2计算各层层底的几何中心
对于各层层底，必然存在一点，使得这一点到其它散点的距离之和最小，则认为该点即是其所属层底的几何中心点。
这一问题实际上是一个最优化点的搜索问题，针对这一问题的求解，有两种方法：
一、利用LINGO软件导入数据后，以求得到其它散点的距离之和最小的点构造目标函数，从而得到几何中心点；二、使用MATLAB中的直接搜索工具，寻找以到其它散点的距离之和为函数的最小值，从而得出几何中心。
方法一LINGO软件求解
以86年第一层层底数据为例，构造其目标函数如下：
（2）
约束条件为：
式中，min表示求各个散点到中心点的距离之和的最小值，d表示各个散点到中心点的距离，i表示第几个散点，Ui1表示第i个点的x坐标，Ui2表示第i个点的y坐标，Ui3表示第i个点的z坐标，x、y、z分别表示中心点的x、y、z坐标。
程序代码详见附录，代入第一层层底数据后，程序输出结果如下：
图5-1：
第1层层底中心坐标LINGO输出图
输出结果中A
（1）、A
（2）、A（3）分别表示中心点的X、Y、Z坐标，得中心点坐标为（566.6647,522.7108,1.787362），该点在散点图中的位置以“*”标识，如下图：
图5-2：
第1层层底几何中心点与其他散点空间图
使用方法二利用MATLAB工具箱求解[3]验证LINGO方法
MATLAB中包含一个专门设计的遗传算法与直接搜索工具箱（GeneticAlforithnmandDirectSearchToolbox），这个工具箱集成有遗传算法工具和直接搜索工具。
直接搜索（DirectSearch）是一种求解优化问题的方法，它不要求任何目标函数梯度的信息。
与使用梯度或高阶倒数信息来搜索优化点的较传统的优化算法相反，直接搜索算法搜索当前周围的一系列点，寻找目标函数值低于当前点值的那些点。
我们可以使用直接搜索算法来求解那些目标函数不可微、甚至不连续的问题。
直接搜索工具实现一类特殊的直接搜索算法，称为模式搜索（PatteernSearch）算法。
模式搜索算法确定一个点的序列，这个点序列呈现越来越接近理想点的趋势。
在每一步，该算法搜索在当前周围的一系列点，称为网格（Mesh）。
当前点是指算法在前一步计算出来的点。
该算法通过把当前点与一个称为模式（Pattern）的固定向量集的标量倍数相加来构成网格。
如果算法在网格中找到一个新点，且在该点比在当前点使目标函数得到改善，则该算法在下一步就将新点作为当前点。
在使用工具箱之前，需要构造目标函数my_fum，目标函数所求的是各散点到一位置点的距离之和，构造的目标函数的代码详见附录。
工具箱的具体操作步骤如下：
（1）健入psearchtool，打开模式搜索工具。
（2）在模式搜索工具的“Objectivefunction”文本框输入@my_fum
（3）在“Startpoint”文本框输入[123]
图6-1：
MATLAB模式搜索工具操作图
（4）勾选“Options”中“Plotfunctions”选项卡里的“Bestfunctionvalue”与“Meshsize”，单击“Start”按钮，运行模式搜索，得到结果如下：
图6-2：
模式搜索工具迭代过程图
图6-3：
模式搜索工具结果输出图
两种方法的比较
LINGO方法与MATLAB方法的结果比较如下表：
表2-1：
LINGO与MATLAB方法结果比较表
方法
X坐标
Y坐标
Z坐标
LINGO
566.6647
522.7108
1.787362
MATLAB
566.665
522.712
1.787
通过上表可以看出，两种方法计算的结果十分接近，相互证明了方法的准确性，验证了模型的合理性。
由于两种方法所计算的结果相差不大，在实际计算中为了操作简便采用LIN
升级会员 