混合型非線性規(guī)劃的MATLAB實(shí)現(xiàn)方法

沈文勝 熊方方 金升平 余后強(qiáng)

（武漢理工大學(xué)理學(xué)院 武漢 430070）

0 引 言

最優(yōu)化問(wèn)題廣泛應(yīng)用于工程應(yīng)用的各部門中.在實(shí)際應(yīng)用中經(jīng)常需要求解混合規(guī)劃問(wèn)題，對(duì)于整數(shù)線性規(guī)劃問(wèn)題，文獻(xiàn)［1］給出了其 MATLAB（本文均指 MATLAB7.0）的實(shí)現(xiàn)，而對(duì)于求解一般的非線性混合整數(shù)規(guī)劃的問(wèn)題時(shí)，在MATLAB中可以調(diào)用一個(gè)基于分枝定界法編寫(xiě)的現(xiàn)成函數(shù)bnb20（）［2］.但是對(duì)于有既含有整數(shù)又含有離散型的混合非線性規(guī)劃問(wèn)題，例如在機(jī)械工程的齒輪設(shè)計(jì)中涉及到的齒輪的模數(shù)以及齒數(shù)時(shí)，模數(shù)就是統(tǒng)一的標(biāo)準(zhǔn)數(shù)值（如系列（GB1357－87））而齒數(shù)則必須是整數(shù)的問(wèn)題，現(xiàn)在的一些數(shù)學(xué)軟件也沒(méi)有可以解決此類問(wèn)題的方法，本文給出一個(gè)在MATLAB中求解此類問(wèn)題的實(shí)現(xiàn)方法.

考慮一般的混合型非線性規(guī)劃問(wèn)題如式（1）所示.

式中：x ＝（x1，x2，…，xn）T為自變量；UB，LB 為向量x的上下界；Aeq，beq分別為線性約束條件中等式約束條件的系數(shù)矩陣和常數(shù)項(xiàng)；A，b分別為線性約束條件中不等式約束條件的系數(shù)矩陣和常數(shù)項(xiàng)；Ceq（x）和C（x）分別為非線性約束條件中等式約束條件和不等式約束條件；xi1，xi2，…，xir為離散型變量，取值范圍分別為Di1，Di2，…，Dir，xir＋1，…，xik為整數(shù)變量；其余變量均為連續(xù)變量.

分枝定界算法是一種在問(wèn)題的解空間樹(shù)（二叉樹(shù)）上搜索問(wèn)題的解的方法.本文利用分枝定界算法對(duì)問(wèn)題的解空間樹(shù)進(jìn)行搜索，所采用的分枝定界法的思想是：將原始問(wèn)題分解，產(chǎn)生一組子問(wèn)題，分枝（branching）是將一組解分為幾組子解，定界（bounding）是建立這些子組解的目標(biāo)函數(shù)的邊界，如果某一子組的解在這些邊界之外，就將這一子組舍棄（剪枝（prune））.

1 MATLAB實(shí)現(xiàn)中的技術(shù)處理

搜索算法按搜索的方式分有兩類：一類是深度優(yōu)先搜索（depth first search，DFS），其耗費(fèi)的時(shí)間取決于所采用的存儲(chǔ)結(jié)構(gòu).當(dāng)用二維數(shù)組表示鄰接矩陣作圖的存儲(chǔ)結(jié)構(gòu)時(shí)，查找每個(gè)頂點(diǎn)的鄰接點(diǎn)所需時(shí)間為O（n2），其中n為圖中頂點(diǎn)數(shù).而當(dāng)以鄰接表作圖的存儲(chǔ)結(jié)構(gòu)時(shí)，查找鄰接點(diǎn)所需的時(shí)間為O（e），其中e為無(wú)向圖中邊的數(shù)或有向圖中弧的數(shù).由此，當(dāng)以鄰接表作存儲(chǔ)結(jié)構(gòu)時(shí)，深度優(yōu)先搜索遍歷圖的時(shí)間復(fù)雜度O（n＋e）；另一類是廣度優(yōu)先搜索（breadth first searchm，BFS）是樹(shù)的按層次遍歷的推廣，遍歷圖的過(guò)程實(shí)質(zhì)上是通過(guò)邊或弧找鄰接點(diǎn)過(guò)程，因此廣度優(yōu)先搜索遍歷圖的時(shí)間復(fù)雜性和深度優(yōu)先搜索遍歷相同，兩者不同之處僅僅在于對(duì)頂點(diǎn)訪問(wèn)的順序不同［3］.

