Monserrat Ríos-Hernández

Juan Manuel

Adriana Herlinda Vilchis-González

Nabil Zemiti
email: nabil.zemiti@lirmm.fr

Miguel A Padilla-Castañeda

Virtual lumbar puncture simulators: where are we today?

The lumbar puncture is an important procedure used in different specialties in medicine, such as neurology and anesthesiology. This procedure is considered complex, and the medical students believe they have not enough experience and confidence to perform it. That is the reason why some researchers are interested in proposing to create solutions to provide a better training of this procedure. This review aims to present a comparison between 10 virtual lumbar puncture simulators developed since 2000, where nine of them use a haptic device to represent a realistic procedure. In addition, the physical and virtual components are discussed, including the haptic device, the haptic algorithm, the virtual environment and the evaluations of the simulators. The contribution of this review is to recognize the limitations and strengths to identify some optimal requirements for the design, development, or selection of the software/hardware for the next generation of virtual lumbar puncture simulators.

I. INTRODUCTION

The lumbar puncture (LP) is an invasive procedure where the physician inserts a thin needle into in the lumbar region, passing through different tissues, piercing the dura and accessing to the subarachnoid space [START_REF] Wright | Cerebrospinal fluid and lumbar puncture: a practical review[END_REF]. This procedure is used for treatment, such as drug delivery or cerebrospinal fluid (CSF) removal to decrease intracranial pressure; and for diagnosis through the CSF sampling [START_REF] Wright | Cerebrospinal fluid and lumbar puncture: a practical review[END_REF], [START_REF] Sargsyan | Augmented reality application to enhance training of lumbar puncture procedure: Preliminary results[END_REF].

LP is considered as a complex and difficult procedure, whose success is related to the knowledge of the contraindications, the anatomy of the zone, the steps of the procedure, the technique and the skills of the physician [START_REF] Moghtaderi | Lumbar puncture: techniques, complications and csf analyses[END_REF], it means, the success of the procedure is correlated with the physician experience. Following an LP, the patient can experience some adverse events, such as post procedural headache, hemorrhage, infections, nerve injury, CSF leak, cerebral herniation, backache, hematoma, implantation of epidermoid tumor, radicular pain or seizures [START_REF] Kroll | Trends in lumbar puncture over 2 decades: A dramatic shift to radiology[END_REF], [START_REF] Moghtaderi | Lumbar puncture: techniques, complications and csf analyses[END_REF].

The LP procedure consists of different steps that involved technical and non-technical skills. In general it is described as the insertion of the needle in an intervertebral space between the L2 and L5 vertebra. The main characteristic regarding this procedure is that the biomechanical behavior of the tissues helped the physician to distinguish the layers when the needle pass trough. These layers are the skin, subcutaneous fat, supraspinous and interspinous ligaments, ligamentum flavum, epidural space, and dura mater [START_REF] Mirbagheri | Introducing a robotic lumbar puncture simulator with force feedback: Lp sim[END_REF]. The biomechanical behavior can be described as a "pop" or loss of resistance when the needle pierces the dura mater and the needle is inside of the subarachnoid space, where the CSF is [START_REF] Färber | Virtual reality simulator for the training of lumbar punctures[END_REF].

Because this procedure in general has a negative perception from the population, the physician has the task of explaining LP and showing a positive and trusting attitude towards it, however some medical students perceived the LP as an important but complicated procedure, with not enough experience and confidence [START_REF] Cranach | Medical students' attitudes toward lumbar puncture-and how to change[END_REF].

Medical students can learn the LP with different techniques, and the lack of experience can generate error in the detection of each layer causing a wrong LP [START_REF] Mirbagheri | Introducing a robotic lumbar puncture simulator with force feedback: Lp sim[END_REF]. Some of the students learn this procedure through observation sessions as in the apprenticeship method, where the student with the guide of an expert, performs the LP directly on a patient in a clinical environment [START_REF] Färber | Virtual reality simulator for the training of lumbar punctures[END_REF]. This method has some disadvantages because the student is under stress for the possibility of make a mistake and cause pain or injury to the patient [START_REF] Färber | Virtual reality simulator for the training of lumbar punctures[END_REF].

Other technique is the use of mannequins or training models built with plastic or rubber, without mechanical or electronic parts. They are portable, easy to use and to set up, with anatomical landmarks for palpation, with them the students can repeat the procedure and can feel the loss of resistance of some tissues, however, they present some limitations such as the limited number of uses before replacement of the fake tissues (consumables), limitation in mimicking the exact biomechanical behavior of tissues, the assessment of the training and the limitation represented by pathological and anatomical variances of the lumbar anatomy [START_REF] Mirbagheri | Introducing a robotic lumbar puncture simulator with force feedback: Lp sim[END_REF], [START_REF] Färber | Virtual reality simulator for the training of lumbar punctures[END_REF], [START_REF] Vaughan | A review of epidural simulators: where are we today?[END_REF].

Virtual simulators have been developed in order to eliminate these limitations using technologies such as haptic devices and computer controls or virtual reality (VR) [START_REF] Vaughan | A review of epidural simulators: where are we today?[END_REF]. In general, the use of this technology allows the use of medical images for more precise anatomical models, and also for the representation of different anatomies or some physiological differences.

VR can be defined as the use of computational models that allows the user interaction with artificial environments and objects, it simulates the reality through interactive devices, and sends or receives information from the virtual environment (VE). The virtual reality is classified according to the level of immersion or presence [START_REF] Riener | Virtual reality in medicine[END_REF], [START_REF] Teklemariam | A case study of phantom omni force feedback device for virtual product design[END_REF]:

• Desktop VR: It is not immersive, but it is the most common and cheapest because it uses a monitor as the visualization system, it can also include a haptic device. • Semi-immersive VR: It transmits to the user a sense of immersion within the VE on a visualization system larger than the monitor. • Fully immersed VR: The user is fully immersed in the VE and isolated from the real world, this is achieved through devices such as the head mounted display (HMD). On the other hand, a haptic device is an electromechanical device that enables communication between the user and computers through the sensation of object manipulation. These devices provide haptic or tactile feedback in the form of force, vibration or with a propioceptive response, following the interaction with a virtual object [START_REF] Jafari | Haptics to improve task performance in people with disabilities: A review of previous studies and a guide to future research with children with disabilities[END_REF].

It is worth mentioning that in the literature there are review articles that analyze virtual medical simulators in a general way like in [START_REF] Wu | Virtual simulation in undergraduate medical education: A scoping review of recent practice[END_REF], or for specific areas, like in [START_REF] Vaughan | A review of epidural simulators: where are we today?[END_REF] where epidural simulators based on mannequins or virtual are presented. However, this review will focus exclusively on those simulators that were developed to simulate LP.

The contribution of this review is to identify the requirements that a virtual LP simulator should cover to improve the skills and training of this procedure. Moreover, this paper also provides a comparison and analysis for researchers interested in developing a LP simulator, showing the necessary characteristics for the design of the virtual environment, and the selection or design of the haptic device. In addition, 10 simulators developed since 2000 are mentioned, each simulator includes its VE description, the technologies implemented (haptic device and visualization system), the haptic algorithm, and if the simulators were tested for usability or they were involved in a user experiment.

