20250412-代码改进-拟蒙特卡洛

---

前言

马尔科夫链蒙特卡罗(MCMC)由于其简单有效性而得到了广泛的应用。在MCMC 中,需要用到独立同分布的随机数样本。经验表明，将拟随机数点列应用 到蒙特卡罗模拟(MC)的拟蒙特卡罗模拟（QMC）能够得到更好的模拟结果，提高 模拟的精度。

资料引用：

拟蒙特卡洛方法（Quasi-Monte Carlo, QMC）
scipy.stats手册

---

一、拟蒙特卡洛（Halton 低差异序列）

代码分析

    # 创建独立的 Halton 生成器（确保 depot 和 node 的序列不重叠）
    depot_halton = qmc.Halton(d=2, scramble=True)
    node_halton = qmc.Halton(d=2, scramble=True)

• 创建独立的 Halton 生成器（确保 depot 和 node 的序列不重叠）
• 使用 scipy 的 qmc 模块生成高度均匀的确定性序列 ，通过 Scrambling 引入随机性
• 通过 scramble=True 对 Halton 序列进行随机置乱，避免高维时的周期性相关性，同时允许统计误差估计。

depot_xy = torch.tensor(
    depot_halton.random(n=batch_size), 
    dtype=torch.float32  # 确保数据类型匹配模型预期
).reshape(-1, 1, 2)

node_xy = torch.tensor(
    node_halton.random(n=batch_size * problem_size),
    dtype=torch.float32
).reshape(batch_size, problem_size, 2)

• 生成depot_xy 和node_xy ，注：为向量PyTorch Tensor

def generate_node_demand(batch_size, problem_size):
    """
    生成节点需求 (node_demand)，形状为 (batch_size, problem_size)
    规则:
        - problem_size=20 → demand_scaler=30
        - problem_size=50 → demand_scaler=40
        - problem_size=100 → demand_scaler=50
    整数范围: [1, 9]（缩放前）
    输出值范围: [1/demand_scaler, 9/demand_scaler]
    """
    # 根据问题规模选择缩放因子
    if problem_size == 20:
        demand_scaler = 30
    elif problem_size == 50:
        demand_scaler = 40
    elif problem_size == 100:
        demand_scaler = 50
    else:
        raise NotImplementedError(f"Unsupported problem_size: {problem_size}")

    # 生成整数 [1, 9] 并缩放 (high=10 确保生成最大值为9)
    node_demand = torch.randint(
        low=1, 
        high=10,  # [1, 10) → 实际生成 1~9
        size=(batch_size, problem_size),
        dtype=torch.float32
    ) / demand_scaler

    return node_demand

• 生成符合形状和规则的 node_demand
• 添加 dtype=torch.float32 直接生成浮点张量，避免后续除法时的隐式类型转换。

二、手动实现 Halton 序列

def van_der_corput(n, base):
    vdc, denom = 0.0, 1.0
    while n != 0:
        denom *= base
        remainder = n % base
        n = n // base
        vdc += remainder / denom
    return vdc

# 生成基为 2 的前 5 个点
for i in range(1, 6):
    print(van_der_corput(i, 2))

• 仅为思路具体待实现

---

代码（一、拟蒙特卡洛（Halton 低差异序列））

# 使用拟蒙特卡洛（Halton 低差异序列）生成随机问题
import torch
import numpy as np
from scipy.stats import qmc


def get_random_problems(batch_size, problem_size):

    # 生成 Halton 序列
    # 使用 scipy 的 qmc 模块生成高度均匀的确定性序列 ，通过 Scrambling 引入随机性
    # 通过 scramble=True 对 Halton 序列进行随机置乱，避免高维时的周期性相关性，同时允许统计误差估计。
    # depot_xy_halton = qmc.Halton(d=2, scramble=True)  # scramble=True 添加随机扰码
    # node_xy_halton = qmc.Halton(d=2, scramble=True)
    # # depot_xy的 shape: (batch, 1, 2)
    # depot_xy = halton.random(size=(batch_size, 1, 2))
   
    # # node_xy的 shape: (batch, problem, 2)
    # node_xy = halton.random(size=(batch_size, problem_size, 2))
    # 创建独立的 Halton 生成器（确保 depot 和 node 的序列不重叠）
    depot_halton = qmc.Halton(d=2, scramble=True)
    node_halton = qmc.Halton(d=2, scramble=True)

    # 生成 depot_xy: (batch_size, 1, 2)
    # 生成并转换数据为 Tensor
    depot_xy = torch.tensor(
        depot_halton.random(n=batch_size),
        dtype=torch.float32
    ).reshape(-1, 1, 2)
    

    # 生成 node_xy: (batch_size, problem_size, 2)
    node_xy = torch.tensor(
        node_halton.random(n=batch_size * problem_size),
        dtype=torch.float32
    ).reshape(batch_size, problem_size, 2)    

    if problem_size == 20:
        demand_scaler = 30
    elif problem_size == 50:
        demand_scaler = 40
    elif problem_size == 100:
        demand_scaler = 50
    else:
        raise NotImplementedError

    node_demand = generate_node_demand(batch_size, problem_size)
    # node_demand的 shape: (batch, problem)

    return depot_xy, node_xy, node_demand


def augment_xy_data_by_8_fold(xy_data):
    # xy_data.shape: (batch, N, 2)

    x = xy_data[:, :, [0]]
    y = xy_data[:, :, [1]]
    # x,y shape: (batch, N, 1)

    dat1 = torch.cat((x, y), dim=2)
    dat2 = torch.cat((1 - x, y), dim=2)
    dat3 = torch.cat((x, 1 - y), dim=2)
    dat4 = torch.cat((1 - x, 1 - y), dim=2)
    dat5 = torch.cat((y, x), dim=2)
    dat6 = torch.cat((1 - y, x), dim=2)
    dat7 = torch.cat((y, 1 - x), dim=2)
    dat8 = torch.cat((1 - y, 1 - x), dim=2)

    aug_xy_data = torch.cat((dat1, dat2, dat3, dat4, dat5, dat6, dat7, dat8), dim=0)
    # shape: (8*batch, N, 2)

    return aug_xy_data



def generate_node_demand(batch_size, problem_size):
    """
    生成节点需求 (node_demand)，形状为 (batch_size, problem_size)
    规则:
        - problem_size=20 → demand_scaler=30
        - problem_size=50 → demand_scaler=40
        - problem_size=100 → demand_scaler=50
    整数范围: [1, 9]（缩放前）
    输出值范围: [1/demand_scaler, 9/demand_scaler]
    """
    # 根据问题规模选择缩放因子
    if problem_size == 20:
        demand_scaler = 30
    elif problem_size == 50:
        demand_scaler = 40
    elif problem_size == 100:
        demand_scaler = 50
    else:
        raise NotImplementedError(f"Unsupported problem_size: {problem_size}")

    # 生成整数 [1, 9] 并缩放 (high=10 确保生成最大值为9)
    node_demand = torch.randint(
        low=1, 
        high=10,  # [1, 10) → 实际生成 1~9
        size=(batch_size, problem_size),
        dtype=torch.float32
    ) / demand_scaler

    return node_demand

# # 测试代码-调试
# depot_xy, node_xy, node_demand = get_random_problems(2,20)
# print('depot_xy',depot_xy)
# print('depot_xy.shape',depot_xy.shape)

# print('node_xy',node_xy)
# print('node_xy.shape',node_xy.shape)

# print('node_demand',node_demand)
# print('node_demand.shape',node_demand.shape)