To get our models to flip"right themselves", we want to add actuation control to them, control their actuation and tweak the this control until we obtain our desired flip a flipping motion. We have an optimization algorithm take care of the tweaking for us. The algorithm will tweak During each step of the optimization, the algorithm tweaks the control in some intelligent way, then we 'll run a forward dynamic simulation with the given tweaked controls, and we 'll see finally check if the cat flipshas flipped. The algorithm tweaking continues to tweak the control until we decide that the cat model has achieved a flip is achieved.

This procedure of algorithmically tweaking control inputs, running a forward dynamic simulation, then judging how well the model/simulation achieved some objective is called dynamic optimization. It is the way that biomechanists would prefer to study human movement (see Sky Higher and Miller, 2012 for just two of many exampleexamples) , and is commonly viewed as the alternative to tracking simulations. The models used for studying such motions are complex, and so it takes a long time to run a dynamic optimization (days to weeks). For cat flipping, however, an optimization takes only an hour or two. Thus, it's an ideal problem for getting one's feet wet with first learning dynamic optimization.

Code Organization

...

First, we present an overview of the 2 source files we use for the optimization. The file optimize.cpp creates an executable that (mostly just) creates the necessary objects , and runs the optimization. The code block below shows only 4 displays the four most important lines from this file, the most important lines, that which show how we connect the various pieces of the optimization. Look at the next code block, a skeleton of the second file, FlippinFelinesOptimizerSystem.h, to see . It is where we define the FlippinFelinesOptimizerTool and FlippinFelinesOptimizerSystem classes.

...

To run any optimization, it's necessary to use SimTK::Optimizer and a class like FlippinFelinesOptimizerSystem. However, a class like FlippinFelinesOptimizerTool is not necessary; it's just something we've done to make it easier to run these optimizations. (TODO so we'll focus on the first 2).

