2016全国大学生数学建模大赛国家一等奖论文A题.pdf

基于反演算法的系泊系统设计探究 摘要 本文探究的是近浅海观测网的系泊系统设计问题。通过平衡条件， 对系泊系统 的整体和各局部进行受力分析，得到系统各要素之间的递推公式。由此，本文采用 反演、 多目标优化等方法分析和解决本问题。 对于问题一，首先对题目中所给出的传输节点示意图建立直角坐标系，进行 受力分析，得出系统各部分的递推公式，从而根据递推-验证的反演算法建立模 型。再将未知量吃水深度设为已知量，以题目所给水深18为约束条件，通过 遍历算法运算出结果。接着对异常数据进行作零处理，分析锚链与海床接触的情 况，推导得出钢管倾角、锚链形状等。最后通过系统是否受力平衡对结果进行检 验。结果显示风速为12/时吃水深度为0.6870，第1-4根钢管的倾斜角(度)依次 为1.120、 1.127、 1.135、 1.142， 钢桶倾斜角(度)为1.150， 游动区域为半径14.34的旋 转面， 24/ 时吃水深度为0.7015/， 第1-4根钢管的倾斜角(度)依次为4.258、 4.285、 4.313、 4.341， 钢桶倾斜角(度)为4.370， 游动区域为半径17.46的旋转面。 对于问题二， 首先考虑第一个子问题， 将风速36/直接代入问题一的模型中， 得出此条件下的吃水深度为0.723， 各钢管倾斜角度(度)依次为8.960、 9.014、 9.068 、 9.123， 钢桶倾斜角(度)为9.179， 锚链链接处的切线方向与海床的夹角(度)为18.414， 游动区域半径为18.80。 发现此条件下， 水声通讯系统设备的工作效果较差,且锚被 拖行。随后考虑第二个子问题，对重物球重量进行调整，以钢桶倾斜角最小和吃水 深度最小为目标函数， 以水深18、 系统各部分递推关系式和钢桶倾斜角小于5度为 约束条件，以重物球重量、吃水深度为决策变量，建立多目标优化模型，再将其转 化为单目标优化。之后，通过分层递进法求出达到要求角度的重物球重力临界值， 求得重物球重力最小为2011。最后对可能存在的缺失值通过代入法进行了添补。 对于问题三，考虑到安全因素，将系统设定为处于较极端环境下，即水深、风 速、 海水速度分别设为16， 36/， 1.5/。在风速与水流速度同向的情况下， 采 用控制变量法对钢桶倾斜角和锚链末端与锚链接处的切线方向与海床的夹角随风 速、海水深度、海水速度影响的变化规律进行分析。从而可以确定使角度满足题意 的变量临界值。 然后， 以钢桶倾斜角最小和吃水深度最小为目标函数， 以水深16、 系统各部分递推关系式和钢桶倾斜角小于5度为约束条件，以重物球重量、吃水深 度和锚链长度为决策变量， 以风速36/、 海水速度1.5/ 为初始条件， 得出不同 锚链类型和长度的最优解。最后， 选取4组可行解， 并对所得系泊系统进行分析。 最后， 对模型作出了客观的评价及推广。 关键词：反演算法;多目标优化;遍历算法;递推公式;系泊系统。 1 一、 问题的重述 1.1 问题的背景 为满足海上资源开发利用的需要， 海上作业平台、 浮游码头及海上浮桥的应用 日渐增加， 这些海上漂浮结构物的海上定位需要锚泊系统的约束来实现。 锚泊系统在在海洋海岸工程中的形式多样，应用十分广泛。其中，近浅海观测 网的传输节点由浮标系统、 系泊系统和水声通讯系统组成， 浮筒式单点系泊系的设 计与应用对海上科学发展有重要意义。 1.2 问题的提出 已知某近浅海传输节点(如图1所示)， 将浮标视作底面直径2为、 高为2、 质量 为1000的圆柱体， 锚的质量为600， 钢管共4节， 每节长度为1， 直径为50， 每节钢管的质量为10。水声通讯系统安装在一个长为1、外径为30的密封圆 柱形钢桶内， 设备和钢桶总质量为100。 图 1: 传输节点示意图 锚链与锚的链接处的切线方向与海床的夹角不超过16度， 否则锚会被拖动。钢 桶竖直时， 水声通讯的效果最佳。 为了控制钢桶的倾斜角度， 链接处可悬挂重物球。 要求通过数学建模来完成以下任务： (1) 某传输节点选用22.05的II型锚链和1200的重物球， 布于深18m、 海床平 坦、 密度为1.025103/3的海域。 若海水静止， 计算海面风速为12/和24/ 时钢桶和各节钢管的倾斜角度、 锚链形状、 浮标的吃水深度和游动区域。 (2) 在问题1的假设下，计算海面风速为36/时钢桶和各节钢管的倾斜角度、 锚链形状和浮标的游动区域。 调节重物球的质量， 使钢桶的倾斜角不超过5度， 锚链 在锚点与海床的夹角不超过16度。 (3) 因潮汐等影响，布放海域的水深实测在16 20之间。布放点的海水速 度可达1.5/、 风速可达36/。 考虑风力、 水流力和水深， 请设计系泊系统， 分析 不同情况下钢桶、 钢管的倾斜角度、 锚链形状、 浮标的吃水深度和游动区域。 2 二、 问题分析 2.1 问题一的分析 问题一要求计算不同海面风速时钢桶和各节钢管的倾斜角度、 锚链形状、 浮标 的吃水深度和游动区域。要解决此问题， 首先对题目中所给出的传输节点示意图建 立直角坐标系， 对传输节点进行先整体后局部的受力分析， 得出系泊系统各部分的 递推公式， 建立递推-校正的反演模型， 得到相关相关方程组。 为求解此方程组，将题目所给的水深18作为约束条件，给吃水深度赋以初始 值，遍历所有可能的吃水深度值，记录对应的水下物体竖直方向的高度总和，取 与18 最接近时所对应的吃水深度为最优解， 推导得相应的倾斜角、 游动区域等。 2.2 问题二的分析 问题二的第一个子问题要求在问题一的基础上， 将风速设定为36/。 因此， 先 将其直接代入问题一的模型中，求得此条件下的吃水深度、各钢管倾斜角度、锚链 形状等。此时， 发现钢桶倾斜角和锚与锚链链接处的切线方向与海床的夹角均不符 合题目要求， 即水声通讯系统设备的工作效果较差,且锚被拖行， 因此无法求解游动 区域。 之后，对于第二个子问题，考虑调整重物球的重量，以使倾斜角度满足题目要 求。首先分析重物球对各个角度的影响，发现重物球越重角度越小。据此设立重物 球重力范围，在满足系统各部分递推关系式和钢桶与竖直方向夹角小于5度的条件 下，将多目标优化模型转化为单目标优化模型。然后遍历所有吃水深度的可取值， 寻找可行的情况， 再通过分层递进法求出恰好达到要求角度的重物球重力临界值。 2.3 问题三的分析 问题三要求考虑海水速度，并给出不同情况下钢桶、钢管的倾斜角度、锚链形 状、 浮标的吃水深度和游动区域。 以风速方向和水流速度方向相同时为例， 为得出各变量与角度的关系， 考虑采 用控制变量法对风速、 海水深度、 海水速度的变化分别对钢桶的倾斜角度和锚链末 端与锚链接处的切线方向与海床的夹角的影响进行分析， 再进一步确定使角度满足 题意的变量临界值。以水深16、系统各部分递推关系式和钢桶与竖直方向夹角小 于5为约束条件， 将多目标优化转化为单目标优化。 通过调节决策变量中锚链的型 号和长度，设定重物球重力的范围，对所有吃水深度的可取值进行遍历，判断对应 的钢管倾斜角度、海水深度是否小于临界值，小于则取其作为可行情况，否则重新 设定锚链的型号和长度，继续运算。最终设计出几种可行的系泊系统方案，并对其 进行分析。再分析在风速与水流速度不在同一方向时的情况， 对该模型进行优化。 3 三、 模型假设 1. 假设锚链的浮力忽略不计； 2. 假设浮标不倾斜； 3. 假设海床水平； 4. 假设水流方向与海床平行； 5. 假设重力加速度为10/。 四、 符号说明 符号意义 吃水深度 浮标所受浮力 浮标的质量 浮标的重力 浮标直径 浮标高度 锚所受的海床弹力 锚的重力 风力 钢管所受浮力 钢桶及以上各物的浮力和 锚链的浮力和 钢桶及以上各物的浮力和 锚链的重力和 锚链长度 链环长度 重球的重力 * 其余符号在文中详细说明 4 五、 模型的建立与求解 5.1 问题一模型的建立与求解 5.1.1 问题一模型的建立 (1) 模型的准备 假设吃水深度已知， 并令其为， 浮标质量为， 浮标直径为， 浮标高度为， 海水密度为。 浮标重力 = (1) 浮标所受浮力 = ( 2 )2 (2) 结合近海风荷载公式， 得到风力 = 0.6252= 0.625( )2(3) (2) 钢管和钢桶的倾斜角度模型 浮标的受力分析 图 2: 受力分析图 其中， 1为第一根钢管对浮标的拉力。 = + 11 = 11 (4) 其中, = 0.625( )2 = 5 根据以上公式可推导出1， 且其与第一根钢管与水平方向的夹角为两平行线间 的内错角， 所以相等。 钢管和钢桶的受力分析 图 3: 受力分析图 其中， 2为第一根钢管对第二根钢管的拉力。 对受力分析图进行进一步地分析， 建立如下函数关系式： 1= = ( 2 )2 1= 1 (5) 此外， 可得到： 2= + 1 1 1= 1 (6) 根据反三角关系， 得： 2= 2(7) 由此可推导得 +1= + =1 =1 , = 1,2,3,4(8) 其中， 当=1， 2， 3时， 为第+1根钢管与水平方向的夹角； 当=4 时， 为 钢桶与水平方向的夹角。 锚的受力分析 6 图 4: 受力分析图 其中， 为摩擦力， 为锚链对锚水平方向的拉力。 图 5: 受力分析图 其中， 为钢桶和重物球及以上各物体的浮力和， 为钢桶和重物球及以上 各物体的重力和，为锚链对钢桶的拉力，即锚的拉力、锚链的重力和浮力的合 力。 其中， 设重力球的质量为， 密度为， 重力球的浮力计算如下： = (9) = (10) 式中取钢的密度， 查阅文献的为7.85/3。 = = (11) = ( 1) , = 1,2, ,210(12) 7 式中为锚链长度， 为链环长度， 为单节链环的重力。 其中， = + 4+ + = + 4+ + = (13) 由此， 可计算出和 (3) 游动区域的计算 定义： 在海平面上半径为的旋转面。 = 5 =1 + 210 =1 其中， 为第个装置(包括各钢管、 钢桶)的长度， 为第个链环的长度。 吃水深度h成立的判断条件： = 5 =1 + 210 =1 = 18 其中， 为海域水深。 5.1.2 问题一模型的求解 Step1: 遍历求解 令吃水深度的初始值为0.1， 以0.0005为单位逐步增加至2。 ( 浮标高度为2， 完全浸没时吃水深度则为2 )，记录对应的数据，选取水下物体竖直方向高度和 与海域水深最接近的组别， 进一步进行计算， 结果如下表所示(具体程序见附录)： 表 1: 不同风速的相关结果表 以风速24/的情况为例， 绘制游动区域图： 8 图 6: 游动区域旋转面图 图表分析 由表可见， 风速为24/时的吃水深度大于12/时的， 且游动区域 半径也较12/时的大， 切合人们实际生活经验。 Step2: 异常数据的处理 在上述操作中， 发现部分负值， 不符合预期情况， 如表所示： 表 2: 异常数据表 经分析，认为异常数据为部分锚链接触海床而产生的，因此，下面将锚链接触 海床情况进行求解。 = ( 1) , = 1,2, ,210 当逐渐增大，锚链受到的竖直向下方向的合力与支持力之差先逐渐接近于0， 再等于0， 直至小于0。当合力小于0时， 锚链以海床接触， 此时海床提供向上的支持 力，其大小与相等。因此可将小于0 的值都作零处理，故锚链接触海床时， 为0。 Step3: 优化绘图 9 在遍历求解的基础上， 考虑锚链是否存在接触海床的情况， 通过将异常数据作 零处理的方式， 绘制下图(具体程序见附录)： 图 7: 不同风速下的锚链形状图 图表分析 由图可见， 风速为12/时锚链与海床接触明显， 风速为24/时仅 有少部分锚链与海床接触。 根据程序计算结果得到， 12/时有64节链环与海床接触， 24/时有1节链环 与海床接触。这与常识相符， 且锚链弧度符合客观实际。 Step4: 受力平衡检验 + + = = 210 其中，为浮标及水下各物的重力总和，水下物体的浮力总和为已知 量， 为锚链对锚竖直方向的拉力，浮标的浮力可表示为关于吃水深度的关系 式。 若计算结果满足受力平衡条件， 则可证明模型的正确性。 5.2 问题二模型的建立与求解 问题二分为两个子问题， 分别对两个子问题进行模型的建立与求解， 并对模型 进行了检验。 5.2.1 第一问模型的建立 由于第一问建立在问题一的假设下， 且仅改变风力， 所以直接应用问题一所建 模型， 详见”5.1.1 问题一模型的建立” 。 5.2.2 第一问模型的求解 由于第一问建立在问题一模型的基础上，将风速= 36代入，得到结果如下 10 表所示(具体程序见附录)： 表 3: 风速为36m/s时的相关结果表 锚链形状如下图所示(具体程序见附录)： 图 8: 风速为36m/s时的锚链形状图 发现锚链末端与锚的链接处的切线方向与海床的夹角超过16度， 锚被拖动， 故 无法计算此时浮标的游动区域。 5.2.3 第二问模型的建立 (1) 目标函数的建立 5(14) (15) (2) 约束条件的建立 11 . = 18 85 5 90 5= (+ 4 =1 4 =1 )/ 210= (+ 4 =1 4 =1 )/ (3) 决策变量的建立 , 5.2.4 第二问模型的求解 Step1: 绘制流程图 图 9: 流程图 Step2: 分层递进求解 为寻找使钢桶的倾斜角度1不超过5度， 锚链在锚点与海床的夹角2不超过16度 时的重物球的最小重力，需选择一个重物球重力取值的较大范围，初步将其设定 为12000,30000， 并取200作为粗略单位进行试运算， 得到此范围内满足该条件的 重物球重力如表2所示(具体程序见附录)： 表 4: 重物球重力部分取值表 12 图表分析 重物球最小重力粗略估计为20200， 且大于该值的任意重物球皆满 足题目要求。 为进一步确定重物球最小重力，将取值范围缩小为20000,20200，得到此范围 内满足该条件的重物球重力为： = 20110 结果解释 该值为满足题目要求下的重物球最小重力，重力大于等于该值的重 物球为可行值。 Step3: 缺失数据的处理 实际运行过程中发现在某一重物球重力范围内， 由于阈值选取不够精确， 存在 部分可行值被排除的可能性。因此， 通过代入法进行了添补。 5.3 问题三模型的建立与求解 5.3.1 问题三模型的建立 问题三可能出现风向与海水流向相同、 风向与海水流向不同两种情况， 考虑将 其分别讨论， 并进行以下分析。 (1) 变量分析 通过控制变量法绘制如下三种因素与角度的关系图(具体程序见附录)： 图 10: 角度与风速关系图 图表分析 由图可知， 钢桶的倾斜角度随风速的增加而增加， 逐渐接近5 度， 且 变化速率随之增大。此外，在风速小于12/时，锚链末端与锚链接处的切线方向 与海床的夹角始终为0， 当风速大于等于12/时， 夹角随风速的增加而增加， 且变 化速率逐渐减小。 13 图 11: 角度与水深关系图 图表分析 由图可知， 钢桶的倾斜角度随水深的增加而缓慢减小， 逐渐接近0度。 此外， 锚链末端与锚链接处的切线方向与海床的夹角随水深的增加而增加。当水深 为16 时， 钢桶的倾斜角度恰好满足5度， 且将呈下降趋势， 因此水深在此基础上继 续加大， 角度依然满足题目要求。 图 12: 角度与海水速度关系图 图表分析 由图可知，钢桶的倾斜角度随海水速度的改变有轻微增大，且锚链 末端与锚链接处的切线方向与海床的夹角随海水速度的增加而呈现不均匀的增加。 (2) 多目标优化模型的建立 初始条件 = 5 =1 + =1 = 16 = 36 = 1.5 (16) 目标函数 5 14 约束条件 . = 16 85 5 90 0 16 5= (+ 4 =1 4 =1 )/(+ 4 ) = (+ 4 =1 4 =1 )/(+ 4 + ) = 3742 = 其中，为锚链末端与锚的链接处的切线方向与海床的夹角，为海水在 水平方向上的作用力,为钢桶和重力球在水平方向上受到的力。 决策变量, , 5.3.2 问题三模型的求解 (1) 临界条件的选取 通过上述图表分析， 得出以下结论： 1.临界风速为36/ 2.临界海水深度为16/ 3.临界海水速度为1.5/ (2) 求解优化模型 图 13: 流程图 15 通过遍历求解， 得出了不同情况下的大量可行解， 筛选四组系泊系统的设计方 案，求解临界条件下相应的吃水深度、各倾斜角度，绘制如下表格(具体程序见附 录)： 表 5: 部分组的相关结果表 表 6: 钢管与竖直方向的夹角结果表 绘制各编号对应的锚链形状图(具体程序见附录)： 图 14: 编号1的锚链形状图图 15: 编号2的锚链形状图 图 16: 编号3的锚链形状图 图 17: 编号4的锚链形状图 5.3.3 问题三模型的优化 考虑到风向有可能不同， 进行如下受力分析： 16 图 18: 受力分析图 由上图可知，风向不同条件下近海风荷载为风力在水平方向上的分力与海水 作用力的合力，即= + ,小于同向时的合力。由角度与风速的关 系(图11)可知， 合力越小， 速度越小， 倾角越小， 因此风向相同时系泊系统的设计方 案满足风向不同时的情况。 六、 模型评价与推广 6.1 模型的优点 (1) 在问题一中对模型结果采用了受力平衡检验法，从而验证了所得分析结果 的合理性和可靠性。 (2) 在问题二中采用了分层递进算法，在不影响准确度的情况下有效地缩短了 程序的运行时间， 提高了求解效率。 (3) 在问题三中提出了临界值概念，考虑极端条件下系泊系统的设计方案，简 化思路， 使结果适用于任何情况。 6.2 模型的缺点 (1) 重物球材料未知， 密度为人为假设， 可能与实际情况存在一定误差。 (2) 题目所给锚链资料有限， 无法计算锚链所受浮力， 对模型精确度存在影响。 6.3 模型的推广 (1) 此文所建立的浮筒式单点系泊系统的相应模型可进一步分析，建立多点系 泊系统的数学模型。 (2) 此文是建立在系泊系统布于近浅海域的基础上， 可结合实际海上作业情况， 将范围扩展至深海。 参考文献 1姜启源,谢金星,叶俊.数学模型M.北京:高等教育出版社,2003. 17 2周凯,宋军全,邬学军.数学建模竞赛入门与提高M.浙江:浙江大学出版社, 2003. 3潘斌,高捷,陈小红,陈家鼎.浮标系泊系统静力计算J.重庆交通学院学报, 1997,16(1):68- 73. 4于定勇.水下锚泊系统计算J.青岛海洋大学学报, 1995:100- 105. 18 附录 一.问题一 matlab 程序 (1)钢桶和各节钢管的倾斜角度、浮标的吃水深度和游动区域的求解： shang=; i=0; r=12;%风速 ml=; H=; for h=0.1:0.0005:2%吃水深度 i=i+1; shang(i,1)=atand(1025*3.14*h*10-10000)/(0.625*r*r*2*(2-h);%第 1 根钢管与水平方向 夹角 shang(i,2)=atand(1025*3.14*h*10-10000-79.884375)/(0.625*r*r*2*(2-h); shang(i,3)=atand(1025*3.14*h*10-10000-79.884375*2)/(0.625*r*r*2*(2-h); shang(i,4)=atand(1025*3.14*h*10-10000-79.884375*3)/(0.625*r*r*2*(2-h); shang(i,5)=atand(1025*3.14*h*10-10000-79.884375*4)/(0.625*r*r*2*(2-h);% 钢 桶 与 水 平方向夹角 x1=sind(shang(i,1)+sind(shang(i,2)+sind(shang(i,3)+sind(shang(i,4)+sind(shang(i,5); y1=cosd(shang(i,1)+cosd(shang(i,2)+cosd(shang(i,3)+cosd(shang(i,4)+cosd(shang(i,5); k=0; for j=1:210 k=k+1; ml(k,i)=atand(1025*3.14*h*10-10000-79.884375*4-275.8375-12000+1566.88-7. 35*(k-1)/(0.625*r*r*2*(2-h);%第 j 节链环与水平方向夹角 if ml(k,i)=0 ml(k,i)=0; end end x2=sum(0.105*sind(ml(:,i); y2=sum(0.105*cosd(ml(:,i); H(i,1)=x1+x2+h;%系统水下部分竖直方向高度 H(i,2)=y1+y2;%浮标游动半径 end (2)锚链形状图像的求解： r=12; k=0; ml=; g=; h=0.687; for j=1:210 k=k+1; ml(k,1)=atand(1025*3.14*h*10-10000-79.884375*4-275.8375-12000+1566.88-7. 35*(k-1)/(0.625*r*r*2*(2-h); if ml(k,1)=0 ml(k,1)=0; end end mll=flipud(ml); g(1,1)=cosd(mll(1,1)*0.105; g(1,2)=sind(mll(1,1)*0.105; for i=1:209 g(i+1,1)=cosd(mll(i+1,1)*0.105+g(i,1); g(i+1,2)=sind(mll(i+1,1)*0.105+g(i,2); end plot(g(:,1),g(:,2) hold on r=24; k=0; ml=; g=; h=0.7015; for j=1:210 k=k+1; ml(k,1)=atand(1025*3.14*h*10-10000-79.884375*4-275.8375-12000+1566.88-7. 35*(k-1)/(0.625*r*r*2*(2-h); if ml(k,1)=0 ml(k,1)=0; end end mll=flipud(ml); g(1,1)=cosd(mll(1,1)*0.105; g(1,2)=sind(mll(1,1)*0.105; for i=1:209 g(i+1,1)=cosd(mll(i+1,1)*0.105+g(i,1); g(i+1,2)=sind(mll(i+1,1)*0.105+g(i,2); end plot(g(:,1),g(:,2) 二.问题二 matlab 程序 (1)问题一假设下钢桶和各节钢管的倾斜角度、浮标的吃水深度和游动区域的求解： shang=; i=0; r=36; ml=; H=; for h=0.1:0.0005:2 i=i+1; shang(i,1)=atand(1025*3.14*h*10-10000)/(0.625*r*r*2*(2-h); shang(i,2)=atand(1025*3.14*h*10-10000-79.884375)/(0.625*r*r*2*(2-h); shang(i,3)=atand(1025*3.14*h*10-10000-79.884375*2)/(0.625*r*r*2*(2-h); shang(i,4)=atand(1025*3.14*h*10-10000-79.884375*3)/(0.625*r*r*2*(2-h); shang(i,5)=atand(1025*3.14*h*10-10000-79.884375*4)/(0.625*r*r*2*(2-h); x1=sind(shang(i,1)+sind(shang(i,2)+sind(shang(i,3)+sind(shang(i,4)+sind(shang(i,5); y1=cosd(shang(i,1)+cosd(shang(i,2)+cosd(shang(i,3)+cosd(shang(i,4)+cosd(shang(i,5); k=0; for j=1:210 k=k+1; ml(k,i)=atand(1025*3.14*h*10-10000-79.884375*4-275.8375-12000+1566.88-7. 35*(k-1)/(0.625*r*r*2*(2-h); if ml(k,i)=0 ml(k,i)=0; end end x2=sum(0.105*sind(ml(:,i); y2=sum(0.105*cosd(ml(:,i); H(i,1)=x1+x2+h; H(i,2)=y1+y2; end (2)问题一假设下锚链形状图像的求解： r=36; k=0; ml=; g=; h=0.7655; for j=1:210 k=k+1; ml(k,1)=atand(1025*3.14*h*10-10000-79.884375*4-275.8375-12000+1566.88-7.35*(k-1)/(0.62 5*r*r*2*(2-h); if ml(k,1)=0 ml(k,1)=0; end end mll=flipud(ml); g(1,1)=cosd(mll(1,1)*0.105; g(1,2)=sind(mll(1,1)*0.105; for i=1:209 g(i+1,1)=cosd(mll(i+1,1)*0.105+g(i,1); g(i+1,2)=sind(mll(i+1,1)*0.105+g(i,2); end plot(g(:,1),g(:,2) （3）重物球的调整： shang=; r=36; ml=; H=; daan=; m=0; for g=20000:10:20200 i=0; for h=0.1:0.005:2 i=i+1; shang(i,1)=atand(1025*3.14*h*10-10000)/(0.625*r*r*2*(2-h); shang(i,2)=atand(1025*3.14*h*10-10000-79.884375)/(0.625*r*r*2*(2-h); shang(i,3)=atand(1025*3.14*h*10-10000-79.884375*2)/(0.625*r*r*2*(2-h); shang(i,4)=atand(1025*3.14*h*10-10000-79.884375*3)/(0.625*r*r*2*(2-h); shang(i,5)=atand(1025*3.14*h*10-10000-79.884375*4)/(0.625*r*r*2*(2-h); x1=sind(shang(i,1)+sind(shang(i,2)+sind(shang(i,3)+sind(shang(i,4)+sind(shang(i,5); k=0; for j=1:210 k=k+1; ml(k,i)=atand(1025*3.14*h*10-10000-79.884375*4-275.8375-g+g/7850*1025-7.35*(k-1)/(0.62 5*r*r*2*(2-h); if ml(k,i)=0 ml(k,i)=0; end end x2=sum(0.105*sind(ml(:,i); H(i)=x1+x2+h; if abs(H(i)-18)85 daan(m,1)=g; end end end 三.问题三 matlab 程序 (1)风速、海水速度、水深临界值的确定（以风速为例） shang=; i=0; daan=; ml=; H=; g=20110; v=1.5; q=0; for r=0:2:36 q=q+1; for h=0.1:0.0005:2 i=i+1; shang(i,1)=atand(1025*3.14*h*10-10000)/(0.625*r*r*2*(2-h)+2*h*v*v*374); a=0.05*h*v*v*374*sind(shang(i,1); shang(i,2)=atand(1025*3.14*h*10-10000-79.884375)/(0.625*r*r*2*(2-h)+a); b=0.05*v*v*374*sind(shang(i,2); shang(i,3)=atand(1025*3.14*h*10-10000-79.884375*2)/(0.625*r*r*2*(2-h)+a+b); c=0.05*v*v*374*sind(shang(i,3); shang(i,4)=atand(1025*3.14*h*10-10000-79.884375*3)/(0.625*r*r*2*(2-h)+a+b+c); d=0.05*v*v*374*sind(shang(i,4); shang(i,5)=atand(1025*3.14*h*10-10000-79.884375*4)/(0.625*r*r*2*(2-h)+a+b+c+d); e=0.3*v*v*374*sind(shang(i,5); x1=sind(shang(i,1)+sind(shang(i,2)+sind(shang(i,3)+sind(shang(i,4)+sind(shang(i,5); y1=cosd(shang(i,1)+cosd(shang(i,2)+cosd(shang(i,3)+cosd(shang(i,4)+cosd(shang(i,5); k=0; for j=1:210 k=k+1; f=0.42*v*v*374; ml(k,i)=atand(1025*3.14*h*10-10000-79.884375*4-275.8375-g+g/7850*1025-7.35*(k-1)/(0.62 5*r*r*2*(2-h)+a+b+c+d+e+f); if ml(k,i)=0 ml(k,i)=0; end (2)重物球重力和浮标吃水深度的求解 shang=; daan=; ml=; H=; v=1.5; r=36; q=0; for g=12000:500:24000 i=0; for h=0.7:0.001:2 i=i+1; shang(i,1)=atand(1025*3.14*h*10-10000)/(0.625*r*r*2*(2-h)+2*h*v*v*374); a=0.05*h*v*v*374*sind(shang(i,1); shang(i,2)=atand(1025*3.14*h*10-10000-79.884375)/(0.625*r*r*2*(2-h)+a); b=0.05*v*v*374*sind(shang(i,2); shang(i,3)=atand(1025*3.14*h*10-10000-79.884375*2)/(0.625*r*r*2*(2-h)+a+b); c=0.05*v*v*374*sind(shang(i,3); shang(i,4)=atand(1025*3.14*h*10-10000-79.884375*3)/(0.625*r*r*2*(2-h)+a+b+c); d=0.05*v*v*374*sind(shang(i,4); shang(i,5)=atand(1025*3.14*h*10-10000-79.884375*4)/(0.625*r*r*2*(2-h)+a+b+c+d); e=0.3*v*v*374*sind(shang(i,5); x1=sind(shang(i,1)+sind(shang(i,2)+sind(shang(i,3)+sind(shang(i,4)+sind(shang(i,5); y1=cosd(shang(i,1)+cosd(shang(i,2)+cosd(shang(i,3)+cosd(shang(i,4)+cosd(shang(i,5); k=0; for j=1:160 k=k+1; f=0.42*v*v*374; ml(k,i)=atand(1025*3.14*h*10-10000-79.884375*4-275.8375-g+g/7850*1025-15*(k-1 )/