A novel approximation method for the solution of weakly nonlinear coupled systems

In this paper, we describe and illustrate the application of a novel approximation technique for coupled, nonlinear dynamic systems. The technique begins by obtaining the analytical (or approximate analytical) solutions to the uncoupled system. Then, these solutions are used to approximate particular terms in the fully-coupled, nonlinear system in such a way that the target system is amenable to (approximate) analytical solution algorithms. This work forms part of a larger effort to develop robust control systems for large-scale industrial manipulators. To this end, the final example examined in this work considers the FutureForge manipulator: a state-of-the-art manipulator which forms part of a next-generation forging platform under development at the Advanced Forming Research Centre in Glasgow. To show the breadth of applications of our approach, we also apply it to more widely-recognised models like the Rayleigh and Van der Pol oscillators. In both of these cases, we consider a system of two oscillators each having dynamic behaviour described by Rayleigh/Van der Pol oscillators and coupled together through the resulting damping matrices.


Introduction
In the context of control systems design for highly nonlinear machines with high degrees of coordinate coupling, obtaining a suitably linear model while minimizing uncertainty in the process is challenging. Specifically, in the context of industrial machinery, the dynamics and control of hybrid manipulators are topics which have received much attention in recent years. This paper puts forward a novel approximation technique that helps to linearise and decouple the systems it is applied to. Our intention behind the development of this algorithm is to apply it in the control of an industrial hybrid manipulator with reduced uncertainty by comparison to more commonly-used approaches. Through this work, our approach to achieving insights into this algorithm's strengths, weaknesses, and general applicability is based on empirical observations. In future work, we will endeavour to calcify these observations and extend the algorithm to more generalised problems.
In the existing literature, many authors choose to adopt techniques that render analytical models of the plant unnecessary; or they solve a portion of the model, ensuring that the controller can accommodate the remaining error and uncertainty. Choi et al. elected to consider the nonlinear damping in their manipulator as a disturbance to be compensated by the controller [1]. This was also a consideration of Sun et al. who produced a decentralized controller with a so-called "extended high-gain observer" to help the system withstand a similarly defined disturbance [2]. Dogan et al. and Okur et al. devised a means of expressing the controller's error signal as a function of the Jacobian in the forward dynamics problem [3], [4]. In Dogan et al.'s work, this resulted in the design of a learningbased control system to accurately track the desired end-effector placement. In Okur et al.'s work, this similar approach led to the use of a backstepping control strategy -enabled by the kinematic configuration of their tendon-driven manipulator under scrutiny. This is an approach more recently interpreted by Gul et al. who used a PID controller alongside an adaptive estimator which uses barrier Lyapunov functions to maintain a bound on the error driving the control signal [5]. In an adaptive control context, Guo et al. adopted a simplified model of the manipulator dynamics from Craig et al.'s work in 1987 [6]. Based on a measurement of the rate of parameter variations, constant approximations were substituted in place of the systems' inertia; and, the stiffness and restoring forces were approximated as linear in time [7]. Working with the state-space models, Saab and Ghanem approached this similarly by approximating coupling terms and nonlinearities with constants [8]. Andreev and Peregudova's approach is indicative of many others adopting adaptive control for robotic manipulators whereby they split the inertia, stiffness, and restoring force matrices into known estimates and uncertain parts, thus providing a robust predictor-corrector structuring for the controller [9]. Several other notable works have also treated nonlinearities and coupling terms as "uncertainties" for the control system to adapt in accordance with [10], [11], [12], [13].
Leite and Lizarralde also took this approach of using an adaptive controller to resolve the complexities of a nonlinear, coupled system. Through their cascade controller design, they proposed to control a manipulator with visual feedback and a predictor function based on the "approximate analytical Jacobian" [14]. Friedrich and Buss adopted feedback linearisation as a control approach to cancel the nonlinearities of the governing dynamics, while the coupling was taken care of via their double-Youla approach [15]. Also known as the dual Youla-Kucera approach, this is a means of finding stabilizing parameters for a given plant and control sys-tem. Helwa et al presented an adaptive multiloop control system in which the outer-loop featured a PD controller and a Gaussian Process regression model [16]. With these in the outer-loop and an inverse dynamics model in the inner-loop, the authors demonstrated their approach as applicable "to any Lagrangian system for which [linearisation techniques] can be used to convert the nonlinear dynamics of the system to a set of decoupled double integrators". Meng et al. designed a robust control system while only considering manipulator models with constant inertia, stiffness, and restoring forces by designing a switching function that allowed the controller to select the most appropriate model given the operating point of the structure [17].
This paper is structured in the following manner. Section 2 outlines the algorithm and demonstrates its application to a range of appropriate systems. Section 2.1 deals with the example of a coupled system of Rayleigh oscillators and Sect. 2.2 deals with the example of a coupled system of Van der Pol oscillators. We explain why these systems are considered in advance of Sect. 2.1. We then illustrate the application of our approach to a state-of-the-art hybrid manipulator which draws on findings from previous work [18,19]. In this last example, we compare the result of our approximation technique with that of a method considered by some of the authors cited in our literature review. Finally, Sect. 3 presents the main conclusions of this work.

