|
| 1 | +# ----------------------------------------------------------------------- # |
| 2 | +# The OpenSim API is a toolkit for musculoskeletal modeling and # |
| 3 | +# simulation. See http://opensim.stanford.edu and the NOTICE file # |
| 4 | +# for more information. OpenSim is developed at Stanford University # |
| 5 | +# and supported by the US National Institutes of Health (U54 GM072970, # |
| 6 | +# R24 HD065690) and by DARPA through the Warrior Web program. # |
| 7 | +# # |
| 8 | +# Copyright (c) 2005-2019 Stanford University and the Authors # |
| 9 | +# Author(s): Thomas Uchida, Akshaykumar Patel, James Dunne # |
| 10 | +# Contributor(s): Andrew Miller, Jeff Reinbolt, Ajay Seth, Dan Jacobs, # |
| 11 | +# Chris Dembia, Ayman Habib, Carmichael Ong # |
| 12 | +# # |
| 13 | +# Licensed under the Apache License, Version 2.0 (the "License"); # |
| 14 | +# you may not use this file except in compliance with the License. # |
| 15 | +# You may obtain a copy of the License at # |
| 16 | +# http://www.apache.org/licenses/LICENSE-2.0. # |
| 17 | +# # |
| 18 | +# Unless required by applicable law or agreed to in writing, software # |
| 19 | +# distributed under the License is distributed on an "AS IS" BASIS, # |
| 20 | +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or # |
| 21 | +# implied. See the License for the specific language governing # |
| 22 | +# permissions and limitations under the License. # |
| 23 | +# ----------------------------------------------------------------------- # |
| 24 | + |
| 25 | +# This example demonstrates how to run an optimization in Python using the |
| 26 | +# cma package (see References, below). The model is a dynamic walker (see |
| 27 | +# the analogous Matlab example files for details); the optimizer adjusts the |
| 28 | +# initial forward velocity of the pelvis and the initial angles and angular |
| 29 | +# velocities of the knees and hips to maximize travel distance in 10 seconds. |
| 30 | +# Note that the OpenSim model 'dynamic_walker_example_model.osim' must be in |
| 31 | +# the same folder as this script and the cma package must have been installed. |
| 32 | +# As configured, the optimizer takes about 3 minutes to run on a standard |
| 33 | +# laptop (Intel i7, Windows 10, Python 2.7). For demonstration purposes only. |
| 34 | +# |
| 35 | +# References for cma package: |
| 36 | +# [1] Project page -- https://pypi.org/project/cma/ |
| 37 | +# [2] Implementation -- https://github.com/CMA-ES/pycma |
| 38 | +# [3] Documentation -- http://cma.gforge.inria.fr/apidocs-pycma/cma.html |
| 39 | + |
| 40 | +import opensim as osim |
| 41 | +import time |
| 42 | +import cma |
| 43 | + |
| 44 | + |
| 45 | +# OBJECTIVE FUNCTION |
| 46 | +# Runs a forward simulation using the initial conditions specified in the |
| 47 | +# candidate solution vector (candsol) and computes the corresponding |
| 48 | +# objective function value (i.e. the final location of the pelvis). |
| 49 | +def walker_simulation_objective_function(candsol): |
| 50 | + global model, initial_state, all_distances, all_candsols |
| 51 | + |
| 52 | + # Set the initial hip and knee angles. |
| 53 | + initial_state.updQ()[3] = candsol[0] # left hip |
| 54 | + initial_state.updQ()[4] = candsol[1] # right hip |
| 55 | + initial_state.updQ()[5] = candsol[2] # left knee |
| 56 | + initial_state.updQ()[6] = candsol[3] # right knee |
| 57 | + |
| 58 | + # Set the initial forward velocity of the pelvis. |
| 59 | + vx0 = candsol[4] # needed to compute the objective function, below |
| 60 | + initial_state.updU()[1] = vx0 |
| 61 | + |
| 62 | + # Set the initial hip and knee angular velocities. |
| 63 | + initial_state.updU()[3] = candsol[5] # left hip |
| 64 | + initial_state.updU()[4] = candsol[6] # right hip |
| 65 | + initial_state.updU()[5] = candsol[7] # left knee |
| 66 | + initial_state.updU()[6] = candsol[8] # right knee |
| 67 | + |
| 68 | + # Simulate. |
| 69 | + manager = osim.Manager(model) |
| 70 | + manager.initialize(initial_state) |
| 71 | + manager.integrate(10.0) |
| 72 | + |
| 73 | + # Get the final location of the pelvis in the X direction. |
| 74 | + st = manager.getStatesTable() |
| 75 | + dc = st.getDependentColumn('/jointset/PelvisToPlatform/Pelvis_tx/value') |
| 76 | + x = dc[ dc.nrow()-1 ] |
| 77 | + # Store the candidate solution and the distance traveled. |
| 78 | + all_candsols.append(candsol) |
| 79 | + all_distances.append(x) |
| 80 | + print('Distance traveled: %f meters' % (x)) |
| 81 | + |
| 82 | + # To maximize distance, minimize its negative. Also penalize candidate |
| 83 | + # solutions that increase the initial pelvis velocity beyond 2 m/s (to |
| 84 | + # avoid simply launching the model forward). This could also be done by |
| 85 | + # adding a hard constraint. |
| 86 | + k = 10.0 |
| 87 | + penalty1 = k*(vx0**2.0) if vx0 < 0 else 0.0 # lower bound: 0 m/s |
| 88 | + penalty2 = k*((vx0-2.0)**2.0) if vx0 > 2 else 0.0 # upper bound: 2 m/s |
| 89 | + J = -x + penalty1 + penalty2 |
| 90 | + return (J) |
| 91 | + |
| 92 | +# MAIN |
| 93 | +# Perform an optimization using cma with the above objective function. The |
| 94 | +# final model will be saved as 'dynamic_walker_example_model_optimized.osim'. |
| 95 | +global model, initial_state, all_distances, all_candsols |
| 96 | +all_distances = [] |
| 97 | +all_candsols = [] |
| 98 | + |
| 99 | +# Load OpenSim model. |
| 100 | +model = osim.Model('dynamic_walker_example_model.osim') |
| 101 | + |
| 102 | +# Create the underlying computational system. Note that we reuse the State |
| 103 | +# in the objective function to avoid the high computational cost of calling |
| 104 | +# initSystem() before every simulation. |
| 105 | +initial_state = model.initSystem() |
| 106 | + |
| 107 | +# Create candidate solution vector. The (arbitrary) initial guess for the |
| 108 | +# initial forward velocity of the pelvis is 0.1. If the optimizer is working |
| 109 | +# correctly, it should increase this value to 2.0 (the upper bound); see the |
| 110 | +# penalty calculation in the objective function. |
| 111 | +candsol = [] |
| 112 | +candsol.append(model.getCoordinateSet().get('LHip_rz').getDefaultValue()) |
| 113 | +candsol.append(model.getCoordinateSet().get('RHip_rz').getDefaultValue()) |
| 114 | +candsol.append(model.getCoordinateSet().get('LKnee_rz').getDefaultValue()) |
| 115 | +candsol.append(model.getCoordinateSet().get('RKnee_rz').getDefaultValue()) |
| 116 | +candsol.append(0.1) |
| 117 | +candsol.append(model.getCoordinateSet().get('LHip_rz').getDefaultSpeedValue()) |
| 118 | +candsol.append(model.getCoordinateSet().get('RHip_rz').getDefaultSpeedValue()) |
| 119 | +candsol.append(model.getCoordinateSet().get('LKnee_rz').getDefaultSpeedValue()) |
| 120 | +candsol.append(model.getCoordinateSet().get('RKnee_rz').getDefaultSpeedValue()) |
| 121 | + |
| 122 | +# Optimize. |
| 123 | +t_start = time.time() |
| 124 | +# For a description of arguments to fmin(), run cma.CMAOptions() or see |
| 125 | +# http://cma.gforge.inria.fr/apidocs-pycma/cma.evolution_strategy.html#fmin |
| 126 | +result = cma.fmin(walker_simulation_objective_function, candsol, 0.5, |
| 127 | + options = {'popsize':20, 'tolfun':1e-3, 'tolx':1e-3, |
| 128 | + 'maxfevals':100}) |
| 129 | +t_elapsed = time.time() - t_start |
| 130 | +print('Elapsed time: %f seconds' % (t_elapsed)) |
| 131 | + |
| 132 | +# Find the best solution. |
| 133 | +max_distance = max(all_distances) |
| 134 | +print('Best distance: %f meters' % (max_distance)) |
| 135 | +idx = all_distances.index(max_distance) |
| 136 | +bestsol = all_candsols[idx] |
| 137 | +print(bestsol) |
| 138 | + |
| 139 | +# Assign best solution to model and save. |
| 140 | +model.getCoordinateSet().get('LHip_rz').setDefaultValue(bestsol[0]) |
| 141 | +model.getCoordinateSet().get('RHip_rz').setDefaultValue(bestsol[1]) |
| 142 | +model.getCoordinateSet().get('LKnee_rz').setDefaultValue(bestsol[2]) |
| 143 | +model.getCoordinateSet().get('RKnee_rz').setDefaultValue(bestsol[3]) |
| 144 | +model.getCoordinateSet().get('Pelvis_tx').setDefaultSpeedValue(bestsol[4]) |
| 145 | +model.getCoordinateSet().get('LHip_rz').setDefaultSpeedValue(bestsol[5]) |
| 146 | +model.getCoordinateSet().get('RHip_rz').setDefaultSpeedValue(bestsol[6]) |
| 147 | +model.getCoordinateSet().get('LKnee_rz').setDefaultSpeedValue(bestsol[7]) |
| 148 | +model.getCoordinateSet().get('RKnee_rz').setDefaultSpeedValue(bestsol[8]) |
| 149 | +model.printToXML('dynamic_walker_example_model_optimized.osim') |
0 commit comments