A small question about a property called "actuator_actearly" in mujoco model

[Yingfan99327]

Intro

Hi!

I am a student majoring in Robotics, I use MuJoCo for my research on human robot interactions.

My setup

I use Mujoco with Python API.

My question

And I am really interested about the property called "actuator_actearly" in MjModel, I read the source codes in the engine_forward.c and find there is a if-else module here:

if (m->actuator_actearly[i]) {
    act = nextActivation(m, d, i, act_adr, d->act_dot[act_adr]);
} else {
    act = d->act[act_adr];
}

The definition for this property is "step activation before force". I did a small experiment that I set the initial ctrl, act both to zero and incremented the excitation signals a small step for each loop and print the value for this parameter and it shows "0", in other words, it will not go to the first branch in this if-else.
So my question for this is as follows:

1) For the muscle model in Mujoco, this value is set to 0 defaultly? which means that the activation used to calculate the active muscle force is the current activation instead of the next activation after integration, right?

2) If not, I guess this parameter would be compiled when I load the model, so is there a way to modify this parameter in the xml file, I didn't find too much information about this.

Thanks very much for your help and wishing you all the best.

Minimal model and/or code that explain my question

If you encountered the issue in a complex model, please simplify it as much as possible (while still reproducing the issue).

Model:

minimal XML

<mujoco>
  <worldbody>
    <light pos="0 0 1"/>
    <geom type="sphere" size="1" rgba="1 0 0 1"/>
  </worldbody>
</mujoco>

Code:

import mujoco
import mediapy as media
model = mujoco.MjModel.from_xml_string(xml)
data = mujoco.MjData(model)

with mujoco.Renderer(model) as renderer:
  mujoco.mj_forward(model, data)
  renderer.update_scene(data)

  media.show_image(renderer.render())

Confirmations

• I searched the latest documentation thoroughly before posting.
• I searched previous Issues and Discussions, I am certain this has not been raised before.

[Yingfan99327]

Sorry that the model and codes in my discussion are not what I used in the experiment. I didn't pay attention to it. Please just check the text info, thx.

[yuvaltassa]

Read this and this

[Yingfan99327]

Thanks for your help, really appreciate.

[Yingfan99327]

May I ask one more question, yuval. I am trying to implement my own muscle dynamics model to calculate the muscle force and hope to directly use the muscle force or joint torque resulted from the muscle force as the control signal to control the muscle. I am not sure whether it is possbile to realize in mujoco, because I think in mujoco, if I define a muscle as an actuator, then the control signal, data.ctrl should be the excitation signals. Is there a way in mujoco to use the muscle force (namely, actuation_force) or joint torque(qfrc_actuator) as the input to control a singal muscle?

I tried to create a single muscle model xml file as follows:

<actuator>
        <general name="DELT1" tendon="DELT1_tendon" biasprm="0.33707 1.53316 1051.85 1 0 1.8 10 0.518315 1.4 0" class="upperarm_muscle" gainprm="0.33707 1.53316 1051.85 1 0 1.8 10 0.518315 1.4 0" lengthrange="0.123363 0.240922"/>
        <motor name="shoulder_joint" joint="shoulder_rot"  gear="1" />
</actuator>

that I want to use the joint torque to control the joint.

[yuvaltassa]

Yes. You are doing the right thing. A motor sends force/torque directly to a tendon/joint.

[Yingfan99327]

Hi Yuval,

Thanks for your reply and sorry for my insist on this problem.

I have anothe two questions about the experiments I hope to do.

Question1:

I met a problem when I implemented the muscle model on my own. The first step I did is to use the muscle model you implemented in the "engine_forward.c" and I just wrote them in Python because I tried to directly use the torque to drive the joint instead of the excitation signal. And I tested in a single muscle model.

Here is the main structure of the codes I use to do the test:

last_act = 0
next_act = 0
gain_prm = model.muscle_gainprm
bias_prm = model.muscle_biasprm
dyn_prm = model.muscle_dynprm

# return to the zero configuration first
data.qpos = 0.0
data.qvel = 0.0
data.qacc = 0.0
data.act = 0.0
data.ctrl = np.zeros(2)
mujoco.mj_inverse(model, data)
mujoco.mj_comPos(model, data)

for i in range(epochs)

    ctrl = i/10
   
    # the function to calculate the activation dynamics which is exactly the same as implemented in Mujoco
    act_dot = cal_activation_dot(ctrl, last_act, dyn_prm)
    next_act = calc_next_activation(last_act, act_dot)

    mujoco.mj_tendon(model, data)
    actuator_length = data.actuator_length
    actuator_velocity = data.actuator_velocity
    lengthrange = model.actuator_lengthrange
    
    # the function used to calculate the muscle actuation force 
    F = single_muscle.calc_muscle_force(last_act,
                                        next_act,
                                        actuator_length,
                                        actuator_velocity,
                                        lengthrange,
                                        gain_prm,
                                        bias_prm)

    muslce_moment_arm = data.actuator_moment[0]
    tau_actuator = F * muslce_moment_arm

    tau_bias = data.qfrc_bias
    tau_passive = data.qfrc_passive
    tau = tau_actuator - tau_passive

    data.ctrl[1] = tau
    mujoco.mj_step(model, data)
    last_act = next_act 

