error in objective function of fmincon and solver stopped prematurely

this is my coding

function [x,fval,flag] = minim

    x0 = [0.5, 0.5, -0.5]; 
    lb = [-1, -1, -1]; 
    ub = [1, 1, 1]; 

    visc = 0.001;
    dens = 1000;
    w = 0.5;
opts = optimset('Display','iter','Algorithm','interior-point');
[x,fval,flag] = fmincon(@phenol,x0,[],[],[],[],lb,ub,@constraints,opts);

 function J = phenol(x)

 % x1 = flowrate, x2 = voidage, x3 = particles diameter
 f1 = (0.778 + 0.0433*(x(1)) + 0.0196*(x(2)) - 0.1111*(x(3))); %efficiency
 f2 = ((150*visc*((1-(x(2)))^2)*(x(1))/(((x(2))^3)*((x(3))^2))) + (1.75*(1-(x(2)))*dens*((x(1))^2)/(((x(2))^3)*(x(3))))); %pressure drop using Ergun equation

 J = w*f1 + (w-1)*f2; %weightage
 end % end of phenol

 function [c, ceq] = constraints (x)
        Pdrop = ((150*visc*((1-(x(2)))^2)*(x(1))/(((x(2))^3)*((x(3))^2))) + (1.75*(1-          (x(2)))*dens*((x(1))^2)/(((x(2))^3)*(x(3)))));
        c = Pdrop - 58.9; %pressure must < 58.9 bar
        ceq = []; % no values assigned for linear equalities
 end % end of constraints
end  % end of minim

when i tried run the data it iterate but end with a message saying

"Solver stopped prematurely. fmincon stopped because it exceeded the function evaluation limit, options.MaxFunEvals = 3000 (the default value)"

and before this when the x0 = [-0.5, 0, 0] it says error in objective function. please kindly help :)

Answer 1

[Star Strider]

To solve the ‘MaxFunEvals’ problem, change your optimset arguments to:

opts = optimset('Display','iter','Algorithm','interior-point', 'MaxIter', 10000, 'MaxFunEvals', 10000);

I didn’t run your code, but since x(2) (and x(3) if I’ reading your code correctly) appear in the denominator of some of your statements in ‘phenol’, they cannot allowed to be zero. Instead of: x0=[-5,0,0] change the zeros to small numbers, for example 1E-10, so:

x0 = [-5, 1E-10, 1E-10];

should solve that problem.

[arina ariff]

I've tried on changing the

    x0 = [-5, 1E-10, 1E-10];

and it said that the values assigned is exceeding the lower boundary value. I've manage to change it to value within the boundary limit and the coding manage to run but the message of "solver stopped prematurely" still appears in the end of the run. I've received some guideline saying that i need to use real value rather than the coded one. can u help?

[Star Strider]

I’ll do my best!

I didn’t look for the bounds you set on your optimisation.

First what do you mean by ‘real value rather than the coded one’? That doesn’t make sense to me. PLEASE include the entire error message (red type in the Command Window), with line numbers, function or script references, and everything else.

[arina ariff]

here

    Solver stopped prematurely.

    fmincon stopped because it exceeded the function evaluation limit,
    options.MaxFunEvals = 10000 (the selected value).

    ans =

       -0.1689    0.6603    1.0000

the x(1), x(2), and x(3) from the previous coding refers to flow rate, voidage and particle diameter respectively.

from what i have retrieved from my lecturer (he's on leave for now so disturbing him is a big no no..hehe). he said to i need to establish real values for the flowrate, voidage and particle diameter as the pressure drop is in real value. so i must change for example

voidage = baseline_voidage + x(??)*dvoidage; dvoidage is the magnitude of voidage change.

all of the coding is posted in the previous post. please help me :(

btw i've change

    x0 = [0.5, 0.5, -0.5]

[Star Strider]

I can help you with the MATLAB code as much as my expertise permits, but I don’t recognise what you are doing (chemical engineering? civil engineering?) and have no expertise in either of those disciplines.

One thing I note is that of your parameter limits, all of them are allowed to be negative. I assume a negative flow rate is reverse flow, I have no idea what a negative voidage implies, but a negative particle size may be a serious problem. (It brings to mind a sort of Schrödinger’s particle that is both there and not there at the same time!) Check those limits to be sure they’re realistic in physical terms and in the context of the system you are modeling.

I suggest you ask someone else in your department to be sure you have your equations coded correctly and your parameter bounds set appropriately. There may be contributors here on MATLAB Answers with the requisite expertise in your area to look over your code and find whatever errors you may have made, but I am not one.

There may also be other optimisation routines that are better suited to your particular problem than fmincon. You might also investigate those possibilities.

[Matt J]

One thing I note is that of your parameter limits, all of them are allowed to be negative. I assume a negative flow rate is reverse flow, I have no idea what a negative voidage implies, but a negative particle size may be a serious problem.

In addition to voidage<0 and diameter<0 being non-physical, allowing x(2)<0 or x(3)<0 means that the feasible region is a subset of x(2)*x(3)~=0 and is split into 4 disjoint quadrants. fmincon is designed to search connected feasible sets only.

[Star Strider]

Thanks, Matt.

I didn’t pick up on the effects of poorly-defined constraints on fmincon performance.

[Matt J]

I get successful termination with the following. It wasn't enough just to modify x0 and the lower bounds. I also had to set MaxFunEvals and MaxIter to Inf. Apparently, it's just a poorly conditioned problem with slow asymptotic convergence.

     x0 = [0.5, 0.5, 0.5]; 
     lb = [0,0,0]; 
      ub = [1, 1, 1]; 
      visc = 0.001;
      dens = 1000;
      w = 0.5;

   opts = optimset('Display','iter','Algorithm','sqp',...
                  'MaxFunEval',inf,'MaxIter',Inf);

[x,fval,flag] = fmincon(@phenol,x0,[],[],[],[],lb,ub,@constraints,opts);

Answer 2

[Alan Weiss]

I wonder if you have an error in this line:

 J = w*f1 + (w-1)*f2; %weightage

Did you mean this instead?

 J = w*f1 + (1-w)*f2; %weightage

Also, the reason that you had trouble with the initial condition

x0 = [-0.5, 0, 0]

is this line from Pdrop:

 ...dens*((x(1))^2)/(((x(2))^3)*(x(3)))));

If x(2) and x(3) are zero, then you get a division by zero here.

Alan Weiss

MATLAB mathematical toolbox documentation