Understanding the transition of flows from steady to chaotic regimes is crucial for many engineering applications. Mathematically, these transitions are characterized by the so-called bifurcating phenomena, through which the system exhibits different qualitative properties. Thus, a comprehensive bifurcation analysis is often fundamental to provide valuable insights on the system behavior, but such in-depth description is usually unfeasible due to the high computational costs needed to obtain accurate enough discretizations. Toward this goal, Reduced Order Models (ROMs) represent a computationally cheap class of strategies for reducing the inherent dimensionality and efficiently investigating bifurcating phenomena for several configurations of the system. As a benchmark study case, we consider the steady and unsteady fluidic pinball problem that exhibits a particularly rich dynamical behavior, and consists of three cylinders with centers placed on the vertices of an equilateral triangle pointing upstream [1]. We investigate flow regimes ranging from steady to chaotic states as the system undergoes multiple branching phenomena, including Hopf and Pitchfork bifurcations, even at moderate Reynolds numbers. Moreover, varying the Reynolds number, in conjunction with introducing a parametrized symmetric forcing through rotating the two rear cylinders [2], we identify diverse flow regimes and track the evolution of the bifurcation points. Detection of bifurcations is accomplished through linear stability analysis of the steady state solution at different parameters, tracking the eigenvalues as they cross the y-axis. To create a low dimensional model that can accurately describe this complex behavior and reconstruct the bifurcations diagrams, we exploit a ROM for the steady Navier-Stokes equations using a Galerkin projection on the modes obtained via Proper Orthogonal Decomposition (POD). To enhance accuracy and efficiency, especially in the more complex multi-parameter framework, we also develop a local ROM approach to perform branch-wise reduction [3, 4]. In particular, during the offline stage we cluster the solution space via an enhanced k-means algorithm, establishing distinct local ROMs, while the selection of the local model in the online phase is done according to the Euclidean distance of the initial state in the solution space.