3.2 非線(xiàn)性規(guī)劃L I NGO程序設(shè)計(jì)基礎(chǔ)

線(xiàn)性規(guī)劃的應(yīng)用范圍十分廣泛，但仍存在較大局限性，對(duì)許多實(shí)際問(wèn)題不能處理。非線(xiàn)性規(guī)劃比線(xiàn)性規(guī)劃有著更強(qiáng)的適用性。事實(shí)上，客觀世界中的問(wèn)題多是非線(xiàn)性的，給予線(xiàn)性處理大多是近似的，是在做了科學(xué)的假設(shè)和簡(jiǎn)化后得到的。在實(shí)際問(wèn)題中，有一些不能進(jìn)行線(xiàn)性化處理，否則將嚴(yán)重影響模型對(duì)實(shí)際問(wèn)題近似的可依賴(lài)性。但是，非線(xiàn)性規(guī)劃問(wèn)題在計(jì)算上通常比較困難，理論上的討論，也不能像線(xiàn)性規(guī)劃的問(wèn)題那樣，給出簡(jiǎn)捷的形式和透徹全面的結(jié)論。

例3.3 板材切割問(wèn)題

某鋼管零售商從鋼管廠進(jìn)貨，將鋼管按照顧客的要求切割后售出。從鋼管廠進(jìn)貨時(shí)得到的都是長(zhǎng)度為19m原料鋼管。現(xiàn)有一個(gè)客戶(hù)需要長(zhǎng)度為4m的鋼管50根、長(zhǎng)度為6m的鋼管20根和長(zhǎng)度為8m的鋼管15根。應(yīng)如何下料最節(jié)?。?/p>

零售商如果采用的不同切割模式太多，將會(huì)導(dǎo)致生產(chǎn)過(guò)程的復(fù)雜性，從而增加生產(chǎn)和管理成本。所以，該零售商規(guī)定采用的不同切割模式不能超過(guò)3種。此外，該客戶(hù)除需要上述中的三種鋼管外，還需要長(zhǎng)度為5m的鋼管10根，應(yīng)如何下料最節(jié)省？

模型設(shè)計(jì)

問(wèn)題一

首先，應(yīng)當(dāng)確定哪些切割模式是可行的，所謂一個(gè)切割模式，是指按照客戶(hù)需要在原料鋼管上安排切割的一種組合。例如，我們可以將長(zhǎng)度為19m的鋼管切割成3根長(zhǎng)度為4m的鋼管，余料為7m；或者將長(zhǎng)度為19m的鋼管切割成長(zhǎng)度為4m、6m和8m的鋼管各1根，余料為1m。顯然，可行的切割模式可以有很多。

其次，應(yīng)當(dāng)確定哪些切割模式是合理的。通常，假設(shè)一個(gè)合理的切割模式的余料不應(yīng)該大于或等于客戶(hù)需要的鋼管的最小尺寸。例如，將19m長(zhǎng)的鋼管切割成3根4m的鋼管是可行的，但余料為7m，可以進(jìn)一步將7m的余料切割成4m鋼管（余料為3m），或者將7m的余料切割成6m鋼管（余料為1m）。在這種合理性假設(shè)下，切割模式一共有7種，如表3-6所示。

問(wèn)題化為在滿(mǎn)足客戶(hù)需要的條件下，按照哪些合理的模式，切割多少根原料鋼管最為節(jié)省。所謂節(jié)省可以有兩種標(biāo)準(zhǔn)：一是切割后剩余的總余料最小，二是切割原料鋼管的總根數(shù)最少。下面將對(duì)這兩個(gè)目標(biāo)分別討論。

決策變量：用xi表示按照第i種模式（i=1,2, …,7）切割的原料鋼管的根數(shù)，顯然它們應(yīng)當(dāng)是非負(fù)數(shù)。決策目標(biāo)：以切割后剩余的總余量最小為目標(biāo)，由表3-6可得到：

minZ1=3x1+x2+3x3+3x4+x5+x6+3x7

以切割原料鋼管的總根數(shù)最少為目標(biāo)，則有：

minZ2=x1+x2+x3+x4+x5+x6+x7

約束條件為滿(mǎn)足客戶(hù)的要求，則有：

