3.9 Iterative Solution of the Loop Closure Equations

In the previous sections we were able to solve the loop closure stepwise or in closed form. Usually such a solution can be found if in each loop closure equation there are two unknowns. In cases where the loop closure equations are coupled i.e. there are more than two unknowns in each loop closure equation, simultaneous solution of the loop closure equations will be required, in which case a search for a stepwise or closed form solution is not practical and iterative solutions are used.

If we look at the loop closure equations, there is a set of nonlinear algebraic equations in the form:

Fj(x) = 0	j = 1, 2, 3, …, n	(1)
	x = [x1 x2 x3 … xn]T

Where in mechanism analysis F_k(x) are the loop equations written in complex numbers or in terms of their real and imaginary parts and x is the vector formed by the position variables such as θ_1j or s_ij. When the set of equations and their derivatives are continuous, a simple numerical iterative method known as Newton-Raphson method may be used to find a solution of such an equation set. There are other more elegant numerical methods which can be used but this is beyond the topic (any textbook on numerical methods will contain a chapter on the solution of nonlinear equations).

Let us consider a single equation F(x) = 0, that is nonlinear in the scalar variable x. Let x = x* be the solution of this equation and let x^(k) be an approximation to x*. The Taylor series expansion of F(x) about x = x^(k) is:

F(x) = F(x(k)) + Fx(x(k))(x − x(k)) + (higjer order terms)

(2)

where F_x  is the derivative of the function with respect to x, F_x(x) = dF(x)/dx.

Let x = x^(k+1) be an improved approximation that is to be determined from the condition:

F(x(k+1)) ≈ F(x(k)) + Fx(x(k))(x(k+1) − x(k)) = 0

(3)

If x^(k+1) − x^(k) is small, higher order terms can be neglected. If F_x(x^(k)) is not equal to zero, equation (3) can be solved for x^(k+1) as:

x(k+1) = x(k) − F(x(k))/Fx(x(k))

(4)

Equation 4 is the basis for Newton-Raphson method (this is also known as the recursive formula). Since higher order terms in Taylor series has been neglected, x^(k+1) will not be the exact solution. If the problem is well behaving , x^(k+1) will be a better approximation, and when this equation is applied repetitively, at an mth step,  F(x^(m)) < ε, where ε is a very small quantity which may be assumed to be equal to zero for the engineering application. Hence Newton-Raphson algorithm can be summarised as:

• Make an estimate x⁽⁰⁾ for the solution of the equation F(x) = 0.
• In iterations k = 0, 1, 2, 3, …, m, evaluate F(x^(k+1)) and if F(x^(k+1)) < ε_e or if x^(k+1) − x^(k) < ε_s , then terminate and use x^(k+1) as the solution of the equation (ε_e is the equation error tolerance and ε_s is the solution error tolerance). If F_x(x^(k+1)) = 0 and F(x^(k+1)) not equal to 0 then go to step 1 for a new estimate. Otherwise go to step 3.
• Replace k+1 in place of k and go to step 2. If the step is repeated m times and the conditions F(x^(k+1)) < ε_e or x^(k+1) − x^(k) < ε_s are not satisfied, stop without a solution.

The sequence of approximate solutions generated by Newton-Raphson algorithm, in many problems will approach a solution F(x) = 0 in a reasonable number of steps. However if x⁽⁰⁾ is not adequately close to x*,  if F(x) is not a well behaving function or if ε_e, or ε_s are selected very small for the problem at hand, a solution may not be obtained. In certain other cases, the solution obtained may not be the solution x* sought.

    

The geometry of Newton-Raphson method is illustrated in the above sketches. For a normal case, the convergence is quadratic (a). In (b) shows that the method may diverge if the solution F(x) = 0 is an inflection point (i.e. the curvature changes at F(x) = 0). If there are multiple roots, the solution found will be a function of the initial estimate x⁽⁰⁾ and you can not be sure whether the root found is the actual root you are looking for (c). If the initial estimate is near a local maximum or minimum (at a minimum or maximum F_x(x) = 0), the algorithm may also fail (d)).

Let us now expand Newton-Raphson algorithm for a set of equations. Consider n nonlinear algebraic equations in n unknowns:

F1(x) = F1(x1, x2, x3, …., xn) = 0 F2(x) = F2(x1, x2, x3, …, xn) = 0 F3(x) = F3(x1, x2, x3, …, xn) = 0 …… Fn(x) = Fn(x1, x2, x3, …, xn) = 0

