FD returns "Unknown exception" when using wrapping cylinder

[Sina Porsa]

Hi all.
I am using a model which uses wrapping cylinders for some of the muscles. I am using this model to solve a dynamic optimization problem and I have two problems. First, let's have a look at my algorithem:

- Load the Osim File
- for i=1:1000
--- create a random muscle excitation pattern for each muscle.
--- create a muscle controller based on excitation pattern and add it to the model.
--- perform Forward Dynamics and save the results
--- remove the controller
- end

1st problem: If I use a model which does not have any wrapping cylinders, everything works, but when I use a model with wrapping cylinders, FD returns this error: "Unknown exception".

2nd problem: To increase the speed, I prefer to run the for loop in parallel. Again, I will not have any problems if I use the model without wrapping objects. But when I use the model with wrapping surfaces, it leads to huge increase in memory usage (see the attached JPEG file) which usually ends with this error: "bad allocation"

Thanks.
Sina

P.S: Here is the code that I am using:

Code: Select all

#include <OpenSim/OpenSim.h>
#include <OpenSim/Common/PiecewiseConstantFunction.h>
#include <omp.h>
using namespace OpenSim;
using namespace SimTK;
using namespace std;

const int NO_OF_SIMULATIONS = 32;
const int NO_OF_MUSCLES = 42;
const double FINAL_TIME = 0.5;

Model OsimModel_00 ("Wrapping_Included.osim");
Model OsimModel_01 ("Wrapping_Included.osim");
Model OsimModel_02 ("Wrapping_Included.osim");
Model OsimModel_03 ("Wrapping_Included.osim");
Model OsimModel_04 ("Wrapping_Included.osim");
Model OsimModel_05 ("Wrapping_Included.osim");
Model OsimModel_06 ("Wrapping_Included.osim");
Model OsimModel_07 ("Wrapping_Included.osim");
Model OsimModel_08 ("Wrapping_Included.osim");
Model OsimModel_09 ("Wrapping_Included.osim");
Model OsimModel_10 ("Wrapping_Included.osim");
Model OsimModel_11 ("Wrapping_Included.osim");
Model OsimModel_12 ("Wrapping_Included.osim");
Model OsimModel_13 ("Wrapping_Included.osim");
Model OsimModel_14 ("Wrapping_Included.osim");
Model OsimModel_15 ("Wrapping_Included.osim");

/*Model OsimModel_00 ("Wrapping_Excluded.osim");
Model OsimModel_01 ("Wrapping_Excluded.osim");
Model OsimModel_02 ("Wrapping_Excluded.osim");
Model OsimModel_03 ("Wrapping_Excluded.osim");
Model OsimModel_04 ("Wrapping_Excluded.osim");
Model OsimModel_05 ("Wrapping_Excluded.osim");
Model OsimModel_06 ("Wrapping_Excluded.osim");
Model OsimModel_07 ("Wrapping_Excluded.osim");
Model OsimModel_08 ("Wrapping_Excluded.osim");
Model OsimModel_09 ("Wrapping_Excluded.osim");
Model OsimModel_10 ("Wrapping_Excluded.osim");
Model OsimModel_11 ("Wrapping_Excluded.osim");
Model OsimModel_12 ("Wrapping_Excluded.osim");
Model OsimModel_13 ("Wrapping_Excluded.osim");
Model OsimModel_14 ("Wrapping_Excluded.osim");
Model OsimModel_15 ("Wrapping_Excluded.osim");*/

double calcCost(const double excVec[], Model & OsimModel) {
	double cost;
	try{
		PrescribedController *muscleController= new PrescribedController();
		muscleController->setActuators(OsimModel.updActuators());
		double Time [2] = {0.0, FINAL_TIME};
		double dummyExct [2];
		for (int muscleIndex=0; muscleIndex<NO_OF_MUSCLES; muscleIndex++){
			for (int col=0; col<2; col++){
				dummyExct[col] = excVec[muscleIndex];}
			muscleController->prescribeControlForActuator(muscleIndex, new PiecewiseConstantFunction(2, Time, dummyExct));}

		OsimModel.addController(muscleController);  //Add the controller
		SimTK::State& system = OsimModel.initSystem(); // Initialize the system and get the default state

		// Define non-zero (defaults are 0) states for the free joint
		CoordinateSet& modelCoordinateSet = OsimModel.updCoordinateSet();
		modelCoordinateSet[2].setValue(system,1.0); //pelvis y

		OsimModel.computeEquilibriumForAuxiliaryStates(system);
		SimTK::RungeKuttaMersonIntegrator integrator(OsimModel.getMultibodySystem());  // Create the integrator
		integrator.setAccuracy(1.0e-4);
		Manager manager(OsimModel,  integrator);		// Create the manager
		manager.setInitialTime(0.0);
		manager.setFinalTime(FINAL_TIME);
		manager.setPerformAnalyses(0);
		manager.integrate(system);		// Integrate from initial time to final time

		cost = 10.0;
		OsimModel.removeController(muscleController);
	}
	catch (OpenSim::Exception ex){
		std::cout << ex.getMessage() << std::endl;
        return 1;}
    catch (std::exception ex){
        std::cout << ex.what() << std::endl;
        return 1;}
    catch (...){
        std::cout << "UNRECOGNIZED EXCEPTION" << std::endl;
        return 1;}
	cout << "DONE\n";
	return(cost);
}
double objectiveFunc(const double excMatrix[][NO_OF_MUSCLES],int whichRow) {
	double cost;
	double muscleExcitations [NO_OF_MUSCLES];
	for (int i=0; i<NO_OF_MUSCLES;i++){
		muscleExcitations [i] = excMatrix[whichRow][i];}

	if (omp_get_thread_num()==0)
		cost = calcCost(muscleExcitations, OsimModel_00);
	else if (omp_get_thread_num()==1)
		cost = calcCost(muscleExcitations, OsimModel_01);
	else if (omp_get_thread_num()==2)
		cost = calcCost(muscleExcitations, OsimModel_02);
	else if (omp_get_thread_num()==3)
		cost = calcCost(muscleExcitations, OsimModel_03);
	else if (omp_get_thread_num()==4)
		cost = calcCost(muscleExcitations, OsimModel_04);
	else if (omp_get_thread_num()==5)
		cost = calcCost(muscleExcitations, OsimModel_05);
	else if (omp_get_thread_num()==6)
		cost = calcCost(muscleExcitations, OsimModel_06);
	else if (omp_get_thread_num()==7)
		cost = calcCost(muscleExcitations, OsimModel_07);
	else if (omp_get_thread_num()==8)
		cost = calcCost(muscleExcitations, OsimModel_08);
	else if (omp_get_thread_num()==9)
		cost = calcCost(muscleExcitations, OsimModel_09);
	else if (omp_get_thread_num()==10)
		cost = calcCost(muscleExcitations, OsimModel_10);
	else if (omp_get_thread_num()==11)
		cost = calcCost(muscleExcitations, OsimModel_11);
	else if (omp_get_thread_num()==12)
		cost = calcCost(muscleExcitations, OsimModel_12);
	else if (omp_get_thread_num()==13)
		cost = calcCost(muscleExcitations, OsimModel_13);
	else if (omp_get_thread_num()==14)
		cost = calcCost(muscleExcitations, OsimModel_14);
	else
		cost = calcCost(muscleExcitations, OsimModel_15);

	return cost;
}