For this review a bibliographic search on the electronic Google Scholar database using the date intervals 2000 to 2022 was performed with the next keywords "Lumbar puncture (LP) simulator, virtual simulator, haptic simulator", where only simulators with a virtual environment that represents a LP were selected. In addition, the papers must describe in general terms the features of this simulator, such as the VE and what devices were implemented for allowing the user interaction with the simulator. Narrative data were extracted from each paper using the categories presented in Table I of Section IV. It is worth mentioning that papers developed for other puncture procedures were excluded from this review, for instance, as the epidural anesthesia or analgesia.

The reminder of this paper is organized as follows, Section II presents the non-commercial LP simulators in chronological order. Section III presents the commercial LP simulator. Section IV analyzes and discusses the advantages and disadvantages of the simulators. And, Section V shows the conclusions.

II. NON-COMMERCIAL LP SIMULATORS

Gorman et al. [START_REF] Gorman | A prototype haptic lumbar puncture simulator[END_REF] describe the development of an LP simulator, the VE includes 3D models of the spine from L1 to the sacrum, the spinal canal area and a virtual needle. The simulation includes six tissues layers: skin, dermis or subcutaneous tissue, muscles, ligamentum flavum and dura. This system consists of a Phantom 1.5 desktop (from Sensable Technologies) as the haptic device, a mannequin, and a monitor. Authors attached a needle in the final effector of the haptic device, in order to improve the interaction between the user and the simulator (see Fig. 1). The haptic device provided force feedback according to the virtual tissue penetrated by the virtual needle. This force was calculated according to the angle, deep of the needle and penetrated tissue. Each tissue was modeled as a haptic box, allowing the user to feel the friction of the surface, the force of the puncture and the damping effect. The values of each haptic box can be modified by the user, allowing to feel when the needle pass through the skin or collides with a bone. This simulator provides the possibility to see through the tissues and has the option to record the training and play it more times. The authors wanted to add a system to measure and track performance, but there is no information on whether it was done or not. The evaluation of this simulator is not presented in the paper.

In 2009, Färber et al. [START_REF] Färber | Virtual reality simulator for the training of lumbar punctures[END_REF] developed a VR simulator with a VE with 2D and 3D scenes. The scenes (see Fig. 2) showed a virtual needle and the surfaces of anatomical structures with different colors and opacities for each organ and the anatomic landmarks can be identified by the trainee, additionally the VE shows annotations related to the anatomy. The scene was obtained from computed tomography data. In general, the scene could be rotated, panned, and zoomed with the mouse. This simulator integrates the 3DoF (degree of freedom) Sensable Phantom Premium 1.5 as haptic device, which gives force feedback according to a haptic algorithm. For the force Fig. 1. Haptic device, mannequin, needle and virtual environment of the simulator developed in [START_REF] Gorman | A prototype haptic lumbar puncture simulator[END_REF].

feedback authors integrated a proxy-based haptic rendering, that consists in a virtual spring damper that connects the end effector of the haptic device and a virtual proxy; and a proxy-based haptic volume rendering, in general the proxy moves perpendicularly to the direction of the image gradient vector, it means that the force can be calculated according to the information of density of the image. In addition, the simulator has three modes according to the experience of the user (beginner, advanced and professional), and provides an evaluation module that records the movement of the needle, and its path is compared with an optimal path, also the piercing of bones or structures at risk are evaluated in a score from 0 to 100. After the evaluation, the system could give feedback related to user performance and advised about how to improve the procedure. On the other hand, the system had some limitations such as the haptic device used cannot provide real force feedback for the collision with bones, also the intuitiveness of the system was limited because they use a monitor causing that the user is not immersed, and the end-effector of the haptic device did not represent a needle. 42 participants evaluated the system in two groups, the control and the training group. Both groups followed the same study protocol, where they studied didactic material, then an introduction to the LP was given to then start a tutorial for the use of the haptic device, then the training group trained with the LP simulator and the control group practiced other procedure; the next day both of them trained with the LP simulator. The participants belonging to the training group obtained better results and a higher success rate than the other group. Jiang et al. [START_REF] Jiang | Remote haptic collaboration for virtual training of lumbar puncture[END_REF] presented another solution for the LP training with remote haptic collaboration (see Fig. 3). In this paper the VE is not described, but the system worked with two haptic devices for a collaboration between a medical teacher and a student, both manipulated a local Phantom Desktop (haptic device), i. e. when the teacher manipulates his Phantom the student receives the information as force feedback with the second haptic device by means of a local network. This haptic device has 6DoF and provides force feedback in 3DoF, this feedback is calculated according to the tissue penetrated. For the haptic algorithm, the system used the mass spring system for the interaction between the needle and viscoelastic tissues, also used a friction model proposed by Karnopp and the clamping force. The user experience of this simulator was evaluated by 20 experienced physicians (10 pairs), at the end all of them succesfully perform the LP and felt the haptic force feedback in the ligament. In [START_REF] Day | Lumbar puncture trainer[END_REF], Day presented an LP trainer that consisted of a haptic device and a virtual environment for training simulation. For the VE, authors used the Unity3D game engine for the creation of a model created with computed tomography data, and consisted of the spine, pelvis, and torso (See Fig. 4). This simulator used the Novint Falcon with 3DoF (X, Y, Z directions) that provides force feedback proportional to the force used by the surgeon during the insertion of a needle. Authors highlighted the need of validate this simulator for training, however this evaluation is not mentioned.

In [START_REF] Wang | A surgical training system for four medical punctures based on virtual reality and haptic feedback[END_REF] Wang et al. developed a Surgical Training System that includes the LP and other three more punctures (thoracentesis, bone marrow puncture, abdominal paracentesis). Each type of puncture had its own models for the surgical instruments (See Fig. 5), and a specific human model with organs and tissue related to the type of puncture. The system considers several models like clinician and nurses, and several surgical scenes were designed. These scenes included operating rooms and a patient on the table, chair, or bed. The surgical trainer was developed in the Unity3D game engine with a server program developed in PHP, and an HMD was implemented as visualization tool for the surgery scene. Authors did not mention the haptic device used by the system, but they mentioned the use of the finite element method and Fig. 4. Virtual environment described in [START_REF] Day | Lumbar puncture trainer[END_REF]. a mass spring system for soft tissue modeling. The system records the number of completed procedures, the competition rate, average finished time, and the procedures. No evaluation was carried out for this simulator. In [START_REF] Jose | Bimanual haptic simulator for training hand palpation and lumbar puncture[END_REF] José et al. developed a bi-manual haptic simulator for LP and palpation was presented. The system consists of two haptic devices and a VE which consists of a virtual needle, a virtual hand and a virtual body with the spinal column from L1 to sacrum (see Fig. 6). One of the haptic devices was used for palpation implementing a Novint Falcon, this device controlled a 3D model of a hand; the other one was a Phantom used to control the movements of a 3D virtual needle. These devices were connected to a computer with a visual display, where the user could see the virtual environment that display a hand, a needle, and a human body model. For the palpation, authors implemented a lateral haptic model consisting of lateral surface force feedback, whose forces acted in a normal direction of the thumb; on the other hand for the puncture, they used a layered haptic model, where each layer is modeled as a haptic box, and the forces were calculated according to the characteristics of the material of the layer, the needle position, and orientation. This simulator was not been evaluated.