Development of the method with example applications
This section outlines our proposed method for obtaining approximate analytical solutions for nonlinear coupled dynamics problems. Each example explores the application of the technique to a different dynamic system. In general terms, the algorithm distils to this: use solutions to the uncoupled system behaviour to approximate terms in the fully-coupled problem as necessary to render it analytically solvable. In the subsequent conclusion, we will note general observations on this approach. Derived initially to describe the acoustic vibration of a clarinet reed, the Rayleigh oscillator has since been used in many applications [20] ranging from robotic bipedal motion [21] to oscillatory phenomena observed in chemical reactions [22]. Given its breadth of applications and the depth of analysis devoted to it, we use this as an illustrative case to demonstrate our approximation technique. For similar reasons we will also consider the Van der Pol oscillator as a further example. These two oscillator models can be derived from one another [23] and they both feature heavily in the study of limit cycle motion and self-sustained vibrations. The Van der Pol oscillator has also appeared in the modelling of vortex shedding patterns in the wake regions of slender bluff bodies [24] as well as in the study of voltage oscillations in electrical circuits (the original context which led to its derivation) [23].
A key feature of these oscillators is their relatively weak damping properties (for small values of ε < 1) which can either feed energy to the systems or dissipate energy from them according to their nonlinear damping coefficient [25]. As will be seen in the examples that follow the accuracy of our approximation technique is partially informed by the accuracy of an approximate analytical solution of the uncoupled problems considered. To explain briefly: in the case of the Van der Pol oscillator system we use the Lindstedt-Poincaré approach to solve a system of uncoupled oscillators first, and then use this approximate solution to inform our approximation of the coupled system. Thus, the low damping exhibited in the Rayleigh and Van der Pol oscillators allows these to serve as good illustrative examples for our purposes.
Both the Rayleigh and Van der Pol systems studied here have relatively weak nonlinear and coupling effects due to the values of the parameter ε in both cases. Discussion of these problems then only gives us insight into similar problems: those where the nonlinear and coupled damping effects are weak relative to the inertia and restoring effects. The FutureForge manipulator system presented in Sect. 2.3 is not one of these systems and so presents an extension of the method's observed applicability. In future work, we will look to expand the range of systems considered. However, for the problems in this current work, the results presented here are similar for all those with ε < 0.5.

Problem specification
In this section, we consider the following system which is in an analogous form to the Rayleigh Differential Equation. This can be written in the simplified form where P(Ẏ ; ε) is the damping matrix including the small parameter ε.
This method aims to find suitably accurate approximate analytical solutions for y 1 and y 2 . Our simulations of this system will involve the following generalised initial conditions where a, b, c, and d are constants.

