This paper concerns the control of a time fractional diffusion system defined in the Riemann–Liouville sense. It is assumed that the system is subject to hysteresis nonlinearity at its input, where the hysteresis is mathematically modeled with the Duhem operator. To compensate the effects of hysteresis nonlinearity, a fractional order Proportional+Integral+Derivative (PID) controller is designed by minimizing integral square error. For numerical computation, the Riemann–Liouville fractional derivative is approximated by the Grünwald–Letnikov approach. A set of algebraic equations arises from this approximation, which can be solved numerically. Performance of the fractional order controllers are analyzed in comparison with integer order controllers by simulation results, and it is shown that the fractional order controllers are more advantageous than the integer ones.