Non Linear optimization problem in MATLAB

I have following problem which I am using MATLAB to solve it seems to me pertty Please check Is this correct or NOT?

Question

objective function (objfun2.m)

function f = objfun2(x)
    f = 7.9 * 0.000000001 * 300 * pi * (x(1)^2 - x(2)^2);
end

Non linear constraint (confun2.m)

function [c, ceq] = confun2(x)
%Nonlinear inequality constraints
 sqrt1=sqrt((1.84*100000* (pi/4)*(x(1)^4-x(2)^4))/(7.9*0.000000001*pi*(x(1)^2-x(2)^2)));
 sqrt2= sqrt(1+(x(3)*300^3)/(3 * 1.84 * 100000 *(pi/4)*(x(1)^4-x(2)^4)));
 c = -(((1.875^2)/(300^2 * 2*pi))*sqrt1*sqrt2)+230 ; % since it should have x<m form
% Nonlinear equality constraints

ceq = [];

Main Code

LB=[9.2,6,0.1];
UB=[15,9.8,18];
A=[];
b=[];
Aeq=[];
beq=[];

x0=[10,9,0.15];

% Make a starting guess at the solution
%options = optimoptions(@fmincon,'Algorithm','sqp');
[x,fval] = fmincon(@objfun2,x0,A,b,Aeq,beq,LB,UB,@confun2);

My Solution (Not sure it is correct or NOT)

 x= 9.7603    9.7603   16.8680

How would I check this?

Answer 1

[John D'Errico]

As Matt points out, the solution you found may have issues.

You can see the source of the problem by looking at your nonlinear constraint. There we see you have the difference (x1^2 - x2^2) in the denominator. So when x1 is approximately equal to x2, things go to hell.

You can fix part of that, by factoring the terms in the first radical.

There we see a fraction with

(x1^4 - x2^4)/(x1^2 - x2^2)

this reduces to (x1^2 + x2^2), so no 0/0 there.

The second radical does have a divide by zero, when x1==x2. So again, this will be a problem. You will be able to get almost any value you want for that term, by TINY changes to x1 and x2.

Did the optimizer converge? Hard to tell. I'start the optimization at a different location though, since you still have a singular constraint.

If your goal is to minimize (x1^2-x2^2) you can do so by maximizing x2, as well as minimize x1. (Drop the constants out front! A minimization does not care about them, and at best, they may just introduce numerical problems.)

So your start point should be roughly:

[9.2 9.8 , ?]

I put a question mark in there for x3, because x3 does not appear in the objective. So we can choose a start point for x3 by trying to satisfy the nonlinear constraint.

So now lets look more closely at the constraint. We know that x3 lies in [0.1, 18].

First, lets simplify the constraint.

K1 = 1.875^2/(300^2*2*pi)*sqrt(1.84e5/4/7.9e-9)
K1 =
        15.0018747945489

K2 = 300^3/(3*1.84e5*pi/4)
K2 =
        62.2780212098721

So we can write the entire nonlinear constraint in the form of a function handle:

CON = @(x123) K1*sqrt(x123(1)^2 + x123(2)^2)*sqrt(1 + K2*x123(3)/(x123(1)^4 - x123(2)^4)) - 230;

Will this function ever be greater than zero, when x1 and x2 are at their best locations?

CON([9.2 9.8 .1])
ans =
        -28.653992101973

CON([9.2 9.8 18])
ans =
       -93.8657086630407

ezplot(@(x3) CON([9.2 9.8 x3]),[.1 18])

So at what I chose as the best point for x1 and x2, no value of x3 will allow this constraint to be positive. But if I choose other values for x1 and x2, as long as x1>x2, there are solutions that allow con to be positive.

ezplot(@(x3) CON([9.8 9.7 x3]),[.1 18])
grid on

Things look very different when x1 is slightly less than x2 however.

ezplot(@(x3) CON([9.6 9.7 x3]),[.1 18])
grid on

In fact, it appears there may be no solution at all when x1 is less than x2. I'd need to look more assiduously at the problem to know if that is true.

But it is clearly true that [9.8 9.7 2] is a feasible point. So what does fmincon find?

Define CONFUN as:

function [c, ceq]  = confun(x123)
K1 = 1.875^2/(300^2*2*pi)*sqrt(1.84e5/4/7.9e-9);
K2 = 300^3/(3*1.84e5*pi/4);
c = K1*sqrt(x123(1)^2+x123(2)^2)*sqrt(1+K2*x123(3)/(x123(1)^4-x123(2)^4))-230;
c = -c;
ceq = [];

I negated CON in there.

LB=[9.2,6,0.1];
UB=[15,9.8,18];
A=[];
b=[];
Aeq=[];
beq=[];
x0=[9.8,9.7,2];

[x,fval,exitflag] = fmincon(@(x123) x123(1)^2 - x123(2)^2,x0,A,b,Aeq,beq,LB,UB,@CONFUN)
Local minimum possible. Constraints satisfied.

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

<stopping criteria details>

x =
          9.53501396574518          9.53501396574518          2.80228092041589
fval =
       2.8421709430404e-14
exitflag =
     2

Fmincon wants to push to the boundary where x1 == x2. The constraint is satisfied at x.

CON(x)
ans =
        1144008644.36073

But, clearly, the solution is a bit of a singular one, since it lies right along the edge where x1==x2, so we have an effective divide by zero.

[Naseer Khan]

Really really thanx for such hardcore explanations........and thanks for your so much time for this problem...........regards

[John D'Errico]

I almost forgot to point this out:

At the minimizer, x1==x2. So your objective is zero. But there are infinitely many points where x1==x2. As long as the constraint is satisfied, you will have a solution. That tends to imply the solution is non-unique.

So now lets look at the locus of points: [x1,x1,x3]. Thus we will set x1==x2.

x(1) - x(2)
ans =
      1.77635683940025e-15

The difference was as small as fmincon could make it. But there will always be infinitely many points that lie essentially along that line with a zero objective value, and that satisfy the constraint. So your problem does not have a unique solution.

Answer 2

[Matt J]

You should call fmincon with more output arguments, so that you can assess the quality of the solution

[x,fval,exitflag,output] = fmincon(___)

However, I would not trust the solution that you've shown above. It looks like fmincon was trying to reach a solution where x1=x2, which is outside the domain of your nonlinear constraint function.

Was your nonlinear constraint active at the solution? We cannot test that because you've only shown the solution to 4 decimal places. If all the constraints were inactive at the final x, then this is definitely not a solution, because the gradient is non-zero there (so the KKT conditions are not satisfied).

You might try re-solving the problem with the demoninators in your nonlinear constraints cleared, i.e., instead of

   num/den-230>=0

express the constraints as

   num-230*den>=0

If you reach a solution where the objective function is significantly less than 0 (and the constraints are satisfied), then you'll know that solution is better.

You could also try evaluating your objective function and constraints on a coarsely sampled grid covering the area within your LB, UB bounds. This would be a way to obtain a better initial guess and a better idea of what the globally minimum value is. Because there are only 3 variables, this should be computationally tractable.

[John D'Errico]

+1