GA, Could not find a feasible initial point

I'm triying to make a portfolio optimization using genetic algorithm.

This is the model that I made, with some constraints, but when I run it, I can't reach a solution, just appear "Could not find a feasible initial point".

Someone could help me fixing my code or helping me to find the mistake.

Thanks.

data = xlsread('Returns.xlsx');

nAssets = size(data, 2);

%Returns and covariance
mu = mean(data);
sigma = cov(data);

%formulate the problem/optimization
r_target = 0.10; %r_target is the required return
f = zeros(nAssets, 1);    %there is no constant
A = [-mu; -eye(nAssets)]; %besides the returns we forbid short selling
b = [-r_target; zeros(nAssets, 1)]; % required return and weights greater/eqauls 0
Aeq = ones(1, nAssets);   %All weights should sum up...
beq = 1;                  %... to one (1)

%solve the optimization
w = ga(@MaxReturn,nAssets,A,b,Aeq,beq);
w

%print the solution
fprintf(2, 'Risk: %.3f%%\n', sqrt(w'*sigma*w)*100);
fprintf(2, 'Ret: %.3f%%\n', w'*mu'*100);


function f = MaxReturn(w,mu);
f = mtimes(w'*mu');
end

[John D'Errico]

I'm not sure what you think this does?

f = mtimes(w'*mu');

Well, besides generate an error, since mtimes is not a monadic operator.

w'*mu' is a valid operation. But then to try to pass that to mtimes, it will fail.

Anyway, without your data, we cannot even test your code.

[Santiago]

Hi John, thanks for your reply.

With that line I want multiple the GA solution with the means of assets, so that result have to be a number (1x1), not a matrix. Also, that is my objetive function of the algorithm, so I wanto to maximize it.

Oh sorry, I already uploaded the data

[John D'Errico]

My point is, IF you want to multiply w' with mu', then just write

f = w'*mu';

The mtimes there as you wrote it is redundant, thus mtimes(w'*mu'), and in fact, will generate a runtime error.

[Walter Roberson]

ga() always passes in row vectors for the unknowns. If you transpose that you would get a column vector, and matrix multiply of a column vector on the left is not going to result in a scalar as required for ga.

The corrected code is in my Answer: f = w*mu'

Answer 1

[Walter Roberson]

format long g
filename = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/835520/Returns.xlsx';

data = readmatrix(filename);

nAssets = size(data, 2);

%Returns and covariance
mu = mean(data);
mu.'

ans = 8×1

0.0017001795070217 0.00238687944682926 0.00534690780710744 0.00276663664554018 7.24577440764558e-06 0.00477387349434592 0.00149865368592419 0.00130429237016939

sigma = cov(data);

%formulate the problem/optimization
r_target = 0.10; %r_target is the required return
f = zeros(nAssets, 1);    %there is no constant
A = [-mu; -eye(nAssets)] %besides the returns we forbid short selling

A = 9×8

-0.0017001795070217 -0.00238687944682926 -0.00534690780710744 -0.00276663664554018 -7.24577440764558e-06 -0.00477387349434592 -0.00149865368592419 -0.00130429237016939 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1

b = [-r_target; zeros(nAssets, 1)] % required return and weights greater/eqauls 0

b = 9×1

-0.1 0 0 0 0 0 0 0 0

Aeq = ones(1, nAssets)   %All weights should sum up...

Aeq = 1×8

1 1 1 1 1 1 1 1

beq = 1                  %... to one (1)

beq = 
     1

%solve the optimization
fcn = @(w)MaxReturn(w,mu);
[w, fval, flag, output] = ga(fcn, nAssets, A, b, Aeq, beq)

Could not find a feasible initial point.

w =

     []


fval =

     []

flag = 
    -2

output = struct with fields:
      problemtype: 'linearconstraints'
         rngstate: [1×1 struct]
      generations: 0
        funccount: 0
          message: 'Could not find a feasible initial point.'
    maxconstraint: []

if isempty(w)
   warning('could not find any solution')
else
    %print the solution
    fprintf(2, 'Risk: %.3f%%\n', sqrt(w'*sigma*w)*100);
    fprintf(2, 'Ret: %.3f%%\n', w'*mu'*100);
end

Warning: could not find any solution

[wF, fvalF, flagF, outputF] = fmincon(fcn, .05*ones(1,nAssets), A, b, Aeq, beq)

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.


Consider enabling the interior point method feasibility mode.

wF = 1×8

-0.000345128325127368 -0.000280137675695389 1.00268187780506 -0.000244196698363124 -0.000505350943988705 -5.42278068495475e-05 -0.000364201130047243 -0.000382595871917903

fvalF = 
       0.00535800915538925

flagF = 
    -2

outputF = struct with fields:
         iterations: 306
          funcCount: 2767
    constrviolation: 0.0946419908446108
           stepsize: 3.97580668930409e-25
          algorithm: 'interior-point'
      firstorderopt: 0.00284530672744689
       cgiterations: 343
            message: '↵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>↵↵Optimization stopped because the relative changes in all elements of x are↵less than options.StepTolerance = 1.000000e-10, but the relative maximum constraint↵violation, 9.464199e-02, exceeds options.ConstraintTolerance = 1.000000e-06.↵↵'
       bestfeasible: []

function f = MaxReturn(w,mu)
  f = w * mu';
end

What's happening is that ga() decides that the constraints cannot be satisfied by any sets of values, and so never calls the function at all. fmincon() cannot find a feasible point either.

[Santiago]

Hi Walter, thank you for your reply.

So, I typed the constraints wrong. May you help to to write then correctly?

I have to maximize the portfolio return, what's the multiplication between w, that are the percentages of assets in the portfolio, and mu, that are the average of asset returns. Also, w is the main result that I need from the genetic algorithm.

The constraints are:

• The sum of all the percentages of the assets have to add up to 1 (100%).
• The percentage of each asset must be greater than 0 and less than 1.

I tried my best, but I'm not that good with MatLab. If you may help me it could be amazing.

[Walter Roberson]

The sum of the coefficients must be 1, right? Now suppose that you put all of that into a single coefficient and made the others zero, then what is the maximum w * mu sum that you could have? It would occur when you put the entire "1" into the largest element of mu and let the other elements be 0. So let us check, is max(mu) >= r_target ? NO! max(mu) is about 0.005 but r_target is 0.1 .

The configuration is impossible.

format long g
filename = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/835520/Returns.xlsx';

data = readmatrix(filename);

nAssets = size(data, 2);

%Returns and covariance
mu = mean(data);
mu.'

ans = 8×1

0.0017001795070217 0.00238687944682926 0.00534690780710744 0.00276663664554018 7.24577440764558e-06 0.00477387349434592 0.00149865368592419 0.00130429237016939

sigma = cov(data);

%formulate the problem/optimization
r_target = 0.10; %r_target is the required return
f = zeros(nAssets, 1);    %there is no constant
A = [-mu; -eye(nAssets)] %besides the returns we forbid short selling

A = 9×8

-0.0017001795070217 -0.00238687944682926 -0.00534690780710744 -0.00276663664554018 -7.24577440764558e-06 -0.00477387349434592 -0.00149865368592419 -0.00130429237016939 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1

b = [-r_target; zeros(nAssets, 1)] % required return and weights greater/eqauls 0

b = 9×1

-0.1 0 0 0 0 0 0 0 0

Aeq = ones(1, nAssets)   %All weights should sum up...

Aeq = 1×8

1 1 1 1 1 1 1 1

beq = 1                  %... to one (1)

beq = 
     1

%solve the optimization
fcn = @(w)MaxReturn(w,mu);

%{
[w, fval, flag, output] = ga(fcn, nAssets, A, b, Aeq, beq)

if isempty(w)
   warning('could not find any solution')
else
    %print the solution
    fprintf(2, 'Risk: %.3f%%\n', sqrt(w'*sigma*w)*100);
    fprintf(2, 'Ret: %.3f%%\n', w'*mu'*100);
end
%}

N = 100000;
x0 = rand(N,nAssets);
x0 = x0 ./ sum(x0,2);
Alhs = mu * x0';
mask = Alhs >= r_target;
x0 = x0(mask,:);
if isempty(x0)
    fprintf('failed all %d initial conditions\n', N);
    fprintf('greatest was %g but need %g\n', max(Alhs), r_target);
    fprintf('mu max is %g\n', max(mu));
else
    fprintf('Succeeded on %d initial conditions\n', size(x0,1));
    
    [wF, fvalF, flagF, outputF] = fmincon(fcn, x0(1,:), A, b, Aeq, beq)

end

failed all 100000 initial conditions

greatest was 0.00428976 but need 0.1

mu max is 0.00534691

[wF, fvalF, flagF, outputF] = fmincon(fcn, .05*ones(1,nAssets), A, b, Aeq, beq)

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.


Consider enabling the interior point method feasibility mode.

wF = 1×8

-0.000345128325127368 -0.000280137675695389 1.00268187780506 -0.000244196698363124 -0.000505350943988705 -5.42278068495475e-05 -0.000364201130047243 -0.000382595871917903

fvalF = 
       0.00535800915538925

flagF = 
    -2

outputF = struct with fields:
         iterations: 306
          funcCount: 2767
    constrviolation: 0.0946419908446108
           stepsize: 3.97580668930409e-25
          algorithm: 'interior-point'
      firstorderopt: 0.00284530672744689
       cgiterations: 343
            message: '↵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>↵↵Optimization stopped because the relative changes in all elements of x are↵less than options.StepTolerance = 1.000000e-10, but the relative maximum constraint↵violation, 9.464199e-02, exceeds options.ConstraintTolerance = 1.000000e-06.↵↵'
       bestfeasible: []

function f = MaxReturn(w,mu)
  f = w * mu';
end

[Walter Roberson]

In the below I lower the required return to 0.004 (a little smaller than the largest mu, so the target becomes do-able.)

You can see that ga() succeeds.

I also create a whole lot of random starting points for fmincon and extract ones that meet the constraints and use one of those as a starting point. You can see that if you start with random points that it is rare but possible to get a feasible starting point (and thus to be sure that fmincon will find some solution.)

We know from the experiments I posted above that fmincon with a starting point that is not feasible might not be able to find a solution.

format long g
filename = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/835520/Returns.xlsx';

data = readmatrix(filename);

nAssets = size(data, 2);

%Returns and covariance
mu = mean(data);
mu.'

ans = 8×1

0.0017001795070217 0.00238687944682926 0.00534690780710744 0.00276663664554018 7.24577440764558e-06 0.00477387349434592 0.00149865368592419 0.00130429237016939

sigma = cov(data);

%formulate the problem/optimization
r_target = 0.004; %r_target is the required return
f = zeros(nAssets, 1);    %there is no constant
A = [-mu; -eye(nAssets)] %besides the returns we forbid short selling

A = 9×8

-0.0017001795070217 -0.00238687944682926 -0.00534690780710744 -0.00276663664554018 -7.24577440764558e-06 -0.00477387349434592 -0.00149865368592419 -0.00130429237016939 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 -1

b = [-r_target; zeros(nAssets, 1)] % required return and weights greater/eqauls 0

b = 9×1

-0.004 0 0 0 0 0 0 0 0

Aeq = ones(1, nAssets)   %All weights should sum up...

Aeq = 1×8

1 1 1 1 1 1 1 1

beq = 1                  %... to one (1)

beq = 
     1

%solve the optimization
fcn = @(w)MaxReturn(w,mu);

[w, fval, flag, output] = ga(fcn, nAssets, A, b, Aeq, beq)

Optimization terminated: average change in the fitness value less than options.FunctionTolerance.

w = 1×8

0.0411595351719939 0.270953348474745 0.326503404668753 0.0742770958832299 0.12114868646004 0.0289719281976104 0.069848533983895 0.0675788103078042

fval = 
       0.00300000034134275

flag = 
     1

output = struct with fields:
      problemtype: 'linearconstraints'
         rngstate: [1×1 struct]
      generations: 61
        funccount: 11800
          message: 'Optimization terminated: average change in the fitness value less than options.FunctionTolerance.'
    maxconstraint: 0.000999999658657245
       hybridflag: []

if isempty(w)
   warning('could not find any solution')
else
    %print the solution
    fprintf(2, 'Risk: %.3f%%\n', sqrt(w*sigma*w')*100);
    fprintf(2, 'Ret: %.3f%%\n', w*mu'*100);
end

Risk: 4.686%

Ret: 0.300%

N = 100000;
x0 = rand(N,nAssets);
x0 = x0 ./ sum(x0,2);
Alhs = mu * x0';
mask = Alhs >= r_target;
x0 = x0(mask,:);
if isempty(x0)
    fprintf('failed all %d initial conditions\n', N);
    fprintf('greatest was %g but need %g\n', max(Alhs), r_target);
    fprintf('mu max is %g\n', max(mu));
else
    fprintf('Succeeded on %d initial conditions\n', size(x0,1));
    
    [wF, fvalF, flagF, outputF] = fmincon(fcn, x0(1,:), A, b, Aeq, beq)

end

Succeeded on 6 initial conditions

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

wF = 1×8

0.0557059863568776 0.0633078972843957 0.409241811954887 0.0747711685077102 0.0386402642299857 0.254060498556273 0.0532410265768362 0.0510313465330336

fvalF = 
        0.0040003437047322

flagF = 
     1

outputF = struct with fields:
         iterations: 18
          funcCount: 171
    constrviolation: 0
           stepsize: 0.00179880774756221
          algorithm: 'interior-point'
      firstorderopt: 3.43704732204755e-07
       cgiterations: 0
            message: '↵Local minimum found that satisfies the constraints.↵↵Optimization completed because the objective function is non-decreasing in ↵feasible directions, to within the value of the optimality tolerance,↵and constraints are satisfied to within the value of the constraint tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 3.437047e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.↵↵'
       bestfeasible: [1×1 struct]

[wF, fvalF, flagF, outputF] = fmincon(fcn, .05*ones(1,nAssets), A, b, Aeq, beq)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

wF = 1×8

0.0564935222995476 0.0687639857461034 0.442177416664919 0.0730119927010965 0.0402106863233072 0.215233398172256 0.0533706189672584 0.0507383791255117

fvalF = 
       0.00400041013681353

flagF = 
     1

outputF = struct with fields:
         iterations: 18
          funcCount: 173
    constrviolation: 0
           stepsize: 0.000479311697708811
          algorithm: 'interior-point'
      firstorderopt: 4.10136813534627e-07
       cgiterations: 0
            message: '↵Local minimum found that satisfies the constraints.↵↵Optimization completed because the objective function is non-decreasing in ↵feasible directions, to within the value of the optimality tolerance,↵and constraints are satisfied to within the value of the constraint tolerance.↵↵<stopping criteria details>↵↵Optimization completed: The relative first-order optimality measure, 4.101368e-07,↵is less than options.OptimalityTolerance = 1.000000e-06, and the relative maximum constraint↵violation, 0.000000e+00, is less than options.ConstraintTolerance = 1.000000e-06.↵↵'
       bestfeasible: [1×1 struct]

function f = MaxReturn(w,mu)
  f = w * mu';
end