# Analytical Synthesis of Regulators for Nonlinear Systems with a Terminal State Method on Examples of Motion Control of a Wheeled Robot and a Vessel

Authors: E. A. Shushlyapin; A. E. Bezuglaya
Journal: Journal of Applied Mathematics
 (2018)

---

## Abstract

The paper is devoted to several examples of control algorithm development for two-wheeled double-track robot and low-tonnage vessel-catamaran with two Azipods that show practical aspects of the application of one nonlinear system control method — terminal state method. This method, developed by the authors of the present paper, belongs to the class of methods for inverse dynamics problem solving. Mathematical models of control objects in the form of normal systems of third-order nonlinear differential equations for the wheeled robot and seventh-order ones for the vessel are presented. Design formulas of the method in general form for terminal and stabilizing controls are shown. A routine of obtaining calculation expressions for control actions is shown. Results of computer simulation of bringing the robot to a given point in a given time, as well as bringing the vessel to a given course during a “strong” maneuver, are described.


---

## Body

## 1. Introduction

Currently, ultimate nonlinear system modes often arise in control problems. In such cases, classical methods based on the linear automatic control theory are unusable. Due to it, more and more attention is being paid to the development and improvement of nonlinear system control methods. We can single out several approaches that gave rise to most nonlinear system control methods. They include methods of solving the inverse mechanics (dynamics) problem [1–4]; differential-geometric methods [5]; and methods based on artificial neural networks [6] and fuzzy logic [7]. These approaches do not require preliminary linearization in the classical sense. Therefore, there is no need to specify reference modes for the system operation and subsequently pursue a composite control. One of the methods of the first group is a “terminal state method” (TSM) [4]. Both the control objective and the system model are formulated via the “terminal state variables” (TSV) in it. These variables depend on two times (terminal t
f and current t) and on the current state x
t. The TSV sense is the finite state of the uncontrolled system which is at the current moment t in the state x
t. The ability to solve control problems via the TSV is related to its obvious property x
-
t
,
t
,
x
t
=
x
t, which is valid for any moments t
∈
t
0
,
t
f of the system operation interval. The terminal state method is developed for various types of nonlinear systems: differential systems with smooth right sides; discrete in time and mixed; with state independent variable time delays; multi-index (discrete analogs of the distributed systems). The initial application area of the TSM was terminal problems on a finite time interval. Consecutive application of the TSM at short intervals allows solving adaptive control problems, including problems with reference models. In general, the calculation of TSM control requires significant computational costs. They are associated with the multiple numerical integration of differential equations to evaluate the TSV and nonlinear (according to V.M. Alekseev) transition matrices as functions of the second argument (in the continuous case) or solutions of nonlinear systems of finite type equations (in discrete cases). There are, however, problems in which the control evaluation is performed by sufficiently simple finite type expressions that are easily realized in practice.The paper is devoted to the questions of practical application of the method under discussion in cases when TSM-control calculations are rather simple. So are the terminal problem of bringing a two-wheeled double-track mobile robot to a given point and the problem of bringing the vessel at a given heading angle. Herewith, it is required to ensure the possibility of pursuing the so-called “strong” vessel maneuvers. These are, for instance, maneuvers with a large angle of deviation of the steering components (in this paper, Azipod).Two-wheeled two-track robots are employed as amphibious platforms for observation and measurement. For example, a similar robot with crawler is proposed in [8] to be used, in particular, for monitoring underwater cables.A mathematical model of a two-wheeled double-track robot is given in paper [9], and an algorithm for its terminal control based on the TSM is constructed in [4]. Because of the robot mathematical model specificity, obtained algorithm is a set of algebraic relations. Control was implemented on real equipment — programmable robot Sparki by ArcBotics. It showed the algorithm efficiency in conditions when discrete time step in the control calculation is not less than 0.1s. Such step, combined with the algorithm relative complexity and the need to evaluate robot current position by the numbering method using its mathematical model, prompts the use of controls calculated at discrete moments [10]. Herewith, the current state can be calculated with a small step, and the control with a large step.The second problem being solved is a part of a larger problem, which is to develop a crew-free operation mode for a small-sized robotic research vessel designed for environmental monitoring of the marine coastal area. In this case, a version of TSM with stabilizing (tracking) type criterion was used. To keep the vessel on a previously specified route, it is necessary to control it in automatic mode so that in case of deviation from the route the vessel returns to the route. In this case, there may be situations (sudden obstacle avoidance, moving in narrow areas, etc.) when the vessel maneuver is performed at a large rudder angle or Azimuthing Electric Propulsion Drive (Azipod) rotation. Then vessel dynamics becomes nonlinear and therefore the appropriate control methods should be applied. In navigation such vessel maneuvers are called “strong”.The specifics of the problem to keep the vessel on a route and the vessel mathematical model limit the range of possible control methods appropriate for obtaining a solution. In particular, it refers to the nonlinear vessel model we have obtained [11]. To solve the problem of vessel control during “strong” maneuvers, when model linearization is difficult, it is proposed, for example, to use methods [12, 13]. The essence of these inventions is the following. Firstly, control is pursued according to information from satellites and onboard measuring devices. Secondly, it is calculated with the use of “engineering” control signals generating methods, i.e., intuitively understandable algorithms that do not require vessel mathematical model using. Thus, for example, rudder deflection angle is calculated by numerical integration of the following differential equation [12]:(1)
d
δ
d
t
=
k
1
K
-
K
∗
+
k
2
ω
-
k
3
δ
,where δ
,
K
,
K
∗
,
ω are, respectively, rudder angle, current course, given course, angular velocity of the ship’s rotation around the center of mass; k
1
,
k
2
,
k
3 are undetermined coefficients.For bow maneuvering devices, it is proposed to use additionally the algorithm in the form of differential equation relative to the transverse control force [13]:(2)
d
P
d
t
=
k
1
β
+
k
2
∫
β
d
t
-
k
3
P
.Herewith algorithms of calculating anticipatory points are used. From these points precalculated program controls are pursued to turn the vessel to a new course. Conditions of the stabilization along a given courseK
∗ mode reverse activation are used as well.A common drawback of these methods, as well as others based on the use of linear feedback from the state coordinates [14, 15], is that they contain indeterminate coefficients that are unknown in advance and may depend on the operating conditions. To adapt to the operation conditions, principles of adaptive and intelligent control are used [16, 17]. To a lesser extent, such drawbacks are inherent in methods based on the use of control object mathematical models [18, 19].Often during the nonlinear system synthesis, methods of solving the inverse dynamics problem are used. These include, in particular, the input-output linearization method [1, 2], Batenko’s method [3], and the terminal state method (TSM) [4]. Method [1, 2] is developed for solving the stabilization or tracking problems and is suitable for keeping the vessel on a given route. The other two methods [3, 4] are terminal; i.e., they are designed to bring the system to a specified state at a given time. Their application to the vessel track-keeping control is possible when the system operation time is divided into small intervals, and then the terminal method is applied on each of these small intervals.


## 2. Terminal State Method

### 2.1. Terminal State Method for Terminal Differential Systems

Terminal state method was developed to solve terminal control problems for nonlinear affine systems in the form of (3) with smooth criterion functions and dynamic constraints of various kinds (differential, delay differential, finite difference, etc.).For differential systems with smooth constraints, criterion problem is formulated according to (3), (4).(3)
J
=
J
z
t
f
=
J
z
-
t
f
,
t
f
,
z
t
f
→
J
∗
,
(4)
d
z
d
t
=
Φ
t
,
z
+
B
t
,
z
u
,
t
∈
t
0
,
t
f
,
z
t
0
=
z
0
.According to this method, control (TSM control) is calculated by the following expressions:(5)
u
t
,
z
t
)
=
G
+
·
f
J
,
G
+
=
G
T
G
G
T
-
1
,
G
=
∂
J
z
∂
z
z
=
z
-
t
f
,
t
,
z
t
·
W
t
f
,
t
,
z
t
B
t
,
z
t
,
f
J
=
J
∗
-
J
z
-
t
f
,
t
,
z
t
T
uwhere f
J is the right side of the differential equation d
J
/
d
t
=
(
J
∗
-
J
)
/
T
u, defining a criterion function exponential approximation with time constant T
u to the desired value J
∗, which occurs at a sufficiently small value of T
u
<
<
t
f. Criterion function J
z
-
t
f
,
t
,
z
t is a function of time t. Other forms for the criterion function curve may be specified. Thus, with a variable during a period T
u
=
t
f
-
t shape of the line, along which the criterion function converges to J
∗, is no longer an exponent, but a straight line.Obtained controls may be constrained by amplitude limitations, if necessary. Such constraints in the active state affect deviation from the given (e.g., exponential) form of the criterion function. At the same time, control goal—to achieve a given terminal state—is commonly accomplished. Sometimes, for the process of control to be sustainable, artificial imposition of restrictions is preferable.In (5), terminal state variables (TSV) z
- and nonlinear transition matrix according to V.M. Alekseev (TM) W are used. They are defined by a joint system of differential equations [20]:(6)
d
z
-
ϑ
,
t
,
z
t
d
ϑ
=
Φ
ϑ
,
z
-
,
d
W
ϑ
,
t
,
z
t
d
ϑ
=
∂
Φ
ϑ
,
z
∂
z
z
=
z
-
ϑ
,
t
,
z
t
·
W
,
ϑ
∈
t
,
t
f
,
z
-
t
,
t
,
z
t
=
z
t
,
W
t
,
t
,
z
t
=
I
,where I is identity matrix.It is important to realize that for calculations in (5), TSV and TM are required as functions of the second argument t, and they are defined as functions of the first argument ϑ for each of the t values by (6).


### 2.2. Terminal State Method for Differential System Stabilization Tasks

For stabilization problems, the algorithm formed from (5), (6) is used at each of the short time intervals t
i
+
1
-
t
i
=
h
=
c
o
n
s
t
.Then the criterion transforms to(7)
J
=
J
z
t
i
+
1
=
J
z
-
t
i
+
1
,
t
i
+
1
,
z
t
i
+
1
→
J
∗
,and instead of (6), simplified expressions obtained by linearizing (6) are used:(8)
z
-
t
i
+
1
,
t
i
,
z
t
i
=
z
t
i
+
Φ
t
i
,
z
t
i
·
h
,
W
t
i
+
1
,
t
i
,
z
t
i
=
I
+
∂
Φ
t
i
,
z
∂
z
z
=
z
t
i
·
h
,
i
=
1,2
,
…
,
N
.


## 2.1. Terminal State Method for Terminal Differential Systems