Code Block
language cpp title FlippinFelinesOptimizerSystem.h (incomplete snippet)
namespace OpenSim { // Manages inputs to the optimization via OpenSim's XML abilities. class FlippinFelinesOptimizerTool : public Object { // ********** PART B ********** }; } // end namespace OpenSim // Finds a control input that achieves flip, as determined by objective function. class FlippinFelinesOptimizerSystem : public SimTK::OptimizerSystem { // ... FlippinFelinesOptimizerSystem(OpenSim::FlippinFelinesOptimizerTool & tool) : _tool(tool), // ... { // ... } // ... // The meat of this class; this is called by the Optimizer that has this system. int objectiveFunc(const SimTK::Vector & parameters, bool new_parameters, SimTK::Real & f) const { // ... f = ...; // ... } // ... // Reference to the FlippinFelinesOptimizerTool. OpenSim::FlippinFelinesOptimizerTool & _tool; // ... };

As you can see, our optimizer system must be a subclass of SimTK::OptimizerSystem. We override its objectiveFunc() function, which is then called by the SimTK::Optimizer as part of its algorithm. It's in objectiveFunc() that we (1) convert the optimizer's tweaks to optimization parameters into changes in the actuation control of the model, (2) run a forward dynamic simulation of the cat's fall, and (3) evaluate the simulation by assigning a value for our objective function, f, which is what we calledon the main Flippin' Felines page. See SimTK::OptimizerSystem reference and SimTK::Optimizer reference for more information.

A:

...

optimize.cpp:

...

the executable that ties the pieces together

Here, we show the complete optimize.cpp. The user must provide one command-line argument, which is the name of an XML file that defines an a FlippinFelinesOptimizerTool object (we'll get to this in the next section). A serialization is a representation of an object or objects in data files (which, in our case, is a plain text file in XML format). See this page for more information about C++ command-line arguments.

#include <iostream>
#include <OpenSim/OpenSim.h>
#include "FlippinFelinesOptimizerSystem.h"

int main(int argc, char * argv[])
{
    // argc is the number of command line inputs, INCLUDING the name of the
    //      exectuable. Thus, it'll always be greater than/equal to 1.
    // argv is an array of space-delimited command line inputs, the first one
    //      necessarily being the name of the executable (e.g., "optimize
    //      optimize_input_template.xml").

    // Get the filename of the FlippinFelinesOptimizerTool serialization.
    if (argc == 2)
    { // Correct number of inputs.

        // Set-up
        // --------------------------------------------------------------------

        // Parse inputs using our Tool class.
        std::string toolSetupFile = argv[1];
        OpenSim::FlippinFelinesOptimizerTool tool(toolSetupFile);

...

We create an optimizer that will optimize our FlippinFelinesOptimizerSystem, and tell the optimizer to figure out the best algorithm to use for the optimization. We played around with the settings a littlebit. We set the maximum number of iterations to be large just because we don't want so that the optimizer to does not quit early.

        // Create the optimizer with our system and the "Best Available"
        // algorithm.
        SimTK::Optimizer opt(sys, SimTK::BestAvailable);

        // Set optimizer settings.
        opt.setConvergenceTolerance(0.001);
        opt.useNumericalGradient(true);
        opt.setMaxIterations(100000);
        opt.setLimitedMemoryHistory(500);

When we call optimize(), we need to give the optimizer a place to start in the parameter space of the parameters. Typically, one would define these initial parameters right here, in the same function where we call optimize(). However, we did not want to hard-code our initial parameters in the source code; we wanted to change this them between optimizations easily. So it's actually something we can specify , the initial parameters are actually specified in the toolSetupFile. The FlippinFelinesOptimizerSystem parses that input , and gives us the initial parameters in a Vector if we call FlippinFelinesOptimizerSystem::initialParameters(). We then pass this to the Optimizer.

The function optimize() returns when the optimization finishes. Once that happens, we print out the model (that should now be able to flip), as well as the actuation controls that produced the flip. In the case that anything goes wrong (an exception is thrown), we still want to see what the optimizer ended up finding (so as not to waste the failed optimization).

        // Initialize parameters for the optimization as those determined
        // by our OptimizerSystem.
        SimTK::Vector initParameters = sys.initialParameters();

		// And we're off! Running the optimization
        // --------------------------------------------------------------------
        try
        {
            double f = opt.optimize(initParameters);
            // Print the optimized model so we can explore the resulting motion.
            sys.printModel(name + "_optimized.osim");
            // Print the control splines so we can explore the resulting actuation.
            sys.printPrescribedControllerFunctionSet(
                    name + "_optimized_parameters.xml");
            
            // Print out the final value of the objective function.
            std::cout << "Done with " << name << "! f = " << f << std::endl;
        }
        catch (...)
        {
            // Print the last model/controls (not optimized) so we have something
            // to look at.
            sys.printModel(name + "_last.osim");
            sys.printPrescribedControllerFunctionSet(
                    name + "_last_parameters.xml");
            std::cout << "Exception thrown; optimization not achieved." << std::endl;
            // Don't want to give the appearance of normal operation.
            throw;
        }
    }
    else
    { // Too few/many inputs, etc.
        std::cout << "\nIncorrect input provided. "
            "Must specify the name of a FlippinFelinesOptimizerTool "
            "serialization (setup/input file).\n\nExamples:\n\t"
            "optimize optimize_input_template.xml\n" << std::endl;
        return 0;
    }
    return EXIT_SUCCESS;
};

...

B: FlippinFelinesOptimizerTool: optimization settings and serialization

Although this part of the Below is a skeleton of the FlippinFelinesOptimizerTool definition. There are three main parts to the definition; we'll address each individually. Although this code is not vital for setting up optimizations, we discuss it now , since the rest of the our code depends on it. Below is a skeleton of the FlippinFelinesOptimizerTool definition. There are three main parts to the definition, which we'll get to shortly.This class

The FlippinFelinesOptimizerTool class operates similarly to OpenSim's tools (e.g., for the  Forward Tool), which can be run via the command line by providing an XML setup file. The way weWe'd like to use our optimizer is via the command like via commands like in a similar way: for example, optimize optimize_setup_file.xml. This Of course, this requires the ability to parse that the provided XML file. That's what this FlippinFelinesOptimizerTool the class is all about.

If we can define a class whose member variables are the settings for our optimization, and we can serialize this class into a text file, then we have a record of the settings we used for a given optimization. Furthermore, if we can deserialize a text file to create an instance of the class we originally serialized, then we have a really nice way to read that the optimize_setup_file.xml input file to set the settings within for our optimization.

Fortunately, all OpenSim classes have the ability to serialize themselves, themselves and to deserialize a file into an instance of an OpenSim class. To create a class that has these features, all we need to do is subclass from OpenSim::Object, and  and define our member variables with a specific syntax. OpenSim uses its ability to de/serialize all the time (all XML files produced by OpenSim or given to OpenSim are serializations done in this way), but we can use this feature for our own purposes as well. That's what we do here; we call our subclass FlippinFelinesOptimizerTool, since  since the class functions similarly (but is different!) from OpenSim Tool and AbstractTool classes.

...

namespace OpenSim {

class FlippinFelinesOptimizerTool : public Object {

OpenSim_DECLARE_CONCRETE_OBJECT(FlippinFelinesOptimizerTool, Object);

public:

    // Declare properties of the class.
    // ------------------------------------------------------------------------
    // ********** PART B.1 **********

    // Constructors
    // ------------------------------------------------------------------------
    // ********** PART B.2 **********

private:

    void setNull() { }

    void constructProperties()
    {
	    // ********** PART B.3 **********
    }
};

B.1: Declare properties of the class

We create member variables using the OpenSim_DECLARE_PROPERTY and OpenSim_DECLARE_OPTIONAL_PROPERTY macros, which are described here. For simplicity, we've omitted many of the properties we have in the actual source code. For a member variable to be de/serializable, it must be delcared via this macro. Some of the properties help with general setup (e.g., results_directory, model_filename), while others modify the computation of the objective function value. These (e.g., desired configurations and weightings). The macros define getters/setters for us; we'll see the usage of the getters in other parts of the code.

// General properties.
OpenSim_DECLARE_PROPERTY(results_directory, std::string,
        "Directory in which to save optimization log and results");
OpenSim_DECLARE_PROPERTY(model_filename, std::string,
        "Specifies path to model file, WITH .osim extension");
OpenSim_DECLARE_PROPERTY(num_optim_spline_points, int,
        "Number of points being optimized in each spline function. "
        "Constant across all splines. If an initial_parameters_filename is "
        "provided, the functions specified in that file must have the "
        "correct number of points. We do not error-check for this.");
// Properties related to the objective function.
OpenSim_DECLARE_PROPERTY(anterior_legs_down_weight, double,
        "Adds terms to the objective to minimize final value of "
        "(roll - Pi) and related speeds");
OpenSim_DECLARE_PROPERTY(posterior_legs_down_weight, double,
        "Adds terms to the objective to minimize final value of "
        "(twist - 0) and related speeds");
OpenSim_DECLARE_PROPERTY(sagittal_symmetry_weight, double,
        "Adds a term to the objective to minimize final value of "
        "(hunch + 2 * pitch)");
// Modifying the model before optimizing.
OpenSim_DECLARE_PROPERTY(use_coordinate_limit_forces, bool,
        "TRUE: use coordinate limit forces, "
        "FALSE: ignore coordinate limit forces");
// Setting initial parameters for the optimization.
OpenSim_DECLARE_OPTIONAL_PROPERTY(initial_parameters_filename, std::string,
        "File containing FunctionSet of SimmSpline's used to initialize "
        "optimization parameters. If not provided, initial parameters are "
        "all 0.0, and this element must be DELETED from the XML file "
        "(cannot just leave it blank). The name of each function must be "
        "identical to that of the actuator it is for. x values are ignored. "
        "The time values that are actually used in the simulation are "
        "equally spaced from t = 0 to t = 1.0 s, and there should be "
        "as many points in each function as given by the num_optim_spline_points "
        "property. y values should be nondimensional and between -1 and 1 "
        "(negative values normalized by minControl if minControl is "
        "negative; otherwise the value is normalized by maxControl). "
        "NOTE that the output optimized splines are NOT NONDIMENSIONAL. "
        "Be careful, we do not do any error checking.");

B.2: Constructors

Classes are serialized by calling their print() method. To be able to deserialize a file, the associated class must have a constructor that takes in a filename , and passes it to the constructor of Object. The properties in the class are overwritten by the values in the serialization by calling the method Object::updateFromXMLDocument().

...

Before we can do that though, we need to construct the properties we've declared.

B.3: Constructing the properties we've declared

These constructProperty functions, generated by the macros in part B.1, require an initial/default value for the property.

void constructProperties()
{
    constructProperty_results_directory("results");
    constructProperty_model_filename("flippinfelines_*FILL THIS IN*.osim");
    constructProperty_num_optim_spline_points(20);
    constructProperty_anterior_legs_down_weight(1.0);
    constructProperty_posterior_legs_down_weight(1.0);
    constructProperty_sagittal_symmetry_weight(1.0);
    constructProperty_use_coordinate_limit_forces(true);
    constructProperty_initial_parameters_filename("");
}

...

Sample Input

Below is what the XML serialization of this class looks like. This is exactly what we meant by the optimizer_setup_file.xml above. This is a simplified version of what's generated by the executable optimize_input_template, which is generated from optimize_input_template.cpp. We've given you this file, but we do not discuss it in this the tutorial.

<?xml version="1.0" encoding="UTF-8" ?>
<OpenSimDocument Version="30000">
    <FlippinFelinesOptimizerTool>
        <!--Directory in which to save optimization log and results-->
        <results_directory>results</results_directory>
        <!--Specifies path to model file, WITH .osim extension-->
        <model_filename>flippinfelines_*FILL THIS IN*.osim</model_filename>
        <!--Number of points being optimized in each spline function. Constant across all splines. If an initial_parameters_filename is provided, the functions specified in that file must have the correct number of points. We do not error-check for this.-->
        <num_optim_spline_points>20</num_optim_spline_points>
        <!--Adds terms to the objective to minimize final value of (roll - Pi) and related speeds-->
        <anterior_legs_down_weight>1</anterior_legs_down_weight>
        <!--Adds terms to the objective to minimize final value of (twist - 0) and related speeds-->
        <posterior_legs_down_weight>1</posterior_legs_down_weight>
        <!--Adds a term to the objective to minimize final value of (hunch + 2 * pitch)-->
        <sagittal_symmetry_weight>1</sagittal_symmetry_weight>
        <!--TRUE: use coordinate limit forces, FALSE: ignore coordinate limit forces-->
        <use_coordinate_limit_forces>true</use_coordinate_limit_forces>
        <!--File containing FunctionSet of SimmSpline's used to initialize optimization parameters. If not provided, initial parameters are all 0.0, and this element must be DELETED from the XML file (cannot just leave it blank). The name of each function must be identical to that of the actuator it is for. x values are ignored. The time values that are actually used in the simulation are equally spaced from t = 0 to t = 1.0 s, and there should be as many points in each function as given by the num_optim_spline_points property. y values should be nondimensional and between -1 and 1 (negative values normalized by minControl if minControl is negative; otherwise the value is normalized by maxControl). NOTE that the output optimized splines are NOT NONDIMENSIONAL. Be careful, we do not do any error checking.-->
        <initial_parameters_filename>initial_parameters.xml</initial_parameters_filename>
    </FlippinFelinesOptimizerTool>
</OpenSimDocument> 

...

C: FlippinFelinesOptimizerSystem

Here's a skeleton of the FlippinFelinesOptimizerSystem class. The two functions that are the most important are the constructor and objectiveFunc(). We focus on these. The definition of an optimization problem consists primarily of (1) parameters/inputs, and (2) the objective function. The parts that deal with the parameters are C.2, C.3, C.4, and C.5. Parts C.6 and C.7 deal with the objective function.

 /**
 * Finds a control input history that achieves certain desired features
 * of a cat's flipping maneuver, as determined by the objective function.
 * The control of the system is performed via a PrescribedController that
 * uses a spline for all actuators.
 * */
class FlippinFelinesOptimizerSystem : public SimTK::OptimizerSystem {
public:

    FlippinFelinesOptimizerSystem(OpenSim::FlippinFelinesOptimizerTool & tool) :
        _tool(tool),
        _duration(1.0),
        _objectiveCalls(0),
        _objectiveFcnValueBestYet(SimTK::Infinity),
    {
        // Parse inputs & prepare output
        // --------------------------------------------------------------------
	    // ********** PART C.1 **********

        // Add a controller to the model to control its actuators
        // --------------------------------------------------------------------
        // ********** PART C.2 **********

		// Set parameter limits/bounds for the optimization.
        // --------------------------------------------------------------------
        // ********** PART C.3 **********

    }

    SimTK::Vector initialParameters()
    {
        // ********** PART C.4 **********
    }
	// ...

    /**
     * Defines desired features for the cat model's flip. Each call to this
     * function updates the model's control, runs a forward dynamic simulation
     * with this control, computes the objective value resulting from the
     * simulation, and saves the results to output.
     * */
    int objectiveFunc(const SimTK::Vector & parameters, bool new_parameters,
            SimTK::Real & f) const
    {
        // ...
        // Unpack parameters into the model (i.e., update spline points)
        // --------------------------------------------------------------------
        // ********** PART C.5 **********

		// Run a forward dynamic simulation
        // --------------------------------------------------------------------
        // ********** PART C.6 **********

        // Construct the objective function value, term by term
        // --------------------------------------------------------------------
        // ********** PART C.7 **********

        // Update the log and outputs with the objective function value
        // --------------------------------------------------------------------
        // ********** PART C.8 **********

    }
    // Miscellaneous functions (used, but uninteresting)
    // ------------------------------------------------------------------------
    // ********** PART C.9 **********

private:

    // Objective function terms
    // ------------------------------------------------------------------------
    // ********** PART C.7 **********

    	// Member variables
    // ------------------------------------------------------------------------
    	// ********** PART C.10 **********

};

Throughout, we'll be making use of member variables we've defined. To see their declarations, see part Cjump to part C.10.

Note that objectiveFunc() is necessarily a const function. That's This is because we wouldndon't want the state of the system to change changing as we are optimizing optimize it. For example, we wouldn't want to change relative masses of the cat while trying to find the optimum control input; the solution would be a moving target.

C.1: Parse inputs

...

& 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 during 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();

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

We will continue to update a log file, and print out various other files (e.g., best-yet controls), as the optimization runs. Here's the 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. More specifically, we want a controller that will specify beforehand what the actuator signal (torque, in our case) will be throughout the simulation (NOTE: one alternative would be to use a feedback law). Such a controller is a PrescribedController. With a PrescribedController, we specify for each actuator a function that is the torque that each CoordinateActuator applies as a function of simulation time. We choose to use splines (SimmSpline) as the functions.

We want our optimization to solve for the control input (prescribed control) that flips the cat. So, 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. As seen above, the number of spline points is specified in the setup file; the corresponding X-values (times) are preset as evenly-spaced over the simulation.

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]);
} 

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 this function might 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.4: Initial parameters

