# Nested Control Co-design of a Spar Buoy Horizontal-axis Floating Offshore
Wind Turbine

Saeid Bayat  Yong Hoon Lee  James T. Allison

###### Abstract

Floating offshore wind turbine (FOWT) systems involve several coupled physical
analysis disciplines, including aeroelasticity, multi-body structural
dynamics, hydrodynamics, and controls. Conventionally, physical structure
(plant) and control design decisions are treated as two separate problems, and
generally, control design is performed after the plant design is complete.
However, this sequential design approach cannot fully capitalize upon the
synergy between plant and control design decisions. These conventional design
practices produce suboptimal designs, especially in cases with strong coupling
between plant and control design decisions. Control co-design (CCD) is a
holistic design approach that accounts fully for plant-control design coupling
by optimizing these decisions simultaneously. CCD is especially advantageous
for system design problems with complex interactions between physics
disciplines, which is the case for FOWT systems. This paper presents and
demonstrates a nested CCD approach using open-loop optimal control (OLOC) for
a simplified reduced-order model that simulates FOWT dynamic behavior. This
simplified model is helpful for optimization studies due to its computational
efficiency, but is still sufficiently rich enough to capture important
multidisciplinary physics couplings and plant-control design coupling
associated with a horizontal-axis FOWT system with a spar buoy floating
platform. The CCD result shows an improvement in the objective function,
annual energy production (AEP), compared to the baseline design by more than
eleven percent. Optimization studies at this fidelity level can provide system
design engineers with insights into design directions that leverage design
coupling to improve performance. These studies also provide a template for
future more detailed turbine CCD optimization studies that utilize higher
fidelity models and design representations.

###### keywords:

floating offshore wind turbine (FOWT) , spar buoy platform , control co-design
(CCD) , open-loop optimal control (OLOC) , reduced-order modeling (ROM)

[author-ise] department = Industrial and Enterprise Systems Engineering,,
organization = University of Illinois at Urbana-Champaign, city=Urbana,
state=IL, country=USA [author-uofm] department = Mechanical Engineering,,
organization = The University of Memphis, city=Memphis, state=TN, country=USA

##  1 Introduction

Over the last decade, offshore wind energy has gained significant attention as
an environmentally friendly and cost-effective energy source. Offshore wind
resources are stronger and more consistent than their onshore counterparts,
while exhibiting fewer issues, such as terrain-induced shear effects, acoustic
noise, visual impact, size limits, and interference with humans and wildlife.
However, many ideal offshore wind locations have water depths beyond 45
meters, rendering monopile foundations installed directly on the sea bed
unsuitable [1]. Therefore, there is a growing need for technical advancements
in floating offshore wind turbine (FOWT) systems to harness the profound
energy potential in deep water regions.

FOWTs are still in their early stages of development and have yet to be
commercially deployed. To facilitate large-scale commercialization of the
FOWTs, several challenges must be addressed, which can be mitigated through
innovative design discoveries and optimization strategies that leverage
inherent complex design couplings and significantly reduce energy costs. For
example, advanced control strategies can help enhance energy production while
simultaneously protecting the turbine system from motions induced by
hydrodynamic interactions with the floating platform. Furthermore, FOWT
systems consist of numerous components with multidisciplinary interactions,
including hydrostatics, hydrodynamics, mooring dynamics, structural dynamics,
aerodynamics, and controls, complicating the design process. Thus, advancing
FOWT performance requires advanced multidisciplinary design methodologies.

This article presents a study of the integrated physical (plant) and control
design of a spar buoy-based horizontal-axis FOWT system using a model
comprised of a set of simplified governing physics formulations. The FOWT
model development is intended for use with integrated physical and control
design optimization (often referred to as a control co-design, CCD), which
utilizes coupling between plant and control design variables to better exploit
the properties of the combined plant-control design space. The design problem
is simple enough to directly demonstrate how FOWT CCD can be formulated but
contains enough richness to explore important physics and design couplings.

Direct optimal control methods have proven particularly effective in the
flexible exploration of new control strategies while accounting for the
realities of physical system design, such as design coupling and failure modes
[2, 3, 4, 5, 6, 7]. Direct transcription is a class of direct optimal control
methods that overcomes many of the shortcomings of conventional shooting-based
direct optimal control methods [8] but requires that optimization solvers have
direct access to the time derivative function. Instead of treating system
dynamics as a black box, direct transcription capitalizes on the problem
structure inherent to continuous dynamic system optimization problems to
enhance numerical robustness and efficiency. Here, in this study, we use
pseudospectral-based direct transcription method [9, 10, 11] for obtaining
optimal control solutions.

Engineered systems are often operated using active control systems, and
designing these systems requires engineers to make both physical system
(plant) and control system design decisions. Plant-control design coupling
exists if changing plant design influences what control decisions are best for
overall system performance, and if control design changes influence what plant
design decisions are best. Despite the coupling mentioned above, a sequential
procedure has largely been used to design these systems in practice (control
design performed after completing plant design). Sequential approaches,
however, leave potential performance gains on the table by not fully
exploiting design coupling to produce system-optimal designs [12].

Here, an integrated design approach is used that fully accounts for bi-
directional plant-control design coupling and produces system-optimal designs,
generally referred to as contro co-design (CCD). The CCD optimization is a
class of integrated design methods that account explicitly for plant-control
coupling to support the discovery of non-obvious system designs that realize
new performance levels [2, 5, 12, 13, 14, 15]. Design coupling is accounted
for when plant and control design decisions are considered simultaneously in
the system optimization problem. Variations on the simultaneous CCD
formulation, such as nested or distributed optimization, have been studied
that provide numerical advantages over the simultaneous approach while being
mathematically equivalent [2, 5, 14, 15, 16].

The importance of CCD in wind energy research has gained recent attention
[17]. A holistic consideration of control and physical systems in a CCD
framework provides significant possibilities of effective and cost-efficient
wind turbine system design discoveries. The potential performance gains are
more pronounced in FOWTs than in land-based turbines due to more complex
couplings that arise when wind turbines are installed on floating platforms.
Various multidisciplinary design optimization (MDO) architectures, including
monolithic methods [18, 19] or nested methods [20] are implemented to explore
potential design advancements in wind energy systems. The results of these
studies demonstrated that the consideration for design coupling provides the
large advantages over designing for each disciplinary domain sequentially.

While previous studies have explored CCD for wind turbines, they have yet to
investigate the potential value of applying the CCD to FOWTs. Fairly
comprehensive CCD studies have been performed specifically for land-based wind
turbines but often relied on computationally-expensive models. Due to the
increased computational expenses beyond land-based models, many previous CCD
studies on FOWTs are limited to non-holistic approaches. While a handful of
articles present integrated design approaches for FOWTs, most utilize
sequential design approaches with open-loop optimal control (OLOC) or model
predictive control (MPC). This study fills a critical gap in integrated FOWT
system design and provides an initial foundation for future studies in this
area.

In this study, we demonstrate a CCD FOWT implementation that utilizes a
reduced-order model (ROM) that emulates the dynamic behaviors of a horizontal-
axis FOWT installed on a spar buoy floating platform. The baseline design
adapts specifications of the NREL 5 MW reference wind turbine [21] installed
on the OC3-Hywind spar buoy platform [22]. We use a nested CCD approach that
accounts for plant-control design coupling while capitalizing on specialized
solution strategies for direct optimal control. The outer-loop optimization
problem navigates the plant design space using the covariant matrix adaptation
evolution strategy (CMA-ES), which is a gradient-free population-based
optimization algorithm that utilizes a covariant matrix. The covariant matrix
adaptation procedure utilizes an estimation of the inverse Hessian for a
quadratic function to further evolve sampling shape toward the descent
direction [23, 24]. The inner-loop OLOC problem finds the optimal state and
control trajectories for a candidate plant design specified by the outer loop.
The OLOC problem is discretized using the pseudospectral direct optimal
control solver, GPOPS-II [10], and numerically solved using the interior-point
method NLP solver IPOPT [25].

Our unique contributions presented in this study include (1) a demonstration
of the holistic nested CCD design procedure for the FOWT system plant and
control design variables, (2) a computationally efficient FOWT model based on
reduced states and degrees-of-freedom (DOF) for system dynamics, including
neural network-based mooring line dynamics, (3) implementation of CMA-ES, a
gradient-free method, in the outer-loop plant design optimization to tackle
rugged numerical responses induced by the inner-loop OLOC solution, and (4)
analyses of the design solution that provide in-depth knowledge and insights
into the coupled effects of the design parameters and the FOWT system
behavior. Section 2 presents the problem definition, including the FOWT model,
design variables, and objective and constraint functions. Section 3 presents
the results and analysis of two CCD case studies: (1) optimization of tower
and control design variables and (2) optimization of the tower, blades, and
control design variables. Finally, Section 4 presents the concluding remarks
of this study.

##  2 Problem Definition

###  2.1 FOWT Dynamic Model and Plant Design Variables

FOWTs can utilize different floating platform types, including, but not
limited to, spar buoy, semi-submersible, tension leg platforms (TLPs), and
barge platforms. The studies presented here is based on a spar buoy platform.
Because many different disciplines and components are interacting with each
other in the system, deriving detailed dynamical equations is complicated. In
addition, as we increase the level of complexity, more advanced models are
needed, and as a result, computational time can increase significantly.
Simulating dynamic behaviors with such high complexity can be achieved using
sophisticated aero-hydro-servo-elastic simulation tools, such as OpenFAST [26,
27]. Working to solve CCD optimization problems by linking optimization
solvers directly to simulation tools such as OpenFAST can enhance design
problem fidelity, but presents several difficult challenges. For example,
direct transcription may not be possible to use with simulation tools that can
only be run as a black box, necessitating the use of the more limited single
shooting method. Recent work in the use of adaptive surrogate models for time
derivative functions has made possible the use of black-box simulation tools
with direct transcription, but is very limited in problem dimension [28]. A
dynamic system model is needed that provides direct access for the
optimization solver to the time derivative function. This need is addressed
here by constructing a reduced-order FOWT model that is compatible with CCD
methods based on direct transcription, while providing sufficient richness and
accuracy to support useful CCD optimization studies.

Several previous efforts have produced reduced-order FOWT models. For example,
Lemmer [29] developed a simplified model for FOWTs with semi-submersible
platforms, and Al-Solihat [30] developed a similar dynamic model for spar buoy
FOWTs. In this work, we have adapted the model published by Al-Solihat [30] to
obtain a spar buoy FOWT model that meets CCD optimization study requirements.
This implementation utilizes analytical dynamics to derive the equations of
motion in terms of platform quasi-coordinates. A two-dimensional dynamic
equation with an assumption of a rigid tower is introduced in Eq. (1):

$\displaystyle\begin{split}&\begin{bmatrix}M_{\text{sys}}+\bar{A}\end{bmatrix}_{4\times
4}\begin{bmatrix}\dot{v}_{x}\\\ \dot{v}_{z}\\\ \dot{\omega}_{y}\\\
\dot{\Omega}\end{bmatrix}_{4\times
1}+\begin{bmatrix}\tilde{S}M_{\text{sys}}+\bar{C}_{A}\end{bmatrix}_{4\times
4}\begin{bmatrix}v_{x}\\\ v_{z}\\\ \omega_{y}\\\ \Omega\end{bmatrix}_{4\times
1}\\\ &=\begin{bmatrix}F^{m}_{2\times 1}\\\ M^{m}_{1\times 1}\\\
0\end{bmatrix}_{4\times 1}+\begin{bmatrix}F_{2\times 1}\\\ M_{1\times 1}\\\
\tau_{\text{a}}-\tau_{\text{g}}\end{bmatrix}_{4\times 1}\end{split}$ (1)

where $M_{\text{sys}}$ is the mass matrix, $\bar{A}$ is the hydrodynamic added
mass matrix, $\bar{C}_{A}$ is the hydrodynamic Coriolis and centripetal
matrix, $v_{x}$ and $v_{z}$ are the platform velocities in the surge and heave
directions, respectively, $\omega_{y}$ is the platform pitch rate, $\Omega$ is
rotor rotational speed, $F^{m}$ is the gravitational force in the surge and
heave directions, $M^{m}$ is the gravitational moment in the pitch direction,
$F$ quantifies the external forces in the surge and heave directions, $M$ is
the external moment in the pitch direction, $\tau_{\text{a}}$ is the
aerodynamic torque, and $\tau_{\text{g}}$ is the generator torque. $\tilde{S}$
is also a matrix that is a function of $v_{x}$, $v_{z}$, and $\omega_{y}$, and
is presented in Eq. (2):

$\displaystyle\tilde{S}=\begin{bmatrix}0&\omega_{y}&0&0\\\
-\omega_{y}&0&0&0\\\ v_{z}&-v_{x}&0&0\\\ 0&0&0&0\end{bmatrix}$ (2)

$M_{\text{sys}}$ is defined in Eq. (3):

$\displaystyle M_{\text{sys}}=\begin{bmatrix}m_{\text{T}}&0&M_{13}&0\\\
0&m_{\text{T}}&-M_{26}&0\\\ M_{13}&-M_{26}&M_{55}&0\\\
0&0&0&I_{\text{r}x}\end{bmatrix}$ (3)

