利用MATLAB求解离散动态规划问题

14 / 28
利用MATLAB求解离散动态规划问题
三、动态规划求解函数-算例1求解过程分析
1. 阶段变量k： 将问题按工程分为4个阶段，即k = 1, 2, 3, 4 2. 状态变量sk ： 状 态 变 量sk 表 示 第k个 工 程 至 第4个 工 程 可 用 的 投 资 总 额 ， 显 然 有s1 = 60
穷举所有可能的结果，找出其中最优值及对应的策略
3 / 28
利用MATLAB求解离散动态规划问题
二、穷举法：算例1-算法与程序
对 于上 述 问 题 ， 在 穷 举 求 解 时 应 注 意 投 资 额 不 多 于60万 元 的 约 束。稍加分析便知道60万元必须全部用完没有剩余，否则结果一 定不会是最优解。因此可以写出以下程序: 输入利润函数： V 1 = [0, 20, 50, 65, 80, 85, 85]; V 2 = [0, 20, 40, 50, 55, 60, 65]; V 3 = [0, 25, 60, 85, 100, 110, 115]; V 4 = [0, 25, 40, 50, 60, 65, 70];
10 / 28
利用MATLAB求解离散动态规划问题
二、穷举法：算例2-算法与程序（续）
n = size(R, 1); Path = zeros(50, 5); Path(1 : 50, 1) = 1; Path(1 : 50, 5) = n; D = inf ∗ ones(50, 1); m = 1; for i = 1 : n; if R(1, i) ∼= inf for j = 1 : n if R(i, j) ∼= inf for k = 1 : n if R(j, k) ∼= inf Path(m, 2 : 4) = [i, j, k]; D(m) = 0; D(m) = D(m)+ R(1, i)+ R(i, j)+ R(j, k)+ R(k, n); m = m + 1; end ...
利用MATLAB求解离散动态规划问题
三、动态规划求解函数-算例1求解过程分析（续）
5. 阶段值Vk ： 第k阶段的阶段值表示投资给第k个工程uk 时的利润 因此可写出以下阶段值函数： function u) V = ObjFun(k, s, 0 0 0 0 20 20 25 25 50 40 60 40 w= 65 50 85 50 80 55 100 60 85 60 110 65 85 65 115 70 w = −w; V = w(u/10 + 1, k);
输入参数： x是状态变量，一列代表一个阶段的状态 函数DecisFun(k, x, u)描述可行决策集 ObjFun(k, s, u)描述阶段值 TransFun(k, s, u)描述状态转移方程 输出变量： p opt =[序号组，状态变量组，最优决策组，阶段值组] fval也是一个列向量，表示最优值
15 / 28
利用MATLAB求解离散动态规划问题
三、动态规划求解函数-算例1求解过程分析（续）
3. 决策变量uk ： 允许决策集：0 ≤ uk ≤ sk ，特别地，当k = 4，uk = sk 因此可写出以下可行决策集函数： function u = DecisFun(k, s, u) if k == 4 u = s; else u = 0 : 10 : s; end
19 / 28
利用MATLAB求解离散动态规划问题
三、动态规划求解函数-算例1求解过程分析（续）
调用动态规划求解函数： [ p opt, fval ] = dynprog(s, DecisFun , ObjFun , TransFun ) 1 2 p opt = 3 4 fval = −160 对于fval，因为之前把阶段值转变为负值，这个结果是-160，而对 于我们的问题中应是160。 运用动态规划求解函数求解总共耗时0.017秒，相比穷举法耗时更 多，因为在阶段数和状态数较少时，求解函数不能体现出优势， 调用函数的过程涉及各种参数的传递和运算，这是相对比较耗时 的。 60 40 40 10 20 0 30 10 −50 0 −85 −25
4 / 28
利用MATLAB求解离散动态规划问题
二、穷举法：算例1-算法与程序（续）
m = 1; J = zeros(84, 1); U = zeros(84, 4); for i = 1 : 7 for j = 1 : 7 for k = 1 : 7 for l = 1 : 7 if i + k + j + l == 10 J (m) = V 1(i) + V 2(j) + V 3(k) + V 4(l); U (m, 1) = i; U (m, 2) = j; U (m, 3) = k; U (m, 4) = l; m = m + 1; else continue; end ...
[inf , 5, 3, inf , inf , inf , inf , inf , inf , inf , inf inf , inf , inf , 1, 3, 6, inf , inf , inf , inf , inf inf , inf , inf , inf , 8, 7, 6, inf , inf , inf , inf inf , inf , inf , inf , inf , inf , inf , 6, 8, inf , inf inf , inf , inf , inf , inf , inf , inf , 3, 5, inf , inf R = inf , inf , inf , inf , inf , inf , inf , inf , 3, 3, inf inf , inf , inf , inf , inf , inf , inf , inf , 8, 4, inf inf , inf , inf , inf , inf , inf , inf , inf , inf , inf , 3 inf , inf , inf , inf , inf , inf , inf , inf , inf , inf , 2 inf , inf , inf , inf , inf , inf , inf , inf , inf , inf , 2 inf , inf , inf , inf , inf , inf , inf , inf , inf , 态规划问题
二、穷举法：算例2
例 6.1 最短路问题。如下图所示，其中小圆圈称为点，两点间的连线称 为弧，弧上的数字称为弧长，箭头表示前进方向。要寻找一条从 起点A到终点E的连通弧，使其总弧长最短。
7 / 28
利用MATLAB求解离散动态规划问题
二、穷举法：算例2
试图利用MATLAB并用最简单易懂的方法解决讲义中所有例题与 练习
2 / 28
利用MATLAB求解离散动态规划问题
二、穷举法：算例1
例 2.1 动 态 规 划 讲 义 第35页 习 题2， 设 国 家 拨 给60万 元 投 资 ， 供 四 个 工 程扩建使用，每个工程扩建后的利润与投资额的大小有关，投资 后的利润函数如下表所示。问如何投资使得总利润最大化？ 利润\投资 v1 (x) v2 (x) v3 (x) v4 (x) 0 0 0 0 0 10 20 20 25 25 20 50 40 60 40 30 65 50 85 50 40 80 55 100 60 50 85 60 110 65 60 85 65 115 70
12 / 28
利用MATLAB求解离散动态规划问题
二、穷举法：优缺点
优点： 程序简单，容易理解；在阶段数量与状态数量并不很大时运算速 度快
缺点： 并未运用动态规划思想，对不同问题通用性低，程序需要较大改 动
13 / 28
利用MATLAB求解离散动态规划问题
三、动态规划求解函数
function [ p opt, fval ] = dynprog(x, DecisFun, ObjFun, TransFun)
18 / 28
利用MATLAB求解离散动态规划问题
三、动态规划求解函数-算例1求解过程分析（续）
6. 各阶段状态变量可能取值： 由已知： s1 = {60} s2 = {0, 10, 20, 30, 40, 50, 60} s3 = {0, 10, 20, 30, 40, 50, 60} s4 = {0, 10, 20, 30, 40, 50, 60} 因此状态变量可能的取值可以描述为： s = nan ∗ ones(7, 4); s(1, 1) = 60; s(:, 2) = [0 10 20 30 40 50 60] ; s(:, 3) = [0 10 20 30 40 50 60] ; s(:, 4) = [0 10 20 30 40 50 60] ;
16 / 28
利用MATLAB求解离散动态规划问题
三、动态规划求解函数-算例1求解过程分析（续）
4. 状态转移方程Tk ： 显然有：sk+1 = sk − u 因此可写出以下状态转移函数： function s next = TransFun(k, s, u) s next = s − u;
17 / 28
在 开 始 编 程 求 解 之 前 需 要 先 简 单 分 析 一下 题 目的 情 景 ， 以 便 简 化求解流程。稍加分析就能发现，其实这个问题只需要做三次决 策 ， 分别 在A、B和C阶 段 ， 行 进 到D阶 段 以 后 下一 步 自 然 就 只 能 到达E点，也就是终点。 根据以上分析以及最短路问题的特点，将各结点用数字标明如下 图：
9 / 28
利用MATLAB求解离散动态规划问题
二、穷举法：算例2-算法与程序
在上面的带权的状态转移矩阵中，若两个结点间不能相互到达或 按着题给示意图的方向到达则留空，若能够到达，则表格中的数 值代表弧长。注意这里的到达关系是单项的，只能顺着题给示意 图的方向行进。 对应上面的状态转移矩阵，在MATLAB中对应初始化，其中inf代 表无穷大数，表示两个结点之间不能相互到达或不能按题给示意 图方向到达：
11 / 28
利用MATLAB求解离散动态规划问题
二、穷举法：算例2-运行结果
最后照例在穷举结果中找最优值： Opt D = min(D(1 : end)) Opt Path = Path(find(D == Opt D), :) 运行结果： Opt D = 14 Opt Path = [1, 2, 5, 8, 11] 翻译过来就是 A—B1—C2—D1—E， 这与我们传统求解方法得到 的结果相一致。整个过程依旧非常迅速，耗时约0.008秒。最后还 顺带得到所有可能的路径为12条。
利用MATLAB求解离散动态规划问题
利用MATLAB求解离散动态规划问题
陈俊先
中山大学 岭南学院
2016年8月30日
1 / 28
利用MATLAB求解离散动态规划问题
一、简介
目前MATLAB内置的最优化工具包不能很好的处理多期最优化问 题，在求解动态规划问题上并不非常适用
目标：利用MATLAB实现离散动态规划问题的求解，寻求并设计 编写一种通用的求解算法和相应的函数
5 / 28
利用MATLAB求解离散动态规划问题
二、穷举法：算例1-运行结果
得到最优值与相应的最优策略： Opt V = max(J (1 : end)) Opt U = (U (find(J == Opt V ), :) − 1) ∗ 10 运行结果： Opt V = 160 Opt U = [20, 0, 30, 10] 由于计算量很小，整个过程的运算速度非常快，从MATLAB自带 的计时器来看只用了大约0.009秒。
8 / 28
利用MATLAB求解离散动态规划问题
二、穷举法：算例2
首先根据最短路问题中状态转移的特点，写出带权的状态转移矩 阵如下：
1 1 2 3 4 5 6 7 8 9 10 11 2 5 3 3 4 1 5 3 8 6 6 7 7 8 9 10 11
6 6 3 8 5 3 8
3 4 3 3 2