Mechanical arm self-learning compliant control method based on sensor online calibration

Abstract

The invention belongs to the technical field of mechanical arm control, and particularly relates to a mechanical arm self-learning compliant control method based on sensor online calibration. Comprising the following steps: 1, denoising data acquired by a six-dimensional force sensor by using a Kalman filter; step 2, acquiring actual data of contact force by utilizing online gravity compensation; 3, the environment interaction model of the robot polishing process of the complex surface is analyzed; 4, identifying initial parameters of the environment through RLS; and 5, by means of the obtained data, a mechanical arm self-learning compliant control method based on sensor online calibration is provided, a reference trajectory is corrected, and zero-error steady-state tracking is achieved. The force tracking effect of the mechanical arm in a complex unknown environment can be improved, and then the grinding precision of the robot is improved to a large extent.

Description

Technical Field

[0001] The invention belongs to the technical field of robot arm control, in particular to a robot arm self-learning compliant control method based on online sensor calibration. Background Art

[0002] In recent years, robots have been widely used in the manufacturing industry, especially in the polishing process. Traditional manual polishing requires high technical requirements and the smoke and dust generated are harmful to health. In addition, the quality of polishing depends on the skills of the operator, and it is difficult to ensure product quality after long working hours. In contrast, robots can maintain high efficiency continuously and accurately position through preset instructions to improve production efficiency. However, robots also face challenges in the polishing process, especially for high-precision, multi-curvature, and easily deformed workpieces. Stable force control tracking becomes the key to achieving high-quality polishing.

[0003] Impedance control is to build a dynamic model of force and position so that the end of the robot can build a dynamic model between position and force on the workpiece surface, which can achieve synchronous control of robot motion and contact force. By adjusting the impedance parameters of the robot, constant force tracking on the workpiece surface can be achieved. However, when the environment and stiffness parameters are uncertain, the steady-state error of the force between the robot and the environment is not zero. Therefore, there is a great need to study how to improve grinding and polishing accuracy by optimizing impedance control. Summary of the invention

[0004] To solve the above problems, the present invention provides a self-learning compliant control method for a robotic arm based on online sensor calibration, which can improve the force tracking effect of the robotic arm in a complex and unknown environment, thereby greatly improving the robot's grinding accuracy.

[0005] The technical solution of the present invention is described as follows in conjunction with the accompanying drawings:

[0006] A method for self-learning compliant control of a robot arm based on online sensor calibration comprises the following steps:

[0007] Step 1: Use Kalman filter to denoise the data collected by the six-dimensional force sensor;

[0008] Step 2: Use online gravity compensation to obtain actual data of contact force;

[0009] Step 3: Analyze the environmental interaction model of the robot polishing process of complex surfaces;

[0010] Step 4: Identify initial parameters of the environment by RLS;

[0011] Step 5: Using the acquired data, a self-learning compliant control method for the robot arm based on online sensor calibration is proposed to correct the reference trajectory and achieve zero-error steady-state tracking.

[0012] Furthermore, the specific method of step one is as follows:

[0013] 11) Select Kalman filtering to filter the data collected by the six-dimensional force sensor, as shown below:

[0014] x k =A·x k-1 +B·u k +w k-1

[0015] z k =H·x k +v k

[0016] Where A is the state transfer matrix; B is the control input matrix; u k is the control vector; H is the state observation matrix; w k-1 is the process noise; v k is the measured noise, which is Gaussian white noise; w k-1 and v k Both are normally distributed; corner labels k and k-1 are time steps;

[0017] 12) The state is predicted as follows:

[0018]

[0019] In the formula, is the prior state estimate, estimated forward from time step k-1 to k;

[0020] 13) Update the status as follows:

[0021]

[0022] In the formula, is a posterior state estimate;

[0023] 14) Perform Kalman gain calculation as follows:

[0024]

[0025] Where K is the prediction error in the optimal state estimation divided by the measurement error;

[0026] 15) Perform covariance update as follows:

[0027]

[0028] Where P k is the estimated mean square error; is the prior mean square error matrix;

[0029] 16) Update the process noise covariance as follows:

[0030]

[0031] In the formula, is the prediction covariance matrix; Q is the process noise covariance matrix.

[0032] Furthermore, the specific method of step 2 is as follows:

[0033] 21) Obtain the force and torque data measured by the six-dimensional force sensor installed at the end of the robot. The data consists of three parts, as follows:

[0034] a) The sensor's own system error is that the reading under no-load condition is not zero. Assuming that the sensor's own system error is a constant, it is called "zero point". Because the sensor has a "zero drift" phenomenon, the value is not constant, expressed as: [F ox F oy F oz M ox M oy M ox ] T ;

[0035] In the formula, F ox is the error force in the x direction; F oy is the error force in the y direction; F oz is the error force in the z direction; M ox is the error torque rotating around the x direction; M oy is the error torque rotating around the y direction; M ox is the error torque rotating around the z direction;

[0036] b) The gravity of the end effector. The gravity of the end effector connected to the sensor acts on the gravity of the six-dimensional force sensor, which changes with the posture of the end effector and is expressed as:

[0037] [F Gx F Gy F Gz M Gx M Gy M Gx ] T ;

[0038] In the formula, F Gx is the gravity of the end effector in the x direction; F Gy is the gravity of the end effector in the y direction; F Gz is the gravity of the end effector in the z direction; M Gx is the torque of the end effector in the x direction; MGy is the torque of the end effector in the y direction; M Gx is the torque of the end effector in the z direction;

[0039] c) The external force on the robot end effector is the relative output of the force detected by the six-dimensional force sensor and is expressed as:

[0040] [F ex F ey F ez M ex M ey M ez ] T ;

[0041] In the formula, F ex is the force on the end effector in the x direction; F ey is the force on the end effector in the x direction; F ez is the force on the end effector in the y direction; M ex is the torque of the end effector in the z direction; M ey is the torque of the end effector in the y direction; M ez is the torque of the end effector in the z direction;

[0042] 22) Obtain the force and torque detected by the six-dimensional force sensor as follows:

[0043]

[0044] In the formula, F x F is the force in the x direction detected by the six-dimensional force sensor; y F is the force in the y direction detected by the six-dimensional force sensor; z M is the force in the z direction detected by the six-dimensional force sensor; x M is the force in the x direction detected by the six-axis force sensor; y M is the force in the y direction detected by the six-dimensional force sensor; z The force in the z direction detected by the six-dimensional force sensor;

