Error using fmincon, the function cannot continue

Hi, everyone.

I'm trying to optimize a function with 8 optimization variables. I set up initial point and costraints but appears the following error:

"Error using optim.problemdef.OptimizationProblem/solve

Objective function is undefined at initial point. Fmincon cannot continue.

Error in curvatura_V7 (line 75)

[sol,fval,exitflag,output] = solve(prob,x0);"

the code is as follow:

aymax = 5;
V = 70;
deltay = 60;
deltax = 3;
intorno = sqrt(deltax^2+deltay^2);
dmin = 3; 

type objfun

rof1 = optimvar('rof1');
rof2 = optimvar('rof2');
rof3 = optimvar('rof3');
rof4 = optimvar('rof4');
sf1 = optimvar('sf1');
sf2 = optimvar('sf2');
sf3 = optimvar('sf3');
sf4 = optimvar('sf4');

obj = fcn2optimexpr(@objfun,rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4);
prob = optimproblem('Objective',obj);

constr1 = (-aymax / (V/3.6)^2) <= rof1;
prob.Constraints.curvatura1 = constr1;
constr2 = (aymax / (V/3.6)^2) >= rof1;
prob.Constraints.curvatura1 = constr2;
constr3 = (-aymax / (V/3.6)^2) <= rof2;
prob.Constraints.curvatura2 = constr3;
constr4 = (aymax / (V/3.6)^2) >= rof2;
prob.Constraints.curvatura2 = constr4;
constr5 = (-aymax / (V/3.6)^2) <= rof3;
prob.Constraints.curvatura3 = constr5;
constr6 = (aymax / (V/3.6)^2) >= rof3;
prob.Constraints.curvatura3 = constr6;
constr7 = (-aymax / (V/3.6)^2) <= rof4;
prob.Constraints.curvatura4 = constr7;
constr8 = (aymax / (V/3.6)^2) >= rof4;
prob.Constraints.curvatura4 = constr8;

constr9 = sf1 >= 0 ;
prob.Constraints.ascissa1 = constr9;
constr10 = sf1 <= sf2;
prob.Constraints.ascissa1 = constr10;
constr11 = sf2 >= sf1;
prob.Constraints.ascissa2 = constr11;
constr12 = sf2 <= sf3;
prob.Constraints.ascissa2 = constr12;
constr13 = sf3 >= sf1+sf2;
prob.Constraints.ascissa3 = constr13;
constr14 = sf3 <= sf4;
prob.Constraints.ascissa3 = constr14;
constr15 = sf4 >= sf1+sf2+sf3;
prob.Constraints.ascissa4 = constr15;
constr16 = sf4 <= 2*intorno;
prob.Constraints.ascissa4 = constr16;

show(prob)

x0.rof1 = 0;
x0.rof2 = 0;
x0.rof3 = 0;
x0.rof4 = 0;
x0.sf1 = 0;
x0.sf2 = 0;
x0.sf3 = 0;
x0.sf4 = 0;

[sol,fval,exitflag,output] = solve(prob,x0);

Answer 1

[Walter Roberson]

aymax = 5;
V = 70;
deltay = 60;
deltax = 3;
intorno = sqrt(deltax^2+deltay^2);
dmin = 3; 

rof1 = optimvar('rof1');
rof2 = optimvar('rof2');
rof3 = optimvar('rof3');
rof4 = optimvar('rof4');
sf1 = optimvar('sf1');
sf2 = optimvar('sf2');
sf3 = optimvar('sf3');
sf4 = optimvar('sf4');

obj = fcn2optimexpr(@objfun,rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4);

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

prob = optimproblem('Objective',obj);

constr1 = (-aymax / (V/3.6)^2) <= rof1;
prob.Constraints.curvatura1 = constr1;
constr2 = (aymax / (V/3.6)^2) >= rof1;
prob.Constraints.curvatura1 = constr2;
constr3 = (-aymax / (V/3.6)^2) <= rof2;
prob.Constraints.curvatura2 = constr3;
constr4 = (aymax / (V/3.6)^2) >= rof2;
prob.Constraints.curvatura2 = constr4;
constr5 = (-aymax / (V/3.6)^2) <= rof3;
prob.Constraints.curvatura3 = constr5;
constr6 = (aymax / (V/3.6)^2) >= rof3;
prob.Constraints.curvatura3 = constr6;
constr7 = (-aymax / (V/3.6)^2) <= rof4;
prob.Constraints.curvatura4 = constr7;
constr8 = (aymax / (V/3.6)^2) >= rof4;
prob.Constraints.curvatura4 = constr8;

