Transition perdiod coordinate limit force

[pauline neerpelt]

Hi, I have a question about the transition period of a coordinate limit force. How is the transition from zero to constant stiffness described? I want to know the stiffness at a specific coordinate.
Any help is much appreciated.

[Thomas Uchida]

Please see the description and plot in the doxygen documentation here: https://simtk.org/api_docs/opensim/api_ ... ml#details. The calculations can be found in the CoordinateLimitForce::calcLimitForce() method (https://github.com/opensim-org/opensim- ... e.cpp#L260).

[pauline neerpelt]

Thank you very much for your reply!

I have attempted to use the source code to write a MATLAB script for determining the CLF at a prescribed angle (attached below).
I get the following results:

CLF.png (36.13 KiB) Viewed 278 times

I'm uncertain whether I have correclty implemented the simmspline function and the CLF transitions and would very much appreciate any feedback.

Kind regards,
Pauline

Code: Select all

import org.opensim.modeling.*;

%CLF parameters for the ARMS Lab Hand and Wrist model (McFarland et al. 2023)
%FLEXION-EXTENSION DOF
%Units: deg and Nm/deg
upperLimit = 60;              
lowerLimit = -60;             
upperStiffness = 200;         
lowerStiffness = 135.3232;    
transitionRange = 92.916;     

% Calculate limit force at angle = -10 degrees
limitForce = calcLimitForce(upperLimit, lowerLimit, upperStiffness, lowerStiffness, transitionRange);
disp(['CLF at angle = -10: ', num2str(limitForce)]);

%GRAPH
angleRange = -200:0.5:200;
limitForces = zeros(size(angleRange));

%Simmsplinbe in transition regions for lower and upper stiffness
%upper
upperSpline = SimmSpline();
upperSpline.addPoint(upperLimit, 0);
upperSpline.addPoint(upperLimit + transitionRange, upperStiffness);
%lower
lowerSpline = SimmSpline();
lowerSpline.addPoint(lowerLimit - transitionRange, lowerStiffness);
lowerSpline.addPoint(lowerLimit, 0);

% Calculate limit forces for each angle in the range
for i = 1:length(angleRange)
    angle = angleRange(i);
    
    % Determine the current stiffness using SimmSpline functions
    upperK = calculateUpperStiffness(angle, upperSpline, upperLimit, upperStiffness, transitionRange);
    lowerK = calculateLowerStiffness(angle, lowerSpline, lowerLimit, lowerStiffness, transitionRange);
    
    % Calculate the total limit force
    upperForce = -upperK * (angle - upperLimit);
    lowerForce = lowerK * (lowerLimit - angle);
    limitForces(i) = upperForce + lowerForce;
end

% Plot the results
figure;
plot(angleRange, limitForces, 'k', 'LineWidth', 2);
grid on;
xlabel('Coordinate Position (deg)');
ylabel('Coordinate Limit Force (Nm)');
title('Coordinate Limit Force for Extension-Flexion DOF');

% Mark the important lines
% ROM for flexon-extension
xline(-70, '--', 'LineWidth', 1);
xline(70, '--', 'LineWidth', 1);
% show upper and lower limit
xline(lowerLimit, 'r', 'LineWidth', 1);
xline(upperLimit, 'r', 'LineWidth', 1);
% show angle at which upper and lower stiffness is reached
xline(lowerLimit - transitionRange, 'b', 'LineWidth', 1);
xline(upperLimit + transitionRange, 'b', 'LineWidth', 1);

%Limits in flexion-extension = 70-70
limitForceAtNeg70 = interp1(angleRange, limitForces, -70);
limitForceAt70 = interp1(angleRange, limitForces, 70);
text(-70, limitForceAtNeg70, sprintf('%.2f', limitForceAtNeg70), ...
    'HorizontalAlignment', 'right', 'VerticalAlignment', 'bottom', 'FontSize', 12);
text(70, limitForceAt70, sprintf('%.2f', limitForceAt70), ...
    'HorizontalAlignment', 'left', 'VerticalAlignment', 'bottom', 'FontSize', 12);


% Calculate CLF at a specific angle
function limitForce = calcLimitForce(upperLimit, lowerLimit, upperStiffness, lowerStiffness, transitionRange)
    import org.opensim.modeling.*;
    angle = -10;
    
    upperSpline = SimmSpline();
    upperSpline.addPoint(upperLimit, 0);
    upperSpline.addPoint(upperLimit + transitionRange, upperStiffness);

    lowerSpline = SimmSpline();
    lowerSpline.addPoint(lowerLimit - transitionRange, lowerStiffness);
    lowerSpline.addPoint(lowerLimit, 0);
    
    %Stiffness
    upperK = calculateUpperStiffness(angle, upperSpline, upperLimit, upperStiffness, transitionRange);
    lowerK = calculateLowerStiffness(angle, lowerSpline, lowerLimit, lowerStiffness, transitionRange);
    
    %Force
    upperForce = -upperK * (angle - upperLimit);
    lowerForce = lowerK * (lowerLimit - angle);
    limitForce = upperForce + lowerForce;
end

% Calculate upper stiffness using SimmSpline with constant behavior beyond the transition zone
function upperK = calculateUpperStiffness(angle, upperSpline, upperLimit, upperStiffness, transitionRange)
    import org.opensim.modeling.*;
    if angle < upperLimit
        upperK = 0; % Zero stiffness below the upper limit
    elseif angle > (upperLimit + transitionRange)
        upperK = upperStiffness; % Constant stiffness beyond the transition range
    else
        upperK = upperSpline.calcValue(Vector(1, angle));
    end
end

% Calculate lower stiffness using SimmSpline with constant behavior beyond the transition zone
function lowerK = calculateLowerStiffness(angle, lowerSpline, lowerLimit, lowerStiffness, transitionRange)
    import org.opensim.modeling.*;
    if angle > lowerLimit
        lowerK = 0; % Zero stiffness above the lower limit
    elseif angle < (lowerLimit - transitionRange)
        lowerK = lowerStiffness; % Constant stiffness beyond the transition range
    else
        lowerK = lowerSpline.calcValue(Vector(1, angle));
    end
end