BUG: Last frequency component is wrong

[cmichelenstrofer]

Describe the bug
The last frequency component of the results is not continues but has a significant jump.

To Reproduce
Run a case with irregular waves where the spectrum has significant ammount of energy in the last frequency (e.g. cut off the spectrum early).

Expected behavior
Smooth tail.

Observed behavior
Large jump at last frequency.

System:
All.

Additional information
We suspect this is related to the DFT matrix.

[cmichelenstrofer]

@michaelcdevin can you add more details?

[michaelcdevin]

See attached zipped file for a basic example to replicate this behavior using WecOptTool v2.4.0 using a similar setup to Tutorial 1, part 1: issue225_test.zip

The highest frequency will give a different optimal excitation force each time as well. Since there is a large amount of energy in that frequency, this will substantially alter the optimal power:

Run 1:

Optimal power: -391.20254096338795 W

Run 2:

Optimal power: -398.47711760183483 W

[michaelcdevin]

At least part of the issue here is that the wot.td_to_fd should not be multiplying the DC and top frequency by 2, since the two-sided wave spectrum would not have Hermitian symmetry at these frequencies for purely real signals (since we're removing the sine component of the top frequency). This is clear when looking at the effort in the frequency domain when calculating power by the top shared frequency being half the magnitude (blue is using WecOptTool v2.1.0, blue is WecOptTool v2.4.0, both using the same hard-coded states except for the added last frequency for the v2.1.0 version):

Fixing this corrects the effort, but the top frequency in the solution is still wrong, so I'm still troubleshooting this.

[michaelcdevin]

The root causes for this behavior have been identified:

1. The truncations of the last row/columns from Remove sine component of 2-point wave #191 now means the last row of the derivative matrix is all zeros, making the evaluated velocity (i.e. from np.dot(derivative_mat, pos)) at the top frequency be zero no matter what, which should not be the case. This makes the last element of the WEC dynamics state vector have no effect on the solution. This causes:

• The last element of the dynamics state vector to be random, which is why the results at the top frequency are not smooth.
• The WEC inertia to be different (since the last element of the dynamics state vector is random), which is why the optimal power converges to a different solution. Note that the control state vector is never multiplied by the derivative matrix for unstructured controllers, so the PTO force will still result in a consistent value, which is why the solver still converges at all.

2. Not an explicit cause, but certainly bears mentioning: versions of WecOptTool prior to v2.2.0 (i.e. before Remove sine component of 2-point wave #191) had a bug with the time vector creation that was masking this issue. This was due to there being one too many components in the time series making the time intervals incorrect, thus the top frequency was not actually the two-point wave. See the example code below using WecOptTool v2.1.0. If the time vector is corrected in v2.1.0, the solution blows up since the last component of the control state vector has no effect on the PTO force (since the sine term at the last time step will always evaluate to zero) and nothing is constraining the velocity, so the last element of the WEC dynamic state vector exponentially increases at the last time step.

Example of bug from cause (2) -- use WecOptTool v2.1.0:

import numpy as np
import wecopttool as wot

nfreq = 10
f1 = 0.01
t = wot.time(f1, nfreq) # should be multiples of 5
w = wot.frequency(f1, nfreq, False)
wt = np.outer(t, w)
top_sin = np.sin(wt)[:, -1] # sine component of two-point wave, should be all zeros
print(f'Time vector: {t}')
print(f'Sine component of 2pt wave: {top_sin}')

# Time vector: [ 0.          4.76190476  9.52380952 14.28571429 19.04761905 23.80952381
# 28.57142857 33.33333333 38.0952381  42.85714286 47.61904762 52.38095238
# 57.14285714 61.9047619  66.66666667 71.42857143 76.19047619 80.95238095
# 85.71428571 90.47619048 95.23809524]
# Sine component of 2pt wave: [ 0.          0.14904227 -0.29475517  0.43388374 -0.56332006  0.68017274
# -0.78183148  0.8660254  -0.93087375  0.97492791 -0.9972038   0.9972038
# -0.97492791  0.93087375 -0.8660254   0.78183148 -0.68017274  0.56332006
# -0.43388374  0.29475517 -0.14904227]

Point regarding the above comment:
Given this new context, the correction of the DC and top frequencies not being multiplied by 2 is still correct in theory and should be fixed, but is not a reason for the top frequency being wrong.

[cmichelenstrofer]

Another bug, but not sure if it fixes it: The second derivative operator is not the 1st derivative squared because of the nyquist frequency. @michaelcdevin lets talk about this next week

[michaelcdevin]

After further discussion, we determined the velocity being zero at the top frequency is actually correct behavior. Since $\omega t = n \pi$ for the Nyquist frequency, $v(f_{Nyquist}) = x'(f_{Nyquist}) \rightarrow \frac{\mathrm{d} }{\mathrm{d} t}(cos(\omega t)) = -\omega sin(\omega t) = 0$. Therefore, for unstructured controllers, the last frequency being arbitrary is the expected response, since the objective function of average power only considers velocity, not position or acceleration.

However, something is still not right here. I am still getting the same behavior when using a PI or PID controller, where position/acceleration should affect the objective function since it impacts the force.

Related bugs mentioned in this issue regarding the incorrect DFT matrix behavior and the 2nd derivative matrix affecting the acceleration have been addressed in #232.

[cmichelenstrofer]

One clarification. The last component of velocity (last component in the state vectors is $cos(N \omega t)$ ) depends on the $sin(N\omega t)$ amplitude of the position which we do not include in our state (arbitrary since the $sin$ is always zero). So we need to make an assumption... and that assumption is that the amplitude of the sine component of position is zero (pure cosine).

[michaelcdevin]

A bit more info on the behavior with a structured controller.

I'm testing convergence with a consistent x_wec_0 and x_opt_0 with a PI controller and 16 frequencies. Rerunning the same problem several times, the optimal x_opt has its first two elements different each time, and x_wec is completely different each time. Except for the first two elements, x_opt will converge to the same values even if x_opt_0 and/or x_wec_0 are changed.

The post processed results are the same as above, with an arbitrary last frequency and a different optimal average power each time.

[michaelcdevin]

Disregard my earlier comment, I forgot to switch my script back to using the same waves each time...

The behavior I'm seeing with the PI controller w/ 16 frequencies is that x_wec and x_opt both converge to consistent values except for the last element of x_opt, which is unchanged from the last element in x_opt_0. Again, the PI controller means that the position of the Nyquist frequency should have an impact, but it isn't.

[cmichelenstrofer]

Through this issue we

• fixed the DFT matrix
• use a different second derivative operator

The last issue we were discussing is not a bug and is actually the expected behavior. See the following writeup showing how for the linear problem with average mechanical power as the objective and an unstructured controller the DC and Nyquist components of the WEC position and PTO force are arbitrary.

Power_last_frequency.pdf

Same as PDF above, but in Markdown, some formatting issues...

Objective function: Average mechanical power

The states are X_wec and F_pto which
are the WEC position and PTO force time series in the frequency domain,
using an unstructured controller for the force. The WEC velocity can be
obtained as the derivative of the position.

$$
X_{wec} = \begin{bmatrix} x_{dc} \\ x_{a1} \\ x_{b1} \\ x_{a2}\end{bmatrix} \qquad
\dot{X}_{wec} = \begin{bmatrix} 0 \\ x_{b1}\omega_1 \\ -x_{a1}\omega_1 \\ 0
\end{bmatrix} \qquad
\ddot{X}_{wec} = \begin{bmatrix} 0 \\ -x_{a1}\omega_1^2 \\ -x_{b1}\omega_1^2 \\
-x_{a2} \omega_2^2
\end{bmatrix} \qquad
F_{pto} = \begin{bmatrix} f_{dc} \\ f_{a1} \\ f_{b1} \\ f_{a2}\end{bmatrix}$$

Note that the second derivative state can not be obtained by applying
the same derivative matrix to the velocity as was applied to the
position to obtain velocity. That is, the second derivative operator
(matrix) is not the first derivative operator squared; these differ in
the treatment of the last frequency component.

In the time domain, these are:

$$\begin{split}
x_{wec}(t) &= x_{dc} + x_{a1}\cos(\omega_1 t) + x_{b1}\sin(\omega_1 t) + x_{a2}\cos(2\omega_1 t) \\
\dot{x}_{wec}(t) &= x_{b1}\omega_1\cos(\omega_1 t) + -x_{a1}\omega_1\sin(\omega_1 t) \\
f_{pto}(t) &= f_{dc} + f_{a1}\cos(\omega_1 t) + f_{b1}\sin(\omega_1 t) + f_{a2}\cos(2\omega_1 t)
\end{split}$$

The power in the time domain is the product of velocity and force. This
is equivalent to a convolution integral in the frequency domain.

$$
\begin{split}
p(t) = &\dot{x}_{wec}(t)\cdot f_{pto}(t) \\
= &\quad f_{dc}x_{b1}\omega_1\cos(\omega_1 t) \\
&+ f_{a1}\cos(\omega_1 t)x_{b1}\omega_1\cos(\omega_1 t)\\
&+ f_{b1}\sin(\omega_1 t)x_{b1}\omega_1\cos(\omega_1 t) \\
&+ f_{a2}\cos(2\omega_1 t)x_{b1}\omega_1\cos(\omega_1 t) \\
&- f_{dc}x_{a1}\omega_1\sin(\omega_1 t) \\
&- f_{a1}\cos(\omega_1 t)x_{a1}\omega_1\sin(\omega_1 t)\\
&- f_{b1}\sin(\omega_1 t)x_{a1}\omega_1\sin(\omega_1 t) \\
&- f_{a2}\cos(2\omega_1 t)x_{a1}\omega_1\sin(\omega_1 t) \\
\end{split}$$

The average power is the integral of power in the time domain, or the DC
component of power in the frequency domain. Using orthogonality of
sinusoids, many of these terms are zero. Similarly, the integral of a
sine or cosine over an integer number of periods is zero.

$$\begin{split}
\bar{P}_{mech} = &\int_0^{2\pi/\omega_1} p(t) \mathop{}\!\mathrm{d}{t} \\
= &\int_0^{2\pi/\omega_1} \biggl[ \\
% &\qquad (f_{dc}x_{b1}\omega_1)\cos(\omega_1 t)\\
&\quad+ (f_{a1}x_{b1}\omega_1)\cos(\omega_1 t)\cos(\omega_1 t) \\
% &\quad- (f_{dc}x_{a1}\omega_1)\sin(\omega_1 t) \\
&\quad- (f_{b1}x_{a1}\omega_1)\sin(\omega_1 t)\sin(\omega_1 t) \\
&\biggr]\mathop{}\!\mathrm{d}{t}\\
= &\pi(f_{a1}x_{b1} - f_{b1}x_{a1})
\end{split}$$

If we generalize to n frequencies:

$$
\bar{P}_{mech} = \pi\left( \sum_{k=1}^{n-1} k \left(f_{ak}x_{bk} - f_{bk}x_{ak} \right) \right)$$

Neither the DC component nor the last frequency component of the PTO
force or the WEC position have any effect on mechanical power.

Dynamics

The linear dynamics can be expressed in the frequency domain, where each
frequency component satisfies the equation of motion independently from
the other frequency components.

$$(A_i + m)\ddot{X}_i + B_i\dot{X}_i + KX_i = F_{e,i} + F_{pto,i}$$

where
X_i = x_a, i + j**x_b, i
is the complex amplitude of position, A_i and
B_i are the added mass and radiation damping coefficients,
and F_e, i is the complex amplitude of excitation force
for frequency ω_i = i * ω₁. The
hydrostatic stiffness K is constant. The complex amplitude of the PTO
force is given as
F_pto, i = f_a, i + j**f_b, i,
where j is the imaginary unit.

DC frequency

For the first (DC) frequency, there is no radiation or excitation.
K**x_d**c = f_d**c
If f_d**c is arbitrary, as is the case with mechanical
power as the objective function and linear dynamics, then the mean (DC)
position is related to this arbitrary DC force by the hydrostatic
stiffness. In practice the mean position can be driven to zero by adding
a max symmetrical PTO force constraint, or by considering electric
power with losses.

Nyquist frequency

For the last frequency n, the EOM is

$$
\begin{split}
-x_{a,n}\omega_n^2(A_n+m) +Kx_{a,n} &= f_{ex,a,n} + f_{pto,a,n}\\
x_{a,n} = \frac{f_{ex,a,n} + f_{pto,a,n}}{K-n^2\omega_1^2(A_n+m)}
\end{split}$$

If f_pto is arbitrary (as above, with
objective=average mechanical power), then the last frequency component
of the position is related to that arbitrary value by the relation
above.

Form observations, the Nyquist component can also be driven to zero by a
force constraint (and probably by an electric energy loss).

Intermediate frequencies

For intermediate frequencies, position, velocity, and acceleration
amplitudes are all non-zero. The velocity and acceleration complex
amplitudes can be written in terms of the position, as
Ẋ_i = j**ω_iX_i and
$\ddot{X}_i = -\omega_i^2 X_i$. The relation between the PTO force
and the position is given as

$$X_i = \frac{F_{e,i} + F_{pto,i}}{-\omega_i^2(A_i+m)+j\omega_iB+K}$$

In general all these intermediate frequency components of the PTO force
and WEC state affect the objective function and are therefore not
arbitrary.

Structured Controller

A structured controller changes how the the PTO force is computed but
the result from the first section
still shows no effect of the DC and Nyquist components on the objective
function. However, since PTO force is no longer a free variable, the
dynamic equations have a specific solution, (even at DC and Nyquist
component). For a PI controller,
equation [eq:dynamics_DC]
is satisfied if either K_p = K (the integral gain equals
the hydrostatic stiffness) which is likely suboptimal, or the DC
component of position is zero.
Equation [eq:dynamics_Ny]
becomes
$$x_{a,n} = \frac{f_{ex,a,n}}{K-n^2\omega_1^2(A_n+m)-K_p}$$

Similiatly, for a PID
equation [eq:dynamics_DC]
becomes

$$\begin{split}
Kx_{dc} &= K_px_{dc} - K_dn^2\omega_1^2x_{dc} \\
&= (K_p-n^2\omega_1^2K_d)x_{dc}
\end{split}$$

which implies either the relation in parenthesis is satisfied (likely
suboptimal) or x_d**c is zero.
Equation [eq:dynamics_Ny]
becomes

$$x_{a,n} = \frac{f_{ex,a,n}}{K-K_p-n^2\omega_1^2(A_n+m-K_d)}$$

Electrical Power

PTO force and velocity are converted to current and voltage using the
impedance matrix. The impedance matrix can capture electrical losses. To
compute the current and voltage we use the PTO impedance in transmission
matrix form, which contains complex Fourier coefficients,

$$
\begin{bmatrix}
I \\
V
\end{bmatrix}
=
\begin{bmatrix}
A & B \\
C & D\\
\end{bmatrix}
\begin{bmatrix}
\dot{X} \\
F_p
\end{bmatrix}$$

Considering that all components in
[eq:pto_transmission]
are complex the intermediate components of the real Fourier coefficients
become significantly longer and we won’t state them right now, since we
want to focus on DC and Nyquist components,

$$
I = \begin{bmatrix}
B_{dc}f_{dc} \\
\cdots \\
\cdots \\
B_{a2}f_{a2}
\end{bmatrix}
\quad
V = \begin{bmatrix}
D_{dc}f_{dc} \\
\cdots \\
\cdots \\
D_{a2}f_{a2}
\end{bmatrix}
\quad$$

However, to re-emphasize, the intermediate components do impact the
average electrical power. The key difference in the electrical states
[eq:elec_states]
vs. the mechanical states
[eq:mech_states]
is that the DC and Nyquist component of the flow variable (current) are
non-zero. We recall, that we manually create a the DC component of the
PTO transmission matrix for ω = 0, assuming that the user will define
the PTO impedance matrix with the same frequency vector that they used
for the BEM code, i.e. no zero component.

The electrical power in the time domain includes extra terms, when
compared to the mechanical power
(equation [eq:mechanicalpower]),
that depend on the DC and Nyquist components, because
I_dc and I_n, a are non-zero and they
also do not go to zero when integrated over the full period, including
the two below:

$$\begin{split}
p_e(t) = &i(t)\cdot v(t) \\
= &\quad\cdots \\
&+ I_{dc}V_{DC} + I_{n,a}V_{n,a}\cos^2(\omega_nt)
\end{split}$$

This results in the average electrical power depending on the DC and
Nyquist compononents of both voltage and current. These in turn depend
on the DC and Nyquist components of the PTO force. The dynamic equations
([eq:dynamics_DC]-[eq:dynamics_Ny])
relate the forces to the position. Since the DC and Nyquist components
of force are no longer arbitrary neither are the DC and Nyquist
components of position.

The DC and Nyquist components of the WEC position and PTO force have
an effect on the electrical average power and are therefor not
arbitrary.

Note: The current DC and Nyquist components I_dc and
I_n, a are only non-zero if the PTO transmission matrix
has a defined B, this generally requires a gyrating energy
transduction within the PTO, such as a generator.

Conclusion

For average mechanical power the DC and Nyquist components of WEC
position and PTO force are arbitrary (have no effect on the objective
function) but ar related to each other through
equations [eq:dynamics_DC]
and
[eq:dynamics_Ny]
(to satify the dynamics).

For average electrical power using an impedance matrix, the DC and
Nyquist components of the WEC position and PTO force have an effect on
the objective function and are therefor not arbitrary.

Appendix

Electrical states

$$\begin{aligned}
\begin{split}
A\dot{X} &= (A_a+{i\mkern 1mu}A_b)(\dot{X}_a + {i\mkern 1mu}\dot{X}_b) \\
&= A_a\dot{X}_a - A_b\dot{X}_b + {i\mkern 1mu}(A_a\dot{X}_b + A_b\dot{X}_c)
\end{split} \\
\begin{split}
BF_p &= (B_a+{i\mkern 1mu}B_b)(f_a + {i\mkern 1mu}f_b) \\
&= B_a f_a - B_bf_b + {i\mkern 1mu}(B_af_b + B_bf_c)
\end{split}\\
\begin{split}
C\dot{X} &= (C_a+{i\mkern 1mu}C_b)(\dot{X}_a + {i\mkern 1mu}\dot{X}_b) \\
&= C_a\dot{X}_a - C_b\dot{X}_b + {i\mkern 1mu}(C_a\dot{X}_b + C_b\dot{X}_c)
\end{split} \\
\begin{split}
DF_p &= (D_a+{i\mkern 1mu}D_b)(f_a + {i\mkern 1mu}f_b) \\
&= D_a f_a - D_bf_b + {i\mkern 1mu}(D_af_b + D_bf_c)
\end{split}\end{aligned}$$
$$\begin{aligned}
\begin{split}
I = \begin{bmatrix}
I_a \\
I_b
\end{bmatrix}
=
\begin{bmatrix}
A_a\dot{X}_a - A_b\dot{X}_b + B_a f_a - B_bf_b \\
A_a\dot{X}_b + A_b\dot{X}_c + B_af_b + B_bf_c
\end{bmatrix}
\end{split} \\
\begin{split}
V = \begin{bmatrix}
V_a \\
V_b
\end{bmatrix}
=
\begin{bmatrix}
C_a\dot{X}_a - C_b\dot{X}_b + D_a f_a - D_bf_b \\
C_a\dot{X}_b + C_b\dot{X}_c + D_a f_a - D_bf_b
\end{bmatrix}
\end{split} \end{aligned}$$

Objective function: Average mechanical power

The states are X_wec and F_pto which
are the WEC position and PTO force time series in the frequency domain,
using an unstructured controller for the force. The WEC velocity can be
obtained as the derivative of the position.

$$
X_{wec} = \begin{bmatrix} x_{dc} \\ x_{a1} \\ x_{b1} \\ x_{a2}\end{bmatrix} \qquad
\dot{X}_{wec} = \begin{bmatrix} 0 \\ x_{b1}\omega_1 \\ -x_{a1}\omega_1 \\ 0
\end{bmatrix} \qquad
\ddot{X}_{wec} = \begin{bmatrix} 0 \\ -x_{a1}\omega_1^2 \\ -x_{b1}\omega_1^2 \\
-x_{a2} \omega_2^2
\end{bmatrix} \qquad
F_{pto} = \begin{bmatrix} f_{dc} \\ f_{a1} \\ f_{b1} \\ f_{a2}\end{bmatrix}$$

Note that the second derivative state can not be obtained by applying
the same derivative matrix to the velocity as was applied to the
position to obtain velocity. That is, the second derivative operator
(matrix) is not the first derivative operator squared; these differ in
the treatment of the last frequency component.

In the time domain, these are:

$$\begin{split}
x_{wec}(t) &= x_{dc} + x_{a1}\cos(\omega_1 t) + x_{b1}\sin(\omega_1 t) + x_{a2}\cos(2\omega_1 t) \\
\dot{x}_{wec}(t) &= x_{b1}\omega_1\cos(\omega_1 t) + -x_{a1}\omega_1\sin(\omega_1 t) \\
f_{pto}(t) &= f_{dc} + f_{a1}\cos(\omega_1 t) + f_{b1}\sin(\omega_1 t) + f_{a2}\cos(2\omega_1 t)
\end{split}$$

The power in the time domain is the product of velocity and force. This
is equivalent to a convolution integral in the frequency domain.

$$
\begin{split}
p(t) = &\dot{x}_{wec}(t)\cdot f_{pto}(t) \\
= &\quad f_{dc}x_{b1}\omega_1\cos(\omega_1 t) \\
&+ f_{a1}\cos(\omega_1 t)x_{b1}\omega_1\cos(\omega_1 t)\\
&+ f_{b1}\sin(\omega_1 t)x_{b1}\omega_1\cos(\omega_1 t) \\
&+ f_{a2}\cos(2\omega_1 t)x_{b1}\omega_1\cos(\omega_1 t) \\
&- f_{dc}x_{a1}\omega_1\sin(\omega_1 t) \\
&- f_{a1}\cos(\omega_1 t)x_{a1}\omega_1\sin(\omega_1 t)\\
&- f_{b1}\sin(\omega_1 t)x_{a1}\omega_1\sin(\omega_1 t) \\
&- f_{a2}\cos(2\omega_1 t)x_{a1}\omega_1\sin(\omega_1 t) \\
\end{split}$$

The average power is the integral of power in the time domain, or the DC
component of power in the frequency domain. Using orthogonality of
sinusoids, many of these terms are zero. Similarly, the integral of a
sine or cosine over an integer number of periods is zero.

$$\begin{split}
\bar{P}_{mech} = &\int_0^{2\pi/\omega_1} p(t) \mathop{}\!\mathrm{d}{t} \\
= &\int_0^{2\pi/\omega_1} \biggl[ \\
% &\qquad (f_{dc}x_{b1}\omega_1)\cos(\omega_1 t)\\
&\quad+ (f_{a1}x_{b1}\omega_1)\cos(\omega_1 t)\cos(\omega_1 t) \\
% &\quad- (f_{dc}x_{a1}\omega_1)\sin(\omega_1 t) \\
&\quad- (f_{b1}x_{a1}\omega_1)\sin(\omega_1 t)\sin(\omega_1 t) \\
&\biggr]\mathop{}\!\mathrm{d}{t}\\
= &\pi(f_{a1}x_{b1} - f_{b1}x_{a1})
\end{split}$$

If we generalize to n frequencies:

$$
\bar{P}_{mech} = \pi\left( \sum_{k=1}^{n-1} k \left(f_{ak}x_{bk} - f_{bk}x_{ak} \right) \right)$$

Neither the DC component nor the last frequency component of the PTO
force or the WEC position have any effect on mechanical power.

Dynamics

The linear dynamics can be expressed in the frequency domain, where each
frequency component satisfies the equation of motion independently from
the other frequency components.

$$(A_i + m)\ddot{X}_i + B_i\dot{X}_i + KX_i = F_{e,i} + F_{pto,i}$$

where
X_i = x_a, i + j**x_b, i
is the complex amplitude of position, A_i and
B_i are the added mass and radiation damping coefficients,
and F_e, i is the complex amplitude of excitation force
for frequency ω_i = i * ω₁. The
hydrostatic stiffness K is constant. The complex amplitude of the PTO
force is given as
F_pto, i = f_a, i + j**f_b, i,
where j is the imaginary unit.

DC frequency

For the first (DC) frequency, there is no radiation or excitation.
K**x_d**c = f_d**c
If f_d**c is arbitrary, as is the case with mechanical
power as the objective function and linear dynamics, then the mean (DC)
position is related to this arbitrary DC force by the hydrostatic
stiffness. In practice the mean position can be driven to zero by adding
a max symmetrical PTO force constraint, or by considering electric
power with losses.

Nyquist frequency

For the last frequency n, the EOM is

$$
\begin{split}
-x_{a,n}\omega_n^2(A_n+m) +Kx_{a,n} &= f_{ex,a,n} + f_{pto,a,n}\\
x_{a,n} = \frac{f_{ex,a,n} + f_{pto,a,n}}{K-n^2\omega_1^2(A_n+m)}
\end{split}$$

If f_pto is arbitrary (as above, with
objective=average mechanical power), then the last frequency component
of the position is related to that arbitrary value by the relation
above.

Form observations, the Nyquist component can also be driven to zero by a
force constraint (and probably by an electric energy loss).

Intermediate frequencies

For intermediate frequencies, position, velocity, and acceleration
amplitudes are all non-zero. The velocity and acceleration complex
amplitudes can be written in terms of the position, as
Ẋ_i = j**ω_iX_i and
$\ddot{X}_i = -\omega_i^2 X_i$. The relation between the PTO force
and the position is given as

$$X_i = \frac{F_{e,i} + F_{pto,i}}{-\omega_i^2(A_i+m)+j\omega_iB+K}$$

In general all these intermediate frequency components of the PTO force
and WEC state affect the objective function and are therefore not
arbitrary.

Structured Controller

A structured controller changes how the the PTO force is computed but
the result from the first section
still shows no effect of the DC and Nyquist components on the objective
function. However, since PTO force is no longer a free variable, the
dynamic equations have a specific solution, (even at DC and Nyquist
component). For a PI controller,
equation [eq:dynamics_DC]
is satisfied if either K_p = K (the integral gain equals
the hydrostatic stiffness) which is likely suboptimal, or the DC
component of position is zero.
Equation [eq:dynamics_Ny]
becomes
$$x_{a,n} = \frac{f_{ex,a,n}}{K-n^2\omega_1^2(A_n+m)-K_p}$$

Similiatly, for a PID
equation [eq:dynamics_DC]
becomes

$$\begin{split}
Kx_{dc} &= K_px_{dc} - K_dn^2\omega_1^2x_{dc} \\
&= (K_p-n^2\omega_1^2K_d)x_{dc}
\end{split}$$

which implies either the relation in parenthesis is satisfied (likely
suboptimal) or x_d**c is zero.
Equation [eq:dynamics_Ny]
becomes

$$x_{a,n} = \frac{f_{ex,a,n}}{K-K_p-n^2\omega_1^2(A_n+m-K_d)}$$

Electrical Power

PTO force and velocity are converted to current and voltage using the
impedance matrix. The impedance matrix can capture electrical losses. To
compute the current and voltage we use the PTO impedance in transmission
matrix form, which contains complex Fourier coefficients,

$$
\begin{bmatrix}
I \\
V
\end{bmatrix}
=
\begin{bmatrix}
A & B \\
C & D\\
\end{bmatrix}
\begin{bmatrix}
\dot{X} \\
F_p
\end{bmatrix}$$

Considering that all components in
[eq:pto_transmission]
are complex the intermediate components of the real Fourier coefficients
become significantly longer and we won’t state them right now, since we
want to focus on DC and Nyquist components,

$$
I = \begin{bmatrix}
B_{dc}f_{dc} \\
\cdots \\
\cdots \\
B_{a2}f_{a2}
\end{bmatrix}
\quad
V = \begin{bmatrix}
D_{dc}f_{dc} \\
\cdots \\
\cdots \\
D_{a2}f_{a2}
\end{bmatrix}
\quad$$

However, to re-emphasize, the intermediate components do impact the
average electrical power. The key difference in the electrical states
[eq:elec_states]
vs. the mechanical states
[eq:mech_states]
is that the DC and Nyquist component of the flow variable (current) are
non-zero. We recall, that we manually create a the DC component of the
PTO transmission matrix for ω = 0, assuming that the user will define
the PTO impedance matrix with the same frequency vector that they used
for the BEM code, i.e. no zero component.

The electrical power in the time domain includes extra terms, when
compared to the mechanical power
(equation [eq:mechanicalpower]),
that depend on the DC and Nyquist components, because
I_dc and I_n, a are non-zero and they
also do not go to zero when integrated over the full period, including
the two below:

$$\begin{split}
p_e(t) = &i(t)\cdot v(t) \\
= &\quad\cdots \\
&+ I_{dc}V_{DC} + I_{n,a}V_{n,a}\cos^2(\omega_nt)
\end{split}$$

This results in the average electrical power depending on the DC and
Nyquist compononents of both voltage and current. These in turn depend
on the DC and Nyquist components of the PTO force. The dynamic equations
([eq:dynamics_DC]-[eq:dynamics_Ny])
relate the forces to the position. Since the DC and Nyquist components
of force are no longer arbitrary neither are the DC and Nyquist
components of position.

The DC and Nyquist components of the WEC position and PTO force have
an effect on the electrical average power and are therefor not
arbitrary.

Note: The current DC and Nyquist components I_dc and
I_n, a are only non-zero if the PTO transmission matrix
has a defined B, this generally requires a gyrating energy
transduction within the PTO, such as a generator.

Conclusion

For average mechanical power the DC and Nyquist components of WEC
position and PTO force are arbitrary (have no effect on the objective
function) but ar related to each other through
equations [eq:dynamics_DC]
and
[eq:dynamics_Ny]
(to satify the dynamics).

For average electrical power using an impedance matrix, the DC and
Nyquist components of the WEC position and PTO force have an effect on
the objective function and are therefor not arbitrary.

Appendix

Electrical states

$$\begin{aligned}
\begin{split}
A\dot{X} &= (A_a+{i\mkern 1mu}A_b)(\dot{X}_a + {i\mkern 1mu}\dot{X}_b) \\
&= A_a\dot{X}_a - A_b\dot{X}_b + {i\mkern 1mu}(A_a\dot{X}_b + A_b\dot{X}_c)
\end{split} \\
\begin{split}
BF_p &= (B_a+{i\mkern 1mu}B_b)(f_a + {i\mkern 1mu}f_b) \\
&= B_a f_a - B_bf_b + {i\mkern 1mu}(B_af_b + B_bf_c)
\end{split}\\
\begin{split}
C\dot{X} &= (C_a+{i\mkern 1mu}C_b)(\dot{X}_a + {i\mkern 1mu}\dot{X}_b) \\
&= C_a\dot{X}_a - C_b\dot{X}_b + {i\mkern 1mu}(C_a\dot{X}_b + C_b\dot{X}_c)
\end{split} \\
\begin{split}
DF_p &= (D_a+{i\mkern 1mu}D_b)(f_a + {i\mkern 1mu}f_b) \\
&= D_a f_a - D_bf_b + {i\mkern 1mu}(D_af_b + D_bf_c)
\end{split}\end{aligned}$$
$$\begin{aligned}
\begin{split}
I = \begin{bmatrix}
I_a \\
I_b
\end{bmatrix}
=
\begin{bmatrix}
A_a\dot{X}_a - A_b\dot{X}_b + B_a f_a - B_bf_b \\
A_a\dot{X}_b + A_b\dot{X}_c + B_af_b + B_bf_c
\end{bmatrix}
\end{split} \\
\begin{split}
V = \begin{bmatrix}
V_a \\
V_b
\end{bmatrix}
=
\begin{bmatrix}
C_a\dot{X}_a - C_b\dot{X}_b + D_a f_a - D_bf_b \\
C_a\dot{X}_b + C_b\dot{X}_c + D_a f_a - D_bf_b
\end{bmatrix}
\end{split} \end{aligned}$$