Rehabilitation (Rehab) exercise can benefit cardiac patients as it can promote the recovery and improve the heart wellness. However, heart failure (HF) patients can only take mild exercise, since excessive exercise may lead to fatal events. It is important to control the exercise intensity at a desired level to maximize exercise benefit. Heart Rate (HR) is an essential factor for measuring exercise intensity. Mathematical models of HR can be used to study exercise physiology. However, HR models involve model uncertainty, resulting from model calibration or variability in patients. It is important to quantify the effect of uncertainty on HR prediction for optimizing exercise intensity, such as treadmill speed. A probabilistic model-based control design is presented in this work to obtain an optimal treadmill speed for Rehab exercise in the presence of uncertainty. To obtain a computationally tractable formulation, the generalized polynomial chaos (gPC) theory is used to propagate uncertainty via a model to HR predictions, and predict slow-acting responses such as peripheral local metabolism that can be used to evaluate exercise outcome for individual patients. The speed control of treadmill is formulated as an optimization problem that can maximize the exercise outcome, while minimizing the slow-acting effects. The effectiveness of the proposed control design was experimentally verified with simulations, showing potentials in the exercise control of individual patients.