Page 1 of 1
Sensitivity Analysis
Posted: Tue Jan 03, 2023 7:25 am
by gent
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
Re: Sensitivity Analysis
Posted: Tue Feb 07, 2023 4:24 am
by gent
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 350 times
Fig2:
- imposed activations.jpg (62.58 KiB) Viewed 350 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')
Re: Sensitivity Analysis
Posted: Tue Feb 07, 2023 8:13 am
by gent
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