How to constrain the every optimized vector component to be nonnegative?

Question

I am building a gradient descent model based on portfolio optimization. Currently, I have finished the model and am able to run it smoothly without any problem. However, there's one issue that I couldn't fix: the optimal vector components have negative values. As my portfolio does not support short-selling stocks, this optimal solution may not be a convincing solution to my portfolio.

Model Description

My cost function for the portfolio problem is defined as

$$F(x)=\frac{\beta_1}{2}(x^T\boldsymbol{\Sigma}x)-\mu^Tx+\frac{\gamma}{2}(e^Tx-1)^2+\lambda(e^Tx-1)+\frac{\beta_2}{2}\Vert{x}\Vert^2_2+\rho\sum_{i=1}^n\max(0,-x_i)^2$$

where $x$ and $\boldsymbol{\Sigma}$ are my minimizer vector and positive definite covariance matrix. $\frac{\gamma}{2}(e^Tx-1)^2+\lambda(e^Tx-1)$ is the augmented lagrangian method for the constraint $e^Tx=1$ and $\rho\sum_{i=1}^n\max(0,-x_i)^2$ is the penalty imposing $x_i>0]$ for all components of $x$.

As can be seen here, my objective is to find the optimum solution for $x$ using the Steepest Descent method. Hence, I start to build my Python code as below:

import numpy as np
from numpy.linalg import norm

def penalty_f(v):
    return  np.sum([np.power(max(0., -xi),2) for xi in v])

def penalty_df(v):
    return -2.*np.array([max(0, -xi) for xi in v])

def aug_lag(lam, v):
    new_lam = lam + gam*(np.sum(v)-1.)
    return new_lam

def lipschitz(v):
    lips = beta_1*np.sqrt(np.trace(cov@cov)) + beta_2*np.sqrt(len(v)) + gam*len(v)
    return min(1,1/lips)

def project_f(v, lam_k):
    func = 0.5*beta_1*v.T@cov@v - b*mean.T@v + gam/2*(np.sum(v)-1.)**2 + lam_k*(np.sum(v)-1.) + beta_2/2*v@v + rho*penalty_f(v)
    return func

def project_df(v, lam_k):
    par_func = beta_1*cov@v - b*mean + gam*np.ones_like(v)*(np.sum(v)-1.) + lam_k*np.ones_like(v) + beta_2*v + rho*penalty_df(v)
    return par_func

def gradient_descent(f, df, ini_v, tolerance, lam, MAX_ITER=10000, output_fname='output.txt'):
    vec = ini_v
    negative = 0
    aug_lam = lam
    with open(output_fname, 'w') as out:
        for i in range(MAX_ITER):
            alpha = lipschitz(vec)
            f_value = f(vec, aug_lam)
            gradient = df(vec, aug_lam)
            direction = np.negative(gradient)
            vec += alpha * direction
            aug_lam = aug_lag(aug_lam, vec)
            print("No. of zeros: ", len(vec)-np.count_nonzero(vec))
            print(norm(gradient,2), "," , f_value)
            msg = f'{f_value},{norm(gradient,2)}, {str(list(vec))}\n'
            out.write(msg)

            # stopping criteria
            if norm(gradient,2) < tolerance:
                print ("The optimum vector for", {df}, " is at ", vec,"at iteration ", i+1)
                for i in vec:
                    if i < 0:
                        negative += 1
                print("No of negative: ", negative)
                print("Vector sum: ",np.sum(vec))
                print("Gradient: ", norm(gradient,2))
                break
            
            if i == MAX_ITER:
                print ('Higher no. of iterations is needed for', {df})
                print ("Vector: ", vec)
                print("Vector sum: ",np.sum(vec))
                print("Gradient: ", norm(gradient,2))
                
    return vec

Here is the result I obtained using 33 stocks:

The optimum vector for {<function project_df at 0x000001FA50FF2D30>}  is at  [-0.03442667 -0.01839447  0.03305939 -0.04319195  0.00353058  0.13003235
  0.09405175  0.00022958  0.08010633 -0.02831581 -0.04066535 -0.02875569
  0.04953858 -0.04019486  0.22866438 -0.03474848 -0.00029308  0.09792178
  0.11669689 -0.08569448  0.00405406  0.05544077  0.12273198  0.09615974
 -0.04374827  0.16994986  0.00670517  0.04772238  0.01729612  0.05218545
  0.03484764 -0.03438397 -0.00811171] at iteration  1230
No of negative:  13
Vector sum:  1.0000000001114249
Gradient:  9.961948573544883e-07

When I assess 30 or more stocks, negative components start to appear. Is there anything I can improve or take note of when performing portfolio optimization?

P.S. The data to the mean and covariance is in the link https://www.dropbox.com/scl/fi/rycj948t4bnq5u60m13ow/Covariance.xlsx?dl=0&rlkey=d8u18ntuxk7wjcl1eup8gmxoa

import pandas as pd

df1 = pd.read_excel('Covariance.xlsx', sheet_name=0, header=None)
df2 = pd.read_excel('Covariance.xlsx', sheet_name=1, header=None)
np_cov_1 = df1.values
mean_1 = df2.values.reshape(len(df2))
ini_vec_1 = np.array([1. / (len(df2)) for i in range(len(df1))])