Maxwell Points of Dynamical Control Systems Based on Vertical Rolling Disc—Numerical Solutions

: We study two nilpotent afﬁne control systems derived from the dynamic and control of a vertical rolling disc that is a simpliﬁcation of a differential drive wheeled mobile robot. For both systems, their controllable Lie algebras are calculated and optimal control problems are formulated, and their Hamiltonian systems of ODEs are derived using the Pontryagin maximum principle. These optimal control problems completely determine the energetically optimal trajectories between two states. Then, a novel numerical algorithm based on optimisation for ﬁnding the Maxwell points is presented and tested on these control systems. The results show that the use of such numerical methods can be beneﬁcial in cases where common analytical approaches fail or are impractical.


Introduction
The apparatus of Lie groups is commonly used in the control theory of planar nonholonomics mechanisms, such as kinematic cars, robotic snakes, etc. [1][2][3][4][5][6][7][8]. The control Lie group derives a control Lie algebra as a tangent linear space of the group with operation Lie bracket. The main advantage of the control on Lie groups is that everything derived for one point of the mechanism's configuration space can be easily generalised for all points of the configuration space thanks to the group properties [9][10][11]. For computational reasons, it is preferred that the corresponding Lie algebra be nilpotent. This property generally does not apply, but the Lie algebra can be approximated, usually by an algorithm called a nilpotent approximation or Bellaïche's algorithm [12,13]. This paper is motivated by [14], where the authors study the Maxwell points of a kinematics model of vertical rolling disk using numerical optimisation methods. The Maxwell point concept is related to optimal control problems [11,15]. At the initial point, the local solution of optimal control problems is the optimal trajectory (referred to as a geodesic), that goes from the initial point and the corresponding control. Maxwell points are points where the geodesics intersect each other for the first time and the corresponding geodesic segments have equal length. For example, Maxwell points on the Heisenberg group lie on the line (see [16]). This implies that Maxwell points are points where geodesics lose optimality. In optimal control problems, there can be other points at which geodesics lose their optimality, for example, conjugate points (see [11,15]).
The problem of existence and finding Maxwell points is an open problem [11,17]. There are only a few concrete situations where the problem is solvable analytically (e.g., the kinematics of a vertical rolling disc) [16,18]. The main goal of this paper is to design a numerical algorithm to find the Maxwell points of the dynamics of a vertical rolling disc. The rolling disk serves as a case study for the methods under development to be tested on, but it may be applicable to a much larger group of complex nonlinear dynamical systems, real-life robots, or may include considerations or models of imperfections of a real dynamical system, as demonstrated in [19].
Using numerical approaches to solve different nonlinear tasks was proven useful in [20,21]. In Section 2 we describe a vertical rolling disc, including the derivation of differential kinematics and dynamics. In particular, we derive the Lie algebra for the differential kinematics and we indicate the shape of the Lie algebra for the dynamics.
In Sections 3 and 4 we study two different approximations of the Lie algebra for the dynamics of the disc. Later, in Section 3, we use the Taylor polynomial and in Section 4 we use Bellaïche's algorithm for the approximation of nonlinear expressions, while both approximations give a nilpotent Lie algebra. For both approximated dynamic systems, the optimal control problem is formulated, where the goal is to get an energetically optimal trajectory between two points of the configuration space. We use the Pontryagin maximum principle to derive the ODEs whose solutions give the optimal trajectories.
In Section 5 we show the description of the numerical algorithm developed for finding the Maxwell points of the system. The input of the algorithm is the set of ODEs given by the Pontryagin maximum principle.
Finally, in Section 6 we test the algorithm on both approximated dynamics systems derived in Sections 3 and 4.

