The project on your CV,
reverse engineered.

That last sentence is the hook. It invites exactly the follow-up you most want, and the answer to it is the section on the MILP below. Deliver the pitch, then stop talking and let them choose where to dig.

# src/processing/activity.py
result = df.sort_values(["vessel_id", "timestamp"]).copy()
time_diff = result.groupby("vessel_id")["timestamp"].diff()
result["activity_hours"] = time_diff.dt.total_seconds() / 3600
result["activity_hours"] = result["activity_hours"].fillna(0)

Three lines. Every modelling assumption in the project that is not stated anywhere lives in them.

conditions = [
    in_anchorage & (speed < 1),
    in_port & (speed > 1),
    ~in_port & (speed >= 1),
]
choices = ["Anchorage", "Maneuver", "Transit"]

result["operating_mode"] = pd.Series("Drifting", index=result.index)
for condition, choice in zip(conditions, choices, strict=False):
    result.loc[condition, "operating_mode"] = choice

# src/processing/load_factor.py
result["max_speed"] = MAX_SPEED_MULTIPLIER * result["vref"]   # 1.066 x reference speed
result["load_factor"] = (result["speed_knots"] / result["max_speed"]) ** 3
result["load_factor"] = result["load_factor"].round(2)

# under way but barely moving still burns something
transit_or_maneuver = result["operating_mode"].isin(["Transit", "Maneuver"])
below_floor = result["load_factor"] < LOAD_FACTOR_FLOOR   # 0.02
result.loc[transit_or_maneuver & below_floor, "load_factor"] = LOAD_FACTOR_FLOOR

The .round(2) is not cosmetic. Downstream, the low-load adjustment table is indexed by integer percent, and rounding to two decimals here is what guarantees load_factor * 100 lands on a key that exists. It is an implicit contract between two modules that nothing in the code states, which is the kind of coupling worth admitting to.

Drag the speed. Watch how little of the range does most of the burning, and where the 2 percent floor takes over.

# src/emissions/calculator.py
in_range = (lf_percent >= 2) & (lf_percent <= 20)
lf_clamped = lf_percent.where(in_range)          # out of range becomes NaN

llaf_co2 = lf_clamped.map(llaf_lookup["CO2"]).fillna(1.0)   # NaN maps to NaN, then to 1.0

The same missing fuel type fails two different ways

Both modules look a fuel type up in a dictionary built from the same spreadsheet sheet. They disagree about what happens when the fuel type is not there.

The second is the dangerous one. A silent zero in a cost term does not crash, does not warn, and biases the optimum toward exactly the vessels whose data you were missing. The fix is not to pick one default, it is to validate the join up front: assert that every fuel type appearing in the movements data has a row in the factors workbook, and fail loudly at load time if it does not.

You pay P today for a ship you will use for n years, and you want the equivalent level annual payment A. The present value of that annuity has to equal P:

P = A * sum(1 / (1+r)^t for t in 1..n)
  = A * (1 - (1+r)^-n) / r

# invert for A, multiply top and bottom by (1+r)^n
A = P * r(1+r)^n / ((1+r)^n - 1)
              ^^^^^^^^^^^^^^^^^^^^^^^^ this is the capital recovery factor

# src/cost/calculator.py
def _calculate_crf() -> float:
    r = DISCOUNT_RATE          # 0.08
    n = SHIP_LIFETIME_YEARS    # 30
    factor = (1 + r) ** n
    return r * factor / (factor - 1)

s = SALVAGE_RATE * p                        # you get 10% back at year 30
annual = (p - s) * crf + DISCOUNT_RATE * s  # amortise the rest, plus carry on the salvage
monthly = annual / 12

# decision variable: x_i in {0, 1}, one per vessel

minimise    sum_i  c_i * x_i                        # adjusted monthly cost

subject to  sum_i  dwt_i * x_i  >=  4,576,667       # the cargo has to fit
            sum_i  (s_i - T) * x_i  >=  0           # average safety >= T  (see below)
            sum_{i in type f}  x_i  >=  1   for each fuel type f

# what you want to say, which is not linear:
sum_i s_i * x_i
---------------  >=  T
  sum_i x_i

