377

Sensors and Materials, Vol. 35, No. 2 (2023) 377–389
MYU Tokyo
S & M 3174

Development and Simulation of Operational Scenarios
for Underwater Cable-driven Parallel Robot
Katutoshi Kodama,* Akihiro Morinaga, and Ikuo Yamamoto
Graduate School of Engineering, Nagasaki University, 1-14 Bunkyo, Nagasaki 852-8521, Japan(Received December 12, 2022; accepted February 1, 2023)

Keywords: cable-driven, parallel robot, underwater robot, operation scenario, feedback control

A cable-driven parallel robot (CDPR) is a mechanism that controls the position and attitude ofan object using multiple cables. We propose an underwater cable-driven parallel robot (UCDPR),
which is a CDPR composed of multiple surface robots. A UCDPR is a type of mobile cabledriven parallel robot (MCDPR) that is composed of multiple mobile robots and is an underwater
application of MCDPR, which has been used only on land and in the air. We describe the detailsof the operational scenario of the UCDPR, from landing on the water to executing a task. We
also simulated numerically the stabilization phase required after landing on water and the
trajectory tracking control phase required for underwater exploration and other tasks. In thisnumerical simulation, we used a high-gain feedback controller as the trajectory tracking
controller to add robustness to the control system. As a result, a trajectory tracking control was
realized within a tolerance tracking error range of 10−2 m underwater in the presence of a
maximum current velocity of 0.8 m/s (≅ 1.5 knots), which is a control requirement.

1.

Introduction

A cable-driven parallel robot (CDPR) is a type of parallel mechanism that uses multiple
cables to constrain an object and control its position and attitude by controlling the cable
lengths.(1,2) When lightweight and flexible cables are used, the mechanism is lightweight andcan achieve a large payload capacity and function in a large workspace.
A CDPR consisting of multiple mobile robots is called a mobile cable-driven parallel robot(MCDPR).(3) An MCDPR is different from a typical CDPR in that it has kinematic redundancy
in its mechanism because the cable feed position can be actively changed by the controller. By
appropriately using kinematic redundancy, an MCDPR can do sub-tasks such as obstacle
avoidance and singularity avoidance in addition to the main task of trajectory tracking control.
The MCDPR has been applied to FASTKIT(4) and MoPICK (Mobile Pick-and-Place),(5) which do
pick-and-place operations in factories, and to mechanisms that coordinate multiple aerial robots
to do transportation and rescue operations.(6,7) As described above, MCDPR has been considered
for many applications on land and in the air, but not much thought has been given to its
application underwater. One of the reasons for this is that the drag force and water currents on*Corresponding author: e-mail: bb52221102@ms.nagasaki-u.ac.jp

https://doi.org/10.18494/SAM4280

ISSN 0914-4935 © MYU K.K.
https://myukk.org/

378

Sensors and Materials, Vol. 35, No. 2 (2023)

the object underwater are large and deform the cable, making it difficult to control the positionof the suspended object.
To solve the problem of underwater application of MCDPR, we proposed an underwater
cable-driven parallel robot (UCDPR), which is composed of multiple robots that work on thesurface of the water [autonomous surface vehicles (ASVs), unmanned surface vehicles (USVs),
and others],(8) and developed a modeling method for the mechanism that takes into account cable
dynamics underwater(9,10) and for a basic control system configuration.(10,11)
To operate the UCDPR in an actual environment, it is necessary to consider operational
scenarios(12,13) from the landing of the robot on the water to the execution of the objective task
(or the recovery of the robot), just as one must consider with general underwater robots [such as
remotely operated vehicles (ROVs) and autonomous underwater vehicle (AUVs)], and to design
an appropriate controller for each scenario. Therefore, we have developed a scenario based on
the operation of a UCDPR in an actual environment and have performed control simulations
based on this operational scenario. To simplify the discussion, we assumed that the verticalmovement of the water surface due to waves on the water may be disregarded.

2.

Overview of a UCDPR

2.1

Mechanism

Figure 1 shows a conceptual diagram of our proposed UCDPR. This UCDPR is an underwater
application of an MCDPR using multiple surface robots. A UCDPR with the configuration
shown in Fig. 1 can control the position and attitude (six degrees of freedom) of a suspended
object in three-dimensional space.

Fig. 1.

(Color online) Image of a UCDPR.

