# NOAI2025教学测试第4题：求解部分数据丢失的单摆运动

简介：此为NOAI2024第二轮试题的第4题，也是NOAI2025教学测试试题

## 一. 题目概述

提供一个单摆运动的数据集，存储在.csv文件中，包含两个变量：

t：时间，单位：秒(s)，单位为标准单位，在解题过程中不用考虑单位换算；

theta：摆角，带有正负，表示方向，单位为弧度，即使用sin函数运算时不用考虑单位换算。

请基于以上的数据集，使用PyTorch回归出来带有阻力的单摆运动的有关参数，补齐残缺数据并进行预测。

 

## 二. 数据集

（1）训练集：pendulum_train.csv ，训练集地址：[训练集](https://bohrium.dp.tech/competitions/1723157880?tab=datasets)；

（2）测试集A：pendulum_testA.csv  比赛过程中选手无法直接下载测试集A；

（3）测试集B：pendulum_testB.csv  比赛过程中选手无法直接下载测试集B；

其中训练集的数据全部可见，可以帮助选手生成求解参数的方法；测试集A的数据选手不可见，但是会在A榜中显示，可以帮助选手验证求解参数的方法是否正确；测试集B的数据选手不可见，最终在B榜计算分数时使用。

注意：训练集、测试集A和测试集B的微分方程参数均不相同，选手需要提交的是一种可以针对不同数据的求解方法。将不同数据集导入时，即训练集、测试集A和测试集B导入时均能求解出来正确参数。

 

## 三. 任务

​	在机器学习中我们有时会遇到需要借助少量数据完成对数据中的规律进行提取的场景。在这种场景下，如何充分利用先验知识（公式）处理数据、设计模型并成功求解未知参数是我们解决问题的关键所在。下面考虑一个单摆场景。

![image-20241020133718907](https://bohrium.oss-cn-zhangjiakou.aliyuncs.com/competition/35/005.jpeg)

​	如上图所示，现有一可以视为质点的**质量为1（$m=1$）**的小球通过一个长度为$l$的轻绳悬挂于一个固定点$O$，设$\theta (t)$ 为$t$时刻绳与过$O$铅垂线的夹角，称作摆角，将绳绷直并以初始的摆角${{\theta }_{0}}<\pi /2$（弧度制，不懂弧度制可以询问大语言模型）。在$t=0$时刻无初速释放小球，则小球可以在轻绳和铅垂线对应的平面进行运动。这里以小球的初始摆角${{\theta }_{0}}$为正，并在铅垂线左侧时，摆角$\theta (t)$的符号记为正，在铅垂线右侧时，摆角$\theta (t)$的符号记为负。设重力加速度为$g=9.8$（整个题目求解过程中可以不用考虑单位换算，直接对数值进行运算即可），考虑大小与速度成正比、方向与速度相反的空气阻力，其中空气阻力大小$\mu $在运动中保持恒定。

​	现有一传感器能够实现对小球摆角的精确记录，但是在实验过程中传感器出现了问题，它在记录时某时刻突然中断了若干秒后（$\ge 1s$）才被重启，且在小球没有停止运动时就已经停止工作。**已知在中断的时间段内的某一时刻${{t}_{Fput}}$开始小球受到了竖直向下大小恒定为$F$的外力，之后一直持续到运动结束**。请尝试通过已经记录的数据对小球的情况进行推断和预测。数据在文件.csv中，文件包括两列，第一列记录了小球运动的时间戳$t$，用变量t表示，第二列记录了小球的摆角信息$\theta $，用变量theta表示。我们将训练集数据中$\theta (t)$绘制如下：

![image-20241020133725296](https://bohrium.oss-cn-zhangjiakou.aliyuncs.com/competition/35/003.jpeg)

​	为了求解这个问题，你需要有微分方程的先验知识，在整个运动过程中。角速度$\omega (t)$指的是$\theta (t)$的瞬时变化率（存在正负），角加速度$a(t)$指的是角速度$\omega (t)$的瞬时变化率（存在正负），即$\omega (t)=\frac{d\theta (t)}{dt}$，$a(t)=\frac{d\omega (t)}{dt}$（不懂导数可以询问大语言模型）。则上述单摆问题中，$a(t)$与$\omega (t)$和$\theta (t)$之间应该满足如下微分方程：

​    $a(t)=-\alpha \cdot \omega (t)-\beta \sin \left( \theta (t) \right)$，其中$\alpha ,\beta $为参数，根据牛顿第二定律，$\alpha =\frac{\mu }{m}$，

​    当$0\le t<{{t}_{Fput}}$时，没有施加外力，此时，$\beta ={{\beta }_{1}}=\frac{g}{l}$；

​    当$t\ge {{t}_{Fput}}$时，施加了外力$F$，此时，$\beta ={{\beta }_{2}}=\frac{g}{l}+\frac{F}{ml}$.

​	利用微分方程，如果$\alpha ,{{\beta }_{1}},{{\beta }_{2}}$已知，当知道$t$时刻的角速度$\omega (t)$和摆角$\theta (t)$，就可以根据微分方程求出来$a(t)$。但是本题的问题是$\alpha ,{{\beta }_{1}},{{\beta }_{2}}$未知，只有记录在每个时刻的摆角数据$\theta (t)$，你需要根据提供的数据，**“回归”**出来$\alpha ,{{\beta }_{1}},{{\beta }_{2}}$，这是你在本题的核心任务，在得到$\alpha ,{{\beta }_{1}},{{\beta }_{2}}$后，你就相当于求得了整个的微分方程，就可以利用微分方程对小球的摆角数据进行补全和预测。

​	具体而言：你需要根据.csv中记录的时间t和theta数据，求解出以下参数：**

​	（1）绳子长度：$l$；

​	（2）空气阻力：$\mu $；

​	（3）中间施加的外力大小：$F$；

​	（4）预测在施加外力后，传感器探测结束后（数据记录结束后），下一次摆角$\theta (t)=0$的时刻：${{t}_{nextzerotheta}}$；

​	（5）中间施加外力的时间：${{t}_{Fput}}$。

## 四. 提交

​	请提交submission.ipynb文件，其中包含训练模型的全部过程和参数求解全部过程。

​    submission.ipynb需要能够生成三个文件：

​    （1）训练集的参数求解结果存储在submission_train.csv中；

​    （2）测试集A的参数求解结果存储在submissionA.csv中；

​    （3）测试集B的参数求解结果存储在submissionB.csv中。

​    参数的存储格式如下表，如果没有求解出部分参数，请填一个默认参数，例如1，不要空白：

| l   | miu | F   | t_nextzerotheta | t_Fput |
| --- | --- | --- | --------------- | ------ |
| 1   | 1   | 10  | 1               | 1      |

​    可以在baseline.ipynb中查询提交格式的参考。

​    baseline.ipynb地址：[NOAI2025教学测试第4题_baseline](https://bohrium.dp.tech/notebooks/17599263382)



## 五. 评分

1.最终的评分矩阵如下：其中$X\_pre$为选手对参数$X$预测值，$X\_real$为参数$X$的真实值；

​    （1）绳长$l$求解评分：${{S}_{1}}=\exp (-10\left| l\_pre-l\_real \right|)$；

​    （2）空气阻力$\mu $求解评分：${{S}_{2}}=\exp (-10\left| \mu \_pre-\mu \_real \right|)$；

​    （3）外力$F$求解评分：${{S}_{3}}=\exp (-\left| F\_pre-F\_real \right|)$；

​    （4）下一次摆角$\theta (t)=0$的时刻${{t}_{nextzerotheta}}$求解评分：${{S}_{4}}=\exp (-10\left| {{t}_{nextzerotheta}}\_pre-{{t}_{nextzerotheta}}\_real \right|)$；

​    （5）外力施加时刻${{t}_{Fput}}$求解评分：${{S}_{5}}=\exp (-10\left| {{t}_{Fput}}\_pre-{{t}_{Fput}}\_real \right|)$；

​    （6）最终得分：$Score=\frac{{{S}_{1}}+{{S}_{2}}+2{{S}_{3}}+2{{S}_{4}}+2{{S}_{5}}}{8}$。

2.未按照格式提交：0分。

 


**六. 附录（提示）：小爱同学使用PyTorch求解带有阻力的微分方程的案例**

​    现在有一个可以看作质点的小球进行带有阻力的直线运动。

​	![image-20241020133535721](https://bohrium.oss-cn-zhangjiakou.aliyuncs.com/competition/35/001.jpeg)

其中位移$s(t)$**，**速度为$v(t)=\frac{\text{d}s(t)}{\text{d}t}$**，**加速度为$a(t)=\frac{\text{d}v(t)}{\text{d}t}$**。**小球受到两个阻力，一个阻力和当前速度成正比，另一个阻力和当前的位移成正比。即小球的运动满足如下微分方程：

​	$a(t)=-\alpha v(t)-\beta s(t)$,

其中$\alpha $是与速度成正比的阻力参数，其中$\beta $是与位移成正比的阻力参数，

​    小爱同学将该运动过程的数据记录在了sv_data.csv中，包括两列数据：

​    t：时间，单位：s，即为标准单位，不用考虑单位换算

​    s：位移，单位，m，即为标准单位，不用考虑单位换算    

​	小爱同学使用PyTorch撰写了代码进行“回归”，求解了其中的$\alpha,\beta$确定了微分方程，并预测了再过5秒之后的小球的位移和速度。他的求解过程的代码如下。


```python
import numpy as np
import torch
from torch import nn
from scipy.integrate import odeint
import pandas as pd
import matplotlib.pyplot as plt
import math

###定义求导/差分运算###
def derivative(t, s):  #使用差分代替导数，可以用记录的位移值代替计算速度
    dsdt = (s[1:] - s[:-1]) / (t[1:] - t[:-1])  #这就意味着每次使用差分求一次导，数组的大小就要-1，s记录n个数据，则v只能记录n-1个数据
    return dsdt 

def load_data(csv_path):  #读取位移和时间数据 
    data = pd.read_csv(csv_path)   
    t = np.array(data["t"])  #读取时间
    s = np.array(data["s"])  #读取位移
    return t, s

###读取数据###
csv_path = "sv_data.csv"   
t, s = load_data(csv_path)  #读取数据
v = derivative(t, s) #计算各个时刻的速度
#print(v)

###定义求解（回归出）微分方程参数的网络###
class StraightLine(nn.Module):  #使用Pytorch搭建微分方程模型
    def __init__(self):
        super().__init__()
        self.alpha = nn.Parameter(torch.randn((1,), requires_grad=True))  #需要记录下来参数，requires_grad = True可以一直不用改，为默认参数
        self.beta = nn.Parameter(torch.randn((1,), requires_grad=True)) 
        
    def forward(self, theta, omega, add_force=False):
        a = - self.alpha * omega - self.beta * theta  #微分方程
        return a

###位移s，速度v，加速度a的数组大小统一，s需要删掉最后两个值，v需要删掉最后一个值，a不需要删###
#为了能够将数据导入到PyTorch中，需要将np.array和list转为torch.tensor，代码如下#
a = torch.from_numpy(derivative(t[:-1], v))  #每次使用差分求一次导，数组的大小就要-1
#print(a)
_v = torch.from_numpy(v[:-1]) #每次使用差分求一次导，数组的大小就要-1,s记录n个数据，则v只能记录n-1个数据,a只能记录n-2个数据
_s = torch.from_numpy(s[:-2]) #每次使用差分求一次导，数组的大小就要-1,s记录n个数据，则v只能记录n-1个数据,a只能记录n-2个数据

###开始进行训练###
num_epoch = 20000  #训练20000轮就差不多了
model = StraightLine()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)

###提取模型求解出来的参数并打印####
alpha = (model.alpha).detach().item() #打印模型参数，与v相关的阻力系数
beta = (model.beta).detach().item() #打印模型参数,与s相关的阻力系数
print("alpha:{}".format(alpha))
print("beta:{}".format(beta))

####使用odient库求解微分方程的方法
def v_func(sv, t, alpha, beta):  #定义一个微分方程，输入分别为sv（一个数对，分别为位移s和速度v；t为时间，直接这么写就行，alpha，beta分别是与v有关的阻力系数和与s有关的阻力系数）
    s,v = sv 
    a = - alpha * v - beta * s #微分方程
    return np.array([v,a])
sv_expect = odeint(v_func, [s[-1],v[-1]], [0,5], args=(alpha,beta))  #打印五秒后的位移和速度，其中s[-1]和v[-1]提取最后一个数，表示提取当前位置的数，0表示当前时刻，5表示5秒后
print("再过5秒后的位移s为",sv_expect[1,0]) #第一列为位移，第一行为当前时刻的位移
print("再过5秒后的速度v为",sv_expect[1,1]) #第二列为速度，第二行为当前时刻的速度
```