where $m_{\text{T}}$ is the total turbine mass and $I_{\text{r}x}$ is rotor
inertia. Additional model parameters are defined as:

$\displaystyle\begin{array}[]{lcl}m_{\text{T}}=m_{\text{p}}+m_{\text{t}}+m_{\text{nc}}+m_{\text{r}}\\\
M_{13}=D_{\text{r}}(m_{\text{r}}+m_{\text{nc}})+m_{\text{t}}D_{\text{t}}\\\
M_{26}=d_{\text{nc}}m_{\text{nc}}-d_{\text{r}}m_{\text{r}}\\\
M_{55}=I_{\text{T}_{y}}+m_{\text{r}}(D_{\text{r}}^{2}+d_{\text{r}}^{2})+m_{\text{nc}}(D_{\text{r}}^{2}+d_{\text{nc}}^{2})+D_{\text{t}}^{2}m_{\text{t}}\\\
I_{\text{T}_{y}}=I_{\text{p}y}+I_{\text{t}y}+I_{\text{nc}y}+I_{\text{r}y},\end{array}$
(9)

where $m_{\text{p}}$ is platform mass, $m_{\text{t}}$ is tower mass,
$m_{\text{nc}}$ is nacelle mass, $m_{\text{r}}$ is rotor mass, $I_{\text{p}y}$
is platform inertia, $I_{\text{t}y}$ is tower inertia, $I_{\text{nc}y}$ is
nacelle inertia, and $I_{\text{r}y}$ is rotor inertia. The dynamic equation
given in Eq. (1) is vectorized to efficiently calculate the time-domain the
simulation. $\bar{A}$ and $\bar{C}_{A}$ in Eq. (1) are defined as:

$\displaystyle\bar{A}$ $\displaystyle=\begin{bmatrix}A_{11}&0&A_{13}&0\\\
0&A_{22}&0&0\\\ A_{13}&0&A_{33}&0\\\ 0&0&0&0\end{bmatrix}$ (10)
$\displaystyle\bar{C}_{A}$ $\displaystyle=\begin{bmatrix}0&0&A_{11}v_{z}&0\\\
0&0&-A_{11}v_{x}&0\\\ -A_{11}v_{z}&A_{11}v_{x}&0&0\end{bmatrix}$ (11)

The parameters used in Eqs. (10) and (11) are defined as:

$\displaystyle\begin{array}[]{lcl}A_{11}=C_{\text{am}}\rho V_{\text{d}}\\\
A_{22}=C_{\text{am}}\left(\frac{1}{12}\rho\pi d_{1}^{3}\right)\\\
A_{13}=A_{11}(a_{\text{cv}}-a_{\text{pf}})\\\
A_{33}=C_{\text{am}}I_{\text{add}}\end{array}$ (16)

where $C_{\text{am}}$ is the hydrodynamic added mass coefficient, $\rho$ is
water density, $V_{\text{d}}$ is submerged volume, $d_{1}$ is the platform
base diameter, $a_{\text{cv}}$ is the distance between the platform base and
the submerged volume center, $a_{\text{pf}}$ is distance between the platform
Center of Gravity (COG) and the platform base, and $I_{\text{add}}$ is the
inertia of the displaced volume.

Gravitational forces and moment are defined in Eq. (17), where
$M^{\mathrm{m}}_{1\times
1}=g\sin{\theta_{p}}\left(D_{r}(m_{nc}+m_{r})+D_{t}m_{t}\right)+g\cos{\theta_{p}}(d_{nc}m_{nc}-d_{r}m_{r})$:

$\displaystyle\begin{bmatrix}F^{\mathrm{m}}_{2\times 1}\\\
M^{\mathrm{m}}_{1\times 1}\\\ 0\end{bmatrix}_{4\times 1}$
$\displaystyle=\begin{bmatrix}m_{T}g\sin{\theta_{p}}\\\
-m_{T}g\cos{\theta_{p}}\\\ M^{\mathrm{m}}_{1\times 1}\\\ 0\end{bmatrix}$ (17)

and other external forces and moments are shown in Eq. (18):

$\displaystyle\begin{bmatrix}F_{2\times 1}\\\ M_{1\times 1}\\\
\tau_{\mathrm{a}}-\tau_{\mathrm{g}}\end{bmatrix}_{4\times 1}$
$\displaystyle=\begin{bmatrix}F_{2\times 1}^{\mathrm{hs}}+F_{2\times
1}^{\mathrm{a}}+F_{2\times 1}^{\mathrm{moor}}+F_{2\times 1}^{\mathrm{hd}}\\\
M_{1\times 1}^{\mathrm{hs}}+M_{1\times 1}^{\mathrm{a}}+M_{1\times
1}^{\mathrm{moor}}+M_{1\times 1}^{\mathrm{hd}}\\\
\tau_{\mathrm{a}}-\tau_{\mathrm{g}}\end{bmatrix}$ (18)

where $\theta_{p}$ is platform pitch, $F^{\mathrm{hs}}$ is the hydrostatic
force, $F^{\mathrm{a}}$ is the aerodynamic force, $F^{\mathrm{moor}}$ is the
mooring force, $F^{\mathrm{hd}}$ is the hydrodynamic force, $M^{\mathrm{hs}}$
is hydrostatic moment, $M^{\mathrm{a}}$ is the aerodynamic moment,
$M^{\mathrm{moor}}$ is the mooring moment, and $M^{\mathrm{hd}}$ is the
hydrodynamic moment. These forces are calculated based on Al-Solihat [30] and
are not detailed here for brevity.

Figure 1: Schematic of the complete spar buoy-based FOWT system, including
spar platform, mooring system, tower, generator, and rotor.

(a) (b) (c) (d)

Figure 2: Schematics of system components, including: (a) tower structure, (b)
spar buoy platform, and (c) mooring line and a depiction its computational
model (d) that uses a neural network (NN) to model the horizontal and vertical
forces exerted by the mooring line.

(a) (b)

Figure 3: Schematics illustrating plant design variables: (a) tower design
variables; (b) blade design variables. Distributed blade design is
parameterized using twist ($\phi$) and chord length ($\zeta$) values at each
indicated ‘optimizing node’. Figure 4: XDSM of the nested CCD problem
depicting the solution process for the overall FOWT CCD problem. The inner-
loop problem solves for the state and control trajectories, while the outer-
loop problem solves for the plant design variables.

Figure 1 shows a schematic of the complete spar buoy FOWT system, including
all the system components considered in this study. In the figure,
$d_{\text{r}}$ is the distance between the rotor COG and tower center line,
$d_{\text{nc}}$ is the distance between the nacelle COG and the tower center
line, $D_{\text{t}}$ is the distance between the tower and platform COGs, and
$D_{\text{r}}$ is the distance between nacelle and platform COGs.

The complete system model presented in Fig. 1 comprised of the tower, spar
buoy floating platform, slack mooring lines, rotor, and generator models.
Figure 2 illustrates the internal structure of the FOWT components considered
in this study. Figure 2(a) depicts how the FOWT tower geometry is modeled as a
tapered hollow cylinder (aspect ratio is exaggerated for visualization), and
Fig. 2(b) shows the internal fixed and variable (water) ballast structures of
the spar buoy platform. Figures 2(c) and (d) depict the configuration and
computing algorithm for the mooring system. In this study, we implemented and
trained a multilayer perceptron neural network (NN) using newff and train
MATLAB commands to calculate mooring line dynamics, instead of directly
solving costly nonlinear equations. In Fig. 2(c), $l$, and $h$ are horizontal
and vertical distances between the anchor point and the fairlead,
$F_{\text{moor},H}$ is the horizontal mooring force, $F_{\text{moor},V}$ is
the vertical mooring force, and $W_{\text{moor},L}$ is the line weight. It is
required to use a different set of nonlinear equations to obtain the mooring
forces depending on whether the mooring line is suspended or resting on the
seabed [30]. When vertical force is greater than mooring line weight, the
mooring line is suspended; otherwise, it rests on the seabed. The neural
network is trained for each distinct case. Inputs of the neural network
include the relative position between the fairlead and the anchor point, and
the outputs are the vertical and the horizontal mooring forces.

Figure 3 presents schematics that help define plant design variables. In this
study the tower and blade designs are the extent of physical system design
decisions. Figures 3(a) and (b) illustrate the tower and blade design
variables, respectively. We defined four plant design variables for the tower:
tower thickness at the base $t_{\mathrm{base}}$, tower thickness at the tip
$t_{\mathrm{tip}}$, tower outer diameter at the tip $d_{\mathrm{tip}}$, and
tower length $l$. We also defined ten blade design variables. The blade is
divided into 17 nodes in the blade length-wise direction. Airfoil shape is
predetermined for each blade node, and is not subject to design in this study.
Distributed blade geometry may be adjusted by changing the values of blade
twist ($\phi$) and chord length ($\zeta$) for select control nodes. These
select nodes are depicted as ‘optimizing nodes’ in Fig. 3(b)), shown as black
circles. The first three nodes (gray nodes in Fig. 3(b)) are assigned fixed
values for twist and chord length. Twist and chord length values for
interpolated nodes (white nodes in Fig. 3(b)) are calculated by interpolating
functions based on a Bézier curve obtained by fixed and optimized node values.
Coordinates of the referenced airfoil shape profiles (CYL, DTU, and NACA) in
Fig. 3(b) can be found in the definition of NREL 5 MW floating offshore wind
turbine [21].

Figure 4 shows an extended design structure matrix (XDSM) of the nested CCD
optimization process. The objective function of the overall problem is
maximizing energy production, given as $E_{\text{in}}^{*}$. The tower and the
blade design variables are represented as
$\mathbf{x}_{\text{t}}=[t_{\text{tip}},\;d_{\text{tip}},\;t_{\text{base}},\;l]^{\textsf{T}}$
and $\mathbf{x}_{\text{b}}=[\phi_{5\times 1}^{\textsf{T}},\;\zeta_{5\times
1}^{\textsf{T}}]^{\textsf{T}}$, respectively. The superscripts $\dagger$ and
$*$ represent the updated and optimal values, respectively. In the outer loop,
the CMA-ES algorithm is used as an optimizer to iteratively create and update
the plant design variables. This algorithm was found to handle the noisy
response of the inner-loop problem efficiently and effectively. Power and
thrust coefficients ($C_{P}$ and $C_{T}$) were calculated using blade element
momentum (BEM) theory [31]. For each plant design the spar buoy platform
ballast mass required to achieve equilibrium in the heave direction is
calculated.

The inner-loop optimization problem is solved using the composition of these
plant-dependent parameters ($C_{P}$ and $C_{T}$), the dynamic model defined
above, and the mooring force NN model. For each candidate plant design
specified by the outer-loop plant design problem, the inner-loop problem is
solved, and the resulting optimal control trajectories and the inner-loop
objective function values for varied average wind speed points are sent back
to the outer-loop optimizer.

###  2.2 Objective and Constraint Functions

The objective function used here is to maximize the AEP. To achieve this goal
at the system-level, both the plant and control design variables need to be
concurrently optimized. In this study, the plant design variables influence
tower and blade geometry, and the control design variables are the time-
dependent generator torque and blade pitch trajectories. The goal is to
maximize the energy generated over a one-year long time horizon in the inner-
loop optimization problem. The annual power generation is calculated based on
a weighted sum of the sample simulations using a probabilistic wind
distribution. Each sample simulation consists of 100 seconds of a predefined
wind speed profile at varied average wind speed values following the
probabilistic distribution by setting the initial time, $t_{i}$, as 0 and the
final time, $t_{f}$, as 100 seconds in the inner-loop optimal control
problems. Using this weighted sum method, the optimal design (in terms of both
plant and control designs) converges to one that extracts the maximum average
power for wind speeds with higher probability of occurrence. In contrast, the
goal of maximum power can be relaxed for wind speeds with lower probability of
occurrence.

The AEP can be obtained using Eq. (19):

$\displaystyle\begin{split}\text{AEP}=8,760\int_{0}^{\infty}P_{\text{out}}(U)f(\bar{U})\,\mathrm{d}\bar{U}\\\
\approx
8,760\sum_{j=1}^{N_{\text{B}}}P_{\text{out}}(u_{j})P(\bar{u}_{j})\end{split}$
(19)

where 8,760 is the total number of hours in a year, $\bar{U}$ is the nominal
free wind velocity perpendicular to the rotor swept area, $U$ is the free wind
speed at the hub height, $P_{\mathrm{out}}$ is the time-averaged wind turbine
output power, and $f(U)$ is the Weibull probability density function. The free
wind speed $U$ is also a function of tower length since the height of the
tower determines the hub height. The nominal free wind speed, $\bar{U}$, is
based on the free stream velocity data at the hub height (90 meters) for the
NREL 5 MW reference turbine. The integral that calculates the AEP can also be
transcribed to a summation over $N_{\text{B}}$ wind speed bins. $\bar{u}_{j}$
and $u_{j}$ are $j$-th discretized quantities of $\bar{U}$ and $U$ among
$N_{\text{B}}$ wind speed bins, respectively.

