# Mini MILP "benchmark"

a guest
Sep 15th, 2021
83
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. using JuMP
2. import GLPK
3. import Cbc
4.
5. function example_urban_plan(factory)
6. model = Model(factory)
7. set_silent(model)
8. # x is indexed by row and column
9. @variable(model, 0 <= x[1:5, 1:5] <= 1, Int)
10. # y is indexed by R or C, the points, and an index in 1:5. Note how JuMP
11. # allows indexing on arbitrary sets.
12. rowcol = ["R", "C"]
13. points = [5, 4, 3, -3, -4, -5]
14. @variable(model, 0 <= y[rowcol, points, 1:5] <= 1, Int)
15. # Objective - combine the positive and negative parts
16. @objective(
17. model,
18. Max,
19. sum(
20. 3 * (y["R", 3, i] + y["C", 3, i]) +
21. 1 * (y["R", 4, i] + y["C", 4, i]) +
22. 1 * (y["R", 5, i] + y["C", 5, i]) -
23. 3 * (y["R", -3, i] + y["C", -3, i]) -
24. 1 * (y["R", -4, i] + y["C", -4, i]) -
25. 1 * (y["R", -5, i] + y["C", -5, i]) for i in 1:5
26. )
27. )
28. # Constrain the number of residential lots
29. @constraint(model, sum(x) == 12)
30. # Add the constraints that link the auxiliary y variables to the x variables
31. for i in 1:5
32. @constraints(model, begin
33. # Rows
34. y["R", 5, i] <= 1 / 5 * sum(x[i, :]) # sum = 5
35. y["R", 4, i] <= 1 / 4 * sum(x[i, :]) # sum = 4
36. y["R", 3, i] <= 1 / 3 * sum(x[i, :]) # sum = 3
37. y["R", -3, i] >= 1 - 1 / 3 * sum(x[i, :]) # sum = 2
38. y["R", -4, i] >= 1 - 1 / 2 * sum(x[i, :]) # sum = 1
39. y["R", -5, i] >= 1 - 1 / 1 * sum(x[i, :]) # sum = 0
40. # Columns
41. y["C", 5, i] <= 1 / 5 * sum(x[:, i]) # sum = 5
42. y["C", 4, i] <= 1 / 4 * sum(x[:, i]) # sum = 4
43. y["C", 3, i] <= 1 / 3 * sum(x[:, i]) # sum = 3
44. y["C", -3, i] >= 1 - 1 / 3 * sum(x[:, i]) # sum = 2
45. y["C", -4, i] >= 1 - 1 / 2 * sum(x[:, i]) # sum = 1
46. y["C", -5, i] >= 1 - 1 / 1 * sum(x[:, i]) # sum = 0
47. end)
48. end
49. # Solve it
50. optimize!(model)
51. if termination_status(model) == MOI.OPTIMAL
52. println("Urban_plan ",solver_name(model)," ", solve_time(model))
53. else
54. println("Urban_plan ",solver_name(model)," failed")
55. end
56. return
57. end
58.
59. function example_transp(factory)
60. ORIG = ["GARY", "CLEV", "PITT"]
61. DEST = ["FRA", "DET", "LAN", "WIN", "STL", "FRE", "LAF"]
62. supply = [1_400, 2_600, 2_900]
63. demand = [900, 1_200, 600, 400, 1_700, 1_100, 1_000]
64. cost = [
65. 39 14 11 14 16 82 8
66. 27 9 12 9 26 95 17
67. 24 14 17 13 28 99 20
68. ]
69. model = Model(factory)
70. set_silent(model)
71. @variable(model, trans[1:length(ORIG), 1:length(DEST)] >= 0)
72. @objective(
73. model,
74. Min,
75. sum(
76. cost[i, j] * trans[i, j] for i in 1:length(ORIG),
77. j in 1:length(DEST)
78. )
79. )
80. @constraints(
81. model,
82. begin
83. [i in 1:length(ORIG)], sum(trans[i, :]) == supply[i]
84. [j in 1:length(DEST)], sum(trans[:, j]) == demand[j]
85. end
86. )
87. optimize!(model)
88. if termination_status(model) == MOI.OPTIMAL
89. println("transp ",solver_name(model)," ", solve_time(model))
90. else
91. println("transp ",solver_name(model)," failed")
92. end
93. return
94. end
95.
96. function example_steelT3(factory; verbose = false)
97. T = 4
98. prod = ["bands", "coils"]
99. area = Dict(
100. "bands" => ("east", "north"),
101. "coils" => ("east", "west", "export"),
102. )
103. avail = [40, 40, 32, 40]
104. rate = Dict("bands" => 200, "coils" => 140)
105. inv0 = Dict("bands" => 10, "coils" => 0)
106. prodcost = Dict("bands" => 10, "coils" => 11)
107. invcost = Dict("bands" => 2.5, "coils" => 3)
108. revenue = Dict(
109. "bands" => Dict(
110. "east" => [25.0, 26.0, 27.0, 27.0],
111. "north" => [26.5, 27.5, 28.0, 28.5],
112. ),
113. "coils" => Dict(
114. "east" => [30, 35, 37, 39],
115. "west" => [29, 32, 33, 35],
116. "export" => [25, 25, 25, 28],
117. ),
118. )
119. market = Dict(
120. "bands" => Dict(
121. "east" => [2000, 2000, 1500, 2000],
122. "north" => [4000, 4000, 2500, 4500],
123. ),
124. "coils" => Dict(
125. "east" => [1000, 800, 1000, 1100],
126. "west" => [2000, 1200, 2000, 2300],
127. "export" => [1000, 500, 500, 800],
128. ),
129. )
130. # Model
131. model = Model(factory)
132. set_silent(model)
133. # Decision Variables
134. @variables(
135. model,
136. begin
137. make[p in prod, t in 1:T] >= 0
138. inventory[p in prod, t in 0:T] >= 0
139. 0 <= sell[p in prod, a in area[p], t in 1:T] <= market[p][a][t]
140. end
141. )
142. @constraints(
143. model,
144. begin
145. [p = prod, a = area[p], t = 1:T], sell[p, a, t] <= market[p][a][t]
146. # Total of hours used by all products may not exceed hours available,
147. # in each week
148. [t in 1:T], sum(1 / rate[p] * make[p, t] for p in prod) <= avail[t]
149. # Initial inventory must equal given value
150. [p in prod], inventory[p, 0] == inv0[p]
151. # Tons produced and taken from inventory must equal tons sold and put
152. # into inventory.
153. [p in prod, t in 1:T],
154. make[p, t] + inventory[p, t-1] ==
155. sum(sell[p, a, t] for a in area[p]) + inventory[p, t]
156. end
157. )
158. # Maximize total profit: total revenue less costs for all products in all
159. # weeks.
160. @objective(
161. model,
162. Max,
163. sum(
164. revenue[p][a][t] * sell[p, a, t] - prodcost[p] * make[p, t] -
165. invcost[p] * inventory[p, t] for p in prod, a in area[p], t in 1:T
166. )
167. )
168. optimize!(model)
169. if termination_status(model) == MOI.OPTIMAL
170. println("SteelT3 ",solver_name(model)," ", solve_time(model))
171. else
172. println("SteelT3 ",solver_name(model)," failed")
173. end
174.
175. return
176. end
177.
178. function example_prod(factory; verbose = false)
179. # PRODUCTION SETS AND PARAMETERS
180. prd = ["18REG" "24REG" "24PRO"]
181. # Members of the product group
182. numprd = length(prd)
183. pt = [1.194, 1.509, 1.509]
184. # Crew-hours to produce 1000 units
185. pc = [2304, 2920, 2910]
186. # Nominal production cost per 1000, used
187. # to compute inventory and shortage costs
188. #
189. # TIME PERIOD SETS AND PARAMETERS
190. firstperiod = 1
191. # Index of first production period to be modeled
192. lastperiod = 13
193. # Index of last production period to be modeled
194. numperiods = firstperiod:lastperiod
195. # 'planning horizon' := first..last;
196. # EMPLOYMENT PARAMETERS
197. # Workers per crew
198. cs = 18
199. # Regular-time hours per shift
200. sl = 8
201. # Wage per hour for regular-time labor
202. rtr = 16.00
203. # Wage per hour for overtime labor
204. otr = 43.85
205. # Crews employed at start of first period
206. iw = 8
207. # Regular working days in a production period
208. dpp = [19.5, 19, 20, 19, 19.5, 19, 19, 20, 19, 20, 20, 18, 18]
209. # Maximum crew-hours of overtime in a period
210. ol = [96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96, 96]
211. # Lower limit on average employment in a period
212. cmin = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
213. # Upper limit on average employment in a period
214. cmax = [8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8]
215. # Penalty cost of hiring a crew
216. hc = [
217. 7500,
218. 7500,
219. 7500,
220. 7500,
221. 15000,
222. 15000,
223. 15000,
224. 15000,
225. 15000,
226. 15000,
227. 7500,
228. 7500,
229. 7500,
230. ]
231. # Penalty cost of laying off a crew
232. lc = [
233. 7500,
234. 7500,
235. 7500,
236. 7500,
237. 15000,
238. 15000,
239. 15000,
240. 15000,
241. 15000,
242. 15000,
243. 7500,
244. 7500,
245. 7500,
246. ]
247. # DEMAND PARAMETERS
248. d18REG = [
249. 63.8,
250. 76,
251. 88.4,
252. 913.8,
253. 115,
254. 133.8,
255. 79.6,
256. 111,
257. 121.6,
258. 470,
259. 78.4,
260. 99.4,
261. 140.4,
262. 63.8,
263. ]
264. d24REG = [
265. 1212,
266. 306.2,
267. 319,
268. 208.4,
269. 298,
270. 328.2,
271. 959.6,
272. 257.6,
273. 335.6,
274. 118,
275. 284.8,
276. 970,
277. 343.8,
278. 1212,
279. ]
280. d24PRO = [0, 0, 0, 0, 0, 0, 0, 0, 0, 1102, 0, 0, 0, 0]
281. # Requirements (in 1000s) to be met from current production and inventory
282. dem = Array[d18REG, d24REG, d24PRO]
283. # true if product will be the subject of a special promotion in the period
284. pro = Array[
285. [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
286. [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 1],
287. [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
288. ]
289. # INVENTORY AND SHORTAGE PARAMETERS
290. # Proportion of non-promoted demand that must be in inventory the previous
291. # period
292. rir = 0.75
293. # Proportion of promoted demand that must be in inventory the previous
294. # period
295. pir = 0.80
296. # Upper limit on number of periods that any product may sit in inventory
297. life = 2
298. # Inventory cost per 1000 units is cri times nominal production cost
299. cri = [0.015, 0.015, 0.015]
300. # Shortage cost per 1000 units is crs times nominal production cost
301. crs = [1.1, 1.1, 1.1]
302. # Inventory at start of first period; age unknown
303. iinv = [82, 792.2, 0]
304. # Initial inventory still available for allocation at end of period t
305. iil = [
306. [
307. max(0, iinv[p] - sum(dem[p][v] for v in firstperiod:t)) for
308. t in numperiods
309. ] for p in 1:numprd
310. ]
311. # Lower limit on inventory at end of period t
312. function checkpro(
313. product,
314. timeperiod,
315. production,
316. promotionalrate,
317. regularrate,
318. )
319. if production[product][timeperiod+1] == 1
320. return promotionalrate
321. else
322. return regularrate
323. end
324. end
325. minv = [
326. [dem[p][t+1] * checkpro(p, t, pro, pir, rir) for t in numperiods]
327. for p in 1:numprd
328. ]
329. # DEFINE MODEL
330. prod = Model(factory)
331. set_silent(prod)
332. # VARIABLES
333. # Average number of crews employed in each period
334. @variable(prod, Crews[0:lastperiod] >= 0)
335. # Crews hired from previous to current period
336. @variable(prod, Hire[numperiods] >= 0)
337. # Crews laid off from previous to current period
338. @variable(prod, Layoff[numperiods] >= 0)
339. # Production using regular-time labor, in 1000s
340. @variable(prod, Rprd[1:numprd, numperiods] >= 0)
341. # Production using overtime labor, in 1000s
342. @variable(prod, Oprd[1:numprd, numperiods] >= 0)
343. # a numperiods old -- produced in period (t+1)-a --
344. # and still in storage at the end of period t
345. @variable(prod, Inv[1:numprd, numperiods, 1:life] >= 0)
346. # Accumulated unsatisfied demand at the end of period t
347. @variable(prod, Short[1:numprd, numperiods] >= 0)
348. # CONSTRAINTS
349. # Hours needed to accomplish all regular-time production in a period must
350. # not exceed hours available on all shifts
351. @constraint(
352. prod,
353. [t = numperiods],
354. sum(pt[p] * Rprd[p, t] for p in 1:numprd) <= sl * dpp[t] * Crews[t]
355. )
356. # Hours needed to accomplish all overtime production in a period must not
357. # exceed the specified overtime limit
358. @constraint(
359. prod,
360. [t = numperiods],
361. sum(pt[p] * Oprd[p, t] for p in 1:numprd) <= ol[t]
362. )
363. # Use given initial workforce
364. @constraint(prod, Crews[firstperiod-1] == iw)
365. # Workforce changes by hiring or layoffs
366. @constraint(
367. prod,
368. [t in numperiods],
369. Crews[t] == Crews[t-1] + Hire[t] - Layoff[t]
370. )
371. # Workforce must remain within specified bounds
372. @constraint(prod, [t in numperiods], cmin[t] <= Crews[t])
373. @constraint(prod, [t in numperiods], Crews[t] <= cmax[t])
374. # 'first demand requirement
375. @constraint(
376. prod,
377. [p in 1:numprd],
378. Rprd[p, firstperiod] + Oprd[p, firstperiod] + Short[p, firstperiod] -
379. Inv[p, firstperiod, 1] == max(0, dem[p][firstperiod] - iinv[p])
380. )
381. # Production plus increase in shortage plus decrease in inventory must
382. # equal demand
383. for t in (firstperiod+1):lastperiod
384. @constraint(
385. prod,
386. [p in 1:numprd],
387. Rprd[p, t] + Oprd[p, t] + Short[p, t] - Short[p, t-1] +
388. sum(Inv[p, t-1, a] - Inv[p, t, a] for a in 1:life) ==
389. max(0, dem[p][t] - iil[p][t-1])
390. )
391. end
392. # Inventory in storage at end of period t must meet specified minimum
393. @constraint(
394. prod,
395. [p in 1:numprd, t in numperiods],
396. sum(Inv[p, t, a] + iil[p][t] for a in 1:life) >= minv[p][t]
397. )
398. # In the vth period (starting from first) no inventory may be more than v
399. # numperiods old (initial inventories are handled separately)
400. @constraint(
401. prod,
402. [p in 1:numprd, v in 1:(life-1), a in (v+1):life],
403. Inv[p, firstperiod+v-1, a] == 0
404. )
405. # New inventory cannot exceed production in the most recent period
406. @constraint(
407. prod,
408. [p in 1:numprd, t in numperiods],
409. Inv[p, t, 1] <= Rprd[p, t] + Oprd[p, t]
410. )
411. # Inventory left from period (t+1)-p can only decrease as time goes on
412. secondperiod = firstperiod + 1
413. @constraint(
414. prod,
415. [p in 1:numprd, t in 2:lastperiod, a in 2:life],
416. Inv[p, t, a] <= Inv[p, t-1, a-1]
417. )
418. # OBJECTIVE
419. # Full regular wages for all crews employed, plus penalties for hiring and
420. # layoffs, plus wages for any overtime worked, plus inventory and shortage
421. # costs. (All other production costs are assumed to depend on initial
422. # inventory and on demands, and so are not included explicitly.)
423. @objective(
424. prod,
425. Min,
426. sum(
427. rtr * sl * dpp[t] * cs * Crews[t] +
428. hc[t] * Hire[t] +
429. lc[t] * Layoff[t] +
430. sum(
431. otr * cs * pt[p] * Oprd[p, t] +
432. sum(cri[p] * pc[p] * Inv[p, t, a] for a in 1:life) +
433. crs[p] * pc[p] * Short[p, t] for p in 1:numprd
434. ) for t in numperiods
435. )
436. )
437. # Obtain solution
438. optimize!(prod)
439. model = prod
440. if termination_status(model) == MOI.OPTIMAL
441. println("prod ",solver_name(model)," ", solve_time(model))
442. else
443. println("prod ",solver_name(model)," failed")
444. end
445.
446. return
447. end
448.
449. function example_knapsack(factory; verbose = false)
450. profit = [5, 3, 2, 7, 4]
451. weight = [2, 8, 4, 2, 5]
452. capacity = 10
453. model = Model(factory)
454. set_silent(model)
455. @variable(model, x[1:5], Bin)
456. # Objective: maximize profit
457. @objective(model, Max, profit' * x)
458. # Constraint: can carry all
459. @constraint(model, weight' * x <= capacity)
460. # Solve problem using MIP solver
461. optimize!(model)
462. if termination_status(model) == MOI.OPTIMAL
463. println("knapsack ",solver_name(model)," ", solve_time(model))
464. else
465. println("knapsack ",solver_name(model)," failed")
466. end
467. return
468. end
469.
470. function example_factory_schedule(factory)
471. # Sets in the problem:
472. months, factories = 1:12, [:A, :B]
473. # This function takes a matrix and converts it to a JuMP container so we can
474. # refer to elements such as `d_max_cap[1, :A]`.
475. containerize(A::Matrix) = Containers.DenseAxisArray(A, months, factories)
476. # Maximum production capacity in (month, factory) [units/month]:
477. d_max_cap = containerize(
478. [
479. 100000 50000
480. 110000 55000
481. 120000 60000
482. 145000 100000
483. 160000 0
484. 140000 70000
485. 155000 60000
486. 200000 100000
487. 210000 100000
488. 197000 100000
489. 80000 120000
490. 150000 150000
491. ],
492. )
493. # Minimum production capacity in (month, factory) [units/month]:
494. d_min_cap = containerize(
495. [
496. 20000 20000
497. 20000 20000
498. 20000 20000
499. 20000 20000
500. 20000 0
501. 20000 20000
502. 20000 20000
503. 20000 20000
504. 20000 20000
505. 20000 20000
506. 20000 20000
507. 20000 20000
508. ],
509. )
510. # Variable cost of production in (month, factory) [\$/unit]:
511. d_var_cost = containerize([
512. 10 5
513. 11 4
514. 12 3
515. 9 5
516. 8 0
517. 8 6
518. 5 4
519. 7 6
520. 9 8
521. 10 11
522. 8 10
523. 8 12
524. ])
525. # Fixed cost of production in (month, factory) # [\$/month]:
526. d_fixed_cost = containerize(
527. [
528. 500 600
529. 500 600
530. 500 600
531. 500 600
532. 500 0
533. 500 600
534. 500 600
535. 500 600
536. 500 600
537. 500 600
538. 500 600
539. 500 600
540. ],
541. )
542. # Demand in each month [units/month]:
543. d_demand = [
544. 120_000,
545. 100_000,
546. 130_000,
547. 130_000,
548. 140_000,
549. 130_000,
550. 150_000,
551. 170_000,
552. 200_000,
553. 190_000,
554. 140_000,
555. 100_000,
556. ]
557. # The model!
558. model = Model(factory)
559. set_silent(model)
560. # Decision variables
561. @variables(model, begin
562. status[m in months, f in factories], Bin
563. production[m in months, f in factories], Int
564. end)
565. # The production cannot be less than minimum capacity.
566. @constraint(
567. model,
568. [m in months, f in factories],
569. production[m, f] >= d_min_cap[m, f] * status[m, f],
570. )
571. # The production cannot be more that maximum capacity.
572. @constraint(
573. model,
574. [m in months, f in factories],
575. production[m, f] <= d_max_cap[m, f] * status[m, f],
576. )
577. # The production must equal demand in a given month.
578. @constraint(model, [m in months], sum(production[m, :]) == d_demand[m])
579. # Factory B is shut down during month 5, so production and status are both
580. # zero.
581. fix(status[5, :B], 0.0)
582. fix(production[5, :B], 0.0)
583. # The objective is to minimize the cost of production across all time
584. ##periods.
585. @objective(
586. model,
587. Min,
588. sum(
589. d_fixed_cost[m, f] * status[m, f] +
590. d_var_cost[m, f] * production[m, f] for m in months, f in factories
591. )
592. )
593. # Optimize the problem
594. optimize!(model)
595. if termination_status(model) == MOI.OPTIMAL
596. println("factory_schedule ",solver_name(model)," ", solve_time(model))
597. else
598. println("factory_schedule ",solver_name(model)," failed")
599. end
600. return
601. end
602.
603. for p in [example_factory_schedule, example_knapsack, example_prod, example_steelT3, example_transp, example_urban_plan ]
604. for s in [GLPK.Optimizer, with_optimizer(Cbc.Optimizer, logLevel = 0)]
605. for trial in 1:5
606. p(s)
607. end
608. end
609. end
610.
RAW Paste Data