week6 常见数据规划-Santa服务调度

规划问题

常见规划问题

线性规划是运筹学的一个重要分支，广泛应用于军事作战、经济分析、经营管理和工程技术等方面。为合理地利用有限的人力、物力、财力等资源作出的最优决策，提供科学的依据。

• LP：Linear Programming 线性规划

•   研究线性约束条件下线性目标函数的极值问题
  

• ILP：Integer Linear Programming 整数线性规划

•   全部决策变量必须为整数
  

• MIP：Mixed Integer Programming 混合整数规划

•   混合整数规划是LP的一种，其中部分的决策变量是整数（不要求全部都是整数）
  

• VRP：Vehicle Routing Problem 车辆路径问题

规划求解步骤

• Step1，列出约束条件及目标函数
• Step2，画出约束条件所表示的可行域
• Step3，在可行域内求目标函数的最优解及最优值

规划工具

• pulp

•   只用于线性模型，包括如整数规划、01规划，还是混合整数线性规划 MILP
  

• ortools

•   开源软件。可以解决车辆路径、流程、整数和线性规划等问题。提供了C++,Python,Java,.NET接口
  

pulp

• LpProblem类，用来构造LP问题实例

# Name，指定问题名，输出信息用
# Sense，LpMinimize或LpMaximize，代表目标是极大值还是极小值
lpp = LpProblem(name='NoName', sense=LpMinimize)
# 在对LpProblem添加完约束条件后，调用solve进行求解
lpp.solve()
# 在对LpProblem添加完约束条件后，调用solve进行求解
lpSum(vector) 

• LpVariable类 ，用来构造LP问题中的变量

# name指定变量名
# lowBound(默认负无穷)和upBound(默认正无穷)是下界和上界，
# cat用来指定变量是离散(Integer,Binary)还是连续(Continuous) 

LpVariable(name, lowBound=None, upBound=None, cat='Continuous', e=None) 

Ortools

线性规划，默认使用GLOP
整数规划，默认使用CBC（Coin-or branch and cut），还包括SCIP、GLPK、Gurobi等

执行流程

# 1. Solver创建
solver = pywraplp.Solver.CreateSolver('SCIP')
# 2. 变量设置
solver.NumVar：创建普通变量
solver.IntVar：创建整数变量

infinity = solver.infinity() # 正无穷
x = solver.IntVar(0.0, infinity, 'x')
print('变量数量：', solver.NumVariables())


# 3. 添加约束条件
solver.Add(x + 7 * y <= 17.5)
print('约束的数量：', solver.NumConstraints())
# 4. Solve求解
# 求解最大值问题
solver.Maximize(x + 10 * y)
status = solver.Solve()
# 5. Solve的结果
print('目标值 =', solver.Objective().Value())
print('x =', x.solution_value())
print('y =', y.solution_value())

实例操作：Santa服务调度

• Step1, 数据加载

data = pd.read_csv('./data/family_data.csv',index_col='family_id')

• Step2，数据预处理
  • 1）计算Perference Cost矩阵 pcost_mat
  • 2）计算Accounting Cost矩阵 acost_mat
  • 3）计算每个家庭的人数 FAMILY_SIZE
  • 4）每个家庭的倾向选择（choice_） DESIRED

N_DAYS = 100 # 安排的天数
N_FAMILY = 5000 #家庭ID的个数
MIN_OCCUPANCY = 125 # 最小承载量
MAX_OCCUPANCY = 300 # 最大承载量
# 计算pcost_mat，每个家庭在什么时候（day0-99）访问时的penalty
# 大小5000*100的矩阵

# 1. 计算Perference Cost矩阵 pcost_mat
pcost_mat = np.full(shape=(N_FAMILY,100),fill_value=999999)
for f in range(N_FAMILY):
    # 家庭成员数
    f_num = data.loc[f,'n_people']
    # 对于第f个家庭，初始化pcost = other choice下的penalty
    pcost_mat[f,:] = get_penalty(f_num,10) #初始值最大惩罚
    # 计算choice 0-9 的penalty
    for choice in range(10):# choice 0-9
        temp = data.loc[f][choice] #choice的天数
        penalty = get_penalty(f_num,choice) # 得到对应choice的惩罚
        pcost_mat[f,temp-1] = penalty # 因为下标是从0开始，所以要在天数基础上-1才是下标值