【思考題】 請(qǐng)思考兩種目標(biāo)函數(shù)下，取得的最優(yōu)方案是否相同？為什么？如果兩種最優(yōu)方案不同，請(qǐng)思考何種更為合理？

問(wèn)題二

按照解問(wèn)題一的思路，可以通過(guò)枚舉法首先確定哪些切割模式是可行的，但由于需求的鋼管規(guī)格增加到4種，所以枚舉法的工作量較大。

通過(guò)枚舉法可以確定16種可行且合理的切割模式。決策變量：用xi表示按照第i種模式（i=1,2, …,16）切割的原料鋼管的根數(shù)，顯然它們應(yīng)當(dāng)是非負(fù)數(shù)。決策目標(biāo)：以切割原料鋼管的總根數(shù)最少為目標(biāo)，則有：

分析第二問(wèn)的難點(diǎn)在于：該零售商規(guī)定采用的不同切割模式不能超過(guò)3種。為了實(shí)現(xiàn)該要求，引入額外的決策變量：用ri 表示采用第i種模式（i=1,2, …,16）。

因此，對(duì)切割原料鋼管的總根數(shù)最少的目標(biāo)需要修正如下：

以上目標(biāo)函數(shù)可以理解為：所采用的切割模式下所用的原料鋼管總數(shù)。

約束條件為滿(mǎn)足客戶(hù)的要求，則有：

其中，dj表示用戶(hù)對(duì)于第j種類(lèi)型鋼管的需求數(shù)，j=1,2,3,4分別表示長(zhǎng)度為4m、5m、6m、8m的鋼管。yij表示表3-7中第i種模式、第j種類(lèi)型鋼管的根數(shù)。

下面介紹一種可帶有普遍性、可以同時(shí)確定切割模式和切割計(jì)劃的方法。同問(wèn)題（1）類(lèi)似，一個(gè)合理的切割模式的余料不應(yīng)該大于或等于客戶(hù)需要的鋼管的最小尺寸（本題為4m）。切割計(jì)劃中只使用合理的切割模式，而由于本題中參數(shù)都是整數(shù)，所以合理的切割模式的余量不能大于3m。此外，這里僅選擇總根數(shù)最少為目標(biāo)進(jìn)行求解。

決策變量 由于不同切割模式不能超過(guò)3種，可以用xi表示按照第i 種模式（i=1,2, 3）切割的原料鋼管的根數(shù)，顯然它們應(yīng)當(dāng)是非負(fù)整數(shù)。設(shè)所使用的第i種切割模式下每根原料鋼管生產(chǎn)長(zhǎng)度為4m、5 m、6 m和8m的鋼管數(shù)量分別為r1i, r2i, r3i, r4i（非負(fù)整數(shù)）。

決策目標(biāo) 以切割原料鋼管的總根數(shù)最少為目標(biāo)，即目標(biāo)為：

min x1+x2+x3

約束條件 為滿(mǎn)足客戶(hù)的需求，應(yīng)有

每一種切割模式必須可行、合理，所以每根原料鋼管的成品量不能超過(guò)19m，也不能小于16m（余料不能大于3m），于是：

與LINDO相比，LINGO軟件主要具有兩大優(yōu)點(diǎn)：除具有LINDO的全部功能外，還可用于求解非線(xiàn)性規(guī)劃問(wèn)題，包括非線(xiàn)性整數(shù)規(guī)劃問(wèn)題；內(nèi)置建模語(yǔ)言，允許以簡(jiǎn)練、直觀的方式描述較大規(guī)模的優(yōu)化問(wèn)題，所需的數(shù)據(jù)可以以一定格式保存在獨(dú)立的文件中。

LINGO不詢(xún)問(wèn)是否進(jìn)行敏感性分析，敏感性分析需要將來(lái)通過(guò)修改系統(tǒng)選項(xiàng)啟動(dòng)敏感性分析后，再調(diào)用“REPORT|RANGE”菜單命令來(lái)實(shí)現(xiàn)。現(xiàn)在同樣可以把模型和結(jié)果報(bào)告保存在文件中。

一般來(lái)說(shuō)，LINGO中建立的優(yōu)化模型可以由五個(gè)部分組成，或稱(chēng)為五“段”（SECTION）：

