# Global Vision Localization of Indoor Service Robot Based on Improved Iterative Extended Kalman Particle Filter Algorithm

Authors: Bingshan Hu; Qingxiao Yu; Hongliu Yu
Journal: Journal of Sensors
 (2021)

---

## Abstract

Localization is the primary problem of mobile robot navigation. Monte Carlo localization based on particle filter has better accuracy and is easier to implement, but there is also the problem of particle degradation. In this paper, the iterative extended Kalman filter is optimized by the Levenberg-Marquardt optimization method. An improved particle filter algorithm based on the upon optimized iterative Kalman filter is proposed, and the importance probability density function of the particle filter is generated by the maximum posterior probability estimation of the improved iterative Kalman filter. Simulation results of the improved particle filter algorithm show that the algorithm can approximate the state posterior probability distribution more closely with fewer sampled particles under the premise of ensuring sufficient state estimation accuracy. Meanwhile, the computation is reduced and the real-time performance is enhanced. Finally, the algorithm is validated on the indoor mobile service robot. The experimental results show that the localization algorithm’s accuracy meets requirement for real-time localizing of the restaurant service robot.


---

## Body

## 1. Introduction

Localization is the primary problem of mobile robot navigation. Accurate, fast. and stable localization is the premise for mobile robot to perform navigation tasks correctly [1]. The global visual localization uses the global camera to collect the image information of the environment where the mobile robot is located. Through the computer vision algorithm, the pose of the mobile robot in the environment is obtained to realize the localization of the mobile robot [2]. In the process of global localizing by ceiling camera, the implementation of the localizing algorithm does not depend on the mobile robot’s own processor but can run on the upper computer. Meanwhile, localization multiple robots in the same space only requires a set of ceiling cameras, so it is cheaper [3]. This paper studies the localization method of mobile service robots for food delivery in restaurants. Multiple food delivery robots work simultaneously in restaurants and the working environment is relatively fixed, so it is very suitable for global visual localization using a ceiling camera.Common global localization methods for service robots include Kalman filtering localization method [4], Markov localization method [5], and Monte Carlo localization based on particle filter algorithm (MCL) [6]. Kalman filtering localization algorithm used the recursive method to estimate the state of the linear dynamic system with Gaussian noise [7]. However, in the practical application environment, the observation equation is generally nonlinear, so it is necessary to extend the application of Kalman filter through model linearization; thus, an extended Kalman filter (EKF) is generated to solve the state estimation problem in some nonlinear environments [8]. The EKF localization algorithm is simple in structure and has certain precision, so it has been widely used. Due to model errors and noise statistical errors caused by linearization of nonlinear systems, EKF is faced with problems such as filter divergence and insufficient estimation accuracy [9]. Iterative extended Kalman filtering (IEKF) algorithm introduces the iterative filtering theory into the extended Kalman filtering method, updates the linearized observation equation by iterative calculation, compensates the truncation error caused by system linearization effectively, and makes full use of observation information to approximate the optimal estimate of state parameters [10]. IEKF localization method approximates the predictive value of robot state by using an iterative sequence, sacrificing a certain amount of computing time, so as to obtain a good performance of robot state estimation [11].By representing the possible distribution of robot pose with grid points in the state space, Markov localization can realize the global localization of mobile robot by constantly updating the probability distribution of state space points [12]. Compared with EKF localization, Markov localization has better robustness, while EKF localization has better localization accuracy [12]. However, due to the need to calculate the reliability of all grids repeatedly, Markov localization consumes a lot of time, which leads to the poor real-time performance, and it is not conducive to practical application [13].Monte Carlo localization is an extension of Markov method, and the basic principle of this method is to use particle filter (PF) to track the probability distribution of the robot’s possible pose. Monte Carlo localization greatly reduces memory consumption and makes effective use of robot resources. It has better localization accuracy and is easier to implement. Montemerlo et al. proposed a fast simultaneous localization and mapping (SLAM) algorithm based on Rao-Blackwellized particle filter in 2003 [14]. However, when the environmental observation noise of mobile robots is low (i.e., the environmental observation sensor is too accurate, such as the laser sensor), particle set degradation is easy to occur, which makes the SLAM algorithm diverge [15]. It is a simple and widely used method of selecting the state transfer probability density function with prior property as the importance probability density function to solve the phenomenon of particle degradation [15]. Montemerlo et al. proposed the Fast SLAM2.0 algorithm, which uses EKF to integrate the current robot’s environmental observation information into the proposal distribution’s design of the particle filter, so that the particles are concentrated in the high observation likelihood region [16]. Fast SLAM 2.0 can create a better importance density at the cost of introducing inaccuracies due to linearization, as well as making a Gaussian assumption on the form of the posterior density. Ullah et al. designs importance density function by using unscented Kalman filter and proposes an unscented particle Filter (PF-UKF) algorithm. It works well at the cost of large computation burden and when the state estimation error covariance is close to the process noise covariance [17]. The proposal distribution selected by the above algorithm does not contain the latest observation data, which also may result in the performance degradation of the algorithm. The iterated extended Kalman filter can make efficient use of the latest observation and improve the performance of particle filter. Wu et al. use the IEKF to generate proposal distribution in particle filtering framework, and an iterative extended Kalman particle filter localization algorithm (PF-IEKF) is proposed [18]. When the nonlinear system satisfies the condition of local linearization, the IEKF iterative updating sequence is equivalent to the Gauss-Newton (G-N) iterative estimation [19]. However, the G-N method has some problems in practical application, such as convergence and precision of calculation results. In order to improve the convergence and precision of the G-N algorithm, the Levenberg-Marquardt (LM) optimization method is usually used to optimize the IEKF iterative sequence to improve the convergence and precision of IEKF iterative sequence [20]. Li et al. proposed a robot vision positioning method based on iterative Kalman particle filter. Experimental results show that the accuracy and real-time performance of robot mobile navigation can meet the requirements [21].Inspired by the Gauss-Newton optimization algorithm, this paper improves iterative extended Kalman filtering using the Levenberg-Marquardt optimization method to improve the convergence and calculation accuracy of IEKF iterative sequence. The maximum posterior probability estimation of improved IEFK is used to generate the importance probability density function of particle filter, and an improved iterative extended Kalman particle filter (PF-UIEKF) localization algorithm is proposed. Then, the effectiveness of the improved particle filter algorithm is verified by simulation. At last, it is applied to the global visual localization on the restaurant service robot and verified by experiment.


