Optimizing Parameters in SCRIMMAGE¶

Mission Setup¶

Here we investigate a predator-prey scenario and discuss how to use SCRIMMAGE to optimize model parameters. There are many ways to do this but here we focus on optimization using Gaussian Processes. For background on Gaussian Processes see [Rasmussen]. For a similar discussion on surrogate modeling with an additional discussion of sampling plans see [Forrester]. When the parameter space is larger than the current application it is recommended to start with a latin hypercube sampling technique. We do not discuss this here but note that it is available elsewhere in SCRIMMAGE (see Multiple Runs on a Local Computer for details).

First let’s run the scenario:

scrimmage missions/predator_prey_boids.xml

The code for the optimization can be found in scripts/optimization_tutorial.py. To run it you will need to install the python dependences:

sudo apt install python3-sklearn python3-pandas python-numpy
sudo pip3 install bayesian-optimization==0.6.0

You should see about 50 blue vehicles (the prey) and 1 red vehicle (the predator). The blue vehicles are running a version of the Reynolds-Boids model [Boids]. The predators attempt to get close enough to the prey to capture them and in the process get points. There are a couple of parameters that affect the efficacy of a predator, namely speed and capture radius. Both factors generally yield better performance when they are larger. To create a more interesting scenario we assume that the two factors are coupled. In particular, we assume that as the speed is increased there is a proportionate decrease in capture radius. We want to know what the optimal speed and capture radius are under this constraint.

To do this note that the predator_prey_boids.xml file uses metrics plugin called SimpleCaptureMetrics that is responsible for writing to a summary.csv file in the log directory. See Create a Metrics Plugin for details of how metrics work. We will be changing the max_speed of the predator autonomy model in order to maximize this output.

To make the problem somewhat more interesting we will introduce noise into the simulation by randomizing the starting locations (see the variance_x, variance_y, and variance_z tags). This means that we need to repeat simulations to get an estimate of the expected value of the output. Gaussian Processes can interpolate with noise by adding a constant to the correlation matrix. We could go with this approach since we can estimate the variance of an output given the repeated runs. However, it is simpler to just let the library handle this for us. This is the purpose of the WhiteKernel below:

def main():
    repeats = 16
    cores = 8
    mission = find_mission('predator_prey_boids.xml')
    nominal_capture_range = 5
    nominal_speed = 35
    kappa = 5       # higher is more exploration, less exploitation

    num_samples = 20
    low_speed = 10
    high_speed = 200
    num_explore_points = 10

    def _run(max_speed):
        num = len(glob.glob('.optimize*'))
        return run(repeats, cores, mission, num,
                   nominal_capture_range, nominal_speed, max_speed)

    pbounds = {'max_speed': (10, 200)}
    bo = bayes_opt.BayesianOptimization(_run, pbounds)

    init_explore_points = \
        {'max_speed':  np.linspace(low_speed, high_speed, num_explore_points)}
    bo.explore(init_explore_points)

    gp_params = {'kernel': 1.0 * RBF() + WhiteKernel()}
    bo.maximize(init_points=1, n_iter=num_samples - num_explore_points,
                kappa=kappa, **gp_params)


if __name__ == '__main__':
    main()

This code essentially provides the function to optimize _run to the Gaussian Process. _run captures some variables so that it is only a function of the optimization variable and calls run:

def run(repeats, cores, mission, num,
        nominal_capture_range, nominal_speed, max_speed):
    """Runs the missions and aggragate the data from each summary.csv"""
    out_dir, out_mission = \
        create_mission(mission, num, nominal_capture_range,
                       nominal_speed, max_speed)

    parallel(repeats, out_mission, cores)

    # Finds and aggragates the score
    files = [os.path.expanduser(os.path.join(out_dir, d, 'summary.csv'))
             for d in os.listdir(os.path.expanduser(out_dir))]

    scores = []
    for f in files:
        try:
            if not os.path.exists(f):
                print("{} does not exists".format(f))
            scores.append(pd.read_csv(f)['score'].sum())
        except (OSError, IndexError):
            scores.append(0)
    score = np.array(scores).mean()

    return score

