Nonlinear optical responses in second harmonic generation (SHG) of van der Waals heterobilayers, Janus MoSSe/MoS2, are theoretically optimized as a function of strain and stacking order by adopting an exchange-correlation hybrid functional and a real-time approach in first-principles calculation. We find that the calculated nonlinear susceptibility, χ(2), in AA stacking (550 pm/V) becomes three times as large as AB stacking (170 pm/V) due to the broken inversion symmetry in the AA stacking. The present theoretical prediction is compared with the observed SHG spectra of Janus MoSSe/MoS2 heterobilayers, in which the peak SHG intensity of AA stacking becomes four times as large as AB stacking. Furthermore, a relatively large, two-dimensional strain (4%) that breaks the C3v point group symmetry of the MoSSe/MoS2, enhances calculated χ(2) values for both AA (900 pm/V) and AB (300 pm/V) stackings 1.6 times as large as that without strain.