Lingo案例分析-供參考

1、最小費(fèi)用運(yùn)輸問(wèn)題單位 銷(xiāo)地運(yùn)價(jià)產(chǎn)地B1B2B3B4B5B6B7B8產(chǎn)量A16267425960A24953858255A35219743351A47673927143A52395726541A65522814352銷(xiāo)量3537223241324338280/302i.e:model:!6發(fā)點(diǎn)8收點(diǎn)運(yùn)輸問(wèn)題;sets: warehouses/wh1.wh6/: capacity; vendors/v1.v8/: demand; links(warehouses,vendors): cost, volume;endsets!目標(biāo)函數(shù); min=sum(links: cost*volume);!需求

2、約束; for(vendors(J): sum(warehouses(I): volume(I,J)=demand(J);!產(chǎn)量約束; for(warehouses(I): sum(vendors(J): volume(I,J)=0, FLOOR returns the largest integer, I, such that I =X.LGM( X)This returns the natural (base e) logarithm of the gamma function of X (i.e., log of (X - 1)!). It is extended to noninteg

3、er values of X by linear interpolation.Note:伽瑪方程表達(dá)式為：(x)=e(-t)*t(x-1)dt (積分的下限式0,上限式+) 利用分部積分法(integration by parts)我們可以得到 (x)=（x-1)*(x-1) LOG( X)This returns the natural logarithm of X.MOD( X,Y)This returns the value of X modulo Y, or, in other words, the remainder of an integer divide of X by Y.PO

4、W( X,Y) This returns the value of X rasied to the Y power.SIGN( X)This returns -1 if X =Staff required today, for each day of the weekThe right-hand side of this expression, Staff required today, is easy to calculate. It is simply the quantity REQUIRED( I). The left-hand side, Staff on duty today, i

5、s a bit trickier to compute. Given that all employees are on a five day on, two day off schedule, the number of employees working today is:Number working today = Number starting today + Number starting 1 day ago + Number starting 2 days ago + Number starting 3 days ago + Number starting 4 days ago.I

6、n other words, to compute the number of employees working today, we sum up the number of people starting today plus those starting over the previous four days. The number of employees starting five and six days back dont count because they are on their days off. So, using mathematical notation, what

7、 one might consider doing is adding the constraint:= j-4, j STARTi ?REQUIRED j , for j?DAYS Translating into LINGO notation, we can write this as:FOR( DAYS( J): SUM( DAYS( I) | I #LE# 5: START( J - I + 1) = REQUIRED( J);In words, the LINGO statement says, for each day of the week, the sum of the emp

8、loyees starting over the five day period beginning four days ago and ending today must be greater-than-or-equal-to the required number of staff for the day. This sounds correct, but there is a slight problem. If we try to solve our model with this constraint we get the error message:To see why we ge

9、t this error message, consider what happens on Thursday. Thursday has an index of 4 in our set DAYS. As written, the staffing constraint for Thursday will be:START( 4 - 1 + 1) + START( 4 - 2 + 1) + START( 4 - 3 + 1) + START( 4 - 4 + 1) + START( 4 - 5 + 1) = REQUIRED( 4);Simplifying, we get:START( 4)

10、 + START( 3) + START( 2) + START( 1) + START( 0) = REQUIRED( 4);The START( 0) term is the root of our problem. START is defined for days 1 through 7. START( 0) does not exist. An index of 0 on START is considered out of range. We would like to have any indices less-than-or-equal-to 0 wrap around to

11、the end of the week. Specifically, 0 would correspond to Sunday (7), -1 to Saturday (6), and so on. LINGO has a function that does just this called WRAP. The WRAP function takes two argumentscall them INDEX and LIMIT. Formally speaking, WRAP returns J such that J = INDEX - K * LIMIT, where K is an i

12、nteger such that J is in the interval 1, LIMIT. Informally speaking, WRAP will subtract or add LIMIT to INDEX until it falls in the range 1 to LIMIT. Therefore, this is just what we need to wrap around an index in multiperiod planning models. Incorporating the WRAP function, we get the corrected, fi

13、nal version of our staffing constraint:FOR( DAYS( J): SUM( DAYS( I) | I #LE# 5: START( WRAP( J - I + 1, 7) = REQUIRED( J);The Solution：Below is our staffing model in its entirety:SETS: DAYS / MON TUE WED THU FRI SAT SUN/: REQUIRED, START;ENDSETSDATA: REQUIRED = 20 16 13 16 19 14 12;ENDDATAMIN = SUM(

14、 DAYS( I): START( I);FOR( DAYS( J): SUM( DAYS( I) | I #LE# 5: START( WRAP( J - I + 1, 7) = REQUIRED( J);Model: STAFFDEMSolving the model, we get the solution report:Optimal solution found at step: 8Objective value: 22.00000 Variable Value Reduced CostSTART( MON) 8.000000 0.0000000START( TUE) 2.00000

15、0 0.0000000START( WED) 0.0000000 0.0000000START( THU) 6.000000 0.0000000START( FRI) 3.000000 0.0000000START( SAT) 3.000000 0.0000000START( SUN) 0.0000000 0.0000000 Row Slack or Surplus Dual Price 1 22.00000 1.000000 2 0.0000000 -0.2000000 3 0.0000000 -0.2000000 4 0.0000000 -0.2000000 5 0.0000000 -0.

16、2000000 6 0.0000000 -0.2000000 7 0.0000000 -0.2000000 8 0.0000000 -0.2000000Solution to STAFFDEMThe objective value of 22 means we need to hire 22 workers. We start our workers according to the schedule:MonTueWedThuFriSatSunStart8206330If we look at the surpluses on our staffing requirement rows (ro

17、ws 2 - 7), we see that the slack values are 0 on all of the days. This means there are no more workers than required and we just meet staffing requirements on every day. Even though this is a small model, trying to come up with a solution this efficient by hand would be a difficult task.職員時(shí)序安排模型 一項(xiàng)工

18、作一周7天都需要有人（比如護(hù)士工作），每天（周一至周日）所需的最少職員數(shù)為20、16、13、16、19、14和12，并要求每個(gè)職員一周連續(xù)工作5天，試求每周所需最少職員數(shù)，并給出安排。注意這里我們考慮穩(wěn)定后的情況。model:sets: days/mon.sun/: required,start;endsetsdata: !每天所需的最少職員數(shù); required = 20 16 13 16 19 14 12; enddata!最小化每周所需職員數(shù); min=sum(days: start); for(days(J): sum(days(I) | I #le# 5:start(wrap(J-I

19、+1,7) = required(J);end計(jì)算的部分結(jié)果為Global optimal solution found at iteration: 0 Objective value: 22.00000 Variable Value Reduced Cost REQUIRED( MON) 20.00000 0.000000 REQUIRED( TUE) 16.00000 0.000000 REQUIRED( WED) 13.00000 0.000000 REQUIRED( THU) 16.00000 0.000000 REQUIRED( FRI) 19.00000 0.000000 REQU

20、IRED( SAT) 14.00000 0.000000 REQUIRED( SUN) 12.00000 0.000000 START( MON) 8.000000 0.000000 START( TUE) 2.000000 0.000000 START( WED) 0.000000 0.3333333 START( THU) 6.000000 0.000000 START( FRI) 3.000000 0.000000 START( SAT) 3.000000 0.000000 START( SUN) 0.000000 0.000000每周最少:22個(gè)職員，周一:8人，周二:2人，周三:0人

21、，周四:6人，周五和周六:3人，周日:0人。用Matlab編程為：先建立一個(gè)M文件：fssp.m然后再建立：%fssp.mf=1,1,1,1,1,1,1; %A=-1,0,0,-1,-1,-1,-1; -1,-1,0,0,-1,-1,-1; -1,-1,-1,0,0,-1,-1; -1,-1,-1,-1,0,0,-1; -1,-1,-1,-1,-1,0,0; 0,-1,-1,-1,-1,-1,0; 0,0,-1,-1,-1,-1,-1;b=-20,-16,-13,-16,-19,-14,-12;l=0,0,0,0,0,0,0;xo,fo,exitflag=linprog(f,A,b,l)結(jié)果為

22、：xo = 8.0000 2.0000 0.0000 6.0000 3.0000 3.0000 0.0000fo = 22.0000exitflag = 13.7集循環(huán)函數(shù)（Set Looping Functions）Set looping functions operate over an entire set and, with the exception of the FOR function, produce a single result. The syntax for a set looping function is:function( setname ( set_index_l

23、ist) | conditional_qualifier : expression_list);function corresponds to one of the set looping functions listed below. setname is the name of the set you want to loop over. The set_index_list is optional and is used to create a list of indices, which correspond to the parent primitive sets that form

24、 the setname set. As LINGO loops through the members of the setname set, it will set the values of the indices in the set_index_list to correspond to the current member of the setname set. The conditional_qualifier is optional and may be used to limit the scope of the set looping function. When LING

25、O is looping over each member of the setname set, it evaluates the conditional_qualifier. If the conditional_qualifier evaluates to true, then function is performed for the set member. Otherwise, it is skipped. The expression_list is a list of expressions to be applied to each member of the setname

26、set. When using the FOR function, the expression_list may contain multiple expressions, separated by semicolons. These expressions will be added as constraints to the model. When using the remaining three set looping functions (SUM, MAX, and MIN), the expression_list must contain one expression only

27、. If the set_index_list is omitted, all attributes referenced in the expression_list must be defined on the setname set.The available set looping functions are listed below:集循環(huán)函數(shù)遍歷整個(gè)集進(jìn)行操作。其語(yǔ)法為function(setname(set_index_list)|conditional_qualifier:expression_list);function相應(yīng)于下面羅列的四個(gè)集循環(huán)函數(shù)之一；setname是要遍

28、歷的集；set_ index_list是集索引列表；conditional_qualifier是用來(lái)限制集循環(huán)函數(shù)的范圍，當(dāng)集循環(huán)函數(shù)遍歷集的每個(gè)成員時(shí)，LINGO都要對(duì)conditional_qualifier進(jìn)行評(píng)價(jià)，若結(jié)果為真，則對(duì)該成員執(zhí)行function操作，否則跳過(guò)，繼續(xù)執(zhí)行下一次循環(huán)。expression_list是被應(yīng)用到每個(gè)集成員的表達(dá)式列表，當(dāng)用的是for函數(shù)時(shí)，expression_list可以包含多個(gè)表達(dá)式，其間用逗號(hào)隔開(kāi)。這些表達(dá)式將被作為約束加到模型中。當(dāng)使用其余的三個(gè)集循環(huán)函數(shù)時(shí)，expression_list只能有一個(gè)表達(dá)式。如果省略set_index_list

29、，那么在expression_list中引用的所有屬性的類(lèi)型都是setname集。FOR( setname ( set_index_list) | cond_qualifier: exp_list)This generates the expressions contained in exp_list for all members of the setname set. FOR is used to generate constraints over all the members of a set, and is one of the most powerful features of L

30、INGO.MAX( setname ( set_index_list) | cond_qualifier: expression) This returns the maximum value of expression taken over the setname set.MIN( setname ( set_index_list) | cond_qualifier: expression) This returns the minimum value of expression taken over the setname set.PROD( setname ( set_index_lis

31、t) | cond_qualifier: expression) This returns the product of an expression over the setname set.SUM( setname ( set_index_list) | cond_qualifier: expression)This returns the sum of an expression over the setname set.Set looping functions are discussed in more detail in section Set Looping Functions i

32、n Chapter 2.3.8 輔助功能函數(shù)（Miscellaneous Functions）IF( logical_condition, true_result, false_result)The IF function evaluates logical_condition and, if true, returns true_result, otherwise it returns false_result. For example, consider the following simple model that uses IF to compute fixed production

33、costs:MIN = COST;COST = XCOST + YCOST;XCOST = IF( X #GT# 0, 100, 0) + 2 * X;YCOST = IF( Y #GT# 0, 60, 0) + 3 * Y;X + Y = 30;Model: IFCOSTWe produce two products梄 and Y. We want to minimize total cost, subject to producing at least 30 total units of X and Y. If we produce X, there is a fixed charge o

34、f 100 along with a variable cost of 2. Similarly, for Y, these respective values are 60 and 3. We use the IF function to determine if either of the products are being produced in order to apply the relevant fixed cost. This is accomplished by testing to see if their production levels are greater tha

35、n 0. If so, we return the fixed cost value, otherwise, we return zero.Experienced modelers know that, without the benefit of an IF function, modeling fixed costs requires invoking some tricks using binary integer variables. The resulting models are not as intuitive as models constructed using IF. Th

36、e caveat, however, is that the IF function is not a linear function. At best, the graph of an IF function will be piecewise linear. In our current example, the IF functions are piecewise linear with a discontinuous break at the origin. As we discuss in the chapter On Mathematical Modeling, it is alw

37、ays best to try and keep a model linear. Barring this, it is best for all functions in a nonlinear model to be continuous. Clearly, then, the IF function is a problem in that it violates both these conditions. Thus, models containing IF functions may be tough to solve to global optimality. Fortunate

38、ly, LINGO has two options that can help overcome the difficult nature of models containing IF functions條inearization and global optimization.To illustrate the difficulty in solving models with discontinuous functions such as IF, we will solve our example model with both linearization and global opti

39、mization disabled. When we do this, we get the following solution:Local optimal solution found at iteration: 42Objective value: 160.0000 Variable Value COST 160.0000 XCOST 160.0000 YCOST 0.000000 X 30.00000 Y 0.000000This solution involves producing only X at a total cost of 160. Clearly, this is me

40、rely a locally optimal point, given that producing only Y and not X will result in a lower total cost of 150. In order to find the globally optimal point we must resort to either the linearization or global optimization features in LINGO. Briefly, linearization seeks to reformulate a nonlinear model

41、 into a mathematically equivalent linear model. This is desirable for two reasons. First, and most important, linear models can always be solved to global optimality. Secondly, linear models will tend to solve much faster than equivalent nonlinear models. Unfortunately, linearization can抰 always tra

42、nsform a model into an equivalent linear state, in which case, it may be of no benefit. Fortunately, our sample model can be entirely linearized. To enable the linearization option, run the LINGO|Options command and set the Linearization Degree to High on the General Solver tab. Global optimization

43、breaks a model down into a series of smaller, local models. Once this series of local models has been solved, a globally optimal solution can be determined. To enable global optimization, run the LINGO|Options command, select the Global Solver tab, then click on the Global Solver checkbox. Note that

44、 the global solver is an add-on option to LINGO. The global solver feature will not be enabled for some installations. Run the Help|About LINGO command to determine if your installation has the global solver capability enabled.Regardless, whether using the linearization option or the global solver,

45、LINGO obtains the true, global solution:Global optimal solution found at iteration: 6Objective value: 150.0000 Variable Value COST 150.0000 XCOST 0.000000 YCOST 150.0000 X 0.000000 Y 30.00000Note:Starting with release 9.0, the false branch of the IF function may contain arithmetic errors without cau

46、sing the solver to trigger an error. This makes the IF function useful in avoiding problems when the solver strays into areas where certain functions become undefined. For instance, if your model involves division by a variable, you might use IF as follows: IF( X #GT# 1.E-10, 1/X, 1.E10).WARN( text,

47、 logical_condition)This function displays the message 憈ext?if the logical_condition is met. This feature is useful for verifying the validity of a models data. In the following example, if the user has entered a negative interest rate, the message INVALID營(yíng)NTEREST燫ATE is displayed:! A model of a home

48、 mortgage;DATA:! Prompt the user for the interest rate, years, and value of mortgage. We will compute the monthly payment; YRATE = ?; YEARS = ?; LUMP = ?;ENDDATA! Number of monthly payment; MONTHS = YEARS * 12;! Monthly interest rate; ( 1 + MRATE) 12 = 1 + YRATE;! Solve next line for monthly payment

49、; LUMP = PAYMENT * FPA( MRATE, MONTHS);! Warn them if interest rate is negative WARN( INVALID INTEREST RATE, YRATE #LT# 0);USER( user_determined_arguments)The user can supply this in an external DLL or object code file. For a detailed example on the use of USER, see User Defined Functions.4.靈敏度分析i.e

50、：某家具公司制造書(shū)桌、餐桌和椅子，所用的資源有三種：木料、木工和漆工。生產(chǎn)數(shù)據(jù)如下表所示：每個(gè)書(shū)桌每個(gè)餐桌每個(gè)椅子現(xiàn)有資源總數(shù)木料8單位6單位1單位48單位漆工4單位2單位1.5單位20單位木工2單位1.5單位0.5單位8單位成品單價(jià)60單位30單位20單位若要求桌子的生產(chǎn)量不超過(guò)5件，如何安排三種產(chǎn)品的生產(chǎn)可使利潤(rùn)最大？max=60*desks+30*tables+20*chairs;8*desks+6*tables+chairs=48;4*desks+2*tables+1.5*chairs=20;2*desks+1.5*tables+.5*chairs=8;tables=5; Global

51、 optimal solution found. Objective value: 280.0000 Total solver iterations: 3 Variable Value Reduced Cost DESKS 2.000000 0.000000 TABLES 0.000000 5.000000 CHAIRS 8.000000 0.000000 Row Slack or Surplus Dual Price 1 280.0000 1.000000 2 24.00000 0.000000 3 0.000000 10.00000 4 0.000000 10.00000 5 5.0000

52、00 0.000000“Slack or Surplus”給出松馳變量的值：第1行松馳變量 =280（模型第一行表示目標(biāo)函數(shù)，所以第二行對(duì)應(yīng)第一個(gè)約束）第2行松馳變量 =24第3行松馳變量 =0第4行松馳變量 =0第5行松馳變量 =5“Reduced Cost”列出最優(yōu)單純形表中判別數(shù)所在行的變量的系數(shù)，表示當(dāng)變量有微小變動(dòng)時(shí), 目標(biāo)函數(shù)的變化率。其中基變量的reduced cost值應(yīng)為0， 對(duì)于非基變量 Xj, 相應(yīng)的 reduced cost值表示當(dāng)某個(gè)變量Xj 增加一個(gè)單位時(shí)目標(biāo)函數(shù)減少的量( max型問(wèn)題)。本例中：變量tables對(duì)應(yīng)的reduced cost值為5，表示當(dāng)非基變量

53、tables的值從0變?yōu)?1時(shí)（此時(shí)假定其他非基變量保持不變,但為了滿足約束條件，基變量顯然會(huì)發(fā)生變化），最優(yōu)的目標(biāo)函數(shù)值 = 280 - 5 = 275?！癉UAL PRICE”（對(duì)偶價(jià)格）表示當(dāng)對(duì)應(yīng)約束有微小變動(dòng)時(shí), 目標(biāo)函數(shù)的變化率。輸出結(jié)果中對(duì)應(yīng)于每一個(gè)約束有一個(gè)對(duì)偶價(jià)格。 若其數(shù)值為p， 表示對(duì)應(yīng)約束中不等式右端項(xiàng)若增加1 個(gè)單位，目標(biāo)函數(shù)將增加p個(gè)單位（max型問(wèn)題）。顯然，如果在最優(yōu)解處約束正好取等號(hào)（也就是“緊約束”，也稱為有效約束或起作用約束），對(duì)偶價(jià)格值才可能不是0。本例中：第3、4行是緊約束，對(duì)應(yīng)的對(duì)偶價(jià)格值為10，表示當(dāng)緊約束 3) 4 DESKS + 2 TABLES

54、 + 1.5 CHAIRS = 20 變?yōu)?3) 4 DESKS + 2 TABLES + 1.5 CHAIRS 1的正整數(shù)）：N點(diǎn)求解 Barrier: 障礙法 (即內(nèi)點(diǎn)法)6.綜合案例旅行售貨員問(wèn)題（Traveling Salesman Problem）有一個(gè)推銷(xiāo)員，從城市1出發(fā)，要遍訪城市2，3，n各一次，最后返回城市1。已知從城市i到j(luò)的旅費(fèi)為，問(wèn)他應(yīng)按怎樣的次序訪問(wèn)這些城市，使得總旅費(fèi)最少？可以用多種方法把TSP表示成整數(shù)規(guī)劃模型。這里介紹的一種建立模型的方法，是把該問(wèn)題的每個(gè)解（不一定是最優(yōu)的）看作是一次“巡回”。在下述意義下，引入一些0-1整數(shù)變量：其目標(biāo)只是使為最小。這里有兩個(gè)

55、明顯的必須滿足的條件：訪問(wèn)城市i后必須要有一個(gè)即將訪問(wèn)的確切城市；訪問(wèn)城市j前必須要有一個(gè)剛剛訪問(wèn)過(guò)的確切城市。用下面的兩組約束分別實(shí)現(xiàn)上面的兩個(gè)條件。123456到此我們得到了一個(gè)模型，它是一個(gè)指派問(wèn)題的整數(shù)規(guī)劃模型。但以上兩個(gè)條件對(duì)于TSP來(lái)說(shuō)并不充分，僅僅是必要條件。例如：以上兩個(gè)條件都滿足，但它顯然不是TSP的解，它存在兩個(gè)子巡回。這里，我們將敘述一種在原模型上附加充分的約束條件以避免產(chǎn)生子巡回的方法。把額外變量附加到問(wèn)題中?？砂堰@些變量看作是連續(xù)的（最然這些變量在最優(yōu)解中取普通的整數(shù)值）?，F(xiàn)在附加下面形式的約束條件。 為了證明該約束條件有預(yù)期的效果，必須證明：（1）任何含子巡回的路線

56、都不滿足該約束條件；（2）全部巡回都滿足該約束條件。首先證明（1），用反證法。假設(shè)還存在子巡回，也就是說(shuō)至少有兩個(gè)子巡回。那么至少存在一個(gè)子巡回中不含城市1。把該子巡回記為，則必有把這k個(gè)式子相加，有，矛盾！故假設(shè)不正確，結(jié)論（1）得證。下面證明（2），采用構(gòu)造法。對(duì)于任意的總巡回，可取訪問(wèn)城市i的順序數(shù)，取值范圍為。因此，。下面來(lái)證明總巡回滿足該約束條件。()總巡回上的邊（）非總巡回上的邊從而結(jié)論（2）得證。這樣我們把TSP轉(zhuǎn)化成了一個(gè)混合整數(shù)線性規(guī)劃問(wèn)題。顯然，當(dāng)城市個(gè)數(shù)較大（大于30）時(shí)，該混合整數(shù)線性規(guī)劃問(wèn)題的規(guī)模會(huì)很大，從而給求解帶來(lái)很大問(wèn)題。TSP已被證明是NP難問(wèn)題，目前還沒(méi)有發(fā)

57、現(xiàn)多項(xiàng)式時(shí)間的算法。對(duì)于小規(guī)模問(wèn)題，我們求解這個(gè)混合整數(shù)線性規(guī)劃問(wèn)題的方式還是有效的。!旅行售貨員問(wèn)題;model:sets: city / 1. 5/: u; link( city, city): dist, ! 距離矩陣; x;endsets n = size( city);data: !距離矩陣，它并不需要是對(duì)稱的; dist = 0 2 6 7 4 2 0 5 3 8 6 5 0 9 7 7 3 9 0 9 4 8 7 9 0 ; !隨機(jī)產(chǎn)生，這里可改為你要解決的問(wèn)題的數(shù)據(jù);enddata !目標(biāo)函數(shù); min = sum( link: dist * x); FOR( city( K)

58、: !進(jìn)入城市K; sum( city( I)| I #ne# K: x( I, K) = 1; !離開(kāi)城市K; sum( city( J)| J #ne# K: x( K, J) = 1; ); !保證不出現(xiàn)子圈; for(city(I)|I #gt# 1: for( city( J)| J#gt#1 #and# I #ne# J: u(I)-u(J)+n*x(I,J)=n-1); ); !限制u的范圍以加速模型的求解，保證所加限制并不排除掉TSP問(wèn)題的最優(yōu)解; for(city(I) | I #gt# 1: u(I)=n-2 ); !定義X為01變量; for( link: bin( x)

59、;endGlobal optimal solution found. Objective value: 25.00000 Extended solver steps: 0 Total solver iterations: 19 Variable Value Reduced Cost N 5.000000 0.000000 U( 1) 0.000000 0.000000 U( 2) 0.000000 0.000000 U( 3) 2.000000 0.000000 U( 4) 1.000000 0.000000 U( 5) 3.000000 0.000000 DIST( 1, 1) 0.000000 0.000000 DIST( 1, 2) 2.000000 0.000000 DIST( 1, 3) 6.000000 0.000000 DIST( 1, 4) 7.000000 0.000000 DIST( 1, 5) 4.000000 0.0000