Guest User

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