Phase-contrast magnetic resonance imaging (PC-MRI) allows us to acquire biofluid flow velocity maps, whereas MRI data is restricted by spatiotemporal resolution limitations and contains theoretically inevitable errors. Although various approaches to estimating actual velocity from MR velocity maps using the mass and momentum conservation laws have been proposed, practically reasonable methodologies are still not well established. This study investigates a practical strategy for estimating physically consistent velocities from MR velocity maps based on variational optimal boundary control through examples of the 2D steady Stokes flow as an incompressible viscous fluid. We defined a minimization problem of the sum of squared residuals between MR and the estimated velocity at all pixels (voxels) considering the image data structure with respect to the Dirichlet boundary velocity condition subject to flow governing equations based on variational formulations. This optimization problem is treated as an unconstrained optimization problem by deriving the Lagrange functional, including the cost function, regularization term, and constraint conditions. The optimality condition is computed using the adjoint variable method in a finite element manner. The boundary velocity profile is iteratively updated by the optimality condition using gradient-based optimization until convergence. Numerical examples for 2D Poiseuille flow with noise-free and noisy reference data demonstrated good convergencies of the cost function minimization. The estimated flow velocities were in excellent agreement with reference data. Finally, we demonstrated the viability of the velocity estimation using the actual MR velocity of the cerebrospinal fluid flow. The proposed approach with further considerations specialized for the MRI may be feasible in providing physically consistent velocity profiles in a versatile target of the biofluid flow.