And the model file is as follows:

<mujoco model="Single_muscle">

    <include file="basic_scene.xml"/>

    <compiler inertiafromgeom="auto" balanceinertia="true" boundmass="0.001" boundinertia=".0001" angle="radian" meshdir="" autolimits="true"/>

    <size njmax='1000' nconmax='1000'/> 

    <visual>
      <global offwidth="1000" offheight="1000"/>
        <headlight diffuse="1 1 1"/>
    <scale actuatorwidth=".075"/>
    </visual>
    
    <option integrator="implicitfast" impratio="10"/>

    <default class="main">

        <geom contype="0" conaffinity="0"/>

        <tendon rgba="0.95 0.3 0.3 41" width="0.001"/>

        <default class="arm">

            <joint limited="true" damping="1.05" armature="1e-4" stiffness="0"/>
            <geom margin="0.001" rgba="0.8 0.8 0.8 1" contype="0" conaffinity="0"/>
            <site size="0.0025" group="2"/>
            <tendon rgba="0.95 0.3 0.3 41" width="0.003"/>

            <default class="muscle">
                <general biasprm="0.75 1.05 -1 200 0.5 1.6 1.5 1.3 1.2 0" biastype="muscle" ctrllimited="true" ctrlrange="0 1" dynprm="0.01 0.04 0 0 0 0 0 0 0 0" dyntype="muscle" gainprm="0.75 1.05 -1 200 0.5 1.6 1.5 1.3 1.2 0" gaintype="muscle"/>
            </default>

        </default>

        <default class="upperarm_muscle">
            <general biasprm="0.75 1.05 -1 200 0.5 1.6 1.5 1.3 1.2 0" biastype="muscle"
                     ctrllimited="true" ctrlrange="0 1" dynprm="0.01 0.04 0 0 0 0 0 0 0 0"
                     dyntype="muscle" gainprm="0.75 1.05 -1 200 0.5 1.6 1.5 1.3 1.2 0" gaintype="muscle"/>
            <tendon rgba="0.95 0.3 0.3 41" width="0.007"/>
        </default>

    </default>

    <asset>

        <mesh file="arm/clavicle.stl" name="clavicle_geom_1_clavicle"/>
        <mesh file="arm/scapula.stl" name="scapula_geom_1_scapula"/>
        <mesh file="arm/humerus.stl" name="humerus_geom_1_humerus"/>

        <mesh name="clavicle" file="arm/clavicle.stl" scale="1 1 1"/>
        <mesh name="scapula" file="arm/scapula.stl" scale="1 1 1"/>
        <mesh name="humerus" file="arm/humerus.stl" scale="1 1 1"/>

        <texture name="texplane" type="2d" builtin="checker" rgb1=".2 .19 .2" rgb2=".1 0.11 0.11" width="50" height="50"/>
        <texture name="texgeom" type="cube" builtin="flat" mark="cross" width="127" height="1278" rgb1="0.7 0.7 0.7" rgb2="0.9 0.9 0.9" markrgb="1 1 1" random="0.01"/>
        <material name="MatPlane" reflectance="0.5" texture="texplane" texrepeat="4 4" texuniform="true"/>
        <material name="geom" texture="texgeom" texuniform="true"/>
        <material name="MatSkin" specular="0.2" shininess=".1" reflectance="0.5" rgba="0.9 .7 .5 1"/>
    </asset>

    <tendon>
        <spatial name="DELT1_tendon" class="upperarm_muscle">
            <site site="DELT1_DELT1-P1"/>
            <site site="DELT1_DELT1-P2"/>
            <site site="DELT1_DELT1-P3"/>
            <site site="DELT1_DELT1-P4"/>
        </spatial>

    </tendon>

    <actuator>
        <general name="DELT1" tendon="DELT1_tendon" biasprm="0.33707 1.53316 0 1 0 1.8 10 0.518315 1.4 0" class="upperarm_muscle" gainprm="0.33707 1.53316 0 1 0 1.8 10 0.518315 1.4 0" lengthrange="0.123363 0.240922"/> I set the force term in gainprm and biasprm to 0 to delete the torque resulted in the muscle actuator, if not, I found the qacc would be twice.
        <general name="shoulder_joint" joint="shoulder_rot" biasprm="0 0 0 0 0 0 0 0 0 0" gainprm="1 0 0 0 0 0 0 0 0 0"/>
    </actuator>
    
    <worldbody>

        <body name="clavicle"  pos="0.006325 -0.02547 0.00693" childclass="arm" euler="1.57 -1.57 0">

            <inertial diaginertia="0.001 0.001 0.001" mass="0.156" pos="-0.011096 0.0063723 0.054168" quat="0.906292 0.0334036 0.164622 -0.387838"/>

            <geom mesh="clavicle_geom_1_clavicle" name="clavicle_geom_1" type="mesh" />
            <site name="DELT1_DELT1-P4" pos="-0.014 0.01106 0.08021"/>

            <body name="clavphant"  pos="-0.01433 0.02007 0.1355">

                <inertial diaginertia="1 1 1" mass="0.001" pos="0 0 0"/>

                <body name="scapula"  pos="0 0 0">

                    <inertial diaginertia="0.00199576 0.00103993 0.001" mass="0.70396" pos="-0.054694 -0.035032 -0.043734" quat="0.453021 0.77921 0.0744129 0.426692"/>
                    <geom mesh="scapula_geom_1_scapula" name="scapula_geom_1" type="mesh"/>

                    <site name="DELT1_DELT1-P3" pos="0.04347 -0.03202 0.00499"/>

                    <body name="scapphant"  pos="-0.00955 -0.034 0.009">

                        <inertial diaginertia="1 1 1" mass="0.001" pos="0 0 0"/>
                        <body name="humphant" pos="0 0 0">

                            <inertial diaginertia="1 1 1" mass="0.001" pos="0 0 0"/>

                            <body name="humphant1"  pos="0 0 0">

                                <inertial diaginertia="1 1 1" mass="0.001" pos="0 0 0"/>

                                <body name="humerus"  pos="0 0 0">

                                    <inertial diaginertia="0.0128576 0.0121576 0.00239265" mass="1.99757" pos="0.018064 -0.140141 -0.012746" quat="0.407851 0.343714 -0.567755 0.627034"/>
                                    
                                    <joint axis="0 0 1" name="shoulder_rot" pos="0 0 0" range="-1.571 2.094" user="0.161391"/>

                                    <geom mesh="humerus_geom_1_humerus" name="humerus_geom_1" type="mesh" />
                                    <site name="DELT1hh_sidesite" pos="0.01226 -0.01367 -0.009599"/>

                                    <site name="DELT1_DELT1-P1" pos="0.00896 -0.11883 0.00585"/>
                                    <site name="DELT1_DELT1-P2" pos="0.01623 -0.1033 -0.00412"/>

                                </body>
                            </body>
                        </body>
                    </body>
                </body>
            </body>
        </body>
