龙门吊问题

这次出的题是一个龙门吊运动过程的角度优化问题, 主要是多目标规划和微分方程, 感觉这个多目标规划的目标函数很难找到一个满意的, 几个目标的无量纲很难找到一个客观的评判标准, 一些参数不能自己主观设置, 但是又没有客观的方法, 确实很麻烦.

然后貌似是陈建业从网上找来的, 直接找到了原题= =不过没有结果, 然而我google搜了一下gantry crane….居然找到了mathworks的一个mathematical modeling的包….里面就包含了这个模型的仿真结果, 不过觉得第一次训练还是好好写写主要不是结果还是要磨合一下队伍就没有用那个结果, 也没和队友说. 自己重新做了一个..下面是用matlab做的一个结果, 感觉还行.

首先是一个随机生成的结果



这里用的是随机优化, 就是个random search….用的是matlab的fmincon..感觉这种题基本很难找到全局最优解, 每一步都有一个微分方程要解, 估计不可导, 就是这个原因弄的我第一天弄了一下午的时间改mathematica的NMinimize…最后还是自己写了个模拟退火…反而可以了..无语

解决过程

第一天和队友讨论了一下基本目标就确定了, 其实要不是这个目标函数的不可导性导致的NMinimize的不稳定性感觉这个题目第一天就基本做完了, 我在优化求解那里花了太多时间了, 后面的时间也基本花在了GUI上, 基本也没有和写作沟通, 这点以后训练要改正, 要把自己的想法和写作说清楚.

最后这个题目其实挺简单的, 不想多讲了, 花了这么多时间也算是一个教训了, 下面直接把代码贴上来吧..= = 这次倒是从华中赛吸取了教训, 把代码都封装了, 写的比较规范, 最后代码都留着在, 感觉REPL还是不要用的太多了, 最好程序还是用编辑器写好了在编译运行, 重复利用率高一些, 之后模型检验也好用.

好了不多说..

第一题角度优化

