How can I solve this fincon problem?

% I have to solve this problem

min x*sigma*x

h-xDx ≤0

where h is:

h=linspace(h_min,h_max,50)

close all;clear all;clc;
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);
h;
%% 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);
opts = optimoptions('fmincon','Display','final-detailed','Algorithm','sqp',...
    'MaxFunctionEvaluations', 5e4,'MaxIterations',5e3,'ConstraintTolerance',...
    1e-6);
x_rao=zeros(n,length(h))
risk_rao=zeros(1,length(h))
[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)
for i=1:length(h)
c=h(i)-(x*D*x');
ceq=[];
[x_rao(:,i), risk_rao(i),exitflag,output]=fmincon(fun,x0_rao,A_rao,b_rao,Aeq_rao,...
    beq_rao,l_b_rao,u_b_rao,constr,opts)
end
end

[Stephan]

New problem? Is it different to your other?

[Ylenia Placella]

%I have a vector h , so I used a cycle for but it doesn't work.
%I have to find a vector x with dimension 1xn for each h, with h=linspace(h_min,h_max,50);
% So x won't be a vector but a matrix with dimension nxh
%I used this cose:
close all;clear all;clc;
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);
h;
%% 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);
opts = optimoptions('fmincon','Display','final-detailed','Algorithm','sqp',...
    'MaxFunctionEvaluations', 5e4,'MaxIterations',5e3,'ConstraintTolerance',...
    1e-6);
x_rao=zeros(n,length(h))
risk_rao=zeros(1,length(h))
risk_rao(1)=var_min
function[c,ceq]= nlcon(x,D,h)
for i=1:length(h)
c(i)=h(i)-(x*D*x');
ceq=[];
[x_rao(:,i), risk_rao(i),exitflag,output]=fmincon(fun,x0_rao,A_rao,b_rao,Aeq_rao,...
    beq_rao,l_b_rao,u_b_rao,constr,opts)
end
end

Answer 1

[Walter Roberson]

constr=@(x)nlcon(x,D,h);

That is in the script

risk_rao(1)=var_min

The main part of the script ends with that line.

function[c,ceq]= nlcon(x,D,h)
for i=1:length(h)
c(i)=h(i)-(x*D*x');
ceq=[];
[x_rao(:,i), risk_rao(i),exitflag,output]=fmincon(fun,x0_rao,A_rao,b_rao,Aeq_rao,...
    beq_rao,l_b_rao,u_b_rao,constr,opts)
end
end

Inside the for loop that is inside nlcon, you try to call fmincon(), calculating some outputs. None of the outputs are returned from nlcon, though, so provided the fmincon calls do not error, then there is no point in calling fmincon there.

The call to fmincon needs variables x0_rao, A_rao, b_rao, Aeq_rao, beq_rao, l_b_rao, u_b_rao, constr, and opts, none of which are defined inside the workspace of nlcon .

It would seem to make more sense to have the for loop outside of nlcon .

Note: x*D*x' does not depend upon h inside nlcon, so instead of that for loop you might as well code

c = h - (x*D*x');

with no loop.

On the other hand, since h is strictly increasing, the maximum value can be used instead of the vector:

c = h(end) - (x*D*x');

[Stephan]

We already discussed this point here

[Ylenia Placella]

I have to make the loop because c depend upon h.

For heach h I will have a different c, then for each c I have to solve the problem.

For this reason I have to use the loop.

But I don't know how.

[Walter Roberson]

I have to make the loop because c depend upon h.

You have

c(i)=h(i)-(x*D*x');

as you go through the for i loop, (x*D*x') does not change, so you can move that outside the loop:

xDx = X*D*x';
for i=1:length(h)
  c(i) = h(i) - xDx;
end

h is a vector, not a function. So you can vectorize that loop:

xDx = X*D*x';
c = h - xDx;

Now, h is

h=linspace(h_min,h_max,50);

Unless h_max somehow became less than h_min, then the smallest value of h will be h(1) and the largest value of h will be h(end)

Suppose h is 1 2 3 4 5 and xDx is 1/2 . Then

c(1) = h(1) - 1/2 -> 1 - 1/2 -> 1/2
c(2) = h(2) - 1/2 -> 2 - 1/2 -> 1 1/2
c(3) = h(3) - 1/2 -> 3 - 1/2 -> 2 1/2
c(4) = h(4) - 1/2 -> 4 - 1/2 -> 3 1/2
c(5) = h(5) - 1/2 -> 5 - 1/2 -> 4 1/2

And compare that to the vectorized version

c = h - 1/2 -> c = [1 2 3 4 5] - 1/2 -> [1/2, 1 1/2, 2 1/2, 3 1/2, 4 1/2]

which is exactly the same.

Therefore you can vectorize.

In order for the vector of constraints to be satisfied for nonlinear constraint purposes, h - xDx must be <= 0 for all entries. Which can be rewritten as h <= xDx for all entries. But that can only be true if max(h) <= xDx . The opposite is not the case -- min(h) might be <= xDx but some other value in h might be greater than xDx (suppose xdx was 3, then [1 2 3 4 5] <= 3 is [1 2 3 4 5] - 3 <= 0 is [-2 -1 0 1 2] <= 0 which is true for the first fiew entries but not for all of them. But if it is true for the largest of the entries in h, then it will be true for all of the entries in h.

Therefore, instead of using

c = h - xDx

you can simplify to

c = h(end) - xDx 

No loop needed.

Inside that for i loop you currently have

[x_rao(:,i), risk_rao(i),exitflag,output]=fmincon(fun,x0_rao,A_rao,b_rao,Aeq_rao,...
    beq_rao,l_b_rao,u_b_rao,constr,opts)

and I showed that most of those variables are not defined inside your nonlinear constraint function. If they were defined inside nlcon so that the call could go ahead, notice that none of the output variables x_rao, risk_rao, exitflag, or output, affect c or ceq . Therefore calling fmincon() inside the nlcon is just a waste of time.