CONTROL SYSTEM AND METHOD FOR ENERGY CAPTURE SYSTEM

Abstract

A computer-implemented method of controlling a power take-off (PTO) of an energy converter apparatus having at least one body which receives energy from their environment and whose energy is absorbed by said PTO is described. The method involves: (a) receiving an input describing the motion of the at least one body making up the apparatus; (b) predicting the excitation force F.sub.ex which will be incident on the body over a prediction horizon H.sub.p, by approximating said excitation force using generalised truncated half-range Chebyshev-Fourier (HRCF) basis functions; (c) solving an optimal control problem, defined in terms of optimising a cost function J representing the energy absorbed by the PTO over the prediction horizon H.sub.p, using a HRCF pseudo spectral optimal control wherein the state and control variables are approximated by their truncated half-range Chebyshev Fourier series to generate an optimal reference trajectory; (d) providing as an output to said PTO a control signal adapted to cause the PTO to approximate said optimal reference trajectory; and (e) repeating steps (a) to (d) after a calculation interval Tc where Tc<Hp.

Claims

1. A computer-implemented method of controlling a power take-off (PTO) of an energy converter apparatus having at least one body which receives energy from their environment and whose energy is absorbed by said PTO, the method comprising the steps of: (a) receiving an input describing the motion of the at least one body making up the apparatus; (b) predicting the excitation force F.sub.ex which will be incident on the body over a prediction horizon H.sub.p, by approximating said excitation force using generalised truncated half-range Chebyshev-Fourier (HRCF) basis functions; (c) solving an optimal control problem, defined in terms of optimising a cost function J representing the energy absorbed by the PTO over the prediction horizon H.sub.p, using a HRCF pseudospectral optimal control wherein the state and control variables are approximated by their truncated half-range Chebyshev Fourier series to generate an optimal reference trajectory; (d) providing as an output to said PTO a control signal adapted to cause the PTO to approximate said optimal reference trajectory; and (e) repeating steps (a) to (d) after a calculation interval Tc where Tc<Hp.

2. A computer-implemented method as claimed in claim 1, wherein step (c) comprises: (vi) expressing the motion of the body in terms of a projection matrix X1 defining the body's position, and a projection matrix X2 defining the body's velocity; (vii) defining a control variable vector and its projection matrix U.sup.T;

1. A computer-implemented method of controlling a power take-off (PTO) of an energy converter apparatus having at least one body which receives energy from their environment and whose energy is absorbed by said PTO, the method comprising the steps of: (a) receiving at a processor an input describing the motion of the at least one body making up the apparatus; (b) predicting at said processor the excitation force F.sub.ex which will be incident on the body over a prediction horizon H.sub.p, by approximating said excitation force using generalised truncated half-range Chebyshev-Fourier (HRCF) basis functions; (c) solving an optimal control problem at said processor, defined in terms of optimising a cost function J representing the energy absorbed by the PTO over the prediction horizon H.sub.p, using a HRCF pseudospectral optimal control wherein the state and control variables are approximated by their truncated half-range Chebyshev Fourier series to generate an optimal reference trajectory; (d) providing as an output from said processor to said PTO a control signal adapted to cause the PTO to approximate said optimal reference trajectory; and (e) repeating steps (a) to (d) after a calculation interval Tc where Tc<Hp.

2. The computer-implemented method as claimed in claim 1, wherein step (c) comprises: (vi) the motion of the body in terms of a projection matrix XI defining the body's position, and a projection matrix X2 defining the body's velocity; (vii) defining a control variable vector and its projection matrix U.sup.T; (viii) defining a cost function J=U.sup.TX.sub.2 (ix) determining a set of 2N+2 variables defining the N+1 components of XI and the N+1 components of X2 for a given projection vector U of the control variable by cancellation of residuals at a series of N+1 collocation points within the prediction horizon interval; and (x) solving the cost function J=U.sup.TX.sub.2 to determine an optimal trajectory [XI, X2, U] to optimise energy absorption while respecting predefined constraints on the position, velocity and control force.

3. The computer-implemented method as claimed in claim 1, further comprising the step of tracking the trajectory of said body using a real time controller.

4. The computer-implemented method as claimed in claim 1, wherein said energy conversion apparatus is a wave energy converter system.

5. The computer-implemented method as claimed in claim 4, wherein said step of predicting the excitation force comprises observing incident wave motion on the apparatus and generating a prediction from said observed motion.

6. The computer-implemented method as claimed in any preceding claim, wherein said prediction horizon H.sub.p is in the range of 2 to 100 seconds, preferably 5 to 50 seconds.

7. The computer-implemented method as claimed in claim 1, wherein the state variables are truncated as a series of N basis functions where N is between 5 and 100, more preferably between 10 and 50.

8. A processor programmed to implement the computer-implemented method of any preceding claim.

9. An energy conversion apparatus comprising at least one body which receives energy from its environment, a power take-off (PTO) configured to absorb energy from said body, and a processor having computer executable code configured to perform the steps of: (a) receiving at said processor an input describing the motion of the at least one body making up the apparatus; (b) predicting at said processor the excitation force F.sub.ex which will be incident on the body over a prediction horizon H.sub.p, by approximating said excitation force using generalised truncated half-range Chebyshev-Fourier (HRCF) basis functions; (c) solving an optimal control problem at said processor, defined in terms of optimising a cost function J representing the energy absorbed by the PTO over the prediction horizon H.sub.p, using a HRCF pseudospectral optimal control wherein the state and control variables are approximated by their truncated half-range Chebyshev Fourier series to generate an optimal reference trajectory; (d) providing as an output from said processor to said PTO a control signal adapted to cause the PTO to approximate said optimal reference trajectory; and (e) repeating steps (a) to (d) after a calculation interval Tc where Tc<Hp.

10. The energy conversion apparatus of claim 9, wherein step (c) comprises: (vi) the motion of the body in terms of a projection matrix XI defining the body's position, and a projection matrix X2 defining the body's velocity; (vii) defining a control variable vector and its projection matrix U.sup.T; (viii) defining a cost function J=U.sup.TX.sub.2 (ix) determining a set of 2N+2 variables defining the N+1 components of XI and the N+1 components of X2 for a given projection vector U of the control variable by cancellation of residuals at a series of N+1 collocation points within the prediction horizon interval; and (x) solving the cost function J=U.sup.TX.sub.2 to determine an optimal trajectory [XI, X2, U] to optimise energy absorption while respecting predefined constraints on the position, velocity and control force.

