Raphael A C Viera
email: raphael.viera@emse.fr

Jean-Max Dutertre
email: dutertre@emse.fr

Philippe Maurine
email: philippe.maurine@lirmm.fr

Rodrigo Possamai Bastos
email: rodrigo.bastos@univ-grenoble-alpes.fr

Standard CAD Tool-Based Method for Simulation of Laser-Induced Faults in Large-Scale Circuits

Designing secure integrated systems requires methods and tools dedicated to simulating -at early design stages-the effects of laserinduced transient faults maliciously injected by attackers. Existing methods for simulation of laser-induced transient faults do not take into account IR drop effects that are able to cause timing failures, abnormal reset, and sram flipping. This paper proposes a novel standard cad tool-based method allowing to simulate laser-induced faults in large-scale circuits. Thanks to a power-grid network modeled by a commercial IR drop cad tool, an additional transient current component causing laser-induced IR drop is taken into consideration. This current component flows from v DD to g ND and may have a significant effect on the fault injection process. The method provides fault sensitivity maps that enable a quick assessment of laser-induced fault effects on the circuit under analysis. As shown in the results, the number of induced faults is underestimated by a factor as large as 3.1 if laser-induced IR drop is ignored. This may lead to incorrect estimations of the fault injection threshold, which is especially relevant for the design of countermeasure techniques for secure integrated systems. Simulation times regarding four different circuits are also presented in the results section.

INTRODUCTION

Lasers have been used since the 1960s in order to emulate the effects caused by radiation on semiconductors [START_REF] Habing | The Use of Lasers to Simulate Radiation-Induced Transients in Semiconductor Devices and Circuits[END_REF]. In the early 2000s, [START_REF] Skorobogatov | Optical Fault Induction Attacks[END_REF] reported the use of laser illumination to induce faults Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from permissions@acm.org. into secure integrated circuits, e.g., a bit-flip into a sram cell. This created a need for designing robust circuits against laser fault injection, consequently generating a demand for simulation tools capable of simulating the effects of laser shots on ics. Although fault simulations can be performed at different abstraction levels of the design flow, i.e. transistor, gate, rtl and software, low abstraction levels provide the highest accuracy. At the electrical level, a double exponential current source has been demonstrated efficient for modeling a laser shot [START_REF] Jordan | Photoeffect on diffused P-N junctions with integral field gradients[END_REF][START_REF] Wirth | The Transient Response of Transistors and Diodes to Ionizing Radiation[END_REF]. This current source is added to the netlists of cells illuminated by the laser. Then an electrical level simulation, which takes into account the effects of the laser attack, can be performed.

The idea commonly accepted is that a laser shot generates parasitic currents [START_REF] Johnston | Charge generation and collection in p-n junctions excited with pulsed infrared lasers[END_REF]. These currents generate an undesired transient voltage that propagates through the logic toward the inputs of registers (D-type Flip Flops) and, if it is still present when the clock edge occurs, bits may be inverted, producing soft errors (SE). Due to the increasing transistor density, a laser shot will affect multiple gates at the same time. Thus, laser illumination also induces, in addition to the well known photoelectric effect, an IR drop phenomenon with a significant effect on the fault injection process that has to be taken into account while simulating laser fault injection [START_REF] Viera | Role of Laser-Induced IR Drops in the Occurrence of Faults: Assessment and Simulation[END_REF]. These effects must be simulated at low abstraction levels taking into account the layout topology to better represent physical phenomenon in the scope of a whole system, i.e., the simulation must be performed in complex circuits and not just in one (or few) cells.

To the best of our knowledge, among the existent fault simulators [START_REF] Hungse Cha | A gatelevel simulation environment for alpha-particle-induced transient faults[END_REF][START_REF] Greenstein | E-PROOFS: A CMOS bridging fault simulator[END_REF][START_REF] Lu | Laser-Induced Fault Simulation[END_REF][START_REF] Meyer | Active timing multilevel fault-simulation with switch-level accuracy[END_REF][START_REF] Santos | Defect-oriented mixed-level fault simulation of digital systems-on-a-chip using HDL[END_REF], the most recent one is [START_REF] Lu | Layout-aware laser fault injection simulation and modeling: From physical level to gate level[END_REF], which is based on the open-source Lifting [START_REF] Bosio | LIFTING: A Flexible Open-Source Fault Simulator[END_REF]. The major issue with these fault simulators is that they rely on electrical models [START_REF] Douin | Electrical modeling for laser testing with different pulse durations[END_REF][START_REF] Godlewski | Electrical modeling of the effect of beam profile for pulsed laser fault injection[END_REF][START_REF] Sarafianos | Building the electrical model of the pulsed photoelectric laser stimulation of an NMOS transistor in 90nm technology[END_REF] that are technology dependent. For instance, in [START_REF] Hériveaux | Electrical Modeling of the Effect of Photoelectric Laser Fault Injection on Bulk CMOS Design[END_REF], the authors proposed a model that includes the vertical parasitic bipolar junctions inherent to mosfets in the fault injection process that may lead to IR drop effects. However, they did not extended their work beyond the scope of a single inverter. In fact, dimensioning the rc network of power/ground rails is a difficult task, since the rc values depend on the technology, the size of cells, the position of voltage taps on the rails, the rc parasitics, etc.

The issue being that, as far as we know, there is no tool capable to simulate laser-induced IR drop and its propagation in a large circuit. Thus, the first and main objective of this work is to introduce the devised methodology to simulate at the electrical level the effect of IR drop on the fault injection sensitivity using standard cad tools; the second objective is to illustrate, on simulation grounds, that laser-induced IR drop has to be considered since it may result in underestimating the risk of fault injection.

STATE OF THE ART 2.1 Modeling laser effects on ICs

2.1.1 Laser Induced Transient Currents. ics are known to be sensitive to induced transient currents. These currents may be caused by laser shots passing through the device, creating electronhole pairs along the path of the laser beam [START_REF] Johnston | Charge generation and collection in p-n junctions excited with pulsed infrared lasers[END_REF]. The induced charge carriers recombine without any significant effect, unless they reach the strong electric field found in the vicinity of reverse biased PN junctions. In this case, the electrical field puts these charges into motion and a transient current appears as well as a transient fault. The nature of this fault is similar to the ionization effect generated by energetic particles [START_REF] Habing | The Use of Lasers to Simulate Radiation-Induced Transients in Semiconductor Devices and Circuits[END_REF].

As an example (the cross section of an inverter), Fig. 1 illustrates where laser shots may generate parasitic currents. In case the inverter input is in low state ('0') the most laser-sensitive part of the inverter is the drain of the nmos transistor since there is a reverse biased PN junction between the drain and the P subst r at e . Thus, an induced transient current (I P h ) flows from the drain of the nmos to the P subst r at e biasing contact. A similar reasoning can be made when the inverter input is high ('1'). In that case, the susceptible part of the inverter is the drain of the pmos transistor. In case of Fig. 1, a part of the induced photocurrent (I P h ) charges the inverter output capacitance. As a result the inverter output undergoes a voltage transient.

Another transient current component flowing from v DD to g ND that may have a significant effect on the fault injection mechanism is taken into consideration by the model of Fig. 1 [START_REF] Viera | Role of Laser-Induced IR Drops in the Occurrence of Faults: Assessment and Simulation[END_REF]. This transient current is induced in the reversed biased Psub-Nwell junction that surrounds every Nwell. Even if the laser beam is directed towards a sensitive nmos, the laser beam also induces charge carriers that will be sufficiently close to a Psub-Nwell junction to induce a transient current in it flowing from v DD to g ND .

The Psub-Nwell junction is always reversed biased and has a larger area than that of a transistor drain (the parameter S in (1)). Thus, it is no surprising that the authors of [START_REF] Dutertre | Improving the ability of Bulk Built-In Current Sensors to detect Single Event Effects by using triple-well CMOS[END_REF] reported on experimental basis that the transient current component flowing directly from v DD to g ND (IP P sub_nwell in Fig. 1) may be more than an order of magnitude greater than those flowing in the drains of the sensitive transistors (I P h in Fig. 1). This transient v DD to g ND current may thus have a significant influence on the laser fault injection mechanism since it will produce a temporary supply voltage drop (IR drop) [START_REF] Dutertre | Improving the ability of Bulk Built-In Current Sensors to detect Single Event Effects by using triple-well CMOS[END_REF][START_REF] Hériveaux | Electrical Modeling of the Effect of Photoelectric Laser Fault Injection on Bulk CMOS Design[END_REF][START_REF] Viera | Role of Laser-Induced IR Drops in the Occurrence of Faults: Assessment and Simulation[END_REF].

