Optimization Studies in Python

The optimization study introduced in the preceding lesson used AnyBody’s builtin facilities for optimizing. Sometimes that is not enough, either because the objective functions depends on data that is not easily included in AnyBody, or because a different solver is needed.

In this tutorial we use an external optimizer together with AnyBody. The example is based on the model from the previous lesson but uses an optimizer from the Scipy python library. The same could also have been achived with other optimization frameworks (like NLopt, or languages like MatLab). The example is written and run in Jupyter Lab, but the Python script can also be run from other environments.

Example Python Script

First we show the full Python script used in this tutorial. Afterwards, we will go through it and explain the different sections in the file. The first part below is the Python script, followed by the corresponding results from Jupyter Lab.

import math

import scipy
from anypytools import AnyPyProcess
from anypytools.macro_commands import Load, OperationRun, Dump, SetValue


def run_model(saddle_height, saddle_pos, silent=False):
    """ Run the AnyBody model and return the metabolism results """
    
    macro = [
        Load("BikeModel2D.main.any"),
        SetValue("Main.BikeParameters.SaddleHeight", saddle_height),
        SetValue("Main.BikeParameters.SaddlePos", saddle_pos),
        OperationRun("Main.Study.InverseDynamics"),
        Dump("Main.Study.Output.Pmet_total"),
        Dump("Main.Study.Output.Abscissa.t"),
    ]
    app = AnyPyProcess(silent=silent)
    results = app.start_macro(macro)
    return results[0]


def objfun(x):
    """ Calculate the objective function value """
    
    saddle_height = x[0]
    saddle_pos = x[1]
    result = run_model(saddle_height, saddle_pos, silent=True)

    if "ERROR" in result:
        raise ValueError("Failed to run model")
    
    # Integrate Pmet_total
    pmet = scipy.integrate.trapezoid(result["Pmet_total"], result["Abscissa.t"])

    return float(pmet)


def seat_distance_constraint(designvars):
    """ Compute contraint value which must be larger than zero """
    
    return math.sqrt(designvars[0] ** 2 + designvars[1] ** 2) - 0.66


CONSTRAINT = {"type": "ineq", "fun": seat_distance_constraint}
BOUNDS = [(0.61, 0.69), (-0.20, -0.05)]
INITIAL_GUESS = (0.68, -0.15)

solution = scipy.optimize.minimize(
    objfun, INITIAL_GUESS, constraints=CONSTRAINT, bounds=BOUNDS, method="SLSQP"
)

print(solution)

 message: Optimization terminated successfully
 success: True
  status: 0
     fun: 503.5763438548318
       x: [ 6.501e-01 -1.139e-01]
     nit: 12
     jac: [ 3.944e+01 -6.962e+00]
    nfev: 45
    njev: 12

Getting Started: Development Environment Setup

You can download the code for this tutorial as an interactive Python (Jupyter) notebook. Download and extract this zip-file.

Running the notebook requires a Python environment with the necessary libraries installed (jupyter, anypytools, scipy).

See the following section for easy instructions on how to easily install these libraries.

Installing Python, Jupyter, and AnyPyTools

The absolute easiest way to get the notebook running is with the Pixi package manager.

Simply run the following command in a powershell terminal to install Pixi (or see here for altinative installation methods):

iwr -useb https://pixi.sh/install.ps1 | iex

Then open a command prompt in the folder you extracted from the zip file above. Then run the following command to start the notebook:

C:\path-to-folder> pixi run jupyter lab lesson3.ipynb

pixi will automatically download and install Python along with all the necessary libararies in virtual environment. Then it will start Jupyter Lab in your web browser.

Importing the Necessary Libraries

The first part of the code is the import statements. They include the libraries which is used by the code:

import math
import scipy

from anypytools import AnyPyProcess
from anypytools.macro_commands import Load, OperationRun, Dump, SetValue

Running a Model From Python

For the external optimizers to work, we need a way to run AnyBody models from Python and record the results of the simulations. We will create a function called run_model() to do this. There are more information on how to do this in the documentaion for AnyPyTools. So here we will just show how the code for this function looks and not discuss its details. Executing the code cell in Jupyter Lab should give the following output:

def run_model(saddle_height, saddle_pos, silent=False):
    #Run the AnyBody model and return the metabolism results
    macro = [
        Load("BikeModel2D.main.any"),
        SetValue("Main.BikeParameters.SaddleHeight", saddle_height),
        SetValue("Main.BikeParameters.SaddlePos", saddle_pos),
        OperationRun("Main.Study.InverseDynamics"),
        Dump("Main.Study.Output.Pmet_total"),
        Dump("Main.Study.Output.Abscissa.t"),
    ]
    app = AnyPyProcess(silent=silent)
    results = app.start_macro(macro)
    return results[0]