(5)

which can be written as :

F(x) = 0

where the solution x* = [x*₁, x*₂, x*₃, …, x*_n]^T. A first-order Taylor series expansion about an estimate x^(k) of the solution can be applied to equation (5). Neglecting higher order terms:

F(x) ≈ F(x(k)) + Fx(x(k))(x − x(k))

(6)

where F_x(x^(k)) is the Jacobian matrix evaluated at x^(k). If  x = x^(k+1) is an improved approximation, it can be determined from the condition:

F(x^(k+1)) ≈ F(x^(k)) + F_x(x^(k))(x^(k+1) − x^(k)) = 0

---

If the Jabobian matrix is nonsingular (it is of rank n), this equation can be solved for x^(k+1) from the equation :

Fx(x(k))Dx(k)= −F(x(k))

(7)

where Dx^(k)= x^(k+1) − x^(k) or x^(k+1)= x^(k) + Dx^(k).

Newton-Raphson algorithm for a system of equations will now be:

1. Make an estimate x⁽⁰⁾ for the solution of the equation F(x) = 0.
2. In iterations k= 0, 1, 2, 3, …, m,  evaluate Dx^(k), x^(k+1) and F(x^(k+1)). If F(x^(k+1)) < ε_e, or if x^(k+1) − x^(k) < ε_s, then terminate and use x^(k+1) as the solution of the equation (ε_e is the equation error tolerance and ε_s is the solution error tolerance). If F_x(x^(k+1))) = 0 and F(x^(k+1)) not equal to 0 then go to step 1 for a new estimate. Otherwise go to step 3.
3. Replace k+1 in place of k and go to step 2. If the step is repeated m times and the conditions F(x^(k+1)) < ε_e or x^(k+1) − x^(k) < ε_s are not satisfied,  stop without a solution.

Example:

In order to show the application of Newton-Raphson method in mechanism analysis consider the mechanism shown. The closed form solution is not known. The loop closure equation in complex numbers yield:

a₂e^iθ₁₂ = b₁ − ic₁ + a₅e^iθ₁₅ + s₄e^iθ₁₄

a₅e^iθ₁₅ + a₄e^iθ₁₄ = −b₁ + ia₁ + s₆

which can be written in the form:

b₁ − ic₁ − a₂e^iθ₁₂ + a₅e^iθ₁₅ + s₄e^iθ₁₄ = 0

b₁ − ia₁ − s₆ + a₅e^iθ₁₅ + a₄e^iθ₁₄ = 0

or as scalar equations:

b₁ − a₂cosθ₁₂ + a₅cosθ₁₅ + s₄cosθ₁₄ = 0

−c₁ − a₂sinθ₁₂ + a₅sinθ₁₅ + s₄sinθ₁₄ = 0

b₁ − s₆ + a₅cosθ₁₅ + a₄cosθ₁₄ = 0

−a₁ + a₅sinθ₁₅ + a₄sinθ₁₄ = 0

which is a set of four equations in four unknowns in the form (input variable θ₁₂ is not considered as an unknown):

F₁(θ₁₄, θ₁₅, s₄) = 0

F₂(θ₁₄, θ₁₅, s₄) = 0

F₃(θ₁₄, θ₁₅, s₆) = 0

F₄(θ₁₄, θ₁₅) = 0

For a given value of input crank angle θ₁₂, we would like to solve the unknown position variables θ₁₄, θ₁₅, s₄, s₆. If θ₁₄⁰, θ₁₅⁰, s₄⁰, s₆⁰, are the initial estimates to the solution, θ₁₄¹, θ₁₅¹, s₄¹, s₆¹ are improved approximations and ds₄= s₄¹ − s₄⁰, ds₆= s₆¹ − s₆⁰, δθ₁₄= θ₁₄¹ − θ₁₄⁰, δθ₁₅ = θ₁₅¹ − θ₁₅⁰ are the difference between the two approximations, then this difference can be solved from equation (7) as:

\displaystyle \left[ {\begin{array}{cccc} {\text{i}{{\text{s}}_{4}}{{\text{e}}^{{\text{i}{{\text{θ}}_{{14}}}^{\text{k}}}}}} & {\text{i}{{\text{a}}_{5}}{{\text{e}}^{{\text{i}{{\text{θ}}_{{15}}}^{\text{k}}}}}} & {{{\text{e}}^{{\text{i}{{\text{θ}}_{{14}}}^{\text{k}}}}}} & 0 \\ {\text{i}{{\text{a}}_{4}}{{\text{e}}^{{\text{i}{{\text{θ}}_{{14}}}^{\text{k}}}}}} & {\text{i}{{\text{a}}_{5}}{{\text{e}}^{{\text{i}{{\text{θ}}_{{15}}}^{\text{k}}}}}} & 0 & {-1} \\ {-\text{i}{{\text{s}}_{4}}{{\text{e}}^{{-\text{i}{{\text{θ}}_{{14}}}^{\text{k}}}}}} & {-\text{i}{{\text{a}}_{5}}{{\text{e}}^{{-\text{i}{{\text{θ}}_{{15}}}^{\text{k}}}}}} & {{{\text{e}}^{{-\text{i}{{\text{θ}}_{{14}}}^{\text{k}}}}}} & 0 \\ {-\text{i}{{\text{a}}_{4}}{{\text{e}}^{{-\text{i}{{\text{θ}}_{{14}}}^{\text{k}}}}}} & {-\text{i}{{\text{a}}_{5}}{{\text{e}}^{{-\text{i}{{\text{θ}}_{{15}}}^{\text{k}}}}}} & 0 & {-1} \end{array}} \right]\left[ {\begin{array}{c} {\text{δ} {{\text{θ}}_{{14}}}} \\ {\text{δ} {{\text{θ}}_{{15}}}} \\ {\text{δ} {{\text{s}}_{4}}} \\ {\text{δ} {{\text{s}}_{6}}} \end{array}} \right]=-\left[ {\begin{array}{c} {{{\text{b}}_{1}}-\text{i}{{\text{c}}_{1}}-{{\text{a}}_{2}}{{\text{e}}^{{\text{i}{{\text{θ}}_{{12}}}^{\text{k}}}}}+{{\text{a}}_{5}}{{\text{e}}^{{\text{i}{{\text{θ}}_{{15}}}^{\text{k}}}}}+{{\text{s}}_{4}}{{\text{e}}^{{\text{i}{{\text{θ}}_{{14}}}^{\text{k}}}}}} \\ {{{\text{b}}_{1}}-\text{i}{{\text{a}}_{1}}-{{\text{s}}_{6}}^{\text{k}}+{{\text{a}}_{5}}{{\text{e}}^{{\text{i}{{\text{θ}}_{{15}}}^{\text{k}}}}}+{{\text{a}}_{4}}{{\text{e}}^{{\text{i}{{\text{θ}}_{{14}}}^{\text{k}}}}}} \\ {{{\text{b}}_{1}}+\text{i}{{\text{c}}_{1}}-{{\text{a}}_{2}}{{\text{e}}^{{-\text{i}{{\text{θ}}_{{12}}}^{\text{k}}}}}+{{\text{a}}_{5}}{{\text{e}}^{{-\text{i}{{\text{θ}}_{{15}}}^{\text{k}}}}}+{{\text{s}}_{4}}{{\text{e}}^{{-\text{i}{{\text{θ}}_{{14}}}^{\text{k}}}}}} \\ {{{\text{b}}_{1}}+\text{i}{{\text{a}}_{1}}-{{\text{s}}_{6}}^{\text{k}}+{{\text{a}}_{5}}{{\text{e}}^{{-\text{i}{{\text{θ}}_{{15}}}^{\text{k}}}}}+{{\text{a}}_{4}}{{\text{e}}^{{-\text{i}{{\text{θ}}_{{14}}}^{\text{k}}}}}} \end{array}} \right]

or in real numbers:

\displaystyle \left[ {\begin{array}{cccc} {-{{\text{s}}_{4}}\sin {{\text{θ}}_{{14}}}^{k}}&{-{{\text{a}}_{5}}\sin {{\text{θ}}_{{15}}}^{\text{k}}}&{\cos {{\text{θ}}_{{14}}}^{\text{k}}}&0 \\ {{{\text{s}}_{4}}\cos {{\text{θ}}_{{14}}}^{\text{k}}}&{{{\text{a}}_{5}}\cos {{\text{θ}}_{{15}}}^{\text{k}}}&{\sin {{\text{θ}}_{{14}}}^{\text{k}}0}&0 \\ {-{{\text{a}}_{4}}\sin {{\text{θ}}_{{14}}}^{\text{k}}}&{{{\text{a}}_{5}}\sin {{\text{θ}}_{{15}}}^{\text{k}}}&0&1 \\ {{{\text{a}}_{4}}\cos {{\text{θ}}_{{14}}}^{\text{k}}}&{{{\text{a}}_{5}}\cos {{\text{θ}}_{{15}}}^{\text{k}}}&0&0 \end{array}} \right] \left[ {\begin{array}{c} {\text{δ} {{\text{θ}}_{{14}}}} \\ {\text{δ} {{\text{θ}}_{{15}}}} \\ {\text{δ} {{\text{s}}_{4}}} \\ {\text{δ} {{\text{s}}_{6}}} \end{array}} \right]=-\left[ {\begin{array}{c} {{{\text{b}}_{1}}-{{\text{a}}_{2}}\cos {{\text{θ}}_{{12}}}^{\text{k}}+{{\text{a}}_{5}}\cos {{\text{θ}}_{{14}}}^{\text{k}}+{{\text{s}}_{4}}\cos {{\text{θ}}_{{14}}}^{\text{k}}} \\ {-{{\text{c}}_{1}}-{{\text{a}}_{2}}\sin {{\text{θ}}_{{12}}}^{\text{k}}+{{\text{a}}_{5}}\sin {{\text{θ}}_{{14}}}^{\text{k}}+{{\text{s}}_{4}}\sin {{\text{θ}}_{{14}}}^{\text{k}}} \\ {{{\text{b}}_{1}}-{{\text{s}}_{6}}^{\text{k}}+{{\text{a}}_{5}}\cos {{\text{θ}}_{{15}}}^{\text{k}}+{{\text{a}}_{4}}\cos {{\text{θ}}_{{14}}}^{\text{k}}} \\ {-{{\text{a}}_{1}}+{{\text{a}}_{5}}\sin {{\text{θ}}_{{15}}}^{\text{k}}+{{\text{a}}_{4}}\sin {{\text{θ}}_{{14}}}^{\text{k}}} \end{array}} \right]

which are in the form:

Aδ = E

where A is the “coefficient matrix”  in either complex or real form. E is the error in the loop closure equations with the initial estimates of the position variables (i.e. if the estimates are correct then E = 0). The solution for the matrix equation is:

δ = A^-1E

After solving for the differences δθ₁₄, δθ₁₅, δs₄, δs₅, the next improved condition is:

s₄^(k+1) = s₄^k + δs₄

s₅^(k+1) = s₅^k + δs₅

θ₁₄^(k+1) = δ₁₄^k + δθ₁₄

θ₁₅^(k+1) = θ₁₅^k + δθ₁₅

We repeat the steps until F(x^(k+1)) < e_e, or x^(k+1) − x^(k) <  e_s. In mechanisms the solution usually converges within less then 10 steps if the initial guess is not too far away from the exact solution.

If the position analysis is to be performed for the whole cycle of the mechanism, one must make an initial guess for the dependent variables for the initial value of the crank angle and then solve for the correct values of the dependent variables using Newton-Raphson method. For the next position, when we increment the crank angle by a certain amount (usually less than 10 degrees), we can use the values of the position variables of the previous step as the initial guess of this step and perform the iterative method for this new position. If the step size is reasonably small, then the solution will converge  to the next position of the mechanism. Since we are solving a nonlinear problem using stepwise linearization, if the step size is too large there is always the possibility of finding the solution for which the mechanism has to be assembled in a different configuration (i.e. if one has started with an open configuration the next solution may give a crossed configuration of the mechanism). Also, you may expect some numerical solution difficulties when the mechanism is at or near critical positions where the Jacobian matrix is zero or very small.

Example:

The six link mechanism shown is used as a rapier drive mechanism in textile industry. Link lengths are:

|C₀A₀|_y = a₁ = 340 mm; |A₀E|_y = b₁ = 520 mm; |C₀A₀|_x = d₁ = 623 mm; |A₀A| = a₂ = 250 mm; |AB| = a₃ = 1225 mm; |CB| = a₄ = 350 mm; |C₀C| = a₅ = 600 mm; |BD| = b₄ = 600 mm and|DE| = a₆ = 1100 mm. You are to determne s₁₆ when θ₁₂ = 60°.

As an initial guess for the position variables assume θ₁₃ = 200°; θ₁₅ = 180°; θ₁₄=80° and s₁₆ = 400 mm when θ₁₂=60°. We shall use MathCAD root finding routine, which utilises Newton-Raphson method (with some modifications). MathCAD output is as follows: You can see that when you write the loop equations and give the initial guess for the position variables, the Given and Find build in functions gives you the result in just a few steps while performing all the necessary iteration.