11. The energy conversion apparatus of claim 9, said computer executable code being further configured to perform the step of tracking the trajectory of said body using a real time controller.

12. The energy conversion apparatus of claim 9, wherein said energy conversion apparatus is a wave energy converter system.

13. The energy conversion apparatus of claim 12, wherein said step of predicting the excitation force comprises observing incident wave motion on the apparatus and generating a prediction from said observed motion.

14. The energy conversion apparatus of claim 9, wherein said prediction horizon H.sub.p is in the range of 2 to 100 seconds, preferably 5 to 50 seconds.

15. The energy conversion apparatus of claim 9, wherein the state variables are truncated as a series of N basis functions where N is between 5 and 100, more preferably between 10 and 50.

Description

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

[0031] FIG. 1 shows a flap-type wave energy controller (WEC) geometry. This generic WEC is based on the same working principle as the Oyster developed by Aquamarine Power Ltd (see M. Folley, T. Whittaker, and J. Vant Hoff, The design of small seabed-mounted bottom-hinged wave energy converters, in Proceedings of the 7th European Wave and Tidal Energy Conference, Porto, Portugal, vol. 455, 2007), recovering wave energy from the oscillatory surge motion of a bottom-hinged vertical panel.

[0032] The density of the sea water is fixed at .sub.water=1025 kg/m.sup.3. The geometric characteristics are W=30 m, T=1 m, h=H=15 m, with the body density is assumed to be uniform and equal to .sub.body=.sub.water/8. The total mass of the device is given by m=.sub.bodyWTL5810.sup.3 kg. The device inertia along the x-axis, expressed at the centre of gravity, is I.sub.G=m(T.sup.2+H.sup.2)/12 and is expressed at the centre of the pivot linkage, I.sub.0=I.sub.G+m(H/2).sup.2. The infinite frequency asymptote of the added inertia of the device, I.sub.=9.8110.sup.7 kgm.sup.2, and the kernel function K of the convolution product used to model the radiation force, are both determined using the potential code Aquaplus (described in G. Delhommeau, Seakeeping codes aquadyn and aquaplus, 19th WEGEMT School Numerical Simulation of Hydrodynamics: Ships and Offshore Structures, 1993), with K presented in FIG. 2.

[0033] Irregular incident waves are considered, and the elevation of the free surface is determined using a Bretschneider spectrum (C. Bretschneider, Wave variability and wave spectra for wind generated gravity waves, Tech. Memo, No. 118. Beach Erosion Board, US Army Corps of Eng., Washington D.C., 1959), with a bandwidth of 210.sup.2 Hz to 510.sup.1 Hz and a frequency step df=510.sup.3 Hz.

[0034] The excitation force is determined using the BEM code Aquaplus, based on the free surface elevation with randomly selected phases for the sinusoidal wave components.

[0035] The linear model is simulated in real time using a Runge-Kutta method, with a time step of t=10.sup.2 s.

[0036] A control algorithm, as described further below, is applied to control the PTO force, with a prediction horizon, under which the optimization process is done, chosen as T=15 s. The number of basis functions used to approximate the solution is N=15. The sequential quadratic programming algorithm is run for a maximum of 15 iterations, leading to an acceptable ratio between convergence quality and computational time, as shown in FIG. 3.

[0037] FIG. 4 shows the control algorithm structure in block diagram form. The structure comprises a hierarchically-organised upper and lower loop, which, respectively, generate and track a reference trajectory.

[0038] The upper loop predicts the excitation force over the prediction horizon, 12, then determines the reference trajectory, 14, using a pseudospectral method with half-range Chebyshev-Fourier (HRCF) basis functions, solved by a nonlinear programming method, such as sequential quadratic programming algorithm (SQP) (see M. J. Powell, A fast algorithm for nonlinearly constrained optimization calculations, in Numerical analysis. Springer, 1978, pp. 144-157) implemented in MATLAB via the fmincon function.

[0039] The lower loop facilitates tracking of the reference trajectory in real time. Trajectory generation is updated at a regular interval T.sub.s and tracked 16 using a standard backstepping method (see P. V. Kokolovie, The joy of feedback: nonlinear and adaptive,1992). A control force command is sent to the power take-off or PTO, 19, which imposes the control force on the WEC, 20. The resulting position and velocity as measured at the WEC are fed back to the trajectory tracking function 18 in the lower loop and to both the excitation force prediction function 12 and the trajectory generation function 14 in the upper loop.

[0040] In order to explain the operation of this algorithm in more detail, the mathematical basis for the algorithm is first established.

[0041] Problem Specification

[0042] The optimal control solution for renewable energy devices typically entails an energy production maximization over a given time horizon and under dynamical constraints dictating the allowable system behaviour. Consideration of physical limitations are also essential in the design of a realistic control algorithm and imply additional inequality constraints in the problem formulation. The cost function is defined by equation (2) and depends on the state and control functions, respectively x and u, of the considered system. The time horizon under which the absorbed energy is maximized is [t.sub.0,t.sub.0+T], where t.sub.0 R.sup.+ corresponds to the initial time and T R.sup.+ is the prediction time horizon. This prediction time horizon corresponds to the time window where the energy resource can be forecasted with an acceptable accuracy for use in a pseudospectral optimal control problem formulation. The interval [t.sub.0, t.sub.0+T] is mapped, for simplicity, into [1,1], using the affine transformation (1).

[00001]$\begin{array}{cc}=\frac{2}{T}.Math.\left(t-t_0\right)-1& \left(1\right)\end{array}$

[0043] where t [t.sub.0,t.sub.0+T] and [1,1]. We consider the following finite-horizon and constrained optimal control problem,

[00002]$\begin{array}{cc}\left\{\begin{array}{cc}\mathrm{Maximize}& J={}_{-1}^1.Math.F\left(x\left(\right),u\left(\right),\right).Math.d.Math..Math.\\ \mathrm{Subject}.Math..Math.\mathrm{to}& \begin{array}{c}\begin{array}{c}\stackrel{.}{x}\left(\right)=f\left(x,u,\right),\left[-1,1\right]\\ x\left(-1\right)=x_0\end{array}\\ H\left(x,u\right)0,\left[-1,1\right]\end{array}\end{array}\right\}& \left(2\right)\end{array}$

[0044] where F corresponds to the absorbed power, H corresponds to inequality constraints and f represents the device dynamics.

