### 1. Introduction

### 2. Preliminary

### 2.1 Ship Dynamic Model

*ψ*as follows: where

*η*=[

*x*,

*y*,

*ψ*]

^{⊤}and

*ν*=[

*u*,

*v*,

*r*]

^{⊤}are the position and velocity vectors, respectively, and

*J*(

*ψ*) is the rotation matrix that transforms the velocities from the body-fixed coordinate to the earth-fixed coordinate.

*M*represents the inertia matrix,

*C*(

*ν*) stands for the Coriolis and centripetal matrix,

*D*(

*ν*) signifies the damping matrix, and

*τ*is the vector of control force. The matrices

*M*and

*C*(

*ν*) can be depicted as the sum of the rigid body and added mass components, as outlined below:

### 2.2 Actuator Model

*T*

_{1}represents the force produced by the port stern propeller,

*T*

_{2}corresponds to the force originating from the stern starboard propeller, while

*T*

_{3}designates the bow tunnel propeller forces. Additionally, the rudder angles at the port and starboard sides are denoted by

*δ*

_{1}and

*δ*

_{2}, respectively. Given these, the control force can be formulated with the control input [

*T*

_{1},

*T*

_{2},

*δ*

_{1},

*δ*

_{2},

*T*

_{3}]

^{⊤}, as detailed below:

*D*

_{1}and

*D*

_{2}represent the drag forces, while

*L*

_{1}and

*L*

_{2}represent the lift forces for the port and starboard rudder, respectively. The terms

*l*

_{x1},

*l*

_{x2}and

*l*refer to the distances from the center of gravity to each propeller. The force generated by the propeller can be defined by the rotational speed

_{y}*n*, and the thruster coefficient

_{i}*K*with advance ratio

_{T}*J*, as follows: where

_{i}*ρ*denotes the density of water and

*D*stands for the propeller's diameter. The tunnel thruster is also modeled using the rotational speed

*n*

_{3}, the propeller's diameter

*D*

_{3}, and the thruster coefficient

*K*

_{T3}. However, it is important to note that the tunnel thruster can only generate sufficient force when the ship moves at low speeds. Consequently, we employ the tunnel thruster at lower speeds and model it with a constant thruster coefficient, independent of the advance ratio, as follows:

### 3. Berthing Trajectory Optimization

### 3.1 Two-phase Approach

### 3.2 Two-point Boundary Value Problem

*x*=[

*x*,

*y*,

*ψ*,

*u*,

*v*,

*r*]

^{⊤}, as detailed below:

*T*

_{f}_{1}and (

*T*

_{f}_{2}−

*T*

_{f}_{1}) represent the time durations for each phase, respectively. Additionally,

*J*denotes the time penalty cost function, which is defined as follows:

_{T}*ρ*is the weighting scalar.

_{t}*J*(

_{a}*t*),

*Q*,

_{a}*P*, and

_{a}*R*denote the weight matrices that penalize the velocity, control input, and rate of change of control input, respectively.

_{a}*η*is the final position vector. Additionally,

_{f}*J*(

_{t}*t*) denotes the cost function for the terminal phase, given as follows:

*Q*,

_{t}*P*, and

_{t}*R*represent the weight matrices that penalize the velocity, control input, and rate of change of control input, respectively. The optimization problem formulated above also encompasses several constraints, which are as follows:

_{t}*f*and

_{a}*f*represent the dynamic models of the approaching and terminal phases, respectively, while

_{t}*x*

_{0}and

*x*denote the initial and target berthing states.

_{f}*and 𝕊*

_{b}*indicate the safety boundary of the ship and state constraints, as shown in Fig. 4.*

_{s}