Figure 5: Free wind velocity values and the corresponding frequency of
occurrence based on the Weibull probability distribution.

Here we use Weibull probability distribution as a wind speed probabilistic
occurrence model over one year. The probability density function of the
Weibull distribution is shown in Eq. (20):

$\displaystyle
f\left(\bar{U}\right)=\frac{k}{c}\left(\frac{\bar{U}}{c}\right)^{k-1}\exp\left(-\left(\frac{\bar{U}}{c}\right)^{k}\right)$
(20)

where $k>0$ is the shape parameter (here, $k=2$), and $c>0$ is the scale
parameter of the distribution (here $c=13.44$). This wind speed probabilistic
occurrence model is discretized with an interval of 1 m/s over a 3 to 25 m/s
range and adjusted to have a sum of 100% for its discretized range, resulting
in $N_{B}=23$ wind speed bins . The probability density function is
illustrated in Fig. 5. The wind load applied at the center of the tower is
also compensated by the power law to account for wind shear, shown in Eq.
(21):

$\displaystyle
u=u\left(\bar{u},l\right)=\bar{u}\left(\frac{l+z_{\text{base}}}{l_{\text{baseline}}+z_{\text{base}}}\right)^{0.2}$
(21)

This accounts for the tower height $l$ variation, where $l_{\text{baseline}}$
is the height of the NREL 5 MW baseline wind turbine tower (76 meters) and
$z_{\text{base}}$ is the distance from sea mean water level to tower base,
which is 12.4 meters. Figure 6(a) shows the wind velocity profiles for select
$\bar{u}$ values with the baseline tower length, and Fig. 6(b) shows the wind
velocity profiles as a function of tower length with a specific $\bar{u}$
value of 12 m/s. Increased tower length tends to exhibit increased AEP in
general because the free wind velocity at the higher altitude is larger, as
shown in Eq. (21). However, this trend only holds up to a specific limit in
practice because taller towers are also subject to increased thrust force,
which leads to an increased platform pitching motion. Increased platform pitch
amplitudes will eventually reach the pitch limits enforced by path constraints
in the inner loop optimal control problem. Excessive platform pitch can be
curtailed in this CCD optimization problem through less aggressive generator
torque or blade pitch control. This limits AEP gains with tower height. In
other words, at some value for tower height, platform pitch limits AEP gains.
Broader design changes, such as those to the platform or mooring subsystems,
could help increase the tower height at which platform pitch begins to
constrain AEP.

(a) (b)

Figure 6: Wind profiles in 100 seconds of time horizon for varied mean wind
velocity ($\bar{u}$) and varied tower height ($l$). (a) Wind profiles for mean
wind velocity values ranging from 6 to 16 m/s; (b) Wind profile for mean wind
velocity of 12 m/s and tower height values ranging from 60 to 100 m.

Optimization problems are often posed as minimization problems for
compatibility with established optimization solvers and by convention. The
outer loop objective function is defined as:

$\displaystyle
J_{\text{out}}=-8,760\sum_{j=1}^{N_{B}}P\left(\bar{u}\right)P_{\text{out}}\left(u\right)=-8,760\cdot
E_{\text{in}}$ (22) $\displaystyle\text{where:
}E_{\text{in}}=\sum_{j=1}^{N_{B}}P\left(\bar{u}\right)P_{\text{out}}\left(u\right)$
(23)

$E_{\text{in}}$ is obtained by solving the inner-loop problem across 23 wind
speeds with their corresponding probabilities, as shown in Fig. 5. The
objective of each inner-loop optimization problem is given in Eq. (24):

$J_{\text{in}}=-(t_{f}-t_{i})P_{\text{out}}\left(u\right)$ (24)

The time-averaged output power $P_{\text{out}}$ for each OLOC problem is
computed by Eq. (25):

$P_{\text{out}}\left(u\right)=\frac{1}{t_{f}-t_{i}}\int_{t_{i}}^{t_{f}}\left(P_{a}-10^{-7}\dot{\tau}_{g}^{2}-10^{7}\dot{\theta}_{b}^{2}\right)dt$
(25)

where $P_{\text{a}}$ is aerodynamic power, $\dot{\tau}_{\text{g}}$ is
generator torque rate, and $\dot{\theta}_{\text{b}}$ is blade pitch rate. The
generator torque and blade pitch rates are control trajectories that are
optimized in the inner-loop problem. The inner-loop objective function
includes control cost terms (i.e., small penalties on the torque and blade
pitch rates), to produce a non-singular optimal control problem and prevent
bang-bang or singular arc control trajectories in the solution. Also, as
discussed in Gros [32], the aerodynamic power is used as an integrand instead
of the generator power to prevent turnpike, which is embodied in turbine
problems as the tendency to maximize generator torque at the termination of a
finite time horizon. Parameter bounds and path constraint functions for the
inner-loop optimization problem are defined in Tables 1 and 2. The results
given in Section 3 indicate that some of these constraints are active in
certain regions of operation. Once all the inner-loop problems (in this study,
23 inner-loop OLOC problems) are solved, $E_{\text{in}}$ is updated by using
Eq. (23) for computing the outer-loop objective function, given in Eq. (22).

Table 1: Parameter bounds for the inner-loop OLOC problem.

|  |  | Lower | Upper
---|---|---|---|---
Name | Explanation | Unit | bound | bound
$\omega_{r}$ | rotor rotational velocity | rad/s | $0.00$ | 1.51
$\theta_{p}$ | platform pitch | deg | $-6.3$ | 6.3
$\theta_{b}$ | blade pitch | rad | $0.0$ | 40.0
$\tau_{g}$ | generator torque |  $10^{6}$ N-m | $0.00$ | 4.18
$\dot{\theta}_{b}$ | blade pitch rate | rad/s | $-0.57$ | 0.57
$\dot{\tau}_{g}$ | generator torque rate | MW | $-0.1$ | 0.1
$\phi$ | blade twist | deg | $-0.001$ | 15.970
$\zeta$ | blade chord length | m | $0.014$ | 5.580

Table 2: Path constraints for the inner-loop OLOC problem.

Name | Explanation | Unit | Constraint
---|---|---|---
$\sigma$ | Tower maximum static stress | MPa |  $\sigma\leq 45.0$ or $90.0$
$P_{u}$ | Generator power | MW | $P_{u}\leq 5.0$

In the nested CCD strategy introduced previously in Fig. 4, the outer loop
optimizer updates the plant design. In this study, we use CMA-ES to generate a
set of design points (population) to evaluate, and the covariant matrix
adaptation procedure finds the direction toward the optimal solution. For the
design points generated by the CMA-ES, the inner-loop OLOC problem is solved.
The Weibull probabilistic distribution is utilized to divide wind speeds into
23 bins with center $\bar{u}\in\mathcal{U}$, where $\bar{u}$ represents the
average wind velocity. As defined in Eq. (21), the wind speed for defining
wind load $u$ is a function of $\bar{u}$ and $l$, and this relationship
depends on tower height. The outer-loop objective function is equivalent to
the weighted sum of the $N_{B}$ inner-loop objective functions, which is
required for stable convergence of the nested CCD problem. When the outer-loop
optimization problem converges, the CMA-ES algorithm provides the final
optimal solution and will be terminated.

##  3 Results and Discussion

This section comprises several studies aimed at showcasing the results of the
CCD approach and gaining insights from different scenarios. In Section 3.1,
two cases are discussed: 1) CCD with four tower parameters as plant design
variables, and 2) CCD with four tower parameters and ten blade parameters as
plant design variables. This section presents the obtained optimal plant
variables and objective values. It also includes an illustration of the CMAES
population, and a sensitivity study is performed. Section 3.2 investigates the
impact of changing the maximum allowable tower stress constraint on power
generation and Annual Energy Production (AEP). It should be noted that here,
stress resulting from fore-aft bending was only considered. Throughout this
paper, when discussing stress upper bound, we specifically mean the upper
bound of fore-aft bending stress. This constraint is used here just as a means
to observe the coupling between plant and controller, and to see how the
design changes as we modify the upper bound. In future work, additional
constraints and higher fidelity models will be employed. Section 3.3 explores
the effect of increasing tower mass on the time series trajectories within the
inner loop. Lastly, Section 3.4 demonstrates the inner loop trajectories in
the presence of waves and includes a fatigue study. These sections
collectively offer a comprehensive analysis of CCD results and provide
valuable insights into various aspects of the design and optimization process
of floating offshore wind turbines.

Table 3: Optimal plant design variable values and corresponding merit function
values for optimization with varied targets. The tower-only (T) CCD problem
solves for tower design variables in conjunction with control trajectories.
The tower and blades (T&B) CCD problem solves for the tower design parameters,
blade design parameters, and control trajectories. Values with gray shading
represent parameters are not optimized and held constant.

|  |  | Optimization target
---|---|---|---
Group | Variable | Unit |  None (baseline) |  Tower only |  Tower & blades
Tower parameters | $t_{\text{base}}$ | m | 0.027 | 0.042 | 0.042
$t_{\text{tip}}$ | m | 0.019 | 0.012 | 0.012
$d_{\text{tip}}$ | m | 3.870 | 4.95 | 5.00
$l$ | m | 77.600 | 83.04 | 81.38
Blade parameters | $\phi_{4}$ | deg | 13.31 | 13.31 | 12.09
$\phi_{6}$ | deg | 11.48 | 11.48 | 9.23
$\phi_{9}$ | deg | 6.54 | 6.54 | 3.96
$\phi_{12}$ | deg | 1.53 | 1.53 | 1.52
$\phi_{17}$ | deg | 0.11 | 0.11 | 0.06
$\zeta_{4}$ | m | 4.557 | 4.557 | 5.580
$\zeta_{6}$ | m | 4.007 | 4.007 | 5.056
$\zeta_{9}$ | m | 3.502 | 3.502 | 3.136
$\zeta_{12}$ | m | 2.764 | 2.764 | 2.458
$\zeta_{17}$ | m | 1.419 | 1.419 | 1.458
Merit functions | $-J^{*}_{\text{out}}$ | GWh | 26.24 | 29.22 | 29.34
$m_{\text{tower}}$ | tonne | 249.6 | 344.4 | 342.0

(a) (b)

Figure 7: Optimized tower and blade designs for the CCD problem with both
tower and blade (T&B) variables optimized. (a) Optimized tower shape compared
to the baseline design; (b) Optimized blade twist ($\phi$) and chord length
($\zeta$) over normalized blade radial position ($r/R$). The tower wall
thickness shapes are exaggerated for visual comparison. Readers are referred
to ‘Tower & blades’ column in Table 3 for accurate wall thickness values.

(a) (b)

Figure 8: Optimized blade power and thrust coefficients as functions of tip
speed ratio ($\lambda$) and blade pitch angle ($\theta_{b}$). (a) Power
coefficient surface; (b) thrust coefficient surface.

###  3.1 Control Co-Design With Tower Only and Tower and Blade Plant Design
Variables

This section presents the results of two CCD studies and compares them with
the optimal control result of the baseline NREL 5 MW wind turbine with the
same model presented in Section 2. Unless mentioned otherwise, the maximum
tower stress constraint value used in this study is 45 MPa. The summary of the
problems we compare in this study are listed in Table 3, where $\phi$ is blade
twist, $\zeta$ is blade chord length, and indices represent blade node values.
Readers are also referred to Fig. 3 for the nodal information where design
variables are located. Gray-shaded cells in the table show which design
variables are held fixed and are not optimized during the CCD procedure. The
‘None (baseline)’ column represents the NREL 5 MW baseline result. All plant
design variables are shaded (fixed parameters), and the inner-loop OLOC
problem is the only design problem that is solved in this case to generate the
optimal AEP, $J_{\text{out}}^{\ast}$. The next ‘Tower only’ column presents
the results from the first CCD problem, which optimizes the tower design
variables and the control trajectories. The last ‘Tower & blades’ column
corresponds to the second CCD problem, which optimizes both the tower and
blade design variables, as well as the control trajectories.

As we see from the overall AEP (objective function) values,
$-J_{\text{out}}^{*}$, from Table 3, increased design flexibility with more
plant design degrees of freedom help improve the objective function value
significantly. Compared to the baseline design, optimizing the tower and blade
design variables and the control trajectories together increases AEP by 7.17%.
Figure 7(a) depicts the optimized tower design for the CCD case with both the
tower and blade design variables optimized. The converged design solution
indicates that the tower height has increased to a particular height (here,
83.12 m) that is larger than the baseline tower height (77.60 m) design. In
addition, the tower thickness has also increased to satisfy tower stress
constraints. The optimized blade is illustrated in Fig. 7(b), and the
resulting power and thrust coefficients corresponding to this optimal blade
design are shown in Fig. 8. The power and thrust coefficients are functions of
blade pitch $\theta_{\text{b}}$ and tip speed ratio $\lambda$.

(a) (b)

