fmincon optimization not running properly, stops beacuase of StepTolerance

I'm designing a strain wave gear for my bsc diploma, but I cannot run this optimization, and I ran out of people I can ask.

This is my code so far:

%Input paramaters:
beta = 62.31*pi/180;
r_0 = 24.519; %eq 3
e = 0.44;
n1 = 1;
n2 = 2;
n3 = 3;

%Input equalites and inequalites from the research:
ineq1_1 = (4*n1^2-1)*cos(2*n1*beta); %eq 8
ineq1_2 = (4*n2^2-1)*cos(2*n2*beta); %eq 8
ineq1_3 = (4*n3^2-1)*cos(2*n3*beta); %eq 8
ineq2_1 = 4*n1^2*(4*n1^2-1)*cos(n1*pi); %eq 9
ineq2_2 = 4*n2^2*(4*n2^2-1)*cos(n2*pi); %eq 9
ineq2_3 = 4*n3^2*(4*n3^2-1)*cos(n3*pi); %eq 9
lambda1 = (4*n1^2-1)*cos(2*n1*beta);
lambda2 = (4*n2^2-1)*cos(2*n2*beta);
lambda3 = (4*n3^2-1)*cos(2*n3*beta);
coef1 = (1/r_0^2)*lambda1; %eq 6
coef2 = (1/r_0^2)*lambda2; %eq 6
coef3 = (1/r_0^2)*lambda3; %eq 6



%Constarints:
x02 = [10,5,1];
lb = [];
ub = [];
nonlcon = [];
Aeq2 = [1,1,1];
beq2 = e;
A2 = [ineq1_1,ineq1_2,ineq1_3;ineq2_1,ineq2_2,ineq2_3];
b2 = [r_0,0];
fun2 = @(a)coef1*a(1)+coef2*a(2)+coef3*a(3);

options = optimoptions('fmincon','EnableFeasibilityMode',true,'Algorithm','sqp','StepTolerance',1e-2);

x2 = fmincon(fun2,x02,A2,b2,Aeq2,beq2,lb,ub,nonlcon,options);

It stops with this message:

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than

the value of the step size tolerance and constraints are

satisfied to within the value of the constraint tolerance.

<stopping criteria details>

Changing the StepTolerance didn't fix it. It should give back 0.450; -0,005; -0,005.

Before this code, I wrote one with only two unknowns, where I calculated the ineqs, lambdas and coefs by hand, and I even forgot to convert degrees to radian and my calculator was in degress instead of radian, but this gave me the good results (0,0463; -0,023):

x02 = [100,100];
Aeq2 = [1,1];
beq2 = 0.44;
A2 = [-1.704,-5.317;11.98,238.55];
b2 = [24.519,0];
fun2 = @(a)-2.835*10^-3*a(1)-8.844*10^-3*a(2);
x2 = fmincon(fun2,x02,A2,b2,Aeq2,beq2);

But this one expanded to have 3 unknowns didn't work either.

I attached the research where the whole process is described and I worked from this.

Answer 1

[Walter Roberson]

fmincon is a general minimizer. It does not stop when the function value is zero: it is fine with going pretty negative.

As such, the only natural way for fmincon to be sure it has reached the best minimum possible would be if the function returned negative infinity.

In all other cases, fmincon has been coded with a number of different artificial ways to know when to give up. For example it knows to give up if the number of function evaluations reaches a configured maximum.

One of the artificial stopping criteria is a combination of circumstances: every time it changes any of the model parameters a little, the function value gets worse; and in trying to pinpoint the best possible spot near where it is, the size of the steps it is taking has reached a configured minimum.

In other words, as far as fmincon can tell, you are in a minimum (the gradient curls upwards in all directions) and maybe you could get a better position by adjusting the position half an iota of a smidgen, but you are darn close to the minimum near the current position.

This does not establish that you have reached a global minima. You cannot tell whether you have reached a global minima without a mathematical analysis of the function, which is not something that the technology can potentially do (you would need symbolic analysis to have a hope of that)

If you have not reached the model parameters that you expected, and those values are more than a trifle different, then adjusting the minimum step size is not going to help. You might have reached a local minima that is too deep for the algorithm to escape. Or you might have programmed the equations incorrectly.

[Levente Varga]