# 2.计算Accounting Cost矩阵 acost_mat
acost_mat = np.zeros(shape=(MAX_OCCUPANCY+1,MAX_OCCUPANCY+1),dtype=np.float64)
for i in range(acost_mat.shape[0]):# 当天安排的人数
    for j in range(acost_mat.shape[1]):# 前一天安排的人数
        diff = abs(i-j)
        acost_mat[i,j] = max(0,(i - 125) / 400 * i ** (0.5 + diff/50.0))

# 3.计算每个家庭的人数 FAMILY_SIZE
FAMILY_SIZE = data['n_people'].values

# 4.每个家庭的倾向选择（choice_） DESIRED
DESIRED = data.values[:,:-1]-1

• Step3，使用LP和MIP求解 规划方案
  • 1）先使用LP 对绝大部分家庭进行规划
  • 2）再使用MIP 对剩余家庭进行规划
  • 3）汇总两边的结果 => 最终规划方案

#进行 使用整数规划求解
def solveIP(families,min_occupancy,max_occupancy):
    # 创建求解器
    solver = pywraplp.Solver('AssignmentProblem', pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)
    # 需要安排的家庭
    n_families = len(families)
    x = {} # family_id在第j天是否参观
    candidates = [[] for x in range(N_DAYS)] # 定义了len为100的list
    for i in families:#family_id
        for j in DESIRED[i,:]:# family_id的choice
            #print(j)
            candidates[j].append(i) # 在第j天，有第i个family参观
            x[i,j] = solver.BoolVar('x[%i,%i]' %(i,j)) # x[%i,%i]中的i代表integer类型
  
    daily_occupancy = [solver.Sum([x[i,j] * FAMILY_SIZE[i] for i in candidates[j]])\
                        for j in range(N_DAYS)] # j代表1-100天
    
    family_presence = [solver.Sum(x[i,j]  for j in DESIRED[i,:])\
                        for i in families]

    preference_cost = solver.Sum([pcost_mat[i,j] * x[i,j] for i in families\
                                for j in DESIRED[i,:]])
    
    # 满足preference_cost最小
    solver.Minimize(preference_cost)    
    
    # 每个家庭都在10个choice中出现一次
    for i in range(n_families):
        solver.Add(family_presence[i]==1)
    # 每天访问人数约束
    for j in range(N_DAYS):
        solver.Add(daily_occupancy[j]>=min_occupancy[j])
        solver.Add(daily_occupancy[j]<=max_occupancy[j])
    result = solver.Solve()
    
    temp = [(i,j) for i in families\
             for j in DESIRED[i,:] if x[i,j].solution_value()>0]
    # 计算剩余家庭的安排
    df = pd.DataFrame(temp,columns=['family_id','day'])
    return df

• Step4, 结果评估
  按照evaluation标准，计算
  Score = preference cost + accounting penalty

def cost_function(prediction):
    penalty,daily_occupancy = pcost(prediction) #统计preference cost和每天承载数量
    accounting_cost,num_out_of_range = acost(daily_occupancy) # 根据每天承载数量计算accounting cost
    final_score = penalty + accounting_cost + num_out_of_range * 9999999
    return final_score

• Step5，方案优化
  通过更换family_id的选择，查找更好的score
  每次更换后，都对方案进行评估，选择更小的score方案

# 寻找更好的替代方案
def find_better(pred):
    fobs = np.argsort(FAMILY_SIZE) # 返回数组从小到大的索引
    #print(np.sort(FAMILY_SIZE)) # 对FAMILY_SIZE按从小到大的顺序排序
    score = cost_function(pred)
    original_score = np.inf #打擂法 正无穷初始值
    
    # 如果找不到更新则退出
    while score < original_score:
        original_score = score
        for family_id in fobs:
            for pick in range(10):
                # 得到family_id在choice pick的日期day
                day = DESIRED[family_id, pick]
                # 该family的原有日期oldvalue
                oldvalue = pred[family_id]
                # 将原有oldvalue替换为现在的choice pick
                pred[family_id] = day
                # 重新计算分数
                new_score = cost_function(pred)
                # 如果比原来分数小，更新choice成功
                if new_score < score:
                    score = new_score
                else:# 设置为原来的oldvalue
                    pred[family_id] = oldvalue
        print(score,end='\r')
    print(score)