simulation_analysis.capacity

  1"""
  2This module provides functions to analyze the capacity of a queueing system,
  3calculate the ultimate capacity, and solve differential equations related to queue dynamics.
  4"""
  5import numpy as np
  6from scipy.optimize import minimize
  7from scipy.integrate import solve_ivp
  8from scipy.optimize import minimize
  9import numpy as np
 10import matplotlib.pyplot as plt
 11import sys
 12import os
 13import glob
 14
 15
 16# parent_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))
 17
 18# # Add parent directory to sys.path
 19# sys.path.append(parent_dir)
 20
 21from constants import SCENARIO_NAME, NUM_MONTHS, ARRIVAL_INCREASE_FACTOR
 22
 23
 24
 25def calculate_mean_wait_time(arrival_rates, queue_lengths, input_mean_wait_time, logfilename):
 26    """
 27    Calculate the mean wait time for a queueing system given arrival rates and queue lengths.
 28
 29    Args:
 30        arrival_rates (list): List of 3 arrival rates for the queues.
 31        queue_lengths (list): List of 3 queue lengths for the queues.
 32        input_mean_wait_time (float): The input mean wait time to compare against.
 33
 34    Returns:
 35        dict: A dictionary with calculated service rates, mean wait time, and relative error.
 36    """
 37    print("\nCapacity Analysis")
 38    print("Arrival Rates: ", arrival_rates)
 39    print("Queue Lengths: ", queue_lengths)
 40    with open(logfilename, 'a') as f:
 41        f.write("\nCapacity Analysis\n")
 42        f.write("Arrival Rates: " + str(arrival_rates) + "\n")
 43        f.write("Queue Lengths: " + str(queue_lengths) + "\n")
 44    arrival_rates = np.array(arrival_rates)
 45    queue_lengths = np.array(queue_lengths)
 46    eps = 1e-6
 47
 48    # Define the fitting function
 49    def fitting_function(mu, arrival_rates, queue_lengths):
 50        numerator = np.sum(arrival_rates / mu**2)
 51        denominator = 1 - (arrival_rates / mu).sum()
 52        return (numerator / denominator) * arrival_rates - queue_lengths
 53
 54    # Objective function for minimization
 55    def objective_function(mu, arrival_rates, queue_lengths):
 56        return np.sum(fitting_function(mu, arrival_rates, queue_lengths)**2)
 57
 58    # Constraint to ensure the stability condition
 59    def stability_constraint(mu, arrival_rates):
 60        return 1 - (arrival_rates / mu).sum()
 61
 62    # Initial guess for service rate
 63    initial_guess = [1.1 * float(arrival_rates.sum())]
 64
 65    # Define constraints
 66    constraints = {
 67        'type': 'ineq',
 68        'fun': stability_constraint,
 69        'args': (arrival_rates,)
 70    }
 71
 72    # Perform optimization
 73    result = minimize(
 74        objective_function,
 75        initial_guess,
 76        method='SLSQP',
 77        args=(arrival_rates, queue_lengths),
 78        constraints=constraints,
 79        options={'ftol': 1e-12, 'maxiter': 1000}
 80    )
 81
 82    if result.success:
 83        mu = result.x[0]
 84        print(f"Calculated service rate (mu): {mu:.2f}")
 85        with open(logfilename, 'a') as f:
 86            f.write(f"Calculated service rate (mu): {mu:.2f}\n")
 87    else:
 88        raise ValueError("Optimization failed: " + result.message)
 89
 90    # Calculate mean wait time
 91    numerator_W_q = np.sum(arrival_rates / mu**2)
 92    denominator_W_q = 1 - np.sum(arrival_rates / mu)
 93    calculated_mean_wait_time = numerator_W_q / denominator_W_q
 94
 95    print(f"Calculated mean wait time: {calculated_mean_wait_time:.2f}")
 96    print(f"Input mean wait time: {input_mean_wait_time:.2f}")
 97    with open(logfilename, 'a') as f:
 98        f.write(
 99            f"Calculated mean wait time: {calculated_mean_wait_time:.2f}\n")