# the denominator is strictly positive whenever any vessel is chosen,
# so multiplying both sides by it preserves the direction of the inequality:
sum_i s_i * x_i  >=  T * sum_i x_i

# collect, and it is linear in x:
sum_i (s_i - T) * x_i  >=  0

src/optimization/solver.py, the constraint named min_avg_safety.

Read the final form as an accounting identity and it stops feeling like a trick. Each vessel carries a signed surplus s_i - T. Vessels above the threshold bank surplus, vessels below it spend surplus, and the fleet is admissible exactly when the books balance. That reading is also what tells you how to extend it: any constraint on a weighted mean linearises the same way, and a constraint on a ratio of two general sums is the Charnes-Cooper transformation, which is the same move applied to the whole program rather than one row.

Toggle vessels into the fleet. Watch the true average and the linear surplus cross zero at exactly the same moment, which is the whole claim.

prob = pulp.LpProblem("Fleet_Selection", pulp.LpMinimize)
x = {vid: pulp.LpVariable(f"x_{vid}", cat=pulp.LpBinary) for vid in vessel_ids}

prob += pulp.lpSum(cost_lookup[vid] * x[vid] for vid in vessel_ids)

prob += (pulp.lpSum(dwt_lookup[vid] * x[vid] for vid in vessel_ids) >= params.min_dwt, "min_dwt")

prob += (pulp.lpSum((safety_lookup[vid] - params.min_avg_safety) * x[vid]
                    for vid in vessel_ids) >= 0, "min_avg_safety")

Three things about the PuLP idiom are worth being able to explain. The first += with a bare expression sets the objective; every later one with an inequality adds a constraint. Naming each constraint is not decoration, it is what lets you read the LP file and inspect duals afterwards. And "each vessel at most once" is not written anywhere because it is already implied by the variable being binary, which is the kind of thing an interviewer will check you noticed.

The method itself is textbook and correctly applied. To trace a cost-against-safety frontier you do not put both objectives in one function with weights, you fix one as a constraint at a series of levels and minimise the other. Here the safety threshold sweeps 3.0 to 5.0 in steps of 0.1, twenty-one solves, each producing a point.

Weighted sums would have been the wrong choice, and knowing why is a good thing to have ready: a weighted sum can only ever find points on the convex hull of the frontier, so with integer decisions it silently skips every solution sitting in a non-convex dent. The epsilon-constraint method finds those.

The Shapley value of a player is its average marginal contribution across all orderings in which the coalition could have been assembled. With a fleet of any size the number of orderings is factorial, so the code samples: one thousand random permutations, add vessels one at a time, record the change in coalition value at each insertion. Sampling for Shapley is the standard approach and it is implemented correctly.

The problem is the value function. An infeasible coalition is assigned a penalty of twice the total cost of every vessel in the dataset, and a feasible one is assigned its actual cost. Contribution is recorded as previous - current.

# every permutation looks like this:

step        coalition state        value                contribution recorded
----        ---------------        -----                ---------------------
1..k-1      infeasible             PENALTY              PENALTY - PENALTY = 0
k           becomes feasible       real cost            PENALTY - cost   = enormous, positive
k+1..n      feasible, growing      real cost + more     small, NEGATIVE

The intent is good and the structure is right. Start from the optimal fleet, propose flipping one vessel in or out, accept downhill moves always and uphill moves with probability exp(-beta * delta), and count how often each vessel appears. Vessels present in nearly every near-optimal fleet are structurally essential; vessels that come and go are substitutable. That is a Metropolis sampler and it is a genuinely nice way to ask a robustness question that a single optimum cannot answer.

beta = 0.0001                                    # src/optimization/mcmc.py, default
cost_diff = proposed_cost - current_cost      # costs are monthly USD, order 1e6
acceptance_prob = 1.0 if cost_diff <= 0 else math.exp(-beta * cost_diff)

Put the numbers in. A modest uphill move of one hundred thousand dollars gives exp(-0.0001 * 100000) = exp(-10), about four in one hundred thousand. A million-dollar move gives exp(-100), which is zero for every practical purpose.

The default beta is marked on the slider. Move the cost difference into the range this problem actually produces and watch the acceptance probability disappear.

The five-minute refresh