ForwardTool failing with PathActuator

[Caitlin Romanczyk]

Hello,

I've been using the ForwardTool with the MoBL-ARMS Dynamic Upper Limb Model to generate a simulation of tremor, and it's been running fine. However, I've just begun developing an actuator framework to model an exoskeleton, and now with the additional actuators in the model, the ForwardTool fails (I'm just setting the controls to zero for simplicity right now). I'm not seeing any error messages in the out log--it just stops before completing the simulation. I'm very new to OpenSim, so it may be something fairly basic that I'm missing--if anyone has an idea of what it might be, I'd appreciate it!

I'd also appreciate if anyone could give me advice on whether to create separate controllers for the muscles and the actuators in the model, or if it's better to use a single controller that's broken down with if/else statements.

My code is copied below for reference.

Thanks!
Caitlin

Code: Select all

#include <OpenSim\OpenSim.h>
#include <iostream>
#include <string>
#include <iostream>
#include <fstream>
#include <string>
#include "PathActuator.h"

using namespace OpenSim;
using namespace SimTK;
using namespace std;


int main()
{
	try{

		// load existing OpenSim model
		Model wristTremorModel("wristModel.osim");
		cout << "Done loading existing model" << endl << endl;
		// get existing objects
		const BodySet &bodySet = wristTremorModel.getBodySet();
		int numBodies = wristTremorModel.getNumBodies();
		const JointSet &jointSet = wristTremorModel.getJointSet();
		int numJoints = wristTremorModel.getNumJoints();
		const CoordinateSet &coordinateSet = wristTremorModel.getCoordinateSet();
		const ForceSet &forceSet = wristTremorModel.getForceSet();
		const Set<Muscle> &muscleSet = wristTremorModel.getMuscles();
		const int numMuscles = muscleSet.getSize();
		cout << "No. of muscles: " << numMuscles << endl << endl;

		/////////////////////////////////
		// ADD ACTUATORS TO OSIM MODEL //
		/////////////////////////////////
		
		//Get reference to bodies for construction of path actuators
		OpenSim::Body& proximal_anchor = wristTremorModel.updBodySet().get("proximal_anchor");
		OpenSim::Body& distal_anchor = wristTremorModel.updBodySet().get("distal_anchor");
		OpenSim::Body& proximal_row = wristTremorModel.updBodySet().get("proximal_row");
		
		//Get reference to WrapCylinder
		WrapObject *wristCylinderRef = proximal_row.getWrapObject("wristCylinder");

		//Construct agonist path actuator and define path geometry
		PathActuator *agonistAct = new PathActuator();
		agonistAct->setName("agonistAct");
		agonistAct->setAllPropertiesUseDefault(true);
		agonistAct->setOptimalForce(10.0);
		agonistAct->setMinControl(0.0);
		agonistAct->setMaxControl(10.0);
		agonistAct->addNewPathPoint("distal_anchor_P1", distal_anchor, Vec3(0, 0, 0.015));
		agonistAct->addNewPathPoint("proximal_anchor_P1", proximal_anchor, Vec3(0, 0, 0.04));
		wristTremorModel.addForce(agonistAct);
		cout << "Added agonist actuator to model" << endl;
		agonistAct->updGeometryPath().addPathWrap(*wristCylinderRef);
		cout << "Added path wrap" << endl;

		//Construct antagonist path actuator and define path geometry
		PathActuator *antagonistAct = new PathActuator();
		antagonistAct->setName("antagonistAct");
		antagonistAct->setAllPropertiesUseDefault(true);
		antagonistAct->setOptimalForce(10.0);
		antagonistAct->setMinControl(0.0);
		antagonistAct->setMaxControl(10.0);
		antagonistAct->addNewPathPoint("distal_anchor_P1", distal_anchor, Vec3(0, 0, -0.015));
		antagonistAct->addNewPathPoint("proximal_anchor_P1", proximal_anchor, Vec3(0, 0, -0.04));
		wristTremorModel.addForce(antagonistAct);
		cout << "Added antagonist actuator to model" << endl;
		antagonistAct->updGeometryPath().addPathWrap(*wristCylinderRef);
		cout << "Added path wrap" << endl;

		
		// Initialize the system and get the states representing the state system
		State &states = wristTremorModel.initSystem();
		// initialize position of upper arm to 0 deg
		for (int i = 0; i < coordinateSet.getSize(); i++)
		{
			cout << "Coordinate: " << coordinateSet.get(i).getName() << ", Value: "
				<< wristTremorModel.updCoordinateSet()[i].getValue(states) << endl;

			if ((coordinateSet.get(i).getName()) == "shoulder_elv")
			{
				coordinateSet.get(i).setDefaultValue(0);
			}
			cout << "Coordinate: " << coordinateSet.get(i).getName() << ", Value: "
				<< wristTremorModel.updCoordinateSet()[i].getValue(states) << endl;
		}
		
		// introducing variable - muscle (dynamic array of pointers to muscle)
		Schutte1993Muscle_Deprecated** muscle = new Schutte1993Muscle_Deprecated*[numMuscles];

		

		///////////////////////////////////////////////////
		// ADD CONTROLLER TO INTRODUCE TREMOR IN MUSCLES //
		///////////////////////////////////////////////////

		// remove pre-existing controller
		double numExistingControllers = wristTremorModel.getControllerSet().getSize();
		if (numExistingControllers > 0)
		{
			Controller &existingController = wristTremorModel.getControllerSet().get(0);
			wristTremorModel.removeController(&existingController);
		}

		// define controller to introduce tremor in muscles
		PrescribedController *tremorController = new PrescribedController();
		tremorController->setName("tremor");
		tremorController->setActuators(wristTremorModel.updActuators());
		cout << "setup controller" << endl << endl;

		// define arrays for tremor function data
		const int ARRAY_SIZE = 15000;
		double timeArray[ARRAY_SIZE];
		double activationArray1[ARRAY_SIZE];
		double activationArray2[ARRAY_SIZE];
		int count;

		// get timeArray
		count = 0;
		ifstream timeFile;
		timeFile.open("C:\\Users\\Katy\\Documents\\MATLAB\\TimeArray.txt");
		if (timeFile.is_open())
		{
			cout << "timeFile is open" << "\n";
			while (count < ARRAY_SIZE && timeFile >> timeArray[count])
			{
				count++;
			}
			timeFile.close();
			if (!(timeFile.is_open()))
			{
				cout << "timeFile has been closed" << "\n";
			}
		}
		else
			std::cout << "timeFile is not open" << "\n";

		// get activationArray1
		count = 0;
		ifstream agonistFile;
		agonistFile.open("C:\\Users\\Katy\\Documents\\MATLAB\\AgonistActivation.txt");
		if (agonistFile.is_open())
		{
			cout << "agonistFile is open" << "\n";
			while (count < ARRAY_SIZE && agonistFile >> activationArray1[count])
			{
				count++;
			}
			agonistFile.close();
			if (!(agonistFile.is_open()))
			{
				cout << "agonistFile has been closed" << "\n";
			}
		}
		else
			std::cout << "agonistFile is not open" << "\n";

		// get activationArray2
		count = 0;
		ifstream antagonistFile;
		antagonistFile.open("C:\\Users\\Katy\\Documents\\MATLAB\\AntagonistActivation.txt");
		if (antagonistFile.is_open())
		{
			cout << "antagonistFile is open" << "\n";
			while (count < ARRAY_SIZE && antagonistFile >> activationArray2[count])
			{
				count++;
			}
			antagonistFile.close();
			if (!(antagonistFile.is_open()))
			{
				cout << "antagonistFile has been closed" << "\n";
			}
		}
		else
			std::cout << "antagonistFile is not open" << "\n";

		// muscle activation times
		double activation = 0.0001, deactivation = 1.0;

		// interpolate activation arrays using piecewise linear functions
		PiecewiseLinearFunction *TF1 = new PiecewiseLinearFunction(15000, timeArray, activationArray1);
		PiecewiseLinearFunction *TF2 = new PiecewiseLinearFunction(15000, timeArray, activationArray2);

		//Constant *TF1 = new Constant(activationArray1[0]);
		//Constant *TF2 = new Constant(activationArray2[0]);
		cout << "Done initializing controller functions" << endl << endl;

		// set initial activation
		for (int i = 0; i < numMuscles; i++)
		{
			cout << i << endl;
			// downcast muscles to type Schuette1993model if they were not originally of that type
			muscle[i] = Schutte1993Muscle_Deprecated::safeDownCast(&muscleSet[i]);
			cout << "Muscle Name: " << muscle[i]->getName() << endl
				<< "Fiber Lengths: " << muscle[i]->getFiberLength(states) << endl
				<< "Optimal Fiber Length: " << muscle[i]->getOptimalFiberLength() << endl
				<< "Muscle Activation: " << muscle[i]->getDefaultActivation() << endl << endl;
		}
		cout << "muscle downcast to Schutte1993Muscle_Deprecated" << endl << endl;

		// wristTremorModel.equilibrateMuscles(states);
		// cout << "muscles in equilibrium" << endl << endl;

		// assigning muscles with the tremor neural activation
		for (int i = 0; i < numMuscles; i++)
		{
			cout << "Controller No. " << i << endl;
			// flexion muscles of wrist
			if (muscle[i]->getName() == "FCR" || muscle[i]->getName() == "FCU")
			{
				//muscle[i]->setActivationTimeConstant(activation);
				//muscle[i]->setDeactivationTimeConstant(deactivation);
				muscle[i]->setDefaultActivation(activationArray1[0]);
				muscle[i]->setDefaultFiberLength(0.1);
				tremorController->prescribeControlForActuator(muscle[i]->getName(), TF1);
				cout << "Muscle Name: " << muscle[i]->getName() << endl
					<< "Fiber Lengths: " << muscle[i]->getFiberLength(states) << endl
					<< "Optimal Fiber Length: " << muscle[i]->getOptimalFiberLength() << endl
					<< "Muscle Activation: " << muscle[i]->getDefaultActivation() << endl << endl;
			}
			// extension muscles of wrist
			else if (muscle[i]->getName() == "ECRL" || muscle[i]->getName() == "ECRB")
			{
				//muscle[i]->setActivationTimeConstant(activation);
				//muscle[i]->setDeactivationTimeConstant(deactivation);
				muscle[i]->setDefaultActivation(activationArray2[0]);
				muscle[i]->setDefaultFiberLength(0.1);
				cout << "Muscle Name: " << muscle[i]->getName() << endl
					<< "Fiber Lengths: " << muscle[i]->getFiberLength(states) << endl
					<< "Muscle Activation: " << muscle[i]->getDefaultActivation() << endl << endl;
				tremorController->prescribeControlForActuator(muscle[i]->getName(), TF2);
			}
			// set all other muscles to no activation
			else
			{
				muscle[i]->setDefaultActivation(0.0);
				//muscle[i]->setActivationTimeConstant(activation);
				//muscle[i]->setDeactivationTimeConstant(deactivation);
				tremorController->prescribeControlForActuator(muscle[i]->getName(), new Constant(0.0));
			}
		}
		cout << "Controller for individual muscles" << endl << endl;

		// add controller
		wristTremorModel.addController(tremorController);
		cout << "Done adding controller" << endl << endl;
		// save model to file
		wristTremorModel.print("wristModel.osim");
		cout << "Done printing to file" << endl << endl;

		/////////////////////////////////////
		// ADD CONTROLLER TO PATH ACTUATOR //
		/////////////////////////////////////
		/*
		// define controller to introduce tremor in muscles
		PrescribedController *pathController = new PrescribedController();
		pathController->setName("pathController");
		pathController->setActuators(agonistAct,antagonistAct);
		cout << "setup path controller" << endl << endl;
		pathController->prescribeControlForActuator("agonistAct", new Constant(0.0));
		pathController->prescribeControlForActuator("antagonistAct", new Constant(0.0));
		// add controller
		wristTremorModel.addController(pathController);
		cout << "Done adding controller" << endl << endl;
		// save model to file
		wristTremorModel.print("wristModel.osim");
		cout << "Done printing to file" << endl << endl;
		*/

		////////////////
		// SIMULATION //
		////////////////

		// set visualizer
		wristTremorModel.setUseVisualizer(false);

		// compute initial conditions for muscles
		wristTremorModel.equilibrateMuscles(states);

		// start and finish times for simulation
		double startTime = 0.0, finishTime = 5.0;

		//forward dynamics
		ForwardTool *tool = new ForwardTool();
		tool->setModel(wristTremorModel);
		tool->setStartTime(startTime);
		tool->setFinalTime(finishTime);
		tool->setSolveForEquilibrium(true);
		tool->setName("wristModel_Trial");
		tool->run();

[Thomas Uchida]

I read through your code until I found an issue:

Code: Select all

// This line
coordinateSet.get(i).setDefaultValue(0);
// Should be
coordinateSet.get(i).setValue(states, 0);

Setting the default_value property after calling initSystem() will, at best, have no effect. To debug further, you may want to add only one muscle since typos between "agonist" and "antagonist" might be difficult to catch. A good general debugging strategy is to reduce the code down to the smallest example that fails. You might also try inspecting the final .osim file (e.g., in the GUI) to ensure the model has been built correctly, and perhaps try passing the .osim filename to the ForwardTool rather than the in-memory Model object.

[Caitlin Romanczyk]

Thomas,

Thank you for your reply! I've fixed the error you caught and am in the process of troubleshooting the rest of the code.

Can I take a step back and ask you a more general question--if I'm trying to model the forward dynamics of a musculoskeletal system with exoskeleton assistance, does it make sense to include the exoskeleton actuators in the same osim model as the musculoskeletal model, or should I be creating a separate file of external forces? I've been reading through the OpenSim documentation, but there are still some things I'm not understanding about how forward dynamics simulations work and how muscles and actuators are treated differently within them.

Really appreciate your help!

Caitlin

[Thomas Uchida]

if I'm trying to model the forward dynamics of a musculoskeletal system with exoskeleton assistance, does it make sense to include the exoskeleton actuators in the same osim model as the musculoskeletal model, or should I be creating a separate file of external forces?

If you know the forces as functions of time, you could apply them to the model as external forces. Alternatively, you could add actuators to the model and design controllers that would compute the exoskeleton forces/torques during the simulation.

[Caitlin Romanczyk]

Thanks again for your reply!

Alternatively, you could add actuators to the model and design controllers that would compute the exoskeleton forces/torques during the simulation.

This is exactly what I'm trying to do--I'm just a little confused on a few points.

1. Do I need to specify states for the actuators, and if so, how? I'm interested in length and lengthening speed (these will be the inputs for my controller), and I know there are functions to obtain these (getLength() and getLengtheningSpeed() ), but I'm not sure how to set them as states for the forward dynamics simulation. As of now, I can see all the states that are generated for the muscles, but none are generated for the actuators.

2. Can you explain the practical differences between using a single controller (with if/else statements) to prescribe control for both the muscles and the actuators, vs. using two controllers (one for muscles and one for actuators) and somehow linking them in the forward dynamics simulation?

3. Is there a benefit to using the ForwardTool, as opposed to setting up an integrator and manager manually as in the TugOfWar example?

I apologize if I'm asking basic questions--I clearly still have some gaps in my understanding, and I really appreciate your help with clarifying! If you can point me toward any examples that use the API to implement exoskeletons with custom actuators and controllers, that would be helpful as well.

Thanks