I'm sure I don't even understand what I'm doing. In the research you can see in table 2 what kind of results I expect. I didn't know how fmincon or the math behind this works sadly. I need this for the optimization of an ellipse shaped wavegenerator of a strain wave gear. The ellipse should be modified with these values as terms of a Fourier expansion.

Answer 2

[Zaikun Zhang]

I conducted a quick experiment using the function and data in the question. The solvers are fmincon of MATLAB and PRIMA, which is available at https://github.com/libprima/prima and https://www.mathworks.com/matlabcentral/fileexchange/127983-prima .

Here is the code I used (the defition of the problem is essnetially the same as that provided by @Levente Varga).

%Input paramaters:
beta = 62.31*pi/180;
r_0 = 24.519; %eq 3
e = 0.44;
n1 = 1;
n2 = 2;
n3 = 3;
%Input equalites and inequalites from the research:
ineq1_1 = (4*n1^2-1)*cos(2*n1*beta); %eq 8
ineq1_2 = (4*n2^2-1)*cos(2*n2*beta); %eq 8
ineq1_3 = (4*n3^2-1)*cos(2*n3*beta); %eq 8
ineq2_1 = 4*n1^2*(4*n1^2-1)*cos(n1*pi); %eq 9
ineq2_2 = 4*n2^2*(4*n2^2-1)*cos(n2*pi); %eq 9
ineq2_3 = 4*n3^2*(4*n3^2-1)*cos(n3*pi); %eq 9
lambda1 = (4*n1^2-1)*cos(2*n1*beta);
lambda2 = (4*n2^2-1)*cos(2*n2*beta);
lambda3 = (4*n3^2-1)*cos(2*n3*beta);
coef1 = (1/r_0^2)*lambda1; %eq 6
coef2 = (1/r_0^2)*lambda2; %eq 6
coef3 = (1/r_0^2)*lambda3; %eq 6
%Constarints:
x02 = [10,5,1]';
lb = [];
ub = [];
nonlcon = [];
Aeq2 = [1,1,1];
beq2 = e;
A2 = [ineq1_1,ineq1_2,ineq1_3;ineq2_1,ineq2_2,ineq2_3];
b2 = [r_0,0]';  % N.B.: Mathematically speaking, b2 should be a column.
fun2 = @(a)coef1*a(1)+coef2*a(2)+coef3*a(3);

fprintf('\nTesting fmincon:\n')
options = optimoptions('fmincon','EnableFeasibilityMode',true,'Algorithm','sqp','StepTolerance',1e-2);
x_fmincon = fmincon(fun2,x02,A2,b2,Aeq2,beq2,lb,ub,nonlcon,options);
fprintf('Solution provided by fmincon: x = [%g; %g; %g]\n', x_fmincon(1), x_fmincon(2), x_fmincon(3));
fprintf('f(x) = %.16f\n', fun2(x_fmincon));
fprintf('Constraint violation: %.16f\n\n', max([0; A2*x_fmincon-b2; abs(Aeq2*x_fmincon-beq2)]));

fprintf('Testing PRIMA:\n')
warning('off', 'prima:ReviseX0');
x_prima = prima(fun2,x02,A2,b2,Aeq2,beq2,lb,ub,nonlcon);
fprintf('Solution provided by PRIMA: x = [%g; %g; %g]\n', x_prima(1), x_prima(2), x_prima(3));
fprintf('f(x) = %.16f\n', fun2(x_prima));
fprintf('Constraint violation: %.16f\n\n', max([0; A2*x_prima-b2; abs(Aeq2*x_prima-beq2)]));

x = [0.450; -0.005; -0.005];
fprintf('Reference solution: x = [0.450; -0.005; -0.005]\n');
fprintf('f(x) = %.16f\n', fun2(x));
fprintf('Constraint violation: %.16f\n\n', max([0; A2*x-b2; abs(Aeq2*x-beq2)]));

Here is the result I obtained.

Testing fmincon:

Local minimum possible. Constraints satisfied.

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

Solution provided by fmincon: x = [3.48821; -2.53259; -0.515619]
f(x) = -0.0166358927551636
Constraint violation: 0.0000000000000001

Testing PRIMA:
Solution provided by PRIMA: x = [8.81918; -6.96796; -1.41122]
f(x) = -0.0431461784076583
Constraint violation: 0.0000000000000005

Reference solution: x = [0.450; -0.005; -0.005]
f(x) = -0.0015141812204616
Constraint violation: 0.0000000000000000