[0045] Control Solution Outline

[0046] Pseudospectral method: The method of weighted residuals (MWR) is widely used in various engineering applications to solve boundary-value or eigenvalue problems by expressing the approximated solution in a finite subspace of a Hilbert space. The MWR can be also used for optimal control problems involving the minimization, or maximization, of a particular cost function while insuring that the state and control variables respect a specific system of dynamical equations and inequality constraints. The MWR derives an approximate solution for optimal control problems, insuring equality and inequality constraints by annulment of the residuals. Projections of the residuals on a finite dimensional space, spanned by a particular set of orthogonal functions, called trial functions, are cancelled. Unlike MPC, which effectively uses local zero-order holder (ZOH) functions to approximate the optimal solution, spectral methods are generally based on global functions defined over the complete control horizon. Various basis functions can be chosen to obtain an approximation of the optimal solution, and the choice of a particular basis is discussed in section II-B2. The choice of trial functions dictates the type of spectral method use, for example the Galerkin method, or collocation methods.

[0047] The cost function minimization, or maximization, is realized using nonlinear programming methods (see Y. Nesterov, A. Nemirovskii, and Y. Ye, Interior-point polynomial algorithms in convex programming. SIAM, 1994, vol. 13), such as sequential quadratic program algorithms, leading to interesting computational aspects since the solution is sought in a finite dimensional vector space.

[0048] Function approximation: Pseudospectral methods are based on an approximation of the state and control variables into a N-dimensional vector space E generated by an orthogonal basis of real functions, B={.sub.1, . . . , .sub.N}. This approximation is usually realized using interpolation methods or truncated generalised Fourier series. Control and state variables are commonly rewritten as:

[00003]$\begin{array}{cc}{x}_i{x}_i^N=\underset{k=1}{\overset{N}{.Math.}}.Math.{}_k.Math.x_{\mathrm{ik}}=.Math..Math.X_i& \left(3\right)\\ u{u}^N=\underset{k=1}{\overset{N}{.Math.}}.Math.{}_k.Math.u_k=.Math..Math.U& \left(4\right)\end{array}$

[0049] with =[.sub.1, . . . , .sub.N], X.sub.i=[x.sub.i1, . . . , x.sub.iN].sup.T and U=[u.sub.1, . . . , u.sub.N].sup.T. A wide variety of basis functions can be used to approximate the state and control variables and the choice of a particular orthogonal base is mainly dictated by the specific requirements of the control problem.

[0050] Orthogonal wavelets are a large family of functions used for a wide variety of applications. For instance, Haar wavelets can be employed to solve nonlinear optimal control problems, Legendre wavelets constitute appropriate candidates for pseudospectral method basis sets for the resolution of various boundary value problems, and other wavelet families, such as Morlet wavelet, are widely used in signal processing in diversified domains. Orthogonal wavelet bases are typically generated through a scaling and shifting process leading to a significant number of derived basis functions.

[0051] Truncated Fourier series give satisfying approximations for relatively smooth functions and have been employed in pseudospectral optimal control, including in the wave energy field (see G. Bacelli and J. V. Ringwood, Numerical optimal control of wave energy converters, Sustainable Energy, IEEE Transactions on, vol. 6, no. 2, pp. 294-302, 2015). Fourier approximations are particularly well adapted to wave signal approximations and constitute a good first choice for pseudospectral resolution. For finite time horizon control, the Fourier basis requires periodicity of the approximated functions in order to avoid the Gibbs phenomenon on the boundaries. Since the Fourier basis generates only periodic functions, different boundary values of the approximated functions lead to discontinuities and higher frequency harmonics are needed to obtain a correct approximation. In F. Huankun and C. Liu, A buffer fourier spectral method for non-periodic pde, International Journal of numerical analysis and modeling, vol. 9, no. 2, pp. 460-478, 2012, a solution has been proposed to avoid such discontinuities by adding a buffer polynomial to construct an extended periodic function after which a standard Fourier pseudospectral method is applied. Lagrange polynomials are commonly used in Legendre methods, Chebyshev methods, or more generally in Jacobi methods, to interpolate or approximate the control and state variables. A certain amount of precaution has to be taken in the choice of the collocation points in order to avoid the Runge phenomenon during the interpolation.

[0052] In the present invention, to deal with realtime control and finite time control horizons requiring approximations of non-periodic functions in order to simulate both the transient and steady state responses of the device, we employ a basis employing HRCF functions, presented in more detail below.

[0053] By way of example, FIG. 5 shows a comparison between the approximations of a non-periodic function defined on [1,1] using different sets of basis functions. As an illustrative example, the approximated function presented in FIG. 6 is a sum of N.sub.s sine functions with random amplitudes, phases and periods uniformly distributed between, respectively, the following intervals [0,1],[0,2] and, [0.5,2]; N.sub.s is arbitrarily set to 10.sup.3.

[0054] FIG. 6 shows the approximation obtained with 25 functions such as ZOH functions or Haar wavelets, truncated Fourier series, Legendre polynomials and HRCF functions. While the legend in the original version of FIG. 6 is colour coded, the Legendre polynomials and HRCF functions give almost identical results and are indistinguishable from the initial function (heavy grey line). The truncated Fourier series is distinguishable due to the boundary discontinuities and the emergence of the Gibbs phenomenon. Haar wavelets, or ZOH functions both need a large number of basis functions to attain the same level of polynomial or HRCF approximation error. Thus for the smaller number of functions graphed in FIG. 6 for ZOH functions, these are readily distinguishable from the initial function, showing as a blocky, stepwise progression compared to the smooth initial function. As a result, a large number of ZOH functions have to be utilized for adequate performance, increasing the computational time for an optimal control problem using a standard MPC algorithm.

[0055] FIG. 6 shows how the error term drops away rapidly for HRCF representations of the FIG. 5 function with modest numbers of basis functions. Legendre representations fare reasonably well, but not as well as HRCF. Both Fourier and ZOH representations are unable to match the level of approximation of the HRCF representation with such modest numbers of basis functions.

[0056] 3) Solution route outline: As indicated when introducing FIG. 4, the method proposed herein comprises a hierarchically-organised upper and lower loop, which, respectively, generate and track a reference trajectory. The upper loop determines the reference trajectory using a pseudospectral method with HRCF basis functions, solved by a nonlinear programming method, while the lower loop facilitates tracking of the reference trajectory in real time. Accordingly, the definition and application of the HRCF basis for pseudospectral methods is now described.