2.1.2 Spatial Distribution of Laser Beam Energy. The beam diameter is one of the most important propagation attribute of a laser beam in a class of commonly measured parameters (beam diameter, spatial intensity distribution, beam quality factor etc.). A commonly used definition of the laser beam diameter is derived from the bivariate normal distribution of its intensity leading to measure the beam diameter at 86.5% of its maximum value [START_REF] Buchner | Pulsed-Laser Testing for Single-Event Effects Investigations[END_REF], or a drop of 1 e 2 from its peak value. The effects of a Near Infrared laser beam have been modeled in [START_REF] Messenger | Collection of Charge on Junction Nodes from Ion Tracks[END_REF] and later in [START_REF] Sarafianos | Building the electrical model of the pulsed photoelectric laser stimulation of an NMOS transistor in 90nm technology[END_REF]. In the latter work, it is shown that the induced photocurrent, which is spatially distributed as a bivariate normal distribution, has a peak amplitude I ph_peak that follows the empirical equation:

I P h_peak = (a × V + b) × α дauss (x,y ) × Pulse w × S ( 1
)
where V is the reverse-biased voltage of the exposed PN junction, a and b are constants that depend on the laser power, α дauss (x,y ) is a term related to the bivariate distribution of the laser beam amplitude in space, Pulse w is a term allowing to take into account the laser pulse duration and S is the area of the PN junction. One can refer to [START_REF] Sarafianos | Building the electrical model of the pulsed photoelectric laser stimulation of an NMOS transistor in 90nm technology[END_REF] for additional details of the above parameters. By way of illustration, Fig. 2 shows a three-dimensional view of the normalized amplitude of a laser spot. Beam intensity at a given coordinate (x,y) represents the amount of power delivered by the laser source at this specific point. 2.1.3 Electrical Model of a Cell Under Laser Illumination. Fig. 3a introduces, in case of an inverter, the classical model showing that the effect of a laser is modeled by a current source placed between the drain and the source of the laser-sensitive transistor (pmos transistor in this example). Fig. 3b shows, in case of an inverter, the enhanced electrical model taking into account the laserinduced IP P sub_nwell current. Without the power-grid model (i.e., considering v DD and g ND ideals), it would be impossible to take the current IP P sub_nwell into account. Consequently, the laserinduced IR drop contribution also would not be taken into account during simulations. This work proposes in its flow the use of an Electromigration/IR drop (EMIR) cad tool to automatically provide the power-grid model for each cell in the circuit.

The current sources in Fig. 3a and Fig. 3b have a profile of a double exponential, such as the one illustrated in Fig. 4. The currents have a peak amplitude defined by [START_REF] Bosio | LIFTING: A Flexible Open-Source Fault Simulator[END_REF]. Since the parameter S (area of the PN junction) corresponds to the cell's Nwell area, thus, the current component IP P sub_nwell is larger than that induced at a sensitive transistor drain (IPh) since the drain area is smaller than the Nwell's area (see [START_REF] Dutertre | Improving the ability of Bulk Built-In Current Sensors to detect Single Event Effects by using triple-well CMOS[END_REF] for an experimental assessment). The IP P sub_nwell current source is attached to the biasing contacts of the Nwell and the P subst r at e (for standard cells without embedded biasing contacts, the current source is connected to the closest). The various IP P sub_nwell currents add up and flow from v DD to g ND trough the power and ground networks of the device under attack. Because the power grid exhibits both resistive and capacitive electrical behaviors, a local voltage drop and ground bounce occurs thus reducing the voltage swing seen by standard cells in the close vicinity of the laser spot. Considering the above, this paper provides a method based on standard cad tools to take at chip level the effect of laser-induced IR drops into account.

C Load '1' '0' >> '1' IPh (a) C Load '1' '0' >> '1' IPpsub_nwel IPh Power-grid Model Power-grid Model (b)

Previous Works on Laser Fault Simulation

Laser fault injection may be anticipated or studied by using simulation tools at different abstraction levels: physical, electrical or logical. In this section, previous works that proposed laser fault simulation tools are reviewed in order to justify the need for the methodology presented in this work.

Physical Level. Based on Technology Computer Aided

Design (TCAD), the authors in [START_REF] Llido | Characterization and TCAD simulation of 90 nm technology transistors under continous photoelectric laser stimulation for failure analysis improvement[END_REF] characterize and analyze photoelectric effects induced by static 1064 nm wavelength laser on a 90 nm technology nmos transistor. In [START_REF] Fleetwood | Using TCAD Modeling to Compare Heavy-Ion and Laser-Induced Single Event Transients in SiGe HBTs[END_REF], Silicon-Germanium Heterojunction Bipolar Transistor (SiGe HBT) models are used in TCAD to investigate single event transients induced by heavy-ion broadbeam and pulsed-laser sources. Although TCAD is the ultimate tool to simulate laser effects on ics, this simulator is extremely CPU consuming and can only be applied to individual transistors or small circuit areas.

Logic

Level. The authors of [START_REF] Papadimitriou | A multiple fault injection methodology based on cone partitioning towards RTL modeling of laser attacks[END_REF] proposed a methodology for multiple fault injection at the Register Transfer Level (rtl). The methodology would reduce the fault space of laser fault injection campaigns by using the locality characteristic of laser fault, and through a partitioning of the rtl description of the circuit. Their efforts involve the development of an rtl fault injection approach more representative of laser attacks than random multi bits fault injection. Unfortunately, as a rtl fault simulator, the fault model is defined as a logic pulse with different widths, which is not sufficient to take into account neither the laser parameters nor IR drop effects.

Electrical Level.

Laser fault simulation at the electrical level is a good tradeoff between speed (logic level) and accuracy (physical level). Therefore, it is possible to represent the laser physical phenomenon in the scope of a whole system. Although the simulation time might be an issue, today's electrical simulators are up to 100x faster than baseline spice simulators without loss of accuracy. Furthermore when large circuits are simulated, it is possible to profit by the use of hybrid simulation in which only the affected zone of the IC is simulated with spice accuracy while the non affected cells are simulated with gate level accuracy.

To the extent of our knowledge, the most recent fault simulator at the electrical level was proposed by [START_REF] Lu | Layout-aware laser fault injection simulation and modeling: From physical level to gate level[END_REF]. Their simulator is based on the open-source Lifting [START_REF] Bosio | LIFTING: A Flexible Open-Source Fault Simulator[END_REF], which allows both 0-delay and delayannotated simulations of digital circuits using layout information to derive the laser spot location. They also use multi-level simulation, trading of speed for accuracy. The major issue with these fault simulators is that they to rely on electrical models [START_REF] Douin | Electrical modeling for laser testing with different pulse durations[END_REF][START_REF] Godlewski | Electrical modeling of the effect of beam profile for pulsed laser fault injection[END_REF][START_REF] Sarafianos | Building the electrical model of the pulsed photoelectric laser stimulation of an NMOS transistor in 90nm technology[END_REF] that are technology dependent. Even though it is possible to dimension these models, it is hard to obtain accurate results when dealing with new technologies.

For instance, the contribution of IR drop effects play a significant role in the fault injection process as reported in [START_REF] Viera | Role of Laser-Induced IR Drops in the Occurrence of Faults: Assessment and Simulation[END_REF]. The authors of [START_REF] Hériveaux | Electrical Modeling of the Effect of Photoelectric Laser Fault Injection on Bulk CMOS Design[END_REF] modeled a rc network in the power/ground rails to demonstrated the significant contribution of the current induced by vertical parasitic bipolar junctions inherent to mosfets in the fault injection process. However, they did not study the effect of the IR drop induced by laser shots, i.e., its impact in the fault injection mechanism. They also did not extended their work beyond the scope of a single inverter since they manually dimensioned the values of the rc components, which would be a difficult task to do for a whole circuit.

