1、本頁只是說明，論文從第二頁開始，下載后請刪除本頁即可:論文內(nèi)容：關(guān)于合理開挖土地問題的數(shù)學(xué)建模競賽論文（含Matlab源程序）特別申明：本論文版權(quán)歸百度文庫賬號dxzsk同學(xué)所有，僅限個人下載學(xué)習(xí)使用，其他人不得轉(zhuǎn)載分享，侵權(quán)必究。以下是本論文原始題目：合理開挖土地問題：A市是一個山區(qū)城市,向山要地是A市發(fā)展的一個必然的選擇,但是如何在一片山地之中選擇合適的方位與開挖深度,從而使總的土石方量最小,就是一個十分有意義的課題.A市某工廠為了在一片長度為1500米,寬度為900米的山地之中,開挖出一個800米600米平坦連續(xù)的長方形地塊作為工廠的廠房地基,前期已經(jīng)在本塊土地上測量出長、寬每隔30米的

2、網(wǎng)格的對應(yīng)網(wǎng)格點(diǎn)的海拔高度（詳細(xì)數(shù)據(jù)見附件）.請你考慮以下幾個問題：問題（1）：用附件中的數(shù)據(jù)畫出工廠的這片土地的三維圖形與等高線圖；問題（2）：從什么地方,什么海拔高度平整一塊800米600米的連片土地能使總的土石方量最小?問題（3）：如果允許平整出來的土地為二層的臺階狀地塊,要求各地塊的長、寬不少于60米,又將從什么地方、什么海拔高度分別開挖,能使總的土石方量最小?提示：在平整土地的過程中,有些地方是要挖山的,但有些地方是要填土的,假設(shè)填土的每立方米所需的費(fèi)用為挖山的每立方米土石方所需費(fèi)用的1/3.2013*大學(xué)金水節(jié)第五屆研究生數(shù)學(xué)建模競賽承 諾 書我們仔細(xì)閱讀了中國大學(xué)生數(shù)學(xué)建模競賽的

3、競賽規(guī)則。我們完全明白，在競賽開始后參賽隊(duì)員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢等）與隊(duì)外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問題。我們知道，抄襲別人的成果是違反競賽規(guī)則的, 如果引用別人的成果或其他公開的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻(xiàn)的表述方式在正文引用處和參考文獻(xiàn)中明確列出。我們鄭重承諾，嚴(yán)格遵守競賽規(guī)則，以保證競賽的公正、公平性。如有違反競賽規(guī)則的行為，我們將受到嚴(yán)肅處理。我們參賽選擇的題號是（從A/B中選擇一項(xiàng)填寫）： A 我們的參賽報(bào)名號為（如果賽區(qū)設(shè)置報(bào)名號的話）： 所屬學(xué)校（請?zhí)顚懲暾娜?參賽隊(duì)員 (打印并簽名) ：1. 2. 3.

4、指導(dǎo)教師或指導(dǎo)教師組負(fù)責(zé)人 (打印并簽名)： 日期： 2013 年 11 月 01 日賽區(qū)評閱編號（由賽區(qū)組委會評閱前進(jìn)行編號）：2013年*大學(xué)金水節(jié)第五屆研究生數(shù)學(xué)建模競賽編 號 專 用 頁賽區(qū)評閱編號（由賽區(qū)組委會評閱前進(jìn)行編號）：賽區(qū)評閱記錄（可供賽區(qū)評閱時(shí)使用）：評閱人評分備注全國統(tǒng)一編號（由賽區(qū)組委會送交全國前編號）：全國評閱編號（由全國組委會評閱前進(jìn)行編號）： 關(guān)于合理開挖土地問題的數(shù)學(xué)模型求解摘要隨著山區(qū)城市的發(fā)展，向山要地開發(fā)成為了一種必然的選擇。為了更好的確定平整地的起始點(diǎn)，做到節(jié)約成本的目的，我們根據(jù)幾個決定性因素利用網(wǎng)格優(yōu)化建立了相關(guān)的模型。主要討論如何在一片山地之中選

