The Industry 4.0 paradigm has emerged to enhance the efficiency, productivity, and adaptability of industrial processes. Within this context, process optimization plays a crucial role. This work introduces a novel bi-objective nonlinear mathematical programming (NLP) model aimed at optimizing the operation of a physicochemical phosphorus removal system of an industrial wastewater treatment plant (WWTP). The model seeks to simultaneously minimize both the process costs and the amount of phosphorus discharged into water bodies, addressing significant environmental concerns. A simulation model was developed using commercial software to estimate the phosphorus concentration and pH levels of the treated effluent and an approximation function was derived to be used in the NLP model. Computational experiments confirmed the competitiveness of the proposed model which is solved using two different bi-objective methods.