如何填充MATLAB反作用轮倒立摆动力学函数RW_pendulum_dynamics输出状态向量导数

前置假设 #

• 摆杆为均质细杆，无弹性形变，转动惯量取绕根部铰接点的计算值
• 反作用轮为均质薄圆盘，安装在摆杆的自由端，转轴与摆杆摆动平面垂直
• 系统所有摩擦、间隙损耗忽略不计
• 状态向量中反作用轮角位置alpha为相对于摆杆的相对角位置，角速度为相对角速度

转动惯量表达式 #

• I1为均质细杆绕一端铰接点的转动惯量：I1 = (1/3)*m1*l1^2;
• I2为均质薄圆盘绕中心转轴的转动惯量：I2 = (1/2)*m2*wr^2;

状态导数分量计算 #

状态向量定义：y = [alpha; alpha_dot; theta; theta_dot]，因此导数分量满足：

1. 第一个分量为反作用轮角位置的导数，即相对角速度：dy(1,1) = y(2);
2. 第三个分量为摆杆角位置的导数，即摆杆角速度：dy(3,1) = y(4);
3. 第二个分量为反作用轮相对角加速度、第四个分量为摆杆角加速度，通过拉格朗日动力学方程联立求解得到。

完整填充后的函数 #

function dy = RW_pendulum_dynamics(y, m1, m2, l1, wr, g, u)
  % 转动惯量计算
  I1 = (1/3) * m1 * l1^2;
  I2 = (1/2) * m2 * wr^2;
  % 运动学导数
  dy(1,1) = y(2);
  dy(3,1) = y(4);
  % 动力学求解角加速度
  tau_g = g * l1 * sin(y(3)) * (m1/2 + m2);
  I_theta = I1 + m2 * l1^2;
  dy(4,1) = (tau_g - 2*u) / I_theta;
  dy(2,1) = u/I2 - dy(4,1);
end
% 若使用Octave环境，将末尾的end替换为endfunction即可

内容的提问来源于stack exchange，提问作者Harsha