Description of the Mechanism
This paper discusses a differential drive wheeled mobile robot's motion where there are two independent motors and one/two auxiliary wheels (like TurtleBot 2 (Figure 1b), Khepera robot, etc.). Let the mechanism be the homogeneous rolling disc of mass m = 1 and radius r = 1 in plane xy. Assume that the disc moves without tilting and slipping. It can change the acceleration in the forward or backward direction and the angular acceleration in the vertical axis. The position of the disc can be described by the vector q = (x, y, θ, v, ω) T , where θ is its orientation, v its velocity in the forward direction, and ω the angular velocity ( Figure 1a). It gives 5-dimensional configuration space C = {q = (x, y, θ, v, ω) T }.

Differential Kinematics of the Mechanism
In this section we want to describe differential kinematics of the mechanism [22][23][24]. The non-slip condition gives us following equation, referred to as the Pfaffian constraint: This equation holds ifẋ = u 1 cos θ, where u 1 = v, u 2 = ω are control parameters for the kinematics. Note that any movement back and forth can be expressed byẋ = u 1 cos θ,ẏ = u 1 sin θ, andθ = 0 and any change of the orientation can be expressed asẋ = 0,ẏ = 0, andθ = u 2 . It is obvious that we obtain two vector fields V 1 = cos θ∂ x + sin θ∂ y , yielding a first-order control systemq where q = (x, y, θ) T . The vector fields V 1 and V 2 generate 3-dimensional Lie algebra, where The operation Lie bracket on the algebra is shown in multiplicative Table 1. Table 1. Lie bracket of Lie algebra for differential kinematics.

Dynamics of the Mechanism
The dynamics of the mechanism with respect to movement without tilting and slipping is m is the mass of the disk, I = 1 4 mr 2 = 1 4 is the moment of inertia about the vertical axis, r is the radius, F is the force that causes forward and backward movement, M 2 is the moment that causes the rotation around the vertical axis, and λ is the friction force [25].
For the first derivation of the position the following holds: v is the velocity in the direction of the movement and ω is the angular velocity around the vertical axis. It givesq =ṠV + SV, and it can be derived thatẋ = v cos θ, The configuration space is C, dim C = 5. The system (5) can be written aṡ where u 1 = M 1 , u 2 = M 2 are control parameters for the dynamics. The vector fields V 0 , V 1 , V 2 generate infinitely dimensional Lie algebra, so for a discussion of controllability, the vector fields need to be approximated appropriately.

Approximation of Vector Fields Using Taylor Polynomial
An approximation by the Taylor polynomial is offered as one of the possibilities. Using a Taylor polynomial of the first order at point θ = 0 we get cos θ ≈ 1, sin θ ≈ θ. Then the dynamics of the mechanism iṡ The vector fieldsṼ 0 ,Ṽ 1 ,Ṽ 2 generate 8-dimensional nilpotent Lie algebra, wherẽ The operation Lie bracket on the algebra is shown in multiplicative Table 2. Table 2. Lie algebra of approximated vector fields using Taylor polynomial.
We can see that the Lie algebra is nilpotent but the dimension of the algebra is higher than the dimension of the configuration space C. However, at point q 0 = (0, 0, 0, 0, 0) T we get nilpotent sub-algebra of dimension 5, so the mechanism is locally controllable. In Section 4 we compute the approximation of vector fields V 0 , V 1 , V 3 using Bellaïche's algorithm to get nilpotent algebra of dimension 5.

Optimal Control Problem
The dynamics (7) with conditions and optimality condition gives the optimal control problem whose solution gives an energetically optimal trajectory between two points q 0 , q T ∈ C, where u 1 , u 2 are control parameters and T is the end time.
The Hamiltonian of the optimal control problem is Using PMP we get the optimality conditions ∂H ∂u 1 = 0, ∂H ∂u 2 = 0. This gives so the Hamiltonian is in the form The adjoint system and the dynamics of the mechanism (7) with respect to (10) gives the following system of differential equations: where The solution of the system with respect to (8) gives the optimal curve between points q 0 , q T ∈ C. The corresponding control u 1 , u 2 can be found from (10).

Nilpotent Approximation Using Bellaïche's Algorithm
In this section we would like to find 5-dimensional nilpotent algebra that approximates the algebra generated by vector fields V 0 , V 1 , V 2 derived in (6). Let us introduce some notions. Define ∆ 1 to be the linear subspace generated by X 1 , . . . , X m , In the construction of a nilpotent approximation of the given Lie algebra we will first proceed with the construction of privileged coordinates according to Bellaïche's algorithm [12].