[0057] Half-Range Chebyshev Fourier Functions

[0058] A. Half-Range Chebyshev Polynomials

[0059] In D. Huybrechs, On the fourier extension of nonperiodic functions, SIAM Journal on Numerical Analysis, vol. 47, no. 6, pp. 4326-4355, 2010, HRCF functions were introduced. An orthogonal basis for non-periodic functions is determined, based on half-range Chebyshev polynomials of the first and second kind. Half-range Chebyshev polynomials of the first and second kind of order k, T.sub.k.sup.h and U.sub.k.sup.h respectively, are orthogonal with lower order monomials with respect to the weights 1/(1y.sup.2), for the first kind, and (1y.sup.2) for the second kind, on the interval [0,1]. Definitions of T.sub.k.sup.h and U.sub.k.sup.h are given in B. Orel and A. Perne, Chebyshev-fourier spectral methods for nonperiodic boundary value problems, Journal of Applied Mathematics, vol. 2014, 2014, based on the Huybrechs 2010 work.

[0060] Definition 1. Let T.sub.k.sup.h(y) be the unique normalized sequence of orthogonal polynomials satisfying

[00004]$\begin{array}{cc}{}_0^1.Math.T_k^h\left(y\right).Math.y^l.Math.\frac{1}{\sqrt{1-y^2}}.Math.\mathrm{dy}=0,l=0,\mathrm{.Math.}.Math.,k-1& \left(5\right)\\ \frac{4}{}.Math.{}_0^1.Math.{T_k^h\left(y\right)}^2.Math.\frac{1}{\sqrt{1-y^2}}.Math.\mathrm{dy}=1& \left(6\right)\end{array}$

[0061] The set {T.sub.k.sup.h(y)}.sub.k=0.sup. is a set of half-range Chebyshev polynomials of the first kind.

[0062] Definition 2. Let U.sub.k.sup.h(y) be the unique normalized sequence of orthogonal polynomials satisfying

[00005]$\begin{array}{cc}{}_0^1.Math.U_k^h\left(y\right).Math.y^l.Math.\sqrt{1-y^2}.Math.\mathrm{dy}=0,l=0,\mathrm{.Math.}.Math.,k-1& \left(7\right)\\ \frac{4}{}.Math.{}_0^1.Math.{U_k^h\left(y\right)}^2.Math.\sqrt{1-y^2}.Math.\mathrm{dy}=1& \left(8\right)\end{array}$

[0063] The set {U.sub.k.sup.h(y)}.sub.k=0.sup. is a set of half-range Chebyshev polynomials of the second kind.

[0064] FIGS. 7A and 7B show half-range Chebyshev polynomials of the first and second kinds, respectively.

[0065] B. Solution to the Approximation Function Problem

[0066] One way to obtain the Fourier series of any nonperiodic function f L.sub.[1,1].sup.2 is to extend f to a periodic function g defined in a larger interval, in this case [2,2]. Determination of the Fourier extension of f can be stated as an optimization problem.

[0067] Problem 1. Let G.sub.n be the space of 4-periodic functions of the form

[00006]$\begin{array}{cc}g{G}_n.Math.\text{:}.Math..Math.g\left(\right)=\frac{a_0}{2}+\underset{k=1}{\overset{n}{.Math.}}.Math.a_k.Math.\cos .Math..Math.\frac{}{2}.Math.k.Math..Math..Math..Math.\mathrm{.Math.}+b_k.Math.\sin .Math..Math.\frac{}{2}.Math.k.Math..Math.& \left(9\right)\end{array}$

[0068] The Fourier extension of f to the interval [2,2] is the solution to the optimization problem

[00007]$\begin{array}{cc}{g}_n:=\arg .Math..Math.\underset{g{G}_n}{\min }.Math.{.Math.f-g.Math.}_{L_{\left[-1,1\right]}^2}& \left(10\right)\end{array}$

[0069] Huybrechs (supra) proved the existence and uniqueness of the solution of the problem 1, based on orthogonal polynomials called half-range Chebyshev Fourier polynomials. Two sets

[00008]${\left\{T_k^h\left(\cos .Math.\frac{}{2}.Math.\right)\right\}}_{k=0}^n.Math..Math.\mathrm{and}.Math..Math.{\left\{U_k^h\left(\cos .Math.\frac{}{2}.Math.\right).Math.\sin .Math.\frac{}{2}.Math.\right\}}_{k=0}^{n-1}$

are introduced, constituting an orthonormal basis for, respectively, the (n+1) and n dimensional spaces spanned by the 4-periodic, respectively, cosine and sine functions. The exact solution of Problem 1 is thus directly found by an orthogonal projection of f, as expressed in equation (11).

[00009]$\begin{array}{cc}{g}_n\left(\right)=\underset{k=0}{\overset{n}{.Math.}}.Math.a_k.Math.T_k^h\left(\cos .Math..Math.\frac{}{2}.Math.\right).Math..Math.\mathrm{.Math.}+\underset{k=0}{\overset{n-1}{.Math.}}.Math.b_k.Math.U_k^h\left(\cos .Math..Math.\frac{}{2}.Math.\right).Math.\sin .Math..Math.\frac{}{2}.Math..Math..Math.\mathrm{where},& \left(11\right)\\ a_k={}_{-1}^1.Math.f\left(\right).Math.T_k^h\left(\cos .Math..Math.\frac{}{2}.Math.\right).Math.d.Math..Math..Math..Math.\mathrm{and},& \left(12\right)\\ b_k={}_{-1}^1.Math.f\left(\right).Math.U_k^h\left(\cos .Math..Math.\frac{}{2}.Math.\right).Math.\sin .Math..Math.\frac{}{2}.Math..Math..Math.d.Math..Math.& \left(13\right)\end{array}$

[0070] As was shown by B. Orel, the sets of basis functions

[00010]${\left\{T_k^h\left(\cos .Math..Math.\frac{}{2}.Math.\right)\right\}}_{k=0}^n.Math..Math.\mathrm{and}.Math..Math.{\left\{U_k^h\left(\cos .Math..Math.\frac{}{2}.Math.\right).Math.\sin .Math..Math.\frac{}{2}.Math.\right\}}_{k=0}^n$

