docs: def optimize(self, n_particles=50, n_iterations=250, restart=1):
  
          """
          the best result of all optimizations will be returned.
          total number of lens models sovled: n_particles*n_iterations
  
          :param n_particles: number of particle swarm particles
          :param n_iterations: number of particle swarm iternations
          :param restart: number of times to execute the optimization;
  
          :return: lens model keywords, [optimized source position], best fit image positions
          """
  
          if restart < 0:
              raise ValueError("parameter 'restart' must be integer of value > 0")
  
          # particle swarm optimization
          penalties, parameters, src_pen_best = [],[], []
  
          for run in range(0, restart):
  
              penalty, params = self._single_optimization(n_particles, n_iterations)
              penalties.append(penalty)
              parameters.append(params)
              src_pen_best.append(self._optimizer.src_pen_best)
  
          # select the best optimization
          best_index = np.argmin(penalties)
  
          # combine the optimized parameters with the parameters kept fixed during the optimization to obtain full kwargs_lens
          kwargs_varied = self._params.argstovary_todictionary(parameters[best_index])
          kwargs_lens_final = kwargs_varied + self._params.argsfixed_todictionary()
  
          # solve for the optimized image positions
          srcx, srcy = self._optimizer.lensing._ray_shooting_fast(kwargs_varied)
          source_x, source_y = np.mean(srcx), np.mean(srcy)
  
          # if we have a good enough solution, no point in recomputing the image positions since this can be quite slow
          # and will give the same answer
          if src_pen_best[best_index] < self._tol_src_penalty:
              x_image, y_image = self.x_pos, self.y_pos
          else:
              # Here, the solver has the instance of "lensing_class" or "LensModel" for multiplane/singleplane respectively.
              print('Warning: possibly a bad fit.')
              x_image, y_image = self.solver.findBrightImage(source_x, source_y, kwargs_lens_final, arrival_time_sort=False)
              #x_image, y_image = self.solver.image_position_from_source(source_x, source_y, kwargs_lens_final, arrival_time_sort = False)
          if self._verbose:
              print('optimization done.')
              print('Recovered source position: ', (srcx, srcy))
  
          return kwargs_lens_final, [source_x, source_y], [x_image, y_image]