Solution of the uncoupled system
First, we define the uncoupled system by where we have used φ k = φ k (t) to denote the uncoupled equivalents of y k = y k (t). By approaching this problem using a technique such as the Lindstedt-Poincaré method, we can obtain asymptotic approximations of its solutions where τ k is the stretched coordinate, n is the same as in Eq. (2.6), and the constant values of ω j are selected so as to remove secularity from the resulting solutions. By substituting Eq. (2.6) and its first two time-derivatives into the uncoupled system (2.5), we can obtain the O(ε 0 ) perturbation equations for φ 1 (t), and for φ 2 (t), with higher-order perturbation equations being gathered in the usual fashion. Solving these gives us analytical expressions for the asymptotic approximations of φ 1 and φ 2 .

Approximation and solution of the fully-coupled system
Our approach proposes that we now substitute the uncoupled solutions φ k into the fully-coupled system in a way that renders it solvable through analytical means. One way of doing this would be to substitute φ 1 and φ 2 in place of y 1 and y 2 respectively (and their derivatives) This simplifies the original system (2.2) into a system of two uncoupled simple harmonic oscillators with excitations applied to both.
Finding the complementary functions for this system is very simple. For generality, we use the forms and where ψ 1 (t) = ψ 1 and ψ 2 (t) = ψ 2 are the bases of the solutions. Using ψ 1 = cos(t) and ψ 2 = sin(t), we can use the method of variation of parameters to compute the general form of the particular integral. Finding the Wronskian we can find the particular integrals as and By computing these, and summing them with their corresponding complementary function, we obtain the general solutions of the system. We can obtain the particular solutions by applying the generalised initial conditions.

Implementation and simulation
To illustrate the method and demonstrate its practical implementation, we consider the scenario of with ε = 0.01, and we will use only the O(ε 0 ) asymptotic approximation to determine our coupled approximation. The numerical solution, which we will use as our reference for the success (or otherwise) of our approach, yields the following plots of the generalised coordinates ( Fig. 1). Solving the uncoupled system (2.5) by means of the perturbation method described previously, under the conditions in (2.17), we obtain the generating solutions and and By computing the particular integrals and summing them with the corresponding complementary functions, we arrive at as the general solution for y 1 , and as the general solution for y 2 . By applying the initial conditions to these two general solutions and their first derivatives, we arrive at the following particular solutions.  of the uncoupled problem. Although these errors are clearly produced by different techniques and different problems, this graph illustrates that the accuracy of our method degrades at a slightly faster rate than this widely-known method for the uncoupled problem.

Problem specification
We now consider the case of a coupled system of two Van der Pol oscillators.
We note that this can be generalised into the same form as the system of Rayleigh oscillators in Sect. 2.1, Eq. (2.2), although, this system's damping matrix is dependent on Y rather thanẎ .
For generality, we use the same notation for the generalised initial conditions as were used in Sect. 2.1.

Solution of the uncoupled system
The uncoupled system, in this example, is very similar to the previous example.
Despite this slight difference to the previous uncoupled system, Eq (2.5), the perturbation method is identical and leads to identical uncoupled solutions at the O(ε 0 )level. This is because the difference between the two systems is limited to the stiffness matrix, P, which is O(ε 1 ) in its contribution to the overall system dynamics. Thus, the generating problems for these solutions are Eqs. (2.8) and (2.9).

Approximation and solution of the fully-coupled system
While the definition of the damping, P(Y ; ε), has no impact on the uncoupled solution by comparison with the previous example, it contributes significantly to the form of our non-homogeneous term in the approximation of the coupled system.
Using this, the coupled system is approximated by a simple harmonic oscillator with an excitation described by Eq. (2.11).

Implementation and simulation
To illustrate the implementation in this scenario, we use the same initial conditions and ε-value as in Sect. 2.1.
The reason for doing this is that it gives us an insight into the impact of the differing complexities of the stiffness matrix.
As should be expected, Fig. 2a and b also represent the uncoupled solution of the Van der Pol system discussed in this section. However, there is a significant reduction in the approximation error when we consider the approximate analytical solutions produced by our method. As can be seen in Fig. 5b, the maximum error is less than 2.5% for this example while it was close to 7% in the previous example (see Fig. 3b).