Figure 9: Explored designs in the two-dimensional space of AEP and mass of
FOWT components for the CCD problem with tower-only (T) plant optimization.
Square markers represents the optimal plant and control design points. Diamond
markers represent a selected suboptimal design point with a tower design that
is heavier than the optimal design. Circular dots represent explored design
points. Here, platform design is not optimized, but ballast mass is adjusted
to keep the mean water level of the floating platform consistent. (a) Explored
designs in AEP and tower mass space; (b) Explored designs in AEP and platform
mass space. Figure 10: Resulting tower shape for the two design points
highlighted in Fig.9. The selected suboptimal design point has a larger tower
mass, but the AEP for the optimal design point (29.22 GWh) is about 9% larger
than the AEP for the suboptimal design point (26.77 GWh).

Candidate designs that were explored during the CMA-ES optimization procedure
are shown in Fig. 9, depicted in two distinct two-dimensional spaces. The
corresponding tower shapes for the optimal design and a selected suboptimal
design are shown in Fig. 10. In Fig. 9, the optimal design point is indicated
with a green square, and a sample suboptimal design point with a slightly
larger tower mass is marked with a blue diamond. The indicated optimal design
point exhibits the highest AEP value among all explored design points. The
zoomed-in view of the plot shows that the optimal point has one of the largest
objective function values among all of its neighbors. Across all of the
designs explored in this study, the platform shape is fixed. As tower mass
changes, the variable ballast (ballast water mass in practice) is adjusted to
satisfy the force equilibrium requirement in the heave direction. Therefore,
the platform capital cost for all these design points will be the same, but
the tower cost varies with tower design changes.

Suppose we consider a multiobjective optimization problem for maximizing the
AEP and minimizing the tower mass. This quantifies the best-possible tradeoffs
between AEP and mass. In this case, all the design points located on the right
side of the hypothetical vertical line crossing through the optimal design
point are inferior to the optimal design since both performance indices (the
AEP and the tower mass) are dominated by the optimal design. However, on the
left side of the vertical line, finding good nondominated solutions in terms
of superior tower mass value by sacrificing the AEP value is possible. Thus,
we have identified a strategy for expanding this CCD problem into a
multiobjective optimization study. Additional objective functions that
consider additional performance and cost indices beyond AEP and mass may be
used.

Table 4: Sensitivity study results with $\pm 5\%$ perturbations of tower
design variable values obtained from the tower-only CCD optimization result.

Tower design | $\Delta_{\mathrm{rel}}$ | Value | $J_{\text{out}}$ | $\Delta J_{\mathrm{rel}}$ | $\Delta J_{\mathrm{rel}}/\Delta_{\mathrm{rel}}$
---|---|---|---|---|---
parameter | [%] | [m] | [GWh] | [%] | [-]
$t_{\text{base}}^{+}$ | $+5.0$ | 0.0445 | $-28.7839$ | $-1.49$ | $-0.298$
$t_{\text{base}}^{-}$ | $-5.0$ | 0.0403 | $-29.2094$ | $-0.04$ | $-0.008$
$t_{\text{tip}}^{+}$ | $+5.0$ | 0.0120 | $-29.2083$ | $-0.04$ | $-0.008$
$t_{\text{tip}}^{-}$ | $-5.0$ | 0.0108 | $-26.9892$ | $-7.63$ | $-1.526$
$d_{\text{tip}}^{+}$ | $+5.0$ | 5.353 | $-29.2076$ | $-0.04$ | $-0.008$
$d_{\text{tip}}^{-}$ | $-5.0$ | 4.843 | $-29.2134$ | $-0.02$ | $-0.004$
$l^{+}$ | $+5.0$ | 86.751 | $-29.1929$ | $-0.09$ | $-0.018$
$l^{-}$ | $-5.0$ | 78.489 | $-29.2023$ | $-0.06$ | $-0.012$

A sensitivity analysis has been performed with a $\pm 5\%$ perturbation of
tower design variables based on the optimal design resulting from the tower-
only CCD problem. The sensitivity results confirm the validity of our design
solution by demonstrating performance reductions. Table 4 presents a summary
of the sensitivity analysis results. In this table, a relative perturbation
$\Delta_{\text{rel}}$ for an arbitrary quantity $\chi$ is defined as
$\Delta_{\text{rel}}=\left(100\%\right)(\chi-\chi^{\ast})/\chi^{\ast}$.
Similarly, a relative change in AEP $\Delta J_{\text{rel}}$ is defined as
$\Delta
J_{\text{rel}}=\left(100\%\right)(\left|J_{\text{out}}\right|-\left|J^{\ast}_{\text{out}}\right|)/\left|J^{\ast}_{\text{out}}\right|$.
As observed in the results, varying any tower design variable by 5% in any
direction from the optimal design decreases the AEP. In addition, the table
indicates that the tower tip thickness and tower base thickness are the most
sensitive parameters that affect the objective function value.

###  3.2 Comparison of Increased Maximum Allowable Tower Stress

Figure 11: Optimal tower shapes for two distinct allowable tower stress
values: 45 MPa and 90 MPa

(a)

(b)

(c)

(d)

(e)

(f)

Figure 12: Generated power curves under various scenarios. Sub-figures refer
to power curves of steady and varied wind cases for (a)
$x_{p}^{\ast}|_{\sigma\leq 90},\,\sigma\leq 90$; (b)
$x_{p}^{\ast}|_{\sigma\leq 45},\,\sigma\leq 45$; (c)
$x_{p}^{\ast}|_{\sigma\leq 90},\,\sigma\leq 45$; Comparison of steady wind
power curves with varied tower stress constraints in the control problem, when
optimal plant design is obtained at (d) $\sigma\leq 90$ MPa; (e) $\sigma\leq
45$ MPa; and (f) Comparison of AEP trends for cases illustrated in (d) and
(e).

To quantify the effect of failure mode constraints on the CCD result, a new
study is presented in Figs. 11–12. In this study, the maximum allowable tower
stress is increased from 45 MPa to 90 MPa. This corresponds to a different
(more expensive) tower material selection. Figure 11 shows the optimal tower
shapes obtained for these two cases. Increased allowable stress results in
optimal designs with increased tower heights, increased wind load, and a
higher AEP.

The study involves comparing power curves generated using two different stress
upper bounds. Additionally, examining the changes in the power curve when the
plant design is obtained based on one case, but the OLOC is solved assuming a
maximum stress constraint from the other case would provide insightful
results. This concept is common in many engineering situations, where plants
are designed either using CCD or sequential approaches with certain
constraints, but the maximum constraint may change in reality. This study
demonstrates that such changes can degrade system response and performance. In
this context, ${x^{\ast}_{p}}|_{\sigma\leq\sigma_{i}}$ represents the optimal
plant obtained using the CCD approach with a maximum stress constraint of
$\sigma_{i}$ MPa. Subsequently, $\sigma\leq\sigma_{j}$ indicates that the
power curve is generated using the optimal plant from the previous step, while
the maximum stress in the OLOC simulation is set to $\sigma_{j}$ MPa.

Figure 12(a) illustrates the power curve obtained when the maximum tower
stress is set to 90 MPa. In this case, ${x^{\ast}_{p}}|_{\sigma=90}$ refers to
the optimal plant obtained using the CCD approach with a maximum stress
constraint of 90 MPa. Subsequently, $\sigma=90$ indicates that the power curve
is generated using the optimal plant from the previous step, while the maximum
stress in the OLOC is also set to 90 MPa. In part (a) of the figure, both the
CCD and OLOC utilize the same maximum stress constraint to generate the power
curve, indicating an expected favorable outcome. The power curve is generated
using a constant wind speed ranging from 3 to 25 m/s in 1 m/s increments. The
simulation is conducted over a duration of 10 minutes, with the average power
value computed. Additionally, to examine the power time series for non-
constant wind speeds, the wind profile depicted in Fig. 6 is employed,
incorporating average values between 3 and 25 m/s. As depicted in Fig. 12(a),
the power curve exhibits an increase roughly proportional to the cube of the
wind speed in region 2 (below rated wind speed). In region 3 (above rated wind
speed), the power remains relatively constant (assuming power constancy
instead of torque constancy). Moreover, the power time series aligns closely
with the average values derived from the constant wind speeds.

Figure 12(b) depicts a similar plot, but in this case, both the CCD and OLOC
simulations employ a constant stress of 45 MPa. Figure 12(c) illustrates the
results when the plant is obtained from the CCD case with a maximum stress of
90 MPa, but in the subsequent OLOC simulation, the stress constraint is
reduced to 45 MPa. As observed, the power curves show a decline, which is
reasonable considering the disparity between the maximum stress assumptions in
the CCD and OLOC simulations.

In Fig. 12 (d), the power curve is presented for various cases. For the first
six cases, the plant design is obtained from ${x^{\ast}_{p}}|_{\sigma=90}$,
but the maximum stress in the OLOC simulation is varied from 45 to 90 MPa. It
can be observed that decreasing the maximum stress in the OLOC simulations
causes a downward shift in the power curve. Furthermore, as the maximum stress
is reduced further, the power decrease becomes more pronounced. The last case
in Fig. 12 (d) corresponds to ${x^{\ast}_{p}}|_{\sigma=45}$, and
interestingly, the power curve for this case is nearly identical to that of
${x^{\ast}_{p}}|_{\sigma=90}$. The only distinction between these two cases
lies in the assumption regarding the maximum stress in the CCD and OLOC
simulations. In the last case, the maximum stress remains the same in both
simulations, while in the previous case, the plant obtained from
${x^{\ast}_{p}}|_{\sigma=90}$ is designed to tolerate stresses ranging from 45
to 90 MPa. However, enforcing a maximum stress of 45 MPa in the OLOC
simulation leads to performance degradation, requiring a compromise in power
generation to meet the stress constraints.

Figure 12 (e) illustrates a similar concept, but this time the plant design is
obtained from ${x^{\ast}_{p}}|_{\sigma=45}$ for the first six cases. As
observed, increasing the maximum stress in the OLOC simulation does not alter
the power curve because the plant design is already based on an assumption of
a maximum stress of 45 MPa. Therefore, increasing this constraint in the OLOC
simulation has no effect on the power curve.

Figure 12 (f) displays the Annual Energy Production (AEP) derived from parts
(d) and (e) of the analysis. For ${x^{\ast}_{p}}|_{\sigma=90}$, reducing the
maximum stress constraint results in a degradation of the power curve and
consequently a decrease in AEP. However, for ${x^{\ast}_{p}}|_{\sigma=45}$,
increasing the maximum stress bound does not alter the power curve
significantly, and the AEP remains almost the same. Upon examining this
figure, one may question the purpose of employing a higher maximum stress
constraint in the CCD if the power curve obtained from a 45 MPa assumption is
not considerably lower than the 90 MPa case. The response to this query is
two-fold:

  1. 1.

As illustrated in Fig. 10 the tower designed for ${x^{\ast}_{p}}|_{\sigma=45}$
is thicker, potentially increasing the tower mass and subsequently the total
cost. In this study, for instance, the tower mass increased from 244.9 tonnes
in ${x^{\ast}_{p}}|_{\sigma=90}$ to 344.4 tonnes in
${x^{\ast}_{p}}|_{\sigma=45}$.

  2. 2.

Even a small enhancement in AEP can yield significant monetary benefits. In
the present case, the difference amounts to approximately 0.5 GWh. Assuming an
electricity cost of $\$0.13$ per kWh, this translates to an annual savings of
approximately $\$65,000$ USD.

Overall, the consideration of factors such as tower mass and cost, as well as
the potential financial gains from even slight improvements in AEP, justifies
the exploration of higher maximum stress constraints in the CCD optimization
process.

###  3.3 Comparison Studies Showing Tower Mass Effect on AEP

Table 5: Comparisons of parameters and results for OLOC studies with optimal tower mass (Case 1), increased tower mass (Case 2), and forward ordinary differential equation (ODE) simulation result using the control trajectories of Case 1 with the tower design of Case 2 (Case 3). $P_{\text{out}}$ is reported for the case where wind speed is equal to 14 m/s. Parameters | Case 1 | Case 2 | Case 3
---|---|---|---
$t_{\text{base}}\text{ [m]}$ | 0.042 | 0.042 | 0.042
$t_{\text{tip}}\text{ [m]}$ | 0.012 | 0.028 | 0.028
$d_{\text{tip}}\text{ [m]}$ | 4.95 | 6.50 | 6.50
$l\text{ [m]}$ | 83.04 | 83.04 | 83.04
$m_{\text{t}}\text{ [kton]}$ | 0.345 | 0.505 | 0.505
$m_{\text{p}}\text{ [kton]}$ | 7.472 | 7.311 | 7.311
$P_{\mathrm{out}}\text{[MW]}$ | 4.93 | 4.83 | 4.83
AEP [GWh] | 29.22 | 28.73 |
Method | OLOC | OLOC | Simulation
Figure 13: Trajectories of control inputs and model outputs that correspond to
the three cases listed in Table 5. Dotted horizontal lines indicate the upper
bounds for corresponding quantities. (a) Generator torque; (b) blade pitch;
(c) Rotor rotational speed; (d) Generator power; (e) Power coefficient.

