how to code epsilon constraint model in Xpress

Question

I want to used the ε-constrained method to solve a problem a posteriori and plot the set of Pareto optimal solutions and the Pareto front.

I think I've managed to do this as my results look similar to what I achieved with the weighted sum method, although I need to discuss the efficiency of each method on obtaining the Pareto front and am completely stuck on what to say.

Any help would be really appreciated!

When I use very large values of epsilon it seems to look "normal" however when I restrict it to any smaller than the range of values for one of my function it becomes completely different.

m_e = xp.problem(name='Prob1')

### Create decision variables and add them to the model ###
x = np.array([xp.var('x1',lb=-20), xp.var('x2',ub=20)], dtype=xp.npvar)
m_e.addVariable(x)



# Define objective vector
f1 = 2 + (x[0]-2)**2 + (x[1]-1)**2

f2 = 5*x[0] + (x[1]-5)**2

### Define constraints and add them to the model

m_e.addConstraint(x[0] -(3*x[1]) + 10 <= 0)
m_e.addConstraint(x[0]**2 + x[1]**2 <= 225)
        ###
    
### Set paramaters and create arrays for storing solutions


# Array containing the pareto front
front = np.zeros((iterations,2))
# Array containing the pareto optimal points
pareto = np.zeros((iterations,2))
# Compute all the weights
         ###
# Create plots  
fig, ax = plt.subplots(1,2)
iterations = 25

# loop through the number of required itterations
for i in range(iterations):
    m_e.addConstraint(5*x[0] + (x[1]-5)**2 <= epsilon[i])
    # Update objective value 
    m_e.setObjective(f1, sense = xp.minimize)
    # Solve the model
    m_e.solve()
    
    # Store the pareto optimal point
    pareto[i] = m_e.getSolution(x)
   
    # Find the corresponding solutions in objective space
    front[i,0] = xp.evaluate(f1,problem=m_e)
    front[i,1] = xp.evaluate(f2,problem=m_e)

    # Plot results
    ax[0].scatter(pareto[i,0],pareto[i,1],color='blue')
    ax[1].scatter(front[i,0],front[i,1],color='red')
    

### Set plot title ect ### 
ax[0].set_ylabel(r'$x_2$')
ax[0].set_xlabel(r'$x_1$')
ax[0].set_title('Decision space space')

ax[1].set_ylabel(r'$f_2(x)$')
ax[1].set_xlabel(r'$f_1(x)$')
ax[1].set_title('Objective space')

fig.tight_layout(pad=1.0)```