int main (int argc, char **argv)	{
	cout << "========\n========\n========\n\n\n";

	unsigned int seed = 0;
	for(int ii=1; ii<argc; ii++) {
		if(strcmp(argv[ii++],"seed") == 0) {
			seed = atoi(argv[ii]);}}
	srand((unsigned)time(0));
	
	double excitation[NO_OF_SIMULATIONS][NO_OF_MUSCLES];
	for (int row=0; row<NO_OF_SIMULATIONS; row++){
		for (int col =0; col<NO_OF_MUSCLES; col++){
			excitation[row][col]= (rand()%100)/100.0;}}

	double COST;

	#pragma omp parallel for num_threads(16)
	for (int i=0; i<NO_OF_SIMULATIONS; i++){
		COST = objectiveFunc(excitation,i);
	}
	
	return 0;
}

[Ayman Habib]

Hi Sina,

The likely problem may have to do with the initialization of muscles hwn they are associated with wrap objects. Please provide a simplified API program to reproduce the problem and we'll help you troubleshoot.

Best regards,
-Ayman

[Sina Porsa]

Hi Ayman and Thank you,
Here is the simple version. If DefaultModel (line 8) contains wapping cylinders, this code will generate an error, but it will work if the loaded model on line 8 does not have wraping objects.
thanks
Sina

P.S: How can I send you the osim files that I am using? I cannot upload them here, as it says the Osim extention is not allowed to be uploaded here.

Code: Select all

#include <OpenSim/OpenSim.h>
using namespace OpenSim;
using namespace SimTK;
using namespace std;
const int NO_OF_SIMULATIONS = 20;
const int NO_OF_MUSCLES = 42;
const double FINAL_TIME = 0.5;
Model DefaultModel ("Wrapping_Included.osim");
void runFD() {
	Model & OsimModel = DefaultModel;
	//Creating random excitation for muscles
	double excitation[NO_OF_MUSCLES][2];
	for (int i =0; i<NO_OF_MUSCLES; i++){
		for (int j=0; j<2; j++){
			excitation[i][j]= (rand()%100)/100.0;}}
	PrescribedController *muscleController= new PrescribedController();
	muscleController->setActuators(OsimModel.updActuators());
	double Time [2] = {0.0, FINAL_TIME};
	for (int muscleIndex=0; muscleIndex<NO_OF_MUSCLES; muscleIndex++){
		muscleController->prescribeControlForActuator(muscleIndex, new PiecewiseLinearFunction(2, Time, excitation[muscleIndex]));}

	OsimModel.addController(muscleController);  //Add the controller
	SimTK::State& system = OsimModel.initSystem(); // Initialize the system and get the default state
	OsimModel.computeEquilibriumForAuxiliaryStates(system);
	SimTK::RungeKuttaMersonIntegrator integrator(OsimModel.getMultibodySystem());  // Create the integrator
	Manager manager(OsimModel,  integrator);		// Create the manager
	manager.setInitialTime(0.0);
	manager.setFinalTime(FINAL_TIME);
	manager.integrate(system);		// Integrate from initial time to final time
}
int main (int argc, char **argv)	{
	unsigned int seed = 0;
	for(int ii=1; ii<argc; ii++) {
		if(strcmp(argv[ii++],"seed") == 0) {
			seed = atoi(argv[ii]);}}
	srand((unsigned)time(0));

	for (int i=0; i<NO_OF_SIMULATIONS; i++){
		runFD();}
	return 0;
}

[Ayman Habib]

Hi Sina,

You can change the extension to .xml and upload it to the forum if it's not too big.

-Ayman

[Sina Porsa]

Hi Ayman,
the Osim files are attached
Tnx