基于修正MD-H模型对机器人进行运行学建模，存在几何参数有a，α，d，θ和β。当这些参数存在微小误差时，机器人的实际相邻连杆之间的变换关系和理论相邻连杆之间变换关系会存在一定的偏差，导致最后实际和理论的末端位姿坐标也存在误差，分别用 Δa、Δα、 Δd,、 Δθ;和 Δβ;来表示MD-H模型中的五个几何参数误差。利用微分变换原理将机器人各个连杆机构之间的微小原始偏差合成积累到末端位姿的误差视为各个连杆机构进行微分变换综合作用导致的结果，基于MD-H运动学模型建立误差模型，由于各个连杆机构都存在几何参数的误差，机器人的相邻连杆之间的变换矩阵也存在着微小偏差，根据微分运动变换原理，连杆之间的实际变换矩阵和理论变换矩阵存在一定关系。 帮我用MATLAB实现结合我做建立的机器人模型和DH参数，建立误差模型。并且举例我输入关节角的值能够得到误差值。clear all; clc; du = pi/180; a = [0+0.001, 185+0.0079, 0+0.005, 120+0.12]; alpha = [pi/2+0.003, 0+0.001, pi/2+0.005, pi/2]; d = [0+0.001, 0+0.0079, 90+0.005, 0+0.12]; theta = [90*du+0.02, 0, 0.023, 0.08]; beta = zeros(1, 4)+0; L1(1) = Link('d', d(1), 'a', a(1), 'alpha', alpha(1), 'qlim', [180*du, 365*du], 'modified'); L1(2) = Link('d', d(2), 'a', a(2), 'alpha', alpha(2), 'qlim', [3*du, 63*du], 'modified'); L1(3) = Link('d', d(3), 'a', a(3), 'alpha', alpha(3), 'qlim', [60*du, 120*du], 'modified'); L1(4) = Link('d', d(4), 'a', a(4), 'alpha', alpha(4), 'qlim', [230*du, 326*du], 'modified'); Needle = SerialLink(L1, 'name', 'Needle'); T1 = DH(1, a(1), alpha(1), d(1), theta(1)+beta(1)); T2 = DH(2, a(2), alpha(2), d(2), theta(2)+beta(2)); T3 = DH(3, a(3), alpha(3), d(3), theta(3)+beta(3)); T4 = DH(4, a(4), alpha(4), d(4), theta(4)+beta(4)); T = T1 * T2 * T3 * T4; delta_a = 0.001; delta_T = zeros(4, 4);帮我续写代码保证能够正确运行

% 定义误差模型
syms a1 a2 a3 a4 alpha1 alpha2 alpha3 alpha4 d1 d2 d3 d4 theta1 theta2 theta3 theta4 beta1 beta2 beta3 beta4 delta_a delta_alpha delta_d delta_theta delta_beta
delta_T = [1, -delta_a, delta_d*sin(alpha1), theta1*delta_d*cos(alpha1);
delta_a, 1, -delta_d*cos(alpha1), -theta1*delta_d*sin(alpha1);
0, delta_alpha, 1, delta_a;
0, 0, 0, 1];
delta_T = delta_T * [1, 0, 0, 0;
0, cos(alpha2), -sin(alpha2), -d2*sin(alpha2);
0, sin(alpha2), cos(alpha2), d2*cos(alpha2);
0, 0, 0, 1];
delta_T = delta_T * [cos(theta2+beta2), -sin(theta2+beta2), 0, a2;
sin(theta2+beta2), cos(theta2+beta2), 0, 0;
0, 0, 1, d3;
0, 0, 0, 1];
delta_T = delta_T * [cos(theta3+beta3), -sin(theta3+beta3), 0, a3;
sin(theta3+beta3), cos(theta3+beta3), 0, 0;
0, 0, 1, d3;
0, 0, 0, 1];
delta_T = delta_T * [cos(theta4+beta4), -sin(theta4+beta4), 0, a4;
sin(theta4+beta4), cos(theta4+beta4), 0, 0;
0, 0, 1, d4;
0, 0, 0, 1];

% 给定关节角，计算误差
q = [theta1, theta2, theta3, theta4];
q_delta = [delta_theta, delta_beta, delta_alpha, delta_a, delta_d];
J = jacobian(T, q);
J_inv = inv(J);
J_delta = jacobian(delta_T, q_delta);
T_delta = J_delta * q_delta';
T_actual = T * T_delta;
position_error = simplify(T_actual(1:3, 4) - T(1:3, 4));
orientation_error = simplify(tr2rpy(T_actual(1:3, 1:3)) - tr2rpy(T(1:3, 1:3)));

% 给定关节角的值，计算误差
q_val = [90*du, 0, 0.023, 0.08];
q_delta_val = [0, 0, 0, 0, 0.001];
position_error_val = double(subs(position_error, [theta1, theta2, theta3, theta4, delta_a, delta_alpha, delta_d, delta_theta, delta_beta], [q_val, q_delta_val]));
orientation_error_val = double(subs(orientation_error, [theta1, theta2, theta3, theta4, delta_a, delta_alpha, delta_d, delta_theta, delta_beta], [q_val, q_delta_val]));
fprintf("Position error: [%.4f, %.4f, %.4f]\n", position_error_val(1), position_error_val(2), position_error_val(3));
fprintf("Orientation error: [%.4f, %.4f, %.4f]\n", orientation_error_val(1), orientation_error_val(2), orientation_error_val(3));