Sensors and Materials, Vol. 35, No. 2 (2023)

379

In our UCDPR, the surface robot’s position is measured by the Global Positioning System(GPS), and the cable lengths are measured by a rotary encoder. For a typical CDPR, the position
and attitude of the suspended object can be easily estimated from the geometric relationships to
the CDPR. However, underwater, it is difficult to accurately estimate the position and attitude ofa suspended object because of the effect of currents. Therefore, the position of the suspended
object was measured by acoustic positioning, and its attitude was measured by a gyro sensor.
2.2

Experimental device

Figure 2 shows the exterior and system configuration of the experimental device. Figure 2(a)shows the experimental device, which is composed of two trolleys at the top of the frame. Two
cables are played out from each trolley to control the position of the suspended object on a twodimensional surface. The two trolleys are attached to independent timing belts, which move left
and right along the guide rails. These timing belts are actuated by DC motors. Winches areattached to the top of the trolleys to control the cable lengths. These winches are actuated by DC
servo motors. All motors have hall encoders to measure the displacement of the trolley and thereeling length of the cables from the angle of rotation of the motors.
The position of the suspended object is measured by acoustic positioning in actual
environments. However, the experimental device shown in Fig. 2(a) is too small to have an
acoustic positioning system installed. Therefore, image processing (color tracking) using
augmented reality (AR) markers attached to the four corners of the experimental device andcolor tape attached to the suspended object is used to obtain the position of the suspended object
in real time.
The experimental device is also attached to a commercially available water pump to simulate
water currents. Although analysis of the velocity distribution is currently underway, theexperimental device enables control experiments that simulate actual conditions at sea.
Figure 2(b) shows a diagram of the control system of this experimental device. The control
system is composed of a Python-based control system on a computer (host system) and a control
system on a control board (Arduino Mega) attached to the experimental device (low-ordersystem). The host system and the low-order system exchange control command values and(a)Fig. 2.

(b)(Color online) (a) Experimental device and (b) system configuration of the experimental device.

380

Sensors and Materials, Vol. 35, No. 2 (2023)

sensor values via serial communication. The host system (computer) acquires the position of thesuspended object by image processing of the video captured by the USB camera. Using theposition of the suspended object, the cable lengths, and the position of the trolley sent from theexperimental device, the computer calculates the manipulated variable and sends the pulse widthmodulation (PWM) value to the low-order system (Arduino Mega). The low-order system
controls the actuator according to this PWM value and sends the rotation angle of the motor tothe host system (computer). The control period of the experimental device is determined by theframes per second (FPS) of the USB camera and the time required for image processing. The
control period of the control system used was 0.1 s.

3.

Development of Operational Scenarios

The operational scenario assumed the basic operation of the UCDPR in an actual environment(especially in the sea). A UCDPR realizes a chain of motion control from UCDPR landing on thewater to task execution by switching controllers according to each scenario. Although Fig. 3
shows only the motion of the UCDPR on a two-dimensional surface for the simplification of the
explanation, the same argument can be made in the actual environment (three-dimensional
space).
3.1

Scenario A: Calibration of mechanism and positioning control

Because the UCDPR is a mechanism composed of multiple surface robots and a suspended
object, it was difficult to accurately pinpoint the landing site on water of all components when
landing the UCDPR from the ground or a mother ship, even if the landing position was plannedin advance. In addition, there was a concern that the suspended object may sway underwater
because of the currents. Therefore, in Scenario A, calibration of the initial positions of thesurface robots and other components and positioning control of the suspended object underwater
was carried out to stabilize the mechanism. At the same time, we also checked the sensors andcommunication system.

Fig. 3.

(Color online) Diagram of operational scenario.

Sensors and Materials, Vol. 35, No. 2 (2023)

3.2

381

Scenario B: Guidance and control

In Scenario B, the UCDPR was guided from its initial position, which was held in place by
Scenario A, to the target position (e.g., the start of the reference path in Scenario C or the target
position in the work area in Scenario D). Guidance control of the UCDPR was achieved by
calculating the guidance path (trajectory) that smoothly connected the initial position and thetarget position in real time and controlling UCDPR autonomously to track this guidance path.
3.3

Scenario C: Trajectory tracking control

In Scenario C, the trajectory tracking control of the UCDPR was carried out according to a
preplanned path (trajectory). The UCDPR trajectory tracking control system must be robust
against parameter errors, modeling errors, and variations in current velocity.
3.4

