diff --git a/examples/ep_load_balancer/README.md b/examples/ep_load_balancer/README.md
index aedc79943a81..eaa277c4a180 100644
--- a/examples/ep_load_balancer/README.md
+++ b/examples/ep_load_balancer/README.md
@@ -2,17 +2,16 @@
 
 Effective load balancing is crucial when leveraging large-scale expert parallelism. As described in the [DeepSeek-V3 paper](https://arxiv.org/abs/2412.19437), redundant experts can be introduced to rebalance the workload across GPUs. This mechanism is known as the Expert Parallelism Load Balancer ([EPLB](https://github.com/deepseek-ai/EPLB)).
 
-> **Note:** Currently, only the offline EP load balancer is supported.
-
 ## Offline EP Load Balancer
 
 ### Step 1: Run Inference and Collect Statistics
 
-To generate the necessary statistics for load balancing, run your model on a target dataset (e.g., GSM8K) while counting the routed expert IDs during inference. Once counting is complete, the statistics will be saved for further processing.
+To generate the necessary statistics for load rebalancing, run your model on a target dataset and count the routed expert IDs during inference. Once the counting is complete, the statistics will be saved for further processing. In this example, we use `deepseek-ai/DeepSeek-R1`.
 
 Set up some environment variables:
 
 ```bash
+export MODEL_NAME=deepseek-ai/DeepSeek-R1
 export MODEL_PATH=<YOUR_MODEL_PATH>
 # Set the expert statistic data path
 export EXPERT_STATISTIC_PATH=./expert_statistic
@@ -20,55 +19,191 @@ export EXPERT_STATISTIC_PATH=./expert_statistic
 export EXPERT_STATISTIC_ITER_RANGE=100-200
 ```
 
-Prepare a configuration file and run inference on GSM8K:
+Prepare a dataset following the [benchmarking documentation](https://github.com/NVIDIA/TensorRT-LLM/blob/main/docs/source/performance/perf-benchmarking.md#preparing-a-dataset) and save it as `./dataset.json`.
+
+Run 32-way expert parallelism inference on the prepared dataset. Please refer to the [LLM API MGMN example](https://github.com/NVIDIA/TensorRT-LLM/blob/main/examples/llm-api/llm_mgmn_trtllm_bench.sh) for details on running `trtllm-bench` on Slurm.
 
 ```bash
 cat > ./extra_llm_api_options.yaml <<EOF
 enable_attention_dp: true
 use_cuda_graph: true
+moe_max_num_tokens: 8192
 EOF
 
-trtllm-eval --model $MODEL_PATH \
-    --tp_size 8 \
-    --ep_size 8 \
+trtllm-llmapi-launch \
+trtllm-bench --model ${MODEL_NAME} \
+    --model_path ${MODEL_PATH} \
+    throughput \
+    --tp 32 \
+    --ep 32 \
     --extra_llm_api_options ./extra_llm_api_options.yaml \
-    --backend pytorch gsm8k
+    --kv_cache_free_gpu_mem_fraction 0.75 \
+    --backend pytorch \
+    --dataset ./dataset.json \
+    --warmup 0 \
+    --eos_id -1
+```
+
+After inference, review the dumped statistic files in `$EXPERT_STATISTIC_PATH`. For each layer and iteration, the load imbalance can be measured using simple metrics such as the standard deviation or the imbalance ratio. Given the routed token counts for all ranks, the imbalance ratio is defined as $(max - mean) / mean$, which represents the excessive workload received by the hottest rank. A perfectly balanced load would have an imbalance ratio of 0. Run the [`report_load_statistics.py`](./report_load_statistics.py) script:
+
+```bash
+python report_load_statistics.py --expert_statistic_path $EXPERT_STATISTIC_PATH
 ```
 
-After inference, review the dumped statistic files in `$EXPERT_STATISTIC_PATH`.
+The output would look like:
+
+```txt
+Load statistics:
+           mean         std  imbalance-ratio
+3        1024.0  187.955200         0.498043
+4        1024.0  202.728516         0.537602
+5        1024.0  209.339981         0.458676
+...
+58       1024.0  570.880676         2.461014
+59       1024.0  341.339447         0.717498
+60       1024.0  381.045471         1.119648
+average  1024.0  491.651199         1.564272
+```
+
+The metrics are computed at each iteration and then averaged. The load imbalance is significant — on average, the hottest rank receives 1.56 times more routed tokens than the mean.
 
 ### Step 2: Generate the EPLB Configuration
 
-Use the provided [`generate_eplb_config.py`](./generate_eplb_config.py) script to convert the collected statistics into an EPLB configuration file. Specify the expert parallelism size (`--ep_size`) and the total number of slots (`--num_slots`) that will be used for deployment:
+Use the provided [`generate_eplb_config.py`](./generate_eplb_config.py) script to convert the collected statistics into an EPLB configuration file. Specify the target expert parallelism size (`--ep_size`) and the total number of slots (`--num_slots`) that will be used for deployment. One potential strategy is to maintain 8 expert slots per rank while increasing expert parallelism to 36 ways. This results in 32 redundant experts and 288 expert slots in total.
 
 ```bash
 python generate_eplb_config.py \
-    --ep_size 8 \
-    --num_slots 320 \
+    --ep_size 36 \
+    --num_slots 288 \
     --expert_statistic_path $EXPERT_STATISTIC_PATH \
     --output_path ./moe_load_balancer.yaml
 ```
 
+The `./moe_load_balancer.yaml` file would look like:
+
+```yaml
+initial_global_assignments:
+  3: [138, 81, 60, ..., 69, 250, 77]
+  4: [24, 243, 72, ..., 90, 251, 52]
+  5: [120, 162, 246, ..., 14, 192, 171]
+  ...
+  58: [67, 70, 160, ..., 212, 103, 125]
+  59: [45, 142, 152, ..., 99, 205, 49]
+  60: [34, 162, 119, ..., 234, 26, 129]
+num_slots: 288
+layer_updates_per_iter: 0
+```
+
+`layer_updates_per_iter` is the number of layers of which the MoE weights are updated per iteration; `layer_updates_per_iter` of 0 means MoE weights are not updated during inference, so it is static EP Load Balancer.
+
+`initial_global_assignments` is a dict that maps MoE layer index to a list of length 288 (`num_slots`); at layer `i`, the `j`-th expert slot is assigned with expert ID `initial_global_assignments[i][j]`. For each layer, every successive 8 expert slots are assigned to a rank.
+
+
 ### Step 3: Run Inference with the EPLB Configuration
 
-Disable the expert ID counting by unsetting the environment variable:
+Set up some environment variables:
+
+```bash
+# Set a new expert statistic data path
+export EXPERT_STATISTIC_PATH=./expert_statistic_eplb
+# Enable counting of routed expert IDs from iteration 100 to iteration 200
+export EXPERT_STATISTIC_ITER_RANGE=100-200
+```
+
+Run 36-way expert parallelism inference with the EPLB configuration incorporated:
+
+```bash
+cat > ./extra_llm_api_options_eplb.yaml <<EOF
+enable_attention_dp: true
+use_cuda_graph: true
+moe_max_num_tokens: 9216
+moe_load_balancer: ./moe_load_balancer.yaml
+EOF
+
+trtllm-llmapi-launch \
+trtllm-bench --model ${MODEL_NAME} \
+    --model_path ${MODEL_PATH} \
+    throughput \
+    --tp 36 \
+    --ep 36 \
+    --extra_llm_api_options ./extra_llm_api_options_eplb.yaml \
+    --kv_cache_free_gpu_mem_fraction 0.75 \
+    --backend pytorch \
+    --dataset ./dataset.json \
+    --warmup 0 \
+    --eos_id -1
+```
+
+Run the [`report_load_statistics.py`](./report_load_statistics.py) script again:
 
 ```bash
-unset EXPERT_STATISTIC_ITER_RANGE
+python report_load_statistics.py --expert_statistic_path $EXPERT_STATISTIC_PATH
+```
+
+The output would look like:
+
+```txt
+Load statistics:
+           mean        std  imbalance-ratio
+3        1024.0  37.612328         0.081947
+4        1024.0  42.367714         0.093256
+5        1024.0  42.623219         0.092623
+...
+58       1024.0  49.167507         0.113420
+59       1024.0  44.529514         0.092314
+60       1024.0  48.408348         0.101029
+average  1024.0  53.976442         0.115378
 ```
 
-Prepare a new configuration file that incorporates the generated EPLB configuration, then run inference on GSM8K:
+Clearly, the load is much more balanced now — on average, the hottest rank receives only about 0.11 times more routed tokens than the mean.
+
+> **Note:** The expert ID counting could significantly hurt performance, so remember to disable it by unsetting `EXPERT_STATISTIC_ITER_RANGE` when running inference for benchmarking or production purposes.
+
+
+## Online EP Load Balancer
+
+Online EP Load Balancer is more suitable for production deployment needs to react timely to the online traffic changes. We still use 8 expert slots per rank and 36-way expert parallelism.
+
+Prepare the EPLB configuration file:
+
+```bash
+cat > ./moe_load_balancer.yaml <<EOF
+num_slots: 288
+layer_updates_per_iter: 2
+EOF
+```
+
+`layer_updates_per_iter` of 2 means that at each iteration, the MoE weights of 2 layers are updated dynamically according to the online statistics. Different from offline EP Load Balancer, `initial_global_assignments` is not important anymore, since the expert assignments will be properly and regularly updated during the inference. Hence, `initial_global_assignments` can be omitted in the configuration.
+
+Run 36-way expert parallelism inference with the EPLB configuration incorporated:
 
 ```bash
 cat > ./extra_llm_api_options_eplb.yaml <<EOF
 enable_attention_dp: true
 use_cuda_graph: true
+moe_max_num_tokens: 9216
 moe_load_balancer: ./moe_load_balancer.yaml
 EOF
 
-trtllm-eval --model $MODEL_PATH \
-    --tp_size 8 \
-    --ep_size 8 \
+trtllm-llmapi-launch \
+trtllm-bench --model ${MODEL_NAME} \
+    --model_path ${MODEL_PATH} \
+    throughput \
+    --tp 36 \
+    --ep 36 \
     --extra_llm_api_options ./extra_llm_api_options_eplb.yaml \
-    --backend pytorch gsm8k
+    --kv_cache_free_gpu_mem_fraction 0.75 \
+    --backend pytorch \
+    --dataset ./dataset.json \
+    --warmup 0 \
+    --eos_id -1
 ```
+
+> **Note:** Similar to offline EP Load Balancer, you can enable expert ID counting to verify the effectiveness of EPLB, but remember to disable it when running inference for benchmarking or production purposes.
+
+> **Explanation on moe_max_num_tokens:** For Large Scale EP, there can be extreme conditions that all ranks send tokens to a single rank since they all want that expert.
+In that case, that rank will have too many tokens to compute. In order not to make the hot rank OOM, there is one strategy that chunk the tokens if there are too much.
+`moe_max_num_tokens` is the parameter that controls the max chunk size. However, this may have performance penalty if there is enough since batch size is smaller.
+So by default, it is set to some value that all tokens can complete in one wave. However, if EP size is large, we may need to trade off that in order not to OOM or got other runtime errors due to lack of memory.
+One good point is that if memory is OK, we can set `moe_max_num_tokens` to `max_batch_size * ep_size` to make all generation requests can be processed in one chunk.
+For example, if `ep_size` is 36 and `max_batch_size` is 256, we may set `moe_max_num_tokens` to 9216.
diff --git a/examples/ep_load_balancer/generate_eplb_config.py b/examples/ep_load_balancer/generate_eplb_config.py
index 22d982e2aaa0..c7b110dde666 100644
--- a/examples/ep_load_balancer/generate_eplb_config.py
+++ b/examples/ep_load_balancer/generate_eplb_config.py
@@ -1,46 +1,13 @@
 import argparse
-import glob
-import json
-import os
 
-import pandas as pd
-import safetensors
 import torch
 import yaml
+from utils import load_expert_statistic
 
 from tensorrt_llm.bindings.internal.runtime import (MoeLoadBalanceMetaInfo,
                                                     MoePlacementCpuInfo,
                                                     do_placement,
                                                     do_replication)
-from tensorrt_llm.logger import logger
-
-logger.set_level("info")
-
-
-def calculate_load_statistics(load_iters: torch.Tensor):
-    # sum the loads over iterations and calculate the statistics
-    load_total = load_iters.sum(dim=0)
-    mean = load_total.mean().item()
-    std = load_total.std().item()
-    imbalance_ratio = (load_total.max().item() - load_total.min().item()) / mean
-    stats = {
-        "mean-total": mean,
-        "std-total": std,
-        "imbalance-ratio-total": imbalance_ratio
-    }
-
-    # calculate the statistics for each iteration and average over iterations
-    mean = load_iters.mean(dim=-1).mean().item()
-    std = load_iters.std(dim=-1).mean().item()
-    imbalance_ratio = (load_iters.max(dim=-1).values -
-                       load_iters.min(dim=-1).values) / load_iters.mean(dim=-1)
-    imbalance_ratio = imbalance_ratio.mean().item()
-    stats.update({
-        "mean-iter": mean,
-        "std-iter": std,
-        "imbalance-ratio-iter": imbalance_ratio
-    })
-    return stats
 
 
 def save_eplb_config(config: dict, path: str):
@@ -58,116 +25,84 @@ def represent_list_inline(dumper, data):
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser()
-    parser.add_argument("--expert_statistic_path",
-                        type=str,
-                        default=os.environ.get("EXPERT_STATISTIC_PATH",
-                                               "expert_statistic"))
-    parser.add_argument("--iter_start", type=int, default=None)
-    parser.add_argument("--iter_stop", type=int, default=None)
+    parser.add_argument(
+        "--expert_statistic_path",
+        type=str,
+        required=True,
+        help="The directory path to the expert statistic files.")
+    parser.add_argument("--iter_start",
+                        type=int,
+                        default=None,
+                        help="The start iteration of used iterations.")
+    parser.add_argument("--iter_stop",
+                        type=int,
+                        default=None,
+                        help="The end iteration of used iterations.")
     parser.add_argument("--output_path",
                         type=str,
-                        default="moe_load_balancer.yaml")
-    parser.add_argument("--ep_size", type=int, default=8)
-    parser.add_argument("--num_slots", type=int, default=320)
-    parser.add_argument("--layer_updates_per_iter", type=int, default=0)
+                        required=True,
+                        help="The output path to the eplb config file.")
+    parser.add_argument(
+        "--ep_size",
+        type=int,
+        default=None,
+        help="The expert parallelism size after load rebalance.")
+    parser.add_argument(
+        "--num_slots",
+        type=int,
+        default=None,
+        help="The total number of expert slots after load rebalance.")
+    parser.add_argument("--layer_updates_per_iter",
+                        type=int,
+                        default=0,
+                        help="The number of layers to update per iteration.")
     args = parser.parse_args()
 
-    with open(f"{args.expert_statistic_path}/meta_info.json", "r") as f:
-        meta_info = json.load(f)
+    meta_info, statistic = load_expert_statistic(args.expert_statistic_path)
     num_experts = meta_info["num_experts"]
     num_experts_per_token = meta_info["num_experts_per_token"]
 
-    statistic = {}
-    for statistic_file in glob.glob(
-            f"{args.expert_statistic_path}/rank*.safetensors"):
-        rank_statistic = safetensors.torch.load_file(statistic_file)
-        for key, data in rank_statistic.items():
-            if key not in statistic:
-                statistic[key] = torch.zeros_like(data)
-            statistic[key] += data
-
-    def parse_key(key: str) -> tuple[int, int]:
-        iter_idx, layer_idx = key.split("_")
-        return int(iter_idx), int(layer_idx)
-
-    statistic = {parse_key(key): data for key, data in statistic.items()}
-
-    iters = sorted(list(set(iter_idx for iter_idx, _ in statistic)))
-    layers = sorted(list(set(layer_idx for _, layer_idx in statistic)))
-    num_iters = len(iters)
-    num_layers = len(layers)
-    assert iters[-1] + 1 - iters[0] == num_iters
-    assert len(statistic) == num_iters * num_layers
-
+    if args.ep_size is None:
+        args.ep_size = meta_info["ep_size"]
+    if args.num_slots is None:
+        args.num_slots = num_experts
     if args.iter_start is None:
-        args.iter_start = iters[0]
+        args.iter_start = meta_info["iter_start"]
     if args.iter_stop is None:
-        args.iter_stop = iters[-1] + 1
-    logger.info(f"Statistic iterations: {iters[0]} - {iters[-1] + 1}")
-    logger.info(f"Used iterations: {args.iter_start} - {args.iter_stop}")
-    logger.info(f"Statistic layers: {layers}")
+        args.iter_stop = meta_info["iter_stop"]
+    num_iters = args.iter_stop - args.iter_start
 
     num_local_slots = args.num_slots // args.ep_size
     initial_global_assignments = {}
-    load_stats = {}
-    load_stats_rebalanced = {}
-
-    for layer_idx in layers:
+    for layer_idx in meta_info["layers"]:
         expert_token_count_iters = [
             data for key, data in statistic.items() if
             args.iter_start <= key[0] < args.iter_stop and key[1] == layer_idx
         ]
         expert_token_count_iters = torch.stack(expert_token_count_iters, dim=0)
-        assert expert_token_count_iters.size(
-            0) == args.iter_stop - args.iter_start
+        assert expert_token_count_iters.size(0) == num_iters
         expert_load_factor = expert_token_count_iters.sum(dim=0).float()
 
-        meta_info = MoeLoadBalanceMetaInfo(expert_count=num_experts,
-                                           top_k=num_experts_per_token,
-                                           ep_rank=0,
-                                           ep_size=args.ep_size,
-                                           slot_count_per_rank=num_local_slots)
+        moelb_info = MoeLoadBalanceMetaInfo(expert_count=num_experts,
+                                            top_k=num_experts_per_token,
+                                            ep_rank=0,
+                                            ep_size=args.ep_size,
+                                            slot_count_per_rank=num_local_slots)
         placement_info = MoePlacementCpuInfo()
         placement_info.expert_replica_count = [0] * num_experts
         placement_info.rank_expert_ids = [[0] * num_local_slots
                                           for _ in range(args.ep_size)]
 
-        do_replication(meta_info, expert_load_factor.tolist(), placement_info)
-        do_placement(meta_info, expert_load_factor.tolist(), placement_info)
+        do_replication(moelb_info, expert_load_factor.tolist(), placement_info)
+        do_placement(moelb_info, expert_load_factor.tolist(), placement_info)
 
         initial_global_assignments[layer_idx] = []
         for local_expert_ids in placement_info.rank_expert_ids:
             initial_global_assignments[layer_idx].extend(local_expert_ids)
 
-        # Report load statistics
-        rank_load_iters = expert_token_count_iters.reshape(
-            expert_token_count_iters.size(0), args.ep_size, -1).sum(dim=-1)
-        load_stats[layer_idx] = calculate_load_statistics(
-            rank_load_iters.float())
-
-        token_load_iters_rebalanced = expert_token_count_iters / torch.tensor(
-            placement_info.expert_replica_count)
-        token_load_iters_rebalanced = token_load_iters_rebalanced[:,
-                                                                  torch.
-                                                                  tensor(initial_global_assignments[
-                                                                      layer_idx]
-                                                                         )]
-        rank_load_iters_rebalanced = token_load_iters_rebalanced.reshape(
-            expert_token_count_iters.size(0), args.ep_size, -1).sum(dim=-1)
-        load_stats_rebalanced[layer_idx] = calculate_load_statistics(
-            rank_load_iters_rebalanced.float())
-
     eplb_config = {
         "num_slots": args.num_slots,
         "initial_global_assignments": initial_global_assignments,
         "layer_updates_per_iter": args.layer_updates_per_iter,
     }
     save_eplb_config(eplb_config, args.output_path)
-
-    load_stats = pd.DataFrame(load_stats).T
-    logger.info(f"Load statistics:\n{load_stats}")
-
-    load_stats_rebalanced = pd.DataFrame(load_stats_rebalanced).T
-    logger.info(
-        f"Load statistics after rebalance (estimated):\n{load_stats_rebalanced}"
-    )
diff --git a/examples/ep_load_balancer/report_load_statistics.py b/examples/ep_load_balancer/report_load_statistics.py
new file mode 100644
index 000000000000..bc26597c0674
--- /dev/null
+++ b/examples/ep_load_balancer/report_load_statistics.py
@@ -0,0 +1,68 @@
+import argparse
+
+import pandas as pd
+import torch
+from utils import load_expert_statistic
+
+
+def calculate_load_statistics(load_iters: torch.Tensor):
+    # calculate the statistics for each iteration and average over iterations
+    mean = load_iters.mean(dim=-1).mean().item()
+    std = load_iters.std(dim=-1).mean().item()
+    imbalance_ratio = load_iters.max(dim=-1).values / load_iters.mean(
+        dim=-1) - 1
+    imbalance_ratio = imbalance_ratio.mean().item()
+    return {"mean": mean, "std": std, "imbalance-ratio": imbalance_ratio}
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--expert_statistic_path",
+        type=str,
+        required=True,
+        help="The directory path to the expert statistic files.")
+    parser.add_argument("--iter_start",
+                        type=int,
+                        default=None,
+                        help="The start iteration of used iterations.")
+    parser.add_argument("--iter_stop",
+                        type=int,
+                        default=None,
+                        help="The end iteration of used iterations.")
+    parser.add_argument("--per_expert",
+                        default=False,
+                        action="store_true",
+                        help="Report the load statistics per expert.")
+    args = parser.parse_args()
+
+    meta_info, statistic = load_expert_statistic(args.expert_statistic_path)
+    num_experts = meta_info["num_experts"]
+    num_experts_per_token = meta_info["num_experts_per_token"]
+
+    if args.iter_start is None:
+        args.iter_start = meta_info["iter_start"]
+    if args.iter_stop is None:
+        args.iter_stop = meta_info["iter_stop"]
+    num_iters = args.iter_stop - args.iter_start
+
+    load_stats = {}
+    for layer_idx in meta_info["layers"]:
+        expert_token_count_iters = [
+            data for key, data in statistic.items() if
+            args.iter_start <= key[0] < args.iter_stop and key[1] == layer_idx
+        ]
+        expert_token_count_iters = torch.stack(expert_token_count_iters, dim=0)
+        assert expert_token_count_iters.size(0) == num_iters
+
+        if args.per_expert:
+            load_iters = expert_token_count_iters
+        else:
+            load_iters = expert_token_count_iters.reshape(
+                num_iters, meta_info["ep_size"], -1).sum(dim=-1)
+        load_stats[layer_idx] = calculate_load_statistics(load_iters.float())
+
+    load_stats = pd.DataFrame(load_stats)
+    load_stats["average"] = load_stats.mean(axis=1)
+    load_stats = load_stats.T
+    print(f"Load statistics:\n{load_stats}")
diff --git a/examples/ep_load_balancer/utils.py b/examples/ep_load_balancer/utils.py
new file mode 100644
index 000000000000..b4c59a34978e
--- /dev/null
+++ b/examples/ep_load_balancer/utils.py
@@ -0,0 +1,37 @@
+import glob
+import json
+
+import safetensors.torch
+import torch
+
+
+def load_expert_statistic(path: str):
+    with open(f"{path}/meta_info.json", "r") as f:
+        meta_info = json.load(f)
+
+    statistic_files = glob.glob(f"{path}/rank*.safetensors")
+    statistic = {}
+    for statistic_file in statistic_files:
+        rank_statistic = safetensors.torch.load_file(statistic_file)
+        for key, data in rank_statistic.items():
+            if key not in statistic:
+                statistic[key] = torch.zeros_like(data)
+            statistic[key] += data
+
+    def parse_key(key: str) -> tuple[int, int]:
+        iter_idx, layer_idx = key.split("_")
+        return int(iter_idx), int(layer_idx)
+
+    statistic = {parse_key(key): data for key, data in statistic.items()}
+
+    iters = sorted(list(set(iter_idx for iter_idx, _ in statistic)))
+    layers = sorted(list(set(layer_idx for _, layer_idx in statistic)))
+    num_iters = len(iters)
+    num_layers = len(layers)
+    assert iters[-1] + 1 - iters[0] == num_iters
+    assert len(statistic) == num_iters * num_layers
+    meta_info["ep_size"] = len(statistic_files)
+    meta_info["iter_start"] = iters[0]
+    meta_info["iter_stop"] = iters[-1] + 1
+    meta_info["layers"] = layers
+    return meta_info, statistic