2.2.4 Summary. What has been observed so far is that there is a great improvement of laser fault models. However the models were developed at the level of a single gate, ignoring thus the effects of laser-induced IR drops at chip level. Regarding laser fault simulators, they usually use the simple fault model in which current sources are attached to the drain and bulk of laser sensitive transistors [START_REF] Jordan | Photoeffect on diffused P-N junctions with integral field gradients[END_REF][START_REF] Wirth | The Transient Response of Transistors and Diodes to Ionizing Radiation[END_REF]. Unfortunately, this fault model was created at a time when laser sources with 1 µm to 5 µm spot diameter were used to target only one sensitive PN junction at the same time. For advanced technologies this model is questionable. For a 28 nm technology, the standard cells have a height value of about 1.2µm, meaning that even lasers with 1 µm spot diameter will also illuminate the Psub-Nwell junction (see Fig. 1) and thus induce significant IR drop in the area surrounding the laser spot.

In order to use a fault model that takes into account the IR drop contribution induced by the current component created between the Psub-Nwell junction, it is necessary to model by a rc network the power/ground rails. Modeling the rc network of a large circuit is not a task to be performed manually. In view of this limitation, i.e., that current laser fault simulators do not use complete and accurate fault models, we propose a fault simulation methodology that uses an EMIR cad tool to automatically provide the rc network of the power/ground rails for a given design. It also provides the transient voltage that propagates along the power rails as a result of the IP P sub_nwell current. The methodology can be used for any circuit designed in any technology supported by the standard cad tools. Next section presents in details the proposed methodology.

PROPOSED METHODOLOGY FOR LASER FAULT SIMULATION

The diagram presented in Fig. 5 proposes a step by step simulation methodology that makes it possible to simulate laser fault injection in large scale circuits. This methodology, which is based on standard CAD tools (Cadence Voltus T M [START_REF]Voltus IC Power Integrity Solution[END_REF] for EMIR simulation and Cadence Spectre XPS [START_REF]Spectre eXtensive Partitioning Simulator[END_REF] for the electrical and hybrid simulation), allows to analyze the impact of laser shots on complex circuits by drawing laser-induced fault sensitivity maps. The methodology can be easily adapted to provide other set of results besides the ones reported in this work. As far as we know, this is the first methodology able to simulate laser effects on ics that takes into account laser-induced IR drop effects. Although Cadence tools were used, any other tools that are able to perform IR drop analysis and spice like simulations can be used. Fig. 5 is subdivided in numbers that represent each step described in the following sections.

Step 1: defining simulation parameters

In the first step, a shell script file (main.scr) defines parameters characterizing the laser shot. Among them, one can find: the laser beam diameter, the duration of the laser shot, the time at which begins the laser shot with regard to the operation of the ic, the (X,Y) displacement step of the laser spot when one aims to draw fault sensitivity maps (detailed in step 5), etc. This file is also responsible for calling the necessary tools and scripts for the correct execution of the simulation flow.

Step 2: data preparation for the EMIR CAD tool

Most of the inputs that are inside the dashed rectangle "EMIR CAD Tool" of Fig. 5 are files that were automatically generated by the design CAD tool (Cadence Innovus [START_REF]Innovus Implementation System[END_REF]). Other files were obtained from the design kit of the technology. It is out of scope of this work to explain each of these files in detail. It suffices to say that they are necessary to model the rc network in the power/ground rails and perform IR drop analysis in Cadence Voltus T M , both necessary for the accomplishment of the proposed methodology.

Step 3: spatial location of the laser spot

In this step it is necessary to know the dimension of the design and the number of simulated laser shots that are going to be applied over the circuit. For this work, an ARM 7 with a 110 µm × 70 µm area was used (more details are provided in Section 4). If a displacement step of x, y: 5µm is set, then, in order to sweep the whole circuit, beginning at x, y: (0, 0) and ending at x, y: (110, 70), it would demand 345 laser shots as illustrate in Fig. 6. This step allows to know where the laser spot illuminates the ic during each simulation. Next step shows which cells are illuminated by the laser spot for each x, y position and at which intensity.

X axis (µm) Y axis (µm)

Step 4: laser-induced fault injection

Faults induced by laser illumination can be simulated by specifying the current amplitude of the current sources that compose the laserinduced transient fault model (Fig. 3) of each standard cell in the circuit. Therefore, to simulate a circuit being attacked by means of laser fault injection it is necessary to know which cells will be affected by the laser.

Several ways can be adopted in order to discover the values to be assigned to the current sources in the fault model (Fig. 3). This methodology benefits from a feature present in Cadence Voltus T M . This tool allows to apply an amount of current to a defined region, in this way, several small rectangular regions are defined and the current amplitude of that region follows the spatial distribution of the laser-induced photocurrent defined by [START_REF] Bosio | LIFTING: A Flexible Open-Source Fault Simulator[END_REF]. Fig. 7 illustrates how the rectangular regions can be used in order to apply the laser power (current induced by the laser) to each rectangle. The following code example represents the characterization of a rectangle (current region) located at the center of the laser spot (Fig. 7). Therefore its I P h_peak is maximum, 100% or 1 mA for this example. The double exponential has a step size of 5 ps, the peak is thus found at it apex, i.e., 1.510ns, considering 10 ps of rise time and fault starting at 1.500 ns. Other parameters such as capacitances are extracted from .lib and .spi files of the technology for each affected cell. The resolution of each rectangle is 250 nm as shown by the last parameter: -region "x1 y1 x2 y2". The dimension of the rectangle can be changed according to the precision needed to model the laser spot.

μm

c r e a t e _ c u r r e n t _ r e g i o nc u r r e n t { 1 . 5 0 0 ns 0 . 0 0 0mA 1 . 5 0 5 ns 0 . 8 2 0mA 1 . 5 1 0 ns 1 . 0 0 0mA 1 . 5 1 5 ns 0 . 9 5 0mA . . . 1 . 8 0 0 ns 0 . 0 0 0mA}l a y e r M2i n t r i n s i c _ c a p C l o a d i n g _ c a p Cr e g i o n " 1 . 5 0 1 . 5 0 1 . 7 5 1 . 7 5 "

Step 5: mapping the circuit

In this step, Cadence Voltus T M is used with the purpose to perform laser-induced IR drop simulations for the different laser spot locations calculated during step 3. All other simulation parameters being kept constant (spot diameter, intensity, etc.).

Clarifying, IR drop can be defined as the power supply noise induced by currents flowing through the resistive parasitic elements of the power distribution network. In this work, the laser-induced IR drop is also considered, meaning that the laser-induced current will accumulate with the dynamic current of a cell, thus increasing its IR drop while the laser is active (IP P sub_nwell ∅).

For each iteration of this step, a table containing the evolution in time of each cell's voltage swing amplitude (v DD -g ND ) is saved for future analysis since different cells are affected by the laser shot. It is also possible to save a table with the dynamic current in time, which translates directly to the amplitude of the current IP P sub_nwell for each cell in the circuit. Table 1 illustrates for three different cells the remaining voltage swing when the laser effect reaches its apex (peak of the double exponential transient current from Fig. 4). Note in this example that, for the laser spot position 1 (cf. Table 1) the cells are more affected (lower voltage swing) as the epicenter of the laser spot is closer to these three cells. For laser spot positions 2 and 10, the cells are less affected since the epicenter of the laser spot is far away from the cells listed in the table.

Step 6: inserting IP Psub_nwell

The IP P sub_nwell current component induces voltage drops in the power/ground rails. This effect is captured thanks to Cadence Voltus T M in the previous steps. In this step a shell script replaces the ideal v DD and g ND in the original spice netlist by waveforms saved in step 5 for each cell in the circuit.

Step 7: inserting I ph

A shell script is used in order to add a current source between the drain and bulk of pmos and nmos transistors. It models the laserinduced currents that may turn into faults. Note that only one of these current sources are activated depending on which drain's PN junction is reversely polarized. For this, it was necessary to run a fault free electrical simulation and save a golden table with all inputs and outputs of each cell as a function of time.

Knowing that the IP P sub_nwell current is defined as a f actor × I ph because of the parameter S in (1), it is possible to compute the f actor value to be applied to each cell by analyzing the .lef file of each cell and to estimate the area of the affected PN junction of the transistor's drain as well as the Nwell of the same cell.

Step 8: electrical/hybrid fault simulation

