系泊系统设计

来源：步遥情感网

2016高教社杯全国大学生数学建模竞赛

承诺书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.

我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是（从A/B/C/D中选择一项填写）：

我们的参赛报名号为（如果赛区设置报名号的话）：

所属学校（请填写完整的全名）：参赛队员 (打印并签名) ：1. 2. 3. 指导教师或指导教师组负责人 (打印并签名)：

日期：年月日

赛区评阅编号（由赛区组委会评阅前进行编号）：

2016高教社杯全国大学生数学建模竞赛

编号专用页

评阅人评分备注赛区评阅编号（由赛区组委会评阅前进行编号）：

赛区评阅记录（可供赛区评阅时使用）：

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）：

系泊系统的设计

摘要

系泊系统的设计问题是一个具有悠久历史的问题，本文根据系泊系统的设计要求，建立优化模型，求解得到使浮标的吃水深度、游动距离以及钢桶的倾斜程度尽可能小的最优解。

针对问题一，本文首先对系泊系统各部分进行受力分析，得到了各部分关于受力平衡和力矩平衡的表达式，初步建立了系泊系统的力学模型，但是由于力学方程组的个数多且复杂，难以直接求解出答案，因此本文运用搜索算法对方程进行简化得到了问题一的结果，在风速为12m/s时浮标的吃水深度为0.735m，游动区域半径为14.184m，钢桶的倾斜程度为1.008°，钢管一到四的倾斜角度分别为0.977°，0.983°，0.988°以及0.994°，在风速为24m/s时浮标的吃水深度为0.751m，游动区域半径为16.919m，钢桶的倾斜程度为3.825°，钢管一到四的倾斜角度分别为3.712°，3.733°，3.755°以及3.776°。

针对问题二，本文首先计算了风速为36m/s时系泊系统各部分参数，发现所得结果不满足题设稳定条件，因此本文考虑在改变重物球重量的情况下寻找使浮标的吃水深度、游动距离以及钢桶的倾斜程度尽可能小的系统系泊最优设计参数，在分析了系泊系统的优化条件和约束条件后，本文成功建立了基于重物球质量调节的系泊系统优化设计模型，并运用变步长的搜索算法得到了问题二的结果，在风速为36m/s的情况下系泊系统的最优设计参数如下重物球的质量为5090kg,浮标的吃水深度为1.940m，游动区域半径为12.300m，钢桶的倾斜程度为0.109°，钢管一到四的倾斜角度分别为0.108°，0.109°，0.109°以及0.109°。

针对问题三，问题三是在问题二基础上的一个延伸，本文首先考虑水流力的影响，对系泊系统各部分重新进行受力分析，得到各部分关于受力平衡和力矩平衡的表达式，之后基于问题二的优化模型，建立了基于重物球质量变初始条件最优系泊系统模型，最后本文同样运用变步长的搜索方法，对本文自己设计的风速为24m/s，水流速度为1.5m/s，水深变化范围为16~20m的情况给出了相对最优解。

最后本文对本文所建立的模型进行了客观的评价，总结了模型的优缺点。

【关键词】刚体力学搜索算法最优设计多目标优化

一、问题重述

1.1 问题背景

系泊系统的设计问题是一个历史很悠久的问题，它的要求就是要确定锚链的型号、长度和重物球的质量，使得浮标的吃水深度和游动区域及钢桶的倾斜角度尽可能小，基于此要求我们依次解决了以下三个问题。 1.2 目标任务

问题一：某型传输节点选用II型电焊锚链22.05m，选用的重物球的质量为1200kg。现将该型传输节点布放在水深18m、海床平坦、海水密度为1.025×103kg/m3的海域。若海水静止，分别计算海面风速为12m/s和24m/s时钢桶和各节钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域。

问题二：在问题1的假设下，计算海面风速为36m/s时钢桶和各节钢管的倾斜角度、锚链形状和浮标的游动区域。请调节重物球的质量，使得钢桶的倾斜角度不超过5度，锚链在锚点与海床的夹角不超过16度。

问题三：由于潮汐等因素的影响，布放海域的实测水深介于16m~20m之间。布放点的海水速度最大可达到1.5m/s、风速最大可达到36m/s。请给出考虑风力、水流力和水深情况下的系泊系统设计，分析不同情况下钢桶、钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域。

二、问题分析

系泊系统的设计问题是一个具有悠久历史的问题，其主要要求为确定锚链的型号、长度和重物球的质量，使得浮标的吃水深度和游动区域及钢桶的倾斜角度尽可能小。

对于问题一，首先本文依次对浮标、钢桶、钢管、重物球、锚链和锚进行受力分析，之后根据受力分析，并考虑锚链是否沉底，得出浮标的游动范围的表达式，最后本文依据之前分析得出的力学平衡方程，使用MATLAB联立求解，得出钢桶和各节钢管的倾斜角度、锚链形状以及浮标的游动区域半径。

对于问题二，问题二是一个多目标优化问题，首先我们根据题设的要求逐步分析待优化的条件和约束条件，之后构造目标函数整体地表示系泊系统的优化程度，最后使用变步长的搜索算法得到了系泊系统相对最优的一组设计参数。

对于问题三，问题三是在问题二基础上的进一步延伸,它和问题二的区别主要在于，问题三考虑了水流力，因此我们要在要考虑水流力影响的情况下对系泊系统各部分重新进行受力分析，基于问题一和问题二的步骤逐步建立基于重物球质量变初始条件最优系泊系统模型，之后同样采用变步长的搜索算法得到了系泊系统相对最优的一组设计结果。不可避免地由于在第三问中有4个待搜索变量，计算量大大增加，因此得到结果的准确度可能会受到一定的影响。

