*Article* **Cytokinin Modulates Cellular Trafficking and the Cytoskeleton, Enhancing Defense Responses**

**Lorena Pizarro 1,\*, Daniela Munoz 1, Iftah Marash 2,3, Rupali Gupta 2 , Gautam Anand 2, Meirav Leibman-Markus 2 and Maya Bar 2,\***

**Abstract:** The plant hormone cytokinin (CK) plays central roles in plant development and throughout plant life. The perception of CKs initiating their signaling cascade is mediated by histidine kinase receptors (AHKs). Traditionally thought to be perceived mostly at the endoplasmic reticulum (ER) due to receptor localization, CK was recently reported to be perceived at the plasma membrane (PM), with CK and its AHK receptors being trafficked between the PM and the ER. Some of the downstream mechanisms CK employs to regulate developmental processes are unknown. A seminal report in this field demonstrated that CK regulates auxin-mediated lateral root organogenesis by regulating the endocytic recycling of the auxin carrier PIN1, but since then, few works have addressed this issue. Modulation of the cellular cytoskeleton and trafficking could potentially be a mechanism executing responses downstream of CK signaling. We recently reported that CK affects the trafficking of the pattern recognition receptor LeEIX2, influencing the resultant defense output. We have also recently found that CK affects cellular trafficking and the actin cytoskeleton in fungi. In this work, we take an in-depth look at the effects of CK on cellular trafficking and on the actin cytoskeleton in plant cells. We find that CK influences the actin cytoskeleton and endomembrane compartments, both in the context of defense signaling—where CK acts to amplify the signal—as well as in steady state. We show that CK affects the distribution of FLS2, increasing its presence in the plasma membrane. Furthermore, CK enhances the cellular response to flg22, and flg22 sensing activates the CK response. Our results are in agreement with what we previously reported for fungi, suggesting a fundamental role for CK in regulating cellular integrity and trafficking as a mechanism for controlling and executing CK-mediated processes.

**Keywords:** cytokinin; endocytosis; cytoskeleton; actin; plant immunity; induced resistance

### **1. Introduction**

The plant hormone cytokinin (CK) is central to plant life, regulating many processes including embryogenesis, cell division and differentiation, stem cell maintenance, growth and branching of roots and shoots, leaf senescence, nutrient balance, and stress tolerance. CHASE domain-containing AHK (Arabidopsis histidine kinase) receptors, which are located primarily in the endoplasmic reticulum (ER), perceive CK and commence the resultant signaling cascades, which mediate CK function [1–3]. Based on the subcellular localization of AHKs, the endoplasmic reticulum (ER) has been considered to be the primary CK perception site [4,5]. Upon perception through AHKs, AHPs (Arabidopsis histidine phosphotransferases) transfer a phosphate from AHKs to downstream type B ARRs (Arabidopsis response regulators), initiating transcriptional reprogramming as a

**Citation:** Pizarro, L.; Munoz, D.; Marash, I.; Gupta, R.; Anand, G.; Leibman-Markus, M.; Bar, M. Cytokinin Modulates Cellular Trafficking and the Cytoskeleton, Enhancing Defense Responses. *Cells* **2021**, *10*, 1634. https://doi.org/ 10.3390/cells10071634

Academic Editor: Suleyman Allakhverdiev

Received: 31 May 2021 Accepted: 23 June 2021 Published: 29 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

result of CK perception [6]. CK signaling can also start at the plasma membrane (PM) [7]. Extracellular CKs were shown to bind to cell-surface receptors and elicit signaling [8]. ER-localized CK receptors were also demonstrated to be transported to the PM [9]. In that work, AHK4 was shown to reside in the PM and in vesicles, in addition to the ER. Moreover, CK was shown to accumulate in BFA compartments, attesting to AHKs being trafficked between the ER and PM on endocytic vesicles. Thus, initiation sites of cytokinin perception occur at both the PM and the ER, and CK receptors undergo endosomal trafficking between these locales. It has been suggested that the dual CK response locales may contribute to the flexibility of the CK response [8].

CK regulates some developmental processes through crosstalk with the auxin pathway. CK was previously shown to regulate the transcription of auxin pathway genes [10–12]. Importantly, CK was also demonstrated to modulate auxin activity by regulating the endocytic trafficking of the auxin carrier PIN1 [13]. CK regulates the recycling of the auxin efflux carrier PIN1 to the plasma membrane, by directing it to the vacuole for degradation, in a rapid regulatory pathway that does not require transcriptional reprogramming [13]. Thus, through the endocytic redirection of PIN1, CK controls auxin fluxes that are required for lateral root organogenesis, and it was postulated that this endocytic regulation may also be involved in other CK-regulated developmental processes. However, CK was found not to influence the stability of two additional membranal auxin carrier proteins, prompting the authors to suggest that CK affects trafficking in a protein-specific manner [13].

CKs were shown to mediate disease resistance through induction of host immunity. Transgenic plants with high CK levels had increased resistance to *Pseudomonas syringae* [14], while transgenes with low levels of CK or Arabidopsis histidine kinase (AHK) receptor mutants displayed enhanced pathogen susceptibility [14,15]. CKs were also found to mediate enhanced resistance in tobacco [16]. In Arabidopsis, it was suggested that CKmediated resistance depends on salicylic acid (SA) [14], and that CK signaling enhances SA-mediated immunity [17]. In tobacco, an SA-independent, phytoalexin-dependent mechanism was reported [16]. We recently reported that CK induces systemic immunity in tomato, promoting resistance to fungal and bacterial pathogens in an SA- and ethylene (ET)-dependent mechanism [18,19]. Until recently, other than the evidence that SA is required for CK-induced immunity [14,15,17], and that phytoalexins are involved [16], no additional cellular mechanisms related to CK-mediated immunity had been reported.

Pattern recognition receptors (PRRs) are the first line of defense and immune activation in plant cells. Since CK was shown to modulate endocytic trafficking of PIN1 in its regulation of auxin [13], in our recent work, we examined whether CK might also affect the trafficking of immune receptors as a possible mechanism for promoting immune responses and disease resistance. We were able to show that CK modulates the cellular trafficking of the PRR LeEIX2, which mediates immune responses to Xyn11 family xylanases [18]. We found that CK enhances both the endosomal presence and the vesicular size of LeEIX2 endosomes, without affecting the total cellular content of the proteins [18]. Furthermore, asking whether this CK-mediated enhancement of PRR endosomal presence acts as a mechanism for increased disease resistance, we used a SlPRA1A-overexpressing line, which has a decreased presence of receptor-like protein (RLP)-type PRRs in the cell plasma membrane, due to receptor degradation [20]. Indeed, we found that in plants overexpressing SlPRA1A, which have decreased levels of RLP-type PRRs, CK-mediated disease resistance to *Botrytis cinerea* is compromised [18].

Recently, we have found that CK directly affects the growth, development, and virulence of fungal plant pathogens [21]. *Botrytis cinerea* (*Bc*) growth, sporulation, and spore germination were all inhibited by CK, with some fungal developmental processes being inhibited by plant physiological CK concentrations [21]. We found similar effects in a variety of plant pathogenic fungi. We also found that CK affects both budding and fission yeast in a similar manner. Transcriptome profiling revealed that cell cycle and DNA replication, cytoskeleton integrity, and endocytosis are all inhibited in *Bc* by CK [21]. We directly confirmed these results, demonstrating that CK affects the cell

cycle and DNA replication, cytoskeleton distribution, and cellular trafficking in fungi and yeasts [21]. CK interfered with actin localization in growing hyphae and lowered endocytic rates and endomembrane compartment sizes, likely underlying the reduced growth and development we observed [21].

In this work, we take an in-depth look at the effects of CK on cellular trafficking and on the cytoskeleton in plant cells, both in the context of plant immunity and in general. We have found that in addition to the xylanase receptor-like protein (RLP) LeEIX2, CK also affects the distribution of the flagellin receptor-like kinase (RLK) flagellin-sensing 2 (FLS2), increasing its presence in the plasma membrane. FLS2, first characterized in Arabidopsis, acts as the PRR for the bacterial PAMP (pathogen-associated molecular pattern) flagellin [22] in several plant species [23–25]. Furthermore, CK enhances the cellular response to the 22-amino-acid, flagellin-derived peptide flg22, and flg22 sensing activates the CK response. Examining cellular trafficking compartments and the cytoskeleton in general, we also show that CK affects endosome distribution and increases the amount of endomembrane compartments. CK also caused disorganization and reduction in actin filaments, but not in tubulin. Our results are in agreement with what we previously reported for fungi, suggesting a fundamental role for CK in regulating cellular integrity and trafficking as a mechanism for controlling and executing CK-mediated processes.

### **2. Results**

### *2.1. CK Affects PRR Trafficking*

We previously reported that CK modulates the trafficking of the PRR LeEIX2 and enhances its defense signaling, and that PRR-RLPs are required to achieve CK-mediated disease resistance [18]. To examine whether this is a general effect on PRRs, we tested whether CK influences the signaling of the LRR-receptor-like kinase (RLK)-type PRR FLS2. FLS2 is very different from LeEIX2 in its structure, specificity, and regulation [20,22,23]. Figure 1 demonstrates that CK affects cellular distribution of FLS2, increasing its presence at the PM.

### *2.2. CK Response Is Activated during Immunity Elicitation*

We and others have previously demonstrated that pathogenesis processes can activate the CK pathway [17–19], raising the possibility that CK pathway activation might be required for pathogen resistance. Interestingly, we previously found that *B. cinerea* infection activates the synthetic CK response promoter TCSv2 [18]. To test whether defense elicitation activates the CK pathway, we examined whether TCSv2 could respond to flg22 elicitation. As seen in Figure 2, flg22 treatment in Arabidopsis roots significantly enhanced TCSv2 driven Venus expression, indicating that flg22 activates the CK response.

Our previous work has shown that in addition to affecting LeEIX2 trafficking, CK enhances defense signals elicited by its ligand, EIX [18]. Since we observed that CK also affected FLS2 cellular distribution, increasing its presence at the PM, we tested whether CK also enhanced defense signaling elicited by flg22. As seen in Figure 3, CK strongly enhances flg22 elicited reactive oxygen species (ROS) production. CK also has a small but significant inductive effect on ROS production on its own (Figure 3B). The CK trans-zeatin was previously reported to induce ROS production in Arabidopsis guard cells [26].

