Optimisation with SimPy: Finding the Best Configuration

import simpy
import random
import numpy as np

def run_simulation(num_servers, arrival_rate, service_rate, duration=1000):
    """Run simulation and return mean wait time."""
    random.seed(42)
    env = simpy.Environment()
    server = simpy.Resource(env, capacity=num_servers)
    wait_times = []

    def customer():
        arrival = env.now
        with server.request() as req:
            yield req
            wait_times.append(env.now - arrival)
            yield env.timeout(random.expovariate(service_rate))

    def arrivals():
        while True:
            yield env.timeout(random.expovariate(arrival_rate))
            env.process(customer())

    env.process(arrivals())
    env.run(until=duration)

    return np.mean(wait_times) if wait_times else float('inf')

# Grid search
results = []
for servers in range(1, 6):
    for arrival_rate in [0.5, 1.0, 1.5, 2.0]:
        wait = run_simulation(servers, arrival_rate, service_rate=0.5)
        results.append({
            'servers': servers,
            'arrival_rate': arrival_rate,
            'mean_wait': wait
        })

# Find best
best = min(results, key=lambda x: x['mean_wait'])
print(f"Best config: {best}")

def total_cost(num_servers, arrival_rate, service_rate):
    """Calculate total cost including wait time cost."""
    SERVER_COST = 100  # Cost per server per hour
    WAIT_COST = 50     # Cost per customer per minute wait

    # Run simulation
    mean_wait = run_simulation(num_servers, arrival_rate, service_rate)
    customers_per_hour = arrival_rate * 60

    # Calculate costs
    server_cost = num_servers * SERVER_COST
    wait_cost = customers_per_hour * mean_wait * WAIT_COST

    return server_cost + wait_cost

# Search for minimum cost
for servers in range(1, 10):
    cost = total_cost(servers, arrival_rate=2, service_rate=0.5)
    print(f"Servers: {servers}, Cost: {cost:.2f}")

from scipy.optimize import minimize_scalar, minimize

def objective(num_servers):
    """Objective function to minimize."""
    num_servers = int(round(num_servers))
    if num_servers < 1:
        return float('inf')

    mean_wait = run_simulation(num_servers, arrival_rate=2, service_rate=0.5)
    server_cost = num_servers * 100
    wait_cost = mean_wait * 50 * 120  # 2 arrivals/min * 60 min

    return server_cost + wait_cost

# Optimize
result = minimize_scalar(objective, bounds=(1, 20), method='bounded')
optimal_servers = int(round(result.x))
print(f"Optimal servers: {optimal_servers}")
print(f"Minimum cost: {result.fun:.2f}")

def multi_objective(params):
    """Optimize multiple parameters."""
    num_servers = int(round(params[0]))
    buffer_size = int(round(params[1]))

    if num_servers < 1 or buffer_size < 1:
        return float('inf')

    # Run simulation with these parameters
    result = run_simulation_complex(num_servers, buffer_size)

    # Weighted objective
    return (result['wait'] * 50 +
            result['abandonment_rate'] * 1000 +
            num_servers * 100)

# Optimize
result = minimize(
    multi_objective,
    x0=[3, 10],  # Initial guess
    method='Nelder-Mead',
    bounds=[(1, 10), (5, 50)]
)

print(f"Optimal: servers={result.x[0]:.0f}, buffer={result.x[1]:.0f}")

def meets_target(num_servers, target_wait=5, target_percentile=0.9):
    """Check if configuration meets service target."""
    # Run multiple replications for confidence
    waits = []
    for seed in range(10):
        wait = run_simulation(num_servers, arrival_rate=2, service_rate=0.5, seed=seed)
        waits.append(wait)

    percentile_wait = np.percentile(waits, target_percentile * 100)
    return percentile_wait <= target_wait

# Binary search for minimum servers
low, high = 1, 20
while low < high:
    mid = (low + high) // 2
    if meets_target(mid):
        high = mid
    else:
        low = mid + 1

print(f"Minimum servers to meet target: {low}")

# Note: Requires deap library
from deap import base, creator, tools, algorithms
import random

# Define fitness (minimize wait, minimize cost)
creator.create("FitnessMin", base.Fitness, weights=(-1.0,))
creator.create("Individual", list, fitness=creator.FitnessMin)

toolbox = base.Toolbox()

# Genes: [servers, buffer_size]
toolbox.register("servers", random.randint, 1, 10)
toolbox.register("buffer", random.randint, 5, 50)
toolbox.register("individual", tools.initCycle, creator.Individual,
                 (toolbox.servers, toolbox.buffer), n=1)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)

def evaluate(individual):
    servers, buffer = individual[0], individual[1]
    cost = run_simulation_with_cost(servers, buffer)
    return (cost,)

toolbox.register("evaluate", evaluate)
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutUniformInt, low=[1, 5], up=[10, 50], indpb=0.2)
toolbox.register("select", tools.selTournament, tournsize=3)

# Run GA
population = toolbox.population(n=50)
result, log = algorithms.eaSimple(population, toolbox, cxpb=0.5, mutpb=0.2,
                                   ngen=20, verbose=True)

best = tools.selBest(result, k=1)[0]
print(f"Best solution: servers={best[0]}, buffer={best[1]}")

class SimulationOptimizer:
    def __init__(self, simulation_func, param_bounds):
        self.sim_func = simulation_func
        self.bounds = param_bounds
        self.history = []

    def objective(self, params):
        """Run simulation and return objective value."""
        result = self.sim_func(**dict(zip(self.bounds.keys(), params)))
        self.history.append({'params': params, 'result': result})
        return result

    def optimize_grid(self, grid_points=5):
        """Grid search optimization."""
        import itertools

        grids = []
        for name, (low, high) in self.bounds.items():
            grids.append(np.linspace(low, high, grid_points))

        best_params = None
        best_value = float('inf')

        for combo in itertools.product(*grids):
            value = self.objective(combo)
            if value < best_value:
                best_value = value
                best_params = combo

        return dict(zip(self.bounds.keys(), best_params)), best_value

    def optimize_scipy(self, x0=None):
        """scipy.optimize-based search."""
        if x0 is None:
            x0 = [(b[0] + b[1]) / 2 for b in self.bounds.values()]

        result = minimize(
            self.objective,
            x0,
            method='Nelder-Mead'
        )

        return dict(zip(self.bounds.keys(), result.x)), result.fun

# Usage
optimizer = SimulationOptimizer(
    simulation_func=run_simulation_with_params,
    param_bounds={
        'servers': (1, 10),
        'buffer': (5, 50),
        'threshold': (0.5, 0.9)
    }
)

best_params, best_value = optimizer.optimize_scipy()
print(f"Optimal: {best_params}, Value: {best_value}")

def robust_objective(params, n_replications=10):
    """Average over multiple replications."""
    results = []
    for seed in range(n_replications):
        result = run_simulation(*params, seed=seed)
        results.append(result)

    return np.mean(results)  # Or use worst-case, percentile, etc.

def percentile_objective(params, n_replications=10, percentile=90):
    """Optimize for percentile (robust to variability)."""
    results = []
    for seed in range(n_replications):
        result = run_simulation(*params, seed=seed)
        results.append(result)

    return np.percentile(results, percentile)