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]