System of Nonlinear Equations exceeds function evaluation limit

9 views (last 30 days)
Attempting to solve a system of three nonlinear equaitons
g = 9.81;
Az_mech = 0.028*g;
Az_dyn = 0.5*g;
Ax_mech = 0.092*g;
Ax_dyn = 2*g;
m=4;
mu = 0.2;
f = @(x) [70*(-m*g+m*(Az_mech+Az_dyn))+55*m*(Ax_mech+Ax_dyn)+140*(x(2)*cos(x(3))+mu*x(2)*sin(x(3)));
-mu*x(1)*cos(x(3))+x(1)*sin(x(3))+m*(Ax_mech+Ax_dyn)-x(2)*sin(x(3))+mu*x(2)*cos(x(3));
x(1)*cos(x(3))+mu*x(1)*sin(x(3))+x(2)*cos(x(3))+mu*x(2)*sin(x(3))+m*(Az_mech+Az_dyn)-m*g];
options = optimoptions(@fsolve, 'MaxFunctionEvaluations', 10000, 'MaxIterations', 10000);
x1 = fsolve(f, [10 10 0])
f(x1)
x2 = fsolve(f, [5 5 0])
f(x2)
x3 = fsolve(f, [0 0 0])
f(x3)
The outputs are:
Solver stopped prematurely.
fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 300 (the default value).
x1 =
9.9967 19.7237 -3.0559
ans =
420.0052
80.9845
-48.6414
Solver stopped prematurely.
fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 300 (the default value).
x2 =
4.9978 17.7964 -2.8936
ans =
680.8757
82.7503
-41.7372
Solver stopped prematurely.
fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 300 (the default value).
x3 =
0.0024 -17.9956 0.3792
ans =
691.5529
85.4092
-36.5683
I try a few different starting points and each gives different solutions. When I try to input the solution into the original function, the result is not zero, as I would expect for a correct solution.
There is no noticable different in runtime, even if i increase the max iterations and evaluations to 1000000 from 300.

Accepted Answer

Star Strider
Star Strider on 5 Oct 2022
You need to tell fsolve that the options sturcture exists —
g = 9.81;
Az_mech = 0.028*g;
Az_dyn = 0.5*g;
Ax_mech = 0.092*g;
Ax_dyn = 2*g;
m=4;
mu = 0.2;
f = @(x) [70*(-m*g+m*(Az_mech+Az_dyn))+55*m*(Ax_mech+Ax_dyn)+140*(x(2)*cos(x(3))+mu*x(2)*sin(x(3)));
-mu*x(1)*cos(x(3))+x(1)*sin(x(3))+m*(Ax_mech+Ax_dyn)-x(2)*sin(x(3))+mu*x(2)*cos(x(3));
x(1)*cos(x(3))+mu*x(1)*sin(x(3))+x(2)*cos(x(3))+mu*x(2)*sin(x(3))+m*(Az_mech+Az_dyn)-m*g];
options = optimoptions(@fsolve, 'MaxFunctionEvaluations', 10000, 'MaxIterations', 10000);
x1 = fsolve(f, [10 10 0],options)
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
x1 = 1×3
-65.8834 36.4872 -3.8490
f(x1)
ans = 3×1
1.0e-12 * -0.9095 0.0178 0.0071
x2 = fsolve(f, [5 5 0],options)
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
x2 = 1×3
-65.8834 36.4872 -3.8490
f(x2)
ans = 3×1
1.0e-11 * 0.5002 -0.0002 0.0036
x3 = fsolve(f, [0 0 0],options)
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
x3 = 1×3
65.8834 -36.4872 -0.7074
f(x3)
ans = 3×1
1.0e-14 * 0 0.7994 0.7105
There are likely myriad (fi not an infinity) of solutions because one of the parameters (‘x(3)’) is an arguments to trigonometric funcitons.
.
  2 Comments