In 2018 [START_REF] Ali | Virtual simulation in enhancing procedural training for fluoroscopy-guided lumbar puncture: A pilot study[END_REF] a virtual fluoroscopy-guided lumbar puncture Fig. 6. Bi-manual haptic virtual lumbar puncture simulator [START_REF] Jose | Bimanual haptic simulator for training hand palpation and lumbar puncture[END_REF].

module was developed using the ImmersiveTouch platform (see Fig. 7). The VE consists of a 3D lower torso model (lower spine, ligamentum flavum, posterior paraspinal muscle groups, subcutaneous fat and skin) created from a segmentation of anonymized CT datasets, and a virtual spinal needle. This system integrates a haptic device and a hand tracking tool that enable interaction between the user and the module, also there are some foot pedals that control the simulation of fluoroscopic images, and with a button they can rotate the image and obtain oblique images. The authors did not mention in what consists the haptic algorithm, but they provide haptic force feedback so the users can feel the resistance of the skin, the "pop" when they pierced it, when the dura is penetrated and the impenetrability bones. This module was evaluated in a pilot study with five radiology trainers and a radiology nurse practitioner, the participants use the simulator for 30 minutes, with at least one LP with fluoroscopy, then they were asked to answer a survey, where the participants mentioned that the module replicated in a realistic way the anatomy and the fluoroscopy, five participants mentioned that also the anatomy was replicated with the haptic feedback, and finally the training overall obtained a score of 4.3 in a range of 3 to 5.

Fluoroscopic images as well as 3D images available

Foot pedals Haptic device Immersive Touch Fig. 7. Lumbar puncture module with fluoroscopy of Immersive Touch [START_REF] Ali | Virtual simulation in enhancing procedural training for fluoroscopy-guided lumbar puncture: A pilot study[END_REF].

In 2020 Mirbagheri et al. [START_REF] Mirbagheri | Introducing a robotic lumbar puncture simulator with force feedback: Lp sim[END_REF] developed the LP Sim, which consists of a software with a force module and a VE developed in Unity. For the VE authors created a geometrical 3D model of the lumbar area, through magnetic resonance and computer tomography, to obtain the soft tissues and the vertebraes (L2-L5). Its hardware consists of a Novint falcon with a modified end effector, a processor unit, a monitor and a base platform. The profile force obtained with the simulation was compared with a real procedure. The force obtained was according to the movement, if there is no movement, there is no force, also the initial point and direction of the penetration is fixed, and all the applied force will be in that direction. After the simulator was integrated, five anesthesiologists evaluate some subjective parameters such as graphical quality, needle motion tracking quality, and similarity with a real patient, the results showed that the graphical component and the needle movement quality were near of excellent. Fig. 8 shows all the components of LP Sim. In [START_REF] Carneiro | Serious game para treinamento de punC ¸Ão lombar utilizando realidade virtual[END_REF] Serra et al developed a simulator, as a serious game, using Unity3D and blender to create the VE through the SCRUM methodology. The VE consists in four phases, the objective of the first one, the presentation phase, is to familiarize the user with the simulator; the second one, the formation phase, shows passively anatomical and technical knowledge; the third phase, the learning evaluation, the formation is evaluated and the simulator indicates if the user is qualified to practice the LP; and finally, the practical phase, the user trains the LP. The authors modeled the vertebraes, a patient, a needle and a surgical table (see Fig. 9). The LP procedure is not explained in detail, and the paper does not mention if the simulator needs another tools like a haptic device or a controller. An advantage of this simulator is that an expert defines the parameter of the training for a personalized and incremental training, however, authors did not mention what parameters can be defined. On the other hand, this simulator has not been tested.

III. COMMERCIAL SIMULATORS

Nowadays there is only one commercial virtual LP simulator, and it is a mixed reality approach, using virtual reality and augmented reality. The name of this simulator is Sim&Care Lumbar Puncture Simulator and is manufactured by InSimo company [START_REF]InSimo. Lumbar puncture simulator[END_REF]. Its development was done in collaboration with the Rheumatology Department at Strasbourg University Hospital. Sim&Care represents different anatomies and simulates the resistance of the different tissues through a haptic device, but the company does not mention the haptic algorithm that the simulator uses. In the last version, the company incorporate the capacity to monitor the patient's emotions and reactions according to the procedure. According to the company's web page, there are different plans to purchase the simulator, it can be used with the institution equipment, can be rented or totally bought. The equipment of this simulator consists of a PC, an optical device, augmented reality glasses and an audio device (see Fig. 10). Although the characteristics of the required devices are not specified, in the demonstrations they use the Phanton Omni. In addition, the simulator includes metrics for a defined assessment, defined with the help of doctors, physicians to observe the user's performance, although it is not specified in what these metrics consist of. IV. DISCUSSION This section will compare the simulators presented in Sections II and III, mentioning the advantages and disadvantages of these systems. Table I shows the main characteristics of the virtual LP simulators like the type of haptic device, the DoF to provide force feedback, the models used in the VE and number of patient variations.

In order to achieve a more realistic experience, the use of haptic devices has been implemented and this can be observed in nine of the 10 simulators presented in this article, where [START_REF] Carneiro | Serious game para treinamento de punC ¸Ão lombar utilizando realidade virtual[END_REF] is the only article where the use of a haptic device is not mentioned. The haptic device used in each system varies since all of them uses different commercial devices. In four of the 10 simulators [START_REF] Day | Lumbar puncture trainer[END_REF], [START_REF] Mirbagheri | Introducing a robotic lumbar puncture simulator with force feedback: Lp sim[END_REF], [START_REF]InSimo. Lumbar puncture simulator[END_REF] the devices monitor and provide force feedback in 3DoF, which means that the translation of the virtual needle can be modified according to the translation of the haptic device, and the user receives the force in X, Y an Z axes. On the other hand, two simulators [START_REF] Gorman | A prototype haptic lumbar puncture simulator[END_REF], [START_REF] Jiang | Remote haptic collaboration for virtual training of lumbar puncture[END_REF] provide force feedback in 3DoF and allow 6DoF movements, providing the possibility of observing the current rotation of the needle according to the rotation of the haptic device and also to modify and represent the angle of the insertion of the needle inside the body, but this angle can be modified at any time since there is no force to fix it. Finally, one simulator [START_REF] Färber | Virtual reality simulator for the training of lumbar punctures[END_REF] uses a haptic device that provides force feedback and allows 6DoF movements, whose benefits are the possibility of fixed the angle of insertion of the needle and provide constraint forces in the end effector of the haptic device which will be reflected in fixing the angle of insertion and eliminate the possibility of change it when the needle is inside the body.

In general the use of a haptic device allows to provide an accurate force feedback, but this force will depend on the characteristics of these devices. One of the limitations mentioned in [START_REF] Färber | Virtual reality simulator for the training of lumbar punctures[END_REF] is the force rendered when the user collides with a bone, since the haptic device provides a maximum of 8.5 N as a peak, and the continuity of the interaction with a bone can not be simulated. However, the use of haptic devices provides the possibility of calculating the interaction forces between the patient and the needle, according to the biomechanical behavior of the tissues, or to imitate a force profile obtained from a real LP.

