Can we optimize parameter of FiberForceLengthCurve with MOCO

[Sepehr Ramezani]

Hi,
I'm working with MOCOParameter and interest in optimizing some parameters (e.g. "strain_at_one_norm_force") of the FiberForceLengthCurve of Millard2012EquilibriumMuscle. But can not add parameter to problem. I used following code and did not work.

Code: Select all

param = MocoParameter('strain_at_norm', '/forceset/bflh_r', 'strain_at_one_norm_force', MocoBounds(0.5,0.9));
problem.addParameter(param);

I would appreciate it if anyone can help me.
Best,

[Nicholas Bianco]

Hi Sepehr,

It looks like the issue is that provide the path to the Millard2012EquilibriumMuscle object for the component path argument, but the 'strain_at_one_norm_force' is a property of the FiberForceLengthCurve object, which is a subcomponent of the muscle. Provide the full path to that object and it should work.

-Nick

[Sepehr Ramezani]

Hi Nick,
That was smart, I have tried below code but did not work, Could you please explain it how can I provide full path?
Thank you,

Code: Select all

param = MocoParameter('strain_at_one_norm_force', '/forceset/bflh_r/bflh_r_FiberForceLengthCurve', 'strain_at_one_norm_force', MocoBounds(0.5,0.9));

[Nicholas Bianco]

Hi Sepehr,

That looks like the correct path. You can try calling the following method on Model to print out model paths:

Code: Select all

model.printSubcomponentInfo();

What error message are you seeing?

-Nick

[Sepehr Ramezani]

Hi
I tried this

Code: Select all

frcset = model.getComponent('/forceset/bflh_r')
frcset.printSubcomponentInfo()

and got this result:

Code: Select all

Millard2012EquilibriumMuscle  /forceset/bflh_r
            [MuscleFixedWidthPennationModel]          /penMdl
    [MuscleFirstOrderActivationDynamicModel]          /actMdl
                              [GeometryPath]          /geometrypath
                                 [PathPoint]              /bflh_r-P1
                                   [Station]                  /station
                                 [PathPoint]              /bflh_r-P2
                                   [Station]                  /station
                                 [PathPoint]              /bflh_r-P3
                                   [Station]                  /station

I can add any variables under these components. e.g.

Code: Select all

param = MocoParameter('Test', '/forceset/bflh_r/geometrypath/bflh_r-P1/','location', MocoBounds(0,1),1);

However I could not find subcomponent with name of bflh_r_FiberForceLengthCurve while I have it in xml (below code) file. Do you have any idea?

Code: Select all

