I’m revisiting the topic of soft time window and hope this is not too late.
Regarding the calculation/estimation of epsilon here, I have some idea:
First we will need a state updater which is similar to UpdateVehicleDependentPracticalTimeWindows. In this state updater, for each activity j and each vehicle v, we do not store only one latestArrivalTime. Instead, we store an array of latestArrivalTimes, one for each of activity j, all subsequent activities and vehicle v latest end.
For example, assume activity j has one subsequent activity m, then, for activity j and vehicle v, we have an array of three latestArrivalTimes [t1, t2, t3], which is the latest arrival time at j such that we will not be late for each of the activities (i.e., j, m and route end), respectively.
Now, assume we have a penalty function for being late for the soft time window, g(delta_t), then we can construct a penalty function for activity j, vehicle v and arrival time t:
G_j,v(t) = g(t-t1)+g(t-t2)+g(t-t3)
Note that G will be a piecewise function, because g is one.
Then this function G can be used in the constraint to calculate epsilon:
epsilon = G_j,v(new arrival time t') - G_j,v(old arrival time t)
where the value for the old arrival time can be passed via state updater (when newAct is not DeliverShipment).
What do you think?