SDP Problem: No Correct Solutions Obtained using cvxpy in Python

Question

I'm working on an optimization problem using cvxpy in Python. The optimization problem is based on a model presented in a IEEE Transactions on Automatic Control journal article：On Maximizing the Second Smallest Eigenvalue of a State-Dependent Graph Laplacian.The objective is to maximize the variable gamma while satisfying a set of constraints. However, the code I have doesn't seem to produce any solutions. I'm wondering if there's an issue with the formulation of the constraints or the way I'm setting up the problem.

Optimization Problem Model

I have implemented the code and attempted to solve the optimization problem, but the solver doesn't return any feasible solutions.

---

import cvxpy as cp
import numpy as np
from matplotlib import pyplot as plt
from scipy.linalg import null_space


def solve_optimization_problem(n, x_init):
    # Weight model parameters
    epsilon = 0.1
    rho1 = 1
    rho2 = 1.5

    # Variables in the cvxpy optimization problem
    gamma = cp.Variable()  # Optimization objective
    x = cp.Variable((n, 2))  # Position matrix
    D = cp.Variable((n, n))  # Euclidean distance matrix
    Aij = cp.Variable((n, n))  # Adjacency matrix, used here to calculate the graph Laplacian matrix
    L_G = cp.Variable((n, n))  # Graph Laplacian matrix

    # Calculate the Euclidean distance matrix and weight matrix in the initial state
    D_initial = np.zeros((n, n))
    W_initial = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            if i != j:
                D_initial[i, j] = (np.linalg.norm(x_init[i, :] - x_init[j, :])) ** 2
                W_initial[i, j] = epsilon ** ((rho1 - D_initial[i, j]) / (rho1 - rho2))
            else:
                D_initial[i, j] = 0
                W_initial[i, j] = 0

    # Create constraint list
    constraints = []

    # Add constraints
    for i in range(n):
        for j in range(n):
            if i != j:
                # Linearized collision avoidance constraint
                constraints.append(
                    D[i, j] == 2 * (x[i, :] - x[j, :]).T @ (x_init[i, :] - x_init[j, :]) - D_initial[i, j])
                constraints.append(D[i, j] >= 1)

                # Update adjacency matrix
                constraints.append(
                    Aij[i, j] == W_initial[i, j] - epsilon ** ((rho1 - D_initial[i, j]) / (rho1 - rho2)) * (
                            D[i, j] - D_initial[i, j]))
            else:
                constraints.append(Aij[i, j] == 0)
                constraints.append(D[i, j] == 0)

    # Graph Laplacian matrix
    constraints.append(L_G == cp.diag(cp.sum(Aij, axis=1)) - Aij)

    # Constraint condition
    P = null_space(np.ones((1, n)))
    constraints.append(P.T @ L_G @ P >> gamma * np.eye(n - 1))

    # Euclidean matrix constraint
    J = np.eye(n) - np.ones((n, n)) / n
    constraints.append(J @ D @ J << 0)

    # Define optimization problem
    problem = cp.Problem(cp.Maximize(gamma), constraints)

    # Solve optimization problem
    # problem.solve(cp.MOSEK,verbose=True)
    problem.solve()

    # Return optimization results
    if problem.status in ['optimal', 'optimal_inaccurate']:
        return gamma.value, x.value, D.value, L_G.value
    else:
        print("Optimization failed with status", problem.status)
        return None, None, None, None


if __name__ == "__main__":
    # Number of robots
    n = 6

    # Initial position
    x_init = np.array([[0, 0], [0, 1], [0, 2], [0, 3], [0, 4], [0, 5]])

    # Store positions for each iteration
    positions = [x_init]

    # Number of iterations
    iterations = 100

    # Optimize position in each iteration
    for _ in range(iterations):
        result = solve_optimization_problem(n, positions[-1])
        if result[0] is not None:
            gamma, x, D, L_G = result
            print("Optimal gamma:", gamma)
            print("Optimal positions:", x)
            positions.append(x)
        else:
            print("Failed to solve the optimization problem.")
            break

    # Plot the trajectory of each robot
    for i in range(n):
        x_coords = [pos[i, 0] for pos in positions]
        y_coords = [pos[i, 1] for pos in positions]
        plt.plot(x_coords, y_coords, marker='o')

    plt.show()

Could someone please review the code and provide insights into what might be causing the issue? Any suggestions on how to modify the code to make it work would be greatly appreciated.