Solver stopped prematurely. fmincon stopped because it exceeded the function evaluation limit, options.Ma​xFunctionE​valuations = 3.000000e+03.

%constraints c, ceq
% size(D)=[69,69]
%size(sigma)=[69,69] 
%length(h)=50
function[c,ceq]= nlcon (x)
load ('workspacefmincon.mat','n','sigma','D','h','h_min','h_max');
n=length(D);
x=zeros(1,n);
h=linspace(h_min,h_max,50);
for i=1:length(h)
c(i)=h(i)-(x*D*x');
ceq=[];
end

%optimization
x=zeros(1,n);
h=linspace(h_min,h_max,50);
fun=@(x) x*sigma*x';
x_rao=zeros(1,n);
x0_rao=zeros(1,n);
A_rao=[];
b_rao=[];
Aeq_rao=ones(length(h),n);
beq_rao=ones(length(h),1);
l_b_rao=zeros(1,n);
u_b_rao=ones(1,n);
risk_rao = zeros(1,n);
risk_rao = var_min;
constr=@nlcon;
[x_rao, risk_rao]=fmincon(fun,x0_rao,A_rao,b_rao,Aeq_rao,beq_rao,l_b_rao,u_b_rao,constr);
% If I add options, the problem is the same
% options = optimoptions('fmincon','Display','iter','Algorithm','sqp');

Answer 1

[Stephan]

You have to give the options to the solver and to insert the MaxFunctionEvaluations Option into the optimoptions struct:

MyValue = 10e4;
options = optimoptions('fmincon','Display','iter','Algorithm','sqp', 'MaxFunctionEvaluations',MyValue);
[x_rao, risk_rao]=fmincon(fun,x0_rao,A_rao,b_rao,Aeq_rao,beq_rao,l_b_rao,u_b_rao,constr,options);

[Ylenia Placella]

%after your suggestion it gives me this error

Converged to an infeasible point.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance but constraints are not
satisfied to within the value of the constraint tolerance.

<stopping criteria details>
>> 

[Stephan]

Are you sure your problem has a solution under the given resitrictions? We can not try, because you did not provide the workspace data.

[Stephan]

Ylenias "answer" moved here:

this is the workspace

[Stephan]

You override x in every call of your nonlinear constraint function. Also the result of:

x*D*x'

with size of x 1x69 and D size of 69x69 is a scalar. So why do you loop like this?:

for i=1:length(h)
c(i)=h(i)-(x*D*x')
ceq=[];
end

Thats confusing - Please show the constraints you want to write in a mathematical representation and/or tell us more about what you want to do here exactly.

[Stephan]

Ylenias "answer" moved here:

I have to solve this optimization problem:

min x*sigma*x

h-xDx ≤0

where h is:

h=linspace(h_min,h_max,50)

[Stephan]

You mean instead i think. But since is a scalar the condtion is always met, when h_max is used - If it is true for all h - isnt it?

The condtion is met by your bounds,

And is true when Aeq and beq are

Aeq_rao=ones(1,n);
beq_rao=1;

instead of

Aeq_rao=ones(length(h),n);
beq_rao=ones(length(h),1);

So le me assume that my thoughts are correct, then we end up here:

load('workspacefmincon.mat','n','sigma','D','h_max');

%optimization
fun=@(x) x*sigma*x';
x0_rao=rand(1,n);
A_rao=[];
b_rao=[];
Aeq_rao=ones(1,n);
beq_rao=1;
l_b_rao=zeros(1,n);
u_b_rao=ones(1,n);
constr=@(x)nlcon(x,D,h_max);

options = optimoptions('fmincon','Display','final-detailed','Algorithm','sqp',...
    'MaxFunctionEvaluations', 5e5,'MaxIterations',5e4,'ConstraintTolerance',...
    1e-4);
[x_rao, risk_rao,exitflag,output]=fmincon(fun,x0_rao,A_rao,b_rao,Aeq_rao,...
    beq_rao,l_b_rao,u_b_rao,constr,options)

function[c,ceq]= nlcon(x,D,h_max)
    c=h_max-(x*D*x');
    ceq=[];
end

Running this code gives still no feasible points. I tried using ga, which also did not find a feasible solution. The sum of x_rao is:

>> sum(x_rao)

ans =

    1.1263

So the bounding condition is not met. I dont think that the problem can be solved under the given constraints - or there is still an error somewhere, that can not be found without insight in the problem.

[Ylenia Placella]

the problem persist.

I have to find the vector x that solve the problem.

i find the h with this code:

index=xlsread('C:\Users\desyp\Desktop\Tesi finass\NASDAQ100.xlsx')
P = index
RR = diff(P)./P(1:end-1,:);
sigma=cov(RR)
rho=corrcoef(RR)
mu=mean(RR)
options=optimset('algorithm','interior-point-convex','MaxIter',1.e7,'TolFun',1.e-10,'TolX',1.e-10)
n=length(mu)
D=1-rho
H=2*D
f=zeros(n,1)
Aeq=ones(1,n)
beq=1
l_b=zeros(1,n)
[x_h_min,var_min]=quadprog(H,f,[],[],Aeq,beq,l_b,[],[],options)
sum_xhmin=sum(x_h_min)
x_hmin_unit=x_h_min/sum_xhmin
h_min=x_hmin_unit'*D*x_hmin_unit
[x_h_max,var_max]=quadprog(-H,f,[],[],Aeq,beq,l_b,[],[],options)
sum_xhmax=sum(x_h_max)
x_hmax_unit=x_h_max/sum_xhmax
h_max=x_hmax_unit'*D*x_hmax_unit
h=linspace(h_min,h_max,50)

[Stephan]

This appears to work - i ran it as 1 script:

index=xlsread('C:\Users\desyp\Desktop\Tesi finass\NASDAQ100.xlsx');
P = index;
RR = diff(P)./P(1:end-1,:);
sigma=cov(RR);
rho=corrcoef(RR);
mu=mean(RR);
n=length(mu);
D=1-rho;
H=2*D;
f=zeros(n,1);
Aeq=ones(1,n);
beq=1;
l_b=zeros(1,n);
x0 = rand(1,n);

options=optimoptions('quadprog','algorithm','active-set','MaxIter',1.e7,...
    'TolFun',1.e-10,'TolX',1.e-10);

[x_h_min,var_min]=quadprog(H,f,[],[],Aeq,beq,l_b,[],x0,options);
sum_xhmin=sum(x_h_min);
x_hmin_unit=x_h_min/sum_xhmin;
h_min=x_hmin_unit'*D*x_hmin_unit;

[x_h_max,var_max]=quadprog(-H,f,[],[],Aeq,beq,l_b,[],x0,options);
sum_xhmax=sum(x_h_max);
x_hmax_unit=x_h_max/sum_xhmax;

h_max=x_hmax_unit'*D*x_hmax_unit;
h=linspace(h_min,h_max,50);


%% optimization
fun=@(x) x*sigma*x';
x0_rao=rand(1,n);
A_rao=[];
b_rao=[];
Aeq_rao=ones(1,n);
beq_rao=1;
l_b_rao=-zeros(1,n);
u_b_rao=ones(1,n);
constr=@(x)nlcon(x,D,h_max);

opts = optimoptions('fmincon','Display','final-detailed','Algorithm','sqp',...
    'MaxFunctionEvaluations', 5e4,'MaxIterations',5e3,'ConstraintTolerance',...
    1e-6);
[x_rao, risk_rao,exitflag,output]=fmincon(fun,x0_rao,A_rao,b_rao,Aeq_rao,...
    beq_rao,l_b_rao,u_b_rao,constr,opts)


function[c,ceq]= nlcon(x,D,h_max)
c=h_max-(x*D*x');
ceq=[];
end

But there is still a warning on the first call of quadprog, that the problem is non-convex. Also note that i assumed that my conclusion regarding the nonlinear constraint function is correct. You have to check if that can be correct.

[Ylenia Placella]

%if now I want to calculate X and var_RAO for every h. How can I do?
%I used this code. Does not mark error, but X and var_RAO are matrices with all zero.
load ("work_rao_completo.mat","sigma","risk_rao","x_rao","n","D","x0_rao","h_max","h_min")
h=linspace(h_min,h_max,50);
X_rao=x_rao'
func=@(x) X(:,i)'*sigma*X(:,i)
X = zeros(n,length(h))
var_RAO = zeros(1,length(h))
var_RAO(1)=risk_rao
vincoli=@(x)con(x,D)
function[c_i,ceq]= con(x,D)
for i=1:length(h)
    c_i= h(i)-X(:,i)'*D*X(:,i)
ceq=[];
[X(:,i),var_RAO(i),exitflag,output]=fmincon(func,x0_rao,A_rao,b_rao,Aeq_rao,...
    beq_rao,l_b_rao,u_b_rao,vincoli,opts)
end
end

Answer 2

[Ylenia Placella]

%if now I want to calculate X and var_RAO for every h. How can I do?
%I used this code. Does not mark error, but X and var_RAO are matrices with all zero.

load ("work_rao_completo.mat","sigma","risk_rao","x_rao","n","D","x0_rao","h_max","h_min")
h=linspace(h_min,h_max,50);
X_rao=x_rao'
func=@(x) X(:,i)'*sigma*X(:,i)
X = zeros(n,length(h))
var_RAO = zeros(1,length(h))
var_RAO(1)=risk_rao
vincoli=@(x)con(x,D)

function[c_i,ceq]= con(x,D)
for i=1:length(h)
    c_i= h(i)-X(:,i)'*D*X(:,i)
ceq=[];
[X(:,i),var_RAO(i),exitflag,output]=fmincon(func,x0_rao,A_rao,b_rao,Aeq_rao,...
    beq_rao,l_b_rao,u_b_rao,vincoli,opts)
end
end