fmincon: Complex optimization, output error

Hello, I am trying to implement a minimization problem provided in this paper (see "Main_Article" in the attachment). More specifically I am trying to implement eq. (7) of the paper for the case vp = 1. In other words this means that the problem should be convex and therefore I can solve it with fmincon. However, I get te following error on fmincon:

Converged to an infeasible point.

fmincon stopped because the size of the current step is less than
the default value of the step size tolerance but constraints are not
satisfied to within the default value of the constraint tolerance.

This is my code (also in the attachment):

main.m

clear all

phi0 = [0.7455,0.2545,0.9475,0.0525];

A = [];
b = [];
Aeq = [];
beq = [];
lb = 1 * ones(4);
up = 5 * ones(4);


phi = fmincon('obj', phi0, A, b, Aeq, beq, lb, up, 'nlcon');

obj.m

function [J] = obj(phi)

    phi = reshape(phi,2,2);

    %q
    q1 = 0.1481;
    q2 = 0.8519;
    q = [q1, q2];

    %v
    vb = 0;
    vp = 1;

    %fb
    fb1 = 0;
    fb2 = 0;
    fb = [fb1, fb2];

    %fp
    fp = vp * phi;

    %alpha
    a11 = 0.0011;
    a12 = 0.0055;
    a21 = 0.003;
    a22 = 0.0043;
    a = [a11, a12; a21, a22];

    %beta
    b11 = 2.5485;
    b12 = 2.1929;
    b21 = 1.5363;
    b22 = 2.2671;
    b = [b11, b12; b21, b22];

    %fs
    fs = fb + fp;

    %OBJECTIVE J
    l = a .* fs + b;    %Delay function l

    J1 = fs .* l;       %Social cost J1
    J1 = sum(J1,2);

    J = q * J1;       %Expected social cost J

nlcon.m

function [c, ceq] = nlcon(phi)
    
    phi = reshape(phi,2,2);
    
    %Define variables
    n_links = 2;
    n_states = 2;

    %q
    q1 = 0.1481;
    q2 = 0.8519;
    q = [q1, q2];

    %v
    vb = 0;
    vp = 1;

    %fb
    fb1 = 0;
    fb2 = 0;
    fb = [fb1, fb2];

    %fp
    fp = vp * phi;

    %alpha
    a11 = 0.0011;
    a12 = 0.0055;
    a21 = 0.003;
    a22 = 0.0043;
    a = [a11, a12; a21, a22];

    %beta
    b11 = 2.5485;
    b12 = 2.1929;
    b21 = 1.5363;
    b22 = 2.2671;
    b = [b11, b12; b21, b22];

    %fs
    fs = fb + fp;
    
    
    %Constraint (1)
    ceq_i = 1 - sum(phi,2);
    ceq(1) = ceq_i(1);
    ceq(2) = ceq_i(2);

    %Constraint (2)
    index_c = 0;

    for j=1:n_links
        for k=1:n_links
            if k ~= j
                index_c = 1 + index_c;
                A1 = q * ( a(:,j) .* phi(:,j) .* fs(:,j) );
                A2 = q * ( b(:,j) .* phi(:,j));

                B1 = q * ( a(:,k) .* phi(:,j) .* fs(:,k) );
                B2 = q * ( b(:,k) .* phi(:,j));
                c_i(index_c) = A1 + A2 -(B1 + B2);
            end
        end
    end
    c(1) = c_i(1);
    c(2) = c_i(2);