Control Co-Design for Buoyancy-Controlled MHK Turbine: A Nested Optimization of Geometry and Spatial-Temporal Path Planning

This work focuses on optimizing an OCT system that was initially designed with a rated power of 700 $\mathrm{kW}$ for harnessing power from the Gulf Stream off Florida’s East Coast following an OCT prototype from IHI Corp. [22], shown in Fig. 1. The turbine was designed such that it will operate at a depth of 50 $\mathrm{m}$ in a homogeneous current speed of 1.6 $\mathrm{m/s}$ assuming 50% filled buoyancy tanks. A schematic of the investigated OCT is shown in Fig. 2, consisting of a body tethered to an anchor via a 607 $\mathrm{m}$ mooring cable and a variable buoyancy tank to control its vertical motion [8].

To formulate the OCT’s kinematics, this system is described through five coordinate frames, including the body-fixed $(\mathcal{T}_{B})$, the inertial $(\mathcal{T}_{I})$, the momentum mesh $(\mathcal{T}_{M})$, the shaft $(\mathcal{T}_{S})$, and multiple rotor blade $(\mathcal{T}_{R})$ frames.

Equations of Motion: Seven degree-of-freedom (DOF) equations of motion are created to describe the dynamic model of the OCT system as suggested in [24], with all forces and moments are shifted to the $\mathcal{T}_{\mathrm{B}}$ frame. The system is represented by 14 states, including the position of OCT body $\mathcal{P}=[x~{}y~{}z]$, linear velocity of the OCT body $\mathcal{\dot{P}}=[u~{}v~{}w]$, Euler angles of the OCT body $\Theta=[\phi~{}\theta~{}\psi]$, angular velocity of the OCT body $\dot{\Theta}=[p_{\mathrm{b}}~{}q~{}r]$, angular velocity of the rotor $p_{\mathrm{r}}$, and rotation angle of the rotor blade $\phi_{\mathrm{r}}$. These equations couple a 6-DOF motion of the OCT body,

$\displaystyle\begin{bmatrix}\mathcal{\ddot{P}}\\ \ddot{\Theta}\end{bmatrix}=\mathcal{M}^{-1}\mathcal{F}$

(1)

where

$\mathcal{M}=\begin{bmatrix}m&0&0&0&m_{\mathrm{b}}z_{\mathrm{cg_{b}}}&0\\[8.61108pt] 0&m&0&-m_{\mathrm{b}}z_{\mathrm{cg_{b}}}&0&mx_{\mathrm{cg}}\\[8.61108pt] 0&0&m&0&-mx_{\mathrm{cg}}&0\\[8.61108pt] 0&-m_{\mathrm{b}}z_{\mathrm{cg_{b}}}&0&I_{x_{\mathrm{{b}}}}&0&-I_{xz_{\mathrm{b}}}\\[8.61108pt] m_{\mathrm{{b}}}z_{\mathrm{cg_{b}}}&0&-mx_{\mathrm{cg}}&0&I_{y}&0\\[8.61108pt] 0&mx_{\mathrm{{cg}}}&0&-I_{xz_{\mathrm{{b}}}}&0&I_{z}\end{bmatrix}$

(2)