三、模型假设

1.风力和水流力平行于海平面，且浮标不会倾斜。 2.锚链整体可视为一条可微的曲线，且质量均匀。 3.忽略锚链和重物球的浮力影响。

4.水流力对钢管和钢桶的作用面积近似为钢管和钢桶的纵截面积。

四、符号说明

符号 f F h G T  L  H s S 含义各物体所受浮力风海风荷载浮标的吃水深度重力物体所受到的拉力系统各部分所受的拉力与竖直方向的夹角钢管的长度钢管与竖直方向的夹角系统各部分垂直高度锚链长度系统各部分的水平长度单位牛顿牛顿米牛顿牛顿度米度米米米

五、模型的建立与求解

5.1 问题一模型的建立与求解 5.1.1 系泊系统各部分受力分析 5.1.1.1 浮标受力分析

首先我们考虑实际情况，假设浮标始终保持竖直状态不会倾斜，并且风始终和海平面平行，在此基础上我们进行浮标的受力分析。浮标一共受到4个力的作用，分别是风力Fb，浮力fi，钢管的拉力T1，浮标自身的重力Gb，如图一所示

图一浮标受力分析图

要使浮标受力平衡，我们应该让浮标的水平方向和竖直方向都保持平衡，我们可以得到以下表达式

fb=T1cos1GbFT1sin1其中根据阿基米德定律可知

fbgV

V是浮标浸入水面的体积易得为h，其中h是浮标的吃水深度。

风力由题目已知条件可得，v为风速，S为物体在风向法平面的投影面积。

F0.625Sv2

5.1.1.2 钢管受力分析

对于钢管的受力分析，由于4节钢管的受力分析类似，我们在此就分析其中一节钢管，其它刚管依次类推。我们假设4节钢管从上到下分别为钢管1，钢管2，钢管3，钢管4，对于第i节钢管，它所受到的右侧的拉力是Ti，它所受到的左侧的拉力为Ti-1，其中AB为钢管，如图二所示

图二钢管受力分析图

根据受力平衡和力矩平衡，我们可以得到以下表达式

Gg1T2cos2fg1T1cos1 T1sin1T2sin20.5TLsin()0.5TLsin()1g112g21其中L是钢管的长度，1是钢管与竖直方向的夹角。

依次类推我们可以依次得到钢管2，钢管3，钢管4的表达式。

钢管2的受力平衡和力矩平衡表达式为：

Gg2T3cos3fg2T2cos2 T2sin2T3sin30.5TLsin()0.5TLsin()2g223g32 钢管3的受力平衡和力矩平衡表达式为：

Gg3T4cos4fg3T3cos3 T3sin3T4sin40.5TLsin()0.5TLsin()3g334g43钢管4的受力平衡和力矩平衡表达式为：

Gg4T5cos5fg4T4cos4 T4sin4T5sin50.5TLsin()0.5TLsin()4g445g545.1.1.3 钢桶与重物球受力分析

如图三所示我们对钢桶和小球进行受力分析，图中DC就为钢桶，小球看成一个点处在D的位置。

图三钢桶和重物球受力分析图

如图所示，钢桶和小球一共受到5个力的作用，钢桶左侧受到的拉力T5,右侧受到的拉力T6，钢桶受到的浮力ft，以及它的重力Gt和小球对钢桶的拉力Tzwq，根据受力平衡分析,可以得到以下表达式：

GtTzwq+T6cos6ftT5cos5T5sin5T6sin

然后考虑钢桶与重物球的力矩平衡，由于钢桶的重力以及浮力作用在质心所以不予考虑，我们只要考虑钢桶与重物球左侧受到的拉力T5，右侧受到的拉力T6，以及小球对钢桶的拉力Tzwq，根据力矩平衡的公式可以得到：

0.5T5Ltsin(55)0.5TzwqLtcos(905)0.5T6Ltsin(65)

其中重物球对钢桶的拉力可以由图四得到

图4 重物球受力分析图

根据竖直方向的受力分析我们可以很简单的得到

TzwqfzwqGzwq

5.1.1.4 锚链的受力分析

由于锚链是无档普通链环，我们可以将其看作无弹性悬垂线，因此我们可以用微元的方法来分析其受力状态。

如图五所示，考虑ds这样一段任意的小弧长

图五锚链受力分析图

这样一段小弧长一共受到三个力的作用，分别是上段锚链对该小弧长的拉力T,下段锚链对该小弧长的拉力T+Dt，以及小弧长自身所受到的重力G,要使这一段小弧长受力平衡，就要使该小弧长在水平方向上和竖直方向上受力平衡，可以得到如下表达式：

TcosTdTcosd TsinTdTsind+dsg将该表达式展开可以得到

TcosTdTcoscosdsinsind TsinTdT(sincosdcossind)+dsg 因为ds是一段任意的小弧长，可以认为ds近似于趋向于0，当ds趋向于0时，d也趋向于0，因此可知cosd趋向于1，sind趋向于d，由此上式可以化简为

dTcosTdsindTdsin0 TdcosdTsindTdcosdsg0 又因为dTd是高阶无穷小，可以忽略不计，这样就可以化简得到锚链的最终表达式。

gcosddsT dTgsinds 并且由问题所给的已知条件可以得到该方程的初始条件为

06T0T6

5.1.1.5 锚的受力平衡分析

如图六所示，锚一共受到4个力的作用，分别是海底给锚的支持力，锚自身的重力，锚受到的海底给它的摩擦力，以及锚链给锚的拉力。

图6 锚受力分析图

由图可知要使锚受力平衡，应该要使锚的水平方向和竖直方向都保持平衡，可以得到如下表达式：

TmcosFfTmsinFNGm

