板凳龙–记录第一次参加数模国赛遇到的难题，哈哈
第一小问
这一部分就是完全使用matlab计算出每一个板凳的前后把手的位置，特点为点的位置由阿基米德螺旋线的方程用一个角度变量来代替繁琐的xy坐标运算，计算两点距离使用余弦定理计算，使计算公式只与变量theta有关，其中误差大小可以调节theta的取样精度以及距离的可接受范围delta_x来实现 以下是matlab的实现代码，如果有使用需求，可以在文末链接处下载
clc; clear; close all;
% 预分配储存所有时刻所有板凳前把手的角度的数组theta_all_chairs_head
theta_all_chairs_head = cell(224, 1);
% 运动时经历的时间，每秒计数一次
time = 0 : 1 : 300;
pitch = 55;
beta = pitch / (2*pi);
% 计算龙头位置
first_chair_head = [];
for each_time = 1 : length(time)
first_chair_head = [first_chair_head; get_ori_theta(time(each_time), pitch)]; %#ok<AGROW>
end
theta_all_chairs_head{1} = first_chair_head;
% 由于步长不同，先将长为286的龙头与第一节龙身的把手距离先计算
second_chair_head = [];
for each_time = 1 : length(time)
second_chair_head = [second_chair_head; get_next_head(theta_all_chairs_head{1}(each_time), 286, 32*pi, pitch)]; %#ok<AGROW>
end
theta_all_chairs_head{2} = second_chair_head;
% 对后面每一个龙身以及最后的龙尾后部，每次进行步长为165的计算
for each_bodies221 = 1 : 222
next_chair_head = [];
for each_time = 1 : length(time)
next_chair_head = [next_chair_head; get_next_head(theta_all_chairs_head{each_bodies221 + 1}(each_time), 165, 32*pi, pitch)]; %#ok<AGROW>
end
theta_all_chairs_head{each_bodies221 + 2} = next_chair_head;
end
% 将theta_all_chairs_head转到二维数组，方便直观得到答案
theta_all_chairs_head_array = zeros(224, length(time));
for i = 1 : 224
theta_all_chairs_head_array(i, :) = theta_all_chairs_head{i, 1}';
end
x_all = beta .* theta_all_chairs_head_array .* cos(theta_all_chairs_head_array);
y_all = beta .* theta_all_chairs_head_array .* sin(theta_all_chairs_head_array);
x300 = x_all(1 : 181, length(time));
y300 = y_all(1 : 181, length(time));
figure;
scatter(x300, y300, '.'); axis equal;
hold on; plot(x300, y300);
% 将位置结果保存到result1.xlsx中
%filename = 'result1.xlsx';
%for i = 1 : 224
% rangeStr1 = sprintf('B%d', 2*i);
% rangeStr2 = sprintf('B%d', 2*i+1);
% writematrix(x_all(i, :), filename, 'Range', rangeStr1);
% writematrix(y_all(i, :), filename, 'Range', rangeStr2);
%end
第二小问
通过一定的矩阵变换，详情看论文中叙述，通过把手的点坐标计算出板凳对应四个点的坐标，然后判断龙头的顶点有没有与最近一圈的板凳内层碰撞，碰撞原理较为简单，就是点到直线的距离小于一个阈值，记录碰撞时的时间，调用自定义函数get_theta_all_chairs_head(time, theta_end_max, pitch)求出对应时刻的各个把手对应的角度
clc; clear;
res = inf;
time = 300 : 500;
pitch = 55;
beta = pitch / (2*pi);
for i = 1 : length(time)
theta_head = get_ori_theta(time(i), pitch);
theta_end_max = theta_head + 2.5 * pi;
theta_all_chairs_head = get_theta_all_chairs_head(time(i), theta_end_max, pitch);
for j = 1 : 224 - 1
if (theta_all_chairs_head{j} < theta_head + 2*pi) && (theta_head + 2*pi < theta_all_chairs_head{j+1})
theta_detect = [theta_all_chairs_head(j), theta_all_chairs_head(j+1)];
E_detect_theta = theta_detect{1};
F_detect_theta = theta_detect{2};
E_head_theta = theta_head;
F_head_theta = theta_all_chairs_head{2};
C = beta .* [1.1667*E_head_theta*cos(E_head_theta) - 0.1667*F_head_theta*cos(F_head_theta) + ...
0.0915*F_head_theta*sin(F_head_theta) - 0.0915*E_head_theta*sin(E_head_theta), ...
1.1667*E_head_theta*sin(E_head_theta) - 0.1667*F_head_theta*sin(F_head_theta) + ...
0.0915*E_head_theta*cos(E_head_theta) - 0.0915*E_head_theta*cos(E_head_theta)];
B = beta .* [1.1667*E_detect_theta*cos(E_detect_theta) - 0.1667*F_detect_theta*cos(F_detect_theta) + ...
0.0915*E_detect_theta*sin(E_detect_theta) - 0.0915*F_detect_theta*sin(F_detect_theta), ...
1.1667*E_detect_theta*sin(E_detect_theta) - 0.1667*F_detect_theta*sin(F_detect_theta) + ...
0.0915*F_detect_theta*cos(F_detect_theta) - 0.0915*E_detect_theta*cos(E_detect_theta)];
A = beta .* [1.1667*F_detect_theta*cos(F_detect_theta) - 0.1667*E_detect_theta*cos(E_detect_theta) + ...
0.0915*E_detect_theta*sin(E_detect_theta) - 0.0915*F_detect_theta*sin(F_detect_theta), ...
1.1667*F_detect_theta*sin(F_detect_theta) - 0.1667*E_detect_theta*sin(E_detect_theta) + ...
0.0915*F_detect_theta*sin(F_detect_theta) - 0.0915*E_detect_theta*sin(E_detect_theta)];
% 计算C到直线AB的距离
% 计算向量AB和AC
AB = B - A;
AC = C - A;
% 计算叉乘结果的模长，即|AB x AC|
crossProductMagnitude = abs(AB(1)*AC(2) - AB(2)*AC(1));
% 计算向量AB的长度
AB_length = norm(AB);
% 点C到直线AB的距离
distance = crossProductMagnitude / AB_length;
disp(distance);
if distance < 0.1
res = time(i);
disp('Time Found!!!!!!!!!!!!');
disp('终止时刻为');
disp(res);
break;
end
end
end
end
if res > 1e9
disp('没有找到终止时刻！！！！');
end
res2_theta = get_theta_all_chairs_head(res, 32*pi, 55);
res2_theta_array = zeros(224, 1);
for i = 1 : 224
res2_theta_array(i) = res2_theta{i};
end
res2_x = beta .* res2_theta_array .* cos(res2_theta_array);
res2_y = beta .* res2_theta_array .* sin(res2_theta_array);
% 将位置结果保存到result2.xlsx中
filename = 'result2.xlsx';
for i = 1 : 224
rangeStr1 = sprintf('B%d', i+1);
writematrix(res2_x(i, :), filename, 'Range', rangeStr1);
rangeStr2 = sprintf('C%d', i+1);
writematrix(res2_y(i, :), filename, 'Range', rangeStr2);
end
第三小问
在调用第二问代码的基础上（哈哈，没有封装第二问代码哈，水平有限），配合题目要求，大致计算出螺距取值范围，再进一步改变范围精确数值
clc; clear;
pitch = 55 : 2 : 80;
delta_time = [];
for ii = 1 : length(pitch)
time_theory = get_time(900*pi / pitch(ii), pitch(ii));
res = inf;
time = 300 : 500;
beta = pitch(ii) / (2*pi);
for i = 1 : length(time)
theta_head = get_ori_theta(time(i), pitch(ii));
theta_end_max = theta_head + 2.5 * pi;
theta_all_chairs_head = get_theta_all_chairs_head(time(i), theta_end_max, pitch(ii));
for j = 1 : 224 - 1
if (theta_all_chairs_head{j} < theta_head + 2*pi) && (theta_head + 2*pi < theta_all_chairs_head{j+1})
theta_detect = [theta_all_chairs_head(j), theta_all_chairs_head(j+1)];
E_detect_theta = theta_detect{1};
F_detect_theta = theta_detect{2};
E_head_theta = theta_head;
F_head_theta = theta_all_chairs_head{2};
C = beta .* [1.1667*E_head_theta*cos(E_head_theta) - 0.1667*F_head_theta*cos(F_head_theta) + ...
0.0915*F_head_theta*sin(F_head_theta) - 0.0915*E_head_theta*sin(E_head_theta), ...
1.1667*E_head_theta*sin(E_head_theta) - 0.1667*F_head_theta*sin(F_head_theta) + ...
0.0915*E_head_theta*cos(E_head_theta) - 0.0915*E_head_theta*cos(E_head_theta)];
B = beta .* [1.1667*E_detect_theta*cos(E_detect_theta) - 0.1667*F_detect_theta*cos(F_detect_theta) + ...
0.0915*E_detect_theta*sin(E_detect_theta) - 0.0915*F_detect_theta*sin(F_detect_theta), ...
1.1667*E_detect_theta*sin(E_detect_theta) - 0.1667*F_detect_theta*sin(F_detect_theta) + ...
0.0915*F_detect_theta*cos(F_detect_theta) - 0.0915*E_detect_theta*cos(E_detect_theta)];
A = beta .* [1.1667*F_detect_theta*cos(F_detect_theta) - 0.1667*E_detect_theta*cos(E_detect_theta) + ...
0.0915*E_detect_theta*sin(E_detect_theta) - 0.0915*F_detect_theta*sin(F_detect_theta), ...
1.1667*F_detect_theta*sin(F_detect_theta) - 0.1667*E_detect_theta*sin(E_detect_theta) + ...
0.0915*F_detect_theta*sin(F_detect_theta) - 0.0915*E_detect_theta*sin(E_detect_theta)];
% 计算C到直线AB的距离
% 计算向量AB和AC
AB = B - A;
AC = C - A;
% 计算叉乘结果的模长，即|AB x AC|
crossProductMagnitude = abs(AB(1)*AC(2) - AB(2)*AC(1));
% 计算向量AB的长度
AB_length = norm(AB);
% 点C到直线AB的距离
distance = crossProductMagnitude / AB_length;
disp(distance);
if distance < 1
res = time(i);
delta_time = [delta_time; abs(res - time_theory)];
disp('Time Found!!!!!!!!!!!!');
disp('终止时刻为');
disp(res);
break;
end
end
end
end
if res > 1e9
disp('没有找到终止时刻！！！！');
delta_time = [delta_time; abs(res - time_theory)];
end
end
function res = get_time(theta_head, pitch)
% 用func1表示每个d(theta)对应的微分量
func1 = @(theta) ((pitch/(2*pi)*theta).^2 + (pitch/(2*pi))^2).^(0.5);
% 对func1从x_theta积分积到第16圈结束，计算路程
temp = integral(func1, theta_head, 16*2*pi);
% t = s / v
res = temp / 100;
end
第四与第五问没有完成，这一次比赛以失败告终，有了这一次的经验，下一次一定会更加优秀
以上代码中get开头的函数为自定义函数，如下所示
get_next_head(theta_start, interval, theta_end_max, pitch)
实现输入上一把手对应角度theta_start, 有效板长interval, 限制最大搜索范围theta_end_max, 阿基米德螺旋线的螺距pitch，寻找下一个把手点对应的角度
function res = get_next_head(theta_start, interval, theta_end_max, pitch)
% 防止res在theta_end超出范围后没有返回值，设定它的初始值为无穷大inf
res = inf;
% beta为螺距与2pi的比值，用于阿基米德螺线方程相关量的计算
beta = pitch / (2*pi);
% 从theta_start开始设置检测用的theta_detect_end的数值，向后检测1圈
% 这里的1圈，是在考虑没有运动到如此内层，以至于1圈不够用
theta_detect_end = theta_start + 6*pi;
theta_detect = theta_start : 0.000001 : theta_detect_end;
for i = 1 : length(theta_detect)
% 将每一次的检测角度赋值给theta_end
theta_end = theta_detect(i);
% 使用余弦定理来计算theta_end对应的“板凳长”
len1 = beta * theta_start;
len2 = beta * theta_end;
len_detect = (len1^2 + len2^2 - 2*len1*len2*cos(theta_end - theta_start))^(0.5);
% 确保theta_end的需求范围
if theta_end < theta_end_max
try
% 如果满足误差范围，则将theta_end赋值给结果res，替代inf
if abs(len_detect - interval) <= 0.03
res = theta_end;
return;
end
% 如果没有计算出数值，找出错误，跳出循环
catch
break;
end
else
% theta_end超出32pi，也就是16圈，回到主函数
return;
end
end
end
get_ori_theta(time, pitch)
实现输入时间time和螺距pitch，返回龙头前把手对应的角度res_theta
function res_theta = get_ori_theta(time, pitch)
% 计算时间time对应的龙头所走过的路程s = v * t
s = 100 * time;
syms theta x_theta
% 用func1表示每个d(theta)对应的微分量
func1 = ((pitch/(2*pi)*theta).^2 + (pitch/(2*pi))^2).^(0.5);
% 对func1从x_theta积分积到第16圈结束，计算路程
temp = int(func1, theta, x_theta, 16*2*pi);
temp_numeric = matlabFunction(temp);
% 两个路程数值相同计算出对应角度
func2 = @(x_theta) s - temp_numeric(x_theta);
res_theta = fsolve(func2, s/(10*pitch), optimset('Display', 'off'));
end
get_theta_all_chairs_head(time, theta_end_max, pitch)
这其实是一个封装的功能，用于计算所有把手节点对应的角度，返回一个cell数组theta_all_chairs_head
function res = get_theta_all_chairs_head(time, theta_end_max, pitch)
% 预分配储存所有时刻所有板凳前把手的角度的数组theta_all_chairs_head
theta_all_chairs_head = cell(224, 1);
% 计算龙头位置
first_chair_head = [];
for each_time = 1 : length(time)
first_chair_head = [first_chair_head; get_ori_theta(time(each_time), pitch)]; %#ok<AGROW>
end
theta_all_chairs_head{1} = first_chair_head;
% 由于步长不同，先将长为286的龙头与第一节龙身的把手距离先计算
second_chair_head = [];
for each_time = 1 : length(time)
second_chair_head = [second_chair_head; get_next_head(theta_all_chairs_head{1}(each_time), 286, theta_end_max, pitch)]; %#ok<AGROW>
end
theta_all_chairs_head{2} = second_chair_head;
% 对后面每一个龙身以及最后的龙尾后部，每次进行步长为165的计算
for each_bodies221 = 1 : 222
next_chair_head = [];
for each_time = 1 : length(time)
next_chair_head = [next_chair_head; get_next_head(theta_all_chairs_head{each_bodies221 + 1}(each_time), 165, theta_end_max, pitch)]; %#ok<AGROW>
end
theta_all_chairs_head{each_bodies221 + 2} = next_chair_head;
end
res = theta_all_chairs_head;
end
以下是资料下载地址
点击这里下载2024A题原题 A.pdf
点击这里下载代码及相关数据 附件_代码数据等.zip
点击这里下载论文 基于阿基米德螺旋线的“板凳龙”模型.pdf