## 2. Modeling of the Indoor Service Robot

### 2.1. Hardware of the Robot

The research object of this paper is the restaurant service robot, whose task is to pick up cooked dishes and send them to the designated table autonomously. Therefore, the restaurant service robot must have good autonomous mobility and can move freely in any direction in the restaurant environment through a motion mechanism. At the same time, in order to make the restaurant service robot take the place of a human to complete the task of picking up and delivering dishes independently, the service robot must also be configured with corresponding manipulating mechanism.As shown in Figure1, the restaurant service robot is composed of three main parts: the robot arm, the head, and the motion mechanism. The robot motion mechanism adopts the differential drive with two driving wheels and is also equipped with two universal wheels which play a supporting role. Two robotic arms are installed on the left and right sides of the robot torso. Each robotic arm has three degrees of freedom to complete the operation of grasping the dish. To achieve the global localization of the service robot, a global camera is installed on the ceiling of the restaurant. The global vision controller receives video signals from the ceiling global camera over a wireless network. The vision processing algorithm is used to identify and measure the pose of the camera relative to the color mark installed at the top of the service robot head. At last, the global localization algorithm is run to determine the robot’s coordinates in the global coordinate system.Figure 1
Structure diagram of the service robot.


### 2.2. Kinematic Model of the Robot

In order to analyze the kinematic model of the service robot, the following assumptions are made. The two driving wheels are rigid bodies with the same size and no deformation. The line of the two wheels’ axis is perpendicular to the front and rear motion direction of the platform. The wheel face of the driving wheel is perpendicular to the ground and maintains point contact. The drive wheel only does pure rolling and will not produce axial movement. The influence of the drive wheel thickness on the movement of the moving platform is ignored.Based on the above assumptions, as shown in Figure2, the pose of the restaurant service robot can be expressed as M=XM,YM,θT, where XM,YM is the projection coordinate of the central point M of the driving wheel axis in the motion plane, and θ is the course angle of the robot. Suppose that the radius of the two driving wheels is r, the axial distance is d, and the rotational angular velocity of left and right wheel is ωL and ωR, respectively, then the localization formula of the robot based on the kinematics model can be obtained as follows:(1)Xi+1=Xi+r2cosθi∫titi+1ωLt+ωRtdt,Yi+1=Yi+r2sinθi∫titi+1ωLt+ωRtdt,θi+1=θi+rd∫titi+1ωRt−ωLtdt.Figure 2
The diagram of the two-wheeled differential driven robot.By discretizing Equation (1) and adding the model noise of the system, the state equation of the restaurant service robot can be obtained as follows:(2)xk=Xk−1+Trcosθk−1ωRk+ωLk2Yk−1+Trsinθk−1ωRk+ωLk2θk−1+TrωRk−ωLkd+wk−1.In Equation (2), xk=Xk,Yk,θkT.wk is the model error, which is the Gaussian white noise with zero mean value. T is the sampling period.


### 2.3. Localization of the Robot through the Global Ceiling Camera

As shown in Figure3, the global coordinate system (GCS), robot coordinate system (RCS), and camera coordinate system (CCS) are established for the service robot firstly.Figure 3
Establishment of the robot coordinate system.In Figure4, 0,0 denotes the origin of the camera coordinate system, HC denotes the height of the ceiling camera relative to the floor; HR denotes the height of the robot; XRC,YRC is the coordinates of the robot in CCS; XHC,YHC is color mark’s coordinates in CCS. According to the geometric relationship, the following equation can be obtained:(3)aa+b=HRHC,ba+b=YRCYHC,ba+b=XRCXHC.Figure 4
The globe vision camera model of the robot.Then, the following results can be calculated as,(4)XRC=HC−HRHC×XHC,YRC=HC−HRHC×YHC.From Equation (4), it is obvious that the coordinates of the robot in CCS XRC,YRC is only dependent on XHC,YHC, which is the color mark’s coordinates in CCS. To get the (XHC, YHC), the image of the monocular camera mounted on the ceiling is acquired and transmitted to the robot global vision controller through the wireless network. The image is firstly processed to remove the image distortion. Then, the color segmentation method is used to locate the robot in the image and distinguish the color mark on the head in the hue-saturation-value (HSV) color space. After distinguishing and extracting the color mark, the coordinates of the color mark in CCS can be obtained [22].At last, the coordinates of the robot in GCS can be solved according to coordinate transformation [22], and the coordinates of the robot in GCS can be calculated as(5)X′Y′=TXRCXRC=XRCcosθ+YRCsinθ+mcosθ+nsinθ−XRCsinθ+YRCcosθ−msinθ+ncosθ.In Equation (5), θ is the rotation angle of the camera coordinate system relative to the global coordinate system; m and n are the offset distances of the camera coordinate system relative to the global coordinate system. In order to determine the current pose of the robot, the current course angle of the service robot is also required. Gyroscope is installed in the robot, and the initial course angle of the service robot is obtained according to the gyroscope firstly. In order to overcome the influence of gyroscope cumulative error on the course angle when the robot moves for a long time or a long distance, it is necessary to correct the course angle using global vision periodically.If the global coordinates of the service robot at two consecutive moments are obtained by the global visual localization method, and they areXt−1,Yt−1 and Xt,Yt, then the current course angle of the robot is(6)θt=arctanXt−Xt−1Yt−Yt−1.By combining Equations (5) and (6), the observation equation of service robot based on global vision after discretization can be obtained as follows:(7)yk=Xk+mcosθk+Yk+nsinθk−Xk+msinθk+Yk+ncosθkarctanXk−Xk−1Yk−Yk−1+vk.In Equation (7), yk=(XRCk, YRCk, θRCk)T, and v(k) is the observed noise of the global ceiling camera of the robot system.


## 2.1. Hardware of the Robot

The research object of this paper is the restaurant service robot, whose task is to pick up cooked dishes and send them to the designated table autonomously. Therefore, the restaurant service robot must have good autonomous mobility and can move freely in any direction in the restaurant environment through a motion mechanism. At the same time, in order to make the restaurant service robot take the place of a human to complete the task of picking up and delivering dishes independently, the service robot must also be configured with corresponding manipulating mechanism.As shown in Figure1, the restaurant service robot is composed of three main parts: the robot arm, the head, and the motion mechanism. The robot motion mechanism adopts the differential drive with two driving wheels and is also equipped with two universal wheels which play a supporting role. Two robotic arms are installed on the left and right sides of the robot torso. Each robotic arm has three degrees of freedom to complete the operation of grasping the dish. To achieve the global localization of the service robot, a global camera is installed on the ceiling of the restaurant. The global vision controller receives video signals from the ceiling global camera over a wireless network. The vision processing algorithm is used to identify and measure the pose of the camera relative to the color mark installed at the top of the service robot head. At last, the global localization algorithm is run to determine the robot’s coordinates in the global coordinate system.Figure 1
Structure diagram of the service robot.