We wanted to make it really easy for the user to specify the initial parameters for the optimization. Since they these parameters end up as Y-coordinates of splinesthe spline points, we thought this would that SimmSpline's would be the easiest way for a user to specify the initial parametersform of input. However, we need the initial parameters as a Vector. So, we wrote a method to convert the Y-coordinates from a FunctionSet of SimmSpline's into a Vector.

SimTK::Vector initialParameters()
{
    SimTK::Vector initParams(getNumParameters(), 0.0);
    if (_tool.get_initial_parameters_filename() != "")
    { // A file is specified.
        // Deserialize the initial parameters XML file.
        OpenSim::FunctionSet initFcns(_tool.get_initial_parameters_filename());
        // Write a copy of the FunctionSet to the results directory.
        initFcns.print(_name + "/" + _name + "_initial_parameters.xml");
        OpenSim::Array<std::string> initNames;
        initFcns.getNames(initNames);
        
        // This loop is set up so that the initFcns do not need to be in the same
        // order as the optimization parameters. Also, the number of initFcns
        // specified by the initial parameters file can be greater than the
        // number of actuators in the model (assuming that the initial parameters
        // file AT LEAST contains the model's actuators).
        using OpenSim::SimmSpline;
        for (int iAct = 0; iAct < _numActuators; iAct++)
        {
            // Get index of the initFcn corresponding to the actuator.
            int iFcn = initFcns.getIndex(_cat.getActuators().get(iAct).getName());
            // Get access to the spline's methods.
            SimmSpline * fcn =
                dynamic_cast<SimmSpline *>(&initFcns.get(iFcn));
            for (int iPts = 0; iPts < _numOptimSplinePoints; iPts++)
            {
                // Find the right index in the optimization parameters.
                int paramIndex = iAct * _numOptimSplinePoints + iPts;
                // Transfer y value from input to initial parameters.
                initParams[paramIndex] = fcn->getY(iPts);
            }
        }
    }
    return initParams;
}

