关于一种Lato-lato小球模型的搭建和相关动力学方程的迭代拟合

Question

lato-lato是一种流行于东南亚的玩具，笔者在此对其进行模型抽象和尝试拟合出它的动力学方程

如图，将两个可视为质点的球分别连接到绳子的两端并且将绳子的中点固定到一个支点上，外力作用于支点，带动整个系统上下振动。当支点在竖直方向振动时，球将会不断相撞。为了更加简化模型，笔者假设在无重力环境下研究整个系统，并且角速度一直足够大，小球相对于支点在圆上运动且轻绳始终紧绷
在无重力条件下，一个支点对称连接着两根轻绳，轻绳与竖直向下方向的夹角为θ，支点到单个小球中心的距离为l。轻绳另一端分别连接着两个完全一样的小球，小球的半径为r。支点受到外力作用，不断上下振动，带动小球运动。由于小球的初速度足够大，小球相对于支点在一个圆周上运动。小球沿圆周运动，到最高点和最低点时两小球相撞，能量发生损失，速度变为原来的e倍（-1<e<0）。每次碰撞突变，系统实际上由于e损失了一些动能。F在不断为系统输入能量。F输入的能量恰好等于一个周期内两次碰撞损失的动能，于是一个周期结束后系统的能量没有发生变化 。
我们进行理论力学的分析并得到以下条件：
1.1. y’’= lθ’^2cosθ 其中y是系统质心的位置
1.2.F=2m lθ’^2cosθ 其中F是拉动支点的外力
1.3. ∫_t^(t+T)▒〖F(t)y'(t)dt〗=2e^2(mlθ’_1^2 mlθ’_2^2) 此处两个角速度分别指碰前和碰后的角速度，这个积分式代表着外力F输入的能量与一个周期内上下两次碰撞损失的动能恰好抵消
此处有三个自由度：F(t),θ(t),y(t)
笔者解这个方程组的思路是：首先给定一个系数已知的θ(t)，代入1.1，1.2得到y(t)和F(t)，然后将其代入1.3这个积分式的左边，检查其是否等于右边，如果不满足1.3的等式，则对θ(t)的系数进行微调，重复之前过程，如此迭代直至得到符合1.3积分式的参数。
之所以将θ(t)作为待定的函数，是因为θ’在θ=0或pi时会发生突变，根据这个直接相关的特征找θ(t)比较方便
根据实际情况和直觉，笔者认为θ(t)和θ‘(t)的图像大致如图

θ(t)先是从0增长到pi，这期间由于外力输入系统能量，角速度不断变大，θ(t)的导数也不断增大。然后在θ=pi时发生碰撞，θ’发生突变，变为原来的e倍。然后角速度又不断增大，对应θ(t)导数为负数，不断减小，直到θ到达0，完成一个周期。
由于突变，θ(t)是一个分段函数，在(0,t1)上θ(t)=θ1(t)，(t1,t2)上θ(t)=θ2(t)
由于我们进行了简化，不考虑重力，那么θ1(t)和θ2(t)应该有某种对称。具体地说，t1=1/2t2=1/2T，θ1(t)=pi-θ2(t-t1)。于是根据对称性，我们只用研究半个周期内，即θ1(t)的函数
θ(t)和θ‘(t)应该有一些特点，反应到迭代的过程中就应该有一些约束：
2.1.θ1(t1)=pi, θ1(0)=0 这个式子意味着θ1(t)从0到pi
2.2.-eθ'1(t1)=θ'1(0) 这个式子意味着突变时刻碰撞后的速度变为原来的e倍
2.3.θ'1(t)>0,θ’’1(t)>0, 这个式子象征着实际情况：角速度在两次碰撞之间方向不变且由于外部能量输入一直在变大
于是笔者提出以下六次多项式拟合θ1(t)：
θ_1 (t)=∑_1^6▒〖c_nπ(t/t_1 )^n 〗
一阶和二阶导数分别为：
〖θ'〗_1 (t)=∑_1^6▒〖nc_nπ*(t/t_1 )^(n-1) 〗
〖θ''〗_1 (t)=∑_2^6▒〖n(n-1)c_nπ*(t/t_1 )^(n-2) 〗
这个多项式将2.1与2.2转化为：
3.1.∑_1^6▒c_n =1
3.2.c1=-e(∑_1^6▒〖n*c_n 〗)
以方便计算。现在笔者根据预实验得到的数据，给出一些合理的初始数据认为
3.3.θ‘1(0)=5pi
以及e=-0.7，t1=0.2，t2=T=0.4。将结合3.1、3.2、3.3，我们可以得到一组特殊解：[1, -0.71, 0.71, 0, 0, 0]
经过程序画出函数图，这个特殊解实际上满足了约束2.1，2.2，但略微违背了约束2.3.已经与理想的需要的函数较为接近。