**Figure 1.** Cytokinin increases the membranal presence of the PRR FLS2. *N. benthamiana* epidermal cells transiently expressing FLS2-GFP were mock treated, or treated with cytokinin (CK, 100 μM 6 benzylaminopurine), for 4 h, and subsequently treated with flg22 (30 min). Membranal plane images (1 μm) were captured in 3 experiments, and images were analyzed using Fiji-ImageJ. (**A**) Representative confocal microscopy images. Contrast was uniformly adjusted; scale bar = 10 μm. (**B**) Twelve images were analyzed. Data are presented as boxplots with inner quartile ranges (boxes), outer quartile ranges (whiskers), and medians (lines in boxes), all points displayed. Letters represent statistical significance in one-way ANOVA with Dunnett's post-hoc test (*p* < 0.035).

**Figure 2.** flg22 activates cytokinin response. Transgenic *A. thaliana* roots stably expressing the cytokinin (CK) response marker *pTCSv2::3XVENUS* were treated with mock (1 mM NaOH), CK (100 μM 6-benzylaminopurine), or flg22 (1 μM), for 12 h. Images were captured in 3 experiments using at least 15 plants per treatment, and analyzed using Fiji-ImageJ; scale bar = 100 μm. (**A**) Representative confocal microscopy images. Contrast was uniformly adjusted; scale bar = 100 μm. (**B**) Images from at least 15 plants were analyzed. Mean fluorescence intensity for each nucleus in each image was quantified, N > 770 nuclei. Data are presented as violin plots—solid gray lines indicate medians, and dotted lines indicate quartile ranges. Letters represent statistical significance in a Kruskal–Wallis test with Dunn's multiple comparisons post-hoc test (*p* < 0.0001).

**Figure 3.** Cytokinin enhances flg22-mediated defense. *S. lycopersicum* cv. M82 leaves were treated for 12 h with mock or with cytokinin (CK, 100 μM 6-benzylaminopurine). After 12 h, leaf discs were harvested and treated with flg22 (1 μM). Reactive oxygen species (ROS) production was measured every 2.5 min for 35 min, using the HRP-luminol method. (**A**) Kinetics of ROS accumulation are plotted; RLU = relative light units; N = 24. (**B**) Total ROS production, area under the graph in A. Average ± SEM is shown, N = 24. Different letters represent statistically significant differences in one-way ANOVA, *p* < 0.0001, with a Dunnett's post-hoc test, *p* < 0.044.

### *2.3. CK Treatment Modulates Cellular Trafficking*

As a developmental hormone, CK regulates various processes that require growth and rapid membrane modeling. These processes could be controlled through cellular trafficking. We have recently found that in fungi, CK attenuates cellular trafficking and causes mislocalization in the actin cytoskeleton [21]. Interestingly, investigation of the effect of CK on the cellular trafficking of the auxin carrier PIN1 showed that CK does not affect the cellular localization of two additional auxin carriers [13].

To test whether the effect of CK on cellular trafficking is general or specific to defense receptors, we examined the distribution of ARA6 and FM-4-64 endosomes following CK treatment. As seen in Figure 4, CK treatment significantly increases the density and size of both FM4-64 compartments (Figure 4A–C) and ARA6 endosomes (Figure 4A,D,E). ARA6 is considered to be a late endosome marker [27], while FM4-64 is present in all endomembrane compartments at the timepoint tested [28].

**Figure 4.** Cytokinin increases endosomal size and density. Transgenic *A. thaliana* root roots stably expressing the endosomal marker ARA6-GFP were grown vertically on 0.5X MS-agar plates, mock or supplemented with cytokinin (CK, 6-benzylaminopurine) at indicated concentrations, and stained with 5 μM FM-4-64 for 5 min on ice. Images were captured in 4 experiments using at least 10 plants per treatment, and analyzed using Fiji-ImageJ. (**A**) Representative confocal microscopy images. Contrast was uniformly adjusted; scale bar = 10 μm. (**B**–**D**) A total of at least 40 cells from at least 10 images per treatment were captured in 4 separate experiments and analyzed. Compartment size and density for FM4-64 (**B**,**C**) and ARA6-GFP (**D**,**E**) are presented as mean ± SEM. Letters represent statistically significant differences in a one-way ANOVA with Tukey's post-hoc test. (**B**) N > 425, *p* < 0.0076. (**C**) N > 40, *p* < 0.0001. (**D**) N > 500, *p* < 0.0037. (**E**) N > 60, *p* < 0.0001.

### *2.4. CK Treatment Affects the Cellular Cytoskeleton*

Growth requires changes to both the composition and the orientation of the cytoskeleton [29]. Furthermore, observed effects on vesicular trafficking can relate to the cytoskeleton, since trafficking relies on an intact cytoskeleton [30]. We previously demonstrated that CK can affect F-actin distribution in fungi [21]. We therefore examined whether CK treatment could affect the plant cytoskeleton. Figure 5 shows that CK affects actin distribution, reducing the quantity and organization of actin in the cell (Figure 5A–C). CK did not affect tubulin (Figure 5D–F).

**Figure 5.** Cytokinin decreases actin organization and density. *Nicotiana benthamiana* epidermal cells transiently expressing Actin-mCherry (**A**–**C**) or mCherry-Tubulin (**D**–**F**) were mock treated, or treated with cytokinin (CK, 100 μM 6-benzylaminopurine) for 4 h. Membranal plane images (1 μm) were captured in 3 experiments, and images were analyzed using Fiji-ImageJ. (**A**,**D**) Representative confocal microscopy images. Contrast was uniformly adjusted; scale bar = 10 μm. (**B**,**C**,**E**,**F**) 25 images per treatment were analyzed. Box plots represent inner quartile ranges (boxes), outer quartile ranges (whiskers), and medians (lines in boxes), all points shown. Asterisks represent statistically significant differences in an unpaired two-tailed *t*-test with Welch's correction, N = 25. (**B**) \*\* *p* = 0.0031. (**C**) \* *p* < 0.028. (**E**,**F**) ns = not significant.

### **3. Discussion**

CK mediates plant immunity and disease resistance [14–19]. Traditionally a plant growth hormone with various developmental roles, CK supports juvenility and delays senescence [31,32], indicating that perhaps commonalities beyond the execution of the cell death program may exist between age-related cell death and pathogen-related cell death, both of which CK can delay or prevent.

In addition to its role in controlling auxin distribution by regulating its trafficking [13], we recently reported that CK can influence cellular trafficking in plants—where it enhances defense signaling of EIX through modulation of the trafficking of its PRR LeEIX2 [18]—and in fungi, where CK treatment lowers growth and development by attenuating trafficking [21]. Furthermore, we observed that CK-mediated disease resistance was lost in a transgenic line with reduced PRR expression [18]. These reports suggest that control of

cellular trafficking could be a general mechanism by which CK exerts its varied effects on development and defense.

In this work, we affirmed the role we previously reported for CK in controlling receptor trafficking, demonstrating that CK enhances flagellin signaling (Figure 2), possibly by increasing the presence of its receptor—FLS2—in the membrane (Figure 1), from which it may transmit part of its defense signals [33]. Interestingly, though ligand-activated FLS2 is endocytosed after about 80–90 min [34], ROS responses occur within 5 min (Figure 3), suggesting that an increase in FLS2 presence at the membrane could underlie the increase in ROS observed upon CK treatment. Furthermore, we observed that flagellin activates CK response (Figure 2). Flagellin was previously reported to increase expression of the CK receptor AHK3 and the CK response regulator ARR2 in Arabidopsis guard cells [26]. We also previously observed a similar mode of action—activation of the CK pathway—in the pathogenesis of *B. cinerea* [18] and bacterial pathogens [19].

The fact that CK enhances signals emanating from both the fungal response RLP LeEIX2 [18] and the bacterial response RLK FLS2 (Figures 1 and 2), and augments defense elicited by both of their respective ligands, EIX [18]—which is derived from *Trichoderma*, which mostly signals through the jasmonic acid pathway [35], and flg22 (Figure 2), which mostly signals through the salicylic acid pathway [36]—indicates that the modulation of the cellular distribution of PRRs is likely a general mechanism by which CK affects cellular defense.

Our recent report [21] observed attenuation of cellular trafficking and disorganization of the actin cytoskeleton in response to CK treatment in *B. cinerea*, suggesting that perhaps control of cellular trafficking and cytoskeletal integrity could be a general mechanism employed to execute CK-mediated responses in plants as well. We examined endocytosis and actin distribution in response to CK treatment in Arabidopsis and *N. benthamiana*. We found that CK significantly increased the size and density of endomembrane compartments, as observed with both the general FM-4-64 dye and the late-endosome-specific ARA6 (Figure 4). Interestingly, the effects of CK on endomembrane compartments did not appear to be dose-dependent in most examined cases, with a cutoff of 1 μM observed, though this effect could be specific to 6-BAP in Arabidopsis, and testing of additional CK concentrations is required in order to definitively make this assertion.

Taken together, our results indicate that one of the early mechanisms for the activation of CK-mediated plant immunity is likely through cellular trafficking. CK may enhance defense signals, both by specifically increasing PRR occupancy in the cellular compartments required for signaling, and by increasing cellular trafficking compartments in general. flg22 is known to induce the entry of FLS2 to endosomes. flg22-induced FLS2-GFP endosomes were shown to co-localize to ARA6 compartments, though there was no difference in the amount of ARA6 endosomes at time points between 30 min and 105 min after flg22 application [37]. This could support the notion that the increased immunity observed upon flg22 and CK co-treatment (Figure 3) stems both from the increased presence of FLS2 in signaling compartments due to flg22 activation, and from a general increase in the available amount of endomembrane compartments, due to CK-mediated activity.

We found that CK inhibits the actin cytoskeleton, though not the tubulin cytoskeleton, in plant cells (Figure 5). We observed similar inhibition of the actin cytoskeleton in *B. cinerea*, underlying the mechanism for growth and development inhibition in fungi [21]. Actin inhibition likely also underlies CK-mediated growth inhibition observed in plants, as exemplified in the inhibition of root growth [38], and also perhaps in processes that maintain plant juvenility, e.g., inhibition of changes to cellular trafficking or degradation of the cytoskeleton required for senescence and cell death. The cytoskeleton is known to be an important component of cellular trafficking. The effect of CK on the actin cytoskeleton could result in altered trafficking rates and paths for endomembrane compartments which could, in turn, affect the distribution of defense-related cargo and result in altered defense signaling.

