Optimization Studies in Python

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.58, 0.64), (-0.21, -0.13)]
INITIAL_GUESS = (0.64, -0.19)

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

print(solution)

c:\Users\jha\AppData\Local\miniforge3\Lib\site-packages\tqdm\auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

 message: Optimization terminated successfully
 success: True
  status: 0
     fun: 516.4476449033699
       x: [ 6.338e-01 -1.841e-01]
     nit: 11
     jac: [ 6.100e+01 -1.719e+01]
    nfev: 43
    njev: 11

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.64, -0.19)
print(result["Main.Study.Output.Pmet_total"])

  0%|          | 0/1 [00:00<?, ?it/s]

100%|██████████| 1/1 [00:05<00:00,  5.36s/it]

Completed: 1
[  19.69795204   48.38931068  114.62742585  190.30259508  272.23849904
  365.13428457  459.07706787  551.23440015  641.61902196  728.58957213
  811.17757193  890.19447134  963.94613992 1031.80854769 1093.40272558
 1148.40311428 1196.52136611 1237.49814426 1271.10123307 1297.12804346
 1315.41061895 1325.82141813 1328.27839346 1322.74816782 1309.24641078
 1287.83484906 1258.61473836 1221.71711994 1177.29084426 1125.49022102
 1066.46528969 1000.35907566  927.31766095  847.52009187  761.23537781
  669.6079932   575.73062944  476.68498989  373.58217557  269.32737802
  167.96959607   71.77796079   25.64607135   31.61062858   38.82607851
   42.10641735   41.28790794   41.41400061   37.74086717   30.50034078
   19.80630727   48.41664074  114.68667561  190.4530785   272.46399335
  365.34221337  459.27237225  551.40625066  641.77964993  728.74310402
  811.32732136  890.34595419  964.11540493 1031.99409392 1093.60250066
 1148.61467311 1196.74199203 1237.72493786 1271.33117373 1297.35802881
 1315.63748708 1326.04195733 1328.48934435 1322.94622267 1309.42821477
 1287.99700834 1258.75383915 1221.82976248 1177.37369405 1125.5400779
 1066.47917232 1000.33431496  927.25199105  847.41172239  761.08302833
  669.29994554  575.43083185  476.42218384  373.3431837   269.10122481
  167.74479086   71.53701528   25.36818809   31.21084533   38.27516228
   41.49378618   41.17405609   41.29571624   37.62119175   30.38308883
   19.69810172]

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.64, -0.19])
print(result_pmet)

517.5231021676813

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.58, 0.64), (-0.21, -0.13)]
initial_guess = (0.64, -0.19)

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: 516.4476449033699
       x: [ 6.338e-01 -1.841e-01]
     nit: 11
     jac: [ 6.100e+01 -1.719e+01]
    nfev: 43
    njev: 11

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.6338 and saddle_pos = -0.1841 resulting in a value of 516.45 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.

Contents