MATLAB在水资源动态规划中的运用与改进

骆欢　刘培　李思晨

【摘 要】水资源是人类赖以生存的宝贵资源，搞好水资源优化配置对实现水资源的可持续利用尤为重要。本文结合一个实例，详细介绍了用MATLAB解决水资源动态规划的思路与方法，同时与传统的逆序法对比，突出了MATLAB的优越性。并且，参照了之前学者所编写的程序，给出了自己的程序优化方案，使程序运行速度、可读性等方面得到了提高。

【关键词】动态规划；水资源优化配置；MATLAB

0.引言

动态规划是解决多阶段决策过程最优化问题的一种方法，该方法是由美国数学家贝尔曼（R.Bellman）等人在20世纪50年代提出。他们针对多阶段决策问题的特点，提出了解决这类问题的最优化原理，并成功解决了生产管理、资源分配等方面的许多实际问题，从而建立了运筹学的分支——动态规划[1]。MATLAB是一个功能强大的用于矩阵运算的数值计算软件，是决策系统优化计算和设计的有力工具，将MATLAB运用到动态规划中可以达到计算简便的目的[2]。现在利用MATLAB解决动态规划问题的程序不少，但是，如何充分利用MATLAB的特点来优化程序确实一个值得深思的问题。

1.MATLAB程序设计

MATLAB有2种工作方式：一是交互式的命令工作方式，直接在Command Window中输入指令求解。另一种是M脚本文件的工作方式，这种工作方式用来解决指令较多和所用指令结构较为复杂的问题[3]。

由于动态规划问题的特殊性，其涉及的实际问题均需要反复的计算求解，所以在使用MATLAB语言编程中采用最多的是循环结构。MATLAB提供了2种实现循环结构的语句——for语句和while语句[4]。

2.应用实例

有一引水渠，其设计最大流量为6m3/s，为甲、乙、丙、丁4个地区供水。据统计，在给四个地区供水时，不同的供水量产生的效益见表1（效益单位为万元）。那么，如何分配供水量可以使产生的总效益最大。（本实例取自文献[5]）

表1 供水量与增产效益的关系

2.1逆序法计算

采用逆序法表格计算，可以得到最优分配方案有3个，分别是（1，2，1，2），（1，2，2，1），（2，2，1，1），最大效益为19万元。

2.2 MATLAB编程计算

2.2.1思路分析

由于本例较为简单，所以不采用M脚本文件的工作方式，直接在Command Window中输入指令即可。

Step1.输入参数

建立A、B、C、D 4个数组，分别来储存供给甲、乙、丙、丁4个地区不同水量时产生的效益。考虑到在动态规划时，某个城市的供水量可以为0，对应的效益也为0，所以A、B、C、D四个数组都应该加入一个元素0。这里采用单下标表示法，分别用i、j、k、

l表示数组A、B、C、D中元素的下标。A（3）=2，表示当i=3时，A中元素为2，即向A供水2m3/s。Step2.找出所有可行解，储存其分配方法和对应的效益，要找出所有的可行解，最简便的的方法是采用循环结构枚举。

①确定循环条件：

循环的限制条件应该是最大设计流量。由于最大设计流量是6 m3/s，现在任意假设一例分配方法以确定循环条件。假设所有的水量全部供给了丁，即l=7，那么此时可有i=j=k=1，得i+j+k+l=10。由于无论是那种分配方法都要满足最大设计流量的限制，所以i+j+k+l=10即为循环条件。由于循环次数是已知的，所以采用for语句循环。

②储存分配方法：

建立一个二维数组E，采用全下标的方法将分配方法存入其中，每一行为一种分配方式。

③储存效益：

建立一个数组d，常用单下标的方式将所有分配方法对应的效益存入其中。

Step3.确定最大效益的值

所有方案的效益都存储在数组d中，所以传统的方法就是采用循环结构遍历数组d中的所有元素以找到最大值。但是在编程中应该尽量避免使用循环结构。因为在循环结构中，要反复进行存储变量间的“存放”操作和算符调用操作，消耗计算时间。MATLAB中为用户提供了大量的现成函数，尽可能多的采用现成函数编程可以使所编程序更可靠、更快速、可读性更好[6]。

①采用循环结构确定最大效益

程序如下：