Interestingly, our findings on the effects of CK on cellular trafficking differ somewhat between plants (Figure 4) and fungi [21]. While in both cases we observed a decrease in both the content and organization of actin filaments, in fungi—in contrast to plants—we observed a decrease in endomembrane compartments upon CK treatment. These differences could relate to differences in both the assay methodologies and the cell biology attributes of different species. Notably, the endocytosis pathway was significantly downregulated upon CK treatment in *B. cinerea* [21].

In summary, we found that CK regulates actin distribution, endocytosis, and PRR trafficking. A model summarizing this work is provided in Figure 6. This is an interesting starting point for future work, which will elucidate how the dual role of CK in growth and defense is regulated, and which additional proteins and feedback circuits are involved, perhaps also elucidating how some of these functionalities evolved in both plants and fungi.

**Figure 6.** Model for the effects of cytokinin on plant cell compartments and defense signaling. Cytokinin (CK) increases endosomal compartment size and density, and decreases actin content and organization. CK promotes defense signaling mediated by both FLS2 and LeEIX2, resulting in increases in reactive oxygen species (ROS) in response to flg22 and EIX, as well as increases in ethylene and electrolyte leakage in response to EIX. CK also primes ROS on its own, and induces gene expression. LeEIX2 and FLS2 elicitor activation induces CK signaling pathways, increasing CK levels and/or activating CK-dependent gene expression. Model created with BioRender.com, accessed on 31 May 2021.

### **4. Materials and Methods**

### *4.1. Plant Growth Conditions*

*Nicotiana benthamiana* and *Solanum lycopersicum* cv. M82 were grown from seeds in soil (Green Mix; Even Ari, Ashdod, Israel) in a growth chamber at 24 ◦C, under long-day conditions (16 h:8 h, light:dark).

*A. thaliana* cv. Columbia transgenic plants overexpressing *pTCS::3XVENUS* [39] or *pARA6::ARA6-Venus* [27] were germinated on 1/2 MS (Duchefa M0225) agar plates containing 50 ug/mL kanamycin, and grown upright in a growth chamber at 22 ◦C under short-day conditions (8 h:16 h, light:dark).

### *4.2. Transient Expression*