constr9 = sf1 >= 0 ;
prob.Constraints.ascissa1 = constr9;
constr10 = sf1 <= sf2;
prob.Constraints.ascissa1 = constr10;
constr11 = sf2 >= sf1;
prob.Constraints.ascissa2 = constr11;
constr12 = sf2 <= sf3;
prob.Constraints.ascissa2 = constr12;
constr13 = sf3 >= sf1+sf2;
prob.Constraints.ascissa3 = constr13;
constr14 = sf3 <= sf4;
prob.Constraints.ascissa3 = constr14;
constr15 = sf4 >= sf1+sf2+sf3;
prob.Constraints.ascissa4 = constr15;
constr16 = sf4 <= 2*intorno;
prob.Constraints.ascissa4 = constr16;

show(prob)

  OptimizationProblem : 

	Solve for:
       rof1, rof2, rof3, rof4, sf1, sf2, sf3, sf4

	minimize :
       objfun(rof1, rof2, rof3, rof4, sf1, sf2, sf3, sf4)


	subject to curvatura1:
       0.013224 >= rof1

	subject to curvatura2:
       0.013224 >= rof2

	subject to curvatura3:
       0.013224 >= rof3

	subject to curvatura4:
       0.013224 >= rof4

	subject to ascissa1:
       sf1 - sf2 <= 0

	subject to ascissa2:
       sf2 - sf3 <= 0

	subject to ascissa3:
       sf3 - sf4 <= 0

	subject to ascissa4:
       sf4 <= 120.1499
     

x0.rof1 = 0;
x0.rof2 = 0;
x0.rof3 = 0;
x0.rof4 = 0;
x0.sf1 = 0;
x0.sf2 = 0;
x0.sf3 = 0;
x0.sf4 = 0;

%[sol,fval,exitflag,output] = solve(prob,x0);

result = objfun(x0.rof1, x0.rof2, x0.rof3, x0.rof4, x0.sf1, x0.sf2, x0.sf3, x0.sf4)

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

Warning: Inf or NaN value encountered.

result = NaN

function erroredistanza = objfun(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4)
    deltay = 60;
    deltax = 3;
    tau0 = 90;
    tau0 = deg2rad(tau0);
    dmin = 3;
    roi1=0;
    si1=0;
    
    PSI1 = @(roi1,rof1,si1,sf1,s)  tau0 + ((rof1*si1-roi1*sf1)/(si1-sf1))*s + ((roi1-rof1)/(si1-sf1))*(s.^2)/2;
    PSI2 = @(rof1,rof2,sf1,sf2,s) ((rof2*sf1-rof1*sf2)/(sf1-sf2))*s + ((rof1-rof2)/(sf1-sf2))*(s.^2)/2;
    PSI3 = @(rof2,rof3,sf2,sf3,s) ((rof3*sf2-rof2*sf3)/(sf2-sf3))*s + ((rof2-rof3)/(sf2-sf3))*(s.^2)/2;
    PSI4 = @(rof3,rof4,sf3,sf4,s) ((rof4*sf3-rof3*sf4)/(sf3-sf4))*s + ((rof3-rof4)/(sf3-sf4))*(s.^2)/2;
    
    funx1 = @(roi1,rof1,si1,sf1,s) integral(@(s) cos(PSI1(roi1,rof1,si1,sf1,s)),0,sf1);
    funx2 = @(rof1,rof2,sf1,sf2,s) integral(@(s) cos(PSI2(rof1,rof2,sf1,sf2,s)),sf1,sf2);
    funx3 = @(rof2,rof3,sf2,sf3,s) integral(@(s) cos(PSI3(rof2,rof3,sf2,sf3,s)),sf2,sf3);
    funx4 = @(rof3,rof4,sf3,sf4,s) integral(@(s) cos(PSI4(rof3,rof4,sf3,sf4,s)),sf3,sf4);
    
    funy1 = @(roi1,rof1,si1,sf1,s) integral(@(s) sin(PSI1(roi1,rof1,si1,sf1,s)),0,sf1);
    funy2 = @(rof1,rof2,sf1,sf2,s) integral(@(s) sin(PSI2(rof1,rof2,sf1,sf2,s)),sf1,sf2);
    funy3 = @(rof2,rof3,sf2,sf3,s) integral(@(s) sin(PSI3(rof2,rof3,sf2,sf3,s)),sf2,sf3);
    funy4 = @(rof3,rof4,sf3,sf4,s) integral(@(s) sin(PSI4(rof3,rof4,sf3,sf4,s)),sf3,sf4);
    
    funx = @(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) funx1(roi1,rof1,si1,sf1) + funx2(rof1,rof2,sf1,sf2) + funx3(rof2,rof3,sf2,sf3) + funx4(rof3,rof4,sf3,sf4);
    funy = @(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) funy1(roi1,rof1,si1,sf1) + funy2(rof1,rof2,sf1,sf2) + funy3(rof2,rof3,sf2,sf3) + funy4(rof3,rof4,sf3,sf4);
    
    fun2x =@(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) funx(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) - deltax;
    fun2y =@(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) funy(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) - deltay;
    fun3x =@(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) fun2x(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4).^2;
    fun3y =@(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) fun2y(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4).^2;
    d =@(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) sqrt(fun3x(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) + fun3y(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4));
    erroredistanza = d(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) - dmin;