<Millard2012EquilibriumMuscle name="bflh_r">
					<!--Minimum allowed value for control signal. Used primarily when solving for control values.-->
					<min_control>0.01</min_control>
					<!--Maximum allowed value for control signal. Used primarily when solving for control values.-->
					<max_control>0.01</max_control>
					<!--The set of points defining the path of the actuator.-->
					<GeometryPath name="geometrypath">
						<!--The set of points defining the path-->
						<PathPointSet>
							<objects>
								<PathPoint name="bflh_r-P1">
									<!--Path to a Component that satisfies the Socket 'parent_frame' of type PhysicalFrame (description: The frame in which this path point is defined.).-->
									<socket_parent_frame>/bodyset/pelvis</socket_parent_frame>
									<!--The fixed location of the path point expressed in its parent frame.-->
									<location>-0.129442 -0.119299 0.070165199999999997</location>
								</PathPoint>
								<PathPoint name="bflh_r-P2">
									<!--Path to a Component that satisfies the Socket 'parent_frame' of type PhysicalFrame (description: The frame in which this path point is defined.).-->
									<socket_parent_frame>/bodyset/tibia_r</socket_parent_frame>
									<!--The fixed location of the path point expressed in its parent frame.-->
									<location>-0.0387213 -0.040215000000000001 0.0290697</location>
								</PathPoint>
								<PathPoint name="bflh_r-P3">
									<!--Path to a Component that satisfies the Socket 'parent_frame' of type PhysicalFrame (description: The frame in which this path point is defined.).-->
									<socket_parent_frame>/bodyset/tibia_r</socket_parent_frame>
									<!--The fixed location of the path point expressed in its parent frame.-->
									<location>-0.0329763 -0.0522795 0.034814699999999997</location>
								</PathPoint>
							</objects>
							<groups />
						</PathPointSet>
						<!--The wrap objects that are associated with this path-->
						<PathWrapSet>
							<objects />
							<groups />
						</PathWrapSet>
						<!--Default appearance attributes for this GeometryPath-->
						<Appearance>
							<!--The color, (red, green, blue), [0, 1], used to display the geometry. -->
							<color>0.80000000000000004 0.10000000000000001 0.10000000000000001</color>
						</Appearance>
					</GeometryPath>
					<!--Maximum isometric force that the fibers can generate-->
					<max_isometric_force>3000</max_isometric_force>
					<!--Optimal length of the muscle fibers-->
					<optimal_fiber_length>0.095423233460067203</optimal_fiber_length>
					<!--Resting length of the tendon-->
					<tendon_slack_length>0.31797313845694503</tendon_slack_length>
					<!--Angle between tendon and fibers at optimal fiber length expressed in radians-->
					<pennation_angle_at_optimal>0.17591823000000001</pennation_angle_at_optimal>
					<!--Maximum contraction velocity of the fibers, in optimal fiberlengths/second-->
					<max_contraction_velocity>15</max_contraction_velocity>
					<!--Compute muscle dynamics ignoring tendon compliance. Tendon is assumed to be rigid.-->
					<ignore_tendon_compliance>true</ignore_tendon_compliance>
					<!--Compute muscle dynamics ignoring activation dynamics. Activation is equivalent to excitation.-->
					<ignore_activation_dynamics>true</ignore_activation_dynamics>
					<!--The linear damping of the fiber.-->
					<fiber_damping>0.01</fiber_damping>
					<!--Assumed initial activation level if none is assigned.-->
					<default_activation>0.01</default_activation>
					<!--Activation lower bound.-->
					<minimum_activation>0.01</minimum_activation>
					<!--Active-force-length curve.-->
					<ActiveForceLengthCurve name="bflh_r_ActiveForceLengthCurve">
						<!--Normalized fiber length where the steep ascending limb starts-->
						<min_norm_active_fiber_length>0.25</min_norm_active_fiber_length>
						<!--Normalized fiber length where the steep ascending limb transitions to the shallow ascending limb-->
						<transition_norm_fiber_length>0.77000000000000002</transition_norm_fiber_length>
						<!--Normalized fiber length where the descending limb ends-->
						<max_norm_active_fiber_length>1.8999999999999999</max_norm_active_fiber_length>
						<!--Slope of the shallow ascending limb-->
						<shallow_ascending_slope>0.75</shallow_ascending_slope>
						<!--Minimum value of the active-force-length curve-->
						<minimum_value>0</minimum_value>
					</ActiveForceLengthCurve>
					<!--Force-velocity curve.-->
					<ForceVelocityCurve name="bflh_r_ForceVelocityCurve">
						<!--Curve slope at the maximum normalized concentric (shortening) velocity (normalized velocity of -1)-->
						<concentric_slope_at_vmax>0</concentric_slope_at_vmax>
						<!--Curve slope just before reaching concentric_slope_at_vmax-->
						<concentric_slope_near_vmax>0.25</concentric_slope_near_vmax>
						<!--Curve slope at isometric (normalized velocity of 0)-->
						<isometric_slope>5</isometric_slope>
						<!--Curve slope at the maximum normalized eccentric (lengthening) velocity (normalized velocity of 1)-->
						<eccentric_slope_at_vmax>0</eccentric_slope_at_vmax>
						<!--Curve slope just before reaching eccentric_slope_at_vmax-->
						<eccentric_slope_near_vmax>0.14999999999999999</eccentric_slope_near_vmax>
						<!--Curve value at the maximum normalized eccentric contraction velocity-->
						<max_eccentric_velocity_force_multiplier>1.3999999999999999</max_eccentric_velocity_force_multiplier>
					</ForceVelocityCurve>
					<!--Passive-force-length curve.-->
					<FiberForceLengthCurve name="bflh_r_FiberForceLengthCurve">
						<!--Fiber strain at zero force-->
						<strain_at_zero_force>0</strain_at_zero_force>
						<!--Fiber strain at a tension of 1 normalized force-->
						<strain_at_one_norm_force>0.69999999999999996</strain_at_one_norm_force>
					</FiberForceLengthCurve>
					<!--Tendon-force-length curve.-->
					<TendonForceLengthCurve name="bflh_r_TendonForceLengthCurve">
						<!--All properties of this object have their default values.-->
					</TendonForceLengthCurve>
				</Millard2012EquilibriumMuscle>

