熊文濤， 雍龍泉

(1.湖北工程學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院， 湖北 孝感 432000；2.陜西理工學(xué)院 數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院， 陜西 漢中 723001)

整數(shù)規(guī)劃模型是運(yùn)籌學(xué)中一類常見的數(shù)學(xué)模型，在現(xiàn)實(shí)生活中有著廣泛的應(yīng)用。在許多仿真實(shí)驗(yàn)中，求解整數(shù)規(guī)劃最常見的工具有Matlab和Lingo軟件。然而，Matlab軟件并沒有提供直接求解一般整數(shù)規(guī)劃的內(nèi)部函數(shù)。盡管它提供了求解0-1規(guī)劃的bintprog函數(shù)，但其輸入的約束條件要求采用矩陣乘積形式，這對于人們習(xí)慣用求和或任意形式構(gòu)造成的模型，不便于轉(zhuǎn)換成標(biāo)準(zhǔn)的矩陣乘積形式，求解起來非常費(fèi)時(shí)。雖然Lingo軟件能根據(jù)人們計(jì)算習(xí)慣將求解一般整數(shù)規(guī)劃問題的模型轉(zhuǎn)換成程序語句，但無法自動(dòng)多次求解整數(shù)規(guī)劃，僅限于每次求解一個(gè)優(yōu)化模型。為方便整數(shù)優(yōu)化模型求解，本文將介紹求解最優(yōu)化模型的Yalmip工具箱。Yalmip是在Matlab軟件的仿真平臺(tái)上，采用Matlab的程序設(shè)計(jì)語言開發(fā)的一個(gè)工具箱。它不僅可以利用Matlab強(qiáng)大的計(jì)算能力，調(diào)用Matlab軟件自帶的優(yōu)化工具箱和Matlab語言編寫的其他工具箱中的函數(shù)，而且還能將多種優(yōu)化求解器(如Cplex和Gurobi等工具)結(jié)合在一起使用。本文闡述了利用Yalmip工具箱求解整數(shù)規(guī)劃的一般方法，并通過一個(gè)具體的實(shí)例說明使用Yalmip工具箱在求解整數(shù)規(guī)劃模型的一般方法。

1 Yalmip工具箱簡介

Yalmip是由Lofberg開發(fā)的一種免費(fèi)的Matlab工具箱[1]，其最初目的是為求解控制與系統(tǒng)理論中的半定規(guī)劃和線性矩陣不等式而設(shè)計(jì)的。后來經(jīng)過不斷發(fā)展，可以求解更多的優(yōu)化模型，如二階錐規(guī)劃、混合整數(shù)規(guī)劃、正項(xiàng)幾何規(guī)劃、雙線性矩陣不等式以及多參數(shù)線性規(guī)劃和二次規(guī)劃等。Yalmip最大的特色在于能集成許多外部的最優(yōu)化求解器，無論這些求解器是否采用Matlab程序語言編寫的。通常，這些求解器把優(yōu)化問題描述成一個(gè)非常緊湊的格式，學(xué)習(xí)、使用起來比較耗時(shí)，并且編程時(shí)容易出錯(cuò)。為了克服這一方面的不足，Yalmip提供了一種合適的建模語言和編程接口，以方便不同求解器之間的集成。Yalmip也能調(diào)用Matlab自帶的優(yōu)化工具箱，特別是對于未知變量較多的線性規(guī)劃問題，調(diào)用Linprog函數(shù)不需要轉(zhuǎn)換成矩陣的形式，使用起來非常方便。Yalmip也包含內(nèi)置的整數(shù)規(guī)劃求解器，可以求解一些小型的整數(shù)規(guī)劃問題，并且與外部的整數(shù)規(guī)劃求解器具有較好的兼容性。

Yalmip作為一個(gè)Matlab的工具箱，與一般的外部工具箱安裝相同。使用者可以通過該工具箱的網(wǎng)址http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Tutorials.Installation下載安裝包并解壓后，可通過如下步驟將對應(yīng)的文件夾放到Matlab路徑里：首先點(diǎn)擊File中的Set Path菜單，然后點(diǎn)擊Add with subfolders菜單，找到解壓好的文件夾Yalmip，最后確定、保存即可。對于外部的一些優(yōu)化求解器，可采用類似的方法放入到Matlab路徑中。

2 整數(shù)規(guī)劃模型及其Yalmip工具箱的程序?qū)崿F(xiàn)

一般的整數(shù)規(guī)劃可表示為：

(1)