At this point, electrical simulations are performed for each laser shot with different locations as defined on step 3. Electrical simulations are time consuming depending mainly on the circuit's size and available computing resources. To circumvent this drawback, a hybrid simulation has to be performed. This simulation defines a region of the circuit where only the most affected cells are simulated with spice accuracy. For the hybrid simulation, Cadence Spectre XPS simulator is used. To define the cells that are going to be simulated at logical level, a threshold voltage is defined based on the v DD -g ND (ir drop + ground bounce) values provided by Table 2. If a cell's power/ground voltage is close to the nominal v DD and g ND , it is considered that this cell is not affected by the laser shot, since it is far away from the epicenter of the laser spot. For example, if a threshold voltage of 10% of the nominal v DD = 1 V is defined, then all cells with ir drop + ground bounce lower than 100 mV are simulated at the logic abstraction level.

Table 2 shows the number of cells simulated with the logic abstraction level for different threshold voltages and different spot locations. The spot locations were selected by chance with the purpose to show that the number of affected cells changed depending on the location where the laser shot was applied.

LASER FAULT SIMULATION RESULTS

Testbench

In order to simulate the effects of laser-induced faults on complex systems, simulations were performed for different circuits, however only results for an ARM 7 processor (dut) are shown in details. All circuits were synthesized using a 28 nm technology. The core nominal voltage of the dut is 1 V and the clock period is 1 ns. The dut has an area equal to 110 µm × 70 µm.

Circuit Inventory.

The evaluated design is composed by 5.21 k cells, 5.34 k nets and 90 k nodes. The power-grid model generated by Cadence Voltus T M has 100 k resistors and 90 k capacitors.

4.1.2 Laser Spot Diameter. Typical laser sources used to produce faults are characterized by a beam diameter equal to 1 µm, 5 µm or 20 µm and a wavelength of 1064 nm. Although the minimum diameter of a laser spot is 1 µm, given the laws of optic its effect area extends far beyond [START_REF] Darracq | Backside SEU laser testing for commercial off-the-shelf SRAMs[END_REF][START_REF] Roscian | Fault Model Analysis of Laser-Induced Faults in SRAM Memory Cells[END_REF]. Consequently, a laser spot does not induce a single transient current in a single cell, but several transient currents at different sensitive nodes of the target. Without loss of generality, a spot diameter of 1 µm has been chosen for the experiments reported below.

Simulation Performance

The performance of the simulation depends directly on the available computing resources and the complexity of the simulated circuit. The available processor used to perform simulations was an Intel Xeon E5630 @ 2.53 GHz with two cores and 16 GB of RAM. Since the proposed method deals with the simulation of laser-induced fault injection, other factors should be also taken into account, such as the laser spot diameter, its power and the duration of the laser shot. Considering the simulation performed using only Spectre accuracy, the simulation takes more time to be performed when comparing to the simulation of the same circuit in a fault free scenario. This happens as the cells no longer have ideal v DD and g ND , thus the simulator has to decrease the simulation step to account with laserinduced transient currents, which are in the ps range. Therefore, since the diameter of the laser spot determines how many cells are affected, it influences on the time required by the simulator to perform necessary calculations. When using hybrid simulation, it is possible to decrease the amount of cells simulated with Spectre accuracy, thus reducing simulation time.

Table 3 shows simulation times for different circuits using Spectre XPS (hybrid simulation). Simulation times for other simulators (Spectre accuracy only) are not shown as they take at least 22 times more to simulate. Simulations were also performed using Spectre and Spectre APS with the intention to compare results regarding the accuracy of Spectre XPS. In all cases the results were the same, i.e., the same sensitivity maps presented in the next section were obtained. In fact, for this kind of analysis there is no need to have the same precision as simulations for RF designs, in which the Spectre RF simulator is recommended.

Laser Propagation on the Circuit Surface

To illustrate how the IR drop propagates in the circuit, refer to Fig. 8a and Fig. 8b. In Fig. 8a, for which no laser effect is considered, the IR drop across the rails reach the maximum of 50 mV. In this figure, the voltage drop is uniquely due to normal switching activity. Even though not fully uniform, the IR drop affects almost the whole circuit. Fig. 8b (obtained in step 5 of the proposed method) illustrates how the laser effect propagates in the circuit. In presence of a single laser shot with a spot diameter of 1 µm at coordinates x=60 µm, y=35 µm, the effect area extends along the X axis of the powergrid main metal lines for more than 100 µm (the effect area has a shape that is stretched horizontally along the power supply rails as they provide a propagation path to the laser-induced IR drop and ground bounce). Whereas its extension along the Y axis is only approximately 3 µm, i.e., three times the laser spot diameter. The peak value of the induced voltage transient in the power lines is 400 mV (Fig. 8b). At this time, the voltage swing is reduced to 600 mV. This value is far below the nominal core voltage of 1 V. Thus laser-induced IR drop may induce faults in the circuit, such as timing errors or even data disruption. There are hundreds of standard cells inside the area affected by the laser when considering a 28 nm technology, meaning that the cells inside the affected area will absorb the laser-induced current according to the surface distribution of the laser beam given by 1.

Simulated Scenarios

We report a total of 4 simulated scenarios among the ones studied. They are illustrated in Fig. 9 that shows in the first line the clock signal waveform used as a time reference. The two other lines give the typical evolutions observed during simulations, of the signal Q x , the output of the cell 'x' of the design under illumination, in two different situations. These two situations represent the behavior when a laser pulse with 250 ps duration starts at 1.5 ns and 1.7 ns respectively. These times are thus closer and closer to the next rising clock edge that occurs at 2 ns. The second line of Fig. 9 gives these evolutions when only the I Ph current sources with a double exponential shape are considered to model laser effects. In the third line, the curve has a smoother double exponential waveform when comparing with the profile of double exponential current pulse (c.f. second line) proposed by [START_REF] Messenger | Collection of Charge on Junction Nodes from Ion Tracks[END_REF] due to the filtering effect (RC effect) of the supply voltage network. In fact, the profile proposed by [START_REF] Messenger | Collection of Charge on Junction Nodes from Ion Tracks[END_REF] aims to model alpha-particle hits, which does not exactly correspond to charge generation and collection in PN junctions excited with pulsed infrared lasers as analyzed in [START_REF] Johnston | Charge generation and collection in p-n junctions excited with pulsed infrared lasers[END_REF]. report fault maps for which I P h , IP P sub_nwell and the power-grid model are considered (scenarios 3 and 4). By comparison to the first line, it reveals that the fault areas are larger than expected for the considered laser shot times. It also unveiled an extension of the laser sensitivity in time, in which the number of faults are increased respectively by a factor of 2.6 and 3.1 for the laser applied at 1.5 ns and 1.7 ns. This demonstrates that IR drops induced by laser shots play an important role in the occurrence of faults. Not taking this effect into account leads to over optimistic results regarding the threshold of fault injection.

Time (ns)

CONCLUSIONS

This paper presented a new method that allows to simulate laserinduced faults at the electrical level in large-scale circuits by using standard CAD tools. Its main intent is to take into account the IR drop effects induced by laser shots: a key parameter in the fault injection process. For each cell in the circuit, a high accuracy electrical fault model that includes the voltage drop effects in the power and ground rails was used thanks to the use of an EMIR CAD tool. The method was applied to a test-chip in order to demonstrate how fault sensitivity maps can be drawn with the purpose of assessing the contribution of laser-induced faults into the circuit.

Figure 1 :

 1 Figure 1: Laser-induced current components. Cross-section of a cmos inverter.

Figure 2 :

 2 Figure 2: Three-dimensional view of a laser beam in terms of intensity per area. 100% of laser beam intensity represents the epicenter of the laser spot.

Figure 3 :Figure 4 :

 34 Figure 3: Laser-induced transient fault model applied to an inverter with its input biased at v DD . (a) Classical model. (b)Improved model including the IR drop and ground bounce contribution induced by IP P sub_nwell for a given power-grid model[START_REF] Viera | Role of Laser-Induced IR Drops in the Occurrence of Faults: Assessment and Simulation[END_REF].

7 8Figure 5 :

 75 Figure 5: Procedure used to draw laser-induced fault sensitivity maps using the proposed methodology.

Figure 6 :

 6 Figure 6: Spatial location of the laser spots. Each point corresponds to a laser shot at different positions (each point corresponds to a simulation).

Figure 7 :

 7 Figure 7: Laser-induced current regions applied over standard cells of a CMOS 28 nm technology. The current amplitude of each region is defined by (1).

Table 1 :

 1 List of cells of the circuit with their voltage swing at the apex of the laser shot.