Conclusions:

1. Both fmincon and PRIMA achieve feasiblity up to the machine precision.
2. Both fmincon and PRIMA provide feasible solutions much better than the reference point [0.450; -0.005; -0.005] in terms of the objective function value.
3. PRIMA achives a much smaller function value than fmincon.

I hope this is helpful.

[Walter Roberson]

format long g

%Input paramaters:
beta = 62.31*pi/180;
r_0 = 24.519; %eq 3
e = 0.44;
n1 = 1;
n2 = 2;
n3 = 3;
%Input equalites and inequalites from the research:
ineq1_1 = (4*n1^2-1)*cos(2*n1*beta); %eq 8
ineq1_2 = (4*n2^2-1)*cos(2*n2*beta); %eq 8
ineq1_3 = (4*n3^2-1)*cos(2*n3*beta); %eq 8
ineq2_1 = 4*n1^2*(4*n1^2-1)*cos(n1*pi); %eq 9
ineq2_2 = 4*n2^2*(4*n2^2-1)*cos(n2*pi); %eq 9
ineq2_3 = 4*n3^2*(4*n3^2-1)*cos(n3*pi); %eq 9
lambda1 = (4*n1^2-1)*cos(2*n1*beta);
lambda2 = (4*n2^2-1)*cos(2*n2*beta);
lambda3 = (4*n3^2-1)*cos(2*n3*beta);
coef1 = (1/r_0^2)*lambda1; %eq 6
coef2 = (1/r_0^2)*lambda2; %eq 6
coef3 = (1/r_0^2)*lambda3; %eq 6
%Constarints:
x02 = [10,5,1]';
lb = [];
ub = [];
nonlcon = [];
Aeq2 = [1,1,1];
beq2 = e;
A2 = [ineq1_1,ineq1_2,ineq1_3;ineq2_1,ineq2_2,ineq2_3];
b2 = [r_0,0]';  % N.B.: Mathematically speaking, b2 should be a column.
fun2 = @(a)coef1*a(1)+coef2*a(2)+coef3*a(3);

options  = optimoptions('ga', 'MaxGenerations', 1e5);

[x_ga, f_ga] = ga(fun2, length(x02), A2, b2, Aeq2, beq2, lb, ub, nonlcon, options)
%double check
fun2(x_ga)

Unfortunately, most of the time this takes more than 55 second run-time limit for the Answers facility.

With coordiantes of [14144.5938655468 -11767.9316586277 -2376.22120927818] the result is -70.3386331774086

Answer 3

[Matt J]

It should give back 0.450; -0,005; -0,005.

If you have that foreknowledge, you should apply lb and ub bounds as I have done below. Also, StepTolerance=1e-2 is way too large.

%Constarints:
x02 = [10,5,1];
lb = [-1,-1,-1];
ub = [1,1,1];
nonlcon = [];
Aeq2 = [1,1,1];
beq2 = e;
A2 = [ineq1_1,ineq1_2,ineq1_3;ineq2_1,ineq2_2,ineq2_3];
b2 = [r_0,0];
fun2 = @(a)coef1*a(1)+coef2*a(2)+coef3*a(3);

options = optimoptions('fmincon','EnableFeasibilityMode',true,...
                      'Algorithm','sqp','StepTolerance',1e-10);

[x2,fval,flag] = fmincon(fun2,x02,A2,b2,Aeq2,beq2,lb,ub,nonlcon,options)

Local minimum found that satisfies the constraints.

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

x2 = 1×3

1.0000 -0.4624 -0.0976

fval = -0.0043

flag = 1

These are not the x2 values you expect, but if we evaluate the objective fun2() at the expected point, we get a worse value than what the code above is giving you. So, either your expectations are wrong, or your problem definition is.

fun2([0.450; -0.005; -0.005])

ans = -0.0015

ans =

-0.0015

[Levente Varga]

There is a good chance my definition is not right, I'm not good with math and MATLAB. Can you maybe look at the research paper I attached? The constrain and objective functions are defined there, as well as the results of the researchers. I try to reproducate their results first, and then applying the minimizer to my parameters.

[Matt J]

The paper isn't going to be an effective way to communicate the optimization problem you are trying to solve to an unspecialized audience. It is advisable for you to write the problem in self-contained mathematical notation using the rich text capabilities of the Answers editor, e.g.,