# Multi-Feeder Restoration using Multi-Microgrid Formation and Management

Valliappan Muthukaruppan, Rongxing Hu, Ashwin Shirsat,
Mesut Baran, Ning Lu, Wenyuan Tang, and David Lubkeman  The authors are with
Department of Electrical and Computer Engineering, North Carolina State
University, Raleigh, NC 27695 USA. This material is based upon work supported
by U.S. Depart- ment of Energy’s Office of Energy Efficiency and Renewable
Energy (EERE) under Solar Energy Technologies Office Award Number DE-
EE0008770.

###### Abstract

This papers highlights the benefit of coordinating resources on mulitple
active distribution feeders during severe long duration outages through multi-
microgrid formation. A graph-theory based multi-microgrid formation algorithm
is developed which is agnostic of the underlying energy management scheme of
the microgrids and solved in a rolling horizon fashion. The algorithm is then
enhanced to handle multiple feeders where formation of long laterals needs to
be avoided due to potential voltage control issues in distribution systems.
The algorithm is evaluated on a synthetic two feeder system derived from
interconnecting two IEEE 123 node system. The results indicate increased
service to loads in the system and better utilization of renewable resources.

###### Index Terms:

multi-feeder restoration, multi-microgrid formation, microgrid management,
active distribution system

##  I Introduction

Forming a microgrid within a distribution system following a significant
weather-related outage is gaining prominence as a practical approach to
enhance resilience at the distribution level. This approach becomes
particularly appealing when the distribution system incorporates substantial
amount of distributed photovoltaic (PV) generation [1]. Nevertheless, the
management of these microgrids at the feeder level or involving multiple
feeders presents several challenges. These challenges encompass constrained
resources for rapid deployment on the feeders and limited capabilities for
real-time monitoring and control within the distribution system [2].

In the presence of significant amount of renewable resources like rooftop
solar and large service territory, it is challenging for a microgrid with low
inertial resources to maintain power balance and support service over a long
period of time [3]. Hence, it is beneficial to manage multiple microgrids
together. When it comes to multi-microgrid management during power system
restoration two methods have found prominence in literature, networked
microgrids [4] and dynamic multi-microgrid formation [5].

There are significant challenges with managing networked microgrids due to
complicated control requirements [6]. Hence, dynamic multi-microgrid is a good
solution especially for Utilities managing microgrids owned by different
stakeholders during extreme weather events. The control architecture is simple
and the boundary of the microgrids can be easily controlled by Utilities. In
existing literature, this problem is solved as a critical service restoration
problem where critical loads are routed to existing microgrid resources
through numerous circuit switches [7]. The issue with this type of restoration
approach is that the topology of microgrid remains same throughout the
restoration process once the critical loads are energized, which may not be
optimal for all operating conditions.

Another approach with multi-microgrid formation is to use mobile devices and
dispatch them from one location to another during the restoration as the
system condition changes [8]. The focus here is the routing of mobile devices
to designated load zones which is solved as a planning problem. Even though
this seems like a viable solution, the time taken to interconnect mobile
devices with distribution system is ignored which can hinder timely
resotration of loads and increase their downtime as well.

Current works in literature do not consider the challenging and realistic
conditions of a distribution system such as limited microgrid resources,
limited controllability, and significantly high solar penetration in the
problem formulation. Furthermore, none of the existing work have considered a
multi-feeder problem and the inherent issues in the setup.

This paper aims at developing a comprehensive restoration strategy for
multiple feeders using multi-microgrid formation. The main contribution of the
proposed method are listed below:

  * •

Multi-feeder restoration through formation of multiple-microgrids has a risk
of leading to long laterals which can lead to voltage issues. We consider a
special formulation of multi-microgrid problem that addresses this issue by
providing DSO with a control on length of laterals that can be formed during
the microgrid formation. Furthermore, the proposed formulation is agnostic of
underlying energy management scheme of the microgrids and significantly
improves the restoration process by coordinating the resources in different
feeders.

  * •

Realistic distribution sytem operating conditions are considered: limited load
and PV visibility and controllability, and limited controllable swtiches on
the feeder. The load and PV are assumed to be controlled only in zones formed
by the existing circuit switches in the network.

The rest of the paper is organized as follows: sec-II introduces the multi-
feeder multi-microgrid formation problem and the interface with the microgrid
energy management schemes, sec-III illustrates the performance of proposed
scheme with a case study using synthetic multi-feeder sytem derived from two
IEEE 123 node feeders and multi-day operating conditions based on field data.

##  II Problem Formulation

###  II-A Microgrid Formation

The underlying multi-microgrid formation problem is a form of graph splitting
problem which can be solved using the single commodity flow method from graph
theory [9]. Let $G=(V,E)$ be a connected and undirected graph with set of
vertices and edges denoted by $V$ and $E$ respectively representing the
distribution network under outage. Let the binary decision variables $y_{ij}$
represent the status of circuit switches connecting load groups $i$,$j$. If
the switch is open, $y_{ij}=0$; otherwise, $y_{ij}=1$. We assume that there
will be only one grid forming source per microgrid, thus the number of
microgrids is equal to the number of master control units. Which means,
$\mid\Pi\mid$ number of microgrids will be formed where $\Pi$ is the set of
grid forming resources in the network.

Objective function is given in (1) where the first term minimizes the total
load shedding in the network during microgrid formation where $D_{j}^{max}$
indicates the total connected load in load zone $j$ and $\mathcal{N}^{d}$
indicates the total number of load zones in the network. The second term
minimizes the fictitious flow $F_{ij}$ in the lines with higher priority to
lines connected directly to critical loads. The larger weight on incoming
flows to critical loads ensures that they are closer to the grid forming
sources of the microgrid thereby ensuring higher probability of service during
the energy management phase.

$\min_{y}\sum_{j=1}^{\mid\mathcal{N}^{d}\mid}(D_{j}^{max}-D_{j})+\sum_{ij\in
E}w_{j}\mid F_{ij}\mid$ (1)

Let, $\mathcal{N}_{p}$ be the index set of number of microgrids to be formed.
The binary variable $x_{i,k}=1$ indicates that load group $i$ belongs to
microrgrid $k$. Equation (2) ensures that each load group belongs to only one
microgrid. While equations (3-4) ensures that nodes $(f,k)$ connected to a
closed switch $l$ belong to the same microgrid $k$.

$\displaystyle\sum_{i\in V}x_{i,k}$ $\displaystyle=1$
$\displaystyle\quad\forall k\in 1,\ldots,\mathcal{N}_{p}$ (2) $\displaystyle
y_{l}$ $\displaystyle=\sum_{k\in 1,\ldots,\mathcal{N}_{p}}z_{l,k}\quad$
$\displaystyle\forall l\in E$ (3) $\displaystyle z_{l,k}$
$\displaystyle=x_{f,k}x_{t,k}\quad$ $\displaystyle\forall l\in E,(f,t)\in
E(l)$ (4)

For all $l\in E$ and $(f,t)\in E(l)$, equation (5) is the McCormick
Linearization of (4).

$\displaystyle z_{l,k}$ $\displaystyle\leq x_{f,k};$ $\displaystyle z_{l,k}$
$\displaystyle\leq x_{t,k};$ (5) $\displaystyle z_{l,k}$ $\displaystyle\leq
x_{f,k}+x_{l,k}-1$

For all $(f(l),t(l))\in E(l)$, $l\in E$, and $i\in V$, (6)-(10) define the DC
power flow equations with limits on the line flow, power generation, and
demand. Where, $T_{l}$ is the net load flowing through line $l$ connecting
load zones $f(l)\text{ and }t(l)$. $P$ and $D$ indicate the total PV
generation and load in the individual load zones. $y_{l}$ indicates the status
of switch $l$.

$\displaystyle\sum_{l:f(l)\rightarrow t(l)}T_{l}-\sum_{l:t(l)\rightarrow
f(l)}T_{l}$ $\displaystyle=P_{f(l)}-D_{f(l)}$ (6) $\displaystyle-
T_{l}^{min}y_{l}$ $\displaystyle\leq T_{l}\leq T_{l}^{max}y_{l}$ (7)
$\displaystyle-M(1-y_{l})$ $\displaystyle\leq T_{l}\leq M(1-y_{l})$ (8)
$\displaystyle P_{i}^{min}$ $\displaystyle\leq P_{i}\leq P_{i}^{max}$ (9)
$\displaystyle D_{i}^{min}$ $\displaystyle\leq D_{i}\leq D_{i}^{max}$ (10)

Equation (11) ensures radiality in the network by determining the total number
of closed switches possible. Here, $\mid V\mid$ is the number of nodes,
$\mid\Pi\mid$ is the number of microgrid resources available in the system,
and $\mid\mathcal{R}\mid$ is the number of load islands. Load islands are load
groups that cannot be connected to any microgrid resources due to existing
faults in the system.

$\sum_{ij\in E}y_{ij}=\mid V\mid-\mid\Pi\mid-\mid\mathcal{R}\mid$ (11)

To ensure the individual microgrid connectivity via mathematical programming
formulations, the single commodity flow method is employed as shown in
equations (12) - (17). The equations are similar to the DC power flow
equations but all nodes are assumed to inject a load of value 1 and the
sources are indicated by $W$

$\displaystyle\sum_{s\in\delta(j)}F_{js}-\sum_{i\in\pi(j)}F_{ij}$
$\displaystyle=1$ $\displaystyle\quad\forall j\in V/\Pi$ (12)
$\displaystyle\sum_{s\in\delta(j)}F_{js}-\sum_{i\in\pi(j)}F_{ij}$
$\displaystyle=W_{j}$ $\displaystyle\quad\forall j\in\Pi$ (13) $\displaystyle-
My_{ij}\leq F_{ij}\leq$ $\displaystyle My_{ij}$ $\displaystyle\quad\forall
ij\in E$ (14) $\displaystyle-M(2-y_{ij})\leq F_{ij}\leq$ $\displaystyle
M(2-y_{ij})$ $\displaystyle\quad\forall ij\in E$ (15) $\displaystyle
W_{j}\leq$ $\displaystyle 1$ $\displaystyle\quad\forall j\in\Pi$ (16)
$\displaystyle F_{ij}\geq$ $\displaystyle n^{min}$ $\displaystyle\quad\forall
i\in\mathcal{N}_{GFM}$ (17)

Consider the simplified graphical representation of a simple multi-feeder
restoration example as shown in fig. 1. Nodes 1 through 5 belong to feeder-1
and 6 through 10 belong to feeder-2. There are two normally open switches
interconnecting the nearby feeders. When restoring multiple feeders using
multiple microgrids, the length of the radial network inside each microgrid
needs to be controlled since the underlying distribution system voltage
regulation may not be capable of handling the different microgrid topologies.
Hence, it is critical to limit the boundary of the microgrids.

Figure 1: Graphical representation of a two feeder restoration problem

To avoid creating long radial networks, we modify the single commodity flow
method such that only leaf nodes (feeder-1: nodes 2 and 5; feeder-2: nodes 6
and 10) will be exchanged between the microgrids during restoration. Leaf
nodes are the nodes connected to the normally open switches in the network.
This strategy is implemented using constraint (17) where we control the
fictious flow on the lines directly connected to the grid forming sources
(root of the graph) in each microgrid.

Figure 2: Fixed and Flexible Topology Classification to avoid voltage issues
caused by long radial networks because of picking up too many load groups from
adjacent feeder.

For example, in figure 2 the left lateral of node-1 with grid forming
resources (GFM) in MG-1 is restricted to two such that load groups 3 and 4 are
always connected to MG-1 while node-5 can be interchanged with the neighboring
microgrids. Since node-2 is directly connected to the GFM in MG-1 the fictious
flow on the right branch of GFM is set to zero. This approach can be extended
any number of feeders with interconnecting switches and different topologies.

###  II-B Energy Management

The advantage of the proposed microgrid formation problem is that it is
agnostic of the underlying energy management scheme (EMS), each microgrid can
use its own EMS implemented through its microgrid controller. No modification
to the underlying control architecture is necessary. In this paper we use our
previously proposed rolling horizon based two-stage hierarchical EMS for each
estabilished microgrid [10].

The decision variables from multi-microgrid formation problem is the optimal
topology of individual microgrids which will be relayed to the microgrid
controllers. The microgrid controllers will then manage the resources and load
groups in their assigned topology until next topology change by multi-
microgrid formation module.

The horizon, time step, and coordination of the various modules in the multi-
microgrid formation and management problem is shown in figure 3. The multi-
microgrid formation module is also formulated as a rolling horizon problem to
avoid frequent topology changes in the microgrids. Note that equations in
section-II-A are defined for a single time step for ease of explanation but
the problem is solved for a whole horizon with multiple timesteps.

Figure 3: Timeline of Different problems and their coordination

Eventhough the multi-microgrid fomration module determines the optimal
topology of each microgrid, the decision to supply the load groups within each
microgrid is still determined by the individual microgrid EMS system depending
on the availability of the resources.

##  III Results

The test system considered for evaluating the proposed algorithm is shown in
figure 4. We have used two IEEE 123 node system which are assumed to be
supplied from same substation to develop the multi-feeder test case. The
location of microgrid resources on feeder-1 is at the substation while in
feeder-2 they are towards the end of the feeder. These locations are pre-
determined and fixed throughout the restoration. Two interconnecting switches
which are normally open at nodes 250 and 350 interconnect the two IEEE 123
feeders. All load groups defined by the shaded polygons have significant
penetration of distributed behind-the-meter PV. There are 3 critical loads per
feeder highlighted as purple stars in the figure. The graphical representation
of this system is shown in figure 1. The rating of the resources are
highlighted in table I. The horizon and time step of the EMS and MMG formation
module are highlighted in table II.

Figure 4: Multi-feeder test system developed from two IEEE 123 node feeders TABLE I: Rating of resources in the multi-feeder test case Resources | Feeder-1 | Feeder-2
---|---|---
Mobile Energy Storage | 3 MW/ 12 MWh | 2 MW/ 8 MWh
Diesel Generator | 4 MW | 4 MW
Total PV | 4 MW | 4 MW
Peak Load (Day-1) | 3.5 MW | 3 MW
Peak Load (Day-2) | 3 MW | 2 MW
TABLE II: Horizon and Time step of different problems in MMG restoration Problem | Horzion | Time Step
---|---|---
Multi-Microgrid Formation | 24 hour | 3 hour
Scheduling Stage-1 | 24 hour | 30 minutes
Dispatching Stage-2 | 30 minutes | 5 minute

