First choice :
Convert the Laplace transform of your process to the Z-domain using the ZOH method as it models your DAC. In your case, your DAC is a PWM.
Second choice :
Work in the Laplace domain but take in account the delay due to your PWM, which should be T/2 where T is your PWM period. To model the delay, you can use the Padé approximation.
I personnaly used the second technique. It was easier, more intuitive. But both techniques will work.
For the numerical issues : The order of your transfer function might get large when including the controllers, I used minreal too, but you can convert your tranfer functions to ZPK, this solves many of the numerical issues you might have.
Edit : As Tim Wescott mentioned in the comments. The ZOH modelization is accurate when using a real DAC. Not so for a PWM. It's better off to work in the s-domain and try to model the delays (sampling, PWM, computation, etc.) as accurately as possible.