动态规划——几个经典问题的实现（Python）
动态规划（dunamic programming，DP）是运筹学的一个分支，是将复杂的多段决策问题分解为若干相互关联的子决策问题以获得最优决策序列（决策链）的方法。由美国数学家贝尔曼(R.E.Bellman)于20世纪50年代提出，其基本原理是贝尔曼在《动态规划》(1957年)一书中所提出的最优化原理，即多段决策过程的最优决策序列具有这样的性质:不论初始状态和初始决策如何，对于前面决策所造成某一状态而言，余下的所有决策总构成一个最优策略。动态规划应用及其广泛，包括工程技术、经济、工业产业、军事以及自动化控制等领域，并在背包问题、生产经营问题、资金管理问题、资源分配问题、最短路径问题和复杂系统可靠性问题等中取得了显著的效果。
一、序贯决策—多阶段决策
在现实生活中，有一类活动的过程，由于他的特殊性，可将过程分成若干个互相联系的阶段，在他的每个阶段都需要作出决策，从而使整个过程达到最好的活动效果。因此各个阶段决策的选取不能任意确定，他依赖与当前面临的状态，又影响以后的发展。当各个阶段决策确定后，就组成一个决策序列，因而也就确定了整个过程的一条活动路线，这种把一个问题看作是一个前后关联具有链状结构的多阶段过程就称为多阶段决策过程，这种问题称为多阶段决策问题。在多阶段决策问题中，各个阶段采取的决策一般来说是与时间有关的，决策依赖于当前状态，又随即引起状态的转移，一个决策序列就是在变化的状态中产生出来的，故有“动态”的含义，称这种解决多阶段决策最优化的过程为动态规划方法。
多阶段决策问题中，各个阶段采用的决策，一般来说是与时间有关的，决策依赖于当前状态，又随即引起状态的转移，一个决策序列就是在变化的状态中产生出来的，故有动态的含义，称这种解决多阶段决策最优化问题的方法为动态规划方法。
在解决问题的过程中，需要经历多个阶段。每个决策阶段产生一组状态，最终通过某组决策序列，得出最优解。使用动态规划的问题必须满足最优化原理和无后效性。
1.1 无后效性
无后效性也叫马尔科夫性，是将各阶段按照一定的次序列排列好之后，对于某个给定的阶段状态，它以前各阶段的状态无法直接影响它未来的决策而只能通过当前的这个状态。换句话说，每个状态都是过去历史的一个完整总结，这就是无后向性，无后效性。
1.2 自相似性
如果一个问题最优解包括规模更小相同问题的最优解，就可以用动态规划来解决。动态规划的自相似性是阶段划分的依据，每一个阶段或者后续子过程都与整体具有模式上的相似性，本质是动态规划问题可分性、可拆分子过程求解的基础。
1.3 子问题的重叠性
动态规划算法的关键在于解决冗余，这是动态规划算法的根本目的。动态规划实质上是一种以空间换时间的技术，他在实现的过程中，不得不存储产生过程中的各种状态，所以他的空间复杂度要大于其他的算法。选择动态规划算法是因为动态规划算法在空间上可以承受，而搜索算法在时间上无法承受，所以我们舍空间而取时间。
动态规划的数学模型，根据决策过程的演变是确定性的还是随机性的，可分为确定性决策过程和随机性决策过程；按照时间参量是离散的或是连续的变量，可分为离散决策过程和连续决策过程。组合起来就有离散确定性、离散随机性、连续确定性、连续随机性四种决策过程模型。对于确定性的决策过程，问题中的下一段的状态已由当前段的状态及决策完全确定。对于随机性决策过程，它与确定性决策过程的区别在于下一段的状态并不能由当前段的状态及决策完全确定，而是按某种概率分布来决定下一段的状态。这种概率分布是由当前段的状态和策略完全确定。
二、动态规划基本思想
动态规划算法通常用于求解具有某种最优性质的问题。在这类问题中可能会有许多可行解。每一个解都对应于一个值，我们希望找到具有最优值的解。动态规划算法与分治法类似，其基本思想也是将待求解问题分解为若干个子问题，先求解子问题，然后从这些子问题的解得到原问题的解。与分治法不同的是，适合于用动态规划求解的问题，经分解得到子问题往往不是互相独立的。若用分治法来解这类问题，则分解得到的子问题数目太多，有些子问题被重复计算了很多次。如果我们能够保存已解决的子问题的答案，而在需要时再找出已求得的答案，这样就可以避免大量的重复计算，节省时间。我们可以用一个表来记录所有已解决的子问题的答案。不管该子问题以后是否被用到，只要他被计算过，就将其结果填入表中，这就是动态规划的基本思路。
最优化原理可以这样阐述：一个最优化策略具有这样的性质，不论过去状态和决策如何，对前面的决策所形成的状态而言，余下的诸决策必须构成最优策略，简言之，一个最优化策略的子策略总是最优的。一个问题满足最优化原理又称其具有最优子结构性质。
三、动态规划的术语
多阶段决策问题：如果一类活动过程可以分为若干个互相联系的阶段，在每一个阶段都需做出决策（采取措施），一个阶段的决策确定以后，常常影响到下一个阶段的决策，从而就完全确定了一个过程的活动路线，则称为多阶段决策问题。各个阶段的决策构成一个决策序列，称为一个策略。每一个阶段都有若干个决策可供选择，因而有许多策略供我们选取，对应于一个策略可以确定活动的效果，这个效果可以用数量来确定，策略不同，效果也不同，多阶段策略问题就是要在可以选择那些策略中间选取一个最优策略，使在预定的标准下达到最好的效果。
阶段：把所给求解问题的过程恰当地分成若干个相互联系的阶段，以便求解，过程不同，阶段数就可能不同，描述阶段的变量称为阶段变量。在多数情况下，阶段变量是离散的，用k表示，此外，也有阶段变量是连续的情形。如果过程可以在任何时刻作出决策，且在任意两个不同时刻之间允许有无穷多个决策时，阶段变量就是连续的。
状态：状态表示每个阶段开始面临的自然状态或客观条件，它不以人们的主观意志为转移，也称为不可控因素。在上面的例子中状态就是某阶段的出发位置，它既是该阶段某路的起点，同时又是前一阶段某支路的终点。
决策：一个阶段的状态给定以后，从该状态演变到下一阶段某个状态的一种选择（行为）称为决策。在最优控制中，也称为控制。在许多问题中，决策可以自然而然地表示为一个数或一组数。不同的决策对应着不同的数值，描述决策的变量称决策变量，因状态满足无后效性，故在每一个阶段选择决策时只需考虑当前的状态而无需考虑过程的历史。决策变量的范围称为允许决策集合。
策略：由每个阶段的决策组成的序列称为策略。对于每一个实际的多阶段策略过程，可供选取的策略有一顶的范围限制，这个范围称为允许策略集合，允许策略集合中达到最优效果的策略称为最优策略。
状态转移方程：给定k阶段状态变量x(k)的值后，如果这一阶段的决策变量一经确定，第k+1阶段的状态变量x(k+1)也就完全确定，即x(k+1)的值随x(k)和第k阶段的决策u(k)的值变化而变化，那么可以把这一关系看成x(k)，u(k)与x(k+1)确定的对应关系，用x(k+1)=Tk(x(k),u(k))表示。这是从k阶段到k+1阶段的状态转移规律，称为状态转移方程。
最优化策略函数：作为整个过程的最优策略函数，它满足：相对前面决策所形成的状态而言，余下的子决策必然构成“最优子策略”。最优性原理实际上要求问题的最优策略的子策略也是最优。
四、经典的动态规划问题
4.1 背包问题
一位旅行者携带背包去登山,已知他所能承受的背包重量限度为\(a(kg)\)，现有\(n\)种物品可供他选择装入背包,第\(i\)种物品的单件重量为\(a_i(kg)\)，其价值(可以是表明本物品对登山的重要性的数量指标)是携带数量\(x_i\)的函数\(v(x_i）(i=1，2，...，n)\)，问旅行者应如何选择携带各种物品的件数，以使总价值最大？
#背包问题4件商品
weights = [2, 3, 4, 5] # 每个物品的重量
values = [3, 4, 5, 6] # 每个物品的价值
max_weight = 8 # 背包的最大承重
def knapsack(weights, values, max_weight):
n = len(weights) # 获取物品的数量
dp = [0] * (max_weight + 1) # 初始化动态规划数组，dp[j]表示承重为j时的最大总价值
for i in range(n): # 遍历每个物品
for j in range(max_weight, weights[i]-1, -1): # 逆序遍历每个承重
if j >= weights[i]: # 如果当前物品的重量不超过当前承重j
dp[j] = max(dp[j], dp[j-weights[i]] + values[i]) # 选择放入物品i或不放入物品i
else: # 如果当前物品的重量超过当前承重j
dp[j] = dp[j] # 不放入物品i
return dp[max_weight] # 返回最大总价值
max_value = knapsack(weights, values, max_weight) # 求解最大总价值
print(max_value)
背包装的最大价值为10。
4.2 资源配置问题
某集团公司拟进行设备升级换代，购进某种高效率机器5台，分配给下属甲、乙、丙三厂，各工厂使用这些设备以后可以产生的效益如下表。问这5台机器如何分配，集团所获总收益为最大？
| 机器数 | 甲 | 乙 | 丙 |
|---|---|---|---|
| 0 | 0 | 0 | 0 |
| 1 | 3 | 5 | 4 |
| 2 | 7 | 10 | 6 |
| 3 | 9 | 11 | 11 |
| 4 | 12 | 11 | 12 |
| 5 | 13 | 11 | 12 |
import numpy as np
#资源离散分配问题
mat = np.array([[0,0,0,],
[3,5,4],
[7,10,6],
[9,11,11],
[12,11,12],
[13,11,12]])
mat_new = mat.copy()
m = mat.shape[0]
n = mat.shape[1]
bestChoice = np.zeros((m,n)).astype(int)
for k in range(n-1):
for i in range(m):
max = 0
for j in range(i+1):
if mat[j][n-2-k] + mat_new[i-j][n-1-k] > max:
max = mat[j][n-2-k] + mat_new[i-j][n-1-k]
bestChoice[i][n-2-k] = j
bestChoice[i][n-1-k] = i-j
mat_new[i][n-2-k] = max #根据算出来的每一阶段的最优价格不断更新价格表与最佳选择表
print("最大盈利为：")
print(mat_new[m-1][0])
print("最佳分配方案为：")
print(bestChoice[m-1][0])
for i in range(n-2):
print(bestChoice[m-1][i+1] - bestChoice[m-1][i+2])
print(bestChoice[-1][-1])
分配方案：甲0台，乙2台，丙3台。
4.3 最短路径问题
已知6个城市之间的机票票价(单位：百元)如矩阵所示，求城市 A 到其它城市的票价最便宜的路径及票价，其中若值为0代表两地之间无直航，若值不为零，表示票价。
| A | B | C | D | E | F |
|---|---|---|---|---|---|
| A | 0 | 50 | 0 | 40 | 25 |
| B | 50 | 0 | 15 | 20 | 0 |
| C | 0 | 15 | 0 | 10 | 20 |
| D | 40 | 20 | 10 | 0 | 10 |
| E | 25 | 0 | 20 | 10 | 0 |
import networkx as nx
# 定义机票票价矩阵
costs = [[0, 50, 0, 40, 25, 10],
[50, 0, 15, 20, 0, 25],
[0, 15, 0, 10, 20, 0],
[40, 20, 10, 0, 10, 25],
[25, 0, 20, 10, 0, 55],
[10, 25, 0, 25, 55, 0]]
# 将矩阵转换成带权有向图
G = nx.DiGraph()
for i in range(len(costs)):
for j in range(len(costs[0])):
if costs[i][j] != 0:
G.add_edge(chr(65+i), chr(65+j), weight=costs[i][j])
# 计算城市A到其它城市的最短路径及票价
shortest_path = nx.shortest_path(G, source='A', weight='weight')
shortest_path_cost = nx.shortest_path_length(G, source='A', weight='weight')
for dest, cost in shortest_path_cost.items():
print(f"从A到{dest}的最便宜路径为{shortest_path[dest]}, 票价为{cost}")
从A到A的最便宜路径为['A'], 票价为0
从A到F的最便宜路径为['A', 'F'], 票价为10
从A到E的最便宜路径为['A', 'E'], 票价为25
从A到B的最便宜路径为['A', 'F', 'B'], 票价为35
从A到D的最便宜路径为['A', 'F', 'D'], 票价为35
从A到C的最便宜路径为['A', 'E', 'C'], 票价为45
4.4 机器负荷问题
某机器可以在高、低两种不同的负荷下进行生产。高负荷下生产时，年产量\(g=8u\)，其中\(u\)为投入生产的机器数量，机器的年折损率为0.7，即年初完好的机器数量为\(u\)，年终就只剩下0.7\(u\)台是完好的，其余均需维修或报废。在低负荷下生产，年产量\(g=5u\)，其中\(u\)为投入生产的机器数量，机器的年折损率为0.9，要求制定一个五年计划，在每年开始时决定如何重新分配好机器在两种不同负荷下工作的数量，使得五年的总产量最高。
from mip import Model, xsum, maximize, INTEGER
def optModel(m, n, g, h, a, b):
model = Model("machine_load_distribution")
x1 = [model.add_var(var_type=INTEGER) for i in range(n)] #高负荷运行的机器数量
x2 = [model.add_var(var_type=INTEGER) for i in range(n)] #低负荷运行的机器数量
model.objective = maximize(xsum(x1[i] * g + x2[i] * h for i in range(n)))
model += x1[0] + x2[0] == m
for i in range(1,n):
model += x1[i] + x2[i] == x1[i-1] * a + x2[i-1] * b
model.verbose = 0
model.optimize()
print(model.objective_value)
lis1 = []
lis2 = []
for i in range(n):
lis1.append(x1[i].x)
lis2.append(x2[i].x)
print(lis1)
print(lis2)
if __name__ == '__main__':
m, n, g, h, a, b = 1000, 5, 8, 5, 0.7, 0.9
optModel(m, n, g, h, a, b)
五年的总收益为23687.0，其中各年初高负荷的机器为[0.0, 0.0, 795.0, 570.0, 399.0]；低负荷的机器为[1000.0, 900.0, 15.0, 0.0, 0.0]。注意：答案和好多书上有些出入，这是因为考虑了分配高低负荷的机器数为整数，不考虑该因素的总收益应为23691.2，好多教材答案有误。还可以考量终端状态固定的情形，如要求在第5年末完好的机器数量是500台，计算其总收益过程略。
4.5 设备更新问题
某种工程设备的役龄为 4 年, 每年年初都面临着是否更新的问题: 若卖旧买新, 就要支付一定的购置费用; 若继续使用, 则要支付更多的维护费用, 且使用年限越长维护费用越多。役龄期内每年的年初购置价格, 当年维护费用及年末剩余净值如下表所示。为该设备制定一个 4 年役龄期内的更新计划, 使总的支付费用最少。
| 年份 | 1 | 2 | 3 | 4 |
|---|---|---|---|---|
| 年初购置价格 (万元） | 25 | 26 | 28 | 31 |
| 当年维护费用 (万元） | 10 | 14 | 18 | 26 |
| 年末剩余净值 (万元） | 20 | 16 | 13 | 11 |
可以把这个问题化为图论中的最短路问题。
构造赋权有向图\(D=(V, A, W)\)，其中顶点集 \(V=\{v_{1}, v_{2}, \cdots, v_{5}\}\)，这里 $v_{i} (i=1,2,3,4) $表示第 \(i\) 年年初的时刻，\(v_{5}\)表示第 4 年年末的时刻；\(A\)为边集；邻接矩阵\(W=(w_{i j})(5 \times 5)\)，这里$ w_{ij}$ 为第 \(i\)年年初至第\(j\)年年初 (或\(j−1\)年末) 期间所支付的费用，计算公式为
其中\(p_{i}\)为第\(i\)年年初的购置价格，\(a_{k}\)为使用到第 \(k\)年当年的维护费用，\(r_{i}\)为使用\(i\)年旧设备的出售价格，邻接矩阵
则制定总的支付费用最小的设备更新计划，就是求有向图\(D\)从顶点\(v_{1}\)到顶点\(v_{5}\)的最短路。
import numpy as np
import networkx as nx
import pylab as plt
p=[25,26,28,31]; a=[10,14,18,26]; r=[20,16,13,11]
b=np.zeros((5,5)); #初始化邻接矩阵
for i in range(5):
for j in range(i+1,5):
b[i,j]=p[i]+np.sum(a[0:j-i])-r[j-i-1];
print(b)
G=nx.DiGraph(b)
p=nx.dijkstra_path(G, source=0, target=4, weight='weight') #求最短路径；
print("最短路径为:",np.array(p)+1) #python下标从0开始
d=nx.dijkstra_path_length(G, 0, 4, weight='weight') #求最短距离
print("所求的费用最小值为：",d)
s=dict(zip(range(5),range(1,6))) #构造用于顶点标注的标号字典
plt.rc('font',size=16)
pos=nx.shell_layout(G) #设置布局
w=nx.get_edge_attributes(G,'weight')
nx.draw(G,pos,font_weight='bold',labels=s,node_color='r')
nx.draw_networkx_edge_labels(G,pos,edge_labels=w)
path_edges=list(zip(p,p[1:]))
nx.draw_networkx_edges(G,pos,edgelist=path_edges,edge_color='r',width=3) #画出最短路径
#plt.savefig("figure10_9.png",pdi=500);
plt.show()
[[ 0. 15. 33. 54. 82.]
[ 0. 0. 16. 34. 55.]
[ 0. 0. 0. 18. 36.]
[ 0. 0. 0. 0. 21.]
[ 0. 0. 0. 0. 0.]]
最短路径为: [1 2 3 5]
所求的费用最小值为： 67.0
4.6 生产库存问题
工厂生产某种产品，每单位（千件）的成本为 1（千元），每次开工的固定成本为 3 （千元），工厂每季度的最大生产能力为 6（千件）。经调查，市场对该产品的需求量第一、二、三、四季度分别为 2，3，2，4（千件）。如果工厂在第一、二季度将全年的需求都生产出来，自然可以降低成本（少付固定成本费），但是对于第三、四季度才能上市的产品需付存储费，每季每千件的存储费为 0.5（千元）。还规定年初和年末这种产品均无库存。试制定一个生产计划，即安排每个季度的产量，使一年的总费用（生产成本和存储费）最少。
'''
这一类生产计划问题（Production planning problem），阶段按计划时间自然划分，状态定义为每阶段开始时的储存量x(k)，决策为每个阶段的产量u(k)，记每个阶段的需求量（已知量）为d(k)，则状态转移方程为:
x(k+1) = x(k) + u(k) − d(k) ( x(k)≥0, k=1,2,...,n )
设每阶段开工的固定成本费为a，生产单位数量产品的成本费为b，每阶段单位数量产品的储存费为c，阶段指标为阶段的生产成本和储存费之和，即
if u(k)>0
v(x(k), u(k)) = cx(k) + (a + bu(k))
if u(k)=0
v(x(k), u(k)) = cx(k) + 0
最优值函数f(x(k))为从第k段的状态x(k)出发到过程终结的最小费用，满足 :
f(x(k)) = min(v(x(k), u(k)) + f(x(k+1)))
'''
# 每阶段单位数量产品的储存费为 C
C = 0.5
# 每阶段开工的固定成本费为 A
A = 3
# 生产单位数量产品的成本费为 B
B = 1
# 市场对该产品的需求量数组
demand = [2, 3, 2, 4]
# 计算市场对该产品的需求量
def demand_function(k):
return demand[k - 1]
# 生产目标，所有初始化为 None，在每次遍历的时候，需要更新这个值
production_target = [None] * len(demand)
def production(k):
return production_target[k - 1]
# 计算每阶段开始时的储存量
def storage(k):
if k <= 1:
return 0
else:
return storage(k - 1) + production(k - 1) - demand_function(k - 1)
# 计算阶段的生产成本和储存费之和
def cost(k):
if production(k) == 0:
return C * storage(k)
else:
return C * storage(k) + (A + B * production(k))
# 计算总成本
def total_cost(k):
if k == len(demand) + 1:
return 0
else:
return cost(k) + total_cost(k + 1)
# 生成所有可能的生产计划组合
possible_plans = []
for s1 in range(sum(demand) + 1):
for s2 in range(sum(demand) + 1):
for s3 in range(sum(demand) + 1):
for s4 in range(sum(demand) + 1):
if sum([s1, s2, s3, s4]) == sum(demand):
possible_plans.append([s1, s2, s3, s4])
# 筛选出符合条件的生产计划组合
for i in range(len(demand)):
possible_plans = [plan for plan in possible_plans if sum(plan[:i + 1]) >= sum(demand[:i + 1])]
# 计算每种生产计划组合的总成本
total_costs = []
for plan in possible_plans:
production_target = plan
total_costs.append(total_cost(1))
# 找到最小总成本和对应的最优生产计划
min_cost = min(total_costs)
min_index = total_costs.index(min_cost)
optimal_plan = possible_plans[min_index]
# 输出结果
print("最小总成本:", min_cost)
print("最优生产计划:", optimal_plan)
4.7 钢条切割问题
一家钢条回收厂推出了不同长度的钢条对应不同价格的回收策略，参看下表。如果小鱼现在有\(n\)米长的钢条，应该怎么切割使得销售后收益\(R_n\)最大？
| 长度\(n\) | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| 价格\(p\) | 1 | 5 | 8 | 9 | 10 | 17 | 17 | 20 | 24 | 30 |
对于\(n\)长度的钢条，每个空隙都可以选择切和不切，一共\(n-1\)个空隙，所以最后有$ 2^n-1$ 种切割方法，然后再考虑每种切割方法的价格，最后进行比较。别的不说，光是第一步时间复杂度就已经达到\(O(2^n)\)了。
#钢条切割问题：自顶向下（由大到小）
#
#获得最大值函数
def max(a,b):
temp = a;
if temp < b:
temp = b;
return temp
# 备忘机制的CUT-ROD
def MEMOIZED_CUT_ROD_AUX(p, n, r): # 实际上参数是输入的数字+1
n = n - 1
# 此时的n是输入的数
if r[n] >= 0:
return r[n]
if n<=10:
if n == 0:
q = 0
else:
q = -9999
for i in range(1, n):
q = max(q, int(p[i]) + int(MEMOIZED_CUT_ROD_AUX(p, n - i, r)))
else:
# 对于n>10的情况，同样要进行拆分成两个子问题再比较最大值，比如n=12应该还是得考虑max（1+11，2+10，3+9...）
q = -9999
for i in range(0, n):
q = max(q, int(r[i+1]) + int(MEMOIZED_CUT_ROD_AUX(p, n - i, r)))
r[n] = q
print("n=",n)
print("r=",r)
return q
def MEMOIZED_CUT_ROD(p,n):
r = {}
for i in range(0,n):
r[i] = -9999
r[1] = p[0]
return MEMOIZED_CUT_ROD_AUX(p,n,r)
if __name__ == '__main__':
p = [1,5,8,9,10,17,17,20,24,30]
#长度 i 0 1 2 3 4 5 6 7 8 9 10
#价格 pi 0 1 5 8 9 10 17 17 20 24 30
print("输入钢条长度n：")
n = int(input())
if n == 1:
print("最大的收益：",p[0])
print("最大的收益：",MEMOIZED_CUT_ROD(p,n+1))
总结
动态规划算法（Dynamic Programming）是一种高效的算法思想，其核心思想是将复杂的问题拆分成多个子问题，然后通过已经求解过的子问题的解来推导出当前问题的解。动态规划对于解决多阶段决策问题的效果是明显的，但是动态规划也有一定的局限性，它没有统一的处理方式，必须根据问题的各种性质并结合一定的技巧来处理；另外当变量的维度增大时，总的计算量及存储量急剧增大，这就是“维数障碍”。动态规划在经济管理、生产调度、工程技术和最优控制、语音识别、图像拼接等方面都得到了广泛应用，例如最短路线、库存管理、资源分配、设备更新、排序、装载等问题，用动态规划方法比用其他方法求解更方便。
动态规划求解步骤如下：
（1）划分：按照问题的特征，把问题分为若干阶段。注意：划分后的阶段一定是有序的或者可排序的，初始状态→│决策1│→│决策2│→…→│决策n│→结束状态
（2）确定状态和状态变量：将问题发展到各个阶段时所处的各种不同的客观情况表现出来。状态的选择要满足无后续性
（3）确定决策并写出状态转移方程：状态转移就是根据上一阶段的决策和状态来导出本阶段的状态。根据相邻两个阶段状态之间的联系来确定决策方法和状态转移方程
（4）边界条件：状态转移方程是一个递推式，因此需要找到递推终止的条件
参考文献
