From 0df29235231444595bf08c49830880e715bf815b Mon Sep 17 00:00:00 2001 From: chloetsetsekas Date: Mon, 11 May 2026 17:50:14 -0400 Subject: [PATCH 1/3] Add Call Center baby staffing simulation model --- docs/source/models/ccbaby.rts | 194 +++++++++++ docs/source/models/tollnw.rst | 164 ++++++++++ notebooks/demo_problem.py | 4 +- simopt/models/ccbaby.py | 589 ++++++++++++++++++++++++++++++++++ 4 files changed, 949 insertions(+), 2 deletions(-) create mode 100644 docs/source/models/ccbaby.rts create mode 100644 docs/source/models/tollnw.rst create mode 100644 simopt/models/ccbaby.py diff --git a/docs/source/models/ccbaby.rts b/docs/source/models/ccbaby.rts new file mode 100644 index 000000000..24f283f1c --- /dev/null +++ b/docs/source/models/ccbaby.rts @@ -0,0 +1,194 @@ +Call Center Staffing +==================== + +See the :mod:`simopt.models.callcenter` module for API details. + +Model: Call Center Queueing and Staffing (CALLCENTER) +----------------------------------------------------- + +Description +^^^^^^^^^^^ + +This model simulates a call center over a 16-hour operating day. Calls arrive +according to a non-homogeneous Poisson process with rate function +:math:`\lambda(t)` for :math:`0 \leq t \leq 16`, where time is measured in +hours. Each incoming call enters the system and waits in queue if there are no +available agents. Then, it either begins service with an available agent or abandons +after its patience time expires. The system has :math:`T` trunk lines, so at most +:math:`T` calls can be in the system at once, including both calls in service and +calls waiting on hold. Calls that arrive when the system is full receive a busy +signal and do not call back. + +Call handling times or service times are gamma distributed with mean :math:`\mu_h` +minutes and variance :math:`\sigma_h^2` minutes^2. Patience times are also gamma +distributed with mean :math:`\mu_p` minutes and variance :math:`\sigma_p^2` minutes^2. +The arrival process, handle times, and patience times are mutually independent. + +Agents work 8-hour shifts structured as follows: they work for :math:`x` +hours, take a 30-minute lunch break, and then work for :math:`8-x` hours, +where :math:`x \in \{3, 3.5, 4, 4.5\}`. Agents can begin work every half hour, +starting from :math:`t=0` through :math:`t=8`. Thus the allowed shift start +times are :math:`0, 0.5, 1, \ldots, 8`. Agents starting at :math:`t=8` finish at +:math:`t=16.5`. At the end of a shift, agents finish any call they are currently +handling and then leave, except if an agent is working the last shift that ends at +:math: `t=16.5`. In that case, at the end of the day, the call center stops receiving +calls at :math:`t=16`, but any calls still in the system, active or queued remain +until they are served or abandon before agents leave, even if this means the shift +extends into overtime, ie. beyond :math: `t=16.5`. Agents are paid a flat rate of +$:math:`r` per hour. This holds even when agents in the last shift work past +:math: `t=16.5`. :math: `c_j` = :math:`8r` for all :math:`J` in :math:`j`. is the +cost of assigning one agent to start at time :math: `j/2`. + +Performance is measured separately in each hour :math:`i=1,\ldots,16`. +Let :math:`N_i` be the number of calls received in hour :math:`i`, including +blocked calls. Let :math:`S_i` be the number of calls answered within 20 +seconds in hour :math:`i`. Calls that abandon or are blocked do not count in +:math:`S_i`. The service level requirement for each hour is + +:math:`E[S_i] \geq 0.8E[N_i], \quad i=1,\ldots,16.` + +The expected number of arrivals :math:`E[N_i]` can be computed analytically, +but :math:`E[S_i]` must be estimated via simulation. + +Sources of Randomness +^^^^^^^^^^^^^^^^^^^^^ + +There are 3 main sources of randomness in this model: + +1. The call arrival process, modeled as a non-homogeneous Poisson process + with rate :math:`\lambda(t)` over :math:`0 \leq t \leq 16`. +2. The call handling times, modeled as gamma random variables with mean + :math:`\mu_h` and variance :math:`\sigma_h^2`. +3. The customer patience times, modeled as gamma random variables with mean + :math:`\mu_p` and variance :math:`\sigma_p^2`. + +The recommended arrival-rate function is + +:math:`\lambda(t) = 500 + 500\sin\left(\frac{3\pi t - 16\pi}{32}\right), \quad 0 \leq t \leq 16.` + +Model Factors +^^^^^^^^^^^^^ + +* arrival_rate_function: + * Default: :math:`500 + 500\sin\left(\frac{3\pi t - 16\pi}{32}\right)` +* time_open: + * Default: 16 +* mean_handle_time: + * Default: 6 +* var_handle_time: + * Default: 2 +* mean_patience_time: + * Default: 2 +* var_patience_time: + * Default: 1 +* trunk_lines: + * Default: 150 +* hourly_wage: + * Default: 18 +* service_level_target: + * Default: 0.8 +* service_level_window_seconds: + * Default: 20 +* constraint_tolerance: + * Default: 0.5 +* allowed_shift_starts: + * Default: [0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8] +* shift_break_options: + * Default: [3, 3.5, 4, 4.5] + +Responses +^^^^^^^^^ + +* total_cost: Total wage cost +* arrivals_by_hour: Number of arrivals per hour +* answered_within_20_seconds_by_hour: Calls meeting service target per hour +* blocked_calls: Calls rejected due to capacity +* abandoned_calls: Calls that leave queue +* service_level_by_hour: :math:`S_i / N_i` +* constraint_lhs_by_hour: :math:`E[S_i] - 0.8E[N_i]` + +References +^^^^^^^^^^ + +Pasupathy, R., & Henderson, S. G. (2007). Call Center Example. + +Optimization Problem: Minimize Staffing Cost Subject to Service Levels (CALLCENTER-1) +-------------------------------------------------------------------------------------- + +Decision Variables +^^^^^^^^^^^^^^^^^^ + +* :math:`x_{j,k}`, number of agents starting at time :math:`j/2`, + and take lunch :math: `h_k` hours after their start time, where + + * :math:`j=0,\ldots,16` + * :math:`h_k \in \{3, 3.5, 4, 4.5\}` + +The staffing decision is a :math:`\times 4` integer matrix with the row determining the +shift start time and the column determinging the lunch start time. + +Objective +^^^^^^^^^^ + +Minimize total labor cost: + +:math:`\sum_{j=0}^{16} c_j x_{j,k}` + + +:math:`c_j` is calculated in the simulation as follows: + TotalCost = :math:`r * np.sum(Agents[:, 3] - Agents[:, 0])` + + +Constraints +^^^^^^^^^^^ + +* For each hour :math:`i=1,\ldots,16`: + + :math:`E[S_i] - 0.8E[N_i] \geq 0` + +* With tolerance: + + :math:`E[S_i] - 0.8E[N_i] \geq -0.5` + +* :math:`x_j \geq 0` and integer + +Problem Factors +^^^^^^^^^^^^^^^ + +* Budget: + * Default: 1000 + * Suggested: 10000, 100000 + +Fixed Model Factors +^^^^^^^^^^^^^^^^^^^ + +* arrival_rate_function +* time_open = 16 +* mean_handle_time = 6 +* var_handle_time = 2 +* mean_patience_time = 2 +* var_patience_time = 1 +* trunk_lines = 150 +* hourly_wage = 18 +* service_level_target = 0.8 +* constraint_tolerance = 0.5 + +Starting Solution +^^^^^^^^^^^^^^^^^ + +* :math:`x_i = 8` for all :math:`i` + +* :math:`x_j` is the number of agents starting their shift at a certain time, + not the number of agents actually working at that time. + +Each :math:`x_j \sim \text{Uniform}\{0,\ldots,10\}` + +Optimal Solution +^^^^^^^^^^^^^^^^ + +unknown + +Optimal Objective Function Value +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +unknown \ No newline at end of file diff --git a/docs/source/models/tollnw.rst b/docs/source/models/tollnw.rst new file mode 100644 index 000000000..938d8169b --- /dev/null +++ b/docs/source/models/tollnw.rst @@ -0,0 +1,164 @@ +Toll Networks +======================= +See the :mod:'simopt.models.tollnw' module for API details. + +Model: Toll Road Improvements in a Newtork +---------------------------------------- + +Description +^^^^^^^^^^^ + +This model represents a connected directed graph :math:`G= (V,E)` describing a network +of toll roads. Customers wishing to travel from point :math:`i` to point :math:`j` +arrive according to a Poisson process with rate :math:`lambda_{ij}`. + +Each directed edge represents a toll road with price :math:`p_{ij}`. Travel times on each +road are triangularly distributed with route-specific parameters. The network operates as +an open Jackson network with interconnected queues. + +Each road has the capacity for one vehicle at a time. If a vehicle is currently traveling +on a road, arriving vehicles form a queue and wait until the road becomes available. + +The operator can choose to invest in road improvements. If :math:`x_{ij}` hundred dollars +is spent improving road :math:`(i,j)`, the resulting mode of the triangular travel-time +distribution becomes: + +:math:`a_{ij} + (b_{ij} - a_{ij}) \exp(-x_{ij})` + +where: +- :math:`a_{ij}` is the minimum travel time +- :math:`b_{ij}` is the maximum travel time + + +If no investment is made (i.e. :math:`x_{ij} = 0`), the mode equals :math:`b_{ij}`. + +The operator seeks to maximize expected profit over a finite horizon :math:`T`: + +:math:`\mathbb{E}[\text{Profit}] = \sum_{(i,j) \in E} p_{ij} \mathbb{E}[N_{ij}(T)] - \sum_{(i,j) \in E} x_{ij}` + +where: +:math:`N_{ij}(T)` is the number of trips completed on road :math:`(i,j)` before time :math:`T`. + + + +Sources of Randomness +^^^^^^^^^^^^^^^^^^^^^ + +There are three source of randomness: + +1. External arrivals modeled as Poisson processes with rates :math:`lambda_{ij}` +2. Internal routing decisions governed by routing matrix :math:`R` +3. Travel times following triangular distribution + + +Model Factors +^^^^^^^^^^^^^ + +- adjacency_matrix: graph structure describing road connectivity +- lambda_ij: external Poisson arrival rates + - default: :math:`\lambda_{ij} = i + j` for :math:`i \neg j`, :math:`i, j \in \{1,2,3\}` +- routing_matrix: internal routing matrix with exit probabilities +- pij: toll/amount charged per completed trips + - default: 1 for all roads +- aij: minimum travel time of triangular distribution + - default: 0 +- bij: maximum travel time of triangular distribution + - default: 1 +- investment_levels: investment decisions :math:`x_{ij}` +- T: length of operating horizon + - default: 96 (12 hours) +- time_unit: one unit equals 10 minutes + + +Recommended Parameter Settings +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Let :math:`G` be a complete graph with three vertices: + +..math:: + + A = + \begin{bmatrix} + 0 & 1 & 1 \\ + 1 & 0 & 1 \\ + 1 & 1 & 0 + \end{bmatrix} + +Routing Matrix: + +..math:: + + R = + \begin{bmatrix} + 0 & 1/3 & 1/3 & 1/3 + 1/2 & 0 & 1/4 & 1/4 + 1/5 & 1/5 & 0 & 3/5 + \end{bmatrix} + +The last column represents exiting the network. + +Responses +^^^^^^^^^ + +- total_profit: expected profit over horizon :math:`T` +- total_trips_completed: total completed trips across all roads +- road_trip_counts: vector of :math:`N_{ij}(T)` values +- average_queue_lengths: time-averaged queue lengths per road + +References +^^^^^^^^^^ + +Eckman, D. (2017). Toll Road Improvements in a Network. + +Optimization Problem: Maximize Expected Profit (TOLLNETWORK-1) +-------------------------------------------------------------- + +Decision Variables +^^^^^^^^^^^^^^^^^^ + +- investment_levels: :math:`x_{ij}` (continuous, :math:`\ge 0`) + +Objective +^^^^^^^^^^ + +Maximize expected profit over horizon :math:`T`. + +Constraints +^^^^^^^^^^^ + +- :math:`x_{ij} \ge 0` for all roads +- Continuous decision Variables xij>=0 + +Problem Factors +^^^^^^^^^^^^^^^ + +- Budget: maximum number of replications for solver + - default: 1000 + +Fixed Model Factors +^^^^^^^^^^^^^^^^^^^ + +- aij = 0 +- bij = 1 +- pij = 1 +- T = 96 + +Starting Solution +^^^^^^^^^^^^^^^^^ + +- :math:`x_{ij} = 1` for all roads + +Random Solutions +^^^^^^^^^^^^^^^^ + +Generate continuous nonnegative investment vectors over all roads + +Optimal Solution +^^^^^^^^^^^^^^^^ + +Unknown + +Optimal Objective Function Value +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +Unknown \ No newline at end of file diff --git a/notebooks/demo_problem.py b/notebooks/demo_problem.py index f33615332..93ceaa47f 100644 --- a/notebooks/demo_problem.py +++ b/notebooks/demo_problem.py @@ -62,7 +62,7 @@ fixed_factors = {"initial_solution": (2,), "budget": 500} myproblem = CntNVMaxProfit(fixed_factors=fixed_factors) -x = (3,) +x = (1,) # ----------------------------------------------- @@ -76,7 +76,7 @@ # mysolution = Solution(x, myproblem) # ----------------------------------------------- -num_macroreps = 10 +num_macroreps = 1500 # %% # Create and attach rngs to solution diff --git a/simopt/models/ccbaby.py b/simopt/models/ccbaby.py new file mode 100644 index 000000000..2fe9f5a35 --- /dev/null +++ b/simopt/models/ccbaby.py @@ -0,0 +1,589 @@ +#imports +import math +import numpy as np +import matplotlib.pyplot as plt + +from simopt.base import ( + ConstraintType, + Model, + Objective, + Problem, + RepResult, + VariableType, +) + +#model simulation +def ccbaby(x, runlength, seed, other=None): + """ + Python translation of the MATLAB function CCBaby. + + Parameters + ---------- + x : array-like + Vector of the number of agents starting shifts at times j/2. + runlength : int + Number of days to simulate. + seed : int + Positive integer random seed. + other : unused + Included only to match original signature. + + Returns + ------- + fn : float + Mean total cost across simulated days. + FnVar : float + Variance of total cost divided by number of days. + FnGrad : float + NaN (not used, matches MATLAB code). + FnGradCov : float + NaN (not used, matches MATLAB code). + constraint : np.ndarray + Mean hourly performance measure across days. + ConstraintCov : np.ndarray + Covariance matrix of hourly performance measures. + ConstraintGrad : float + NaN (not used, matches MATLAB code). + ConstraintGradCov : float + NaN (not used, matches MATLAB code). + """ + + FnGrad = np.nan + FnGradCov = np.nan + constraint = np.nan + ConstraintCov = np.nan + ConstraintGrad = np.nan + ConstraintGradCov = np.nan + + #C sets maximum number of agents allowed to start at any time slot + #x is your staffing decision vector, here it is converted into a NumPy array of integers + C = 100 + x = np.asarray(x, dtype=int) + hk = np.array([3.0, 3.5, 4.0, 4.5]) + + #if anything about the inputs is invalid, stop the simulation and return NaN + #checks if staffing value is negative, if any entry exceeds the max number of agents allowed to start, + #if the simulation runs for at least one day, random seed must be positive and integer + if ( + x.ndim != 2 + or x.shape[1] != 4 + or np.any(x < 0) + or np.any(x > C) + or runlength <= 0 + or seed <= 0 + or int(seed) != seed + ): + print("x must be a 2D nonneg integer array with 4 columns: x[j,k]") + fn = np.nan + FnVar = np.nan + return ( + fn, + FnVar, + FnGrad, + FnGradCov, + constraint, + ConstraintCov, + ConstraintGrad, + ConstraintGradCov, + ) + + # Parameters + lambdaMax = 1000 # Arrival upper bound + serviceShape = 18 # Gamma shape for service times + serviceScale = 1 / 3 / 60 # Gamma scale for service times (hours) + patienceShape = 4 # Gamma shape for patience times + patienceScale = 1 / 2 / 60 # Gamma scale for patience times (hours) + TrunkMax = 150 # Number of trunk lines + Salary = 18 # Agent hourly wage + nDays = int(runlength) + tMax = 16.5 # Length of a day + m = 16 # Number of constraints (not otherwise used) + + #converts staffing plan into a working variable, then sets variable for total agents in system + xjk = x + nAgents = int(np.sum(xjk)) + + print("x matrix:") + print(xjk) + print("x shape:", xjk.shape) + print("total agents =", nAgents) + print("agents by start time =", np.sum(xjk, axis=1)) + print("agents by lunch choice =", np.sum(xjk, axis=0)) + + #storage arrays to store (# of hours in a day, # of days), and a 1D array of length nDays + PerformanceMeasure = np.zeros((math.ceil(tMax), nDays)) + TotalCost = np.zeros(nDays) + + #create 5 separate random number generators to mimic separate MATLAB streams + seed_seq = np.random.SeedSequence(seed) + child_seeds = seed_seq.spawn(4) + + dummy_rng = np.random.default_rng(child_seeds[0]) + service_rng = np.random.default_rng(child_seeds[1]) + patience_rng = np.random.default_rng(child_seeds[2]) + ar_rng = np.random.default_rng(child_seeds[3]) + + #generates exponentially distributed time gapes between customer call arrivals + Dummy = dummy_rng.exponential(scale=1 / lambdaMax, size=(nDays, 100000)) + + #pre-generates random values for service times, customer patience, and arrival acceptance + sTime = service_rng.gamma(shape=serviceShape, scale=serviceScale, size=(nDays, 100000)) + pTime = patience_rng.gamma(shape=patienceShape, scale=patienceScale, size=(nDays, 100000)) + + AR = ar_rng.random(size=(nDays, 100000)) + + #counters for each day + D = np.ones(nDays, dtype=int) + S = np.ones(nDays, dtype=int) + P = np.ones(nDays, dtype=int) + A = np.ones(nDays, dtype=int) + + #convert MATLAB 1-based counters to Python 0-based by subtracting 1 when indexing + #simulate one full day of call center operations, then repeat for the next day + for k in range(nDays): + served_immediate = 0 + served_wait = 0 + abandoned = 0 + blocked = 0 + #stores current state of every agent for that day, including shift start time, + #lunch start time, next time agent is available, departure/off-duty time, + #whether lunch has happened already + Agents = np.zeros((nAgents, 5)) + #stores call information, including call arrival time, time service starts, + #time call leaves the system + Calls = np.zeros((3, 1)) + + #create start times and break times for agents + #converts decision variable into actual agents + #for each start time, and for each lunch option, create the specificed number + #of agents + num = 0 + max_j = int((tMax - 8.5) * 2) + + for j in range(max_j + 1): + for k_idx in range(4): + num_agents_starting = xjk[j, k_idx] if j < xjk.shape[0] else 0 + for _ in range(num_agents_starting): + #each agent gets assigned these 3 things + start_time = j/2 + lunch_time = start_time + hk[k_idx] + + Agents[num, 0] = start_time + Agents[num, 1] = lunch_time + Agents[num, 2] = start_time + num += 1 + + if k == 0: + print("\nFirst 10 agents (start time, lunch time):") + print(Agents[:10, 0:2]) + + #first arrival via acceptance-rejection + #finds first customer of the day + nextCall = 0.0 + dummy = 0.0 #running timeof candidate arrivals + #generates candidate arrivals and accepts only some + while nextCall == 0: + dummy += Dummy[k, D[k] - 1] + D[k] += 1 + #time-varying arrival rate at that moment so that call intensity depends on time of day + funValue = 500 + 500 * math.sin((3 * math.pi * dummy - 16 * math.pi) / 32) + a = AR[k, A[k] - 1] + A[k] += 1 + #random acceptance-rejection number + if a <= funValue / 1000: + nextCall = dummy + nCalls = 1 + Calls[0, 0] = nextCall + + trunks = np.zeros(TrunkMax) + + #as long as the next call arrives before the end of the day, process it + while nextCall < tMax: + #if there isn't enough room in Calls var, add another column + if Calls.shape[1] < nCalls: + Calls = np.hstack([Calls, np.zeros((3, 1))]) + + #stores call's arrival time + Calls[0, nCalls - 1] = nextCall + minTime = 25.0 + minAgent = -1 + callTaken = False + + #check trunk availability + minTrunkIndex = np.argmin(trunks) + minTrunk = trunks[minTrunkIndex] + + #a blocked call immediately leaves the system at arrival time, this is consistent with "busy signal" + if minTrunk >= nextCall: + Calls[2, nCalls - 1] = nextCall + callTaken = True + blocked += 1 + + #assigns random charactertistics to current caller, including + #how long service will take if answered, and caller's patience time + serviceTime = sTime[k, S[k] - 1] + patienceTime = pTime[k, P[k] - 1] + S[k] += 1 + P[k] += 1 + + #checks agent statuses to decide who can take the call + for i in range(nAgents): + #update agent's status if available by nextCall and not already left + if Agents[i, 2] <= nextCall and Agents[i, 3] == 0 and not callTaken: + #check if Agent left for lunch break + if Agents[i, 1] <= nextCall and Agents[i, 4] == 0: + Agents[i, 2] = nextCall + 0.5 + Agents[i, 4] = 1 + + #check if Agent left job already/already off shift + if Agents[i, 0] + 8.5 <= nextCall and Agents[i, 0] + 8.5 < 16.5: + Agents[i, 3] = Agents[i, 0] + 8.5 + + #if agent available after those two checks, they take the call + if Agents[i, 2] <= nextCall and Agents[i, 3] == 0 and not callTaken: + Calls[1, nCalls - 1] = nextCall + #update agent's next available time, mark and record taken call + Agents[i, 2] = nextCall + serviceTime + callTaken = True + served_immediate += 1 + Calls[2, nCalls - 1] = nextCall + serviceTime + + #check whether lunch happens during the call + #if so, then agent starts lunch after call is over + if Agents[i, 2] >= Agents[i, 1] and Agents[i, 4] == 0: + Agents[i, 2] += 0.5 + Agents[i, 4] = 1 + + #check whether shift end happens during call + if Agents[i, 2] >= Agents[i, 0] + 8.5 and Agents[i, 0] + 8.5 < 16.5: + Agents[i, 3] = Agents[i, 2] + + break + + #track next available agent + if Agents[i, 2] <= minTime and Agents[i, 3] == 0: + minTime = Agents[i, 2] + minAgent = i + + #no immediately available agents + if minTime < nextCall + patienceTime and not callTaken and minAgent != -1: + Calls[1, nCalls - 1] = Agents[minAgent, 2] + Agents[minAgent, 2] = Agents[minAgent, 2] + serviceTime + Calls[2, nCalls - 1] = Agents[minAgent, 2] + served_wait += 1 + + if Agents[minAgent, 2] >= Agents[minAgent, 1] and Agents[minAgent, 4] == 0: + Agents[minAgent, 2] += 0.5 + Agents[minAgent, 4] = 1 + + if ( + Agents[minAgent, 2] >= Agents[minAgent, 0] + 8.5 + and Agents[minAgent, 0] + 8.5 < 16.5 + ): + Agents[minAgent, 3] = Agents[minAgent, 2] + elif not callTaken: + Calls[2, nCalls - 1] = nextCall + patienceTime + abandoned += 1 + + #generate next arrival/call + #advances simulation clock to the next real customer arrival + #so the loop can process the next calll + success = False + while not success: + dummy += Dummy[k, D[k] - 1] + D[k] += 1 + funValue = 500 + 500 * math.sin((3 * math.pi * dummy - 16 * math.pi) / 32) + a = AR[k, A[k] - 1] + A[k] += 1 + if a <= funValue / 1000: + nextCall = dummy + success = True + + trunks[minTrunkIndex] = Calls[2, nCalls - 1] + nCalls += 1 + + #update departure times for last agents + #some agents may not yet have a recorded departure time- this fills that in + for i in range(nAgents): + if Agents[i, 3] == 0: + if Agents[i, 2] > 16.5: + Agents[i, 3] = Agents[i, 2] + else: + Agents[i, 3] = Agents[i, 0] + 8.5 + + #compute days outputs + callsPerHour = np.zeros(math.ceil(tMax)) + under20sPerHour = np.zeros(math.ceil(tMax)) + + for i in range(nCalls - 1): + hour = int(math.floor(Calls[0, i])) + callsPerHour[hour] += 1 + wait_time = Calls[1, i] - Calls[0, i] + if 0 <= wait_time <= 0.0056: + under20sPerHour[hour] += 1 + print("callsPerHour:", callsPerHour.astype(int)) + print("under20sPerHour:", under20sPerHour.astype(int)) + print("under20 <= calls each hour:", np.all(under20sPerHour <= callsPerHour)) + + print(f"\n--- Day {k+1} ---") + print("total calls:", nCalls - 1) + print("served immediately:", served_immediate) + print("served after waiting:", served_wait) + print("abandoned:", abandoned) + print("blocked:", blocked) + print("sum check:", served_immediate + served_wait + abandoned + blocked) + #cost calculation + TotalCost[k] = Salary * np.sum(Agents[:, 3] - Agents[:, 0]) + PerformanceMeasure[:, k] = 0.8 * callsPerHour - under20sPerHour + print("daily cost:", TotalCost[k]) + + fn = np.mean(TotalCost) + FnVar = np.var(TotalCost, ddof=1) / nDays if nDays > 1 else 0.0 + + Performance = np.zeros((math.ceil(tMax), 2)) + Performance[:, 0] = np.mean(PerformanceMeasure, axis=1) + if nDays > 1: + Performance[:, 1] = 1.96 * np.std(PerformanceMeasure, axis=1, ddof=1) / np.sqrt(nDays) + else: + Performance[:, 1] = 0.0 + + constraint = Performance[:, 0] + ConstraintCov = np.cov(PerformanceMeasure) + + return ( + fn, + FnVar, + FnGrad, + FnGradCov, + constraint, + ConstraintCov, + ConstraintGrad, + ConstraintGradCov, + ) + +#model class +class CallCenter(Model): + """SimOpt wrapper for the ccbaby call center simulation.""" + + class_name_abbr = "CALLCENTER" + class_name = "Call Center Staffing" + n_rngs = 1 + n_responses = 2 + + def __init__(self, fixed_factors=None): + super().__init__(fixed_factors) + + def replicate(self): + x = self.factors["staffing_matrix"] + runlength = self.factors.get("runlength", 1) + seed = self.factors.get("seed", 123) + + ( + fn, + FnVar, + FnGrad, + FnGradCov, + constraint, + ConstraintCov, + ConstraintGrad, + ConstraintGradCov, + ) = ccbaby(x=x, runlength=runlength, seed=seed) + + responses = { + "total_cost": fn, + "constraint_lhs_by_hour": constraint, + } + + return responses, {} + +#optimization problem +class CallCenterMinCost(Problem): + """Simulation-optimization problem for call center staffing.""" + + class_name_abbr = "CALLCENTER-1" + class_name = "Minimize Staffing Cost Subject to Service Levels" + + model_class = CallCenter + n_objectives = 1 + n_stochastic_constraints = 16 + + minmax = (-1,) + constraint_type = ConstraintType.STOCHASTIC + variable_type = VariableType.DISCRETE + gradient_available = False + + optimal_value = None + optimal_solution = None + + model_default_factors = {} + model_decision_factors = {"staffing_matrix"} + + @property + def dim(self): + #17 start times: 0, 0.5, ..., 8 + #4 lunch options: 3, 3.5, 4, 4.5 + return 17 * 4 + + @property + def lower_bounds(self): + return (0,) * self.dim + + @property + def upper_bounds(self): + return (100,) * self.dim + + def vector_to_factor_dict(self, vector): + """ + Convert flat optimization vector into 17 x 4 staffing matrix. + """ + x_matrix = np.array(vector, dtype=int).reshape((17, 4)) + return { + "staffing_matrix": x_matrix, + } + + def factor_dict_to_vector(self, factor_dict): + """ + Convert 17 x 4 staffing matrix back into flat vector. + """ + return tuple(np.array(factor_dict["staffing_matrix"]).flatten()) + + def replicate(self, _x): + """ + Run one simulation replication at solution x. + """ + responses, _ = self.model.replicate() + + objective = responses["total_cost"] + + stochastic_constraints = responses["constraint_lhs_by_hour"] + + return RepResult( + objectives=[Objective(stochastic=objective)], + stochastic_constraints=[ + Objective(stochastic=value) for value in stochastic_constraints + ], + ) + + #checks simple bounds and integrality + def check_deterministic_constraints(self, x): + + if not super().check_deterministic_constraints(x): + return False + + x_matrix = np.array(x).reshape((17, 4)) + + # Optional: total staffing cap + if np.sum(x_matrix) > 1000: + return False + + return True + + #generate random feasible staffing plan + def get_random_solution(self, rand_sol_rng): + + return tuple( + rand_sol_rng.integer_random_vector_from_simplex( + n_elements=self.dim, + summation=200, + with_zero=True, + ) + ) + +#function making number of agents dependent on demand +#this makes an alternative ("better") staffing plan than the random feasible one above +def make_demand_based_x(): + hk = np.array([3.0, 3.5, 4.0, 4.5]) + start_times = np.arange(0, 8.5, 0.5) + + x = np.zeros((17, 4), dtype=int) + + #estimated number of agents needed by half-hour + times = np.arange(0, 16.5, 0.5) + demand = 500 + 500 * np.sin((3 * np.pi * times - 16 * np.pi) / 32) + + #each agent handles about 10 calls/hour if avg service time is 6 minutes + required_agents = np.ceil(demand / 10).astype(int) + + #add buffer + required_agents = np.ceil(required_agents * 1.25).astype(int) + + #put more shift starts before the peak + starts_by_time = np.array([ + 2, 4, 6, 8, 10, 12, 14, 16, 18, + 18, 16, 14, 12, 10, 8, 6, 4 + ]) + + #split each start time across the 4 lunch options + for j in range(17): + base = starts_by_time[j] // 4 + remainder = starts_by_time[j] % 4 + x[j, :] = base + x[j, :remainder] += 1 + + return x + +#plot staffing vs. demand to test +def plot_staffing_vs_demand(x): + x = np.asarray(x, dtype=int) + + hk = np.array([3.0, 3.5, 4.0, 4.5]) + times = np.arange(0, 16.5, 0.5) + + working_agents = [] + + for t in times: + count = 0 + + for j in range(17): + start_time = j / 2 + + for k in range(4): + lunch_start = start_time + hk[k] + lunch_end = lunch_start + 0.5 + shift_end = start_time + 8.5 + + #agent is working if in shift, but not at lunch + is_on_shift = start_time <= t < shift_end + is_at_lunch = lunch_start <= t < lunch_end + + if is_on_shift and not is_at_lunch: + count += x[j, k] + + working_agents.append(count) + + working_agents = np.array(working_agents) + + #arrival-rate demand function from simulation + demand = 500 + 500 * np.sin((3 * np.pi * times - 16 * np.pi) / 32) + + #plot agents working vs. demand + plt.figure(figsize=(10, 5)) + plt.plot(times, working_agents, marker="o", label="Agents working") + plt.plot(times, demand, marker="x", label="Demand / arrival rate") + + plt.xlabel("Time of day") + plt.ylabel("Count") + plt.title("Agents Working vs. Demand by Half-Hour") + plt.legend() + plt.grid(True) + plt.show() + +#test +if __name__ == "__main__": + + x = make_demand_based_x() + + print("x matrix:") + print(x) + print("total agents:", np.sum(x)) + + plot_staffing_vs_demand(x) + + fn, FnVar, FnGrad, FnGradCov, constraint, ConstraintCov, ConstraintGrad, ConstraintGradCov = ccbaby( + x=x, + runlength=2, + seed=123 + ) + + print("Objective value:", fn) + print("Constraints:", constraint) \ No newline at end of file From 24737a7626aa65fa9bc4a29398373b3666bed98d Mon Sep 17 00:00:00 2001 From: chloetsetsekas Date: Mon, 11 May 2026 18:16:17 -0400 Subject: [PATCH 2/3] Remove unrelated files from PR --- docs/source/models/tollnw.rst | 164 ---------------------------------- notebooks/demo_problem.py | 140 ----------------------------- 2 files changed, 304 deletions(-) delete mode 100644 docs/source/models/tollnw.rst delete mode 100644 notebooks/demo_problem.py diff --git a/docs/source/models/tollnw.rst b/docs/source/models/tollnw.rst deleted file mode 100644 index 938d8169b..000000000 --- a/docs/source/models/tollnw.rst +++ /dev/null @@ -1,164 +0,0 @@ -Toll Networks -======================= -See the :mod:'simopt.models.tollnw' module for API details. - -Model: Toll Road Improvements in a Newtork ----------------------------------------- - -Description -^^^^^^^^^^^ - -This model represents a connected directed graph :math:`G= (V,E)` describing a network -of toll roads. Customers wishing to travel from point :math:`i` to point :math:`j` -arrive according to a Poisson process with rate :math:`lambda_{ij}`. - -Each directed edge represents a toll road with price :math:`p_{ij}`. Travel times on each -road are triangularly distributed with route-specific parameters. The network operates as -an open Jackson network with interconnected queues. - -Each road has the capacity for one vehicle at a time. If a vehicle is currently traveling -on a road, arriving vehicles form a queue and wait until the road becomes available. - -The operator can choose to invest in road improvements. If :math:`x_{ij}` hundred dollars -is spent improving road :math:`(i,j)`, the resulting mode of the triangular travel-time -distribution becomes: - -:math:`a_{ij} + (b_{ij} - a_{ij}) \exp(-x_{ij})` - -where: -- :math:`a_{ij}` is the minimum travel time -- :math:`b_{ij}` is the maximum travel time - - -If no investment is made (i.e. :math:`x_{ij} = 0`), the mode equals :math:`b_{ij}`. - -The operator seeks to maximize expected profit over a finite horizon :math:`T`: - -:math:`\mathbb{E}[\text{Profit}] = \sum_{(i,j) \in E} p_{ij} \mathbb{E}[N_{ij}(T)] - \sum_{(i,j) \in E} x_{ij}` - -where: -:math:`N_{ij}(T)` is the number of trips completed on road :math:`(i,j)` before time :math:`T`. - - - -Sources of Randomness -^^^^^^^^^^^^^^^^^^^^^ - -There are three source of randomness: - -1. External arrivals modeled as Poisson processes with rates :math:`lambda_{ij}` -2. Internal routing decisions governed by routing matrix :math:`R` -3. Travel times following triangular distribution - - -Model Factors -^^^^^^^^^^^^^ - -- adjacency_matrix: graph structure describing road connectivity -- lambda_ij: external Poisson arrival rates - - default: :math:`\lambda_{ij} = i + j` for :math:`i \neg j`, :math:`i, j \in \{1,2,3\}` -- routing_matrix: internal routing matrix with exit probabilities -- pij: toll/amount charged per completed trips - - default: 1 for all roads -- aij: minimum travel time of triangular distribution - - default: 0 -- bij: maximum travel time of triangular distribution - - default: 1 -- investment_levels: investment decisions :math:`x_{ij}` -- T: length of operating horizon - - default: 96 (12 hours) -- time_unit: one unit equals 10 minutes - - -Recommended Parameter Settings -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Let :math:`G` be a complete graph with three vertices: - -..math:: - - A = - \begin{bmatrix} - 0 & 1 & 1 \\ - 1 & 0 & 1 \\ - 1 & 1 & 0 - \end{bmatrix} - -Routing Matrix: - -..math:: - - R = - \begin{bmatrix} - 0 & 1/3 & 1/3 & 1/3 - 1/2 & 0 & 1/4 & 1/4 - 1/5 & 1/5 & 0 & 3/5 - \end{bmatrix} - -The last column represents exiting the network. - -Responses -^^^^^^^^^ - -- total_profit: expected profit over horizon :math:`T` -- total_trips_completed: total completed trips across all roads -- road_trip_counts: vector of :math:`N_{ij}(T)` values -- average_queue_lengths: time-averaged queue lengths per road - -References -^^^^^^^^^^ - -Eckman, D. (2017). Toll Road Improvements in a Network. - -Optimization Problem: Maximize Expected Profit (TOLLNETWORK-1) --------------------------------------------------------------- - -Decision Variables -^^^^^^^^^^^^^^^^^^ - -- investment_levels: :math:`x_{ij}` (continuous, :math:`\ge 0`) - -Objective -^^^^^^^^^^ - -Maximize expected profit over horizon :math:`T`. - -Constraints -^^^^^^^^^^^ - -- :math:`x_{ij} \ge 0` for all roads -- Continuous decision Variables xij>=0 - -Problem Factors -^^^^^^^^^^^^^^^ - -- Budget: maximum number of replications for solver - - default: 1000 - -Fixed Model Factors -^^^^^^^^^^^^^^^^^^^ - -- aij = 0 -- bij = 1 -- pij = 1 -- T = 96 - -Starting Solution -^^^^^^^^^^^^^^^^^ - -- :math:`x_{ij} = 1` for all roads - -Random Solutions -^^^^^^^^^^^^^^^^ - -Generate continuous nonnegative investment vectors over all roads - -Optimal Solution -^^^^^^^^^^^^^^^^ - -Unknown - -Optimal Objective Function Value -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Unknown \ No newline at end of file diff --git a/notebooks/demo_problem.py b/notebooks/demo_problem.py deleted file mode 100644 index 93ceaa47f..000000000 --- a/notebooks/demo_problem.py +++ /dev/null @@ -1,140 +0,0 @@ -# --- -# jupyter: -# jupytext: -# formats: ipynb,py:percent -# text_representation: -# extension: .py -# format_name: percent -# format_version: '1.3' -# jupytext_version: 1.17.3 -# kernelspec: -# display_name: simopt -# language: python -# name: python3 -# --- - -# %% [markdown] -# # Demo script for the Problem class. -# -# This script is intended to help with debugging a problem. -# -# It imports a problem, initializes a problem object with given factors, sets up pseudorandom number generators, and runs multiple replications at a given solution. - -# %% [markdown] -# ## Append SimOpt Path -# -# Since the notebook is stored in simopt/notebooks, we need to append the parent simopt directory to the system path to import the necessary modules later on. - -# %% -import sys -from pathlib import Path - -# Take the current directory, find the parent, and add it to the system path -sys.path.append(str(Path.cwd().parent)) - -# %% [markdown] -# ## Configuration Parameters -# -# This section defines the core parameters for the demo. -# -# To query model/problem/solver names, run `python scripts/list_directories.py` - -# %% -# Import problem. -# from models. import -# Replace with name of .py file containing problem class. -# Replace with name of problem class. -# Fix factors of problem. Specify a dictionary of factors. -# fixed_factors = {} # Resort to all default values. -# Look at Problem class definition to get names of factors. -# Initialize an instance of the specified problem class. -# myproblem = (fixed_factors=fixed_factors) -# Replace with name of problem class. -# Initialize a solution x corresponding to the problem. -# x = (,) -# Look at the Problem class definition to identify the decision variables. -# x will be a tuple consisting of the decision variables. -# The following line does not need to be changed. -# mysolution = Solution(x, myproblem) -# Working example for CntNVMaxProfit problem. -# ----------------------------------------------- -from simopt.models.cntnv import CntNVMaxProfit - -fixed_factors = {"initial_solution": (2,), "budget": 500} -myproblem = CntNVMaxProfit(fixed_factors=fixed_factors) -x = (1,) - -# ----------------------------------------------- - -# Another working example for FacilitySizingTotalCost problem. (Commented out) -# This example has stochastic constraints. -# ----------------------------------------------- -# from models.facilitysizing import FacilitySizingTotalCost -# fixed_factors = {"epsilon": 0.1} -# myproblem = FacilitySizingTotalCost(fixed_factors=fixed_factors) -# x = (200, 200, 200) -# mysolution = Solution(x, myproblem) -# ----------------------------------------------- - -num_macroreps = 1500 - -# %% -# Create and attach rngs to solution -from mrg32k3a.mrg32k3a import MRG32k3a -from simopt.base import Solution - -rng_list = [MRG32k3a(s_ss_sss_index=[0, ss, 0]) for ss in range(myproblem.model.n_rngs)] -mysolution = Solution(x, myproblem) - -mysolution.attach_rngs(rng_list, copy=False) - -# %% -# Simulate a fixed number of replications at the solution x. -myproblem.simulate(mysolution, num_macroreps=num_macroreps) - -# %% -# Print results to console. -print( - f"Ran {num_macroreps} replications of the {myproblem.name} problem " - f"at solution x = {x}.\n" -) -obj_mean = round(mysolution.objectives_mean[0], 4) -obj_stderr = round(mysolution.objectives_stderr[0], 4) -print(f"The mean objective estimate was {obj_mean} with standard error {obj_stderr}.") -print("The individual observations of the objective were:") -for idx in range(num_macroreps): - print(f"\t {round(mysolution.objectives[idx][0], 4)}") -if myproblem.gradient_available: - print("\nThe individual observations of the gradients of the objective were:") - for idx in range(num_macroreps): - grads = mysolution.objectives_gradients[idx][0] - print(f"\t {[round(g, 4) for g in grads]}") -else: - print("\nThis problem has no known gradients.") -if myproblem.n_stochastic_constraints > 0: - print( - f"\nThis problem has {myproblem.n_stochastic_constraints} stochastic " - "constraints of the form E[LHS] <= 0." - ) - for stc_idx in range(myproblem.n_stochastic_constraints): - stoch_const_mean = mysolution.stoch_constraints_mean[stc_idx] - stoch_const_mean_round = round(stoch_const_mean, 4) - stoch_const_stderr = mysolution.stoch_constraints_stderr[stc_idx] - stoch_const_stderr_round = round(stoch_const_stderr, 4) - print( - f"\tFor stochastic constraint #{stc_idx + 1}, " - f"the mean of the LHS was {stoch_const_mean_round} " - f"with standard error {stoch_const_stderr_round}." - ) - print("\tThe observations of the LHSs were:") - for idx in range(num_macroreps): - # Quick check to make sure the stochastic constraints are not None. - if mysolution.stoch_constraints is None: - error_msg = "Stochastic constraints should not be None." - raise ValueError(error_msg) - # Print out the current stochastic constraint value. - stoch_const = mysolution.stoch_constraints[stc_idx][idx] - stoch_const_round = round(stoch_const, 4) - print(f"\t\t {stoch_const_round}") -else: - print("\nThis problem has no stochastic constraints.") From a916f86300ea4ec9ec1508197635d89061bc73f0 Mon Sep 17 00:00:00 2001 From: chloetsetsekas Date: Mon, 11 May 2026 18:22:29 -0400 Subject: [PATCH 3/3] Revert demo_problem changes --- notebooks/demo_problem.py | 140 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 140 insertions(+) create mode 100644 notebooks/demo_problem.py diff --git a/notebooks/demo_problem.py b/notebooks/demo_problem.py new file mode 100644 index 000000000..f33615332 --- /dev/null +++ b/notebooks/demo_problem.py @@ -0,0 +1,140 @@ +# --- +# jupyter: +# jupytext: +# formats: ipynb,py:percent +# text_representation: +# extension: .py +# format_name: percent +# format_version: '1.3' +# jupytext_version: 1.17.3 +# kernelspec: +# display_name: simopt +# language: python +# name: python3 +# --- + +# %% [markdown] +# # Demo script for the Problem class. +# +# This script is intended to help with debugging a problem. +# +# It imports a problem, initializes a problem object with given factors, sets up pseudorandom number generators, and runs multiple replications at a given solution. + +# %% [markdown] +# ## Append SimOpt Path +# +# Since the notebook is stored in simopt/notebooks, we need to append the parent simopt directory to the system path to import the necessary modules later on. + +# %% +import sys +from pathlib import Path + +# Take the current directory, find the parent, and add it to the system path +sys.path.append(str(Path.cwd().parent)) + +# %% [markdown] +# ## Configuration Parameters +# +# This section defines the core parameters for the demo. +# +# To query model/problem/solver names, run `python scripts/list_directories.py` + +# %% +# Import problem. +# from models. import +# Replace with name of .py file containing problem class. +# Replace with name of problem class. +# Fix factors of problem. Specify a dictionary of factors. +# fixed_factors = {} # Resort to all default values. +# Look at Problem class definition to get names of factors. +# Initialize an instance of the specified problem class. +# myproblem = (fixed_factors=fixed_factors) +# Replace with name of problem class. +# Initialize a solution x corresponding to the problem. +# x = (,) +# Look at the Problem class definition to identify the decision variables. +# x will be a tuple consisting of the decision variables. +# The following line does not need to be changed. +# mysolution = Solution(x, myproblem) +# Working example for CntNVMaxProfit problem. +# ----------------------------------------------- +from simopt.models.cntnv import CntNVMaxProfit + +fixed_factors = {"initial_solution": (2,), "budget": 500} +myproblem = CntNVMaxProfit(fixed_factors=fixed_factors) +x = (3,) + +# ----------------------------------------------- + +# Another working example for FacilitySizingTotalCost problem. (Commented out) +# This example has stochastic constraints. +# ----------------------------------------------- +# from models.facilitysizing import FacilitySizingTotalCost +# fixed_factors = {"epsilon": 0.1} +# myproblem = FacilitySizingTotalCost(fixed_factors=fixed_factors) +# x = (200, 200, 200) +# mysolution = Solution(x, myproblem) +# ----------------------------------------------- + +num_macroreps = 10 + +# %% +# Create and attach rngs to solution +from mrg32k3a.mrg32k3a import MRG32k3a +from simopt.base import Solution + +rng_list = [MRG32k3a(s_ss_sss_index=[0, ss, 0]) for ss in range(myproblem.model.n_rngs)] +mysolution = Solution(x, myproblem) + +mysolution.attach_rngs(rng_list, copy=False) + +# %% +# Simulate a fixed number of replications at the solution x. +myproblem.simulate(mysolution, num_macroreps=num_macroreps) + +# %% +# Print results to console. +print( + f"Ran {num_macroreps} replications of the {myproblem.name} problem " + f"at solution x = {x}.\n" +) +obj_mean = round(mysolution.objectives_mean[0], 4) +obj_stderr = round(mysolution.objectives_stderr[0], 4) +print(f"The mean objective estimate was {obj_mean} with standard error {obj_stderr}.") +print("The individual observations of the objective were:") +for idx in range(num_macroreps): + print(f"\t {round(mysolution.objectives[idx][0], 4)}") +if myproblem.gradient_available: + print("\nThe individual observations of the gradients of the objective were:") + for idx in range(num_macroreps): + grads = mysolution.objectives_gradients[idx][0] + print(f"\t {[round(g, 4) for g in grads]}") +else: + print("\nThis problem has no known gradients.") +if myproblem.n_stochastic_constraints > 0: + print( + f"\nThis problem has {myproblem.n_stochastic_constraints} stochastic " + "constraints of the form E[LHS] <= 0." + ) + for stc_idx in range(myproblem.n_stochastic_constraints): + stoch_const_mean = mysolution.stoch_constraints_mean[stc_idx] + stoch_const_mean_round = round(stoch_const_mean, 4) + stoch_const_stderr = mysolution.stoch_constraints_stderr[stc_idx] + stoch_const_stderr_round = round(stoch_const_stderr, 4) + print( + f"\tFor stochastic constraint #{stc_idx + 1}, " + f"the mean of the LHS was {stoch_const_mean_round} " + f"with standard error {stoch_const_stderr_round}." + ) + print("\tThe observations of the LHSs were:") + for idx in range(num_macroreps): + # Quick check to make sure the stochastic constraints are not None. + if mysolution.stoch_constraints is None: + error_msg = "Stochastic constraints should not be None." + raise ValueError(error_msg) + # Print out the current stochastic constraint value. + stoch_const = mysolution.stoch_constraints[stc_idx][idx] + stoch_const_round = round(stoch_const, 4) + print(f"\t\t {stoch_const_round}") +else: + print("\nThis problem has no stochastic constraints.")