</worldbody>

</mujoco>

The ground-truth I use to compare is the qpos, qvel and qacc obtained after mj_step by using the same excitation signals as the data.ctrl singal. There is only one difference in the ground-truth model file which is the actuator,

<actuator>
        <general name="DELT1" tendon="DELT1_tendon" biasprm="0.33707 1.53316 1051.85 1 0 1.8 10 0.518315 1.4 0" class="upperarm_muscle" gainprm="0.33707 1.53316 1051.85 1 0 1.8 10 0.518315 1.4 0" lengthrange="0.123363 0.240922"/>
</actuator>

So, the experiment I did is to compare the qpos, qvel, qacc of the joint between using excitation signals in data.ctrl and joint torque in data.ctrl because I want to directly use the joint torque to control the bone.

Suppose the muscle actuation function is implemented correctly by me in Python, which has also been checked by calling many functions before and after mj_step, for instance, I checked gain, bias, actuation_force etc. And I just read the codes in engine_forward.c, but the problem is that

When data.ctrl=0.0, the results are good.
When data.ctrl=0.1, the qfrc_actuator before the mj_step are correct but after the mj_step are different. 
qfrc_actuator before mj_step:
a) for my own implementation: [3.64496526]
b) for mujoco: [3.64496534]
qfrc_actuator after mj_step:
a) for my own implementation: [3.64496526]
b) for mujoco:  [3.64101303]
qacc_smooth after mj_step:
a) for my own implementation: [60.07820937]
b) for mujoco: [60.1534012]

I thought the problem would be the tau I use to directly drive the joint is wrong with a very small difference. The tau I set when ctrl=0.1 is still [3.6450] but I tested that if I set the tau=[3.6408542602], then I would get the same results as shown by using Mujoco.

I don't know where this minor difference comes from, and I am not clear what I missed when I calculate the right joint torques to drive the joint. I checked the codes of calculating the muscle force and I use a very simple model without any constraints and also compared with the source codes, but I just don't know where this minor difference comes from and the difference in qacc would be enlarged a lot if I iterate for more epochs, for instance, 1000000. Would numerical precision matters in this experiment. Could you give me some advice or some notes I should pay attention to if I calculate the joint torques?

Question2:

Because of the problem above, due to the potential numerical precision problem in calculation I don't know (because I don't know whether the numerical precision does count in mj_solveM function etc), I was thinking about whethere there are other options to avoid this and also help me achieve my goal. Instead of calculating the joint torques, is it possible by using muscle actuation force, which means I only need to implemented the mju_muscleGain and mju_muscleBias function on my own and then set them by calling some function. Really appreciate that if you could give me some advice on my experiments.

Thanks very much and wishing you all the best.