The impact of higher-order approximations of the uncoupled system
An intuitive means of increasing the accuracy of the approximations of the uncoupled system may be to increase the accuracy of the solutions to the uncoupled system. This idea has some intuitive support. Finding the O(ε 0 ) and O(ε 1 ) uncoupled solutions and then substituting both of these into the fully-coupled system, as we have shown previously, we can see the following reduction in the relative errors (Fig. 6).
There is clearly some small benefit in using the O(ε 1 ) approximations for the uncoupled problems. However, there is a significant increase in the computational complexity of then attempting to solve the approximation of the coupled system analytically. Given the motivation behind using this approach in the first instance, we don't find a significant advantage in using the O(ε 1 ) approximation despite its minor improvement in accuracy.
In Fig. 7a and b we illustrate the error growth arising from our approximation in comparison to that arising from the standard Lindstedt-Poincaré solution of the uncoupled Van der Pol oscillator. Notably, our approximation retains higher accuracy than the O(ε 0 ) uncoupled solution until t = 10 for y 1 (t) and around t = 5 for y 2 (t). Additionally, this error growth indicates the presence of secular terms in our approximate solution which do not occur in the Lindtedt-Poincaré approach. In further work, we will look to examine a means of accounting for such behaviour appropriately.

Problem specification
The FutureForge manipulator is a large-scale industrial manipulator designed to manoeuvre metallic workpieces in, out, and around a state-of-the-art forging environment. In particular, one of the critical areas of performance for this manipulator should be its ability to place accurately and manoeuvre metallic specimens through treatment in the hydraulic press. The wider scope of our work is to develop a control system for this manipulator (Fig. 8).
Considering the structure of this manipulator, rather than fitting into the traditional categories of either series or parallel manipulators, this is a hybrid structure as a series of parallelograms. The advantage to this is that it has the enhanced workspace reach of a serial manipulator with the rigidity of a parallel manipulator.

Solution of the uncoupled system
In [18] the solution of uncoupled actuation scenarios of this manipulator was investigated. Between this and [19], it is shown that the system can be written as This schematic is also shown previously in [19] and in [18]  where F v and F h are the remainders of the excitation terms when the components counteracting gravitational restoring forces are removed. It was then subsequently shown, in [18], that an approximate analytical solution of this system can be obtained via a perturbation method and introducing a small parameter, ε, via its nonlinear terms. The O(ε 0 ) problems for this method were given to be in which α 0 and β 0 are the leading-order terms in the asymptotic approximations of α and β. These can be solved to give the O(ε 0 ) approximate solutions in the forms where the K -terms are constants determined by the specific problem parameters. In the general case, these can be written as follows where u v0 and v v0 are the initial position and velocity conditions respectively for Eq. (2.30); u h0 and v h0 are the initial position and velocity conditions respectively for Eq. (2.31); and, F v and F h are the actuating torques in these two problems.

Approximation and solution of the fully-coupled system
To approximate the fully-coupled system, we use a similar approach as in the previous two examples of the Van der Pol and Rayleigh oscillators. In one of the system's equations, we substitute α for the approximate uncoupled solution Eq. (2.39). We denote these approximations α(t) ≈ f (t) and β(t) ≈ g(t). The system can then be written as which is in a similar form to the uncoupled actuation problem. For readability through the rest of this work, we denote (2.51) The key difference between this system and the uncoupled one is the structure of the excitation terms which can present a challenge for analytical solution methods. As a result, we further approximate the system by finding Taylor Series expansions of the excitation terms. With regards to the order of expansions required, we notice that higher-order terms do not significantly increase computational complexity for these cases (merely the number of terms that must be computed) since the difficulty that they allow us to overcome is the integration of the transcendental terms in the excitations of our approximated model. In practice, the order of expansion used here is at the practitioners' discretion. Through trial and error, we notice an acceptable loss of accuracy when using a 10 th -order Taylor Series approximation in both equations as will be illustrated in the next section. In the following equations, we use the function h(·) to indicate that (·) has been approximated in this way.
With the system approximated in this manner, we proceed with an analytical solution via the method of variation of parameters. The nonlinear terms occurring on the left-hand sides of Eqs. (2.49) and (2.50), i.e. on the diagonals of M and P, mean we can benefit from reusing a perturbation method as in the uncoupled case. As such, we introduce the small perturbation parameter via a 2 = εã 2 , a 6 = εã 6 , b 2 = εb 2 , and b 6 = εb 6 to give a 1α + a 5α + ε ã 2 cos(2α)α +ã 6 sin(2α)α 2 As such, when we introduce the generalised asymptotic approximations for any positive integer value of n, and substitute them into Eqs. (2.52) and (2.52), we find the following perturbation equations at the generating level.
Solving these equations analytically is not overly troublesome with determining the particular integral made simple via the method of variation of parameters. Due to the high order of Taylor Series approximation used for h, we refrain from writing the full form of the solution explicitly in this work.