Total load and PV on each feeder is highlighted in figure 5. Since, these
feeders are adjacent to each other, the overall PV profile looks similar, but
the load profiles are different to highlight the difference in characteristics
of these feeders. The overall load is considerable lower on day-2 in feeder-2.
With 4 MW of PV in each feeder, it would be challenging to absorb all the PV
in feeder-2 with the smaller 2 MW battery.

Figure 5: Total Load and PV on feeder-1 and feeder-2

Table III highlights the topology changes on day-1 and day-2 of the
restoration. The change in topology on day-1 is minimal since the load and PV
are quite similar between the two microgrids. Initially, more load groups are
connected to MG1 since it has higher sized battery. During PV, load group 10
which is a critical load group is connected to MG1 in exchange for load group
5 which is a non-critical load gorup.

On day-2 we see a lot of topology changes that happen towards the end of
restoration which is due to the limited resources towards the end of
restoration. Also, more load groups are shifted from MG1 to MG2 eventhough it
has a smaller battery size because the peak load on day-2 in MG2 is much lower
than MG1, so MG1 load groups have more chance of being served when connected
to MG2.

TABLE III: Optimal Topology Changes on Day-1 and Day-2 Duration | MG1 Load Groups | MG2 Load Groups | Comments
---|---|---|---
Day-1
00:00 to 12:00 | 1, 2, 3, 4, 5, 6 | 7, 8, 9, 10 | LG-6 connected to MG1 due to bigger battery size.
12:00 to 15:00 | 1, 2, 3, 4, 10 | 5, 6, 7, 8, 9 |  | Since 2MW battery in MG2 cannot absorb all PV,
---
might lead to outage of LG-10.
15:00 to 24:00 | 1, 2, 3, 4, 5, 6 | 7, 8, 9, 10 | LG-6 connected back to MG1 due to bigger battery size.
Day-2
00:00 to 12:00 | 1, 2, 3, 4, 5, 6 | 7, 8, 9, 10 |  | More chance of load being served in MG1 since it has
---
higher battery size. Hence, LG6 connected to MG1.
12:00 to 15:00 | 1, 2, 3, 4, 5 | 6, 7, 8, 9, 10 | Default topology of feeder is maintained.
15:00 to 18:00 | 1, 2, 3, 4, 5, 6 | 7, 8, 9, 10 | More chance of service to load groups in MG1.
18:00 to 21:00 | 1, 2, 3, 4 | 5, 6, 7, 8, 9, 10 |  More loads connected to MG2 even though lower battery size because loading is much lower in MG2 compared to MG1 on day-2.
21:00 to 24:00 | 1, 3, 4 |  | 2, 5, 6, 7, 8,
---
9, 10

The SoC of mobile energy storage devices and the fuel of diesel generators in
both microgrids are highlighted in figure 6. Both battery and DG have been
appropriately used and completely depleted by the end of the restoration. The
devices were kept within limits even with significant number of topology
changes.

Figure 6: (top) State of Charge of Mobile Energy Storage Devices and (bottom)
Fuel usage of Diesel Generators during the two day restoration period.

The load group connectivity status is highlighted in figure 7 which shows how
different load groups are served by the two microgrids in the system. There
are instances where load groups are not served by the microgrids even though
they are assigned to the microgrid by multi-microgrid algorithm. The downside
of multi-microgrid algorithm is that when a load group is assigned to a
particular microgrid there is no guarantee that the load group will be served
by that microgrid which is subject to availability of resources and handled by
the individual energy management scheme. Nevertheless, we do minimize such
outage of critical loads which can be communicated to the customer in advance
and the customer can prepare back up local generation during these scheduled
outage hours.

Figure 7: Load Group connectivity status highlighted using different colors to
indicate which microgrid was utilized to energize the load groups. Two
topology important topology changes to improve critical load connectivity is
highlighted in red with corresponding topologies.

Critical load groups 4 and 10 have significant duration of unserved periods by
both microgrids and this is because under any topology configuration the
distance from load group 4 and 10 to microgrid resources 1 and 7 is long.
Hence, to serve these load groups all the other load groups between the
resources and these groups need to be served which is challenging with limited
resources and during peak load duration. Hence, load groups 4 and 10 have
maximum number of unserved duration.

To analyze the performance of our proposed multi-microgrid algorithm we
compare the results against base case where the topology of the microgrids are
fixed throught the restoration period. Nodes 1 through 5 are always connected
to MG1 and 6 through 10 to MG2 through out the 2 day restoration.

Percentage service to each load group is compared between proposed scheme and
base case in figure 8. It can be seen that service to critical loads is
significantly increased with proposed scheme and service to leaf nodes like
load groups 5 and 6 is also considerably increased. Due to efficient
management of the leaf nodes in flexible topology, other nodes pertaining to
fixed topology like load group 4 and 8 also have improved service. This
highlights the advantage of coordinating multi-microgrids during restoration
rather than using a fixed topology. The overall PV utilization with proposed
scheme is 82% which is higher than the utilization by base case with 79.12%.

Figure 8: Percentage of Load Groups served during restoration period between
base case (fixed topology) and proposed scheme (flexible topology). Critical
Load Groups: 2, 3, 4, 7, 9, 10.

##  IV Conclusion

This paper highlights the advantages of optimizing resource coordination
across multiple feeders during prolonged outages through dynamic multi-
microgrid formation. We introduce a graph-theory based multi-microgrid
formation problem, extending it to multiple feeders with adjustable lateral
lengths to mitigate potential voltage issues. A case study, using a synthetic
two-feeder system with varied load characteristics derived from a realistic
IEEE-123 node system, demostrates substantial enhancements in service to
critical and non-critical loads. Coordinating microgrid resources across
multiple feeders also results in improved utilization of photovoltaic (PV)
systems compared to a fixed microgrid topology.

## References

  * [1] D. Labonte, “Microgrids Can Provide Resiliency During Extreme Weather Events.” [Online]. Available: https://www.powermag.com/microgrids-can-provide-resiliency-during-extreme-weather-events/
  * [2] M. H. Saeed, W. Fangzong, B. A. Kalwar, and S. Iqbal, “A Review on Microgrids’ Challenges & Perspectives,” vol. 9, pp. 166 502–166 517.
  * [3] A. Pratt, “Addressing Challenges for Single Microgrids and Networked Microgrids at Large Scales.” [Online]. Available: https://www.osti.gov/biblio/1778193
  * [4] A. Arif and Z. Wang, “Service restoration in resilient power distribution systems with networked microgrid,” in _2016 IEEE Power and Energy Society General Meeting (PESGM)_ , pp. 1–5.
  * [5] M. A. Igder and X. Liang, “Service Restoration Using Deep Reinforcement Learning and Dynamic Microgrid Formation in Distribution Networks,” vol. 59, no. 5, pp. 5453–5472.
  * [6] M. Islam, F. Yang, and M. Amin, “Control and optimisation of networked microgrids: A review,” vol. 15, no. 6, pp. 1133–1148.
  * [7] Y. Xu, C.-C. Liu, K. P. Schneider, F. K. Tuffner, and D. T. Ton, “Microgrids for Service Restoration to Critical Load in a Resilient Distribution System,” vol. 9, no. 1, pp. 426–437.
  * [8] S. Lei, J. Wang, C. Chen, and Y. Hou, “Mobile Emergency Generator Pre-Positioning and Real-Time Allocation for Resilient Response to Natural Disasters,” vol. 9, no. 3, pp. 2030–2041.
  * [9] T. Ding, Y. Lin, Z. Bie, and C. Chen, “A resilient microgrid formation strategy for load restoration considering master-slave distributed generators and topology reconfiguration,” vol. 199, pp. 205–216.
  * [10] V. Muthukaruppan, A. Shirsat, R. Hu, V. Paduani, B. Xu, Y. Li, M. Baran, N. Lu, D. Lubkeman, and W. Tang. Feeder Microgrid Management on an Active Distribution System during a Severe Outage. arXiv.org. [Online]. Available: https://arxiv.org/abs/2208.10712v1


# The energy equation for nonholonomic systems with nonlinear kinematic
restrictions: a special category of constraints

F. Talamucci  DIMAI, Dipartimento di Matematica e Informatica “Ulisse Dini”,
Università degli Studi di Firenze, Italy
e-mail<EMAIL_ADDRESS>
###### Abstract

The main topic of this work concerns the formulation of the equations of
motion and the consequent energy balance that they imply for this type of
systems, In particular, the analytical development that we will carry out on
the equations of motion has as its objective the energy balance of the system.
the delicate question of defining the displacements admitted by the system
leads, as we shall see, to a non-univocal definition of the energy of the
system, which finds coherence and unity for a particular class of nonholonomic
constraints.

2010 Mathematics Subject Classification: 37J60, 70F25, 70H03.

37J60: Nonholonomic dynamical systems 70F25: Nonholonomic systems 70H03:
Lagrange’s equations 70G10 Generalized coordinates; event, impulse-energy,
configuration, state, or phase space for problems in mechanics

Keywords:

Nonholonomic mechanical systems - Linear and nonlinear kinematic constraints
-Lagrangian equations of motion - Voronec’s equations of motion.

##  1 Introduction

The study of mechanical systems subject to kinematic constraints has
undoubtedly registered a growing interest, starting from the historical works
that laid the foundations of the theory, among which we mention [1], [5],
[14]. A text of fundamental importance which takes stock of the situation up
to that moment and which does the groundwork for many future developments is
undoubtedly [21].

We can state that the treatment of non-holonomic systems is still today at an
incomplete stage in which there is no general theory, comparable to that of
holonomic systems, nor a satisfactory understanding of fundamental questions
such as the link between symmetries and conservation laws or the extension of
Noether’s theorem.

Truthfully, the category of nonholonomic constraints is wide: the definition
and characterization of holonomic and nonholonomic constraints appears in
numerous texts on mechanics and kinematics, based at times on non-unique
choices and criteria. Moreover, in textbooks the treatment of the mechanics of
nonholonomic systems is often reduced to the study of particular systems,
developed with ad hoc techniques. In this regard it should be noted that since
the first examples and models studied (for example [2], [31]), typically rigid
bodies that slip or roll without sliding on surfaces, find the ideal treatment
in the use of the so-called pseudo-coordinates, selected in a targeted way for
the problem ([13] [3]), starting from the theoretical foundations exposed in
[15].

It must be said that there is a remarkable conceptual and methodological
difference between the case of linear kinematic constraints and that of non-
linear constraints: if on one hand the first type of constraints can be
treated convincingly as a plain extension of the holonomic case, on the other
the case of nonlinear conditions on the kinetic variables does not seem to
have a definitive theoretical arrangement.

In both cases of linear and nonlinear nonholonomic problems we can broadly
outline more than one approach to nonholonomic systems, possibly proposing a
mixed way:

  * •

analytical methods, following and extending the Lagrangian formalism which is
familiar and consolidated in the case of geometric constraints (among others
[32], [33], [4], [18]),

  * •

an approach via the formalism of Lagrangian multipliers, where the derivatives
with respect to the kinetic variables of the constraint equations replace
those with respect to the geometric coordinates of the holonomic case ([34],
[24], [22], [19]),

  * •

methods based on a geometrical support, extending the holonomic formulation
via differentiable manifolds to nonholonomic systems, by involving more
complex aspects of the theory,as the jet manifold approach ([20], [17], [26]),

  * •

a variational approach, essentially investigating the possibility of modifying
the principle in order to incorporate the case of nonholonomic constraints
([23], [25]); the analysis developed in [10] and in [12] shows that the
problem is delicate and still open.

The mathematical investigation presented just after substantially favours the
first approach.

In the context of what we have introduced, the present work develops through
the following selection of the modelling areas and the following steps. We
will consider nonholonomic constraints which concern (also in a non-linear
way) restrictions on the kinematic state (both geometric coordinates and
kinetic variables), possibly depending on time explicitly (rheonomic
constraints). To the same extent that holonomic restrictions allow to identify
the essential positional parameters, we use the kinematic type restrictions to
give rise to a subset of independent generalized velocities, in order to write
the equations of motion and develop the theory according to the precepts of
the Lagrangian formalism. This setting based on the selection among kinematic
parameters is inspired by the dated but fundalmental work [30], where the
theory is for linear kinematic restrictions (and still present in the linear
version in [21]), here extended to the general case of nonlinear kinematic
constraints.

As it is known, for holonomic systems an energy-type equation can be deduced
from the equations of motion and the concept of energy of the system can be
opportunely introduced. If we implement the same formal procedure to the
nonoholonomic case we will detect some differences if the starting point is
the set of equations in the essential parameters or the set of equations with
the lagrangian multipliers (where the independent velocities do not need to be
selected). The discrepancy disappears for a remarkable class of kinematic
constraints, for which even the notion of virtual displacements theorized to
infer equations of motion (${\rm{\check{C}}}$hetaev condition, [7], [8], [9])
turns out to be natural and in step with the holonomic case.

###  1.1 The equations of motion

The model layout consists in a mechanical system whose configurations are
uniquely determined by the set of Lagrangian coordinates $q_{1},\dots,q_{n}$.
Besides that, we assume that the system is subject to the following
restrictions, involving the generalized velocities and possible the tiem $t$:

$\begin{cases}\phi_{1}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{n},t)=0\\\
\qquad\qquad\dots\dots\dots\\\ \qquad\qquad\dots\dots\dots\\\
\phi_{k}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{n},t)=0\end{cases}$
(1)

In order to write the equations of motion A direct approach consists in
writing the Lagrangian equations of motion

$\dfrac{d}{dt}\dfrac{\partial T}{\partial{\dot{q}}_{i}}-\dfrac{\partial
T}{\partial q_{i}}={\cal F}^{(i)}+{\cal R}^{(i)},\qquad i=1,\dots,n$ (2)

where $T(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{n},t)$ is the
kinetic energy of the system and ${\cal F}^{(q_{i})}$, ${\cal R}^{(q_{i})}$
are the $i$–th lagrangian components of the acting forces and of the
constraint forces, respectively. The structure of the constraints can be
outlined by means of the Lagrangian multipliers by setting

