...

C.1: Parse inputs and prepare output

Here's where we load the model. If the setup file specifies that we turn off the coordinate limit forces, we turn them off. We do this in the constructor because it's not something we would need to do in each call to the objective function.

// Parse inputs
// --------------------------------------------------------------------
_name = _tool.get_results_directory();
_cat = OpenSim::Model(_tool.get_model_filename());

// Disable coordinate limit forces?
if (!_tool.get_use_coordinate_limit_forces())
{
    SimTK::State& initState = _cat.initSystem();
    // Loop over all forces in model.
    using OpenSim::CoordinateLimitForce;
    for (int iFor = 0; iFor < _cat.getForceSet().getSize(); iFor++)
    {
        CoordinateLimitForce * LF = dynamic_cast<CoordinateLimitForce *> (&_cat.updForceSet().get(iFor));
        if (LF)
        { // If it is a coordinate limit force, disable it.
            _cat.updForceSet().get(iFor).setDisabled(initState, true);
        }
    }
}

_numOptimSplinePoints = _tool.get_num_optim_spline_points();

We will continue to update a log file, and print out various other files, as the optimization runs. Here's some code that creates a folder for our outputs and opens the log file.

// Prepare output
// --------------------------------------------------------------------

// Create a results directory for all the output files.
#if defined(_WIN32)
    _mkdir(_name.c_str());
#else
    mkdir(_name.c_str(), 0777);
#endif

// Serialize what we just deserialized, in the results directory, to keep
// everything in one organized place.
_tool.print(_name + "/" + _name + "_setup.xml");

// Create a log for the objective function value and the separate terms
// that go into it.
_optLog.open((_name + "/" + _name + "_log.txt").c_str(), std::ofstream::out);
_optLog << "Flippin Felines optimization log: " << _name << std::endl;
time_t rawtime; time(&rawtime);
_optLog << ctime(&rawtime);
_optLog << "Model file name: " << _tool.get_model_filename() << std::endl;

We will continue to update a log file, and print out various other files, as the optimization runs. Here's some code that creates a folder for our outputs and opens the log file.

// Prepare output
// --------------------------------------------------------------------

// Create a results directory for all the output files.
#if defined(_WIN32)
    _mkdir(_name.c_str());
#else
    mkdir(_name.c_str(), 0777);
#endif

// Serialize what we just deserialized, in the results directory, to keep
// everything in one organized place.
_tool.print(_name + "/" + _name + "_setup.xml");

// Create a log for the objective function value and the separate terms
// that go into it.
_optLog.open((_name + "/" + _name + "_log.txt").c_str(), std::ofstream::out);
_optLog << "Flippin Felines optimization log: " << _name << std::endl;
time_t rawtime; time(&rawtime);
_optLog << ctime(&rawtime);
_optLog << "Model file name: " << _tool.get_model_filename() << std::endl;

C.2: Add a controller to the model to control its actuators

...

C.2: Add a controller to the model to control its actuators

Right now, our models have CoordinateActuator's, but we have no way of controlling them. We must add a controller to the model to do this. In general, a controller can manage any number of the model's actuators. We choose to have one controller that controls all of the cat's actuators separately. We want a controller that will specify beforehand what the actuator signal (torque in our case) is throughout the simulation. One alternative would be to use some sort of feedback law. So in our case, we want to use a PrescribedController. With a PrescribedController, we specify for each actuator a function that is the torque that each CoordinateActuator applies as a function of time in the simulation. We choose to use SimmSpline's as the functions.

Now, we want our optimization to solve for the control input (prescribed control) that achieves the flip. So that means our optimization parameters need to be related to our SimmSpline's. The simplest way to relate them is by choosing the optimization parameters to be the Y-values of the spline points. The corresponding X-values (times) are preset as evenly-spaced, and the number of spline points is a setting.

Parameters are ordered by actuator, then by spline point index for a given actuator. Parameters are nondimensionalized by the minimum or maximum control values for the associated actuator.