在MATLAB實(shí)現(xiàn)過(guò)程中不需要保留樹(shù)的結(jié)構(gòu)，所以本文利用了兩種常見(jiàn)的數(shù)據(jù)結(jié)構(gòu)來(lái)解決變量存儲(chǔ)問(wèn)題：對(duì)于深度優(yōu)先用棧（stack）數(shù)據(jù)結(jié)構(gòu)來(lái)做其存儲(chǔ)結(jié)構(gòu)，這是利用棧是按照后進(jìn)先出的原則存儲(chǔ)數(shù)據(jù)的性質(zhì)；對(duì)于廣度優(yōu)先法則用隊(duì)列（queue）來(lái)表示，這是利用它只允許在表的前端（front）進(jìn)行刪除操作，而在表的后端（rear）進(jìn)行插入操作的特性.再結(jié)合MATLAB的特點(diǎn)，在實(shí)現(xiàn)過(guò)程中來(lái)節(jié)約時(shí)間和空間開(kāi)銷.文獻(xiàn)［4］中例子的搜索樹(shù)如圖1所示，以此樹(shù)為例介紹DFS和BFS兩種遍歷過(guò)程中變量的存儲(chǔ)情況.首先介紹DFS這種情況見(jiàn)表1.其初始狀態(tài)為P0，但是P0中的要左分枝即第1步；P1滿足約束條件是葉節(jié)點(diǎn)需要剪枝，在表1中就是刪除P1，剪枝后P0右分枝得到第2步；P2中的解不滿足約束條件，再左分枝得到P3表1第3步；P3中的解不滿足約束條件，繼續(xù)左分得到P4表1第4步；P4滿足約束條件，P4小于P1，所以P4是葉節(jié)點(diǎn)需要剪枝，剪枝后返回對(duì)P3右分枝得到P5（表1第5步），其他的以此類推.

圖1 例題的搜索樹(shù)

表1 DFS遍歷過(guò)程

BFS具體的實(shí)際情況見(jiàn)表2.初始狀態(tài)是相同，但是分之后存儲(chǔ)的方式不同，每次分枝后其原來(lái)節(jié)點(diǎn)就被刪除，被兩個(gè)分枝節(jié)點(diǎn)所代替，即P0分枝后得到表2第1步，P1葉節(jié)點(diǎn)需要?jiǎng)h除，得到第2步；P2分枝得到第3步；P3分枝后得到第4步.其他的以此類推.

表2 BFS遍歷過(guò)程

由此可見(jiàn)經(jīng)過(guò)這種處理，極大的減少了計(jì)算機(jī)的存儲(chǔ)空間.尤其對(duì)于運(yùn)用在大規(guī)模計(jì)算中，其優(yōu)點(diǎn)是顯然的.

對(duì)于線性問(wèn)題，文獻(xiàn)［4］中認(rèn)為DFS方法優(yōu)于BFS方法，并且給出了3點(diǎn)理由.基于類似的原因，對(duì)于非線性問(wèn)題，建議使用DFS方法，第二部分的設(shè)計(jì)方法就是以DFS為基礎(chǔ)編制的，也可以用BFS來(lái)編寫(xiě)，原則上只是訪問(wèn)的順序不同.

2 基于MATLAB的設(shè)計(jì)方法

用一個(gè)3×k的數(shù)組xstatus來(lái)設(shè)定離散或者整數(shù)變量的狀態(tài)，其中xstatus（1；：）是自變量中離散或整數(shù)變量的下標(biāo)；xstatus（2；：）中的每個(gè)元素取值為1或2，1表示xstatus（1；：）中對(duì)應(yīng)列的自變量為整數(shù)型變量，2表示xstatus（1；：）中對(duì)應(yīng)列的變量為離散變量；xstatus（3；：）是xstatus（1；：）中對(duì)應(yīng)列的離散變量的取值范圍的序號(hào)，若變量為整數(shù)則為零.于是xstatus就可以表示為

構(gòu)造節(jié)點(diǎn)的數(shù)據(jù)結(jié)構(gòu)如下.