$\mathcal{F}=\begin{bmatrix}f_{x}+m(vr-wq)+mx_{\mathrm{cg}}\left(q^{2}+r^{2}\right)-m_{\mathrm{b}}z_{cg_{b}}p_{b}r\\[8.61108pt] f_{y}-mur+w\left(m_{\mathrm{b}}p_{\mathrm{b}}+m_{\mathrm{r}}p_{\mathrm{r}}\right)-m_{\mathrm{b}}z_{\mathrm{cg_{b}}}qr\\ -m_{\mathrm{b}}x_{\mathrm{cg_{b}}}qp_{\mathrm{b}}-m_{\mathrm{r}}x_{\mathrm{cg_{r}}}qp_{\mathrm{r}}\\[8.61108pt] f_{z}+muq-v\left(m_{\mathrm{b}}p_{\mathrm{b}}+m_{\mathrm{r}}p_{\mathrm{r}}\right)+m_{\mathrm{b}}z_{\mathrm{cg_{b}}}\left(p_{\mathrm{b}}^{2}+q^{2}\right)\\ -m_{\mathrm{b}}x_{\mathrm{cg_{b}}}rp_{\mathrm{b}}-m_{\mathrm{r}}x_{\mathrm{cg_{r}}}rp_{\mathrm{r}}\\[8.61108pt] M_{x_{\mathrm{b}}}+\tau_{\mathrm{em}}-qr\left(I_{z_{\mathrm{b}}}-I_{y_{\mathrm{b}}}\right)+I_{xz_{\mathrm{b}}}p_{\mathrm{b}}q\\ -m_{\mathrm{b}}z_{\mathrm{cg_{b}}}\left(wp_{\mathrm{b}}-ur\right)\\[8.61108pt] M_{y}-rp_{\mathrm{b}}\left(I_{x_{\mathrm{b}}}-I_{z_{\mathrm{b}}}\right)-rp_{\mathrm{r}}\left(I_{x_{\mathrm{r}}}-I_{z_{\mathrm{r}}}\right)-I_{xz_{\mathrm{b}}}\left(p_{\mathrm{b}}^{2}-r^{2}\right)\\ +m_{\mathrm{b}}z_{\mathrm{cg_{b}}}(vr-wq)-mx_{\mathrm{cg}}uq+m_{\mathrm{b}}x_{\mathrm{cg_{b}}}vp_{\mathrm{b}}+m_{\mathrm{r}}x_{\mathrm{cg_{r}}}vp_{\mathrm{r}}\\[8.61108pt] M_{z}-qp_{\mathrm{b}}\left(I_{y_{\mathrm{b}}}-I_{x_{\mathrm{b}}}\right)-qp_{\mathrm{r}}\left(I_{y_{\mathrm{r}}}-I_{x_{\mathrm{r}}}\right)-I_{xz_{\mathrm{b}}}rq\\ -mx_{\mathrm{cg}}\text{ u }r+m_{\mathrm{b}}x_{\mathrm{cg_{b}}}wp_{\mathrm{b}}+m_{\mathrm{r}}^{v}x_{\mathrm{cg_{r}}}wp_{\mathrm{r}}\end{bmatrix}$

(3)

and a 1-DOF rotation of the rotor about the x-axis of the body frame:

$$\dot{p}_{\mathrm{r}}=\frac{M_{x_{\mathrm{r}}}-\tau_{\mathrm{em}}-qr(I_{z_{\mathrm{r}}}-I_{y_{\mathrm{r}}})}{I_{x_{\mathrm{r}}}}$$

(4)

where, $f_{(.)}$ denotes the force, and $M_{(.)}$ is the moment. Note that the mass $m_{(.)}$, the moment of inertia $I_{(.)}$, and the center of gravity ${(.)}_{\mathrm{cg}}$ include both the actual inertial properties and added inertial properties of the OCT (denoted as virtual) as suggested in [24]. $(.)_{x}$, $(.)_{y}$, and $(.)_{z}$ denote the portion $(.)$ about $x-$, $y-$, and $z-$ axes. $\tau_{\mathrm{em}}$ denotes the electromechanical shaft torque. Eventually, $(.)_{\mathrm{r}}$ and $(.)_{\mathrm{b}}$ refer to the rotor and the body.

The total forcing on the OCT system is modeled through the gravitational and buoyancy force $f_{\mathrm{gb}}$, rotor force $f_{\mathrm{r}}$, body force $f_{\mathrm{b}}$, and cable force $f_{\mathrm{c}}$, namely $f=f_{\mathrm{gb}}+f_{\mathrm{r}}+f_{\mathrm{b}}+f_{\mathrm{c}}$; the total moments are accordingly calculated as fully presented in [24]. To avoid complexities arise when applying the nonlinear dynamic model, the OCT model is linearized for controller development, as detailed in the followings.

Linear Model: To form the linear state space model for the OCT system, the equations of motion are averaged around the nominal condition (equilibrium point), thereby formulating the linear model as follows:

$\displaystyle\delta\dot{x}_{n\times 1}$

$\displaystyle=A_{n\times n}\delta x_{n\times 1}+B_{n\times m}\delta u_{m\times 1}$

(5)