Implementation and simulation
In testing the algorithm, we undertake a range of manoeuvres as detailed in Table 1. These particular manoeuvres and torques are chosen because they allow us to observe how our approximation performs with different excitations and different initial configurations. For the manoeuvres with durations of t = 2 seconds, we will also inspect the shorter timescale performance from the same set of results as necessary.
Throughout this section, we compare our approximate analytical solution to one that is commonly used in the control of manipulators such as this one. This method essentially regards all nonlinearities as disturbances/uncertainties in the system model. So, we can represent it here as the solutions to the uncoupled problem in Eq. (2.36).
For Case 1, Fig. 9a-d, the error of our approach is initially lower than that of the uncoupled approximation before it grows rapidly and extremely high around t = 1 second in both the vertical and horizontal case. In Fig. 10a and b, we restrict the displayed time interval as well as the percentage error to highlight the comparative performance over shorter timescales. For Case 2, Fig. 11a-d, we see very similar behaviour to that observed in Case 1. The difference between these scenarios is the value of the excitation torque F h . With the lower value of F h in Case 2, we observe that the rapid growth of error in our approximation of α is delayed slightly: more notably occurring just after t = 1 (comparing Figs. 10a and 11b). Additionally, in this case, the error growth in β seems slower before it eventually surpasses that of the uncoupled approximation at approximately the same time as in the previous case (comparing Figs. 10b and 12).
For Case 3 and those that follow, we default to displaying a restricted interval in the time domain. We do so with the proviso that the rapid growth in error seen in Cases 1 and 2 can also be observed similarly in all cases. The error in this simulation of α is seen to be lower than in Case 1 and similar to that in Case 2 (over the restricted domain in t). Equally, the error in the Case 3 simulation of β is seen to be lower than in Case 1 and slightly lower than that in Case 2 (Fig. 13).
In Case 4, the errors computed for both α and β are lower than in Cases 1-3 (Fig. 14).
In Case 5, we see our approach perform more desirably than for Cases 1 and 2. However, it exhibits a significant loss in accuracy in comparison to Case 4 (Figs. 15,16,17).
In the interest of comparing the accuracy of our approach versus the uncoupled approximation in Cases 1-7, we show Fig. 18a and b. The curves on these graphs are found by subtracting the error of our approach from the error of the uncoupled approximation. Therefore, when a curve is greater than 0 on the vertical axis, our approximation is superior at that point.
For interest, we highlight that Case 2 is the only simulation (of these seven) that sees our approximation initially produce a less accurate result. The curve for this particular simulation, however, soon surpasses the reference approach. In Cases 1-7, for α(t), Case 3 exhibits the greatest improvement in accuracy while Case 4 shows the longest retention of its accuracy. Case 7 appears to show the lowest improvement in the accuracy for this set of simulations. And, Case 6 is the first to become less valuable than the reference approach. For β(t), Case 2 shows the highest improvement in accuracy (although closely followed by Case 1) while Case 7 is the only result in this set of simulations (for either coordinate) to show no improvement by comparison to the reference method. Additionally, for β(t), Case 6 retains superior accuracy for the longest duration by contrast to the aforementioned Case 7.
In Cases 8-7, we see some similarities in comparison to the previous cases. However, there are also some key differences which will be reflected upon in the subsequent section.
In this set of simulations, for α(t), Case 12 shows the greatest improvement in accuracy while Case 14 shows the least. Case 13 retains its superiority for the longest duration, and Case 9 is the first to degrade below the reference method's accuracy. Note that Case 14 initially dips below the 0% improvement line, but rises up beyond this and retains its benefit for longer than Cases 9 or 11. Regarding β(t) in Cases 8-14, Case 12 exhibits the greatest improvement in accuracy as Case 13 exhibits the least (showing largely inferior performance to the uncoupled approximation). Additionally, it is also Case 12 that retains its accuracy for the longest duration while Case 13 is also the first to degrade.