1．集合段（SETS）：以“SETS”開(kāi)始，“ENDSETS”結(jié)束，定義必要的集合變量（SET）及其元素（MEMBER，含義類(lèi)似于數(shù)組的下標(biāo)）和屬性（ATTRIBUTE，含義類(lèi)似于數(shù)組）。

2．目標(biāo)與約束段：目標(biāo)函數(shù)、約束條件等，沒(méi)有“段”的開(kāi)始和結(jié)束標(biāo)記，因此實(shí)際上就是除其他四個(gè)段（都有明確的段標(biāo)記）外的LINGO模型。這里一般要用到LINGO的內(nèi)部函數(shù)，尤其是與集合相關(guān)的求和函數(shù)@SUM和循環(huán)函數(shù)@FOR等。

3．?dāng)?shù)據(jù)段（DATA）：以“DATA”開(kāi)始，“ENDDATA”結(jié)束，對(duì)集合的屬性（數(shù)組）輸入必要的常數(shù)數(shù)據(jù)。格式為：“attribute（屬性）=value_list（常數(shù)列表）; ”常數(shù)列表（value_list）中數(shù)據(jù)之間可以用逗號(hào)“, ”分開(kāi)，也可以用空格分開(kāi)（回車(chē)等價(jià)于一個(gè)空格）。

4．初始段（INIT）：以“INIT”開(kāi)始，“ENDINIT”結(jié)束，對(duì)集合的屬性（數(shù)組）定義初值（因?yàn)榍蠼馑惴ㄒ话闶堑惴?，所以用?hù)如果能給出一個(gè)比較好的迭代初值，對(duì)提高算法的計(jì)算效果是有益的）。如果有一個(gè)接近最優(yōu)解的初值，對(duì)LINGO求解模型是有幫助的。定義初值的格式為：“attribute（屬性）=value_list（常數(shù)列表）; ”

5．計(jì)算段（CALC）：以“CALC”開(kāi)始，“ENDCALC”結(jié)束，對(duì)一些原始數(shù)據(jù)進(jìn)行計(jì)算處理。在實(shí)際問(wèn)題中，輸入的數(shù)據(jù)通常是原始數(shù)據(jù)，不一定能在模型中直接使用，可以在這個(gè)段對(duì)這些原始數(shù)據(jù)進(jìn)行一定的“預(yù)處理”，得到模型中真正需要的數(shù)據(jù)。

普通LINGO程序設(shè)計(jì)：

    min=x1+x2+x3；
    x1*r1 1+x2*r1 2+x3*r1 3＞=5 0；
    x1*r2 1+x2*r22+x3*r23＞=1 0；
    x1*r3 1+x2*r32+x3*r33＞=2 0；
    x1*r41+x2*r42+x3*r43＞=15；
    4*r11+5*r21+6*r31+8*r41＜=19；
    4*r12+5*r22+6*r32+8*r42＜=19；
    4*r13+5*r23+6*r33+8*r43＜=19；
    4*r1 1+5*r2 1+6*r3 1+8*r4 1＞=1 6；
    4*r12+5*r22+6*r32+8*r42＞=16；
    4*r13+5*r23+6*r33+8*r43＞=16；
    x1+x2+x3＞=26；
    x1+x2+x3＜=3 1；
    @gin（x1）; @gin（x2）; @gin（x3）; @gin（r11）; @gin（r12）; @gin（r13）; @gin（r21）; @gin
（r22）; @gin（r23）; @gin（r31）; @gin（r32）; @gin（r33）; @gin（r41）; @gin（r42）; @gin
（r43）；
    End

運(yùn)行結(jié)果一：

    Local optimal solution found.
    Objective value:                        28.00000
    Extended solver steps:                       29 1
    Total solver iterations:                     9681

建?；Z(yǔ)言程序設(shè)計(jì)：

    Model：
     sets：
    number/1..3/:x；
    modes/1..4/:a, b；
    mode（number, modes）:r；

    endsets
    data：
    a=50,10,20,15；
    b=4,5,6,8；
    enddata
    min=@sum（number∶x）；
    @for（modes（i）:@sum（number（j）:x（j）*r（j, i））＞=a（i））；
    @for（number（i）:@sum（modes（j）:b（j）*r（i, j））＜=19）；
    @for（number（i）:@sum（modes（j）:b（j）*r（i, j））＞=1 6）；
    @sum（number∶x）＞=26；
    @sum（number∶x）＜=3 1；
    @for（number:@gin（x））；
    @for（mode:@gin（r））；
    End