Here, the states $\delta x$ and control inputs $\delta u$ are defined as deviation from the nominal condition. $\delta x\in\mathbb{R}^{n}$ and $\delta u\in\mathbb{R}^{m}$ with $n$ and $m$ denoting the number of states and control inputs used by the linear model ($n=13$ and $m=3$ in our OCT system), namely:

$$\delta x=[\delta u~{}\delta v~{}\delta w~{}\delta p_{b}~{}\delta p_{r}~{}\delta q~{}\delta r~{}\delta x~{}\delta y~{}\delta z~{}\delta\phi~{}\delta\theta~{}\delta\psi]$$

(6)

$$\delta u=[\delta B_{\mathrm{f}}~{}\delta B_{\mathrm{a}}~{}\delta\tau_{\mathrm{em}}]$$

(7)

where the control inputs are defined by forward buoyancy tank fill fraction $B_{\mathrm{f}}$, aft buoyancy tank fill fraction $B_{\mathrm{a}}$, and $\tau_{\mathrm{em}}$. The fill fractions are limited by the ratio between the buoyancy tank size $\nu_{\mathrm{B}}$ and the base buoyancy tank size $\nu_{\mathrm{B}}^{\mathrm{b}}$ (as formulated in (21)). It is noted that the linear model has one less state than the non-linear model as the rotor rotation angle state, $\phi_{\mathrm{r}}$, is eliminated during the linearization process as described in [25]. Given that the nominal condition is defined by an averaged homogeneous flow speed of $v_{\mathrm{eq}}=1.6\mathrm{m/s}$, the nominal states and control inputs are characterized by $x_{\mathrm{eq}}=[0~{}0~{}0~{}0~{}1.49~{}0~{}0~{}554.50~{}0.38~{}50~{}0.01~{}0.00~{}3.14]$ and $u_{\mathrm{eq}}=[0.4677~{}0.4677~{}-188280]$.

Note that among all the states the major focus here is on the vertical position (depth) $z$, e.g., vertical path planning, which is primarily controlled through the fill fractions $B_{(.)}$. To further reduce the complexity of the model, this work leverages a linear equation relating $z$ to $B_{(.)}$ based on the authors’ previous work [26]. To control the operating depth, the OCT system should be able to change and hold the optimal depth. Hence, the fill fractions must be set such that the OCT system is placed in an optimal depth and maintains that depth when the flow velocities change. The fill fraction at time instant $k$ is calculated by:

$$B_{\mathrm{(.)}}(k+1)=B_{\mathrm{(.)}}(k)+\Delta B_{\mathrm{(.)}}^{\mathrm{HD}}+\Delta B_{\mathrm{(.)}}^{\mathrm{CD}}$$

(8)

where $\Delta B_{\mathrm{(.)}}^{\mathrm{HD}}$ and $\Delta B_{\mathrm{(.)}}^{\mathrm{CD}}$ are the changes in the fill fractions to hold the operating depth and change the depth, respectively. These fill fraction changes are calculated assuming the linear equations between fill fraction changes and depth changes $\frac{dF_{\text{F}}}{dz}$ as well as flow speed changes and depth changes $\frac{dv}{dz}$, namely:

$$\Delta B_{\mathrm{(.)}}^{\mathrm{HD}}=\frac{dF_{\text{F}}}{dz}\frac{dz}{dv}\Delta v=\kappa_{1}\Delta v$$

(9)

$$\Delta B_{\mathrm{(.)}}^{\mathrm{CD}}=\frac{dF_{\text{F}}}{dz}\Delta z=\kappa_{2}\Delta z$$

(10)

where $\kappa_{1}=\frac{dF_{\text{F}}}{dz}\frac{dz}{dv}$ and $\kappa_{2}=\frac{dF_{\text{F}}}{dz}$ are the constant coefficients, which are presented in Table I. Each fill fraction, $B_{(.)}$, is limited between 0 (chamber filled with water) and 1 (chamber filled with air).

OCT Output Power Model: The net harnessed power from the OCT system has been previously formulated in [26], consisting of three terms of (i) power generation $P_{\mathrm{OCT}}$; (ii) power consumption to hold the operating depth $P_{\mathrm{HD}}$; and (iii) power consumption to change the operating depth and relocate the OCT $P_{\mathrm{CD}}$, namely:

$$P_{\mathrm{net}}=P_{\mathrm{OCT}}-P_{\mathrm{HD}}-P_{\mathrm{CD}}$$