result = run_model(0.66, -0.16)
print(result["Main.Study.Output.Pmet_total"])

Completed: 1
[  17.20903341   73.49291834  209.58490241  379.67002659  559.57715608
  736.92126247  901.88875426 1045.75303378 1162.65470516 1248.32088806
 1299.79539032 1315.38241529 1294.6947524  1238.68684947 1149.59584772
 1030.78784505  886.60667952  722.43408728  547.1840971   368.64175002
  198.07134668   53.41928909   25.84379129   30.60376508   23.17442367
   24.30809055  139.3209062   292.35610808  469.73382854  649.02576552
  821.74094457  977.02863522 1108.05435136 1209.79739513 1278.65973442
 1312.31195028 1309.70451022 1271.1212895  1198.17227557 1093.68215448
  961.51890402  806.51623776  634.74029158  458.00117565  280.40563034
  121.30841553   21.54859903   28.97200722   26.82989147   17.2090334 ]

In generel terms the function run_model() takes the saddle_height and saddle_pos as input and returns the total metabolism as output Pmet_total from the AnyBody model. It begins by loading the AnyBody model file BikeModel2D.main.any and sets the saddle_height and saddle_pos parameters. The function then runs the inverse dynamics study and finally returns the Pmet_total values for each time step.

As we expected the above output contains the total metabolism of the muscles, the Main.Study.Output.Pmet_total value, for each timestep in our model.

Defining the Objective Function

The next step is to define the objective function called objfun(). The objective function should take a list of design values as input and return the objective function value. In Lesson 2 the objective function was the time integral of the metabolism variable. So we will define it similarly here using Scipy’s numpy.trapz() function.

def objfun(x):
    #Calculate the objective function value
    saddle_height = x[0]
    saddle_pos = x[1]
    result = run_model(saddle_height, saddle_pos, silent=True)

    if "ERROR" in result:
        raise ValueError("Failed to run model")
    
    # Integrate Pmet_total
    pmet = scipy.integrate.trapezoid(result["Pmet_total"], result["Abscissa.t"])

    return float(pmet)


result_pmet = objfun([0.66, -0.16])
print(result_pmet)

505.3293995327723

Note

The function also checks the results for errors reported by AnyBody and raises a ValueError exception if that happens. There could be ways to handle errors without failing but it is important to handle model failures, otherwise they may go unnoticed or mess with the optimization results.

Executing the code cell returns the time integral of the Pmet_total variable as a single value, and we are now ready to define the optimization process. The interpretation of this number is the total metabolism of the muscles over the entire simulation time for the given values of saddle_height and saddle_pos.

Setting Up the Optimization Study

We wrap things up by creating a function, similar to what we did in AnyBody, as well as defining the constraints, bounds and initial guess for the design variables.

def seat_distance_constraint(x):
     #Compute contraint value which must be larger than zero
     return (math.sqrt(x[0] ** 2 + x[1] ** 2) - 0.66)


constraints = {"type": "ineq", "fun": seat_distance_constraint}
bounds = [(0.61, 0.69), (-0.20, -0.05)]
initial_guess = (0.68, -0.15)

Finally, we call the scipy.optimize.minimize function and run the optimizer. This function minimizes the objective function while respecting the constraints and bounds. In this case we used the SLSQP algorithm.

Note

The documentation scipy.optimize.minimize() has more information on bounds, contraints, tolerances, methods etc.

solution = scipy.optimize.minimize(
    objfun, initial_guess, constraints=constraints, bounds=bounds, method="SLSQP"
)

print(solution)

 message: Optimization terminated successfully
 success: True
  status: 0
     fun: 503.5763438548318
       x: [ 6.501e-01 -1.139e-01]
     nit: 12
     jac: [ 3.944e+01 -6.962e+00]
    nfev: 45
    njev: 12

Executing the code cell should give the above output, indicating that the optimization study was successfully completed. The optimal values for the design variables are determined to be saddle_height = 0.6501 and saddle_pos = -0.1139 resulting in a value of 503.58 for the objective function.

And there we have it! You have now used the optimizer from the Scipy Python library to optimize the AnyBody model. You can now take advantage of the many different algorithms and settings available for scipy.optimize.minimize(). It is also possible to use a different package or customize our own algorithm, constraints etc.

The possibilities are practically endless.

For more information regarding the AnyPyTools python package follow this link. or check out this webcast.