運(yùn)行結(jié)果二：

    Local optimal solution found.
    Objective value:                        28.00000
    Objective bound:                        28.00000
    Infeasibilities:                        0.000000
    Extended solver steps:                   166
    Total solver iterations:                  11142

例3.4 太陽(yáng)影子定位問(wèn)題

如何確定視頻的拍攝地點(diǎn)和拍攝日期是視頻數(shù)據(jù)分析的重要方面，太陽(yáng)影子定位技術(shù)就是通過(guò)分析視頻中物體的太陽(yáng)影子變化，確定視頻拍攝的地點(diǎn)和日期的一種方法。

1．建立影子長(zhǎng)度變化的數(shù)學(xué)模型，分析影子長(zhǎng)度關(guān)于各個(gè)參數(shù)的變化規(guī)律，并應(yīng)用你們建立的模型畫(huà)出2015年10月22日北京時(shí)間9:00—15:00之間天安門(mén)廣場(chǎng)（北緯39°54′26″，東經(jīng)116°23′29″）3米高的直桿的太陽(yáng)影子長(zhǎng)度的變化曲線(xiàn)。

2．根據(jù)某固定直桿在水平地面上的太陽(yáng)影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，建立數(shù)學(xué)模型確定直桿所處的地點(diǎn)。將你們的模型應(yīng)用于表3-8的影子頂點(diǎn)坐標(biāo)數(shù)據(jù)，給出若干個(gè)可能的地點(diǎn)。

問(wèn)題分析

問(wèn)題一：本小題要求建立影子長(zhǎng)度變化模型，并將該模型應(yīng)用于實(shí)際地理背景，進(jìn)行太陽(yáng)影子的仿真計(jì)算。考慮引起影子長(zhǎng)度變化的直接因素是太陽(yáng)高度角和桿子測(cè)量高度，據(jù)此分析影響太陽(yáng)高度角變化的因子，最終得到影響影長(zhǎng)變化的因子有太陽(yáng)高度角、太陽(yáng)赤緯角、太陽(yáng)方位角和時(shí)角等參數(shù)。將模型應(yīng)用到題目所給時(shí)間和位置上進(jìn)行求解，繪制不同時(shí)刻下直桿太陽(yáng)影子長(zhǎng)度的變化曲線(xiàn)。

問(wèn)題二：本小題要求根據(jù)已知某日期和時(shí)刻所對(duì)應(yīng)的位置坐標(biāo)，建立模型并應(yīng)用到實(shí)際數(shù)據(jù)中，找出滿(mǎn)足要求的地理位置。考慮所給數(shù)據(jù)中坐標(biāo)系朝向未知，分析造成坐標(biāo)變化的因素，考慮軌跡隨時(shí)間的變化，基于坐標(biāo)變化思想，對(duì)任意時(shí)刻的方位角求解。用最小二乘法使真實(shí)軌跡和計(jì)算軌跡偏差最小。

模型設(shè)計(jì)

問(wèn)題一

考慮任意物體影子長(zhǎng)度都是由太陽(yáng)光照射引起的，且物體高度和太陽(yáng)高度角是影響影長(zhǎng)的兩個(gè)直接因素。某一地點(diǎn)太陽(yáng)高度角在不同時(shí)刻是變化的，由于地球存在公轉(zhuǎn)，且地球地心一直在赤道平面上，故太陽(yáng)直射點(diǎn)在南北回歸線(xiàn)之間移動(dòng)，進(jìn)而要考慮物體所在地理緯度?？紤]地球繞太陽(yáng)公轉(zhuǎn)時(shí)，地軸指向不變，導(dǎo)致太陽(yáng)和地心連線(xiàn)與赤道平面夾角發(fā)生變化，故需考慮太陽(yáng)赤緯角和太陽(yáng)方位角?？紤]地球自身在自轉(zhuǎn)，故在探討影子長(zhǎng)度幾何關(guān)系時(shí)還需要考慮時(shí)角等因素。