${\cal
R}^{(i)}=\sum\limits_{\nu=1}^{k}\lambda_{\nu}\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{i}},\qquad
i=1,\dots,n$ (3)

so that the systems of $n+k$ equations (2), (1) contain properly the $n+k$
unknown functions $q_{i}$, $i=1,\dots,n$, $\lambda_{\nu}$, $nu=1,\dots,k$.

Although the generalized velocities ${\dot{q}}_{i}$, $i=1,\dots,n$ are not
independent, adopting the point of view (2) allows to not explicitly declare
the virtual displacements and simply considering the constraints of the type
(3) as smooth.

A quite different formal path consists in explicitly declaring the virtual
displacements and directly ignoring the presence of the constraint forces: the
first step is to use the restrictions (1) in order to construct the set of
virtual displacements.

The set of conditions (1) can be made explicit if the following condition
holds (without losing generality, except for re-enumerating the variables):

$det\,\left(\dfrac{\partial\phi_{i}}{\partial{\dot{q}}_{m+j}}\right)_{i,j=1,\dots,k}\not=0,\quad
m=n-k$ (4)

In this case it is possible to determine the functions

$\begin{cases}{\dot{q}}_{m+1}=\alpha_{1}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)\\\
\dots\\\ \dots\\\
{\dot{q}}_{n}=\alpha_{k}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)\end{cases}$
(5)

with $m=n-k$. The kinetic parameters ${\dot{q}}_{r}$, $=1,\dots,m$, play the
role of independent velocities that is, fixing a $m$–uple
${\dot{q}}_{1},\dots,{\dot{q}}_{m}$ establishes the entire kinematic state of
the system in the position $q_{1}$, $\dots$, $q_{n}$.

In the presence of the setting (5), the formulation of the equations of motion
can be based on the assumption that the virtual displacements of the system
are

$\delta
q_{m+\nu}=\sum\limits_{r=1}^{m}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}\delta
q_{r},\qquad\nu=1,\dots,k$ (6)

where the virtual variations $\delta q_{1}$, $\dots$, $\delta q_{r}$ are
independent. By assuming (8), the corresponding equations are

$\dfrac{d}{dt}\dfrac{\partial T}{\partial{\dot{q}}_{r}}-\dfrac{\partial
T}{\partial
q_{r}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}\left(\dfrac{d}{dt}\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}-\dfrac{\partial T}{\partial
q_{m+\nu}}\right)={\cal
F}^{(r)}+\sum\limits_{\nu=1}^{k}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}{\cal
F}^{(m+\nu)},\qquad r=1,\dots,m.$ (7)

###### Remark 1

If ${\bf r}_{1}$, $\dots$, ${\bf r}_{N}$ are the vectors which locate the
points of the system, assumption (6) corresponds to

$\delta{\bf r}_{j}=\sum\limits_{r=1}^{m}\dfrac{\partial{\dot{\bf
r}}_{j}}{\partial{\dot{q}}_{r}}\delta q_{r},\qquad j=1,\dots,N$ (8)

Actually, since

${\dot{\bf r}_{j}}=\sum\limits_{r=1}^{m}\dfrac{\partial{\bf r}_{j}}{\partial
q_{r}}{\dot{q}}_{r}+\sum\limits_{\nu=1}^{k}\dfrac{\partial{\bf
r}_{j}}{\partial q_{m+\nu}}\alpha_{\nu}+\dfrac{\partial{\bf r}_{j}}{\partial
t},\qquad j=1,\dots,N$ (9)

we see that

$\dfrac{\partial{\dot{\bf r}}_{j}}{\partial{\dot{q}}_{r}}=\dfrac{\partial{\bf
r}_{j}}{\partial q_{r}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial{\bf
r}_{j}}{\partial
q_{m+\nu}}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}},\qquad
j=1,\dots,N,\quad r=1,\dots,m$

hence (8) takes the form

$\delta{\bf r}_{j}=\sum\limits_{r=1}^{m}\left(\dfrac{\partial{\bf
r}_{j}}{\partial
q_{r}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}\dfrac{\partial{\bf
r}_{j}}{\partial q_{m+\nu}}\right)\delta q_{r},\qquad j=1,\dots,N$ (10)

which entails the same information as (6) for what concerns the virtual
displacements. In terms of virtual velocities the assumption (8) can be
formally written as

${\widehat{\dot{\bf r}}_{j}}=\sum\limits_{r=1}^{m}\dfrac{\partial{\dot{\bf
r}}_{j}}{\partial{\dot{q}}_{r}}{\dot{q}}_{r}=\sum\limits_{r=1}^{m}\dfrac{\partial{\bf
r}_{j}}{\partial q_{r}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial{\bf
r}_{j}}{\partial
q_{m+\nu}}\sum\limits_{r=1}^{m}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}{\dot{q}}_{r},\qquad
j=1,\dots,N$ (11)

which does not match the set of velocities consistent with the instantaneous
configuration of the system (i. e. at a blocked time $t$), obtained by
removing the term $\frac{\partial{\bf r}_{j}}{\partial t}$ in (9).

Coming back to the equations (7), we point out that they can be appropriately
written in terms of the only independent generalized velocities
${\dot{q}}_{r}$, $r=1,\dots,m$: indeed, by defining the kinetic energy $T^{*}$
restricted to the constraints (5) as

$\begin{array}[]{l}T^{*}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)\\\
\\\
=T(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},\alpha_{1}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{n},t),\dots,\alpha_{k}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{n},t),t)\end{array}$
(12)

it can be seen that the equations of motion (7) take the form (see [28] for
details)

$\dfrac{d}{dt}\dfrac{\partial T^{*}}{\partial{\dot{q}}_{r}}-\dfrac{\partial
T^{*}}{\partial q_{r}}-\sum\limits_{\nu=1}^{k}\dfrac{\partial T^{*}}{\partial
q_{m+\nu}}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}_{r}}}-\sum\limits_{\nu=1}^{k}B_{r}^{\nu}\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}={\cal
F}^{(q_{r})}+\sum\limits_{\nu=1}^{k}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}{\cal
F}^{(q_{m+\nu})},\qquad r=1,\dots,m$ (13)

where

$B_{r}^{\nu}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)=\dfrac{d}{dt}\left(\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}\right)-\dfrac{\partial\alpha_{\nu}}{\partial
q_{r}}-\sum\limits_{\mu=1}^{k}\dfrac{\partial\alpha_{\mu}}{\partial{\dot{q}}_{r}}\dfrac{\partial\alpha_{\nu}}{\partial
q_{m+\mu}},\quad r=1,\dots,m,\;\;\nu=1,\dots,k$ (14)

Equations (13) contain the unknown functions $q_{1}$, $\dots$, $q_{n}$ but
only the derivatives ${\dot{q}}_{r}$, ${\mathop{q}\limits^{..}}_{r}$
$r=1,\dots,m$ are present, owing to (5); the variables
$({\dot{q}}_{k+1},\dots,{\dot{q}}_{n})$ that are present in $\frac{\partial
T}{\partial{\dot{q}}_{m+\nu}}$ of (13) have to be expressed in terms of
$(q_{1},\dots,q_{n},$ ${\dot{q}}_{1},\dots,{\dot{q}}_{m},t)$ by making use of
(5). The set of equations (13) joined to (5) represents a system of $m+k=n$
equations in the $n$ unknown functions $(q_{1},\dots,q_{n})$. Such a system
represents an extension to the nonlinear case of the so-called Voronec
equations, developed in [21] for the linear constraints

$\left\\{\begin{array}[]{l}\sum\limits_{i=1}^{n}\sigma_{1,i}(q_{1},\dots,q_{n},t){\dot{q}}_{i}+\zeta_{1}(q_{1},\dots,\dots,q_{n},t)=0\\\
\qquad\qquad\dots\dots\dots\\\ \qquad\qquad\dots\dots\dots\\\
\sum\limits_{i=1}^{n}\sigma_{k,i}(q_{1},\dots,q_{n},t){\dot{q}}_{i}+\zeta_{k}(q_{1},\dots,\dots,q_{n},t)=0\end{array}\right.$
(15)

associated to the explicit form (5), whenever the $k\times n$ matrix
$(\sigma)_{\nu,i}$ has full rank $k$, according to (4),

$\begin{cases}{\dot{q}}_{m+1}=\sum\limits_{r=1}^{m}\alpha_{1,r}(q_{1},\dots,q_{n},t){\dot{q}}_{r}+\beta_{1}(q_{1},\dots,q_{n},t)\\\
\dots\\\ \dots\\\
{\dot{q}}_{n}=\sum\limits_{r=1}^{m}\alpha_{k,r}(q_{1},\dots,q_{n},t){\dot{q}}_{r}+\beta_{k}(q_{1},\dots,q_{n},t)\end{cases}$
(16)

for suitable coefficients $\alpha_{nu,r}$ and $\beta_{r}$, $\nu=1,\dots,k$,
$r=1,\dots,m$. In this simpler context, (6) is

$\delta
q_{m+\nu}=\sum\limits_{r=1}^{m}\alpha_{\nu,r}(q_{1},\dots,q_{n},t)\delta
q_{r},\qquad\nu=1,\dots,k.$

###### Remark 2

The case of merely holonomic constraints corresponds to the absence of the
conditions (1) and by deleting all the terms containing the $\alpha_{\nu}$
functions in (13) we find the usual Euler-Lagrange equations of motion with
$T^{*}=T$.

###### Remark 3

Although the set of equations (13) appears to depend on the selection of the
independent velocities, hence the $\alpha_{\nu}$ functions, the particular
dependence on the latter ones only through the derivatives allows us to write
the coefficients in terms of the constraint functions: this occurs using the
relations that hold for any index $\nu=1,\dots,k$:

$\begin{array}[]{ll}\dfrac{\partial\phi_{\nu}}{\partial
q_{i}}+\sum\limits_{\mu=1}^{k}\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{m+\mu}}\dfrac{\partial\alpha_{\mu}}{\partial
q_{i}}=0,&i=1,\dots,n\\\ \\\
\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{r}}+\sum\limits_{\mu=1}^{k}\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{m+\mu}}\dfrac{\partial\alpha_{\mu}}{\partial{\dot{q}}_{r}}=0,&r=1,\dots,m.\end{array}$
(17)

For a fixed $i=1,\dots,n$ and writing the first line for each $j=1,\dots,k$,
we obtain a non–singular linear system of $k$ equations in the $k$ unknowns
$\dfrac{\partial\alpha_{\nu}}{\partial q_{i}}$, $\nu=1,\dots,k$. The same
consideration applies to the second set of conditions, where by fixing
$r=1,\dots,m$ and letting assume for $j$ the values $1,\dots,k$ we get a
non–singular linear system of $k$ equations for the $k$ unknowns
$\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}$, $\nu=1,\dots,k$.

###### Remark 4

A point of contact between the two approaches (2) (together with (3)) and (13)
can be found in the so-called ${\rm{\check{C}}}$hetaev condition ([7], [8],
[9])

$\sum\limits_{i=1}^{n}\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{i}}\delta
q_{i}=0,\qquad i=1,\dots,\nu$ (18)

($\phi_{\nu}$ are the constraints functions (1)). Obviously, (18) implies that
the virtual work of the constraint forces is null, provided that they are
formulated as in (3):

$\sum\limits_{i=1}^{n}{\cal R}^{(i)}\delta
q_{i}=\sum\limits_{\nu=1}^{k}\sum\limits_{i=1}^{n}\lambda_{\nu}\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{i}}\delta
q_{i}=0$

On the other hand, assumption (10) provides a set of virtual displacements to
(3), which makes the constraints ideal: actually

$\sum\limits_{i=1}^{n}{\cal R}^{(i)}\delta
q_{i}=\sum\limits_{r=1}^{m}\left({\cal R}^{(r)}+\sum\limits_{\nu=1}^{k}{\cal
R}^{(m+\nu)}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}\right)\delta
q_{r}=\left(\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{r}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{m+\nu}}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}\right)\delta
q_{r}=0$ (19)

owing to the relations (17), second group.

#### The Lagrangian function

If the generalized forces are such that

${\mathcal{F}}^{(i)}=\dfrac{\partial U}{\partial
q_{i}}(q_{1},\dots,q_{n},t),\quad i=1,\dots,n\quad$ (20)

for a suitable function $U$, then it is possible to define the Lagrangian
function ${\cal L}=T+U$ and (2) (with assumption (3)) and (7) can be written,
respectively, as

$\displaystyle\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{i}}-\dfrac{\partial{\cal L}}{\partial
q_{i}}=\sum\limits_{\nu=1}^{k}\lambda_{\nu}\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{i}}$
$\displaystyle i=1,\dots,n$ (21)
$\displaystyle\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{r}}-\dfrac{\partial{\cal L}}{\partial
q_{r}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}\left(\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}-\dfrac{\partial{\cal L}}{\partial
q_{m+\nu}}\right)=0$ $\displaystyle r=1,\dots,m$

In turn, the restricted kinetic energy (12) makes us define

${\cal
L}^{*}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)=T^{*}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)+U(q_{1},\dots,q_{n},t).$
(22)

In terms of the function (22) the equations of motion (13) become

$\dfrac{d}{dt}\dfrac{\partial{\cal
L}^{*}}{\partial{\dot{q}}_{r}}-\dfrac{\partial{\cal L}^{*}}{\partial
q_{r}}-\sum\limits_{\nu=1}^{k}\dfrac{\partial{\cal L}^{*}}{\partial
q_{m+\nu}}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}_{r}}}-\sum\limits_{\nu=1}^{k}B_{r}^{\nu}\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}=0,\qquad r=1,\dots,m$ (23)

where $B_{r}^{\nu}$ are the same functions (14).

##  2 The energy equation

The next aim is to deduce from (23) an information concerning the energy of
the system: to this end, let us briefly recall the familiar case of holonomic
systems (HC) (which can be identified with $m=n$ and $T^{*}=T$ in (12)): the
lagrangian equations of motion imply the following formula:

$\dfrac{d}{dt}\dfrac{\partial T}{\partial{\dot{q}}_{i}}-\dfrac{\partial
T}{\partial
q_{i}}={\mathcal{F}}^{(i)},\;i=1,\dots,n\;\;\Rightarrow\;\;\dfrac{d}{dt}\left(\sum\limits_{i=1}^{n}{\dot{q}}_{i}\dfrac{\partial
T}{\partial{\dot{q}}_{i}}-T\right)=\sum\limits_{i=1}^{n}{\mathcal{F}}^{(q_{i})}{\dot{q}}_{i}-\dfrac{\partial
T}{\partial t}$ (24)

Furthermore, if (20) holds then
$\sum\limits_{i=1}^{n}{\dot{q}}_{i}{\mathcal{F}}^{(i)}=\dfrac{dU}{dt}-\dfrac{\partial
U}{\partial t}$ and (24) is simply, by defining with ${\cal L}=T+U$,

$\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{i}}-\dfrac{\partial{\cal L}}{\partial
q_{i}}=0\;\;\Rightarrow\;\;\dfrac{d}{dt}\left(\sum\limits_{i=1}^{n}{\dot{q}}_{i}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{i}}-{\cal L}\right)=-\dfrac{\partial{\cal L}}{\partial
t}$ (25)

showing the total variation of the hamiltonian function (in round brackets,
associable to the energy of the system) with respect to the possible explicit
dependence on time of the Lagrangian function; synthetically we recall that
(24) is achieved by multiplying each equation of motion by ${\dot{q}}_{i}$, by
summing the equations and by opportunely rearranging the terms.

We examine now the case of nonlinear nonholonomic systems (NNC). For our
purposes it is reasonable to consider the case of forces (20), since the
presence of forces that cannot be traced back to a potential does not affect
the meaning of the results. The equations of motion to be taken into
consideration are therefore either (21) or (23).

In order to better delineate the definition of energy in nonholonomic systems,
let us start with the setting (21), first set of equations, recurrent in
literature (, ). Retracing what was done before for (HC) systems, that is
multiplying the $i$–th equation by ${\dot{q}}_{i}$ and summing up from $1$ to
$n$, after rearranginng the terms one gets

$\dfrac{d}{dt}\left(\sum\limits_{i=1}^{n}{\dot{q}}_{i}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{i}}-{\cal L}\right)=-\dfrac{\partial{\cal L}}{\partial
t}+\sum\limits_{\nu=1}^{k}\sum\limits_{i=1}^{n}\lambda_{\nu}\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{i}}{\dot{q}}_{i}.$
(26)

With respect to (25), the additional terms of the double sum provide the
contribution of the kinematic constraints, by means of the power
$\sum\limits_{i=1}^{n}{\cal R}^{(i)}{\dot{q}}_{i}$ of the of the constraint
forces (3). We stress that the function ${\cal L}$ without the asterisk
indicates that no replacement (5) has been made) and we do not expect the
power of the constraint forces to be zero.

The main aspect that we are going to examine more deeply is the clarifying
role of the explicit functions $\alpha_{\nu}$, $\nu=1,\dots,k$, in order to
examine the energy balance of the system. Let us start with the following

###### Proposition 1

The constraint forces (3) verify

$\sum\limits_{i=1}^{n}{\cal
R}^{(i)}{\dot{q}}_{i}=\sum\limits_{\nu=1}^{k}(\alpha_{\nu}-{\overline{\alpha}}_{\nu})\left(\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}-\dfrac{\partial{\cal L}}{\partial
q_{m+\nu}}\right)$ (27)

where $\alpha_{\nu}$, $\nu=1,\dots,k$ are the functions (5) and

${\overline{\alpha}}_{\nu}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)=\sum\limits_{r=1}^{m}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}{\dot{q}}_{r}.$
(28)

Proof. The comparison of (21) considered for $i=1,\dots,m$ with the
immediately following equations implies

$\sum\limits_{\nu=1}^{k}\lambda_{\nu}\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{r}}=-\sum\limits_{\nu=1}^{k}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}\left(\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}-\dfrac{\partial{\cal L}}{\partial
q_{m+\nu}}\right),\qquad r=1,\dots,m$

thus

$\sum\limits_{\nu=1}^{k}\sum\limits_{r=1}^{m}\lambda_{\nu}\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{r}}{\dot{q}}_{r}=-\sum\limits_{\nu=1}^{k}\overbrace{\sum\limits_{r=1}^{m}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}{\dot{q}}_{r}}^{={\overline{\alpha}}_{\nu}}\left(\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}-\dfrac{\partial{\cal L}}{\partial
q_{m+\nu}}\right).$

On the other hand, considering $i=m+1,\dots,n$ in (21) one easily gets

$\sum\limits_{\nu,\mu=1}^{k}\lambda_{\mu}\dfrac{\partial\phi_{\mu}}{\partial{\dot{q}}_{m+\nu}}{\dot{q}}_{m+\nu}=\sum\limits_{\nu=1}^{k}\overbrace{{\dot{q}}_{m+\nu}}^{=\alpha_{\nu}}\left(\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}-\dfrac{\partial{\cal L}}{\partial
q_{m+\nu}}\right).$

Putting together the last two formulas we overall obtain

$\sum\limits_{\nu=1}^{k}\sum\limits_{i=1}^{n}\lambda_{\nu}\dfrac{\partial\phi_{\nu}}{\partial{\dot{q}}_{i}}{\dot{q}}_{i}=\sum\limits_{\nu=1}^{k}(-{\overline{\alpha}}_{\nu}+\alpha_{\nu})\left(\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}-\dfrac{\partial{\cal L}}{\partial
q_{m+\nu}}\right)$

that is (27). $\quad\square$

The property (27) allows us to write (26) in the form

$\dfrac{d}{dt}\left(\sum\limits_{i=1}^{n}{\dot{q}}_{i}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{i}}-{\cal L}\right)=-\dfrac{\partial{\cal L}}{\partial
t}-\sum\limits_{\nu=1}^{k}({\overline{\alpha}}_{\nu}-\alpha_{\nu})\left(\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}-\dfrac{\partial{\cal L}}{\partial
q_{m+\nu}}\right)$ (29)

The latter formula shows that the difference with the holonomic case (25) lies
exclusively in the possible difference between $\alpha_{\nu}$ and the function
${\overline{\alpha}}_{\nu}$ defined in (28).

It is natural at this point to apply the same technique to the equations of
motion (23) in order to obtain an information analogous to (26) involving this
time the restricted function ${\cal L}^{*}$: to this end we prove the
following

###### Proposition 2

The equations of motion (23) imply the relation

$\dfrac{d}{dt}\left(\sum\limits_{r=1}^{m}{\dot{q}}_{r}\dfrac{\partial{\cal
L}^{*}}{\partial{\dot{q}}_{r}}-{\cal L}^{*}\right)=-\dfrac{\partial{\cal
L}^{*}}{\partial
t}+\sum\limits_{\nu=1}^{k}({\overline{\alpha}}_{\nu}-\alpha_{\nu})\dfrac{\partial{\cal
L}^{*}}{\partial
q_{m+\nu}}+\sum\limits_{\nu=1}^{k}{\overline{B}}_{\nu}\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}$ (30)

where ${\overline{\alpha}}_{\nu}$ is defined in (28) and

${\overline{B}}_{\nu}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)=\sum\limits_{r=1}^{m}B_{r}^{\nu}{\dot{q}}_{r},\qquad\nu=1,\dots,k$
(31)

with $B_{r}^{\nu}$ tbe same coefficients (14).

Proof. As before, it is sufficient to multiply each equation (23) by
${\dot{q}}_{r}$ and sum up with respect to the index $r$: bearing in mind the
following formula of differentiation

$\frac{d}{dt}{\cal
L}^{*}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)=\sum\limits_{r=1}^{m}\left(\frac{\partial{\cal
L}^{*}}{\partial q_{r}}{\dot{q}}_{r}+\frac{\partial{\cal
L}^{*}}{\partial{\dot{q}}_{r}}{\mathop{q}\limits^{..}}_{r}\right)+\sum\limits_{\nu=1}^{k}\frac{\partial{\cal
L}^{*}}{\partial q_{m+\nu}}\alpha_{\nu}+\frac{\partial{\cal L}^{*}}{\partial
t}$ (32)

we see that

$\sum\limits_{r=1}^{m}{\dot{q}}_{r}\left(\frac{d}{dt}\dfrac{\partial{\cal
L}^{*}}{\partial{\dot{q}}_{r}}-\dfrac{\partial{\cal L}^{*}}{\partial
q_{r}}\right)-\sum\limits_{\nu=1}^{k}\dfrac{\partial{\cal L}^{*}}{\partial
q_{m+\nu}}\sum\limits_{r=1}^{m}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}_{r}}}{\dot{q}}_{r}-\sum\limits_{\nu=1}^{k}\sum\limits_{r=1}^{m}B_{r}^{\nu}{\dot{q}}_{r}\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}$
$=\frac{d}{dt}\left(\sum\limits_{r=1}^{m}{\dot{q}}_{r}\dfrac{\partial{\cal
L}^{*}}{\partial{\dot{q}}_{r}}\right)-\overbrace{\sum\limits_{r=1}^{m}\left(\dfrac{\partial{\cal
L}^{*}}{\partial q_{r}}{\dot{q}}_{r}+\dfrac{\partial{\cal
L}^{*}}{\partial{\dot{q}}_{r}}{\mathop{q}\limits^{..}}_{r}\right)}^{=\frac{d{\cal
L}^{*}}{dt}-\frac{\partial{\cal L}^{*}}{\partial
t}-\sum\limits_{\nu=1}^{k}\frac{\partial{\cal L}^{*}}{\partial
q_{m+\nu}}\alpha_{\nu}}-\sum\limits_{\nu=1}^{k}\dfrac{\partial{\cal
L}^{*}}{\partial
q_{m+\nu}}{\overline{\alpha}}_{\nu}-\sum\limits_{\nu=1}^{k}{\overline{B}}_{r}^{\nu}\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}$

from which (30) immediately follows. $\quad\square$

The formula (30) clearly the nonholonomic contribution originating from the
functions $\alpha_{\nu}$ and the comparison with the holonomic case (25) is
evident: the additional terms appear only if the system undergoes nonholonomic
constraints, in the absence of which ${\cal L}^{*}={\cal L}$ and the usual
energy balance for holonomic systems is recovered.

###### Remark 5

It is worth quoting the following formulation of (31)

${\overline{B}}_{\nu}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)=\dfrac{d}{dt}({\overline{\alpha}}_{\nu}-\alpha_{\nu})-\sum\limits_{\mu=1}^{k}\dfrac{\partial\alpha_{\nu}}{\partial
q_{m+\mu}}({\overline{\alpha}}_{\mu}-\alpha_{\mu})+\dfrac{\partial\alpha_{\nu}}{\partial
t}$ (33)

the check of which is based on a calculation analogous to (32). The version
(33) shows more clearly the dependence of the coefficients on the deviation of
${\overline{\alpha}}$ from $\alpha$ and on the possible mobility of the
constraints, marked by the presence of the variable $t$ in the functions (1)
(rheonomic constraints).

###### Remark 6

If the forces of the system do not come from a potential then a procedure
analogous to the one developed leads to the equation of energy type

$\dfrac{d}{dt}\left(\sum\limits_{r=1}^{m}{\dot{q}}_{r}\dfrac{\partial
T^{*}}{\partial{\dot{q}}_{r}}-T^{*}\right)-\sum\limits_{\nu=1}^{k}({\overline{\alpha}}_{\nu}-\alpha_{\nu})\dfrac{\partial
T^{*}}{\partial
q_{m+\nu}}-\sum\limits_{\nu=1}^{k}{\overline{B}}_{\nu}\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}=-\dfrac{\partial T^{*}}{\partial
t}\sum\limits_{r=1}^{m}{\mathcal{F}}^{(r)}{\dot{q}}_{r}+\sum\limits_{\nu=1}^{k}{\mathcal{F}}^{(m+\nu)}{\overline{\alpha}}_{\nu}.$
(34)

which replaces (30); the functions ${\overline{\alpha}}$ and
${\overline{B}}_{\nu}$ are the same as (28) and (31). The mixed case of forces
partly attributable to a potential and partly not is clear.

###  2.1 The energy function

Let us introduce now the main question of this analysis. The energy-type
information (25) is formulated using only the independent variables
$q_{1},dots,q_{n}$. On the other hand, the expression for nonholonomic systems
(26) comes from equations (21) , first line, where the kinetic variables are
not independent and it could not be otherwise, since the deduced relationships
(5) are ignored. The result (29) makes us define the quantity in round
brackets in the left of the equality as the energy of the system. On the other
hand, the same comment can be done for (30) in order to define the energy in
terms of the restricted function ${\cal L}^{*}$ and to refer to the quantity
in round brackets as the energy of the system.

It is therefore worthwhile to examine and compare the two expressions that
emerge from (26) and (29), namely

${\cal
E}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{n},t)=\sum\limits_{i=1}^{n}{\dot{q}}_{i}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{i}}-{\cal L}$ (35)

and

${\cal
E}^{*}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)=\sum\limits_{r=1}^{m}{\dot{q}}_{r}\dfrac{\partial{\cal
L}^{*}}{\partial{\dot{q}}_{r}}-{\cal L}^{*}.$ (36)

in order to deal with the question of appropriately defining the energy of the
system and with the existence of possible first integrals, which originate
uniquely from the balance equations we have introduced above.

Let us start by pointing out that in general the two functions (35) and (36)
do not coincide: they can actually be compared once they are expressed in the
same variables, i. e. setting

${\cal
E}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)=\left.\left(\sum\limits_{i=1}^{n}{\dot{q}}_{i}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{i}}-{\cal
L}\right)\right|_{{\dot{q}}_{m+\nu}=\alpha_{\nu},\nu=1,\dots,k}$ (37)

that is replacing the dependent kinetic variables by means of (5); we indicate
again with ${\cal E}$ the recalculated function for simplicity. The function
${\cal E}^{*}$ calculated by the restricted function ${\cal L}^{*}$ does not
generally coincide with (37), whose calculation requires first the use of
${\cal L}$ and then the restriction (5): let us show in the following standard
situation the two calculations.

###### Example 1

Assume that the kinetic energy is of the standard form

$T(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{n},t)=\dfrac{1}{2}\sum\limits_{i,j=1}^{n}g_{i,j}{\dot{q}}_{i}{\dot{q}}_{j}+\sum\limits_{i=1}^{n}b_{i}{\dot{q}}_{i}+c$
(38)

