Inclusion-driven structure learning of Bayesian networks, or idlBNs, converges to the generative structure as the sample size grows large and as long as that structure is an acyclic digraph (DAG) over the observed random variables. Because Markov equivalence of Bayesian networks organizes the search space of DAGs in equivalence classes, an obvious choice for such an approach is the greedy equivalence search (GES) algorithm, which carefully traverses the space of essential graphs, the canonical elements of those equivalence classes, following an inclusion path. GES is adapted to data produced by multiple intervention experiments in the greedy interventional equivalence search (GIES) algorithm. The algorithmic complexity of both GES and GIES is in the worst case exponential in the number of vertices, but it can be reduced to polynomial by bounding the vertex degree during the search, albeit at the cost of losing the large-sample optimality guarantee. Inclusion-driven structure learning can also be implemented in the search space of DAGs, as in the hill-climber Monte Carlo (HCMC) algorithm, whose stochastic nature confers the advantage of a polynomial-time bounded algorithmic complexity. Here, we introduce the interventional HCMC (iHCMC) algorithm, an inclusion-driven structure learning algorithm for interventional data in DAG-space. Using synthetic Gaussian data, we verify that iHCMC preserves the large-sample optimality for interventional data with polynomial-time complexity independent of the sparsity of the generative structure.