Constraint optimization problem, fmincon, converged to an infeasible point.-?

Hi

I'm having trouble with an exercise. I don't know how to use the starting guess value x0. When I run it I get this:

_Converged to an infeasible point.

fmincon stopped because the predicted change in the objective function is less than the default value of the function tolerance but constraints are not satisfied to within the default value of the constraint tolerance.

< stopping criteria details>_

Although I get final values I am not sure they are right because of this message I get. These are the values: yopt=6.1011; zopt=0.1098; fval=3.9996.

First I have the cost function, then the constraints function, and the script file to obtain the optimized value:

function [ J ] = objfun( x )

%Objective function

P=25.0231; E=2; h=2;

J=((P*h)/E)*(1/(x(1)+sqrt(2)*x(2)));

end

function [ c,ceq ]=confun( x )

%Non-linear constraints

%Constraint functions

S=0.0092; P=25.0231;

c= [(P*(x(2)+x(1)*sqrt(2))/(x(1)^2*sqrt(2)+2*x(1)*x(2)))-17.5;

      (P/(x(1)+x(2)*sqrt(2)))-17.5*(1+0.1*S);

      (P/x(1)^2*sqrt(2)+2*x(1)*x(2))-1.2;

      0.2-x(1);

      x(1)-6*(1+0.2*S);

      0.2*(1-0.1*S)-x(2);

      x(2)-6];

ceq=[];

end

S=0.0092
P=25*(1+0.1*S);

x0=[3 3]; %Starting guess

options=optimset('Algorithm','active-set');

[x,fval]=fmincon(@objfun,x0,[],[],[],[],[],[],@confun,options);

%Solution x and function value fval

yopt=x(1);zopt=x(2);

P.s. I would also like to know how to use the command Global in this exercise, for P and S, as they repeat in both function files and script file; and it looks much better.

[Walter Roberson]

Do not use "global" for that purpose; see http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[George Schwarzenneger]

Thank you

Answer 1

[Walter Roberson]

Your last 4 constraints, [0.2-x(1); x(1)-6*(1+0.2*S); 0.2*(1-0.1*S)-x(2); x(2)-6]

are all linear constraints of one variable each. You should convert them to upperbound and lowerbound constraints instead, which will make the search much more efficient.

A = []; b = []; 
S=0.0092;
lb = [0.2, 0.2*(1-0.1*S)];
ub = [6*(1+0.2*S), 6];
[x, fval] = fmincon(@objfun, x0, A, b, Aeq, beq, lb, ub, @confun, options);

Now look at your second constraint, (P/(x(1)+x(2)*sqrt(2)))-17.5*(1+0.1*S)) . With a brief examination you can see that since x(1) and x(2) are both positive, x(1)+x(2)*sqrt(2) will never be 0, so you can rewrite this as a linear constraint

(-(0.1) * (17.5) * S - 17.5) * x1 + (-sqrt(2) * 0.1 * 17.5 * S - sqrt(2) * 17.5) * x2 + P <= 0

which can be coded in the linear constraints as

A = [-(0.1) * (17.5) * S - 17.5), (-sqrt(2) * 0.1 * 17.5 * S - sqrt(2) * 17.5)];
b = [-P];

and you are now down to two non-linear constraints.

Putting these together,

E=2;
h=2;

S = 0.0092;
P = 25*(1+0.1*S);

A = [-(0.1) * (17.5) * S - 17.5), (-sqrt(2) * 0.1 * 17.5 * S - sqrt(2) * 17.5)];
b = [-P];
Aeq = []; beq = [];
lb = [0.2, 0.2*(1-0.1*S)];
ub = [6*(1+0.2*S), 6];

x0 = [3 3]; %Starting guess

options = optimset('Algorithm','active-set');

objfun = @(x) ((P*h)/E)*(1/(x(1)+sqrt(2)*x(2)));

[x, fval] = fmincon(objfun, x0, A, b, Aeq, beq, lb, ub, @(x) confun(x, S, P), options);

%Solution x and function value fval

yopt=x(1); zopt=x(2);

with no separate file for objfun, and with

function [c, ceq] = confun(x, S, P)
%Non-linear constraints
%all the other constraints moved to bounds constraints or linear constraints

    c = [(P * (x(2)+x(1)*sqrt(2)) / (x(1)^2*sqrt(2) + 2*x(1)*x(2))) - 17.5;
         (P/x(1)^2*sqrt(2) + 2*x(1)*x(2)) - 1.2];

    ceq=[];