[0045] 23) Obtain the installation posture of the robot; the installation posture of the robot represents the rotation matrix of the robot base coordinate system B relative to the world coordinate system G; the world coordinate system is defined as 0X o Y o Z0, let the Z0 axis direction be vertically upward, the robot base coordinate system is defined as 00X1Y1Z1, the origin of B and G coincide; the world coordinate system G is transformed to the base coordinate system B by rotating V around the Y0 axis first, and then rotating U around the current X0 axis. The posture transformation matrix is ​​expressed as:

[0046]

[0047] Where U is the angle of rotation around the current X0 axis to the base coordinate system B;

[0048] V is the angle of rotation of the world coordinate system G around the Y0 axis;

[0049] 24) Obtain the posture of the robot's end effector; the posture of the robot's end effector represents the rotation matrix of the robot's end effector coordinate system T relative to the base coordinate system B; the industrial robot's connecting rod length and offset are obtained by measurement, and the joint angle is directly obtained by the robot control system. Assuming that the posture of the robot's sixth joint is consistent with the posture of the end effector, the posture of the end effector is solved by the forward kinematics formula;

[0050] 25) According to the composite transformation relationship of the conversion matrix, the rotation matrix of the sensor coordinate system S relative to the world coordinate system G is as follows:

[0051]

[0052] In the formula, is the rotation matrix of the base coordinate system B relative to the world coordinate system G;

[0053] is the rotation matrix of coordinate system T relative to coordinate system B;

[0054] is the rotation matrix of coordinate system S relative to coordinate system T;

[0055] 26) According to the Cartesian transformation relationship of the robot static force, the transformation relationship between force and torque in the G and S coordinate systems is obtained as follows:

[0056]

[0057] In the formula, S F S and S M S are the force and torque components of the gravity of the end effector in the sensor coordinate system S, respectively; is the rotation matrix of coordinate system G relative to coordinate system S; G F G , G M G is the gravity G of the end effector expressed in the world coordinate system G; S P G × is the coordinate of the center of gravity of the end effector in the sensor coordinate system, as shown below:

[0058]

[0059] Where P zIt represents the z-axis of the end effector in the sensor coordinate system;

[0060] P x is the x-axis representation of the end effector in the sensor coordinate system;

[0061] P y is the y-axis representation of the end effector in the sensor coordinate system;

[0062] Simplified to:

[0063]

[0064] Further simplified to:

[0065] S M S = S P G × S F G

[0066]

[0067] Combined

[0068]

[0069] Thus, the following formula is obtained:

[0070]

[0071] In the formula, k1, k2, and k3 are intermediate variables;

[0072] The definitions of intermediate variables k1, k2, and k3 are as follows:

[0073]

[0074] The above formula is simplified to:

[0075] M S =F·P

[0076] Where M S , F represents the data stacking amount of sensor force and sensor torque obtained by sampling N (N≥3) times of posture transformation of the robot end effector; P is the coordinate of the center of gravity of the end effector;

[0077] The calculation method of the center of gravity coordinate P of the end effector is as follows:

[0078] P=(F(N)· T F(N) -1 ·F(N)· T M S (N)

[0079] Where F(N) is the data stacking amount of the sensor force; M S (N) is the data stacking amount of the sensor torque;

[0080] Thus, the center of gravity coordinates (P x , P y , P z ) and intermediate variables k1, k2, k3;

[0081] By formula:

[0082]

[0083]

[0084] Simplified to:

[0085]

[0086] In the formula, k4, k5, and k6 are intermediate variables;

[0087] Among them, the intermediate variables k4, k5, and k6 are defined as follows:

[0088]

[0089] Further simplified to:

[0090] F s =R·A

[0091] Where Fs and R represent the data stack of sensor force and robot joint angle obtained by sampling N (N ≥ 3) times of the robot end effector posture transformation; A is the representation of the gravity of the end effector in the sensor coordinate system.

[0092] The calculation method of the gravity of the end effector in the sensor coordinate system is as follows: A = (R (N) T R(N) -1 ·R(N) T F S (N)

[0093] Where, Fs and R represent the data stacking of sensor force and robot joint angle obtained by sampling N (N≥3) times of the robot end effector posture transformation; A is identified by the least squares method, and k4, k5, and k6 are calculated;

[0094] 27) Calculate the gravity G of the end effector and the installation angle (U, V) of the robot base as:

[0095]

[0096]

[0097] Where G is the gravity of the end tool after the robot body and the end tool are installed; the installation angle (U, V) of the robot base is a constant value;

[0098] 28) Perform online compensation for the gravity and zero point of the end effector; the online compensation divides the six-dimensional force sensor calibration stage into the sensor zero point calibration period 0~T1 and the end effector gravity and zero point compensation period T1~T2 for processing: first, at a certain sampling time t within the period 0~T1, the robot is controlled to move freely, and the joint variables obtained by the robot controller are [θ1(t)θ2(t)θ3(t)θ4(t)θ5(t)θ6(t)] T , R(t) is calculated in real time, and the force and torque data sampled by the six-dimensional force sensor are [F(t)M(t)] T , gravity compensation value of end effector [F G (t)M G (t)] T for:

[0099]

[0100]

[0101] In the formula, S F G (k) is the representation of the gravity compensation value of the end effector in the S coordinate system; S M G (S) is the representation of the gravity compensation value of the end effector in the S coordinate system; is the rotation matrix of coordinate system G relative to coordinate system S; S P G is the center of gravity coordinate of the end effector in the sensor coordinate system. ; is the rotation matrix of coordinate system G relative to coordinate system S; G F G is the gravity compensation value of the end effector; G M G is the gravity compensation value of the end effector; F0(t) and M0(t) are the zero-point data of the sensor after each power-on, which are constant values ​​after a period of fluctuation; F(t) and M(t) are the force and torque data sampled by the six-dimensional force sensor; F G (t), M G (t) is the gravity compensation value of the end effector;

[0102] Define the zero point data of the sensor as [F0(t)M0(t)] T , then:

[0103]

[0104] In the formula, S F G (t) is the representation of the gravity compensation value of the end effector in the S coordinate system; S M G (t) is the representation of the gravity compensation value of the end effector in the S coordinate system; is the rotation matrix of coordinate system G relative to coordinate system S; S P G is the center of gravity coordinate of the end effector in the sensor coordinate system; is the rotation matrix of coordinate system G relative to coordinate system S;