Terminal state method was developed to solve terminal control problems for nonlinear affine systems in the form of (3) with smooth criterion functions and dynamic constraints of various kinds (differential, delay differential, finite difference, etc.).For differential systems with smooth constraints, criterion problem is formulated according to (3), (4).(3)
J
=
J
z
t
f
=
J
z
-
t
f
,
t
f
,
z
t
f
→
J
∗
,
(4)
d
z
d
t
=
Φ
t
,
z
+
B
t
,
z
u
,
t
∈
t
0
,
t
f
,
z
t
0
=
z
0
.According to this method, control (TSM control) is calculated by the following expressions:(5)
u
t
,
z
t
)
=
G
+
·
f
J
,
G
+
=
G
T
G
G
T
-
1
,
G
=
∂
J
z
∂
z
z
=
z
-
t
f
,
t
,
z
t
·
W
t
f
,
t
,
z
t
B
t
,
z
t
,
f
J
=
J
∗
-
J
z
-
t
f
,
t
,
z
t
T
uwhere f
J is the right side of the differential equation d
J
/
d
t
=
(
J
∗
-
J
)
/
T
u, defining a criterion function exponential approximation with time constant T
u to the desired value J
∗, which occurs at a sufficiently small value of T
u
<
<
t
f. Criterion function J
z
-
t
f
,
t
,
z
t is a function of time t. Other forms for the criterion function curve may be specified. Thus, with a variable during a period T
u
=
t
f
-
t shape of the line, along which the criterion function converges to J
∗, is no longer an exponent, but a straight line.Obtained controls may be constrained by amplitude limitations, if necessary. Such constraints in the active state affect deviation from the given (e.g., exponential) form of the criterion function. At the same time, control goal—to achieve a given terminal state—is commonly accomplished. Sometimes, for the process of control to be sustainable, artificial imposition of restrictions is preferable.In (5), terminal state variables (TSV) z
- and nonlinear transition matrix according to V.M. Alekseev (TM) W are used. They are defined by a joint system of differential equations [20]:(6)
d
z
-
ϑ
,
t
,
z
t
d
ϑ
=
Φ
ϑ
,
z
-
,
d
W
ϑ
,
t
,
z
t
d
ϑ
=
∂
Φ
ϑ
,
z
∂
z
z
=
z
-
ϑ
,
t
,
z
t
·
W
,
ϑ
∈
t
,
t
f
,
z
-
t
,
t
,
z
t
=
z
t
,
W
t
,
t
,
z
t
=
I
,where I is identity matrix.It is important to realize that for calculations in (5), TSV and TM are required as functions of the second argument t, and they are defined as functions of the first argument ϑ for each of the t values by (6).


## 2.2. Terminal State Method for Differential System Stabilization Tasks

For stabilization problems, the algorithm formed from (5), (6) is used at each of the short time intervals t
i
+
1
-
t
i
=
h
=
c
o
n
s
t
.Then the criterion transforms to(7)
J
=
J
z
t
i
+
1
=
J
z
-
t
i
+
1
,
t
i
+
1
,
z
t
i
+
1
→
J
∗
,and instead of (6), simplified expressions obtained by linearizing (6) are used:(8)
z
-
t
i
+
1
,
t
i
,
z
t
i
=
z
t
i
+
Φ
t
i
,
z
t
i
·
h
,
W
t
i
+
1
,
t
i
,
z
t
i
=
I
+
∂
Φ
t
i
,
z
∂
z
z
=
z
t
i
·
h
,
i
=
1,2
,
…
,
N
.


## 3. Model Formulation

### 3.1. A Mathematical Model of a Two-Wheeled Double-Track Mobile Robot

The robot model excluding the wheel drives dynamics description, i.e., its kinematic model, has the following form [10]:(9)
d
x
1
d
t
=
D
w
ω
1
+
ω
2
4
cos
⁡
ϕ
,
d
x
2
d
t
=
D
w
ω
1
+
ω
2
4
sin
⁡
ϕ
,
d
ϕ
d
t
=
D
w
ω
1
-
ω
2
2
L
r
,where x
1
,
x
2 are linear movements of the robot along the coordinate axes; D
w is the diameter of the first and second wheels; ω
1
,
ω
2 are the first and second wheels angular speeds of rotation (control actions); L
r is the distance between wheels.The kinematic model defined by (9) actually coincides with the model from [16]. The only difference between the models is that control actions in (9) have the form of robot wheels angular velocities, and in [16] the form of robot wheels linear velocities V
L
,
V
R. Replacing linear velocities with expressions V
L
=
(
D
w
/
2
)
ω
1, V
R
=
(
D
w
/
2
)
ω
2 makes the model from [16] coincide with (9).First of all, it should be noted that models of this type take place in the case of using stepper driving motors, as well as when using crawler propulsions.Some features of this model are as follows: the right sides of its differential equations are independent from the state vector; controlsω
1
,
ω
2 are included linearly in the right sides with state dependent coefficients. The latter feature allows us to classify this model as affine (see below).


### 3.2. A Mathematical Model of a Vessel

#### 3.2.1. A Mathematical Model of a Vessel for Integrating Drives

Mathematical model was obtained in [11] by the methods described in [21–24] assuming that the Azipod drives are of integrating type. In accordance with the problem formulation, vessel is controlled by rotating two Azipod units at the same angle δ with given thrusts equal on both units: T
v
1, T
v
2, i.e., T
v
1
=
T
v
2
=
T
v.Mathematical model has the form of (4), where the vessel state vector consists of seven components:(10)
z
=
c
o
l
o
n
V
x
,
V
y
,
ω
,
x
1
,
y
1
,
K
,
δ
.Herewith,δ
,
K
,
ω are, respectively, Azipod rotation angle, heading angle, vessel angular velocity relative to its center of mass; V
x
,
V
y
,
x
1
,
x
2 are, respectively, speeds and positions of the vessel center of mass along its longitudinal and transverse axes. The seventh component δ
≡
z
7 that is a control action essentially, is included in vessel model equations not additively (see below). To transform the system to an affine form, its state vector is extended by one component z
7
≡
δ with a differential equation for the new state:(11)
d
δ
d
t
=
uwhere u is the new control action. Equation (11) describes a drive of integrating type. The transformed system of the seventh order, in accordance with the form (4) and extended state vector, has the following components.(12)
Φ
t
,
z
≔
2
·
T
v
·
cos
⁡
z
7
-
R
x
z
-
m
y
·
z
3
·
z
2
m
x
2
·
T
v
·
sin
⁡
z
7
-
m
x
·
z
3
·
z
1
m
y
L
·
T
v
·
sin
⁡
z
7
+
M
r
e
z
z
+
M
D
z
-
z
1
·
z
2
·
m
x
-
m
y
J
w
V
z
·
cos
⁡
z
6
-
β
z
V
z
·
sin
⁡
z
6
-
β
z
z
3
0
(13)
B
t
,
z
≔
c
o
l
o
n
0
0
0
0
0
0
1State functionsV
,
R
x
,
R
y
,
M
K
,
M
D
,
β are defined in terms of variables z by the expressions [21]:(14)
V
=
z
1
2
+
z
2
2
,
β
=
arctg
⁡
z
2
z
1
,
(15)
R
x
V
,
z
1
,
z
2
=
2000
·
R
K
·
V
·
cos
⁡
β
,
(16)
R
y
V
,
z
1
,
z
2
=
2000
·
R
K
·
V
·
sin
⁡
β
,
(17)
R
K
V
=
8
,
1
-
8
,
436
·
V
+
2
,
536
·
V
2
,
(18)
M
K
V
,
z
1
,
z
2
=
R
y
V
,
z
1
,
z
2
·
l
K
,
l
K
=
0
,
4
·
L
·
0.5
-
β
π
,
(19)
M
D
V
,
z
3
=
-
0
,
148
ρ
2
V
·
24
,
092
·
L
2
z
3
.In (12)-(19), V is the vessel resulting speed in the coordinate system associated with the vessel; R
x
,
R
y are the vessel resistance forces in the longitudinal and transverse directions considering the two catamaran hulls and conversion of the force dimension from kilonewtons to newtons; R
K is approximated by the polynomial vessel resistance in kilonewtons dependent on the speed V (in m
·
s
-
1), obtained by simulating vessel in the OpenFOAM software; L is the vessel length in meters; m
x
,
m
y are the vessel weights, taking into account water added masses in longitudinal and transverse movements; ρ
=
1025
k
g
·
m
-
3 is the sea water density.Masses and moment of inertia from (12) are determined by the relations [21]:(20)
m
x
=
1
+
k
11
D
,
m
y
=
1
+
k
22
D
,
J
ω
=
1
+
k
66
J
Z
,where J
Z is the moment of inertia of the vessel; D is the weight draught (mass) of the vessel in kilograms. Coefficients of the added masses and moment are selected from the following ranges: k
11
=
0.05
÷
0.15
;
k
22
=
0.7
÷
1.2
;
k
66
=
0.6
÷
1.3. For the considered vessel of small tonnage D
=
155000
k
g, the lower bounds of the mentioned coefficients were taken.


#### 3.2.2. A Mathematical Model of a Vessel Considering Drives Inertance

To obtain control formula considering Azipod drives dynamics instead of (11), we should use the following drive dynamics equation in the form of first order inertial link of (4):(21)
d
z
7
d
t
=
-
1
T
δ
z
7
+
k
δ
T
δ
vwith a coefficient k
δ, time constant T
δ, and a control action v.Equating right sides of (21) with the control u for the case when drives dynamics is not taken into account, results in(22)
-
1
T
δ
z
7
+
k
δ
T
δ
v
=
u
.Subsequently we obtain expression demonstrating dependence between the controls with and without assumption of the drive dynamics:(23)
v
=
1
k
δ
z
7
+
T
δ
k
δ
u
.


## 3.1. A Mathematical Model of a Two-Wheeled Double-Track Mobile Robot

The robot model excluding the wheel drives dynamics description, i.e., its kinematic model, has the following form [10]:(9)
d
x
1
d
t
=
D
w
ω
1
+
ω
2
4
cos
⁡
ϕ
,
d
x
2
d
t
=
D
w
ω
1
+
ω
2
4
sin
⁡
ϕ
,
d
ϕ
d
t
=
D
w
ω
1
-
ω
2
2
L
r
,where x
1
,
x
2 are linear movements of the robot along the coordinate axes; D
w is the diameter of the first and second wheels; ω
1
,
ω
2 are the first and second wheels angular speeds of rotation (control actions); L
r is the distance between wheels.The kinematic model defined by (9) actually coincides with the model from [16]. The only difference between the models is that control actions in (9) have the form of robot wheels angular velocities, and in [16] the form of robot wheels linear velocities V
L
,
V
R. Replacing linear velocities with expressions V
L
=
(
D
w
/
2
)
ω
1, V
R
=
(
D
w
/
2
)
ω
2 makes the model from [16] coincide with (9).First of all, it should be noted that models of this type take place in the case of using stepper driving motors, as well as when using crawler propulsions.Some features of this model are as follows: the right sides of its differential equations are independent from the state vector; controlsω
1
,
ω
2 are included linearly in the right sides with state dependent coefficients. The latter feature allows us to classify this model as affine (see below).


## 3.2. A Mathematical Model of a Vessel

### 3.2.1. A Mathematical Model of a Vessel for Integrating Drives

