PARAFAC2 is a useful algorithm for decomposing tensors that do not have low-rank variation such as e.g. PARAFAC requires. It has been applied in analyzing different types of multi-way data, such as GC-MS data. Since the optimization in fitting the loss function of PARAFAC2 is non-convex, the PARAFAC2 model suffers from local minima. In this paper, we investigated the local minima issue in the PARAFAC2 decomposition. We observed that the decomposed loading matrices corresponding to local minima and global minima are in general very different even though the loss function values can be very similar. We also observed that overfitting always led to higher fraction of local minima. For coping with local minima in PARAFAC2, we investigated the effect of non-negativity constraints on avoiding local minima, proposed some PARAFAC2 algorithms with optimized initializations, and implemented a simple local minima avoidance procedure in the general PARAFAC2 algorithm. These remedies are evaluated and illustrated by the use of different types of datasets such as simulated data, GC-MS data and EEG data. It is finally concluded that imposing non-negativity constraints on the PARAFAC2 model decreased the fraction of local minima significantly, using the proposed optimized initialization PARAFAC2 algorithm decreased the average fraction of local minima from the highest (84%) to the lowest (0%), and using the proposed local minima avoidance PARAFAC2 algorithm eliminated all the local minima in our case. In case of modeling data with high complexity, it will be promising to use all the remedies at the same time in order to avoid local minima PARAFAC2 models successfully. To the best of our knowledge, this is the first paper that investigate the local minima issue in the context of PARAFAC2 for irregular tensor.