can be employed in the resolution of a non-periodic boundary value problem using pseudospectral methods (see B. Orel and A. Perne, Chebyshev-fourier spectral methods for nonperiodic boundary value problems, Journal of Applied Mathematics, vol. 2014, 2014).
C. Computation with Half-Range Chebyshev Polynomials

Differentiation Matrix:

[0071] Based on the HRCF functions, B. Orel developed in an efficient method of calculation of the derivatives of truncated HRCF series (see B. Orel and A. Perne Computations with half-range chebyshev polynomials, Journal of Computational and Applied Mathematics, vol. 236, no. 7, pp. 1753-1765, 2012). Let g.sub.n be the truncated HRCF series of a function f L.sub.[1,1].sup.2, so that

[00011]$\begin{array}{cc}f\left(\right)g_n\left(\right)=\underset{k=0}{\overset{n}{.Math.}}.Math.a_k.Math.T_k^h\left(\cos .Math..Math.\frac{}{2}.Math.\right).Math..Math.\mathrm{.Math.}+\underset{k=0}{\overset{n-1}{.Math.}}.Math.b_k.Math.U_k^n\left(\cos .Math..Math.\frac{}{2}.Math.\right).Math.\sin .Math..Math.\frac{}{2}.Math.& \left(14\right)\end{array}$

[0072] The derivative of g.sub.n is then defined by equation (15).

[00012]$\begin{array}{cc}\frac{d.Math..Math.g_n}{d.Math..Math.}=\underset{k=0}{\overset{n}{.Math.}}.Math.a_k^{}.Math.T_k^h\left(\cos .Math..Math.\frac{}{2}.Math.\right).Math..Math.\mathrm{.Math.}+\underset{k=0}{\overset{n-1}{.Math.}}.Math.b_k^{}.Math.U_k^n\left(\cos .Math..Math.\frac{}{2}.Math.\right).Math.\sin .Math..Math.\frac{}{2}.Math.& \left(15\right)\end{array}$

[0073] In matrix form, D R.sup.(n+1)(2n+1) is defined such that f=Df, where f=[a.sub.0, . . . , a.sub.n,b.sub.0, . . . , b.sub.n1].sup.T and f=[a.sub.0, . . . , a.sub.n,b.sub.0, . . . , b.sub.n1].sup.T (See Computations with half-range chebyshev polynomials, Journal of Computational and Applied Mathematics, vol. 236, no. 7, pp. 1753-1765, 2012.) The differentiation matrix D is written in the following form,

[00013]$\begin{array}{cc}D=\frac{}{2}.Math.\left(\begin{array}{cc}0& H_1\\ H_2& 0\end{array}\right)& \left(16\right)\end{array}$

where the matrices H.sub.1 R.sup.(n+1)n and H.sub.2 R.sup.n(n+1) are not expanded here for brevity; a complete exposition is to be found in B. Orel and A. Perne Computations with half-range chebyshev polynomials, Journal of Computational and Applied Mathematics, vol. 236, no. 7, pp. 1753-1765, 2012.

Function Multiplication:

[0074] In the same manner, the Orel & Perne (2012) paper presents a matrix formulation to determine the truncated series coefficients of the product p.sub.n=g.sub.n,h.sub.n of truncated series of an arbitrary function f L.sub.[1,1].sup.2 and a known function h L.sub.[1,1].sup.2. In matrix form,

a matrix F R.sup.(2n+1)(2n+1) is defined such that p=Ff, where p=[a.sub.0, . . . , a.sub.n,b.sub.0, . . . , b.sub.n1].sup.T and f=[a.sub.0, . . . , a.sub.n,b.sub.0, . . . , b.sub.n1].sup.T, with

[00014]$\begin{array}{cc}F=\left(\begin{array}{cc}{G}_1& G_2\\ G_3& G_4\end{array}\right),& \left(17\right)\end{array}$

where the matrix G.sub.1 R.sub.(n+1)(n+1), G.sub.2 R.sub.(n+1)n, G.sub.3 R.sub.n(n+1) and G.sub.4 R.sub.nn are functions of the coefficients of the truncated series g.sub.n that approximates the known function h, and are defined in the Orel & Perne (2012) paper.

Convolution Product:

[0075] The convolution product between an arbitrary causal function f L.sup.2(R) and a known causal function h L.sup.2(R) is defined by equation (18).

(f*h)(t)=.sub.0.sup.tf(s).Math.h(ts)ds, t .sup.+(18)

[0076] Convolution products can occur in dynamical equations describing physical systems and an illustration is given in the chosen wave energy converter application below where the fluid-structure interaction force, more specifically the radiation force, involves a convolution product with the radiation kernel function and the velocity of the floating body.

[0077] The determination of the radiation force is generally achieved using various approximation methods, generally leading to a state-space model of the convolution product, see T. Perez and T. I. Fossen, Time-vs. frequency-domain identification of parametric radiation force models for marine structures at zero speed, Modeling, Identification and Control, vol. 29, no. 1, pp. 1-19, 2008. For example, Prony's method (G. Prony, Essai experimental et analytique, etc, J. de l'Ecole Polytechnique, vol. 1, no. 2, pp. 24-76, 1975) approximates the radiation kernel function by a sum of complex exponentials, creating new state variables with their respective partial differential equations. Using such approximations will result in the unwelcome increase in the dimension of the function f(x,u), taking into account new radiation state variables.

[0078] In the present embodiment, a direct approximation of the radiation kernel function is achieved using the HRCF functions. The coefficients of the truncated HRCF series of the convolution product c=[a.sub.0, . . . , a.sub.n, b.sub.0, . . . , b.sub.n1].sup.T between an arbitrary function f L.sub.[1,1].sup.2 and the known kernel function is given by:

c=Pf, (19)

where f=[a.sub.0, . . . , a.sub.n, b.sub.0, . . . , b.sub.n1].sup.T is the vector of coefficients of the truncated HRCF series of f, and P R.sup.(2n+1)(2n+1) is the convolution matrix depending on the kernel function of the convolution product. A direct application on hydrodynamic radiation convolution product is presented in an appendix.

The Wave Energy Device Case

[0079] As discussed previously, real-time wave energy control algorithms deals with nonperiodic signals, since the prediction horizon is limited, and wave surface elevation, fluid-structure interaction forces and body motion are usually described using spectral forms. The choice of HRCF functions appears to fit all the requirements needed to solve the WEC optimal control problem with a pseudospectral method.

