使用Moco工具生成人体蹲起运动并分析辅助器械作用

本文链接：https://blog.youkuaiyun.com/weixin_50084948/article/details/141717936

1. 思路和解释

2. 代码逻辑

3. 代码注释

4. 关节扭矩驱动与肌肉驱动的模型处理方法

1. 思路和解释

Moco可以通过轨迹优化生成预测的运动，并用来进行优化模型参数。本篇文章以蹲起运动为例子进行说明。优化目标为最小化肌肉激活度平方和(权重：1.0)与动作持续时间（权重：1.0）。初始状态被设定为蹲姿[74]，最终状态被设定为直立站姿。没有使用动捕的运动数据。该模型包含了躯干和9块肌肉驱动的单右腿，其中6块肌肉设置有顺应性肌腱动力学（腘绳肌、股直肌、股外侧肌、腓肠肌、比目鱼肌和胫骨前肌）。为了简化将模拟双腿合为一腿，力量加倍。足部与地面设定为焊接接触，进行了简化。模型中包含了先前mocoInverse相同的髌骨的运动学约束，初始猜测值是变量界限的中点，肌肉激励和激活的初始猜测值是0.05。起立过程中，臀大肌、腘绳肌和股外侧肌等伸肌表现出了显著的活动性。

关节扭矩驱动与肌肉驱动的模型处理方法中，关节扭矩驱动会移除模型中的肌肉并添加致动器，肌肉驱动中，会对肌肉参数进行配置。

2. 代码逻辑

study(study.solve)_关节扭矩驱动的轨迹优化，或者说是运动预测

1. 导入.osim肌肉骨骼模型，创建MocoStudy和problem

2. 设置状态量和时间的边界条件，addGoal设置最小力为cost

3. 配置solver参数

4. study.solve()求解结果，并.write保存.sto运动轨迹数据predictSolution.sto，包含以下关节角度和速度状态量等条目/jointset/ankle_r/ankle_angle_r/value /jointset/knee_r/knee_angle_r/value /jointset/patellofemoral_r/knee_angle_r_beta/value /jointset/hip_r/hip_flexion_r/value /jointset/ankle_r/ankle_angle_r/speed /jointset/knee_r/knee_angle_r/speed /jointset/patellofemoral_r/knee_angle_r_beta/speed /jointset/hip_r/hip_flexion_r/speed /tau_hip_flexion_r /tau_knee_angle_r /tau_ankle_angle_r lambda_cid4_p0 gamma_cid4_p0

track(MocoStateTrackingGoal)_关节扭矩驱动的追踪问题，以上面生成运动轨迹为参考

1.TableProcessor加载predictSolution.sto数据并滤波

2. MocoStateTrackingGoal载入TableProcessor数据并添加到addGoal，.setWeight设置cost项的权重

3. 配置study的solver，使用setGuessFile()载入生成的predictSolution.sto运动轨迹数据，作为初始猜测值，生成trackingSolution.sto数据。包含条目与predictSolution相同，都是关节坐标状态。

inverse(.MocoInverse)_肌肉驱动的逆动力问题，求解肌肉力量

1. 创建MocoInverse实例化

2. ModelProcessor加载模型，模型的自由度提供额外的驱动力，并加载到inverse中

3. setKinematics加载TableProcessor 作为运动参考

4. 配置inverse时间点，网格数等边界条件，以及肌肉激活度为目标

5. solve，并将运动数据保存到inverseSolution.sto中，包含以下条目内容# activation，# normalized_tendon_force，# joint angle，joint speed，# muscle force,reserve jointset force,# implicitderiv_normalized_tendon_force，包含状态量，tendon_force，muscle_force。inverseSolution.sto条目如下