Figure 8 :

 8 Figure 8: ARM 7 layout with 5k+ instances: (a) Maximum voltage drop (IR-Drop + ground bounce) in normal operation condition. (b) Maximum voltage drop in presence of a laser shot with spot diameter equal 1µm.

Figure 9 :

 9 Figure 9: waveforms observed during simulations at the output of gates illuminated by a laser beam. Line 1: clock signal. Line 2: waveforms observed when considering I P h contribution only. Line 3: waveforms observed when considering I P h + IP P sub_nwell contributions.

4. 5 . 2

 52 Contribution of I ph and IP Psub_nwell . Fig.10cand Fig.10d

Figure 10 :

 10 Figure 10: Maps of laser-induced faults for the simulated scenarios: (a-b) laser applied at 1.5 ns and 1.7 ns respectively, considering I Ph contribution only. (c-d) laser applied at 1.5 ns and 1.7 ns respectively, considering I P h + IP P sub_nwell contributions.

Table 2 :

 2 Number of cells simulated with the logic abstraction level for different threshold voltages and different spot locations. (5.21k cells in the circuit.)

	Threshold	No. of cells No. of cells
	(ir drop + bounce) (spot loc. 1) (spot loc. 2)
	10%	2535	2625
	15%	4510	4585
	20%	4641	4620

