Fmincon error and solver stopped prematurely?

Hi, I'm using Fmincon for optimization but when I try to run the script with the initial conditions that I give it, it gives me the following error after performing a single iter:

Converged to an infeasible point.

fmincon stopped because the size of the current search direction is less than
twice the default value of the step size tolerance but constraints are not
satisfied to within the default value of the constraint tolerance.

<stopping criteria details>

when trying to relax the restrictions the error is next after 1 iter

Local minimum found that satisfies the constraints.

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

<stopping criteria details>

Any idea why this problem???.

Answer 1

[Walter Roberson]

"Local minimum found that satisfies the constraints." on the first iteration means that all of the directions permitted by the constraints increase the function value relative to the starting point. The starting point is at the bottom of a local minimum.

One of the ways this could happen is if the constraints allow the starting point but not other points. Or you might just have gotten lucky with the starting point.

Remember, fmincon is strictly a local minimizer.

[alan ghi]

thanks but this happens to me with all the starting points that I put. and when I remove constraints of my problem.

[Walter Roberson]

Please post your code

[Walter Roberson]

Your code has

    objfunc = @(x)tiempo( x0,P,option );

which ignores x and always passes in x0. The output of that is going to be constant, so it is not possible to optimize it.

[alan ghi]

I did not understand what you're talking about. you say he's taking x0 always to solve.

[Walter Roberson]

Change

      objfunc = @(x)tiempo( x0,P,option );

to

      objfunc = @(x)tiempo( x,P,option );

Answer 2

[alan ghi]

I'm sorry my code is very long and it's very messy.

thanks for your help.

optimainaleron.m

clear all
clc
load 'resultscirculo.mat'
option = 'cerrado';
% P=P_opt;
% R=R_opt;
% clear P_opt R_opt
axmax = inf*ones(max(size(P)),1);
axmin = -inf*ones(max(size(P)),1);
aymax = inf*ones(max(size(P)),1);
aymin = -inf*ones(max(size(P)),1);
% ffrz = 0.*ones(max(size(P)),1);
% fflz = 0.*ones(max(size(P)),1);
% frrz = 0.*ones(max(size(P)),1);
% % frlz = 0.*ones(max(size(P)),1);
% Frmax=[ffrz;fflz;frrz;frlz];

vx0min = -100;
% ffrz = -inf.*ones(max(size(P)),1);
% fflz = -inf.*ones(max(size(P)),1);
% frrz = -inf.*ones(max(size(P)),1);
% frlz = -inf.*ones(max(size(P)),1);
% Frmin=[ffrz;fflz;frrz;frlz];

vx0max = 100;
vy0min = -inf;
vy0max = inf;
 v0min = -inf;
 v0max = inf;
theta0min=-18*ones(max(size(P)),1);
theta0max=40*ones(max(size(P)),1);
 lb = [axmin;aymin;theta0min;vx0min;vy0min];
 ub = [axmax;aymax;theta0max;vx0max;vy0max];
%   lb = [theta0min;vx0min;vy0min];
%  ub = [theta0max;vx0max;vy0max];
theta0 =linspace(-18,45,50);
theta0=theta0';
% theta0 = -50.*ones(max(size(P)/2),1);
% theta0 = [theta0;-50.*ones(max(size(P)/2),1)];
% ax0 = 1.*ones(max(size(P)),1);
% ay0 = 1.*ones(max(size(P)),1);
ax0 =linspace(0,0,50);
ay0 =linspace(0,0,50);

ffrz0 = -100.*ones(max(size(P)),1);
fflz0 = -100.*ones(max(size(P)),1);
frrz0 = -200.*ones(max(size(P)),1);
frlz0 = -200.*ones(max(size(P)),1);
Fr0=[ffrz0;fflz0;frrz0;frlz0];
v1=50;
vx1=10;
vy1=10;
x0=[ax0';ay0';theta0;vx1;vy1];

noncon=@(x)moddim(x0,R,P,option);
% x = fzero(noncon,x0)
objfunc = @(x)tiempo( x0,P,option );
options.Display = 'iter';
options.PlotFcns = {@optimplotfval,@optimplotfirstorderopt};
options.MaxFunEvals =10000000;
options.TolFun = 1e-5;
options.Algorithm = 'active-set';
options.UseParallel = 0;
for i = 1:1
%     [x0,fval0,exitflag0,output0,lambda0,grad0,hessian0]= fmincon(objfunc, x0, [],[],[],[],lb,ub,noncon,options);
    [x0,fval0,exitflag0,output0,lambda0,grad0,hessian0]= fmincon(objfunc, x0, [],[],[],[],[],[],noncon,options);