/forceset/bifemsh_r/activation /forceset/bifemsh_r/normalized_tendon_force /forceset/med_gas_r/activation /forceset/med_gas_r/normalized_tendon_force /forceset/glut_max2_r/activation /forceset/glut_max2_r/normalized_tendon_force /forceset/psoas_r/activation /forceset/psoas_r/normalized_tendon_force /forceset/rect_fem_r/activation /forceset/rect_fem_r/normalized_tendon_force /forceset/semimem_r/activation /forceset/semimem_r/normalized_tendon_force /forceset/soleus_r/activation /forceset/soleus_r/normalized_tendon_force /forceset/tib_ant_r/activation /forceset/tib_ant_r/normalized_tendon_force /forceset/vas_int_r/activation /forceset/vas_int_r/normalized_tendon_force /jointset/ankle_r/ankle_angle_r/value /jointset/ankle_r/ankle_angle_r/speed /jointset/knee_r/knee_angle_r/value /jointset/knee_r/knee_angle_r/speed /jointset/patellofemoral_r/knee_angle_r_beta/value /jointset/patellofemoral_r/knee_angle_r_beta/speed /jointset/hip_r/hip_flexion_r/value /jointset/hip_r/hip_flexion_r/speed /forceset/bifemsh_r /forceset/med_gas_r /forceset/glut_max2_r /forceset/psoas_r /forceset/rect_fem_r /forceset/semimem_r /forceset/soleus_r /forceset/tib_ant_r /forceset/vas_int_r /forceset/reserve_jointset_knee_r_knee_angle_r /forceset/reserve_jointset_ankle_r_ankle_angle_r /forceset/reserve_jointset_hip_r_hip_flexion_r lambda_cid4_p0 /forceset/bifemsh_r/implicitderiv_normalized_tendon_force /forceset/med_gas_r/implicitderiv_normalized_tendon_force /forceset/glut_max2_r/implicitderiv_normalized_tendon_force /forceset/psoas_r/implicitderiv_normalized_tendon_force /forceset/rect_fem_r/implicitderiv_normalized_tendon_force /forceset/semimem_r/implicitderiv_normalized_tendon_force /forceset/soleus_r/implicitderiv_normalized_tendon_force /forceset/tib_ant_r/implicitderiv_normalized_tendon_force /forceset/vas_int_r/implicitderiv_normalized_tendon_force

6. 使用STOFileAdapter工具，将inverseSolution导出到muscleOutputs.sto。muscleOutputs.sto文件包含以下每条肌肉的fiber_length和passive_force条目，单纯与肌肉相关。

/forceset/bifemsh_r|normalized_fiber_length /forceset/bifemsh_r|passive_force_multiplier /forceset/med_gas_r|normalized_fiber_length /forceset/med_gas_r|passive_force_multiplier /forceset/glut_max2_r|normalized_fiber_length /forceset/glut_max2_r|passive_force_multiplier /forceset/psoas_r|normalized_fiber_length /forceset/psoas_r|passive_force_multiplier /forceset/rect_fem_r|normalized_fiber_length /forceset/rect_fem_r|passive_force_multiplier /forceset/semimem_r|normalized_fiber_length /forceset/semimem_r|passive_force_multiplier /forceset/soleus_r|normalized_fiber_length /forceset/soleus_r|passive_force_multiplier /forceset/tib_ant_r|normalized_fiber_length /forceset/tib_ant_r|passive_force_multiplier /forceset/vas_int_r|normalized_fiber_length /forceset/vas_int_r|passive_force_multiplier

inverse_器械辅助下的肌肉驱动逆问题，求解器械辅助力

主要是修改模型，在膝盖坐标添加一个SpringGeneralizedForce弹簧阻尼模型。

1. 在knee_angle_r关节添弹簧元件，并添加到模型中

2. 模型中添加1个reserve actuators

3. inverse.solve求解，并保存inverseDeviceSolution.sto运动数据结果

4. getObjective获取最终cost，并画图，inverseSolution对比inverseDeviceSolution，观察肌肉激活度，肌腱力tendon_force的对比

# 第0部分：加载Moco库和opensim内置模型
import opensim as osim
import exampleSquatToStand_helpers as helpers
import mocoPlotTrajectory as plot
import os
import numpy as np
# 导入扭矩驱动模型和肌肉驱动模型
torqueDrivenModel = helpers.getTorqueDrivenModel()
muscleDrivenModel = helpers.getMuscleDrivenModel()