[0105] At a certain time k in the end effector gravity and zero point compensation period T1~T2, the robot is controlled to move freely, and the joint variables of the robot are [θ1(k)θ2(k)θ3(k)θ4(k)θ5(k)θ6(k)] T , calculate R(t) in real time, and calculate the gravity compensation data of the end effector online as [F G (k)M G (k)] T , according to the sensor zero point is a constant value [F0 M0] T ,have:

[0106]

[0107] Where, F(t) is the force data sampled by the six-dimensional force sensor; M(t) is the torque data sampled by the six-dimensional force sensor;

[0108] According to the online compensation results, the external force / torque is accurately detected. The external force and torque detected by the six-dimensional force sensor at the k sampling time is defined as [F e (k)M e (k)] T :

[0109]

[0110] In the formula, F e (k) is the force compensation error of the six-dimensional force sensor; M e (k) is the torque compensation error of the six-dimensional force sensor; F(k) is the torque data sampled by the six-dimensional force sensor; M(s) is the torque data sampled by the six-dimensional force sensor; F G (k) is the gravity compensation value of the end effector; M G (k) is the gravity compensation value of the end effector; F0 is the constant value of the force at the zero point of the sensor; M0 is the constant value of the torque at the zero point of the sensor.

[0111] Furthermore, the specific method of step three is as follows:

[0112] 31) Establish the impedance model as follows:

[0113]

[0114] Where ΔF is the force error and K is the stiffness X c is the end position X r is the reference trajectory, B is the damping, M is the inertia term, X c is the command position trajectory; X r is the reference position trajectory; is the derivative of position is the reference speed of the sixth axis of the robot; is the second derivative of position is the acceleration of the sixth axis of the robot;

[0115] 32) Assuming that the position control method is ideal, the command position trajectory X c =Measurement position trajectory x m , the position error caused by the contact force is expressed as:

[0116] E=X c -X r =X m -X e ;

[0117] Where, X r is the reference trajectory X m Command trajectory X e For environmental location

[0118] Assuming the environment is rigid, use the stiffness K e The spring is represented by x m =x c ;

[0119] 33) The contact force between the robot and the environment is expressed as:

[0120] F e =K e (X c -X r )=K e (X m -X e );

[0121] In the formula, F e is the contact force; K e is the environmental stiffness;

[0122] 34) The stiffness and damping of the environmental parameters are estimated by least squares. In the environmental model, the environmental model in a specific direction is considered. The environmental model is a scalar expression, b e , k e ,x,f e B e , K e , X, F e Damping, stiffness, displacement and force in a certain direction in:

[0123]

[0124] In the formula, is the derivative of the position error; Δx is the position error;

[0125] The discrete changes are:

[0126]

[0127] Where T is the sampling time; z is a complex variable;

[0128] is equivalent to:

[0129]

[0130] Furthermore, the specific method of step 4 is as follows:

[0131] 41) The RLS algorithm gives θ(k) as follows:

[0132]

[0133]

[0134] In the formula, G(k) is the gain factor; λ is the forget factor, which is 0.9; p(k) is the parameter estimation covariance matrix, and the recursive calculation formula is:

[0135]

[0136] 42) When the change of θ(k) value is less than the set threshold, it is considered that θ(k) value is stable. After obtaining the θ(k) value, the calculation formula of a1 and a2 is used to obtain (b e , k e ) are estimated to be:

[0137]

[0138]

[0139] Furthermore, the specific method of step five is as follows:

[0140] 51) Perform Laplace transform on the impedance model as follows:

[0141]

[0142] Where ΔF is the force error; M is the mass; E is the position error; B is the damping; K is the stiffness;

[0143] 52) When it is a one-dimensional vector, due to the expected contact force f d and the reference position trajectory x r is a constant, so The following conditions are met:

[0144]

[0145] f e =k e (x m -x e )

[0146]

[0147]

[0148] 53) The steady-state error of the force tracking function is as follows:

[0149]

[0150] In the formula, e ss is the steady-state error; k e is the environmental stiffness; k is the stiffness; x e is the environmental position; x r is the reference trajectory; f r is the reference force;

[0151] 54) In order to achieve zero-error steady-state tracking, make:

[0152]

[0153] 55) Estimate command position trajectory x c :

[0154]

[0155] In the formula, is the estimated value of the environmental position; f d For the power of expectation; is the estimated value of the environmental stiffness;

[0156] 56) Calculate the actual contact force f e , as shown below:

[0157]

[0158] In the formula, x m For end position

[0159] Define k x =k e ·x e ,

[0160] In the formula, is the estimated value of the actual contact force;

[0161]

[0162] 57) Define the positive definite symmetric matrix Q1 as:

[0163]

[0164] 58) The Lyapunov function is defined as:

[0165] V=Q1 T PQ1

[0166]

[0167] Where V is the Lyapunov function; P is a diagonal positive definite matrix; q1 and q2 are positive constants. When , it tends to be stable;

[0168] Compute the derivative of the Lyapunov function as follows:

[0169]

[0170] In the formula, is the derivative of the Lyapunov function;

[0171] 59) It is deduced that the derivative of the Lyapunov function must be less than or equal to zero, so we have for:

[0173]

[0174]

[0175]

[0176]

[0177] The beneficial effects of the present invention are:

[0178] 1) In terms of robot force control, the vibration problem of the workpiece when the robot is polishing the workpiece surface makes the traditional impedance control feedback untimely. To address this problem, the present invention proposes a robot arm self-learning compliant control method based on online sensor calibration; Kalman filtering is performed on the grinding force data to obtain the actual grinding force data, and then gravity compensation is performed on the grinding force data to obtain the actual grinding force, which solves the problem of large force / position tracking error when the robot is polishing the workpiece surface. The reference trajectory of the impedance control model is dynamically compensated by designing an adaptive law, which effectively constrains the contact force during workpiece polishing, reduces contact jitter, and enhances the robustness of the control system;

[0179] 2) The present invention enables the robot arm to better polish in unknown environments and provides more reliable system stability; it enables the robot arm to better adapt to complex environments; it improves the flexibility of the robot tool end, fully proving that it has practicality and universality. BRIEF DESCRIPTION OF THE DRAWINGS

