❮❮❮ Multicommodity network flow model of a human resource allocation problem considering time periods ❮❮❮ ❯❯❯ Disassembly line optimization with reinforcement learning ❯❯❯

MILP models of a patient transportation problem

Tags: OR Algebraic-Modelling Reading-Notes Cejor

Reading time: 13 min.

Metadata

Tóth, Martin; Hajba, Tamás; Horváth, Adrián MILP models of a patient transportation problem. Cent. Eur. J. Oper. Res. 32, No. 4, 903-922 (2024).

Links: doi zbMATH

Authors:

corresponding: Tamás Hajba hajbat@math.sze.hu
math department
hp

Summary

Questions

modelling and benchmarking
- value of $M$
open implementation
- sourcecode
- data files

Notes

Glossary

DARP = Dial a ride problem
TSP = travelling salesman problem
MTZ = classic formulation of the TSP

Parameters

$N\ge{}1$ by (3)!!
$K\ge{}2$ by (7) and (8)
travel distance in $D$ does not have to be symmetric
there is an implicit constraint of $K\ge{}2$ (because there is a last round)
there is no unboarding time at all!
there should be $N\ge{}1$, otherwise no problem
introducing an index set of patient types would unite a lot of parameters, decision variables and constraints

Decision variables

$pn$ and $dpn$ are not continuous but in $\mathbb{N}$ by (12)-(17)

Pickup locations vs physical locations

pickup locations are not 1:1 physical locations
there can be several pickups at the same physical location of different groups of patients with different properties
- distance $D$ does not change, but all other $i$ indexed parameters
- in particular traveling time $T$ (depending on time of the day and traffic etc)

Closing round

the closing round could be eliminated
then depot closure time would be $B_0 - D_{N+1,0}$
the objective function shifts by the constant $D_{N+1, 0}$
only $K\ge{}1$ necessary
- mmh, if $N=0$, then still $K\ge{}1$?
(11) could be eliminated
indices in (8) an (9) could be simplified

Opening round

if the vehicle would travel depot -> hospital directly in round 1, the formulation would be much easier
but this would change the problem

Leaving out the depot

possible at the end, as the last leg is hospital -> depot
- this eleminates the closing round in model 1
- this would need a Hamiltonian path instead of a Hamiltonian cycle in model 2
but because of pickup of first patient impossible to eliminate in the beginning

Constraint indices in model 1

these jump between $i$ and $j$ without reason

Data in model 2

is the same as in model 1, but needs to be adapted
see reformulation of model 2

TSP in model 2

the indices constraints (37)-(40) are off by one
(39) is already part of the domain bound of $u$ (and off by one)
see the interlude model below

Pickup-Dropoff order in model 2

Error in paper: not (50) and (51), but (40) and (41) force dropoffs after pickups
There is no proof for this order in the paper
A combinatorial argument:
- $p_i$ and $h_i$ (patient pickup and hospital dropoff of patient group $i$)
- Route is string of letters
- Invariants: maximal $p$ substring $p_{i_1}\dotso{}p_{i_k}$ is directly followed by $h_{i_1}\dotso{}h_{i_k}$
- Case $N=0$: $$
- Case $N=1$: $p_1 h_1$
- Case $N=2$: $p_1 h_1 p_2 h_2 d$ and $d h_1 p_2 h_1 h_2$ and their permutations of the patient groups
- Case $N\ge{}2$:
  - By recursion
  - Insert $p_N$ first
    - Can not be inside $h$ substring => this would break invariant of corresponding $p$ substring to the left
    - Inside or adjoining a $p$ substring => insert $h_N$ at same relative position in corresponding $h$ substring to the right

Model strengthening

(60) is not a valid constraint for model 2 as written
- it should be $x_{ij}=x_{ji}=0$ (as there are no explicit rounds in the graph reformulation)

Big M

in (16)-(17), (20)-(23) and (30)-(31) there is a common, unspecified $M$
same in (47)-(48) and (51)-(53)
Probably part of a big-M formulation
no value given
should be taylored to constraint and index!
- otherwise really bad for benchmarking!

