2019年高教社杯全国大学生数学建模竞赛B题代码和思路

同心鼓运动模型的理论分析与策略研究

摘要

“同心协力”是一项团队协作能力拓展项目，有较广泛的普及面。本文对同心鼓的颠球策略进行研究，分析同心鼓与排球的碰撞过程，综合考虑队员牵拉同心鼓的实际操作、边界和初始条件，在队员水平牵拉绳的假设下，建立球鼓碰撞模型和同心鼓运动模型，该模型可用于计算在特定情况下同心鼓的运动参数，并在鼓面发生倾斜时，给出相应调整策略。

对于问题一：通过假设队员的手仅在水平方向上移动，将同心鼓的三维运动模型进行简化。考虑队员施加在绳子上的力大小恒定，建立基于牛顿第二定律的二阶常系数微分方程，确定初边值条件，得到描述同心鼓的高度及速度的方程。基于动量守恒定理与能量守恒定理，考虑到同心鼓及排球碰撞时发生的能量损失，建立起队员拉绳所消耗的功率的数学描述。由于不同参赛队员的用力偏好差异性极大，在给定队员所偏好的用力大小及频率下，可求得相对应的最优方案。在本文中，我们以参赛队员的耗能应尽可能地小为目标给出团队的最佳协作策略。

对于问题二：基于欧拉旋转定理，鼓的运动可描述为沿 Z 轴方向的平动及绕 X轴、Y轴转动。由于同心鼓的转角 α、β都为小量，我们讨论并分析了小角度转动与转动顺序无关的条件，将模型三维转动求解转化为两个自由度上的独立求解问题。基于牛顿第二定律及转动定律，得到0.1s 时，九种工况下鼓面的倾斜角度分别为：0.2124°0.3885°0.1649°0.2270°0.4123°0.1749°0.2175°0.1359°0.3401°。 从理论和结果进行分析，得到转动顺序对结果影响占次要因素，模型可分解为两个自由度的独立转动。

对于问题三：由问题二建立的模型可知，在非理想状态下，发力大小与时机的不精准控制将导致同心鼓发生小角度倾斜，因此我们针对问题一所建立的模型存在一定的不完善性，需要进行相应的调整。我们分别从调整参赛队员的拉绳拉力大小、调整鼓面倾斜状态、调整队员的水平位置三个方面整体地对模型进行优化，分析并讨论了各个方面的调整对最终结果的影响，以给出在现实情形中适应性能更强的团队协作策略。

对于问题四：本文针对给定排球偏斜竖直方向运动的特殊情况，分析调整鼓面倾角与倾向的解决方法，并以此建立多目标规划的最佳策略模型，基于模拟退火算法的相关思想，通过分析约束条件及目标函数，针对团队协作策略的有效性与协调性进行全局优化，并对优化结果进行分析处理，得到了在题设情况下队员的拉绳策略。并以此分析在现实情形中这种调整策略的实施效果。

综上所述，本文综合分析了同心鼓在不同情况下的运动姿态，得到了有效地调整碰撞过程中所产生的偏差的方法，并且该方法给出了一种适应不同情况下团队的协作策略，具有较好的普适性。

关键词：同心鼓，欧拉刚体旋转定理，模拟退火算法，鼓面倾角调整，最优协作策略

1 问题重述

“同心协力”(又称“同心鼓”) 是一项团队协作能力拓展项目。该项目用一个排球(质量0.27kg)和一只牛皮双面鼓(质量3.6kg，鼓面直径为40cm) 作为道具，并且在同心鼓鼓身(高度为 22cm)中部固定有多根长度相同且沿圆周均匀分布的绳子。由不少于 8名的参赛队员每人牵拉一根绳子，使鼓面在空间内作以竖直方向运动为主，水平面内的横向运动及水平面外的转动为辅的运动。通过同心鼓的主运动将排球颠至高度 40cm以上，否则项目停止。重复上下颠球过程，要求队员最终连续颠球次数最多。

结合现实中的亲身游戏体验以及网上优秀团队的颠球视频，笔者发现该项目的难度主要在于不同参赛队员需要在不同方向对同心鼓同时进行合理有效的控制，进而易出现颠球高度和方向难以把握的问题。这与不同参赛队员的发力大小、方向、时机、手法以及走位等因素相关。现实中物理过程极其复杂，如何简化过并建立起合理的数学模型是笔者研究的重点。

