Creating an Optimization

The steps to creating an optimization are:

Overview

In this section, we will write a main program to perform an optimization study using OpenSim. We will build it up in pieces, starting by programmatically loading an existing OpenSim model. The model will be a simple arm model named Arm26.osim, consisting of 2 degrees of freedom and 6 muscles. We will then define an optimization problem that finds a set of muscle controls to maximize the forward velocity of the forearm/hand segment mass center. The resulting source code and associated files for this example come with the OpenSim 2.0 distribution under the directory:

C:\Program Files\OpenSim 2.0\sdk\APIExamples\OptimizationExample_Arm26

As in Performing a Simulation, the following sections explain the steps to create your own main program. Additionally, we will be extending existing optimizer classes in the OpenSim API. For more information on the OptimizerSystem class, see the SimTKmath User's Guide available on the SimTK project site (https://simtk.org/home/simtkcore under the "Documents" tab).

worddavdb3909b4fe9323a260e6d6b21c562f23.png [imported from a Word document]

Extending the OptimizerSystem class

Before we get into extending the class, we need to include the proper header files and define a few global variables.

#include <OpenSim/OpenSim.h>
using namespace OpenSim;
using namespace SimTK;
int stepCount = 0; 
 
// Global variables to define integration time window, optimizer step 
// count, the best solution.
double initialTime = 0.0;
double finalTime = 0.25;
double bestSoFar = Infinity;

In OpenSim, optimization problems are set up within an OptimizerSystem, which uses the SimTK-level algorithms to determine a solution. To set up our optimization problem, we need to create our own OptimizerSystem, called ExampleOptimizationSystem, by extending the existing base OptimizerSystem class.

class ExampleOptimizationSystem : public OptimizerSystem {
  public: 
    /* Constructor class. Parameters accessed in objectiveFunc() class */
    ExampleOptimizationSystem(int numParameters, State& s, Model& aModel): 
    numControls(numParameters), OptimizerSystem(numParameters), si(s), osimModel(aModel){} 

  /* The objectiveFunc() class will go here. */
  private:
    int numControls;
    State& si;
    Model& osimModel;
};

Writing the main()

We can perform an optimization by creating our own main program that will invoke our OptimizerSystem.

int main()
{
  try {
    // Create a new OpenSim model
    Model osimModel("Arm26_Optimize.osim"); 

    /* The guts of your main() will go here */
  }

  catch (std::exception ex)
  {
    std::cout << ex.what() << std::endl;
    return 1;
  } 

  // End of main() routine.
  return 0;
}

Defining controls and initializing muscle states

As in Chapter 2, we need to define the controls (i.e., muscle excitations), initial activation, and fiber length of each muscle. Moreover, we initialize the states for each muscle after setting the states.

// Define the initial and final controls
ControlLinear *control_TRIlong = new ControlLinear();
ControlLinear *control_TRIlat = new ControlLinear();
ControlLinear *control_TRImed = new ControlLinear();
ControlLinear *control_BIClong = new ControlLinear();
ControlLinear *control_BICshort = new ControlLinear();
ControlLinear *control_BRA = new ControlLinear(); 
 
/* NOTE: Each contoller must be set to the corresponding 
* muscle name in the model file. */
control_TRIlong->setName("TRIlong"); control_TRIlat->setName("TRIlat"); 
control_TRImed->setName("TRImed"); control_BIClong->setName("BIClong");
control_BICshort->setName("BICshort"); control_BRA->setName("BRA"); 
 
ControlSet *muscleControls = new ControlSet();
muscleControls->append(control_TRIlong);
muscleControls->append(control_TRIlat);
muscleControls->append(control_TRImed);
muscleControls->append(control_BIClong);
muscleControls->append(control_BICshort);
muscleControls->append(control_BRA); 
 
ControlSetController *muscleController = new ControlSetController();
muscleController->setControlSet(muscleControls);
muscleController->setName("MuscleController"); 
 
// Add the controller to the model
osimModel.addController(muscleController); 
 
// Initialize the system and get the state 
State& si = osimModel.initSystem(); 
 
// Define the initial muscle states
const OpenSim::Set<OpenSim::Actuator> &muscleSet = 
osimModel.getActuators();
 
for(int i=0; i< muscleSet.getSize(); i++ ){
  Actuator* act = &muscleSet.get;
  OpenSim::Muscle* mus = dynamic_cast<Muscle*>(act);
  mus->setDefaultActivation(0.5);
  mus->setDefaultFiberLength(0.1);
} 
 
// Make sure the muscles states are in equilibrium
osimModel.equilibrateMuscles(si);

Define the optimizer

