抽样平均近似方法 Sample Average Approximation （SAA）

心态与习惯

已于 2023-09-13 18:01:22 修改

阅读量1.7w

点赞数 22

于 2019-05-18 19:23:55 首次发布

本文链接：https://blog.csdn.net/robert_chen1988/article/details/90319257

版权

数学优化同时被 2 个专栏收录

145 篇文章

订阅专栏

随机规划

11 篇文章

订阅专栏

多次听说或看到抽样平均方法 Sample Average Approximation ，但不是真正理解这个方法到底怎么用的。决定写一篇博客，将每次看到的心得体会记录下来。

对于一个随机规划问题，

$\min\quad f(x, \xi)$

其中 $\xi$ 为一个随机变量，而 $x$ 是一个求解变量。抽样平均方法 SAA 的大致含义是用抽样的方法将随机变量用样本表示，从而将随机规划问题转化为确定性问题。假设随机变量 $\xi$ 的样本分别为： $\xi_1, \xi_2,\dots, \xi_N$ . 则原问题转化为：

$\min\quad \frac{1}{N}\sum_{i=1}^{N} f(x, \xi_i)$

而在每个样本下， $\min f(x, \xi)$ 可以通过线性规划求解。

该方法有一定的适用范围，这一点还没弄清楚。

1. 单阶段时

例如一个单阶段报童模型：

$\max\quad\pi= p\min\{Q, d\}-cQ+\gamma\max\{Q-d, 0\}$

假如随机需求抽样了 $N$ 个样本，每个样本值 $d^i$ ，则表达式可以转化为：
$\max\quad\pi= \frac{1}{N}\sum_{i=1}^Np\min\{Q, d^i\}-cQ+\gamma\max\{Q-d^i, 0\}$

进一步转化为可以用规划软件求解的混合整数规划模型，首先引入两个辅助变量表示目标函数中的 min ，max 表达式：

$\begin{aligned} \max\quad&\pi= \frac{1}{N}\sum_{i=1}^N(px^i+\gamma y^i)-cQ\\ s.t.\quad&\forall i=1, 2, \dots, N\\ &x^i=\min\{Q, d^i\}\\ &y^i=\max\{Q-d^i, 0\}\\ & Q\geq 0 \end{aligned}$

可以再进一步引入辅助变量 $\delta$ 将约束条件中的 min 或 max 转化为线性表达式，不过很多规划软件已经可以求解了。

使用 python 调用 gurobi 实现单阶段报童模型的例子：

# -*- coding: utf-8 -*-
"""
Created on Tue Feb 23 15:08:06 2021

@author: zhen chen

MIT Licence.

Python version: 3.7


Description: SAA for single period newsvendor
"""
    
import numpy as np
import scipy.stats as st
from gurobipy import *


# 使用拉丁超立方抽样生成样本
def generate_sample(sample_num, trunQuantile, mu):
    samples = [0 for i in range(sample_num)]
    for i in range(sample_num):
        rand_p = np.random.uniform(trunQuantile*i/sample_num, trunQuantile*(i+1)/sample_num)
        samples[i] = st.poisson.ppf(rand_p, mu)
    return samples


price  =  6
vari_cost = 2
sal_value = 1
mean_demands = 5
sample_num = 10
trunQuantile = 0.999

samples = generate_sample(sample_num, trunQuantile, mean_demands)
   

try:
    # Create a new model
    m = Model("saa-mip")
    
    # Create variables
    Q = m.addVar(vtype = GRB.CONTINUOUS) # only one value for all scenarios           
    realized_demand = [m.addVar(vtype = GRB.CONTINUOUS) for s in range(sample_num)] # auxiliary variable
    remnant_demand = [m.addVar(vtype = GRB.CONTINUOUS) for s in range(sample_num)] # auxiliary variable
    profit = LinExpr() 
    
    profit = price * sum(realized_demand) / sample_num - vari_cost * Q + sal_value * sum(remnant_demand) / sample_num
    
    m.update()
    m.setObjective(profit, GRB.MAXIMIZE)
    
    # Add constraints
    for s in range(sample_num):
        m.addGenConstrMin(realized_demand[s], [Q], samples[s])
        m.addConstr(remnant_demand[s] == Q - realized_demand[s])
        
    m.update()
    m.optimize()
    
    print('ordering quantity in the first period is:  %.2f ' % Q.X)
    
except GurobiError as e:
        print('Error code ' + str(e.errno) + ": " + str(e))    
        
except AttributeError:
        print('Encountered an attribute error')