Sensitivity Analysis

[Jeremy Genter]

Dear Community,

I want to do a sensitivity analysis in OpenSim 3.3. What I want to find out is how sensible the location of the insertion sites are. I want to presribe the muscle forces (not the muscle activations) in a forward dynamics manner. Then repeat this for for a grid of different locations of the insertion sites. How can I implement this in an efficient way through the MATLAB API?

Thank you for any suggestions.

Kind regards,
JG

[Jeremy Genter]

After some time of digging I coded the script below (snipped of the code). I use Wu's Shoulder model https://simtk.org/projects/wu-shoulder to calculate the muscle forces. I try to do a forward simulation to see if I get the same motion if I impose the muscle forces directly (Fig. 1). I get a bit better results if I impose the activation (Fig.2) I expect a -30degrees (as i did Wu's muscle optimization with this abduction) of abduction and other DOF should remain 0. Can somebody help me to find out what I did wrong? I appreciate any help.

Kind Regards,
JG

Fig1:
[

imposed forces.jpg (66.89 KiB) Viewed 349 times

Fig2:

imposed activations.jpg (62.58 KiB) Viewed 349 times

Code: Select all

import org.opensim.modeling.*
% initialize model, parameters, etc,
my_model = Model('Shoulder_Simulator.osim')
dt = 0.2;
strart_t = 0.5;
final_t = 6.5;
t = (strart_t: dt : final_t)';
data = readmatrix([folders.results_folder_muscle_opt,...
                  	'simulator_data\abduction_simulator_muscle_force.csv']); % muscle forces calculated with https://simtk.org/projects/wu-shoulder
act_m = interp1(data(:,1), data(:,2:end), t);

state = my_model.initSystem();

% Get the model's joint set
jointSet = my_model.getJointSet();
joint = jointSet.get(0);  
coordSet = joint.getCoordinateSet();
% Fix some degrees of freedom (DOF)
DOF2fix = {'Model_Positioning','Model_Positioning_yRotation','Model_Positioning_zRotation'...
           'Model_Positioning_xTranslation', 'Model_Positioning_yTranslation', 'Model_Positioning_zTranslation'};
for j = 1:length(DOF2fix)
    coordSet.get(DOF2fix{j}).setLocked(state, true);
end
joint = jointSet.get(1);  
coordSet = joint.getCoordinateSet();
coordSet.get('Scapula_rotation').setLocked(state, true);

% Get the ForceSet
forces = my_model.getForceSet();

% Run the forward simulation
manager = Manager(my_model);
% integrator = SimTK.RungeKuttaMersonIntegrator();
% manager.setIntegrator(integrator);
% manager.initialize(state,);
tic
GH_mot = zeros(length(t), 3);

    for i = 1:length(t)
        % Set the muscle activations for this time step

        for j = 1:nMusc
            muscle = Thelen2003Muscle.safeDownCast(my_musc_set.get(j-1));
            muscle.setIgnoreTendonCompliance(state, true);
%             muscle.setActivation(state, act_m(i,j))
            muscle.overrideForce(state, true)
            muscle.setForce(state, act_m(i,j));            
        end
        my_model.equilibrateMuscles(state);
        % Integrate the model for this time step
        manager.integrate(state, dt);
        
        jointSet = my_model.getJointSet();
        joint = jointSet.get(2); % get the GH joint
        for j = 1:3
            coordinateSet = joint.getCoordinateSet(); % get the set of coordinates for the joint
            coordinate = coordinateSet.get(j-1); % get the coordinates (abduction, flexion, int. rot.)
            GH_mot(i,j) = coordinate.getValue(state); % get the value of the coordinate from the state
        end
    end

toc
figure
plot(t-t(1), GH_mot*180/pi)
xlabel('time (s)')
ylabel('joint motion (derees)')
legend('abduction', 'flexion', 'internal rotation')

[Jeremy Genter]

I found an error:
I should have used

Code: Select all

muscle.setOverrideForce(state, act_m(i,j));

instead of

Code: Select all

muscle.setForce(state, act_m(i,j));

Further, in muscle optimization I accidendently solved for an adduction motion and not an abduction.
Still, I get oscillating solution.
How do I define the integrator method and step size?
I tried to set

Code: Select all

 manger.setIntegratorMethod(RungeKuttaMersonIntegrator())

it didn't work in MATLAB