Mathematical model was obtained in [11] by the methods described in [21–24] assuming that the Azipod drives are of integrating type. In accordance with the problem formulation, vessel is controlled by rotating two Azipod units at the same angle δ with given thrusts equal on both units: T
v
1, T
v
2, i.e., T
v
1
=
T
v
2
=
T
v.Mathematical model has the form of (4), where the vessel state vector consists of seven components:(10)
z
=
c
o
l
o
n
V
x
,
V
y
,
ω
,
x
1
,
y
1
,
K
,
δ
.Herewith,δ
,
K
,
ω are, respectively, Azipod rotation angle, heading angle, vessel angular velocity relative to its center of mass; V
x
,
V
y
,
x
1
,
x
2 are, respectively, speeds and positions of the vessel center of mass along its longitudinal and transverse axes. The seventh component δ
≡
z
7 that is a control action essentially, is included in vessel model equations not additively (see below). To transform the system to an affine form, its state vector is extended by one component z
7
≡
δ with a differential equation for the new state:(11)
d
δ
d
t
=
uwhere u is the new control action. Equation (11) describes a drive of integrating type. The transformed system of the seventh order, in accordance with the form (4) and extended state vector, has the following components.(12)
Φ
t
,
z
≔
2
·
T
v
·
cos
⁡
z
7
-
R
x
z
-
m
y
·
z
3
·
z
2
m
x
2
·
T
v
·
sin
⁡
z
7
-
m
x
·
z
3
·
z
1
m
y
L
·
T
v
·
sin
⁡
z
7
+
M
r
e
z
z
+
M
D
z
-
z
1
·
z
2
·
m
x
-
m
y
J
w
V
z
·
cos
⁡
z
6
-
β
z
V
z
·
sin
⁡
z
6
-
β
z
z
3
0
(13)
B
t
,
z
≔
c
o
l
o
n
0
0
0
0
0
0
1State functionsV
,
R
x
,
R
y
,
M
K
,
M
D
,
β are defined in terms of variables z by the expressions [21]:(14)
V
=
z
1
2
+
z
2
2
,
β
=
arctg
⁡
z
2
z
1
,
(15)
R
x
V
,
z
1
,
z
2
=
2000
·
R
K
·
V
·
cos
⁡
β
,
(16)
R
y
V
,
z
1
,
z
2
=
2000
·
R
K
·
V
·
sin
⁡
β
,
(17)
R
K
V
=
8
,
1
-
8
,
436
·
V
+
2
,
536
·
V
2
,
(18)
M
K
V
,
z
1
,
z
2
=
R
y
V
,
z
1
,
z
2
·
l
K
,
l
K
=
0
,
4
·
L
·
0.5
-
β
π
,
(19)
M
D
V
,
z
3
=
-
0
,
148
ρ
2
V
·
24
,
092
·
L
2
z
3
.In (12)-(19), V is the vessel resulting speed in the coordinate system associated with the vessel; R
x
,
R
y are the vessel resistance forces in the longitudinal and transverse directions considering the two catamaran hulls and conversion of the force dimension from kilonewtons to newtons; R
K is approximated by the polynomial vessel resistance in kilonewtons dependent on the speed V (in m
·
s
-
1), obtained by simulating vessel in the OpenFOAM software; L is the vessel length in meters; m
x
,
m
y are the vessel weights, taking into account water added masses in longitudinal and transverse movements; ρ
=
1025
k
g
·
m
-
3 is the sea water density.Masses and moment of inertia from (12) are determined by the relations [21]:(20)
m
x
=
1
+
k
11
D
,
m
y
=
1
+
k
22
D
,
J
ω
=
1
+
k
66
J
Z
,where J
Z is the moment of inertia of the vessel; D is the weight draught (mass) of the vessel in kilograms. Coefficients of the added masses and moment are selected from the following ranges: k
11
=
0.05
÷
0.15
;
k
22
=
0.7
÷
1.2
;
k
66
=
0.6
÷
1.3. For the considered vessel of small tonnage D
=
155000
k
g, the lower bounds of the mentioned coefficients were taken.


### 3.2.2. A Mathematical Model of a Vessel Considering Drives Inertance

To obtain control formula considering Azipod drives dynamics instead of (11), we should use the following drive dynamics equation in the form of first order inertial link of (4):(21)
d
z
7
d
t
=
-
1
T
δ
z
7
+
k
δ
T
δ
vwith a coefficient k
δ, time constant T
δ, and a control action v.Equating right sides of (21) with the control u for the case when drives dynamics is not taken into account, results in(22)
-
1
T
δ
z
7
+
k
δ
T
δ
v
=
u
.Subsequently we obtain expression demonstrating dependence between the controls with and without assumption of the drive dynamics:(23)
v
=
1
k
δ
z
7
+
T
δ
k
δ
u
.


## 3.2.1. A Mathematical Model of a Vessel for Integrating Drives

Mathematical model was obtained in [11] by the methods described in [21–24] assuming that the Azipod drives are of integrating type. In accordance with the problem formulation, vessel is controlled by rotating two Azipod units at the same angle δ with given thrusts equal on both units: T
v
1, T
v
2, i.e., T
v
1
=
T
v
2
=
T
v.Mathematical model has the form of (4), where the vessel state vector consists of seven components:(10)
z
=
c
o
l
o
n
V
x
,
V
y
,
ω
,
x
1
,
y
1
,
K
,
δ
.Herewith,δ
,
K
,
ω are, respectively, Azipod rotation angle, heading angle, vessel angular velocity relative to its center of mass; V
x
,
V
y
,
x
1
,
x
2 are, respectively, speeds and positions of the vessel center of mass along its longitudinal and transverse axes. The seventh component δ
≡
z
7 that is a control action essentially, is included in vessel model equations not additively (see below). To transform the system to an affine form, its state vector is extended by one component z
7
≡
δ with a differential equation for the new state:(11)
d
δ
d
t
=
uwhere u is the new control action. Equation (11) describes a drive of integrating type. The transformed system of the seventh order, in accordance with the form (4) and extended state vector, has the following components.(12)
Φ
t
,
z
≔
2
·
T
v
·
cos
⁡
z
7
-
R
x
z
-
m
y
·
z
3
·
z
2
m
x
2
·
T
v
·
sin
⁡
z
7
-
m
x
·
z
3
·
z
1
m
y
L
·
T
v
·
sin
⁡
z
7
+
M
r
e
z
z
+
M
D
z
-
z
1
·
z
2
·
m
x
-
m
y
J
w
V
z
·
cos
⁡
z
6
-
β
z
V
z
·
sin
⁡
z
6
-
β
z
z
3
0
(13)
B
t
,
z
≔
c
o
l
o
n
0
0
0
0
0
0
1State functionsV
,
R
x
,
R
y
,
M
K
,
M
D
,
β are defined in terms of variables z by the expressions [21]:(14)
V
=
z
1
2
+
z
2
2
,
β
=
arctg
⁡
z
2
z
1
,
(15)
R
x
V
,
z
1
,
z
2
=
2000
·
R
K
·
V
·
cos
⁡
β
,
(16)
R
y
V
,
z
1
,
z
2
=
2000
·
R
K
·
V
·
sin
⁡
β
,
(17)
R
K
V
=
8
,
1
-
8
,
436
·
V
+
2
,
536
·
V
2
,
(18)
M
K
V
,
z
1
,
z
2
=
R
y
V
,
z
1
,
z
2
·
l
K
,
l
K
=
0
,
4
·
L
·
0.5
-
β
π
,
(19)
M
D
V
,
z
3
=
-
0
,
148
ρ
2
V
·
24
,
092
·
L
2
z
3
.In (12)-(19), V is the vessel resulting speed in the coordinate system associated with the vessel; R
x
,
R
y are the vessel resistance forces in the longitudinal and transverse directions considering the two catamaran hulls and conversion of the force dimension from kilonewtons to newtons; R
K is approximated by the polynomial vessel resistance in kilonewtons dependent on the speed V (in m
·
s
-
1), obtained by simulating vessel in the OpenFOAM software; L is the vessel length in meters; m
x
,
m
y are the vessel weights, taking into account water added masses in longitudinal and transverse movements; ρ
=
1025
k
g
·
m
-
3 is the sea water density.Masses and moment of inertia from (12) are determined by the relations [21]:(20)
m
x
=
1
+
k
11
D
,
m
y
=
1
+
k
22
D
,
J
ω
=
1
+
k
66
J
Z
,where J
Z is the moment of inertia of the vessel; D is the weight draught (mass) of the vessel in kilograms. Coefficients of the added masses and moment are selected from the following ranges: k
11
=
0.05
÷
0.15
;
k
22
=
0.7
÷
1.2
;
k
66
=
0.6
÷
1.3. For the considered vessel of small tonnage D
=
155000
k
g, the lower bounds of the mentioned coefficients were taken.


## 3.2.2. A Mathematical Model of a Vessel Considering Drives Inertance

To obtain control formula considering Azipod drives dynamics instead of (11), we should use the following drive dynamics equation in the form of first order inertial link of (4):(21)
d
z
7
d
t
=
-
1
T
δ
z
7
+
k
δ
T
δ
vwith a coefficient k
δ, time constant T
δ, and a control action v.Equating right sides of (21) with the control u for the case when drives dynamics is not taken into account, results in(22)
-
1
T
δ
z
7
+
k
δ
T
δ
v
=
u
.Subsequently we obtain expression demonstrating dependence between the controls with and without assumption of the drive dynamics:(23)
v
=
1
k
δ
z
7
+
T
δ
k
δ
u
.


## 4. Terminal State Method Applications

### 4.1. Robot Control