Clear["Global`*"]
max1[pts_] := Module[{
   time = 100, g = 9.8, m = 6000, r = 15},
  piece[t_] := 
   Piecewise[{{1/2, 0 <= t="" <="" pts[[1]]},="" {-1="" 2,="">= (pts[[1]] + pts[[2]]) && 
       t < (pts[[1]] + pts[[2]] + pts[[3]])}}];
  eq = g*m*r*Sin[\[Theta][t]] + m*r*Cos[\[Theta][t]]*piece[t] + 
    m*r^2*\[Theta]''[t];
  sol = NDSolve[{eq == 0, \[Theta][0] == 0, \[Theta]'[0] == 
      0}, \[Theta][t], {t, 0, time}];
  intep = 
   Cases[Plot[(\[Theta][t] /. sol)*180/Pi, {t, 0, time}], _Line, 
     Infinity][[1, 1]];
  Max[#[[2]] & /@ 
    Select[intep, #[[1]] > (pts[[1]] + pts[[2]] + pts[[3]]) &]]
  ]
random[] := 
 Module[{t1, t2, t3},
  t1 = RandomReal[{1, 10}];
  t3 = RandomReal[{1, t1}];
  t2 = (240 - t1^2 - 2*t1*t3 + t3^2)/(2*t1);
  {t1, t2, t3}]
T = 30000; SeedRandom[1];
at = 0.9999;(* 降温系数, 在0.8~1 *)
lowestTemp = .01; acceptOrder = 
 Table[random[], {i, 10}]
acceptDis = Min[max1[#] & /@ acceptOrder]
bestOrder = {5.880219584586618`, 15.724289580055993`, 
   2.12808475083862`};
bestDis = max1[bestOrder];
While[T > lowestTemp,
  neworder = Table[random[], {i, 10}];
  newdistance = Min[max1[#] & /@ neworder];
  delF = newdistance - bestDis;
  If[newdistance == -Infinity, Continue[]];
  If[
   delF < 0 ,
   bestOrder = neworder[[Ordering[max1[#] & /@ neworder, 1][[1]]]];
   bestDis = newdistance;
   Print[T];
   Print[bestDis];
   Print[bestOrder];
   ];
  T = T *at;
  ];

结果是0.27度

第二题约束优化

  Clear["Global`*"]
max1[pts_] := Module[{
   time = 100, g = 9.8, m = 6000, r = 15},
  piece[t_] := 
   Piecewise[{{1/2, 0 <= t="" <="" pts[[1]]},="" {-1="" 2,="">= (pts[[1]] + pts[[2]]) && 
       t < (pts[[1]] + pts[[2]] + pts[[3]])}}];
  eq = g*m*r*Sin[\[Theta][t]] + m*r*Cos[\[Theta][t]]*piece[t] + 
    m*r^2*\[Theta]''[t];
  sol = NDSolve[{eq == 0, \[Theta][0] == 0, \[Theta]'[0] == 
      0}, \[Theta][t], {t, 0, time}];
  intep = 
   Cases[Plot[(\[Theta][t] /. sol)*180/Pi, {t, 0, time}], _Line, 
     Infinity][[1, 1]];
  Max[#[[2]] & /@ 
    Select[intep, #[[1]] > (pts[[1]] + pts[[2]] + pts[[3]]) &]]
  ]
random[] := 
 Module[{t1, t2, t3},
  t1 = RandomReal[{1, 10}];
  t3 = RandomReal[{1, t1}];
  t2 = (240 - t1^2 - 2*t1*t3 + t3^2)/(2*t1);
  {t1, t2, t3}]
counttime[a_] := Total[a] + 10/((a[[1]] - a[[3]])/2);
count[pts_] := Module[{
    time = 100, g = 9.8, m = 6000, r = 15},
   piece[t_] := 
    Piecewise[{{1/2, 0 <= t="" <="" pts[[1]]},="" {-1="" 2,="">= (pts[[1]] + pts[[2]]) && 
        t < (pts[[1]] + pts[[2]] + pts[[3]])}}];
   eq = g*m*r*Sin[\[Theta][t]] + m*r*Cos[\[Theta][t]]*piece[t] + 
     m*r^2*\[Theta]''[t];
   sol = NDSolve[{eq == 0, \[Theta][0] == 0, \[Theta]'[0] == 
       0}, \[Theta][t], {t, 0, time}];
   intep = 
    Interpolation[
     Cases[Plot[(\[Theta][t] /. sol), {t, 0, time}], _Line, 
       Infinity][[1, 1]]];
   f = Interpolation[
     Cases[Plot[
        Evaluate[D[(\[Theta][t] /. sol), t]], {t, 0, time}], _Line, 
       Infinity][[1, 1]]];
   Abs[pts[[1]] - pts[[3]] - 
     2*r*f[counttime[pts]]*Cos[intep[counttime[pts]]]]];
T = 30000; SeedRandom[28];
at = 0.9999;(* 降温系数, 在0.8~1 *)
lowestTemp = .01; acceptOrder = 
 random[];
bestOrder = acceptOrder
acceptDis = counttime[acceptOrder];
bestDis = counttime[bestOrder]
Quiet[While[T > lowestTemp,
   neworder = random[];
   newdistance = counttime[neworder];
   flag = count[neworder];
   If[flag > 1, Continue[]];
   If[max1[neworder] > 2.0, Continue[]];
   delF = newdistance - bestDis;
   If[
    delF < 0 ,
    bestOrder = neworder;
    bestDis = newdistance;
    Print[bestDis];
    Print[bestOrder];
    ];
   T = T *at;
   ]];

结果忘记了..懒得贴了..

第三题

积分优化

Clear["Global`*"]
r = 15;
g = 9.8;
counttime[a_] := Total[Drop[a, -1]] + 10/((a[[1]] - a[[3]])*a[[4]]);
max1[pts_] := Module[{
   time = 100, g = 9.8, m = 6000, r = 15},
  piece[t_] := 
   Piecewise[{{pts[[4]], 0 <= t="" <="" pts[[1]]},="" {-pts[[4]],="">= (pts[[1]] + pts[[2]]) && 
       t < (pts[[1]] + pts[[2]] + pts[[3]])}}];
  eq = g*m*r*Sin[\[Theta][t]] + m*r*Cos[\[Theta][t]]*piece[t] + 
    m*r^2*\[Theta]''[t];
  sol = NDSolve[{eq == 0, \[Theta][0] == 0, \[Theta]'[0] == 
      0}, \[Theta][t], {t, 0, time}];
  intep = 
   Interpolation[
    Cases[Plot[(\[Theta][t] /. sol)*180/Pi, {t, 0, time}], _Line, 
      Infinity][[1, 1]]];
  NIntegrate[Abs[intep[t]], {t, 0, counttime[pts]}]
  ]
random[] := 
 Module[{t1, t2, t3, a},
  a = RandomReal[];
  t1 = RandomReal[{6, 10}];
  t3 = RandomReal[{4, t1}];
  t2 = (120 - a t1^2 - 2 a t1 t3 + a t3^2)/(2 a t1);
  {t1, t2, t3, a}]
fun1[x_, y_, z_, z1_] := 
 If[(Sin[z] - 0) > 0.001, (x*y - x*r*z1)/Sin[z], 0]
fun2[x_, z_] := x*g*Cos[z]
count[pts_] := Module[{
    time = 100, g = 9.8, r = 15, m = 6000},
   piece[t_] := 
    Piecewise[{{pts[[4]], 0 <= t="" <="" pts[[1]]},="" {-pts[[4]],="">= (pts[[1]] + pts[[2]]) && 
        t < (pts[[1]] + pts[[2]] + pts[[3]])}}];
   eq = g*m*r*Sin[\[Theta][t]] + m*r*Cos[\[Theta][t]]*piece[t] + 
     m*r^2*\[Theta]''[t];
   sol = NDSolve[{eq == 0, \[Theta][0] == 0, \[Theta]'[0] == 
       0}, \[Theta][t], {t, 0, time}];
   intep = 
    Interpolation[
     Cases[Plot[(\[Theta][t] /. sol), {t, 0, time}], _Line, 
       Infinity][[1, 1]]];
   f1p = Cases[
      Plot[Evaluate[D[(\[Theta][t] /. sol), t]], {t, 0, time}], _Line,
       Infinity][[1, 1]];
   f1 = Interpolation[f1p]; 
   f2p = Cases[
      Plot[Evaluate[D[(\[Theta][t] /. sol), {t, 2}]], {t, 0, 
        time}], _Line, Infinity][[1, 1]];
   f2 = Interpolation[f1p];
   x1 = Abs[
     pts[[1]] - pts[[3]] - 
      2*r*f1[counttime[pts]]*Cos[intep[counttime[pts]]]];
   x2 = FindMaximum[{fun1[m, pts[[4]], intep[zz], f2[zz]], 
       0 < zz < pts[[1]] || (pts[[1]] + pts[[2]]) < 
         zz < (pts[[1]] + pts[[2]] + pts[[3]])}, {zz, pts[[1]]/2}][[
     1]];
   x3 = FindMaximum[{fun2[m, intep[zx]],
       pts[[1]] < 
         zx < (pts[[1]] + pts[[2]]) || (pts[[1]] + pts[[2]]) < 
         zx < (pts[[1]] + pts[[2]] + pts[[3]])}, {zx, (pts[[1]] + 
          pts[[2]])/2}][[1]];
   x1 <= 1 && x2 <= 196000 && x1 <= 196000;
   {x1, x2, x3}
   ];
T = 30000; SeedRandom[28];
at = 0.9999;(* 降温系数, 在0.8~1 *)
lowestTemp = .01; acceptOrder = 
 random[];
bestOrder = acceptOrder
acceptDis = max1[acceptOrder];
bestDis = 100
Quiet[While[T > lowestTemp,
   neworder = random[];
   newdistance = max1[neworder];
   flag = count[neworder];
   Print[flag];
   If[flag == False, Continue[]];
   delF = newdistance - bestDis;
   If[
    delF < 0 || 
     RandomVariate[UniformDistribution[]] < Exp[-(delF/T)],
    acceptOrder = neworder;
    acceptDis = newdistance;
    ];
   If[
    delF < 0 ,
    bestOrder = neworder;
    bestDis = newdistance;
    Print[bestDis];
    Print[bestOrder];
    ];
   T = T *at;
   ]];

Mmax求解

Clear["Global`*"]
r = 15;
g = 9.8;
counttime[a_] := Total[Drop[a, -2]] + 10/((a[[1]] - a[[3]])*a[[4]]);
fun1[x_, y_, z_, z1_] := 
 If[(Sin[z] - 0) > 0.001, (x*y - x*r*z1)/Sin[z], 0]
fun2[x_, z_] := x*g*Cos[z]
count[pts_] := Module[{
    time = 100, g = 9.8, r = 15, m = pts[[5]]},
   piece[t_] := 
    Piecewise[{{pts[[4]], 0 <= t="" <="" pts[[1]]},="" {-pts[[4]],="">= (pts[[1]] + pts[[2]]) && 
        t < (pts[[1]] + pts[[2]] + pts[[3]])}}];
   eq = g*m*r*Sin[\[Theta][t]] + m*r*Cos[\[Theta][t]]*piece[t] + 
     m*r^2*\[Theta]''[t];
   sol = NDSolve[{eq == 0, \[Theta][0] == 0, \[Theta]'[0] == 
       0}, \[Theta][t], {t, 0, time}];
   intep = 
    Interpolation[
     Cases[Plot[(\[Theta][t] /. sol), {t, 0, time}], _Line, 
       Infinity][[1, 1]]];
   f1p = Cases[
      Plot[Evaluate[D[(\[Theta][t] /. sol), t]], {t, 0, time}], _Line,
       Infinity][[1, 1]];
   f1 = Interpolation[f1p]; 
   f2p = Cases[
      Plot[Evaluate[D[(\[Theta][t] /. sol), {t, 2}]], {t, 0, 
        time}], _Line, Infinity][[1, 1]];
   f2 = Interpolation[f1p];
   x1 = Abs[
     pts[[1]] - pts[[3]] - 
      2*r*f1[counttime[pts]]*Cos[intep[counttime[pts]]]];
   x2 = FindMaximum[{fun1[m, pts[[4]], intep[zz], f2[zz]], 
       0 < zz < pts[[1]] || (pts[[1]] + pts[[2]]) < 
         zz < (pts[[1]] + pts[[2]] + pts[[3]])}, {zz, pts[[1]]/2}][[
     1]];
   x3 = FindMaximum[{fun2[m, intep[zx]],
       pts[[1]] < 
         zx < (pts[[1]] + pts[[2]]) || (pts[[1]] + pts[[2]]) < 
         zx < (pts[[1]] + pts[[2]] + pts[[3]])}, {zx, (pts[[1]] + 
          pts[[2]])/2}][[1]];
   {x2, x3}
   ];
count[{8.718908223457277`, 5.708759029792792`, 7.1725028248090155`, 
  0.4815493490832745`, 11380}]
NestWhileList[# + 
   0.01 &, 11383.1, (count[{8.718908223457277`, 5.708759029792792`, 
        7.1725028248090155`, 0.4815493490832745`, #}][[1]] - 196000) <
    0 &]

第四题仿真

2D仿真

Clear["Global`*"]
x = {8.718908223457277`, 5.708759029792792`, 7.1725028248090155`, 
   0.4815493490832745`};
counttime[a_] := Total[Drop[a, -1]] + 10/((a[[1]] - a[[3]])*a[[4]]);
c1 = x[[1]];
c2 = x[[2]];
c3 = x[[3]];
a = x[[4]];
time = counttime[x];
aa[t_] := 
  Piecewise[{{a, t < c1}, {-a, t >= (c1 + c2) && t < (c1 + c2 + c3)}}];
int1[t_] := 
  Piecewise[{{t*a, 0 <= t="" <="" c1},="" {c1*a,="">= c1 && t <= (c1="" +="" c2)},="" {c1*a="" -="" (t="" c1="" c2)*a,="" t=""> (c1 + c2) && t <= (c1="" +="" c2="" c3)},="" {c1*a="" -="" c3*a,="" t=""> (c1 + c2 + c3)}}];
int2[t1_] := Integrate[int1[t], {t, 0, t1}];
m = 6000;
r = 15;
g = 9.8;
eq1 = g*m*r*Sin[\[Theta][t]] + m*r*Cos[\[Theta][t]]*aa[t] + 
   m*r^2*\[Theta]''[t];
sol = NDSolve[{eq1 == 0, \[Theta][0] == 0, \[Theta]'[0] == 
     0}, \[Theta][t], {t, 0, time}];
int = Interpolation[
   Cases[Plot[(\[Theta][t] /. sol)*180/Pi, {t, 0, time}, 
      PlotRange -> All], _Line, Infinity][[1, 1]]];
Manipulate[
 With[{\[Theta] = int[t], x = (-13 - 11.5)/2 + int2[t], r = 6},
  Graphics[{
    Style[
     Text[Row[{Style["\[Theta]", Italic], " = ", 
        NumberForm[\[Theta], {3, 2}], "度"}], {5, 15.5}], 12], 
    Style[Text[
      Row[{Style["x", Italic], " = ", 
        NumberForm[x - (-13 - 11.5)/2, {3, 2}], "米"}], {5, 18.5}], 12],
    Text[A, {-10, -25}],
    Text[B, {(-13 - 11.5)/2 + 60, -25}],
    Text[C, {(-13 - 11.5)/2 + 70, -25}],
    Black, Rectangle[{-17, 2}, {80, 1.75}], Black, 
    Rectangle[{x - 0.75, 2}, {x + 0.75, 3.5}], 
    Line[{{x - r*Sin[\[Theta]*Degree]*2, -r*2*Cos[\[Theta]*Degree] + 
        2}, {x, 2}}], {Blue,
     Rotate[
      Rectangle[{x - r*2*Sin[\[Theta]*Degree] - 0.75 - 
         0.3, -r*2*Cos[\[Theta]*Degree] + 0.5 - 0.3}, {x - 
         r*2*Sin[\[Theta]*Degree] + 0.75 + 
         0.3, -r*2*Cos[\[Theta]*Degree] + 2 + 0.3}], -\[Theta]*
       Degree]}, Red, Rectangle[{-16, -20}, {-10, -14}], 
    Rectangle[{(-13 - 11.5)/2 + 60, -20}, {(-13 - 11.5)/2 + 
       100, -14}], Dashed, 
    Line[{{(-13 - 11.5)/2 + 70, -14}, {(-13 - 11.5)/2 + 70, 2}}], 
    Line[{{(-13 - 11.5)/2, -14}, {(-13 - 11.5)/2, 2}}], 
    Line[{{(-13 - 11.5)/2 + 60, -14}, {(-13 - 11.5)/2 + 60, 2}}]}, 
   PlotRange -> Automatic, ImageSize -> {500, 350}]],
 {{t, 0., "time"}, 0., time, 0.01, Appearance -> "Labeled"},
 {{m, 6000, "货物质量"}, 4000, 9000},
 {{M, 6000, "吊车质量"}, 4000, 9000},
 {{\[Mu], .02, "滑动摩擦系数"}, 0.01, 0.2},
 {{t1, 3, "减速时间"}, 1, 4},
 {{t3, 3, "加速时间"}, 1, 4},
 {{a, 0.5, "加速度"}, 0, 1},
 {{D1, 60, "船到码头距离"}, 50, 70},
 {{d2, 10, "码头到终点距离"}, 5, 15},
 {{l, 15, "绳长"}, 10, 20}, TrackedSymbols :> {t, m1, m2, m3}]
```

gui界面如下
![](http:\/\/ooo.0o0.ooo\/2016\/07\/25\/579649c76c0a5.png)


#### 3D仿真
```mathematica
Clear["Global`*"]
x = {8.56607806024076`, 5.485675399190234`, 7.7076876866414015`, 
   0.5};
counttime[a_] := Total[Drop[a, -1]] + 10/((a[[1]] - a[[3]])*a[[4]]);
c1 = x[[1]];
c2 = x[[2]];
c3 = x[[3]];
a = x[[4]];
time = counttime[x];
aa[t_] := 
  Piecewise[{{a, t < c1}, {-a, t >= (c1 + c2) && t < (c1 + c2 + c3)}}];
int1[t_] := 
  Piecewise[{{t*a, 0 <= t="" <="" c1},="" {c1*a,="">= c1 && t <= (c1="" +="" c2)},="" {c1*a="" -="" (t="" c1="" c2)*a,="" t=""> (c1 + c2) && t <= (c1="" +="" c2="" c3)},="" {c1*a="" -="" c3*a,="" t=""> (c1 + c2 + c3)}}];
int2[t1_] := Integrate[int1[t], {t, 0, t1}];
m = 6000;
r = 15;
g = 9.8;
eq1 = g*m*r*Sin[\[Theta][t]] + m*r*Cos[\[Theta][t]]*aa[t] + 
   m*r^2*\[Theta]''[t];
sol = NDSolve[{eq1 == 0, \[Theta][0] == 0, \[Theta]'[0] == 
     0}, \[Theta][t], {t, 0, time}];
int = Interpolation[
   Cases[Plot[(\[Theta][t] /. sol)*180/Pi, {t, 0, time}, 
      PlotRange -> All], _Line, Infinity][[1, 1]]];
Manipulate[
 Show[
  Graphics3D[Cylinder[{{-18, 0 + bf, 9}, {80, 0 + bf, 9}}, 1]],
  Graphics3D[Cylinder[{{-18, -10, 9}, {-18, 10, 9}}, 1]],
  Graphics3D[Cylinder[{{80, -10, 9}, {80, 10, 9}}, 1]],
  With[{\[Theta] = int[t], x = (-13 - 11.5)/2 + int2[t], r = 6},
   Graphics3D[{
     Style[
      Text[Row[{Style["\[Theta]", Italic], " = ", 
         NumberForm[\[Theta], {3, 2}], "度"}], {-5, 0, -15.5}], 12], 
     Style[Text[
       Row[{Style["x", Italic], " = ", 
         NumberForm[x - (-13 - 11.5)/2, {3, 2}], "米"}], {-5, 
        0, -18.5}], 12],
     Cuboid[{-2 + b + x - (-13 - 11.5)/2, -1 + bf, 
       7}, {2 + b + x - (-13 - 11.5)/2, 1 + bf, 9}], 
     Cylinder[{{b + x - (-13 - 11.5)/2 - 2*r*Sin[\[Theta]*Degree], 
        0 + bf, -r*2*Cos[\[Theta]*Degree] + 0.75}, {b + 
         x - (-13 - 11.5)/2, bf, 7}}, 0.18],
     Cuboid[{x - r*2*Sin[\[Theta]*Degree] - 0.75 + 
        b - (-13 - 11.5)/2, -1 + bf, -r*2*Cos[\[Theta]*Degree] + 
        0.5}, {x - r*2*Sin[\[Theta]*Degree] + 0.75 + 
        b - (-13 - 11.5)/2, 1 + bf, -r*2*Cos[\[Theta]*Degree] + 2}]}, 
    PlotRange -> Automatic, ImageSize -> {500, 350}]]],
 "龙门吊前后移动",
 {{bf, 0, ""}, -15, 9},
 "龙门吊左右移动",
 {{b, (-13 - 11.5)/2, ""}, -15, 8},
 "绳长",
 {{ud, -5, ""}, -7, 5},
 {{m, 6000, "货物质量"}, 4000, 9000},
 {{M, 6000, "吊车质量"}, 4000, 9000},
 {{\[Mu], .02, "滑动摩擦系数"}, 0.01, 0.2},
 {{t1, 3, "减速时间"}, 1, 4},
 {{t3, 3, "加速时间"}, 1, 4},
 {{a, 0.5, "加速度"}, 0, 1},
 {{D1, 60, "船到码头距离"}, 50, 70},
 {{d2, 10, "码头到终点距离"}, 5, 15},
 {{t, 0., "time"}, 0., time, 0.01, Appearance -> "Labeled"},
 ControlPlacement -> Left]
```

3D gui如下
![](http:\/\/ooo.0o0.ooo\/2016\/07\/25\/579649c635755.png)

### 附录
最后是整体的一个测试函数, 为了验证不同输入的验证程序....显示了加速度,速度,位移和角度还有一个简单的GUI实现..

```mathematica
Clear["Global`*"]
x = {3.5052860093210825`, 30.743884624269622`, 3.1792840069469843`, 
   0.5};
counttime[a_] := Total[Drop[a, -1]] + 10/((a[[1]] - a[[3]])*a[[4]]);
c1 = x[[1]];
c2 = x[[2]];
c3 = x[[3]];
a = x[[4]];
time = counttime[x];
aa[t_] := 
  Piecewise[{{a, t < c1}, {-a, t >= (c1 + c2) && t < (c1 + c2 + c3)}}];
int1[t_] := 
  Piecewise[{{t*a, 0 <= t="" <="" c1},="" {c1*a,="">= c1 && t <= (c1="" +="" c2)},="" {c1*a="" -="" (t="" c1="" c2)*a,="" t=""> (c1 + c2) && t <= (c1="" +="" c2="" c3)},="" {c1*a="" -="" c3*a,="" t=""> (c1 + c2 + c3)}}];
int2[t1_] := Integrate[int1[t], {t, 0, t1}];
m = 6000;
r = 15;
g = 9.8;
eq1 = g*m*r*Sin[\[Theta][t]] + m*r*Cos[\[Theta][t]]*aa[t] + 
   m*r^2*\[Theta]''[t];
sol = NDSolve[{eq1 == 0, \[Theta][0] == 0, \[Theta]'[0] == 
     0}, \[Theta][t], {t, 0, time}];
Plot[aa[t], {t, 0, time}, PlotLabel -> "a-time"]
Plot[int1[t], {t, 0, time}, PlotLabel -> "u-time"]
Plot[int2[t3], {t3, 0, time}, PlotLabel -> "x-time"]
intep = Cases[
    Plot[(\[Theta][t] /. sol)*180/Pi, {t, 0, time}, 
     PlotRange -> All], _Line, Infinity][[1, 1]];
int = Interpolation[intep];
Print["t1的时间为: ", c1]
Print["t2的时间为: ", c2]
Print["t3的时间为: ", c3]
Print["t4的时间为: ", time - Total[Drop[x, -1]]]
Print["a为: ", x[[4]]]
Print["总的时间为: ", time]
Print["B到C的最大角度为", 
  Max[#[[2]] & /@ Select[intep, #[[1]] > (c1 + c2 + c3) &]]];
Plot[(\[Theta][t] /. sol)*180/Pi, {t, 0, time}, PlotRange -> All, 
 PlotLabel -> "\[Theta]-time"]
Manipulate[
 With[{\[Theta] = int[t], x = (-13 - 11.5)/2 + int2[t], r = 6},
  Graphics[{
    Style[
     Text[Row[{Style["\[Theta]", Italic], " = ", 
        NumberForm[\[Theta], {3, 2}], "度"}], {-5, -15.5}], 12], 
    Style[Text[
      Row[{Style["x", Italic], " = ", 
        NumberForm[x - (-13 - 11.5)/2, {3, 2}], "米"}], {-5, -18.5}], 
     12],
    Black, Rectangle[{-15, 2}, {80, 1.75}], Black, 
    Rectangle[{x - 0.75, 2}, {x + 0.75, 3.5}], 
    Line[{{x - r*Sin[\[Theta]*Degree]*2, -r*2*Cos[\[Theta]*Degree] + 
        2}, {x, 2}}],
    Rotate[
     Rectangle[{x - r*2*Sin[\[Theta]*Degree] - 
        0.75, -r*2*Cos[\[Theta]*Degree] + 0.5}, {x - 
        r*2*Sin[\[Theta]*Degree] + 0.75, -r*2*Cos[\[Theta]*Degree] + 
        2}], -\[Theta]*Degree]}, PlotRange -> Automatic, 
   ImageSize -> {500, 350}]],
 {{t, 0., "time"}, 0., time, 0.01, Appearance -> "Labeled"},
 Style["\nmasses", Bold], TrackedSymbols :> {t, m1, m2, m3}]

下面是角度和gui界面, 还有一些测试用的图像就不截图了…


最后…

昨天算是结束了数模第一阶段的培训, 感觉还是学了点东西, 这次训练感觉和华中赛时候犯的错误类似, 三个队友间的沟通太少, 自己只顾着编程了, 没有去考虑队友那边的进度, 三个人的分工太明确, 导致沟通太少, 每个人对对方做的工作都不太了解, 最后导致的就是论文问题非常的大, 基本没有把我程序的意思, 建模的思路写进去, 这点本来在华中赛就显示的非常严重了, 但是这次还是出现了.

不过考虑到这次还是我们组第一次开始, 还有5次训练, 之后一定要吸取教训了. 明天打算回家去了, 本来今天和王衍田和王浩然准备这两天去做一下毛概的调查报告, 今天下午过去把我热死了= =感觉今天也没做啥, 下午讨论了一下觉得王浩然那个题目主题太不明确了, 最后做的效果可能不太好. 所以还是先回去算了, 实践的话在自己家那边找个养老院或者社区做一下调查好了, 感觉比大学生调查之类的靠谱一些.

回家后最好能把借的那么多书看一下…现在倒是觉得R比python都好用了..两个都学下吧, julia就算了, 感觉没啥用…还不如用matlab, 把simulink用熟点感觉大三上课做作业之类的还是有点用的…

最后大三基本就是最后一年了, 保持现在这个状态吧, 不管最后结果怎么样, 还是要把自己该做的做了, 至于保研还是出国….感觉还是不太确定, 托福感觉口语太差了..两次都跪在口语了orz不知道怎么练习了, 不管怎么样最好能在下个学期把GRE考完了, 然后在寒假结束前结束托福, 之后安心刷绩点, 大三刷下题, 最好大三暑假找个实习或者去哪个学校当个RA, 经历也好用来申请之类的…谷歌那个夏令营感觉可以试试, 不过还不知道选什么主题哎..我现在算法感觉都没刷太多, 本来老是想着暑假刷, 结果到了暑假根本没动力啊…毕竟上学的时候每天上课, 现在放假了根本没动力啊摔!…科三挂了两次…过两天再考一次! 一定要过了…免得又浪费时间..醉了= =

恩, 总结完毕..睡觉!