Four-week-old *N. benthamiana* plant leaves were abaxially infiltrated with a needleless syringe, with an *Agrobacterium tumefaciens* (strain GV3101) harboring the following previously described constructs: *pAtFLS2::AtFLS2-3xmyc-GFP* [40] (Addgene plasmid#86157; http://n2t.net/addgene:86157, accessed on 31 May 2021; RRID:Addgene\_86157), *pCMU-ACTLr* (AtUBQ10::F-actin LifeAct-mCherry; Addgene plasmid #61193; http://n2t.net/ addgene:61193, accessed on 31 May 2021; RRID:Addgene\_61193) [41], and *pCMU-MTUBr* (AtUBQ10::mCherry-MAP4-MBD; Addgene plasmid #61196; http://n2t.net/addgene, accessed on 31 May 2021; RRID:Addgene\_61196) [41]. Treatments were applied 40 h after infiltration. Transgenic lines and constructs used in this study are detailed in Table 1.

**Table 1.** Transgenic lines and constructs used in this study.

### *4.3. Chemical Treatments*

The CK 6-benzylaminopurine (6-BAP, Benzyladenine, sigma) was dissolved in 10 mM NaOH. *pARA6:ARA6-Venus* Arabidopsis plants were grown in 0.5X MS-agar media and transferred 3 days after germination to 0.5X MS-agar media supplemented with 6-BAP at indicated concentrations. *pTCS::3XVENUS* Arabidopsis plants, 5 days after germination, were incubated for 16 h in 0.5X MS media supplemented with 6-BAP at indicated concentrations (Figures 2 and 4). In the transient expression assays, 6-BAP water solutions were infiltrated into *N. benthamiana* leaves, and images were captured 6 h after infiltration (Figures 1 and 5). flg22 (1 μM, PhytoTech Labs #P6622) was dissolved in DMSO and applied as described for CK.

FM-4-64 (Invitrogen # T13320) was dissolved in water and applied to Arabidopsis roots in solution at 5 μM. Arabidopsis roots were incubated with this solution for 5 min on ice prior to live-cell imaging.

### *4.4. Confocal Microscopy*

Confocal microscope images were acquired using a Zeiss LSM780 confocal microscope with Objective C-Apochromat 40 ×/1.2 W Corr M27 (Figure 1, Figure 2, and Figure 5) or Objective C-Apochromat 63×/1.2 W Corr (Figure 4). GFP was excited using a laser of 488 nm (5% power), and emission was collected in the range of 493–535 nm. mCherry fluorescence was excited using a laser of 561 nm (3% power), and emission was collected in the range of 588–641 nm. FM4-64 was excited using a laser of 514 nm (3% power), and emission was collected in the range of 650–750 nm. Images of 8 bits were acquired using a pixel dwell time of 1.27 μs, pixel averaging of 4, and pinhole of 1 airy unit. Image analysis was performed using Fiji-ImageJ with the raw images [42]. Endosome density and size measurements were conducted with the 3D Object counter tool, pixel intensity was

measured using the measurement analysis tool, cytoskeleton organization was measured using the skew tool, and occupancy was assessed using the 3D Object tool [42].

### *4.5. Reactive Oxygen Species Burst Measurement*

ROS measurement was performed as previously described [43]. Leaf disks of 0.5 cm in diameter were harvested from leaves 4 to 6 of 5–6-week-old M82 plants, treated for 12 h with mock (1 mM NaOH) or CK (100 μM 6-benzylaminopurine). Disks were floated in a white 96-well plate (SPL Life Sciences, Korea) containing 250 μL distilled water for 4–6 h at room temperature. After incubation, water was removed, and an ROS measurement reaction containing 1 μM flg22 was added. Light emission was immediately measured using a luminometer (Spark, Teacan, SL).

### *4.6. Data Analysis*

Differences between two groups were analyzed for statistical significance using an unpaired two-tailed *t*-test, with Welch's correction where samples were found to have unequal variances, or with Holm–Sidak correction where multiple *t*-tests were applied. Differences among three groups or more were analyzed for statistical significance with a one-way ANOVA, or with a Kruskal–Wallis test. For ANOVA analyses, regular ANOVA was used for groups with equal variances, and Welch's ANOVA for groups with unequal variances. When a significant result for a group was returned, differences between the means of different samples in the group were assessed using a post-hoc test. Tukey's test was employed for ANOVA conducted on samples with equal variances, Dunnett's test was employed for ANOVA conducted on samples with unequal variances, and Dunn's multiple comparisons post-hoc test was employed for Kruskal–Wallis tests. All statistical analyses were conducted using Prism9.

**Author Contributions:** Conceptualization, L.P., M.L.-M., and M.B.; methodology, L.P. and M.B.; validation, L.P., D.M., I.M., R.G., G.A., and M.L.-M.; formal analysis, L.P., D.M., I.M., M.L.-M., and M.B.; writing—original draft preparation, M.B.; writing—review and editing, L.P. and M.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was partially supported by grant PAI 77190027 to LP. The APC was funded by M.B.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The authors declare that the data supporting the findings of this study are available within the paper. Raw data are available from the corresponding author upon reasonable request.

**Acknowledgments:** We thank Takashi Ueda for providing seeds of *A. thaliana* col. *pARA6:ARA6- Venus,* Silke Robatzek for *pAtFLS2:FLS2-3xmyc-GFP*, Maria Harrison for *pCMU-ACTLr* and *pCMU-MTUBr*, and Adi Avni for assistance with microscopy. We thank all members of the Pizarro and Bar groups for helpful discussions and support.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


*Article* **Optimal Damping Concept Implementation for Marine Vessels' Tracking Control**

**Evgeny I. Veremey**

Computer Applications and Systems Department, Saint-Petersburg State University, Universitetskii Prospekt 35, Petergof, 198504 Saint Petersburg, Russia<EMAIL_ADDRESS>
**Abstract:** This work presents the results of studies related to the design of stabilizing feedback connections for marine vessels moving along initially given trajectories. As is known, in mathematical formalization, this question leads to a problem of tracking control synthesis for nonlinear and nonautonomous plants. To provide desirable stability and performance features of the closed-loop systems to be synthesized, it is appropriate to use an optimization approach. Unlike the known synthesis methods, which are usually used within the framework of this approach, it is proposed to implement the optimal damping concept first developed by V.I. Zubov in the early 60s of the last century. Modern interpretation of this concept allows constructing numerically effective procedures of control law synthesis taking into account its applicability in a real-time regime. Central attention is focused on the questions connected with practical adaptation of the optimal damping methods for marine control systems. The operability and effectiveness of the proposed approach are illustrated by a practical example of tracking control design.

**Keywords:** marine vessel; tracking controller; stability; functional; optimal damping

#### **1. Introduction**

The nonlinear tracking control of modern marine vessels is one of the most practically significant and theoretically considerable problems in the area of automatic control analysis and design. In particular, tracking control systems are widely used in different branches such as hydrography, inspection of marine constructions, wreck investigation, underwater cable laying, and so on [1,2]. Central theoretical and practical background of tracking control for various moving plants is presented in [3–5] and other fundamental works.

Various issues associated with the design of tracking controllers for marine surface vessels have already been extensively researched and presented in numerous publications (for example, [1,2,6–16]). To evaluate the state of the art in marine tracking control, let us address some modern works presenting this direction of research.

Currently, it is possible to use various ideas to design nonlinear tracking control laws that are reflected in numerous publications, for example, [7,8,13–16]. However, let us note that the mentioned works are not directly oriented to the application of the optimization technique. This makes it difficult to provide the desired dynamic features of the closedloop connections. Now, it seems to be quite evident that the most effective analytical and numerical tool for feedback connections design is the optimization approach. Several aspects of nonlinear tracking control optimization technique are presented in multitudinous scientific publications, including such popular monographs as [4,5,17–21]. As for the simplest stabilization problem, the autopilots with multipurpose structures of optimized control laws are discussed in detail in [22–26].

The sliding mode control technique for marine control applications is discussed in [1,2,6,8]. This direction seems to be quite constructive, but poorly applicable, since it leads to intensive wear of the actuators.

**Citation:** Veremey, E.I. Optimal Damping Concept Implementation for Marine Vessels' Tracking Control. *J. Mar. Sci. Eng.* **2021**, *9*, 45. https:// doi.org/10.3390/jmse9010045

Received: 26 November 2020 Accepted: 29 December 2020 Published: 4 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

As for the model predictive control (MPC) approach [9], its most significant disadvantage is the large dimension of the minimization problem that is solved at each step of the control process.

Notably, the complexity of this problem is vast because of the many dynamic requirements, restrictions, and conditions that must be satisfied by the chosen control actions.

It should be noted that many scientific works devoted to the tracking control for marine vessels use linear time invariant models of their motion. However, such models are not quite adequate for the problems of deep maneuvering control in angle and positional dynamic variables. Respectively, one of the most important practical difficulties requiring consideration in the design process is the account of nonlinearity and non-autonomy of the control plant model. In most cases, this problem is a source of dynamic instability and poor performance for various systems that were designed based only on linear approximations.

As for the aforementioned optimization approach, its advantages are determined by the flexibility and convenience of modern optimization methods with respect to the relevant practical demands for control design implementation. Certain analytical and numerical methods are used now to compute the optimal controllers for nonlinear and non-autonomous systems subject to various given performance indices. Nevertheless, there is no saying that the optimization approach is recognized overall as a universal instrument to be put into practice for marine tracking controllers design. This can be explained by the presence of some disadvantages connected with computational troubles. Therefore, there exists a vital necessity to develop persistently analytical and numerical methods of control laws design based on optimization ideology adapting to the specific problems for various marine applications.

At present time, numerous approaches are used for a practical solution of these problems [1–12,17,18]. Usually, they are based on Pontryagin's maximum principle, on Bellman's dynamic programming principle, on finite-dimensional approximation in the range of model predictive control (MPC) technique, etc. Unfortunately, all these approaches are connected with the huge extent of calculations that essentially impedes their implementation both for laboratory design activity and real-time regimes of control.

The existence of numerical difficulties motivates us to use other approaches that allow avoiding the aforementioned shortcomings. This work focuses on a different concept that can be applied to design tracking controllers using the theory of optimal damping (OD). This theory, which was first proposed and developed by V.I. Zubov in his works [19–21], provides effective analytical and numerical methods for control calculations with essentially reduced computational consumptions with respect to classical techniques. We believe that this theory was ahead of its time and was undeservedly underutilized for practical control problem solving. This work is one of the attempts to overcome this omission, taking into account the impressive development of modern computer technologies.

In this article, special attention is paid to the control of marine vessels in terms of the forward speed and heading angle. We are considering the regime of the acceleration in order to achieve the specified forward speed with one-time turn along the heading. To achieve desirable stability and performance features of the reference motion, the correspondent tracking controllers were designed based on the OD technique.

The main contribution of this paper is determined by the following statements. First, we propose to use the OD concept to design tracking controllers for marine vessel speed and heading. This has not been the case before. Second, we discuss a new methodology for selecting the functional to be damped, taking into account the desirable features of the closed-loop system in the range of the optimization technique. Hereby, the choice of this functional as the basis is argued by the guarantee asymptotic stability and the desired quality of control processes. Third, we point to the possibility of applying the OD approach to a wide class of nonlinearities in the mathematical model of the vessel. It is noted that this approach can be implemented in real-time regime of a ship's motion. The practical applicability and effectiveness of the proposed technique is illustrated by a controller design for a transport marine ship.

The novelty of the proposed approach with respect to other works lies in the universality and flexibility of proposed nonlinear non-autonomous control laws based on OD computational procedure, which can be implemented in a real-time regime of functioning for marine control plants.

In general, the present study is an extension of the multipurpose approach proposed in [22–27] and developed in [28] with respect to the marine autopilot control laws with the novel structure, taking into account actuators' time delays.

This article is organized as follows. In Section 2, the optimal damping concept for control law synthesis for nonlinear non-autonomous systems is discussed, taking into account certain specific stability and performance requirements for marine control applications. The known background is presented, and the novel ways are proposed to provide tracing controllers synthesis. Section 3.1 is devoted to the OD synthesis problem statement for the forward speed tracking controller and for the tracking autopilot. Central attention is paid to the presentation of mathematical models of the control plant and dynamic requirements for the quality of the closed-loop connection. Section 3.2 presents an exhaustive novel solution for the mentioned synthesis problem based on the optimal damping concept. In Section 3.3, a practical example of tracking controller synthesis is presented to illustrate the applicability and effectiveness of the proposed approach. Finally, Section 4 concludes the article by discussing the overall results of the investigation and indicates how these results can be further developed.

#### **2. Materials and Methods**

As mentioned above, the essence of this paper involves developing an optimal damping technique of tracking control law synthesis for marine vessels with nonlinear and non-autonomous models. In this section, let us first consider the background and some essential features of the OD approach that define the methodological basis of the study.

First of all, let us introduce a commonly used nonlinear robot-like model of the control plant, which represents marine vessel motion for various regimes of its operation [1,2,6]:

$$\begin{cases} \mathbf{M}\dot{\mathbf{v}} + \mathbf{C}(\mathbf{v})\mathbf{v} + \mathbf{D}(\mathbf{v})\mathbf{v} + \mathbf{g}(\eta) = \mathbf{G}\_u \pi + \mathbf{d}, \\ \dot{\eta} = \mathbf{J}(\eta)\mathbf{v}, \end{cases} \tag{1}$$

where vector ν ∈ *Rn* presents velocities defined in a plant-fixed frame and vector η ∈ *Rn* contains position dynamical parameters (displacements and angles) in an Earth-fixed frame. External disturbances and controls are presented by the vectors **d** ∈ *Rn* and τ ∈ *Rm*, respectively. Let us accept that the inertia matrix is positive definite: **M** = **M***T* > 0, the matrix of Coriolis-centripetal terms is skew-symmetrical: **C**(ν) = −**C***T*(ν), and the damping matrix **D**(ν) > 0 is positive definite but non-symmetrical. Vector **g**(η) represents gravitational and buoyancy forces and moments, **J**(η) is the matrix of rotations, and the matrix **G***u* with the constant components reflects controls allocation.

Let us provide a transformation of the body-fixed frame representation (1) to the Earth-fixed one with respect to the vector η. Following [1], this can be done using the following notations:

$$\begin{array}{lcl}\mathbf{M}\_{\eta}(\eta) := \mathbf{J}^{-T}(\eta)\mathbf{M}\mathbf{J}^{-1}(\eta),\\ \mathbf{C}\_{\eta}(\mathbf{v},\eta) := \mathbf{J}^{-T}(\eta)\left[\mathbf{C}(\mathbf{v}) - \mathbf{M}\mathbf{J}^{-1}(\eta)\dot{\mathbf{J}}(\eta)\right]\mathbf{J}^{-1}(\eta),\\ \mathbf{D}\_{\eta}(\mathbf{v},\eta) := \mathbf{J}^{-T}(\eta)\mathbf{D}(\mathbf{v})\mathbf{J}^{-1}(\eta), \; \mathbf{g}\_{\eta}(\eta) := \mathbf{J}^{-T}(\eta)\mathbf{g}(\eta),\\ \mathbf{\tau}\_{\eta} := \mathbf{J}^{-T}(\eta)\mathbf{G}\_{\boldsymbol{\mu}}\mathbf{\tau}, \; \mathbf{d}\_{\eta} = \mathbf{J}^{-T}(\eta)\mathbf{d}.\end{array} \tag{2}$$

In accordance with (2), initial model (1) of the plant takes the form:

$$\ddot{\boldsymbol{\eta}} = -\mathbf{M}\_{\eta}^{-1}(\boldsymbol{\eta}) \Big( (\mathbf{C}\_{\eta}(\mathbf{v}, \boldsymbol{\eta}) + \mathbf{D}\_{\eta}(\mathbf{v}, \boldsymbol{\eta})) \dot{\boldsymbol{\eta}} + \mathbf{g}\_{\eta}(\boldsymbol{\eta}) \ -\boldsymbol{\pi}\_{\eta} - \mathbf{d}\_{\eta} \Big), \ \mathbf{v} = \mathbf{J}^{-1}(\boldsymbol{\eta}) \dot{\boldsymbol{\eta}}.\tag{3}$$

The essence of the tracking control problem is to provide given desirable motion η = η*d*(*t*) of the vessel, using the following state feedback

$$
\mathfrak{r} = \mathfrak{r}(\mathfrak{v}\_{\prime}, \mathfrak{v}, \mathfrak{n}\_{d}(t)),
\tag{4}
$$

which is a nonlinear non-autonomous tracking controller.

Within mathematical formalization, controller (4) must be implemented to provide the zero equilibrium with respect to the tracking error **e**(*t*) := η(*t*) − η*d*(*t*) for the closed-loop system (3), (4), where **d**(*t*) ≡ 0. Naturally, this equilibrium point must be asymptotically stable to guarantee that **e**(*t*) → 0 as *t* → ∞. Let us especially note that the mentioned closed-loop system is nonautonomous, if we have no constant reference motion η*d*(*t*). This gives reasons for us to require the uniform asymptotic stability in global (UGAS) or local (UAS) form. An additional requirement is that the controller (4) provides the desired dynamical features for the closed-loop system (3), (4) under the action of an admissible control τ ∈ *Tu*.

To set the perform of the controller (4) synthesis, we assume that the vector functions η*d*(*t*), ν(*t*) := **J**−1(η*d*(*t*)) **.** η*d*(*t*), and the corresponding τ*d*(*t*) are given. These functions satisfy Equations (1) or (3), i.e., we have:

$$\begin{array}{l} \mathbf{M}\dot{\mathbf{v}}\_{d}(t) + \mathbf{C}(\mathbf{v}\_{d}(t))\mathbf{v}\_{d}(t) + \mathbf{D}(\mathbf{v}\_{d}(t))\mathbf{v}\_{d}(t) + \mathbf{g}(\eta\_{d}(t)) \equiv \mathbf{G}\_{u}\mathbf{r}\_{d}(t),\\ \dot{\mathbf{\eta}}\_{d}(t) \equiv \mathbf{J}(\eta\_{d}(t))\mathbf{v}\_{d}(t). \end{array} \tag{5}$$

Let introduce the following additional notations:

**.**

$$\tilde{\mathbf{x}} := \begin{pmatrix} \tilde{\mathbf{x}}\_1 \\ \mathbf{x}\_2 \end{pmatrix} = \left( \frac{\mathbf{v}}{\eta} \right), \mathbf{f}(\tilde{\mathbf{x}}) := \left( \frac{-\mathbf{M}^{-1}}{\eta} \mathbf{\tilde{L}}^{-1} \underline{\mathbf{C}}(\underline{\mathbf{v}}) + \frac{\mathbf{D}(\mathbf{v})}{\mathbf{J}(\eta)\mathbf{v}} \underline{\mathbf{v}} - \underline{\mathbf{M}}^{-1} \underline{\mathbf{g}}(\underline{\eta}) \right), \mathbf{B} := \left( \frac{\mathbf{M}^{-1} \mathbf{G}\_{\tilde{\mathbf{u}}}}{\mathbf{0}} \underline{\mathbf{v}} \right)$$

which allows us to present Equations (1) and (5) as

$$
\tilde{\mathbf{x}} = \mathbf{f}(\tilde{\mathbf{x}}) + \mathbf{B}\boldsymbol{\tau}\dot{\mathbf{x}}\_d \equiv \mathbf{f}(\mathbf{x}\_d) + \mathbf{B}\boldsymbol{\tau}\_d \tag{6}
$$

supposing that **d**(*t*) ≡ 0.

Let us also consider deflections

$$\mathbf{x} := \tilde{\mathbf{x}} - \mathbf{x}\_d = \begin{pmatrix} \mathbf{x}\_1 \\ \mathbf{x}\_2 \end{pmatrix} := \begin{pmatrix} \mathbf{e}\_\mathbf{v} \\ \mathbf{e} \end{pmatrix} := \begin{pmatrix} \mathbf{e}\_\mathbf{v} \\ \mathbf{e} \end{pmatrix} := \begin{pmatrix} \mathbf{v} - \mathbf{v}\_d \\ \boldsymbol{\eta} - \mathbf{\eta}\_d \end{pmatrix}, \mathbf{u} := \mathbf{r} - \mathbf{r}\_d \tag{7}$$

of the vessel dynamical parameters from the desirable motion.

Then, we can present equations of the vessel in the deflections from the desirable motion. Using notations (7) on the base of (6), we obtain

$$
\dot{\mathbf{x}} = \mathbf{a}(t, \mathbf{x}) + \mathbf{B}\mathbf{u},\tag{8}
$$

where

$$\mathbf{a}(t,\mathbf{x}) = \mathbf{a}(t,\mathbf{e},\mathbf{e}\_{\mathbf{v}}) := \begin{pmatrix} -\mathbf{M}^{-1}[\mathbf{C}(\mathbf{e}\_{\mathbf{v}} + \mathbf{v}\_{d}(t)) + \mathbf{D}(\mathbf{e}\_{\mathbf{v}} + \mathbf{v}\_{d}(t))](\mathbf{e}\_{\mathbf{v}} + \mathbf{v}\_{d}(t)) - \mathbf{M}^{-1}\mathbf{g}(\mathbf{e} + \mathbf{v}\_{d}(t))} \\ -\mathbf{M}^{-1}[\mathbf{C}(\mathbf{e}\_{d}(t)) + \mathbf{D}(\mathbf{v}\_{d}(t))]\mathbf{v}\_{d}(t) - \mathbf{M}^{-1}\mathbf{g}(\mathbf{n}\_{d}(t))} \\ -\mathbf{M}^{-1}[\mathbf{C}(\mathbf{v}\_{d}(t)) + \mathbf{D}(\mathbf{v}\_{d}(t))]\mathbf{v}\_{d}(t) - \mathbf{M}^{-1}\mathbf{g}(\mathbf{n}\_{d}(t))} \end{pmatrix} \tag{9}$$