5、擇合適的方位與開挖深度，從而使總的土石方量最小，使得在實(shí)際工程中達(dá)到最可靠，最節(jié)約成本。采用MATLAB軟件模擬功能，畫出了該工廠的這片土地的三維圖形與等高線圖，并對山地地形作出分析。采用二重積分的集合意義，即二重積分是曲頂柱體的體積，利用這個思想，進(jìn)行合理假設(shè)，對曲頂柱體進(jìn)行拆分，分成多個細(xì)小的小區(qū)頂柱體，再對這些小曲頂柱體的體積求和，并考慮挖土，和填土的費(fèi)用代價(jià)，求出土石方量最小時(shí)的海拔高度。在不同的底面投影位置下，分別求出對應(yīng)的費(fèi)用消耗，求出最小費(fèi)用，從而決定從什么地方，什么海拔高度平整一塊固定大小的長方形連片土地，能使總的土石方量最小，使工程施工經(jīng)濟(jì)上達(dá)到最優(yōu)。關(guān)鍵詞： 合理開挖土地

6、最優(yōu) 土石方量最小一、問題提出及其分析1. 提出問題A市是一個山區(qū)城市，向山要地是A市發(fā)展的一個必然的選擇，但是如何在一片山地之中選擇合適的方位與開挖深度，從而使總的土石方量最小，就是一個十分有意義的課題.A市某工廠為了在一片長度為1500米，寬度為900米的山地之中，開挖出一個800米600米平坦連續(xù)的長方形地塊作為工廠的廠房地基，前期已經(jīng)在本塊土地上測量出長、寬每隔30米的網(wǎng)格的對應(yīng)網(wǎng)格點(diǎn)的海拔高度（詳細(xì)數(shù)據(jù)見附件）。請你考慮以下幾個問題：問題一：用附件中的數(shù)據(jù)畫出工廠的這片土地的三維圖形與等高線圖；問題二：從什么地方，什么海拔高度平整一塊800米600米的連片土地能使總的土石方量最?。?

7、. 問題分析1) 問題一：我們可以利用MATLAB軟件直接調(diào)用surf函數(shù)和contour函數(shù)分別 畫出這片土地的三維圖形和等高線圖。2) 問題二：考慮要平整一塊800米600米的連片土地，且要使土石方量最小， 我們需要確定這塊800米600米的連片土地的底面投影區(qū)域，以及海拔高度，使得其土石方量最小。 連片土地底面投影區(qū)域的確定：分析一：對于指定底面是800米600米大小區(qū)域內(nèi)，取海拔的高度做方差，如果方差最小，說明海拔都比較接近平均值，就是說明這片區(qū)域比較平坦，挖土或填土的土石方量近似最小。我們可以通過這個思路解決問題，即為求出滿足海拔方差值最小的底面位置的問題。分析二：讓底面投影區(qū)域在長

8、度為1500米，寬度為900米內(nèi)枚舉，此時(shí)計(jì)算出對應(yīng)的土石方量最小的是，在這些所有的最小值中取得最小值的那塊投影區(qū)域即為所求。我們采用此種方法求解。 連片土地海拔高度的確定：在平整土地的過程中，凸出的地方是要挖土的，凹下的地方是可以填土的，這樣挖下來的土可以填充到凹下去的部分。通過計(jì)算體積的方式計(jì)算挖土和填土大小，求出最小值對應(yīng)的海拔高度，即可確定平整連片土地的具體高度。二、問題的假設(shè)1、在平整土地的過程中，有些地方是要挖山的，但有些地方是要填土的，假設(shè)填土的每立方米所需的費(fèi)用為挖山的每立方米土石方所需費(fèi)用的1/3。2、假設(shè)除了挖土和填土以外，在平整土地的過程中其他作業(yè)（如登高）產(chǎn)生的費(fèi)用都與