1.
Choose an adapted frame Y 1 , . . . , Y n at p. In our case we choose frame (Y 1 : Choose coordinates (y 1 , . . . , y n ) centered at p such that ∂ y i | p = Y i (p). For our adapted frame we choose the coordinates with |α| = α 1 + · · · + α n . First notice that for coordinates of degree w i < 3 the transformation z i = h(y 1 , . . . , y i ) degenerates [12]. Thus, we will focus on the evaluation of z 5 . The sum yields We can see that for the choice of θ = 0 the transformation degenerates, and hence from now on let us introduce a substitution C = 1 4 sin(θ)| p . Now, express the base of tangent bundle of old coordinates in terms of new coordinates Finally, let us express original vector fields in algebraic privileged coordinates.
The second step of nilpotent approximation is the approximation of the vector fields' elements (from a function viewpoint). First let us define a weighted degree of monomial. Given a sequence of integers α = (α 1 , . . . , α n ), we define the weighted degree of the monomial z α = z α 1 1 . . . z α n n to be w(α) = w 1 α 1 + . . . + w n α n and the weighted degree of the monomial vector field z α ∂ z j to be w(α) − w j .
For vector field X with a Taylor expansion the order of X is the least weighted degree of a monomial vector field having a nonzero coefficient in the Taylor series. Grouping together the monomial vector fields of the same weighted degree we can express each vector field as a series of monomial vector fields. A monomial vector field with the least weighted degree gives us nilpotent approximation of the Lie algebra [13], thus The operation Lie bracket on the algebra is shown in multiplicative Table 3. Table 3. Lie bracket of nilpotent Lie algebra found by Bellaïche's algorithm.
The dynamics of the mechanisṁ and conditions (8) and (9) define a different approximation of the same optimal control problem as in Section 3. Like in Section 3 we get that for the control the following hold: The Hamiltonian is in the form and the system of differential equations is in the forṁ The solution of the system with respect to (8) gives the optimal curve between points q 0 , q T ∈ C.

Numerical Analysis
As stated in the previous sections, obtaining the Maxwell points of the system is crucial if our goal is to solve the optimal control problem. However, finding them analytically is not always possible, and even if an analytical solution is technically reachable it may not be a practically feasible one.
In practice, a viable alternative to analytical methods is a numerical solution to the problem. Therefore, we developed a numerical apparatus to help discover the approximate location of these Maxwell points of a dynamical system, and to be able to validate any Maxwell point candidates.
The numerical approach proposed throughout this section utilises the Maxwell point definition in the form: a Maxwell point (MP) is any state of an autonomous dynamical system such that it is the end point of (at least) two different trajectories equal in both length in the state space and time, starting at any given state. Generally, there is no guarantee that Maxwell points exist in any given case. However, when presented with a specific system, we can attempt to find them using this definition, and if the attempt fails it can be presumed that the Maxwell points do not exist.
To describe the methods used in further subsections a basic nomenclature will be described here. First, we assume that we deal with an autonomous system of ODEs-the control system, which was created by expanding the original system using the PMP as described in Sections 3.1 and 4.1, resulting in the explicit form (15) where y ∈ R 2n is the state space of the control system, q ∈ R n is the state space of the original system, and h ∈ R n is the extended space generated by the PMP with respect to a minimisation functional L(y, t) described by Equation (9). These spaces are related as Additionally, we denote a specific state at a given time instant as q(t) = q t , for example, q(0) = q 0 representing the initial state or q(T) = q T the final state, if t ∈ (0; T). If more than one state trajectory is being described, the respective states are marked q i (t) = q it , with i being the trajectory index. The same notation is used for vectors y, q, and h.
Given this formulation, the task of looking for a Maxwell point in the q space corresponds to solving the boundary value problem for a system described by (15) with the boundary values given by any q 0 and q T = q M representing the Maxwell point. Note that we only trace the position of Maxwell points in the state space of the original system q, not in the whole state space of the control system y.
Further, given an initial condition y 0 = [q T 0 , h T 0 ] = [q 10 . . . h n0 ] T we can simulate the system numerically and produce a trajectory using numerical integration methods such as Runge-Kutta [26]. This allows us to acquire trajectory end points. Different system trajectories can be generated using different initial conditions. These properties allow us to effectively convert the BVP solution, which is generally a very difficult task, to an IVP coupled with an optimisation algorithm (commonly referred to as the shooting method) [27].
As stated earlier, if a point is a Maxwell point of the system, then two different trajectories from the same initial point q 0 will intersect at time T at this point. Based on this assumption we can now reformulate our problem as a shooting problem, where we have a fixed initial state of the system q 0 and we are looking for two different vectors h 10 and h 20 which generate trajectories leading to the same point q 1T = q 2T .
Finally, viewing the trajectory end point position in the state space as a function of its initial condition q T = q T (q 0 , h 0 ) will allow us to formulate an optimisation task for both finding and validating the MPs. To solve these tasks we implement an interior point search algorithm. This constrained nonlinear optimisation algorithm is used to solve problems of the form arg min where J(β) is the cost function, φ(β) is the equality constraint function (which is meant to match the length of multiple trajectories if necessary), γ(β) is the inequality constraint function (which is meant to penalise multiple similar trajectories), and β is the vector of parameters that are being optimised. The optimisation algorithm itself is well-described in the literature (e.g., [28]), and we employ it using a standardised implementation in the Matlab environment [29]. The optimisation tasks utilised in later sections follow the same formulation.

Finding Trajectories Leading to Known MPs
If the position of a Maxwell point is known we may want to verify that this point is truly the Maxwell point of the system. This may seem counter-intuitive but proves to be useful when we only assume the property and need to validate it. This method will be used in the following sections to validate our assumptions made about a subset of the state space which we had assumed the MPs to lay on.
As per our previous definition, assuming a known position of an MP, we need to find two trajectories of the same length which intersect at the point at time instant T. We can look for the trajectories one at a time, adding constraints to prevent the optimisation algorithm from converging to the same solution repeatedly.

1.
Find the first trajectory Assuming the MP position is invariant to the actual trajectory length in time, we can fix the initial configuration of the system to the origin (q 10 = 0) and choose a final time T ∈ R + as any positive real number (affecting only how fast the system passes through the optimal trajectory). We can now optimise the end point of the trajectory to be equal to the suspected Maxwell point q 1T ≈ q M . This means we choose such h 10 that drives the system to q M in time T. This is a straightforward convex minimisation task deduced from (17), since q 1T now depends only on the choice of h 10 and no constraints are required. The cost function will take the form of a euclidean distance result in an optimisation task defined as where h * 10 is the optimal solution. The whole situation is depicted in Figure 2. Figure 2. Illustration of the iterations for the optimisation task described by Equation (18). q 0 depicts the initial condition, q M the assumed Maxwell point, and r distance to the actual iteration final point q T .

Find the second trajectory
Now we need to find h * 20 for a second trajectory, which is different from h * 10 and drives the system towards the same point q M , as shown in Figure 3. We can formulate our constrained minimisation problem as (19) with constraints (20) and (21). The equality constraint (20) is meant to satisfy the same functional (9) for both trajectories.
The parameter ρ is a tuning parameter which specifies the minimal distance between the vectors h 10 and h 20 , which are normalized so that the choice of this parameter is independent of their respective scales.

Review results
There are multiple termination conditions to the interior point algorithm. These include: • First-order optimality; • Iteration step size; • Number of iterations.
Any of these can cause the algorithm to stop, but the main concern is the case where the termination occurs due to a local minimum (first-order optimality reached). If this happens, we may need to return to step 2 and restart the algorithm, giving it a different initial guess.
To decide whether our result is valid, or if we need to restart the optimisation with a different initial guess, we can compare the value of the optimisation cost function J(h 20 ) with a numerical tolerance parameter. If the value is smaller than our chosen numerical tolerance, the result is valid and the Maxwell point has been confirmed.

Finding MPs
In the previous section we described how to verify assumed Maxwell points, now we will describe how to find them using a similar approach. We experimented with various combinations of optimisation and random shooting and found that typical approaches such as Monte Carlo or grid search are not viable in this case because it is quite difficult to "hit" the Maxwell point exactly using shooting from random initial guesses. To overcome this issue, the task needs to be formulated as a convex optimisation problem which can converge with arbitrary precision if left iterating.
The proposed process itself is very similar to the one described in the previous section. However, in this case the location of the end point where the two trajectories intersect is unknown, so we need to optimise both h 10 and h 20 at the same time to make the end points q 1T and q 2T converge to a single point. This means that the dimensionality of the problem increases twofold because the number of unknown parameters is 2n instead of n. The whole process is shown in Figure 4, and can be described as looking for two different trajectories which start and end at the same points and have the same lengths. Based on the MP definition described earlier, if the end points do converge, we can take the end point for a Maxwell point. Formally, the optimisation task is then defined as (22), also with the constraints (20) and (21). Similar to the previous section, we also need to perform an evaluation of the results through the cost function value and check whether the found minimum is truly the Maxwell point or if it is a local minimum.
r c(y 1 (t);y 2 (t)) After a Maxwell point is found we can restart the method with a different initial guess and repeat the process to obtain more Maxwell points. Initial guesses can be chosen randomly. It is important to note that the search process appears stochastic, and not all initial guesses are guaranteed to converge to an MP, even if they do exist. Naturally, the more initial guesses that are tried, the higher the chance of finding an MP.
When a sufficient number of Maxwell points is found, we can formulate a hypothesis about the set they form. Generally, the MPs form a subset of the full space, often a hyperplane in linear cases. The form of the set can be analysed using standard data analysis tools such as principal component analysis or proper orthogonal decomposition (PCA) [30]. Afterwards, our hypothesis can be verified using the algorithm from the previous section applied on a new MP predicted by the hypothesis.

Numerical Experiments
In Sections 3 and 4 we described the optimal control problems developed from two different approximations of the rolling disc dynamics, and in Section 5 we proposed numerical tools based on previous research [14] for analysing these control problems and finding the Maxwell points of the systems. In this section, we apply the numerical tools on both systems as case studies.

Rolling Disc-Nilpotent Approximation with Drift
Using the approximation from Section 4, based on (14), in the form of (15), the optimal control problem gives us the system of ODEs of the forṁ This is the starting point of our numerical algorithm. We can use one more simplification based on our system of equations. The optimality functional (9) is dependent on the optimal control inputs (13). However, by calculating this we see that the resulting value of the functional is only dependent on h 1 , h 2 , which are decoupled states (states which are not dependent on any other states, only on time). This means that it is not necessary to compute the value of the resulting integral (9). If the norm of our states h 10 , h 20 is the same for both trajectories, the resulting integral value will be the same. This simplifies the constraint condition (20) into (24) The algorithm was able to find some results, but the found points were not the true Maxwell points of the system. The system (23) has too many decoupled states and so there are infinitely many solutions to the minimisation problem (22) that satisfy the constraintstwo trajectories have different initial vectors h(0) and converge to the same point, but the trajectories themselves are the same. An example of such trajectories is shown in Figure 5.

Rolling Disc-Approximation Based on Taylor Linearisation
The second approximation leads us to a different set of equations (based on (11)) of the form (15) Deriving a simplification to the calculation of the optimality functional as seen in the previous section is not trivial in the case of the system (25). The final value is dependent on h 4 , h 5 , which are not decoupled, and is also dependent on other states. Therefore, in this case there is no choice but to compute the final value of the optimality functional (9) at each iteration, which increases the computational complexity.
This time, after running the optimisation, we obtain different trajectories leading to the same points (example trajectory shown in Figure 6). In Figure 7 we can clearly see that the found points span a subspace of the original configuration space. Some dimensions in the data are numerically insignificant, forming planes or lines in specific views, see Figure 8.
To further analyse the form of the subset formed by these points, PCA analysis was performed on the data to discover the principal components in the data along with their specific latent variables which specify the amount of variance each specific component encapsulates. We can then calculated the cumulative sum of this variance for each specific number to see how many of them are needed to describe our data. This is shown in Figure 9. The calculated principal componentsQ = [q 1 T . . .q 5 T ] form an orthonormal basis which can be further reduced using our latent values and cumulative sum. In our case, the first 3 components describe 99.9 % of the variance in our data. So we take only the first 3 columns of the matrixQ which forms a basis for our set of Maxwell points.
Note that this basis is valid only for the range of the space where the data points are available. Hypothesis based on extrapolating the data is possible but needs to be verified further.

Conclusions
In this paper, a numerical algorithm for finding and validating the Maxwell points of an affine control system was proposed in Section 5. The algorithm was applied to two different affine control systems, which were derived by different approximations of dynamics of the vertical rolling disc, as described in Section 6. Both of the approximations serve two purposes, that is, as a case study for the proposed numerical approach and to evaluate the effect of different approximations on the whole process of optimal control design. These numerical algorithms are the result of continuing research based on preliminary results published in [14], which proved that the Maxwell points found using numerical optimisation comply with analytical solutions.
As can be seen from the results, there are advantages and disadvantages of both derived systems. In the nilpotent-based approximation the numerical algorithm collapses to a trivial solution, as the state space trajectories are not dependent on some of the expanded states, so the optimality of the control system is lost in points different from the Maxwells ones. The problem of finding these points is much more complicated and will be part of our next research.
In opposition to this result, the common linearisation based on Taylor approximation, while still being nilpotent, yields quite different results. In this case, the numerical algorithm is able to locate a set of Maxwell points suitable for further system analyses and optimal control design. Unfortunately, the dimension of the control Lie algebra is higher than the dimension of the configuration space.
The proposed algorithms allow us, in theory, to find or validate assumed Maxwell points on any affine control system. This applies mainly in situations when analytical solutions are not available, with the only limitations being the computation power and numerical accuracy of the optimisation task. This paper only demonstrates its function on relatively simple examples, but all of the algorithm's steps should extend well to more complicated cases, as it does not impose any insurmountable conditions on the systems being analysed. The optimisation algorithm can reach arbitrary accuracy (within reason), limited mainly by the numerical accuracy of 64-bit double precision data type. However, increasing the accuracy slows down the computation, and obtaining even a reasonable number of points can take a very long time. The task is by definition suitable for parallelisation on multiple cores of the system to reduce the computation time if needed.
Inaccuracies can also be introduced into the algorithm via integration steps of the numerical ODE solver. This is especially so if the system is stiff. For linear systems we can overcome this issue using an analytical solution to the system of ODEs, but with nonlinear systems only choosing a specialised solver or reducing the step size can solve this problem, which again trades off the computation time.
As seen in Section 6.1 there is a risk of false positives-either the system is ill defined and the constraints specified in (22) are insufficient, or a local minimum occurs when converging towards the Maxwell point. Both of these instances result in false positives and can be eliminated by checking that the point we obtained is truly a Maxwell point of the system. In the future research this problem will be addressed more concretely.
While the numerical algorithm presented in this paper is not perfect and requires understanding and interaction with the user, the results are promising and form a basis for future research. Funding: The work presented in this article has also been supported by the Czech Republic Ministry of Defence (University of Defence development program "Advanced Automated Command and Control System II -PASVR II") as well as by the Faculty of Mechanical Engineering, Brno University of Technology under the projects FSI-S-20-6407 "Research and development of methods for simulation, modelling a machine learning in mechatronics" and FSI-S-20-6187 "Modern methods in applied mathematics". Very special thanks to Associate Professor Jaroslav Hrdina for his understanding and supports for this research.