沈文勝 熊方方 金升平 余后強

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

0 引 言

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

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

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

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

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

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

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

圖1 例題的搜索樹

表1 DFS遍歷過程

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

表2 BFS遍歷過程

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

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

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

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

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

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

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

3 應(yīng)用舉例

3.1 一個實例

給定一個求解某齒輪的問題，自變量為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分別為選取本問題的初始值為

在MATLAB中運用DFS編制的程序來計算，計算過程如下：在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后，運行本方法所編軟件，經(jīng)過約4 min的計算，得出解向量和目標(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)用于文獻［5］問題

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

對于文獻［5］中的目標(biāo)函數(shù)及其約束條件，運用本文方法和所開發(fā)的程序進行求解，設(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.

文獻［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)過實際驗算，還發(fā)現(xiàn)文獻［5－7］圓整后的解向量不滿足第一和第二個約束條件，是不可行解，說明其圓整的方法對所述問題是行不通的.而本文的方法勿需"圓整"，簡單易行，且求出的解為最優(yōu)解.

4 結(jié)束語

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

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

［2］薛定宇，陳陽泉.高等應(yīng)用數(shù)學(xué)問題的 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］崔樹平，趙 亮，崔 涵.基于 MATLAB最小體積齒輪減速器的優(yōu)化設(shè)計［J］.機械管理開發(fā)，2008，23（6）：54－55.

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

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