This function adjusts the predator_prey_boids.xml file so that the constraint is satisfied. In particular, it sets the max_speed for the predator and makes the corresponding adjustment to capture_radius. It then calls a helper function parallel that runs SCRIMMAGE instances locally (see [Parallel] and Multiple Runs on a Local Computer for more details). For larger problems where grid engine is available for a cluster, batch jobs can instead call the scrimmage.qsub and scrimmage.wait_for_job functions. Here is the create_mission function:

def create_node(tag, text):
    """Create an xml node."""
    el = ET.Element(tag)
    el.text = "{:.2f}".format(text) if isinstance(text, float) else str(text)
    return el


def create_mission(mission, num, nominal_capture_range, nominal_speed, max_speed):
    """Modify the mission xml with custom parameters"""
    tree = ET.parse(mission)
    root = tree.getroot()

    # Removes the seed for each run
    seed_node = root.find('seed')
    if seed_node != None:
        root.remove(seed_node)

    run_node = root.find('run')
    run_node.attrib['enable_gui'] = "false"
    run_node.attrib['time_warp'] = "0"

    log_dir_node = root.find('log_dir')
    out_dir = os.path.join(log_dir_node.text, 'optimize' + str(num))
    log_dir_node.text = out_dir

    ratio = nominal_speed / max_speed
    capture_range = nominal_capture_range * ratio

    # Applies the max_speed and capture_range attributes to the Predator and SimpleCapture
    for entity_node in root.findall('entity'):
        autonomy_node = entity_node.find('autonomy')
        if autonomy_node.text == 'Predator':
            autonomy_node.attrib['max_speed'] = str(max_speed)
            autonomy_node.attrib['capture_range'] = str(capture_range)

    for interaction_node in root.findall('entity_interaction'):
        if interaction_node.text == 'SimpleCapture':
            interaction_node.attrib['capture_range'] = str(capture_range)

    out_mission = '.optimize' + str(num) + '.xml'
    tree.write(out_mission)

    return out_dir, out_mission

We can now run this file and get the following:

Initialization
-------------------------------------------
 Step |   Time |      Value |   max_speed |
| 05m05s |    0.52941 |     10.0000 |
| 04m15s |    1.50000 |     31.1111 |
| 04m04s |    6.88235 |     52.2222 |
| 04m09s |    6.20000 |     73.3333 |
| 04m06s |    6.11765 |     94.4444 |
| 04m11s |    5.52941 |    115.5556 |
| 04m08s |    6.29412 |    136.6667 |
| 04m06s |    5.11765 |    157.7778 |
| 04m06s |    6.05882 |    178.8889 |
| 04m10s |    4.17647 |    200.0000 |
| 04m07s |    5.88235 |     36.7489 |
Bayesian Optimization
-------------------------------------------
 Step |   Time |      Value |   max_speed |
| 04m08s |    6.62500 |    114.1056 |
| 04m09s |    5.81250 |    111.5516 |
| 04m08s |    6.64706 |    110.5069 |
| 04m06s |    7.76471 |     82.7220 |
| 04m08s |    7.35294 |     79.4004 |
| 04m09s |    6.41176 |     78.4301 |
| 04m07s |    6.11765 |     79.6248 |
| 04m08s |    7.05882 |     81.1437 |
| 04m06s |    6.58824 |     80.4494 |
| 04m09s |    6.35294 |     80.4652 |

The best speed found so far is 65.6019 (note that we could have had more exploration by setting kappa to something higher). We can either continue the search process with more or use a function minimizer to minimize the Gaussian Process to get a final estimate of optimal value.

Boids: Reynolds, Craig W. “Flocks, herds and schools: A distributed behavioral model.” ACM SIGGRAPH computer graphics. Vol. 21. No. 4. ACM, 1987.
Rasmussen: Rasmussen, Carl Edward. “Gaussian processes in machine learning.” Advanced lectures on machine learning. Springer, Berlin, Heidelberg, 2004. 63-71.
Forrester: Forrester, Alexander, and Andy Keane. Engineering design via surrogate modelling: a practical guide. John Wiley & Sons, 2008.
Parallel: Tange, Ole. “Gnu parallel-the command-line power tool.” The USENIX Magazine 36.1 (2011): 42-47.