Another set of comparison studies was performed to show the effect of changing
tower mass on AEP while maintaining a consistent tower height. Table 5
contains the summary of parameter values and the results of this comparison
study. Case 1 is the OLOC problem with the optimal tower design parameters
obtained via the CCD problem presented above. Case 2 is another OLOC problem
with the same tower height but with increased tower mass. Case 3 is based on
solving the dynamic system model (ordinary differential equations) using a
forward simulation and the control trajectories obtained from Case 1 and the
tower design of Case 2. All three problems are solved using the same wind
profile with an average wind speed $\bar{u}$ of 14 m/s. In summary, Cases 1
and 2 are solved for obtaining the optimal control trajectory designs, while
Case 3 is solved using a forward simulation with provided control
trajectories.

Comparisons of results for these three cases are shown in Figs. 13–14. Figure
13(a) and (b) shows the generator torque and blade pitch trajectories, which
are controlled quantities. Actual control signals used in numerical
optimization are rates of these quantities. The control trajectories for Case
3 are extracted from the OLOC result of Case 1. However, due to a difference
in numerical integration schemes, generator torque trajectories exhibit a
slight difference. Case 1 uses the Gaussian quadrature method for the
numerical integration, while Case 3 uses the Dormand and Prince [33] method, a
family of higher-order Runge-Kutta methods. Notably, the optimal control
trajectories for Cases 1 and 2 are significantly different since the tower
designs for these cases are different. This difference indicates that the
control strategy should be individually established for a specific plant
design, which could be attainable using the CCD strategy presented in this
study. Feedback control laws could be synthesized based on the OLOC CCD
results. It is possible that for some CCD problems, no feedback control exists
that can mimic the behavior of the OLOC trajectories. Methods for
systematically informing CCD optimization studies based on feedback control
law limitations is a topic of ongoing work.

Figure 14: Trajectories of model dynamic outputs, including motion
trajectories, stress, and velocities that correspond to the three cases listed
in Table 5. Dotted horizontal lines indicate the upper bounds for
corresponding quantities. (a) Platform surge motion; (b) platform heave
motion; (c) platform pitch motion; (d) tower stress; (e) rotor hub fore-aft
velocity; (f) relative wind velocity.

Figure 13(c)–(e) shows the trajectories of rotor rotational velocity,
generator power, and power coefficient curves, and Fig. 14 shows the
trajectories of model dynamic responses, including platform motions in
(a)–(c), tower stress in (d), and velocities in (e)–(f). The plant design for
Case 1 is selected from the ‘Tower only’ CCD result, which is already
optimized and is given in Tab. 3. In contrast, the plant design for Case 2
deviates from the optimal design since tower mass was intentionally increased.
Thus, in Fig. 13(d), the area under the generator power curve $P_{\text{u}}$
for Case 1 is significantly larger than for Case 2. Also, Case 1 exhibits a
much longer period of constraint saturation (active path constraint indicated
by the dotted horizontal line) compared to Case 2, as shown in Fig. 13(d). In
addition, because the control trajectories used in Case 3 are not optimized
with respect to its unique plant design parameters, we see that the Case 3
simulation result violates bounds and path constraints, as shown in Figs.
13(c)–(d) and 14(c).

From the results given in Figs. 13–14, increased tower mass produces a larger
tower mass moment of inertia with respect to the rotational center for the
platform pitch DOF. Thus, Case 3 exhibits a larger value in the peak platform
pitch motion, as shown in Fig. 14(c). As Cases 2 and 3 have common plant
designs, the OLOC result from Case 2 can provide a mitigating strategy for the
increased platform pitch motion from Case 3. By comparing the control
trajectories of Cases 2 and 3, we observe that the blade pitch motion in Case
2 is significantly increased, resulting in reduced aerodynamic thrust applied
to the turbine. With this optimized control strategy, the platform pitch
motion can be mitigated to satisfy the constraint with the increased tower
mass in Case 2, but with a sacrifice in power generation.

We used the same absolute wind profiles at the same mean wind velocity value
($\bar{u}=14$ m/s) for fair comparisons across these three cases. One item to
note is that the relative wind velocity profiles for these three cases are
different from each other, as shown in Fig. 14(f). The reason behind this
difference is that the thrust load exerted to the wind turbine is based on the
relative wind velocity considering the rotor hub fore-aft motion. Thus, the
relative wind velocity is calculated by subtracting the rotor hub fore-aft
velocity from the absolute wind velocity. Due to significant platform pitch
motion (unlike land-based wind turbines), the rotor hub fore-aft velocity is
significant. As a result, the relative wind velocity profile can vary
depending on plant design, and this variation impacts the resulting optimal
control strategy and power generation.

It should also be noted that the platform shape across all studies presented
here is fixed. The only variable here that can change the platform properties
is the variable ballast mass. When the tower design parameters are modified,
the variable ballast mass is automatically modified to satisfy the mass and
buoyancy equilibrium requirement in the platform heave DOF. However, there are
a few other ways to achieve mass and buoyancy equilibrium: (1) adjust the
submerged volume of the platform by changing shape, or (2) adjust equilibrium
draught height. Each of these options may affect optimization results
significantly. However, to maintain simplicity, we kept the platform shape and
the equilibrium draught height as fixed constants. Instead, we defined the
variable ballast mass as a variable that depends on the tower mass. More
flexible design studies may reveal new insights and potentially higher
performance.

This study highlights an important finding that increasing tower mass does not
necessarily result in an increase in generated power and Annual Energy
Production (AEP). In fact, it was observed that for a fixed tower length,
increasing the tower diameter and tip thickness led to a worse outcome in
terms of power generation. This underscores the need for an integrated
approach in obtaining optimal plant design variables, considering the coupling
between the design variables and the controller. By taking full advantage of
this coupling, it becomes possible to optimize both the structural aspects and
the controller, leading to improved overall performance and energy production.

###  3.4 Wave and Fatigue Study

Figure 15: Optimal control performed on the “tower only” plant design given in
Tab. 3 without and with additional wave loading. (a) Wave amplitudes; (b) Wave
velocities in surge and heave directions; (c) Platform heave motions; (d)
Platform pitch motions; (e) Optimal generator torque trajectories; (f) Tower
stress trajectories.

In the previous sections, the dynamics of the FOWT system took into account
the interaction between the water and the platform motion, but wave loads were
excluded for simplicity. However, wave loading is an essential external
disturbance that can cause FOWT system failures. Wave loading leads an
additional computational burden in solving the CCD problem of the FOWT system.
Thus, instead of solving all CCD cases with wave loading, an optimal result
obtained for the tower-only case with wave loading is compared to the result
of the same problem without wave loading obtained in the previous section.
Fatigue analysis is also performed on the tower-only CCD cases (with and
without wave loading).

Figure 15 shows the optimal control results for two cases: without and with
wave loading. The wave is obtained from Al-Solihat [30] and is a irregular
signal with a period ($T_{p}$) of $10$ s and significant wave height ($H_{s}$)
of $6$ m. Including wave loading results in higher frequency motion responses,
requiring higher frequency control inputs. This may accelerate fatigue damage.
To see this effect, Mlife [34], a fatigue prediction code, is used to compute
fatigue life of the tower. The lifetime is considered 20 years, and the
ultimate strength is set in such a way that 20 years life time is satisfied.
For the case that wave loading is not included, 20 years of lifespan is
satisfied with the tower design strength of $45.08$ MPa, which is close to the
previous CCD result without considering fatigue lifespan ($45$ MPa). However,
when the wave loading is implemented, the required tower design strength needs
to be increased to $61.78$ MPa to satisfy the 20 years of lifespan. It should
be noted that this is an initial study utilizing low-order models to decrease
computation time. More detailed models need to be incorporated to verify
fatigue lifespan of the entire FOWT system design.

##  4 Conclusion

In this article, a nested CCD method using OLOC is implemented to solve an
integrated FOWT design optimization problem. The reduced order FOWT dynamic
model is based on a spar buoy floating platform and the NREL 5-MW reference
turbine. Tower and blade designs are parameterized using 14 plant design
variables for the outer-loop problem. The time-domain OLOC problem is
formulated using the pseudospectral method for the inner-loop problem. The
objective function for this nested CCD problem is to maximize the AEP, subject
to constraints defined to enforce stress, motion, and other various operating
limits. Neural network-based mooring line dynamic models are trained for
rested and suspended mooring line conditions to facilitate efficient
computation for use with optimization studies. The simplified reduced-order
modeling approach enabled efficient CCD studies that are practical to solve
using standard computational resources. At the same time, the model fully
captures all of the essential physics behind the dynamic behaviors of the
FOWT, especially those that give rise to key design coupling relationships.

The CCD result shows an increase in the AEP compared to the baseline design by
more than eleven (11) percent. In addition, we demonstrated a possibility to
expand this CCD problem with various additional objective functions, including
tower mass, to form a multiobjective optimization problem to help reduce the
cost of the FOWT system further. Influences of the design variables on the
objective function are investigated using sensitivity analysis. By comparing
the inner-loop OLOC solutions for increased tower mass design, we identified a
dynamic mechanism that leads to a certain optimal tower mass that produces
maximal power. Due to platform pitch motion, there exists a sweet spot for the
optimal tower mass and height that maintains structural integrity and limits
platform motion while supporting large wind speeds at high rotor hub
elevations.

Furthermore, as explored in Section 3.2, it was demonstrated that the power
curve and AEP are influenced by the maximum tower stress constraint used in
CCD. It was observed that if the maximum stress in the OLOC simulation is
lower than the assumption made in CCD, it results in a degradation of the
power curve and AEP. Conversely, if the maximum stress in the simulation is
higher than the assumption made in CCD, the power curve remains relatively
unchanged, but it leads to an increase in tower mass and consequently higher
capital costs. Additionally, even when the power curves are closely aligned,
even a small difference between them can translate into substantial financial
discrepancies. Therefore, these types of constraints in CCD should be set in a
manner that accurately represents real-world conditions, with an added margin
of safety.

The nested CCD implementation presented here (based on a simplified reduced-
order FOWT dynamic model) is appropriate for early-stage FOWT design studies
in practice. CCD studies based on this implementation can reveal overall
preferred design directions for system designers to explore. Also, this
implementation can be used as a template for developing more advanced CCD
approaches that incorporate higher fidelity FOWT system models, such as
OpenFAST [27] simulation code, as well as enhanced design
flexibility/representation fidelity. An additional ongoing effort to improve
the practicality of similar CCD implementations is to adapt closed-loop
control (CLC) or MPC design capabilities in the inner-loop optimal control
problem. In this way, the ideal OLOC design results can inform promising
directions for engineers to explore for more practical CLC or MPC-based
controllers.

## Funding

The information, data, or work presented herein was funded in part by the
Advanced Research Projects Agency-Energy (ARPA-E), U.S. Department of Energy,
under a project titled: “WEIS: A Tool Set to Enable Controls Co-Design of
Floating Offshore Wind Energy Systems” associated with award no DE-
AC36-08GO28308.

## Declaration of Competing Interest

The authors declare that they have no competing interest.

## Data Availability

All data required to replicate the results can be generated by the MATLAB
optimization code. A MATLAB optimization code for a select problem
demonstrated in the manuscript are available upon request to the first or
corresponding authors. Readers may need a commercial pseudospectral optimal
control solver, GPOPS-II, to run the provided MATLAB codes.

