In this paper, multiple virtual synchronous machines (VISMAs) with fuzzy proportional integral derivative (FPID) controllers optimized by differential evolution (DE) are proposed to maintain frequency stability in the grid in the presence of renewable penetration, such as wind and solar photovoltaic (PV) systems, residential loads, and industrial loads, by reducing the area control error in the objective function. Simulations are conducted using MATLAB/Simulink, and in the optimization process, the integral of the time-weighted absolute error (ITAE) is used as the objective function. In the work to obtain optimized values of renewable energy sources (RESs), fuzzy membership functions, controller gain parameters, and loads for system modeling, differential evolution and genetic algorithm (GA) methods are applied and the results were compared. It was shown that better results were achieved while FPID controllers were optimized by DE in the presence of multiple VISMAs than DE in the presence of single VISMAs and GA in multiple VISMAs. Moreover, the study is compared to integral control methods in which, compared to all controllers, the proposed controller reduces undershoot by 0.0674 Hz more than a single VISMAs, in which it is improved approximately by 97.82%. Similarly, the proposed controller improves the system settling time, rise time, and overshoot by more than 99.5% compared to the classical integral controller. To examine the robust operation of the system under the proposed controller, the system was run under a wide range of disturbances and uncertainties using random load perturbation of ± 20%, in which the proposed controller retains the system frequency by reducing or damping the system oscillation.