end

... Look at the code. You divide by (si1-sf1) and (sf1-sf2) and (sf2-sf3) and (sf3-sf4) . si1 is 0, and sf1 to sf4 are passed in. Your initial values for all of the variables is 0, so you are dividing by (0-0) in multiple places.

If you divide by something in your objective, it is your responsibility to ensure that your constraints exclude the possibility that you are going to be dividing by 0, and that your initial point is not going to trigger a division by 0.

[Andrea Mazzetto]

Ok fixed, thank you!

Answer 2

[Tarunbir Gambhir]

The error message says the objective function is not defined at initial point. This means objfun(xo) returns a 'NaN' or errors out (where "objfun" is objective function and "xo" is the initial point). It could be that the initial point does not satisfy all the provided constraints for the objective function. Since you did not provide the 'objfun', I would not be able to help any further.

Anyway, to avoid this error you can simply try choosing a better initial point.

[Andrea Mazzetto]

this is my objfun

function erroredistanza = objfun(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4)
deltay = 60;
deltax = 3;
tau0 = 90;
tau0 = deg2rad(tau0);
dmin = 3;
roi1=0;
si1=0;

PSI1 = @(roi1,rof1,si1,sf1,s)  tau0 + ((rof1*si1-roi1*sf1)/(si1-sf1))*s + ((roi1-rof1)/(si1-sf1))*(s.^2)/2;
PSI2 = @(rof1,rof2,sf1,sf2,s) ((rof2*sf1-rof1*sf2)/(sf1-sf2))*s + ((rof1-rof2)/(sf1-sf2))*(s.^2)/2;
PSI3 = @(rof2,rof3,sf2,sf3,s) ((rof3*sf2-rof2*sf3)/(sf2-sf3))*s + ((rof2-rof3)/(sf2-sf3))*(s.^2)/2;
PSI4 = @(rof3,rof4,sf3,sf4,s) ((rof4*sf3-rof3*sf4)/(sf3-sf4))*s + ((rof3-rof4)/(sf3-sf4))*(s.^2)/2;

funx1 = @(roi1,rof1,si1,sf1,s) integral(@(s) cos(PSI1(roi1,rof1,si1,sf1,s)),0,sf1);
funx2 = @(rof1,rof2,sf1,sf2,s) integral(@(s) cos(PSI2(rof1,rof2,sf1,sf2,s)),sf1,sf2);
funx3 = @(rof2,rof3,sf2,sf3,s) integral(@(s) cos(PSI3(rof2,rof3,sf2,sf3,s)),sf2,sf3);
funx4 = @(rof3,rof4,sf3,sf4,s) integral(@(s) cos(PSI4(rof3,rof4,sf3,sf4,s)),sf3,sf4);

funy1 = @(roi1,rof1,si1,sf1,s) integral(@(s) sin(PSI1(roi1,rof1,si1,sf1,s)),0,sf1);
funy2 = @(rof1,rof2,sf1,sf2,s) integral(@(s) sin(PSI2(rof1,rof2,sf1,sf2,s)),sf1,sf2);
funy3 = @(rof2,rof3,sf2,sf3,s) integral(@(s) sin(PSI3(rof2,rof3,sf2,sf3,s)),sf2,sf3);
funy4 = @(rof3,rof4,sf3,sf4,s) integral(@(s) sin(PSI4(rof3,rof4,sf3,sf4,s)),sf3,sf4);

funx = @(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) funx1(roi1,rof1,si1,sf1) + funx2(rof1,rof2,sf1,sf2) + funx3(rof2,rof3,sf2,sf3) + funx4(rof3,rof4,sf3,sf4);
funy = @(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) funy1(roi1,rof1,si1,sf1) + funy2(rof1,rof2,sf1,sf2) + funy3(rof2,rof3,sf2,sf3) + funy4(rof3,rof4,sf3,sf4);

fun2x =@(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) funx(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) - deltax;
fun2y =@(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) funy(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) - deltay;
fun3x =@(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) fun2x(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4).^2;
fun3y =@(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) fun2y(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4).^2;
d =@(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) sqrt(fun3x(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) + fun3y(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4));
erroredistanza = d(rof1,rof2,rof3,rof4,sf1,sf2,sf3,sf4) - dmin;