MATLAB在水資源動態(tài)規(guī)劃中的運用與改進

駱歡　劉培　李思晨

【摘 要】水資源是人類賴以生存的寶貴資源，搞好水資源優(yōu)化配置對實現水資源的可持續(xù)利用尤為重要。本文結合一個實例，詳細介紹了用MATLAB解決水資源動態(tài)規(guī)劃的思路與方法，同時與傳統(tǒng)的逆序法對比，突出了MATLAB的優(yōu)越性。并且，參照了之前學者所編寫的程序，給出了自己的程序優(yōu)化方案，使程序運行速度、可讀性等方面得到了提高。

【關鍵詞】動態(tài)規(guī)劃；水資源優(yōu)化配置；MATLAB

0.引言

動態(tài)規(guī)劃是解決多階段決策過程最優(yōu)化問題的一種方法，該方法是由美國數學家貝爾曼（R.Bellman）等人在20世紀50年代提出。他們針對多階段決策問題的特點，提出了解決這類問題的最優(yōu)化原理，并成功解決了生產管理、資源分配等方面的許多實際問題，從而建立了運籌學的分支——動態(tài)規(guī)劃[1]。MATLAB是一個功能強大的用于矩陣運算的數值計算軟件，是決策系統(tǒng)優(yōu)化計算和設計的有力工具，將MATLAB運用到動態(tài)規(guī)劃中可以達到計算簡便的目的[2]?，F在利用MATLAB解決動態(tài)規(guī)劃問題的程序不少，但是，如何充分利用MATLAB的特點來優(yōu)化程序確實一個值得深思的問題。

1.MATLAB程序設計

MATLAB有2種工作方式：一是交互式的命令工作方式，直接在Command Window中輸入指令求解。另一種是M腳本文件的工作方式，這種工作方式用來解決指令較多和所用指令結構較為復雜的問題[3]。

由于動態(tài)規(guī)劃問題的特殊性，其涉及的實際問題均需要反復的計算求解，所以在使用MATLAB語言編程中采用最多的是循環(huán)結構。MATLAB提供了2種實現循環(huán)結構的語句——for語句和while語句[4]。

2.應用實例

有一引水渠，其設計最大流量為6m3/s，為甲、乙、丙、丁4個地區(qū)供水。據統(tǒng)計，在給四個地區(qū)供水時，不同的供水量產生的效益見表1（效益單位為萬元）。那么，如何分配供水量可以使產生的總效益最大。（本實例取自文獻[5]）

表1 供水量與增產效益的關系

2.1逆序法計算

采用逆序法表格計算，可以得到最優(yōu)分配方案有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個地區(qū)不同水量時產生的效益?？紤]到在動態(tài)規(guī)劃時，某個城市的供水量可以為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.找出所有可行解，儲存其分配方法和對應的效益，要找出所有的可行解，最簡便的的方法是采用循環(huán)結構枚舉。

①確定循環(huán)條件：

循環(huán)的限制條件應該是最大設計流量。由于最大設計流量是6 m3/s，現在任意假設一例分配方法以確定循環(huán)條件。假設所有的水量全部供給了丁，即l=7，那么此時可有i=j=k=1，得i+j+k+l=10。由于無論是那種分配方法都要滿足最大設計流量的限制，所以i+j+k+l=10即為循環(huán)條件。由于循環(huán)次數是已知的，所以采用for語句循環(huán)。

②儲存分配方法：

建立一個二維數組E，采用全下標的方法將分配方法存入其中，每一行為一種分配方式。

③儲存效益：

建立一個數組d，常用單下標的方式將所有分配方法對應的效益存入其中。

Step3.確定最大效益的值

所有方案的效益都存儲在數組d中，所以傳統(tǒng)的方法就是采用循環(huán)結構遍歷數組d中的所有元素以找到最大值。但是在編程中應該盡量避免使用循環(huán)結構。因為在循環(huán)結構中，要反復進行存儲變量間的“存放”操作和算符調用操作，消耗計算時間。MATLAB中為用戶提供了大量的現成函數，盡可能多的采用現成函數編程可以使所編程序更可靠、更快速、可讀性更好[6]。

①采用循環(huán)結構確定最大效益

程序如下：

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函數比用循環(huán)結構更好，不僅使程序的輸入得到了簡化，同時也使程序的可讀性更好。

Step4.確定最大效益所對應的分配方法

數組E中的橫坐標和數組d中的下標是一一對應的，即數組E中橫坐標為a的行產生的總效益儲存在數組d中下標為a的元素中。所以可以采用循環(huán)遍歷的方法，找到數組d中所有值等于MAX的元素的下標n，在把數組E中對應的第n行輸出，即可得到最優(yōu)分配方法。當然，我們也可以利用matlab中現成的函數find來代替這個循環(huán)結構。

①采用循環(huán)結構確定最大值

程序如下：

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語言編程在動態(tài)規(guī)劃中的優(yōu)越性，這里分別在程序的開頭和結尾加入tic和ti=toc，以記錄下程序運行需要的時間。但應當注意的是，ti的值受到電腦配置、matlab版本和該程序是否首次運行等因素的影響。

優(yōu)化后的完整程序如下：

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

運行結果表示有三個最優(yōu)水資源分配方案，即（1，2，1，2）、（1，2，2，1）、（2，2，1，1），最大效益為19萬元，程序運行時間為0.0255 s。

2.2.3小結

從本例可以看出，運用MATLAB解決動態(tài)規(guī)劃問題比逆序法更具有可操作性。因為實際工作中要解決的問題往往比本例復雜得多，運用MATLAB編程可以避免繁瑣的人工計算。同時，相比之前學者的程序而言，優(yōu)化后的程序盡量回避了循環(huán)結構，充分利用了MATLAB中現成的函數，不僅使程序可讀性更好，輸入更簡便，還大大提高了程序運行速度。經筆者測試，若是采用沒有優(yōu)化的方案編程，程序運行時間為0.0401s，相比之下，優(yōu)化后的程序運行速度提高了36%，程序代碼也由原來的36行減少到了27行。在實際工作的運用中，特別是在復雜的實例中，這種優(yōu)點將得到更大的體現。

3.結語

本文結合一簡單實例，詳細說明了MATLAB編程的思路與方法，并給出了程序優(yōu)化方案。從本例的計算結果可以看出，運用MAYLAB編程解決動態(tài)規(guī)劃問題是實際可行的，而優(yōu)化后的程序在運算速度、可讀性等方面都有提高，這為解決實際工作中的動態(tài)規(guī)劃問題提供了新的思路。但是，如何利用MATLAB的向量化計算特點來做到更深層次的優(yōu)化還值得更多的思考與探討。 [科]

【參考文獻】

[1]孫趙杰.運籌學[M].北京：機械工業(yè)出版社，2006.

[2]趙云，王希云.基于MATLAB的動態(tài)規(guī)劃常用算法的實現[J].太原師范大學學報（自然科學版），2008，7（4）：26-30.

[3]劉衛(wèi)國.科學計算與MATLAB語言[M].北京：中國鐵道出版社，2000.

[4]樓順天.MATLAB程序與語言[M].西安：西安電子科技大學出版社，1997.

[5]左兼金.水利水電施工組織管理與系統(tǒng)分析[M].北京：水利水電出版社，1993.

[6]張志涌，楊祖櫻.MATLAB教程R2012a[M].北京：北京航空航天大學出版社，2013，129.