5.1.2 系泊系统垂直高度分析

根据上文受力分析，我们可以简单得到系泊系统各部分的高度，在此我们定义，系泊系统的总高度为H，浮标的水中垂直高度为Hb，钢管的水中垂直高度为Hg，钢桶的水中垂直高度为Ht，锚链的水中垂直高度为Hl。

首先浮标的垂直高度就为浮标的浸水深度，因此

Hbh

其次钢管的水中垂直高度为4节钢管的长度总和，每节钢管长度是它们各自的长度乘以它们的倾斜程度，因此

Hgli14gcosi

然后钢管的水中垂直高度等于它的长度乘以它的倾斜程度

Htltcos5

最后由于锚链可以近似的看成无弹性的悬垂线，因此它的垂直高度可以用微元法积分得到

Hl622.50dssind

将该式拆开得到

Hl6sindscosdsd

022.5因此我们可以得到系泊系统的总高度H为

HHbHgHtHl

5.1.4 浮标游动距离分析

如图七所示首先我们定义浮标的游动距离为锚链到浮标中点的距离。

图7 浮标游动距离示意图

和系泊系统垂直高度分析类似,根据上文的受力分析我们也可以得到系泊系统各部分的水平距离,从而得到浮标的游动距离, 在此我们定义，浮标的游动距离为S，钢管的水平距离为Sg，钢桶的水平距离为St，锚链的水平距离为Sl。

用与系泊系统垂直高度类似的分析方法我们可以得到

Sglgsinii14Stltsin5Sl6

dcosd22.50因此浮标的游动距离就为这三部分的水平距离之和

SSgStSl

5.1.4 模型的求解算法 5.1.4.1 浮标部分求解

首先我们假设浮标的吃水深度h已知，由此进行分析，由上文对浮标的分析我们已经知道

fbgV

2 F0.625Sv

在h已知的条件下我们就可以算出F与fb，将这两个值代入下式

fb=T1cos1GbFT1sin1

就可以解出T1和1。 5.1.4.2 钢管部分求解

在已经得到T1和1后，我们分析钢管部分，根据上文钢管一的受力分析表达式

Gg1T2cos2fg1T1cos1 T1sin1T2sin20.5TLsin()0.5TLsin()1g111g218

三个方程三个未知数，我们就可以算出T2和2和1的大小，依次类推将T2和2代入钢管二的受力分析表达式

Gg2T3cos3fg2T2cos2 T2sin2T3sin30.5TLsin()0.5TLsin()2g222g32我们就可以算出T3和3，再将T3和3代入钢管三的受力分析表达式，再将求出的结果代入钢管四的受力分析表达式，就可以得到T4和4，T5和5。 5.1.4.3 钢桶部分求解

在已经得到T5和5的后我们分析钢桶部分，根据上文钢桶的受力分析表达式

GtTzwq+T6cos6ftT5cos5T5sin5T6sin60.5T5Ltsin(55)0.5TzwqLtcos(905)0.5T6Ltsin(65)

三个方程三个未知数，可以计算出T6和6和5。 5.1.4.4 锚链部分求解

在已知T6和6后，我们分析锚链部分，由锚链部分的表达式和初值条件

gcosddsT dTgsinds

06T0T6

我们就可以用matlab画出锚链的图形[2]并计算出锚链对锚的拉力的大小Tm。 5.1.4.5 锚部分求解

在已经知道锚链对锚的拉力的大小Tm后，并且对系统的整体分析可以知道锚受到的摩擦力等于风力的大小，由此根据锚的受力平衡表达式

TmcosFfTmsinFNGm

我们就可以算出的大小。

5.1.4.6 模型具体求解过程

我们已经知道在已经知道浮标的吃水深度以后，我们可以一步步推算出系泊系统的各个参数，因此我们考虑用搜索算法[2]对模型进行求解，具体求解过程如下 1．取一个最初的浮标吃水深度h,由它一步步推算得到H和。

2．判断得到的H和是否同时满足H与18的差值小于一个十分小的量（比如说这个量取0.0.1），并且的角度小于16度，如果满足这两个条件就输出结果，否则就让h增加一个h的长度重复步骤一，直到输出结果。

3．其中当浮标只受到重力和浮力的影响时，浮标的吃水深度最小，经计算得为0.31，因此我们把0.31作为浮标最初的吃水深度。 5.1.5 模型一问题求解结果展示

5.1.5.1风速为12 m/s时求解结果展示

根据上文的求解过程我们可以得到风速为12m/s时系泊系统的状态如下表格所示

表1 风速为12 m/s 时的系泊系统参数表

得到的锚链形状如图8所示

图8 风速为12 m/s 时的锚链形状图

根据图8和参数表我们可以知道在风速为12m/s时锚链沉底，即锚链和锚的夹角为0度。

5.1.5.2风速为24 m/s时求解结果展示

根据上文的求解过程我们可以得到风速为24m/s时系泊系统的状态如下表格所示

表2 风速为24 m/s 时的系泊系统参数表

得到的锚链形状如图9所示

图9 风速为24 m/s 时的锚链形状图

根据图9和参数表我们可以知道在风速为24m/s时锚链不沉底，锚链和锚的夹角为3.248度。

5.2 问题二模型的建立与求解

5.2.1 风速为36 m/s时系泊系统状态分析

当风速为36m/s时，我们可以根据问题一的分析方法分析系泊系统的状态，结果如下。

表3 风速为36 m/s 时的系泊系统参数表

锚链的形状如图10所示

图10 风速为36 m/s时锚链形状图