## 2.2. Kinematic Model of the Robot

In order to analyze the kinematic model of the service robot, the following assumptions are made. The two driving wheels are rigid bodies with the same size and no deformation. The line of the two wheels’ axis is perpendicular to the front and rear motion direction of the platform. The wheel face of the driving wheel is perpendicular to the ground and maintains point contact. The drive wheel only does pure rolling and will not produce axial movement. The influence of the drive wheel thickness on the movement of the moving platform is ignored.Based on the above assumptions, as shown in Figure2, the pose of the restaurant service robot can be expressed as M=XM,YM,θT, where XM,YM is the projection coordinate of the central point M of the driving wheel axis in the motion plane, and θ is the course angle of the robot. Suppose that the radius of the two driving wheels is r, the axial distance is d, and the rotational angular velocity of left and right wheel is ωL and ωR, respectively, then the localization formula of the robot based on the kinematics model can be obtained as follows:(1)Xi+1=Xi+r2cosθi∫titi+1ωLt+ωRtdt,Yi+1=Yi+r2sinθi∫titi+1ωLt+ωRtdt,θi+1=θi+rd∫titi+1ωRt−ωLtdt.Figure 2
The diagram of the two-wheeled differential driven robot.By discretizing Equation (1) and adding the model noise of the system, the state equation of the restaurant service robot can be obtained as follows:(2)xk=Xk−1+Trcosθk−1ωRk+ωLk2Yk−1+Trsinθk−1ωRk+ωLk2θk−1+TrωRk−ωLkd+wk−1.In Equation (2), xk=Xk,Yk,θkT.wk is the model error, which is the Gaussian white noise with zero mean value. T is the sampling period.


## 2.3. Localization of the Robot through the Global Ceiling Camera

As shown in Figure3, the global coordinate system (GCS), robot coordinate system (RCS), and camera coordinate system (CCS) are established for the service robot firstly.Figure 3
Establishment of the robot coordinate system.In Figure4, 0,0 denotes the origin of the camera coordinate system, HC denotes the height of the ceiling camera relative to the floor; HR denotes the height of the robot; XRC,YRC is the coordinates of the robot in CCS; XHC,YHC is color mark’s coordinates in CCS. According to the geometric relationship, the following equation can be obtained:(3)aa+b=HRHC,ba+b=YRCYHC,ba+b=XRCXHC.Figure 4
The globe vision camera model of the robot.Then, the following results can be calculated as,(4)XRC=HC−HRHC×XHC,YRC=HC−HRHC×YHC.From Equation (4), it is obvious that the coordinates of the robot in CCS XRC,YRC is only dependent on XHC,YHC, which is the color mark’s coordinates in CCS. To get the (XHC, YHC), the image of the monocular camera mounted on the ceiling is acquired and transmitted to the robot global vision controller through the wireless network. The image is firstly processed to remove the image distortion. Then, the color segmentation method is used to locate the robot in the image and distinguish the color mark on the head in the hue-saturation-value (HSV) color space. After distinguishing and extracting the color mark, the coordinates of the color mark in CCS can be obtained [22].At last, the coordinates of the robot in GCS can be solved according to coordinate transformation [22], and the coordinates of the robot in GCS can be calculated as(5)X′Y′=TXRCXRC=XRCcosθ+YRCsinθ+mcosθ+nsinθ−XRCsinθ+YRCcosθ−msinθ+ncosθ.In Equation (5), θ is the rotation angle of the camera coordinate system relative to the global coordinate system; m and n are the offset distances of the camera coordinate system relative to the global coordinate system. In order to determine the current pose of the robot, the current course angle of the service robot is also required. Gyroscope is installed in the robot, and the initial course angle of the service robot is obtained according to the gyroscope firstly. In order to overcome the influence of gyroscope cumulative error on the course angle when the robot moves for a long time or a long distance, it is necessary to correct the course angle using global vision periodically.If the global coordinates of the service robot at two consecutive moments are obtained by the global visual localization method, and they areXt−1,Yt−1 and Xt,Yt, then the current course angle of the robot is(6)θt=arctanXt−Xt−1Yt−Yt−1.By combining Equations (5) and (6), the observation equation of service robot based on global vision after discretization can be obtained as follows:(7)yk=Xk+mcosθk+Yk+nsinθk−Xk+msinθk+Yk+ncosθkarctanXk−Xk−1Yk−Yk−1+vk.In Equation (7), yk=(XRCk, YRCk, θRCk)T, and v(k) is the observed noise of the global ceiling camera of the robot system.


## 3. Improved Iterative Extended Kalman Particle Filter Localization Algorithm

### 3.1. Improved Iterative Extended Kalman Filtering

