The paper presents a novel batch-to-batch control approach for crystallization processes, which can be used for designing the shape of the crystal size distribution (CSD) and the habit of the crystals, to robustly achieve desired product properties. The method is implemented in a hierarchical structure. On the lower level a supersaturation control approach is used that drives the system in the phase diagram according to a concentration versus temperature trajectory, providing the within batch control. On the higher level a robust model-based optimization algorithm adapts the setpoint of the supersaturation controller, after each batch, using the adapted model parameters by applying a batch-to-batch moving horizon estimation (MHE) approach, to counteract the effects of changing operating conditions, and parameter uncertainties. The process is modelled using the two dimensional population balance equation (PBE), which is solved using the method of characteristics (MOC). The control approach is corroborated through simulations studies. The results illustrate the benefits of the robust iterative learning approach which is able to control both the entire CSD and the habit of the product crystals, reducing significantly the variability in the product properties after only four batches.