根据计算结果我们可以发现当风速为36m/s时，锚链与锚的夹角为21.363度，不满足题目所给的锚链与锚的夹角要小于16度的条件，同时钢桶的倾斜角度为7.975度不满足题目所给的钢桶倾斜程度小于5度的条件，因此我们考虑调整重物球的重量，使得钢桶的倾斜角度不超过5度，锚链在锚点与海床的夹角不超过16度。 5.2.2 基于重物球质量调节的系泊系统优化设计模型的建立

通过对风速为36m/s时的系泊系统分析我们已经知道当风速较大时，钢桶的倾斜程度和锚链与锚的夹角将不满足题目条件，因此我们要调整重物球的重量，来使得钢桶的倾斜角度不超过5度，锚链在锚点与海床的夹角不超过16度，但是在调整重物球重量的过程中，满足这两个条件的系统有很多个，哪个才是我们所要寻找的设计最优的系泊系统呢？为了寻找设计最优的系泊系统，我们考虑建立以下三个目标函数来确定浮标的吃水深度和游动区域及钢桶的倾斜角度最小的系泊系统。

5.2.2.1 目标函数分析 5.2.2.1.1 浮标吃水深度

要让系泊系统设计达到最优，该系统根据题设条件应该满足浮标吃水深度最小的条件，根据问题一的分析，即要满足

minh

5.2.2.1.2 浮标游动区域

要让系泊系统设计达到最优，该系统根据题设条件应该满足浮标游动区域最小的条件，根据问题一的分析，即要满足

minS

其中

SSgStSlSglgsinii14Stltsin5Sl6

22.50dcosd5.2.2.1.3 钢桶倾斜角度

要让系泊系统设计达到最优，该系统根据题设条件应该满足钢桶倾斜角度最小的条件，根据问题一的分析，即要满足

min5

其中5可以由钢桶的力矩平衡和受力平衡表达式求得。

5.2.2.1.4 目标函数的确定

在分析了浮标的吃水深度，浮标游动区域和钢桶倾斜程度后,我们可以确定总的目标函数M,令

M5hS

当M最小时即可达到相对最优的系泊系统设计，但是以上表达式没有统一量纲，

5和h对M的影响就可以忽略不计了，在此种情况下假如S相对于5和h特别大，

因此我们考虑对目标函数M进行优化，让三个变量各自除以它们在系统中可能达到的

最大值以统一三个变量的量纲[5]。

M55maxhhmaxSSmax

5.2.2.2约束条件分析

5.2.2.2.1系泊系统总高度分析

首先无论重物球质量怎么变化，系泊系统的总高度是不变的一直等于海水的深度为18m，即

HHbHgHtHl

其中

HbhHglgcosii14

Htltcos5Hl622.50dssind5.2.2.2.2 锚链与锚夹角分析

根据题目的约束条件锚链和锚的夹角不能大于16度，因此

16

综上所述，我们就建立了基于重物球质量调节的系泊系统优化设计模型。在保证以下这两个约束条件成立的情况下

HHbHgHtHl 16我们要寻找

minM

这样的一个设计最优化的系泊系统。 5.2.3 问题二模型的求解

1．问题二的具体求解思路和问题一类似还是采用搜索算法，但是由于现在搜索的变量变成了两个，计算量大大增加，因此我们考虑采用变步长的方法来求得近似最优解[4]。 2．和问题一类似，在参数的有效范围内给定一个较大的步长，h,G,计算得到每一组数据的H和。

3．首先判断得到的H和是否同时满足H与18的差值小于一个十分小的量（比如说这个量取0.0.1），并且的角度小于16度，在满足这两个条件并满足上文约束条件的数据中寻找M最小的数据，记录该数据的G与h

4．以在上步得到的G与h为中心，取一个较小的范围，选定一个较小的步长h,G再次计算H和，然后重复步骤2

直到某一次循环找到的最小的M大于上一次循环找到的最小的M，或者达到了规定的循环次数后，跳出，输出结果 5.2.4 问题二求解结果展示

根据以上思路以及问题二的已知条件，我们得到了相对最优的系泊模型设计，结果如下

表4 风速为36 m/s 时的系泊系统最优参数表

在该情况下锚链的形状如下

图11 风速为36m/s时最优系泊设计锚链形状图

5.3 问题三模型的建立与求解

问题一和问题二都是在水深，风速恒定，水流静止的特殊情况下建立的，那么到了问题三，题目给定水深，风速和水流的范围，考虑实际系泊系统中水流速度的影响，在水速，风速，水深不定的情况下，选择合适的锚链型号，长度以及重物球的质量，来得到最优的系泊系统设计。

根据上文所述，问题三的决策目标和问题二是相同的，都是选择合适的锚链型号，长度以及重物球的质量，不同的是由于水流速度的加入，我们要重新分析系泊系统的条件，在实际情况中，水平风速和水流力方向夹角是不定的，但是在此为了简化计算的复杂度，下文仅在水流力和水平风速同向的情况下进行分析。 5.3.1 当水流力和水平风速同向时系泊系统受力分析 5.3.1.1 浮标受力分析

对于浮标的受力分析，首先我们考虑实际情况，依旧假设浮标始终保持竖直状态不会倾斜，在此基础上我们进行浮标的受力分析。浮标一共受到5个力的作用，分别是风力F，浮力fb，钢管的拉力T1，浮标自身的重力Gb，水流力Wb如图11所示：

图12 浮标受力分析图

根据理论力学的知识，要使浮标受力平衡，应该让浮标在水平方向和竖直方向上都保持平衡，可得到如下表达式

fb=T1cos1GbFWbT1sin1其中根据阿基米德定律可知

fbgV

V是浮标浸入水面的体积易得为h，其中h是浮标的吃水深度。

风力由题目已知条件可得，vF为风速，SF为物体在风向法平面的投影面积可由简单的数学计算得到