%      [x,fval,exitflag,output] =  patternsearch(objfunc, x0, [],[],[],[],lb,ub,noncon,options);

%     x0 = fmincon(objfunc, x0, [],[],[],[],[],[],noncon,options);

      x_opt = x0;
  end

tiempo.m

function [ velp ] = tiempo( x,P,option )

% siz = length(x)-3;
% % ay = x(1:(siz/7));
% % ax = x((siz/7)+1:(((2*siz)/7)));
% theta=x(((1)):(end-2));
% % Fr=x(((3*siz)/7)+1:end-3);
% % Fr=reshape(Fr,[],4);
% vx1=x(end-1);
% vy1=x(end);
siz = length(x)-2;
ax = x(1:(siz/3));
ay = x((siz/3)+1:(((2*siz)/3)));
theta=x(((2*siz/3)+1):(end-2));
% Fr=x(((3*siz)/7)+1:end-3);
% Fr=reshape(Fr,[],4);
vx1=x(end-1);
vy1=x(end);
np=max(size(P));
[ D, s,sx,sy ] = distancia( P,option );
[ v,vx,vy ] = velcalculator( ax,ay,vx1,vy1, s,sx,sy)

v=(vx'.^2+vy'.^2).^(1/2);

     p00 =      0.1248  ;%(-1.512, 1.762)
         p10 =   -0.008334  ;%(-0.2887, 0.2721)
         p01 =     0.03139  ;%(-0.084, 0.1468)
         p20 =   -0.002709  ;%(-0.01576, 0.01034)
         p11 =   -0.001824  ;%(-0.01569, 0.01204)
         p02 =   -0.002159  ;%(-0.008204, 0.003885)
         p30 =  -3.832e-06  ;%(-0.0002225, 0.0002148)
         p21 =  -0.0002478  ;%(-0.0008435, 0.0003478)
         p12 =  -3.538e-05  ;%(-0.0003742, 0.0003035)
         p03 =  -0.0001129  ;%(-0.0004093, 0.0001834)
         p40 =   5.228e-08  ;%(-1.143e-06, 1.247e-06)
         p31 =   2.519e-07  ;%(-9.385e-06, 9.889e-06)
         p22 =  -1.056e-05  ;%(-1.759e-05, -3.52e-06)
         p13 =   6.798e-06  ;%(-5.695e-07, 1.417e-05)
         p04 =   8.076e-06  ;%(-3.413e-06, 1.956e-05)
         p41 =  -2.174e-10  ;%(-5.245e-08, 5.202e-08)
         p32 =  -1.048e-08  ;%(-5.829e-08, 3.734e-08)
         p23 =  -4.085e-08  ;%(-1.007e-07, 1.898e-08)
         p14 =  -1.177e-07  ;%(-2.057e-07, -2.973e-08)
         p05 =  -1.094e-07  ;
         v=v';
   Fd =( p00 + p10.*theta + p01.*v + p20.*theta.^2 + p11.*theta.*v + p02.*v.^2 + p30.*theta.^3 + p21.*theta.^2.*v + p12.*theta.*v.^2 + p03.*v.^3 + p40.*theta.^4 + p31.*theta.^3.*v + p22.*theta.^2.*v.^2 + p13.*theta.*v.^3 + p04.*v.^4 + p41.*theta.^4.*v + p32.*theta.^3.*v.^2 + p23.*theta.^2.*v.^3 + p14.*theta.*v.^4 + p05.*v.^5).*9.81;

a=((ax+(( p00 + p10.*theta + p01.*v + p20.*theta.^2 + p11.*theta.*v + p02.*v.^2 + p30.*theta.^3 + p21.*theta.^2.*v + p12.*theta.*v.^2 + p03.*v.^3 + p40.*theta.^4 + p31.*theta.^3.*v + p22.*theta.^2.*v.^2 + p13.*theta.*v.^3 + p04.*v.^4 + p41.*theta.^4.*v + p32.*theta.^3.*v.^2 + p23.*theta.^2.*v.^3 + p14.*theta.*v.^4 + p05.*v.^5).*9.81)/1000).^2+ay.^2).^(1/2);
time(1)=0;

   vel = -((2.*(a(1:end-1)+a(2:end))).*(((s(end-1)+v(1:end-1).^2)))).^(0.5);
   vellast = -((2.*(a(end)+a(1))).*(((s(end)+v(end).^2)))).^(0.5);
  vel=[vel;vellast]
   velp=sum(vel)/50
  for i=1:49
  time(i+1)= (2.*s(i)./(a(i)+a(i+1)))+time(i).^2.^(0.5);
  end
  timelast = (((s(1)./(a(end)+a(1)))+time(end-1)).^(0.5));
  time =  time(end)+timelast ;
  % timev = sum(time);
  % 