VE is the other component of a VR simulator, and it provides different tools and advantages. First, the use of a VE allows modeling of different tissues or areas of the human body, and the simulators presented in this article use a similar anatomical model, all of them focused in the lumbar area with the skin, bones and some soft tissues. Furthermore, four of the anatomy models [START_REF] Färber | Virtual reality simulator for the training of lumbar punctures[END_REF], [START_REF] Day | Lumbar puncture trainer[END_REF], [START_REF] Mirbagheri | Introducing a robotic lumbar puncture simulator with force feedback: Lp sim[END_REF], [START_REF] Ali | Virtual simulation in enhancing procedural training for fluoroscopy-guided lumbar puncture: A pilot study[END_REF] were created with medical images for a realistic representation of the anatomy, but the disadvantage is that these models represent an specific anatomy, because the data were collected from an specific person. Second, the use of virtual models offers the opportunity to modify specific zones or tissues, providing the capacity of represent different anatomies and conditions, like elderly, calcification of the ligamentum flavum, and obesity. In the comparison of the simulators, only two of them have the option of choose between two or more patients and therefore the possibility of training with two anatomical variances. Finally, another advantage in using a VE is the possibility of evaluating the skills of the user with quantitative data and not to rely in the observation of the experts. For instance, the number of successful punctures, the time of the training, the trajectory, the number of unwanted collisions (bones).

One of the advantages of the mannequins over the virtual simulators is the possibility to train the palpation, which is an important step in the LP. This step has not been considered for most simulators, and just one [START_REF] Jose | Bimanual haptic simulator for training hand palpation and lumbar puncture[END_REF] tried to replicate it with a haptic device, but in [START_REF] Gorman | A prototype haptic lumbar puncture simulator[END_REF] authors used a mannequin, although it is not clearly stated that palpation is its objective.

The evaluation of a simulator is important to validate its performance and its capacity of training medical students. The usability of some simulators [START_REF] Jiang | Remote haptic collaboration for virtual training of lumbar puncture[END_REF], [START_REF] Ali | Virtual simulation in enhancing procedural training for fluoroscopy-guided lumbar puncture: A pilot study[END_REF], [START_REF] Mirbagheri | Introducing a robotic lumbar puncture simulator with force feedback: Lp sim[END_REF] was evaluated, and from this evaluation authors could detect some improvements, such as enhance the haptic algorithm, refined the fluoroscopy images, and improve user experience. On the other hand, just one article [START_REF] Färber | Virtual reality simulator for the training of lumbar punctures[END_REF] describes the evaluation of the simulator to validate its functionality for training proposes. Finally, one of the 10 analyzed simulators is commercially available, Sim&Care [START_REF]InSimo. Lumbar puncture simulator[END_REF] by InSimo, with different prices and plans to acquire it.

V. CONCLUSIONS

The development and use of virtual simulators look for a realistic training, in order to overcome the advantages offered by the mannequins, such as the palpation, the possibility of leak CSF if the dura is pierced, the possibility of remove the skin and look inside, and the physical similarity with a patient [START_REF] Vaughan | A review of epidural simulators: where are we today?[END_REF]. The LP virtual simulators mentioned in this article have proposed solutions to meet these characteristics and to provide more, such as immersion and different biomechanical behavior in tissues.

The simulators presented in this article have shown different characteristics that try to simulate a realistic situation for a LP procedure, but exists the potential of improve these systems, to provide a realistic and functional training. If these characteristics are put together in a simulator it can result in a system with benefits for the patients and the physicians. The physicians could train in a safe and controlled environment, and it will be reflected in a successful procedure on patient.

After the comparison between these simulators, some characteristics were detected, and can be combined in order to develop a more realistic LP virtual simulator and provide a better training to the user in a controlled and realistic environment, resulting in a successful and safe procedure. An important point to consider and probably, the most noticeable difference with the mannequins is the palpation, this characteristic must be considered in training, not necessarily within the VE but as a complement to the simulator. The 3D models of the VE has to consider the lumbar zone with the landmarks and its tissues, these models have to be anatomically correct. In addition, the simulator should has the possibility of choose between patient variations, or the possibility of change the biomechanical behavior of the tissues in order to represent different situations. Also, the skin can change its opacity to visualize the inner part of the area, and a representation of the CSF leak should be integrated. Moreover, the use of a haptic device is required to provide a realistic experience, ideally a device that provides feedback in 6DoF should be used to represent the angle of insertion.

Finally, a user experience evaluation has to be performed to confirm the usability of the system, also, an user study must be carried out to confirm that the simulator works properly as a training tool. As mentioned above, researchers interested in developing a LP simulator can identify the best software and hardware for the design, development or selection of the virtual environment and the haptic device.

Fig. 2 .

 2 Fig. 2. Virtual lumbar puncture simulator with stereographic visualization [6].

Fig. 3 .

 3 Fig. 3. Virtual environment with remote haptic collaboration [14].

Fig. 5 .

 5 Fig. 5. Virtual environment with 3D surgical tools [16].

Fig. 8 .

 8 Fig. 8. Components of the simulator LP Sim [5].

Fig. 9 .

 9 Fig. 9. Serious game for the LP training [19].

Fig. 10 .

 10 Fig. 10. Lumbar puncture simulator Sim & Care [20].

TABLE I A

 I COMPARISON OF TEN VIRTUAL LUMBAR PUNCTURE SIMULATORS

	Simulator	Haptic Device	Modeled	Patient	3D Models	Visualiza-	Evaluation	Palpation	Adjustable	Patient	User ex-	Evalua-
		(DoF force	anatomy	model		tion	module		needle	varia-	perience	tion as a
		feedback/DoF)							insertion	tions	evalua-	training
									point		tion	tool
	Gorman et	Phantom 1.5	Spinal canal	Layers	Virtual needle	NM	Playback of			1	X	X
	al. (2000)	desktop (3/6)	area and	modeled	and patient		the training
	[13]		column, six		anatomy		record
			tissue layers
	Farber et	Phantom	Skin, bones	Created	Virtual needle	Monitor	Recording,	X		2	X
	al. (2009)	Premium 1.5	and soft	from CT	and patient		comparison
	[6]	(6/6)	tissues	Data	anatomy		with optimal
							path, score
							from 0 to
							100.
	Jiang et al.	Phantom	Skin,	NM	NM	Monitor	Immediately	X		1		X
	(2013)	desktop device	muscles,				feedback
	[14]	(3/6)	ligaments				give it by an
							expert
	Day	Novint Falcon	Spine, pelvis,	Created	Patient anatomy	NM	NM	X	X	1	X	X
	(2014)	(3/3)	and torso	from CT
	[15]			data
	Wang et	NM	Organs	Modeled	Patient, nurses,	Monitor	Number of	X	NM	1	X	X
	al. (2017)			according	physician,	and HMD	trains,
	[16]			to the	surgical		completion
				puncture	instruments,		rate, average
					operating rooms.		finished time
	José et al.	Phantom Omni	Spinal	NM	Human hand and	Monitor	X			1	X	X
	(2018)	/Novint Falcon	column (L1		lumbar needle
	[17]	(3/3)	to the
			sacrum)
	Fluoroscopy	ImmersiveTouch	Lower torso	Created	Virtual needle	Stereostopic	NM	X		NM		X
	Guided	Platform (NM)	model	from CT	and patient	interaction
	Lumbar			datasets	anatomy
	Puncture
	(2018)
	[18]
	LP Sim	Novint Falcon	Lumbar area	Created	Virtual needle an	Monitor	NM	X	X	1		X
	(2020) [5]	(3/3)		from CT	patient anatomy
				and MR
	Sim&Care	Phantom Omni	Skin, bones,	Not	Patient anatomy	HMD,	NM	X		Minimum	NM	NM
	(2020)	(3/3)	and tissues	mentioned		Aug-				4
	[20]		where the			mented
			needle passes			reality
	Serra et al.	NA	Vertebraes	NM	Chirurgical table	NM	X	X	NM	1	X	X
	(2021)				and lumbar
	[19]				needle
	NM = Not mentioned, NA = Not available