where the coefficients $g_{i,j}(q_{1},\dots,q_{n},t)$ are the entries of a
$n\times n$ positive definite matrix and the functions
$b_{i}(q_{1},\dots,q_{n},t)$, $c(q_{1},\dots,q_{n},t)$ appear only in the case
where the holonomic constraints depend on time explicitly. Then (12)takes the
form

$T^{*}=\frac{1}{2}\left(\sum\limits_{r,s=1}^{m}g_{r,s}{\dot{q}}_{r}{\dot{q}}_{s}+\sum\limits_{\nu,\mu=1}^{k}g_{m+\nu,m+\mu}\alpha_{\nu}\alpha_{\mu}\right)+\sum\limits_{r=1}^{m}\sum\limits_{\nu=1}^{k}g_{r,m+\nu}{\dot{q}}_{r}\alpha_{\nu}+\sum\limits_{r=1}^{m}b_{r}{\dot{q}}_{r}+\sum\limits_{\nu=1}^{k}b_{m+\nu}\alpha_{\nu}+c$
(39)

and the computation of (36), in the validity of (20), leads to

$\displaystyle{\cal E}^{*}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum\limits_{r,s=1}^{m}g_{r,s}{\dot{q}}_{r}{\dot{q}}_{s}+\sum\limits_{\nu,\mu=1}^{k}g_{m+\nu,m+\mu}\alpha_{\mu}({\overline{\alpha}}_{\nu}-\frac{1}{2}\alpha_{\nu})$
$\displaystyle+$
$\displaystyle\sum\limits_{r=1}^{m}\sum\limits_{\nu=1}^{k}g_{r,m+\nu}{\overline{\alpha}}_{\nu}{\dot{q}}_{r}+\sum\limits_{\nu=1}^{k}b_{m+\nu}({\overline{\alpha}}_{\nu}-\alpha_{\nu})-c-U.$

On the contrary, the computation (35) provides

${\cal
E}=\frac{1}{2}\sum\limits_{r,s}^{m}g_{r,s}{\dot{q}}_{r}{\dot{q}}_{s}+\dfrac{1}{2}\sum\limits_{\nu,\mu=1}^{k}g_{m+\nu,m+\mu}\alpha_{\nu}\alpha_{\mu}+\sum\limits_{r=1}^{m}\sum\limits_{\nu=1}^{k}g_{r,m+\nu}\alpha_{\nu}{\dot{q}}_{r}-c-U$
(41)

The Example suggests that the two expressions coincides if
${\overline{\alpha}}_{\nu}=\alpha_{\nu}$, $\nu=1,\dots,k$: as a matter of
facts, the following property shows the deviation between (35) and (36).

###### Proposition 3

The functions (35) and (36) are related by the formula

${\cal E}^{*}={\cal
E}+\sum\limits_{\nu=1}^{k}({\overline{\alpha}}_{\nu}-\alpha_{\nu})\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}.$ (42)

Proof. We write

$\sum\limits_{i=1}^{n}{\dot{q}}_{i}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{i}}=\sum\limits_{r=1}^{m}{\dot{q}}_{r}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}\alpha_{\nu}.$

On the other hand, the chain differentiation

$\dfrac{\partial{\cal L}^{*}}{\partial{\dot{q}}_{r}}=\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{r}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}\qquad\textrm{for
any}\;\;r=1,\dots,m$

entails, by multiplying each relation by ${\dot{q}}_{r}$ and summing up,

$\sum\limits_{r=1}^{m}{\dot{q}}_{r}\dfrac{\partial{\cal
L}^{*}}{\partial{\dot{q}}_{r}}=\sum\limits_{r=1}^{m}{\dot{q}}_{r}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{r}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}\sum\limits_{r=1}^{m}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}{\dot{q}}_{r}=\sum\limits_{r=1}^{m}{\dot{q}}_{r}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{r}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}{\overline{\alpha}}_{\nu}$

owing to definition (28). Hence

${\cal E}^{*}=\sum\limits_{r=1}^{m}{\dot{q}}_{r}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{r}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}{\overline{\alpha}}_{\nu}-{\cal
L}^{*},\qquad{\cal E}=\sum\limits_{r=1}^{m}{\dot{q}}_{r}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{r}}+\sum\limits_{\nu=1}^{k}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}\alpha_{\nu}-\left.{\cal
L}\right|_{{\dot{q}}_{m+\nu}=\alpha_{\nu},\nu=1,\dots,k}.$

Since $\left.{\cal
L}^{*}(q_{1},\dots,q_{m},{\dot{q}}_{1},\dots,{\dot{q}}_{n},t)={\cal
L}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{n},t)\right|_{{\dot{q}}_{m+\nu}=\alpha_{\nu},\nu=1,\dots,k}$,
we finally deduce (42). $\quad\square$

###### Corollary 1

In the non-singularity situation $\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}\not=0$, $\nu=1,\dots,k$, ${\cal E}^{*}={\cal E}$
if and only if ${\overline{\alpha}}_{\nu}=\alpha_{\nu}$ for any
$\nu=1,\dots,k$.

###### Remark 7

The relation

$\dfrac{d{\cal E}^{*}}{dt}=\dfrac{d{\cal
E}}{dt}+\dfrac{d}{dt}\left(\sum\limits_{\nu=1}^{k}({\overline{\alpha}}_{\nu}-\alpha_{\nu})\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}\right)$

which is obtained by differentiating (42) leads back to (30) having in mind
(57), and vice versa: the check is based on (33) and on the rule

$\dfrac{\partial{\cal L}^{*}}{\partial y}=\dfrac{\partial{\cal L}}{\partial
y}+\sum\limits_{\nu=1}^{k}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}\dfrac{\partial\alpha_{\nu}}{\partial
y}\qquad\textrm{for}\;\;\textrm{any of the
variables}\;\;\;q_{m+1},\dots,q_{m+k},t$

The attention paid to the dependence of the various formulas on (28) shows
that the circustance $\alpha_{\nu}={\overline{\alpha}}_{\nu}$, $\nu=1,\dots,k$
is a special case to be this to be studied separately; Before doing this we
point out some notable situations, regarding either the structure of the
Lagrangian function or that of the constraint functions.

#### The stationary case

We identify the stationary case with systems verifying the assumptions (going
back to the general case (1))

$\dfrac{\partial{\cal L}}{\partial
t}=0,\qquad\dfrac{\partial\phi_{\nu}}{\partial t}=0,\quad\nu=1,\dots,k$ (43)

From the latter condition it follows that also the explicit functions
$\alpha_{\nu}$ do not depend on $t$, that is
$\dfrac{\partial\alpha_{\nu}}{\partial t}=0$, hence $\dfrac{\partial{\cal
L}^{*}}{\partial t}=0$, because of (22). The energy equations (29) and (30)
take the form

$\begin{array}[]{l}\dfrac{d{\cal
E}}{dt}=\sum\limits_{\nu=1}^{k}(\alpha_{\nu}-{\overline{\alpha}}_{\nu})\left(\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}}-\dfrac{\partial{\cal L}}{\partial
q_{m+\nu}}\right),\\\ \\\ \dfrac{d{\cal
E}^{*}}{dt}=-\sum\limits_{\nu=1}^{k}({\overline{\alpha}}_{\nu}-\alpha_{\nu})\dfrac{\partial{\cal
L}^{*}}{\partial
q_{m+\nu}}+\sum\limits_{\nu=1}^{k}\left(\dfrac{d}{dt}({\overline{\alpha}}_{\nu}-\alpha_{\nu})-\sum\limits_{\mu=1}^{k}\dfrac{\partial\alpha_{\nu}}{\partial
q_{m+\mu}}({\overline{\alpha}}_{\mu}-\alpha_{\mu})\right)\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}\end{array}$ (44)

As (25) recalls, for holonomic systems the energy is conserved in the
stationary case: the property cannot be extended to nonholomic systems, unless
the property ${\bar{\alpha}}_{u}=\alpha_{\nu}$ for each $\nu=1,\dots,k$ holds.

In the context of stationary systems a special case recurring in the early
examples of nonholonomic constraints concerns the additional assumptions

$\begin{array}[]{ll}\phi_{\nu}=\phi_{\nu}(q_{1},\dots,q_{m},{\dot{q}}_{1},\dots,{\dot{q}}_{m})&\textrm{for
each}\;\nu=1,\dots,k\;\textrm{and}\;j=1,\dots,m\\\ \\\ \ {\cal L}={\cal
L}(q_{1},\dots,q_{m},{\dot{q}}_{1},\dots,{\dot{q}}_{n}),&\end{array}$ (45)

that is the lagrangian coordinates $(q_{m+1},\dots,q_{n})$ corresponding to
the dependent velocities are absent. Equations (13) reduce to

$\dfrac{d}{dt}\dfrac{\partial{\cal
L}^{*}}{\partial{\dot{q}}_{r}}-\dfrac{\partial{\cal L}^{*}}{\partial
q_{r}}-\sum\limits_{\nu=1}^{k}\left(\dfrac{d}{dt}\left(\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}\right)-\dfrac{\partial\alpha_{\nu}}{\partial\dot{q}_{r}}\right)\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}=0\quad r=1,\dots,m$ (46)

and they extend the linear $\check{\rm C}$aplygin’s equations (see ([5]) and
[21]) to the nonlinear case. The evident advantage is that (46) contains only
the unknown functions $q_{1}$, $\dots$, $q_{m}$ and it is disentangled from
the constraints equations (16).

In the same context, assumptions (45) make (44) of the simplified form

$\begin{array}[]{ll}\dfrac{d{\cal
E}}{dt}=\sum\limits_{\nu=1}^{k}(\alpha_{\nu}-{\overline{\alpha}}_{\nu})\dfrac{d}{dt}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{m+\nu}},&\dfrac{d{\cal
E}^{*}}{dt}=\sum\limits_{\nu=1}^{k}\dfrac{d}{dt}({\overline{\alpha}}_{\nu}-\alpha_{\nu})\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}\end{array}$ (47)

which points up the contribution of each kinematic constraint in the variation
of the energy.

A further category with particular treatment due to the special form of the
Lagrangian function is represented by those systems where
$T=\dfrac{1}{2}\sum\limits_{i=1}^{n}M^{(i)}{\dot{q}}_{i}^{2}$ with $M^{(i)}$
positive constant (for instance mass of the point associated with the
$i$–coordinate) with the corresponding restricted form (12)

$T^{*}=\dfrac{1}{2}\sum\limits_{r=1}^{m}M^{(i)}{\dot{q}}_{r}^{2}+\dfrac{1}{2}\sum\limits_{\nu=1}^{k}M^{(\nu)}\alpha_{\nu}^{2}$
(48)

where $\alpha_{\nu}$, $\nu=1,\dots,k$ are the generic functions (5). It is
simple to verify that the equations of motion (13) are

$M^{(r)}{\mathop{q}\limits\limits^{..}}_{r}+\sum\limits_{\nu=1}^{k}M^{(m+\nu)}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}_{r}}}\dfrac{d\alpha_{\nu}}{dt}={\cal
F}^{(q_{r})}+\sum\limits_{\nu=1}^{k}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}{\cal
F}^{(q_{m+\nu})}\quad r=1,\dots,m$ (49)

with the obvious modification in case that (20) holds. In the latter case,
(35) and (36) are

${\cal
E}=\dfrac{1}{2}\sum\limits_{r=1}^{m}M^{(r)}{\dot{q}}_{r}^{2}+\dfrac{1}{2}\sum\limits_{\nu=1}^{k}M^{(\nu)}\alpha_{\nu}^{2}-U,\qquad{\cal
E}^{*}=\frac{1}{2}\sum\limits_{r=1}^{m}M^{(r)}{\dot{q}}_{r}^{2}+\sum\limits_{\nu=1}^{k}M^{(\nu)}\alpha_{\nu}\left({\overline{\alpha}}_{\nu}-\frac{1}{2}\alpha_{\nu}\right)-U,$

and an expressive way to write (29) or (30) is

$\dfrac{d}{dt}\left(\frac{1}{2}\sum\limits_{r=1}^{m}M^{(r)}{\dot{q}}_{r}^{2}-U\right)+\sum\limits_{\nu=1}^{k}M^{(\nu)}{\overline{\alpha}}_{\nu}\dfrac{d\alpha_{\nu}}{dt}=\sum\limits_{\nu=1}^{k}({\overline{\alpha}}_{\nu}-\alpha_{\nu})\dfrac{\partial
U}{\partial q_{m+\nu}}-\dfrac{\partial U}{\partial t}.$ (50)

The latter form draws attention to the energy contribution of the $\nu$–th
constraint on the applied force, absent if
$\overline{\alpha}_{\nu}=\alpha_{\nu}$.

In the next Example the function ${\overline{\alpha}}_{\nu}$ does not coincide
with $\alpha_{\nu}$.

###### Example 2

A peculiar instance of (48) concerns a point $P$ moving in the space and
constrained to mantain the same length of the velocity:

$|{\dot{P}}|=C(t)$ (51)

with $C$ nonnegative function of time. The example is analyzed by [29], [27])
for $C$ constant and by [27], [16] in the case $C(t)=1/\sqrt{t}$. In a
coordinate system $(q_{1},q_{2},q_{3})=(x,y,z)$ (51) is
${\dot{q}}_{1}^{2}+{\dot{q}}_{2}^{2}+{\dot{q}}_{3}^{2}-C(t)=0$, so $k=1$ $m=2$
and (5) is

${\dot{q}}_{3}=\alpha_{1}({\dot{q}}_{1},{\dot{q}}_{2},t)=\pm\sqrt{C^{2}(t)-{\dot{q}}_{1}^{2}-{\dot{q}}_{2}^{2}}.$

Assuming for the forces the potential $U(q_{1},q_{2},q_{3})$, evidently ${\cal
L}=\dfrac{1}{2}M({\dot{q}}_{1}^{2}+{\dot{q}}_{2}^{2}+{\dot{q}}_{3}^{2})+U$ and
(48) is $T^{*}=\frac{1}{2}MC^{2}(t)$, hence ${\cal L}^{*}=T^{*}+U={\cal
L}^{*}(q_{1},q_{2},q_{3},t)$. the equations of motion (49) are