[0180] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings required for use in the embodiments are briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without creative work.

[0181] Figure 1 It is a schematic diagram of the principle of the self-learning compliant control method of the robot arm based on the online calibration of the sensor;

[0182] Figure 2 Schematic diagram for modeling contact forces;

[0183] Figure 3 This is the robot force tracking effect under the plane;

[0184] Figure 4 This is a schematic diagram of the robot force tracking effect on an inclined plane;

[0185] Figure 5 Schematic diagram of the robot force tracking effect under the curved surface. DETAILED DESCRIPTION

[0186] The present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments. It is to be understood that the specific embodiments described herein are only used to explain the present invention, rather than to limit the present invention. It should also be noted that, for ease of description, only parts related to the present invention, rather than all structures, are shown in the accompanying drawings.

[0187] Embodiment 1

[0188] This embodiment provides a method for self-learning compliance control of a robot arm based on online sensor calibration, comprising the following steps:

[0189] Step 1: Use the Kalman filter to denoise the data collected by the six-dimensional force sensor, as follows:

[0190] 11) During the polishing process, the force signal received by the force sensor is mixed with Gaussian noise. Kalman filtering is selected to filter the data collected by the six-dimensional force sensor, as shown below:

[0191] x k =A·x k-1 +B·u k +w k-1

[0192] z k =H·x k +v k

[0193] Where A is the state transfer matrix; B is the control input matrix; u k is the control vector; H is the state observation matrix; w k-1 is the process noise; v k is the measured noise, which is Gaussian white noise; w k-1 and v k Both obey normal distribution; Gaussian white noise with process covariance Q P(w)~N(0,Q); Gaussian white noise with measurement noise covariance P(v)~N(0,R);

[0194] Corner labels k+1, k and k-1 are all time steps;

[0195] 12) The state is predicted as follows:

[0196]

[0197] In the formula, is the prior state estimate, estimated forward from time step k-1 to k;

[0198] 13) Update the status as follows:

[0199]

[0200] In the formula, is a posterior state estimate;

[0201] 14) Perform Kalman gain calculation as follows:

[0202]

[0203] Where K is the prediction error in the optimal state estimation divided by the measurement error;

[0204] 15) Perform covariance update as follows:

[0205]

[0206] Where P k is the estimated mean square error; is the prior mean square error matrix;

[0207] 16) Update the process noise covariance as follows:

[0208]

[0209] In the formula, is the prediction covariance matrix; Q is the process noise covariance matrix.

[0210] Step 2: Use online gravity compensation to obtain actual contact force data, as follows:

[0211] 21) Obtain the force and torque data measured by the six-dimensional force sensor installed at the end of the robot. The data consists of three parts, as follows:

[0212] a) The sensor's own system error is that the reading under no-load condition is not zero. Assuming that the sensor's own system error is a constant, it is called "zero point". Because the sensor has a "zero drift" phenomenon, the value is not constant, expressed as: [F ox F oy F oz M ox M oy M ox ] T ;

[0213] In the formula, F ox is the error force in the x direction; F oy is the error force in the y direction; F oz is the error force in the z direction; M ox is the error torque rotating around the x direction; M oy is the error torque rotating around the y direction; M ox is the error torque rotating around the z direction;

[0214] b) The gravity of the end effector. The gravity of the end effector connected to the sensor acts on the gravity of the six-dimensional force sensor, which changes with the posture of the end effector and is expressed as:

[0215] [F Gx F Gy F Gz M GxM Gy M Gx ] T ;

[0216] In the formula, F Gx is the gravity of the end effector in the x direction; F Gy is the gravity of the end effector in the y direction; F Gz is the gravity of the end effector in the z direction; M Gx is the torque of the end effector in the x direction; M Gy is the torque of the end effector in the y direction; M Gx is the torque of the end effector in the z direction;

[0217] c) The external force on the robot end effector is the relative output of the force detected by the six-dimensional force sensor and is expressed as:

[0218] [F ex F ey F ez M ex M ey M ex ] T ;

[0219] In the formula, F ex is the force on the end effector in the x direction; F ey is the force on the end effector in the x direction; F ez is the force on the end effector in the y direction; M ex is the torque of the end effector in the z direction; M ey is the torque of the end effector in the y direction; M ez is the torque of the end effector in the z direction;

[0220] 22) Obtain the force and torque detected by the six-dimensional force sensor as follows:

[0221]

[0222] In the formula, F x F is the force in the x direction detected by the six-axis force sensor; y F is the force in the y direction detected by the six-dimensional force sensor; z M is the force in the z direction detected by the six-dimensional force sensor; x M is the force in the x direction detected by the six-dimensional force sensor; y M is the force in the y direction detected by the six-dimensional force sensor; z The force in the z direction detected by the six-dimensional force sensor;

[0223] 23) Obtain the installation posture of the robot; the installation posture of the robot represents the rotation matrix of the robot base coordinate system B relative to the world coordinate system G; the world coordinate system is defined as 0X o Y o Z0, let the Z0 axis direction be vertically upward, the robot base coordinate system is defined as 00X1Y1Z1, the origin of B and G coincide; the world coordinate system G is transformed to the base coordinate system B by rotating V around the Y0 axis first, and then rotating U around the current X0 axis. The posture transformation matrix is ​​expressed as:

[0224]

[0225] Where U is the angle of rotation around the current X0 axis to the base coordinate system B;

[0226] 24) Obtain the posture of the robot's end effector; the posture of the robot's end effector represents the rotation matrix of the robot's end effector coordinate system T relative to the base coordinate system B; the industrial robot's connecting rod length and offset are obtained by measurement, and the joint angle is directly obtained by the robot control system. Assuming that the posture of the robot's sixth joint is consistent with the posture of the end effector, the posture of the end effector is solved by the forward kinematics formula;

[0227] 25) According to the composite transformation relationship of the conversion matrix, the rotation matrix of the sensor coordinate system S relative to the world coordinate system G is as follows:

[0228]

[0229] In the formula, is the rotation matrix of the base coordinate system B relative to the world coordinate system G;

[0230] is the rotation matrix of coordinate system T relative to coordinate system B;

[0231] is the rotation matrix of coordinate system S relative to coordinate system T;