// Compute the number of optimization parameters for the model.
_numActuators = _cat.getActuators().getSize();
int numParameters = _numActuators * _numOptimSplinePoints;
setNumParameters(numParameters);

// Size the vector of best-yet splines.
_splinesBestYet.resize(_numActuators);

// Add a PrescribedController to the model.
using OpenSim::PrescribedController;
PrescribedController * flipController = new PrescribedController();
flipController->setName("flip_controller");

// Set the PrescribedController to control all actuators.
flipController->setActuators(_cat.getActuators());

// Add the PrescribedController to the model.
_cat.addController(flipController);

// Create SimmSpline's for each actuator and push into vector.
double indexToTimeInSeconds =
    _duration / (double)(_numOptimSplinePoints - 1);
for (int i = 0; i < _numActuators; i++)
{
    // Create a spline function for the current actuator.
    _splines.push_back(new OpenSim::SimmSpline());
    // Name this spline with the name of the current actuator.
    _splines[i]->setName(_cat.getActuators().get(i).getName());
    // Add the correct number of points to this spline.
    for (int j = 0; j < _numOptimSplinePoints; j++)
    {
        // Scale the time points appropriately.
        double thisTime = (double)j * indexToTimeInSeconds;
        _splines[i]->addPoint(thisTime, 0.0);
    }
    // Tell the controller about this function.
    flipController->prescribeControlForActuator(i, _splines[i]);
} 

...

 about this function.
    flipController->prescribeControlForActuator(i, _splines[i]);
} 

C.3: Set parameter limits/bounds for the optimization

We tell the optimizer that there are bounds on the values of the parameters. We've chosen to nondimensionalize our parameters, so our bounds are simple. The OptimizerSystem::setParameterLimits() function is public, so sometimes you might see that this function is be called in optimize.cpp instead.

// Set (nondimensional) parameter limits/bounds for the actuators.
SimTK::Vector lowerLimits(getNumParameters(), -1.0);
SimTK::Vector upperLimits(getNumParameters(), 1.0);
setParameterLimits(lowerLimits, upperLimits);

...

C.5: Unpack parameters into the model

The Optimizer has given us new parameters to try  out. We need to use these parameters to change the Y-values of the spline. First, we must dimensionalize the parameters.

Note that the number of calls to the objective function (_objectiveCalls) is NOT equal to the number of iterations of the optimization algorithm.

// Increment the count of calls to this function.
_objectiveCalls++;

// Unpack parameters into the model (i.e., update spline points)
// --------------------------------------------------------------------
for (int iAct = 0; iAct < _numActuators; iAct++)
{
    // Max/min for all spline points for the i-th actuator.
    double minControl = _cat.getActuators().get(iAct).getMinControl();
    double maxControl = _cat.getActuators().get(iAct).getMaxControl();

    // Visit each point in this spline.
    for (int iPts = 0; iPts < _numOptimSplinePoints; iPts++)
    {
        int paramIndex = iAct * _numOptimSplinePoints + iPts;

        // Dimensionalize the control value.
        double dimControlValue;
        if  (parameters[paramIndex] < 0 && minControl < 0)
        { // parameter and minControl are both neg; must negate again.
            dimControlValue = -minControl * parameters[paramIndex];
        }
        else
        {
            dimControlValue = maxControl * parameters[paramIndex];
        }

        // Edit the spline with the dimensional control value.
        _splines[iAct]->setY(iPts, dimControlValue);
    }
}

C.6: Run a forward dynamic simulation6: Run a forward dynamic simulation

Now that we've updated the control functions, we can run a simulation and see how the cat ends up.

SimTK::State& initState = _cat.initSystem();

// Construct an integrator.
SimTK::RungeKuttaMersonIntegrator integrator(_cat.getMultibodySystem());
integrator.setAccuracy(1.0e-6);