100        f.write(f"Input mean wait time: {input_mean_wait_time:.2f}\n")
101
102    # Calculate relative error
103    relative_error = (calculated_mean_wait_time -
104                      input_mean_wait_time) / input_mean_wait_time * 100
105    print(f"Relative error (%): {relative_error:.2f}")
106
107    # Return results as a dictionary
108    output = {
109        'Service Rate': mu,
110        'Calculated Mean Wait Time': calculated_mean_wait_time,
111        'Input Mean Wait Time': input_mean_wait_time,
112        'Relative Error (%)': relative_error
113    }
114
115    print("Capacity Analysis Complete")
116    print("Relative Error: ", relative_error)
117    with open(logfilename, 'a') as f:
118        f.write("Relative Error: " + str(relative_error) + "\n")
119    return output
120
121
122def dN_dlambda(lambda_, T, N, Cu, theta):
123    """
124    Differential equation for N with respect to lambda.
125    Args:
126        lambda_ (float): The current value of lambda.
127        T (float): Total time in hours.
128        N (float): Current value of N.
129        Cu (float): Ultimate capacity.
130        theta (float): Degree of congestion.
131    Returns:
132        float: The derivative of N with respect to lambda.
133    """
134    return T * (1 - (N / (T * Cu))**theta)
135
136
137def solve_N(lambda_eval, Cu, T, theta):
138    """
139    Solve the ODE for N using the initial condition N(0) = 0.
140    Args:
141        lambda_eval (array): Array of lambda values for evaluation.
142        Cu (float): Ultimate capacity.
143        T (float): Total time in hours.
144        theta (float): Degree of congestion.
145    Returns:
146        array: Solution for N at the given lambda values.
147    """
148    sol = solve_ivp(
149        fun=lambda lam, N: dN_dlambda(lam, T, N, Cu, theta),
150        t_span=(0, lambda_eval[-1]),
151        y0=[0],
152        t_eval=lambda_eval,
153        method='RK45',
154        rtol=1e-8,
155        atol=1e-8
156    )
157    return sol.y[0]
158
159
160def objective(params, constant_args):
161    """
162    Objective function to minimize the difference between predicted and actual N values.
163    Args:
164        params (list): List containing Cu and theta.
165        constant_args (list): List containing a, b, c, d, and T.
166    Returns:
167        float: Sum of squared differences between predicted and actual N values.
168    """
169    Cu, theta = params
170    a, b, c, d, T = constant_args
171    lambda_vals = np.array([a, c])
172    N_vals = np.array([b, d])
173    N_pred = solve_N(lambda_vals, Cu, T, theta)
174    return np.sum((N_pred - N_vals)**2)
175
176def calculate_ultimate_capacity() -> tuple[float, float]:
177    """
178    Auto‐discover all scenario logs (BASE + others), but only read
179    'Operating capacity:' from the BASE file. Fit dN/dλ = T*(1 - (N/(Cu*T))**θ),
180    write results, and plot every sample labeled.
181    """
182    # 1) Setup
183    T = NUM_MONTHS * 30 * 24
184    LOG_DIR = 'bottleneckAnalysis/logs'
185
186    # 2) Identify files
187    base_file = os.path.join(LOG_DIR, f'{SCENARIO_NAME}_BASE.txt')
188    pattern   = os.path.join(LOG_DIR, f'{SCENARIO_NAME}_*.txt')
189    all_files = sorted(glob.glob(pattern))
190    high_files = [f for f in all_files if f != base_file]
191
192    # 3) Read BASE file
193    lambdas = []
194    Ns      = []
195    Ns_actual = []
196    op_capacity = None
197
198
199    with open(base_file) as f:
200        lam0 = N0 = None
201        for line in f:
202            if line.startswith('Operating capacity:'):
203                op_capacity = float(line.split()[2])
204            elif line.startswith('Arrival rate:'):
205                lam0 = float(line.split()[2])
206            elif line.startswith('Exited ships:'):
207                N0 = float(line.split()[2])
208        if lam0 is None or N0 is None:
209            raise RuntimeError(f"Could not parse λ or N from BASE file {base_file}")
210        lambdas.append(lam0)
211        Ns.append(N0)
212        Ns_actual.append(N0)
213
214    # 4) Read every non-BASE file (only λ and N)
215    for fn in high_files:
216        lam = Ni = None
217        with open(fn) as f:
218            for line in f:
219                if line.startswith('Arrival rate:'):
220                    lam = float(line.split()[2])
221                elif line.startswith('Exited ships:'):
222                    Ni = float(line.split()[2])
223        if lam is not None and Ni is not None:
224            lambdas.append(lam)
225            Ns.append(Ni)
226            Ns_actual.append(Ni)
227        else:
228            print(f"Skipping {fn}: missing Arrival rate or Exited ships")
229
230        
231    lambdas = np.array(lambdas)
232    Ns      = np.array(Ns)
233
234    sort_idx = np.argsort(lambdas)
235    lambdas = lambdas[sort_idx]
236    Ns = Ns[sort_idx]
237
238    print("Discovered sample points:")
239    for i, (λ, N) in enumerate(zip(lambdas, Ns), start=1):
240        print(f"  sample point {i}: λ = {λ:.2f},  N = {N:.2f}")
241
242    # 5) Define multi‐point objective
243    def multi_obj(x, lam_arr, N_arr, T):
244        Cu, θ = x
245        N_pred = solve_N(lam_arr, Cu, T, θ)
246        return np.sum((N_pred - N_arr) ** 2)
247
248    # 6) Optimize
249    initial = [max(Ns)/T, 5.0]
250    bounds  = [(0.1, max(lambdas)), (0.1, 20.0)] #TODO: adjust bounds if needed
251    res = minimize(multi_obj, x0=initial, args=(lambdas, Ns, T), bounds=bounds)
252    Cu_opt, theta_opt = res.x
253    print(f"\nOptimal Cu = {Cu_opt:.4f},  θ = {theta_opt:.4f}")
254
255    # 7) Append to results.txt
256    out_path = 'bottleneckAnalysis/results.txt'
257    with open(out_path, 'a') as f:
258        f.write(f"===== Capacity Analysis for {SCENARIO_NAME} =====\n")
259        f.write(f"Months: {NUM_MONTHS} (T = {T} hours)\n")
260        f.write("Sample points (λ_i, N_i):\n")
261        for i, (λ, N) in enumerate(zip(lambdas, Ns), start=1):
262            f.write(f"  Point {i}: λ = {λ:.4f}, N = {N:.4f}\n")
263        if op_capacity is not None:
264            f.write(f"Operating capacity (Co): {op_capacity:.4f} vessels/hr\n")
265        f.write(f"Fitted ultimate capacity (Cu): {Cu_opt:.4f}\n")
266        f.write(f"Fitted exponent (θ): {theta_opt:.4f}\n\n")
267
268    # 8) Plot everything
269    λ_range = np.linspace(0, max(lambdas)*1.1, 300)
270    N_fit   = solve_N(λ_range, Cu_opt, T, theta_opt)
271
272    plt.figure(figsize=(6, 5))
273    plt.plot(λ_range, N_fit, label='ODE fit', linewidth=2, color='orange')
274
275    for i, (λ, N) in enumerate(zip(lambdas, Ns), start=1):
276        if i == 1:
277            plt.scatter(λ, N, s=20, label=f'Simulation runs', c='blue')
278        else:
279            plt.scatter(λ, N, s=20, c='blue')
280        # plt.annotate(f"sample point {i}",
281        #              (λ, N),
282        #              textcoords="offset points",
283        #              xytext=(5, 5),
284        #              fontsize=9)
285
286    plt.axhline(y=Cu_opt*T, linestyle='--', label='Ultimate capacity', color = 'magenta')
287    # find y corresponding to operating capacity
288
289    try:
290        y = solve_N([op_capacity], Cu_opt, T, theta_opt) if op_capacity is not None else None
291        plt.axvline(x=op_capacity, linestyle='--', label='Operating capacity', color ='green')
292        plt.axhline(y=y, linestyle='--', color='green')
293
294    except Exception as e:
295        print(f"Error occurred while finding y for operating capacity: {e}")
296
297    # set limits on y axis from 0 to 6500
298    plt.ylim(0, 7000)
299    plt.xlabel(r'$\lambda$', fontsize=14)
300    plt.ylabel(r'$N^{exits}(\lambda)$', fontsize=14)
301    plt.xticks(fontsize=12)
302    plt.yticks(fontsize=12)
303    plt.legend()
304    plt.tight_layout()
305    plt.savefig(f"bottleneckAnalysis/{SCENARIO_NAME}_capacity_analysis.pdf")
306    plt.close()
307
308    return Cu_opt, theta_opt
309# calculate_ultimate_capacity()