接下来我们回到1.1，1.2，1.3三个式子。我们将得到的系数特殊解代入，得到一个θ1(t)。然后利用该θ1(t)得到对应的y、F函数，并算出1.3积分式左边部分，检验是否等于右边部分。不等于则对c_n进行迭代。因为迭代程序会破坏函数的合理性，因此我们应在迭代程序中加入以下约束条件：
4.1.∑_1^6▒c_n=1
4.2.c1=-e(∑_1^6▒〖n*c_n 〗)
4.3. θ'1(t)>0,θ’’1(t)>0
4.4. θ1，θ‘1，θ’‘1之间符合导数关系

接下来是笔者的困惑：我的编程能力较弱，在gpt的帮助下完成了前面的部分。但当涉及到迭代的部分，迭代出的效果一直不甚理想，即总是会破坏迭代的约束条件。笔者目前已经尝试了遗传算法、约束优化方法，迭代求积分法，效果均不佳。希望有人能指教笔者，在当前问题的背景下，什么方法能更好地在不破坏4.1,4.2,4.3,4.4四个约束条件的前提下迭代出满足1.3积分式的θ1(t)的参数

Answer 1

以下回复参考：皆我百晓生、券券喵儿等免费微信小程序作答:

对于这个问题，确实需要一种有效的迭代方法来求解满足所有约束条件的θ1(t)。考虑到问题的复杂性，可以尝试以下几种方法来改进迭代过程：

1. 采用带有约束的优化算法： 针对带有约束的优化问题，可以使用如序列二次规划（Sequential Quadratic Programming）或内点法（Interior Point Method）等优化算法。这些算法可以在满足约束条件的情况下寻找最优解。

2. 拉格朗日乘数法： 对于包含等式和不等式约束的优化问题，拉格朗日乘数法是一种常用的数学规划技术。可以尝试用这种方法将约束条件融入目标函数中，然后寻找满足条件的解。

3. 粒子群优化算法（Particle Swarm Optimization）： 这是一种基于群体智能的优化算法，通过模拟鸟群、鱼群等生物的社会行为来寻找最优解。这种方法对初始值不敏感，可能能够更好地处理您遇到的问题中的非凸性。

4. 差分进化算法（Differential Evolution）： 该算法是一种简单、高效、鲁棒性强的全局优化算法，适用于多维参数优化问题。它可以处理复杂的非线性问题，并且在迭代过程中能够保持解的多样性。

5. 模拟退火技术： 模拟退火是一种用于解决组合优化问题的随机搜索算法。它可以避免陷入局部最优解，从而有可能找到全局最优解。通过调整温度参数和冷却速度，可以在一定程度上处理约束条件。

在进行迭代的过程中，确保遵循以下步骤来提高效率：

• 选择合适的初始值：选择合理的初始参数值可以提高算法的收敛速度。可以尝试根据问题的特点，选择一个接近真实解的初始值。
• 动态调整约束权重：在迭代过程中，根据解的进展动态调整约束条件的权重，以帮助算法更快地找到满足条件的解。
• 利用问题特性设计算法：针对问题中的物理意义（如角速度的变化、能量的输入与损失等），设计更为贴近实际物理过程的迭代策略。
• 使用适当的评估函数：设计一个评估函数来评估解的优劣，并根据这个函数值来更新解。确保评估函数能够反映所有约束条件以及目标函数的要求。

最后，编程实现时需要注意代码的质量和效率，确保迭代过程稳定并且能够快速收敛到满足所有约束条件的解。如果可能的话，寻求一个具有优化算法经验的程序员或者专家帮助来实现这个过程。