F0.625SFvF2

水流力也可由题目已知条件得到，vw为水流的速度，Swater为物体在水流速度法平面的投影面积，经过计算得到Swater与物体在水中的横截面积的偏角小于2度，可以忽略不计，因此我们就用物体在水中的横截面积作为Swater。

F374SWvwater2

5.3.1.2钢桶受力分析

如图12所示我们对钢桶和小球进行受力分析，图中DC就为钢桶，小球看成一个点处在D的位置。

图13 钢桶受力分析图

如图所示，钢桶和小球一共受到6个力的作用，钢桶左侧受到的拉力T5,右侧受到的拉力T6，钢桶受到的浮力ft，以及它的重力Gt和小球对钢桶的拉力Tzwq和水流对它的作用力Wt根据受力平衡分析,可以得到以下表达式：

GtTzwq+T6cos6ftT5cos5T5sin5+WtT6sin6

然后考虑钢桶与重物球的力矩平衡，由于钢桶的重力以及浮力还有水流力作用在质心所以不予考虑，我们只要考虑钢桶与重物球左侧受到的拉力T5,右侧受到的拉力T6，以及小球对钢桶的拉力Tzwq，根据力矩平衡的公式可以得到：

0.5T5Ltsin(55)0.5TzwqLtcos(905)0.5T6Ltsin(65)

其中重物球对钢桶的拉力可以由图四得到

图14 重物球受力分析图

根据竖直方向的受力分析我们可以很简单的得到

TzwqfzwqGzwq

5.3.1.3钢管受力分析

对于钢管的受力分析，和问题一类似，4节钢管的受力分析类似，我们在此就分析其中一节钢管，其它刚管依次类推。钢管14如图所示，一共受到来自它左侧的拉力，来自它右侧的拉力，它自身的重力，浮力以及水流力4个力。

图15 钢管受力分析图

根据钢管的力矩平衡和受力平衡我们可以得到以下表达式：

Gg1T2cos2fg1T1cos1 T1sin1Wg1T2sin20.5T2Lgsin(11)0.5T1Lgsin(21)依次类推我们可以得到钢管2的表达式

Gg2T3cos3fg2T2cos2 T2sin2Wg2T3sin30.5T3Lgsin(22)0.5T2Lgsin(32)同理可得钢管3的表达式

Gg3T4cos4fg3T3cos3 T3sin3Wg3T4sin40.5T4Lgsin(33)0.5T3Lgsin(43)同样可以得到钢管4的表达式

Gg4T5cos5fg4T4cos4 T4sin4Wg4T5sin50.5T5Lgsin(44)0.5T4Lgsin(54)5.3.2 基于重物球质量变初始条件最优系泊系统模型的确立与求解 5.3.2.1 模型的确立

根据以上的受力分析运用于问题二所类似的分析方法我们就可以建立同样的目标函数M

M约束条件为

55maxhhmaxSSmax

HHbHgHtHl 16需要强调的是问题三与问题二唯一的区别就在于推导出上述目标函数和约束条件的基本表达式由于考虑了水流力的影响发生了改变。 5.3.2.2 模型的具体求解思路

问题三和问题二求解思路基本相同，仍然采用搜索算法，只是我们需要考虑的变量新增了锚链的型号和锚链的个数，变量从两个变成了四个，首先我们假设锚链的型号为P，锚链链环的个数为N，具体的求解过程如下

1.和问题二类似，在参数的有效范围内给定一个较大的步长，h,G,P,N(其中锚链的型号从1到4循环)，计算得到每一组数据的H和。

2．首先判断得到的H和同时满足H与18的差值小于一个十分小的量（比如说这个量取0.01），并且的角度小于16度，在满足这两个条件并满足上文约束条件的数据中我们改变一开始随意设定的海水深度（在题目所给的16~20m的海水深度范围内用随机函数任意取100个高度）测试这些数据是否依旧满足上面的两个条件，从筛选出来的数据里，寻找M最小的数据，记录该数据的G,P,N与h

3.以在上步得到的G,P,N与h为中心，取一个较小的范围，选定一个较小的步长

h,G,P,N再次计算H和，然后重复步骤2。

4.直到某一次循环找到的最小的M大于上一次循环找到的最小的M，或者达到了规定的循环次数后，跳出，输出结果。 5.3.4 模型求解结果

首先我们假设风速为24m/s，水速为1.5m/s，在此基础上我们按照上文的思路进行求解得到了最优的系统设计为重物球的质量为3900kg，锚链型号是3号，使用的锚链总数是250个，接下来给出水深为16,18,20米的情况来验证该最优设计的正确性。 5.3.4.1 水深为16m时详解

水深为16m时系泊系统的具体参数如下

表5 水深16m，风速24 m/s，水速1.5m/s时的系泊系统优化参数表

由此可以得到锚链的形状

图16 水深16m，风速24 m/s，水速1.5m/s时的锚链形状图

5.3.4.2 水深为18m时详解

水深为18m时系泊系统的具体参数如下

表6 水深18m，风速24 m/s，水速1.5m/s时的系泊系统优化参数表

由此可以得到锚链的形状

图17 水深18m，风速24 m/s，水速1.5m/s时的锚链形状图

5.3.4.3 水深为20m时详解

水深为20m时系泊系统的具体参数如下

表7 水深20m，风速24 m/s，水速1.5m/s时的系泊系统优化参数表

由此可以得到锚链的形状

图18 水深20m，风速24 m/s，水速1.5m/s时的锚链形状图

六、模型的评价

6.1 模型的优点

1．通过对系泊系统各个部分严谨的受力分析，得出了比较完美的受力平衡关系式，进而建立了比较完善的模型，思路严晰，结构严谨。

