System of Nonlinear Equations exceeds function evaluation limit

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.

Answer 1

[Star Strider]

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.

.

[Brian]

Thank you! is there a way to bound the outputs? e.g. must > 0

[Star Strider]

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.)

.

Answer 2

[Torsten]

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