为了在理论上进一步分析同心鼓项目的玩法，并得到合理且可操作性较高的颠球技巧与合作策略，笔者将综合分析参赛队员、同心鼓及排球的运动过程，抓住主要因素，忽略次要因素，进而建立合理的力学简化模型，并利用该模型对以下几个问题进行研究：

·在每位参赛人员均能精准控制用力方向、时机及力度大小的前提下，通过理论分析并与实际经验相结合的办法，给出最佳团队协作的颠球策略，以及对应的颠球高度。

·结合相关文献，利用图解的方法发现并展示物体在三维空间中旋转的不可交换性问题，确定在小角度的情况下，这种空间旋转交换理论的误差值极小，可满足精度要求。且需提出合理可行的状态量，用于描述空间小角度的倾斜状态。

·在非理想状态下，参赛队员的力度和发力时机均可能发生偏差，基于上述研究成果与物理简化模型，计算九种不同工况下的鼓面倾斜角度。

·在非理想的状态下，对团队颠球策略进行理论上的合理性分析以及现实情况中的可行性分析，并对其进行合理的修正。

·若排球因鼓面倾斜而发生偏移后，给出队员发力大小及时机的合理方案，使得排球重新回到竖直的直线运动。

综合上述问题，查找相关文献及数据，简化同心鼓的颠球过程，设计合理的物理模型，分析并得到合理的团队协作策略，得到小角度倾斜情况下用于描述的三个状态量，探究团队协作中作用力大小与发力时机出现偏差对鼓面的倾角与倾向的物理关系，并针对实际情况下出现的排球偏移问题，给出参赛团队应采取的合理应对措施。结合以上分析与研究成果，撰写一篇研究论文。

2 模型假设

1. 在分析研究鼓体及排球运动过程中的受力时，忽略排球和鼓的运动以及牵引绳振动中所产生的空气阻力的影响。

2. 假设牵引同心鼓运动的绳子是理想绳，即绳子不可伸长，且绳子为轻质绳，绳所提供的拉力始终沿着绳。

3. 假设队员每次拉绳的力恒定，且在其发力时机时瞬时将拉绳的力变化为新的受力大小，且力是瞬间施加在鼓上的。

4. 假设参与拉绳的队员均匀对称地分布在以同心鼓鼓面中心为圆心的圆周上，鼓的中心始终在过圆心的竖直轴上运动，鼓的水平运动由队员整体的水平运动提供，不考虑鼓相对于队员的水平运动。

5. 假设鼓与排球碰撞时，排球始终保持平动状态，且忽略排球与鼓面非对心碰撞时所产生的瞬时冲量矩。

6. 忽略鼓与球碰撞之后，队员恢复初始拉力时鼓的复杂运动，且假设每次鼓与球碰撞后鼓都能够恢复到初始水平状态。

7. 为保证比赛的稳定且有序的进行，假设队员在调整好颠球的状态之后，更倾向于有规律的周期性颠球，即其每两次发力之间的时间间隔与排球每次被颠起并落下的时间相同，并称该时间为一个运动周期。

3 符号约定

符号	说明
M	“同心鼓”的质量
9	重力加速度， 这里取g=9.8m/s²
n	每次参与“同心鼓”比赛队员的人数
L	牵引“同心鼓”运动的绳子的长度
θ₀	初始时刻绳子与水平方向的夹角
θ	任意时刻绳子与水平方向的夹角
x	同心鼓鼓心相对于拉绳平面的垂直距离
a	同心鼓的加速度
H	击鼓后排球的反弹高度
tcrash	排球与同心鼓碰撞的时刻
Ucrash	排球与同心鼓碰撞时同心鼓的速度
e	弹性恢复系数
α	xoy 平面绕 y 轴转动角度(方向满足左手螺旋法则)
β	xoy 平面绕 x 轴转动角度(方向满足右手螺旋法则)
γ	鼓面发生倾斜时与竖直方向夹角大小

4 团队最佳协作策略之水平拉绳模型

4.1 模型分析

题目要求我们在精确控制用力方向、时机和力度的理想状态下，给出颠球的最优协作策略及在该策略下的颠球高度。