In SimTK and OpenSim, an Optimizer operates on an OptimizationSystem, which we will initialize as an ExampleOptimizerSystem . We then define the bounds for the parameters of the problem, the optimizer tolerance, and the numerical gradient flag before finally invoking the optimizer.

// Number of controls will equal the number of muscles in the model
int numControls = 6; 
 
// Initialize the optimizer system we've defined.
ExampleOptimizationSystem sys(numControls, si, osimModel);
Real f = NaN; 
 
/* Define and set bounds for the parameter we will optimize */
Vector lower_bounds(numControls);
Vector upper_bounds(numControls); 

for(int i=0;i<numControls;i++) {
  lower_bounds[i] = 0.01;
  upper_bounds[i] = 1.0;
}
 
sys.setParameterLimits( lower_bounds, upper_bounds ); 
 
// Set the initial values (guesses) for the controls
Vector controls(numControls);
controls[0] = 0.01; controls[1] = 0.01;
controls[2] = 0.01; controls[3] = 0.01;
controls[4] = 0.01; controls[5] = 0.01; 

try {
  // Create an optimizer. Pass in our OptimizerSystem
  // and the name of the optimization algorithm.
  Optimizer opt(sys, SimTK::LBFGSB); 

  // Specify settings for the optimizer
  opt.setConvergenceTolerance(0.05);
  opt.useNumericalGradient(true);
  opt.useNumericalJacobian(true); 

  // Optimize it!
  f = opt.optimize(controls);
}

catch(const std::exception& e) {
  std::cout << "Caught exception :" << std::endl;
  std::cout << e.what() << std::endl;
}

Writing the objective function

Within ExampleOptimizationSystem, we need to define our objective function as a public member of the class. This member function will take the parameters of the muscle controls that we want to vary and will return a real number about the performance we want to optimize (i.e., forward velocity of the hand) , which will then be minimized (Note: To maximize a value, just multiply it by -1). In this case, the parameters we want to vary are the muscle control values, and we will return a real number determined by our objective function (i.e., the forearm/hand mass center velocity). Additionally, if we had an analytical gradient or Jacobian function for our system, they could also be defined as member functions of the ExampleOptimizationSystem.

int objectiveFunc(const Vector &newControls, const bool new_coefficients, Real& f ) const { 

 // Make a copy of out initial states
:State s = si; 

// Access the controller set of the model and update the control values
((ControlSetController *)(&osimModel.updControllerSet()[0]))>updControlSet()>setControlValues(initialTime, &newControls[0]);
((ControlSetController *)(&osimModel.updControllerSet()[0]))>updControlSet()>setControlValues(finalTime, &newControls[0]); 

// Create the integrator for the simulation.
RungeKuttaMersonIntegrator integrator(osimModel.getMultibodySystem());
integrator.setAccuracy(1.0e-4); 

// Create a manager to run the simulation
Manager manager(osimModel, osimModel.getMultibodySystem(), integrator); 

// Integrate from initial time to final time
manager.setInitialTime(initialTime);
manager.setFinalTime(finalTime);
osimModel.getMultibodySystem().realize(s, Stage::Acceleration);
manager.integrate(s); 
 
/* Calculate the scalar quantity for the optimizer to minimize 
* In this case, we're maximizing forward velocity of the 
* forearm/hand mass center so compute the velocity and
just multiply it by -1./
Vec3 massCenter;
Vec3 velocity;
osimModel.getBodySet().get("r_ulna_radius_hand").getMassCenter(massCenter);
osimModel.getMultibodySystem().realize(s, Stage::Velocity);
osimModel.getSimbodyEngine().getVelocity(s,osimModel.getBodySet().get("r_ulna_radius_hand"),
massCenter, velocity); 
f = -velocity[0];
stepCount++; 
 
/* Use an if statement to only store and print the results of an 
* optimization step if it is better than a previous result.
*/
if( f < bestSoFar){
  Storage statesDegrees(manager.getStateStorage());
  osimModel.updSimbodyEngine().convertRadiansToDegrees(statesDegrees);
  statesDegrees.print("bestSoFar_states_degrees.sto");
  bestSoFar = f;
  std::cout << "\nOptimization Step #: " << stepCount << " controls = " << newControls << " bestSoFar = " << f << std::endl;
} 
 
return(0);
}

Now you can build and run your main program, and then load the model and results into OpenSim to visualize the optimized control pattern and resulting kinematics.

worddav9399476b35521358079e8139e597acb6.png [imported from a Word document]

worddav03be4cb724694f6296e41cdd78709fb1.png [imported from a Word document]

Next: Creating a Customized Actuator

Home: Adding New Functionality | Developer's Guide