It is a matter of simple calculations to check that equation (6) has zero equilibrium position, which must be stabilized by the choice of the controller **u** = **u**(*t*, **x**). If this controller is found, then the tracking feedback (4) can be presented as

$$\boldsymbol{\pi} = \boldsymbol{\pi}(t, \mathbf{v}, \boldsymbol{\eta}) = \boldsymbol{\pi}(\boldsymbol{\eta}, \mathbf{v}, \boldsymbol{\eta}\_d) = \mathbf{u}(t, \mathbf{x}) + \boldsymbol{\pi}\_d(t), \; \mathbf{x} := \begin{pmatrix} \mathbf{v} - \mathbf{v}\_d(t) \\ \overline{\mathbf{v}} - \overline{\mathbf{v}\_d(t)} \end{pmatrix}. \tag{10}$$

As for the desirable dynamic features of the closed-loop connection, the most widely used formalized approach is based on minimizing the following integral functional:

$$J = J(t, \mathbf{x}, \mathbf{u}) = \int\_{t\_0}^{\infty} f(\mathbf{r}, \mathbf{x}, \mathbf{u}) \, d\mathbf{r} \,\omega$$

that determines the quality of control processes for a closed-loop system. Here, subintegral function *f* is positive definite, i.e., *f*(*t*, **x**, **u**) ≥ 0 ∀*t* ≥ *t*0, ∀**x**, ∀**u**, *f* = 0 ⇔ **x** = 0, **u** = 0. However, as is well known, there are certain difficulties in directly implementing this approach, which are determined by the huge extent of calculations that essentially impedes their practical implementation.

We propose to overcome the mentioned obstacles using a novel technique based on the OD concept, connected to the OD problem for the synthesis of the control **u**. To state this problem, firstly, let us introduce the functional to be damped as follows:

$$L = L(t, \mathbf{x}, \mathbf{u}) = V(t, \mathbf{x}) + \int\_{t\_0}^{t} f\_0(\mathbf{r}, \mathbf{x}, \mathbf{u}) \, d\mathbf{r},\tag{11}$$

where *V* = *V*(*t*, **x**) is a Lyapunov function candidate, and *f*0 is a positively defined function. Let us additionally accept that the function *V* satisfies the following conditions:

$$
\alpha\_1(\|\mathbf{x}\|) \le V(t, \mathbf{x}) \le \alpha\_2(\|\mathbf{x}\|) \forall \mathbf{x} \in B\_r \subset E^n, \forall t \ge t\_0. \tag{12}
$$

where α1, α2 ∈ *K* are comparison functions [4].

The essence of the OD approach consists of the control generation as a function from the current values of variables *t*, **x** in the form

$$\mathbf{u}\_0 = \mathbf{u}\_0(t, \mathbf{x}) = \operatorname\*{argmin}\_{\mathbf{u} \in \mathcal{U}} \mathcal{W}(t, \mathbf{x}, \mathbf{u}) \tag{13}$$

where *U* ⊂ *Em* is the metric compact set, and *W* is a rate of the functional *L* change along the motions of the plant (8):

$$\mathcal{W}(t, \mathbf{x}, \mathbf{u}) := \left. \frac{dL}{dt} \right|\_{\left(8\right)} = \left. \frac{dV}{dt} \right|\_{\left(8\right)} + f\_0(t, \mathbf{x}, \mathbf{u}) = \frac{\partial V}{\partial t}(t, \mathbf{x}) + \frac{\partial V}{\partial \mathbf{x}}(t, \mathbf{x}) \, \mathbf{a}(t, \mathbf{x}) + \frac{\partial V}{\partial \mathbf{x}}(t, \mathbf{x}) \, \mathbf{B} \mathbf{u} + f\_0(t, \mathbf{x}, \mathbf{u}).$$

In other words, it is necessary to find OD controller (13), using an admissible set *U* ⊂ *Em* such that ∀**u** ∈ *U* : τ*d*(*t*) + **u**(*t*, **x**) ∈ *Tu*, ∀*t* ≥ 0.

$$\left. \frac{d\mathcal{W}}{d\mathbf{u}} \right|\_{\mathbf{u}=\mathbf{u}\_{0}(t,\mathbf{x})} = 0 \quad \Rightarrow \quad \left[ \frac{\partial V}{\partial \mathbf{x}}(t,\mathbf{x}) \, \mathbf{B} + \frac{\partial F}{\partial \mathbf{u}}(t,\mathbf{x},\mathbf{u}) \right]\_{\mathbf{u}=\mathbf{u}\_{0}(t,\mathbf{x})} = 0. \tag{14}$$

Using the necessary condition (14), one can solve the following nonlinear system

$$\mathbf{a}(t, \mathbf{x}, \mathbf{u}) = \mathbf{b}(t, \mathbf{x})$$

for any point (*t*, **x**) with respect to the vector **u**, where **a**(*t*, **x**, **u**) = col[*ai*(*t*, **x**)], **b**(*t*, **x**) = col[*bi*(*t*, **x**)],

$$a\_i(t, \mathbf{x}, \mathbf{u}) := \left[\frac{\partial F}{\partial \mathbf{u}}(t, \mathbf{x}, \mathbf{u})\right]\_{i} \\ \text{ } b\_i(t, \mathbf{x}) := -\left[\frac{\partial V}{\partial \mathbf{x}}(t, \mathbf{x}) \, \mathbf{B}\right]\_{i} \\ \text{ } i = \overline{1, n}.$$

Based on [4,5], it can be shown that if the function *V* = *V*(*t*, **x**)is such that *W*(*t*, **x**, **u**0(*t*, **x**)) ≤ −α3(**x**) ∀**x** ∈ *Br*, ∀*t* ≥ *t*0, where α3 ∈ *K*, then the function *V* is control Lyapunov function (CLF) for the plant (8), and the zero equilibrium for the closed-loop system (8), (13) is UAS.

#### **3. Results**

This section is devoted to a practical implementation of the proposed OD approach to nonlinear tracking controllers design for marine vessels. Particular attention is given to the forward speed tracking control law with initially given reference signal.

#### *3.1. Tracking Control Problem for Marine Vessels*

To consider the problems of tracking control for marine applications, let us accept the following widely used [1,2] nonlinear dynamical model of a marine vessel:

$$\begin{array}{ll}\frac{d\overline{V}\_{\mathbf{x}}}{dt} = T\_{\overline{v}}(\overline{V}\_{\mathbf{x}}, \overline{\mathfrak{L}}, \overline{u}\_{1}) + T\_{h}(\overline{V}\_{\mathbf{x}}, \overline{\mathfrak{L}}),\\\frac{d\overline{V}\_{\mathbf{y}}}{dt} = h\_{2}(\overline{V}\_{\mathbf{x}}, \overline{\mathfrak{L}}), \ \frac{d\overline{w}}{dt} = h\_{3}(\overline{V}\_{\mathbf{x}}, \overline{\mathfrak{L}}),\\\frac{d\overline{w}}{dt} = \overline{w}, \ \frac{d\overline{s}}{dt} = \overline{u}\_{2}.\end{array} \tag{15}$$

