function dotq = open_link_PID_ode_param ...
(t,q, l1,lc1,m1,J1, l2,lc2,m2,J2, l3,lc3,m3,J3, l4,lc4,m4,J4, g, ...
 x1,y1, x3,y3, alpha, ...
 theta1d, theta3d, Kp1, Ki1, Kd1, Kp3, Ki3, Kd3)
% •ÂƒŠƒ“ƒN‹@�\ PID�§Œä •W�€Œ`
    theta1 = q(1); theta2 = q(2); theta3 = q(3); theta4 = q(4);
    omega1 = q(5); omega2 = q(6); omega3 = q(7); omega4 = q(8);
    xi1 = q(9); xi3 = q(10);
    
    Hleft  = open_link_inertia_matrix ( theta1, theta2, l1,lc1,m1,J1, l2,lc2,m2,J2 );
    Hright = open_link_inertia_matrix ( theta3, theta4, l3,lc3,m3,J3, l4,lc4,m4,J4 );
    
    [ lleft,  uleft  ] = open_link_lower_upper_terms ...
        (t, theta1,theta2,omega1,omega2, l1,lc1,m1,J1, l2,lc2,m2,J2, g);
    [ lright, uright ] = open_link_lower_upper_terms ...
        (t, theta3,theta4,omega3,omega4, l3,lc3,m3,J3, l4,lc4,m4,J4, g);
    
    [ gx, gy, Cx, Cy ] = closed_link_CSM ...
        ( theta1, theta2, theta3, theta4, q(5:8), ...
          l1, l2, l3, l4, x1, y1, x3, y3, alpha );
    
    tauleft  = [ lleft  - Kp1*(theta1-theta1d) - Kd1*omega1 - Ki1*xi1; ...
                 uleft ];
    tauright = [ lright - Kp3*(theta3-theta3d) - Kd3*omega3 - Ki3*xi3; ...
                 uright ];
    
    Hlr = [ Hleft, zeros(2,2); zeros(2,2), Hright ];
    A = [gx, gy];
    H = [ Hlr, -A;
          -A', zeros(2,2) ];
    tau = [ tauleft; tauright; Cx; Cy ];
    dotomega = H\tau;
    dotq = [ omega1; omega2; omega3; omega4; dotomega(1:4); theta1-theta1d; theta3-theta3d ];
end
