运筹系列87：julia求解随机动态规划问题入门

原创

已于 2024-01-10 11:11:45 修改 · 1.7k 阅读

24 ·

CC 4.0 BY-SA版权

文章标签：

#julia #动态规划 #代理模式

于 2023-12-18 11:17:35 首次发布

本文介绍了如何使用Julia的SDDP算法包解决随机动态规划问题，包括线性策略图示例、报童模型的随机需求优化、水力和火力发电调度，以及探索迷宫和牛奶制造的决策制定，展示了SDDP在处理随机环境下的决策问题中的应用。

随机动态规划问题的特点是：

有多个阶段，每个阶段的随机性互不相关，且有有限个实现值 (finite realizations)
具有马尔可夫性质，即每个阶段只受上一个阶段影响，可以用状态转移方程来描述阶段与阶段之间的变化过程。

我们使用julia的SDDP算法包来求解随机动态规划问题。

1. 入门案例：LinearPolicyGraph

看一个简单的数值优化的例子：
在这里插入图片描述
我们将其建立为一个N阶段的问题：

初始值为M。
使用SDDP.jl进行求解：

using SDDP
import Ipopt

M, N = 5, 3

model = SDDP.LinearPolicyGraph(
    stages = N,
    lower_bound = 0.0,
    optimizer = Ipopt.Optimizer,
) do subproblem, node
    @variable(subproblem, s >= 0, SDDP.State, initial_value = M)
    @variable(subproblem, x >= 0)
    @stageobjective(subproblem, x^2)
    @constraint(subproblem, x <= s.in)
    @constraint(subproblem, s.out == s.in - x)
    if node == N
        fix(s.out, 0.0; force = true)
    end
    return
end

SDDP.train(model)
println(SDDP.calculate_bound(model))
simulations = SDDP.simulate(model, 1, [:x])
for data in simulations[1]
    println("x_$(data[:node_index]) = $(data[:x])")
end

结果为

8.333333297473812
x_1 = 1.6666655984419778
x_2 = 1.6666670256548375
x_3 = 1.6666673693365108

非常接近理论最优值。

2 报童模型：加入随机变量

报童每天早上进一批报纸，需求有个大致的数，但是并不确定。那么报童要进多少报纸合适？
变量如下：

报纸成本c，零售价v
残值s，缺货损失p
需求为随机变量x，概率密度为f
决策变量为订货量Q
目标函数为max(利润-成本)，假如x>Q，则利润=vQ，成本=cQ+p(x-Q)
最佳订货量Q满足到Q的概率积分=(v+p-c)/(v+p-s)，推倒过程这里忽略

.我们缺货无损失，残值为0，零售价30.5，订货价20.5，则最佳订货量满足积分=0.328。若f为均值300，方差50的正态分布函数，则最佳订货量 = -0.445*50+300 = 278
我们使用Julia进行求解，首先要将正态分布离散化：

using Distributions,StatsPlots
D = Distributions.Normal(300,50)
N = 1000
d = rand(D, N)
P = fill(1 / N, N)
StatsPlots.histogram(d; bins = 80, label = "", xlabel = "Demand")

在这里插入图片描述
接下来是一个trick，这里是先决策再实现随机变量，因此把问题建成两阶段的模型，第一阶段决策为购买量u_make，第二阶段实现随机变量ω，随后决策为销售量u_sell。
模型如下：

using SDDP, HiGHS
model = SDDP.LinearPolicyGraph(;
    stages = 2,
    sense = :Max,
    upper_bound = 20 * maximum(d),  # The `M` in θ <= M
    optimizer = HiGHS.Optimizer,
) do subproblem::JuMP.Model, stage::Int
    @variable(subproblem, x >= 0, SD