$\left\\{\begin{array}[]{l}\dfrac{M}{R}\left((C^{2}-{\dot{q}}_{2}^{2}){\mathop{q}\limits^{..}}_{1}+{\dot{q}}_{1}({\dot{q}}_{2}{\mathop{q}\limits^{..}}_{2}-C{\dot{C}})\right)=\dfrac{\partial
U}{\partial q_{1}}\mp\dfrac{{\dot{q}}_{1}}{\sqrt{R}}\dfrac{\partial
U}{\partial q_{1}}\\\ \\\
\dfrac{M}{R}\left((C^{2}-{\dot{q}}_{1}^{2}){\mathop{q}\limits^{..}}_{2}+{\dot{q}}_{2}({\dot{q}}_{1}{\mathop{q}\limits^{..}}_{1}-C{\dot{C}})\right)=\dfrac{\partial
U}{\partial q_{2}}\mp\dfrac{{\dot{q}}_{2}}{\sqrt{R}}\dfrac{\partial
U}{\partial q_{2}}\end{array}\right.$

where we set $M=M^{(1)}$ and
$R({\dot{q}}_{1},{\dot{q}}_{2},t)=C^{2}(t)-{\dot{q}_{1}^{2}-{\dot{q}}_{2}^{2}}$.

Since ${\overline{\alpha}}_{1}=\mp\dfrac{C^{2}-R}{\sqrt{R}}$,
${\overline{\alpha}}_{1}-\alpha_{1}=\mp\dfrac{C^{2}}{\sqrt{R}}$,
${\overline{\alpha}}_{1}\dfrac{d\alpha_{1}}{dt}=-\dfrac{(C^{2}-R){\dot{R}}}{2R}$,
(50) ca be written as

$\dfrac{d}{dt}\left(\frac{M}{2}(C^{2}-R)-U\right)-\dfrac{M}{2}(C^{2}-R)\dfrac{\dot{R}}{R}=\mp\dfrac{C^{2}}{R}\dfrac{\partial
U}{\partial q_{3}}-\dfrac{\partial U}{\partial t}$

where $C^{2}-R={\dot{q}}_{1}^{2}+{\dot{q}}_{2}^{2}$. In the simpler case
$U=U(q_{1},q_{2})$, the expression we have written may help in order to solve
the equations of motion or give indications on the existence of possible first
integrals of the motion. In particular, for $C>0$ constant (50) reduces to
$\dfrac{d}{dt}(\frac{M}{2}y-U)-\frac{M}{2}\frac{y}{C^{2}-y}{\dot{y}}=0$, with
$y={\dot{q}}_{1}^{2}+{\dot{q}}_{2}^{2}$.

In a more general context, (51) can be considered part of the following
kinematic condition, involving each singular component:

$a{\dot{x}}^{2}+b{\dot{y}}^{2}+c{\dot{z}}^{2}=C^{2}(t)$

The case $a=b$, $c=-1$, $C(t)=0$ is examined in [11].

#### Linear kinematic constraints

Many examples and concrete applications show linear kinematic constraints.
indeed, most of the literature and theory only consider this case. In fact,
the most common examples of kinematic constraints such as the pure rolling of
a disk or a sphere on a plane or a body that slides frictionlessly with a
knife edge that slides only longitudinally can be formulated using linear
equations in the velocities and have been the first studied examples ([6],
[2], [31]).

In the case of linear constraints (15), we see in the explicit expressions
(16) that
$\alpha_{\nu}=\sum\limits_{r=1}^{m}\alpha_{\nu,r}{\dot{q}}_{r}+\beta\nu$ for
each $\nu=1,\dots,k$ so that the coefficients (28) and (31) are

${\overline{\alpha}}_{\nu}=\sum\limits_{r=1}^{m}\alpha_{\nu,r}{\dot{q}}_{r}=\alpha_{\nu}-\beta_{\nu},\qquad{\overline{B}}_{\nu}=\sum\limits_{j=1}^{m}\left(\dfrac{\partial\alpha_{\nu,j}}{\partial
t}-\dfrac{\partial\beta_{\nu}}{\partial
t}+\sum\limits_{\mu=1}^{k}\left(\dfrac{\partial\alpha_{\nu,j}}{\partial
q_{m+\mu}}\beta_{\mu}-\dfrac{\partial\beta_{\nu}}{\partial
q_{m+\mu}}\alpha_{\mu,j}\right)\right)$

so that (30) is (a similar observation may be done regarding (29))

$\dfrac{d{\cal E^{*}}}{dt}=-\dfrac{\partial{\cal L}^{*}}{\partial
t}-\sum\limits_{\nu=1}^{k}\beta_{\nu}\dfrac{\partial{\cal L}^{*}}{\partial
q_{m+\nu}}+\sum\limits_{\nu=1}^{k}\sum\limits_{j=1}^{m}\left(\dfrac{\partial\alpha_{\nu,j}}{\partial
t}-\dfrac{\partial\beta_{\nu}}{\partial
t}+\sum\limits_{\mu=1}^{k}\left(\dfrac{\partial\alpha_{\nu,j}}{\partial
q_{m+\mu}}\beta_{\mu}-\dfrac{\partial\beta_{\nu}}{\partial
q_{m+\mu}}\alpha_{\mu,j}\right)\right)\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}.$

The latter equation shows that, except for the possible explicit dependence on
$t$ of the terms, the presence of the functions $\zeta_{\nu}$ in (15), hence
of $\beta_{\nu}$ in (16), is the only responsible for the variation of (36).
In other words, if ${\cal L}$ does not depend explicitly on time and (15) are
linear homogeneous functions of ${\dot{q}}_{1}$, $\dots$, ${\dot{q}}_{n}$,
then the energy is conserved.

A special situation concerns the affine constraints of the form

$\alpha_{\nu}=\sum\limits_{j=1}^{m}a_{\nu,j}(q_{1},\dots,q_{n},t){\dot{q}}_{j}+c_{\nu}(t),\quad\nu=1,\dots,k$
(52)

with $c_{\nu}$ non zero function. In that case
${\bar{\alpha}}_{\nu}=\alpha_{\nu}-c_{\nu}$, so that
${\overline{B}}_{\nu}=\dfrac{\partial{\overline{\alpha}}_{\nu}}{\partial t}$
and (30) is

$\dfrac{d{\cal E^{*}}}{dt}=-\dfrac{\partial{\cal L}^{*}}{\partial
t}-c_{\nu}\dfrac{\partial{\cal L}^{*}}{\partial
q_{m+\nu}}+\sum\limits_{\nu=1}^{k}\left(\sum\limits_{j=1}^{m}\dfrac{\partial\alpha_{\nu,j}}{\partial
t}{\dot{q}}_{j}+{\dot{c}}_{\nu}\right)\dfrac{\partial
T}{\partial{\dot{q}}_{m+\nu}}.$ (53)

###### Example 3

For a particle in ${\mathbb{R}}^{3}$ with mass $M$ and submitted to the
constraint $ax{\dot{y}}+b{\dot{x}}y+c-{\dot{z}}=0$, $c\not=0$, we set
$(q_{1},q_{2},q_{3})=(x,y,z)$ so that (54) writes
${\dot{q}}_{3}=aq_{1}{\dot{q}}_{2}+bq_{2}{\dot{q}}_{1}+c$ (in this case $m=2$
and $k=1$); moreover

${\cal
L}^{*}=\frac{1}{2}M[(1+a^{2}q_{1}^{2}){\dot{q}}_{1}^{2}+(1+b^{2}q_{2}^{2}){\dot{q}}_{2}^{2}+2abq_{1}q_{2}{\dot{q}}_{1}{\dot{q}}_{2}+2c(aq_{1}{\dot{q}}_{2}+bq_{2}{\dot{q}}_{1})]+U(q_{1},q_{2},q_{3}).$

Whenever $U=U(q_{1},q_{2})$, (53) supplies the conservation of the quantity

$\frac{1}{2}M({\dot{q}}_{1}^{2}+{\dot{q}}_{2}^{2})-\frac{1}{2}M(aq_{1}{\dot{q}}_{2}+bq_{2}{\dot{q}}_{1}+c)(aq_{1}{\dot{q}}_{2}+bq_{2}{\dot{q}}_{1}-c).$

A necessary and sufficient condition (in terms of geometrical properties of
the constraint manifold) in order that the energy integral exists is discussed
and proved in [11].

###### Remark 8

More generally, for a positive integer $p$ the set of conditions

$\sum\limits_{j=1}^{n}\sigma_{\nu,j}(q_{1},\dots,q_{n},t){\dot{q}}_{j}^{p}+\zeta_{\nu}(q_{1},\dots,q_{n},t)=0,\quad\nu=1,\dots,k$
(54)

defines affine nonholonomic constraints of degree $p$; the linear case $p=1$
corresponds to (15). The explicit form (5) which generalizes (16) is

${\dot{q}}_{m+\nu}^{p}=(\pm
1)^{p+1}\left(\sum\limits_{j=1}^{m}\alpha_{\nu,j}(q_{1},\dots,q_{n},t){\dot{q}}_{j}^{p}+\beta_{\nu}(q_{1},\dots,q_{n},t)\right)^{1/p},\qquad\nu=1,\dots,k$
(55)

for suitable coefficients $\alpha_{\nu,j}$ and $\beta_{\nu}$. We deduce from
(55):

$\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{i}}=\dfrac{1}{p}\alpha_{\nu}\dfrac{\alpha_{\nu,i}{\dot{q}}_{i}^{p-1}}{\sum\limits_{j=1}^{m}\alpha_{\nu,j}{\dot{q}}_{j}^{p}+\beta_{\nu}}\qquad
i=1,\dots,m,\quad\nu=1,\dots,k$

hence

${\overline{\alpha}}_{\nu}=\dfrac{1}{p}\alpha_{\nu}\dfrac{\sum\limits_{i=1}^{m}\alpha_{\nu,i}{\dot{q}}_{i}^{p}}{\sum\limits_{j=1}^{m}\alpha_{\nu,j}{\dot{q}}_{j}^{p}+\beta_{\nu}}\qquad\nu=1,\dots,k$

so that (56) is valid if and only if

$(1-p)\sum\limits_{j=1}^{m}\alpha_{\nu,j}{\dot{q}}_{j}^{p}=p\beta_{\nu}.$

An evident solution is the case of linear homogeneous constraints ($p=1$,
$\beta_{\nu}=0$ for each $\nu=1,\dots,k$).

##  3 A relevant class of nonholonomic constraints

As it emerged, it is worth dwelling on the special case

${\overline{\alpha}}_{\nu}=\sum\limits_{r=1}^{m}{\dot{q}}_{i}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{r}}=\alpha_{\nu}\quad\textrm{for
each}\quad\nu=1,\dots,k$ (56)

It is evident that (34) and (30) deserve a distinctive treatment whenever (56)
holds; such an assumption plays the crucial role for the conservation of the
energy of the system. According to what observed in (17), the condition
depends only on the property of the functions $\phi_{\nu}$, $nu=1,\dots,k$,
and not on the choiche of the explicit functions $\alpha_{\nu}$.

We start by remarking that, whenever (56) holds, then the inconsistency
(asserted in Remark 1) between the virtual displacements and the velocities
compatible with a blocked configuration is overcome: indeed, in this case the
expression (11) matches up with the velocity (9) where the last term is
dropped. The property (56 )therefore extends to nonlinear nonholonomic
constraints what happens for holonomic systems and for nonholonomic systems
with linear homogeneous constraints, that is (15) with $\zeta_{\nu}=0$,
$\nu=1,\dots,k$. At the same time, (27) shows that the virtual power of the
constraints forces is null, as one can see also from (19).

Regarding the energy of systems fulfilling (56), we state the following
property, which pertains to (35), (36) and the corresponding energy equations.

###### Proposition 4

Assume that (56) holds for each $\nu=1,\dots,k$. Then
$\sum\limits_{i=1}^{n}{\dot{q}}_{i}\dfrac{\partial{\cal
L}}{\partial{\dot{q}}_{i}}=\sum\limits_{r=1}^{m}{\dot{q}}_{r}\dfrac{\partial{\cal
L}^{*}}{\partial{\dot{q}}_{r}}$ so that ${\cal E}^{*}={\cal E}$ and (29), (30)
reduce to

$\dfrac{d{\cal E}}{dt}=-\dfrac{\partial{\cal L}}{\partial
t}\;\;\textrm{or}\;\;\dfrac{d{\cal E}^{*}}{dt}=-\dfrac{\partial{\cal
L}^{*}}{\partial
t}+\sum\limits_{\nu=1}^{k}\dfrac{\partial\alpha_{\nu}}{\partial
t}\dfrac{\partial T}{\partial{\dot{q}}_{m+\nu}}=-\dfrac{\partial{\cal
L}}{\partial t}.$ (57)

In the stationary case (43), the quantity ${\cal E}={\cal E}^{*}$ is conserved
during the motion.

The statements are straightforward consequences (42), (33) (which implies
${\overline{B}}_{\nu}=\dfrac{\partial\alpha_{\nu}}{\partial t}$) and (44).

An remarkable consequence of (57) is the following

###### Corollary 2

Let $\alpha_{\nu}$, $\nu=1,\dots,k$ be functions satisfying (56). Then, the
quantity ${\cal E}={\cal E}^{*}$ is a first integral for the motion (13) if
and only if the Lagrangian function ${\cal L}$ does not depend explicitly on
time, that is $\dfrac{\partial{\cal L}}{\partial t}=0$.

If one refers to the expression (38) for $T$, the Corollary requires that
$g_{i,j}$, $b_{i}$ and $c$, $i,j=1,\dots,n$, do not depend on $t$ explicitly
and the same goes for $U=U(q_{1},\dots,q_{n})$; the conserved quantity is the
one calculated in (41):

$\frac{1}{2}\sum\limits_{r,s}^{m}g_{r,s}{\dot{q}}_{r}{\dot{q}}_{s}+\frac{1}{2}\sum\limits_{\nu,\mu=1}^{k}g_{m+\nu,m+\mu}\alpha_{\nu}\alpha_{\mu}+\sum\limits_{r=1}^{m}\sum\limits_{\nu=1}^{k}g_{r,m+\nu}{\dot{q}}_{r}\alpha_{\nu}-U-c=T^{*}_{2}-U-c$
(58)