...

 initParams;
}

C.5: Unpack parameters into the model

At this point, the Optimizer has given us new parameters to try. We need to use these parameters to change the Y-values of the splines. We must first 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 simulation

Now that we've updated the controls, 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);

...

(_duration);
manager.integrate(initState);

C.7: Construct the objective function value

Here's the brain of the optimization. We assign to f the value of the 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; it provides a convenient way to (in one step) update the log and weight 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;
}

The following functions must be 'const' because  because objectiveFunc() is is const; we cannot call member functions within a a const member  member function unless those functions are also also const. These functions are also defined inline (i.e., definition is inside the class definition) for efficiency: the compiler will 'paste' the function body into the place where the function is called.

double anteriorLegsDownTerm(
        double roll, double rollRate, double rollAccel) const
{
    if (_tool.get_anterior_legs_down_weight() != 0.0)
    {
        return processObjTerm("anterior_legs_down", _tool.get_anterior_legs_down_weight(),
                pow(roll - SimTK::Pi, 2) + _relw * pow(rollRate, 2) + _relw * pow(rollAccel, 2));
    }
    return 0.0;
}
double posteriorLegsDownTerm(
        double twist, double twistRate, double twistAccel) const
{
    if (_tool.get_posterior_legs_down_weight() != 0.0)
    {
        return processObjTerm("posterior_legs_down", _tool.get_posterior_legs_down_weight(),
                pow(twist - 0.0, 2) + _relw * pow(twistRate, 2) + _relw * pow(twistAccel, 2));
    }
    return 0.0;
} 

