import numpy as np t_star = 8 + 50/60 tas = np.linspace(7.0, 10.0, 10) tts = np.array([15.0, 17.0, 17.0, 21.0, 23.0, 29.0, 27.0, 23.0, 19.0, 15.0]) if t_star not in tas: t_star_id = np.argmax(tas > t_star) tas = np.insert(tas, t_star_id, t_star) t_star_tt = (tts[t_star_id-1] + tts[t_star_id]) / 2 tts = np.insert(tts, t_star_id, t_star_tt) alpha = 15 beta = 5 gamma = 20 schedule_costs = ( beta * (t_star - tas) * (t_star > tas) + gamma * (tas - t_star) * (tas > t_star) ) costs = alpha * tts / 60 + schedule_costs mu = 1 exp_costs = np.exp(-costs / mu) probs = list() for i in range(len(tas)-1): m = min(exp_costs[i], exp_costs[i+1]) M = max(exp_costs[i], exp_costs[i+1]) prob = (tas[i+1] - tas[i]) * (m + (M-m) / 2) probs.append(prob) probs = np.array(probs) / np.sum(probs) probs = np.insert(probs, 0, 0) cum_probs = np.cumsum(probs)