## References

  * Bhattacharya [2014] S. Bhattacharya,  Challenges in design of foundations for offshore wind turbines,  Engineering & Technology Reference (2014). doi:10.1049/etr.2014.0041.
  * Allison et al. [2014] J. T. Allison, T. Guo, Z. Han,  Co-design of an active suspension using simultaneous dynamic optimization,  Journal of Mechanical Design 136 (2014) 081003. doi:10.1115/1.4027335.
  * Allison and Herber [2014] J. T. Allison, D. R. Herber,  Multidisciplinary design optimization of dynamic engineering systems,  AIAA Journal 52 (2014) 691–710. doi:10.2514/1.J052182.
  * Böhme and Frank [2017] T. J. Böhme, B. Frank,  Indirect methods for optimal control,  in: Hybrid Systems, Optimal Control and Hybrid Vehicles, Advances in Industrial Control, Springer, Cham, 2017, pp. 215–231. doi:10.1007/978-3-319-51317-1\\_7.
  * Herber and Allison [2019] D. R. Herber, J. T. Allison,  Nested and simultaneous solution strategies for general combined plant and control design problems,  Journal of Mechanical Design 141 (2019) 011402. doi:10.1115/1.4040705.
  * Bayat and Allison [2023a] S. Bayat, J. T. Allison,  SS-MPC: A user-friendly software based on single shooting optimization to solve model predictive control problems,  Software Impacts 17 (2023a) 100566.
  * Bayat and Allison [2023b] S. Bayat, J. T. Allison,  Control Co-Design with varying available information applied to vehicle suspensions,  in: International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American Society of Mechanical Engineers, 2023b.
  * Betts [1998] J. T. Betts,  Survey of numerical methods for trajectory optimization,  Journal of Guidance, Control, and Dynamics 21 (1998) 193–207. doi:10.2514/2.4231.
  * Ross and Karpenko [2012] I. M. Ross, M. Karpenko,  A review of pseudospectral optimal control: From theory to flight,  Annual Reviews in Control 36 (2012) 182–197. doi:10.1016/j.arcontrol.2012.09.002.
  * Patterson and Rao [2014] M. A. Patterson, A. V. Rao,  GPOPS-II: A MATLAB software for solving multiple-phase optimal control problems using hp-adaptive Gaussian quadrature collocation methods and sparse nonlinear programming,  ACM Transactions on Mathematical Software 41 (2014) 1–37. doi:10.1145/2558904.
  * Garg et al. [2009] D. Garg, M. A. Patterson, W. W. Hager, A. V. Rao, D. A. Benson, G. T. Huntington,  An overview of three pseudospectral methods for the numerical solution of optimal control problems,  in: A. V. Rao, T. A. Lovell, F. K. Chan, L. A. Cangahuala (Eds.), Advances in the Astronautical Sciences: AAS/AIAA Astrodynamics Conference, American Astronautical Society, Pittsburgh, PA, 2009, pp. 1–17.
  * Fathy et al. [2001] H. K. Fathy, J. A. Reyer, P. Y. Papalambros, A. G. Ulsoy,  On the coupling between the plant and controller optimization problems,  in: Proceedings of the 2001 American Control Conference, Institute of Electrical and Electronics Engineers, Arlington, VA, 2001, pp. 1864–1869. doi:10.1109/ACC.2001.946008.
  * Peters et al. [2009] D. L. Peters, P. Y. Papalambros, A. G. Ulsoy,  On measures of coupling between the artifact and controller optimal design problems,  in: International Design Engineering Technical Conferences, DETC2009-86868, American Society for Mechanical Engineers, San Diego, CA, 2009, pp. 1363–1372. doi:10.1115/DETC2009-86868.
  * Allison and Nazari [2010] J. T. Allison, S. Nazari,  Combined plant and controller design using decomposition-based design optimization and the minimum principle,  in: International Design Engineering Technical Conferences, DETC2010-28887, American Society for Mechanical Engineers, Montreal, Canada, 2010, pp. 765–774. doi:10.1115/DETC2010-28887.
  * Lee et al. [2019] Y. H. Lee, Vedant, R. H. Ewoldt, J. T. Allison,  Strain-actuated solar arrays for spacecraft attitude control assisted by viscoelastic damping,  in: Advances in Structural and Multidisciplinary Optimization, Proceedings of the 13th World Congress of Structural and Multidisciplinary Optimization (Beijing), Dalian University of Technology Electronic & Audio-Visual Press, Dalian, China, 2019, pp. 149–155.
  * Reyer and Papalambros [2000] J. A. Reyer, P. Y. Papalambros,  An investigation into modeling and solution strategies for optimal design and control,  in: Design Engineering Technical Conferences, DETC2000/DAC-14253, American Society for Mechanical Engineers, Baltimore, MD, 2000, pp. 283–292. doi:10.1115/DETC2000/DAC-14253.
  * Garcia-Sanz [2019] M. Garcia-Sanz,  Control co-design: An engineering game changer,  Advanced Control for Applications 1 (2019) e18. doi:10.1002/adc2.18.
  * Ashuri et al. [2014] T. Ashuri, M. B. Zaaijer, J. R. R. A. Martins, G. J. W. Van Bussel, G. A. M. Van Kuik,  Multidisciplinary design optimization of offshore wind turbines for minimum levelized cost of energy,  Renewable Energy 68 (2014) 893–905. doi:10.1016/j.renene.2014.02.045.
  * Du et al. [2020] X. Du, L. Burlion, O. Bilgen,  Control co-design for rotor blades of floating offshore wind turbines,  in: International Mechanical Engineering Congress and Exposition, IMECE2020-24605, American Society for Mechanical Engineers, 2020, p. V07AT07A052. doi:10.1115/IMECE2020-24605.
  * Deshmukh and Allison [2016] A. P. Deshmukh, J. T. Allison,  Multidisciplinary dynamic optimization of horizontal axis wind turbine design,  Structural and Multidisciplinary Optimization 53 (2016) 15–27. doi:10.1007/s00158-015-1308-y.
  * Jonkman et al. [2009] J. Jonkman, S. Butterfield, W. Musial, G. Scott, Definition of a 5-MW reference wind turbine for offshore system development, Technical Report NREL/TP-500-38060, National Renewable Energy Laboratory, Golden, CO, 2009.
  * Jonkman [2010] J. Jonkman, Definition of the floating system for phase IV of OC3, Technical Report NREL/TP-500-47535, National Renewable Energy Laboratory, Golden, CO, 2010.
  * Hansen et al. [2003] N. Hansen, S. D. M uller, P. Koumoutsakos,  Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES),  Evolutionary Computation 11 (2003) 1–18. doi:10.1162/106365603321828970.
  * Hansen [2016] N. Hansen,  The CMA evolution strategy: a tutorial,  arXiv CoRR abs/1604.00772 (2016). URL: http://{}arxiv.org{}/abs{}/1604.00772.
  * Wächter and Biegler [2006] A. Wächter, L. T. Biegler,  On the implementation of a primal-dual interior point filter line search algorithm for large-scale nonlinear programming,  Mathematical Programming 106 (2006) 25–57. doi:10.1007/s10107-004-0559-y.
  * Jonkman [2013] J. Jonkman,  The new modularization framework for the FAST wind turbine CAE tool,  in: 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, AIAA 2013-0202, American Institute of Aeronautics and Astronautics, American, Grapevine, TX, 2013, pp. 1–26. doi:10.2514/6.2013-202.
  * National Renewable Energy Laboratory [2021] National Renewable Energy Laboratory, Openfast, https://github.com/OpenFAST/openfast/releases/tag/v2.6.0, 2021. Online, accessed on June 7, 2021.
  * Deshmukh and Allison [2017] A. P. Deshmukh, J. T. Allison,  Design of dynamic systems using surrogate models of derivative functions,  ASME Journal of Mechanical Design 139 (2017) 101402–101402–12. doi:10.1115/1.4037407.
  * Lemmer [2018] F. Lemmer, Low-order modeling, controller design and optimization of floating offshore wind turbines, Ph.D. thesis, University of Stuttgart, Stuttgart, Germany, 2018. doi:10.18419/opus-10526.
  * Al-Solihat [2017] M. K. Al-Solihat, Dynamics modeling, simulation and analysis of a floating offshore wind turbine, Ph.D. thesis, McGill University, Montreal, Canada, 2017.
  * Hansen [2015] M. O. L. Hansen, Aerodynamics of Wind Turbines, 3 ed., Routledge, London, UK, 2015.
  * Gros [2013] S. Gros,  An economic NMPC formulation for wind turbine control,  in: 52nd IEEE Conference on Decision and Control, International Electrical and Electronics Engineers, Florence, Italy, 2013, pp. 1001–1006. doi:10.1109/CDC.2013.6760013.
  * Dormand and Prince [1980] J. R. Dormand, P. J. Prince,  A family of embedded runge-kutta formulae,  Journal of Computational and Applied Mathematics 6 (1980) 19–26. doi:10.1016/0771-050X(80)90013-3.
  * Hayman [2012] G. Hayman,  Mlife theory manual for version 1.00,  National Renewable Energy Laboratory, Golden, CO 74 (2012) 106.


*mps*

# Path Integral Method for Barrier Option Pricing Under Vasicek Model

Qi Chen 1<EMAIL_ADDRESS>Chao Guo 2<EMAIL_ADDRESS>1 School of
Economics and Management, Langfang Normal University, Langfang 065000, China
2 School of Science, Langfang Normal University, Langfang 065000, China

###### Abstract

Path integral method in quantum theory provides a new thinking for time
dependent option pricing. For barrier options, the option price changing
process is similar to the infinite high barrier scattering problem in quantum
mechanics; for double barrier options, the option price changing process is
analogous to a particle moving in a infinite square potential well. Using path
integral method, the expressions of pricing kernel and option price under
Vasicek stochastic interest rate model could be derived. Numerical results of
options price as functions of underlying prices are also shown.

##  I Introduction

In 1973, Black and Scholes derived a closed-form solution for European options
Black . Since then the field of derivative pricing has grown greatly Amin ;
Rubinstein . In resent years, path depend options have become increasingly
popular in financial markets, and barrier options are considered to be the
simplest types of path dependent options. Owing to the confinement of the
barrier, barrier options give investors more flexibility to transactions.
Snyder was the first to discuss down-and-out barrier options Snyder . Merton
derived a closed-form solution for the down-and-out barrier call option
Merton:1973 . Chiara et al priced barrier options by a Mellin transform method
Chiara . Karatzas and Wang gave closed-form expressions for the prices of
American up-and-out put options Karatzas . Baaquie et al discussed barrier
options and double barrier options by path integral method, and derived the
corresponding analytical expressions Baaquie . Kunitomo and Ikeda studied the
barrier options with floating barriers Kunitomo . Chen et al gave the integral
expressions for floating barrier option prices by path integral methods Chen .
All the papers mentioned above assumed that the model parameters such as
interest rate and volatility are constants. However, as we know, interest rate
and volatility change with time, and the Black-Scholes model should be revised
to be close to the reality. Lo et al give a closed-form expression for the
price of barrier options with time dependent parameters Lo . Lemmens et al
discussed Stochastic interest rate and volatility by means of path integral
method Lemmens . More about the discussions on option with time dependent
parameters could be studied in Hull ; Wiggins ; Roberts .

In this paper, we assume that the interest rate $r$ obeys the Vasicek model
Vasicek . Using $\Delta$-arbitrage strategy, the partial differential equation
(PDE) of option price $V(S,r,t)$ could be derived, which could be rewirtten
into a Schrödinger-type equation. Using path integral method, the pricing
kernel could be written out by the Hamiltonian of the option system Baaquie ;
Chen . The pricing kernel represents all the information for option price
evolution, which is in correspondence to the propagator in quantum theory. If
we can derive the expression of the pricing kernel, the expression of option
pricing could also be derived.

Our work is organized as follows. In Section 2, we will review the derivation
of the option price PDE under Vasicek model. In Section 3, we will derive the
integral expressions for the barrier option price by path integral method. In
Section 4, we will derive the integral expressions for the double barrier
option price by path integral method. Numerical results for option price as a
function of underlying asset price will be discussed in Section 5. We
summarize our main results in Section 6.

##  II Option Pricing Under the Vasicek Model

Assuming that the model of stock price $S_{t}$ obeys Brownian movement,

$\frac{{\rm d}S_{t}}{S_{t}}=r_{t}{\rm d}t+\sigma_{1}{\rm d}W_{t}^{1}$ (1)

where $\sigma_{1}$ is the volatility of $S_{t}$, and $r_{t}$ is the stochastic
interest rate which obeys Vasicek model Vasicek

${\rm d}r_{t}=a(\theta-r_{t}){\rm d}t+\sigma_{2}{\rm d}W_{t}^{2}$ (2)

where $\sigma_{2}$ is the volatility of $r_{t}$, $a$ and $\theta$ indicate
regression rate and long-term mean, respectively. $W_{t}^{1}$ and $W_{t}^{2}$
are standard Brownian movements, with the covariance

${\rm Cov}({\rm d}W_{t}^{1},{\rm d}W_{t}^{2})=\rho{\rm d}t,\ \ |\rho|\leq 1$
(3)

Using $\Delta$-arbitrage strategy, the option price $V(S,r,t)$ satisfies the
follow partial differential equation (PDE)

$\frac{\partial V}{\partial
t}+\frac{1}{2}\sigma_{1}^{2}S^{2}\frac{\partial^{2}V}{\partial
S^{2}}+\sigma_{1}\sigma_{2}\rho S\frac{\partial^{2}V}{\partial S\partial
r}+\frac{1}{2}\sigma_{2}^{2}\frac{\partial^{2}V}{\partial
r^{2}}+rS\frac{\partial V}{\partial S}+a(\theta-r)\frac{\partial V}{\partial
r}-rV=0$ (4)

for a call option, the final condition at $t=\tau$ is

$V(S,r,\tau)=(S-K)^{+}$ (5)

where $\tau$ is the maturity, and $K$ is the exercise price. By means of
variable substitutions

$y=\frac{S}{P(r,t;\tau)},\ \ \ \hat{V}(y,t)=\frac{V(S,r,t)}{P(r,t;\tau)}$ (6)

the PDE (4) could be simplified into

$\frac{\partial\hat{V}}{\partial
t}+\frac{1}{2}\hat{\sigma}^{2}(t)y^{2}\frac{\partial^{2}\hat{V}}{\partial
y^{2}}=0$ (7)

where $P(r,t;\tau)$ is the zero-coupon bond price. For Vasicek model,
$P(r,t;\tau)$ obeys the following PDE and the final condition

$\begin{split}&\frac{\partial P}{\partial
t}+\frac{\sigma^{2}}{2}\frac{\partial^{2}P}{\partial
r^{2}}+a(\theta-r)\frac{\partial P}{\partial r}-rP=0\\\ &P(r,T)=1\end{split}$
(8)

(LABEL:VasicekPDE) has the following explicit solution