% timev =  0 ;

end

moddim.m

function [c,ceq]=moddim(x,R,P,option);
siz = length(x)-2;
ax = x(1:(siz/3));
ay = x((siz/3)+1:(((2*siz)/3)));
theta=x(((2*siz/3)+1):(end-2));
% Fr=x(((3*siz)/7)+1:end-3);
% Fr=reshape(Fr,[],4);
vx1=x(end-1);
vy1=x(end);
np=max(size(P));
[ D, s,sx,sy ] = distancia( P,option );
vx=zeros(np,1);
vy=zeros(np,1);
% for i=1:np
%   [A]=calcacce(R(i),theta(i),vx1,vy1);
%   ax(i)=A(1);
%   ay(i)=A(2);
%   ar(i)=A(3);
%   Fr(:,i)=A(4:end);
[ v,vx,vy ] = velcalculator( ax,ay,vx1,vy1, s,sx,sy)
%     vx1=vx(i);
%     vy1=vy(i);
% end
% ax=ax';
% ay=ay';
% ar=ar';
% 
% vx=vx';
% vy=vy';
 a=(ax.^2+ay.^2).^(1/2);
 ar=a./R;
m=1000;d=0;b=1.7;tf=0.95;tr=0.95;h=0.48;Ixz=160;Izz=2800;Krr=(500*17.858);Krl=(500*17.858);Kfr=(700*17.858);Kfl=(700*17.858);da=1.2;l=2.6;g=9.81;
v=(vx.^2+vy.^2).^(1/2);

          p00 =      -22.07 ; %
         p10 =       -9.15  ;%(-17.42, -0.8807)
         p01 =       2.109  ;%(-17.98, 22.2)
         p20 =      0.4215 ; %(-0.01162, 0.8547)
         p11 =      0.3237 ; %(-0.6698, 1.317)
         p02 =     -0.6287 ; %(-1.564, 0.3064)
         p30 =     0.03312 ; %(0.01188, 0.05435)
         p21 =   -0.001976 ; %(-0.02626, 0.02231)
         p12 =   -0.008085 ; %(-0.05077, 0.0346)
         p03 =    6.67e-05 ; %(-0.0156, 0.01573)
         p40 =    -0.00181 ; %(-0.002633, -0.0009865)
         p31 =   -0.001228 ; %(-0.001756, -0.0007002)
         p22 =   0.0009839 ; %(0.0004797, 0.001488)
         p13 =  -1.491e-05 ; %(-0.0007055, 0.0006757)
         p04 =  -3.561e-06 ; %(-8.921e-05, 8.208e-05)
         p50 =   2.152e-05 ; %(1.04e-05, 3.264e-05)
         p41 =   2.599e-05 ; %(1.969e-05, 3.23e-05)
         p32 =  -1.592e-05 ; %(-2.021e-05, -1.163e-05)
         p23 =   4.868e-07 ; %(-2.94e-06, 3.913e-06)
         p14 =   6.517e-08 ; %(-3.678e-06, 3.808e-06)
  Fl=  p00 + p10.*theta + p01.*v + p20.*theta.^2 + p11.*theta.*v + p02.*v.^2 + p30.*theta.^3 + p21.*theta.^2.*v + p12.*theta.*v.^2 + p03.*v.^3 + p40.*theta.^4 + p31.*theta.^3.*v + p22.*theta.^2.*v.^2 + p13.*theta.*v.^3 + p04.*v.^4 + p50.*theta.^5 + p41.*theta.^4.*v + p32.*theta.^3.*v.^2 + p23.*theta.^2.*v.^3 + p14.*theta.*v.^4;
  deltaang=l./R;
  r=(vx./l).*(deltaang);
  betaang=atan(vy./abs(vx));
  alfaang1= atan(v.*(sin(betaang)))./(v.*cos((betaang))+tf)-deltaang;
  alfaang2= atan(v.*(sin(betaang)))./(v.*cos((betaang))-tf)-deltaang;
  alfaang3= atan(v.*(sin(betaang)))./(v.*cos((betaang))+tr);
  alfaang4= atan(v.*(sin(betaang)))./(v.*cos((betaang))-tr);
  alfaang = [alfaang1;alfaang2;alfaang3;alfaang4];
  g=9.81;
  for i=1:siz/3
   Fr(:,i)=calcforce(vx(i),vy(i),r(i),ax(i),ay(i),ar(i),Fl(i));
  end
   Fr1=Fr';

k=((0.45./(0.45-abs(Fr)/2.1582e+05))-1);

