function optimization(data, n, k, ρ)
    data_mean = mean(data)
    data_std = std(data)
    norm_data = (data .- data_mean) ./ data_std

    m = Model(solver=GurobiSolver(OutputFlag=0))
    @variable(m, d)
    @variable(m, x1[1:n], Bin)
    @variable(m, x2[1:n], Bin)
    mu1 = 1 / k * sum(norm_data .* x1)
    mu2 = 1 / k * sum(norm_data .* x2)
    sigma1 = 1 / k * sum((norm_data .^ 2) .* x1)
    sigma2 = 1 / k * sum((norm_data .^ 2) .* x2)
    @constraint(m, d >= mu1 - mu2 + ρ*sigma1 - ρ*sigma2)
    @constraint(m, d >= mu1 - mu2 + ρ*sigma2 - ρ*sigma1)
    @constraint(m, d >= mu2 - mu1 + ρ*sigma1 - ρ*sigma2)
    @constraint(m, d >= mu2 - mu1 + ρ*sigma2 - ρ*sigma1)
    @constraint(m, sum(x1) == k)
    @constraint(m, sum(x2) == k)
    @constraint(m, x1 .+ x2 .== 1)
    @constraint(m, x1[1] == 0)
    @objective(m, Min, d)
    @show status=solve(m)
    x1 = getvalue(x1)
    x2 = getvalue(x2)

    first_group, second_group = zeros(0), zeros(0)
    for i = 1:n
        if x1[i] == 1
            append!(first_group, data[i])
        else
            append!(second_group, data[i])
        end
    end
    return first_group, second_group
end