本题目的重点是要将球撞击同心鼓并弹起这一过程进行物理模型的简化，该碰撞过程为接近弹性碰撞的非弹性碰撞模型，每次都会将体系中小部分能量通过碰撞的形式耗散，而队员拉鼓为整个体系补充消耗的能量。为了达到最佳的协作策略，队员需在使得连续颠球的次数尽可能多的前提下，减少能量的消耗，获得更佳的舒适度和游戏体验感。

现实中人通过拉动绳子而使同心鼓做上升运动的过程十分复杂，根据主次分析法，笔者将人手的运动轨迹理想化为水平方向的运动，且每位队员手的高度位于同一水平面上。为了更好地利用能量，需假设在鼓速度最大时与排球发生碰撞，且碰撞后小球仅发生速度方向的改变，而速度大小不发生改变。此时，小球弹起的高度与队员所施加的力以及力持续的时间可以建立起物理关系，从而得到最佳的协作策略。

4.2  模型的建立与求解

为了更好地控制绳的方向，以便于团队协作，在考虑到现实情况的基础上，笔者假设每位队员的手仅在水平方向上移动，而无竖直方向的运动，且每位队员手的高度位于同一水平面上。此时以队员手所在的平面为基准面，向下方向作为坐标的正方向，建立同心鼓的位置坐标系。假设每次参赛的队员人数为n，每根绳子的长度为 L。初始静止状态时，n位队员对绳施加恒力 F₀ 用以平衡同心鼓的重力G = Mg，则有初始平衡方程：

 n⋅F₀sinθ₀=Mg

       其中，初始位置绳与水平方向的夹角满足  sinθ0=H0L,代入上式可得  H0=MgLnF0。

到达一定时刻时，队员通过手在同一水平面内的移动，在瞬间对绳增大拉力，此阶段时间可忽略不计。同心鼓在拉力的作用下，离开初始位置向上运动，绳与水平方向的夹角θ随着同心鼓的位置坐标的改变而改变，即有 θ=sin-1xL,则“同心鼓”的运动状态可以由下图所表示：

为了更合理地简化碰撞模型，假设队员从拉绳到同心鼓与排球发生碰撞前，施加在绳子上的力的大小是恒定的，且不考虑绳子的重力，则由牛顿第二定律 ∑F=Ma可得：

 -n⋅F₀sinθ+Mg=Ma

即

 d2xdt2=a=-nFsinθ+MgM=g-nFsinθM

由此，可以建立同心鼓的坐标x关于时间t的二阶微分方程，易知初始时刻同心鼓的坐标为 H₀,速度

为v₀=0, 即

 d2xdt2=9-nFMLx,x0=H0=M2n0,x0=v0=0.

上述含有的初值的二阶常系数微分方程具有其理论解，可解得同心鼓的位置随时间的关系函数x(t)为

                          xt=MgLnF+MgLnFF0F-F0⋅cosnFMLt                       (1)

对上式求导，可得同心鼓的运动速度随时间的关系函数 v(t) 为

                         vt=dxdt=-gF0⋅nFMLF-F0⋅sinnFMLt                     (2)

其同心鼓运动随时间的关系图可以由下图表示

考虑到实际情况，若在同心鼓的速度最大时，同心鼓与排球发生碰撞，队员可在用力最小的情况下使球达到尽可能高的高度。因此，笔者假设在同心鼓的速度最大时，鼓与球发生碰撞，即碰撞时刻为 tcrash=π2nFML,此时

 vcrash=nFML⋅gF0⋅F-F0

当同心鼓与排球发生碰撞时，考虑到同心鼓与排球碰撞时会发生能量损失，需引入恢复系数的概念，碰撞恢复系数是碰撞前后两物体接触点的法向相对分离速度与法向相对接近速度之比，即

 e=|v'1-v'2v1-v2|

根据同心鼓及排球的材料性质，碰撞十分接近弹性碰撞，可假设e=0.95。同心鼓的质量为 M，碰撞前后的速度分别为 v₁、v₁';排球的质量为 m，碰撞前后的速度分别为 v₂、v₂，故可联立方程：

版权归全国大学生数学建模竞赛所有，如有侵权请联系删除

 v1'-v2'=-e⋅v1-v2,Mv1+mv2=Mv1'+mv2'.

易得到碰撞后同心鼓和排球的速度 v'₁、v'₂分别如下

 ν1'=m1-em2ν1+1+em2ν2m1+m2,ν2'=m2-em12+1+em1v1m1+m2.