[0232] 26) According to the Cartesian transformation relationship of the robot static force, the transformation relationship between force and torque in the G and S coordinate systems is obtained as follows:

[0233]

[0234] In the formula, S F S and S M S are the force and torque components of the gravity of the end effector in the sensor coordinate system S, respectively; is the rotation matrix of coordinate system G relative to coordinate system S; G F G , G MG is the gravity G of the end effector expressed in the world coordinate system G; S P G × is the coordinate of the center of gravity of the end effector in the sensor coordinate system, as shown below:

[0235]

[0236] In the formula,

[0237] P z is the z-axis representation of the end effector in the sensor coordinate system; P x is the x-axis of the end effector in the sensor coordinate system; P y is the y-axis representation of the end effector in the sensor coordinate system; simplified as:

[0238]

[0239] Further simplified to:

[0240]

[0241] Combined

[0242]

[0243] Thus, the following formula is obtained:

[0244]

[0245] In the formula, k1, k2, and k3 are intermediate variables;

[0246] The definitions of intermediate variables k1, k2, and k3 are as follows:

[0247]

[0248] The above formula is simplified to:

[0249] M S =F·P

[0250] Where M S , F represents the data stacking amount of sensor force and sensor torque obtained by sampling N (N≥3) times of posture transformation of the robot end effector; P is the coordinate of the center of gravity of the end effector;

[0251] The calculation method of the center of gravity coordinate P of the end effector is as follows:

[0252] P=(F(N)· T F(N) -1 ·F(N)· T M S (N)

[0253] Where F(N) is the data stacking amount of the sensor force; M S (N) is the data stacking amount of the sensor torque;

[0254] Thus, the center of gravity coordinates (P x , P y , P z ) and intermediate variables k1, k2, k3;

[0255] By formula

[0256]

[0257] It turns out that:

[0258]

[0259]

[0260] Simplified to:

[0261]

[0262] In the formula, k4, k5, and k6 are intermediate variables;

[0263] Among them, the intermediate variables k4, k5, and k6 are defined as follows:

[0264]

[0265] Further simplified to:

[0266] F s =R·A

[0267] Where Fs and R represent the data stack of sensor force and robot joint angle obtained by N (N≥3) samplings of the robot end effector posture transformation; A represents the gravity of the end effector in the sensor coordinate system;

[0268] The calculation method of the gravity of the end effector in the sensor coordinate system is as follows:

[0269] A=(R(N) T R(N) -1 ·R(N) T F S (N)

[0270] Where, Fs and R represent the data stacking of sensor force and robot joint angle obtained by sampling N (N≥3) times of the robot end effector posture transformation; A is identified by the least squares method, and k4, k5, and k6 are calculated;

[0271] Thus k4, k5, and k6 are calculated;

[0272] 27) Calculate the gravity G of the end effector and the installation angle (U, V) of the robot base as:

[0273]

[0274]

[0275] Where G is the gravity G of the end tool after the robot body and the end tool are installed, and the center of gravity (P X ,P y ,P z ); If the installation angle (U, V) of the robot base is a constant value, online compensation of the gravity of the end effector and the sensor zero point is required;

[0276] 28) Perform online compensation for the gravity and zero point of the end effector; the online compensation divides the six-dimensional force sensor calibration stage into the sensor zero point calibration period 0~T1 and the end effector gravity and zero point compensation period T1~T2 for processing: first, at a certain sampling time t within the period 0~T1, the robot is controlled to move freely, and the joint variables obtained by the robot controller are [θ1(t)θ2(t)θ3(t)θ4(t)θ5(t)θ6(t)] T , R(t) is calculated in real time, and the force and torque data sampled by the six-dimensional force sensor are [F(t)M(t)] T , gravity compensation value of end effector [F G (t)M G (t)] T for:

[0277]

[0278]

[0279] In the formula, S F G (k) is the representation of the gravity compensation value of the end effector in the S coordinate system; S M G (S) is the representation of the gravity compensation value of the end effector in the S coordinate system; is the rotation matrix of coordinate system G relative to coordinate system S; S P G is the center of gravity coordinate of the end effector in the sensor coordinate system. ; is the rotation matrix of coordinate system G relative to coordinate system S; G F G is the gravity compensation value of the end effector; G M Gis the gravity compensation value of the end effector; F0(t) and M0(t) are the zero-point data of the sensor after each power-on, which are constant values ​​after a period of fluctuation; F(t) and M(t) are the force and torque data sampled by the six-dimensional force sensor; F G (t), M G (t) is the gravity compensation value of the end effector. The zero point data of the sensor is defined as [F0(t)M0(t)] T .

[0280] Define the zero point data of the sensor as [F0(t)M0(t)] T , then:

[0281]

[0282] In the formula, S F G (t) is the representation of the gravity compensation value of the end effector in the S coordinate system; S M G (t) is the representation of the gravity compensation value of the end effector in the S coordinate system; is the rotation matrix of coordinate system G relative to coordinate system S; S P G is the center of gravity coordinate of the end effector in the sensor coordinate system; is the rotation matrix of coordinate system G relative to coordinate system S;

[0283] At a certain time k in the end effector gravity and zero point compensation period T1~T2, the robot is controlled to move freely, and the joint variables of the robot are [θ1(k)θ2(k)θ3(k)θ4(k)θ5(k)θ6(k)] T , calculate R(t) in real time, and calculate the gravity compensation data of the end effector online as [F G (k)M G (k)] T , according to the sensor zero point is a constant value [F0 M0] T ,have:

[0284]

[0285] Where, F(t) is the force data sampled by the six-dimensional force sensor; M(t) is the torque data sampled by the six-dimensional force sensor;

[0286] According to the online compensation results, the external force / torque is accurately detected. The external force and torque detected by the six-dimensional force sensor at the k sampling time is defined as [F e (k)M e (k)] T :

[0287]

[0288] In the formula, F e (k) is the force compensation error of the six-dimensional force sensor; M e (k) is the torque compensation error of the six-dimensional force sensor; F(k) is the torque data sampled by the six-dimensional force sensor; M(s) is the torque data sampled by the six-dimensional force sensor; F G (k) is the gravity compensation value of the end effector; M G (k) is the gravity compensation value of the end effector; F0 is the constant value of the force at the zero point of the sensor; M0 is the constant value of the torque at the zero point of the sensor.