綜上分析，在研究影子長(zhǎng)度相關(guān)的幾何參數(shù)時(shí)，考慮對(duì)太陽(yáng)方位角、太陽(yáng)赤緯角、太陽(yáng)高度角和時(shí)角等參數(shù)進(jìn)行分析定位，各參數(shù)的幾何意義分別為：

太陽(yáng)赤緯角δ：太陽(yáng)中心和地球中心的連線(xiàn)與地球赤道平面的夾角。其中，在春秋分時(shí)刻夾角最小為0°，在夏至和冬至?xí)r角度達(dá)到最大為±23°26′36″；

太陽(yáng)高度角h：地球表面上任意一點(diǎn)和太陽(yáng)的連線(xiàn)與地平線(xiàn)的夾角；

時(shí)角t：?jiǎn)挝粫r(shí)間內(nèi)地球自轉(zhuǎn)的角度，定義正午時(shí)角為0°，記上午時(shí)角為負(fù)值，下午時(shí)角為正值。

太陽(yáng)赤緯角是地球赤道面與日地中心連線(xiàn)的夾角，黃赤交角為黃道面與赤道面的交角，是赤緯角的最大值。查閱相關(guān)文獻(xiàn)，目前地球的黃赤交角約為23°26′，在日地二體系統(tǒng)中，選定地球?yàn)閰⒄瘴?，假設(shè)太陽(yáng)沿赤道面繞地球做勻速圓周運(yùn)動(dòng)。定義積日零點(diǎn)為2015年1月1日，以2015年春分日為3月21日，與2015年1月1日相差79天。由于春分時(shí)δ=0°，因此以春分日為基準(zhǔn)，則得到赤緯角的計(jì)算式為：

考慮太陽(yáng)的運(yùn)行軌跡，計(jì)算太陽(yáng)對(duì)地球上某一物體的影子長(zhǎng)度變化軌跡，查閱文獻(xiàn)得到主要由當(dāng)?shù)氐牡乩砭暥群蜁r(shí)間這兩個(gè)因素決定。假設(shè)物體所在地理緯度記為φ，地理經(jīng)度記為α，北京時(shí)間記為t0，北京時(shí)間所在經(jīng)度記為α0，各參數(shù)的意義如下：

時(shí)角t：考慮任一經(jīng)度位置的物體在計(jì)算時(shí)角時(shí)都是以北京時(shí)間為基準(zhǔn)，且規(guī)定正午時(shí)角為0°。故得到時(shí)角計(jì)算式為：

太陽(yáng)高度角h：入射至地表某點(diǎn)的太陽(yáng)光線(xiàn)與該點(diǎn)切平面夾角，計(jì)算公式為：

h=sin-1（sinφsinδ+cosφcosδcost）, h∈[-90°,90°]

太陽(yáng)赤緯角δ：太陽(yáng)中心和地球中心的連線(xiàn)與地球赤道平面的夾角，由赤緯角的理論推導(dǎo)得到公式為：

綜上分析得到，影子長(zhǎng)度函數(shù)式為：

l=H×coth=H×cot（sin-1（sinφsinδ+cosφcosδcost））

考慮在影長(zhǎng)模型中根據(jù)經(jīng)緯度和時(shí)間對(duì)各個(gè)時(shí)刻的太陽(yáng)高度角求解，進(jìn)而求得影子長(zhǎng)度。現(xiàn)將該模型應(yīng)用到實(shí)際地理背景中，已知拍攝地點(diǎn)為天安門(mén)廣場(chǎng)（北緯39°54′26″，東經(jīng)116°23′29″）時(shí)，且拍攝時(shí)間為2015年10月22日北京時(shí)間9:00—15:00。此時(shí)得到桿長(zhǎng)為3m的物體在各整點(diǎn)時(shí)的影子長(zhǎng)度如表3-9所示。

分析表3-9中數(shù)據(jù)可得，在該時(shí)段內(nèi)呈現(xiàn)的影子長(zhǎng)度先減小后增大，且總大于桿子實(shí)際長(zhǎng)度。繪制影子長(zhǎng)度隨時(shí)間連續(xù)變化曲線(xiàn)如圖3-6所示。