double sagittalSymmetryTerm(
        double hunch, double pitch, double pitchRate, double pitchAccel) const
{
    if (_tool.get_sagittal_symmetry_weight() != 0.0)
    {
        return processObjTerm("sagittal_symmetry", _tool.get_sagittal_symmetry_weight(),
                pow(hunch + 2 * pitch, 2) + _relw * pow(pitchRate, 2) + _relw * pow(pitchAccel, 2));
    }
    return 0.0;
}

C.8: Update the log and outputs

We now return to objectiveFunc(), right after the code where we've assembled  f.  We update the log with the new objective function value , and might print an updated model file and updated control functions/controls.

_lastCallWasBestYet = _thisCallIsBestYet;
_thisCallIsBestYet = f <= _objectiveFcnValueBestYet;
if (_thisCallIsBestYet) _objectiveFcnValueBestYet = f;
if (_objectiveCalls % _outputPeriod == 0)
    _optLog << " objfcn " << f << " objfcn_best_yet " << _objectiveFcnValueBestYet << std::endl;
// If this is the best yet (i.e., smallest objective function value),
// save a copy of the splines.
if (_thisCallIsBestYet)
{
    for (unsigned int iFcn = 0; iFcn < _splinesBestYet.size(); iFcn++)
        _splinesBestYet[iFcn] = *_splines[iFcn];
}
// If this is worse, print out the best-yet splines and current model
// (NOTE: current model is NOT "best yet", but the idea is that it will
// be close). Also print out nondimensionalized splines so that they can
// be used directly as input for a subsequent optimization.
if (_lastCallWasBestYet && !_thisCallIsBestYet)
{
    printBestYetPrescribedControllerFunctionSet(
            _name + "_best_yet_parameters.xml");
    printBestYetPrescribedControllerFunctionSet(
            _name + "_best_yet_parameters_nomdim.xml", true);
    printModel(_name + "_best_yet.osim");
}
// Print out to the terminal/console every so often.
if (_objectiveCalls % _outputPeriod == 0)
    std::cout << "Objective call # " << _objectiveCalls << std::endl;