9、800米600米的連片土地所處的位置方向和海拔高度均無關(guān)。這樣我們只需考慮挖土和填土的土石方量及費(fèi)用，據(jù)此來考慮連片土地所處的位置方向和海拔高度。3、假設(shè)計(jì)算體積的過程中，分割的小曲頂柱體不能達(dá)到無窮小，而產(chǎn)生的誤差忽略不計(jì)。實(shí)際計(jì)算中我們?nèi)∫粋€很小的步長去劃分，使其劃分盡可能的小，產(chǎn)生的誤差忽略不計(jì)。三、 符號說明符號說明L平整連片土地塊的長W平整連片土地塊的寬D平整連片土地塊的底面所在閉區(qū)域平整連片土地塊底面劃分的子區(qū)域子區(qū)域的面積平整土地所需的挖土量平整土地所需的填土量C平整土地挖土和填土所需要的總費(fèi)用四、模型的建立與求解4.1. 問題一求解1利用MATLAB軟件的三維繪圖功能，畫出工廠

10、這片土地的三維圖形如下：圖表4.1.1山地的三維圖形2，這片土地的等高線圖如下：圖表4.1.2山地的三維等高線圖圖表4.1.3山地的二維等高線圖4.2. 問題二模型建立與求解4.2.1. 平整塊的海拔高度的確定通過土石方量費(fèi)用最小的原則，確定平整土地海拔的開挖高度。平整土地快所對應(yīng)的的底面區(qū)域記為D，顯然D是一個800米600米的矩形區(qū)域。設(shè)頂部海拔高度是一個非負(fù)連續(xù)函數(shù)，。假設(shè) D的位置已經(jīng)確定，下面我們來討論在什么樣的海拔高度下，土石方量費(fèi)用最小?？紤]底面區(qū)域D對應(yīng)的山體是一個曲頂柱體，把區(qū)域D分成個小區(qū)域，為了簡化運(yùn)算，我們分成個等大的區(qū)域個等大的區(qū)域，每個區(qū)域都是邊長為的正方形。再分別

11、以這些小區(qū)域的邊界為準(zhǔn)線，作母線平行于軸的柱面，這些柱面將原曲頂柱體分成個細(xì)的曲頂柱體。每個區(qū)域上任取一點(diǎn)。當(dāng)很大時(shí)，所分的區(qū)域很小，邊長的值很小，每個區(qū)域?qū)?yīng)的海拔高度z可以近似為一個相同的值，記所有小區(qū)域的最小海拔高度為，最大海拔高度為設(shè)在平整塊連片土地建在海拔高度為處。當(dāng)小區(qū)域所對應(yīng)的海拔高度時(shí)，說明這塊區(qū)域需要挖土，挖土量為，這里表示的面積。時(shí)，說明這塊區(qū)域需要填土，填土量為，總的挖土量，其中為滿足條件的小區(qū)域的個數(shù)?？偟奶钔亮?， ，其中為滿足條件的小區(qū)域的個數(shù)。當(dāng)，挖土量大于填土量，這樣挖出來的土足夠填充凹下的山地?？偟耐潦搅抠M(fèi)用為：當(dāng)，挖土量小于填土量，這樣挖出來的土不足以填充凹

12、下的山地。需要從別的地方挖土來填充。總的土石方量費(fèi)用為：因此如果讓平整塊的底面投影區(qū)域取遍所有可能的值，再根據(jù)上述思想求解此底面投影區(qū)域?qū)?yīng)土石方量費(fèi)用最小的開挖海拔高度。按照最小的對應(yīng)的底部和開挖海拔高度來開挖，所得平整的連片土地能使總的土石方量最小。4.2.2. 平整塊的底面投影區(qū)域的確定將平整塊區(qū)域投影到底面（XOY平面），確定開挖方向，即要確定平整塊的底面位置。顯然底面是800米600米的矩形，假設(shè)矩形的四個頂點(diǎn)分別為A、B、C、D。山地是長度為1500米，寬度為900米的矩形區(qū)域，因此平整塊底面投影矩形ABCD可以在1500900的區(qū)域范圍內(nèi)任意移動，四個頂點(diǎn)不能超出這個范圍我們要取

