Fmincon converged to an infeasible point, I don't know why

Hi guys,

I'm trying to minimize a function using fmincon, the program runs correctly but fmincon doesn't give me the answer I'd like. I don't know why fmincon stops and gives me the message 'Converged to an infeasible point'.

Here you can find my Matlab code I use.

Thank you for your help, I'm a beginner in coding.

% Optimization model

fun = @(d) parameterfun(d,L_mean,L_std);
L_mean = 2000*10^3;
L_std = 200*10^3;

% Component 1
c1_mean = 3500*10^3;
c1_std = 350*10^3;
pf1 = 9.920*10^-5;
R1 = 1 - pf1;

% Component 2
c2_mean = 3200*10^3;
c2_std = 260*10^3;
pf2 = 1.2696*10^-4;
R2 = 1-pf2;

% Component 3
c3_mean = 4000*10^3;
c3_std = 400*10^3;
pf3 = 3.87*10^-6;
R3 = 1-pf3;

S1_mean = L_mean+norminv(R1)*sqrt(c1_std^2+L_std^2);
S2_mean = L_mean+norminv(R2)*sqrt(c2_std^2+L_std^2);
S3_mean = L_mean+norminv(R3)*sqrt(c3_std^2+L_std^2);

L_mean = [S1_mean; S2_mean; S3_mean];
L_std = [c1_std^2 L_std^2 L_std^2; L_std^2 c2_std^2 L_std^2; L_std^2 L_std^2 c3_std^2];

lb = [0,0];
ub = [1,1];
A = [];
B = [];
Aeq = [];
Beq = [];
d0 = (lb + ub)/2;
const = @(d) nonlcon(d,L_mean,L_std);
options = optimoptions('fmincon','Display','iter');
[d,fval] = fmincon(fun,d0,A,B,Aeq,Beq,lb,ub,const,options);

function Rs = parameterfun(d,L_mean,L_std)

% System load
L_mean = 2000*10^3;
L_std = 200*10^3;

% Component 1
c1_mean = 3500*10^3;
c1_std = 350*10^3;
pf1 = 9.920*10^-5;
R1 = 1 - pf1;

% Component 2
c2_mean = 3200*10^3;
c2_std = 260*10^3;
pf2 = 1.2696*10^-4;
R2 = 1-pf2;

% Component 3
c3_mean = 4000*10^3;
c3_std = 400*10^3;
pf3 = 3.87*10^-6;
R3 = 1-pf3;

S1_mean = L_mean+norminv(R1)*sqrt(c1_std^2+L_std^2);
S2_mean = L_mean+norminv(R2)*sqrt(c2_std^2+L_std^2);
S3_mean = L_mean+norminv(R3)*sqrt(c3_std^2+L_std^2);

L_mean = [S1_mean; S2_mean; S3_mean];
L_std = [c1_std^2 L_std^2 L_std^2; L_std^2 c2_std^2 L_std^2; L_std^2 L_std^2 c3_std^2];

d = zeros(size(L_mean));

Rs = mvncdf(d, L_mean, L_std);

end


function [c, ceq] = nonlcon(d,L_mean,L_std)

% System load
L_mean = 2000*10^3;
L_std = 200*10^3;

% Component 1
c1_mean = 3500*10^3;
c1_std = 350*10^3;
ns1_min = 1.5;
ns1_max = 2.5;
c1_min = 0.08;
c1_max = 0.2;
pf1 = 9.920*10^-5;
R1 = 1 - pf1;

% Component 2
c2_mean = 3200*10^3;
c2_std = 260*10^3;
ns2_min = 1.5;
ns2_max = 2.5;
c2_min = 0.08;
c2_max = 0.2;
pf2 = 1.2696*10^-4;
R2 = 1-pf2;

% Component 3
c3_mean = 4000*10^3;
c3_std = 400*10^3;
ns3_min = 1.5;
ns3_max = 2.5;
c3_min = 0.08;
c3_max = 0.2;
pf3 = 3.87*10^-6;
R3 = 1-pf3;

% Define inequality constraints which use parameters L_mean and L_std
c(1) = ns1_min - (L_mean+norminv(R1)*sqrt(c1_std^2+L_std^2)/L_mean);
c(2) = ns2_min - (L_mean+norminv(R2)*sqrt(c2_std^2+L_std^2)/L_mean);
c(3) = ns3_min - (L_mean+norminv(R3)*sqrt(c3_std^2+L_std^2)/L_mean);
c(4) = (L_mean+norminv(R1)*sqrt(c1_std^2+L_std^2)/L_mean) - ns1_max;
c(5) = (L_mean+norminv(R2)*sqrt(c2_std^2+L_std^2)/L_mean) - ns2_max;
c(6) = (L_mean+norminv(R3)*sqrt(c3_std^2+L_std^2)/L_mean) - ns3_max;
c(7) = c1_min - (c1_std/(L_mean+norminv(R1)*sqrt(c1_std^2+L_std^2)));
c(8) = c2_min - (c2_std/(L_mean+norminv(R2)*sqrt(c2_std^2+L_std^2)));
c(9) = c3_min - (c3_std/(L_mean+norminv(R3)*sqrt(c3_std^2+L_std^2)));
c(10) = (c1_std/(L_mean+norminv(R1)*sqrt(c1_std^2+L_std^2))) - c1_max;
c(11) = (c2_std/(L_mean+norminv(R2)*sqrt(c2_std^2+L_std^2))) - c2_max;
c(12) = (c3_std/(L_mean+norminv(R3)*sqrt(c3_std^2+L_std^2))) - c3_max;