[0289] Step 3: Analyze the environmental interaction model of the robot polishing process of complex surfaces to facilitate the transformation of forces in various coordinate systems, as follows:

[0290] 31) Establish the impedance model as follows:

[0291]

[0292] Where ΔF is the force error and K is the stiffness X c is the end position X r is the reference trajectory, B is the damping, M is the inertia term, X c is the command position trajectory; X r is the reference position trajectory; is the derivative of position is the reference speed of the sixth axis of the robot; is the second derivative of position is the acceleration of the sixth axis of the robot;

[0293] 32) Assuming that the position control method is ideal, the command position trajectory X c =Measurement position trajectory x m , the position error caused by the contact force is expressed as:

[0294] E=X c -X r =X m -X e ;

[0295] In the formula, X r is the reference trajectory X m Command trajectory X e for environmental location;

[0296] Assuming the environment is rigid, use the stiffness K e The spring is represented by x m =x c ;

[0297] 33) The contact force between the robot and the environment is expressed as:

[0298] F e =K e (X c -X r )=K e (X m -X e );

[0299] In the formula, F e is the contact force; K e is the environmental stiffness;

[0300] 34) The stiffness and damping of the environmental parameters are estimated by least squares. In the environmental model, the environmental model in a specific direction is considered. The environmental model is a scalar expression, b e , k e ,x,f e B e , K e , X, F e Damping, stiffness, displacement and force in a certain direction in:

[0301]

[0302] In the formula, is the derivative of the position error; Δx is the position error;

[0303] The discrete changes are:

[0304]

[0305] Where T is the sampling time; z is a complex variable;

[0306] is equivalent to:

[0307]

[0308] Step 4: RLS identifies the initial parameters of the environment, as follows:

[0309] 41) The RLS algorithm gives θ(k) as follows:

[0310]

[0311]

[0312] In the formula, G(k) is the gain factor; λ is the forget factor, which is 0.9; p(k) is the parameter estimation covariance matrix, and the recursive calculation formula is:

[0313]

[0314] 42) When the change of θ(k) value is less than the set threshold, it is considered that θ(k) value is stable. After obtaining the θ(k) value, the calculation formula of a1 and a2 is used to obtain (b e , k e ) are estimated to be:

[0315]

[0316]

[0317] Step 5: Using the acquired data, a self-learning compliant control method for the robot arm based on online sensor calibration is proposed to correct the reference trajectory and achieve zero-error steady-state tracking, as follows:

[0318] 51) Perform Laplace transform on the impedance model as follows:

[0319]

[0320] Where ΔF is the force error; M is the mass; E is the position error; B is the damping; K is the stiffness;

[0321] 52) When it is a one-dimensional vector, due to the expected contact force f d and the reference position trajectory x r is a constant, so The following conditions are met:

[0322]

[0323] f e =k e (x m -x e )

[0324]

[0325]

[0326] 53) The steady-state error of the force tracking function is as follows:

[0327]

[0328] In the formula, e ss is the steady-state error; k e is the environmental stiffness; k is the stiffness; x e is the environmental position; x r is the reference trajectory; f r is the reference force;

[0329] 54) In order to achieve zero-error steady-state tracking, make:

[0330]

[0331] 55) Estimate command position trajectory x c :

[0332]

[0333] In the formula, is the estimated value of the environmental position; f d For the power of expectation; is the estimated value of the environmental stiffness;

[0334] 56) Calculate the actual contact force f e , as shown below:

[0335]

[0336] In the formula, x m For x m is the end position;

[0337] Define k x =k e ·x e ,

[0338] In the formula, is the estimated value of the actual contact force;

[0339]

[0340] 57) Define Q1 as:

[0341]

[0342] 58) The Lyapunov function is defined as:

[0343] V=Q1 T PQ1

[0344]

[0345] Where V is the Lyapunov function; P is a diagonal positive definite matrix; q1 and q2 are positive constants. When , it tends to be stable;

[0346] Compute the derivative of the Lyapunov function as follows:

[0347]

[0348] In the formula, is the derivative of the Lyapunov function;

[0349] 59) It is deduced that the derivative of the Lyapunov function must be less than or equal to zero, so we have for:

[0350]

[0351]

[0352]

[0353]

[0354] Finally, the gravity of the end effector is compensated through sensor calibration, and the Kalman filter is used for noise reduction to correct the reference trajectory, thus completing active compliant control.

[0355] Embodiment 2

[0356] This embodiment conducts experimental verification on a six-degree-of-freedom robotic arm.

[0357] When the robot end contacts the plane, the force tracking effect is as follows: Figure 3 shown.

[0358] When the robot end contacts the inclined surface, the force tracking effect is as follows: Figure 4 shown.

[0359] When the robot end contacts the surface, the force tracking effect is as follows: Figure 5 shown.

[0360] The above experimental results show that the robot arm self-learning compliant control method based on online sensor calibration of the present invention can effectively adapt to the environment of complex curved workpieces and has a good tracking effect on position and force.

[0361] Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. A self-learning compliant control method for a robotic arm based on online sensor calibration, characterized in that: The following steps are involved: Step 1: Use Kalman filter to denoise the data collected by the six-dimensional force sensor; Step 2: Use online gravity compensation to obtain actual data of contact force; Step 3: Analyze the environmental interaction model of the robot polishing process of complex surfaces; Step 4: Identify initial parameters of the environment by RLS; Step 5: Using the acquired data, a self-learning compliant control method for the robot arm based on online sensor calibration is proposed to correct the reference trajectory and achieve zero-error steady-state tracking.

2. The method for self-learning compliant control of a robot arm based on online sensor calibration according to claim 1 is characterized in that: The specific method of step one is as follows: 11) Select Kalman filtering to filter the data collected by the six-dimensional force sensor, as shown below: x k =A·x k-1 +B·u k +w k-1 z k =H·x k +v k Where A is the state transfer matrix; B is the control input matrix; u k is the control vector; H is the state observation matrix; w k-1 is the process noise; v k is the measured noise, which is Gaussian white noise; w k-1 and v k Both are normally distributed; k and k-1 are time steps; 12) The state is predicted as follows: In the formula, is the prior state estimate, estimated forward from time step k-1 to k; 13) Update the status as follows: In the formula, is a posterior state estimate; 14) Perform Kalman gain calculation as follows: Where K is the prediction error in the optimal state estimation divided by the measurement error; 15) Perform covariance update as follows: Where P k is the estimated mean square error; is the prior mean square error matrix; 16) Update the process noise covariance as follows: In the formula, is the prediction covariance matrix; Q is the process noise covariance matrix.