where $T^{*}_{2}$ (whose definition is clear in (58)) collects the quadratic
terms. The absence of the terms containing the coefficients $b_{i}$,
$i=1,\dots,n$, is evidently in line with the Jacobi integral for holonomic
systems. In case of linear kinematic constraints (15) the quantity (58) is
(see also (16))

$\sum\limits_{r,s=1}^{m}\left(\frac{1}{2}g_{r,s}+\frac{1}{2}\sum\limits_{\nu,\mu=1}^{k}g_{m+\nu,m+\mu}\alpha_{\nu,r}\alpha_{\mu,s}+\sum\limits_{\nu=1}^{k}g_{r,m+\nu}\alpha_{\nu,s}\right){\dot{q}}_{r}{\dot{q}}_{s}-U-c.$
(59)

###### Example 4

Consider a system of two points $P_{1}$ and $P_{2}$ of mass $M_{1}$ and
$M_{2}$ constrained on a plane and whose velocities are perpendicular and the
veloocity of $P_{2}$ is parallel to the straight line joining $P_{1}$ with
$P_{2}$ (or, equivalently, the velocity of $P_{1}$ is perpendicular to the
joining line). The constraints equations are

${\dot{x}}_{1}{\dot{x}}_{2}+{\dot{y}}_{1}{\dot{y}}_{2}=0,\qquad(y_{1}-y_{2}){\dot{x}}_{2}-(x_{1}-x_{2}){\dot{y}}_{2}=0$

and the explicit formulation is, setting $q_{1}=x_{1}$, $q_{2}=x_{2}$,
$q_{3}=y_{1}$, $q_{4}=y_{2}$,

${\dot{q}}_{3}=-\dfrac{q_{1}-q_{2}}{q_{3}-q_{4}}{\dot{q}}_{1},\qquad{\dot{q}}_{4}=\dfrac{q_{3}-q_{4}}{q_{1}-q_{2}}{\dot{q}_{2}}$

Having in mind (16) with $k=2$ and $m=2$, we see that
$\alpha_{1,1}=-(q_{1}-q_{2})/(q_{3}-q_{4})$, $\alpha_{1,2}=\alpha_{2,1}=0$,
$\alpha_{2,2}=(q_{3}-q_{4})/(q_{1}-q_{2})$, $\beta_{1}=\beta_{2}=0$ and the
function (12) is

$T^{*}(q_{1},q_{2},q_{3},q_{4},{\dot{q}}_{1},{\dot{q}}_{2})=\dfrac{1}{2}M_{1}\left(1+\alpha_{1,1}^{2}\right){\dot{q}}_{1}^{2}+\dfrac{1}{2}\left(1+\alpha_{2,2}^{2}\right){\dot{q}}^{2}_{2}$

where $M_{1}$ and $M_{2}$ are the masses. Assuming that the two points are
connected by a spring exerting the force $-\kappa(P_{1}-P_{2})$ on $P_{1}$ and
the opposite one on $P_{2}$ ($\kappa$ positive constant) and including also
the gravitational force directed in the direction of decreasing $y$, the
Lagrangian function is

${\cal
L}=\frac{1}{2}M_{1}({\dot{q}}_{1}^{2}+{\dot{q}}_{3}^{2})+\frac{1}{2}M_{2}({\dot{q}}_{2}^{2}+{\dot{q}}_{2}^{2})-\frac{\kappa}{2}\left((q_{1}-q_{2})^{2}+(q_{3}-q_{4})^{2}\right)-g(M_{1}q_{3}+M_{2}q_{4})$

showing $\dfrac{\partial{\cal L}}{\partial t}=0$; the constant of motion (58)
is

$I(q_{1},q_{2},q_{3},q_{4},{\dot{q}}_{1},{\dot{q}}_{2})=T^{*}+\frac{\kappa}{2}\left((q_{1}-q_{2})^{2}+(q_{3}-q_{4})^{2}\right)+g(M_{1}q_{3}+M_{2}q_{4}).$

We see that the circumstance (56) brings the energy balance back to that of
holonomic systems (25), provided that the function ${\cal L}$ or ${\cal
L}^{*}$ is proper positioned in (57).

###### Example 5

We consider a system of two material points $(P_{1},M_{1})$ and
$(P_{2},M_{2})$ on a plane whose velocities are both orthogonal to the
straight line joining $P_{1}$ with $P_{2}$. In the present case it is $m=k=2$:
we set $(q_{1},q_{2},q_{3},q_{4})=(x_{1},x_{2},y_{1},y_{2})$ and write the
constraint equations (16) as
${\dot{q}}_{3}=\frac{q_{2}-q_{1}}{q_{3}-q_{4}}{\dot{q}}_{1}$,
${\dot{q}}_{4}=\frac{q_{2}-q_{1}}{q_{3}-q_{4}}{\dot{q}}_{2}$. Assuming that
the active forces verify (20), the Lagrangian function is ${\cal
L}=\frac{1}{2}M_{1}(q_{1}^{2}+q_{3}^{2})+\frac{1}{2}M_{2}(q_{2}^{2}+q_{4}^{2})+U(q_{1},q_{2},q_{3},q_{4})$
not depending on time and with respect to (59) one has
$g_{1,1}=g_{3,3}=M_{1}$, $g_{2,2}=g_{4,4}=M_{2}$, $g_{i,j}=0$ for $i\not=j$,
$\alpha_{1,1}=\alpha_{2,2}=\frac{q_{2}-q_{1}}{q_{3}-q_{4}}$,
$\alpha_{2,1}=\alpha_{2,2}=0$ so that the conserved quantity is

$I(q_{1},q_{2},q_{3},q_{4},{\dot{q}}_{1},{\dot{q}}_{2})=\frac{1}{2}(m_{1}{\dot{q}}_{1}^{2}+m_{2}{\dot{q}}_{2}^{2})\left(1+\left(\frac{q_{2}-q_{1}}{q_{3}-q_{4}}\right)^{2}\right)-U(q_{1},q_{2},q_{3},q_{4}).$

Let us comment the time dependence of the constraints (1): the linked
equalities in (57) make us realize that a first integral can be conferred upon
the system even though the nonholonomic constraints are not stationary, that
is the functions (5) depend explicitly on $t$ and so in turn
$\dfrac{\partial\alpha_{\nu}}{\partial t}\not=0$. Actually, (57) ensures that
the quantity (35) is conserved if and only if ${\cal L}$ does not depend on
time explicitly, even though the constraints (1) do: an example in this sense
can be sketched as follows.

###### Example 6

In the three-dimensional space a point $Q\equiv(x_{Q},y_{Q},z_{Q})$ is moving
according to the given relations $x_{Q}=x_{Q}(t)$, $y_{Q}=y_{Q}(t)$,
$z_{Q}=z_{Q}(t)$. A point $P\equiv(x,y,z)$ can move only pursuing $Q$, that
means its velocity has at any time the direction of the straight line joining
$P$ and $Q$: the condition
$({\dot{x}},{\dot{y}},{\dot{z}})\wedge(x_{Q}(t)-x,y_{Q}(t)-y,z_{Q}(t)-z)={\bf
0}$ gives the two independent constraints