% Define equlaity constraints
ceq = [];

end

[Torsten]

You don't use the parameters "d" from "fmincon" in your calculations - neither in "parameterfun" nor in "nonlcon". So how can you expect that the optimizer could converge to something that makes sense if you repeatedly return wrong information ?

[John D'Errico]

As I show in my answer, not only are the constraints not a function of d, some constraints return positive numbers, some are negative. And that means the solution will NEVER return a feasible result.

Answer 1

[John D'Errico]

I'm sorry, but this code is a bit strange. Ok, wildly strange.

Why are you using fmincon? You have no linear inequality constratins, you have no linear inequality constraints. So then I looked at nonlcon.

Just for kicks, look at the nonlinear ineqality constraints. You pass in L_mean, and L_std, but then you immediately overwrite them. I think you need what functions do and how to use them, but that is just a general comment. But now look carefully at nonlcon.

nonlcon constrains the parameters. You hve TWO unknown parameters, the elements of d. Are ANY of the values of the output of nonlcon a function of d? NO. Let me repeat that. NO. So your constraints are completely independent of the problem parameters.

Next, just for kicks, I evaluated the nonlinear constraints you created.

I'll create the vector d, as NaNs. This way, in case ANYTHING uses d that I did not notice, we will see a NaN in the resulting vector c.

d = [nan, nan];

L_mean = 2000*10^3;
L_std = 200*10^3;
% System load
L_mean = 2000*10^3;
L_std = 200*10^3;
% Component 1
c1_mean = 3500*10^3;
c1_std = 350*10^3;
ns1_min = 1.5;
ns1_max = 2.5;
c1_min = 0.08;
c1_max = 0.2;
pf1 = 9.920*10^-5;
R1 = 1 - pf1;
% Component 2
c2_mean = 3200*10^3;
c2_std = 260*10^3;
ns2_min = 1.5;
ns2_max = 2.5;
c2_min = 0.08;
c2_max = 0.2;
pf2 = 1.2696*10^-4;
R2 = 1-pf2;
% Component 3
c3_mean = 4000*10^3;
c3_std = 400*10^3;
ns3_min = 1.5;
ns3_max = 2.5;
c3_min = 0.08;
c3_max = 0.2;
pf3 = 3.87*10^-6;
R3 = 1-pf3;
% Define inequality constraints which use parameters L_mean and L_std
c(1) = ns1_min - (L_mean+norminv(R1)*sqrt(c1_std^2+L_std^2)/L_mean);
c(2) = ns2_min - (L_mean+norminv(R2)*sqrt(c2_std^2+L_std^2)/L_mean);
c(3) = ns3_min - (L_mean+norminv(R3)*sqrt(c3_std^2+L_std^2)/L_mean);
c(4) = (L_mean+norminv(R1)*sqrt(c1_std^2+L_std^2)/L_mean) - ns1_max;
c(5) = (L_mean+norminv(R2)*sqrt(c2_std^2+L_std^2)/L_mean) - ns2_max;
c(6) = (L_mean+norminv(R3)*sqrt(c3_std^2+L_std^2)/L_mean) - ns3_max;
c(7) = c1_min - (c1_std/(L_mean+norminv(R1)*sqrt(c1_std^2+L_std^2)));
c(8) = c2_min - (c2_std/(L_mean+norminv(R2)*sqrt(c2_std^2+L_std^2)));
c(9) = c3_min - (c3_std/(L_mean+norminv(R3)*sqrt(c3_std^2+L_std^2)));
c(10) = (c1_std/(L_mean+norminv(R1)*sqrt(c1_std^2+L_std^2))) - c1_max;
c(11) = (c2_std/(L_mean+norminv(R2)*sqrt(c2_std^2+L_std^2))) - c2_max;
c(12) = (c3_std/(L_mean+norminv(R3)*sqrt(c3_std^2+L_std^2))) - c3_max;

Now what is the return vector c?

format long
c

c = 1×12

1.0e+06 * -1.999999250000617 -1.999999100001466 -1.999999500026034 1.999998250000617 1.999998100001466 1.999998500026035 -0.000000019999965 -0.000000001249926 -0.000000019998698 -0.000000100000035 -0.000000118750074 -0.000000100001302

Soem elements of c are positive, some are negative. None of them are NaNs. your constraints are FIXED, and completely independent of the parameters in the vector d, thus the unknowns.

What does fmincon do with the constraints? It tries to find a solution that satisfies all of the constraints, and still minimizes the objective. BUT YOUR CONSTRAINTS ARE INDEPENDENT OF d. Nothing that fmincon will ever do will satisfy the constraints.

Now go back and read the message you got:

'Converged to an infeasible point'.

Do you understand why it gives you that message? A feasible solution is one that satisfies all of the constraints. For it to be feasible, EVERY element of c must be <= 0, to within a tolerance. Your constraints are constant values, with some positive, some negative. What did you expect fmincon to tell you?

[Paul AGAMENNONE]

Hello John, thank you for your constructive answer. I haven't found the solution yet but I noticed I made a big mistakes when writing my inequality constraints. I still have some to learn for writing my code but your answer help me a lot.

Thank you

[John D'Errico]

If it helped you to solve this problem, then please click on the accept button, so that someone else need not try to solve it again. While I might guess there still may be issues for you once you change the constraints to something valid, this specific question should now be marked as resolved.