Best

[Nicholas Bianco]

Hi Sepehr,

After looking into this more closely, the FiberForceLengthCurve is stored in the Millard muscle object as a property. It is not a subcomponent of the Millard muscle, that's why it's not showing up on the path. So unfortunately, I don't think you can optimize the this curve with the current MocoParameter interface.

If you use the DeGrooteFregly2016Muscle, it has a property called 'passive_fiber_strain_at_one_norm_force'. If you're trying to optimize passive forces and don't mind switching to that muscle, that could be viable alternative.

-Nick

[Sepehr Ramezani]

Hi Nick,
Thank you so much.
I'll try it, last time when I compared passive force of Millard and Degroot on my study, there was a large differences (up to 45%). I believe Degroot's parameter are optimized for walking or behavior close to walking (in the same range of motion) and the more you get far from that situation the more error you get. I think it happens in my study because I cover a larger range of motion.
I will share the result with you once I try it.
Thanks!

[Sepehr Ramezani]

If you use the DeGrooteFregly2016Muscle, it has a property called 'passive_fiber_strain_at_one_norm_force'. If you're trying to optimize passive forces and don't mind switching to that muscle, that could be viable alternative.

Hi,
quick update for this post it might be useful. I have used DegrooteMuscle and successfully optimized 2 parameters of tendon_slack_length and passive_fiber_strain_at_one_norm_force by the error of 12%. There were some issues related to this model which needs to be more careful using this model. here and here (I believe scaling issues are already solved and would be added to 4.4, this problem also caused a huge passive fiber force for tfl muscle after scaling and I had to remove this muscle).
I decided to add another parameter (tendon_strain_at_one_norm_force), simulation converged however, regardless of using different data set, Moco always gives me the same value for this parameter for all muscles which is average of the lower and higher bound!
Is it related to the CasADiSolver? appreciate any suggestions.
Best,
Sepehr

[Nicholas Bianco]

Hi Sepehr,

What is your problem formulation? It could be that the initial value is close to optimal or that the objective gradient with respect to this parameter is close to zero.

Best,
Nick

[Sepehr Ramezani]

Hi Nick,
My case is: parameter estimation of knee's flexors and extensors during Isokinetic knee flexion (0deg to 90deg).
I'm using MOCOTrack, stateTracking for motion, plus control tracking (external force).
estimation of these 2 parameters was successful:
tendon_slack_length
passive_fiber_strain_at_one_norm_force
However, when I try to estimate tendon_strain_at_one_norm_force and MOCO gives a fix number which is equal to the average of the lower and upper of boundary as the output of estimation. I don't think it's because the value is close to optimal or objective gradient with respect to this parameter is close to zero. I tried different scenarios and all the times it gives me same exact number. it seems MOCO cant not play with this parameter during optimization and just gives the initial value which is the average of boundaries.
Best,
Sepehr