Regardless of the huge advances made in the design and fabrication of mid-infrared and terahertz quantum cascade lasers, success in accessing the ~3-4 mm region of the electromagnetic spectrum has remained limited. This fact has brought about the need to exploit resonant intersubband transitions as powerful nonlinear oscillators, consequently enabling the occurrence of large nonlinear optical susceptibilities as a means of reaching desired wavelengths. In this work, we present a computational model developed for the optimization of second-order optical nonlinearities in In0.53Ga0.47As/Al0.48In0.52As quantum cascade laser structures based on the implementation of the Genetic algorithm. The carrier transport and the power output of the structure were calculated by self-consistent solutions to the system of rate equations for carriers and photons. Both stimulated and simultaneous double-photon absorption processes occurring between the second harmonic generation-relevant levels are incorporated into rate equations and the material-dependent effective mass and band non-parabolicity are taken into account, as well. The developed method is quite general and can be applied to any higher order effect which requires the inclusion of the photon density equation. [Projekat Ministarstva nauke Republike Srbije, br. III 45010]