(11)

where

$$P_{\mathrm{OCT}}=\min(\frac{1}{2}\rho\pi(\frac{d}{2})^{2}v^{3}c_{\mathrm{p}},P^{\mathrm{r}}_{\mathrm{g}})$$

(12)

$$P_{\mathrm{HD}}=\begin{cases}0,&\text{$\Delta v<0$}\\ \frac{\zeta(\kappa_{1}\Delta v)}{\Delta t},&\text{$\Delta v>0$}\\ \end{cases}$$

(13)

$$P_{\mathrm{CD}}=\begin{cases}0,&\text{$\Delta z>0$}\\ \frac{\zeta(\kappa_{2}\Delta z)}{\Delta t},&\text{$\Delta z<0$}\\ \end{cases}$$

(14)

where $\rho$ is the water density, $d$ is the rotor diameter, $v$ is the ocean velocity, $c_{\mathrm{p}}$ denotes the average power coefficient, $P^{\mathrm{r}}_{\mathrm{g}}$ is the rated power of the OCT, $\Delta t$ is the sampling time, and $\zeta$ denotes the constant coefficient.

A prominent design of the OCT system may entail several prerequisites: (i) harnessing the maximum power while tethered to the anchor; (ii) conducting a real-time search to find the vertical position with the maximum power subject to consuming the minimal power to relocate the OCT system; and (iii) acquiring the maximum power with a minimum OCT weight. To address these conditions, the dominant design parameters of the OCT system should be pointed out:

Rotor: The rotor is a key parameter in the OCT system design, directly affecting the whole system design and power generation. To highlight the importance of the rotor, a complete study on the hydrodynamic design of the OCT’s rotor (airfoil size and shape with the corresponding power and thrust coefficients) has been conducted in [27]. On account of the results obtained through [27], the main focus in the current study is on rotor diameter. Larger rotor diameters increase the power production from OCT systems operating in a homogeneous flow field in the order of the square of the rotor diameter (see (12)), but with increased hydrodynamic system drag/thrust (also a function of the square of the diameter) and increased weight. The increased drag/thrust force on the rotor requires increased buoyancy forces to maintain the same elevation, which can be achieved by maintaining less water in the buoyancy tanks. However, for higher flow speeds, this increases the minimum achievable operating depth (i.e., depth when buoyancy tanks are completely filled with air [8]), resulting in a decrease in the available energy density as available energy is typically strongest near the sea surface.

Buoyancy Tank: In our investigated OCT system, the buoyancy tank is another important design parameter that keeps the system afloat. Accordingly, the operating depth of the tethered OCT plant is initially determined through the amount of water in the tank, thereby imposing constraints on the vertical forcing and movement based on tank size. The results obtained in the authors’ previous work [8] suggest that the buoyancy tank size was not optimally designed since the OCT system can operate in a wide range of depth (e.g., $z=27$ $\mathrm{m}$ to 216 $\mathrm{m}$ for a flow speed of 2.5 $\mathrm{m/s}$) and can hit the depth of 50 $\mathrm{m}$ even with a strong flow speed of 2.5 $\mathrm{m/s}$. Note that in the buoyancy tank, a pump drives water while setting the tank’s pressure at vacuum pressure ($\cong 0$ $\mathrm{kPa}$). Therefore, the key design parameter is the tank size (denoted as $\nu_{\mathrm{B}}$ for the volume of each of the two buoyancy tanks). A larger tank size allows for the larger movement of the OCT, but increases the system’s weight and cost.

Generator: To design a reliable and efficient OCT plant, the generator configuration plays an important role, limiting the power produced from the ocean currents (again, see (12)). In the current control co-design framework, the primary concentration is on the generator rated power. Picking a generator with a large rated power allows an electrical power output that is directly related to available hydrodynamic power from the ocean currents, however, at the expense of increased generator size, weight, and cost. Therefore, optimizing generator power rating to maximize the power-to-weight ratio (a surrogate for the cost of electricity) requires site specific energy resource data (commonly measured through the “power density” [28]).

