""" robust_gne_experiment.py Reproducible numerical experiment accompanying: "Robust Generalized Nash Equilibria under Uncertain Shared Constraints: A Variational Approach to Decentralized Resource Allocation" Soumaya Et-tryky, Khalid El Hajioui Computes the robust variational equilibrium path for the paper's five-player, two-resource quadratic potential game under Bertsimas-Sim budgeted uncertainty (Section 5), certifies in-set feasibility (Corollary 4.5), runs the diagnostic out-of-set Beta(2,5) stress test (Section 5.3), and regenerates Tables 1-2 and Figures 1-3. Usage: python robust_gne_experiment.py Requires: numpy, scipy, matplotlib (tested with numpy>=1.24, scipy>=1.10) Outputs (written relative to the current working directory): results/equilibria.csv Full equilibrium path: Gamma, x1..x5, Total, Phi, Pi results/stress_test.csv Stress diagnostics: Gamma, Total, Pi, Violation_rate_pct results/model_parameters.json Q, r, A_bar, D, b, grid, seed, eigenvalue figures/fig1_equilibrium_path.{pdf,png} figures/fig2_robustness_premium.{pdf,png} figures/fig3_stress_test.{pdf,png} """ import json import os import numpy as np from scipy.optimize import minimize import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt # --------------------------------------------------------------------------- # 1. Model parameters (paper, Section 5.1) # --------------------------------------------------------------------------- N_PLAYERS = 5 N_RESOURCES = 2 X_UB = 3.0 # local bound 0 <= x_i <= 3 STRESS_SEED = 20261007 # fixed seed for the Beta(2,5) diagnostic stress test N_INSET_SAMPLES = 20_000 # budget-feasible certification draws per Gamma N_STRESS_SAMPLES = 50_000 # out-of-set diagnostic draws per Gamma VIOLATION_TOL = 5e-7 # feasibility tolerance (Section 5.2) GAMMA_GRID = np.round(np.arange(0.0, 3.0 + 1e-9, 0.25), 2) # Figures 1-3 resolution TABLE_GAMMAS = [0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0] # Tables 1-2 resolution Q = np.array([ [2.80, 0.25, 0.18, 0.12, 0.15], [0.25, 2.40, 0.20, 0.14, 0.10], [0.18, 0.20, 2.60, 0.22, 0.13], [0.12, 0.14, 0.22, 2.20, 0.18], [0.15, 0.10, 0.13, 0.18, 2.50], ]) r = np.array([5.4, 4.9, 5.1, 4.7, 5.0]) A_BAR = np.array([ [1.0, 0.8, 1.1, 0.9, 1.0], [0.6, 1.2, 0.7, 1.1, 0.9], ]) D = np.array([ [0.22, 0.18, 0.25, 0.20, 0.23], [0.15, 0.28, 0.18, 0.24, 0.21], ]) B_CAP = np.array([7.0, 5.8]) n, m = N_PLAYERS, N_RESOURCES # --------------------------------------------------------------------------- # 2. Potential, robust counterpart, and solver (Section 3.3, 4.3, 5.2) # --------------------------------------------------------------------------- def potential(x): """Phi(x) = 0.5 x^T Q x - r^T x (Eq. 5.3).""" return 0.5 * x @ Q @ x - r @ x def solve_robust_equilibrium(Gamma, x0=None): """ Solve the lifted convex program (Eq. 5.4): min_{x,p,q} 0.5 x^T Q x - r^T x s.t. 0 <= x_i <= X_UB a_bar_r^T x + Gamma*p_r + sum_i q_ri <= b_r (3.9a) p_r + q_ri >= d_ri * x_i, p_r >= 0, q_ri >= 0 (3.9b) Returns the unique equilibrium allocation x (length n). """ nvar = n + m + m * n def unpack(z): return z[:n], z[n:n + m], z[n + m:].reshape(m, n) def obj(z): x, _, _ = unpack(z) return potential(x) def obj_grad(z): x, _, _ = unpack(z) return np.concatenate([Q @ x - r, np.zeros(m), np.zeros(m * n)]) cons = [] for rr in range(m): def cap_con(z, rr=rr): x, p, q = unpack(z) return B_CAP[rr] - (A_BAR[rr] @ x + Gamma * p[rr] + np.sum(q[rr])) cons.append({'type': 'ineq', 'fun': cap_con}) for i in range(n): def dual_con(z, rr=rr, i=i): x, p, q = unpack(z) return p[rr] + q[rr, i] - D[rr, i] * x[i] cons.append({'type': 'ineq', 'fun': dual_con}) bounds = [(0, X_UB)] * n + [(0, None)] * m + [(0, None)] * (m * n) if x0 is None: x0 = np.zeros(nvar) x0[:n] = 1.0 res = minimize(obj, x0, jac=obj_grad, bounds=bounds, constraints=cons, method='SLSQP', options={'maxiter': 2000, 'ftol': 1e-14}) if not res.success: raise RuntimeError(f"SLSQP failed to converge at Gamma={Gamma}: {res.message}") x, _, _ = unpack(res.x) return x def sigma_support(dx, Gamma): """ Exact Bertsimas-Sim support function (Eq. 3.7), closed form: sort d_ri*x_i descending, sum the top floor(Gamma) values, plus the fractional part of Gamma times the next-largest value. """ s = np.sort(dx)[::-1] k = int(np.floor(Gamma)) frac = Gamma - k total = s[:k].sum() if k > 0 else 0.0 if k < len(s): total += frac * s[k] return total def robust_lhs(x, Gamma): """Exact robust left-hand side a_bar_r^T x + sigma_r(x;Gamma) for each resource.""" return np.array([A_BAR[rr] @ x + sigma_support(D[rr] * x, Gamma) for rr in range(m)]) # --------------------------------------------------------------------------- # 3. In-set certification and out-of-set stress test (Section 5.3, Cor. 4.5) # --------------------------------------------------------------------------- def sample_inset_disturbances(rng, Gamma, n_samples): """ Sample budget-feasible disturbances in U_r(Gamma) = {0<=z_i<=1, sum z_i<=Gamma}: draw uniformly in [0,1]^n, then radially scale down any sample whose sum exceeds Gamma, as described in Section 5.3. """ z = rng.uniform(0.0, 1.0, size=(n_samples, n)) sums = z.sum(axis=1, keepdims=True) scale = np.minimum(1.0, Gamma / np.maximum(sums, 1e-12)) return z * scale def certify_inset_feasibility(rng, x, Gamma, n_samples=N_INSET_SAMPLES): """Corollary 4.5 check: budget-feasible disturbances must never exceed capacity.""" max_violation = np.zeros(m) for rr in range(m): z = sample_inset_disturbances(rng, Gamma, n_samples) realized = A_BAR[rr] @ x + z @ (D[rr] * x) max_violation[rr] = np.max(realized - B_CAP[rr]) return max_violation def stress_test(rng, x, n_samples=N_STRESS_SAMPLES): """Diagnostic out-of-set test: z_ri ~ Beta(2,5). Returns violation rate (%).""" z = rng.beta(2, 5, size=(n_samples, m, n)) a_part = A_BAR @ x d_part = D * x[None, :] sample_lhs = a_part[None, :] + np.einsum('smi,mi->sm', z, d_part) violation = np.any(sample_lhs > B_CAP[None, :], axis=1) return 100.0 * violation.mean() # --------------------------------------------------------------------------- # 4. Main experiment # --------------------------------------------------------------------------- def main(): here = os.path.dirname(os.path.abspath(__file__)) root = os.path.dirname(here) # one level up from code/ results_dir = os.path.join(root, "results") figures_dir = os.path.join(root, "figures") os.makedirs(results_dir, exist_ok=True) os.makedirs(figures_dir, exist_ok=True) eigvals = np.linalg.eigvalsh(Q) print(f"Eigenvalues of Q: {np.round(eigvals, 6)}") print(f"Smallest eigenvalue (strong-monotonicity constant mu): {eigvals.min():.3f}") # -- Equilibrium path over the full Gamma grid -------------------------- # Each Gamma is solved independently from the same fixed starting point; # this is slightly slower than warm-starting but more numerically robust # for SLSQP across the full grid (warm-starting occasionally triggered a # "positive directional derivative" line-search failure for this model). X_path = [solve_robust_equilibrium(Gamma) for Gamma in GAMMA_GRID] X_path = np.array(X_path) Phi_path = np.array([potential(x) for x in X_path]) Pi_path = Phi_path - Phi_path[0] Total_path = X_path.sum(axis=1) # -- In-set feasibility certification (Corollary 4.5) -------------------- rng_certify = np.random.default_rng(1) max_violations = np.array([certify_inset_feasibility(rng_certify, x, Gamma) for x, Gamma in zip(X_path, GAMMA_GRID)]) print(f"Max in-set violation over all Gamma/resources: {max_violations.max():.3e} " f"(tolerance {VIOLATION_TOL:.0e})") assert max_violations.max() <= VIOLATION_TOL, "In-set feasibility certification failed!" # -- Out-of-set diagnostic stress test (fixed seed) ----------------------- rng_stress = np.random.default_rng(STRESS_SEED) stress_rates = np.array([stress_test(rng_stress, x) for x in X_path]) # -- results/equilibria.csv ----------------------------------------------- eq_path = os.path.join(results_dir, "equilibria.csv") with open(eq_path, "w") as f: f.write("Gamma," + ",".join(f"x{i + 1}" for i in range(n)) + ",Total,Phi,Pi\n") for Gamma, x, total, Phi, Pi in zip(GAMMA_GRID, X_path, Total_path, Phi_path, Pi_path): f.write(f"{Gamma:.2f}," + ",".join(f"{xi:.6f}" for xi in x) + f",{total:.6f},{Phi:.6f},{Pi:.6f}\n") print(f"Wrote {eq_path}") # -- results/stress_test.csv ----------------------------------------------- stress_path = os.path.join(results_dir, "stress_test.csv") with open(stress_path, "w") as f: f.write("Gamma,Total_allocation,Robustness_premium,Violation_rate_pct\n") for Gamma, total, Pi, rate in zip(GAMMA_GRID, Total_path, Pi_path, stress_rates): f.write(f"{Gamma:.2f},{total:.6f},{Pi:.6f},{rate:.6f}\n") print(f"Wrote {stress_path}") # -- results/model_parameters.json ------------------------------------- params = { "n_players": n, "n_resources": m, "Q": Q.tolist(), "r": r.tolist(), "A_bar": A_BAR.tolist(), "D": D.tolist(), "b": B_CAP.tolist(), "x_upper_bound": X_UB, "gamma_grid": GAMMA_GRID.tolist(), "stress_seed": STRESS_SEED, "n_inset_samples": N_INSET_SAMPLES, "n_stress_samples": N_STRESS_SAMPLES, "smallest_eigenvalue_Q": float(eigvals.min()), } params_path = os.path.join(results_dir, "model_parameters.json") with open(params_path, "w") as f: json.dump(params, f, indent=2) print(f"Wrote {params_path}") # -- Console summary matching Table 1 / Table 2 -------------------------- print("\nTable 1 -- Robust equilibrium path (selected Gamma):") header = f"{'Gamma':>6}" + "".join(f"{'x' + str(i + 1):>9}" for i in range(n)) \ + f"{'Total':>9}{'Phi':>10}{'Pi':>8}" print(header) for Gamma in TABLE_GAMMAS: idx = int(np.argmin(np.abs(GAMMA_GRID - Gamma))) x = X_path[idx] print(f"{Gamma:6.1f}" + "".join(f"{xi:9.3f}" for xi in x) + f"{Total_path[idx]:9.3f}{Phi_path[idx]:10.3f}{Pi_path[idx]:8.3f}") print("\nTable 2 -- Out-of-set stress-test results (selected Gamma):") print(f"{'Gamma':>6}{'Total':>10}{'Premium':>10}{'Violation%':>12}") for Gamma in TABLE_GAMMAS: idx = int(np.argmin(np.abs(GAMMA_GRID - Gamma))) print(f"{Gamma:6.1f}{Total_path[idx]:10.3f}{Pi_path[idx]:10.3f}" f"{stress_rates[idx]:12.3f}") # -- Figures ---------------------------------------------------------------- plt.rcParams.update({ "font.size": 11, "axes.labelsize": 11, "legend.fontsize": 9.5, "xtick.labelsize": 10, "ytick.labelsize": 10, "figure.autolayout": True, }) fig, ax = plt.subplots(figsize=(5.6, 3.6)) colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd'] markers = ['o', 's', '^', 'v', 'D'] for i in range(n): ax.plot(GAMMA_GRID, X_path[:, i], marker=markers[i], ms=4, lw=1.4, color=colors[i], label=fr'$x_{i + 1}$') ax.set_xlabel(r'Uncertainty budget $\Gamma$') ax.set_ylabel('Equilibrium allocation') ax.legend(ncol=n, loc='upper center', bbox_to_anchor=(0.5, 1.18), frameon=False) ax.grid(alpha=0.25) for ext in ("pdf", "png"): fig.savefig(os.path.join(figures_dir, f"fig1_equilibrium_path.{ext}"), dpi=300) plt.close(fig) fig, ax = plt.subplots(figsize=(5.6, 3.6)) ax.plot(GAMMA_GRID, Pi_path, marker='o', ms=4, lw=1.6, color='#1f77b4') ax.set_xlabel(r'Uncertainty budget $\Gamma$') ax.set_ylabel(r'Robustness premium $\Phi(x^{\Gamma})-\Phi(x^{0})$') ax.grid(alpha=0.25) for ext in ("pdf", "png"): fig.savefig(os.path.join(figures_dir, f"fig2_robustness_premium.{ext}"), dpi=300) plt.close(fig) fig, ax = plt.subplots(figsize=(5.6, 3.6)) ax.plot(GAMMA_GRID, stress_rates, marker='o', ms=4, lw=1.6, color='#1f77b4') ax.set_xlabel(r'Uncertainty budget $\Gamma$') ax.set_ylabel('Stress-test violation rate (%)') ax.grid(alpha=0.25) for ext in ("pdf", "png"): fig.savefig(os.path.join(figures_dir, f"fig3_stress_test.{ext}"), dpi=300) plt.close(fig) print(f"\nWrote {figures_dir}/fig1_equilibrium_path.{{pdf,png}}") print(f"Wrote {figures_dir}/fig2_robustness_premium.{{pdf,png}}") print(f"Wrote {figures_dir}/fig3_stress_test.{{pdf,png}}") if __name__ == "__main__": main()