Here, the following notations are used: *Vx*, *Vy*, and ω are the projections of the velocity vectors on the axes of a vessel-fixed frame *Oxyz* (Figure 1); ξ := *Vy* ω ϕ δ *T* is the auxiliary vector of dynamical parameters; ϕ is the heading angle, and δ is the vertical rudder deflection. The functions *Tv* and *Th* represent hydrodynamical forces, which are produced by the ship's engine and the water resistance correspondingly.

**Figure 1.** Earth-fixed and vessel-fixed coordinate frames.

The control signals *u*1 and *u*2 must be composed by the automatic control system to be designed. The first of them determines a reference surge speed of the vessel, and the second one presents a reference speed of the rudders' deflections.

Assuming that the number of the screw rotations is proportional to a reference surge speed, the function *Tv*(*Vx*, ξ, *u*1) determining a trust force of the screw can be presented in the form

$$T\_v(\nabla\_x, \mathsf{E}, \pi\_1) \equiv \alpha \pi\_1^2 + \beta (\nabla\_x, \mathsf{E}) \pi\_1 + \widetilde{\gamma} (\nabla\_x, \mathsf{E}),\tag{16}$$

We accept here that α ≡ *const*, and that the variables β and *γ*1 are treated as known functions of the dynamical parameters *Vx*, ξ.

Let us introduce new vector variables **h**(*Vx*, ξ) := *h*2(*Vx*, ξ) *h*1(*Vx*, ξ) ξ3 0 *T* , **b** = 0001 *T* , and, taking into account (16), we can rewrite nonlinear model (1) of the vessel dynamic as follows:

$$\begin{array}{l}\frac{d\nabla\_{\mathbf{x}}}{dt} = \alpha\overline{u}\_{1}^{2} + \beta(\nabla\_{\mathbf{x}}, \overline{\mathbf{\xi}})\overline{u}\_{1} + \gamma(\overline{\nabla}\_{\mathbf{x}}, \overline{\mathbf{\xi}}),\\\frac{d\overline{\mathbf{\xi}}}{dt} = \mathbf{h}(\nabla\_{\mathbf{x}}, \overline{\mathbf{\xi}}) + \mathbf{b}\overline{u}\_{2}.\end{array}\tag{17}$$

Now we specify a certain reference motion *Vx*ρ(*t*), ξρ(*t*), *u*1ρ(*t*), and *u*2ρ(*t*) of the plant (17), satisfying the identities

$$\begin{array}{l}\frac{dV\_{x\rho}}{dt} \equiv \mathfrak{a}u\_{1\rho}^{2} + \mathfrak{k}(V\_{x\rho\prime}\mathfrak{k}\_{\rho})u\_{1\rho} + \gamma(V\_{x\rho\prime}\mathfrak{k}\_{\rho}),\\\frac{d\underline{k}\_{\rho}}{dt} \equiv \mathfrak{h}(V\_{x\rho\prime}\mathfrak{k}\_{\rho}) + \mathfrak{b}u\_{2\rho}.\end{array} \tag{18}$$

Using systems (17) and (18), we can present equations of vessel dynamic in deflections from the desirable reference motion of the form

$$\begin{split} \frac{dV\_x}{dt} + \frac{dV\_{\rm{sp}}}{dt} &= \alpha \left( u\_1 + u\_{1\rho} \right)^2 + \beta \left( V\_x + V\_{x\rho}, \xi + \xi\_\rho \right) \left( u\_1 + u\_{1\rho} \right) + \\ &+ \gamma \left( V\_x + V\_{x\rho}, \xi + \xi\_\rho \right), \\ \frac{d\xi}{dt} + \frac{d\xi\_\rho}{dt} &= \mathbf{h} \left( V\_x + V\_{x\rho}, \xi + \xi\_\rho \right) + \mathbf{b} \left( u\_2 + u\_{2\rho} \right), \end{split} \tag{19}$$

where

$$V\_{\mathbf{x}} := \overline{V}\_{\mathbf{x}} - V\_{\mathbf{x}\boldsymbol{\rho}\_{\prime}} \ \xi := \overline{\xi} - \xi\_{\boldsymbol{\rho}\_{\prime}} \ u\_{1} := \overline{u}\_{1} - u\_{1\boldsymbol{\rho}\_{\prime}} \ u\_{2} := \overline{u}\_{2} - u\_{2\boldsymbol{\rho}\_{\prime}}.\tag{20}$$

One can see that *Vx*ρ = *Vx*ρ(*t*), ξρ = ξρ(*t*), *u*1ρ = *u*1ρ(*t*), and *u*2ρ = *u*2ρ(*t*) are known functions of *t*: using new correspondent notations β1, *γ*1, **h**1, we can rewrite (19) as follows:

$$\begin{array}{l} \frac{dV\_{\mathbf{x}}}{dt} = \alpha u\_1^2 + \beta\_1 (t, V\_{\mathbf{x}\prime} \mathfrak{E}) u\_1 + \gamma\_1 (t, V\_{\mathbf{x}\prime} \mathfrak{E}),\\ \frac{d\mathfrak{E}}{dt} = \mathbf{h}\_1 (t, V\_{\mathbf{x}\prime} \mathfrak{E}) + \mathbf{b} u\_2 \end{array} \tag{21}$$

stating that the resulting system (21) has a zero-equilibrium position.

Let us especially note that, unlike (17), this system is not only nonlinear, but also non-autonomous. This is due to the explicit introduction of the time-dependent reference signals into the vessel dynamics equations.

The purpose of the control design procedure is to construct the following stabilizing feedback controllers in deviations:

$$
\mu\_1 = \mu\_1(t, V\_\mathbf{x}, \mathfrak{E}), \; \mu\_2 = \mu\_2(t, V\_\mathbf{x}, \mathfrak{E}), \tag{22}
$$

such that the zero equilibrium of the closed-loop connection (21), (22) is asymptotically stable. This means that if the motion of the initial plant (15) takes place under the action of tracking controllers of the form

$$
\overline{u}\_1 = u\_{1\rho}(t) + u\_1(t, V\_\mathbf{x}, \mathfrak{E}), \\
\overline{u}\_2 = u\_{2\rho}(t) + u\_2(t, V\_\mathbf{x}, \mathfrak{E}), \tag{23}
$$

starting in some neighborhood of the point {*Vx*ρ(0), ξρ(0)}, then this motion tends to the reference one as *t* becomes infinite.

As for the performance of control processes, they are usually formalized mathematically using certain functionals, which are given on the motions of the closed-loop system (21), (22). Currently, the commonly used approach to design stabilizing controllers (22) is setting and solving the following optimization problem

$$J = J(u\_1, u\_2) \to \min\_{\{u\_1, u\_2\} \in \Omega} \tag{24}$$

based on the integral functional

$$J = J(\mu\_1, \mu\_2) = \int\_{t\_0}^{\infty} f^\*(t, V\_x(t), \xi(t), \mu\_1(t), \mu\_2(t)) \, dt,\tag{25}$$

where *U*1 is the set of admissible pairs {*u*1, *u*2} and subintegral function *f* ∗ is positive definite for all its arguments.

In contrast to the problems (24), (25) with traditional methods of its solving, as mentioned above, we propose to use novel approach based on the OD concept.

Let us especially note that both the solution of the problem (24) and the solution of the OD problem significantly depend on the initially given mathematical model (21) of a marine vessel. Naturally, this is evidence that this model cannot be formed accurately, which raises very important questions about the robust features of the closed-loop system to be designed. This problem is quite solvable for various types of marine vessels: the paper [29] is an example. However, the study of the robust features of tracking OD controllers presents an independent problem to be addressed in future studies.

#### *3.2. Forward Speed OD Tracking Controller Design*

In general, the synthesis of the two stabilizing controllers (22) using OD technique can be carried out simultaneously. Nevertheless, one can also apply a combined approach to the feedback design for the considered plant (21) with two control channels. In the range of this approach, the first control is formed based on the OD problem, and the second one can be designed in any other way, providing desirable performance features. However, a joint closed-loop system must have zero equilibrium, and this motion must be asymptotically stable.

Realizing this idea, let us accept the dynamic controller for rudders (marine autopilot) in the following form:

$$\begin{array}{l} \frac{d\mathbf{z}}{dt} = \mathbf{y}\_z(t, \mathbf{z}, V\_{\mathbf{x}}, \mathbf{\xi}),\\ \mu\_2 = \mathbf{g}\_z(t, \mathbf{z}, V\_{\mathbf{x}}, \mathbf{\xi}), \end{array} \tag{26}$$

where **z** ∈ *E*ν is the state vector of this controller.

In particular, the feedback (26) may have a multipurpose structure, which is presented in detail in [22,27] for linear time-invariant case. Its implementation taking into account control time-delay is investigated in the paper [28].

The equations of the control plant now take the form

$$\begin{array}{l}\frac{dV\_{\mathbf{r}}}{dt} = \alpha u\_{1}^{2} + \beta\_{1}(t, V\_{\mathbf{x}}, \xi)u\_{1} + \gamma\_{1}(t, V\_{\mathbf{x}}, \xi),\\\frac{d\xi}{dt} = \mathbf{h}\_{1}(t, V\_{\mathbf{x}}, \xi) + \mathbf{b}\varrho\_{z}(t, \mathbf{z}, V\_{\mathbf{x}}, \xi),\\\frac{d\mathbf{z}}{dt} = \mathbf{y}\_{z}(t, \mathbf{z}, V\_{\mathbf{x}}, \xi).\end{array} \tag{27}$$

In order to synthesize a feedback for the first control channel, i.e., design the OD forward speed stabilizing controller, let us introduce the functional to be damped as

$$L = L(t, V\_{\mathbf{x}}, \xi, \mathbf{z}, u\_1) = V(V\_{\mathbf{x}}, \xi, \mathbf{z}) + \int\_{t\_0}^{t} \lambda\_1^2 u\_1^2 \, d\mathbf{\tau}.\tag{28}$$

Let us take as a Lyapunov function candidate the following sum of quadratic forms:

$$V = V(V\_{\mathbf{x}}, \xi, \mathbf{z}) = \lambda^2 V\_{\mathbf{x}}^2 + \xi^T \mathbf{Q} \xi + \mathbf{z}^T \mathbf{Q}\_1 \mathbf{z}, \ \mathbf{Q} > 0, \ \mathbf{Q}\_1 > 0. \tag{29}$$

Based on (28), (29), we can state the correspondent OD problem of the form