## Part 1: 求解关节扭矩驱动的预测问题，即由起始位置生成中间运动
# Part 1a: 创建MocoStudy类
study = osim.MocoStudy()

# Part 1b: 获取problem引用并加载模型
problem = study.updProblem()
problem.setModel(torqueDrivenModel)

# Part 1c: 边界设置
#
# 时间边界输入参数格式problem.setTimeBounds(initial_bounds, final_bounds)
# 状态量边界输出参数格式problem.setStateInfo(path, trajectory_bounds, inital_bounds, final_bounds)
#
# 所有边界参数可以设置为一个范围, [lower upper]； 或一个数值，表示上下界相等，或空的符号[]，表示默认值；
# 可以同时设置多个状态setStateInfoPattern()，格式如下:
# problem.setStateInfoPattern(pattern, trajectory_bounds, inital_bounds, ...
#       final_bounds)
#
# setStateInfoPattern函数支持“pattern”参数中的正则表达式；
# 使用'.*'匹配状态/控制路径的任何子字符串，例如，以下将设置所有坐标值状态信息：
# problem.setStateInfoPattern('/path/to/states/.*/value', ...)

# 时间边界设置[0, 1]，1秒之内从蹲到站起
problem.setTimeBounds(0, 1)

# 位置约束: 模型需要从蹲姿开始并以站姿结束
# 设定了了三个关节角度，
# 髋关节屈曲角度，膝关节角度，踝关节角度
# 关节角度上下界，初始位置上下界，和结束位置的上下界
# 髋关节屈曲角度，关节角度上下界[-114.59°, 28.65°]，初始值-114.59°，结束0°
problem.setStateInfo('/jointset/hip_r/hip_flexion_r/value', 
    [-2, 0.5], -2, 0)
problem.setStateInfo('/jointset/knee_r/knee_angle_r/value', 
    [-2, 0], -2, 0)
problem.setStateInfo('/jointset/ankle_r/ankle_angle_r/value', 
    [-0.5, 0.7], -0.5, -0.1)

# 速度界限：所有坐标 在开始和结束是都是静止
problem.setStateInfoPattern('/jointset/.*/speed', [], 0, 0)

# 添加控制量作为cost项，命名为myeffort
problem.addGoal(osim.MocoControlGoal('myeffort'))

# 配置求解器
# 前两行是必须的
solver = study.initCasADiSolver()
solver.set_num_mesh_intervals(25)
# 后面的选择项，可以不设置
solver.set_optim_convergence_tolerance(1e-4)
solver.set_optim_constraint_tolerance(1e-4)

# 如果'predictSolution.sto'这个文件不存在，执行求解并保存新的
if not os.path.isfile('predictSolution.sto'):
    # Part 1f: Solve! Write the solution to file, and visualize.
    predictSolution = study.solve()
    predictSolution.write('predictSolution.sto')
    # study.visualize(predictSolution)


## Part 2: 扭矩驱动的追踪问题，即给出运动轨迹参考
# 使用上述预测结果，使用TableProcessor加载数据，并低通滤波
# 将滤波数据，追加到tableProcessor中
tableProcessor = osim.TableProcessor('predictSolution.sto')
tableProcessor.append(osim.TabOpLowPassFilter(6))

# 使用MocoStateTrackingGoal将状态量作为Cost
# 生成的预测的TableProcessor作为参考轨迹
# 启用setAllowUnusedReference忽略不需要的列数据
tracking = osim.MocoStateTrackingGoal()
tracking.setName('mytracking')
tracking.setReference(tableProcessor)
tracking.setAllowUnusedReferences(True)
problem.addGoal(tracking)

# Part 2c: 使用track时，为了追踪轨迹的误差更小，一般会减小控制量惩罚项的权重
problem.updGoal('myeffort').setWeight(0.001)

# Part 2d: 使用predictSolution.sto作为初始猜测值
solver.setGuessFile('predictSolution.sto')
# 减小收敛公差，以确保平稳控制。
solver.set_optim_convergence_tolerance(1e-6)