其中，f(x)為目標(biāo)函數(shù)，gi(x)和hj(x)是定義在Rn上的多元實(shí)值函數(shù)，決策變量x至少有一個(gè)分量為整數(shù)。當(dāng)x全部為整數(shù)時(shí)，模型(1)為純整數(shù)規(guī)劃。

特別地，當(dāng)f(x),gi(x),hj(x)均為線性函數(shù)時(shí)，此模型(1)稱為整數(shù)線性規(guī)劃，模型(1)可重新改寫為如下的矩陣表示形式：

(2)

其中c是n維列向量，A和Aeq分別是m×n矩陣和l×n矩陣，b是m維列向量，beq是l維列向量。

對于許多實(shí)際問題，人們習(xí)慣于用求和形式與任意的符號(hào)形式進(jìn)行表示，而不是矩陣的乘積形式，于是模型(2)又可寫為模型(3)的形式，即：

(3)

其中ci，bi，beqs分別向量c，b，beq的分量，aij，aeqsj則分別A,Aeq是矩陣中元素。

為了方便，假設(shè)模型(3)中所有未知變量 全部為整數(shù)。若部分變量不是整數(shù)，則只需改變申明變量類型的函數(shù)即可。表1給出了求解整數(shù)線性規(guī)劃的一般程序。

表1 求解整數(shù)線性規(guī)劃的一般程序模式

3 實(shí)例求解

3.1 問題描述

以文獻(xiàn)[2]中一個(gè)電力生產(chǎn)問題作為實(shí)例，介紹使用Yalmip工具箱求解整數(shù)規(guī)劃模型的求解方法。為滿足每日電力需求(單位為MW，見表2)，有4種不同類型的發(fā)電機(jī)可供選用。所有發(fā)電機(jī)都存在一個(gè)啟動(dòng)成本(在每個(gè)時(shí)段開始時(shí)才允許啟動(dòng)或關(guān)閉)，以及工作于最小功率狀態(tài)時(shí)的固定成本，并且如果功率高于最小功率，則超出部分的功率每兆瓦每小時(shí)還存在一個(gè)邊際成本。另外，每種發(fā)電機(jī)都有一個(gè)最大發(fā)電能力，當(dāng)接入電網(wǎng)時(shí)，其輸出功率不應(yīng)低于某一最小輸出功率，具體數(shù)據(jù)見表3。與啟動(dòng)發(fā)電機(jī)不同，關(guān)閉發(fā)電機(jī)不需要付出任何代價(jià)。在任意時(shí)刻，正在工作的發(fā)電機(jī)組必須留出20%的發(fā)電能力余量，以防用電量突然上升。每個(gè)時(shí)段應(yīng)分別使用哪些發(fā)電機(jī)才能使每天的總成本最??？

表2 每日用電需求(MW)

表3 發(fā)電機(jī)情況

3.2 模型建立與求解

假設(shè)P為發(fā)電機(jī)型號(hào)集合，T={1，2，…，NT}為每天的時(shí)段集合，Lenj,Demj分別表示第j個(gè)時(shí)段的長度(小時(shí)數(shù))和電網(wǎng)的用電量需求，Pmini，Pmaxi,Availi,Cstarti,Cmini,Caddi為i型發(fā)電機(jī)的最小功率、最大功率、可用的發(fā)電機(jī)數(shù)目、啟動(dòng)成本、最小功率每小時(shí)工作成本、以及每兆每小時(shí)的邊際成本。未知變量startij，workij及paddij分別表示在時(shí)段j開始工作的型號(hào)i的發(fā)電機(jī)數(shù)量，在時(shí)段j內(nèi)工作型號(hào)i的發(fā)電機(jī)數(shù)量，以及型號(hào)i的發(fā)電機(jī)在時(shí)段j內(nèi)功率超出其最低功率之上的量。根據(jù)上述問題描述可以得到如下的整數(shù)規(guī)劃模型。

(4)

