问题三

最小螺距单目标优化模型的建立

问题二考虑了在螺距固定的条件下计算舞龙队盘入的终止时间，问题三在第二问的基础提出了改变螺距的要求，即求解在螺距最小为多少时，龙头前把手能够沿着相应的螺线盘入到调头空间的边界。故可将其转换为一个单目标优化问题，目标函数为求解最小的螺距，将板凳龙模型中具体物理要求与龙头前把手能够沿着相应的螺线盘入到调头空间的边界作为约束条件，由于该目标函数非线性且较为复杂，故本文采取了带罚函数的双重模拟退火算法与变步长搜索算法结合进行求解。

目标函数与约束条件

由于此时螺距的取值受到约束条件的限制，故将螺距$a$的表达修改为变量$A$，可写出目标函数的表达为：
$$\min A$$
其中约束条件如下如下：

（1）. 根据题意，龙头能够进入调头空间的边界，即此时结束位置的极径应该小于4.5m，此时螺线方程为$r=\frac{A}{2\pi }\cdot \theta$，故可列约束条件为：
$$\frac{A}{2\pi }\cdot \theta \le 4.5(m)$$
（2）.根据题意，在龙头能够进入调头空间前，板凳之间不能相撞，即此时的决策变量应该满足的约束条件为：
$$\sum _{i=1}^{223}\sum _{j=1,}^{223}{\delta }_{ij}=0$$
（3）.板凳的宽度为15cm，螺距的宽度能够支持两个板凳并肩移动，即此时螺距应该满足的约束条件为：
$$A\ge 30(cm)$$
（4）.调头空间意味着舞龙队的总长度不能长于螺线的总弧长，此时应该满足的约束条件为：
$$A{\int }_0^{32\pi }\sqrt{1+{\theta }^2}d\theta \ge 220\times 222+341$$
（5）.此时碰撞时间的求解模型基于问题一中对位置的迭代求解模型，给出问题一中的迭代方程为：
$$\{\begin{array}{ll}(x_{i+1}-x_i{)}^2+(y_{i+1}-y_i{)}^2=165^2& \\ x_{i+1}=\frac{55}{2\pi }\cdot {\theta }_{i+1}\cdot \cos ({\theta }_{i+1})& \\ y_{i+1}=\frac{55}{2\pi }\cdot {\theta }_{i+1}\cdot \sin ({\theta }_{i+1})& \end{array}$$
此处$(x_{i+1},y_{i+1})$为下一龙身前孔的坐标，$(x_i,y_i)$为上一龙身前孔的坐标，${\theta }_{i+1}$为下一龙身前孔此时的弧度。

综上所述，给出问题三最小螺距单目标优化模型如下：
$$\begin{gathered} \min A \\ \text{s.t.}\{\begin{array}{ll}\frac{A}{2\pi }\cdot \theta \le 4.5(m)& \\ {\sum }_{i=1}^{223}{\sum }_{j=1,}^{223}{\delta }_{ij}=0& \\ A\ge 30(cm)& \\ A{\int }_0^{32\pi }\sqrt{1+{\theta }^2}d\theta \ge 220\times 222+341& \\ (x_{i+1}-x_i{)}^2+(y_{i+1}-y_i{)}^2=165^2& \\ x_{i+1}=\frac{55}{2\pi }\cdot {\theta }_{i+1}\cdot \cos ({\theta }_{i+1})& \\ y_{i+1}=\frac{55}{2\pi }\cdot {\theta }_{i+1}\cdot \sin ({\theta }_{i+1})& \end{array} \end{gathered}$$

模型的求解

本题的方程比较复杂且非线性，但是最终结果为一维且搜索空间较小，于是本文在求解的时候采取了带罚函数的模拟退火算法与定步长遍历算法结合的办法，启发式算法可以很好的在解空间内探查到最优解，而变步长遍历算法则可以在解空间较小的情况下遍历空间得到整个函数的分布情况进而查找到最优解。下面是两种算法的求解过程：

带罚函数的双重模拟退火算法

双重模拟退火算法(Dual Annealing)是一种启发式算法，适合用于求解复杂的非线性优化问题。在该问题中，由于目标函数和约束条件都比较复杂，使用双重模拟退火算法可以有效搜索到接近全局最优解的结果。为了处理约束条件，本文引入罚函数，将约束条件违背的程度引入到目标函数中，确保求解过程中尽量满足约束条件。

Step1 设定初始解，初始温度和罚函数

设定初始螺距$A_0=30cm$，初始温度$T_0=100$，设定罚函数因子${\lambda }_1$，${\lambda }_2$与${\lambda }_3$，其大小为某远大于目标值的数，定义罚函数为：
$$P(A)={\lambda }_1\cdot \max (0,\frac{A}{2\pi }\cdot \theta -4.5)+{\lambda }_2\cdot \left(\sum _{i=1}^{223}\sum _{j=1,}^{223}{\delta }_{ij}\right)+{\lambda }_3\cdot \max (0,220\times 222+341-A{\int }_0^{32\pi }\sqrt{1+{\theta }^2}d\theta )$$
由此实现对约束条件的处理

Step2 全局搜索

在当前解$x_i$的邻域内生成新解$x_{i+1}$，计算新解的函数值$f(x_{i+1})$，如果$f(x_{i+1})<f(x_i)$，则接受新解作为当前解。若不满足该条件，则以概率$P=\exp \left(\frac{f(x_i)-f(x_{i+1})}{T}\right)$接受新解，随着温度的降低，接受劣解的概率逐渐减小。

Step3 局部搜索

在接受新解后，使用局部优化算法（如最速下降法或梯度下降法）对当前解进行局部优化，以快速收敛到邻域内的局部最优解。

Step4 温度下降

更新温度 $T_{i+1}=\alpha \cdot T_i$,逐步降低温度。

Step5 搜索结束，输出结果

当最终结果收敛时，将搜索得到的结果输出，获得此时的螺距长度。

变步长遍历算法

变步长求解是一种结合了大步长的快速搜索和小步长的精细搜索的优化策略。在开始阶段使用较大的步长快速缩小解的范围，后面逐步缩小步长，通过类似二分查找的方式，精确找到最优解。

Step1 设定初始解和初始步长

设定初始螺距$A_0=30cm$，最大螺距为在通过观察模拟退火算法求解过程基础上确定的边界解$A_{max}=50cm$，初始步长为$\Delta A=2cm$

Step2 大步长快速搜索

从初始螺距$A_0$开始搜索，以$\Delta A=2cm$为步长开始遍历，快速搜索解空间直达搜索到$A_{max}$，计算每个搜索点的目标函数值，若当前螺距下的目标函数值优于前面搜索到的最优值，则更新最优值，重复上述过程，直到搜索完初始步长的范围。

Step3 二分精细搜索

经过大步长搜索找到最优螺距$A_{best}$，在其左右邻域$[A_{best}-\Delta A,A_{best}+\Delta A]$内进行搜索，此时设定步长为$\frac{\Delta A}{2}$。重复Step2与Step3的步骤。

Step4 导出搜索结果

当搜索步长缩小到$\frac{\Delta A}{1024}$时，其值以达到误差容限内，输出此时的目标函数值$A_{best}$。

问题四

模型的建立

问题四设置盘入螺线的螺距为1.7m，盘入螺线与盘出螺线为中心对称关系，并且舞龙队需要在调头空间内实现调头，其中调头位置在盘入螺线上，调出位置在盘出螺线，即现在确定了调头空间的调入点与调出点，以及调入圆弧与调出圆弧之间的数量关系，在未知圆弧进入点与圆弧出去点的情况下对圆弧曲线的长度进行优化，求解出最短的调头曲线。为了求解这个问题，首先需要求解调头螺线的数学方程，进而通过几何关系建立起调头圆弧的数学模型，根据圆弧的弧度和调头的位置确定整个调头曲线长度的表达式，作为模型的目标函数，然后通过数学模型以及实际板凳的物理模型给出目标函数的约束条件，

调头弧线模型的建立

首先已知盘入螺线的螺距为 1.7 m,盘出螺线与盘入螺线关于螺线中心呈中心对称，故可先求解出盘入螺线$C_1$与盘出螺线$C_2$的表达式分别为：
$$\{\begin{array}{l}C1:r=\frac{170}{2\pi }\cdot \theta \\ C2:r=-\frac{170}{2\pi }\cdot \theta \end{array}$$
由于目前只可推出进入调头空间的点与离开调头空间的点，由于两点分别位于中心对称的两螺线与调头区域圆的交点上，所以两点坐标存在中心对称关系,可设进入点为$A(x_A,y_A)$，离开点为$B(x_B,y_B)$，则两点满足的方程为：
$$\{\begin{array}{l}{x}_A+x_B=0\\ y_A+y_B=0\\ x^2+y^2=450^2\end{array}$$
由于开始调头的点在盘入螺线和前段圆弧上，且经过该点后龙头前把手运动方程从螺线方程变为圆方程，结束调头的点则相反，故可设开始调头的点为点$C$，坐标为$(x_C,y_C)$，结束调头的点为点$D$，坐标为$(x_D,y_D)$。其示意图如下图所示：

故可表示调头弧线的总长度$L$可表示为：
$$L=\stackrel{⌢}{\mathrm{AC}}+\stackrel{⌢}{\mathrm{CD}}+\stackrel{⌢}{\mathrm{DB}}$$
为确定两圆的圆心位置，需计算出点$C$处螺线的导数，将其带入曲率半径的公式即可计算出前段圆弧的半径，进而确定后段圆弧半径，即可确定两圆心的位置以及调头弧线的长度。

由于点$C$位于螺线上，故可先求解螺线的$\frac{dx}{d\theta }$与$\frac{dy}{d\theta }$，再通过链式法则即可求解得到$\frac{dy}{dx}$，计算过程如下：
$$\begin{gathered} \frac{dx}{d\theta }=\frac{170}{2\pi }\cdot (\cos (\theta )-\theta \cdot \sin (\theta )) \\ \frac{dy}{d\theta }=\frac{170}{2\pi }\cdot (\sin (\theta )+\theta \cdot \cos (\theta )) \\ \frac{dy}{dx}=\frac{\frac{dy}{d\theta }}{\frac{dx}{d\theta }}=\frac{\sin (\theta )+\theta \cdot \cos (\theta )}{\cos (\theta )-\theta \cdot \sin (\theta )} \end{gathered}$$
故点$C$处的导数为：
$$y_C^{\prime }(x)=\frac{d{y}_C}{d{x}_C}=\frac{\frac{d{y}_C}{d{\theta }_C}}{\frac{d{x}_C}{d{\theta }_C}}=\frac{\sin ({\theta }_C)+{\theta }_C\cdot \cos ({\theta }_C)}{\cos ({\theta }_C)-\theta \cdot \sin ({\theta }_C)}$$
根据互相垂直的直线之间存在$k_1\cdot k_2=-1$可知$C$点切线的垂线斜率$k_{w1}$为：
$$k_{w1}=\frac{\frac{d{y}_C}{d{\theta }_C}}{\frac{d{x}_C}{d{\theta }_C}}=-\frac{\cos ({\theta }_C)-\theta \cdot \sin ({\theta }_C))}{\sin ({\theta }_C)+{\theta }_C\cdot \cos ({\theta }_C}$$
同理，$D$点切线的垂线斜率$k_{w2}$为：
$$k_{w2}=\frac{\frac{d{y}_D}{d{\theta }_D}}{\frac{d{x}_D}{d{\theta }_D}}=-\frac{\cos ({\theta }_D)-\theta \cdot \sin ({\theta }_D))}{\sin ({\theta }_D)+{\theta }_D\cdot \cos ({\theta }_D)}$$
作出切线图后，调头空间中存在如图的几何约束，示意图如下图所示：

其中$r_1$为前段圆的半径，$r_2$为后段圆的半径，设两圆圆心与坐标分别为$O_1(x_{O1},y_{O1})$与$O_2(x_{O2},y_{O2})$，此时可依图列出方程组为：
$$\{\begin{array}{ll}{x}_{O1}=x_C+\frac{1}{\sqrt{k_{w1}^2+1}}\cdot r_1& \\ y_{O1}=y_C+\frac{k_{w1}}{\sqrt{k_{w1}^2+1}}\cdot r_1& \\ x_{O2}=x_D+\frac{1}{\sqrt{k_{w2}^2+1}}\cdot r_2& \\ y_{O2}=y_D+\frac{k_{w2}}{\sqrt{k_{w2}^2+1}}\cdot r_2& \\ r_1=2r_2& \\ \sqrt{(x_{O1}-x_{O2}{)}^2+(y_{O1}-y_{O2}{)}^2}=r_1+r_2& \end{array}$$
求解上述方程组，即可计算得到两圆圆心的坐标。

得到圆心坐标后，则需要再计算出两圆的切点$E$的坐标，由于已知两圆的圆心坐标和圆心到$E$点的直线的斜率，故需要通过简单的比例计算就可得到$E$点的坐标如下：
$$\{\begin{array}{l}{x}_E=x_{O1}+\frac{r1}{r1+r2}\cdot (x_{O2}-x_{O1})\\ y_E=y_{O1}+\frac{r1}{r1+r2}\cdot (y_{O2}-y_{O1})\end{array}$$
下面依据$E$点的坐标，注意到前段圆弧均为顺时针方向，后段圆弧为逆时针方向，在导出公式中将角度的顺序进行调整，可导出两圆转过的弧度${\varphi }_1$与${\varphi }_2$为：
$$\{\begin{array}{l}{\varphi }_1=\arctan (\frac{y_C-y_{O1}}{x_C-x_{O1}})-\arctan (\frac{y_E-y_{O1}}{x_E-x_{O1}})\\ {\varphi }_2=\arctan (\frac{y_D-y_{O2}}{x_D-x_{O2}})-\arctan (\frac{y_E-y_{O2}}{x_E-x_{O2}})\end{array}$$
即可求解出$L$的表达式为：
$$L=\frac{170}{2\pi }{\int }_{{\theta }_C}^{{\theta }_A}\sqrt{1+{\theta }^2}d\theta +r_1\cdot {\varphi }_1+r_2\cdot {\varphi }_2+\frac{170}{2\pi }{\int }_{{\theta }_D}^{{\theta }_B}\sqrt{1+{\theta }^2}d\theta$$

建立最短弦长单目标优化模型

前文已经求解出$L$的表达式，根据题意，需要通过调整圆弧，保持各部分相切在条件下使得调头曲线变短，故可基于此构建最短调头曲线的优化模型，其目标函数为最短调头曲线弧长，约束条件则基于调头区域内需要保持各部分相切以及圆弧之间的几何关系约束。

目标函数与约束条件

目标函数为最短调头曲线弧长，由于前文的推导过程均基于$C$点与$D$点的不同坐标，即${\theta }_C与{\theta }_D$，故可视$L$为${\theta }_C$与${\theta }_D$的函数，故目标函数数学表达式为：
$$\min L({\theta }_C,{\theta }_D)=\frac{170}{2\pi }{\int }_{{\theta }_C}^{{\theta }_A}\sqrt{1+{\theta }^2}d\theta +r_1\cdot {\varphi }_1+r_2\cdot {\varphi }_2+\frac{170}{2\pi }{\int }_{{\theta }_D}^{{\theta }_B}\sqrt{1+{\theta }^2}d\theta$$
其中$r_1,r_2,{\varphi }_1,{\varphi }_2$均可由${\theta }_C$与${\theta }_D$确定，${\theta }_A$与${\theta }_B$为定值。

约束条件在前文调头弧线模型的建立已经进行了详细的推导，得出了方程组(44)，(45)，(46)，（47）与(48)，方程组本问中需要的约束条件。

综上所述，最短弦长单目标优化模型建立如下：
$$\begin{gathered} \underset{{\theta }_C,{\theta }_D}{\min }L({\theta }_C,{\theta }_D)=\frac{170}{2\pi }{\int }_{{\theta }_C}^{{\theta }_A}\sqrt{1+{\theta }^2}d\theta +r_1\cdot {\varphi }_1+r_2\cdot {\varphi }_2+\frac{170}{2\pi }{\int }_{{\theta }_D}^{{\theta }_B}\sqrt{1+{\theta }^2}d\theta \\ \text{s.t.}\{\begin{array}{ll}{k}_{w1}=-\frac{\cos ({\theta }_C)-\theta \cdot \sin ({\theta }_C))}{\sin ({\theta }_C)+{\theta }_C\cdot \cos ({\theta }_C}& \\ k_{w2}=-\frac{\cos ({\theta }_D)-\theta \cdot \sin ({\theta }_D))}{\sin ({\theta }_D)+{\theta }_D\cdot \cos ({\theta }_D)}& \\ x_{O1}=x_C+\frac{1}{\sqrt{k_{w1}^2+1}}\cdot r_1& \\ y_{O1}=y_C+\frac{k_{w1}}{\sqrt{k_{w1}^2+1}}\cdot r_1& \\ x_{O2}=x_D+\frac{1}{\sqrt{k_{w2}^2+1}}\cdot r_2& \\ y_{O2}=y_D+\frac{k_{w2}}{\sqrt{k_{w2}^2+1}}\cdot r_2& \\ r_1=2r_2& \\ \sqrt{(x_{O1}-x_{O2}{)}^2+(y_{O1}-y_{O2}{)}^2}=r_1+r_2& \\ x_E=x_{O1}+\frac{r1}{r1+r2}\cdot (x_{O2}-x_{O1})& \\ y_E=y_{O1}+\frac{r1}{r1+r2}\cdot (y_{O2}-y_{O1})& \\ {\varphi }_1=\arctan (\frac{y_C-y_{O1}}{x_C-x_{O1}})-\arctan (\frac{y_E-y_{O1}}{x_E-x_{O1}})& \\ {\varphi }_2=\arctan (\frac{y_D-y_{O2}}{x_D-x_{O2}})-\arctan (\frac{y_E-y_{O2}}{x_E-x_{O2}})& \end{array} \end{gathered}$$

位置速度计算模型基于调头弧线的优化

在问题一与问题二中，求解板凳龙各个圆孔的位置与速度，按照本文的思路需要计算螺线的积分，但是第四问中路线修改为螺线与双圆弧进行拼接，无法直接套用问题一中对于位置的求解模型，于是在第一问的基础上进行了优化，以$C$点处的弧度${\theta }_C$与$D$点弧度${\theta }_D$作为盘入圆弧与盘出圆弧的范围限制，在从$[{\theta }_D,{\theta }_C]$的弧度范围内，设定映射关系$g(\theta )$表述双圆弧段的板凳位置计算，则整个板凳龙所处位置与弧度的关系可表示为函数$f(x,y)$如下：
$$f(x,y)=\{\begin{array}{ll}r(\theta ),& \text{if}\theta \in [{\theta }_C，32\pi ]\\ g(\theta ),& \text{if}\theta \in [{\theta }_D,{\theta }_C]\\ -r(\theta ),& \text{if}\theta \in [-32\pi ，{\theta }_D]\end{array}$$
其中$r(\theta )$与$-r(\theta )$分别对应前段螺线和后段螺线的位置计算函数，可沿用第一问模型，对于$g(\theta )$，求解过程如下：

由于驶入圆弧前与盘出圆弧后的位置$C(x_C,y_C)$与$D（x_D,y_D)$可根据最短路径$L$求解导出，由此可导出对应的${\theta }_C$与${\theta }_D$，进而求解出圆弧部分总长度$S={\varphi }_1\cdot r_1+{\varphi }_2\cdot r_2$.

其中在整个圆弧路线中，盘入点弧度${\theta }_C$为其中弧度的最大值，即${\theta }_{max}$，同理盘出点${\theta }_D$为其中弧度的最小值，即${\theta }_{min}$，故可定义映射系数$k$，表示螺旋线开始到当前角度 $\theta$ 的路径长度相对于整个双圆弧路径长度的比例，其公式为：
$$k=\frac{\theta -{\theta }_{min}}{{\theta }_{max}-{\theta }_{min}}$$
基于比例系数与映射关系，可计算出此时在圆弧部分前进路程与后段圆弧长度的差值$\Delta x$及其含义为：
$$\Delta x=k\cdot S-{\varphi }_2\cdot r_2$$
其中：

• 如果 $\Delta x>0$，则当前点在前段圆弧上；
• 如果 $\Delta x\le 0$，则当前点在后段圆弧上。

若该点在前段圆弧上

首先需求解出该点在前段圆弧上的角度占比$\eta$的表达式为：
$$\eta =\frac{kS-r_2{\varphi }_2}{2\pi r_1}$$
最终导出该点的极坐标下的极角${\theta }_0$表达为：
$${\theta }_0=\arctan (\frac{y_E-y_{O1}}{x_E-x_{O1}})+2\pi \eta$$
最终将${\theta }_0$与圆心$O_1(x_{O1}，y_{O1})$的坐标结合即可推导出此时该点的位置坐标为：
$$f(x,y)=\{\begin{array}{l}x=x_{O1}+r_1\cos {\theta }_0\\ y=y_{O1}+r_1\sin {\theta }_0\end{array}$$
若该点在后段圆弧上

推导过程与前段圆弧类似，只是将其中相对的位置计算点改为$O_2(x_2,y_2)$，并且由于其在后段圆弧上，投影范围从$[{\theta }_D,{\theta }_C]$，则其圆弧长度小于后段圆弧的总长，故公式表达修改如下：
$$\eta =\frac{r_2{\varphi }_2-kS}{2\pi r_1}$$

$${\theta }_0=\arctan (\frac{y_E-y_{O2}}{x_E-x_{O2}})+2\pi \eta$$

$$f(x,y)=\{\begin{array}{l}x=x_{O2}+r_2\cos {\theta }_0\\ y=y_{O2}+r_2\sin {\theta }_0\end{array}$$

最终导出位置计算的全曲线表达式为：
$$\begin{array}{r}f(x,y)=\{\begin{array}{llc}r(\theta )& \mathrm{if}\theta \in [{\theta }_C,32\pi ]& \\ \{\begin{array}{l}x=x_{O1}+r_1\cos ({\theta }_0),\\ y=y_{O1}+r_1\sin ({\theta }_0),\end{array}& \mathrm{if}\theta \in [{\theta }_D,{\theta }_C]\mathrm{and}\Delta x>0& \\ \{\begin{array}{l}x=x_{O2}+r_2\cos ({\theta }_0),\\ y=y_{O2}+r_2\sin ({\theta }_0),\end{array}& \mathrm{if}\theta \in [{\theta }_D,{\theta }_C]\mathrm{and}\Delta x\le 0& \\ -r(\theta )& \mathrm{if}\theta \in [-32\pi ,{\theta }_D]& \end{array}\end{array}$$
由于数学公式过于繁琐，故在本问中速度的求解模型的建立只采取了问题一中以极小时间步长$\Delta t$进行数值模拟的方法，具体过程将置于模型的求解中进行。

模型的求解

本题中需要求解的部分是单目标优化模型最优值的求解与速度计算小步长数值模拟的求解。对于本题最短弦长单目标优化模型的求解，由于主要决策变量仅为${\theta }_C$与${\theta }_D$，计算量较小，故采取定步长搜索策略，具体步骤如下：

Step1 设定初始解向量和初始步长

设定初始解$[{\theta }_{C0},{\theta }_{D0}]$为$[\pi ,\pi ]$,步长均为$\frac{2\pi }{360}$，即角度制中的$1°$。

Step2 遍历解空间

采用双层循环的方式遍历解空间，通过固定其中一个变量，改变另一个变量，并计算每一对$[{\theta }_C,{\theta }_D]$的目标函数值，若比之前的最优值更好，则更新最优解和最优值。

Step3 得到最优解

将最终保存的最优解与最优值导出。

对于问题四速度的求解，仍沿用问题一中提出的求解办法，具体步骤参照问题一模型的求解部分。