In the present study, we investigated the possibilities and limitations of computer-assisted method development (CAMD) for the HILIC separation optimization of a mixture of 13 isomeric hydroxy- and aminobenzoic acids on a ZIC-HILIC column. The isocratically obtained Neue and Kuss retention parameters enabled the accurate gradient retention modeling for peaks eluting well within the gradient (mean error of 2.7 %). The prediction errors for peaks eluting at the end of the gradient could be reduced from 8.8 to 6.1 % by implementing the isocratic regime after the gradient into the expression for the gradient retention factor. The prediction of the corresponding peak widths improved significantly for certain compounds and gradient profiles using individual gradient N values for each compound compared to employing a single N value for all compounds and gradient profiles. Two gradient optimization strategies (constructing the Rs map based on individual retention modeling and predictive elution stretching and shifting, PEWS2) resulted in a reasonable separation of the challenging mixture of 13 isomeric hydroxy- and aminobenzoic acids on the ZIC-HILIC column. Overall, the optimization was limited by the steep decrease in N (dropping to the isocratic N value) and corresponding increase in peak width when increasing the gradient time. The discrimination factors d0 were used to assess the resolution between peaks varying widely in height. The best separation was found to be obtained via the PEWS2 approach. Both the individual retention modeling and PEWS2 strategies corresponded to a total instrument time less than 12 h (including equilibration). Finally, it was found that the salt concentration had a significant effect on both the retention and the peak shape of the compounds, resulting in a small “solution domain” at 10 mM. Coupled columns with higher efficiencies are suggested to improve the resolution and robustness of the separation.