$$\mathcal{W}(t, V\_{\mathbf{x}}, \xi, \mathbf{z}, u\_1) \to \min\_{u\_1 \in \mathcal{U}\_1} \tag{30}$$

where the rate *W* of the functional *L* change is determined as follows:

$$\begin{split} \mathcal{W} &= \mathcal{W}(t, V\_{\mathbf{x}}, \boldsymbol{\xi}, \mathbf{z}, \boldsymbol{u}\_{1}) := \left. \frac{d\mathcal{L}}{dt} \right|\_{\begin{pmatrix} 2\varPi \end{pmatrix}} = \left. \frac{dV}{dt} \right|\_{\begin{pmatrix} 2\varGamma \end{pmatrix}} + \lambda\_{1}^{2} \boldsymbol{u}\_{1}^{2} = \left. \frac{\partial V}{\partial V\_{\mathbf{x}}} \frac{dV\_{\mathbf{x}}}{dt} + \frac{\partial V}{\partial \mathbf{k}} \frac{d\mathbf{k}}{dt} + \frac{\partial V}{\partial \mathbf{z}} \frac{d\mathbf{z}}{dt} + \lambda\_{1}^{2} \boldsymbol{u}\_{1}^{2} \\ &= 2\lambda^{2} V\_{\mathbf{x}} \Big[ \alpha \mathbf{u}\_{1}^{2} + \beta\_{1} (t, V\_{\mathbf{x}}, \boldsymbol{\xi}) \boldsymbol{u}\_{1} + \gamma\_{1} (t, V\_{\mathbf{x}}, \boldsymbol{\xi}) \Big] + 2\xi^{T} \mathbf{Q} \frac{d\mathcal{L}}{dt} + 2\mathbf{z}^{T} \mathbf{Q}\_{1} \frac{d\mathbf{z}}{dt} + \lambda\_{1}^{2} \boldsymbol{u}\_{1}^{2}. \end{split} \tag{31}$$

Let us solve the problem (30), taking into account (31) and assuming that the extreme is achieved at the inner point of the set *U*1: with necessity we obtain

$$\frac{dW}{du\_1} = 2\lambda^2 V\_\mathbf{x} (2\pi u\_1 + \beta\_1 (t, V\_\mathbf{x}, \xi)) + 2\lambda\_1^2 u\_1 = 0$$

that determines the following controller

$$u\_1 = u\_{10}(t, V\_{\mathbf{x}'} \mathfrak{E}) = -\frac{\lambda^2 V\_{\mathbf{x}} \mathfrak{E}\_1(t, V\_{\mathbf{x}'} \mathfrak{E})}{2\lambda^2 V\_{\mathbf{x}} \mathfrak{a} + \lambda\_1^2}.\tag{32}$$

Since, according to formulae (15)–(21), we have

$$\begin{aligned} \beta\_1 &= \beta\_1(t, V\_\ge, \xi) := 2\alpha u\_{1\rho}(t) + \beta \left(V\_\ge + V\_{\ge \rho}(t), \xi + \xi\_\rho(t)\right), \\ \beta(V\_\ge + V\_{\ge \rho}(t), \xi + \xi\_\rho(t)) &= \beta\_0 \sqrt{\left(V\_\ge + V\_{\ge \rho}(t)\right)^2 + \left(V\_\ge + V\_{y\rho}(t)\right)^2}, \ \beta\_0 = \text{const} \end{aligned}$$

we arrive from (32) at the following OD controller:

$$u\_{10}(t, V\_x, \xi) = -\frac{\lambda^2 V\_x \left(2\alpha u\_{1\rho}(t) + \beta\_0 \sqrt{\left(V\_x + V\_{x\rho}(t)\right)^2 + \left(V\_y + V\_{y\rho}(t)\right)^2}\right)}{2\lambda^2 V\_x \alpha + \lambda\_1^2},\tag{33}$$

where *Vy* = 1000 ξ, *Vy*ρ(*t*) = 1000 ξρ(*t*).

It is necessary to note that if the value *u*10(*t*, *Vx*, ξ) is out of an inner part of the admissible set *U*1, we must use the control signal *u*1 = 0.01*P*ρ*u*1ρ(*t*) · sign(*u*10(*t*, *Vx*, ξ)) instead of (33).

It is necessary to note that practical problems involve situations where the right part of the equation for forward speed has a more complex structure than for the system (17). In general, this equation can be presented as

$$\frac{d\nabla\_{\mathbf{x}}}{dt} = F\_{\mathbf{x}}(\overline{V}\_{\mathbf{x}}, \overline{\mathbf{\xi}}\_{\mathbf{\prime}} \overline{u}\_{1}) .$$

which results in the correspondent equation for the system (21) in deflections:

$$\frac{dV\_{\mathbf{x}}}{dt} = G\_{\mathbf{x}}(t\_{\prime}V\_{\mathbf{x}\prime}\boldsymbol{\xi}\_{\prime}\boldsymbol{u}\_{1})\_{\prime\prime}$$

where *Gx*(*t*, *Vx*, ξ, *u*1) := *Fx Vx* + *Vx*ρ(*t*), ξ + ξρ(*t*), *u*1 + *u*1ρ(*t*) − *Fx Vx*ρ(*t*), ξρ(*t*), *u*1ρ(*t*) . In this case, instead of (31) we have

$$\mathcal{W}(t, V\_{\mathbf{x}}, \xi, \mathbf{z}, \boldsymbol{\mu}\_{1}) = 2\lambda^{2} V\_{\mathbf{x}} G\_{\mathbf{x}}(t, V\_{\mathbf{x}}, \xi, \boldsymbol{\mu}\_{1}) + 2\xi^{T} \mathbf{Q} \frac{d\xi}{dt} + 2\mathbf{z}^{T} \mathbf{Q}\_{1} \frac{d\mathbf{z}}{dt} + \lambda\_{1}^{2} \boldsymbol{\mu}\_{1}^{2}.$$

and the analytical search for stationary points of the function *W* becomes problematic. Nevertheless, it is always possible to consider a question about the numerical solution of OD problem (30) for each aggregate of fixed parameters *t*, *Vx*, ξ. Moreover, it is always possible to pose a finite dimensional minimization problem

$$u\_1 = u\_{10}^n(t, V\_{\ge}, \xi) = \arg\min\_{u\_1 \in \mathcal{U}\_{1n}} \left( \mathcal{g}\_1(t, V\_{\ge}, \xi, u\_1) + \lambda\_1^2 u\_1^2 \right),\tag{34}$$

on the finite net *U*1*n* ⊂ *U*1 for any compact set *U*1. This problem should be solved at the time *t* for the fixed parameters *t*, *Vx*(*t*), ξ(*t*). Obviously, with a sufficiently large number of points for the set *U*1*n*, the solution of the problem (34) tends to solution for the same problem on the set *U*1.

Finally, in accordance with formula (23), we can compose the tracking controller of the form

$$
\overline{u}\_1 = \overline{u}\_{01}(t, V\_{\ge \prime} \mathfrak{E}) := u\_{1\rho}(t) + u\_{10}(t, V\_{\ge \prime} \mathfrak{E})\_{\prime},
$$

where the stabilizing part *u*10(*t*, *Vx*, ξ) is determined by equality (33). One can see that the resulting representation for OD forward speed tracking controller is as follows:

$$\overline{u}\_{1} = \overline{u}\_{01}(t, \overline{V}\_{x}, \overline{V}\_{y}) = \frac{\lambda\_{1}^{2} u\_{1\rho}(t) - \lambda^{2} (\overline{V}\_{x} - V\_{x\rho}(t)) \left\| \beta\_{0} \sqrt{\overline{V}\_{x}^{2} + \overline{V}\_{y}^{2}} \right. \tag{35}$$

The current values of the dynamic parameters *Vx*(*t*), *Vy*(*t*) must be measured to implement the controller (35).

#### *3.3. Numerical Example of OD Tracking Controller Synthesis*

To illustrate a practical implementation of the proposed OD approach, let us consider a practical example of forward speed tracking controller design for the transport ship with a displacement of 3500 tons, a length of 110 m, a width of 14 m, and an immersion of 5 m. As a mathematical model of the plant, let us accept equations (17), which are presented in [1,2] with the given functions β(*Vx*, ξ), *γ*(*Vx*, ξ), **h**(*Vx*, ξ) and parameter α.

First, we define the reference motion of the ship as a partial solution of system (17) with the initial conditions *Vx*(0) = 2 *m*/*s*, ξ(0) = 0. Let us accept the reference forward speed control as follows:

$$u\_{1\rho}(t) = \begin{cases} 2 + 0.08t, & \text{if } 0 \le t \le 100s, \\ 10, & \text{if } t \ge 100s. \end{cases} \tag{36}$$

To determine the reference control for the rudders, let us initially construct a model feedback as a controller

$$u\_m = k\_1 \mathbf{x}\_{m1} + k\_2 \mathbf{x}\_{m2} + k\_3 (\mathbf{x}\_{m3} - \boldsymbol{\varphi}\_z) + k\_4 \mathbf{x}\_{m4} \\ := \mathbf{k}\_m \mathbf{x}\_m - k\_3 \boldsymbol{\varphi}\_z. \tag{37}$$

fulfilling the command-heading signal ϕ*z*. The basis item **k***m***x***m* in (37) stabilizes Linear Time Invariant (LTI) plant

$$\frac{d\mathbf{x}\_m}{dt} = \mathbf{A}\_m \mathbf{x}\_m + \mathbf{b}\_m u\_m \tag{38}$$

which is a result of the plant (17) linearization in the neighborhood of the origin for a fixed forward speed *Vx*. Let us especially notice that the linear model (38) is used only for the synthesis of controller (37) but is not used for simulation and for performance indices computation: controller (37) closes full initial model (17) of the vessel.

Accepting *Vx* = 10 m/s , we obtain

$$\mathbf{A}\_{m} = \begin{pmatrix} a\_{11} & a\_{12} & 0 & b\_{1} \\ a\_{21} & a\_{22} & 0 & b\_{2} \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \ \mathbf{b}\_{m} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 1 \end{pmatrix} \\ a\_{11} = -0.0936, \quad a\_{12} = -6.34, \quad b\_{1} = -0.190, \\ b\_{21} = -0.00480, \quad a\_{22} = -0.717, \quad b\_{2} = 0.0160.$$

We design the basis control *um* = **k***m***x***m* for plant (38) with aforementioned parameters as the Linear Quadratic Regulator (LQR) optimal controller with respect to the functional

$$J\_m = J(u\_m) = \int\_0^\infty \left( \mathbf{x}\_m^T \mathbf{Q}\_m \mathbf{x}\_m + \lambda\_m u\_m^2 \right) dt\_\prime \tag{39}$$

where **Q***m* = diag 0 1.975 0.0250 0 , *λm* = 60. Let us especially remark that the choice of the presented parameters for the functional (39) is determined by the initially given requirements to a dynamic quality of the control processes for the nonlinear time varying system. These requirements provide desirable settling time and overshoot for the closed-loop connection.

As a result of computations, we get the constant vector **k***m* = *k*1 *k*2 *k*3 *k*4 , where

$$k\_1 = 0.00234, \ k\_2 = -0.0495, \ k\_3 = -0.0204, \ k\_4 = -0.0497.$$

Let us substitute the reference heading control *u*1ρ(*t*) and the rudders feedback control

$$
\overline{u}\_2 = k\_1 \overline{\xi}\_1 + k\_2 \overline{\xi}\_2 + k\_3 (\overline{\xi}\_3 - \varphi\_z) + k\_4 \overline{\xi}\_4,\\
\varphi\_z = 90^\circ,\tag{40}
$$

into Equation (17). After integration, we obtain the functions *u*2ρ = *u*2ρ(*t*), *Vx*ρ = *Vx*ρ(*t*), and ξρ = ξρ(*t*) = *Vy*ρ(*t*) ωρ(*t*) ϕρ(*t*) δρ(*t*) *T* , which determine given reference motion of the vessel. The graphs of aforementioned functions are presented in Figures 2 and 3.

**Figure 2.** Graphs of the reference functions *u*1ρ(*t*), *u*2ρ(*t*), *Vx*ρ(*t*), and δρ(*t*).

**Figure 3.** Graphs of the reference functions ϕρ(*t*), *Vy*ρ(*t*), and ωρ(*t*).

To provide simulation process, instead of (23), (26), we accept the following simple feedback control law for a marine autopilot:

$$
\overline{u}\_2 = k\_1 V\_y + k\_2 \overline{w} + k\_3 (\overline{\varphi} - \varphi\_z) + k\_4 \overline{\delta}, \tag{41}
$$

where ϕ*z* is the command-heading signal. This controller, corresponding to (40), can be directly used for actuators of the vertical rudders.

Thus, one can see that, as a result, initial control plant (17) is closed by the tracking controllers (35) and (41). The current values of the dynamic parameters *Vx*(*t*), *Vy*(*t*), ω(*t*), ϕ(*t*), and δ(*t*) must be measured by the corresponding sensors for the actual implementation of these controllers.

For simulation of the closed-loop system dynamics, the following parameters values are accepted: *λ*2 = 100, *λ*2 1 = 1, α = 0.00462, β0 = −0.00322. In addition, let us take into account the restrictions δ(*t*) ≤ *dm* = 35◦ and |*u*2(*t*)| ≤ *um* = 3 ◦/c.

To illustrate the practical applicability of the proposed approach, let us simulate the control processes for the closed-loop connection. The aim is to make the proposed approach comparable to other methods. This determined the choice of design parameters and regimes of vessel's motion. These regimes represent the most popular options for movement on quiet water and under the action of sea waves.

The results of simulation are presented in Figures 4–6 as the graphs of corresponding functions, which reflect control signals and the vessel's state variables for the transient process. This process is determined by the aforementioned reference motion, which is realized with the help of the designed tracking controllers. Initial conditions for all variables are zero with the exception of forward speed and heading angle. By these variables, the initial conditions *Vx*(0) = 4 m/s and ϕ(0) = −10◦ are accepted to distinguish them from the reference motion.

**Figure 4.** First control signals *u*1 = *u*1*b*(*t*), *u*1ρ(*t*), and *u*1(*t*) = *u*1(*t*) − *u*1ρ(*t*).

**Figure 5.** Forward speed *Vx* = *Vxb*(*t*), *Vx*ρ(*t*), and *Vx*(*t*) = *Vx*(*t*) − *Vx*ρ(*t*).

**Figure 6.** Rudders deflections δ = δ*b*(*t*), δρ(*t*), and δ(*t*) = δ(*t*) − δρ(*t*).

Let us note that the dynamical quality of the presented transient process seems to be quite satisfactory. In addition, it is suitable to illustrate the dynamics of the closed-loop connection under the action of sea wave external disturbance. Figure 7 shows the same control process for the forward speed *Vx*, taking into account presence the approximate representation of sea waves with an intensity of 5 on the Beaufort scale.

**Figure 7.** Forward speed *Vx* = *Vxb*(*t*), *Vx*ρ(*t*), and *Vx*(*t*) = *Vx*(*t*) − *Vx*ρ(*t*) under sea waves action.

Finally, let us analyze the stability of the closed-loop system (17), (35), (41) with synthesized tracking controllers. One can easily see that in order to provide aforementioned analysis, it is sufficient to consider the zero-equilibrium stability for the closed-loop system (27), (33) presented in deflections from the reference motion.

It is obvious that systems (27), (33) have zero equilibrium, and to investigate its stability features, let us introduces the following Lyapunov function candidate:

$$V = V(\mathbf{x}) = \lambda^2 V\_x^2 + \xi^T \mathbf{S} \boldsymbol{\xi} = \mathbf{x}^T \mathbf{S}\_1 \mathbf{x}, \; \mathbf{S}\_1 = \left(\cdot, \frac{\lambda^2}{0}, \frac{\lambda}{\cdot} \frac{0}{\mathbf{S}} \cdot \right), \; \mathbf{x} := \begin{pmatrix} V\_x \\ \xi \end{pmatrix}. \tag{42}$$

Here, the symmetrical matrix **S** > 0 is a solution of the algebraic Riccati equation, which is used in the range of LQR controller (37) synthesis:

$$\mathbf{S} = \begin{pmatrix} 0.0107 & -0.243 & -0.102 & -0.141 \\ -0.243 & 7.04 & 2.59 & 2.97 \\ -0.102 & 2.59 & 1.27 & 1.22 \\ -0.141 & 2.97 & 1.22 & 2.98 \end{pmatrix}.$$

Let us notice that the introduced function *V*(**x**) satisfy the relationships

$$
\alpha\_1(\|\mathbf{x}\|) \le V(\mathbf{x}) \le \alpha\_2(\|\mathbf{x}\|), \forall \mathbf{x} \in E^5,\tag{43}
$$

with the following *K*-class functions:

$$\alpha\_1(\left\|\mathbf{x}\right\|) = \sigma\_{\min}(\mathbf{S}\_1) \cdot \left\|\mathbf{x}\right\|^2, \ \alpha\_2(\left\|\mathbf{x}\right\|) = \sigma\_{\max}(\mathbf{S}\_1) \cdot \left\|\mathbf{x}\right\|^2.$$

where *σ*min = 0.00115 and *σ*max = 100 are the minimum and maximum eigenvalues of the matrix **S**1, respectively.

Using the additional function α3(**x**) ≡ **x** ∈ *K*, one can check if the following inequality is valid

$$\mathcal{W}(t, \mathbf{x}) = \frac{dV}{dt}\Big|\_{(2\mathcal{T}),(3\mathbf{3})} \le -\alpha\_{\mathbb{B}}(\|\mathbf{x}\|), \forall t \ge 0, \ \forall \mathbf{x} \in B\_{\mathbf{x}}.\tag{44}$$

where *Bx* is a box, determined by the relationships:

$$B\_{\mathbf{x}}(t) = \left\{ \mathbf{x} \in E^5 : |\mathbf{x}\_i| \le \mathbf{x}\_{\text{im}}(t), \ i = \overline{1,5} \right\},$$

*x*1*m*(*t*) = 0.01*P*ρ|*Vx*ρ(*t*)|, *x*2*m*(*t*) = 0.01*P*ρ *Vy*ρ(*t*) ,*x*3*m*(*t*) = 0.01*P*ρ|ωρ(*t*)|, *x*4*m*(*t*) = 0.01*P*ρ|ϕρ(*t*)|,

$$\mathbf{x}\_{\mathsf{Sm}}(t) = 0.01 P\_{\mathsf{P}} |\boldsymbol{\delta}\_{\boldsymbol{\mathcal{P}}}(t)|\ldots$$

Here, the variable *P*ρ determines the relative width of the box *Bx* compared to the current values of the reference signals. This variable was increased to such a value that the condition (44) was met. The obtained value *P*ρ = 25% seems to be still admissible in the range of (43), (44) [4,5], determining the region of the local uniform asymptotic stability for the reference motion, which is realized by tracking controllers (35), (41).

#### **4. Discussion**

The main goal of this work was to propose constructive methods for marine tracking controllers' design taking into account the real conditions of a vessel's motions. We focused our main attention on a situation where the rudders' deflections and the forward speed are presented by initially given reference signals to be realized using tracking feedback controls.

This problem can be solved using different popular optimization approaches (Bellman's theory, MPC technique, sliding mode control, etc.). Nevertheless, we propose a new specific method for tracking controllers' design, which ensures the desirable reference motion of the vessel along the forward speed and heading angle.

This method is based on the optimal damping concept, which has certain advantages related to the practical requirements for the dynamic features of a closed-loop connection. The main advantage of the aforementioned approach is that the numerical solution of the optimization problems is essentially simplified. In contrast to well-known approaches [1–4], we applied OD tracking feedback as a control law with special features that allow it to be adjusted and implemented in real-time regime of motion.

The main result of this study is the development of the optimal damping concept to ensure its practical applicability and effectiveness that is illustrated by a controller design for a transport ship.

The investigations presented above could be further developed to consider the robust features of the tracking control laws, information about the measurement noise, and the presence of transport delays [24]. The results of the executed research could also be implemented to provide desirable reference motion of the dynamic positioning systems for marine vessels [26,27]. Certain attention may be given to the multipurpose control laws applications [22–25]. As for direct development of this study, it is possible to also use OD controller for the vertical rudders. The extension of the proposed approach to various robotic systems is also of considerable interest. The scope of the proposed approach may additionally include remotely operated vehicles [11] and offshore structures [12].

**Funding:** This research was funded by the Russian Foundation for Basic Research (RFBR) controlled by the Government of Russian Federation, research project number 20-07-00531.

**Acknowledgments:** The author expresses his sincere gratitude to A.P. Zhabko for the interest to the topic of this study and for the high-level professional suggestions.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**