13、遍所有可能的位置，必須枚舉所有的可能值。設(shè)矩形ABCD的邊長為AB=W=600，AC=L=800，且各點(diǎn)坐標(biāo)位置也已經(jīng)確定，設(shè)各頂點(diǎn)的坐標(biāo)分別為A（，）、B（，）、C（，）、D（，）。設(shè)AB邊的中點(diǎn)為E，CD邊的中點(diǎn)為F，那么唯一確定了線段EF。OxyADBCEF由中點(diǎn)公式，中點(diǎn)E，F(xiàn)的坐標(biāo)分別可以表示為E（，），F(xiàn)（，），且EF的邊長為EF=AC=800,因此只要矩形ABCD的位置和大小確定，那么EF的大小和位置就唯一確定。反過來如果EF的位置和大小確定了（E,F分別為短邊的中點(diǎn)），那么矩形ABCD的位置和大小也唯一確定。下面我們分別說明矩形四個頂點(diǎn)坐標(biāo)，以及四個邊的直線方程計(jì)算方法。設(shè)EF

14、=800，E、F兩點(diǎn)的坐標(biāo)分別為E（，）、F（，），為了使計(jì)算不重復(fù)，保證EF的唯一性，我們不妨設(shè) ，且當(dāng)時(shí), ,若EF所在直線的斜率存在，設(shè)為。當(dāng)時(shí)，斜率存在，分以下四種情況討論：(1)當(dāng),即EF平行于y軸時(shí)此時(shí)矩形如下圖所示：xyOADBCEF此時(shí)，四條邊對應(yīng)的直線方程分別為AB：，BC：，CD：，AD：，四個頂點(diǎn)的坐標(biāo)分別為：A（，）, B（，），C（，）, D（，）矩形ABCD的充分必要條件是：(2)當(dāng)時(shí)，xyOADBCEF此時(shí)，四條邊對應(yīng)的直線方程分別為AB: ,即： CD: ,即：AD：BC: 任意兩條直線的交點(diǎn)即為頂點(diǎn)，A、B、C、D四個頂點(diǎn)的坐標(biāo)分別由以下方程組解出：矩形ABC

15、D的充分必要條件是：(3)當(dāng)，即EF垂直于x軸時(shí)xyOADBCEF此時(shí)，四條邊對應(yīng)的直線方程分別為AB：，BC：，CD：，AD：，四個頂點(diǎn)的坐標(biāo)分別為：A（，）, B（，），C（，）, D（，）矩形ABCD的充分必要條件是：(4)當(dāng)，xyODABCEF此時(shí)，四條邊對應(yīng)的直線方程分別為AB: ,即： CD: ,即：BC：AD: 任意兩天直線的交點(diǎn)即為頂點(diǎn)，A、B、C、D四個頂點(diǎn)的坐標(biāo)分別由以下方程組解出：矩形ABCD的充分必要條件是：如果讓EF在平面內(nèi)取遍所有可能的值，再根據(jù)上述的思想求解可以確定平整塊的底面投影區(qū)域的位置，再根據(jù)4.2.1的介紹求解海拔位置。為了讓EF能取遍所有可能值，先在山地

16、海底平面上任意取一點(diǎn)作為E點(diǎn)，再根據(jù)E點(diǎn)確定F點(diǎn)的位置。這樣只需考慮E點(diǎn)取遍所有值的可能情況，根據(jù)EF的長度，確定F點(diǎn)可能的有效位置。平面區(qū)域內(nèi)任意取一點(diǎn)E， ，以E為圓心，R=L=800為半徑作圓，顯然F可能取得圓周上的所有點(diǎn)。OxyEFnF1F3F2為了避免重復(fù)，令 ，且當(dāng)時(shí), ,這樣F的取值為上圖所示的圓弧上。這樣線段EF就可以確定下來。4.2.3. 編程求解具體開挖位置根據(jù)4.2.1和4.2.2介紹的方法，編寫MATLAB程序求解。值得說明的是，在計(jì)算過程中，我們根據(jù)已有的檢測數(shù)據(jù)，用了插值的方法，多維插值擬合計(jì)算某一點(diǎn)的海拔高度，如線性插值，三次樣條插值等，計(jì)算思想在數(shù)值分析中有說明