A. Isidori’s linearizing by input/output method [1] in Wei Wu’s modification [2] and TSM [4] were used in [9] to bring robot to a given point at a given time. Because of the specificity of the model described by (4), both methods led to the same algorithm:(24)
u
t
=
G
+
J
∗
-
J
z
t
T
,
G
+
=
G
T
G
G
T
-
1
,
G
=
∂
J
z
t
∂
x
B
t
,
z
t
,where J
t
≡
J
z
t is the criterion function that defines desired exponential trajectory with a time constant T
u. Differential equation for this function is the following:(25)
d
J
t
d
t
=
J
∗
-
J
t
T
u
.(24), (25) result from (5), (6) and the fact that TM W in this case is identity matrix on the entire period of the system’s operation, since the uncontrolled system vector Φ is equal to zero according to (9).In (24) state variables z
=
z
i
,
i
=
1
,
…
,
n and control variables u
=
u
j
,
j
=
1
,
…
,
r, are redefined according to(26)
z
=
x
1
x
2
φ
,
u
=
ω
1
ω
2
.For the criterion function, the following quadratic form is used:(27)
J
z
t
f
=
F
1
z
1
t
f
-
x
1
∗
2
+
F
2
z
2
t
f
-
x
2
∗
2
+
F
3
z
3
t
f
-
ω
∗
2
→
J
∗
=
0
,where z
∗
=
c
o
l
o
n
x
1
∗
,
x
2
∗
,
ω
∗ are the desired values of the robot coordinates at the finite time t
f.The sense of criterion by (27) is that when the current state z
t converges to a given point z
∗, the value of the criterion function tends to the given value J
∗, which must obviously be equal to zero. This target may be unattainable by the control (24) that results in deviation of the actual behavior by (27) from the desired trajectory (25). The main reason for unattainability is that the scalar G
G
T nulls at some points or intervals of the time axis. Obviously, this fact causes unattainability, as control target does not meet the way it is being achieved along the trajectory according to (25). Another possible reason for the deviation of the actual trajectory from the given one is the amplitude constraints on the control actions that cannot be taken into account in the methods used. Thus, rotational speeds of the robot wheels are limited in amplitude by the ω
m value:(28)
ω
i
=
-
ω
m
,
u
i
<
-
ω
m
,
u
i
,
-
ω
m
≤
u
i
≤
ω
m
,
ω
m
,
u
i
>
-
ω
m
,
i
=
1,2
.The output of the controls evaluated by (24) with constraints by (28) also distorts the form of the criterion function defined by (25). But still the criterion function decreases, converging to some minimally attainable value J
∗. Herewith attained value J
∗ determines the accuracy of bringing state vector components to given values.System of (9) is a nonlinear affine system of (4) with two control actions ω
1
,
ω
2, where Φ is zero vector, and the components of the matrix B
=
B
i
j
,
i
=
1,2
,
3
,
j
=
1,2, are the following:(29)
B
11
=
B
12
=
D
w
cos
⁡
ϕ
4
,
B
21
=
B
22
=
D
w
sin
⁡
ϕ
4
,
B
31
=
-
B
32
=
D
w
2
L
r
.In [9], the results of Mathcad modeling of the system of (9) with control by (24), criterion by (27), and amplitude constraints by (28) are applied to the ArcBotics Sparki robot with the following parameters: D
w = 0,05 m; L
r = 0.085 m; F
1 = 10; F
2 = 0; F
3 = 1; x
1
∗ = 0.6 m; x
2
∗ = 0.3 m; x
3
∗ = 0.8 m; ω
m = 1.5s
-
1; h = 0.1 s; Tu = 0.01 s; t
f = 40 s.The same algorithm is implemented in the onboard C++ version of a real robot. As follows from the kinematic scheme, which is the basis for (9), there is a nonlinear dependence between the states of the system of (9).(30)
tg
⁡
z
3
=
z
2
z
1Consequently, equations (9) are dependent. It means that during control by (24) implementation, specified values can be reached simultaneously for a pair of variables x
1
,
φ
z
1
,
z
3 or for a pair x
2
,
φ
z
2
,
z
3. As for the pair x
1
,
x
2
z
1
,
z
2, this target is generally unattainable.This is clear from (9) with other control actions:(31)
v
1
=
ω
1
+
ω
2
2
,
v
2
=
ω
1
-
ω
2
2
.The resulting system is(32)
d
z
1
d
t
=
D
w
2
v
1
cos
⁡
ϕ
,
d
z
2
d
t
=
D
w
2
v
1
sin
⁡
ϕ
,
d
z
3
d
t
=
D
w
L
r
v
2
.From (32) it is clear that control v
2 affects only z
3, and control v
1 affects simultaneously z
1 and z
2. If a certain value of z
3 is reached at a selected value v
2, then v
1 may be chosen for the known z
3, providing either a given z
1 or given z
2. A pair of variables z
1
,
z
2 can simultaneously reach specified values x
1
∗ and x
2
∗ only if those values are agreed.However, it is more interesting to bring the robot to a given pointz
1
,
z
2 at a given time from a practical point of view. Therefore, it is necessary to obtain a robot control variant, which takes (30) into account.For this purpose, the control problem should be modified by assigning the other criterion function instead of (27).(33)
J
z
t
=
F
1
z
1
t
-
x
1
∗
2
+
F
3
z
3
t
-
arctg
⁡
x
2
∗
t
-
z
2
t
x
1
∗
-
z
1
t
2This function reflects the following control objective: bringing the coordinatez
1 to a value x
1
∗, and the coordinate φ
≡
z
3 to such a value that will provide a value x
2
∗ according to (30). As movement trajectory of the robot is curvilinear, at each current moment the angle of convergence to the target is variable and is reckoned from the attained state. This explains the use of the differences under the sign of arctangent in (33). Algorithm (24) modification for the case of transforms to new expressions for J
z
t (33) instead of (27) and for the vector components:(34)
∂
J
z
∂
z
=
∂
J
z
∂
z
1
∂
J
z
∂
z
2
∂
J
z
∂
x
3
.The latter are now much more complicated than in the case of criterion in the form of (27). This may be especially critical when implementing on the onboard robot controller, whose software does not include inverse trigonometric functions. Expressions for the components have the following form (we omit argument t for the sake of brevity):(35)
∂
J
z
∂
z
1
=
2
F
1
z
1
-
x
1
∗
+
2
F
3
z
2
-
x
2
∗
z
3
-
arctg
⁡
z
2
-
x
2
∗
/
z
1
-
x
1
∗
x
2
∗
z
1
-
x
1
∗
2
z
2
-
x
2
∗
2
/
z
1
-
x
1
∗
2
+
1
,
(36)
∂
J
z
∂
z
2
=
2
F
3
z
3
-
arctg
⁡
z
2
-
x
2
∗
/
z
1
-
x
1
∗
x
2
∗
z
1
-
x
1
∗
z
2
-
x
2
∗
2
/
z
1
-
x
1
∗
2
+
1
,
(37)
∂
J
z
∂
z
3
=
2
F
3
z
3
-
arctg
⁡
z
2
-
x
2
∗
z
1
-
x
1
∗
.Control by (24) with criterion by (33) was verified by Mathcad modeling, in particular, for the case of x
1
∗
=
1.5 m, x
2
∗
=
0.8 m, t
f
= 60 s, integration step h
= 0.001 s, F
1
=
10, F
3
=
1. The specified point is reached with a sufficiently high accuracy (1.488 m versus 1.5 m for x
1, and 0.784 m versus 0.8 m for x
2). In other tried cases, the accuracy is approximately the same. Figures 1–3 show time dependencies graphs of coordinates and trajectory of the robot.Figure 1
Time dependence of the robot left wheel coordinate.Figure 2
Time dependence of the robot right wheel coordinate.Figure 3
Robot trajectory.Figures4 and 5 show time dependencies graphs of the robot control actions.Figure 4
Time dependence of the robot left wheel angular velocity (control action).Figure 5
Time dependence of the robot right wheel angular velocity (control action).


### 4.2. Vessel Motion Control

#### 4.2.1. Heading Angle Control by Simultaneous Rotation of Two Azipod Units at the Same Angle

For the denoted problem, criterion by (3) is specified in the form:(38)
J
=
F
6
z
6
t
i
+
1
-
z
6
∗
2
+
F
7
z
7
t
i
+
1
-
z
7
∗
2
→
J
∗
.It should be noted that we had to add a term containingz
7 (Azipod rotation angle) to the criterion by (38), because without this term the TSM control cannot be obtained.The matrices from the system of (4) have the form of (12), (13). Applying (5), (6) to this case, we obtain the following:(39)
∂
J
∂
z
-
W
=
0
0
0
0
0
2
F
6
z
-
6
-
z
6
∗
2
F
7
z
-
7
-
z
7
∗
·
W
11
…
W
17
…
…
…
W
71
…
W
77
=
2
a
1
a
2
a
3
a
4
a
5
a
6
a
7
,
a
i
=
W
6
i
F
6
z
-
6
-
z
6
∗
+
W
7
i
F
7
z
-
7
-
z
7
∗
,
i
=
1
,
…
,
7
,where(40)
∂
J
∂
z
-
=
∂
J
∂
z
z
=
z
-
t
i
+
1
,
t
i
,
z
t
i
,
z
-
=
z
-
t
i
+
1
,
t
i
,
z
t
i
.For the case when the drive (integrating) inertia is not taken into account, we have(41)
G
=
2
a
1
a
2
a
3
a
4
a
5
a
6
a
7
·
0
0
0
0
0
0
1
=
2
a
7
=
2
W
67
F
6
z
-
6
-
z
6
∗
+
2
W
77
F
7
z
-
7
-
z
7
∗
=
2
F
7
W
77
z
-
7
-
z
7
∗
.Herewith, according to (23),(42)
W
67
=
0
+
∂
Φ
6
∂
z
7
·
h
=
∂
z
3
∂
z
7
·
h
=
0
,
W
77
=
1
+
∂
Φ
7
∂
z
7
·
h
,where Φ
7
=
0 and W
77
=
1.The expression for the pseudoinverse matrix is the following:(43)
G
+
=
G
T
G
G
T
-
1
=
1
2
F
7
·
W
77
·
z
-
7
-
z
7
∗
=
1
2
F
7
·
z
-
7
-
z
7
∗and the expression for TSV is(44)
z
-
6
t
i
+
1
,
t
i
,
z
t
i
=
z
6
t
i
+
Φ
6
t
i
,
z
·
h
=
z
6
t
i
+
h
·
z
3
t
i
≡
z
6
+
h
·
z
3
,
z
-
7
t
i
+
1
,
t
i
,
z
t
i
=
z
7
t
i
+
Φ
7
t
i
,
z
·
h
=
z
7
t
i
.As a result, we obtain the following expression for the control action at the integrating drive input:(45)
u
t
,
z
=
J
∗
-
F
6
z
6
+
h
·
z
3
-
z
6
∗
2
-
F
7
z
7
-
z
7
∗
2
2
T
u
F
7
z
7
-
z
7
∗
,where z
i
≡
z
i
t are the state variables. Weights and parameter T
u are chosen depending on the goals of the system.In the algorithm of (45), h is the integration step of the differential equations; F
6
,
F
7 are the weight coefficients, and F
7
<
<
F
6 so that the control target (bringing the vessel path z
6 to a given value z
6
∗) is not distorted. Since the value of Azipod rotation angle is not of interest, z
7
∗ can actually have any value. However, during the modeling process it was revealed that the algorithm stability depends on F
7 and z
7
∗ values. The values for which algorithm of (45) stable operation is achieved in a wide range of operating conditions are F
7
=
0
,
01 and z
7
∗
=
0
,
01
·
sign
⁡
K
∗. Because of the small value of the weight coefficient F
7, criterion by (38) in case of J
∗
=
0 displays approximately the mentioned initial control target. To prevent the division by zero in (45), a small number ε can be added to the denominator. A control amplitude limit can be set for this purpose as well.Control by (23), where u is substituted by (45), leads to a given Azipod rotation angle in case of drive dynamics being considered.


#### 4.2.2. Angular Velocity and Heading Angle Control by Simultaneous Rotation of Two Azipod Units to the Same Angle