Model size comparison

the model size comparison in table 1 tells a lot
- model 1 is just to sparse
but setting $K=N+1$ really hurts model 1
- an a priori heuristic solution should give a much tigher bound than that
- e.g. a greedy heuristic by picking up patients by arrival time
- get some estimate for $K$, try a solution => timeout, then double $K$

Benchmarking results

The Big-M issue makes it unclear what was really benchmarked
- different M for each equation at index
- different M for each equation
- different M for each model
- same M globally for both models
there seems to be only a single benchmark run per scenario
- no truncated means/shifted geometric means etc
data only available on request
no public benchmarking code
why CPLEX and not something better?

Models

In ASCI-printable Symplex notation.

Model 1 as in paper

$N$ and $K$ have correct domains

input N in {1, 2, ...}
input K in {2, 3, ...}
input D[0..(N + 1)][0..(N + 1)] in [0, infinity[
input T[0..(N + 1)][0..(N + 1)] in [0, infinity[
input C1 in {0, 1, ...}
input C2 in {0, 1, ...}
input p[i in 1..N] in {0, 1, ...}
input dp[i in 1..N] in {0, 1, ...}
input MaxCap[PatientTypes] in {0, 1, ...}
input NumPatientsAt[Location][PatientTypes] in {0, 1, ...}
input I[1..N] in {0, 1, ...}
input a[1..N] in [0, infinity[
input tr[1..N] in [0, infinity[
input ta[1..N] in [0, infinity[
input A0 in [0, infinity[
input B0 in [0, infinity[

decision s[1..K] in [0, infinity[
decision e[1..K] in [0, infinity[
decision pn[1..(N + 1)] in {0, 1, ...}
decision dpn[1..(N + 1)] in {0, 1, ...}
decision w[1..N] in [0, infinity[
decision m[1..N] in [0, infinity[
decision x[0..(N + 1)][0..(N + 1)][1..K] in {0, 1}

constraint EQ_1[j in 1..N] is
  sum[i in 0..(N + 1)][k in 1..K] x[i][j][k] = 1
constraint EQ_2[i in 1..N][k in 1..K] is
  sum[j in 0..(N + 1)] x[i][j][k] = sum[j in 0..(N + 1)] x[j][i][k]
constraint EQ_3 is
  sum[i in 1..N] x[0][i][1] = 1
constraint EQ_4[i in 1..(N + 1)] is
  x[i][0][1] = 0
constraint EQ_5 is
  sum[i in 1..N] x[i][N + 1][1] = 1
constraint EQ_6[i in 0..N] is
  x[N + 1][i][1] = 0
constraint EQ_7[k in 2..K] is
  sum[i in 1..(N + 1)] x[0][i][k] = 0
constraint EQ_8[k in 1..(K - 1)] is
  sum[i in 1..(N + 1)] x[0][i][k] = 0

constraint EQ_9[k in 2..(K - 1)] is
  sum[i in 0..N] x[N + 1][i][k] = sum[i in 0..N] x[i][N + 1][k]
constraint EQ_10[k in 2..K] is
  sum[i in 0..N] x[N + 1][i][k] <= 1>
constraint EQ_11 is
  x[N + 1][0][K] = 1
constraint EQ_12 is
  pn[0] = 0
constraint EQ_13 is
  pn[N + 1] = 0
constraint EQ_14 is
  dpn[0] = 0
constraint EQ_15 is
  dpn[N + 1] = 0
constraint EQ_16[i in 0..(N + 1)][j in 1..N][k in 1..K] is
  pn[i] + p[j] + M * (1 - x[i][j][k]) <= pn[j]
constraint EQ_17[i in 0..(N + 1)][j in 1..N][k in 1..K] is
  dpn[i] + dp[j] + M * (1 - x[i][j][k]) <= dpn[j]
constraint EQ_18[i in 1..N] is
  dp[i] <= C1
constraint EQ_19[i in 1..N] is
  dpn[i] <= C2

constraint EQ_20[j in 1..N] is
  s[1] + T[0][j] + M * (1 - x[0][j][1]) <= m[j]
constraint EQ_21[j in 1..N][k in 2..(K - 1)] is
  s[k] + T[N + 1][j] + M * (1 - x[N + 1][j][k]) <= m[j]
constraint EQ_22[i,j in 1..N][k in 1..(K - 1)] is
  m[i] + T[i][j] + M * (1 - x[i][j][k]) <= m[j]
constraint EQ_23[i in 1..N][k in 1..(K - 1)] is
  m[i] + T[i][N + 1] + M * (1 - x[i][N + 1][k]) <= e[k]
constraint EQ_24 is
  s[K] + T[N + 1][0] = e[K]
constraint EQ_25[i in 1..N] is
  tr[i] <= m[i]
constraint EQ_26 is
  A0 <= s[1]
constraint EQ_27 is
  e[K] <= B0
constraint EQ_28[k in 1..K] is
  s[k] <= e[k]
constraint EQ_29[k in 1..(K - 1)] is
  e[k] <= s[k + 1]
constraint EQ_30[i in 1..N] is
  a[i] + T[i][N + 1] - M * (1 - sum[k in 1..(K - 1)] x[i][N + 1][k]) <= w[i]

constraint EQ_31[i, j in 1..N] is
  m[j] - m[i] + w[j] - M * (1 - sum[k in 1..(K - 1)] x[i][j][k]) <= w[i]
constraint EQ_32[i in 1..N] is
  w[i] <= I[i]
constraint EQ_33[i in 1..N] is
  m[i] + w[i] <= ta[i]

minimize sum[i,j in 0..(N + 1)][k in 1..K] x[i][j][k] * D[i][j]

Model 1 rewritten

introduces a patient type to unify elements
eliminates last round … this adapts the number of rounds etc
data checks
makes locations abstract

input PatientLocations finitesubsetof Symbols
local Depot : Symbol
check Depot not in PatientLocations
local Hospital : Symbol
check Hospital not in PatientLocations
local AllLocations := PatientLocations union {Depot, Hospital}

input Distance[AllLocations][AllLocations] in [0, infinity[

input NumRounds in {1, 2, ...}
local Rounds := {1, 2, ..., NumRounds} 

input DepotOpens in [0, infinity[
input DepotCloses in [0, infinity[
check DepotOpens <= DepotCloses

local PatientTypes := {Normal, Disabled}

input TravelTime[PatientLocations][PatientLocations] in [0, infinity[
input MaxTravelTime[PatientLocations] in [0, DepotCloses - DepotOpens]
input BoardingDelay[PatientLocations] in [0, DepotCloses - DepotOpens]
input EarliestBoarding[PatientLocations] in [DepotOpens, DepotCloses]
input LatestArrival[PatientLocations] in [DepotOpens, DepotCloses]
check forall[l in PatientLocations] EarliestBoarding[l] <= LatestArrival[l]

input ShuttleCap[PatientTypes] in {0, 1, ...}
input NumPatientsAt[PatientLocations][PatientTypes] in {0, 1, ...}

decision RoundStart[Rounds] in [DepotOpens, DepotCloses]
decision RoundEnd[Rounds] in [DepotOpens, DepotCloses]
decision NumInShuttle[AllLocations][PatientTypes] in {0, 1, ...}
decision TravelTime[PatientLocations] in [0, DepotCloses - DepotOpens]
decision ArrivalTime[PatientLocations] in [DepotOpens, DepotCloses]
decision Runs[AllLocations][AllLocations][Rounds] in {0, 1}

constraint EQ_1[l in PatientLocations] is
  sum[ll in AllLocations][r in Rounds] Runs[ll][l][r] = 1
constraint EQ_2[l in PatientLocations][r in Rounds] is
  sum[ll in AllLocations] Runs[ll][l][r] = sum[ll in AllLocations] Runs[l][ll][r]
constraint EQ_3 is
  sum[l in PatientLocations] Runs[Depot][l][1] = 1
constraint EQ_4[l in PatientLocations union {Hospital}] is
  Runs[l][Depot][1] = 0
constraint EQ_5 is
  sum[l in PatientLocations] Runs[l][Hospital][1] = 1
constraint EQ_6[l in PatientLocations union {Depot}] is
  Runs[Hospital][l][1] = 0
constraint NeverLeavesDepotAgain[k in Rounds \ {1}] is
  Runs[Depot][Hospital][r] + sum[l in PatienLocations] Runs[Depot][l][r] = 0

constraint EQ_9[k in 2..(K - 1)] is
  sum[i in 0..N] x[N + 1][i][k] = sum[i in 0..N] x[i][N + 1][k]
constraint EQ_10[k in 2..K] is
  sum[i in 0..N] x[N + 1][i][k] <= 1>
constraint EQ_11 is
  x[N + 1][0][K] = 1
constraint EQ_12[pt in PatientTypes] is
  NumInShuttle[Depot][pt] = 0
constraint EQ_12[pt in PatientTypes] is
  NumInShuttle[Hospital][pt] = 0
constraint EQ_16[i in 0..(N + 1)][j in 1..N][k in 1..K] is
  pn[i] + p[j] + M * (1 - x[i][j][k]) <= pn[j]
constraint EQ_17[i in 0..(N + 1)][j in 1..N][k in 1..K] is
  dpn[i] + dp[j] + M * (1 - x[i][j][k]) <= dpn[j]
constraint EQ_18[i in 1..N] is
  dp[i] <= C1
constraint EQ_19[i in 1..N] is
  dpn[i] <= C2

constraint EQ_20[j in 1..N] is
  s[1] + T[0][j] + M * (1 - x[0][j][1]) <= m[j]
constraint EQ_21[j in 1..N][k in 2..(K - 1)] is
  s[k] + T[N + 1][j] + M * (1 - x[N + 1][j][k]) <= m[j]
constraint EQ_22[i,j in 1..N][k in 1..(K - 1)] is
  m[i] + T[i][j] + M * (1 - x[i][j][k]) <= m[j]
constraint EQ_23[i in 1..N][k in 1..(K - 1)] is
  m[i] + T[i][N + 1] + M * (1 - x[i][N + 1][k]) <= e[k]
constraint EQ_24 is
  s[K] + T[N + 1][0] = e[K]
constraint EQ_25[i in 1..N] is
  tr[i] <= m[i]
constraint EQ_26 is
  A0 <= s[1]
constraint EQ_27 is
  e[K] <= B0
constraint EQ_28[k in 1..K] is
  s[k] <= e[k]
constraint EQ_29[k in 1..(K - 1)] is
  e[k] <= s[k + 1]
constraint EQ_30[i in 1..N] is
  a[i] + T[i][N + 1] - M * (1 - sum[k in 1..(K - 1)] x[i][N + 1][k]) <= w[i]

constraint EQ_31[i, j in 1..N] is
  m[j] - m[i] + w[j] - M * (1 - sum[k in 1..(K - 1)] x[i][j][k]) <= w[i]
constraint EQ_32[i in 1..N] is
  w[i] <= I[i]
constraint EQ_33[i in 1..N] is
  m[i] + w[i] <= ta[i]

constraint EQ_34[] is

minimize sum[l1 in AllLocations][l2 in AllLocations][r in Rounds] Runs[l1][l2]][r] * Distance[l1][l2]

Extension 1S as in paper

NYW

Model 2 as in paper

NYW

Interlude: travelling salesman

Generic form of (35) - (40)

input V finitesubsetof Symbols
local N := |V|
check N >= 1
input Distance[V][V] in ]-inf, inf[
check forall[v in V]: Distance[v][v] = 0
// start node ... implies N >= 1
decision Ordering[V] in {1, ..., N}
decision TakeEdge[V][V] in {0, 1}

constraint EntersOnce[v in V] is
  sum[w in V] TakeEdge[w][v] = 1
constraint LeavesOnce[v in V] is 
  sum[w in V] TakeEdge[v][w]] = 1

constraint OrderStartFirst is
  Ordering[Depot] = 1
// This one is optional, if OrderByTakenEdge ranges over V^2
constraint OrderOthers[v in V setminus {s}] is
  Ordering[v] >= 2
// This eliminates subtours
constraint OrderByTakenEdge[v, w in V setminus {s}] is
  Ordering[v] + 1 <= Ordering[w] + N * (1 - TakeEdge[v][w])

Model 2 rewritten

joined patient types
nodes have semantic types
adaption of input data from model 1 (not spelled out in paper)
corrected subtour elimination

input PatientLocations finitesubsetof Symbols
local NumPatientLocations := |PatientLocations|
local Depot : Symbol
check Depot not in PatientLocations
local Hospital : Symbol
check Hospital not in PatientLocations
local AllLocations := PatientLocations union {Depot, Hospital}

input Distance[AllLocations][AllLocations] in [0, infinity[

input DepotOpens in [0, infinity[
input DepotCloses in [0, infinity[
check DepotOpens <= DepotCloses

local PatientTypes := {Normal, Disabled}

input TravelTime[AllLocations][AllLocations] in [0, infinity[
input MaxTravelTime[PatientLocations] in [0, DepotCloses - DepotOpens]
input BoardingDelay[PatientLocations] in [0, DepotCloses - DepotOpens]
input EarliestBoarding[PatientLocations] in [DepotOpens, DepotCloses]
input LatestArrival[PatientLocations] in [DepotOpens, DepotCloses]
check forall[l in PatientLocations] EarliestBoarding[l] <= LatestArrival[l]

input ShuttleCap[PatientTypes] in {0, 1, ...}
input NumPatientsAt[PatientLocations][pt in PatientTypes] in {0, 1, ..., ShuttleCap[pt]}

local PatientNodes := PatientLocations x {Pickup, Dropoff}
local Nodes := PatientNodes union {Depot}

local NodeTravelTime[n, m in Nodes] :=
  | (Depot, Depot) -> 0
  | (Depot, (l, PickUp)) -> TravelTime[Depot][l]
  | (Depot, (_, Dropoff)) -> TravelTime[Depot][Hospital]
  | ((l, Pickup), (ll, Pickup)) -> TravelTime[l][ll]
  | ((l, Pickup), (_, Dropoff)) -> TravelTime[l][Hospital]
  | ((l, Pickup), Depot) -> TravelTime[l][Depot]
  | ((_, Dropoff), (ll, Pickup)) -> TravelTime[Hospital][ll]
  | ((_, Dropoff), (_, Dropoff)) -> 0
  | ((_, Dropoff), Depot) -> TravelTime[Hospital][Depot]

local NodeDistance[n, m in Nodes] :=
  | (Depot, Depot) -> 0
  | (Depot, (l, PickUp)) -> NodeDistance[Depot][l]
  | (Depot, (_, Dropoff)) -> NodeDistance[Depot][Hospital]
  | ((l, Pickup), (ll, Pickup)) -> NodeDistance[l][ll]
  | ((l, Pickup), (_, Dropoff)) -> NodeDistance[l][Hospital]
  | ((l, Pickup), Depot) -> NodeDistance[l][Depot]
  | ((_, Dropoff), (ll, Pickup)) -> NodeDistance[Hospital][ll]
  | ((_, Dropoff), (_, Dropoff)) -> 0
  | ((_, Dropoff), Depot) -> NodeDistance[Hospital][Depot]

decision NumInShuttle[Nodes][PatientTypes] in {0, 1, ...}
decision ArrivalTime[PatientLocations] in [DepotOpens, DepotCloses]
decision Ordering[Nodes] in {0, 1, ..., 2 * NumPatientLocations}
decision TakeEdge[Nodes][Nodes] in {0, 1}

// (35)
constraint EntersOnce[n in Nodes] is
  sum[m in Nodes] TakeEdge[m][n] = 1
// (36)
constraint LeavesOnce[n in Nodes] is 
  sum[m in Nodes] TakeEdge[n][m]] = 1

// (37)
constraint OrderDepotFirst is
  Ordering[Depot] = 0
// (38)
constraint OrderPatientNode[n in PatientNodes] is
  Ordering[n] >= 1

// (39) is handled by domain of Ordering

// (40)
constraint OrderByTakenEdge[n, m in PatientNodes] is
  Ordering[n] + 1 <= Ordering[m] + 2 * NumPatientLocations * (1 - TakeEdge[n][m])

// (41)
constraint PickupBeforeDropoff[l in PatientLocations] is
  Ordering[(l, Pickup)] + 1 <= Ordering[(l, Dropoff)]
// (42)
constraint PickupDropoffOrderCoupling[l, ll in PatientLocations] is
  TakeEdge[(l, Pickup), (ll, Pickup)] = TakeEdge[(l, Dropoff), (ll, Dropoff)]

// (45) and (47)
constraint NumPatientsAtDepot[pt in PatientTypes] is
  NumInShuttle[Depot][pt] = 0

// (44) and (46)
constraint NumPatientsDroppedOff[l in PatientLocations][pt in PatientTypes] is
  NumInShuttle[(l, Dropoff)][pt] = 0

// (47) and (48)
constraint NumPatientsTransported[l in PatientLocations][n in Nodes][pt in PatientTypes] is
  NumInShuttle[n][pt] + NumPatientsAt[l][pt] - M[l][pt] * (1 - TakeEdge[n][(l, Pickup)]) <= NumInShuttle[(l, Pickup)][pt]
  with
  M[l][pt] := NumPatientsAt[l][pt] + ShuttleCap[pt]

// (51)
constraint TimeDepotToPatientNode[n in PatientNodes] is
  ArrivalTime[Depot] + NodeTravelTime[Depot][n] - M[n] * (1 - TakeEdge[Depot][n]) <= ArrivalTime[n]
  with
  M[n] := (DepotCloses - DepotOpens) + NodeTravelTime[Depot][n]

// (52)
constraint TimePickupToPatientNode[l in PatientLocations][n in PatientNodes] is
  ArrivalTime[(l, Pickup)] + NodeTravelTime[(l, Pickup)][n] + BoardingDelay[l] - M * (1- TakeEdge[(l, Pickup)][n]) <= ArrivalTime[n]
  with
  M[l][n] := (DepotCloses - DepotOpens) + NodeTravelTime[(l, Pickup)][n] + BoardingDelay[l]

// (53)
constraint TimeHospitalToPatientNode[l in PatientLocations][n in PatientNodes] is
  ArrivalTime[h] + NodeTravelTime[h][n] - M * (1- TakeEdge[h][n]) <= ArrivalTime[n]
  with
  h := (l, Dropoff)
  M[l][n] := (DepotCloses - DepotOpens) + NodeTravelTime[h][n]

// (54)
constraint DoNotTooArriveEarly[l in PatientLocations] is
  EarliestBoarding[l] <= ArrivalTime[(l, Pickup)]

// (55)
constraint TravelTimeBound[l in PatientLocations] is
  ArrivalTime[(l, DropOff)] - ArrivalTime[(l, Pickup)] <= MaxTravelTime[l]

// (56)
constraint ArriveAtHospitalInTime[l in PatientLocations] is
  ArrivalTime[(l, DropOff)] <= LatestArrival[l]

// (57) is handled by bounds in domain of ArrivalTime

// (58)
constraint FinishAtDepotBeforeClosing[l in PatientLocations] is
  ArrivalTime[(l, DropOff)] + TravelTime[Hospital][Depot] <= DepotCloses

minimize sum[n, m in Nodes] NodeDistance[n][m] * TakeEdge[n][m]

Extension 2S as in paper

NYW