Extrusion bioprinting is a popular biofabrication technique for creating viable biological constructs with applications ranging from regenerative medicine and cancer research to therapeutics screening. A primary challenge within extrusion bioprinting lies in preserving cell viability and function amidst the shear forces generated during extrusion. Therefore, the formulation of shear-thinning bioinks with tailored flow properties is essential to minimize shear stress accumulation during printing. However, this task is complex, as the bioink’s flow properties are contingent upon numerous factors, including but not limited to the polymer concentration, spatial cell density in the bioink, molecular weight of the polymer, and the nozzle’s geometry. The current methods for adjusting flow properties are heuristic, posing challenges to the rapid development models and automation of the bioprinting workflow. In this study, we introduce a Gaussian Process-based multi-objective Bayesian optimization framework that rapidly identifies Pareto-optimal flow property values to minimize hydrodynamic shear stresses during extrusion. The optimization focuses on four flow parameters, minimum viscosity, maximum viscosity, consistency index, and flow behavior index from the power-law model of shear-thinning fluids alongside the nozzle outlet diameter. Following an initial set of 24 experimental runs, our sequential design approach required 21 iterations to achieve convergence on a set of flow parameters that meet the Pareto optimality criteria. The presented strategy, which integrates physics-based numerical models with data-driven optimization models, emerges as a tool for the rapid fine-tuning of bioink flow behaviors. The multi-objective framework also sets the stage for future exploration into cell viability and functional outcomes post-bioprinting.