17、，在這里不做細(xì)說，直接調(diào)用MATLAB中的多維插值函數(shù)interp2()求解即可。根據(jù)算法的基本思想，作出程序流程圖見附錄1，MATLAB程序見附錄2.程序的運(yùn)性結(jié)果為，A，B，C，C，海拔高度，即在以這樣的ABCD構(gòu)成的矩形區(qū)域?yàn)榈酌妫诤０蚊滋幤秸粔K800米600米的連片土地，所所需要的土石方量費(fèi)用最低。五、模型的評價(jià)與推廣模型最大的優(yōu)點(diǎn)在于對原始數(shù)據(jù)進(jìn)行擬合時(shí)，使用MATLAB軟件描述山地的三維圖形和等高線圖形，在此基礎(chǔ)上可以得到具有很高的擬合度和適度性。對模型的進(jìn)一步討論便可以得到一系列可靠而實(shí)用的信息，所得結(jié)論與客觀事實(shí)很好地吻合，從而進(jìn)一步說明模型是合理的。在解決合理開挖土地的問

18、題中，通過土石方量費(fèi)用最小的原則，確定平整土地海拔的開挖高度。模型的缺點(diǎn)在于在建立模型時(shí)，由于要提高精度，必須減小循環(huán)迭代步長，從而使循環(huán)次數(shù)增大，增大了程序的運(yùn)算量，運(yùn)行數(shù)據(jù)量比較大，會消耗大量的運(yùn)算時(shí)間，所以造成程序運(yùn)行過程中結(jié)果出來比較慢，算法效率比較低。本模型雖然討論的是以通過土石方量費(fèi)用最小的原則來確定平整土地海拔的開挖高度，但是可以推廣到各種合理開挖土地的問題中，用類似的方法進(jìn)行建模求解。解決該問題時(shí)，我們首先利用MATLAB軟件的三維繪圖功能，畫出工廠這片土地的三維圖形，再確定平整塊開挖的方向。本文所建立的模型可以推廣到其他各種模型。參考文獻(xiàn)1 陳紀(jì)修，於崇華，金路，數(shù)學(xué)分析 第

19、二版，M，北京：高等教育出版社，1999. 2 李慶揚(yáng)，王能超，易大義，數(shù)值分析 第五版，M，北京：清華大學(xué)出版社，2008. 3 韓明，張積林，李林，林杰等，數(shù)學(xué)建模案例，M，北京：同濟(jì)大學(xué)出版社，2012. 4 侯進(jìn)軍，肖艷清，譚敏，高明柯，數(shù)學(xué)建模方法與應(yīng)用，M，南京：東南大學(xué)出版社，2012. 5 宋世德，郭滿才，數(shù)學(xué)實(shí)驗(yàn)，M，北京：高等教育出版社，2002. 6 姜啟源等，數(shù)學(xué)模型，M，北京：高等教育出版社，2003.7 盛聚等，概率論與數(shù)理統(tǒng)計(jì)，M，北京：高等教育出版社，1983.8 薛定宇，陳陽泉，高等應(yīng)用數(shù)學(xué)問題的MATLAB求解，M，北京：清華大學(xué)出版社，2004。9 何正風(fēng)

20、，MATLAB在數(shù)學(xué)方面的應(yīng)用，M，北京：清華大學(xué)出版社，2012.附錄文件1 流程圖問題二解決方案的程序流程圖如下圖所示：2. 求解過程中的程序（已經(jīng)在MATLAB Version 7.10.0.499(R2010a)環(huán)境下調(diào)試通過）。%腳本文件drawpicture.m%用附件中的數(shù)據(jù)畫出工廠的這片土地的三維圖形與等高線x=0:30:900;y=0:30:1500;X,Y=meshgrid(x,y);Z=xlsread(data,1,b3:az33);surf(X,Y,Z); %三維圖xlabel(橫坐標(biāo)),ylabel(縱坐標(biāo)),zlabel(海拔高度);title(山地三維圖);fig

