[![INFORMS Journal on Computing Logo](https://INFORMSJoC.github.io/logos/INFORMS_Journal_on_Computing_Header.jpg)](https://pubsonline.informs.org/journal/ijoc)

# Machine Learning-Assisted Stochastic Kriging for Offline Simulation Online Application

This archive is distributed in association with the [INFORMS Journal on
Computing](https://pubsonline.informs.org/journal/ijoc) under the [MIT License](LICENSE).

The software and data in this repository are a snapshot of the software and data
that were used in the research reported on in the paper 
[Machine Learning-Assisted Stochastic Kriging for Offline Simulation Online Application](https://doi.org/10.1287/ijoc.2023.0130) by Haoting Zhang, Jinghai He, Rhonda Righter, Zuo-Jun Max Shen, Zeyu Zheng.
The snapshot is based on 
[this SHA](https://github.com/INFORMSJoC/2023.0130/commit/main) 
in the development repository. 


## Cite

To cite the contents of this repository, please cite both the paper and this repo, using their respective DOIs.

https://doi.org/10.1287/ijoc.2023.0130

https://doi.org/10.1287/ijoc.2023.0130.cd

Below is the BibTex for citing this snapshot of the repository.

```
@misc{ML-SK-OSOA,
  author =        {Haoting Zhang and Jinghai He and Rhonda Righter and Zuo-Jun Max Shen and Zeyu Zheng},
  publisher =     {INFORMS Journal on Computing},
  title =         {{Machine Learning-Assisted Stochastic Kriging for Offline Simulation Online Application}},
  year =          {2024},
  doi =           {10.1287/ijoc.2023.0130.cd},
  url =           {https://github.com/INFORMSJoC/2023.0130},
  note =          {Available for download at https://github.com/INFORMSJoC/2023.0130},
}  
```

## Description

The goal of this software is to implement the Machine Learning-Assisted Stochastic Kriging (ML-SK) approach for Offline Simulation Online Application (OSOA) tasks. The ML-SK approach addresses the challenge where expensive simulation experiments must be conducted offline to learn system behavior, and then fast predictions are needed online for newly observed high-dimensional covariates.

The implementation includes:
- Core ML-SK algorithms with Linear Regression, Polynomial Regression, and Neural Network models
- Sequential and Joint training procedures as described in the paper
- Uncertainty quantification through bootstrap methods
- Three experimental scenarios: Modified Ackley function, Multi-server queueing, and Assortment optimization

All instances used in the paper are generated by the scripts provided.

## Prerequisites

To run the ML-SK implementation, you need Python 3.7+ with the following packages:

```
pip install numpy scipy torch scikit-learn matplotlib
```

For GPU acceleration (optional), ensure CUDA is properly installed for PyTorch.

## Usage

### Quick Start

To run the basic ML-SK example:

```python
python example_use_case.py
```

This demonstrates a complete supply chain optimization workflow using ML-SK.

### Running Individual Experiments

To replicate the results from the paper experiments:

**Experiment 1: Modified Ackley Function**
```python
python experiment_1_ackley.py
```

**Experiment 2: Multi-server Queueing System**
```python
python experiment_2_queue.py
```

**Experiment 3: Assortment Optimization**
```python
python experiment_3_assortment.py
```

### Comprehensive Testing

To run the full test suite validating all algorithms:

```python
python test_ml_sk_complete.py
```

## Core Implementation

The main ML-SK implementation is in `ml_sk_implementation.py` with the following key components:

### MLSKModel Class

```python
from ml_sk_implementation import MLSKModel, sparse_grid_design

# Create model with choice of ML type and kernel
model = MLSKModel(
    design_points_x=x_design,
    ml_type='neural_network',  # 'linear', 'polynomial', or 'neural_network'
    kernel_type='gaussian',    # 'gaussian' or 'matern'
    basis_degree=1
)

# Train using sequential or joint methods
model.train_sequential(covariate_samples, simulation_outputs, simulation_variances)
# OR
model.train_joint(covariate_samples, simulation_outputs, simulation_variances, alpha=1.0)

# Make fast online predictions
prediction, uncertainty = model.predict(x_test, theta_test)
```

### Device Selection

For neural network models, the system automatically detects GPU availability and prompts for CPU/GPU selection:

```
CUDA is available. Use GPU? (y/n, default=y): y
Using GPU: NVIDIA GeForce RTX 3090
```

## Results

The implementation demonstrates the effectiveness of ML-SK across three experimental scenarios:

### Experiment 1: Modified Ackley Function Results
- **Problem**: 2D decision variables, 6D covariates, closed-form test function
- **Key Finding**: Neural network ML-SK achieves lowest RMSE with efficient online prediction
- **Computational Advantage**: O(K³) vs O(N²K²) complexity reduction

### Experiment 2: Multi-server Queueing Results  
- **Problem**: Complex stochastic simulation with 31D covariates
- **Key Finding**: ML-SK handles high-dimensional covariates where classical SK fails
- **Performance**: Maintains accuracy while scaling to realistic problem sizes

### Experiment 3: Assortment Optimization Results
- **Problem**: Discrete choice modeling with market-dependent preferences  
- **Key Finding**: Effective encoding of discrete decisions in continuous ML-SK framework
- **Business Impact**: Adaptive inventory strategies based on market conditions

## Replicating Paper Results

To replicate the main experimental results from the paper:

**All Experiments with Performance Comparison:**
```bash
python test_ml_sk_complete.py
```

**Individual Experiment Replication:**
```bash
# Ackley function experiment (Figure 1 equivalent)
python experiment_1_ackley.py

# Queueing experiment (Figure 2 equivalent)  
python experiment_2_queue.py

# Assortment experiment (Figure 3 equivalent)
python experiment_3_assortment.py
```

**Supply Chain Use Case (Practical Application):**
```bash
python example_use_case.py
```

## Key Algorithmic Features

### Training Methods
- **Sequential Training**: Equations 156-161 from paper
- **Joint Training**: Equation 128 from paper with hyperparameter α

### Machine Learning Models
- **Linear Regression**: Fast, interpretable baseline
- **Polynomial Regression**: Nonlinear relationships with feature expansion
- **Neural Networks**: Maximum flexibility with GPU acceleration

### Kernel Functions
- **Gaussian Kernel**: R_M(x-x'; γ) = exp(-γ||x-x'||²)
- **Matérn Kernel**: Full implementation of Equation 27

### Uncertainty Quantification
- **Bootstrap Method**: Algorithm 1 from paper
- **Credible Intervals**: Accounts for both learning and kriging uncertainty