Scenario D: Execution of objective tasks

In Scenario D, the UCDPR carried out an arbitrary task at a preplanned target point by
positioning the surface robots and controlling the cable lengths. The objective task was done byremote control from the ground or mother ship, or by autonomous control of the UCDPR.
Concrete tasks were expected to include underwater fixed-point observation and installation ofunderwater structures.

4.

Numerical Simulation Based on Operational Scenarios

We numerically simulated a simple operational example composed of Scenario A andScenario C for the experimental device shown in Fig. 2(a).
4.1Mathematical models

First, the dynamics model of the entire UCDPR to be used as the control object of thenumerical simulation is given by Eq. (1). For a detailed derivation of this equation, refer toRef. 10.
 M (q ) −G T (q )   q   F   =  
0  λ   γ 
 G (q )(1)

Here, q is the state vector of the position and attitude of the entire UCDPR, M(q) is the mass
matrix, and F(q) is a vector of generalized forces. The terms G(q) and γ are defined by theholonomic constraint Φ(q) = 0 for the entire mechanism as382

Sensors and Materials, Vol. 35, No. 2 (2023)∂Φ=(q ) =q G (q )q ,
Φ
∂q

(2)

 (q ) G (q )q + G (q )q ,
Φ
=

(3)

γ =
−G (q )q − 2α G (q )q − β 2 Φ (q ) ,

(4)

and

where Φ(q) = 0 is a position constraint, Eq. (2) is a velocity constraint, and Eq. (3) is anacceleration constraint. Equation (4) is a Baumgarte stabilization term that reduces numericalerrors and has α > 0 and β ≥ 0.
Next, we describe the inverse kinematics problem around a suspended object of UCDPR,
which is necessary for the controller design. The inverse kinematics problem for UCDPR isdifferent from that for CDPR because it requires solving the inverse dynamics problem in thecalculation. Therefore, the inverse kinematics problem of UCDPR is an inverse kinematics
problem in the broadest sense. For a detailed discussion of this problem, refer to Ref. 10. The
solution of the inverse kinematics around the suspended object is given byκ = J obj (qobj , λobj )qobj  J obj (qobj , λobj )u

(5)

+T
=
λobj Gobj( M obj qobj − Fobj ) + hobj kobj ,

(6)

and

where qobj is the state vector of the position and attitude of the suspended object, and κ is avector of the position of the trolleys and the cable lengths, which are manipulated variables of the
UCDPR. The term u is a virtual control input in the kinematics controller, and by giving anappropriate feedback input, the manipulated variables κ can be calculated while reducing theeffects of disturbances. The term Jobj is a matrix that expresses the velocity relationship between
qobj and κ, and is represented by the manipulation force (tension) λobj on the suspended object,which is the inverse dynamics solution of Eq. (6). In Eq. (6), Mobj and Fobj are the mass matrix
and the vector of the generalized force of the suspended object. hobj kobj is the internal force
caused by the kinematic redundancy of the UCDPR, and hobj satisfies Gobj hobj = 0. The term kobj
is an arbitrary vector expressing kinematic redundancy, which can be used to subtask themechanism if appropriately determined by solving optimization problems related to the tension
distribution of the UCDPR, trolley placement, and other features.(11)Sensors and Materials, Vol. 35, No. 2 (2023)

4.2

383

Control system

The trajectory tracking control system of the UCDPR is configured as shown in Fig. 4. In Fig.
ref
ref
ref
, qobjand qobj
are the reference trajectories of the suspended object. These reference4, qobj
trajectories are defined in Sect. 4.2.1 together with the UCDPR control objectives. About thevirtual input u, the positioning controller for Scenario A and the trajectory tracking controller
for Scenario C are defined in Sect. 4.2.2. The optimization problem of calculating the kinematic
redundancy parameter kobj in Eq. (6) is formulated in Sect. 4.2.3.
4.2.1 Control objectiveFirst, we describe the design of the reference trajectory of the suspended object. As shown ina related study,(14) considering the actual control of trolleys and winches, if the manipulated
variable κ and its first-order derivative (i.e., velocity) κκ are not continuous functions, thesuspended object will vibrate. As shown in Eq. (5), κ includes the inverse kinematic solution λobj
of Eq. (6). Therefore, qobj should be continuous until jerked, i.e., a third-order differentiable. In
this study, the reference path and reference trajectory of the suspended object were designed
independently, using the parameter s(t) (0 ≤ s ≤ 1), which explicitly includes the time t (0 ≤ t ≤ T).
init
and the coordinates of the
Denoting the coordinates of the starting position of the path as qobj
fin
ending position as qobj , the reference trajectory of the suspended object is obtained using s(t) as
follows:
ref
init
fin
qobj
+ s (t )qobj( s(t ) ) =(1 − s(t ) ) qobj
.