[0080] Furthermore, the exponential rate of convergence of the Fourier series of the extended periodic function offers a better rate of convergence than extensions realized with the use of cut-off functions (see D. Huybrechs, On the fourier extension of nonperiodic functions, SIAM Journal on Numerical Analysis, vol. 47, no. 6, pp. 4326-4355, 2010; and F. Huankun and C. Liu, A buffer fourier spectral method for non-periodic pde, International Journal of numerical analysis and modeling, vol. 9, no. 2, pp. 460-478, 2012).

A. WEC Model

[0081] For clarity, we consider a wave energy device with only one degree of freedom. The fluid is assumed to be inviscid and the flow incompressible, allowing the use of potential theory to determine fluid-structure interactions. The body displacement and the amplitude of the wave field are considered small enough to use linearised potential theory, leading to a linear equation of motion of the device, or Cummin's equation (see W. Cummins, The impulse response function and ship motions, DTIC Document, Tech. Rep., 1962), defined by:

(I+I.sub.){umlaut over (x)}+.sub.0.sup.tK(ts){dot over (x)}(s)ds+S.sub.hx=F.sub.ex+F.sub.c, (20)

where, x,{dot over (x)},{umlaut over (x)} are the body displacement and its first and second derivatives, F.sub.c is the control force, m is the mass or inertia of the body (depending on the type of degree of freedom considered), .sub. corresponds to the infinite frequency added mass of the device, K is the kernel function for the radiation convolution product, used to determine radiation forces from the velocity of the body, S.sub.h is the linearised hydrostatic stiffness and F.sub.ex corresponds to the wave excitation force, or Froude-Krylov force, experienced by the device. Recent control approaches for nonlinear wave energy converter models based on pseudospectral method are presented in G. Li, Nonlinear model predictive control of a wave energy converter based on differential flatness parameterisation, International Journal of Control, pp. 1-10, 2015.

[0082] From (20) and using the affine transformation

[00015]$g.Math.\text{:}.Math.\left[-1,1\right].fwdarw.\left[t_0,t_0+T\right],\frac{T}{2}.Math.\left(+1\right)+t_0$

we define the scaled excitation force f.sub.ex=F.sub.ex g, the kernel function of the radiation convolution product k=K g, position x.sub.1=x g, velocity x.sub.2={dot over (x)} g and control force u=F.sub.c g, leading to the system of differential equations (21):