其中：‘Left＿or＿Right’表示左右節(jié)點(diǎn)，‘branch＿variable＿index’表示分枝變量的下標(biāo)，‘xlowe’和‘xupper’分別表示分枝變量的下界和上界，‘solution’表示解向量，‘obj＿value’表示目標(biāo)函數(shù)的值，‘LB’和‘UB’分別表示該節(jié)點(diǎn)中變量的下界和上界.分枝的過(guò)程就是不斷地修改‘LB’和‘UB’的過(guò)程.

對(duì)于此類問(wèn)題所采取的策略是：（1）設(shè)定當(dāng)前最好解 MIDP（1）.設(shè)定目標(biāo)函數(shù) MIDP（1）.obj＿value＝inf；（2）利用函數(shù)FMINCON（）計(jì)算當(dāng)前節(jié)點(diǎn)松弛后的解.使用MATLAB自帶函數(shù)FMINCON（）進(jìn)行求解，若EXITFLAG的返回值小于等于零則求解出現(xiàn)錯(cuò)誤，該節(jié)點(diǎn)需要剪枝.若當(dāng)前的目標(biāo)函數(shù)值大于 MIDP（1）.obj＿value，則該節(jié)點(diǎn)需剪枝；（3）判斷是否符合離散和整數(shù)（checkIntDisc）.若是整數(shù)變量則執(zhí)行判斷該變量是否在某個(gè)整數(shù)的δ（δ是計(jì)算時(shí)允許的誤差限，根據(jù)實(shí)際問(wèn)題選取一較小的值）鄰域內(nèi)，若在該鄰域內(nèi)則說(shuō)明該值是一可行解，返回到（3）繼續(xù)判斷下一個(gè)變量；若不在該鄰域內(nèi)，則對(duì)該變量做向上舍入（ceil）和向下舍入（floor），且將舍入后的值記為該變量的下一次分枝的區(qū)間的上界和下界，返回該值.結(jié)束此次判斷.若是離散變量則將該變量的所有取值范圍與此時(shí)的值進(jìn)行判斷，判斷此變量是否在其取值范圍中的某個(gè)值的δ鄰域內(nèi)，若在該鄰域內(nèi)則說(shuō)明該值是一可行解，返回到（3）繼續(xù)判斷下一個(gè)變量.否則得到該變量的離散取值范圍中的上下值作為下一次分枝的上界和下界；（4）對(duì)（3）的結(jié)果進(jìn)行判斷.若符合整數(shù)和離散的要求，則更新 MIDP（1），此時(shí)該節(jié)點(diǎn)需剪枝；否則執(zhí)行（5）；（5）對(duì)于節(jié)點(diǎn)（node）進(jìn)行左分枝（left）或者右分枝（right）.設(shè)定分枝開(kāi)關(guān)，當(dāng) MIDP（index）.Left＿or＿right＝1時(shí)向左進(jìn)行分枝，當(dāng)MIDP（index）.Left＿or＿right＝2時(shí)進(jìn)行右分枝.具體方法如下：設(shè)定MIDP（index）.Left＿or＿right＝1時(shí)更新自變量的取值區(qū)間，將上一節(jié)點(diǎn)取值區(qū)間的下界記為這一節(jié)點(diǎn)的取值區(qū)間的上界；否則 MIDP（index）.Left＿or＿right＝2向右進(jìn)行分枝，將上一節(jié)點(diǎn)取值區(qū)間的上界記為這一節(jié)點(diǎn)的取值區(qū)間的下界，返回（2）.

3 應(yīng)用舉例

3.1 一個(gè)實(shí)例

給定一個(gè)求解某齒輪的問(wèn)題，自變量為x＝［x1，x2，x3，x4，x5，x6，x7，x8，x9，x10，x11，x12，x13，x14］.式中：x1，x2，x3為模數(shù)且取值只能為D1中的數(shù)；x4，x5，x6，x7，x8，x9為齒輪的齒數(shù)，故取值只能為整數(shù)；x10，x11為螺旋角的弧度；x12，x13，x14為嚙合齒齒寬.并且x的上下界和D1分別為選取本問(wèn)題的初始值為

