Solve the optimization problem by value iteration. See the
        `solve` method.

def value_iteration(self, v_init=None, epsilon=None, max_iter=None):
        """
        Solve the optimization problem by value iteration. See the
        `solve` method.

        """
        if self.beta == 1:
            raise NotImplementedError(self._error_msg_no_discounting)

        if max_iter is None:
            max_iter = self.max_iter
        if epsilon is None:
            epsilon = self.epsilon

        try:
            tol = epsilon * (1-self.beta) / (2*self.beta)
        except ZeroDivisionError:  # Raised if beta = 0
            tol = np.inf

        v = np.empty(self.num_states)
        if v_init is None:
            self.s_wise_max(self.R, out=v)
        else:
            v[:] = v_init

        # Storage array for self.bellman_operator
        Tv = np.empty(self.num_states)

        num_iter = self.operator_iteration(T=self.bellman_operator,
                                           v=v, max_iter=max_iter, tol=tol,
                                           Tv=Tv)
        sigma = self.compute_greedy(v)

        res = DPSolveResult(v=v,
                            sigma=sigma,
                            num_iter=num_iter,
                            mc=self.controlled_mc(sigma),
                            method='value iteration',
                            epsilon=epsilon,
                            max_iter=max_iter)

        return res