muflx=(ax./abs(ax)).*(1.48.*sin(1.37*atan(18.22.*k(1,:)'+0.46.*(18.22.*k(1,:)'-atan(18.22.*k(1,:)'))))).*(cos(1.1231.*atan(alfaang1.*(13.476./(1+11.354^2.*k(1,:)'.^2)))));
mufrx=(ax./abs(ax)).*(1.48.*sin(1.37*atan(18.22.*k(2,:)'+0.46.*(18.22.*k(2,:)'-atan(18.22.*k(2,:)'))))).*(cos(1.1231.*atan(alfaang2.*(13.476./(1+11.354^2.*k(2,:)'.^2)))));
murlx=(ax./abs(ax)).*(1.48.*sin(1.37*atan(18.22.*k(3,:)'+0.46.*(18.22.*k(3,:)'-atan(18.22.*k(3,:)'))))).*(cos(1.1231.*atan(alfaang3.*(13.476./(1+11.354^2.*k(3,:)'.^2)))));
murrx=(ax./abs(ax)).*(1.48.*sin(1.37*atan(18.22.*k(4,:)'+0.46.*(18.22.*k(4,:)'-atan(18.22.*k(4,:)'))))).*(cos(1.1231.*atan(alfaang4.*(13.476./(1+11.354^2.*k(4,:)'.^2)))));

mufly=(1.22.*sin(1.25*atan(17.8.*alfaang1-0.02.*(17.8.*alfaang1-atan(17.8.*alfaang1))))).*(cos(1.0533.*atan(k(1,:)'.*(7.7856./(1+8.1697^2.*alfaang1.^2)))));
mufry=(1.22.*sin(1.25*atan(17.8.*alfaang2-0.02.*(17.8.*alfaang2-atan(17.8.*alfaang2))))).*(cos(1.0533.*atan(k(2,:)'.*(7.7856./(1+8.1697^2.*alfaang2.^2)))));
murly=(1.22.*sin(1.25*atan(17.8.*alfaang3-0.02.*(17.8.*alfaang3-atan(17.8.*alfaang3))))).*(cos(1.0533.*atan(k(3,:)'.*(7.7856./(1+8.1697^2.*alfaang3.^2)))));
murry=(1.22.*sin(1.25*atan(17.8.*alfaang4-0.02.*(17.8.*alfaang4-atan(17.8.*alfaang4))))).*(cos(1.0533.*atan(k(4,:)'.*(7.7856./(1+8.1697^2.*alfaang4.^2)))));

Fflz=Fr1(:,1);
Ffrz=Fr1(:,2);
Frlz=Fr1(:,3);
Frrz=Fr1(:,4);
Fr=[Fflz;Ffrz;Frlz;Frrz];
Ax=(ax.*m-ar.*d.*m+Fflz.*muflx+Ffrz.*mufrx+Frlz.*murlx+Fl.*murrx+Frrz.*murrx-(b.*r+vy).*m.*r);
Ay=(ay.*m+ar.*b.*m+Fflz.*mufly+Ffrz.*mufry+Frlz.*murly+Frrz.*murrx+Fl.*murry+m.*r.*(-d.*r+vx));
My=(ay.*b.*m-ax.*d.*m+ar.*(Izz+b.*h.*m)+Ffrz.*(l.*mufly-muflx.*tf)+Fflz.*(l.*mufly+muflx.*tf)+Frlz.*murlx.*tr-Frrz.*murrx.*tr+m.*r.*(b.*vx+d.*vy));
Az=(-Fflz-Ffrz-Fl-Frlz-Frrz-g.*m);
Mp=(-d.*g.*m+ay.*h.*m+ar.*(Ixz+b.*h.*m)+Fflz.*tf-Ffrz.*tf+Frlz.*tr-Frrz.*tr+h.*m.*r.*(-d.*r+vx));
Mr=(Fflz.*l+Ffrz.*l+b.*g.*m-ax.*h.*m+ar.*d.*h.*m+(Ixz.*+b.*h.*m).*r.^2+h.*m.*r.*vy);
Hp=(-Frlz.*Krl.*tf+Frrz.*Krr.*tf+Fflz.*Kfl.*tr-Ffrz.*Kfr.*tr);
%   c=[Ax;Ay;My;Fr];
    c=[];
% %  ceq=[Ax;Ay;My;Mp;Mr;Fr;Az;Hp];
   ceq=[];
%   ceq=[Az;Hp];
%  ceq=[Hp];

end

[Walter Roberson]

To test we will need your resultscirculo.mat and distancia.m and velcalculator.m and calcforce.m (I might have overlooked something)

[alan ghi]

thank you very much for the help I found the error in the code

Answer 3

[alan ghi]

hello walter I send you the scripts thank you very much for your help