2．以浮标吃水深度最小，钢桶倾斜角度最小和游动范围最小为优化目标，建立多目标优化模型。可以根据工程需要，通过调整权重来衡量每个目标的重要性，得到更满足需求的结果。

3．模型求解时，给出一系列锚链型号，锚链长度，重物球重力和浮标吃水深度，通过

搜索算法找到最优解，简化了求解的难度。 6.2 模型的缺点

1．未知量较多，模型求解运算量太大。

2．当风力和水流力同时作用在浮标上时，没有考虑浮标的倾斜的情况。

七、参考文献

[1] 卓金武.MATLAB在数学建模中的应用（第二版）.北京:北京航空航天大学出版社，

2014.9

[2] 姜启源.数学模型（第四版）.北京:高等教育出版社，2011.1 [3] 司守奎.数学建模算法与应用 .北京: 国防工业出版社，2011.8

[4] 张平. MATLAB基础与应用 .北京:北京航空航天大学出版社，2007.5 [5] 刘来福.数学建模方法与分析（第三版）.北京:机械工业出版社，2009.5 [6] 姚泽清，郑旭东，赵颖.全国与优秀论文评析.北京:国防工业出版社，2012.5 大学生数学建模竞赛

八、附录

程序源码

%Q1求解程序 clc clear

%-------------全局变量定义---------------------- global ro g v Gb Gg fg Gt ft Tzwq yita lianjie len; ro = 1025; g = 9.8; v = 36;

Gb = 1000 * g; Gg = 10 * g;

fg = 0.025 * 0.025 * pi * ro * g; Gt = 100 * g;

ft = 0.15 * 0.15 * pi *ro * g; %Tzwq = 1200 * g; Tzwq = 38220; yita = 7;

lianjie = 0.120; len = 30;

%------------------------------------------------ num = len / lianjie;

h = 1.660000000000000; for i = 1:length(h)

T1(i) = sqrt((Gb - ro * g * pi * h(i)).^2 + (0.625 * (4 - 2 * h(i)) * v .^ 2 ).^2); alpha1(i) = atan(-(0.625 * (4 - 2 * h(i)) * v .^ 2)/(Gb - ro * g * pi * h(i)));

T2(i) = sqrt((Gg - fg - T1(i) * cos(alpha1(i))).^2 + (T1(i) * sin(alpha1(i))).^2); alpha2(i) = asin(T1(i) * sin(alpha1(i)) / T2(i));

sita1(i) = atan((T1(i) * sin(alpha1(i))+T2(i)*sin(alpha2(i)))/(T1(i) * cos(alpha1(i))+T2(i)*cos(alpha2(i))));

T3(i) = sqrt((Gg - fg - T2(i) * cos(alpha2(i))).^2 + (T2(i) * sin(alpha2(i))).^2); alpha3(i) = asin(T2(i) * sin(alpha2(i)) / T3(i));

sita2(i) = atan((T2(i) * sin(alpha2(i))+T3(i)*sin(alpha3(i)))/(T2(i) * cos(alpha2(i))+T3(i)*cos(alpha3(i))));

T4(i) = sqrt((Gg - fg - T3(i) * cos(alpha3(i))).^2 + (T3(i) * sin(alpha3(i))).^2); alpha4(i) = asin(T3(i) * sin(alpha3(i)) / T4(i));

sita3(i) = atan((T3(i) * sin(alpha3(i))+T4(i)*sin(alpha4(i)))/(T3(i) * cos(alpha3(i))+T4(i)*cos(alpha4(i))));

T5(i) = sqrt((Gg - fg - T4(i) * cos(alpha4(i))).^2 + (T4(i) * sin(alpha4(i))).^2); alpha5(i) = asin(T4(i) * sin(alpha4(i)) / T5(i));

sita4(i) = atan((T4(i) * sin(alpha4(i))+T5(i)*sin(alpha5(i)))/(T4(i) * cos(alpha4(i))+T5(i)*cos(alpha5(i))));

T6(i) = sqrt((Gt + Tzwq - ft - T5(i) * cos(alpha5(i))).^2 + (T5(i) * sin(alpha5(i))).^2); alpha6(i) = asin(T5(i) * sin(alpha5(i)) / T6(i));

sita5(i) = atan((T5(i) * sin(alpha5(i))+T6(i)*sin(alpha6(i)))/(T5(i) * cos(alpha5(i))+T6(i)*cos(alpha6(i)) + Tzwq ));

for n = 1:num if n == 1

T7y(i,n) = T6(i) * cos(alpha6(i)); T7x(i,n) = T6(i) * sin(alpha6(i)); alpha7(i,n) = atan(T7y(i,n)/T7x(i,n)); tx(i,n) = lianjie * cos(alpha7(i,n)); ty(i,n) = lianjie * sin(alpha7(i,n)); continue end

T7y(i,n) = T7y(i,n - 1) - lianjie * 12.5 * g; T7x(i,n) = T7x(i,n - 1);

alpha7(i,n) = atan(T7y(i,n)/T7x(i,n)); if lianjie * sin(alpha7(i,n)) <= 0 tx(i,n) = tx(i, n - 1) + lianjie; ty(i,n) = ty(i, n - 1); continue end

tx(i,n) = lianjie * cos(alpha7(i,n)) + tx(i, n - 1); ty(i,n) = lianjie * sin(alpha7(i,n)) + ty(i, n - 1); end

sx(i) = tx(i,n) + sin(sita5(i)) + sin(sita4(i)) + sin(sita3(i)) + sin(sita2(i)) + sin(sita1(i)); sy(i) = ty(i,n) + cos(sita5(i)) + cos(sita4(i)) + cos(sita3(i)) + cos(sita2(i)) + cos(sita1(i)); end