The robot state Equation (2) and observation Equation (7) above can be simplified as(8)xk=fxk−1,uk+wkyk=hxk,uk+vk.In Equation (8), f∗ and h∗ are the state transfer function and observation function of the robot system, respectively. xk and yk are the state value and observed value of the system at time k, respectively. wk and vk are the state noise and observation noise of the system, and their covariances are Q and R, respectively. uk is the known input. When EKF performs nonlinear transformation, it uses Taylor series expansion and ignores higher-order terms, which will introduce truncation error and make the performance of state estimation worse. In order to reduce nonlinear error and improve the performance of nonlinear filtering, IEKF algorithm introduces the iterative filtering theory into the EKF algorithm. In the state update stage of IEKF, multistep iteration is adopted, and the observation equation is linearized for many times. The whole filtering process is repeated to achieve the optimal state estimation finally [23].The main steps of IEKF algorithm are as follows:(1)
State value prediction. Suppose the estimated state value at timek−1 is x^k−1∣k−1, and the known input state uk and the system state equation are used to obtain the estimated state value x^k−1∣k−1 at time k. Then, the Jacobian matrix Fk−1 of the system state equation and the error covariance matrix Pk∣k−1 of the system state are calculated, respectively, to complete the state prediction process(2)
Observation value acquisition. The state estimation valuex^k∣k−1 at time k and the system observation equation were used to obtain the predicted observation value y^k∣k−1 at time k(3)
Updating the predicted value iteratively using the observation value. Assuming that there areN iterations, the Jacobian matrix Hki and the filtered gain matrix Kki of the system observation equation at the ith iteration are calculated, respectively. Then use Hki and Kki to update the error covariance matrix. The difference value yk−yki of the observed variable is calculated and the status is updatedIt can be seen that if the updatedx^k∣k−1 after observation is taken as the initial value of iteration and the observation equation is linearized again, the obtained state estimation value will be more accurate. The IEKF algorithm is formed by carrying out this iterative process several times.Bell and Cathey analyzed the relationship between the iterative sequence of IEKF and the nonlinear least squares Gauss-Newton method. When the nonlinear system satisfies the condition of local linearization, the iterative update sequence of iteratively extended Kalman filter is equivalent to Gauss-Newton iterative estimation [19]. Due to the introduction of errors such as linearization, the system’s state space equation may not be completely consistent with the actual observed data, and solutions using the Gauss-Newton method usually do not yield stable results. It also cannot guarantee the convergence of the state observation updated estimation results, and the estimated value of the covariance matrix is lower than the true value, thus affecting the effective use of the observation information.A modified Levenberg-Marquardt method can correct the original iterative sequence by introducing a damping factor, which improves the convergence of the algorithm and the accuracy of the calculation results [24]. Therefore, in this paper, the Levenberg-Marquardt method is used to adjust the covariance matrix in the IEKF to ensure the convergence of the estimated error. The core of this method is that the damping factor λ is used to modify the prediction error covariance matrix during each iteration, that is, to adjust the error covariance matrix as follows:(9)P~k/k−1=I−Pk/k−1Pk/k−1+λ−1I−1Pk/k−1.Then, an iterative observation update is performed with a modifiedP~k/k−1. The implementation steps of the improved IEKF algorithm based on the Levenberg-Marquardt optimization method are as follows:(1)
Initial value setting. The calculation formula of the initial value is(10)k=0,x^0=Ex0,P0=Ex^0−x0x∧0−x0T,wherex^0 is the estimated value of initial state and P0 is the initial state error covariance matrix(2)
Prediction of state. Based on the estimated state value at timeK−1 and known input uk, the predicted state value x^k/k−1, Jacobian matrix Fk−1, and error covariance matrix Pk∣k−1 of the system at time K are calculated as follows:(11)x^k/k−1=fx^k−1/k−1,uk,Fk−1=∂f∂xx^k−1/k−1,uk,Pk∣k−1=Fk−1Pk−1/k−1Fk−1T+Q(3)
Obtain observation information. Based onx^k/k−1 and the observation equation of the system, the predicted observation value of the system can be calculated as follows:(12)y^k/k−1=hx^k/k−1,uk(4)
Iterative updates of status. Assume thatN iterations are performed. For the ith iteration 1≤i≤N, the following calculation is performed successivelyThe Jacobian matrix of the system state is calculated as(13)Fk−1i=∂f∂xx^k−1/k−1i,uk.The predicted state and error covariance matrix of the system are calculated as(14)x^k/k−1i=fx^k−1/k−1i,uk,Pk/k−1i=Fk−1iPk−1/k−1Fk−1iT+Q.The modified error covariance matrix is optimized using the Levenberg-Marquardt method(15)P~k/k−1i=I−Pk/k−1iPk/k−1i+λ−1I−1Pk/k−1i.The Jacobian matrix of the observation equation is calculated as(16)Hki=∂h∂xx^k/ki,uk.The filter gain matrix is calculated as(17)Kki=P~k/k−1iHkiHkiP~k/k−1iHkiT+R−1.The updated error covariance matrix is,(18)Pk/ki=I−KkiHkiP~k/k−1i.The predicted value of the observed value of the system is(19)y^ki=hx^k/k−1i,uk.The system status estimated value is updated as(20)x^k/ki=x^k/k−1i+Kkiyk−y^ki.When the number of iterations reaches the set maximum number of iterationsN or the error between two consecutive iterations is less than the set minimum error value, the iteration process is stopped. The iteration result is taken as the status update value at time k,(21)x^k/k=x^k/ki,Pk/k=Pk/ki.In the above equation, “^” represents the estimated state. The subscript k/k−1 represents the prior estimate. The subscript k/k is the posterior estimate.


### 3.2. Improved Iterative Extended Kalman Particle Filter Algorithm

Compared with the iteratively extended Kalman particle filter algorithm, the difference of the improved particle filter algorithm proposed in this paper is that, the importance probability density function of particle filter is generated by iterative Kalman filter based on L-M optimization. The importance probability density function reduces the linearization error of system. It is more consistent with the actual posterior probability distribution of the state variable. The implementation steps of the improved iteratively extended Kalman particle filter algorithm (PF-UIEKF) are as follows:(1)
Initializing the particle set. The initial state and prior probability densitypx0 of the robot were sampled Ns times, and Ns particle points were obtained(2)
Fork=1,2,⋯⋯N, N is the number of samplesFori=1,2,⋯⋯Ns,(i)
Calculating the Jacobian matrixFki of the state equation using Equation (13)(ii)
The improved IEKF algorithm is used to update the particle set using Equation (14)Forj=1,2,⋯⋯c(a)
The Jacobian matrixHkji of the observation equation of the system is calculated using Equation (16)(b)
The state covariancePkji is updated by Equation (18)(c)
Equation (17) is used to update the state gain Kkji(d)
The state estimationx∧kji is updated according to Equation (20)End for(iii)
Construct the importance probability density function(iv)
(22)xki∼qxkixk−1i,yk=Nx^kci,PkciIn Equation (22), qxkixk−1i,yk is importance sampling density function.(v)
Ns times of sampling are conducted to generate state prediction particle set(vi)
Calculating the particle weightωki of the sample set using the following equation(23)ωki∝ωk−1ipyk∣xkipxki∣xk−1iqxki∣xk−1i,ykIn equation (23), pyk∣xki is the likelihood function of the observed value yk.End for(3)
The sample point weight is normalized as(24)ωki=ωki∑j=1Nωkj(4)
To represent the degree of particle degradation, the effective particle set indexNeff in particle filter is introduced [25](25)Neff=1∑i=1Nωki2Preset a threshold valueNth, and Nth=αN0<α<1. If Neff>Nth, step 6 is entered, otherwise step 5 is entered.(5)
Sequence importance resampling. Polynomial resampling is performed on the original particle set. Particles with larger sample weights are propagated, while particles with smaller sample weights are eliminated, and a new set of particles is obtained(26)S^t=x^t1,x^t2,⋯,x^tNsThe sample weights of the new particle set are reallocated,ωki=1/Ns. At last, step 7 is entered.(6)
The state estimation value of the system at timek is output according to the particle setSt=xt1,xt2,⋯,xtNs(7)
Makek=k+1 and go to step 2, end for


