{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# An Introduction to Optimisation via Simulation. \n",
"# How to choose the best setup for a system.\n",
"## A tutorial for the Simulation Workshop 2021\n",
"\n",
"**Dr Christine S.M Currie and Tom Monks**\n",
"\n",
"This tutorial introduces Ranking and Selection procedures for stochastic computer simulation. We focus on in difference zone and optimal computer budget allocation procedures. The procedures are:\n",
"\n",
"* KN and KN++\n",
"* OCBA and OCBA-m\n",
"\n",
"> The code implementing the Ranking and Selection procedures explored in this tutorial is available online with an MIT license https://github.com/TomMonks/sim-tools"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# run this first to install the OvS procedures so we can use them in the tutorial.\n",
"# !pip install sim-tools"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"from sim_tools.ovs.toy_models import (custom_gaussian_model,\n",
" gaussian_sequence_model,\n",
" random_gaussian_model,\n",
" ManualOptimiser)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Simulation Models\n",
"\n",
"To keep the focus of this tutorial on Ranking and Selection procedures and what they do, and how they perform in practice, the models we introduce are deliberately simple. The following sections describe how to create a a set of simulation models where the outputs are independent guassian distributions. There are three options:\n",
"\n",
"* the means of these outputs distributions follow a sequence (e.g. 1, 2, 3, 4, 5) variance is 1.0. \n",
"* the means and variances are user specified.\n",
"* the means and variances are randomly generated."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Creating a simple sequence of normal distributions\n",
"\n",
"A simple model can be created using the `guassian_sequance_model()` function. This creates a simulation model where the output of each design follows a $N - (\\mu_{i}, 1.0)$ :\n",
"\n",
"The function accepts three keyword arguments:\n",
"\n",
"* **start** - int, the first mean in the sequence (inclusive)\n",
"* **end** - int, the last mean int the sequence (inclusive)\n",
"* **step** - int, the difference between mean i and mean i + 1.\n",
"\n",
"For example, the following code creates a simulation model with 10 designs with means 1 to 10 and unit variance.\n",
"\n",
"```python\n",
"model = guassian_sequence_model(1, 10)\n",
"```\n",
"\n",
"To create a simulation model with 5 designs where $\\mu_{i+1} - \\mu_i = 2 $ :\n",
"\n",
"```python\n",
"model = guassian_sequence_model(1, 10, step=2)\n",
"```\n",
"\n",
"### Creating a custom model with known designs\n",
"Instead of a sequence of normal distributions with unit variance, it is possible to create a custom set of designs with varying variances. Use the `custom_guassian_model` function for this task. For example to create a custom set of designs: \n",
"\n",
"```python\n",
"means = [5, 8, 1, 2, 1, 7]\n",
"variances = [0.1, 1.2, 1.4, 0.3, 0.8]\n",
"\n",
"custom_model = custom_guassian_model(means, variances)\n",
"```\n",
"\n",
"The following code demonstrates how to create a sequence of 100 designs with variances that are 10% of the mean. \n",
"\n",
"```python\n",
"n_designs = 100\n",
"means = [i for i in range(n_designs+1)]\n",
"variances = [j*0.1 for j in range(n_designs+1)]\n",
"\n",
"custom_model = custom_guassian_model(means, variances)\n",
"```\n",
"\n",
"\n",
"### Creating a model with randomly sampled designs\n",
"\n",
"To create a model with a set of unknown designs (within a specified mean and variance tolerance) use \n",
"\n",
"```python\n",
"mean_low = 1.0\n",
"mean_high = 15.0\n",
"var_low = 0.1\n",
"var_high = 2.0\n",
"n_designs = 15\n",
"\n",
"model = random_guassian_model(mean_low, mean_high, var_low, var_high, n_designs)\n",
"```\n",
"\n",
"Where:\n",
"* **mean_low** - float, a lower bound on the means of the output distributions\n",
"* **mean_high** - float, an upper bound on the means\n",
"* **var_low** - float, a lower bound on the variance of the output distributions\n",
"* **var_high** - float, an upper bound on the variances.\n",
"* **n_designs** - int, the number of designs to create.\n",
"\n",
"---\n",
"\n",
"## A manual optimisation framework\n",
"\n",
"Before using the **Optimisation via Simulation** (OvS) procedures, it is recommended that you get a feel for the framework in which the OvS procedures operate. To do this we will create some models and explore them using a `ManualOptimiser`. This allows the user to run independent and multiple replications of the model yourself independent of any algorithm. The `ManualOptimiser` keeps track of the means, variances and number of replications run for each design.\n",
"\n",
"A `ManualOptimiser` object requires two parameters when it is created. \n",
"\n",
"* model - object, e.g. a model that is a sequence of normal distributions\n",
"* n_designs - the number of designs to be considered."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"manual_opt = ManualOptimiser(model=gaussian_sequence_model(1, 10),\n",
" n_designs=10)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* We can print the optimiser object to help us remember what parameters we set."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ManualOptimiser(model=BanditCasino(), n_designs=10, verbose=False)\n"
]
}
],
"source": [
"print(manual_opt)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Follow Law and Kelton's advice and run 3 initial replications of each design."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"manual_opt.simulate_designs(replications=3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The optimiser keeps track of the allocation of replications across each design. For efficiency it doesn't store each individual observation, but it does compute a running mean and variance.\n",
"\n",
"* Let's have a look at the replication allocation between and results of each design."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([3, 3, 3, 3, 3, 3, 3, 3, 3, 3], dtype=int32)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"manual_opt.allocations"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.7, 2.4, 4. , 3.7, 5. , 6. , 6.9, 8.8, 8.4, 9.8])"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.set_printoptions(precision=1) # this is a hack to view at 1 decimal place.\n",
"manual_opt.means"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.5, 1.9, 0.6, 0.8, 0.2, 0.9, 3.1, 0.5, 1. , 3.3])"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"manual_opt.vars"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's run 1 additional replication of the top 3 designs. \n",
"Note - in python arrays are **zero indexed**. This means that design 1 has index 0."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"manual_opt.simulate_designs(design_indexes=[7, 8, 9], replications=1)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([3, 3, 3, 3, 3, 3, 3, 4, 4, 4], dtype=int32)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"manual_opt.allocations"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.7, 2.4, 4. , 3.7, 5. , 6. , 6.9, 8.7, 8.7, 9.5])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"manual_opt.means"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.5, 1.9, 0.6, 0.8, 0.2, 0.9, 3.1, 0.4, 1. , 2.5])"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"manual_opt.vars"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Manual Optimisation of a model with unknown means."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Now have a go yourself**. This time create a model with random designs. Run as many replications of each design as you think is neccessary to make a decision about which design is best. Here we define best as the design with the largest mean value."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"#create random problem\n",
"rnd_model = random_gaussian_model(mean_low=5, mean_high=25, \n",
" var_low=0.5, var_high=2.5,\n",
" n_designs=10)\n",
"\n",
"#create manual optimiser for you to use.\n",
"manual_opt = ManualOptimiser(model=rnd_model, n_designs=10)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"#insert your code here...\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"# Indifference Zone Ranking and Selection"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Procedure **KN**\n",
"The first Ranking and Selection (R&S) algorithm we will explore is **Procedure KN** developed by Kim and Nelson. This is an *Indifference Zone* (IZ) approach to R&S. IZ procedures exploit the the idea that the performance of one or more of the sub-optimal designs are in fact so close to the best performing system that a user does not care which one is chosen (the decision maker is said to be *indifferent* to this choice). To do this a user must specify a quantity $\\delta$. Procedure KN provides a theorectical guarantee to select the best system or a system within $\\delta$ of the best with probability of $1 - \\alpha$.\n",
"\n",
"A key feature of KN is that it only estimates the variance of each design once. This happens after an initial number of replications (specified by the user) $n_0$. \n",
"\n",
"To run Kim and Nelson's R&S procedure KN, create an instance of `ovs.indifference_zone.KN`\n",
"\n",
"An object of type KN takes the following parameters:\n",
"\n",
"* **model** - a simulation model\n",
"* **n_designs** - int, the number of competing designs to compare\n",
"* **delta** - float, the indifference zone\n",
"* **alpha** - float, $PCS = 1-\\alpha$ (default=0.05)\n",
"* **n_0** - int, $n_0$ the number of initial replications (default=2)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"from sim_tools.ovs.indifference_zone import KN, KNPlusPlus"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Example 1. A simple sequence of normal distributions.\n",
"\n",
"The first problem we will test KN against is selecting the largest mean from a sequence of 10 normal distributions (means raning from 1 to 10) with unit variance. We want a $PCS = 0.95$ and are indifferent to designs that are within an average of 1.0 of the best. \n",
"\n",
"For this simple problem, we will follow Law's longstanding advice of setting $n_0 = 5$ i.e 5 initial replications of each design.\n",
"\n",
"#### Setting up KN"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"N_DESIGNS = 10\n",
"DELTA = 1.0\n",
"ALPHA = 0.05\n",
"N_0 = 5\n",
"\n",
"#first we create the simulation model. \n",
"model = gaussian_sequence_model(1, N_DESIGNS)\n",
"\n",
"#then we create the KN R&S object and pass in our parameters.\n",
"#note that this doesn't run the optimisation yet.\n",
"kn = KN(model=model, \n",
" n_designs=N_DESIGNS, \n",
" delta=DELTA, \n",
" alpha=ALPHA, \n",
" n_0=N_0)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"KN(n_designs=10, delta=1.0, alpha=0.05, n_0=5, obj=max)\n"
]
}
],
"source": [
"#quickly check if KN is parameterised as expected\n",
"print(kn)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Running KN\n",
"Now that we have created the KN object we can run the solver at any time we choose. To do this call the `KN.solve()` method. This method returns a design within the indifference zone $(1 - alpha$) of the time."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"best design\t[9]\n",
"allocations\t[ 5 5 5 5 7 5 5 5 13 14]\n",
"total reps\t69\n"
]
}
],
"source": [
"#this runs KN (you can run this cell multiple times if you wish)\n",
"best_design = kn.solve()\n",
"\n",
"#print out the results\n",
"print(f'best design\\t{best_design}')\n",
"print(f'allocations\\t{kn._allocations}')\n",
"print(f'total reps\\t{kn._allocations.sum()}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In general you should see that KN simulates the top 2-3 designs the most with the lower performing designs eliminated early on (possibly after the initial stage of replication)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Procedure **KN++**\n",
"\n",
"Kim and Nelson's KN++ procedure is an enhancement on KN. KN++ introduces an **update** step that recalculates the variance of the differences between designs. \n",
"\n",
"To run procedure KN++, create an instance of `ovs.indifference_zone.KNPlusPlus`\n",
"\n",
"### Example 2. Solving the simple problem using KN++\n",
"\n",
"We will run KN++ on the same simulation model as KN. In general you should see that KN++ uses less replication to find a design within the indifference zone."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"best design\t[9]\n",
"allocations\t[5 5 5 6 5 6 6 7 6 8]\n",
"total reps\t59\n"
]
}
],
"source": [
"knpp = KNPlusPlus(model=model, \n",
" n_designs=N_DESIGNS, \n",
" delta=DELTA, \n",
" alpha=ALPHA, \n",
" n_0=N_0)\n",
"\n",
"best_design = knpp.solve()\n",
"print(f'best design\\t{best_design}')\n",
"print(f'allocations\\t{knpp._allocations}')\n",
"print(f'total reps\\t{knpp._allocations.sum()}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise 1. Comparing KN and KN++ on a larger design.\n",
"\n",
"Your task is to compare KN and KN++ optimising a simulation model with 100 competing designs. These designs are a sequence of 100 normal distributions with unit variance. The objective is to find the design with the maximum mean.\n",
"\n",
"What do you notice about the replications required from each method?"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"#the code to set up the simulation model has been provided below.\n",
"N_DESIGNS = 100\n",
"DELTA = 1.0\n",
"ALPHA = 0.05\n",
"N_0 = 5\n",
"\n",
"#first we create the simulation model. \n",
"model = gaussian_sequence_model(1, N_DESIGNS)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"#your code goes here...\n",
"# hint 1. you need to create KN and KNPlusPlus objects, pass in the simulation model and run them. \n",
"# hint 2. try running the models a few times. \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise 2: A more difficult optimisation problem\n",
"\n",
"We will now explore how the IZ procedures perform on a slightly harder optimisation problem where the mean of the designs are closer together and have differing variances.\n",
"\n",
"The problem is: \n",
"\n",
"```python\n",
"means = [4, 4.1, 4.2, 4, 4.1, 4.3, 4, 4.1, 4.2, 4.2]\n",
"variances = [1, 1, 1, 0.1, 0.1, 10, 10, 10, 10, 0.1]\n",
"```\n",
"\n",
"The 'optimal' mean design is at index 5 (4.3). It is a very noisy design. We will set a $\\delta$ so that there are multiple designs within the IZ.\n",
"\n",
"Your task is to compare the performance of KN and KN++ given the following parameters:\n",
"\n",
"```python\n",
"N_DESIGNS = 10\n",
"N_0 = 10\n",
"DELTA = 0.15\n",
"ALPHA = 0.05\n",
"```\n",
"\n",
"We are also setting a random seed before running each procedure. This ensures that the first stage of replications returns the same samples.\n",
"\n",
"Questions: \n",
"* Do the methods return the 'optimal' design or another design within the indifference zone?\n",
"* Which procedure require the most replication effort? Does it vary?\n",
"* Are there any differences in designs selected and the number of replications needed? Try different values of `SEED`.\n",
"* What happens if you try difference values of $n_0$ e.g. 20 and 50?\n",
"* What happens if you change $\\alpha$ to 0.1? \n",
"\n",
"\n",
"\n",
"RUNTIME WARNING! This problem requires many more replications. It will take a few seconds for the procedures to return a design.
\n"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[2 5 8 9]\n"
]
}
],
"source": [
"# set problem parameters and create simulation model.\n",
"\n",
"N_DESIGNS = 10\n",
"N_0 = 20\n",
"DELTA = 0.15\n",
"ALPHA = 0.05\n",
"\n",
"#change the value of seed.\n",
"SEED = 999\n",
"\n",
"means = [4, 4.1, 4.2, 4, 4.1, 4.3, 4, 4.1, 4.2, 4.2]\n",
"variances = [1, 1, 1, 0.1, 0.1, 10, 10, 10, 10, 0.1]\n",
"\n",
"guass_model = custom_gaussian_model(means, variances)\n",
"\n",
"#print out which design indexes are in the indifference zone\n",
"print(np.where(4.3 - np.array(means) <= DELTA)[0])"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"best design\t[2]\n",
"allocations\t[ 415 314 41678 20 20 29 6697 41677 23238 20]\n",
"total reps\t114108\n"
]
}
],
"source": [
"# create an instance of KN and solve.\n",
"\n",
"kn = KN(model=guass_model, \n",
" n_designs=N_DESIGNS, \n",
" delta=DELTA, \n",
" alpha=ALPHA, \n",
" n_0=N_0)\n",
"\n",
"np.random.seed(SEED)\n",
"best_design = kn.solve()\n",
"\n",
"#print out the results\n",
"print(f'best design\\t{best_design}')\n",
"print(f'allocations\\t{kn._allocations}')\n",
"print(f'total reps\\t{kn._allocations.sum()}')"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"best design\t[5]\n",
"allocations\t[ 190 200 107 20 20 28167 8476 7204 28166 15383]\n",
"total reps\t87933\n"
]
}
],
"source": [
"# use KN++ on the same problem\n",
"knpp = KNPlusPlus(model=guass_model, \n",
" n_designs=N_DESIGNS, \n",
" delta=DELTA, \n",
" alpha=ALPHA, \n",
" n_0=N_0)\n",
"\n",
"# we use the same random seed to help with comparison.\n",
"np.random.seed(SEED)\n",
"best_design = knpp.solve()\n",
"\n",
"# print out the results\n",
"print(f'best design\\t{best_design}')\n",
"print(f'allocations\\t{knpp._allocations}')\n",
"print(f'total reps\\t{knpp._allocations.sum()}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Optimal Computing Budget Allocation (OCBA)\n",
"\n",
"We will first consider the scenario where the objective is find the optimal design. The distinguishing feature of OCBA procedures from IZ procedures a fixed budget formulation of the optimisation problem. I.e. achieve the best $PCS$ with the budget available.\n",
"\n",
"Single solution OCBA object can be created by importing `ovs.fixed_budget.OCBA`\n",
"\n",
"An object of type `OCBA` takes the following parameters:\n",
"\n",
"* **model** - a simulation model\n",
"* **n_designs** - int, the number of competing designs to compare\n",
"* **budget** - int, the total number of replications to allocate across designs\n",
"* **delta** - int, the incremental amount of replications to allocate at each round\n",
"* **n_0** - int, $n_0$ the number of initial replications (default=5)\n",
"* **obj** - str, 'min' if minimisation; 'max' if maximisation (default='min')"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"from sim_tools.ovs.fixed_budget import OCBA, OCBAM"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 3. Applying OCBA to the simple problem.\n",
"\n",
"We will first apply OCBA to the simple sequence of normal random variables that we used in Example 1. \n",
"\n",
"Initially we will set a simulation budget of 1000 replications. This is overkill, but makes it easier to see what OCBA is doing."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [],
"source": [
"N_DESIGNS = 10\n",
"\n",
"#remember that delta is not the same as in IZ!\n",
"#it is the amount of replication to allocate at each stage.\n",
"DELTA = 5\n",
"\n",
"#the intial stage will run 5 reps of each design. \n",
"N_0 = 5\n",
"\n",
"BUDGET = 1000\n",
"\n",
"#first we create the simulation model. \n",
"model = gaussian_sequence_model(1, N_DESIGNS)\n",
"\n",
"ocba = OCBA(model=model, \n",
" n_designs=N_DESIGNS, \n",
" budget=BUDGET, \n",
" delta=DELTA, \n",
" n_0=N_0, \n",
" obj='max')"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OCBA(n_designs=10, budget=1000, delta=5, n_0=5, obj=max)\n"
]
}
],
"source": [
"print(ocba)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"call the `solve()` method to run the optimisation. If needed, run OCBA a few times."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"best design\t9\n",
"allocations\t[ 5 11 5 5 15 5 5 101 426 422]\n",
"total reps\t1000\n"
]
}
],
"source": [
"best_design = ocba.solve()\n",
"\n",
"#print out the results\n",
"print(f'best design\\t{best_design}')\n",
"print(f'allocations\\t{ocba._allocations}')\n",
"print(f'total reps\\t{ocba._allocations.sum()}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"* Run the cell above a few times. You should see that OCBA is allocating most of the budget across the top 3 designs."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise 3. Exploring OCBA and $PCS$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"OCBA conceptualises optimisation as a fixed budget problem. Let's run an experiment where we estimate what $PCS$ is actually being achieved under a given budget. KN++ was reporting a range of 70-80 replications. Let's try 75.\n",
"\n",
"To help with this test you can import an Experiment object: `ovs.evaluation.Experiment`\n",
"\n",
"An `Experiment` object takes the following parameters:\n",
"\n",
"* **env** - the simulated environment aka the model.\n",
"* **procedure** - the method to test e.g. KNPlusPlus or OCBA\n",
"* **best_index** - the index of the best design in the simulation model (zero based).\n",
"* **objective** - 'max' or 'min'\n",
"* **replications** - the number of repeats of the experiment to run.\n",
"\n",
"**Your task is to run the experiment with different budgets. E.g. 60, 70 and 75. Remember that you cannot go lower than N_0 + DELTA**\n"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"from sim_tools.ovs.evaluation import Experiment"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"RUNTIME WARNING! The experiment will will take a few seconds to complete.
"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PCS=0.999\n"
]
}
],
"source": [
"#this is the budget you can vary.\n",
"BUDGET = 75\n",
"\n",
"#you could also vary delta. What happens if delta = 1?\n",
"DELTA = 5\n",
"\n",
"N_DESIGNS = 10\n",
"N_0 = 5\n",
"\n",
"ocba = OCBA(model=model, \n",
" n_designs=N_DESIGNS, \n",
" budget=BUDGET, \n",
" delta=DELTA, \n",
" n_0=N_0, \n",
" obj='max')\n",
"\n",
"#create an experiment and repeat 1000 times.\n",
"exp = Experiment(env=model, procedure=ocba, best_index=9, objective='max', \n",
" replications=1000)\n",
"\n",
"results = exp.execute()\n",
"\n",
"#print the probability of correct selection achieved.\n",
"print(f'PCS={results.p_correct_selections}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise 4. Calculating the Expected Opportunity Cost.\n",
"\n",
"Another way to evaluate a procedures performance is to calculate the Expected Opportunity Cost. This is the average difference between the best mean and the selected mean in each repeat of the experiment.\n",
"\n",
"**Task: try varying the budget, delta and n_0. What is the effect?**\n",
"\n",
"We will only run a small number of repeats of the experiments to save runtime!"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"PCS=0.7\n",
"selected designs: [5 5 5 8 5 5 5 9 2 5]\n",
"EOC=0.029999999999999895\n"
]
}
],
"source": [
"# KN++ needed a large budget to reach its desired precision. \n",
"# We will try a smaller budgets here\n",
"BUDGET = 20000\n",
"DELTA = 5\n",
"N_DESIGNS = 10\n",
"N_0 = 5\n",
"\n",
"means = [4, 4.1, 4.2, 4, 4.1, 4.3, 4, 4.1, 4.2, 4.2]\n",
"variances = [1, 1, 1, 0.1, 0.1, 10, 10, 10, 10, 0.1]\n",
"\n",
"guass_model = custom_gaussian_model(means, variances)\n",
"\n",
"ocba = OCBA(model=guass_model, \n",
" n_designs=N_DESIGNS, \n",
" budget=BUDGET, \n",
" delta=5, \n",
" n_0=50, \n",
" obj='max')\n",
"\n",
"# create an experiment and repeat 10 times (in practice we would use more!)\n",
"exp = Experiment(env=guass_model, procedure=ocba, best_index=5, objective='max', \n",
" replications=10)\n",
"\n",
"results = exp.execute()\n",
"\n",
"# the probability of correct selection achieved.\n",
"print(f'PCS={results.p_correct_selections}')\n",
"\n",
"# IF YOU REMEMBER FROM Example 2 THE DESIGNS IN THE IZ WERE [2 5 8 9]\n",
"# Let's print out the selections made...\n",
"print(f'selected designs: {results.selections}')\n",
"\n",
"# the average of the best_mean - selected_design_mean for each rep of the exp.\n",
"print(f'EOC={results.expected_opportunity_cost}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Optimal Computing Budget Allocation Top M (OCBA-m)\n",
"\n",
"In many situations it is not possible or desirable to model all of the detail in a simulation model. A procedure that offers a top performing group of designs is more useful in these circumstances. It allows a decision maker to choose between good designs while taking other factors into account.\n",
"\n",
"So far we have only considered selecting the best system. OCBA-m extended OCBA to identify the **top m designs**. \n",
"\n",
"A top-m solution OCBA object can be created by importing `ovs.fixed_budget.OCBAM`\n",
"\n",
"An object of type `OCBAM` takes the following parameters:\n",
"\n",
"* **model** - a simulation model\n",
"* **n_designs** - int, the number of competing designs to compare\n",
"* **budget** - int, the total number of replications to allocate across designs\n",
"* **delta** - int, the incremental amount of replications to allocate at each round\n",
"* **n_0** - int, $n_0$ the number of initial replications (default=5)\n",
"* **m** - int, $m$ the number of top designs to return. m must be >=2. For m = 1 see OCBA.\n",
"* **obj** - str, 'min' if minimisation; 'max' if maximisation (default='min')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 4. Selecting the top 2 designs."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"N_DESIGNS = 10\n",
"BUDGET = 1000\n",
"DELTA = 5\n",
"N_0 = 5\n",
"\n",
"#you can vary this parameter to change the number selected.\n",
"M = 3\n",
"\n",
"guass_model = gaussian_sequence_model(1, N_DESIGNS)\n",
"\n",
"ocbam = OCBAM(model=guass_model, \n",
" n_designs=N_DESIGNS, \n",
" budget=BUDGET, \n",
" delta=DELTA, \n",
" n_0=N_0, \n",
" m=M,\n",
" obj='max')"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"OCBA(n_designs=10, m=3, budget=1000, delta=5, n_0=5, obj=max)\n"
]
}
],
"source": [
"print(ocbam)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"best designs\t[7 8 9]\n",
"allocations\t[ 20 23 21 38 39 75 79 303 302 100]\n",
"total reps\t1000\n"
]
}
],
"source": [
"#print out the results\n",
"best_designs = ocbam.solve()\n",
"print(f'best designs\\t{best_designs}')\n",
"print(f'allocations\\t{ocbam._allocations}')\n",
"print(f'total reps\\t{ocbam._allocations.sum()}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise 5. Select the top 3 designs out of 100 designs.\n",
"\n",
"Task: Use OCBA-m to select the **smallest** 3 designs from the example problem. What impact does varying the budget, delta and $n_0$ have on performance?"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"N_DESIGNS = 100\n",
"BUDGET = 1000\n",
"DELTA = 1\n",
"N_0 = 5\n",
"\n",
"#you can vary this parameter to change the number selected.\n",
"M = 3\n",
"\n",
"guass_model = gaussian_sequence_model(1, N_DESIGNS)\n",
"\n",
"\n",
"#insert your code here\n",
"#ocbam = ?"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"best designs\t[7 8 9]\n",
"allocations\t[ 20 23 21 38 39 75 79 303 302 100]\n",
"total reps\t1000\n"
]
}
],
"source": [
"#print out the results\n",
"best_designs = ocbam.solve()\n",
"print(f'best designs\\t{best_designs}')\n",
"print(f'allocations\\t{ocbam._allocations}')\n",
"print(f'total reps\\t{ocbam._allocations.sum()}')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# THE END\n",
"\n",
"What did you think of the procedures?"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}