  The time-evolution or equations of motions for many systems are usually
described by a set of first-order ordinary differential equations, and for a
variety of physical observables, the long-time limit or steady state solution
is desired. In the case of open quantum systems, the time-evolution of the
reduced density matrix is described by the Liouvillian. For inverse design or
optimal control of such systems, the common approaches are based on brute-force
search strategies. Here, we present a novel methodology, based on automatic
differentiation, capable of differentiating the steady state solution with
respect to any parameter of the Liouvillian. Our approach has a low memory
cost, and is agnostic to the exact algorithm for computing the steady state. We
illustrate the advantage of this method by inverse designing the parameters of
a quantum heat transfer device that maximizes the heat current and the
rectification coefficient. We also optimize the parameters of various Lindblad
operators used in the simulation of energy transfer under natural incoherent
light.
