Abstract Mixed‐integer linear programming is widely used in hydropower reservoirs operation optimisation because of its modelling flexibility, solution stability, and global search capability. A major obstacle to applying mixed‐integer linear programming to hydropower reservoirs operation optimisation, especially for long‐term operations of large‐scale hydropower systems, is the non‐linear phenomenon while converting water energy into electricity, which is normally expressed as a non‐linear bivariate function, that is, hydropower production function. To cross this technological barrier, a novel linearisation method for hydropower production function is proposed. The method uses rectangular meshing techniques to approximate hydropower production function. Then Special Ordered Sets of type 2 (SOS2)constraints are adopted for sub‐rectangle selection, which are formulated based on the binary branching schemes generated by binary reflected grey code. The proposed method can linearise the hydropower production function with only logarithmic‐sized binary variables of the traditional methods and maintain high accuracy. The method is then applied to optimising the operation of the Lancang River Basin hydropower system, which has 13 head‐dependent reservoirs with a deterministic model and a two‐stage stochastic model. The results show that the authors' method significantly outperforms conventional linearisation methods and can obtain accurate solutions within an acceptable time for deterministic and stochastic large‐scale hydropower reservoirs operation optimisation problems.