Star Strider
Star Strider on 5 Oct 2022
As always, my pleasure!
The fsolve function does not allow constraints. The fsolve funciton is a root-finder, however if the minima are at the zero-crossings, fmincon could be an option —
g = 9.81;
Az_mech = 0.028*g;
Az_dyn = 0.5*g;
Ax_mech = 0.092*g;
Ax_dyn = 2*g;
m=4;
mu = 0.2;
f = @(x) [70*(-m*g+m*(Az_mech+Az_dyn))+55*m*(Ax_mech+Ax_dyn)+140*(x(2)*cos(x(3))+mu*x(2)*sin(x(3)));
-mu*x(1)*cos(x(3))+x(1)*sin(x(3))+m*(Ax_mech+Ax_dyn)-x(2)*sin(x(3))+mu*x(2)*cos(x(3));
x(1)*cos(x(3))+mu*x(1)*sin(x(3))+x(2)*cos(x(3))+mu*x(2)*sin(x(3))+m*(Az_mech+Az_dyn)-m*g];
options = optimoptions(@fmincon, 'MaxFunctionEvaluations', 10000, 'MaxIterations', 10000);
format long
x1 = fmincon(@(x)norm(f(x)), [10 10 0],[],[],[],[],zeros(1,3),[],[],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.
x1 = 1×3
0.000001231879232 0.000000000218967 1.016044901310691
f(x1)
ans = 3×1
1.0e+03 * 3.218464802001204 0.082090176498341 -0.018521199567835
x2 = fmincon(@(x)norm(f(x)), [5 5 0],[],[],[],[],zeros(1,3),[],[],options)
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.
x2 = 1×3
0.371864025910980 84.129346717274629 39.738544078784273
f(x2)
ans = 3×1
0.301147983536339 -0.203252665739122 -41.609768616988219
x3 = fmincon(@(x)norm(f(x)), [0 0 0],[],[],[],[],zeros(1,3),[],[],options)
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.
x3 = 1×3
0.014702062575239 83.706997747634830 33.456762357334291
f(x3)
ans = 3×1
0.299182857958385 -0.107027744891436 -41.512214619537083
These results do not appear to be the function roots.
I believe this is the only constrained optimisation function available in the Optimization Toolbox. (Other functions exist in the Global Optimization Toolbox. They are optimisers, not root-finders, so are likely to produce the same results.)
.

Sign in to comment.

More Answers (1)

Torsten
Torsten on 5 Oct 2022
Edited: Torsten on 5 Oct 2022
Substitute
a = x(2)*cos(x(3)), b = x(2)*sin(x(3)), c= x(1)*cos(x(3)), d=x(1)*sin(x(3))
in your equations.
You get an underdetermined linear system of equations.
The general solution is
v = lambda*n + sol
with a parameter lambda, n the nullspace of the system matrix and sol a special solution.
Incorporating the condition b/a = d/c gives
lambda = vpasolve(v(2)/v(1)-v(4)/v(3)==0);
Now
x3 = atan(v(2)/v(1))
x2 = v(1)/cos(x3)
x1 = v(3)/cos(x3)
follow.
The calculation of x3 gives you degrees of freedom for the solution because x3 + n*pi will also satisfy the equation. But it turns out that - although there are infinitely many solutions - none of them satisfies that x1, x2 and x3 are simultaneously positive.
Code is as follows:
g = 9.81;
Az_mech = 0.028*g;
Az_dyn = 0.5*g;
Ax_mech = 0.092*g;
Ax_dyn = 2*g;
m=4;
mu = 0.2;
M = [140 140*mu 0 0; mu -1 -mu 1; 1 mu 1 mu];
n = null(M)
n = 4×1
-0.1387 0.6934 -0.1387 0.6934
rhs = [-55*m*(Ax_mech+Ax_dyn)-70*(-m*g+m*(Az_mech+Az_dyn));-m*(Ax_mech+Ax_dyn);-m*(Az_mech+Az_dyn)+m*g];
sol = M\rhs;
syms lambda
v = lambda*n + sol;
s = vpasolve(v(2)/v(1)-v(4)/v(3)==0);
v = double(subs(v,lambda,s));
x3 = atan(v(2)/v(1))
x3 = -0.7074
x3s = atan(v(4)/v(3))
x3s = -0.7074
x1 = v(4)/sin(x3)
x1 = 65.8834
x1s = v(3)/cos(x3)
x1s = 65.8834
x2 = v(1)/cos(x3)
x2 = -36.4872
x2s = v(2)/sin(x3)
x2s = -36.4872

Categories

Find more on Systems of Nonlinear Equations in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by