# 如果'trackingSolution.sto'这个文件不存在，执行求解并保存新的
if not os.path.isfile('trackingSolution.sto'):
    # Part 2e: Solve! Write the solution to file, and visualize.
    trackingSolution = study.solve()
    trackingSolution.write('trackingSolution.sto')
    # study.visualize(trackingSolution)


## Part 3: 对比 Predictive与Tracking结果
# 使用绘图函数mocoPlotTrajectory.m，加载.sto运动数据文件
plot.mocoPlotTrajectory('predictSolution.sto', 'trackingSolution.sto', 
    'predict', 'track')

## Part 4: 肌肉驱动的逆动力问题，求解肌肉力量
# MocoInverse实例化
inverse = osim.MocoInverse()

# 使用ModelProcessor加载.osim肌肉骨骼模型，可以通过.append操作符来修改模型
modelProcessor = osim.ModelProcessor(muscleDrivenModel)
# 将备用执行器添加到模型中,给模型的自由度提供额外的驱动力或控制
# 添加model,operation,AddReserves
# 2为optimalForce，该自由度最大的力
modelProcessor.append(osim.ModOpAddReserves(2))
# 将配置完成的模型选项给到MocoInverse实例中	
inverse.setModel(modelProcessor)

# 使用TableProcessor 作为运动参考
inverse.setKinematics(tableProcessor)

# Part 4c: Set the time range, mesh interval, and convergence tolerance.
# 初始时间0s,终止时间1s,间隔为0.05s.
inverse.set_initial_time(0)
inverse.set_final_time(1)
inverse.set_mesh_interval(0.05)
#后面两行可以不设置
inverse.set_convergence_tolerance(1e-4)
inverse.set_constraint_tolerance(1e-4)

# 默认情况下，运动学数据包含额外的列，Moco会报错。
# 在这里，告诉Moco允许（并忽略）这些额外的列。
inverse.set_kinematics_allow_extra_columns(True)
# 设置肌肉激活程度最小
inverse.set_minimize_sum_squared_activations(True)

# 设定特定肌肉结果的关键字属性，设置muscleOutputs.sto文件的表头
inverse.append_output_paths('.*normalized_fiber_length')
inverse.append_output_paths('.*passive_force_multiplier')

# Part 4d: 求解! 将结果进行保存到 inverseSolution.sto文件中
# 内容是如下，每列个表头为一列数据
# activation，
# normalized_tendon_force，
# joint angle，joint speed，
# muscle force,reserve jointset force, 
# implicitderiv_normalized_tendon_force
inverseSolution = inverse.solve()
solution = inverseSolution.getMocoSolution()
solution.write('inverseSolution.sto')

# Part 4e: 将inverse结果写入到 muscleOutputs.sto中
# 内容是每个肌肉的fiber_length和passive_force_multiplier
inverseOutputs = inverseSolution.getOutputs()
sto = osim.STOFileAdapter()
sto.write(inverseOutputs, 'muscleOutputs.sto')

## Part 5: 器械辅助下的肌肉驱动逆问题，求解器械辅助力
# Part 5a: 修改肌肉驱动模型，膝盖坐标添加一个SpringGeneralizedForce。

# 在knee_angle_r广义坐标上添加一个弹簧阻尼组件，输入坐标名称
device = osim.SpringGeneralizedForce('knee_angle_r')
device.setStiffness(50)
device.setRestLength(0)
device.setViscosity(0)
# 在模型中添加力组件
muscleDrivenModel.addForce(device)

# Part 5b: Create a ModelProcessor similar to the previous one, using the same
# reserve actuator strength so we can compare muscle activity accurately.
# 修改后的模型，加载到ModelProcessor，添加储备执行器强度，比较是否器械情况下的肌肉情况
modelProcessor = osim.ModelProcessor(muscleDrivenModel)
modelProcessor.append(osim.ModOpAddReserves(2))
inverse.setModel(modelProcessor)

# Part 5c:固定格式，求解! 获取结果并保存
inverseDeviceSolution = inverse.solve()
deviceSolution = inverseDeviceSolution.getMocoSolution()
deviceSolution.write('inverseDeviceSolution.sto')

