Problem

The classical Markowitz problem is as follows:
$$\underset{x}{\max }{r}^{T}x-\lambda x^{T}Qx$$
The quadprog solver addresses this quadratic programming problem. Some cardinality constraints would be added. $x$ is the vector of asset allocation fractions, with $0\le x(i)\le 1$ for each $i$. The indicator variables $v$ such that $v(i)=0$ when $x(i)=0$, and $v(i)=1$ when $x(i)>0$. To get variables that satisfy this restriction, set the $v$ vectror to be a binary variable, and impose the linear constraints.
$$v(i)f_{\min }\le x(i)\le v(i)f_{\max }$$
and
$$m\le \sum v(i)\le M$$

Objective and Successive linear Approximations

The objective function is
$$\underset{x}{\min }\lambda x^{T}Qx-r^{T}x$$
However, the intlinprog MILP solver requires a linear objective function. There is a standard technique to reformulate this problem into one with linear objective and nonlinear constraints. Introduce a slack variable $z$ to represent the quadratic term.
$$\begin{gathered} \underset{x,z}{\min }\lambda z-r^{T}x \\ \{\begin{array}{l}{x}^{T}Qx-z\le 0\\ z\ge 0\end{array} \end{gathered}$$
To solve MILP approximations iteratively, we could approximate the nonlinear constraint locally near the current point
$$x^{T}Qx-z=x_0^{T}Q{x}_0+2x_0^{T}Q\delta -z+O(|\delta {|}^2)$$
Replacing $\delta$ by $x-x_0$
$$x^{T}Qx-z=-x_0^{T}Q{x}_0+2x_0^{T}Qx-z+O(|x-x_0{|}^2)$$
For each intermediate solution $x_k$ you introduce a new linear constraint in $x$ and $z$ as the linear part of the expression above:
$$-x_k^{T}Q{x}_k+2x_k^{T}Qx-z\le 0$$

Code

%% integer programming
clc;
clear all;
close all;

%%
load port5;
r = mean_return;
Q = Correlation .* (stdDev_return * stdDev_return');
N = length(r);
% indexes
xvars = 1:N;
vvars = N+1:2*N;
zvar = 2*N+1;

lb = zeros(2*N+1, 1);
ub = ones(2*N+1, 1);
ub(zvar) = Inf;

% set the number of assets in the solution to between 100 and 150.
M = 150;
m = 100;
A = zeros(1, 2*N+1);
A(vvars) = 1; % A*x represents the sum of v(i)
A = [A; -A];
b = zeros(2, 1);
b(1) = M;
b(2) = -m;

% set the minimal nonzero fraction of assets to be 0.05 and the maximal to
% be 0.20
fmin = 0.001;
fmax = 0.05;

% inequalities fmin(i)*v(i)<=x(i)<= fmax(i)*v(i)
Atemp = eye(N);
Amax = horzcat(Atemp, -Atemp*fmax, zeros(N, 1));
A = [A; Amax];
b = [b; zeros(N, 1)];
Amin = horzcat(-Atemp, Atemp*fmin, zeros(N, 1));
A = [A; Amin];
b = [b; zeros(N, 1)];

% constraint sum(x_i)=1
Aeq = zeros(1, 2*N+1);
Aeq(xvars) = 1;
beq = 1;

% risk-aversion coefficient lambda
lambda = 100;

% objective function
f = [-r; zeros(N, 1); lambda];

%% solve
ops = optimoptions(@intlinprog, 'Display', 'off');
[xLinInt, fval, exitFlag, output] = intlinprog(f, vvars, A, b, Aeq, beq, lb, ub, ops);

%% stopping condition
thediff = 1e-4;
iter = 1;
assets = xLinInt(xvars); 
trueVar = assets'*Q*assets;
zslack = xLinInt(zvar);

cvx code

%% MLIP
n = 10;
k_high = 8;
k_low = 2;
fmax = 0.30;
fmin = 0.05;
mu = abs(randn(n, 1));
Sigma = randn(n, n)/10;
Sigma = Sigma'*Sigma;

cvx_solver mosek
cvx_begin
    variable w(n);
    variable bi(n) binary;
    ret = mu'*w;
    risk = quad_form(w, Sigma);
    maximize(ret-risk);
    subject to 
        ones(1, n)*w==1;
        w>=0;
        w<=bi;
        % k_high>=nnz(w)>=k_low;
        % norm(w, Inf)<=fmax;
        fmin*bi<=w;
        fmax*bi>=w;
        sum(bi)>=k_low;
        sum(bi)<=k_high;
cvx_end

w=value(w);

Reference

Mathlab MILP