Equilibrium computation strategies can optimize the wave-based dispatch of quantum circuits in hybrid HPC-quantum systems to minimize resource consumption.

Applying Nash equilibrium or correlated equilibrium computation strategies to schedule quantum circuit dispatch in hybrid HPC-quantum systems will reduce aggregate resource consumption (qubit-hours, classical CPU cycles, and memory bandwidth) by ≥15% compared to baseline first-come-first-served (FCFS) or greedy dispatch policies, measurable under controlled workload conditions with ≥3 distinct quantum circuit types and ≥2 QPU backends, with statistical significance p < 0.05 across ≥30 independent trials.

PHASE 1 — Simulation Baseline (Days 1–15): Implement a discrete-event simulation of a hybrid HPC-quantum system using SimPy or custom event engine. Establish FCFS, round-robin, and shortest-job-first baselines across 3 workload profiles.

PHASE 2 — Equilibrium Solver Integration (Days 16–35): Implement Nash equilibrium solver (support enumeration + Lemke-Howson) and correlated equilibrium solver (linear programming) as scheduling decision engines. Integrate with wave-dispatch queue model.

PHASE 3 — Controlled Comparison (Days 36–55): Run 30 independent trials per configuration (6 dispatch policies × 4 workload intensities × 3 circuit type mixes = 720 total simulation runs). Measure qubit-hours, CPU cycles, queue wait time, and throughput.

PHASE 4 — Emulation on Real Hardware (Days 56–80): Deploy top-2 performing strategies on IBM Quantum Network or AWS Braket hybrid simulator with real QPU access. Validate simulation predictions against hardware measurements on ≥100 real circuit executions.