Table 3 :

 3 Simulation performance for different circuits regarding one point of the fault sensitivity map (1 simulation out of 345 simulations to create a complete map).

	Circuit	No. of cells Simulation time
	Arm7	5.210	1min 02s
	S38584 (ISCAS'89)	20.705	1min 20s
	B18 (ITC'99)	52.601	3min 05s
	B19 (ITC'99)	105.344	6min 35s


Wisama Khalil
email: wisama.khalil@ircyn.ec-nantes.fr

Sébastien Besnard

Philippe Lemoine

COMPARISON STUDY OF THE GEOMETRIC PARAMETER CALIBRATION METHODS

Keywords: Calibration, robots, static accuracy, identification, geometric parameters, identifiable parameters

come

Introduction

The absolute accuracy of a robot depends to a large extent on the accuracy of the numerical values of the geometric parameters used in the direct and inverse kinematic models. Geometric calibration of robots is the process by which the geometric parameters are precisely identified. In general, the geometric calibration is carried out by solving a system of linear or nonlinear equations, which are a function of the geometric parameters, the joint positions, and the location (position and orientation) of the end-effector frame.

A large variety of methods have been proposed for the calibration of the geometric parameters. They vary by the type of endpoint sensing, the endpoint constraints, and by the number of equations. In this paper the following methods have been used in the comparison:

• Classical Open Loop Method: this method is based on using a set of configurations for which the joint positions and the corresponding Cartesian coordinates of the terminal frame position or location are given [1 4]. • Relative Location Measure: the calibration is carried out using the joint positions and the relative Cartesian location between each couple of configurations.

• Distance Measure: this method uses the joint positions of a set of configurations and the distance between the terminal points of each couple of configurations. • Frame (or position) Link: this autonomous method needs only the joint positions: no external sensor is needed. It can be applied for redundant robots or nonredundant robots, which can achieve the same location (or position in the case of position link) of the terminal frame by multiple configurations. The data used in the calibration correspond to the joint variables of some sets of robot configurations, giving the same Cartesian location (or position) of the terminal link of the robot [5 8]. • Plane Link: in this case, the calibration is carried out using the joint positions of a set of configuration whose terminal points are in the same plane [9 13]. This paper focuses primarily on the identifiable parameters of these methods and on their convergence rate. The geometric parameters are estimated using a linearized model, which is solved iteratively using a least squares criterion, and by updating the identified parameters and the observation matrix after each iteration. Simulations show that even in the case of large errors in the parameters, such techniques work very well.

The paper is organized as follows: the parameters defining the robot are presented in Section 2. Section 3 describes the calibration methods. Section 4 presents the simulation of the calibration methods on the Puma and Stanford robots as well as the comparison of the different methods. Section 5 is the conclusion.

Description of the Geometric Parameters

The geometric parameters are the constant parameters and the offsets of the joint positions that are used to calculate the location of the end-effector of the robot with respect to the fixed world frame.

We consider serial robots consisting of n joints and n + 1 links. Link 0 is the base and link n is the terminal link, frame j is defined fixed on link j, we denote:

• frame -1: the fixed reference frame • frame n + 1: the end-effector frame The end-effector location can be computed with respect to the reference frame by the direct geometric model:

-1 T n+1 = -1 T 0 . 0 T n (q). n T n+1 (1)
where i T j is the 4×4 transformation matrix defining frame j relative to frame i and q is the vector of joint variables.

The definition of the link frames used here follows the notations of the modified Denavit and Hartenberg method proposed by Khalil and Kleinfinger [START_REF] Khalil | A new geometric notation for open and closed loop robots[END_REF][START_REF] Khalil | Mod elisation identi cation et commande des robots[END_REF]. Frame j is defined such that z j is along the axis of joint j, and x j is perpendicular to z j and z j+1 . Frame j is defined relative to frame j -1 by the matrix j-1 T j , which is a function of the four parameters (α j , d j , θ j and r j ) as shown in Fig. 1. The joint variable q j is equal to θ j if joint j is revolute and r j if joint j is translational. Although frames -1 T 0 and n T n+1 can be arbitrarily defined, it has been shown that the calculation of -1 T n+1 can be obtained as a function of the geometric parameters of n + 2 frames represented by four parameters for each frame (α j , d j , θ j , and r j ) for j = 0, . . . , n + 1, except for the first frame for which α 0 = 0, d 0 = 0 [START_REF] Khalil | Mod elisation identi cation et commande des robots[END_REF][START_REF] Khalil | Identifiable parameters and optimum configurations for robots calibrations[END_REF].

In the calibration process, we have to identify the deviation of the real parameters from the nominal values, thus to identify: ∆θ 0 , ∆r 0 and ∆α j , ∆d j , ∆θ j , ∆r j (for j = 1, . . . , n + 1).

If the axis of joint j is parallel to the axis of joint j -1, an additional parameter ∆β j , which is called the Hayati parameter, must be considered [START_REF] Hayati | Robot arm geometric link calibration[END_REF]. The nominal value of β j is equal to zero such that: j-1 T j = Rot(y, β j ).Trans(x, d j ).Rot(x, α-j).Trans(z, r j ).Rot(z, θ j )

(2) Using the same procedure, one can also identify the errors on the joint angle gains giving the joint positions from the position sensor readings. They will be denoted by ∆K j (for j = 1, . . . , n) [START_REF] Hollerbach | A review of kinematic calibration[END_REF].

The geometric parameters of the Puma manipulator are given in Table 1 and those of the Stanford arm are listed in Table 2. We have supposed the same -1 T 0 and n T n+1 for both robots. Beside the parameters β j , α j , d j , θ j , and r j we see the following parameters:

• σ j : which defines the type of the frame

σ j = 0 if joint j is revolute σ j = 1 if joint j is prismatic σ j = 2 if frame j is fixed
• K j : the joint angle gain of joint j , for j = 1, . . . , n Table 1 The Geometric Parameters of the Puma (6R) Robot Units are given in meters for the distances and radians for the angles.

j σ j α j d j θ j r j β j K j 0 2 0 0 π/2 0.5 0 1 0 0.1 0 θ 1 0 0 2 0 -π/
j σ j α j d j θ j r j β j K j 0 2 0 0 π/2 0.5 0 0 1 0 0.1 0 θ 1 0 0 1 2 0 -π/2 0 θ 2 0.
1 r 6 is not equal to zero because of the supposed terminal frame.

Calibration Methods

The calibration methods studied in this paper differ according to the variables used in the calibration model and the type of measuring device. A unified approach can be formulated using the framework suggested in [START_REF] Wampler | An implicit loop method for kinematic calibration and its application to closedchain mechanisms[END_REF]. According to this formulation, the kinematic calibration equation can be written in the general form:

0 = f (q, x, η r ) (3)
with • x as the Cartesian variables giving the position and orientation of the terminal frame • q as the (n × 1) joint variables vector • η r as the (N p × 1) vector of the real (unknown) values of the geometric parameters 2

Each method can be classified by a calibration index giving the number of equations of the function f [START_REF] Hollerbach | The calibration index and taxonomy of kinematic calibration methods[END_REF].

The nonlinear calibration equation ( 3) can be linearized to get the differential equation: ∆y(q, x, η) = Ψ(q, η).∆η

with

• ∆η = η r -η, defining the vector of the errors on the geometric parameters • η as the vector of the nominal values of the geometric parameters To estimate ∆η, (4) will be applied using a sufficient number of configurations Q = q 1 , . . . q m to obtain an overconstrained system of equations:

∆Y (Q, X, η) = W(Q, η).∆η (5)
with

∆Y =      ∆y 1 (q 1 , x 1 , η) . . . ∆y m (q m , x m , η)      , W =      Ψ(q 1 , η) . . . Ψ(q m , η)     
where W is the (r × N p) observation matrix with r N p, X = x 1 , . . . , x m .

It can be seen that if some columns of W are dependent, and W is of rank b < N p, then relation ( 5) can be reduced to [START_REF] Khalil | Calculation of the identifiable parameters for robots calibration[END_REF]:

∆Y (Q, X, η) = W b (Q, η).∆η b (6)
where W b contains arbitrarily b independent columns of W, the corresponding parameters are called the identifiable (or base) parameters. They will be denoted by the vector ∆η b .

The determination of the identifiable parameters ∆η b must be done before the identification process. They can be obtained numerically using the QR decomposition of a matrix W similar to that defined in ( 5), but obtained using random configurations satisfying the constraints of the calibration method. The outline of this method is given in Appendix A.

Equation ( 6) will be solved to get the least squares error solution to the current parameter estimate. This procedure is iterated until the vector ∆η b is sufficiently small. After each iteration, the geometric parameters will be updated in W b and ∆Y . In the following we will use the standard pseudo-inverse to estimate ∆η b . Other methods using statistical solutions are given in [START_REF] Wampler | An implicit loop method for kinematic calibration and its application to closedchain mechanisms[END_REF]. The condition number of W b gives an index measure for the observability of the parameters in the calibration system of equations. To get good results the condition number must be close to one [START_REF] Khalil | Identifiable parameters and optimum configurations for robots calibrations[END_REF]. Other indices of the observability using the product of the singular values can be also used [START_REF] Borm | Determination of optimal measurement configurations for robot calibration based on observability measure[END_REF]. It has been shown that the condition number is more sensitive [START_REF] Hollerbach | Closed-loop kinematic calibration of the RSI 6-DOF hand controller[END_REF].

The elements of ∆η b have different units: meters (for distances) radians (for angles) or even no units (gain coefficients). The effect of this heterogeneity can be reduced by a column normalization of the matrix W b [START_REF] Hollerbach | The calibration index and taxonomy of kinematic calibration methods[END_REF][START_REF] Schr 'oer | Theory of kinematic modelling and numerical procedures for robot calibration[END_REF]. Thus, (6) becomes:

∆Y = b j=1 w j w j
.∆η j . w j [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF] where w j is the jth column of the matrix W b . It is to be noted that since the columns of W b correspond only to the identifiable parameters, w j = 0.

Classical Calibration Method

This method needs external sensors to measure the location of the terminal frame with respect to the world reference frame.

The nonlinear equation of the calibration is given as:

-1 T n+1 (q, η r ) --1 T n+1 (x) = 0 (8)
The linear differential model defining the deviation of the end-effector location due to the differential error in the geometric parameters can be obtained as [START_REF] Khalil | Calibration of the geometric parameters of robots without external sensors[END_REF][START_REF] Khalil | Identifiable parameters and optimum configurations for robots calibrations[END_REF]:

∆x(q, η) = J(q, η).∆η (9)
with • ∆x represents the (6 × 1) vector of the position and orientation errors (representing the difference between measured and computed -1 T n+1 ) • J is the (6 × N p) Jacobian matrix of frame (n + 1) with respect to the geometric parameters. Its columns can be calculated directly using the elements of the transformation matrices -1 T j as given in Appendix B.

The calibration index of this method is equal to 6. If the external sensor gives only the position of the endpoint, only the first three equations of (9) will be used and the calibration index will be equal to 3.

Calibration Using Relative Location Measure

In this method we need an external sensor giving the relative location of the end-effector frame when moving from configuration q a to q b . Let the (4 × 4) matrix a F b denote this measure, thus:

[ -1 T n+1 (q a , η r )] -1 . -1 T n+1 (q b , η r ) = a F b (10)
The nonlinear calibration relation is given by:

-1 T n+1 (q b , η r ) --1 T n+1 (q a , η r ). a F b = 0 (11)
The differential model of the first-order development is given as:

[J(q b , η) -J2(q a , η, a F b )].∆η = ∆x(q a , q b , η, a F b (12)
J(q a , η) is obtained as given in Appendix B, while J2(q a , η, a F b can be obtained from the expressions of J(q a , η) after replacing the vector -1 P n+1 (q a , η) by the displacement vector of the transformation -1 T n+1 (q a , η). a F b . The right-hand side of ( 12) is the differential error vector between the frames -1 T n+1 (q b , η) and -1 T n+1 (q a , η). a F b .

The calibration index of the method is equal to 6. 3

Calibration Using Distance Measure

This method handles the calibration of the geometric parameters of robots if the Cartesian distance of the endpoint can be measured using an external sensor when moving from configuration q a to q b . Suppose two configurations, q a and q b , where the Cartesian distance between the position of their endpoints is Dr, the nonlinear calibration equation is given as:

[P x(q b , η r ) -P x(q a , η r )] 2 + [P y(q b , η r ) -P y(q a , η r )] 2

+[P z(q b , η r ) -P z(q a , η r )] 2 = Dr 2 (13)
The first-order differential model is given by:

2.[P x(q b ) -P x(q a )].[Jx(q b ) -Jx(q a )]
+2.[P y(q b ) -P y(q a )].[Jy(q b ) -Jy(q a )]

+2.[P z(q b ) -P z(q a )].[Jz(q b ) -Jz(q a )] ∆η = Dr 2 -D 2 (14) where:

• D is the distance between the endpoints of configurations q a and q b , using the nominal parameters and the geometric model • Jx, Jy, and Jz denote respectively, the first, second, and third row of the Jacobian matrix defined in [START_REF] Zhong | A new method for autonomous robot calibration[END_REF]. The calibration index of this method is equal to 1.

Frame Link and Position Link Calibration Methods

The main problem of the previous methods is the need to have an accurate, fast, and inexpensive external sensor to measure the Cartesian variables. The frame link (or position link) method can be used for robots (redundant or not), where for a given location (or position) of the terminal effector, multiple configurations can be obtained from the inverse kinematic model. Let q a and q b represent two configurations giving the same location of the terminal link, then the nonlinear calibration model is given as:

-1 T n+1 (q a , η r ) = -1 T n+1 (q b , η r ) (15)
The first-order differential model gives:

J(q b , η) -J(q a , η) .∆η = ∆x(q a , q b , η)

where J is the Jacobian matrix of frame (n+1) with respect to the geometric parameter variations, as defined in [START_REF] Zhong | A new method for autonomous robot calibration[END_REF]. ∆x is the differential position and orientation vector between the tool locations -1 T n+1 (q a , η) and -1 T n+1 (q b , η).

The calculation of the identifiable parameters can be carried out, as demonstrated in Appendix A, by studying the QR decomposition of an observation matrix W calculated from (16) using a sufficient random couple of configurations, which give the same tool location. Such couple of configurations are obtained by supposing a random configuration q a , then computing the terminal frame location -1 T n+1 (q a ). Finally q b is obtained using the solution of the inverse geometric model for -1 T n+1 (q a ) and such that q b = q a .

In the case of the position link method, we take q a and q b giving the same terminal position such that -1 P n+1 (q a , η r ) --1 P n+1 (q b , η r ). Thus, we use only the first three equations of [START_REF] Khalil | Identifiable parameters and optimum configurations for robots calibrations[END_REF].

The calibration index, in the case of frame link, is equal to 6. It will be reduced to 3 in the case of position link.

Planar Calibration Methods

In this case the calibration will be carried out using the values of the joint positions of a set of configurations of the robot whose endpoints are in the same plane. Different methods based on this technique have been proposed previously [9 13].

Two methods are used in our study.

The First Method: Calibration Using Plane Equation

The general nonlinear equation of the calibration is:

a.P x(q, η r ) + b.P y(q, η r ) + c.P z(q, η r ) + 1 = 0

where a, b, and c represent the plane coefficients P x, P y, and P z represent the position coordinates of the terminal point in the world frame Using a first-order development, the differential model is obtained as:

[P x(q, η) P y(q, η) P z(q, η) a.Jx(q, η) +b.Jy(q, η) + c.Jz(q, η)] .

        ∆a ∆b ∆c ∆η        
= -a.P x(q, η) -b.P y(q, η) -c.P z(q, η) -1 [START_REF] Wampler | An implicit loop method for kinematic calibration and its application to closedchain mechanisms[END_REF] where:

Jx, Jy, and Jz are respectively, the first, second, and third row of the Jacobian matrix defined in [START_REF] Zhong | A new method for autonomous robot calibration[END_REF]. P u(q, η r ), for u = x, y, z represent the u coordinate of the terminal point position at configuration q.

The coefficients of the plane are initialized by calculating the equation of the nearest plane to the terminal points of the given configurations.

If the coefficients of the plane are known, the parameters a, b, and c will be not identified. The corresponding columns and unknowns in (18) will be eliminated.

The calibration index of this method is equal to 1.

Second Method: Calibration using Normal Coordinates

In this method, we make use of the fact that the scalar product of the coordinates of the vector normal to the plane and of any vector between two points (i and j) in the plane is equal to zero [START_REF] Zhong | A new method for autonomous robot calibration[END_REF][START_REF] Maurine | D'eveloppement et mise en oeuvre de m ethodologies d' etalonnage de robots manipulateurs industriels[END_REF]. The main advantage of this procedure is that the coordinates of the normal can be obtained easily using inclinometers.

The system equation is thus:

a. P x(q j , η r ) -P x(q i , η)r) +b. P y(q j , η r ) -P y(q i , η r ) +c. P z(q j , η r ) -P z(q i , η r ) = 0

Using a first-order development for η and assuming that the normal coordinates are known, we obtain: a. Jx(q j ) -Jx(q i ) + b. Jy(q j ) -Jy(q i ) +c. Jz(q j ) -Jz(q i ) .∆η = -a. P x(q j ) -P x(q i ) -b. P y(q j ) -P y(q i ) -c. P z(q j ) -P z(q i ) (

The calibration index of this method is equal to 1.

Evaluation of the Calibration Methods

In this section we study the previous methods on both the Puma and Stanford robots, and obtain some general conclusions.

The comparison focuses on: the identifiable parameters, the condition number of random data, the convergence to the real parameters, and the number of iterations to achieve convergence.

Identifiable Parameters

Tables 3 and4 give the identifiable parameters of the two robots. The parameters indicated by '0' are not identifiable and they have no effect on the identification model. The parameters indicated by 'n' are not identified because they have been regrouped to some other parameters. The regrouped relations can be obtained as given in Appendix A. From Tables 3 and4, the following general remarks are deduced:

1. The location measure method can identify the maximum number of parameters (36 parameters for a Puma robot and 34 for a Stanford robot), which corresponds to the following general relation, in agreement with that given in [START_REF] Everett | Kinematic modeling for robot calibration[END_REF][START_REF] Goswami | Identifying robot parameters using partial pose information[END_REF]: Number of Identifiable Parameters = 4(n + 1) + 2 -2n P + n. This number can be interpreted as: (a) 4 parameters for frames 1, . . . , n + 1 (b) 2 parameters for frame 0 (c) -2n P because two parameters are not identifiable for each prismatic joint (d) n parameters of joint angle gains 2. The parameters of frames 0 and 1 have no effect on the model and cannot be identified when using the following methods: relative location measure, distance measure, point link and frame link. 3. Most of the parameters of frame n and n + 1 are not identifiable in the frame link method. 4. Most of the parameters of frames 0, 1, and 7 are not identifiable in the planar methods. Some of them have been regrouped to other parameters. 5. The effect of the errors on the non-identifiable parameters will be calibrated through the parameters on which they have been regrouped. 6. The parameter β j is not identifiable when α j = 0, that is to say, when the axis j is not parallel to the axis j -1. 7. The offsets of the joint variables 2, . . . , n -1 and all the gains K j are identifiable in all the methods. q n , however, is not identifiable in frame link method, whereas q 1 is not identifiable for the following methods: relative location measure, distance measure, point link, frame link, and unknown planar methods. 8. The parameter r 6 is not identifiable in the point link and plane link methods for the two robots; it represents the scale factor of those closed-loop methods. In the frame link method, the parameter r 6 has no effect and the scale factor is given by r 4 for the Puma manipulator and r 2 for the Stanford arm. It is to be noted that in the case of Stanford arm, the scale factor could be the prismatic variable r 3 (instead of r 6 or r 2 ) if we suppose that the joint gain K 3 is known and has not been identified. 9. The parameters α 7 , θ 7 , and r 7 are not identifiable in the case of the position measure method because the effector is just a point defined by r 6 , θ 6 , and d 7 . For the same reason they are not identifiable for point link and plane link methods.

Relative

Plane Link 1 Plane Link n: non-identifiable parameter, its effect is regrouped on some other parameters 0: non-identifiable parameter having no effect on the model * in this method the coefficient c of the normal to the plane is not identifiable

Simulation Results

The following errors are supposed on the parameters of the Puma robot: These errors give about 30 cm errors on the endpoint position for each robot. This error is very large with respect to the performance of industrial robots. Even in this extreme case, the standard pseudo-inverse works very well.

It is to be noted that these errors are identifiable in all the methods. The unknown vector ∆η contains all of the identifiable parameters of the corresponding method.

The calculation of the precision of the robots after the calibration is done with a set of 200 configurations, randomly distributed in the workspace of the robot, which are not used in the calibration process.

To get reliable results, each method has been tested using 10 different random sets of data. The number of configurations used for each method is taken so that we obtain 120 equations for the calibration, which is about four times the number of parameters.

Simulation without Noise

In this case we assumed exact measuring of q, and X and exact constraint links. The estimation process has been converged to the exact values of the geometric parameters in all the methods in about 5-8 iterations. The accuracy of the robot after calibration is seen to be 10 -14 meters for the position and 10 -14 radians for the orientation, which correspond to the numerical precision.

Simulation with Noise

In order to simulate a realistic calibration process, noisy data have been considered. The noise on the joint positions is randomly generated as equal to half a step of the encoders: we consider an encoder with 100,000 points for each joint, which corresponds to 6.6E-5 radians for a step. In the same way, we add normally zero mean distributed noise on the Cartesian measures and on all the virtual links in the autonomous methods with a variance of 0.1 millimeter on the position and 0.1 milliradians on the orientation.

It is to be noted that in the case of exact geometric parameters, the error on the position of the endpoint due to the imprecision on the joint position values is about 0.3 mm for the Puma robot and 0.2 mm for the Stanford robot.

Tables 5 and6 give the identification results for the Puma and Stanford arms respectively. In all the methods, the convergence of the estimation of the geometric parameters has taken place within 6-13 iterations.

Remarks

From the foregoing simulations, the following remarks are given:

• In general, random data give good condition numbers for all methods, except plane link methods. Therefore, for planar points the user may need to select the location of the plane and the robot configurations by an optimization procedure. • The normalization of the columns of the matrix W b reduces its condition number. • The condition number range is larger for the Stanford robot than for the Puma robot in the case of frame link. This can be explained by the fact that the inverse kinematic model of the Stanford arm has only four solutions, which limits the distribution of the configurations in the workspace. • In general, the number of iterations increases when the condition number increases. This can be especially seen with the plane link methods, which have slower convergence than the other methods. • In the case of joint angle noise, bias errors may result on the identified parameters, but the accuracy of the robot after calibration will be in the same order as the errors. For instance, in the case of the Puma robot, the joint noises give an error on the terminal point of about 0.3, which is in the same range as the errors due to the errors on the identified parameters after calibration in all the methods except in the case of some simulations concerning the plane link methods. • We can observe that the final calibration accuracy of the robot using plane link identification methods with noisy data is less efficient than the other methods. • When the level of noise is as given here, which is higher than typical values on industrial robots, the pseudoinverse method works very well. If needed, statistical methods can take as initial conditions the solution obtained by the standard pseudo-inverse method.

Conclusion

This paper presents a comparison study on the geometric parameter calibration methods of robots. The study is carried out on the Puma and Stanford robots. Comparison of the methods focused on the identifiable parameters, the condition number, the convergence rate, and the precision. Noisy data have been also considered. The identifiable parameters are determined for all of the methods using a numerical method based on QR decomposition, while the estimation of the geometric parameters is carried out using linearized iterative techniques. General conclusions have been obtained concerning the identifiable parameters of the different methods. It has been shown that the iterative-pseudo inverse method with column normalization works very well even in the case of noisy data and large initial errors in the geometric parameters.

Figure 1 .

 1 Figure 1. The geometric parameters of a serial robot.

Figure 2 .

 2 Figure 2. Puma manipulator.

Figure 3 .

 3 Figure 3. Stanford manipulator.

∆α 2 :

 2 -0.05 ∆α 5 : +0.05 ∆d 2 : +0.05 ∆d 3 : -0.05 ∆r 2 : -0.05 ∆r 5 : +0.05 ∆θ 2 : +0.05 ∆θ 4 : +0.05 ∆θ 5 : -0.05 ∆K 1 : +0.05 ∆K 3 : -0.05 ∆K 6 : +0.05 The following errors are supposed on the parameters of the Stanford robot: ∆α 2 : +0.05 ∆α 5 : -0.05 ∆d 2 : +0.05 ∆d 3 : +0.05 ∆r 3 : +0.05 ∆r 5 : +0.05 ∆θ 2 : -0.05 ∆θ 4 : +0.05 ∆θ 5 : -0.05 ∆K -1 : +0.1 ∆K 3 : -0.1 ∆K 6 : -0.05

  1 Plane Link 2 Plane Link 2

		Position Location Location Distance Point Frame	(Known	(Unknown	(Known	(Unknown
	Parameter Measure Measure Measure Measure Link Link	Plane)	Plane)	Plane)	Plane)
	θ7	0			0	0	0	0	0	0	0
	r7	n			n	n	0	n	n	n	n
	β7	n	n	n	n	n	0	n	n	n	n
	K1
	K2
	K3
	K4
	K5
	K6
	Total	31	34	28	25	24	21	27	24	26	24 *

Table 5

 5 Calibration of the Geometric Parameters of the Puma Robot with Noise

		Interval of	Interval of
		maximum	maximum
		position error	orientation error Convergence:	Interval of	Interval of
		on the terminal	on the terminal	Interval of non-normalized normalized
		point (m)	point (rad)	number	condition	condition
		After calibration	After calibration of interations	number	number
	Position
	Measure	1.62E-04 : 3.19E-04 4.01E-04 : 5.90E-04	6 : 6	10.8 : 13.2	6.6 : 9.0
	Location
	Measure	1.69E-04 : 2.74E-04 1.70E-04 : 2.64E-04	6 : 6	10.5 : 14.9	8.3 : 10.1
	Relative
	Location Measure 1.08E-04 : 2.44E-04 1.33E-04 : 2.48E-04	6 : 6	9.9 : 12.3	6.9 : 10.7
	Distance
	Measure	1.72E-04 : 3.22E-04 3.30E-04 : 6.03E-04	6 : 7	17.1 : 23.6	6.9 : 9.1
	Point Link	1.26E-04 : 2.63E-04 3.19E-04 : 6.29E-04	6 : 6	9.9 : 12.5	6.1 : 8.0
	Frame Link	6.89E-05 : 3.92E-04 1.23E-04 : 4.98E-04	6 : 6	9.6 : 16.4	6.9 : 10.7
	Plane Link 1
	(known plane)	2.61E-04 : 7.81E-04 4.48E-04 : 1.07E-03	6 : 7	12.5 : 197.1	8.7 : 31.3
	Plane Link 1
	(unknown plane) 2.56E-04 : 7.22E-04 4.29E-04 : 1.12E-03	7 : 9	18.6 : 472.4	9.4 : 37.9
	Plane Link 2
	(known normal) 3.14E-04 : 1.19E-03 5.79E-04 : 1.65E-03	6 : 7	14.2 : 41.7	9.0 : 27.2
	Plane Link 2
	(unknown normal) 2.73E-04 : 1.21E-03 4.99E-04 : 1.71E-03	7 : 13	32.4 : 4272.2	8.9 : 27.9

Table 6

 6 Calibration of the Geometric Parameters of the Stanford Robot with Noise

		Interval of	Interval of
		maximum	maximum
		position error	orientation error Convergence:	Interval of	Interval of
		on the terminal	on the terminal	Interval of non-normalized normalized
		point (m)	point (rad)	number	condition	condition
		After calibration	After calibration of interations	number	number
	Position
	Measure	1.69E-04 : 3.06E-04 3.60E-04 : 7.14E-04	6 : 7	11.0 : 14.0	6.6 : 8.2
	Location
	Measure	1.46E-04 : 3.51E-04 1.27E-04 : 2.59E-04	6 : 6	14.0 : 20.6	10.6 : 16.6
	Relative
	Location Measure 1.01E-04 : 2.52E-04 6.36E-05 : 1.82E-04	6 : 6	14.6 : 22.4	10.1 : 16.2
	Distance
	Measure	1.79E-04 : 3.80E-04 3.25E-04 : 7.07E-04	6 : 7	19.0 : 26.1	9.6 : 13.7
	Point Link	1.10E-04 : 2.14E-04 1.76E-04 : 5.74E-04	6 : 6	9.2 : 11.6	6.2 : 7.9
	Frame Link	1.75E-04 : 7.82E-04 1.26E-04 : 3.63E-04	6 : 7	50.6 : 98.6	15.6 : 36.4
	Plane Link 1
	(known plane)	2.26E-04 : 6.26E-03 3.69E-04 : 6.66E-03	6 : 9	20.8 : 685.7	10.1 : 60.0
	Plane Link 1
	(unknown plane) 2.22E-04 : 1.85E-03 3.80E-04 : 2.16E-03	6 : 10	18.8 : 204.8	9.4 : 59.5
	Plane Link 2
	(known normal) 3.48E-04 : 1.94E-02 6.56E-04 : 1.73E-02	7 : 11	17.6 : 871.1	9.8 : 45.7
	Plane Link 2
	(unknown normal) 2.31E-04 : 1.60E-03 5.93E-04 : 1.54E-03	7 : 11	21.8 : 316.8	9.7 : 46.0

Let us consider the following overconstrained system of linear equations:

where W is (r × c) matrix with r c. If b is the rank of W, we can write:

where:

We can write:

with β as the (b × (c -b)) matrix with constant elements. Using ( 23), ( 22) can be written as:

with

The solution of the full rank system (24) will yield ∆η b , which is called the identifiable (or base) parameter vector. The matrix β is not needed in the identification process.

Numerically, the study of the identifiable parameters is equivalent to the study of the space spanned by the columns of the matrix W.

The QR decomposition of the matrix W is given as [START_REF] Dongarra | LINPACK user's guide[END_REF][START_REF] Lawson | Solving least squares problems[END_REF]:

where:

Q is an (r × r) orthogonal matrix R is a (c × c) upper triangular matrix 0 i×j is the (i × j) matrix of zeros Theoretically, the non-identifiable parameters are those whose corresponding elements on the diagonal of the matrix R are equal to zero. Assuming τ is the numerical zero, if |R ii | representing the absolute value of the (i, i) element of R is less than τ , the corresponding parameter ∆η i is not identifiable. The numerical zero τ can be taken as [START_REF] Dongarra | LINPACK user's guide[END_REF]:

where is the machine precision.

To calculate the identifiable (base) parameters as a function of the standard parameters, let us permute the columns of W such that the first b columns are independent:

where:

• P is a permutation matrix

The QR decomposition of (W.P) gives:

where R 1 is a (b × b) regular matrix. Then it becomes:

Thus from ( 23) and (29) we get:

Appendix B: Calculation of the Jacobian Matrix

Assuming the (4 × 4) transformation matrix defining frame j relative to the fixed frame as:

The calculation of the columns of the matrix J can be done as follows [START_REF] Hollerbach | A review of kinematic calibration[END_REF][START_REF] Khalil | Mod elisation identi cation et commande des robots[END_REF][START_REF] Khalil | Identifiable parameters and optimum configurations for robots calibrations[END_REF]:

∂q ∂k j where:

• × denotes the vector product

• L j,n+1 is the (3 × 1) displacement vector between the origin of frame j and the origin of frame n + 1 equal to P n+1 -P j

• 0 (3×1) is the (3 × 1) zero vector

• Jq j is equal to Jθ j if joint j is revolute or Jr j if joint j is prismatic All the vectors of (32) refer to the measuring fixed frame.

Biographies