In contrast to the problem solved in the previous paragraph, now let us set the problem of controlling the vessel course by means of an implicit criterion:(46)
J
=
z
3
t
→
z
3
∗
t
,where z
3
∗
t
=
ω
∗
t is the specified behavior of the angular velocity. From the sixth equation of the system of (4), (12) for the heading angle d
z
6
/
d
t
=
z
3, at a constant angular velocity z
3
t
=
c
o
n
s
t
=
z
3
∗
≡
ω
∗, it follows that(47)
z
6
∗
=
z
6
t
0
+
∫
t
0
t
1
ω
∗
d
t
=
z
6
t
0
+
ω
∗
·
t
∗
-
t
0
,
t
∗
=
t
0
+
z
6
∗
-
z
6
t
0
ω
∗
,
z
3
∗
t
=
if
t
<
t
∗
,
ω
∗
,
0
.The second expression of (47) determines the moment t
∗ when the vessel reaches the given route z
6
∗. Therefore, the time-varying master control z
3
∗
t is evaluated using a conditional expression written in the third line of (47) according to the rule: if the condition t
<
t
∗ is true, then z
3
∗
t
=
ω
∗; otherwise z
3
∗
t
=
0.To evaluate the control, we apply (5), (6) and get(48)
G
=
0
0
1
0
0
0
0
·
W
11
…
W
17
⋮
…
⋮
W
71
…
W
77
·
0
0
0
0
0
0
1
=
W
37
,
G
+
=
1
G
=
1
W
37
.In accordance with (8) and taking into account the definition of Φ
3 from (12):(49)
Φ
3
t
,
z
=
1
J
ω
·
L
·
T
v
·
sin
⁡
z
7
+
M
K
z
1
,
z
2
+
M
D
z
1
,
z
2
-
m
x
-
m
y
·
z
1
·
z
2
,the element W
37 of transition matrix W is evaluated by the following expression:(50)
W
37
t
i
+
1
,
t
i
,
z
t
i
=
0
+
h
·
∂
Φ
3
t
i
,
z
∂
z
7
z
7
=
z
7
t
i
=
h
·
L
·
T
v
·
cos
⁡
z
7
J
ω
.Then, forf
J
=
z
3
∗
-
z
-
3
t
i
+
1
,
t
i
,
z
t
i
/
T
u, we get the following expression for control:(51)
u
t
,
z
=
z
3
∗
-
z
3
¯
t
T
u
W
37
t
,
z
·
,
z
3
¯
t
=
z
3
+
Φ
3
t
,
z
·
h
.Simulation showed that, in any case, for the vessel in question, instead of the latter formula, the simplified one can be used. This follows from the fact thatΦ
3 (angular acceleration) has sufficiently small value and from the second relation of (51):(52)
u
t
,
z
=
z
3
∗
-
z
3
t
T
u
W
37
t
,
z
.In (51), (52) h is a time step of the control values calculation; T
u is a given time constant of the exponential motion of a variable z
3
¯
t to a given value z
3
∗ according to the equation: (53)
d
z
3
¯
d
t
=
-
z
3
∗
-
z
3
¯
T
u
.Control by (23), where u is substituted by (51) or (52), depending on used control (original or simplified), will result in a given Azipod rotation angle in case of drive dynamics being considered.Thus, two algorithms for bringing the vessel to a given course were obtained: by means of control by (45) and by means of controls by (47)-(50), (51) or controls by (47)-(50), (52). The first approach is simpler in terms of the amount of calculations. Still it has low accuracy because its criterion in the form of (38) reflects an artificial target—bringing to the given values both the course and Azipod rotation angle. The second approach is more difficult to implement, but its target in the form of criterion by (46) and (47) is exactly the attainment of the given course K
∗
≡
z
6
∗.To achieve both simplicity and accuracy of computations, a composite algorithm is offered. This algorithm includes a “strong” maneuver through the control by (45). A “slight” maneuver for static error compensation is included through the control by (50), (52). In the latter case transition function W
37 is calculated by the simplified formula for small Azipod rotation angles, when cos
⁡
z
7
≈
1.(54)
W
37
t
,
z
=
h
·
L
·
T
v
J
ωApplying composite algorithm makes sense in the case of using onboard controllers, intended to perform simple computational operations. This is of course not typical for the catamaran vessel in question. It is more relevant to ordinary small vessels intended for everyday purposes, for example, for cleaning leaves and relatively large debris from the water surface.Simulation results in the form of the timet dependence of the vessel model coordinates are obtained by Mathcad modeling. The results are obtained for the model defined by (4), (12)–(19) with correction by (21) using the control by (47), (50), (52) for the given values ω
∗ (constant angular velocity of the vessel rotation towards the given course) and z
6
∗ (the given heading angle). Calculations were made for the following vessel and control algorithm parameters: D
=
155000
k
g; k
11
=
0.05; k
66
=
0.6; k
22
=
0.7; ρ
=
1025
k
g
·
m
-
3; L
=
22.3
m; J
Z
=
0.05
·
ρ
·
D
·
L
2; T
v
1
=
T
v
2
=
20400
N; T
u
=
0.1
s; k
δ
=
1; T
δ
=
2
s; h
=
0.1
s, ω
m
a
x
=
20.It is possible to regulate both trajectory and time of bringing the vessel to the specified course by choosing parametersz
3
∗
,
z
6
∗
,
T
u.Figures6 and 7 show the vessel trajectory and heading angle during the maneuver for the case of z
3
∗
=
0.005
r
a
d
·
s
-
1, z
6
∗
=
3.14
r
a
d, corresponding to a relatively slow vessel rotation towards the return course.Figure 6
Vessel trajectory during rotation towards the return course forz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d.Figure 7
Time dependence of the vessel heading angle during rotation towards the return course forz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d.Figures8 and 9 show the vessel trajectory and heading angle during the maneuver for the case of z
3
∗
=
0.005
r
a
d
·
s
-
1 and z
3
∗
=
0.05
r
a
d
·
s
-
1, corresponding to a slowed and accelerated vessel rotation at a given heading angle of 57.3 degrees (z
6
∗
=
1
r
a
d).Figure 8
Vessel trajectory during slowed rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.Figure 9
Time dependence of the vessel heading angle during slowed rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.A certain influence on the trajectory is also provided by the value of the free parameterT
u, on which transient period for angular velocity depends. Figure 9 shows how by increasing parameter T
u by a factor of twenty, a smoother vessel turn trajectory can be attained (compare with the trajectory in Figure 6 for T
u
=
0.1
s).Graphs shown in Figures6–13 illustrate how the free parameters z
3
∗ (vessel angular velocity), z
6
∗ (specified heading angle), and T
u (speed parameter) affect the shape of the vessel trajectory.Figure 10
Vessel trajectory during accelerated rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.05
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.Figure 11
Time dependence of the vessel heading angle during accelerated rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.05
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.Figure 12
Vessel trajectory during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d; T
u
=
2
s.Figure 13
Time dependence of the heading angle during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d; T
u
=
2
s.We note that in all the cases presented, the control error is close to zero (in the simulation, the first significant digit is in fourth place after the decimal point). By control error, we mean the difference between specified and attained values of the course.Graphs in Figures14 and 15 show the results of vessel motion simulation under the influence of composite control by (45), (52), (54).Figure 14
Vessel trajectory during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.4
r
a
d; T
u
=
0.1
s, F
6
=
10, F
7
=
0.01, z
7
∗
=
0.01
sign
⁡
z
6
∗.Figure 15
Time dependence of the heading angle during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.4
r
a
d; T
u
=
0.1
s, F
6
=
10, F
7
=
0.01, z
7
∗
=
0.01
⋅
sign
⁡
z
6
∗.Simulation was performed for two consecutive time intervals with a duration of500
s. During the first time interval, the vessel was controlled by algorithm of (45) with all zero initial conditions except speed V
x
0
=
4
m
·
s
-
1. During the second time interval, the vessel was controlled by algorithm of (47), (52), (54). The initial conditions for the second interval were the values of all coordinates at the end of the first interval, except for z
7
0
=
0. From a technical point of view, this means that Azipod should be brought to the zero (neutral) position at the beginning of second interval. This peculiarity of composite control is connected with the fact that the algorithm of (47), (52), (54) does not contain a feedback from z
7. Therefore, the nonzero initial condition of this coordinate is ignored by the control algorithm. If the initial condition z
7
0 for the second interval is nonzero, an uncompensated static error is formed at the end of the process.As can be seen from the graphs of Figures14 and 15, vessel reached the heading angle of 3.051 radians under the influence of control by (45). Then the vessel made a turn of 3.14 radians to the set course under the influence of control by (52), (54). Herewith, the actual time of the turn is very short—it is about one minute (more precisely, 1193 steps with a duration of 0.05 seconds, that is, 59.65 seconds).It is important to note that the control by (47), (52), (54) explicitly depends on the mass-dimensional and propulsion parameters of the vessel, as well as on parameters h
,
T
u
,
z
3
∗
≡
ω
∗ affecting the shape of the vessel track to the given course. This fact makes it possible to eliminate the problem of smoothing the trajectory (calculating anticipatory points, etc., proposed in [12]).It should be noted that the results concerning vessel control adduced in this paper represent development of methods and algorithms published in [25], regarding a definite implementation of the composite algorithm and refinement of the methods of analytical synthesis of controls as well. A Mathcad program code simulating the vessel motion with a model in the form of (12)-(19) and control by (52) is given in the mentioned paper. This program is easy to modify with (23) and composite control by (45), (52), (54) for the case of the Azipod inertial drive.


## 4.1. Robot Control

