Bio-inspired soft robots are robotic systems made of materials that have the similar moduli order of natural organisms () . Different from the conventional rigid-link robotic systems with impedance control  or stiffness modulation , compliance is an intrinsic property of soft robots and is achieved through the material mechanical property (for example, the softness of silicone) or morphological control . Soft robots enable safe interaction with external objects and exhibit capability to manipulate inside confined environments or achieve complicated tasks. These robots have been used to grasp objects , to adapt to environments , to estimate external shapes , to locomote on rough terrain , to regain motor skills in rehabilitation applications , to support heart function  and many more [11, 12, 13].
Soft robots are usually actuated in two forms, namely the variable length tendons and pneumatic/fluidic actuation 
. Pneumatically actuated soft robots offer certain advantages due to the benefit of high power-to-weight ratio and easy implementation. Pneumatic muscle actuators (PMA) or McKibben actuator, a simple one degree of freedom (DoF) pneumatic soft actuator fabricated by a rubber tube and fiber sleeves, was proposed for artificial limb research since 1950’s and commercialized in the 1980s by Bridgestone Company. The history of pneumatically actuated soft robotic system dates back to 1992 when Suzumori et al. developed a 3-DoF (pitch, yaw, and extension) microactuator . Thanks to the widespread availability of 3D printing technique and affordable fabrication materials, the design and fabrication of pneumatically actuated soft robotics has been continuously growing in the past 10 years. For instance, a pneumatically-actuated soft bending actuator can be fabricated by (1) fabricating the outer molds of the soft actuator, (2) pouring the silicone to the molds, (3) leaving the silicone to cure, and (4) demolding the actuator .
Forward and inverse kinematic modeling of pneumatic soft robot are required to understand the robot inherent characteristic and enable accurate closed-loop control. In practice, each segment of a soft robot can be regarded as a continuum section when it is powered. Several approaches have been reported to study soft or continuum robot modeling. Piecewise constant-curvature assumption has been widely used in different soft/continuum robot forward kinematics modelings, such as the D-H parameter approach , Frenet-Serret frames approach , exponential coordinates approach  (see the detailed continuum robot kinematic modeling in ). This assumption provides a simple method to obtain the mapping from robot configuration space to task space. However, the piecewise constant curve assumption is not valid for all the conditions (see Figure 4 of  or Figure 4 in this paper), especially when the shape of a soft or continuum robot is significantly changed due to the large internal joint space pressure input, requiring more accurate modeling and control. A more general approach is a modal representation proposed by Chirikjian et al. in 1994 , which described the kinematics of a hyper-redundant snake robot with both constant and variable curvature. Wang and Simaan [23, 24] adapted this modal representation for multi-backbone continuum robots and presented kinematic error propagation and geometric calibration to address equilibrium shape deviation. Similar to the rigid serial robot, the inverse kinematics of soft continuum robot is more challenging compared to its forward kinematics. Neppalli et al. developed a closed-form geometric approach by modeling each segment of a continuum robot as a spherical joint and rigid link and using the conventional analytical method to solve the inverse kinematics . The Jacobian-based approach is able to search the arc parameters to achieve the desired configuration [19, 26]. Garbin et al. recently proposed a design that mechanically couples a pneumatic parallel bellows actuator (PBA) with a multi-backbone continuum user interface (CUI), and presented a unified kinematics for both the PBA (soft robot) and the CUI (continuum) .
With the aforementioned kinematic modeling approaches, the soft robots could be accurately controlled in the free space environment. However, due to the compliance of the soft robot, an external wrench disturbance (e.g. force or torque) in a realistic environment could potentially lead to significant robot configuration variation, resulting in inaccurate control. To estimate the external wrench applied to the soft robot body, Ozel et al. proposed a soft robot curvature sensing technique by integrating a resistive flex sensor or a magnetic curvature sensor . Kramer et al. proposed a polydimethylsiloxane (PDMS) based curvature sensing technique in 2011 . However, the soft robot kinematic model variation caused by the external contact remains unsolved.
Inspired by the work done by Bajo [30, 31], this paper aims (i) to develop a new kinematic modeling approach for soft robot with generalized configurations and (ii) to enable accurate contact detection during motion with feasible sensing feedback. We will systematically study the kinematic modeling of soft robot with contact, contact detection and contact localization. The contact force magnitude estimation is out of the scope of current work. We use an 1-DoF pneumatic bending actuator  to experimentally validate the modeling method. This paper is a step toward our overarching goal of robust soft robot control method for actual applications when the soft robot is in contact with the external objects. The rest of the paper is arranged as follows. Section II details the research problems and the basic modeling assumptions. Section III describes the theoretical kinematic modeling and contact detection algorithm. Section IV presents the simulation results followed by Conclusions in section V.
Ii Problem Statement
Figure 1 shows the soft robot prototype, which consists of a pneumatic bellow soft actuator and a manual pneumatic pump. The actuator fabrication process with fiber reinforcement technique was detailed in . The cross-section of pneumatic bellow actuator is semicircle shape, which has shown to have the smallest bending resistance compared to the rectangular shape and circle shape. A strain-limiting layer was added to its radial base to achieve the uni-directional bending motion and to prevent the axial extending and contracting motion. The circumferential reinforcement wires prevent the pneumatic bellow from significant radial expansion at high joint space pressure input. The maximum bending angle is occurred when the input pressure is 25 psi.
In this paper, we aim to address the following three problems:
Problem 1: Derive the generalized kinematic model of the 1-DoF pneumatic soft bellow actuator. Find the direct kinematics and instantaneous kinematics that maps from robot actuation space to task space.
Problem 2: Given the soft bellow end effector position and velocity, determine whether the external contact is applied to the bellow or not.
Problem 3: Given the soft bellow end effector position and velocity, find the contact location along the backbone of the bellow if the external contact is applied to the robot.
When discussing these problems, we use the following modeling assumptions:
The pneumatic bellow bends in the plane.
The curve of pneumatic bellow and its first-order derivative are continuous.
The pneumatic bellow end effector position and orientation can be obtained from an extrinsic sensor (e.g. electromagnetic tracker).
The portion of a soft bellow from base to contact position maintains its shape after contact and the unconstrained portion continues to bend as if it was a shorter bellow.
Iii Materials and Methods
Iii-a Forward Kinematics: Actuator Space to Task Space
The kinematics of continuum robots can be described using two concatenated mappings according to [18, 23]. The first one is a mapping from actuator space (input air pressure in this paper) to configuration space (tangential angle along the bending curve at the arc length of and the given input pressure ). The second one is from the configuration space (tangential angle ) to task space (soft actuator tip position and orientation). With the piecewise constant-curvature assumption , the robot end effector position in the inertial frame can be written as . The end effector orientation can be achieved by rotating about the axis for degree (see Figure 2). Therefore, the homogeneous transformation from robot based to end effector can be written in equation 1
where is the arc length and is the curvature. However, this forward kinematics expression is problematic when the curvature is zero. As can be seen from figure 2, zero curvature indicates that the continuum robot is a straight line, while the last column of homogeneous transformation matrix is undefined.
To address the above expression singularity and model the soft robot when the piecewise constant-curvature assumption is insufficiently accurate, we use the modal approach to describe the pneumatic bending characteristics from the actuator space to configuration space. Once the soft pneumatic bellow is powered and at the static position, the geometrically exact configuration only varies with respect to the input pressure and is determined by the minimal energy solution. Experimental calibration is performed with the hardware platform in figure 1 to obtain a family of shapes under the specific pneumatic pressure input. The mapping from actuator space to configuration space can be approximated in the following equation :
where and are described in the following modal representations
The bending actuator backbone can be decomposed into discrete points along the strain limiting layer. Each point has a specific corresponding tangential angle at different input pressure. Then equation 2 can be expressed in the following matrix form:
where , is the number of discrete points along the bending actuator backbone and is the number of input pressure. and are written as follows
Equation 5 can be reshaped into the following form:
where is the Kronecker’s product, and
are the vectorized form ofand respectively, which can be written as:
Solving matrix could get the mapping from the soft actuator joint space to the configuration space. The mapping from configuration space to task space can be explained by using the definition of curvature (see figure 2), as described in the following equation:
where is the point position expressed in the inertial frame at the arc length of and is the corresponding tangential angle (figure 2). The forward kinematics that maps from actuator space and task space is obtained by substituting equation 2 into equation 11 and integrate the tangent vector along the arc length.
Iii-B Instantaneous Kinematics
This section describes the instantaneous kinematics modeling that maps the actuator speed to the end effector twist . The end effector twist consists of the linear velocity and angular velocity . However, since the actuator bends in the plane, the linear velocity in direction and angular velocity about and axis are zero. We seek to find the Jacobian that satisfies the following equation:
The Jacobian the relates the joint space motion rates
with the end effector twist can be derived by using the differential chain rule:
where and can be written as follows:
Combining equation 13, 14 and 15, we can get the Jacobian that relates the actuator space speed to the robot end effector twist (linear and angular velocities). Given the specific arc length , the end effector twist is the function of input pressure and pressure change rate . The Jacobian was also validated numerically by using the resolved rates method , as can be seen in the following equation:
where is the positive scaling factor.
Iii-C Contact Detection
This section describes the contact detection algorithms as initially discussed in . The joint force deviation (JFD) method calculates the external wrench value based on the deviation of generalized joint space forces with respect to the nominal model. It can be described in the following equation:
where is the gradient of elastic energy with respect to the configuration perturbation, and are the joint force and external wrench respectively. When there is no external wrench applied on the soft robot, then . Once the wrench applied on the robot (), which leads to the variation of the joint force if the pre-determined detection threshold is . Based on the deviation value, we can obtain the external wrench according to the equation below:
where is the deviation of joint force. In the real applications, the contact detection threshold is typically not equal to 0, which is caused by the friction force in the pull-wire continuum robots and the compressibility characteristic of pneumatic supplyin the current robot. To implement this algorithm, one needs to have a accurate joint space pressure sensor to get the generalized joint force feedback. However, the major limitation of this algorithm is that it does not provide the contact location .
The fixed centrode deviation (FCD) method calculates the deviation of theoretical and actual loci of the fixed centrode of end effector to detect the contact position. The loci of theoretical fixed centrode (or instantaneous screw axis, ISA) can be expressed in the following equation :
where and are the end effector linear velocity and angular velocity respectively, which is calculated from the instantaneous kinematics model described in section III-B.
is the skew-symmetric angular velocity matrix defined as. Based on the modeling assumption (3) in section II, the actual end effector position and orientation can be captured through a tracking system (e.g. EM tracker or motion camera). The loci of actual fixed centrode is expressed as following:
where the subscript in equation 20 indicates the sensor data. Once the deviation of fixed centrode occurs, there exists an external wrench along the pneumatic bellow soft robot. Similarly, the detectability threshold also exists in the FCD method and it is affected by the tracking system resolution.
Iii-D Contact Location Estimation
This section presents the contact location estimation algorithm by using the kinematics modeling and contact detection method described above. Once the contact occurs at a given location , the robot kinematic model can be decomposed into two sections (see figure 3. Note that, without future explanation, position and position are expressed in the unit of pixel, which is 0.2959 mm/pixel). Due to the inherent compliance characteristic of the soft robot, the soft robot base to contact position forms the constrained portion. The contact position to end effector forms the unconstrained portion, which continues to bend with the increase joint space pressure. The unconstrained portion kinematic model is similar to the model described in section III-A and III-B, which has a shorter arc length and different base position  (which is the contact location). The modified kinematic model can be written as:
where is the input pressure when contact occurs, is the contact location, and is the tangential angles along the arc when contact occurs.
Having this relationship, the pneumatic bellow actuator instantaneous kinematics can be obtained readily by changing the actual arc length.
where is the Jacobian of the unconstrained portion by substituting the arc length into equation 13, 14 and 15. Then the loci of the fixed center can be derived by combining equations 11, 12, 19, 21, 22, and 23.
The contact location is estimated by solving a least square optimization problem that minimizes the Euclidian distance between the actual loci of fixed centrode obtained from sensor data and theoretical loci of the fixed centrode obtained from modeling results. The objective function is given by
where and is the weighting matrix. The gradient of fixed centrode in terms of contact location needs to be derived to implement this algorithm. According to equation 19, the gradient can be written in the following form:
With the gradient calculated above, the least square optimization problem will be solved iteratively according to .
This section describes the results of pneumatic bellow calibration performance, forward/instantaneous kinematics, motion simulation after contact, and contact position estimation.
Iv-a Pneumatic Bellow Calibration
This section presents the pneumatic bellow calibration results. The experimental setup includes a soft bending bellow, manual pneumatic pump, a high resolution camera (C920 Webcam, Logitech, Swiss) and custom designed graph user interface (GUI). The manual pneumatic pump has the resolution of 1Psi. High resolution input pressure control can be achieved by closing the pressure control loop with a digital pressure sensor. The pneumatic bellow with different air pressure inputs was captured by the camera system and stored for post processing. In this experiment, the input pressure samples 0 Psi, 6 psi, 10 Psi, 15 Psi and 21 Psi respectively (Figure 4).
The strain limiting layer on the pneumatic bellow was manually segmented to represent the pneumatic bellow shape. We annotated ten points along the each backbone and find the corresponding positions and tangential angles. By solving equation 8, we can obtain the model-dependent matrix. The segmented curve and calibrated model can be seen in Figure 5. The result shows that the calibrated model closely matches the experimental results. The difference between the calibrated model and the actual shape increases with the increase of input pressure. However, the maximum error occurred at 21 Psi input pressure is less than 2.1 mm. Increasing the number of point could increase the calibration accuracy at the expense of increasing the computation load.
Iv-B Instantaneous Kinematics Evaluation
This section validates the instantaneous model of the pneumatic bellow actuator as presented in equation (16). Given the actuator initial position and desired position, the resolved rates algorithm is applied to find the corresponding joint space pressure input. Note that the instantaneous kinematics model is evaluated in the scenario of no contact occurs. Figure 6 shows the iterative process when the pneumatic bellow converges to the desired position. The pneumatic bellow actuator end effector space error is 0.0326 mm (0.0978 pixels) and the joint space error is 0.016 Psi, which validates the accuracy of the proposed instantaneous kinematics model.
Iv-C Contact Detection and Estimation
This section presents the contact detection and contact location estimation. In the experimental study, the contact occurs at . The input pressure was increased from 5 Psi to 20 Psi with the step size of 0.05 Psi. The loci of fixed centrode throughout the trajectory was calculated simultaneously. We also provided the robot motion without any contact for comparative study. The simulation results are shown in Figure 7, which illustrates that the ISA difference between these two scenarios increases steadily due to the growth of input pressure. This is mainly due to the significant end effector position difference as the pressure goes up, see equation (19). Maximum ISA difference is observed when the input pressure is 20 Psi.
We also simulate the relation between ISA difference vs. the contact location. In the simulation, the contact position was chosen at to 400 from the actuator base with 50 increment per step (Figure 8). The performance index is defined as follows:
where is the fixed centrode when robot is in contact and is the fixed centrode when robot is free of contact. The ISA difference serves the sensitivity performance for the FCD method and it rises as a result of the contact location moves further from the base. Intuitively, contacts can be detected easily when they are closer to end effector due to the large difference in both end effector twist and position when such contacts occur. In contrast, it is difficult to detect contacts when they are close to the base. The ISA difference equals to 0 when (indicates that the contact point is located at the fixed base of the soft bellow actuator), which also represents that the contact does not affect the robot shape and the ISA.
Next we present the contact position estimation result by using the method discussed in section III-D. The ISA with a given contact location () is obtained from the simulation results and used as the ground truth. The estimation process starts with an initial guess () and then calculates the ISA error as well as the . Using the nonlinear least square optimization method, the solution iteratively converge to the ground truth. Another initial guess () was also tested. Both cases studies were able to find the contact position with the acceptable accuracy (figure 9 and 10).
This paper presents the a mathematical framework for kinematic modeling, contact detection, and contact location estimation of a 1-DoF soft bellow actuator. A modal representation approach was used to derive the forward and instantaneous kinematics model. We presented the simulation study of the pneumatic bellow actuator with contact by using a modified kinematic model. The fixed centrode deviation was used as the performance index for contact detection. Contact position was estimated by solving a nonlinear least square optimization problem. Simulation results show that the proposed method could accurately estimate the contact location with sub-pixel error.
The results in this paper show for the first time how the mapping from joint space and task space can be generated through the modal approach and how the contact detection could be estimated through a least square optimization manner. This will inspire future development of complex soft robot modeling and contact control to accomplish the desired tasks.
-  D. Rus and M. T. Tolley, “Design, fabrication and control of soft robots,” Nature, vol. 521, no. 7553, p. 467, 2015.
-  S. Part, “Impedance control: An approach to manipulation,” Journal of dynamic systems, measurement, and control, vol. 107, p. 17, 1985.
-  N. Simaan and M. Shoham, “Geometric interpretation of the derivatives of parallel robots’ jacobian matrix with application to stiffness control,” Journal of Mechanical Design, vol. 125, no. 1, pp. 33–42, 2003.
-  C. Laschi and M. Cianchetti, “Soft robotics: new perspectives for robot bodyware and control,” Frontiers in bioengineering and biotechnology, vol. 2, p. 3, 2014.
-  F. Ilievski, A. D. Mazzeo, R. F. Shepherd, X. Chen, and G. M. Whitesides, “Soft robotics for chemists,” Angewandte Chemie, vol. 123, no. 8, pp. 1930–1935, 2011.
-  S. A. Morin, R. F. Shepherd, S. W. Kwok, A. A. Stokes, A. Nemiroski, and G. M. Whitesides, “Camouflage and display for soft machines,” Science, vol. 337, no. 6096, pp. 828–832, 2012.
-  K. C. Galloway, K. P. Becker, B. Phillips, J. Kirby, S. Licht, D. Tchernov, R. J. Wood, and D. F. Gruber, “Soft robotic grippers for biological sampling on deep reefs,” Soft Robotics, vol. 3, no. 1, pp. 23–33, 2016, pMID: 27625917. [Online]. Available: https://doi.org/10.1089/soro.2015.0019
-  M. T. Tolley, R. F. Shepherd, B. Mosadegh, K. C. Galloway, M. Wehner, M. Karpelson, R. J. Wood, and G. M. Whitesides, “A resilient, untethered soft robot,” Soft Robotics, vol. 1, no. 3, pp. 213–223, 2014.
-  P. Polygerinos, Z. Wang, K. C. Galloway, R. J. Wood, and C. J. Walsh, “Soft robotic glove for combined assistance and at-home rehabilitation,” Robotics and Autonomous Systems, vol. 73, pp. 135–143, 2015.
-  E. T. Roche, M. A. Horvath, I. Wamala, A. Alazmani, S.-E. Song, W. Whyte, Z. Machaidze, C. J. Payne, J. C. Weaver, G. Fishbein et al., “Soft robotic sleeve supports heart function,” Science translational medicine, vol. 9, no. 373, p. eaaf3925, 2017.
-  C. Majidi, “Soft robotics: a perspective—current trends and prospects for the future,” Soft Robotics, vol. 1, no. 1, pp. 5–11, 2014.
-  S. Kim, C. Laschi, and B. Trimmer, “Soft robotics: a bioinspired evolution in robotics,” Trends in biotechnology, vol. 31, no. 5, pp. 287–294, 2013.
-  D. Trivedi, C. D. Rahn, W. M. Kier, and I. D. Walker, “Soft robotics: Biological inspiration, state of the art, and future research,” Applied Bionics and Biomechanics, vol. 5, no. 3, pp. 99–117, 2008.
-  H. Schulte Jr, “The characteristics of the mckibben artificial muscle (1961) the application of external power in prosthetics and orthotics,” National Academy of Sciences-National Research Council, Washington DC, Appendix H, pp. 94–115.
-  K. Suzumori, S. Iikura, and H. Tanaka, “Applying a flexible microactuator to robotic mechanisms,” IEEE control systems, vol. 12, no. 1, pp. 21–27, 1992.
-  “softroboticstoolkit,” https://softroboticstoolkit.com/, accessed: 2019-06-10.
-  M. W. Hannan and I. D. Walker, “Kinematics and the implementation of an elephant’s trunk manipulator and other continuum style robots,” Journal of Field Robotics, vol. 20, no. 2, pp. 45–63, 2003.
-  R. J. Webster III and B. A. Jones, “Design and kinematic modeling of constant curvature continuum robots: A review,” The International Journal of Robotics Research, vol. 29, no. 13, pp. 1661–1683, 2010.
-  P. Sears and P. Dupont, “A steerable needle technology using curved concentric tubes,” in Intelligent Robots and Systems, 2006 IEEE/RSJ International Conference on. IEEE, 2006, pp. 2850–2856.
-  R. J. Webster III, J. S. Kim, N. J. Cowan, G. S. Chirikjian, and A. M. Okamura, “Nonholonomic modeling of needle steering,” The International Journal of Robotics Research, vol. 25, no. 5-6, pp. 509–525, 2006.
-  S. Wakimoto, K. Ogura, K. Suzumori, and Y. Nishioka, “Miniature soft hand with curling rubber pneumatic actuators,” in Robotics and Automation, 2009. ICRA’09. IEEE International Conference on. IEEE, 2009, pp. 556–561.
-  G. S. Chirikjian and J. W. Burdick, “A modal approach to hyper-redundant manipulator kinematics,” IEEE Transactions on Robotics and Automation, vol. 10, no. 3, pp. 343–354, 1994.
-  L. Wang and N. Simaan, “Geometric calibration of continuum robots: Joint space and equilibrium shape deviations,” IEEE Transactions on Robotics, vol. 35, no. 2, pp. 387–402, 2019.
-  ——, “Investigation of error propagation in multi-backbone continuum robots,” in Advances in Robot Kinematics. Springer, 2014, pp. 385–394.
-  S. Neppalli, M. A. Csencsits, B. A. Jones, and I. D. Walker, “Closed-form inverse kinematics for continuum manipulators,” Advanced Robotics, vol. 23, no. 15, pp. 2077–2091, 2009.
-  R. Webster, J. Swensen, J. Romano, and N. Cowan, “Closed-form differential kinematics for concentric-tube continuum robots with application to visual servoing,” in Experimental Robotics. Springer, 2009, pp. 485–494.
-  N. Garbin, L. Wang, J. H. Chandler, K. L. Obstein, N. Simaan, and P. Valdastri, “Dual-continuum design approach for intuitive and low-cost upper gastrointestinal endoscopy,” IEEE Transactions on Biomedical Engineering, vol. 66, no. 7, pp. 1963–1974, July 2019.
-  S. Ozel, E. H. Skorina, M. Luo, W. Tao, F. Chen, Y. Pan, and C. D. Onal, “A composite soft bending actuation module with integrated curvature sensing,” in Robotics and Automation (ICRA), 2016 IEEE International Conference on. IEEE, 2016, pp. 4963–4968.
-  R. K. Kramer, C. Majidi, R. Sahai, and R. J. Wood, “Soft curvature sensors for joint angle proprioception,” in Intelligent Robots and Systems (IROS), 2011 IEEE/RSJ International Conference on. IEEE, 2011, pp. 1919–1926.
-  A. Bajo and N. Simaan, “Finding lost wrenches: Using continuum robots for contact detection and estimation of contact location,” in Robotics and Automation (ICRA), 2010 IEEE International Conference on. IEEE, 2010, pp. 3666–3673.
-  ——, “Kinematics-based detection and localization of contacts along multisegment continuum robots,” IEEE Transactions on Robotics, vol. 28, no. 2, pp. 291–302, 2012.
-  Z. Wang, P. Polygerinos, J. Overvelde, K. Galloway, K. Bertoldi, and C. Walsh, “Interaction forces of soft fiber reinforced bending actuators,” IEEE/ASME Transactions on Mechatronics, 2016.
-  P. Polygerinos, Z. Wang, J. T. Overvelde, K. C. Galloway, R. J. Wood, K. Bertoldi, and C. J. Walsh, “Modeling of soft fiber-reinforced bending actuators,” IEEE Transactions on Robotics, vol. 31, no. 3, pp. 778–789, 2015.
-  B. A. Jones and I. D. Walker, “Kinematics for multisection continuum robots,” IEEE Transactions on Robotics, vol. 22, no. 1, pp. 43–55, 2006.
-  J. Zhang and N. Simaan, “Design of underactuated steerable electrode arrays for optimal insertions,” Journal of Mechanisms and Robotics, vol. 5, no. 1, p. 011008, 2013.
-  D. E. Whitney, “Resolved motion rate control of manipulators and human prostheses,” IEEE Transactions on man-machine systems, vol. 10, no. 2, pp. 47–53, 1969.
-  J. Angeles, Fundamentals of robotic mechanical systems: theory, methods, and algorithms. Springer Science & Business Media, 2013, vol. 124.
-  R. M. Murray, Z. Li, S. S. Sastry, and S. S. Sastry, A mathematical introduction to robotic manipulation. CRC press, 1994.
-  P. Lancaster and M. Tismenetsky, The theory of matrices: with applications. Elsevier, 1985.