[00016]$\begin{array}{cc}\{\begin{array}{c}{\stackrel{.}{x}}_2=\frac{T}{2}.Math.(f_{\mathrm{ex}}+u.Math..Math.\mathrm{.Math.}\\ -{}_{g^{-1}\left(0\right)}^{}.Math.k\left(-s+g^{-1}\left(0\right)\right).Math.x_2\left(s\right).Math.\mathrm{ds}.Math..Math.\mathrm{.Math.}\\ -S_h.Math.x_1)/\left(I+I_{}\right)\\ {\stackrel{.}{x}}_1=\frac{T}{2}.Math.x_2\end{array}& \left(21\right)\end{array}$

where g.sup.1(0)=12t.sub.0/T is the antecedent of 0 from the affine transformation g. The system of differential equations (21) defines the equality constraints in the optimal control problem (2), specifying the function f of the state variable x=[x.sub.1,x.sub.2].sup.T and the control variable u via the following equation,

[00017]$\begin{array}{cc}\stackrel{.}{x}=\left[\begin{array}{c}{\stackrel{.}{x}}_1\\ {\stackrel{.}{x}}_2\end{array}\right]=f\left(x,u,\right),\left[-1,1\right]& \left(22\right)\end{array}$

[0083] Practical limitations are considered on the body motion, control force and power output, to ensure feasibility of the control, and the safety and durability of the wave energy device in real sea conditions. For all [1,1],

[00018]$\begin{array}{cc}\{\begin{array}{c}.Math.x_1\left(\right).Math.X_{\max }\\ .Math.x_2\left(\right).Math.V_{\max }\\ .Math.u\left(\right).Math.U_{\max }\end{array},& \left(23\right)\end{array}$

where X.sub.max, V.sub.max, U.sub.max are real positive constants corresponding to position, velocity and control force limitations determined in advance to ensure feasibility and safety. Inequality constraints are rewritten in a matrix form using approximation functions, and [1,1],

[00019]$\begin{array}{cc}H\left(x,u\right)=\left(\left[\begin{array}{c}\left(\right)\\ -\left(\right)\end{array}\right].Math.I_3\right).Math.X-\left[\begin{array}{c}1\\ 1\end{array}\right].Math.X_{\max }0& \left(24\right)\end{array}$

where X=[X.sub.1.sup.T,X.sub.2.sup.2,U.sup.T].sup.T represents the projections of the state and control variables, X.sub.max=[X.sub.max,V.sub.maxU.sub.max].sup.T represents position, velocity and control force limitations, I.sub.3 the 3-dimension identity matrix and .Math. denotes the Kronecker product. Equation (24) specifies the linear function H from the optimal control formulation in (2).

[0084] Finally, the cost function, representing the absorbed energy, is written using the approximation functions, in the following form,

J=.sub.1.sup.1U.sup.T.sup.T()()X.sub.2dt (25)

[0085] The cost function J represents the recovered energy by the control through the control force u, determined by integrating instantaneous absorbed power over the control horizon. In the case of orthonormal approximation functions, with respect to the inner product, .sub.i,.sub.j=.sub.1.sup.1.sub.i().sub.j()dt, the cost function described in equation (25) is simplified into the form J=U.sup.TX.sub.2.

B. Cancellation of the Residuals

[0086] Replacing state variable derivatives terms and convolution terms coming from the determination of the radiation force in the expression of f from equation (21), one can express the two residual terms r.sub.1 and r.sub.2 as

[00020]$\begin{array}{cc}\{\begin{array}{c}{r}_1\left(\right)=\left(\right).Math.(\frac{2.Math.\left(I+I_{}\right)}{T}.Math.{\mathrm{DX}}_2-{\mathrm{PX}}_2.Math..Math.\mathrm{.Math.}\\ -S_h.Math.X_1-U)-f_{\mathrm{ex}}\left(\right)-c_0\left(\right)\\ r_2\left(\right)=\left(\right).Math.\left(X_2-\frac{2}{T}.Math.{\mathrm{DX}}_1\right)\end{array}& \left(26\right)\end{array}$

[0087] Trial functions used by pseudospectral methods are Dirac distributions, leading to the cancellation of residual terms at particular instants, or collocation points. The interval [1,1] is covered by N+1 collocation points i.e. Chebyshev points of the second kind.

[00021]$\begin{array}{cc}{}_i=-\cos \left(\frac{.Math..Math.i}{N}\right),i=0,1,\mathrm{.Math.}.Math.,N.& \left(27\right)\end{array}$

[0088] Cancellation of the first residual term r.sub.1 at the collocation points leads to N+1 equations. The cancellation of the second residual term r.sub.2 can be simplified, leading to N+1 equations linking the velocity and position projection vectors X.sub.1 and X.sub.2, as

[00022]$\begin{array}{cc}\{\begin{array}{cc}{r}_i\left({}_i\right)=0,& i=0,1,\mathrm{.Math.}.Math.,N.\\ X_2-\frac{2}{T}.Math.{\mathrm{DX}}_1=0& \end{array}& \left(28\right)\end{array}$

[0089] This leads to 2N+2 dynamical equations determining, for a given projection vector U of the control variable u, the 2N+2 variables composed of the N+1 components of X.sub.1 and X.sub.2.

C. Trajectory Tracking

[0090] Various control algorithms, such as PID (see Y. Choi and W. K. Chung, PID trajectory tracking control for mechanical systems. Springer, 2004, vol. 298), sliding control (see V. Utkin, J. Guldner, and J. Shi, Sliding mode control in electromechanical systems. CRC press, 2009, vol. 34) or backstepping control (see S. M. Rozali, M. Rahmat, and A. R. Husain, Backstepping design for position tracking control of nonlinear system, in Control System, Computing and Engineering (ICCSCE), 2012 IEEE International Conference on. IEEE, 2012, pp. 77-82), can be used to realise the tracking of the reference trajectory.

[0091] For the wave energy device case study, the chosen lower control loop is based on a backstepping method and is described now. The variables x.sub.1 and x.sub.2 refer, respectively, to the position and the velocity of the device, and the reference position and velocity trajectories are denoted, respectively, by x.sub.1,ref and x.sub.2,ref. x.sub.1,ref and x.sub.2,ref are known and determined previously by the upper control loop.

[0092] We can define the function V.sub.1 depending on the error e.sub.1=x.sub.1x.sub.1,ref, as

V.sub.11/2e.sub.1.sup.2 (29)

[0093] In order for V.sub.1 to be Lyapounov stable, we define the desired velocity xd=x.sub.2,ref.sub.1e.sub.1 and the function V.sub.2, as

V.sub.2=V.sub.1+1/2e.sub.2.sup.2, (30)

where e.sub.2=x.sub.2x.sub.2,ref. The derivative of V.sub.2 is

{dot over (V)}.sub.2=.sub.1e.sub.2.sup.2+e.sub.2(e.sub.1+{umlaut over (x)}{umlaut over (x)}.sub.d) (31)

[0094] Replacing the expression for the acceleration {umlaut over (x)} from Cummin's equation (20) in the derivative function V.sub.2, we obtain

[00023]$\begin{array}{cc}{\stackrel{.}{V}}_2=-{}_1.Math.e_1^2+e_2\left(e_1+\frac{1}{I+I_{}}.Math.\left(F_{\mathrm{ex}}+F_c-K*\stackrel{.}{x}-S_h.Math.x\right)-{\stackrel{\mathrm{.Math.}}{x}}_d\right)& \left(32\right)\end{array}$

where * denotes the convolution product. In order for V.sub.2 to be Lyapounov stable, the control force u has to be of the following form,

F.sub.c=(I+I.sub.)(e.sub.1{umlaut over (x)}.sub.d+.sub.2e.sub.2)F.sub.ex+S.sub.hx+K*{dot over (x)}(33)

Results of Simulation

[0095] Reverting now to the geometry shown in FIG. 1 and the control solution outlined in FIG. 4, the process described above in mathematical form is now outlined in flowchart form in FIG. 7.

[0096] In FIG. 7, the process is divided into two halves. On the left side 30, a number of pre-computed parameters are indicated. These are specific to the particular physical system being modelled, and include the computation 32 of the coefficients of the HRCF basis functions (equations 5 and 6 for HRCF polynomials of the first kind and equations 7 and 8 for the second kind), the computation 34 of the differentiation matrix D (equation 16) and convolution matrix P (equation 19 and appendix), and the definition 36 of the optimal control problem. The latter step involves generating the cost function to be used (equation 25), the generation of the residuals (equation 26), and of the path constraints matrix (equation 24).

[0097] In the right hand side of the flowchart 40, a real time loop is implemented. The excitation force is estimated 42 based on sensor data (e.g. a wave buoy). Then in a moving window with period T.sub.p, the current velocity and position of the moving body are set as the initial conditions, allowing definition of the convolution function c.sub.0(t), step 44.

[0098] In step 46, the excitation force is approximated using generalised truncated half-range Chebyshev-Fourier (HRCF) basis functions, allowing the previously described mathematical treatment to be applied to the system.

[0099] In step 48, the control problem is solved, taking the pre-computed parameters from the left-hand side 30 of the flowchart, cancelling the residuals and determining the maximum energy achievable over the future window, based on the relationship J=U.sup.TX, and the dynamical equations describing the relationship between the projection vector U of the control variable u, and the components of X.sub.1 and X.sub.2. The number of equations is (2N+1) for the dynamical equation linking U, X.sub.1 and X.sub.2. This corresponds to the first residual term R.sub.1 in equation 28, that we cancel at (2N+1) collocation points. There are also (2N+1) equations linking X.sub.1 and X.sub.2, which corresponds to the second residual term R.sub.2 in equation (28). Thus, the vectors U,X.sub.1 and X.sub.2 lead to 3*(2N+1) variables. Cancellation of R.sub.1 and R.sub.2, results in 2*(2N+1) equations. The algorithm has thus (2N+1) degrees of freedom to choose the optimal coefficients of the control force U.

[0100] The output 50 is a reference trajectory and reference control force which are passed to the real time controller 16 (FIG. 4). The real time controller uses these references to generate a control force command to the PTO, and in step 50, the real-time controller tracks the trajectory using the backstepping method.

[0101] The results are presented for the system of FIG. 1, for a significant wave height of H.sub.s=1 m and a peak wave period of T.sub.w=8 s, using a Bretschneider wave model. The computational time to solve one optimal control problem over the finite prediction horizon is less than 0.6 s for the complete simulation horizon, allowing its real-time application.

[0102] FIG. 8 shows the optimal trajectory for the position (FIG. 8(a)) and velocity (FIG. 8(b)) of the device, and the optimal control force (FIG. 8(c)), while the dashed curves represent the results from the real-time receding horizon control algorithm described above, for a 250 s simulation. The optimal trajectory can be theoretically calculated from the exact cancellation of the reactive terms in the equation of motion (20) by the optimal control force, i.e. stiffness, mass and added mass terms (see J. Falnes, Ocean waves and oscillating systems: linear interactions including wave-energy extraction. Cambridge university press, 2002). In the optimal case, the device is brought into resonance for all frequencies and thus a small variation of the control force can induce large displacement differences. The velocity and control force determined using the pseudospectral method are seen to be close to the optimal solution from FIG. 8. However, some differences between the optimal and calculated position of the wave energy converter, due to restoring force terms, are also evident. Nevertheless, the differences between the optimal trajectory and the actual position of the device under control do not significantly affect energy absorption, since the variations are of lower frequency than the wave excitation, and relate to (stored) potential energy that is not directly involved in the energy absorption.

[0103] FIG. 9 shows the optimal unconstrained trajectories (solid line) and actual trajectories of the device under control (dashed line). The horizontal dotted lines show the position under constraints (FIG. 9(a)) and control force under constraints (FIG. 9(b)) for a 250 s simulation.

[0104] FIG. 10 presents the optimal absorbed power and the maximal energy absorbed by the wave energy device, and under control. Despite variations from the optimal trajectory, specifically for the position and the control force, the control gives a large absorption rate, and the total efficiency reaches 98.2%, meaning that almost all the available energy is recovered.

[0105] It can be seen from these results that the application of pseudospectral methods for the real-time receding-horizon energy maximisation problem provides an excellent solution due to the choice of a basis function that can cater for both transient and steady-state response components. In particular, to approximate non-periodic functions, specific adjustments need to be made, for example using cut-off functions to ensure periodicity. HRCF functions complement Fourier analysis and satisfy the non-periodicity requirement without introducing bias, via the direct transformation of the input variable, and some additional computation, since HRCF functions are directly implementable in the pseudospectral approach.

[0106] The method described herein is adjustable in computational complexity in order to match real-time requirements. The number of basis functions employed and the number of iterations in the optimization process constitute variables that have to be adapted to the control hardware limitations. The number of basis functions corresponds to the dimension of the space where the optimization algorithm searches for a solution, hence, the computational time can be reduced if the solution depends only on few variables. Pseudospectral methods present a pragmatic approach to solve optimal control problems, and allows path limitations that ensure safety and durability of the wave energy device to be taken into account.

[0107] The solution according to the invention suits renewable energy systems, where the objective is to maximise converted energy, subject to the retention of stability and adherence to physical system constraints. In particular, the HRCF basis function set is well suited to the wave energy control problem, where the excitation force and system variables are well modelled using harmonic signals. The example case study shows that the developed real-time controller can achieve energy capture rates approaching the theoretical optimum, for a flap-type wave energy converter. The same approach can be used for other energy converter apparatuses, by following the same steps of parametrising the model and the excitation force in HRCF terms, approximating the state and control functions using truncated HRCF series, and solving the resulting optimal control problem by cancellation of the residuals.

Appendix: Convolution Product

[0108] Cummins' equation (20) involves the determination of a convolution product in order to estimate part of the radiation force acting on the device hull. The convolution product between the scaled velocity x.sub.2 of the device and the scaled radiation kernel function k, is split into two terms c.sub.0 and c,

[00024]$\begin{array}{cc}L\left(\right)={}_{g^{\left(-1\right)}\left(0\right)}^{}.Math.k\left(-s+g^{-1}\left(0\right)\right).Math.x_2\left(s\right).Math.\mathrm{ds}=c_0\left(\right)+c\left(\right)& \left(34\right)\\ \mathrm{where}& \\ \{\begin{array}{c}{c}_0\left(\right)={}_{g^{-1}\left(0\right)}^{-1}.Math.k\left(-s+g^{-1}\left(0\right)\right).Math.x_2\left(s\right).Math.\mathrm{ds}\\ c\left(\right)={}_{-1}^{}.Math.k\left(-s+g^{-1}\left(0\right)\right).Math.x_2\left(s\right).Math.\mathrm{ds}\end{array}& \left(35\right)\end{array}$

[0109] For s [g.sup.1(0),1], x.sub.2 corresponds to past values of the velocity, and is thus known. The first part of the convolution product, c.sub.0, is directly computed based on past velocity values. Replacing the state variable x.sub.2 in the second term c, by its HRCF truncated series, defined in equation (36), get

[00025]$\begin{array}{cc}{x}_2\left(\right)=\underset{k=0}{\overset{N}{.Math.}}.Math..Math.{}_k\left(\right).Math.x_{2.Math.k}=\left(t\right).Math.X_2,& \left(36\right)\\ \mathrm{and},& \\ c\left(\right)=\underset{k=0}{\overset{N}{.Math.}}.Math..Math.c_k\left(\right).Math.x_{2.Math.k},& \left(37\right)\\ \mathrm{with}& \\ c_k\left(\right)={}_{-1}^{}.Math.k\left(-s+g^{-1}\left(0\right)\right).Math.{}_k\left(s\right).Math.\mathrm{ds}.& \left(38\right)\end{array}$

[0110] The orthogonal projection of c.sub.k over each basis function .sub.j is p.sub.kj and defined by equation

p.sub.kj=.sub.1.sup.1c.sub.k().sub.j()d(39).

[0111] Finally, the matrix P=[p.sub.kj].sub.kj is used to define the radiation force, so that

[00026]$\begin{array}{c}L\left(\right)=c_0\left(\right)+\underset{k=0}{\overset{N}{.Math.}}.Math..Math.\underset{j=0}{\overset{N}{.Math.}}.Math..Math.{}_j\left(t\right).Math.p_{\mathrm{kj}}.Math.x_{2.Math.k}.Math..Math.\left(40\right)\\ =c_0\left(\right)+\left(\right).Math.{\mathrm{PX}}_2.Math..Math.\left(41\right)\end{array}$

[0112] The expression for the convolution term is replaced in the expression of the equality constraints (21), i.e. inside the function f.