A. Isidori’s linearizing by input/output method [1] in Wei Wu’s modification [2] and TSM [4] were used in [9] to bring robot to a given point at a given time. Because of the specificity of the model described by (4), both methods led to the same algorithm:(24)
u
t
=
G
+
J
∗
-
J
z
t
T
,
G
+
=
G
T
G
G
T
-
1
,
G
=
∂
J
z
t
∂
x
B
t
,
z
t
,where J
t
≡
J
z
t is the criterion function that defines desired exponential trajectory with a time constant T
u. Differential equation for this function is the following:(25)
d
J
t
d
t
=
J
∗
-
J
t
T
u
.(24), (25) result from (5), (6) and the fact that TM W in this case is identity matrix on the entire period of the system’s operation, since the uncontrolled system vector Φ is equal to zero according to (9).In (24) state variables z
=
z
i
,
i
=
1
,
…
,
n and control variables u
=
u
j
,
j
=
1
,
…
,
r, are redefined according to(26)
z
=
x
1
x
2
φ
,
u
=
ω
1
ω
2
.For the criterion function, the following quadratic form is used:(27)
J
z
t
f
=
F
1
z
1
t
f
-
x
1
∗
2
+
F
2
z
2
t
f
-
x
2
∗
2
+
F
3
z
3
t
f
-
ω
∗
2
→
J
∗
=
0
,where z
∗
=
c
o
l
o
n
x
1
∗
,
x
2
∗
,
ω
∗ are the desired values of the robot coordinates at the finite time t
f.The sense of criterion by (27) is that when the current state z
t converges to a given point z
∗, the value of the criterion function tends to the given value J
∗, which must obviously be equal to zero. This target may be unattainable by the control (24) that results in deviation of the actual behavior by (27) from the desired trajectory (25). The main reason for unattainability is that the scalar G
G
T nulls at some points or intervals of the time axis. Obviously, this fact causes unattainability, as control target does not meet the way it is being achieved along the trajectory according to (25). Another possible reason for the deviation of the actual trajectory from the given one is the amplitude constraints on the control actions that cannot be taken into account in the methods used. Thus, rotational speeds of the robot wheels are limited in amplitude by the ω
m value:(28)
ω
i
=
-
ω
m
,
u
i
<
-
ω
m
,
u
i
,
-
ω
m
≤
u
i
≤
ω
m
,
ω
m
,
u
i
>
-
ω
m
,
i
=
1,2
.The output of the controls evaluated by (24) with constraints by (28) also distorts the form of the criterion function defined by (25). But still the criterion function decreases, converging to some minimally attainable value J
∗. Herewith attained value J
∗ determines the accuracy of bringing state vector components to given values.System of (9) is a nonlinear affine system of (4) with two control actions ω
1
,
ω
2, where Φ is zero vector, and the components of the matrix B
=
B
i
j
,
i
=
1,2
,
3
,
j
=
1,2, are the following:(29)
B
11
=
B
12
=
D
w
cos
⁡
ϕ
4
,
B
21
=
B
22
=
D
w
sin
⁡
ϕ
4
,
B
31
=
-
B
32
=
D
w
2
L
r
.In [9], the results of Mathcad modeling of the system of (9) with control by (24), criterion by (27), and amplitude constraints by (28) are applied to the ArcBotics Sparki robot with the following parameters: D
w = 0,05 m; L
r = 0.085 m; F
1 = 10; F
2 = 0; F
3 = 1; x
1
∗ = 0.6 m; x
2
∗ = 0.3 m; x
3
∗ = 0.8 m; ω
m = 1.5s
-
1; h = 0.1 s; Tu = 0.01 s; t
f = 40 s.The same algorithm is implemented in the onboard C++ version of a real robot. As follows from the kinematic scheme, which is the basis for (9), there is a nonlinear dependence between the states of the system of (9).(30)
tg
⁡
z
3
=
z
2
z
1Consequently, equations (9) are dependent. It means that during control by (24) implementation, specified values can be reached simultaneously for a pair of variables x
1
,
φ
z
1
,
z
3 or for a pair x
2
,
φ
z
2
,
z
3. As for the pair x
1
,
x
2
z
1
,
z
2, this target is generally unattainable.This is clear from (9) with other control actions:(31)
v
1
=
ω
1
+
ω
2
2
,
v
2
=
ω
1
-
ω
2
2
.The resulting system is(32)
d
z
1
d
t
=
D
w
2
v
1
cos
⁡
ϕ
,
d
z
2
d
t
=
D
w
2
v
1
sin
⁡
ϕ
,
d
z
3
d
t
=
D
w
L
r
v
2
.From (32) it is clear that control v
2 affects only z
3, and control v
1 affects simultaneously z
1 and z
2. If a certain value of z
3 is reached at a selected value v
2, then v
1 may be chosen for the known z
3, providing either a given z
1 or given z
2. A pair of variables z
1
,
z
2 can simultaneously reach specified values x
1
∗ and x
2
∗ only if those values are agreed.However, it is more interesting to bring the robot to a given pointz
1
,
z
2 at a given time from a practical point of view. Therefore, it is necessary to obtain a robot control variant, which takes (30) into account.For this purpose, the control problem should be modified by assigning the other criterion function instead of (27).(33)
J
z
t
=
F
1
z
1
t
-
x
1
∗
2
+
F
3
z
3
t
-
arctg
⁡
x
2
∗
t
-
z
2
t
x
1
∗
-
z
1
t
2This function reflects the following control objective: bringing the coordinatez
1 to a value x
1
∗, and the coordinate φ
≡
z
3 to such a value that will provide a value x
2
∗ according to (30). As movement trajectory of the robot is curvilinear, at each current moment the angle of convergence to the target is variable and is reckoned from the attained state. This explains the use of the differences under the sign of arctangent in (33). Algorithm (24) modification for the case of transforms to new expressions for J
z
t (33) instead of (27) and for the vector components:(34)
∂
J
z
∂
z
=
∂
J
z
∂
z
1
∂
J
z
∂
z
2
∂
J
z
∂
x
3
.The latter are now much more complicated than in the case of criterion in the form of (27). This may be especially critical when implementing on the onboard robot controller, whose software does not include inverse trigonometric functions. Expressions for the components have the following form (we omit argument t for the sake of brevity):(35)
∂
J
z
∂
z
1
=
2
F
1
z
1
-
x
1
∗
+
2
F
3
z
2
-
x
2
∗
z
3
-
arctg
⁡
z
2
-
x
2
∗
/
z
1
-
x
1
∗
x
2
∗
z
1
-
x
1
∗
2
z
2
-
x
2
∗
2
/
z
1
-
x
1
∗
2
+
1
,
(36)
∂
J
z
∂
z
2
=
2
F
3
z
3
-
arctg
⁡
z
2
-
x
2
∗
/
z
1
-
x
1
∗
x
2
∗
z
1
-
x
1
∗
z
2
-
x
2
∗
2
/
z
1
-
x
1
∗
2
+
1
,
(37)
∂
J
z
∂
z
3
=
2
F
3
z
3
-
arctg
⁡
z
2
-
x
2
∗
z
1
-
x
1
∗
.Control by (24) with criterion by (33) was verified by Mathcad modeling, in particular, for the case of x
1
∗
=
1.5 m, x
2
∗
=
0.8 m, t
f
= 60 s, integration step h
= 0.001 s, F
1
=
10, F
3
=
1. The specified point is reached with a sufficiently high accuracy (1.488 m versus 1.5 m for x
1, and 0.784 m versus 0.8 m for x
2). In other tried cases, the accuracy is approximately the same. Figures 1–3 show time dependencies graphs of coordinates and trajectory of the robot.Figure 1
Time dependence of the robot left wheel coordinate.Figure 2
Time dependence of the robot right wheel coordinate.Figure 3
Robot trajectory.Figures4 and 5 show time dependencies graphs of the robot control actions.Figure 4
Time dependence of the robot left wheel angular velocity (control action).Figure 5
Time dependence of the robot right wheel angular velocity (control action).


## 4.2. Vessel Motion Control

### 4.2.1. Heading Angle Control by Simultaneous Rotation of Two Azipod Units at the Same Angle

For the denoted problem, criterion by (3) is specified in the form:(38)
J
=
F
6
z
6
t
i
+
1
-
z
6
∗
2
+
F
7
z
7
t
i
+
1
-
z
7
∗
2
→
J
∗
.It should be noted that we had to add a term containingz
7 (Azipod rotation angle) to the criterion by (38), because without this term the TSM control cannot be obtained.The matrices from the system of (4) have the form of (12), (13). Applying (5), (6) to this case, we obtain the following:(39)
∂
J
∂
z
-
W
=
0
0
0
0
0
2
F
6
z
-
6
-
z
6
∗
2
F
7
z
-
7
-
z
7
∗
·
W
11
…
W
17
…
…
…
W
71
…
W
77
=
2
a
1
a
2
a
3
a
4
a
5
a
6
a
7
,
a
i
=
W
6
i
F
6
z
-
6
-
z
6
∗
+
W
7
i
F
7
z
-
7
-
z
7
∗
,
i
=
1
,
…
,
7
,where(40)
∂
J
∂
z
-
=
∂
J
∂
z
z
=
z
-
t
i
+
1
,
t
i
,
z
t
i
,
z
-
=
z
-
t
i
+
1
,
t
i
,
z
t
i
.For the case when the drive (integrating) inertia is not taken into account, we have(41)
G
=
2
a
1
a
2
a
3
a
4
a
5
a
6
a
7
·
0
0
0
0
0
0
1
=
2
a
7
=
2
W
67
F
6
z
-
6
-
z
6
∗
+
2
W
77
F
7
z
-
7
-
z
7
∗
=
2
F
7
W
77
z
-
7
-
z
7
∗
.Herewith, according to (23),(42)
W
67
=
0
+
∂
Φ
6
∂
z
7
·
h
=
∂
z
3
∂
z
7
·
h
=
0
,
W
77
=
1
+
∂
Φ
7
∂
z
7
·
h
,where Φ
7
=
0 and W
77
=
1.The expression for the pseudoinverse matrix is the following:(43)
G
+
=
G
T
G
G
T
-
1
=
1
2
F
7
·
W
77
·
z
-
7
-
z
7
∗
=
1
2
F
7
·
z
-
7
-
z
7
∗and the expression for TSV is(44)
z
-
6
t
i
+
1
,
t
i
,
z
t
i
=
z
6
t
i
+
Φ
6
t
i
,
z
·
h
=
z
6
t
i
+
h
·
z
3
t
i
≡
z
6
+
h
·
z
3
,
z
-
7
t
i
+
1
,
t
i
,
z
t
i
=
z
7
t
i
+
Φ
7
t
i
,
z
·
h
=
z
7
t
i
.As a result, we obtain the following expression for the control action at the integrating drive input:(45)
u
t
,
z
=
J
∗
-
F
6
z
6
+
h
·
z
3
-
z
6
∗
2
-
F
7
z
7
-
z
7
∗
2
2
T
u
F
7
z
7
-
z
7
∗
,where z
i
≡
z
i
t are the state variables. Weights and parameter T
u are chosen depending on the goals of the system.In the algorithm of (45), h is the integration step of the differential equations; F
6
,
F
7 are the weight coefficients, and F
7
<
<
F
6 so that the control target (bringing the vessel path z
6 to a given value z
6
∗) is not distorted. Since the value of Azipod rotation angle is not of interest, z
7
∗ can actually have any value. However, during the modeling process it was revealed that the algorithm stability depends on F
7 and z
7
∗ values. The values for which algorithm of (45) stable operation is achieved in a wide range of operating conditions are F
7
=
0
,
01 and z
7
∗
=
0
,
01
·
sign
⁡
K
∗. Because of the small value of the weight coefficient F
7, criterion by (38) in case of J
∗
=
0 displays approximately the mentioned initial control target. To prevent the division by zero in (45), a small number ε can be added to the denominator. A control amplitude limit can be set for this purpose as well.Control by (23), where u is substituted by (45), leads to a given Azipod rotation angle in case of drive dynamics being considered.


### 4.2.2. Angular Velocity and Heading Angle Control by Simultaneous Rotation of Two Azipod Units to the Same Angle