笔者认为队员对同心鼓施加作用力的过程可简化以下四个阶段：第一阶段( O-t₁为使同心鼓加速上升阶段，队员施加的恒力 F； 第二阶段( t₁-t₂为碰撞完成后队员对鼓施加非线性变化的力； 第三阶段 t₂-t₂+△t为同心鼓到达最低点后，队员为了使同心鼓停止运动的突加力，其作用时间 △t 极短，这一阶段可忽略不计； 第四阶段( t₂-t₃为同心鼓面位于最低处时，为了平衡同心鼓的重力而施加的恒力 F₀，故整个阶段如下图示

因此，为了达到最佳的协作策略，需保证每个队员在用力较小的情况下，尽可能减小周期运动的频率，使队员达到一个较为舒适的状态。已知排球与同心鼓发生碰撞时， v1=2gH,v2=vcrash。为了得到稳定的颠球高度，需保证碰撞过程仅改变排球的速度方向，而不改变其大小，即 v₁'=-v₁。进而可建立 H、F₀、F 之间的关系 H = H(F₀,F), 即

 1-2m11+em22gH=MLnF⋅gF0F-F0

接下来我们分析一下该情形下的最优策略，我们希望在每个碰撞周期中所发生的总能量的损失应尽可能地小，参考文献[1]，我们可以令 E 表示一次碰撞同心鼓与排球体系总能量的损失，则

 E=12m1v1'2+12m2v2'2-12m1v12-12m2v22=m1m22m1+m221-e2v1-v22

即有

 E=12mu21-e2

在颠球上下运动的一个周期内，假设参赛队员所提供的能量全部用以同心鼓与排球发生碰撞时所消耗的能量。因此，为了达到最佳的协作策略，需保证每个队员消耗的能量较小，笔者认为将每位队员在一个周期内所用的功率 P 作为评判指标较为合理可行，即

 P=En⋅T

由于现实中，不同参赛队员的用力偏好差异性极大，即一些参赛队员可能偏向于每次对同心鼓作用较小的力，但是在此情况下会使用力的频率相应增大；但是另外一些参赛队员可能更偏向于每个作用周期内，用力的频率较小，但是在此情况下每次需要对同心鼓作用较大的力。由此可见，队员的每次对同心鼓作用力的大小及频率为相互制约的关系，但是对于给定的用力的大小(或者用力的频率) 的情况，可通过上文所描述的公式求得相对应的用力的频率(或者用力的大小)。题设要求排球每次被颠起的高度应离开鼓面40cm以上，通过逻辑推理易知，在给定的参赛队员所偏好的用力的大小及用力的频率下，颠球高度越低，参赛队员所需的能量越小，即在该策略下，颠球高度 40cm为佳。

5 鼓面小角度倾斜求解之刚体定点旋转模型

5.1  模型分析

在现实情形中，队员发力时机和力度都存在一定误差，鼓面可能会出现倾斜而使球的位置发生偏移，题目要求我们不同队员出现失误的情况下，求出鼓面的倾斜角度。

假设同心鼓为刚体，本题目的重点是计算得到不同大小与不同时机的力作用下，同心鼓在空间中的姿态。刚体在空间中的大角度转动不具有交换性，但是空间小角度转动与次序无关，因此本题可将同心鼓的转动分解成多个方向的转动，再进行叠加求和。每位队员用力相等时，同心鼓上下运动而不会发生偏转；当队员在同心鼓上作用较大的力时，则所有的力对转动点的矩不再能够平衡，鼓面将发生转动；当队员作用力提早 0.1s 施加，则前0.1s内，队员对同心鼓的作用可视作提早施力的队员作用较大的力，而其他队员仍保持用以平衡同心鼓重力的恒力F0，使同心鼓在零时刻后以一个初始偏转角进行运动。因此，在第一问所建立的模型基础上，可通过同心鼓的运动微分方程，得到不同工况下 0.1s时同心鼓鼓面的倾斜角度。

5.2  模型建立与求解

根据欧拉提出欧拉旋转定理( Euler’ srotationtheorem)：刚体绕定点的任意有限转动可由绕过定点的某根轴的一次有限转动实现，即刚体绕过定点某根轴的一次有限转动可由绕定点的任意有限转动实现。因此，在第一问所建立的模型基础上，假设同心鼓仅有三个自由度，即沿 Z 轴方向的平动与绕 X轴、Y轴转动。笔者将同心鼓两个方向上的转动设想为刚体的两次相继有限转动，即同心鼓从起始位置(O- xoyozo) 出发, 先绕 -Y 轴转动α角到达 O-x₁y₁z₁,再绕 X 轴转动β角到达(O- xyz) 位置,则由三维空间坐标系变换，可得到同心鼓各次转动方向余弦矩阵分别为

 Iα=cosα0sinα010-sinα0cosα,Iβ=10sinβ0-sinβcosβ

在空间坐标系中，刚体连续作两次有限转动后的方向余弦矩阵可由各次转动方向余弦矩阵的乘法运算确定。

由此，可利用矢量的加法计算两次转动的结果。省略二阶小量，可将同心鼓转动的余弦矩阵[A] 简化为

 A=Iβ⋅Iα=Iα⋅Iβ=cosα0-sinβsinαsinβcosαcosβ

采用前面所述的水平拉绳团队协作策略，当参赛队员握住绳端并水平拉动时，恒力 F 会随着同心鼓鼓面上升而愈发趋于水平方向，在该合作策略下，结合上述刚体在微小转动时的空间转动可交换的性能，通过α、β两个状态量来描述鼓面的法相量的倾角γ

由此解得α、β，由上述旋转变换公式，可得同心鼓与水平面的夹角γ为

 γ=arcsinsin2α+sin2βsin2α+sin2β+cos2αcos2β≈arcsinsin2α+sin2β

对同心鼓整体做受力分析，其中同心鼓所受的总弯矩可由 M0=∑r×F1得到

因此，同心鼓在 X、Y、Z三个自由度方向上的运动可由微分方程表示为

 d2zdt2=-∑FsMd,dααdt2-1y≥2d24β2≤1J∑z,

同心鼓的绕x 和y轴的转动惯量的算法，由于同心鼓的对称性，则有： Jx=Jy。将同心鼓的外形抽象成为一个空心的圆柱，且圆柱壁的厚度相对于同心鼓的半径来说相对较小，因此鼓的质量集中在一个圆柱薄壁上，厚度不计。则质量密度为

 σ=M2πrH

则建立如上图示坐标系，则转动惯量 Jₓ=∫σl²dA,l 为对应单位质量点到转轴的距离。即

 Jx=σ-H2H2 dz02π r2sin2θ+z2rdθ=12Mr2+112MH2

对于所求得的J₂，我们可以理解为在 x轴方向所产生的力矩对应的过鼓心的转动惯量，即同心鼓被视作在为绕该轴进行的定轴旋转，同理，对于 Jy 和y轴也是如此，我们可以做出如下侧视图与立体图加以理解。

随后，我们对于已经列出的二次微分方程组，利用MATLAB 中 ode45 命令对上述微分方程组进行求解，由前面所论述的在小倾角的情况下，我可以将空间的任意转动看成由两个相对独立的转角α和β表示。故上述方程式并不耦合，可以分别求解也可以一起求解。

引入向量的方法，可以不用分析由于九种工况带来的空间几何上分析的复杂度，仅有有向向量完成空间方向的描述，进而大大减少了运算量，简化模型的同时，也让该模型算法可以应对多种发力时机与用力大小的组合情况。将九种情况统一求解后得到以下表格的结果。

 Table 1： 9种不同的情况下， 0.1s时鼓面的倾斜角度一览表

序号	用力参数	1	2	3	4	5	6	7	8	鼓面倾角 (度)
1	发力时机用力大小	0 90	0 80	0 80	0 80	0 80	0 80	0 80	0 80	0.2124°
2	发力时机用力大小	0 90	0 90	0 80	0 80	0 80	0 80	0 80	0 80	0.3855°
3	发力时机用力大小	0 90	0 80	0 80	0 90	0 80	0 80	0 80	0 80	0.1649°
4	发力时机用力大小	-0.1 80	0 80	0 80	0 80	0 80	0 80	0 80	0 80	0.2247°
5	发力时机用力大小	-0.1 80	-0.1 80	0 80	0 80	0 80	0 80	0 80	0 80	0.4094°
6	发力时机用力大小	-0.1 80	0 80	0 80	-0.1 80	0 80	0 80	0 80	0 80	0.1738°
7	发力时机用力大小	-0.1 90	0 80	0 80	0 80	0 80	0 80	0 80	0 80	0.2157°
8	发力时机用力大小	0 90	-0.1 80	0 80	0 90	-0.1 80	0 80	0 80	0 80	0.1354°
9	发力时机用力大小	0 90	0 80	0 80	0 90	-0.1 80	0 80	0 80	-0.1 80	0.3356°

5.3  求解结果分析

为了方便论述与理解，我们对8 根绳子按逆时针顺序进行标号，以1号绳为平面坐标系的 X 轴，其具体示意图如下所示：

基于上图所规定的各物理量的正方向，利用 MATLAB的作图工具，我们可以得到 9种情况下同心鼓的各运动参数(x，α，β)随时间的变化曲线，如下所示：

 Figure9： 9 种情况下的同心鼓各运动参数随时间的变化示意图

代码

 function

[ alpha   best, beta   best, gamma   best,FO  best,FF  best, tt   best]= find   opt(F00,FFO,tt0,T

k,TO, Lf, alpha)

 alpha   final= asin( sin(0.5* pi/180)* cos(12* pi/180))*180/ pi;

 beta   final= asin( sin(0.5* pi/180)* sin(12* pi/180))*180/ pi;

 gamma   final= asin( sqrt( sin( beta   final* pi/180)^2+ sin( alpha   final* pi/180)^2))*180/p

i;

 if nargin==3

 Tk=1000;

TO=1;

 Lf=2000;

 alpha=0.5;

 end

FO  current=F00;

FF  current=FFO;

 tt   current=tt0;

FO  best=FOO;

FF  best=FFO;

 tt   best=tt0;

 pur   current= optfun([F00 FF0 tt0]);

 pur   best= optfun([F00 FF0 tt0]);

 result=[];

 while Tk>TO

 for i=1: Lf

F0=disturbF0(F0  current);

FF=disturbFF(FF  current);

 tt= disturbtt( tt   current);

 fprintf('已执行第%d次扰动 \n',i);

X=[F0 FF tt];

 pur= optfun(X);

 if pur< pur   current

 pur   current= pur;

F0  current=F0;

FF  current=FF;

 tt   current= tt;

 if pur< pur   best

 pur   best= pur;

F0  best=F0;

FF  best=FF;

 tt   best= tt;

 end

 else

 if rand(1,1)< exp(( pur   current- pur)/ Tk)

 pur   current= pur;

F0  current=F0;

FF  current=FF;

 tt   current= tt;

 else

 end

 end

 fprintf('扰动之后的目标函数参数为%1.4f \n', pur);

 fprintf(当前最优的目标函数参数为%1.4f \n ’, pur   current);

 end

 Tk= alpha* Tk;

 fprintf('完成一次~此时目的函数参数为%1.4f \n', pur   best);

 result=[ result, pur   best];

 end

[~, alpha   best, beta   best, gamma   best]= gamma  F t(FO  best, FF  best, tt   best);

 fprintf(经过退火算法优化后,同心鼓与水平面的最优夹角 gamma   best

为: %4.4f° \n', gamma   best);

 fprintf('同心鼓与水平面的目标的夹角 gamma   final为: %4.4f°\n ', gamma   final);

 figure; hold on; grid on;

 plot( result,’r—’,’LineWidth’,2);

 xlabel('{zɪkj..火迭代次数n ',' Fontsize',18,' Fontname', ' Times');

 ylabel('目标::数指标值',' Fontsize',18,' Fontname',' Times');

 title(模拟退火算法搜索全局最优分配方案’,’ Fontsize,20,’ Fontname’,’ Times’);

 end

 function F0  dis=disturbF0(F0)

F0  dis=F0+0.5* randn();

 end

 function FF  dis=disturbFF(FF)

 for i=1:10

FF(i)=FF(i)+0.5* randn();

 end

ind1=0; ind2=0;

 while(ind1==ind2)

ind1= ceil( rand.*10);

ind2= ceil( rand.*10);

 end

FF  dis=FF;

FF  dis(ind1)=FF(ind2);

FF  dis(ind2)=FF(ind1);

 end

 function tt   dis= disturbtt( tt)

 for i=1:10

 tt(i)= tt(i)+0.001* randn();

 end

 tt( tt>0)=0;

ind1=0; ind2=0;

 while(ind1==ind2)

ind1= ceil( rand.*10);

ind2= ceil( rand.*10);

 end

 tt   dis= tt;

 tt   dis(ind1)= tt(ind2);

 tt   dis(ind2)= tt(ind1);

 end