For a graph $G= (V,E)$, a double Roman dominating function (DRDF) is a function $f : V \to \{0,1,2,3\}$ having the property that if $f (v) = 0$, then vertex $v$ must have at least two neighbors assigned $2$ under $f$ or {at least} one neighbor $u$ with $f (u) = 3$, and if $f (v) = 1$, then vertex $v$ must have at least one neighbor $u$ with $f (u) \ge 2$. In this paper, we consider the double Roman domination problem, which is an optimization problem of finding the DRDF $f$ such that $\sum_{v\in V} f (v)$ is minimum. We propose {five integer linear programming (ILP) formulations and one mixed integer linear programming formulation with polynomial number of constraints for this problem. Some additional valid inequalities and bounds are also proposed for some of these formulations.} Further, we prove that {the first four models indeed solve the double Roman domination problem, and the last two models} are equivalent to the others regardless of the variable relaxation or usage of a smaller number of constraints and variables. Additionally, we use one ILP formulation to give an $H(2(Δ+1))$-approximation algorithm. All proposed formulations and approximation algorithm are evaluated on randomly generated graphs to compare the performance.