Wind as a Stochastic Process

Where the gradient-estimate noise comes from, derived from the inflow spectrum rather than measured

Aykut C. Satici

Robot Control Lab · Systems Engineering · UT Dallas

The number we are trying to compute

Extremum seeking perturbs the torque gain and correlates the measured power against the perturbation. One period of that correlation returns one number, \(\hat g\), an estimate of the slope of the objective. Because the wind is random, \(\hat g\) is random.

\[\sigma(\hat g) := \text{standard deviation of one period's gradient estimate}\]

\(\sigma(\hat g)\) is the quantity this deck computes. Three design choices turn on it: the probe amplitude, how long to average, and how widely the gain scatters about its optimum once the adaptation has settled. Each of those dependences is derived in the closing section, once the machinery to do so exists.

Every value quoted for \(\sigma(\hat g)\) so far has been measured: read off a simulation, or bounded from a published table. None was predicted from the turbulence.

Why it is not a one-line calculation

\[\hat g = \frac{1}{T}\int_0^T \ln P(t)\,\frac{2}{a}\sin\omega t\;dt\]

\(P\) depends on the wind, which is a random process, not a random number. So \(\hat g\) is an integral of a random process against a deterministic weight, and its variance depends on how the wind’s randomness is distributed in time. A gust lasting an hour and a gust lasting a second contribute quite differently to the same integral.

Three things are therefore needed, and none of them is optional:

  1. a way to describe how a random process distributes its variance over timescales: the spectral density;
  2. a model of which spectrum turbulence actually has: Kaimal;
  3. a rule for pushing a process through a linear integral: Wiener-Khinchin together with Fubini.

The route

Four boxes left to right with arrows between them. A: reduce the randomness, y minus its mean equals three epsilon, needs P, C_P and lambda. B: make it a linear functional, g-hat noise equals six over aT times Z, needs the estimator. C: variance of a linear functional, sigma squared Z from S and w-hat, needs Wiener-Khinchin. D: put in the measured frequency content of turbulence, a number. Below, a banner giving the goal.

Four links. A shows that although log-power depends on the wind three separate ways, at the optimum only one survives. B uses that to write the estimator’s noise as a single integral of one random process. C is the general rule for the variance of such an integral. D puts in the measured frequency content of atmospheric turbulence and produces a number. Each is stated, then checked against simulation at the end.

What we need from the turbine

A wind turbine below rated wind speed. Write \(\bar V\) for the mean wind speed, \(R\) for the rotor radius, \(\Omega\) for the rotor angular speed, \(\rho\) for air density and \(A = \pi R^2\) for the swept area. Define the tip-speed ratio and the power coefficient:

\[\lambda := \frac{R\Omega}{V}, \qquad P = \tfrac12\rho A V^3\,C_P(\lambda)\]

\(C_P\) is a property of the blades. For the rotor studied here it has a single maximum \(C_P^{\max}\) at \(\lambda = \lambda^\star\); where that curve comes from, and what is physics versus measurement in it, is Appendix A.0b. The controller’s only actuator is the torque gain \(u\) in the generator-torque law \(\tau_g = u\Omega^2\); setting \(\dot\Omega = 0\) makes the rotor rest where \(C_P(\lambda)/\lambda^3 = c\,u\) for a rotor constant \(c\) (derived in Appendix A.0), so \(u\) selects an equilibrium \(\lambda_{\text{eq}}(u)\).

Define the objective \(J(u) := \ln C_P\big(\lambda_{\text{eq}}(u)\big)\), and write \(J'\), \(J''\) for its first two derivatives in \(u\). The optimum \(u^\star\) is where \(J'(u^\star) = 0\), equivalently \(\lambda_{\text{eq}}(u^\star) = \lambda^\star\). Numerically, for the NREL 5 MW rotor at \(\bar V = 8\) m/s, \(J'' = -1.41\times10^{-13}\) in SI units.

The decomposition

Split the wind into a mean and a fluctuation about it:

\[V(t) = \bar V\big(1 + \varepsilon(t)\big), \qquad \varepsilon(t) := \frac{\tilde v(t)}{\bar V}, \qquad \mathrm{TI} := \frac{\sigma_{\tilde v}}{\bar V} = \sigma_\varepsilon\]

\(\tilde v\) is the fluctuating part in m/s, \(\varepsilon\) its dimensionless counterpart, and \(\mathrm{TI}\) the turbulence intensity, \(0.10\) in the case studied here. For a random variable \(x\), \(\sigma_x := \sqrt{\mathbb{E}[x^2] - (\mathbb{E}[x])^2}\) is its standard deviation. For a process \(x(t)\) it is evaluated at a fixed \(t\); the assumption that makes it independent of \(t\) is stated on Link A, executed. Since \(\bar V := \mathbb{E}[V]\), both \(\tilde v\) and \(\varepsilon\) have mean zero, so \(\sigma_{\tilde v}^2 = \mathbb{E}[\tilde v^2]\) and \(\sigma_\varepsilon^2 = \mathbb{E}[\varepsilon^2]\).

Remark. The standard form is additive, \(V = \bar V + \tilde v\), with the statistics written for \(\tilde v\) [IEC]. Dividing by the constant \(\bar V\) gives the form above: \(\varepsilon = \tilde v/\bar V\) and \(\sigma_\varepsilon = \sigma_{\tilde v}/\bar V = \mathrm{TI}\).

\(\varepsilon\) is small

\(\sigma_\varepsilon = 0.10\), so a series in \(\varepsilon\) converges quickly and second-order terms contribute about \(1\%\).

\(\varepsilon\) is dimensionless

Every statistic of \(\varepsilon\) is a pure number, \(\mathrm{TI}\) among them. From here on the wind enters only through \(\varepsilon\), and \(\bar V\) appears explicitly wherever it does.

What we need from the estimator

Extremum seeking perturbs the gain with a sinusoid of amplitude \(a\) and period \(T\). Time is cut into dither periods \([nT,(n+1)T)\), \(n = 0,1,2,\dots\); \(u_n\) is the gain about which the \(n\)-th period dithers, and the measured log-power is correlated against the sinusoid over exactly that period (written below with the period’s own clock, \(t \in [0,T]\)):

\[u(t) = u_n + a\sin\omega t, \qquad \omega := \frac{2\pi}{T}, \qquad f_1 := \frac{1}{T} = \frac{\omega}{2\pi}\]

\[y(t) := \ln P(t), \qquad \hat g_n := \frac{1}{T}\int_0^T y(t)\,\frac{2}{a}\sin\omega t\;dt\]

Why that integral is a gradient estimate, and why the factor is \(2/a\), is Appendix A.1.

\(\hat g_n\) is an estimate of \(J'(u_n)\). Its mean error is the finite-amplitude bias \(\tfrac{a^2}{8}J'''\), derived on The bias grows as \(a^2\); its scatter, which is what this deck computes, comes entirely from the wind. Published values for the NREL 5 MW case [CLR19]: \(T = 150\) s and \(a = 13.6\%\) of \(u^\star\).

The power spectral density

\(\varepsilon\) is strictly stationary (Link A, executed) with finite variance, hence wide-sense stationary: \(\mathbb{E}[\varepsilon(t)] = 0\) for all \(t\), by construction of \(\varepsilon\), and \(\mathbb{E}[\varepsilon(t)\varepsilon(s)]\) depends on \(t - s\) alone. Define the autocovariance and, as its Fourier transform, the power spectral density; \(\Gamma\) is even, since \(\Gamma(-\tau) = \mathbb{E}[\varepsilon(t)\varepsilon(t-\tau)] = \Gamma(\tau)\) by stationarity (shift \(t\) by \(\tau\)):

\[\Gamma(\tau) := \mathbb{E}\big[\varepsilon(t)\,\varepsilon(t+\tau)\big], \qquad S^{2s}(f) := \int_{-\infty}^{\infty}\! \Gamma(\tau)e^{-i2\pi f\tau}d\tau, \qquad \Gamma(\tau) = \int_{-\infty}^{\infty}\! S^{2s}(f)e^{i2\pi f\tau}df\]

The pair is the Wiener-Khinchin theorem: for \(\Gamma \in L^1(\mathbb{R})\) the transform \(S^{2s}\) exists, is real, even and nonnegative, and the inversion formula holds. Setting \(\tau = 0\) in the inversion gives \(\Gamma(0) = \sigma_\varepsilon^2 = \int S^{2s}df\): the variance is the integral of \(S^{2s}\), which is why it is called a density. How the variance is weighted across frequency by a given integral is exactly what Link C computes, with \(|\hat w|^2\) as the weight.

Two conventions, and a factor of two

\(\Gamma\) is real and even, so \(S^{2s}\) is real and even. Tables therefore fold the negative frequencies onto the positive ones and quote the one-sided density:

\[S(f) := 2\,S^{2s}(f)\quad (f>0), \qquad\text{so}\qquad \sigma_\varepsilon^2 = \int_{-\infty}^{\infty}\!S^{2s}(f)\,df = \int_0^{\infty}\! S(f)\,df\]

Everything quoted from a standard is one-sided; every integral over all of \(\mathbb{R}\) below needs two-sided. The conversion \(S^{2s} = \tfrac12 S\) is applied once, explicitly, at the point of use.

Units follow from the definition: \(\varepsilon\) is dimensionless, so \([\Gamma] = 1\), \([S] = \text{Hz}^{-1}\), and \(\int S\,df\) is dimensionless as it must be.

How \(\varepsilon\) is synthesized, and its distribution

The validation section needs sample paths of \(\varepsilon\) with a prescribed one-sided density \(S\). For a record of length \(D\) sampled every \(\Delta t\):

  • Definition. \(N := D/(2\Delta t)\) bins, \(f_k := k/D\), \(\Delta f := 1/D\), and \[\varepsilon(t) := \sum_{k=1}^{N} A_k\cos(2\pi f_k t + \phi_k), \qquad A_k := \sqrt{2\,S(f_k)\,\Delta f}, \qquad \phi_k \ \text{independent, uniform on } [0,2\pi)\] In the runs here \(D = 36\,600\) s, \(\Delta t = 0.5\) s, \(N = 36\,600\).
  • Mean. \(\mathbb{E}\cos(2\pi f_k t + \phi_k) = 0\) for uniform \(\phi_k\), so \(\mathbb{E}[\varepsilon(t)] = 0\).
  • Autocovariance. Independent phases remove every cross term, so \[\Gamma(\tau) = \sum_{k=1}^{N} \tfrac12 A_k^2\cos(2\pi f_k\tau) = \sum_{k=1}^{N} S(f_k)\cos(2\pi f_k\tau)\,\Delta f \;\xrightarrow{\ \Delta f\to0\ }\; \int_0^\infty S(f)\cos(2\pi f\tau)\,df\] the inversion on The power spectral density for real even \(\Gamma\). At \(\tau = 0\), \(\operatorname{Var}\varepsilon = \sum_k S(f_k)\Delta f \to \int_0^\infty S\,df = \sigma_\varepsilon^2\).

How \(\varepsilon\) is synthesized, continued

  • Strict stationarity. Shifting \(t\) by \(\tau\) shifts \(\phi_k\) to \(\phi_k + 2\pi f_k\tau\); a uniform phase shifted modulo \(2\pi\) is uniform, so the joint law of \((\phi_1,\dots,\phi_N)\), and with it of the whole path, is shift invariant. This is assumption 2 of Link A, standing assumptions.
  • Distribution at fixed \(t\). \(\varepsilon(t)\) is a sum of \(N\) independent terms, each bounded by \(A_k\) with variance \(\tfrac12 A_k^2\). As \(N \to \infty\) with \(\max_k A_k \to 0\), the central limit theorem gives \(\varepsilon(t) \to \mathcal{N}(0,\sigma_\varepsilon^2)\). This matches IEC 61400-1, whose normal turbulence model takes the fluctuation at a fixed time to be Gaussian [IEC].

Every result in this deck uses \(\varepsilon\) through \(\mathbb{E}[\varepsilon] = 0\), \(\Gamma\), and \(\mathbb{E}[\varepsilon^4] < \infty\) only, so it holds for the synthesized process and for any other process with the same \(\Gamma\).

A shortcut that does not apply here

The squared magnitude of the window transform against frequency in units of the dither frequency, overlaid on the Kaimal spectrum normalized at the dither frequency. The window's main lobe runs from zero to twice the dither frequency; across that range the spectrum falls by a factor of twelve.

If \(S^{2s}\) were constant across the support of \(|\hat w|^2\) it could be taken outside, and Parseval would finish the job:

\[\operatorname{Var}(Z) \;\approx\; S^{2s}(f_1)\int_{\mathbb{R}}|\hat w|^2df = S^{2s}(f_1)\,\frac{T}{2} \;=\; S_\varepsilon(f_1)\,\frac{T}{4} \qquad\text{(valid iff } S \text{ is flat over the lobe)}\]

The condition fails here. \(w\) is one period of a sine, so \(|\hat w|^2\) has its main lobe on \([0, 2f_1]\), a \(100\%\) relative bandwidth, and across that interval \(S_\varepsilon\) falls by a factor of \(S_\varepsilon(0)/S_\varepsilon(2f_1) = 11.8\) at \(T = 150\) s.

The integral length scale

One more definition before the model, because it is the only free parameter in it. The normalized autocovariance and the integral time and length scales are

\[\varrho(\tau) := \frac{\Gamma(\tau)}{\Gamma(0)}, \qquad \mathcal{T} := \int_0^{\infty}\!\varrho(\tau)\,d\tau, \qquad L := \bar V\,\mathcal{T}\]

\(\mathcal{T}\) is the lag over which the signal retains memory of itself, and \(L = \bar V\mathcal{T}\) converts it to a length by Taylor’s frozen-turbulence hypothesis: eddies of size \(L\) are carried past the rotor at speed \(\bar V\).

\(L\) is fixed by the spectrum at zero frequency. Since \(\Gamma\) is real and even, \(e^{-i2\pi f\tau} = \cos 2\pi f\tau - i\sin 2\pi f\tau\) contributes only its cosine, and the two half-lines contribute equally: \(S(f) = 2S^{2s}(f) = 2\int_{\mathbb{R}}\Gamma(\tau)\cos(2\pi f\tau)\,d\tau = 4\int_0^\infty\Gamma(\tau)\cos(2\pi f\tau)\,d\tau\). Hence \[S(0) = 4\int_0^{\infty}\! \Gamma(\tau)\,d\tau = 4\sigma_\varepsilon^2\mathcal{T} = \frac{4L}{\bar V}\,\sigma_\varepsilon^2\] so any spectrum with integral scale \(L\) must have that value at \(f = 0\). That is where the prefactor in the next slide comes from.

The Kaimal spectrum

A specific one-parameter model for \(S\), fitted by Kaimal, Wyngaard, Izumi and Coté [KWIC72] to the 1968 Kansas boundary-layer measurements and adopted by IEC 61400-1 [IEC] as its normal turbulence model. In terms of \(\varepsilon\):

\[\boxed{\;S_\varepsilon(f) = \mathrm{TI}^2\;\frac{4L/\bar V}{\big(1+6fL/\bar V\big)^{5/3}}\;} \qquad f \ge 0,\ \text{one-sided}\]

Every piece is forced or measured, none is free:

\(4L/\bar V\)

The value at \(f = 0\) that the previous slide requires of any spectrum with integral scale \(L\).

exponent \(5/3\)

Kolmogorov’s inertial-subrange law, \(S\propto f^{-5/3}\), which the data follow above the corner.

the \(6\)

The one empirical number: it places the corner, at \(f = \bar V/6L\).

It is correctly normalized: substituting \(x = 6fL/\bar V\), \(\int_0^{\infty}\! S_\varepsilon df = \tfrac{2}{3}\int_0^\infty (1+x)^{-5/3}dx = \tfrac23\cdot\tfrac32\,\mathrm{TI}^2 = \mathrm{TI}^2\), as the definition of \(\mathrm{TI}\) demands.

What it looks like, and where the dither sits

Two panels. Left, the Kaimal one-sided density on log-log axes: flat below a corner frequency, then falling as f to the minus five thirds, with the 150 and 600 second dither frequencies marked on it. Right, f times S against log frequency, a single broad peak at U over 4L, showing where the variance actually accumulates.

[IEC] fixes \(L\) from the hub height \(z\) alone: \(L = 8.1\Lambda_1\) with \(\Lambda_1 = 0.7\min(z,60)\) m. For the NREL 5 MW rotor, \(z = 90\) m gives \(\Lambda_1 = 42\) m and \(L = 340\) m.

At \(\bar V = 8\) m/s that puts the corner at \(\bar V/6L = 3.9\times10^{-3}\) Hz. The right panel plots \(fS_\varepsilon\) against \(\ln f\), whose area is the variance (\(S\,df = fS\,d\ln f\)); it peaks at \(\bar V/4L = 5.9\times10^{-3}\) Hz (both shown in A.7). The published dither, \(f_1 = 1/150\,\text{s} = 6.7\times10^{-3}\) Hz, sits directly on that peak. The last part of this deck returns to that.

The answer, and what the shortcut would have cost

Evaluating the exact integral numerically against the shortcut, and both against a direct simulation of the turbine under a synthesized Kaimal inflow:

\(T\) exact \(\operatorname{Var}(Z)\) shortcut \(\sigma(\hat g)\) exact shortcut simulated exact/sim
\(150\) s \(15.65\) \(12.18\) \(5.21{\times}10^{-7}\) \(4.59{\times}10^{-7}\) \(5.48{\times}10^{-7}\) \(0.95\)
\(300\) s \(53.08\) \(45.73\) \(4.80{\times}10^{-7}\) \(4.45{\times}10^{-7}\) \(5.09{\times}10^{-7}\) \(0.94\)
\(600\) s \(152.9\) \(141.3\) \(4.07{\times}10^{-7}\) \(3.91{\times}10^{-7}\) \(3.83{\times}10^{-7}\) \(1.06\)

\[\boxed{\;\sigma(\hat g) = \frac{6}{aT}\sqrt{\int_{\mathbb{R}}\tfrac12 S_\varepsilon(f)\,|\hat w(f)|^2\,df}\;}\] Agreement with simulation is within \(\mathbf{6\%}\) across a fourfold range of \(T\), with no fitted parameter. The shortcut under-predicts by \(11\) to \(13\%\), because it samples \(S\) at \(f_1\) while the window also collects the larger \(S\) below \(f_1\).

Does it hold?

Each link checked separately against a simulation of the full nonlinear turbine, so a failure can be attributed to one link rather than to the chain.

How the check is set up

The turbine is integrated with an adaptive solver under a synthesized Kaimal inflow (How \(\varepsilon\) is synthesized), \(240\) periods of \(T = 150\) s at \(\bar V = 8\) m/s, realized \(\mathrm{TI} = 9.86\%\). The gain is held at \(u^\star\), with no adaptation between periods, so what is measured is the estimator alone.

testing A

Record \(y(t)\) and the \(\varepsilon(t)\) that drove it. Regress one on the other: the slope should be \(3\).

testing C

Form \(Z = \int\varepsilon w\,dt\) each period, take its sample variance, compare with \(\int S^{2s}|\hat w|^2 df\).

testing D

Take the sample standard deviation of \(\hat g\) itself and compare with \(\tfrac{6}{aT}\sqrt{\operatorname{Var}Z}\).

Link B is algebra given A: \(\hat g_{\text{noise}} = \tfrac{6}{aT}Z\) follows from \(y = y_0 + 3\varepsilon + r\) by linearity of the estimator, so it is tested through D, which compares the scatter of \(\hat g\) itself with the prediction.

What it decides

With \(\sigma(\hat g)\) in hand, the three design dependences can be derived. Each uses only what the four links established.

Bias and standard deviation: the definitions

Link A gave \(y - \bar y = 3\varepsilon + O(\varepsilon^2)\) at the optimum, and Link B gave \(\hat g_{\text{noise}} = \tfrac{6}{aT}Z\). So the estimator splits into a part that would be present with no wind and a part the wind adds:

\[\hat g = \underbrace{\frac{2}{aT}\int_0^T J\big(u + a\sin\omega t\big)\sin\omega t\,dt}_{\hat g_{\text{det}},\ \text{deterministic}} \;+\; \underbrace{\frac{6}{aT}Z}_{\hat g_{\text{noise}},\ \text{random}}\]

\[\operatorname{Bias}(\hat g) := \mathbb{E}[\hat g] - J'(u), \qquad \sigma(\hat g) := \sqrt{\operatorname{Var}(\hat g)} = \sqrt{\mathbb{E}\big[(\hat g - \mathbb{E}\hat g)^2\big]}\]

\(\hat g_{\text{det}}\) is a number, and \(\mathbb{E}[Z] = \int_0^T\mathbb{E}[\varepsilon]\,w\,dt = 0\) because \(\varepsilon\) is zero-mean. So the two quantities separate: \[\operatorname{Bias}(\hat g) = \hat g_{\text{det}} - J'(u), \qquad \operatorname{Var}(\hat g) = \Big(\frac{6}{aT}\Big)^2\operatorname{Var}(Z)\] The bias comes only from the probe; the scatter comes only from the wind.

The bias grows as \(a^2\)

Assume \(J\) is four times differentiable near \(u\) and \(a\) is small enough for the Taylor series to be used. Write \(s := \sin\omega t\) and \(\langle\cdot\rangle := \tfrac1T\int_0^T(\cdot)\,dt\). The moments of a sine over a whole period are

\[\langle s\rangle = 0,\qquad \langle s^2\rangle = \tfrac12,\qquad \langle s^3\rangle = 0,\qquad \langle s^4\rangle = \tfrac38,\qquad \langle s^5\rangle = 0\]

Expand \(J(u + as) = J + J'as + \tfrac{J''}{2}a^2s^2 + \tfrac{J'''}{6}a^3s^3 + O(a^4)\) and multiply by \(\tfrac{2}{a}s\) before averaging:

\[\hat g_{\text{det}} = \frac2a\Big[J\langle s\rangle + J'a\langle s^2\rangle + \tfrac{J''}{2}a^2\langle s^3\rangle + \tfrac{J'''}{6}a^3\langle s^4\rangle\Big] + O(a^4) = J' + \frac{a^2}{8}J''' + O(a^4)\]

\[\operatorname{Bias}(\hat g) = \frac{a^2}{8}\,J'''(u) + O(a^4)\] Every even power of \(a\) is killed by an odd moment of \(s\), which is why the leading term is \(a^2\) and not \(a\). Checked numerically on a cubic \(J\), where \(J'''\) is exact: agreement to six digits at \(a = 0.02\), \(0.05\), \(0.10\).

The \(\tfrac18\) is specific to a sinusoidal dither, through \(\langle s^4\rangle = \tfrac38\). A square-wave dither gives \(\tfrac16\).

The standard deviation falls as \(1/a\)

From the separation, \(\operatorname{Var}(\hat g) = (6/aT)^2\operatorname{Var}(Z)\) with \(Z = \int_0^T\varepsilon\sin\omega t\,dt\).

\(Z\) is built from the wind fluctuation and the fixed weight \(\sin\omega t\) on \([0,T]\). The amplitude \(a\) does not appear in it. Therefore \[\sigma(\hat g) = \frac{6}{aT}\sqrt{\operatorname{Var}(Z)} \;\propto\; \frac1a\] exactly for \(\hat g_{\text{noise}}\). The full scatter also carries \(\hat g_r\) from the remainder \(r\) of Link A, executed, concluded, with \(\operatorname{sd}(r - \mathbb{E}r) = O(\sigma_\varepsilon^2)\); the validation measures both together.

Why \(1/a\): the wind’s contribution to log-power is \(3\varepsilon(t)\) whatever the gain does, and the estimator divides by \(a\) to turn a power change into a slope. A fixed disturbance divided by a smaller number is a larger result.

The two scalings pull opposite ways. Bias wants \(a\) small; scatter wants \(a\) large. Balancing them needs the constant in the second relation, which is what Links C and D supplied.

Averaging over \(N\) periods

Average \(N\) consecutive estimates, \(\bar{\hat g}_N := \tfrac1N\sum_{n=1}^N \hat g_n\). Each has the same variance \(\sigma^2 := \operatorname{Var}(\hat g)\); if they were independent, \(\operatorname{Var}(\bar{\hat g}_N) = \sigma^2/N\). They are not quite, because the same wind spans adjacent periods. The exact statement uses the correlation between estimates \(k\) periods apart, which Link C also gives with the window shifted by \(kT\) (derived in A.8), and the variance of a mean of correlated samples (derived in A.9):

\[\rho_k := \frac{\operatorname{Cov}(Z_0,Z_k)}{\operatorname{Var}(Z)} = \frac{\int S^{2s}(f)\,|\hat w(f)|^2\cos(2\pi fkT)\,df}{\int S^{2s}(f)\,|\hat w(f)|^2\,df}, \qquad \operatorname{Var}(\bar{\hat g}_N) = \frac{\sigma^2}{N}\Big[1 + 2\sum_{k=1}^{N-1}\big(1-\tfrac kN\big)\rho_k\Big]\]

\(T\) \(\rho_1\) \(\rho_2\) sd reduction at \(N=3\) at \(N=8\) at \(N=24\) if independent
\(150\) s \(-0.099\) \(-0.008\) \(1.87\times\) \(3.14\times\) \(5.51\times\) \(1.73\), \(2.83\), \(4.90\)
\(600\) s \(-0.036\) \(-0.001\) \(1.78\times\) \(2.93\times\) \(5.09\times\) same

Adjacent estimates are slightly negatively correlated, because the half-window that ends one period begins the next with the opposite sign, so averaging does a little better than \(\sqrt N\). Measured lag-1 correlation of \(\hat g\) over \(300\) simulated periods at \(T = 600\) s: \(-0.024\), against the predicted \(-0.036\).

The terminal error of the closed loop

Close the loop with a constant gain \(\kappa > 0\), \(u_{n+1} = u_n + \kappa\hat g_n\) (a step up the gradient). Near the optimum \(J'(u_n) = J''\,\tilde u_n + O(\tilde u_n^2)\) with \(\tilde u_n := u_n - u^\star\), and \(\hat g_n = J'(u_n) + \tfrac{a^2}{8}J''' + \hat g_{\text{noise},n}\) by the bias and noise results. The bias is a constant: it shifts the fixed point to \(\tilde u = -a^2J'''/(8J'')\) and leaves the variance unchanged, so it is dropped here. Then, since \(J'' < 0\):

\[\tilde u_{n+1} = (1-K)\,\tilde u_n + \kappa\,\hat g_{\text{noise},n}, \qquad K := \kappa|J''|\]

a first-order autoregression. Its stationary variance, taking the noise uncorrelated across periods (the \(\rho_k\) above are small):

\[\operatorname{Var}(\tilde u) = \frac{\kappa^2\sigma^2(\hat g)}{1-(1-K)^2} = \frac{K}{2-K}\cdot\frac{\sigma^2(\hat g)}{J''^2} \qquad\Longrightarrow\qquad \sigma(\tilde u) = \sqrt{\frac{K}{2-K}}\;\frac{\sigma(\hat g)}{|J''|}\]

This is the sense in which the terminal error is “\(\sigma(\hat g)\) divided by the curvature”: the ratio \(\sigma(\hat g)/|J''|\), scaled by a factor that depends only on the loop gain \(K\). It is a stationary distribution: running longer leaves its width unchanged.

Checked in closed loop at \(K = 0.015\), \(T = 600\) s: predicted \(\sigma(\tilde u) = 0.106\,u^\star\), measured \(0.082\), ratio \(0.777\). That shortfall is the record, not the formula: with \(\alpha = 1-K = 0.985\) the loop’s correlation time is \(67\) periods, the run’s \(270\) post-transient samples span only four of them, and the sample variance of such a series is biased low by an expected factor of \(0.79\) in standard deviation (Appendix A.6). Prediction and measurement agree once that is accounted for.

Consequences

One design rule, and one discrepancy the derivation does not remove.

The dither period has a worst case

Gradient noise and second-harmonic curvature bias against dither period on a logarithmic axis. Noise rises to a broad maximum near 110 seconds and falls slowly for longer periods; the curvature bias rises from minus fifty per cent toward zero as the period lengthens. The published 150 second period sits one per cent below the noise maximum.

For the trend only, take \(\operatorname{Var}(Z) \approx S_\varepsilon(f_1)\,T/4\) with \(f_1 = 1/T\). Above the corner (short \(T\)), \(S_\varepsilon \propto f^{-5/3} \propto T^{5/3}\), so \(\operatorname{Var}(Z) \propto T^{8/3}\) and \(\sigma(\hat g) = \tfrac{6}{aT}\sqrt{\operatorname{Var}Z} \propto T^{1/3}\): rising with \(T\). Below the corner (long \(T\)), \(S_\varepsilon\) is flat, so \(\operatorname{Var}(Z) \propto T\) and \(\sigma(\hat g) \propto T^{-1/2}\): falling. The turnover is therefore a maximum. The exact integral (figure) puts it at \(T \approx 110\) s, with \(\sigma(\hat g)\) within \(1\%\) of that maximum for \(80 \le T \le 150\) s.

The dither period, continued What a longer period buys

The published \(T = 150\) s sits \(1\%\) below that maximum. Lengthening the period reduces the noise, and it also helps a second estimator this loop can run: \(\hat H\), the curvature read from the second harmonic \(2\omega\) of the same dither (the subject of the matrix-gain deck). The rotor is a low-pass filter with time constant \(\tau_{\text{rotor}} = 6.4\) s at \(8\) m/s (A.0c), and it attenuates that harmonic, so

\[\frac{\hat H}{J''} - 1 \;=\; \frac{1}{1+(2\omega\tau_{\text{rotor}})^2} - 1 \;<\; 0\]

At \(T = 600\) s the gradient noise is \(22\%\) lower (\(4.07\) against \(5.21\times10^{-7}\), the table on The answer) and this bias moves from \(-22\%\) to \(-2\%\).

What the derivation does not fix

Feed the published configuration [CLR19] through the closed-loop result: \(T = 150\) s, \(a = 13.6\%\) of \(u^\star\), \(J'' = -1.41\times10^{-13}\), and a loop that settles in about three periods, which for \(\tilde u_{n+1} = (1-K)\tilde u_n\) means \((1-K)^3 \approx 1/8\), \(K = 0.5\), so \(\sqrt{K/(2-K)} = 1/\sqrt3\). With \(\sigma(\tilde u) = \sigma(\hat g)/(\sqrt3\,|J''|)\) and \(d\ln\lambda/d\ln u = -1/3\) (Link A under the dither):

\[\sigma(u) = 0.96\,u^\star \qquad\Longrightarrow\qquad \sigma_\lambda = 2.4 \qquad\text{against}\qquad \sigma_\lambda \lesssim 0.14 \ \text{from the published table}\]

inflow spectrum: eliminated

Reconciling would need \(S_\varepsilon\) smaller by \(292\times\), i.e. \(\mathrm{TI} = 0.58\%\) against the stated \(10\%\).

power signal: eliminated

Demodulating generator power \(u\Omega^3\) rather than aerodynamic \(\tau_{\text{aero}}\Omega\) changes \(\sigma(\hat g)\) by \(2\%\) in simulation.

The factor of \(17\) is now analytic, independent of the solver, and two explanations are gone. What remains: \(|J''|\) wrong by that factor, an effective averaging far longer than the reported settling time permits, or an inflow genuinely far quieter at \(f_1\) than IEC Kaimal at \(10\%\). Figure 3 of [CLR19] plots the inflow PSD; its value at \(6.7\times10^{-3}\) Hz decides among them.

Where this leaves the program

established

\(\sigma(\hat g)\) from the inflow spectrum to \(6\%\), nothing fitted. The \(\lambda\) channel drops out at the optimum and only there.

actionable

\(T = 150\) s is near the worst available dither period. Lengthening it reduces gradient noise and curvature bias together.

open

A factor of \(17\) against the published scatter, with two of four candidate explanations eliminated.

Appendices

The derivations the main line leaned on without stopping to do.

A.0 Where the equilibrium relation comes from

The rotor is a rigid body of inertia \(I\), driven by aerodynamic torque and braked by the generator: \(I\dot\Omega = \tau_{\text{aero}} - \tau_g\) with \(\tau_g = u\Omega^2\). Power is torque times angular speed, and \(\Omega = \lambda V/R\), so with \(A = \pi R^2\)

\[\tau_{\text{aero}} = \frac{P}{\Omega} = \frac{\tfrac12\rho\pi R^2V^3C_P(\lambda)}{\lambda V/R} = \tfrac12\rho\pi R^3V^2\,\frac{C_P(\lambda)}{\lambda}\]

At rest \(\dot\Omega = 0\), so \(\tau_{\text{aero}} = \tau_g = u\lambda^2V^2/R^2\). The \(V^2\) on the two sides cancels, and dividing through:

\[\boxed{\;\frac{C_P(\lambda)}{\lambda^3} = \frac{2}{\rho\pi R^5}\,u =: c\,u\;}\]

A curve set by the blades meets a line set by the gain at \(\lambda_{\text{eq}}(u)\), and \(V\) is gone. That cancellation is the reason for the \(\Omega^2\) law: any other power of \(\Omega\) leaves a factor of \(V\) behind.

\(C_P/\lambda^3\) is not monotone; it peaks near \(\lambda \approx 4\). The physical root is on the falling branch, the other is unstable, and \(\lambda_{\text{eq}}(u)\) means the falling-branch root throughout.

A.0b Where the shape of \(C_P(\lambda)\) comes from

Three layers, from what is a theorem down to what is a fit.

physics fixes the envelope

\(C_P \ge 0\), and \(C_P \le 16/27 \approx 0.593\) for any rotor: the Betz limit, from momentum conservation on the stream tube alone. At \(\lambda \to 0\) the rotor is stopped and extracts nothing; at large \(\lambda\) the blades run faster than the wind can feed them and drag dominates. So \(C_P\) must rise from zero and fall back, and has at least one maximum.

blade design fixes the curve

For a given blade, \(C_P(\lambda)\) is computed by blade-element momentum theory: slice the blade spanwise, look up each section’s lift and drag against its local angle of attack from wind-tunnel airfoil tables, integrate. The result is empirical in its inputs and deterministic in its arithmetic. For a well-designed blade at fixed pitch it has one peak.

what this deck uses

Heier’s closed-form fit to typical BEM curves [H98], remapped so its peak sits at \((\lambda^\star, C_P^{\max}) = (7.5,\,0.49)\), the design values of the NREL 5 MW rotor from [J09]. A stand-in with the right peak and the right asymmetry, steep on the stall side and gentle above.

So “single maximum” is an empirical property of good blades, and the specific curve here is a fit. Two consequences the main line depends on survive any such curve: \(C_P' = 0\) at the peak, which is what Link A uses, and \(C_P/\lambda^3\) being non-monotone, which is why A.0 must pick a branch.

A.0c The rotor time constant

Write the rotor equation of A.0 at fixed \(u\) as \(I\dot\Omega = F(\Omega,V) := \tau_{\text{aero}}(\Omega,V) - u\Omega^2\), with equilibrium \(F(\Omega_{\text{eq}},\bar V) = 0\). For \(\delta\Omega := \Omega - \Omega_{\text{eq}}\) and \(\delta V := V - \bar V = \bar V\varepsilon\), to first order

\[I\,\dot{\delta\Omega} = F_\Omega\,\delta\Omega + F_V\,\delta V, \qquad F_\Omega := \partial F/\partial\Omega,\ F_V := \partial F/\partial V \ \text{at the equilibrium.}\]

From A.0, \(\tau_{\text{aero}} = \tfrac12\rho\pi R^3V^2\,C_P(\lambda)/\lambda\) with \(\lambda = R\Omega/V\), so \(\partial\lambda/\partial\Omega = R/V\) and

\[\frac{\partial\tau_{\text{aero}}}{\partial\Omega} = \tfrac12\rho\pi R^3V^2\,\frac{\lambda C_P'(\lambda) - C_P(\lambda)}{\lambda^2}\,\frac{R}{V} = \tfrac12\rho\pi R^4V\,\frac{\lambda C_P' - C_P}{\lambda^2}, \qquad 2u\Omega_{\text{eq}} = \frac{2\tau_{\text{aero}}}{\Omega_{\text{eq}}} = \rho\pi R^4V\,\frac{C_P}{\lambda^2}\]

the second using \(u\Omega_{\text{eq}}^2 = \tau_{\text{aero}}\) at equilibrium. Therefore

\[F_\Omega = -\tfrac12\rho\pi R^4\bar V\,\frac{3C_P - \lambda C_P'}{\lambda^2} < 0, \qquad \tau_{\text{rotor}} := -\frac{I}{F_\Omega} = \frac{2I\lambda^2}{\rho\pi R^4\bar V\,\big(3C_P(\lambda) - \lambda C_P'(\lambda)\big)}, \qquad b := \frac{F_V\bar V}{I}\]

At the peak \(C_P' = 0\), so \(\tau_{\text{rotor}} = 2I\lambda^{\star2}/(3\rho\pi R^4\bar V C_P^{\max}) = 6.4\) s for \(I = 4.05\times10^{7}\) kg m\(^2\) (rotor \(3.54\times10^{7}\) plus \(97^2\times534\) from the generator through the gearbox, [J09]) and \(\bar V = 8\) m/s. Dividing the linearized equation by \(I\) gives \(\dot{\delta\Omega} + \delta\Omega/\tau_{\text{rotor}} = b\,\varepsilon(t)\), the form on Link A, standing assumptions. The full equation is stable about the same equilibrium and forgets \(\Omega(0)\) in the same way; the linearization makes the rate explicit.

A.0c, continued The rotor as a low-pass filter

For a single-frequency input \(\varepsilon(t) = \cos(2\pi f t)\) the steady-state solution of \(\dot{\delta\Omega} + \delta\Omega/\tau_{\text{rotor}} = b\,\varepsilon\) is

\[\delta\Omega(t) = \frac{b\,\tau_{\text{rotor}}}{\sqrt{1+(2\pi f\tau_{\text{rotor}})^2}} \cos\big(2\pi f t - \varphi\big), \qquad \varphi = \arctan(2\pi f\tau_{\text{rotor}})\]

a first-order low-pass with corner frequency \(1/(2\pi\tau_{\text{rotor}}) = 0.025\) Hz.

Below the corner, \(\Omega\) follows the wind and \(\lambda\) stays near \(\lambda_{\text{eq}}(u)\); above it, \(\Omega\) barely moves and \(\lambda\) follows \(1/V\). These are the two limits bracketing \(\delta\lambda\) on Link A, executed, continued. A harmonic of the dither at frequency \(f\) is attenuated by \(1/\sqrt{1+(2\pi f\tau_{\text{rotor}})^2}\) and delayed by \(\varphi\); demodulating it in phase with the dither projects onto \(\cos\varphi = 1/\sqrt{1+(2\pi f\tau_{\text{rotor}})^2}\) a second time. The product at the second harmonic \(2\omega\) is \(1/(1+(2\omega\tau_{\text{rotor}})^2)\), the factor on The dither period has a worst case.

A.0d Why \(\bar\lambda - \lambda^\star\) is second order

Under the standing assumptions \(\lambda(t)\) is a functional of the wind path, and the next slide derives its expansion in powers of \(\varepsilon\) about \(\varepsilon \equiv 0\), where \(\lambda = \lambda_{\text{eq}}(u)\):

\[\lambda(t) = \lambda_{\text{eq}}(u) + \Lambda_1[\varepsilon](t) + \Lambda_2[\varepsilon](t) + O(\varepsilon^3)\]

\[\Lambda_1[\varepsilon](t) = \int_0^\infty h(s)\,\varepsilon(t-s)\,ds, \qquad \Lambda_2[\varepsilon](t) = \iint_0^\infty k(s_1,s_2)\,\varepsilon(t-s_1)\,\varepsilon(t-s_2)\,ds_1\,ds_2\]

with \(h\) and \(k\) fixed kernels set by the rotor equation. Taking expectations,

\[\mathbb{E}\,\Lambda_1 = \int_0^\infty h(s)\,\mathbb{E}[\varepsilon(t-s)]\,ds = 0, \qquad \big|\mathbb{E}\,\Lambda_2\big| = \Big|\iint k(s_1,s_2)\,\Gamma(s_1-s_2)\,ds_1ds_2\Big| \le \|k\|_{1}\,\Gamma(0) = \|k\|_{1}\,\sigma_\varepsilon^2\]

the first because \(\mathbb{E}[\varepsilon] = 0\), the second because \(|\Gamma(\tau)| \le \Gamma(0)\) (Cauchy-Schwarz). At \(u = u^\star\), \(\lambda_{\text{eq}}(u^\star) = \lambda^\star\), so \(\bar\lambda - \lambda^\star = \mathbb{E}\,\Lambda_2 + O(\mathbb{E}|\varepsilon|^3) = O(\sigma_\varepsilon^2)\).

In the two rotor limits. Frozen \(\Omega\): \(\lambda = \lambda^\star/(1+\varepsilon) = \lambda^\star(1 - \varepsilon + \varepsilon^2 - \dots)\), so \(\Lambda_1 = -\lambda^\star\varepsilon\) has mean zero and \(\Lambda_2 = \lambda^\star\varepsilon^2\) has mean \(\lambda^\star\sigma_\varepsilon^2 = 0.01\,\lambda^\star\). Tracking \(\Omega\): \(\lambda = \lambda^\star\) and the difference is zero. With the linearized rotor of A.0c, \(\lambda = \lambda^\star(1 + \delta\Omega/\Omega_{\text{eq}})(1 - \varepsilon + \varepsilon^2 - \dots)\); the first-order part \(\delta\Omega/\Omega_{\text{eq}} - \varepsilon\) has mean zero, and the second-order part \(\varepsilon^2 - \varepsilon\,\delta\Omega/\Omega_{\text{eq}}\) has mean \(\sigma_\varepsilon^2 - \mathbb{E}[\varepsilon\,\delta\Omega]/\Omega_{\text{eq}}\) with \(|\mathbb{E}[\varepsilon\,\delta\Omega]| = \big|b\int_0^\infty e^{-s/\tau_{\text{rotor}}}\Gamma(s)\,ds\big| \le b\,\tau_{\text{rotor}}\,\sigma_\varepsilon^2\).

A.0d, continued Where the expansion comes from

\(\Omega(t)\) solves \(I\dot\Omega = F(\Omega,V)\), \(F := \tau_{\text{aero}} - u\Omega^2\), driven by \(V = \bar V(1+\varepsilon)\), with \(F\) smooth and the equilibrium stable (\(F_\Omega < 0\), A.0c). Seek the solution ordered by powers of \(\varepsilon\), \(\Omega = \Omega_{\text{eq}} + \Omega_1 + \Omega_2 + \dots\) with \(\Omega_n\) of order \(n\), expand \(F\) about \((\Omega_{\text{eq}},\bar V)\), and collect equal orders:

  • Order 1. \(I\dot\Omega_1 = F_\Omega\Omega_1 + F_V\bar V\varepsilon\), the linearized equation of A.0c, so \[\Omega_1(t) = b\int_0^\infty e^{-s/\tau_{\text{rotor}}}\,\varepsilon(t-s)\,ds\] linear in the path, with kernel \(b\,e^{-s/\tau_{\text{rotor}}}\).
  • Order 2. \(I\dot\Omega_2 = F_\Omega\Omega_2 + \tfrac12F_{\Omega\Omega}\Omega_1^2 + F_{\Omega V}\bar V\,\Omega_1\varepsilon + \tfrac12F_{VV}\bar V^2\varepsilon^2\): the same stable operator, forced by products of two first-order quantities, so \[\Omega_2(t) = \frac1I\int_0^\infty e^{-s/\tau_{\text{rotor}}} \Big[\tfrac12F_{\Omega\Omega}\Omega_1^2 + F_{\Omega V}\bar V\,\Omega_1\varepsilon + \tfrac12F_{VV}\bar V^2\varepsilon^2\Big](t-s)\,ds\] With \(\Omega_1\) itself an integral of \(\varepsilon\), this is a double integral of \(\varepsilon(t-s_1)\varepsilon(t-s_2)\) against a fixed kernel \(k(s_1,s_2)\).
  • Order 3 and higher are forced by cubic and higher products: \(O(\varepsilon^3)\).

Time invariance of \(F\) makes the kernels depend on lags only; causality puts the integrals on \(s \ge 0\); stability makes \(e^{-s/\tau_{\text{rotor}}}\) integrable, so each \(\Omega_n\) is bounded by a constant times \(\max|\varepsilon|^n\) over the filter’s memory.

A.0d, continued From \(\Omega\) to \(\lambda\)

With \(R\Omega_{\text{eq}}/\bar V = \lambda_{\text{eq}}(u)\) and \(1/(1+\varepsilon) = 1 - \varepsilon + \varepsilon^2 - \dots\), \[\lambda = \frac{R\Omega}{\bar V(1+\varepsilon)} = \lambda_{\text{eq}}\Big(1 + \frac{\Omega_1}{\Omega_{\text{eq}}} + \frac{\Omega_2}{\Omega_{\text{eq}}} + \dots\Big)\big(1 - \varepsilon + \varepsilon^2 - \dots\big)\] \[\Lambda_1 = \lambda_{\text{eq}}\Big(\frac{\Omega_1}{\Omega_{\text{eq}}} - \varepsilon\Big), \qquad \Lambda_2 = \lambda_{\text{eq}}\Big(\frac{\Omega_2}{\Omega_{\text{eq}}} - \varepsilon\frac{\Omega_1}{\Omega_{\text{eq}}} + \varepsilon^2\Big)\]

These are the functionals \(\Lambda_1\), \(\Lambda_2\) of A.0d: \(\Lambda_1\) is linear in the path because \(\Omega_1\) is; \(\Lambda_2\) is quadratic because \(\Omega_2\), \(\varepsilon\Omega_1\) and \(\varepsilon^2\) are. The frozen rotor is the case \(\Omega_1 = \Omega_2 = 0\), which gives \(\Lambda_1 = -\lambda_{\text{eq}}\varepsilon\) and \(\Lambda_2 = \lambda_{\text{eq}}\varepsilon^2\) as on A.0d.

A.1 The estimator is a Fourier coefficient

Let \(y\) be \(T\)-periodic with Fourier series \(y = \tfrac{a_0}{2} + \sum_{k\ge1}(a_k\cos k\omega t + b_k\sin k\omega t)\), \(\omega := 2\pi/T\). Multiply by \(\sin\omega t\), integrate over one period, and use three product-to-sum identities:

\[\int_0^T\!\sin\omega t\,dt = 0, \qquad \int_0^T\!\cos k\omega t\,\sin\omega t\,dt = \tfrac12\!\int_0^T\!\big[\sin(k{+}1)\omega t - \sin(k{-}1)\omega t\big]dt = 0,\] \[\int_0^T\!\sin k\omega t\,\sin\omega t\,dt = \tfrac12\!\int_0^T\!\big[\cos(k{-}1)\omega t - \cos(k{+}1)\omega t\big]dt = \begin{cases} T/2, & k = 1\\ 0, & k \ne 1\end{cases}\]

because a non-constant sinusoid integrates to zero over whole periods, while for \(k = 1\) the first cosine is \(\cos0 = 1\). Every term dies except one:

\[\int_0^T y(t)\sin\omega t\,dt = b_1\cdot\frac{T}{2} \qquad\Longrightarrow\qquad b_1 = \frac{2}{T}\int_0^T y(t)\sin\omega t\,dt\]

Now the estimator, with the constants pulled out of the integral: \[\hat g = \frac{1}{T}\int_0^T y(t)\,\frac{2}{a}\sin\omega t\,dt = \frac{1}{a}\cdot\underbrace{\frac{2}{T}\int_0^T y(t)\sin\omega t\,dt}_{=\;b_1} = \frac{b_1}{a}\] The estimator is the first sine coefficient of the measured log-power, divided by the probe amplitude. Its \(2/a\) is the Fourier \(2/T\) times the \(1/a\) that turns a power swing into a slope.

A.1, continued From log-power to \(J\)

Two idealizations, both stated: \(\varepsilon \equiv 0\), so \(V = \bar V\) is constant; and the rotor sits at its equilibrium at every instant, so \(\lambda(t) = \lambda_{\text{eq}}(u(t))\) with \(u(t) = u + a\sin\omega t\). From What we need from the turbine, \(P = \tfrac12\rho A\bar V^3C_P(\lambda)\) and \(J(u) := \ln C_P(\lambda_{\text{eq}}(u))\). Therefore

\[y(t) = \ln P(t) = \ln\big(\tfrac12\rho A\bar V^3\big) + \ln C_P\big(\lambda_{\text{eq}}(u(t))\big) = \underbrace{\ln\big(\tfrac12\rho A\bar V^3\big)}_{\text{a constant}} + J\big(u + a\sin\omega t\big)\]

The constant contributes nothing to \(b_1\), because \(\int_0^T\sin\omega t\,dt = 0\) (A.1). So for the purpose of computing \(b_1\), \(y\) may be replaced by \(J(u + a\sin\omega t)\), which is what the next slide expands. Both idealizations are relaxed in the main line: the first by Links A through D, the second by the rotor-dynamics correction on The dither period has a worst case.

A.1, continued What that coefficient is

Taylor-expand \(J(u + a\sin\omega t)\) in \(a\):

\[J(u + a\sin\omega t) = J + J'a\sin\omega t + \frac{J''}{2}a^2\sin^2\omega t + \frac{J'''}{6}a^3\sin^3\omega t + O(a^4)\]

To read off \(b_1\), each power of \(\sin\omega t\) must be rewritten as a sum of harmonics. Two identities, the second from \(\sin3\theta = 3\sin\theta - 4\sin^3\theta\) rearranged:

\[\sin^2\theta = \tfrac12 - \tfrac12\cos2\theta, \qquad \sin^3\theta = \tfrac34\sin\theta - \tfrac14\sin3\theta\]

Now collect the coefficient of \(\sin\omega t\) term by term:

term its \(\sin\omega t\) content contribution to \(b_1\)
\(J\) none \(0\)
\(J'a\sin\omega t\) itself \(J'a\)
\(\tfrac{J''}{2}a^2\sin^2\omega t\) constant and \(\cos2\omega t\) only \(0\)
\(\tfrac{J'''}{6}a^3\sin^3\omega t\) \(\tfrac34\sin\omega t\) \(\tfrac{J'''}{6}\cdot\tfrac34\,a^3 = \tfrac{a^3}{8}J'''\)

\[b_1 = aJ' + \frac{a^3}{8}J''' + O(a^5) \qquad\Longrightarrow\qquad \hat g = \frac{b_1}{a} = J' + \frac{a^2}{8}J''' + O(a^4)\] Every even power of \(\sin\) is a constant plus cosines (A.2), so only odd powers of \(a\) reach \(b_1\), and after the division only even powers reach \(\hat g\). Same \(a^2/8\) as on The bias grows as \(a^2\), by a different route.

A.2 Moments of a sine over a whole period

With \(s = \sin\omega t\) and \(\langle\cdot\rangle = \tfrac1T\int_0^T(\cdot)\,dt\) over one period:

\[\langle s^{2m+1}\rangle = 0 \quad\text{for every } m \ge 0\]

because \(\sin\omega t\) is odd about \(t = T/2\) and the interval is symmetric about that point, so the two halves cancel. For the even moments,

\[\langle s^2\rangle = \Big\langle\tfrac{1-\cos2\omega t}{2}\Big\rangle = \tfrac12, \qquad \langle s^4\rangle = \Big\langle\tfrac{3 - 4\cos2\omega t + \cos4\omega t}{8}\Big\rangle = \tfrac38\]

using \(\sin^4\theta = \tfrac18(3 - 4\cos2\theta + \cos4\theta)\) and that every non-constant cosine averages to zero over a whole period. These are the only moments the bias derivation uses.

A.3 The window’s transform, in closed form

The window is \(w(t) = \sin\omega_1 t\) on \([0,T]\) with \(\omega_1 = 2\pi/T\), and zero elsewhere. Write \(\beta := 2\pi f\) and \(\sin\omega_1 t = (e^{i\omega_1 t} - e^{-i\omega_1 t})/2i\):

\[\hat w(f) = \int_0^T \sin\omega_1 t\,e^{-i\beta t}\,dt = \frac{1}{2i}\left[\frac{e^{i(\omega_1-\beta)T}-1}{i(\omega_1-\beta)} + \frac{e^{-i(\omega_1+\beta)T}-1}{i(\omega_1+\beta)}\right]\]

Since \(\omega_1T = 2\pi\), both exponentials equal \(e^{-i\beta T}\) and the fractions combine:

\[\boxed{\;\hat w(f) = \frac{\omega_1\big(1 - e^{-i2\pi fT}\big)}{\omega_1^2 - (2\pi f)^2}\;} \qquad\Longrightarrow\qquad |\hat w(f)|^2 = \frac{4\omega_1^2\sin^2(\pi fT)}{\big(\omega_1^2 - 4\pi^2f^2\big)^2}\]

At \(f = f_1 = 1/T\) numerator and denominator both vanish; the limit is \(|\hat w(f_1)|^2 = T^2/4\). The main lobe lies between the first two zeros of \(\sin^2(\pi fT)\), at \(f = 0\) and \(f = 2f_1\), which is the \([0, 2f_1]\) quoted on A shortcut that does not apply here. Checked against a numerical transform at four frequencies to all printed digits.

A.4 Parseval for the window

Parseval’s identity for square-integrable \(w\) states \(\int_{\mathbb{R}}|\hat w(f)|^2\,df = \int_{\mathbb{R}}|w(t)|^2\,dt\). The right side is elementary:

\[\int_{\mathbb{R}}|\hat w(f)|^2\,df = \int_0^T\sin^2\omega_1 t\,dt = \int_0^T\frac{1-\cos2\omega_1 t}{2}\,dt = \frac{T}{2}\]

the cosine integrating to zero over the two full cycles it completes on \([0,T]\). This fixes the total weight the window places on the spectrum, and it is what the narrowband shortcut in Link C uses when it pulls \(S\) outside the integral.

A.5 Stationary variance of a first-order autoregression

Let \(x_{n+1} = \alpha x_n + e_n\) with \(|\alpha| < 1\), \(e_n\) zero-mean with variance \(\sigma_e^2\) and uncorrelated with \(x_n\). If \(x\) is stationary then \(\operatorname{Var}(x_{n+1}) = \operatorname{Var}(x_n) =: \Sigma\), and taking variances of both sides,

\[\Sigma = \alpha^2\Sigma + \sigma_e^2 \qquad\Longrightarrow\qquad \Sigma = \frac{\sigma_e^2}{1-\alpha^2}\]

On The terminal error of the closed loop, \(\alpha = 1 - K\) and \(e_n = \kappa\hat g_{\text{noise},n}\), so \(\sigma_e^2 = \kappa^2\sigma^2(\hat g)\) and \(1-\alpha^2 = 1-(1-K)^2 = K(2-K)\). Substituting \(\kappa = K/|J''|\) gives the boxed result there. Stationarity requires \(|1-K| < 1\), that is \(0 < K < 2\), which is the loop’s stability condition.

A.6 Why a short record under-measures a persistent process

Take \(N\) consecutive samples of a stationary process with variance \(\sigma^2\) and autocorrelations \(\rho_k\), and form the sample variance \(s^2 := \tfrac1N\sum_{n=1}^N(x_n - \bar x)^2\) with \(\bar x := \tfrac1N\sum_n x_n\) the sample mean. Take the process zero-mean (subtract the mean otherwise). Expanding the square, \(\sum_n(x_n-\bar x)^2 = \sum_n x_n^2 - 2\bar x\sum_n x_n + N\bar x^2 = \sum_n x_n^2 - N\bar x^2\), so

\[\mathbb{E}[s^2] = \frac1N\sum_{n=1}^N\mathbb{E}[x_n^2] - \mathbb{E}[\bar x^2] = \sigma^2 - \operatorname{Var}(\bar x), \qquad \operatorname{Var}(\bar x) = \frac{\sigma^2}{N}\Big[1 + 2\sum_{k=1}^{N-1}\big(1-\tfrac kN\big)\rho_k\Big]\ \text{(A.9)}\]

For an AR(1), \(\rho_k = \alpha^k\). With \(\alpha = 0.985\) and \(N = 270\) the bracket is \(170\), so \(\mathbb{E}[s^2] = 0.628\,\sigma^2\) and the expected sample standard deviation is \(0.79\,\sigma\). A persistent series sampled for only a few correlation times has not had room to wander as far as it eventually will, and the sample variance records that. The measured/predicted ratio of \(0.777\) on the terminal-error slide is this effect, not a failure of the formula.

A.7 Variance per decade, and the peak of \(fS_\varepsilon\)

Variance per decade. With \(u := \ln f\), \(df = f\,du\), so \[\sigma_\varepsilon^2 = \int_0^\infty S(f)\,df = \int_{-\infty}^{\infty} f\,S(f)\,d(\ln f)\] On a logarithmic frequency axis the area under \(fS\) is the variance; the contribution of one decade is the area of \(fS\) over that decade.

The peak. Write \(c := 6L/\bar V\), so \(S_\varepsilon = \mathrm{TI}^2(4L/\bar V)(1+cf)^{-5/3}\) and \(fS_\varepsilon \propto f(1+cf)^{-5/3}\). Differentiate: \[\frac{d}{df}\Big[f(1+cf)^{-5/3}\Big] = (1+cf)^{-5/3} - \tfrac53\,cf\,(1+cf)^{-8/3} = (1+cf)^{-8/3}\Big[(1+cf) - \tfrac53 cf\Big] = (1+cf)^{-8/3}\Big[1 - \tfrac23 cf\Big]\] which vanishes at \(cf = 3/2\), that is \(f = \tfrac{3}{2c} = \tfrac{3\bar V}{12L} = \tfrac{\bar V}{4L}\), and changes sign from positive to negative there: a maximum. At \(\bar V = 8\) m/s, \(L = 340\) m this is \(5.9\times10^{-3}\) Hz, against the corner \(\bar V/6L = 3.9\times10^{-3}\) Hz.

A.8 Correlation between estimates \(k\) periods apart

The estimate from period \(k\) uses the window shifted by \(kT\), \(w_k(t) := w(t-kT)\), so \(Z_k := \int_{\mathbb{R}}\varepsilon(t)\,w_k(t)\,dt\). Substituting \(t' = t - kT\), the shift becomes a phase on the transform: \(\hat w_k(f) = \int w(t-kT)\,e^{-i2\pi ft}dt = e^{-i2\pi fkT}\int w(t')\,e^{-i2\pi ft'}dt' = e^{-i2\pi fkT}\,\hat w(f)\). Repeating Link C’s calculation with the two windows \(w\) and \(w_k\),

\[\operatorname{Cov}(Z_0,Z_k) = \iint w(t)\,w(s-kT)\,\Gamma(t-s)\,dt\,ds = \int_{\mathbb{R}} S^{2s}(f)\,\overline{\hat w(f)}\,e^{-i2\pi fkT}\hat w(f)\,df = \int_{\mathbb{R}} S^{2s}(f)\,|\hat w(f)|^2\,e^{-i2\pi fkT}\,df\]

\(S^{2s}|\hat w|^2\) is even in \(f\), so the sine part of \(e^{-i2\pi fkT}\) integrates to zero and \[\rho_k := \frac{\operatorname{Cov}(Z_0,Z_k)}{\operatorname{Var}(Z)} = \frac{\int S^{2s}(f)\,|\hat w(f)|^2\cos(2\pi fkT)\,df}{\int S^{2s}(f)\,|\hat w(f)|^2\,df}\] the formula on Averaging over \(N\) periods. Dividing by \(6/aT\) on both sides leaves the ratio unchanged, so \(\rho_k\) is also the correlation of \(\hat g_{\text{noise}}\) across periods.

A.9 Variance of the mean of \(N\) correlated samples

Let \(x_1,\dots,x_N\) have common variance \(\sigma^2\) and correlations \(\operatorname{Corr}(x_n,x_m) = \rho_{|n-m|}\), and \(\bar x := \tfrac1N\sum_n x_n\). Variance of a sum is the sum of all covariances:

\[\operatorname{Var}(\bar x) = \frac1{N^2}\sum_{n=1}^N\sum_{m=1}^N\operatorname{Cov}(x_n,x_m) = \frac{\sigma^2}{N^2}\sum_{n=1}^N\sum_{m=1}^N\rho_{|n-m|}\]

Count the pairs by lag \(k := |n-m|\): there are \(N\) pairs with \(k = 0\) and \(2(N-k)\) pairs with lag \(k \ge 1\) (each unordered pair twice). Hence

\[\operatorname{Var}(\bar x) = \frac{\sigma^2}{N^2}\Big[N + 2\sum_{k=1}^{N-1}(N-k)\rho_k\Big] = \frac{\sigma^2}{N}\Big[1 + 2\sum_{k=1}^{N-1}\Big(1-\frac kN\Big)\rho_k\Big]\]

With all \(\rho_k = 0\) the bracket is \(1\) and \(\operatorname{Var}(\bar x) = \sigma^2/N\), the independent case. Applied to \(\hat g_n\) with the \(\rho_k\) of A.8 this is the formula on Averaging over \(N\) periods; applied to the closed-loop iterates with \(\rho_k = \alpha^k\) it is the bracket in A.6.

References

[R17]  M. A. Rotea, Logarithmic power feedback for extremum seeking control of wind turbines, IFAC-PapersOnLine 50(1) (2017) 4504–4509. doi:10.1016/j.ifacol.2017.08.381

[CLR19]  U. Ciri, S. Leonardi & M. A. Rotea, Evaluation of log-of-power extremum seeking control for wind turbines using large eddy simulations, Wind Energy 22 (2019) 992–1002. doi:10.1002/we.2336

[KR22]  D. Kumar & M. A. Rotea, Wind turbine power maximization using log-power proportional-integral extremum seeking, Energies 15 (2022) 1004. doi:10.3390/en15031004

[KR24]  D. Kumar & M. A. Rotea, Optimal tip-speed ratio for degraded blades, Wind Energy Science 9 (2024) 2133–2146. doi:10.5194/wes-9-2133-2024

[RKAJ24]  M. A. Rotea, D. Kumar, E. J. Aju & Y. Jin, Multi-row extremum seeking for wind farm power maximization, J. Phys. Conf. Ser. 2767 (2024) 032043. doi:10.1088/1742-6596/2767/3/032043

[MGR24]  S. P. Mulders, N. Gallo & M. A. Rotea, Analysis of extremum seeking control for wind turbine torque controller optimization by aerodynamic and generator power objectives, ACC (2024). arXiv:2407.08059

[GKN12]  A. Ghaffari, M. Krstić & D. Nešić, Multivariable Newton-based extremum seeking, Automatica 48 (2012) 1759–1767. doi:10.1016/j.automatica.2012.05.059

[S00]  J. C. Spall, Adaptive stochastic approximation by the simultaneous perturbation method, IEEE Trans. Automat. Contr. 45 (2000) 1839–1853. PDF

[LM23]  C. Lauand & S. Meyn, Quasi-stochastic approximation: design principles with applications to extremum seeking control, IEEE Control Systems Magazine 43 (2023). doi:10.1109/MCS.2023.3291884

[A22]  N. J. Abbas et al., A reference open-source controller for fixed and floating offshore wind turbines, Wind Energy Science 7 (2022) 53–73. doi:10.5194/wes-7-53-2022

[H98]  S. Heier, Grid Integration of Wind Energy Conversion Systems, Wiley (1998). Source of the analytic \(C_P(\lambda,\beta)\) fit.

[J09]  J. Jonkman, S. Butterfield, W. Musial & G. Scott, Definition of a 5-MW reference wind turbine, NREL/TP-500-38060 (2009). PDF  ·  [IEC]  IEC 61400-1 ed. 3, Wind turbines, Part 1: Design requirements, normal turbulence model (2005). IEC

[KWIC72]  J. C. Kaimal, J. C. Wyngaard, Y. Izumi & O. R. Coté, Spectral characteristics of surface-layer turbulence, Q. J. R. Meteorol. Soc. 98, 563–589 (1972). DOI

Companion to Extremum Seeking Control of Wind Turbines.