function energie = U(c, I_moyen, beta, gamma, S, R)

    N = length(I_moyen);
    
    Ui = 1 - 2 ./ ( 1 + exp( - gamma * (I_moyen / S - 1) ) );

    delta = 0;
    for i=1:N
        for j=i+1:N
            if norm(c(j,:) - c(i,:)) <= sqrt(2) * R
                delta = delta + 1;
            end
        end
    end

    energie = sum(Ui) + beta * delta;

end