We present AquaNutriOpt, a user-friendly Python package designed to tackle a complex combinatorial optimization problem aimed at optimizing nutrient management for the control of harmful algal blooms. This optimization process involves the identification of optimal Best Management Practices (BMPs) and Treatment Technologies (TTs). AquaNutriOpt is constructed based on a novel integer programming model, which we present in this paper. The package can accommodate various user inputs, automatically transforming them into an optimization model, and then solving it using a free solver. To demonstrate AquaNutriOpt’s efficacy, we conduct a series of experiments on two watersheds around Lake Okeechobee in Florida, USA. These experiments illustrate that the optimal BMPs/TTs obtained by AquaNutriOpt can significantly reduce Phosphorus loads into the lake across various budget scenarios. We validate the results by running simulations with the process-based Watershed Assessment Model (WAM), confirming that the estimated percentage reductions closely align with the reports from WAM.