Optimization of an implicit function using fmincon

I am trying to carry out multi-objective constrained optimization using fmincon in MATLAB. However, one of my objective function can not be brought in an explicit form depending on variables. For example, I have an equation like:

x^a + x^b = y

and I want to minimize x given a range of values for a, b and y respectively.

How can I do that?

Answer 1

[Walter Roberson]

You cannot do multi-objective optimization using fmincon(). fmincon() is restricted to single objectives.

Which of a, b, x, y are to be varied during a single call to fmincon() ?

If a, b, y are constant for this call, then use nonlinear equality

ceq = x.^a + x.^b - y

[Harshit Agrawal]

Let us take single objective only.

Even in that case, since I can not express x in terms of a,b and y, what should I take as my objective function?

a, b and y are variables. The constraints are for example: 0<a<=100; 0<b<=59; 0<y<=675

[Walter Roberson]

If you can constrain x to be positive and real (mathematically, there are negative and complex x that would solve the equation) then you can calculate upper and lower bounds on x based upon the equation. Then add it as an additional optimization variable, using those upper and lower bounds, and adding ceq = x.^a + x.^b - y as a nonlinear equality. Be sure to use an initial value that satisfies the constraint or MATLAB can spend a long time trying to find something that works.

[Harshit Agrawal]

Thanks for your answer.

Are you suggesting to define objective function as 'x' and keep an equality constraint as x.^a + x.^b - y = 0 along with other inqeuality constraints for a, b and y for their permissible ranges?

[Walter Roberson]

That could work.

[Harshit Agrawal]

That generally works but result varies with varying initial point. Could you please suggest a robust way to find global minima?

[Walter Roberson]

Numeric analysis

x = eps(realmin) 
a = 1
b = 1
y = 2*x

The only smaller non-negative x possible in floating point is x = 0. But 0 to a power is 0 unless the power is 0, but a and b are not permitted to be exactly 0. So x = 0, the left would be 0+0=0. However, y = 0 is not permitted.

We have thus proven that 0 exactly is not permitted. So let x be the smallest nonzero value. Can we make that work? Yes, as indicated above.

You can do better on y though. Let x = eps(realmin) and a = 1 and b>= 1+eps. Then x^b < x because b>1 and x<1 and since there is no smaller representable number, x^b underflows to 0. x^1 = x so x^a stays eps(realmin) and we add the 0 from x^b to that, getting eps(realmin) for y, which is within the permitted range. This is the smallest possible y since y=0 is not permitted and this is the next representable number.

[Harshit Agrawal]

Thanks. Numerical analysis does work in the case of this generalised problem. However, the problem I am trying to solve is a more complex one.

I am trying to solve following equation for minimum D_shell.

Q - (10536339*D_ti*H_coil*pi*((2778046668940015*D_coil^2)/281474976710656 + p^2)^(1/2))/(50000000000000*p*((D_to*log(D_to/D_ti))/90000 + (180143985094819840*D_ti*(1000/D_ti)^(9/200)*(D_coil/D_ti)^(27/125)*(D_shell/D_ti)^(3/250)*(H_coil/D_ti)^(3/100))/(2895380937030248271*(D_v/D_ti)^(3/125)*(((D_v - 1000)^2 + D_shell^2)^(1/2)/D_ti)^(13/1000)*(p/D_ti)^(11/1000)*(4702891291781029703125000/(4329144886923138451*D_ti*pi))^(137/200)) + (2251799813685248000*D_ti*(1000/D_to)^(103/100)*(D_shell/D_to)^(41/50)*(1000*pi*(D_shell^2 - 222784) - (pi*D_to^2*H_coil*((2778046668940015*D_coil^2)/281474976710656 + p^2)^(1/2))/p))/(486758140305134317113*D_to*pi*(D_coil/D_to)^(189/500)*(D_v/D_to)^(139/250)*(H_coil/D_to)^(43/1000)*(((D_v - 1000)^2 + D_shell^2)^(1/2)/D_to)^(561/1000)*(p/D_to)^(69/500)*((H_coil*((2778046668940015*D_coil^2)/281474976710656 + p^2)^(1/2))/p + 1000)*(D_to + D_shell + 472)*(-(7576000000*D_to)/(723*pi*((D_coil + D_to)^2 - D_shell^2 - (D_coil - D_to)^2 + 222784)))^(723/1000)))) == 0

Q is a constant.

D_ti,D_to,H_shell,H_coil,D_coil,p,D_v are all variables with a range of positive values.

I have tried many algorithms like fmincon, fminsearch but nothing seems to work for global minima.

Please help.

[Walter Roberson]

Numeric algorithms can never promise global minima, even given a really good starting point. Too much round-off.

[Walter Roberson]

Are there upper limits for the variables?

[Harshit Agrawal]

Yes

[Harshit Agrawal]

Any idea on how to solve the problem?

[Walter Roberson]

You will need to do numeric analysis. It will probably take you a number of days to carry out, possibly weeks. The expression is not very tractable, so it will be difficult.

Remember that the only way to prove global minima is by theoretical analysis. The best that you can do with fmincon or related tools is to find a local minima, and it might even be a great looking local minima, but numeric tools such as fmincon cannot rule out the possibility that the expression has a narrow unstable global minima (balancing a pin on its point sort of thing.)

[Harshit Agrawal]

Will geneteic algorithm help in finding globala minima?

[Walter Roberson]

Genetic algorithms explore more of the space, and can help get out of ruts. However:

• they never know (and cannot know) whether they have found a global minima. All they can ever know is that they have stalled finding better points -- but that could just mean that more iterations are needed
• they are not at all efficient at finding the local minima within whichever catch-basin they are in. They might be very close to a great minima and never locate it
• they can get misled and waste a lot of time. For example if you have a asymtopic descent around a central dip, they might spend their time chasing the ever-lower value further and further from the center, only ending up in the center through chance mutations

Mathworks claims that patternsearch() is typically more efficient than ga()

[Walter Roberson]

Saswath Ghosh comments to @Harshit Agrawal

Did this approach work? to pose implicit function as objective