The simulation of nonlinear optical experiments in the condensed phase often requires including relaxation between different states of matter due to interaction with an external bath. This is commonly done using second-order Redfield or Lindblad equations. Here, we derive closed expressions for relaxation equations to any order in the coupling to a harmonic bath. This is done using a compact superoperator notation combined with reduced equations of motions based on Time Convolutionless theory (TCL), partial ordering prescription.