## Part 6: 比较无辅助和辅助情况下的逆问题结果
print('Cost without device: ', solution.getObjective())
print('Cost with device: ', deviceSolution.getObjective())

# 使用此函数进行方便的结果对比画图
helpers.compareInverseSolutions(inverseSolution, inverseDeviceSolution)

4. 关节扭矩驱动与肌肉驱动的模型处理方法

模型的预处理封装在exampleSquatToStand_helpers.py中供直接使用，此处只截取模型处理部分进行注释说明。


def addCoordinateActuator(model, coordName, optForce):
    """向给定的生物力学模型添加坐标致动器
    参数：
    - model:生物力学模型对象。
    - coordName:要添加致动器的坐标名称，要与.osim文件"Coordinate name=**"中一致
    - optForce:最大扭矩值。
    """
    # 获取模型的关节坐标集
    coordSet = model.updCoordinateSet()
    # 创建一个坐标致动器
    actu = osim.CoordinateActuator()
    # 设置致动器的名称
    actu.setName('tau_' + coordName)
    # setCoordinate将致动器关联到特定的坐标
    actu.setCoordinate(coordSet.get(coordName))
    # setOptimalForce设置致动器的最大扭矩
    actu.setOptimalForce(optForce)
    # 输入的control的最小/小控制值，一般设置inf,.osim中有对应属性
    actu.setMinControl(-1)
    actu.setMaxControl(1)
    # 将致动器添加到模型中
    model.addComponent(actu)

def getTorqueDrivenModel():
    # 加载基础模型
    model = osim.Model('squatToStand_3dof9musc.osim')

    # 移除模型中的肌肉
    model.updForceSet().clearAndDestroy()
    model.initSystem()

    # 向模型的关节自由度上添加坐标致动器
    addCoordinateActuator(model, 'hip_flexion_r', 150)
    addCoordinateActuator(model, 'knee_angle_r', 300)
    addCoordinateActuator(model, 'ankle_angle_r', 150)

    return model

def getMuscleDrivenModel():

    # 加载基础模型
    model = osim.Model('squatToStand_3dof9musc.osim')
    # 进行重载，完成模型的连接
    model.finalizeConnections()

    # 使用DeGrooteFregly提出的肌肉替换原模型，
    # 它们的特性曲线针对Direct Collocation法进行了优化（即无间断点、二次可微等）
    osim.DeGrooteFregly2016Muscle().replaceMuscles(model)

    # 通过强化模型和加宽主动力-长度曲线，使问题更容易求解
    # 遍历模型中的所有肌肉，对于每个肌肉有如下操作
    for m in np.arange(model.getMuscles().getSize()):
        musc = model.updMuscles().get(int(m))
        # 设置最小控制值为 0
        musc.setMinControl(0.0)
        # 设置不忽略激活动力学，即使用肌肉激活动力学和肌腱顺应性
        musc.set_ignore_activation_dynamics(False)
        musc.set_ignore_tendon_compliance(False)
        # 将最大等长力设置为原来的两倍，例子中只考虑身体一半的右腿
        musc.set_max_isometric_force(2.0 * musc.get_max_isometric_force())
        # 肌肉转换为
        dgf = osim.DeGrooteFregly2016Muscle.safeDownCast(musc)
        # 肌肉主动力-长度曲线（Muscle Active Force-Length Curve）中
        # 设置肌肉主动力宽度为 1.5倍
        dgf.set_active_force_width_scale(1.5)
        # 肌腱的顺应性动力学模型设置为隐式形式
        dgf.set_tendon_compliance_dynamics_mode('implicit')

        # 忽略比目鱼肌“soleus_r”被动纤维力，被动纤维力对应具有刚性
        # 因为比目鱼肌肌腱很长，比目鱼肌肌腱不具有被动纤维力，即没有刚性，为纯柔性
        # 否则，因为纤维长度，会产生过多的被动纤维力
        # 设置后，soleus_r_passive_force列的值都为0
        if str(musc.getName()) == 'soleus_r':
            dgf.set_ignore_passive_fiber_force(True)

    return model