ACKNOWLEDGMENT

The authors would like to thank to the Secretary for Research and Advanced Studies office of the Universidad Autónoma del Estado de México, for financing the project no. 6482/2022CIB. M. R. H thanks to Consejo Nacional de Ciencia y Tecnología (CONACYT) for the master scholarship (CVU 1078220).


Alessandro Colotti
email: alessandro.colotti@ls2n.fr

Alexandre Goldsztejn
email: alexandre.goldsztejn@ls2n.fr

Practical stability and attractors of systems with bounded perturbations

The classical Lyapunov analysis of stable fixed points is extended to perturbed dynamical systems that may not have any fixed point due to perturbations. Practical stability is meant here to assess the convergence of such systems. This is achieved by investigating a parametric optimization problem encoding some worst-case Lie derivative. Key properties of this parametric optimization problem are formulated. The proposed framework is finally applied to a class of perturbed linear systems tracking a highly nonlinear reference.

I. INTRODUCTION

Lyapunov's method is one of the most widely used tools to assess the stability of dynamical systems [START_REF] Hafstein | Review on computational methods for Lyapunov functions[END_REF], [START_REF] Khalil | Nonlinear systems[END_REF]. Broadly applied for both linear and non-linear systems, it can be generalized to study the stability of a larger class of dynamical systems, such as uncertain and non-continuous systems [START_REF] Cortes | Discontinuous Dynamical Systems: A tutorial on solutions, nonsmooth analysis, and stability[END_REF]. Lyapunov's theory is intrinsically local in its formulation and, while it allows for the analytical assessment of global stability in general, evaluating the exact Region of Attraction (ROA) of an equilibrium is a hard problem, even when the dynamics are known. The ROA's estimation problem has been extensively studied in the literature [START_REF] Blanchini | Set invariance in control[END_REF] and we can find several well known approaches to the problem. Classically, we can estimate the ROA by evaluating the invariant sub-level sets of a given fixed Lyapunov function. In the case of polynomial systems, the evaluation is performed solving optimization problems with linear matrix constraints, exploiting the sum-of-squares relaxation in polynomial optimization [START_REF] Chesi | Solving quadratic distance problems: an LMI-based approach[END_REF], [START_REF] Hachicho | Estimating domains of attraction of a class of nonlinear dynamical systems with LMI methods based on the theory of moments[END_REF], [START_REF] Chesi | Domain of Attraction: Analysis and Control Via SOS Programming[END_REF]. It's possible to generalize this approach to allow for parametric Lyapunov functions [START_REF] Tibken | Computing the domain of attraction for polynomial systems via BMI optimization method[END_REF], [START_REF] Chesi | Domain of Attraction: Analysis and Control Via SOS Programming[END_REF]. These techniques can be generalized to compute the robust ROA of perturbed systems under bounded disturbances [START_REF] Goldsztejn | Estimating the robust domain of attraction for non-smooth systems using an interval Lyapunov equation[END_REF], [START_REF] Iannelli | Robust estimations of the Region of Attraction using invariant sets[END_REF].

All these methods require that there exists a stable and robust (i.e. perturbation-independent) equilibrium for the system. This is a major drawback, since in a large class of real-world problems this hypothesis is not satisfied. An important example of this is the problem of tracking a moving target with unknown dynamics, which is a common problem that cannot, a priori, allow for a robust equilibrium point. In these problems it could still be possible however to discuss the practical stability of the system, i.e. to assess whether the system, while technically unstable in theory, can be considered to be stable in practice. Consider, for instance, the case where the tracking system quickly reaches a small neighborhood of the target while never actually converging to it. Motivated by this, we want to relax the requirement of a robust equilibrium by using instead the concept of attractor. An attractor is a set of states which is positively invariant and attracting, i.e. such that the system converges to it (possibly in infinite time) when its initial state is close enough to the attractor. In the tracking system example discussed above, the attractor would be the small neighborhood of the target; this characterization, however, holds in more general settings.

The main contribution of this paper is to characterize the practical stability of perturbed systems with bounded perturbations by exploiting a given candidate Lyapunov function. The proposed definition of value function associated to the system and candidate Lyapunov function allows characterizing the system's attractors, hence extending the classical Lyapunov analysis of attracting fixed points. Our formulation does not require any underlying structure for the dynamics of the system, relaxing the classical hypothesis of polynomial dynamics and allowing for the study of a large class of nonlinear, perturbed systems.

The paper is organized as follows. In Section II we generalize the classical ROA's estimation problem to the identification of an attractor and the estimation of its ROA, and we show that this problem can be framed as a parametric optimization problem, using the candidate Lyapunov function's level sets as parameter. In Section III we study the continuity of the parametric optimization problem's solution, showing that, under common assumptions, we can expect it to be well-behaved. Finally, in Section IV we apply our results to discuss the practical stability of a class of perturbed linear systems.

II. PRACTICAL STABILITY OF PERTURBED SYSTEMS

We consider a perturbed system ẋ = f (x, w), where w ∈ W represents a time-dependent, unknown perturbation, with f continuous1 with respect to both x and w and W a compact set. We also define a candidate Lyapunov function, that we denote with V (x), to be a continuously differentiable function satisfying:

V (x) ≥ 0 and V (x 0 ) = 0 (1)
for some x 0 . For practical purposes, one can think of V to be a standard Lyapunov function (see, e.g., [START_REF] Khalil | Nonlinear systems[END_REF]) for the nominal system ẋ = f (x, 0). If f (x, 0) is a feedback system, we can then see our method as a characterization of the perturbation's effects on the nominal controller. Still, this is not a strict requirement for our analysis and the weaker hypotheses (1) suffice. Given f and V , we propose to define the function m(c) that associates to a target value c of the candidate Lyapunov function the worst-case Lie derivative on the corresponding level set:

m(c) = sup x∈Vc Φ(x) (2a) Φ(x) = sup w∈W ∇V (x) T f (x, w), (2b)
where V c = {x ∈ X : V (x) = c} is the level set with value c. We call this function the value function2 associated to f and V , or simply the value function when the associated system and candidate Lyapunov function are clear from the context. Sufficient conditions for the two supremum to be maximum are investigated in the next section. The obvious usefulness of the value function is summarized in the following two properties, valid under the typical assumption that the level sets of the candidate Lyapunov function involved in these properties are compact: 1) If m(c) < 0 inside the interval (c, c], for some 0 < c < c, then the sublevel set V ≤c is an attractor and the sublevel set V ≤c is inside its region of attraction, i.e. V ≤c is an estimated region of attraction for the attractor. 2) If m(c) < 0 inside (c, +∞) for some 0 < c, then the sublevel set V c is a global attractor. As discussed in the introduction, for a practically stable system we typically expect that m(c) > 0 inside (0, c * ) and m(c) < 0 inside (c * , +∞): provided that V is proper, i.e., has bounded level sets, all trajectories will converge inside the positively invariant set V c * . The usefulness of the value function defined here is of course balanced by the difficulty of computing it. Its formal evaluation is carried out for a class of perturbed linear system in Section IV.

In the general case, the quest of finding an interval where the value function is negative provides a strong temptation: computing its roots is of course a simpler and appealing approach. The roots of the value function can be characterized by using the Karush-Kuhn-Tucker (KKT) conditions applied to the optimization problem 3 (2): where ∇ 2 V denotes the Hessian of V , f x , f w the Jacobians of f w.r.t. x and w, g(w) represents the constraints on w (i.e. W = {w : g(w) ≤ 0}) and λ, µ ∈ R are the KKT multipliers. In addition to the usual KKT conditions, we consider c as a variable and we add the constraint V (x, w) = 0 in order to find a solution compatible with m(c) = 0. However, this simple intuition hides a trap: if the value function is not continuous, we might incur in a change of sign without passing through zero, rendering this approach useless. Thus, before we succumb to temptation, It is crucial to characterize the properties of f , V that guarantee the value function's continuity.

∇ 2 V (x)f (x, w) + f x (x, w) T + λI ∇V (x) = 0 (3a) f w (x, w) T + µI ∇g(w) = 0 (3b) V (x) -c = 0 (3c) µg(w) = 0 (3d) ∇V (x) T f (x, w) = 0, (3e)

III. CONTINUITY OF THE VALUE FUNCTION

As seen in Section II, the continuity of the value function plays an important role in choosing the strategy to evaluate it. We want to study the continuity of m(c) by investigating its upper semicontinuity and lower semicontinuity separately. Intuitively, the continuity property can be "split" into these two sub-properties, and we find that a function is continuous when both hold at the same time. For a scalar function (like m(c)), the following definition applies [START_REF] Aliprantis | Infinite dimensional analysis : a hitchhiker's guide[END_REF].

Definition 1: Let h : H ⊂ R → R ∪ {-∞, +∞}. Then:
1) The function h is called upper semicontinuous at c if, for any > 0, there exists δ > 0 such that:

|c -c| < δ =⇒ h(c) ≤ h(c) + (4)
2) The function h is called lower semicontinuous at c if, for any > 0, there exists δ > 0 such that:

|c -c| < δ =⇒ h(c) ≥ h(c) - (5)
A function h which is upper (lower) semicontinuous on its whole domain is simply called upper (lower) semicontinuous.

If it is both upper and lower semicontinuous, then it is a continuous function.

An example of the difference between upper and lower semicontinuity can be seen in Figure 1. While, technically, in this definition we ask the codomain of h to be the extended reals R∪{-∞, +∞}, in practice, for m(c), we can consider to have the set of real numbers R as codomain. Moreover, due to the positive definiteness of V , its domain is R + . We are ready to discuss our results. First of all, we show that, under typical assumptions on f and V , the function Φ(x) is continuous.

Proposition 1: Let V : X → R be a function of class C 1 , and let f : X × W → R n be continuous in x and w, where W is a compact set. Then, Φ(x) is a continuous function and the maximum is attained, i.e. Φ(x) can be written as a max instead of sup.

Proof: Being both ∇V (x) and f (x, w) continuous, the product ∇ T V (x)f (x, w) is continuous as well. Thus, by considering W as a constant (and so continuous) setvalued function with nonempty compact values, we can directly apply Berge's maximum theorem [1, Th. 17.31], which proves that Φ(x) is continuous. While promising, the continuity of Φ(x) is not sufficient to guarantee the continuity of the value function. Consider, for instance, the unperturbed, scalar system defined as:

ẋ = f (x) = -10x (0.1x 2 +5) 2 if x > 0 -2x if x ≤ 0 , (6)
and let:

V (x) = x 2 0.1x 2 +5 if x > 0 x 2 if x ≤ 0 (7)
be its candidate Lyapunov function. The plot of V is shown in Figure 2.

It is easy to verify the continuity and continuous differentiability of f and V , respectively. Moreover, we can see that f (x) = -V (x), meaning that the Lie derivative is V (x) = -(V (x)) 2 , which is negative for all x ∈ R\{0}. Even with a scalar and globally asymptotically stable system, we can see, in Figure 2, that the value function is not upper semicontinuous: this is due to a sudden change of the worstcase Lie derivative at the level of the horizontal asymptote.

The lower semicontinuity of the value function is not guaranteed either. Consider, for instance, the system:

ẋ = f (x) = -4x(x 2 -x -2), (8)
with candidate Lyapunov function:

V (x) = x 4 - 4 3 x 3 -4x 2 + 32 3 . (9)
Similarly to the previous example, we have a system of the form f (x) = -V (x). As we can see in Figure 3, in this case as well, the value function is discontinuous, in particular it is not lower semicontinuous: this is due to sudden changes of the worst-case Lie derivative at the levels of the local minimizer and local maximizer. Motivated by this, we want, in the following theorem, to independently study the upper and lower semicontinuity of On the left, plot of V (x) defined in [START_REF] Hachicho | Estimating domains of attraction of a class of nonlinear dynamical systems with LMI methods based on the theory of moments[END_REF]. On the right, the corresponding value function, where the solid black points denote the values of m at the discontinuities. the value function. The theorem's proof relies on the same foundations of the already introduced Berge's maximum theorem: as we will see, the upper and lower semicontinuity of m(c) depend on the continuity of Φ(x) and the hemicontinuity 4 of V c , seen as a set-valued function V c : R + ⇒ X . We have the following.

Theorem 1: Let V , f satisfy the same hypotheses of Prop. 1. Then, the following holds true:

1) If V is radially unbounded, then m(c) is upper semicontinuous and the maximum is attained, i.e. m(c) can be written as a max instead of sup. 2) If, for all x ∈ V c, ∇V (x) = 0, then m(c) is lower semicontinuous at c. Remark 1: Case 1 rules out situations illustrated by ( 6)-( 7), while Case 2 rules out situations illustrated by ( 8)- [START_REF] Hachicho | Estimating domains of attraction of a class of nonlinear dynamical systems with LMI methods based on the theory of moments[END_REF].

Proof: Due to Proposition 1, Φ(x) is continuous in the theorem's hypotheses. Thus, we focus in this proof on the upper and lower hemicontinuity of V c .