// Construct a manager to run the integration.
OpenSim::Manager manager(_cat, integrator);
_cat.getMultibodySystem().realize(initState, SimTK::Stage::Acceleration);

// Integrate from initial time to final time
manager.setInitialTime(0.0);
manager.setFinalTime(_duration);
manager.integrate(initState);

C.7: The objective function value

Here's the brain of the optimization. We assign to f the value of objective function. Here, we have 3 terms in our objective function, which are described on the main Flippin' Felines page. The source code contains a few more terms; we've cut them out for simplicity.

// Will be writing to a log while constructing the objective function.
if (_objectiveCalls % _outputPeriod == 0)
    _optLog << _objectiveCalls;

// Copy the model's state. We need a consistent state to access the
// model's configuration and related speeds.
SimTK::State aState = initState;
_cat.getMultibodySystem().realize(aState, SimTK::Stage::Acceleration);
const OpenSim::CoordinateSet& coordinates = _cat.getCoordinateSet();
double roll = coordinates.get("roll").getValue(aState);
double rollRate = coordinates.get("roll").getSpeedValue(aState);
double rollAccel = coordinates.get("roll").getAccelerationValue(aState);
double twist = coordinates.get("twist").getValue(aState);
double twistRate = coordinates.get("twist").getSpeedValue(aState);
double twistAccel = coordinates.get("twist").getAccelerationValue(aState);
double hunch = coordinates.get("hunch").getValue(aState);
double pitch = coordinates.get("pitch").getValue(aState);
double pitchRate = coordinates.get("pitch").getSpeedValue(aState);
double pitchAccel = coordinates.get("pitch").getAccelerationValue(aState);

// --------------------------------------------------------------------
f = 0;
f += anteriorLegsDownTerm(roll, rollRate, rollAccel);
f += posteriorLegsDownTerm(twist, twistRate, twistAccel);
f += sagittalSymmetryTerm(hunch, pitch, pitchRate, pitchAccel);
// --------------------------------------------------------------------

We use some helper functions to keep this code clean. This first helper function is called by each of the "*Term()" functions above, and provides a convenient way of, in one step, updating the log and weighting a term in the function.

/**
 * Helper function for writing data to the log.
 * */
double processObjTerm(std::string name, double weight, double term) const
{
    if (_objectiveCalls % _outputPeriod == 0)
        _optLog << " " << name << " " << weight * term;
    return weight * term;
}

...

There are some functions that we use, but that aren't that too important for learning how optimization works. They are used to print either the model or control functions.

...

// Reference to the FlippinFelinesOptimizerTool.
OpenSim::FlippinFelinesOptimizerTool & _tool;

// Simulation runtime (seconds).
double _duration;

// Model containing the attributes described in this class' description
mutable OpenSim::Model _cat;

// Number of actuators in the model.
int _numActuators;

// Number of points used to define splines for actuator control.
int _numOptimSplinePoints;

// Vector of splines used in the PrescribedController.
std::vector<OpenSim::SimmSpline *> _splines;

// Vector of splines that gave the best objective value yet.
mutable std::vector<OpenSim::SimmSpline> _splinesBestYet;

// Number of calls to the objective function.
mutable int _objectiveCalls;

// Name of results directory.
std::string _name;

// Log for recording data during optimization
mutable std::ofstream _optLog;

// Period for how often objective terms are printed to the log.
static const int _outputPeriod = 100;

// Best (lowest) value of the objective function, for logging.
mutable double _objectiveFcnValueBestYet;

// To aid with conservative printing of best yet actuation:
mutable bool _lastCallWasBestYet;
mutable bool _thisCallIsBestYet; 

TODO MUTABLES

TODO talk about controls

TODO talk about how the parameters work.

C.2 is the part that would be custom between different pplz

CAN BE ANYTHING

iterations is not number of objective function calls.

There's an additional file, optimizer_tool_template.cpp, whose code we do not show here. It generates a template XML file TODO .... TODO paste a version of it into this page.

Getting started with optimization and reproducing our results

...