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.