サッカーボールを強く蹴るにはどの筋肉を鍛えたらよいか？（２）

import opensim as osim
import numpy as np
import time
import cma


def getSpeed():
    data = np.genfromtxt('kick/leg6dof9musc_knee_stop_Kinematics_u.sto', dtype='float', names=True, autostrip=True, skip_header=10)
    return np.max(data['ball_tx'])


# OBJECTIVE FUNCTION
# Runs a forward simulation using the initial conditions specified in the
# candidate solution vector (candsol) and computes the corresponding
# objective function value (i.e. the maximum speed of the ball).
def succer_kicking_objective_function(candsol):
    global all_max_speeds, all_candsols

    model = osim.Model('SoccerKickingModel.osim')

    # Set the initial max isometric forces
    model.getMuscles().get('bifemlh_r').setMaxIsometricForce(candsol[0])
    model.getMuscles().get('rect_fem_r').setMaxIsometricForce(candsol[1])
    model.getMuscles().get('soleus_r').setMaxIsometricForce(candsol[2])
    model.getMuscles().get('tib_ant_r').setMaxIsometricForce(candsol[3])

    # Simulate
    tool = osim.ForwardTool('runFD.xml', True, False)
    tool.setModel(model)
    tool.run()

    x = getSpeed()
    # Store the candidate solution and the distance traveled.
    all_candsols.append(candsol)
    all_max_speeds.append(x)
    print(candsol)
    print(f'Max speed: {x} m/s')

    # To maximize max speed, minimize its negative.
    J = -x
    return (J)

# MAIN
# Perform an optimization using cma with the above objective function.
global all_max_speeds, all_candsols
all_max_speeds = []
all_candsols = []

# Load OpenSim model.
model = osim.Model('SoccerKickingModel.osim')

# Create candidate solution vector. If the optimizer is working
# correctly, it should increase this value to 9000.0 (the upper bound); see the
# penalty calculation in the objective function.
candsol = []
candsol.append(model.getMuscles().get('bifemlh_r').getMaxIsometricForce())
candsol.append(model.getMuscles().get('rect_fem_r').getMaxIsometricForce())
candsol.append(model.getMuscles().get('soleus_r').getMaxIsometricForce())
candsol.append(model.getMuscles().get('tib_ant_r').getMaxIsometricForce())

# Optimize.
t_start = time.time()
# For a description of arguments to fmin(), run cma.CMAOptions() or see
# http://cma.gforge.inria.fr/apidocs-pycma/cma.evolution_strategy.html#fmin
x, es = cma.fmin(succer_kicking_objective_function, candsol, 200.0,
                 options = {'popsize':20, 'tolfun':1e-3, 'tolx':1e-3,
                            'bounds': [0.0, 9000.0]})

t_elapsed = time.time() - t_start
print(f'Elapsed time: {t_elapsed} seconds')

# Find the best solution.
max_speed = max(all_max_speeds)
print(f'Best max speed: {max_speed} m/s')
idx = all_max_speeds.index(max_speed)
bestsol = all_candsols[idx]
print(bestsol)