for i = 1 : length(h)

sita1(i) = sita1(i) / pi * 180; sita2(i) = sita2(i) / pi * 180; sita3(i) = sita3(i) / pi * 180; sita4(i) = sita4(i) / pi * 180; sita5(i) = sita5(i) / pi * 180; end

%Q2求解程序 clc clear

%-------------全局变量定义----------------------

global ro g v Gb Gg fg Gt ft yita lianjie len trange mrange; ro = 1025; g = 9.8; v = 36;

Gb = 1000 * g; Gg = 10 * g;

fg = 0.025 * 0.025 * pi * ro * g; Gt = 100 * g;

ft = 0.15 * 0.15 * pi *ro * g; trange = 5 / 180 * pi; mrange = 16 / 180 * pi; yita = 7;

lianjie = 0.105; len = 22.05;

%------------------------------------------------ num = len / lianjie; Tzwqs = 2800:10:5100; Tzwqs = Tzwqs * g; h = 0.1:0.01:2; s = 1;

for t = 1:length(Tzwqs) Tzwq = Tzwqs(t); for i = 1:length(h)

T1(i) = sqrt((Gb - ro * g * pi * h(i)).^2 + (0.625 * (4 - 2 * h(i)) * v .^ 2 ).^2); alpha1(i) = atan(-(0.625 * (4 - 2 * h(i)) * v .^ 2)/(Gb - ro * g * pi * h(i)));

T2(i) = sqrt((Gg - fg - T1(i) * cos(alpha1(i))).^2 + (T1(i) * sin(alpha1(i))).^2); alpha2(i) = asin(T1(i) * sin(alpha1(i)) / T2(i));

sita1(i) = atan((T1(i) * sin(alpha1(i))+T2(i)*sin(alpha2(i)))/(T1(i) * cos(alpha1(i))+T2(i)*cos(alpha2(i))));

T3(i) = sqrt((Gg - fg - T2(i) * cos(alpha2(i))).^2 + (T2(i) * sin(alpha2(i))).^2);

alpha3(i) = asin(T2(i) * sin(alpha2(i)) / T3(i));

sita2(i) = atan((T2(i) * sin(alpha2(i))+T3(i)*sin(alpha3(i)))/(T2(i) * cos(alpha2(i))+T3(i)*cos(alpha3(i))));

T4(i) = sqrt((Gg - fg - T3(i) * cos(alpha3(i))).^2 + (T3(i) * sin(alpha3(i))).^2); alpha4(i) = asin(T3(i) * sin(alpha3(i)) / T4(i));

sita3(i) = atan((T3(i) * sin(alpha3(i))+T4(i)*sin(alpha4(i)))/(T3(i) * cos(alpha3(i))+T4(i)*cos(alpha4(i))));

T5(i) = sqrt((Gg - fg - T4(i) * cos(alpha4(i))).^2 + (T4(i) * sin(alpha4(i))).^2); alpha5(i) = asin(T4(i) * sin(alpha4(i)) / T5(i));

sita4(i) = atan((T4(i) * sin(alpha4(i))+T5(i)*sin(alpha5(i)))/(T4(i) * cos(alpha4(i))+T5(i)*cos(alpha5(i))));

T6(i) = sqrt((Gt + Tzwq - ft - T5(i) * cos(alpha5(i))).^2 + (T5(i) * sin(alpha5(i))).^2); alpha6(i) = asin(T5(i) * sin(alpha5(i)) / T6(i));

sita5(i) = atan((T5(i) * sin(alpha5(i))+T6(i)*sin(alpha6(i)))/(T5(i) * cos(alpha5(i))+T6(i)*cos(alpha6(i)) + Tzwq ));

for n = 1:num if n == 1

T7y(i,n) = T6(i) * cos(alpha6(i)); T7x(i,n) = T6(i) * sin(alpha6(i)); alpha7(i,n) = atan(T7y(i,n)/T7x(i,n)); tx(i,n) = lianjie * cos(alpha7(i,n)); ty(i,n) = lianjie * sin(alpha7(i,n)); continue end

T7y(i,n) = T7y(i,n - 1) - lianjie * 7 * g; T7x(i,n) = T7x(i,n - 1);

alpha7(i,n) = atan(T7y(i,n)/T7x(i,n)); if lianjie * sin(alpha7(i,n)) <= 0 tx(i,n) = tx(i, n - 1) + lianjie; ty(i,n) = ty(i, n - 1); continue end

tx(i,n) = lianjie * cos(alpha7(i,n)) + tx(i, n - 1); ty(i,n) = lianjie * sin(alpha7(i,n)) + ty(i, n - 1); end

sx(i) = tx(i,n) + sin(sita5(i)) + sin(sita4(i)) + sin(sita3(i)) + sin(sita2(i)) + sin(sita1(i)); sy(i) = ty(i,n) + cos(sita5(i)) + cos(sita4(i)) + cos(sita3(i)) + cos(sita2(i)) + cos(sita1(i));

if abs(sy(i) - 18) < 0.5 && alpha7(i,n) < mrange && sita5(i) < trange he(s,1) = Tzwq; he(s,2) = h(i); if s == 1

youhua = sx(i) + h(i) + sita5(i);

he1 = [Tzwq, h(i)];

else if sx(i) + h(i) + sita5(i) < youhua youhua = sx(i) + h(i) + sita5(i); he1 = [Tzwq, h(i)]; end end

s = s + 1; end end end

%Q3求解程序

%-------------全局变量定义----------------------

global ro g v Gb Gg fg Gt ft yita lianjie len shui trange mrange; ro = 1025; g = 9.8; v = 36;