$\begin{split}P(r,t;\tau)&=A(t)e^{-rB(t)}\\\ A(t)&={\rm
exp}\bigg{[}\frac{1}{a^{2}}[B^{2}(t)-(\tau-t)]\bigg{(}a^{2}\theta-\frac{\sigma^{2}}{2}-\frac{\sigma^{2}}{4a}B^{2}(t)\bigg{)}\bigg{]}\\\
B(t)&=\frac{1}{a}(1-e^{-a(\tau-t)})\end{split}$ (9)

and

$\hat{\sigma}(t)=\sqrt{\sigma_{1}^{2}+2\rho\sigma_{1}\sigma_{2}B(t)+\sigma_{2}^{2}B^{(}t)}$
(10)

the final condition for (7) is

$\hat{V}(y,\tau)=\frac{V(S,r,\tau)}{P(r,\tau;\tau)}=(y-K)^{+}$ (11)

Considering the variable substitution

$y=e^{x},\ \ \ -\infty<x<+\infty$ (12)

(7) could be changed into a Schrödinger-type equation

$\frac{\partial\hat{V}}{\partial t}=H\hat{V}$ (13)

with the Hamiltonian $H$ given by

$H=-\frac{1}{2}\hat{\sigma}^{2}(t)\frac{\partial^{2}}{\partial
x^{2}}+\frac{1}{2}\hat{\sigma}^{2}(t)\frac{\partial}{\partial x}$ (14)

which embodies the dynamic characteristic of the option system.

##  III Path Integral Method for Barrier Option Pricing

The barrier option Hamiltonian is

$H=-\frac{1}{2}\hat{\sigma}^{2}(t)\frac{\partial^{2}}{\partial
x^{2}}+\frac{1}{2}\hat{\sigma}^{2}(t)\frac{\partial}{\partial
x}=e^{\frac{1}{2}x}H_{\rm eff}e^{-\frac{1}{2}x}+U(x)$ (15)

where

$H_{\rm eff}=-\frac{1}{2}\hat{\sigma}^{2}(t)\frac{\partial^{2}}{\partial
x^{2}}+\frac{1}{8}\hat{\sigma}^{2}(t)$ (16)

and the potential $U(x)$ is