(7)

init
fin
and qobj
and isThe path of the suspended object shown by Eq. (7) is the line connecting qobj
clearly a third-order differentiable. At the same time, s(t) should also be a third-orderdifferentiable, so s(t) is given here as a seventh-order polynomial. Eight boundary conditions are

Fig. 4.

(Color online) Block diagram of the trajectory tracking system.

384

Sensors and Materials, Vol. 35, No. 2 (2023)

needed to determine the coefficients of this polynomial. For these boundary conditions, s(0) = 0,
init
fin
at the start point of the path, and s(T) = 1, so that qobj (T ) = qobj
so Eq. (7) becomes qobj ( 0 ) = qobj
at the endpoint of the path. Also, s=( 0 ) s=(T ) s=( 0 ) s=(T ) s=( 0 ) s=(T ) 0 , so that thesuspended object stops at the start and end of the path (i.e., the velocity, acceleration, and jerk of
the suspended object are all set to 0). On the basis of this analysis, the path parameter s(t) is
given by
4

5

6

7

t 
t 
t 
t 
s (t ) = 35   − 84   + 70   − 20   .
T 
T 
T 
T 

(8)

On the basis of this discussion, the reference trajectory of the suspended object is given in
T
T
init
fin
= [ 0.10, 0.30] and qobj
= [ 0.40, 0.15] . The 0 to 10 s phase is the positioning
Fig. 5, where qobj
control phase of Scenario A, the 10 to 20 s phase is the trajectory tracking control phase of
Scenario C, and the 20 to 30 s phase is for observing the residual vibration of the suspended
object.
Next, we describe the control objectives of the numerical simulations. We designed thecontrol system to stabilize the mechanism within 10 s after landing on water (Scenario A) and to
track the reference trajectory of the suspended object shown in Fig. 5 within ± 10−2 m of position
error (Scenario C) in the sea area with a maximum current velocity of 0.8 m/s ( ≅ 1.5 knots).
4.2.2 Design of the controllers
The feedback controller in the following equation is the positioning controller used in
Scenario A:

Fig. 5.

(Color online) Reference trajectory of the suspended object.

Sensors and Materials, Vol. 35, No. 2 (2023)

ref=uA k A (qobj− qobj ) .

385

(9)

Here, kA is the feedback gain. The feed-forward + feedback controller in the following equationis the trajectory tracking controller used in Scenario C:
ref
ref
uC =
qobj
+ kC (qobj− q) .

(10)

Here, the first term on the right-hand side is the feedforward controller for trajectory tracking,and the second term is the feedback controller for reducing disturbances and modeling errors.
The feedback gain kC is given a high gain to reduce the effect of current and modeling errors
caused by the dynamics of the cable. In general, this is called a high-gain feedback control (15)and is a simple method to improve the robustness of the control system. In this study, kC is
adjusted by trial and error for the motion of the UCDPR at the maximum current velocity of
0.8 m/s specified in the control objectives, and its effectiveness is verified in Sect. 4.3.
4.2.3 Optimization problem for redundancy parameter
The optimizer is designed to give the redundancy parameter kobj, which indicates the
kinematic redundancy. We considered the optimization problem of the tension distribution ofcables, which is often considered a subtask of the typical MCDPR.(3,11) The optimization
problem of tension distribution in UCDPR is defined asT
min V = λobj
λobj
kobj

+T
s.t.
=
λobj Gobj( M obj qobj − Fobj ) + hobj kobjκ = J obj (qobj , λobj )uκ = ∫ κ dt(11)

0 ≤ xvi ≤ 0.7 (i =1, 2)0.10 ≤ xv 2 − xv1 .

In this nonlinear optimization problem, xvi is the displacement of the trolley, and theinequality constraints are the range of motion of the two trolleys and the requirement thatthe two trolleys do not collide.
4.3