$\left\\{\begin{array}[]{l}(y_{Q}(t)-y){\dot{x}}-(x_{Q}(t)-x){\dot{y}}=0,\\\
\\\ (z_{Q}(t)-z){\dot{x}}-(x_{Q}(t)-x){\dot{z}}=0\end{array}\right.$ (60)

The system consists of only the point $P$ and the lagrangian coordinates can
be $q_{1}=x$, $q_{2}=y$, $q_{3}=z$; the explicit form (16) is, wherever
$x\not=x_{Q}$,

$\left\\{\begin{array}[]{l}{\dot{q}}_{2}=\dfrac{y_{Q}(t)-q_{2}}{x_{Q}(t)-q_{1}}{\dot{q}}_{1}=\alpha_{1,1}(q_{1},q_{2},t){\dot{q}}_{1},\\\
\\\
{\dot{q}}_{3}=\dfrac{z_{Q}(t)-q_{3}}{x_{Q}(t)-q_{1}}{\dot{q}}_{1}=\alpha_{2,1}(q_{1},q_{3},t){\dot{q}}_{1}\end{array}\right.$
(61)

The kinetic energy of the system (formed only by the point $P$) is
$T=\frac{1}{2}M({\dot{q}}_{1}^{2}+{\dot{q}_{2}}^{2}+{\dot{q}}_{3}^{2})$, where
$M$ is the mass, and the form (12) is, owing to (61),

$T^{*}(q_{1},q_{2},q_{3},{\dot{q}}_{1},t)=\dfrac{1}{2}M{\dot{q}}_{1}^{2}\left(1+\dfrac{(q_{2}-\eta(t))^{2}+(q_{3}-\zeta(t))^{2}}{(q_{1}-\xi(t))^{2}}\right).$
(62)

It is immediate to check that the constraints (61 verify the assumption (56)).
Assuming (without losing generality) that there are no applied forces, we have
${\cal L}=T$ hence $\dfrac{\partial T}{\partial t}=0$ and the quantity (35)
(or (36)) is conserved even though the constraints are mobile, by virtue of
(57). This quantity is equal to ${\dot{q}}_{1}\frac{\partial
T^{*}}{\partial{\dot{q}}_{1}}=T^{*}$, hence

${\dot{P}}^{2}=\frac{2}{M}T^{*}={\dot{q}}_{1}^{2}\left(1+\dfrac{(q_{2}-\eta)^{2}+(q_{3}-\zeta)^{2}}{(q_{1}-\xi)^{2}}\right)$

is constant during the motion, that is the magnitude of the speed does not
vary over time. Such an information can support the resolution of the equation
of motion (13)

$\left(1+\dfrac{(q_{2}-\eta)^{2}+(q_{3}-\zeta)^{2}}{(q_{1}-\xi)^{2}}\right){\mathop{q}\limits^{..}}_{1}=-\dfrac{{\dot{q}}_{1}{\dot{\xi}}}{q_{1}-\xi}\dfrac{(q_{2}-\eta)^{2}+(q_{3}-\zeta)^{2}}{(q_{1}-\xi)^{2}}+\dfrac{{\dot{q}}_{1}}{(q_{1}-\xi)^{2}}\left({\dot{\eta}}(q_{2}-\eta)+{\dot{\zeta}}(q_{3}-\zeta)\right).$
(63)

We finally remark that for systems with $T^{*}$ as in (48), the validity of
(56) eliminates in (50) the extra-contribution of the applied forces (with
potential $U$); at the same time in that case
$\sum\limits_{\nu=1}^{k}M^{(\nu)}{\overline{\alpha}}_{\nu}\dfrac{d\alpha_{\nu}}{dt}=\dfrac{d}{dt}\left(\dfrac{1}{2}\sum\limits_{\nu=1}^{k}M^{(\nu)}\alpha_{\nu}^{2}\right)$,
hence (50) assumes the ordinary form

$\dfrac{d}{dt}\left(\frac{1}{2}\sum\limits_{i=1}^{m}M^{(i)}{\dot{q}}_{i}^{2}+\frac{1}{2}\sum\limits_{\nu=1}^{k}M^{(\nu)}\alpha^{2}_{\nu}-U\right)=\dfrac{d}{dt}(T^{*}-U)=-\dfrac{\partial
U}{\partial t}$ (64)

analogous to that pertinent to holonomic systems.

###  3.1 Some special categories of kinematic constraints

In order to intercept the categories of constraints that have this particular
treatment (57) it is essential, from the mathematical point of view, to
specify which are the functions fulfilling condition (56): a first simple but
considerable first step is the following

###### Proposition 5

The function
${\overline{\alpha}}_{\nu}=\sum\limits_{i=1}^{m}{\dot{q}}_{i}\dfrac{\partial\alpha_{\nu}}{\partial{\dot{q}}_{i}}$
coincides with $\alpha_{\nu}$ if and only if $\alpha_{\nu}$ is a positive
homogeneous function of degree $1$ w. r. t. ${\dot{q}}_{1}$, $\dots$,
${\dot{q}}_{m}$, that is

$\alpha_{\nu}(q_{1},\dots,q_{n},\lambda{\dot{q}}_{1},\dots,\lambda{\dot{q}}_{m},t)=\lambda\alpha_{\nu}(q_{1},\dots,q_{n},{\dot{q}}_{1},\dots,{\dot{q}}_{m},t)\quad\textrm{for
any}\;\;\lambda>0.$ (65)

Proof. It is merely the Euler’s homogeneous function theorem: the positively
homogeneous functions $f(x_{1},\dots,x_{N})$ of degree are exactly the
solutions of the partial differential equation
$\sum\limits_{i=1}^{N}x_{i}\dfrac{\partial f}{\partial x_{i}}=kf$; in the
present case, the variables $x_{1},\dots,x_{N}$ are the kinetic variables
${\dot{q}}_{1}$, $\dots$, ${\dot{q}}_{m}$ and $f$ is $\alpha_{\nu}$, for each
$\nu=1,\dots,k$. $\quad\square$

Shifting now the attention to the assigned functions (1), it is certainly
appropriate to ask how the functions $\phi_{\nu}$ must be so that the explicit
expressions verify (65): we can certainly identify at least the following
types of constraints:

(1)

Linear homogeneous kinematic constraints, that is (15) with $\zeta_{j}=0$:

$\sum\limits_{i=1}^{n}\sigma_{\nu,i}(q_{1},\dots,q_{n},t){\dot{q}}_{i}=0\qquad\nu=1,\dots,k$
(66)

(2)

Nonlinear kinematic constraints of type

$\sum\limits_{i,j=1}^{n}a_{i,j}^{(\nu)}(q_{1},\dots,q_{n},t){\dot{q}}_{i}{\dot{q}}_{j}=0,\qquad\nu=1,\dots,k$
(67)

that is homogeneous quadratic functions with respect to the kinetic variables.

The check of (65) for the linear case (66) is immediate; as for the case (2),
it can be seen by means of the implicit function theorem that the functions
(5) defined by (67) are actually homogeneous functions of degree $1$ with
respect tothe variables ${\dot{q}}_{1}$, $\dots$, ${\dot{q}}_{r}$.

Examples of the typology (67) are standard kinematic restrictions which are
frequently considered in applications and in the literature (we limit
ourselves to systems of only two points for ease of writing):

$\begin{array}[]{ll}{\dot{x}}_{1}^{2}+{\dot{y}}_{1}^{2}+{\dot{z}}_{1}^{2}-\left({\dot{x}}_{2}^{2}+{\dot{y}}_{2}^{2}+{\dot{z}}_{2}^{2}\right)=0&\textrm{same
magnitude of the velocities}\\\
{\dot{x}}_{1}{\dot{y}}_{2}-{\dot{x}}_{2}{\dot{y}}_{1}=0,\quad{\dot{x}}_{1}{\dot{z}}_{2}-{\dot{x}}_{2}{\dot{z}}_{1}=0&\textrm{parallel
velocities}\\\
{\dot{x}}_{1}{\dot{x}}_{2}+{\dot{y}}_{1}{\dot{y}}_{2}+{\dot{z}}_{1}{\dot{z}}_{2}=0&\textrm{perpendicular
velocities}\end{array}$

The mentioned systems have access to the balance equation (57) and the absence
of $t$ from the Lagrangian function ${\cal L}$ will provide the first integral
of motion (36).

For the examples listed above the coefficients $a_{i,j}^{(\nu)}$ in (67) are
constant: in the following final example, known as a nonholonomic pendulum,
the coefficients depend on the Lagrangian coordinates.

###### Example 7

Two points $P_{1}$ and $P_{2}$ are moving on a plane in a way that the
straight lines perpendicular to the velocities of the points ${\dot{P}}_{1}$,
${\dot{P}}_{2}$ intersect in one of the points of a given curve $\gamma$ lying
on the plane.

If $(x_{i},y_{i},0)$ are the coordinates of $P_{i}$, $i=1,2$, the two straight
lines perpendicular to the velocities intersect at the point of coordinates

${\bar{x}}=\dfrac{{\dot{y}}_{1}{\dot{y}}_{2}(y_{1}-y_{2})-x_{2}{\dot{y}}_{1}{\dot{x}}_{2}+x_{1}{\dot{x}}_{1}{\dot{y}}_{2}}{{\dot{x}}_{1}{\dot{y}}_{2}-{\dot{x}}_{2}{\dot{y}}_{1}},\;\;{\bar{y}}=\dfrac{{\dot{x}}_{1}{\dot{x}}_{2}(x_{2}-x_{1})-y_{1}{\dot{y}}_{1}{\dot{x}}_{2}+y_{2}{\dot{x}}_{1}{\dot{y}}_{2}}{{\dot{x}}_{1}{\dot{y}}_{2}-{\dot{x}}_{2}{\dot{y}}_{1}}$

so that, if the curve $\gamma$ is the graph of $y=g(x)$, the kinematic
constraint is formulated by

${\bar{y}}(x_{1},y_{1},x_{2},y_{2},{\dot{x}}_{1},{\dot{y}}_{1},{\dot{x}}_{2},{\dot{y}}_{2})=g({\bar{x}}(x_{1},y_{1},x_{2},y_{2},{\dot{x}}_{1},{\dot{y}}_{1},{\dot{x}}_{2},{\dot{y}}_{2}))$

Let us assume for simplicity that $\gamma$ is the straight line $ax+by=0$: the
constraint is

$a\left({\dot{y}}_{1}{\dot{y}}_{2}(y_{1}-y_{2})-x_{2}{\dot{y}}_{1}{\dot{x}}_{2}+x_{1}{\dot{x}}_{1}{\dot{y}}_{2}\right)+b\left({\dot{x}}_{1}{\dot{x}}_{2}(x_{2}-x_{1})-y_{1}{\dot{y}}_{1}{\dot{x}}_{2}+y_{2}{\dot{x}}_{1}{\dot{y}}_{2}\right)=0.$
(68)

Wih respect to (1), we have $n=4$ and $k=1$ (hence $m=3$): setting
$q_{1}=x_{1}$, $q_{2}=y_{1}$, $x_{2}=q_{3}=x_{2}$, $q_{4}=y_{2}$, the explicit
form of (68) is

${\dot{q}}_{4}=\dfrac{b(q_{1}-q_{3}){\dot{q}}_{1}+(bq_{2}+aq_{3}){\dot{q}}_{2}}{(aq_{1}+bq_{4}){\dot{q}}_{1}+a(q_{2}-q_{4}){\dot{q}}_{2}}{\dot{q}}_{3}=\alpha_{1}(q_{1},q_{2},q_{3},q_{4},{\dot{q}}_{1},{\dot{q}}_{2},{\dot{q}}_{3})$
(69)

The constraint (68) is of (67) type and we see that the explicit form (69) is
a homogeneous function of degree $1$ with respect to the kinematic independent
variables ${\dot{q}}_{1}$, ${\dot{q}}_{2}$, ${\dot{q}}_{3}$. Assuming that the
forces give rise to the potential $U(q_{1},q_{2},q_{3},q_{4})$, neither the
function

${\cal
L}^{*}=\dfrac{1}{2}M_{1}({\dot{q}}_{1}^{2}+{\dot{q}}_{2}^{2})+\dfrac{1}{2}M_{2}{\dot{q}}_{3}^{2}\left(1+\left(\dfrac{b(q_{1}-q_{3}){\dot{q}}_{1}+(bq_{2}+aq_{3}){\dot{q}}_{2}}{(aq_{1}+bq_{4}){\dot{q}}_{1}+a(q_{2}-q_{4}){\dot{q}}_{2}}\right)^{2}\right)+U(q_{1},q_{2},q_{3},q_{4})$

nor the constraint (69) depend explicitly on time $t$. The conserved quantity
coming from (57) is obtained by inverting the sign before $U$ in ${\cal
L}^{*}$.

##  4 Conclusion and upcoming research

The central theme of the first part of the work is the condition 18) that
although it dates back to over seventy years ago, it is accepted directly as
an axiom (anyway in accordance with physical measurements) rather than having
a provenance from the laws of mechanics.

if one follows the Lagrange multiplier procedure, this hypothesis ensures that
the constraint forces do not perform virtual work. If instead we turn to a
subset of independent velocities by means of (5, the same hypothesis provides
the set of virtual displacements useful for writing the equations of motion
(13).

The fundamental difference that brings systems with nonlinear constraints to a
more complicated level than holonomic systems and those with linear
constraints lies in the fact that for the latter ones the notion of
displacement overlaps with that of virtual velocity. On the contrary, the
formula (11) does not trace a virtual velocity seeing as the components
${\dot{q}}_{m+\nu}=\alpha_{\nu}$, $\nu=1,\dots,k$, are replaced by the
functions ${\overline{\alpha}}_{\nu}$ defined in (28). Such a discrepancy can
be approached to the analysis carried out in [23] (even though in a much more
elaborate formal context), where the concept of admissible or virtual
displacements are explored and discriminated.

Inevitably the standard procedure in the Lagrangian formalism of producing
energy-type information from the equations of motion, which we dealt with in
the second part of the work, entails also in this case a difference with the
holonomic or linear nonholonomic situation for the reason, essentially, that
the powers of forces cannot directly emerge in the energy equation (34) and in
the similar ones.

It is therefore legitimate to wonder how much the Hamiltonian type function
(36) built with only the independent velocities can be assumed as the energy
of the system: the discrepancy that we first reported means that this function
does not necessarily correspond to the energy (35) calculated with all the
variables.

The care we took in highlighting the dependence of the various formula (firt
of all (42)) on the differences ${\overline{\alpha}}_{\nu}-\alpha_{\nu}$ leads
us to conclude that the hypothesis (56)is essential and indicative for
treating the case of nonlinear nonholonomic constraints as a plain extension
of the linear case. Under the same assumption, the discrimination between
“virtual” and “admissible”, to say, does not exist.

The conclusion leads to further topics to investigate: which types of
constraints (1) make the hypothesis (56) hold? The mathematical conjecture is
the following: if all the functions $\phi_{n}u$ in (1), $\nu=1,\dots,k$, are
homogeneous functions (even if of different degree) with respect to the
kinetic variables ${\dot{q}}_{1}$, $\dots$, ${\dot{q}}_{n}$, then the explicit
functions (5) are homogeneous functions (with respect the independent
velocities) of degree one.

On the other hand, one may suppose that the zero set defined by (1) can be
always realized by means of homogeneous functions (it is known that the same
constraints, even kinematic ones, can be realized by different sets of
conditions).

In conclusion, the topic to be explored is whether the disentangling condition
(28) is a prerogative for systems with homogeneous nonholonomic constraints,
clearing the way for claiming that only this type of systems can show (under
further assumptions) the conservation of energy.

## References

  * [1] Appell P., Sur une forme nouvelle des èquations de la dynamique, J. Reine Angew. Math. , 121 , 301–319, 1899.
  * [2] Appell P., Les mouvements de roulement en dynamique, Paris: Scientia Phys. Math., 4, 1899.
  * [3] Bloch, A. M., Mardsen, J. E. and Zenkov, D. V., Quasivelocities and Symmetries in Non-Holonomic Systems, Dynamical Systems, 24, 187–222, 2009.
  * [4] Cameron, J. , M. , Book, W., Modeling Mechanisms with Nonholonomic Joints Using the Botzmann-Hamel Equations, The International Journal of Robotics Research 16(1), 47–59, 1997
  * [5] ${\check{\rm C}}$aplygin, S. A., On the theory of motion of nonholonomic systems. Theorems on the reducing multiplier, Mat. Sb. 28, 303–314, 1911.
  * [6] ${\check{\rm C}}$aplygin, S. A., On a rolling sphere on a horizontal plane, Mat. Sb. 24, 139-–168, 1903.
  * [7] ${\check{\rm C}}$etaev, N. G., Modification of the Gauss principle, Prikl. Mat. i Mekh. 5 n. 1 11–-12, 1941 (In Russian)
  * [8] ${\check{\rm C}}$etaev, N. G., On the Gauss Principles, Papers on Analytical Mechanics, 323, Science Academy, 1962.
  * [9] Rumyantsev, V. V., On the Chetaev principl Dokl. Akad. Nauk SSSR , 210 n. 4 787–790, 1973 (In Russian)
  * [10] Cronström, C., On the Compatibility of Nonholonomic Systems and Related Variational Systems, Acta Physica Universitatis comenianae, Vol. L–LI, 1 and 2, 25–36, 2010.
  * [11] Fassò, F. , Sansonetto, N., Conservation of Energy and Momenta in Nonholonomic Systems with Affine Constraints, Regular and Chaotic Dynamics, 20 no. 4, 449–462, 2015.
  * [12] Flannery M. R., The enigma of nonholonomic constraints, _American Journal of Physics_ 73, 265, 2005.
  * [13] Gantmacher, F. R., Lectures in analytical mechanics, Moskow: Mir Publisher, 1970.
  * [14] Hamel, G., Über nichtholonome Systeme, Math. Ann. 92, 31–41, 1924.
  * [15] Hamel, G., Die Lagrange–Eulersche Gleichungen der Mechanik, Z. Math. Phys. , 50, 1–57, Fortschritte 34, 757, 1904.
  * [16] Krupkovà, O., Mechanical systems with nonholonomic constraints, J. Math. Phys., 38, 5098–5126, 1997.
  * [17] de León, M., de Diego, D. M. , On the geometry of non-holonomic Lagrangian systems, J. Math. Phys. 37, 3389-–3414, 1996.
  * [18] Li, S. M. , Berakadar, J., A generalization of the Chetaev condition for nonlinear nonholonomic constraints: The velocity-determined virtual displacement approach, Reports on Mathematical Physics 63, Issue 2, 179–189, 2009.
  * [19] Lurie, A. I., Analytical Mechanics, Springer–Verlag, Berlin Heidelberg, 2002.
  * [20] Massa, E., Pagani, E., Classical dynamics of nonholonomic systems: a geometric approach, Ann. Inst. H. Poincaré Phys. Théor. 55 n. 1, 511–544, 1991.
  * [21] Ne${\check{\rm i}}$mark Ju. I., Fufaev N. A., Dynamics of Nonholonomic Systems, American Mathematical Society, Providence, 1972.
  * [22] Papastavridis, J. G., Analytical Mechanics: a Comprehensive Treatise on the Dynamics of Constrained Systems, World Scientific, 2014.
  * [23] Papastravidis, J. G., Time–integral Variational Principles for Nonlinear Nonholonomic Systems, J. Appl. Math., 64, 985–991, 1997.
  * [24] Cariñena, J. F., Rañada, M. F., Lagrangian systems with constraints: a geometric approach to the method of Lagrange multipliers, J. Phy2. A: Math. Gen. 26, 1335-1351, 1993.
  * [25] Rumyantsev, V. V., Forms of Hamilton’s principle for nonholonomic systems, Mechanics, Automatic and Robotics, 2, n. 10, 1035–1048, 2000.
  * [26] Salehani, M. K., A Jet Bundle Approach to the Variational Structure of Nonholonomic Mechanical Systems, Reports on Mathematical Physics 83, Issue 3, 373–385, 2019.
  * [27] Swaczyna, M., Several examples of nonholonomic mechanical systems, Communications in Mathematics 19, 27–56, The University of Ostrava, 2011.
  * [28] Talamucci, F., Rheonomic Systems with Nonlinear Nonholonomic Constraints: The Voronec Equations, _Regular and Chaotic Dynamics_ , 25 no. 6, 662–673, 2020.
  * [29] Virga, E., Un’osservazione sui vincoli anolonomi non perfetti, Riv. Mat. Univ. Parma, 13, 379–384, 1987.
  * [30] Voronec, P. V., On the equations of motion for nonholonomic systems, Mat. Sb. 22, 659–686, 1901.
  * [31] Voronec, P. V., On the equations of motion of a heavy rigid body rolling without sliding on a horizontal plane, Kiev. Univ. Izv. no. 11, 1–17, 1901.
  * [32] Zeković, D. N., Dynamics of mechanical systems with nonlinear nonholonomic constraints – II Differential equations of motion, Z. Angew. Math. Mech, 91 no. 11, 899–922, 2011.
  * [33] Zeković, D. N., Dynamics of mechanical systems with nonlinear nonholonomic constraints – III Analysis of motion, Z. Angew. Math. Mech, 93 no. 8, 550–574, 2013.
  * [34] Zeković, D. N., Energy integrals for mechanical ssystems with nonlinear nonholonomic constraints, Int. J. Non–Linear Mech. 73, 105–113, 2015.