### 3.3. Global Vision Localization Based on PF-UIEKF Algorithm

In order to reduce the navigation cost of mobile service robot in restaurant environment, the robot global localization is realized based on the restaurant environment image collected by the global ceiling camera in this paper. The localization method based on global vision has been introduced in Section2.3 in detail, and it can realize the global localization of the robot theoretically. However, due to the influence of noise in the restaurant environment, the measured value of global vision often contains large non-Gaussian noise. It seriously affects the localization accuracy and fails to meet the mobile navigation requirements. Here, a robot global visual localization algorithm based on PF-UIEKF is proposed.According to Equations (2) and (7), the robot system’s state model and observation model is obtained. Equations (17), (18), and (20) are used to obtain the state estimation and covariance iterative update equation of the robot’s position using IEKF based on L-M optimization. Then, the importance probability density function is generated by using the maximum posterior probability estimation of the obtained robot state, and the improved particle filter algorithm is used to obtain more accurate global coordinates of the robot. Figure 5 shows the implementation steps of the global visual localization algorithm of the restaurant service robot based on the PF-UIEKF.Figure 5
The improved iterative extended Kalman particle filter localization algorithm.


## 3.1. Improved Iterative Extended Kalman Filtering

The robot state Equation (2) and observation Equation (7) above can be simplified as(8)xk=fxk−1,uk+wkyk=hxk,uk+vk.In Equation (8), f∗ and h∗ are the state transfer function and observation function of the robot system, respectively. xk and yk are the state value and observed value of the system at time k, respectively. wk and vk are the state noise and observation noise of the system, and their covariances are Q and R, respectively. uk is the known input. When EKF performs nonlinear transformation, it uses Taylor series expansion and ignores higher-order terms, which will introduce truncation error and make the performance of state estimation worse. In order to reduce nonlinear error and improve the performance of nonlinear filtering, IEKF algorithm introduces the iterative filtering theory into the EKF algorithm. In the state update stage of IEKF, multistep iteration is adopted, and the observation equation is linearized for many times. The whole filtering process is repeated to achieve the optimal state estimation finally [23].The main steps of IEKF algorithm are as follows:(1)
State value prediction. Suppose the estimated state value at timek−1 is x^k−1∣k−1, and the known input state uk and the system state equation are used to obtain the estimated state value x^k−1∣k−1 at time k. Then, the Jacobian matrix Fk−1 of the system state equation and the error covariance matrix Pk∣k−1 of the system state are calculated, respectively, to complete the state prediction process(2)
Observation value acquisition. The state estimation valuex^k∣k−1 at time k and the system observation equation were used to obtain the predicted observation value y^k∣k−1 at time k(3)
Updating the predicted value iteratively using the observation value. Assuming that there areN iterations, the Jacobian matrix Hki and the filtered gain matrix Kki of the system observation equation at the ith iteration are calculated, respectively. Then use Hki and Kki to update the error covariance matrix. The difference value yk−yki of the observed variable is calculated and the status is updatedIt can be seen that if the updatedx^k∣k−1 after observation is taken as the initial value of iteration and the observation equation is linearized again, the obtained state estimation value will be more accurate. The IEKF algorithm is formed by carrying out this iterative process several times.Bell and Cathey analyzed the relationship between the iterative sequence of IEKF and the nonlinear least squares Gauss-Newton method. When the nonlinear system satisfies the condition of local linearization, the iterative update sequence of iteratively extended Kalman filter is equivalent to Gauss-Newton iterative estimation [19]. Due to the introduction of errors such as linearization, the system’s state space equation may not be completely consistent with the actual observed data, and solutions using the Gauss-Newton method usually do not yield stable results. It also cannot guarantee the convergence of the state observation updated estimation results, and the estimated value of the covariance matrix is lower than the true value, thus affecting the effective use of the observation information.A modified Levenberg-Marquardt method can correct the original iterative sequence by introducing a damping factor, which improves the convergence of the algorithm and the accuracy of the calculation results [24]. Therefore, in this paper, the Levenberg-Marquardt method is used to adjust the covariance matrix in the IEKF to ensure the convergence of the estimated error. The core of this method is that the damping factor λ is used to modify the prediction error covariance matrix during each iteration, that is, to adjust the error covariance matrix as follows:(9)P~k/k−1=I−Pk/k−1Pk/k−1+λ−1I−1Pk/k−1.Then, an iterative observation update is performed with a modifiedP~k/k−1. The implementation steps of the improved IEKF algorithm based on the Levenberg-Marquardt optimization method are as follows:(1)
Initial value setting. The calculation formula of the initial value is(10)k=0,x^0=Ex0,P0=Ex^0−x0x∧0−x0T,wherex^0 is the estimated value of initial state and P0 is the initial state error covariance matrix(2)
Prediction of state. Based on the estimated state value at timeK−1 and known input uk, the predicted state value x^k/k−1, Jacobian matrix Fk−1, and error covariance matrix Pk∣k−1 of the system at time K are calculated as follows:(11)x^k/k−1=fx^k−1/k−1,uk,Fk−1=∂f∂xx^k−1/k−1,uk,Pk∣k−1=Fk−1Pk−1/k−1Fk−1T+Q(3)
Obtain observation information. Based onx^k/k−1 and the observation equation of the system, the predicted observation value of the system can be calculated as follows:(12)y^k/k−1=hx^k/k−1,uk(4)
Iterative updates of status. Assume thatN iterations are performed. For the ith iteration 1≤i≤N, the following calculation is performed successivelyThe Jacobian matrix of the system state is calculated as(13)Fk−1i=∂f∂xx^k−1/k−1i,uk.The predicted state and error covariance matrix of the system are calculated as(14)x^k/k−1i=fx^k−1/k−1i,uk,Pk/k−1i=Fk−1iPk−1/k−1Fk−1iT+Q.The modified error covariance matrix is optimized using the Levenberg-Marquardt method(15)P~k/k−1i=I−Pk/k−1iPk/k−1i+λ−1I−1Pk/k−1i.The Jacobian matrix of the observation equation is calculated as(16)Hki=∂h∂xx^k/ki,uk.The filter gain matrix is calculated as(17)Kki=P~k/k−1iHkiHkiP~k/k−1iHkiT+R−1.The updated error covariance matrix is,(18)Pk/ki=I−KkiHkiP~k/k−1i.The predicted value of the observed value of the system is(19)y^ki=hx^k/k−1i,uk.The system status estimated value is updated as(20)x^k/ki=x^k/k−1i+Kkiyk−y^ki.When the number of iterations reaches the set maximum number of iterationsN or the error between two consecutive iterations is less than the set minimum error value, the iteration process is stopped. The iteration result is taken as the status update value at time k,(21)x^k/k=x^k/ki,Pk/k=Pk/ki.In the above equation, “^” represents the estimated state. The subscript k/k−1 represents the prior estimate. The subscript k/k is the posterior estimate.