Simulation

Numerical simulations were performed using the control system described in Sect. 4.2. In
this numerical simulation, the simulation environment was constructed using MATLAB/
Simulink.

386

Sensors and Materials, Vol. 35, No. 2 (2023)

4.3.1 Variation of trajectory trackability in the case of uniform flow
Figure 6 shows the simulation results when the UCDPR tracked the reference trajectory
shown in Fig. 5 underwater at a current velocity of 0.8 m/s. Figures 6(a)–6(c) show the results ofthe positioning control of Scenario A, and Figs. 6(d)–6(f) show the results of the trajectory
tracking control of Scenario C. As shown in Fig. 6(a), the cable lengths and the position of thetrolleys, which are determined from the geometric relationship of the UCDPR in static water, aregiven as appropriate initial values for the manipulated variables κ in the numerical simulation. In
Figs. 6(a) and 6(b), the UCDPR swung from the initial position of the path due to the current, but
by executing Scenario A, the state of the mechanism was stabilized as shown in Fig. 6(c) while

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 6. (Color online) Simulation results for scenario A and scenario C (current velocity: −0.8 m/s). (a) t = 0 s,(b) t = 1.0 s, (c) t = 3.0 s, (d) t = 10.0 s, (e) t = 15.0 s, and (f) t = 20.0 s.

Sensors and Materials, Vol. 35, No. 2 (2023)

387

the effect of the current was reduced. In Figs. 6(d)–6(f), the UCDPR shifts to Scenario C andaccurately tracks the preplanned reference path.
Figure 7 shows the variations of the tracking error for the numerical simulations for several
current velocities. In Scenario A, the effect of the current velocity on the UCDPR was able to be
reduced within a predefined time for velocities up to 0.8 m/s. After shifting to Scenario C, theposition of the suspended object was able to be accurately tracked to the reference trajectory
within the tolerance range. The maximum tracking error was observed around 15 s for all
current velocities when the velocity of the reference trajectory was at its maximum.
Although the maximum current velocity was set at 0.8 m/s in the control objective in Sect. 4.1,
the effect of the current velocity was reduced in a given time even when a disturbance of 1.2 m/s,
which is outside the design range, was applied, and the tracking error was also kept within thetolerance range.
The numerical simulations in this study did not include considerations of constraints such as
the maximum reeling speed of the winches and the speed of the trolley. In the future, it will benecessary to develop controllers that consider these constraints of the mechanism.
4.3.2 Variation of trajectory trackability with variation of flow velocity
Figure 8 shows the variation in the behavior of a suspended object and the tracking error
when the current velocity was changed stepwise from 0.7 to 0.8 m/s (t = 13 s). As shown in Fig.
8, the suspended object deviated strongly from the reference path due to the step change incurrent velocity during Scenario C and then converged at the end of the reference path withdamped oscillations caused by the recovery motion. However, it takes a considerable amount oftime for the suspended object to recover to the reference path.
These results indicate that a high-gain feedback controller improves the stability and the
trajectory tracking performance of the control system but deteriorates other control features thathave a trade-off relationship with the stability. In particular, the gain tuning has been done bytrial and error for the motion of the UCDPR at a maximum current velocity of 0.8 m/s, which isdefined as the control objective. Therefore, excessive robustness may have been given to thecontrol system, resulting in a conservative system.

Fig. 7.

(Color online) Comparisons of tracking errors.

388

Sensors and Materials, Vol. 35, No. 2 (2023)(a)Fig. 8.

(b)(Color online) Comparison of tracking error: (a) t = 0 s and (b) t = 1.0 s.

In general, the control parameter of high-gain feedback is reported to depend on unknown
parameters such as the modeling error of the control model.(16) Excessive robustness was shownto cause degradation of control performance, so the high gain was developed into an adaptive
control system that performs an online search based on control errors and other variables. In thisway, it may be possible to achieve both optimal stability and control performance at all times
while suppressing the effects of unknown parameters.

5.

Conclusions