在MATLAB中運(yùn)用DFS編制的程序來(lái)計(jì)算，計(jì)算過(guò)程如下：在MATLAB中輸入A＝［］；b＝［］；Aeq＝［］；beq＝［］；xstatus＝［1，2，3，4，5，6，7，8，9；2，2，2，1，1，1，1，1，1；1，1，1，0，0，0，0，0，0］；以及LB，UB，D1，X0后，運(yùn)行本方法所編軟件，經(jīng)過(guò)約4 min的計(jì)算，得出解向量和目標(biāo)函數(shù)值分別為：［2.75，2.5，6，28，85，11，74，12，53，0.254 8，0.254 8，31.052 5，31.157 8，74.778 7］和1.157 8×107.

3.2 應(yīng)用于文獻(xiàn)［5］問(wèn)題

對(duì)于優(yōu)化問(wèn)題，工程人員常用的方法是把先忽略自變量的離散或整數(shù)約束要求，當(dāng)作連續(xù)型的非線性規(guī)劃問(wèn)題進(jìn)行求解，然后將得到的結(jié)果圓整［5］，即將自變量選取和它最近的離散值或整數(shù)值.這在數(shù)學(xué)中是沒(méi)有理論依據(jù)的，進(jìn)行圓整所得到的結(jié)果很有可能不滿足約束條件.

對(duì)于文獻(xiàn)［5］中的目標(biāo)函數(shù)及其約束條件，運(yùn)用本文方法和所開(kāi)發(fā)的程序進(jìn)行求解，設(shè)置參數(shù)如下：D1＝［2，2.5，3，4，5］；D2＝［3，4，5，6］；xstatus＝［1，2，3，4；2，2，1，1；1，2，0，0］.可得到如下解向量和目標(biāo)函數(shù)的值：解向量為［2，4，19，17，5.816 9，0.990 3］；目標(biāo)函數(shù)值為351.054 7.

文獻(xiàn)［5］的解向量和目標(biāo)函數(shù)的值分別是［2.005 0，4.004 2，14.001 7，16.376 6，5.8，0.99］和346.259，其圓整后結(jié)果解向量為［2，4，15，17，6，0.990 3］，目標(biāo)函數(shù)值是350.

經(jīng)過(guò)實(shí)際驗(yàn)算，還發(fā)現(xiàn)文獻(xiàn)［5－7］圓整后的解向量不滿足第一和第二個(gè)約束條件，是不可行解，說(shuō)明其圓整的方法對(duì)所述問(wèn)題是行不通的.而本文的方法勿需"圓整"，簡(jiǎn)單易行，且求出的解為最優(yōu)解.

4 結(jié)束語(yǔ)

對(duì)于變量中既含離散約束、又含整數(shù)約束的混合非線性規(guī)劃問(wèn)題，本文依據(jù)分枝定界原理給出了MATLAB中的實(shí)現(xiàn)方法，在其實(shí)現(xiàn)過(guò)程中，在樹(shù)的遍歷上運(yùn)用了兩種存儲(chǔ)技巧，減少了運(yùn)行時(shí)間和存儲(chǔ)空間，開(kāi)發(fā)了相應(yīng)的軟件，實(shí)例表明方法的正確性和軟件的可行性.

［1］潘 君.整數(shù)規(guī)劃的分枝定界法及其MATLAB實(shí)現(xiàn)［J］.科技信息，2008，（7）：167－168.

［2］薛定宇，陳陽(yáng)泉.高等應(yīng)用數(shù)學(xué)問(wèn)題的 MATALAB求解［M］.北京：清華大學(xué)出版社，2004.

［3］嚴(yán)蔚敏，吳偉民.數(shù)據(jù)結(jié)構(gòu)［M］.北京：清華大學(xué)出版社，1997.

［4］Vanderbei R J.Linear programming：foundations and extensions［M］.Boston：Kluwer Academic Publisher，2001.

［5］崔樹(shù)平，趙 亮，崔 涵.基于 MATLAB最小體積齒輪減速器的優(yōu)化設(shè)計(jì)［J］.機(jī)械管理開(kāi)發(fā)，2008，23（6）：54－55.

［6］金全意，孟 航.基于 MATLAB的圓柱齒輪減速器優(yōu)化設(shè)計(jì)［J］.遼寧工程技術(shù)大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，29（z1）：55－59.

［7］張慧鵬.基于MATLAB的二級(jí)圓柱齒輪減速器優(yōu)化設(shè)計(jì) ［J］.Machinery Design ＆ Manufacture，2010（4）：101－106.