## 3.2. Improved Iterative Extended Kalman Particle Filter Algorithm

Compared with the iteratively extended Kalman particle filter algorithm, the difference of the improved particle filter algorithm proposed in this paper is that, the importance probability density function of particle filter is generated by iterative Kalman filter based on L-M optimization. The importance probability density function reduces the linearization error of system. It is more consistent with the actual posterior probability distribution of the state variable. The implementation steps of the improved iteratively extended Kalman particle filter algorithm (PF-UIEKF) are as follows:(1)
Initializing the particle set. The initial state and prior probability densitypx0 of the robot were sampled Ns times, and Ns particle points were obtained(2)
Fork=1,2,⋯⋯N, N is the number of samplesFori=1,2,⋯⋯Ns,(i)
Calculating the Jacobian matrixFki of the state equation using Equation (13)(ii)
The improved IEKF algorithm is used to update the particle set using Equation (14)Forj=1,2,⋯⋯c(a)
The Jacobian matrixHkji of the observation equation of the system is calculated using Equation (16)(b)
The state covariancePkji is updated by Equation (18)(c)
Equation (17) is used to update the state gain Kkji(d)
The state estimationx∧kji is updated according to Equation (20)End for(iii)
Construct the importance probability density function(iv)
(22)xki∼qxkixk−1i,yk=Nx^kci,PkciIn Equation (22), qxkixk−1i,yk is importance sampling density function.(v)
Ns times of sampling are conducted to generate state prediction particle set(vi)
Calculating the particle weightωki of the sample set using the following equation(23)ωki∝ωk−1ipyk∣xkipxki∣xk−1iqxki∣xk−1i,ykIn equation (23), pyk∣xki is the likelihood function of the observed value yk.End for(3)
The sample point weight is normalized as(24)ωki=ωki∑j=1Nωkj(4)
To represent the degree of particle degradation, the effective particle set indexNeff in particle filter is introduced [25](25)Neff=1∑i=1Nωki2Preset a threshold valueNth, and Nth=αN0<α<1. If Neff>Nth, step 6 is entered, otherwise step 5 is entered.(5)
Sequence importance resampling. Polynomial resampling is performed on the original particle set. Particles with larger sample weights are propagated, while particles with smaller sample weights are eliminated, and a new set of particles is obtained(26)S^t=x^t1,x^t2,⋯,x^tNsThe sample weights of the new particle set are reallocated,ωki=1/Ns. At last, step 7 is entered.(6)
The state estimation value of the system at timek is output according to the particle setSt=xt1,xt2,⋯,xtNs(7)
Makek=k+1 and go to step 2, end for


## 3.3. Global Vision Localization Based on PF-UIEKF Algorithm

In order to reduce the navigation cost of mobile service robot in restaurant environment, the robot global localization is realized based on the restaurant environment image collected by the global ceiling camera in this paper. The localization method based on global vision has been introduced in Section2.3 in detail, and it can realize the global localization of the robot theoretically. However, due to the influence of noise in the restaurant environment, the measured value of global vision often contains large non-Gaussian noise. It seriously affects the localization accuracy and fails to meet the mobile navigation requirements. Here, a robot global visual localization algorithm based on PF-UIEKF is proposed.According to Equations (2) and (7), the robot system’s state model and observation model is obtained. Equations (17), (18), and (20) are used to obtain the state estimation and covariance iterative update equation of the robot’s position using IEKF based on L-M optimization. Then, the importance probability density function is generated by using the maximum posterior probability estimation of the obtained robot state, and the improved particle filter algorithm is used to obtain more accurate global coordinates of the robot. Figure 5 shows the implementation steps of the global visual localization algorithm of the restaurant service robot based on the PF-UIEKF.Figure 5
The improved iterative extended Kalman particle filter localization algorithm.


## 4. Algorithm Verification and Experiment

### 4.1. Simulation and Analysis