In contrast to the problem solved in the previous paragraph, now let us set the problem of controlling the vessel course by means of an implicit criterion:(46)
J
=
z
3
t
→
z
3
∗
t
,where z
3
∗
t
=
ω
∗
t is the specified behavior of the angular velocity. From the sixth equation of the system of (4), (12) for the heading angle d
z
6
/
d
t
=
z
3, at a constant angular velocity z
3
t
=
c
o
n
s
t
=
z
3
∗
≡
ω
∗, it follows that(47)
z
6
∗
=
z
6
t
0
+
∫
t
0
t
1
ω
∗
d
t
=
z
6
t
0
+
ω
∗
·
t
∗
-
t
0
,
t
∗
=
t
0
+
z
6
∗
-
z
6
t
0
ω
∗
,
z
3
∗
t
=
if
t
<
t
∗
,
ω
∗
,
0
.The second expression of (47) determines the moment t
∗ when the vessel reaches the given route z
6
∗. Therefore, the time-varying master control z
3
∗
t is evaluated using a conditional expression written in the third line of (47) according to the rule: if the condition t
<
t
∗ is true, then z
3
∗
t
=
ω
∗; otherwise z
3
∗
t
=
0.To evaluate the control, we apply (5), (6) and get(48)
G
=
0
0
1
0
0
0
0
·
W
11
…
W
17
⋮
…
⋮
W
71
…
W
77
·
0
0
0
0
0
0
1
=
W
37
,
G
+
=
1
G
=
1
W
37
.In accordance with (8) and taking into account the definition of Φ
3 from (12):(49)
Φ
3
t
,
z
=
1
J
ω
·
L
·
T
v
·
sin
⁡
z
7
+
M
K
z
1
,
z
2
+
M
D
z
1
,
z
2
-
m
x
-
m
y
·
z
1
·
z
2
,the element W
37 of transition matrix W is evaluated by the following expression:(50)
W
37
t
i
+
1
,
t
i
,
z
t
i
=
0
+
h
·
∂
Φ
3
t
i
,
z
∂
z
7
z
7
=
z
7
t
i
=
h
·
L
·
T
v
·
cos
⁡
z
7
J
ω
.Then, forf
J
=
z
3
∗
-
z
-
3
t
i
+
1
,
t
i
,
z
t
i
/
T
u, we get the following expression for control:(51)
u
t
,
z
=
z
3
∗
-
z
3
¯
t
T
u
W
37
t
,
z
·
,
z
3
¯
t
=
z
3
+
Φ
3
t
,
z
·
h
.Simulation showed that, in any case, for the vessel in question, instead of the latter formula, the simplified one can be used. This follows from the fact thatΦ
3 (angular acceleration) has sufficiently small value and from the second relation of (51):(52)
u
t
,
z
=
z
3
∗
-
z
3
t
T
u
W
37
t
,
z
.In (51), (52) h is a time step of the control values calculation; T
u is a given time constant of the exponential motion of a variable z
3
¯
t to a given value z
3
∗ according to the equation: (53)
d
z
3
¯
d
t
=
-
z
3
∗
-
z
3
¯
T
u
.Control by (23), where u is substituted by (51) or (52), depending on used control (original or simplified), will result in a given Azipod rotation angle in case of drive dynamics being considered.Thus, two algorithms for bringing the vessel to a given course were obtained: by means of control by (45) and by means of controls by (47)-(50), (51) or controls by (47)-(50), (52). The first approach is simpler in terms of the amount of calculations. Still it has low accuracy because its criterion in the form of (38) reflects an artificial target—bringing to the given values both the course and Azipod rotation angle. The second approach is more difficult to implement, but its target in the form of criterion by (46) and (47) is exactly the attainment of the given course K
∗
≡
z
6
∗.To achieve both simplicity and accuracy of computations, a composite algorithm is offered. This algorithm includes a “strong” maneuver through the control by (45). A “slight” maneuver for static error compensation is included through the control by (50), (52). In the latter case transition function W
37 is calculated by the simplified formula for small Azipod rotation angles, when cos
⁡
z
7
≈
1.(54)
W
37
t
,
z
=
h
·
L
·
T
v
J
ωApplying composite algorithm makes sense in the case of using onboard controllers, intended to perform simple computational operations. This is of course not typical for the catamaran vessel in question. It is more relevant to ordinary small vessels intended for everyday purposes, for example, for cleaning leaves and relatively large debris from the water surface.Simulation results in the form of the timet dependence of the vessel model coordinates are obtained by Mathcad modeling. The results are obtained for the model defined by (4), (12)–(19) with correction by (21) using the control by (47), (50), (52) for the given values ω
∗ (constant angular velocity of the vessel rotation towards the given course) and z
6
∗ (the given heading angle). Calculations were made for the following vessel and control algorithm parameters: D
=
155000
k
g; k
11
=
0.05; k
66
=
0.6; k
22
=
0.7; ρ
=
1025
k
g
·
m
-
3; L
=
22.3
m; J
Z
=
0.05
·
ρ
·
D
·
L
2; T
v
1
=
T
v
2
=
20400
N; T
u
=
0.1
s; k
δ
=
1; T
δ
=
2
s; h
=
0.1
s, ω
m
a
x
=
20.It is possible to regulate both trajectory and time of bringing the vessel to the specified course by choosing parametersz
3
∗
,
z
6
∗
,
T
u.Figures6 and 7 show the vessel trajectory and heading angle during the maneuver for the case of z
3
∗
=
0.005
r
a
d
·
s
-
1, z
6
∗
=
3.14
r
a
d, corresponding to a relatively slow vessel rotation towards the return course.Figure 6
Vessel trajectory during rotation towards the return course forz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d.Figure 7
Time dependence of the vessel heading angle during rotation towards the return course forz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d.Figures8 and 9 show the vessel trajectory and heading angle during the maneuver for the case of z
3
∗
=
0.005
r
a
d
·
s
-
1 and z
3
∗
=
0.05
r
a
d
·
s
-
1, corresponding to a slowed and accelerated vessel rotation at a given heading angle of 57.3 degrees (z
6
∗
=
1
r
a
d).Figure 8
Vessel trajectory during slowed rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.Figure 9
Time dependence of the vessel heading angle during slowed rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.A certain influence on the trajectory is also provided by the value of the free parameterT
u, on which transient period for angular velocity depends. Figure 9 shows how by increasing parameter T
u by a factor of twenty, a smoother vessel turn trajectory can be attained (compare with the trajectory in Figure 6 for T
u
=
0.1
s).Graphs shown in Figures6–13 illustrate how the free parameters z
3
∗ (vessel angular velocity), z
6
∗ (specified heading angle), and T
u (speed parameter) affect the shape of the vessel trajectory.Figure 10
Vessel trajectory during accelerated rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.05
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.Figure 11
Time dependence of the vessel heading angle during accelerated rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.05
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.Figure 12
Vessel trajectory during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d; T
u
=
2
s.Figure 13
Time dependence of the heading angle during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d; T
u
=
2
s.We note that in all the cases presented, the control error is close to zero (in the simulation, the first significant digit is in fourth place after the decimal point). By control error, we mean the difference between specified and attained values of the course.Graphs in Figures14 and 15 show the results of vessel motion simulation under the influence of composite control by (45), (52), (54).Figure 14
Vessel trajectory during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.4
r
a
d; T
u
=
0.1
s, F
6
=
10, F
7
=
0.01, z
7
∗
=
0.01
sign
⁡
z
6
∗.Figure 15
Time dependence of the heading angle during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.4
r
a
d; T
u
=
0.1
s, F
6
=
10, F
7
=
0.01, z
7
∗
=
0.01
⋅
sign
⁡
z
6
∗.Simulation was performed for two consecutive time intervals with a duration of500
s. During the first time interval, the vessel was controlled by algorithm of (45) with all zero initial conditions except speed V
x
0
=
4
m
·
s
-
1. During the second time interval, the vessel was controlled by algorithm of (47), (52), (54). The initial conditions for the second interval were the values of all coordinates at the end of the first interval, except for z
7
0
=
0. From a technical point of view, this means that Azipod should be brought to the zero (neutral) position at the beginning of second interval. This peculiarity of composite control is connected with the fact that the algorithm of (47), (52), (54) does not contain a feedback from z
7. Therefore, the nonzero initial condition of this coordinate is ignored by the control algorithm. If the initial condition z
7
0 for the second interval is nonzero, an uncompensated static error is formed at the end of the process.As can be seen from the graphs of Figures14 and 15, vessel reached the heading angle of 3.051 radians under the influence of control by (45). Then the vessel made a turn of 3.14 radians to the set course under the influence of control by (52), (54). Herewith, the actual time of the turn is very short—it is about one minute (more precisely, 1193 steps with a duration of 0.05 seconds, that is, 59.65 seconds).It is important to note that the control by (47), (52), (54) explicitly depends on the mass-dimensional and propulsion parameters of the vessel, as well as on parameters h
,
T
u
,
z
3
∗
≡
ω
∗ affecting the shape of the vessel track to the given course. This fact makes it possible to eliminate the problem of smoothing the trajectory (calculating anticipatory points, etc., proposed in [12]).It should be noted that the results concerning vessel control adduced in this paper represent development of methods and algorithms published in [25], regarding a definite implementation of the composite algorithm and refinement of the methods of analytical synthesis of controls as well. A Mathcad program code simulating the vessel motion with a model in the form of (12)-(19) and control by (52) is given in the mentioned paper. This program is easy to modify with (23) and composite control by (45), (52), (54) for the case of the Azipod inertial drive.


## 4.2.1. Heading Angle Control by Simultaneous Rotation of Two Azipod Units at the Same Angle

For the denoted problem, criterion by (3) is specified in the form:(38)
J
=
F
6
z
6
t
i
+
1
-
z
6
∗
2
+
F
7
z
7
t
i
+
1
-
z
7
∗
2
→
J
∗
.It should be noted that we had to add a term containingz
7 (Azipod rotation angle) to the criterion by (38), because without this term the TSM control cannot be obtained.The matrices from the system of (4) have the form of (12), (13). Applying (5), (6) to this case, we obtain the following:(39)
∂
J
∂
z
-
W
=
0
0
0
0
0
2
F
6
z
-
6
-
z
6
∗
2
F
7
z
-
7
-
z
7
∗
·
W
11
…
W
17
…
…
…
W
71
…
W
77
=
2
a
1
a
2
a
3
a
4
a
5
a
6
a
7
,
a
i
=
W
6
i
F
6
z
-
6
-
z
6
∗
+
W
7
i
F
7
z
-
7
-
z
7
∗
,
i
=
1
,
…
,
7
,where(40)
∂
J
∂
z
-
=
∂
J
∂
z
z
=
z
-
t
i
+
1
,
t
i
,
z
t
i
,
z
-
=
z
-
t
i
+
1
,
t
i
,
z
t
i
.For the case when the drive (integrating) inertia is not taken into account, we have(41)
G
=
2
a
1
a
2
a
3
a
4
a
5
a
6
a
7
·
0
0
0
0
0
0
1
=
2
a
7
=
2
W
67
F
6
z
-
6
-
z
6
∗
+
2
W
77
F
7
z
-
7
-
z
7
∗
=
2
F
7
W
77
z
-
7
-
z
7
∗
.Herewith, according to (23),(42)
W
67
=
0
+
∂
Φ
6
∂
z
7
·
h
=
∂
z
3
∂
z
7
·
h
=
0
,
W
77
=
1
+
∂
Φ
7
∂
z
7
·
h
,where Φ
7
=
0 and W
77
=
1.The expression for the pseudoinverse matrix is the following:(43)
G
+
=
G
T
G
G
T
-
1
=
1
2
F
7
·
W
77
·
z
-
7
-
z
7
∗
=
1
2
F
7
·
z
-
7
-
z
7
∗and the expression for TSV is(44)
z
-
6
t
i
+
1
,
t
i
,
z
t
i
=
z
6
t
i
+
Φ
6
t
i
,
z
·
h
=
z
6
t
i
+
h
·
z
3
t
i
≡
z
6
+
h
·
z
3
,
z
-
7
t
i
+
1
,
t
i
,
z
t
i
=
z
7
t
i
+
Φ
7
t
i
,
z
·
h
=
z
7
t
i
.As a result, we obtain the following expression for the control action at the integrating drive input:(45)
u
t
,
z
=
J
∗
-
F
6
z
6
+
h
·
z
3
-
z
6
∗
2
-
F
7
z
7
-
z
7
∗
2
2
T
u
F
7
z
7
-
z
7
∗
,where z
i
≡
z
i
t are the state variables. Weights and parameter T
u are chosen depending on the goals of the system.In the algorithm of (45), h is the integration step of the differential equations; F
6
,
F
7 are the weight coefficients, and F
7
<
<
F
6 so that the control target (bringing the vessel path z
6 to a given value z
6
∗) is not distorted. Since the value of Azipod rotation angle is not of interest, z
7
∗ can actually have any value. However, during the modeling process it was revealed that the algorithm stability depends on F
7 and z
7
∗ values. The values for which algorithm of (45) stable operation is achieved in a wide range of operating conditions are F
7
=
0
,
01 and z
7
∗
=
0
,
01
·
sign
⁡
K
∗. Because of the small value of the weight coefficient F
7, criterion by (38) in case of J
∗
=
0 displays approximately the mentioned initial control target. To prevent the division by zero in (45), a small number ε can be added to the denominator. A control amplitude limit can be set for this purpose as well.Control by (23), where u is substituted by (45), leads to a given Azipod rotation angle in case of drive dynamics being considered.