PHASE 5 — Statistical Analysis and Reporting (Days 81–90): Apply ANOVA + post-hoc Tukey HSD. Compute effect sizes (Cohen's d). Generate reproducibility package.

# ARCHITECTURE OVERVIEW: Equilibrium Wave Dispatcher (EWD)
# ============================================================

# --- Core Data Structures ---
@dataclass
class QuantumJob:
    job_id: str
    circuit_depth: int        # gate count
    n_qubits: int
    circuit_type: str         # 'variational', 'clifford', 'qft'
    arrival_time: float
    priority: int = 1

@dataclass
class QPUNode:
    node_id: str
    n_qubits: int
    gate_error_rate: float    # e.g., 0.001
    t1_us: float              # coherence time microseconds
    current_job: Optional[QuantumJob] = None
    available_at: float = 0.0

# --- Wave Buffer ---
class WaveBuffer:
    def __init__(self, wave_size: int = 8):
        self.wave_size = wave_size
        self.queue: List[QuantumJob] = []
    
    def is_ready(self) -> bool:
        return len(self.queue) >= self.wave_size
    
    def get_wave(self) -> List[QuantumJob]:
        wave = self.queue[:self.wave_size]
        self.queue = self.queue[self.wave_size:]
        return wave

# --- Payoff Matrix Construction ---
def build_payoff_matrix(jobs: List[QuantumJob], 
                         qpus: List[QPUNode]) -> np.ndarray:
    """
    Payoff[i][j] = negative resource cost of assigning job_class_i to QPU_j
    Resource cost = qubit_hours + error_penalty + wait_time_penalty
    """
    n_job_classes = len(set(j.circuit_type for j in jobs))
    n_qpus = len(qpus)
    payoff = np.zeros((n_job_classes, n_qpus))
    
    job_classes = list(set(j.circuit_type for j in jobs))
    for i, jc in enumerate(job_classes):
        class_jobs = [j for j in jobs if j.circuit_type == jc]
        avg_depth = np.mean([j.circuit_depth for j in class_jobs])
        avg_qubits = np.mean([j.n_qubits for j in class_jobs])
        
        for k, qpu in enumerate(qpus):
            exec_time = estimate_execution_time(avg_depth, qpu)
            error_penalty = avg_depth * qpu.gate_error_rate * 100
            qubit_hours = (avg_qubits * exec_time) / 3600.0
            wait_penalty = max(0, qpu.available_at - current_time()) * 0.1
            
            # Negative because we minimize cost (maximize negative cost)
            payoff[i][k] = -(qubit_hours + error_penalty + wait_penalty)
    
    return payoff

# --- Nash Equilibrium Solver ---
def compute_nash_equilibrium(payoff: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    """
    For 2-player case: use Nashpy support enumeration
    For N-player: use Gambit via subprocess call
    Returns: (row_strategy, col_strategy) as probability distributions
    """
    import nashpy as nash
    
    if payoff.shape[0] <= 4:  # Small game: direct solver
        game = nash.Game(payoff, -payoff.T)  # Zero-sum approximation
        equilibria = list(game.support_enumeration())
        if equilibria:
            return equilibria[0]  # Return first NE found
    
    # Fallback: correlated equilibrium via LP
    return compute_correlated_equilibrium(payoff)

def compute_correlated_equilibrium(payoff: np.ndarray) -> np.ndarray:
    """
    Solve for correlated equilibrium using CVXPY LP
    Maximize: sum(p[i,j] * payoff[i,j])
    Subject to: incentive compatibility, p >= 0, sum(p) = 1
    """
    import cvxpy as cp
    
    n_rows, n_cols = payoff.shape
    p = cp.Variable((n_rows, n_cols), nonneg=True)
    
    objective = cp.Maximize(cp.sum(cp.multiply(p, payoff)))
    constraints = [cp.sum(p) == 1]
    
    # Incentive compatibility: no player wants to deviate
    for i in range(n_rows):
        for i_prime in range(n_rows):
            if i != i_prime:
                constraints.append(
                    cp.sum(cp.multiply(p[i,:], payoff[i,:] - payoff[i_prime,:])) >= 0
                )
    
    prob = cp.Problem(objective, constraints)
    prob.solve(solver=cp.GLPK, verbose=False)
    
    return p.value if p.value is not None else np.ones((n_rows, n_cols)) / (n_rows * n_cols)

# --- Wave Dispatcher ---
class EquilibriumWaveDispatcher:
    def __init__(self, qpus: List[QPUNode], wave_size: int = 8):
        self.qpus = qpus
        self.buffer = WaveBuffer(wave_size)
        self.metrics = ResourceMetrics()
        self.solver_times = []
    
    def dispatch_wave(self, wave: List[QuantumJob]) -> Dict[str, str]:
        """Returns mapping: job_id -> qpu_id"""
        t_start = time.perf_counter()
        
        payoff = build_payoff_matrix(wave, self.qpus)
        eq_strategy = compute_nash_equilibrium(payoff)
        
        t_solve = time.perf_counter() - t_start
        self.solver_times.append(t_solve)
        
        # Sample assignments from equilibrium distribution
        assignments = {}
        job_classes = list(set(j.circuit_type for j in wave))
        
        for job in wave:
            class_idx = job_classes.index(job.circuit_type)
            row_probs = eq_strategy[0] if isinstance(eq_strategy, tuple) else eq_strategy[class_idx]
            
            # Normalize and sample QPU assignment
            col_probs = np.abs(row_probs) / np.sum(np.abs(row_probs))
            qpu_idx = np.random.choice(len(self.qpus), p=col_probs)
            assignments[job.job_id] = self.qpus[qpu_idx].node_id
        
        return assignments
    
    def run(self, job_stream: Iterator[QuantumJob]):
        for job in job_stream:
            self.buffer.queue.append(job)
            if self.buffer.is_ready():
                wave = self.buffer.get_wave()
                assignments = self.dispatch_wave(wave)
                self.execute_wave(wave, assignments)
                self.metrics.record_wave(wave, assignments, self.qpus)

# --- Metrics Collection ---
@dataclass
class ResourceMetrics:
    total_qubit_hours: float = 0.0
    total_wait_time: float = 0.0
    throughput_jobs: int = 0
    qpu_utilization: Dict[str, float] = field(default_factory=dict)
    
    def record_wave(self, wave, assignments, qpus):
        for job in wave:
            qpu = next(q for q in qpus if q.node_id == assignments[job.job_id])
            exec_time = estimate_execution_time(job.circuit_depth, qpu)
            self.total_qubit_hours += (job.n_qubits * exec_time) / 3600.0
            self.throughput_jobs += 1

# --- Simulation Entry Point ---
def run_experiment(n_qpus: int, arrival_rate: float, 
                   policy: str, n_trials: int = 30) -> pd.DataFrame:
    results = []
    for trial in range(n_trials):
        np.random.seed(trial * 137)
        qpus = [QPUNode(f"qpu_{i}", n_qubits=27, 
                        gate_error_rate=0.001, t1_us=100.0) 
                for i in range(n_qpus)]
        
        dispatcher = EquilibriumWaveDispatcher(qpus) if policy == 'equilibrium' \
                     else BaselineDispatcher(qpus, policy=policy)
        
        job_stream = PoissonJobGenerator(rate=arrival_rate, duration=3600)
        dispatcher.run(job_stream)
        
        results.append({
            'trial': trial, 'policy': policy, 'n_qpus': n_qpus,
            'arrival_rate': arrival_rate,
            'qubit_hours': dispatcher.metrics.total_qubit_hours,
            'throughput': dispatcher.metrics.throughput_jobs,
            'mean_wait': dispatcher.metrics.total_wait_time / max(1, dispatcher.metrics.throughput_jobs)
        })
    
    return pd.DataFrame(results)

# --- Statistical Analysis ---
def analyze_results(df: pd.DataFrame) -> Dict:
    from scipy import stats
    import pingouin as pg
    
    baseline = df[df['policy'] == 'fcfs']['qubit_hours']
    equilibrium = df[df['policy'] == 'equilibrium']['qubit_hours']
    
    t_stat, p_value = stats.ttest_ind(baseline, equilibrium)
    cohens_d = (baseline.mean() - equilibrium.mean()) / np.sqrt(
        (baseline.std()**2 + equilibrium.std()**2) / 2)
    pct_improvement = (baseline.mean() - equilibrium.mean()) / baseline.mean() * 100
    
    return {
        'pct_improvement': pct_improvement,
        'p_value': p_value,
        'cohens_d': cohens_d,
        'significant': p_value < 0.05 and pct_improvement >= 15.0
    }