To verify the effect of the improved iterative extended Kalman particle filter algorithm, the following standard nonlinear system model is adopted in simulation experiments, and the standard nonlinear system state equation is(27)xk+1=1+sinωπk+ϕ1xk+νk.And its observation equation is(28)zk=ϕ2xk2+nk,t≤30,ϕ3xk−2+nk,t>30.The system process noisevk obeys the gamma distribution, and vk∼ζα3,2. The observed noise obeys Gaussian distribution and nk∼N0,0.00001. The following are given parameters in a nonlinear system model: ω=0.04, Φ1=0.5, Φ2=0.2, and Φ3=0.5. The number of particles used in the simulation experiment is N=20, and the observation time is T=60s. The polynomial resampling method is used to resample the particle set. State variables were initialized again before each simulation experiment, and 100 independent simulation experiments were repeated. The output of the algorithm is the mean value of the particle set, and then, the state estimation of the nonlinear system can be expressed as(29)x¯k=1N∑i=1Nxki.In order to explain the state estimation performance of the PF-UIEKF algorithm more intuitively, simulation experiments of various particle filtering algorithms such as standard particle filter (PF), extended Kalman particle Filter (PF-EKF), and iterative extended Kalman particle filter (PF-IEKF) were also carried out for this nonlinear model, and the performance of state estimation was compared and analyzed. Figure6 shows posterior means estimated by four different nonlinear filtering algorithms (PF, PF-EKF, PF-IEKF, and PF-UIEKF) in an independent experiment.Figure 6
Posterior means estimated by four different nonlinear filtering algorithms.In this paper, the root mean square error (RMSE) of state estimation value is used to reflect the state estimation accuracy of each particle filter algorithm, and its value can be calculated from the following formula.(30)RMSE=1n∑k=1Nx¯k−xk2.Table1 summarizes the state estimation performance data for each particle filter. The mean and variance of the RMSE of state estimation value of each particle filter algorithm are given, and the execution time of each algorithm is also given. It can be clearly seen from Table 1 that the state estimation accuracy using the PF-UIEKF algorithm is the best. Although PF-UIEKF particle filter algorithm has a longer execution time than others, the difference is not significant.Table 1
RMSE and execution time of different algorithms (20 particles).
AlgorithmRMSEExecution time (s)MeanVariancePF0.815610.0320850.69404PF-EKF0.393390.0262610.70998PF-IEKF0.157230.0119120.71868PF-UIEKF0.0122761.1699e-0050.77623In order to analyze the state estimation performance of the improved particle filter algorithm furtherly, the number of particles used in the simulation experiment was increased to 50 and 200. Other parameters remain the same. 100 independent simulation experiments were performed on all particle filter algorithms again, and initialization was required before each simulation experiment. Tables2 and 3 summarize the state estimation performance data for each particle filter when the number of particles is 50 and 200.Table 2
RMSE and execution time of different algorithms (50 particles).
AlgorithmRMSEExecution time (s)MeanVariancePF0.656810.0361780.70529PF-EKF0.291220.0102150.83844PF-IEKF0.0308930.00282031.1196PF-UIEKF0.00635662.5751e−0061.6613Table 3
RMSE and execution time of different algorithms (200 particles).
AlgorithmRMSEExecution time (s)MeanVariancePF0.459810.0578050.79223PF-EKF0.289790.0164742.6039PF-IEKF0.00621121.5201e−0063.5898PF-UIEKF0.00511362.3769e−0066.0593As can be seen from Table2, when the number of particles rises from 20 to 50, the state estimation performance of each particle filter is greatly improved with the increase of the number of particles. The state estimation accuracy using the PF-UIEKF algorithm is also the best. However, the increase in the number of particles also led to the increase in the execution time of the PF-UIEKF algorithm, which reached 1.6613 s. As can be seen from Table 3, with the number of particles increasing to 200, the state estimation performance of each particle filter is still improved. However, the performance improvement of PF-UIEKF is not obvious compared with other particle filtering algorithms. At this time, PF-IEKF nearly has the same state estimation accuracy, mean, and variance of the RMSE of state estimation value as the PF-UIEKF particle filter. The estimation accuracy of these two particle filtering algorithms is the best and far better than that of other filtering algorithms. However, at this time, the execution time of the PF-UIEKF algorithm is greatly increased, which is nearly twice that of the PF-IEKF algorithm.Compared with Table1–3, PF-UIEKF can still get more accurate state estimation under the condition of fewer particles (N=20), and its estimation performance is much better than that of PF-IEKF algorithm when N=50. Moreover, the execution time of the PF-UIEKF algorithm is shorter than that of PF-IEKF when N=50. The improved particle filter algorithm can approach the posterior probability distribution of the state more closely with fewer sampled particles under the premise of ensuring sufficient state estimation accuracy, which reduces the computation and enhances the real-time performance.


### 4.2. Experiments

In this paper, the following localization experiment is designed: the robot starting from the global coordinate150,150, moving in a counterclockwise direction for one turn, and finally returning to the vicinity of the starting point. During the movement of the robot, the localization algorithm based on PF-UIEKF proposed in the paper is used to carry out on-line for real-time localization of the robot.Firstly, the global image in Figure7 is obtained after correcting the original image acquired by the global ceiling camera. The upper left corner of Figure 7 is the origin of the coordinate system CCS, and the origin of the coordinate system GCS is located in the upper right corner of the carpet area in Figure 7. According to the color characteristics of the service robot, an appropriate threshold is set in the HSV color space of the image by using the color segmentation method to detect the robot region in the image and the color mark above the robot head [22]. After localizing the position of the color mark of the robot’s head in the global image, the global coordinates of the service robot can be calculated according to the method in Section 2.3. Since the height of the global camera and the robot are both unchanged and known in advance (Table 4 lists the given parameters required for the global visual localization).Figure 7
Global vision image processing.Table 4
The used parameters in this experiment.
HeightWidthXHC,YHCHCHR576768(356.5, 72.67)350150The current coordinates of the robot are updated every 100 ms, and the results were recorded. In order to display the experimental effect more intuitively, 55 sample points of the actual motion trajectory were extracted at certain intervals, and the robot localization results and the actual coordinate positions were displayed in Figure8.Figure 8
The movement trajectory of the service robot.The experimental results show that the maximum localization error of the robot is 4.8 cm. The mean localization error was 3.74 cm. That is, the localization accuracy of the proposed global visual localizing algorithm based on PF-UIEKF can be kept within 5 cm, which fully meets the requirements of mobile navigation and localization of the restaurant service robots.


## 4.1. Simulation and Analysis