## 4.2.2. Angular Velocity and Heading Angle Control by Simultaneous Rotation of Two Azipod Units to the Same Angle

In contrast to the problem solved in the previous paragraph, now let us set the problem of controlling the vessel course by means of an implicit criterion:(46)
J
=
z
3
t
→
z
3
∗
t
,where z
3
∗
t
=
ω
∗
t is the specified behavior of the angular velocity. From the sixth equation of the system of (4), (12) for the heading angle d
z
6
/
d
t
=
z
3, at a constant angular velocity z
3
t
=
c
o
n
s
t
=
z
3
∗
≡
ω
∗, it follows that(47)
z
6
∗
=
z
6
t
0
+
∫
t
0
t
1
ω
∗
d
t
=
z
6
t
0
+
ω
∗
·
t
∗
-
t
0
,
t
∗
=
t
0
+
z
6
∗
-
z
6
t
0
ω
∗
,
z
3
∗
t
=
if
t
<
t
∗
,
ω
∗
,
0
.The second expression of (47) determines the moment t
∗ when the vessel reaches the given route z
6
∗. Therefore, the time-varying master control z
3
∗
t is evaluated using a conditional expression written in the third line of (47) according to the rule: if the condition t
<
t
∗ is true, then z
3
∗
t
=
ω
∗; otherwise z
3
∗
t
=
0.To evaluate the control, we apply (5), (6) and get(48)
G
=
0
0
1
0
0
0
0
·
W
11
…
W
17
⋮
…
⋮
W
71
…
W
77
·
0
0
0
0
0
0
1
=
W
37
,
G
+
=
1
G
=
1
W
37
.In accordance with (8) and taking into account the definition of Φ
3 from (12):(49)
Φ
3
t
,
z
=
1
J
ω
·
L
·
T
v
·
sin
⁡
z
7
+
M
K
z
1
,
z
2
+
M
D
z
1
,
z
2
-
m
x
-
m
y
·
z
1
·
z
2
,the element W
37 of transition matrix W is evaluated by the following expression:(50)
W
37
t
i
+
1
,
t
i
,
z
t
i
=
0
+
h
·
∂
Φ
3
t
i
,
z
∂
z
7
z
7
=
z
7
t
i
=
h
·
L
·
T
v
·
cos
⁡
z
7
J
ω
.Then, forf
J
=
z
3
∗
-
z
-
3
t
i
+
1
,
t
i
,
z
t
i
/
T
u, we get the following expression for control:(51)
u
t
,
z
=
z
3
∗
-
z
3
¯
t
T
u
W
37
t
,
z
·
,
z
3
¯
t
=
z
3
+
Φ
3
t
,
z
·
h
.Simulation showed that, in any case, for the vessel in question, instead of the latter formula, the simplified one can be used. This follows from the fact thatΦ
3 (angular acceleration) has sufficiently small value and from the second relation of (51):(52)
u
t
,
z
=
z
3
∗
-
z
3
t
T
u
W
37
t
,
z
.In (51), (52) h is a time step of the control values calculation; T
u is a given time constant of the exponential motion of a variable z
3
¯
t to a given value z
3
∗ according to the equation: (53)
d
z
3
¯
d
t
=
-
z
3
∗
-
z
3
¯
T
u
.Control by (23), where u is substituted by (51) or (52), depending on used control (original or simplified), will result in a given Azipod rotation angle in case of drive dynamics being considered.Thus, two algorithms for bringing the vessel to a given course were obtained: by means of control by (45) and by means of controls by (47)-(50), (51) or controls by (47)-(50), (52). The first approach is simpler in terms of the amount of calculations. Still it has low accuracy because its criterion in the form of (38) reflects an artificial target—bringing to the given values both the course and Azipod rotation angle. The second approach is more difficult to implement, but its target in the form of criterion by (46) and (47) is exactly the attainment of the given course K
∗
≡
z
6
∗.To achieve both simplicity and accuracy of computations, a composite algorithm is offered. This algorithm includes a “strong” maneuver through the control by (45). A “slight” maneuver for static error compensation is included through the control by (50), (52). In the latter case transition function W
37 is calculated by the simplified formula for small Azipod rotation angles, when cos
⁡
z
7
≈
1.(54)
W
37
t
,
z
=
h
·
L
·
T
v
J
ωApplying composite algorithm makes sense in the case of using onboard controllers, intended to perform simple computational operations. This is of course not typical for the catamaran vessel in question. It is more relevant to ordinary small vessels intended for everyday purposes, for example, for cleaning leaves and relatively large debris from the water surface.Simulation results in the form of the timet dependence of the vessel model coordinates are obtained by Mathcad modeling. The results are obtained for the model defined by (4), (12)–(19) with correction by (21) using the control by (47), (50), (52) for the given values ω
∗ (constant angular velocity of the vessel rotation towards the given course) and z
6
∗ (the given heading angle). Calculations were made for the following vessel and control algorithm parameters: D
=
155000
k
g; k
11
=
0.05; k
66
=
0.6; k
22
=
0.7; ρ
=
1025
k
g
·
m
-
3; L
=
22.3
m; J
Z
=
0.05
·
ρ
·
D
·
L
2; T
v
1
=
T
v
2
=
20400
N; T
u
=
0.1
s; k
δ
=
1; T
δ
=
2
s; h
=
0.1
s, ω
m
a
x
=
20.It is possible to regulate both trajectory and time of bringing the vessel to the specified course by choosing parametersz
3
∗
,
z
6
∗
,
T
u.Figures6 and 7 show the vessel trajectory and heading angle during the maneuver for the case of z
3
∗
=
0.005
r
a
d
·
s
-
1, z
6
∗
=
3.14
r
a
d, corresponding to a relatively slow vessel rotation towards the return course.Figure 6
Vessel trajectory during rotation towards the return course forz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d.Figure 7
Time dependence of the vessel heading angle during rotation towards the return course forz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d.Figures8 and 9 show the vessel trajectory and heading angle during the maneuver for the case of z
3
∗
=
0.005
r
a
d
·
s
-
1 and z
3
∗
=
0.05
r
a
d
·
s
-
1, corresponding to a slowed and accelerated vessel rotation at a given heading angle of 57.3 degrees (z
6
∗
=
1
r
a
d).Figure 8
Vessel trajectory during slowed rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.Figure 9
Time dependence of the vessel heading angle during slowed rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.A certain influence on the trajectory is also provided by the value of the free parameterT
u, on which transient period for angular velocity depends. Figure 9 shows how by increasing parameter T
u by a factor of twenty, a smoother vessel turn trajectory can be attained (compare with the trajectory in Figure 6 for T
u
=
0.1
s).Graphs shown in Figures6–13 illustrate how the free parameters z
3
∗ (vessel angular velocity), z
6
∗ (specified heading angle), and T
u (speed parameter) affect the shape of the vessel trajectory.Figure 10
Vessel trajectory during accelerated rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.05
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.Figure 11
Time dependence of the vessel heading angle during accelerated rotation at a given heading angle of 57.3 degrees for valuesz
3
∗
=
0.05
r
a
d
⋅
s
-
1, z
6
∗
=
1
r
a
d.Figure 12
Vessel trajectory during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d; T
u
=
2
s.Figure 13
Time dependence of the heading angle during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.14
r
a
d; T
u
=
2
s.We note that in all the cases presented, the control error is close to zero (in the simulation, the first significant digit is in fourth place after the decimal point). By control error, we mean the difference between specified and attained values of the course.Graphs in Figures14 and 15 show the results of vessel motion simulation under the influence of composite control by (45), (52), (54).Figure 14
Vessel trajectory during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.4
r
a
d; T
u
=
0.1
s, F
6
=
10, F
7
=
0.01, z
7
∗
=
0.01
sign
⁡
z
6
∗.Figure 15
Time dependence of the heading angle during rotation towards the return course for valuesz
3
∗
=
0.005
r
a
d
⋅
s
-
1, z
6
∗
=
3.4
r
a
d; T
u
=
0.1
s, F
6
=
10, F
7
=
0.01, z
7
∗
=
0.01
⋅
sign
⁡
z
6
∗.Simulation was performed for two consecutive time intervals with a duration of500
s. During the first time interval, the vessel was controlled by algorithm of (45) with all zero initial conditions except speed V
x
0
=
4
m
·
s
-
1. During the second time interval, the vessel was controlled by algorithm of (47), (52), (54). The initial conditions for the second interval were the values of all coordinates at the end of the first interval, except for z
7
0
=
0. From a technical point of view, this means that Azipod should be brought to the zero (neutral) position at the beginning of second interval. This peculiarity of composite control is connected with the fact that the algorithm of (47), (52), (54) does not contain a feedback from z
7. Therefore, the nonzero initial condition of this coordinate is ignored by the control algorithm. If the initial condition z
7
0 for the second interval is nonzero, an uncompensated static error is formed at the end of the process.As can be seen from the graphs of Figures14 and 15, vessel reached the heading angle of 3.051 radians under the influence of control by (45). Then the vessel made a turn of 3.14 radians to the set course under the influence of control by (52), (54). Herewith, the actual time of the turn is very short—it is about one minute (more precisely, 1193 steps with a duration of 0.05 seconds, that is, 59.65 seconds).It is important to note that the control by (47), (52), (54) explicitly depends on the mass-dimensional and propulsion parameters of the vessel, as well as on parameters h
,
T
u
,
z
3
∗
≡
ω
∗ affecting the shape of the vessel track to the given course. This fact makes it possible to eliminate the problem of smoothing the trajectory (calculating anticipatory points, etc., proposed in [12]).It should be noted that the results concerning vessel control adduced in this paper represent development of methods and algorithms published in [25], regarding a definite implementation of the composite algorithm and refinement of the methods of analytical synthesis of controls as well. A Mathcad program code simulating the vessel motion with a model in the form of (12)-(19) and control by (52) is given in the mentioned paper. This program is easy to modify with (23) and composite control by (45), (52), (54) for the case of the Azipod inertial drive.


## 5. Conclusion

The practical aspects of using the author’s nonlinear systems control method—terminal state method (TSM)—are discussed in the paper. Expressions for control calculating according to TSM with respect to terminal and stabilization (tracking) problems are presented. A special feature of the method is the use of terminal state variables (TSV), whose sense is the forecast of uncontrolled system final state, dependent on the current time and current state. Two examples of the method application are provided. The first one is a terminal state control of a two-wheeled double-track mobile robot. The second one is a control of the course of a small-capacity vessel-catamaran with Azipods during the “strong” maneuvers performing, i.e., maneuvers at large rudder angles. Two TSM algorithms are proposed for the course control. The first one is based on a quadratic criterion for deviations from a given course. The second one is based on the criterion for tracking the vessel angular velocity in relation to its center of mass. The first algorithm is proven to be simpler from a computational point of view, but gives a noticeable static error. The second algorithm is more complicated but has a zero static error. A composite algorithm, consisting of a consecutive implementation of the first and second algorithms, combines the advantages of both algorithms. Mathematical models of two control objects and expressions for control are presented. A routine of obtaining expressions for control actions is described. Results of computer simulation in the form of graphs of time dependencies for robot and vessel current coordinates, their trajectories, and robot control actions are shown.



---

*Source: 4868791-2018-12-19.xml*