$U(x)=\left\\{\begin{aligned} 0&,&x<B,\\\ \infty&,&x\geq
B.\end{aligned}\right.$ (17)

The Schrödinger equation for $H_{\rm eff}$ is

$-\frac{1}{2}\hat{\sigma}^{2}(t)\frac{\partial^{2}\phi}{\partial\phi^{2}}+\frac{1}{8}\hat{\sigma}^{2}(t)\phi=E\phi,\
\ \ x<B$ (18)

where $\phi$ is the scattering state wave function of the option, which
describes the option price at time $t$, and $E$ is corresponding to the
“option energy”. The solution of (18) is

$\phi(x,t)=e^{ip(t)(x-B)}-e^{-ip(t)(x-B)},\ \ \ x<B$ (19)

where $B$ is the barrier level, and

$p^{2}(t)=\frac{2E}{\hat{\sigma}^{2}(t)}-\frac{1}{4}$ (20)

for $x\geq B$, the wave function $\phi(x,t)=0$.

Owing to $\hat{\sigma}(t)$ changes with time, we discretize $\tau$ so that
there are $N$ steps to maturity, with each time step $\epsilon=\tau/N$.
$\hat{\sigma}(t)$ tends to be a constant during each small enough time step.
The pricing kernel could be denoted as

$\begin{split}p_{\rm BO}(x,x^{\prime};\tau)&=\braket{x}{e^{-\tau
H}}{x^{\prime}}\\\ &=\lim_{\epsilon\to 0}\int_{-\infty}^{B}{\rm
d}x_{1}\int_{-\infty}^{B}{\rm d}x_{2}...\int_{-\infty}^{B}{\rm
d}x_{N-1}\braket{x}{e^{-\tau H}}{x_{1}}\braket{x_{1}}{e^{-\tau
H}}{x_{2}}...\braket{x_{N-1}}{e^{-\tau H}}{x^{\prime}}\end{split}$ (21)

where the completeness condition has been used. The $j$th matrix element is

$\begin{split}\braket{x_{j}}{e^{-\epsilon
H}}{x_{j+1}}=&e^{\alpha(x_{j}-x_{j+1})}e^{-\epsilon\gamma_{j}}\int_{0}^{+\infty}\frac{{\rm{d}}p_{j}}{2\pi}e^{-\frac{1}{2}\epsilon\sigma_{j}^{2}p_{j}^{2}}\big{[}e^{ip_{j}(x_{j}-B)}-e^{-ip_{j}(x_{j}-B)}\big{]}\big{[}e^{-ip_{j+1}(x_{j}-B)}-e^{ip_{j+1}(x_{j}-B)}\big{]}\\\
&=e^{\alpha(x_{j}-x_{j+1})}e^{-\epsilon\gamma_{j}}\int_{-\infty}^{+\infty}\frac{{\rm{d}}p_{j}}{2\pi}e^{-\frac{1}{2}\epsilon\sigma_{j}^{2}p_{j}^{2}}\big{[}e^{ip_{j}(x_{j}-x_{j+1})}-e^{ip_{j}(x_{j}+x_{j+1}-2B)}\big{]}\end{split}$
(22)

Using Gaussian integral formula, the integral in (22) could be calculated as

$\begin{split}&\int_{-\infty}^{+\infty}\frac{{\rm{d}}p_{j}}{2\pi}e^{-\frac{1}{2}\epsilon\sigma_{j}^{2}p_{j}^{2}}\big{[}e^{ip_{j}(x_{j}-x_{j+1})}-e^{ip_{j}(x_{j}+x_{j+1}-2B)}\big{]}\\\
=&\int_{-\infty}^{+\infty}\frac{{\rm{d}}p_{j}}{2\pi}e^{-\frac{1}{2}\epsilon\sigma_{j}^{2}\big{[}p_{j}-\frac{i(x_{j}-x_{j+1})}{\epsilon\sigma_{j}^{2}}\big{]}^{2}}e^{-\frac{(x_{j}-x_{j+1})^{2}}{2\epsilon\sigma_{j}^{2}}}-\int_{-\infty}^{+\infty}\frac{{\rm{d}}p_{j}}{2\pi}e^{-\frac{1}{2}\epsilon\sigma_{j}^{2}\big{[}p_{j}-\frac{i(x_{j}+x_{j+1}-2B)}{\epsilon\sigma_{j}^{2}}\big{]}^{2}}e^{-\frac{(x_{j}+x_{j+1}-2B)^{2}}{2\epsilon\sigma_{j}^{2}}}\\\
=&\frac{1}{\sqrt{2\pi\epsilon\sigma_{j}^{2}}}\bigg{[}e^{-\frac{(x_{j}-x_{j+1})^{2}}{2\epsilon\sigma_{j}^{2}}}-e^{-\frac{(x_{j}+x_{j+1}-2B)^{2}}{2\epsilon\sigma_{j}^{2}}}\bigg{]}\end{split}$
(23)

Similarly,

$\begin{split}\braket{x_{j-1}}{e^{-2\epsilon
H}}{x_{j+1}}=&\int_{-\infty}^{B}{\rm d}x_{j}\braket{x_{j-1}}{e^{-\epsilon
H}}{x_{j}}\braket{x_{j}}{e^{-\epsilon H}}{x_{j+1}}\\\ =&\int_{-\infty}^{B}{\rm
d}x_{j}e^{\frac{1}{2}(x_{j-1}-x_{j+1})}e^{-\epsilon(\gamma_{j-1}+\gamma_{j})}\frac{1}{\sqrt{2\pi\epsilon\sigma_{j-1}^{2}}}\frac{1}{\sqrt{2\pi\epsilon\sigma_{j}^{2}}}\times\\\
&\bigg{[}e^{-\frac{1}{2\epsilon\sigma_{j-1}^{2}}(x_{j-1}-x_{j})^{2}}-e^{-\frac{1}{2\epsilon\sigma_{j-1}^{2}}(x_{j-1}+x_{j}-2B)^{2}}\bigg{]}\bigg{[}e^{-\frac{1}{2\epsilon\sigma_{j}^{2}}(x_{j}-x_{j+1})^{2}}-e^{-\frac{1}{2\epsilon\sigma_{j}^{2}}(x_{j}+x_{j+1}-2B)^{2}}\bigg{]}\end{split}$
(24)

now we calculate the four integrals in (24)

$\begin{split}&\int_{-\infty}^{B}{\rm
d}x_{j}e^{-\frac{1}{2\epsilon\sigma_{j-1}^{2}}(x_{j-1}-x_{j})^{2}-\frac{1}{2\epsilon\sigma_{j}^{2}}(x_{j}-x_{j+1})^{2}}\\\
=&\int_{-\infty}^{B}{\rm
d}x_{j}e^{-\frac{\sigma_{j-1}^{2}+\sigma_{j}^{2}}{2\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}\big{[}x_{j}-\frac{x_{j-1}\sigma_{j}^{2}+x_{j+1}\sigma_{j-1}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}\big{]}^{2}}e^{-\frac{(x_{j-1}-x_{j+1})^{2}}{2\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\\\
=&\bigg{(}\frac{1}{2}\sqrt{\frac{2\pi\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}}+\sqrt{\frac{\pi\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}}\
{\rm
erf}\bigg{[}\frac{B\sigma_{j}^{2}+B\sigma_{j-1}^{2}-x_{j-1}\sigma_{j}^{2}-x_{j+1}\sigma_{j-1}^{2}}{\sqrt{2\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\bigg{]}\bigg{)}e^{-\frac{(x_{j-1}-x_{j+1})^{2}}{2\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\end{split}$
(25) $\begin{split}&\int_{-\infty}^{B}{\rm
d}x_{j}e^{-\frac{1}{2\epsilon\sigma_{j-1}^{2}}(x_{j-1}+x_{j}-2B)^{2}-\frac{1}{2\epsilon\sigma_{j}^{2}}(x_{j}+x_{j+1}-2B)^{2}}\\\
=&\int_{-\infty}^{B}{\rm
d}x_{j}e^{-\frac{\sigma_{j-1}^{2}+\sigma_{j}^{2}}{2\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}\big{[}x_{j}+\frac{(x_{j-1}-2B)\sigma_{j}^{2}+(x_{j+1}-2B)\sigma_{j-1}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}\big{]}^{2}}e^{-\frac{(x_{j-1}-x_{j+1})^{2}}{2\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\\\
=&\bigg{(}\frac{1}{2}\sqrt{\frac{2\pi\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}}+\sqrt{\frac{\pi\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}}\
{\rm
erf}\bigg{[}\frac{-B\sigma_{j}^{2}-B\sigma_{j-1}^{2}+x_{j-1}\sigma_{j}^{2}+x_{j+1}\sigma_{j-1}^{2}}{\sqrt{2\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\bigg{]}\bigg{)}e^{-\frac{(x_{j-1}-x_{j+1})^{2}}{2\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\end{split}$
(26) $\begin{split}&\int_{-\infty}^{B}{\rm
d}x_{j}e^{-\frac{1}{2\epsilon\sigma_{j-1}^{2}}(x_{j-1}-x_{j})^{2}-\frac{1}{2\epsilon\sigma_{j}^{2}}(x_{j}+x_{j+1}-2B)^{2}}\\\
=&\int_{-\infty}^{B}e^{-\frac{\sigma_{j-1}^{2}+\sigma_{j}^{2}}{2\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}\big{[}x_{j}-\frac{x_{j-1}\sigma_{j}^{2}-(x_{j+1}-2B)\sigma_{j-1}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}\big{]}^{2}}e^{-\frac{(x_{j-1}+x_{j+1}-2B)^{2}}{2\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\\\
=&\bigg{(}\frac{1}{2}\sqrt{\frac{2\pi\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}}+\sqrt{\frac{\pi\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}}\
{\rm
erf}\bigg{[}\frac{B\sigma_{j}^{2}-B\sigma_{j-1}^{2}-x_{j-1}\sigma_{j}^{2}+x_{j+1}\sigma_{j-1}^{2}}{\sqrt{2\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\bigg{]}\bigg{)}e^{-\frac{(x_{j-1}+x_{j+1}-2B)^{2}}{2\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\end{split}$
(27) $\begin{split}&\int_{-\infty}^{B}{\rm
d}x_{j}e^{-\frac{1}{2\epsilon\sigma_{j-1}^{2}}(x_{j-1}+x_{j}-2B)^{2}-\frac{1}{2\epsilon\sigma_{j}^{2}}(x_{j}-x_{j+1})^{2}}\\\
=&\int_{-\infty}^{B}e^{-\frac{\sigma_{j-1}^{2}+\sigma_{j}^{2}}{2\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}\big{[}x_{j}+\frac{(x_{j-1}-2B)\sigma_{j}^{2}-x_{j+1}\sigma_{j-1}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}\big{]}^{2}}e^{-\frac{(x_{j-1}+x_{j+1}-2B)^{2}}{2\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\\\
=&\bigg{(}\frac{1}{2}\sqrt{\frac{2\pi\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}}+\sqrt{\frac{\pi\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}}{\sigma_{j-1}^{2}+\sigma_{j}^{2}}}\
{\rm
erf}\bigg{[}\frac{-B\sigma_{j}^{2}+B\sigma_{j-1}^{2}+x_{j-1}\sigma_{j}^{2}-x_{j+1}\sigma_{j-1}^{2}}{\sqrt{2\epsilon\sigma_{j-1}^{2}\sigma_{j}^{2}(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\bigg{]}\bigg{)}e^{-\frac{(x_{j-1}+x_{j+1}-2B)^{2}}{2\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\end{split}$
(28)

where

${\rm erf}(x)=\int_{0}^{x}e^{-\eta^{2}}{\rm d}\eta$ (29)

is the error function. Considering the error function is an odd function, (24)
could be simplified into

$\braket{x_{j-1}}{e^{-2\epsilon
H}}{x_{j+1}}=\frac{e^{\frac{1}{2}(x_{j-1}-x_{j+1})}e^{-\frac{1}{8}\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}{\sqrt{2\pi\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\bigg{[}e^{-\frac{(x_{j-1}-x_{j+1})^{2}}{2\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}-e^{-\frac{(x_{j-1}+x_{j+1}-2B)^{2}}{2\epsilon(\sigma_{j-1}^{2}+\sigma_{j}^{2})}}\bigg{]}$
(30)

Repeating the above calculation, the pricing kernel is

$\braket{x}{e^{-\tau
H}}{x^{\prime}}=\frac{e^{\frac{1}{2}(x-x^{\prime})}e^{-\frac{1}{8}\epsilon\sum_{j=1}^{N}\sigma_{j}^{2}}}{\sqrt{2\pi\epsilon\sum_{j=1}^{N}\sigma_{j}^{2}}}\bigg{[}e^{-\frac{(x-x^{\prime})^{2}}{2\epsilon\sum_{j=1}^{N}\sigma_{j}^{2}}}-e^{-\frac{(x+x^{\prime}-2B)^{2}}{2\epsilon\sum_{j=1}^{N}\sigma_{j}^{2}}}\bigg{]}$
(31)

where

$\begin{split}\epsilon\sum_{j=1}^{N}\sigma_{j}^{2}&=\epsilon\sum_{j=1}^{N}[\sigma_{1}^{2}+2\rho\sigma_{1}\sigma_{2}B(t)+\sigma_{2}^{2}B(t)^{2}]\\\
&=(\sigma_{1}^{2}+\frac{2\rho}{a}\sigma_{1}\sigma_{2}+\frac{\sigma_{2}^{2}}{a^{2}})\tau-\frac{2\sigma_{2}}{a^{2}}(\rho\sigma_{1}+\frac{\sigma_{2}}{a})(1-e^{-a\tau})+\frac{\sigma_{2}^{2}}{2a^{3}}(1-e^{-2a\tau})\end{split}$
(32)

and the option price could be denoted as

$V(x,r;\tau)=P(r;t,\tau)\int_{{\rm ln}K}^{B}\braket{x}{e^{-\tau
H}}{x^{\prime}}(e^{x^{\prime}}-K)$ (33)

##  IV Path Integral Method for Double Barrier Option Pricing

The double barrier option Hamiltonian is

$\begin{split}H&=e^{\frac{1}{2}x}H_{\rm eff}e^{-\frac{1}{2}x}+U(x)\\\
&=e^{\frac{1}{2}x}\bigg{(}-\frac{1}{2}\hat{\sigma}^{2}(t)\frac{\partial^{2}}{\partial
x^{2}}+\frac{1}{8}\hat{\sigma}^{2}(t)\bigg{)}e^{-\frac{1}{2}x}+U(x)\end{split}$
(34)

the potential $U(x)$ is

$U(x)=\left\\{\begin{aligned} \infty&,&x\leq B_{1},\\\ 0&,&B_{1}<x<B_{2},\\\
\infty&,&x\geq B_{2}.\end{aligned}\right.$ (35)

where $B_{1}$ and $B_{2}$ are the lower and higher barrier levels,
respectively. The Schrödinger equation for $H_{\rm eff}$ is

$-\frac{1}{2}\hat{\sigma}^{2}(t)\frac{\partial^{2}\phi}{\partial\phi^{2}}+\frac{1}{8}\hat{\sigma}^{2}(t)\phi=E\phi,\
\ \ B_{1}<x<B_{2}$ (36)

with the eigenstate

$\phi_{n}(x)=\left\\{\begin{aligned}
\sqrt{\frac{n\pi}{B_{2}-B_{1}}}\sin&[p_{n}(x-B_{1})],&B_{1}<x<B_{2},\\\
0&,&x\leq B_{1},\ x\geq B_{2}.\end{aligned}\right.$ (37)

where

$p_{n}=\frac{n\pi}{B_{2}-B_{1}},\ \
E_{n}=\frac{1}{2}\hat{\sigma}^{2}(t)p_{n}^{2},\ \ n=1,2,3,...$ (38)

Owing to $\hat{\sigma}(t)$ changes with time, we also discretize $\tau$ to $N$
steps, with each step $\epsilon=\tau/N$. The $j$th matrix element is

$\begin{split}\braket{x_{j}}{e^{-\epsilon
H}}{x_{j+1}}&=e^{\frac{1}{2}(x_{j}-x_{j+1})}e^{-\frac{1}{8}\epsilon\sigma_{j}^{2}}\sum_{n=1}^{\infty}\braket{x_{j}}{e^{-\frac{1}{2}\epsilon\hat{\sigma}^{2}(t)\hat{p}^{2}}}{n}\braket{n}{x_{j+1}}\\\
&=e^{\frac{1}{2}(x_{j}-x_{j+1})}e^{-\frac{1}{8}\epsilon\sigma_{j}^{2}}\sum_{n=1}^{\infty}e^{-\frac{1}{2}\epsilon
p_{n}^{2}\sigma_{j}^{2}}\phi_{n}(x_{j})\phi_{n}(x_{j+1})\end{split}$ (39)

where $\hat{p}=-i\frac{\partial}{\partial x}$ is the momentum operator, and

$\braket{x_{j}}{n}=\phi_{n}(x_{j})=\sqrt{\frac{2}{B_{2}-B_{1}}}\sin{p_{n}(x_{j}-B_{1})}$
(40)

similarly,

$\begin{split}\braket{x_{j}}{e^{-2\epsilon
H}}{x_{j+2}}&=\int_{B_{1}}^{B_{2}}{\rm d}x_{j+1}\braket{x_{j}}{e^{-\epsilon
H}}{x_{j+1}}\braket{x_{j+1}}{e^{-\epsilon H}}{x_{j+2}}\\\
&=e^{\frac{1}{2}(x_{j}-x_{j+2})}e^{-\frac{1}{8}\epsilon(\sigma_{j}^{2}+\sigma_{j+1}^{2})}\sum_{n=1}^{\infty}e^{-\frac{1}{2}\epsilon
p_{n}^{2}(\sigma_{j}^{2}+\sigma_{j+1}^{2})}\phi_{n}(x_{j})\phi_{n}(x_{j+2})\end{split}$
(41)

where the orthonormalization condition

$\int_{B_{1}}^{B_{2}}{\rm d}x\
\phi_{n}(x)\phi_{n^{\prime}}(x)=\delta_{nn^{\prime}}=\left\\{\begin{aligned}
0&,&n\neq n^{\prime},\\\ 1&,&n=n^{\prime}.\end{aligned}\right.$ (42)

Repeating the above calculation, the pricing kernel for double barrier option
is

$\braket{x}{e^{-\tau
H}}{x^{\prime}}=e^{\frac{1}{2}(x-x^{\prime})}e^{-\frac{1}{8}\epsilon\sum_{j=1}^{N}\sigma_{j}^{2}}\sum_{n=1}^{\infty}e^{-\frac{1}{2}\epsilon
p_{n}^{2}\sum_{j=1}^{N}\sigma_{j}^{2}}\phi_{n}(x)\phi_{n}(x^{\prime})$ (43)

and the option price is

$V(x,r;\tau)=P(r;t,\tau)\int_{{\rm ln}K}^{B_{2}}\braket{x}{e^{-\tau
H}}{x^{\prime}}(e^{x^{\prime}}-K)$ (44)

##  V Numerical Results

In Fig. 1-Fig. 3, we show the up-and-out barrier call price (left) and the up-
and-out double barrier call price (right) as functions of underlying price.
Without loss of generality, we set $\sigma_{1}=\sigma_{2}=0.3$. In Fig. 1,
solid lines of different colors represent different price curves under
different regression rates $a$. It is shown that the option prices increase
with the increasing of $a$. In Fig. 2, solid lines of different colors
represent different price curves under different long-term mean $\theta$. It
is shown that the option prices decrease with the increasing of $\theta$. In
Fig. 3, dashed lines represent price curves in the case of $\rho<0$, chain
lines represent price curves in the case of $\rho>0$, and the solid line is
corresponding to the case of $\rho=0$.

Figure 1: Barrier call price (left) and double-barrier call price (right) as
functions of underlying price for different regression rates, respectively.
Parameters: $B_{1}={\rm ln}100=4.5$, $B=B_{2}={\rm ln}130=4.867$, $K=100$,
$\rho=0.5$, $\theta=0.04$, $r_{0}=0.05$, $\sigma_{1}=\sigma_{2}=0.3$,
$\tau=1$.

Figure 2: Barrier call price (left) and double-barrier call price (right) as
functions of underlying price for different long-term means, respectively.
Parameters: $B_{1}={\rm ln}100=4.5$, $B=B_{2}={\rm ln}130=4.867$, $K=100$,
$\rho=0.5$, $a=1$, $r_{0}=0.05$, $\sigma_{1}=\sigma_{2}=0.3$, $\tau=1$.

Figure 3: Barrier call price (left) and double-barrier call price (right) as
functions of underlying price for different correlation coefficients,
respectively. Parameters: $B_{1}={\rm ln}100=4.5$, $B=B_{2}={\rm
ln}130=4.867$, $K=100$, $a=1$, $\theta=0.04$, $r_{0}=0.05$,
$\sigma_{1}=\sigma_{2}=0.3$, $\tau=1$.

##  VI Conclusion

Barrier option price changing with time could be analogous to a particle
moving under some special potential. For stochastic interest rates, the
maturity should be split into $N$ steps, and each matrix element during a tiny
step is calculated by Gaussian integral. After $N$ times integral calculation,
the pricing kernel could be written in the form of series operation. Path
integral is an effective method linking option price changing to a particle
moving under some potential in the space. The pricing pf other barrier options
such as step options could be studied by defining appropriate potentials.

## References

  * (1) F. Black and M. Scholes, _The Pricing of Options and Corporate Liabilities_ , Journal of political economy 81, 637 (1973).
  * (2) K. Amin and A. J. Morton, _Implied Volatility Functions in Arbitrage-Free Term Structure Models_ , Journal of Financial Economics 35, 141 (1994).
  * (3) M. Rubinstein, _Displaced Diffusion Option Pricing_ , Journal of Finance 38, 213 (1983).
  * (4) G. L .Snyder, _Alternative Forms of Options_ , Financial Analysts Journal 25, 93-99 (1969).
  * (5) R. C. Merton, _Theory of Rational Option Pricing_ , The Rand Journal of Economics 4, 141-183 (1973).
  * (6) G. Chiara, M. R. Rodrigo and S. Simona _A Mellin Transform Approach to Barrier Option Pricing_ , Ima Journal of Management Mathematics 1 (2018).
  * (7) I. Karatzas and H. Wang _A Barrier Option of American Type_ , Applied Mathematics and Optimization 42: 259-279 (2000).
  * (8) B. E. Baaquie, C. Coriano and M. Srikant, _Hamiltonian and Potentials in Derivative Pricing Models: Exact Results and Lattice Simulations_ , Physica A 334, 531-57 (2004).
  * (9) N. Kunitomo and M. Ikeda, _Pricing Options With Curved Boundaries_ , Mathematical Finance 2(4) (1992).
  * (10) Q. Chen, H. T. Wang and C. Guo, _Path Integral Method for Floating Barrier Option Pricing_ , [arXiv:2209.12542].
  * (11) C. F. Lo, H. C. Lee and C. H. Hui _A Simple Approach for Pricing Barrier Options with Time-Dependent Parameters_ , Quantitative Finance 3(2), (2003), 98-107.
  * (12) D. Lemmens, M. Wouters, J. Tempere and S. Foulon _A Path Integral Approach to Closed-Form Option Pricing Formulas with Applications to Stochasitic Volatility and Interest Rate Models_ , Phys. Rev. E 78(1), (2008), 016101.
  * (13) J. C. Hull and W. Alan, _The Pricing of Options on Assets with Stochastic Volatilities_ , Journal of Finance 42, (1987),281-300.
  * (14) J. Wiggins, _Option Values under Stochastic Volatility: Theory and Empirical Estimates_ , Journal of Financial Economics 19, (1987),351-372.
  * (15) G. O. Roberts and C. F. Shortland _Pricing Barrier Options with Time-Dependent Coefficients_ , Mathematical Finance 7(1), (1997), 83-93.
  * (16) O. A. Vasicek _An Equilibrium Characterization of the Term Structure_ , Journal of Financial Economics 5: 177-88 (2004).