docs: def robust_rule_simple(self, P_init=None, max_iter=80, tol=1e-8):
          """
          A simple algorithm for computing the robust policy F and the
          corresponding value function P, based around straightforward
          iteration with the robust Bellman operator.  This function is
          easier to understand but one or two orders of magnitude slower
          than self.robust_rule().  For more information see the docstring
          of that method.
  
          Parameters
          ----------
          P_init : array_like(float, ndim=2), optional(default=None)
              The initial guess for the value function matrix.  It will
              be a matrix of zeros if no guess is given
          max_iter : scalar(int), optional(default=80)
              The maximum number of iterations that are allowed
          tol : scalar(float), optional(default=1e-8)
              The tolerance for convergence
  
          Returns
          -------
          F : array_like(float, ndim=2)
              The optimal control matrix from above
          P : array_like(float, ndim=2)
              The positive semi-definite matrix defining the value
              function
          K : array_like(float, ndim=2)
              the worst-case shock matrix K, where
              :math:`w_{t+1} = K x_t` is the worst case shock
  
          """
          # == Simplify names == #
          A, B, C, Q, R = self.A, self.B, self.C, self.Q, self.R
          beta, theta = self.beta, self.theta
          # == Set up loop == #
          P = np.zeros((self.n, self.n)) if P_init is None else P_init
          iterate, e = 0, tol + 1
          while iterate < max_iter and e > tol:
              F, new_P = self.b_operator(self.d_operator(P))
              e = np.sqrt(np.sum((new_P - P)**2))
              iterate += 1
              P = new_P
          I = np.identity(self.j)
          S1 = P.dot(C)
          S2 = C.T.dot(S1)
          K = inv(theta * I - S2).dot(S1.T).dot(A - B.dot(F))
  
          return F, K, P