分析圖3-6中影子長(zhǎng)度變化軌跡得到，影子長(zhǎng)度在北京時(shí)間12:14時(shí)刻影子長(zhǎng)度最短為3.7992m，即此時(shí)刻太陽(yáng)高度角最大為38.296°。影子長(zhǎng)度變化率先減小后增大，即越靠近北京時(shí)間12:14時(shí)刻，太陽(yáng)高度角和影子長(zhǎng)度單位時(shí)間的變化率越小，且曲線(xiàn)是連續(xù)變化的。聯(lián)系實(shí)際背景可得，在北京時(shí)間10月22日，太陽(yáng)直射點(diǎn)在南半球，且往南回歸線(xiàn)方向移動(dòng)，故相對(duì)應(yīng)天安門(mén)這一地點(diǎn)的影子長(zhǎng)度均大于實(shí)際長(zhǎng)度是合理的，一定程度上也說(shuō)明了模型的可靠性。

問(wèn)題二

考慮對(duì)給定太陽(yáng)頂點(diǎn)影子坐標(biāo)數(shù)據(jù)，建立模型確定所觀測(cè)點(diǎn)的地理位置。建立確定直桿所處位置與時(shí)間的數(shù)學(xué)模型，由于描述太陽(yáng)狀態(tài)的角度對(duì)時(shí)間十分敏感，無(wú)法直接使用北京時(shí)間，為了計(jì)算直桿所處地點(diǎn)，需要轉(zhuǎn)換為真太陽(yáng)時(shí)，即當(dāng)?shù)氐牡胤綍r(shí)間，建立優(yōu)化模型，運(yùn)用最小二乘法使得預(yù)測(cè)影長(zhǎng)軌跡與實(shí)際影長(zhǎng)軌跡之間誤差最小，進(jìn)而得到可能的觀測(cè)位置。

由問(wèn)題一模型幾何求解得到影長(zhǎng)與太陽(yáng)高度角之間函數(shù)關(guān)系為：

l=H×coth

其中，h為太陽(yáng)高度角，h=arcsin（sinφsinδ+cosφcosδcost）。

考慮建立平面直角坐標(biāo)系xOy, x軸正方向?yàn)檎龞|，y軸正方向?yàn)檎?，則影長(zhǎng)在x軸和y 軸的分量分別為：

其中，A為太陽(yáng)方位角。

考慮表3-8中所給出的x和y 坐標(biāo)數(shù)據(jù)沒(méi)有規(guī)定坐標(biāo)軸方向，由于任意兩個(gè)時(shí)刻的坐標(biāo)系角度是變化的，故需要統(tǒng)一坐標(biāo)系。再用最小二乘法計(jì)算求解，將得到的影長(zhǎng)坐標(biāo)數(shù)據(jù)與真實(shí)坐標(biāo)數(shù)據(jù)進(jìn)行比較。坐標(biāo)變換示意圖如圖3-7所示。

將各時(shí)刻坐標(biāo)系旋轉(zhuǎn)β角后，由圖3-7轉(zhuǎn)換示意圖分析得到新舊坐標(biāo)變換公式為：

其中，β為兩平面直角坐標(biāo)系變化前后旋轉(zhuǎn)角度；x和y 為變化前坐標(biāo)；x′和y′為變化后坐標(biāo)。

求解的目標(biāo)是得到若干個(gè)可能的直桿所在地點(diǎn)，為使得陰影的長(zhǎng)度和位置即軌跡與附件所給盡可能相近，利用最小二乘法建立如下優(yōu)化模型。目標(biāo)函數(shù)為：

其中，xi和yi為真實(shí)的影子橫縱坐標(biāo)；x′i和y′i為預(yù)測(cè)得到的影子橫縱坐標(biāo)。

約束條件為：

參數(shù)關(guān)系為：

對(duì)上述幾何關(guān)系式中未知參數(shù)在全范圍遍歷使得計(jì)算得到軌跡與真實(shí)軌跡之間相對(duì)誤差盡可能小，并運(yùn)用最小二乘法，通過(guò)MATLAB對(duì)優(yōu)化模型進(jìn)行計(jì)算求解。通過(guò)對(duì)優(yōu)化模型求解，得到滿(mǎn)足影子軌跡最優(yōu)和影子長(zhǎng)度最優(yōu)目標(biāo)下的位置信息如表3-10所示?？紤]優(yōu)化模型求解后，得到通過(guò)影長(zhǎng)優(yōu)化的地理位置位于海南省。