In this paper, we described an operational scenario for a UCDPR composed of multiple
surface robots and assuming an actual environment. We simulated the control of the UCDPR inan operation scenario that included the control phase of positioning after landing on water(Scenario A) and the control phase of tracking along a predefined path (Scenario C).
Because the main purpose of this paper was to develop an operational scenario, a high-gain
feedback controller was applied in the numerical simulation of Scenario C as a robust controllerthat is relatively easy to implement. As a result, it was found that while sufficient trajectory
tracking was achieved for different current velocities in a uniform current, trajectory tracking
deteriorated when the current velocity changed during the trajectory tracking.
The control system design and simulation in this study focused mainly on the current velocity
as the only factor that changes the control performance of the UCDPR. However, in actual
environments, other factors may have a strong effect on the performance, such as the effect ofwaves on surface robots, the stiffness of cables, and parameter errors of suspended objects.
The trajectory tracking controller in Eq. (10) is based on the kinematics around the suspended
object. Therefore, even if the trajectory tracking of the suspended object is sufficient for a givenreference path, there is no guarantee that it will be sufficient for a different reference trajectory
because the inertia and other dynamic characteristics of the suspended object will change. In the
future, we will develop an online retuning method of feedback gain based on adaptive control tosolve these problems as well as that of the conservativeness of robust control and to verify amore active robust control method.

Sensors and Materials, Vol. 35, No. 2 (2023)

389

In addition, we discussed only control methods for Scenarios A and C. In the future, we willdevelop theories for online trajectory generation and guidance control methods for Scenario B
and for determining the workspace for Scenario D.

References
1 H. Kino and S. Kawamura: J. Robot. Mechatron. 27 (2015) 599. https://doi.org/10.20965/jrm.2015.p0599
2 H. Tachiya: Parallel Mechnism (Morikita-Syuppan, Tokyo, 2019) 1st ed., Chap. 8 (in Japanese).
3 Collaborative Mobile Cable-Driven Parallel Robots–Archive ouverte HAL: https://theses.hal.science/tel02938097 (accessed January 2023).
4 N. Pedemonte, T. Rasheed, D. Marquez-Gamez, P. Long, E. Hocquard, F. Babin, C. Fouche, G. Caverot, A.
Girin, and S. Caro: Advances in Robotics Research: From Lab to Market (Springer, 2020). https://doi.
org/10.1007/978-3-030-22327-4_8
5 T. Rasheed, P. Long, D. Marquez-Gamez, and S. Caro: Proc. 4th Int. Conf. Cable-Driven Parallel Robots(Springer, 2019) 257. https://doi.org/10.1007/978-3-030-20751-9_22
6 N. Michael, J. Fink, and V. Kumar: Auton. Rob. 30 (2011) 73. https://doi.org/10.1007/s10514-010-9205-0
7 G. A. Cardona, D. Tellez-Castro, and E. Mojica-Nava: IFAC-PapersOnLine 52 (2019) 145. https://doi.
org/10.1016/j.ifacol.2019.12.149
8 K. Kodama, A. Morinaga, and I. Yamamoto: Proc. SICE SI-2021 (SICE, 2021) 1488 (in Japanese).
9 K. Kodama, A. Morinaga, and I. Yamamoto: Proc. Robomech2022 (JSME, 2022) 2A1-B07 (in Japanese).
10 K. Kodama, A. Morinaga, and I. Yamamoto: J. Jpn. Soc. Nav. Archit. Ocean Eng. (JASNAOE) 36 (2022) 53 (inJapanese).
11 K. Kodama, A. Morinaga, and I. Yamamoto: Proc. SICE SI-2022 (2022) 929 (in Japanese).
12 M. Nakamura, T. Hyakudome, H. Yoshida, and T. Aoki: J. Jpn. Soc. Nav. Archit. Ocean Eng. 9 (2009) 139 (inJapanese). https://doi.org/10.2534/jjasnaoe.9.139
13 M. Nakamura, S. Ishibashi, T. Hyakudome, H. Yoshida, and T. Aoki: J. Jpn. Soc. Nav. Archit. Ocean Eng. 14(2011) 143 (in Japanese). https://doi.org/10.2534/jjasnaoe.14.143
14 M. Yamamoto, N. Yanai, and A. Mohri: JRSJ 17 (1999) 1037 (in Japanese). https://doi.org/10.7210/jrsj.17.1037
15 H. Maeda: J. SICE 26 (1987) 381 (in Japanese). https://doi.org/10.11499/sicejl1962.26.381
16 S. Masuda: J. SICE 40 (2001) 717 (in Japanese). https://doi.org/10.11499/sicejl1962.40.717