21、ure,contour(X,Y,Z); %二維等高線title(二維等高線圖);xlabel(橫坐標(biāo)),ylabel(縱坐標(biāo)),zlabel(海拔高度);figure,contour3(Z); %三維等高線title(三維等高線圖);xlabel(橫坐標(biāo)),ylabel(縱坐標(biāo)),zlabel(海拔高度);%函數(shù)文件getypos.m%功能：根據(jù)已知E點(diǎn)，及F點(diǎn)的x坐標(biāo)，取得滿足EF=800且yfye所對應(yīng)的F點(diǎn)的坐標(biāo)function yf=getypos(xe,ye,xf)%輸入:E點(diǎn)的坐標(biāo)(xe,ye),F點(diǎn)的x坐標(biāo)xf%輸出:F的y坐標(biāo)yfyf=sqrt(800*800-(xe-xf)2

22、) + ye;%函數(shù)文件isvalidpoint.m%功能：判斷P是否是有效的點(diǎn)function ret = isvalidpoint( p )% p是向量，長度為2，分別存放x,y軸的坐標(biāo)if (p(1) 0) & (p(1) 0) & (p(2) 0 A1 = 1 k;k -1; %AB和AD的直線方程系數(shù) b1 = pe(1)+k*pe(2),(-sqrt(1+k*k)*width/2-pe(2)+k*pe(1); %AB和AD的直線方程對應(yīng)的常數(shù)項(xiàng) pa = A1b1; A2 = 1 k; k (-1); %AB和BC的直線方程系數(shù) b2 = (pe(1)+k*pe(2),(sqrt(

23、1+k2)*width/2-pe(2)+k*pe(1); %AB和BC的直線方程對應(yīng)的常數(shù)項(xiàng) pb = A2b2; A3 = 1 k; k -1; %CD和BC的直線方程系數(shù) b3 = pf(1)+k*pf(2),sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和BC的直線方程對應(yīng)的常數(shù)項(xiàng) pc = A3b3; A4 = 1 k;k -1; %CD和AD的直線方程系數(shù) b4 = pf(1)+k*pf(2),(-sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和AD的直線方程對應(yīng)的常數(shù)項(xiàng) pd = A4b4; else if k 0 A1 =

24、1 k;k -1; %AB和BC的直線方程系數(shù) b1 = pe(1)+k*pe(2),(-sqrt(1+k2)*width/2-pe(2)+k*pe(1); %AB和BC的直線方程對應(yīng)的常數(shù)項(xiàng) pb = A1b1; A2 = 1 k; k (-1); %AB和AD的直線方程系數(shù) b2 = (pe(1)+k*pe(2),(sqrt(1+k2)*width/2-pe(2)+k*pe(1); %AB和AD的直線方程對應(yīng)的常數(shù)項(xiàng) pa = A2b2; A3 = 1 k; k -1; %CD和AD的直線方程系數(shù) b3 = pf(1)+k*pf(2),sqrt(1+k2)*width/2-pe(2)+k*

25、pe(1); %CD和AD的直線方程對應(yīng)的常數(shù)項(xiàng) pd = A3b3; A4 = 1 k;k -1; %CD和BC的直線方程系數(shù) b4 = pf(1)+k*pf(2),(-sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和BC的直線方程對應(yīng)的常數(shù)項(xiàng) pc = A4b4; end end endend A2 = 1 k; k (-1); %AB和AD的直線方程系數(shù) b2 = (pe(1)+k*pe(2),(sqrt(1+k2)*width/2-pe(2)+k*pe(1); %AB和AD的直線方程對應(yīng)的常數(shù)項(xiàng) pa = A2b2; A3 = 1 k; k -1; %CD和

26、AD的直線方程系數(shù) b3 = pf(1)+k*pf(2),sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和AD的直線方程對應(yīng)的常數(shù)項(xiàng) pd = A3b3; A4 = 1 k;k -1; %CD和BC的直線方程系數(shù) b4 = pf(1)+k*pf(2),(-sqrt(1+k2)*width/2-pe(2)+k*pe(1); %CD和BC的直線方程對應(yīng)的常數(shù)項(xiàng) pc = A4b4; end end endend%函數(shù)文件isinregion.m%功能：判斷點(diǎn)pt是否在矩形區(qū)域內(nèi)function ret = isinregion( pt , pe, pf)% pe, pf

27、分別表示E,F點(diǎn)的坐標(biāo)% pt表示待判斷的點(diǎn)，是否在矩形區(qū)域內(nèi)ret = 1;% 0表示不再區(qū)域內(nèi)，1表示不再區(qū)域內(nèi)long = 800;width = 600;x0 = pt(1);y0 = pt(2);if pe(2) = pf(2) %EF平行于x軸 if (x0 pf(1) | (y0 (pf(2) + width/2) ret = 0; endelse if x0 = pf(1) %EF平行于x軸 if (x0 (pf(1) - long/2) | (y0 pf(2) ) ret = 0; end else k = (pf(2) - pe(2)/(pf(1) - pf(1); if

28、k 0 ab = x0 + k*y0 - pe(1) - k * pe(2); cd = x0 + k*y0 - pf(1) - k * pf(2); ad = -k*x0 + y0 - sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); bc = -k*x0 + y0 + sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); if (ab 0) | (bc 0) | (cd 0) ret =0; end if (ab 0) | (bc 0) | (cd 0) ret =0; end else if k 0 ab = x0 + k*

29、y0 - pe(1) - k * pe(2); cd = x0 + k*y0 - pf(1) - k * pf(2); bc = -k*x0 + y0 - sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); ad = -k*x0 + y0 + sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); if (ab 0) | (bc 0) | (cd 0) ret =0; end if (ab 0) | (ad 0) | (cd 0) ret =0; end end end endendend ad = -k*x0 + y0 - sqrt

30、(1 + k2)*width/2 - pe(2) + k * pe(1); bc = -k*x0 + y0 + sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); if (ab 0) | (bc 0) | (cd 0) ret =0; end if (ab 0) | (bc 0) | (cd 0) ret =0; end else if k 0 ab = x0 + k*y0 - pe(1) - k * pe(2); cd = x0 + k*y0 - pf(1) - k * pf(2); bc = -k*x0 + y0 - sqrt(1 + k2)*width/2

31、 - pe(2) + k * pe(1); ad = -k*x0 + y0 + sqrt(1 + k2)*width/2 - pe(2) + k * pe(1); if (ab 0) | (bc 0) | (cd 0) ret =0; end if (ab 0) | (ad 0) | (cd 0) ret =0; end end end endendend%函數(shù)文件getmaxminpoint.m%功能：取得矩形區(qū)域內(nèi)的點(diǎn)對應(yīng)的橫縱坐標(biāo)的最大值和最小值function xmin, xmax, ymin, ymax = getmaxminpoint( pa, pb, pc, pd )xs=pa(

32、1), pb(1), pc(1), pd(1);ys=pa(2), pb(2), pc(2), pd(2);xmin = min(xs);xmax = max(xs);ymin = min(ys);ymax = max(ys);end%腳本文件main.m%計(jì)算并輸出最小石方量對應(yīng)的底面矩形坐標(biāo)和海拔高度x=0:30:900;y=0:30:1500;long = 800;%平整塊矩形區(qū)域的長width = 600;%平整塊矩形區(qū)域的寬X,Y=meshgrid(x,y);Z=xlsread(data,1,b3:az33);%disp(Z);dp=30;%步長dr=5;%分區(qū)域用dh=0.2;%海

33、拔高度步長direct=zeros(4,2);%開挖方向，存放底面投影的四個頂點(diǎn)seaheight=0;%開挖海拔高度minmincost=0;%最小費(fèi)用代價(jià)pexi=0;% E的橫坐標(biāo)while pexi 900 % 循環(huán)E的橫坐標(biāo) peyi=0;% E的縱坐標(biāo) fprintf(pexi); disp(pexi); while peyi 1500 % 循環(huán)E的縱坐標(biāo) % E(pexi, peyi) pfxj=max(pexi-long,0);%F的橫坐標(biāo)初始值 fprintf(peyi); disp(peyi); maxpfxj=min(pexi + long),900);%小于900范圍內(nèi)

34、取值 while (pfxj maxpfxj) % 循環(huán)F的橫坐標(biāo)，窮舉F的所有可能情況，由距離確定F的縱坐標(biāo) pfyj=getypos(pexi, peyi, pfxj);%F的縱向坐標(biāo) % F(pfxj, pfyj) pe = pexi peyi;%E點(diǎn)的坐標(biāo) ret=isvalidpoint(pe);%F判斷是否是有效范圍內(nèi)的點(diǎn) if ret = 0 pfxj = pfxj + dp; continue; end %取得ABCD四點(diǎn)坐標(biāo)，并判斷是否超出總的范圍 pf = pfxj, pfyj;%F點(diǎn)的坐標(biāo) pa, pb, pc, pd = getpoints( pe, pf );%取得A

35、BCD四點(diǎn)坐標(biāo) %至少有一個點(diǎn)超出范圍 if (isvalidpoint(pa) = 0) | (isvalidpoint(pb) = 0) | (isvalidpoint(pc) = 0) | (isvalidpoint(pd) = 0) pfxj = pfxj + dp; continue; end % 固定矩形區(qū)域內(nèi)求土石方量費(fèi)用最小 % mincost = 0;%最小代價(jià) height = 0;%最小代價(jià)，對應(yīng)的海拔高度 MP=zeros(ceil(long+1)*(width+1)/dp/dp), 3);%分小區(qū)域數(shù)long*width/d/d，記錄小區(qū)域的任意一點(diǎn) xmin, xm

36、ax, ymin, ymax = getmaxminpoint( pa, pb, pc, pd );%找到區(qū)域內(nèi)的最高點(diǎn)和最低點(diǎn)的橫縱坐標(biāo) px = xmin; index=1; while px xmax %fprintf(px); %disp(px); py = ymin; while py ymax %判斷點(diǎn)P(px,py)是否在矩形區(qū)域內(nèi) 在此寫一個函數(shù)去判斷.預(yù)留 pt = px, py; ret = isinregion( pt , pe, pf); if ret=0 py = py + dr; continue; end %如果該點(diǎn)在矩形區(qū)域內(nèi)，記錄該點(diǎn).為了節(jié)約時(shí)間，以空間換時(shí)

37、間。 MP(index, 1) = px; MP(index, 2) = py; MP(index, 3) = interp2(x, y, Z, px, py);%多維插值，取值 index = index + 1; py = py + dr; end px = px + dr; end %如果在區(qū)域內(nèi)在此區(qū)域內(nèi)求出體積最大值 mh = MP(:,3); minheight = min(mh); maxheight = max(mh); zh = minheight; while zh maxheight % 循環(huán)所有高度（從最低點(diǎn)的高度到最高點(diǎn)的高度） S = 0;% 記錄挖土量 T = 0

38、;% 記錄填土量 index=1; len = size(mh, 1); while (index 0 S = S + temp*dr*dr;% 挖土量 else T = T - temp*dr*dr;% 填土量 end index = index + 1; end if (S T) cost = 4*T/3; else cost = T/3 + S; end if (cost mincost) | (mincost = 0) mincost = cost; height = zh; end zh = zh+dh; end % 固定矩形區(qū)域 % if (minmincost mincost) | (minmincost = 0) minmincost = mincost; seaheight = height; direct = pa, pb, pc, pd; end pfxj = pfxj + dp; end peyi = peyi + dp; end pexi = pexi + dp;endfprintf(the final direct);disp(direct);fprintf(the final seaheight);disp(seaheight); %至少有一個點(diǎn)