3. The method for self-learning compliant control of a robot arm based on online sensor calibration according to claim 1 is characterized in that: The specific method of step 2 is as follows: 21) Obtain the force and torque data measured by the six-dimensional force sensor installed at the end of the robot. The data consists of three parts, as follows: a) The sensor's own system error is that the reading under no-load condition is not zero. Assuming that the sensor's own system error is a constant, it is called "zero point". Because the sensor has a "zero drift" phenomenon, the value is not constant, which is expressed as: [F ox F oy F oz M ox M oy M ox ] T ; In the formula, F ox is the error force in the x direction; F oy is the error force in the y direction; F oz is the error force in the z direction; M ox is the error torque rotating around the x direction; M oy is the error torque rotating around the y direction; M ox is the error torque rotating around the z direction; b) The gravity of the end effector. The gravity of the end effector connected to the sensor acts on the gravity of the six-dimensional force sensor, which changes with the posture of the end effector and is expressed as: [F Gx F Gy F Gz M Gx M Gy M Gx ] T ; In the formula, F Gx is the gravity of the end effector in the x direction; F Gy is the gravity of the end effector in the y direction; F Gz is the gravity of the end effector in the z direction; M Gx is the torque of the end effector in the x direction; M Gy is the torque of the end effector in the y direction; M Gx is the torque of the end effector in the z direction; c) The external force on the robot end effector is the relative output of the force detected by the six-dimensional force sensor and is expressed as: [F ex F ey F ez M ex M ey M ex ] T ; In the formula, F ex is the force on the end effector in the x direction; F ey is the force on the end effector in the x direction; F ez is the force on the end effector in the y direction; M ex is the torque of the end effector in the z direction; M ey is the torque of the end effector in the y direction; M ez is the torque of the end effector in the z direction; 22) Obtain the force and torque detected by the six-dimensional force sensor as follows: In the formula, F x F is the force in the x direction detected by the six-axis force sensor; y F is the force in the y direction detected by the six-dimensional force sensor; z M is the force in the z direction detected by the six-dimensional force sensor; x M is the force in the x direction detected by the six-dimensional force sensor; y M is the force in the y direction detected by the six-dimensional force sensor; z The force in the z direction detected by the six-dimensional force sensor; 23) Obtain the installation posture of the robot; the installation posture of the robot represents the rotation matrix of the robot base coordinate system B relative to the world coordinate system G; the world coordinate system is defined as 0X o Y o Z0, let the Z0 axis direction be vertically upward, the robot base coordinate system is defined as 00X1Y1Z1, the origin of B and G coincide; the world coordinate system G is transformed to the base coordinate system B by rotating V around the Y0 axis first, and then rotating U around the current X0 axis. The posture transformation matrix is ​​expressed as: Where U is the angle of rotation around the current X0 axis to the base coordinate system B; V is the angle of rotation of the world coordinate system G around the Y0 axis; 24) Obtain the posture of the robot's end effector; the posture of the robot's end effector represents the rotation matrix of the robot's end effector coordinate system T relative to the base coordinate system B; the industrial robot's connecting rod length and offset are obtained by measurement, and the joint angle is directly obtained by the robot control system. Assuming that the posture of the robot's sixth joint is consistent with the posture of the end effector, the posture of the end effector is solved by the forward kinematics formula; 25) According to the composite transformation relationship of the conversion matrix, the rotation matrix of the sensor coordinate system S relative to the world coordinate system G is as follows: In the formula, is the rotation matrix of the base coordinate system B relative to the world coordinate system G; is the rotation matrix of coordinate system T relative to coordinate system B; is the rotation matrix of coordinate system S relative to coordinate system T; 26) According to the Cartesian transformation relationship of the robot static force, the transformation relationship between force and torque in the G and S coordinate systems is obtained as follows: In the formula, S F S and S M S are the force and torque components of the gravity of the end effector in the sensor coordinate system S, respectively; is the rotation matrix of coordinate system G relative to coordinate system S; G F G , G M G is the gravity G of the end effector expressed in the world coordinate system G; S P G × is the coordinate of the center of gravity of the end effector in the sensor coordinate system, as shown below: Where P z It represents the z-axis of the end effector in the sensor coordinate system; P x is the x-axis representation of the end effector in the sensor coordinate system; P y is the y-axis representation of the end effector in the sensor coordinate system; Simplified to: Further simplified to: S M S = S P G × S F G Combined Thus, the following formula is obtained: In the formula, k1, k2, and k3 are intermediate variables; The definitions of intermediate variables k1, k2, and k3 are as follows: The above formula is simplified to: M S =F·P Where M S , F represents the data stacking amount of sensor force and sensor torque obtained by sampling N (N≥3) times of posture transformation of the robot end effector; P is the coordinate of the center of gravity of the end effector; The calculation method of the center of gravity coordinate P of the end effector is as follows: P=(F(N)· T F(N)) -1 ·F(N)· T M S (N) Where F(N) is the data stacking amount of the sensor force; M S (N) is the data stacking amount of the sensor torque; Thus, the center of gravity coordinates (P x , P y , P z ) and intermediate variables k1, k2, k3; By formula It turns out that: Simplified to: In the formula, k4, k5, and k6 are intermediate variables; Among them, the intermediate variables k4, k5, and k6 are defined as follows: Further simplified to: F s =R·A Where Fs and R represent the data stack of sensor force and robot joint angle obtained by sampling N (N ≥ 3) times of the robot end effector posture transformation; A is the representation of the gravity of the end effector in the sensor coordinate system. The calculation method of the gravity of the end effector in the sensor coordinate system is as follows: A=(R(N) T R(N)) -1 ·R(N) T F S (N) Where, Fs and R represent the data stacking of sensor force and robot joint angle obtained by sampling N (N≥3) times of the robot end effector posture transformation; A is identified by the least squares method, and k4, k5, and k6 are calculated; 27) Calculate the gravity G of the end effector and the installation angle (U, V) of the robot base as: Where G is the gravity of the end tool after the robot body and the end tool are installed; the installation angle (U, V) of the robot base is a constant value; 28) Perform online compensation for the gravity and zero point of the end effector; the online compensation divides the six-dimensional force sensor calibration stage into the sensor zero point calibration period 0~T1 and the end effector gravity and zero point compensation period T1~T2 for processing: first, at a certain sampling time t within the period 0~T1, the robot is controlled to move freely, and the joint variables obtained by the robot controller are [θ1(t)θ2(t)θ3(t)θ4(t)θ5(t)θ6(t)] T , R(t) is calculated in real time, and the force and torque data sampled by the six-dimensional force sensor are [F(t)M(t)] T , gravity compensation value of end effector [F G (t)M G (t)] T for: In the formula, S F G (k) is the representation of the gravity compensation value of the end effector in the S coordinate system; S M G (S) is the representation of the gravity compensation value of the end effector in the S coordinate system; is the rotation matrix of coordinate system G relative to coordinate system S; S P G is the center of gravity coordinate of the end effector in the sensor coordinate system. ; is the rotation matrix of coordinate system G relative to coordinate system S; G F G is the gravity compensation value of the end effector; G M G is the gravity compensation value of the end effector; F0(t) and M0(t) are the zero-point data of the sensor after each power-on, which are constant values ​​after a period of fluctuation; F(t) and M(t) are the force and torque data sampled by the six-dimensional force sensor; F G (t), M G (t) is the gravity compensation value of the end effector. Define the zero point data of the sensor as [F0(t)M0(t)] T , then: In the formula, S F G (t) is the representation of the gravity compensation value of the end effector in the S coordinate system; S M G (t) is the representation of the gravity compensation value of the end effector in the S coordinate system; is the rotation matrix of coordinate system G relative to coordinate system S; S P G is the center of gravity coordinate of the end effector in the sensor coordinate system; is the rotation matrix of coordinate system G relative to coordinate system S; At a certain time k in the end effector gravity and zero point compensation period T1~T2, the robot is controlled to move freely, and the joint variables of the robot are [θ1(k)θ2(k)θ3(k)θ4(k)θ5(k)θ6(k)] T , calculate R(t) in real time, and calculate the gravity compensation data of the end effector online as [F G (k)M G (k)] T , according to the sensor zero point is a constant value [F0 M0] T ,have: Where, F(t) is the force data sampled by the six-dimensional force sensor; M(t) is the torque data sampled by the six-dimensional force sensor; According to the online compensation results, the external force / torque is accurately detected. The external force and torque detected by the six-dimensional force sensor at the k sampling time is defined as [F e (k)M e (k)] T : In the formula, F e (k) is the force compensation error of the six-dimensional force sensor; M e (k) is the torque compensation error of the six-dimensional force sensor; F(k) is the torque data sampled by the six-dimensional force sensor; M(s) is the torque data sampled by the six-dimensional force sensor; F G (k) is the gravity compensation value of the end effector; M G (k) is the gravity compensation value of the end effector; F0 is the constant value of the force at the zero point of the sensor; M0 is the constant value of the torque at the zero point of the sensor.