We start with the first statement. Being V continuous, then the preimage of a closed set through V is closed as well, meaning that V c is closed. Moreover, the radial unboundness of V implies that V c is bounded, for all c ∈ R + . Thus, by [3, Cor. 21], V c is upper hemicontinuous, which means that, by [START_REF] Aliprantis | Infinite dimensional analysis : a hitchhiker's guide[END_REF]Lemma 17.30], m(c) is a max and it is upper semicontinuous. The second statement's proof is proved in a similar fashion, albeit more involved. Firstly, we show that V c is inner hemicontinuous at c, which is a property equivalent to the lower hemicontinuity [START_REF] Kim | Introduction to Correspondences[END_REF]Prop. 23]. By its definition, we have that V c is inner hemicontinuous at c if, for all x ∈ V c, it holds true that:

∀{c k } → c, ∃x k ∈ V c k : x k → x (10)
We prove the inner hemicontinuity of V c by contradiction. Suppose that V c is not inner hemicontinuous at c. That means that there exists x ∈ V c such that:

∃{c k } → c : ∀x k ∈ V c k , x k → x. (11)
Consider the restriction of V , centered in x and along the direction of its gradient evaluated at x, i.e. the scalar function V (h) = V (x + ∇hV (x)). Notice that, by the continuity of 4 Hemicontinuity generalized the concept of continuity for set-valued functions. Intuitively, a set-valued function is not lower hemicontinuous if we have a sudden "explosion" of its value (e.g. a new connected component appearing), and similarly it is not upper hemicontinuous if we have a sudden "implosion". For a formal introduction see, e.g., [START_REF] Aliprantis | Infinite dimensional analysis : a hitchhiker's guide[END_REF].

V , V is continuous. Moreover, there exists an open interval H ⊂ R, with 0 ∈ H, such that V is strictly increasing over H. This is due to the fact that the derivative D V (h) evaluated at zero is:

D V (h)| h=0 = ∇V (x) 2 > 0, (12)
and, by the continuity of ∇V , there exists a neighbourhood of the origin (i.e. an open interval H containing zero) where ( 12) holds true, meaning that V is strictly increasing on H. Thus, V admits a continuous inverse V -1 on H. Its continuity means that:

∀{c k } → c, V -1 (c k ) → V -1 (c) = 0 ( 13
)
Let

h k = V -1 (c k )
. By defining the sequence {x k } k∈N as

x k = x + h k ∇V (x)
, it follows from Eq. ( 13) that:

∀{c k } → c, x k → x ( 14
)
which is in contradiction with Eq. 11, proving the inner (and, thus, lower) hemicontinuity of V c at c. By [1, Lemma 17.29], the lower hemicontinuity of V c at c implies that m(c) is lower semicontinuous at c, concluding the proof.

IV. APPLICATION TO LINEAR TRACKING SYSTEMS WITH BOUNDED VELOCITY TARGET

A. Upper bound of the value function for linear systems with bounded uncertain right hand side

We consider the linear system

ż = Az + w, ( 15
)
where A is a stable matrix (i.e. all eigenvalues have negative real part) and the uncertainty w is bounded inside W = {w ∈ R n : w ≤ M }. We consider a quadratic Lyapunov function V (z) = z T P z for the nominal system with P symmetric positive-definite (SPD). The Lie derivative for the system with nominal w = 0 is given by z T (P A + A T P )z, which is supposed to be negative-definite. Such a matrix P is usually obtained by choosing an arbitrary SPD matrix Q, typically Q = I, and solving the Lyapunov equation

P A + A T P + Q = 0. ( 16
)
The following theorem provides an upper bound for the value function associated to the system and the Lyapunov function defined above.

Theorem 2: Let P , Q be positive-definite matrices that Lyapunov's equation ( 16) holds. Then, m(c) ≤ m(c) with

m(c) = σ 1 c + 2M λ 1 √ c, (17)
where σ 1 is the greatest eigenvalue of -P -1 Q and λ 1 is the greatest eigenvalue of P . Furthermore, σ 1 < 0.

Proof: With these hypotheses, m(c) can be written as:

m(c) = max z∈Vc Φ(z). ( 18
)
Here we have Using Cauchy-Schwarz inequality and the bound on the norm of w, we obtain z T P w ≤ P z w ≤ M P z , with the equality holding when w is parallel to P z. Since W is a ball of radius M , such a w parallel to P z with norm M exists and we have Φ(z) = -z T Qz + 2M P z . This leads us to:

Φ(z) = max w∈W -z T Qz + 2z T P w (19) = -z T Qz + 2 max w∈W z T P w. (20)
m(c) = max z∈Vc -z T Qz + 2M P z (21a) ≤ max z∈Vc -z T Qz + 2M max z∈Vc P z , (21b)
where the optimization problem is split into two independent subproblems that can be solved independently.

For the second subproblem, we solve max z∈Vc P z 2 . Lagrange's first-order conditions yield P 2 z -µP z = 0, which must hold for some µ ∈ R. Being P invertible, this is equivalent to solving P z = µz, leading to max z T P z=c P z 2 = max z T P z=c µ∈Λ(P )

µz T P z = λ 1 c, (22)
where Λ(P ) is the spectrum of P and λ 1 = max {Λ(P )} is the greatest eigenvalue of P . The second subproblem maximum is therefore √ λ 1 c. For the first subproblem, Lagrange's first-order conditions of max z∈Vc -z T Qz lead to -Qz -µP z = 0 for some µ ∈ R. We notice that the values of µ that satisfy this equation are the generalized eigenvalues [START_REF] Beresford | The Symmetric Eigenvalue Problem[END_REF] of the symmetric matrix pencil (-Q, P ), whose set we denote with Λ(-Q, P ). Being P positive-definite, we know that we have n real generalized eigenvalues (see [START_REF] Beresford | The Symmetric Eigenvalue Problem[END_REF]Th. 15.3.3]), that we denote, w.l.o.g., with σ 1 ≥ • • • ≥ σ n . Thus, by substitution, we have that:

max z T P z=c -z T Qz = max z T P z=c µ∈Λ(-Q,P ) µz T P z = σ 1 c, (23)
which allows us to write:

m(c) ≤ σ 1 c + 2M λ 1 √ c. (24)
We conclude the proof by noticing that, being P invertible, the generalized eigenvalues of the pencil (-Q, P ) are equal to the eigenvalues of -P -1 Q, and that, being -z T Qz < 0 for all z = 0, we have that σ 1 c < 0 for c = 0, proving that σ 1 is indeed negative.

c * = 4M 2 λ 1 σ 2 1 , (25)
and is negative for greater values of c. The positive definite quadratic Lyapunov function being radially unbounded, the sublevel set V ≤c * is proved to be a global attractor of the uncertain system (15). In the typical case where Q = I we have σ 1 = -1 λ1 and the root's expression simplifies to

c * = 4M 2 λ 3 1 . (26)

B. Linear tracking systems with bounded velocity target

We consider a simple linear system with known input matrix gain A of the form:

ẋ = Au, (27)
where x ∈ R n is the system's state, u ∈ R n its input and A ∈ R n×n is stable. The target reference signal r(t) ∈ R n is measured at all times, and we assume the knowledge of an upper bound M on ṙ(t) . The proportional control u(t) = -k(x(t) -r(t)), k > 0, leads to the following closed loop system:

ẋ(t) = kA x(t) -r(t) . ( 28
)
The dynamic of the tracking error z(t) = x(t) -r(t) is then

ż(t) = kAz(t) -ṙ(t). ( 29
)
We use Theorem 2 to find a globally attracting sublevel set V ≤c * . In order to compare the sublevel sets for different gains, we consider the Lyapunov function V (z) = z T P z with P A+A T P +I = 0 independently of the gain k. Therefore, the matrix Q that satisfies the Lyapunov equation ( 16) is Q = kI. Since Q is proportional to I, the sublevel value (25) now simplifies to

c * = 4M 2 λ 3 1 k 2 . ( 30
)

C. Numerical application

We consider the following highly oscillating target, r(t) = cos(2t) + 2 sin(2t) + 0.5 sin(3t) sin(t) -cos(3t) , .

(32)

The largest eigenvalue of P is λ 1 = 1 28 (15 + √ 85) ≈ 0.86. The globally attractive sublevel values obtained using Theorem 2 for different values of the gain k are given in Table I. The trajectories of the system and the globally attracting sublevel sets obtained by Theorem 2 are shown in Figure 6, while in Figure 7 we can see the value of the Lyapunov function with respect to time.

D. Application of Theorem 1

The positive definite quadratic Lyapunov function satisfies the hypothesis of Theorem 1, which proves that the exact value function m(c) is continuous. We can therefore infer the sign of m(c) by computing its roots. The system (3) that characterizes its roots is now -Qx + P w + λP x = 0 (33a) P x + µw = 0 (33b)

x T P x = c (33c) µ(w T w -M 2 ) = 0 (33d) -x T Qx + 2x T P w = 0. (33e)
It is polynomial of degree three 5 and can be solved using formal computations. For each gain value, we find four solutions showing two distinct values of c, given in Table II.

The solutions with c 2 correspond exactly to the solutions previously computed, while c 1 < c 2 . We therefore expect the latter to be local maximizers. This is confirmed by evaluating m(c 1 ) and m(c 2 ) as follows: for these values of c, we solve the system (33) with the last equation removed. This is a square system of equations6 that encodes Lagrange first order conditions for the Lie derivative, its solutions therefore include the global maximizer. The global maximizers obtained this way are shown in Table II and confirm that c 2 is indeed the unique positive root of m(c). Surprisingly, both m(c) computed here and m(c) computed in the previous section have the same unique positive root (compare Table I and Table II), while the upper bound (21) should entail some overestimation. This coincidence is left for future investigations.

V. CONCLUSION

In this paper, we introduced the concept of value function associated to a perturbed dynamical system and a candidate Lyapunov function. We have shown that by analyzing this function it is possible to identify and study the system's attractors, even in cases where the system lacks robust equilibria and when its vector field is nonlinear. In applications, this translates in the ability to assess the practical stability of a system, providing quantitative information on the trajectories' proximity to the unperturbed equilibrium and allowing for the evaluation of a region of attraction's estimate. On the downside, the value function is difficult to compute, requiring ad-hoc, non-trivial theoretical results even in the case of linear systems with a simple perturbation structure.

However, Theorem 1 opens the way to using root finding algorithms for locating the roots of m(c), deducing intervals of negative values and identifying the corresponding attractors. The use of numerical solvers and computer-assisted proof methods to perform this evaluation will be the subject of future studies. In particular, the use of interval analysis [START_REF] Goldsztejn | Estimating the robust domain of attraction for non-smooth systems using an interval Lyapunov equation[END_REF], which do not require any particular structure for the system's vector field, to find the value function's root in a validated way seems promising and is currently under investigation.

Even though we specialize the value function's definition for perturbed systems, it can clearly be generalized to allow for a generic differential inclusion. Thus, it appears logical to try to extend this approach to a larger class of systems, naturally defined as differential inclusions (for instance, systems with discontinuous vector fields). Moreover, by using the notion of set-valued Lie derivative, it would be possible to relax the continuous differentiability of V (x), allowing for Lipschitz continuous candidate Lyapunov functions, which are necessary in several applications (e.g. in sliding mode control [START_REF] Colotti | New convergence conditions for the super twisting algorithm with uncertain input gain[END_REF]).

Fig. 1 .

 1 Fig. 1. Example of upper (on the left) and lower (on the right) semicontinuous functions, where the solid black point denotes h(1). In practice, the only difference is whether h(1) is at the top or bottom portion of the discontinuity.

Fig. 2 .

 2 Fig.2. On the left, plot of V (x) defined in[START_REF] Cortes | Discontinuous Dynamical Systems: A tutorial on solutions, nonsmooth analysis, and stability[END_REF]. The dotted line represents the horizontal asymptote for x → +∞. On the right, the corresponding value function, where the solid black point denotes the value of m at the discontinuity.

Fig. 3 .

 3 Fig. 3.On the left, plot of V (x) defined in[START_REF] Hachicho | Estimating domains of attraction of a class of nonlinear dynamical systems with LMI methods based on the theory of moments[END_REF]. On the right, the corresponding value function, where the solid black points denote the values of m at the discontinuities.

Fig. 4 .

 4 Fig. 4. Graph of -c + √ c.

Fig. 5 .

 5 Fig. 5. Left: r(t); right: r (t) with a circle of radius 6.

Fig. 6 .

 6 Fig. 6. Tracking error of the proportional linear tracking system for the same initial condition, whose initial error is denoted with a black point, together with the globally attracting ellipsoids, for different values of the gain k ∈ {5, 10, 50, 100}, with western reading direction.

Fig. 7 .

 7 Fig. 7. Lyapunov function's value along the trajectories of the proportional linear tracking system (in blue) together with the globally attracting level set's value (in orange), for different values of the gain k ∈ {5, 10, 50, 100}, with western reading direction.

TABLE I

 I Values of c * obtained using Theorem 2 for different gain values.

TABLE II

 II Values of c for the extremal values leading to a zero Lie derivative, with the corresponding values of the value function.

	k	5	10	50	100
	c 1	0.0507	0.0127	0.000507 0.000127
	c 2	3.73	0.932	0.0373	0.00931
	m(c 1 )	2.21	1.10	0.221	0.110
	m(c 2 )	0	0	0	0

Lipschitz continuity with respect to x would enforce uniqueness of trajectories, but is not required for Lyapunov-based stability analysis.

The theoretical machinery, i.e., Berge's maximum theorem and the hemicontinuity of set valued functions, used throughout the paper comes mainly from economics, the field in which the name value function has been defined.

By substituting the definition of Φ(x) into (2a), we can see m(c) as the solution of a single optimization problem in x and w.

The only degree three monomial in the system is the complementarity constraint (33d), which is actually a simple alternative. The system can therefore be solved as two degree two systems.

With respect to the square system (33), the variable c is now fixed, and one equation has been removed, hence leading to a new square system.

VI. ACKNOWLEDGMENTS

This work was supported partially by the French Agence National de la Recherche through the grants ANR-16-CE33-0024 and ANR-18-CE33-0011.