To verify the effect of the improved iterative extended Kalman particle filter algorithm, the following standard nonlinear system model is adopted in simulation experiments, and the standard nonlinear system state equation is(27)xk+1=1+sinωπk+ϕ1xk+νk.And its observation equation is(28)zk=ϕ2xk2+nk,t≤30,ϕ3xk−2+nk,t>30.The system process noisevk obeys the gamma distribution, and vk∼ζα3,2. The observed noise obeys Gaussian distribution and nk∼N0,0.00001. The following are given parameters in a nonlinear system model: ω=0.04, Φ1=0.5, Φ2=0.2, and Φ3=0.5. The number of particles used in the simulation experiment is N=20, and the observation time is T=60s. The polynomial resampling method is used to resample the particle set. State variables were initialized again before each simulation experiment, and 100 independent simulation experiments were repeated. The output of the algorithm is the mean value of the particle set, and then, the state estimation of the nonlinear system can be expressed as(29)x¯k=1N∑i=1Nxki.In order to explain the state estimation performance of the PF-UIEKF algorithm more intuitively, simulation experiments of various particle filtering algorithms such as standard particle filter (PF), extended Kalman particle Filter (PF-EKF), and iterative extended Kalman particle filter (PF-IEKF) were also carried out for this nonlinear model, and the performance of state estimation was compared and analyzed. Figure6 shows posterior means estimated by four different nonlinear filtering algorithms (PF, PF-EKF, PF-IEKF, and PF-UIEKF) in an independent experiment.Figure 6
Posterior means estimated by four different nonlinear filtering algorithms.In this paper, the root mean square error (RMSE) of state estimation value is used to reflect the state estimation accuracy of each particle filter algorithm, and its value can be calculated from the following formula.(30)RMSE=1n∑k=1Nx¯k−xk2.Table1 summarizes the state estimation performance data for each particle filter. The mean and variance of the RMSE of state estimation value of each particle filter algorithm are given, and the execution time of each algorithm is also given. It can be clearly seen from Table 1 that the state estimation accuracy using the PF-UIEKF algorithm is the best. Although PF-UIEKF particle filter algorithm has a longer execution time than others, the difference is not significant.Table 1
RMSE and execution time of different algorithms (20 particles).
AlgorithmRMSEExecution time (s)MeanVariancePF0.815610.0320850.69404PF-EKF0.393390.0262610.70998PF-IEKF0.157230.0119120.71868PF-UIEKF0.0122761.1699e-0050.77623In order to analyze the state estimation performance of the improved particle filter algorithm furtherly, the number of particles used in the simulation experiment was increased to 50 and 200. Other parameters remain the same. 100 independent simulation experiments were performed on all particle filter algorithms again, and initialization was required before each simulation experiment. Tables2 and 3 summarize the state estimation performance data for each particle filter when the number of particles is 50 and 200.Table 2
RMSE and execution time of different algorithms (50 particles).
AlgorithmRMSEExecution time (s)MeanVariancePF0.656810.0361780.70529PF-EKF0.291220.0102150.83844PF-IEKF0.0308930.00282031.1196PF-UIEKF0.00635662.5751e−0061.6613Table 3
RMSE and execution time of different algorithms (200 particles).
AlgorithmRMSEExecution time (s)MeanVariancePF0.459810.0578050.79223PF-EKF0.289790.0164742.6039PF-IEKF0.00621121.5201e−0063.5898PF-UIEKF0.00511362.3769e−0066.0593As can be seen from Table2, when the number of particles rises from 20 to 50, the state estimation performance of each particle filter is greatly improved with the increase of the number of particles. The state estimation accuracy using the PF-UIEKF algorithm is also the best. However, the increase in the number of particles also led to the increase in the execution time of the PF-UIEKF algorithm, which reached 1.6613 s. As can be seen from Table 3, with the number of particles increasing to 200, the state estimation performance of each particle filter is still improved. However, the performance improvement of PF-UIEKF is not obvious compared with other particle filtering algorithms. At this time, PF-IEKF nearly has the same state estimation accuracy, mean, and variance of the RMSE of state estimation value as the PF-UIEKF particle filter. The estimation accuracy of these two particle filtering algorithms is the best and far better than that of other filtering algorithms. However, at this time, the execution time of the PF-UIEKF algorithm is greatly increased, which is nearly twice that of the PF-IEKF algorithm.Compared with Table1–3, PF-UIEKF can still get more accurate state estimation under the condition of fewer particles (N=20), and its estimation performance is much better than that of PF-IEKF algorithm when N=50. Moreover, the execution time of the PF-UIEKF algorithm is shorter than that of PF-IEKF when N=50. The improved particle filter algorithm can approach the posterior probability distribution of the state more closely with fewer sampled particles under the premise of ensuring sufficient state estimation accuracy, which reduces the computation and enhances the real-time performance.


## 4.2. Experiments

In this paper, the following localization experiment is designed: the robot starting from the global coordinate150,150, moving in a counterclockwise direction for one turn, and finally returning to the vicinity of the starting point. During the movement of the robot, the localization algorithm based on PF-UIEKF proposed in the paper is used to carry out on-line for real-time localization of the robot.Firstly, the global image in Figure7 is obtained after correcting the original image acquired by the global ceiling camera. The upper left corner of Figure 7 is the origin of the coordinate system CCS, and the origin of the coordinate system GCS is located in the upper right corner of the carpet area in Figure 7. According to the color characteristics of the service robot, an appropriate threshold is set in the HSV color space of the image by using the color segmentation method to detect the robot region in the image and the color mark above the robot head [22]. After localizing the position of the color mark of the robot’s head in the global image, the global coordinates of the service robot can be calculated according to the method in Section 2.3. Since the height of the global camera and the robot are both unchanged and known in advance (Table 4 lists the given parameters required for the global visual localization).Figure 7
Global vision image processing.Table 4
The used parameters in this experiment.
HeightWidthXHC,YHCHCHR576768(356.5, 72.67)350150The current coordinates of the robot are updated every 100 ms, and the results were recorded. In order to display the experimental effect more intuitively, 55 sample points of the actual motion trajectory were extracted at certain intervals, and the robot localization results and the actual coordinate positions were displayed in Figure8.Figure 8
The movement trajectory of the service robot.The experimental results show that the maximum localization error of the robot is 4.8 cm. The mean localization error was 3.74 cm. That is, the localization accuracy of the proposed global visual localizing algorithm based on PF-UIEKF can be kept within 5 cm, which fully meets the requirements of mobile navigation and localization of the restaurant service robots.


## 5. Conclusions

According to the global localization requirement of indoor restaurant service robot which is driven by two differential wheels, the kinematics model of the service robot and the localization model through the global ceiling camera are established in this paper firstly. Then, the iterative extended Kalman filter is optimized and modified by the Levenberg-Marquardt optimization method. A particle filter algorithm based on an improved iterative Kalman filter is proposed, and the importance probability density function of the particle filter is generated by the maximum posterior probability estimation of the improved iterative Kalman filter. A typical nonlinear system model is used to simulate the improved particle filter algorithm. The simulation results show that the particle filter algorithm can obtain higher estimation accuracy and improve the real-time performance of state estimation by using fewer particles. Finally, a method of global localization of indoor restaurant service robot based on improved iterative extended Kalman particle filter algorithm and global visual localization is proposed and verified by experiments. The experimental results show that the localizing accuracy of the proposed global localizing algorithm can be kept within 5 cm, which meets the requirements of mobile navigation and localizing of restaurant service robots.



---

*Source: 8819917-2021-07-09.xml*