The efficacy of the proposed control co-design framework is tested on a sample buoyancy-controlled OCT detailed in Section II-B. The base design parameters suggested in [8] are presented in Table I. The simulations are fully conducted in Python on a machine equipped with a 2.3 $\mathrm{\,GHz}$ CPU and 32 $\mathrm{\,GB}$ of RAM. For the MPC-based spatial-temporal power optimization algorithm, the time step is set as 1 hour, the prediction horizon is set as $T=2$ hours, and the number of discrete depths is 17 within the allowable depth range (i.e., between $z_{\mathrm{min}}$ and $z_{\mathrm{max}}$ as presented in Table I). The design parameters are assumed to change within $\upsilon^{\min}=0.1\times\upsilon^{b}$ and $\upsilon^{\max}=1.1\times\upsilon^{b}$ to help ensure that the compatibility is maintained with the OCT components that are satisfied at their baseline values. Also, to introduce spatial-temporal uncertainties of the ocean currents into our model, we use a set of real measured ocean data at a latitude of $26.09^{\circ}N$ and longitude of $-79.80^{\circ}E$ by a 75 kHz acoustic Doppler current profiler (ADCP) with a vertical resolution of 5 $\mathrm{m}$ within 400 $\mathrm{m}$ [28].

The following three cases are simulated and compared:

• •

  Baseline OCT design: The baseline OCT design investigated in this study is discussed in Section II, and its major design parameters are presented in Table I. The mechanical and structural design of this system was detailed in [8], which lacks an analysis to account for the structural design and its effect on the harnessed power. Here, the spatial-temporal optimization is performed for the baseline OCT system, and the maximum harnessed power along with the power-to-weight ratio for the baseline are found.

• •

  Co-design for a single OCT design parameter: The co-design framework solves the optimization problem of maximizing the power-to-weight ratio for a single OCT design parameter. Three separate optimizations are performed for the buoyancy tank, generator, and rotor as major design parameters.

• •

  Co-design for multiple OCT design parameters: The proposed co-design framework (as shown in Fig. 3) evaluates multiple design parameters at the same time in relation to the maximum power-to-weight for the OCT system when operating in the spatial-temporal uncertain ocean environment.

Fig. 4 compares the obtained power-to-weight ratios for changes applied in each design parameter ($\pm 10$ % from its base value). From this figure, the rotor diameter is the dominant design parameter, which significantly affects the power-to-weight ratio. The design parameters of the baseline OCT, the optimal design parameters obtained through co-design for the rotor, and co-design for multiple parameters simultaneously are presented in Table II. The hourly power-to-weight ratios for baseline design and optimal designs for the other two cases, along with the corresponding ocean current velocity (i.e., the turbine operating environment), are shown in Fig. 5. The optimal OCT designs for case II and case III generate the same power (e.g., $P(\upsilon)=256.99$ $\mathrm{kW}$) but different total mass. The rotor diameter dominates other design parameters due to its significant mass and effect on the power generation ($P_{\mathrm{OCT}}$). From all these three cases, we see that the generated powers $P(\upsilon)$ are smaller than the baseline design of 700 $\mathrm{kW}$, which means that the baseline design was away from the optimal. The reason is that the average ocean current speed that we used in this paper is around 0.67 $\mathrm{m/s}$, much smaller than the speed considered for formulating the baseline design 1.015 $\mathrm{m/s}$. Through control co-design, the optimal generator rated power in case III is largely decreased from 700 $\mathrm{kW}$ in baseline design to 495 $\mathrm{kW}$, much closer to the harvesting power 256.99 $\mathrm{kW}$. These results highlight the importance of applying the geometric/spatial-temporal control co-design approach to design the optimal OCT in dynamic ocean environments.

This work is the first attempt to control co-design for a buoyancy-controlled OCT system. The efficacy of the proposed approach, which maximizes the power-to-weight ratio and finds the optimal values of three key design parameters, was validated through simulations. The current study can be expanded on different aspects. The cost minimization objective should be added to the control co-design framework in addition to the optimization of the power-to-weight ratio. Furthermore, the spatial-temporal optimization in the current format ignores the complete linear model of the OCT and, in a broader view, a dynamic model of the OCT as well as the coupling between all OCT states. Although the presented model takes into account measured ocean currents, the increased forcing on different components is neglected, which is bounded by defining an upper limit for changing the design parameters. Finally, the ultimate goal of control co-design is to optimize the whole OCT design parameters with detailed modeling of the corresponding couplings.