4. The method for self-learning compliant control of a robot arm based on online sensor calibration according to claim 1, characterized in that: The specific method of step three is as follows: 31) Establish the impedance model as follows: Where ΔF is the force error and K is the stiffness X c is the end position X r is the reference trajectory, B is the damping, M is the inertia term, X c is the command position trajectory; X r is the reference position trajectory; is the derivative of position is the reference speed of the sixth axis of the robot; is the second derivative of position is the acceleration of the sixth axis of the robot; 32) Assuming that the position control method is ideal, the command position trajectory X c =Measurement position trajectory x m , the position error caused by the contact force is expressed as: E=X c -X r =X m -X e ; Where, X r is the reference trajectory X m Command trajectory X e for environmental location; Assuming the environment is rigid, use the stiffness K e The spring is represented by x m =x c ; 33) The contact force between the robot and the environment is expressed as: F e =K e (X c -X r )=K e (X m -X e ); In the formula, F e is the contact force; K e is the environmental stiffness; 34) The stiffness and damping of the environmental parameters are estimated by least squares. In the environmental model, the environmental model in a specific direction is considered. The environmental model is a scalar expression, b e , k e ,x,f e B e , K e , X, F e Damping, stiffness, displacement and force in a certain direction in: In the formula, is the derivative of the position error; Δx is the position error; The discrete changes are: Where T is the sampling time; z is a complex variable; is equivalent to:

5. The method for self-learning compliant control of a robot arm based on online sensor calibration according to claim 1, characterized in that: The specific method of step 4 is as follows: 41) Obtain θ by RLS algorithm (k) as follows: In the formula, G(k) is the gain factor; λ is the forget factor, which is 0.9; p(k) is the parameter estimation covariance matrix, and the recursive calculation formula is: 42) When the change of θ(k) value is less than the set threshold, it is considered that θ(k) value is stable. After obtaining the θ(k) value, the calculation formula of a1 and a2 is used to obtain (b e , k e ) are estimated to be:

6. The method for self-learning compliant control of a robot arm based on online sensor calibration according to claim 1, characterized in that: The specific method of step five is as follows: 51) Perform Laplace transform on the impedance model as follows: Where ΔF is the force error; M is the mass; E is the position error; B is the damping; K is the stiffness; 52) When it is a one-dimensional vector, due to the expected contact force f d and the reference position trajectory x r is a constant, so The following conditions are met: f e =k e (x m -x e ) 53) The steady-state error of the force tracking function is as follows: In the formula, e ss is the steady-state error; k e is the environmental stiffness; k is the stiffness; x e is the environmental position; x r is the reference trajectory; f r is the reference force; 54) In order to achieve zero-error steady-state tracking, make: 55) Estimate command position trajectory x c : In the formula, is the estimated value of the environmental position; f d For the power of expectation; is the estimated value of the environmental stiffness; 56) Calculate the actual contact force f e , as shown below: In the formula, x m For x m is the end position; Define k x =k e ·x e , In the formula, is an estimate of the actual contact force; 57) Define Q1 as: 58) The Lyapunov function is defined as: Where V is the Lyapunov function; P is a diagonal positive definite matrix; q1 and q2 are positive constants. When, stable; Compute the derivative of the Lyapunov function as follows: In the formula, is the derivative of the Lyapunov function; 59) It is deduced that the derivative of the Lyapunov function must be less than or equal to zero, so we have for:

Images