其中，約束條件(4.1)表示對于每個(gè)類型和每個(gè)時(shí)段的發(fā)電機(jī)(型號(hào)為i)，其額外功率都不能超過Pmaxi-Pmini；約束條件(4.2)表示每個(gè)時(shí)段內(nèi)發(fā)電量都應(yīng)滿足用電需求；約束條件(4.3)則表示在不啟動(dòng)新的發(fā)電機(jī)的前提下，當(dāng)前正在工作的發(fā)電機(jī)的發(fā)電能力必須高于實(shí)際用電量的20%；另外，根據(jù)需求，當(dāng)不關(guān)閉發(fā)電機(jī)時(shí)，在開始時(shí)段時(shí)，投入開始工作的發(fā)電機(jī)數(shù)目應(yīng)等于時(shí)段j和j-1間工作的發(fā)電機(jī)數(shù)量之差。如果有些發(fā)電機(jī)可能會(huì)被關(guān)閉，那么開始工作的發(fā)電機(jī)數(shù)量應(yīng)至少等于此差值，于是得到約束條件(4.4)；約束條件 (4.5)表示當(dāng)日午夜和次日凌晨這兩個(gè)時(shí)段的情形。

針對此模型，可以通過表4的程序語句計(jì)算出其最優(yōu)解和最優(yōu)值，得到發(fā)電機(jī)的使用計(jì)劃，其結(jié)果見表5。在Windows xp操作系統(tǒng)，3.20 GHz Intel Core i5 CPU和 4 GB RAM的PC上求解此整數(shù)規(guī)劃，時(shí)間僅用了0.3 s。表5的結(jié)果表明：對于型號(hào)1的發(fā)電機(jī)，時(shí)段1時(shí)有3臺(tái)在工作，時(shí)段2開始時(shí)啟動(dòng)1臺(tái)，時(shí)段4開始時(shí)啟動(dòng)3臺(tái)，時(shí)段4之后關(guān)閉4臺(tái),4臺(tái)型號(hào)2的發(fā)電機(jī)都將連續(xù)工作；對于型號(hào)3的發(fā)電機(jī)，時(shí)段1時(shí)有2臺(tái)在工作，時(shí)段2 開始時(shí)啟動(dòng)6臺(tái)，時(shí)段6結(jié)束后關(guān)閉4臺(tái)，一天結(jié)束時(shí)再關(guān)閉2臺(tái)；對于型號(hào)4的發(fā)電機(jī)，時(shí)段2開始時(shí)啟動(dòng)3臺(tái)，到下一時(shí)段時(shí)關(guān)閉2臺(tái)，然后在下面的幾個(gè)時(shí)段這2臺(tái)反復(fù)開閉，

表4 求解模型(4)的Matlab程序

直到一天結(jié)束時(shí)，所有發(fā)電機(jī)都關(guān)閉。每種型號(hào)的發(fā)電機(jī)在每個(gè)時(shí)段的總功率可以通過發(fā)電機(jī)使用數(shù)量、最小輸出功率和額外功率計(jì)算，如型號(hào)1的發(fā)電機(jī)在第2時(shí)段的總功率為4600 MW。

4 結(jié)語

整數(shù)規(guī)劃模型是運(yùn)籌學(xué)中一類常見的優(yōu)化求解問題，在實(shí)際應(yīng)用中有著非常廣泛的應(yīng)用。由于Matlab軟件并沒有提供直接求解一般整數(shù)規(guī)劃的內(nèi)部函數(shù)，不便于人們直接對一般整數(shù)規(guī)劃優(yōu)化求解問題。為方便整數(shù)規(guī)劃問題求解，本文利用Matlab平臺(tái)，介紹了一種使用Yalmip工具箱求解整數(shù)規(guī)劃問題的一般方法。該方法能利用Matlab的強(qiáng)大計(jì)算能力，以Matlab為平臺(tái)調(diào)用Yalmip工具箱對優(yōu)化問題進(jìn)行求解，編程方式具有靈活、直觀等優(yōu)點(diǎn)。而且利用Yalmip能夠集成多種優(yōu)化求解器的優(yōu)點(diǎn)，該方法能方便求解其他類型的優(yōu)化問題，如線性規(guī)劃、二次規(guī)劃和一般的非線性規(guī)劃，便于對多個(gè)求解器計(jì)算的結(jié)果進(jìn)行比較和選擇。以電力生產(chǎn)問題作為實(shí)例，說明了使用Yalmip工具箱求解整數(shù)規(guī)劃模型的過程。計(jì)算結(jié)果表明了本文方法的有效性。

表5 發(fā)電機(jī)的使用計(jì)劃

[參 考 文 獻(xiàn)]

[1] Lofberg J. YALMIP: A toolbox for modeling and optimization in MATLAB[C]//Proceedings of 2004 IEEE International Symposium on Computer Aided Control Systems Design, 2004:284-289.

[2] Gueret C,Prins C,Sevaux M.運(yùn)籌學(xué)案例[M].北京: 北京林森科技發(fā)展有限公司,2006.