The relay optimization expresses quite a challenge for smooth and optimal operation of power system networks. The relay optimization is formulated as a mixed integer non-linear problem and is highly constrained. Furthermore, a reliable relaying system must be able to detect and isolate the faulted portion in a timely manner. Therefore, it is necessary to find optimal parameters for relay settings to be able to respond in a timely way to the encountered fault and at the same time keep in consideration the operational and coordination constraints. This paper proposes modified Harris hawk optimization (MHHO), which is based on the intelligent preying tactics of Harris hawks and the improvement of intended modifications, crowding distance and roulette wheel selection. The proposed algorithm has been tested on IEEE 8 and 15-bus systems, using MATLAB programming. The test systems are the distribution networks covering the medium level voltage for consideration. The simulation results verified the success of MHHO to find optimal settings for the relays. For IEEE 8-bus system, MHHO was able to give 35.45% improvement in the results in comparison to other algorithms. Furthermore, for the IEEE 15-bus system, MHHO showed 24.09% improvement on average. The comparison of the results obtained by MHHO with the other state-of-the-art algorithms proved that it is the strong candidate for optimization of the relay coordination problem.