Discussion regarding the FutureForge manipulator simulations
In the simulations detailed in Table 1, we have explored the following question. Is the accuracy of our approach affected by changing actuation forces, initial conditions, and/or manoeuvre durations? Clearly, yes. A question remains though as to how precisely these changes impact the accuracy of our approximations.
In the first set of simulations, Cases 1-7, the greatest improvement offered by our approximation was found in Cases 1-4. In these cases, we are varying the excitations applied to the system while maintaining the initial conditions across them. Cases 5-7, in which we varied the initial conditions while maintaining the excitations, generally performed worse than Cases 1-4. That being said, a small improvement in accuracy was still exhibited in these results. In the second set of simulations, Cases 8-14, the improvement offered by our approximation technique was generally lower than in Cases 1-7 although still present. Inevitably, the accuracy of our method will stem from how reasonable it is to use Fig. 21 A comparison of our approximation method with the uncoupled approximation (used as a measure of success in the FutureForge manipulator example) applied to a the Rayleigh oscillator system and b the Van der Pol oscillator system the uncoupled system dynamics to approximate aspects of the coupled system's behaviour.
In the context of the control of hybrid manipulators such as this, accurate system prediction over short timescales can be arguably more beneficial than accurate prediction over longer durations. In model predictive control systems, the prediction horizon is the parameter specifying how far into the future the prediction model forecasts. If this is less than 1 second, then our approximation technique is superior to the standard uncoupled approximation as it accounts for a higher proportion of model uncertainty. However, for prediction horizons greater than this timescale, preference may still lie in the uncoupled approach.

Conclusions
In this work, we have described a novel method for the approximation and resulting analytical solution of coupled, nonlinear dynamic systems. We have also illus-trated its application in three different contexts. The first of these was a system of coupled Rayleigh oscillators; the second was a system of coupled Van der Pol oscillators; and, the third was a next-generation industrial hybrid manipulator. In each of these contexts, we showed the application of our approximation technique before showing the advantage offered by this approach.
In the cases of the Rayleigh and Van der Pol systems, our method performs well by comparison to the O(ε 0 ) solution of the uncoupled problems (see Figs. 4, 7a, b). For the manipulator example, we compared our approach to a widely-used approximation in the context of predictive control systems: considering all nonlinearities and coupling effects as uncertainties. If we compare our approximations of the Rayleigh and Van der Pol systems to this same technique, then we see our method outperforms this reference approach in both cases (see Fig. 21a, b).
By comparison to the uncoupled approximation of the system, we see that our method generally achieves greater accuracy over an initial interval in the time domain before becoming inferior after this interval is exceeded. This seems to be a reasonable trade-off: our method is not a universal cure-all for the approximation and analytical solution of coupled, nonlinear systems. Rather, it offers practitioners (particularly in the field of control) an alternative means of approximation which has enhanced accuracy in the short term. We appreciate that there is much work to be carried out to fully understand the error growth, optimal means of application, and potential pitfalls of this approximation technique.

Comments for further work
In future work, we will explore the generalisation of this method in two senses. Firstly, we will conduct a parameter study for the Rayleigh and Van der Pol systems described here. This is intended to give an understanding of the parameter space in which our approach is applicable and affords some benefit to the user. Secondly, we will look to extend this-from a parameter study based on the specific problems studied in this paper-to a study of the parameter space of more generalised ordinary differential equations. Additionally, we look to explore the application of this technique in the context of designing a control system for the FutureForge manipulator. However, for now, we hope that this broad presentation of the concept is interesting and finds some utility in the realm of the linearisation of nonlinear systems.