MAX=d（1）；

for n=1：length（d）

if d（n）>MAX

MAX=d（n）；

end

end

②采用max函数确定最大效益

由于max（X）求的是数组X中各列的最大值，所以应该先对数组d转置。

程序如下：

MAX=max（d'）

可以清楚的看到，采用matlab中现成的max函数比用循环结构更好，不仅使程序的输入得到了简化，同时也使程序的可读性更好。

Step4.确定最大效益所对应的分配方法

数组E中的横坐标和数组d中的下标是一一对应的，即数组E中横坐标为a的行产生的总效益储存在数组d中下标为a的元素中。所以可以采用循环遍历的方法，找到数组d中所有值等于MAX的元素的下标n，在把数组E中对应的第n行输出，即可得到最优分配方法。当然，我们也可以利用matlab中现成的函数find来代替这个循环结构。

①采用循环结构确定最大值

程序如下：

for n=1：length（d）

if d（n）==MAX

E（n，：）-1

end

end

②采用find函数确定最大值

程序如下：

n=find（d==MAX）；

E（n，：）-1

2.2.2完整MATLAB程序

为了突出MTALAB语言编程在动态规划中的优越性，这里分别在程序的开头和结尾加入tic和ti=toc，以记录下程序运行需要的时间。但应当注意的是，ti的值受到电脑配置、matlab版本和该程序是否首次运行等因素的影响。

优化后的完整程序如下：

clear all

tic

m=1；

A=[0 3 5.5 7.5 9 10 10]；

B=[0 3 6 8 9.5 10.5 11]；

C=[0 4 6.5 8.5 9 9 9]；

D=[0 3.5 6 7.5 8.5 9 9]；

for i=1：7

for j=1：7

for k=1：7

for l=1：7

if i+j+k+l==10

d（m）=A（i）+B（j）+C（k）+D（l）；

E（m，1）=i；

E（m，2）=j；

E（m，3）=k；

E（m，4）=l；

m=m+1；

end

end

end

end

end

MAX=max（d'）

n=find（d==MAX）；

E（n，：）-1

ti=toc

按回车后可以得到以下结果：

MAX =

19

ans =

1 2 1 2

1 2 2 1

2 2 1 1

ti =

0.0255

运行结果表示有三个最优水资源分配方案，即（1，2，1，2）、（1，2，2，1）、（2，2，1，1），最大效益为19万元，程序运行时间为0.0255 s。

2.2.3小结

从本例可以看出，运用MATLAB解决动态规划问题比逆序法更具有可操作性。因为实际工作中要解决的问题往往比本例复杂得多，运用MATLAB编程可以避免繁琐的人工计算。同时，相比之前学者的程序而言，优化后的程序尽量回避了循环结构，充分利用了MATLAB中现成的函数，不仅使程序可读性更好，输入更简便，还大大提高了程序运行速度。经笔者测试，若是采用没有优化的方案编程，程序运行时间为0.0401s，相比之下，优化后的程序运行速度提高了36%，程序代码也由原来的36行减少到了27行。在实际工作的运用中，特别是在复杂的实例中，这种优点将得到更大的体现。

3.结语

本文结合一简单实例，详细说明了MATLAB编程的思路与方法，并给出了程序优化方案。从本例的计算结果可以看出，运用MAYLAB编程解决动态规划问题是实际可行的，而优化后的程序在运算速度、可读性等方面都有提高，这为解决实际工作中的动态规划问题提供了新的思路。但是，如何利用MATLAB的向量化计算特点来做到更深层次的优化还值得更多的思考与探讨。 [科]

【参考文献】

[1]孙赵杰.运筹学[M].北京：机械工业出版社，2006.

[2]赵云，王希云.基于MATLAB的动态规划常用算法的实现[J].太原师范大学学报（自然科学版），2008，7（4）：26-30.

[3]刘卫国.科学计算与MATLAB语言[M].北京：中国铁道出版社，2000.

[4]楼顺天.MATLAB程序与语言[M].西安：西安电子科技大学出版社，1997.

[5]左兼金.水利水电施工组织管理与系统分析[M].北京：水利水电出版社，1993.

[6]张志涌，杨祖樱.MATLAB教程R2012a[M].北京：北京航空航天大学出版社，2013，129.