fmincon optimization not running properly, stops beacuase of StepTolerance

2 次查看（过去 30 天）
Levente Varga
Levente Varga 2023-4-18
评论： Matt J 2023-4-21
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.

请先登录，再回答此问题。

回答（3 个）

Walter Roberson
Walter Roberson 2023-4-18
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.
  1 个评论
Levente Varga
Levente Varga 2023-4-20
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.

请先登录，再进行评论。

Zaikun Zhang
Zaikun Zhang 2023-4-18
编辑：Zaikun Zhang 2023-4-19
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.
  1 个评论
Walter Roberson
Walter Roberson 2023-4-18
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

请先登录，再进行评论。

Matt J
Matt J 2023-4-18
编辑：Matt J 2023-4-18
Ran in:
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
  2 个评论
Levente Varga
Levente Varga 2023-4-20
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
Matt J 2023-4-21
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.,

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

产品


版本

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by