return 0;

C.9: Miscellaneous functions

There Here 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.

/**
 * @brief Prints the current model to the given filename.  The file will
 * be located in the directory containing the log file for this
 * optimization run.
 * */
void printModel(std::string filename) const
{
    _cat.print(_name + "/" + filename);
}
/**
 * @brief Serializes the current set of functions used for the actuators.
 * The file will be located in the directory containing the log file
 * for this optimization run.
 * */
void printPrescribedControllerFunctionSet(std::string filename) const
{
    // Create the FunctionSet that we'll then serialize.
    OpenSim::FunctionSet fset;
    fset.setSize(_splines.size());
    for (unsigned int iFcn = 0; iFcn < _splines.size(); iFcn++)
    {
        fset.insert(iFcn, *_splines[iFcn]);
    }
    fset.print(_name + "/" + filename);
}
/**
 * @brief Serializes the set of functions associated with the best-yet
 * value of the objective function. The file will be located in the
 * directory containing the log file for this optimization run.
 *
 * @param (nondimensionalize) Divide spline values by minControl for the
 * actuator, if value and minControl are negative, and by maxControl, if
 * value is positive. Once nondimensionalized, this FunctionSet can be used
 * as initial parameters.
 * */
void printBestYetPrescribedControllerFunctionSet(std::string filename,
        bool nondimensionalize=false) const
{
    using OpenSim::SimmSpline;
    // Create the FunctionSet that we will then serialize.
    OpenSim::FunctionSet fset;
    fset.setSize(_splinesBestYet.size());
    for (unsigned int iFcn = 0; iFcn < _splinesBestYet.size(); iFcn++)
    {
        fset.insert(iFcn, _splinesBestYet[iFcn]);
        if (nondimensionalize)
        { // Go through each y-value and nondimensionalize based on its sign.
            // Max/min for all spline points for the i-th actuator.
            double minControl = _cat.getActuators().get(iFcn).getMinControl();
            double maxControl = _cat.getActuators().get(iFcn).getMaxControl();
            for (int iPts = 0; iPts < _numOptimSplinePoints; iPts++)
            {
                SimmSpline * fcn = dynamic_cast<SimmSpline *>(&fset.get(iFcn));
                double dimValue = fcn->getY(iPts);
                double nonDimValue = 0;
                if (dimValue < 0 && minControl < 0)
                {
                    nonDimValue = -dimValue / minControl;
                }
                else
                {
                    nonDimValue = dimValue / maxControl;
                }
                fcn->setY(iPts, nonDimValue);
            }
        }
    }
    fset.print(_name + "/" + filename);
}

C.10: Member variables

The mutable keyword lets us modify a member variable inside const member functions, such as objectiveFunc() above.

// 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.

 printing of best yet actuation:
mutable bool _lastCallWasBestYet;
mutable bool _thisCallIsBestYet; 

Getting started with optimization and reproducing our results

We've given you some setup/input files to for the optimization to get you started and to reproduce reproducing our results. Enter a terminal or command prompt, navigate to a the directory containing the optimize executable and the counterrotation_setup.xml file. The following will start an optimization. It should take under 2 hours to complete.

...

Compare your results to those in the 'References' folder of the 'Tutorial code' codedirectory. TODO. You can do the same repeat this for the second optimization, if counterrotation_variableinertia_setup.xml and counterrotation_variableinertia_initial_parameters.xml are in the directory.

...

Remember that it's the setup file that refers to the initial parameters file. This second optimization should take about the same amount of time to complete.

Enjoy playing around with the code. We hope you've learned something useful!