Gb = 1000 * g; Gg = 10 * g;

fg = 0.025 * 0.025 * pi * ro * g; Gt = 100 * g;

ft = 0.15 * 0.15 * pi *ro * g; yita = 7;

lianjie = 0.120; len = 30;

shui = 374 * 1.5 ^ 2; trange = 5 / 180 * pi; mrange = 16 / 180 * pi;

%------------------------------------------------ num = len / lianjie;

Tzwqs = 1200:100:8000; Tzwqs = Tzwqs * g; h = 0.1:0.01:2; s = 1;

for t = 1:length(Tzwqs) Tzwq = Tzwqs(t); for i = 1:length(h)

T1(i) = sqrt((Gb - ro * g * pi * h(i)).^2 + (0.625 * (4 - 2 * h(i) ) * v .^ 2 + shui * h(i) * 2).^2); alpha1(i) = atan(-(0.625 * (4 - 2 * h(i)) * v .^ 2 + shui * h(i) * 2)/(Gb - ro * g * pi * h(i)));

T2(i) = sqrt((Gg - fg - T1(i) * cos(alpha1(i))).^2 + (T1(i) * sin(alpha1(i)) + 0.05 * shui).^2); alpha2(i) = asin((T1(i) * sin(alpha1(i)) + 0.05 * shui) / T2(i));

sita1(i) = atan((T1(i) * sin(alpha1(i))+T2(i)*sin(alpha2(i)))/(T1(i) * cos(alpha1(i))+T2(i)*cos(alpha2(i))));

T3(i) = sqrt((Gg - fg - T2(i) * cos(alpha2(i))).^2 + (T2(i) * sin(alpha2(i)) + 0.05 * shui).^2); alpha3(i) = asin((T2(i) * sin(alpha2(i)) + 0.05 * shui) / T3(i));

sita2(i) = atan((T2(i) * sin(alpha2(i))+T3(i)*sin(alpha3(i)))/(T2(i) * cos(alpha2(i))+T3(i)*cos(alpha3(i))));

T4(i) = sqrt((Gg - fg - T3(i) * cos(alpha3(i))).^2 + (T3(i) * sin(alpha3(i)) + 0.05 * shui).^2); alpha4(i) = asin((T3(i) * sin(alpha3(i)) + 0.05 * shui) / T4(i));

sita3(i) = atan((T3(i) * sin(alpha3(i))+T4(i)*sin(alpha4(i)))/(T3(i) * cos(alpha3(i))+T4(i)*cos(alpha4(i))));

T5(i) = sqrt((Gg - fg - T4(i) * cos(alpha4(i))).^2 + (T4(i) * sin(alpha4(i)) + 0.05 * shui).^2); alpha5(i) = asin((T4(i) * sin(alpha4(i)) + 0.05 * shui) / T5(i));

sita4(i) = atan((T4(i) * sin(alpha4(i))+T5(i)*sin(alpha5(i)))/(T4(i) * cos(alpha4(i))+T5(i)*cos(alpha5(i))));

T6(i) = sqrt((Gt + Tzwq - ft - T5(i) * cos(alpha5(i))).^2 + (T5(i) * sin(alpha5(i)) + shui * 0.3).^2); alpha6(i) = asin((T5(i) * sin(alpha5(i)) + shui * 0.3) / T6(i));

sita5(i) = atan((T5(i) * sin(alpha5(i))+T6(i)*sin(alpha6(i)))/(T5(i) * cos(alpha5(i))+T6(i)*cos(alpha6(i)) + Tzwq ));

for n = 1:num if n == 1

T7y(i,n) = T6(i) * cos(alpha6(i)); T7x(i,n) = T6(i) * sin(alpha6(i)); alpha7(i,n) = atan(T7y(i,n)/T7x(i,n)); tx(i,n) = lianjie * cos(alpha7(i,n)); ty(i,n) = lianjie * sin(alpha7(i,n)); continue end

T7y(i,n) = T7y(i,n - 1) - lianjie * 12.5 * g; T7x(i,n) = T7x(i,n - 1);

alpha7(i,n) = atan(T7y(i,n)/T7x(i,n)); if lianjie * sin(alpha7(i,n)) <= 0 tx(i,n) = tx(i, n - 1) + lianjie; ty(i,n) = ty(i, n - 1); continue end

tx(i,n) = lianjie * cos(alpha7(i,n)) + tx(i, n - 1); ty(i,n) = lianjie * sin(alpha7(i,n)) + ty(i, n - 1); end

sx(i) = tx(i,n) + sin(sita5(i)) + sin(sita4(i)) + sin(sita3(i)) + sin(sita2(i)) + sin(sita1(i)); sy(i) = ty(i,n) + cos(sita5(i)) + cos(sita4(i)) + cos(sita3(i)) + cos(sita2(i)) + cos(sita1(i)); if abs(sy(i) - 16) < 0.5 && alpha7(i,n) < mrange && sita5(i) < trange he2(s,1) = Tzwq; he2(s,2) = h(i); s = s + 1;

end end end

%Q3求解程序

%-------------全局变量定义----------------------

global ro g v Gb Gg fg Gt ft yita lianjie len shui trange mrange; ro = 1025; g = 9.8; v = 36;

Gb = 1000 * g; Gg = 10 * g;

fg = 0.025 * 0.025 * pi * ro * g; Gt = 100 * g;

ft = 0.15 * 0.15 * pi *ro * g; yita = 7;

lianjie = 0.120; len = 30;

shui = 374 * 1.5 ^ 2; trange = 5 / 180 * pi; mrange = 16 / 180 * pi;

%------------------------------------------------ num = len / lianjie; h = 1.65; Tzwq = 38220;

for i = 1:length(h)