end

I make no promises that this will allow the equation to be solved. I think I have found there is no solution with those constraints, but the proof of that will wait for another posting.

[Walter Roberson]

Your third constraint together with your bounds constraint rule out any possible solutions.

Your 4th constraint 0.2-x(1) forces x(1) to be positive (0.2 exactly is the minimum)

Your 6th constraint 0.2*(1-0.1*S)-x(2) <= 0 forces x(2) to be positive (just under 0.2 is the minimum.)

Your third constraint is (P/x(1)^2*sqrt(2)+2*x(1)*x(2))-1.2 <= 0 . Rewrite to make it easier to read:

P/x1^2 * sqrt(2) + 2*x1*x2 - 6/5 <= 0

x2 occurs only in the +2*x1*x2 term. We know that x1 and x2 are positive from the above analysis, so +2*x1*x2 is contributing a positive value, making the <= 0 harder to satisfy. The constraint will be easiest to satisfy when x2 is smallest, 0.2*(1-0.1*S), or 24977/125000 expressed as a fraction. This, together with the rational equivalent of P, leads to

(25023/1000) * sqrt(2)/ x1^2 + (24977/62500)*x1 - 6/5 <= 0

Examining, when x1 is very small, the first term is large and positive, certainly greater than 6/5, so near 0 the left expression is non-negative. When x1 approaches infinity, the first term is very small but the second term is large, certainly greater than 6/5, so near infinity the left expression is non-negative. So there might not be any roots . Are there? Differentiate the left side with respect to x1, solve for 0; it is a cubic so there are 3 roots, two of which are imaginary. The real root is at (5/24977)*15610611787167^(1/3)*2^(1/6) which is about x1 = 5.616, so the only real-valued critical point is near x1 = 5.616 . And at that x1, the left expression has a value of (3/25000)*15610611787167^(1/3)*2^(1/6)-6/5 which is about 2.1664 and is a minima. Therefore there are no positive x1 for which (25023/1000) * sqrt(2)/ x1^2 + (24977/62500)*x1 - 6/5 <= 0, and that was the third constraint with your x2 set to make it easiest to satisfy.

Therefore your third constraint cannot be satisfied in the presence of your 4th constraint and 6th constraint (positive x(1) and positive x(2))

[George Schwarzenneger]

Thank you very much for your time!

Answer 2

[George Schwarzenneger]

Thank you for your time. I'm gonna try that way too. Anyways, I made some changes and I think it works as expected:

1st: in the objective function (objfun) I changed the inputs, adding P, h and E. I did the same for the constraints function (confun), including P and S in the input.

2nd: I rewrote the constraints equations, adding more brackets to make it more clear.

3rd: x0 has been changed. from [3 3] to [-2 2].

Here is what I have done:

function [ J ] = objfun(x,P,h,E)

%Objective function

J=((P*h)/E)*(1/(x(1)+sqrt(2)*x(2)));

end

%

function [ c,ceq ]=confun(x,P,S)

%Non-linear constraints

%Constraint functions

c=[P*(((x(2))+(sqrt(2))*(x(1)))/(((sqrt(2))*(x(1)).^2)+(2*(x(1))*(x(2)))))-17.5; P*(((x(1))+(sqrt(2)*x(2))).^-1)-17.5*(1+0.1*S); P*((x(2))/(((sqrt(2))*(x(1)).^2)+(2*(x(1))*(x(2)))))-1.2; 0.2-(x(1)); (x(1))-6*(1+0.2*S); 0.2*(1-0.1*S)-(x(2)); (x(2))-6];

ceq=[];

end

%

S=0.0092;%'S' parameter

P=25*(1+0.1*S);E=2;h=2;

%

x0=[-2 2]; %Starting guess

options=optimset('Algorithm','active-set');

[x,fval]=fmincon(@(x)objfun(x,P,h,E),x0,[],[],[],[],[],[],@(x)confun(x,P,S),options);

%Solution x and function value fval

yopt=x(1);

zopt=x(2);

fprintf('yopt is %5.2f \n',x(1));

fprintf('zopt is %5.2f \n',x(2));

fprintf('fmin is %5.2f \n',fval);

%

And this is the result:

Local minimum found that satisfies the constraints.

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

< stopping criteria details>

Active inequalities (to within options.TolCon = 1e-06): lower upper ineqlin ineqnonlin 3 5

yopt is 6.01

zopt is 5.79

fmin is 1.76