function [x_opt] = Numerical_opt(lambda,B,c_f,c_delta)
% Numerical optimization
% x_opt are the optimal hit probabilities
 
N = length(lambda);  % Number of variables
h_init = ones(N, 1) * B / N;  % Initial guess for h_i, evenly distributed
x_init=h_init;
objective1 = @(x) sum(lambda .*( (1-x)*c_f+x.^2*c_delta/2));  % Objective function
constraint1 = @(x) sum(x) - B; % cache size constraint
% constraint on hit prob. 0<=h_i<=1
nonlcon1 = @(x) deal(constraint1(x), []);  % 
lb = zeros(N, 1);  % Lower bounds
ub = ones(N, 1);     % Upper bounds
options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'sqp');
[x_opt, fval_h] = fmincon(objective1, x_init, [], [], [], [], lb, ub, nonlcon1, options);
end